经典分子动力学方法 中国科学院固体物理研究所 计算材料科学研究室 范巍 分子动力学基本原理 ? 一个体系有N个原子 ? 体系的状态 由这N个原子的位置{r i }和动量{p i } 或速度{v i }来标志。 ? 体系的能量为H({r i ,p i }) ? 体系的运动方程为 H r p t i i ? ? ?= ? ? H p r t i i ? ? = ? ? 分子动力学的主要目的是解上面的方程求得体系状态 相空间演化的轨迹 {r i p i } t0 ,{r i p i } t1 ,{r i p i } t2 ,{r i p i } t3 ,……。 进而可计算我们感兴趣的物理量的值 Q({r i ,p i })。 在实际的应用中，我们把上面的哈密顿方程化为下面的 牛顿方程，并且用位置 r i 和速度 v i 做为描述体系的参量。 })({ 2 1 1 2 i N i ii rVvmH += ∑ = })({ 2 2 i i i dt d m rV r r i ? ? ?= V({r i })是原子间相互作用势，通过解上面的方 程我们可以得到体系在相空间得由轨迹，进而 求得物理量得平均值 [t(1),t(2),t(3),…t(M)] ),...,1....()},({ 1 )( )( MmvrQ M Q mt mtii == ∑ 积分牛顿方程的方法(I) 1.Verlet法则 ))......((})({ )()(2)( 4 2 AhOrV rh htrtrhtr m i i iii i + ? ? ?= ?+?+ ))......(( 2 )()( )( 3 BhO h htrhtr tv ii i + ??+ = *由前两个时刻的位置，根据方程 (A)推得下一个时刻得位置 速度由方程 (B)计算速度。 *需要连续记录两个时刻得位置。 积分牛顿方程的方法(II) 2.Verlet 速度法则 ).........( 2 )( )( )()( Ah m tF tv h trhtr i i ii += ?+ )........( 2 )()( )()( Bh m tFhtF tvhtv i ii ++ +=+ 需要知道上一个时刻得位置，速度和力，首先由方程(A)计算 新得位置，然后计算新得力F(t+h)，再由由方程(B)计算新时 刻得速度，需要储存前一个时刻的位置，速度和力。 积分牛顿方程的方法(III) 3.Leap Frog法则 )().........( )()( 2 1 Bhtv h trhtr i ii += ?+ )().........( )()( 2 1 2 1 AtF h htvhtv m i ii i = ??+ )........( 2 )()( )( 2 1 2 1 C htvhtv tv ?++ = 需要知道上一个时刻得位置 r i (t)，中间时刻的速度 v i (t-(1/2)h)和力F i (t)，首先由方程(A)计算 v i (t+(1/2)h),然后由方程(B)计算新时刻得位置 r i (t+h)，再由方程(C)计算 v i (t)。 积分牛顿方程的方法(IV) 4.Gear Predictor-Corrector法则 ).........( 2 )()()()( B tFhtF h tvhtv m iiii i ++ = ?+ )().........( 2 )()( Atv h htrhtr i i p i = ??+ ).........( 2 )()()()( C tvhtv h trhtr iiii ++ = ?+ 首先由方程(A)初步预测新的位置r i p (t+h),并且计算 力F i (t+h)。然后根据方程(B)计算新的速度v i (t+h). 再由方程(C)计算最后新的位置r i (t+h) 原子间相互作用势 V({r i }) ? 两体势 (Lennard-Jones,Morse) 惰性气体、简单金属 。 ? 原子内嵌势 (Embed Atom Methods) 简单金属，过渡金属 ? 类原子内嵌势 Johnson, Mei potential, Glue Potential, Finnis-Sutton Potential ? 紧束缚势 (Tight Binding) ? Tersoff,Brennerd,Stillinger-Weber Potential 半导体 两体势 Lennard-Jones Morse M.J.Weins,Surf. Sci. 31, 138(1972) L.A.Girifalco and G.V.Weizer, Phys.Rev, 114,687(1959) ])()([4 612 ijij rr ji V σσ ε ?= ∑ < ]2[ )1()1(2 0 0 0 0 ?? < ?? ∑ ?= r ij r r ij r eeV ji ρρ ε 吸引势 冲击势 c r c a 阻止原子靠的过近）部分排斥冲击 ()( 原子间平衡距离 )(约束原子不散开吸引势部分 部分组成。势不同，基本上由这两其它原子势尽管数学形 原子内嵌势 (EAM) EAM,Johnson,Finnis-Sutton,Glue,Mei ? 嵌到电子气中原子所受 到的力，是原子所处位 置电子密度的函数 F(ρ i ) 是 吸引势。 ∑ ∑ ≠ ?= iji iij FrE )()( 2 1 ρφ ∑ ≠ = ij iji rf )(ρ ? 两体势 ? 一般情况下，带正 电的原子实间相互 作用，是 冲击势 EAM. S.W.Daw and M.I.Baskes, PRL50,1285(1983);PRB29,6443(1984) S.M.Foiles, M.I.Bakes and S.W.Daw, PRB33,7983(1986) Johnson. PRB37,3924(1988);PRB39,12554(1989). Finnis-Sutton. Philos. Mag.50,45(1984),Philos.Mag.Lett,61,139(1990) Glue. F.Ercolessi and J.B.Adams, EuroPhys.Lett.26,583(1994) Mei. J.Mei etal Phys.Rev.B43,4653(1991) e ∑ ∑ ≠ ?= iji iij FrE )()( 2 1 ρφ ⊕ ⊕ 紧束缚势 ∑∑∑ ???? ?= ij q ji p r ij r r ij r eeAE 2 1 00 }{ )1(2 2 , )1( αβ αβ αβ αβ αβαβ ξ 紧束缚势与 EAM势形式上完全一样，第一项代表离子 间的排斥项，第二项为电子气对原子的相互势，是吸 引项。再进行合金的模拟中，需要确定不同原子间的 紧束缚参数。 PRB48,22(1993) Tersoff&Brenner势 (I) ∑ < += ji ijAijijRijijC rfbrfarfE )]()()[( )exp( 2 rBf A λ??= 吸引势 )exp( 1 rAf R λ?= 冲击势 ? ? ? ? ? +> +<<?? ?< = ? DRr DRrDR DRr rf D Rr C ,0 ],sin[ ,1 )( )( 22 1 2 1 π Tersoff?PRB37,6991(1988),PRB39,5566(1989) Tersoff&Brenner势 (II) 三体效应包含在 a ij 和 b ij 中 n n ij n ij b 2 1 )1( ? += ζβ n n ij n ij a 2 1 )1( ? += ηα ∑ ≠ ?= jik ikijijkijCij rrgrf , 33 3 ])(exp[)()( λθ? i j k θ ijk ∑ ≠ ?= jik ikijikCij rrrf , 33 3 ])(exp[)( λη ])cos([ 22 2 2 2 1)( ijk hd c d c ijk g θ θ ?? ?+= Brenner C-H势与 Tersoff势基本相同， 情参阅 ,Brenner?PRB42,9458(1990) 分子动力学中的系综 1.微正则系综 2.等温分子动力学 3.等压分子动力学 4.非平衡态分子动力学 5.最速下降 (Steep Desend) Nose-Hover等温分子动力学 rsv & rr = )ln()1( 2 2 1 0 sTkfsQHH B +?+= & ∑ + ??= i s Tkf B srrmsQ )1( & r & r && s rs ms F r & r & r && r 2 2 ?= 引入了新的自由度s,定义与其相关的动能和势能。 与老系统的耦合体现在对速度的重新标定。 f为系统的自由度数，Q为Nose质量。 新的自由度的引入相当与引入了一个恒温热库，系统与 其耦合趋于和热库热平衡。系统的温度恒温热库相同而 保持恒温 Nose, Mol. Phys. 52, 255-68(1984) Anderson等压分子动力学 )1( 2 2 1 0 ??+= VPVQHH ex & )4( 3 2 3 1 ??= ss V V mV F & r && r & r )5( )( ?= ? Q PP ex V && h )2(?= shr rr )3(?= shv & rr ∑ ∑ > += iij ijij m pp V FrP i ii , 3 1 )( α αα αα 3 hV = 一个正立方体中的原子系统，立方体边长为h，s为原子相立 方体的坐标，现在把(V,s)均做为变量。和原先系统(r) 比自由度增加1.(1)中右边2,3项为与V相关的能量。可以 看出等压过程是通过调节立方体的体积与速度来实现的。 Q为Anderson质量，调节Q可以控制压力涨落的大小. Andeson, H.C. J.Chem.Phys. 72, 2384-93(1980) 有限体系的等压分子动力学 0=+≡ PPP ext ) ext P L P })({ 0 ijext rVPLL += 0 L ∑ ?= N i i m p rL i i })({ 2 0 2 φ 0)..( 2 3 1 =????= ∑ N i eiiii V VPrrvmP φ ) ∑ ??= N i iii V rvmP ).( 2 3 1 φ ext N i i V P PVr e ?=? ∑ . 3 为常压的系统也为常数则为常数显然 0 ,,. LPVr N i i∑ ? 是常数从而是常数则的齐次函数是若选取 PnVVrrVrVrV N i i n ,.),()( ∑ =?=λλ 常数则 ≡== ∑ e N i r i PPV ij ,)( 3 23 4 γ π = 不需要引入质量参数 D.Y.Sun and X.G.Gong, J.Phys.Condens.Matter 14 (2002)L487-L493 Parrinello-Rahman等压分子动力学 sTr rr ?=VPTQHH ex ?+= ∑ αβ αβ 2 2 1 0 & ),,( 321 tttT rrr = sGmGFHsm & & r && r 11 ?? ?= 321 tttV rrr ×?= T ex HVIPPTQ )()( 1? ?= && TTG T = 与 Anderson方法不同之处在于 ,不但能改变 元胞的大小又可以改变元胞的形状，附加 的自由度是 6个 Parrinello and Rahman Phys.Rev. Lett. 45, 1196-99(1980); J. appl. Phys. 52, 7182-90(1981) 非平衡态分子动力学 (NEMD) )1()(}){},({ ?Γ?Α+= tprHH ii ne r rr )2()( ?Γ?+= tAr r m p r & r r )3()( ?Γ??= tAFp p rr & r AA rr ?= AA PP ?= 在外场的作用下，系统内由粒子输运和热输运。系统 一般处于非平衡态下。可用含时哈密顿量来表示 (1)。 运动方程由式 (2)(3)表示。 A是 3Nx3N阶矩阵， Γ是 3N维 的矢量 . M.P.Allen and D.J.Tildesley, Computer Simulation of Liquids, Oxford University Press(1987) Langevin分子动力学 )()( tGrmrVrm iiiii +???= & r && r γ Gauess白噪声 )()()( '' ttCtGtG ?>=< δ 通过调节噪声强度 C来控制系统的温度 K.W.Jacabsen, et.al Phys.Rev. B35, 7423(1987) 最速下降 (Steep Descend) 1 v r 0 1 =v r 2 2 1 112 ttvrr m F δδ r rrr ++= MD 0 1 =v r 0 2 =v r 0 4 =v r 稳定构形 F r 0 3 =v r 2 2 1 12 trr m F δ r rr += Steep Descend VF ??= r 由于力是势能的梯度，方向垂直与势能面，并且沿着能量减小的方向， 因此沿着力的方向运动，势能总是减小，最后趋向一个极小点 物理量的计算 A.扩散系数的计算 ∑∑∑ === ?+=?+= M m N i mimi NM s N i ii N sRstRsRstRtR 11 2 11 1 2 1 2 )]()([)]()([)( rrrr t tR t D 6 )( 2 lim ?? ∞→ = .代表对初试点的平均 s 为了使均方位移与初始点 (s)的选取无关，通常选取 多个初始点 (M个 )，然后再取平均。 ,6)(, 2 DttR =子系统对于处于平衡态的多原 α DttR =)( 2 处于非平衡态的系统 )(1 ),(sup,1 onsubdiffusi nerdiffusio 是亚扩散 是超扩散 < > α α 也可以称为扩散系数慢的标志 做为扩散快但我们仍然可以用 , ,1 D≠α 分布距离是飞行的特征，其跳跃的 其扩散轨迹具有对于团簇的扩散 Levy Levy, 1＜α 1＝α 1＞α 3',1 )(,)( ' << == ?? μμ μμ CttPCllP taillongLevy ?分布 玻尔兹曼分布 B.速度关联函数 ∑∑ = ∑ ∑ + = = = == N i svsv svstv N N i vv vtv N M m mimi M m mimi ii ii tC 1 )()( )()( 1 1 )0()0( )0()( 1 1 1 )( ∑∑ = ∑ ∑ + = = = == N i svsv svstv N N i vv vtv N M m mimi M m mimi ii ii tC 1 )()( )()( 1 1 )0()0( )0()( 1 1 ,, 1 ,, ,, ,, )( βα βα βα βα αβ 速度关联函数是相同原子不同时间的关联 ,N是你所感兴趣的原子数， 可以是一个原子，或多个原子。我们也可以研究，速度 不同分量间的关联。 α=x, β=y ∫ = dtttCC P )cos()()( 2 ωω π 声子普 (声子态密度） ∑ ≠ ??= ji ji rrrrg |)|()( 0 δρ 1)(lim = ∞→ rg r 0)(lim 0 = → rg r 0)(lim = ∞→ rg r 0)(lim 0 = → rg r )(4 2 rgrRCF π= )(rg r 1)( =rg )(rg r C.对相关函数 无限系统 有限系统 径向分布函数 RCF是实验可以直接测量的量 (中子衍射， X射线衍射） D.结构因子 s(k) ∑ ??? = ji rrki N ji ekS , )( 1 2 )( rr r r )(kS r 1 有序程度 需要选择一系列倒格失 k进行计算或选择某一特殊的 k 0 3 2 1 ............. 3 3 2 2 1 1 ρ=→ → → → ∞ ∞ ∞ V N V N V N V N r r r r ? ? ? ? ? ? ? ? 11 Vr ?→ 22 Vr ?→ 33 Vr ?→ 1 N? 2 N? 3 N? ∑ ??= ji ji N rrtrtrG , 1 )|)0()((|),( 2 δ E.Van Hove 函数 相同原子间的关联，不同原子间的关联，时间演化效应 ∑ ??= i ii N rrtrtrG )|)0()((|),( 2 1 δ Van Hove自相关函数 ∑ ≠ ??= ji ji N rrtrtrG )|)0()((|),( 2 1 δ 随时间演化的对关联函数 )(rg 对时间求和得到对关联函数 ∫ = c r kr kr drtrGtkF 0 )sin( ),(),( 中子非弹性散射 The Intermediate Scattering Function 角分布函数 ∑∑ =< ? ??Θ?Θ= N kji rr rr kjki N kjki kjki c rrrrrA 1 1 )(cos)()(),(cos rr θδθ ∑∑ =< ?Θ?Θ= N kji kjkiC rrrrN 1 )()( 2 1 ),(cos rA θ ),(cos rA θ 0 2 1 θcos θcos 2 1 ? 1?0 分布函数 P(Q) ∫ ∑ == = dQQQPQQ step N i i T )( 1 1 ∫ =1)( dQQP 如果分子动力学运行 N step 步，那么它在相空间采样 N step 个点 ,其中 P(Q)是这 Nstep个采样点在区间 Q+dQ中出现的概率。其中 Q可以是 能量 或其它你感兴趣 的物理量。可以从比较少的采样中获得系统的有用信息 ∑ = =? N i ii vmK 3 1 2 2 1 动能 )(EP B Nk K T 3 2 =?温度 ∑∑ >= ?+=? ji ijij N i ii V FrvmP )( 1 2 3 1 r rr 压力是 E 1 E 2 E 一个假想系统的能量分布函 数，能量在 E 1 ,E 2 附近有两个 峰值，预示系统有两个理想 的状态 2 2 2 )( 2 3 2 3 )1( TkN KK C Nk B B ? =?? υ 比热 分 子 动 力 学 的 程 序 结 构 Verlet列表技术 f r skin r 1 d 2 d 0 t 0 t ?+ 0 t ?+ 0 t ).(, 2121 dddd >中两个最大的是对初始时刻位置偏离 有变化。 本身没列表重新计算，但列表的变化大于 离原子相对于中心原子距没有变化。例如图中的 列表有可能尽管列表就全都重新计算， 整个对原子该条件被破坏，序设计中，只要任何一 在实际的程列表不需要重新计算。都小于 间的距离变化也就是说任何两个原子 。如果变化不会大于任何两个原子间的距离 Verletr c VerletVerlet Verltr rdd dd skin skin skin , , , 21 21 <+ + c c c r )(rV r 算量大大增加。来的复杂性有可能使计会更加复杂性，由次带 计算列表的条件算列表的次数，但判断我们可以尽一步减小计 周期性边界条件 )int( )int( )int( boz z ijij boy y ijij box x ijij ij ij ij anbozzz anboyyy anboxxx ×?= ×?= ×?= )|(|)( )|(|)( )|(|)( 2 1 2 1 2 1 ??= ??= ??= ijijijij ijijijij ijijijij zzsigzz yysigyy xxsigxx θ θ θ 是质心 ccc ijijij boz zz ij boy yy ij box xx ij zyx zyx z y x cij cij cij ,, ),,( 2 1 2 1 ?≥≥ ? ? ? ? ? ? 2 3 2 1 2 1 || 1 1 ≤ +=?< ?=> x xxx xxx 要求 则如果 则如果 期性边界条件。 耗时比较少的周 0 2 1 2 1 ? box boy θ t tt δ+ Cii N i i N C Rrr rR r rr r r ?= = ∑ = 1 1 1 Cii N i i N C vvv vv rrr rr ?= = ∑ = 1 1 1 JI vrmJ N i iii r r rr r ?= ×= ? = ∑ 1 1 )( ω trrrr rvvv iiii iiii δω ω ?×+=? ×+=? )( 112 112 rrrrr rrrr 质心移至原点 相对于质心速度 整体转动角速度 以后的位置和速度 扣除平动和转动 与以削除。 由度的影响可暂时体平动和转动等外部自 自由度感兴趣，整质时，我们只对其内部 究团簇等体系的性统的平动和转动。在研 差也有可能引起系分运动方程所带来的误 动和平动。另外积统始终会有一个整体转 ，那么系速度不为果初始状态的角速度和 角动量均守恒，如在孤立系统中，动量、 0 ∑ ?+++= N i iiiiii xxxxxmI ])1()[( 2 3 2 2 2 1 βα δ αβαβ αβ δ ? ? ? ? ? ? ? ? ? ? ?? ?? ?? 333231 232221 131211 III III III 扣除整体平动与转动