def get_delta_plus_epsilon(nat, i_now, dr_now, dr_old):
XIJ = get_XIJ(nat, i_now, dr_now, dr_old)<EOL>YIJ = get_YIJ(nat, i_now, dr_old)<EOL>YIJ_invert = array_inverse(YIJ)<EOL>epsilon = np.sum(XIJ.reshape(-<NUM_LIT:1>,<NUM_LIT:3>,<NUM_LIT:1>,<NUM_LIT:3>)*YIJ_invert.reshape(-<NUM_LIT:1>,<NUM_LIT:1>,<NUM_LIT:3>,<NUM_LIT:3>), axis=<NUM_LIT:3>)<EOL>return epsilon<EOL>
Calculate delta_ij+epsilon_ij, i.e. the deformation gradient matrix