【3.2.2】gromacs教程--蛋白质配体
我们必须下载将要使用的蛋白质结构文件。在本教程中，我们将使用T4溶菌酶L99A/M102Q（PDB代码 3HTB）。转到RCSB网站，并下载晶体结构的PDB文本。
source /opt/fastone/softwares/hpc-kits/hpc-kits.sh
source /data/software/gromacs/gromacs-2020.3/bin/GMXRC
cd /data/user/sam/project/gromacs/test5
1.1 生成拓扑结构
下载完结构后，您可以使用查看程序（例如VMD，Chimera，PyMOL等）将其可视化。查看了分子后，您将要去除结晶水,PO4，和BME。注意，这种方法并非普遍适用（即，结合了活性位点水分子的情况）。就我们的意图而言，我们不需要结晶水或其他物质，它们只是结晶助溶剂。相反，我们将专注于称为“ JZ4”的配体，即2-丙基苯酚。
如果您希望使用.pdb文件的干净版本来检查工作，则可以在此处下载。我们现在面临的问题是，在GROMACS随附的任何力场中，JZ4配体都不是公认的实体，因此，如果您尝试使该文件通过它，则pdb2gmx将给出致命错误。如果在力场的.rtp（残基拓扑）文件中存在构建块的条目，则拓扑只能自动组装。既然不是这种情况，我们将分两步准备系统拓扑：
- 使用pdb2gmx准备蛋白质拓扑
- 使用外部工具准备配体拓扑
由于我们将分别准备这两种拓扑，因此必须将蛋白质和JZ4配体保存到单独的坐标文件中。像这样保存JZ4坐标：
grep JZ4 3HTB_clean.pdb > jz4.pdb
然后只需从3HTB_clean.pdb中删除JZ4行。 在这一点上，准备蛋白质拓扑结构是微不足道的。 我们将在本教程中使用的力场是CHARMM36，可从MacKerell实验室网站获得。 在此处，下载最新的CHARMM36力场压缩包和“ cgenff_charmm2gmx.py”转换脚本，我们将在以后使用。 下载与您安装的Python版本相对应的转换脚本版本（Python 2.x或3.x）。
在工作目录中取消归档力场压缩包：
tar -zxvf charmm36-mar2019.ff.tgz
现在，您的工作目录中应该有一个“ charmm36-mar2019.ff”子目录。 用pdb2gmx编写T4溶菌酶的拓扑：
grep -v 'JZ4' 3HTB_clean.pdb > 3HTB_clean_remove_jz4.pdb
gmx pdb2gmx -f 3HTB_clean_remove_jz4.pdb -o 3HTB_processed.gro
出现提示时，选择默认的水模型。 该结构将由pdb2gmx处理，并且将提示您选择一个力场：
Select the Force Field:
From current directory:
1: CHARMM36 all-atom force field (July 2017)
From '/usr/local/gromacs/share/gromacs/top':
2: AMBER03 protein, nucleic AMBER94 (Duan et al., J. Comp. Chem. 24, 1999-2012, 2003)
3: AMBER94 force field (Cornell et al., JACS 117, 5179-5197, 1995)
4: AMBER96 protein, nucleic AMBER94 (Kollman et al., Acc. Chem. Res. 29, 461-469, 1996)
5: AMBER99 protein, nucleic AMBER94 (Wang et al., J. Comp. Chem. 21, 1049-1074, 2000)
6: AMBER99SB protein, nucleic AMBER94 (Hornak et al., Proteins 65, 712-725, 2006)
7: AMBER99SB-ILDN protein, nucleic AMBER94 (Lindorff-Larsen et al., Proteins 78, 1950-58, 2010)
8: AMBERGS force field (Garcia & Sanbonmatsu, PNAS 99, 2782-2787, 2002)
9: CHARMM27 all-atom force field (CHARM22 plus CMAP for proteins)
10: GROMOS96 43a1 force field
11: GROMOS96 43a2 force field (improved alkane dihedrals)
12: GROMOS96 45a3 force field (Schuler JCC 2001 22 1205)
13: GROMOS96 53a5 force field (JCC 2004 vol 25 pag 1656)
14: GROMOS96 53a6 force field (JCC 2004 vol 25 pag 1656)
15: GROMOS96 54a7 force field (Eur. Biophys. J. (2011), 40,, 843-856, DOI: 10.1007/s00249-011-0700-9)
16: OPLS-AA/L all-atom force field (2001 aminoacid dihedrals)
对于本教程，选择CHARMM36强制字段（选项1），该字段首先在列表中的“来自当前目录”下列出。
1.2 Prepare the Ligand Topology 准备配体拓扑
我们现在必须处理配体。但是，如何为力场无法自动识别的某些物种提供参数呢？正确处理配体是分子模拟中最具挑战性的任务之一。力场作者花了数年的时间来获得自洽的力场，将新物种引入该框架并非易事。任何新物种的力场参数必须以与原始力场一致的方式导出和验证。
对于OPLS，AMBER和CHARMM力场，此推导通常采用各种量子力学计算的形式。这些力场的主要文献描述了所需的过程。对于GROMOS力场，参数化方法尚不清楚，它依赖于凝聚相行为的经验拟合。也就是说，为每种原子类型计算一些初始电荷和Lennard-Jones参数，对其准确性进行评估并进行完善。尽管最终结果非常令人满意，即流体类似于现实世界中的流体，但推导过程可能很费力且令人沮丧。
由于这个原因，非常需要自动化工具。对于每个力场，都有一些方法或软件程序旨在提供与各种力场兼容的参数。并非所有这些都同样准确。有关一些示例，请参考下表：
表格来源： http://www.mdtutorials.com/gmx/complex/02_topology.html
1.2.1 将氢原子添加到JZ4 Add Hydrogen Atoms to JZ4
在本教程中，我们将使用CGenFF服务器生成JZ4拓扑。 单击上表中的链接访问该网站。 注册一个（免费）帐户并激活它。 CGenFF需要一个Sybyl .mol2文件作为其输入，以收集基本原子类型信息和绑定连接。 CHARMM也是一个全原子力场，这意味着所有H原子都被明确表示。 晶体结构通常不分配H坐标，因此必须将其内置。要生成.mol2文件并添加H原子，请使用Avogadro程序。 在Avogadro中打开jz4.pdb，然后从“构建”菜单中选择“添加氢”。 Avogadro会将所有H原子构建到JZ4配体上。 保存一个名为“ jz4.mol2”的.mol2文件（文件->另存为…，然后从下拉菜单中选择Sybyl Mol2）。
在使用jz4.mol2之前，必须对其进行一些更正。 在您喜欢的纯文本编辑器中打开jz4.mol2，您将找到：
@<TRIPOS>MOLECULE
*****
22 22 0 0 0
SMALL
GASTEIGER
@<TRIPOS>ATOM
1 C4 24.2940 -24.1240 -0.0710 C.3 167 JZ4167 -0.0650
2 C7 21.5530 -27.2140 -4.1120 C.ar 167 JZ4167 -0.0613
3 C8 22.0680 -26.7470 -5.3310 C.ar 167 JZ4167 -0.0583
4 C9 22.6710 -25.5120 -5.4480 C.ar 167 JZ4167 -0.0199
5 C10 22.7690 -24.7300 -4.2950 C.ar 167 JZ4167 0.1200
6 C11 21.6930 -26.4590 -2.9540 C.ar 167 JZ4167 -0.0551
7 C12 22.2940 -25.1870 -3.0750 C.ar 167 JZ4167 -0.0060
8 C13 22.4630 -24.4140 -1.8080 C.3 167 JZ4167 -0.0245
9 C14 23.9250 -24.7040 -1.3940 C.3 167 JZ4167 -0.0518
10 OAB 23.4120 -23.5360 -4.3420 O.3 167 JZ4167 -0.5065
11 H 25.3133 -24.3619 0.1509 H 1 UNL1 0.0230
12 H 23.6591 -24.5327 0.6872 H 1 UNL1 0.0230
13 H 24.1744 -23.0611 -0.1016 H 1 UNL1 0.0230
14 H 21.0673 -28.1238 -4.0754 H 1 UNL1 0.0618
15 H 21.9931 -27.3472 -6.1672 H 1 UNL1 0.0619
16 H 23.0361 -25.1783 -6.3537 H 1 UNL1 0.0654
17 H 21.3701 -26.8143 -2.0405 H 1 UNL1 0.0621
18 H 21.7794 -24.7551 -1.0588 H 1 UNL1 0.0314
19 H 22.2659 -23.3694 -1.9301 H 1 UNL1 0.0314
20 H 24.5755 -24.2929 -2.1375 H 1 UNL1 0.0266
21 H 24.0241 -25.7662 -1.3110 H 1 UNL1 0.0266
22 H 23.7394 -23.2120 -5.1580 H 1 UNL1 0.2921
@<TRIPOS>BOND
1 4 3 ar
2 4 5 ar
3 3 2 ar
4 10 5 1
5 5 7 ar
6 2 6 ar
7 7 6 ar
8 7 8 1
9 8 9 1
10 9 1 1
11 1 11 1
12 1 12 1
13 1 13 1
14 2 14 1
15 3 15 1
16 4 16 1
17 6 17 1
18 8 18 1
19 8 19 1
20 9 20 1
21 9 21 1
22 10 22 1
需要进行的第一个更改是在MOLECULE标题中。 将“ *****”替换为“ JZ4”，例如：
@ <TRIPOS>分子
JZ4
接下来，固定残基名称和编号，使它们都相同。 该.mol2文件仅包含一个分子，因此应仅指定一个残基名称和编号。 编辑后，jz4.mol2的更正的ATOM部分应显示为：
@<TRIPOS>ATOM
1 C4 24.2940 -24.1240 -0.0710 C.3 1 JZ4 -0.0650
2 C7 21.5530 -27.2140 -4.1120 C.ar 1 JZ4 -0.0613
3 C8 22.0680 -26.7470 -5.3310 C.ar 1 JZ4 -0.0583
4 C9 22.6710 -25.5120 -5.4480 C.ar 1 JZ4 -0.0199
5 C10 22.7690 -24.7300 -4.2950 C.ar 1 JZ4 0.1200
6 C11 21.6930 -26.4590 -2.9540 C.ar 1 JZ4 -0.0551
7 C12 22.2940 -25.1870 -3.0750 C.ar 1 JZ4 -0.0060
8 C13 22.4630 -24.4140 -1.8080 C.3 1 JZ4 -0.0245
9 C14 23.9250 -24.7040 -1.3940 C.3 1 JZ4 -0.0518
10 OAB 23.4120 -23.5360 -4.3420 O.3 1 JZ4 -0.5065
11 H 25.3133 -24.3619 0.1509 H 1 JZ4 0.0230
12 H 23.6591 -24.5327 0.6872 H 1 JZ4 0.0230
13 H 24.1744 -23.0611 -0.1016 H 1 JZ4 0.0230
14 H 21.0673 -28.1238 -4.0754 H 1 JZ4 0.0618
15 H 21.9931 -27.3472 -6.1672 H 1 JZ4 0.0619
16 H 23.0361 -25.1783 -6.3537 H 1 JZ4 0.0654
17 H 21.3701 -26.8143 -2.0405 H 1 JZ4 0.0621
18 H 21.7794 -24.7551 -1.0588 H 1 JZ4 0.0314
19 H 22.2659 -23.3694 -1.9301 H 1 JZ4 0.0314
20 H 24.5755 -24.2929 -2.1375 H 1 JZ4 0.0266
21 H 24.0241 -25.7662 -1.3110 H 1 JZ4 0.0266
22 H 23.7394 -23.2120 -5.1580 H 1 JZ4 0.2921
最后，请注意@
perl sort_mol2_bonds.pl jz4.mol2 jz4_fix.mol2
在下一步中使用“ jz4_fix.mol2”。
1.2.2 使用CGenFF生成JZ4拓扑
现在可以使用jz4_fix.mol2文件生成拓扑了。访问CGenFF服务器，登录到您的帐户，然后单击页面顶部的“上载分子”。上传jz4_fix.mol2，CGenFF服务器将以CHARMM“流”文件（扩展名为.str）的形式快速返回拓扑。将其内容从Web浏览器保存到名为“ jz4.str”的文件中。您也可以在这里下载此文件的副本。
CHARMM流文件包含所有拓扑信息-原子类型，电荷和绑定连接。它还有一些附加的粘结参数部分，这些参数是通过类比为力场未涵盖的任何内部相互作用生成的。 CGenFF还为每个参数提供罚分，即评估分配的参数的可靠性。低于10的任何东西都被认为可以立即使用。介于10到50之间的值表示必须对拓扑进行一些验证，并且任何大于50的惩罚通常需要手动重新参数化。罚分是CGenFF服务器最重要的功能之一。许多其他服务器生成拓扑，它们是“黑匣子”，用户只需要隐式信任即可。在没有验证的情况下将整个研究项目放在自动程序上非常危险。较差的配体拓扑会导致大量的时间浪费和不可靠的结果。始终验证新参数化物种的拓扑！至少应对照力场中现有的部分检查电荷的大小和分配给配体的原子类型。
检查jz4.str的内容，并查看charge和新的二面体参数的惩罚。它们都非常低，表明该拓扑的质量非常好，可以直接用于我们的仿真。
CHARMM格式的JZ4拓扑很好，但是如果我们尝试在GROMACS中运行仿真，它就没有用。使用我们先前从MacKerell网站下载的cgenff_charmm2gmx.py脚本。脚本名称将包含版本_py2或_py3，但为简单起见，此处省略了此信息，但是您将需要使用为Python版本下载的脚本的实际名称。使用以下命令执行转换：
activate3
python3 cgenff_charmm2gmx_py3_nx2.py JZ4 jz4_fix.mol2 jz4.str charmm36-mar2019.ff
请注意：此脚本需要NetworkX软件包，版本1.11-与NetworkX 2.x不兼容，并且将退出并显示明显错误。
成功转换结束时，请注意以下屏幕输出：
============ DONE ============
Conversion complete.
The molecule topology has been written to jz4.itp
Additional parameters needed by the molecule are written to jz4.prm, which needs to be included in the system .top
============ DONE ============
该配体引入了新的结合参数，这些参数不是现有力场的一部分，并且将这些参数写入到名为“ jz4.prm”的文件中，该文件采用GROMACS .itp文件的格式。 我们将在短期内处理此文件，但重要的是要注意它的存在。 现在将配体拓扑写入“ jz4.itp”，其中包含配体[分子类型]定义。
1.2.3 建立综合体 Build the Complex
在pdb2gmx中，我们有一个名为“ 3HTB_processed.gro”的文件，其中包含蛋白质的已加工，与力场相符的结构。 我们还从cgenff_charmm2gmx.py获得了“ jz4_ini.pdb”，它具有所有必需的H原子，并且与JZ4拓扑中的原子名称匹配。 使用editconf将此.pdb文件转换为.gro格式：
gmx editconf -f jz4_ini.pdb -o jz4.gro
将3HTB_processed.gro复制到要操作的新文件中，例如，将其称为“ complex.gro”，因为在蛋白质中添加JZ4会生成我们的蛋白质-配体复合物。 接下来，只需复制jz4.gro的坐标部分并将其粘贴到complex.gro中，在蛋白质原子的最后一行之下，在框矢量之前，如下所示：
163ASN C 1691 0.621 -0.740 -0.126
163ASN O1 1692 0.624 -0.616 -0.140
163ASN O2 1693 0.683 -0.703 -0.011
5.99500 5.19182 9.66100 0.00000 0.00000 -2.99750 0.00000 0.00000 0.00000
添加以后
163ASN C 1691 0.621 -0.740 -0.126
163ASN O1 1692 0.624 -0.616 -0.140
163ASN O2 1693 0.683 -0.703 -0.011
1JZ4 C4 1 2.429 -2.412 -0.007
1JZ4 C7 2 2.155 -2.721 -0.411
1JZ4 C8 3 2.207 -2.675 -0.533
1JZ4 C9 4 2.267 -2.551 -0.545
1JZ4 C10 5 2.277 -2.473 -0.430
1JZ4 C11 6 2.169 -2.646 -0.295
1JZ4 C12 7 2.229 -2.519 -0.308
1JZ4 C13 8 2.246 -2.441 -0.181
1JZ4 C14 9 2.392 -2.470 -0.139
1JZ4 OAB 10 2.341 -2.354 -0.434
1JZ4 H1 11 2.531 -2.436 0.015
1JZ4 H2 12 2.366 -2.453 0.069
1JZ4 H3 13 2.417 -2.306 -0.010
1JZ4 H4 14 2.107 -2.812 -0.407
1JZ4 H5 15 2.199 -2.735 -0.617
1JZ4 H6 16 2.304 -2.518 -0.635
1JZ4 H7 17 2.137 -2.681 -0.204
1JZ4 H8 18 2.178 -2.476 -0.106
1JZ4 H9 19 2.227 -2.337 -0.193
1JZ4 H10 20 2.458 -2.429 -0.214
1JZ4 H11 21 2.402 -2.577 -0.131
1JZ4 H12 22 2.374 -2.321 -0.516
5.99500 5.19182 9.66100 0.00000 0.00000 -2.99750 0.00000 0.00000 0.00000
由于我们已经在.gro文件中添加了22个以上的原子，因此请增加complex.gro的第二行以反映此更改。 现在，坐标文件中应该有2636个原子。
1.2.4 建立拓扑（修改topol.top）
在系统拓扑中包括JZ4配体的参数非常容易。 在包含位置限制文件后，只需在topol.top中插入一行，显示#include“ jz4.itp”。 位置限制的包含表明“蛋白质”分子类型部分的结束。
; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif
; Include water topology
#include "./charmm36-mar2019.ff/tip3p.itp"
变成
; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif
; Include ligand topology
#include "jz4.itp"
; Include water topology
#include "./charmm36-mar2019.ff/tip3p.itp"
配体引入了新的二面体参数，这些参数由cgenff_charmm2gmx.py脚本写入“ jz4.prm”。 在topol.top的TOP处，插入#include语句以添加以下参数：
; Include forcefield parameters
#include "./charmm36-mar2019.ff/forcefield.itp"
[ moleculetype ]
; Name nrexcl
Protein_chain_A 3
变成：
; Include forcefield parameters
#include "./charmm36-mar2019.ff/forcefield.itp"
; Include ligand parameters
#include "jz4.prm"
[ moleculetype ]
; Name nrexcl
Protein_chain_A 3
#include语句的位置很关键-它必须出现在任何[Moleculartype]条目之前，因为必须在构造任何分子之前定义所有参数。 它也必须出现在父力场的#include语句之后，因为必须引入所有原子类型，然后才能引入利用它们的键合参数。
最后要做的调整是在[分子]指令中。 考虑到complex.gro中存在一个新分子这一事实，我们必须在此处添加它，如下所示：
[ molecules ]
; Compound #mols
Protein_chain_A 1
变成：
[ molecules ]
; Compound #mols
Protein_chain_A 1
JZ4 1
现在，拓扑和坐标文件在系统内容方面是一致的。
1.3 Defining the Unit Cell & Adding Solvent 定义晶胞并添加溶剂
至此，工作流程就像其他任何MD模拟一样。 我们将定义单位晶胞并填充水。
gmx editconf -f complex.gro -o newbox.gro -bt dodecahedron -d 1.0
gmx solvate -cp newbox.gro -cs spc216.gro -p topol.top -o solv.gro
可视化solv.gro时，您可能会想知道editconf为什么不产生所需的十二面体晶胞形状，因为系统看起来像这样：
GROMACS程序始终使用最数值有效的坐标表示形式，将所有内容重新包装为三斜晶胞。 mdrun执行的物理计算可以使用不同的坐标换行等效地进行，因此最有效的方法是首选。 生成.tpr文件后，可以稍后恢复所需的晶胞形状。
1.4 Adding Ions 添加离子
现在，我们有了包含带电蛋白质的溶剂化系统。 pdb2gmx的输出告诉我们，该蛋白质的净电荷为+ 6e（基于其氨基酸组成）。 如果您在pdb2gmx输出中错过了此信息，请查看topol.top中[atoms]指令的最后一行； 它应该（部分）显示为“ qtot 6”。 由于生命不存在净电荷，因此我们必须向系统中添加离子。
使用grompp可以使用任何.mdp文件来汇编.tpr文件。 我使用.mdp文件来运行能量最小化，因为它们需要最少的参数，因此最容易维护。 例如这个。
gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr
现在，我们将.tpr文件传递给genion：
gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral
选择的15
Will try to add 0 NA ions and 6 CL ions.
Select a continuous group of solvent molecules
Group 0 ( System) has 33518 elements
Group 1 ( Protein) has 2614 elements
Group 2 ( Protein-H) has 1301 elements
Group 3 ( C-alpha) has 163 elements
Group 4 ( Backbone) has 489 elements
Group 5 ( MainChain) has 651 elements
Group 6 ( MainChain+Cb) has 803 elements
Group 7 ( MainChain+H) has 813 elements
Group 8 ( SideChain) has 1801 elements
Group 9 ( SideChain-H) has 650 elements
Group 10 ( Prot-Masses) has 2614 elements
Group 11 ( non-Protein) has 30904 elements
Group 12 ( Other) has 22 elements
Group 13 ( JZ4) has 22 elements
Group 14 ( Water) has 30882 elements
Group 15 ( SOL) has 30882 elements
Group 16 ( non-Water) has 2636 elements
Select a group: 15
Selected 15: 'SOL'
用-pname和-nname指定的离子名称在以前的GROMACS版本中是特定于力场的，但在4.5版中已标准化。 指定的原子名称始终是所有大写字母中的元素符号，以及[moleculetype]。 残基名称可能会也可能不会附加charge符号（+/-）。 如果遇到困难，请参考ions.itp以获得适当的命名。
您的topol.top [molecules]指令现在应如下所示：
[ molecules ]
; Compound #mols
Protein_chain_A 1
JZ4 1
SOL 10228
CL 6
1.5 Energy Minimization 能源最小化
现在已经组装好系统，使用grompp使用此输入参数文件创建二进制输入：
gmx grompp -f em.mdp -c solv_ions.gro -p topol.top -o em.tpr
确保在运行genbox和genion时已更新了topol.top文件，否则您将收到很多讨厌的错误消息（“坐标文件中的坐标数与拓扑不匹配”等）。
现在，我们准备调用mdrun来执行EM：
gmx mdrun -v -deffnm em
对我来说，该系统收敛比较快(16min， g01)：
Steepest Descents converged to Fmax < 1000 in 143 steps
Potential Energy = -4.9014547e+05
Maximum force = 8.7411469e+02 on atom 27
Norm of force = 5.6676244e+01
如溶菌酶教程中一样，可以使用能量模块监视势能的各种成分。 我将不在这里说明这些原则。 玩得开心。
现在我们的系统处于最低能耗状态，我们可以开始真正的动力了。
1.6 Equilibration 平衡
平衡我们的蛋白质-配体复合物就像平衡任何其他在水中含有蛋白质的系统一样。 在这种情况下，有一些特殊的注意事项：
- 对配体施加约束
- 温度偶合基团的处理
1.6.1 限制配体 Restraining the Ligand
为了约束配体，我们将需要为其生成位置约束拓扑。 首先，为JZ4创建一个仅包含其非氢原子的索引组：
gmx make_ndx -f jz4.gro -o index_jz4.ndx
...
> 0 & ! a H*
> q
然后，执行genrestr模块并选择此新创建的索引组（在index_jz4.ndx文件中将是组3）：
gmx genrestr -f jz4.gro -n index_jz4.ndx -o posre_jz4.itp -fc 1000 1000 1000
现在，我们需要将此信息包括在拓扑中(topol.top)。 我们可以根据希望使用的条件以几种方式进行此操作。 如果我们只是想在蛋白质也被限制时就限制配体，请在指定的位置在拓扑中添加以下几行：
; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif
; Include ligand topology
#include "jz4.itp"
; Include water topology
#include "./charmm36-mar2019.ff/tip3p.itp"
变成：
; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif
; Include ligand topology
#include "jz4.itp"
; Ligand position restraints
#ifdef POSRES
#include "posre_jz4.itp"
#endif
; Include water topology
#include "./charmm36-mar2019.ff/tip3p.itp"
位置很重要！ 您必须按照指示将对posre_jz4.itp的调用放入拓扑中。 jz4.itp中的参数为我们的配体定义了[moleculetype]指令。 分子类型以水拓扑结构（tip3p.itp）结尾。 将调用放置在posre_jz4.itp的任何其他位置将尝试将位置限制参数应用于错误的分子类型。
如果您希望在平衡过程中获得更多控制，即独立地限制蛋白质和配体，则可以控制将配体位置限制文件包含在其他#ifdef块中，如下所示：
; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif
; Include ligand topology
#include "jz4.itp"
; Include water topology
#include "./charmm36-mar2019.ff/tip3p.itp"
变成
; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif
; Include ligand topology
#include "jz4.itp"
; Ligand position restraints
#ifdef POSRES_LIG
#include "posre_jz4.itp"
#endif
; Include water topology
#include "./charmm36-mar2019.ff/tip3p.itp"
在后一种情况下，要限制蛋白质和配体，我们需要在.mdp文件中指定define = -DPOSRES -DPOSRES_LIG。 您要如何对待系统取决于您。 这些示例仅用于说明GROMACS提供的灵活性。 对于标准的平衡程序，同时限制蛋白质和配体可能就足够了。 您自己的需求可能会有所不同。
1.6.2 温控器 Thermostats
正确控制温度耦合是一个敏感的问题。 将每种分子类型耦合到其自身的恒温基团是一个坏主意。 例如，如果您执行以下操作：
tc-grps = Protein JZ4 SOL CL
您的系统可能会崩溃，因为温度耦合算法的稳定性不足以控制具有几个原子（即JZ4和CL）的基团产生的动能波动。 不要单独耦合系统中的每个物种。
典型的方法是设置tc-grps = Protein Non-Protein 并继续进行。 不幸的是，“非蛋白质”组也包括JZ4。 由于JZ4和蛋白质在物理上紧密相连，因此最好将它们视为单个实体。 也就是说，出于温度偶联的目的，JZ4与蛋白质分组在一起。 同样，我们插入的少量Cl-被认为是溶剂的一部分。 为此，我们需要一个特殊的索引组来合并蛋白质和JZ4。 我们使用make_ndx模块完成此任务，提供完整系统的任何坐标文件。 在这里，我使用em.gro，这是系统的输出（最小化）结构：
gmx make_ndx -f em.gro -o index.ndx
...
0 System : 33506 atoms
1 Protein : 2614 atoms
2 Protein-H : 1301 atoms
3 C-alpha : 163 atoms
4 Backbone : 489 atoms
5 MainChain : 651 atoms
6 MainChain+Cb : 803 atoms
7 MainChain+H : 813 atoms
8 SideChain : 1801 atoms
9 SideChain-H : 650 atoms
10 Prot-Masses : 2614 atoms
11 non-Protein : 30892 atoms
12 Other : 22 atoms
13 JZ4 : 22 atoms
14 CL : 6 atoms
15 Water : 30864 atoms
16 SOL : 30864 atoms
17 non-Water : 2642 atoms
18 Ion : 6 atoms
19 JZ4 : 22 atoms
20 CL : 6 atoms
21 Water_and_ions : 30870 atoms
将“ Protein”和“ JZ4”组与以下内容合并，其中“>”表示make_ndx提示：
> 1 | 13
> q
现在，我们可以设置tc-grps = Protein_JZ4 Water_and_ions来实现所需的“非蛋白质蛋白质”效果。
使用此.mdp文件进行NVT平衡。
gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -n index.ndx -o nvt.tpr
gmx mdrun -deffnm nvt
NVT模拟完成后，请使用以下.mdp文件继续进行NPT：
gmx grompp -f npt.mdp -c nvt.gro -t nvt.cpt -r nvt.gro -p topol.top -n index.ndx -o npt.tpr
gmx mdrun -deffnm npt
1.8 Production MD 生产MD
在完成两个平衡阶段后，系统现在已在所需的温度和压力下达到了良好的平衡。 现在，我们准备释放位置限制并运行生产MD进行数据收集。 这个过程就像我们之前看到的一样，因为我们将把检查点文件（在这种情况下现在包含保留压力耦合信息）用于grompp。 我们将运行10 ns的MD模拟，有关脚本，请参见此处。
gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -n index.ndx -o md_0_10.tpr
gmx mdrun -deffnm md_0_10
1.9 数据分析
1.9.1 重新排列和重新包装坐标 Recentering and Rewrapping Coordinates
就像在周期性边界条件下进行的任何模拟中一样，分子可能会在盒子中“破碎”或来回跳跃。 要更新蛋白质并重新包裹单位细胞内的分子以恢复所需的菱形十二面体形状，请调用trjconv：
gmx trjconv -s md_0_10.tpr -f md_0_10.xtc -o md_0_10_center.xtc -center -pbc mol -ur compact
选择“蛋白质”作为居中，选择“系统”作为输出。 请注意，对于较长的模拟（涉及跨周期性边界的多次跳跃），centering复合物（蛋白质-配体，蛋白质-蛋白质）可能比较困难。 在那些情况下（尤其是在蛋白质-蛋白质复合物中），可能有必要创建一个自定义索引组以用于居中，它对应于复合物中一种蛋白质的活性位点或一种单体的界面残基。
要提取轨迹的第一帧（t = 0 ns），请对最近的轨迹使用trjconv -dump：
gmx trjconv -s md_0_10.tpr -f md_0_10_center.xtc -o start.pdb -dump 0
为了获得更平滑的可视化效果，执行旋转和平移拟合可能是有益的。 执行trjconv，如下所示：
gmx trjconv -s md_0_10.tpr -f md_0_10_center.xtc -o md_0_10_fit.xtc -fit rot+trans
选择“主干”对蛋白质骨架进行最小二乘拟合，选择“系统”进行输出。 请注意，同时进行PBC重新包装和坐标拟合在数学上是不兼容的。 如果您希望执行拟合（这对于可视化很有用，但对于大多数分析例程而言不是必需的），请分别执行这些坐标操作，如此处所示。
1.9.2 分析蛋白质-配体相互作用和配体动力学
本教程无法涵盖您可能希望执行的所有分析方法。 这里将说明一些基本操作。
2-丙基苯酚配体可以与Gln102侧链形成氢键。 GROMACS hbond模块可轻松用于计算任何原子组之间的氢键数量，但在我们的情况下，唯一值为1或0。有关配体如何与Gln102相互作用的更详细的研究，我们 将计算JZ4的羟基与Gln102的羰基O原子之间的距离。 对于要形成的氢键，通常的标准是供体和受体原子之间的距离应≤3.5Å（0.35 nm）。 使用距离模块，使用命令行选择语法（请参见gmx帮助选择以获取示例和更多语法），来计算轨迹上的距离。
gmx distance -s md_0_10.tpr -f md_0_10_center.xtc -select 'resname "JZ4" and name OAB plus resid 102 and name OE1' -oall
平均距离为0.31±0.05 nm，与氢键的形成一致。
通常用于确定氢键存在的第二个标准是供体，氢和受体原子之间形成的夹角。 有不同的惯例来计算角度。 在GROMACS hbond模块中，该角度定义为氢供体受体，该角度应≤30°。 要执行此分析，请首先为供体原子（必须同时包括供体重原子和键合氢）和受体原子创建索引组：
gmx make_ndx -f em.gro -n index.ndx
...
> 13 & a OAB | a H12
(creates group 23)
> 1 & r 102 & a OE1
(creates group 24)
> 23 | 24
> q
使用角度模块分析这三个原子形成的角度：
gmx angle -f md_0_10_center.xtc -n index.ndx -ov angle.xvg
请注意，与大多数其他GROMACS分析模块不同，angle不带-s参数。 角度计算不需要质量或周期性信息，原子名称等，因此无需.tpr或坐标文件即可处理轨迹。 该命令返回的值应约为23±9°。 由于我们构建的索引组的顺序为OAB，H12，OE1，因此该结果可能有些出乎意料，这将对应于供体-氢-受体距离，我们预计该距离接近线性（〜180°）。 检查索引文件的内容，您会发现：
[ JZ4_&_OAB_H12_Protein_&_r_102_&_OE1 ]
1616 2624 2636
make_ndx自动将原子序数从低到高排序，因此计算的结果是受体-供体-氢角，该角与hbond模块所计算的角相同。 因此，该结果与氢键的形成一致，因为它≤30°。 为了获得所需的供体氢受体角度，我们将必须在文本文件中手动编辑索引组以重新排列原子序号（2624 2636 1616）。 使用该索引组重新运行角度计算将得出平均值147±11°。
最后，我们可能对量化模拟过程中配体结合姿势发生了多少变化感兴趣。 要计算仅JZ4的重原子RMSD，请为其创建一个新的索引组：
gmx make_ndx -f em.gro -n index.ndx
...
> 13 & ! a H*
> name 26 JZ4_Heavy
> q
执行rms模块，选择“ Backbone”进行最小二乘拟合，选择“ JZ4_Heavy”进行RMSD计算。 这样，通过拟合去除了蛋白质的整体旋转和翻译，并且报告的RMSD是JZ4位置相对于蛋白质变化了多少，这很好地表明了在模拟过程中结合姿势得以保留的程度。
gmx rms -s em.tpr -f md_0_10_center.xtc -n index.ndx -tu ns -o rmsd_jz4.xvg
计算得出的RMSD应该约为0.1 nm（1Å），表明配体位置仅发生很小的变化。
1.9.2 蛋白质-配体相互作用能 Protein-Ligand Interaction Energy
为了量化JZ4和T4溶菌酶之间相互作用的强度，计算这两个物种之间的非键相互作用能可能是有用的。 GROMACS能够分解任意数量的已定义基团之间的短距离非结合能。重要的是要注意，该数量不是自由能或结合能。实际上，大多数力场的参数设置方式实际上并没有实际意义。 CHARMM的参数设置为专门针对与水的量子力学相互作用能，因此它在本质上与有意义的量保持平衡，因此相互作用能非常有用。
交互能量的计算是通过.mdp文件中的energygrps关键字进行的。尽管是.mdp关键字，但不应将交互能量计算视为正常模拟的一部分。短程能量的分解与在GPU上运行不兼容，并且也不必要地减慢了计算速度。 mdrun模块不需要执行额外的工作即可执行有效的模拟。因此，仅将相互作用能计算为分析的一部分，而不是动力学。从定义了energygrps = Protein JZ4的.mdp文件中创建一个新的.tpr文件，如下所示：
gmx grompp -f ie.mdp -c npt.gro -t npt.cpt -p topol.top -n index.ndx -o ie.tpr
接下来，使用-rerun选项调用mdrun以从现有的仿真轨迹重新计算能量：
gmx mdrun -deffnm即-rerun md_0_10.xtc -nb cpu
注意，使用-deffnm读取ie.tpr并将所有输出文件写入ie。*作为其文件名。 -rerun选项采用您要为其重新计算能量的轨迹的名称，而-nb cpu告诉mdrun仅尝试在CPU硬件上运行，而忽略任何可能可用的GPU。 如上所述，这种类型的计算不能在GPU上执行。 重新运行应该非常快，只需几分钟即可完成。
通过能量模块提取感兴趣的能量项。 我们感兴趣的术语是Coul-SR：Protein-JZ4和LJ-SR：Protein-JZ4。
gmx energy -f ie.edr -o interaction_energy.xvg
平均短程库伦相互作用能为-20.5±7.4 kJ mol-1，短程Lennard-Jones能量为-99.1±7.2 kJ mol-1。 从这些数量的相对大小得出结论可能很诱人，但是即使CHARMM的参数是针对显式的水相互作用能，但相互作用能进一步分解为这些成分的过程也不一定是真实的。 无法通过实验验证这些数量，因此无法知道它们是否有意义。 然而，在这种情况下，总的相互作用能是有用的。 该值（在按照标准公式将两个数量相加的误差传播后）为-119.6±10.3 kJ mol-1。
参考资料
