gmx_mmpbsa使用说明

• 2019-07-31 10:14:58 初稿
• 2019-09-19 10:07:02 修订

缘起

很多年后, 当我再次更新这个工具, 就会记起, 很多年前, 我开始研究GROMACS中MMPBSA实现方法的情境. 那时的我就已经发现已有了两种解决方式: GMXPBSA脚本与g_mmpbsa程序. 我粗略地测试一下, 还算可用, 但未能尽美, 可改进之处也颇多. 因为我不经常使用, 便也就无所谓了, 暴虎冯河, 用过即忘, 倒也合适.

时间流逝, 形势在变, 心境也在变. 当我觉得真正需要测试一下MMPBSA方法时, 又回想起这两种解决方式, 重新捡拾起先前的操作, 才发现颇不顺手. GMXPBSA对GROMACS版本有要求, 脚本也颇得意大利人的风骨, 一如Italic的中文, 让我厌烦. 出起错来却又如意大利面, 纠纠结结, 纠缠不清, 让人无可措手, 只想烈马快刀, 一斩了之. g_mmpbsa呢, 印度人的作品, 传至东土已有些时日, 受者众多, 却仍有水土不服之虞. 安装麻烦, 二进制包Windows下也无法使用. 外加只能使用特定版本的GROMACS和APBS, 而无法兼容最新版本. 我本想好好参阅代码真经, 希冀修饬一二, 传与他人, 也算得我功德一件. 可总不得五蕴皆空, 六根清净. 既已如此, 又何须执著? 倒不如我也来写一段, 看看传得传不得.

简介

gmx_mmpbsa用于计算GROMACS轨迹的MMPBSA结合能, 并进行能量分解. 计算时, MM部分由脚本自行完成, PBSA部分借助APBS程序完成. 因此, 在使用gmx_mmpbsa前, 必需安装好GROMACS和APBS程序. 如果你已经需要使用这个脚本了, 那说明你已经得到了GROMACS的轨迹, GROMACS自然是已经安装好了的. 而APBS的安装也很简单, 官方网站提供了二进制版本, Linux和Windows下的都有, 拖过来就可以用.

gmx_mmpbsa计算MM和PBSA所需的原子参数(电荷, 半径, LJ参数等)来自GROMACS的.tpr文件, 因此, 需要调用gmx dump程序处理.tpr文件, 以便抽取所需的参数. 采用这种方法, 可以避免版本不兼容的问题, 因此也就可以支持任意版本的GROMACS. 对于APBS也是一样, 只要APBS的输入文件格式没有变化, 那gmx_mmpbsa生成的APBS输入文件就可以用于任意版本的APBS. 这样我们就可以安装GROMACS和APBS的最新版本而不用担心兼容性问题.

gmx_mmpbsa还使用了gmx trjconv程序, 用于处理轨迹的周期性叠合问题. 当然, 这一步你可以自己完成, 而无需借助脚本, 特别是在需要对周期性进行特殊处理的情况下.

在Windows下使用时, 还需要一个bash环境, 你可以使用msys2, Git Bash, 或cygwin.

使用流程

1. 获得输入文件

1. 生成tpr文件, 预平衡
2. 运行成品模拟, 得到xtc轨迹文件
3. 生成ndx文件, 其中需要定义三个组: 复合物(com), 蛋白(pro), 配体(lig). 名字虽然是pro, lig, 但其实可以代表任意分子, 比如两个有机小分子, 而不一定就是蛋白和配体.
4. 处理轨迹. 如果存在二聚体, 团簇等情况, 确保组成的原子间没有分离. 其他情况脚本可自动处理.

2. 设定gmx_mmpbsa计算参数

程序路径

• apbs='c:/apbs1.5/bin/apbs.exe' # APBS: APBS可执行文件的完整路径. Linux系统下如: apbs='/home/users/APBS/bin/apbs'.

基本参数

• ff=AMBER: 力场类型, AMBER, OPLS, CHARMM
• trj=1EBZ.xtc: 轨迹文件
• tpr=1EBZ.tpr: tpr文件
• ndx=index.ndx: 索引文件
• com=System: 复合物索引组
• pro=Protein: 蛋白索引组
• lig=BEC: 配体索引组

将trj, tpr, ndx等变量指向相应文件, 将com, pro, lig指向索引文件中相应的索引组.

其他参数主要是APBS计算所用的极性参数, 非极性参数, 网格参数, 一般无需更改太多.

3. 运行脚本

1. 脚本首先处理轨迹: 1. 完整化; 2. 居中叠合. 之后将构型输出到pdb文件.
2. 脚本其次抽取tpr中的原子信息, 存放在qrv文件中. 主要是复合物中每个原子的电荷, 半径, LJ参数以及残基信息.
3. 脚本根据pdb文件中原子的坐标获取APBS网格参数, 并将每帧构型输出到APBS所需的pqr文件, 同时生成APBS输入文件*.apbs. 然后调用APBS计算每帧构型对应的apbs文件, 并计算极性PB, 非极性SA部分的贡献, 再计算MM贡献, 同时进行残基分解, 输出结果.

由于脚本在计算时是分步进行的, 因此你可以先将所有apbs文件都产生出来然后并行计算它们, 最后再统计结果. 当然这只有在计算构型较多的情况下才值得尝试.

文件说明

用户反馈

与GMXPBSA相比

• GMXPBSA只能给出总的结果, 而gmx_mmpbsa可给出MM, PB, SA总的结果和残基分解结果;
• GMXPBSA使用时若体系中存在自定义拓扑或残基, 则需要提供相应的拓扑或力场文件;
• GMXPBSA原有脚本生成并使用的Mm.mdp文件较老, 需要手动修改print_files.dat脚本;
• GMXPBSA出错时需要进入子目录查看错误原因, 重新运行时又要删除子目录, 否则不能运行, 因而调试较为繁琐;
• GMXPBSA需先后运行三个脚本, 而gmx_mmpbsa只需要运行一次;
• GMXPBSA可实现甘氨酸扫描(CAS), gmx_mmpbsa尚未实现.

与g_mmpbsa相比

• g_mmpbsa需要安装, 而gmx_mmpbsa无需安装, 简单修改后可直接使用;
• g_mmpbsa需手动分别运行计算MM, PB, SA的能量, 而gmx_mmpbsa只需要运行一次.
• g_mmpbsa计算速度较慢, 而gmx_mmpbsa脚本计算速度较快, 因可以使用最新版本的APBS程序, 支持并行.
• g_mmpbsa可实现任意残基突变, gmx_mmpbsa尚未实现.
• g_mmpbsa支持其他非极性模型, 如SAV, WCA, gmx_mmpbsa尚未实现.

有待改进

• 结果的误差统计
• GB方法
• 其他非极性模型
• 调用其他PB程序, 如AMBER的PBSA
• 熵的贡献
• 形变能的贡献
• 甘氨酸扫描, 任意残基突变
• APBS mg-para方法处理超大蛋白

测试: g_mmpbsa的1EBZ示例

MM能量

计算方法: 直接对势累加, 考虑溶质(蛋白)的介电常数pdie, 不使用周期性, 不使用截断, 不对构型进行优化.

注意: GMXPBSA未考虑pdie, 且使用周期性, 使用截断, 对构型进行了优化, 但所得结果差距不大.

MM能量分解

计算方法: 直接累加

PB(极性)能量, Polar solvation enrgy

计算方法: 水相能量减去真空中的能量, PB=Esol-Evac

与g_mmpbsa比较, PB能量有差距, 原因在于

• 所用网格不同, g_mmpbsa对每个分子构型使用基于其自身的网格, 而gmx_mmpbsa使用基于所有构型的网格, 所有APBS计算使用相同的网格.

• 有些原子的半径不同, g_mmpbsa根据蛋白原子名称推断半径的方式有些怪异, 与论文所给数据不一致. 对配体, 由于原子名称可以任意取, 所以没有办法直接根据名称推断半径. 下面是g_mmpbsa所给蛋白中原子的半径, 可以看到有很多不一致之处.

    CA      1.7
    CB      1.7
    CD      1.7
    CD1     1.7
    CD1     1.77
    CD2     1.7
    CD2     1.77
    CE      1.7
    CE1     1.77
    CE2     1.77
    CE3     1.77
    CG      1.7
    CG      1.77
    CG1     1.7
    CG2     1.7
    CH2     1.77
    CZ      1.7
    CZ      1.77
    CZ2     1.77
    CZ3     1.77
    H       1.2
    H1      1.2
    H2      1.2
    HA      1.2
    HA1     1.2
    HA2     1.2
    HB      1.2
    HB1     1.2
    HB2     1.2
    HB3     1.2
    HD1     1.
    HD1     1.2
    HD11    1.2
    HD12    1.2
    HD13    1.2
    HD2     1.
    HD2     1.2
    HD21    1.2
    HD22    1.2
    HD23    1.2
    HD3     1.2
    HE      1.2
    HE1     1.
    HE1     1.2
    HE2     1.
    HE2     1.2
    HE21    1.2
    HE22    1.2
    HE3     1.
    HE3     1.2
    HG      1.2
    HG1     1.2
    HG11    1.2
    HG12    1.2
    HG13    1.2
    HG2     1.2
    HG21    1.2
    HG22    1.2
    HG23    1.2
    HH      1.2
    HH11    1.2
    HH12    1.2
    HH2     1.
    HH21    1.2
    HH22    1.2
    HZ      1.
    HZ1     1.2
    HZ2     1.
    HZ2     1.2
    HZ3     1.
    HZ3     1.2
  

• 换用与g_mmpbsa相同的网格与半径后, 所得PB作用能差距很小.
• GMXPBSA所用网格更大, 且使用了来自LJ参数的半径, 并对构型进行了优化, 所得结果也更大.

PB能量分解

计算方法: 残基处于复合物中时的能量减去处于其所属分子(蛋白/配体)中的能量, G(res@com)-G(res@pro)

SA(非极性)能量, APolar solvation energy

计算方法: 每个原子的溶剂可及表面积乘以表面张力, 加上常量

SA能量分解

计算方法: 直接累加

能量分解图示

收敛性测试

对于分解后的各个能量, 只有PB能量对选项的设置比较敏感, 使用默认参数得到的值未必是收敛的. 因此最好进行收敛性测试, 确定每个选项的最佳设置.

计算PB能量时, 有三个选项影响很大, cfac, fadd, df. 下面列出了这三个选项取不同值时, 所得的PB能量.

可见, 网格大小和格点间距对最终结果的影响很大. 对此, 引用一下我咨询专家时的答复.

邮件1

有一条分子动力学模拟得到的轨迹, 是蛋白和配体的, 含有很多帧. 如果计算蛋白和配体之间的PB作用能, 我需要分别计算蛋白+配体, 蛋白, 配体三个体系的PB能量, 然后做差值获得PB相互作用能. 在计算每个体系的时候都要设置网格大小, 而且使用APBS的话, 还需要设置两套网格, 粗糙网格和细密网格.

我的第一个问题是, 只算一帧构型的PB相互作用能的话, 对这三个体系是不是 必须 使用相同的网格大小和格点间距, 还是可以对每个体系使用独立的, 只与其自身坐标有关的网格大小和格点间距? 不考虑计算效率的话, 哪种方法更正确?

如果前一个问题的答案是, 对每帧的三个体系必须使用相同网格大小和格点间距, 那我下面的问题是, 如果我要算很多帧构型的PB相互作用能, 那所有这些帧是不是也要使用相同的网格大小和格点间距, 也就是说我们需要对整条轨迹使用相同的网格大小和格点间距? 是否这样得到的结果更自洽?

我没有很细致的检验过APBS的细节，不过我觉得/感觉，如果 网格取得充分大（粗网格能cover两倍的分子大小以上，对于溶剂化能的计算细网格能cover整个分子那我认为应该是够了 —- 而有的只关心反应活性部位的电场计算，那细网格只要包含比活性部位大一点就行了），格点间距取得足够小（比如0.2 A 应该够了），并且如果你每次算的是 分子的 solvation energy （这样它内部已经解了两次 方程），那么我觉的每次计算（三个体系）应该不需要使用同样的网格，就是说可以各自独立的采用粗网格和细网格，因为这时我们“认为”单个体系APBS计算是准确的。 你可以计算试试看。 这样的话，那么对一帧的计算还是整条轨迹多帧的计算应该都没关系，可独立计算。不过，感觉上对多帧的同一个体系比如复合物，最好采用一样的网格精度，以减少网格精度上到带来互相间的计算差别。

当然，对于一帧中的三个体系的计算，我觉得，如果你每次计算都采用同一套网格来计算，并且网格大小满足 我上面讲的覆盖关系，那么这样当然会略为减少网格选取的不同（坐标系平移，转动，网格尺寸大小）带来数值误差，我估计这个偏差较小（如果像上面说的计算正确的话），但这样就是计算量大（因为你计算单体和复合体时用的是一样大的网格系统）。

这些建议你要设计实验来验证一下

邮件2

选择了一个含200个氨基酸残基的蛋白, 含90个原子的配体的体系进行了一些测试, 结果如下:

(见上表)

大致可以看到, 影响最大的是格点间距, 其他粗细网格的增加值影响小, 如您所言, 只要能覆盖即可. 此外, 使用整体一致的网格还是单独的网格对这个体系影响也很小,.

对太大的蛋白, 我们没法使用很小的网格, 所以我检验了一下, 在网格大小固定, 格点间距不太大的情况下, PB相互作用能大致与格点间距的三次方成正比. 我想问下, 这个关系是否正确, 还是理论上有其他更正确的关系式? 如果有的话, 这种关系式对格点间距适用范围的最大值如何确定? 如果这些问题都能解决的话, 我们就可以做外推得到格点间距为零时的PB相互作用能, 这样就可以避免格点间距设置的影响了.

你前面的观察应该合理的。

后半部分你说 “PB相互作用能大致与格点间距的三次方成正比”，这个应该不对，因为如果这样，格点为零，则PB能为0，当然不是。可能的话也是“PB相互作用能的误差大致与格点间距的三次方成正比”，但这个我的印象里也没有看到有人这么说（当然我关于APBS的文章读的很少）。一般来说是这样，一个软件的计算误差由它的算法决定，对于普通线性有限元来说，它的“解的误差”（L2 norm）与网格尺寸（大致可认为类似于格点间距）成平方关系（二阶），有限差分方法或有限体方法的误差也与具体算法格式（方法）有关，一般的方法对一般光滑问题的解的误差也能达到二阶精度，达到三阶需要更精细一点的方法，但有些问题的有限差分方法的解随着格点间距的减小并不能总是保证一致的精度阶，有时候随着格点间距减小误差甚至会变得更大。 你用的APBS的方法应该是用其中的有限差分方法，其实可能应该说有限体方法更恰当。APBS的计算里面还有奇异电荷分配的问题，还有粗细网格同时用的问题，所以具体的误差阶我没有研究过，不知道它是几阶，大致的三阶也是有可能的。不过就算是这样，你要用这个关系也挺麻烦，因为每次计算你想得到精确值就要多计算若干次得到准确的三阶估计才能用。

一些资料

• MMPBSA
• 计算MMPBSA的新脚本——gmxmmpbsa使用说明