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

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

摘要: 介绍了分子动力学模拟的基本原理及常用的原子间相互作用势, 如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(r

1

,r

2

,r

N

),式中，m

i

，r

i

分别为第i个原子的质量和位置。▽

i

=-∂/∂r

i

，

V(r

1,r

2

,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时刻的速度