docs: 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