本文摘要(由AI生成):
这篇文章介绍了如何使用 PWmat 与 PLUMED 的接口实现 Metadynamics 功能。Metadynamics 是一种增强采样方法,可以帮助我们克服分子动力学模拟中的能垒,实现对自由能面的更全面探索。文章首先介绍了 Metadynamics 的基本原理,然后详细介绍了 PLUMED 这个开源库,它包含了增强采样方法、自由能方法以及多种分析 MD 轨迹的工具。接下来,文章通过实例演示了如何使用 PWmat 与 PLUMED 的接口进行 Metadynamics 计算,以模拟带电的 Cl(-)取代 CH(3)Cl 中的另一个 Cl(-)的亲核取代反应。最后,文章提供了一些 PLUMED 的使用技巧,如利用 sum_hills 功能预估自由能面等。
本文旨在介绍利用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接口使用示例,可在文末点击阅读原文查看文档)
$ 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处两者交叉并交换了位置,表明亲核反应的发生。