Could we add an offset term to r in the NonbondedForce?

[zhang-ivy]

We (@dominicrufa, @jchodera ) are exploring implementations for constructing alchemical systems for use in the lab's free energy calculation software Perses.

Currently we are using CustomNonbondedForces with complicated energy expressions so that we can softcore the electrostatics/sterics of alchemical atoms. However, this current implementation is not conducive to running enhanced sampling during the alchemical protocol, so we are trying to re-write the implementation. It would simplify things tremendously for us if we could avoid the need to use softcores and instead, add an offset term to r (let's call this w) that would allow us to prevent singularities.

Here is our proposed feature:
In the openmm.NonbondedForce , we define a potential u_{nb} = u_{elec}(r; {charge}) + u_{steric}(r; {sigma}, {epsilon}). Would it be possible to, when we add a particle as nbf.addParticle(charge, sigma, epsilon) to instead add a particle as nbf.addParticle(charge, sigma, epsilon, w) s.t. u_{nb} = u_{elec}(r + w; {charge}) + u_{steric}(r + w; {sigma}, {epsilon}). Specifically, can we define a mixing rule identical to the epsilon parameter s.t. w = sqrt(w1*w2) ? Can we also support this in addException , addParameterOffset and addExceptionParameterOffset for this ‘lifting’ term w?

[zhang-ivy]

tagging @peastman to get your feedback on this!

[peastman]

One challenge with doing this is how it affects the nonbonded method. PME and LJPME are only derived for potentials of the form C/r^n. They don't apply to potentials of the form C/(r+w)^n.

Possibly one could work around that by using a potential that smoothly interpolates between the modified form at short distances and the standard form at the cutoff. The reciprocal space part would still get computed with the standard form, but we could then correct for it in the direct space part. In principle I think this could work, although the details would need some thought.

[dominicrufa]

Possibly one could work around that by using a potential that smoothly interpolates between the modified form at short distances and the standard form at the cutoff

yeah i think that would work.

The reciprocal space part would still get computed with the standard form, but we could then correct for it in the direct space part. In principle I think this could work

I think that for our purposes, we don't have to have to correct for the reciprocal space part.

[peastman]

For Coulomb, you might be able to get away without correcting. For LJPME the correction would be huge. At short distances, the dispersion energy diverges as 1/r^6. The reciprocal space term is scaled by erf(r), which at short distances is approximately r. So the reciprocal energy will diverge as 1/r^5. I'm not sure it's even possible to correct for that.

Reaction field would also need to be modified. Otherwise, you would get a big discontinuity at the cutoff.

[jchodera]

Reaction field would also need to be modified. Otherwise, you would get a big discontinuity at the cutoff.

We've been working on a reaction field variant that is much better behaved, following this paper. We're doing some preliminary testing right now, but it's likely we may want to revisit the reaction field treatment we support in OpenMM.

Essentially, we're looking for the ability to control the r_eff used in direct-space treatment(s) so that we can "lift" some interactions into the fourth dimension a bit as a function of global parameter(s). We would not need to "lift" reciprocal space interactions, but it does help if we can scale or modify groups of charges and epsilons.

@zhang-ivy and @dominicrufa are currently prototyping something similar in Custom*Forces, but it is likely worth brainstorming what could be supported at the C++ API layer to greatly simplify free energy codes.

[peastman]

Is there a paper describing the method you want to implement? Also, just to make sure I understand,

It would simplify things tremendously for us if we could avoid the need to use softcores and instead, add an offset term to r (let's call this w) that would allow us to prevent singularities.

This is still a softcore potential, right? You just mean that you want to switch to a softcore potential with a different functional form?

Essentially, we're looking for the ability to control the r_eff used in direct-space treatment(s) so that we can "lift" some interactions into the fourth dimension a bit as a function of global parameter(s).

I assume you're speaking figuratively when you talk about the fourth dimension? If you want to literally apply an offset in a fourth dimension, you would replace r by sqrt(r^2+w^2). (I also talked to Yutong about this a while back when he was interested in literally running four dimensional simulations. I don't know what he's done with that since.)

[zhang-ivy]

Is there a paper describing the method you want to implement?

Yes! Here it is: https://aip.scitation.org/doi/10.1063/1.1946750 (I've also attached a PDF in case you don't have access)

Re your last 2 questions, the answer is yes to both. I think Yutong has been running the fourth dimension simulations using his timemachine engine.

rodinger2005.pdf

[peastman]

Thanks! That gives me a much clearer idea of what you want to do. So basically you're going to integrate w from 0 (the solute is fully coupled to the rest of the system) to infinity (it's fully decoupled)?

This means that w is not just a small correction to avoid singularities. You'll be running simulations with large values of w that dramatically change the interaction strength all the way out to the cutoff and even beyond?

I gather you'll need to compute derivatives of the energy with respect to w?

Is a single value of w really sufficient? For example, would you ever want to do calculations of transforming one molecule into a different one, so you simultaneously add some atoms and remove others?

We would not need to "lift" reciprocal space interactions

Can you explain why? If you don't modify the reciprocal space part, the molecules will remain strongly interacting even at w = infinity.

[zhang-ivy]

Thanks! That gives me a much clearer idea of what you want to do. So basically you're going to integrate w from 0 (the solute is fully coupled to the rest of the system) to infinity (it's fully decoupled)?

I don't think we would follow the paper exactly -- I think would interpolate between w = 0 and w = a small value (e.g. 1), so that w is just a small correction to avoid singularities.

I gather you'll need to compute derivatives of the energy with respect to w?

I don't think we'll need to compute the derivative of the energy w.r.t. w for our relative free energy calcs, but that may be something useful for other users?

Is a single value of w really sufficient? For example, would you ever want to do calculations of transforming one molecule into a different one, so you simultaneously add some atoms and remove others?

We do want to do calculations involving transforming one molecule into another, but I think we could add offsets to control the value of w for different sets of atoms.

Can you explain why? If you don't modify the reciprocal space part, the molecules will remain strongly interacting even at w = infinity.

Tagging @jchodera or @dominicrufa to answer this one

[dominicrufa]

So basically you're going to integrate w from 0 (the solute is fully coupled to the rest of the system) to infinity (it's fully decoupled)?

So we don't want to do this per se. really, we just want to take w to some value >0 temporarily while we push q -> 0 and eps->0 (or in the opposite direction) so that we can avoid energies blowing up at r=0.

I gather you'll need to compute derivatives of the energy with respect to w?

strictly speaking, no, but it might be useful if it isn't too difficult.

Is a single value of w really sufficient? For example, would you ever want to do calculations of transforming one molecule into a different one, so you simultaneously add some atoms and remove others?

a single value of w should be sufficient if we can add them as ParameterOffsets and tether them to Global Parameters, i think (as ivy mentions)

[peastman]

So you're essentially using the same protocol you've been using? All you're changing is the particular form of the softcore interaction. It will go something like this:

• Change w from 0 to 1.
• Change q and eps from their initial values to their final values.
• Change w from 1 back to 0.

strictly speaking, no, but it might be useful if it isn't too difficult.

If w is a global parameter, we can compute derivatives with respect to it. If it's specified per-particle, we don't have any way to calculate it.

[dominicrufa]

So you're essentially using the same protocol you've been using? All you're changing is the particular form of the softcore interaction. It will go something like this:

yes, that's about right.

If w is a global parameter, we can compute derivatives with respect to it. If it's specified per-particle, we don't have any way to calculate it.

it might be better in our case if w is a per particle parameter and we can add an offset to it with a ParameterOffset with a global parameter multiplier.

[zhang-ivy]

@peastman Thanks for the brainstorming call today! We talked about a lot possible ways to make our alchemical factory implementation (that allows for rest scaling) easier, but just to summarize -- the main features we think will be feasible:

1. Separate the direct space from the reciprocal space expressions in the NonbondedForce
2. Add a w term offset to r (as a per particle parameter, potentially also multiplied by a global parameter lambda) in the direct space part of the NonbondedForce

Does this sound right to you? If you are able to think of alternate ways to make our alchemical factory implementation (with rest scaling) easier, feel free to share them here and we can discuss more.

[peastman]

Item 2 would be done with a CustomNonbondedForce. I'll also look into features that could make that easier, like being able to completely disable the direct space part of NonbondedForce so you don't need to worry about specifying the right groups in every call to getState().

[zhang-ivy]

Sounds great! Any updates on whether it'll be possible to separate the direct space from the reciprocal space in the NonbondedForce?

[peastman]

It's definitely possible to separate the direct space and reciprocal space parts. Right now you do it by assigning them to different force groups. Allowing you to completely disable one of them will make it easier, though, since you won't have to worry about passing the right flags to getState().

[zhang-ivy]

Right sorry, I meant to ask whether you were able to figure out a way to completely disable one of them to make it easier, but no worries if its still in progress

[zhang-ivy]

Dominic and I have started drafting the implementation that we discussed (using the CustomNonbondedForce) in the meeting -- I will share it sometime in the next couple days, but just wanted to check in on the status of disabling the direct space in the NonbondedForce!

[peastman]

I haven't started on that yet. I have some other things I need to deal with first.

[peastman]

When we met to discuss this, we decided not to put it into NonbondedForce but instead to use a CustomNonbondedForce for the direct space calculation along with the reciprocal space calculation from NonbondedForce. I added a feature to make that slightly easier in #3173. Are there any other code changes still needed to support this, or can the issue be closed?

[zhang-ivy]

I've implemented a version of the HybridTopologyFactory class that uses a CustomNonbondedForce for the direct space. I will chat with Dominic later this week to review my implementation and will tag you for review afterwards. We can close this issue for now as we work on this and I can open another issue if we are still having problems.