# Subtraction and Periodic Boundary Conditions

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

$d^{2} = \sum_i \left(\vec{r}_{i} \ominus_{\vec{L}} \vec{s}_{i} - \vec{\delta}\right)^2$

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

$d^2 = \frac{1}{N}\sum_{\left\langle i,j\right\rangle }\left(\vec{r}_{ij} \ominus_{\vec{L}} \vec{s}_{ij}\right)^2$

The reason this math works (and the $$\vec \delta$$ drops out!) is talked about a bit more in this post. The Python function for this is as follows: