vdW energy mismatch between Amber and GROMACS/OpenMM

[mattwthompson]

In [23]: ff14sb = ForceField("ff14sb_off_impropers_0.0.3.offxml")

In [24]: out = Interchange.from_smirnoff(
    ...:     force_field=ff14sb,
    ...:     topology=[Molecule.from_polymer_pdb(get_data_file_path('proteins/MainChain_HIE.pdb'))],
    ...: )

In [25]: out.box = [10, 10, 10]

In [26]: get_summary_data(out)
Out[26]:
             Bond      Angle    Torsion       vdW  Electrostatics
OpenMM   2.041475  80.760040  54.966743  5.605402     -211.167846
Amber    2.041374  80.759986  54.966882  6.151735     -211.158530
GROMACS  2.041476  80.760101  54.673031  5.606108     -211.196062

[mattwthompson]

Bookmarking for later: https://gromacs.bioexcel.eu/t/mdp-parameters-for-non-bonded-interaction-using-the-amber-force-field/6458/5

[mattwthompson]

Okay, the problem is that fswitch does not apply a switching function. It applies a shift to the entire potential. The interaction energy at the cutoff is simply subtracted out of the potential for the entire range, disregarding where a shifting function is meant to be applied.

If I turn off the switching function (vdWHandler.switch_width *= 0.0), nothing in particular happens at 8 A and there's a hard, discontinuous cut-off at 9 A:

(The pseudo-periodic noise is probably just due to float rounding; Amber reports only to something like 0.001 kcal/mol by default.)

If use OpenFF's current default (vdWHandler.switch_width = 8 * unit.angstrom) you can see that the energies reported by OpenMM and GROMACS code paths are unaffected below 8 A and follow a continuous transition to smoothly decay to zero at 9 A. What's reported through our Amber code paths (including passing fswitch=8 through to sander) is a complete shift of the potential over the entire range:

I'm not sure where the value 8 ends up; nothing happens at 8 A so maybe it'll behave the same if I use 7 or 6.

From the Amber22 manual (emphasis mine):

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.

Something is clearly going wrong here. I can't find any mention of an option that adds a potential shift (like vdw-modifier = Potential-shift in GROMACS or pair_modify shift yes in LAMMPS) to sander.

For good measure, here's the second plot zoomed out a bit - clearly what's happening is applied the same everywhere below 8 A:

import matplotlib.pyplot as plt
import numpy as np
from openff.interchange.drivers import *
from openff.toolkit import *

from openff.interchange._tests import get_test_file_path

ion_ff = ForceField(get_test_file_path("ions.offxml"))
# ion_ff['vdW'].switch_width = Quantity(0.0, "angstrom")
ions = Topology.from_molecules([
    Molecule.from_smiles("[#3+1]"),
    Molecule.from_smiles("[#17-1]"),
])
ions.box_vectors = Quantity([4, 4, 4] , "nanometer")

out = ion_ff.create_interchange(ions)

fig, (ax1, ax2) = plt.subplots(2)

r = np.linspace(0.3, 1.0, 101)

omm, gmx, amb = list(), list(), list()

for d in r:
    out.positions = Quantity(
        np.array(
            [
                [0.0, 0.0, 0.0],
                [d, 0.0, 0.0],
            ]
        ),
        "nanometer",
    )

    omm.append(get_openmm_energies(out, combine_nonbonded_forces=False).energies['vdW'].m)
    gmx.append(get_gromacs_energies(out).energies['vdW'].m)
    amb.append(get_amber_energies(out).energies['vdW'].m_as(unit.kilojoule / unit.mol))

ax1.plot(r, omm, 'k-', label='OpenMM')
ax1.plot(r, gmx, 'c.', label='GROMACS')
ax1.plot(r, amb, 'm.-', label='Amber')
ax2.plot(r, np.asarray(omm) - np.asarray(gmx), 'c.', label='OpenMM - GROMACS')
ax2.plot(r, np.asarray(omm) - np.asarray(amb), 'm.', label='OpenMM - Amber')

ax1.set_ylabel(r'$U(r)$')
ax2.set_ylabel(r'$U(r) - U_{OpenMM}(r)$')

ax1.set_xlabel("r, nm")
ax1.legend(loc=0)
ax2.legend(loc=0)
fig.show()

Here are the files used under the hood, with and without fswitch meant to be applied:

single-point energy
&cntrl
imin=1,
maxcyc=0,
ntb=1,
fswitch=8.0,
ntc=1,
ntf=1,
cut=9.0,
/
&ewald
order=4
skinnb=1.0
//

single-point energy
&cntrl
imin=1,
maxcyc=0,
ntb=1,
fswitch=-1.0,
ntc=1,
ntf=1,
cut=9.0,
/
&ewald
order=4
skinnb=1.0
//

[mrshirts]

One bigger picture thing that might be applicable logic to use here - I decided with InterMol that "success" was that all the parameters were translated correctly; if the functional forms were not the same, that wasn't the converter's problem. i.e. if I could get the energies to match with 9 A abrupt cutoffs in both engines, then that was matching. Lack of similar implementation between programs was the program's problem, not mine.

[mrshirts]

fswitch is I think short for force-switching, which is not the same multiplying the potential by a scaling function. I don't think there is a way to multiply the potential by a switching function, yes.

[mattwthompson]

I agree that getting parameters to being written correctly to the data files is a good target; the problem is that OpenFF's force fields use a switching function and there's no better place to try and get that translated that here. People will get different results if the switching function is turned into a hard cutoff. How much? Probably less than other error in most cases, I don't know.

Unfortunately it seems the best we can do in this case is document this gap in functionality so that users are aware of what non-bonded settings are included in the force field that can't precisely be used in their workflow.

[mrshirts]

Unfortunately it seems the best we can do in this case is document this gap in functionality so that users are aware of what non-bonded settings are included in the force field that can't precisely be used in their workflow.

Exactly. I would not change it to a hard cutoff, that was just the setting that was implementable in all version of the codes InterMol used, and therefore could be used to show that in THAT mode, the energies were the same. I think using fswitch for AMBER outputs and WARNING the user the change that was made is 100% a great way to go. We can't make simulation engines more cross-compatible.

[mattwthompson]

Agree on warning users - but my intuition is to default to not using fswitch at all for these cases. I think shifting the force is probably better for practical use but it's more disruptive (shifting the potential over the entire range introduces more numerical error, even if the physics are usually fine, compared to having an incorrect cut-off and some error accumulating between 8 and 9 A) compared to what's "supposed to be used." I have this implemented in #914, including the same warning for LAMMPS since we weren't able to figure out that case.

The scope here is internal energy comparisons and an API point that only ever writes out a script to set up a 0-step "simulation." In both cases a warning will be emitted. I expect power users will have their own infrastructure for production runs and aren't as worried about the tail correction, but I do see it as our responsibility to communicate out the settings that are encoded in force field.

[mattwthompson]

In #914 I have added a warning - until Amber implements this particular type of tail correction, there isn't anything we can do