分子动力学模拟方法概述

SILICONLLEY■■【高新技术产业发展分子动力学模拟方法概述周晓平田壮壮忽晓伟(郑州大学西亚斯国际学院河南郑州451100)摘要:主要介绍分子动力学模拟的基本原理,阐述分子动力学方法的运动方程、数值解法、势函数、边界条件、适用系综以及体系相关性质的计算。最后指出分子动力学模拟方法的优势和发展方向。关键词:分子动力学;势函数;边界条件中图分类号:0414文献标识码:A文章编号:1671-7597(2012)1210040-021分子动力模拟计算的基本原理利用蛙跳法计算仅需储存(-)与r()两种资料,可节分子动力计算的基本原理,即为利用牛顿运动定律。在分子动力模拟中,体系原子的一系列位移是通过对牛顿运动方程省储存空间。其缺点是位置与速度不同步。这意味着在位置的积分得到的,结果是一条运动轨迹,它表明了系统内原子的定时,不可能同时计算动能对总动能的贡献位置与速度是如何随时间而发生变化2.2Gea算法[1先由系统中各分子位置计算系统的势能,按照经典力学Gear所提出的一种利用数值解的方法,称为校正预测法系统中任一原子i所受的力为势能的梯度( predictor- corrector method)。时间t+时的位置、速度等可由时间t的泰勒展开式预测得到(1)r"(+6)=()+8n()+12()+18b(+由牛顿第二定律可得i原子的加速度为F(+8)=()+n()+182b()(t+)=a()+b(1)将牛顿运动定律方程式对时间积分,可预测i原子经过时间t后的速度与位置:b(+)=6()+式中的(、a()、b()为()的1次、2次、3次微分。式(7)dt所产生的速度、加速度不是由牛顿运动方程解得的,所以并非(3)完全正确。可由所预测的位置产(+)计算所受的力及正确的加速度d(+)。设正确的加速度与预测的加速度之间的误差为r=r +p t+△a(+8)=a(+8)-a(+8)(8)式中,及"分别是粒子i的位置与速度,上标“0”为各可得各量的校正式为:物理量的初始值[1]r(+8)=r"(+8)+c△a(+2牛顿运动方程的数值解法为了得到原子的运动轨迹,必须解式(3)的牛顿运动方(+8)=(+6)+c△(+)(9)程,可采用有限差分法。有限差分法的基本思想就是将积分分a(+8)=a(+6)+c,△(+8)成很多小步,每一小步的时间固定为c。常用的有以下儿种算法:① Verlet算法;② Velocity- Verlet算法:③leapb(+)=b(+)+c3△a(+)frog算法(蛙跳算法);④ Beeman算法:⑤Gear算法式中,C0、C1c2、C3均为常数。以上为Gear的一次预测1 eap-frog算法和Gear算法由于使用简便,准确性及稳定性校正法,也可将此计算推展至更高次的校正高,节省储存空间等作者: photon优点,已被广泛采用3势函数2.1leap-frog算法势函数表明了原子间的相互作用。针对不同的计算物质Leap-frog算法速度与位置的数学式为:不同的模拟目的,势函数有不同的形式。计算结果的可靠性与势函数密切相关。在分子动力学发展初期,主要采用对势。随(+-=v(r--S)+sat着模拟体系的复杂性,逐渐出现了多体势,以弥补对势的不r(t+)=r()+Cp(t+-)3.1对势为了执行1eap-frog算法,必须首先由t-0.5时刻的速主要是 Lennard- Jones势(L-J势),又叫12-6势能,它度与t时刻的加速度计算出(t+8),然后由方程的数学表达式是(r+)=(+16(+)+5b()-16(t-06)(5)U(r)=4中国煤化工(10)计算出位置r(+6t)。时间为t时的速度可由式(6)算出,即式中,r为CNMHG能参数能中,r-2项是排斥项,r项是吸引项。当r很大时,L-J势能p()=(+)+p(-O)(6)趋近于零,表示当原子对相距很远时,彼此之间已经没有非键SILICON【高新技术产业发展】 VALLE结作用了。这个形式的势函数表达的作用力较弱,适合描述惰定义为性气体的固体和液体[2]多体势Daw和 Baskes首次提出了嵌入原子势(EAM),其体系的总式中,P为系统密度,dN为距离参考分子中心由r→r+ah势能可以表示为:球壳内的平均粒子数。径向分布函数可以理解为所模拟体系的区域密度与平均密U=∑F(p)+∑()(11)度的比[1]。当r值小时,距离参考分子很近,区域密度不同于平均密度;当r值大时,距离参考分子很远,区域密度与系统的式中,第一项F是嵌入能,表示原子核镶嵌在电子云背平均密度相同,径向分布函数接近于1景中的嵌入能;第二项是对势能,表示位于晶格点阵上的原子2配位数核之间的相互作用。多体势大都用来模拟金属的微观性质配位数( coordination number)是指某个原子的最近邻[2]。后来 Baskes等人修正了嵌入原子法,提出MEA势,可以原子个数,它显示了该原子周围的原子分布的密度大小,是微描述共价键材料观结构分析的辅助物理量[4]。其定义为:4周期性边界条件与最近镜像(13)执行分子动力计算时,通常选取一定数目N的分子或原n=4mJ”g()2bh子,将其放在一个立方体的盒子中,该盒子即为模拟系统。而式中,P是对应点的数密度,积分的上下限由径向分布函计算机最多只能模拟儿百到几千个粒子的系统,为了能够用较数g(r)得到,当万=0,n=rm时,得到第一配位圈的配位少的粒子数目来模拟真实的宏观体系,引入三维周期性边界条数件。在三维体系中,每一个单胞会被其它26个单胞所包围。当6.3扩散系数模拟的单胞中一个粒子由于力的作用离开这个单胞的时候,就分子动力学模拟中,采用两种方法来计算体系分子的扩散会有另一个和它对应的粒子运动到这个单胞中来,这样,模拟系数D,分别是 Einstein法和reen-Kubo法。 Einstein法是利的整个体系的粒子数就保持不变,密度也不变,符合实际的要对均方位移求斜率来求扩散系数,计算公式为求[3]。D=imaF()-7(0)(14)计算系统中分子间的作用力时,采取最近镜像的方法。计算分子A与B的作用力,是取与分子A和其最近的距离镜像分子Green-Kubo法是通过对速度自相关函数("ACF)的积分获B,而非计算系统中的分子A与分子B。因为在计算中利用了最近得扩散系数,计算公式为镜像的概念,因此必须采用截断半径的方法计算非键结的远程(15)作用力。定义为截断半径,则当>时,势能值趋近于零表示分子间的范德瓦耳斯作用力可忽略不计。截断半径最大不其中,D表示粒子i的扩散系数,r()、r(0)分别表示粒能超过盒子边长的一半,即r