（十二）体系平衡

1.设定对配体的限制：
$gmx/genrestr -f jz4.gro -o posre_jz4.itp -fc 1000 1000 1000

** -f jz4.gro 输入配体的gro文件
** -o posre_jz4.itp 输出加了限制的配体itp文件
** -fc 1000 1000 1000 限制Force constants (kJ/mol nm^2)

2. 将限制的配体itp文件加入到topol.top中并保存：

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

; Include ligand topology
#include "drg.itp"

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

; Include water topology
#include "gromos43a1.ff/spc.itp"

3.
$gmx/make_ndx -f em.gro -o index.ndx

系统显示：
 0 System              : 33037 atoms
 1 Protein             :  1693 atoms
 2 Protein-H           :  1301 atoms
 3 C-alpha             :   163 atoms
 4 Backbone            :   489 atoms
 5 MainChain           :   653 atoms
 6 MainChain+Cb        :   805 atoms
 7 MainChain+H         :   815 atoms
 8 SideChain           :   878 atoms
 9 SideChain-H         :   648 atoms
10 Prot-Masses         :  1693 atoms
11 non-Protein         : 31344 atoms
12 Other               :    15 atoms
13 JZ4                 :    15 atoms
14 CL                  :     6 atoms
15 Water               : 31323 atoms
16 SOL                 : 31323 atoms
17 non-Water           :  1714 atoms
18 Ion                 :     6 atoms
19 JZ4                 :    15 atoms
20 CL                  :     6 atoms
21 Water_and_ions      : 31329 atoms

此时我们将蛋白（1 Protein）与配体（13 JZ4）融合：
> 1 | 13
> q  #退出make_ndx


4. 设定NVT系综下体系平衡的配置文件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
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 JZ4
; Bond parameters
continuation    = no            ; first dynamics run
constraint_algorithm = lincs    ; holonomic constraints
constraints     = all-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_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
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

5.grompp预处理——NVT系综的体系平衡
$gmx/grompp -f nvt.mdp -c em.gro -p topol.top -n index.ndx -o nvt.tpr

6. NVT系综运行体系平衡
$gmx/mdrun -deffnm nvt

7. 设定NPT系综下体系平衡的配置文件npt.mdp

用记事本编辑如下参数，保存为npt.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
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 JZ4
; Bond parameters
continuation    = yes           ; first dynamics run
constraint_algorithm = lincs    ; holonomic constraints
constraints     = all-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_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
refcoord_scaling    = com
; Periodic boundary conditions
pbc         = xyz       ; 3-D PBC
; Dispersion correction
DispCorr    = EnerPres  ; account for cut-off vdW scheme
; Velocity generation
gen_vel     = no        ; velocity generation off after NVT

8.grompp预处理——NPT系综的体系平衡
$gmx/grompp -f npt.mdp -c nvt.gro -t nvt.cpt -p topol.top -n index.ndx -o npt.tpr

9.NPT系综运行体系平衡
$gmx/mdrun -deffnm npt


至此，体系平衡已经大功告成，（温度，压力均已确定）

（十三）动力学模拟

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

$gmx/mdrun -deffnm md_0_1