5.2.2. gromacs.setup – Setting up a Gromacs MD run

Individual steps such as solvating a structure or energy minimization are set up in individual directories. For energy minimization one should supply appropriate mdp run input files; otherwise example templates are used.

Warning

You must check all simulation parameters for yourself. Do not rely on any defaults provided here. The scripts provided here are provided under the assumption that you know what you are doing and you just want to automate the boring parts of the process.

5.2.2.1. User functions

The individual steps of setting up a simple MD simulation are broken down in a sequence of functions that depend on the previous step(s):

topology()

generate initial topology file (limited functionality, might require manual setup)

solvate()

solvate globular protein and add ions to neutralize

energy_minimize()

set up energy minimization and run it (using mdrun_d)

em_schedule()

set up and run multiple energy minimizations one after another (as an alternative to the simple single energy minimization provided by energy_minimize())

MD_restrained()

set up restrained MD

MD()

set up equilibrium MD

Each function uses its own working directory (set with the dirname keyword argument, but it should be safe and convenient to use the defaults). Other arguments assume the default locations so typically not much should have to be set manually.

One can supply non-standard itp files in the topology directory. In some cases one does not use the topology() function at all but sets up the topology manually. In this case it is safest to call the topology directory top and make sure that it contains all relevant top, itp, and pdb files.

5.2.2.2. Example

Run a single protein in a dodecahedral box of SPC water molecules and use the GROMOS96 G43a1 force field. We start with the structure in protein.pdb:

from gromacs.setup import *
f1 = topology(protein='MyProtein', struct='protein.pdb', ff='G43a1', water='spc', force=True, ignh=True)

Each function returns “interesting” new files in a dictionary in such a away that it can often be used as input for the next function in the chain (although in most cases one can get away with the defaults of the keyword arguments):

f2 = solvate(**f1)
f3 = energy_minimize(**f2)

Now prepare input for a MD run with restraints on the protein:

MD_restrained(**f3)

Use the files in the directory to run the simulation locally or on a cluster. You can provide your own template for a queuing system submission script; see the source code for details.

Once the restraint run has completed, use the last frame as input for the equilibrium MD:

MD(struct='MD_POSRES/md.gro', runtime=1e5)

Run the resulting tpr file on a cluster.

5.2.2.3. User functions

The following functions are provided for the user:

gromacs.setup.topology(struct=None, protein='protein', top='system.top', dirname='top', posres='posres.itp', ff='oplsaa', water='tip4p', **pdb2gmx_args)

Build Gromacs topology files from pdb.

Keywords:

struct

input structure (required)

protein

name of the output files

top

name of the topology file

dirname

directory in which the new topology will be stored

ff

force field (string understood by pdb2gmx); default “oplsaa”

water

water model (string), default “tip4p”

pdb2gmxargs

other arguments for pdb2gmx

Note

At the moment this function simply runs pdb2gmx and uses the resulting topology file directly. If you want to create more complicated topologies and maybe also use additional itp files or make a protein itp file then you will have to do this manually.

gromacs.setup.solvate(struct='top/protein.pdb', top='top/system.top', distance=0.9, boxtype='dodecahedron', concentration=0, cation='NA', anion='CL', water='tip4p', solvent_name='SOL', with_membrane=False, ndx='main.ndx', mainselection='"Protein"', dirname='solvate', **kwargs)

Put protein into box, add water, add counter-ions.

Currently this really only supports solutes in water. If you need to embedd a protein in a membrane then you will require more sophisticated approaches.

However, you can supply a protein already inserted in a bilayer. In this case you will probably want to set distance = None and also enable with_membrane = True (using extra big vdw radii for typical lipids).

Note

The defaults are suitable for solvating a globular protein in a fairly tight (increase distance!) dodecahedral box.

Arguments:

structfilename

pdb or gro input structure

topfilename

Gromacs topology

distancefloat

When solvating with water, make the box big enough so that at least distance nm water are between the solute struct and the box boundary. Set boxtype to None in order to use a box size in the input file (gro or pdb).

boxtype or bt: string

Any of the box types supported by Editconf (triclinic, cubic, dodecahedron, octahedron). Set the box dimensions either with distance or the box and angle keywords.

If set to None it will ignore distance and use the box inside the struct file.

bt overrides the value of boxtype.

box

List of three box lengths [A,B,C] that are used by Editconf in combination with boxtype (bt in editconf) and angles. Setting box overrides distance.

angles

List of three angles (only necessary for triclinic boxes).

concentrationfloat

Concentration of the free ions in mol/l. Note that counter ions are added in excess of this concentration.

cation and anionstring

Molecule names of the ions. This depends on the chosen force field.

waterstring

Name of the water model; one of “spc”, “spce”, “tip3p”, “tip4p”. This should be appropriate for the chosen force field. If an alternative solvent is required, simply supply the path to a box with solvent molecules (used by genbox()’s cs argument) and also supply the molecule name via solvent_name.

solvent_name

Name of the molecules that make up the solvent (as set in the itp/top). Typically needs to be changed when using non-standard/non-water solvents. [“SOL”]

with_membranebool

True: use special vdwradii.dat with 0.1 nm-increased radii on lipids. Default is False.

ndxfilename

How to name the index file that is produced by this function.

mainselectionstring

A string that is fed to Make_ndx and which should select the solute.

dirnamedirectory name

Name of the directory in which all files for the solvation stage are stored.

includes

List of additional directories to add to the mdp include path

kwargs

Additional arguments are passed on to Editconf or are interpreted as parameters to be changed in the mdp file.

gromacs.setup.energy_minimize(dirname='em', mdp='/home/docs/checkouts/readthedocs.org/user_builds/gromacswrapper/checkouts/latest/gromacs/templates/em.mdp', struct='solvate/ionized.gro', top='top/system.top', output='em.pdb', deffnm='em', mdrunner=None, mdrun_args=None, **kwargs)

Energy minimize the system.

This sets up the system (creates run input files) and also runs mdrun_d. Thus it can take a while.

Additional itp files should be in the same directory as the top file.

Many of the keyword arguments below already have sensible values.

Keywords:

dirname

set up under directory dirname [em]

struct

input structure (gro, pdb, …) [solvate/ionized.gro]

output

output structure (will be put under dirname) [em.pdb]

deffnm

default name for mdrun-related files [em]

top

topology file [top/system.top]

mdp

mdp file (or use the template) [templates/em.mdp]

includes

additional directories to search for itp files

mdrunner

gromacs.run.MDrunner instance; by default we just try gromacs.mdrun_d() and gromacs.mdrun() but a MDrunner instance gives the user the ability to run mpi jobs etc. [None]

mdrun_args

arguments for mdrunner (as a dict), e.g. {'nt': 2}; empty by default

New in version 0.7.0.

kwargs

remaining key/value pairs that should be changed in the template mdp file, eg nstxtcout=250, nstfout=250.

Note

If mdrun_d() is not found, the function falls back to mdrun() instead.

gromacs.setup.em_schedule(**kwargs)

Run multiple energy minimizations one after each other.

Keywords:

integrators

list of integrators (from ‘l-bfgs’, ‘cg’, ‘steep’) [[‘bfgs’, ‘steep’]]

nsteps

list of maximum number of steps; one for each integrator in in the integrators list [[100,1000]]

kwargs

mostly passed to gromacs.setup.energy_minimize()

Returns:

dictionary with paths to final structure (‘struct’) and other files

Example:

Conduct three minimizations:

1. low memory Broyden-Goldfarb-Fletcher-Shannon (BFGS) for 30 steps

2. steepest descent for 200 steps

3. finish with BFGS for another 30 steps

We also do a multi-processor minimization when possible (i.e. for steep (and conjugate gradient) by using a gromacs.run.MDrunner class for a mdrun executable compiled for OpenMP in 64 bit (see gromacs.run for details):

import gromacs.run
gromacs.setup.em_schedule(struct='solvate/ionized.gro',
          mdrunner=gromacs.run.MDrunnerOpenMP64,
          integrators=['l-bfgs', 'steep', 'l-bfgs'],
          nsteps=[50,200, 50])

Note

You might have to prepare the mdp file carefully because at the moment one can only modify the nsteps parameter on a per-minimizer basis.

gromacs.setup.MD_restrained(dirname='MD_POSRES', **kwargs)

Set up MD with position restraints.

Additional itp files should be in the same directory as the top file.

Many of the keyword arguments below already have sensible values. Note that setting mainselection = None will disable many of the automated choices and is often recommended when using your own mdp file.

Keywords:

dirname

set up under directory dirname [MD_POSRES]

struct

input structure (gro, pdb, …) [em/em.pdb]

top

topology file [top/system.top]

mdp

mdp file (or use the template) [templates/md.mdp]

ndx

index file (supply when using a custom mdp)

includes

additional directories to search for itp files

mainselection

make_ndx selection to select main group [“Protein”] (If None then no canonical index file is generated and it is the user’s responsibility to set tc_grps, tau_t, and ref_t as keyword arguments, or provide the mdp template with all parameter pre-set in mdp and probably also your own ndx index file.)

deffnm

default filename for Gromacs run [md]

runtime

total length of the simulation in ps [1000]

dt

integration time step in ps [0.002]

qscript

script to submit to the queuing system; by default uses the template gromacs.config.qscript_template, which can be manually set to another template from gromacs.config.templates; can also be a list of template names.

qname

name to be used for the job in the queuing system [PR_GMX]

mdrun_opts

option flags for the mdrun command in the queuing system scripts such as “-stepout 100”. [“”]

kwargs

remaining key/value pairs that should be changed in the template mdp file, eg nstxtcout=250, nstfout=250 or command line options for grompp` such as ``maxwarn=1.

In particular one can also set define and activate whichever position restraints have been coded into the itp and top file. For instance one could have

define = “-DPOSRES_MainChain -DPOSRES_LIGAND”

if these preprocessor constructs exist. Note that there must not be any space between “-D” and the value.

By default define is set to “-DPOSRES”.

Returns:

a dict that can be fed into gromacs.setup.MD() (but check, just in case, especially if you want to change the define parameter in the mdp file)

Note

The output frequency is drastically reduced for position restraint runs by default. Set the corresponding nst* variables if you require more output. The pressure coupling option refcoord_scaling is set to “com” by default (but can be changed via kwargs) and the pressure coupling algorithm itself is set to Pcoupl = “Berendsen” to run a stable simulation.

gromacs.setup.MD(dirname='MD', **kwargs)

Set up equilibrium MD.

Additional itp files should be in the same directory as the top file.

Many of the keyword arguments below already have sensible values. Note that setting mainselection = None will disable many of the automated choices and is often recommended when using your own mdp file.

Keywords:

dirname

set up under directory dirname [MD]

struct

input structure (gro, pdb, …) [MD_POSRES/md_posres.pdb]

top

topology file [top/system.top]

mdp

mdp file (or use the template) [templates/md.mdp]

ndx

index file (supply when using a custom mdp)

includes

additional directories to search for itp files

mainselection

make_ndx selection to select main group [“Protein”] (If None then no canonical index file is generated and it is the user’s responsibility to set tc_grps, tau_t, and ref_t as keyword arguments, or provide the mdp template with all parameter pre-set in mdp and probably also your own ndx index file.)

deffnm

default filename for Gromacs run [md]

runtime

total length of the simulation in ps [1000]

dt

integration time step in ps [0.002]

qscript

script to submit to the queuing system; by default uses the template gromacs.config.qscript_template, which can be manually set to another template from gromacs.config.templates; can also be a list of template names.

qname

name to be used for the job in the queuing system [MD_GMX]

mdrun_opts

option flags for the mdrun command in the queuing system scripts such as “-stepout 100 -dgdl”. [“”]

kwargs

remaining key/value pairs that should be changed in the template mdp file, e.g. nstxtcout=250, nstfout=250 or command line options for :program`grompp` such as maxwarn=1.

Returns:

a dict that can be fed into gromacs.setup.MD() (but check, just in case, especially if you want to change the define parameter in the mdp file)

5.2.2.4. Helper functions

The following functions are used under the hood and are mainly useful when writing extensions to the module.

gromacs.setup.make_main_index(struct, selection='"Protein"', ndx='main.ndx', oldndx=None)

Make index file with the special groups.

This routine adds the group __main__ and the group __environment__ to the end of the index file. __main__ contains what the user defines as the central and most important parts of the system. __environment__ is everything else.

The template mdp file, for instance, uses these two groups for T-coupling.

These groups are mainly useful if the default groups “Protein” and “Non-Protein” are not appropriate. By using symbolic names such as __main__ one can keep scripts more general.

Returns:

groups is a list of dictionaries that describe the index groups. See gromacs.cbook.parse_ndxlist() for details.

Arguments:

structfilename

structure (tpr, pdb, gro)

selectionstring

is a make_ndx command such as "Protein" or r DRG which determines what is considered the main group for centering etc. It is passed directly to make_ndx.

ndxstring

name of the final index file

oldndxstring

name of index file that should be used as a basis; if None then the make_ndx default groups are used.

This routine is very dumb at the moment; maybe some heuristics will be added later as could be other symbolic groups such as __membrane__.

gromacs.setup.check_mdpargs(d)

Check if any arguments remain in dict d.

gromacs.setup.get_lipid_vdwradii(outdir='.', libdir=None)

Find vdwradii.dat and add special entries for lipids.

See gromacs.setup.vdw_lipid_resnames for lipid resnames. Add more if necessary.

gromacs.setup._setup_MD(dirname, deffnm='md', mdp='/home/docs/checkouts/readthedocs.org/user_builds/gromacswrapper/checkouts/latest/gromacs/templates/md_OPLSAA.mdp', struct=None, top='top/system.top', ndx=None, mainselection='"Protein"', qscript='/home/docs/checkouts/readthedocs.org/user_builds/gromacswrapper/checkouts/latest/gromacs/templates/local.sh', qname=None, startdir=None, mdrun_opts='', budget=None, walltime=0.3333333333333333, dt=0.002, runtime=1000.0, **mdp_kwargs)

Generic function to set up a mdrun MD simulation.

See the user functions for usage.

Defined constants:

gromacs.setup.CONC_WATER = 55.345

Concentration of water at standard conditions in mol/L. Density at 25 degrees C and 1 atmosphere pressure: rho = 997.0480 g/L. Molecular weight: M = 18.015 g/mol. c = n/V = m/(V*M) = rho/M = 55.345 mol/L.

gromacs.setup.vdw_lipid_resnames = ['POPC', 'POPE', 'POPG', 'DOPC', 'DPPC', 'DLPC', 'DMPC', 'DPPG']

Hard-coded lipid residue names for a vdwradii.dat file. Use together with vdw_lipid_atom_radii in get_lipid_vdwradii().

gromacs.setup.vdw_lipid_atom_radii = {'C': 0.25, 'H': 0.09, 'N': 0.16, 'O': 0.155}

Increased atom radii for lipid atoms; these are simply the standard values from GMXLIB/vdwradii.dat increased by 0.1 nm (C) or 0.05 nm (N, O, H).