GROMACS-A Molecular Dynamic Approach

The level of detail in this tutorial will be focused on membrane protein-specific considerations, and will not provide an exhaustive explanation of every step, as with the Lysozyme tutorial.

gmx pdb2gmx -f KALP-15_princ.pdb -o KALP-15_processed.gro -ignh -ter -water spc

To use the parameters in lipid.itp, we will have to make some changes to our pre-packaged GROMOS96 53A6 force field files (in $GMXLIB/gromos53a6.ff). Make a copy of this directory in your working directory called “gromos53a6_lipid.ff” (assuming you have GROMACS installed in /usr/local/gromacs):
cp -r /usr/local/gromacs/share/gromacs/top/gromos53a6.ff/ ./gromos53a6_lipid.ff

Inside this directory, you will find the following files:
aminoacids.c.tdb
aminoacids.hdb
aminoacids.n.tdb
aminoacids.r2b
aminoacids.rtp
aminoacids.vsd
atomtypes.atp
ff_dum.itp
ffbonded.itp
ffnonbonded.itp
forcefield.doc
forcefield.itp
ions.itp
spc.itp
spce.itp
tip3p.itp
tip4p.itp
watermodels.dat

Next, modify forcefield.doc to update the description of the force field parameters in it based on the modifications we will make in this tutorial. Mine contains something like:
GROMOS96 53A6 force field, extended to include Berger lipid parameter.s

The files contained in the gromos53a6_lipid.ff directory constitute the complete description of the force field. The files serve the following purposes:
aminoacids.c.tdb and aminoacids.n.tdb – These files are read by pdb2gmx and list the available patching actions that can be applied to protein chains
aminoacids.hdb – This file is read by pdb2gmx and contains a database for constructing hydrogen atoms
aminoacids.r2b – This file is read by pdb2gmx and contains translations between conventional residue names and any force field-specific building block names
aminoacids.rtp – This file is read by pdb2gmx and contains all the available protein residues (and some cofactors, water, and ions – typical contents of a coordinate file passed to pdb2gmx)
atomtypes.atp – This file is read by pdb2gmx and contains all the available atom types in the force field.

ff_dum.itp – This file contains constants that are used in constructing virtual sites (more on this in a different tutorial, and not relevant for our purposes here
ffbonded.itp – All the parameters for bonded interactions in the force field
ffnonbonded.itp – All the nonbonded (Lennard-Jones) parameters in the force field
ions.itp – Contains [moleculetype] definitions for all monoatomic ions in the force field
spc.itp, spce.itp, tip3p.itp, and tip4p.itp – Water model topologies (the standard water model for GROMOS96 simulations is SPC)

watermodels.dat – A text file read by pdb2gmx for interactively selecting the water model to be used in the system topology

Now, to add the lipid parameters into the parent force field, we will need to copy and paste the entries in the [atomtypes], [nonbond_params], and [pairtypes] sections from lipid.itp into the corresponding headings within ffnonbonded.itp. You will find that the lipid.itp [atomtypes] section lacks atomic numbers (the at.num column), so add these in. The newly-modified lines should be:
LO 8 15.9994 0.000 A 2.36400e-03 1.59000e-06 ;carbonyl O, OPLS
LOM 8 15.9994 0.000 A 2.36400e-03 1.59000e-06 ;carboxyl O, OPLS
LNL 7 14.0067 0.000 A 3.35300e-03 3.95100e-06 ;Nitrogen, OPLS
LC 6 12.0110 0.000 A 4.88800e-03 1.35900e-05 ;Carbonyl C, OPLS
LH1 6 13.0190 0.000 A 4.03100e-03 1.21400e-05 ;CH1, OPLS
LH2 6 14.0270 0.000 A 7.00200e-03 2.48300e-05 ;CH2, OPLS
LP 15 30.9738 0.000 A 9.16000e-03 2.50700e-05 ;phosphor, OPLS
LOS 8 15.9994 0.000 A 2.56300e-03 1.86800e-06 ;ester oxygen, OPLS
LP2 6 14.0270 0.000 A 5.87400e-03 2.26500e-05 ;RB CH2, Bergers LJ
LP3 6 15.0350 0.000 A 8.77700e-03 3.38500e-05 ;RB CH3, Bergers LJ
LC3 6 15.0350 0.000 A 9.35700e-03 3.60900e-05 ;CH3, OPLS
LC2 6 14.0270 0.000 A 5.94700e-03 1.79000e-05 ;CH2, OPLS

In the [nonbond_params] section, you will find the line “;; parameters for lipid-GROMOS interactions.” Delete this line and all of the subsequent lines in this section (preserving the section that starts with “;; lipid-SPC/SPCE interactions”). The protein-specific nonbonded combinations are specific to the deprecated “ffgmx” (a modified GROMOS87) force field. Removing these lines allows the interactions between the protein and the lipids to be generated by the standard combination rules of GROMOS96 53A6. Non-bonded interactions involving atom type HW are also present; since these are all zero you can delete these lines as well, or otherwise rename HW as H to be consistent with the GROMOS96 53A6 naming convention. If you do not rename these lines or remove them, grompp will later fail with a fatal error.
Append the contents of the [ dihedraltypes ] to the corresponding section of ffbonded.itp. Do not be concerned that these lines look a bit different. They are Ryckaert-Bellemans dihedrals, which differ in form from the standard periodic dihedrals used in the default GROMOS96 53A6 force field.
To recap, you must make all of the following changes for the force field to be functional:
Copy [atomtypes] from lipid.itp to ffnonbonded.itp and add a column for atomic number
Copy [nonbond_params] from lipid.itp to ffnonbonded.itp
Remove the ;; parameters for lipid-GROMOS interactions and all subsequent lines in the [nonbond_params] section in ffnonbonded.itp
Remove all lines containing “HW” in [nonbond_params] or otherwise rename them to “H”
Copy [pairtypes] from lipid.itp to ffnonbonded.itp
Copy [dihedraltypes] from lipid.itp to ffbonded.itp
Finally, change the #include statement in your topol.top from:
#include “gromos53a6.ff/forcefield.itp”

to:
#include “gromos53a6_lipid.ff/forcefield.itp”

2.gmx grompp -f minim.mdp -c dppc128.pdb -p topol_dppc.top -o dppc.tpr
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 ; 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 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 ; Short-range electrostatic cut-off rvdw = 1.2 ; Short-range Van der Waals cut-off pbc = xyz ; Periodic Boundary Conditions3. gmx trjconv -s dppc.tpr -f dppc128.pdb -o dppc128_whole.gro -pbc mol -ur compact 4. gmx editconf -f KALP-15_processed.gro -o KALP_newbox.gro -c -box 6.41840 6.44350 6.59650 5. gmx editconf -f protein.gro -o protein_newbox.gro -box (membrane box vectors) -center x y z 6.cat KALP_newbox.gro dppc128_whole.gro > system.gro
7.Now, generate this new position restraint file using genrestr:
gmx genrestr -f KALP_newbox.gro -o strong_posre.itp -fc 100000 100000 100000
8.perl inflategro.pl system.gro 4 DPPC 14 system_inflated.gro 5 area.dat
9.gmx grompp -f minim_inflategro.mdp -c system_inflated.gro -p topol.top -r system_inflated.gro -o system_inflated_em.tpr

10.gmx mdrun -deffnm system_inflated_em
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
define = -DSTRONG_POSRES ; Prevent protein from moving
integrator = steep ; Algorithm (steep = steepest descent minimization)
emtol = 1000.0 ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm 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 ; Buffered neighbor searching ns_type = grid ; Method to determine neighbor list (simple, grid) rlist = 1.2 ; Cut-off for making neighbor list (short range forces) coulombtype = cutoff ; Treatment of long range electrostatic interactions rcoulomb = 1.2 ; Short-range electrostatic cut-off rvdw = 1.2 ; Short-range Van der Waals cut-off pbc = xyz ; Periodic Boundary Conditions11.gmx solvate -cp system_shrink26_em.gro -cs spc216.gro -o system_solv.gro -p topol.top 12.perl water_deletor.pl -in system_solv.gro -out system_solv_fix.gro -ref O33 -middle C50 -nwater 3 13.gmx grompp -f ions.mdp -c system_solv_fix.gro -p topol.top -o ions.tpr 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) rlist = 1.2 ; Cut-off for making neighbor list (short range forces) coulombtype = cutoff ; Treatment of long range electrostatic interactions rcoulomb = 1.2 ; Short-range electrostatic cut-off rvdw = 1.2 ; Short-range Van der Waals cut-off pbc = xyz ; Periodic Boundary Conditions in all 3 dimensions14.gmx genion -s ions.tpr -o system_solv_ions.gro -p topol.top -pname NA -nname CL -neutral15.gmx grompp -f minim.mdp -c system_solv_ions.gro -p topol.top -o em.tpr 16.gmx mdrun -v -deffnm em 17.gmx make_ndx -f em.gro -o index.ndx ... > 1 | 13
> q
18.gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -n index.ndx -o nvt.tpr
Nvt.mdp
title = NVT equilibration for KALP15-DPPC
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 = 100 ; save coordinates every 0.2 ps
nstvout = 100 ; save velocities every 0.2 ps
nstenergy = 100 ; save energies every 0.2 ps
nstlog = 100 ; update log file every 0.2 ps
; 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
ns_type = grid ; search neighboring grid cels
nstlist = 5 ; 10 fs
rlist = 1.2 ; short-range neighborlist cutoff (in nm)
rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)
rvdw = 1.2 ; 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_DPPC Water_and_ions ; two coupling groups – more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 323 323 ; 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
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = yes ; assign velocities from Maxwell distribution
gen_temp = 323 ; temperature for Maxwell distribution
gen_seed = -1 ; generate a random seed
; COM motion removal
; These options remove motion of the protein/bilayer relative to the solvent/ions
nstcomm = 1
comm-mode = Linear
comm-grps = Protein_DPPC Water_and_ions

19.gmx mdrun -deffnm nvt
20. gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -n index.ndx -o npt.tpr (copy script from online for npt.mdp)

21.gmx mdrun -deffnm npt
22.gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -n index.ndx -o md_0_1.tpr (copy script from online for md.mdp)
23.gmx mdrun -deffnm md_0_1