Proposal to add an in-plane angle parameter to the DivalentLonePair virtual site type

[tom-potter-cresset]

No description provided.

[j-wags]

As a start for discussion - I think a minimum requirement for approving things should be whether we can implement them in OpenMM and this seems to pass that check. On a technical level I feel reasonably confident this will generally be possible to export to OpenMM without any nasty edge cases/implementation details (though @mattwthompson can correct me if I'm overlooking something), and @tomdpotter mentions that this type of vsites is also supported in GROMACS.

Conceptually, the thing that gives me pause here is the... possible asymmetry (for lack of a better term?) of these sites, in combination with possible unexpected symmetries of molecules that might come in and the (admittedly byzantine) deduplication logic for when multiple vsites affect the same atom(s). For example, a seemingly straightforward vsite for hydroxyls (like is present in the current version of this PR) would create asymmetric water, since the ordering of the hydrogens is arbitrary and the match=once keyword is specified:

<VirtualSite
    type="DivalentLonePair"
    name="VS1"
    smirks="[#8X2:1](~[*:2])~[#1:3]"
    distance="0.4 * angstrom"
    charge_increment1="0.1*elementary_charge"
    outOfPlaneAngle="60.0 * degree"
    inPlaneAngle="20.0 * degree"
    match="once" >
</VirtualSite>
<VirtualSite
    type="DivalentLonePair"
    name="VS2"
    smirks="[#8X2:1](~[*:2])~[#1:3]"
    distance="0.4 * angstrom"
    charge_increment1="0.1*elementary_charge"
    outOfPlaneAngle="-60.0 * degree"
    inPlaneAngle="20.0 * degree"
    match="once" >
</VirtualSite>

I don't know off-hand of a general way to police against this/warn users when this happens. On one hand, we could say that FF developers are free to do risky stuff like this and don't need guardrails. On the other hand, maybe there are some common-sense things that we can do to structurally avoid these cases?

The first thought that comes to mind is to evaluate whether in_plane_angle is nonzero and create two virtual particles if so, but as @trevorgokey humbly summarized after months of thinking about vsites, "The toolkit never compares the physical numbers to determine equality as this can lead to instability during e.g. parameter fitting.". So I don't think I'd approve a solution that looks at whether in_plane_angle is nonzero to decide how many virtual sites/particles to make. Maybe we should be thinking about creating an entirely new type of vsite, like AsymmetricDivalentLonePair, that by construction expects nonzero in_plane_angle values and will always make pairs of virtual sites/particles.

[tom-potter-cresset]

Good point about the hydroxyl example. The smirks should really be something like "[#8X2:1](~[!#1:2])~[#1:3]" to prevent that issue, as we don't intend to cover water with that site definition.

I assume that if an AsymmetricDivalentLonePair type was added, the sulfur sigma hole use case would still need to be covered by an extended DivalentLonePair definition? For those, we're looking at something like this dimethylsulfide example:

[lilyminium]

I support extending the spec to allow for both these functionalities. I agree with Jeff though that how to do it might require a bit of workshopping. Just on the technical side --

The first thought that comes to mind is to evaluate whether in_plane_angle is nonzero and create two virtual particles if so, but as @trevorgokey humbly summarized after months of thinking about vsites, "The toolkit never compares the physical numbers to determine equality as this can lead to instability during e.g. parameter fitting.". So I don't think I'd approve a solution that looks at whether in_plane_angle is nonzero to decide how many virtual sites/particles to make. Maybe we should be thinking about creating an entirely new type of vsite, like AsymmetricDivalentLonePair, that by construction expects nonzero in_plane_angle values and will always make pairs of virtual sites/particles.

It's worth noting that the toolkit already compares whether the out_of_plane_angle is nonzero to determine whether match="once" is supported upon loading the FF to avoid asymmetries, so this might not be out of the question, and naively I would expect you could apply this to in_plane_angle too. So one option is potentially to only allow match="once" for the extended DivalentLonePair if both in_plane_angle and out_of_plane_angle are 0, and force users to define their SMIRKS uniquely if they specify an angle but only want a single match. If I'm reading the code correctly I think that's how the current DivalentLonePair works in practice.

[davidlmobley]

In general I agree with this approach and think this would be a worthwhile change.

In terms of backwards compatibility I think this makes sense:

Alternatively, backwards compatibility could be kept by making inPlaneAngle an optional parameter, with a default value of 0.

as we probably have a few test force fields in the wild using this and probably don't want to break them unless there is a good reason to do so (which this doesn't seem to be).

I also agree that this could cause issues if match="once" is allowed, so I think I agree with Lily's proposed approach:

It's worth noting that the toolkit already compares whether the out_of_plane_angle is nonzero to determine whether match="once" is supported upon loading the FF to avoid asymmetries, so this might not be out of the question, and naively I would expect you could apply this to in_plane_angle too. So one option is potentially to only allow match="once" for the extended DivalentLonePair if both in_plane_angle and out_of_plane_angle are 0, and force users to define their SMIRKS uniquely if they specify an angle but only want a single match. If I'm reading the code correctly I think that's how the current DivalentLonePair works in practice.

Specifically, matching once should only be allowed if both are exactly zero. As noted, this could cause slight issues in fitting, though, e.g. if a fit starts with match=once and then tries to adjust the value away from zero. (Lily's point, I suppose, is that we already have this issue.)

I think we should tackle problems one at a time, and so deal with this enhancement in the current framework (allowing only if zero) and then raise an issue on the tracker to deal with later wherein we might address the numerical issue, or at least just document it. (Alternatively we could just document in the spec and call it a day?)

Sorry for the delay. I'm OK with this going ahead to a formal spec change when ready!

[davidlmobley]

I think we are OK to proceed with this when you are ready, @mattwthompson ?

[mattwthompson]

If the next step is morphing this into a formal EP - yes, please, don't let me stand in the way of that!

I haven't thought through (much less tested) the implementation details I'm likely to be responsible for, but IIUC that's a future problem