分子动力学简介及其在热传导方面的应用

分子动力学简介及其应用

摘要：综述了分子动力学的发展历史，基本步骤，作用势与动力学计算，时间步长与约束动力学，及其应用，最后还指出了分子动力学的进一步研究方向。

关键字：分子动力学，作用势，时间步长，约束动力学

1.引言

分子动力学是一套分子模拟方法，该方法主要是依靠牛顿力学来模拟分子体系的运动，以在由分子体系的不同状态构成的系综中抽取样本，从而计算体系的构型积分，并以构型积分的结果为基础进一步计算体系的热力学量和其他宏观性质。

1957年，Alder和Wainwright[1]首先在硬球模型下，采用分子动力学研究气体和液体的状态方程，从而开创了利用分子动力学模拟方法研究物质宏观性质的先例。后来，人们对这一方法作了许多改进，并运用它对固体及其缺陷以及液体作了大量的研究。但由于受计算机速度及内存的限制，早期模拟的空间尺度和时间尺度都受到很大限制。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；做界面体系，Dl_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 the

Synthesis 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 Nanostructured Materials，1998，10(2)：205～215

[8] Wen Yuhua，Zhou puxin．Liu Yuewu．Molecular dynamics simulations of Microstructure of nanocrystalline Copper．Chinese Physics Letters，2001，31(3)：411—413

[9] Wen Yuhua，Zhou Fuxin，Liu Ynewu. The influence of grain size and temperature on the mechanical deformation of nanocrystalline materials：molecular dynamics simulation Chinese Physics，2001，10(5)：407～412

[10] 孙祉伟经典流体的计算机模拟试验——蒙特卡洛法和分子动力学法力学与实践，1983，5（6)：59～62

[11] Zhu Ruzeng，Wen Yuhun and Qian Jin．The virial theorem in refined Thomas-Fermi-Dirac theory for the inner of the atoms in a solid Chinese Physics，2002，1l(11)：1193～1195