Skip to content

MD with GROMACS for SMALL molecules

I am interested in MD simulation of "small" molecules, and I found little about doing that with GROMACS (and a lot to deal with huge biological molecules). Here I try to describe a little bit the main steps to do that. I wrote this page for myself, but since it might be help to someone I decided to published it. This is not a tutorial and also does not deal with the basics of MD simulations. Is a cookbook where I gathered webpages and scripts since I keep forgetting where to find things ;-). You can follow the following links if you want to do some of the following:

 


Installing and simulating with GROMACS

First of all, it would be a good ide to download GROMACS and INSTALL it. You might also install it using

  • "sudo get-apt install gromacs"

Is easier but you cannot control it!

Input files needed:

See the GROMACS manual for a description of all files

  • A file with the molecular structure such as a PDB, mol2, or gro file.
  • A file with the force-field: the *.top in GROMACS can contain all force fields needed for the simulation (VW parameters, dihedrals and so on). But it can also contain a refererence to a *.itp which contains the force-field for a particular molecule. This is to recycle the ForceFiles that you have already generated for another molecules
  • A file that tells GROMACS what to do, i.e. number of steps, wether to use an NVE, NVT, NPT ensemble and so on. This is *.mdp file.

All these files are summarized in a *.tpr file (usually called topol.tpr) builded using the grompp program. This file containts EVERYTHING needed to start the simulations. Unfortunatelly this file is not intended for human eyes!!!

Goal: get a force field to perform an MD simulation:

  1. Search on internet a PDB file. If you don't find one, you might want to generate a raw PDB using, for example, avogadro. Once you have this you have different forcefields to do the MD simulation
  2. Generate a topology file using a force-field:

 

CHARMM

An excellent option is to use the automatic generator from SWISSPARAM.

Another option is to use CGENF, but then you will have to use this   python script to translate the output to a GROMACS compatibel topology file

AMBER (GAFF)

You might use antechamber and tleap to prepare your simulations and then this program to translate the amber file to gromacs format. For conversion between GAFF and GROMACS see below.

GROMOS

You might want to search in the ATB data base  your molecule or build one with the tools in this web page

You can also use the PRODGR program to generate your topology file

OPLS

You can usethe TPP software. This web page will explain you how to get an OPLS force field using this program. Morover, you will have to add a residue to your GROMACS distribution. Information about how to do that can be found here. You can also click here to read the recipe to do that have generated with the help of Nebojsa Zec.

 

From AMBER to GROMACS

Here you can find an explanation to translate a force field from amber to GAFF (amber) to gromacs format, summarizing.

For non-bonded:

  • εgroGAFF·4.184
  • σgro=RGAFF·25/6  * 1/10 = 0.17818·RGAFF

For bonded:

  • Bonds: b0gro = b0GAFF / 10   ;    kb,gro = kb,GAFF /10 *2 *4.184
  • Angles: Egro = Emin/2 cth = EGAFF *2 * 4.184
  • Dihedrals: kdgro = PKGAFF * 4.184 / IDIVFGAFF

 

Goal: Produce a stable molecule to simulate

 

 Well, you should now have a PDB and an ITP. We will now put one molecule in a big box to be sure that the molecule is stable. If it expoltoes it is not worth to go any further. We will thus generate a *.gro file which is the kind of molecular files to feed GROMACS. The syntaxis is:

  • gmx editconf -f molecule.pdb -o molecule.gro -box Lx Ly Lz

where molecule.pdb is the input molecule and molecule_eq.gro is the equilibrated molecule. Lx, Ly and Lz are the lengths of the box. Now we should have the molecule in a box (molecule.gro) a force field file (molecule.top) and we have to add a file that tells gromacs what to do (an mdp file). You can find one here but be sure to change the integrator md (used to do an MD simulation) to steep to minimize the energy.We will call this file min.mdp Then with:

  • gmx grompp -f min.mdp -c molecule.gro -p molecule.top

The files should be ready to use. To actually do the simulation use

  • gmx mdrun

And everything should work...

Goal: to generate a box with the right number of molecules

For the case of a crystal structure click here

Now we will produce the actual box to simulate. Fist of all you should have to decide a size of your box that is compatible with the density of your system. The problem will be to fit the molecules in the box. For this reason is better to start with a big box (a lower density) and then increase the density at small steps. Let's assume that you want to simulate 125 molecules in a box of 50 AA with a density of 1000g/l. Be aware that unit in GROMACS are Nm! First we will change the size of the box with one molecule so that it is big enough to fit all molecules you want to add. "Big enough" means that the program does not complains because he couldn't put all molecules into the box.:

  • gmx editconf -f mol_eq.gro -o mol_eqb.gro -box 5 5 5

We will then add 124 molecules to the box with:

  • gmx insert-molecules -f mol_eqb.gro -ci mol_eqb.gro -o box.gro -nmol 124 -rot xyz -try 100

Since we want to insert the same molecule (mol_eqb.gro) to the existing file mol_eqb.gro the configurations to be inserted and inserted are the same. We will also insert molecules randomly rotated with -rot, and will try 100 times to insert the molecules. After that you minimize the energy... if it blows up, you should start with an even smaller density.

  • gmx grompp -f min.mdp -c box_noeq.gro -p topol.top
  • gmx mdrun
  • mv confout.gro box_eq.gro

I added a last line to change the name of the output *.gro file so that it says it is the equilibrated box.

O.K., let's assume you were able to minimize. Now you increase the density and minimize again.... and repeat until you arrive to the desired density.

  • gmx editconf -f box_eq.gro -density 750 -o box_noeq.gro
  • gmx grompp -f min.mdp -c box_noeq.gro -p topol.top 
  • gmx mdrun
  • mv confout.gro box_eq.gro

You must repeat these lines until you have converged... and that's all (if you manage at the very first time, the you are really very lucky or very good!!!).

It is maybe better to do this in an automatic way using a bash file like this one.

By now you shoud have produced:

  • A confout.gro file with the final configuration
  • A traj.trr with the trajectory of your molecules
  • An ener.edr with the energy

Let's have a look at these files:

ANALYSIS OF GROMACS RESULTS