最近在做蛋白与小分子复合物的动力学模拟，学习了Gromacs的用法，以此记录一下。这里贴出Gromacs的官方文档http://www.mdtutorials.com/gmx/complex/index.html，有兴趣的同学可以复现一遍。

0. Gromacs在Windows下的安装

打开http://sobereva.com/458网址，找到sobereva老师编译的版本，按需下载。这里我下载的是2020.6 CUDA加速版。下载后将.bin目录添加至环境变量.path中即可。具体步骤如下：

Windows编译版

1. 搜索框搜索编辑系统环境变量；

   
   
   搜索编辑系统变量

2. 打开环境变量；

   
   
   打开环境变量
3. 点击path，再点击编辑
   
   编辑.path
4. 复制下载好的Gromacs文件夹中.bin路径，粘贴至path中，确定关闭对话框。至此，环境变量设置完毕。
   
   添加.bin路径至path中
5. 测试Gromacs能否正常运行。地址栏输入cmd打开终端后，输入gmx，若返回如图所示，则代表安装成果。
   
   gmx测试

1.预处理受体蛋白、配体分子

将蛋白另存为3HTB.pdb，选择分子对接后的排名第一的配体小分子并输出.pdb文件后，用pymol或其他软件另存为jz4.mol2。

2. 准备蛋白拓扑

gmx pdb2gmx -f 3HTB.pdb -o 3HTB_processed.gro -water spc

1. 提示选择立场，这里选择Amber99sb-ildn力场，输入6。

   
   
   选择Amber99sb-ildn力场

2. 输出三个文件，如下所示：

   
   
   输出文件

3.处理小分子生成.gro和.top文件

Gromacs无法对小分子配体进行处理生成拓扑文件，这里使用Sob老师开发的处理小分子软件Sobtop，实用简单。

1. .mol2文件的预处理：使用文本编辑器打开.mol2文件，将图中位置改成配体名称，并保存；
   

   将图中位置的名称改为配体的名称.png

2. 打开Sobtop，将配体分子jz4.mol2拖进去，并回车；

3. 选择1，生成top文件；

   
   

4. 选择3，尽可能使用GAFF力场；

   
   

5. 选择0，进入下一步

   
   

6. 选择4，如果可能，预先构建成键参数，任意猜测缺少的选项

   
   

7. 回车，生成的top文件生成在sobtop软件根目录下

8. 回车，生成的itp位置限制文件在sobtop软件根目录下

9. 选择2，生成gro文件

   
   

10. 回车，生成的gro文件在sobtop软件根目录下

11. 回车，退出sobtop软件

处理过后得到jz4.itp/.gro/.top文件。将这三个文件剪切至工作目录中。

4. 构建复合体及拓扑

1. 构建复合体： 将3HTB_processed.gro复制后粘贴，重命名为complex.gro。将jz4.gro中的内容插入到3HTB_processed.gro中的蛋白原子之后，盒向量之前，并修改第二行的原子总数（修改为蛋白和小分子的原子总和）
   
   绿色部分是加进来的部分.png
   
   
   修改第二行的总和数
2. 构建拓扑：将JZ4的拓扑包含到体系拓扑中去，在topol.top中引入位置限定文件的后面插入，

; Include ligand topology
#include "jz4.itp"

如图所示

插入拓扑

5. 定义盒子、添加溶剂和离子、能量最小化

1. 定义盒子并向盒子充满水。

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，使之成为单独一行

2. 添加离子。由于生命体系中不存在净电荷，必须在体系中添加离子。可随意用一个.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

3. 能量最小化：
   现在系统已组装完毕，创建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. 限制配体
   平衡蛋白配体复合物与平衡水中只包含蛋白的体系非常类似，值得注意的是：
   1、对配体施加限定。
   2、对于温度耦合群的处理。

为了限定配体，需要为他生成一个位置限定拓扑。首先，为JZ4创建一个只包含非氢原子的index group。

gmx make_ndx -f jz4.gro -o index_jz4.ndx
 > 0 & ! a H*
> q

生成index group

选择新创建的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

创建包含蛋白质和配体的组

2. 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文件作少许改动，否则会报错名字不符。

3. 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

结果文件

可以看到结合能非常强！

结合能