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

本文旨在介绍利用PWmat与PLUMED的接口实现Metadynamics的功能。

什么是Metadynamics？

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

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

详细的原理可以参考doi: 10.1073/pnas.202427399

PLUMED

PLUMED是一个开源的方法库，里面包含了增强采样方法、自由能方法以及多种分析分子动力学（MD）轨迹的工具。

PWmat在MD的过程中嵌入了PLUMED的API接口，无需额外的安装即可实现Metadynamics的计算。

利用PWmat与PLUMED的接口探索CH3Cl的亲核取代反应

本例为CH3Cl的亲核取代反应，以PLUMED与PWmat的接口利用Metadynamics的功能来模拟一个带电的 Cl-取代CH3Cl中的另一个 Cl-。

结构文件

atom.config

结构构建了一个较大的胞来减少与相邻周期的相互作用，从而模拟孤立的分子反应。

输入文件

etot.input

在输入文件中，设置300K下NVT模拟50000步步长为1fs的分子动力学。在这里为了可以设置较长的步长，我们需要将H原子的质量增加为3.016 a.u.，方法是在赝势H.SG15.PBE.UPF中，在下图位置处添加

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

IN.MDOPT

plumed.input

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

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

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

下图节选了MOVEMENT轨迹中转变的阶段21100fs-22200fs，C-Cl断裂后与另一个Cl相连接的过程

利用PLUMED自带的sum_hills功能可以预估得到的自由能面（sum_hills功能需要自行安装PLUMED软件，参考Module73-plumed_pwmat接口使用示例，查看文档：https://www.pwmat.com/modulefiles/pwmat-resource/module-download7/pdf/guide_plumed_20240320.pdf）

$ 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处两者交叉并交换了位置，表明亲核反应的发生。