OFF-EP 0005b

[jchodera]

An alternative proposal to #30 to fix #29.

Other modifications that may be possible:

• We could use solvent_dielectric to specify both reaction field dielectric constants and the Ewald boundary conditions dielectric constant (currently specified as part of the verbose Ewald3D-ConductingBoundary keyword)
• We could migrate the specification of CODATA 2018 physical constants to a higher-level attribute
• We could streamline or specify which physical constants are available in the potential expressions (here, pi and epsilon0 are used)
• We could specify some common reaction field models (such as the OpenMM reaction field model used in SireMol) in this OFF-EP rather than relying on users to provide the complete functional form

[jchodera]

I haven't yet updated the smirnoff.md from #30, but please take a look at the proposed alternative OFF-EP 0005.

[davidlmobley]

I think I'm fine with this as well; I'd like to hear from @mattwthompson and @j-wags but it still seems relatively compact.

[mattwthompson]

I don't have a seat on the committee and but since I was asked -

#30 solves, in my view, the problems we have identified as short-term blockers for version 0.11.0 of the toolkit (and a 0.2.0 release of Interchange, which will be paired with it) as I am trying to prevent #29 from creeping into Interchange. I intended for it to be unambiguously discrete in scope to the changes I need, both because I have not done research cycles on other topics and because I hoped it would be easier to review a more narrow proposal; the assumption is that it is easier to gain consensus on one move part than several. (This has been my hope with the OFF-EPs I have written, that targeted, narrow proposals would be easier to review, much like PRs to source code - even at the cost of having a larger number of proposals creating extra chatter.)

I think specifying physical constants would be an excellent improvement to SMIRNOFF, but I'm not sure it's in scope here and think it would be better split out into a separate proposal. Two questions that immediately come to mind are

• How should implementations handle exports to engines that use different versions of CODATA (or a different system altogether)? Forbidding export to to i.e. GROMACS 2019 (which uses CODATA 2010) would be onerous, but allowing implementations to skip this check altogether implies it's not necessary - or should be communicated as a recommendation, not a requirement.
• Assuming one doesn't wish SMIRNOFF to be on CODATA 2018 forever, what is the process for applying updates?

PME (and other methods) are permissible approximations to Ewald as long as they are controlled.

These needs clarification, presumably with coming changes to smirnoff.md. What does "controlled" mean here? What approximations are permissible, i.e. is it fine that Amber and OpenMM use different B-splines of different order?

Should the default dielectric be 78.3 (water) or 78.5 (consistency with Amber)? Does this difference matter? Separately from this proposal, it would be worth seeing if OpenFF's fitting/benchmarking pipelines are consistent - both in terms of this constant and any other changes to the electrostatics definitions.

I have no idea if the reaction field expression given in the example is supported by engines other than OpenMM; since it's not been an operational priority, I haven't dug into it and it's not a technique I'm familiar with. So I'm not able to given feedback on it except to repeat the same concern as above - if it's a requirement, it might be ignored when interfacing with slightly different implementations.

I also don't know how different engines compute electrostatics of intramolecular interactions. To be honest, it's only recently come to my attention that this is a detail that engines have opinions about.

I just need the specification revised to support the most common use case of our force fields: PME (with details possibly left unspecific) for periodic systems and something-close-to-exact-Coulomb for non-periodic systems. This proposal seems like it only adds content, if I'm reading things accurately. Assuming that updates to smirnoff.md in a final version of this proposal don't modify this, that need is met. I'm unable to provide useful feedback on some other points, though, so I will defer to the committee's judgement unless asked again.

[j-wags]

Thanks @jchodera for adding your ideas to this proposal. I made several comments but none of them are necessary to resolve for me to approve this EP. I think the proposed changes are great and would approve this EP once the changes are propagated through to the smirnoff.md file.

You and @mattwthompson have done most of the writing on this, so I'd be happy to edit the smirnoff.md file with the proposed changes, and the "Backward compatibility" section of this EP with a more thorough recommendation for up-conversion (I'll have to write an up-converter anyway, so it couldn't hurt to specify exactly what it will do). Please tell me if you'd like me to make these edits and I can do so on Monday!

[j-wags]

(not blocking) I would approve an EP modifying the text of the 0.3 SMIRNOFF specification to say that 2018 CODATA values are assumed. I would approve that change in this EP or a different one, but given the asynchronous nature of our org I'd prefer that go in another EP.

I'd also approve an EP to add a top-level attribute to a SMIRNOFF FF specifying which CODATA should be used. But that should definitely be a different EP.

I share Matt's concerns that ensuring use of a specific CODATA will be very difficult when crossing the the boundaries of software ecosystems (like, when we export to another package, what if they use different physical constants? What if that other package updates its constants in a new release - How would we know which version we're exporting to?). But the "different physical constants" problem already exists. This step would change it from a "the answer is not defined" issue into a "the implementation doesn't follow the spec" issue, and that's an improvement to me. I don't know when I'll have time to do a deep audit of the physical constants in our current implementation, but a bug in the tracker is better than undefined behavior.

[jchodera]

I share Matt's concerns that ensuring use of a specific CODATA will be very difficult when crossing the the boundaries of software ecosystems (like, when we export to another package, what if they use different physical constants?

As I note below, if the constants defined in the force field and simulation package are explicitly specified/known, the parameters can be correctly converted to reproduce the intended energies.

What if that other package updates its constants in a new release - How would we know which version we're exporting to?).

This question answers itself: Export requires the specific release of the package to be specified in order to correctly export the parameters in a manner that ensures the energies are correctly computed.

[j-wags]

(not blocking) Thankfully, I don't think this is a concern. Since this change will be accompanied by an Electrostatics section version change, implementations will know which keywords are expected from each version. My plan for the Toolkit would be to just automatically up-convert any Electrostatics section that is read from 0.3 to 0.4 (though if there are strong objections I could try to implement round trip support for both versions, but that would be expensive, brittle, and complex).

[j-wags]

(not blocking) To other reviewers - I have no strong opinion about this default value, but you should be aware that 78.5 matches the default solvent_dielectric value in the GBSA tag.

[SimonBoothroyd]

It should likely be none unless a reaction-field / custom expression that makes use of it is used?

[j-wags]

PME (and other methods) are permissible approximations to Ewald as long as they are controlled.

These needs clarification, presumably with coming changes to smirnoff.md. What does "controlled" mean here? What approximations are permissible, i.e. is it fine that Amber and OpenMM use different B-splines of different order?

(not blocking) This would be good to answer. I don't see it as a hard blocker for my approval, and I don't know enough about PME implementations to write this in smirnoff.md myself. But I'd basically trust any reasonable scheme agreed upon by more knowledgeable committee members.

[jchodera]

@mattwthompson @j-wags : Thanks so much for your detailed feedback!

Scope: My assumed scope differed from @mattwthompson's #29 in that I considered the scope of this OFF-EP to be "eliminate ambiguities in the definition of the Electrostatics tag", rather than "eliminate one ambiguity in the definition but leave a boatload of other ambiguities". I appreciate that there is tension between the desire to make narrowly-scoped, easily reviewable changes and the amount of additional work required to review, accept, and then implement each change, especially if this requires backward- (and possibly forward-)compatibility. I'm totally OK with implementing a series of narrowly-scoped changes if the software scientists are fine with implementing all of them, but it seems like this has the potential to create a great deal more work.

I think specifying physical constants would be an excellent improvement to SMIRNOFF, but I'm not sure it's in scope here and think it would be better split out into a separate proposal.

With my assumed scope of eliminating ambiguities, this was an ambiguity in the definition of how the potential was calculated that we need to explicitly address in this OFF-EP.

How should implementations handle exports to engines that use different versions of CODATA (or a different system altogether)? Forbidding export to to i.e. GROMACS 2019 (which uses CODATA 2010) would be onerous, but allowing implementations to skip this check altogether implies it's not necessary - or should be communicated as a recommendation, not a requirement.

If we know which constants a package uses internally to compute the potential energy from our parameters, it is trivial to adjust the parameters on export to ensure the energies match their intended values. This is only possible if we specify what physical constants we intend the spec to use. If we do not, we will just get different energies in each package.

Assuming one doesn't wish SMIRNOFF to be on CODATA 2018 forever, what is the process for applying updates?

As noted above, a future OFF-EP could explicitly specify the physical constants used at a higher level---perhaps at top level---and remove the specification of which physical constants are used from the Electrostatics tag. We could then update the physical constants as desired. This change was outside of my assumed scope, however---which was confined to just the Electrostatics tag.

These needs clarification, presumably with coming changes to smirnoff.md. What does "controlled" mean here? What approximations are permissible, i.e. is it fine that Amber and OpenMM use different B-splines of different order?

We need to agree on a criteria for this, but my assumption here is that we will first need to see some data before we can construct a criteria. One possible approach, for example, would be to curate a set of test systems for which we could provide Interchange parameters and reference energies for a true Ewald implementation, and benchmark the exported energies in other packages with a variety of settings (e.g. different pme-tolerance values in OpenMM with PME). We could then specify "your settings should reproduce these energies to within X kcal/mol maximum deviation on the provided systems and snapshots". This is but one idea, but it would provide a concrete way of determining whether the force field is correctly implemented in a package with particular settings.

Should the default dielectric be 78.3 (water) or 78.5 (consistency with Amber)? Does this difference matter? Separately from this proposal, it would be worth seeing if OpenFF's fitting/benchmarking pipelines are consistent - both in terms of this constant and any other changes to the electrostatics definitions.

We can choose not to define a default dielectric in this spec.

I have no idea if the reaction field expression given in the example is supported by engines other than OpenMM; since it's not been an operational priority, I haven't dug into it and it's not a technique I'm familiar with. So I'm not able to given feedback on it except to repeat the same concern as above - if it's a requirement, it might be ignored when interfacing with slightly different implementations.

I did not define any keywords for reaction field electrostatics because I haven't had time to review which functional forms are implemented in different packages. It will require some effort to make sure implementations are standardized or can be implemented or approximated in multiple packages. For now, the safest thing to do is either exclude them (which is inconvenient, since we know we need to support Cresset with an optimized reaction field form) or explicitly specify the functional form (which likely for Cresset with their OpenMM implementation). Some refinement of this part of the proposal may be necessary, however, if we want to easily be able to render this to NonbondedForce instead of CustomNonbondedForce in Interchange.

I also don't know how different engines compute electrostatics of intramolecular interactions. To be honest, it's only recently come to my attention that this is a detail that engines have opinions about.

These details matter for the energies they produce, so we cannot ignore them.

You and @mattwthompson have done most of the writing on this, so I'd be happy to edit the smirnoff.md file with the proposed changes, and the "Backward compatibility" section of this EP with a more thorough recommendation for up-conversion (I'll have to write an up-converter anyway, so it couldn't hurt to specify exactly what it will do). Please tell me if you'd like me to make these edits and I can do so on Monday!

@j-wags : Please feel free to make further changes to the branch directly! My goal here was to provide the basis of an alternative OFF-EP that more thoroughly addressed all identified ambiguities in the Electrostatics tag, but I am not committed to any of the specific solutions proposed here.

I'll respond to the other comments inline.

[mattwthompson]

The remaining loose ends (PME details, reaction field implementation, intramolecular interactions) are each improved or not made worse by this proposal, so I don't see any reason to reject these proposals (a. I think we're in a good position to tackle them in the future.

[j-wags]

I'm updating the SMIRNOFF spec text and realized that I don't understand this - The existing Electrostatics tag already has a cutoff attribute, which defaults to 9.0 * unit.angstrom. Is the proposal here to split the existing cutoff attribute into cutoff_periodic and cutoff_nonperiodic?

[mattwthompson]

I think (especially if I wrote this, I don't know) that line 88 should refer to "cutoff_potential if applicable ..."and not allowed/ignored for Ewald/Ewald-X methods that do not have a concept of a cutoff.

Another cup of coffee later, I guess we do have to consider the possibility of encoding separate cutoffs for periodic and non-periodic methods? I think reaction field has the concept of a cutoff? If the periodic method is reaction field and the non-periodic method is ... I don't know, some sort of hard cutoff, I guess they could differ?

OpenFF's common use case of PME + no-cutoff would result in both being ignored. I'm fine kicking the can down the road on this, but I think the most thorough solution suitable for the long term would be to split cutoff into two new attributes.

[j-wags]

I've assumed that the answer to my last question is "yes, we need to split cutoff into periodic_cutoff and nonperiodic_cutoff", then realized that we need to do the same with switch_width, and gone ahead with:

• extending the Reverse Compatibility section of the EP with suggestions for upconverting to the new spec
• propagated the proposed changes into the SMIRNOFF spec text

I approve merging this EP in its current form.

[jchodera]

Thanks for the changes, @j-wags!

Do we expect that (1) we will want to use a cutoff for non-periodic systems (which are primarily intended to match QM calculations) at some point in the future, and (2) we will want to at some point enable the cutoffs to be distinct in these cases? If so, this seems like a good future-proof change. I just couldn't envision a use case for either.

[j-wags]

Do we expect that (1) we will want to use a cutoff for non-periodic systems (which are primarily intended to match QM calculations) at some point in the future, and (2) we will want to at some point enable the cutoffs to be distinct in these cases?

No opinion - my experience with cutoff potentials is limited. @SimonBoothroyd or @davidlmobley, do you have an opinion?

[mattwthompson]

Do we expect that (1) we will want to use a cutoff for non-periodic systems

The Shirts, 2017 way of testing for energy equivalence (which I would like to be able to do, #6) uses a hard cut-off for electrostatics. But I don't know if this is something that must be supported for both periodic and non-periodic systems. I can't think of a use case outside of that (which would be strictly testing-only and never touch production).

If there is any friction associated with officially supporting this feature I'd just modify a system after preparation to do this, which is only a tiny detour away from official SMIRNOFF.

(2) we will want to at some point enable the cutoffs to be distinct in these cases?

I can't think of any cases; a cutoff being used in either periodic or non-periodic treatments should be rare and the possibility of shipping a force field using both should be vanishingly rare ... but I guess possible?

[davidlmobley]

Re the question directed at me: I personally don't expect a hard cutoff for nonperiodic systems. And I don't see a use case for separate cutoffs in the two cases.

I also read everything else here and am fine with where this has ended up, so I approve this EP as written in its current form. @SimonBoothroyd ?

[SimonBoothroyd]

Thanks @jchodera these changes would indeed be a very welcome clarification / improvement to the spec.

I also can't think of a reason to use a hard cut-off with a plain Coulomb potential, and would suggest this be disallowed by the spec. If a user really wants this behaviour they can just provide the Coulomb potential as a custom expression and set the cut-off.

I think the two documents need another pass through to ensure consistency as the EP and SMIRNOFF spec seem to be mismatched in places.

[j-wags]

Thanks for the close read, @SimonBoothroyd and everyone else! I'll take another swing this today - basically

• reverting the {non,}periodic_{cutoff,switch_width} fields to what they were before, and clarify that they only apply to periodic potentials and custom functional forms
• Fixing the inconsistencies that Simon pointed out

Some of the defaults that Simon suggested changing to none will require a bit more thought. I remember there's a mess of complexity when a default value has to change, but I can't verbalize it well now. I'll think this through and write it up with my edits.

[j-wags]

Ok, I think this is in good shape now. Thanks again @SimonBoothroyd for catching the inconsistencies.

The notable changes now (beyond just fixing inconsistencies) are:

• There's now just one cutoff and switch_width (I've un-split it from having separate periodic_ and nonperiodic_ values)
• The upconversion table has the full formula for reaction-field
• switch_width, cutoff, and solvent_dielectric now default to "none" in the 0.4 Electrostatics section spec
• The upconversion section recommends that, when upconverting 0.3 Electrostatics tags containing method=PME and non-none cutoff or switch_width, that those values should be upconverted to none (This will prevent massive section-compatibility headaches in workflows like "load sage and then append this other FF", trust me!)
• Because cutoff/switch_width/solvent_dielectric would be required to be non-none if ANY of the three possibly-custom potentials included them, I think it'd be futile/brittle to check for all three terms in all three potentials, so I've softened the statement

Only cutoff="none" and switch_width="none" is allowed for Ewald methods---only finite-ranged methods should set these to a non-none value.

to

It is possible to define an Electrostatics section where no potential uses cutoff, switch_width, or solvent_dielectric. In these cases it is strongly recommended that these values be set to none to avoid ambiguity.

[SimonBoothroyd]

Just a few last inconsistencies in the text but otherwise looks good to me.

[SimonBoothroyd]

Suggested change

	* `Coulomb` denotes that the standard Coulomb potential should be used with specified cutoff `cutoff` and optionally switch width `switch_width`
	* `Coulomb` denotes that the standard Coulomb potential should be used with no cutoff or reaction-field attenuation