賈志剛,余 敏,張素英

Nose-Poincare恒溫分子動(dòng)力學(xué)模擬中廣義質(zhì)量Q的研究

賈志剛1,余 敏2,張素英1

(1.山西大學(xué)物理電子工程學(xué)院理論物理研究所,山西太原030006;2.海軍蚌埠士官學(xué)校基礎(chǔ)部,安徽蚌埠233012)

研究了恒溫分子動(dòng)力學(xué)模擬中系統(tǒng)能量的演化問(wèn)題,數(shù)值結(jié)果顯示廣義質(zhì)量Q對(duì)應(yīng)了系統(tǒng)對(duì)外界作用的敏感程度,是一種系統(tǒng)慣量.Q值減小,系統(tǒng)慣量就減小,系統(tǒng)達(dá)到平衡狀態(tài)所需的時(shí)間就更長(zhǎng),系統(tǒng)能量的相對(duì)波動(dòng)也更大,外界對(duì)系統(tǒng)的作用更明顯.

Nose-Poincare哈密爾頓量;恒溫分子動(dòng)力學(xué);廣義質(zhì)量

0 引言

分子動(dòng)力學(xué)模擬經(jīng)過(guò)多年的發(fā)展,已經(jīng)成為一種重要的研究工具,在物理學(xué)、化學(xué)及生物學(xué)研究中得到了廣泛的應(yīng)用.傳統(tǒng)分子動(dòng)力學(xué)模擬的是微正則(NVE)系綜,但是在實(shí)際應(yīng)用當(dāng)中,恒溫或者恒壓的實(shí)驗(yàn)更為容易實(shí)現(xiàn).于是,Anderson提出了用一種擴(kuò)展Hamilton量來(lái)研究恒壓(NPH)分子動(dòng)力學(xué)[1].而后,Nose又提出了Nose恒溫箱,用來(lái)模擬恒溫的正則系綜[2].此后,擴(kuò)展Hamilton量方法逐漸發(fā)展成為一種重要的分子動(dòng)力學(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提出來(lái)的擴(kuò)展的Hamilton量[2]:

以及Laird等人對(duì)其所做的Poincare變換[7]:

式中,s是人為引入的時(shí)間標(biāo)度因子,我們可以將其看作一個(gè)廣義的正則坐標(biāo),ps是與s共軛的正則動(dòng)量,mi代表分子質(zhì)量,pi是經(jīng)過(guò)標(biāo)度的各個(gè)分子的動(dòng)量,qij表示分子間的距離,g是系統(tǒng)的自由度,k是波爾茲曼常量,T代表外界溫度,H0代表H在初始時(shí)刻的值.特別值得注意的是:參數(shù)Q代表廣義質(zhì)量,本文主要研究Q在系統(tǒng)演化過(guò)程中所起的重要作用,并對(duì)其進(jìn)行詳細(xì)討論.

1 Nose-Poincare恒溫分子動(dòng)力學(xué)

擴(kuò)展Hamilton量的方法已經(jīng)成為分子動(dòng)力學(xué)模擬的一種標(biāo)準(zhǔn)方法,在不同的系綜下,可以構(gòu)造不同的Hamilton量來(lái)模擬實(shí)際的實(shí)驗(yàn).擴(kuò)展Hamilton量的基本思想就是:引入一個(gè)標(biāo)度因子,重新標(biāo)度時(shí)間或者空間,并將引入的標(biāo)度因子作為一個(gè)新的廣義坐標(biāo),并引入與其相對(duì)應(yīng)的新的廣義動(dòng)量和廣義質(zhì)量,從而構(gòu)造出一個(gè)新的Hamilton量.

本文從擴(kuò)展的Nose-Poincare Hamilton量(2)出發(fā)進(jìn)行模擬,研究其中引入的廣義質(zhì)量Q對(duì)系統(tǒng)演化過(guò)程的影響.這里,要注意一點(diǎn),為了保證這個(gè)系統(tǒng)能夠產(chǎn)生正則分布平均值,式(1)中最后一項(xiàng)的g+1在(2)中變成了g.可以證明:在等概率假設(shè)和各態(tài)遍歷假設(shè)的前提下,以上系統(tǒng)會(huì)產(chǎn)生正則分布平均值.

物理量A(p,q)的統(tǒng)計(jì)平均值[1]為:

由擴(kuò)展的Hamilton量(2),可以寫出運(yùn)動(dòng)方程:

下一節(jié),我們將對(duì)運(yùn)動(dòng)方程進(jìn)行離散,數(shù)值模擬其演化過(guò)程.

2 數(shù)值模擬

在分子動(dòng)力學(xué)模擬中,常見的算法主要有四種[8]:Verlet算法,速度Verlet算法,leap-frog算法,預(yù)校正算法(詳見附錄).這里,我們采用修正的leap-frog算法,這種算法是一種辛算法,能夠保持系統(tǒng)長(zhǎng)時(shí)間演化的穩(wěn)定性.下面就是經(jīng)過(guò)離散化的運(yùn)動(dòng)方程:

具體地,本文選擇Lennard-Jones流體進(jìn)行數(shù)值模擬,Lennard-Jones勢(shì)的表達(dá)式如下:

其中,ε為勢(shì)阱深度,σ為平衡長(zhǎng)度,即排斥勢(shì)與吸引勢(shì)相等的分子間距.初始時(shí)刻,我們將四千個(gè)粒子放在簡(jiǎn)單立方晶格中,關(guān)于勢(shì)能的計(jì)算,采用周期邊界條件和最小映像原理,而截?cái)嚅L(zhǎng)度的選取則是一項(xiàng)很有技巧性的工作,對(duì)于不同系統(tǒng),截?cái)嚅L(zhǎng)度的選取方法與選取的截?cái)嚅L(zhǎng)度是不同的,這里我們所選模擬的是一個(gè)典型的Lennard-Jones流體,其最合適的截?cái)嚅L(zhǎng)度為2.5σ[7,8].計(jì)算過(guò)程以及計(jì)算結(jié)果的給出,均運(yùn)用對(duì)比單位,對(duì)于Lennard-Jones流體,選擇m,ε,σ作為質(zhì)量、能量、長(zhǎng)度的基本單位,溫度在對(duì)比單位下就是T′=由于σ為長(zhǎng)度的基本單位,這時(shí)截?cái)嚅L(zhǎng)度即2.5.

各物理量的統(tǒng)計(jì)平均值由(3)給出.特別地,能量關(guān)于pi,qi的表達(dá)式由(1)給出.我們計(jì)算了Q從0.5到100不同取值下系統(tǒng)的演化情況,這里,我們只列出Q=0.5,1,10,100這四種典型的情況來(lái)加以研究.

圖1(P221)描述了在四種不同的情況下,系統(tǒng)能量的演化過(guò)程,從中可以看出,Q值越大,達(dá)到平衡所經(jīng)過(guò)的時(shí)間越短,并且達(dá)到平衡后,系統(tǒng)能量的演化更為平穩(wěn).另外,從中可以看出,Q取值越大,系統(tǒng)能量的波動(dòng)越小,這一點(diǎn)在圖2(P221)中可以更清楚地看出.

圖2所描述的是系統(tǒng)能量的相對(duì)波動(dòng),即系統(tǒng)能量波動(dòng)與系統(tǒng)能量平均值的比值:

圖1 能量演化過(guò)程Fig.1 Evolvement of energy

(8)式也可以稱作能量的相對(duì)誤差,但這種誤差不是由于算法的原因而產(chǎn)生的,而是由于系統(tǒng)本身與外界相互作用而形成的一種物理現(xiàn)象.從圖中可以看出,Q的值越大,相對(duì)波動(dòng)越小,系統(tǒng)的能量更不容易偏離平均值.也就是說(shuō),Q值較大的系統(tǒng)具有更強(qiáng)的穩(wěn)定性,不會(huì)由于與外界的能量交換而產(chǎn)生較大波動(dòng).

圖2 能量相對(duì)波動(dòng)Fig.2 Relative fluctuation of energy

圖3 (P222)中,六角星圖線是系統(tǒng)達(dá)到平衡時(shí)的動(dòng)量分布,實(shí)線表示經(jīng)典的波爾茲曼分布曲線,可以看出,模擬所得的結(jié)果與理論基本符合,也從另一個(gè)側(cè)面說(shuō)明模擬的有效性與正確性.

3 結(jié)論

本文討論了在不同的Q值情況下,系統(tǒng)的演化過(guò)程和相對(duì)能量波動(dòng),并考察了不同Q值的平衡狀態(tài)下動(dòng)量分布情況.從中我們可以看出:廣義質(zhì)量Q較大時(shí),(1)系統(tǒng)總能量演化更為平穩(wěn),受外界影響較小;(2)系統(tǒng)達(dá)到平衡狀態(tài)后,系統(tǒng)能量相對(duì)波動(dòng)的幅度更小.這就說(shuō)明了,Q的取值越大,對(duì)應(yīng)的系統(tǒng)漲落越小,相反,Q取值越小,系統(tǒng)漲落越大.

圖3 動(dòng)量分布Fig.3 The distribution of momentum

綜上所述,我們可以得出:同Nose-Hoover恒溫分子動(dòng)力學(xué)一樣,在Nose-Poincare恒溫分子動(dòng)力學(xué)中, Q對(duì)應(yīng)于系統(tǒng)的慣量,Q取值越大,系統(tǒng)慣量越大,外部作用所產(chǎn)生的波動(dòng)就越小.反之,Q越小,系統(tǒng)慣量就越小,外部作用所產(chǎn)生的波動(dòng)就越大,系統(tǒng)對(duì)外界的作用越敏感.

附錄:四種常見的分子動(dòng)力學(xué)算法

在分子動(dòng)力學(xué)模擬中,算法的好壞直接決定了模擬過(guò)程的計(jì)算效率與模擬結(jié)果的有效性.這里我們介紹四種常用的算法:Verlet算法,速度Verlet算法,Leap-frog算法,預(yù)校正算法.

Verlet算法:

速度Verlet算法:

Leap-frog算法:

預(yù)校正算法:

從(9)-(12)式可以看出,Verlet算能夠保持系統(tǒng)的時(shí)間反演不變性,leap-frog算法是一種辛格式,但是位置與速度不同步,而速度Verlet算法與預(yù)校正算法位置與速度同步,但是卻破壞了時(shí)間反演不變性.本文中,我們采用一種修正的leap-frog算法構(gòu)造了一種辛格式[3],不僅能夠保持系統(tǒng)長(zhǎng)時(shí)間演化的穩(wěn)定性,同時(shí)也實(shí)現(xiàn)了位置與速度的同步計(jì)算,是一種有效的分子動(dòng)力學(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

國(guó)家自然科學(xué)基金(10972125)

賈志剛(1983-),男,山西太原人,碩士生,研究方向:分子動(dòng)力學(xué)模擬.E-mail:firefox_me@yahoo.cn