分子动力学

第六章 分子动力学方法 Molecular Dynamics Simulations 第一节 引言 计算机模拟分类: (1)随机模拟方法。 优点: 随机模拟方法计算的程序简单，占内存少，但 是该方法难于处理非平衡态的问题。 (2)分子动力学方法（Molecular Dynamics或简称MD） 。 可以处理非平衡态问题。但是使用该方法的程序较复 杂，计算量大，占内存也多。 分子动力学方法利用牛顿古典力学来计算许多分子在 相空间中的轨迹。 分子动力学简介 1. 分子动力学是在原子、分子水平上求解多体问题的 重要的计算机模拟方法。 2. 通过求解所有粒子的运动方程，分子动力学方法可 以用于模拟与粒子运动路径相关的基本过程。 3. 在分子动力学中，粒子的运动行为是通过经典的 Newton运动方程所描述。 第二节 粒子运动方程的数值求解 粒子体系的运动方程⎯⎯Lagrangian方程 1. 定义Lagrangian函数为 则运动的Lagrangian方程为 0 d d =∂ ∂−∂ ∂ kk qqt LL & 粒子的运动方程-Lagrangian方程 k k k k k k k k k q p p q pq q p ∂ ∂−=∂ ∂= −= ∂ ∂= ∑ HH LH L & && & ),(),( qqqp 单原子的牛顿运动方程 V iiii i i i m m rfr pr −∇== = && & Verlet algorithm r(t+∆t) = r(t) + v(t)∆t + (1/2)a(t)∆t2 (1) r(t-∆t) = r(t) – v(t)∆t + (1/2)a(t)∆t2 (2) 将上面两式相加，得到： r(t+∆t) = 2r(t) – r(t- ∆t) + a(t)∆t2 (3) v(t+∆t) = v(t) + a(t)∆t + (1/2)b(t)∆t2 (4) a(t+∆t) = a(t) + b(t)∆t (5) 将 (5) 式的b(t) 代入 (4) ，得到： v(t+∆t) = v(t) + (1/2)[a(t) + a(t+∆t)] ∆t (6) 其他求解算法： • Leap-frog algorithm r(t+∆t) = r(t) + v(t+(1/2)∆t) ∆t v(t+(1/2)∆t) = v(t-(1/2)∆t) + a(t) ∆t • Beeman’s algorithm r(t+∆t) = r(t) + v(t)∆t + (2/3)a(t)∆t2 – (1/6)a(t-∆t)∆t2 v(t+∆t) = v(t) + v(t)∆t + (1/3)a(t)∆t + (5/6)a(t)∆t–(1/6)a(t∆t)∆t 第三节 分子动力学模拟的基本步骤 分子动力学模拟的实际步骤可以划分为四步：设定模拟 所采用的模型；给定初始条件；趋于平衡的计算过程；宏 观物理量的计算。 1．模拟模型的设定 硬球势 Lennard-Jones型势 . , σ σ ≥ < r r 如果 如果( ) ⎩⎨ ⎧ ∞+= ,0 , rU ( ) ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ ⎟⎠ ⎞⎜⎝ ⎛−⎟⎠ ⎞⎜⎝ ⎛= 612 4 rr rU σσε 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 -3 -2 -1 0 1 2 3 4 位势V(r) 力F(r) 吸引力排斥力 ™ 根据经典物理学的规律我们就可以知道在系综模拟 中的守恒量。 ™ 微正则系综的模拟中能量、动量和角动量均为守恒 量。在此系综中他们分别表示为： 元胞 周期性边界条件, ( ) ( )∑ ⎥⎦⎤⎢⎣⎡ += i ii rVrmE r&r 221 ∑= i ipP rr ∑ ×= i ii prM rrr 分子动力学模拟的最小像力约定示意图 ( )Lnrrr jiij rrr +−= min最小像力约定, 2.给定初始条件 给定粒子的初始位置和速度的数值： （1）令初始位置在差分网格格子上，初始速度从玻尔兹曼分布随机 抽样得到。 （2）令初始位置随机地偏离差分网格格子，初始速度为零。 （3）令初始位置随机地偏离差分网格格子，初始速度从玻尔兹曼分 布随机抽样得到。 3．趋于平衡 使系统达到平衡，模拟中需要一个趋衡过程。在这个过程中，我 们增加或从系统中移出能量，直到系统具有所要求的能量。 4．宏观物理量的计算 => ( ) ( ) ( )( )∫ ′∞→′ −′= tt NNt prAdttA 0 }{},{1lim )()(0 τττ rr ( ) ( )( )∫ ′∞→′ −′= tt Nktk pEdttE 0 }{1lim )(0 ττ r ∑∑ > =−= n n N i i k m p nn E 0 1 )(2 0 2 )(1 μ μ ( )∑∑ > <−= n n ji ijrunn U 0 )( 0 1 μ μ 第四节 平衡态分子动力学模拟 1.微正则系综的分子动力学模拟 粒子数恒定、体积恒定、能量恒定、整个系统的总动 量恒等于零。 分子动力学模拟步骤如下（分子动力学模拟步骤如下（VerletVerlet算法）算法）： （1）给定初始空间位置。 （2）计算在第步时粒子所受的力。 （3）利用如下公式，计算在第步时所有粒子所处的空 间位置。 （4）计算第步的速度。 （5）返回到步骤（2），开始下一步的模拟计算。 mhFrrr ni n i n i n i /2 2)()1()()1( rrrr +−= −+ ( ) hrrv ninini 2/)1()1()( −+ −= rrr Verlet算法的速度形式： （1） 给定初始空间位置 。 （2） 给定初始速度 。 （3） 利用公式： 计算在第n+1步时所有粒子所处的空间位置 。 （4） 计算在第n+1步时所有粒子的速度： （5） 返回到步骤（3），开始第n+2步的模拟计算。 { })1(ivr { })1(irr mhFvhrr ni n i n i n i 2/ 2)()()()1( rrrr ++=+ { })1( +nirr ( ) mFFhvv nininini 2/)()1()()1( rrrr ++= ++ 速度标度因子: 2.正则系综的分子动力学模拟 速度标度因子: 正则系综分子动力学的模拟具体步骤: (Verlet算法的速度形式) 2/1 2 * 16 )1( ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ −= ∑ i iv NTβ 2/1 2 )43( ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ −= ∑ i imv kTNβ （1）给定初始空间位置，， （2）给定初始速度， （3）利用公式： 计算在第n+1步时所有粒子所处的空间位置， （4） 计算在第步时所有粒子的速度： 动能和速度标度因子： （5） 计算将速度乘以标度因子的值，并让该值作为下一次计算时， 第n+1步粒子的速度： 。 （6） 返回到步骤（3），开始第n+2步的模拟计算。 mhFvhrr ni n i n i n i 2/ 2)()()()1( rrrr ++=+ ( ){ }mFFhvv nininini 2/)()1()()1( rrrr ++= ++ ( )2)1( 2 1 ∑ += i n ik vmE 2/1 2)1( )( )43( ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ −= ∑ + i n ivm kTNβ { } { }11 ++ → nini vv rr β Thank You!