本帖最后由 jackshoulder 于 2025-5-15 20:18 编辑
基于gmx_MMPBSA程序的结合自由能计算-MMGBSA方法
1 前言
结合自由能，在模拟计算中，用以量化受体-配体之间的结合强度的一种指标。尤其是对于药物设计，蛋白-蛋白互作，蛋白-核酸，核酸-小分子等研究领域，均会涉猎分子动力学模拟的工作。在众多工作中，分子动力学数据后处理中，直接用以量化研究对象之间的识别强度和结合能力的指标较少，所以结合自由能在研究受体-配体间结合强度的尤为重要。但对于初学者，众多动力学程序及计算结合自由能相关程序，难以入手，特别是没有耐心看用户手册的新手们（毕竟大部分都是英文，畏难情绪更加大 ）。于是，作为一个踩过坑的新手我来说，再次分享一下有关结合自由能计算的经验。
由于本人刚入坑时，也是一头雾水，不经意间逛到强大的计算化学公社，从此几乎将其视为圣经。所以根据流行程度，选择了Gromacs程序来研究分子动力学。为什么选择，可以参考sob老师的帖子，《2021年计算化学公社论坛“你最常用的计算化学程序和DFT泛函”投票结果统计》（http://bbs.keinsci.com/thread-23482-1-1.html）和《2024年计算化学公社举办的计算化学程序和DFT泛函的流行程度投票结果》（http://bbs.keinsci.com/thread-45272-1-1.html）毕竟这是真适合新手入门的程序（对于我而言）。于是做完分子动力学之后，为了我课题研究的完整性，我要通过计算的方法来量化受体-配体间的结合强度，因此又要挑选相关程序来进行计算。起初，本人也是使用了李老师的脚本来进行结合自由能计算，也许是本人能力有限，每次算出来的值要么为0，要么为正值，这明显的不对（因为动力学轨迹来看，我的蛋白-核酸结合良好），于是就不停查文献找其他程序，后里找到了gmx_MMPBSA这个程序。专门为Gromacs开发的用以计算结合自由能程序。
gmx_MMPBSA，该文献发表于2021年，各位可以进入其官方网址，进行详细了解（https://valdes-tresanco-ms.github.io/gmx_MMPBSA/dev/）。对其相关推荐和讨论，可以参考sob老师的《GROMACS做MM/PBSA结合自由能计算最理想的工具：gmx_MMPBSA》（http://bbs.keinsci.com/thread-23634-1-1.html）。正如其官网介绍所说，“gmx_MMPBSA 是基于 AMBER 的 MMPBSA.py 的新工具,旨在使用 GROMACS 文件执行自由能量计算”。因此，在这里，大佬们不要抬杠，毕竟 AMBER我不会，会AMBER的可以直接用MMPBSA.py计算。从我开始使用这个程序以来，已经有一年之久，从其1.6.2版本到现在的1.6.4版本，使用其计算过较多体系：蛋白-蛋白，蛋白-小分子，蛋白-核酸，核酸-小分子，核酸-核酸，核酸-多肽，蛋白-多肽等体系。也常与使用此程序的一些友友们交流经常。本人对此程序稍有了解和使用经验，突发奇想来论坛写一篇帖子，互相交流和分享经验。
![]()
2 安装 本人使用的是ubuntu系统，从最开始的20，到22再到现在的24，均自己安装配置部署，也可以交流分享ubuntu系统部署经验。 安装方法，我这里不详细介绍，毕竟安装方法，官网介绍的非常详细（https://valdes-tresanco-ms.github.io/gmx_MMPBSA/dev/installation/） 。我是1.6.2版本入坑的，最开始我纯属小白，按照官网的安装，包括env.file 始终会卡在git和pip某些程序安装不了，于是找人安装（100元子，大概23年12月左右）。反正也是用着把毕业论文做完了。但由于后续维护，以及后面我继续在组里读博，程序安装始终都要解决。又捣鼓了很久，自己在装1.6.3版本的时候，装成功了。总结经验：有些安装包，可能是在当前的conda源没有，也有可能是网络问题，导致有些环境配置过程中断; 因为，那个时候流行用chatgpt，都知道的，要**上网嘛，于是我在组里的机器（双系统，我都装了个工具），使用env.file 安装的时候，非常顺畅；还有一种情况就是不要用校园网，毕竟有些校园网拒绝访问很多网址。 现在安装相对简单一行命令就可以搞定，其实就是官网给的方法 。只要你网络顺畅，conda源update一下，基本上都能安装成功。- conda env create --file env.yml
由于在使用env.file时，里边有gromacs的安装，如果你的系统中安装了gromacs程序，你可以将env.file中“- gromacs<=2024.3”所在行删除掉。后续在使用的时候，程序会自动读取你当前系统中，bashrc文件中写入的gromacs版本。 3 使用方法 gmx_MMPBSA 和 gmx_MMPBSA_ana两部分主要程序 gmx_MMPBSA 负责计算, gmx_MMPBSA_ana负责绘图 本人先声明一点，我分享的经验不一定全面，所使用的代码和方法，以及命令参数，均是前言中提到的体系中，成功跑通过，如果你按照我的方法，没有跑通，表明缺失了相应的文件，或者中间有些步骤没有正确处理。我的处理方式和方法不一定最优，但可以畅通无阻。一般体系分别包括“蛋白-小分子”，“蛋白-蛋白”，“蛋白-核酸”，“核酸-小分子”为什么不提多肽，其实多肽可以分为两种情况，一种是短肽，还有一种是长肽，短肽你可以把其视为“小分子”处理，这样就跟“蛋白-小分子”体系类似；另一种是长肽，你可以将其为一种较短的蛋白质处理，这样就跟“蛋白-蛋白”体系类似。为什么要提到“核酸-小分子”这种小众的体系，其实这并不小众，可能大家不怎么了解，现在在检测领域，核酸又称适配体，现在更流行的叫法，功能核酸，在抗生素残留和肿瘤标志物检测，会做核酸-小分子的动力学，并且适配体几乎都是单链的DNA和RNA。做这种体系的模拟比较少。还有一个原因是，本人就是做适配体的，所以将其作为一个例子分享一下，因为适配体-小分子的模拟比较小众，对很多刚入门新手，几乎找不到很好的方法和策略。于是借此分享一点个人经验。 后文提到的方法中，均是使用MMGBSA方法进行的计算，MMPBSA和RISM后续再专门写帖子介绍 强烈建议大家，有时间，多去看看官网给的使用手册以及例子。
(1) 准备工作
(2) 输入文件
参数文件：mmpbsa.in 结构文件：tpr（是你正式md时的tpr，一般会命名为md.tpr或者md_0_1.tpr) 索引文件 index.ndx 轨迹文件: xtc或者trr top文件(可选,但被要求) 参考结构文件 pdb (可选) 加粗均是必须提供的。
这里我提供一个基于MMGBSA方法的mmpbsa.in文件，大家可以直接用到自己的计算中，稍作修改即可，主要是更改力场，温度，离子浓度即可。通用性参数文件；
mmpbsa.in
(5.29 KB, 下载次数 Times of downloads: 94)
对相关参数做一定的解释，注：我提供的mmpbsa.in文件中，在接下来没有进行解释的，都是按照官网推荐默认值。
sys_name = "SAs16-1 with SMP" # 用户自定义；
startframe = 1 # 起始帧，建议轨迹中的第一帧，但也可以自己定义从自己轨迹的第多少帧开始；
endframe = 9999999 # 结束帧可以是自己轨迹的最后一帧所对应的帧号，也可以写9999999，这样的目的是使轨迹中的帧尽可能算完，因为你的轨迹也不可能超过9999999帧
forcefields = "leaprc.protein.ff14SB/leaprc.gaff2" # 最关键的一点来了，力场参数，“/”左边是受体力场，右边是配体力场。（这里“/”也可以是“，”两种方法，程序都是能识别的。力场的写法至关重要，见官网对力场的写法要求“https://valdes-tresanco-ms.github.io/gmx_MMPBSA/dev/input_file/#general-namelist-variables”中“forcefields”部分的介绍。这里力场尽可能与你动力学轨迹中所用力场一致。
ions_parameters = 1 # 这个参数跟“forcefields ”字段一样。
PBRadii = 4 # 这个参数，我改过，因为我是用MMGBSA中的“igb=8”所以这个地方就必须写4 ，具体这个值写多少，必须根据你用的哪种“igb”方法；具体见官网“https://valdes-tresanco-ms.github.io/gmx_MMPBSA/dev/input_file/#general-namelist-variables”中“PBRadii”字段。
temperature = 298.15 # 温度，不用过多赘述，这个跟你跑动力学时，保持一致即可。
gmx_path = "" # 这个地方，我常会看见有网友去进行定义，其实这里不需要，如果你是使用的ubuntu系统，只要你在bashrc中添加了gromacs的环境变量，那么这个地方空着不写，程序会自动去检查的系统gromacs版本并使用的。
igb = 8 # 这个地方，根据大家的选择，官网默认用的5 如果用的5 那么记得把PBRadii 改成3，这个地方我习惯性使用8，具体哪个好哪个差，我判别不了。
intdiel = 5 #内部介电常数，我一般写的5，因为我平时算核酸和小分子体系比较都，因为核酸的带电性质。
extdiel = 78.5 # 外部介电常数，一般就默认，78.5
saltcon = 0.15 # 盐浓度，这个就跟你动力学模拟时一致吧，基本上要么为0，要么为0.15，因为现在动力学大部分基本都是按照0.15生理条件在设置。
idecomp = 2 # 残基能量分解，这里我习惯用2 将每一个残基都进行分解到1-4 EEL added to EEL and 1-4 VDW added to VDW potential terms，具体看个，只是分解的不同项归到哪个类别里边。“https://valdes-tresanco-ms.github.io/gmx_MMPBSA/dev/input_file/#decomp-namelist-variables”
dec_verbose = 3 # 这里也是看个人，因为我有时候会去看Complex, Receptor, Ligand,的各项值的分解，以及主链和侧链对能量的贡献，因此我会选择3。
print_res = "within 1000" # 这里，很多伙伴会写6或者5，有或者更小，我一般是选择尽可能的大，这样程序会自动去判断你受体与配体之间多远的距离的残基，去做能量分解，写的越大，基本上会把你体系中，每一个残基都进行考虑，这个地方，其实，建议考虑所有残基，通常，无限远的残基，结合自由能分解出来，也是0，这里不会增加计算的时间，这里只是做统计。因为你写的小了，你也不能完全保证，在整个轨迹中，哪些残基之间是保持在这些范围内的，因为这里在分解的时候，是根据你提供的初始结构去判断的，但你也不能保证初始结构按照within 6或者更短判断出来的残基在动力学过程中，是否有超出6或者其他值，也有可能在轨迹后期，有其它残基距离在6之内，但你在提供做残基分解的时候，最开始没有考虑到这一残基对，这样就不会分解到此残基上。索性，这里就写很大，尽肯能考虑所有残基。
|
(3)具体实现
假设你已经跑完某个体系的动力学模拟。备注：以下动力学代码，来自sob老师培训班中，学习后，总结的。
1）首先对md结束的轨迹进行周期性矫正
- gmx trjconv -s md_0_1.tpr -f md_0_1.xtc -o md1.xtc -pbc nojump -n index.ndx
2）做平动与转动的消除
- gmx trjconv -f md1.xtc -s md_0_1.tpr -fit rot+trans -o fit.xtc -n index.ndx
3） 新建一个文件夹“energy”（名字自定义），将相应的mmpbsa.in文件拷贝进去；并修改其中的力场，温度，离子浓度，内部介电常数
注意1：“蛋白-蛋白”体系，你需要将mmpbsa.in文件中的forcefields = "leaprc.DNA.bsc1,leaprc.gaff2"的力场改为自己md运行时所用的力场，intdiel = 1，内部介电常数改为1。由于我跑蛋白体系，力场一般都会用amber14sb，所以这里提供一个“蛋白-蛋白”体系，我常用的参数文件：
mmpbsa.in
(5.31 KB, 下载次数 Times of downloads: 53)
注意2：“蛋白-小分子”体系，直接提供参数文件了，改动主要是力场的写法。
mmpbsa.in
(5.3 KB, 下载次数 Times of downloads: 93)
注意3：“蛋白-核酸”体系，核酸我一般会用bsc1那个力场，由于我们经常做适配体（含有单链核酸的体系，我曾试过解析的晶体结构，bsc1，OL15，以及OL3等，柑橘bsc1在多个适配体晶体结构，模拟效果较佳）
mmpbsa.in
(5.3 KB, 下载次数 Times of downloads: 12)
注意4：“核酸-小分子”体系，内部介电常数改成5。
mmpbsa.in
(5.29 KB, 下载次数 Times of downloads: 9)
4） gmx_MMPBSA运行
若你用env.file成功且正确安装，那么你应该也正确安装了MPI并行化程序，这样就可以直接用下边的命令：
- mpirun --allow-run-as-root --oversubscribe -np 64 gmx_MMPBSA MPI -O -i mmpbsa.in -cs ../md_0_1.tpr -ct ../fit.xtc -ci ../index.ndx -cg 1 12 -cp ../topol.top -o FINAL_RESULTS_MMPBSA.dat -eo FINAL_RESULTS_MMPBSA.csv
--oversubscribe： 这个命令，没记错应该是我当时报超线程那个错误，加了这个就不报错了。
-np 64 ：这里是系统中的CPU核心数，这里的数字一定要小于你轨迹中的帧数，倘若-np后边的值大于你轨迹中的帧数也会报错，如果报这种错，那就把-np降一降。
-cg 1 12 ：这里是index.ndx中，你受体，配体所在的组号。这里也是与forcefields = "leaprc.DNA.bsc1,leaprc.gaff2"字段对应的，先写受体的，再写配体的。
-nogui ： 如果你在命令上加上这个参数，代表你只进行计算，不自动调用gmx_MMPBSA_ana做图形化的数据可视化分析。 后续可以手动可视化分析，调用gmx_MMPBSA_ana。直接运行就行了
注意：由于我是把gmx_MMPBSA写到了bashrc中，所以就可以直接调用，如果你们严格按照env.file安装的，那你应该要先激活环境
所有的运行就到这里结束了！
总结：此程序对于gromaca用户真的很有好，大家总是容易报错，主要原因还是没有好好的看官网的介绍，其实此程序官网写的是真的非常好，该做的提醒都写到了，只是大家没有仔细认真的去看完，导致处处碰壁报错。
目前只是基于此程序做了MMGBSA模型的结合自由能计算，正如题目所写，后期空闲时，再来总结MMPBSA，以及IE，3D-RISM，还有ALA突变给总结出来。至于本贴提到的一些参数，也仅仅是一个初探，具体设置为多少，没有一个特好的标准，如比：内部介电常数，在不同体系，可能就不太一样。
结语：祝大家工作顺利，科研畅通。祝计算化学公社千千万万同仁，一帆风顺！
| 
发表于 Post on 2025-5-14 15:26:47
收藏 Add to favorites
楼主 Author