数值模拟导论 -第十七、十八讲
分子动态学
分子动态学
分子动态学是计算经典的多体系统的平衡性和非平
衡性的一种技术
*连续粒子的原子运动遵循经典的机械定律（牛顿定
律）
参考资料：
1.液体计算机模拟， M.P. Allen & D.J. Tildesley,
Clarendon, Oxford, 1987.
2.了解分子模拟：从计算到应用， D. Frenkel and B.
Smit, Academic Press, 1997.
3.Moldy手册
Moldy
z免费的分子动态学模拟资料压缩包可以在
CCP5程序库
（ http://www.ccp5.ac.uk/librar.shtml） 的
Moldy中下载，该网页中还有各种非常有用的
信息和分子模拟工具。
zMoldy 简单易用，并附有一个很好的手册可
参考帮助，我希望布置的家庭作业可以借助
Moldy来完成。
我们为什么要研究分子动态学？
nF= =数字密度， 流量
类似真实实验
1.让我们去研究和理解材料性能以便我们建立模型
2.告诉我们答案是什么，什么时候模型不存在
质量守恒：
实例：不收敛方程
极限：
z该方程只有给定n和F的关系才能求解
z实验或分子性能表明在不同条件下
方程不收敛！
Kn 1
Kn 1
0.1 Kn 1
线性梯度连续定律
z大梯度
z不平衡气态流动
－冲击波
－小量流动 （高努森数字流动）
－稀薄流动（低密度） （高努森数字流动）
高努森数字流动（气体）
zKn为分子平均自由路径比率，表示纵向尺度
z分子平均自由路径分子两次碰撞间隔间的平均运动距离
z碰撞趋向于回复平衡
时，粒子运动是扩散的（接近平衡）
时，粒子运动是扩散的（远离平衡）
时，物理性能变化不定（难于建立模型）
的分解
z当
z当
z当
其它高努森数字流动
z小尺度流动（在空气压力接近60纳米时
的平均自由路径）
z真空技术（低压状态）
z在高空气流密度低（分子碰撞几率小）
z较长的平均自由路径
z典型的高努森数字流动实例
实例：
再进入交通工具的高空飞机制动
SMA-HPC ?2003 MIT
简要介绍统计力学
统计力学为物质的微观描述（如分子的位置和速度）
和诸如肉眼可见的压力、密度、温度等宏观描述之间提
供了理论联系。
这种联系时通过统计方法和总体均值的概念获得的，
一个整体是由很多个同一的的系统（ M）在相同的宏观
条件不同的微观初始条件下逐渐演变而成的集合体
令 为系统在状态 下的分子，
i则 可以理解为在状态 下找到全体分子的可能性
k
z宏观特性（可见的）为加权平均值
或者为其连续极限
z统计力学的一个最基本的结果就是一个系统在固定温度 T
时，能量为 E的系统达到平衡状态的可能性，
确定，其中
z分子方法类似于试验方法，直接测量 即可，不用计算ρ (Γ )。
为玻耳兹曼常数。它由
z对于非平衡系统，问题的求解则变成计算ρ (Γ ).
表示一个给定的全体系统，任一可见的
－由于系统的初始状态和其它时间状态有一一对应的关
系，多个系统初始状态和一个系统不同时间状态是等同的。
可通过多个系统
值的加权平均来测定。
但是，实际上我们不需要用很多个系统来做实验。我们通
常是观察一个系统的不同时间段来测定，这是因为我们采用
了各态历经的假设。
－各态历经假设允许我们将求多个系统的均值转化为求同
一系统不同时间段的均值，这看起来是一个在大多数的情况
下都非常好用的假设。
一个简化的 MD程序结构
tt=
最后
tt>
平衡
z初始化：
－读入运行参数（初始参数，时间步数，密度，粒子数，
时间步长）
－生成或读取初始配置（粒子的位置和速度）
－计算粒子受力
－对运动方程积分
－若
z结果输出
，系统达到平衡状态
z循环时间步增量，直到
运动方程
牛顿定律
因而
其中 包含了三个实体的交互影响
Lennard-Jones势能
ε
σ
常用的最简单的势能为 Lennard-Jones 势能
典型的应用是模拟简单液体和固体
为井深（能量）
为交互长度比例
时强烈排斥
时电势最小
时弱相吸
Lennard-Jones势能（ U）和受力（ F）为间距（ r）的函数
约化单元
zLennard-Jones粒子的运动方程积分的合适
的时间步长是多少？
zLennard-Jones分子组成的液体 /气体的密度
是多少？
zLennard-Jones液体的汽化温度是多少？
32
10
?
分子密度
温度
压强
时间
z对于所有 Lennard-Jones分子，结果有效
必定有误！） z点误差放宽：（
z在这些单元中，分子具有物理意义
积分算法
积分算法必须要求：
a)快速，消耗内存较少
b)允许采用大时间步
c)尽可能精确地复制经典的轨道线
d)满足冲量和能量转化并且可逆
e)简单易于程序编写
积分算法的讨论
a)不是特别重要，但是模拟中耗费最大的是计算
力
b)很重要，因为对于一个给定的模拟，时间步越
长力需要赋值的次数越少
c）不太重要，因为没有哪个数值计算能长时间提
供精确的解（在轨迹线附近程指数偏离），由于
MD是一种在给定宏观常数下获取 “所有初始状态 ”
的平均时间的方法，因此冲量和能量转变更重要
d）非常重要（见 C）
e）重要，在没有增速情况下无需复杂计算
Verlet算法
在 MD的最初几十年内 Verlet算法是最受欢迎的算法之一
公式推导：泰勒级数在 t时刻展开
优点：
1） .程序非常的简单紧凑。
2） .非常好的能量守恒工具（得益于时间 ——可逆性）
3）时间可逆性
4）局部误差
1.对速度的处理非常笨拙
之前需要解
缺点：
a)在求得之前需要解
2.或许对切断误差比较敏感，这是因为在公式
中一个小数与微分方程中的两个大数相加。
b)速度的误差
改进的Verlet 算法
Beeman运算法则
与 Verlet算法有相同的适应性
速度的计算比 Verlet算法更加精确
预校正算法
a)预知
虽然这些方法都已经很精确，但是 MD仿真自身的属性
决定了这些方法并不适应他们。这是因为任何的预言 -校正
算法，如果没有把邻居考虑在内都是不可靠的。
基本构成：
时刻的位置，速度，加速度。
b)从新的位置估算速度和加速度（如果力由速度来决定）
c)应用新的加速度来校正位置，速度，加速度
d)返回（ a）步
这里我们可以用改进的 Beeman算法来证明这种方法用来处
理速度依赖于力（以后讨论）的情况是可以的。
如果出现速度不依赖于力的情况将超出 Beeman算法的讨论范围。
周期性边界条件
?周期性边界条件非常普
遍：简化表面作用
本部分内容可以参考：流体的仿真
计算，作者 M.P. Allen 和 D.J.
Tildesley,，牛津科学出版社，
1987。
?今天的计算机可以很容易的处理 N>1000的周期性边界条件
受限制的小系统
?平衡条件不受影响
?长波是不可能的
?在程序的运行中，粒子与其最相近的粒子的镜像会相互影响
宏观量一般都被定义为微观量的平均值
如果系统处于非平衡状态，在长度（空间，时间）
变化很小的时候我们可以通过求其平均值来将这一工具
定义为空间和时间的函数。
（宏观速度）
仿真的开始
?需要知道系统中所有粒子的位置速度的初始值
?初始化已知系统的密度，温度，数量
?由于高阶非线性粒子的相互作用，要想在完全的任
意状态下开始仿真几乎不可能成功
?我们可以用平衡分布的方法非常安全的初始化速度
的自由度
如果粒子的位置任意初始化，有很 大的可能性至少有一对粒子
离的很近，这将导致能量的畸变。
由于动能这一附加量 :
?典型的初始化可以用在流体仿真上
导致对每一个粒子速度都有相对独立的概率分布
- 一个近似的随机方针数据结构避免了粒子的重
叠（作为一个例子我们可以参考 moldy手册）
- 晶体在平衡点融化是系统仿真的一部分
平衡
由于在不恰当的状态我们使用了典型系统来模拟，势
能往往不是转化成热能，就是热能随着平衡的到来而消
失。为了阻止温度的浮动，我们需要用到下面的方法：
1）.变速度尺度法
这里Td是需要的速度并且
是瞬时速度。
在仿真过程中，这是一种最简单也是最直接的控制
温度的一种方法。
使用温控器也是一种保持温度稳定性的一种方法，
它是通过把温度恒定作为一个约束加入到运动方程中去
的。温控器经常被描述成为 “动力系统约束 ”
2） .温控器。
大范围力、切断误差和相邻粒子列表
?Lennard-Jones势能在
处衰减的非常的快
?尽管如此，当在 时，相邻粒子之间相互作用的数量将会增加。
?在某一距离 处由于节省计算耗时，相互作用被忽略
?对于大范围力来说，在计算大多数的敏感量（表面张力）
时往往将
近似为
。
?典型的流体力学的计算经常将 近似为
?静电的相互作用经常需要特殊的方法（多重扩展， Ewald求和
法） -参考 moldy手册
?虽然系统的特征（如：状态方程，传送系数，潜热，弹性不变量）受
如果需要，我们需要测量新的
?最简单的切断法就是分段求解。
由于 是间断的，所以这个方法没有什么优势。并没有保存能量
的影响，
。
?下面是根据切断量和潜在的改变进行修正后的方法
?即便是一个很小的切断值，计算耗时也会与
—— Verlet列表法
这是因为我们需要检测所有的成对的分离量。
成比例，
?我们可以使用 “相邻粒子列表 ”来简少计算耗时。
——单元索引法
Verlet相邻粒子列表法
?保持一个扩展的相邻粒子列表
计算出相邻的例子之间的相互作用
?我们可以通过测试 10-20个时间步来选择
。
可以计算
本部分的内容来源于流体的计算机
仿真。作者 M.P. Allen 和
D.J. Tildesley，牛津科学出版社，
1987。
(需要很多的存储空间 )
?
从而不需要在每一个时间步
单元索引法
?在每个方向上将仿真系统分解成 个亚单元（这里 ）
?这里我们只是在分离点考察亚单元（保守算法）
例子：
如果亚单元的尺寸大与分离点那么，如果我们要研究单元 13那么我
们只是需要考察单元 7， 8， 9， 10， 12， 13， 14， 17， 18， 19。
?计算两个维度的耗时是 4.5
由原来了 变成了现在的
。
?计算三个维度的耗时为 13.5 ，而不是原来的
（当 被恰当的重新定义）
。
约束法
?牛顿的分子动力学保存了系统的能量，体积和
粒子的数量。
?典型的工程物理学系统就是要在一定的压力，
温度和频繁的与外界交换中（也就是开放系统）
进行处理。
?对这样的系统进行仿真被提了出来。
*在系统处于平衡状态时这种方法可以对其提供
一种准确的统计学描述。在系统远离平衡时，我
么还没有找到确切的物理学的定理来对系统进行
描述。因此，我们需要转换，尤其需要像变尺度
法这样自然的方法
恒温仿真
?牛顿方程所保存的能量和温度是可变的。
-变速度尺度法 （未经加工）
?三种主要方法的类型
?同样我们可以以次来考察压力。我们也重视希望在
多变的系统中来进行恒定的压力计算。
?在大多数实际影响的计算中我们总是愿意去描述温度。（在
现实中，大容量的系统，例如大气层，系统的各种影响相互
作用并且保持温度恒定）
-扩展系统法 （添加一个或者更多的自由度到大容量系统中）
运动方程：
其中 f是动力学可变量并且我们可以通过以下公式获得：
g=系统的自由度
Q=大容量系统的集合（可调整的参数控制平衡动力学）
这种方法被称为 Nosé-Hoover温控器
-约束法 （运动方程被强制放在相空间的一个平面上）
在这种情况下，约束是
这产生了下面的运动方程。
是一个拉格朗日乘数（不是一个动力学变量）
在 Moldy手册中我们可以看到一个关于本方法和
压力约束法的介绍
纪录“速度 -力”之间的关系。
作业讨论
?应用（或者修改）控制文档 control.argon。 你可以用这
一文件来进行简短系统的分析模拟。
? Control.argon调用了系统的 argon.in系统文件。将粒子的
数量增加到至少 1000个（减小了小系统的影响和噪音）
?在 Control.argon文件中密度的单位是
?我们可以通过增加粒子的数量和 /或平均时间间隔来减少
噪音。
?确保温度没有发生明显的漂移。你或许会经常处于正在重
新调节或者停止调节的过程中。确保在重新调节的过程中
没有进行取样
?通过设置 ，你忽略了放射性波动函数的计
算使的输出文件更加的易读。
采样控制文件
#
#采样控制文件。不同于代码产生的文件
#
#密度 =1335kg/m
3
#初始温度 =120k
#每 2个时间步长对速度进行一次调解
#在时间步 5000处停止改变速度的尺度
#总共运行 25000个时间步长
#每 500个时间步长打印一次结果
#每 500个时间步长进行一次循环。
#在第 5001时间步开始求平均值
#在 10000时间步求平均值
#每一事件步 0.01ps
#亚单元法的亚单元尺寸
#应用精密级来计算力
#3 * sigma
#在第 2501步开始放射分区函数的计算
