Gromacs with PLUMED Demo
Prepare the molecule for MD
- Acquire a PDB file frome the protein database.
- Load gromacs/plumed module on ACISS
>>> module load gromacs/4.6.1p
>>> pdb2gmx_mpi_d -f .pdb -o _processed.gro -water spce
- Choose a force-field (usually 5)
Select the Force Field:
From '/usr/local/gromacs/share/gromacs/top':
1: AMBER03 force field (Duan et al., J. Comp. Chem. 24, 1999-2012, 2003)
2: AMBER94 force field (Cornell et al., JACS 117, 5179-5197, 1995)
3: AMBER96 force field (Kollman et al., Acc. Chem. Res. 29, 461-469, 1996)
4: AMBER99 force field (Wang et al., J. Comp. Chem. 21, 1049-1074, 2000)
5: AMBER99SB force field (Hornak et al., Proteins 65, 712-725, 2006)
6: AMBER99SB-ILDN force field (Lindorff-Larsen et al., Proteins 78, 1950-58, 2010)
7: AMBERGS force field (Garcia & Sanbonmatsu, PNAS 99, 2782-2787, 2002)
8: CHARMM27 all-atom force field (with CMAP) - version 2.0
9: GROMOS96 43a1 force field
10: GROMOS96 43a2 force field (improved alkane dihedrals)
11: GROMOS96 45a3 force field (Schuler JCC 2001 22 1205)
12: GROMOS96 53a5 force field (JCC 2004 vol 25 pag 1656)
13: GROMOS96 53a6 force field (JCC 2004 vol 25 pag 1656)
14: OPLS-AA/L all-atom force field (2001 aminoacid dihedrals)
15: [DEPRECATED] Encad all-atom force field, using full solvent charges
16: [DEPRECATED] Encad all-atom force field, using scaled-down vacuum charges
17: [DEPRECATED] Gromacs force field (see manual)
18: [DEPRECATED] Gromacs force field with hydrogens for NMR
Solvation
Two steps to defining a box and filling it with solvent.
- Define the box dimensions.
- Fill the box with water.
>>> editconf_mpi_d -f _processed.gro -o _newbox.gro -c -d 1.0 -bt cubic
The above command centers the protein in the box (-c), and places it at least 1.0 nm from the box edge (-d 1.0). The box type is defined as a cube (-bt cubic).
>>> genbox_mpi_d -cp _newbox.gro -cs spc216.gro -o _solv.gro -p topol.top
- The configuration of the protein (-cp) is contained in the output of the previous editconf step, and the configuration of the solvent (-cs) is part of the standard GROMACS installation.
- The output is called < name >_solv.gro, and we tell genbox the name of the topology file (topol.top) so it can be modified.
Add ions
We now have a solvated system that contains a charged protein. Look at the last line of your [ atoms ] directive in topol.top; it should read (in part) "qtot #" (# stands for the total charge).
grompp processes the coordinate file and topology (which describes the molecules) to generate an atomic-level input (.tpr). The .tpr file contains all the parameters for all of the atoms in the system.
To produce a .tpr file with grompp, we will need an additional input file, with the extension .mdp (molecular dynamics parameter file); grompp will assemble the parameters specified in the .mdp file with the coordinates and topology information to generate a .tpr file.
(An example .mdp file can be downloaded here)
>>> grompp_mpi_d -f ions.mdp -c _solv.gro -p topol.top -o ions.tpr
Generate ions in our box.
When prompted, choose group 13 "SOL" for embedding ions. You do not want to replace parts of your protein with ions.
>>> genion_mpi_d -s ions.tpr -o _solv_ions.gro -p topol.top -pname NA -nname CL -nn 8
In the genion command, we provide the structure/state file (-s) as input, generate a .gro file as output (-o), process the topology (-p) to reflect the removal of water molecules and addition of ions, define positive and negative ion names (-pname and -nname, respectively), and tell genion to add only the ions necessary to neutralize the net charge on the protein by adding the correct number of negative ions (-nn 8).
Energy Minimization
Use an energy minimization .mdp file for this pre-processing. (found here)
>>> grompp_mpi_d -f minim.mdp -c _solv_ions.gro -p topol.top -o em.tpr
Now run the energy minimization.
>>> mdrun_mpi_d -v -deffnm em
The -v flag is for the impatient: it makes mdrun verbose, such that it prints its progress to the screen at every step. The -deffnm flag will define the file names of the input and output. So, if you did not name your grompp output "em.tpr," you will have to explicitly specify its name with the mdrun -s flag. In our case, we will get the following files:
- em.log: ASCII-text log file of the EM process
- em.edr: Binary energy file
- em.trr: Binary full-precision trajectory
- em.gro: Energy-minimized structure
NOTE: To check if the energy is minimized
Type 10 0 at the prompt after entering the following on the command line. Then plot the xvg file.
>>> g_energy_d -f em.edr -o potential.xvg
Equilibration
We need to make sure the structure/geometry is reasonable. Use the .itp file from above to make sure the protein's structure doesn't change.
Two phases to equilibrate the system:
1. NVT equilibration (isothermal-isocharic or canonical) to bring temperature to logical choice
Use the .mdp file here
>>> grompp_d -f nvt.mdp -c em.gro -p topol.top -o nvt.tpr
>>> mdrun_d -deffnm nvt
Check the temperature for equilibration
Type 15 0 and plot file
>>> g_energy_d -f nvt.edr
2. NPT equilibration ("isothermal-isobaric")
Use the .mdp file here
>>> grompp_d -f npt.mdp -c nvt.gro -t nvt.cpt -p topol.top -o npt.tpr
>>> mdrun_d -deffnm npt
Analyze the pressure
Type 16 0 and plot
>>> g_energy_d -f npt.edr -o pressure.xvg
Useful commands in GROMACS
Use the following command for checking out the RMSD of the system. This outputs a file called rmsd.xvg that can be viewed on xmgrace
>>> g_rms_mpi_d -s .tpr -f full .xtc