谈谈分子模拟中的能量最小化，弛豫和平衡态

目前分子模拟在很多领域得到了广泛地应用。这主要得益于分子模拟理论的逐渐成熟和模拟工具的开发及完善。Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) 是目前材料领域应用广泛的分子模拟工具包；而Gromacs则是生物体系模拟的一把利器。根据统计热力学的要求，只有对平衡后结构做分析才能得到合理的结果。[1] 平衡态分子模拟的一般运行流程是能量最小化，弛豫、平衡模拟及采样。下面我们将详细介绍这三个流程。

1        能量最小化

分子模拟的初始构型往往偏离平衡态较远。如果初始结构中有些原子靠的比较近，这样系统的能量就会过高而使得模拟崩溃。另外如果原子靠的很近，相互作用力就是非常大，可能会导致原子在一个时间步内移出盒子边界从而出现原子丢失的现象。一般分子动力学软件包都提供了能量最小化(minimize)工具，能够在一定程度上优化结构以减少原子间的坏的接触。除了能量最小化工具，LAMMPS还提供了一些独特的手段处理比较差的初始构形：

①  pair_style soft :使用一个余弦对势将重叠的原子排开，这个势能确保原子重叠时，体系的能量不会爆炸。进行一段时间的模拟后，可能会得到比较合理的构型，从而可以用正常的力场参数继续进行模拟。

②  fix nve/limit 0.1 : 限制一个原子在一个时间步内的最大位移，从而避免原子丢失。

当然能量最小化工具只能在一定程度上缓解(alleviate)坏的构型带来的影响。为了顺利地模拟，在搭建构型时就应注意避开原子靠近的情况。周期性带来的原子重叠问题往往容易被忽略，因此需格外小心。

2        弛豫

弛豫(relaxation)指的是将非平衡的系统慢慢演变为平衡态。[2] 模拟刚开始时，系统的温度，压力等性质可能偏离设定值很远，因此需要经过一段时间的模拟趋于平衡。通常针对体系的不同，可以分几步弛豫不同的物理量。但是最终都需要一个预平衡过程中将系统弛豫到平衡。弛豫通常采用什么系综和控温、控压方法呢？

1. 1       系综的选择

①   如果需要在NPT 系综的平衡态下抽样，可以考虑先将体系在NVT下模拟，把温度T弛豫到稳态，再转NPT系综进行预平衡，从而把压强P弛豫到稳态。因为温度相较于压强容易控制，在NVT系综下模拟能够在短时间内能进入平衡态。

②   如果需要在NVT 系综的平衡态下抽样，可以先在NPT系综下将系统密度弛豫到稳态，然后再转NVT进行预平衡。

③   特殊体系处理： (A) 两相体系[1,3,4]：在构建初始构型前，可以对两相分别进行平衡，再组合成新的结构。在NVT系综下，限制固相的位置，让液体弛豫以适应固相的位置。然后放开限制，转到NPT系综下进行预平衡。 (B) 不同分子的尺寸相差较大，如蛋白质在溶液中[1,5]：在NVT系综下，限制蛋白的位置，让溶液弛豫以适应蛋白的位置。然后放开限制，转到NPT系综下进行预平衡。(C)对于平衡时间长的体系，如离子液体[1]：此类体系所需的平衡时间较长，可以考虑在NVT下高温跑一段时间，使体系空间分布变得均匀，然后退火到目标温度并弛豫温度，再转到NPT系综下进行预平衡。

1. 2       控温方法的选择

通常的速度控制方法有

①  速度标度(Velocity scaling):简单速度标度以及Berendsen 热浴

速度标度方法简单，但是调节的速度不能严格符合玻尔兹曼(Boltzmann)分布，也就是不属于严格的正则系综。Berendsen是一种弱耦合热浴，在系统远离平衡态时，对温度的调节较好，相较于Nosé-Hoover温度的震荡较小。因此可以先用Berendsen将温度弛豫到设定值，再转用Nosé-Hoover进行预平衡。

②  随机力或随机速度热浴[6]：郎之万(Langevin) 热浴， Andersen 热浴

速度中等，精确度中等。Andersen和郎之万热浴都引入了非物理的随机力，干扰了系统正常的时间演化，粒子之间的速度相关性减弱，因此算出的动力学性质是不可靠的，但是可以用来计算热力学性质。

③  恒温扩展法：Nosé-Hoover热浴[6,7]

比较复杂，计算的速度相对慢了些。严格遵守正则系综，体系可以时间反演，通常用于平衡采样。但是在系统远离平衡态时，可能温度涨落巨大，难以收敛，因此可以考虑先使用Berendsen、Langevin热浴先将系统温度弛豫到合理的值，再转Nosé-Hoover热浴进行预平衡、平衡采样。[8]

LAMMPS 中实现Berendsen和Langevin的命令分别为：fix temp/berendsen 和 fix langevin。但是这两种方式只会调节速度或力，并不会积分位置，因此需要配合 fix NVE 更新位置，因此其实系统是NVT正则系综。而实现Nosé-Hoover控温只需要一个fix nvt/npt就可以更新速度、力以及位置。

1. 3       控压方法的选择

控压方法的选择原理与控温方法基本一致。Berendsen 控压适合用于最初的压力弛豫，不适用于平衡采样。而基于Nosé-Hoover算法的Nosé-Hoover以及Parrinello-Rahman 恒压器(barostat)适用于平衡态控压[9]。前者常用于LAMMPS，而后者是Gromacs中常用的控温器。

对于各向异性的体系，比如双脂膜，均匀压力缩放(isotropic)并不适合。双脂层在xy平面无限大，而在z方向尺寸有限。因此系统在x，y方向上的压力调节与z方向的调节应该脱钩。在Gromacs中通过设置pcoupltype=semiisotropic实现[9]，而在LAMMPS中可使用fix npt 命令中分别指定x，y，z方向的参数并指定couple xy 实现。[10]

Nosé-Hoover控温和控压算法并不是万能的。Nosé-Hoover算法中温度和压力的标度并不是瞬时的，而是在一个关联时间tau的周期内弛豫温度或者压力。LAMMPS推荐的控温关联时间Tdamp为100个时间步(100.0*dt)，而控压关联时间Pdamp推荐值为1000个时间步。[10]由于关联时间的存在，在计算当前压力时会有几个步骤用到前几个时间步的压力。因此对于需要获得准确压力波动，或者系统的温度、压力关联时间比较短的情况，Nosé-Hoover算法有可能是不合适的。[9]。

3        平衡及采样

对于平衡态模拟，采样必须在确认系统进入平衡态之后才能进行。预平衡和平衡态模拟的参数(系综，控温控压方式)基本一致，但也可以通过减少预平衡过程数据输出(轨迹，日志等)以加快平衡的速度。平衡并不是意味着数值稳定。温度、压力波动的幅度都非常大，尤其是压力变化可能会出现符号的变化。压力为正说明系统有膨胀的趋势，而压力为负系统有缩小的趋势。因此判断系统是不是处于平衡态需要考察所关注的物理量的平均值是不是收敛在设定值附近。通常可以通过考察系统的势能(potential  energy)、温度(T)、压力、均方根波动(RMSD) [11] 等物理量是不是都达到了收敛从而判定系统是否进入了平衡态。再次强调，数据需要从到达平衡态的时间点开始收集。

采样的一个原则是样本数目足够多[12]。样本量越大，结果的可信度就越高。另一个原则是样本抽样要尽量无关联。合理地选取收集样本的间隔能够降低样本的关联性。[12] 一个经验做法是通过计算速度自相关函数确定大致的弛豫时间(tau)[13]，采样的时间尺度至少跨越100*tau，或者采集至少5000个样本。两相体系以及非平衡态模拟的采样请参考教程[3,12]。

平衡态模拟所采用的系综往往是按照自己的模拟需求选用的。一般可以考虑选用NPT系综进行平衡模拟，这样系综的密度能够弛豫到与实验接近的状态。而扩散系数的模拟一般是在NVT系综下，因为NPT系综可能会产生原子坐标标度的问题。如果算一些物性，比如热容，压缩系数之类的，需要在公式对应系综下抽样。[1,14,15]

定容比热容Cv定义为恒容条件下能量对温度的偏导数，这就决定了必须在定容的系综下抽样。根据热力学公式可以分别推导得到Cv在不同定容系综下的公式[14]：

其中< >代表的是系综平均。

可以看到NVT是计算定容比热容Cv最简单的系综，因为温度得到了控制，因此只需要得到势能的涨落即可计算得到Cv。

计算定压热容Cp可以在NPT系综下进行，需要得到H+PV的涨落。

参考资料：

[1]公众号“模拟小窝”：平衡态MD的经典套路

[2]公众号“模拟小窝”：什么是MD中说的弛豫？

[3]公众号“模拟小窝”：两相平衡体系模拟的建议

[4] GROMACS教程：构建双相体系：http://jerkwin.github.io/GMX/GMXtut-4/

[5] GROMACS教程之水中的溶菌酶：http://jerkwin.github.io/GMX/GMXtut-1

[6] 小木虫：关于几个温度控制方法的优缺点的讨论：http://muchong.com/html/200811/1074798.html

[7] 苑世领，分子模拟：理论与实验[M]

[8] GROMACS教程：DPPC膜中的KALP15：http://jerkwin.github.io/GMX/GMXtut-2

[9] http://manual.Gromacs.org/documentation/current/user-guide/mdp-options.html

[10] https://LAMMPS.sandia.gov/doc/fix_nh.html

[11]分子动力学模拟教程-多肽-蛋白相互作用:

http://jerkwin.github.io/2017/10/20/GROMACS分子动力学模拟教程-多肽-蛋白相互作用

[12]公众号“模拟小窝”：初学者数据处理注意事项

[13] http://bbs.keinsci.com/thread-14133-1-7.html

[14]Materials Modeling at the Molecular Level:

 http://www.carc.unm.edu/~acpineda/CS471/HTML/frink3/LECT3.html

[15]Lecture 3: The canonical ensemble:

https://chem.libretexts.org/Courses/New_York_University/CHEM-UA_652%3A_Thermodynamics_and_Kinetics/Lecture_03%3A_The_canonical_ensemble