[SciPy-User] Distance matrix with periodic boundary conditions

[Steve Schmerler]

On Nov 17 11:17 -0600, Dan Lussier wrote:
> Hi Ramon,
> 
> I've run into a similar experience before - I wasn't able to find
> another scipy module that did the job so I ended up doing a somewhat
> crude implementation of my own for doing post-processing of molecular
> dynamics simulation data.
> 
> It is a variation on a cell-linked-list decomposition method that is
> commonly seen in molecular simulation and relies on the particle
> interaction being cutoff at some small or intermediate value relative
> to the total size of the domain.

I came across the same problem while calculating the radial pair
distribution function (RPDF) from molecular dynamics data. If you do not
have too many atoms, you might get away with a simple numpy approach
(which will require somewhat more memory though.) 

Say you have a MD trajectory `coords` as a 3d array of shape (N, 3, T),
where N = number of atoms, T = number of time steps. Instead of a double
loop:

    sij = coords[:,None,...] - coords[None,...] 

This gives you an array (N, N, 3, T) of distance vectors for each time
step. I always have to consult [1] for this kind of stuff. 

Introducing PBC means applying the Minimum Image Convention [2] to get
nearest neighbor separations.

    sij[sij > 0.5] -= 1.0
    sij[sij < -0.5] += 1.0

Note that this assumes that you have coords in fractional (reduced,
crystal) coordinates.

To get the distances, transform to cartesian coords. For an arbitrarily
shaped simulation cell with the cell basis vectors as rows of a 3x3
array `cell`

    sij = sij.reshape(N**2, 3, T)
    rij = np.dot(sij.swapaxes(-1,-2), cell).swapaxes(-1,-2)
    dists = np.sqrt((rij**2.0).sum(axis=1))

Note that the sketched algorithm double-counts each separation
(sij[i,j,...] == -sij[j,i,...]), which seems not to be very clever.
However, one usually deals with different atom selections (e.g. the RPDF
between two different types of atoms, sij = coords1[..] - coords2[..].
Then, sij = sij.reshape(N1*N2, 3, nstep). For only one atom type, one
could select the upper (or lower) "triangle": ind =
np.triu_indices(N,k=1); sij=sij[ind[0], ind[1],...].

[1] http://www.scipy.org/EricsBroadcastingDoc
[2] M. P. Allen, D. J. Tildesley, Computer Simulation of Liquids

best,
Steve