MD模拟统计力学

1、MD模拟统计力学 粒子数N，温度T，体积V都相同的热力学体系组成的系综称为正则系综，正则系综的热力学体系必须处于刚性容器之中，没有任何体积变化，与环境之间也没有物质交换，但是，如果正则系综热力学体系与外界没有能量交换，则热力学体系的温度将其组成粒子的动能与势能之间的相互转换化而发生波动。为了保证正则系综热力学体系的温度恒定，每个学体系必须与一个热容巨大、温度为T的恒温热浴接触，同时，为了保证热力学体系与热浴随时处于热平衡状态，它们之间的热传导速度必须达到无穷大。 因此，正则系综热力学体系的总能量是变化的、不是固定的。1. 非Hamilton体系统计理论1.1 Liouvile方程对于任何经典力

2、学体系,给定体系的Hamilton函数,可以得到体系的Hamilton运动方程 Hamilton运动方程具有重要性质,1)Hamilton运动方程对时间反演可逆,当对运动方程的时间变量作t至-t变换时,运动方程不变,由于运动方程对时间反演可逆,对应的微观过程也对时间反演可逆,与时间方向无关;2)在体系随时间的演化过程中,体系的Hamilton函数守恒。 由于体系的Hamilton函数对应体系的总能量，它的守恒与能量守恒等价。 引入新的符号x(q,p)=,用于统一表达并处理体系的广义坐标和广义动量。根据统计系统的概念，x表示2f维相空间中的一个矢量，对应相空间中的一个点，即代表点。同时，组成统计

3、系综的任何一个经典力学体系，都有与空间中的一个代表点对应，而空间中的全部点的集合代表了统计系综的所有体系，在统计系综理论中，一个系综完全由系综分布函数确定，系综分布函数满足Liouville方程，式中，表示2f维相空间中的梯度。Liouville方程是系综分布函数守恒的直接结果，表明任意相空间体积中相点的变化等于流经该相体积边界的相点数，系综分布函数守恒也表明相空间度量守恒，即体积元是不变的，根据系综分布函数，可以计算任意力学量的系综平均，2. 非Hamilton体系统计力学 假设，某动力学体系的广义坐标和广义动量的演化不符合Hamilton运动方程，但遵循下列运动方程， 式中为体系的广义力，

4、显含时间由于体系的演化不遵循Hamilton运动方程，该动力学体系是非Hamilton体系。定义相空间的压缩率 根据统计力学理论，Hamilton体系相空间不可压缩，压缩率，相空间体积元为不变量，相反，非Hamilton体系相空间可压缩，压缩率，相空间体积元不再是不变量。 对于该非Hamilton体系，如果0时刻体系处于初始相点，t时刻体系演化到相点,则演化前后的两个相点可以通过Jacobi变换矩阵联系起来。式中，=1，随时间的演化由下列方程给出、由上式可知，只有压缩率恒为零的Hamilton体系，Jacobi矩阵才恒等于1。相反，非Hamilton体系的相空间度量或体积元按下式变换，仅当时，

5、当，.在Hamilton体系中，体积元是不变量，但在非Hamilton体系统计理论中，不变量取如下形式： 与Hamilton体系的Liouville方程对应，非Hamilton体系概率分布函数满足广义Liouville方程，在没有外界驱动力或同时显示相关的作用力的条件下，非Hamilton体系微正则系统可以通过不变量定义，如果动力系统存在M个守恒量满足则微正则系综的分布函数为：对应分配函数为3. 扩展Hamilton体系的MD模拟3.1 Nose算法 受Andersen在恒压MD模拟中通过引入广义变量扩展Hamilton函数启发，1984年Nose提出了在恒温MD模拟中通过引入额外变量扩展Ha

6、milton函数的方法，实现模拟体系与热浴之间的耦合。具体方法为引入额外的广义坐标及其对应的动量作为体系的额个自由度，利用与广义坐标对应的广义力修正体系中各粒子的速度，实现体系与热浴之间的耦合，Nose扩展体系的Hamilton函数为：扩展体系的运动方程为： Nose方法的最大贡献是通过扩展体系Hamilton函数的方法，在MD模拟中实现正则分布，成为MD模拟理论的基础。但是，Nose方法是通过对虚拟时间的等距采样来实现正则分布，但在真实时间上不能等距采样，给后期计算和处理带来困难。同时，Nose的扩展Hamilton函数不满足辛几何结构，无法采用当前在效率和稳定性上最好的辛算法，对简单体系的

7、模拟也不满足准各态历经假设。3.2 Nose-Hoover算法 为了克服Nose方法的缺陷，Hoover发展了Nose的扩展体系MD模拟方法，实现了正则系综的MD模拟，Hoover的扩展体系运动方程具有如下形式：可以证明，Nose-Hoover扩展体系中下列函数守恒，根据相空间压缩率的定义式代入Nose-Hoover方程得到得到Jacobi矩阵，相空间度量为：体系分配函数，利用函数的性质，对广义坐标积分时，仅当积分才不为零，得到现正则分布一致。3.2 Nose-Hoover算法 它为正则系综MD模拟Nose-Hoover算法的发展，通过使体系与M个广义坐标，广义动量为，广义质量为的热浴耦合的方

8、法调控温度，实现正则系综MD模拟，相应地，扩展体系运动方程具有如下形式：可以证明下列量守恒：相这空间压缩率对应的相空间度量为：3.3 对元胞体积的各向同性调整实现NPT系综 NPT系综是比正则系综更难实现的系综，在模拟过程中不但要调控温度，还必须通过调整体系的体积实现对压力的调控。因此，实现NPT系综MD模拟的关键是把元胞体积作为动力学变量，实现对压力的调控。在下列运动方程中，通过对元胞体积各同性调整实现NPT系综。式中,为与元胞体积的对数关联的广义动量;W为恒压器的广义质量,Q分别为与热浴对应的广义坐标、广义动量、广义质量，为施加的外压，为体系的内压，按照下面公式计算可以证明下列量守恒：得到

9、空间压缩率，对应的Jacobi矩阵为：相空间度量为：3.4 演化算符与差分格式 利用差分法求解经典力学体系运动方程式，随着差分过程的不断推进，差分轨迹并不收敛于实现轨迹，而是离开实际轨迹越来越远。虽然差分轨迹的误差随着时间步长的缩短而降低，但缩短时间步长需要更多的差分步为代价，才能实现相同的实际演化时间。因此，在MD模拟中需要在可以容忍误差的前提下，尽可能地延长时间步长，以减少需要进行的差分步数。在MD模拟发展的早期，普遍采用Taylor展开法设计差分格式，把坐标和速度在处展开成时间步长的幂级数，导出差分格式 。但是，采用这种方法设计的差分格式一般只能精确到时间步长的两阶，而更高阶的差分格式不

10、可避免地要求计算受力的空间导数，消耗大量计算时间。本文以Liouville算符表述的经典统计学，并在此基础上引入演化算符，用以系统地设计MD模拟的差分格式。1）Liouville算符与演化算符 在经典统计力学中，Liouville算符被子定义为：Liouville方程形式可以表达为：体系的广义坐标和广义动量随着时间的演化服从如果已知体系的初始条件，则Liouville方程的形式解为由于算符称为经典传播子或经典演化算符，简称传播子或演化算符，相应地，体系的广义坐标和广义动量的演化服从如果体系不是从0时刻开始演化，而是从时刻演化到时刻，则其演化算符写成2） Trotter定理 由于演化算符决定了体

11、系状态随时间演化的规律，因此，任何经典力学问题都归结为从Liouville算符求演化算符，为了便于计算，把Liouville算符写成两项之和但由于Liouville算符不具有对易性演化算符不能因子化 因此，无法直接利用演化算符推导MD模拟差分格式。根据Trotter定理，可用于设计差分格式，对有限的M值，得到令为单步演化的时间步长，有由此得到MD模拟中单个时间步的近似演化算符近似演化算符是酉算符，满足时间反演对称性条件，保证微观动力学过程的时间可逆性，同时，算符具有两阶精度，精确到。3） Hamilton体系的差分格式 利用近似演化算符式可以系统地设计MD模拟差分格式，具有重要意义：这时可以得

12、到利用恒等式以及近似演化算符的性质可以得到上述两式与速度Verlet差分格式一致，证明了这是满足辛对称性和时间反演对称差分格式。4） 多重时间步长差分格式 将分子体系总势能写成快速变化的分子内势能和慢速变化的分子势能之和。同时，相互作用力也可以写成分子间和分子内相互作用之和。一般地，分子内相互作用和分子间相互作用分别对应高频运动和低频运动，这时体系算符可以写成：通过把总体Liouville算符写成没有分子间作用参考态和分子间相互作用校正项之和，利用Trotter定理，演化算符可以成 这样得到了具有两种时间步长的传播子，短的时间步长对应快速变化的作用力，长的时间步长对应慢速变化的作用力。在MD模

13、拟过程中，每更新快速作用力N次，才更新慢速变化作用力一次，使计算精度在两种具有不同变化速率的作用力之间达到平衡。3.5 非Hamilton体系差分格式1. 正则系综的差分格式 利用Nose-Hoover链算法的扩展Hamilton函数，可以得到体系的Liouville算符，其中热浴对应的广义力为：把Liouville算符写成三项和的形式，再利用Trotter定理由此可推导出MD模拟的差分算法。2. NPT系综的差分格式（各向同向性） 利用NPT系综的扩展Hamilton函数,可以得到体系的Liouville算符。式中与正则系综的差分格式定义相同，但体系的运动方程可以通过下列演化算符计算，由此可

14、推导出MD模拟的差分算法。 量子力学：它是研究微观粒子的运动规律的物理学分支学科，它主要研究原子、分子、凝聚态物质，以及原子核和基本粒子的结构、性质的基础理论，它与相对论一起构成了现代物理学的理论基础。量子力学不仅是现代物理学的基础理论之一，而且在化学等有关学科和许多近代技术中也得到了广泛的应用。量子力学是非常小的领域亚原子粒子中的主要物理学理论 。该理论形成于20世纪早期，彻底改变了科学家对物质组成成分的观点。在量子世界，粒子并非是台球，而是嗡嗡跳跃的概率云，它们并不只存在一个位置，也不会从点A通过一条单一路径到达点B 。根据量子理论，粒子的行为常常像波，用于描述粒子行为

15、的“波函数”预测一个粒子可能的特性，诸如它的位置和速度，而非实际的特性。物理学中有些怪异的想法，诸如纠缠和不确定性原理，就源于量子力学 。 波函数：是量子力学中用来描述粒子的德布罗意波的函数。量子力学中描写微观系统状态的函数。在经典力学中，用质点的位置和动量（或速度）来描写宏观质点的状态，这是质点状态的经典描述方式，它突出了质点的粒子性。由于微观粒子具有波粒二象性，粒子的位置和动量不能同时有确定值（见测不准关系），因而质点状态的经典描述方式不适用于对微观粒子状态的描述。为了定量地描述微观粒子的状态，量子力学中引入了波函数，并用表示。一般来讲，波函数是空间和时间的函数，并且是复函数，即

16、=(x，y，z，t)。将爱因斯坦的“鬼场”和光子存在的概率之间的关系加以推广，玻恩假定就是粒子的概率密度，即在时刻t，在点(x,y,z)附近单位体积内发现粒子的概率。波函数因此就称为概率幅。不确定性原理：是量子力学的一个基本原理，由德国物理学家海森堡（Werner Heisenberg）于1927年提出。本身为傅立叶变换导出的基本关系：若复函数f(x)与F(k)构成傅立叶变换对，且已由其幅度的平方归一化（即f*(x)f(x)相当于x的概率密度；F*(k)F(k)/2相当于k的概率密度，*表示复共轭），则无论f(x)的形式如何，x与k标准差的乘积xk不会小于某个常数（该常数的具体形式与f(x)的

17、形式有关）。 薛定谔方程又称薛定谔波动方程，是由奥地利物理学家薛定谔提出的量子力学中的一个基本方程，也是量子力学的一个基本假定，其正确性只能靠实验来检验。它是将物质波的概念和波动方程相结合建立的二阶偏微分方程，可描述微观粒子的运动，每个微观系统都有一个相应的薛定谔方程式，通过解方程可得到波函数的具体形式以及对应的能量，从而了解微观系统的性质。）在量子力学中，体系的状态不能用力学量（例如x）的值来确定，而是要用力学量的函数（x,t），即波函数（又称概率幅，态函数）来确定，因此波函数成为量子力学研究的主要对象。力学量取值的概率分布如何，这个分布随时间如何变化，这些问题都可以通过求解波函数的薛定谔方程得到解答。这个方程是奥地利物理学家薛定谔于1926年提出的，它是量子力学最基本的方程之一，在量子力学中的地位与牛顿方程在经典力学中的地位相当。第一性原理：根据原子核和电子互相作用的原理及其基本运动规律，运用量子力学原理，从具体要求出发，经过一些近似处理后直接求解薛定谔方程的算法，习惯上称为第一性原理。第一性原理通常是跟计算联系在一起的，是指在进行计算的时候除了告诉程序你所使用的原子和他们的位置外，没有其它的实验的，经验的或者半经验的参量，且具有很好的移植性。作为评价事物的依据，第一性原理和经验参数是两个极端。第一性原理是某些硬性规定