第四讲 运动方程的积分算法

节省内 存？？
能量守 恒？？
积分运动方程的注意点
要想获得允许使用长时间步的算法，必须将 信息存储在粒子坐标的较高阶导数中。结果 是这需要更多的内存。对于一个通常的模拟 来说，这并不是一个严重的缺点。因为除非 处理很大的体系，存储这些导数所需的内存 与即便是在一个通常的工作站上所能获取的 总量相比来说仍是很小。
二阶对称算子
➢可以看出原子应力具有“能量密度”的量纲。 ➢包括了原子动量流和原子间作用力的贡献。 ➢是原子的一种力学“活性能”，反应了原子产生运动的潜在能 力。
原子应力
Voronoi几何构形的数学描述
原子应力
Lammps原子应力的定义
原子对势相互作用
原子键作用
键角作用
面弯曲作用
约束作用
二面角作用
能量守恒是一个重要的判据，但实际上需要 区分两种能量守恒，即短时间的和长时间的。 复杂的高阶算法通常在短时间内（如在几个 时间步内）有很好的能量守恒性，然而该方 法通常会有所不期望的特征，即时间较长时 总能量漂移。（不一定最好）
积分运动方程的注意点
最好有一个算法能同时在短时间和长时间准确地预
测所有粒子的轨迹。
第四讲 运动方程的积分算法
黄敏生
积分运动方程的注意点
• 何为一个好的运动方程积分算法？
计算速 度？？
但这并不是很贴切，因为花在积分运动方程上 的时间分率（相对于计算相互作用来说）很小， 至少对原子或简单分子体系是这样。
准确度对较大的时间步长来说更重要， 因为所能使用的时间步长越长，单位时 间内力的计算量越少。因此，这意味着 准确度？？ 采用允许用长时间步长的成熟算法是有 利的。（然而，保证准确度，不允许发 散）
计算 不 同 热力学状态的总能E和温度，可以得到内 能一温度曲线E(T )。这对于监测相变的发生非常有 用，相变发生时，该曲线会有跳跃。
原子应变
• 应变 是 表 示物体变形大小的测度，可以 有不同的定义方式。
• 应变是一个相对量，反映了物体相对初始 构型的变形程度。
• 应变又是纯几何量，与空间尺度变化无关。 • 因此宏观、纳观尺度下应变可以采用相同
RMS (Reservoir Modeling System保守系统的模拟) Energy deviation
(log/log scale)
时间步长
• 时间步长与研究对象、系统温度、所采用 的数值积分算法及势能函数有关。不存在 一个通用的时间步长值。
• 一般情况下，体系的温度越低，允许采用 的时间步长越大；而模拟较高温度时必须 采取较小的时间步长。
正加速度ac(t+ t)。定义预测误差
❖第三步根据加速度的预测误差对各预测量进行修 正：
Gear预测校正法
• 预测－校正法允许的时间步长比其它算法 长两倍以上
• 每个积分步内要计算两次体系势能，以得 到原子间相互作用力。
• 该算法的稳定差，能量波动较大，较verlet 算法占用更多的内存。
时间步长
• 一般来说，在探讨其力学行为时，我们需要研究 一些宏观力学概念的物理量，比如温度、能量、 压强、应力状态等，这就需要我们对分子动力学 计算出来的粒子数据进行分析。
• 统计力学是连接微观层次的物理量和宏观概念物 理量的桥梁，利用统计力学原理我们可以从系统 中单个粒子的运动学状态得到整个粒子系统的一 些性质。
• 一般认为，时间步长应小于原子振动周期 的十分之一，而通常原子振动周期的数量 级为0.1皮秒(10-12s)，即时间步长应选择在 飞秒级(10-15s)。
宏观物理量的计算
• 分子 动 力 学计算最终得到的是系统各个时刻的 相空间轨道(the phase-space trajectory)，包括任 意时刻所有原子的坐标和速度，这些都是微观原 子层次的物理量。
Too short - computation needlessly slow Too long - errors result from approximations
Just right - errors acceptable, maximum speed
时间步长
过长的时间步
时间步长过大，原子的作用 力急剧改变。 误差逐渐累计，导致结果发 散.
两个 Argon (氩)原子 在两个不同时间步长 ts的模拟。 图中画出的是计算模拟值与理论值的差。
误差与时间步长的关系
Circles: Verlet Squares: Gear 4th order Triangles: Gear 5th order Diamonds: Gear 6th order
Gear预测校正法
预测校正法是分子动力学模拟中的常用算法之一， 其基本思想是Taylor展开，这种算法包含三个部 分。 第一步，利用泰勒展开预测下一时刻的位置及其一 阶、二阶、三阶导数：
式中v、a、b、c分别是位置矢量r的一阶、二阶、三阶和四阶导数。
Gear预测校正法
❖第二步根据新的原子位子rp，计算受力以及修
• 对标准Verlet算法进行改进而得到的蛙跳 算 法 (Leapfrog methods) 。 相 比 标 准 Verlet算法，它有两个优点。
• 一、包含显式速度项。 • 二、计算量稍小。 • 缺点
原子的位置与速度计算不同步，这就意 味着在确定位置时，不能同时计算体系的 动能，给模拟过程带来不便。
的定义方式
原子应力
宏观应力
反映了单位面积上作用力的 大小，是关于面积的强度量。
原子应力
离散原子系统的原子应力是 关于体积的强度量。实际上， 原子应力只是形式地沿用了 应力的概念，具有与宏观应 力完全不同的特征。
原子 应 力 表明一个原子与周围原子相互作用的强弱程度
原子应力
动量 势能
Hamilton 自由能
谢谢
标准Verlet算法
• Verlet提出的Verlet算法在分子动力学中 应用最为广泛。
由Taylor公式展开有：
+
位移
标准Verlet算法
由Taylor公式展开有：
-
速度
不出现在 算法中。
在典型分子动力学模拟中，只有原子的初始位置和初 始速度是给定的，而在verlet算法中，计算下一步的 原子的位置需要前两步的信息。
Leap-frog算法
速度Verlet算法
速度Verlet算法
Velocity-Verlet算法不仅可以获得相同精度的原 子位置和速度量，给出了显式的速度项，而且在 每步积分中只需要存储一个时刻的状态变量，模 拟稳定性好，允许采用较大的时间步长，计算量 适中，因而在分子动力学方程的积分算法中得到 了最广泛的应用。
那么如何采用verlet算法计算第二个时间步的原子位 置呢？
标准Verlet算法
• 解决方法一是，在第二个时间步中，把该 步力看成常量，且使用普通运动学方程:
ri (t) ri (0) ri (0)t
• 解决ห้องสมุดไป่ตู้法二是
• 流程
标准Verlet算法
标准Verlet算法
加速度项
Leap-frog算法
不存在
MD模拟所研究的所有体系，体系的轨迹穿过相空
间（即对于用由所有粒子坐标和动量所跨过的6N维
空间）时敏感地依赖于初始条件。
这意味着两个初始靠得很近的轨迹随时间的演绎
将会显著分开。
李雅普诺夫不稳定性
MD获得的轨迹在某种意义上与真实的轨迹相接近。 MD的目标并不是精确地预测一个已知初始条件的体系轨迹将会发生什么（卫 星轨道预测）。而对统计预测感兴趣。在MD中，统计预测是足够精确的。
宏观物理量的计算
• 系统 的 物 理性质是系统中粒子坐标和速度 的函数，对于任意一个时刻宏观概念物理 量A，定义为
统计平均
宏观物理量的计算
能量
• 势能 部 分 可以按势函数计算，对于不同的原子势函数表 达式有不同的计算公式。
• 动 能计 算 公式为
温度
温度T直接与粒子动能相关，即著名的均匀分
布公式，每个自由度赋予kBT/2的能量, N个 粒子的总自由度为3N，故动能为
Lammps中原子应力没有除于体 积，是能量的量纲。（e.v）
原子应力
原子应力与宏观应力的关系
Lammps中原子应力
• 宏观应力实际上是原子应力在一定程度上的统计平均。 • 如果要计算体系中某个区域（由region定义，可以是整个
模拟盒）所围成的“块”的应力，在Lammps只需将该区 域里的所有原子的单原子应力值加起来，再除以这个区域 的体积即可，无须进行单个原子体积的计算。