计算物理第六章

1、第六章动力学方法Molecular Dynamics Simulations第一节引言计算机模拟分类:(1)随机模拟方法。优点:随机模拟方法计算的程序简单，占内存少，但是该方法难于处理非平衡态的问题。该方法可以被用到 没有任何内禀动力学模型体系的模拟上。(2)动力学方法（MolecularDynamics或简称MD）。可以处理非平衡态问题。但是使用该方法的程序较复 杂，计算量大，占内存也多。动力学方法利用相空间中的轨迹。素。古典力学来计算许多在动力学方法中不存在任何随机因动力学简介v动力学方法（MolecularDynamics或简称MD）是、水平上求解多体问题的重要的计算机模拟方法。v动力学

2、方法可以用于模拟与粒子运动路径相关的基本过程。v在动力学中，粒子的运动行为是利用古典力学来计算许多原子、在相空间中的轨迹。系统的动力学机制决定运动方程的形式在动力学方法处理过程中，方程组的建立是通过对物理体系的微观数学描述给出的。每个都各自服从经典的力学。这种方法可以处理与时间有关的过程，因而可以处理非平衡态问题。但是使用该方法的程序较复杂，计算量大，占内存也多。适用范围广泛原则上，动力学方法所适用的微观物理体系并无什么限制。这个方法适用的体系既可以是少体系统，也可以是多体系统；既可以是点粒子体系，也可以是 具有内部结构的体系；处理的微观客体既可以是也可以是其它的微观粒子。，实际使用的限制实际

3、上，动力学模拟方法和随机模拟方法一样都着两个基本限制：(1)有限观测时间的限制；(2)有限系统大小的限制。通常人们感的是体系在热力学极限下（即粒子数目趋于无穷时）的性质。但是计算机模拟的体系大小要比热力学极限小得多，因此可能会出现有限 效应。为了减小有限效应，人们往往引入周期性、全反射、漫反射等边界条件。当然边界条件的引入显然 会影响体系的某些性质。计算机平均值之计算1rrt()( )()òA = lim, p t=dtA r tAt ®¥ t0温度之计算1 mv2= 1 k TB22第二节动力学基础知识粒子运动方程及其数值求解：1.定义Lagrangian函数为

4、则运动的Lagrangian方程为d ¶ L¶ L= 0t d单原子的运动方程p=r&iim=if i=-Ñmi &r&iVr i数值求解时的离散化方法对体系的运动方程组采用计算机进行数值求解时，需要将运动方程离散化为有限差分方-程。常用的求解方法有欧、法等。数值计算的误差阶数显然取决于所采用的数值求解方法的近似阶数。原则上， 只要计算机计算速度足够大，内存足够多，我 们可以使计算误差足够小。Verlgorithmr(t+t) = r(t) + v(t)t + (1/2)a(t)t2(1)r(t-t)= r(t) v(t)t+ (1/2)a

5、(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 a(t+t) = a(t) + b(t)(4)将 (5) 式的b(t) 代入 (4) ，得到：v(t+t) = v(t) + (1/2)a(t) + a(t+t) t(6)其他求解算法：Leap-frog algorithmr(t+t)= r(t) + v(t+(1/2)t) tv(t+(1/2)t) = v(t-(1/2)t) + a(t) tBeemans algorithmr(t+t) = r(t) + v(t)

6、t + (2/3)a(t)t2 (1/6)a(t-t)t2v(t+t) = v(t) + v(t)t + (1/3)a(t)t +(5/6)a(t)t(1/6)a(tt)t第三节动力学模拟的基本步骤动力学模拟的实际步骤可以划分为四步：设定模拟所采用的模型；给定初始条件；趋于平衡的计算过程；宏 观物理量的计算。1模拟模型的设定硬球势ì +¥ ,0 ,如果r < s ,U (r ) =í如果r ³î.Lennard-Jones型势位势V(r)43力F(r)2éæ s ö12æ s ö6 

7、49;1U (r ) = 4e êç- ç÷÷úúû0êëè r øè r ø-1-2排斥力吸引力 -30.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6v根据经典物理学的规律我们就可以知道在系综模拟中的守恒量。微正则系综的模拟中能量、动量和角动量均为守恒 量。在此系综中他们分别表示为：vm(r ) + V (r )ùé1P = å piM = å r ´ på

8、2&E =êúiiiië2ûiii元胞周期性边界条件,该边界条件的具体实现是这样操作的：当有一个粒子穿过基本动力学元胞的六方体表面时，就让这个粒子以相同的速度穿过此表面对面的表面重新进入动力学元胞内。rc不同动力学元胞盒子内粒子间的相互作用：对于不同动力学元胞盒子内粒子间的相互作用，如果相互作用是短程力，我们可以在长度rc处截断。这里V(rc)必须要足够小，以使截断显著地影响模拟结果。型力学元胞尺度L通常选得比rc大很多。我们往往选择元胞尺度满足不等式条件L/2> rc ，使得距离大于L/2的粒子的相互作用可以忽略，以避免效应。通常L的数值

9、应当选得很大。有限rr - rr+ nrLr= min粒子间相互作用的最小像力约定,ijij动力学模拟的最小像力约定示意图采用最小像力约定后，元胞内第 i 个粒子与周围粒子的相互作用势和相互作用力为Ui ( )=ri ¹ j( )åRu rijj =1, Nr (r)=i ¹ jå F ( )FiRrijrij.j =1, NR = rr , rr ,.rr,表示元胞内所有粒子的坐标。 rij是沿12Nrj - ri方向的矢量采用最小像力约定会使得在截断处粒子的受力有一个d -函数的奇异性，这会给模拟计算带来误差。为减小这种误差，我们总可以将相互作用势能

10、移到V (r ) - V (rc )，以保证在截断处相互作用接近零。2.给定初始条件给定粒子的初始位置和速度的数值：（1） 令初始位置在差分网格格子上，初始速度从玻尔兹曼分布随机抽样得到。（2） 令初始位置随机地偏离差分网格格子，初始速度为零。（3） 令初始位置随机地偏离差分网格格子，初始速度从玻尔兹曼分 布随机抽样得到。3趋于平衡使系统达到平衡，模拟中需要一趋衡过程。在这个过程中，我们增加或从系统中移出能量，直到系统具有所要求的能量。4宏观物理量的计算A(rr ( N ) ( ), pr ( N ) ( ) 1(t ¢ - t0 )¢tòA = limdt

11、62;®¥t0dtE ( p (N ) (t )t¢r1(t ¢ - t0 )òtE= limkkt¢®¥0( p 2 )( m )1nNååiE=>kn - n2mm >n0 i=10u()n=1åå( m )Urijn - nm >n0 i< j0第四节 平衡态动力学模拟1.微正则系综的动力学模拟粒子数恒定、体积恒定、能量恒定、整个系统的总动 量恒等于零。(N,V,E)动力学模拟步骤如下（Verlet算法）：（1）给定初始空间位置。（2）计算在第n

12、步时粒子所受的力。（3）计算在第n+1步时所有粒子所处的空间位置。r (n)2r (n+1) = 2r (n) - r (n-1)+ Fih/ miii（4）计算第n步的速度。=( n +1)- r( n -1)v ( n)r/ 2hiii（5）返回到步骤（2），开始的模拟计算。Verlet算法的速度形式：r(1)r（1）（2）给定初始空间位置。ir(1)v给定初始速度。i（3）利用公式：rr ( n+1) = r ( n) + hv ( n)+ Fih / 2m( n)2iii( n +1)计算在第n+1步时所有粒子所处的空间位置r。i（4）计算在第n+1步时所有粒子的速度：+ h(+ Fi

13、)/ 2mrrr( n +( n += v1)( n )1)( n )vFiii（5）返回到步骤（3），开始第n+2步的模拟计算。速度标度因子:ù1 / 2éb = êT (N - 1) ú*êúå i216vêúëûi2.正则系综的动力学模拟ù1/ 2é速度标度因子:b = ê(3N - 4)kT úêúåi2mvêúëiû正则系综动力学的模拟具体步骤:(Verlet算法的速度形式)（1）给定初始空间位置，（2）给定初始速度，（3）利用公式：rrr(n+= r+(n) + F (n) h 21)( n)rhv/ 2miiii计算在第n+1步时所有粒子所处的空间位置，（4） 计在时所有粒子的速度：+ h()/ 2mrr( n +( n += v+1)( n )1)( n )vFFiiii动能和速度标度因子：ù1 / 2&#