分子动力学分子动力学方法是一种计算机模拟实验方法，是研究凝聚态系统的有力工具。

该技术不仅可以得到原子的运动轨迹，还可以观察到原子运动过程中各种微观细节。

它是对理论计算和实验的有力补充。

分子动力学总是假定原子的运动服从某种确定的描述，这种描叙可以牛顿方程、拉格朗日方程或哈密顿方程所确定的描述，也就是说原子的运动和确定的轨迹联系在一起。

在忽略核子的量子效应和Born-Oppenheimer绝热近似下，分子动力学的这一种假设是可行的[1]。

所谓绝热近似也就是要求在分子动力学过程中的每一瞬间电子都处于原子结构的基态。

要进行分子动力学模拟就必须知道原子间的相互作用势。

在分子动力学模拟中，我们一般采用经验势来代替原子间的相互作用势，如Lennard-Jones势、Mores势、EAM原子嵌入势、F-S多体势。

然而采用经验势必然丢失了局域电子结构之间存在的强相关作用信息，即不能得到原子动力学过程中的电子性质[1]。

事实上，分子动力学就是模拟原子系统的趋衡过程。

实际上，分子动力学方法就是确定某一描述与初始条件、边值关系的数值解。

我们假定系统经过M步长之后达到稳定，而这一稳定状态正是我们所求的。

１、分子动力学的算法分析首先，我们假定我们研究的系统服从 Newton 方程所确定的描述，即：)(1)(..t F mt r =（1） 式中r(t)表征原子在t 时刻的位置矢量F(t)表征原子在t 时刻所受到的力，它与所有原子的位置矢有关m 表征原子的质量。

如果我们给定初始条件，即方程（１）的定解条件r(0)和v(0)，那么方程（１）的解就可以确定。

６０年代中期发展了大量的分子动力学算法，如两步差分算法[2]、预测-校正算法[3]、中心差分算法[4]、蛙跳算法[5]等等。

为了方便导出它们，我们以Euler 一步法[6]来讨论之。

我们令)()(..t r t v =（表征粒子的速度），则有：)()()(1)()(....t v t r t F m t r t v === （２）记⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎦⎤⎢⎣⎡=)()(1)()()()(.t v t F m t f t r t v t w （３）则有)()(.t f t w = ？？？？？？ （４） 欧拉一步法就是用向前差商来替代一阶导数，即：)()()1(.t w hk w k w =-+，其中h 是时间步长，将之代入（４）则有：)()()1(t hf k w k w =-+ （５）即：⎥⎥⎦⎤⎢⎢⎣⎡=⎥⎦⎤⎢⎣⎡-+-+)()(1)()1()()1(k v k F m h k r k r k v k v )()()1()(1)()1(k hv k r k r k F mhk v k v +=++=+ （６） 对于（６）式,因为给定了r(0)和v(0),故r(k+1) 和v(k+1)可以确定。

根据上述思路我们可以运用欧拉梯形法[6]、龙格－库塔法[6]、Adams 预测校正法[6]推导多种分子动力学算法，如Verlet 算法[7]、Rahman’s 预测与校正格式[8]、Rahman’s 和Verlet 预测校正格式[9]等。

其实从物理意义上考虑，Verlet 算法的推导只须用到Taylor 展开，其过程如下：)()(!31)(!21)()()(432h o t rh t r h t rh t r h t r ++++=+ （7） )()(!31)(!21)()()(432h o t rh t r h t r h t r h t r +-+-=- （7*） 上式中h 代表时间步长，将以上两式相加我们得到：)()()()(2)(42h o t rh h t r t r h t r ++--=+ （8） hh t r h t r t v 2)()()(--+=（8*）该方法的优点是只须二维变量而得到四阶精度，但是不能同时得到同一步的速度和位置。

下面列出了一些格式仅供参考： １． Verlet 格式[8]mk F k a k a k a h k v k v k a h k hv k r k r )()(2)]()1([)()1()(2)()()1(2=+++=+++=+ （9） ２．2)]()()()[(r C v C a E r P [9]预测：［midpoint method predictor ］)(2)1()1(k hv k r k r +-=+ (10) 校正：[the second-order moucton’s corrector]2)]()1([)()1(2)]()1([)()1(k v k v h k r k r k a k a h k v k v +++=++++=+ (11)3．n r C v C a E r P )]()()()[([9]预测格式：[the third-order Adams-Bashforth perdictor]12)]2(5)1(16)(23[()()1(-+--+=+k v k v k v h k r k r (12)校正格式：[the third-order Adams-Moucton’scorrector]12)]1()(8)1(5[)()1(12)]1()(8)1(5[)()1(--+++=+--+++=+k v k v k v h k r k r k a k a k a h k v k v (13) 4．Rahman 和Verlet 预测校正格式[9]预测格式：∑-=+-++=+112)1()()()1(q p p p k a b hk hv k r k r (14)校正格式：∑∑-=-=+-+-+=++-++=+112112)2()()1()1()2()()()1(q p p q p pp k a d h h k r k r k v p k a chk hv k r k r (15)系数p p p d c b ,,可以从Taylor 展开中得到。

对于q=3,4,5情形，其系数如下表表一Rahman 和Verlet 预测校正格式的系数２、分子动力学中的势模型（分子模拟为何研究原子势能？？？）为了计算力F(t),一般的考虑是从原子之间的相互作用势Φ（interatomic interaction potential ）着手，然后根据Hartree 定理求相互作用力，即F(t)=-▽Φ.从而将求力转化从求势。

作用势的选取是计算模拟的一个重要环节。

合理选择作用势的一般要求有：（1）具有足够的精度，在实验范围内，计算结果能较好的和实验现象相吻合，在实验范围外，能够提供合理的预测结果；（2）在保证计算精度的前提下尽量简单，使计算机处理起来简易且速度快。

为了计算势，现我们将势在原子平衡位置附近展开[10]： Φt o t(r 1r 2…r n )=∑∑∙∙∙++)3()2(!31!21v v (16) 其中)2(v 表示二位体势（pair-potential ）,只与原子的相对位置有关。

)3(v 表示三位体势如果只考虑前两项，历史上已存在对势模型有如下几种： 1、 L ennard-Jones 势[11-13]两原子间的势模型为：Φ(r)=])()[(4612rrσσε-(17)式中第一项为排斥项，第二项为吸引项。

ε，σ为根据实验数据拟合的势参数， r 为原子间的距离。

2、 M orse 势[14] 两原子间的势模型为：)]}(exp[2)](2{exp[)(00r r A r r A B r -----=φ （18）式中A,B,r 0为根据实验数据确定的经验参数，r 0为原子间的距离。

3、 B orn 和Mayer 根据Van de waals 理论提出了一种对势模型。

即在Coulomb 势上附加一Tail [15]。

其模型为：)e x p ()(02112222rr r bre z z rr -++=φ （19）其中z 1 , z 2 为核电荷数， r 1 ,r 2 是Goldsmidt 半径， b 12，r 0根据实验数据确定的可调参数。

实际上，多粒子系统的相互作用是相互关联的，在计算两原子间的相互作用势时应该考虑周围原子的位形，此时就不能简单的应用以上对势模型，如处理d,f 轨道没有被充满的过渡金属原子原子，近年来也发展了多种势模型以处理此类情形。

4、F-S 多体势[16-18] 其模型为：])(21[121∑∑=≠-=ni ij ij tot s R s s E jv ρ （20）式中 第一项代表核-核排斥作用的能量，第二项代表聚合能的多体吸引部分。

其中 s jρ=)(ij ij R s s ji ∑≠φ 下标S i ,S j 代表原子的种类，R ij 代表原子i 与原子j 之间的距离。

v 和φ是根据实验确定的经验短程势函数，一般采用下列立体仿样函数：)()()()()()(321613R s s H R s s s s R s s R s s H R ss s s R s s r r A r r a v j i jijijij i jijijikkk k k kkk --=--=∑∑==φ上式中，r k ,R k 为选用的节点，分别代表V 和φ的切断半径。

系数a 1……a 6,A 1,A 2通过各自平衡位置的晶格参数。

H(x)是Heaviside 函数。

铜的各拟合参数如下表所示:表一Parameters of potentials for Copper.For potentials of pure elements,a k ,A k are in ev and r k ,R k are the corresponding[18]5、嵌入原子势（EAM ）[19,20]嵌入原子势（Embedded Atom Method,简称EMA ）是基于电子密度泛函理论提出的，因为只考虑电子浓度的影响，而无须提供电子的原子种类，因而特别适合合金体系的计算。

其模型为：E tot =∑∑∑=≠=+n i ni j i ni ij F R 11)]()([21φ (21)式中φ(R i j )为相距R i j 的原子i 和j 之间的核－核的相互作用势；F(ρi )为嵌入函数，可以理解为将一个原子嵌入到系统中的其它原子在I 位置处产生的局部电子云中所需要的能量；ρi 为系统中所有其它原子在i 位置处所产生的电子密度：ρI =∑≠ni j ij R )(ρ.其中电子密度分布ρ采用Thomas-Fermi 屏蔽函数的形式。

６、用ＣＧＥ方法与陈氏反演法推导势模型ＣＧＥ方法与陈氏反演法的共同特征是从晶体的结合能曲线出发计算原子的相互作用势。

◐ 关于ＣＧＥ方法[10]理想晶体中的结合能Ｅ（Ｘ）可以表示从所有原子对势的迭加，即：∑∑∞===1)(21)(21)(n n n R x s w r x E φφ （22）式中 r是原子的位置矢量；)(x φ是相距为x 的两原子之间的相互作用势；n w 是相同近邻的原子个数；x s n 是第n 近邻原子的距离。