第40卷第5期 当 代 化 工 Vol.40，No.5 2011年5月 Contemporary Chemical Industry May，2011
分子模拟的方法及其应用
李春艳，刘 华，刘波涛
（中北大学物理系，山西 太原030001）
摘 要： 综述了分子模拟的方法及其在材料科学中的应用，介绍了分子模拟的基本原理、牛顿运动方程及其有限差分算法、势函数的发展、平衡态系综及其温度和压力的控制，最后还指出了分子模拟的应用及其进一步的研究方向。
关 键 词： 分子模拟; 有限差分法; 原子间相互作用势; 材料科学
中图分类号： O 642 文献标识码： A 文章编号： 1671-0460（2011）05-0494-03
Method of Molecular Simulation and Its Application
LI Chun-yan, LIU Hua, LIU Bo-tao
（North University of China, Shanxi Taiyuan 030001, China）
Abstract: Method of molecular simulation and its application were reviewed. The principle of molecular simulation, Newtonian equation of motion were introduced as well as its finite difference technique, progress of potential function, realization of equilibrium ensembles and the temperature and pressure control. Directions of researches and applications of molecular simulation were put forward.
Key words: Molecular simulation; Finite difference technique; Interatomic potential; Material science
在高分子领域中计算机模拟的应用从20世纪60年代至今已发展到了一个崭新的阶段。这个新阶段的特点是，计算机模拟不仅能提供定性的描述，而且能模拟出高分子材料结构与性能的定量的结果；模拟分子体系算法的发展，主要是从描述简单的非真实分子体系的算法到复杂的真实分子体系的算法。该发展克服了用蒙特卡洛法仅能够描述不同温度下分子结构的特征，但不能描述随温度变化而引起体系宏观物理性质的变化中分子结构的变化过程。
[1-2]
1.2 分子动力学方法(简称MD)的基本原理
分子动力学方法(MD)是另一种主要的计算机模拟方法，在材料科学、物理、化学等学科的各个领域得到广泛应用。模拟过程是在一定体系及分子势能函数条件下，从计算分子间作用力着手，求解牛顿运动方程，得到体系中各分子微观状态随时间的变化，再将分子的位置和动量组成的微观状态对时间平均，即可求出体系的压力、能量、粘度等宏观性质以及组成粒子的空间分布等微观结构。
[8]
1 分子模拟方法
1.1 蒙特卡洛方法(简称MC)的基本原理
MC法早在20世纪50年代就被应用到高分子科学的研究之中。先驱性的工作是Wall在为研究高分子链的排除体积问题时所做的蒙特卡洛模拟
[3-5]
2 分子模拟的计算
2.1 分子动力学的运动方程
分子动力学模拟的出发点是假定粒子的运动可以用经典动力学来处理, 对一个由N 个粒子构成的孤立体系, 粒子的运动由牛顿运动方程决定, 也就是
d 2r
m =-∇V (r , r , r ) （1）
i 12N
d t 式中：m i ，r i —第i 个原子的质量和位置；
▽i = -d/dri ,V (r 1, r 2, r N ) —体系所处的势。
。基
本模拟过程是在一定系统条件下，将系统内粒子进行随机的位移、转动，或粒子在两相间转移位置。根据给定的分子势能函数，进行粒子间内能的加和。采用Metropolis取样方法
[6-7]
，生成一系列体系的微观
粒子随机构型，从而逐渐趋近于平衡时的Boltzmann分布。目前蒙特卡洛的应用几乎涉及到高分子科学中的各个领域，为现代高分子科学理论基础的建立和发展起到了十分重要的推动作用。
2.2 分子动力学运动方程的求解
在分子动力学中，为了得到原子的运动，一般采用以下几种有限差分法来求解运动方程。常用的有：
基金项目：国家自然科学基金（批准号：50730009）。
作者简介：李春艳，女，硕士研究生，主要从事分子力学模拟研究。E-mail:[email protected]。
第40卷第5期 李春燕，等：分子模拟的方法及其应用 495
Verlet算法是在60年代后期出现的, 对扩散分子的质心运动的积分是最稳定的也是最常用的数值方法。Swope提出的Velocity–Verlet算法可以同时给出位置、速度与加速度，并且不牺牲精度。这种算法的优点是给出了显速度项，并且计算量适中，应用比较广泛。Hockney提出的“蛙跳”Leap–forg 算法是Verlet算法的变化，这种算法与Verlet算法相比有两个优点:（1）包括显速度项;（2）收敛速度快, 计算量小。这种算法明显的缺陷是位置和速度不同步。Beeman提出的Beeman算法运用了更精确的速度表达式。它的表示式比Verlet算法复杂得多，因此计算量较大。Gear的预测——校正算法分为三步来完成:首先, 根据Taylor展开, 预测新的位置、速度和加速度。然后, 根据新的计算的力计算加速度。这个加
[13]
[12]
[11]
[10]
[9]
模型主要有EAM势、FS势、TB势等。 3.2.1 EAM势
Daw和Baskes在1984年首次提出了嵌入原子法（EAM）。EAM势的基本思想是把晶体的总势能分成两部分：一部分是晶格点阵上的原子核之间的相互作用对势，另一部分是原子核镶嵌在电子云中的嵌入能，它代表多体相互作用。在嵌入原子法中，系统的总势能表示为：
E =F ρ+1φr
（6）∑i (i )∑ij (i , j )tot
i
[20]
2
ij
式中：F 是嵌入能：第2项是对势项，根据需
要可以取不同的形式。ρi 是除第i 个原子以外的所有其它原子的核外电子在第i 个原子处产生的电子云
ρi (r ij )，密度之和，可以表示为：ρi =∑其中，ρi （r ij ）j ≠i
是第j 个原子的核外电子在第i 个原子处贡献的电荷
速度再由与Taylor级数展开式中的加速度进行比较。
密度，r ij 是第i 个原子与第j 个原子之间的距离。对于
这种方法的缺点就是占有计算机的内存大。
不同的金属，嵌入能函数和对势函数需要通过拟合金属的宏观参数来确定。 3 分子间的作用势 3.1 对势
在分子动力学模拟的初期，人们经常采用的是对势模型有以下几种：
（1）在1957年Alder和Wainwright使用间断对势: 当r ≤r 1时,E =∞;当r ＞r 1时, E =0。
（2）连续对势
Lennard-Jones(L-J )势是最老的原子间作用势之一，是在1924年由J.E.Lennard-Jone提出来的，
[15]
[14]
3.2.2 FS势
Finnis等根据金属能带的紧束缚理论，形成了
[21]
一种在数学上等同于EAM的势函数，并给出了多体势的函数形式，把嵌入能函数设为平方根形式，即
ρi (r ij )项分别为：FS势。其中的多体项及对势ρi =∑j ≠i
ρ(r ij )=∑A k (R k −r ij )H (R K −r ij )
3
k =1
6
2
3
σij (r ij )=∑a k (r k −r ij )H (r k −r ij ) （7）
k =1
L-J 势的基本形式为：
φ(R )=4ε⎡⎛σ⎞−⎛σ⎞⎤ （2）
⎢⎜⎟⎜⎟⎥LJ
⎢⎝R ⎠⎝R ⎠⎥⎣⎦
式中：σ、ε—经验常数。
12
6
式中：当x >0时，H (x )=0；当x
基于EAM势的势函数还有很多种的非球形对称。于是Baskes等
[24-26]
[22-23]
（3）Born-Lande(B-J )势的形式如下： 2
Z i Z j b
φij (r ij )=e +m （3）
r ij 4πεr ij （4）Morse势：
根据双原子分子的振动谱，提出了指数形式的相互作用势：
φ(r )=Ae (−αr )−Be (−βr ) （4）
它有4个参数A , B , α和β，与L-J 势的普适形式相类似。
（5）Johnson势的形式如下：
φij (r ij )=−A n (r ij −B n )+C n r ij −D n （5）
3[19][18]
[17]
[16]
。为了把
将EAM势推广到共价键材料中去，需考虑到电子云
提出了修正型嵌入
原子核法(MEAM)。总的电子密度为下式：
3
(ρ)2=t (h )(ρ(h ))2 （8）
i
∑
h =0
i i
给出。其中t (h ) i 是常系数。在MEAM里，镶嵌能项采取了如下形式：F (ρ)=AEρln ρ,其中A ，E 是系数，因物质不同而不同。
3.2 多体势
多体势于20世纪80年代初期开始出现。多体势于20世纪80年代初期开始出现。常用的多体势函数
4 分子动力学模拟中平衡系综的控制
方法
平衡态分子动力学模拟是在一定的系综下进行的, 经常用的平衡系综是NVT或NPT系综。在这两种系综中, 牵涉到控制温度和压力的几种技术,分别介
496 当 代 化 工 2011年5月
绍如下： 4.1 控温方法
在平衡系综中，经常要把温度调整到期望值。而体系的有效温度是由动能决定的：
k =∑m i v i / 2=fk B T /2
式中：f 为体系的自由度，K B 为波尔兹曼常数。 控温常用的方法有: （1）速度标定法
系统的温度与动能存在如下的关系：
E k =∑m i v i 2/2=(3N −N C )K B T /2 （9）
i N
[27]
2
与他的控温方法类似, 如下:
μ=[1+△t (P -P 0) γ/Г]1/3
式中γ是一个可调参数, P 0为所控制的压力。然
后, 在每一次的分子动力学模拟迭代中, 粒子坐标
x , y , z 均用μ相乘, 得到新的坐标, 即可实现对压力P 的控制。
4.2.2 Andersen方法
这种方法是通过改变体系的Lagrangian来实现的, 新的Lagrangian定义如下：
i i N N ii i
L =m s G h /2−V (r )+WT h ' h /2−p Ω（12）
[32]
∑
i i
∑∑
i =1j >1
ij r
式中：N —原子数；
N C —约束数；
K B —波尔兹曼常数； v i —原子i 的速度。
其中：h =﹛a , b , c ﹜—模拟单胞的基矢；
G =h ’h ,原子坐标r i 为r i =hs ； V —原子之间的相互作用势； Ω—单胞体积； T r — 矩阵的迹；
这种方法的基本思想是：如果t 时刻的温度是T （t ）, 速度乘以因子λ后, 温度的变化为：
2
ΔT =（λ-1）T （t ）
T 为所控制体系的温度。 （2） Berendsen热浴法法
这种也称恒温槽方法，是假设所模拟的体系与一个恒温槽连在一起, 那么两者之间就可以通过热交换而使所模拟的体系达到恒温目的, 方法如下:
λ=,式中,参数Г表征系统与恒温槽
[28]
W — 一个可调参量这种方法的优点是模拟单胞不仅大小可变,
形状也可变。
5 分子模拟在材料科学中的应用
分子动力学模拟计算液态金属的结构及热力学性质在国外得到了广泛地应用。其中波兰的Janusz 等利用多体相互作用，借助分子动力学计算了Ag、Au、Cu、Ni等面心立方金属的热力学性质。Holzman 等采用EAM 势计算了面心立方金属液-气界面的特殊性，得到的密度、内能、结构因子等结果与实验值相一致。N.Tajima 等采用EAM 势模拟了面心立方金属的点缺陷。K.Kadau等运用分子动力学模拟了烧结Fe-Ni 纳米颗粒的马氏体转变。Tsai 采用第一性原理的分子动力学方法，研究了GaN 表面气态Ga 和N 原子的吸附特征，模拟结果显示：引入的Ga 原子相对不活泼，而引入的N 原子却很活泼；Ga在GaN 表面的吸附系数要比N 原子高很多。
[37]
[36]
[35]
[34]
[33]
之间的热交换速率,Δt 为分子动力学模拟的时间步
长,通过v new =λv old 校正即可保持体系的温度在T 0附近振动, 而参数Г可用于控制这个振动幅度。
（3）Nose-Hoover方法法
[29-30]
Nose-Hoover热浴法是通过改变模拟体系的哈密顿量来实现控温的,其基本思想就是在哈密顿量里加入一个假想的项来代表一个恒温源, 具体做法如下：
22
H =∑m i v i /2+V (r )+Q ζ/2+gKT ln S （10） 式中：S 和ζ分别是假想项的坐标和动量。
近几年，在研究液态金属方面，中国科学院冶这样体系的微分方程就变为：
[38]
/r +m i v i z )/m i ， 金研究所陈柳等对晶体的生长过程进行了分子动v i =d r i /d r ，a i =-(d V /d
[39]
力学模拟研究。王金照等利用分子动力学方法模2
（11）
/t =(∑m i v i -gK B T )/Q d z /d
式中：g — 体系的自由度；
Q — 一个可调参量, 表征着假想项的质量, T 为温度。
拟研究了铜铝合金的比热。戎咏华等对铁锰合金层错位能进行了计算。张建民等采用改进嵌入原
[41]
[42]
[40]
4.2 压力控制
分子动力学模拟通常要求模拟体系的压力保持不变，而体积可以随温度的变化而改变，即进行等压模拟。这可以通过改变模拟原胞的3个方向或一个方向的尺寸来实现体积的变化。通常有： 4.2.1 Berendsen方法
这种控压方法的基本思路与基本控制方程均
[31]
子法计算了铜晶体的表面能。张朝阳等采用
Discover／Materials Studio 中的MD（NVT）及MM （COMPASS力场）方法，通过对两种氟聚合物在TATB 不同晶面上吸附行为的动力学模拟和结构优化发现：二者相互作用时,TATB原子位置变化的大小次序为：（010）>（100）>（001）；而氟聚合物的展开程度按以上次序增加。宋小艳等利用三
[43]
第40卷第5期 李春燕，等：分子模拟的方法及其应用 497
维MC 技术模拟了较完整的单相材料正常晶粒长大过程，得到了晶粒长大动力学和拓扑学的全面信息。莫春立等以单相铁素体不锈钢为研究对象，利用基于实验数据的模型，对焊接影响区晶粒长大过程进行了MC 动态模拟。清华大学的郭靖原等运用MC 方法模拟了金属凝固与晶体生长过程。孙雷等
[47]
[46]
[45]
[44]
[20] Daw M S,Baskes M I. Embedded atom method derivation and application
to impurities, surfaces, and other defects in metals[J]. Physical Review B, 1984,29（12）: 8486-8495.
[21] Finnis M W，Sinclair J E．A simple empirical n —body potential for
transition—metal[J]．Philosophic Magazine A，1984，50（1）：45. [22] Foiles S M, Baskes M I,Daw M S. Embedded-atmo-method function for
the fcc metals Cu, Ag, Au,Ni, Pd, Pt and their slloys[J]. Physical Review B, 1986, 33（12）: 7983-7991.
[23] Manninen M,Johnson R A. Interatomic interactions in solids[J]. Physical
Review B,1986, 34（12）:8486-8495.
[24] Baskes M I. Modified embedded-atom potentials for cubic materials and
impurities[J]. Physical Review B, 1992, 46（5）:2727-2742. [25] Baskes M I, Johnson R A.Modified embedded-atom potentials for hcp
metals. Modelling Simul Mater Sci Eng[J]. Phys.Rev.1994,2: 147-163. [26] Baskes M I.Calculation of the behaviour of Si ad-dimers on
Si(001).Modelling Simul Mater Sci Eng[J]. Phys.Rev., 1997,5: 149-158. [27] Hofh曲nn K H，Schreiber~L C0mputat．oml Physics．Berlm Heidelberg：
Spring—verlag[J].Phys.Rev.,1996:268.
[28] Berendsen H J C，Postma J P M，Gunsteren W F V，et a1.[J]．chdTl Phys，
1984，81：3684.
[29] N05e S J.[J] Ch1 Phys，1984，8l：5. [30] H00ver W G．[J].Ph”RA，1985，3l：1695.
[31] Berendsen H J C，Postma J P M，Gunsteren W F V，et a1．[J].J chdTl
Phys，1984，81：3684.
[32] AnderseIl H C.Chem Phys[J]．Journal of Computational
Physics,1980:721-2384.
[33] Anuse J,Holender L M.[J]. Phys Rev B, 1990,141（12）: 8054. [34] Holzman L M,Adams J B,et al.[J]. Mater Res,1991,6（2）:343.
[35] Tajima N,Takai O,Kogure Y, et al. Computer simulation of point defects in
fcc metals using EAM potentials[J]. Computational Materials Science, 1999, 14:152
[36] Kadau K, Entel P, Lomdahl P S. Molecular-dynamics study of martensitic
transformations in sintered Fe-Ni nanoparticles[J]. Computer Physics Communications, 2002, 147: 126.
[37] Tsai M H．[J].Computer Physics,2002，147：130.
[38] 陈柳，程兆年，唐鼎元.晶体生长过程的分子动力学模拟研究[J].人
体晶体学报，1999，28（2）：109.
[39] 王金照，杨春，陈民,等.铜银合金比热的分子动力学模拟[J].工程热
物理报，2002,23（3）：274.
[40] 戎咏华，孟庆平，何刚,等.Fe-Mn 合金层错能的嵌入原子法计算[J].
上海交通大学学报，2003,37（2）：171.
[41] 张建民，徐可为，马飞.用改进嵌入原子法计算Cu 晶体的表面能[J].
半导体学报,2003,52（8）：1993.
[42] 张朝阳，舒远杰，黄辉，等．爆炸与冲击[J].金属学报，2004，24
（6）：563.
[43] 宋小艳．基于计算机模拟的正常晶粒长大的三维特征二维截面表征
[J]．材料研究学报，1998，12（3）：245
[44] 莫春立．铁素体不锈钢焊接热影响区晶粒长大过程的模拟[J]．金属
学报，2001，37（3）：307.
[45] 郭靖原，洪泽恺，陈民，等．金属凝固与晶体生长过程的蒙特卡罗
模拟[J]．工程热物理学报，2001，22（6）：725.
[46] 孙雷，杜刚，刘小彦，等．蒙特卡罗方法模拟金属一半导体基础的
直接隧穿效应（英文）[J]．半导体学报，2001，22（11）：1364. [47] 姜寿文，赵国群，关小军，等．冷轧薄钢板退火过程组织演变的
Monte Carlo模拟[J]．钢铁研究学报，2002，14（5）：53.
模拟了金属-半导体接触的直接隧穿效应。姜寿文等对冷轧薄钢板退火过程的组织演变进行了MC 模拟。
参考文献：
[1] E.B.Wilson, J.C.Decius, P.C.Cross. Molecular Vibrations, McGraw-Hill[M].
New York:1995:1-51.
[2] T.Shimanochi.Vibration Spectroscopy and Its Chemical application[M].
Tokyo:Univ.of Tokyo, 1997:105-406.
[3] F.T.Wall.Strain fluctuations and elastic constants [J].Chem.Phys, 1953,21
（19142）:98-235.
[4] F.T.Wall,L.A.Hiller, Jr., D.J.Wheeler. A new molecular dynamics
method[J].Chem. Phys., 1954,22（1036）:76-343.
[5] F.T.Wall,et al.. Interatomic interactions in solids[J].Chem.Phys.,1955:
23-231.
[6] Frenkel D，Smit B．Understanding molecular simulation：From algorithms
to applications[M]．SanDiego：Academic Press，1996:21-196. [7] 分子模拟：从算法到应用[M]．汪文川，等译．北京：化学工业出版
社，2002.
[8] Heermann D W.ComDuter simulation methods in theoretical physics.Berlin
Heidelberg：Springer—Verlag[J].Phys.Rev.,1990:81-345.
[9] L.Verlet.Computer ‘experiments’ classical fluids.I.Thermodynamical
properties of Lennard-Jonesmolecules[J].Phys.Rev.,1967（159）:98-103. [10] 孙伟，常明，杨保和.分子动力学模拟纳米晶体铜的结构与性质[J].
物理学报，1998,47（4）：591-597. [11]
R.W.Honeycutt.Thepotential
calculation
and
some
applications[J].Mcthods in Computational Physics,1970,9:136-211. [12] Beeman D.Some multistep methods for use in molecular dynamics
calculations.Journal of Computational Physics[J]. Phys.Rev.,1976,20: 130-139.
[13] Gear CW.Numerical imitial value problems in ordinary differential
equation[M]. Englewood Cliffs:Prentice--Hall,1971.
[14] Alder B J, Wainwright T E. Phase transition for a hardsphere
system[J] .The Journal of Chemical Physics,1957,27:1208-1209. [15] 钱学森.物理学讲义[M].北京：科学出版社，1962:223-227. [16] 冯端,等.金属物理学第一卷：结构与缺陷[M].北京：科学出版社，
2000:31-76.
[17] Cotterill R M J,Doyama M.In: Hasiguti R eds.Energier and Atomic
Configuration of Line Defects and Plane Defects.Lattice Defects and their Interactions[M]. New York: Gordon and Breach Science Publishers Inc.,1976:62-75.
[18] Morse P M. Diatomic molecules according to the ware mechanics
Vibrational Levels[J] .Phys.Rev.,1929（34）.
[19] Johnson R A, Wilson W D. In: Gehlen P C, et al, eds. Defect Calculations
for Fcc and Bcc Metals. Interatomic Potentials and Simulation of Lattice Defects[M]. New York:Plenum,1971:301-305.
第40卷第5期 当 代 化 工 Vol.40，No.5 2011年5月 Contemporary Chemical Industry May，2011
分子模拟的方法及其应用
李春艳，刘 华，刘波涛
（中北大学物理系，山西 太原030001）
摘 要： 综述了分子模拟的方法及其在材料科学中的应用，介绍了分子模拟的基本原理、牛顿运动方程及其有限差分算法、势函数的发展、平衡态系综及其温度和压力的控制，最后还指出了分子模拟的应用及其进一步的研究方向。
关 键 词： 分子模拟; 有限差分法; 原子间相互作用势; 材料科学
中图分类号： O 642 文献标识码： A 文章编号： 1671-0460（2011）05-0494-03
Method of Molecular Simulation and Its Application
LI Chun-yan, LIU Hua, LIU Bo-tao
（North University of China, Shanxi Taiyuan 030001, China）
Abstract: Method of molecular simulation and its application were reviewed. The principle of molecular simulation, Newtonian equation of motion were introduced as well as its finite difference technique, progress of potential function, realization of equilibrium ensembles and the temperature and pressure control. Directions of researches and applications of molecular simulation were put forward.
Key words: Molecular simulation; Finite difference technique; Interatomic potential; Material science
在高分子领域中计算机模拟的应用从20世纪60年代至今已发展到了一个崭新的阶段。这个新阶段的特点是，计算机模拟不仅能提供定性的描述，而且能模拟出高分子材料结构与性能的定量的结果；模拟分子体系算法的发展，主要是从描述简单的非真实分子体系的算法到复杂的真实分子体系的算法。该发展克服了用蒙特卡洛法仅能够描述不同温度下分子结构的特征，但不能描述随温度变化而引起体系宏观物理性质的变化中分子结构的变化过程。
[1-2]
1.2 分子动力学方法(简称MD)的基本原理
分子动力学方法(MD)是另一种主要的计算机模拟方法，在材料科学、物理、化学等学科的各个领域得到广泛应用。模拟过程是在一定体系及分子势能函数条件下，从计算分子间作用力着手，求解牛顿运动方程，得到体系中各分子微观状态随时间的变化，再将分子的位置和动量组成的微观状态对时间平均，即可求出体系的压力、能量、粘度等宏观性质以及组成粒子的空间分布等微观结构。
[8]
1 分子模拟方法
1.1 蒙特卡洛方法(简称MC)的基本原理
MC法早在20世纪50年代就被应用到高分子科学的研究之中。先驱性的工作是Wall在为研究高分子链的排除体积问题时所做的蒙特卡洛模拟
[3-5]
2 分子模拟的计算
2.1 分子动力学的运动方程
分子动力学模拟的出发点是假定粒子的运动可以用经典动力学来处理, 对一个由N 个粒子构成的孤立体系, 粒子的运动由牛顿运动方程决定, 也就是
d 2r
m =-∇V (r , r , r ) （1）
i 12N
d t 式中：m i ，r i —第i 个原子的质量和位置；
▽i = -d/dri ,V (r 1, r 2, r N ) —体系所处的势。
。基
本模拟过程是在一定系统条件下，将系统内粒子进行随机的位移、转动，或粒子在两相间转移位置。根据给定的分子势能函数，进行粒子间内能的加和。采用Metropolis取样方法
[6-7]
，生成一系列体系的微观
粒子随机构型，从而逐渐趋近于平衡时的Boltzmann分布。目前蒙特卡洛的应用几乎涉及到高分子科学中的各个领域，为现代高分子科学理论基础的建立和发展起到了十分重要的推动作用。
2.2 分子动力学运动方程的求解
在分子动力学中，为了得到原子的运动，一般采用以下几种有限差分法来求解运动方程。常用的有：
基金项目：国家自然科学基金（批准号：50730009）。
作者简介：李春艳，女，硕士研究生，主要从事分子力学模拟研究。E-mail:[email protected]。
第40卷第5期 李春燕，等：分子模拟的方法及其应用 495
Verlet算法是在60年代后期出现的, 对扩散分子的质心运动的积分是最稳定的也是最常用的数值方法。Swope提出的Velocity–Verlet算法可以同时给出位置、速度与加速度，并且不牺牲精度。这种算法的优点是给出了显速度项，并且计算量适中，应用比较广泛。Hockney提出的“蛙跳”Leap–forg 算法是Verlet算法的变化，这种算法与Verlet算法相比有两个优点:（1）包括显速度项;（2）收敛速度快, 计算量小。这种算法明显的缺陷是位置和速度不同步。Beeman提出的Beeman算法运用了更精确的速度表达式。它的表示式比Verlet算法复杂得多，因此计算量较大。Gear的预测——校正算法分为三步来完成:首先, 根据Taylor展开, 预测新的位置、速度和加速度。然后, 根据新的计算的力计算加速度。这个加
[13]
[12]
[11]
[10]
[9]
模型主要有EAM势、FS势、TB势等。 3.2.1 EAM势
Daw和Baskes在1984年首次提出了嵌入原子法（EAM）。EAM势的基本思想是把晶体的总势能分成两部分：一部分是晶格点阵上的原子核之间的相互作用对势，另一部分是原子核镶嵌在电子云中的嵌入能，它代表多体相互作用。在嵌入原子法中，系统的总势能表示为：
E =F ρ+1φr
（6）∑i (i )∑ij (i , j )tot
i
[20]
2
ij
式中：F 是嵌入能：第2项是对势项，根据需
要可以取不同的形式。ρi 是除第i 个原子以外的所有其它原子的核外电子在第i 个原子处产生的电子云
ρi (r ij )，密度之和，可以表示为：ρi =∑其中，ρi （r ij ）j ≠i
是第j 个原子的核外电子在第i 个原子处贡献的电荷
速度再由与Taylor级数展开式中的加速度进行比较。
密度，r ij 是第i 个原子与第j 个原子之间的距离。对于
这种方法的缺点就是占有计算机的内存大。
不同的金属，嵌入能函数和对势函数需要通过拟合金属的宏观参数来确定。 3 分子间的作用势 3.1 对势
在分子动力学模拟的初期，人们经常采用的是对势模型有以下几种：
（1）在1957年Alder和Wainwright使用间断对势: 当r ≤r 1时,E =∞;当r ＞r 1时, E =0。
（2）连续对势
Lennard-Jones(L-J )势是最老的原子间作用势之一，是在1924年由J.E.Lennard-Jone提出来的，
[15]
[14]
3.2.2 FS势
Finnis等根据金属能带的紧束缚理论，形成了
[21]
一种在数学上等同于EAM的势函数，并给出了多体势的函数形式，把嵌入能函数设为平方根形式，即
ρi (r ij )项分别为：FS势。其中的多体项及对势ρi =∑j ≠i
ρ(r ij )=∑A k (R k −r ij )H (R K −r ij )
3
k =1
6
2
3
σij (r ij )=∑a k (r k −r ij )H (r k −r ij ) （7）
k =1
L-J 势的基本形式为：
φ(R )=4ε⎡⎛σ⎞−⎛σ⎞⎤ （2）
⎢⎜⎟⎜⎟⎥LJ
⎢⎝R ⎠⎝R ⎠⎥⎣⎦
式中：σ、ε—经验常数。
12
6
式中：当x >0时，H (x )=0；当x
基于EAM势的势函数还有很多种的非球形对称。于是Baskes等
[24-26]
[22-23]
（3）Born-Lande(B-J )势的形式如下： 2
Z i Z j b
φij (r ij )=e +m （3）
r ij 4πεr ij （4）Morse势：
根据双原子分子的振动谱，提出了指数形式的相互作用势：
φ(r )=Ae (−αr )−Be (−βr ) （4）
它有4个参数A , B , α和β，与L-J 势的普适形式相类似。
（5）Johnson势的形式如下：
φij (r ij )=−A n (r ij −B n )+C n r ij −D n （5）
3[19][18]
[17]
[16]
。为了把
将EAM势推广到共价键材料中去，需考虑到电子云
提出了修正型嵌入
原子核法(MEAM)。总的电子密度为下式：
3
(ρ)2=t (h )(ρ(h ))2 （8）
i
∑
h =0
i i
给出。其中t (h ) i 是常系数。在MEAM里，镶嵌能项采取了如下形式：F (ρ)=AEρln ρ,其中A ，E 是系数，因物质不同而不同。
3.2 多体势
多体势于20世纪80年代初期开始出现。多体势于20世纪80年代初期开始出现。常用的多体势函数
4 分子动力学模拟中平衡系综的控制
方法
平衡态分子动力学模拟是在一定的系综下进行的, 经常用的平衡系综是NVT或NPT系综。在这两种系综中, 牵涉到控制温度和压力的几种技术,分别介
496 当 代 化 工 2011年5月
绍如下： 4.1 控温方法
在平衡系综中，经常要把温度调整到期望值。而体系的有效温度是由动能决定的：
k =∑m i v i / 2=fk B T /2
式中：f 为体系的自由度，K B 为波尔兹曼常数。 控温常用的方法有: （1）速度标定法
系统的温度与动能存在如下的关系：
E k =∑m i v i 2/2=(3N −N C )K B T /2 （9）
i N
[27]
2
与他的控温方法类似, 如下:
μ=[1+△t (P -P 0) γ/Г]1/3
式中γ是一个可调参数, P 0为所控制的压力。然
后, 在每一次的分子动力学模拟迭代中, 粒子坐标
x , y , z 均用μ相乘, 得到新的坐标, 即可实现对压力P 的控制。
4.2.2 Andersen方法
这种方法是通过改变体系的Lagrangian来实现的, 新的Lagrangian定义如下：
i i N N ii i
L =m s G h /2−V (r )+WT h ' h /2−p Ω（12）
[32]
∑
i i
∑∑
i =1j >1
ij r
式中：N —原子数；
N C —约束数；
K B —波尔兹曼常数； v i —原子i 的速度。
其中：h =﹛a , b , c ﹜—模拟单胞的基矢；
G =h ’h ,原子坐标r i 为r i =hs ； V —原子之间的相互作用势； Ω—单胞体积； T r — 矩阵的迹；
这种方法的基本思想是：如果t 时刻的温度是T （t ）, 速度乘以因子λ后, 温度的变化为：
2
ΔT =（λ-1）T （t ）
T 为所控制体系的温度。 （2） Berendsen热浴法法
这种也称恒温槽方法，是假设所模拟的体系与一个恒温槽连在一起, 那么两者之间就可以通过热交换而使所模拟的体系达到恒温目的, 方法如下:
λ=,式中,参数Г表征系统与恒温槽
[28]
W — 一个可调参量这种方法的优点是模拟单胞不仅大小可变,
形状也可变。
5 分子模拟在材料科学中的应用
分子动力学模拟计算液态金属的结构及热力学性质在国外得到了广泛地应用。其中波兰的Janusz 等利用多体相互作用，借助分子动力学计算了Ag、Au、Cu、Ni等面心立方金属的热力学性质。Holzman 等采用EAM 势计算了面心立方金属液-气界面的特殊性，得到的密度、内能、结构因子等结果与实验值相一致。N.Tajima 等采用EAM 势模拟了面心立方金属的点缺陷。K.Kadau等运用分子动力学模拟了烧结Fe-Ni 纳米颗粒的马氏体转变。Tsai 采用第一性原理的分子动力学方法，研究了GaN 表面气态Ga 和N 原子的吸附特征，模拟结果显示：引入的Ga 原子相对不活泼，而引入的N 原子却很活泼；Ga在GaN 表面的吸附系数要比N 原子高很多。
[37]
[36]
[35]
[34]
[33]
之间的热交换速率,Δt 为分子动力学模拟的时间步
长,通过v new =λv old 校正即可保持体系的温度在T 0附近振动, 而参数Г可用于控制这个振动幅度。
（3）Nose-Hoover方法法
[29-30]
Nose-Hoover热浴法是通过改变模拟体系的哈密顿量来实现控温的,其基本思想就是在哈密顿量里加入一个假想的项来代表一个恒温源, 具体做法如下：
22
H =∑m i v i /2+V (r )+Q ζ/2+gKT ln S （10） 式中：S 和ζ分别是假想项的坐标和动量。
近几年，在研究液态金属方面，中国科学院冶这样体系的微分方程就变为：
[38]
/r +m i v i z )/m i ， 金研究所陈柳等对晶体的生长过程进行了分子动v i =d r i /d r ，a i =-(d V /d
[39]
力学模拟研究。王金照等利用分子动力学方法模2
（11）
/t =(∑m i v i -gK B T )/Q d z /d
式中：g — 体系的自由度；
Q — 一个可调参量, 表征着假想项的质量, T 为温度。
拟研究了铜铝合金的比热。戎咏华等对铁锰合金层错位能进行了计算。张建民等采用改进嵌入原
[41]
[42]
[40]
4.2 压力控制
分子动力学模拟通常要求模拟体系的压力保持不变，而体积可以随温度的变化而改变，即进行等压模拟。这可以通过改变模拟原胞的3个方向或一个方向的尺寸来实现体积的变化。通常有： 4.2.1 Berendsen方法
这种控压方法的基本思路与基本控制方程均
[31]
子法计算了铜晶体的表面能。张朝阳等采用
Discover／Materials Studio 中的MD（NVT）及MM （COMPASS力场）方法，通过对两种氟聚合物在TATB 不同晶面上吸附行为的动力学模拟和结构优化发现：二者相互作用时,TATB原子位置变化的大小次序为：（010）>（100）>（001）；而氟聚合物的展开程度按以上次序增加。宋小艳等利用三
[43]
第40卷第5期 李春燕，等：分子模拟的方法及其应用 497
维MC 技术模拟了较完整的单相材料正常晶粒长大过程，得到了晶粒长大动力学和拓扑学的全面信息。莫春立等以单相铁素体不锈钢为研究对象，利用基于实验数据的模型，对焊接影响区晶粒长大过程进行了MC 动态模拟。清华大学的郭靖原等运用MC 方法模拟了金属凝固与晶体生长过程。孙雷等
[47]
[46]
[45]
[44]
[20] Daw M S,Baskes M I. Embedded atom method derivation and application
to impurities, surfaces, and other defects in metals[J]. Physical Review B, 1984,29（12）: 8486-8495.
[21] Finnis M W，Sinclair J E．A simple empirical n —body potential for
transition—metal[J]．Philosophic Magazine A，1984，50（1）：45. [22] Foiles S M, Baskes M I,Daw M S. Embedded-atmo-method function for
the fcc metals Cu, Ag, Au,Ni, Pd, Pt and their slloys[J]. Physical Review B, 1986, 33（12）: 7983-7991.
[23] Manninen M,Johnson R A. Interatomic interactions in solids[J]. Physical
Review B,1986, 34（12）:8486-8495.
[24] Baskes M I. Modified embedded-atom potentials for cubic materials and
impurities[J]. Physical Review B, 1992, 46（5）:2727-2742. [25] Baskes M I, Johnson R A.Modified embedded-atom potentials for hcp
metals. Modelling Simul Mater Sci Eng[J]. Phys.Rev.1994,2: 147-163. [26] Baskes M I.Calculation of the behaviour of Si ad-dimers on
Si(001).Modelling Simul Mater Sci Eng[J]. Phys.Rev., 1997,5: 149-158. [27] Hofh曲nn K H，Schreiber~L C0mputat．oml Physics．Berlm Heidelberg：
Spring—verlag[J].Phys.Rev.,1996:268.
[28] Berendsen H J C，Postma J P M，Gunsteren W F V，et a1.[J]．chdTl Phys，
1984，81：3684.
[29] N05e S J.[J] Ch1 Phys，1984，8l：5. [30] H00ver W G．[J].Ph”RA，1985，3l：1695.
[31] Berendsen H J C，Postma J P M，Gunsteren W F V，et a1．[J].J chdTl
Phys，1984，81：3684.
[32] AnderseIl H C.Chem Phys[J]．Journal of Computational
Physics,1980:721-2384.
[33] Anuse J,Holender L M.[J]. Phys Rev B, 1990,141（12）: 8054. [34] Holzman L M,Adams J B,et al.[J]. Mater Res,1991,6（2）:343.
[35] Tajima N,Takai O,Kogure Y, et al. Computer simulation of point defects in
fcc metals using EAM potentials[J]. Computational Materials Science, 1999, 14:152
[36] Kadau K, Entel P, Lomdahl P S. Molecular-dynamics study of martensitic
transformations in sintered Fe-Ni nanoparticles[J]. Computer Physics Communications, 2002, 147: 126.
[37] Tsai M H．[J].Computer Physics,2002，147：130.
[38] 陈柳，程兆年，唐鼎元.晶体生长过程的分子动力学模拟研究[J].人
体晶体学报，1999，28（2）：109.
[39] 王金照，杨春，陈民,等.铜银合金比热的分子动力学模拟[J].工程热
物理报，2002,23（3）：274.
[40] 戎咏华，孟庆平，何刚,等.Fe-Mn 合金层错能的嵌入原子法计算[J].
上海交通大学学报，2003,37（2）：171.
[41] 张建民，徐可为，马飞.用改进嵌入原子法计算Cu 晶体的表面能[J].
半导体学报,2003,52（8）：1993.
[42] 张朝阳，舒远杰，黄辉，等．爆炸与冲击[J].金属学报，2004，24
（6）：563.
[43] 宋小艳．基于计算机模拟的正常晶粒长大的三维特征二维截面表征
[J]．材料研究学报，1998，12（3）：245
[44] 莫春立．铁素体不锈钢焊接热影响区晶粒长大过程的模拟[J]．金属
学报，2001，37（3）：307.
[45] 郭靖原，洪泽恺，陈民，等．金属凝固与晶体生长过程的蒙特卡罗
模拟[J]．工程热物理学报，2001，22（6）：725.
[46] 孙雷，杜刚，刘小彦，等．蒙特卡罗方法模拟金属一半导体基础的
直接隧穿效应（英文）[J]．半导体学报，2001，22（11）：1364. [47] 姜寿文，赵国群，关小军，等．冷轧薄钢板退火过程组织演变的
Monte Carlo模拟[J]．钢铁研究学报，2002，14（5）：53.
模拟了金属-半导体接触的直接隧穿效应。姜寿文等对冷轧薄钢板退火过程的组织演变进行了MC 模拟。
参考文献：
[1] E.B.Wilson, J.C.Decius, P.C.Cross. Molecular Vibrations, McGraw-Hill[M].
New York:1995:1-51.
[2] T.Shimanochi.Vibration Spectroscopy and Its Chemical application[M].
Tokyo:Univ.of Tokyo, 1997:105-406.
[3] F.T.Wall.Strain fluctuations and elastic constants [J].Chem.Phys, 1953,21
（19142）:98-235.
[4] F.T.Wall,L.A.Hiller, Jr., D.J.Wheeler. A new molecular dynamics
method[J].Chem. Phys., 1954,22（1036）:76-343.
[5] F.T.Wall,et al.. Interatomic interactions in solids[J].Chem.Phys.,1955:
23-231.
[6] Frenkel D，Smit B．Understanding molecular simulation：From algorithms
to applications[M]．SanDiego：Academic Press，1996:21-196. [7] 分子模拟：从算法到应用[M]．汪文川，等译．北京：化学工业出版
社，2002.
[8] Heermann D W.ComDuter simulation methods in theoretical physics.Berlin
Heidelberg：Springer—Verlag[J].Phys.Rev.,1990:81-345.
[9] L.Verlet.Computer ‘experiments’ classical fluids.I.Thermodynamical
properties of Lennard-Jonesmolecules[J].Phys.Rev.,1967（159）:98-103. [10] 孙伟，常明，杨保和.分子动力学模拟纳米晶体铜的结构与性质[J].
物理学报，1998,47（4）：591-597. [11]
R.W.Honeycutt.Thepotential
calculation
and
some
applications[J].Mcthods in Computational Physics,1970,9:136-211. [12] Beeman D.Some multistep methods for use in molecular dynamics
calculations.Journal of Computational Physics[J]. Phys.Rev.,1976,20: 130-139.
[13] Gear CW.Numerical imitial value problems in ordinary differential
equation[M]. Englewood Cliffs:Prentice--Hall,1971.
[14] Alder B J, Wainwright T E. Phase transition for a hardsphere
system[J] .The Journal of Chemical Physics,1957,27:1208-1209. [15] 钱学森.物理学讲义[M].北京：科学出版社，1962:223-227. [16] 冯端,等.金属物理学第一卷：结构与缺陷[M].北京：科学出版社，
2000:31-76.
[17] Cotterill R M J,Doyama M.In: Hasiguti R eds.Energier and Atomic
Configuration of Line Defects and Plane Defects.Lattice Defects and their Interactions[M]. New York: Gordon and Breach Science Publishers Inc.,1976:62-75.
[18] Morse P M. Diatomic molecules according to the ware mechanics
Vibrational Levels[J] .Phys.Rev.,1929（34）.
[19] Johnson R A, Wilson W D. In: Gehlen P C, et al, eds. Defect Calculations
for Fcc and Bcc Metals. Interatomic Potentials and Simulation of Lattice Defects[M]. New York:Plenum,1971:301-305.