Return the squared distance between all combinations of
    a set of points, in periodic space.

    Parameters
    ----------
    r: shape (n, d) for n points in d dimensions.
        Set of points
    L: float array, shape (d,)
        System lengths.

    Returns
    -------
    d_sq: float array, shape (n, n, d)
        Squared distances

def pdist_sq_periodic(r, L):
    """Return the squared distance between all combinations of
    a set of points, in periodic space.

    Parameters
    ----------
    r: shape (n, d) for n points in d dimensions.
        Set of points
    L: float array, shape (d,)
        System lengths.

    Returns
    -------
    d_sq: float array, shape (n, n, d)
        Squared distances
    """
    d = csep_periodic(r, r, L)
    d[np.identity(len(r), dtype=np.bool)] = np.inf
    d_sq = np.sum(np.square(d), axis=-1)
    return d_sq