GROMACS全称为Groningen Machine for Chemical Simulations,是用于研究生物分子体系的分子动力学程序包。它可以用分子动力学、随机动力学或者路径积分方法模拟溶液或晶体中的任意分子，进行分子能量的最小化，分析构象等。它的模拟程序包包含GROMACS力场(蛋白质、核苷酸、糖等)，研究的范围可以包括玻璃和液晶、到聚合物、晶体和生物分子溶液。GROMACS是一个功能强大的分子动力学的模拟软件，其在模拟大量分子系统的牛顿运动方面具有极大的优势。
GROMACS几乎支持所有分子动力学的算法，也有很多优势，如：精度高、速度快；支持GPU加速；开源、免费等。
本教程所使用的GROMACS版本为GROMACS 2020.2，不同版本的使用方法大同小异。在前面的推文Gromacs中基于MM/PBSA方法进行结合自由能计算并进行能量分解中，我们为大家介绍了MM/PBSA方法以及使用李继存老师的脚本进行结合自由能计算和能量分解的步骤。今天我们为大家介绍另外一款非常方便的基于Gromacs模拟后产生的轨迹文件进行结合自由能计算（包括MM-PB(GB)SA方法计算结合自由能、熵计算、能量分解等）的工具--gmx_MMPBSA。
gmx_MMPBSA是基于Amber中的MMPBSA.py进行的结合自由能计算，因此需要事先安装AmberTool20或者更高版本。AmberTool20可由官方网站免费下载，大家可根据官方手册说明自行安装。
1.安装gmx_MMPBSA
首先需要安装gmx_MMPBSA依赖包PyQt5，在终端输入以下命令：
amber.python -m pip install PyQt5
随后输入以下命令安装gmx_MMPBSA:
amber.python -m pip install gmx_MMPBSA
如果已安装gmx_MMPBSA可输入以下命令更新版本：
amber.python -m pip install gmx_MMPBSA --upgrade
安装完毕后，将gmx_MMPBSA(*/amber20/miniconda/bin)写入环境变量：
export PATH="/path_to_amber_install/amber20/miniconda/bin:$PATH"
2.MM-PB(GB)SA方法计算结合自由能
这里我们默认已使用Gromacs完成了完整的分子动力学流程。在此基础上使用gmx_MMPBSA进行结合自由能计算就很方便了。我们只需提供一个计算的参数文件--mmpbsa.in。
如果需要进行GB计算，mmpbsa.in文件的内容可参考如下，复制内容保存为mmpbsa.in文件即可。
Sample input file for GB calculation
This input file is meant to show only that gmx_MMPBSA works. Althought,
we tried to used the input files as recommended in the Amber manual,
some parameters have been changed to perform more expensive calculations
in a reasonable amount of time. Feel free to change the parameters
according to what is better for your system.
&general
startframe=5, endframe=20, verbose=2, interval=1,
protein_forcefield="oldff/leaprc.ff99SB", ligand_forcefield="leaprc.gaff"
/
&gb
igb=5, saltcon=0.15,
/
可根据自己的需求更改文件中的参数。
在“&general”模块中，startframe指定从完整轨迹中提取计算结合自由能的初始轨迹帧数；endframe指定计算的最后轨迹帧数；verbose指定在输出文件中打印多少输出的变量（verbose=2表示如果使用一个轨迹，也输出详细内容）；interval指定计算的轨迹间隔；protein_forcefield和ligand_forcefield设置蛋白和配体的力场，需根据自己动力学时使用的力场确定。
在“&gb”模块中，igb设定Generalized Born方法，推荐5；saltcon设定离子浓度，一般为0.15。
准备好mmpbsa.in文件后，可新建文件夹把需要的文件均放置在该文件夹中，打开终端，输入以下命令（注意文件名称的对应）：
gmx_MMPBSA -O -i mmpbsa.in -cs com.tpr -ci index.ndx -cg 1 13 -ct com_traj.xtc -lm ligand.mol2
-i：指定参数文件；-cs：指定动力学模拟的.tpr文件；-ci：指定.ndx索引文件；-cg：指定在.ndx文件中受体和配体对应的组；-ct：指定轨迹文件；-lm：动力学模拟的配体文件
等待计算完成后，计算结果输出至“FINAL_RESULTS_MMPBSA.dat”文件中，该文件记录了结合自由能计算结果及各个能量分项。
Differences (Complex - Receptor - Ligand):
Energy Component Average Std. Dev. Std. Err. of Mean
-------------------------------------------------------------------------------
BOND 0.0000 0.0000 0.0000
ANGLE -0.0000 0.0000 0.0000
DIHED -0.0002 0.0036 0.0008
VDWAALS -48.7779 2.5404 0.5544
EEL -71.2140 12.1375 2.6486
1-4 VDW -0.0000 0.0000 0.0000
1-4 EEL -0.0000 0.0000 0.0000
EGB 78.7078 8.3950 1.8319
ESURF -7.5909 0.1825 0.0398
DELTA G gas -119.9922 11.6497 2.5422
DELTA G solv 71.1169 8.3927 1.8314
DELTA TOTAL -48.8753 4.5584 0.9947
其中DELTA TOTAL为总的结合自由能，其他我们比较关注的如：VDWAALS为范德华能量（van der Waals energy）；EEL为静电作用能（Electrostatic energy）；EGB为极性溶剂化能（Polar solvation energy）；ESURF为非极性溶剂化能（Non-polar solvation energy）；DELTA G gas为总气相自由能（Total gas phase free energy）；DELTA G solv为总溶剂化自由能（Total solvation free energy）。
如果同时需要计算PB，mmpbsa.in文件的内容可参考如下，复制内容保存为mmpbsa.in文件即可。
Sample input file for GB and PB calculation
This input file is meant to show only that gmx_MMPBSA works. Althought,
we tried to used the input files as recommended in the Amber manual,
some parameters have been changed to perform more expensive calculations
in a reasonable amount of time. Feel free to change the parameters
according to what is better for your system.
&general
startframe=5, endframe=20, verbose=2, interval=1,
protein_forcefield="oldff/leaprc.ff99SB", ligand_forcefield="leaprc.gaff"
/
&gb
igb=5, saltcon=0.15,
/
&PB
istrng=0.15, fillratio=4.0
/
同样的，输入同样的命令执行计算即可。
gmx_MMPBSA -O -i mmpbsa.in -cs com.tpr -ci index.ndx -cg 1 13 -ct com_traj.xtc -lm ligand.mol2
3.熵效应计算
熵效应的准确计算一般有下面几种方法：正则模分析 （normal-mode analysis）、准协波分析 （ quasi-harmonic analysis）以及准高斯方法 （quasi-Gaussian approach）。前面两种方法可能比较适合于生物体系，它们的原理基本相同，不同之处在于准协波分析的原子波动矩阵不来自于正交分析的计算，而来自于分子动力学的样本。考虑到正则模分析已经被成功用于多个生物体系，因此，在一般情况下会采用正则模分析来计算体系熵效应对结合自由能的影响。
--引自：Egtok，知乎,https://www.zhihu.com/question/28399313/answer/2168578717
如我们使用正则模分析的方法计算熵效应，mmpbsa.in文件内容可参考如下，一般计算注意更改计算的起始和结束帧数即可。
Sample input file for entropy calculations (nmode)
This input file is meant to show only that gmx_MMPBSA works. Althought,
we tried to used the input files as recommended in the Amber manual,
some parameters have been changed to perform more expensive calculations
in a reasonable amount of time. Feel free to change the parameters
according to what is better for your system.
&general
#
startframe=5, endframe=21, verbose=2, interval=1,
protein_forcefield="oldff/leaprc.ff99SB", ligand_forcefield="leaprc.gaff"
/
&gb
igb=2, saltcon=0.150,
/
Note that nmode will use only a fraction of the no. of frames selected in
&general variable (21-5/1=16 in this case). This way, nmode will only
process 2 frames (15th and 16th frames) #note also that some parameters
have been changed to perform the calculation faster (maxcyc=5, drms=100).
The typical values for these parameters are (maxcyc=50000, drms=0.0001)
&nmode
nmstartframe=15, nmendframe=16, nminterval=1,
maxcyc=5, drms=100,
/
此外，在“&general”模块中，可添加参数“qh_entropy=1”，使用准协波分析（quasi-harmonic entropy ，QH)来进行熵效应计算。
需要注意的是，使用正则模分析或准协波分析的方法计算熵效应，会耗费较大的计算资源，如果我们对计算准确度要求不高的话，可以使用Interaction Entropy (IE)方法来计算近似值，准确度虽不如正则模分析或准协波分析方法，但计算速度较快。关于该方法的介绍，可参考文献：
DOI: 10.1021/jacs.6b02682
mmpbsa.in文件内容可参考如下：
Sample input file for entropy calculations (IE)
This input file is meant to show only that gmx_MMPBSA works. Althought,
we tried to used the input files as recommended in the Amber manual,
some parameters have been changed to perform more expensive calculations
in a reasonable amount of time. Feel free to change the parameters
according to what is better for your system.
&general
#
startframe=5, endframe=21, verbose=2, interval=1,
protein_forcefield="oldff/leaprc.ff99SB",
#entropy variable control whether to perform a quasi-harmonic entropy (QH)
# approximation or the Interaction Entropy approximation
# (https://pubs.acs.org/doi/abs/10.1021/jacs.6b02682)
entropy=2, entropy_seg=25, temperature=298
/
&gb
igb=2, saltcon=0.150,
/
输入同样的命令计算即可。计算完成后，“_GMXMMPBSA_iteraction_entropy.dat”文件中记录了每帧的熵效应（-TΔS）计算结果。
Calculation for last 13 frames:
Interaction Entropy (-TΔS): 19.8648 +/- 0.0494
Interaction Entropy per-frame:
Frame # | IE value
4000 0.00
4020 4.11
4040 9.27
4060 9.12
4080 8.98
4100 14.09
4120 14.01
4140 13.94
4160 13.87
4180 13.81
4200 13.99
4220 13.94
4240 20.60
4260 20.55
4280 20.51
4300 20.47
4320 20.44
4340 20.40
4360 20.37
4380 20.34
4400 20.31
4420 20.29
4440 20.26
4460 20.23
4480 20.21
4500 20.19
4520 20.17
4540 20.14
4560 20.12
4580 20.10
4600 20.08
4620 20.06
4640 20.05
4660 20.03
4680 20.01
4700 19.99
4720 19.98
4740 19.96
4760 19.95
4780 19.93
4800 19.92
4820 19.90
4840 19.89
4860 19.88
4880 19.86
4900 19.85
4920 19.84
4940 19.82
4960 19.81
4980 19.80
5000 19.79
当然，“FINAL_RESULTS_MMPBSA.dat”文件中记录了结果汇总，可以查看。
ENTROPY RESULTS (INTERACTION ENTROPY):
Entropy Term Average Std. Dev. Std. Err. of Mean
-------------------------------------------------------------------------------
IE 19.8648 0.0494 0.0137
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Differences (Complex - Receptor - Ligand):
Energy Component Average Std. Dev. Std. Err. of Mean
-------------------------------------------------------------------------------
BOND 0.0000 0.0001 0.0000
ANGLE -0.0000 0.0001 0.0000
DIHED -0.0025 0.0056 0.0008
VDWAALS -48.6267 3.2625 0.4568
EEL -71.7993 12.0651 1.6894
1-4 VDW -0.0000 0.0001 0.0000
1-4 EEL -0.0000 0.0000 0.0000
EPB 89.4338 9.4081 1.3174
ENPOLAR -4.8158 0.1152 0.0161
EDISPER 0.0000 0.0000 0.0000
DELTA G gas -120.4285 11.5944 1.6235
DELTA G solv 84.6180 9.3624 1.3110
DELTA TOTAL -35.8105 5.1304 0.7184
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Using Interaction Entropy Approximation:
DELTA G binding = -15.9457 +/- 5.1307
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
4.残基能量分解
如果需要对结合自由能进行能量分解，获得每个残基的自由能贡献，可将“&decomp”参数添加至mmpbsa.in文件中，运行命令即可。
&decomp
idecomp=2, dec_verbose=3,
print_res="within 4"，csv_format=0,
/
idecomp=1 & 2 时输出结果为单个残基的结合自由能，其中=1将1-4非键相互作用能（1-4 EEL & 1-4 VDW）归为内部势能，而=2则将1-4非键相互作用各自归为静电势能和范德华势能上。idecomp=3 & 4 时输出结果为残基对的结合自由能。
dec_verbose可指定输出的结果信息，其中dec_verbose=0输出总体的ΔG；dec_verbose=1则输出侧链、骨架、总体的ΔG；dec_verbose=2输出配体、受体、复合物的能量值以及总体的ΔG；dec_verbose=3则输出配体、受体、复合物的能量值以及侧链、骨架、总体的ΔG。一般情况下，我们关注总体的ΔG，因此，定义dec_verbose=0即可。
print_res可指定计算能量的残基，print_res="within 4"表示输出配体周围4埃范围内氨基酸残基的能量贡献，也可自己定义氨基酸残基，如：
print_res="A/1,3-10,15,100 B/25"
准备好mmpbsa.in文件后，在终端运行同样的命令即可。计算完成后，“FINAL_DECOMP_MMPBSA.dat”文件是分解自由能计算结果的矩阵。
DELTAS:
Total Energy Decomposition:
Residue,Location,Internal,,,van der Waals,,,Electrostatic,,,Polar Solvation,,,Non-Polar Solv.,,,TOTAL,,
,,Avg.,Std. Dev.,Std. Err. of Mean,Avg.,Std. Dev.,Std. Err. of Mean,Avg.,Std. Dev.,Std. Err. of Mean,Avg.,Std. Dev.,Std. Err. of Mean,Avg.,Std. Dev.,Std. Err. of Mean
ASN 37,R ASN 37,0.0,0.0,0.0,-1.4309411764705884,0.7815249399948214,0.10943538086459338,-8.071313725490196,1.7747180001362721,0.24851022703568404,7.57750980392157,1.0792264726717995,0.15112193358381198,0.0,0.0,0.0,-1.924745098039215,1.4974694188914808,0.20968766037147377
ALA 38,R ALA 38,0.0,0.0,0.0,-0.6670392156862747,0.15301761538787756,0.021426751933310288,-0.2385098039215686,0.12379733161763956,0.01733509379199917,0.07815686274509803,0.15482505909409813,0.02167984467579253,0.0,0.0,0.0,-0.8273921568627449,0.20207745345136235,0.02829650334990399
ASP 40,R ASP 40,0.0,0.0,0.0,0.5048431372549019,0.9310411286496514,0.13037183498590668,-14.873078431372548,2.175897307293652,0.3046865664293475,16.748372549019603,3.394088932944574,0.4752675136222166,0.0,0.0,0.0,2.3801372549019617,2.4134204432843998,0.3379464581117196
ALA 41,R ALA 41,0.0,0.0,0.0,-1.184450980392157,0.3139400663249864,0.04396040224533109,0.37013725490196075,0.15912271988820859,0.022281637557582615,-0.4641372549019608,0.14922427122246623,0.02089557750464135,0.0,0.0,0.0,-1.2784509803921567,0.3584043238603364,0.050186643673115386
LYS 44,R LYS 44,0.0,0.0,0.0,-1.4620392156862747,0.37853473748354693,0.05300546540109817,0.07815686274509868,3.073378312200416,0.4303590441257861,0.8505098039215686,3.4635862994626607,0.48499909144498377,0.0,0.0,0.0,-0.5333725490196067,1.0737390773857807,0.15035354455059186
ASP 79,R ASP 79,0.0,0.0,0.0,0.5180980392156861,0.6115706957323156,0.08563702652091153,-10.518588235294118,1.8278457713338905,0.25594960302736625,14.04972549019608,2.3990752227982193,0.3359377254569407,0.0,0.0,0.0,4.049235294117646,1.767683868332279,0.2475252515683281
ILE 82,R ILE 82,0.0,0.0,0.0,-1.6432156862745095,0.371377156937685,0.05200320365227849,-0.1511960784313727,0.1667027435002305,0.023343053167620254,-0.010352941176470587,0.2542339376251115,0.03559987195404629,0.0,0.0,0.0,-1.804764705882353,0.4480776625714585,0.06274342269966361
GLY 83,R GLY 83,0.0,0.0,0.0,-0.4262745098039217,0.431580318947923,0.06043333252812328,-0.5668039215686272,2.1696541788199557,0.303812353582973,0.9104117647058821,0.9514163059946851,0.1332249304903856,0.0,0.0,0.0,-0.08266666666666672,1.464094495860994,0.20501423636891794
MET 84,R MET 84,0.0,0.0,0.0,-1.80943137254902,0.3511398947462806,0.04916942013209083,-0.3033725490196079,0.2769473625066412,0.03878038760421352,0.32299999999999995,0.18925167391183462,0.026500534984770528,0.0,0.0,0.0,-1.7898039215686272,0.43970159308900986,0.06157053837179285
ASP 88,R ASP 88,0.0,0.0,0.0,-0.26523529411764696,0.09961013969870688,0.013948209478734706,1.5999803921568623,1.0961847649147567,0.15349656947250595,-0.9814901960784325,1.171803175143519,0.16408526485542552,0.0,0.0,0.0,0.3532549019607833,0.5732447666363546,0.08027032299941714
ASN 92,R ASN 92,0.0,0.0,0.0,-2.3814901960784307,0.578206536562406,0.08096510976031503,-1.9740000000000009,1.158507288390928,0.16222346831351103,2.50086274509804,0.9613472168562026,0.13461553615994967,0.0,0.0,0.0,-1.8546274509803922,0.8680187709361051,0.1215469397504299
LEU 93,R LEU 93,0.0,0.0,0.0,-1.7162352941176469,0.22421318533890133,0.03139612580065445,0.2904509803921567,0.23354055059785345,0.032702218181510985,0.020333333333333335,0.24585462519958037,0.034426533523314325,0.0,0.0,0.0,-1.4054509803921569,0.39100074287705877,0.054751055309090205
GLY 121,R GLY 121,0.0,0.0,0.0,-0.6724705882352944,0.14501485426837224,0.020306141232022585,0.20601960784313722,0.5256335296179933,0.07360341630214229,0.8367058823529411,0.54270333838554,0.07599366762767854,0.0,0.0,0.0,0.37025490196078437,0.3063805683656345,0.042901860801558
VAL 122,R VAL 122,0.0,0.0,0.0,-1.2274313725490196,0.19551783980606385,0.027377973715261003,-0.37601960784313737,0.45419420093672747,0.06359990938527155,0.14117647058823538,0.20136862168497838,0.0281972470493646,0.0,0.0,0.0,-1.4622745098039218,0.5899417744094633,0.08260837174417164
GLY 123,R GLY 123,0.0,0.0,0.0,-0.36786274509803935,0.13370585503595472,0.018722564592478355,0.7118627450980394,0.210910818399517,0.029533421851089196,-0.5353921568627452,0.14991044103277784,0.020991660496605415,0.0,0.0,0.0,-0.19139215686274522,0.22407544934086104,0.03137683890316349
PHE 124,R PHE 124,0.0,0.0,0.0,-1.9531372549019599,0.39950655935167956,0.05594210784987109,-0.6548039215686273,0.6978348874774749,0.09771642948746352,0.8308823529411767,0.36048538340599207,0.050478050296661185,0.0,0.0,0.0,-1.7770588235294111,0.9089074115740859,0.1272724946652642
TYR 125,R TYR 125,0.0,0.0,0.0,-0.6450784313725489,0.10533752345920999,0.0147502036200765,-0.30015686274509773,0.16802245094489657,0.02352784917277119,0.4708039215686274,0.24950366570515092,0.034937501397890515,0.0,0.0,0.0,-0.4744313725490196,0.2787563596321543,0.03903369786890551
THR 170,R THR 170,0.0,0.0,0.0,-1.2008235294117646,0.39603662466218725,0.055456219806013694,0.24270588235294105,0.45519693191864,0.06374031980763255,0.4794509803921571,0.46218913511319776,0.06471942409531417,0.0,0.0,0.0,-0.47866666666666713,0.5471678118996813,0.07661881896243079
VAL 172,R VAL 172,0.0,0.0,0.0,-0.4652941176470587,0.09206717371138316,0.012891982974079699,-0.06023529411764676,0.07367245693736547,0.010316207419080405,0.08878431372549023,0.06113447527922584,0.008560538818100264,0.0,0.0,0.0,-0.43674509803921513,0.12673019714388395,0.017745777110552405
LIG 211,L LIG 1,-0.002470588235293654,0.00568573357805688,0.0007961619492442365,-24.313254901960782,1.6312588527971321,0.22842192834662028,-35.89976470588234,6.032459377460879,0.8447132723966475,44.93198039215687,4.087689488275192,0.5723910180122398,0.0,0.0,0.0,-15.283509803921566,3.1477585886241064,0.4407743660978453
这里获得的数据较多，我们可以将数据复制到Excel中，然后以逗号为分隔符进行分列，方便作图和数据分析。
数据结果中主要看average列，可具体查看哪个氨基酸的贡献比较大，可能是相互作用的关键氨基酸，之后根据自己需要进行作图即可，如十分方便的柱状图。
转自：“叮当科研”微信公众号
如有侵权，请联系本站删除！
