Gromacs模拟教程：配体-双链蛋白质复合物体系准备
Gromacs模拟教程：配体-双链蛋白质复合物体系准备
Gromacs是一款功能强大的分子动力学模拟软件，广泛应用于生物大分子体系的模拟研究。本文将详细介绍如何使用Gromacs准备一个蛋白质-配体复合物体系，包括蛋白质和配体的预处理、拓扑文件生成、溶剂化和离子添加等关键步骤。
1. 蛋白质的准备
首先需要从RCSB蛋白质数据库下载所需的蛋白晶体结构（例如3htb）。使用文本编辑器（如Notepad）或分子可视化软件去除结构中的非蛋白分子或离子。
这里采用的是一个经过分子对接后的蛋白质PDB文件和配体小分子的PDB文件。配体是2-丙基苯酚（JZ4），这是一个非肽类分子，同时存在于3htb蛋白中。可以使用Notepad将其提取为一个单独的PDB文件。如果是一个短肽类分子，可以利用Gromacs自带的蛋白力场为其生成力场参数。
接下来，使用以下命令对蛋白质PDB文件进行处理，生成top文件和gro坐标文件。在弹出的力场选项中选择第一个CHARMM36力场：
gmx pdb2gmx -f 3HTB_clean.pdb -o 3HTB_processed.gro -ter
选择TIP3P模型作为水分子，选择NH3+和COO-作为蛋白链的末端原子。
2. 配体分子的转换
使用Discovery Studio为配体分子加氢并添加力场，然后输出为mol2格式文件。需要下载相关的CHARMM36力场文件（charmm36-jul2022.ff.tgz）和一个Python脚本（cenff_charmm2gmx.py）。下载链接为：http://mackerell.umaryland.edu/charmm_ff.shtml#gromacs
将下载的tgz文件解压到存放蛋白和配体PDB文件的工作目录下，成为一个子目录。
3. 为配体mol2文件生成拓扑文件
在CGenFF服务器注册账号，然后使用该网站完成拓扑数据的生成。如果无法注册，可以使用Sob老师提供的Sobtop程序生成Gromacs兼容的拓扑文件。下载界面的下方有11个例子参考，可以根据例子2完成对小分子拓扑文件的生成。
操作步骤如下：
## 打开sobtop.exe
sm.mol2 # 输入小分子的mol2路径
2// # 生成gro文件
enter # 回车键在当前路径下生成gro
1// # 生成top文件
4// # 使用sobtop自带的力场数据，缺少的就靠猜
enter # 回车键在当前路径下生成top
enter # 回车键在当前路径下生成itp
4. 建立复合物的gro文件
用protein.gro拷贝成一个新文件，修改名为complex.gro。将ligand.gro里的原子坐标放到complex里，再在盒子数据上方空一行：
这里推荐使用Notepad3进行编辑，Win10自带的笔记本确实不行，不能看清楚数据之间的空格情况，容易报错。
这里不用留空行，留了会报错。
主要的问题应该还是空格没有处理好。
蛋白有4892个原子，配体是33个原子，所以在complex.gro的第2行（表示原子数）修改为4825。配体里的原子数从1又重新排序了，这个应该不伤大雅。
5. 建立复合物的拓扑文件
同目录下有个topol，是一个包含蛋白质拓扑信息的文件（在使用pdb2gmx时生成）。这个蛋白是一个双链蛋白，因为活性口袋在两个蛋白双螺旋之间。同目录下还有两个分开的独立链的top文件。
在topol里，用#include语句导入蛋白的top信息，如果蛋白是单链的话，数据会直接导入到topol中。根据前面导入力场，导入蛋白链拓扑的命令的格式，插入导入配体拓扑的命令：
6. 避开报错，关于[atomtypes]
利用Sob老师生成的小分子itp文件中，有[atomtypes]这个部分，如果不处理，会报错。确实在官方manual里的top文件模板中，不会出现[atomtypes]这一项。
在manual里，top文件的格式如下，是没有[atomtypes]这个部分的。
;
; Example topology file
;
[ defaults ]
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
1 1 no 1.0 1.0
; The force field files to be included
#include "rt41c5.itp"
[ moleculetype ]
; name nrexcl
Urea 3
[ atoms ]
; nr type resnr residu atom cgnr charge
1 C 1 UREA C1 1 0.683
2 O 1 UREA O2 1 -0.683
3 NT 1 UREA N3 2 -0.622
[ bonds ]
; ai aj funct c0 c1
3 4 1 1.000000e-01 3.744680e+05
3 5 1 1.000000e-01 3.744680e+05
6 7 1 1.000000e-01 3.744680e+05
[ pairs ]
; ai aj funct c0 c1
2 4 1 0.000000e+00 0.000000e+00
2 5 1 0.000000e+00 0.000000e+00
2 7 1 0.000000e+00 0.000000e+00
[ angles ]
; ai aj ak funct c0 c1
1 3 4 1 1.200000e+02 2.928800e+02
1 3 5 1 1.200000e+02 2.928800e+02
4 3 5 1 1.200000e+02 3.347200e+02
[ dihedrals ]
; ai aj ak al funct c0 c1 c2
2 1 3 4 1 1.800000e+02 3.347200e+01 2.000000e+00
6 1 3 4 1 1.800000e+02 3.347200e+01 2.000000e+00
2 1 3 5 1 1.800000e+02 3.347200e+01 2.000000e+00
[ dihedrals ]
; ai aj ak al funct c0 c1
3 4 5 1 2 0.000000e+00 1.673600e+02
; Include SPC water topology
#include "spc.itp"
[ system ]
Urea in Water
[ molecules ]
Urea 1
SOL 1000
我的做法是将这些原子添加到力场文件中。如图，根据topol里提供的信息，找到Gromacs提供的所用力场的文件地址：
（PS：这里图截错了，实际应该是#include "amber99sb-ildn.ff/forcefield.itp"）
打开forcefield.itp文件，会发现其还是用include命令调用了两个itp文件。将这三个itp文件拷贝到我的工作目录（存放了protein.gro、ligand.gro、topol什么的文件夹）下的子文件夹ff里：
在ffnonbonded.itp文件里有原子类型信息（用笔记本打开查看），可以将配体itp里的atomtype的数据拷贝过去，最后保存。
之后，修改topol里调用forcefield.itp的命令，将其路径改为调用修改后的forcefield.itp文件存放的位置ff文件夹：
到这里，提示拓扑文件里有[atomtypes]的fatal报错就会消失。
7. 定义盒子大小和添加溶剂分子
定义盒子：
gmx editconf -f complex.gro -o newbox.gro -bt dodecahedron -d 3.0
# 要查看每个参数是什么意思，可以用gmx editconf -h
# 这里修改了d参数为3，原先是1，这个距离在我的体系里原子挨太近了，EM过程报错
添加溶剂水：
gmx solvate -cp newbox.gro -cs spc216.gro -p topol.top -o solv.gro
# 没有tip3p.gro
# 用spc216的gro代替，不过拓扑信息还是用tip3p的水
我觉得不一定要用tip3p.gro，不会报错就不错了。
8. 规避报错，分子数量的报错
这里添加了水到复合物的gro文件里，但是topol文件里的[melocules]还是没有更新的状态，所以有下面的报错：
这里提示了坐标文件中原子数和拓扑文件中的原子数不同。Fatal error: number of coordinates in coordinate file (solv.gro, 29510) does not match topology...-CSDN博客
这篇博客中老师的方法给了我提醒。要在topol文件的molecules里添加水的分子数，这里多出75549个原子，就是25183个水分子（除以3）。
这里为什么是SOL， 其实，这个记录在每个分子的itp文件的[moleculetypes]部分里。如下图，在[dir]\gmx\share\gromacs\top\amber99sb-ildn.ff文件夹下的tip3p.itp文件里的[moleculetype]：
所以，这个水的分子名是SOL。
9. 添加离子
写一个mdp文件，重命名为ions.mdp，真不懂为什么加离子就要用到这种文件。
# 生成一个tpr文件：
gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr
# 添加离子，生成一个新的gro文件：
gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral
# 运行最后一个添加离子的命令后，会出现如下的选项：
Reading file ions.tpr, VERSION 2022.2 (single precision)
......
Group 8 ( SideChain) has 3352 elements
Group 9 ( SideChain-H) has 1174 elements
Group 10 ( Prot-Masses) has 4892 elements
Group 11 ( non-Protein) has 75582 elements
Group 12 ( Other) has 33 elements
Group 13 ( MOL) has 33 elements
Group 14 ( Water) has 75549 elements
Group 15 ( SOL) has 75549 elements
Group 16 ( non-Water) has 4925 elements
Select a group:
# 这里输入15
# 因为15是溶剂，这里会将水分子的位置替换为离子的位置。
输入15，因为我的topol里的溶剂分子类型写的就是SOL，查看topol文件，添加的离子文件已经更新到文件里：
采用VMD或是Pymol检查配体和蛋白之间是否存在挨太近等不合理的状况：
参考教程：
Protein-Ligand Complex