正则系综
平衡系综
等温等压系综
等压等焓系综
系统原子数N，压强P，焓值H=E+PV保持不变。 在模拟中较少见。
系综的控温
• 温度调控机制可以使系统的温度维持在给 定值，也可以根据外界环境的温度使系统 温度发生计 系综，即调温后各粒子位形发生的概率可 以满足统计力学法则。 直接速度标定法 Berendsen温控机制 Gaussian温控机制 Nose-Hoover温控机制
A.原子速度的初始化
• 1.为使模拟尽快达到平衡，分子初始速度的 分布应该尽量接近真实情形。 • 2. 采用近似的 Maxwell-Boltzmann 统计分布 来赋予原子的初始速度是比较合理的。能 够使得系统尽快弛豫。
x,[0,1]范围内的随机数
0，2Pi随机数
速度的高斯分布
A初始条件二
• 1.在满足以上温度的条件下，必须保证系统 净总动量为零。 • 2.另一种获得初始条件的方法是选取模拟过 程某一时刻的原子坐标和速度。 • 3.分子动力学模拟经常分成不同的物理阶段 进行，上一个模拟过程结束时的原子位置 和速度就可以作为下一次模拟的初始条件。
Verlet近邻表
网格近邻表
近邻表
Verlet近邻表
• 最早由Verlet提出
首先通过建立链表记录各原子周围( 首先通过建立链表记录各原子周围 (截断 半径:r 半径:rc)区域内的所有原子， 区域内的所有原子，然后通过链表 计算各分子所受到的作用力。 计算各分子所受到的作用力。 为了计算原子1 为了计算原子1 所受的作用力。 所受的作用力。假定一个 球形的区域， 球形的区域 ， 半径为r 半径为 r1 ， 且 r1>rc 。 (rc 为截 断半径) 断半径 ) ， 并给截断半径r 并给截断半径 rc 区域内的分子建 立一个链表。 立一个链表。 当截断半径r 当截断半径rc范围内的分子没有离开r 范围内的分子没有离开r1球 形区域时， 形区域时，只需要根据链表中的分子， 只需要根据链表中的分子，即可 计算出分子1 计算出分子1所受到的作用力。 所受到的作用力。 如果截断半径r 如果截断半径rc范围内的分子离开了r 范围内的分子离开了r1， 球形区域， 球形区域，则需要建立新的链表， 则需要建立新的链表，即在计算 过程中， 过程中，每隔一定时间， 每隔一定时间，更新列表。 更新列表。
分子动力学模拟的基本过程
黄敏生
基本步骤
A.原子位置的初始化
• 建立分子动力学模拟过程的首要问题和第 一步是确定分子体系的初始条件。 • 两种方式，一是采用实验数据，二是借助 各种理论模型得到分子结构的几何参数， 如面心立方(FCC)模型等。
A.原子位置的初始化
• 1.无论采取哪种方法，给定分子结构的空间坐 标都不一定处在分子力场最稳定的位置，即 各原子并非处在平衡态，造成体系的能量比 较高。 • 2.要进行一个不施加载荷的弛豫过程，使得系 统达到稳定的平衡状态（共轭梯度法）。 • 3.在这个过程中，系统从人为的初始构形转变 成真实初始构形，势能减小并达到稳定。 • 4.初始条件最好与真实构形类似，Fcc BCC, 固体结果影响较大，气体影响较小。
Nose-Hoover温控机制
• 根据以上Hamiltian推导真实系统的运动方程：
摩擦项
摩擦系数率
控压技术
• 系统常压条件 压力诱导下的相变 • 对于各向同性材料，瞬时压力定义为：
系统的压力可以通过改变体积来调节。 直接体积标定法 Berendesen控压机制 Aandersen控压机制
直接体积标定法 • 在一定时间步 内将 系统体 积自 动 重 分 配 以 达到目标压力值。 • 计算中引入材料的体积模量 B，并在每个时 间步内计算系统应力，由体积模量和应 力，得到系统的体积变化为
Time RList Update interval 5.78 12.50 26.32 43.48 83.33 100.00 N=256 3.33 2.24 2.17 2.28 2.47 2.89 N=500 10.00 4.93 4.55 4.51 4.79 5.86
占用了一定量的计算机内存。分子数较 时，此方法具有一定的优势
Berendsen温控机制
Gaussian温控机制 • 在Gaussian法中，给每个原子施加一个摩擦 力，摩擦力的系数依赖于当前系统温度和热 浴温度之间的差值。 • 若摩擦力系数为正，则表示系统温度高于热 浴温度，需对系统进行冷却；若系数为负， 则表示加热系统。
摩擦力系数
摩擦力项
Gaussian温控机制 • 保持温度不变就是保持体系的总动能不变。 也 就是体系 内部 力 所做功 （ 功率 为零(Ff).v=0） • 根据上式，可以得到摩擦力系数为：
力场的截断
在分子动力学中， 在分子动力学中，出于计算上的考虑， 出于计算上的考虑，力场的截 即在某一范围内力场是有效的， 断是必须的， 断是必须的 ，即在某一范围内力场是有效的 ，因 此会导致一些计算上的困难。 此会导致一些计算上的困难。
势函数直接截断： 势函数直接截断：
V (rij ) - V c V (rij ) = 0
非周期边界条件
• 分子体系的模拟并不是都使用周期性边界条件，在 很多情况下，如溶液中沉淀的分子团簇、蛋白质分 子、病毒分子、材料的表面等并不需要周期性边界 条件。 • 可以根据分子体系所处的外界环境对非周期边界上 的粒子施加一定的限制。 • 例如，边界上的原子设计为位置固定的，就可以形 成刚性边界。（原子始终不动） • 对边界上的原子施加一定荷载或考虑边界上原子与 外界环境之间的作用力，就形成阻尼边界
近邻表
网格近邻表
• 在计算过程中，每次更新原子位置时，将跨越本 网格边界的原子从网格中删除，将其插入到相应 的邻近网格中。 • 与邻域列表法相比，此方法不占用多余内存，在 进行大规模分子系统模拟时，此方法可以明显的 减少计算量。
近邻表
网格近邻表
单位的无量纲化
• 在模拟中涉及很多浮点和指数运算，为提 高计算效率，往往将温度、密度、压力等 类似量表示成无量纲的形式。
Berendsen温控机制 • 又称Berendsen外部热浴法。 • 其基本思想是假设系统和一个恒温的外部 基本思想是假设系统和一个恒温的外部 热浴耦合在一起， 热浴耦合在一起，通过热浴吸收和释放能 量来调节系统的温度， 量来调节系统的温度，使之与恒温热浴保 持一致。 持一致。 • 对速度每一步进行标定， 对速度每一步进行标定 ， 以保持温度的变 化率与热浴和系统的温差(Tbath-T(t))成比例。 成比例。 每一步温度变化： 每一步温度变化：
• 所得期望温度为体系初始温度。
程序编码简单的优点，但消除局域的相关运动。
Nose-Hoover温控机制 • 该方法可以把任何数量的原子与一个热浴 耦合起来，可以消除局域的相关运动，而 且可以模拟宏观系统的温度涨落现象。 • 是等温系综统计力学研究领域的一个里程 碑。 基本思想
引入一个反映真实系统与热浴相互作用的 广义变量s(无量纲量)，将真实系统与热 浴作为统一的扩展系统来考虑。扩展系统 的哈密顿量为：
式中Li是计算元胞在i方向上尺寸，Si为i方向上 的应力，n是一个控制参数。调节元胞趋向合适 的尺寸，即可使压力达到期望值。
Berendesen控压机制
Aandersen控压机制 • 假定系统与外界“活塞”耦合，当外部压强不 能补偿系统内部压强时，“活塞”运动引起系 统均匀地膨胀或收缩，最终使得系统压强 等于外部压强。 • Andersen方法具有重要的意义，后来的各 种压力控制方法基本都是基于Andersen思 想发展起来的。
直接速度标定法 • 引入速度标定因子 λ 或 • Ndf为系统的总自由度数。Tc为目标温度；K 为标度前体系的总动能。标度时
实际分子动力学模拟中，并不需要对每一步的速度都 进行标定，而是每隔一定的积分步，对速度进行周期 性的标定，从而使系统温度在目标值附近小幅波动。
直接速度标定法 • 直 接速度标定法的 优点 是原理 简单 ， 易 于 程序编制. • 缺点 是模拟系统无法和 任 何一个统计力学 的系综对应起来。 • 突然的速度标定引起体系能量的突然改 变， 致 使模拟系统和真实结构的平衡态 相 差较远。
直接速度标定法
• 系统温度和粒子的速度直接相关，可以通过调整 粒子的速度使系统温度维持在目标值。
m v (t ) T (t ) = ∑ i =1 k B N f
• N系统粒子的个数。Nf=3*N-3为体系的自由度数（总 动量固定）。波尔兹曼常数kB= 1.38×10-23 J/K。
N
2 i i
绝对温度T与体系的总动能密切相关。 • 温度是一个统计平均的量。对于单个的原子， 其温度？？？？
边界条件
• 分子动力学模拟中，只有足够的粒子数量， 才能准确的描述材料的宏观性能。 • 为了减小计算规模，人们引入了周期性、固 定、全反射等边界条件。目前常用的边界条 件包括周期性边界条件、对称边界条件和固 定边界条件。
1.周期性边界(Periodic boundary condition)
像胞元
No 2.60 2.70 2.90 3.10 3.43 3.50
近邻表
Verlet近邻表算法
近邻表
Verlet近邻表—周期性
方法（a）
方法（b）
近邻表
网格近邻表
这种方法的思想是将研究对象看成一个方盒，将这一方盒划分 为M×M×M个单元（cell），每个单元的边长必须大于势函数 的截断半径。 与单元13内的原子间距小于截断半径的其它原子必然在单元13 与其邻近单元7，8，9，12，13，14，17，18，19共9个单元内。 由于每个单元内原子数为Nc=N/M2，因此对每个原子只要计算 9Nc个原子间距，对整个原子系统就要计算9NNc个原子间距。 对三维结构则要计算27NNc个原子间距（Nc=N/M3）。 关于原子间距的计算量就与微结构的尺度即原子系统的原子数 N成正比。