摘? ?要： 采用分子動力學(xué)方法，對SPC、SPCE、TIP3P、TIP4P、TIP5P五種勢能模型下水分子的熱力學(xué)性質(zhì)、徑向分布函數(shù)、自擴散系數(shù)進行計算，優(yōu)選最適合于水分子動力學(xué)性質(zhì)研究的勢能模型。研究發(fā)現(xiàn)：1）采用SPCE模型通過模擬得到的內(nèi)能、汽化焓最為準確，其他模型受到了模型固有性質(zhì)和系綜（大量宏觀上完全相同的體系的抽象集合）的影響，模擬結(jié)果不甚理想;2）采用SPCE、TIP4P、TIP5P模型可以較好地反映水分子短程有序、長程無序的微觀結(jié)構(gòu);3）采用SPCE模型通過模擬得到的自擴散系數(shù)相對誤差僅為0.39%，遠好于其他模型。結(jié)果表明，SPCE模型能夠較真實地反映水分子的熱力學(xué)性質(zhì)、微觀結(jié)構(gòu)和自擴散特性等動力學(xué)性質(zhì)。

關(guān)鍵詞： 水分子;勢能模型;動力學(xué);徑向分布函數(shù);自擴散系數(shù)

中圖分類號：O645? ? 文獻標識碼：A? ? 文章編號：2095-8412 （2020） 05-085-05

工業(yè)技術(shù)創(chuàng)新 URL： http：//gyjs.cbpt.cnki.net? ? DOI： 10.14103/j.issn.2095-8412.2020.05.016

引言

水是生命的源泉，水分子的研究對水的開發(fā)利用具有非常重要的意義。分子動力學(xué)（MD）方法是模擬水分子動力學(xué)性質(zhì)的一種非常有效的方法，一直受到國內(nèi)外研究者的關(guān)注[1-3]。

單個水分子看似結(jié)構(gòu)比較簡單，但氫鍵的作用使得水分子相互聚合，形成多分子聚合物，水分子的微觀結(jié)構(gòu)發(fā)生了復(fù)雜化。采用分子動力學(xué)方法研究水分子，最重要的也是選擇分子間的相互作用勢。在水分子動力學(xué)計算中，可以選擇多種勢能模型，包括量子力學(xué)從頭計算模型（如MCY模型[4]）、半經(jīng)驗?zāi)Ｐ停ㄈ鏢T2、SPC、SPCE、TIP3P、TIP4P、TIP5P）等[5-8]。MCY模型結(jié)構(gòu)較為復(fù)雜，計算時間成本較高，因此在實際應(yīng)用中更多采用半經(jīng)驗?zāi)Ｐ汀?/p>

本文對水分子的SPC、SPCE、TIP3P、TIP4P、TIP5P等勢能模型進行計算，考察水分子的熱力學(xué)性質(zhì)、徑向分布函數(shù)、自擴散系數(shù)，優(yōu)選最適合于水分子動力學(xué)性質(zhì)研究的勢能模型。

1? 勢能模型

SPC、SPCE、TIP3P、TIP4P、TIP5P這幾種模型都是剛體固定點電荷模型，計算時選用非極性作用Lennard-Jones勢模，分布為有效電荷表征電子云。

勢能函的公式為

式（1）由長程靜電相互作用和短程Lennard-Jones作用兩個部分構(gòu)成。為分子對作用勢能，、分別表示為、原子對、作循環(huán)，為原子所帶電荷，為原子所帶電荷，為原子間距離，為兩個分子的氧原子間作用距離，、為氧原子Lennard-Jones作用參數(shù)。此外，設(shè)為O-H鍵鍵長，為H-O-H鍵角，代表玻爾茲曼常數(shù)，是熱力學(xué)溫度單位。具體參數(shù)[9-10]如表1所示。

勢能的截斷，會使計算時能量在邊界處發(fā)生斷裂，導(dǎo)致哈密頓函數(shù)不守恒。采用上移的勢能函數(shù)[7]來彌補缺陷，計算公式為

其中，表示上移的位能，表示勢能的截斷半徑。

2? 模擬方法

模擬的水分子總數(shù)為300個，模擬溫度為298 K、壓強為1 atm（常溫常壓），模擬系綜（系綜表示大量宏觀上完全相同的體系的抽象集合）為正則系綜（NVT），水的起始密度設(shè)為1.0 g/cm3。

計算模型設(shè)為面心立方晶格，計算開始時水分子取向隨機，各分子的起始速度按Maxwell-Boltzmann分布取樣，對溫度的控制采用速度標定法，以確保體系總動量為零。

計算過程中采用最小鏡像準則和兩體有效勢基本近似，周期性邊界條件采用立方結(jié)構(gòu)[11]，計算盒子尺寸由密度確定。采用球形截去法，設(shè)置截斷半徑為0.9 nm。采用長程校正法校正截斷半徑之外的分子間相互作用能，長程靜電相互作用采用Ewald加和法[12]進行計算。

幾何構(gòu)型的固定采用SHAKE算法，采用5階Gear預(yù)估—校正法[13]求解運動方程。計算時所取時間步長為2 fs，每次模擬總時間為300 ps，前200 ps使體系達到平衡，后100 ps用于進行統(tǒng)計計算各種物理性質(zhì)。

3? 結(jié)果與討論

3.1? 熱力學(xué)性質(zhì)

計算常溫常壓下，不同勢能模型下水的內(nèi)能、汽化焓，并把所得結(jié)果與實驗值和其他文獻值進行比較。

汽化焓計算公式為

其中，H表示焓值，E代表體系內(nèi)能，P是壓強，V是體積，T是溫度，R是普適氣體常量。

模擬結(jié)果與比較如表2所示，其中MC代表量子力學(xué)從頭計算模型。選用SPCE模型模擬得到的水的內(nèi)能、汽化焓與實驗值和其他文獻值較接近，而其他模型可能受到了模型固有性質(zhì)和系綜的影響，模擬結(jié)果不太理想。總之，采用SPCE模型計算水分子的熱力學(xué)性質(zhì)最為合適。

3.2? 徑向分布函數(shù)

徑向分布函數(shù)是反映流體微觀結(jié)構(gòu)的物理量，它表明了與某個原子距離為的單位體積元內(nèi)出現(xiàn)的分子數(shù)，表達式為

其中，為體系的體積，為分子數(shù)，為兩個分子的質(zhì)心的距離，為函數(shù)。

計算常溫常壓下水分子在不同勢能模型下的gOO（O-O徑向分布）、gHO（H-O徑向分布）、gHH（H-H徑向分布），并與實驗值[12]進行比較。

圖1a所示為O-O徑向分布函數(shù)與實驗值?？梢钥闯觯煌瑒菽苣Ｐ偷腛-O徑向分布函數(shù)與實驗值的峰值變化趨勢都相符，其中SPCE、TIP4P和TIP5P三個模型與實驗值一樣出現(xiàn)了三個峰值，SPC、TIP3P只有一個峰值。綜合來看，在研究水的O-O徑向分布函數(shù)時，選用TIP5P模型更為合適。此外，SPCE、TIP4P和TIP5P模型的第二峰值出現(xiàn)在0.45 nm，第三峰值出現(xiàn)在0.68 nm，這與文獻[15]使用蒙特卡羅方法所得結(jié)果相一致。

圖1b所示為H-O徑向分布函數(shù)與實驗值。比較不同模型的峰值出現(xiàn)的位置，出現(xiàn)第一峰值和第二峰值的位置與實驗值相吻合，第一峰值的出現(xiàn)是氫原子與氧原子形成的共價鍵產(chǎn)生的結(jié)果。

圖1c所示為H-H徑向分布函數(shù)與實驗值。觀察第一峰值的峰位和第二峰值的峰位，其分別反映的是中心水分子與近臨水分子、第二配位圈水分子間的H-H距離。SPC和TIP3P模型的H-H徑向分布函數(shù)與實驗值吻合得不太理想，其他模型的模擬均與實驗值吻合的較好。

綜合分析徑向分布函數(shù)，可以發(fā)現(xiàn)隨著r的增加，函數(shù)趨于緩和，這種現(xiàn)象反映了水分子短程有序、長程無序的微觀結(jié)構(gòu)特征。綜上所述，選用SPCE、TIP4P和TIP5P三種勢能模型可以較好地反映出水分子的微觀結(jié)構(gòu)。

3.3? 自擴散系數(shù)

可以用兩種方法計算自擴散系數(shù)，一種是對均方位移（Mean Square Displacement，簡稱MSD）求斜率的Einstein方法[12]，另一種是對速度自相關(guān)函數(shù)求積分的Green-Kubo方法。Einstein方法的公式為

其中，為水分子的質(zhì)量，為玻爾茲曼常數(shù)，為體系的溫度，、分別為0時刻和時刻粒子的速度。

采用以上兩種方法分別計算水的自擴散系數(shù)。圖2所示采用Einstein方法為不同模型計算的均方位移?？梢钥闯鼍轿灰菩甭视尚〉酱蠓謩e對應(yīng)SPCE、TIP5P、TIP4P、SPC和TIP3P模型，表明自擴散系數(shù)也依次增大。

表3給出了采用兩種方法計算的自擴散系數(shù)與其他文獻[17]所得結(jié)果及實驗值[18]的比較?？梢钥闯?，采用Einstein方法計算的自擴散系數(shù)要比Green-Kubo方法更為理想，所得結(jié)果比文獻[17]的值更符合實驗值。采用Einstein方法計算所得結(jié)果比Green-Kubo方法小，這與速度自相關(guān)函數(shù)的積分區(qū)間選取有關(guān)，從而導(dǎo)致了松弛時間的差異，即當(dāng)速度自相關(guān)函數(shù)從1突降到0附近后，一直在0附近保持振蕩，這表明此后粒子的速度與時間不再相關(guān)。松弛時間一般很難準確確定，這將會影響到自擴散系數(shù)積分區(qū)間的選取。SPCE模型得到的結(jié)果與實驗值最接近，相對誤差僅為0.39%。TIP5P模擬結(jié)果次之，相對誤差為13.87%，但比SPC、TIP3P、TIP4P模擬結(jié)果要好得多。在TIP3P模型中，兩種方法得到的結(jié)果與實驗值相差最大，分析其原因可能是此模型下氫氧有效電荷數(shù)相差較大，使得氧原子束縛了氫原子。SPC和TIP4P計算結(jié)果也與實驗值相差較大，但是相比文獻[17]要好得多。

4? 結(jié)束語

勢能模型的選取是分子動力學(xué)計算的關(guān)鍵。本文采用分子動力學(xué)方法計算了五個勢能模型下水分子的內(nèi)能、汽化焓、徑向分布函數(shù)及自擴散系數(shù)。

綜合分析發(fā)現(xiàn)，SPCE模型最適用于水分子動力學(xué)性質(zhì)的研究，能夠較真實地反映水分子的熱力學(xué)性質(zhì)、微觀結(jié)構(gòu)和自擴散特性，值得學(xué)者廣泛關(guān)注。

基金項目

教育后勤疫情防控專項課題（課題編號：ZXYBKT202022）

參考文獻

[1] Alder B J， Wainwright T E. Phase Transition for a Hard Sphere System[J]. Journal of Chemical Physics， 1957， 27（5）： 1208-1209.

[2] Rahman A. Correlations in the Motion of Atoms in Liquid Argon[J]. Physical Review， 1964， 136（2）： 405-411.

[3] Cheung P S Y， Powles J G. The properties of liquid nitrogen. IV. A computer simulation[J]. Molecular Physics， 1975， 30： 921-949.

[4] Matsuoka O， Clementi E， Yoshimine M. CI study of the water dimer potential surface[J]. Journal of Chemical Physics， 1976， 64（4）： 1351-1361.

[5] Stilliger F H， Rahman A. Improved simulation of liquid water by molecular dynamics[J]. Journal of Chemical Physics， 1974， 60（4）： 1545-1557.

[6] Berendsen H J C， Postma J P M， Gunsteren W F， et al. Intermolecular Forces[M]. Holland： D. Reidel Publishing Company， 1981： 331.

[7] Berendsen H J C， Grigena J R， Straatsma T P. The missing term in effective pair potentials[J]. J Journal of Physical Chemistry， 1987， 91（24）： 6269-6271.

[8] Mahoney M W， Jorgensen W L. A five-five model for liquid water and the reproduction of the density anomaly by rigid， nonpolarizable potential functions[J]. Journal of Chemical Physics， 2000， 112（20）： 8910-8922.

[9] 李印實， 何雅玲， 孫杰， 等. 不同水分子模型凝結(jié)系數(shù)的分子動力學(xué)對比研究[J]. 西安交通大學(xué)學(xué)報， 2006， 40（11）： 1272-1275.

[10] Mills R. Self-diffusion in normal and heavy water in the range 1-45.deg[J]. Journal of Physical Chemistry， 1973， 77（5）： 685-688.

[11] Allen M P， Tildesley D J. Computer Simulation of liquids[M]. Oxford： Clareendon Press， 1987.

[12] De Leeuw S W， Perram J W， Smith E R. Simulation of electrostatic systems in periodic boundary conditions. I. Lattice sums and dielectric constant[C]// Proceedings of the Royal Society A， 1980， 373： 27-56.

[13] Gear C W. Numerical Integration of Ordinary Differential Equations[M]. New Jersey： Prentice Hill， 1971.

[14] Soper A K， Phillips M G. A new determination of the structure of water at 25℃[J]. Chemical Physics， 1986， 107（11）： 47-60.

[15] 吳熊武， 騰藤， 李以圭， 等. 水—甲醇體系的Monte Carlo分子動力學(xué)模擬[J]. 化學(xué)學(xué)報， 1992， 50： 543-548.

[16] Krynicki K， Green C D， Sawyer D W. Pressure and temperature dependence of self-diffusion in water[J]. Faraday Discussions of the Chemical Society， 1978， 66： 199-208.

[17] Jorgensen W L， Jenson C. Temperature Dependence of TIP3P， SPC， and TIP4P Water from NPT Monte Carlo Simulations： Seeking Temperatures of Maximum Density[J]. Journal of Computational Chemistry， 1998， 19（10）： 1179-1186.

[18] 劉娟芳， 曾丹苓， 蔡智勇， 等. 擴散系數(shù)的分子動力學(xué)模擬[M]. 工程熱物理學(xué)報， 2006， 27（3）： 373-375.

作者簡介：

俞聯(lián)夢（1983—），通信作者，女，云南大理人，碩士研究生，講師，從事高校教學(xué)工作。主要研究方向：分子動力學(xué)。

E-mail： 175365186@qq.com

（收稿日期：2020-06-18）