Doc string clarification request

[davidlmobley]

trafficstars

In HybridTopologyFactory, there's an option which seems important but for which we can't fully understand the doc string. Specifically:

 interpolate_old_and_new_14s : bool, default False
            Whether to turn off interactions for new exceptions (not just
            1,4s) at lambda = 0 and old exceptions at lambda = 1; if False,
            they are present in the nonbonded force.

Issues: a) What does this have to do with 1,4s? (Implied by the name) b) Why does the doc string say it doesn't apply to just 1,4s? c) What are "old" vs "new"? Are we talking about the two states in the perturbation, e.g. A vs B state? d) Why is this referring to exceptions? e) What's the interpolation implied by the name? ("interpolate" vs "whether to turn off" seem to be two very different things)

Our best guess as to what this is doing/claiming to do is to control whether nonbonded interactions are turned off for dummy atoms. We ran into a case where our dummy atoms were forcing the attached physical/real atoms to stay in a particularly odd conformation (relative to the conformation expected if the dummy atoms are not present) and turning this option to True resulted in the physical atoms returning to the "expected" conformation. My student's interpretation of this was that this is because the nonbonded interactions are otherwise left on (if we don't set this to True) and that results in weird forces/conformations. However, then I had all the questions above and she couldn't answer them.

I'm guessing the fix here is to clarify the doc string and possibly change the name of the option, if it doesn't specifally have to do with interpolation and 1-4s?

cc @Yoshanuikabundi

[IAlibay]

@davidlmobley thanks for opening up this issue.

This is a bit of an obscure part of the HTF, so it might take us a little bit to come up with a complete answer to all your questions. I do suspect that the docstring is at best misleading here.

because the nonbonded interactions are otherwise left on

Apologies, I'm not sure I fully understand this bit. Here do you mean "all nonbonded interactions involving dummy atoms"? Happy to jump on a call to discuss this further if it's easier.

[davidlmobley]

I agree my statement may be unclear, but I doubt it's necessary for me to clarify, since the whole point here is that I don't understand what the doc string is saying so my interpretation of it is also likely incorrect, because I'm having to invoke my imagination and make wild guesses. If you NEED me to get on a call I'm of course happy to do so, but my point is that I'm confused and uncertain, yet this setting is doing something clearly important that we need to understand.

[IAlibay]

@davidlmobley this is a particularly complicated topic, apologies for the rather lengthy response. The OpenFE team did quite the deep dive on this (particularly thanks to @hannahbaumann). Happy to go into more details on any parts of the topic as necessary.

I’ll start by rephrasing some of your questions into the following:

1. What do the references to “old” and “new” mean?
2. Why are we talking about exceptions when dealing with 1,4 interactions?
3. What does this method do?
4. Why is it not limited to just 1,4s?
5. Should we be using interpolate_old_and_new_14s?
6. What are the changes that OpenFE is looking to do here?

What do the references to “old” and “new” mean?

As you correctly guessed, this refers to the parameterized systems of the end states, e.g. old == state A system, new == state B system.

Apologies for the confusion here, it’s some leftover terminology from Perses which we need to clean up.

Why are we talking about exceptions when dealing with 1,4 interactions?

In OpenMM, nonbonded exceptions are the means by which one defines special direct space pair interactions (e.g. not defined by the parameters of the individual particles). Specifically, this is how you set 1-2, 1-3, and 1-4 interactions.

Exceptions are not limited to 1,N interactions, but are also used for other things like virtual sites.

What does this method do?

Definitions

Unique atoms Atoms which are unique to a given end state and alchemically disappear to dummy atoms in the alternate end state

Core atoms Atoms which are “mapped” between the two end states and are alchemically interpolated

Environment atoms Atoms which are the same between the two end states and are not alchemically interpolated

When interpolate_old_and_new_14s is False

For the following parts of the system, the following happens to the Coulombic and LJ Exceptions within the NonbondedForce:

1. Between “unique atoms of different states” Interactions are fully excluded, e.g. there are no direct space interactions between unique state A atoms and unique state B atoms.
2. “Unique atoms of the same state” Exceptions for those atoms are taken from the corresponding endstate system
3. “Unique atoms interacting with shared core atoms” Exceptions are taken from the relevant unique atom end state system
4. Between core-core and core-environment atoms Exceptions for given atom pairs are taken from each end state system and set to a midpoint value.
5. Between environment atoms Exceptions are taken directly from the state A endstate system

In all these cases, the exceptions are kept constant throughout the transformation. This specifically means that in the case number 3, 1-4 interactions are still applied between core and unique atoms even when the unique atoms should not be interacting with the rest of the system. This is probably what is leading to the odd conformations you are observing.

Note: this is only affecting direct space interactions - reciprocal space coulombics are not affected by exceptions and charges are scaled through lambda as appropriate.

When interpolate_old_and_new_14s is True

In this case, the Exceptions are handled by both a standard NonbondedForce and a custom BondForce.

The same exceptions are extracted as in the interpolate_old_and_new_14s=False case, however the “unique-unique of the same state” and “unique-core” exceptions are lambda scaled through the BondForce.

Why is it not limited to just 1,4s?

Exceptions don’t have tags which let you know what they come from, so the HTF is a bit heavy handed and assumes that all exceptions need handling. This is fine in the case of 1-2 and 1-3 interactions which are set to zero, but may be problematic for custom exceptions.

Should we be using interpolate_old_and_new_14s=True?

This is a set of methods inherited from Perses and as far as we know haven’t been particularly well validated.

Specifically:

• We haven’t run full benchmarks using this option - this is something we will try to do now.
• We aren’t 100% sure about the reciprocal space correction, which we need to discuss with the Perses team.

edit: It would however seem that, if working correctly, it would provide a better behaved transformations.

What changes are needed here?

There appears to be quite a lot of work to do here. Some initial thoughts include;

• We should update documentation
  • The HTF docstring
  • The settings name & description (maybe switch from interpolate to annhilate? or switch off?)
• Validate the setting
  • We’ll re-run some benchmark systems to check its use
  • Investigate with relevant parties any concerns about reciprocal space corrections
• Change our settings default as necessary

[hannahbaumann]

We renamed interpolate_old_and_new_14 to turn_off_core_unique_exceptions and changed the doc string in the relative free energy protocol setting https://github.com/OpenFreeEnergy/openfe/blob/56816d347939e604f23cf8a245946a0618eba00e/openfe/protocols/openmm_rfe/equil_rfe_settings.py#L85-L90

We haven't changed the name in the hybrid topology factory itself yet, but will discuss doing so soon! Please let us know if this is still unclear!