分子动力学模拟基本原理和主要技术.pdf

1、第1 8 卷第1 期2 0 0 5 年3 月聊城大学学报(自然科学版)J o u r n a lo fL i a o c h e n gU n i v e r s i t y(N a t S c i)V o I 1 8 N o 1M a r 2 0 0 5分子动力学模拟基本原理和主要技术崔守鑫D 胡海泉D肖效光D黄海军2(”聊城大学物理科学与信息工程学院，山东聊城2 5 2 0 5 9；中国工程物理研究院流体物理研究所，四川绵阳6 2 1 9 0 0)摘要综述了分子动力学模拟技术的发展，从原子(或分子)运动的角度出发，简要地介绍了分子动力学方法的基本理论以及对一些感兴趣量的提取并且还简要地说明

2、了分子动力学的用途以及与其它方法的结合最后还指出了分子动力学模拟方法本身进一步的研究方向关键词分子动力学，有限差分法，相互作用势，结构分析技术中图分类号0 4 1 1 3文献标识码A文章编号1 6 7 2 6 6 3 4(2 0 0 5)0 1 0 0 3 0 0 5凝聚态物质的宏观现象是由极大量分子、原子的运动和相互作用的反映，在本质上都带有统计性质随着人们对凝聚态物质的微观结构的不断深入了解，人们愈加需要从统计角度探讨凝聚态物质的宏观现象的基本规律，寻找内在的联系现在凝聚态理论已经深入到知其所以然的阶段，在量子力学和统计力学基础上探讨宏观现象和微观结构的本质统计力学是研究凝聚态物质的微观结

3、构和宏观现象的桥梁，它的理论已取得了很大的成就，但凝聚态物理问题非常复杂，纯粹用解析理论来描述目前尚难以实现随着计算机的软硬件和计算技术的长足发展，利用计算机对物质的微观结构和运动进行数值模拟计算得到了迅速发展，并在此基础上发展了用数值运算统计求和的方法：分子动力学(M D)方法和蒙特卡罗(M c)方法这两种方法可以按照原子(或分子)的排列和运动的模拟以实现宏观物理量的计算，不仅可以直接模拟许多物质的宏观凝聚特性，得出与实验结果相符合或可比拟的计算结果，而且可以提供微观结构、粒子运动以及它们和物质宏观性质关系的明确图像，有利于从中提取新的概念和理论这两种方法又可视为“计算机实验”，利用粒子之间

4、的真实相互作用势，也可用半经验势计算凝聚态物质的结构和热力学性质，从而很方便的对相互作用与宏观性质之间的关系进行考察，据此提出有关物理现象的理论从而，这种“计算机实验”建立了理论与实验的桥梁，对新的理论和新的实验的建立将起愈来愈大的作用1 9 5 7 年，A l d e r 和W a i n w r i g h t 妇首先在硬球模型下，采用经典M D 方法研究气体和液体的状态方程，取得了尝试性的成功，由于计算机速度和内存的限制，早期模拟的空间尺度和时间尺度都受到很大的限制2 1 世纪8 0 年代后期，由于计算机技术的飞速发展和多体势函数的提出与发展，为分子动力学模拟技术注入了新的活力1 9 7

5、 2年，L e e s 和E d w a r d s 2 首先将分子动力学模拟应用到了非平衡态的研究，进一步扩展了M D 方法的应用范围目前，许多在实际实验中无法获得的微观细节，而在分子动力学模拟中可以很方便的得到，这种优点使得分子动力学模拟在新材料的研究和开发、预测材料的介观和宏观性质等领域显得非常重要，由此改变了过去理论计算与实验的关系1 分子动力学的运动方程分子动力学模拟的出发点是假定粒子的运动可以用经典动力学来处理，对一个由个粒子构成的孤立体系，粒子的运动由牛顿运动方程决定，也就是：m i d 2，r 出2=一V r V(r。，r。，r N)，式中，m，r，分别为第i国家自然科学基金重

6、点项目(0 2 9 9 0 4 0)、中国工程物理研究院科学技术基金资助课题(2 0 0 0 1 0 3)收稿日期：2 0 0 4 1 2 1 5 万方数据万方数据第1 期崔守鑫等：分子动力学模拟基本原理和主要技术3 1个原子的质量和位置V，一一a O r i，V(r，。，r N)为体系所处的势2 运动方程的数值积分计算机模拟方法的基点是利用现代计算机高速和精确的优点，对几百个以至上千个分子的运动方程进行数值积分有许多不同的积分方法，它们的效率和方便程度各异问题基本上就是用有限差分法来对二阶常微分方程进行积分常用的有以下几种方法2 1V e r l e t 算法 3 1V e r l e t

7、算法是在6 0 年代后期出现的，对扩散分子的质心运动的积分是最稳定的也是最常用的数值方法它运用t 时刻的位置和加速度以及t 时刻的位置来预测f+位置其积分方案，以三阶T a y l o r 展开为基础，由以下方程给出r(t+&)=2 r(t)一r(t 一&)+2 口(f)，(1)这里，为简单计，省略了i 速度可以按微分的基本法则得出：y()一，(t+艿t)一r(t 一8 t)2 8 t 这种算法的优点是占有计算机的内存小，并且很容易编程，但它的缺点是位置r(f+&)要通过小项(3 t 2)与非常大的两项2 r(t)和r(f 一乩)的差相加得到(见(1)式)，这容易造成精度损失并且从式(1)可以

8、看出，这种算法不是一个自启动算法，新位置必须由t 和f 一&时刻的位置得到2 2“蛙跳”(L e a p f r o g)算法 4 H o c k e y 提出的L e a p 一 r o g 算法是V e r l e t 算法的变化，这种方法涉及半时间间隔的速度，即：r(t+3 t)=r()+艿f 口(+3 t 2)，v(t+艿t 2)=v(t 一3 t 2)+3 t a()t 时刻的速度由下式给出：口(f)一 口(+艿t 2)+v(t 一3 t 2)1 2 这种算法与V e r l e t 算法相比有两个优点：(1)包括显速度项；(2)收敛速度快，计算量小这种算法明显的缺陷是位置和速度不同

9、步2 3V e l o c i t y V e r l e t 算法 5 1这种算法可以同时给出位置、速度和加速度，并且不牺牲精度，给出了显速度项，计算量适中，目前应用比较广泛2 4G e a r 的预测一校正算法这种算法分为三步来完成：首先，根据T a y l o r 展开，预测新的位置、速度和加速度然后，根据新的计算的力计算加速度这个加速度再由与T a y l o r 级数展开式中的加速度进行比较两者之差在校正步里用来校正位置和速度项这种方法的缺点就是占有计算机的内存大除了上述提及的几种方法外，还有B e e m a n 算法 6。、R a h m a n 刀等3 周期性边界条件和长程力即

10、使是使用现代的巨型计算机，M D 方法还是只能用于粒子数大约是几百到几千的系统这就引起一个问题：用这样少量的粒子，如何来模拟宏观体系?为了解决这个问题，引入了周期性边界条件 8 采用这种方法，模拟体系实际上是由基本单元(也称为模拟计算元胞)在各个方向上重复叠合而成但在模拟中只需保留基本单元，所有其它单元与基本单元由平移对称性关联在处理粒子之间的相互作用时，通常采用“最小影像”约定这个约定是在由无穷重复的M D 基本模拟计算元胞中，一个粒子只与它所在的基本元胞内的另外N 一1 个(设在此元胞内有N 个粒子)中的每个粒子或其最邻近影像粒子发生相互作用实际上，这个约定就是通过满足不等式条件r,L 2

11、 来截断位势(n 为截断半径，是元胞的边长)通常的数值应当选得很大，以避免有限尺寸效应，但这样会增大计算量，通常采用对相互作用势的修正来近似处理4 势函数M D 模拟结果准确与否的关键在于对系统内的原子之间相互作用势函数的选取，总的来说，原子(或分子)之间的相互作用势的研究进展一直很缓慢，在一定程度上制约了M D 方法在实际研究中的应用原子间的势函数的发展经历了从对势到多体势的过程，对势认为原子之间的相互作用是两两之间的作用，与 万方数据万方数据3 2聊城大学学报(自然科学版)第1 8 卷其它原子的位置无关，而实际上，在多原子体系中，一个原子的位置不同，将影响其它原子间的有效相互作用，所以，多

12、体势能更准确地表示多原子体系势函数4 1 对势在分子动力学模拟的初期，经常采用的就是对势。对势可以分为间断对势 1 和连续对势，而连续对势主要有以下几种：L e n n a r d J o n e s(L J)势、B o r n l a n d e(B L)势、M o r s e 势和J o h n s o n 势等，其中，L J 势是为描述惰性气体分子之间的相互作用而建立的，因此它表达的力比较弱，描述材料的行为也比较柔韧，也可以用来描述过渡金属原子之间的相互作用B L 势是用来描述离子晶体离子之间的相互作用的M o r s e 势和J o h n s o n 势多用于描述金属原子之间的相互作

13、用对势虽然简单，得到的结果往往也符合某些宏观物理规律，但其缺点是必然导致C a u c h y 关系，所以，对势实际上不能准确地描述晶体的弹性性质4 2 多体势多体势是在2 0 世纪8 0 年代初期开始出现的，1 9 8 4 年D a w 和B a s k e s 首次提出了原子嵌入(E A M)势9 此势的基本思想是：把晶体的总势能分成位于晶格点阵上的原子核之间的相互作用对势和原子核镶嵌在电子云背景中的嵌入能(多体相互作用)两部分，其中，对势和多体作用势的函数形式往往根据经验选取基于E A M 势的势函数还有很多种 10 n ，这些多体势大多用于金属的微观模拟，此外，还有许多形式的多体势函数

14、形式，1 9 8 7 年J a c o b s e n 等人 1 2 在等效介质原理的基础上提出的另一种多体势函数形式，由于其简单、有效，也得到了广泛地应用5 分子动力学模拟的系综平衡态分子动力学模拟是在一定的系综下进行的，经常用的平衡系综是N V T 或N P T 系综在这两种系综中，牵涉到控制温度和压力的几种技术，分别介绍如下5 1 控温方法N(1)速度标定 1 引系统的温度与动能存在如下的关系：E t 一芝：m iI 口中z 一(3 N 一G)K B T 2，式_中，是原子数，c 是约束数，K B 是B o l t z m a n n 常数，口i 是原子i 的速度由于系统的温度和动能存在

15、这样的关系，所以一种最简单和最直观的方法是直接对速度进行标定，这种方法的基本思想是：如果t 时刻的温度是丁()，速度乘以因子A 后，温度的变化为：丁一(A 2 1)丁(t)，其中，A 一丁丁(f)，丁为所控制体系的温度(2)B e r e n d s e n 恒温槽方法 1 4 3 这种控温方法假设所模拟的体系与一个恒温槽连在一起，则两者之间J _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _-。-_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 一就可以通过热交换而使模拟体系达到恒温的目的，方法如下，定义

16、一个参数A：A 一1+f(T T。)r，式中，参数r 表征系统与恒温槽之间的热交换速率，为M D 的时间步长，那么通过7-。，=A 口o-a 校正即可保持体系的温度在丁。附近振动，而参数r(通常取为0 1 o 4 加)则可用于控制这个振动幅度(3)N o s 6 一H o o v e r 方法 i s,1 6 这种方法是通过改变模拟体系的H a m i l t o n i a n 来实现控温的，因而有更强的物理意义，其基本思想就是在H a m i l t o n i a n 加入一个假想的项来代表一个恒温源，具体做法如下：日=芝：m，口；2+y(r)+Q f 2 2+g K T l n S，式

17、中，s 和f 分别是假想项的坐标和动量这样体系的微分方程就变为：口i=d r,d r，口r 一一(d V d r+m m f)m i，d f d t 一(：m r 口i g K 丑丁)Q 式中，g 为体系的自由度，Q 为一个可调参量，表征着假想项的质量，丁为温度在这三种控温方法中，速度标定是一种非真实的物理效应，但这种方法可以使系统很快达到平衡，在经典分子动力学方法中，这是一种比较常用的控温方法B e r e n d s e n 控温方法是通过系统和恒温槽进行热交换来控制温度，此方法的优点是非常简单，应用起来非常方便N o s 6 一H o o v e r 控温方法是基于统计力学而提出来的一种

18、控温方法，当系统和恒温槽进行热交换时，在系统中粒子出现的几率遵从统计力学规律，所以这是一种真实的物理效应5 2 控压方法用M D 方法研究压力的诱导相变和结构重构，在等压下模拟比在等体积下更容易实现在等压模拟 万方数据万方数据第1 期崔守鑫等：分子动力学模拟基本原理和主要技术3 3下，可以通过改变模拟元胞的三个方向或一个方向的尺寸来实现体积的变化，类似于温度控制的方法，也有许多方法用于压力控制，在实际应用中，有以下两种较为常用的压力控制方法(1)B e r e n d s e n 方法m 这种控压方法的基本思路与基本控制方程均与他的控温方法类似，如下：卢一 14-出(P P。)r z-3“3，

19、式中，y 是一个可调参数，尸。为所控制的压力然后，在每一次的M D 迭代中，粒子坐标X，Y，z 均用卢相乘，得到新的坐标，即可实现对压力P 的控制该方法的缺点是可能导致在很长的模拟时间步长内，保持着起伏(2)P a r r i n e l l o R a h m a n 方法 15。这种方法是通过改变体系的L a g r a n g i a n 来实现的，新的L a g r a n g i a nNN定义如下：L 一m，羹G 矗z 一V(I i，)-4-W T r h 左z 一户力，其中，一 口，6，c)是模拟单胞的基矢，一一i=1 两G=h h，原子坐标n 为r，=h s，V 为原子之间的相

20、互作用势，0 为单胞体积，n 表示矩阵的迹，是一个可调参量这种方法的优点是模拟单胞不仅大小可变，形状也可变6M D 模拟的初始条件M D 模拟进入对系统微分方程组做数值求解的过程时，需要知道粒子的初始位置和速度的数值，不同的算法要求不同的初始条件常用的初始条件可以选择为：(1)令初始位置在差分划分网格的格子上，初始速度则从玻尔兹曼分布随机抽样得到(2)令初始位置随机地偏离差分划分网格的格子，初始速度为零(3)令初始位置随机地偏离差分划分网格的格子，初始速度从玻尔兹曼分布随机抽样得到7 热力学性质的计算经典力学的一个基本前提是，只要知道了物理体系在相空间的运动，就可以导出其所有的宏观性质在分子动

21、力学中，对模拟得到的位形进行结构分析是很重要的，比如在研究晶体的熔化时，需要及时分析它的结构变化就结构分析而言，目前常用的方法有：径向分布函数、静态结构因子和配位数等方法7 1径向分布函数(r a d i a ld i s t r i b u t i o nf u n c t i o n(R D F)【1 8 径向分布函数可以由下式来定义：g(，)=y 1 l(r)4 1 t r 2 A r，式中，扎(r)表示距离原点r 到r 4-A r 之间的平均粒子数，y 是模拟的体积R D F 表示的物理含义是：在空间位置r 点周围的体积元中单位体积内发现另一个粒子的几率由此看出，R D F 表征着结构

22、的无序化程度7 2 静态结构因子(s t a t i cs t r u c t u r ef a c t o r(S S F)1 妇N静态结构因子的表示式为：5(忌)=I 乏：e x p(i K，)l，其中，N 是总原子数，K 是倒格矢，r 为原子。j。=l的位置矢量对理想晶体而言，S S F 为1，而对理想流体，则为0 7 3 配位数(c o o r d i n a t i o nn u m b e r(C N)所谓配位数，就是原子的第一最近邻原子的个数配位数的多少，表示了该原子周围的原子分布的密度大小，它经常在一些分析中作为结构分析的一种辅助手段8M D 方法的进一步研究方向经典分子动力学

23、方法存在两个缺陷：(1)元胞的体积和形状保持不变，限制了M D 方法的应用；(2)不适合含有自由电子的系统，因而对金属等系统的计算的结果不理想为了克服经典分子动力学方法的局限性，2 0 世纪8 0 年代初，A n d e r s o n z o ,P a r r i n e l l o 等人E z l 发展了可变元胞分子动力学方法由于能带论中的密度泛函(D F)方法是研究凝聚态物质电子结构的有效方法，1 9 8 5 年C a r 等人 2 2 提出了第一原理的M D 方法但是，进一步拓宽分子动力学方法的应用面显然是一个热门的发展方向对M D 方法本身而言，最重要的两个要素是初始位形的给定和相互

24、作用势函数的选取M D 模拟时间足够长，初始条件的选择不会影响计算的结果，但是选择合理的初始条件可以加快系统趋于平衡，可以节省机时影响M D 计算结果精确程度的最主要的因素是相互作用势的精确性，人们对提高势函数的精确性进行了大量的研究，提出了许多势 万方数据万方数据3 4聊城大学学报(自然科学版)第1 8 卷函数的形式，但是，对大多数原子而言，努力找到即准确而形式又不太复杂的势函数，仍然是一个很有挑战性的课题我们认为，把第一原理、量子化学分析和参数拟合相结合是势函数研究的最好的方法，N a n x i a nC h e n 等人3 就用此方法研究了氯化钠晶体离子之间的相互作用势，并研究了微结构

25、的变化，这是一个值得进一步研究的方向参考文献 1 B J A l d e r，T E W a i n w r i g h t P h a s et r a n s i t i o n f o rah a r ds p h e r es y s t e m J J C h e m P h y s，1 9 5 7，2 7 1 12 0 8 I2 0 9 2 A W L e e s，S F E d w a r d s T h ec o m p u t e rs t u d yo ft r a n s p o r tp r o c e s su n d e re x t r e m ec o n d

26、i t i o n s J 1 J P h y s C 1 9 7 5，C 5：19 2 1 19 2 9 3 L V e r l e t C o m p u t e r e x p e r i m e n t s c l a s s i c a lf l u i d s I T h e r m o d y n a m i c a lp r o p e r t i e so fL e n n a r d J o n e a m o l e c u l e s J P h y s R e v，1 9 6 7，1 5 9：9 8 1 0 3 4 R W H o n e y c u t t T h

27、ep o t e n t i a lc a l c u l a t i o na n ds o m ea p p l i c a t i o n s J M e t h o d si nC o m p u t a t i o n a lP h y s i c s，1 9 7 0 9：1 3 6 2 1 1 5 w c S w o p e，H C A n d e r s o n，P H B e r e n s，K R W i l s o n Ac o m p u t e rs i m u l a t i o nm e t h o df o rt h ec a l c u l a t i o no

28、 fe q u i l i b r i u mc o n s t a n t sf o rt h ef o r m a t i o no fp h y s i c a lc l u s t e r so fm o l e c u l e s：a p p l i c a t i o n st os m a l lw a t e rc l u s t e r s J J C h e m P h y s，1 9 8 2，7 6：6 3 7 6 4 9 6 D B e e m a n S o m em u l t i s t e pm e t h o d sf o ru s ei nm o l e c

29、 u l a rd y n a m i c sc a l c u l a t i o n s J J C o m p u t P h y s，1 9 7 6，2 0 1 1 3 0 1 3 9 7 A R a h m a n C o r r e l a t i o n si nt h em o t i o no fa t o m si nl i q u i d sa r g o n J P h y s R e v A。1 9 6 4。1 3 6：4 0 5 4 1 1 8 N M e t r o p o l i s，A W R o s e n b l u t h，M N R o s e n b

30、 l u t h，A H T e l l e r，E T e l l e r E q u a t i o no fs t a t ec a l c u l a t i o n sb yf a s tc o m p u t i n gm a c h i n e s口 J C h e m P h y s，1 9 5 3，2 11 10 8 7 10 9 2 9 M S D a w，M I B e s k e s E m b e d d e da t o mm e t h o dd e r i v a t i o na n da p p l i c a t i o nt oi m p u r i t

31、 i e s，s u r f a c e s，a n do t h e rd e f e c t si nm e t a l s J P h y s R e v B，1 9 8 4 2 9：84 8 6 84 9 5 1 0 s M F o i l e s，M I B a s k e s，M S D a w E m b e d d e d a t o m m e t h o df u n c t i o nf o rt h ef c cm e t a l sC u，A g，A u，N i P d，P ta n dt h e i ra l l o y s J P h y s R e v B，1

32、9 8 6，3 3 t 79 8 3 79 9 1 1 1 3M M a n n i n e n，R A J o h n s o n I n t e r a t o m i ci n t e r a c t i o n si ns o l i d s J P h y s R e v B，1 9 8 6，3 4：84 8 6 84 9 5 1 2 3K W J a c o b s e n，J K N o r s k o v，M J P u s k a I n t e r a t o m i c i n t e r a c t i o n s i n e f f e c t i v e-m e d

33、 i u m t h e o r y J P h y s R e v B，19 8 7，3 5：7 4 2 3 74 4 2 1 3 K H H o f f m a n n，M S c h r e i b e r C o m p u t a t i o n a lP h y s i c s M B e r l i nH e i d e l b e r g：S p r i n g e r V e r l a g，1 9 8 6，2 6 8 3 2 6 1 4 H J C B e r e n d s e n，J P M P o s t m a，W F V G u n s t e r e n M o l e c u l a rd y n a m i c sw i t hc o u p l i n gt oa ne x t e r n a lb a t h J J C h e m P h y s，