Skip to content

Commit

Permalink
guarantee that total energy is unchanged when not updating exchange-c…
Browse files Browse the repository at this point in the history 
…orrelation energy
  • Loading branch information
@JWilhelm @fstein93
JWilhelm authored and fstein93 committed Jun 29, 2020
1 parent 84d4f61 commit 4a622d8
Showing 1 changed file with 63 additions and 24 deletions.
87 changes: 63 additions & 24 deletions src/mp2.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -470,6 +470,7 @@ SUBROUTINE mp2_main(qs_env, calc_forces)
IF (mp2_env%ri_rpa%do_ri_g0w0) THEN
CALL compute_vec_Sigma_x_minus_vxc_gw(qs_env, mp2_env, mos_mp2, E_ex_from_GW, t3, unit_nr)
END IF
IF (free_hfx_buffer) THEN
CALL timeset(routineN//"_free_hfx", handle2)
CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
Expand Down Expand Up @@ -689,6 +690,7 @@ SUBROUTINE mp2_main(qs_env, calc_forces)
CASE DEFAULT
CPABORT("")
END SELECT
t2 = m_walltime()
IF (unit_nr > 0) WRITE (unit_nr, *)
IF (mp2_env%method .NE. ri_rpa_method_gpw) THEN
Expand Down Expand Up @@ -723,9 +725,9 @@ SUBROUTINE mp2_main(qs_env, calc_forces)
ELSE
IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.6)') 'Total RI-RPA Time=', t2 - t1
update_xc_energy = .TRUE.
IF(mp2_env%ri_rpa%do_ri_g0w0 .AND. .NOT. mp2_env%ri_g0w0%update_xc_energy) update_xc_energy = .FALSE.
IF(.NOT. update_xc_energy) Emp2 = 0.0_dp
update_xc_energy = .TRUE.
IF (mp2_env%ri_rpa%do_ri_g0w0 .AND. .NOT. mp2_env%ri_g0w0%update_xc_energy) update_xc_energy = .FALSE.
IF (.NOT. update_xc_energy) Emp2 = 0.0_dp
IF (unit_nr > 0 .AND. update_xc_energy) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'RI-RPA energy = ', Emp2
IF (mp2_env%ri_rpa%do_ri_axk) THEN
Expand Down Expand Up @@ -1206,10 +1208,14 @@ SUBROUTINE calculate_exx(qs_env, unit_nr, do_gw, do_admm, E_ex_from_GW, t3)
! take the exact exchange energy from GW or calculate it
IF (do_gw) THEN
IF(qs_env%mp2_env%ri_g0w0%update_xc_energy) THEN
IF (qs_env%mp2_env%ri_g0w0%update_xc_energy) THEN
CALL remove_exc_energy(energy)
energy%total = energy%total + E_ex_from_GW
energy%ex = E_ex_from_GW
t2 = m_walltime()
IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.6)') 'Total EXX Time=', t2 - t1 + t3
IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'EXX energy = ', energy%ex
END IF
ELSE
Expand All @@ -1236,12 +1242,13 @@ SUBROUTINE calculate_exx(qs_env, unit_nr, do_gw, do_admm, E_ex_from_GW, t3)
! include the EXX contribution to the total energy
energy%ex = ehfx
energy%total = energy%total + energy%ex
END IF
t2 = m_walltime()
t2 = m_walltime()
IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.6)') 'Total EXX Time=', t2 - t1 + t3
IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'EXX energy = ', energy%ex
IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.6)') 'Total EXX Time=', t2 - t1 + t3
IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'EXX energy = ', energy%ex
END IF
! reset to zero the Hartree-Fock energy
energy%ex = 0.0_dp
Expand All @@ -1250,25 +1257,31 @@ SUBROUTINE calculate_exx(qs_env, unit_nr, do_gw, do_admm, E_ex_from_GW, t3)
END SUBROUTINE calculate_exx
! **************************************************************************************************
!> \brief ...
!> \param energy ...
! **************************************************************************************************
SUBROUTINE remove_exc_energy(energy)
TYPE(qs_energy_type), POINTER :: energy
CHARACTER(len=*), PARAMETER :: routineN = 'remove_exc_energy', routineP = moduleN//':'//routineN
CHARACTER(len=*), PARAMETER :: routineN = 'remove_exc_energy', &
routineP = moduleN//':'//routineN
INTEGER :: handle
CALL timeset(routineN, handle)
! Remove the Exchange-correlation energy contributions from the total energy
energy%total = energy%total - (energy%exc + energy%exc1 + energy%ex + &
energy%exc_aux_fit + energy%exc1_aux_fit)
! Remove the Exchange-correlation energy contributions from the total energy
energy%total = energy%total - (energy%exc + energy%exc1 + energy%ex + &
energy%exc_aux_fit + energy%exc1_aux_fit)
energy%exc = 0.0_dp
energy%exc1 = 0.0_dp
energy%exc_aux_fit = 0.0_dp
energy%exc1_aux_fit = 0.0_dp
energy%ex = 0.0_dp
energy%exc = 0.0_dp
energy%exc1 = 0.0_dp
energy%exc_aux_fit = 0.0_dp
energy%exc1_aux_fit = 0.0_dp
energy%ex = 0.0_dp
CALL timestop(handle)
CALL timestop(handle)
END SUBROUTINE
Expand Down Expand Up @@ -1302,8 +1315,8 @@ SUBROUTINE compute_vec_Sigma_x_minus_vxc_gw(qs_env, mp2_env, mos_mp2, energy_ex,
LOGICAL :: charge_constrain_tmp, do_admm_rpa, &
do_hfx, do_kpoints_cubic_RPA, &
do_ri_Sigma_x, really_read_line
REAL(KIND=dp) :: eh1, ehfx, exx_minus_vxc, hfx_fraction, &
t1, t2, tmp
REAL(KIND=dp) :: eh1, ehfx, energy_exc, energy_exc1, energy_exc1_aux_fit, &
energy_exc_aux_fit, energy_total, exx_minus_vxc, hfx_fraction, t1, t2, tmp
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: vec_Sigma_x_minus_vxc_gw, &
vec_Sigma_x_minus_vxc_gw_im
TYPE(admm_type), POINTER :: admm_env
Expand All @@ -1316,6 +1329,7 @@ SUBROUTINE compute_vec_Sigma_x_minus_vxc_gw(qs_env, mp2_env, mos_mp2, energy_ex,
TYPE(dbcsr_type) :: matrix_tmp, matrix_tmp_2, mo_coeff_b
TYPE(dft_control_type), POINTER :: dft_control
TYPE(kpoint_type), POINTER :: kpoints
TYPE(qs_energy_type), POINTER :: energy
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(qs_rho_type), POINTER :: rho
TYPE(section_vals_type), POINTER :: hfx_sections, input, xc_section, &
Expand Down Expand Up @@ -1349,7 +1363,8 @@ SUBROUTINE compute_vec_Sigma_x_minus_vxc_gw(qs_env, mp2_env, mos_mp2, energy_ex,
input=input, &
dft_control=dft_control, &
para_env=para_env, &
ks_env=ks_env)
ks_env=ks_env, &
energy=energy)
! RPA/GW with ADMM for EXX or the exchange self-energy only implemented for
! ADMM_PURIFICATION_METHOD NONE
Expand All @@ -1369,7 +1384,8 @@ SUBROUTINE compute_vec_Sigma_x_minus_vxc_gw(qs_env, mp2_env, mos_mp2, energy_ex,
dft_control=dft_control, &
para_env=para_env, &
kpoints=kpoints, &
ks_env=ks_env)
ks_env=ks_env, &
energy=energy)
nkp = kpoints%nkp
Expand All @@ -1381,7 +1397,8 @@ SUBROUTINE compute_vec_Sigma_x_minus_vxc_gw(qs_env, mp2_env, mos_mp2, energy_ex,
input=input, &
dft_control=dft_control, &
para_env=para_env, &
ks_env=ks_env)
ks_env=ks_env, &
energy=energy)
nkp = 1
Expand Down Expand Up @@ -1507,10 +1524,32 @@ SUBROUTINE compute_vec_Sigma_x_minus_vxc_gw(qs_env, mp2_env, mos_mp2, energy_ex,
END DO
END IF
IF (.NOT. mp2_env%ri_g0w0%update_xc_energy) THEN
energy_total = energy%total
energy_exc = energy%exc
energy_exc1 = energy%exc1
energy_exc_aux_fit = energy%ex
energy_exc1_aux_fit = energy%exc_aux_fit
energy_ex = energy%exc1_aux_fit
END IF
! Remove the Exchange-correlation energy contributions from the total energy
energy%total = energy%total - (energy%exc + energy%exc1 + energy%ex + &
energy%exc_aux_fit + energy%exc1_aux_fit)
! calculate KS-matrix without XC and without HF
CALL qs_ks_build_kohn_sham_matrix(qs_env=qs_env, calculate_forces=.FALSE., &
just_energy=.FALSE.)
IF (.NOT. mp2_env%ri_g0w0%update_xc_energy) THEN
energy%exc = energy_exc
energy%exc1 = energy_exc1
energy%exc_aux_fit = energy_ex
energy%exc1_aux_fit = energy_exc_aux_fit
energy%ex = energy_exc1_aux_fit
energy%total = energy_total
END IF
! set the DFT functional and HF fraction back
CALL section_vals_val_set(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", &
i_val=myfun)
Expand Down

0 comments on commit 4a622d8

Please sign in to comment.