docs: def subtract_vectors(a, b, periodic):
      '''Returns the difference of the points vec_a - vec_b subject 
         to the periodic boundary conditions.
  
      '''
      r = a - b
      delta = np.abs(r)
      sign = np.sign(r)
      return np.where(delta > 0.5 * periodic, sign * (periodic - delta), r)