OpenMM 7: Rapid development of high performance algorithms for molecular dynamics — TEI-tagged view

Peter EastmanpersName, Jason SwailspersName, John D. ChoderapersName, Robert T. McGibbonpersName, Yutong ZhaopersName, Kyle A. BeauchamppersName, Lee-Ping WangpersName, Andrew C. SimmonettpersName, Matthew P. HarriganpersName, Chaya D. SternpersName, Rafal P. WiewiorapersName, Bernard R. BrookspersName, Vijay S. PandepersName (2017date). “OpenMM 7: Rapid development of high performance algorithms for molecular dynamics.” PLOS Computational Biologytitle·j.

Object of study — OpenMMterm 7 architecture, extensibilityterm mechanisms, hardware portability, advanced molecular models, precision modesterm·var, and performance

Research question — how to make new molecular-dynamics methods easy to implement, portable across hardware, performant, and independently distributable

Contribution — a layered molecular-dynamics toolkit combining Python scripting, expression-defined forces and integratorsterm, runtime pluginsterm, and multiple input pipelinesterm

Full text — tagged

The reference list is not part of the text we hold for this paper. Our copy was extracted from the publisher's page and the bibliography did not survive that extraction, so the article below ends before its references. Nothing has been removed by us.

OpenMMterm 7: Rapid development of high performance algorithms for molecular dynamics Abstract OpenMMterm is a molecular dynamics simulation toolkitterm with a unique focus on extensibilityterm. It allows users to easily add new features, including forces with novel functional formsterm, new integration algorithmsterm, and new simulation protocolsterm. Those features automatically work on all supported hardware typesterm·var (including both CPUsterm and GPUsterm) and perform well on all of them. In many cases they require minimal coding, just a mathematical description of the desired function. They also require no modification to OpenMMterm itself and can be distributed independently of OpenMMterm. This makes it an ideal tool for researchers developing new simulation methods, and also allows those new methods to be immediately available to the larger community.

This is a PLOS Computational Biology Software paper.

Introduction Background Molecular dynamics simulationterm is a rapidly advancing field. Many aspects of it are subjects of current research and development. Some of the more important examples include the development of new force fieldsterm [ 1 , 2 ], sometimes involving novel functional forms for the interactions [ 3 , 4 ]; new integration algorithmsterm [ 5 – 7 ]; new sampling methodsterm [ 8 – 10 ]; and support for new types of hardware [ 11 – 13 ].

There are many popular software packages for conducting molecular dynamics simulationsterm. They vary considerably in their capabilities and feature sets. This is especially true when it comes to cutting edge, recently developed simulation techniques. The inventor of a new method will typically implement it in a single package, whichever one they are most comfortable working with. From that point, it may take years to be implemented in other packages, depending on the interests and priorities of the development team behind each one. In many cases, it may never get implemented. Even the initial implementation may not be accepted into an official release of the package it was created in. Or it may have limited usefulness, for example because it executes slowly or cannot be used on advanced hardware such as graphics processing units (GPUsterm).

The main reason for this problem is that most molecular dynamics packages were not designed with extensibilityterm in mind. Adding new features, even very simple ones, is often labor intensive and requires a deep understanding of the code. Once a prototype implementation is complete, it may be even more difficult to turn that into a well optimized version that works on all hardware typesterm·var. In most cases the simulation enginesterm are also monolithic, so the only way to add features to them is to directly modify their source code. There is no plugin interfaceterm or other mechanism for allowing new features to be implemented and distributed independently. This turns the core development teamterm into gatekeepers, restricting what features can be added to the package.

A complete molecular dynamics package is, of course, much more than just a simulation engineterm. Each one typically has its own collection of tools for preparing molecular systems to simulate, its own file formats, and sometimes even its own force fieldsterm. This makes it difficult for users to switch back and forth between them, or to combine features from different packages.

OpenMMterm OpenMMterm is a molecular dynamics package designed to address these problems. It began as simply a library for performing certain types of calculationsterm on GPUsterm, but in recent versions has grown into a complete simulation package with unique and powerful features. This article describes OpenMMterm 7.0, which is the latest release at the time of writing. An earlier version (OpenMMterm 4.1) was described in a previous publication [ 14 ]. This article focuses primarily on what has changed since that version, but for completeness there is some overlap between the two.

OpenMMterm is based on a layered architectureterm which (see Fig 1 ), to the best of our knowledge, is unique among molecular dynamics packages. This allows it to be used in several different ways by users with varying needs and interests. Depending on how a particular user chooses to interact with it, OpenMMterm can act as: A high-performance libraryterm , callable from other programs, for performing a wide range of calculationsterm used in molecular modelling and simulation on a range of advanced hardware platforms (both CPUsterm and GPUsterm).

A domain specific languageterm for easily implementing new algorithms for molecular modelling and simulation.

A complete packageterm for running molecular simulations.

Architecture of OpenMMterm.

Extensibilityterm is built into every layer of the architecture as a fundamental design goal. A guiding philosophy is that users should be able to implement new features as easily as possible, by writing as little code as possible. Those features should then work on all types of hardware, including both CPUsterm and GPUsterm, and have good performance on all of them. Finally, the developer of a feature should be able to package and distribute their code independently, without needing the approval or participation of the core OpenMMterm development team.

The highest layer of the architecture is based on the Python scripting languageterm. Users can easily extend it by writing their own Python code to implement the algorithms of their choice. A wide range of simulation protocols, sampling methods, etc. can be implemented in this way, often with only a few lines of code.

The next layer down defines the calculationsterm that are tied together through Python scripting. This layer includes many classes for creating “custom” forces and integratorsterm. These classes provide a simple but powerful mechanism for extensibilityterm. The user provides one or more mathematical expressionsterm to describe the calculationterm to be done. For example, they might give an expression for the interaction energyterm of a pair of particles as a function of the distance between them. The expression is parsed and analyzed, and just-in-time compilationterm is used to generate an efficient implementationterm of the code for calculating that interaction [ 15 ]. This allows users to easily define a huge variety of interactions and integration algorithms. They can then be used on any supported type of hardware, and involve little or no loss in performance.

At the lowest layer, OpenMMterm is based on a pluginterm mechanism. Calculationsterm are defined by “computational kernelsterm”. A pluginterm may define new kernels for doing new types of calculationsterm, or alternatively it may provide new implementations of existing kernels, for example to support a new type of hardware. Pluginsterm are dynamically discovered and loaded at runtime. Each one is packaged as a file that can be distributed separately from the rest of OpenMMterm and installed by any user.

Another unique feature of OpenMMterm is its support for multiple input pipelinesterm. Before a molecular system can be simulated, it first must be modelled. This is sometimes a complex process involving such steps as combining multiple molecules into a single file, building missing loops, selecting a force fieldterm, and parametrizing small molecules. Typically, each simulation package provides its own tools for doing this. They often differ in significant ways, such as what force fieldsterm are available.

OpenMMterm does include modelling tools, but it also can directly read the file formats used by Amberterm [ 16 ], CHARMMterm [ 17 ], Gromacsterm [ 18 ], and Desmondterm [ 19 ]. A user can prepare their system with the tools from any of those packages, or with other tools that are designed to work with them, then simulate it in OpenMMterm. This gives great flexibility, since the user can use whatever tools are best suited to the system they want to simulate. It also lets OpenMMterm easily fit into their existing workflow. A user who is accustomed to a particular tool can continue to use it, but still run their simulation in OpenMMterm.

Design and implementation Extensibilityterm Python scripting The highest level of the architecture consists of a set of Python classes and functions. They may be chained together to create simple scripts that run simulations, or more complicated ones that implement a variety of advanced algorithms. These are some of the functions provided by OpenMMterm that may be used by Python scripts: Reading input files, including standard formats like PDB or PDBx/mmCIF, as well as the proprietary formats used by applications such as Amberterm, CHARMMterm, Gromacsterm, and Desmondterm.

Editing molecular models, such as by combining molecules together, adding or deleting atoms, building solvent boxes, etc.

Defining the forces acting on a molecular system, either by specifying them explicitly or by loading a force fieldterm definition from a file.

Computing forces and energiesterm.

Running simulations.

Outputting results.

An example of a script to run a simulation is shown in Listing 1. It loads a PDB file, models the forces with the AMBER99SB-ILDN force fieldterm [ 20 ] and TIP-3P water model [ 21 ], performs a local energy minimization to eliminate clashes, and then simulates 1 million steps of Verlet dynamics. Every 1000 steps, it writes the current structure to a DCD file, and the current time, potential energy, and temperature to a log file.

Listing 1: A Python script executing a complete molecular simulation from a PDB file.

from simtk.openmmterm.app import * from simtk.openmmterm import * from simtk.unit import * pdb = PDBFile('input.pdb') forcefield = ForceField('amber99sbildn.xml', 'tip3p.xml') system = forcefield.createSystem(pdb.topology, nonbondedMethod = PME, nonbondedCutoff = 1*nanometer, constraints = HBonds) integrator = VerletIntegrator(0.002*picoseconds) simulation = Simulation(pdb.topology, system, integrator) simulation.context.setPositions(pdb.positions) simulation.context.setVelocitiesToTemperature(300*kelvin) simulation.minimizeEnergy() simulation.reporters.append(DCDReporter('output.pdb', 1000)) simulation.reporters.append(StateDataReporter('output.log', 1000, time = True, potentialEnergy = True, temperature = True)) simulation.step(1000000) This script runs a simple simulation, much like one might run in any molecular dynamics package. More sophisticated or exotic algorithms and protocols can be implemented in exactly the same way. Here are some examples of features or applications that use or extend OpenMMterm’s Python scripting features: A class for simulated tempering [ 10 ], an accelerated sampling method that varies the temperature of a simulation to accelerate barrier crossings. The entire algorithm was implemented in roughly 200 lines of code.

PDBFixerterm, an application for cleaning up molecular models in preparation for simulating them. It includes such features as building missing loops, replacing nonstandard amino acids with standard ones, adding hydrogens, and building solvent boxes. By using the features provided by OpenMMterm, all of these algorithms were implemented in only about 1000 lines of code.

YANKterm [ 22 ] a sophisticated application and toolkit for alchemical free energy calculationsterm. It implements Hamiltonian exchange molecular dynamics simulationsterm to efficiently sample multiple alchemical states, and utilizes the “custom” forces provided by OpenMMterm to allow exploration of many different alchemical intermediate functional forms.

Although these tools are written in Python, all expensive calculationsterm are done by OpenMMterm and take full advantage of the available hardware, including GPUsterm and multicore CPUsterm. Because they interact with OpenMMterm only through well-defined public interfaces, they can be packaged and distributed independently. No changes to OpenMMterm itself are required to use them.

A set of advanced examples is included in the supporting information. They demonstrate more complex simulation techniques, and show how to use OpenMMterm in combination with other programs.

Custom forcesterm In addition to the standard forces provided by OpenMMterm (such as Lennard-Jones forces, PME and reaction field electrostatics, and generalized Born models), custom forcesterm are a mechanism for creating interactionsterm between particles with entirely novel functional forms. There are many different custom forceterm classes, each supporting a particular category of interactions. They are listed in Table 1 . Since OpenMMterm 4.1 was described in an earlier publication, several new custom forceterm classes have been added, including CustomCompoundBondForce (added in OpenMMterm 5.0), CustomManyParticleForce (added in OpenMMterm 6.2), and CustomCentroidBondForce (added in OpenMMterm 7.0).

Custom forcesterm supported by OpenMMterm 7.0.

In each case, the user provides an algebraic expression for the interaction energyterm as a function of the relevant variables. OpenMMterm analytically differentiates the expression to determine the corresponding force, then uses just-in-time compilationterm to generate machine codeterm for efficiently computing the force and energy on the current hardware (either CPUterm or GPUterm).

As an example, Listing 2 defines a harmonic restraint that can be applied to the angles formed by triplets of atoms. It specifies that the energy of each triplet is given by k ( θ - θ 0 ) 2 . It also specifies that k and θ 0 are per-angle parameters: each triplet can have different values for them.

Listing 2: Implementation of a harmonic angle restraint using a CustomAngleForce.

force = CustomAngleForce('k*(theta-theta0)^2') force.addPerAngleParameter('k') force.addPerAngleParameter('theta0') Custom forcesterm are designed to achieve several goals that usually conflict with each other. First, it should be exceptionally easy to implement completely new functional forms for interactions. As seen in Listing 2, it often requires no more than a few lines of Python code. Second, a single implementation should work on all types of hardware. The exact same code can be used whether the program is being run on a CPUterm or GPUterm. Third, the user should not need to sacrifice performance. Because the expression is converted to machine codeterm before the simulation is run, there often is little or no difference in speed between a custom forceterm and a hand-written implementation of the same interaction.

Custom integratorsterm Just as custom forcesterm allow users to implement novel interactions, custom integratorsterm allow them to implement novel integration algorithms. The algorithm is defined by a sequence of operations, each defining a calculationterm to be done. Various types of operations are supported. Examples include: Evaluating a mathematical expressionterm for each degree of freedom, then assigning the result to a variable for each one.

Evaluating a mathematical expressionterm once and assigning the result to a global variable.

Summing an expression over all degrees of freedom and assigning the result to a global variable.

Applying constraints to positions or velocities.

In the simplest case, all operations are executed in order to take a single integration time step. In addition, OpenMMterm 7.0 added support for more complex flow control through if and while blocks.

Listing 3 shows Python code that uses a custom integratorterm to implement the leapfrog Verlet algorithm. The function addPerDofVariable() defines a new variable that has a different value for each degree of freedom. The function addComputePerDof() defines a calculationterm to be performed independently for each degree of freedom. In the absence of constraints, each time step of this algorithm consists of the operations v ← v + d t · f / m x ← x + d t · v where x is the position at time t , v is the velocity at time t - dt /2, dt is the step size, f is the force, and m is the particle mass. When constraints are present, the positions must then be adjusted to satisfy them, and finally the velocities are recalculated as v ← ( x − x 0 ) / d t where x 0 is the position at the start of the step.

Listing 3: Leapfrog Verlet algorithm implemented as a CustomIntegrator.

integrator = CustomIntegrator(dt) integrator.addPerDofVariable('x0', 0) integrator.addComputePerDof('x0', 'x') integrator.addComputePerDof('v', 'v+dt*f/m') integrator.addComputePerDof('x', 'x+dt*v') integrator.addConstrainPositions() integrator.addComputePerDof('v', '(x-x0)/dt') Far more complicated and sophisticated algorithms can be implemented in the same way. Here are some examples of integrators that have been created with this mechanism.

The rRESPA multiple time step integration algorithm [ 23 ].

The aMD accelerated sampling algorithm [ 9 ].

Metropolis-Hastings Monte Carlo [ 24 ] with Gaussian displacement proposals.

Hybrid Monte Carlo and variants, such as Generalized hybrid Monte Carlo (GHMC) [ 25 ], a Metropolized form of Langevin dynamics.

Nonequilibrium candidate Monte Carlo (NCMC) [ 26 ], where an external field is changed during the course of dynamics and the resulting nonequilibrium proposal accepted or rejected to preserve the equilibrium distribution.

As with custom forcesterm, a single implementation works on all types of hardware. Because just-in-time compilationterm is used to generate efficient machine codeterm for the algorithm, there usually is little or no performance cost relative to using hand-written GPUterm code.

Pluginsterm The lowest layer of the OpenMMterm architecture is based around pluginsterm. This allows new features to be packaged as libraries, distributed independently, and loaded dynamically at runtime. For example, a pluginterm can implement a new type of interaction or a new integration algorithmterm, or it can add support for a new type of hardware. In fact, many of the core features of OpenMMterm are actually implemented as pluginsterm, including its implementations of the AMOEBA force fieldterm [ 27 ], ring polymer molecular dynamics (RPMD) [ 28 ], and polarizable Drude particlesterm [ 29 ].

This provides nearly unlimited extensibilityterm, allowing users to implement any feature they might want. Writing a pluginterm involves far more work than the other extensibilityterm features described above. For example, it is up to the user to write whatever code is necessary to make it work on each type of hardware, such as CUDA or OpenCL code for GPUsterm. When possible, it is therefore usually preferable to use one of the other mechanisms. Nevertheless, pluginsterm are an important option when extreme extensibilityterm and performance is needed. Like the other mechanisms, they allow a developer to create an extension and distribute it directly to users. No modifications to OpenMMterm itself are needed.

Advanced features OpenMMterm has many other features beyond those discussed above, some of which are themselves unique or noteworthy. The following are some of the more significant ones.

AMOEBA OpenMMterm has an implementation of the AMOEBA polarizable force fieldterm which is, to the best of our knowledge, the fastest available in any code [ 3 ]. AMOEBA is designed to transcend the limitations of conventional point charge force fieldsterm and achieve much higher accuracy in force and energy computations. It uses two main mechanisms to achieve this. First, instead of approximating atoms as point charges, it assigns each one a multipole momentterm up to the level of quadrupoles. Second, it explicitly models atomic polarization by assigning an induced dipoleterm to each atom. Because the induced dipolesterm interact with each other, they must be computed at each time step using an iterative self-consistent field calculationterm. Both of these features make AMOEBA far slower to simulate than conventional force fieldsterm.

Much research has been done recently on ways to reduce this cost, and new versions of OpenMMterm have incorporated several of the most recent algorithms. Interactions between multipoles are computed using spherical harmonics in a quasi-internal coordinate system [ 30 , 31 ] (added in OpenMMterm 7.0). The iterative solver for induced dipolesterm uses the Direct Inversion in the Iterative Subspace (DIIS) algorithm [ 32 ] (added in OpenMMterm 6.1). Alternatively, it can use the recently developed extrapolated polarization approximationterm [ 33 ] (added in OpenMMterm 7.0). In this method, only a few iterations are performed, and then an analytic approximation is used to extrapolate to the limit of infinite iterations. This can give a large improvement in speedterm with only a very small loss in accuracy.

Drude particlesterm A new feature introduced in OpenMMterm 5.2 is support for Drude particlesterm [ 29 ] as an alternative way of modelling polarizability. In this method, each polarizable atom is modelled as a pair of charges connected by an anisotropic harmonic force. When an electric field is applied, the two particles are displaced from each other, creating a dipole moment. The strength of the force connecting them determines the atomic polarizability. The particle positions can be determined using a self-consistent field calculationterm or, more commonly, a dual-thermostat Langevin integratorterm that couples the center of mass of each pair to a high temperature heat bath (e.g. 300K), but the internal motion of each pair to a low temperature heat bath (e.g. 1K).

Polarizable force fieldsterm based on Drude particlesterm are included with OpenMMterm. This includes the SWM4-NDP water model [ 34 ], and the CHARMMterm polarizable force fieldterm for proteins [ 4 ]. They aim to incorporate some of the same physical effects as AMOEBA at a lower computational cost.

Virtual sitesterm Virtual sitesterm are interaction sitesterm within a molecule whose positions are not integrated directly. Instead, they are calculated at each time step based on the positions of other particles. They are often used to provide a more detailed charge distribution than would be possible using only a single point charge for each atom. For example, they appear in many multisite water models (such as TIP-4P and TIP-5P), and also in the CHARMMterm polarizable protein force fieldterm.

There are a multitude of possible ways a virtual siteterm position can be specified based on the positions of other mobile atoms. Typically, a simulation package will provide a limited choice of rules, covering only those cases needed for the particular force fieldsterm that package supports. For example, Gromacsterm 5 offers a choice of four methods for calculating a virtual siteterm position based on the positions of three atoms. Each one covers one very specific case that occurs in a supported force fieldterm. One can easily imagine other cases that would be impossible to construct with any of the current rules, and would therefore require adding a fifth rule.

OpenMMterm also offers a few specialized rules for positioning virtual sitesterm, but in addition, it has a very general method designed to cover all cases that are ever likely to occur in which a virtual siteterm depends on the positions of three atoms. This method, added in OpenMMterm 6.1, can reproduce all four of the rules provided by Gromacsterm, as well as supporting many other situations they could not.

In this method, three vectors are first calculated as weighted averages of the positions of the three atoms: o = w 1 o r 1 + w 2 o r 2 + w 3 o r 3 d x = w 1 x r 1 + w 2 x r 2 + w 3 x r 3 d y = w 1 y r 1 + w 2 y r 2 + w 3 y r 3 where r 1 , r 2 , and r 3 are the atom positions, and the coefficients are user-defined. They are then used to construct a set of orthonormal coordinate axes ( x ^ , y ^ , z ^ ): d z = d x × d y x ^ = d x / | d x | z ^ = d z / | d z | y ^ = z ^ × x ^ Finally, the virtual siteterm position is set to an arbitrary user-defined location within this coordinate system: r = o + p 1 x ^ + p 2 y ^ + p 3 z ^ This method is another example of how flexibility and extensibilityterm are core design goals of OpenMMterm. Instead of supporting only a limited set of specialized virtual siteterm types, it tries to provide a very general type that can cover as wide a range of cases as possible, thus giving maximum flexibility to users in designing their models and force fieldsterm.

Triclinic periodic boxes Earlier versions of OpenMMterm supported only rectangular periodic boxes. In OpenMMterm 6.3 it was extended to support triclinic boxesterm as well: ones formed by combinations of three arbitrary lattice vectors. It can be shown that this formulation is extremely general; all standard periodic box shapes, including the popular rhombic dodecahedronterm and truncated octahedron, can be represented as triclinic boxesterm [ 35 ].

This feature serves two important functions. First, it allows one to simulate crystals, which very often have non-rectangular unit cells. Second, it is useful when simulating freely rotating molecules in solvent. One needs to include a certain amount of padding around the molecule to ensure that no two periodic copies ever come too close together. Because the molecule can freely rotate, the same padding is required along all directions, so one wants the periodic box to be as close as possible to spherical. For a given padding distance, the rhombic dodecahedronterm has only about 71% the volume of a rectangular unit cell. It therefore requires less solvent and reduces the cost of the simulationterm.

Multiple precision modesterm·var Many aspects of a molecular dynamics code involve tradeoffs between speed and accuracy. This is especially true when executing on a GPUterm, since they often have very poor double precisionterm performance. To optimize execution speed, it is preferable to use single precisionterm whenever possible, resorting to double precisionterm only when absolutely necessary. Unfortunately, there is no unique standard for when it is "necessary". The minimum acceptable level of error can vary widely depending on the details of a simulation and the type of information one wishes to obtain from it.

OpenMMterm 5.0 introduced support for multiple precision modesterm·var. When running on a GPUterm, the user has a choice of three modes: Single Precisionterm : Nearly all calculationsterm are done in single precisionterm. Double is used in only a handful of places where it has negligible impact on performance and is most important for accuracy.

Mixed Precisionterm : Forces are computed in single precisionterm, but integration and energy accumulation are done in double precisionterm. This gives a large improvement in the accuracy of some quantitiesterm, while only having a small impact in performance.

Double Precisionterm : All calculationsterm are done in double precisionterm. This gives the best accuracy, but often has a very large effect on performance.

The effects of the different modes are illustrated below in Results.

Regardless of the precision modeterm·var, forces are accumulated as 64 bit fixed point values. This improves accuracy when working in single or mixed precisionterm modes, and ensures that force accumulation is deterministic. It also allows force accumulation to be done with integer atomic operations, which substantially improves performance. OpenMMterm has used this method since version 4.0, released in January 2012. Since that time, it has found its way into other GPUterm accelerated MD codes, such as AMBERterm [ 36 ].

Results Performance To evaluate the speed of OpenMMterm, we benchmarked its performance with three molecular systems of varying size: Dihydrofolate reductase (DHFRterm), a 2489 atom protein solvated with 7023 water molecules to give a total of 23,558 atoms.

Abl kinase (ABL1term), a 4067 atom protein solvated with 13,692 water molecules to give a total of 45,143 atoms The mechanistic target of rapamycin (MTORterm), a 19,019 atom protein solvated with 56,733 water molecules to give a total of 189,218 atoms.

Benchmarks were run on the following types of hardware: An NVIDIA Titan X Pascal GPUterm.

An NVIDIA Tesla K80 GPUterm.

A 4 core, 3.5 GHz Intel Core i7-2700K CPUterm.

All GPUterm simulations used CUDA 7.5. The K80 consists of two independent GPUsterm on a single board. OpenMMterm can parallelize a single simulation across multiple GPUsterm, or alternatively run a different simulation on each one at the same time. We therefore included benchmarks using only one of the GPUsterm (thus leaving the other free for a different simulation), as well as ones using both GPUsterm for a single simulation.

In the discussion below, we summarize the most important parameters for each set of simulations. Full details can be found in the scripts used to run the simulations, which are included in the Supplemental Information.

Amberterm We first benchmarked the performance using the AMBER99SB-ILDN force fieldterm and TIP3P water model. All simulations used a Langevin integrator with a temperature of 300 K and a friction coefficient of 1 ps -1 . Long range Coulomb interactions were computed with the Particle Mesh Ewald (PME) method.

Simulations were run using integration time steps of both 2 fs and 5 fs. For the 2 fs simulations, covalent bonds involving a hydrogen atom were modelled as rigid constraints. For the 5 fs simulations, all covalent bonds were modelled as rigid constraints and hydrogen mass repartitioning was used to increase the mass of hydrogen atoms to 4 amu (while decreasing the masses of the atoms they were bonded to so as to keep the total system mass constant). In all cases, water molecules were kept rigid. All of the GPUterm simulations used single precisionterm.

The results are shown in Table 2 .

Benchmark results for various protein systems in explicit solvent simulated with PME.

Depending on the molecule and settingsterm·var, using two GPUsterm is anywhere from 22% to 62% faster than a single GPUterm, with the larger molecules generally having the higher speedups. If the goal is to run a single simulation as quickly as possible, using multiple GPUsterm is therefore quite useful. On the other hand, if the goal is to generate as much total simulation time as possible, it is more efficient to run a separate simulation on each one.

AMOEBA We next benchmarked performance using the AMOEBA2013 force fieldterm. The AMOEBA water model is designed to be flexible rather than rigid, which requires a smaller step size. We therefore used a rRESPA multiple time step integrator, in which bonded forces were evaluated every 1 fs and nonbonded forces every 2 fs. No degrees of freedom were constrained. As above, we used PME for long range Coulomb interactions and single precisionterm.

Simulations were run with two different methods of calculating the induced dipolesterm: Full mutual polarization with a tolerance of 10 −5 for the induced dipolesterm.

The extrapolated polarization approximationterm.

The CPUterm implementation of AMOEBA in OpenMMterm is not well optimized, so we only ran benchmarks on GPUsterm. Using multiple GPUsterm for a single simulation is not supported with AMOEBA. Because AMOEBA is a very expensive force fieldterm and is normally only used for modest sized systems, we only ran benchmarks for DHFRterm.

The results are shown in Table 3 .

Benchmark results for DHFRterm in explicit solvent using AMOEBA2013.

Effect of precision When running on a GPUterm, OpenMMterm gives a choice of three precision modesterm·var: single, mixed, and double. To measure the effect of this choice on performance, we repeated the 2 fs time step DHFRterm simulations in mixed and double modes. The results are shown in Table 4 .

Effect of precision model on performance.

The speed difference between single and mixed precisionterm is quite small, whereas double precisionterm is much slower. This is especially true on the Titan X, a GPUterm primarily targeted at consumers that has very poor double precisionterm performance. The Tesla K80, which is targeted at high performance computing, does much better, although there is still a large decrease in performanceterm. Overall, the Titan X is far faster in single or mixed precisionterm modes, while the K80 is faster in double precisionterm mode.

To see the benefits of higher precision, we performed additional simulations of DHFRterm. Because a thermostat tends to mask the effect of error, these simulations used a leapfrog Verlet integrator to simulate constant energy. All simulations used the AMBER99SB-ILDN force fieldterm, a 2 fs time step, rigid water, and constraints on bonds involving hydrogen.

Each simulation was 1 ns in length. The total energy was recorded every 1 ps, and a linear regression was used to estimate the rate of energy changeterm·var. Ten independent simulations were performed for each precision modeterm·var, giving ten estimates of the rate. Table 5 reports the mean and standard error of those ten rates for each mode.

Energy driftterm for different precision models for the DHFRterm explicit solvent system.

The energy driftterm in single precisionterm is more than two orders of magnitude larger than in mixed or double precisionterm. When accurate energy conservation is important, using mixed precisionterm has a very large benefit at low cost. The average drift rates in mixed and double precisionterm are not significantly different from each other, indicating that numeric precision is no longer the dominant source of error. In other cases, such as when using a smaller step size or when simulating a larger molecule, statistically significant differences between them might emerge.

Input pipelinesterm A key feature of OpenMMterm is its support for multiple input pipelinesterm. This allows users to prepare molecular systems with the tools of their choice, then simulate them in OpenMMterm. Support for Gromacsterm input files was added in OpenMMterm 5.1, Desmondterm file support was added in OpenMMterm 6.0, and CHARMMterm file support was added in OpenMMterm 6.1.

The code in Listing 1 began from a PDB file and force fieldterm definition, using those to construct a description of the molecular system. Listing 4 shows the changes needed to instead construct it from an Amberterm prmtop file, as created by the AmberTools suite of software. More complete examples of using Amberterm and CHARMMterm input files are included in the supporting information.

Listing 4: Loading a system from Amberterm prmtop/inpcrd files.

prmtop = AmberPrmtopFile('input.prmtop') system = prmtop.createSystem(nonbondedMethod = PME, nonbondedCutoff = 1*nanometer, constraints = HBonds) To validate the accuracy of the input pipelinesterm, we constructed systems using the setup tools from other packages, then loaded them into those packages and into OpenMMterm and compared the forces and energiesterm. We performed these tests on two systems: DHFRterm, a 159 residue protein, and 2KOC, a 14-mer hairpin RNA. Comparisons were made to Amberterm 16, Gromacsterm 4.6.5, and CHARMM-LITE c40b1. For Amberterm, we performed comparisons in both explicit solvent and OBC1 implicit solvent. For Gromacsterm and CHARMMterm, we compared only explicit solvent.

To create the Amberterm input files, ParmEd [ 37 ] was used to download PDB files 4M6J (DHFRterm) and 2KOC (RNA hairpin), which were then stripped of water molecules and adjusted to standard amino acids. In the case of 2KOC, the first model was used. Cofactors and phosphates were deleted. Amberterm prmtop and inpcrd files were created with LEaP from the AmberTools 16 distribution. For implicit solvent simulations, mbondi3 GB radii were used; for explicit solvent simulations, the system was solvated with TIP3P waters in an octahedral box with 15 Å of clearance, and 20 Na+ and 20 Cl- counterions were added.

To create the Gromacsterm input files, ParmEd was used to convert the Amberterm prmtop and inpcrd files into Gromacsterm top and gro files.

To create the CHARMMterm input files, CHARMM-GUI [ 38 ] was used to download the 4M6J and 2KOC PDB files. Crystallographic water molecules were deleted, and the system was solvated with a rectangular water box with 15 Å of padding. Default values were accepted for all other options, including replacing nonstandard amino acids and patching terminal residues.

All of the input files, as well as scripts needed to run the comparisons, are included in the Supplemental Information.

Results are shown for Amberterm in Tables 6 and 7 , for Gromacsterm in Tables 8 and 9 , and for CHARMMterm in Tables 10 and 11 . In all cases the agreementterm is excellent, with all energy components matching to at least four significant digits. In systems that use PME, the nonbonded energies have somewhat larger differences than other energy components. This is partly because of the larger magnitude of this interaction, and partly because of the fact that different applications compute nonbonded interactions in slightly different ways. For example, Amberterm uses 4th order splines for charge spreading, while OpenMMterm uses 5th order splines. Nonetheless, they both compute the forces and energy to similar overall accuracy.

Comparison of energy components, as calculated by Amberterm and OpenMMterm.

Comparison of forces as computed by Amberterm and OpenMMterm.

Comparison of energy components, as calculated by Gromacsterm and OpenMMterm.

Comparison of forces as computed by Gromacsterm and OpenMMterm.

Comparison of energy components, as calculated by CHARMMterm and OpenMMterm.

Comparison of forces as computed by CHARMMterm and OpenMMterm.

Availability and future directions OpenMMterm is available from http://openmmterm.org . Ongoing development is conducted through the Github community at https://github.com/pandegroup/openmmterm . Detailed instructions on how to compile it from source are found in the OpenMMterm User Guide at http://openmmterm.org/documentation.html .

Supporting information (ZIP) Click here for additional data file.

(ZIP) Click here for additional data file.

(DOCX) Click here for additional data file.