物理问题的计算机模拟方法(1)—分子动力学

硕士研究生课程
《物理问题的计算机模拟方法》讲义
适用专业:凝聚态物理、材料物理与化学、理论物理、光学工程
学时：30—40学时
参考教材：
1．[德]D.W.Heermann 著，秦克诚译，理论物理中的计算机模拟方法，北京大学出版社，1996。

2．[荷] Frenkel & Smit 著，汪文川等译，分子模拟—从算法到应用，化学工业出版社，2002。

3．M.P.Allen and D.J.Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1989.
4. A.R.Leach, Molecular Modelling: Principles and Applications, Addison Wesley
Longman, England, 1996.
5. [德] D.罗伯著，计算材料学，化学工业出版社，2002。

6. [英] B. Chopard & Michel Droz 著，物理系统的元胞自动机模拟，祝玉学，赵学
龙译，清华大学出版社，2003。

目录
第一章计算机模拟方法概论
1.1 序言
1.2 热力学系统物理量的统计平均
1.3 分子动力学方法模拟的基本思想
1.4 蒙特卡罗方法模拟的基本思想
1.5 元胞自动机模拟的基本思想
1.5.1 简要的发展历程
1.5.2 简单元胞自动机：奇偶规则
1.5.3 元胞自动机的一般定义
第二章确定性模拟方法—分子动力学方法（MD）
2.1 分子动力学方法
2.2 微正则系综分子动力学方法
2.3 正则系综分子动力学方法
2.4 等温等压系综分子动力学方法
第三章随机性模拟方法—蒙特卡罗方法（MC）
3.1 预备知识
3.2 布朗动力学（BD）
3.3 蒙特卡罗方法
3.4 微正则系综蒙特卡罗方法
3.5 正则系综蒙特卡罗方法
3.6 等温等压系综蒙特卡罗方法
3.7 巨正则系综蒙特卡罗方法
第四章离散性模拟方法—原胞自动机（CA）
4.1 引言
4.2 元胞自动机模拟
*4.3元胞自动机模拟的应用
第一章 计算机模拟方法概论
§ 1.1 序言
1．什么是计算机模拟？ Simulation Modelling 2．为什么要进行计算机模拟？ 3．常用的计算机模拟方法
确定性模拟方法：MD 模拟 Molecular Dynamics 随机性模拟方法：MC 模拟 Monte Carlo 离散性模拟方法：CA 模拟 Cellular Automata
§ 1.2 热力学系统物理量的统计平均
描述系统的坐标（自由度）：x(t)={x 1(t),x 2(t),…x N (t)} 系统的物理量：A (x (t)) 1．时间平均
dt t A t t A t
t t ⎰
-=
))((1
x ← 分子动力学（MD ）模拟 (1-1)
2．系综平均
⎰⎰
Ω
Ω
=>=
<x
x x x
x x d A d H f A Z
A )()( ))(()(1ρ ← 蒙特卡罗(MC )模拟 (1-2)
))((1
)(x x H f Z
=
ρ— 分布函数（几率密度函数） (1-3) ⎰Ω
=x x d H f Z ))(( — 配分函数 (1-4)
Ω—相空间
H (x )—系统的哈密顿函数 对于处于平衡态的系统，可以证明： ∞>=<A A 对于实际的有限时间内的平均，则有 A A >≈<
实际模拟的系统大小也是有限的：有限的粒子数N 或有限的系统限度L 对统计平均结果有影响。

§ 1.3 分子动力学(MD)方法模拟的基本思想
1．基本原理
系统：N 个粒子，体积V ，粒子质量为m
描述一个粒子运动状态的自由度：（r i , p i ） （p i =m v i ） 相空间：6N 维，相空间中的一点的坐标 X N =[r N , (m v N )] r N =(r 1, r 2, …, r N ), v N =(v 1, v 2, …, v N ) 粒子间的相互作用势：U(r N
)=U(r 1, r 2, …, r N )=∑<N
j i ij u )(r
决定系统相轨迹X N （t ）的运动方程：
⎪⎪
⎪⎪⎩
⎪
⎪⎪⎪⎨⎧==-∇== ( )0(N ..., 2, ,1)
( , 0）加上边界条件（周期性初始条件））（N N N i i i i X X i U dt d m dt
d r v
v r (1-5)
物理量A 的宏观值，由A(X N ) 的时间平均获得，即
⎰''=t
t d t A t t A 0
)]([1
)(X （离散情况：∑==k i i A k t A 11)(）
对于平衡态：)(lim t A A t ∞
→=
实际模拟时间总是有限的，模拟时间的长短可通过判断时间的增加对平均值的影响来确定，当继续增加时间带来的平均值得变化在允许的误差范围之内时，即可认为模拟足够长了。

2．计算步骤 运动方程：
)( , N i i i i U dt
d m dt d r v
v r -∇== 即
i N i i
U dt
d m F r r =-∇=)(22 (1-6)
或 m dt
d i i
/2
2F r = (1-7) 数值求解：用差分近似表示微分 （采用不同的差分格式，可得到不同的算法）。

用显示中心差分格式，将（7）式写为
2
2
2)/()]()(2)([t t t t t t dt
d i i i i ∆∆-+-∆+≈r r r r (1-8) 由（7）和（8）式可得：
m t t t t t t i i i i /)()()(2)(2F r r r ∆+∆--=∆+ (1-9)
第一步：由（9）式计算第i 个粒子在t +Δt 时刻的位置坐标。

要启动计算，我们必须要知道最初两点r i (0)和r i (Δt ) 第二步：对不同时刻
t = Δt , 2Δt , 3Δt , ……, L Δt (t 0 = 0) 计算物理量
A(r 1(l Δt ), r 2(l Δt ), ……, r N (l Δt )) (l = 1, 2, ……, N)
第三步：计算物理量A 的平均值
∑=∞→∆⋯⋯∆∆>=<L
l N L t l t l t l L A 1
21))( , ),( ),(A(1lim r r r
L 的大小由继续增大L 而<A>不变（或变化在误差范围内）来确定。

§ 1.4 蒙特卡罗（MC ）方法模拟的基本思想
1．基本原理
以正则系综（T, V , N ）为例
正则分布：s
E s e Z
βρ-=
1 正则配分函数：)( !1
3N N E N
d d d d e
h N s
p r =ΩΩ=
Z ⎰-β
系统能量：)(2)(2
N i
i
i
N p s U m U E E r p r +=+=∑
物理量：A(r N ) = A(r 1, r 1, …, r N ) 系综平均：
)(1 )(!1 )(!1 )(!1/)(/2/)(3/332
⎰⎰⎰⎰
⎰--
--=∑==Ω>=
<N kT
U N N
N
kT
m N
kT U N N N
N kT E N N s N
N d e A Q d e d e A Z h N d d e A Z h N d A h N A N i i i
N s r r p r r p
r r r r p r ρ （1-10）
⎰-=N kT
U N d e Q N
r r
/)( （位形积分） （1-11）
用MC 方法计算上述多维积分。

2．计算步骤
（1） 划分原胞
N 个粒子—3N 个空间自由度，3N 维空间
划分成s 个相等的原胞（s>>1）
注意：由于积分中不含动量，所以我们只需要在位置空间积分，而不需要在相空间中积分。

当系统的代表点落入第i 个原胞时，则认为系统处在状态i ，因此，s 为系统可能的微观状态数目。

于是，积分（10）和（11）可近似表述为 ∑->=
<i
kT
U i N
i e A Q A /1
(1-12) ∑-=i
kT U N i e Q / (1-13)
（2） 建立马尔可夫(Ma ρko в)过程（链）
将s 个状态看作一组随机事件
马尔可夫链：从状态i →j 状态j (εi →εj )的概率p ij ，只与εi 和εj 有关。

1=∑j
ij
p
，i =1, 2, …, s
若εi 经历n 步到达εj ，其概率表示为)
(n ij p ，存在极限概率
j n ij n u p =∞
→)
(lim ( j=1, 2, …, s) ∑=>j
j j u u 1 ,0
u j 为系统处在状态j 的概率。

于是，沿无限长的马尔可夫链，物理量A 的平均值可写为
)(∑∑∑=>=<i
j
ij i i
i i p A u A A (1-14)
选取 kT
U N
i i e
Q u /1-=，则（14）式为A 的正则系综平均值。

（3）抽样方法
采用怎样的抽样方法所构成的马尔可夫链能得到上述平均值？
粒子位置坐标：)
(i x αγ
粒子编号：γ =1, 2, …N 坐标的三个方向：α = 1, 2, 3 系统状态：i = 1, 2, …, s
给定粒子位置坐标的变化量δ (小于系统体积的限度)
给定系统的初态i ，随机选定4个随机数，其中三个ξα （α =1， 2，3），且 –1 ≤ ξα ≤1，一个表示粒子编号 γ =1, 2, …N ，由此随机确定粒子位置的变化：
δξααγαγαγ+=→)()()(i j i x x x （ 确保 δαγαγ≤-)()(j i x x ）
若i j U U ≤，则 γ 运动到新的位置，即系统由状态i 过渡到状态j ； 若i j U U >，则再选一个随机数ξ 4（0 ≤ ξ 4 ≤1），
若kT
U U i j e /)(4-->ξ，则粒子保留在原位置，不发生i →j 的跃迁； 若kT
U U i j e
/)(4--<ξ，则发生i → j 的跃迁。

由此进行下去，则形成一个马尔可夫链（或过程），此链的长度L (即粒子行走的步数，远大于s)，由所计算的物理量的平均值
∑=∞→>=<L
i i L A L A 1
1lim (1-15)
不再随链的加长而改变来确定。

由此得到的平均值即正则平均值。

一般来说，L 与N , V , T 有关，比如，N =32~108, L =3000~5000。

归纳起来，计算系统物理量的正则系综平均值的具体步骤如下： 第一步：给定系统的初始状态(粒子的初始位置)r i 和每一步的改变量δ；
第二步：选择四个随机数，其中一个代表粒子的编号i ( 1≤ i ≤ N )；另外三个表示粒子空间坐标的改变δx , δy , δz ( -δ ≤ δα ≤δ , α=1, 2, 3);
第三步：计算粒子i 的新位置i r '
),,(),,(z i y i x i i i i i z y x z y x δδδ+++='''='r 第四步：计算粒子在新旧两个位置系统的能量之差 ),...,,...,,(),...,,...,,(2121N i N i U U U r r r r r r r r -'=∆ 第五步：由U ∆的大小判断粒子i 是否从r i 运动到i r ': 若0<∆U ，则r i →i r '；
若0>∆U ，则再选一个随机数R(0 < R < 1), 如果 kT U e R /_∆<，则r i →i r '；
如果 kT U e R /_∆>，则r i 不变，返回第二步。

第六步：计算),...,,...,,(21N i A r r r r '
第七步：重复上述各个步骤，直到完成L 步为止，最后利用公式（15）计算A 的平均值。

3．粒子间相互作用势模型的选取 最简单的两种模型： （1）硬球模型
⎩
⎨
⎧<∞≥=)()
(0)(σσr r r u (σ 为硬球的直径) （2）L —J 势
126)(r
B r A r u -=
∑<=j
i ij N r u U )(),...,,(21r r r
§ 1.5 元胞自动机(CA)模拟的基本思想
元胞自动机：时间和空间都离散、物理参量只取有限数值集的物理系统的理想化模型
cellular automata 或 cellular automaton — CA
1.5.1 简要的发展历程
1．自繁殖系统
20世纪40年代，V on Neumann, 构造能解决非常复杂问题的计算机，设想模仿人脑的行为——寻求与生物过程无关的情况下自繁殖机理的逻辑抽象。

根据S. Ulam的建议，V on Neumann 在由元胞构成的完全离域的框架下处理这个问题，构造了一个完全离散的动力学系统——元胞自动机。

第一个自复制元胞自动机——由二维方形网络组成，有数千个基本元胞构成的自繁殖结构。

（1）一般机器只能构造比自己简单的客体，而采用自复制元胞机，可获得一种能产生新的、具有同样复杂性和功能的“机器”；
（2）V on Neumann 的元胞自动机规则具有所谓通用计算的性质，这意味着，存在一种元胞自动机的初始构型，该元胞自动机能产生任何计算机算法的解。

通用计算的性质指：用元胞自动机演化规则能够模拟任何计算机流程（逻辑选择器开关）。

2．生命游戏机
1970年，数学家John Conway 生命游戏机的概念，寻找能导致复杂行为的简单规则。

设想一个类似于棋盘的二维方形网格，每个元胞可能的状态是活（状态1）或死（状态0），其更新规则是：有三个活元胞包围的一个死元胞恢复为活元胞；由两个以下或三个以上活元胞包围的活元胞因孤立或拥挤而死亡。

结果表明，生命游戏机有出乎意料的丰富行为，从原“汤”中显示出来的复杂结构，演变发展成为某些特殊的技艺，例如，可能形成所谓的滑翔机——紧邻元胞的特殊排列，这些元胞具有沿直线弹道穿越空间运动的特性。

生命游戏机也是具有计算通用性的元胞自动机。

3．模拟物理系统
（1）20世纪70年代，Hardy, Pomeau 和Pazzis 建立了所谓的HPP格子气体模型，用以在质量和动量守恒的情况下在方形网格上模拟粒子的碰撞行为。

（2）1986年，Frisch, Hasslacher 和Pomeau 提出了著名的FHP模型，这是在六边形网格上模拟二维流体动力学的第一个严格模型——全离散计算机模型替代风洞试验。

HPP和FHP 通常称之为格子气自动机（LGA—Lattice Gas Automata）
（3）Ising 自旋动力学模型，20世纪80年代末，Vichniac 提出Q2R 规则。

（4）格子Boltzmann 方法与多粒子模型
格子Boltzmann 方法或模型（LBM ）：网格上定义的物理模型，在这个网格上与每个格位相关联的变量是平均粒子数，或具有一定速度粒子出现的概率。

该模型可以用均化或因子分解方法由元胞自动机动力学推导出来，或者自定义，而与特定的实现无关。

格子Boltzmann 模型保持了元胞自动机方法的微观水平解释，但忽略了多体的相关函数，但这种方法已经成为目前模拟物理系统中最有前途的方法之一。

在严格的元胞自动机方法与较灵活的格子Boltzmann 模型之间，有一种目前处于发展中的模型——多粒子模型。

这种模型保留量化状态的概念，但接受无限数值集，因此既保证了数值稳定性（与LBM 相反），又考虑了多体相关性。

在模拟物理系统时，大量的可能状态提供更多的灵活性，并产生小的统计噪声。

但多粒子动力学更难设计，且在数值计算上比格子Boltzmann 方法更慢。

1.5.2 简单元胞自动机：奇偶规则
简单元胞自动机的演化规则（20世纪70年代，Edward Fredkin 提出，定义在二维方形网上）：
网格的每一个格位是一个元胞，以其位置r =(i , j ) 来标记，其中i 和j 为行和列的标号。

函数)(r t ψ描述每个元胞在时间t 的状态，其值为0或1。

从时间 t = 0 的初始条件及网格上给定的构形值)(0r ψ开始，t = 1 时状态按下列步骤求得：
（1）对于每个格位r ，都计算出其位于东、南、西、北4个最近邻格位r '的)(0r 'ψ值之和。

应使系统在i 和 j 两个方向循环（如同在环面上），从而确定出所有格位的)(0r 'ψ计算值。

（2）如果这个和值为偶数，则新状态)(1r ψ为0（白色），否则为1（黑色）。

重复上述步骤，得到t = 2, 3, 4, …的状态。

这个元胞自动机的奇偶规则可表示为：
)1,()1,(),1(),1(),(1-⊕+⊕-⊕+=+j i j i j i j i j i t t t t t ψψψψψ
式中，符号⊕代表“异或”逻辑运算，也即模2和：1⊕1=0⊕0=0，1⊕0=0⊕1=1。

反复迭代这个规则时，可得到非常精美的几何图形。

1.5.3 元胞自动机的一般定义
定义：
（1）规整的元胞网格覆盖d 维空间的一部分；
（2）归属于网格的每个格位r的一组布尔变量Φ(r, t) = {Φ1(r, t), Φ2(r, t), …, Φm(r, t)}给出每个元胞在时间t = 0, 1, 2, …的局部状态；
（3）演化规则R = { R1, R2, …, R m}按下列方式指定状态Φ(r, t)的时间演化过程：
Φj(r, t+1) = R j[Φ(r, t), Φ(r + δ1, t), Φ(r + δ2, t),…, Φ(r + δq, t)]
式中r + δq指定从属于元胞r的给定邻居。

邻居：
二维元胞自动机的两种邻居：（1）V on Neumann邻居，有一个中心元胞（要演化的元胞）和四个位于其近邻东西南北方位的元胞组成；（2）Moore邻居，除了前面涉及的最近邻元胞外，还包括次近邻的4个元胞，共九个元胞。

还有一种有用的邻居称为Margolus邻居，将空间划分成2⨯2元胞的邻接单元块，这个规则对位于单元块内的位置即左上、右上和左下、右下很敏感。

三种邻居如下图所示。

边界条件：
（1）周期性边界条件；
（2）固定边界条件；
（3）绝热边界条件；
（4）映射边界条件。

备注：
确定型元胞自动机：演化规则确定，给定的初始状态将始终演化出同样的式
样。

概率型自动元胞机：演化规则包含一定的随机性，给定的初始状态可能演化
出不同的式样。

第二章 确定性模拟方法—分子动力学方法（MD ）
§ 2.1 分子动力学方法
用分子动力学方法模拟气体、液体或固体系统
1．分子动力学的描述形式： 哈密顿描述 拉格朗日描述 牛顿运动方程描述 2．划分MD 元胞
选取适当大小的 V=L 3，太大耗费计算时间，太小不能准确反映系统的性质。

目的：维持系统恒定的密度。

对平衡态系统，液体和气体原胞的形状无关紧要，但对晶体则不然。

3．周期性边界条件
)()(L A A n r r +=
目的：减少表面效应
4．相互作用势、最小影像约定和切断距离
∑∑∑<<+-+=n
j
i j i j
i ij N L u r u U )()(),...,,(21n r r r r r n = (n 1, n 2, n 3 )
⇓ ⇓
元胞内的粒子间 元胞内的粒子与影像
的相互作用 粒子的相互作用
最小影像约定(minimum image convention): 最小影像 L r j i ij n r r +-=min
元胞中的一个粒子只与元胞中的另外N-1个粒子中的每一个或其最近邻影像发生相互作用。

切断距离（cutoff distance ）：r c ≤ L /2 （切断半径）
上述约定相当于用r c 来切断位势，即r ij > r c 的相互作用忽略不计，代价是忽略了背景。

5．积分格式
∑∑<<=
∇-=-∇=j
i ij
ij j
i i N i i r r u U )()()()()(F r F
于是，运动方程可写为
∑<==j
i ij i i i r dt d m dt d )()( , F v
v r (2-1) （1）递推公式
在第一章中求解运动方程时，我们是直接求解的关于粒子位置坐标的二阶微分方程，得到的递推公式需要知道最初的两点位置才能启动计算，但在实际计算中，我们常常是给出最初的位置和速度，于是，我们可通过选取一定的差分格式，有运动方程(16)得到关于位置和速度的递推公式。

⎩⎨⎧→→++)
1()()
1()(n i
n i n i n i v v r r 或 ⎩
⎨
⎧+→+→)()()
()(h t v t v h t r t r 已知⎩⎨⎧==00)0()0(v v r r 递推公式的具体形式取决于差分各式的选取。

（2）时间步长
选择时间步长的原则：在保证计算精度的前提下，尽量节省计算时间。

实例：Ar 原子系统，L —J 势
时间步长取为 s h 14210~)(10~--约化时间步长(=100ps=10fs) （3）约束条件
保证能量、动量或角动量守恒。

（4）减小计算误差的技巧
数值计算不可避免有误差，与误差有关的因素主要有： 差分格式 时间步长 切断半径 最小影像约定 等等 6. 计算热力学量
dt t V t t A t t A t t N t ))(),(,)(()(lim 00v r ⎰'
∞
→'-'=
（1） 若系统的NVE 恒定，则
N V E A A >=< （微正则系综平均—microcannonical ensemble average ）
（2） 若系统的NVT 恒定，则
N V T A A >=< （正则系综平均—cannonical ensemble average ）
在实际计算中，往往还需要要求系统的总动量P =0或恒定，因整个系统不受外界的作用。

（3） 温度与能量均分定理
)3(2
1
c k N N kT E -=
N —总粒子数，N c —约束条件个数, 一般情况下，N >> N c P =0 ⇒ N c = 3 （4） 位势切断误差与修正
g (r ) — 对关联函数（pair correlation function ）
V
N
=
ρ — 粒子数密度 r r d g )(ρ—当原点位置上一个粒子时，在r 附近d r 内发现有一个粒子的几
率。

对气体或液体，)()(r ρρ=r （各向同性） 平均每个粒子的能量
⎰⎰⎰⎰⎰⎰≈+==Ω==∞
∞
c c
c
r r r V
V dr
r r g r u N
U
dr r r g r u dr r r g r u dr r r g r u d dr r r g r u d r g r u N U 0220
2
2
2)()(2)()(2)()(2)()(2 )()(21
)()(21ρπρπρπρπρρr
切断误差修正为：
⎰∞
=c
r c dr r r g r u U 2)()(2ρπ
6．计算机模拟的组织
（1） 初始化——给定粒子的初始位置、初始速度等；
（2） 趋衡——从初始转台开始，按运动方程要求的规律，从非平衡态趋向平
衡态；
（3） 投产——计算物理量的统计平均值。

§ 2.2 微正则系综分子动力学方法
微正则系综：NVE 恒定，且P = 0 （分子动力学方法的自然选择） 1．Verlet 算法
运动方程的显示中心差分格式
m t t t t t t i i i i /)()()(2)(2F r r r ∆+∆--=∆+ t t t t t t i i i ∆∆--∆+=2/))()(()(r r v
由可得到如下递推公式
⎩⎨⎧-=+-=-+-+ 2/)(/21
1211h r r v m
h F r r r n i n i n i n i n i n i n i （2-2） 此即Verlet 算法（A2）, 其特点是：
（1） 要启动此算法，必须已知两点10,i i r r ，所以又称为二步法； （2） n i r 和n i v 不能同步算出，第n+1步算出的第n 步的速度。

2．Verlet 算法的速度形式
2
2
2
2)(21)()( )()
(21)()()(t m
t t t t dt t d t dt t d t t t i
i i i i i i ∆+∆+=∆+∆+≈∆+F v r r r r r
2)()()()( )()()( )()()(t m t t t t t t t m t t t t t t m t t t t ∆+∆++=∆+⇒⎪⎩
⎪⎨⎧
∆∆++=∆+∆+=∆+F F v v F v v F v v （向后差分）
（向前差分）
由此可得如下递推公式：
⎩⎨⎧++=++=+++
2/)(2/1121m F F h v v m
F h hv r r n
i n i n i n i n i n i n i n i （2-3） 此即Verlet 算法的速度形式（A3）, 其特点是：
（1） 启动此算法，必须已知两点00,i i v r ；
（2）1+n i r 和1+n i v 同步计算。

（1+n i F 的计算需要知道1+n i r ） 此算法比A2算法更稳定。

3．趋恒阶段的能量调整
由于系统初始状态的能量离我们要模拟系统的能量（或温度）有一定差异，这就需
要在系统趋于平衡的阶段进行调整。

最常用的方法之一就是进行速度标度：
11++←n i n i v v β 标度因子：
∑-=2/12]/)33[(i ref mv kT N β （2-4）
T ref 为设定的系统的温度值。

能量设定值： r e f r e f kT N E )33(2
1
-= 能量参考值： ∑=22
1
i mv E
2β=E
E r e f
2)(2
1∑=i
r e f v m E β 要注意的问题：
（1） 一般不需要每一步都进行标度，可每隔若干步（比如50步）调整一次； （2） 当能量达到所设定的值后，停止标度，在以下的模拟时间内能量保持不
变。

4．实例
氩原子系统， N=256 L —J 势:
⎥⎥⎦
⎤⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛-⎪⎪⎭⎫
⎝⎛=6124)(ij ij
ij r r r u σσ
ε 作用力： ⎥⎥⎦
⎤⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛-⎪⎪⎭⎫
⎝⎛-⎪⎭
⎫
⎝⎛=∂∂-
=8
14221)(48)()(ij ij j i i
ij ij x r r x x x r u r F σσ
σε
取时间单位：s m 122
/12
103~48-⨯⎪⎪⎭
⎫ ⎝
⎛ε
σ
取长度单位：σ = 0.3405nm
取质量单位：m=6.63382⨯10-26 kg(氩原子质量)
⎥⎥⎦
⎤⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛-⎪⎪⎭⎫ ⎝⎛⎪⎪⎭⎫ ⎝⎛-=⎥⎥⎦
⎤⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛-⎪⎪⎭⎫ ⎝⎛⎪⎪⎭⎫ ⎝⎛-⎪⎭⎫ ⎝⎛=8
1428
1422148 2148)(ij ij j i ij ij j
i ij x r r x x m m r r x x r F σσσσσσεσσσσσσε
使用上述单位，可得到约化单位下的作用力表达式：
⎥⎥⎦⎤⎢⎢⎣
⎡⎪⎪⎭⎫ ⎝⎛-⎪⎪⎭⎫ ⎝⎛-=8*14*****11)()(ij ij j i ij x r r x x r F 标度因子：
2
/12
**2
/12
*2
/12
*2
2
2/12]16/)1[( ]16/)
1[( ]48/)1[3( ]/)33[(∑∑∑∑-=-=-=-=i i ref
i ref i ref v T N v kT N v m m kT N m v kT N ε
σε
σβ
约化温度： ⎪⎭⎫
⎝⎛==k T T k
T r e f r e f εε/*
K k
8.119=ε
取时间步长h = 2⨯10-14s = 20fs, 相当于约化时间0.064 约化温度：T* = 2.53 → 303K, T* = 0.722 → 86.5K 约化数密度：ρ* = N /L 3, N = 256
ρ* = 0.636 → L = 7.83 , ρ* = 0.83134 → L = 6.75
当 N = 64 时，ρ* = 0.83134 → L = 4.25, 所以，取r c = 2.5 是错的，因为要求 r c > L/2 = 2.13。

（习题3.5）
计算源程序见p129, 程序PL1
§ 2.3 正则系综分子动力学方法
正则系综：NVT 恒定，且P = 0 如何保持温度T 恒定？
1．速度标度方法 11++←n i n i v v β
∑-=2/12]/)33[(i ref mv kT N β
算法A4:
(1) 规定初始位置0i r 和初始速度0i v ；
(2) 计算⎩⎨⎧++=++=+++ 2/)(2/1121m F F h v v m
F h hv r r n
i n i n i n i n i n i n i n i ； (3) 计算∑-=2/12]/)33[(i ref mv kT N β； (4) 对速度进行标度：11++←n i n i v v β，返回（2）。

说明：
（1） 在前面的微正则系综的模拟中也对速度进行了标度，但其目的是使总能
量达到所需要得值，而且不必每步进行速度标度，可隔若干步标度一次，一旦总能量达到所设定的值，即可停止标度；而在正则系综的模拟中，每步都必须进行标度，以始终保持温度恒定。

（2） 通过一个广义位势引进能量涨落，连同细致耦合的一种特殊选择，会导
致速度标度机制。

（见书中的证明，公式3.59有误）
（3） 由于控制温度的反馈环内的时间延迟，温度将在一定程度上涨落。

2．随机方法（Anderson 热浴）
体系与一强加了温度的热浴相耦合，与热浴的耦合由偶尔作用于随机选择的粒子上的脉冲力表示。

在开始模拟前首先要选择系统与热浴的耦合强度，这一耦合强度由随机碰撞的频率ν 决定。

模拟步骤：
（1）规定初始位置0i r 和初始速度0i v ；
（2）计算⎩⎨⎧++=++=+++
2/)(2/1121m F F h v v m F h hv r r n
i n i n i n i n i n i n i n i ； （3）粒子i 在时间间隔h 内发生碰撞的几率为h ν，产生随机数ξ，如果ξ<
h ν, 则粒子i 经历一次随即碰撞，否则，返回（2）；
（4）粒子i 从温度为T 的麦克斯韦—玻尔兹曼分布中获得一个新的速度，返回（2）。

§ 2.4 等温等压系综分子动力学方法
1．等压等焓系综
N P H 恒定，P = 0 (P :压强，P ：总动量) V P E H E += (P E 为外压强，平衡时，P E =P ) V P E H E ∆+∆=∆
若系统与外界绝热，则V P E E ∆-=∆ ，即0=∆H （等焓过程）
要保持压强不变，体积则必须变化，为此，将体积V 作为一个新的动力学变量，构造一个新的拉氏函数：
V P V M U m U T U T V V E i i V
V -+-=-+-=∑2
22
1)(21 ),,,( r r r r L 引入标度坐标：
L r V r i i i //3/1==ρ 或 i i L r ρ=
i i
i i
i i L L L m
p v r
ρρρ ≈+===)( （0≈L ，体积变化缓慢） ⎪⎩
⎪⎨⎧∂∂=→→i i i i i m L r
m r ρρρ L 2 ， V M V → （共轭动量） V P V M L U m L V V E i
i -+-=∑2
2221)(21),,,( ρρρL ρ
系统的哈密顿量为：
V P V M L U m L V V H E i
i +++=∑22221)(21),,,( ρρρρ （2-5）
运动方程为：
i
i H m L dt d ρρ
∂∂-=)(2 （对粒子） （2-6）
V H V M dt d ∂∂-=)( （对活塞） （2-7）
将（2-5）代入（2-6）得
i
i
i
i i i i LF r r U U m L L m L =∂∂∂∂-
=∂∂-=+ρρρρ
22
L
L Lm F i i i ρ
ρ2-= V V dt dV dV dV dt dV dt dL L 3/23/13/13
1-==== ⎪⎪⎭
⎫ ⎝⎛==-V V V V L L
31313/13/2 所以，
⎪⎪⎭
⎫ ⎝⎛-=V V Lm F dt d i i i ρρ3222 （2-8） 将（2-5）代入（2-7）得
E
E P P P V
U V M -=-∂∂-= 所以，
M P P dt
V d E
-=2
2 （2-9） 由维里定理
∑∑∑<<+=
∙+=j
i ij ij i i
j
i ij
ij
F r m L K PV 312132 3
1322
2ρ F r
∑∑<+=
j
i ij
ij i
i F
r V
m L P 31312
ρ (2-10)
归纳起来：
⎪⎪⎪⎩
⎪⎪⎪⎨⎧+=-=⎪⎪⎭⎫ ⎝⎛-=∑∑<j i ij ij i i E i i i F r V m L P M P P dt V d V V Lm F dt
d 3131 3222222ρρρ 这些方程给出一个恒定的平均压强。

下面推出递推公式，利用泰勒展开保留到二阶项得：
⎪⎪⎩
⎪⎪⎨⎧++=++=++222122212121dt V d h dt dV h V V dt d h dt d h n n n n n n n n ρρρρ 将（2-8）和（2-9）式代入上式得：
⎪⎪⎩
⎪⎪⎨⎧-++=-++=++ )(21312121221E n n n n n n n n n n n n P P h M V h V V V V h m L F h h ρρρρ （2-11） 定义偏速度：
⎪⎩
⎪⎨⎧-+='-+='++ )(21312111E n n n n n n n n n n P P h M V V V V h mL F h ρρρ 递推公式可写为：
⎩⎨⎧'+='+=++++1111n n n n n n V
h V V h ρρρ 为了计算第n +1步的压强，需要知道第n +1的速度，为了解决这个问题，利用偏速度近似有：
⎪⎪⎩
⎪⎪⎨⎧-+'='-+'=++++++++++ )(2131211211111211211E n n n n n n n n n n P P h M V h V h V V h m L F h h h ρρρ （没有严格的证明） 压强的递推公式为：
∑∑<+++++++=j i n ij n ij i
n n i n n F r V m L P 111211131)(31ρ （2-12） 计算时，用偏速度1+'n ρ
代替1+n ρ 。

于是，可把上述算法表述为如下形式：
算法A5（NPH 系综）
（1） 规定初始位置和初始速度：0000,,i i i i r
r ρρ →； （2） 规定一个同所要求的密度不矛盾的0V ；
（3） 规定0V （比如，00
=V ）； （4） 由递推公式计算11,++n n V ρ；
（5） 计算偏速度11,++''n n V ρ
； （6） 计算111,+++n n ij n F r F ；
（7） 利用偏速度计算1+n P ；
（8） 计算1+n V
； （9） 计算1+n ρ。

将上述算法同速度标度结合起来，就可得到等温等压系综算法，即
算法A6（NPT 系综）
（1）-（9）和算法A5相同；
（10）速度标度11++←n n ρβρ。