分子动力学与原子多体势

本文格式为Word版，下载可任意编辑——分子动力学与原子多体势

一．分子动力学简介

随着纳米科技的到来，大量新的学科产生了，例如纳米电子学、纳米生物学、纳米材料学、纳米机械学等。人们的注意力逐渐从宏观物体转向小尺度及相应的器件，其中微机械系统(mieromachine)或称微型机电系统(mieroe一eetro一meeh耐ealsystem，MEMs)特别取得了成功，并正被拓展应用于各种工业过程。由图1可知，分子动力学正是处于nm尺度下的研究方法。

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

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

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

在纳米尺度下，材料由离散的原子排列而成，由于比表面积大、表面效应明显，材料的力学性能和力学行为将与宏观材料迥异。基于连续性假设的宏观连续介质理论在研究材料的损伤演化、失效过程时，往往在时间和空间上将原子尺度的缺陷进行平均化处理，但这种处理仅适用于大量缺陷分布在材料中计算区域的情形，而对大量微弱观材料和力学试验观测到的现象都无法解释，如疲乏与蠕变过程中的位错模式、塑性变形的不均匀性、脆性断裂的统计本质、尺寸效应等。因此，连续介质理论显然难以确凿求解纳米尺度的力学问题。同时，假使直接从第一原理出发进行计算，除了类氢原子以外其他材料的薛定愕方程求解难度都太大，而且局域密度泛函近似理论并不是总能满足实际问题的需要。另一方面，材料本身在空间、时间和能量等方面存在藕合和脱祸现象[3,4]，直接从头开始的量子力学计算难以很好地应用到几百个原子以下的计算规模中，无法达到一般纳米材料和器件的模拟要求。此外，由于试验条件控制的困难和合成、制备方式不同，各种纳米材料力学性能的有关试验结果分散性较大甚至相反，以至于目前难以通过纳米力学试验得到普遍适用的定量力学规律。

鉴于理论和试验上的困难，通过分子动力学方法模拟纳米尺度的力学性能和行为来摸索纳米尺度的一般规律，是进行纳米力学研究的有效方法。分子动力学最早用于热动力学和物理化学，计算不同物理系统如固体、气体、液体的整体或平均热化学性能。1957年Alder首次提出并采用分子动力学方法分析刚性球系统的固液相变问题取得成功，此后，分子动力

1

学开始逐渐应用于材料领域。随着上世纪80年代计算机硬件水平的提高和各种描述原子间作用的势函数的提出，分子动力学模拟日益活跃。通过分子动力学模拟不仅可以深入了解繁杂的机制，发现本质上崭新的现象，而且可以定量模拟真实固体中所发生的过程，是诸如表面结构和扩散中的动力学和稳定性等大量结果的唯一来源[5-9]。在EAM理论逐渐成熟和Baskes试验室[10-13]、Ackland试验室[14-18]等确切测出大量常用材料的EAM参数以后，分子动力学方法在模拟材料的物理性能和现象方面逐渐显示出强大的计算能力和较高的精度，大量的模拟特别是固体结构、位错运动、表面界面现象、力学性能、变形机制和流体中的电泳、电渗透流等都得到了理想的结果。只要能将基于物理的模型建立起来，通过分子动力学计算就可以透露出物理现象的本质，逐渐被广泛应用于凝聚态物理、材料学、力学、生物学、微电子学和微纳米加工等领域。

晶体是由大量的原子有序排列而成，材料的强度来源于原子间的相互作用，塑性来源于原子间的相互运动。因此，直接从原子尺度上对材料的微观力学行为进行研究显得十分有必要。分子动力学模拟技术既能得到原子的运动轨迹，还能像做试验一样观测。对于平衡系统，可以在一个分子动力学观测时间(ObservationTime)内做时间平均来计算一个物理量的统计平均值，对于一个非平衡系统过程，只要发生在一个分子动力学观测时间内(一般为1一10ps)的物理现象也可以用分子动力学计算进行直接模拟。特别是大量与原子有关的微观细节，在实际试验中无法获得，而在计算机模拟中可以便利的得到。这种优点使得分子动力学方法广泛运用于材料科学与工程中，如材料设计、断裂分析等。

近年来，利用计算机模拟技术研究材料的力学性能日益成为人们感兴趣的课题。由于计算机处理速度的迅速提高，计算机模拟已经和试验观测、理论分析并列成为本世纪科学研究的三种方法[19]。计算机模拟的数据(从模型中得来)可以用来比较、验证各种近似理论;同时，计算机模拟方法还可用来对试验和模型进行比较，从而提供了评估一个模型正确与否的手段。计算机模拟方法还有一个优点就是可以沟通理论和试验。某些量或行为可能是无法或难以在试验中测量的，而用计算机模拟方法，这些量可以被确切的计算出来[20]。分子动力学模拟方法更以其建模简单、模拟结果确凿的特征而倍受研究者的关注。

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

二．分子动力学基本原理

分子动力学方法是纳米计算力学的主要手段。它对多经典粒子系统的运动方程组进行数值积分，得到相空间轨道，进而研究系统的平衡热力学性质、结构动力学性质和非平衡传输性质等。在分子动力学模拟中，将原子视为质点或忽略内部结构的固体球。首先建立粒子系统的几何模型，通过描述原子间相互作用的势函数，根据给定的边界条件、初始条件求出系统中每一时刻单个粒子或原子的能量和所受到的力，代入牛顿动力学方程组求解原子的位置和速度，得到系统在相空间的运动轨迹。对足够大数量的粒子在足够长时间的结果进行统平均，则可以得到类似于宏观意义上的物理量和力学量。

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

a.粒子是相互作用的质点，运动由位置矢量和速度矢量来描述。粒子间的相互作用取决于彼此的空间位置。

b.系统无质量交换。即模拟过程中系统的原子数不变。

保守系综分子动力学常用于模拟能量守恒的孤立绝热系统。但是要确凿模拟实际纳米材料和器件的表面、界面等现象，真实反映粒子系统受外界约束状况下的物理行为，以及进行跨尺度的祸合模拟，保守系综分子动力学模拟明显是远远不够的。在此基础上，研究者提出了耗散系统与周边介质进行能量交换的不同理论和算法，对粒子系统进行温度、压力、粒子数、体积等不同控制，适合于不同系综的分子动力学模拟。

2

分子动力学模拟中假设系统所有粒子的运动遵循经典牛顿运动定律，且忽略原子背景电子云的量子效应，原子间的相互作用满足叠加原理。可见，分子动力学方法是一种对广义牛顿运动方程的近似的数值积分方法。

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

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

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

?为质点位置对时间的导数。在笛式中q为广义坐标，指定质点的空间位置。q卡儿坐标系下，由N个原子构成的模拟系统的拉格朗日方程可写为：

式中

子系统，可取拉格朗日量：

，为原子i的位置矢量，在三维空间n=N/3。对于互不影响的粒

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

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

式右端的两项分别表示系统的动能和势能。代入拉格朗日方程可得系统牛顿运动方程:

式中Fi即为原了i所受的内力，即由于系统中其他原子的作用而在原子i上表达出的合力。由牛顿运动方程建立线性微分方程组，给定初始条件(初始位置、初始速度)，求解该封闭方程，可得到确定解，求出原子运动的轨迹，即任意时刻原子的位置ri(t)和速度ri(t)。分子动力学模拟中的经典拉格朗日方程常用来计算原子的运动，得到单个原子的运动轨迹，描述原子系统的运动过程，以及反映在原子系统整体特征下的位移、变形、缺陷等。(2)哈密顿(Hamiltonian)运动方程：

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

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

3

.

式中Pi，为粒子i的动量，Pi=mii，哈密顿正则方程为:

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

在分子动力学模拟的基本方程中，拉格朗日运动方程更适合于求解原子系统运动的过程(求解原子速度和位置)，并能施加外部荷载如外力、约束、边界条件等，而保守统的哈密顿正则方程更适合于求解系统的动力学演化过程和热力学状态，如温度、热滚动。在实际的分子动力学模拟中，拉格朗日方程寻常与外部约束和边界条件一起，构成特定原子系统模型的初值问题;而Hamiltonian运动方程往往在保守系统的基础上，通过改变系统状态变量，与外部环境进行能量交换，构成不同系综的分子动力学模拟。

2.分子动力学原子间作用势函数分子动力学方法是通过原子间的相互作用势，依照经典牛顿运动定律求出原子运动轨迹及其演化过程。分子动力学计算的关键是原子间势函数的选取，它决定着计算的工作量以及计算模型与真实系统的近似程度，直接影响到模拟结果的成功与否。由于物质系统的繁杂性以及原子间相互作用类型的不同，很难得到满足各种不同体系和物质的一般性而又精度较高的势函数。所以针对不同的物质体系人们不断发展了大量的经验和半经验的势函数[22]。

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

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

r4

式中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。ε表征原子间的吸引/排斥强度，对原子中的一个从平衡位置r0移到无穷远处所作的功，对于不同物质的原子,ε各不一致，寻常为10-19J或ev量级。能量最低时刻的r值，即r0，表示原子间的平衡位置，此时原子间作用力为零。对上述方程求极值，可得r0=62??1.22σ

将上式代入式模型中原了间作用力为：

，可得Lconard一JoneS

2.Morse势函数

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