Implementation of nonbondend interaction force between two groups for enhanced sampling algorithm

[Fabian Paul]

Dear Peter, dear community,

I'm thinking of porting an enhanced sampling scheme, that has been successfully demonstrated for protein-peptide binding [1], to OpenMM. I would like to ask for advice.

The enhanced sampling scheme can be best described a "Replica exchange between differently boosted Lenneard-Jones protein-peptide interaction Hamiltonians". Originally, I had implemented it as a plugin for ACEMD [2].

I would like to have full open-source version of the algorithm. That's why I'm considering OpenMM as a platform.

There are a couple of problems that would have to be solved for an implementation:

(a) The nonbonded force that is biased is between two molecules (or more generally between two selections of atoms). Usually these selections do not involve the solvent. So a very efficient implementation should be possible. (In fact the original implementation of the scheme runs as a CPU plugin to ACEMD without slowing down ACEMD much).
(b) The boosting potential requires two "passes" over the atoms. A boost potential can be seen as an adaptive scaling factor that is applied to the force F (chain rule!). The adaptive scaling factor itself is a function of the energy that corresponds to the unbiased F. So one pass over the atoms is needed to compute the unbiased energy and a second pass over the atoms is needed to compute the forces.
(c) Replica exchange requires that the hyper-parameters can be exchanged quickly between parallel simulations.

I would already be very happy, if (a) can be solved. What I would like to achieve is to have some degree of generality in the implementation of (a). Not only is it conceivable to apply such an algorithm to the Lennard-Jones interaction, but also to the direct space short-range part of the Coulomb interaction. Or both of them. Also switching/shifting functions should be respected, if the force field uses them. I would like to have a platform for some experimentation, where I can play with different nonbonded forces.

To implement (b), is it perhaps possible to use the aMD integrator? However let me add that I'm not interested to use boost potentials to compute kinetics with aMD. What I like instead about boost potentials is that they leave the potential unbiased if it falls below a given threshold. This increases thermodynamic overlap and produces Boltzmann reweighting factors that are = 1 for many conformations. I presume it is possible to implement the two passes inside the same force object.

Probably the implementation would require to write an custom CUDA kernel. I would be fine to do that. Could you provide me with some intial pointers to the OpenMM internals documentation? Could you give me some hints how to implement the atom selection and an efficient loop over the two selections, which also uses distance information like a pair list?

To re-iterate: just solving (a) and leaving out the more difficult problems for later would be very nice.

Thanks,
Fabian

[1] https://doi.org/10.1038/s41467-017-01163-6

[2] Suppl. Note 3.1
https://static-content.springer.com/esm ... M1_ESM.pdf

[Peter Eastman]

The first thing I'd suggest is looking at the implementations of the various aMD integrators: https://github.com/pandegroup/openmm/bl ... nmm/amd.py. They're all based on CustomIntegrator so they're very simple, just a few lines of Python code. There's a good chance you can implement the boosting for your own algorithm in the same way.

For your nonbonded force between molecules, you may be able to implement it with CustomNonbondedForce. You can define "interaction groups" to restrict it to a small subset of interactions. See http://docs.openmm.org/latest/api-pytho ... ondedForce.

If you decide you do need to write a plugin, the place to start is the developer guide: http://docs.openmm.org/latest/developerguide/index.html. It describes the sort of internal information you're looking for.