def get_delta_plus_epsilon(nat, i_now, dr_now, dr_old):
    """
    Calculate delta_ij+epsilon_ij, i.e. the deformation gradient matrix
    """
    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
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 ) epsilon = np . sum ( XIJ . reshape ( - 1 , 3 , 1 , 3 ) * YIJ_invert . reshape ( - 1 , 1 , 3 , 3 ) , axis = 3 ) return epsilon
Calculate delta_ij+epsilon_ij, i.e. the deformation gradient matrix