hfx: some refactoring of integrate_four_center (#488)

cp2k · Jul 24, 2019 · 8139c36 · 8139c36
1 parent 4adbcf0
commit 8139c36
Show file tree

Hide file tree

Showing 5 changed files with 26 additions and 34 deletions.
diff --git a/src/hfx_admm_utils.F b/src/hfx_admm_utils.F
@@ -34,6 +34,7 @@ MODULE hfx_admm_utils
                                               dbcsr_type
    USE hfx_derivatives,                 ONLY: derivatives_four_center
    USE hfx_energy_potential,            ONLY: integrate_four_center
+   USE hfx_types,                       ONLY: hfx_type
    USE input_constants,                 ONLY: &
         do_admm_aux_exch_func_bee, do_admm_aux_exch_func_default, do_admm_aux_exch_func_none, &
         do_admm_aux_exch_func_opt, do_admm_aux_exch_func_pbex, do_potential_coulomb, &
@@ -227,6 +228,7 @@ SUBROUTINE hfx_ks_matrix(qs_env, matrix_ks, rho, energy, calculate_forces, &
                                                             rho_ao_1d
       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_h, matrix_ks_orb, rho_ao_orb
       TYPE(dft_control_type), POINTER                    :: dft_control
+      TYPE(hfx_type), DIMENSION(:, :), POINTER           :: x_data
       TYPE(pw_env_type), POINTER                         :: pw_env
       TYPE(pw_poisson_type), POINTER                     :: poisson_env
       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
@@ -255,7 +257,8 @@ SUBROUTINE hfx_ks_matrix(qs_env, matrix_ks, rho, energy, calculate_forces, &
                       matrix_ks_aux_fit=matrix_ks_aux_fit, &
                       matrix_ks_aux_fit_im=matrix_ks_aux_fit_im, &
                       matrix_ks_aux_fit_hfx=matrix_ks_aux_fit_hfx, &
-                      s_mstruct_changed=s_mstruct_changed)
+                      s_mstruct_changed=s_mstruct_changed, &
+                      x_data=x_data)
 
       nspins = dft_control%nspins
       nimages = dft_control%nimages
@@ -316,7 +319,7 @@ SUBROUTINE hfx_ks_matrix(qs_env, matrix_ks, rho, energy, calculate_forces, &
          ! Finally the real hfx calulation
          ehfx = 0.0_dp
          DO ispin = 1, mspin
-            CALL integrate_four_center(qs_env, matrix_ks_orb, eh1, rho_ao_orb, hfx_sections, &
+            CALL integrate_four_center(qs_env, x_data, matrix_ks_orb, eh1, rho_ao_orb, hfx_sections, &
                                        para_env, s_mstruct_changed, irep, distribute_fock_matrix, &
                                        ispin=ispin)
             ehfx = ehfx+eh1
@@ -365,7 +368,7 @@ SUBROUTINE hfx_ks_matrix(qs_env, matrix_ks, rho, energy, calculate_forces, &
 
             ehfxrt = 0.0_dp
             DO ispin = 1, mspin
-               CALL integrate_four_center(qs_env, matrix_ks_orb, eh1, rho_ao_orb, hfx_sections, &
+               CALL integrate_four_center(qs_env, x_data, matrix_ks_orb, eh1, rho_ao_orb, hfx_sections, &
                                           para_env, .FALSE., irep, distribute_fock_matrix, &
                                           ispin=ispin)
                ehfxrt = ehfxrt+eh1
@@ -977,6 +980,7 @@ SUBROUTINE tddft_hfx_matrix(matrix_ks, rho_ao, qs_env)
       TYPE(cp_para_env_type), POINTER                    :: para_env
       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_ks_kp, rho_ao_kp
       TYPE(dft_control_type), POINTER                    :: dft_control
+      TYPE(hfx_type), DIMENSION(:, :), POINTER           :: x_data
       TYPE(qs_energy_type), POINTER                      :: energy
       TYPE(section_vals_type), POINTER                   :: hfx_sections, input
 
@@ -989,7 +993,8 @@ SUBROUTINE tddft_hfx_matrix(matrix_ks, rho_ao, qs_env)
                       energy=energy, &
                       input=input, &
                       para_env=para_env, &
-                      s_mstruct_changed=s_mstruct_changed)
+                      s_mstruct_changed=s_mstruct_changed, &
+                      x_data=x_data)
 
       hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
       CALL section_vals_get(hfx_sections, explicit=do_hfx)
@@ -1015,7 +1020,7 @@ SUBROUTINE tddft_hfx_matrix(matrix_ks, rho_ao, qs_env)
             ! the real hfx calulation
             ehfx = 0.0_dp
             DO ispin = 1, mspin
-               CALL integrate_four_center(qs_env, matrix_ks_kp, eh1, rho_ao_kp, hfx_sections, para_env, &
+               CALL integrate_four_center(qs_env, x_data, matrix_ks_kp, eh1, rho_ao_kp, hfx_sections, para_env, &
                                           s_mstruct_changed, irep, distribute_fock_matrix, ispin=ispin)
                ehfx = ehfx+eh1
             END DO

diff --git a/src/hfx_energy_potential.F b/src/hfx_energy_potential.F
@@ -110,6 +110,7 @@ MODULE hfx_energy_potential
 !> \brief computes four center integrals for a full basis set and updates the
 !>      Kohn-Sham-Matrix and energy. Uses all 8 eri symmetries
 !> \param qs_env ...
+!> \param x_data ...
 !> \param ks_matrix ...
 !> \param ehfx energy calculated with the updated HFX matrix
 !> \param rho_ao density matrix in ao basis
@@ -120,8 +121,6 @@ MODULE hfx_energy_potential
 !> \param distribute_fock_matrix Flag that indicates whether to communicate the
 !>        new fock matrix or not
 !> \param ispin ...
-!> \param do_exx ...
-!> \param axk_x_data ...
 !> \par History
 !>      06.2007 created [Manuel Guidon]
 !>      08.2007 optimized load balance [Manuel Guidon]
@@ -132,11 +131,12 @@ MODULE hfx_energy_potential
 !>              [Tobias Binninger + Valery Weber]
 !> \author Manuel Guidon
 ! **************************************************************************************************
-   SUBROUTINE integrate_four_center(qs_env, ks_matrix, ehfx, rho_ao, hfx_section, para_env, &
+   SUBROUTINE integrate_four_center(qs_env, x_data, ks_matrix, ehfx, rho_ao, hfx_section, para_env, &
                                     geometry_did_change, irep, distribute_fock_matrix, &
-                                    ispin, do_exx, axk_x_data)
+                                    ispin)
 
       TYPE(qs_environment_type), POINTER                 :: qs_env
+      TYPE(hfx_type), DIMENSION(:, :), POINTER           :: x_data
       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: ks_matrix
       REAL(KIND=dp), INTENT(OUT)                         :: ehfx
       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: rho_ao
@@ -146,8 +146,6 @@ SUBROUTINE integrate_four_center(qs_env, ks_matrix, ehfx, rho_ao, hfx_section, p
       INTEGER                                            :: irep
       LOGICAL, INTENT(IN)                                :: distribute_fock_matrix
       INTEGER, INTENT(IN)                                :: ispin
-      LOGICAL, OPTIONAL                                  :: do_exx
-      TYPE(hfx_type), DIMENSION(:, :), OPTIONAL, POINTER :: axk_x_data
 
       CHARACTER(LEN=*), PARAMETER :: routineN = 'integrate_four_center', &
          routineP = moduleN//':'//routineN
@@ -190,7 +188,7 @@ SUBROUTINE integrate_four_center(qs_env, ks_matrix, ehfx, rho_ao, hfx_section, p
       INTEGER, SAVE                                      :: shm_number_of_p_entries
       LOGICAL :: bins_left, buffer_overflow, do_disk_storage, do_dynamic_load_balancing, do_it, &
          do_kpoints, do_p_screening, do_periodic, do_print_load_balance_info, is_anti_symmetric, &
-         ks_fully_occ, my_do_exx, my_geo_change, treat_lsd_in_core, use_disk_storage
+         ks_fully_occ, my_geo_change, treat_lsd_in_core, use_disk_storage
       LOGICAL, DIMENSION(:, :), POINTER                  :: shm_atomic_pair_list
       REAL(dp) :: afac, bintime_start, bintime_stop, cartesian_estimate, compression_factor, &
          compression_factor_disk, ene_x_aa, ene_x_aa_diag, ene_x_bb, ene_x_bb_diag, eps_schwarz, &
@@ -241,7 +239,6 @@ SUBROUTINE integrate_four_center(qs_env, ks_matrix, ehfx, rho_ao, hfx_section, p
          DIMENSION(:, :, :, :, :, :), POINTER            :: radii_pgf, screen_coeffs_pgf
       TYPE(hfx_screening_type)                           :: screening_parameter
       TYPE(hfx_task_list_type), DIMENSION(:), POINTER    :: shm_task_list, tmp_task_list
-      TYPE(hfx_type), DIMENSION(:, :), POINTER           :: x_data
       TYPE(hfx_type), POINTER                            :: actual_x_data, shm_master_x_data
       TYPE(kpoint_type), POINTER                         :: kpoints
       TYPE(pair_list_type)                               :: list_ij, list_kl
@@ -292,19 +289,6 @@ SUBROUTINE integrate_four_center(qs_env, ks_matrix, ehfx, rho_ao, hfx_section, p
          CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
       END IF
 
-      my_do_exx = .FALSE.
-      IF (PRESENT(do_exx)) my_do_exx = do_exx
-      NULLIFY (x_data)
-      IF (PRESENT(axk_x_data)) THEN ! In case of AXK, x_data is a subroutine argument
-         x_data => axk_x_data
-      ELSE
-         IF (my_do_exx) THEN
-            x_data => qs_env%mp2_env%ri_rpa%x_data
-         ELSE
-            x_data => qs_env%x_data
-         END IF ! my_do_exx
-      END IF ! AXK_x_data
-
       !! Calculate l_max used in fgamma , because init_md_ftable is definitely not thread safe
       nkind = SIZE(atomic_kind_set, 1)
       l_max = 0

diff --git a/src/mp2.F b/src/mp2.F
@@ -1210,9 +1210,10 @@ SUBROUTINE calculate_exx(qs_env, unit_nr, do_gw, do_admm, E_ex_from_GW, t3)
             ELSE
                matrix_ks_2d(1:ns, 1:1) => matrix_ks(1:ns)
             END IF
-            CALL integrate_four_center(qs_env, matrix_ks_2d, eh1, rho_ao_2d, hfx_sections, &
+            CALL integrate_four_center(qs_env, qs_env%mp2_env%ri_rpa%x_data, matrix_ks_2d, eh1, &
+                                       rho_ao_2d, hfx_sections, &
                                        para_env, .TRUE., irep, .TRUE., &
-                                       ispin=1, do_exx=.TRUE.)
+                                       ispin=1)
             ehfx = ehfx+eh1
          END DO
 
@@ -1547,9 +1548,10 @@ SUBROUTINE compute_vec_Sigma_x_minus_vxc_gw(qs_env, mp2_env, mos_mp2, energy_ex,
                   matrix_ks_2d(1:ns, 1:1) => matrix_ks(1:ns)
                END IF
 
-               CALL integrate_four_center(qs_env, matrix_ks_2d, eh1, rho_ao_2d, hfx_sections, &
+               CALL integrate_four_center(qs_env, qs_env%mp2_env%ri_rpa%x_data, matrix_ks_2d, eh1, &
+                                          rho_ao_2d, hfx_sections, &
                                           para_env, .TRUE., irep, .TRUE., &
-                                          ispin=1, do_exx=.TRUE.)
+                                          ispin=1)
                ehfx = ehfx+eh1
             END DO
          END IF

diff --git a/src/mp2_cphf.F b/src/mp2_cphf.F
@@ -964,6 +964,7 @@ SUBROUTINE cphf_like_update(qs_env, para_env, homo, virtual, dimen, &
       REAL(KIND=dp)                                      :: ehfx, total_rho
       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho_work_ao
       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mat_2d, rho_ao_2d
+      TYPE(hfx_type), DIMENSION(:, :), POINTER           :: x_data
       TYPE(qs_ks_env_type), POINTER                      :: ks_env
 
       alpha_beta = .FALSE.
@@ -975,7 +976,7 @@ SUBROUTINE cphf_like_update(qs_env, para_env, homo, virtual, dimen, &
       my_recalc_hfx_integrals = .FALSE.
       IF (PRESENT(recalc_hfx_integrals)) my_recalc_hfx_integrals = recalc_hfx_integrals
 
-      CALL get_qs_env(qs_env, ks_env=ks_env)
+      CALL get_qs_env(qs_env, ks_env=ks_env, x_data=x_data)
       ! perform back transformation
       SELECT CASE (transf_type_in)
       CASE (1)
@@ -1233,7 +1234,7 @@ SUBROUTINE cphf_like_update(qs_env, para_env, homo, virtual, dimen, &
             rho_ao_2d(1:ns, 1:1) => rho_work_ao(1:ns)
             ns = SIZE(mat_mu_nu)
             mat_2d(1:ns, 1:1) => mat_mu_nu(1:ns)
-            CALL integrate_four_center(qs_env, mat_2d, ehfx, rho_ao_2d, hfx_sections, &
+            CALL integrate_four_center(qs_env, x_data, mat_2d, ehfx, rho_ao_2d, hfx_sections, &
                                        para_env, my_recalc_hfx_integrals, irep, .TRUE., &
                                        ispin=1)
          END DO

diff --git a/src/rpa_axk.F b/src/rpa_axk.F
@@ -452,9 +452,9 @@ SUBROUTINE integrate_exchange(qs_env, dbcsr_Gamma_munu_P, mat_munu, para_env_sub
             rho_ao_2d(1:ns, 1:1) => rho_work_ao(1:ns)
 
             CALL dbcsr_set(mat_2d(1, 1)%matrix, 0.0_dp)
-            CALL integrate_four_center(qs_env, mat_2d, ehfx, rho_ao_2d, hfx_sections, &
+            CALL integrate_four_center(qs_env, x_data, mat_2d, ehfx, rho_ao_2d, hfx_sections, &
                                        para_env_sub, my_recalc_hfx_integrals, irep, .TRUE., &
-                                       ispin=1, axk_x_data=x_data)
+                                       ispin=1)
          END DO
 
          my_recalc_hfx_integrals = .FALSE.