Gromacs软件进行蛋白-配体体系MD模拟

1、进行能量最小化

下载官方提供的参数文件：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

 mdp文件里分3个区块，标题、模拟的步数怎么跑和原子之间相互作用力设置的有关参数。

例如模拟的步数怎么跑里的设置：能量下降的积分算法用最陡下降最小化。体系中最大能量达到1000停止最小化。每步下降能量大小是0.01以及最大的步数设置为50000。

gmx grompp -f em.mdp -c solv_ions.gro -p topol.top -o em.tpr

先编译执行文件tpr，然后调用mdrun执行能量最小化：

gmx mdrun -v -deffnm em
#-v 将详细内容输出到终端
#-deffnm em 指定输出和输入文件名都叫em

 查看输出的结果：

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

输出的结果上看，算法在第143步就实现体系中最大能量小于1000，然后结束最小化。这里的力是指原子在原位上受到的键作用力和非键作用力的总和。

这里的Maximum force=800 on atom 27说明体系中有最大能量是27号原子。

2、完成能量最小化后，建立一个需要限制的配体原子的信息的文件。

我的配体的gro文件是ligand.gro，在命令行输入：

gmx make_ndx -f ligand.gro -o index_ligand.ndx #进入交互界面

 > 0 & ! a H* #选择0的system中全部非H原子
 > q #退出

index_ligand.ndx文件内容入下,就是配体原子的原子数据索引：

[ System ] #整个系统的配体原子序号（索引）
   1    2    3    4    5    6    7    8    9   10   11   12   13   14   15
  16   17   18   19   20   21   22   23   24   25   26   27   28   29   30
  31   32   33
......
[ System_&_!H* ] #选出的要限制的非H配体原子序号（索引）
   1    2    3    4    5    6    7    8    9   10   11   12   14   15   16
  17   18   19   20   21   22   23   24

建立对选出的非H原子带有限制力信息的itp文件（拓扑文件）

gmx genrestr -f ligand.gro -n index_ligand.ndx -o posre_ligand.itp -fc 1000 1000 1000
#生产出itp拓扑文件，-fc指出限制力，对ndx里的选出的非H原子用在笛卡尔坐标中各个轴用1000的力栓住
#该命令会进入交互界面
#选出要限制的原子就是了

 在topol文件里加入导入限制力参数posre_l.itp文件的命令：

3、进行温度耦合

温度耦合就是对原先没有温度概念的系统附加上温度。对每种分子类型和它自己的恒温基团耦合是一个坏主意。Gromacs的恒温器不适合对小分子单独进行耦合

 因为gmx的温度耦合算法不稳定，少量原子组成的组别单独进行单独进行耦合，动能波动会很大，会导致体系的崩溃。例如教程里提到的，如果在nvt.mdp里设置温度耦合的组别是这样：

tc-grps = Protein JZ4 SOL CL

SOL是水，Protein是蛋白，这两的原子数目量多（超过1000），单独进行温度耦合是不会有动能波动大的问题的。但是离子Cl-和配体JZ4，体系里CL只有6个，JZ4是22个（教程中的数据），原子数目完全不是体系中蛋白和水分子在一个量级。不能单独进行温度耦合。

JZ4是配体，和蛋白挨得很近，可以将其视为一个整体进行温度耦合，离子则是苟在水分子群里，可以和所有水分子组成一个整体进行耦合。典型的tc-grps设置是：

tc-grps=Protein Non-Protein  #典型的设置

现在就是要把配体和蛋白整合为Protein，水和离子整合为Non-Protein。

首先，采用make_ndx模块将蛋白质和配体的index抽出，做为一个区块[Protein_JZ4]，这个名字是默认的，输出到文件index.ndx里。这个区块的名字要记住，可以在输出的index.ndx里看到。

gmx make_ndx -f em.gro -o index.ndx
.......  #下面是选项
 0 System              : 33506 atoms
  1 Protein             :  2614 atoms
......
 13 JZ4                 :    22 atoms
......
 21 Water_and_ions      : 30870 atoms

> 1 | 13 #这里选中蛋白和配体， 这里已经存在water_and_ions，就不需要再抽出index。
> q

在我的index.ndx文件里存在water_and_ions:

同时也存在蛋白和配体的组和（我的配体是分子名叫MOL，不是JZ4）：

下载官方的nvt.mdp文件，或者使用下面的：

注意mdp里的这个设置，要根据自己体系的分子名进行调整：

tc-grps=Protein_JZ4 Water_and_ions  #这是这个体系的设置，修改在nvt.mdp里

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模拟：

#生成执行文件tpr：
gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -n index.ndx -o nvt.tpr

#-n参数指定index文件，里面是
gmx mdrun -deffnm nvt

4、进行压力耦合

下载官方的参数文件npt.mdp，注意修改里面的tc-grps参数，然后运行模拟：

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

5、成品模拟：

经过温度耦合和压力耦合后，体系达到预期的室温和常压条件（旧制条件），即是温度在25°C，298.15K，压强在1.01bar，101.325kPa条件。

下载官方提供的md.mpt文件，可以进行成品模拟：

注意在这个md.mpt文件里是没有限制力参数的，没有下面这段参数：

define = -DPOSRES  ; position restrain the protein and ligand

 这个参数是出现在nvt和npt模拟中，表示使用拓扑文件中的这个函数，对配体进行位置约束：

; Ligand position restraints
#ifdef POSRES   #函数名是POSRES，mdp里就是DPOSRES,前面多个D
#include "posre_jz4.itp"
#endif

另外，教程里的拓扑topol里还有水分子的位置限制参数，如下：

#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
;  i funct       fcx        fcy        fcz
   1    1       1000       1000       1000
#endif

但是mdp文件并没有define=-DPOSRES_WATER的参数。 

运行没有约束的成品10ns模拟：

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

 参考：

Protein-Ligand Complex (mdtutorials.com)

关于.mdp文件中的DPOSRES的问题 - 分子模拟 (Molecular Modeling) - 计算化学公社 (keinsci.com)

GROMACS中文手册：第三章　算法|Jerkwin