Nonbonded force between atom groups

[peastman]

This is a commonly requested feature, so I'd like to open a discussion on how it should work. The idea is to create a nonbonded force that, instead of letting every atom interact with every other atom, only considers the interaction between the atoms in one group and the atoms in another group. There are lots of uses for this. Here are some examples:

• For analysis. For example, you might want to compute the interaction energy between solute and solvent, but exclude all interactions between two solute atoms or two solvent atoms.
• For efficiency. Suppose you're running a simulation where everything is held fixed except a small number of atoms around the active site. Most of the interactions (anything between two immobile atoms) will never change, so recomputing them every time step is inefficient. It might be faster if you could only compute the interactions that involve at least one atom that can move.
• Certain types of restraints, such as for NMR refinement, take this form. You have two small clusters of atoms, and you want to compute an interaction between every atom in the first cluster and every atom in the second cluster.

What other uses would it have?

An obvious solution is to create a version of CustomNonbondedForce that lets you specify pairs of atom groups. For each pair of groups you specified, the interaction would be computed between every atom in the first group and every atom in the second group. (But if an atom appeared in both groups, presumably the self interaction should be skipped?) This could be a new custom force, or just a new feature added to CustomNonbondedForce so it could work in either mode (all pair interactions, or only interactions between groups).

[jchodera]

We're extremely interested in this feature for alchemical free energy calculations.

Currently, we create alchemically-modified intermediates by setting the epsilon and charge atom properties of the ligand to zero in NonbondedForce and creating a new CustomNonbondedForce to handle the interactions between the ligand (usually a few tens of atoms) and the rest of the system (tens of thousands of atoms). To ensure only interactions between the ligand and the rest of the system (and, optionally, within the ligand) are computed, we add an additional alchemical property that only computes interactions between dissimilar types.

With the proposed "Nonbonded force between atom groups" (let's call it CustomNonbondedGroupForce for purposes of discussion), we would instead define a force between the ligand atoms in one group and all other atoms in another group.

If some atoms belong to both groups, it does make sense to permit the interactions between non-identical atoms in this group to be computed, though this would require that exclusions can be specified, introducing additional complexity.

The ability to assign difference Force objects to different force groups and extract the total potential energy and the potential energy of each force group individually would also be useful.

[peastman]

It should definitely be possible to have an atom in both groups. Consider the example where you want to skip computing interactions between atoms that don't move. In that case, one group would include only the mobile atoms, while the other group would include all atoms (whether mobile or fixed). Exclusions would work exactly the same as before.

[peastman]

Does anyone else have opinions on this? Thoughts about the API? Particular features it needs to support?

[kyleabeauchamp]

So I've briefly thought about using something like this to implement predictors of chemical shifts. A lot of the terms in chemical shift prediction look similar to nonbonded forces, but I'd need a way to access atom-specific energies to do this.

I don't think this is exactly the same as what other people are envisioning, but I thought I'd discuss this as a possible application.

[peastman]

Could you describe in more detail what you would want to do? What would the calculation look like, and what information would you want to get out?

[kyleabeauchamp]

Actually, right now I'm leaning towards not doing this in OpenMM--I think I'd be better off just doing this in Python. However, there are some things that OpenMM has that could be "useful" in this sort of analysis (creating virtual sites, for example). I'd say you should ignore my idea for now, it's a pretty tangential application...

One example is ring current effect, which looks something like this:

Magnetic anisotropy looks something like this:

My idea was to try to encode such energies as custom forces and use OpenMM to extract them from trajectories of conformations.

This is pretty far fetched, but it's a way to use some of OpenMM to simplify our code.

[peastman]

Nothing particularly far fetched about those functions. But would you compute them for all possible atom pairs? Or only a limited set? Or just one? And would you want to get out the value for every pair, or some sort of sum?

[kyleabeauchamp]

So the ring current would be between a given atom (we'll probably want all atoms though) and the sum over all ring atoms in a protein. The S_{qp} term actually refers to a projection of a given atom onto the plan of a particular ring.

For the anisotropy, again I would be looking at the interaction between a given atom (repeated for all atoms in protein) and summing over all other sites in the protein, where a "site" is probably some virtual site that we place at the center of the peptide bond.

[peastman]

It would be easy enough to implement that calculation with this feature. You'd just ask for the interaction of one atom with every other site of the desired type, then compute the energy. But you'd need to rebuild the Context to repeat the calculation for the next atom, so it would be a really inefficient way to implement it.

[jchodera]

This feature will also be useful for protein or ligand mutations in which two mutually exclusive sets of atoms (say, sidechains of differing identities) are built at the same site and need to be switched on and off via a global parameter.

[peastman]

I guess you would implement that like this: all the atoms in both sidechains would have their nonbonded parameters set to zero. Then you would add two additional custom forces. Each would implement the interaction between one side chain and everything else in the system. And each would scale its energy by a global parameter so you could turn it on and off. Is that the idea?

[jchodera]

@peastman That's exactly right! The two sidechains wouldn't interact via the CustomNonbondedGroupForce because the other "ghost sidechain" would not be in either atom group.

The global parameter scaling the energies of each sidechain version would also probably "soften" the interactions with softcore interactions as we currently do in ligand binding free energy calculations with CustomNonbondedForce.

[peastman]

What about the API? Should this be a new class, or just a new feature of CustomNonbondedForce? (I'm leaning toward the latter, since it's simpler and makes it easier to set up a system as usual and then restrict the computation.) Should NonbondedForce also have this feature? (I'm inclined to say no, since there's no way to make this work with those features of NonbondedForce that can't also be implemented with CustomNonbondedForce, such as PME.) How should we handle the truncation correction when using this feature? (I don't think there's really any meaningful way to define it, so my inclination is to throw an exception if the user tries to enable both this feature and the truncation correction at the same time.)

[jchodera]

It would be cleanest from an API perspective to add this as a feature of CustomNonbondedForce, where the user would only need an additional method call such as:

force.restrictToAtomGroups(group1, group2)

I think there is a meaningful way to implement an analytical truncation correction, but I will have to think about how that should work. For now, throwing an exception at Context creation time if the user tries to enable group restrictions and truncation corrections simultaneously is totally fine.

The cases I'm concerned about where the truncation correction is important are:

1. hydration free energies where the solute is in one group and the solvent in another group
2. protein sidechain mutations where the sidechain atoms are in one group and the solvent (and rest of protein) is in another group

Speaking of sidechain mutations, I'm playing around with ways to get the OpenMM app to set up a System that allows multiple residue choices at a given site, as discussed above. There are a few ways to implement this. I'll give my thoughts in a new issue soon, since I think these can be implemented most easily within (or alongside) ForceField. But if you had any thoughts, please share!

[mrshirts]

We've got some crazy ideas that would involve running hundreds of these, so ideally this would not require a huge amount of overhead to create or run. For example, we would have 100's of groups whose energies with the surrounding molecules would all be calculated.

[jchodera]

@mrshirts You also want the ability to retrieve the energy of each Force separately, correct?

[mrshirts]

Yes, correct, we need the energies of each Force group separately.

[mrshirts]

And forces, as well, if we want to linearly interpolate between groups.

[jchodera]

Can you clarify why you need the forces? If you put a linear scaling parameter lambdaXYZ into each energy expression, the forces will automagically be scaled by that amount and linearly combined in computing the total system force.

[jchodera]

I also wanted to note that this feature will also be very useful for relative free energy calculations in which some atoms are changed from dummy (noninteracting) atoms to real atoms, some atoms are changed from real to dummy, and some are changed from one real to another real type through softcore intermediates. This kind of relative transformation requires three different energy functions that depend on the same global parameter lambda, so implementing it with three separate group forces would be both simpler and more efficient than constructing a complex single CustomNonbonderForce with lots of per-particle parameters to pick out these different cases.

[peastman]

You will be able to create as many CustomNonbondedForce objects as you want. Note, though, that the number of force groups is limited to 32. That's because it uses bit flags to specify particular force groups.

A single CustomNonbondedForce will be able to specify many different particle groups. For NMR refinement, for example, you might have a hundred different groups, where each one involves just a few particles interacting with a few other particles.

@mrshirts Could you provide more information about what you have in mind?

[peastman]

I've created a reference implementation. Take a look at https://github.com/peastman/openmm/blob/9840ca70f11a857eaa90717f41509e4da4dfbb1e/openmmapi/include/openmm/CustomNonbondedForce.h. Does the API look ok?

[jchodera]

Looks good to me!

Are you requesting API use cases at this point, or are you just going to forge ahead?

[mrshirts]

Quick question: For truncation correction, are you referring to the analytical dispersion correction?

One possible solution: There is a contribution to the dispersion term from due to the pairwise terms included in each particular list of pairwise terms. If that quantity can be calculated on a per-pair-of-groups created, then the user can be responsible for figuring out how all the components of the dispersion correction can be added together to get the total system dispersion.

[jchodera]

@mrshirts Yes, but it's more general than a 'dispersion correction' because the energy function being computed may have nothing to do with dispersion forces. We call it a general analytical truncation correction.

[peastman]

Are you requesting API use cases at this point, or are you just going to forge ahead?

The most important thing at this point is to make sure the API is ok. I'm forging ahead with the implementation, so if you see any problems with the API, the sooner you tell me the better.

If that quantity can be calculated on a per-pair-of-groups created, then the user can be responsible for figuring out how all the components of the dispersion correction can be added together to get the total system dispersion.

The issue I see is that this feature can be used in so many different ways for so many different purposes. You might be using it just as an optimization, to avoid recalculating terms you know haven't changed. Or you might be using it as a decomposition, to look at particular components of the energy. Or you might be using it to implement restraints. You might be using it as the only nonbonded force, or in combination with a regular NonbondedForce or another CustomNonbondedForce. For any one of these use cases, you can say what is the "right" way to treat the truncation correction. But that won't be the same for all of them.

[jchodera]

The issue I see is that this feature can be used in so many different ways for so many different purposes. You might be using it just as an optimization, to avoid recalculating terms you know haven't changed. Or you might be using it as a decomposition, to look at particular components of the energy. Or you might be using it to implement restraints. You might be using it as the only nonbonded force, or in combination with a regular NonbondedForce or another CustomNonbondedForce. For any one of these use cases, you can say what is the "right" way to treat the truncation correction. But that won't be the same for all of them.

Starting with just letting the user switch the standard truncation correction 'on' (and having it default to 'off') would cover many cases, if not all. We definitely need the ability to switch the standard dispersion correction 'on' of alchemical free energy calculations, but I agree that it may make sense to have them 'off' so as not to surprise the user.

[peastman]

When you say "the standard truncation correction", do you mean the value computed based on all particles? Or one that is computed based only on the pairs in the interaction groups?

[jchodera]

When you say "the standard truncation correction", do you mean the value computed based on all particles? Or one that is computed based only on the pairs in the interaction groups?

Computed based on only pairs in the interaction group. We assume there is an isotopic distribution of particles in the box and compute the average interaction energy for the components beyond the cutoff, as usual.

Obviously this doesn't make sense for anything but periodic simulations.

@mrshirts does this sound right to you?

[mrshirts]

John, this is exactly what I was attempting to describe, but failing because I hadn't been following the correct terminology.

[Lnaden]

Based on how you extract energies in OpenMM, am I right in understanding if we are wanting only the energy of this CustomNonbondedForce with InteractionGroups we would have to: assign it to its own ForceGroup, get the State from the Context with the correct groups flags, then get the energy? This would mean that we would be limited by only 31 sets CustomNonbondedForces (leaving 1 for the rest of the forces).

If so, will there be a way to get the energy of a single InteractionGroup, or are we only able to access the CustomNonbondedForce with all of its InteractionGroups through the method above?

[peastman]

That is correct. No, there is no way to query the energy of just one interaction group out of several in a single CustomNonbondedForce.

Depending on what you're doing, it might be simplest to just repeatedly call setInteractionGroupParameters() then reinitialize the context. That would be fine for a one time energy evaluation, but obviously not for running a simulation. But if you're running a simulation, be aware that having 31 CustomNonbondedForces each with a single interaction group will be slower than a single CustomNonbondedForce with 31 interaction groups.

[Lnaden]

Since we cannot access the InteractionGroups directly, is it possible to identify which InteractionGroups and which set inside the InteractionGroup a particle belongs to? This would be possible to indirectly affect the evaluations by doing setParticleParameters() and updateParametersInContext(). An alternate method would be returning a list of particles of belonging to an InteractionGroup and which set they are in.

Is this the purpose of the getInteractionGroupParameters() function? Otherwise, I do not know what parameters it is looking at since it requires the index of the InteractionGroup and both sets of the atoms inside the InteractionGroup.

[peastman]

That's exactly what getInteractionGroupParameters() does. You give it the index of an interaction group, and it returns the two sets of particles.

[Lnaden]

This seems as though it would be cumbersome to decompose energies as it stands now.

For instance, if you have a protein with 10 side chains and you have 3 amino acids at each side chain you want test against all other amino acids, you would need to specify 10^3 groups to capture them all. You would then need just as many per-particle parameters to manipulate to actually decompose the energy (e.g. distinguishing from interactions between sites 1 and 2 vs sites 1 and 3) and seems to defeat the purpose of the InteractionGroup. Even if you were looking at a small subset, this quickly becomes unruly. I imagine there are other applications which suffer from the same problem. Is this how you see the energy decomposition or is there another way you could use the InteractionGroups to do this?

Would it be possible to send a command through the Force similar to updateParametersInContext() which would allow either turning InteractionGroups on/off so they are ignored or included in an energy calculation, or possibly adding/removing particles on the fly. This would make decomposing the energies easier and save on the bookkeeping from many per-particle parameters. I suspect the first would be easier since it avoids changing the number of particles in a the Force while its part of the Context.

[jchodera]

Would it be possible to send a command through the Force similar to updateParametersInContext() which would allow either turning InteractionGroups on/off so they are ignored or included in an energy calculation, or possibly adding/removing particles on the fly

I like this idea.

You can mimic some of the function of this by adding a globalParameter multiplying each force group interaction by a controllable scalar, but (1) there may be many such groups, and (2) setting the scalar multiple to zero doesn't save you any work. So there could be utility in turning these groups on or off.

[peastman]

Would it be possible to send a command through the Force similar to updateParametersInContext() which would allow either turning InteractionGroups on/off so they are ignored or included in an energy calculation, or possibly adding/removing particles on the fly.

Call setInteractionGroupParameters() to modify the particles in the interaction group, then call reinitialize() on the Context. Reinitializing the Context is more expensive than calling updateParametersInContext(), but that's because we carefully avoid doing anything in updateParametersInContext() that's too expensive. And modifying the particles in an interaction group requires rebuilding a data structure that's somewhat expensive to create.

[Lnaden]

Call setInteractionGroupParameters() to modify the particles in the interaction group, then call reinitialize() on the Context. Reinitializing the Context is more expensive than calling updateParametersInContext()

This is why I was wondering if it was possible just to efficiently ignore certain InteractionGroups when evaluating since you already have to keep track of them in some way. This would prevent having to modify the particles of a group, and instead just track the indices of the InteractionGroup.

[peastman]

Unfortunately, the implementation is more complicated than that. It breaks each interaction group into small pieces, then filters the pieces based on exclusions and duplicate interactions, and then merges them into larger blocks that can be computed efficiently on the GPU. So by the time it gets to the actual interaction kernel, all information about what interactions came from what groups has been lost.

[rmcgibbo]

This is implemented and now in the docs, so I think we can close this issue. Feel free to reopen if necessary.