def METHOD_NAME(nat, i_now, dr_old):
    """
    Calculates the Y_{ij} matrix
    """
    # 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