GromacsWrapper
GromacsWrapper is a python package that wraps system calls to Gromacs tools into thin classes. This allows for fairly seamless integration of these tools into Python scripts. This is generally superior to shell scripts because of Python’s better error handling and superior data structures. It also allows for modularization and code re-use. In addition, commands, warnings and errors are logged to a file so that there exists a complete history of what has been done.
Availability
Source code is available from the GromacsWrapper github repository under the GNU Public Licence, together with the online documentation.
Example
As an example we show some high-level functions in the gromacs.setup module, which make use of the GromacsWrapper library.
Setting up and running a simulation
We want to 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 topology, solvate, energy_minimize, MD_restrained, MD
filenames = topology(protein='MyProtein', struct='protein.pdb', ff='G43a1', water='spc', force=True, ignh=True)
filenames = solvate(**filenames)
filenames = energy_minimize(**filenames)
MD_restrained(runtime=1000, qscript=["local.sh", "cluster.pbs"], **filenames)
# run position restraint simulation: cd MD_POSRES && qsub cluster.pbs
MD(struct='MD_POSRES/md.gro', runtime=1e5, qscript=["local.sh", "cluster.pbs"])
# run equilibrium MD simulation: cd MD && qsub cluster.pbs
The next section explains all steps in a bit more detail. See also gromacs.setup: Example.
Step-by-step
First import all functions from the module
from gromacs.setup import *
(This saves some typing but of course one can simply say import gromacs.setup and call each function as gromacs.setup.FUNCTION().)
We want to 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 and generate topology with the topology() function:
filenames = topology(protein='MyProtein', struct='protein.pdb', ff='G43a1', water='spc', force=True, ignh=True)
Each function returns “interesting” new files in a dictionary filenames in such a way that it can 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).
Add the solvent to the protein in the simulation box with the solvate() function:
filenames = solvate(**filenames)
Now we run a simple steepest descent energy minimization with energy_minimze()
filenames = energy_minimize(**filenames)
(It is also possible to run different kinds of energy minimizers in series with gromacs.setup.em_schedule().)
Once properly minimized, we prepare a 1000-ps simulation with position restraints on the protein with the help of MD_restrained()
MD_restrained(runtime=1000, qscript=["local.sh", "cluster.pbs"], **filenames)
Because this simulation can take a while you would typically run it on a cluster. If you provide a template for a queuing system script then the function will create an appropriate script for you so that you can directly submit the job:
cd MD_POSRES && qsub cluster.pbs
(Alternatively, run the simulation locally with a command such as mdrun -v -deffnm md -cpi.)
From the final frame of the position restraint run we set up the unrestrained equilibrium MD simulation (using the gromacs.setup.MD() function) – in this example for 100 ns (100,000 ps):
MD(struct='MD_POSRES/md.gro', runtime=1e5, qscript=["local.sh", "cluster.pbs"])
Submit the simulation with
cd MD && qsub cluster.pbs
The gromacs.analysis module contains plugins for a range of analysis tasks.
For further details see the GromacsWrapper documentation
Discuss: “GromacsWrapper”