戴 偉 陳柳楊 鄭利敏 楊明暉,*

(1湖北第二師范學院物理與機電工程學院, 武漢 430205;2中國科學院武漢物理與數(shù)學研究所, 武漢 430071)

應用多中心分區(qū)方法構建H3分子反應勢能面

戴 偉1,2陳柳楊2鄭利敏2楊明暉2,*

(1湖北第二師范學院物理與機電工程學院, 武漢 430205;2中國科學院武漢物理與數(shù)學研究所, 武漢 430071)

勢能面在分子反應動力學的研究中起著非常重要的作用. 本文提出了一種新的勢能面構建方法——多中心分區(qū)法. 通過比較London-Eyring-Polanyi-Sato (LEPS)勢能、多體展開勢能、置換對稱不變多項式三種方法, 確定了H3分子的最佳勢能函數(shù)表達形式, 并應用準經(jīng)典軌線方法分析了勢能面的合理性, 結果表明置換對稱不變多項式能很好地描述H3分子的勢能面特征. 結合置換對稱不變多項式和本文提出的多中心分區(qū)方法,可以有效改善H3分子勢能面的精度并可能推廣到高維反應勢能面.

勢能面擬合; 準經(jīng)典軌線法; 多中心分區(qū)

1 引 言

隨著實驗技術的快速發(fā)展, 分子反應動力學實驗方法1取得了長足進步, 理論研究能夠解釋并預測基元化學反應的動態(tài)結構、反應過程和反應機理,是實驗研究的重要補充手段. 勢能面是理論研究中的一個重要概念, 它描述反應體系中多個原子間的相互作用, 因此需要構造準確的勢能面以盡可能真實表達分子反應的勢能面,2通常構造勢能面的方法可以分為兩類, 一種是插值方法,3–23Collins等3–15在高維插值勢能面方面做了很多系統(tǒng)的工作, 還有應用三次樣條插值、16,17Shepard 插值、18,19Modified Shepard 插值20–22來構建分子反應勢能面的研究工作,通過高維插值構建勢能面, 需要較多的從頭算數(shù)據(jù)點. 另一種是函數(shù)擬合,23–30Corchado等24通過解析表達式擬合構建了H+CH4勢能面, Varandas等25通過雙多體展開(DMBE2001)構建了H3勢能面, Bowman等26通過選擇滿足置換對稱性的多項式作為勢能函數(shù), 構建了置換對稱不變多項式(PIP)勢能面, 以及最近得到廣泛應用的的神經(jīng)網(wǎng)絡方法.27–30選擇合適的勢能函數(shù)形式能夠有效地提高勢能面精度.

在勢能面構建過程中, 基于誤差是按區(qū)域分布這一特點, 我們認為如果按區(qū)域多中心分區(qū)擬合,將會有利于改善勢能面的精度. 基于這一思想, 本文發(fā)展了多中心分區(qū)擬合方法并應用到H3分子反應勢能面的構建中, 相比高維插值方法, 節(jié)約計算資源, 相比傳統(tǒng)的函數(shù)擬合方法, 提高勢能面的精度,該項工作可為研究較大體系的反應動力學提供一個新的思路.

2 H3分子勢能面的解析表達式

2.1 訓練集的數(shù)據(jù)來源及要求

本文H3結構與能量的訓練集數(shù)據(jù)由DMBE2001勢能面25產(chǎn)生. 其中R1, R2, R3取值區(qū)間為0.05–0.5 nm,步長為0.01 nm, 坐標的選取如圖1所示. 剔除能量高于40000 cm–1的能量點, 同時剔除不能構成幾何形狀的數(shù)據(jù)點.

2.2 勢能函數(shù)形式的選取

改善勢能面擬合精度的關鍵在于選擇合適的勢能函數(shù)形式. Connor31指出, 理想的勢能函數(shù)應滿足以下條件:

(1) 能準確描述反應物和產(chǎn)物通道的漸進性質;

(2) 能正確顯示出體系的對稱性;

圖1 H3體系坐標Fig.1 Coordinates for H3system

(3) 在有實驗數(shù)據(jù)或者(非經(jīng)驗的)理論數(shù)據(jù)的反應區(qū)域能準確描述體系的能量;

(4) 在沒有實驗數(shù)據(jù)和理論數(shù)據(jù)的反應區(qū)域能合理描述體系的行為;

(5) 保證各個區(qū)域光滑連接;

(6) 勢能函數(shù)及其微分的代數(shù)形式盡可能簡單; (7) 盡可能少的數(shù)據(jù)點就能達到擬合的精度.本文討論以下三種勢能函數(shù)在擬合H3反應勢能面的精度.

2.2.1 LEPS勢能

Eyring和Polanyi32利用London對三原子體系的量子力學勢能近似構建了London-Eyring-Polanyi (LEP)勢能面, Sato33又在他們的基礎上進行了修正,去掉了勢壘頂端不合理的勢阱, 得到了London-Eyring-Polanyi-Sato (LEPS)勢能面. LEPS也適合于多原子分子, 由于形式比較簡單, 參數(shù)少, 能夠反映出體系的主要特征, 在分子動力學發(fā)展的早期應用特別廣泛. 三原子體系的能量采用London形式表示:

其中, Q為庫能積分, J為交換積分,

式中1Ei為雙原子Morse函數(shù),3Ei為反Morse函數(shù), 表達式如下:

式中Di、和βi分別表示雙原子分子基態(tài)離解能、平衡核間距和光譜常數(shù), 可通過從頭計算確定, Si是針對不同體系引入的可調參數(shù), βi為光譜常數(shù), 可表示為

c為光速, ωe為振動光譜常數(shù), μ為折合質量.

為確定LEPS勢能中的可調參數(shù)Si, 掃描了一條均方根偏差(RMSD)隨Si變化的曲線, 由圖2可以看出, Si取0.10時, 均方根偏差最小. 構建的LEPS勢能面如圖3所示, 勢能中所用參數(shù)如表1所示.

2.2.2 多體展開勢

多體展開法34是構建多原子分子勢能面的常見方法之一, 它由單體項、二體項和多體項構成. 表達形式如下:

其中,

其中RAB為雙原子體系AB的鍵長, DeAB是離解能, ReAB是平衡核間距, αAB是Morse參數(shù). 同理可得AC、BC雙原子系統(tǒng)的Morse勢表達式. H3體系的離解能為0.293526 hartree, 平衡核間距為1.390 Bohr, Morse參數(shù)為8.10036.

式中的為非線性系數(shù), ρBC和 ρAC的表達式與此類似, d是擬合參數(shù), M是展開的階數(shù), 本文中M = 5.

圖2 均方根偏差(RMSD)隨Si參數(shù)的變化關系Fig.2 Relationship between the root mean square deviation (RMSD) and the Siparameter

圖3 H3分子30°的勢能等值圖Fig.3 Contour plots of potential energy surface for H3at 30°

表1 LEPS勢能參數(shù)數(shù)值Table1 Values of LEPS potential parameters

2.2.3 置換對稱不變多項式(PIP)方法

美國Emory大學化學系Joel M. Bowman教授在構建多原子分子勢能面時, 選擇滿足置換對稱性的多項式作為勢能函數(shù), 對勢能面精度有很大改善,其擬合表達式27如下:

式中a為可調參數(shù), 在本文中取1.8 bohr, n1, n2, n3= 1, 2,, 10且n1+ n2+ n3≤ 18, N為擬合表達式中待定參數(shù)總和, 文中取970項.

2.2.4 計算結果

為了進一步研究勢能函數(shù)帶來數(shù)值誤差的分布狀況, 本文統(tǒng)計了三種不同勢能函數(shù)在不同散射角度下的均方根偏差(見表2). 結果表明: (1) PIP方法能準確描述H3分子勢能面的結構特征, 統(tǒng)計誤差最小; (2) 所有勢能函數(shù)在60°散射時, 誤差最大, 分析原因可能是DMBE2001勢能面25提供的數(shù)據(jù)集在這一角度本身不能很好描述H3體系的相互作用, 需要通過從頭計算做進一步的驗證.

圖4、圖5、圖6是三種擬合表達式60°時的勢能面等值圖, 由圖可以看出, 勢能面整體趨勢區(qū)別不大, 第一種擬合表達式近程排斥區(qū)域最大, 第二種擬合表達式近程排斥區(qū)域最小, 在三體解離區(qū), 也存在一定差別, 其它區(qū)域三種形式基本吻合.

表2 不同散射角下勢能面的RMSDTable2 RMSD of potential energy surface under different scattering angles

圖7是三種擬合勢能面及LEPS勢能面的準經(jīng)典軌線計算結果, 總體來說, 與文獻提供的勢能面在

圖4 公式1擬合下60°時的勢能面等值圖Fig.4 Contour plots of potential energy surface for the Function 1 at 60°

圖5 公式2擬合下60°時的勢能面等值圖Fig.5 Contour plots of potential energy surface for the Function 2 at 60°

圖6 公式3擬合下60°時的勢能面等值圖Fig.6 Contour plots of potential energy surface for the Function 3 at 60°

圖7 不同勢能函數(shù)表達式的準經(jīng)典軌線結果比較Fig.7 Comparison of the quasi-classical trajectory results between different potential energy functions

基本趨勢上是一致的, LEPS勢能面參數(shù)少, 使用方便, 但高能區(qū)域偏差較大, 第一種擬合表達勢能面在低能區(qū)域偏差較大, 不適宜用來描述H3分子相互作用, 第二種和第三種勢能面在低能區(qū)與文獻值符合很好, 高能區(qū)域, 第三種勢能面符合更好, 綜合動力學結果分析, PIP方法能更好地描述H3體系的相互作用.

3 應用多中心分區(qū)擬合構建H3分子勢能面

3.1 多中心分區(qū)擬合的思想由來

圖8是本文構建勢能面與數(shù)據(jù)測試集的差分圖(散射角60°), 圖中表示的是采用PIP方法構建的勢能面與測試數(shù)據(jù)的差值, 用來描述勢能面的回歸特性. 由圖可以看出, 誤差呈區(qū)域分布, 大部分區(qū)域吻合較好, 誤差不大, 只有在鞍點區(qū)域到三體離解區(qū)域過渡時, 兩種勢能差別較大. 擬合勢能面與測試數(shù)據(jù)集的吻合度按區(qū)域分布, 有些地方高于測試值,有些地方低于測試值, 鑒于誤差是按區(qū)域分布這一特點, 我們預測如果在構建勢能面時, 按區(qū)域多中心分區(qū)擬合, 對于提高勢能面的精度會有改善.

圖8 置換對稱不變多項式勢能面與測試數(shù)據(jù)的差值等高圖Fig.8 Contour plots of the difference between the permutation-invariant polynomial method and the test data

3.2 多中心展開構建H3勢能面

3.2.1 坐標的選取

為了描述H + H2反應體系, 選取H2為中心原子,通過兩個鍵長(R1, R2)和一個鍵角(θ)表示, 坐標選取如圖1所示.

3.2.2 計算方法

3.2.2.1 訓練集數(shù)據(jù)的選取

采用DMBE2001勢能面25作為測試集, 取R1, R2∈ [0.05 nm, 0.5 nm]均勻100 個點; 均勻取θ ∈[0°, 180°]30個點, 得到各數(shù)據(jù)點的能量值.

3.2.2.2 分區(qū)原則

在H + H2→ H2+ H的最小反應路徑上取7個(R1, R2)(單位: nm), 作為七個中心點, 計算訓練集中的數(shù)據(jù)點(R1, R2)與七個中心的距離, 依據(jù)最近及次近距離將數(shù)據(jù)集分成7個區(qū)域, 當R1= R2時, 存在兩個次小數(shù)據(jù)點與各中心點間距值, 將此點納入相應的兩個分區(qū), 擬合數(shù)據(jù)的分布與誤差統(tǒng)計如表3所示.

3.2.2.3 邊界處理

合并擬合分區(qū), 將各數(shù)據(jù)點根據(jù)公式(18)計算相應的能量值: centq為此數(shù)據(jù)點所在數(shù)據(jù)域的中心點, q = 1, 2, 3, , 7. wk為分區(qū)能量在總能量中的權重, 表示如下:

表3 擬合數(shù)據(jù)分布及誤差統(tǒng)計Table3 Distribution of the fitting data and the error statistics

圖9 多中心分區(qū)法計算準經(jīng)典軌線的結果Fig.9 Results of quasiclassical trajectory by “multi-center partition” method

當此數(shù)據(jù)點僅存在于兩個分區(qū)時, w3= 0, vfitj3= 0, vfitjk表示處在第j個分區(qū)時，第k次近中心點對總勢能面的貢獻. 則多中心分區(qū)能量(單位為cm–1)為:

3.3 多中心展開構建H3勢能面

采用多中心分區(qū)擬合后, 構建勢能面的準經(jīng)典軌跡計算結果與文獻完全一致, 如圖9所示. 由圖9可以看出, 多中心分區(qū)擬合后, 構建勢能面與原勢能面回歸性很好, 相比全域擬合(圖7), 精度進一步提高, 無論是在低能區(qū)域, 還是在高能區(qū)域, 都與文獻值完全吻合.

4 結 論

構建了H3體系的LEPS勢能面, 采用三種不同勢能函數(shù), 構建了H3分子勢能面, 為多中心分區(qū)方法確定了最佳擬合函數(shù)形式, 應用準經(jīng)典軌跡法驗證了勢能面的合理性, 最后采用多中心分區(qū)擬合, 構建了H3分子勢能面. 研究結果表明:

(1) LEPS勢能面參數(shù)少, 除Si參數(shù)需調節(jié)確定,其余參數(shù)都可通過從頭計算得到, 勢能函數(shù)形式簡單, 使用方便, 在定性計算中可以廣泛使用, 但不適宜描述高精度反應勢能面, 動力學驗證誤差較大.

(2) 比較了三種不同的勢能函數(shù)形式. 訓練集采用了213200個構型的能量信息, 截斷能為40000 cm–1,數(shù)據(jù)量大, 能域寬. PIP方法統(tǒng)計誤差最小, 動力學信息最吻合, 可以用來描述H3分子勢能面的能量特征.

(3) 擬合勢能面的精度與表達式展開的項數(shù)、數(shù)據(jù)集的數(shù)據(jù)量大小密切相關. 多中心分區(qū)之后,每一區(qū)域數(shù)據(jù)集變少, 改變了數(shù)據(jù)集和擬合參數(shù)的相對比, 可以降低擬合表達式展開項的要求, 同時有效改善擬合精度. 從計算結果可以看出: 多中心分區(qū)擬合法相比全域擬合, 勢能面回歸性很好, 無論是在低能區(qū)域, 還是在高能區(qū)域, 都與文獻值吻合很好.

(1)Sun, Z. F.; Gao, Z.; Wu, X. K.; Tang, G. Q.; Zhou, X. G.; Liu, S. L. Acta Phys. -Chim. Sin. 2015, 31, 829. [孫中發(fā), 高 治, 吳向坤, 唐國強, 周曉國, 劉世林. 物理化學學報, 2015, 31, 829.] doi: 10.3866/PKU.WHXB201503041

(2)Li, W.; Qu, J. Y.; Zhao, X. S. Acta Phys. -Chim. Sin. 2003, 19 (8), 751. [李 巍, 屈軍艷, 趙新生. 物理化學學報, 2003, 19 (8), 751.] doi: 10.3866/PKU.WHXB20030816

(3)Le, H. A.; Frankcombe, T. J.; Collins, M. A. J. Phys. Chem. A 2010, 114 (40), 10783. doi: 10.1021/jp1060182

(4)Ramazani, S.; Frankcombe, T. J.; Andersson, S.; Collins, M. A. J. Chem. Phys. 2009, 130 (24), 244302. doi: 10.1063/1.3156805

(5)Moyano, G. E.; Jones, S. A.; Collins, M. A. J. Chem. Phys. 2006, 124 (12), 124318. doi: 10.1063/1.2181571

(6)Moyano, G. E.; Collins, M. A. Theor. Chem. Acc. 2005, 113 (4), 225. doi: 10.1007/s00214-004-0626-8

(7)Evenhuis, C. R.; Collins, M. A.; Lin, X.; Zhang, O. H. Amer. Chem. Soc. 2004, 227, 259.

(8)Moyano, G. E.; Pearson, D.; Collins, M. A. J. Chem. Phys. 2004, 121 (24), 12396. doi: 10.1063/1.1810479

(9)Castillo, J. F.; Aoiz, F. J.; Ba?ares, L.; Collins, M. A. J. Phys. Chem. A 2004, 108 (32), 6611. doi: 10.1021/jp048366b

(10)Moyano, G. E.; Collins, M. A. J. Chem. Phys. 2003, 119 (11), 5510. doi: 10.1063/1.1599339

(11)Fuller, R. O.; Bettens, R. P. A.; Collins, M. A. J. Chem. Phys. 2001, 114 (24), 10711. doi: 10.1063/1.1377602

(12)Song, K.; Collins, M. A. Chem. Phys. Lett. 2001, 335 (5), 481.

(13)Collins, M. A.; Zhang, D. H. J. Chem. Phys. 1999, 111 (22), 9924. doi: 10.1063/1.480344

(14)Bettens, R. P. A.; Hansen, T. A.; Collins, M. A. J. Chem. Phys. 1999, 111 (14), 6322. doi: 10.1063/1.479937

(15)Jordan, M. J. T.; Collins, M. A. J. Chem. Phys. 1996, 104 (12), 4600. doi: 10.1063/1.471207

(16)Li, H.; Le, R. R. J. J. Chem. Phys. 2006, 125 (4), 044307. doi: 10.1063/1.2212933

(17)Li, A. Y.; Xie, D. Q. J. Chem. Phys. 2010, 133 (14), 144306. doi: 10.1063/1.3490642

(18)Wu, T.; Manthe, U. J. Chem. Phys. 2003, 119 (1), 14. doi: 10.1063/1.1577328

(19)Ishida, T.; Schatz, G. C. J. Chem. Phys. 1997, 107 (9), 3558. doi: 10.1063/1.474695

(20)Takata, T.; Taketsugu, T.; Hirao, K.; Gordon, M. S. J. Chem. Phys. 1998, 109 (11), 4281. doi: 10.1063/1.477032

(21)Crespos, C.; Collins, M. A. J. Chem. Phys. 2004, 120 (5), 2392. doi: 10.1063/1.1637337

(22)Wang, M.; Sun, X.; Bian, W.; Cai, Z. J. Chem. Phys. 2006, 124 (23), 234311. doi: 10.1063/1.2203610

(23)Dawes, R.; Thompson, D. L.; Guo, Y.; Wagner, A. F.; Minkoff, M. J. Chem. Phys. 2007, 126 (18), 184108. doi: 10.1063/1.2730798

(24)Corchado, J. C.; Bravo, J. L.; Espinosa-Garcia, J. J. Chem. Phys. 2009, 130 (18), 184314. doi: 10.1063/1.3132223

(25)Varandas, A. J. C.; Brown, F. B.; Mead, C. A. J. Chem. Phys. 1987, 86, 6258. doi: 10.1063/1.452463

(26)Bowman, J. M.; Czako, G.; Fu, B. Phys. Chem. Chem. Phys. 2011, 13 (18), 8094. doi: 10.1039/c0cp02722g

(27)Pukrittayakamee, A.; Malshe, M.; Hagan, M.; Raff, L. M.; Narulkar, R.; Bukkapatnum, S.; Komanduri, R. J. Chem. Phys. 2009, 130 (13), 134101. doi: 10.1063/1.3095491

(28)Le, H. M.; Raff, L. M. J. Phys. Chem. A 2010, 114 (1), 45. doi: 10.1021/jp907507z

(29)Sumpter, B. G.; Noid, D. W. Chem. Phys. Lett. 1992, 192 (5), 455.

(30)Blank, T. B.; Brown, S. D.; Calhoun, A. W.; Doren, D. J. J. Chem. Phys. 1995, 103 (10), 4129. doi: 10.1063/1.469597

(31)Connor, J. N. L. Comput. Phys. Commun. 1979, 17, 117. doi: 10.1016/0010-4655(79)90075-4

(32)Eyring, H.; Polanyi, M. Phys. Chem. Abt. B 1931, 12, 279.

(33)Sato, S. J. Chem. Phys. 1955, 23 (12), 2465.

(34)Aguado, A.; Paniagua, M. J. Chem. Phys. 1992, 96 (2), 1265. doi: 10.1063/1.462163

Application of the Multi-Center Partition Method to Construct the Potential Energy Surface of H3

DAI Wei1,2CHEN Liu-Yang2ZHENG Li-Min2YANG Ming-Hui2,*
(1School of Physics and Mechanical & Electrical Engineering, Hubei University of Education, Wuhan 430205, P. R. China;2Wuhan Institute of Physics and Mathematics, Chinese Academy of Sciences, Wuhan 430071, P. R. China)

The potential energy surface plays an important role in studying molecular reaction dynamics. In this work, a new method, namely the “multi-center partition” method, is proposed to construct the potential energy surface of H3. The optimized function is first determined by comparing the London-Eyring-Polanyi-Sato (LEPS) potential, the many-body expansion potential, and the permutation-invariant polynomial potential. This comparison shows that the permutation-invariant polynomial fitting proposed by Bowman is the most efficient method for describing the topology of the H3system. The quasi-classical trajectory method is used to analyze the rationality of those potential energy surfaces. By combining the multi-center partition method with the permutation-invariant polynomial method, the accuracy of the H3molecular potential energy surface is greatly improved and could possibly be used in the fitting of potential energy surfaces in other systems.

Potential energy surface fitting; Quasi-classical trajectory method; Multi-center partition

O641.3

10.3866/PKU.WHXB201509143

Received: July 6, 2015; Revised: September 14, 2015; Published on Web: September 14, 2015.

*Corresponding author. Email: yangmh@wipm.ac.cn; Tel: +86-27-87197783.

This project was supported by the China Scholarship Council (201408420174), Science and Technology Research Program of the Education Department, Hubei Province, China (Q20133005), and Natural Science Foundation of Hubei Province, China (2014CFB428, 2015CFB502).

國家留學基金委項目(201408420174), 湖北省教育廳科學技術項目(Q20133005)和湖北省自然科學基金(2014CFB428, 2015CFB502)資助

?Editorial office of Acta Physico-Chimica Sinica