质心牵引-介绍

当两组之间的距离不断变化时，力会作用在整个体系上，这意味着体系不再处于平衡。管在非常缓慢的拉动极限下，系统再次处于平衡状态，但对于许多系统来说，在合理的计算时间内无法达到这一极限。然而，人们可以使用Jarzynski关系从许多非平衡模拟中获得两个距离之间的平衡自由能差：

其中，所执行的工作是迫使系统沿着从状态A到B的一条路径移动，角括号表示在初始状态A和的规范系综上的平均值。

pull代码对集合变量（有时称为反应坐标）施加力或约束。基本的集体拉力坐标是原子群质心之间的距离、角度和二面角，即所谓的“拉力群”。更复杂的集合变量可以使用变换拉坐标来定义。

拉组可以是一个或多个拉坐标的一部分。此外，坐标还可以对空间中的单个组和绝对参考位置进行操作。一对组之间的距离可以在1、2或3维上确定，也可以沿着用户定义的向量确定。

参考距离可以是恒定的，也可以随时间线性变化。

通常，所有原子都按其质量加权，但也可以使用额外的加权因子。

用伞形采样拉出脂质双层的示意图。Vrup是弹簧缩回的速度，Zlink是弹簧所连接的原子，Zspring是弹簧的位置。

gromacs支持几种不同的拉力类型，即施加拉力的方式，在所有情况下，参考距离可以是恒定的或随时间线性变化的。

Umbrella pulling：在两组质心之间施加一个谐波势。因此，力与位移成正比。

Constraint pulling：约束拉伸两组质心之间的距离受到约束。约束力可以写入文件。该方法使用SHAKE算法，但如果只有两个组受到约束，则只需要1次迭代即可精确。

Constant force pulling：恒拉力在两组质心之间施加恒定力。因此，势能是线性的。在这种情况下，没有拉力的参考距离。

Flat bottom pulling：平底拉动类似于伞形拉动，但对于低于（）或高于（）参考值的坐标值，势能和力为零。

这有助于将两个分子之间的距离限制在某个区域pull-coord?-type = flat-bottom            pull-coord?-type = flat-bottom-high

External potential ：外部势能这从另一个模块获取作用于反应坐标的势能。目前仅支持加速权重直方图方法（见第节AWH自适应偏置），该方法提供拉力坐标的自适应偏置。

此外，还有不同类型的反应坐标，即所谓的拉力几何。这些是用mdp选项设置的。pull-coord?-geometry

定义拉动的方向

最常见的设置是沿着包含两个拉力组的矢量的方向拉力，这是通过选择的。您可能希望拖动某个向量，该向量是用选择的。但这可能会在系统中产生不必要的扭矩，除非你以（几乎）固定的方向拉动参考组，例如沿x/y方向嵌入膜中的膜蛋白，同时沿z方向拉动。

如果你的参考组没有固定的方向，你可能应该使用，见图40。由于势能现在取决于定义方向的另外两个组的坐标，因此扭矩将作用于这两个组。

pull-coord?-geometry = distancepull-coord?-geometry = directionpull-coord?-geometry = direction-relative

图40几何形状的拉力设置。“正常”拉力组为1和2。第3组和第4组定义了拉力方向，从而定义了法向拉力的方向（红色）。这导致组3和4上的反作用力（蓝色）垂直于拉力方向。它们的大小由“正常”拉力乘以第3组和第4组之间的比率和距离给出。

角度和二面角拉力几何形状的定义

需要四个拉力组。与具有两个组的几何体相同，每个连续的组对都定义了一个连接组和的COM的向量。角度定义为两个结果向量之间的角度。例如，mdp选项定义了从组1的COM到组2的COM的向量与从组2的OM到组4的COM的矢量之间的角度。

角度取闭合区间[0，180]度内的值。因为角度是相对于由给出的参考轴定义的，只需要给出两组。

二面角几何形状需要六个拉力组。这些以与上述相同的方式配对，因此定义了三个向量。

二面角定义为由两个第一矢量和最后两个矢量跨越的两个平面之间的角度。等效地，二面角可以看作是当这些矢量投影到垂直于第二矢量（轴矢量）的平面上时，第一和第三矢量之间的角度。例如，考虑一个涉及四组的二面角：1、5、8和9。这里，mdp选项指定了定义二面角的三个向量：第一个向量是从组8到1的COM距离向量，第二个向量是组1到5的COM距离矢量，第三个向量是第5到9的COM距离变量。二面角取值范围为（-180  180]度，具有周期性边界。

pull-coord?-geometry = angle

pull-coord?-groups = 1 2 2 4

pull-coord?-geometry = angle-axis

pull-coord?-vecpull-coord?-groups = 8 1 1 5 5 9

变换拉坐标

这里有两个使用转换坐标设置的mdp输入拉取部分的示例。第一个是接触反作用坐标，在接触时为1，在较大距离时为0：

pull                     = yes

pull-ngroups             = 2

pull-ncoords             = 2

pull-group1-name         = groupA

pull-group2-name         = groupB

pull-coord1-type         = umbrella

pull-coord1-geometry     = distance

pull-coord1-groups       = 1 2

pull-coord1-dim          = Y Y Y

pull-coord1-k            = 0      ; avoid forces working directly on this distance

pull-coord2-type         = umbrella

pull-coord2-geometry     = transformation

pull-coord2-expression   = 1/(1 + exp(50*(x1 - 1.8*0.3)))  ; x1 refers to the value of coord1

pull-coord2-init         = 1      ; this restrains the distance to having the contact

pull-coord2-k            = 100

第二个例子是两个距离的平均值：

pull                     = yes

pull-ngroups             = 4

pull-ncoords             = 3

pull-group1-name         = groupA

pull-group2-name         = groupB

pull-group3-name         = groupC

pull-group4-name         = groupD

pull-coord1-type         = umbrella

pull-coord1-geometry     = distance

pull-coord1-groups       = 1 2

pull-coord1-dim          = Y Y Y

pull-coord1-k            = 0      ; avoid forces working directly on this distance

pull-coord2-type         = umbrella

pull-coord2-geometry     = distance

pull-coord2-groups       = 3 4

pull-coord2-dim          = Y Y Y

pull-coord2-k            = 0      ; avoid forces working directly on this distance

pull-coord3-type         = umbrella

pull-coord3-geometry     = transformation

pull-coord3-expression   = 0.5*(x1 + x2)  ; x1 and x2 refer to the value of coord1 and coord2

pull-coord3-init         = 0.8    ; restrains the average distance to 0.8 nm

pull-coord3-k            = 1000

局限性

有一个理论限制：严格来说，约束力只能在没有通过约束连接到系统其余部分的组之间计算。

如果一个基团包含键长受约束的分子的一部分，则应同时迭代拉力约束和LINCS或SHAKE键约束算法。GROMACS中没有这样做。这意味着，对于mdp文件中的模拟，严格来说，拉取仅限于整个分子或分子组。

在某些情况下，可以通过使用自由能代码来避免这种限制，使用自由能编码计算PMF。

在实践中，当拉力群由大量原子组成和/或拉力很小时，不迭代两种约束算法造成的误差可以忽略不计。在这种情况下，与键长长度相比，拉力组的约束校正位移较小。

constraints = all-bonds