分子动力学应用

1、6. 分子动力学原理与方法 2013年9月4日,计算材料学,关于研究报告,选取一种计算方法和一个研究对象，进行计算获得数据结果，并得出结论，完成专题研究。 课程报告的选题，应该具有一定的意义，能够说明问题。另外，选题应该尽量避免雷同。 根据研究结果编写研究报告，报告一般包括研究目的，研究方法，计算模型，结果和数据，结论。 计算中的原子模型更选用白色背底，通过缩放和旋转，以最佳的角度反应出模型中的结构。 数据画度要调整图形的标度，让最的用的数清楚地显示出来。必要时可以导出数据，采用专业的画图软件画图。 报告应该尽量的条理清楚，逻辑性强，结论可靠。 课程研究应独立完成，或者以2人为小组进行。,本课

2、提纲,本节课将向同学位介绍分子动力学的应用。 系综简介 统计平均值 静态性能分析 动态性能分析 GULP的使用 分子动力学计算范例,1.系综简介,系综（Ensemble)是统计力学的一个概念，它是1901年由吉布斯创立完成的。分子动力学所研究的对象是多粒子体系，统计物理的规律仍然成立，因此计算机模拟的多粒子体系用统计物理的规律来描述。与微观量相对 应的宏观量是在一定的宏观条件下所有可能的运动状态的平均值。 分子动力学模拟方法中包括平衡态和非平衡态模拟。根据研究对象的特性，主要的系综有微正则系综（NVE)、正则系综（NVT)、等温等压系综（NPT)、等焓等压系综 (NPH)等。 采用分子动力学模

3、拟时，必须要在一定的系综下进行,1.系综简介,微正则系综 微正则系综，又称NVE系综，它是孤立的、保守的系统的统计系综。 在这种系综中，体系与外界不交换能量，体系的粒子数守恒，体系的体积也不发生变化，系统沿着相空间中的恒定能量轨道演化。 在分子动力学模拟中，通常用时间平均代替系综平均。 在微正则系综中，轨道(坐标和动量轨迹)在一切具有同一能量的相同体积内经历相同的时间，则轨道平均等于微正则系综平均。,1.系综简介,正则系综 在热力学统计物理中正则系综是一个粒子数为N、体积为V、温度为T和总动量为守恒 量的系综，在这个系综中系统的粒子数（N)、体积(V)和温度（T)都保持不变，并且总动量为 零，

4、因此为称NVT系综。 在恒温下，系统的总能量不是一个守恒量，系统要与外界发生能量交换。保持系统的温度不变，通常运用的方法是让系统与外界的热浴处于热平衡状态。由于温度与系统的动能有直接的关系，通常的做法是把系统的动能固定在一个给定值上，这是对速 度进行标度来实现的。,1.系综简介,等温等压系综 等温等压系综，即NPT系综，就是系统处于等温、等压的外部环境下的系综。 在这种系综下，体系的粒子数(N)、压力(P)和温度(T)都保持不变。 这种系综是最常见的系综，许多分子动力学模拟都要在这个系综下进行。 这时，不仅要保证系统的温度恒定，还要保持它的压力恒 定，温度的恒定和以前一样，是通过调节系统的速度

5、来实现的。 由于系统的压力P与其体积V是共辄量，要调节压力值可以通过标度系统的体积来实现。目前有许多调压的方法都是采用这个原理。,1.系综简介,等温等压系综,在分子动力学计算中实现等温等压的方法,2. 统计平均值,微观领域往往研究单个粒子的行为，宏观性质是大量粒子的综合行为。分子动力学 (MD)方法能够再现宏观行为，同时又存储了大量的微观信息，因此是联系宏观和微观的重要 工具。 利用此方法可以研究由热力学统计物理能够给出的各种性能参数。统计力学将系统的微观量与宏观量通过统计物理联系起来。 在运行分子动力学程序之后，得到了系统的所有粒子的坐标和速度随时间的变化轨迹。 接下来研究我们所感兴趣体系的

6、性质。最简单的是热力学性质，如温度、压力、热容。,2. 统计平均值,这些量可由体系的坐标和动量的统计平均得到，称为静态性能。但有一类热力学性质不能在一次模拟中直接得到。也就是说，这些性质不能表达为体系中所有粒子坐标和动量的一些函数的简单平均，称为动态性能。 物性参量可以根据原子的坐标和速度通过统计处理得出，在统计物理中可以利用系综微观量的统计平均值来计算物性参量值 在分子动力学中，使用了时间平均等于系综平均的各态历经假设。虽然各态历经假设在热力学统计物理中没有证明，但它的正确性已被实验结果证明是正确的。,3.静态性能分析,体系的热力学性质和结构性质都不依赖于体系的时间演化，它们是静态平衡性质。

7、 温度T 在正则系综（NVT)中，体系的温度为一常数；然而在微正则系综中，温度将发生涨落。 温度是体系最基本的热力学量，它直接与系统的动能有关， N为粒子总数，Nc为受限制的自由度数目，通常N=3,3.静态性能分析,能量E 体系的动能EK、势能U、总能量E可以由下式给出 压强P 压力通常通过虚功原理模拟得到。虚功定义为所有粒子坐标与作用在粒子上的力的乘积的和。,3.静态性能分析,径向分布函数 径向分布函数(radial distribution function)是描述系统结构的很有用的方法，是距离一个原子为r时找到另一个原子的概率，表示为g(r)。 径向分布函数在模拟过程中很容易计算，因为所

8、有的距离在计算力时巳完成。它表征着结构的无序化程度。 在晶体中，径向分布函数g(r)有无限个尖锐的峰，位置和高度由晶体结构决定。,4.动太性能分析,关联函数 分子动力学可产生与时间有关的体系的构型，因此可以用来计算与时间有关的性质。 与时间有关的性质通常通过时间关联函数来计算。体系的动态性质必须通过体系的动力学轨迹得到。 假设有两套数据x和y，要确定它们之间在一定条件下的关联。关联函数(关联系数)是提供数据之间关联强度的一个量。可以定义很多关联函数：,4.动太性能分析,时间关联函数 一分子动力学模拟可以提供特定时刻的值，这样使得我们可以计算一个时刻的物理量与同 一时刻或另一时刻（时间t以后）的

9、另一物理量的关联函数，这个值被称为时间关联系数： 如果 和是不同的物理量，则关联函数称为交叉关联函数。 如果 和y是同一量，则关联函数称为自关联函数。 自关联函数就是一个量对先前的值的记忆程度，或者反过来说，就是系统需要多长时间忘记先前 的值。 一 些关联函数可以通过系统内所有粒子求平均得到，而另外一些关联函数是整个系统粒子的函数。,4.动太性能分析,弛豫时间 一个自关联函数例如速度自关联函数初始值为1，随时间的增加，变为0。关联函数从1变为0的时间称为关联时间，或弛豫时间。 两种密度的气体初始速 度关联函数为1，然后随时间衰减到0。对低密 度，速度关联函数逐步衰减到0;而对高密度的情 况，（

10、0越过轴变为负值，然后又变为0。负的 速度关联函数意义就是粒子以与0时刻速度相 反的方向运动。,4.动太性能分析,输运性质 输运性质是指物质从一个区域流动到另一个区域的现象，比如非平衡溶质分布的溶液，溶质原子会发生扩散直到溶质浓度均匀。 如果体系存在温度梯度，就会发生能量输运直到温度达到平衡，动量梯度产生粘滞性。 输运意味着体系处于非平衡态,处理非平衡态的分子动力学方法这里不予讨论, 但可以用平衡态模拟中的微观局域涨落来实现非平衡态性质的计算。 当然也应该意识到用非平衡态分子动力学来计算非平衡性质更有效。,4.动太性能分析,扩散 扩散的通量用Fick第一定律来描述，即 Jz =-D(dN/dz

11、) 。扩散行为随时间的演化由Fick第二定律来描述： Fick第二定律的解为： 利用爱因斯坦关系可以在平衡模拟中计算扩散系数、平均平方位移与时间的曲线：,5. GULP的介绍,Material Studio中有几个专门用于分子动力学的模块，包括Discovery, Forcite, GULP等。 GULP是和款分子动力学的计算程序，也被收入到Material Studio中。 除了在Material Studio中使用之外，还可以从GULP的网站上下载，编译后使用。网址：/gulp/。 GULP对于学术界免费，注册确认后可以下截。但商业用户除外。

12、对于免费用户，GULP不提供支持。直接拿来用就可以了。 对于本课程，我们可以直接在Material Studio中使用。,5. GULP的介绍,GULP的计算功能非常强大，主要表现如下： 能够处理的系统类型包括： 0D：团簇，点缺陷 1D：高分子，线缺陷， 2D：表面，薄层，晶界 3D：块体材料 能量最小化方法，包括： 等体积或等压强 约束条件下的结构优化5. GULP的介绍 Newton/Raphson, conjugate gradients 或 有理函数Optimisers DFP or BFGS 法更新hessian矩阵,5. GULP的介绍,晶体性质计算，包括： 力学性质，如弹性常数

13、，体弹性模量，杨氏模量，泊松比，剪切模量 光学性质上，如静态介电常数，高频介电常数，折射率，反射率，频率有关的介电常数张量 压电常数 声子谱，包括声子频率，DOS，PDOS，色散曲线 热力学性质，如熵，比热，自由能 电学性质，如静电势，电场，电场梯度，波恩有效电荷 分子动力学，包括 NVE,NVT和NPT系综 壳层模式分子动力学 绝热算法的壳层外推,5. GULP的介绍,计算类型的选取： 单点能量 结构优化 分子动力学 力场拟合 表面计算 力场的选取： 有大量的力场可供选用 电荷的设定 电场的设定,GULP的计算参数设置,5. GULP的介绍,分子动力学的计算设置： 选取系综 选取温度 选取压

14、强 选取步长 选取步数 选取输出的间隔 进行分子动力学计算可以看到模型中原子或分子的运动状况。,5. GULP的介绍,性质选项： GULP的性质选项不多，包括晶格，频率，势和电场梯度。 要得到有用的性能，很多时候需要对模型的灵活应用。 事实上，GULP的计算中单点能量或有条件的结构优化可以得到更多的信息。主要在主输出文件*.gout中列出。 大多数的研究是在模型建立上下功夫，通过模型的更变导致的性能的变化，得到有用的结论。,6. 分子动力学计算的应用,聚合物与氧所铝表面的相互作用 构建Al2O3（012)的解理面 驰豫表面结构，内层原子不变，表面层参与驰豫。可以用分子动力学中的几何优化方法。

15、增加表面的面积，并变成3D的薄层结构。 用重复单元构建聚合物。 将聚合物加入到层状结构中。 参数的选取：NVT系统，温度298，平衡时间50ps。 相互作用能计算：总能量减去相互作用能。,6. 分子动力学计算的应用,Cu纳米线、纳米薄膜、单晶材料力学性能的模拟 原子模拟的建立 对于面心立方的晶胞，建立8x8x30的超晶胞，实比分大小为2.89nmx2.89nmx10.83nm。 模型分了上下面端的边界区和中间的驰豫区，边界区3层晶胞，驰豫区24层晶胞。 三种边界条件：x，y方向自由，z方向为周期性边界条件；z方向自由，x，y方向为周期性边界条件；x，y，z方向都为周期性边界条件；分别对应纳米线

16、，纳米薄膜和块体。 模拟过程和参数选取 先沿z方向均匀施加0.3%的拉伸应变，然后驰豫1000步；重复此拉伸，直到材料发生破坏。模拟时采用EAM势(在GULP中可以选取Johnson势，属于EAM势），温度为0K。,6. 分子动力学计算的应用,结果与讨论 上图显示了三种边界条件下铜单晶（或者说纳米线、纳米薄膜、块体）的拉伸应力应变曲线，图中应力纵轴代表铜单晶弛豫区原子的Z向正应力平均值。,拉伸曲线1为纳米线的，开始表现为线性，应变为0.09时应力达到峰值11. 17 GPa后突 然下降，随后发生塑性流动,流动应力基本保持在3.8 GPa附近。 纳米薄膜的应力应变曲线2与曲线1类似，初始为线性，

17、应变为0.1时应力峰值为11. 65 GPa，塑性流动应力在5 GPa左右波动。 铜单晶块体的曲线3表现出完全不同的特征，应变0.05前应力应变曲线基本为线性，此后逐步弯曲；应变超过0.135后，应力发生小的突降，接着继续上升，但斜率不断减小;应 变为0.21时,应力达最大值18.88 GPa;随后应力很快下降到零，未出现类似曲线1和曲线2 的塑性流动。,6. 分子动力学计算的应用,结论 纳米线、纳米薄膜的应力应变曲线及内在变形机理基本类似，但纳米薄膜中位错运动受到更大的阻碍从而强度略高;位错运动是塑性的来源；两者都在应变0.1附近出现位错，随后发生很大的塑性流动。 位错在铜单晶块体中受到极大

18、阻碍，因此块体的强度极高，但破坏时延性相对较小。 纳米器件表面存在明显的表面张力。,6. 分子动力学计算的应用,非晶态形成过程的模拟 由液态急冷而成的非晶态材料具有许多优良性能，这主要是由其特殊的微观结构决定的，而固态微观结构又是由其液态母体微观结构在凝固过程不断转变后得到的。 对非晶态形成过程进行跟踪研究，了解和控制这种微观结构随温度变化的规律，将具有重要的理论和实际意义。 在目前的实验条件下，对于液态金属凝固过程很难实现微观跟踪研究。采用计 算机模拟快速凝固过程，可以获得实验上无法得到的微观结构演变信息，使材料设计和性能预测成为可能。 不同成分的合金可以在什么条件下可形成非晶态？在这方面分

19、子动力学有意想不到的效果，已经模拟了很多单原子液体和合金的急冷过程，得到了非晶材料。 以下以液态金属Ga的快速凝固模拟为例介绍模拟方法。,6. 分子动力学计算的应用,1.模型的建立， 采用分子动力学方法对500个Ga原子的急冷凝固过程中的微观结构组态变化进行模拟研究。 2.参数的选取 采用三维周期性边界条件，原子间相互作用采用有效双体势， 双体势的截止距离为 20 au，运行的时间步长为10-15 s。模拟计算从Ga的熔点(302.80 K) 附近的324 K开始，在该温度下，让系统运行20000次，以确保其处于平衡态。 分别以1. 69x1012 K . s-1和 1.01x1011 K s

20、-1 的冷却速率降温，在分别达到 254 K, 204 K，154 K，134 K, 74 K, 4 K温度点时，再让系统运行4 000步，使其达到平衡态。,6. 分子动力学计算的应用,3.结果分析 (1)平均原子总能量随温度的变化 通过模拟，得到不同冷却速度下的平均原子总能量图。,当冷却速度为1. 69x1012 K s-1的冷却速率降温时，没有发生晶化，E随 T下降线性下降。 在1.01x1011 K _ s1的冷却速度下发生明显晶化现象，在温度约为144 K附近出现能量跳跃下降，说明该过程中体系发生了相转变，使得体系能量突然下降，体系趋 向稳定相结构，该温度为结晶转变温度Tc。,6. 分

21、子动力学计算的应用,(2)双体分布函数随温度的变化,7. 第一性原理分子动力学,经典分子动力学方法的中心思想是利用经验原子间相互作用势计算体系的平衡态和非平衡态的物理性质。 其优点是可用于大块固体、原子团簇、非晶态和液体等物质系统。 对于惰性元素组成的体系，原子间势可用Lennard - Jones势；对于金属体系，提出了原子嵌入的模型势；对于共价晶体有Stillinger-Weber势。 在研究晶体的电子态方面，密度泛函（DF)理论取 得了很大成功。 然而由于计算上的复杂性，基于密度泛函理论的第一性原理计算长期以来是不可能直接用于统计物理的模拟。 如何将这两种方法结合起来，一直是人们的理论追

22、求。,7. 第一性原理分子动力学,电子的质量比原子核的质量轻得多，通常可以假设Born-OppenheimeKBO)绝热近似， 即将电子的运动和原子核的运动之间的耦合忽略，在每个时间步对电子和原子核的运动独立 考虑。 在原子核的运动速度不太快和温度非常低时，假定电子处于基态。“绝热势能面”或 “Born - Oppenheimer面”（BO面）表示在3Natom个坐标空间在Born - Oppenheimer近似下的能量面。在这里必须对BO近似的使用作两点说明。 第一，由轻原子(如氢原子)构成的体系， 在绝热近似下，讨论稳定原子构型的能量差时，原子核的零点振动能变得非常重要，在这种情况下，必须将零点振动能考虑在内。 第二，如果电子不在基态或者原子运动的速度很大时，BO 近似失效。如