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

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

1、精选优质文档-倾情为你奉上分子动力学模拟方法的基本原理与应用摘要: 介绍了分子动力学模拟的基本原理及常用的原子间相互作用势, 如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、

7、Verlet算法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)的差相加得到，这容易造成精度损失。并且从式中可以看出，这种算法不是一个自启动算法，新

8、位置必须由t和t-,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(

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

10、之间的相互作用时，通常采用“最小影像”约定。这个约定是在由无穷重复的MD基本模拟计算元胞中，一个粒子只与它所在的基本元胞内的另外N-1个（设在此元胞内有N个粒子） 中的每个粒子或其最邻近影像粒子发生相互作用。实际上，这个约定就是通过满足不等式条件rc1Vrij+WTrhh/2-p ，其中，h=a,b,c是模拟单胞的基矢，Ghh，原子坐标ri为rihs，V为原子之间的相互作用势，为单胞体积，Tr表示矩阵的迹，W是一个可调参量。这种方法的有点是模拟单胞不仅大小可变，形状也可变。6、热力学性质的计算： 经典力学的一个基本前提是，只要知道了物理体系在相空间的运动，就可以导出其所有的宏观性质。在分子动力

11、学中，对模拟得到的位形进行结构分析是很重要的。比如在研究晶体的熔化时，需要及时分析它的结构变化。就结构分析而言，目前常用的方法有：径向分布函数、静态结构因子和配位数等方法。6.1、径向分布函数： 径向分布函数可以由下式来定义：gr=VNn(r)/4r2r，式中，n(r)表示距离原点r到r+r之间的平均粒子数，V是模拟的体积。RDF表示的物理含义是：再空间位置r点周围的体积元中单位体积内发现另一个粒子的几率。由此看出，RDF表征着结构的无序化程度。6.2、静态结构因子：静态结构因子表示式为：Sk=|j=1NexpiKrj|/N，其中，N是总原子数，K是倒格矢，rj为原子位置矢量。对理想晶体而言，

12、SSF为1，而对理想流体，则为0。6.3、配位数：所谓配位数，就是原子的第一最近邻原子的个数 配位数的多少 , 表示了该原子周围的原子分布的密度大小，它经常在一些分析中作为结构分析的一种辅助手段。7、MD方法的进一步研究方向：经典分子动力学方法存在两个缺陷：（1）元胞体积和形状保持不变，限制MD方法的应用；（2）不适合含有自由电子的系统，对金属等系统的计算的结果不理想。为克服经典分子动力学方法的局限性，20世纪80年代初Anderson、Parrinello等人发展了可变元胞分子动力学方法。由于能带论中的密度泛函（DF）方法是研究凝聚态物质电子结构的有效方法，1985年Car等人提出了第一原理

13、的MD方法。但是，进一步拓宽分子动力学方法的应用面显然是一个热门的发展方向。对MD方法本身而言，最重要的两个要素是初始位形的给定和相互作用势函数的选取。MD模拟时间足够长，初始条件的选择不会影响计算的结果，但是选择合理的初始条件可以加快系统趋于平衡，可以节省机时。影响MD计算结果精确程度的最主要的因素是相互作用势的精确性，人们对提高势函数的精确性进行了大量的研究提出了许多势函数的形式，但是，对大多数原子而言，努力找到即准确而形式又不太复杂的势函数，仍然是一个很有挑战性的课题。我们认为，把第一原理、量子化学分析和参数拟合相结合是势函数研究的最好的方法，NanxianChen等人就用此方法研究了氯

14、化钠晶体离子之间的相互作用势，并研究了微结构的变化，这是一个值 得进一步研究的方向。58、分子动力学模拟方法的应用：8.1、氨的自扩散系数的分子动力学模拟研究：扩散系数是描述传递现象的重要的流体热物理性质之一，较宽温度和压力范围的扩散系数在化工生产中具有重要应用价值。扩散系数的实验测量是一项很困难的工作，而分子动力学（MD）模拟撒获得分子扩散系数的有效方法。氨，作为大气中最重要的碱性气体之一，它在大气化学、物质传输与沉降过程中发挥着重要作用。采用MD模拟方法计算氨在较宽温度和压力范围的自扩散系数，通过将两种不同的控温方式得到的模拟值与大量的实验值进行比较，验证了所采用的模拟方法是合理的。从而可

15、有效地预测极端条件下实验难以测量的扩散系数。MD模拟采用TINKER v5.0程序包，OPLS-AA全原子力场。模拟体系包含300个分子。分子间的作用势采用Lennard-Jones势，体系中的长程作用力采用Ewald加和的形式。在模拟过程中均采用周期性边界条件、Beeman算法求解运动方程。温控方式分别采用Andersen和Berendsen。体系的截断半径为1nm，时间步长为1fs。MD模拟在NVT系综进行。通过模拟发现：从常温到高温，采用Berendsen控温方式得到的氨的自扩散系数的模拟值与实验值吻合的很好。于是，可以采用分子动力学模拟来代替实验，获得高温高压条件下实验难以测量的分子自

16、扩散系数。另一方面，氨的自扩散系数受温度的影响比受压力的影响更明显。68.2、半晶态聚合物的分子动力学模拟：聚合物材料区别于其它材料的一个典型特征是它具有半结晶形态.当从熔点以上温度开始缓慢冷却到室温时，很多聚合物将发生部分区域结晶、部分保留为非晶态的转变，而形成半晶态的结构。实验观测到折叠链结构为聚合物结晶的基本形式，由折叠链片晶可形成多层堆叠片晶、球晶、串晶和枝晶等。针对半晶态聚合物的总体形态，人们提出了多种多样的聚合物结构模型，如缨状微束模型、折叠链模型和折叠链缨状胶粒模型等。 采用Meyer和Muler-Plathe提出的粗粒化聚乙烯醇模型，应用分子动力学方法，模拟熔融态聚合物经过缓慢

17、冷却、局部结晶形成半晶态聚合物的过程，研究冷却过程中聚合物的结晶和凝固行为，揭示半晶态聚合物的结构特征及其形成机制。为半晶态聚合物物理、力学性能的分子模拟研究提供基础。 应用分子动力学方法 , 模拟熔融态聚合物经过缓慢冷却，局部结晶形成半晶态聚合物的过程。静态结构因子的演变显示出，在结晶行为的初期小角散射强度的增大先于布拉格峰的出现，这个现象与小角/大角x射线散射实验结果相一致。模拟得到了半晶态的聚合物结构，其中的微小晶区由折叠链构成，它们随机取向分布在非晶态区中，总体形态可用缨状微束结构模型描述。模拟得到的半晶态聚合物的结晶度为45%，分子链持久长度（即晶区平均厚度）约为6.5nm。通过比较

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

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

20、一步合理设计A高效短肽抑制剂具有非常重要的理论指导意义。88.4、超临界水中醋酸锌水解反应的分子动力学模拟：超临界水具有特殊的溶解度、可变的密度、较低的粘度和表面张力以及较高的扩散特性，利用超临界水热合成法制备纳米ZnO逐渐成为常用的方法之一。在这种方法中，超临界水（SCW，温度T=647K，压力p=22MPa，密度=0.322g.cm-3）同时作为反应介质和反应物，含Zn离子的盐如Zn(NO3)2 、ZnSO4、Zn(CH3COO)2的水溶液在超临界状态下与碱（如KOH）或纯水先发生水解反应生成Zn(OH)2然后Zn(OH)2脱水得到纳米级的ZnO颗粒。利用超临界水热合成还成功制备出不同种类

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

22、OO)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团簇内部的水分子周围静电势的分布发现，水分子周围分布有

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

24、N.Metropolis,A.W.Rosenbluth,M.N.Rosenbluth,A.H.Teller,E.Teller.Equation of state calculations by fast computing 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-8495

25、5 崔守鑫，胡海泉，肖效光，黄海军，分子动力学模拟基本原理和主要技术O.聊城大学学报，2005，18（1）：30-346周昌林，黎多来，李东凯，等，氨的自扩散系数的分子动力学模拟究J. 广东化工，2012，39（13）：1831847段芳莉，颜世铛，半晶态聚合物的分子动力学模拟A. 计算物理，2012，29（5）：7597658董晓燕，都文婕，刘夫锋，多肽抑制剂抑制淀粉质多肽42构象转换的分子动力学模拟和结合自由能计算O，物理化学学报，2012，28（11），2734-27449王晓娟，李志义，刘志军，超临界水中醋酸锌水解反应的分子动力学模拟O，物理化学学报，2013，29（1），23-29专心-专注-专业