Gromacs结合自由能计算教程（四）

以下文章来源于网络，如有侵权,请联系管理员删除!

本教程使用的分析脚本文献：

J Comput Aided Mol Des (2015) 29: 397. 

https://doi.org/10.1007/s10822-015-9840-9

完成了之前的模拟，我们得到了一系列dH/dλ文件，接下来就是进行分析来得到自由能了；分析的方法又有很多种，但总体而言可归为两类——自由能微扰（FEP）和热力学积分法（TI），这一步的理论对于新手来说可能会显得艰深。在这里我强烈建议初学者优先学习热力学积分法，这可以比较直观地理解并利用dH/dλ文件的意义。

· 首先从自由能的定义开始，Q为配分函数

· 对λ求导，得到

· 改写为

· 最后，可以在整个λ范围内进行积分，以得到最终的TI方程

当体积功忽略不计时（计算结合自由能的大部分情况），dU=dH，所以dA=dG。

从以上推导可以看出，实际上的自由能就是对dH/dλ进行积分。而实际上由于选取的λ是离散值，最终的积分也是近似的，亦即：

以上所说的看上去很复杂对不对？不过实际上你完全不懂这些也不要紧，因为已经有人开发了一键分析的python脚本alchemical-analysis，只需要配置好脚本后在dHdl文件夹内运行

alchemical-analysis -t 温度 -s 平衡时间 等参数

即可，自动分析并调用一系列方法（可选）得到误差和最终自由能

如果按照默认的分析参数，以及在mdp中设置了calc-lambda-neighbors = -1，那么会得到以下一系列方法的分析结果

# Command line was: /bin/alchemical_analysis -t 300 -s 0 -u kcal -w -g

------------ --------------------- --------------------- --------------------- --------------------- --------------------- ---------------------

   States           TI (kcal/mol)   TI-CUBIC (kcal/mol)       DEXP (kcal/mol)       IEXP (kcal/mol)        BAR (kcal/mol)       MBAR (kcal/mol)

------------ --------------------- --------------------- --------------------- --------------------- --------------------- ---------------------

   0 -- 1         2.006  +-  0.012      2.009  +-  0.024      1.922  +-  0.024      1.915  +-  0.017      1.923  +-  0.011      1.923  +-  0.011

   1 -- 2         4.961  +-  0.088      5.493  +-  0.101      6.041  +-  0.078      5.963  +-  0.418      5.691  +-  0.058      5.690  +-  0.058

   2 -- 3         3.350  +-  0.094      3.504  +-  0.106      3.144  +-  0.513      1.967  +-  0.192      3.043  +-  0.169      3.043  +-  0.170

   3 -- 4         0.033  +-  0.044     -0.263  +-  0.053     -0.251  +-  0.206     -0.389  +-  0.150     -0.108  +-  0.051     -0.108  +-  0.051

   4 -- 5        -1.261  +-  0.064     -1.243  +-  0.073     -1.046  +-  0.077     -1.560  +-  0.191     -1.227  +-  0.062     -1.252  +-  0.061

   5 -- 6        -1.952  +-  0.059     -2.030  +-  0.068     -1.947  +-  0.136     -1.946  +-  0.097     -1.960  +-  0.049     -1.925  +-  0.045

   6 -- 7        -2.120  +-  0.065     -2.086  +-  0.076     -2.298  +-  0.043     -1.372  +-  0.349     -2.244  +-  0.033     -2.287  +-  0.032

   7 -- 8        -2.406  +-  0.067     -2.403  +-  0.075     -2.436  +-  0.126     -2.499  +-  0.171     -2.515  +-  0.055     -2.608  +-  0.038

   8 -- 9        -1.382  +-  0.013     -1.382  +-  0.014     -1.371  +-  0.021     -1.419  +-  0.019     -1.388  +-  0.013     -1.423  +-  0.010

   9 -- 10       -1.469  +-  0.008     -1.476  +-  0.009     -1.426  +-  0.014     -1.500  +-  0.009     -1.493  +-  0.006     -1.469  +-  0.005

  10 -- 11       -1.507  +-  0.005     -1.507  +-  0.008     -1.532  +-  0.005     -1.471  +-  0.016     -1.520  +-  0.004     -1.492  +-  0.003

  11 -- 12       -0.602  +-  0.003     -0.602  +-  0.003     -0.600  +-  0.004     -0.603  +-  0.004     -0.602  +-  0.003     -0.601  +-  0.001

  12 -- 13       -0.905  +-  0.003     -0.908  +-  0.004     -0.913  +-  0.005     -0.898  +-  0.002     -0.900  +-  0.002     -0.906  +-  0.001

  13 -- 14       -1.370  +-  0.025     -1.440  +-  0.032     -1.157  +-  0.006     -1.340  +-  0.122     -1.194  +-  0.008     -1.116  +-  0.004

  14 -- 15       -2.153  +-  0.025     -2.194  +-  0.028     -2.635  +-  0.096     -1.577  +-  0.069     -1.862  +-  0.014     -1.843  +-  0.005

  15 -- 16       -5.884  +-  0.075     -5.713  +-  0.102     -7.271  +-  0.395     -4.572  +-  0.281     -6.370  +-  0.036     -5.308  +-  0.012

  16 -- 17       -3.943  +-  0.038     -3.946  +-  0.038     -4.185  +-  0.084     -3.171  +-  0.050     -3.335  +-  0.018     -3.723  +-  0.006

  17 -- 18       -4.701  +-  0.025     -4.686  +-  0.026     -5.253  +-  0.021     -4.168  +-  0.169     -5.148  +-  0.011     -5.101  +-  0.009

  18 -- 19      -12.421  +-  0.075    -12.289  +-  0.091    -12.869  +-  0.330    -10.787  +-  0.467    -12.284  +-  0.072    -13.453  +-  0.043

  19 -- 20      -17.461  +-  0.077    -17.403  +-  0.088    -18.187  +-  0.241    -15.696  +-  0.526    -17.415  +-  0.081    -17.454  +-  0.081

------------ --------------------- --------------------- --------------------- --------------------- --------------------- ---------------------

  Coulomb:      -47.933  +-  0.190    -47.670  +-  0.213    -51.556  +-  0.583    -41.310  +-  0.790    -47.609  +-  0.117    -47.997  +-  0.105

  vdWaals:       -3.254  +-  0.265     -2.893  +-  0.264     -2.713  +-  0.595     -3.813  +-  0.658     -3.298  +-  0.213     -3.415  +-  0.234

    TOTAL:      -51.187  +-  0.326    -50.563  +-  0.339    -54.269  +-  0.833    -45.123  +-  1.028    -50.907  +-  0.243    -51.413  +-  0.256

对于一个稳健的模拟，各种分析方法得到的结果一般不会相差太大

另外，你还会得到一个OMBAR邻接矩阵

邻接矩阵反映出某λ值下的模拟与相邻λ值下的模拟的相空间重叠程度，对于一个稳健的模拟，其邻接矩阵的主对角线和相邻两条对角线上的值都不应该为0。

另外一个文件便是以热力学积分形式显示的积分图：

使用离散值近似积分运算，如上图所示；其中范德华力部分和库仑力部分分别积分。

原则上说λ越多积分越准确，但耗时也将大大增加；λ太少，或间距太大，相空间不能很好重叠，会导致积分结果与真实值误差大。因此选取好适宜的λ间距值很重要。眼尖的读者这时候可能看到了以上图片中λ间距值并不等距，这倒并不是说就应该这么取，实际的λ间距值应该考虑你的体系。

简单而言，如果某个λ值附近的曲率很大，那么离散近似带来的误差就会很大，这时候应该在这附近多取几个间距小的λ值以降低误差；反之，如果曲率较小，那么可以少取几个λ值以减少模拟用时。一般而言这种大曲率的情况在范德华力积分中容易出现，在库仑力中较少，因此库仑力对应的总λ值数目也较少。你可能需要进行多次模拟才能找到合适的λ值的分布。

另外，一般而言，为了得到可信的数据，一个模拟需要运行3-5个副本，得到的数据才具有较好的统计学可靠性。