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 publish 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 use these links if you want to do some of the following:

Install and run GROMACS (stay in this page)

Position restraints

Installing and simulating with GROMACS

First of all, it would be a good idea to download GROMACS, compile it, 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:

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
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. Another interesting option is to use the web http://zarbi.chem.yale.edu/ligpargen/ . Thanks to Yitian to make me aware of this!!

Of course more information and more reliable is found in this GROMACS web page.

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:

ε_gro=ε_GAFF·4.184
σ_gro=R_GAFF·2^5/6 * 1/10 = 0.17818·R_GAFF

For bonded:

Bonds: b_0gro = b_0GAFF / 10 ; k_b,gro = k_b,GAFF /10 *2 *4.184
Angles: E_gro = E_min/2 cth = E_GAFF *2 * 4.184
Dihedrals: kd_gro = PK_GAFF * 4.184 / IDIVF_GAFF

You will find more useful information about force fields and units used in MD by GROMACS here.

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

Position Restraints

Sometimes you might want to keep atoms of molecules more or less fixed during the equilibration. Maybe you want to do this to first "relax" molecule structure and then the global structure, or maybe you want to randomize the system (a crystal for example) going to high temperature... butyou don't want the centre of mass of the molecules to leave the crystalline structure. There is an option to do that in GROMACS. We proudly present: the Position Restraint

Again, look at the manual for a long and complete description, here you will find the recipe. You must do four things:

Add the position restraint to your "itp" file
If you want to use "IFDEF" change your topology file
Tell gromacs the file you use for position restraint when "grompp"-ing the MD simulation
Generate the file with the positons you want to restraint to. (can be the same as the input conf. file!)

1.- Add position restraint

Go to your ".itp" file and add, after section atoms (one per molecule):

#ifdef POSRES
[ position_restraints ]
; ai func g r K
1 2 2 0.15 10
.

.
#endif

This will activate position restraints when told by *.mdp. Position Restraint will "fix" atom 1, to the position of atom 1 in the position restraint configuration file. for a description of the options go to the manual... it is soooo well explained

2.- Change your topology file

Now you must add this line to your topology file if you want to actually do a position restraint:

#define POSRES

... and you can switch it off by adding a ; in front of the line: ; #define POSRES

3.- Tell gromacs you want to positionrestrain your simulation

Just write the usual command in grompp... but adding the file with the coordinates you want to be retrained to:

gmx grompp -f min.mdp -c configuration.gro -p topol.top -r configuration_r.gro

Where configuration_r.gro (file with the coordinates you want to restraint to) can be the same (I mean with teh same name) as the configuration file

4.- Generating coordinate file for restraints

Well, usually you don't need to do that since you want to restraint the positions to the initial configration. If you want to do some fancy stuff (like restraining to the centre os mass) you might want to use my lovely, sweety, tiny, winny traductor software.