第三章_经典分子动力学方法资料

《第三章_经典分子动力学方法资料》由会员分享，可在线阅读，更多相关《第三章_经典分子动力学方法资料（39页珍藏版）》请在金锄头文库上搜索。

1、第三章 经典分子动力学方法 3.1 引言 分子动力学 (Molecular Dynamics, 简写为MD) 方法是确 定性模拟方法，这方法是按该体系内部的内禀动力学规 律来计算确定位形的转变。 首先需要建立一组分子的运动方程，然后通过直接对系 统中的每一个原子/分子运动方程进行数值求解，得到每 个时刻每个原子/分子的坐标与动量(速度)，即在相空间 的运动轨迹，再利用统计方法得到多体系统的静态和动 态特性，从而得到系统的宏观性质。 在MD方法的处理过程中，方程组的建立是通过对物理体 系的微观数学描述给出的，在这个微观的物理体系中， 每个原子/分子都各自服从经典的牛顿力学定律。 MD方法是实现玻

2、尔兹曼的统计力学途径，可以处理与时 间有关的过程，因而可以处理非平衡态问题，但是该方 法的计算机程序较复杂，计算量大，占内存也多。 MD方法的发展史 MD方法是20世纪50年代后期由B.J Alder和T.E. Wainwright创造发展的。他们在1957年利用MD 方法，发现了早在1939年根据统计力学预言的“ 刚性球组成的集合系统会发生由其液相到结晶相 的相转变”。 20世纪70年代，产生了刚性体系的动力学方法 被应用于水和氮等分子性溶液体系的处理，取得 了成功。1972年，A.W. Less和S.F. Edwards等 人发展了该方法，并扩展到了存在速度梯度(即 处于非平衡状态)的系统

3、。 之后，此方法被M.J. Gillan等人推广到了具有温 度梯度的非平衡系统，从而构造并形成了非平衡 MD方法体系。 MD方法的发展史 MD方法真正作为材料科学领域的一个重要研究 方法，开始于恒压MD方法(1980)和恒温MD方法 (1984)的建立及在应用方面的成功。 1985年人们又提出了将电子论和分子动力学方 法有机统一起来的所谓CarParrinello方法，即 第一性原理MD方法。它不仅可以处理半导体和 金属的问题，同时还可应用于处理有机物和化学 反应。 1991年有人进一步提出了巨正则系综MD方法， 从而又可适用于吸附问题的处理等，该方法还在 进一步发展之中。 分子动力学方法的主

4、要发展可见表3.1。 年代创立者创造内容工作（MD分类名称） 1957 B. J Alder & T. E. Wainwright 刚性球MD方法 1963A. Rahman质点系MD方法 1971 A. Rahman & F. H. Stillinger 刚性系统MD方法 1972 A. W. Lees & S. F. Edwards 平衡系统MD方法（存在速度梯度） 1977J. P. Rychaert et al.约束系统MD方法 1980 Andersen, Parrinello & Rahman 恒压压MD方法 1983 N. J. Gillan & M. Dixon 非平衡MD方法

5、（存在温度梯度） 1984 S. Nos恒温MD方法 1985 R. Car & M. Parrinello 第一性原理MD方法 (Car-Parrinello方法) 1991Cagin & Pettitt巨正则系统MD方法 表3.1 MD方法的里程碑工作 3.2 MD方法计算初步 在计算机出现以前，作为根据原子间相互作用力 等微观信息了解多原子或分子团的结构、性质的 方法，所采用的是基于统计理论的数学解析法。 然而，原子间相互作用力稍微复杂一些，不用说 求解统计理论严格方程解，就是进行数值求解也 是一件很困难的事。 MD方法就是数值求解多体系统的确定性运动方 程，并根据对所求结果进行统计处理

6、，决定粒子 的轨迹，从而给出物性预测和微观结构信息的一 种模拟方法。 内能 比热容 运动方程 温度、压力 相互作用 原子位置坐标 3维结构 原子坐标、速度 原子运动 热力学性质 动力学性质 光学性质 （输出信息）（二次信息） 扩散系数 粘滞系数 电导率 红外吸收 图3-1 MD方法信息输入输出信息方框图 MD这种方法并不严格。因此，必须根据 情况，检验改变所模拟的基本单元尺寸所 得结果是否会改变，直到所得结果不随基 本单元尺寸变化而变化。通常这样的处理 在很多情况下是有效的。 对于基本单元内的原子、分子运动方程，使用什么样的 形式合适，要具体问题具体分析。若是考虑具有确定的 粒子数N，体积V和

7、能量E的NEV系综(称为微正则系综， Micro-Canonical Ensemble)，则其运动方程可以表达成 式(3-2-1)所示的普通牛顿方程的形式 (3-2-1) 式中mi为所考察的原子质量，ri为原子的位置坐标，Fi为 作用在原子上的原子相互作用的合力，它由下式给出 (3-2-2) 其中，ij是原子和原子j之间的势函数(有时亦称为力场) 例如，由氩原子等组成的稀薄气体，其势 函数可采用LennardJones势， (3-2-3) 式中，r 是原子间距， 是结合强度参数 ， 是表示原子半径的参数。 在t时间内，对系统内的所有 粒子解运动方程 Yes No 启动计算 设定坐标、速度初始值

8、 计算作用在原子上的力Fi 计算要求的物理量，将数据写 入轨迹文件 ttma x 输出计算结果，并结束计算 对(3-2-1)可用数值积 分法求解，其数据处 理流程图见图3-2 图3-2 MD数据处理流程图 MD方法 NEV 能量恒定 NTV 恒温 NHP 恒压 NTP 恒温恒压 VT 巨正则系 VL 恒化学势 弹性力学 （原子分子 ） 质点力学 （原子分子 ） 刚体力学 （分子） 约束力（分 子和晶体） 系） 动力学模型 目标系统 团簇 块体材料 表面 界面 MBE/CVD Berrele 法 Green 法 多重时间 刻度法 数值积分法 边界条件 到目前为止已经确立的MD方法的主要技术体系

9、统计系综 系综是一个巨大的系统，由组成、性质、尺寸和形状完 全一样的全同体系构成数目极多的系统的集合。不同的 系综，MD方法的基本方程有所不同。 目前除微正则系综（NEV系综）外，已完成了正则系综 （NTV系综），等温等压系综（NTP系综），等压等焓 系综（NHP系综），巨正则系综（VL系综），恒定化 学势系综（VT系综）等五个系综的MD方法的基本方程 的确立。 现在已经能够处理许多体系，例如： 孤立宏观团簇的模拟(用NEV或NTV系综) 固体的结构相变，玻璃转变，晶化过程的模拟(用NTP系综) 固体(晶体)表面的原子、分子吸附现象的模拟 (用VT系综) 力学条件 已建立了弹性力学、质点力学、

10、刚体力学、约 束力学等不同力学条件下的四种体系的MD方法 。 弹性力学方法是将所考察的原子分子看作刚性球 来处理，建立完全弹性碰撞方程，借以求解出原 子、分子的运动规律。这种处理可以在液晶的模 拟中使用。 质点力学模型是将原子、分子作为质点处理，粒 子间的相互作用力采用坐标的连续函数。这种力 学体系的应用对象非常多，可以用于处理陶瓷、 金属、半导体等无机化合物材料以及有机高分子 、生物大分子等几乎所有的材料。 力学条件 刚体力学方法是把分子作为刚体处理，建立对于 刚体的欧拉方程和对于质点系的牛顿方程，联立 求解所研究对象的运动问题，以前主要是在处理 像水和四氯化碳那样的低分子量体系，现已用于

11、研究晶体的相变。 约束力学是冻结粒子体系的一部分自由度，进而 求解因此而生成约束条件下的质点力学运动方程 。对于有机分子来讲，因为键长、键角的振动变 化对体系的结构影响不大，将这些自由度冻结是 合适的。另外也有采用固定晶体结构的考虑。 边界条件问题 在处理原子、分子的聚合体问题时，就MD方法 而言，能处理的原子(分子)数目要受到计算机运 行速度和能力的限制。目前报导的最好水平是能 处理109量级的原子数目。这与现实物质含有 1023个原子或分子的差距还很大，导致模拟系统 原子数少于真实系统的所谓“尺寸效应”的问题。 为了减小“尺寸效应”而又不至于使计算工作量过 大，对于平衡态MD模拟采用 “周

12、期性边界条件” 。 周期性边界条件 图3-4 假设现实物质中的 一部分原子(通常为 102105个)，被取 出配置在所谓基本 单元的箱中，由于 基本单元周围的原 子、分子变成表面 ，从而不同于本来 要处理的体状态。 为了防止这种情况 ，在基本单元周围 配置其复制品(图3- 4)。 周期性边界条件 对于在基本单元周围边界的原子、分子所受到的 作用力，不仅有来自基本单元之内部的原子或分 子的作用力的贡献，还要考虑来自其近邻复制单 元的原子或分子的作用力的贡献。这样，把来自 距粒子某一距离（截断距离）内的j粒子的贡献 （这里与粒子j是否处在基本单元内或复制单元 内无关）求和。给出力的方法称为最小镜像

13、距离 法（Minimum Image Distance）。 周期性边界条件 如果在基本单元内的原子的位置为, 周期边界条 件会产生该原子的镜像。它的位置在 这儿a、b、c是基本单元的三个边长，l，m，n 是整数，取值范围为- 到+。 在基本单元内的原子不仅与在本单元内的其他 原子有相互作用，还与在相邻单元内的镜像原子 有作用。 (3-2-4) 基本单元大小的选择 基本单元的大小必须大于2Rcut（Rcut是相互作用 势的截断距离）或Rcut1012 m2/s时适用。 (3-3- 10) (3-3-9) GreenKubo formula 扩散系数D也可用速度自相关函数通过下 面的GreenKu

14、bo 公式表示 （3311 ） GreenKubo公式和Einstein 关系是等价 的。 Partial diffusion 在多种原子组成的体系中（例如AB2）， 可以对某一种原子求和，得到Partial diffusion，从而得到原子A或原子B的扩散 系数DA、DB。 非均匀体系中的扩散 对一非均匀体系，例如一块材料，有体和 晶界部分。可分别对某一类原子做平均， 如对晶界原子，或对体中的原子，得到不 同部分原子不同的扩散。 与温度有关的扩散 在一个系统里，若原子的扩散是通过跳跃进行的 ，可将扩散用Arrhenius关系 来表达，式中是prefactor，是原子迁移率或激 发能(acti

15、vation energy)。对（3312）两边 求对数有 (3-3- 12) (3-3- 13) 与温度有关的扩散 作lnD1/T图(图3.8)，从图中直线的斜率可以得到Ea， 从截距可以得到D0。同理也可以得到表面原子、晶界原 子的扩散系数(图3.9)。 图3.8 图3.9 MD模拟材料的溶解 固态到液态的相变（溶解）的过程发生在几ps时 间，在MD模拟中可通过以下手段知道是否发生 了固-液相变： 观察不同时间段系统的原子位置情况（snapshots） 能量从动能向势能转移 总能量和势能随温度的变化有跳跃 扩散系数有变化 径向分布函数有变化 在做有固定体积的模拟时，压力有突变 MD模拟材料的溶解 溶点对应于固-液共存状态，这时在 computational cell 中应含有固态、液态及 固-液界面。由图310可以得到溶点温度 。 图310 体膨胀系数可以由VT曲线（图311）得到 图311 对不同的压强P做模拟，可以得到 （T，P）。