1E0Q with gromacs

The following is a procedure to run dynamic simulation on a protein with Gromacs 5.1.2, all module files in older version are combine with only one binary file gmx(or gmx_mpi) as a flag, some modules might deceased in new version.

Here is tutorial of gromacs, you can find many help there. This is another tutorial about how to run gromacs in older version, they gave more explaination on the parameter files and results analysis.

Download 1E0Q.pdb from PDB.

minCHARMM.pl 1E0Q.pdb >1e0q.pdb
gmx_mpi pdb2gmx -f 1e0q.pdb -o 1e0q.gro -p 1e0q.top

(output 1e0q.gro/top, posre.itp)

The pdb2gmx command creates input coordinates and topology for Gromacs to use. We choose CHARMM27 force field and the TIP3P water model.

Open 1e0q.top, we can find qtot=0,  system only show protein.

 

Now we set up a box and fill it with solvent.

gmx_mpi editconf -f 1e0q.gro -o 1e0q_pbc.gro -bt cubic -d 1.5

gmx_mpi solvate -cp 1e0q_pbc.gro -cs spc216.gro -o 1e0q_solv.gro -p 1e0q.top

Check 1e0q.top, [molecules] include 1 Protein and 6694 waters. (Note this only apply to water solvent, if solvate with non-water molecules, the topology file won't be updated)

Since net charge of this molecule is 0, we do not need to add ions to neutralize it, but it does not hurt to add ions as a control of ionic strength of the solvent.

gmx_mpi grompp -f ions.mdp -c 1e0q_solv.gro -p 1e0q.top -o ions.tpr

gmx_mpi genion -s ions.tpr -o 1e0q_solv_ions.gro -p 1e0q.top -pname NA -nname CL -conc 0.1

when prompted, choose group 13 SOL. Check 1e0q.top, [molecules] include 1 Protein 6670 waters and 12 sodium ions and 12 chloride ions.

Now we perform energy minimization to remove steric clashes or inappropriate geometry.

gmx_mpi grompp -f minim.mdp -c 1e0q_solv_ions.gro -p 1e0q.top -o em.tpr
gmx_mpi mdrun -v -deffnm em

 gmx_mpi energy -f em.edr -o potential.xvg

1e0q_em_potential

Two-step equilibration procedure - Now we need to relax the solvent and ions while keeping the protein atom positions restrained. This will be accomplished in two phases: 100 ps NVT followed by 100 ps NPT ensembles.

gmx_mpi grompp -f nvt.mdp -c em.gro -p 1e0q.top -o nvt.tpr

gmx_mpi mdrun -deffnm nvt

gmx_mpi energy -f nvt.edr -o temperature.xvg

1e0q_nvt_temperature

gmx_mpi grompp -f npt.mdp -c nvt.gro -t nvt.cpt -p 1e0q.top -o npt.tpr

gmx_mpi mdrun -deffnm npt

gmx_mpi energy -f npt.edr -o pressure.xvg

1e0q_npt_pressure

gmx_mpi energy -f npt.edr -o density.xvg 

1e0q_npt_density

Now the system is well equilibrated with desired temperature and pressure, next step is to run production MD without restraints on heavy atoms.

gmx_mpi grompp -f md.mdp -c npt.gro -t npt.cpt -p 1e0q.top -o md_0_1.tpr

gmx_mpi mdrun -deffnm md_0_1

After the simulation finished, if you load the trajectory file into vmd, we will see the jump of atoms across the boundaries, now we use trjconv to remove such jump.

gmx_mpi trjconv -s md_0_1.tpr -f md_0_1.xtc -o md_0_1_noPBC.xtc -pbc mol -ur compact

Select 0 system.Now we load the trajectory into vmd

vmd md_0_1.gro md_0_1_noPBC.xtc

All analysis are based on this corrected trajectories.

gmx_mpi rms -s md_0_1.tpr -f md_0_1_noPBC.xtc -o rmsd.xvg -tu ns

1e0q_md_rmsd

This is the rmsd relative to minimized equilibrated structure.

gmx_mpi gyrate -s md_0_1.tpr -f md_0_1_noPBC.xtc -o gyrate.xvg

1e0q_md_gyrate