第六章-new分子

1 第六章 分子动力学模拟方法 6.1 引言 分子动力学方法利用牛顿力学基本原理，通过积分运动方程产生体系的连续构型，从而得到所有原子的轨迹（即体系中各粒子的位置和速度依时间而变化的轨迹），进一步从轨迹中计算得到各种性质。 复习牛顿运动定律三要素： 1. 物体总是连续地作匀速直线运动，除非一个外力作用于它。 2. 力等于动量的比率，即

madtdvmdtmvddtdpf)(

3. 对每一个作用有一个相等并相对抗的反作用。 对于一个质量为im位置矢量为iX的原子，其在t时刻的加速度为

22xii

i

dXF

dtm （6.1）

而粒子受力等于势能对坐标的偏导数 xii

dVFdX

由方程6.1得： 1ixi

i

dXdFdtdtm

dtFmdvttxiivvi00

1

)(00,ttmFvvixiii 从而得： 0,0)(iixiivttmFv （新速度） 又 0,02)(iixiivttmFdtXd dtvttmFdXiixii0,0)( dtvdtttmFdXttittixiiXXii000,0,0)( 2

)()(2100,200,ttvttmFXXiixiii 从而得 0,00,20)()(2iiixiiXttvttmFX （新的位置） 经过)(0ttt变化后，力xiF也发生改变。 所以新的力xiF又被计算，从而又导出新的位置和新的速度。 在应用牛顿定律时，由于物体运动时受力的情况不同，因此运动的状态也有所不同。讨论三种不同情况： 第一种情况：最简单情况，在二次碰撞间，无力作用质点上，

即 0xiF )(0tt

所以 0,0,0)(iiixiivvttmFv （常速度）

0,00,0,00,20)()()(iiiiixiiXttvXttvttmFX （位置有iv t改变）

第二种情况，在二次碰撞间，粒子受到一个恒力。比如一个带电粒子在均场中运动，即受力：

EqF 第三种情况，粒子受力，依赖其它粒子的相对位置，所以受力偶的影响。

6.2 利用简单模型的分子动力学 第一个利用简单硬球模型报道分子动力学模拟工作的是Alder and Wainwright 1957完成的。 1. 模型的特点： a. 碰撞间球作均速直线运动； b. 完全弹性碰撞，球中心之间的距离等于碰撞直径; c. 成对势，采用如图6.1中所示的形式： 左图为硬球势模拟

当212,1rrr时发生碰撞。 右图为方井势模拟 当 r ＞σ2 时， V（r）= 0 ; r ＜σ1 时， V（r）→; r 在 [σ1 ，σ2 ] 之间时， V（r）=V0 。

2. 硬球模型计算的具体步骤

a. 确定碰撞发生时,碰撞球的相邻对。 b. 计算在碰撞时间，所有球的位置。 3

c. 测定碰撞后两个碰撞球的新的速度。 d. 重复由a至结束。 3. 新速度计算运用线性动量的守恒原理。 即：

221122112121vmvmvmvmpppp。

简单相互作用模型有许多不足，但它仍提供许多有用的性质。

6.3 具有连续势的分子动力学 在分子间相互作用的实际模型中，质点上的力随位置变化而改变；或随其它质点(与其互相作用)位置改变，力也改变。 利用连续势模拟的第一个模型是氩原子[Rahman 1964]，后来他又完成第一个液体分子（水）的模拟[Rahman和Stillinger 1971]。 因此Rahman等人在动力学的方法学上作出了重要的贡献。 在一个连续势的作用下，所有质点运动被偶连在一起，引起一个多体问题，不能用解析法解。在这种情况下，运动的方程只能利用有限差分方法来积分。

6.3.1有限差值法 有限差值方法的基本思想是整体的折散成许多小的段，每一个段通过固定的时间t被分隔。在时间t构型中每一质点上的总力等于与其它质点相互作用的矢量加和。 采用有限差值方法的分子动力学计算可以按照以下步骤进行： 1、 设定粒子的初始位置和速度 2、 根据粒子的位置计算每个粒子的受力 3、 根据粒子的位置、速度和受力计算粒子的新位置和新速度 4、 更新粒子的位置和速度，然后回到第2步 分子动力学牵涉到的几个关键性步骤，按重要性排列依次是积分方法（第3步）、初始化（第1步）和力的求算（第2步）

利用有限差值方法积分运动方程有许多算法，其中若干的方法常用在分子动力学计算中，所有算法都假定位置和速度、加速度（动力学性质）能够以泰勒级数展开被近似：

)(241)(61)(2

1)()()(432tcttbttatttvtrttr（6.2）

)(61)(2

1)()()(32tcttbtttatvttv（6.3）

c(t)t2

1tb(t)a(t)t)t(a2δδδ（6.4）

)()()(ttctbttb（6.5）

v: 速度，位置相对时间的一级微商。 a: 加速度，二级微商。 b. 三级微商。 一、积分方法 所谓积分方法就是根据粒子的位置、速度和受力计算δt时刻之后粒子的新位置和新速度的算法。 4

下面介绍几种算法： 1． Verlet法又称中心差分法。[Verlet, L., Phys. Rev. 1967, 159, 98; 1968, 165, 201.]

将ittr和ittr在t 时刻作Taylor展开

Verlet算法利用在时间t时的位置和加速度，及前一步长)(tt时的位置)(ttr,计算tt时的新位置。具体表达如下： a(t)t2

1tv(t)(t)rt)t(r2δδδ（6.6）

-)t(at2

1tv(t)-(t)rt)t(r2δδδ（6.7）

相加这两个方程给出： )()()(2)(2tatttrtrttr（6.8）

相减两个方程除以t2得 tttrttrtv2/)]()([)(（6.9）

tt 2tt tt t 3tt

PositionPotentialForceTime itF 

ittF 



2ittF

ittr 

itr ittr 2ittr 

3ittr

iVtr 

iVttr 

2iVttr

Verlet法的算法图示 说明： （1）Verlet法是分子动力学中采用最普遍的算法。

（2）计算精度： 关于有限差分法的计算精度详见附录D。微分方程用差分来求，显然其误差来自Taylor展

开，如（3.1.1-1）式等。Verlet法在（3.1.1-3）式中弃去的项是4tO，故Verlet法的精度为4阶。 （3）Verlet法对计算机存贮的要求属中等。它计算过程中需要存贮的量有itr、

ittr和加速度ita。

（4）Verlet法的缺点： 5

a、从（3.1.1-3）式可见，等式右边为：（大项）－（大项）＋（小项），故不利于准确度。 b、Verlet法中不出现速度项，故如有需要，需另外计算速度itv。

c、从图3.1-1可见，Verlet算法还需要初始条件。如0t时需用itr。解决办法：从（3.1.1-1）式得到00iiittrrv。 （5）后来从Verlet法发展出来几种变种。蛙跳法就是其中的一种。

* 蛙跳法 [Hockrey, R. W. Methods Comput. Phys. 1970, 9, 136.] 蛙跳法也是经常采用的算法。令t为模拟中的时间步长。将速度/2inttv和/2inttv分别在ntt时刻处Taylor展开：

22321222!2inininindtdtttttttdtdtvvvvO （3.1.2-1）

22321222!2inininindtdtttttttdtdtvvvvO （3.1.2-2）

两式相减得到： 322inininittt

ttttmFvvO （3.1.2-3）

ttrttrttv/)]()([)21( （6.10）

同样，再将位置inttr和intr分别在/2nttt时刻处Taylor展开： 2232/2/21222!2inininindttdtttttttttdtdtrrrrO

（3.1.2-4） 2232/2/21222!2inininindttdttttttttdtdtrrrrO

（3.1.2-5） 两式相减得到：

3/2ininintttttttrrvO （3.1.2-6）

所以，整个体系所有粒子的位置、速度都可以用图3.1-2所示的算法求得：