MD模拟统计力学

MD模拟统计力学粒子数N,温度T,体积V都相同的热力学体系组成的系综称为正则系综，正则系综的热力学体系必须处于刚性容器之中，没有任何体积变化，与环境之间也没有物质交换，但是，如果正则系综热力学体系与外界没有能量交换，则热力学体系的温度将其组成粒子的动能与势能之间的相互转换化而发生波动。为了保证正则系综热力学体系的温度恒定，每个学体系必须与一个热容巨大、温度为T的恒温热浴接触，同时，为了保证热力学体系与热浴随时处于热平衡状态，它们之间的热传导速度必须达到无穷大。因此，正则系综热力学体系的总能量是变化的、不是固定的。1.非Hamilton体系统计理论Liouvile方程对于任何经典力学体系，给定体系的Hamilton函数，H(p,q)二H(q,q,q，…q;p,p,p,...,p)二丈 +u(q,q,...,q)TOC\o"1-5"\h\z123f123f 2m 12fi=1 i可以得到体系的Hamilton运动方程dHpq= =—idp mV i idH du(q,q，…q) 「p=-——= 1_2 j=fidq dq ii iHamilton运动方程具有重要性质,1)Hamilton运动方程对时间反演可逆,当对运动方程的时间变量作t至-t变换时，运动方程不变,由于运动方程对时间反演可逆,对应的微观过程也对时间反演可逆，与时间方向无关;2)在体系随时间的演化过程中，体系的Hamilton函数守恒。dH二工(岂q-+岂p)二工(迺匹-廻翌)二0dt dqidpi dqdpdpdqi=1 i i i=1 ii ii由于体系的Hamilton函数对应体系的总能量，它的守恒与能量守恒等价。引入新的符号x三(q,p)=(q,q，…q,p,p，…p)，用于统一表达并处理1 2f1 2f体系的广义坐标和广义动量。根据统计系统的概念，x表示2f维相空间中的一个矢量，对应相空间中的一个点，即代表点。同时，组成统计系综的任何一个经典力学体系，都有与空间中的一个代表点对应，而空间中的全部点的集合代表了统计系综的所有体系，在统计系综理论中，一个系综完全由系综分布函数f(q，p，t)三f(x，t)确定，系综分布函数满足Liouville方程，df(x，t)=df(x，t)+xVf(x,t)=0dt dt式中，v表示2f维相空间中的梯度。Liouville方程是系综分布函数f(x,t)守恒的直接结果，表明任意相空间体积中相点的变化等于流经该相体积边界的相点数，系综分布函数f(x,t)守恒也表明相空间度量守恒，即体积元dr=dx2f=dx=dqfdpf=dqdp是不变的，根据系综分布函数，可以计算任意力学量A(x)的系综平均,Jdxf(x,t)A(x)〈A= J—Jdxf(x,t)2.非Hamilton体系统计力学假设，某动力学体系的广义坐标和广义动量的演化不符合Hamilton运动方程，但遵循下列运动方程，X=£(X,t)式中8(X,t)为体系的广义力，显含时间t由于体系的演化不遵循0Hamilton运动方程，该动力学体系是非Hamilton体系。定义相空间的压缩率k(x,t)=V・X根据统计力学理论，Hamilton体系相空间不可压缩，压缩率k(x,t)三0，相空间体积元亦=如为不变量，相反，非Hamilton体系相空间可压缩，压缩率k(X,t)丰0，相空间体积元dr=dx2f不再是不变量。对于该非Hamilton体系，如果0时刻体系处于初始相点x，t时0刻体系演化到相点X，则演化前后的两个相点可以通过Jacobi变换矩t阵联系起来。)Q(Xl,...X2f)J(X；X)二tt

t0 3(X1,...X2f)00式中，J(X;X)=1，J(X;X)随时间的演化由下列方程给出、0t 0tdJ(x,x) )t0=k(x,t)J(x;x)

dt tt0由上式可知，只有压缩率k(x,t)恒为零的Hamilton体系，Jacobi矩阵J(x；x)才恒等于1。相反，非Hamilton体系的相空间度量或体积元0t按下式变换，dX=J(x;x)dX1 t0 0仅当J(X;X)三1时，dx三dx，当J(X;X)丰1，dx丰dx•在Hamilton体系0t t 0 0t t 0中，体积元dr是不变量，但在非Hamilton体系统计理论中，不变量取如下形式：dP=七g(x,t)dx

vg(x,t)=e-w(x,t)如卫=k(x,t)dt与Hamilton体系的Liouville方程对应，非Hamilton体系概率分布函数f(x,t)满足广义Liouville方程，+V(x話f)=0ct在没有外界驱动力或同时显示相关的作用力的条件下，非Hamilton体系微正则系统可以通过不变量定义，如果动力系统存在M个守恒量K九(x)满足则微正则系综的分布函数为:Mo

f(x)= 5(K(x)-KJ九 九九=1对应分配函数为Q(N,V,K，…K)=Jdx/g(x)M5(K(x)—KJ1 M w 九 人九=13.扩展Hamilton体系的MD模拟Nose算法受Andersen在恒压MD模拟中通过引入广义变量扩展Hamilton函数启发，1984年Nose提出了在恒温MD模拟中通过引入额外变量扩展Hamilton函数的方法，实现模拟体系与热浴之间的耦合。具体

方法为引入额外的广义坐标s及其对应的动量p作为体系的额个自s由度，利用与广义坐标s对应的广义力修正体系中各粒子的速度，实现体系与热浴之间的耦合，Nose扩展体系的Hamilton函数为：H(r,p,s,p)三£ Pi—二u(rr，…r)+ +3NkTlnss 2ms2 1,2N 2Q Bi=1 i扩展体系的运动方程为:dHdpdHdu(r,rdHdu(r,r,.・.r12drdHdHds-3NkT)i=1Nose方法的最大贡献是通过扩展体系Hamilton函数的方法，在MD模拟中实现正则分布，成为MD模拟理论的基础。但是，Nose方法是通过对虚拟时间的等距采样来实现正则分布，但在真实时间上不能等距采样，给后期计算和处理带来困难。同时，Nose的扩展Hamilton函数不满足辛几何结构，无法采用当前在效率和稳定性上最好的辛算法，对简单体系的模拟也不满足准各态历经假设。Nose-Hoover算法为了克服Nose方法的缺陷，Hoover发展了Nose的扩展体系MD模拟方法，实现了正则系综的MD模拟，Hoover的扩展体系运动方程具有如下形式：

p.=f-p.咚iiiQp

q=—qQp仝p-3NkTqmBi=1i可以证明，Nose-Hoover扩展体系中下列函数守恒,H'(r,p,q,p)三兰 +u(rr)+q+3NkTq二H(r,p)+丐+3NkTq二Cq2m1N2Q B 2QBi=1 i根据相空间压缩率的定义式k(xt)=vx氾(dTr+appi)+紬i=1i i代入Nose-Hoover方程得到k(x,t)=-3N得到Jacobi矩阵,J(x,t)二e-3Nn相空间度量为:爲=exp(叫P2P2仝+3Nv;-C)Q(N,V,E)=JdpJdrJdpdqexp(3N)5(H(r,p)+TOC\o"1-5"\h\zq q 2Q利用5函数的性质，对广义坐标耳积分时，仅当i p2q二 (H(r,p)+q-C)3NkT 2QB积分才不为零，得到

13NkTB(C-H(r,p)13NkTB(C-H(r,p)-gQ(N,V,T)B现正则分布一致。Nose-Hoover算法它为正则系综MD模拟Nose-Hoover算法的发展，通过使体系与M个广义坐标耳，广义动量为p，广义质量为q的热浴耦合的方法j 中 '调控温度，实现正则系综MD模拟，相应地，扩展体系运动方程具有如下形式:pr=iimip=f—严piiQi

匕1n= ,j=i,2,...mjQp=(兰佯-3NkT)-厶2p1m BQn1TOC\o"1-5"\h\zi=1i 2p pp.=(心1一kT)一ni+1p,=1,2,...M-1./QbQ帀j-i j+ip' p-1-kTMQ B、 M-1)三迓止+瓦隘)三迓止+瓦隘2mi j=12Qji=1+u(r,...,r)+3NkT1N Br|i相这空间压缩率k(x,t)=—3N—艺rj=2对应的相空间度量为：屈=exp(3N屈=exp(3Nij=23.3对元胞体积的各向同性调整实现NPT系综NPT系综是比正则系综更难实现的系综，在模拟过程中不但要调控温度，还必须通过调整体系的体积实现对压力的调控。因此，实现NPT系综MD模拟的关键是把元胞体积作为动力学变量，实现对压力的调控。在下列运动方程中，通过对元胞体积各同性调整实现NPT系综。r=匕+PrimWiip.=f-(1+N)Wp—牛p.NW1QTOC\o"1-5"\h\z3Vp= £Wp=3V(P-P)+丄迓聲-厶p£ inexti=1 ipQp +伫-(3N+1)kT耳 mW b' i=1 i式中,p为与元胞体积的对数£=!ln(V)关联的广义动量；W为恒压器£ 3V0的广义质量凋,p,Q分别为与热浴对应的广义坐标、广义动量、广义质量，p为施加的外压，p为体系的内压，按照下面公式计算ext inPn三sV忆2m+7f-3V將)1ii=1H(H(r,p,V,p弭,p)仝£ 耳i=1釜+养+S+U(r，V)+(3N+叽PxtVi得到空间压缩率,k(x,t)=Vx=-(3N+1)n对应的Jacobi矩阵为:J(x,t)=exp((-3N+1)n)相空间度量为:y/g=exp((3N+1)n)3.4演化算符与差分格式利用差分法求解经典力学体系运动方程式，随着差分过程的不断推进，差分轨迹并不收敛于实现轨迹，而是离开实际轨迹越来越远。虽然差分轨迹的误差随着时间步长A的缩短而降低，但缩短时间步长At t需要更多的差分步为代价，才能实现相同的实际演化时间。因此，在MD模拟中需要在可以容忍误差的前提下，尽可能地延长时间步长A，以减少需要进行的差分步数。t在MD模拟发展的早期，普遍采用Taylor展开法设计差分格式，把坐标r(t)和速度V(t)在t=t0处展开成时间步长At的幕级数，导出差分格式。但是，采用这种方法设计的差分格式一般只能精确到时间步长At的两阶,而更高阶的差分格式不可避免地要求计算受力的空间导数，消耗大量计算时间。本文以Liouville算符表述的经典统计学，并在此基础上引入演化算符，用以系统地设计MD模拟的差分格式。1)Liouville算符与演化算符在经典统计力学中，Liouville算符被子定义为：

iL仝(i=16Ha6piL仝(i=16Ha6pdqi ii iLiouville方程形式可以表达为:idfidfdt=Lf体系的广义坐标和广义动量随着时间的演化服从x=iLx如果已知体系的初始条件x(0)，则Liouville方程的形式解为x(t)=eiLtx(0)由于算符U(t)=eiLtU(t)称为经典传播子或经典演化算符，简称传播子或演化算符，相应地，体系的广义坐标和广义动量的演化服从x(t)=U(t)x(0)如果体系不是从0时刻开始演化，而是从t时刻演化到t时刻，则其12演化算符写成U(t,t)=U(t-1)=eiL(『2-[)21212)Trotter定理由于演化算符U(t)决定了体系状态随时间演化的规律，因此，任何经典力学问题都归结为从Liouville算符求演化算符u(t)，为了便于计算，把Liouville算符写成两项之和L=L+L12但由于Liouville算符不具有对易性[L,L]三LL—LL丰0121221

演化算符不能因子化eiLt=ei(L]+厶2)t丰eiL^eiLf因此，无法直接利用演化算符推导MD模拟差分格式。根据Trotter定理，ei(人+l2)t=lim(e2meme2m)mMs可用于设计差分格式，对有限的M值，得到ei(L+厶2)〃M沁e2meme2m令At二t/M为单步演化的时间步长，有ei(L]+乙2)AtUejAt/2eiL1AteiL2At/2由此得到MD模拟中单个时间步的近似演化算符U(At)=eiL2At/2eiL1AteiL2At/2近似演化算符是u(At)酉算符，满足时间反演对称性条件，保证微观动力学过程的时间可逆性，同时，u(At)算符具有两阶精度，精确到At2。3)Hamilton体系的差分格式利用近似演化算符式可以系统地设计MD模拟差分格式，具有重要意义:TOC\o"1-5"\h\z5号QHQyd

iL三 = v——1 QpQr iQri=1 ii i=1 iiL三y更丄正fa2 QrQp mQvi=1 ii i=1 ii这时可以得到U(At)=exp(乞y£—)exp(Atyv—)exp(A —)i=12mQv iQr 2mi=1利用恒等式dexp(c )g(q)=g(q+c)oq以及近似演化算符的性质r(At)uU(At)r(0)TOC\o"1-5"\h\zi iv(At)uU(At)v(0)

i i可以得到iv(At)二v(0)+¥(f(0)+f(At))i i 2i i上述两式与速度Verlet差分格式一致，证明了这是满足辛对称性和时间反演对称差分格式。4)多重时间步长差分格式将分子体系总势能写成快速变化的分子内势能u和慢速变化的intra分子势能u之和。inter同时，相互作用力也可以写成分子间和分子内相互作用之和。f=fintra+finterii i一般地，分子内相互作用和分子间相互作用分别对应高频运动和低频运动，这时体系算符可以写成：iL仝(v2+(化+a)2)iOr mmOvi=1 i i i i通过把总体Liouville算符写成没有分子间作用参考态l和分子间相ref互作用校正项AL之和，L=L +ALref

dfintradTOC\o"1-5"\h\ziL 二乙(v +dfintradrf 1dr m dvi=1 i i\iAiAL=Z(化2mdv利用Trotter定理，演化算符可以成U(At)=eiALAt/2eiLr(eftcezwceiALeiALAt/2=exp( 乙iAtyfinterd2m dv这样得到了具有两种时间步长的传播子，短的时间步长81对应快速变化的作用力，长的时间步长At=n81对应慢速变化的作用力。在MD模拟过程中，每更新快速作用力N次，才更新慢速变化作用力一次，使计算精度在两种具有不同变化速率的作用力之间达到平衡。3.5非Hamilton体系差分格式正则系综的差分格式利用Nose-Hoover链算法的扩展Hamilton函数，可以得到体系的Liouville算符，idr mdv 叶1dv n/'dqnidr mdv 叶1dv n/'dqnjn/+1dq Mdvj nMi=1 ii=1ii i=1 i j=1 j j=1其中热浴对应的广义力为:QG=£mv2一3NkTQG=Qv2一kT,j>1把Liouville算符写成三项和的形式,

LL=L+L+L1 2NHC再利用Trotter定理exp(iLAt)沁exp(iL 空)exp(iL空)exp(iL乞)exp(iL竺)exp(iL 乞)nhc2 12 22 12 nhc2由此可推导出MD模拟的差分算法。NPT系综的差分格式(各向同向性)利用NPT系综的扩展Hamilton函数，可以得到体系的Liouville算符。1、yQ QQTOC\o"1-5"\h\ziL=iL一(1+ )乙vv +(G一vv) +v一nhcN£iQv £ £nQv eQsi=1 i £歹/ \Qyf(r)+乙(v+vr) +乙——i£iQr Qvi=1 i i=1 iG=1((1+1)lLmv2+迓rf(r)-3VQu(r'V)-3PV)£VW N ii ii QV exti=1 i=11£=lg(V/V)3 0iL与正则系综的差分格式定义相同，但NHCG=—(ymv2+Wv2-(3N+1)kT)1Q ii £ B1i=1体系的运动方程可以通过下列演化算符计算,iLp-(1+丄疋vvV+(G-vviLpNsivi£gqQviL眾竺亘iL眾竺亘mdvi=1 iiiL=艺(v+vr)icii=1、QQ+viQr£Qsi由此可推导出MD模拟的差分算法。量子力学：它是研究微观粒子的运动规律的物理学分支学科，它主要研究原子、分子、凝聚态物质，以及原子核和基本粒子的结构、性质的基础理论，它与相对论一起构成了现代物理学的理论基础。量子力学不仅是现代物理学的基础理论之一，而且在化学等有关学科和许多近代技术中也得到了广泛的应用。量子力学是非常小的领域一一亚原子粒子中的主要物理学理论。该理论形成于20世纪早期，彻底改变了科学家对物质