Use OpenMM to calculate neighbors list?

[Maxim Imakaev]

Dear Peter,

I'm working with polymer systems, and I often need to calculate contact maps for very large polymer systems (up to a million monomers) for thousands of conformations. I'm currently doing it using an in-house C++ code which works in about O(N^(1.3-1.6)).

I was wondering if it is possible to use OpenMM to get a contact list for a given conformation (i.e. list of particle pairs separated by less than cutoff distance)? I know Openmm does this to calculate NonBonded force, but I can't see a way to extract this information out of OpenMM for now. It would significantly speed up my code, as now it takes tens of seconds for a million particles, and OpenMM perhaps could do it in under a second.

Thanks,
Max

[Yutong Zhao]

Unfortunately this information isn't exposed directly. In fact, we don't even build a canonical neighborlist on the GPU. We find neighbors for a block of atoms instead of a single atom.

The relevant code are located here:

https://github.com/SimTk/openmm/blob/ma ... gBlocks.cu
https://github.com/SimTk/openmm/blob/ma ... lities.cpp

Note that these functions operate on the reordered atoms, since once every 100 steps we reorder the atoms using a Hilbert curve to improve spatial coherence and memory access.

Maybe it's possible to just do your entire simulation in OpenMM with customForces if you'd like to go into a bit more detail.

[Peter Eastman]

We do have a reference implementation, though, which you could call. ReferenceNeighborList::computeNeighborListVoxelHash() computes a neighbor list. It's not very fast, but it does use an O(n) algorithm, so it may be faster than what you're doing right now. In the latest version of the code in Github, we also have CpuNeighborList which is much faster. Perhaps one of those would suite your needs?

Peter

[Maxim Imakaev]

Thanks, Peter!

I'll try these and see how it goes.

M

[Maxim Imakaev]

Hi Peter,

I was able to strip down testNeighborsList and get a c++ code which calculates neighbors list. It is indeed very fast! I will try to incorporate it in our library. This will speed up our contact map calculations by a factor of 100 or so!

I also noted that it returns more particle pairs that are actually in contact (e.g. for 256000 particles it returned 3970540 pairs, but only 2123960 of them were actually less than the cutoff distance apart). This list however included all the actual contacts. Most of the "extra" pairs were at a slightly longer distances.

I'm assuming that this is a consequence of a particular implementation, and I can filter those pairs with a line of code, but I thought I'd let you know about it.

Thank you!
Max

[Maxim Imakaev]

Dear Peter,

I'm trying to wrap this function in Python, but I don't even know how to start. My knowledge of C++ is limited to editing code and writing simple C programs/inserts (here, I just modified the test file and re-run make... and it somehow worked).

Is there a way to access this method now using current Python bindings? If not, is it possible to modify SWIG files and make it accessible?

Thanks,
Max

[Peter Eastman]

Hi Max,

I also noted that it returns more particle pairs that are actually in contact

That's because it operates on blocks of four atoms: it picks four atoms that are close together, then lists other atoms that are within the cutoff distance of any one of those four atoms. It works that way because the code to compute nonbonded interactions is vectorized, so it computes four different interactions at a time.

Is there a way to access this method now using current Python bindings?

No, the Python bindings only expose the public API, not the internal details of each platform. If you want an easy way to invoke the C++ code directly from Python, take a look at Cython or Weave.

Peter

[Maxim Imakaev]

Hi Peter,

I have been using OpenMM to find neighbor pairs in 3D for a couple years now. For large systems (10k+), OpenMM works faster than any other solution I found, including the KD tree in scipy and other methods. This is all despite the fact that piping the data between C++ and Python in my code takes more time than the actual neighbor list calculation.

As for seamlessly integrating it with Python, I found that the code needs to be compiled against the current version of OpenMM, and my library ended up with three different versions of binaries. I also never figured out the way to pass the proper compile arguments to Cython, so I ended up compiling a C++ binary (by replacing TestCpuNeighborsList with my own code) and piping binary data in and out of it.

Although it works for now, I found that sharing it with the community would still be a problem. For example, in my recent stackOverflow post, I could not provide any ready-to-use recipe to use OpenMM to calculate neighbors. http://stackoverflow.com/questions/2128 ... 3#34257083 (note that the test there is for a uniform arrangement of points; for non-uniform systems with high-density regions, OpenMM has a noticeable advantage).

Do you think it would be a nice feature to include in the OpenMM Python bindings? It would then be version-independent, and would not require the user to have GCC installed. How difficult would it be to do? I could try to write it myself and commit it to github, but I'm not sure if I have enough skills to do that.

Thanks,
Max and the Mirnylab.

[Peter Eastman]

Hi Max,

Funny you should ask about this! I just started writing a neighbor list builder for mdtraj (http://mdtraj.org/latest). It's going to be simpler and less optimized than what's in the OpenMM CPU platform. At least to start with, I'm just using Numpy and Cython. But it will still use an O(N) voxel based algorithm, so hopefully it will be fast enough for your purposes. Once I have it written, I can point you to the code to try it out.

Peter

[Maxim Imakaev]

That would be great! I would definitely incorporate it in our lab pipeline.

Also, starting with Scipy 0.17, their cKDTree class has a nice contact finder, cKDTree.query_pairs(cutoff,output_type="ndarray"). It is faster than the previous version because they now output contacts as a numpy array and not as a set of tuples. That contact finder may be hard to beat.

Max