OpenMM Functionality¶
Since ParmEd is a library for reading, converting, and modifying full molecular mechanical descriptions of chemical systems in a wide variety of different families of force fields, supporting molecular simulation directly using the fantastic OpenMM Python API was a natural extension.
This page is not meant as an exhaustive description of the OpenMM library and its usage. Instead, you should visit the OpenMM website for that. However, this page will provide a brief description of OpenMM, what I consider to be its strengths that can make it an invaluable tool in the field of molecular mechanics. It will also present an introduction to using OpenMM with the tools provided by ParmEd through a handful of examples.
What is OpenMM¶
OpenMM is not a program in the traditional sense. There is no OpenMM program that you can run from the terminal like you can with AMBER, CHARMM, or GROMACS, for example. Instead, it is a library of routines written in C++ that allow you to program your own molecular models in C, C++, Fortran, or, as we will demonstrate here, Python. Its basic features include:
Basic valence-term forces, like bonds, angles, torsions, and coupled-torsion correction maps.
Nonbonded potential terms, like the typical electrostatic and Lennard-Jones potentials, Generalized Born implicit solvent models, and more, with accurate long-range potentials.
Various integrators to carry out molecular dynamics, like the Verlet integrator for traditional MD, or stochastic integrators like those for Langevin and Brownian dynamics.
Thermostats and barostats for sampling from various statistical ensembles.
What makes OpenMM so awesome?¶
There are two key features that not only make OpenMM awesome, but make it unlike any other molecular dynamics software in existence.
Stellar performance on GPUs¶
Among the primarily advertised selling-points of OpenMM is the computational performance that is possible by running OpenMM on NVidia and AMD/ATI graphics processing units (GPUs). OpenMM utilizes both CUDA and OpenCL to program molecular models for GPUs, and runs the entire calculation on the accelerator device to eliminate the GPU-to-CPU communications that can significantly limit performance.
In general, the CUDA platform is faster than the OpenCL platform, but OpenCL works on a much wider array of hardware. I will not expound on the details of the GPU performance, as you are directed to look at the various GPU-related publications for such comparisons. Suffice it to say that, currently (February, 2015), GPU performance is ca. one, possibly up to two orders of magnitude faster than running on every core of a standard multicore server processor (8 to 16 cores) with a multithreaded MD application.
Custom Forces¶
In my opinion, what makes OpenMM truly awesome is its custom force capabilities. OpenMM allows you to implement a new potential for bonds, angles, torsions, and even nonbonded pairwise and multi-particle interactions simply by providing an equation for the potential energy.
OpenMM will then analytically differentiate this analytical function to get an analytical, closed-form expression for the gradients (forces), and build a modestly optimized GPU implementation of that term!
For example, suppose you want to model a few bonds in your biomolecule using a typical Morse potential, the OpenMM Python code that will do this is shown below:
force = CustomBondForce("D*(1-exp(-alpha*(r-req)))^2")
And that’s it! Once you add the bonds to that force, with their parameters, OpenMM will build an optimized kernel for the calculation. Not only does this eliminate the development cost of building and optimizing a GPU implementation of a new potential, it also provides a pre-tested, and pre-debugged version of that kernel.
How does ParmEd enhance OpenMM?¶
|
This class acts as a state data reporter for OpenMM simulations, but it is a little more generalized. |
|
NetCDFReporter prints a trajectory in NetCDF format |
|
MdcrdReporter prints a trajectory in ASCII Amber format. |
|
Use a reporter to handle writing restarts at specified intervals. |
|
A class that prints out a progress report of how much MD (or minimization) has been done, how fast the simulation is running, and how much time is left (similar to the mdinfo file in Amber) |
|
Wrapper for parsing OpenMM-serialized objects. |
|
Creates a |
|
This computes the energy of every force group in the given structure and computes the energy for each force group for the given Context. |
ParmEd provides a common framework for constructing and representing fully parametrized force field models for various systems instantiated from a wide range of file formats. In particular, it supports the following file formats:
Standard Amber topology file
Amber Chamber-style topology file*
Amber AMOEBA-style topology file*
CHARMM PSF file
Gromacs topology file
The -marked options are *not available in the OpenMM Python application layer. While the OpenMM Python application layer does support standard Amber topology files, it does not support either old-style topology files or those that define a 10-12 nonbonded potentials between certain pairs of atoms. ParmEd supports both of these. With Gromacs topology files, the OpenMM Python application layer does not correctly handle virtual sites, while the Gromacs topology file parser in ParmEd does.
In addition to these file formats, ParmEd also supports several new reporter classes in addition to the small number provided by ParmEd:
StateDataReporter– This takes extra arguments specifying the units of each of the types of data (like the energy, volume, and time units). The defaults correspond to the AKMA unit system, which is more familiar to Amber and CHARMM users.
NetCDFReporter– This allows you to write a trajectory file in the Amber NetCDF format.
MdcrdReporter– This allows you to write a trajectory file in the Amber ASCII mdcrd format.
RestartReporter– This allows you to periodically write NetCDF or ASCII restart files during the course of the calculation
ProgressReporter– This prints a file during the course of the simulation tracking the runtime speed of the calculation and predicting the amount of time remaining.
The energy_decomposition() function takes as input a Structure instance, OpenMM Context, and an optional
energy unit (nrg) and returns a dictionary of all energy components for the
different force groups. This permits an form of energy decomposition that allows
energy components to be compared between programs more effectively. For
example:
>>> import parmed as pmd
>>> from simtk.openmm import app
>>> from simtk import openmm as mm
>>> # Instantiate the parm and create the system
... parm = pmd.load_file('tip4p.parm7', 'tip4p.rst7')
>>> system = parm.createSystem(nonbondedMethod=app.PME,
... nonbondedCutoff=8*pmd.unit.angstrom)
>>> # Make the context and set the positions
... context = mm.Context(system, mm.VerletIntegrator(0.001))
>>> context.setPositions(parm.positions)
>>> # Find the energy decomposition
... pmd.openmm.energy_decomposition(parm, context)
{'total': -2133.295388974015, 'nonbonded': -2133.2953890231834, 'bond': 1.0126508125518531e-07}
Loading OpenMM Objects¶
The OpenMM Topology object and System object contain the information
stored in Structure. You can use
load_topology() to load an OpenMM Topology and create a
Structure instance from it. If you provide
either a file containing a serialized System in XML format or a System
object directly, parameters will be extracted from the various forces and added
to the generated Structure. You can also
pass coordinates (or any coordinate file, including an OpenMM XML-serialized
State file) with the xyz argument and unit cell dimensions with the box
argument (which will override any unit cell information contained in the input
Topology, System or coordinate file if applicable).
The XmlFile class can parse and return a deserialized object from an
OpenMM-generated XML file. Supported XML files are:
Serialized
System(returns an OpenMMSysteminstance)Serialized
State(returns a container object with attributescoordinates,velocities,forces,energy, andtime)Serialized
Integrator(returns an OpenMMIntegratorsubclass)XML
ForceFieldfile (returns an OpenMMForceFieldinstance)
Examples¶
Right now, the two main pathways to run an OpenMM simulation starting from a fully parametrized molecular mechanical model (click on either option to view the annotated and explained example):
Taking OpenMM Topology and System to a ParmEd Structure¶
While the above sections described how you would generate an OpenMM System
and Topology instance from any of a number of file formats (e.g., Amber
topology file, Gromacs topology file, or CHARMM PSF file), it is also possible
to go in the reverse direction. That is, to start with Topology and
System instances and convert those to a Structure instance.
The function for this is load_openmm, and takes a mandatory Topology
instance, along with either an optional System instance or name of a
serialized System XML file defining an OpenMM System, and returns a
populated Structure from it. This is
particularly useful when you parametrize a system using the OpenMM modelling
capabilities, but want to use that parametrized system in another program, like
Amber or Gromacs.
An example is shown below, using the OpenMM functionality to parametrize a PDB file with the ff99SB force field:
import parmed as chem
import parmed.unit as u
from simtk.openmm import app
from simtk import openmm as mm
pdb = app.PDBFile('input.pdb')
forcefield = app.ForceField('amber99sb.xml', 'tip3p.xml')
system = forcefield.createSystem(pdb.topology, nonbondedMethod=app.PME,
nonbondedCutoff=1*u.nanometer)
struct = chem.openmm.load_topology(pdb.topology, system)
There are some limitations to this functionality, itemized below:
ParmEd does not recognize all of the different OpenMM Systems it can generate, such as any of those features implemented using a
CustomNonbondedForce(e.g., NBFIX, 10-12 potential, 12-6-4 potential, etc.).You should make sure to create the
Systemwith no constraints, since OpenMM may be missing bond or angle terms associated with the constrained degrees of freedom.With the exception of certain
CustomTorsionForcewhich are recognized as quadratic improper torsions, the presence of anyCustomForceinstances prevents ParmEd from recognizing the potential.