search

choderalab / openmmtools

A batteries-included toolkit for the GPU-accelerated OpenMM molecular simulation engine.
http://openmmtools.readthedocs.io
MIT License
250 stars 80 forks source link

Alchemical CustomGBForce Issues #287

Open Lnaden opened 7 years ago

Lnaden commented 7 years ago

There are some issues I have been trying to track down in choderalab/yank#768 and in turn have found 2 issues with the Alchemical Custom GB force.

The first is that there is a term missing from the calculation of the I term here which differs from the energy expression in the OpenMMDocs for a non-alchemical CustomGBForce here. Specifically, the C term is missing from the first line of the computation so it should actually look like this:

... log(L/U)/r+C);"
    "U=r+sr2;"
    "C=2*(1/or1-1/L)*step(sr2-r-or1);"
    ...

If you compare what is currently in the code, with and without the C term relative to the built in GBSAOBCForce with the same particles, you get something that looks like this:

Force Comparison

So although this only affects particles very close to each other, it is still incorrect and with alchemical particles, they can get very close

Secondly, this image shows something else odd which happens, when the particles are 0 distance apart, the CustomGBForce energy and force go to NaN, but the GBSAOBCForce does not, which may account for what I see in choderalab/yank#768, but may also be a red herring in this case.

While I think we should fix the first issue, the second one may need some discussion.

jchodera commented 7 years ago

We definitely need to fix the NaN issue, and it sounds like it is easy to fix the missing terms issue.

Should we add a switching function that sends the force to zero at zero distance?

Lnaden commented 7 years ago

easy to fix the missing terms issue.

That should be easy enough to fix

Should we add a switching function that sends the force to zero at zero distance?

There are a couple issues with this. The first is the force is not actually zero at r=0. Second, I'm not sure a switch would help given how I think Custom*Forces handles a NaNs such that something + 0*NaN is NaN. Something to test. @andrrizzi has also reported that he does not see the NaN issues in implicit solvents that I am seeing with some of the FXR data set, so solving the NaN issue here may not actually fix the NaN problem in the FXR set. I just happened to find this while I was looking into the CustomGBForce