首页/文章/ 详情

Gromacs结合自由能计算教程(四)

1月前浏览1406

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

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

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个副本,得到的数据才具有较好的统计学可靠性。




来源:模拟之家
python理论GROMACS
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2024-09-08
最近编辑:1月前
刘十三613
博士 分子动力学、GROMACS
获赞 133粉丝 95文章 82课程 29
点赞
收藏
作者推荐
未登录
还没有评论
课程
培训
服务
行家
VIP会员 学习 福利任务 兑换礼品
下载APP
联系我们
帮助与反馈