赤字率王延颋 2012年12月12日
10. 分子动力学（Molecular Dynamics ，MD ）算法 分子动力学算法通过计算每个质点上所受的合力并求解多体牛顿方程来模拟多体体系的运动过程。通过第一性计算得到质点间作用力的方法目前只可以模拟到皮秒量级。实际应用中主要还是把粒子看作经典粒子，从而模拟时间可以达到纳秒至微秒量级。 10.1. 基本的MD 算法 10.1.1. 算法流程
力可以从作用势得到：()()d
f r V r dr
=-
，注意r 是在周期边界条件（PBC ）作用之后。
10.1.2. 牛顿运动方程的数值积分
不用考虑计算耗费的时间，关键是精度和长时间稳定性。以下是几种常用的运动方程积分方法。
a ） Verlet 算法：一类通用算法的基础。
泰勒展开
()()()()()()()()()()()()()()()()()()()()()3
243
242
42
2
23!
23!
222f t t r t t r t v t t t r O t m f t t r t t r t v t t t r O t m f t r t t r t t r t t O t m
f t r t t r t r t t t r t r t t a t m
∆+∆=+∆+∆++∆∆-∆=-∆+∆-+∆⇒+∆+-∆=+∆+∆⇒+∆≈--∆+∆=--∆+∆ (10.1)
瞬时速度可以从位置估计出来：
()()()
()22r t t r t t v t O t t
+∆--∆=
+∆∆
(10.2)
b ） Velocity Verlet 算法： 性能非常好的算法，应尽量使用，但是不能直接与Nosé-Hoover 热耦配合使用。
()()()()()()2222
f t a
r t t r t v t t t r t v t t t m +∆=+∆+∆=+∆+∆ (10.3)
()()()()()()()122
f t t f t v t t v t t v t a t t a t t m +∆-+∆=+
∆=++∆+∆⎡⎤⎣⎦ (10.4)
c ） Leap Frog 算法：定义了一半时间的速度，适合于配合某些力场算法使用。 定义
()()()()22r t r t t t v t t r t t r t t v t t --∆∆⎛⎫-≡
⎪∆⎝
⎭+∆-∆⎛⎫+≡
⎪∆⎝
⎭
(10.5)
则可以得到
()()222t r t t r t v t t t t v t v t a t ⎧∆⎛⎫+∆=++⋅∆ ⎪⎪⎪⎝
⎭⎨
∆∆⎛⎫⎛⎫⎪+=-+⋅∆ ⎪ ⎪⎪⎝
⎭⎝⎭⎩ (10.6)
d ） 预估-校正算法：时间不对称。 泰勒展开
()()2233
232!3!
r r t r t r t t r t t t t t ∂∂∆∂∆+∆=+∆+++
∂∂
∂ (10.7)
定义
()()()()()012222
33
33
2!3!x t r t r
x t t t
r t x t t r t x t t ≡∂≡
褐煤干燥∆∂∂∆≡∂∂∆≡∂
p308(10.8)
用t 时刻的x 值预估()n x t t +∆：
()()()()()()()()()()()()()()predicted 00123predicted
1
123predicted
223predicted 33233x t t x t x t x t x t x t t x t x t x t x t t x t x t x t t x t ⎧+∆=+++⎪+∆=++⎪⎨+∆=+⎪
⎪+∆=⎩
(10.9)
因为由牛顿方程可以从()0x t t +∆的值算出实际加速度2
corrected
x ，可以对第二项进行修正
corrected predicted 222x x x ∆=-
(10.10)
其它项可以修正为
corrected predicted 2n n n x x c x =+∆
(10.11)
对于给定的算法阶数，选定的常数n c 可以使得精度和稳定性均衡。
理论上，应该递归计算预估和校正项使结果自洽。实际计算中，因为每步递归都要进行
耗时的力的计算，更好地提高精度的方法是减少步长。
MD 的主要特征：1）模拟真实实验；2）时间平均而非系综平均，验证了Boltzmann 分
布；3）步长有限制；4）可用于非平衡过程；5）计算力和势能。
10.2. 正则系综的MD 算法
基本的MD 算法实现的是微正则系综，为了实现MD 对正则系综的模拟，需要在牛顿
积分后引入热耦（thermostat ）对系统温度进行调节。
10.2.1. Isokinetics Thermostat (Evans & Morriss)
每一步或几步把系统温度重置为环境温度T ：
2B 31
,22scale i i i i i
Nk T m v v v λλ=⇒= =∑ (10.12)寡核苷酸探针
问题：系统温度没有涨落，不是正则系综。可以正确给出只和位置相关的量。
10.2.2. Berendsen Thermostat
对温度的校正并非一步到位，而是给定一个系统与热耦的耦合强度T τ：
12
011T t T T λτ⎡⎤⎛⎫∆=+-⎢⎥ ⎪⎝⎭⎣⎦
(10.13)
等效于
0T
T T dT dt τ-= (10.14)
实际体系在正则系综下的温度是有涨落的。
证明：正则系综下体系的速度分布服从Maxwell-Boltzmann 分布：
()3
22
2B 1exp ,223p f p k T mv m m ββπ⎛⎫⎛⎫
=- = ⎪ ⎪
⎝⎭⎝⎭
(10.15)
其中m 是粒子质量，v 是粒子速度，p mv =是粒子动量。根据上式有
()223m
儒林外史的讽刺艺术p dp p f p β
=⋅=
⎰
(10.16)
()2
4
415m p dp p f p β⎛⎫
=⋅= ⎪⎝⎭
⎰
(10.17)
动能的方差为
2
22
2
242
2
2
2
22315233p
m m p p m p
p
σβββ⎛⎫⎛⎫- ⎪
⎪-⎝⎭⎝⎭=
==⎛⎫
⎪⎝⎭
姚国志
(10.18)
从而温度的方差为
(
)
2
22
2
2
2
42
222
2
2
22
4
22
2
112
3T T T T
N p N N p p N p N
p
p p
N N
p σ-=
+--=
-=
=
(10.19)
10.2.3. Andersen Thermostat
用随机碰撞来模拟热浴对系统的影响。如果对一个粒子的两次连续碰撞不相关，则在t 时间内被碰撞的次数为Poisson 分布：
()();exp P t t λλλ=-
(10.20)
使用Andersen Thermostat 实现MD 模拟正则系综的算法流程为：
1） 计算t ∆的常规MD 运动积分；
2） 随机选择一组粒子与热浴进行碰撞。每个粒子被选中的几率为t λ∆；
3） 如果某粒子被选中，则从温度为T 的Maxwell-Boltzmann 分布中随机产生一个速度赋予
该粒子，其它粒子不变。
可以证明，使用Andersen thermostat 的体系确实是对正则系综的采样。问题在于系统角动量不守恒，模拟的体系随时间旋转，不方便数据处理。
10.2.4. Nosé-Hoover Thermostat
基本思想是扩展系统的拉格朗日量，加入额外的、虚构的位置和速度自由度，使得扩展的体系是微正则系综，而实际体系是正则系综。
拉格朗日量为
k p L E E =-
(10.21)
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议