-
分子动力学模拟及其在金属材料中的研究进展
摘要
本文综述了分子动力学模拟技术
的发展，介绍了分子动力学的分类、运动方程的求解、
初始条件和边界条件的选取、
平衡系综及其控制、
感兴趣量的提取以及分子动力学模拟在金
属材料中的研究进展。
关键词
：分子动力学模拟
平衡态系综
金属材料
感兴趣量
径向分布函数
引言
科学工作者在长期的科学研究实
践中发现，当实验研究方法不能满足研究工作的需求
时，
用计算
机模拟却可以提供实验上尚无法获得或很难获得的重要信息；
尽管计算机模拟不
能完全取代实验，
但可以用来指导实验，
并验证
某些理论假设，
从而促进理论和实验的发展。
特别是在材料形成
过程中许多与原子有关的微观细节，
在实验中基本上是无法获得的，
而在
计算机模拟中即可以方便地得到。
这种优点使分子动力
学模拟在金属材料研究中显得非常有
吸引力。
分子动力学
MD (Molecular Dynamics)
模拟就是用计算机方法来表示统计力学，作为实
验的一个辅助手
段。
MD
模拟就是对于原子核和电子所构成的多体系统，求解运
动方程
(
如牛
顿方程、哈密顿方程或拉
格朗日方程
)
，其中每一个原子核被视为在全部其它原子核和电
子
作用下运动，
通过分析系统中各粒子的受力情况，
用经典或量子的方法求解系统中各粒子在
某时刻的位置和速度，
以确定粒子的运动状态，
进而计算系统的结构和性质。
该模拟技术主
要涉及粒子运动的动力学问题
,
与蒙特卡罗模拟方法
(
简称
MC)
相比
,
分子动力学是一种
“确定
性方法”
,
它
所计算的是时间平均
,
而
MC
进行的是系综平均。
然而按照统计力学各态历经假设
,
时间平均等价于系综平均。因此
,
两种方法严格的比较计算能给出几乎相同的结果。
经典的分
子动力学方法是
Alder
等于
195
7
年提出并首先在“硬球”液体模型下应用，发
现了由
Kirkwood
在
1939
< br>年根据统计力学预言的
“刚性球组成的集合系统会发生有液相到结晶
相的转变”。后来人们称这种相变为
Alder
相变
。
Rahman
于
1963
年采用连续势模型研究了液
体的分子动力学模拟。
1972
年
Less
等发展了该方法并
扩展了存在速度梯度的非平衡系统。
1980
年
Andersen
等创造了恒压分子动力学方法。
19
83
年
Gillan
等将该方法推广到
具有温度梯
度的非平衡系统，
从而形成了非平衡系统分子动力学
方法体系。
1984
年
Nose
等完成了恒温分
子动力学方法的创建。
198
5
年针对势函数模型化比较困难的半导体和金属等，
Car
等提出了
将电子论与分子动力学方法有机统一起来的第一性原理分子
动力学方法。
1991
年
Cagin<
/p>
等
[1]
进一步提出了应用于处理吸附问
题的巨正则系综分子动力学方法。
20
世纪
80
年代后期，
计算
机技术飞速发
展，
加上多体势函数的提出与发展，
使分子动力学模拟技术有了
进一步的发展。
1
分子动力学分类
分子动力学的目标是
研究体系中与时间和温度等有关的性质而不只是静力学模拟中研
究的构型方面。
分子动力学假定原子的运动是由牛顿运动方程决定的，
这意味着原子的运
动
第
1
页
共
9
页
是与特定的轨道联系在一起的。
分子动力学模拟的关键问题是原
子间作用势的确定，
主要是
求解下述牛顿运动方程组。
<
/p>
其中
M
a
为原子
质量，
R
a
为原子空间位置，
t
表示时间，
F
α
为原子间作用力。
确定原子间的
相互作用力
F
α
，也就是确定原子间作用势
E(R
)
。
确定原子间作用势，
必须知道相应的电子基态。
电子基态的计算是一个非常复杂的量子
多体问题，即解多体薛定愕
(Schrodinger)
方程
(
式
(2)):
式中：
E
tot
表示系统的总能量，
r
i
表示第
i
电子的空间坐标，
Ψ
(r
i
，
R
α
)
是系统波函数。系统
的哈密顿算子
(Hamilt
onian)H
可表示为
:
其中
P
、
p
分别表示核和电子的动量算子，
M
、
m
分别表示核和电子的质量，
α
< br>、
β
表示原
子核的序号，
Z
表示电荷数。式
(3)
右端第一、二项分别表示核和电子的动能，
第三项表示电
子
间的相互作用势，
第四项表示核和电子的相互作用势，
第五项表
示核间的相互作用势。
根
据
Born-
Oppenheimer
近似
(
电子云
结构受核运动的影响极小
)
，
系统的薛
定愕方程可分离为原
子核薛定愕方程和电子薛定愕方程。而电子薛定愕方程进一步可写为
:
其中
Ψ
(r
i
，
R<
/p>
α
)
为电子的波函数，
< br>E(R
α
)
的物理意义是核静止
时系统的基态能量，是核
坐标
R
。的函
数，可以理解为原子间作用势。当
F
α
确定时，就可以通过求解牛顿运动方程分
析系统的力学行为。
事实上，
求解薛定愕方程是非常困难的，
因此通常是通过试验拟合或半经验解法得到原
子间作用势，然后求得系统能量。也就是
说，
分子动力学模拟通常是经验或半经验的。
根据
对原子间作用势不同的简化处理方法，
分子动力学可划分为经典分子动力学和
现代分子动力
学。
N
（
1
）经典分子动力学
经典分子动力学
(ClassicalMD)
通过实验结果或经验模型确定原子间作用势，计算量较
小，可以解决较大规模的问题，
但是可移植性
(Transferability)
差。针对不
同的问题，可
能需要确定不同的经验参数。在
20
世纪
80
年代以前，分子动力学模拟一般都采用对势
模型
(Pairpotential)
，该模型仅考虑近邻原子
间的库仑作用力和短程相互作用，并认为系统能
量为各粒子能量总和。对势可以比较好地
描述除金属和半导体以外的几乎所有无机化合物。
比较常用的对势有硬球势、
Lennard-Jones(LJ)
势、
Mor
se
势、
Johnson
势等，它们在
特定的
问题中均有各自的优越性。
实
际上，
在多原子体系中一个原子的位置不同将影响空间一定范围内的电子云分布，
从
而影响其他原子之间的有效相互作用，因此，人们开始考虑粒子间的
多体作用
第
2
页
共
9
页
(Many-
bodyeffects)
，构造出多体势结构。
多体势于
20
世纪
80<
/p>
年代初期开始出现，
Daw
等在
1984
年首次提出了嵌人原子法
(Embed
ded-atommethod,EAM),EAM
势很好地描述了金属原子之间的相互
作用，是描述金属体
系最常用的一种势函数。
对于由共价键结合
的有机分子以及半导体材料并不适用。
为更好描
述各种含有共价
键作用的物质，人们考虑了电子云的非球形对称，将
EAM
势推
广到共价健材
料。
为此，
Baskes
等提出了修正嵌人原子核法
(MEAM)
。
从某种意义上说这个模型是半经验的，
因为它从局域电子密
度观点出发解决全部问题，使用的参数有从实验中获得的数据
(
如晶格
常数、转变能、体积弹性模量、弹性系数等
)
。
（
2
）现代分子动力学
为了克服传统分子动力学可移植性
差这一缺陷，人们考虑直接从量子力学
(Quantummechanics,QM)<
/p>
轨道理论出发获取原子间作用势。基于
QM
的分子动力学称之为现代
分子动力学
(ModernMD)<
/p>
，
也称之为从头分子动力学
(Abini
tiomoleculardynamics,AIMD)
。
密
度泛函分子动力学
(DFMD)
和第一
原理分子动力学
(FPMD)
是比较常用的。
< br>DFMD
是基于量子力学
密度泛函理论
< br>(Densityfunctionaltheory)
，直接从量子力学基本原
理考虑电子云结构，模
拟更为准确，可移植性更好，但计算量大。密度泛函理论是在量子
理论基础上建立起来的，
从波函数出发定义电子的密度，赋予波函数确切的物理意义，通
过求解
Schrodinger
方程，
确定电子的密度，再根据能量与密度的关系给出系统的能量。
第一原理分子动力学
(First-principlesmo
leculardynamics,FPMD)
是利用第一原理法
对电子结构进行计算，
解决材料中各元素间的成键、
结合和相
稳定性，
材料的力学行为与电
子结构和成键性质、
电荷分布的主要方向等。
Smargiassi
等给
出了一种
FPMD
算法，
它把系统
总能量进一步细分为
8
个部分，
并利用各部分现成的显式结果，
只有凡由
LDA<
/p>
得到。
这种方法
最关键的一点是计算中不
需要波函数的显式结构，
使其计算量大大减小。
对于不同的应用
对
象原子间势函数的势参数互不相同，势参数的确定一般有
3<
/p>
种方法
:
①通过实验值
< br>(
如晶格常
数、弹性常数、内聚能和空位形成能等
)
拟合势参数；②通过蒙特卡罗方法确定势参数；③
< br>通过基于量子力学得到的各种微观信息来确定势参数。
2
运动方程的求解
在分子动力学中，
系统中原子的一系列位形是通过牛顿运动方程积分得到的。<
/p>
为了得到
原子的运动，一般采用各种有限差分法来求解运动方程。
常用的几种算法如下
:
(1)Verlet
< br>算法是一种目前应用最广泛的数值积分求解运动方程组的算法。由系统的哈
密顿量
可以推导出牛顿方程形式的运动方程组
:
令
则可以得到如下差分方程组的形式
:
第
3
页
共
9
页
由空间坐标可以算出粒子的运动速度为
:
令
则
(8)
式可以写为
:
利用上式则可以在计算中得到同一时间步上的空间位
置和速度，
并且提高了数值计算的稳定
性。
(2)Swope
提出的
Vel
ocity-Verlet
算法可以同时给出位置、速度与加速度，并且不牺牲
精度。这种算法的优点是给出了显式速度项，并且计算量适中，目前应用比较广泛。
3
分子动力学模拟的初始条件和边界条件
系统的初始位形和初始速度可通过实验数据、
理论模型或两者的结合来确定。
如果模拟
系统无固定晶格结构，则每个原子的位置可用舍选法或
Petropolis
等方法从初始密度分布
n
(r)
得到；
每个原子的初速度可从初始温度分布
T
c
(r)
下的
Maxwell-Boltzmann
分布来随机选
取。
对于流体分子系统，
粒子的初始速度按高斯分布来选取比
较合理。
物体的宏观性质是大
量粒子的统计行为，
模拟系统的粒子数必须非常大才能准确地再现系统的行为。
分子动力学
模拟采用周期性边界条件的假定。
取模拟系统为中心元胞，
一般为立方体，
体系中含有几百
至上万个粒子，<
/p>
使中心元胞在三维空间上重复排列，
于是系统粒子￡的象粒子将在
三维空间
周期性出现，
体系势能包括粒子与象粒子间的相互作用
在内。
应用周期性边界条件后，
整个
体
系就变得粒子数赝无限了。
周期性边界条件是一个近似，
它给体系强加一个实际并不存在的长程序，
如果是带电粒
子系统还会人为造成极化。显然元胞所含粒子数
N
不能太小，否则将因边界条件的影响而极
大地偏离实际情况。
4
分子动力学模拟的系综
在分子动力学模拟方法的应用中，
存在着对两种系统状态的
MD
模拟。
一种是对平衡态的
MD
模拟，
另一种是对非平衡态的
M
D
模拟，
两者的差别在于平衡系统的感兴趣量及边界条件与
时间无关。
为了准确反映实际的物理过程，
分子动力学模拟总是在一定的系综下进行，
具体
有微正则系综，
即体系的原子数
(N)
、体积
(V)<
/p>
和能量
(E)
都保持不变；正则系综，即
体系的
原子数
(N)
、体积
(V)
和温度
(T)
都保
持不变，并且总动量为零
(P=0)
；
等温等压系综，
即系统
的原子数
(N)
、压力
(P)
和温度
< br>(T)
都保持不变；
等压等烩系综，
即保持系统的原子数
(N)
、压
力
(P)
和熔值
(H)
< br>都不变。
5
感兴趣量的提取
求出模拟系统中原子
在相空间的轨迹后，
应用统计物理原理即可得到系统的宏观特性和
结构特点。在此讨论几种常用物理参量的统计形式。
第
4
页
共
9
页
-
-
-
-
-
-
-
-
-
上一篇：全球SCI收录材料期刊影响因子排名
下一篇：美语里最常见的10个问句