Gromacs Quick Start

GROMACS Quick Start

Official Documentation:

https://manual.gromacs.org/current/user-guide/index.html

Gromacs Server and Environment:

Gromacs is pre-installed on a Levich server and user accounts are initilized to include Gromacs path.

In case you want to see how the Gromacs environment was set, Here it is:

In file ~/.profile, the last line points to gromacs environment initialization:

source /usr/local/gromacs/bin/GMXRC

Users are expected to install their own Python and Python modules.

Install Conda Python and Python modules:

wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh

chmod +x Miniconda3-latest-Linux-x86_64.sh

./Miniconda3-latest-Linux-x86_64.sh (answer yes to initialize Minicoda3)

Logout and back in to have conda environment, then run:

conda install numpy scipy matplotlib pandas

Conda can be deactivated or reactivated by command:

conda deactivate

conda activate

Gromacs Command

Main Command: gmk

Examples:

gmx -h (print help)

gmk -v (show version)

gmx help [module] (documentation of a module)

Example 1: Lysozyme

In this example, we are going to use hen egg white lysozyme (PDB ID 1AKI) to run a MD simulation.

Prepare structure and topology files

Create a working directory:

mkdir lysozyme

cd lysozyme

Download a pdb file from Protein Data Bank

getpdb 1aki

Pymol molecule structure viewer is preinstalled on the system. Due the Python dependency, one needs to deactivate conda to run Pymol.

conda deactivate

pymol 1aki.pdb

After verifying structure, quit pymol and activate conda again.

conda activate

Delete crystal water molecules. Water in the pdb file are crystal water molecules. Gromacs has its own way to add solvent water.

grep -v HOH 1aki.pdb > 1aki_clean.pdb

Create topology file and process the structure file

gmx pdb2gmx -f 1aki_clean.pdb -o 1AKI_processed.gro -water spce

When it will prompts to select a force field, select choice 15 - OPLS force field.

The explanation of pdb2gmx module can be found by running:

gmx help pdb2gmx

-f 1aki_clean.pdb: Input file

-o 1AKI_processed.gro: Output file for Gromacs use

-water spce: The model to build solvent water

Other two files are

topol.top Topology file defines atom and bond parameters.
posre.itp Position restraint file

Define the simulation box and solvent

Step 1: Define the box dimensions using the editconf module

gmx editconf -f 1AKI_processed.gro -o 1AKI_newbox.gro -c -d 1.0 -bt cubic

This command use module editconf to center and define a box.

-f 1AKI_processed.gro: Input file
-o 1AKI_newbox.gro: Output file, the box is 0,0,0 and the coordinates at the bottom line
-c: center the box
-d 1.0: Leave at least 1 nm at the edge
-bt cubic: Use cubic box. There are other choices such as rhombic dodecahedron.

Step 2: Fill the box with water using the solvate module

gmx solvate -cp 1AKI_newbox.gro -cs spc216.gro -o 1AKI_solv.gro -p topol.top

This command uses module solvate to add water molecules into the box.

-cp 1AKI_newbox.pro: Configuration of the protein from the named file
-cs spc216.gro: Configuration of the solvent. Spc216.gro is a generic equilibrated 3-point solvent model good for SPC, SPC/E, or TIP3P water.
-o 1AKI_solv.gro: Output file name
-p topol.top: Topology file name. Solvate module will update this file to include both protein molecule and solvate (SOL) line

Add ions

In topology file [ atom ] section, the protein total charge is calculated. The charge at the end of this section is the net charge.

1960   opls_272    129    LEU     O2    682       -0.8    15.9994   ; qtot 8

In this example, the net charge is 8.

In MD simulation, we need to balance the charge with ions so that we have a neutral system. This is a two step procedure.

Step 1: Prepare a run input file (extension .tpr) for genion module

MD parameter file ions.mdp contains instructions for Gromacs Preprocessor module grompp to assemble coordinates and topology into an atomic-level input .tpr file.

Sample ions.mdp file

; ions.mdp - used as input into grompp to generate ions.tpr
; Parameters describing what to do, when to stop and what to save
integrator  = steep         ; Algorithm (steep = steepest descent minimization)
emtol       = 1000.0        ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep      = 0.01          ; Minimization step size
nsteps      = 50000         ; Maximum number of (minimization) steps to perform

; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist         = 1         ; Frequency to update the neighbor list and long range forces
cutoff-scheme	= Verlet    ; Buffered neighbor searching 
ns_type         = grid      ; Method to determine neighbor list (simple, grid)
coulombtype     = cutoff    ; Treatment of long range electrostatic interactions
rcoulomb        = 1.0       ; Short-range electrostatic cut-off
rvdw            = 1.0       ; Short-range Van der Waals cut-off
pbc             = xyz       ; Periodic Boundary Conditions in all 3 dimensions

This mdp file tells Gromacs to run an energy minimization.

gmx grompp -f ions.mdp -c 1AKI_solv.gro -p topol.top -o ions.tpr

Module grompp is a gromacs preprocessor. Its job is to make a .tpr file.
-f ions.mdp: Read instruction from ions.mdp file.
-c 1AKI_solv.gro: Coordinates file
-p topol.top: Topology file
-o ions.tpr: Output file. This is an atomic level input file with coordinates and topology all assembled. It's going to be the input of MD simulation. In this case, it will be the input of ion adding module's input file.

Step 2: Use module genion to replace some water molecules with ions

gmx genion -s ions.tpr -o 1AKI_solv_ions.gro -p topol.top -pname NA -nname CL -neutral

-s ions.tpr: Specify structure file ions.tpr.
-o 1AKI_solv_ions.gro: Write to this output file.
-p topol.top: Update topology file to reflect the removal of water and addition of ions.
-pname NA: Use NA for position ion.
-nname CL: Use CL for negative ion.
-neutral: Neutralize the system. In this case, it will replace 8 waters by CL- to offset the 8 positive net charge.

When prompted, choose option 13 SOL so module genion will replace solvent molecules.

After this step, the topol.top file will include CL in its [ molecules ] section:

[ molecules ]
; Compound #mols
Protein_chain_A 1
SOL 10636
CL 8

Energy minimization:

We have a solvated and charge neutral system by now in coordinates file 1AKI_solv_ions.gro with the molecule toplogy in file topol.top. Before we run production MD, we have a few more prepartion steps:

Energy minimization: Remove structure clashes.
Equilibration: Move the structure from high energy state to equilibrated state.

Similar to add ions step, we need to make a MD include-all atomic level .tpr file, with 3 pieces of information:

.mdp file: MD parameter file that serves as instruction script
.gro file: Gromacs coordinate file
.top file: Topology file

Sample minim.mdp file:

; minim.mdp - used as input into grompp to generate em.tpr
; Parameters describing what to do, when to stop and what to save
integrator  = steep         ; Algorithm (steep = steepest descent minimization)
emtol       = 1000.0        ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep      = 0.01          ; Minimization step size
nsteps      = 50000         ; Maximum number of (minimization) steps to perform

; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist         = 1         ; Frequency to update the neighbor list and long range forces
cutoff-scheme   = Verlet    ; Buffered neighbor searching
ns_type         = grid      ; Method to determine neighbor list (simple, grid)
coulombtype     = PME       ; Treatment of long range electrostatic interactions
rcoulomb        = 1.0       ; Short-range electrostatic cut-off
rvdw            = 1.0       ; Short-range Van der Waals cut-off
pbc             = xyz       ; Periodic Boundary Conditions in all 3 dimensions

This .mdp file is the same as ions.mdp except the coulombtype. PME is Fast smooth Particle-Mesh Ewald (SPME) electrostatics, more accurate than cutoff in ions.mdp.

gmx grompp -f minim.mdp -c 1AKI_solv_ions.gro -p topol.top -o em.tpr

grompp: Gromacs preprocessor to assemble .gro and .top files.
-f minim.mdp: MD parameter file
-c 1AKI_solv_ions.gro: Coordinate file
-p topol.top: Toplogy file
-o em.tpr: Output file

Once the em.tpr is ready, we can pass it to mdrun module to run an energy minimization.

gmx mdrun -v -deffnm em

mdrun: MD simulation module
-v: Verbose mode
-deffnm: Define file names of the input and output. If you did not name your grompp output "em.tpr," you will have to explicitly specify its name with the mdrun -s flag.

Since we passed the model em in, mdrun takes em.tpr as input and writes out 4 files that start with em. Gromacs will detect CPU and uses OpenMP to run on maximum available threads. If you would like control the number of threads manually, use option

-ntomp

For example

gmx mdrun -v -ntomp 8 -deffnm em

Will uses 8 threads.

This step produces output files as following:

(base) jmao@gromacs:~/demo/lysozyme$ ls -lt
total 9044
-rw-rw-r-- 1 jmao jmao 305030 Oct 18 13:32 em.log
-rw-rw-r-- 1 jmao jmao 1524475 Oct 18 13:32 em.gro
-rw-rw-r-- 1 jmao jmao 406632 Oct 18 13:32 em.trr
-rw-rw-r-- 1 jmao jmao 129712 Oct 18 13:32 em.edr
-rw-rw-r-- 1 jmao jmao 848248 Oct 18 13:27 em.tpr

These files are:

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

Gromacs Quick Start