第六章分子动力学方法 6.1引言 对于一个多粒子体系的实验观测物理量的数值可以由总的平均得到。但是由于实验体系又非常 大,我们不可能计算求得所有涉及到的态的物理量数值的总平均。按照产生位形变化的方法,我们有 两类方法对有限的一系列态的物理量做统计平均: 第一类是随机模拟方法。它是实现Gibs的统计力学途径。在此方法中,体系位形的转变是通过 马尔科夫( Markov)过程,由随机性的演化引起的。这里的马尔科夫过程相当于是内禀动力学在概率方 面的对应物。该方法可以被用到没有任何内禀动力学模型体系的模拟上。随机模拟方法计算的程序简 单,占内存少,但是该方法难于处理非平衡态的问题
第六章 分子动力学方法 6.1 引言 对于一个多粒子体系的实验观测物理量的数值可以由总的平均得到。但是由于实验体系又非常 大，我们不可能计算求得所有涉及到的态的物理量数值的总平均。按照产生位形变化的方法，我们有 两类方法对有限的一系列态的物理量做统计平均： 第一类是随机模拟方法。它是实现 Gibbs 的统计力学途径。在此方法中，体系位形的转变是通过 马尔科夫(Markov)过程，由随机性的演化引起的。这里的马尔科夫过程相当于是内禀动力学在概率方 面的对应物。该方法可以被用到没有任何内禀动力学模型体系的模拟上。随机模拟方法计算的程序简 单，占内存少，但是该方法难于处理非平衡态的问题
另一类为确定性模拟方法,即统计物理中的所谓分子动力学方法( Molecular Dynamics Method。 这种方法广泛地用于研究经典的多粒子体系的研究中。该方法是按该体系内部的内禀动力学规律来计 算并确定位形的转变。它首先需要建立一组分子的运动方程,并通过直接对系统中的一个个分子运动 方程进行数值求解,得到每个时刻各个分子的坐标与动量,即在相空间的运动轨迹,再利用统计计算 方法得到多体系统的静态和动态特性,从而得到系统的宏观性质。因此,分子动力学模拟方法可以看 作是体系在一段时间内的发展过程的模拟。在这样的处理过程中我们可以看出:分子动力学方法中不 存在任何随机因素。 系统的动力学机制决定运动方程的形式: 在分子动力学方法处理过程中,方程组的建立是通过对物理体系的微观数学描述给出的。在这个 微观的物理体系中,每个分子都各自服从经典的牛顿力学。每个分子运动的内禀动力学是用理论力学 上的哈密顿量或者拉格朗日量来描述,也可以直接用牛顿运动方程来描述。这种方法可以处理与时间 有关的过程,因而可以处理非平衡态问题。但是使用该方法的程序较复杂,计算量大,占内存也多
另一类为确定性模拟方法，即统计物理中的所谓分子动力学方法（Molecular Dynamics Method）。 这种方法广泛地用于研究经典的多粒子体系的研究中。该方法是按该体系内部的内禀动力学规律来计 算并确定位形的转变。它首先需要建立一组分子的运动方程，并通过直接对系统中的一个个分子运动 方程进行数值求解，得到每个时刻各个分子的坐标与动量，即在相空间的运动轨迹，再利用统计计算 方法得到多体系统的静态和动态特性, 从而得到系统的宏观性质。因此，分子动力学模拟方法可以看 作是体系在一段时间内的发展过程的模拟。在这样的处理过程中我们可以看出：分子动力学方法中不 存在任何随机因素。 系统的动力学机制决定运动方程的形式： 在分子动力学方法处理过程中，方程组的建立是通过对物理体系的微观数学描述给出的。在这个 微观的物理体系中，每个分子都各自服从经典的牛顿力学。每个分子运动的内禀动力学是用理论力学 上的哈密顿量或者拉格朗日量来描述，也可以直接用牛顿运动方程来描述。这种方法可以处理与时间 有关的过程，因而可以处理非平衡态问题。但是使用该方法的程序较复杂，计算量大，占内存也多
适用范围广泛: 原则上,分子动力学方法所适用的微观物理体系并无什么限制。这个方法适用的体系既可以是少 体系统,也可以是多体系统;既可以是点粒子体系,也可以是具有内部结构的体系;处理的微观客体 既可以是分子,也可以是其它的微观粒子。 自五十年代中期开始,分子动力学方法得到了广泛的应用。它与蒙特卡洛方法一起已经成为计算 机模拟的重要方法。应用分子动力学方法取得了许多重要成果,例如气体或液体的状态方程、相变问 题、吸附问题等,以及非平衡过程的研究。其应用已从化学反应、生物学的蛋白质到重离子碰撞等广 泛的学科研究领域
适用范围广泛： 原则上，分子动力学方法所适用的微观物理体系并无什么限制。这个方法适用的体系既可以是少 体系统，也可以是多体系统；既可以是点粒子体系，也可以是具有内部结构的体系；处理的微观客体 既可以是分子，也可以是其它的微观粒子。 自五十年代中期开始，分子动力学方法得到了广泛的应用。它与蒙特卡洛方法一起已经成为计算 机模拟的重要方法。应用分子动力学方法取得了许多重要成果，例如气体或液体的状态方程、相变问 题、吸附问题等，以及非平衡过程的研究。其应用已从化学反应、生物学的蛋白质到重离子碰撞等广 泛的学科研究领域
实际使用的限制: 实际上,分子动力学模拟方法和随机模拟方法一样都面临着两个基本限制:一个是有限观测时间 的限制;另一个是有限系统大小的限制。通常人们感兴趣的是体系在热力学极限下(即粒子数目趋于 无穷时)的性质。但是计算机模拟允许的体系大小要比热力学极限小得多,因此可能会出现有限尺寸 效应。为了减小有限尺寸效应,人们往往引入周期性、全反射、漫反射等边界条件。当然边界条件的 引入显然会影响体系的某些性质。 数值求解时的离散化方法: 对体系的分子运动方程组采用计算机进行数值求解时,需要将运动方程离散化为有限差分方程(参见 第四章)。常用的求解方法有欧拉法、龙格-库塔法、辛普生法等(参见附录D。数值计算的误差阶 数显然取决于所采用的数值求解方法的近似阶数。原则上,只要计算机计算速度足够大,内存足够多, 我们可以使计算误差足够小
实际使用的限制： 实际上，分子动力学模拟方法和随机模拟方法一样都面临着两个基本限制：一个是有限观测时间 的限制；另一个是有限系统大小的限制。通常人们感兴趣的是体系在热力学极限下（即粒子数目趋于 无穷时）的性质。但是计算机模拟允许的体系大小要比热力学极限小得多，因此可能会出现有限尺寸 效应。为了减小有限尺寸效应，人们往往引入周期性、全反射、漫反射等边界条件。当然边界条件的 引入显然会影响体系的某些性质。 数值求解时的离散化方法： 对体系的分子运动方程组采用计算机进行数值求解时，需要将运动方程离散化为有限差分方程(参见 第四章)。常用的求解方法有欧拉法、龙格-库塔法、辛普生法等（参见附录 D）。数值计算的误差阶 数显然取决于所采用的数值求解方法的近似阶数。原则上，只要计算机计算速度足够大，内存足够多， 我们可以使计算误差足够小
6.2分子动力学基础知识 分子运动方程及其数值求解 采用分子动力学方法时,必须对一组分子运动微分方程做数值求解。从计算数学的角度来看,这 是个求一个初值问题的微分方程的解。实际上计算数学为了求解这种问题已经发展了许多的算法,但 是并不是所有的这些算法都可以用来解决物理问题。 例子: 以一个一维谐振子为例,来看一下如何用计算机数值计算方法求解初值问题。一维谐振子的经典 哈密顿量为: H kx 这里的哈密顿量(即能量)为守恒量
6．2 分子动力学基础知识 一、分子运动方程及其数值求解 采用分子动力学方法时，必须对一组分子运动微分方程做数值求解。从计算数学的角度来看，这 是个求一个初值问题的微分方程的解。实际上计算数学为了求解这种问题已经发展了许多的算法，但 是并不是所有的这些算法都可以用来解决物理问题。 例子： 以一个一维谐振子为例，来看一下如何用计算机数值计算方法求解初值问题。一维谐振子的经典 哈密顿量为： 2 2 2 1 2 kx m p H += . (6.2.1) 这里的哈密顿量（即能量）为守恒量
假定初始条件为x0p0),则它的哈密顿方程是对时间的一阶微分方程 dx aH P 中 d t (6.2.2) 计算在相空间中的运动轨迹(p():采用有限差分法,将微分方程变为有限差分方程,以便在计算 机上做数值求解,并得到空间坐标和动量随时间的演化关系。 首先,取差分计算的时间步长为h,采用一阶微分形式的向前差商表示,它是直接运用展开到h的一 阶泰勒展开公式 (+6)=/()+h4+0矿), 即得到 d、f(+h-f() (6.2.3
假定初始条件为 () () px 0,0 ，则它的哈密顿方程是对时间的一阶微分方程 m p p H dt dx = ∂ ∂ = ， kx x H dt dp −= ∂ ∂ −= . (6.2.2) 计算在相空间中的运动轨迹( ( ), (tptx ))：采用有限差分法，将微分方程变为有限差分方程，以便在计算 机上做数值求解，并得到空间坐标和动量随时间的演化关系。 首先，取差分计算的时间步长为 ，采用一阶微分形式的向前差商表示，它是直接运用展开到 的一 阶泰勒展开公式 h h ( ) () )( 2 hO dt df htfhtf ++=+ ， 即得到 ( ) () h tfhtf dt df −+ ≈ ， (6.2.3)
则微分方程6.2.2)可以被改写为差分形式 dx x(+h)-x( p() d t h (6.2.4) 中_p+)p2=-6 (6.2.5 将上面两个公式整理后,我们得到解微分方程(6.2.2)的欧拉( Euler)算法(参见附录C) x(t+h)=x(t (6.2.6) +h hkx (6.2.7) 这是x()p)的一组递推公式。有了初始条件x()p(0),就可以一步一步地使用前一时刻的坐标、动量 值确定下一时刻的坐标、动量值。这个方法是一步法的典型例子
则微分方程(6.2.2)可以被改写为差分形式 ( ) ( ) ( ) m tp h txhtx dt dx = −+ = . (6.2.4) ( ) ( ) ( )tkx h tphtp dt dp −= −+ = . (6.2.5) 将上面两个公式整理后, 我们得到解微分方程(6.2.2)的欧拉(Euler)算法（参见附录 C）： ( ) () ( ) m thp txhtx +=+ . (6.2.6) ( ) =+ ( ) − (thkxtphtp ) . (6.2.7) 这是 ( ), (tptx )的一组递推公式。有了初始条件 ( ) px (0,0 )，就可以一步一步地使用前一时刻的坐标、动量 值确定下一时刻的坐标、动量值。这个方法是一步法的典型例子
由于在实际数值计算时h的大小是有限的,因而在上述算法中微分被离散化为差分形式来计算时 总是有误差的。可以证明一步法的局部离散化误差与总体误差是相等的,都为OH)的量级。在实际应 用中,适当地选择h的大小是十分重要的。h取得太大,得到的结果偏离也大,甚至于连能量都不守 恒;h取得太小,有可能结果仍然不够好。这就要求我们改进计算方法,进一步考虑二步法。 差分计算的二步法 实际上泰勒展开式的一般形式为 +)=10+∑r0)+- 6.2.8) 其中O")表示误差的数量级。前面叙述的欧拉算法就是取n=1
由于在实际数值计算时 的大小是有限的，因而在上述算法中微分被离散化为差分形式来计算时 总是有误差的。可以证明一步法的局部离散化误差与总体误差是相等的，都为 的量级。在实际应 用中，适当地选择 的大小是十分重要的。 取得太大，得到的结果偏离也大，甚至于连能量都不守 恒； 取得太小，有可能结果仍然不够好。这就要求我们改进计算方法，进一步考虑二步法。 h )( 2 hO h h h 差分计算的二步法： 实际上泰勒展开式的一般形式为 ( ) () ( ) ( ) 1 1 )( ! + = +=+ ∑ + i n n i i hOtf ih tfhtf . (6.2.8) 其中 ( ) n+1 hO 表示误差的数量级。前面叙述的欧拉算法就是取n = 1
现在考虑公式(6.2.8)中直到含h的二次项的展开(即取n=2),则得到 f(+h)=f( (6.2.9 f(t-h)=f()-h,+ (6.2.10) 将上面两式相加、减得到含二阶和一阶导数的公式 2=1[/(+b)-2()+y(- (6.2.11) df f(t+h)-f(t-h) d t 2h 6.2
现在考虑公式(6.2.8)中直到含 的二次项的展开（即取 h n = 2 ），则得到 ( ) () ( ) 3 2 22 2 hO dx fdh dt df htfhtf +++=+ . (6.2.9) ( ) () ( ) 3 2 22 2 hO dx fdh dt df htfhtf ++−=− . (6.2.10) 将上面两式相加、减得到含二阶和一阶导数的公式 [ () ()])(2 122 2 htftfhtf hdx fd −+−+= . (6.2.11) ( ) ( ) h htfhtf dt df 2 + − − = . (6.2.12)
令/()=x),利用牛顿第二定律公式F()=mx,公式6.2.1)写为坐标的递推公式 x(+b)=2()-x-b)+h2 (6.2.13) 公式(6.2.12)写为计算动量的公式得到 p()=m()=m()=[x(+b)-x-h) (6.2.14) 这样我们就推导出了一个比(6.2.6)和(6.2.⑦)更精确的递推公式。这是二步法的一种,称为 Verlet 方法。 当然我们还可以建立更高阶的多步算法,然而大部分更高阶的方法所需要的内存比一步法和二步 法所需要的大得多,并且有些更高阶的方法还需要用迭代来解出隐式给定的变量,内存的需求量就更 大。并且当今的计算机都仅仅只有有限的内存,因而并不是所有的高阶算法都适用于物理系统的计算 机计算。 Verlet算法是分子动力学模拟中求解常微分方程最通用的方法
令 ( ) = (txtf )，利用牛顿第二定律公式 ( ) 2 2 dt xd = mtF ，公式(6.2.11)写为坐标的递推公式 ( ) () ( ) ( ) m tF hhtxtxhtx 2 2 +−−=+ . (6.2.13) 公式(6.2.12)写为计算动量的公式得到 () () () [ ] ( )( ) htxhtx hm tmvtxmtp −−+=== 2 & . (6.2.14) 这样我们就推导出了一个比(6.2.6)和(6.2.7)更精确的递推公式。这是二步法的一种, 称为 Verlet 方法。 当然我们还可以建立更高阶的多步算法，然而大部分更高阶的方法所需要的内存比一步法和二步 法所需要的大得多，并且有些更高阶的方法还需要用迭代来解出隐式给定的变量，内存的需求量就更 大。并且当今的计算机都仅仅只有有限的内存，因而并不是所有的高阶算法都适用于物理系统的计算 机计算。Verlet 算法是分子动力学模拟中求解常微分方程最通用的方法