|Home | About | Journals | Submit | Contact Us | Français|
OpenMM is a software toolkit for performing molecular simulations on a range of high performance computing architectures. It is based on a layered architecture: the lower layers function as a reusable library that can be invoked by any application, while the upper layers form a complete environment for running molecular simulations. The library API hides all hardware-specific dependencies and optimizations from the users and developers of simulation programs: they can be run without modification on any hardware on which the API has been implemented. The current implementations of OpenMM include support for graphics processing units using the OpenCL and CUDA frameworks. In addition, OpenMM was designed to be extensible, so new hardware architectures can be accommodated and new functionality (e.g., energy terms and integrators) can be easily added.
Molecular mechanics simulation is a powerful and widely used technique for studying the behavior of biological macromolecules. Many software applications exist for performing these simulations. Although these applications are powerful and feature rich, nearly all of them share some important limitations:
These limitations can create bottlenecks for the molecular simulation community. For a researcher developing new algorithms, prototyping their ideas inside an existing simulation code is typically unnecessarily difficult and time consuming. Once an algorithm is shown to work, it may take years to become available in all the major applications. Moreover, it is frequently impractical (and in many cases, not permitted by the license) for anyone but the core developers of an application to add the new feature and make it available to the community. Even after one application adds the feature, that does not bring any other application closer to having it; their architectures are different enough that most features must be implemented independently for every application.
In this paper, we describe OpenMM, a molecular simulation code designed to address these challenges. It is architected from the ground up to be extensible and reusable, while still permitting very high performance. At its core, it is not just an application; it is a library that can be easily used by any application. It is based around a modular plugin architecture allowing new features to be added without modifying the OpenMM library itself. Finally, it is available under a permissive open source license allowing it to be used for any purpose, including commercial applications.
OpenMM was designed with the following goals in mind:
OpenMM aims to make it very easy for programmers to add molecular simulation features to their applications. This means having a simple application programming interface (API) that is easy to learn and understand, and that requires minimal code to use. It also means that OpenMM should be easy to use correctly. The API should minimize the opportunities for programmers to make mistakes and guide them toward making good choices.
Molecular simulations are commonly run on many types of hardware, including conventional processors, supercomputing clusters, graphics processing units (GPUs)5–8, and dedicated special purpose processors.9 Technology continues to change quickly, and it is difficult to predict what hardware platforms will be important to support as little as five years in the future.
OpenMM acts as a hardware abstraction layer, similar to LAPACK for linear algebra or OpenGL for graphics. It provides a hardware independent API for users to call from their programs. Implementations of that API can be written for many different types of hardware, and any program that uses OpenMM can run efficiently and without modification on any of those hardware platforms—even ones that did not exist yet when the program was written.
OpenMM provides a large number of simulation features (molecular force fields, integration methods, etc.), but it cannot hope to include every possible feature any user might want. Instead, it is designed to be extensible so users can add new features without modifying the OpenMM library itself. This allows users to implement the features they need and, if they wish, distribute them to the molecular simulation community. Any program that uses OpenMM can easily add support for those features without needing to reimplement them.
Whatever other benefits a simulation code offers, it will not become widely used if its performance is not competitive with other applications. OpenMM is designed to permit efficient implementation on a wide variety of hardware platforms. This requires careful thought in the design of the API to avoid making assumptions about data locality, the latency of operations, or other features that vary between platforms. Otherwise, one could easily create an API that was impossible to implement efficiently on certain types of hardware. OpenMM strives to avoid these pitfalls. It includes very fast GPU implementations of all its standard features, giving users instant access to high performance simulation code.
OpenMM is based on a modular architecture. It consists of the following major components:
Each of these components is described in more detail below.
The most recent version of OpenMM at the time of this writing is 4.1.1, and the descriptions contained below correspond to this version. Many of these features were present in earlier versions, but some important features were added only recently. Most importantly, the application layer was first introduced in version 4.0, and custom integrators were added in 4.1. The AMOEBA force field was first added in version 3.0. Custom forces have existed since the very first stable release, but the set of custom force types has grown steadily over time.
The OpenMM core provides the standard covalent and noncovalent interactions used in most biomolecular force fields. Two implicit solvent models are also included in the OpenMM core. The covalent energy terms implemented in OpenMM include the standard harmonic bond and angle energies and the Ryckaert-Bellemans and periodic torsion energies. The CHARMM cross-term energy correction map, CMAP torsion, is also available.10 The CMAP torsion term couples two dihederal angles using a tabulated potential energy surface.
The nonbonded interactions implemented in the OpenMM core include the Lennard-Jones and Coulomb energies. These interactions may be calculated with or without a distance cutoff; if a cutoff is applied, it is the same for both sets of interactions.
The Lorentz-Bertelot combining rule is used for Lennard-Jones interactions. If a distance cutoff and periodic boundary conditions are applied, an analytical, isotropic long-range dispersion correction may be included that approximately represents the Lennard-Jones energy contribution from all interactions beyond the cutoff distance.11
The Coulomb interaction can be calculated with one of five options:
If no distance cutoff is applied, the standard Coulomb energy is included for each particle pair. If a cutoff is applied, then the Coulomb energy is modeled using the reaction-field approximation.14 This approximation is derived by assuming molecules beyond the cutoff distance are solvent molecules and the dielectric constant is uniform.
For the Ewald summation and particle-mesh Ewald (PME) options, the user specifies the direct space cutoff and an error tolerance. The program then determines the values of the internal parameters, including the reciprocal space cutoff for Ewald and the mesh size for PME, required to calculate the forces to the specified tolerance.
Two implicit solvent forces are available in the OpenMM core: the OBC Generalized Born model15, and the GB/VI model developed by Labute.16 The OBC model is comprised of a polarization term to represent the electrostatic interaction between the solute and solvent, and a surface area term to represent the free energy of solvating a neutral molecule.
The GB/VI model includes a polarization term similar to the one in the OBC model, but with the replacement of the exponent in the solute volume integral with the value −6 instead of −4. In addition, the GB/VI model employs a volumetric dispersion term to account for the nonpolar contributions to the solvation free energy.
Five integrators are available in the OpenMM core: three fixed time step integrators (Verlet, Langevin, and Brownian), and two variable time step integrators (Verlet and Langevin). The Verlet and Langevin integrators implement the leap-frog versions of their respective algorithms.
The variable time step integrators are variants of the respective fixed time step integrators that continuously adjust the step size to keep the integration error below a user-specified tolerance.17 This is accomplished by comparing the positions generated by Verlet integration with those that would be generated by an explicit Euler integrator, and taking the difference between them as a conservative estimate of the integration error. The step size is then adjusted to make the estimated error equal the specified error tolerance.
An Andersen thermostat and Monte Carlo barostat are also provided in the OpenMM core. The thermostat couples the system to a heat bath by randomly selecting a subset of particles at the start of each time step, then setting their velocities to new values chosen from a Maxwell-Boltzmann distribution. This represents the effect of random collisions between particles in the system and particles in the heat bath.18 This is only needed when using a Verlet integrator, since the Langevin and Brownian integrators provide their own temperature coupling.
The Monte Carlo barostat models the effect of constant pressure by allowing the size of the periodic box to vary during the simulation.19–20 It does this by periodically rescaling the box, computing the change in energy, then accepting or rejecting the change based on a Metropolis criterion.
In OpenMM, a “platform” is an implementation of all the computations defined by the API. The architecture is modular so that new platforms may be added and the client code can select at run time which platform to use for a simulation. Three platforms are included with OpenMM:
This platform is written in C++, and serves as a reference when writing other platforms. The emphasis is on clarity and correctness, not performance; it intentionally avoids any optimization that would reduce the clarity of the code. This platform is generally too slow to be useful for production simulations, but it serves two important functions. First, it can be used as a starting point for writing other platforms, since it contains complete implementations of all required algorithms. Second, it serves as a standard for “correct” behavior. Results produced by other platforms can be compared to the reference platform to verify their correctness.
This platform uses the OpenCL framework21 for its computational kernels. This allows it to run efficiently on Nvidia and AMD GPUs, as well as on multicore CPUs. It also can parallelize computations across multiple GPUs in a single computer. In most cases, this is the fastest platform, and it is the recommended platform for most production simulations.
This platform uses Nvidia’s CUDA framework22 for running on GPUs. Unlike the OpenCL platform, it only supports Nvidia GPUs. It also is more limited in other respects: it cannot parallelize computations across multiple GPUs, and it does not support as many types of custom forces (described below) as the OpenCL platform.
The CUDA platform support in OpenMM is derived from an older code base which makes it very difficult to support many of the newer features of OpenMM. We are currently in the process of writing a completely new CUDA-based platform modeled after the OpenCL platform. Once it is released in a future version of OpenMM, CUDA will support all the same features as OpenCL, as well as having significantly better performance on recent Nvidia GPUs.
The OpenMM core library is focused on the essential computational tasks of running a molecular simulation: computing forces and energies, and integrating the equations of motion. To actually run a simulation, however, other features are also required: loading and saving files in various standard formats, building models based on force field descriptions, etc. That is what the application layer does. It consists of a set of Python libraries that perform these functions, allowing it to be used as a complete simulation application.
Strictly speaking, of course, a set of libraries does not constitute an application, and the user is responsible for writing a true application that calls them. In practice, they are complete and simple enough that they may collectively be thought of as an “application”, in which the user provides a control script. This is illustrated by Listing 1. It is a complete Python script that loads a PDB file, builds a mathematical model from it using the AMBER99SB force field and TIP3P water model, does a local energy minimization, and simulates it for 1 million time steps, saving a frame in PDB format every 1000 steps. All of this is done in only a dozen lines of code. The script is no more complicated than the control scripts used by many other MD applications, and its logic can be easily understood even by users with no programming experience.
from simtk.pyopenmm.app import * from simtk.openmm import * from simtk.unit import * pdb = PDBFile(‘5dfr_solv-cube_equil.pdb’) ff = ForceField(‘amber99sb.xml’, ‘tip3p.xml’) sys = ff.createSystem(pdb.topology, nonbondedMethod=PME, nonbondedCutoff=1*nanometer, constraints=HBonds) integ = LangevinIntegrator (300*kelvin, 1/picosecond, 0.002*picoseconds) sim = Simulation (pdb.topology, sys, integ) sim.context.setPositions (pdb.positions) sim.minimizeEnergy ( ) sim.reporters.append (PDBReporter (‘output.pdb’, 1000)) sim.step (1000000)
The advantage of this approach is its extreme flexibility. A user who simply wants to run simulations can take an example script such as this and modify it in clearly documented ways to suit their needs. More advanced users have direct access to the entire OpenMM API, as well as the full power of the Python language and libraries. For example, Listing 2 shows a complete simulated annealing algorithm implemented in only three lines:
for i in range(100): integ.setTemperature(3*(100-i)*kelvin) sim.step(1000)
Extensibility is a fundamental feature of OpenMM. We cannot hope to provide every possible feature that any user might want, so it is important that other people can add those features themselves. OpenMM is also designed to be useful for researchers developing new simulation methods. By allowing them to easily prototype and test new ideas, it serves as a powerful tool for algorithm development.
Extensibility is provided through two mechanisms: plugins and custom forces/integrators.
A plugin is a piece of code that adds new features to OpenMM. Typically, each plugin is packaged as a dynamically linked library. At run time, OpenMM can be told to load all the plugins found in a directory and make them available to the running program.
Plugins generally fall into two categories. First, there are plugins that add new features to OpenMM. This is illustrated by the AMOEBA plugin described below. It defines a set of new forces that users can include in their simulations. Plugins can also define new thermostats, barostats, and integrators. Programs can use features defined by plugins just as easily as those defined by the core library. In fact, the lower levels of the architecture do not distinguish at all between core features and those provided by plugins. All features are handled identically, wherever they are defined.
Another type of plugin is ones that implement new platforms. In fact, the OpenCL and CUDA platforms are actually implemented in plugins, not in the core library. Because plugins can add new platforms, and because the choice of what platform to use is made at run time, any program that uses OpenMM has a high level of hardware independence. If a new type of hardware becomes available, a plugin can be written to support it, and all existing programs can take advantage of it with no changes required.
Three plugins that provide additional features beyond those in the core library are bundled with OpenMM: the AMOEBA force field, the free energy plugin, and the Ring Polymer Molecular Dynamics (RPMD) plugin. An outline of each plugin is given below.
AMOEBA (Atomic Multipole Optimized Energetic for Biomolecular Applications) is a polarizable force field23–24 developed in Jay Ponder’s lab. The AMOEBA force field was designed with the goal of obtaining chemical accuracy on the order of 2 kJ/mol for small molecule and protein-ligand interactions. To realize this level of accuracy, both the covalent and noncovalent interactions employed in AMOEBA differ from those employed in more common force fields.
The bond-stretch, angle-bending and bond-angle interactions are based on the corresponding terms in the MM3 force field25; anharmonicity is included in the bond-stretching term up to fourth order and in the angle-bending term up to sixth order. A Wilson-Decius-Cross interaction26 is used at sp2-hybridized trigonal centers to limit out-of-plane bending. A periodic torsion term is included and is based on the typical Fourier expansion. In addition a Bell torsion term27 is applied to dihederal angles involving two trigonal centers. Finally a CMAP torsion interaction coupling two dihedral angles is also included.
The van der Waals interactions are modeled using a buffered 14-7 functional form.28 One significant difference between the AMOEBA implementation and others is that the van der Waals interaction site for hydrogen atoms is not at the hydrogen nucleus, but at a position on the H-X bond, where X is the hydrogen’s covalent partner. The interaction site is typically about 90% of the distance from the X nucleus along the H-X bond. This modification has been shown to improve the fit to water dimer structures based on quantum mechanical calculations.
The major difference between AMOEBA and traditional force fields is in AMOEBA’s treatment of the electrostatic interaction. There are two major differences associated with the electrostatics: for each atom, AMOEBA includes a permanent charge, dipole and quadrupole, in contrast to traditional force fields which typically only include a permanent charge. The permanent electrostatic components are based on ab initio calculations performed on small molecules. The second major difference between conventional force fields and AMOEBA in the treatment of the electrostatic interactions is the inclusion of electronic polarization in AMOEBA. In traditional force fields polarization is only treated in an average way. In AMOEBA polarization is modeled as induced dipoles at each atom. The induced dipoles are determined by the atomic polarizability and the electric field at the site. The atomic polarizability in turn is based on Thole’s damped interaction model.29 Because the electric field at each site includes contributions from the induced dipoles of neighboring sites, the calculation of the induced dipoles must be iterated until the induced dipoles are self-consistent. The iterative process must be carried out at each time step and depending on the accuracy required can be the rate-limiting step in the calculation of the forces.
The AMOEBA plugin uses particle-mesh Ewald summation to calculate the electrostatic nonbonded interaction for simulations with explicit solvent. For implicit solvent simulations, the generalized Kirkwood algorithm is used.30
The free energy plugin allows “alchemical” free energy calculations to be performed. In these types of calculations, subsets of particles are inserted, removed, or transformed into other particle types in stages. The staging is performed using a parameter λ that modulates the interactions of the particles being transformed with the rest of the system: λ=0 results in no interactions with the rest of the system, while λ=1 recovers the standard interaction terms. At intermediate values of λ during particle insertion or deletion, a “soft-core” form of the interactions is employed that allows particles to overlap without generating numerical instabilities.31 The soft-core form of the Lennard-Jones potential is given by32:
However, other choices for free energy pathways33 can be easily implemented using custom force as described below. The free energy plugin also includes modifications to the OBC and GB/VI implicit solvent models to allow alchemical free energy calculations. For the OBC model, the contributions of transformed particles to the cavity energy term are scaled by λ. In addition, for both implicit solvent models, the contribution of transformed particles to the calculation of the Born radii is scaled by λ. Free energies can be computed by evaluating the system energy at additional values of λ and applying the Bennett acceptance ratio (BAR)34 or the multistate Bennett acceptance ratio (MBAR)35. Derivatives of the potential are not possible in this framework, so thermodynamic integration (TI) cannot be used. This is not a major obstacle for most applications, as BAR and MBAR are as efficient as or more efficient than TI for standard free energy calculations36–37.
This plugin supports running simulations with the ring-polymer molecular dynamics (RPMD) method.38 This is a method for incorporating limited quantum mechanical effects, such as tunneling and zero-point energy, into a classical molecular simulation. It does this by simulating many copies of a system at once, which are connected by harmonic springs to form a “ring polymer”. The plugin implements the path integral Langevin equation (PILE) thermostat algorithm to efficiently sample the static and dynamic properties of quantum systems.39
Custom forces are a unique feature of OpenMM that allow a wide range of novel interactions to be implemented with very little effort, yet with high speed of execution. The user provides a mathematical description of the interaction in the form of a simple text string. OpenMM parses the expression, figures out how to compute it efficiently, and then uses it exactly as if it were a standard force. All of this is done at run time, so no compilation step is needed.
Listing 3 gives an example. It uses the CustomBondForce class to implement a Morse potential. This class implements pairwise bonded interactions whose energy is an arbitrary function of the inter-atomic distance r. This example specifies that the energy of each bond is given by
It also specifies that D, a, and r0 are per-bond parameters, so a value of each one will be specified for each bond.
bonds = CustomBondForce(“D*(1-exp(a*(r0-r)))^2”) bonds.addPerBondParameter(“D”) bonds.addPerBondParameter(“a”) bonds.addPerBondParameter(“r0”)
This code is very simple, but it hides a much more complicated computation. OpenMM parses the mathematical expression, analytically differentiates it to get an expression for the force, applies a series of optimizations to both expressions, generates OpenCL code for evaluating both force and energy, synthesizes a complete kernel to loop over bonds and evaluate each one, and sends this kernel to the OpenCL driver for compilation to the GPU’s native instruction set. All of this is done at run time, and takes only a fraction of a second. The custom force can then be used exactly like any other force, and runs at the full speed of the GPU hardware. The details of this process are described elsewhere.40
OpenMM includes several custom force classes. In addition to the pairwise bonded force shown above, there are custom forces that depend on the angle between three atoms, the dihedral angle formed by four atoms, and the pairwise distance between nonbonded atoms. There also are more complicated custom forces. These include CustomHbondForce, which can implement a wide range of hydrogen bonding potentials, and CustomGBForce, which provides a flexible framework for implementing Generalized Born implicit solvent models. Listing 4 shows the essential code required to implement the OBC implicit solvent model using the CustomGBForce. This code may look complicated, but it would take roughly 1000 lines of code to implement the same force by more conventional means. Moreover, since it expresses the nature of the calculation directly in terms of equations (rather than OpenCL or CUDA code), it is generally much more transparent and easier to modify than traditional approaches. We have used CustomGBForce to implement several other solvation models, including the Hawkins-Cramer-Truhlar41, GBn42, and ABSINTH43 models, thus demonstrating its great flexibility.
gb = CustomGBForce( ) gb.addPerParticleParameter(“q”) gb.addPerParticleParameter(“radius”) gb.addPerParticleParameter(“scale”) gb.addGlobalParameter(“solventDielectric”, 78.3) gb.addGlobalParameter(“soluteDielectric”, 1.0) gb.addComputedValue(“I”, “““step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(1/U^2-1/L^2)*(r- sr2*sr2/r)+0.5*log(L/U)/r+C); U=r+sr2; C=2*(1/or1-1/L)*step(sr2-r-or1); L=max(or1, D); D=abs(r-sr2); sr2 = scale2*or2; or1 = radius1-0.009; or2 = radius2-0.009”””, CustomGBForce.ParticlePairNoExclusions) gb.addComputedValue(“B”, “““1/(1/or-tanh(1*psi-0.8*psi^2+4.85*psi^3)/radius); psi=I*or; or=radius-0.009”””, CustomGBForce.SingleParticle) gb.addEnergyTerm(“28.3919551*(radius+0.14)^2*(radius/B)^6- 0.5*138.935456*(1/soluteDielectric-1/solventDielectric)*q^2/B”, CustomGBForce.SingleParticle) gb.addEnergyTerm(“““−138.935456*(1/soluteDielectric- 1/solventDielectric)*q1*q2/f; f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))”””, CustomGBForce.ParticlePairNoExclusions)
Just as custom forces allow arbitrary new interactions to be defined, custom integrators can implement entirely new integration algorithms. An algorithm consists of a series of steps, each defined by a user supplied mathematical expression. A step may calculate a single value, or be executed separately for each degree of freedom. It also may involve arbitrary, user defined variables. This system is flexible enough to support a wide range of deterministic and stochastic algorithms including multiple time step integrators, Metropolized integrators, Generalized Langevin Equation based integrators, and even some simple Monte Carlo methods.
Listing 5 presents a simple example of using a custom integrator to implement the velocity Verlet algorithm. It creates a new CustomIntegrator, specifies a time step of 0.001 ps, and adds three steps to the algorithm to be calculated for each degree of freedom:
Notice that the force F appearing in the first step is different from the force appearing in the third step, since the atom position x has changed in between. The user does not need to specify when the force should be recomputed. OpenMM recognizes when anything changes that might affect it and recomputes it automatically.
integrator = CustomIntegrator(0.001) integrator.addComputePerDof(“v”, “v+0.5*dt*f/m”) integrator.addComputePerDof(“x”, “x+dt*v”) integrator.addComputePerDof(“v”, “v+0.5*dt*f/m”)
In both cases the appropriate OpenMM routines are called within the time-step loop of each program to calculate the energies and forces and perform the time-step integration. The setup (reading of the input files, mapping the force-field parameters to the system particles, etc.) and any run-time analysis are executed by the molecular dynamics program.
The main steps involved in embedding OpenMM into these applications are identical:
In both cases only the source file containing the main time-step loop of the molecular dynamics programs required modification. The changes included one-line calls to execute the four steps outlined above. A new source file was included that implemented these four steps. The last three steps consist of only a few lines of code; for both cases, the setup step is the most involved. In the case of Gromacs, the code in the new, added file used the C++ OpenMM API, whereas the C OpenMM API was used for TINKER which is written in Fortran. By making minimal changes to the original molecular dynamics code and isolating the OpenMM-specific code, the risk of introducing errors into the molecular dynamics program is reduced and maintenance of the program is much more straightforward.
Three types of validation of OpenMM are performed with each major release: unit tests, system tests, and direct comparison with other applications if possible. In the following, a brief description of each type of validation is given.
Unit tests are provided for each major force and integrator class and other auxiliary functions (e.g., the random number generator) and can be run by users to check that the code is working properly for their particular hardware setup. This ability is critical since GPU software and hardware is rapidly evolving and incompatibilities or unexpected results on next generation hardware are not uncommon. The unit tests exercise the basic functionality of each class to probe for problems; a separate collection of unit tests for each force and integrator is available for each of the different platforms. Most of these tests use simple model systems comprised of a small number of particles. However for the nonbonded force and implicit solvent models, some of the unit tests include relatively large systems with ‘extreme’ parameter values to help isolate issues.
Whereas unit tests validate the operation of small pieces of code in isolation, system tests are designed to test the library as a whole. They are performed on biomolecules (proteins, RNA, and DNA) typical of what users actually simulate, and they use realistic force fields. Systems are included that use both implicit and explicit solvent. System sizes range from 75 up to 173,181 atoms. Detailed results for these tests are provided in the OpenMM Users Guide, which the reader should consult for more information. Here we present a summary of the most important results.
The types of tests included in this category are listed below:
A 1 ns simulation is performed. For Verlet integration, it verifies that the magnitude of fluctuations in the total energy are below a cutoff. For Langevin integration, it verifies that the average instantaneous temperature equals the expected value to within a tolerance. In both cases, it also monitors all distance constrains and verifies that none is violated at any point during the simulation.
Direct comparison between energies and forces computed using OpenMM and another application.
For the OpenMM core functions, the comparisons are made with Gromacs. For the OpenMM AMOEBA plugin, the comparisons are made with TINKER. Table 3 shows the average relative difference in forces between the OpenMM Reference platform and Gromacs 4.5.3 (compiled in single precision mode) over a collection of molecules. Numbers are not shown for PME, since differences between how OpenMM and Gromacs handle cutoffs on nonbonded interactions (charge groups, shifting functions, etc.) complicate the comparison and require a more sophisticated analysis.
To benchmark the performance of OpenMM, we used the dihydrofolate reductase (DHFR) models taken from the Joint Amber/Charmm benchmark. This is a 159 residue protein with 2489 atoms. The version used for explicit solvent simulations included 7023 TIP3P water molecules, giving a total of 23,558 atoms.
The first set of simulations used the AMBER99SB force field.46 We used three different methods to model the effect of solvent:
We also used two different sets of constraints within the protein, allowing for different time steps:
Water molecules were fully rigid in all cases, and the constraint error tolerance was set to 10−4. (For H-bonds constraints, the error tolerance has a negligible effect on performance. For H-angles, reducing it to 10−5 causes a measurable but still quite small decrease in performance.) All simulations used an NVT ensemble, with a Langevin thermostat to maintain a temperature of 300K. The results for both the OpenCL and CUDA platforms are shown in Table 1. The CUDA simulations were run on a single Nvidia GTX 580 GPU. For OpenCL, simulations were run on both the GTX 580 and on an AMD Radeon HD 7970.
The OpenCL platform can parallelize an explicit solvent simulation across multiple GPUs. We therefore repeated each of the above benchmarks using 1 to 4 Nvidia C2070 GPUs in parallel. The results are shown in Table 2.
The scaling with number of GPUs is much less than linear. Using up to three GPUs produces a significant speedup, but there is little benefit beyond that. For PME, using four GPUs is actually slightly slower than using three. This is primarily due to the cost of transferring data over the high latency PCIe bus.
To demonstrate the performance of custom forces, we repeated the Explicit-RF, H-bonds benchmark for the OpenCL platform, but instead of using OpenMM’s built-in implementation of Coulomb and Lennard-Jones forces, we reimplemented them using the CustomNonbondedForce class. Running on a single Nvidia C2070 it achieves 23.0 ns/day, compared to 25.9 ns/day for the built-in forces, a speed penalty of only 11%. This shows the great value of custom forces: completely arbitrary functional forms, defined with only a few lines of code, can nearly match the performance of a highly optimized, hand-written implementation.
Benchmarks for the AMOEBA plugin were carried out for villin, a small protein with 596 atoms, and DHFR in both explicit and implicit solvent on an Nvidia GTX 580 GPU. The error tolerance for the calculation of the induced dipoles was 0.01 and the time step was 1 femtosecond. The results are shown in Table 3.
To measure energy conservation, we performed a 10 ns simulation of ubiquitin47, a 1231 atom protein in OBC implicit solvent. The simulation used the OpenCL Platform, a 1 fs time step, no constraints, and no cutoff on the nonbonded interactions, and had an average temperature of 298.6 K. A linear regression to the curve of energy versus time gives an overall drift rate of 2.5 kJ/mole/ns, or 0.00027 kT/ns/dof in the more commonly reported units. Note, however, that these units are not very meaningful since energy drift does not actually scale linearly with time, temperature, or system size. Furthermore, a simple linear regression is not sufficient to fully characterize energy conservation. A thorough discussion of these issues and a more detailed study of energy conservation in OpenMM is published elsewhere.48
OpenMM is a software toolkit that addresses an unmet need in the molecular simulation community. It is designed to be modular, extensible, reusable, and hardware independent. This combination of traits makes it uniquely valuable to several different types of users. Application developers have instant access to a large library of robust, high performance algorithms. Researchers developing new methods can use OpenMM’s plugin architecture and custom forces to implement and test new algorithms. Once written, those algorithms can quickly be incorporated into any program that uses OpenMM. For chemists and biologists interested in running simulations, the Python based application layer is a uniquely flexible and powerful environment offering complete control over nearly every aspect of the simulation protocol.
This work was supported by Simbios via the NIH Roadmap for Medical Research Grant U54 GM072970, and by NIH grant R01-GM062868. JDC acknowledges support from a QB3-Berkeley Distinguished Postdoctoral Fellowship.