Reference implementation energies [WIP]

Objective

Some scientists want to use our force fields but cannot or do not wish to use our software stack. A concern in these cases is that the force fields may be used in ways that are not consistent with how they are intended to be used. (Some may be relatively minor, like using a different tail correction, but other details such as the cutoff distance, combination rule, etc. can have more significant impacts.) A concrete way to communicate the intended use of a force field is to distribute it alongside a specification that describes in detail everything needed to reproduce the potential energy function - the SMIRNOFF spec does this (caveats below). Yet a better way is to also distribute an implementation of this - the OpenFF Toolkit does this (caveat(s) below). In principle, any user can prepare their system and pass it through the toolkit alongside a force field to get out objects (probably an openmm.System and openmm.Topology) that can be evaluated to get reference energies that can be used to validate a different implementation. This can be a fair amount of work for a user, and many probably won’t bother going through this process if it requires installing software they don’t intend to use and/or aren’t familiar with. These energies may also be subject to modifications by the user, intentional or not, resulting from the rabbit hole of different input settings that can impact energy evaluations.

We should make this process easier by curating a data set of systems, described below, and what we think the energy of each should be (including both the potential energy and per-term contributions). This would make it easier for users to validate their implementations and should mostly take control of the details out of their hands. This also provides internal value; it can uncover errors in our own implementations and serve as a natural data set to run regression tests against. The toolkit already runs some regression tests against a non-exhaustive set of molecules in an effort to safeguard against critical bugs, but leaves many use cases uncovered (and could be incorrect itself). Something similar is in the Interchange test suite; tracking reference energies from a canonical data set will reduce duplication of effort and increase the quality and reliability of these tests.

Curation

This data set should probably be distributed as a version-controlled set of files and scripts to process those files, supplemented with a table describing everything in detail. I think the result could look something like this, with potentially many (10s to 100s) of rows:

• Each record should specify, one way or another:

  • Force field (i.e. specific file/DOI)

  • Periodicity

  • Periodic box vectors, if any

  • Atomic coordinates

    • Specified in a file, not generated on the fly

  • Bond graph/connectivity

• Each section of the SMIRNOFF spec should be covered. This includes

  • Multi-term torsions

  • Improper torsions

  • WBO interpolated torsions and bonds

  • AM1BCC charges

  • Library charges

  • Custom charge increments

  • Constrained and non-constrained force fields

  • Anything specific to biopolymers?

  • GBSA?

  • Virtual sites

• Some vaguely broad amount of chemistry should be covered

  • Single ligand(s) in vacuum

    • Molecules with charged (zwitterionic?) groups

    • Molecules that are perceived differently by different aromaticity models

    • Molecules with non-planar impropers (i.e. pyrimidal nitrogens)

  • Molecule dimers in vacuum

  • Box of organic liquids

  • Box of water?

  • Protein(s) in vacuum

  • Protein(s) in saltwater

  • Ligand(s) bound to protein(s)?

  • Other biomolecules

Issues

The toolkit is an imperfect reference implementation

Because the toolkit’s implementation of SMIRNOFF force fields is (necessarily, for historical reasons) OpenMM-centric, there are some combinations of inputs in which it does not strictly follow the specification. Unfortunately, one of them is a highly-trafficked path; vdW interactions are computed with no cutoff if the input topology does not have box vectors specified (as is commonly the case when evaluating single-molecule systems in the gas phase, i.e. to compare to QM). There may continue to be good reasons to allow the toolkit to slightly disobey the SMIRNOFF specification, but it should be made clear the cases in which that is true and what should be used as the strictest reference (i.e. this project).

Stochasticity of AM1-BCC partial charges

It is well-known that the many subtleties in the generation of AM1-BCC partial charge result in small inconsistencies in the results. This makes it difficult to, given a molecular graph with coordinates, ensure to high precision that the electrostatic energy will be any approximate value and impossible to ensure a precise value.

PME is an approximation

It is also well-known that PME, while accurate enough for the vast majority of use cases, necessarily introduces small approximations. Different simulation engines also expose different PME parameters, which makes it practically impossible to ensure high precision comparisons between implementations or engines, even if they were self-consistent across different hardware.

For purposes of ensuring parameters in a force field are applied, other methods of computing the electrostatic energy may be sufficient. The simplest method would be applying a hard cut-off with no long-range interactions, tail corrections, or reaction-field treatment. This should enable comparisons, but it unacceptable for the vast majority of use cases that our force fields find themselves in. Also, not all engines offer this method.

The switching function is never explicitly specified

The SMIRNOFF spec describes (in somewhat unclear terms) the pair distance at which a switching function should be applied. It does not, however, describe what this switching function actually is. Engines each use slightly different forms of a switching function and/or make it unclear what is actually used, if it is supported.