MD模拟统计力学

1MD模拟统计力学粒子数N，温度T，体积V都相同的热力学体系组成的系综称为正则系综，正则系综的热力学体系必须处于刚性容器之中，没有任何体积变化，与环境之间也没有物质交换，但是，如果正则系综热力学体系与外界没有能量交换，则热力学体系的温度将其组成粒子的动能与势能之间的相互转换化而发生波动。为了保证正则系综热力学体系的温度恒定，每个学体系必须与一个热容巨大、温度为T的恒温热浴接触，同时，为了保证热力学体系与热浴随时处于热平衡状态，它们之间的热传导速度必须达到无穷大。因此，正则系综热力学体系的总能量是变化的、不是固定的。1.非Hamilton体系统计理论1.1Liouvile方程对于任何经典力学体系,给定体系的Hamilton函数,),...,,(2),...,,,;,...,,(),(2112321321ffiiiffqqqumpppppqqqqHqpH可以得到体系的Hamilton运动方程iifiiiiiifqqqquqHpmppHq),...,(21...Hamilton运动方程具有重要性质,1)Hamilton运动方程对时间反演可逆,当对运动方程的时间变量作t至-t变换时,运动方程不变,由于运2动方程对时间反演可逆,对应的微观过程也对时间反演可逆,与时间方向无关;2)在体系随时间的演化过程中,体系的Hamilton函数守恒。0)()(1..1.iifiiiiifiiiqHpHpHqHppHqqHdtdH由于体系的Hamilton函数对应体系的总能量，它的守恒与能量守恒等价。引入新的符号x(q,p)=),...,,,...,(2121ffpppqqq,用于统一表达并处理体系的广义坐标和广义动量。根据统计系统的概念，x表示2f维相空间中的一个矢量，对应相空间中的一个点，即代表点。同时，组成统计系综的任何一个经典力学体系，都有与空间中的一个代表点对应，而空间中的全部点的集合代表了统计系综的所有体系，在统计系综理论中，一个系综完全由系综分布函数),(),,(txftpqf确定，系综分布函数满足Liouville方程，0),(),(),(.txfxttxfdttxdf式中，表示2f维相空间中的梯度。Liouville方程是系综分布函数),(txf守恒的直接结果，表明任意相空间体积中相点的变化等于流经该相体积边界的相点数，系综分布函数),(txf守恒也表明相空间度量守恒，即体积元dqdpdpdqdxdxdfff2是不变的，根据系综分布函数，可以计算任意力学量)(xA的系综平均，),()(),(txdxfxAtxdxfA2.非Hamilton体系统计力学3假设，某动力学体系的广义坐标和广义动量的演化不符合Hamilton运动方程，但遵循下列运动方程，),(.txx式中),(tx为体系的广义力，显含时间0t由于体系的演化不遵循Hamilton运动方程，该动力学体系是非Hamilton体系。定义相空间的压缩率.),(xtxk根据统计力学理论，Hamilton体系相空间不可压缩，压缩率0),(txk，相空间体积元fdxd2为不变量，相反，非Hamilton体系相空间可压缩，压缩率0),(txk，相空间体积元fdxd2不再是不变量。对于该非Hamilton体系，如果0时刻体系处于初始相点0x，t时刻体系演化到相点tx,则演化前后的两个相点可以通过Jacobi变换矩阵联系起来。),...(),...();(2010210fftttxxxxxxJ式中，);(0txxJ=1，);(0txxJ随时间的演化由下列方程给出、);(),(),(00xxJtxkdtxxdJttt由上式可知，只有压缩率),(txk恒为零的Hamilton体系，Jacobi矩阵);(0txxJ才恒等于1。相反，非Hamilton体系的相空间度量或体积元按下式变换，001);(dxxxJdxt仅当1);(0txxJ时，0dxdxt，当1);(0txxJ，0dxdxt.在Hamilton体系4中，体积元d是不变量，但在非Hamilton体系统计理论中，不变量取如下形式：dxtxgd),('),(),(txwetxg),(),(txkdttxdw与Hamilton体系的Liouville方程对应，非Hamilton体系概率分布函数),(txf满足广义Liouville方程，0)()(.fgxtfg在没有外界驱动力或同时显示相关的作用力的条件下，非Hamilton体系微正则系统可以通过不变量定义，如果动力系统存在M个守恒量)(xK满足MdtxdK,...,1,0)(则微正则系综的分布函数为：MKxKxf1_))(()(对应分配函数为MMKxKxgdxKKVN1___1))(()(),...,,(3.扩展Hamilton体系的MD模拟3.1Nose算法受Andersen在恒压MD模拟中通过引入广义变量扩展Hamilton函数启发，1984年Nose提出了在恒温MD模拟中通过引入额外变量扩展Hamilton函数的方法，实现模拟体系与热浴之间的耦合。具体5方法为引入额外的广义坐标s及其对应的动量sp作为体系的额个自由度，利用与广义坐标s对应的广义力修正体系中各粒子的速度，实现体系与热浴之间的耦合，Nose扩展体系的Hamilton函数为：sTNkQprrrusmppsprHBsNNiiisln32),...(2),,,(22,1122扩展体系的运动方程为：NiBiisssiNiiiiiiTNksmpssHpQppHsrrrruqHpsmppHr122..21.2.)3(1),...,(Nose方法的最大贡献是通过扩展体系Hamilton函数的方法，在MD模拟中实现正则分布，成为MD模拟理论的基础。但是，Nose方法是通过对虚拟时间的等距采样来实现正则分布，但在真实时间上不能等距采样，给后期计算和处理带来困难。同时，Nose的扩展Hamilton函数不满足辛几何结构，无法采用当前在效率和稳定性上最好的辛算法，对简单体系的模拟也不满足准各态历经假设。3.2Nose-Hoover算法为了克服Nose方法的缺陷，Hoover发展了Nose的扩展体系MD模拟方法，实现了正则系综的MD模拟，Hoover的扩展体系运动方程具有如下形式：6NiBiiiiiiiiTNkmppQpQppfpmpr12....3可以证明，Nose-Hoover扩展体系中下列函数守恒，CTNkQpprHTNkQprrumppprHBBNNiii32),(32),...,(2),,,('22112根据相空间压缩率的定义式....1.)(),(pppprrxtxkiiiNii代入Nose-Hoover方程得到.3),(Ntxk得到Jacobi矩阵，NetxJ3),(相空间度量为：)3exp(Ng体系分配函数，)32),(()3exp(),,(2CTNkQpprHNddpdrdpEVNB利用函数的性质，对广义坐标积分时，仅当)2),((312CQpprHTNkB积分才不为零，得到7),,())2),((1exp(31),,(2TVNQQpprHCTkdpdrdpTNkEVNBB现正则分布一致。3.2Nose-Hoover算法它为正则系综MD模拟Nose-Hoover算法的发展，通过使体系与M个广义坐标j，广义动量为jp，广义质量为iQ的热浴耦合的方法调控温度，实现正则系综MD模拟，相应地，扩展体系运动方程具有如下形式：TkQppMpQpTkQpppQpTNkmppMjQppQpfpmprBMMMjjiBjijMiBiijijiiiiii112.1111.12212.1.11...1,...2,1,)()3(,...2,1,可以证明下列量守恒：MjjBBNMjjjNiiiTkTNkrruQpmppprH21112123),...,(22),,,('相这空间压缩率MjjNtxk2.13),(对应的相空间度量为：)3exp(21MjjNg83.3对元胞体积的各向同性调整实现NPT系综NPT系综是比正则系综更难实现的系综，在模拟过程中不但要调控温度，还必须通过调整体系的体积实现对压力的调控。因此，实现NPT系综MD模拟的关键是把元胞体积作为动力学变量，实现对压力的调控。在下列运动方程中，通过对元胞体积各同性调整实现NPT系综。NiBiiNiiiextiniiiiiiiiTkNWpmppQppQpmpNPPVpWVpVpQppWpNfprWpmpr122..12....)13(1)(33)11(式中,p为与元胞体积的对数)ln(310VV关联的广义动量;W为恒压器的广义质量,,p,Q分别为与热浴对应的广义坐标、广义动量、广义质量，extP为施加的外压，inP为体系的内压，按照下面公式计算)),(32(31112VVruVfrmpVPiNiiNiiiin可以证明下列量守恒：VPTkNVruQpWpmpppVprHextBNiii)13(),(222),,,,,(2122'得到空间压缩率，9..)13(),(Nxtxk对应的Jacobi矩阵为：))13exp((),(NtxJ相空间度量为：))13exp((Ng3.4演化算符与差分格式利用差分法求解经典力学体系运动方程式，随着差分过程的不断推进，差分轨迹并不收敛于实现轨迹，而是离开实际轨迹越来越远。虽然差分轨迹的误差随着时间步长t的缩短而降低，但缩短时间步长t需要更多的差分步为代价，才能实现相同的实际演化时间。因此，在MD模拟中需要在可以容忍误差的前提下，尽可能地延长时间步长t，以减少需要进行的差分步数。在MD模拟发展的早期，普遍采用Taylor展开法设计差分格式，把坐标)(tr和速度)(tV在0tt处展开成时间步长t的幂级数，导出差分格式。但是，采用这种方法设计的差分格式一般只能精确到时间步长t的两阶，而更高阶的差分格式不可避免地要求计算受力的空间导数，消耗大量计算时间。本文以Liouville算符表述的经典统计学，并在此基础上引入演化算符，用以系统地设计MD模拟的差分格式。1）Liouville算符与演化算符在经典统计力学中，Liouville算符被子定义为：10NiiiiipqHqpHiL1)(Liouville方程形式可以表达为：Lftfi体系的广义坐标和广义动量随着时间的演化服从iLxx.如果已知体系的初始条件)0(x，则Liouville方程的形式解为)0()(xetxiLt由于算符iLtetU)()(tU称为经典传播子或经典演化算符，简称传播子或演化算符，相应地，体系的广义坐标和广义动量的演化服从)0()()(xtUtx如果体系不是从0时刻开始演化，而是从1t时刻演化到2t时刻，则其演化算符写成)(121212)(),(ttiLettUttU2）Trotter定理由于演化算符)(tU决定了体系状态随时间演化的规律，因此，任何经典力学问题都归结为从Liouville算符求演化算符)(tU，为了便于计算，把Liouville算符写成两项之和21LLL