-
Notifications
You must be signed in to change notification settings - Fork 3
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
SMIRNOFF: PME electrostatics cannot be applied to non-periodic systems #29
Comments
I agree that the spec should be updated (to at least record the behavior that users have been getting using our FFs), and would be happy with either of these courses of action. If needed I can generate a preference, but I'd approve an EP that adds either change. One alternative that came to mind is that we could expand the spec to have |
@jchodera @davidlmobley @SimonBoothroyd Resolving this question is high-priority. Wearing my "Toolkit developer" hat: This issue is effectively blocking to biopolymer support. The chain of logic is:
Now, putting on my "SMIRNOFF committee member" hat: Per my last message, I think that @mattwthompson's proposed solutions would satisfy everyone. Before he takes the time to draft a formal EP to change the spec, could the other committee members weigh in? |
To take the last first, I think it's important that doing what I had suggested would break a lot of existing code, as you point out, @j-wags , and as @jchodera pointed out elsewhere. That's enough to make me back off somewhat from my earlier position. So, I agree with @mattwthompson 's proposal that we update the SMIRNOFF spec to convey that the behavior users have been getting is appropriate for gas phase systems. So, I like this:
Both of these are fine for me.
I think this is accurate. We kind of assume that the default behavior is to get non-periodic systems with no cutoff, and if we want a cutoff and PME we have to explicitly enable. So it probably makes sense to update the spec to reflect this, even though GROMACS has moved away from it.
Right, so in the gas phase with no periodicity or with a very large box there should not be much of an impact on dynamics, conformational preferences, etc. There COULD be dramatic impacts if doing a free energy calculation, such as for the free energy of turning off a molecule's electrostatics in the gas phase (which can be part of a solvation free energy calculation). But for most other cases no. And, this should be something users could avoid -- I mean, in GROMACS (before their updates the last several years) one could deliberately choose to simulate a system in the gas phase with or without periodicity and with or without PME. The main argument for making electrostatics (PME) part of the force field is, I think, for condensed phase systems/periodic systems. But if we can require it there, while not requiring it when a periodic box is not set, that seems like a perfect solution. |
Thanks for tagging me into this issue, and for clearly laying out the problem and potential solutions that attempt to eliminate ambiguity from the specification! There are a few key principles we need to keep in mind as we determine the best way to clarify the spec and ensure all implementations are conformant:
In light of all of this, we should develop a strategy that addresses the following:
My suggestions (only suggestions!) would be to
Regarding gromacs not supporting cutoff: Presumably, specifying |
@mkgilson points out we should also take this opportunity to clarify how net charged periodic systems are treated. |
Thanks all for your thoughtful replies - the discussion here was integral in stewing on the problem and deciding which path to propose as a solution. Earlier I suggested two ideas:
and I decided to move forward with a revised version of #2, more explicitly @j-wags's suggestion of I wrote up an EP that codifies this change and would like to continue the discussion there: #30 |
What should we do about the fact that "PME" is not an accurate description of the electrostatics model we intend (where the true model is Ewald, and methods like PME that perform controlled numerical approximations are acceptable implementations)? |
It's important but I believe the scope and complexity warrants moving it to a separate EP - one that may be as simple as renaming/clarifying what "PME" means, or more prescriptive with respect to the relevant options and implementation details. I've started discussion here: #15 |
I'm not sure that makes sense. If the issue is that the current definition is ambiguous, shouldn't the crisp definition of how nonbonded electrostatics should be handled the purvey of a single EP to address that issue? |
Hm, I think this is letting perfect be the enemy of good. Neither I nor @mattwthompson have a detailed understanding of long-range electrostatics, we largely see this as the minimal change needed to unblock the integration of Interchange into the Toolkit, and make the biopolymer release. I'd be happy to review an alternative EP from @jchodera that handles this in a more complete way. I just don't think that @mattwthompson or I have the background to craft changes that cleanly handle the mathematical/implementation/historical scope of this area. |
I agree with Jeff here. I don't know how to write what's suggested either. |
I know I am late into the discussion, and it looks like it is being resolved adequately. For those interested, here is just a point of clarification on the "boundary conditions" associated with the Ewald. The Ewald formulation has the concept that the infinite periodic system is embedded in a dielectric at infinity. (Strange concept, right?) The "dielectric at infinity" is usually assumed to be conducting, sometimes also called tinfoil boundary conditions. This corresponds to epsilon=infinity. If the boundary is conducting, charges organize in the conductor to produce a field that cancels any total system dipole that might form. Other choices, however, include vacuum, epsilon=1, and I have seen epsilon=80. For finite values of epsilon, there is an extra term in the electrostatic energy expression that must be included. It is always positive and has the square of the total dipole moment of the system. Including this results in a supression of the transient system dipole moments that form during the simulation. Although all of these values (80, 1 and infinity) are, in some sense, reasonable, the choice can have profound effects on the structure and dynamics of the system. It affects the apparent dielectric constant of water, for instance. In biomolecular codes, I have only ever heard of tinfoil boundary conditions. I would vote for supporting just these and to state so in the documentation. |
I am happy with the EP in #30 and I think it resolves this. @jchodera @SimonBoothroyd there are six days left in the feedback window there, so please review so we can get this wrapped up. |
In openforcefield/openff-toolkit#1084 I did a poor job of keeping things in any sort of scope. This issue splits out one issue into something discrete and actionable and focused on just one topic. To resolve this, it will also be helpful to keep the discussion in scope and split out other issues unless they are strictly blocking progress here. I've also attempted to push as much into footnotes as possible.
The problem: PME electrostatics cannot be applied to non-periodic systems
The default and most commonly-used electrostatics method (PME) is fundamentally not compatible with non-periodic systems [1]. This means a force field with
<Electrostatics method="PME" .../>
cannot be used with gas-phase systems, which is a common use case.The history: OpenFF has always modified PME to "no-cutoff" for gas-phase systems
The OpenFF Toolkit, if provided a force field with PME electrostatics and a non-periodic topology, changes the electrostatics method to OpenMM's
Nocutoff
. This is sensible for historical and practical reasons, but modifying the contents of a force field is not desirable and similar implementation in other simulation engines are difficult [2] or not tractable [3].This code is from the early days of OpenFF and as a result all force fields have been trained using this behavior [4]. This should be of no consequence to the vdW refits, but may have an impact on the valence terms for large/charged molecules with relevant
non-bonded intramolecular interactions [5].
This has all led to a situation in which the software and force fields agree what
<Electrostatics method="PME" .../>
means for non-periodic systems, but it's not what the spec says.The question: Should the SMIRNOFF specification be updated to convey this meaning?
I'd argue the answer is "yes," since I don't want implementations to choose between implementing the text of the specification and having
The idea here would be to, in version 0.4 of the Electrostatics section, either
method="pme"
with a clause along the lines of "by the way, if non-periodic, use full electrostatics with no or a very long cutoff"method
flag that conveys the meaning of "PME if periodic, full no-cutoff electrostatics if non-periodic" and make this the default moving forwardIf the answer is "no," the above discrepancy will continue to be baked into our code base and force fields. It will also make it difficult for other implementations to decide, particularly when exporting to non-OpenMM engines, if the specification or implementation should be followed.
There are some unresolved questions [6] that are out of scope here but I want to gather feedback on this specific question before moving forward. If any of these ideas is workable I will move it/them into a proper OFF-EP for consideration.
Coulomb - direct Coulomb interactions (with no reaction-field attenuation) should be used
) but this leads to a mirror of the same problem - full electrostatics for a condensed-phase system is not tractable and a user will simply switch it to PME.Topology
without setting box vectors. Somebody from the fitting team surely has a better understanding of the code and history here.The text was updated successfully, but these errors were encountered: