SMIRNOFF: The switching function is never explicitly specified

[mattwthompson]

The SMIRNOFF spec makes references to switching functions when discussing the pair distances at which it should be applied, but never states what switching function should be applied. Most simulations implement some form of a switching function, but no engines seem to agree on precisely how to do so, in particular the actual polynomial used (openforcefield/openff-toolkit#896 (comment)). It's not likely that these differences make a huge impact on most scientific results, but strict implementations will always be impossible while these details are unclear.

This causes a significant issue when validating any result relating to non-bonded interactions; mismatches in tail corrections tend to result in different non-bonded energies. This obfuscates whether differences are due to the values of the non-bonded parameters, the implementation of the switching function, other non-bonded settings (cutoff distance, switching distance, non-bonded method, etc.) or some combination thereof. This makes it difficult to compare implementations of SMIRNOFF force fields in different engines (openforcefield/openff-interchange#159).

Related #9

cc: @mrshirts

[davidlmobley]

Ugh. This is a painful one. Probably our reference implementation ought to be OpenMM so we probably ought to look at what we'd do when evaluating energies in OpenMM, then specify that in the spec. It'd be helpful though to cross-compare with GROMACS and see if we can choose something also supported by GROMACS.

[mrshirts]

I am still advocating for separating the parameters and short range functional forms for the force fields for how the nonbonded long range terms are handled. One will never get complete agreement between engines on how to handle these. I've tried. Even supposedly matching functional forms are only accurate up to one part in 10^5, which for a large system can be several kcals/mol. I tried to do it during the course of this: https://link.springer.com/content/pdf/10.1007/s10822-016-9977-1.pdf, it just couldn't be done better.

[mrshirts]

Functional form would be specified, and then there's some set of optional information - if you do a cutoff, you get this, if you do some type of LJ with switch you get this, etc.

I'm just noticing that like half of the issues have to do with trying to match long range interaction choices, and it would be much simpler if that was somehow and optional part of the spec, or intentionally or specifically undefined in the spec.

[mattwthompson]

Definitely explicitly and precisely the function form of the switching function would be an improvement in the specification, but it's not clear they could strictly be used outside of OpenMM if that's what's matched. (Of course, matching another engine's opinion doesn't really move the needle.)

Amber supports a switching function, but I'm not sure what it is. Here's what it says in section 21.7 (page 402) of https://ambermd.org/doc12/Amber22.pdf about fswitch:

When off, fswitch<=0 , uses a truncation cutoff. When on fswitch>0, sets a force switching region
where the force cutoff smoothly approaches 0 between the region of the fswitch value to the cut
value. Force values below the fswitch value follow the standard Lennard-Jones force. Default is
-1. This option is not supported for use with GB (i.e., only igb=0 and ntb>0), nor is it compatible
with the 12-6-4 Lennard-Jones model (lj1264=1). Due to performance regressions (about 20%)
with running with the force switching on, it is recommended that simulations run with fswitch off
unless using a force field that requires or recommends using the force switch.

So there's one option.

GROMACS has three options, one of which might line up with OpenMM's. I'm not sure what functional forms are used. Maybe one matches OpenMM exactly. From this extensive discussion one might match Amber's fswitch in principle, though it's unclear if it matches exactly.

LAMMPS has some options for shifts, but not switching. So that's not ideal.