分子动力学简介

1.4 分子动力学模拟系综(ensemble)

平衡分子动力学模拟，总是在一定的系综下进行的。
微正则系综(micro系统原子数N,体积V,能量E保持不变。是 canonical ensemble, 孤立、保守的系统。和外界没有能量或物 质交换，只有动能和势能之间的转换。 NVE) 研究扩散系数
正则系综(canonical ensemble, NVT)
每个原子的速度按此式随机产生
⑥计算第n+1步的位置
⑦计算第n步的速度
Water Model Geometries
Explicit Water models SPC, SPC/E, TIPnP (n=3,4,5,6), POL5
Water Model Parameters
• SPC, SPC/E (Berendsen) • TIP3P, TIP4P, TIP5P (Jorgensen) • TIP4P/FQ, POL5 (Berne)
初始化
能量优化
升温
长时间平衡模拟
数据分析
根据所有分子的当前坐标计算分子的受力 根据受力更新分子的坐标 在此过程中收据用于计算宏观性质的有关信息
避免局部分子重叠
平衡分子动力学模拟的基本步骤
1.6 数据分析：微观到宏观信息转化

分子动力学模拟可以得到原子显微水平的信息，包括原子位 置和运动速度。要将这些微观信息转化为宏观水平信息，如 压强、能量、热容量等，需要用到统计力学理论。统计力学
是分子动力学模拟生物系统研究的基本理论。
温度: T
1 dNk B
mi vi2
i
Ei / k BT
1 ' ' 2 溶液扩散系数: D lim r (t t ) r (t ) 6t
) N ,T
压强:
ln( i e p k BT V
热容:
2 ( E E ) E CV (k BT ) 2 T V
常用的分子力场
AMBER —Assisted model building with energy refinement CHARMM—Chemistry at Harvard macromolecular mechanics GROMOS—Gronigen molecular simulation
Gsol Gel Gnonel
Gel Poisson Bolzmann equation
Gnonel surface area of solute
静电部分溶解能 (GB模型）
Coulomb field integral:
2.3 增强取样方法


以空间坐标为导向
1.2 势能函数形式
d2 mi 2 ri V ({ri }) dt ri
简单力场的函数形式
第一项表示成键原子间的相互作用 第二项表示分子中所有键角的贡献，用的均为谐振子模型
第三项为扭曲势，表示键旋转时的能量变化 第四项是非键项，在简单力场中，用 Lennard-Jones 势描 述van der Waals作用,库仑势描述静电相互作用。
张军华等，石油地球物理勘探，2010， 45， 918
2.2 显式vs隐式水模型
A: 显式模型：TIP3P，TIP4P，SPC/E… B: 隐式模型：GB/SA (general Born/surface area)
GB/SA隐式模型

通过计算溶质的溶解化能来模拟溶液环境对蛋白质的影响。 将溶解化能分成静电和非静电两个部分。
蛙跳（Leap-frog)算法：半步算法
1 1 ri (t t ) ri (t ) vi (t )t ai (t )t 2 bi (t )t 3 2 6 1 1 ri (t t ) ri (t ) vi (t )t ai (t )t 2 bi (t )t 3 2 6
力场(force field)
• 描述势能的方程以及其中的参数构成了力场。 • 势能函数和结构参数决定了分子力场的性能即它的计 算结果的准确性和可靠性。 • 参数的决定要在大量的热力学、光谱学实验数据的基 础上进行,有时也需要由量子化学计算的结果提供数 据。 • 各类键长、键角的“本征值”一般取自晶体学、电子 衍射或其他的谱学数据，键伸缩和角变力常数主要由 振动光谱数据确定，扭转力常数经常要从分子内旋转 垒来推算。
1 1 vi (t t ) vi (t t ) 2 2 vi (t ) 2
1 2

t-△t/2
t
t+△t/2
t+△t
T+3△t /2
t+2△t
Leap-frog算ຫໍສະໝຸດ Baidu的表述:
算法启动 ①规定初始位置 ②规定初始速度 ③扰动初始速度： vi (t / 2) vi (0) ai (0)t / 2 ④计算第n步的力
伞状取样（umbrella sampling) 引导分子动力学（steered molecular dynamics) 构象填充（conformation flooding) 绝热分子动力学（adiabatic molecular dynamics) Metadynamics
1 算法启动 r (t t ) r (t ) v (t )t a (t )t 2 ①规定初始位置 ②规定初始速度 ③扰动初始位置： ri (t ) ri (0) vi (0)t ④计算第n步的力

i i i i 2
1 bi (t )t 3 6
2 r ( t t ) 2 r ( t ) r ( t t ) a ( t ) t ⑤计算第n+1步的位置：i i i i ri (t t ) ri (t t ) ⑥计算第n步的速度： vi (t ) 2t
周期性边界条件（Periodic boundary conditions）



周期性边界条件使人们能够基 于少量粒子的模拟的同时考虑 周围溶液和溶质分子对研究体 系的影响。 把研究体系看作一个特定形状 空间包围的区域，采用周期性 边界条件，则这个基本单元会 沿着所有的方向进行周期性扩 展以形成一个无限的周期排列。 当基本单元中的粒子离开这个 单元进入一侧的映射单元时， 其映射粒子会从基本单元的另 一侧进入基本单元。

两式微分后相减
vi (t t ) vi (t t ) 2ai (t )t
首先利用当前时刻的加速度，计算半个时间步长后的速度：
1 1 vi (t t ) vi (t t ) ai (t )t 2 2

计算下一步长时刻的位置： ri (t t ) ri (t ) vi (t t )t 计算当前时刻的速度：
粒子位置： ri (t t ) 2ri (t ) ri (t t ) ai (t )t 2 粒子速度：
ri (t t ) ri (t t ) vi (t ) 2t
由前两个时刻的位置,计算获得下一时刻的位置、这一时刻的速度 需要连续记录两个时刻的位置
Verlet算法的表述:
势能面的复杂性导致分子模拟构象取样 的困难
如何加快分子模拟的计算效率是目前我们面 临的重大挑战



计算机运行速度的快速发展 隐式水模型的应用 增强取样方法的发展
2.1计算资源的快速发展
Anton computer, developed by D.E. Shaw Research
Dror et al., Annu. Rev. Biophys. 2012, 41, 429.
方程①中的能量Eel(势能面)仅仅是原子核坐标有关。相应的，方程② 所表示的为在核势能面E(R)上的核运动方程。 直接求解方程①，采用的是从头算或者是半经验，这样的量化计算都 是把电子的波函数和能量处理成原子核坐标的函数。由于量子化学求 解电子波函数和势能面耗时巨大，常常将势能面进行经验性的拟合， 成为力场，由此构成分子力学的基础。
分子动力学依靠牛顿力学来模拟分子体系的运动，以在由分子 体系的不同状态构成的系统中抽取样本，从而计算体系的热力 学量和其他宏观性质。

分子力学中势能面是用解析经验势函数来描述的，体系 的势函数可以分解成多种势函数的共同作用。通常分为键 伸缩能Es 、键弯曲能Eb 、二面角扭曲能Etor、范德华作用 能Evdw 和静电作用能Eelec

⑤计算第n+1/2步的速度：vi (t 1 / 2t ) vi (t 1 / 2t ) ai (t )t ⑥计算第n+1步的位置： ri (t t ) ri (t ) vi (t 1 / 2t )t ⑦计算第n步的速度：
vi (t ) (vi (t 1 / 2t ) vi (t 1 / 2t )) / 2


1. 分子动力学简介
1.1基本原理


一个体系有N个原子 体系的状态由这N个原子的位置{ri}和动量{pi}或速度{vi}来 标志。 体系的能量为：H {ri, pi } 体系的运动方程为：

pi H t ri
ri H t pi
分子动力学的主要目的: 解上面的方程求得体系状态相空间演化的轨迹{ri, pi}t0, {ri, pi}t1, {ri, pi}t2, {ri, pi}t3,……

Verlet算法
1.3 分子动力学的算法
粒子位置的Taylor展开式：
1 1 2 ri (t t ) ri (t ) vi (t )t ai (t )t bi (t )t 3 2 6 1 1 2 ri (t t ) ri (t ) vi (t )t ai (t )t bi (t )t 3 2 6


将方程②用牛顿运动方程代替，势能面采用力场拟合，就构成了分子 动力学的基础。
分子级别的模拟


分子水平的模拟

以分子的运动为主要模拟对象 一般情况下不考虑电子转移效应，因而不能准确模拟化 学成键作用 采用经验性的分子间作用函数模拟微粒之间的作用 分子模拟的最早发展 1950s, Alder, 劳伦斯利物默实验室，分子动力学模拟32 个原子 1950s, Metropolis, 洛斯阿洛莫斯实验室，蒙特卡洛模 拟32个原子
分子动力学模拟：基本原理及 应用介绍
邵强 药物发现与设计中心 中科院上海药物所 2015.09.09
k e 0
电子运动方程： 核运动方程：

k e 0
① ②
ˆ V ) ( R, r ) E ( R, r ) (H el NN e el e ˆ E ) ( R) E ( R) (T el N N

其他的算法还包括： Velocity Verlet算法 预测-校正（predictor-correlation)格式算法
算法表述:
算法启动 ①规定初始位置 ②规定初始速度 ③扰动初始速度 ④计算第n步的力 ⑤计算第n+1/2步的速度

1 3 2 能量均分原理: 2mi vi 2 Nk BT
研究溶液特性 等温等压系综(NPT) 研究相变
系统原子数N,体积V,温度T保持不变,且总 动量保持不变。虚拟热浴耦合→系统动能 固定→保持温度不变→原子速度标度
系统原子数N,压强P,温度T保持不变 。压 强P与体积V共轭，控压可以通过标度系 统的体积来实现。
1.5 分子动力学模拟的基本步骤
读入模型参数、模拟控制参数
自由能:
A k BT ln( P)
分子模拟与统计力学结合，为实验宏观 观测量提供分子机理。 对研究体系的结构、动态运动做出预测, 给实验提供理论指导。
24
2. 分子力学模拟应用局限性及解决方案
例：用100ns/天的计算速度研究半衰期为1毫秒的反应， 通过计算观察到1次反应平均需要超过20年！
在实际的应用中，我们把哈密顿方程化为下面的牛顿方程， 并且用位置ri 和速度vi 做为描述体系的参量。
1 N H(ri ,pi ) mi vi2 V ({ri }) 2 i 1
pi H t ri
d2 mi 2 ri V ({ri }) dt ri
Fi mi ai