分子动力学模拟笔记-GROMACS模拟蛋白质小分子体系（二）

九、限制配体

gmx genrestr -f Ligand.gro -o posre_Ligand.itp  -fc 1000 1000 1000

出现以下信息：

Reading structure file
Select group to position restrain
Group     0 (         System) has    68 elements
Group     1 (          Other) has    68 elements
Group     2 (            MOL) has    68 elements
Select a group:

选哪个都行，其实是一样的。我个人选择了2

得到posre_Ligand.itp

输出配体施加位置限制的itp文件

找到topol.top中的以下部分：

; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif

; Include ligand topology
#include "Ligand.itp"

; Include water topology
#include "amber99sb-ildn.ff/spc.itp"

在配体拓扑的下面、水拓扑的前面插入配体限制：

; Ligand position restraints
#ifdef POSRES
#include "posre_19.itp"
#endif

千万别放错地方了

十、体系整合

gmx make_ndx -f solv_ions.gro -o index.ndx

出现以下信息：

  0 System              : 125394 atoms
  1 Protein             :  3265 atoms
  2 Protein-H           :  1644 atoms
  3 C-alpha             :   202 atoms
  4 Backbone            :   606 atoms
  5 MainChain           :   809 atoms
  6 MainChain+Cb        :   993 atoms
  7 MainChain+H         :  1000 atoms
  8 SideChain           :  2265 atoms
  9 SideChain-H         :   835 atoms
 10 Prot-Masses         :  3265 atoms
 11 non-Protein         : 122129 atoms
 12 Other               :    68 atoms
 13 MOL                 :    68 atoms
 14 NA                  :     3 atoms
 15 Water               : 122058 atoms
 16 SOL                 : 122058 atoms
 17 non-Water           :  3336 atoms
 18 Ion                 :     3 atoms
 19 MOL                 :    68 atoms
 20 NA                  :     3 atoms
 21 Water_and_ions      : 122061 atoms

 nr : group      '!': not  'name' nr name   'splitch' nr    Enter: list groups
 'a': atom       '&': and  'del' nr         'splitres' nr   'l': list residues
 't': atom type  '|': or   'keep' nr        'splitat' nr    'h': help
 'r': residue              'res' nr         'chain' char
 "name": group             'case': case sensitive           'q': save and quit
 'ri': residue index

>

选择你要整合的部分。在本实验中，我将蛋白和小分子绑定在一起。

> 1 | 13
> q

使用的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
nstxout     = 500       ; save coordinates every 1.0 ps
nstvout     = 500       ; save velocities every 1.0 ps
nstenergy   = 500       ; save energies every 1.0 ps
nstlog      = 500       ; update log file every 1.0 ps
energygrps  = Protein MOL
; Bond parameters
continuation    = no            ; first dynamics run
constraint_algorithm = lincs    ; holonomic constraints 
constraints     = h-bonds     ; all bonds (even heavy atom-H bonds) constrained
lincs_iter      = 1             ; accuracy of LINCS
lincs_order     = 4             ; also related to accuracy
; Neighborsearching
cutoff-scheme   = Verlet
ns_type         = grid      ; search neighboring grid cells
nstlist         = 10        ; 20 fs, largely irrelevant with Verlet
rcoulomb        = 1.4       ; short-range electrostatic cutoff (in nm)
rvdw            = 1.4       ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype     = PME       ; Particle Mesh Ewald for long-range electrostatics
pme_order       = 4         ; cubic interpolation
fourierspacing  = 0.16      ; grid spacing for FFT
; Temperature coupling
tcoupl      = V-rescale                     ; modified Berendsen thermostat
tc-grps     = Protein_MOL 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
DispCorr    = EnerPres  ; account for cut-off vdW scheme
; Velocity generation
gen_vel     = yes       ; assign velocities from Maxwell distribution
gen_temp    = 300       ; temperature for Maxwell distribution
gen_seed    = -1        ; generate a random seed

如果有报错，可以根据提示在constraints、energygrps或tc-grps查找问题。一般报错的话就出在这几个地方。

然后我们就即将开始NVT、NPT平衡和正式的分子动力学模拟了。在此之前我想特别提醒一下，强烈建议各位使用北京超算！！！NVT和NPT平衡还好，一两个小时能做出来，但是成品MD模拟总共50万步，用自己的电脑跑真的非常非常非常慢。北京超算注册条件又简单又有200元免费额度，实际上买起来价格也不太贵，RTX3090区是4.5多一点/小时，Tesla性能高的区稍微贵一些，这个看个人需求，总之租3090和Tesla V100都是很好的选择。按照手册登录平台，选Linux系统，把GROMACS安装包和你已经做好的文件上传上去，打开web ssh就可以正常在超算平台上使用GROMACS了，贼给力！

十一、NVT平衡

gmx grompp -f GMXtut-5_nvt.mdp -c em.gro -r em.gro -p topol.top -n index.ndx -o nvt.tpr

如果出现

Fatal error:
Cannot find position restraint file restraint.gro (option -r).
From GROMACS-2018, you need to specify the position restraint coordinate files
explicitly to avoid mistakes, although you can still use the same file as you
specify for the -c option.

则说明你没有用上面那行命令中的-r来指定.gro文件。在早期版本的教程没有关于这个问题的描述，官网最新版教程可以看看。

然后就可以开始运行NVT平衡了！！

gmx mdrun -deffnm nvt -v

敲下回车，开始坐等

完成时出现以下信息：

Writing final coordinates.

               Core t (s)   Wall t (s)        (%)
       Time:    25564.000     6391.000      400.0
                         1h46:31
                 (ns/day)    (hour/ns)
Performance:        1.352       17.752

十二、NPT平衡

gmx grompp -f npt.mdp -c nvt.gro -t nvt.cpt -r nvt.gro -p topol.top -n index.ndx -o npt.tpr

使用的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_MOL 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 mdrun -deffnm npt -v

来进行最后的NPT平衡了！

（坐等again）

Writing final coordinates.
step 50000, remaining wall clock time:     0 s
               Core t (s)   Wall t (s)        (%)
       Time:    10568.000     2642.000      400.0
                         44:02
                 (ns/day)    (hour/ns)
Performance:        3.270        7.339

得到一堆npt开头的文件

至此，两步平衡结束，准备开始最终模拟！

十三、成品MD

gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -n index.ndx -o md_0_1.tpr

得到md_0_1.tpr

gmx mdrun -deffnm md_0_1 -v 

-v可以让你看到现在运行到哪一步

然后就慢慢等呗