KDTree For Chemical Monte Carlo

The Problem:

This was one of those projects you start out thinking “how hard could it be?” and then after some hacking, the fine print comes, finds you, and slaps you up side the head.  So, what were we trying to accomplish anyway?  Well, an old high school buddy of mine works on these chemical simulations that at their core are pretty approachable:

1. Get a pizza box.
2. Cut a hole in the box.
3. Put N particles in the box interacting with only pair-wise central potentials (if its easier, you can imagine our particles are all planets in outer-space only feeling the force of gravity on each other).
4. Allow system to evolve, in the Monte Carlo sense, at finite temperature.  This means sometimes two particles that attract each other move towards each other (minimizing potential energy) and sometimes they move away.  Damn you, heat.
5. Submit like a bazillion jobs, and wait patiently as your thesis project is delayed by slow code.

Slow code, did you say?? Well slow is relative, so why don’t we say that there was a bottleneck inside the thermal simulation.  The problem is was that the pair-wise potentials ma buddy wanted to study are short-ranged but the interparticle distance was calculated for all pairs.  Or rather, that’s the *advantage* we exploited.  I could babble on and on about the science details but lemme cut to the chase.  If you have a system of N particles, you only need to examine those in some small 3D range at any given step.  We don’t need to visit each N particle and ask “Are you close enough to my test particle, that I *actually* need to calculate the force???”.  The hope (and end result) is that we can use some non-trivial geometrically-aware data structure to store our particles that allows us to query it faster than O(N).   And finally…we can’t make any approximations.  We’re here for the science, god damnit!

The Solution:

Turns out, such a data structure exists (no surprise); I implemented a KD-Tree with pretty significant modifications.  Before you, oh avid read, proclaim “dude you said range.  why aren’t you using a range-tree?”  I will shameless admit that I made this code to solve a related but simplified pair-wise interaction that only needed the nearest-neighbor, so at the time KD-Tree looked perfect.  Also, people use OctTrees for these things too, and in theory I could cross-compare, but our speed-up was so significant that we became dominated by other things.  Anyways, if you don’t know much about KD-Trees think of them as the k-dimensional generalization of balanced sort tree.  Also read more here . And check out one for 210 particles in only 2-d…

But I said ours has pretty significant modifications, like what, oh master of big claims?

Handling wrap-around

This was, hands-down, the hardest issue to solve because it breaks the fundamental assumption of KD-Trees.  Yea, that’s right, I hammered on that square peg so damn hard, I got it through the circular hole.  I have no regrets.  The basic problem is that our universe is like a circle not a line: go too far in one direction and you get back to where you started.  Ru-roh.  I’m gonna go more in depth on this wily topic, with smart jargon, simple explanations and a solution for k-dimensions, in another post!

Mutability

So in the standard texts KD-Trees are often pitched as a build-once query many times.  But ours is different, we build it once, query it once, move only one particle, then query it once again, the move that one particle and so on.  I stared at this problem for a really long time, I was sure I could do some kind of smart balancing (a la AVL or Red-Black trees).  In the end, to move a particle I just remove it from the tree, rebalance all the nodes under it, then reinsert it.  Inelegant, I agree, but it might be the theoretical best you can do.  I wish I had a proof to offer you =\

Range Searches

KD-Trees are really optimal for nearest-neighbor searches, but can also be used to find the m nearest neighbors.  That’s not too hard if you push them nodes into a priority queue as you traverse the tree.  Push Push Push; Push Pop Push Pop–we whistle while we work.  If you want to collect the m you just don’t start popping nodes till you have at least m in the queue.  To do range searches, you don’t start popping until the furthest node is greater than your search range.  This is really more about the queue object/search algorithm than the tree itself.

Indexing

Super boring, but I guess I should list it.  Our data structure maintains the ability to perform O(1) index searches, so that finding the third particle you inserted, always takes the same amount (ie 0) time.  This isn’t too tricky, but was crucial for science reasons.  Science.

Wow, you read a page of text. Congrats for getting this far and thanks for reading!  I’ll try to get up