atom reordering with cutoffs

The functionality of OpenMM will (eventually) include everything that one would need to run modern molecular simulation.
POST REPLY
User avatar
Charles Brooks
Posts: 35
Joined: Fri Feb 24, 2012 11:48 am

atom reordering with cutoffs

Post by Charles Brooks » Wed Feb 03, 2016 3:15 pm

Peter et al.,

We are trying to finish up our CPHMD plugin to complement the recent implementation of the GBSW model in the CHARMM/OpenMM interface and are running into an issue that I think I ran across before but cannot find any detail of. This has to do with reordering of atoms to optimize the tile-based computations on the GPU (in cuda) when one has cutoffs. I am pretty sure this happens. We discovered an array (atomIndex) and a cuda platform kernel routine (CudaContext.cpp) where we believe this happens, our problem is 1) we want to gather the energies of interaction between specific sets of atoms and their neighbors at each timestep and 2) we use these energies to propagate the values of the charges on the atoms, and hence need to repopulate the charges with their updated values after the propagation step. It appears that this array, atomIndex, contains the information we need, but applying it to "reorder" our calculation doesn't seem to yield the correct result. I note that if we do no cutoff and don't use this array to reorder atoms, our results agree with those from or reference or CPU platform (i.e., CHARMM). Clearly we are not doing something correct, but we don't seem to be able to figure out where/why. Any help would be appreciated.

Thanks,

Charlie

User avatar
Peter Eastman
Posts: 2543
Joined: Thu Aug 09, 2007 1:25 pm

Re: atom reordering with cutoffs

Post by Peter Eastman » Wed Feb 03, 2016 3:46 pm

Hi Charlie,

Two possible solutions, depending on exactly what you're doing.

First, you may simply be applying the atomIndex array in the wrong way. If atomIndex[ i] = j, that means the i'th atom in the internal ordering (as seen by GPU kernels) corresponds to atom j in the original ordering (as seen by public API methods). Whenever I get confused by this (which I do, often!), I look at the implementation of get/set routines for positions or velocities to make sure of the correct way to do the indexing: https://github.com/pandegroup/openmm/bl ... #L234-L291

The other possibility is simply to prevent any reordering that affects your force. We mostly try to do that, so we don't have to update lots of arrays of force field parameters after every reordering. The only type of reordering that can happen is to swap the indices of identical molecules. Mostly that means we reorder the waters. But since all the force field parameters are the same for every water molecule, that doesn't affect them.

The decision about what is "identical" is determined by the CudaForceInfo classes that all the forces register. They're described in more detail at http://docs.openmm.org/6.3.0/developerg ... -particles. So possibly you just need to define and register a CudaForceInfo. From your description, though, I'm not sure that's actually the right solution for your situation.

Peter

POST REPLY