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

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

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

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

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

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

6、 机技术飞速发展，加上多体势函数的提出与发展，使分子动力学模拟技术有了进一 步的发展。1、分子动力学的运动方程：分子动力学模拟的出发点是假定粒子的运动可以用经典动力学来处理，对一个由N 个粒子构成的孤立体系，粒子的运动由牛顿运动方程决定，也就是：mid2ri/dt2=- r2, rf),式中，nu, “分别为第i个原子的质量和位置。产-0/3“, V(m r2, rf) 为体系所处的势。2、运动方程的数值积分：计算机模拟方法的基点是利用现代计算机高速和精确的优点，对几百个以至上千 个分子的运动方程进行数值积分有许务不同的积分方法，它们的效率和方便程度 各界问题基本上就是用河限差分法来对二阶常微

7、分方程进行积分常用的有以卜几种 方法：2 1、Verlot算法叫Verlet算法是在60年代后期岀现的，对扩散分子的质心运动的积分是最稳定的也 是最常用的数值方法。它运用t时刻的位置和加速度以及t时刻的位置來预测t+8t 位置，其积分方案，以三阶Taylor展开为基础，由以下方程给岀：r(t+8t) =2r(t)-r(t-8t)+8t2a (t)这里，为简单计，省略Ti速度可按微分的基本法则得出:V(t)=r(t+5t)-r(t-8t)/25t0 这种算法的优点是占有计算机的内存小，并且很容易编程。但它的缺点是位Br(t+8t) 要通过小项2)与非常人的两项2r(t)和心&)的差相加得到，这容

8、易造成精度损失。 并且从式中可以看出，这种算法不是一个自启动算法，新位置必须由t和tdt时刻的 位置得到。2.2、 的预侧-校正算法：这种算法分为三步来完成：首先，根据Taylor展开，预测新的位置、速度和加 速度。然后，根据新的计算的力计算加速度。这个加速度再由与Taylor级数展开式 中的加速度进行比较，两者之差在校正步里用來校正位置和速度项。这种方法的缺 点就是占有计篦机的内存大。2.3、 “蛙跳”(Leap-frog)算法：Hockey提出的Leap-frog算法是Verlet算法的变化，这种方法设计半时间间隔 的速度,即：r(t+8t)=r(t)+8tv(t+5t/2), v(t+8

9、t/2)=v(t-8t/2)+5ta(t) t 时刻的速度 由下式给Hl：v（t）=v（t+8t/2）+v（t-5t/2）/2。这种算法与Verier算法相比冇两个优点：（1）包括显速度项；（2）收敛速度快，计算最小。这种算法明显的缺陷是位置和 速度不同步。除了上述提及的儿种方法外，还右 Beeman隽法、Rahman等。3、 周期性边界条件和长程力：即使是便用现代的巨型计SZ机，MD方法还星只能用于粒子数大约是几百到几 千的系统。这就引起一个问题：用这样少量的粒子，如何来模拟宏观体系？为了解 决这个问题，引入了周期性边界条件卩】。采用这种方法，模拟体系实际上是由基本 单元（也称为模拟计算元胞

10、）在各个方向上重复叠合而成。但在模拟中只需保留基 本单元，所有其它单元与基本单元由平移对称性关联。在处理粒子之间的相互作用时，通常采用“最小影像”约定。这个约定是在由 无穷重复的MD基本模拟计算元胞中，一个粒子只与它所在的基本元胞内的另外N-1 个（设在此元胞内有N个粒子）中的每个粒子或其垠邻近影像粒子发生相互作用。 实际上，这个约定就是通过满足不等式条件rc1(中)+ WTrh h/2 pO.其中，/2 = a,b,c是模拟单胞的基矢，G=hh,原子坐标耳为r, =hs, V为原子之间 的相互作用势，Q为单胞体积，Tr表示矩阵的迹，W是一个可调参量。这种方法的 有点是模拟单胞不仅大小可变，形

11、状也可变。6、热力学性质的计算：经典力学的一个基本前提是，只要知道了物理体系在相空间的运动，就可以导出其 所行的宏观性质。在分子动力学中，对模拟得到的位形进行结构分析是很巫要的。 比如在研究晶体的熔化时，需要及时分析它的结构变化。就结构分析而言，目前常 用的方法育：径向分布函数、静态结构因子和配位数等方法。6.1、 径向分布函数：径向分布函数可以由下式来定义：g(r) = n(r)/4irr2Ar,式中，n(r)表示距离原点 N到之间的平均粒子数，是模拟的体积。RDF表示的物理含义是：再空间位置点 周围的体积元中单位体积内发现另一个粒子的几率。由此看出，RDF表征着结构的无序化 程度。6.2、

12、 静态结构因子：静态结构因子表示式为：S(k) = |Mexp(iK7j)|/N,其中，N是总原子数，K 是倒格矢，z；为原子位置矢量。对理想晶体而言，SSF为1,而对理想流体，则为0。6.3、 配位数:所谓配位数，就是原子的第一最近邻原子的个数配位数的多少，表示了该原子周围的原子分布的密度大小，它经常在一些分析中作为结构分析的一种辅助手段。7.MD方法的进一步研究方向：经典分子动力学方法存在两个缺陷：（1）元胞体积和形状保持不变，限制MD 方法的应用：（2）不适合含有自由电子的系统，对金属等系统的计算的结果不理想。为克服经典分子动力学方法的局限性，20世纪80年代初Anderson、Parr

13、inello 等人发展了可变元胞分子动力学方法由于能带论中的密度泛函（DF）方決是研究 凝聚态物质电子结构的有效方法，1985年Cai等人提出了第一原理的MD方法。但是，进一步拓宽分子动力学方法的应用而显然是一个热门的发展方向。对 MD方法本身而言，最重要的两个要素是初始位形的给定和相互作用势函数的选取。 MD模拟时间足够长，初始条件的选择不会影响计算的结果，但是选择合理的初始 条件口J以加快系统趋于平衡，町以节省机时。影响MD计算结果精确程度的垠主要的因索是相互作用坍的精确性，人们对提 高势函数的精确性进行了人量的研究提出了许多势函数的形式，但是，对人多数原 子而言，努力找到即准确而形式又不

14、太复杂的势函数，仍然是一个很有挑战性的课 题。我们认为，把第一原理、量子化学分析和参数拟合相结合是势函数研究的最好 的方法，NaiixianChen等人就用此方法研究了氯化钠晶体离子之间的相互作用势， 并研究了微结构的变化，这是一个值得进一步研究的方向。P8、 分子动力学模拟方法的应用：8.1、氨的自扩散系数的分子动力学模拟研究：扩散系数是描述传递现象的巫要的流体热物理性质么一，较宽温度和压力范圉 的扩散系数在化工生产中具右重要应用价值。扩散系数的实验测量是一项很困难的 工作，而分子动力学（MD）模拟撒获得分子扩故系数的有效方法。氨，作为大气中最重要的碱性气体之一，它在大气化学、物质传输与沉降

15、过程 中发挥着重要作用。采用MD模拟方法计算氨在较宽温度和压力范圉的自扩散系数，通过将两种不 同的控温方式得到的模拟值与大最的实验值进行比较，验证了所采用的模拟方法是 合理的。从而可有效地预测极端条件卜实验难以测量的扩散系数。MD模拟采用TINKERv5.0程序包，OPLSAA全原子力场。模拟体系包含300个分子。分子间的作用势采用LennaidJones势，体系中的长程 作用力采用Ewald加和的形式祚模拟过程中均采用周期性边界条件、Beeman法 求解运动方程。温控方式分別采用Andersen和Berendsen。体系的截断半径为limi, 时间步长为lfs。MD模拟在NVT系综进行。通过

16、模拟发现：从常温到髙温，采用Berendsen控温方式得到的氨的自扩散系数 的模拟值与实验值吻合的很好。于是，可以釆用分子动力学模拟来代替实验，获得 高温高压条件卜实验难以测量的分子自扩散系数。另一方面，氨的自扩散系数受温 度的影响比受压力的影响更明显。8.2、半晶态聚合物的分子动力学模拟：聚合物材料区别于其它材料的一个典型特征是它具有半结晶形态当从熔点以 上温度开始缓慢冷却到室温时，很多聚合物将发生部分区域结晶、部分保留为非晶 态的转变,而形成半晶态的结构。实验观测到折叠链结构为聚合物结晶的基本形式， 由折穂链片晶可形成多层堆穂片晶、球晶、串晶和枝晶等。针对半晶态聚合物的总 体形态，人们提出

17、了多种多样的聚合物结构模型，如缨状微束模型、折亮链模型利 折叠链缨状胶粒模型等。采用Meyer和MuleP4the提出的粗粒化聚乙烯醇模型，应用分子动力学方法，模 拟熔融态聚合物经过缓慢冷却、局部结晶形成半晶态聚合物的过程，研究冷却过程 中聚合物的结晶和凝固行为，揭示半晶态聚合物的结构特征及其形成机制。为半晶 态聚合物物理、力学性能的分子模拟研究提供基础。应用分子动力学方法，模拟熔融态聚合物经过缓慢冷却，局部结晶形成半晶态聚合 物的过程。静态结构因子的演变显示出，在结晶行为的初期小角散射强度的增大先 于布拉格峰的出现，这个现象与小角/大角X射线散射实验结果相一致。模拟得到了 半品态的聚合物结构

18、，其中的微小晶区由折叠链构成，它们随机取向分布在非晶态 区中，总体形态可用缨状微束结构模型描述。模拟得到的半晶态聚合物的结晶度为 45%,分子链持久长度（即晶区平均厚度）约为6.5nm。通过比较结晶度和分子链 持久长度随温度的变化曲线发现在不同的冷却阶段具育不同的仃序结构形成机制C 从随机缠绕的分子链到折叠链晶区，需要两种形式的结构转变：分子链的伸展和伸 直分子链之间的平行排列。在从结晶温度到玻璃化温度的凝固转变过程中，存在上 述两种形式的结构转变；而在玻璃化温度之后，材料的活性只能允许调整伸直分子 链之间的相对排列位置，但由此也导致了聚合物结晶度明显的增大。卩】8.3、 多肽抑制剂抑制淀粉质

19、多肽42构象转换的分子动力学模拟和结合自由能计算:应用分子动力学模拟和结合自由能计算方法研究务肽抑制剂KLVFF、VVIA利 LPFFD抑制淀粉质务肽42（Ag42）构象转换的分子机理。结果表明，三种多肽抑制剂 均能够有效抑制Ap42的二级结构由8螺旋向P折叠的构象转换。另外，多肽抑制剂降低了 Ap42分子内的疏水相互作用，减少了多肽分子内远距离 的接触，有效抑制了 Ap42的疏水塌缩，从而起到稳定其初始构象的作用。这些抑 制剂与Ap42之间的疏水和静电相互作用（包括氢键）均有利于它们抑制A卩42的 构象转换。此外，抑制剂中的带电氨基酸残基町以增强其和AR42之间的静电相互 作用（包括氢键）并

20、降低抑制剂之间的聚集，从而大大增强对Ap42构象转换的抑 制能力。但脯氨酸的引入会破坏多肽的线性结构，从而人人降低其与A042之间的 作用力。上述分子模拟的结果揭示了多肽抑制剂KLVFF. VVIA和LPFFD抑制A|342 构象转换的分子机理，对进一步合理设计AB高效短肽抑制剂具有非常重要的理论 指导意义。同8.4、 超临界水中醋酸锌水解反应的分子动力学模拟：超临界水具冇特殊的溶解度、可变的密度、较低的粘度和表面张力以及较高的 扩散特性，利用超临界水热合成法制备纳米ZnO逐渐成为常用的方法之一。在这种 方法中，超临界水(SCW,温度T二647K,压力p二22MPa,密度p二T?) 同时作为反

21、应介质和反应物，含Zn离子的盐如Zn(NO3)2、ZnS04 Zn(CH3COO)2 的水溶液在超临界状态下与碱(如KOH)或纯水先发生水解反应生成Zn(OH)2然后 Zn(OH)2脱水得到纳米级的颗料利用超临界水热合成还成功制备出不同种类 和形貌的纳米颗粒，如MgO、TiO、ZrO2 SnO2 ZnS和FesO斗等。利用分子动力学模拟Zn(CH3COO)2与SCW反应的能量变化，以及产物CHsCOOH在 SCW的分布特性，为实验上通过控制温度或压力生成纳米ZnO颗粒提供参考。 利用分子动力学模拟研究Zn(CH3COO)2在SCW中的平衡特征以及水解反应过程的 能量变化，定性分析了水分子在Zn

22、(CH3COO)2团簇内部由于静电场的影响而发生解 离的可能性，得出以卜结论：(1)Zn(CH3COO)2在SCW中易聚集成团，1个Zn周围平均结合5个CH3COO和1 个H20满足Zn形成6配位八面体的空间构型的碍要。Zi0在Zn(CH3COO)2团簇内 部和农明的配位情况有所不同，更多的氏0在Zn(CH3COO)2团簇衷而参巧ZM配位。(2)在标准状态，Zn(CH3COO)2水解生成Zn(OH)2和CH3COOH的反应焙变为 60.4kJ.moli,从热力学的角度看难以发生，但在超临界条件下，模拟发现上述水解 反应使得体系势能降低，表明反应容易发生，并且水解反应发生的同时伴随 Zn(CH3

23、COO)2团簇结构的改变。水解产物0H-容易进入Zn(CH3COO)2团簇内部，富 集Zn,而产物CH3COOH扩散进入水相，并倾向分布在水和一气相界面。(3)定性分析处于Zn(CH3COO)2团簇内部的水分子周国静电势的分布发现， 水分子周田分布有很强的负的静电势，容易使水分子极化变形进而解离，促进 Zn(CH3COO)2在SCW中的水解反应。叨参考文献：1LVerlet.Coniputer experiment s classical fluids. I Therniodynamical properties ofLennard-JonesmoleculesJ. Phys. Rev; ,1

24、967,159:98-1032RW. Honeycutt. The potential calculation and some applications J. Methods in Computational Physics,1970,9： 136-2113N.MetropolisAWRosenbluth,M.N.RosenbluthfA.H.Tellei;E.Teller.Equation of state calculations by fast computing niachinesQJ.J.ChenLPliys.,1953,21:1087-10924M.S.Daw;M丄Baskes.Embedded atom method derivation and application to impui itiessuifaces.and other defects in metalsJ.Phys.RevB,1984,29:8485 -84955崔守鑫，胡海泉，肖效光，黄海军，分子动力学模拟基本原理和主要技术0.聊城大学学报，2005. 18 (1)： 30-346周昌林，黎多来，李东凯，等，氨的口扩散系数的分子动力学模拟究J.广东化工，2012, 39 (13)： 183-1847段芳莉，颜世铛，半晶态聚合物的分子动力学模拟A.计算物理，2012,29(5)： 759