分子模拟的方法及其应用
合集下载
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
式中：f为体系的自由度，KB为波尔兹曼常数。 控温常用的方法有: （1）速度标定法[27]
系统的温度与动能存在如下的关系：
N
∑ Ek = mi vi 2 / 2 = (3N − NC ) KBT / 2 （9）
i
式中：N —原子数； NC —约束数； KB —波尔兹曼常数； vi —原子i的速度。
这种方法的基本思想是：如果t 时刻的温度是T
取了如下形式：F(ρ)=AEρlnρ,其中A，E是系数，
因物质不同而不同。
（5）Johnson势的形式如下[19]：
( ) ( ) φij rij = − An rij − Bn 3 + Cnrij − Dn
（5）
3.2 多体势
多体势于20世纪80年代初期开始出现。多体势
于20世纪80年代初期开始出现。常用的多体势函数
模拟过程是在一定体系及分子势能函数条件下从计算分子间作用力着手求解牛顿运动方程得到体系中各分子微观状态随时间的变化再将分子的位置和动量组成的微观状态对时间平均即可求出体系的压力能量粘度等宏观性质以及组成粒子的空间分布等微观结构8
第 40 卷第 5 期 2011 年 5 月
当代化工 Contemporary hemical Industry
当r≤r1时,E=∞;当r＞r1时, E=0。 （2）连续对势
Lennard-Jones(L-J )势是最老的原子间作用势
之一，是在1924年由J.E.Lennard-Jone提出来的[15]，
L-J 势的基本形式为：
φLJ
(
R)
=
4ε
⎡⎢⎢⎣⎛⎜⎝
σ R
⎞12 ⎟⎠
−
⎛ ⎜⎝
σ R
⎞6 ⎟⎠
⎤ ⎥ ⎥⎦
在高分子领域中计算机模拟的应用从 20 世纪 60 年代至今已发展到了一个崭新的阶段[1-2]。这个新阶段 的特点是，计算机模拟不仅能提供定性的描述，而且 能模拟出高分子材料结构与性能的定量的结果；模拟 分子体系算法的发展，主要是从描述简单的非真实分 子体系的算法到复杂的真实分子体系的算法。该发展 克服了用蒙特卡洛法仅能够描述不同温度下分子结 构的特征，但不能描述随温度变化而引起体系宏观物 理性质的变化中分子结构的变化过程。
比有两个优点:（1）包括显速度项;（2）收敛速度快, 系统的总势能表示为：
计算量小。这种算法明显的缺陷是位置和速度不同 步。Beeman提出的Beeman算法[12]运用了更精确的速度 表达式。它的表示式比Verlet算法复杂得多，因此计 算量较大。Gear的预测——校正算法[13]分为三步来完 成:首先, 根据Taylor展开, 预测新的位置、速度和加 速度。然后, 根据新的计算的力计算加速度。这个加 速度再由与Taylor级数展开式中的加速度进行比较。 这种方法的缺点就是占有计算机的内存大。
∑ ∑ ( ) ( ) Etot =
i
Fi
ρi
+1 2
φij
ij
ri, j
（6）
式中：F 是嵌入能：第2项是对势项，根据需
要可以取不同的形式。ρi是除第i个原子以外的所有
其它原子的核外电子在第i个原子处产生的电子云
密度之和，可以表示为：ρi = ∑ ρi (rij ) ，其中，ρi（rij） j≠i
Vol.40，No.5 May，2011
分子模拟的方法及其应用
李春艳，刘 华，刘波涛
（中北大学物理系，山西 太原 030001）
摘
要： 综述了分子模拟的方法及其在材料科学中的应用，介绍了分子模拟的基本原理、牛顿运动方程
及其有限差分算法、势函数的发展、平衡态系综及其温度和压力的控制，最后还指出了分子模拟的应用及其进
（4）
它有4个参数A,B,α和β，与L-J 势的普适形式 相类似。
3.2.2 FS势 Finnis等根据金属能带的紧束缚理论，形成了
一种在数学上等同于EAM的势函数[21]，并给出了多体
势的函数形式，把嵌入能函数设为平方根形式，即
FS势。其中的多体项及对势 ρi = ∑ ρi (rij ) 项分别为： j≠i
（t）, 速度乘以因子λ后, 温度的变化为： ΔT=（λ2-1）T（t）
T为所控制体系的温度。
（2） Berendsen热浴法法[28] 这种也称恒温槽方法，是假设所模拟的体系与
一个恒温槽连在一起, 那么两者之间就可以通过热
交换而使所模拟的体系达到恒温目的, 方法如下:
λ = 1+ ∆t (T −TC ) /τ ,式中,参数Г表征系统与恒温槽
给出位置、速度与加速度，并且不牺牲精度。这种 （EAM）[20]。EAM势的基本思想是把晶体的总势能
算法的优点是给出了显速度项，并且计算量适中， 分成两部分：一部分是晶格点阵上的原子核之间的
应用比较广泛。Hockney提出的“蛙跳”Leap–forg 相互作用对势，另一部分是原子核镶嵌在电子云中
算法[11]是Verlet算法的变化，这种算法与Verlet算法相 的嵌入能，它代表多体相互作用。在嵌入原子法中，
2.2 分子动力学运动方程的求解
（1）
在分子动力学中，为了得到原子的运动，一般采
用以下几种有限差分法来求解运动方程。常用的有：
基金项目：国家自然科学基金（批准号：50730009）。 作者简介：李春艳，女，硕士研究生，主要从事分子力学模拟研究。E-mail:lcy104104@。
第 40 卷第 5 期
（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
3.2.3 TB势
基于EAM势的势函数还有很多种[22-23]。为了把
将EAM势推广到共价键材料中去，需考虑到电子云
的非球形对称。于是Baskes等[24-26]提出了修正型嵌入
原子核法(MEAM)。总的电子密度为下式：
∑ ( ) ( )ρi
2=
3
ti(h)
ρi(h)
2
h=0
（8）
给出。其中t(h)i是常系数。在MEAM里，镶嵌能项采
李春燕，等：分子模拟的方法及其应用
495
Verlet算法[9]是在60年代后期出现的, 对扩散分 模型主要有EAM势、FS势、TB势等。
子的质心运动的积分是最稳定的也是最常用的数值 3.2.1 EAM势
方法。Swope提出的Velocity–Verlet算法[10]可以同时
Daw和Baskes在1984年首次提出了嵌入原子法
2.1 分子动力学的运动方程
分子动力学模拟的出发点是假定粒子的运动
可以用经典动力学来处理, 对一个由N个粒子构成
的孤立体系, 粒子的运动由牛顿运动方程决定, 也
就是
m
d2ri dt 2
=-∇iV (r1,r2,rN
)
式中：mi，ri —第i个原子的质量和位置；
▽i = -d/dri,V(r1,r2,rN)—体系所处的势。
Q — 一个可调参量, 表征着假想项的质量, T为温度。
4.2 压力控制
分子动力学模拟通常要求模拟体系的压力保
持不变，而体积可以随温度的变化而改变，即进行
等压模拟。这可以通过改变模拟原胞的3个方向或一
个方向的尺寸来实现体积的变化。通常有： 4.2.1 Berendsen方法[31]
这种控压方法的基本思路与基本控制方程均
1.2 分子动力学方法(简称MD)的基本原理 分子动力学方法(MD)是另一种主要的计算机模
拟方法，在材料科学、物理、化学等学科的各个领 域得到广泛应用。模拟过程是在一定体系及分子势 能函数条件下，从计算分子间作用力着手，求解牛 顿运动方程，得到体系中各分子微观状态随时间的 变化，再将分子的位置和动量组成的微观状态对时 间平均，即可求出体系的压力、能量、粘度等宏观 性质以及组成粒子的空间分布等微观结构[8]。
1 分子模拟方法
2 分子模拟的计算
1.1 蒙特卡洛方法(简称MC)的基本原理 MC法早在20世纪50年代就被应用到高分子科
学的研究之中。先驱性的工作是Wall在为研究高分 子链的排除体积问题时所做的蒙特卡洛模拟[3-5]。基 本模拟过程是在一定系统条件下，将系统内粒子进 行随机的位移、转动，或粒子在两相间转移位置。 根据给定的分子势能函数，进行粒子间内能的加和。 采用Metropolis取样方法[6-7]，生成一系列体系的微观 粒子随机构型，从而逐渐趋近于平衡时的Boltzmann 分布。目前蒙特卡洛的应用几乎涉及到高分子科 学中的各个领域，为现代高分子科学理论基础的 建立和发展起到了十分重要的推动作用。
( ) ∑ ( ) ( ) ρ
rij
=
2
Ak
Rk − rij
3
H
RK − rij
k =1
( ) ∑ ( ) ( ) σij
rij
=
6
ak
rk − rij
3
H
rk − rij
k =1
（7）
式中：当x>0时，H (x)=0；当x<0时，H(x)=1。
Ak、RK、aK、rK为常数，且有R1＞R2,r1＞r2＞…＞r6。 它们的值随金属物质的不同而有所不同。
一步的研究方向。
关 键 词： 分子模拟; 有限差分法; 原子间相互作用势; 材料科学
中图分类号： O 642
文献标识码： A 文章编号： 1671-0460（2011）05-0494-03
Method of Molecular Simulation and Its Application
LI Chun-yan, LIU Hua, LIU Bo-tao
如下： H=∑miv2i/2+V(r)+Qζ2/2+gKT lnS
式中：S和ζ分别是假想项的坐标和动量。
（10）
这样体系的微分方程就变为：
vi = dri / dr ， ai =-(dV /dr+mivi z ) /mi ，
(∑ ) dz/dt= mivi2 -gKBT /Q
（11）
式中：g — 体系的自由度；
是第j个原子的核外电子在第i个原子处贡献的电荷
密度，rij是第i个原子与第j个原子之间的距离。对于
不同的金属，嵌入能函数和对势函数需要通过拟合
3 分子间的作用势
金属的宏观参数来确定。
3.1 对势
在分子动力学模拟的初期，人们经常采用的是
对势模型有以下几种：
（1）在1957年Alder和Wainwright使用间断对势[14]:
式中：σ、ε—经验常数。
（2）
（3）Born-Lande(B-J )势的形式如下[16]：
( ) φij
rij
= e2 ZiZ j + b 4πε rij rijm
（4）Morse势[17]：
（3）
根据双原子分子的振动谱[18]，提出了指数形式
的相互作用势：
( ) φ r = Ae(−αr) − Be(−βr)
4 分子动力学模拟中平衡系综的控制 方法
平衡态分子动力学模拟是在一定的系综下进行 的, 经常用的平衡系综是NVT或NPT系综。在这两种 系综中, 牵涉到控制温度和压力的几种技术,分别介
496
当代
绍如下：
4.1 控温方法
在平衡系综中，经常要把温度调整到期望值。
而体系的有效温度是由动能决定的： k=∑miv2i/ 2=fkBT/2
之间的热交换速率,Δt为分子动力学模拟的时间步
长,通过vnew=λvold校正即可保持体系的温度在T0附近 振动, 而参数Г可用于控制这个振动幅度。
（3）Nose-Hoover方法法[29-30]
Nose-Hoover热浴法是通过改变模拟体系的哈
密顿量来实现控温的,其基本思想就是在哈密顿量
里加入一个假想的项来代表一个恒温源, 具体做法
化工
2011 年 5 月
与他的控温方法类似, 如下: μ=[1+△t(P-P0)γ/Г]1/3
式中γ是一个可调参数, P0为所控制的压力。然 后, 在每一次的分子动力学模拟迭代中, 粒子坐标
x,y,z均用μ相乘, 得到新的坐标, 即可实现对压力P
系统的温度与动能存在如下的关系：
N
∑ Ek = mi vi 2 / 2 = (3N − NC ) KBT / 2 （9）
i
式中：N —原子数； NC —约束数； KB —波尔兹曼常数； vi —原子i的速度。
这种方法的基本思想是：如果t 时刻的温度是T
取了如下形式：F(ρ)=AEρlnρ,其中A，E是系数，
因物质不同而不同。
（5）Johnson势的形式如下[19]：
( ) ( ) φij rij = − An rij − Bn 3 + Cnrij − Dn
（5）
3.2 多体势
多体势于20世纪80年代初期开始出现。多体势
于20世纪80年代初期开始出现。常用的多体势函数
模拟过程是在一定体系及分子势能函数条件下从计算分子间作用力着手求解牛顿运动方程得到体系中各分子微观状态随时间的变化再将分子的位置和动量组成的微观状态对时间平均即可求出体系的压力能量粘度等宏观性质以及组成粒子的空间分布等微观结构8
第 40 卷第 5 期 2011 年 5 月
当代化工 Contemporary hemical Industry
当r≤r1时,E=∞;当r＞r1时, E=0。 （2）连续对势
Lennard-Jones(L-J )势是最老的原子间作用势
之一，是在1924年由J.E.Lennard-Jone提出来的[15]，
L-J 势的基本形式为：
φLJ
(
R)
=
4ε
⎡⎢⎢⎣⎛⎜⎝
σ R
⎞12 ⎟⎠
−
⎛ ⎜⎝
σ R
⎞6 ⎟⎠
⎤ ⎥ ⎥⎦
在高分子领域中计算机模拟的应用从 20 世纪 60 年代至今已发展到了一个崭新的阶段[1-2]。这个新阶段 的特点是，计算机模拟不仅能提供定性的描述，而且 能模拟出高分子材料结构与性能的定量的结果；模拟 分子体系算法的发展，主要是从描述简单的非真实分 子体系的算法到复杂的真实分子体系的算法。该发展 克服了用蒙特卡洛法仅能够描述不同温度下分子结 构的特征，但不能描述随温度变化而引起体系宏观物 理性质的变化中分子结构的变化过程。
比有两个优点:（1）包括显速度项;（2）收敛速度快, 系统的总势能表示为：
计算量小。这种算法明显的缺陷是位置和速度不同 步。Beeman提出的Beeman算法[12]运用了更精确的速度 表达式。它的表示式比Verlet算法复杂得多，因此计 算量较大。Gear的预测——校正算法[13]分为三步来完 成:首先, 根据Taylor展开, 预测新的位置、速度和加 速度。然后, 根据新的计算的力计算加速度。这个加 速度再由与Taylor级数展开式中的加速度进行比较。 这种方法的缺点就是占有计算机的内存大。
∑ ∑ ( ) ( ) Etot =
i
Fi
ρi
+1 2
φij
ij
ri, j
（6）
式中：F 是嵌入能：第2项是对势项，根据需
要可以取不同的形式。ρi是除第i个原子以外的所有
其它原子的核外电子在第i个原子处产生的电子云
密度之和，可以表示为：ρi = ∑ ρi (rij ) ，其中，ρi（rij） j≠i
Vol.40，No.5 May，2011
分子模拟的方法及其应用
李春艳，刘 华，刘波涛
（中北大学物理系，山西 太原 030001）
摘
要： 综述了分子模拟的方法及其在材料科学中的应用，介绍了分子模拟的基本原理、牛顿运动方程
及其有限差分算法、势函数的发展、平衡态系综及其温度和压力的控制，最后还指出了分子模拟的应用及其进
（4）
它有4个参数A,B,α和β，与L-J 势的普适形式 相类似。
3.2.2 FS势 Finnis等根据金属能带的紧束缚理论，形成了
一种在数学上等同于EAM的势函数[21]，并给出了多体
势的函数形式，把嵌入能函数设为平方根形式，即
FS势。其中的多体项及对势 ρi = ∑ ρi (rij ) 项分别为： j≠i
（t）, 速度乘以因子λ后, 温度的变化为： ΔT=（λ2-1）T（t）
T为所控制体系的温度。
（2） Berendsen热浴法法[28] 这种也称恒温槽方法，是假设所模拟的体系与
一个恒温槽连在一起, 那么两者之间就可以通过热
交换而使所模拟的体系达到恒温目的, 方法如下:
λ = 1+ ∆t (T −TC ) /τ ,式中,参数Г表征系统与恒温槽
给出位置、速度与加速度，并且不牺牲精度。这种 （EAM）[20]。EAM势的基本思想是把晶体的总势能
算法的优点是给出了显速度项，并且计算量适中， 分成两部分：一部分是晶格点阵上的原子核之间的
应用比较广泛。Hockney提出的“蛙跳”Leap–forg 相互作用对势，另一部分是原子核镶嵌在电子云中
算法[11]是Verlet算法的变化，这种算法与Verlet算法相 的嵌入能，它代表多体相互作用。在嵌入原子法中，
2.2 分子动力学运动方程的求解
（1）
在分子动力学中，为了得到原子的运动，一般采
用以下几种有限差分法来求解运动方程。常用的有：
基金项目：国家自然科学基金（批准号：50730009）。 作者简介：李春艳，女，硕士研究生，主要从事分子力学模拟研究。E-mail:lcy104104@。
第 40 卷第 5 期
（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
3.2.3 TB势
基于EAM势的势函数还有很多种[22-23]。为了把
将EAM势推广到共价键材料中去，需考虑到电子云
的非球形对称。于是Baskes等[24-26]提出了修正型嵌入
原子核法(MEAM)。总的电子密度为下式：
∑ ( ) ( )ρi
2=
3
ti(h)
ρi(h)
2
h=0
（8）
给出。其中t(h)i是常系数。在MEAM里，镶嵌能项采
李春燕，等：分子模拟的方法及其应用
495
Verlet算法[9]是在60年代后期出现的, 对扩散分 模型主要有EAM势、FS势、TB势等。
子的质心运动的积分是最稳定的也是最常用的数值 3.2.1 EAM势
方法。Swope提出的Velocity–Verlet算法[10]可以同时
Daw和Baskes在1984年首次提出了嵌入原子法
2.1 分子动力学的运动方程
分子动力学模拟的出发点是假定粒子的运动
可以用经典动力学来处理, 对一个由N个粒子构成
的孤立体系, 粒子的运动由牛顿运动方程决定, 也
就是
m
d2ri dt 2
=-∇iV (r1,r2,rN
)
式中：mi，ri —第i个原子的质量和位置；
▽i = -d/dri,V(r1,r2,rN)—体系所处的势。
Q — 一个可调参量, 表征着假想项的质量, T为温度。
4.2 压力控制
分子动力学模拟通常要求模拟体系的压力保
持不变，而体积可以随温度的变化而改变，即进行
等压模拟。这可以通过改变模拟原胞的3个方向或一
个方向的尺寸来实现体积的变化。通常有： 4.2.1 Berendsen方法[31]
这种控压方法的基本思路与基本控制方程均
1.2 分子动力学方法(简称MD)的基本原理 分子动力学方法(MD)是另一种主要的计算机模
拟方法，在材料科学、物理、化学等学科的各个领 域得到广泛应用。模拟过程是在一定体系及分子势 能函数条件下，从计算分子间作用力着手，求解牛 顿运动方程，得到体系中各分子微观状态随时间的 变化，再将分子的位置和动量组成的微观状态对时 间平均，即可求出体系的压力、能量、粘度等宏观 性质以及组成粒子的空间分布等微观结构[8]。
1 分子模拟方法
2 分子模拟的计算
1.1 蒙特卡洛方法(简称MC)的基本原理 MC法早在20世纪50年代就被应用到高分子科
学的研究之中。先驱性的工作是Wall在为研究高分 子链的排除体积问题时所做的蒙特卡洛模拟[3-5]。基 本模拟过程是在一定系统条件下，将系统内粒子进 行随机的位移、转动，或粒子在两相间转移位置。 根据给定的分子势能函数，进行粒子间内能的加和。 采用Metropolis取样方法[6-7]，生成一系列体系的微观 粒子随机构型，从而逐渐趋近于平衡时的Boltzmann 分布。目前蒙特卡洛的应用几乎涉及到高分子科 学中的各个领域，为现代高分子科学理论基础的 建立和发展起到了十分重要的推动作用。
( ) ∑ ( ) ( ) ρ
rij
=
2
Ak
Rk − rij
3
H
RK − rij
k =1
( ) ∑ ( ) ( ) σij
rij
=
6
ak
rk − rij
3
H
rk − rij
k =1
（7）
式中：当x>0时，H (x)=0；当x<0时，H(x)=1。
Ak、RK、aK、rK为常数，且有R1＞R2,r1＞r2＞…＞r6。 它们的值随金属物质的不同而有所不同。
一步的研究方向。
关 键 词： 分子模拟; 有限差分法; 原子间相互作用势; 材料科学
中图分类号： O 642
文献标识码： A 文章编号： 1671-0460（2011）05-0494-03
Method of Molecular Simulation and Its Application
LI Chun-yan, LIU Hua, LIU Bo-tao
如下： H=∑miv2i/2+V(r)+Qζ2/2+gKT lnS
式中：S和ζ分别是假想项的坐标和动量。
（10）
这样体系的微分方程就变为：
vi = dri / dr ， ai =-(dV /dr+mivi z ) /mi ，
(∑ ) dz/dt= mivi2 -gKBT /Q
（11）
式中：g — 体系的自由度；
是第j个原子的核外电子在第i个原子处贡献的电荷
密度，rij是第i个原子与第j个原子之间的距离。对于
不同的金属，嵌入能函数和对势函数需要通过拟合
3 分子间的作用势
金属的宏观参数来确定。
3.1 对势
在分子动力学模拟的初期，人们经常采用的是
对势模型有以下几种：
（1）在1957年Alder和Wainwright使用间断对势[14]:
式中：σ、ε—经验常数。
（2）
（3）Born-Lande(B-J )势的形式如下[16]：
( ) φij
rij
= e2 ZiZ j + b 4πε rij rijm
（4）Morse势[17]：
（3）
根据双原子分子的振动谱[18]，提出了指数形式
的相互作用势：
( ) φ r = Ae(−αr) − Be(−βr)
4 分子动力学模拟中平衡系综的控制 方法
平衡态分子动力学模拟是在一定的系综下进行 的, 经常用的平衡系综是NVT或NPT系综。在这两种 系综中, 牵涉到控制温度和压力的几种技术,分别介
496
当代
绍如下：
4.1 控温方法
在平衡系综中，经常要把温度调整到期望值。
而体系的有效温度是由动能决定的： k=∑miv2i/ 2=fkBT/2
之间的热交换速率,Δt为分子动力学模拟的时间步
长,通过vnew=λvold校正即可保持体系的温度在T0附近 振动, 而参数Г可用于控制这个振动幅度。
（3）Nose-Hoover方法法[29-30]
Nose-Hoover热浴法是通过改变模拟体系的哈
密顿量来实现控温的,其基本思想就是在哈密顿量
里加入一个假想的项来代表一个恒温源, 具体做法
化工
2011 年 5 月
与他的控温方法类似, 如下: μ=[1+△t(P-P0)γ/Г]1/3
式中γ是一个可调参数, P0为所控制的压力。然 后, 在每一次的分子动力学模拟迭代中, 粒子坐标
x,y,z均用μ相乘, 得到新的坐标, 即可实现对压力P