code: def get_delta_plus_epsilon(nat, i_now, dr_now, dr_old):
      
      XIJ = get_XIJ(nat, i_now, dr_now, dr_old)
      YIJ = get_YIJ(nat, i_now, dr_old)
  
      YIJ_invert = array_inverse(YIJ)
  
      # Perform sum_k X_ik Y_jk^-1
      epsilon = np.sum(XIJ.reshape(-1,3,1,3)*YIJ_invert.reshape(-1,1,3,3), axis=3)
  
      return epsilon
docs: Calculate delta_ij+epsilon_ij, i.e. the deformation gradient matrix