def get_YIJ(nat, i_now, dr_old):
<EOL>dr_dr = dr_old.reshape(-<NUM_LIT:1>,<NUM_LIT:3>,<NUM_LIT:1>)*dr_old.reshape(-<NUM_LIT:1>,<NUM_LIT:1>,<NUM_LIT:3>)<EOL>yij = np.zeros([nat,<NUM_LIT:3>,<NUM_LIT:3>])<EOL>for i in range(<NUM_LIT:3>):<EOL><INDENT>for j in range(<NUM_LIT:3>):<EOL><INDENT>yij[:,i,j] = np.bincount(i_now, weights=dr_dr[:,i,j])<EOL><DEDENT><DEDENT>return yij<EOL>
Calculates the Y_{ij} matrix