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
• Analyze some results
• Generating scattering functions with Sassena
• Generating a F(q) with my code (through g(r))
• Crystals are a bit special. Some things usefull to deal with them:
• Going from a CIF file to a PDB file
• Calculating Bragg peaks from a PDB
• Disordered crystals
• Units and summary of Force Fields (a nightmare)
• Some LINUX command I keep forgetting

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:

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:
   

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:

1. Add the position restraint to your "itp" file
2. If you want to use "IFDEF" change your topology file
3. Tell gromacs the file you use for position restraint when "grompp"-ing the MD simulation
4. 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.