【龙讯module小课堂】Metadynamics：PWmat+PLUMED实现元动力学算例详解

在各个学科领域里分子动力学(MD)都已经产生了深远的影响。然而这个模拟方法的应用范围却一直受到制约，因为该方法所能探索的时间尺度有限，无法遍历所有的结构状态。其中一种情况就是自由能面具有很多个局部的最小值，而在它们之间却有一个较高的能垒难以越过，标准的AIMD就很难越过这个能垒达到另一个亚稳态。这时我们就需要一些特殊的采样方法来帮助我们度过这个能垒。

元动力学(Metadynamics, mMD)是一种在分子动力学模拟中增强采样的方法，它可以将自由能面转变为由几个自己选定的自由度的函数，这些自己选定的自由度通常叫做集 合变量(collective variable, CV)。在元动力学中，基于CV构建了的与历史相关的偏置势能(具体形式如下，由一系列离散的基于CV的高斯函数累加而成)可以帮助某一亚稳态的构型越过到另一个亚稳态的自由能垒，比起AIMD可以更容易观察到想要的过程。

结构文件

atom.config

输入文件

etot.input

除此之外为了调用PLUMED，需要设置IN.MDOPT=T，并准备一个MD的设置文件IN.MDOPT，在其中设置MD_USE_PLUMED=T。而关于Metadynamics的参数需要在plumed.input内设置。

IN.MDOPT

PLUMED的输入和输出的默认单位是能量:kJ/mol, 长度:nm, 时间:ps。如果需要改变可以通过参数UNITS来设置。

在PLUMED的输入文件内，先监控两个C-Cl的键长，即第一个原子和第五/六个原子的距离设为d1,d2，然后将ss=-1*d2+1*d1作为集 合变量(CV)，在此基础上累加高斯势推动集 合变量ss的变化。并且为了防止Cl与CH₃Cl的解离，需要设定一个集 合变量ss的上下限，在集 合变量ss的-0.3nm与0.3nm处添加一个势能墙。

下图为随时间变化的集 合变量ss的变化，可以看到随着时间推移，集 合变量ss在不断地震荡，经历不同的结构，在20ps处发现了转变的过程，可以看到C-Cl断裂后，与另一个Cl相连接的从反应物到产物的过程。

利用PLUMED自带的sum_hills功能可以预估得到的自由能面（sum_hills功能需要自行安装PLUMED软件，参考Module73-plumed_pwmat接口使用示例，可在文末点击阅读原文查看文档）

$ plumed sum_hills --hills HILLS

其中HILLS是储存高斯势的文件

想要观察到自由能面逐渐被填满的过程，可以在sum_hills的命令中添加--stride指令，分别输出高斯势沉积的过程中阶段性的高斯势，用于判断收敛的情况

$ plumed sum_hills --hills HILLS --stride 100

(需要注意的是这里的100指的是100次的沉积，与plumed.input中的METAD中PACE的意义需要做区分，指的是在MD_DETAIL中设置的50000步中每隔100*PACE(50)=5000步记录一次。)

然而这个例子比较特殊，在很多时候我们都很难用一个集 合变量来描述整个反应的变化。所以接下来我们尝试同时使用d1,d2（两个C-Cl的距离）作为集 合变量来做演示。其他参数都与之前相同，仅改变plumed.input中的设置

plumed.input

这里我们把METAD的ARG参数由ss改为d1、d2，注意SIGMA（展宽）参数需要为d1和d2分别赋值。同样为了防止它们之间的解离我们在d1,d2的0.5nm处加上一个势能墙。

在下图中可以看出两个键长在开始时各自在两侧振荡，到了约32ps处两者交叉并交换了位置，表明亲核反应的发生。