-
Notifications
You must be signed in to change notification settings - Fork 92
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Fix #882 #896
Fix #882 #896
Conversation
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
- 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.
- 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
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 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 |
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. |
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. |
Tagging @davidlmobley and @jchodera for feedback, since I'm also less clear on this:
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
Can this instead be handled by a combination of "if
Agreed. We should update the SMIRNOFF spec on this point. |
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. |
force.setEwaldErrorTolerance(1.0e-4) | ||
|
||
force.setSwitchingDistance(self.switch_width) | ||
force.setUseSwitchingFunction(self.switch_width < self.cutoff) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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.
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 |
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 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 |
Thanks for the detailed reply, @jchodera. I think we should to view this on two levels:
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". |
force.setEwaldErrorTolerance(1.0e-4) | ||
|
||
force.setSwitchingDistance(self.switch_width) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@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.
force.setSwitchingDistance(self.switch_width) | |
force.setSwitchingDistance(self.cutoff - self.switch_width) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@j-wags is correct here: For OpenMM, you need to use force.setSwitchingDistance(self.cutoff - self.switch_width)
.
Stipulating that
I believe
|
Sounds good. |
@j-wags if the above plan sounds good to you, I'll update and re-open this to point to the refactor branch? |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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!
Co-authored-by: SimonBoothroyd <[email protected]>
* Port PR #896 to biopolymer-topology-refactor branch Co-authored-by: SimonBoothroyd <[email protected]> * Fix GBSA tests * Fix logic of calling setUseSwitchingFunction * Condense switching function logic, update tests * Update release history * Remove an annotation * Add comments about switching function with LJ-PME Co-authored-by: SimonBoothroyd <[email protected]>
These changes are now in the 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! |
Fixes #882