分子动力学模拟方法的基本原理与应用

分子动力学模拟方法的基本原理与应用摘要:介绍了分子动力学模拟的基本原理及常用的原子间相互作用势,如Lennard-Jones势;论述了几种常用的有限差分算法,如Verlet算法;说明了分子动力学模拟的几种系综及感兴趣的宏观统计量的提取。关键词:分子动力学模拟;原子间相互作用势;有限差分算法;分子动力学是一门结合物理，数学和化学的综合技术。分子动力学是一套分子模拟方法，该方法主要是依靠牛顿力学来模拟分子体系的运动，以在由分子体系的不同状态构成的系统中抽取样本，从而计算体系的构型积分，并以构型积分的结果为基础进一步计算体系的热力学量和其他宏观性质。从统计物理学中衍生出来的分子动力学模拟方法(MolecularDynamicsSimulation,MDS),实践证明是一种描述纳米科技研究对象的有效方法,得到越来越广泛的重视。所谓分子动力学模拟,是指对于原子核和电子所构成的多体系统,用计算机模拟原子核的运动过程,从而计算系统的结构和性质,其中每一个原子核被视为在全部其他原子核和电子所提供的经验势场作用下按牛顿定律运动。它被认为是本世纪以来除理论分析和实验观察之外的第三种科学研究手段,称之为“计算机实验”手段,在物理学、化学、生物学和材料科学等许多领域中得到广泛地应用。科学工作者在长期的科学研究实践中发现，当实验研究方法不能满足研究工作的需求时，用计算机模拟却可以提供实验上尚无法获得或很难获得的重要信息；尽管计算机模拟不能完全取代实验，但可以用来指导实验，并验证某些理论假设，从而促进理论和实验的发展。特别是在材料形成过程中许多与原子有关的微观细节，在实验中基本上是无法获得的，而在计算机模拟中即可以方便地得到。这种优点使分子动力学模拟在材料研究中显得非常有吸引力。分子动力学模拟就是用计算机方法来表示统计力学，作为实验的一个辅助手段。分子模拟就是对于原子核和电子所构成的多体系统，求解运动方程(如牛顿方程、哈密顿方程或拉格朗日方程)，其中每一个原子核被视为在全部其它原子核和电子作用下运动，通过分析系统中各粒子的受力情况，用经典或量子的方法求解系统中各粒子在某时刻的位置和速度，以确定粒子的运动状态，进而计算系统的结构和性质。该模拟技术主要涉及粒子运动的动力学问题,与蒙特卡罗模拟方法(简称MC)相比,分子动力学是一种“确定性方法”,它所计算的是时间平均,而MC进行的是系综平均。然而按照统计力学各态历经假设,时间平均等价于系综平均。因此,两种方法严格的比较计算能给出几乎相同的结果。经典的分子动力学方法是Alder等于1957年提出并首先在“硬球”液体模型下应用，发现了由Kirkwood在1939年根据统计力学预言的“刚性球组成的集合系统会发生有液相到结晶相的转变”。后来人们称这种相变为Alder相变。Rahman于1963年采用连续势模型研究了液体的分子动力学模拟。1972年Less等发展了该方法并扩展了存在速度梯度的非平衡系统。1980年Andersen等创造了恒压分子动力学方法。1983年Gillan等将该方法推广到具有温度梯度的非平衡系统，从而形成了非平衡系统分子动力学方法体系。1984年Nose等完成了恒温分子动力学方法的创建。1985年针对势函数模型化比较困难的半导体和金属等，Car等提出了将电子论与分子动力学方法有机统一起来的第一性原理分子动力学方法。1991年Cagin等进一步提出了应用于处理吸附问题的巨正则系综分子动力学方法。20世纪80年代后期，计算机技术飞速发展，加上多体势函数的提出与发展，使分子动力学模拟技术有了进一步的发展。1、分子动力学的运动方程：分子动力学模拟的出发点是假定粒子的运动可以用经典动力学来处理，对一个由Ｎ个粒子构成的孤立体系，粒子的运动由牛顿运动方程决定，也就是：mid2ri/dt2=-▽iV(r1,r2,rN),式中，mi，ri分别为第i个原子的质量和位置。▽i=-∂/∂ri，V(r1,r2,rN)为体系所处的势。2、运动方程的数值积分：计算机模拟方法的基点是利用现代计算机高速和精确的优点,对几百个以至上千个分子的运动方程进行数值积分有许多不同的积分方法,它们的效率和方便程度各异问题基本上就是用有限差分法来对二阶常微分方程进行积分常用的有以下几种方法：2.1、Verlet算法［1］：Verlet算法是在60年代后期出现的，对扩散分子的质心运动的积分是最稳定的也是最常用的数值方法。它运用t时刻的位置和加速度以及t时刻的位置来预测t＋δt位置，其积分方案，以三阶Taylor展开为基础，由以下方程给出：r(t＋δt)＝2r(t)-r(t-δt)+δt2α(t)这里，为简单计，省略了i。速度可按微分的基本法则得出：V(t)=[r(t+δt)-r(t-δt)]/2δt。这种算法的优点是占有计算机的内存小，并且很容易编程。但它的缺点是位置r(t+δt)要通过小项(δt2)与非常大的两项2r(t)和r(t-δt)的差相加得到，这容易造成精度损失。并且从式中可以看出，这种算法不是一个自启动算法，新位置必须由t和t-,δt时刻的位置得到。2.2、Gear的预侧-校正算法：这种算法分为三步来完成：首先，根据Taylor展开，预测新的位置、速度和加速度。然后，根据新的计算的力计算加速度。这个加速度再由与Taylor级数展开式中的加速度进行比较，两者之差在校正步里用来校正位置和速度项。这种方法的缺点就是占有计算机的内存大。2.3、“蛙跳”（Leap-frog）算法[2]：Hockey提出的Leap-frog算法是Verlet算法的变化，这种方法设计半时间间隔的速度，即：r(t+δt)=r(t)+δtv(t+δt/2)，v(t+δt/2)=v(t-δt/2)+δtα(t)。t时刻的速度由下式给出：v(t)=[v(t+δt/2)+v(t-δt/2)]/2。这种算法与Verlet算法相比有两个优点：（1）包括显速度项；（2）收敛速度快，计算量小。这种算法明显的缺陷是位置和速度不同步。除了上述提及的几种方法外,还有Beeman算法、Rahman等。3、周期性边界条件和长程力：即使是使用现代的巨型计算机，MD方法还是只能用于粒子数大约是几百到几千的系统。这就引起一个问题：用这样少量的粒子，如何来模拟宏观体系？为了解决这个问题，引入了周期性边界条件[3]。采用这种方法，模拟体系实际上是由基本单元（也称为模拟计算元胞）在各个方向上重复叠合而成。但在模拟中只需保留基本单元，所有其它单元与基本单元由平移对称性关联。在处理粒子之间的相互作用时，通常采用“最小影像”约定。这个约定是在由无穷重复的MD基本模拟计算元胞中，一个粒子只与它所在的基本元胞内的另外N-1个（设在此元胞内有N个粒子）中的每个粒子或其最邻近影像粒子发生相互作用。实际上，这个约定就是通过满足不等式条件rcL/2来截断位势（rc为截断半径，L是元胞的边长）通常L的数值应当选得很大，以避免有限尺寸效应，但这样会增大计算量，同城采用对相互作用势的修正来近似处理。4、势函数：MD模拟结果准确与否的关键在于对系统内的原子之间相互作用势函数的选取，总的来说，原子（或分子）之间的相互作用势的研究进展一直很缓慢，在一定程度上制约了MD方法在实际研究中的应用，原子间的势函数的发展经历了从对势到多体势的过程，对势认为原子之间的相互作用是两两之间的作用，与其它原子的位置无关，而实际上，在多原子体系中，一个原子的位置不同，将影响其它原子间的有效相互作用。所以，多体势能更准确地表示多原子体系势函数。4.1、对势：在分子动力学模拟的初期，经常采用的就是对势。对势可以分为间断对势[1]和连续对势，而连续对势主要有以下几种：Lennard-Jones（L-J）势、Born-lande（B-L）势、Morse势和Johnson势等，其中，L-J势是为描述惰性气体分子之间的相互作用而建立的，因此它表达的力比较弱，描述材料的行为也比较柔韧，也可以用来描述过渡金属原子之间的相互作用。B-L势是用来描述离子晶体离子之间的相互作用的。Morse势和Johnson势多用于描述金属原子之间的相互作用。对势虽然简单，得到的结果往往也符合某些宏观物理规律，但其缺点是必然导致Cauchy关系，所以，对势实际上不能准确地描述晶体的弹性性质。4.2、多体势：多体势是在20世纪80年代初期开始出现的，1984年Daw和Baskes首次提出了原子嵌入（EAM）势[4]。此势的基本思想是：把晶体的总势能分成位于晶格点阵上的原子核之间的相互作用对势和原子核镶嵌在电子云背景中的嵌入能（多体相互作用）两部分，其中，对势和多体作用势的函数形式往往根据经验选取。基于EAM势的势函数还有很多种，这些多体势大多用于金属的微观模拟，此外，还有许多形式的多体势函数形式，1987年Jacobsen等人在等效介质原理的基础上提出的另一种多体势函数形式，由于其简单、有效，也得到了广泛地应用。5、分子动力学模拟的系综：平衡态分子动力学模拟是在一定的系综瞎进行的，经常用的平衡系综是NVT或NPT系综。在这两种系综中，牵涉到控制温度和压力的几种技术，分别如下：5.1、控温方法：（1）速度标定：系统的温度与动能存在如下关系：𝐸𝑘=∑𝑚𝑖|𝑣𝑖|22=(3𝑁−𝑁𝑐)𝐾𝐵𝑇2𝑁𝑖，式中，N是原子数，Nc是约束数，KB是Boltzmann常数，vi是原子i的速度。由于系统的温度和动能存在这样的关系，所以一种最简单和最直观的方法是直接对速度进行标定。这种方法的基本思想是：如果t时刻的温度是T(t)，速度乘以因子λ后，温度的变化为ΔT=(λ2-1)T(t)，其中，λ=√()，T为所控制体系的温度。（2）Berendsen恒温槽方法：这种控温方法假设所模拟的体系与一个恒温槽连载一起，则两者之间就可以通过热交换而使模拟体系达到恒温的目的，方法如下，定义一个参数λ：λ=√()，式中τ表征系统与恒温槽之间的热交换速率，Δt为MD的时间步长，那么通过vnew=λvold校正即可保持体系的温度在T0附近震动，而参数（通常取为0.1~0.4ps）则可用于控制这个震动幅度。（3）Nose-Hoover方法：这种方法是通过改变模拟体系的Hamiltonian来实现控温的，因而有更强的物理意义，其基本思想就是在Hamiltonian加入一个假想的项来代表一个恒温源，具体做法如下：𝐻=∑𝑚𝑖𝑣𝑖22𝑉(𝑟)𝑄ζ22𝑔𝐾𝑙𝑛𝑆，式中，S和ζ分别是假想项的坐标和动量。这样体系的微分方程就变为：vi=dri/dr，ai=-(dV/dr+miviζ)/mi，d𝑑ζdt=(∑𝑚𝑖𝑣𝑖2𝑔𝐾𝐵)Q。式中，g为体系的自由度，Q为一个可调参量，表征着假想项的质量，T为温度。在这三种控温方法中，速度标定是一种非真实的物理效应。但这种方法可以使系统很快达到平衡，在经典分子动力学方法中，这是一种比较常用的控温方法。Berendsen控温方法是通过系统和恒温槽进行热交换来控制温度，此方法的优点是非常简单，应用起来非常方便。Nose-Hoover控温方法是基于统计力学而提出来的一种控温方法，当系统和恒温槽进行热交换时，在系统中粒子出现的几率遵从统计力学规律，所以这是一种真实的物理效应。5.2、控压方法:用MD方法研究压力的诱导相变和结构重构，在等压下模拟比在等体积下更容易实现。在等压模拟下，可以通过改变模拟元胞的三个方向或一个方向的尺寸来实现体积的变化，类似于温度控制的方法，也有许多方法用于压力控制，在实际应用中，有以下两种较为常用的压力控制方法。（1）Berendsen方法：这种控压方法的基本思路与基本控制方程均与他的控温方法类似，如下：μ=[Δt(P𝑃0)γ]13，式中γ是一个可调参数，𝑃0为所控制的压力。然后，在每一次的MD迭代中，粒子坐标x，y，z均用μ相乘，得到新的坐标，即可实现对压力P的控制。该方法的缺点是可能导致在很长的模拟时间步长内，保持着起伏。（2）Parrinello-Rahaman方法：这种方法是通过改变体系的Lagran