软件与材料
本文内容使用到如下软件：
- Ambertools 20
- GROMACS 2019.5
- Gaussian 16
如果使用win10操作系统，可以自己编译前两个软件，也可以找别人编译好的版本。比如Ambertools、GROMACS2019.5；在此表示感谢。Gaussian 16 是要付费购买的，在本文中是用于计算配体的RESP电荷，RESP电荷还可以通过Multiwfn得到；配体还可以采用bcc电荷(准确性不如RESP)，可以直接用Ambertools计算得到。
本文还用到一个蛋白protein.pdb和一个从drugbank下载的配体小分子ligand.pdb。蛋白将利用GROMACS的力场处理，而配体将用amber的GAFF力场处理。
配体处理
配体加氢
本次下载得到的ligand.pdb文件里是没有包含非极性氢的，所以需要给它加氢。加氢的方式有很多，比如Ambertools里的reduce、或者pymol等。
reduce ligand.pdb > ligand_h.pdb
reduce 可能更适用于氨基酸组成的分子；不论使用什么软件进行加氢操作，最后一定要检查结果文件中是否包含了加上的氢。还需要检查加氢之后的结构是否合理，不合理的结构无法通过高斯进行计算。
之后通过Ambertools的antechamber处理成Gaussian的输入文件。
antechamber -i ligand_h.pdb -fi pdb -o ligand.gjf -fo gcrt -pf y
# 当然输入的配体文件也可以是mol2等
# antechamber -i liand_h.mol2 -fi mol2 -o ligand.gjf -fo gcrt -pf y
然后利用Gaussian 16 计算：
g16 < ligand.gjf > ligand.out
在Gaussian的输出文件中就有ESP电荷数据；利用antechamber得到包含RESP电荷的mol2文件：
antechamber -i ligand.out -fi gout -c resp -o ligand_resp.mol2 -fo mol2 -pf y
要在Ambertools中对配体应用GAFF力场，我们还需要生成一个输入文件：
parmchk -i ligand_resp.mol2 -f mol2 -o ligand.frcmod
ligand.frcmod中包含小分子的一些必要参数。
如此，我们所需要的的Ambertools输入文件就准备好了。
对配体应用GAFF力场
打开Ambertools中的tleap，输入如下命令：
# 载入力场
source leaprc.gaff
# 载入配体参数化文件
loadamberparams ligand.frcmod
# 载入配体
lig = loadmol2 ligand_resp.mol2
# 检查配体
check lig
# 保存amber输入文件
saveamberparm lig lig.prmtop lig.inpcrd
转换配体结构和拓扑文件到GROMACS格式
接下来我们需要将amber的输入文件转换为GROMACS可用的输入文件，也即top文件和gro文件。
通过Ambertools中的acpype程序进行转换：
# 输出帮助信息
# acpype -h
# 做文件转换，输出文件名默认为Ambertools中unit的名字
acpype -p lig.prmtop -x lig.inpcrd
如此我们就得到了top文件和gro文件；配体的top文件还需要一点修改才可以被GROMACS使用。
构建模拟体系
生成蛋白质拓扑文件
我们已经有了蛋白文件protein.pdb，配体的top和gro文件；接下来首先要做的就是创建蛋白质的拓扑文件：
gmx pdb2gmx -f protein.pdb -o protein.gro -water spc -ignh
# 选择力场
# 这里选择amber99sb力场
然后我们在相应目录下就可以看到生成了protein.gro、topol.top、posre.itp。topol.top里面包含了选择的力场的一些参数以及蛋白质的拓扑结构参数等等。
接下来我们需要做一系列文件更改，来把配体加到我们的体系中。
拓扑文件修改
先说明本例子中各文件的修改步骤：
- 复制protein.gro并重命名为complex.gro
- 将配体的gro文件的坐标部分，复制到complex.gro的坐标部分后面，并相应修改第二行的原子数目
- 将配体的top文件重命名为ligand.itp
- 将pdb2gmx命令生成的topol.top文件重命名为protein.itp
- 新建文件topol.top
- 将protein.itp中
[ moleculetype ]部分之前的内容完全剪切到topol.top - 将protein.itp中
[ dihedrals ]部分之后的内容完全剪切到topol.top，包括各种文件的导入部分、[ system ]以及[ molecules ]。 - 在topol.top中导入蛋白质位置限制文件的前面，导入蛋白质拓扑文件protein.itp
- 对照ligand.itp中最后的
[ molecules ]，在topol.top中[ molecules ]部分添加配体的名字和数目；在topol.top导入水拓扑文件的前面导入配体拓扑文件ligand.itp。 - 删去ligand.itp中最后的
[ system ]和[ molecules ]部分 - 适当修改topol.top中[ system ]的名字（这个可能无所谓）
- 对于ligand.itp中
[ moleculetype ]前面的内容，我们也需要整合到topol.top中，一般包含两部分，[ defaults ]和[ atomtypes ]。对于这两部分内容，直接复制到topol.top文件中导入蛋白质拓扑文件部分的前面即可。 - 现在我们topol.top文件的最开始，是导入力场参数文件，之后就是配体的
[ defaults ]部分，由于力场文件中通常也包含这一部分，所以会出错，将[ defaults ]部分注释掉；然后尝试对整个体系进行溶剂化、电荷平衡操作，如果不报错，可能就可以了，如果报错，那可能还需要进一步修改topol.top的力场部分。比如用配体的[ defaults ]部分置换掉力场的[ defaults ]部分等。 - 有关力场文件的详细信息，可参考jerkwin的博文，主要是5.7.1部分；在此表示感谢。
关于本体系中拓扑文件修改的详细情况，说明如下。
首先，在配体和蛋白质的拓扑文件（.top）中，主要包含如下几个部分的内容：
- [ defaults ]
- [ atomtypes ]
- [ moleculetype ]
- [ atoms ]
- [ bonds ]
- [ pairs ]
- [ angles ]
- [ dihedrals ]
- [ system ]
- [ molecules ]
前两个部分[ defaults ] 和[ atomtypes ] 主要是相应的力场参数的部分，后面的部分则定义了一些分子拓扑结构相关的参数以及系统参数。一般的itp文件的内容则主要包括3-8的部分，所以我们把ligand.top重命名为ligand.itp、pdb2gmx生成的topol.top重命名为protein.itp之后需要去掉多余的部分；但是去掉的部分并不是不需要的，这部分内容需要与体系的topol.top文件进行整合。
首先将protein.itp的1-2、9-10部分剪切到空白的topol.top中，并include蛋白质的itp文件，得到topol.top如下：
; Include forcefield parameters
#include "amber99sb.ff/forcefield.itp"
; Include protein topology
#include "protein.itp"
; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif
; Include water topology
#include "amber99sb.ff/spc.itp"
#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
; i funct fcx fcy fcz
1 1 1000 1000 1000
#endif
; Include topology for ions
#include "amber99sb.ff/ions.itp"
[ system ]
; Name
Protein
[ molecules ]
; Compound #mols
Protein 1然后我们需要将ligand.itp中的内容整合到topol.top中来。
首先是整合[ system ]和[ molecules ]部分，整合后的topol.top的这两部分如下：
[ system ]
; Name
Protein_and_ligand
[ molecules ]
; Compound #mols
Protein 1
MOL 1
;主要就是添加了配体这一行，这里配体名为MOL整合之后记得删除ligand.itp中的这两部分，并在topol.top导入水拓扑文件的前面导入ligand.itp。
然后整合ligand.itp的前两个部分到topol.top中，直接剪切并粘贴到topol.top导入蛋白质拓扑文件前面、导入力场的后面，然后注释掉[ default ]部分。
整合之后的topol.top如下：
; Include forcefield parameters
#include "amber99sb.ff/forcefield.itp"
; [ defaults ]
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
; 1 2 yes 0.5 0.8333
[ atomtypes ]
;name bond_type mass charge ptype sigma epsilon Amb
c2 c2 0.00000 0.00000 A 3.39967e-01 3.59824e-01 ; 1.91 0.0860
c1 c1 0.00000 0.00000 A 3.39967e-01 8.78640e-01 ; 1.91 0.2100
os os 0.00000 0.00000 A 3.00001e-01 7.11280e-01 ; 1.68 0.1700
o o 0.00000 0.00000 A 2.95992e-01 8.78640e-01 ; 1.66 0.2100
; Include protein topology
#include "protein.itp"
; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif
; Include ligand topology
#include "ligand.itp"
; Include water topology
#include "amber99sb.ff/spc.itp"
#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
; i funct fcx fcy fcz
1 1 1000 1000 1000
#endif
; Include topology for ions
#include "amber99sb.ff/ions.itp"
[ system ]
; Name
Protein_and_ligand
[ molecules ]
; Compound #mols
Protein 1
MOL 1后续建模
通过GROMACS进行后续的建模，如果在添加离子等后续步骤中不报错的话，体系一般就是ok的。
gmx editconf -f complex.gro -o newbox.gro -bt cubic -d 1.0
gmx solvate -cp newbox.gro -cs spc216.gro -p topol.top -o solv.gro
gmx grompp -f em.mdp -c solv.gro -r solv.gro -p topol.top -o ions.tpr -maxwarn 2
gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -np 5关于力场
这种体系中，小分子的力场最好选择和蛋白质的力场比较相近的力场；也即同一个体系中，所有的物质，最好都采用同一个力场或者同一类力场。本文中，GAFF力场和amber99sb算是同一类力场吧。
如果采用不同类型的力场（如假设本文中蛋白质采用GROMOS系列的力场），可能会有一些理论上的问题，请咨询相关的资深人员。
关于多配体体系
如果蛋白配体体系中需要有多个位置不同的配体，那可以通过pymol等工具，载入ligand_resp.mol2文件，通过复制移动等功能实现蛋白周围的不同位置配体的生成，然后分别保存每个配体为mol2文件。这样就可以得到多个有resp电荷的配体文件，但要注意：pymol等工具导出的文件中原子类型可能和ligand_resp.mol2中定义的「原子类型」不同了，需要对照ligand_resp.mol2文件对原子类型一列进行修改，才能与ligand.frcmod配合进行后续建模。之后再通过Ambertools处理，acpype转换。最后整合和include到体系中应该就可以了。
关注公众号蓝色火烈鸟，了解更多技术笔记和工具推荐：
