作用的模型，故假定一个孤立的多粒子体系，粒子间的相 互作用位势是球对称的，则系统的哈密顿量为：
H 1
2
i
1 m
pi2
i j
u(rij )
（rij第i个粒子和第j个粒子 之间距离）
时间不显示地出现在哈密顿量中，系统受一下约束： ①系统能量E守恒； ②系统的粒子数恒定； ③系统的体积不变； ④整个系统并未运动，系统的总动量P恒等于零。
2
分子动力学方法的分类
• 按所遵循的运动规律或研究的对象分类：
经典MD
平衡态MD
非平衡态 MD
微正则系综 （NVE）
正则系综 （NVT）
等温等压系 综（NPT）
量子MD
等焓等压系 综（NPH）
3
分子动力学模拟的基本步骤
i • 设定模拟的物理模型，建立运动微分方程 ii • 给定初始条件，数值求解运动微分方程 iii • 趋于平衡的计算过程 iv • 宏观物理量的计算
16
平衡态分子动力学模拟
经典MD
平衡态MD
非平衡态 MD
量子MD
微正则系综 （NVE）
正则系综 （NVT）
等温等压系 综（NPT）
等焓等压系 综（NPH）
• 微正则系综（NVE）
系统的能量、粒子数和体积均被处理成恒定量。
17
微正则系综
• 前提条件
对微正则系综进行MD模拟时首先要确定系统中粒子相互
18
微正则系综
• 从牛顿运动方程出发：
d 2ri (t) dt 2
1 m
i j
Fi (rij
)
i=1,2,…,N
采用Verlet方法（中心差分格式）求上式：
ri (t h) 2ri (t) ri (t h) Fi (t)h2 / m i=1,2,…,N
i j
可知，由t和t-h的位置及t时刻的受力能求出t+h时刻粒子的位置；
由n+1步位置算出n步的速度，可见动能的计算比势能落后一步。
19
微正则系综
• 具体模拟步骤：
①给定初始空间位置：
r (0)
i
，ri(1)
②在n步时计算粒子所受的力： Fi(n) Fi (tn )
③计算粒子第n+1步的位置
r (n1) i
r (n1) i
2ri ( n )
ri(n1)
Fi ( n ) h 2
式中rij为位于ri处的粒子i同ri处的粒子j之间的距离。
8
具体步骤
➢ 对于短程力采用截断半径法
相互作用势的确定十分关键，然而势函数的确定 往往相当困难，并不能精确知道。
9
具体步骤
• 给定初始条件
一般来说MD中初始条件是不知道的，而且初始 条件往往不重要； 常用的初始条件选择：
1、令初始位置处在离散网格的格点上，初始速度则从波 尔兹曼分布随机抽样确定； 2、令初始位置随机地偏离网格的格点，初始速度取为零； 3、初始位置与速度都随机确定：初始位置随机偏离网格 格点，初始速度有波尔兹曼分布随机抽样确定。
分子动力学方法
1
概述
• 首先，分子动力学方法是一种确定性方法。 • 分子动力学方法是按照该体系内部的动力
学规律来确定位置和速度随时间的演化的 一种模拟方法。 • 跟踪系统中每个粒子的个体运动。 • 根据统计物理规律，给出微观量（座标、 速度）与宏观可观测量（温度、压力、比 热、弹性模量等）的关系，从而研究物质 和材料的性能。
/m
④计算第n步的速度，Vi ( n )
(ri(n1)
r (n1) i
)
/
2h
⑤返回步骤2，开始下一次模拟计算。
• 改进：
把N个粒子的初始位置放置在网格的格点上，然后加以扰 动，给出的初始条件是粒子的空间位置和运动速度，可用
如下公式计算粒子位置： ri(1) 2ri(0) hvi0 Fi(0)h2 / 2m 20
10
具体步骤
• 数值求解运动方程和趋衡计算
模拟计算首先从数值求解运动方程出发，方法有： ➢ 一步法：
采用有限差分法，将微分方程变为有限差分方程， 取差分计算的时间步长为h,运用到展开的一阶泰 勒公式：
11
具体步骤
整理后得到微分方程的欧拉算法：
这是一组递推公式。 ➢ Verlet方法（二步法的一种）
微正则系综
• 上述Verlet算法可以作变形改进，使数值计算的稳定性得 到改善，即为Verlet算法的速度形式：
令，
z(n) i
(ri(n1)
ri(n) ) /
将泰勒公式二次展开
12
具体步骤
➢ 蛙跳Leap-Frog算法：
13
具体步骤
14
具体步骤
➢ 另外还有速度Verlet算法 可以同时给出位置、速度与加速度，并且不牺牲 精度。 预测—校正算法
➢ 时间步长h的选取 过大误差太大，过小平衡时间太长，同时增加累 加误差。
15
具体步骤
• 宏观物理量的计算 最终目的 通过相对轨迹求平均得到。
值， σ可以确定长度的单位。
➢ 周期性边界条件
A(x)=A(x+nL) A为任意的可观测量
7
具体步骤
➢ MD中考虑粒子间相互作用时，通常采用最小像力 的约定： 在无穷重复的分子动力学基本元胞中，只考 虑每个粒子只同它所在的基本元胞内的另外N-1个 粒子或其最邻近的影像粒子相互作用；
既有， rij min(| ri rj nL |)
现令， tn nh
r (n) i
ri (tn )
F (n) i
Fi (tn )
上式可以改写为：ri(n1)
2ri(n)
r (n1)
i
Fi(n)h2 / m , i=1,2,…,N
i j
即为计算粒子坐标随时间演化的递推公式；
一组初始空间位置
随后粒子的所有位置；
粒子运动速度：Vi(n) (ri(n1) ri(n1) ) / 2h
4
具体步骤
• 设定模拟物理模型
首要步骤
➢ 对经典MD而言，每个分子遵循牛顿运动规律
mai(t)=Fi(r), i=1,2,…,N 总Fi 和是。体系内N个经典粒子作用在第i个粒子上力的
N
^
Fi (r) Fi (| ri rj |) rij j 1
i j
5
具体步骤
➢分子之间相互作用势：
势函数的选取是模型设定步骤中的关键 M最先被应用的势函数）
∞， 如果r<σ V(r)=
0， 如果r≥ σ 其中σ为相互作用势为0时的两分子之间的距离。
6
具体步骤
2、Lennard-Jones位势（最简单和最常用）
V (r) 4[( )12 ( 6 )]
r
r
其中，-ε是位势的最小值（即势阱最低处的V(r)
值）， ε可确定能量的单位；σ为V(r)等于零时的r