基于gmx_MMPBSA程序的结合自由能计算-MMGBSA方法 - 计算化学公社

结合自由能，在模拟计算中，用以量化受体-配体之间的结合强度的一种指标。尤其是对于药物设计，蛋白-蛋白互作，蛋白-核酸，核酸-小分子等研究领域，均会涉猎分子动力学模拟的工作。在众多工作中，分子动力学数据后处理中，直接用以量化研究对象之间的识别强度和结合能力的指标较少，所以结合自由能在研究受体-配体间结合强度的尤为重要。但对于初学者，众多动力学程序及计算结合自由能相关程序，难以入手，特别是没有耐心看用户手册的新手们（毕竟大部分都是英文，畏难情绪更加大

）。于是，作为一个踩过坑的新手我来说，再次分享一下有关结合自由能计算的经验。

由于本人刚入坑时，也是一头雾水，不经意间逛到强大的计算化学公社，从此几乎将其视为圣经。所以根据流行程度，选择了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版本，使用其计算过较多体系：蛋白-蛋白，蛋白-小分子，蛋白-核酸，核酸-小分子，核酸-核酸，核酸-多肽，蛋白-多肽等体系。也常与使用此程序的一些友友们交流经常。本人对此程序稍有了解和使用经验，突发奇想来论坛写一篇帖子，互相交流和分享经验。

安装方法，我这里不详细介绍，毕竟安装方法，官网介绍的非常详细（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 安装的时候，非常顺畅；还有一种情况就是不要用校园网，毕竟有些校园网拒绝访问很多网址。

本人先声明一点，我分享的经验不一定全面，所使用的代码和方法，以及命令参数，均是前言中提到的体系中，成功跑通过，如果你按照我的方法，没有跑通，表明缺失了相应的文件，或者中间有些步骤没有正确处理。我的处理方式和方法不一定最优，但可以畅通无阻。一般体系分别包括“蛋白-小分子”，“蛋白-蛋白”，“蛋白-核酸”，“核酸-小分子”为什么不提多肽，其实多肽可以分为两种情况，一种是短肽，还有一种是长肽，短肽你可以把其视为“小分子”处理，这样就跟“蛋白-小分子”体系类似；另一种是长肽，你可以将其为一种较短的蛋白质处理，这样就跟“蛋白-蛋白”体系类似。为什么要提到“核酸-小分子”这种小众的体系，其实这并不小众，可能大家不怎么了解，现在在检测领域，核酸又称适配体，现在更流行的叫法，功能核酸，在抗生素残留和肿瘤标志物检测，会做核酸-小分子的动力学，并且适配体几乎都是单链的DNA和RNA。做这种体系的模拟比较少。还有一个原因是，本人就是做适配体的，所以将其作为一个例子分享一下，因为适配体-小分子的模拟比较小众，对很多刚入门新手，几乎找不到很好的方法和策略。于是借此分享一点个人经验。

如工作流程图所示(https://valdes-tresanco-ms.github.io/gmx_MMPBSA/dev/howworks/),需要对提供的complex，recepotor，ligand通过parmed和tleap进行转换为AMBER程序所能识别的文件；注意：大多因为结构文件有问题报错，均会在这个环节出现，如小分子结构扭曲，键连关系不正确，下载小分钟结构时为2d导致的结构等问题均会在这里报错。另外一类常见报错就是index索引超出范围，这里一般均是你提供给程序的结构文件和你的轨迹文件中原子数不一致导致的。

gmx trjconv -s md_0_1.tpr -f md_0_1.xtc -o md1.xtc -pbc nojump -n index.ndx

复制代码

gmx trjconv -f md1.xtc -s md_0_1.tpr -fit rot+trans -o fit.xtc -n index.ndx

复制代码

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

复制代码

conda activate gmxMMPBSA

复制代码

[INFO ] Building Normal Complex Amber topology...
File "/home/wsy/miniconda3/envs/gmxMMPBSA/bin/gmx_MMPBSA", line 7, in <module>
sys.exit(gmxmmpbsa())
File "/home/wsy/miniconda3/envs/gmxMMPBSA/lib/python3.10/site-packages/GMXMMPBSA/app.py", line 98, in gmxmmpbsa
app.make_prmtops()
File "/home/wsy/miniconda3/envs/gmxMMPBSA/lib/python3.10/site-packages/GMXMMPBSA/main.py", line 681, in make_prmtops
self.FILES.mutant_receptor_prmtop, self.FILES.mutant_ligand_prmtop) = maketop.buildTopology()
File "/home/wsy/miniconda3/envs/gmxMMPBSA/lib/python3.10/site-packages/GMXMMPBSA/make_top.py", line 125, in buildTopology
tops = self.gmxtop2prmtop()
File "/home/wsy/miniconda3/envs/gmxMMPBSA/lib/python3.10/site-packages/GMXMMPBSA/make_top.py", line 566, in gmxtop2prmtop
com_top = self.cleantop(self.FILES.complex_top, self.indexes['COM']['COM'])
File "/home/wsy/miniconda3/envs/gmxMMPBSA/lib/python3.10/site-packages/GMXMMPBSA/make_top.py", line 853, in cleantop
rtemp_top = parmed.gromacs.GromacsTopologyFile(temp_top.name)
File "/home/wsy/miniconda3/envs/gmxMMPBSA/lib/python3.10/site-packages/parmed/gromacs/gromacstop.py", line 346, in __init__
self.read(fname, defines, parametrize)
File "/home/wsy/miniconda3/envs/gmxMMPBSA/lib/python3.10/site-packages/parmed/gromacs/gromacstop.py", line 383, in read
for line in f:
File "/home/wsy/miniconda3/envs/gmxMMPBSA/lib/python3.10/site-packages/parmed/gromacs/_gromacsfile.py", line 42, in __iter__
for line in self._handle:
File "/home/wsy/miniconda3/envs/gmxMMPBSA/lib/python3.10/site-packages/parmed/gromacs/_cpp.py", line 174, in __iter__
self._ppcmdmap[cmd](self, args)
File "/home/wsy/miniconda3/envs/gmxMMPBSA/lib/python3.10/site-packages/parmed/gromacs/_cpp.py", line 29, in wrapper
return func(self, args)
File "/home/wsy/miniconda3/envs/gmxMMPBSA/lib/python3.10/site-packages/parmed/gromacs/_cpp.py", line 283, in _pp_include
raise PreProcessorError(f'Could not find {includefile}')
PreProcessorError: Could not find amber14sb_parmbsc1.ff/forcefield.itp
Error occurred on rank 0.
Exiting. All files have been retained.
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD
with errorcode 1.
NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.

复制代码