Default non-bonded cutoff should maybe be 0.9 nm

[mattwthompson]

OpenFF's mainline force fields use a 9 A cutoff for vdW interactions. Code paths involving openmmforcefields do not, by default, use this cutoff and instead fall back to the default of 1 nm:

>>> openmm.NonbondedForce().getCutoffDistance()
Quantity(value=1.0, unit=nanometer)

This can be set in general in OpenMM like

ForceField.createSystem(
    ...,
    nonbondedCutoff=0.9 * openmm.unit.nanometer,
    ...,
)

Unfortunately this is hardcoded in SystemGenerator.create_system which appears to be what this codebase calls. Something like .setCutoffDistance() after system creation should work in either case.

I don't know how much of an impact using 10 A has, but OpenFF's force fields thus far are trained with a 9 A cutoff in condensed-phase fitting and are shipped with that in the non-bonded settings.

Related https://github.com/openmm/openmmforcefields/pull/324

[IAlibay]

Thanks for raising this @mattwthompson !

I don't know how much of an impact using 10 A has, but OpenFF's force fields thus far are trained with a 9 A cutoff in condensed-phase fitting

This, similar to your other issue, is also something that would need testing / validation. I suspect cancellation of error minimizes the impact here but it will need further discussion.

Unfortunately this is hardcoded in SystemGenerator.create_system which appears to be what this codebase calls.

I will have to review this a bit later, but my understanding is that it is possible to set alternative cutoffs via system_generator / the OpenFE Protocol settings (for example see this test: https://github.com/OpenFreeEnergy/openfe/blob/main/openfe/tests/protocols/test_openmm_equil_rfe_protocols.py#L560).

P.S. The link you point to is for DummySystemGenerator, I'm not sure if this is intentional?

[mattwthompson]

Good catch, I'm confused by the different classes and method names. It might get passed through here when called here, not sure if it's set later?

[IAlibay]

We're passing the cutoff on object creation here: https://github.com/OpenFreeEnergy/openfe/blob/main/openfe/protocols/openmm_utils/system_creation.py#L97

Which then gets used within SystemGenerator here: https://github.com/openmm/openmmforcefields/blob/0.12.0/openmmforcefields/generators/system_generators.py#L333

That being said, one thing to note is that we should add some tests for the AFE and MD Protocols that custom cutoffs are used/retained (I'm reasonably sure that they are, but sure != certain).

[mattwthompson]

Ah okay I missed those are passed to that constructor. I'm used to using createSystem

[mattwthompson]

FTR

• This is, I think, accurately captured in the SMIRNOFF example upstream
• Because openmmforcefields does not itself implement ForceField.createSystem itself, users must remember those settings if using that code path
• SystemGenerator itself is an openmmforcefields class, so that package has more control over non-bonded settings there
• The hard-coding in DummySystemGenerator.create_system is something that could be changed upstream (cc: @epretti), maybe requiring logic to remember which force field a generator came from since GAFF need not implement SMIRNOFF non-bonded settings.

I'm not sure if OpenFE is doing something wrong here per se, but it's a sharp edge to be aware of in general

[epretti]

I'm finally getting around to taking a look at this... sorry for the delay.

It looks like the hard-coding in DummySystemGenerator is intentional, or at least not a problem. See the docstring:

    Dummy force field that can add basic parameters to any system for testing purposes.

    * All particles are assigned carbon mass
    * All particles interact with a repulsive potential
    * All bonds have equilibrium length 1 A
    * All angles have equilibrium angle dependent on number of substituents of central atom
        * 2, 3 bonds: 120 degrees
        * 4 bonds: 109.8 degrees
        * 5 or more bonds: 90 degrees
    * Torsions are added with periodicity 3, but no barrier height

In other words, it seems like this shouldn't be used for any real system generation, just basic testing, and I didn't see anywhere in OpenFE that uses this, so I don't think OpenFE is doing anything incorrectly.

I believe that the way that the openmmforcefields API is set up now is that whoever is using SystemGenerator should be passing the cutoff, and any other system creation options, into [non]periodic_forcefield_kwargs. I agree that this is something users have to be aware of and perhaps openmmforcefields could handle it better. However, it's not clear to me what an automatic approach would look like. Users could define their own template generators, pass in arbitrary FFXMLs (in which cutoff information is not stored)—what if there's a conflict, or a user wants control over the cutoff to optimize PME performance, for instance? If you think that some form of automated cutoff selection could be a useful feature, though, please feel free to open an issue on openmmforcefields to discuss more.

[mattwthompson]

I'm like 80% sure that this is not actually an issue for OpenFE as scoped here, or something that would only benefit from testing/documentation improvements - please

Outside of the scope of OpenFE, SMIRNOFF and OpenMM disagree on what constitutes a force field so there will always be some pain at that interface

[IAlibay]

@mattwthompson I'm going to re-open this issue, not for the "nonbonded cutoff isn't being passed" issue, but rather "the default nonbonded cutoff value is 1.0 when it should be 0.9".

As part of the "fast settings", we are going to start encouraging folks to move towards 0.9 nm, and hopefully we'll do a full switch in a 2.0. Until we do, it's best to keep this issue open as a reminder that we have to change that default.