《装备制造技术》 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]。pppd-298
$(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微正则系综
(NVE 微正则系综是孤立的、保守的系统 , 即所模拟的系统中粒子数 N 、体积 V 、能量 E 都不变。此外 , 由于整个系统并未运动 , 所以整个系统的总动量为零。在系统的演化过程中 , 系统沿着相空间中的恒定能量轨道演化。在模拟中不需要对系统进行能量控制 , 系统能量自然保持守恒。 4.2正则系综
(NVT 在正则系综下 , 系统的粒子数 N 、体积 V 和温度 T 都保持不变。在恒温下 , 系统的总能量不是一个守恒量 , 系统要与外界发生能量交换。保持系统的温度不变 , 通常运用的方法是让系统与外界的热浴处于热平衡状态 , 也可以通过对速度进行直接标度来实现恒定的温度。 4.3等温等压系综
(NPT 等温等压系综具有确定的粒子数 N 、恒定的温度 T 和恒定的压力 P 。温度控制和以前一样 , 通过对速度加一些约束来实现。而对压力的调节 , 比较复杂。由于系统的压力与其体积是共轭量 , 要调节压力值可以通过标度系统的体积来实现。 4.4等焓等压系综
(NHP 等焓等压系综是保持系统的粒子数 N 、压力 P 和焓值 H 都不变。模拟时要保持压力与焓值为固定值 , 有一定难度。事实上 , 这种系综在实际的分子动力学模拟中很少见。
5宏观统计量的提取
5.1温度
阻尼系数温度 T 直接与粒子的动能相关 , 即著名的均匀分布公式。
粒子的每个自由度赋予 K B T /2的能量 , N 个粒子的总自由度为
3N , 故系统温度 T 与动能 K 的关系为 [11]:
K=N
i =1%
m i v 2
=2NK T (12
式中 , K B 为 Boltzmann 常数 , v i 为原子 i 速度。
5.2压力
PV =NK B T +1i
%r i ? F i 〉 (13
这是著名的 Virial 方程 [11]。式中 , P 为压力 , V 为元胞体积 ,
T 为系统温度 , N 为系统总的粒子数目 , D 为系统维数 , K B 为Boltzmann 常数 , F i 为原子 i 所受的原子间作用力 , r i 为原子 i
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议