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