《装备制造技术》 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 =Ni =1! j ≠ i!" (r ij , 所有的物理量都是随时间变化的 , 即 A=A (t , 控制方程如下 :m i r ・・i =F i =j ≠ i! fij(1 F i =-#r iU (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-(σr6%(3式中ε和σ是势参数 , r 表示原子间距。
通常为更好地拟合元素已有的实验数据 , 如结合能、晶格常数等 , L-J 势则使用一种普适的形式 [6]。
$(r =4ε(n (σr m -(n (σr n%(42.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 +1i , j! i ≠ j&ij (r ij (6式中第一项为嵌入能项 , 表示原子嵌入到电子密度为ρi处的能量 , 第二项为两体相互作用项 , 而基体电子密度则表示为原子电子密度的线性叠加 , 即ρi =j (≠ i! f ij (r ij(73有限差分算法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=Ni =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的位置 , 〈〉表示对括号内的量在 t 时刻求平均。
5.3系统总能量E =E k +E p =〈 1i %m i ・ v i 2〉+〈 i%u ij ・ (r ij 〉 (14式中右端第一项为系统总动能 , 第二项为系统总势能。
6结束语分子动力学方法是用运动方程来计算系统的性质 , 结果得到的既有系统的静态特性 , 也有动态特性 [12]。
M D 模拟方法面临着两个基本限制 :一个是有限观测时间的限制 ; 另一个是有限系统大小的限制。
自 20世纪 50年代中期开始 , M D 方法得到了广泛的应用。
目前 , 分子动力学方法还在进一步发展中。
参考文献 :[1][2][3][4][5][6][7][8][9][10][11][12]文玉华 , 朱如曾 , 周富信 , 等 . 分子动力学模拟的主要技术 [J]. 力学进展 , 2003, 33(1:65-73.Heermann D.W. Computer Simulation M ethods in Theoretical Physics[M ]. Springer-Verlag Press,1990.Liem S Y , Chen K Y. Simulation study of platinum adsorption on.graphite using the Sutton – Chen potential [J].Surf. Science, 1995, (328 :119-128.温诗铸 . 纳米摩擦学 [M ]. 北京 :清华大学出版社 , 1998:60-80.Lennard-Jones J E. On the determination of molecular fields [J]. Pro-ceedings of the Poyal Society. 1924,A(106 :463-469. Zhen Shu, Davies G J. Calculation of the Lennard-Jones nm potentialenergy parameters for metals[J]. Phys Stat Sol, 1983,A(78 :595-605. M orse P M . Diatomic molecules according to the ware mechanics ⅡVi-brational Levels[J]. Phys Rev, 1929,(34 :57-64. Daw M S, Baskes M I. Embedded atom method:derivation and applica-tion to impurties, surfaces, and other defects in metals [J]. Phys Rev, 1984,B(29 :6443-6453. Verlet L. Computer ‘ experiments ’ on classical fluidsⅠ :Thermody-namical properties of Lennard-Jones molecules [J]. Pysical Review, 1967,(159 :98-103. Honeycutt R W. The potential calculation and some application [J].M ethods in Computational Pysics, 1970,(9 :136-211. 吴恒安 . 纳米尺度下结构和材料力学行为的分子动力学模拟研究[D]. 合肥 :中国科学技术大学 , 2002.张会臣 , 严立 . 纳米尺度润滑理论及应用 [M ]. 北京 :化学工业出版社 ,2005:182-185.[上接第 30页 ]The Overview of Molecular Dynamics SimulationSHEN Hai-lan , ZHAO Jing-song(M echanical Engineering School of Xidian University, Xi ’ an 710071, ChinaAbstract :This article discusses the principle of molecular dynamics, interatomic potentials (for example, Lennard-Jones potential, somecommon finite difference algorithms (such as Verlet algorithm, several molecular dynamics simulation ensembles, and the extraction of useful information.Key words :molecular dynamics ; interatomic potential ; finite difference algorithm ; ensemblesDesign and Development of Dynamic Balance Experiment Device for the Dual-RotorSynchronized SystemLI Xi-ling , ZHANG Zhi-xin , HE Shi-zheng(Institute of Chemistry and M achinery Engineering, Zhejiang University, Hangzhou 310027, ChinaAbstract :This paper designs an experimental machine of dual-rotor for dynamic balance. And the system for getting signal based on Lab-VIEW is proposed. Through the machine designed in this paper the basis and vibration signal can be obtained, and several experiment about dynamic-balance can be taken.Key words :Dual-rotor system ; Dynamic-balance图 , 转鼓材质为无缝钢管 , 钢管端面焊接法兰用于连接端盖。