PLUMED MetaDynamics Simulation 输入文件编写教程

CocoCream

已于 2023-12-14 12:16:29 修改

阅读量982

点赞数 2

分类专栏： GROMACS 文章标签： Gromacs PLUMED

于 2023-12-01 19:14:06 首次发布

原文链接：https://www.plumed.org/doc-v2.4/user-doc/html/_m_e_t_a_d.html

版权

GROMACS 专栏收录该内容

19 篇文章 79 订阅

订阅专栏

本文翻译自 PLUMED 官方教程，旨在指导读者完成MetaDynamics Simulation 输入文件的编写，该文件可用于诸如 GROMACS 的gmx_mpi mdrun -plumed，以完成最终模拟。

Introduction

这是偏置模块的一部分，用于对一个或多个集合变量执行 MetaDynamics。在元动力学模拟中，由间歇性添加的高斯函数组成的历史相关偏置被添加到势能中：
$V(\vec{s},t) = \sum_{ k \tau < t} W(k \tau) \exp\left( -\sum_{i=1}^{d} \frac{(s_i-s_i^{(0)}(k \tau))^2}{2\sigma_i^2} \right).$ 这种势能迫使系统远离势能表面的动力学陷阱，进入能量景观中未探索的部分。关于构成该势的高斯函数的信息被输出到一个名为HILLS的文件，该文件用于重新开始计算和将自由能重建为 CV 的函数。自由能可以从元动力学计算中重建，因为最终偏差由下式给出： $V(\vec{s}) = -F(\vec(s)).$ 在后处理过程中，可以使用sum_hills程序以这种方式计算自由能。

在最简单的元动力学计算实现中，元动力学计算的费用随着模拟的长度而增加，因为在每一步都必须评估越来越多的高斯值。为了避免这个问题，您可以将偏差存储在网格上。这种方法类似于¹中提出的方法，具有网格间距与高斯宽度无关的优点。请注意，对于这种方法，您应该为每个集合变量（GRID_BIN）或所需的网格间距（GRID_SPACING）提供 bins 的数量。如果您同时提供这两个选项，则将对每个维度使用最保守的选择（二者中最大的 bins 数量值）。如果您没有提供任何关于 bins 大小的信息（既不是 GRID_BIN 也不是GRID_SPACING），并且如果高斯宽度是固定的，PLUMED 将使用高斯宽度的 1/5 作为网格间距。对于大多数程序来说，此默认选择应该是合理的。

您可以从HILLS文件和GRID重新启动元动力学，在第二种情况下，可以首先使用GRID_WFILE（和GRID_WSTRIDE）保存 GRID，然后在稍后阶段使用 GRID_RFILE 读取它。

在 PLUMED 中可用的另一个选项是回火元动力学（well-tempered metadynamics）²。在这种元动力学的变化中，高斯山丘的高度在每一步都会重新缩放，因此偏置现在由下式给出： $\sum_{t'=0,\tau_G,2\tau_G,\dots}^{t'<t} W e^{-V({s}({q}(t'),t')/\Delta T} \exp\left( -\sum_{i=1}^{d} \frac{(s_i({q})-s_i({q}(t'))^2}{2\sigma_i^2} \right)$ 这种方法可以确保偏差更平滑地收敛。应该注意的是，在回火元动力学的情况下，在输出中会使用偏置因子重新缩放高斯高度。还要注意的是，对于回火元动力学，HILLS文件不包含偏差，而是包含自由能估计的负值。这种选择的优点是可以使用不同的 ΔT 值重新启动模拟，所施加的偏置也将相应地缩放。

请注意，您也可以在这里使用 flexible gaussian approach，在该方法中，您可以将高斯势调整为变量覆盖的笛卡尔空间的范围，或调整为给定时间内覆盖的 CV 中的空间。在这种情况下，沉积的高斯势的宽度仅由一个值表示，该值是笛卡尔空间（ADAPTIVE=GEOM）或时间（ADAPTIVE=DIFF）。请注意，在这种情况下，应使用沉积高斯的特定积分技术，具体请查看文档中的程序sum_hills。

使用关键字INTERVAL，可以更改元动力学算法以设置偏置力为等于零的外部边界。例如，如果在 CV s 上执行元动力学，并且仅对 s > s_w 的自由能感兴趣，则根据上述方程仍然更新s > s_w 范围的历史相关势，但是对于 s < s_w，元动力学偏置力被设置为零。请注意，如果 s < s_w，也会添加高斯势，因为这些高斯势的尾部会影响相关区域 s > s_w 中的 VG。这样，s > s_w 区域系统上的力来自元动力学和常规力场，s < s_w 区域仅来自后者。这种方法允许获得在稳定估计器周围波动的历史相关偏置势 VG，该势等于离边界足够远的自由能的负值。请注意：

它只适用于一维偏置；
它既可以使用GRID，也可以不使用GRID；
在自由能导数不大的区域中的区间极限 s_w；
如果在限制 s_w 之外的区域中，系统具有最小自由能，则INTERVAL关键字应与 s_w 处的UPPER_WALLS或LOWER_WALLS一起使用。

最后要注意的是，自2.0.2版本以来，当系统在所选间隔之外时，力设置为零，偏置值设置为相应边界处的值。这允许正确计算副本交换（replica exchange）方法的接受度。

也可以使用多个walkers（可以理解成同时使用多个搜索点），请参阅后面的示例。

c(t) 重加权因子也可以使用 ³ 中给出的方程实时计算。用于计算 c(t) 的表达式直接来自于上文的方程3中的方程12，并且给出了比该文中的等效方程13和方程14更平滑的结果。c(t) 由 rct 分量给出，而由c(t) 通过使用关键字REWEIGHTING_NGRID来启用 c(t) 的计算，其中指定了用于计算的网格。此网格的大小应等于或大于GRID_BIN./中给定的网格。默认情况下，c(t) 每50个高斯山丘更新一次，但这可以通过使用REWEIGHTING_NHILLS关键字来更改，此选项只能与回火元动力学模拟一起使用，并且需要使用网格。

请注意，与 PLUMED 1.3 不同的是，现在可以直接应用并行元动力学，这确实可以通过在同一输入文件中多次使用 METAD 操作来取得。

Description of Components

**默认情况下元动力学将计算以下量，**你可以通过在元动力学标签后面跟着一个点和下面列表中量的名称来引用相应的量（例如 meta.bais 表示标签为meta的元动力学的瞬时偏置势）。

bais：偏置势的瞬时值；
work：累计功；

此外，可以通过使用下面列出的关键字来计算以下量：

rbias：使用REWEIGHTING_NGRID关键词。使用 c(t) 重加权因子 [ rbias＝bias-c(t) ] 归一化的偏置势的瞬时值。该分量可用于获得重新加权的直方图；
rct：使用REWEIGHTING_NGRID关键词。重加权因子 c(t) ；
acc：使用ACCELERATION关键词。元动力学加速度因子；
maxbias：使用CALC_MAX_BIAS关键词。元动力学 V(s，t) 的最大值；
transbias：使用CALC_TRANSITION_BIAS关键词。元动力学转变偏置 V*(t)；

以下为必填关键字（括号内为默认值）：

SIGMA：高斯丘陵的宽度；
PACE：斜坡加法的频率；
FILE（=HILLS）：存储已添加山丘列表的文件；

以下为可选关键字（括号内为默认值）：

NUMERICAL_DERIVATIVES ( default=off )：以数字方式计算这些量的导数；
GRID_SPARSE ( default=off )：使用稀疏栅格来存储山丘；
GRID_NOSPLINE ( default=off )：不将样条曲线插值与网格一起使用；
STORE_GRIDS ( default=off )：存储计算生成的所有网格文件。如果此关键字不存在，它们将被删除；
WALKERS_MPI ( default=off )：使用MPI版本的多个步行器。与WALKERS_DIR以外的WALKERS_*选项不兼容；
ACCELERATION ( default=off )：如果要计算元动力学加速度因子，请设置为TRUE；
CALC_MAX_BIAS ( default=off )：如果要计算元动力学 V(s，t) 的最大值，请设置为TRUE；
CALC_TRANSITION_BIAS ( default=off )：如果要计算元动力学转换偏置 V*(t)，请设置为TRUE；

ARG：此关键字的输入是其他标签的标量输出，您需要在ARG关键字后写明你希望使用的标量的标签。如果标签单独出现，则假定原标签计算单个标量值。如果出现*或*.*，则输入标签中所有计算的标量均作为输入；这种情况可能发生在某些标签具有多个输出，但一般输出的每个值都有一个额外的特定标签；例如，一个名为 dist 的DISTANCE标签可能有三个输出值x、y和z，如果你希望只使用x部分，您应该使用dist.x，如果您想同时使用这三个输出值，请使用dist.*。有关该标签的更多信息，请参阅手册的PLUMED入门部分。该标签的输入也可以使用POSIX的正则表达式，详见正则表达式一节。要使用此功能，您必须使用适当的标志编译PLUMED。您可以为该关键字添加多个输入，即 ARG1、ARG2、ARG3；

HEIGHT：高斯山丘的高度。必须选项，除非给出TAU和BIASACTOR或DAMPACTOR；

FMT：指定HILLS文件的格式（有助于减少regtest中的位数）；

TARGET：指向预定义的分布；

RECT：所有副本的偏置因子的列表；

DAMPFACTOR：以exp(-max(V)/(kbT*DAMPFACTOR)进行沉积；

BIASFACTOR：偏置因子，使用回火元动力学时启用的关键字。请注意，您还必须指定TEMP；
TEMP：系统温度，只有进行回火元动力学时才需要；
TAU：在回火元动力学中，设置高斯高度为(kb*DeltaT*pace*timestep)/TAU

TTBIASFACTOR：偏置因子，使用转温元动力学（transition tempered metadynamics）时启用的关键字。请注意，您还必须指定TEMP；
TTBIASTHRESHOLD：偏置阈值，使用转温元动力学时启用的关键字。请注意，您还必须指定TTBIASFACTOR；
TTALPHA：山丘大小衰减指数，使用转温元动力学时启用的关键字。请注意，您还必须指定TTBIASFACTOR；
TRANSITIONWELL：在转温元动力学中使用，此关键字作为TRANSITIONWELLx重复出现，其中x=0,1,2，…，n；每个都指定了一个井的坐标，。必须至少提供一个。您可以使用此关键字的多个实例，即TRANSITIONWELL1、TRANSITIONWILL2、TRANSITION WELL3…

GRID_MIN：网格的最小边界；
GRID_MAX：网格的最大边界；
GRID_BIN：网格的窗口数量；
GRID_SPACING：网格的窗口大小，在效果上等效于GRID_BIN；
REWEIGHTING_NGRID：计算 c(t) 重加权因子，并使用该因子来获得归一化偏置 [ rbias＝bias-c(t) ]。在此，您应该指定每个维度中所需的网格的数量。网格点的数量应等于或大于GRID_BIN中给定的网格数量。此方法与不使用网格的元动力学不兼容。
INTERVAL：一维的下限和上限，超出该阈限系统将不会使用偏置；

REWEIGHTING_NHILLS：计算 c(t) 重新加权因子的高斯山丘间隔。默认情况下，每50个山丘执行一次此操作。
GRID_WSTRIDE：每隔N步将网格写入文件
GRID_WFILE：写入网格的文件
GRID_RFILE：一个网格文件，用于在模拟的初始步骤读取偏差

ACCELERATION_RFILE：在模拟的初始步骤读取加速度的数据文件；

ADAPTIVE：使用基于几何（=GEOM）或扩散（=DIFF）的山丘宽度方案。西格玛（标准差）是一个具有距离单位或时间步长维度的数字；
SIGMA_MAX：使用自适应山丘时SIGMA的上限（以CV为单位）。负数表示无边界
SIGMA_MIN：使用自适应山丘时SIGMA的下限（以CV为单位）。负数表示无边界

WALKERS_ID：walker ID；
WALKERS_N：walker 数量；
WALKERS_DIR：所有不同 walker 产生的 HILLS 文件的共享文件夹；
WALKERS_RSTRIDE：读取 HILLS 文件的步幅；

RESTART：模拟运行标签前的关键词，可设置是否为重新运行（可以理解为续跑）
UPDATE_FROM：仅更新此时间之前的模拟结果；
UPDATE_UNTIL：仅更新此时间之后的模拟结果；

Examples

Examples 1.

以下输入用于使用原子3和5之间的距离以及原子2和4之间的距离作为集合变量（CV）的标准元动力学计算。CV 的值和元动力学偏置势每 100 步写入COLVAR文件：

DISTANCE ATOMS=3,5 LABEL=d1
DISTANCE ATOMS=2,4 LABEL=d2
METAD ARG=d1,d2 SIGMA=0.2,0.2 HEIGHT=0.3 PACE=500 LABEL=restraint
PRINT ARG=d1,d2,restraint.bias STRIDE=100  FILE=COLVAR

Examples 2.

如果你使用扩散方案的自适应高斯，高斯势应该覆盖集体变量中20个时间步长的空间。请注意，在这种情况下，对山丘求和时需要进行直方图校正：

DISTANCE ATOMS=3,5 LABEL=d1
DISTANCE ATOMS=2,4 LABEL=d2
METAD ARG=d1,d2 SIGMA=20 HEIGHT=0.3 PACE=500 LABEL=restraint ADAPTIVE=DIFF
PRINT ARG=d1,d2,restraint.bias STRIDE=100  FILE=COLVAR

如果你使用几何方案的自适应高斯，高势斯应该覆盖笛卡尔空间中0.05nm的空间。请注意，在这种情况下，对山丘求和时需要进行直方图校正：

DISTANCE ATOMS=3,5 LABEL=d1
DISTANCE ATOMS=2,4 LABEL=d2
METAD ARG=d1,d2 SIGMA=0.05 HEIGHT=0.3 PACE=500 LABEL=restraint ADAPTIVE=GEOM
PRINT ARG=d1,d2,restraint.bias STRIDE=100  FILE=COLVAR

使用自适应高斯时，可能需要限制山丘宽度的变化方式。您可以使用SIGMA_MIN和SIGMA_MAX关键字。应根据 CV 指定 sigmas，并使用 CV 的单位。请注意，如果使用负数，则表示不设置限制。还要注意的是，在这种情况下，对山丘求和时需要进行直方图校正：

DISTANCE ATOMS=3,5 LABEL=d1
DISTANCE ATOMS=2,4 LABEL=d2
METAD ...
  ARG=d1,d2 SIGMA=0.05 HEIGHT=0.3 PACE=500 LABEL=restraint ADAPTIVE=GEOM
  SIGMA_MIN=0.2,0.1 SIGMA_MAX=0.5,1.0
... METAD
PRINT ARG=d1,d2,restraint.bias STRIDE=100  FILE=COLVAR

Examples 3.

可以使用多个 walkers，如[75]中所述。通过设置使用的 walkers 数量、解释输入文件的当前 walkers 的id、包含文件的山丘所在的目录以及读取其他助行器时的频率，来启用这些助行器。下面是一个例子：

DISTANCE ATOMS=3,5 LABEL=d1
METAD ...
   ARG=d1 SIGMA=0.05 HEIGHT=0.3 PACE=500 LABEL=restraint
   WALKERS_N=10
   WALKERS_ID=3
   WALKERS_DIR=../
   WALKERS_RSTRIDE=100
... METAD

其中WALKERS_N是 walkers 的总数，WALKERS_ID是当前步行者的ID（从0开始），WALKERS_DIR是所有 walkers 所在的目录。WALKERS_RSTRIDE是一次更新和另一次更新之间的步数。从2.2.5版本开始，hills 文件在WALKERS_RSTRIDE的每一步都会自动刷新。

Examples 4.

c(t) 重加权因子可以使用上述[85]中提出的方程动态计算。这是通过使用关键字REWEIGHTING_NGRID来启用的，其中设置了用于计算的网格。REWEIGHTING_NGRID中给出的格点数量应等于或大于GRID_BIN中给出的网格点数量。

phi: TORSION ATOMS=5,7,9,15
psi: TORSION ATOMS=7,9,15,17
METAD ...
 LABEL=metad
 ARG=phi,psi SIGMA=0.20,0.20 HEIGHT=1.20 BIASFACTOR=5 TEMP=300.0 PACE=500
 GRID_MIN=-pi,-pi GRID_MAX=pi,pi GRID_BIN=150,150
 REWEIGHTING_NGRID=150,150
 REWEIGHTING_NHILLS=20
... METAD

在这里，我们要求使用REWEIGHTING_NHILLS关键字在每20个山丘进行一次计算。如果未给定此关键字，则默认情况下每50个山丘执行一次计算。c(t) 重加权因子将在 rct 分量中给出，而使用 c(t) 重新加权因子归一化的偏置势的瞬时值在 rbias 分量 [rbias=bias-c(t)] 中给出，该分量可用于 HISTOGRAM 分析获得重新加权的直方图或自由能面。

盆地之间转换的动力学也可以像[84]中那样进行动态分析。标签ACCELERATION开启加速度因子的累积，然后可用于确定速率。该方法可与COMMITTOR分析一起使用，以在系统到达目标盆地时停止模拟。它必须与回火元动力学一起使用。

您还可以使用关键字TARGET提供目标分发。TARGET应为包含自由能（即所需目标分布的 -kbT*log）的网格。然后，高斯将按该因子缩放 $e^{\beta(\tilde{F}(s)-\tilde{F}_{max})}$ 这里，F 是在网格上定义的自由能，Fmax 是它的最大值。请注意，我们在这里使用了参考文献[46]中的最大值。这个选择可以避免添加过大的高斯势。然而，它可能会使高斯势过小。在这种情况下，您应该始终仔细选择HEIGHT参数。网格文件应与其他 PLUMED 网格文件类似，因为它应包含目标自由能及其导数。

请注意，如果希望模拟收敛到目标自由能，则应使用DAMPACTOR命令提供全局回火。或者，如果使用BIASACTOR，则模拟将收敛到另一个自由能，该自由能是目标自由能和由原始力场确定的固有自由能的线性组合。

DISTANCE ATOMS=3,5 LABEL=d1
METAD ...
 LABEL=t1
 ARG=d1 SIGMA=0.05 TAU=200 DAMPFACTOR=100 PACE=250
 GRID_MIN=0 GRID_MAX=2 GRID_BIN=200
 TARGET=dist.dat
... METAD

PRINT ARG=d1,t1.bias STRIDE=100 FILE=COLVAR

此计算的dist.dat文件中的标头将为：

#! FIELDS d1 t1.target der_d1
#! SET min_d1 0
#! SET max_d1 2
#! SET nbins_d1  200
#! SET periodic_d1 false

请注意，BIASACTOR也可以选择为等于1。在这种情况下，将执行无偏采样。此时应使用TAU参数，而不是使用HEIGHT。例如：

d: DISTANCE ATOMS=3,5
METAD ARG=d SIGMA=0.1 TAU=4.0 TEMP=300 PACE=100 BIASFACTOR=1.0

获得的HILLS文件仍像 PLUMED 的sum_HILLS同样可用，以便绘制自由能。有意义的情况可能是 RECT 模拟。而关于 RECT 模拟，您也可以使用RECT关键字，以避免使用多个输入文件。例如，单个输入文件的情况可能像这样，其中RECT后的元素个数应等于副本数：

d: DISTANCE ATOMS=3,5
METAD ARG=d SIGMA=0.1 TAU=4.0 TEMP=300 PACE=100 RECT=1.0,1.5,2.0,3.0

在这里插入图片描述

V. Babin, C. Roland, and C. Sagui. Adaptively biased molecular dynamics for free energy calculations. J. Chem. Phys., 128:134101, 2008. ↩︎
A Barducci, G Bussi, and M Parrinello. Well-tempered metadynamics: A smoothly converging and tunable free-energy method. Phys. Rev. Lett., 100(2):020603, Jan 2008. ↩︎
Pratyush Tiwary and Michele Parrinello. A time-independent free energy estimator for metadynamics. The Journal of Physical Chemistry B, 119(3):736–742, Jan 2015. ↩︎