• 之后，此方法被M.J. Gillan等人推广到了具有温 度梯度的非平衡系统，从而构造并形成了非平衡 MD方法体系。
MD方法的发展史
• MD方法真正作为材料科学领域的一个重要研究方 法，开始于恒压MD方法(1980)和恒温MD方法 (1984)的建立及在应用方面的成功。
• 1985年人们又提出了将电子论和分子动力学方法 有机统一起来的所谓Car－Parrinello方法，即第 一性原理MD方法。它不仅可以处理半导体和金属 的问题，同时还可应用于处理有机物和化学反应。
周期性边界条件
假设现实物质中的 一部分原子(通常为 102－105个)，被取 出配置在所谓基本 单元的箱中，由于 基本单元周围的原 子、分子变成表面， 从而不同于本来要 处理的体状态。为 了防止这种情况， 在基本单元周围配 置其复制品(图3-4)。
图3-4
周期性边界条件
• 对于在基本单元周围边界的原子、分子所受到的 作用力，不仅有来自基本单元之内部的原子或分 子的作用力的贡献，还要考虑来自其近邻复制单 元的原子或分子的作用力的贡献。这样，把来自 距粒子某一距离（截断距离）内的j粒子的贡献 （这里与粒子j是否处在基本单元内或复制单元内 无关）求和。给出力的方法称为最小镜像距离法 （Minimum Image Distance）。
• 为了检验不同基本单元大小是否会引入“人为效 应”，必须用不同的基本单元尺寸做计算，若结 果能收敛，则尺寸选择是合适的。
边界条件
• 边界条件的问题比较复杂，因所考虑的具体情况而定，概 括地说，作为处理原子、分子团簇的边界条件可进行推广。 例如，适用于块体状态的三维周期边界条件，可扩展处理 表面重构和表面吸附、单层膜相变的二维周期边界条件， 以及处理异质晶体界面的边界系统等，同时还开发了模拟 非平衡的方法。诸如针对分子束外延（MBE）的边界条件， 人们进行了很多研究。
ij

4ij

ij
r
12



ij
r
6




(3-2-3)
式中，r 是原子间距， ij 是结合强度参数， 是表示ij 原子半径的参数。
对(3-2-1)可用数值积 分法求解，其数据处 理流程图见图3-2
图3-2 MD数据处理流程图
启动计算
设定坐标、速度初始值 计算作用在原子上的力Fi
巨正则系统MD方法
3.2 MD方法计算初步
• 在计算机出现以前，作为根据原子间相互作用力 等微观信息了解多原子或分子团的结构、性质的 方法，所采用的是基于统计理论的数学解析法。 然而，原子间相互作用力稍微复杂一些，不用说 求解统计理论严格方程解，就是进行数值求解也 是一件很困难的事。
• MD方法就是数值求解多体系统的确定性运动方程， 并根据对所求结果进行统计处理，决定粒子的轨 迹，从而给出物性预测和微观结构信息的一种模 拟方法。
基本单元大小的选择
• 基势本的单截元 断的 距大 离小 ）必或须Rc大ut<于1/22R基cu本t（单R元cut是的相大互小作。用这 保证了任何原子只与原子的一个镜像有相互作用， 不与自己的镜像作用。这个条件称为“minimum image criterion”
• 在我们所研究的体系内的任何结构特性的特征尺 寸或任何重要的效应的特征长度必须小于基本单 元的大小。
• 1991年有人进一步提出了巨正则系综MD方法， 从而又可适用于吸附问题的处理等，该方法还在 进一步发展之中。
• 分子动力学方法的主要发展可见表3.1。
年代 1957 1963 1971
1972 1977 1980
1983 1984 1985 1991
表3.1 MD方法的里程碑工作
创立者
创造内容工作（MD分类名称）
周期性边界条件
• 如果在基本单元内的原子的位置为, 周期边界条 件会产生该原子的镜像。它的位置在
uvimage uv v v v ri ri la mb nc
(3-2-4)
这儿a、b、c是基本单元的三个边长，l，m，n是 整数，取值范围为-∞ 到+∞。
• 在基本单元内的原子不仅与在本单元内的其他 原子有相互作用，还与在相邻单元内的镜像原子 有作用。
③边界条件问题
• 在处理原子、分子的聚合体问题时，就MD方法而 言，能处理的原子(分子)数目要受到计算机运行 速度和能力的限制。目前报导的最好水平是能处 理109量级的原子数目。这与现实物质含有1023个 原子或分子的差距还很大，导致模拟系统原子数 少于真实系统的所谓“尺寸效应”的问题。
• 为了减小“尺寸效应”而又不至于使计算工作量 过大，对于平衡态MD模拟采用 “周期性边界条 件”。
约束系统MD方法
Andersen, Parrinello & Rahman
恒压MD方法
N. J. Gillan & M. Dixon
非平衡MD方法（存在温度梯度）
S. Nosé R. Car & M. Parrinello Cagin & Pettitt
恒温MD方法
第一性原理MD方法 (Car-Parrinello方法)
第三章 经典分子动力学方法
3.1 引言
• 分子动力学 (Molecular Dynamics, 简写为MD) 方法是确 定性模拟方法，这方法是按该体系内部的内禀动力学规律 来计算确定位形的转变。
• 首先需要建立一组分子的运动方程，然后通过直接对系统 中的每一个原子/分子运动方程进行数值求解，得到每个 时刻每个原子/分子的坐标与动量(速度)，即在相空间的运 动轨迹，再利用统计方法得到多体系统的静态和动态特性， 从而得到系统的宏观性质。
②力学条件
• 刚体力学方法是把分子作为刚体处理，建立对于 刚体的欧拉方程和对于质点系的牛顿方程，联立 求解所研究对象的运动问题，以前主要是在处理 像水和四氯化碳那样的低分子量体系，现已用于 研究晶体的相变。
• 约束力学是冻结粒子体系的一部分自由度，进而 求解因此而生成约束条件下的质点力学运动方程。 对于有机分子来讲，因为键长、键角的振动变化 对体系的结构影响不大，将这些自由度冻结是合 适的。另外也有采用固定晶体结构的考虑。
(3-2-1)
式中mi为所考察的原子质量，ri为原子的位置坐标，Fi为 作用在原子上的原子相互作用的合力，它由下式给出
uur
N
Fi i ij
j 1
(3-2-2)
其中，Φij是原子和原子j之间的势函数(有时亦称为力场)
例如，由氩原子等组成的稀薄气体，其势 函数可采用Lennard－Jones势，
④数值积分法
• MD方法的基本方程是线性或非线性的二阶常微 分方程。对此人们研究了许多求解方法。例如， 贝鲁勒（Berreele）法，阿达姆斯（Adams）法， 龙格-库塔（Runge-Kuta）法，追赶法等。
• 最近，M.Tuckerman 和B.J.Berne 提出了所谓多 重时间宽度计算方法。即在质量相差很大的多原 子体系，力学常数有很大不同的体系以及含有长、 短程力的体系等情况下，必须结合具有最小振动 周期的自由度选择积分的时间标度（ ）.这t 方法 缩短了计算用的时间，有望成为MD方法中有前 途的数值积分法。
• 现在已经能够处理许多体系，例如：
– 孤立宏观团簇的模拟(用NEV或NTV系综) – 固体的结构相变，玻璃转变，晶化过程的模拟(用NTP系综) – 固体(晶体)表面的原子、分子吸附现象的模拟 (用μVT系综)
②力学条件
已建立了弹性力学、质点力学、刚体力学、约束 力学等不同力学条件下的四种体系的MD方法。
在Δt时间内，对系统内的所有 粒子解运动方程


ri (tn ) ri (tn1)


vi (tn ) vi (tn1)
tn1 tn t
计算要求的物理量，将数据写 入轨迹文件
No
tn1 tmax
t>tmax
Yes
输出计算结果，并结束计算
MD方法
到目前为止已经确立的MD方法的主要技术体系
①统计系综
• 系综是一个巨大的系统，由组成、性质、尺寸和形状完全 一样的全同体系构成数目极多的系统的集合。不同的系综， MD方法的基本方程有所不同。
• 目前除微正则系综（NEV系综）外，已完成了正则系综 （NTV系综），等温等压系综（NTP系综），等压等焓系 综（NHP系综），巨正则系综（μVL系综），恒定化学势 系综（μVT系综）等五个系综的MD方法的基本方程的确 立。
• 以二维边界条件为例。只在x-y平面配置基本单元的复制 品，使用周期性边界条件，同时在z方向不赋予周期性边 界条件，固定两端的数层原子。由于采用了这样的人工边 界条件，所以在z轴方向上的原子层数要有适当的数目， 一般考虑的标准线度为4－5纳米。采用二维边界条件，就 可以用MD在原子尺度上研究界面的结构。
对于基本单元内的原子、分子运动方程，使用什么样的形
式合适，要具体问题具体分析。若是考虑具有确定的粒子 数N，体积V和能量E的NEV系综(称为微正则系综，MicroCanonical Ensemble)，则其运动方程可以表达成式(3-2-1) 所示的普通牛顿方程的形式
r
mi
d 2 ri dt 2
uur Fi
• 通常 t ～10-15s=1fs
3.3 MD模拟的数据分析
MD模拟给出的直接数据是坐标和速度，对 它们进行加工，可得到各种所要的信息。 这些信息大致可分为3类： • ·所模拟系统的平衡态性质 • ·系统处于亚稳态时的结构与性质 • ·系统在远离平衡态时的动力学过程
平衡态

• 一个物理量A随时间变化而趋于平衡态的值A0的过程可用 下式表示
• 在MD方法的处理过程中，方程组的建立是通过对物理体 系的微观数学描述给出的，在这个微观的物理体系中，每 个原子/分子都各自服从经典的牛顿力学定律。