分子动力学模拟方法概述(精)

《装备制造技术》 2007年第 10期

收稿日期 :2007-08-21

作者简介 :申海兰 , 24岁 , 女 , 河北人 , 在读研究生 , 研究方向为微机电系统。

分子动力学模拟方法概述

申海兰 , 赵靖松

(西安电子科技大学机电工程学院 , 陕西西安 710071

摘要 :介绍了分子动力学模拟的基本原理及常用的原子间相互作用势 , 如Lennard-Jones 势 ; 论述了几种常用的有限差分算法 , 如 Verlet 算法 ; 说明了分子动力学模拟的几种系综及感兴趣的宏观统计量的提取。关键词 :分子动力学模拟 ; 原子间相互作用势 ; 有限差分算法 ; 系综中图分类号 :O3

文献标识码 :A

文章编号 :1672-545X(200710-0029-02

从统计物理学中衍生出来的分子动力学模拟方法 (molec-

ular dynamics simulation , M DS , 实践证明是一种描述纳米科技

研究对象的有效方法 , 得到越来越广泛的重视。所谓分子动力学模拟 , 是指对于原子核和电子所构成的多体系统 , 用计算机模拟原子核的运动过程 , 从而计算系统的结构和性质 , 其中每一个原子核被视为在全部其他原子核和电子所提供的经验势场作用下按牛顿定律运动 [1]。它被认为是本世纪以来除理论分析和实验观察之外的第三种科学研究手段 , 称之为“计算机实验” 手段 [2], 在物理学、化学、生物学和材料科学等许多领域中得到广泛地应用。

根据模拟对象的不同 , 将它分为平衡态分子动力学模拟 (EM DS (和非平衡态分子动力学模拟 (NEM DS 。其中 , EM DS 是分子动力学模拟的基础 ; NEM DS 适用于非线性响应系统的模拟 [3]。下面主要介绍 EM DS 。

1分子动力学方法的基本原理

计算中根据以下基本假设 [4]:

(1 所有粒子的运动都遵循经典牛顿力学规律。 (2 粒子之间的相互作用满足叠加原理。

显然这两条忽略了量子效应和多体作用 , 与真实物理系统存在一定差别 , 仍然属于近似计算。

假设 N 为模拟系统的原子数 , 第 i 个原子的质量为 m i , 位置坐标向量为 r i , 速度为 v i =r ・ i , 加速度为 a i =r ・・

i , 受到的作用力为 F i , 原子 i 与原子 j 之间距离为 r ij =r i -r j , 原子 j 对原子 i 的作用力为 f ij , 原子 i 和原子 j 相互作用势能为 ! (r ij , 系统总的势能为 U (r 1, r 2, K r N =

N

i =1! j ≠ i

!

" (r ij , 所有的物理量都是随时

间变化的 , 即 A=A (t , 控制方程如下 :

m i r ・・

i =F i =j ≠ i

! f

ij

(1 F i =-#r i

U (r 1, r 2, K r N

(2

以此建立一个线性的微分方程组 , 给定初始位置和速度 , 方程是封闭的 , 可以得到任意时刻系统中所有原子的位置

r i (t 和速度 v i (t 。

2势函数

2.1Lennard-Jones

(L-J 势 M ie 最先提出了两体势的解析形式 , Lennard-Jones [5]将它

应用到铜的自扩散研究计算中。该两体势可表示为

#(r =4ε(

σr 12-(σr

6

%

(3

式中ε和σ是势参数 , r 表示原子间距。通常为更好地拟合元素已有的实验数据 , 如结合能、晶格常数等 , L-J 势则使用一种普适的形式 [6]。

$(r =4ε(n (

σr m -(n (σr n

%

(4

2.2Morse 势

根据双原子分子的振动谱 [7], 提出了指数形式的相互作用

势

%(r =Aexp (-αr -B exp (-βr (5

它有 4个势参数A , B , α和β, 与 L-J 势的普适形式相类

似。

2.3Embedded Atom Method

(EAM 势 Baskes 和 Daw [8]基于密度函数理论和准原子近似理论 , 导出了嵌入原子理论模型势 , 能量表示为

E to t =i

! F i (ρi +1

i , j

! i ≠ j

&ij (

r ij (6

式中第一项为嵌入能项 , 表示原子嵌入到电子密度为ρi

处的能量 , 第二项为两体相互作用项 , 而基体电子密度则表示为原子电子密度的线性叠加 , 即

ρi =j (≠ i

! f ij (r ij

(7

3有限差分算法

3.1Verlet 算法 [9]

将t +Δt 和 t -Δt 时粒子的位置坐标分别用时刻 t 的位置坐标作泰勒展开有 : [下转第 34页 ]

r (t +Δt =r (t +Δt・ V

(t +(Δt 2a (t +L r (t -Δt =r (t -Δt ・ V

(t +(Δt 2a (t +L ! ####" ####$

(8

由式 (8 可得t +Δt 时刻粒子的位置为:r (t +Δt B2r (t -r (t -Δt +

(Δt 2F (t (9

相应的速度为 :

V

(t B r (t +Δt -r (t -Δt (10 式中 , m , V (t , a (t 和 F (t 分别为原子的质量、速度、加速

度以及所受到的力。

3.2Leap-frog 算法 [10]

Hockney 对 Verlet 算法进行了改进 , 提出 Leap-frog 算法:r (t +Δt =r (t +Δt ・ V (t +Δt V (t +Δt =V (t -Δt +Δt F (t ! #

#

##" #

#

##

$

(11 3.3其它算法

Verlet 和 Leap-frog 算法只能求解线性常微分方程。为了

求解非线性常微分方程 , Gear 提出了基于预测 -校正积分方法的 Gear 算法。 Swope 提出的 Velocity-Verlet 算法可同时得出位置、速度与加速度 , 且不牺牲精度 , 优点是给出了显式速度项 , 计算量适中。 Beeman 提出的 Beeman 算法运用了更精确的速度表达式 , 能更精确地计算系统动能 , 但表达式很复杂 , 计算量很大。

4不同系综的分子动力学模拟

4.1微正则系综