Fast Python implementation of short-range interacting particles under Metroplis algorithm

Question

Can anyone write a Python implementation of a set of particles interacting in 2D according to a short-range particle-particle force and evolving in time under a Metropolis algorithm, which randomly chooses a particle and proposes a random spatial move? I can write such an implementation using the naive approach of updating all the particle-particle interaction energies associated with a given cell that has been proposed to be moved, but would like to see an implementation using a Verlet List or a KD Tree, which I hope can be at least an order of magnitude faster when the number of particles is large (say 100-1000).

An added complication in my particular problem is that the particles replicate and die, but incorporating that into the Verlet/KD Tree approach would be a secondary goal.

Thanks for any help you can offer.

Answer 1

Scipy offers an implementation of KDTree class:

http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.spatial.KDTree.html

Yu should be able to incorporate this into your code.

Answer 2

I have never done this myself, I do mostly MD. Still, your energy is

$E = \frac{1}{2} \sum_i e_i = \frac{1}{2} \sum_i \sum_{j\ne i} V_{ij}$

where the second sum is performed on neighbors only. If you store for every particle $i$ the sum $e_i = \sum_{j\ne i} V_{ij}$ you can just recompute it for the MC update. Explicitly, you compute the difference in energy for the particle that has moved. This idea extends to particle removal and addition easily.

The average cost is $O(N_\textrm{neigh})$ where $N_\textrm{neigh}$ is the average number of neighbors for a particle.