If you have two points, and periodic boundary conditions, and you want to know the closest distance vector between them, that’s easy to write in python:

If you have a set of, say, 10 points, and you want to know the distance between all 55 pairs of them, that’s still not too complicated:

Note that the above function is really only two lines, yet it has some nice properties:

It takes vectors of any length d

Its vectorized, aside from the splitting by dimension (the list comprehension)

The box size can be a single float, or a vector for rectangular-but-not-square boxes

In Matlab

We can do something similar in Matlab:

The matlab code has about 8 lines to Python’s 2, and also can’t handle non-square boxes, but it is also using vectorized operations and any number of dimensions.

Benchmarking

Out of curiosity, I benchmarked these two functions:

Input size (rs)

Matlab

Python

100 x 3

0.7 ms

2.6 ms

200 x 3

2.5 ms

8.3 ms

500 x 3

18.2 ms

60.8 ms

1000 x 3

83.3 ms

288.0 ms

Matlab clearly has the upper hand, by a factor of ~3.

Python Update

We can also do a nice broadcasting method:

Note that we had to add in a call to np.rollaxis here, which converts diffs
from \(N \times d \times N\) to \(N \times N \times d\), and allows us to broadcast
over the adding L in the case where L is not a scalar.

We can also extend this to a pairwise function for comparing two different sets
of particles:

And if our two pairs of particles may have had a center-of-mass motion that we
want to ignore, we can calculate

where \(\ominus_{\vec{L}}\) means “shortest distance given periodic boundary
conditions in a box of shape \(\vec{L}\),” and \(\vec{\delta}\) is the center-of-mass
motion. It turns out that this is equal to