Only print+update wavefunction correlation energy when explicitly req…

…uested in GW. Reason: When doing RI for the exchange-self energy, Fock exchange is not computed since it includes computational overhead.
cp2k · Jun 29, 2020 · 84d4f61 · 84d4f61
1 parent d9a0891
commit 84d4f61
Show file tree

Hide file tree

Showing 4 changed files with 50 additions and 14 deletions.
diff --git a/src/input_cp2k_mp2.F b/src/input_cp2k_mp2.F
@@ -804,6 +804,15 @@ SUBROUTINE create_ri_g0w0(section)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
+      CALL keyword_create(keyword, __LOCATION__, name="UPDATE_XC_ENERGY", &
+                          description="If true, the Hartree-Fock and RPA total energy are printed and the total energy "// &
+                          "is corrected using exact exchange and the RPA correlation energy.", &
+                          usage="UPDATE_XC_ENERGY", &
+                          default_l_val=.FALSE., &
+                          lone_keyword_l_val=.TRUE.)
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
       ! here we generate a subsection for the periodic GW correction
       CALL create_periodic_gw_correction_section(subsection)
       CALL section_add_subsection(section, subsection)

diff --git a/src/mp2.F b/src/mp2.F
@@ -158,7 +158,7 @@ SUBROUTINE mp2_main(qs_env, calc_forces)
       INTEGER(KIND=int_8)                                :: mem
       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: kind_of
       LOGICAL :: calc_ex, do_admm, do_admm_rpa_exx, do_dynamic_load_balancing, do_exx, do_gw, &
-         do_kpoints_cubic_RPA, free_hfx_buffer
+         do_kpoints_cubic_RPA, free_hfx_buffer, update_xc_energy
       REAL(KIND=dp) :: E_ex_from_GW, Emp2, Emp2_AA, Emp2_AA_Cou, Emp2_AA_ex, Emp2_AB, Emp2_AB_Cou, &
          Emp2_AB_ex, Emp2_BB, Emp2_BB_Cou, Emp2_BB_ex, Emp2_Cou, Emp2_ex, Emp2_S, Emp2_T, maxocc, &
          mem_real, t1, t2, t3
@@ -722,7 +722,12 @@ SUBROUTINE mp2_main(qs_env, calc_forces)
          IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'Second order perturbation energy  =   ', Emp2
       ELSE
          IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.6)') 'Total RI-RPA Time=', t2 - t1
-         IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'RI-RPA energy  =   ', Emp2
+
+         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
             IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'RI-RPA-AXK energy=', mp2_env%ri_rpa%ener_axk
          ENDIF
@@ -1198,22 +1203,19 @@ SUBROUTINE calculate_exx(qs_env, unit_nr, do_gw, do_admm, E_ex_from_GW, t3)
          END IF
       END DO
 
-      ! 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
-
       ! take the exact exchange energy from GW or calculate it
       IF (do_gw) THEN
-         energy%total = energy%total + E_ex_from_GW
-         energy%ex = E_ex_from_GW
+
+         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
+         END IF
+
       ELSE
 
+         CALL remove_exc_energy(energy)
+
          ehfx = 0.0_dp
          DO irep = 1, n_rep_hf
             ns = SIZE(rho_ao)
@@ -1248,6 +1250,28 @@ SUBROUTINE calculate_exx(qs_env, unit_nr, do_gw, do_admm, E_ex_from_GW, t3)
 
    END SUBROUTINE calculate_exx
 
+
+   SUBROUTINE remove_exc_energy(energy)
+      TYPE(qs_energy_type), POINTER                      :: energy
+      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)
+
+         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)
+
+   END SUBROUTINE
+
 ! **************************************************************************************************
 !> \brief ...
 !> \param qs_env ...

diff --git a/src/mp2_setup.F b/src/mp2_setup.F
@@ -135,6 +135,8 @@ SUBROUTINE read_mp2_section(input, mp2_env)
                                 r_vals=mp2_env%ri_g0w0%ic_corr_list_beta)
       CALL section_vals_val_get(mp2_section, "RI_RPA%GW%GAMMA_ONLY_SIGMA", &
                                 l_val=mp2_env%ri_g0w0%do_gamma_only_sigma)
+      CALL section_vals_val_get(mp2_section, "RI_RPA%GW%UPDATE_XC_ENERGY", &
+                                l_val=mp2_env%ri_g0w0%update_xc_energy)
 
       CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE", &
                                 l_val=mp2_env%ri_g0w0%do_bse)

diff --git a/src/mp2_types.F b/src/mp2_types.F
@@ -151,6 +151,7 @@ MODULE mp2_types
       REAL(KIND=dp), DIMENSION(:), POINTER               :: ic_corr_list, ic_corr_list_beta
       INTEGER :: print_exx
       LOGICAL :: do_gamma_only_sigma
+      LOGICAL :: update_xc_energy
    END TYPE
 
    TYPE ri_basis_opt