分子动力学与原子多体势.docx

分子动力学与原子多体势.docx

《分子动力学与原子多体势.docx》由会员分享，可在线阅读，更多相关《分子动力学与原子多体势.docx（19页珍藏版）》请在冰点文库上搜索。

分子动力学与原子多体势

一．分子动力学简介

随着纳米科技的到来，许多新的学科产生了，例如纳米电子学、纳米生物学、纳米材料学、纳米机械学等。

人们的注意力逐渐从宏观物体转向小尺度及相应的器件，其中微机械系统（mieromachine）或称微型机电系统（mieroe一eetro一meeh耐ealsystem，MEMs）尤其取得了成功，并正被拓展应用于各种工业过程。

由图1可知，分子动力学正是处于nm尺度下的研究方法。

图1.不同模拟方法所对应的空间和时间尺度

1957年Alder和Wainwright[1]开创了分子动力学（Moleeularnynamies，MD）方法，之后经过多位科学家的努力，拓展了分子动力学方法的理论、技术及应用领域，尤其是在20世纪80年代由Andersen等[2]先后完成的恒温、定压分子动力学方法，标志着分子动力学方法的科研应用进入了一个新阶段。

分子动力学方法是研究纳米尺度物理现象的重要手段。

随着越来越多的材料原子间作用势函数被精确描述并经过实验验证、计算机硬件水平的快速更新以及高效率新算法的提出，分子动力学模拟被广泛应用于纳米尺度力学行为和纳米材料力学性能的研究。

在纳米尺度下，材料由离散的原子排列而成，由于比表面积大、表面效应明显，材料的力学性能和力学行为将与宏观材料迥异。

基于连续性假设的宏观连续介质理论在研究材料的损伤演化、失效过程时，往往在时间和空间上将原子尺度的缺陷进行平均化处理，但这种处理仅适用于大量缺陷分布在材料中计算区域的情形，而对许多细微观材料和力学实验观测到的现象都无法解释，如疲劳与蠕变过程中的位错模式、塑性变形的不均匀性、脆性断裂的统计本质、尺寸效应等。

因此，连续介质理论显然难以准确求解纳米尺度的力学问题。

同时，如果直接从第一原理出发进行计算，除了类氢原子以外其他材料的薛定愕方程求解难度都太大，而且局域密度泛函近似理论并不是总能满足实际问题的需要。

另一方面，材料本身在空间、时间和能量等方面存在藕合和脱祸现象[3,4]，直接从头开始的量子力学计算难以很好地应用到几百个原子以下的计算规模中，无法达到一般纳米材料和器件的模拟要求。

此外，由于实验条件控制的困难和合成、制备方式不同，各种纳米材料力学性能的有关实验结果分散性较大甚至相反，以至于目前难以通过纳米力学实验得到普遍适用的定量力学规律。

鉴于理论和实验上的困难，通过分子动力学方法模拟纳米尺度的力学性能和行为来探索纳米尺度的一般规律，是进行纳米力学研究的有效方法。

分子动力学最早用于热动力学和物理化学，计算不同物理系统如固体、气体、液体的整体或平均热化学性能。

1957年Alder首次提出并采用分子动力学方法分析刚性球系统的固液相变问题取得成功，此后，分子动力学开始逐渐应用于材料领域。

随着上世纪80年代计算机硬件水平的提高和各种描述原子间作用的势函数的提出，分子动力学模拟日益活跃。

通过分子动力学模拟不仅可以深入了解复杂的机制，发现本质上崭新的现象，而且可以定量模拟真实固体中所发生的过程，是诸如表面结构和扩散中的动力学和稳定性等许多结果的唯一来源[5-9]。

在EAM理论逐渐成熟和Baskes实验室[10-13]、Ackland实验室[14-18]等精确测出大量常用材料的EAM参数以后，分子动力学方法在模拟材料的物理性能和现象方面逐渐显示出强大的计算能力和较高的精度，大量的模拟尤其是固体结构、位错运动、表面界面现象、力学性能、变形机制和流体中的电泳、电渗透流等都得到了理想的结果。

只要能将基于物理的模型建立起来，通过分子动力学计算就可以揭示出物理现象的本质，逐渐被广泛应用于凝聚态物理、材料学、力学、生物学、微电子学和微纳米加工等领域。

晶体是由大量的原子有序排列而成，材料的强度来源于原子间的相互作用，塑性来源于原子间的相互运动。

因此，直接从原子尺度上对材料的微观力学行为进行研究显得非常有必要。

分子动力学模拟技术既能得到原子的运动轨迹，还能像做实验一样观察。

对于平衡系统，可以在一个分子动力学观察时间（ObservationTime）内做时间平均来计算一个物理量的统计平均值，对于一个非平衡系统过程，只要发生在一个分子动力学观察时间内（一般为1一10ps）的物理现象也可以用分子动力学计算进行直接模拟。

特别是许多与原子有关的微观细节，在实际实验中无法获得，而在计算机模拟中可以方便的得到。

这种优点使得分子动力学方法广泛运用于材料科学与工程中，如材料设计、断裂分析等。

近年来，利用计算机模拟技术研究材料的力学性能日益成为人们感兴趣的课题。

由于计算机处理速度的迅速提高，计算机模拟已经和实验观察、理论分析并列成为本世纪科学研究的三种方法[19]。

计算机模拟的数据（从模型中得来）可以用来比较、验证各种近似理论;同时，计算机模拟方法还可用来对实验和模型进行比较，从而提供了评估一个模型正确与否的手段。

计算机模拟方法还有一个优点就是可以沟通理论和实验。

某些量或行为可能是无法或难以在实验中测量的，而用计算机模拟方法，这些量可以被精确的计算出来[20]。

分子动力学模拟方法更以其建模简单、模拟结果准确的特征而倍受研究者的关注。

分子动力学模拟（MolecularDynamicsSimulation）用于计算以固体、液体、气体为模型的单个分子的运动情况，是一种联系微观世界与宏观世界的强有力的计算机模拟方法侧。

二．分子动力学基本原理

分子动力学方法是纳米计算力学的主要手段。

它对多经典粒子系统的运动方程组进行数值积分，得到相空间轨道，进而研究系统的平衡热力学性质、结构动力学性质和非平衡传输性质等。

在分子动力学模拟中，将原子视为质点或忽略内部结构的固体球。

首先建立粒子系统的几何模型，通过描述原子间相互作用的势函数，根据给定的边界条件、初始条件求出系统中每一时刻单个粒子或原子的能量和所受到的力，代入牛顿动力学方程组求解原子的位置和速度，得到系统在相空间的运动轨迹。

对足够大数量的粒子在足够长时间的结果进行统平均，则可以得到类似于宏观意义上的物理量和力学量。

早期研究孤立系统的保守系综分子动力学基于两个基本假设[21]:

a.粒子是相互作用的质点，运动由位置矢量和速度矢量来描述。

粒子间的相互作用取决于彼此的空间位置。

b.系统无质量交换。

即模拟过程中系统的原子数不变。

保守系综分子动力学常用于模拟能量守恒的孤立绝热系统。

但是要准确模拟实际纳米材料和器件的表面、界面等现象，真实反映粒子系统受外界约束情况下的物理行为，以及进行跨尺度的祸合模拟，保守系综分子动力学模拟明显是远远不够的。

在此基础上，研究者提出了耗散系统与周围介质进行能量交换的不同理论和算法，对粒子系统进行温度、压力、粒子数、体积等不同控制，适合于不同系综的分子动力学模拟。

分子动力学模拟中假设系统所有粒子的运动遵循经典牛顿运动定律，且忽略原子背景电子云的量子效应，原子间的相互作用满足叠加原理。

可见，分子动力学方法是一种对广义牛顿运动方程的近似的数值积分方法。

1.下面分别描述分子动力学方法的基本方程。

（1）拉格朗日（Lagrangian）运动方程

对于N自由度的由互相作用的质点构成的运动系统，拉格朗日方程为：

式中q为广义坐标，指定质点的空间位置。

为质点位置对时间的导数。

在笛

卡儿坐标系下，由N个原子构成的模拟系统的拉格朗日方程可写为：

式中

，为原子i的位置矢量，在三维空间n=N/3。

对于互不影响的粒子系统，可取拉格朗日量：

式中mi，为粒子i的质量，L即为系统的动能，即系统所有原子的动能之和。

如果考虑到原子间的互相作用，L可改一记为

式右端的两项分别表示系统的动能和势能。

代入拉格朗日方程可得系统牛顿运动方程:

式中Fi即为原了i所受的内力，即由于系统中其他原子的作用而在原子i上体现出的合力。

由牛顿运动方程建立线性微分方程组，给定初始条件（初始位置、初始速度），求解该封闭方程，可得到确定解，求出原子运动的轨迹，即任意时刻原子的位置ri（t）和速度

i（t）。

分子动力学模拟中的经典拉格朗日方程常用来计算原子的运动，得到单个原子的运动轨迹，描述原子系统的运动过程，以及反映在原子系统整体特征下的位移、变形、缺陷等。

（2）哈密顿（Hamiltonian）运动方程：

如果采用广义坐标系和动量的形式来描述粒子系统的运动，则可以得到哈密顿形式的运动方程，求出多粒子系统的状态和演化过程。

在笛卡儿坐标系中，对含n个粒子的保守系统，取系统哈密顿量为：

式中Pi，为粒子i的动量，Pi=mi

i，哈密顿正则方程为:

如果给定系统的初始状态（初始位置和速度），则可以求解线性微分方程组，得到系统原子的运动轨迹，即任意时刻的原子位置和动量，由统计平均得到系统的热力学表征。

在分子动力学模拟的基本方程中，拉格朗日运动方程更适合于求解原子系统运动的过程（求解原子速度和位置），并能施加外部荷载如外力、约束、边界条件等，而保守统的哈密顿正则方程更适合于求解系统的动力学演化过程和热力学状态，如温度、热流动。

在实际的分子动力学模拟中，拉格朗日方程通常与外部约束和边界条件一起，构成特定原子系统模型的初值问题;而Hamiltonian运动方程往往在保守系统的基础上，通过改变系统状态变量，与外部环境进行能量交换，构成不同系综的分子动力学模拟。

2.分子动力学原子间作用势函数

分子动力学方法是通过原子间的相互作用势，按照经典牛顿运动定律求出原子运动轨迹及其演化过程。

分子动力学计算的关键是原子间势函数的选取，它决定着计算的工作量以及计算模型与真实系统的近似程度，直接影响到模拟结果的成功与否。

由于物质系统的复杂性以及原子间相互作用类型的不同，很难得到满足各种不同体系和物质的一般性而又精度较高的势函数。

所以针对不同的物质体系人们陆续发展了大量的经验和半经验的势函数[22]。

从分子动力学模拟的基本方程可以看出，分子动力学模拟即求解偏微分方程，求解的精度关键在于势函数U（r1,r2,…rn）描述原子间相互作用的准确性。

多原子系统的势函数可以表示成:

式中Um为m体势，即原子势能受m阶多体效应的影响，其中U1为系统原子受外力场（如重力）的影响项。

在通常的计算中，由于系统内原子受原子间相互作用的影响远大于单个原子本身受外力场的影响，同时为了减少计算量，一般忽略外力场的影响U1和三阶以上的多体效应Um（m>3），这就是应用于分子动力学模拟的对势模型（PairPotential）。

对不同的固体材料，将其原子三阶以上多体效应的影响作为修正项记入对势模型的二阶作用项中，就形成了各种不同的多体势函数。

（1）对偶式

对偶势理论认为任何两个原子结合键的强度不因为它周围其它原子的存在而受到影响。

1.Lennard一Jones（L一J）势函数

在对势模型的典型代表、应用广泛的Lennard一Jones模型[23-24]中，取原子间作用的势函数方程为:

式中σ—原子与零点势的距离;

ε—最小势能处的能量;

m,n—调整系数，一般（m,n）取值为（12,6）,（10,5）,C8,4）

方程右端的第一项描述原子间的排斥作用，第二项描述原子间的吸引作用。

|ry|=|ri-rj|，即原子间的距离。

当原了间距离rij=σ时，原子势能Ф（rij）=0;rij>σ时，Ф（rij）<0;rij<σ时，Ф（rij）>0。

ε表征原子间的吸引/排斥强度，对原子中的一个从平衡位置r0移到无穷远处所作的功，对于不同物质的原子,ε各不相同，通常为10-19J或ev量级。

能量最低时刻的r值，即r0，表示原子间的平衡位置，此时原子间作用力为零。

对上述方程求极值，可得r0=

1.22σ

将上式代入式

，可得Lconard一JoneS模型中原了间作用力为：

2.Morse势函数

1929年，Morse注意到双原子分子的振动谱的量子力学问题可用指数形式的势函数解析地解决，并发现计算结果与实验一致[22]。

于是他提出如下形式的势函数:

式中D—结合能系数;

α—势能曲线梯度系数;

r0—分子间作用力为零时的原子间距。

原子间作用力分别描述为:

同样，等式右边第一项描述原子核间的短程排斥作用，第二项描述原子核对电

子云的吸引作用。

3.Born一Mayer势函数

Bom和Mayer[25]估计碱金属离子之间的排斥项可用指数形式表示，于是提出如下形式的势函数:

Bom一Mayer排斥项等价于Morse势函数的第一项，第三项表示偶极子一偶极子

相互作用，第四项表示偶极子一四极子相互作用。

Aij反映了泡利不相容效应。

该势

函数考虑了离子之间的库仑力，所以它能很好地应用于碱金属卤化物以及碱土金属

卤化物等各种离子型化合物。

4.Jonson势函数

利用对偶势人们能很好地拟合FCC金属的弹性常数但用于BCC金属却误差较大。

后来，考虑了原子体积效应，Johnson在对势项的基础上提出了用于BCC金属的Johnson势函数[26]:

.

（2）多体势

20世纪80年代陆续发展出许多考虑多体相互作用的新的势函数（关于中心对称的）。

主要有Daw和Bakes在1983年提出的嵌入原子势，Norskov和Lang在1980年发展的有效介质理论，凝胶模型，Finnis一Sinclair多体势等。

其中前三种势函数的理论基础是一致的，为我们在这里主要介绍嵌入原子势（EAM）和F一S多体势。

在模拟固体和复杂的分子结构时，要考虑原子间的多体效应，即势函数方程中的高阶项。

但多体势项中

的引入将导致参数过多，一方面这些参数的求解和实验拟合比较困难，另一方面过多的参数将带来计算上的困难。

从上世纪八十年代开始，科学家开始研究在对势的基础上进行多体效应修正，即考虑三阶以上的多体作用，将其作为一个修正项引入到对势模型中，使势函数能更精确描述固体的原子间作用。

根据这一思路，研究者先后提出了一系列新的多体势函数，如能较精确地描述C，Si原子间作用的Tersoff势[27]，描述C和烃类化合物的Brenne:

势[28]、适于描述体心立方晶格的Finnis一Sinclair[29]:

势等等。

这些势函数的共同点是在对势模型的原子间吸引项上增加一个修正参数，即将势函数写为:

式中

（rij）和

（rij）分别是对势模型中原子距离rij=|rij|的原子i、j之间的排斥作用项和吸引作用项。

εij表征其他的邻近空间原子对当前原子对的多体效应。

1.镶嵌原子势（EAM）

镶嵌原子势（EAM）[10-13]由Daw和BaskeS于1984年率先提出。

它在对势模型的基础上，假设将每个原子镶嵌在由所有其他近邻原子产生的背景电了云中，由电子云密度和将原子“镶嵌”到该电子云中所需的能量来衡量原子所处背景电子云的特征和背景电子石对势能的影响。

由于在电子云密度的计算中采用球面平均，镶嵌原子法特别适用于电子云基本呈球对称分布的金属原子。

这种方法单独求解电子百对原子核的多体效应，而不是在原子间对势吸引项上采取参数修正的拟合方法，结果更为精确。

依据上述原理，在镶嵌原子法中，将系统势能写为：

式中右边第一项即为传统的对势项，，第二项为多体修正项。

（rij）表示原子核i、j之间的排斥作用：

zi（rij）和zj（rij）为两原子核中的有效电荷数。

由于受到屏蔽作用，有效电荷数随着原子间距离rij=|rij|的增大而衰减：

式中Zi为原子的价电子数，α为待定常数，由材料的弹性模量确定。

F（ρi）是把第i个原子核“镶嵌”到密度为ρi的背景电子云中所需的镶嵌能，可以写为：

F（ρi）=DρiInρi

D是由材料决定的参数，ρi即为除第i个原子以外其他所有原子在原子i处所产生的电子云密度的线性叠加：

ρi=

Фj（rij）是离原子i处距离为rij的j原子在原子i处产生的电子云密度。

在EAM模型提出并得到逐步完善后，特别是计算结果得到大量实验验证后，原子间作用的精确描述问题逐渐转移到测试和求解EAM势参数上来。

针对原始的镶嵌原子求解电子云密度时采用球面平均而不适用于极性原子键结构，并对部分材料预测不准，暴露出对堆垛层错计算过大的缺点，Baskes研究组对EAM进行了部分修正，并于上世纪90年代先后测出了元素周期表上常用元素的EAM参数，使得EAM模型得到更广泛和方便的应用。

Doyama等从第一原理出发给出了求解立方晶格金属原子镶嵌势参数的方法，并得到实验验证。

在他们的求解中，将势函数方程中的对势项取为：

电子云密度假设为：

以上两式中rc1和rc2分别为计算原子间作用项和电子云密度项的截断半径。

这样，立方晶格金属晶体势函数的求解就集中在求解由材料物理属性如结合能、玻恩稳定性、弹性常数、空位形成能和堆垛能等决定的5个参数A1、A2、c1、c2、D上。

2.Finnis一Sinelair势[25]

F一S势和EAM势的理论基础不同，它是建立在紧束缚（TB）近似的基础上的。

F一S势函数可以写成如下形式:

V和Φ均为经验拟合的对势函数，式（2一23）中第一项表示核一核之间的相互作用，第二项不再是嵌入能，而表示多体相互作用的关联能，Ф可以理解为二分量紧束缚近似下的积分平方求和。

Finnis一Sinclair模型最初用来描述具有BCC结构的纯金属的原子间相互作用，Ackland等人将它推广到二元合金体系。

对于二元合金体系，可以用下面6个对势函数表示VAA、VBB、VAB、ΦAA、ΦBB、ΦAB。

Ackland等认为VAA、VBB、ΦAA、ΦBB和纯金属的参数相同；函数ΦAB取ΦAA和ΦBB的几何平均，这和紧束缚近似下的积分求和的解释是一致的。

因此在已知纯金属势参数的情况下，仅有对势项屹。

需要通过拟合而得到。

Ackland给出上述模型的函数:

H（x）是阶跃函数。

与EAM相比，F一S多体势形式较为简单，势参数易于拟合。

此为计算量比较少，程序容易收敛，所以应用广泛。

3.分子动力学模拟模型

在给定原子间作用势函数的条件下，可以建立模拟系统的几何模型，根据模拟的实际物理情况施加初始条件和边界条件，选取时间步长，然后按照一定的动力学积分算法，通过求解牛顿运动方程组或系统Hamilton正则运动方程组，计算系统中原子在每一步的状态，包括每个原子的位置、速度和原子势，并由原子势计算每个原子上的作用力，根据原子力给出下一时间步的原子位置和速度。

在模拟过程中对原子系统实施相应的边界条件调节，实现系统和周围介质的能量交换。

对足够长时间的计算结果和足够大系统的原子计算结果进行统计平均，可以得到类似于宏观概念的系统物理量。

（1）几何模型的建立

几何模型的建立就是按照要模拟的真实晶体中原子排列给出原子系统构型。

根据晶体结构分析，绝大多数金属晶体均为面心立方晶格（（FaceCenteredCrystal,FCC，如Au,Ag,Cu,Ni,Al,γ-Fe等）、体心立方晶格（BodyCenteredCrystal,BCC，如Cr,Mo,W,Nb,V,α-Fe等）和密排六方晶格（HexagonalClosedPacked,HCP，如Mg,Cd,Zn,Be等）这三种典型的紧密结构，单个晶胞的原子排列结构如下图。

空间方向上依次排列多个单胞，规定相应的品格尺寸，就构成了要模拟的单晶纳米材料，同时可以根据模拟需要在模型中引入晶格缺陷，如表面、裂纹等。

如果要构造多晶材料，则需要按照一定的分布形式如Voronoi几何构造方法构筑多个单晶晶粒，形成晶界。

图.常见金属晶体结构

（2）分子动力学模拟中的有限差分算法

分子动力学模拟的主要目的是计算任意时刻相互作用的原子系统中各原子的运动轨迹，即求解动力学系统初值问题，理论上可以采用有限差分法求解系统的牛顿运动方程。

分子动力学模拟中，粒子的运动方程是一组常微分方程组，求解的基本原理都是利用有限差分法，以一定的时间步长对方程沿时间轴逐步积分，但分子动力学对计算时间有一些特殊要求。

首先由于计算分子间作用力所占的工作量极大，所以凡是要求在一个时间步中重复计算分子间作用力的算法是不适宜的。

第二方面的要求涉及到解的精度和稳定性，这要求算法的阶数至少应大于3。

适用于分子动力学模拟的算法有Verlet算法、Leap-frog算法、Gear算法、Euler算法等

1）Verlet算法

Verlet方法基于有限差分思想，由初始时刻的状态值及其时间导数得到下一时

刻的状态值，并由此逐步计算出整个物理过程。

在分子动力学计算中，对原子的状

态进行Taylor级数展开:

两式相加可得:

由中两式相减可得原子速度:

根据系统的牛顿运动方程可得原子加速度:

这就是Verlet算法的基本形式。

2）Gear算法

为了求解非线性常微分方程，Gear提出了基于预测一校正积分方法的Gear算法[30]。

Gear0预测一校正法通过Taylor级数展开，由当前时刻原子的位置、速度、加速度以及其他高阶导数计算下一时刻原子的位置、速度、加速度，即预报值，然后根据原子位置的预报值计算原子间作用力和实际加速度值，通过比较实际加速度值和预报加速度值，进行误差校正，得到下一时刻的原子位置、速度、加速度以及其他高阶导数。

将t+∆t时刻的位置、速度和加速度对时间t作泰勒展开可得预测步为：

由新给出的位置rp（t十∆t）可计算t+∆t时刻的力F（t+∆t）和加速度αr（t+∆t）,根据加速度αr评价预测的加速度αp（t+∆t）的误差为:

将此差值与预测项相加得修正项为:

式中C0、C1、C2和C3为修正项系数。

在分子动力学模拟中，除上述三种常用的算法外，还有很多其它一些用于特殊要求情况下的算法。

Swope提出的Velocity-Verlet算法可以同时给出位置、速度与加速度，并且不牺牲精度。

这种算法的优点是给出了显式速度项，计算量适中。

Beeman提出的Beeman算法运用了更精确的速度表达式，所以能更精确地计算系统动能。

然而，由于它的表达式很复杂，所以计算量很大[31]。

Satoh就各种分子动力学算法的稳定性所作的对比研究表明，最优越的算法是Velocity-Verlet算法，之后是Leap-frog算法，Beeman算法等。

Gear算法比这些算法都要差些，因其更复杂，且稳定性和能量波动性较差，因而只是常见于早期的分子动力学模拟以及简单系统的计算。

参考文献：

[1]AlderBJ,WainwrightTE.PhaseTransitionforHardSphereSystem.TheJournalof

ChemicalP坤sics.1957,27（5）:

1208-1209.

[2]StillingerFH,WeberTA.PhysRevB.1985,31:

52-62.

[3]冯瑞，金国均.凝聚态物理中的基本概念.物理学进展.2000,20

（1）:

1-21.

[4」胡壮麟，王鲁红.电子和原子层次材料行为的计算机模拟.材料研究学报.1998,

12

（1）:

1一19.

[5〕梁海戈.纳米铜力学行为的分子动力学模拟〔博士论文].合肥:

中国科技大

学.2001.

[6〕吴恒安，倪向贵，王宇等.金属纳米棒弯曲力学行为的分子动力学模拟.物理学

报.2002,51（7）:

1412-1415.

[7〕吴恒安，王秀喜，梁海戈等.金属纳米丝力学行为研究进展.金属学报.2002,

38（9）:

903-907.

[8〕吴恒安，王秀喜，倪向贵等.金属纳米杆弹性模量的直接原子计算.金属学

报.2002,38（11）:

1219-1222.

[9〕吴恒安.纳米尺度下结构和材料力学行为的分子动力学模拟