GROMACS分子动力学的MM/PBSA结合自由能计算教程
GROMACS不像Amber那样自带MM/PBSA.py结合自由能计算程序【1】。2014年，Kumari等人发布第一个用于GROMACS分子动力学轨迹的MM/PBSA计算程序g_mmpbsa【2】。但该程序安装麻烦，对GROMACS版本有特定要求（非常古老的4.5/4.6/5.0/5.1），现已不再更新。另一个严重的瑕疵是它给出的能量分解结果竟然没有范德华力与静电力项。对于习惯Amber MMPBSA.py程序输出结果的人来说，这是不可接受的。同年，C. Paissoni等人发布GMXPBSA 2.0【3】，次年发布了GMXPBSA 2.1【4】。它在精度和便捷度方面都有所改进，但仍然对GROMACS版本有要求，它给的mdp文件版本太老，修改非常麻烦。更致命的是，计算功能多有欠缺：只能给出总结果，不能给出能量项与能量分解结果。因此，对于要做MMPBSA计算的动力学模拟，我基本不用GROMACS。然而，也有不得不使用的时候。以往的做法是通过MDTraj的mdconvert脚本将GROMACS轨迹转为Amber支持的binpos或NetCDF文件，再利用Amber进行计算，但这需要安装第三方软件。我的做法则更简单，只需通过GROMACS的gmx trjconv命令将轨迹转为通用的PDB文件。这样做，原理上是一样，只是生成的PDB文件巨大，要有较大的存储空间。这个思路是非常自然的。但直到去年，才有人基于Amber开发出一款名为gmx_MMPBSA的自动化程序【5-6】。值得一提的是，早在2019年，李继存老师基于GROMACS的相关模块与APBS程序开发了gmx_mmpbsa脚本【7】，与g_mmpbsa和GMXPBSA相比，有不同方面的改进和修正。使用起来也非常方便，只需输入必要文件，脚本内部已经设定好MM/PBSA有关参数。
出于便捷性、稳健性与完备性等方面的考虑，本文使用gmx_MMPBSA程序来开展计算。本文以GROMACS的官方教程T4溶菌酶与JZ4配体小分子复合物的分子动力学模拟【8】为例，讲述蛋白-小分子体系的结合自由能计算流程。
系统环境Linux系统GROMACS 4.x.x/5.x.x/20xx.xAnaconda 3或Miniconda3OpenMPI或MPICH2（可选）硬盘空间1.3 GB以上
软件安装
既然是GROMACS分子动力学模拟，GROMACS当然是安装好且能直接调用的。除此，gmx_MMPBSA程序还要求安装好AmberTools20或21以及Python3。编译安装AmberTools是个耗时费力的工程。这里提供一个简单的安装方法：使用Anaconda（或Miniconda）联网安装。Anaconda（或Miniconda）的安装很简单，这里不再赘述。尚未安装的读者，请自行百度查询。首先，我们创建一个干净的虚拟环境：conda create -n gmx_mmpbsa python=3.9 -y -q并激活环境：conda activate gmx_mmpbsa
然后执行以下命令安装gmx_MMPBSA的最新版本：conda install -c conda-forge gmx_mmpbsa -y
文件准备
在使用gmx_MMPBSA程序进行MM/PBSA结合自由能计算之前，需要准备好以下文件：
文件含义
complex.topGROMACS拓扑文件
md.tprGROMACS输入文件（二进制）
md.xtcGROMACS轨迹文件（二进制）
index.ndxGROMACS索引文件
本文假定这些文件放在同一个目录里。
经过长时间的动力学模拟，体系很容易出现由PBC引起的问题，包括：跨周期漂移、亚单元分离。这时需要先对轨迹文件做适当处理。这里只介绍最简单的处理方式——整体随蛋白居中；特殊体系（比如，膜蛋白）可能需要额外处理。这方面的工作就留给感兴趣的读者，这里不再赘述。我们使用GROMACS的命令来处理轨迹，得到消除了PBC问题的轨迹文件md-nopbc.xtc：echo 1 0 | gmx trjconv -s md.tpr -f md.xtc -o md-nopbc.xtc -n index.ndx -pbc mol -center -ur compact命令中的1和0分别是Protein组和System组，请根据实际情况做出调整。MMPBSA计算还需要准备一个参数文件mmpbsa.in，可以通过下面的命令即时产生：gmx_MMPBSA --create gb pb decomp该文件列出了GB、PB方法做结合自由能计算和能量分解的所有参数，具体含义请见官方网站【9】。大部分情况下默认参数就够用了，只有少量需要修改。这里给一个精简版的mmpbsa.in文件：# General namelist variables
&general
sys_name = ""
startframe = 1
endframe = 10000
interval = 50
/
# Generalized-Born (GB) namelist variables
&gb
igb = 5
saltcon = 0.150
/
# Poisson-Boltzmann (PB) namelist variables
&pb
istrng = 0.15
npbopt = 0
/
# Energy decomposition namelist variables
&decomp
idecomp = 2
dec_verbose = 3
print_res = "within 6"
/说明：
[*]• startframe和endframe：是用于结合自由能计算的动力学轨迹起始帧和终末帧，请根据实际情况修改。
[*]• interval：采样间隔，即每隔多少帧取一帧进行计算。按照上述例子，最终用于计算的帧数为 (10000 - 1 + 1) / 50 = 200帧。一般而言，200-500帧就足够了。
[*]• igb：GB方法的计算模型，igb = 5即采用GB-OBC2模型，其他模型详见官网或Amber手册。
[*]• saltcon：离子浓度，单位：摩尔浓度（M）。
[*]• istrng：离子强度，单位：摩尔浓度（M）。
[*]• npbopt：采用线性（0）或非线性（1）PB方程。对于一般体系，前者就足够了；而后者更推荐用于高度带电的体系，相应地，需要调整eneopt和cutnb参数。
[*]• idecomp：能量分解方案，数值2表示将按残基分解的1-4 EEL项加到EEL项中，而1-4 VDW项则加至VDW势能项。
[*]• dec_verbose：内容输出程度，数值3表示输出全部数据，包括：复合物、受体、配体的能量以及它们的差值，还有总体、侧链以及骨架的贡献值。
[*]• print_res：指定输出哪些残基的能量分解数据，within 6表示输出距离受体、配体6 Å以内的残基。也可以直接指定残基，例如，"A/1,3-10,15,100 B/25"，表示输出复合物中A链的1、3至10、15、100号残基以及B链的25号残基。
[*]能量计算
[*]下面分别介绍串行和并行两种计算方式：串行计算在包含上述文件的目录里执行以下命令：
[*]gmx_MMPBSA -O -i mmpbsa.in -cs md.tpr -cp complex.top -ci index.ndx -cg 1 13 -ct md-nopbc.xtc -o energy.dat -do decomp.csv -nogui && gmx_MMPBSA --clean
命令中的1和13分别是Protein组（受体）和JZ4组（配体），请根据实际情况做出调整。
并行计算如果安装了MPI，可以使用多核并行计算，这样能极大加快速度，减少等待时间。只需在串行计算的命令基础上加上mpirun -np N（N为核数）以及MPI即可：mpirun -np 2 gmx_MMPBSA MPI -O -i mmpbsa.in -cs md.tpr -cp complex.top -ci index.ndx -cg 1 13 -ct md-nopbc.xtc -o energy.dat -do decomp.dat -nogui && gmx_MMPBSA --clean输出结果成功结束后，总能量输出在energy.dat文件中，使用文本编辑器打开即可见到下图所示内容。其中ΔTOTAL是我们需要的结合自由能，其他以Δ开头的是各个能量项。能量单位是kcal/mol，温度是 298.15 K。SD和SEM分别为标准偏差和平均值标准误差，SD(Prop.)和SEM(Prop.)是由误差传递公式计算得到的值。能量分解数据输出到decomp.csv文件中，可用Excel打开，一般只需看DELTAS下面的数据。
总结
本教程简单介绍了公开报道的几种针对GROMACS动力学轨迹的MM/PBSA结合自由能计算程序脚本，并着重介绍了gmx_MMPBSA程序的安装与使用流程。因系统环境与网络环境千差万别，有的读者可能在安装或使用过程中出现意外情况。因此，我们强烈建议读者使用由计算专家精心打造的殷赋云平台，有浏览器即可访问使用，一键提交计算，更有强大的分析功能。我们的口号是：解放大脑，专注思考，君子善假于物，科研快人一步！
页:
[1]