Replica Exchange Molecular Dynamics Simulation for small peptides/Globular Proteins with Gromacs-4.5 and above

In this method for the peptide or proteins, a number of simulations are performed at different temperatures in parallel, and exchanges of configurations [Amino acids Corodinates in Gromacs] are tried periodically. Even if a trajectory is temporarily trapped in a local minimum, the simulation can escape from this minimum via an exchange with a higher temperature configuration. With this method, one can obtain various thermodynamic quantities as a function of temperature for a wide temperature range from a single simulation run. Moreover, because each replica can be simulated using its own computer processor, the REMD method is well suited for and very efficiently runs on parallel computers.

Step1: pdb2gmx –ignh –f xx.gro –o xx_pro.gro –water spce –ter

Choose 11. popc_ffG45a3 force field

Step2: editconf –f xx_pro.gro –o xx_newbox.gro –c –d 1.0 –bt cubic

Here the box size is 1.0 nm and the box is cubic

Step3: genbox –cp xx_newbox.gro –cs spc216.gro –o xx_solv.gro –p topol.top

Step4: grompp –f ions.mdp –c xx_solv.gro –p topol.top –o ions.tpr –maxwarn 100

Step5: genion –s ions.tpr –o xx_solv_ions.gro –p topol.top –pname NA –nname CL –np x –nn x

Add any number of NA or CL atoms to have equal charge in the system, x indicates the no. of atoms

Step6: grompp –f minim.mdp –c xx_solv_ions.gro –p topol.top –o em.tpr –maxwarn 100

Step7: mdrun –nt 8 –v –deffnm em

Where 8 is the no. of processsors in local machine

Step8: grompp –f nvt.mdp –c em.gro –p topol.top –o nvt.tpr –maxwarn 100

Step9: mdrun –deffnm nvt

Step10: Replica Set-Up

Now Count the no. of water molecules and No. of proteins atoms in the system

Now use the following link to generate the temperature generator for REMD

http://folding.bmc.uu.se/remd/

Keep Exchange Probability: 0.2 (which is supposed to give the best results, refer to the literatures for more information)

Select the Minimum and maximum temperature range for the particular system

For Example:

No. of water molecules: 63012

No. of protein atoms: 186

Low temperature : 298.15K = 25 °C [Room Temperature]

High Temperature: 310.15K = 37 °C [Human Body Temperature]

This will generate 16 different temperatures:

298.15, 298.17, 299.80, 300.62, 301.45, 302.28, 303.11, 303.95, 304.78, 305.61, 306.45, 307.29, 308.13, 308.97, 309.82, 310.66

[Note: Higher temperature should be chosen so that the particular peptide or protein would cross the energy barrier level and will fold properly]

So, now create 16 different folders (or depends on the how many of temperatures we get for the particular system)

For example give the folder name equil and this folder keep the following files: xx.top , npt.gro , equil.mdp

Now will have 16 different folder with name equil and each folder will have 3 above described files and now change the temperature of all 16 equil.mdp’s from different folders

Step11: Equilibration Run

We will run the equilibration run in all different folders altogether

(for dir in equil{1..15}; do cd $dir; grompp –f equil –c npt –maxwarn 200; cd..; done)

mpirun –np 16 mdrun –v –multidir equil{1..15}

It will run equilibration in all different 16 folders

Gromacs version should be MPI enabled. Refer to Gromacs website for more information.

Step12: MD Run

After Equilibration Run make a folder with name sim and and that folder keep the following 4 files:

1. confout.gro (this file will generate after finishing the equilibration run in every equil folders)

2. sim.mdp (can be downloaded from gromacs website)

3. state.cpt

4. topol.top

[Now change all of the file name as same like xx.mdp, xx.gro and xx.top or as you wish]

[Now write the all different 16 temperature in xx.mdp]

Now can type the following commands:

(for dir in sim{1..15}; do cd $dir; grompp –f xx –c xx –t state –maxwarn 200; cd..; done)

mpirun –np 16 mdrun –v –multidir sim{1..15} –s xx –c xx –o xx –g xx –multi 16 –replex 1000 –reseed 173529

Step13: Replica exchange run in Supercomputers {The code might be different from the below}

1. In your working directory

source /pkg/local/gcc/4.6.1/gcc461.sh

2. Prepare the *.tpr file with grompp program

/pkg/chem/gromacs46b1/bin/grompp -f *.mdp -c *.gro -p *.top -o *.tpr -maxwarn 200

where files *.mdp *,gro *.top *.tpr should be replaced with the real file names respectively

3. Prepare the job script file jobname.gmx like this

#!/bin/bash

#BSUB -J jobname

#BSUB -e jobname.err

#BSUB -o jobname.stderr

#BSUB -R 'span[ptile=16]'

#BSUB -q 48cpu

#BSUB -n 16

#BSUB -R "ipathavail==0"

source /pkg/chem/gromacs46b1/setgmx

source /pkg/chem/scripts/setgccmvapich2

cat $LSB_DJOB_HOSTFILE > ./host.$$

echo "Your GROMACS job starts at `date` "

mpirun_rsh -np 16 -hostfile ./host.$$ /pkg/chem/gromacs46b1/bin/mdrun_mpi -s jobname.tpr -multi 16 -replex 1000 -reseed 175329 -o jobname.trr -c jobname_after_md -v -g jobname.log

4. Submit your job

bsub < jobname.gmx

In this job script file jobname.gmx, 16 cores are used, 16 replicas, (i.e., each replica use 1 Core, here multiple no. of cores can also be given like 16*2 = 32, 16*3 = 48)

48cpu is the name of the queue, input files for mdrun_mpi have the filename root as jobname.