分子动力学概述-精品文档.docx

《分子动力学概述-精品文档.docx》由会员分享，可在线阅读，更多相关《分子动力学概述-精品文档.docx（23页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、分子动力学概述分子动力学分子动力学方法是一种计算机模拟实验方法，是研究凝聚态系统的有力工具。该技术不仅能够得到原子的运动轨迹，还能够观察到原子运动经过中各种微观细节。它是对理论计算和实验的有力补充。分子动力学总是假定原子的运动服从某种确定的描绘，这种描叙能够牛顿方程、拉格朗日方程或哈密顿方程所确定的描绘，也就是讲原子的运动和确定的轨迹联络在一起。在忽略核子的量子效应和Born-Oppenheimer绝热近似下，分子动力学的这一种假设是可行的1。所谓绝热近似也就是要求在分子动力学经过中的每一霎时电子都处于原子构造的基态。要进行分子动力学模拟就必须知道原子间的互相作用势。在分子动力学模拟中，我们一

2、般采用经历势来代替原子间的互相作用势，如Lennard-Jones势、Mores势、EAM原子嵌入势、F-S多体势。然而采用经历势必然丢失了局域电子构造之间存在的强相关作用信息，即不能得到原子动力学经过中的电子性质1。事实上，分子动力学就是模拟原子系统的趋衡经过。实际上，分子动力学方法就是确定某一描绘与初始条件、边值关系的数值解。我们假定系统经过M步长之后到达稳定，而这一稳定状态正是我们所求的。、分子动力学的算法分析首先，我们假定我们研究的系统服从Newton方程所确定的描绘，即：)(1)(.tFmtr=1式中r(t)表征原子在t时刻的位置矢量F(t)表征原子在t时刻所遭到的力，它与所有原子的

3、位置矢有关m表征原子的质量。假如我们给定初始条件，即方程的定解条件r(0)和v(0)，那么方程的解就能够确定。年代中期发展了大量的分子动力学算法，如两步差分算法2、预测-校正算法3、中心差分算法4、蛙跳算法5等等。为了方便导出它们，我们以Euler一步法6来讨论之。我们令)()(.trtv=表征粒子的速度，则有：)()()(1)()(.tvtrtFmtrtv=记?=?=)()(1)()()()(.tvtFmtftrtvtw则有)()(.tftw=？欧拉一步法就是用向前差商来替代一阶导数，即：)()()1(.twhkwkw=-+，其中h是时间步长，将之代入则有：)()()1(thfkwkw=-+

4、即：?=?-+-+)()(1)()1()()1(kvkFmhkrkrkvkv)()()1()(1)()1(khvkrkrkFmhkvkv+=+=+对于式,由于给定了r(0)和v(0),故r(k+1)和v(k+1)能够确定。根据上述思路我们能够运用欧拉梯形法6、龙格库塔法6、Adams预测校正法6推导多种分子动力学算法，如Verlet算法7、Rahmans预测与校正格式8、Rahmans和Verlet预测校正格式9等。其实从物理意义上考虑，Verlet算法的推导只须用到Taylor展开，其经过如下：)()(!31)(!21)()()(432hotrhtrhtrhtrhtr+=+7)()(!31)

5、(!21)()()(432hotrhtrhtrhtrhtr+-+-=-7*上式中h代表时间步长，将以上两式相加我们得到：)()()()(2)(42hotrhhtrtrhtr+-=+8hhtrhtrtv2)()()(-+=8*该方法的优点是只须二维变量而得到四阶精度，但是不能同时得到同一步的速度和位置。下面列出了一些格式仅供参考：Verlet格式8mkFkakakahkvkvkahkhvkrkr)()(2)()1()()1()(2)()()1(2=+=+=+92)()()()(rCvCaErP9预测：midpointmethodpredictor)(2)1()1(khvkrkr+-=+(10)校

6、正：thesecond-ordermouctonscorrector2)()1()()1(2)()1()()1(kvkvhkrkrkakahkvkv+=+=+(11)3nrCvCaErP)()()()(9预测格式：thethird-orderAdams-Bashforthperdictor12)2(5)1(16)(23()()1(-+-+=+kvkvkvhkrkr(12)校正格式：thethird-orderAdams-Mouctonscorrector12)1()(8)1(5)()1(12)1()(8)1(5)()1(-+=+-+=+kvkvkvhkrkrkakakahkvkv(13)4Ra

7、hman和Verlet预测校正格式9预测格式：-=+-+=+112)1()()()1(qpppkabhkhvkrkr(14)校正格式：-=-=+-+-+=+-+=+112112)2()()1()1()2()()()1(qppqpppkadhhkrkrkvpkachkhvkrkr(15)系数pppdcb,能够从Taylor展开中得到。对于q=3,4,5情形，其系数如下表表一Rahman和Verlet预测校正格式的系数、分子动力学中的势模型分子模拟为何研究原子势能？为了计算力F(t),一般的考虑是从原子之间的互相作用势interatomicinteractionpotential着手，然后根据Ha

8、rtree定理求互相作用力，即F(t)=-.进而将求力转化从求势。作用势的选取是计算模拟的一个重要环节。合理选择作用势的一般要求有：1具有足够的精度，在实验范围内，计算结果能较好的和实验现象相吻合，在实验范围外，能够提供合理的预测结果；2在保证计算精度的前提下尽量简单，使计算机处理起来简易且速度快。为了计算势，现我们将势在原子平衡位置附近展开10：tot(r1r2rn)=?+)3()2(!31!21vv(16)其中)2(v表示二位体势pair-potential,只与原子的相对位置有关。)3(v表示三位体势假如只考虑前两项，历史上已存在对势模型有如下几种：1、Lennard-Jones势11-

9、13两原子间的势模型为：(r)=)()(4612rr-(17)式中第一项为排挤项，第二项为吸引项。，为根据实验数据拟合的势参数，r为原子间的距离。2、Morse势14两原子间的势模型为：)(exp2)(2exp)(00rrArrABr-=18式中A,B,r0为根据实验数据确定的经历参数，r0为原子间的距离。3、Born和Mayer根据Vandewaals理论提出了一种对势模型。即在Coulomb势上附加一Tail15。其模型为：)exp()(02112222rrrbrezzrr-+=19其中z1,z2为核电荷数，r1,r2是Goldsmidt半径，b12，r0根据实验数据确定的可调参数。实际上

10、，多粒子系统的互相作用是互相关联的，在计算两原子间的互相作用势时应该考虑周围原子的位形，此时就不能简单的应用以上对势模型，如处理d,f轨道没有被充满的过渡金属原子原子，近年来也发展了多种势模型以处理此类情形。4、F-S多体势16-18其模型为：)(21121=-=niijijtotsRssEjv20式中第一项代表核-核排挤作用的能量，第二项代表聚合能的多体吸引部分。其中sj=)(ijijRssji下标Si,Sj代表原子的种类，Rij代表原子i与原子j之间的距离。v和是根据实验确定的经历短程势函数，一般采用下列立体仿样函数：)()()()()()(321613RssHRssssRssRssHRs

11、sssRssrrArravjijijijijijijijikkkkkkkk-=-=上式中，rk,Rk为选用的节点，分别代表V和的切断半径。系数a1a6,A1,A2通过各自平衡位置的晶格参数。H(x)是Heaviside函数。铜的各拟合参数如下表所示:表一ParametersofpotentialsforCopper.Forpotentialsofpureelements,ak,Akareinevandrk,Rkarethecorresponding 185、嵌入原子势EAM19,20嵌入原子势EmbeddedAtomMethod,简称EMA是基于电子密度泛函理论提出的，由于只考虑电子浓度的影响

12、，而无须提供电子的原子种类，因此十分合适合金体系的计算。其模型为：Etot=+ninijiniijFR11)()(21(21)式中(Rij)为相距Rij的原子i和j之间的核核的互相作用势；F(i)为嵌入函数，能够理解为将一个原子嵌入到系统中的其它原子在I位置处产生的局部电子云中所需要的能量；i为系统中所有其它原子在i位置处所产生的电子密度：I=nijijR)(.其中电子密度分布采用Thomas-Fermi屏蔽函数的形式。、用方法与陈氏反演法推导势模型方法与陈氏反演法的共同特征是从晶体的结合能曲线出发计算原子的互相作用势。?关于方法10理想晶体中的结合能能够表示从所有原子对势的迭加，即：=1)(

13、21)(21)(nnnRxswrxE22式中r是原子的位置矢量；)(x是相距为x的两原子之间的互相作用势；nw是一样近邻的原子个数；xsn是第n近邻原子的距离。我们定义nR算符知足)()(xswRRnnn=据此能够改写式2为：=nnrRrRxE)(21)()(21)(?)(2)(1xErR-=?+-=+=-=-=-22312211112111111)()1()(nnnnnnnnRRRRRRRRRR由于：)(1)()()(11111111srwrRxswrR=?=-所以：?+-=-)()(1)()(21221111sxsEwwsxEwxERxnnn经历表明方法的收敛速度很慢。?陈氏反演方法、分子

14、动力学中的基本问题我们在推导出以上分子动力学算法和势模型后，就来利用来研究原子团簇。考虑到算法的收敛与稳定，结果的精度要求，计算速度以及计算机的内存限制，我们常碰到下列问题：、初始条件问题为了启动程序，算法要求我们给出一定的初始条件一般是初始位置和或初始速度。当然对于不同的算法和系统会要求不同的初始条件。一般来讲，我们所模拟系统的初始条件是不知道的晶体构造中的初始条件为未知吗？。事实上准确知道系统的初始条件也无实际意义。由于我们所研究的系统的初始状态具有随机性。这种随机性也知足一定的统计规律。但较好的初始条件会大量节省系统趋于稳定的时间，也就是讲，体系在相空间运动的时间将大大缩短，如下图：初始

15、条件能够有下列几种选择：1、位置在正规格子上，速度随机；2、位置随机的偏离格子，速度为零；3、位置速度都随机。通常我们的做法是假定原子的初始位置位于晶格格点上该格点能够是面心构造、体心构造等各种构造，而初始速度从Bolzmann分布规律或麦克思韦速率分布规律中抽取。假设我们给定的初始条件不在我们所研究的相图上怎么办？能否能够随意给定原子构造？对于详细的问题我们怎么给出它的晶格构造，如碳纳米管，C60等等？、时间步长和模拟步数以及终止条件的选取模拟的时间步长和步数共同决定了我们所研究的系统在相空间的大小。它的选取没有固定的规则，都是根据经历确定。这里有一经历规则是：系统总能量的涨落不超过势能的百

16、分之几。注：在分子动力学模拟计算中，时间步长h的选取是很重要的。为了减小误差，时间步长h必须取小些，但假如太小，趋于平衡所需要的时间就越长。因而，为了更快的趋于平衡，h又应该获得大些，即让粒子在相空间的步伐大一点。一般时间选为10-1410-17s所谓终止条件即研究的系统到达什么样的状态后，我们能够终止程序的运算。一般来讲当系统持续给出一样的平均动能和平均势能，我们就讲我们研究的系统已到达稳定状态。假如在模拟经过中发现系统很难到达稳定状态，能够通过适当的表度来解决这一问题。即适当降低或升高系统的动能。如我们得到的系统温度为T*,而我们所求的系统温度为T，则标度因子为*TT=，我们能够通过给每一原子的速度乘以这一标度因子，进而加快趋衡速度。