MD: Plain run
Introduction
How setup a MD simulation in water for a PPI complex.
Instructions based on the Lysozime in Water GROMACS tutorial.
Setup of the environment
Open a terminal in the server (this instruction depend on the
configuration of the server)
term>module load cuda/cuda-9.1
Setup the environment of the terminal
term> export PATH=/opt/gromacs/bin:$PATH
term> export PATH=/opt/biki/BiKiLifeSciences/vmd-1.9.2/vmd_bin/:$PATH
If a setup file is present it is possible to use the instruction source
source $HOME/gromacs_setup.source
In this case the setup file is located in the HOME directory and is named
gromacs_setup.source
If a modulefile is present
module load ge/gromacs
Setup of the system
The starting structure can be downloaded from the PDB database or can be
the result of PPI building sites.
- Use pdb format
- The protein will be completed using GROMACS tools, thus it is better
to use a file without any hydrogen atoms but GROMACS tools can ignore
hydrogen atoms
Prepare gro file from pdb file.
gmx pdb2gmx -ignh -f model_4_MinChimera.pdb -o model_4.gro
-water spce
Select the FF that will be used by GROMACS
- I usually use AMBER99SB that is used by ACTYPE tool
- In lysozime tutorial is suggested OPLS-AA (a
deeper analysis ìs required)
The topology file is created.
gmx editconfig -f model_4.gro -o model_4_box.gro -c -d 1.0
-bt cubic
The box is prepared and the length of the cubic box is fixed.
gmx solvate -cp model_4_box.gro -cs spc216.gro -o
model_4_solv.gro -p topol.top
The solvated box is prepared.
Now, it is possible to check is the dimension of the box and the number
of water molecule is convenient to avoid long calculation. It is possible
to view the system with VMD and check if the box if too large and then
repeat editconf and solvate to obtain a convenient system. For example,
the box dimensions can be fixed using -box flag followed by the dimensions
in nm.
gmx editconf -f model_4.gro -o model_4_box.gro -c -d 1.0
-box 8.0 8.0 8.0
gmx solvate -cp model_4_box.gro -cs spc216.gro -o model_4_solv.gro -p
topol.top
Finally positive and negative ions have to be added to neutralize the
system
gmx grompp -f ions.mdp -c model_4_solv.gro -p topol.top -o
ions.tpr
gmx genion -s ions.tpr -o model_4_solv_ions.gro -p topol.top -pname NA
-nname CL -neutral
Choose the option 13.
The ions.mdp is a script only used for the preparation of the file to
submit to the ion addition
ions.mdp
; ions.mdp - used as input into grompp to generate ions.tpr
; 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 < 1000.0 kJ/mol/nm
emstep =
0.01 ; Minimization
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 ; Buffered
neighbor searching
ns_type =
grid ; Method to determine neighbor list
(simple, grid)
coulombtype = cutoff ; Treatment
of long range electrostatic interactions
rcoulomb =
1.0 ; Short-range electrostatic
cut-off
rvdw =
1.0 ; Short-range Van der Waals
cut-off
pbc
= xyz ; Periodic Boundary Conditions
in all 3 dimensions
Setup of the MD protocol
Minimization step
gmx grompp -f minim.mdp -c model_4_solv_ions.gro -p
topol.top -o em
gmx mdrun -v -deffnm em -cpt 1 -cpo em_restart.cpt
minim.mdp
; minim.mdp - used as input into grompp to generate em.tpr
; 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 < 1000.0 kJ/mol/nm
emstep =
0.01 ; Minimization
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 =
20 ; Frequency to update
the neighbor list and long range forces
cutoff-scheme = Verlet ; Buffered neighbor
searching
ns_type =
grid ; Method to determine neighbor list
(simple, grid)
coulombtype =
PME ; Treatment of long range
electrostatic interactions
rcoulomb =
1.0 ; Short-range electrostatic
cut-off
rvdw =
1.0 ; Short-range Van der Waals
cut-off
pbc
= xyz ; Periodic Boundary Conditions
in all 3 dimensions
Check the energy.
gmx energy -f em.edr -o potential.xvg
At the prompt, select Potential (10); zero (0) terminates input and then
view the results with xmgrace
xmgrace potential.xvg
Run the MD
Equilibration steps
The first equilibration step (nvt) stabilizes the temperature of the
system.
gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o
nvt
gmx mdrun -deffnm nvt -nb gpu
- the flag nb is used to use gpu in the calculaton of non bond energy
nvt.mdp
title
= VSP9 NVT equilibration
define
= -DPOSRES ; position restrain the protein
; 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
; Bond parameters
continuation
= no ; first dynamics run
constraint_algorithm = lincs ;
holonomic constraints
constraints
= h-bonds ; bonds involving H are constrained
lincs_iter
= 1 ; accuracy of LINCS
lincs_order
= 4 ; also related to
accuracy
; Nonbonded settings
cutoff-scheme
= Verlet ; Buffered neighbor searching
ns_type
= grid ; search neighboring grid cells
nstlist
= 10 ; 20 fs, largely irrelevant
with Verlet
rcoulomb
= 1.0 ; short-range electrostatic
cutoff (in nm)
rvdw
= 1.0 ; short-range van der Waals
cutoff (in nm)
DispCorr
= EnerPres ; account for cut-off vdW scheme
; Electrostatics
coulombtype
= PME ; Particle Mesh Ewald for
long-range electrostatics
pme_order
= 4 ; cubic interpolation
fourierspacing =
0.16 ; grid spacing for FFT
; Temperature coupling is on
tcoupl
=
V-rescale
; modified Berendsen thermostat
tc-grps
= Protein Non-Protein ; 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 is off
pcoupl
= no ; no pressure coupling in
NVT
; Periodic boundary conditions
pbc
= xyz ; 3-D PBC
; Velocity generation
gen_vel
= yes ; assign velocities from Maxwell
distribution
gen_temp
= 300 ; temperature for Maxwell
distribution
gen_seed
= -1 ; generate a random seed
Check the temperature with gmx energy tool.
The second equilibration step (npt) stabilizes the pressure of the
sytem.
gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p
topol.top -o npt
gmx mdrun -deffnm npt -nb gpu
npd.mdp
title
= VSP9 NPT equilibration
define
= -DPOSRES ; position restrain the protein
; 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
; Bond parameters
continuation
= yes ; Restarting after NVT
constraint_algorithm = lincs ;
holonomic constraints
constraints
= h-bonds ; bonds involving H are constrained
lincs_iter
= 1 ; accuracy of LINCS
lincs_order
= 4 ; also related to
accuracy
; Nonbonded settings
cutoff-scheme
= Verlet ; Buffered neighbor searching
ns_type
= grid ; search neighboring grid cells
nstlist
= 10 ; 20 fs, largely
irrelevant with Verlet scheme
rcoulomb
= 1.0 ; short-range electrostatic
cutoff (in nm)
rvdw
= 1.0 ; short-range van der Waals
cutoff (in nm)
DispCorr
= EnerPres ; account for cut-off vdW scheme
; Electrostatics
coulombtype
= PME ; Particle Mesh Ewald for
long-range electrostatics
pme_order
= 4 ; cubic
interpolation
fourierspacing =
0.16 ; grid spacing for FFT
; Temperature coupling is on
tcoupl
=
V-rescale
; modified Berendsen thermostat
tc-grps
= Protein Non-Protein ; 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 is on
pcoupl
= Parrinello-Rahman ; Pressure coupling on in
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
; Velocity generation
gen_vel
= no ; Velocity generation is
off
Check pressure and density with gmx energy tool.
Production
Finally, it is possible to run the production step.
gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o
md_0_1
gmx mdrun -deffnm md_0_1 -nb gpu
md.mdp
title
= VSP9 production
; Run parameters
integrator
= md ; leap-frog integrator
nsteps
= 500000 ; 2 * 500000 = 1000 ps (1 ns)
dt
= 0.002 ; 2 fs
; Output control
nstxout
= 0 ; suppress bulky
.trr file by specifying
nstvout
= 0 ; 0 for output
frequency of nstxout,
nstfout
= 0 ; nstvout, and
nstfout
nstenergy
= 5000 ; save energies every 10.0 ps
nstlog
= 5000 ; update log file every 10.0 ps
nstxout-compressed =
5000 ; save compressed coordinates every
10.0 ps
compressed-x-grps =
System ; save the whole system
; Bond parameters
continuation
= yes ; Restarting after NPT
constraint_algorithm = lincs ;
holonomic constraints
constraints
= h-bonds ; bonds involving H are constrained
lincs_iter
= 1 ; accuracy of LINCS
lincs_order
= 4 ; also related to
accuracy
; Neighborsearching
cutoff-scheme
= Verlet ; Buffered neighbor searching
ns_type
= grid ; search neighboring grid cells
nstlist
= 10 ; 20 fs, largely
irrelevant with Verlet scheme
rcoulomb
= 1.0 ; short-range electrostatic
cutoff (in nm)
rvdw
= 1.0 ; 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 is on
tcoupl
=
V-rescale
; modified Berendsen thermostat
tc-grps
= Protein Non-Protein ; 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 is on
pcoupl
= Parrinello-Rahman ; Pressure coupling on in
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
DispCorr
= EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel
= no ; Velocity generation is
off
The length of the simulation depends on the number of the steps.
- dt
= 0.002 ; 2 fs
- nsteps
= 500000 ; 2 fs * 500000 = 1000 ps (1 ns)