Gromacs-蛋白小分子复合物分子动力学模拟实战（包括软件安装
最近在做蛋白与小分子复合物的动力学模拟，学习了Gromacs的用法，以此记录一下。这里贴出Gromacs的官方文档http://www.mdtutorials.com/gmx/complex/index.html，有兴趣的同学可以复现一遍。
0. Gromacs在Windows下的安装
打开http://sobereva.com/458网址，找到sobereva老师编译的版本，按需下载。这里我下载的是2020.6 CUDA加速版。下载后将.bin目录添加至环境变量.path中即可。具体步骤如下：
Windows编译版
-
搜索框搜索编辑系统环境变量；
搜索编辑系统变量
-
打开环境变量；
打开环境变量
- 点击
path，再点击编辑
编辑.path
- 复制下载好的Gromacs文件夹中
.bin路径，粘贴至path中，确定关闭对话框。至此，环境变量设置完毕。
添加.bin路径至path中
- 测试Gromacs能否正常运行。地址栏输入
cmd打开终端后，输入gmx，若返回如图所示，则代表安装成果。
gmx测试
1.预处理受体蛋白、配体分子
将蛋白另存为3HTB.pdb，选择分子对接后的排名第一的配体小分子并输出.pdb文件后，用pymol或其他软件另存为jz4.mol2。
2. 准备蛋白拓扑
gmx pdb2gmx -f 3HTB.pdb -o 3HTB_processed.gro -water spc
-
提示选择立场，这里选择Amber99sb-ildn力场，输入6。
选择Amber99sb-ildn力场
-
输出三个文件，如下所示：
输出文件
3.处理小分子生成.gro和.top文件
Gromacs无法对小分子配体进行处理生成拓扑文件，这里使用Sob老师开发的处理小分子软件Sobtop，实用简单。
-
.mol2文件的预处理：使用文本编辑器打开.mol2文件，将图中位置改成配体名称，并保存；
将图中位置的名称改为配体的名称.png
-
打开Sobtop，将配体分子jz4.mol2拖进去，并回车；
-
选择1，生成top文件；
-
选择3，尽可能使用GAFF力场；
-
选择0，进入下一步
-
选择4，如果可能，预先构建成键参数，任意猜测缺少的选项
-
回车，生成的top文件生成在sobtop软件根目录下
-
回车，生成的itp位置限制文件在sobtop软件根目录下
-
选择2，生成gro文件
-
回车，生成的gro文件在sobtop软件根目录下
-
回车，退出sobtop软件
处理过后得到jz4.itp/.gro/.top文件。将这三个文件剪切至工作目录中。
4. 构建复合体及拓扑
- 构建复合体： 将
3HTB_processed.gro复制后粘贴，重命名为complex.gro。将jz4.gro中的内容插入到3HTB_processed.gro中的蛋白原子之后，盒向量之前，并修改第二行的原子总数（修改为蛋白和小分子的原子总和）
绿色部分是加进来的部分.png
修改第二行的总和数
- 构建拓扑：将JZ4的拓扑包含到体系拓扑中去，在topol.top中引入位置限定文件的后面插入，
; Include ligand topology
#include "jz4.itp"
如图所示
插入拓扑
在文档末尾
Molecules部分加上JZ4配体名称，保存文件。
5. 定义盒子、添加溶剂和离子、能量最小化
- 定义盒子并向盒子充满水。
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
-bt：盒子的形状为菱形十二面体，
-c：复合物距离盒子边缘至少1.0nm。
-cs spc216.gro：表示所有三点水模型都用spc216.gro填充盒子。
运行完命令后，打开.top文件，将最后一行SOL和后面的数字成为单独一行，保存。
回车SOL，使之成为单独一行
- 添加离子。由于生命体系中不存在净电荷，必须在体系中添加离子。可随意用一个
.mdp文件来生成.tpr文件如官网列出的，因为其参数少、易维护。在工作目录新建一个文本文档，将以下内容粘贴进去，修改名为ions.mdp。
; LINES STARTING WITH ';' ARE COMMENTS
title = Minimization ; Title of run
; Parameters describing what to do, when to stop and what to save
integrator = steep ; Algorithm (steep = steepest descent minimization)
emtol = 1000.0 ; Stop minimization when the maximum force < 10.0 kJ/mol
emstep = 0.01 ; Energy step size
nsteps = 50000 ; Maximum number of (minimization) steps to perform
; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist = 1 ; Frequency to update the neighbor list and long range forces
cutoff-scheme = Verlet
ns_type = grid ; Method to determine neighbor list (simple, grid)
rlist = 1.0 ; Cut-off for making neighbor list (short range forces)
coulombtype = cutoff ; Treatment of long range electrostatic interactions
rcoulomb = 1.0 ; long range electrostatic cut-off
rvdw = 1.0 ; long range Van der Waals cut-off
pbc = xyz ; Periodic Boundary Conditions
gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr
将.tpr文件传递给genion程序，提示选择替换溶剂分子的group时选择“SOL”(15)，得到得到solv_ions.gro文件。
gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral
-pname 阳离子的名称
-nname 阴离子的名称
-np 阳离子个数
-nn 阴离子个数
选择SOL
- 能量最小化：
现在系统已组装完毕，创建em.mdp文件，内容如下：
; LINES STARTING WITH ';' ARE COMMENTS
title = Minimization ; Title of run
; Parameters describing what to do, when to stop and what to save
integrator = steep ; Algorithm (steep = steepest descent minimization)
emtol = 1000.0 ; Stop minimization when the maximum force < 10.0 kJ/mol
emstep = 0.01 ; Energy step size
nsteps = 50000 ; Maximum number of (minimization) steps to perform
; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist = 1 ; Frequency to update the neighbor list and long range forces
cutoff-scheme = Verlet
ns_type = grid ; Method to determine neighbor list (simple, grid)
rlist = 1.2 ; Cut-off for making neighbor list (short range forces)
coulombtype = PME ; Treatment of long range electrostatic interactions
rcoulomb = 1.2 ; long range electrostatic cut-off
vdwtype = cutoff
vdw-modifier = force-switch
rvdw-switch = 1.0
rvdw = 1.2 ; long range Van der Waals cut-off
pbc = xyz ; Periodic Boundary Conditions
DispCorr = no
使用grompp创建二进制输入文件：
gmx grompp -f em.mdp -c solv_ions.gro -p topol.top -o em.tpr
能量最小化：
gmx mdrun -v -deffnm em
运行完毕，可以看到体系能很快收敛，能量最小化成功！
体系很快收敛
6. NVT、NPT平衡
- 限制配体
平衡蛋白配体复合物与平衡水中只包含蛋白的体系非常类似，值得注意的是：
1、对配体施加限定。
2、对于温度耦合群的处理。
为了限定配体，需要为他生成一个位置限定拓扑。首先，为JZ4创建一个只包含非氢原子的index group。
gmx make_ndx -f jz4.gro -o index_jz4.ndx
> 0 & ! a H*
> q
生成index group
执行genrestr模块，选择新创建的index group（在index_jz4.ndx文件中是group 3）。
gmx genrestr -f jz4.gro -n index_jz4.ndx -o posre_jz4.itp -fc 1000 1000 1000
选择新创建的group3
下一步需要把这个信息包括在拓扑文件中。有几种方式，取决于我们希望使用的条件。这里介绍第一种：如果我们想在无论蛋白是否被限定的条件下对配体施加限定，在拓扑文件的Include ligand topology之后中加入如下内容：
; Ligand position restraints
#ifdef POSRES
#include "posre_jz4.itp"
#endif
修改后的.top文件如下所示：
修改后的.top文件
必须要把限定这段话加在Include ligand topology之后，Include water topology之前，加在其他任何地方都会报错。
由于对只有几个原子的组在控制其动能的涨落时, 温度耦合算法不够稳定，不要对体系中的每一单个物种使用独立的耦合。创建一个特殊的组, 其中包含蛋白质和配体。
gmx make_ndx -f em.gro -o index.ndx
> 1 | 13
> q
创建包含蛋白质和配体的组
- NVT平衡
创建nvt.mdp文件，将以下部分复制入文件中：
title = Protein-ligand complex NVT equilibration
define = -DPOSRES ; position restrain the protein and ligand
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 50000 ; 2 * 50000 = 100 ps
dt = 0.002 ; 2 fs
; Output control
nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps
nstxout-compressed = 500 ; save coordinates every 1.0 ps
; Bond parameters
continuation = no ; first dynamics run
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds to H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighbor searching and vdW
cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cells
nstlist = 20 ; largely irrelevant with Verlet
rlist = 1.2
vdwtype = cutoff
vdw-modifier = force-switch
rvdw-switch = 1.0
rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Protein_JZ4 Water_and_ions ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 300 300 ; reference temperature, one for each group, in K
; Pressure coupling
pcoupl = no ; no pressure coupling in NVT
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction is not used for proteins with the C36 additive FF
DispCorr = no
; Velocity generation
gen_vel = yes ; assign velocities from Maxwell distribution
gen_temp = 300 ; temperature for Maxwell distribution
gen_seed = -1 ; generate a random seed
这里需要对nvt.mdp文件作少许改动，否则会报错名字不符。
运行：
gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -n index.ndx -o nvt.tprgmx mdrun -deffnm nvt
- NPT平衡
新建文件nvt.mdp，将以下内容粘贴入文件中。
title = Protein-ligand complex NPT equilibration
define = -DPOSRES ; position restrain the protein and ligand
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 50000 ; 2 * 50000 = 100 ps
dt = 0.002 ; 2 fs
; Output control
nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps
nstxout-compressed = 500 ; save coordinates every 1.0 ps
; Bond parameters
continuation = yes ; continuing from NVT
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds to H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighbor searching and vdW
cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cells
nstlist = 20 ; largely irrelevant with Verlet
rlist = 1.2
vdwtype = cutoff
vdw-modifier = force-switch
rvdw-switch = 1.0
rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
rcoulomb = 1.2
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Protein_JZ4 Water_and_ions ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 300 300 ; reference temperature, one for each group, in K
; Pressure coupling
pcoupl = Berendsen ; pressure coupling is on for NPT
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
refcoord_scaling = com
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction is not used for proteins with the C36 additive FF
DispCorr = no
; Velocity generation
gen_vel = no ; velocity generation off after NVT
这里需要修改少许部分，同上。
运行：
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
7. 运行MD
在温度和压力平衡结束后，体系已经趋于稳定。现在可以释放位置限定，执行MD来收集数据。
创建md.mdp文件，将以下内容粘贴至内容中：
title = Protein-ligand complex MD simulation
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 5000000 ; 2 * 5000000 = 10000 ps (10 ns)
dt = 0.002 ; 2 fs
; Output control
nstenergy = 5000 ; save energies every 10.0 ps
nstlog = 5000 ; update log file every 10.0 ps
nstxout-compressed = 5000 ; save coordinates every 10.0 ps
; Bond parameters
continuation = yes ; continuing from NPT
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds to H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighbor searching and vdW
cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cells
nstlist = 20 ; largely irrelevant with Verlet
rlist = 1.2
vdwtype = cutoff
vdw-modifier = force-switch
rvdw-switch = 1.0
rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
rcoulomb = 1.2
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Protein_JZ4 Water_and_ions ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 300 300 ; reference temperature, one for each group, in K
; Pressure coupling
pcoupl = Parrinello-Rahman ; pressure coupling is on for NPT
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction is not used for proteins with the C36 additive FF
DispCorr = no
; Velocity generation
gen_vel = no ; continuing from NPT equilibration
这里模拟的是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
①如果要使用gpu加速，可在末尾加上 -nb gpu 来使用gpu加速计算，-v 可显示计算结束时间。
gmx mdrun -deffnm md_0_10 -nb gpu -v
②如果模拟过程中因突发状况导致模拟终止，可使用-cpi命令续跑。
在同一个文件夹下运行。
gmx mdrun -s md_0_10.tpr -cpi md_0_10.cpt -deffnm md_0_10
运行后基本分析
提取能量数据观察动力学能量是否达到平衡，如下：
温度：Temperature
gmx energy -f md_0_10.edr -o md_0_10_Temperature.xvg
选择14 0
压力：Pressure
gmx energy -f md_0_10.edr -o md_0_10_Pressure.xvg
选择16 0
密度：Density
gmx energy -f md_0_10.edr -o md_0_10_Density.xvg
选择22 0
势能：Potential
gmx energy -f md_0_10.edr -o md_0_10_Potential.xvg
选择10 0
总能量：Total-Energy
gmx energy -f md_0_10.edr -o md_0_10_TotalEnergy.xvg
选择12 0
结果采用qtgrace可视化，观察中线是否有波动，是否达到平衡。
8. RMSD和RMSF分析
正如任何在周期性边界条件下进行的模拟一样，分子可能会被盒子截断或者反复穿过盒子。为了将蛋白重新置于盒子中央并且将晶胞恢复成菱形十二面体，使用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选择“Protein”，输出选择“System”（0）。
为了更平滑的可视化，进行适当的旋转和平移拟合可能是有益的。按如下方式执行trjconv：
gmx trjconv -s md_0_10.tpr -f md_0_10_center.xtc -o md_0_10_fit.xtc -fit rot+trans
选择“Backbone”进行蛋白质主链的最小二乘拟合，选择“System”进行输出。
由于原始文件太大，我们间隔抽取运动轨迹。
gmx trjconv -s md_0_10.tpr -f md_0_10_fit.xtc -o whole.pdb -dt 500 -b 1000 -e 20000 -n index.ndx
选择“蛋白-配体”进行输出。
RMSD分析
gmx rms -s md_0_10.tpr -f md_0_10_center.xtc -tu ns -o rmsd.xvg
执行rms模块，选择“Backbone”（4）进行最小二乘拟合和RMSD计算，得到蛋白质骨架的前后偏移程度。
也可以选择配体（13），查看配体的前后偏移程度。
RMSF分析
gmx rmsf -s md_0_10.tpr -f md_0_10_center.xtc -o Tanpp_rmsf.xvg -ox Tanpp_avg.pdb -res -oq Tabpp_bfac.pdb
选择1（蛋白质）。
将结果.xvg拖入qtgrace软件可视化。
蛋白与配体氢键数量分析
gmx hbond -s md_0_10.tpr -n index.ndx -f md_0_10.xtc -num md_0_10_HB-Num_Protein_mol.xvg
先选择protein（1），回车后再选择配体（13）
回旋半径分析蛋白的密实度（Radius of gyration）
gmx gyrate -s md_0_10.tpr -n index.ndx -f md_0_10.xtc -o md_0_10_Gyrate-Protein.xvg
选择protein（1）
MMPBSA计算（结合能计算）
采集40-50ns的模拟数据（这里通常采集体系平衡的数据）。
gmx trjconv -f md_0_10_center.xtc -o trj_40-50ns.xtc -b 40000 -e 50000
检查轨迹路径
gmx check -f trj_40-50ns.xtc
检查index.ndx， 配体的索引只能是一个，多的必须删除（ctrl+F搜索配体名称，剩一个配体即可）。删除后保存。
下载 gmx_mmpbsa.bash，网址：https://github.com/Jerkwin/gmxtools/blob/master/gmx_mmpbsa/gmx_mmpbsa.bsh
修改文件：
源文件
修改后
在文件夹下右键-Open Git bash Here，运行bash gmx_mmpbsa.bsh
image.png
打开结果
结果文件
可以看到结合能非常强！
结合能
