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

分子动力学简介及其应用

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

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

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模拟，此过程将调整体系的压强进而使体系密度收敛。