Port plugins from Interchange

[mattwthompson]

Description

A different version of openforcefield/openff-interchange#45, copying over some existing code.

Todos

Notable points that this PR has either accomplished or will accomplish.

• Port damped Buckingham potential
• Port double exponential potential
• Create a base class for plugin collections
• Support on virtual sites (at least charges)
• Investigate differences in scaled energies
• Remove double 1-4 scaling (resolves Store "slot" identifiers as str openff-interchange#51)
• Update tests

Questions

• Existing force fields with an empty vdW handler might cause issues with Interchange's logic around guessing how dispersion interactions are descript - deal with this here or in specific force field repos?

Status

• Ready to go

[jthorton]

hey @mattwthompson, thanks for writing this up made it a lot easier to track down, the clue is that the first energy you report is very close to the expected energy at r_min so it seemed like the distances between the waters were wrong, this is now fixed. I also had to make the box vectors and cutoff larger to get the test to pass however, I tried setting the interchange box vectors to None as I expected this to set the vdW cutoff to NoCutoff but it didn't seem to work. Can we have no cutoff for this test? if not we can always just set the box vectors to be larger than the system.

[mattwthompson]

Thanks for fixing the tests! The energy differences reporting in CI are now small, though not quite to the desired tolerance.

Can we have no cutoff for this test?

Maybe - the issue that concerns me is how vdW and custom vdW handlers define method="cutoff" and producing something with OpenMM's NoCutoff is in conflict with that. If I understand/recall correctly, when the non-bonded interactions are split out into multiple forces, OpenMM throws an exception when trying to mix PME electrostatics with NoCutoff set on the CustomNonbondedForce. I try to catch this specifically since it's easy to run into.

if Interchange is provided a force field that says its vdW(-like) interactions should be cutoff, I don't want to give the user back a force with no cutoff. In spec, there's nothing that covers this case, but out of spec I imagine we can fudge it by allowing a method="nocutoff" via modifying this line. That would provide Interchange enough information to know the user really does want NoCutoff. I'd punt on whether or not a force field fit to a 9 A cutoff can be used with no cutoff, since I don't know if that's a trivial or significant difference.

if not we can always just set the box vectors to be larger than the system.

I'd like this to work but it's already set fairly large (20 nm) and continuing to bump it up might cause tests to run much slower (at least on my machine it does). Messing around with the platform and maybe loosening the tolerance could help here, but it could be fragile and might not be worth the trouble of deviating from the hand-calculated reference energies.

When I bump it to 40 nm and switch to the Reference platform, the double exponential energy test passes locally but the damped Buckingham fails (if barely):

E           assert 329.30501078508985 == 329.30542 ± 1.0e-05
E             comparison failed
E             Obtained: 329.30501078508985
E             Expected: 329.30542 ± 1.0e-05

of course this is probably going to run differently on CI runners, different developers' machines, etc.

Do you have a preference between these options?

[jthorton]

if Interchange is provided a force field that says its vdW(-like) interactions should be cutoff, I don't want to give the user back a force with no cutoff.

What about a second method like create_gas_system or having an argument in the function to declare the user wants a gas phase system with no cutoff so it is very clear? Then it's up to other programs like forcebalance to decide when it wants a gas system vs a system with cutoffs and say PME electrostatics.

I think for this test we can manually change the system to have no cutoff to match the values I calculated by hand with no cutoffs.

Seems like the last issue is the Buckingham water test simulation which has v-sites should we expect that to work yet?

[mattwthompson]

What about a second method like create_gas_system

Do you imagine this as an alternative to Interchange.to_openmm or somewhere else? It's fine if it's OpenMM-specific since other engines either don't support gas-phase simulations or aren't as flexible with non-bonded functional forms.

I think for this test we can manually change the system to have no cutoff to match the values I calculated by hand with no cutoffs.

Sounds good to me.

Seems like the last issue is the Buckingham water test simulation which has v-sites should we expect that to work yet?

I think so - the virtual sites are only on the water and only include charges, right? That should work since that's the same thing that the double exponential model supports, I probably just forgot to copy something over. I'll have a look.

[codecov-commenter]

Codecov Report

Patch coverage: 88.96% and project coverage change: -3.86 ⚠️

Comparison is base (051ab7c) 79.06% compared to head (98685ab) 75.21%.

Additional details and impacted files

@@            Coverage Diff             @@
##             main      openforcefield/openff-interchange#53      +/-   ##
==========================================
- Coverage   79.06%   75.21%   -3.86%     
==========================================
  Files           5        4       -1     
  Lines         301      238      -63     
==========================================
- Hits          238      179      -59     
+ Misses         63       59       -4     

Flag	Coverage Δ
unittests	75.21% <88.96%> (-3.86%)	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

Impacted Files	Coverage Δ
smirnoff_plugins/utilities/openmm.py	52.88% <69.76%> (+5.94%)	⬆️
smirnoff_plugins/collections/nonbonded.py	95.60% <95.60%> (ø)
smirnoff_plugins/handlers/nonbonded.py	85.00% <100.00%> (-8.89%)	⬇️

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

☔ View full report in Codecov by Sentry.
📢 Do you have feedback about the report comment? Let us know in this issue.

[mattwthompson]

Some egg on my face; even though Interchange.positions was properly set to the hand-crafted positions matrix, to_openmm_positions had a bug in it that didn't read those positions if the molecules had conformers. The way this was constructed, the copies of the water molecule each had the same conformers, so all molecules were placed on the same positions. I have tests passing here with a quick hack (removing the molecule's conformers after making the Topology) that's functional but not aesthetic. The next tag or release of Interchange will have the relevant fix in it.

However, now that tests are passing, I think this is finally ready for a review!

[jthorton]

If possible I would like to remove some of the duplicated code in the collections to simplify making new plugins, otherwise, this looks great!

[jthorton]

What is this method intended to do? Why would the creation of a v-site be different?

[mattwthompson]

Well, this one turned out a little weirder than I expected. The short explanation is that nothing in this repo needs to handle processing virtual sites as a special case, so I've just removed that method altogether.

Earlier I put this in because I wanted to get the non-bonded functionals working before making sure the virtual sites logic was also working, so throwing an exception was just a stop-gap. The way Interchange stores virtual site information is a little complex, it does it in three steps: It turns out Interchange doesn't even use these methods anymore! No wonder this wasn't causing any problems. Previously I was thinking of a design in which each non-bonded handler made its own virtual sites, but now I have it set up so that a single class manages the virtual sites' geometry, charges, and vdW parameters. That complexity is tucked away here, and I'm pretty sure it works just the same as if the 4-site water model used 12-6 Lennard Jones; it just creates a non-interacting potential:

In [1]: from openff.toolkit import ForceField, Molecule

In [2]: from rich.pretty import pprint

In [3]: interchange = ForceField(
   ...:     "buckingham-force-field.offxml",
   ...:     load_plugins=True,
   ...: ).create_interchange(
   ...:     Molecule.from_smiles("O").to_topology(),
   ...: )

In [4]: pprint([*interchange["DampedBuckingham68"].potentials.values()])
[
│   Potential(
│   │   parameters={
│   │   │   'a': <Quantity(1600000.0, 'kilojoule_per_mole')>,
│   │   │   'b': <Quantity(42.0, '1 / nanometer')>,
│   │   │   'c6': <Quantity(0.003, 'kilojoule_per_mole * nanometer ** 6')>,
│   │   │   'c8': <Quantity(3e-05, 'kilojoule_per_mole * nanometer ** 8')>
│   │   },
│   │   map_key=None
│   ),
│   Potential(
│   │   parameters={
│   │   │   'a': <Quantity(0.0, 'kilojoule_per_mole')>,
│   │   │   'b': <Quantity(0.0, '1 / nanometer')>,
│   │   │   'c6': <Quantity(0.0, 'kilojoule_per_mole * nanometer ** 6')>,
│   │   │   'c8': <Quantity(0.0, 'kilojoule_per_mole * nanometer ** 8')>
│   │   },
│   │   map_key=None
│   ),
│   Potential(
│   │   parameters={'sigma': <Quantity(0.0, 'angstrom')>, 'epsilon': <Quantity(0.0, 'kilocalorie_per_mole')>},
│   │   map_key=None
│   )
]

These three objects are the Buckingham interactions for oxygen, hydrogen, and the virtual site, respectively.

If you ever want to have custom dispersion interactions on the virtual sites this'll get tricker. I'm honestly not sure if getting this to look correct would be a small fix or large undertaking, and if you don't need it now then I'd prefer leaving this as-is. These changes would probably go into Interchange, not here.

[jthorton]

LGTM, great job!

[mattwthompson]

Thanks!