Fix #882

[SimonBoothroyd]

Fixes #882

• Tag issue being addressed
• Add tests
• Update docstrings/documentation, if applicable
• Lint codebase
• Update changelog

[codecov]

Codecov Report

Merging #896 (d354bde) into master (e54052a) will decrease coverage by 9.97%.
The diff coverage is 100.00%.

[mattwthompson]

1. How much of an impact does patch have on physical properties? As far as I can tell, mainline OpenFF force fields using this change (and a new toolkit version) will use the switching function, as directed by the force field itself, whereas they currently do not because the toolkit currently does not respect this detail of the spec. I have a general sense of how much cutoff treatments matter, but I am not familiar with switching functions.
2. Any guesses why my code snippet is failing? (Is it also failing on yours/others' computers?) The unit tests are passing on my machine and I cycled the installation, so I hope it's just the ghost in my machine, but I'd like confirmation before digging into this more than a quick grep.

# Load a molecule into the OpenFF Molecule object
from openff.toolkit.topology import Molecule
from openff.toolkit.utils import get_data_file_path
sdf_file_path = get_data_file_path('molecules/ethanol.sdf')
molecule = Molecule.from_file(sdf_file_path)

# Create an OpenFF Topology object from the molecule
from openff.toolkit.topology import Topology
topology = Topology.from_molecules(molecule)

# Load the latest OpenFF force field definition
from openff.toolkit.typing.engines.smirnoff import ForceField
forcefield = ForceField('openff-1.2.0.offxml')

# Create an OpenMM system representing the molecule with SMIRNOFF-applied parameters
openmm_system = forcefield.create_openmm_system(topology)


from simtk import openmm, unit

assert forcefield['vdW'].switch_width == 1.0 * unit.angstrom
assert forcefield['vdW'].cutoff == 9.0 * unit.angstrom

for force in openmm_system.getForces():
    if isinstance(force, openmm.NonbondedForce):
        assert force.getUseSwitchingFunction()  # AssertionError
        assert force.getSwitchingDistance() == 1.0 * unit.angstrom

[SimonBoothroyd]

How much of an impact does patch have on physical properties? As far as I can tell, mainline OpenFF force fields using this change (and a new toolkit version) will use the switching function, as directed by the force field itself, whereas they currently do not because the toolkit currently does not respect this detail of the spec. I have a general sense of how much cutoff treatments matter, but I am not familiar with switching functions.

Oh wow I'd mis-read the force fields and had thought they had a switch width of 10A defined, not 1...

I realise also now that the spec is unclear about what a switch_width actually is - I had assumed that it was referring to the 'switch distance' as defined in the OMM docs here but now I'm not at all sure...

I think for this to be safely fixed a SMIRNOFF spec change is needed to include a flag which toggles switching on / off (default = off) and for the spec document to be updated to explicitly define how the switch_width is defined and how it fits in to the vdW energy expression...

[SimonBoothroyd]

as for point 2) the switch width is only set in cases where a cutoff is used, which currently require the input topology to have box vectors.

[mattwthompson]

It's a relief that I'm not alone in my confusion! I assumed that the switch width was a distance away from the cutoff that a treatment was applied, i.e. the SMIRNOFF defaults would have a switching function applied at 8 A with its cutoff at 9 A. Applying anything at 1 A doesn't make much intuitive sense to me - I think the terminology of switch distance used by OpenMM implies it's a value close to the cutoff distance, but switch width as used in the SMIRNOFF spec doesn't make clear to me if it's mean to be the distance at which a switch starts or the distance away from the cutoff - I'd assume the latter, but not confidently.

[j-wags]

Tagging @davidlmobley and @jchodera for feedback, since I'm also less clear on this:

I realise also now that the spec is unclear about what a switch_width actually is

For some reason I had imaged it being like a "skin distance", where "cutoff=9 switch_width=1" means that between 8 and 9, the switch is being applied. So in that case, we'd call NonbondedForce.setSwitchingDistance(self.cutoff - self.switch_width)

I think for this to be safely fixed a SMIRNOFF spec change is needed to include a flag which toggles switching on / off (default = off)

Can this instead be handled by a combination of "if switch_width==0, then turn switching on/off" and "never use a cutoff/switch if the system is nonperiodic"? Or are there cases where we'd want to use a switching function in a nonperiodic system? Genuinely curious, I don't know much of the background here.

and for the spec document to be updated to explicitly define how the switch_width is defined and how it fits in to the vdW energy expression.

Agreed. We should update the SMIRNOFF spec on this point.

[mattwthompson]

Or are there cases where we'd want to use a switching function in a nonperiodic system?

As far as I know, the coupling of periodicity to the cutoff method is a quirk of OpenMM's nonbonded forces and not necessarily fundamental to how cutoffs/corrections/etc. are applied in general. I believe in other engines (I'm familiar with GROMACS and LAMMPS but I would assume similar things for Amber, etc.) the periodicity and cutoff/correction/etc. are generally separated.

[j-wags]

It seems to me that having a larger switch_width than the cutoff should raise an error -- Silently disabling the switching it would likely be the more surprising behavior.

[j-wags]

I realise also now that the spec is unclear about what a switch_width actually is - I had assumed that it was referring to the 'switch distance' as defined in the OMM docs here but now I'm not at all sure...

I think for this to be safely fixed a SMIRNOFF spec change is needed to include a flag which toggles switching on / off (default = off) and for the spec document to be updated to explicitly define how the switch_width is defined and how it fits in to the vdW energy expression...

Given that this will affect the usage of our Sage fitting and release, and that the SMIRNOFF steering committee is unlikely to make a decision before then, let's deal with this as an implementation issue for now.

I think that the most natural interpretation of switch_width is what I described above (nonbondedforce.setSwitchingDistance(self.cutoff - self.switch_width)). but I could see how others would interpret it differently. I'm also not the user in this case, and have a strong track record of making user-unfriendly decisions, so maybe another definition is better for real-world use cases. So I'm happy to put my weight behind @SimonBoothroyd's interpretation (nonbondedforce.setSwitchingDistance(self.switch_width)). Any objection to that @davidlmobley or @mattwthompson ?

[jchodera]

switch_width was originally intended to reflect the distance span over which the switching function would be applied.
OpenMM uses a different "switching distance" to determine at what pair separation the switch function should turn on, but that does not allow the cutoff to be as independently controlled.

Note that we have see very little impact on the use or omission of switching functions for Lennard-Jones in OpenMM use, so I would anticipate the exact choice would have little impact. I recommend using a switching function if possible just to guarantee continuity of derivatives, but small switch widths could cause force artifacts near the cutoff. Switching also makes the long-range dispersion corrections a more complex to compute, and may not be supported by all downstream packages.

For example, it does appear that AMBER supports some manner of switching function for LJ, but this is specified not in the prmtop but in the mdin file, as the pair separation distance (in A) at which the switch turns on, but carries a 20% performance penalty:

gromacs appears to support both force-switch and potential-switch options, but (1) the potential-switch differs from CHARMM, and (2) they note it can create large artifacts at the boundaries, but that this is acceptable for LJ:

So it may be the case that simply including switch_width in the SMIRNOFF spec is insufficient because of the heterogeneity in switching functions supported by different packages. We would have to be explicit about the functional form the force field is parameterized with and work with package developers to ensure the same form is supported in those packages.

[j-wags]

Thanks for the detailed reply, @jchodera. I think we should to view this on two levels:

• Short term: The current behavior is "always turn switching off". Even without deciding on which behavior is correct, we can say that this behavior is definitely incorrect, since our force fields set a nonzero vdW switch_width. Nobody can be expected to match the toolkit's current behavior by reading the SMIRNOFF spec, and I'm thinking that this PR should just pick a plausible interpretation and run with it.
• Long term: Like you posted above, the SMIRNOFF spec is underspecified on how to correctly use the switch_width value in the force field, and it's likely that we'll need to update the specs for the vdW and Electrostatics sections. Our main-line FFs also specify a switch_width of 0 for Electrostatics, and 1 angstrom for vdW, which can not be supported by an OpenMM NonbondedForce. So that's a big thorny knot that we can begin to untangle with future SEPs, toolkit PRs, and FF releases.

John, would you mind if I copied your post above to the openFF standards issue tracker? That will be a good point to nucleate discussion and craft a Smirnoff Enhancement Proposal. In the meantime, I'd like to get this PR merged so that we're at least following some interpretation of the phrase "swtich width".

[j-wags]

@SimonBoothroyd I think this is the only substantive question remaining before I can make this merge-ready. The use of "width" makes me think that this is the more correct interpretation, since the switching function applies between this distance and the cutoff. But since you're more on the physics side, I'll defer this decision to you.

Suggested change

	force.setSwitchingDistance(self.switch_width)
	force.setSwitchingDistance(self.cutoff - self.switch_width)

[davidlmobley]

The width would usually refer to the distance over which it is tapered, e.g. an angstrom or two or some such (e.g. switched between 8 and 10 angstroms) if that helps answer this question.

[jchodera]

@j-wags is correct here: For OpenMM, you need to use force.setSwitchingDistance(self.cutoff - self.switch_width).

[mattwthompson]

Stipulating that

1. These changes update the toolkit to adhere to this particular detail of the SMIRNOFF spec (as it's understood and interpreted here)
2. These changes introduce a non-zero behavior change
3. Behavior changes, when justified, are to be grouped into version 0.11.0

I believe

1. These changes should accepted, targeting v0.11.0
2. The wording of the SMIRNOFF spec should be clarified (SMIRNOFF: Clarifying the meaning of switch_width standards#9), but no change is needed to it.

[davidlmobley]

Sounds good.

[mattwthompson]

@j-wags if the above plan sounds good to you, I'll update and re-open this to point to the refactor branch?

[j-wags]

Sorry for the delay -- yes, I think this one is good to go, and I agree with @mattwthompson's plans for next steps. Thanks for picking back up on this, everyone!

[mattwthompson]

These changes are now in the topology-biopolymer-refactor branch and therefore slated to be a part of the version 0.11.0 release.

I don't think we plan to include this change in the 0.10.x series of releases, based on my comment above, so I'm going to close this. Please re-open if I'm incorrect in my recollection!