Python API, CustomGBForce, TabulatedFunction, OpenCL vs. CPU

[Nicholas Schafer]

Hi Peter -

I'm working on a model that is implemented primarily using the Python API. One of my forces is a CustomGBForce that includes tabulated functions. When running a test simulation, simply changing the platform from OpenCL to CPU (or Reference) gives me different answers, and the cause of the problem seems to be the way the two platforms are reading the tabulated functions. I am using both Discrete2DFunctions and Discrete3DFunctions. If I eliminate the tabulated functions from the energy term expression, the platforms then give similar answers. When the tabulated functions are part of the energy expression, OpenCL is giving me an answer that is closer to what I expect. Can you think of what I might be doing wrong? Is there some particular order that I should call addPerParticleParameter, addGlobalParameter, addTabulatedFunction, addComptuedValue, and addEnergyTerm? What is the difference in the way that the platforms handle tabulated functions?

Thanks,
Nick

[Peter Eastman]

Is it possible for you to post an example that shows the problem? Just changing platforms shouldn't affect the results at all, aside from numerical precision. And tabulated functions should be interpreted exactly the same way on all platforms.

[Nicholas Schafer]

Hi Peter -

I've pared it down to a pretty minimal example. When I run testminimalexample.py, I get:

Code: Select all

CPU
#"Step","Potential Energy (kJ/mole)","Temperature (K)"
1,0.0,0.0
Reference
#"Step","Potential Energy (kJ/mole)","Temperature (K)"
1,0.0,0.0
OpenCL
#"Step","Potential Energy (kJ/mole)","Temperature (K)"
1,1.0,0.0

So, as I mentioned previously, OpenCL seems to be giving me the answer I expect, whereas both CPU and Reference are somehow missing one of the values in the tabulated function(?).

Thanks as always for the help,
Nick

[Peter Eastman]

Sorry for taking a while to get to this! I've been a bit overwhelmed with other things.

I believe the problem is that your force is nonsymmetric. Here's what the API documentation says about energy terms of type ParticlePair:

Note that energy terms are assumed to be symmetric with respect to the two interacting particles, and therefore are evaluated only once per pair. In contrast, expressions for computed values need not be symmetric and therefore are calculated twice for each pair: once when calculating the value for the first particle, and again when calculating the value for the second particle.

Your system has 368 atoms, and you define the tabulated function like this:

Code: Select all

gamma_ij = [0.0]*self.natoms*self.natoms
gamma_ij[32] = 1.0
am_dd.addTabulatedFunction("gamma_ij", Discrete2DFunction(self.natoms, self.natoms, gamma_ij))

So gamma_ij(32, 0) is 1, but gamma_ij(0, 32) (along with all other values) is 0. Here's how you define your energy term:

Code: Select all

am_dd.addEnergyTerm("gamma_ij(index1,index2)", CustomGBForce.ParticlePair)

It isn't symmetric, which means the result is undefined. The interaction between atoms 0 and 32 will only be calculated once, and you can't depend on whether it will end up equal to gamma_ij(0, 32) or gamma_ij(32, 0). Different platforms happen to pick different ones.

[Nicholas Schafer]

Hi Peter -

No trouble at all. Thanks for looking into it. I made everything symmetric and now all of the platforms are giving nearly identical answers. Good catch!

Thanks,
Nick