第六章 分子动力学方法
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
11 10 所以模拟 1000步,弛豫时间只有 s. 3.4宏观量的Biblioteka Baidu量
1014 s,
在分子动力学模拟中 ,宏观物理量的测量是通过对 t 1 N N A lim A ( r ( t ), p (t ))dt 相空间轨迹的平均来得到的: t t t 0 t 例如,动能为 n N 1 1
2
dp H kx dt x
用数值积分计算在相空间中的轨道 ( x(t ), p(t )) 方法是将微分变成差分，引入时间步子 h 作泰勒展 开 df
dt
[ f (t h) f (t )] / h
dx [ x(t h) x(t )] p(t ) 于是，哈密顿方程变为 dt h m dp [ p(t h) p(t )] kx(t ) dt h
3.3趋于平衡 在分子动力学模拟计算中 , 时间步子 $h$ 的选择是 很重要的 . 为了减小误差 ,h 必须取得小些 ; 但如果 h 取 得太小,趋于平衡的时间就要长,因此,为了更快地趋于 平衡 ,h 又应该取得大些 , 因而 , 在相空间中走的步子大 一些.对几百个Ar分子,采用L--J位势,发现取 h 102 就可得到好的相图结果.h是无量纲的,实际时间为
用中心差分法,有 所以
d 2 ri (t ) 2 [ r ( t h ) 2 r ( t ) r ( t h )] / h i i i dt 2
ri (t h) 2ri (t ) ri (t h) Fi (t )h 2 / m
t n nh, ri n ri (t n ), Fi n Fi (t n )
分子动力学方法 （ＭＤ）
1.引言
什么是分子动力学方法(MD)? 确定论的模拟方法实际上在本书第一章介绍的的第 一个计算机模拟工作 , 即费米 -- 巴斯塔 -- 乌拉姆实 验中已经产生了 . 但是在统计物理中所用的分子动 力学方法则产生于 1957 年 . B.J.Alder 等首先用此 方法研究了硬球气体的状态方程.
12 6 例:单原子(Ar)粒子系统的MD计算. u (rij ) 4 r 取256个粒子,位势取为L-J势: r ij ij
为位势的极小值(取 为能量单位),其位置在
r 2(1/ 6)
dx 1000( x t ) 2 2t 例:用 Euler 法解如下初值问题 dt x(0) 0
计算到x(1) , 取 h 10m , m 0,1,2,3, , 结果如下:
可见要取m=3,4,… ,才能得到稳定的结果.在实际数值 计算中,必须十分注意舍入误差与稳定性问题.
2
两式相加相减分别得到
F (t ) x(t h) x(t h) 2 x(t ) h m x(t h) x(t h) v(t ) 2h
于是，通过递推求得数值解，称为Verlet方法。两步法 还有龙格－库塔方法和辛算法（冯康先生的贡献）
考虑二级项后，谐振子的数值解在取h=0.05 时, 迭 代1,000步后能量变为
T
1
Ek
4.微正则系综的分子动力学模拟
pi 1 u (rij ) 一个多粒子系统的哈密顿量为 H 2 i m i j
2
系统受到的限制为
E const. V const.
N const. P0
牛顿方程
Fi (rij ) d 2 r (t )i 2 dt m i j
即
hp(t ) x(t h) x(t ) m p(t h) p(t ) hkx(t )
这是对 x(t),p(t) 的一组递推方程，给定了初 始条件x(0),p(0) ，就可由 t 0 h 2h 3h 一步步计算出粒子在相空间中的轨迹，这个解 法称为 Euler 方法. 由于 h 是有限量，计算总有误差，正确选择 h 是 重要的 . 图左图右分别给出了取 h=0.05 和 h=0.005 的轨迹.可见当h 取得太大时，偏离太大，以至能量 都不守恒了 :E_0=1, 在 1000 步后变为 E=12.14 ！右 图所取 h 小一个数量级，在同样的时间后能量变为 E=1.28 ，轨迹是闭合的，但结果仍不够好 , 必须改 进计算方法，进一步考虑二级项.
取周期性边界条件 , 并采用最小像力约定 . 采用自然单位，长度与 2 1/ 2 时间单位分别为 (m / 48 ) 与 .方程为无量纲形式 ,对Ar原 子,时间单位为 3 1012 s * * ( T , ) (2.53,0.636) 和(0.722,0.83134) 研究相图上的两个点: 相应的MD模拟盒子分别为L=7.38 和L=6.75. 作为初始条件 , 取原子处于面心立方的格子上 , 而速度按麦克斯韦 分布.位势的截断取两个不同的值, rc 2.5 和 rc 3.6 ,以便观察其影响.在平衡化过程中,调节速度的 1/ 2 * T ( N 1) scaling 因子 2
16i vi
重复这个步骤 , 直到达到所希望的能量 ( 或等价的平均温 度).在此MD模拟中,平衡化用了1000个MD步子, 在图中给出了 动能、位能和总能的演化过程.平均速度应当为 v 1.13 T * / 24 0.3668 而模拟的结果为 v 0.3654 ,考虑到只有256 个粒子,这个符合是相当好的.
0
Ek
n n0
2
n0 i
m(vi )
2
d 由动能可根据能量均分定理得出温度 2 Nk B 压强可通过计算在dA面积元上的 F S S 净动量转移 dFn p 的平均值得到: , 或从维里定理得到 : dA 1 N pV Nk BT ri Fi d i
MD与MC的比较 分子动力学方法实现了 Boltzmann 的统计力 学途径 ,而蒙特卡罗方法实现了 Gibbs 的统计力学 途径. 分子动力学方法可处理与时间有关的过程，因 而可处理非平衡态，是更有力的模拟方法，但计算 量大、占内存多、程序也较复杂; 蒙特卡罗方法则具有程序简单、占内存少的优 点,但难以处理非平衡态.
在分子动力学模拟中 , 在考虑粒子间相互作用时 , 通常采用最 小像力约定,如图所示. 在无穷重复的基本单元中,考虑一个 粒子和它最近的粒子或像粒子之间的作用,即
rij min( ri rj nL)
选取L,使得 rij L / 2 , 即大于L/2的 相互作用可忽略. 3.2.初条件： 对位置和速度作初始化, 原则上初始条件并不重要, 只要时间足够长就会忘掉初始条件的 , 但在实际模拟中 , 初始 条件选得好可以更快趋于平衡.可选择：位置在正规格子上, 速度随机；位置随机地偏离格子 ,速度为0；位置速度都随机.
换一个写法:
因此,可写成 ri n1 2ri n ri n1 Fi n h 2 / m 在以上式子中,可以由 ri0 , ri1 得到 ri 2 , ri3 , ,而粒子的速度为
v (ri
n i
n 1
ri
n 1
) / 2h
这个计算方法称为 Verlet 方法. 具体步骤如下: 0 1 (1)指定初始的ri , ri n F i (2)计算在n步粒子受到的力 n 1 n n 1 n 2 n 1 r 2 r r F r i i i h /m i (3)计算在n+1步的位置 :i
经典分子动力学方法就是通过直接数值解 N 个分子 的运动方程得到这些分子每个时刻的坐标与动量， 即相空间中的轨迹，从而用统计方法计算出多体系 统的宏观性质。
分子：微观客体（不一定是分子）. 运动方程：经典牛顿方程或哈密顿方程. 数值积分：欧拉法、龙格--库塔法、辛算法等. 相互作用：Lennard--Jones 势等. 初始条件：有序、无序等. 边界条件：周期性、全反射、漫反射等. 经典MD 平衡态MD 微正则系综（NVE) 正则系综（NVT） 等温等压系综（NPT） 等焓等压系综（NPH） 非平衡态 MD 量子MD（密度泛函分子动力学，第一原理计算）
结果如图示 . 而在取 h=0.005 时 , 迭代 10,000 步后 , 能量变为
E 0.99999956
二阶法解谐振子1000步的结果（h=0.05）
2) 数值计算的误差与稳定性 在数值计算中的误差可分为两类: (a) 截断误差 ( 或方法误差 ). 即由于计算方法引起的 误差 , 这个误差即使计算机有无穷位有效数字也会 存在,例如上面的 Tayler 展开中的O(h3 ) 就是展开到 2项时的截断误差. 要减小截断误差, 必须改进算法. (b) 舍入误差 ( 或计算误差 ). 在一定算法下 , 用计算 机去实现时, 由于机器字长有限,数据在计算机里的 表示会有误差 , 在运算过程中还会产生新的误差 . 随着运算次数的增加, 舍入误差可能会累积得很大. 如果舍入误差引起的后果在 N 时是有限的, 那么 这个算法是稳定的, 否则是不稳定的.不稳定的算法 是不能用的 . 为了减小舍入误差可以采用双精度 , 另外要避免相近数相减,并注意运算次序.
n n 1 n 1 v ( r r ) / 2 h i i i (4)计算速度
必须注意数值计算的稳定性 . 实际上 , 不可能精 确知道对应于给定能量的初始条件,故需要调节系统 能量达到给定值 :先计算动能和位能 , 如果总能量不 等于指定值,就改变速度,即乘以一个标度(scaling) 因子,再反复调节直到系统达到平衡.最后,还要检验 粒子的速度分布是否符合麦克斯韦分布.
在微正则系综模拟中守恒量有: i ) 2 / 2 V (ri ) 能量 E i m(r P p 动量 i i 角动量 M i ri pi 采用周期性边界条件: A( x ) A( x nL), 其中 n (n1 , n2 , n3 ) , ni 为整数.
2.单粒子问题---初值问题的数值解法 2.1.一个例子---一维简谐振子
我们先通过简谐振子的例子看看如何用计算机数值解单粒子问题 , 一维简谐振子的哈密顿量为
p 1 2 H kx 2m 2
初始条件为x(0), p(0) ，哈密顿量（能量）为守恒量。 运动方程可用拉格朗日形式和哈密顿形式，等价，但计算方法有 区别，哈密顿形式对时间一阶导数： dx H p dt p m
Euler法解谐振子的结果 (a)h=0.05, (b)h=0.005.
2.2.初值问题的数值解法和误差问题 1)初值问题的数值解法 初值问题最直接的解法可从 Taylor 展开得出到2级：
dx h 2 x(t h) x(t ) h dt 2 dx h 2 x(t h) x(t ) h dt 2 d 2x 3 O ( h ) 2 dt d 2x 3 O ( h ) 2 dt
E 0.9999459
而在取h=0.005时, 迭代10,000步后, 能量变为
E 0.9998397
为什么在取小的 h 值后结果反而不好了呢 ? 这是因为 计算机精度不够 , 而步数增大后累积的舍入误差增 大了. 改用双精度后 ,取 h=0.05时 , 迭代1,000步后 能量变为
E 0.99995860
3.分子动力学模拟的基本步骤 3.1.模型设定
, r 两体势：最简单的是硬球势: V (r ) 0, r （早期常用） 现在常用的是Lenard-Jones势：
12 6 V (r ) 4 r r
1014 s,
在分子动力学模拟中 ,宏观物理量的测量是通过对 t 1 N N A lim A ( r ( t ), p (t ))dt 相空间轨迹的平均来得到的: t t t 0 t 例如,动能为 n N 1 1
2
dp H kx dt x
用数值积分计算在相空间中的轨道 ( x(t ), p(t )) 方法是将微分变成差分，引入时间步子 h 作泰勒展 开 df
dt
[ f (t h) f (t )] / h
dx [ x(t h) x(t )] p(t ) 于是，哈密顿方程变为 dt h m dp [ p(t h) p(t )] kx(t ) dt h
3.3趋于平衡 在分子动力学模拟计算中 , 时间步子 $h$ 的选择是 很重要的 . 为了减小误差 ,h 必须取得小些 ; 但如果 h 取 得太小,趋于平衡的时间就要长,因此,为了更快地趋于 平衡 ,h 又应该取得大些 , 因而 , 在相空间中走的步子大 一些.对几百个Ar分子,采用L--J位势,发现取 h 102 就可得到好的相图结果.h是无量纲的,实际时间为
用中心差分法,有 所以
d 2 ri (t ) 2 [ r ( t h ) 2 r ( t ) r ( t h )] / h i i i dt 2
ri (t h) 2ri (t ) ri (t h) Fi (t )h 2 / m
t n nh, ri n ri (t n ), Fi n Fi (t n )
分子动力学方法 （ＭＤ）
1.引言
什么是分子动力学方法(MD)? 确定论的模拟方法实际上在本书第一章介绍的的第 一个计算机模拟工作 , 即费米 -- 巴斯塔 -- 乌拉姆实 验中已经产生了 . 但是在统计物理中所用的分子动 力学方法则产生于 1957 年 . B.J.Alder 等首先用此 方法研究了硬球气体的状态方程.
12 6 例:单原子(Ar)粒子系统的MD计算. u (rij ) 4 r 取256个粒子,位势取为L-J势: r ij ij
为位势的极小值(取 为能量单位),其位置在
r 2(1/ 6)
dx 1000( x t ) 2 2t 例:用 Euler 法解如下初值问题 dt x(0) 0
计算到x(1) , 取 h 10m , m 0,1,2,3, , 结果如下:
可见要取m=3,4,… ,才能得到稳定的结果.在实际数值 计算中,必须十分注意舍入误差与稳定性问题.
2
两式相加相减分别得到
F (t ) x(t h) x(t h) 2 x(t ) h m x(t h) x(t h) v(t ) 2h
于是，通过递推求得数值解，称为Verlet方法。两步法 还有龙格－库塔方法和辛算法（冯康先生的贡献）
考虑二级项后，谐振子的数值解在取h=0.05 时, 迭 代1,000步后能量变为
T
1
Ek
4.微正则系综的分子动力学模拟
pi 1 u (rij ) 一个多粒子系统的哈密顿量为 H 2 i m i j
2
系统受到的限制为
E const. V const.
N const. P0
牛顿方程
Fi (rij ) d 2 r (t )i 2 dt m i j
即
hp(t ) x(t h) x(t ) m p(t h) p(t ) hkx(t )
这是对 x(t),p(t) 的一组递推方程，给定了初 始条件x(0),p(0) ，就可由 t 0 h 2h 3h 一步步计算出粒子在相空间中的轨迹，这个解 法称为 Euler 方法. 由于 h 是有限量，计算总有误差，正确选择 h 是 重要的 . 图左图右分别给出了取 h=0.05 和 h=0.005 的轨迹.可见当h 取得太大时，偏离太大，以至能量 都不守恒了 :E_0=1, 在 1000 步后变为 E=12.14 ！右 图所取 h 小一个数量级，在同样的时间后能量变为 E=1.28 ，轨迹是闭合的，但结果仍不够好 , 必须改 进计算方法，进一步考虑二级项.
取周期性边界条件 , 并采用最小像力约定 . 采用自然单位，长度与 2 1/ 2 时间单位分别为 (m / 48 ) 与 .方程为无量纲形式 ,对Ar原 子,时间单位为 3 1012 s * * ( T , ) (2.53,0.636) 和(0.722,0.83134) 研究相图上的两个点: 相应的MD模拟盒子分别为L=7.38 和L=6.75. 作为初始条件 , 取原子处于面心立方的格子上 , 而速度按麦克斯韦 分布.位势的截断取两个不同的值, rc 2.5 和 rc 3.6 ,以便观察其影响.在平衡化过程中,调节速度的 1/ 2 * T ( N 1) scaling 因子 2
16i vi
重复这个步骤 , 直到达到所希望的能量 ( 或等价的平均温 度).在此MD模拟中,平衡化用了1000个MD步子, 在图中给出了 动能、位能和总能的演化过程.平均速度应当为 v 1.13 T * / 24 0.3668 而模拟的结果为 v 0.3654 ,考虑到只有256 个粒子,这个符合是相当好的.
0
Ek
n n0
2
n0 i
m(vi )
2
d 由动能可根据能量均分定理得出温度 2 Nk B 压强可通过计算在dA面积元上的 F S S 净动量转移 dFn p 的平均值得到: , 或从维里定理得到 : dA 1 N pV Nk BT ri Fi d i
MD与MC的比较 分子动力学方法实现了 Boltzmann 的统计力 学途径 ,而蒙特卡罗方法实现了 Gibbs 的统计力学 途径. 分子动力学方法可处理与时间有关的过程，因 而可处理非平衡态，是更有力的模拟方法，但计算 量大、占内存多、程序也较复杂; 蒙特卡罗方法则具有程序简单、占内存少的优 点,但难以处理非平衡态.
在分子动力学模拟中 , 在考虑粒子间相互作用时 , 通常采用最 小像力约定,如图所示. 在无穷重复的基本单元中,考虑一个 粒子和它最近的粒子或像粒子之间的作用,即
rij min( ri rj nL)
选取L,使得 rij L / 2 , 即大于L/2的 相互作用可忽略. 3.2.初条件： 对位置和速度作初始化, 原则上初始条件并不重要, 只要时间足够长就会忘掉初始条件的 , 但在实际模拟中 , 初始 条件选得好可以更快趋于平衡.可选择：位置在正规格子上, 速度随机；位置随机地偏离格子 ,速度为0；位置速度都随机.
换一个写法:
因此,可写成 ri n1 2ri n ri n1 Fi n h 2 / m 在以上式子中,可以由 ri0 , ri1 得到 ri 2 , ri3 , ,而粒子的速度为
v (ri
n i
n 1
ri
n 1
) / 2h
这个计算方法称为 Verlet 方法. 具体步骤如下: 0 1 (1)指定初始的ri , ri n F i (2)计算在n步粒子受到的力 n 1 n n 1 n 2 n 1 r 2 r r F r i i i h /m i (3)计算在n+1步的位置 :i
经典分子动力学方法就是通过直接数值解 N 个分子 的运动方程得到这些分子每个时刻的坐标与动量， 即相空间中的轨迹，从而用统计方法计算出多体系 统的宏观性质。
分子：微观客体（不一定是分子）. 运动方程：经典牛顿方程或哈密顿方程. 数值积分：欧拉法、龙格--库塔法、辛算法等. 相互作用：Lennard--Jones 势等. 初始条件：有序、无序等. 边界条件：周期性、全反射、漫反射等. 经典MD 平衡态MD 微正则系综（NVE) 正则系综（NVT） 等温等压系综（NPT） 等焓等压系综（NPH） 非平衡态 MD 量子MD（密度泛函分子动力学，第一原理计算）
结果如图示 . 而在取 h=0.005 时 , 迭代 10,000 步后 , 能量变为
E 0.99999956
二阶法解谐振子1000步的结果（h=0.05）
2) 数值计算的误差与稳定性 在数值计算中的误差可分为两类: (a) 截断误差 ( 或方法误差 ). 即由于计算方法引起的 误差 , 这个误差即使计算机有无穷位有效数字也会 存在,例如上面的 Tayler 展开中的O(h3 ) 就是展开到 2项时的截断误差. 要减小截断误差, 必须改进算法. (b) 舍入误差 ( 或计算误差 ). 在一定算法下 , 用计算 机去实现时, 由于机器字长有限,数据在计算机里的 表示会有误差 , 在运算过程中还会产生新的误差 . 随着运算次数的增加, 舍入误差可能会累积得很大. 如果舍入误差引起的后果在 N 时是有限的, 那么 这个算法是稳定的, 否则是不稳定的.不稳定的算法 是不能用的 . 为了减小舍入误差可以采用双精度 , 另外要避免相近数相减,并注意运算次序.
n n 1 n 1 v ( r r ) / 2 h i i i (4)计算速度
必须注意数值计算的稳定性 . 实际上 , 不可能精 确知道对应于给定能量的初始条件,故需要调节系统 能量达到给定值 :先计算动能和位能 , 如果总能量不 等于指定值,就改变速度,即乘以一个标度(scaling) 因子,再反复调节直到系统达到平衡.最后,还要检验 粒子的速度分布是否符合麦克斯韦分布.
在微正则系综模拟中守恒量有: i ) 2 / 2 V (ri ) 能量 E i m(r P p 动量 i i 角动量 M i ri pi 采用周期性边界条件: A( x ) A( x nL), 其中 n (n1 , n2 , n3 ) , ni 为整数.
2.单粒子问题---初值问题的数值解法 2.1.一个例子---一维简谐振子
我们先通过简谐振子的例子看看如何用计算机数值解单粒子问题 , 一维简谐振子的哈密顿量为
p 1 2 H kx 2m 2
初始条件为x(0), p(0) ，哈密顿量（能量）为守恒量。 运动方程可用拉格朗日形式和哈密顿形式，等价，但计算方法有 区别，哈密顿形式对时间一阶导数： dx H p dt p m
Euler法解谐振子的结果 (a)h=0.05, (b)h=0.005.
2.2.初值问题的数值解法和误差问题 1)初值问题的数值解法 初值问题最直接的解法可从 Taylor 展开得出到2级：
dx h 2 x(t h) x(t ) h dt 2 dx h 2 x(t h) x(t ) h dt 2 d 2x 3 O ( h ) 2 dt d 2x 3 O ( h ) 2 dt
E 0.9999459
而在取h=0.005时, 迭代10,000步后, 能量变为
E 0.9998397
为什么在取小的 h 值后结果反而不好了呢 ? 这是因为 计算机精度不够 , 而步数增大后累积的舍入误差增 大了. 改用双精度后 ,取 h=0.05时 , 迭代1,000步后 能量变为
E 0.99995860
3.分子动力学模拟的基本步骤 3.1.模型设定
, r 两体势：最简单的是硬球势: V (r ) 0, r （早期常用） 现在常用的是Lenard-Jones势：
12 6 V (r ) 4 r r