Running a simulation in Gromacs

These brief instructions are a general guide to how to run a simple MD simulation using Gromacs. More detailed tutorials for a variety of systems can be found elsewhere. In this case we are simulating the alanine dipeptide using Gromacs 2018 but things should not change too much for other recent versions of the software. You must have a working version installed in your machine.

Get the files

First of all, you will have to download the compressed files in alaTB_files.tar.gz to your computer and extract them typing

tar -xvf ala_TB_files.tar.gz

in your terminal. When you do this, you should be able to find a file called alaTB.gro together with a number of mdp files that we will later use.

Generate the topology

The next step is generating a topology file from the configuration file alaTB.gro. We will do this using the pdb2gmx program. Before we do this we will generate variables in bash with informative names which will be useful for sensible file name structure. This is what we will run in our terminal

proot="alaTB"
ff="ff03"
wat="tip3p"

out="${proot}_${ff}_${wat}"
top="${proot}_${ff}_${wat}"
gmx pdb2gmx -f $proot -o $out -p $out <<EOF
1
1
EOF

Using 1 and 1 as options after the EOF we are choosing specific options for our force field and water model (Amber03 and TIP3P respectively). These will suffice in this illustrative example but it would be advisable to do some research on your system for an educated choice of force field and water model combination.

Editing the box and solvation

Next we will change the size of our simulation box and we will fill it with water molecules. For this we will use two additional Gromacs programs:

inp=$out
out="${proot}_${ff}_${wat}_edit"
gmx editconf -f $inp -o $out -c -d 1 -bt cubic

inp=$out
out="${proot}_${ff}_${wat}_solv"
gmx solvate -p $top -cp $inp -o $out

The first command edits the box so that the molecule is centered (-c) in the cubic (-bt cubic) simulation box and it has 1 nm between molecule and box sides (-d 1). The second command simply solvates the molecule filling the box with the water model geometry of choice.

Energy minimization

onsidering that we are starting from a peptide conformation and a fixed water box configuration that are not particularly meaningful, we will first energy minimize the system using steepest descent.

mdp="mini.mdp"
inp="${proot}_${ff}_${wat}_solv"
out="${proot}_${ff}_${wat}_mini"
gmx grompp -p $top -c $inp -f $mdp -o $out

inp=$out
out=$inp
gmx mdrun -v -s $inp -deffnm $out

As we will do many times later, we first use the Gromacs preprocessor (grompp) and then run the calculation (mdrun). Since we are using the verbose (-v) option you will find that this produces a lot of output, corresponding to the gradual decrease in the energy of the system.

MD in the canonical ensemble with position restraints

We then equilibrate the water molecules while keeping the peptide mostly fixed using position restraints on its atoms. This is not too important in the context of a small peptide like the one we are considering, but can be critical for larger systems.

mdp="nvt_posre.mdp"
inp=$out
out="${proot}_${ff}_${wat}_nvt_posre"
gmx grompp -f $mdp -c $inp -p $top -r $inp -o $out

inp=$out
out=$inp
gmx mdrun -v -s $inp -deffnm $out

MD in the isothermal-isobaric ensemble

The next step is running MD and fixing the pressure to the value of interest (typically 1 bar) with the help of a barostat. Most of the remaining parameters stay the same as were in the NVT simulation.

mdp="npt.mdp"
inp=$out
out="${proot}_${ff}_${wat}_npt"
gmx grompp -f $mdp -c $inp -p $top -r $inp -o $out -t ${inp}.cpt

inp=$out
out=$inp
gmx mdrun -v -s $inp -deffnm $out

Production run

Finally we run a production run, again in the canonical ensemble, but now removing restraints on the peptide so that we can adequately sample conformational space.

mdp="sd_nvt.mdp"
inp=$out
out="${proot}_${ff}_${wat}_nvt"
gmx grompp -f $mdp -c $inp -p $top -r $inp -o $out -t ${inp}.cpt

inp=$out
out=$inp
gmx mdrun -v -s $inp -deffnm $out