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
Calculates the X_{ij} matrix