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

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

摘要: 介绍了分子动力学模拟的基本原理及常用的原子间相互作用势, 如Lennard-Jones势; 论述了几种常用的有限差分算法, 如Verlet算法; 说明了分子动力学模拟的几种系综及感兴趣的宏观统计量的提取。

关键词: 分子动力学模拟; 原子间相互作用势; 有限差分算法;

分子动力学是一门结合物理，数学和化学的综合技术。分子动力学是一套分子模拟方法，该方法主要是依靠牛顿力学来模拟分子体系的运动，以在由分子体系的不同状态构成的系统中抽取样本，从而计算体系的构型积分，并以构型积分的结果为基础进一步计算体系的热力学量和其他宏观性质。

从统计物理学中衍生出来的分子动力学模拟方法(Molecular Dynamics Simulation, 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=-▽i V(r1,r2,r N),式中，m i，r i分别为第i个原子的质量和位置。▽i=-∂/∂r i，V(r1,r2,r N)为体系所处的势。

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个粒子）中的每个粒子或其最邻近影像粒子发生相互作用。实际上，这个约定就是通过满足不等式条件r c

4、势函数：

MD模拟结果准确与否的关键在于对系统的原子之间相互作用势函数的选取，总的来说，原子（或分子）之间的相互作用势的研究进展一直很缓慢，在一定程度上制约了MD方法在实际研究中的应用，原子间的势函数的发展经历了从对势到多体势的过程，对势认为原子之间的相互作用是两两之间的作用，与其它原子的位置无关，而实际上，在多原子体系中，一个原子的位置不同，将影响其它原子间的有效相互作用。所以，多体势能更准确地表示多原子体系势函数。

4.1、对势：