gromacs计算自由能

在 GROMACS 中可以使用“慢增长”方法计算两个分子物种之间的自由能差值。这种不同分子物种之间的自由能差值没有物理意义，但可以通过热力学循环利用它们获得有意义的量。

这种方法要求在模拟过程中，系统的哈密顿量缓慢地从描述一个系统（A）转变到描述另一个系统（B）。变化必须非常缓慢，以保证系统在整个过程中仍保持平衡状态。

如果满足了这个要求，变化就是可逆的，并且从 B 到 A 的慢增长模拟会与从 A 到 B 的慢增长模拟得到相同的结果（但符号相反）。

这是一个有用的核查方法，但用户应该意识到，正向和反向增长结果的相等并不能保证结果的正确性。

GROMACS 也可以在一次模拟中在从 A 到 B 的整个范围内进行积分(TI方法)。但是，如果变化很大，采样可能不充分，用户可能更愿意在一些精心选择的<dG/dλ>的值。

这很容易通过将步长 delta_lambda 设置为零来实现。每次模拟都可以先进行平衡，再λ中间值上精确地确定

根据φH/φλ的涨落得到每个dG/dλ值适当的误差估计。然后通过适当的数值积分方法确定总的自由能变化。

GROMACS 现在也支持利用 Bennett 接受率（BAR，Bennett’s Acceptance Ratio）方法计算从状态 A 转变到状态 B 的 ΔG 值，使用的程序为gmx bar。

利用相同的数据也可以使用 MBAR方法来计算自由能，但目前需要借助于外部 pymbar 程序包https://SimTK.org/home/pymbar 中的工具。

λ设置

所有常见势能类型和约束类型都可以从状态 A(λ = 0) 平滑地插值到状态 B(λ = 1)，反之亦然。所有成键相互作用的插值都是通过对相互作用参数的线性插值来实现的。

非键相互作用可以使用线性插值或软核相互作用。自 GROMACS 4.6 开始，λ是一个向量，允许自由能转变的不同分量以不同速率进行。库仑，Lennard Jones，成键和限制项都可以单独控制。

简谐势：这里给出的例子针对键势能, 在GROMACS中使用的是简谐势. 然而, 这些方程也同样适用于键角势与不当二面角势.

GROMOS–96键和键角：对四次键伸缩势与基于余弦的键角势, 通过对力常数与平衡位置的线性插值进行插值

正常二面角：对于正常二面角，方程稍微复杂一些：

注意：多重度无法参数化，因为函数在 [0, 2𝜋] 区间内应当保持周期性。

表格成键相互作用:对表格成键相互作用，只能对力常数进行插值

库仑相互作用:两个粒子之间的库仑相互作用, 粒子电荷随λ的变化关系如下:

反应场库仑相互作用：对两个粒子间包含反应场的库仑相互作用, 粒子电荷随λ的变化关系如下:

Lennard-Jones相互作用：对两个粒子间的Lennard-Jones相互作用, 原子类型随λ的变化关系可写为:

应该注意的是，从状态 A 到状态 B 的途径也可以使用σ和 ε来表示。看起来改变力场参数σ 和 ε比改变衍生参数C12和C6在物理上更合理。然而，在参数空间中这两种途径之间的差别并不大，而且自由能本身也不依赖于途径，所以我们使用上面的简单公式。

动能:当粒子的质量发生变化时，动能对自由能也有贡献 (注意，我们不能将动量 p写成mv，因为那样会导致dEk/dλ的符号错误

求导后，可以代入p=mv

约束：约束在形式上是哈密顿量的一部分，因此它们对自由能也有贡献。在 GROMACS 中可以使用 LINCS 或SHAKE 算法计算约束。对 LINCS，如果我们有k = 1 … K个约束方程 gk，那么

其中rk为两粒子间的位移向量，dk为两粒子间的约束距离。利用约束距离与 λ的依赖关系：

因此约束方程对λ的依赖关系为：

约束对哈密顿量的（零）贡献 G (使用拉格朗日乘子λk，它们逻辑上并不同于自由能λ) 为：

λ设置-软核势

在自由能计算中，当有粒子凭空出现或消失时，对 Lennard-Jones 势和库仑势使用方程所述的简单线性插值可能收敛性很差。当粒子几乎消失或接近出现时（λ接近 0 或 1），相互作用能会变得非常弱，粒子彼此可能非常靠近，导致得到的 dV/dλ值的波动非常大（因为简单的线性插值同时取决于势能在λ两个端点的值）。为避免这些问题，需要消除势能的奇异性。这可以通过对常规的 Lennard-Jones 势和库仑势进行“软核”修正来实现，软核势限制了λ 值介于 0 和 1 之间时的能量和力，但不会改变λ等于 0 或 1 时的能量和力。在 GROMACS 中，软核势 Vsc是常规势能的移位势，因此势能及其导数的奇点完全不可能出现在r = 0处：

其中VA和VB分别为状态 A(λ= 0) 和状态 B(λ = 1) 中正常的“硬核”范德华势或静电势，α为软核参数（在mdp文件中利用 sc_alpha 进行设置）; p为软核势中λ的次数（使用 sc_power 进行设置），σ为相互作用半径，其值为 (C12/C6)1/6。当 C6或 C12为零时，σ是一个输入参数（由 sc_sigma设置）。

自由能循环 A: 计算抑制剂 I 分别与酶 E，E′ 结合时的自由能差值 ΔG12。B: 计算抑制剂 I，I′ 分别与酶 E 结合时的自由能差值 ΔG12.