分子动力学模拟

系综调节
系综调节主要是指在进行分子动力学计算 过程中，对温度和压力参数的调节
调温技术： Berendsen热浴、速度标度、Gaussian热浴、
Nose-Hoover热浴
调压技术：Berendsen方法、Anderson方法、Parrinello-Rahman方法
三、MD模拟所需条件
MD模拟所需条件 初始条件：模拟对象的起始位置，速度，
任务：通过求解经典牛顿运动方程，计算一个经典多 体体系的平衡和非平衡性质 系统描述：粒子坐标x，速度（动量）v，受力f，时间t 模拟体系大小：几百到上百万个粒子，对应于几个到几十个nm。
MD模拟的一般过程
构建构型 动力学过程模拟 构型性能计算
MS构建 MS 晶胞
所需条件
势函数 系综 初始条件 周期性边 界条件
c p c p c p
c
p
谢谢
2011.12
i≠ j
ρi (rij ) 为第j个原子在i个原子处贡献的电荷密度
长程F-S势函数
对势
c:正的无量纲常数 ε：有能量量纲的参数 α：有长度量纲的参数
多体势 m, n：正整数 对于不同研究体系，5 个参数取值不同
系综简介
系综（Ensemble） 系综（Ensemble）：相空间中具有相同热力学性质的所有
结果分析
RDF,CN等
二、势函数与系综
势函数简介
对势（Pair potential ）：认为原 子间的相互作用是两两之间的作 用 与其他原子的位置无关
硬球势、Lennard-Jones势、 Morse势、Born-Lande势 及Johnson势
原子间作用势 多体势（Many-body effects）： 在多原子体系中 一个原子的位 置不同 将影响其它原子间的有 效相互作用 嵌入原子法（EAM 势）、多体相互作用 势（FS势）、TB势 等
1. Verlet算法: 将粒子位置以泰勒式展开
d 1 d r (t + δ t ) = r (t ) + r ( t )δ t + dt 2 ! dt
2 2
r ( t )( δ t ) 2 + ....
d 1 d 2 r (t − δ t ) = r (t ) − r ( t )δ t + r ( t ) (δ t ) 2 + .... dt 2! dt2
0 0
0
牛顿运动方程
r i vi
d2 d ri = vi = ai 2 dt dt
0
0
vi = v i + a i t
1 2 r i = r i + vi t + a i t 2
0 0
0
给 定
反 复 计 算 得 轨 迹
计 算
t=δt
δt
femto second
r i vi
求解
常用算法： Verlet算法，Leap-frog(蛙跳 算法 Gear 算法 蛙跳)算法 蛙跳 算法,
关于分子动力学模拟的 初步认识
摸金校尉
主要内容
分子动力学模拟概述 MD模拟所需条件 势函数与系综 牛顿运动方程及其求解
一、分子动力学模拟概述
为什么要搞MD模拟
Chemistry is no longer a purely experimental science.
无法获得过程中粒子 微观细节，成本高等
∆ a (t + δ t ) = ∆ a (t + δ t ) − ∆ a (t + δ t )
p
c
引入常数c0, c1, c2, c3修正各变量得
r (t + δ t ) = r (t + δ t ) + c0 ∆ a(t + δ t ) v (t + δ t ) = v (t + δ t ) + c1∆ a(t + δ t ) a (t + δ t ) = a (t + δ t ) + c2 ∆ a(t + δ t ) b (t + δ t ) = b (t + δ t ) + c3 ∆ a (t + δ t )
执行温度，积分步长等值得确定。直接关系 到模拟计算的复杂程度。 初始位置可采用能量最小化的方法取能量最 低的结构为起点。均匀相的液态系统常取其 晶体结构；不知道结构的可以以面心立方为 起点 初始速度由初始温度下的Maxwell-Boltzmann distribution随机选取
MD模拟所需条件 周期性边界条件
实验 方法
指导
计算机 模拟
定义
分子动力学模拟（ Simulation） 分子动力学模拟（Molecular Dynamics Simulation）：通过
计算机对原子核和电子所构成的多体体系中的微观粒子之间相互作用 和运动进行模拟，把每一原子核视为在全部其他的原子核和电子所构 成的经验势场的作用下按照牛顿定律进行运动，进而得到体系中粒子 的运动轨迹，再按照统计物理的方法计算得出物质的结构和性质等宏 观性能。
Lennard-Jones势（LJ)
间距为R的两个原子总势能：
排斥项
吸引项
σ和ε为因原子而异的势能参数 势能最低点为r=21/6σ, σ大小表征原子 间平衡距离。ε为由势能最低点到势 能为0点的差。 L-J势能曲线
EAM势（嵌入原子法）
系统中能量：
嵌入能 对势项
ρ i 是除第i
个原子以外的所有其它原子的核外电子 在第 i 个 原子处产生的电子云密度之和：ρi = ∑ ρi (rij )
1 vi = 2 1 1 v t+ 2 δ t + v t- 2 δ t
3. Gear 算法 （校正预测法-predictor –corrector method)）
将位置函数泰勒展开
1 1 2 r (t + δ t ) = r (t ) + v(t )δ t + a (t )δ t + b(t )δ t 3 + ... 2 6
r ( t+ δ t ) = 2 r ( t ) -r ( t- δ t ) + δ t a ( t )
2
v ( t ) = r ( t+ δ t ) -r ( t- δ t ) / 2 δ t
2. Leap-frog(蛙跳 算法 蛙跳)算法 蛙跳
δt δt v i t+ = v i t + a i ( t )δ t 2 2 δt ri ( t+ δ t ) = ri ( t ) + v i ( t + )δ t 2
点的集合。
微正则系综 (micro canonical ensemble)--NVE系综，al ensemble)--NVT系综，动量为0，封闭体系 巨正则系综 (grand canonical ensemble)-- µVT系综，开放体系 吉布斯系综 (Gibbs ensemble)-- NPT系综 等压等焓系综（constant-pressure, constant-enthalpy ensemble ） --NPH系综
p
1 v (t + δ t ) = v(t ) + a (t )δ t + b(t )δ t 2 + ... 2
p
a (t + δ t ) = a (t ) + b(t )δ t + ... b (t + δ t ) = b(t ) + ...
p
p
v, a, b,为r的1次，2 次，3次微分
由于速度加速度均来自泰勒展开式，而非牛顿运动方程， 故并不完全准确。正确加速度与预测加速度差值为：
总作 用力：
牛顿运动方程
原子i受力：
加速度：
∂ ∂ ∂ Fi = −∇Ei = −i − j − k Ei ∂y ∂z ∂x
F i ai = mi
d2 d ri = vi = ai 2 dt dt
i原子经过t时间后的位置
vi = v i + a i t
1 2 r i = r i + vi t + a i t 2
（periodic boundary condition）
：是为了
解决少数粒子来模拟宏观体系的问题而引入的。模拟体系由 基本单元在各个方向上重复叠合而成，模拟时只需保留基本 单元，其他单元与基本单元由平移对称性关联。
rc<L/2
四、牛顿运动方程及其求解
分子力场（Force Field）
原子 i 在其它原 子的作用势场 Ei (ri) 中运动