5.3 Creating input files of MD simulations using the Gromacs input files
Contents
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.
Preparations
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
6: AMBER99SB-ILDN
Solvation
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.
[INPUT] grotopfile = topol.top grocrdfile = proa_ionized.gro groreffile = proa_ionized.gro [OUTPUT] dcdfile = em.dcd rstfile = em.rst [MINIMIZE] nsteps = 5000 # number of steps eneout_period = 500 # energy output period crdout_period = 500 # coordinates output period rstout_period = 5000 # restart output period [BOUNDARY] 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] [SELECTION] group1 = rno:4-64 and heavy [RESTRAINTS] 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.
3. 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. Protein-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