分子动力学――势函数的种类与发展

1.1原子论模型材料的研究目标之一：确定宏观性能与微观结构之间的关系。现代材料研究从某种意义上来说就是对微结构的研究。微结构，是指横跨埃到米的空间尺度上所有热力学非平衡态晶格缺陷的集合。空间尺度为几个埃～几米。时间尺度为ps～几年。考虑空间的离散化程度，材料模拟分为两类即连续体模型和原子论模型。连续体模型是在考虑了唯象和经验本构方程及平衡性、相容性和守恒定律所附加的约束条件下，建立起通常宏观情况下描述材料响应特性的微分方程，并由此微分方程求出单个原子的平均性质。如果要获得微结构性质更为详细的预测信息，则连续体模型将代之为原子论模型。原子论模型可给出更好的空间分辨率，与连续体模型相比，原子论模型包含有较少的唯象假说。原子论模型的基础是总能模型，即所有分子与凝聚态的原子模拟都必须建立在总能模型的基础上，多原子系统总能量E是原子空间坐标{xi,yi,zi}的函数。有三个层次的总能模型：第一性原理方法（First-principlesmethods），半经验方法（Semi-EmpiricalMethods）和经验势（EmpiricalPotential）或分子力场（ForceField）。这三个层次的总能模型在现代计算机上原子模拟的空间/时间尺度见表1。从表中可以看出，第一性原理方法的时间空间跨度较大，可计算的原子数也有限，所以如果单纯从第一性原理出发，而对材料进行完全的量子力学处理，那么不仅会在计算方法上存在困难，而且更会难以获得全面而准确的计算结果。于是，针对于不同的物质体系，研究人员研究发展了大量的经验和半经验势函数。表1原子模拟的空间/时间尺度First-principlesmethodsSemi-EmpiricalMethodsEmpiricalPotential/ForceFieldNo.ofatomsHundredsthousandsMillionsLengthscales1-5nm10nm50-100nmTimescales100psns100ns1.2对势分子动力学模拟的核心问题是求解牛顿运动方程组，而关键问题是原子间作用势的确定。（这句话也应该放在分子动力学部分而非相互作用势部分）原子间相互作用势是所有原子模拟的基础，可用于描述原子（分子）间相互作用的函数。控制着原子间的相互作用行为，从根本上决定着材料的所有性质。一般可由两个或多个原子之间的相对位置来定量确定出相互作用势，其中可以包括一系列参数诸如电荷，离子极化率。局域电子密度等。原子间相互作用势的精确与否将直接影响着模拟结果的准确性，而计算机模拟所需要的计算机机时则取决于势函数的复杂程度。故原子间相互作势的构建对材料的计算机模拟起着至关重要的作用。早期的原子间作用势多数是一些纯经验拟合势，近年来人们通过对基本电子结构的理论计算，发展了一些半经验的“有效势”。1957年，Alder和Wainwright对硬球模型系统的液固相变过程采用经典分子动力学（MD）方法进行模拟研究，这是人们第一次从原子尺度来研究物质的宏观性质，是分子动力学模拟史上的一座里程碑。（这一块内容应写在分子动力学的发展部分而非原子间相互作用势）对势在早期的材料研究中发挥了极为重要的作用，并且至今仍然活跃在计算机模拟的许多领域。若不考虑外力的作用，某体系包含N个相互作用的粒子，则该体系的总能可表示为：,321nNE（1.1）其中n表示N体相互作用的能量总和。此外，如果该体系为N个无相互作用的粒子，则该体系的总能量可以表示为：1'NE（1.2）二者的总能量之差为'NNEE，则该N个相互作用粒子体系总的相互作用能为,32'nNNEE-（1.3）若将其表示为粒子的坐标函数，则有如下表示：Nrrr,,,21（1.4）jijirrU,22（1.5）kjikjirrrU,,33（1.6）⋮kjinjirrrU,,33（1.7）其中，U2、U3和Un分别表示为两体、三体和n体相互作用，且Φ是可测量的，表示系统的总势能。一般认为Φ的多体展开式收敛很快，因此Φ中的高阶项可以被忽略不计。然而在实际中，对于Φ的n体展开式的计算处理上存在较大的困难，此外，多体原子的互作用也不容易定义。因此，人们在处理多体展开式时通常在三体项之后便进行截断，而其中的截断项贡献则可以通过添加不同的参数在剩余项中进行近似。在早期，人们甚至将展开项中的三体项也忽略了，故仅由两体相互作用之和来表示总势能Φ，即用对势来近似。在早期的材料研究中，对势发挥了极为重要的作用，尤其被广泛的应用于计算机模拟的重多领域，并且取得了相当不错的结果。根据对系统的总能量的贡献，我们可以将对势分为两类：一类为总能量完全由对势函数所决定的，该类对势可以有效的描述vanderWaals相互作用占据主导地位的体系；而另一类则为对势函数仅仅用来描述恒定的材料平均密度下系统的能量随着原子构型的变化，此类对势适用于描述sp-价态金属。经典对势可以写成UrENjiijijtot121（1.8）其中Ω为原子平均体积或材料平均密度，第一类对势U=0,第二类对势U(Ω)为原子凝聚对系统总能量的贡献,是原子平均体积Ω或材料平均密度的函数,同时φij也可以是Ω的函数。原子间作用势与原子间距的关系如图1所示。在每一个基本计算步，原子间作用力可由目标原子在截断半径以内与其周围原子之间的相互作用势的导数求出，即drrdEF，图1蓝色曲线所示。截断半径（cutoffdistance）是指：小于模拟盒子的边长的一半以避免与同一粒子的两个镜像同时作用。图1原子间作用势与原子间距之间的关系。（图1和图2如何命名，可以考虑二者合二为一）1924年，Lennard-Jones发表了著名的负幂函数式的Lennard-Jones势函数的解析式。1929年，Morse发表了指数式的Morse势。1931年，Born和Mayer发表了描述离子晶体中两离子间相互作用的Born-Mayer势函数。1983年至1984年间，Daw和Baskes提出了金属材料中的嵌入原子势(EAM)的概念和算法。几乎同时，Finnis和Sinclair根据密度函数二次矩理论提出形式上与EAM基本相同的经验F-S模型。接下来将简单介绍一下上述势函数形式。1.2.1Lennard-Jones势Lennard-Jones势也称L-J势，或6-12势，最初由JohnEdwardLennard-Jones于1924年提出，势函数采用负幂函数形式，用来描述两个中性原子（分子）之间的相互作用模型，是提出的最早的二体势模型。L-J势是一种典型的对势形式，其具体的表达式为：])()([000mrrnrrmmnBArnmnrr（1.9）其中，φ0表示相互作用的强度；r0表示原子的大小。通过量子力学的二次微扰论中偶极子-偶极子相互作用，可以推导出n=12，而该项则描述了范德瓦耳斯力。后一项表示排斥力，其一方面是来源于原子核之间的库仑斥力，而另一方面则是来源于由泡利不相容原理所产生的电子之间的交叠能。实际应用中，Lennard-Jones势常取作以下形式：])(2)[(601200rrrrr（1.10）有时还会增加一个rc参数，代表势函数的截断距离：当原子间距离超过rc时，认为原子间相互作用为0.从公式中可以很容易看出，L-J势由作用相反（正负号表示）的两部分组成：第一项（12次方项）代表排斥能，第二项（6次方项）代表吸引能。原子间相互作用势与原子间距的关系如图2所示。其中σ为相互作用的势能正好为零时的两体距离，ε为势能阱的深度。则公式（1.10）又可写为])()[(4612rrr（1.11）图2原子间相互作用势与原子间距的关系虽然式（1.10）、（1.11）为Lennard-Jones势最为常见的形式，但人们也经常采用其他形式的排斥项，尤其是指数形式ijijBrArexp。例如，采用指数排斥势和倒六方吸引项所构成的势，通常被称为Buckingham势。式(2.2)形式的Lennard-Jones势经常用来模拟惰性气体的凝聚态。LJ势函数模型很简单，计算量很小，适用于中性原子，尤其是惰性气体原子的描述。也有文献将该LJ势与其他势函数形式结合起来。1.2.2Morse势1929年，物理学家PhilipM.Morse注意到可以采用指数形式的势函数解析式来解决双原子分子振动谱的量子力学问题，并且计算结果与实验值一致。于是他提出了下面的势函数形式，即我们所谓的Morse势：]2[)()(200rrarraeeDr（1.12）或者DeDrrra2)(]1[0（1.13）其中，a和A是关于相互作用距离的维数和能量的常数。r0是两原子之间的平衡距离，r是两原子间距离。式中的参数D、r0和a可以通过对晶体的结合能、平衡点阵常数和体弹模量进行拟合得到。由于φ(r0)=-D，D是Morse势的阱深（势能零点可任意选取，在此将解离极限设为势能零点，即，两核间距趋于无穷远时令体系势能为零）；a则控制了势阱的“宽度”，a越小，势阱越宽。阱深D减去零点能就得到了解离能，在此E(0)为零，解离能为D。对Morse势在r0附近做Taylor展开，得到)(2100rrkr（1.14）其中二阶项中的k为平衡位置处的力常数。由此推导出，a，D和k有如下关系Dka2（1.15）而当a=6时，Morse势与Lennard-Jones势（6-12势）则极为接近，陈难先等人利用晶格反演的方法获得了若干金属元素的Morse势。Morse势常常用来构造各种多体势的对势部分，Morse势是一种针对于双原子分子间势能的简易模型，由于双原子分子势中隐含地包括了键断裂的现象。因此，Morse势对于分子振动微细结构能够得到很好的近似。1.2.3Born-Mayer势Born-Mayer势是由Born和Mayer于1932年提出的一种新型短程势形式，其目的主要是为了描述离子晶体中离子间闭壳层电子所产生的排斥作用，而其中的排斥项是来源于电子云和VanderWaals的相互作用，其具体的形式为：BrAer（1.16）其中，式中的参数A、B是通过平衡态的结合能和晶格常数来确定，主要用来处理离子晶体的模拟计算。然而有一种Born-Mayer-Huggins势模型的具体表达形式为：nijijijrCreAr)(（1.17）其中，式中的参数A、C和n是通过计算或实验值的拟合进行确定的。考虑到两离子间的相互作用可分为静电作用、短程排斥作用以VanderWaals作用。Fumi和Tosi认为两个离子之间的相互作用势是由Coulombic静电作用项ijjireZZ2和Born-Mayer形式的短程排斥项)(ijijreA以及VanderWaals吸引项(仅考虑偶极子—偶极子和偶极子—四极子吸引势)86ijijrDrC-这三部分组成的可叠加的有效两体势。有效的两体势（即Fumi-Tosi势）可写为862)(ijijijijijjitotrDrCreAreZZE-（1.18）其中，rij为i离子和j离子间距，e为单位电子电荷，Zi，Zj为离子带电量，在这里为±1,Αij为i离子和j离子间的软化参数(softnessparameter)，Aij反映了泡利不相容效应。C和D为VanderWaals系数。若(1.6)式中Aij的范围在无限大，则它的排斥作用会根据距离的增加而快速衰减，所以排斥力是短程的。该势考虑了离子之间的库仑力，所以它能很好地应用于碱金属卤化物以及碱土金属卤化物等各种离子型化合物。1.2.4Johnson势利用对势模型人们可以很好地拟合fcc金属的弹性常数，却不适用于bcc金属，此时，C12=C44,而实际上Cauchy压力P=(C12-C44)/2却往往不等于零.后来，Johnson在对势项的基础上添加一