How to handle constraints in a self-consistent way with `SMIRNOFFTemplateGenerator`?

[epretti]

I'm still trying to work on #423. Though it wasn't the original focus of that PR, testing virtual site handling with the SMIRNOFF ports of various water models has turned up some unpleasant edge cases involving constraints. Specifically, I've found broken behavior where waters end up with neither rigid constraints nor harmonic bond/angle terms.

I'm opening this issue primarily as a request for comments to establish what the behavior should be, and the implementation details can follow after. I think there's a mismatch between how OpenMM and OpenFF treat constraints. OpenMM treats them entirely at system creation time, controlled by the constraints, rigidWater, and flexibleConstraints options. What this means is that constraints aren't part of an OpenMM XML (atom typed) force field.¹ Instead, the lengths and angles of bond/angle terms in the force field are used to add constraints that the user requests. SMIRNOFF force fields, on the other hand, contain <Constraints> tags to explicitly specify constraints, and there is no equivalent to the constraints and rigidWater options if you are creating an OpenMM System using the OpenFF Toolkit or Interchange (@j-wags @mattwthompson please correct me if that's wrong). Furthermore, for the SMIRNOFF water force fields with constraints, there aren't any harmonic bond/angle terms present.

What should the behavior be if a user's constraint options to ForceField.createSystem() request that certain bonds or angles be flexible, but there are only rigid constraint parameters in the SMIRNOFF force field for them?

• Set some default force constants for the bonds and angles² silently
• Set some default force constants, and warn about it
• Add the rigid constraints anyway, silently
• Add the rigid constraints anyway, and warn about it
• Stop and raise an error

Since OpenMM treats water as a special case for constraints, we'd need to detect that a user is trying to use SMIRNOFF water (which is relatively easy). Will water be the only case for SMIRNOFF constraints, or do/will other OpenFF force fields use them too?

Footnotes

1. There is one place where XML force fields have information related to constraints: if a residue template for water (which must be named HOH) has a rigidWater attribute with a value other than true, and the rigidWater user option is not specified, this instructs OpenMM to make flexible water molecules. ↩

2. The XML force fields distributed with OpenMM for rigid water models do this in some sense; many contain bond and angle force constants equal to 1106 kcal/mol/Å^2 and 200 kcal/mol/rad^2. Where did these numbers come from? ↩

[peastman]

Where did these numbers come from?

I believe they came from Gromacs. Compare the force constants in https://github.com/openmm/openmmforcefields/blob/main/amber/files/tip3p.xml to https://github.com/gromacs/gromacs/blob/5208be5adafc1023c94db6c0a360e9f3cd41d158/share/top/amber99sb.ff/ffbonded.itp#L16 and https://github.com/gromacs/gromacs/blob/5208be5adafc1023c94db6c0a360e9f3cd41d158/share/top/amber99sb.ff/ffbonded.itp#L124.

[peastman]

A few comments on the use of flexible water. Each water model is parameterized to be either rigid or flexible, and should only be simulated that way. No one should ever try to run a simulation with, for example, flexible TIP4P. That isn't really TIP4P and won't give good results. Likewise, no one should try to simulate rigid AMOEBA water. That just isn't AMOEBA water.

On the other hand, it's really useful to make water flexible for certain purposes. For example, energy minimization is a lot easier and more robust when there are no constraints. So you might minimize your system with flexible water, then switch back to rigid water for the simulation.

The same is true with OpenFF. The water model is designed to be rigid and should only be simulated that way. But it's still useful to be able to make it flexible.

[epretti]

If water is the only case where a SMIRNOFF force field might have constraints, we can have special handling for it. In this case, I think I'd be in favor of adding default spring constants to the generated templates so that users can still get flexible water for, e.g., minimization if desired.

The question is then what to do for molecules other than water to which a SMIRNOFF force field adds <Constraints>. Or, if this shouldn't be the case for any OpenFF force fields now or in the forseeable future, we can just raise a NotImplementedError or something if we find it, leaving open the option to support it down the road if necessary. Right now, SMIRNOFFTemplateGenerator silently deletes the section and ignores the constraints.

[mattwthompson]

Your understanding of the SMIRNOFF specification is correct, below the horizontal line I'm mostly repeating the last two sections from here, if something is inconsistent I'm probably wrong.

If I understand the problem correctly, the hope is to have a design which both

• faithfully implements constraints as defined by parameters in SMIRNOFF force field files(s)
• produces a system with constraints matching the constraints to openmm.app.ForceField.createSystem()

Unfortunately I think these can't be achieved at the same time without errors being thrown for many common use cases. For example, using a default mainline OpenFF force field and using the default of createSystem(..., constraints=None) are in conflict.

Since OpenMM treats water as a special case for constraints

So I understand this correctly, are you talking about the application or library layer? I'm pretty sure in Interchange I'm handling water and non-water constraints the same (just `my_system.addConstraint()1, I think). If I should be doing something different I will open a separate issue.

---

Constraints are part of the force field so openff.toolkit.ForceField.create_openmm_system cannot have a constraints argument as found in openmm.app.ForceField.createSystem. A flexible water model¹ could be implemented by dropping the Constraints parameters and adding in harmonic bond and angle parameters like so. Peter is right about flexible vs. rigid water; the difference is that OpenFF implement of that toggle would involve loading a different file/object in, not toggling an argument to a function.

Also note that the water parameters found in Sage are just TIP3P and all water models shipped as .offxml files are just ports. OpenFF may refit a water model in the future but that's some ways out and likely to be a 3- or 4-site model with tuned parameters, not a whole new style of model altogether.

Will water be the only case for SMIRNOFF constraints, or do/will other OpenFF force fields use them too?

The so-called "mainline" force fields have constraints for all bonds between hydrogen and heavy atoms. (There are "unconstrained" versions of each that strip this parameter out.) This has been the case since Parsley 1.0.0 and I believe all force fields in development are intended to carry this through. Compare here and here

In terms of how this works, Sage 2.0.0 applied to an ethanol molecule has 6 constraints:

>>> from openff.toolkit import ForceField, Molecule
>>> ForceField("openff-2.0.0.offxml").create_openmm_system(Molecule.from_smiles("CCO").to_topology()).getNumConstraints()
6

no more and no less. Same goes for something that combines this with a different force field-topology pairing in a larger system. Nothing's stopping users from removing all of them or adding other constraints, but at that point it becomes something different than Sage 2.0.0 applied to ethanol.

Lastly, people outside of OpenFF could also define their own constraints to chemical systems that I'm not creative enough to list off the top of my head right now. The specification is here, but so far has only been used in the above two cases. Interchange implements every case I can think of, I doubt constraining other types of bonds would end up being so different.

Footnotes

1. Though TIP3P is canonically rigid and I don't know the ultimate source of bond and angle force constants ↩

[epretti]

Since OpenMM treats water as a special case for constraints

So I understand this correctly, are you talking about the application or library layer? I'm pretty sure in Interchange I'm handling water and non-water constraints the same (just my_system.addConstraint(), I think). If I should be doing something different I will open a separate issue.

Just the application layer (specifically ForceField.createSystem()). From the library layer, there's nothing special about them, so what you're doing shouldn't be wrong there.

[j-wags]

Oh, this is a good edge case to discuss. Thanks for writing this up @epretti.

(I'm largely reiteration Matt's points here but posting for the sake of completeness)

SMIRNOFF force fields, on the other hand, contain tags to explicitly specify constraints, and there is no equivalent to the constraints and rigidWater options if you are creating an OpenMM System using the OpenFF Toolkit or Interchange

This is correct

If water is the only case where a SMIRNOFF force field might have constraints

Just being precise here - The most-used version of our flagship force fields (eg openff-2.2.1.offxml) apply constraint parameters to all chemical bonds involving Hs (inheriting the length from bond parameters that match the same chemical bond), as well as special constraints with pre-defined distance for water Hs.

The SMIRNOFF way of deciding what to do with bonds/constraints is:

1. If a (constraint parameter WITH a defined distance) matches the bond, apply only that, even if a harmonic bond parameter also matches. Otherwise...
2. If (a constraint parameter WITHOUT a defined distance) and (a harmonic bond parameter) both match the bond, apply a rigid constraint with the length from the harmonic bond. Otherwise...
3. If only (a harmonic bond parameter) matches the bond, use the harmonic bond parameters. Otherwise...
4. Raise a missing parameter error

Right now, SMIRNOFFTemplateGenerator silently deletes the section and ignores the constraints.

This would run into trouble if a SMIRNOFF force field has constraints defined using case 1 above. This is the case for water models in our mainline FFs.

My two thoughts here are:

• Users could use the unconstrained variant of our mainline FFs (eg. openff_unconstrained-2.2.1.offxml) when calling them through OMMFFs. This still has rigid water, but the non-water bonds to hydrogens would all have harmonic bond parameters, and OMMFFs could decide what to do with those based on which constraints the user requested.
• If OMMFFs gets a rigid water from SMIRNOFFTemplateGenerator, but the user wants it to be flexible, I'd recommend raising an error. While the water model is currently TIP3P and maybe OMMFFs has some precedent for "flexible TIP3P", we may change that at any time, and if someone was trying to experiment with a different water model I'd prefer they get an error rather than silently getting some other parameters substituted in.

[epretti]

Thanks for enumerating those cases. In case 2, will the OpenFF Toolkit still create the harmonic bond terms in an OpenMM System for bonds with rigid constraints? If so, I think that's why openmmforcefields is working for this case presently. Right now, case 1 is what is entirely broken.

I don't know if we can meaningfully detect at template generation time what kinds of constraints a user is requesting, since the generator will register a template in the ForceField that might get used in a later call to createSystem that has different constraint options.

[j-wags]

In case 2, will the OpenFF Toolkit still create the harmonic bond terms in an OpenMM System for bonds with rigid constraints?

No, it wont. Here's an example using ethyne:

from openff.toolkit import ForceField, Molecule
mol = Molecule.from_smiles('C#C') # ethyne
ff = ForceField('openff-2.2.1.offxml')
sys = ff.create_openmm_system(mol.to_topology())
sys.getNumConstraints()

2

This is the two C-H bonds getting constraints

bond_force = sys.getForce(3)
bond_force..bf.getNumBonds()

1

This is the C-C triple bond getting a harmonic bond parameter.

But I think OMMFFs gets harmonic bond parameters by deleting all the constraints . If I do the same thing in a local example it yields harmonic bond parameters for all three bonds as expected:

ff.deregister_parameter_handler("Constraints")
sys = ff.create_openmm_system(mol.to_topology())
bf = sys.getForce(3)
bf.getNumBonds()

3

sys.getNumConstraints()

0

However, as you point out, this is bad for waters, since then they get the bond and angle parameters from the general force field instead of the TIP3P geometry that the constraints encoded.

[mattwthompson]

Jeff's right and I'm not sure the below is necessarily useful (more complexity and cases to enumerate). But still, just in OpenFF world, harmonic bond and angle parameters can be added in even when they're ~pointless due to being overridden by constraints. It's just not the default behavior¹

Building on Jeff's code:

>>> sys = ff.create_interchange(mol.to_topology()).to_openmm_system(add_constrained_forces=True)
>>> sys.getNumConstraints()
2
>>> sys.getForce(3).getNumBonds()
3

i.e.

1. C#C has an (unconstrained) harmonic bond
2. each C-H also has a harmonic bond, but overriden by a constraint (with a distance the same as the equilibrium bond distance)

double-checking:

>>> for index in range(2):
...     sys.getConstraintParameters(index)
...
[0, 2, 0.10709717863209998 nm]
[1, 3, 0.10709717863209998 nm]
>>> for index in range(3):
...     sys.getForce(3).getBondParameters(index)
...
[0, 1, 0.12142555389339998 nm, 972533.4709478018 kJ/(nm**2 mol)]
[0, 2, 0.10709717863209998 nm, 389978.3886748757 kJ/(nm**2 mol)]
[1, 3, 0.10709717863209998 nm, 389978.3886748757 kJ/(nm**2 mol)]

Footnotes

1. Going on memory, Peter told us there wasn't much of a performance hit associated with including these but we also decided there was literally no case in which these parameters should be added in. So we left them out ↩

[epretti]

Thanks for that example demonstrating the actual reason why "case 2" is working at present.

One possible approach could be as follows:

1. Use Interchange with add_constrained_forces=True to create the reference System that a template will get built from.
2. For each constraint that gets added to the force field, look for an associated bond or angle term.
   a. If the geometry (length/angle) from the term matches the constraint, add the term and its force constant to the template.
   b. If it doesn't match, or if there isn't a term, and the molecule is not water, raise an error.
   c. For water, as a special case, when constraints are found without terms, put terms into the template with the appropriate lengths/angles determined by the constraints, and the default force constants from OpenMM's water models.

This approach does special-case water, but OpenMM's ForceField already does. Application of constraints would be controlled by OpenMM as usual; constraints=HBonds would give behavior matching the constrained OpenFF variants and the workflows of existing users not already using the unconstrained force field variants wouldn't break.

I thought there would be no way to raise an error from a TemplateGenerator if a user asks for flexible water. I realize you could embed Python as a string inside a <Script> tag in the XML for the template that gets cached (I think it'd work but haven't tested). This is an ugly solution though, and also would prevent the use of flexible water for minimization.

[j-wags]

Following our discussion earlier this week about the tension between

• The OpenMMForceFields API allows users to request SMIRNOFF water without constraints
• No current OpenFF force field has harmonic bond/angle parameters intended for water (while some harmonic bond/angle SMIRKS will match water, they are always superseded by rigid constraints in the same file)

From OpenFF's side it's quite important that OpenMMForceFields users never think they're getting flexible water parameters from our FFs, as vanilla usage of our FFs can't yield flexible water parameters, and the harmonic parameters that do match are not intended for water.

I think the solution to getting OpenMMForceFields' behavior right might lie in the fact that we release two versions of each force field, eg openff-2.3.0.offxml and openff_unconstrained-2.3.0.offxml. The _unconstrained version of each FF can be thought of "the minimal number of constraints that this FF can be used with". It omits small molecule hydrogen bond constraints (since we actually fit those harmonic bond parameters as part of our training), but still contains constraints to make waters rigid since we don't intend to provide harmonic bond parameters for them.

So a workable solution here might be:

• If the user requests a SMIRNOFF force field, try to find the version with minimal constraints. This could be accomplished by:
  • checking whether a version of the same FF but with _unconstrained in the name exists and switching to it (a little automagic but probably harmless)
  • instead of removing the entire constraints section in OMMFFs, just remove the constraint named c1, which will leave the TIP3P constraints in place (somewhat brittle, but this is in practice how we generate the _unconstrained variants now, and we're likely to stay with this convention for a while, and it's failsafe (if we name the constraint something else then users will get errors instead of incorrect parameters))
  • raise an error telling the user to directly call the _unconstrained version (clearest about demanding user intent but will break basically all current user workflows)
• Use the "unconstrained" SMIRNOFF FF for parameter assignment
• When OMMFFs sets constraints according to its keywords, never let it remove a constraint that was left there by the FF (only ever let it keep or add them)

I think this should give the "correct" behavior in all cases.

[epretti]

If the user requests a SMIRNOFF force field, try to find the version with minimal constraints.

We could do this; since users could pass filenames to custom OFFXML files and aren't restricted to using OpenFF builtin forcefields, I'd prefer to do this directly based on the constraints rather than by looking at the FF name.

When OMMFFs sets constraints according to its keywords, never let it remove a constraint that was left there by the FF (only ever let it keep or add them)

The problem is that this is not under the control of OMMFFs, but instead OpenMM (when the user calls createSystem and gives keyword arguments). As is, the residue templates generated by OMMFFs and given to OpenMM do not contain information about constraints (and OpenMM doesn't provide functionality for them to have this information apart from flexible bond and angle terms). I understand that you don't want OMMFFs to make up harmonic parameters for SMIRNOFF water. In that case, do we think that a change to OpenMM is necessary?

[j-wags]

Ahhh, I think I see. So you're saying that, at template creation time (which happens in OpenMMForceFields), the TemplateGenerator doesn't know which constraint settings the user passed into the createSystem call (which is in OpenMM itself), so it can't do the policing of constraints. And then once the template's made, OpenMM doesn't know which of the bonds would have been constrained by the SMIRNOFF FF because all it sees are harmonic bonds in the template.

I've tried to reason through what the current state of OpenFF water in OpenMMForceFields is, and I think it's some flavor of "not correct" (either waters get no bond parameters because they're overridden by the constraints, or they get hydroxyl parameters that aren't intended for water). So with that as precedent, I think it'd be ok to do something like unconditionally raising an error in this case.

Looking ahead, I think that adopting a behavior like

• switch to the _unconstrained version of the SMIRNOFF FF if one can be found [footnote 1]
• parameterize the molecules with add_constrained_bonds=False
• have SMIRNOFFTemplateGenerator raise an error if a chemical bond in the reference molecule doesn't get a harmonic bond potential assigned (since that would mean it was overwritten by a constraint that was still in the _unconstrained version of the FF, indicating that there were no harmonic parameters suitable for application)

should do the right thing for the cases I can think of. Specifically:

• It would raise errors for waters currently
• If some day we released a flexible water model, this would be reflected by the removal of water constraints in the _unconstrained FF, so the flexible parameters could pass through

[footnote 1] It's tempting to think that removing the c1 constraint from a SMIRNOFF FF would work more easily than looking for the _unconstrained version of the file, but that will hit a problem if we release a FF with flexible water, since the constrained version would need an H-H constraint to fix the angle, and that angle-setting constriant would have a name different than c1

[peastman]

ForceField can also constrain angles based on harmonic bonds and angles. It uses the equilibrium angle to figure out what the constrained distance should be.

I don't think we need to care what the constraints are named. If the forcefield says a distance should be constrained, and there's no corresponding bond definition, throw an exception.

[j-wags]

I don't think we need to care what the constraints are named. If the forcefield says a distance should be constrained, and there's no corresponding bond definition, throw an exception.

Just double checking that this is referring to a future where the _unconstrained version of the FF is permitted and switched to if it's available. Otherwise, no main-line OpenFF FF could be correctly handled by OpenMMForceFields, since

• openff-2.3.0.offxml constrains all bond distances to hydrogen (and if add_constrained_bonds=True is specified to OpenFF during system creation, then small molecule bonds to Hs would get correct harmonic bond parameters, but waters would get nonsense harmonic parameters)
• openff_unconstrained-2.3.0.offxml would be forbidden by this logic.