分子动力学模拟方法

1、 1957年年：基于刚球势的分子動力学法（：基于刚球势的分子動力学法（Alder and Wainwright） 1964年年：利用：利用Lennard-Jone势函数法对液态势函数法对液态氩氩性质的模拟（性质的模拟（Rahman） 1971年年：模拟具有分子：模拟具有分子团簇团簇行为的行为的水水的性质（的性质（Rahman and Stillinger） 1977年年：约束动力学方法（：约束动力学方法（Rychaert, Ciccotti van Gunsteren） 1980年年：恒压条件下的动力学方法（：恒压条件下的动力学方法（Andersen法、法、Parrinello-Rahman法

2、）法） 1983年年：非平衡态动力学方法（：非平衡态动力学方法（Gillan and Dixon） 1984年年: 恒温条件下的动力学方法恒温条件下的动力学方法(Berendsen et al.) 1984年年：恒温条件下的动力学方法（：恒温条件下的动力学方法（Nos-Hoover法）法） 1985年年：第一原理分子動力学法（：第一原理分子動力学法（Car-Parrinello法）法） 1991年年：巨正则系综巨正则系综的分子动力学方法（的分子动力学方法（Cagin and Pettit） 分子动力学简史分子动力学简史 粒子的运动取决于经典力学粒子的运动取决于经典力学 (牛顿定律（牛顿定律（F

3、=ma） 原理：原理： 计算一组分子的相空间轨道，其中每个分子各自服从计算一组分子的相空间轨道，其中每个分子各自服从 牛顿运动定律：牛顿运动定律： 1 111 2 )( 2 1 N i N ij ij N i i i rU m p H ii i ii m dt d mv r p 1 11 1 11 )( )( N i N ij ij N i N ij ij i rU r dt d ij r F p )0( 0iti rr )0( 0it i dt v dr 初始条件：初始条件： 粒子位置的粒子位置的Taylor展开式：展开式： 粒子位置粒子位置 ： 粒子速度粒子速度 ： t (t) t)-(t

4、(t)2 t)(t 2 i arrr iii t2 t)-(tt)(t (t) ii i rr v 粒子加速度：粒子加速度： 3 i 2 ii t (t) 6 1 t (t) 2 1 t (t) (t) t)(tbavrr ii 3 i 2 ii t (t) 6 1 t (t) 2 1 t (t) (t) t)(tbavrr ii (t) (t) i m i i F a 开始运动时需要开始运动时需要r(t-t)： t (0) (0) t)( i vrr + n 算法启动算法启动 u 规定初始位置规定初始位置 u 规定初始速度规定初始速度 u 扰动初始位置：扰动初始位置： u 计算第计算第n步的

5、力步的力 u 计算第计算第n+1步的位置：步的位置： u 计算第计算第n步的速度：步的速度： u 重复重复至至 t (0) (0) t)( i vrr t (t) t)-(t(t)2 t)(t 2 i arrr iii t2 t)-(tt)(t (t) ii i rr v Do 100 I = 1, N RXNEWI = 2.0 * RX(I) RXOLD(I) + DTSQ * AX(I) RYNEWI = 2.0 * RY(I) RYOLD(I) + DTSQ * AY(I) RZNEWI = 2.0 * RZ(I) RZOLD(I) + DTSQ * AZ(I) VXI = ( RXNE

6、WI RXOLD(I) ) / DT2 VYI = ( RYNEWI RYOLD(I) ) / DT2 VZI = ( RZNEWI RZOLD(I) ) / DT2 RXOLD(I) = RX(I) RYOLD(I) = RY(I) RZOLD(I) = RZ(I) RX(I) = RXNEWI RY(I) = RYNEWI RZ(I) = RZNEWI 100 CONTINUE t (t) t) 2 1 -(t t) 2 1 (t i avv ii t t) 2 1 (t(t) t)(t i vrr ii 2 t) 2 1 -(tt) 2 1 (t (t) ii i vv v 1. 首先利

7、用当前时刻的加速度，计算半个时间步长后的速度：首先利用当前时刻的加速度，计算半个时间步长后的速度： 2. 计算下一步长时刻的位置：计算下一步长时刻的位置： 3. 计算当前时刻的速度：计算当前时刻的速度： t-t/2tt+t/2t+tt+3t/2t+2t v r v 开始运动时需要开始运动时需要v(-t/2)： t/2 (0) (0) t/2)( i avv n 算法启动算法启动 u 规定初始位置规定初始位置 u 规定初始速度规定初始速度 u 扰动初始速度：扰动初始速度： u 计算第计算第n步的力步的力 u 计算第计算第n+1/2步的速度：步的速度： u 计算第计算第n+1步的位置：步的位置：

8、u 计算第计算第n步的速度：步的速度： u 重复重复至至 t/2 (0) (0) t/2)( i avv t (t) t) 2 1 -(t t) 2 1 (t i avv ii t t) 2 1 (t(t) t)(t i vrr ii 2 t) 2 1 -(tt) 2 1 (t (t) ii i vv v 2 ii t (t) 2 1 t (t) (t) t)(tavrr ii tt)(t (t) 2 1 (t) t)(t ii aavv ii t(t) 2 1 (t) t) 2 1 (t i avv ii tt)(t 2 1 t) 2 1 (t t)(t i avv ii t t) 2 1

9、(t (t) t)(t i vrr ii 等价于等价于 n 算法启动算法启动 u 规定初始位置规定初始位置 u 规定初始速度规定初始速度 u 计算第计算第n+1步的位置：步的位置： u 计算第计算第n+1步的力步的力 u 计算第计算第n+1步的速度：步的速度： u 重复重复至至 2 ii t (t) 2 1 t (t) (t) t)(tavrr ii tt)(t (t) 2 1 (t) t)(t ii aavv ii Verlet Leap- frog Velocity Verlet 预测预测(Predictor)阶段：其基本思想是阶段：其基本思想是Taylor展开，展开， )()( )()(

10、)( )( 2 1 )()()( )( 6 1 )( 2 1 )()()( 2 32 ttt ttttt ttttttt ttttttttt p p p p bb baa bavv bavrr 根据新的原子位置根据新的原子位置r rp p，可以计算获得校正后的，可以计算获得校正后的a ac c( (t t+ + t t),),定义预测误差定义预测误差: : )()()(tttttt pc aaa )()()( )()()( )()()( )()()( 3 2 1 0 ttcttt ttcttt ttcttt ttcttt pc pc pc pc abb aaa avv arr 利用此预测误差，

11、对预测出的位置、速度、加速度等量进行校正：利用此预测误差，对预测出的位置、速度、加速度等量进行校正： 校正校正(Corrector)阶段：阶段： u预测阶段运动方程的变换：预测阶段运动方程的变换： 定义一组矢量定义一组矢量: 3 3 3 3 2 2 2 2 1 0 6 1 d )(d 2 1 d )(d d )(d )( t t t t t t t t t t r r r r r r rr )( )( )( )( 1000 3100 3210 1111 )( )( )( )( 3 2 1 0 3 2 1 0 t t t t tt tt tt tt p p p p r r r r r r r r

12、 )()( )()()( )( 2 1 )()()( )( 6 1 )( 2 1 )()()( 2 32 ttt ttttt ttttttt ttttttttt p p p p bb baa bavv bavrr u校正阶段运动方程的变换：校正阶段运动方程的变换： r r r r r r r r r )( )( )( )( )( )( )( )( 3 2 1 0 3 2 1 0 3 2 1 0 c c c c tt tt tt tt tt tt tt tt p p p p c c c c r 的形式：的形式：C0， C1， C2， C3的值以及的值以及 C0， 取决于运动方程的阶数。取决于运动

13、方程的阶数。 n 一阶运动方程：一阶运动方程： )(rrf)()( 11 tttt pc rrr n 二阶运动方程之一：二阶运动方程之一： )(rrf )()( 22 tttt pc rrr n 二阶运动方程之二：二阶运动方程之二： ),(rrr f)()( 22 tttt pc rrr 室温下，室温下， t 1 fs (femtosecond 10-15s)，温，温 度越高，度越高，t 应该减小应该减小 n 太长的时间步长会造成分子间的激烈碰撞，体系数据溢出太长的时间步长会造成分子间的激烈碰撞，体系数据溢出; ; n 太短的时间步长会降低模拟过程搜索相空间的能力太短的时间步长会降低模拟过程搜

14、索相空间的能力 u 它是分子动力学方法的最基本系综它是分子动力学方法的最基本系综 u 具有确定的粒子数具有确定的粒子数N，能量，能量E和体积和体积V u 算法：算法： u 规定初始位置和初始速度规定初始位置和初始速度 u 对运动方程积分若干步对运动方程积分若干步 u 计算势能和动能计算势能和动能 u 若能量不等于所需要的值，对速度进行标度若能量不等于所需要的值，对速度进行标度 u 重复重复至，直到系统平衡至，直到系统平衡 微正则系综微正则系综(NVE)MD模拟算法的流程图：模拟算法的流程图： 给定每个分子的初始位置给定每个分子的初始位置ri(0)和速度和速度vi(0) 计算每个分子的受力计算每

15、个分子的受力Fi和加速度和加速度ai 解运动方程并求出每个分子运动一个时间步解运动方程并求出每个分子运动一个时间步 长后到达的位置所具有的速度长后到达的位置所具有的速度 统计系统的热力学性质及其它物理量统计系统的热力学性质及其它物理量 统计性质不变？统计性质不变？ 打印结果，结束打印结果，结束 YesNo 移动所有分子到新的位置并具有当前时刻的移动所有分子到新的位置并具有当前时刻的 速度速度 微正则系综微正则系综MD模拟程序模拟程序F3讲解（讲解（LJ, NVE）：）： 无因次量：无因次量： /*rr 3 * /*EE /*ff /* 3 PP /*kTT 2/12 * )/(m t t 2/

16、1 * )/( v v m MD模拟中几个热力学量的计算：模拟中几个热力学量的计算： 对于由对于由N个单原子组成的系统：个单原子组成的系统： 动能和温度：动能和温度： TNkmE B N i K 2 3 v 2 1 1 2 i 采用对比量：采用对比量：/*kTT 2/1 * )/( v v m /*EE * 1 2 * i * 2 3 v 2 1 NTE N i K 2 1 11 4) ru 2 )(r(dr N rUUUU c r N i N ij ijclrcc 对于对于LJ流体：流体： 39 3 3 1 3 8 cc lrc rr NU 势能：势能： 4)( 612 rr rUc /*r

17、r 3 */*EE 采用对比量：采用对比量： 11 4)( 6 * 12 * * rr rU c 3 * 9 * * 11 3 1 3 8 cc rr NU lrc 内能：内能： K EUE 内能由势能和动能组成：内能由势能和动能组成： 采用对比量：采用对比量： * K EUE 压力：压力： 39 32 3 2 3 16 cc lrc rr P lrc c B P V W TkP /*rr 3 */*EE 采用对比量：采用对比量： /* 3 PP /*ff 1 11 )( 3 1 N i N ij ijcc rWW ij ij ijijijijc rU rW r rfr )( 3 1 )( *

18、 * * * lrc c P V W TP 3 * 9 * 2 * 11 3 2 3 16 cc lrc rr P 1 11 * )( 3 1 N i N ij ijcc rWW * * * )( )( j j jjjjc rU rW i i iiii r rfr 练习：练习： 推导推导LJ流体分子间力的表达式流体分子间力的表达式(fx, fy, fz及其对比量及其对比量)： 4)( 612 rr rU 势能函数形式：势能函数形式： 力：力： r )( )( rU rf 2 612 2 1 84 r r rrr r r U r U f /*rr 采用对比量：采用对比量：/*ff =x, y,

19、z 2 * * 6 * 12 * * 1 2 11 84 r r rr f * )(fr rWc LJ分子间的维里项：分子间的维里项： 2 * * 6 * 12 * * 1 2 11 84 r r rr f =x, y, z 1 2 11 84 )( 6 * 12 * * rr rWc 2 * * * )( r r rWf c t (t) t) 2 1 -(t t) 2 1 (t i m F vv ii t t) 2 1 (t(t) t)(t i vrr ii 2 t) 2 1 -(tt) 2 1 (t (t) ii i vv v (t) i a /*ff 2/12 * )/(m t t 2/

20、1 * )/( v v m 采用对比量：采用对比量： /*rr 2/1 2 * * 2/1 * 2/1 * t )(t )t 2 1 -(t )t 2 1 (t m mmm * i * i * i F vv 2/1 2/1 2 * * * t )(t )t 2 1 -(t )t 2 1 (t mm m * i * i * i F vv * t )(t )t 2 1 -(t )t 2 1 (t * i * i * i Fvv 最终得到：最终得到： * t )t 2 1 (t )(t )t(t * i * i * i vrr 2 )t 2 1 -(t)t 2 1 (t (t) * * i * i

21、* i vv v 同理得到：同理得到： m Tk v B x 2 根据能量均分原理，可知：根据能量均分原理，可知： 2/1 2 * * 2/1 2 2/1 2 2/1 333 v T vm Tk vNm TNk T T reqrefBrefBref 标度因子标度因子: 对比量对比量 速度标度速度标度: oldnew vv * old * new vv 或或 微正则系综微正则系综MD模拟程序模拟程序F3讲解（讲解（LJ, NVE）：）： 初始化：初始化： READ (*,(A) TITLE ! 运行作业题目运行作业题目 READ (*,*) NSTEP ! 运行步数运行步数 READ (*,*)

22、 IPRINT ! 打印步数打印步数 READ (*,(A) CNFILE ! 位型文件位型文件 READ (*,*) DENS ! 对比密度对比密度 READ (*,*) RTEMP ! 对比温度对比温度 READ (*,*) RCUT ! 对比截断半径对比截断半径 READ (*,*) DT ! 对比时间步长对比时间步长 CALL READCN ( CNFILE ) 初始位型：初始位型： u面心立方面心立方 (face-centered cubic, FCC)： 每面中心有一格点每面中心有一格点 u体心立方体心立方 (body-centered cubic, BCC)： u简单立方简单立方

23、 (simple cubic, SC)： XL 初始位型：初始位型：面心立方面心立方 ( FCC) (程序程序F23) NC=(REAL(N)/4.0)*(1.0/3.0) XL = 1.0 / REAL ( NC ) Y = 0.5 * XL R(1)=(0, 0, 0) R(2)=(0, Y, Y) R(3)=(Y, 0, Y) R(4)=(Y, Y, 0) M = 0 DO 10 I = 1, NC DO 10 J = 1, NC DO 10 K = 1, NC DO 11 IJ = 1, 4 RX(IJ+M)=RX(IJ) + XL*(K-1) RY(IJ+M)=RY(IJ) + XL

24、*(J -1) RZ(IJ+M)=RZ(IJ) + XL*(I -1) 11 CONTINUE M = M + 4 10 CONTINUE DO 100 I = 1, N RX(I) = RX(I) - 0.5 RY(I) = RY(I) - 0.5 RZ(I) = RZ(I) - 0.5 100 CONTINUE 将模拟盒子的中心移到原点：将模拟盒子的中心移到原点： 初始速度：初始速度： n 简单的选择：简单的选择： V random(-0.5, 0.5) =x, y, z m Tk v B 2 2/1 * * 2/1 2 * * 2/1 2 2/1 T T v T vm Tk T T re

25、qreqrefBref 标度因子标度因子: * old * new vv 速度标度速度标度: * 2 * Tv FACTOR = SQRT ( RTEMP ) DO 100 I = 1, N VX(I) = FACTOR * ( RANF(DUMMY) - 0.5 ) VY(I) = FACTOR * ( RANF(DUMMY) - 0.5 ) VZ(I) = FACTOR * ( RANF(DUMMY) - 0.5 ) 100 CONTINUE * 2 * Tv 随机安排初始速度：随机安排初始速度： 标度初始速度：标度初始速度： SUMKX = 0.0 SUMKY = 0.0 SUMKZ =

26、 0.0 DO 200 I = 1, N SUMKX = SUMKX + VX(I)*2 SUMKY = SUMKY + VY(I)*2 SUMKZ = SUMKZ + VZ(I)*2 200 CONTINUE BEITAX =SQRT(RTEMP/SUMKX) BEITAY =SQRT(RTEMP/SUMKY) BEITAZ =SQRT(RTEMP/SUMKZ) DO 300 I = 1, N VX(I) = VX(I) * BEITAX VY(I) = VY(I) * BEITAY VZ(I) = VZ(I) * BEITAZ 300 CONTINUE 2/1 * * 2/1 2 * *

27、T T v T reqreq 标度因子标度因子: * old * new vv SUMX = 0.0 SUMY = 0.0 SUMZ = 0.0 DO 200 I = 1, N SUMX = SUMX + VX(I) SUMY = SUMY + VY(I) SUMZ = SUMZ + VZ(I) CONTINUE SUMX = SUMX / REAL ( N ) SUMY = SUMY / REAL ( N ) SUMZ = SUMZ / REAL ( N ) DO 300 I = 1, N VX(I) = VX(I) - SUMX VY(I) = VY(I) - SUMY VZ(I) = V

28、Z(I) - SUMZ 300 CONTINUE 控制体系的总动量为零：控制体系的总动量为零： n 从从Maxwell分布中抽样：分布中抽样： x )(xP 0 x x dx 高斯高斯(Gauss)分布分布: 对于等几率随机试验对于等几率随机试验(Bernoulli试验试验)， 大量的试验结果满足高斯分布大量的试验结果满足高斯分布 2 2 2 )(2 exp )(2 1 x xx x xP 麦克斯韦速度分布定律麦克斯韦速度分布定律: dve kT m Ndvvf Tk mv B 2 2 1 2 2 m Tk v B 2 由于：由于： 2 2 2 )(2 exp )(2 1 x xx x xP

29、=x, y, z 单位体积的分子再每个分量上的速度分布实际上就是单位体积的分子再每个分量上的速度分布实际上就是 高斯分布。高斯分布。 n 从从Maxwell分布中抽样：分布中抽样： 高斯高斯(Gauss)分布的随机数生成方法分布的随机数生成方法: 2 2 2 )(2 exp )(2 1 x xx x xP 生成随机数：生成随机数： i，i=1, 2, , 12 46 12 1 i i R )R a R ) a R ) a )R a Ra ( 1 2 3 2 5 2 7 2 9 SUM = 0.0 DO 10 I = 1, 12 SUM = SUM + RANF ( DUMMY ) 10 CON

30、TINUE R = ( SUM - 6.0 ) / 4.0 R2 = R * R GAUSS = ( A9 * R2 + A7 ) * R2 + A5 ) * R2 + A3 ) * R2 +A1 ) * R 高斯高斯(Gauss)分布的随机数生成分布的随机数生成(程序程序F24) FACTOR = SQRT ( RTEMP ) DO 100 I = 1, N VX(I) = FACTOR * GAUSS(DUMMY) VY(I) = FACTOR * GAUSS(DUMMY) VZ(I) = FACTOR * GAUSS(DUMMY) CONTINUE 控制总动量为零：同前面一样处理。控制总

31、动量为零：同前面一样处理。 从从Maxwell分布中抽样分布中随机安排初始速度：分布中抽样分布中随机安排初始速度： 微正则系综微正则系综MD模拟程序模拟程序F3讲解讲解(LJ, NVE) ： 量纲变换：量纲变换： SIGMA = ( DENS / REAL ( N ) ) * ( 1.0 / 3.0 ) RCUT = RCUT * SIGMA DT DT * SIGMA DENS = DENS / ( SIGMA * 3 ) 333 N V N * 模拟盒子的边长为模拟盒子的边长为1 V N 3 * * L 2/1 * )/(m t t 长程校正：长程校正： 微正则系综微正则系综MD模拟程序模

32、拟程序F3讲解讲解(LJ, NVE) ： SR3 = ( SIGMA / RCUT ) * 3 SR9 = SR3 * 3 SIGCUB = SIGMA * 3 VLRC = ( 8.0 /9.0 ) * PI * DENS * SIGCUB * REAL ( N ) : * ( SR9 - 3.0 * SR3 ) WLRC = ( 16.0 / 9.0 ) * PI * DENS * SIGCUB * REAL ( N ) :* ( 2.0 * SR9 - 3.0 * SR3 ) 39 3* 3 9 8 cc lrc rr NU 39 3* 32 9 16 cc lrclrc rr NVPW

33、 算法：算法启动算法：算法启动 微正则系综微正则系综MD模拟程序模拟程序F3讲解讲解(LJ, NVE) ： CALL FORCE ( -DT, SIGMA, RCUT, NEWV, NEWVC, NEWW ) CALL MOVE ( -DT ) CALL FORCE ( -DT, SIGMA, RCUT, V, VC, W ) CALL FORCE ( DT, SIGMA, RCUT, V, VC, W ) CALL KINET ( OLDK ) CALL MOVE ( DT ) CALL FORCE ( DT, SIGMA, RCUT, NEWV, NEWVC, NEWW ) CALL KI

34、NET ( NEWK ) 算法：差分格式：算法：差分格式： SR2 = SIGSQ / RIJSQ VIJ = 4.0 * ( SR12 - SR6 ) WIJ = 24.0 * ( 2.0 * SR12 - SR6 ) VELIJ = WIJ * DT / RIJSQ DVX = VELIJ * RXIJ . VXI = VXI + DVX . VX(J) = VX(J) DVX V = V + VIJ W = W + WIJ CALL MOVE(DT) 11 4)( 6 * 12 * * c rr rU 1 2 11 84 )( 6 * 12 * * rr rWc 2 * * * )( r

35、 r rWf c * t )(t )t 2 1 -(t )t 2 1 (t * i * i * i Fvv DO 1000 I = 1, N RX(I) = RX(I) + VX(I) * DT RY(I) = RY(I) + VY(I) * DT RZ(I) = RZ(I) + VZ(I) * DT 1000 CONTINUE MOVE(DT): * t )t 2 1 (t )(t )t(t * i * i * i vrr 速度的标定（只用于平衡阶段）速度的标定（只用于平衡阶段） SUMK = 0.0 DO 200 I = 1, N SUMK = SUMKX + VX(I)*2+VY(I)*

36、2 + VZ(I)*2 200 CONTINUE BEITA =SQRT( 3.0 * RTEMP / SUMK ) DO 300 I = 1, N VX(I) = VX(I) * BEITA VY(I) = VY(I) * BEITA VZ(I) = VZ(I) * BEITA 300 CONTINUE 2/1 * * 2/1 2 * * 3 T T v T reqreq * old * new vv u 具有确定的粒子数具有确定的粒子数N，温度，温度T和体积和体积V u速度的直接标度速度的直接标度 u热浴方法热浴方法 (Andersen Thermostat) u约束方法约束方法(阻尼力方

37、法阻尼力方法) u系统扩展方法系统扩展方法(Extended Systems Method) u 问题的关键：温度的约束问题的关键：温度的约束 Nose-Hoover方法方法 一、热浴方法一、热浴方法 (Andersen Thermostat) u 引入一个与虚拟粒子碰撞的随机力引入一个与虚拟粒子碰撞的随机力 u 想象系统浸在热浴当中想象系统浸在热浴当中 u系统和热浴间的相互作用强度由随机碰撞的频率决定系统和热浴间的相互作用强度由随机碰撞的频率决定 u碰撞的几率等于碰撞的几率等于Nudt u如果一个粒子经历碰撞，它的速度将从约束温度下的如果一个粒子经历碰撞，它的速度将从约束温度下的 Maxwe

38、ll分布中随机抽取分布中随机抽取 u 总能量和总动量均不守恒总能量和总动量均不守恒 二、约束方法二、约束方法 u 是等动能是等动能(Iso-Kinetics)分子动力学方法分子动力学方法 u 系统的运动方程为：系统的运动方程为： m/pr pprfp),( u 引入阻尼系数引入阻尼系数 以保证将温度约束在恒定值以保证将温度约束在恒定值 u 根据高斯最小约束原理：根据高斯最小约束原理： i i i ii 2 p fp 三、三、Nose-Hoover扩展方法扩展方法 基本思想：基本思想： 设想原系统与一个耦合系统共同组成一个扩展系统，设想原系统与一个耦合系统共同组成一个扩展系统， 允许热流在原系统

39、和耦合系统之间交换。允许热流在原系统和耦合系统之间交换。 Q：等效质量：等效质量 S: 扩展坐标变量扩展坐标变量 ： 热力学阻尼系数热力学阻尼系数 L： 扩展系统的自由度扩展系统的自由度 u Predictor-corrector algorithm is straightforward u Verlet algorithm is feasible, but tricky to implement 积分方案：积分方案： update of depends on p update of p depends on Nos-Hoover方法正确地描述了方法正确地描述了NVT系综中的动系综中的动 量和位型，而等动能方法只正确地描述了后者。量和位型，而等动能方法只正确地描述了后者。 扩展系统的哈密顿量扩展系统的哈密顿量Hamilt