賈志剛,余 敏,張素英
Nose-Poincare恒溫分子動力學(xué)模擬中廣義質(zhì)量Q的研究
賈志剛1,余 敏2,張素英1
(1.山西大學(xué)物理電子工程學(xué)院理論物理研究所,山西太原030006;2.海軍蚌埠士官學(xué)校基礎(chǔ)部,安徽蚌埠233012)
研究了恒溫分子動力學(xué)模擬中系統(tǒng)能量的演化問題,數(shù)值結(jié)果顯示廣義質(zhì)量Q對應(yīng)了系統(tǒng)對外界作用的敏感程度,是一種系統(tǒng)慣量.Q值減小,系統(tǒng)慣量就減小,系統(tǒng)達(dá)到平衡狀態(tài)所需的時間就更長,系統(tǒng)能量的相對波動也更大,外界對系統(tǒng)的作用更明顯.
Nose-Poincare哈密爾頓量;恒溫分子動力學(xué);廣義質(zhì)量
分子動力學(xué)模擬經(jīng)過多年的發(fā)展,已經(jīng)成為一種重要的研究工具,在物理學(xué)、化學(xué)及生物學(xué)研究中得到了廣泛的應(yīng)用.傳統(tǒng)分子動力學(xué)模擬的是微正則(NVE)系綜,但是在實(shí)際應(yīng)用當(dāng)中,恒溫或者恒壓的實(shí)驗(yàn)更為容易實(shí)現(xiàn).于是,Anderson提出了用一種擴(kuò)展Hamilton量來研究恒壓(NPH)分子動力學(xué)[1].而后,Nose又提出了Nose恒溫箱,用來模擬恒溫的正則系綜[2].此后,擴(kuò)展Hamilton量方法逐漸發(fā)展成為一種重要的分子動力學(xué)模擬方法,被應(yīng)用于現(xiàn)代納米技術(shù)及非平衡系統(tǒng)的研究[3,4].各種不同的方法相互融合,產(chǎn)生了適用于恒溫恒壓(NPT)系綜的擴(kuò)展Hamilton量[5].另一方面,Nose恒溫箱也進(jìn)一步發(fā)展,產(chǎn)生了Nose-Hoover恒溫箱[6]和Nose-Poincare恒溫箱[7].本文,我們所做的模擬就是基于Nose提出來的擴(kuò)展的Hamilton量[2]:
以及Laird等人對其所做的Poincare變換[7]:
式中,s是人為引入的時間標(biāo)度因子,我們可以將其看作一個廣義的正則坐標(biāo),ps是與s共軛的正則動量,mi代表分子質(zhì)量,pi是經(jīng)過標(biāo)度的各個分子的動量,qij表示分子間的距離,g是系統(tǒng)的自由度,k是波爾茲曼常量,T代表外界溫度,H0代表H在初始時刻的值.特別值得注意的是:參數(shù)Q代表廣義質(zhì)量,本文主要研究Q在系統(tǒng)演化過程中所起的重要作用,并對其進(jìn)行詳細(xì)討論.
擴(kuò)展Hamilton量的方法已經(jīng)成為分子動力學(xué)模擬的一種標(biāo)準(zhǔn)方法,在不同的系綜下,可以構(gòu)造不同的Hamilton量來模擬實(shí)際的實(shí)驗(yàn).擴(kuò)展Hamilton量的基本思想就是:引入一個標(biāo)度因子,重新標(biāo)度時間或者空間,并將引入的標(biāo)度因子作為一個新的廣義坐標(biāo),并引入與其相對應(yīng)的新的廣義動量和廣義質(zhì)量,從而構(gòu)造出一個新的Hamilton量.
本文從擴(kuò)展的Nose-Poincare Hamilton量(2)出發(fā)進(jìn)行模擬,研究其中引入的廣義質(zhì)量Q對系統(tǒng)演化過程的影響.這里,要注意一點(diǎn),為了保證這個系統(tǒng)能夠產(chǎn)生正則分布平均值,式(1)中最后一項(xiàng)的g+1在(2)中變成了g.可以證明:在等概率假設(shè)和各態(tài)遍歷假設(shè)的前提下,以上系統(tǒng)會產(chǎn)生正則分布平均值.
物理量A(p,q)的統(tǒng)計平均值[1]為:
由擴(kuò)展的Hamilton量(2),可以寫出運(yùn)動方程:
下一節(jié),我們將對運(yùn)動方程進(jìn)行離散,數(shù)值模擬其演化過程.
在分子動力學(xué)模擬中,常見的算法主要有四種[8]:Verlet算法,速度Verlet算法,leap-frog算法,預(yù)校正算法(詳見附錄).這里,我們采用修正的leap-frog算法,這種算法是一種辛算法,能夠保持系統(tǒng)長時間演化的穩(wěn)定性.下面就是經(jīng)過離散化的運(yùn)動方程:
具體地,本文選擇Lennard-Jones流體進(jìn)行數(shù)值模擬,Lennard-Jones勢的表達(dá)式如下:
其中,ε為勢阱深度,σ為平衡長度,即排斥勢與吸引勢相等的分子間距.初始時刻,我們將四千個粒子放在簡單立方晶格中,關(guān)于勢能的計算,采用周期邊界條件和最小映像原理,而截斷長度的選取則是一項(xiàng)很有技巧性的工作,對于不同系統(tǒng),截斷長度的選取方法與選取的截斷長度是不同的,這里我們所選模擬的是一個典型的Lennard-Jones流體,其最合適的截斷長度為2.5σ[7,8].計算過程以及計算結(jié)果的給出,均運(yùn)用對比單位,對于Lennard-Jones流體,選擇m,ε,σ作為質(zhì)量、能量、長度的基本單位,溫度在對比單位下就是T′=由于σ為長度的基本單位,這時截斷長度即2.5.
各物理量的統(tǒng)計平均值由(3)給出.特別地,能量關(guān)于pi,qi的表達(dá)式由(1)給出.我們計算了Q從0.5到100不同取值下系統(tǒng)的演化情況,這里,我們只列出Q=0.5,1,10,100這四種典型的情況來加以研究.
圖1(P221)描述了在四種不同的情況下,系統(tǒng)能量的演化過程,從中可以看出,Q值越大,達(dá)到平衡所經(jīng)過的時間越短,并且達(dá)到平衡后,系統(tǒng)能量的演化更為平穩(wěn).另外,從中可以看出,Q取值越大,系統(tǒng)能量的波動越小,這一點(diǎn)在圖2(P221)中可以更清楚地看出.
圖2所描述的是系統(tǒng)能量的相對波動,即系統(tǒng)能量波動與系統(tǒng)能量平均值的比值:
圖1 能量演化過程Fig.1 Evolvement of energy
(8)式也可以稱作能量的相對誤差,但這種誤差不是由于算法的原因而產(chǎn)生的,而是由于系統(tǒng)本身與外界相互作用而形成的一種物理現(xiàn)象.從圖中可以看出,Q的值越大,相對波動越小,系統(tǒng)的能量更不容易偏離平均值.也就是說,Q值較大的系統(tǒng)具有更強(qiáng)的穩(wěn)定性,不會由于與外界的能量交換而產(chǎn)生較大波動.
圖2 能量相對波動Fig.2 Relative fluctuation of energy
圖3 (P222)中,六角星圖線是系統(tǒng)達(dá)到平衡時的動量分布,實(shí)線表示經(jīng)典的波爾茲曼分布曲線,可以看出,模擬所得的結(jié)果與理論基本符合,也從另一個側(cè)面說明模擬的有效性與正確性.
本文討論了在不同的Q值情況下,系統(tǒng)的演化過程和相對能量波動,并考察了不同Q值的平衡狀態(tài)下動量分布情況.從中我們可以看出:廣義質(zhì)量Q較大時,(1)系統(tǒng)總能量演化更為平穩(wěn),受外界影響較小;(2)系統(tǒng)達(dá)到平衡狀態(tài)后,系統(tǒng)能量相對波動的幅度更小.這就說明了,Q的取值越大,對應(yīng)的系統(tǒng)漲落越小,相反,Q取值越小,系統(tǒng)漲落越大.
圖3 動量分布Fig.3 The distribution of momentum
綜上所述,我們可以得出:同Nose-Hoover恒溫分子動力學(xué)一樣,在Nose-Poincare恒溫分子動力學(xué)中, Q對應(yīng)于系統(tǒng)的慣量,Q取值越大,系統(tǒng)慣量越大,外部作用所產(chǎn)生的波動就越小.反之,Q越小,系統(tǒng)慣量就越小,外部作用所產(chǎn)生的波動就越大,系統(tǒng)對外界的作用越敏感.
附錄:四種常見的分子動力學(xué)算法
在分子動力學(xué)模擬中,算法的好壞直接決定了模擬過程的計算效率與模擬結(jié)果的有效性.這里我們介紹四種常用的算法:Verlet算法,速度Verlet算法,Leap-frog算法,預(yù)校正算法.
Verlet算法:
速度Verlet算法:
Leap-frog算法:
預(yù)校正算法:
從(9)-(12)式可以看出,Verlet算能夠保持系統(tǒng)的時間反演不變性,leap-frog算法是一種辛格式,但是位置與速度不同步,而速度Verlet算法與預(yù)校正算法位置與速度同步,但是卻破壞了時間反演不變性.本文中,我們采用一種修正的leap-frog算法構(gòu)造了一種辛格式[3],不僅能夠保持系統(tǒng)長時間演化的穩(wěn)定性,同時也實(shí)現(xiàn)了位置與速度的同步計算,是一種有效的分子動力學(xué)算法.
[1] ANDERSEN H C.Molecular Dynamics Simulations at Constant Pressure And/or Temperature[J].J Chem Phys,1980,72 (4):2384-2393.
[2] NOSE S.A Molecular Dynamics Method for Simulations in the Canonical Ensemble[J].Molecular Physics,2002,100(1): 191-198.
[3] AKIMOV A V,NEMU KHIN A V,MOSKOVSKY A A.Molecular Dynamics of Surface-Moving Thermally Driven Nanocars[J].J Chem Theory Comput,2008,4(4):652-656.
[4] CUENDET Michel A.The Jarzynski Identity Derived from General Hamiltonian or Non-Hamiltonian Dynamics Reproducing NVT or NPT Ensembles[J].The J ournal of Chemical Physics,2006,125(14):144109.
[5] CUI Zhi-wei,SUN Yi,QU Jian-min.Constant Pressure Molecular Dynamics Simulation for Ionic System[J].Journal of Computational and Theoretical N anoscience,2008,5(8):1646-1650.
[6] HOOVER W G.Canonical Dynamics:Equilibrium Phase-space Distributions[J].Phys Rev A,1985,31(4):1695-1697.
[7] BOND S D,LEIMKU HL ER B J,LAIRD B B.The Nose-poincare Method for Constant Temperature Molecular Dynamics [J].J Comput Phys,1999,151(1):114-134.
[8] FRENKEL&SMIT.分子模擬——從算法到應(yīng)用[M].汪文川,等譯.北京:化學(xué)工業(yè)出版社,2002.
Generalized Mass Q in Nose-Poincare Constant Temperature Molecular Dynamics
J IA Zhi-gang1,YU Min2,ZHAN G Su-ying1
(1.School ofPhysics and Electronics Engineering Institute ofTheoretical Physics,Shanxi University,Taiyuan030006,China; 2.Basic Courses Department,Bengbu N aval Petty Of f icer Academy,Bengbu233012,China)
We simulate the evolvement of system enegy for constant temperature molecular dynamics,and conclude that the generalized mass Q is the inertia of system which determines the sensitivity of system to external influence.If Q decreases,the inertia of system decreases,it takes long to reach equilibrium,and the relative fluctuation of system energy increanses,and the system is more sensitive to external influence.
Nose-poincare Hamiltonian;constant temperature molecular dynamics;generalized mass
O411.3
A
0253-2395(2010)02-0219-05
2009-09-05;
2009-12-10
國家自然科學(xué)基金(10972125)
賈志剛(1983-),男,山西太原人,碩士生,研究方向:分子動力學(xué)模擬.E-mail:firefox_me@yahoo.cn