基于TraPPE-UA勢能模型的超高壓乙醇分子動力學模擬

周曉平，周 偉，郝江濤

（鄭州大學西亞斯國際學院，河南 鄭州 451150）

0 引 言

分子動力學模擬是依靠牛頓力學來模擬分子體系運動的一套分子模擬方法，不僅可以用于模擬物質(zhì)的宏觀性質(zhì)，還可以得出微觀結(jié)構(gòu)、粒子運動等信息[1]。一般認為壓強超過100MPa就是超高壓，物質(zhì)在高壓作用下，其物理和化學性質(zhì)會發(fā)生巨大的變化，在超高壓作用下，液體的壓縮性顯著增大[2-3]。高壓物理實驗面臨著實驗條件、費用、實驗設(shè)備的限制，而分子動力學計算利用計算機模擬，可以很好地解決這些問題。乙醇分子是常見的有機、極性分子，常用于化工生產(chǎn)中，研究其熱力學性質(zhì)和動力學性質(zhì)對工業(yè)生產(chǎn)有重要意義。且乙醇和水在很多性質(zhì)上都存在相似性，被認為是水的理想替代物[4]，關(guān)于乙醇性質(zhì)的研究也一直是熱點。張陽[5]用Monte Carlo方法研究了超臨界乙醇的結(jié)構(gòu)性質(zhì)和氫鍵，壓強范圍從常壓到高壓70 MPa；Saiz等[6]用分子動力學方法模擬研究了不同熱力學狀態(tài)下液態(tài)乙醇的結(jié)構(gòu)、氫鍵和擴散等動態(tài)性質(zhì)；Markus等[7]研究了溫度和壓強（0.1～35 MPa）對甲醇和乙醇體系中氫鍵的影響。目前對乙醇的研究壓強一般都在100MPa以內(nèi)，更高壓強下乙醇性質(zhì)的文章還少見報道。因此本文對液態(tài)乙醇的研究進行了補充，采用平衡分子動力學方法，模擬溫度298～500K、壓強從10MPa到超高壓300 MPa條件下不同狀態(tài)點乙醇的熱力學性質(zhì)、結(jié)構(gòu)性質(zhì)和動力學性質(zhì)。本文的乙醇分子動力學模擬值可為其傳質(zhì)性質(zhì)應用提供理論參考。

1 模型建立

1.1 勢能模型

勢能模型的選取對模擬結(jié)果的準確與否有著決定性影響。目前文獻中分子模擬常用的乙醇模型主要是全原子模型和聯(lián)合原子模型。鑒于全原子模型所需的計算時間較長，因此大部分文獻采用聯(lián)合原子模型進行模擬。聯(lián)合原子模型有OPLS-UA（optimized potentials for liquid simulations-united atom）模型和TraPPE-UA（transferable potentials for phase equilibria-united atom）模型兩種。OPLS-UA模型由Jorgensen[8]開發(fā)，主要針對液相有機小分子，能夠很好地模擬液態(tài)乙醇的各種動態(tài)性質(zhì)。但是，在Chalaris和Samios[9-10]用OPLS-UA模型對超臨界甲醇進行分子動力學模擬的研究中，壓強模擬值的準確性很差，說明OPLS-UA勢能模型在壓強較高時并不能真實反映分子的動態(tài)性質(zhì)。TraPPE-UA模型是由Siepmann[11]小組開發(fā)，由于其模型參數(shù)已經(jīng)用實驗的相平衡數(shù)據(jù)優(yōu)化過，因此在很大的壓強范圍內(nèi)都有較好的準確度。本文采用TraPPE-UA勢能模型對高壓下乙醇分子團簇的性質(zhì)進行研究。

在TraPPE-UA模型中，乙醇分子由CH3、CH2基團和O原子、H原子合計4個作用位點構(gòu)成。乙醇分子間的勢能計算式為

式中：qi、qj——第i、j個原子或基團上所帶電荷；

rij——原子或離子i與j間的距離；

σij、εij——原子間L-J作用參數(shù)。

式（1）右側(cè)第 1 項表示短程 Lennard-Jones（L-J）勢，第2項表示長程庫侖勢。TraPPE-UA模型和OPLS-UA模型的相關(guān)參數(shù)見表1。

表1 TraPPE-UA模型和OPLS-UA模型參數(shù)列表

1.2 實驗方法

本研究采用Materials Explorer軟件，在溫度為298，400，500 K，壓強分別為 10，100，300 MPa 條件下對乙醇分子體系進行了分子動力學模擬。采用Nose-Hoover熱浴法控制溫度，完成不同壓強條件下的等溫模擬[13]。分子動力學模擬需進行以下假設(shè)：1）兩體有效勢能近似；2）最小鏡像準則[14]。分子動力學模擬的起始構(gòu)型為面心立方晶格，將乙醇分子的起始取向定為隨機分布，采用立方周期性邊界條件。模擬過程中所用的乙醇分子總數(shù)為100個。采用球形截去法進行位能截斷[14]，采用Ewald加和法求解長程靜電相互作用力[15]，采用Shake算法固定乙醇分子的幾何構(gòu)型[16]。分子動力學模擬過程中，每個分子都可視為是剛性體，運用五階Gear預測-校正法求解體系的運動方程[17]。分子動力學模擬過程中，時間步長取10-15s，單次模擬總時間取3×10-10s，采用NVT系綜進行10-10s的初始化，使得體系達到相對平衡，隨后在NPT系綜中完成2×10-10s的模擬，完成各種參數(shù)的計算，模擬過程中每隔10-14s就輸出一次構(gòu)型，用于結(jié)果分析。

2 結(jié)果與分析

2.1 熱力學性質(zhì)

焓是物質(zhì)重要的熱力學狀態(tài)量，體系焓計算式為

式中：E——體系的內(nèi)能，kJ；

P——體系的壓強，MPa；

V——體系的體積，m3。

表2是本文模擬溫度為298，400，500K，壓強為10，100，300 MPa 時乙醇的焓值，并與文獻[18]的模擬值和文獻[19]用狀態(tài)方程推導的焓值進行比較，結(jié)果較為接近，可見本文模擬結(jié)果可靠。

表2 乙醇在不同狀態(tài)點的焓值

從表中可以看出，當壓強相同時，隨著溫度的升高乙醇體系的焓值逐漸增大。這是由于當溫度升高，乙醇分子間的熱運動加劇，體系內(nèi)能增大，焓也跟著增大。相同溫度下，本文模擬的焓值隨壓強的增大而增大，但文獻[18-19]中焓隨壓強的變化關(guān)系不明確，由于現(xiàn)有文獻中尚未見到高壓下乙醇焓值實驗數(shù)據(jù)的相關(guān)報道，此計算結(jié)果有待進一步驗證。

2.2 結(jié)構(gòu)性質(zhì)

徑向分布函數(shù)（RDF)是反映物質(zhì)微觀結(jié)構(gòu)的重要物理量，主要用于描述距離中心粒子為r處出現(xiàn)另一粒子的概率密度與隨機分布概率密度的比值[13]，其計算式為

式中：ρ——系統(tǒng)密度；

dN——距離中心分子（r，r+dr）區(qū)間內(nèi)的分子數(shù)目。

圖1為T=298 K時乙醇分子的氧氫徑向分布函數(shù)隨壓強的變化。圖2為P=10MPa時氧氫徑向分布函數(shù)隨溫度的變化。

從圖1可看出，壓強變化對OH徑向分布函數(shù)第1峰影響較大，隨著壓強的增大，OH徑向分布函數(shù)第1峰分別出現(xiàn)在1.69，1.67，1.65?，峰位逐漸左移，峰值依次減小，第2峰也下降，但不如第1峰明顯，在Petravic[20]和Dellis[21]等的研究中，也發(fā)現(xiàn)了OH徑向分布函數(shù)的第1峰隨壓強增大而減小的現(xiàn)象。計算結(jié)果表明，隨著壓強的增加，第一配位圈內(nèi)乙醇分子間的氫原子與氧原子間距逐漸減小，氫鍵作用逐漸增強。

圖1 298K時乙醇的氧氫徑向分布函數(shù)

圖2 10MPa時乙醇的氧氫徑向分布函數(shù)

從圖2可看出，當壓強相同時，隨著溫度的升高，OH徑向分布函數(shù)第1峰和第2峰逐漸下降，峰谷抬高，峰位右移。計算結(jié)果表明隨著溫度的升高，乙醇分子間的氫原子與氧原子間距逐漸增大，氫鍵作用逐漸減弱。表3為各狀態(tài)點乙醇的OH、OO、HH徑向分布函數(shù)第1特征峰值與峰位。隨著溫度的升高和壓強的增大，OH、OO、HH徑向分布函數(shù)第1峰峰值均下降，但OH徑向分布函數(shù)的變化較OO、HH徑向分布函數(shù)更為明顯。

表3 不同狀態(tài)下乙醇的OH、OO、HH徑向分布函數(shù)第1特征峰值與峰位

圖3、圖4為100MPa時OO、HH徑向分布函數(shù)隨溫度的變化。從圖中可看出，相同壓強下，隨著溫度的升高，OO、HH徑向分布函數(shù)第1峰下降，峰位右移，峰谷抬高，第2峰也下降，且在500K時，第2峰消失。這說明隨著溫度的升高，乙醇的長程有序程度逐漸下降。這與高壓下水的結(jié)構(gòu)隨溫度和壓強變化的趨勢相同[22]。

圖3 100MPa時乙醇的氧氧徑向分布函數(shù)

圖4 100MPa時乙醇的氫氫徑向分布函數(shù)

圖5～圖8為乙醇的MeH、MeO、MeMe徑向分布函數(shù)。隨著溫度的升高和壓強的增大，各徑向分布函數(shù)第1峰變化最明顯，峰高下降，峰谷抬高。500 K時，MeH徑向分布函數(shù)的第2、3、4峰峰高已超過第1峰。MeMe徑向分布函數(shù)第2峰逐漸消失。

通過對比不同原子和基團間徑向分布函數(shù)的第一波谷位置rmin1可發(fā)現(xiàn)，rOH＜rHH＜rOO＜rMeH＜rMeO＜rMeMe。 計算結(jié)果表明：相同溫度和壓強下，第一配位圈內(nèi)親水的O原子與H原子距離最小，而疏水的烷基之間的距離最遠。

2.3 動力學性質(zhì)

圖5 298K時乙醇的甲基氫徑向分布函數(shù)

圖6 100MPa時乙醇的甲基氧徑向分布函數(shù)

圖7 400K時乙醇的甲基甲基徑向分布函數(shù)

結(jié)合國內(nèi)外研究現(xiàn)狀，自擴散系數(shù)D的計算方法主要包括Einstein法和Green-Kubo法。Eintein法通過對均方位移（MSD）求斜率得到擴散系數(shù)。鑒于分子動力系統(tǒng)中的原子在永不停息地做無規(guī)則運動，各原子在不同時刻的位置都不相同。將r→（t）、r→（0）作為粒子在t時刻和初始時刻的位置，則粒子均方位移計算式[13]為

圖8 300MPa時乙醇的甲基甲基徑向分布函數(shù)

式中＜＞表示平均值。

根據(jù)統(tǒng)計學原理，只要分子數(shù)目足夠多，計算時間足夠長，系統(tǒng)的任一瞬間都可以當作時間的零點，所計算的平均值應該相同。根據(jù)愛因斯坦的擴散定律，擴散系數(shù)（diffusion constant）與均方位移的關(guān)系是：

式中D即為粒子的擴散系數(shù)。

當模擬的時間很長時，均方位移對時間曲線的斜率的1/6就是擴散系數(shù)。表4是本文用均方位移計算的乙醇的自擴散系數(shù)。

表4 乙醇在不同狀態(tài)下的自擴散系數(shù)

由表可知，在相同壓強下，隨溫度的升高，乙醇的自擴散系數(shù)增大。這與文獻[23]測定的常壓下乙醇自擴散系數(shù)隨溫度的變化趨勢相同，且隨著壓強的增大，這種變化趨勢逐漸減緩。而相同溫度下，隨著壓強的增加，乙醇自擴散系數(shù)逐漸減小。這是因為壓強增大使乙醇分子之間的距離變小，分子間排列更為緊密，相互碰撞現(xiàn)象明顯增加，從而阻礙了分子的擴散運動；當溫度不斷上升時，乙醇分子間的距離逐漸增加，相互作用力減弱，此時將有利于分子的擴散運動。本文模擬值與文獻[18]的模擬結(jié)果基本一致，但在500K，10MPa時相差較大，由于實驗條件的限制，高壓下乙醇分子團簇的自擴散系數(shù)實驗數(shù)據(jù)匱乏，有待進一步驗證。

3 結(jié)束語

本文模擬了溫度298～500K，壓強從高壓10MPa到超高壓300MPa條件下乙醇分子體系的熱力學性質(zhì)、結(jié)構(gòu)性質(zhì)和動力學性質(zhì)，并得出以下結(jié)論：

1）相同壓強下，乙醇體系的焓值隨溫度的升高而增大，本文模擬值與文獻模擬值接近，本文所用TraPPE-UA模型可靠。

2）通過對徑向分布函數(shù)的研究，發(fā)現(xiàn)在壓強很大的情況下，乙醇體系的結(jié)構(gòu)依然很有規(guī)律。隨著壓強和溫度的升高，各徑向分布函數(shù)的第1峰高度明顯下降，乙醇的長程有序程度下降。

3）當壓強相同時，乙醇的自擴散系數(shù)隨溫度的升高而增大，但壓強增大時，這種變化趨勢逐漸減緩。乙醇的自擴散系數(shù)隨溫度的升高而增大，而壓強對擴散系數(shù)的影響與溫度正好相反。

[1]RAHMAN A.Correlations in the motion of atoms in liquid argon[J].Phys Rev，1964（136）：405-411.

[2]曾勇平，朱曉敏，楊正華.水、甲醇和乙醇液體微結(jié)構(gòu)性質(zhì)的Car-Parrinello分子動力學模擬[J].物理化學學報，2011，27（12）：2779-2785.

[3]曾勇平，時榮，楊正華.Be～（2+）在水、甲醇和乙醇中結(jié)構(gòu)性質(zhì)的從頭算分子動力學模擬[J].物理化學學報，2013，2（10）：2180-2186.

[4]羅輝，王睿，范維玉，等.分子動力學模擬計算超臨界CO的溶解度參數(shù)[J].石油學報（石油加工），2015，31（1）：78-83.

[5]張陽.分子模擬在純超臨界流體及其二元混合物體系中的應用[D].北京：清華大學，2005.

[6]SAIZ L，PADRO J A，GUARDIA E.Structure and dyna mics of liquid ethanol[J].Phys Chem B，1997（101）：78-86.

[7]MARKUS M，HOFFMANN，MARK S C.Are there hydrogen bonds in supercritical methanol and ethanol[J].Phys Chem B，1998（102）：263-271.

[8]JORGENSEN W L.Optimized intermolecular potential functions for liquid alcohols[J].Phys Chem，1986（90）：1276-1284.

[9]CHALARIS M，SAMIOS J.Hydrogen bonding in supercritical methanol A molecular dynamics investigation[J].Phys Chem B，1999（103）：1161-1166.

[10]CHALARIS M，SAMIOS J.Translational and rotational dynamics in supercritical methanol from molecular dynamics simulation[J].Pure and Applied Chemistry，2004（76）：203-213.

[11]CHEN B，POTOFF J J，SIEPMANN J I.Monte carlo calculations for alcohols and their mixtures with alkanes.transferable potentials for phase equilibria.5.United-Atom description of primary secondary and tertiary alcohols[J].Phys Chem B，2001（105）：3093-3104.

[12]AND I Y G，PIOTROVSKAYA E M.Properties of coexisting phases for the ethanol-ethane binary system by computer simulation[J].Journal of Physical Chemistry B，1999，103（36）：7681-7686.

[13]陳正隆，徐為人，湯立達.分子模擬的理論與實踐[M].北京：化學工業(yè)出版社，2007：6.

[14]ALLEN M P，TILDESLEY D J.Computer simulation of liquids[M].Oxford： Clareendon Press，1987：327-341.

[15]DE-LEEUW S W，PERRAM J W，SMITH E R.Simulation ofelectrostatic systemsin periodic boundary conditions I Lattice sums and dielectric constant[J].Proc R Soc Lond A，1980（373）：27-56.

[16]RYCHAERT J P，CICCOTTI G，BERENDSEN H J C.Numercial intergration of the cartesian equations of motion of a system with constraints:molecular dynamics of n-alkanes[J].J Chem Phys，1977（23）：327-341.

[17]GEAR C W.Numerical integration of ordinary differential equations[J].Mathemetics of Computation，1967，21（98）：146-156.

[18]李勇，劉錦超，蘆鵬飛，等.從常溫常壓到超臨界乙醇的分子動力學模擬[J].物理學報，2010，59（7）：4880-4887.

[19]DILLON H E，PENONCELLO S G.A fundamental equation for calculation of the thermodynamic properties ofethanol[J].International Journal of Thermophysics，2004，25（2）：321-335.

[20]PETRAVIC J，DELHOMMELLE J.Influence of temperature，pressure and internal degrees of freedom on hydrogen bonding and diffusion in liquid ethanol[J].Chemical Physics，2003，286（2/3）：303-314.

[21]DELLIS D，MICHALIS C A，SAMIOS J.Pressure and temperature dependence of the hydrogen bonding in supercriticalethanol： A computersimulation study[J].Journal of Physical Chemistry B，2005，109（39）：18575-18590.

[22]周曉平，楊向東，劉錦超.超高壓水的分子動力學模擬[J].高壓物理學報，2009，23（4）：1-5.

[23]KARGER N，VARDAG T.Temperature dependence of self iffusion in compressed monohydric alcohols[J].Journal of Chemical Physics，1990，93（93）：3437-3444.