分子动力学简介及其应用摘要:综述了分子动力学的发展历史,基本步骤,作用势与动力学计算,时间步 长与约束动力学,及其应用,最后还指出了分子动力学的进一步研究方向 关键字:分子动力学,作用势,时间步长,约束动力学1. 引言分子动力学是一套分子模拟方法,该方法主要是依靠牛顿力学来模拟分子体 系的运动,以在由分子体系的不同状态构成的系综中抽取样本,从而计算体系的 构型积分,并以构型积分的结果为基础进一步计算体系的热力学量和其他宏观性 质1957年,Alder和Wainwrights首先在硬球模型下,采用分子动力学研究气 体和液体的状态方程,从而开创了利用分子动力学模拟方法研究物质宏观性质的 先例后来,人们对这一方法作了许多改进,并运用它对固体及其缺陷以及液体 作了大量的研究但由于受计算机速度及内存的限制,早期模拟的空间尺度和时 间尺度都受到很大限制21 世纪 80年代后期,由于计算机技术的飞速发展,加 上多体势函数的提出与发展,为分子动力学模拟技术注入了新的活力分子动力 学模拟不仅能得到原子的运动细节,还能像做实验一样进行各种观察对于平衡 系统,可以用分子动力学模拟作适当的时间平均来计算一个物理量的统计平均 值。
对于非平衡系统,发生在一个分子动力学观察时间内(一般为1〜100ps)的物 理现象也可以用分子动力学计算进行直接模拟特别是许多在实际实验中无法获 得的微观细节,而在分子动力学模拟中都可以方便地观察到这种优点使分子动 力学存物理、化学、材料科学等领域研究中显得非常有吸引力[2〜10]分子动力学模拟虽然不如第一原理模拟精确,但以程序简单,计算量小,可 计算的原子体系大大超过第一原理等方法,而保持有巨大的发展和应用前景2. 分子动力学简史1957年:基于刚球势的分子动力学法(Alder and Wainwright)1964年:利用Lennard-Jones势函数法对液态氩性质的模拟(Rahman)1971年:模拟具有分子团簇行为的水的性质(Rahman and Stillinger)1977 年:约束动力学法(Rychaert,Ciccotti)1980年:恒压条件下的动力学法(Ander法,Parrinello-Rahman法)1983年:非平衡态动力学法(Gillan and Dixon)1984年:恒温条件下的动力学法(Berendsen et al.)1985年:第一原理分子动力学法(Car-Parrinello法)1991年:巨正则系综的分子动力学方法(Cagin and Pettit)3. 基本步骤3.1 确定起始构型进行分子动力学模拟的第一步是确定起始构型,一个能量较低的起始构型是进行分子模拟的基础,一般分子的起始构型主要来自实验数据或量子化学计算。
在确定起始构型之后要赋予构成分子的各个原子速度,这一速度是根据波尔 兹曼分布随机生成的,由于速度的分布符合波尔兹曼统计,因此在这个阶段,体 系的温度是恒定的另外,在随机生成各个原子的运动速度之后须进行调整,使 得体系总体在各个方向上的动量之和为零,即保证体系没有平动位移3.2 进入平衡相由上一步确定的分子组建平衡相,在构建平衡相的时候会对构型、温度等参 数加以监控3.3 进入生产相进入生产相之后体系中的分子和分子中的原子开始根据初始速度运动,可以 想象其间会发生吸引、排斥乃至碰撞,这时就根据牛顿力学和预先给定的粒子间 相互作用势来对各个粒子的运动轨迹进行计算,在这个过程中,体系总能量不变, 但分子内部势能和动能不断相互转化,从而体系的温度也不断变化,在整个过程 中,体系会遍历势能面上的各个点,计算的样本正是在这个过程中抽取的3.4 计算结果用抽样所得体系的各个状态计算当时体系的势能,进而计算构型积分4. 作用势与动力学计算作用势的选择与动力学计算的关系极为密切,选择不同的作用势,体系的势 能面会有不同的形状,动力学计算所得的分子运动 和 分子内部运动的轨迹也会 不同,进而影响到抽样的结果和抽样结果的势能计算,最初的分子动力学计算采 用比较简单的刚球势,现在更多地采用兰纳-琼斯势,后者能够更好的与粒子间 相互作用拟合。
5. 时间步长与约束动力学分子动力学计算的基本思想是赋予分子体系初始运动状态之后利用分子的 自然运动在相空间中抽取样本进行统计计算,时间步长就是抽样的间隔,因而时 间步长的选取对动力学模拟非常重要太长的时间步长会造成分子间的激烈碰 撞,体系数据溢出;太短的时间步长会降低模拟过程搜索相空间的能力,因此一 般选取的时间步长为体系各个自由度中最短运动周期的十分之一但是通常情况下,体系各自由度中运动周期最短的是各个化学键的振动,而 这种运动对计算某些宏观性质并不产生影响,因此就产生了屏蔽分子内部振动或 其他无关运动的约束动力学,约束动力学可以有效地增长分子动力学模拟的时间 步长,提高搜索相空间的能力6. 考虑事项1)软件的选择,这通常和软件主流使用的力场有关,而软件本身就具体一定的偏向性,比如说,做蛋白体系,Gromacs, Amber, Namd均可;做DNA, RNA体系,首选肯定是Amber;做界面体系,D1_POLY比较强大,另外做材料体系,Lammps 会是一个不错的选择2)力场的选择力场是来描述体系中最小单元间的相互作用的,是用量化 等方法计算拟合后生成的经验式常见的有三类力场:全原子力场,联合力场, 粗粒化力场。
3)通过实验数据或者是某些工具得到体系内的每一个分子的初始结构坐标 文件,之后,我们需要按我们的想法把这些分子按照一定的规则或是随机的排列 在一起,从而得到整个系统的初始结构,这也是我们模拟的输入文件4)结构输入文件得到了,我们还需要力场参数输入文件,也就是针对我们 系统的力场文件,这通常由所选用的力场决定,比如键参数和非键参数等势能函 数的输入参数5) 体系的大小通常由你所选用的box大小决定,我们必须对可行性与合理 性做出评估,从而确定体系的大小,这依赖于具体的体系6) 由于初始构象可能会存在两个原子挨的太近的情况(称之为 bad contact), 所以需要在正式模拟开始的第一步进行体系能量最小化,比较常用的能量最小化 有两种,最速下降法和共轭梯度法,最速下降法是快速移除体系内应力的好方法, 但是接近能量极小点时收敛比较慢,而共轭梯度法在能量极小点附近收敛相对效 率高一些,所有我们一般做能量最小化都是在最速下降法优化完之后再用共轭梯 度法优化,这样做能有效的保证后续模拟的进行7) 以平衡态模拟为例,你需要设置适当的模拟参数,并且保证这些参数设 置和力场的产生相一致,举个简单的例子, gromos 力场是用的范德华势双截断 来定范德华参数的,若你也用 gromos 力场的话也应该用双截断来处理范德华相 互作用。
常见的模拟思路是,先在NVT下约束住你的溶质(剂)做限制性模拟, 这是一个升温的过程,当温度达到你的设定后,接着做 NPT 模拟,此过程将调 整体系的压强进而使体系密度收敛经过一段时间的平衡模拟,在确定系统弛豫已经完全消除之后,就可以开始 取数据了如何判断体系达到平衡,这个问题是比较技术性的问题,简单的讲可 以通过以下几种方式,一,看能量(势能,动能和总能)是否收敛;二,看系统 的压强,密度等等是否收敛;三看系统的 RMSD 是否达到你能接受的范围,等 等8)运行足够长时间的模拟以确定我们所感兴趣的现象或是性质能够被观测到,并且务必确保此现象出现的可重复性9)数据拿到手后,很容易通过一些可视化软件得到轨迹动画,但这并不能 拿来发文章真正的工作才刚刚开始——分析数据,你所感兴趣的现象或性质只 是表面,隐含在它们之中的机理才是文章中的主题7. 应用分子动力学的计算过程给定了体系的总能量,因此适用于对微正则系综的模 拟计算,另外由于分子动力学计算过程始终是时间的函数,因此一些与时间有关 的宏观量如扩散系数的模拟必须应用分子动力学另外,在实际应用中,经常把分子动力学方法和蒙特•卡罗方法联合使用。
在近年,多尺度模拟计算已经得到了很多学者的关注在多尺度模拟计算中,分 子动力学方法研究纳米级现象,蒙特•卡罗方法研究微观形态,有限元方法应用 于宏观领域通过多种尺度的多种模拟计算方法的联合应用,令纳观与宏观联结 起来8. 进一步研究方向对于分子动力学模拟而言,最重要的两要素是:初始位形的给定和原子间作 用势的确定系统的恰当的初始位形可必根据实验而得到或者根据物理分析人为 地预先给出,以大大减少计算量影响结果精确程度的最主要因素是原子间作用 势的精确性人们为了提高势函数的精确性做了许多研究,提出了许多势函数的 形式,但是对大多数原子而言,努力找到尽可能精确而形式又不太复杂的势函数, 仍然是极大的挑战把第一原理计算,量子化学分析与参数拟合相结合是势函数 研究的最佳方法值得指出,半经典的 Thomas-Fermi-Dirac Theory 及由此导出 的位力定理[11等] 定律在分析原子同作用势的可能形式时也是有用的对模拟得到的位形进行结构分析的方法虽然已经提出了上述数种,然而不能 认为它们已经构成了凝聚态结构分析方法的完备集合寻找更多的结构分析方 法,以便对模拟得到的位形进行更详尽的结构分析也是分子动力学模拟领域的一 个重要的发展方向。
参考文献[1] Alder B J.Wainwright T E.Phase transition for a hardsphere system.The Journal of Chemical Physics, 1957, 27: 1208〜1209[2] Schiotz J.Tolla F D D, Jacobsen W.Softening of nanocrystalline metals at very small grain size.Nature, 1998, 391: 561〜563[3] Swygenhoven H V, Caro A.Molecular dynamics computer simulation of nanophase Ni: structure and mechanical properties. Nanostructured Materials, 1997,9: 669〜672[4] Phillpot S R, Wolf D, Gleiter H Molecular-dynamics study of theSynthesis and characterization of a fully dense , three-dimensionally nanocrystalline materials.Journal Applied Physics, 1995, 78(2): 847〜861[5] Faux D A.Molecular dynamics studies of sodium diffusion in hydrated Na-Zeolite-4A.Journal of Physics Chemistry B,1998, 102: 10658〜10662[6] Zeng P,Zajac S, Clapp P C, Rifkin J A.Nanoparticle sintering simulations Materials Science and Engineerin A, 1998,252: 301〜306[7] Chen z Y, Ding J Q.Molecular dynamics studies on dislocation in crystallites of nanocrystalline iron Nanostructur。