分子动力学模拟橡胶

《分子动力学模拟橡胶》由会员分享，可在线阅读，更多相关《分子动力学模拟橡胶（46页珍藏版）》请在装配图网上搜索。

1、装订线 二甲醚在橡胶中扩散特性的分子动力学模拟 摘 要随着计算机的性能提高和计算技术的长足发展，利用计算机对物质的微观结构和运动进行数值模拟计算得到了迅速发展，并在此基础上发展了运用数值运算统计求和的方法:分子动力学法。本文通过分子动力学方法，应用Accelrys公司的软件Materials Studio来模拟二甲醚在不同橡胶中的扩散行为，用Einstein法和Green-kubo法求得二甲醚在6中橡胶不同温度下的扩散系数，并对两种算法进行比较。基于COMPASS力场进行建模。为了使模拟体系接近真实体系，模型的构建采用三维周期性边界条件，选取能量最低、体系最稳定的周期边界作为初始构型。然后对模

2、拟体系进行能量优化和预平衡，进行NVT系综的分子动力学模拟以便得到二甲醚在橡胶内的扩散过程和运动轨迹，最后采用两种算法计算二甲醚在橡胶内的扩散系数。通过本研究的计算数据，发现分别通过两种算法得到的二甲醚在橡胶中的扩散系数在数量级上保持一致，除去个别点外，两者的相对误差在一定稳定范围内。Green-Kubo法计算需要积分，较为繁琐，误差较大，Einstein法计算简便，计算结果与实验值比较接近。用Green-Kubo法求得二甲醚在343k温度下的氟橡胶中的模拟扩散系数最大，可推测此温度下在所有研究的橡胶中氟橡胶的耐二甲醚的能力最差。基于本文研究得到的数据，进行综合分析后认为本文所研究的橡胶中丙烯

3、腈含量为20%的氢化丁腈橡胶的耐二甲醚性能最好。关键词： 二甲醚，分子动力学，橡胶，扩散系数 ABSTRACTWith the rapid development of computer performance and computing technology, the use of computers and the movement of material microstructure numerical simulation have been developing rapidly, and the development of statistical summation of nume

4、rical calculation method based on this: Molecular Dynamics. In the article, use Materials Studio of Accelrys Company to simulation DME diffusion behavior in the different rubber, and obtain the diffusion coefficient of DME in 6 types of rubber that has different temperature by using Einstein and Gre

5、en-kubo algorithms, and compare the two algorithms.Create model based on the COMPASS force field .In order to simulate the system close to the real system, build a model of three-dimensional periodic boundary conditions, and select the lowest energy, the most stable system with periodic boundary as

6、the initial configuration. Then the simulation system for energy optimization and pre-equilibrium, molecular dynamics simulations with NVT ensemble in order to get DME diffusion processes and trajectories in the rubber, finally, calculate DME diffusion coefficient in the rubber with two algorithms.

7、After analyzing the data of this study, found that the DME diffusion coefficients in the rubber were obtained by the two algorithms consistent in magnitude, except from a few points, the relative error of the two algorithms within a stable range. Green-Kubo method requires integration, more cumberso

8、me, error is large, Einstein method is simple, and the results are close to the experimental value.Under the 343k ,simulation of DME diffusion coefficient in fluorine rubber obtained with green-kubo method is the biggest of the other rubber in the study ,it can be speculated that the ability resista

9、nting DME of fluorine rubber is the worst of the other rubber in the rubber under the temperature .In this paper ,based on the data obtained,after a comprehensive analysis，i think that the ability resistanting DME of HNBR20 is the best of the other rubber studied in this paper.KEY WORDS: Dimethyl et

10、her, Molecular dynamics, Rubber，Diffusion coefficient 目录第一章 绪论11.1分子模拟方法的分类11.1.1量子力学法11.1.2分子力学法11.1.3蒙特卡洛(MC)法21.1.4分子动力学(MD)法21.2分子动力学的发展31.3本文的研究意义51.4本文研究内容6第二章 分子动力学方法介绍82.1模拟软件Materials Studio82.1.1 Materials Studio软件简介82.1.2 使用模块简介82.2分子动力学基本原理92.2.1牛顿运动方程和数值积分算法92.2.2分子动力学模拟的基本步骤112.3分子动力学模

11、拟的关键技术122.3.1系综122.3.2初始条件122.3.3势函数122.3.4常见力场142.3.5周期性边界条件与最近镜像理论152.3.6截断距离与电荷基团152.3.7长程力的计算162.3.8调温控制172.3.9能量最小化192.4扩散系数的计算20第三章 分子动力学模拟过程223.1建立二甲醚与橡胶分子模型223.2构建二甲醚与橡胶分子的晶格体系243.3进行分子动力学模拟273.4模拟分析计算均方位移323.5模拟分析计算速度自相关函数33第四章 模拟结果计算与分析354.1Einstein方法模拟结果计算与分析354.2Green-Kubo方法模拟结果计算与分析374.

12、3对两种算法进行比较39结论及展望41致谢42参考文献43II 第一章 绪论1.1分子模拟方法的分类 分子模拟方法主要包括量子力学(QM)法、分子力学(MM)法、蒙特卡洛(MC)法和分子动力学(MD)法。1.1.1量子力学法量子力学(Quantum Mechanics)是描写微观物质的一个物理学理论,和相对论一起被认为是现代物理学的两大基本支柱，是研究微观粒子的运动规律的物理学分支学科，它主要研究原子、分子、凝聚态物质，以及原子核和基本粒子的结构、性质的基础理论，它和相对论一起构成了现代物理学的理论基础。量子力学不仅是近代物理学的基础理论之一，而且在化学等有关学科和许多近代技术中也得到了广泛的

13、应用。量子力学的基木原理包括量子态的概念，运动方程、理论概念和观测物理量之间的对应规则和物理原理。在量子力学中，一个物理体系的状态由波函数表示，波函数的任意线性叠加仍然代表体系的一种可能状态。状态随时间的变化遵循一个线性微分方程，这个方程预言体系的行为，物理量由满足一定条件的、代表某种运算的算符表示;测量处于某一状态的物理体系的某一物理量的操作，对应于代表该量的算符对其波函数的作用;测量的可能取值由该算符的木征方程决定，测量的期待值由一个包含该算符的积分方程计算。波函数的平方代表作为其变数的物理量出现的几率。根据这些基本原理并附以其他必要的假设，量子力学能够解释原子和亚原子的各种现象。随着计算

14、机硬件的发展，量子力学法已可用于大分子包括聚合物的低聚物在内的体系，尽管如此，将量子力学法应用在多达数千的聚合物体系仍然是不现实的。1.1.2分子力学法分子力学(MM)法，也叫做Force Field Methed，就是忽略了电子运动而只计算与原子位置相关的体系能量，是一种在原子分子水平上用力场解决问题的非量子力学方法，目前广泛地用于计算分子的构象和能量。分子力学从本质上说上是能量最小值方法，就是在原子间相互作用势的作用下，通过改变粒子分布的几何位型，以能量最小为判据，从而获得体系的最佳结构。分子力学主要应用分子力场方法计算分子的势能，分子力学使用解析经验势能函数来描述分子的势能面。每个分子都

15、有固定的力场。在分子力学实际计算时,将力场分解成不同的组分，使用理论计算和实验拟合的方法建立力场参数，力场参数要有可移植性，它要适应于同类分子的计算，即同一类分子也要有很高的计算精度。分子力学认为分子体系的势能函数是分子体系中原子位置的函数。分子力学将分子体系作为在势能面上运动的力学体系来处理，求解的是经典力学方程,而不是量子力学的薛定愕方程。所以分子力学方法可以求得分子的平衡结构和热力学性质等,但不能得到子体系与电子结构有关的其他性质。现在，分子力学是模拟蛋白质、核酸等生物大分子结构和性质以及配体和受体相互作用的常用方法，在生命科学领域得到了广泛的应用。随着分子图形学的不断发展,分子力学已经

16、广泛应用于分子模型设计。当今优秀的分子设计程序都将分子力学作为初始模型优化的主要方法，分子力学和分子图形学己经充分地揉合在分子设计中，分子模型的构建也是分子力学为主，分子力学方法是九三级辅助分子设计中常用的方法,尤其是在药物设计中，已经离不开分子力学计算和模拟方法。应用分子力学方法，可以非常迅速得到分子的低能构象，通过构象分析可以获得合理的药效构象和药效集团。如果己知受体的三维结构，可以用分子力学模拟药物与受体的相互作用。在分子的定量结构活性关系研究中，也需要用分子力学方法进行计算。1.1.3蒙特卡洛(MC)法蒙特卡洛(Monte Carlo)方法，或称计算机随机模拟方法或随机抽样方法或统计试

17、验方法，属于计算数学的一个分支，是一种基于“随机数”的计算方法。Monte Carlo方法也称为统计模拟方法，是二十世纪四十年代中期由于科学技术的发展和电子计算机的发明，而被提出的一种以概率统计理论为指导的一类非常重要的数值计算方法。Monte Carlo方法可以很好地用来对付维数的灾难，因为这种方法的计算复杂性不再依赖于维数，以前那些本来是无法计算的问题现在也能够计算。科学家们提出了许多“方差缩减”技巧来提高方法的效率。Monte Carlo模拟适用于研究复杂体系，研究具有多得数不清的结构、状态的体系，对此我们可以采用蒙特卡洛模拟，以统计的方法寻找出现几率最高的结构、状态，或相应的有关数据。

18、Monte Carlo方法作为一种概率性统计方法在相空间中形成Markov链，每一尝试步移动的结果只依赖于上一步，也就是相空间中的随机行走。它局限于平衡态热力学量的计算，一般不能预测体系的动力学特性，平衡态物理量通过系综平均得到。 1.1.4分子动力学(MD)法自五十年代中期开始，分子动力学方法得到了广泛的应用。它和蒙特卡洛方法一起己经成为计算机模拟的重要方法。应用分子动力学方法取得了很多重要成果，例如气体或液体的状态方程、相变问题、吸附问题等，以及非平衡过程的研究。它的应用已经从化学反应、生物学的蛋白质到重离子碰撞等广泛的学科研究领域。分子动力学是在原子、分子水平上求解多体问题的重要的计算机

19、模拟方法。通过求解所有粒子的运动方程，分子动力学方法可以用于模拟与粒子运动路径相关的基本过程。在分子动力学中，粒子的运动行为是通过经典的牛顿运动方程所描述。1.2分子动力学的发展1956年，Alder和Wainwright1通过研究硬球的集聚首次报道了正确的分子动力学模拟，采用分子动力学研究气体和液体的状态方程，结果显示这种方法在解决实际问题中能够起到重要的作用，从而开创了利用分子动力学模拟方法研究物质宏观性质的先例。后来，人们对这一方法作了许多改进，并运用它对固体及其缺陷以及液体作了大量的研究。真实材料模型的首次分子动力学模拟在1959年被报道了(发表于1960年)，该研究由位于Brookh

20、aven的Vilneyard2领导的小组开展，模拟了钻晶体的辐射危害。1963年Rahman3采用连续势模型对真实流体(氢)进行了首次分子动力学模拟。1967年Verlet4给出了著名的Verlet算法，就是在分子动力学模拟中对粒子运动的距离、速度和加速度的逐步计算法，这种算法后来被广泛采用，为分子动力学模拟计算做出了很大贡献。1973年Rahman和Stillinger5对水(液态多原子分子)做了第一原理分子动力学模拟计算。1980年，Andersen6做了恒压状态下的分子动力学研究，提出了等压分子动力学模型；在同一年，Hoover7进行了非平衡态分子动力学的研究。1981年，Parrine

21、llo和Rahman8给出了恒定压强的分子动力学模型，他们将等压分子动力学推广到允许元胞的形状可以随其中粒子的运动而发生改变的范围，为之后的发展做出了里程碑的贡献。 1983-1984年间，Daw和Baskes9，10根据准原子概念和密度泛函理论，提出了一个能较好描述过渡金属各种性质的新方法一镶嵌原子法(Embedded Atom Method，简称EAM)。后来经过Johnson11的改进，EAM法在固态材料的诸多领域己经取得了很大成功12-15。 1984年，Nose提出恒温分子动力学方法14。1985年，Car和Parrinello16又提出了共价键系统的分子动力学从头计算法，此算法需在

22、高性能计算机上进行。1986年，Frisch, Hasslacher和Pomeau17提出晶格气体自动学。1989年，Frenkel和Emst提出晶格气体自动学和长时尾的处理方法。因为计算机速度运行速度慢以及内存不够大，早期模拟的空间尺度和时间尺度都受到很大限制。多体势函数的提出和发展，为分子动力学模拟技术注入了新的活力。分子动为学模拟不但能得到原子的活动细节，还能像做实验一样进行各种观察。对于平衡系统，可以用分子动力学模拟作适当的时间平均来计算一个物理量的统计平均值。对于非平衡系统，发生在一个分子动力学观察时间内的物理现象也可以用分子动力学计算进行直接模拟。特别是许多在实际实验中无法获得的微

23、观细节，而在分子动力学模拟中都可以方便地观察得到。随着计算机技术的飞速发展和分子模拟技术的不断进步，利用计算机模拟粒子在分子、原子微观水平上的运动，并对粒子运动轨迹进行分析，可以计算出粒子的扩散系数。用分子模拟方法求取扩散系数，一方面节省了大量的时间和成本，另一方面能够模拟在极端环境，比如在高温高压以及临界条件下的扩散系数。对扩散过程的分子模拟研究，国内外己经有很多报道。周健等18,19运用分子动力学方法模拟计算了甲烷、氧气等气体在常温水中的扩散系数，其模拟计算结果与实验值吻合较好，验证了模型的准确性。他们用同样的方法计算了在不同温度下水的自扩散系数，除了在低温条件下模拟值与实验值差异较大外，

24、其它温度下模拟的扩散系数值与实验值吻合较好。赵长伟等20测定了298.15K328.15K不同温度下葡萄糖水溶液的积分扩散系数，考察了温度对扩散系数的影响，提出一个有温度关联项的半经验模型。朱春英等21以Fick定律和不可逆过程热力学为基础，提出了氨基酸在水溶液中扩散系数模型，并利用该模型对多种氨基酸及尿素在水溶液中的扩散系数进行了计算。孙炜等22采用分子动力学模拟方法，计算了L-J流体氩及氩/氪溶液的自扩散系数和互扩散系数，并研究了扩散系数与温度的关系。在扩散系数的数据处理方面，陈敏伯等23 针对模拟计算的均方位移曲线，在原来Einstein方法的基础上有所改进，提出了微分一曲线变分法以及线

25、性回归-区限变分法，能够有效的剔除无效数据，使计算结果更加准确。Pavel等24通过分子动力学模拟计算氧气和二氧化碳在无定形聚对苯二甲酸乙二酯、芳环聚酯及其共混物中的扩散过程，考察了温度、高分子动态性质、密度和自由体积分布对扩散系数的影响。Karayiannis等25模拟小分子在PEI和PET中的扩散过程，研究发现透过物分子在高分子中的扩散主要取决于能容纳渗透物分子的自由体积及其分布。小分子在聚合物或其它多孔介质内的扩散较之本体扩散要复杂得多，这是当前理论研究的热点之一。有关小分子在聚合物中扩散的模拟研究方面，黄宇等对聚二甲基硅氧烷(PDMS)和聚三甲基硅烷基苯乙烯(PSI)体系进行模拟，发现

26、COMPASS力场对于计算氧气和氮气在PDMS和PSI中的扩散非常接近实验值,氧气运动范围大于氮气，小分子的运动轨迹基本与聚合物的自由体积分布相对应。Hofmann等先后采用PCFF和COMPASS等力场对小分子气体在PDMS、聚酰亚胺(PI)等聚合物中的扩散进行了模拟，并分析了聚合物的结构对扩散的影响，他们发现，由于玻璃态聚合物链段结合紧密而难于形成分子扩散的通道，小分子气体在其内部的扩散系数要小于其在橡胶态聚合物中的扩散系数。Boyd等研究了甲烷在聚对苯二甲酸乙二酯（PET)、聚2,6-萘二甲酸乙二酯(PEN)、聚乙烯(PE)和聚异丁烯(PIB)等聚合物内的扩散，分析了体系温度对扩散的影响

27、；研究发现随着温度的升高，甲烷跳跃频繁，扩散速率也增大；由于PIB的自由体积小于PE，甲烷在PIB内的扩散系数要小于其在PE中的扩散系数。在这些气体扩散的模拟研究中，模拟体系中气体分子的数目是任意给定的，这势必会影响气体扩散性质计算的准确性，因为多孔介质内吸附质的扩散和其在吸附剂内的密度紧密相关，如果任意给定小分子在聚合物微孔中的数目，经常不能反映实际条件下的吸附量而导致得到的模拟结果不可靠。2010年，Farkhondeh等26在很大温度范围内对氩、甲烷、二氧化碳、氮气等多种气体在聚苯乙烯中的渗透和扩散进行了分子动力学模拟，计算得到的扩散系数与实验数据以及以前的模拟结果吻合很好。2010年G

28、arruchet等27应用分子动力学方法对在950 K-1600 K时氧原子在镍中的扩散进行了研究，发现氧气扩散符合空位扩散机理。这些研究都对以后的科学研究提供了宝贵的经验。然而，进一步扩大分子动力学模拟的应用面仍然是一个热门的发展方向。1.3本文的研究意义能源作为重要的战略资源，关乎国家的军事安全、经济发展、乃至国民生活。在当今能源消费结构中，石油因其广泛的使用与不可再生性,成为最重要的能源之一。随着我国汽车保有量的增加和经济持续的高增长,我国能源不足的问题日益严重。尤其近年来，我国石油的对外依存度逐年攀升，这严重威胁到国家的能源安全。随着人民生活水平的提高和车价的不断降低，人们积累的汽车需

29、求不断得到释放。预计我国汽车保有量在今后一段时间内将持续增加，因此国内车用燃料油需求会进一步增大。为了减少对石油的消耗，寻找车用替代燃料势在必行。 二甲醚是一种理想的车用替代燃料。常温常压下，二甲醚是一种淡淡香味的无色无毒气体。同时二甲醚不会腐蚀金属。二甲醚在光化学反应中不会产生甲醛，因此不会破坏臭氧层，是一种环保燃料。天然气、甲醇、液化石油气的十六烷值都不超过10只能用于点燃式发动机。而二甲醚的诸多性质决定了它适用于压燃式发动机，可以做柴油发动机的替代燃料。首先，二甲醚的十六烷值大于55，自燃温度低于柴油；其次，二甲醚的蒸发压力低，容易被液化。二甲醚做为柴油发动机代用燃料的优势。二甲醚的化学

30、结构式为CH3-0-CH3，它是最简单的一种醚结构，其中不含C-C键，只有C-O键和C-H键，氧气质量分数为34.8%，因此燃烧后生成的一氧化碳 、碳氢化合物和微粒比较少，发动机可以承受较高的废气再循环率以降低氮氧化合物的生成和排放。二甲醚的低热值比柴油小，因此为了达到和柴油相当的动力水平，必须增大每循环供油量。二甲醚理论混合气热值是3066.7 KJ/Kg,而二甲醚的理论混合气热值是2911KJ/Kg，因此柴油机燃用二甲醚的升功率会升高。二甲醚的沸点低，这使得它在喷入气缸后就可以气化，它的油束雾化特性将明显优于柴油，有可能在低喷射压力下就能满足燃烧要求。二甲醚作为柴油机代用燃料还可以降低发动

31、机噪音。现有的汽车发动机只要略加改造就可以使用二甲醚燃料。二甲醚是重要的甲醇衍生物，传统上采用气相甲醇催化脱水法生产，这种方法的优点是简单，缺点是产量低，成本高。Haldor Topsoe A /S公司开发了一种由合成气生产二甲醚的低成本技术，省去了中间制甲醇步骤，而这种合成气几乎可以从任何含碳的原料转化而得到，如煤、天然气、木材等，从而大大降低了生产成本。所以，二甲醚是一种可以不依赖于石油的可再生燃料。利用我国丰富的煤炭资源，再结合二甲醚生产技术与工艺的改进和批量的增加，二甲醚的价格将会降低，可以加快二甲醚用作柴油机的代用燃料的推广。城市公交汽车由于行驶里程多，运行又集中于市区，重型大客车和

32、汽油车因汽缸容量大，排气量大，排放的氮氧化物较多，因而城市公交汽车已成为城市机动车尾气污染治理的重中之重。由于二甲醚具有清洁燃料性能，所以二甲醚燃料汽车将有望成为新一代城市客车的替代燃料。二甲醚做为车用燃料有助于缓解我国柴油供应紧张的局面，有利于环境保护，市场前景十分广阔。但DME是一种溶剂，能够溶解多种橡胶材料。传统柴油机如果不经过改造直接使用DME燃料，燃油供给系统中的耐油橡胶密封件会因长期接触DME发生溶胀和老化等现象，力学性能会逐渐变差从而失去密封作用，影响发动机的可靠运行。国内外关于DME作为柴油机代用燃料的研究时间还很短，对橡胶材料耐DME的性能研究还不足。因此研究开发耐DME的橡

33、胶密封材料具有重要的实际意义，有利于DME汽车的推广。通过研究二甲醚在橡胶中的扩散系数可以了解橡胶的耐二甲醚性能。预测扩散系数的方法可以分为理论方法、经验关联和实验方法三类.理论方法主要有Stokes-Einstein方法、Darken方法、Hartley-Crank方法、Eyring方法和各种摩擦系数方法。它们一般不能精确地预测液相的扩散系数。经验关联方法因为实验数据有限，应用范围也很窄。实验方法包括物理实验和计算机实验.由于物理实验代价昂贵，所以计算机实验方法受到广泛的重视。本文运用分子动力学方法对二甲醚在橡胶材料的扩散系数进行研究。随着计算机模拟技术的发展，运用分子动力学方法获得的数据和

34、实际试验数据比较接近。分子动力学模拟通过对研究物系中粒子的运动方程的求解从而得到粒子的运动速度和运动轨迹，通过统计计算即可求得扩散系数，由于是计算机实验，它的成本比较低，并且可以缩短试验时间，因此它是研究扩散系数等传递性质的有效手段。并且运用分子动力学模拟可以看到很多在实际实验中无法获得的微观细节，这可以在橡胶材料的研究和开发中发挥重大作用。1.4本文研究内容本文分别用Einstein方法和Green-Kubo方法计算二甲醚在氯丁橡胶、氟橡胶、丙烯腈含量为20%的丁腈橡胶和丙烯腈含量为15%、20%、30%的氢化丁腈橡胶在三种不同温度下的扩散系数，并对两种算法进行比较。第二章 分子动力学方法介

35、绍2.1模拟软件Materials Studio2.1.1 Materials Studio软件简介美国ACCELRYS公司在2000年初推出了新一代的模拟软件Materials Studio，将高质量的材料模拟带入了个人电脑(PC)的时代。本文计算采用的软件正是Accelrys公司推出的基于PC平台的材料模拟软件Materials Studio。Materials Studio是Accelrys专门为材料科学领域开发的可以在个人电脑上运行的材料计算软件，能够帮助研究人员解决现在化学及材料工业中的许多重要问题。支持Windows, Unix以及Linux等多种操作平台的Materials St

36、udio使化学及材料科学的研究者们能更方便的建立三维分子模型，深入的分析有机、无机晶体、无定形材料以及聚合物。该软件可以生成高质量的图片，同时可以处理各种不同来源的图形、表格和文本。因为Materials Studio可以综合运用多种先进算法，所以它的模拟能力非常强大。不管论是性质预测、聚合物建模还是X射线衍射模拟，我们都可以通过一些简单易学的操作来得到切实可靠的数据。该软件模拟的内容包含了催化剂、聚合物、固体化学、结晶学、晶粉衍射以及材料特性等材料科学研究领域的主要课题。 2.1.2 使用模块简介Materials Studio包含了很多模块，每个模块都有不用的功能。本文研究二甲醚在橡胶中的

37、扩散系数，主要用到其中的Amorphous Cell, Discover和Forcite三个模块。（1)Amorphous Cell 模块Amorphous Cell是一个可以建立复杂无定型系统的代表性模型，并且可以预测它的主要性质的模块。 Amorphous Cell创建结构采用的基于建立的产生包含链分子的疏松无序系统。它的特点是：任意混合体系的建模方法、特殊的产生有序向列型中间相以及层状无定型材料的能力。它通过观察体系结构和性质的关系，对分子的一些重要性质有更深入的了解，从而设计出更好的新物质的配方。它可以预测并研究的性质包括状态方程行为、内聚能密度、链堆砌和局部链运动。(2)Discov

38、er模块Discover是Materials Studio中非常重要的模块，它提供了功能强大的原子模拟方法，适用于大部分的分子和材料。Discover集成的很大范围已经被证明适用于分子设计的动力学方法，它以PCFF、COMPASS、CVFF等多个力场为基础，能够准确地计算出分子的最低能量构像，同时能够给出不同系综下体系结构的动力学轨迹。Discover为Amorphous Cell等模块提供了计算的基础。周期性边界条件的引入使得它可以对固态体系进行研究，如晶体、非晶和溶剂化体系。Discover还提供强大的分析工具，能对模拟结果进行分析，从而得到各类结构参数、力学性质、动力学量和热力学性质等。

39、(3)Forcite模块Forcite是分子力学模块，能够使用经典力学对任意分子和周期系统进行势能和几何优化计算。Forcite支持UFF,、COMPASS和Dreiding力场。因为支持力场的广泛性，Forcite原则上可以处理任何材料。几何优化算法提供了最速下降法、共扼梯度法和牛顿法，并且提供了连续使用这些方法的灵巧方法。这使得研究人员可以准确地进行能量最小化计算。2.2分子动力学基本原理2.2.1牛顿运动方程和数值积分算法分子动力学是对物理系统确定的微观描述，该系统既可以是多体系统，也可以是少体系统。该方法的基本原理是:建立一个粒子系统，对所研究的微观现象进行仿真，各粒子间的相互作用根据

40、量子力学来确定。对于符合经典牛顿力学规律的大量粒子系统，通过粒子运动学方程组的数值求解，得到粒子在相空间的运动规律和轨迹，然后按照统计物理原理得出该系统相应的宏观物理特性。以下是分子动力学模拟的基本原理: 在包括N个粒子的分子动力学模拟系统中，系统的总能量是系统中粒子动能与势能的总和。其中，势能是系统中粒子位置的函数： 。 通常认为势能为各粒子对间势能的加成，即: （2-1） 公式中，表示第i、 j个粒子的势能；是粒子i、 j间的距离。 按照经典力学，系统中任何一个粒子所受的力为势能的梯度，则粒子j对粒子i的作用力为： （2-2） 公式中是粒子j对粒子i的作用力；是粒子i对粒子j的作用力。根据

41、叠加原理和牛顿运动定律,可得到第i个粒子所受的作用力和加速度的表达式为: （2-3） （2-4）根据加速度与粒子的位移和速度的导数关系，可以对加速度积分从而得到第i个粒子的速度和位移: （2-5）公式中， 分别是第i个粒子的位置和速度。然后，一般采用各种有限差分法来求解运动方程，从而得到粒子的运动，实现对粒子的运动规律和运动轨迹的模拟，最终实现对系统的宏观特性模拟计算的目的。有很多可以求解牛顿运动方程的有限差分算，但是分子动力学计算的原子数目比较多，因此对空间和时间的要求都比较高，需要给出适合用于分子动力学的时间积分算法。时间积分算法是基于有限差分法，对时间进行离散，是时间步长，根据t时刻的变

42、量和它的时间导数，可以推演得到t+时刻的物理量，逐步计算就可以得到整个历史过程的物理量。这个过程是近似的，因此必然会有误差，包括截断误差和舍入误差。截断误差主要是因为有限差分法用Taylor展开截断某些高阶项引入的，舍入误差是由于计算机的本身数值精度引入的。两种误差都可以通过减小时间步长而减小。对于比较大的，主要是截断误差起作用，但是随着的减小截断误差会迅速降低。舍入误差随着下降相对来说比较慢，在较小的时起主导作用，采用双精度数据类型（64位）有助于把截断误差降到最低。在分子动力学计算领域主要采用Gear算法和Verlet算法。(1)Verlet算法在1967年被提出，之后有对它进行了一些变化

43、。该算法基于Taylor展开： （2-6）两项相加，可以得到 （2-7） 这就是Verlet算法的基本形式。从此式能得到粒子在相空间的轨迹，而且该式计算的截断误差是。公式（2-7）中没有v(t),可以通过下面的公式计算： （2-8）得出粒子的速度就可以计算系统的动能，从公式（2-8）可以看出速度计算的截断误差是O。在下一步的位置没有得到之前，难以得到速度项。而且，它不是一个自启动算法。在t=0时，只有一组位置，因此需要运用别的方法得到的位置。为了解决这个问题，对Verlet算法进行了一些改变。Verlet算法的速度形式具有的形式如下： （2-9）这种算法的缺点是需要储存的量较多，但这种算法可以

44、同时给出位置、速度和加速度，本文选用这种算法。(2)Gear算法在分子动力学模拟中，在微正则系综中保持能量守恒的前提下，要尽可能用长的时间步长，或在时间步长相同时得到较高的精度，保存前一步的力，并且用预测-校正算法更新位置和速度。这个方法计算过程如下:首先，根据Taylor展开式，预测新的位置、速度和加速度。然后，根据新位置计算得到的力，求出加速度。这个加速度再和Taylor级数展开式预测的加速度进行比较。两个加速度的差在校正步里用来校正位置和速度项。 2.2.2分子动力学模拟的基本步骤MD模拟可以分为以下几个步骤：（1） 确定起始的构型,把由N个粒子构成的系统看成是N个相互作用的质点，这些质

45、点有质量、初始温度、坐标等。然后选择适当的势函数。（2）给定初始条件,对系统中所有粒子给以初始位置和速度，根据目标温度按照波尔兹曼分布随机指定各质点的初速度。（3）趋于平衡计算,我们确定了模拟需要的初始条件和边界条件后就求解运动方程和模拟，但是我们需要系统具有的能量在刚开始时系统并不具备，为了使系统达到平衡态，我们在模拟过程中必须设立一个趋衡过程。在这个过程中我们从系统中移出或者增加能量，一直到系统持续给出确定的能量值。（4）宏观物理量的计算,宏观物理量的计算一般在模拟的最后阶段进行，它是沿着相空间轨迹通过求平均计算得到。2.3分子动力学模拟的关键技术2.3.1系综研究人员对同一个结构研究的目

46、的不同，在进行分子动力学模拟时会使用不同的系综。经常使用的平衡系综有微正则系综(NVE)、正则系综(NVT)、等温等压系综（NPT)和等焓等压系综（NPH)。微正则系综(NVE) ，它是一个孤立的、保守的系统的统计系综。在这种系综中，系统沿着相空间中的恒定能量轨道演化，在分子动力学模拟的过程中，系统中的原子数N，体积V和能量E都保持不变。正则系综(NVT)，在这个系综中，系统的原子数N，体积V和温度T都保持不变，并且总动量为零。一般保持系统的温度不变的方法是让系统与虚拟的热浴处于热平衡状态，由于温度与系统的动能有直接的关系，通常的做法是把系统的动能固定，这可以通过对原子的速度进行标度来实现，也

47、可以对运动方程中出现的力加一个约束力项来实现。本文采用NVT系综。等温等压系综（NPT)，即系统的原子数N，压力P和温度T都保持不变。温度的恒定以前一样，是通过调节系统的速度或加一约束力来实现的。因为系统的压力P与其体积V是共轭量，所以要调节压力值可以通过标度系统的体积来实现，现在有很多调压的方法都是采用的这个原理。等焓等压系综（NPH)，就是保持系统的原子数，压力和焓值都不变，所以在该系综下进行模拟时要保持压力与焓值是固定值，其调节技术的实现的难度较大。2.3.2初始条件对一个系统进行模拟，需要对体系中所有粒子赋予初始位置和速度，选择的粒子位置应该和所模拟的对象相适应。但是，通常系统的初始位

48、置都是未知的。从实际来看，精确选择要研究的系统的初始条件是无意义的，因为模拟的时间足够长时，系统就会忘掉初始条件。但初始条件会对模拟过程产生一些影响，初始条件的不同系统趋于平衡的时间也不同。常用的初始条件有:让初始位置在差分划分网格的格子上，初始速度是从玻尔兹曼分布随机抽样得到；让初始位置随机地偏离差分划分网格的格子，同时初始速度为零；让初始位置随机地偏离差分划分网格的格子，初始速度同样从玻尔兹曼分布随机抽样得到。 2.3.3势函数在分子动力学模拟中，相互作用势是用分子动力学模拟材料性质的基础和关键。原子间相互作用势的选择是决定模拟结果在多大程度上符合实际的关键因素，同时，模拟结果也可以用来检

49、验原子间相互作用势。从分子动力学的发展进程来看，原子之间的相互作用势的研究发展的很慢，这影响了分子动力学在实际研究中的应用。原子间势函数概念已经把电子云对势函数的贡献考虑在里面了。原子间势函数刚开始在研究对势，后来发展到多体势。对势理论认为，原子间的相互作用是两两之间的作用，和其他原子的位置无关。换句话说就是任何两个原子结合键的强度不因为它周围其它原子的存在而受影响。我们可以假设在N个原子中作用能量是相互独立的二体作用之和。它的势能泛函式为: （2-10） 对势能够比较好地描述除金属和半导体以外的几乎所有无机化合物。有的对势是经过一定的理论分析而得到的，但其中一些参数则需要根据宏观试验参数用经

50、验方法来确定，这些宏观实验参数主要有弹性常数、平衡点阵常数以及内聚能、空位形成能和层错能等，这些称为半经验势。后来，人们经过实践又发展了经验势。所谓经验势就是人们为了便于拟合，在选择势函数的形式时，并不一定要求有确切的理论依据，而是出于经验的估计和便于拟合，相对自由的选择势函数形式。下面对常用的几种势函数进行介绍： Lennard-Jones(L-J)势：Lennard-Jones(L-J)势函数的解析表达式如下所示: （2-11）其中，U。表示相互作用的强度，反映了原子的大小。根据量子力学二次微扰论的偶极子一偶极子相互作用可导出=12,这一项描述了范德瓦耳斯力。后面一项是排斥力，其来源之一是

51、原子核之间的库仑斥力，来源之二是电子之间由于泡利不相容原理产生的交叠能。在实际使用中Lennard-Jones势常取: （2-12）上式中，为长度参数，为能量参数。上式是Lennard-Jones势最为常用的形式，但也经常采用别的形式的排斥项，例如指数形式。用指数排斥项和倒六次方吸引项构成的势，通常称为Buckingham势。Morse势：Morse注意到双原子分子的振动谱的量子力学问题可用指数形式的势函数解析地解决，并发现计算结果与实验一致。因此他提出如下形式的势函数： （2-13）Morse势和Lennard-Jones势的曲线形式很接近。Morse势经常被用来构造各种多体势的对势部分。B

52、orn-Mayer势：Born和Mayer认为碱金属离子之间的排斥项可用指数形式表示，因此提出如下形式的势函数： （2-14）Born-Mayer排斥项等价于Morse势的第一项。其中第三项表示偶极子一偶极子相互作用，第四项表示偶极子一四极子相互作用。反映了泡利不相容效应。该势考虑了离子之间的库仑力，所以它能很好地应用于碱金属卤化物和碱土金属卤化物等各种离子型化合物。 2.3.4常见力场目前用于材料研究领域的分子力场主要有CVFF、UFF、PCFF和COMPASS等。CVFF(Consistent Valence Force Field)是一种典型的较早应用的分子力场算法，这个算法包含了部分非

53、简谐作用以及能量藕合作用项，但并不是全部的。适用于高分子小分子、某些晶体和气相物质的计算。它可以处理氨基酸、蛋白质和大量高分子系统。它主要用于研究物质的结构和结合能，并能有效预测振动频率和构象能。UFF(Universal Force Field)力场适用的范围最广泛，它可以模拟的元素几乎覆盖了整个元素周期表。能够较为准确的预测高分子、主族无机分子和金属的结构和构型能。可以用于计算高分子金属系统和其他力场没有参数化的系统。PCFF(Polymer Consistent Force Field)力场是在CFF91力场基础上发展起来的，侧重于对高分子和高分子材料的研究，诸如聚碳酸酯、多聚糖、碳水化

54、合物、脂类、核酸等。COMPASS (Condensed-phase Optimized Molecular Potential for Atomistic Simulation Study)力场是第一个由凝聚态性质以及孤立分子的各种从头算和经验数据经参数化并得到广泛验证的从头算力场，也是第一个将原来分别处理有机体系和无机分子体系的力场统一起来了的力场。COMPASS力场适用的体系广泛，包括常用的高分子和无机分子、一些金属离子、金属氧化物和金属，以及用于研究表面、共混等复杂体系都可以使用这个力场。经过大量的试验表明这个力场能够准确的预报孤立态和凝聚态的结构、构象、振动以及内聚能等热力学性质，可

55、以说COMPASS力场是当前发展得相当完善的“求全”型的分子力场。因为使用COMPASS力场可在很大温度、压力范围内精确地预测出孤立体系或凝聚态体系中各种分子的结构、构象、振动以及热物理性质，因此本文的力场就选择为COMPASS力场。2.3.5周期性边界条件与最近镜像理论在分子动力学模拟中，为了保持系统的密度不变，通常采用周期性边界条件。就是认为在所模拟系统周围排布着和模拟系统内部原子排布和运动相同的镜像系统，并且当有一个原子进入模拟盒子时，模拟系统中一定有一个原子从相反的方向移出模拟盒子。为了进一步说明周期性边界条件，下面给出周期性边界条件的二维图像 如图2-1所示。当系统中的粒子按照周期性

56、边界条件排布时，计算粒子间的作用力就必须遵守最小镜像原理。最小镜像原理如图2-2所示 。就是计算粒子1与粒子3的相互作用力时，只计算粒子1与离其最近的镜像3间的相互作用力。因为在计算中应用最小镜像原理，在实际的模拟计算中必须引入截断半径，只计算截断半径以内的粒子间的相互作用力，一般取，这样就能够有效地避免重复计算分子间的作用力。 图2-1 二维周期性边界条件 图2-2粒子的镜像原理 2.3.6截断距离与电荷基团（1）非键截断距离势函数中非键相互作用的直接计算包括所有的原子对，计算量的大小和系统的大小密切相关，当系统增大时计算量会迅速增大。为了缩短模拟的时间，在模拟比较大的系统时，为了减少计算量

57、，必须采用特殊的算法。在处理高分子系统时，引入非键截断距离是最常用也是最重要的方法。使用该方法时，超过非键截断距离的原子对之间的相互作用会被忽略。在计算的时候，必须慎重的选择截断距离的大小来平衡计算的效率和精度。在计算过程中，如果仅仅采用截断距离会导致势能面不是连续的，进而体系的模拟结果会变的不稳定。为了避免由于截断引起的势函数的不连续性，需要运用一定的处理方法。通常的做法是引入switching函数在一定距离（spline宽度)内能够让能量的变化变的缓慢，同时计算量又不会过大，其作用原理如 2-3图所示 （2-15）公式内的是加权因子(01),的函数形式如下方程所示： （2-16 ) 图2-

58、3 非键阶段距离和switching函数(2)电荷基团通常采用基于基团的截断距离来减小直接采用基于原子截断距离而引起的人为的单极一单极相互作用和偶极分裂。电荷基团是由临近原子组成中性基团。一般将电荷基团中临近几何中心的原子定义为中心原子。如果中心原子在截断距离之内，电荷基团中所有的原子都会包括在相互作用的计算中。电荷基团通常是化学官能团。计算时，引入电荷基团能够减小非键相互作用的计算量，可以提高计算的速度同时计算的精度不会受到影响。2.3.7长程力的计算处理短程力一般用周期性边界条件和截断半径，但是对于静电引力等长程有效作用力，简单的位能截断已经不再适用。需要使用别的处理方法，经常使用的处理方

59、法如下： (1) Tree-Based方法基于树结构的方法是一种被广泛应用的计算算法，这种方法是基于树的方法计算效率非常高，在处理长程作用方面的优势非常明显，但是不经常被用在计算模拟上。 (2)Ewald Sum加和法Ewald Sum加和是比较常用于处理长程力作用的方法，它的作用原理是在空间里处理截断半径内的作用，周期性边界条件下处理截断半径以外的长程力。 (3)反应场方法反应场方法也是一种常见的计算长程力的方法，它基本思想是作用偶极子的场分成两个部分。一部分是由位于截断半径内的短程作用，另一部分由截断半径外的作用贡献，这部分堪称是介电常数的反应场。 (4)PPMM法PPMM法由Eastwo

60、od等人在1980年提出，短程力处理PP(particle-particle)法，与长程力计算处理方法比较相似;长程力处理采用PM(particle-meth)法。PP和PM两种方法结合得到的结果比较精准。2.3.8调温控制因为外部或者摩擦力引起的过热，积分误差和力的截断，这些情况都会引起能量的飘移，所以控制系统的温度就显得非常必要。经常使用的控温方法有: 速度标度法、Gaussian热浴、Nose-Hoover方法和Berendsen弱耦合方法。 （1）速度标度法改变温度最简单和最直观的方法是直接对速度进行标度，如果t时刻的温度是，速度乘以因子后温度的变化是 （2-17） （2-18）控制温度最简单的方法是在每一步乘上标度因子。上面公式中，是期望的参考温度。（2） Gaussian热浴这种方法的思路是在运动方程中加上一项“摩擦力”项，这个力又和系统的速度联系在一起。运动方程如下 (2-19)上式中，用来保持总动能守恒，为了把温度调到期望的温度下，需要把初始的温度设定为期望值。(3)N