分子动力学模拟方法的基本原理与应用

《分子动力学模拟方法的基本原理与应用》由会员分享，可在线阅读，更多相关《分子动力学模拟方法的基本原理与应用（12页珍藏版）》请在装配图网上搜索。

1、word分子动力学模拟方法的根本原理与应用摘要: 介绍了分子动力学模拟的根本原理与常用的原子间相互作用势, 如Lennard-Jones势; 论述了几种常用的有限差分算法, 如Verlet算法; 说明了分子动力学模拟的几种系综与感兴趣的宏观统计量的提取。关键词: 分子动力学模拟; 原子间相互作用势; 有限差分算法; 分子动力学是一门结合物理，数学和化学的综合技术。分子动力学是一套分子模拟方法，该方法主要是依靠牛顿力学来模拟分子体系的运动，以在由分子体系的不同状态构成的系统中抽取样本，从而计算体系的构型积分，并以构型积分的结果为根底进一步计算体系的热力学量和其他宏观性质。从统计物理学中衍生出来的

2、分子动力学模拟方法(Molecular Dynamics Simulation, MDS) , 实践证明是一种描述纳米科技研究对象的有效方法, 得到越来越广泛的重视。所谓分子动力学模拟, 是指对于原子核和电子所构成的多体系统, 用计算机模拟原子核的运动过程, 从而计算系统的结构和性质, 其中每一个原子核被视为在全部其他原子核和电子所提供的经验势场作用下按牛顿定律运动。它被认为是本世纪以来除理论分析和实验观察之外的第三种科学研究手段, 称之为“计算机实验手段, 在物理学、化学、生物学和材料科学等许多领域中得到广泛地应用。科学工作者在长期的科学研究实践中发现，当实验研究方法不能满足研究工作的需求时

3、，用计算机模拟却可以提供实验上尚无法获得或很难获得的重要信息；尽管计算机模拟不能完全取代实验，但可以用来指导实验，并验证某些理论假设，从而促进理论和实验的开展。特别是在材料形成过程中许多与原子有关的微观细节，在实验中根本上是无法获得的，而在计算机模拟中即可以方便地得到。这种优点使分子动力学模拟在材料研究中显得非常有吸引力。分子动力学模拟就是用计算机方法来表示统计力学，作为实验的一个辅助手段。分子模拟就是对于原子核和电子所构成的多体系统，求解运动方程(如牛顿方程、哈密顿方程或拉格朗日方程)，其中每一个原子核被视为在全部其它原子核和电子作用下运动，通过分析系统中各粒子的受力情况，用经典或量子的方法

4、求解系统中各粒子在某时刻的位置和速度，以确定粒子的运动状态，进而计算系统的结构和性质。该模拟技术主要涉与粒子运动的动力学问题,与蒙特卡罗模拟方法(简称MC)相比,分子动力学是一种“确定性方法,它所计算的是时间平均,而MC进展的是系综平均。然而按照统计力学各态历经假设,时间平均等价于系综平均。因此,两种方法严格的比拟计算能给出几乎一样的结果。经典的分子动力学方法是Alder等于1957年提出并首先在“硬球液体模型下应用，发现了由Kirkwood在1939年根据统计力学预言的“刚性球组成的集合系统会发生有液相到结晶相的转变。后来人们称这种相变为Alder相变。Rahman于1963年采用连续势模型

5、研究了液体的分子动力学模拟。1972年Less等开展了该方法并扩展了存在速度梯度的非平衡系统。1980年Andersen等创造了恒压分子动力学方法。1983年Gillan等将该方法推广到具有温度梯度的非平衡系统，从而形成了非平衡系统分子动力学方法体系。1984年Nose等完成了恒温分子动力学方法的创建。1985年针对势函数模型化比拟困难的半导体和金属等，Car等提出了将电子论与分子动力学方法有机统一起来的第一性原理分子动力学方法。1991年Cagin等进一步提出了应用于处理吸附问题的巨正如此系综分子动力学方法。20世纪80年代后期，计算机技术飞速开展，加上多体势函数的提出与开展，使分子动力学模

6、拟技术有了进一步的开展。1、分子动力学的运动方程：分子动力学模拟的出发点是假定粒子的运动可以用经典动力学来处理，对一个由个粒子构成的孤立体系，粒子的运动由牛顿运动方程决定，也就是：mid2ri/dt2=-iV(r1,r2,rN),式中，mi，ri分别为第i个原子的质量和位置。i=-/ri，V(r1,r2,rN)为体系所处的势。2、运动方程的数值积分：计算机模拟方法的基点是利用现代计算机高速和准确的优点 , 对几百个以至上千个分子的运动方程进展数值积分有许多不同的积分方法 , 它们的效率和方便程度各异问题根本上就是用有限差分法来对二阶常微分方程进展积分常用的有以下几种方法：2.1、Verlet算

7、法1：Verlet算法是在60年代后期出现的，对扩散分子的质心运动的积分是最稳定的也是最常用的数值方法。它运用t时刻的位置和加速度以与t时刻的位置来预测tt位置，其积分方案，以三阶Taylor展开为根底，由以下方程给出：r(tt)2r(t)-r(t-t)+t2(t)这里，为简单计，省略了i。速度可按微分的根本法如此得出：V(t)=r(t+t)-r(t-t)/2t。这种算法的优点是占有计算机的存小，并且很容易编程。但它的缺点是位置r(t+t)要通过小项(t2)与非常大的两项2r(t)和r(t-t)的差相加得到，这容易造成精度损失。并且从式中可以看出，这种算法不是一个自启动算法，新位置必须由t和t

8、-,t时刻的位置得到。2.2、Gear的预侧-校正算法：这种算法分为三步来完成：首先，根据Taylor展开，预测新的位置、速度和加速度。然后，根据新的计算的力计算加速度。这个加速度再由与Taylor级数展开式中的加速度进展比拟，两者之差在校正步里用来校正位置和速度项。这种方法的缺点就是占有计算机的存大。2.3、“蛙跳Leap-frog算法2：Hockey提出的Leap-frog算法是Verlet算法的变化，这种方法设计半时间间隔的速度，即：r(t+t)=r(t)+tv(t+t/2)，v(t+t/2)=v(t-t/2)+t(t)。t时刻的速度由下式给出：v(t)=v(t+t/2)+v(t-t/2

9、)/2。这种算法与Verlet算法相比有两个优点：1包括显速度项；2收敛速度快，计算量小。这种算法明显的缺陷是位置和速度不同步。除了上述提与的几种方法外, 还有Beeman算法、Rahman等。3、周期性边界条件和长程力：即使是使用现代的巨型计算机，MD方法还是只能用于粒子数大约是几百到几千的系统。这就引起一个问题：用这样少量的粒子，如何来模拟宏观体系？为了解决这个问题，引入了周期性边界条件3。采用这种方法，模拟体系实际上是由根本单元也称为模拟计算元胞在各个方向上重复叠合而成。但在模拟中只需保存根本单元，所有其它单元与根本单元由平移对称性关联。在处理粒子之间的相互作用时，通常采用“最小影像约定

10、。这个约定是在由无穷重复的MD根本模拟计算元胞中，一个粒子只与它所在的根本元胞的另外N-1个设在此元胞有N个粒子中的每个粒子或其最邻近影像粒子发生相互作用。实际上，这个约定就是通过满足不等式条件rcL/2来截断位势rc为截断半径，L是元胞的边长通常L的数值应当选得很大，以防止有限尺寸效应，但这样会增大计算量，同城采用对相互作用势的修正来近似处理。4、势函数： MD模拟结果准确与否的关键在于对系统的原子之间相互作用势函数的选取，总的来说，原子或分子之间的相互作用势的研究进展一直很缓慢，在一定程度上制约了MD方法在实际研究中的应用，原子间的势函数的开展经历了从对势到多体势的过程，对势认为原子之间的

11、相互作用是两两之间的作用，与其它原子的位置无关，而实际上，在多原子体系中，一个原子的位置不同，将影响其它原子间的有效相互作用。所以，多体势能更准确地表示多原子体系势函数。4.1、对势：在分子动力学模拟的初期，经常采用的就是对势。对势可以分为连续对势1和连续对势，而连续对势主要有以下几种：Lennard-JonesL-J势、Born-landeB-L势、Morse势和Johnson势等，其中，L-J势是为描述惰性气体分子之间的相互作用而建立的，因此它表达的力比拟弱，描述材料的行为也比拟柔韧，也可以用来描述过渡金属原子之间的相互作用。B-L势是用来描述离子晶体离子之间的相互作用的。Morse势和J

12、ohnson势多用于描述金属原子之间的相互作用。对势虽然简单，得到的结果往往也符合某些宏观物理规律，但其缺点是必然导致Cauchy关系，所以，对势实际上不能准确地描述晶体的弹性性质。4.2、多体势：多体势是在20世纪80年代初期开始出现的，1984年Daw和Baskes首次提出了原子嵌入EAM势4。此势的根本思想是：把晶体的总势能分成位于晶格点阵上的原子核之间的相互作用对势和原子核镶嵌在电子云背景中的嵌入能多体相互作用两局部，其中，对势和多体作用势的函数形式往往根据经验选取。基于EAM势的势函数还有很多种，这些多体势大多用于金属的微观模拟，此外，还有许多形式的多体势函数形式，1987年Jaco

13、bsen等人在等效介质原理的根底上提出的另一种多体势函数形式，由于其简单、有效，也得到了广泛地应用。5、分子动力学模拟的系综：平衡态分子动力学模拟是在一定的系综瞎进展的，经常用的平衡系综是NVT或NPT系综。在这两种系综中，牵涉到控制温度和压力的几种技术，分别如下：5.1、控温方法：1速度标定：系统的温度与动能存在如下关系：，式中，N是原子数，Nc是约束数，KB是Boltzmann常数，vi是原子i的速度。由于系统的温度和动能存在这样的关系，所以一种最简单和最直观的方法是直接对速度进展标定。这种方法的根本思想是：如果t时刻的温度是T(t)，速度乘以因子后，温度的变化为T=(2-1)T(t)，其

14、中，T为所控制体系的温度。2Berendsen恒温槽方法：这种控温方法假设所模拟的体系与一个恒温槽连载一起，如此两者之间就可以通过热交换而使模拟体系达到恒温的目的，方法如下，定义一个参数：，式中表征系统与恒温槽之间的热交换速率，t为MD的时间步长，那么通过vnew=vold校正即可保持体系的温度在T0附近震动，而参数ps如此可用于控制这个震动幅度。3Nose-Hoover方法：这种方法是通过改变模拟体系的Hamiltonian来实现控温的，因而有更强的物理意义，其根本思想就是在Hamiltonian参加一个假想的项来代表一个恒温源，具体做法如下：，式中，S和分别是假想项的坐标和动量。这样体系的

15、微分方程就变为：vi=dri/dr，ai=-(dV/dr+mivi)/mi，d。式中，g为体系的自由度，Q为一个可调参量，表征着假想项的质量，T为温度。在这三种控温方法中，速度标定是一种非真实的物理效应。但这种方法可以使系统很快达到平衡，在经典分子动力学方法中，这是一种比拟常用的控温方法。Berendsen控温方法是通过系统和恒温槽进展热交换来控制温度，此方法的优点是非常简单，应用起来非常方便。Nose-Hoover控温方法是基于统计力学而提出来的一种控温方法，当系统和恒温槽进展热交换时，在系统中粒子出现的几率遵从统计力学规律，所以这是一种真实的物理效应。、控压方法: 用MD方法研究压力的诱导

16、相变和结构重构，在等压下模拟比在等体积下更容易实现。在等压模拟下，可以通过改变模拟元胞的三个方向或一个方向的尺寸来实现体积的变化，类似于温度控制的方法，也有许多方法用于压力控制，在实际应用中，有以下两种较为常用的压力控制方法。（1） Berendsen方法：这种控压方法的根本思路与根本控制方程均与他的控温方法类似，如下：，式中是一个可调参数，为所控制的压力。然后，在每一次的MD迭代中，粒子坐标x，y，z均用相乘，得到新的坐标，即可实现对压力P的控制。该方法的缺点是可能导致在很长的模拟时间步长，保持着起伏。 2Parrinello-Rahaman方法：这种方法是通过改变体系的Lagrangian

17、来实现的，新的Lagrangian定义如下：，其中，是模拟单胞的基矢，原子坐标为，V为原子之间的相互作用势，为单胞体积，Tr表示矩阵的迹，W是一个可调参量。这种方法的有点是模拟单胞不仅大小可变，形状也可变。6、热力学性质的计算：经典力学的一个根本前提是，只要知道了物理体系在相空间的运动，就可以导出其所有的宏观性质。在分子动力学中，对模拟得到的位形进展结构分析是很重要的。比如在研究晶体的熔化时，需要与时分析它的结构变化。就结构分析而言，目前常用的方法有：径向分布函数、静态结构因子和配位数等方法。、径向分布函数：径向分布函数可以由下式来定义：，式中，n(r)表示距离原点r到r+r之间的平均粒子数，

18、V是模拟的体积。RDF表示的物理含义是：再空间位置r点周围的体积元中单位体积发现另一个粒子的几率。由此看出，RDF表征着结构的无序化程度。、静态结构因子：静态结构因子表示式为：，其中，N是总原子数，K是倒格矢，为原子位置矢量。对理想晶体而言，SSF为1，而对理想流体，如此为0。、配位数：所谓配位数，就是原子的第一最近邻原子的个数配位数的多少 , 表示了该原子周围的原子分布的密度大小，它经常在一些分析中作为结构分析的一种辅助手段。7、MD方法的进一步研究方向：经典分子动力学方法存在两个缺陷：1元胞体积和形状保持不变，限制MD方法的应用；2不适合含有自由电子的系统，对金属等系统的计算的结果不理想。

19、为克制经典分子动力学方法的局限性，20世纪80年代初Anderson、Parrinello等人开展了可变元胞分子动力学方法。由于能带论中的密度泛函DF方法是研究凝聚态物质电子结构的有效方法，1985年Car等人提出了第一原理的MD方法。但是，进一步拓宽分子动力学方法的应用面显然是一个热门的开展方向。对MD方法本身而言，最重要的两个要素是初始位形的给定和相互作用势函数的选取。MD模拟时间足够长，初始条件的选择不会影响计算的结果，但是选择合理的初始条件可以加快系统趋于平衡，可以节省机时。影响MD计算结果准确程度的最主要的因素是相互作用势的准确性，人们对提高势函数的准确性进展了大量的研究提出了许多势

20、函数的形式，但是，对大多数原子而言，努力找到即准确而形式又不太复杂的势函数，仍然是一个很有挑战性的课题。我们认为，把第一原理、量子化学分析和参数拟合相结合是势函数研究的最好的方法，NanxianChen等人就用此方法研究了氯化钠晶体离子之间的相互作用势，并研究了微结构的变化，这是一个值得进一步研究的方向。58、分子动力学模拟方法的应用：、氨的自扩散系数的分子动力学模拟研究：扩散系数是描述传递现象的重要的流体热物理性质之一，较宽温度和压力围的扩散系数在化工生产中具有重要应用价值。扩散系数的实验测量是一项很困难的工作，而分子动力学MD模拟撒获得分子扩散系数的有效方法。氨，作为大气中最重要的碱性气体

21、之一，它在大气化学、物质传输与沉降过程中发挥着重要作用。采用MD模拟方法计算氨在较宽温度和压力围的自扩散系数，通过将两种不同的控温方式得到的模拟值与大量的实验值进展比拟，验证了所采用的模拟方法是合理的。从而可有效地预测极端条件下实验难以测量的扩散系数。MD模拟采用TINKER v5.0程序包，OPLS-AA全原子力场。模拟体系包含300个分子。分子间的作用势采用Lennard-Jones势，体系中的长程作用力采用Ewald加和的形式。在模拟过程中均采用周期性边界条件、Beeman算法求解运动方程。温控方式分别采用Andersen和Berendsen。体系的截断半径为1nm，时间步长为1fs。M

22、D模拟在NVT系综进展。通过模拟发现：从常温到高温，采用Berendsen控温方式得到的氨的自扩散系数的模拟值与实验值吻合的很好。于是，可以采用分子动力学模拟来代替实验，获得高温高压条件下实验难以测量的分子自扩散系数。另一方面，氨的自扩散系数受温度的影响比受压力的影响更明显。6、半晶态聚合物的分子动力学模拟：，很多聚合物将发生局部区域结晶、局部保存为非晶态的转变，而形成半晶态的结构。实验观测到折叠链结构为聚合物结晶的根本形式，由折叠链片晶可形成多层堆叠片晶、球晶、串晶和枝晶等。针对半晶态聚合物的总体形态，人们提出了多种多样的聚合物结构模型，如缨状微束模型、折叠链模型和折叠链缨状胶粒模型等。采用

23、Meyer和Muler-Plathe提出的粗粒化聚乙烯醇模型，应用分子动力学方法，模拟熔融态聚合物经过缓慢冷却、局部结晶形成半晶态聚合物的过程，研究冷却过程中聚合物的结晶和凝固行为，揭示半晶态聚合物的结构特征与其形成机制。为半晶态聚合物物理、力学性能的分子模拟研究提供根底。应用分子动力学方法 , 模拟熔融态聚合物经过缓慢冷却，局部结晶形成半晶态聚合物的过程。静态结构因子的演变显示出，在结晶行为的初期小角散射强度的增大先于布拉格峰的出现，这个现象与小角/大角x射线散射实验结果相一致。模拟得到了半晶态的聚合物结构，其中的微小晶区由折叠链构成，它们随机取向分布在非晶态区中，总体形态可用缨状微束结构模

24、型描述。模拟得到的半晶态聚合物的结晶度为45%，分子链持久长度即晶区平均厚度约为。通过比拟结晶度和分子链持久长度随温度的变化曲线，发现在不同的冷却阶段具有不同的有序结构形成机制。从随机缠绕的分子链到折叠链晶区，需要两种形式的结构转变：分子链的伸展和伸直分子链之间的平行排列。在从结晶温度到玻璃化温度的凝固转变过程中，存在上述两种形式的结构转变；而在玻璃化温度之后，材料的活性只能允许调整伸直分子链之间的相对排列位置，但由此也导致了聚合物结晶度明显的增大。7、多肽抑制剂抑制淀粉质多肽42构象转换的分子动力学模拟和结合自由能计算：应用分子动力学模拟和结合自由能计算方法研究多肽抑制剂KLVFF、VVIA

25、和LPFFD抑制淀粉质多肽42(A42)构象转换的分子机理。结果明确，三种多肽抑制剂均能够有效抑制A42的二级结构由-螺旋向-折叠的构象转换。另外，多肽抑制剂降低了A42分子的疏水相互作用，减少了多肽分子远距离的接触，有效抑制了A42的疏水塌缩，从而起到稳定其初始构象的作用。这些抑制剂与A42之间的疏水和静电相互作用包括氢键均有利于它们抑制A42的构象转换。此外，抑制剂中的带电氨基酸残基可以增强其和A42之间的静电相互作用包括氢键并降低抑制剂之间的聚集，从而大大增强对A42构象转换的抑制能力。但脯氨酸的引入会破坏多肽的线性结构，从而大大降低其与A42之间的作用力。上述分子模拟的结果揭示了多肽抑

26、制剂KLVFF、VVIA和LPFFD抑制A42构象转换的分子机理，对进一步合理设计A高效短肽抑制剂具有非常重要的理论指导意义。8、超临界水中醋酸锌水解反响的分子动力学模拟：超临界水具有特殊的溶解度、可变的密度、较低的粘度和外表力以与较高的扩散特性，利用超临界水热合成法制备纳米ZnO逐渐成为常用的方法之一。在这种方法中，超临界水SCW，温度T=647K，压力p=22MPa，密度=0.322同时作为反响介质和反响物，含Zn离子的盐如Zn(NO3)2、ZnSO4、Zn(CH3COO)2的水溶液在超临界状态下与碱如KOH或纯水先发生水解反响生成Zn(OH)2然后Zn(OH)2脱水得到纳米级的ZnO颗粒

27、。利用超临界水热合成还成功制备出不同种类和形貌的纳米颗粒，如MgO、TiO、ZrO2、SnO2、ZnS和Fe3O4等。利用分子动力学模拟Zn(CH3COO)2与SCW反响的能量变化，以与产物CH3COOH在SCW的分布特性，为实验上通过控制温度或压力生成纳米ZnO颗粒提供参考。利用分子动力学模拟研究Zn(CH3COO)2在SCW中的平衡特征以与水解反响过程的能量变化，定性分析了水分子在Zn(CH3COO)2团簇部由于静电场的影响而发生解离的可能性，得出以下结论：1Zn(CH3COO)2在SCW中易聚集成团，1个Zn周围平均结合5个CH3COO-和1个H2O满足Zn形成6配位八面体的空间构型的需

28、要。Zn2+在Zn(CH3COO)2团簇部和明确的配位情况有所不同，更多的H2O在Zn(CH3COO)2团簇外表参与Zn2+配位。2在标准状态，Zn(CH3COO)2水解生成Zn(OH)2和CH3COOH的反响焓变为60.4kJmol-1，从热力学的角度看难以发生，但在超临界条件下，模拟发现上述水解反响使得体系势能降低，明确反响容易发生，并且水解反响发生的同时伴随Zn(CH3COO)2团簇结构的改变。水解产物OH-容易进入Zn(CH3COO)2团簇部，富集Zn，而产物CH3COOH扩散进入水相，并倾向分布在水相气相界面。3定性分析处于Zn(CH3COO)2团簇部的水分子周围静电势的分布发现，水

29、分子周围分布有很强的负的静电势，容易使水分子极化变形进而解离，促进Zn(CH3COO)2在SCW中的水解反响。9参考文献：1L.Verlet puterexperimentsclassical fluids. I. Thermodynamical properties of Lennard-JonesmoleculesJ. Phys. Rev. ,1967,159:98-1032R.W.Honeycutt. The potential calculation and some applicationsJ. Methods in putational Physics,1970,9:136-211

30、3N.Metropolis,A.W.Rosenbluth,M.N.Rosenbluth,A.H.Teller,E.Teller.Equation of state calculations by fast puting machinesJ.J.Chem.Phys.,1953,21:1087-10924M.S.Daw,M.I.Baskes.Embedded atom method derivation and application to impurities,surfaces,and other defects in metalsJ.Phys.Rev.B,1984,29:8485-84955 X守鑫，海泉，肖效光，黄海军，分子动力学模拟根本原理和主要技术O.聊城大学学报，2005，181：30-346周昌林，黎多来，东凯，等，氨的自扩散系数的分子动力学模拟究J. 化工，2012，3913：1831847段芳莉，颜世铛，半晶态聚合物的分子动力学模拟A. 计算物理，2012，295：7597658董晓燕，都文婕，夫锋，多肽抑制剂抑制淀粉质多肽42构象转换的分子动力学模拟和结合自由能计算O，物理化学学报，2012，2811，2734-27449王晓娟，志义，志军，超临界水中醋酸锌水解反响的分子动力学模拟O，物理化学学报，2013，291，23-2912 / 12