5.3 Creating input files of MD simulations using the Gromacs input files

GENESIS can read .top and .gro files generated by Gromacs. Here we create input files including these two by using simulation setup programs of Gromacs.


This tutorial was made by using Gromacs 2021.4. Please refer to the Gromacs documentation web page (https://manual.gromacs.org/documentation/) for download and installation of the latest version.

1. Single-chain protein

We use PDB ID: 3MP9 in this example.Let’s make a pdb file containing only a protein chain. First, make proa.pdb according to Tutorial5.2.1 using VMD.

Now we create .gro, .top and .itp files using pdb2gmx of Gromacs.

$ gmx pdb2gmx -f proa.pdb -o proa.gro -water spce
# This will ask you to choose a force filed. Type 6 for this example.
> 6


Define a box with editconf to obtain proa_box.gro

$ gmx editconf -f proa.gro -o proa_box.gro -c -d 1.0 -bt cubic

This command places the protein at the center (-c) of a cubic box (-bt cubic), and at 1.0 nm from the box edge (-d 1.0). Let’s solvate the box with solvate.

$ gmx solvate -cp proa_box.gro -cs spc216.gro -o proa_solv.gro -p topol.top

This will generate proa_solv.gro and rewrite topol.top. -cs spc216.gro defines the configuration of the solvent as a generic equilibrated 3-point water models, such as SPC or TIP3P water.

Add ions

Create .tpr file using grompp, then add ions using genion

$ gmx grompp -f ions.mdp -c proa_solv.gro -p topol.top -o ions.tpr
$ gmx genion -s ions.tpr -o proa_ionized.gro -p topol.top -pname NA -nname CL –neutral

Input on GENESIS

Create a link to the force field parameter file in Gromacs. Make sure that the force field parameter file is located in the same directory as topol.top.

$ ln -s path_to_gromacs/share/top/amber99sb-ildn.ff ./

Here is an example of an input file for energy minimization of this system. Get the box size from the output of editconf (editconf.out) to specify them in the [BOUNDAY] section.

grotopfile      = topol.top
grocrdfile      = proa_ionized.gro
groreffile      = proa_ionized.gro

dcdfile         = em.dcd
rstfile         = em.rst

nsteps          = 5000     # number of steps
eneout_period   = 500      # energy output period
crdout_period   = 500      # coordinates output period
rstout_period   = 5000     # restart output period

type            = PBC      # [PBC,NOBC]
box_size_x      = 70.22     # box size (x) in [PBC]
box_size_y      = 70.22     # box size (y) in [PBC]
box_size_z      = 70.22     # box size (z) in [PBC]

group1          = rno:4-64 and heavy

nfunctions      = 1     # number of functions
function1       = POSI  # restraint function type
constant1       = 10.0  # force constant
select_index1   = 1     # restrained group

Note: GENESIS currently does not support CHARMM force fields with Gromacs input style. Please utilize the psf/top/par file formats when using CHARMM force fields (see Chapter 5.1 for more details).

2. Changing the protonation states of side-chains

pdb2gmx allows you to interactively choose side-chain protonation states of Lys, Arg, Asp, Glu, Gln, His using -lys/arg/asp/glu/gln/his flags. The protonation states of termini can be also interactively chosen with -ter flag.

Adding disulfide bonds

If the two cysteins are within a difined distance by pdb2gmx, you can interactively select disulfide bond formation using -ss flag for pdb2gmx. Alternatively, you can edit the contents of specbond.dat in gromacs and run pdb2gmx, followed by careful energy minimization of the system so that the disulfide bond length is relaxed.

4. P
rotein-DNA complex

In this section, we use PDB ID: 3LNQ. 3LNQ. Make a pdb file that contains only 1 protein and 2 dna chains (pro1_dnabc.pdb). Then run pdb2gmx as follows.

$ gmx pdb2gmx -f proa_dnabc.pdb -o proa_dnabc.gro -water tip3p

# Type 6 for selecting the force field
> 6
6: AMBER99SB-ILDN protein, nucleic AMBER94

Updated by Chigusa Kobayashi @RIKEN R-CCS

October 10, 2023

Written by Ai Niitsu@RIKEN Theoretical molecular science laboratory
April 1, 2022