docs: def get_XIJ(nat, i_now, dr_now, dr_old):
      """
      Calculates the X_{ij} matrix
      """
      # Do an element-wise outer product
      dr_dr = dr_now.reshape(-1,3,1)*dr_old.reshape(-1,1,3)
  
      xij = np.zeros([nat,3,3])
      for i in range(3):
          for j in range(3):
              # For each atom, sum over all neighbors
              xij[:,i,j] = np.bincount(i_now, weights=dr_dr[:,i,j])
  
      return xij