分子模拟教程

N1 (V - V1 ) Pparticle (m → n) = min{1, exp[−δU nm / kT ] } (N - N1 + 1)V1 change
总结与讨论
系综的等效性
采用什么系综来求体系的热力学量是无关紧要的， 采用什么系综来求体系的热力学量是无关紧要的， 其结果实际上都是一样的，这称为不同系综的热力 其结果实际上都是一样的， 学等同性。 学等同性。
GEMC 的配分函数
对于原子系统，位型（构型） 对于原子系统，位型（构型）的配分函数
V1N1 (V − V1 ) N − N1 Q( N ,V1 ,V2 , T ) = ∑ 3 N N1!( N − N1 )! N1 =0 λ
N N ds2 − N1 ∫ ds1N1 exp[−U (s1N1 / k BT )] ∫
气体分子运动论 其它
主要的学习参考书籍： 主要的学习参考书籍：
1. Computer Simulation of Liquids. --- M. P. Allen and D.J. Tildesley, Oxford University (1987) 2. The Art of Molecular Dynamics Simulation. --- D. C. Rapaport, Cambridge University (1995) 3. Understanding Molecular Simulation: From Algorithms to Applications. ---- D. Frenkel and B. Smit, Academic (1996)
exp(µN / k B T )V N PµVT ( s N , N ) ∝ exp[−U ( s N )] λ3 N N !
Metropolis GCMC algorithm
产生巨正则系综的马尔可夫链的过程涉及到典型的、 产生巨正则系综的马尔可夫链的过程涉及到典型的、 三种不同的随机移动： 三种不同的随机移动：
巨正则系综配分函数
对于原子系统，位型（构型）的配分函数： 对于原子系统，位型（构型）的配分函数：
V N exp( Nµ / kT ) Ξ( µ , V , T ) = ∑ ⋅ ⋅ ⋅ ∫ exp[−U (s)/ kT ]dS 3 N ∫ N!λ3 N N =0
∞
其中， s 为标度坐标，r = V1/3 。 其中， 为标度坐标， 概率密度为： 概率密度为：
但是： 但是：
系综选择关系到热力学量的涨落大小； 系综选择关系到热力学量的涨落大小； 系综选择决定于物理现象的特点和关键物理量的计 算方便与否。 算方便与否。
3. MD方法简介 方法简介
概念： 概念：
从系统中各粒子间的相互作用解每个粒子的 牛顿运动方程: 牛顿运动方程: (F=ma) 。 =ma
3.1 MD 的一般算法： 的一般算法：
Attempt to move a particle (just like canonical MC) Attempt to create a particle Attempt to delete a particle
各种随机移动的概率： 各种随机移动的概率：
Pmove (m → n) = min{1, exp[−δU nm / kT ] }
2. 什么是计算机分子模拟方法？ 什么是计算机分子模拟方法？
分子模拟的定义： 分子模拟的定义：
从统计力学基本原理出发，将一定数量的分子输入 从统计力学基本原理出发， 计算机内进行分子微观结构的测定和宏观性质的计算。 计算机内进行分子微观结构的测定和宏观性质的计算。
按照获得微观态的方法不同，分子模拟分为： 按照获得微观态的方法不同，分子模拟分为：
某个特定构型的发生概率为 PNVT(r)
PNVT (r) ∝ exp[−U (r) / k BT ]
典型算法 (Metropolis)
随机选择一个原子，并且计算其当前构型m 随机选择一个原子，并且计算其当前构型m的能量 U(rm). 尝试随机移动此原子，rn = rm + ∆×(ranf-0.5), 并且 ∆×(ranf尝试随机移动此原子， 计算此新构型n 计算此新构型n时的能量 U(rn). 接受此移动的概率为 Pmove
3.4 经典分子动力学的分类： 经典分子动力学的分类：
NVE分子动力学 NVE分子动力学 恒温分子动力学 恒压分子动力学 恒温恒压分子动力学
其它分子动力学
粗粒化分子动力学（Coarse-grained molecular 粗粒化分子动力学（CoarseDynamics, CGMD) 耗散粒子动力学 (Dissipative Particle Dynamics，DPD) Dynamics，
Pmove = min{1, exp[−δU nm / kT ] }
其中，δUmn = U(rn)- U(rm) 其中，
2.2 巨正则 MC simulation (GCMC)
恒定 V, T, 和 µ，体系的粒子数发生波动； 体系的粒子数发生波动； 可用于预测 EOS-type的性质，但主要是用来模 EOS-type的性质， 的性质 拟吸附过程. 拟吸附过程.
N=N1+N2 V=V1+V2 constant T
V-V1 V1
GEMC 模拟算法： 模拟算法： 随机选择一个粒子进行移动 (NVT). 改变每个模拟盒子的体积，但总体 改变每个模拟盒子的体积， 积保持不变 (NTP). 盒子间交换粒子 (µVT)。 VT)。
各种随机移动的概率： 各种随机移动的概率：
计算机分子模拟 (Molecular Simulation)
计算机分子模拟的部分应用领域： 计算机分子模拟的部分应用领域：
• 生物、制药： 生物、制药： • 力学、物理学： 力学、物理学： • 化学、化工: 化学、化工: • 微电子、微机械： 微电子、微机械： 大分子性质、药物设计 大分子性质、 应力与裂纹扩展、团簇研究 应力与裂纹扩展、 溶液理论、吸附、界面化学 溶液理论、吸附、 半导体工艺、微加工、 半导体工艺、微加工、 超薄膜润滑： 超薄膜润滑： • 地质、矿产： 地质、矿产： 地核动力学、熔融盐结构 地核动力学、 及其性质
从上个世纪九十年代初期以来，计算机模拟技 从上个世纪九十年代初期以来， 术得到了飞速Biblioteka Baidu展，主要基于三个方面的发展： 术得到了飞速发展，主要基于三个方面的发展：
分子力场的发展（基石） 分子力场的发展（基石） （Amber，OPLS、Compass） ， 、 ） 原子间的键长、键角、 原子间的键长、键角、分子间的内聚能等
Pmove (m → n) = min{1, exp[−δU nm / kT ] }
(V1n ) N1 (V - V1n ) N - N1 Pvolume (m → n) = min{1, exp[−δU nm / kT ] } m N1 m N - N1 (V1 ) (V - V1 ) change
分类： 分类：
正则系综 (Canonical MC) (NVT) MC， 巨正则系综 (Grand Canonical MC，GCMC) (µVT) Gibbs系综 Gibbs系综 (Gibbs Ensemble MC， GEMC) MC，
2.1 正则系综 正则系综MC (Canonical MC Simulation
< A(r, p) >τ =< A(r, p) > z lim
τ →∞
3.3 解经典牛顿运动方程的算法： 解经典牛顿运动方程的算法：
数值积分 有限差分法 泰勒展开 Verlet算法： Verlet算法： 算法 优点：精确，时间可逆； 优点：精确，时间可逆； 缺点：速度有较大误差，轨迹与速度无关， 缺点：速度有较大误差，轨迹与速度无关，无法与热 浴耦联 Leapfrog形式 Leapfrog形式 优点：轨迹与速度有关， 优点：轨迹与速度有关，可与热浴耦联 缺点：速度近似； Verlet算子多花时间 缺点：速度近似；比Verlet算子多花时间 预测－校正格式 预测－
1) 在某t时刻，体系内的粒子i 在某t时刻，体系内的粒子i具有位置 ri(t), 速度 vi(t). 从ri(t), 计算力fi(t). 利用质量mi, 获得ai(t). 计算力f 利用质量m 获得a 利用r 利用ri(t), vi(t), ai(t) 来计算出下一时刻 (t+dt) 时 的位置r 速度v 的位置ri(t+dt), 速度vi(t+dt). (t+dt)，从第一步开始循环重复， 利用r 利用ri(t+dt), vi(t+dt)，从第一步开始循环重复， 直到我们所需要的模拟时间长度为止。 直到我们所需要的模拟时间长度为止。
Nλ = min{1, exp[−( µ + δU nm ) / kT ] } V
3
Pdelete
Pcreate
V = min{1, 3 exp[( µ − δU nm ) / kT ] } λ ( N + 1)
2.3 Gibbs Ensemble MC (GEMC)
特点: 特点:
尤其适用于研究纯流体或混合物的相平衡问题； 尤其适用于研究纯流体或混合物的相平衡问题； 此方法不能用于涉及到非常稠密流体的相平衡问题； 此方法不能用于涉及到非常稠密流体的相平衡问题； 此方法能同时获得共存相的各自密度及其组成； 此方法能同时获得共存相的各自密度及其组成； 此方法避免了共存相界面的问题。 此方法避免了共存相界面的问题。
2) 3)
4)
< A(r , p ) >=
1
τ
∫ 0 A(r (t ), p (t ))dt
τ
3.2 各态遍历假说： 各态遍历假说：
一个力学体系在长时间的运动中， 一个力学体系在长时间的运动中 ， 它的代 表点可以无限接近能量曲面上的任何点， 表点可以无限接近能量曲面上的任何点，系 综平均等于长时间的时间平均。 综平均等于长时间的时间平均。
(1) 蒙特卡罗方法 (Monte Carlo, MC) (2) 分子动力学方法 (Molecular Dynamics, MD) (3) 混合方法 (hybrid method，HM) ，
计算机分子模拟的发展历史： 计算机分子模拟的发展历史：
1. 蒙特卡罗方法（MC） 蒙特卡罗方法（MC）
1953 Metropolis, Ulam, Rosenbluth and Tell Los Alamos National Lab Monte Carlo simulation of hard sphere.
掌握分子模拟方法的必备知识： 掌握分子模拟方法的必备知识：
编程技能 (Fortran or C/C++) 统计物理学（统计力学） 统计物理学（统计力学）：
统计物理学基础； 统计物理学基础； 系综原理； 系综原理； 非平衡统计力学基础； 涨落理论 非平衡统计力学基础；
分子热力学 ：
分子间相互作用理论； 分布函数理论 分子间相互作用理论；
particles N Volume V Temperature T Initial state r{N}
正则系综的配分函数
对于原子系统，位型（构型）的配分函数： 对于原子系统，位型（构型）的配分函数：
1 Q( N ,V , T ) = ⋅ ⋅ ⋅ ∫ exp[−U (r ) / k B T ]dr 3 N N! N !λ3 N ∫
2. 分子动力学方法（MD） 分子动力学方法（MD）
1957 Alder and Wainwrigth Livermore Lab Molecular dynamics simulation of hard spheres.
计算机分子模拟的发展历史（ 计算机分子模拟的发展历史（续） 分子模拟的发展历史
模拟算法（途径） 模拟算法（途径） 计算机硬件（工具） 计算机硬件（工具）
HPCx
计算机分子模拟的特点： 计算机分子模拟的特点：
1. 原子水平的模拟 2. 计算机实验 3. 检验理论、筛选实验 检验理论、 4. 科学研究中的第三种方法
2. MC方法简介 方法简介
利用马尔可夫链方法产生感兴趣的系综微观态。 利用马尔可夫链方法产生感兴趣的系综微观态 。 对这些微观态进行系综平均来获得热力学性质。 对这些微观态进行系综平均来获得热力学性质 。 只对系统中粒子的位置抽样， 不包含动能部分， 只对系统中粒子的位置抽样 ， 不包含动能部分 ， 提供的是系统的超额性质。 提供的是系统的超额性质。