1、分子动力学模拟方法的基本原理与应用分子动力学模拟方法的基本原理与应用摘要：介绍了分子动力学模拟的基本原理及常用的原子间相互作用势，如Lennard-Jones势；论述 了几种常用的何限差分算法，如Verlet算法；说明了分子动力学模拟的几种系综及感兴趣的宏观 统计量的提取。关键词：分子动力学模拟；原子间相互作用势；何限差分算法；分子动力学是一门结合物理，数学和化学的综合技术。分子动力学是一套分子 模拟方法，该方法主要是依靠牛顿力学来模拟分子体系的运动，以在由分子体系的 不同状态构成的系统中抽取样本，从而计算体系的构型积分，并以构型积分的结果 为基础进一步计算体系的热力学量和其他宏观性质。从统计
2、物理学中衍生出来的分子动力学模拟方法(Molecular Dynamics Simulation, NfDS),实践证明是一种描述纳米科技研究对象的有效方法，得到越来 越广泛的重视。所谓分子动力学模拟，是指对于原子核和电子所构成的多体系统, 用计算机模拟原子核的运动过程，从而计算系统的结构和性质，其中每一个原子核 被视为在全部其他原子核和电子所提供的经验势场作用下按牛顿定律运动。它被认 为是本世纪以来除理论分析和实验观察之外的第三种科学研究手段，称之为“计算 机实验”手段，在物理学、化学、生物学和材料科学等许多领域中得到广泛地应用。科学工作者在长期的科学研究实践中发现，当实验研究方法不能满足研
3、究工作 的需求时，用计算机模拟却可以提供实验上尚无法获得或很难获得的重要信息；尽 管计算机模拟不能完全取代实验，但可以用来指导实验，并验证某些理论假设，从 而促进理论和实验的发展。特别是在材料形成过程中许多与原子有关的微观细节， 在实验中基本上是无法获得的，而在计算机模拟中即可以方便地得到。这种优点使 分子动力学模拟在材料研究中显得非常有吸引力。分子动力学模拟就是用计算机方法来表示统计力学,作为实验的一个辅助手段。 分子模拟就是对于原子核和电子所构成的多体系统，求解运动方程(如牛顿方程、哈 密顿方程或拉格朗日方程)，其中每一个原子核被视为在全部其它原子核和电子作用 下运动，通过分析系统中各粒子
4、的受力情况，用经典或量子的方法求解系统中各粒 子在某时刻的位置和速度，以确定粒子的运动状态，进而计算系统的结构和性质。 该模拟技术主要涉及粒子运动的动力学问题，与蒙特卡罗模拟方法(简称MC)相比, 分子动力学是一种“确定性方法”，它所计算的是时间平均，而MC进行的是系综平均。 然而按照统计力学备态历经假设,时间平均等价于系综平均。因此，两种方法严格的 比较计算能给出几乎相同的结果。经典的分子动力学方法是Alder等于1957年提出并首先在“硬球”液体模型下 应用，发现了由Kirkwood在1939年根据统计力学预言的“刚性球组成的集合系统 会发生有液相到结晶相的转变”。后来人们称这种相变为Al
5、der相变。Rahman于1963 年采用连续势模型研究了液体的分子动力学模拟。1972年Less等发展了该方法并 扩展了存在速度梯度的非平衡系统。1980年Andersen等创造了恒压分子动力学方 法。1983年Gillan等将该方法推广到具有温度梯度的非平衡系统，从而形成了非 平衡系统分子动力学方法体系。1984年Nose等完成了恒温分子动力学方法的创建。 1985年针对势函数模型化比较困难的半导体和金属等，Car等提出了将电子论与分 子动力学方法有机统一起来的第一性原理分子动力学方法。1991年Cagin等进一步 提岀了应用于处理吸附问题的巨正则系综分子动力学方法。20世纪80年代后期，
6、 计算机技术飞速发展，加上多体势函数的提出与发展，使分子动力学模拟技术有了 进一步的发展。1、分子动力学的运动方程：分子动力学模拟的岀发点是假定粒子的运动可以用经典动力学来处理，对一个 由N个粒子构成的孤立体系，粒子的运动由牛顿运动方程决定，也就是： mid2ri/dt2=-VV(rb r2, rx)，式中，nu,分别为第i个原子的质量和位置。 i=-3/5ri, V (ri, r:, rx)为体系所处的势。2、运动方程的数值积分：计算机模拟方法的基点是利用现代计算机高速和精确的优点，对几百个以至 上千个分子的运动方程进行数值积分有许多不同的积分方法，它们的效率和方便 程度各界问题基本上就是用
7、有限差分法来对二阶常微分方程进行积分常用的有以 下几种方法：2. 1、Verlet 算法叫Verlet算法是在60年代后期出现的，对扩散分子的质心运动的积分是最稳定 的也是最常用的数值方法。它运用t时刻的位置和加速度以及t时刻的位置来预测 t +乱位置，其积分方案，以三阶Taylor展开为基础，由以下方程给出：r (t+St) =2r (t)-r (t-6t)+6t a (t)这里，为简单计，省略了i。速度可按微分的基本法则得出:V(t)=r(t+8t)-r(t-8t)/26to 这种算法的优点是占有计算机的存小，并且很容易编程。但它的缺点是位 要通过小项2)与非常大的两项和心的差相加得到，这
8、容易造成精度损失。 并且从式中可以看出，这种算法不是一个自启动算法，新位置必须由时刻的 位置得到。2.2、 Gs啲预侧-校正算法：这种算法分为三步来完成：首先，根据Taybr展开，预测新的位置、速度和加 速度。然后，根据新的计算的力计算加速度。这个加速度再由与Taylor级数展开式 中的加速度进行比较，两者之差在校正步里用来校正位置和速度项。这种方法的缺 点就是占有计算机的存大。2.3、 “蛙跳” (Leap-frog)算法：Hockey提出的Leap-frog算法是Verlet算法的变化，这种方法设计半时间间隔 的速度，即：r(t+8t)=r(t)+ 8tv(t+8t/2), v(t+6t/
9、2)=v(t-8t/2)+8ta(t)o t 时刻的速度 由下式给岀：V（t）= v（t+8t/2）+v（t-8t/2）/2o这种算法与Verlet算法相比有两个优点:（1）包括显速度项；（2）收敛速度快，计算量小。这种算法明显的缺陷是位置和 速度不同步。除了上述提及的几种方法外，还有Beeman算法、Rahman等。3、 周期性边界条件和长程力：即使是使用现代的巨型计算机，MD方法还是只能用于粒子数大约是几百到几 千的系统。这就引起一个问题：用这样少量的粒子，如何来模拟宏观体系？为了解 决这个问题，引入了周期性边界条件。采用这种方法，模拟体系实际上是由基本 单元（也称为模拟计算元胞）在各个方
10、向上重复叠合而成。但在模拟中只需保留基 本单元，所有其它单元与基本单元由平移对称性关联。在处理粒子之间的相互作用时，通常釆用“最小影像”约定。这个约定是在由 无穷重复的MD基本模拟计算元胞中，一个粒子只与它所在的基本元胞的另外N-1 个（设在此元胞有N个粒子）中的每个粒子或其最邻近影像粒子发生相互作用。 实际上，这个约定就是通过满足不等式条件rc Parrinello 等人发展了可变元胞分子动力学方法。由于能带论中的密度泛函（DF）方法是研究 凝聚态物质电子结构的有效方法，1985年Cai等人提出了第一原理的MD方法。但是，进一步拓宽分子动力学方法的应用而显然是一个热门的发展方向。对 MD方法
11、本身而言，最重要的两个要素是初始位形的给定和相互作用势函数的选取。 MD模拟时间足够长，初始条件的选择不会影响计算的结果，但是选择合理的初始 条件可以加快系统趋于平衡，可以节省机时。影响MD计算结果精确程度的最主要的因素是相互作用势的精确性，人们对提 高势函数的精确性进行了大量的研究提出了许多势函数的形式，但是，对大多数原 子而言，努力找到即准确而形式又不太复杂的势函数，仍然是一个很有挑战性的课 题。我们认为，把第一原理、量子化学分析和参数拟合相结合是势函数研究的最好 的方法，NanxianChen等人就用此方法研究了氯化钠晶体离子之间的相互作用势， 并研究了微结构的变化，这是一个值得进一步研
12、究的方向。【58、 分子动力学模拟方法的应用：&1、氨的自扩散系数的分子动力学模拟研究：扩散系数是描述传递现象的重要的流体热物理性质之一，较宽温度和压力围的 扩散系数在化工生产中具有重要应用价值。扩散系数的实验测量是一项很困难的工 作，而分子动力学（MD）模拟撒获得分子扩散系数的有效方法。氨，作为大气中最重要的碱性气体之一，它在大气化学、物质传输与沉降过程 中发挥着重要作用。采用MD模拟方法计算氨在较宽温度和压力围的自扩散系数，通过将两种不同 的控温方式得到的模拟值与大量的实验值进行比较，验证了所釆用的模拟方法是合 理的。从而可有效地预测极端条件下实验难以测量的扩散系数。MD模拟采用TINKE
13、Rv5.0程序包，OPLS-AA全原子力场。模拟体系包含300个分子。分子间的作用势釆用Leiinaid-Jones，体系中的长程 作用力釆用Ewald加和的形式。在模拟过程中均采用周期性边界条件、Beeman算法 求解运动方程。温控方式分别采用Andersen和Beiendsen。体系的截断半径为lnm, 时间步长为lfs。MD模拟在NVT系综进行。通过模拟发现：从常温到高温，釆用Berendsen控温方式得到的氨的自扩散系数 的模拟值与实验值吻合的很好。于是，可以采用分子动力学模拟来代替实验，获得 高温高压条件下实验难以测量的分子自扩散系数。另一方而，氨的自扩散系数受温 度的影响比受压力的
14、影响更明显。8.2、半晶态聚合物的分子动力学模拟：聚合物材料区别于其它材料的一个典型特征是它具有半结晶形态当从熔点以 上温度开始缓慢冷却到室温时，很多聚合物将发生部分区域结晶、部分保留为非晶 态的转变，而形成半晶态的结构。实验观测到折叠链结构为聚合物结品的基本形式, 由折叠链片晶可形成多层堆叠片晶、球晶、串晶和枝晶等。针对半晶态聚合物的总 体形态，人们提岀了多种多样的聚合物结构模型，如缨状微束模型、折叠链模型和 折叠链缨状胶粒模型等。采用Meyer和Muler-Plathe提出的粗粒化聚乙烯醇模型，应用分子动力学方法, 模拟熔融态聚合物经过缓慢冷却、局部结晶形成半晶态聚合物的过程，研究冷却过
15、程中聚合物的结晶和凝固行为，揭示半晶态聚合物的结构特征及其形成机制。为半 晶态聚合物物理、力学性能的分子模拟研究提供基础。应用分子动力学方法，模拟熔融态聚合物经过缓慢冷却，局部结晶形成半晶态 聚合物的过程。静态结构因子的演变显示出，在结晶行为的初期小角散射强度的增 大先于布拉格峰的出现，这个现象与小角/大角x射线散射实验结果相一致。模拟得 到了半晶态的聚合物结构，其中的微小晶区由折叠链构成，它们随机取向分布在非 晶态区中，总体形态可用缨状微束结构模型描述。模拟得到的半晶态聚合物的结晶 度为45%,分子链持久长度（即晶区平均厚度）约为6.5nmo通过比较结晶度和分 子链持久长度随温度的变化I山线
16、，发现在不同的冷却阶段具有不同的有序结构形成 机制。从随机缠绕的分子链到折叠链晶区，需要两种形式的结构转变：分子链的伸 展和伸直分子链之间的平行排列。在从结晶温度到玻璃化温度的凝固转变过程中， 存在上述两种形式的结构转变；而在玻璃化温度之后，材料的活性只能允许调整伸 直分子链之间的相对排列位置，但由此也导致了聚合物结晶度明显的增大。【7】8.3、多肽抑制剂抑制淀粉质多肽42构象转换的分子动力学模拟和结合自由能计算:应用分子动力学模拟和结合自由能计算方法研究多肽抑制剂KLVFF、VVIA和 LPFFD抑制淀粉质多肽42（Ap42）构象转换的分子机理。结果表明，三种多肽抑制剂 均能够有效抑制Ap4
17、2的二级结构由a-螺旋向卜折叠的构象转换。另外，多肽抑制剂降低了 A042分子的疏水相互作用，减少了多肽分子远距离 的接触，有效抑制了 A042的疏水塌缩，从而起到稳定其初始构象的作用。这些抑 制剂与Ap42之间的疏水和静电相互作用（包括氢键）均有利于它们抑制Ap42的 构象转换。此外，抑制剂中的带电氨基酸残基可以增强其和AR42之间的静电相互 作用（包括氢键）并降低抑制剂之间的聚集，从而大大增强对Ap42构象转换的抑 制能力。但脯氨酸的引入会破坏多肽的线性结构，从而大大降低其与A042之间的 作用力。上述分子模拟的结果揭示了多肽抑制剂KLVFF、WIA和LPFFD抑制A042 构象转换的分子
18、机理，对进一步合理设计Ap高效短肽抑制剂具有非常重要的理论 指导意义。8】8.4、超临界水中醋酸锌水解反应的分子动力学模拟：超临界水具有特殊的溶解度、可变的密度、较低的粘度和表面力以及较高的扩 散特性，利用超临界水热合成法制备纳米ZnO逐渐成为常用的方法之一。在这种方 法中，超临界水(SCW,温度T二647K,压力p=22MPa,密度p二0.322匠云！I)同 时作为反应介质和反应物,含Zn离子的盐如Zn(NOs)2、ZnSOa Zn(CHsC00)2的 水溶液在超临界状态下与碱(如KOH)或纯水先发生水解反应生成Zn(OH)2然后 Zn(OH)2脱水得到纳米级的ZnO颗粒。利用超临界水热合成
19、还成功制备岀不同种类 和形貌的纳米颗粒，如MgO. TiO、ZrO2 SnO2 ZnS和FesCU等。利用分子动力学模拟Zn(CH3COO)2与SCW反应的能量变化，以及产物CHsCOOH 在SCW的分布特性，为实验上通过控制温度或压力生成纳米ZnO颗粒提供参考。利用分子动力学模拟研究Zn(CH3COO)2在SCW中的平衡特征以及水解反应过 程的能量变化，定性分析了水分子在Zn(CH3COO)2团簇部由于静电场的影响而发生 解离的可能性，得出以下结论：(1)Zn(CH3COO)2在SCW中易聚集成团，1个Zn周围平均结合5个CHsCOO 和1个H20满足Zn形成6配位八面体的空间构型的需要。Z
20、n2+在Zn(CH3COO)2团 簇部和表明的配位情况有所不同，更多的氏0在Zn(CH3COO)2团簇表而参与Zn2+ 配位。(2)在标准状态，Zn(CH3COO)2水解生成Zn(OH)2和CH3COOH的反应焙变为 60.4kJ.moli,从热力学的角度看难以发生，但在超临界条件下，模拟发现上述水解 反应使得体系势能降低，表明反应容易发生，并且水解反应发生的同时伴随 Zn(CH3COO)2团簇结构的改变。水解产物OH-容易进入Zn(CH3COO)2团簇部，富集 Zn,而产物CHsCOOH扩散进入水相，并倾向分布在水相一气相界而。(3)定性分析处于Zn(CH3COO)2团簇部的水分子周围静电势
21、的分布发现，水 分子周围分布有很强的负的静电势，容易使水分子极化变形进而解离，促进 Zn(CH3COO)2在SCW中的水解反应。【9 参考文献：IlL.Verlet.Computefexperimentsclassical fluids. I. Thermodynamical properties of Lennard-JonesmoleculesJ. Phys. Rev.,1967,159:98-1032 R.W.Honeycutt The potential calculation and some applications J Methods in Computational Physi
22、cs970,9:136211N.MetropolisA.W.RosenbkithM.N.RosenbluthAH.TelleE.Teller.Equation of state calculations by fost computing machinesJJChemPhysJ953,21：108710924M.S.Daw,M丄Baskes.Embedded atom method derivation and application to impurities,surfaces,and other defects in metals J Phys.RcvB1984,29：848584955守鑫，胡海泉，肖效光，黄海军，分子动力学模拟基本原理和主要技术0.聊城 大学学报,2005, 18 (1)： 30-34周昌林，黎多来，东凯，等，氨的自扩散系数的分子动力学模拟究J.化工，2012,39 (13)： 183-1847段芳莉，颜世铛，半晶态聚合物的分子动力学模拟A.计算物理，2012, 29 (5)： 759-7658董晓燕，都文婕，夫锋，多肽抑制剂抑制淀粉质多肽42构象转换的分子动力学 模拟和结合自由能计算0,物理化学学报，2012, 2
