Molecular dynamics of an aqueous system with GROMACS

Preparation of a water box

Example of how to simulate a box of 2 nm full of TIP3P water.

To run a simulation you need a coordinate file and a topology file, and the file with the simulation parameters

To generate the coordinate file (solvated.gro)

gmx solvate -box 2 -cs spc216.gro -o solvated.gro

The topology file (topol.top)

#include "amber03.ff/forcefield.itp"
#include "amber03.ff/tip3p.itp"
#include "amber03.ff/ions.itp"

[system]
water

[molecules]
SOL         [N]

Where N is the number of water molecules in solvated.gro

Preparion of a solvated macromolecule

You need the pdb coordinates of the molecule, with a format which GROMACS can understand

gmx pdb2gmx -f molecule.pdb

generates both the coordinate conf.gro file and the topology file topol.top. During the following steps, further information will be added at the topology file by GROMACS. The last lines of the topology file show the number of solvent molecules and ions.

First you have to specify the box size editing the last line of conf.gro The size units are nm.

Then you add water molecules to your simulation box:

gmx solvate -cp conf.gro -cs spc216.gro -p topol.top -o solvated.gro

And now you have your system ready for minimization

Adding ions to the simulation box

Sometimes the macromolecule is charged. Then you have to add ions to the simulation box in order that the net charge of the system is 0.

Suppose we want to simulate an acidic molecule (negative charge). To make any operation on GROMACS you first have to generate the tpr file:

gmx grompp -f ions.mdp -c solvated.gro -p topol.top -o ions.tpr

Now you can add the ions:

gmx genion -s ions.tpr -o solvated.gro -p topol.top -pname NA -np 1

gmx genion -s ions.tpr -o solvated.gro -p topol.top -pname NA -nname CL -np 18 -nn 6

gmx genion -s ions.tpr -o solvated.gro -p topol.top -pname NA -nname CL -neutral

genion adds the ions to the coordinate and topology files.

You can find an example here

Energy minimization

Once you have the coordinates and the topology, you can prepare the simulation file with the simulation parameters em.mdp

gmx grompp -f em.mdp -c solvated.gro -p topol.top -o em

And this generates the tpr (simulation) file em.tpr

Now you can run the minimization:

gmx mdrun -deffnm em

This generates em.gro, the coordinates of the minimized system.

More information at the page of Luis Carlos Pardo

After minimization, you make the equilibration (keeps fixed the coordinates of the macromolecule while letting move - rearrange the water molecules)

gmx grompp -f pr.mdp -c em.gro -r em.gro  -p topol.top -o pr

gmx mdrun -deffm pr

MD simulation

After equilibration you can finally make the long MD simulation, starting from pr.gro, the coordinates of the equilibrated system

gmx grompp -maxwarn 2 -f run.mdp -c pr.gro -p topol.top -o run

(the optional maxwarn is to ignore the usual warnings)

gmx mdrun -deffnm run

This simulation can be run in any computer, it only needs the tpr file as input. Usually you generate the tpr file of a long simulation in your computer, and then you run the long simulation in a cluster.

The trajectory files are run.trr or run.xtc. From them you can generate the frames of the trajectory in PDB format

gmx trjconv -s run.tpr -f run.trr -o traj.pdb  -pbc nojump -b 100 -e 200

In this example we generated the movie of the trajectory from t=100 ps till t=200 ps, preventing the breaking the macromolecule at the box wall (-pbc nojump)

You can find an example of the parameter for minimization em.mdp, for equilibration pr.mdp and for simulation run.mdp in this document