RPA: Refactor contraction with metric

cp2k · Jul 5, 2019 · 63907af · 63907af
1 parent cf54a08
commit 63907af
Show file tree

Hide file tree

Showing 2 changed files with 55 additions and 22 deletions.
diff --git a/src/rpa_main.F b/src/rpa_main.F
@@ -37,7 +37,6 @@ MODULE rpa_main
                                               cp_fm_set_element,&
                                               cp_fm_to_fm,&
                                               cp_fm_type
-   USE cp_gemm_interface,               ONLY: cp_gemm
    USE cp_para_env,                     ONLY: cp_para_env_create,&
                                               cp_para_env_release
    USE cp_para_types,                   ONLY: cp_para_env_type
@@ -88,6 +87,7 @@ MODULE rpa_main
                                               RPA_postprocessing_start,&
                                               alloc_im_time,&
                                               calc_mat_Q,&
+                                              contract_P_omega_with_mat_L,&
                                               dealloc_im_time,&
                                               get_mat_3c_overl_int_cut
    USE util,                            ONLY: get_limit
@@ -1934,27 +1934,14 @@ SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_s
 
                IF (.NOT. (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma)) THEN
 
-                  ! multiplication with RI metric/Coulomb operator
-                  CALL dbcsr_multiply("N", "T", 1.0_dp, mat_P_omega(jquad, 1)%matrix, mat_L%matrix, &
-                                      0.0_dp, mat_P_global%matrix, filter_eps=eps_filter_im_time)
+                  CALL contract_P_omega_with_mat_L(mat_P_omega(jquad, 1)%matrix, mat_L%matrix, mat_P_global%matrix, &
+                                                   eps_filter_im_time, fm_mat_work, dimen_RI, fm_mat_L(1, 1)%matrix, fm_mat_Q)
 
-                  CALL copy_dbcsr_to_fm(mat_P_global%matrix, fm_mat_work)
-
-                  CALL cp_gemm('N', 'N', dimen_RI, dimen_RI, dimen_RI, 1.0_dp, fm_mat_L(1, 1)%matrix, fm_mat_work, &
-                               0.0_dp, fm_mat_Q)
-
-                  ! For open-shell SOS-MP2 we have two different matrices we have to deal with
+                  ! For open-shell SOS-MP2 we have two different matrices to deal with
                   IF (my_open_shell .AND. do_ri_sos_laplace_mp2) THEN
-                     CALL cp_fm_set_all(matrix=fm_mat_work, alpha=0.0_dp)
-
-                     ! multiplication with RI metric/Coulomb operator
-                     CALL dbcsr_multiply("N", "T", 1.0_dp, mat_P_omega_beta(jquad, 1)%matrix, mat_L%matrix, &
-                                         0.0_dp, mat_P_global%matrix, filter_eps=eps_filter_im_time)
-
-                     CALL copy_dbcsr_to_fm(mat_P_global%matrix, fm_mat_work)
-
-                     CALL cp_gemm('N', 'N', dimen_RI, dimen_RI, dimen_RI, 1.0_dp, fm_mat_L(1, 1)%matrix, fm_mat_work, &
-                                  0.0_dp, fm_mat_Q_beta)
+                     CALL contract_P_omega_with_mat_L(mat_P_omega_beta(jquad, 1)%matrix, mat_L%matrix, mat_P_global%matrix, &
+                                                      eps_filter_im_time, fm_mat_work, dimen_RI, fm_mat_L(1, 1)%matrix, &
+                                                      fm_mat_Q_beta)
                   END IF
                END IF
 

diff --git a/src/rpa_util.F b/src/rpa_util.F
@@ -18,7 +18,8 @@ MODULE rpa_util
                                               cp_cfm_release,&
                                               cp_cfm_set_all,&
                                               cp_cfm_type
-   USE cp_dbcsr_operations,             ONLY: copy_fm_to_dbcsr,&
+   USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
+                                              copy_fm_to_dbcsr,&
                                               dbcsr_allocate_matrix_set,&
                                               dbcsr_deallocate_matrix_set
    USE cp_fm_basic_linalg,              ONLY: cp_fm_syrk,&
@@ -73,7 +74,8 @@ MODULE rpa_util
 
    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_util'
 
-   PUBLIC :: RPA_postprocessing_nokp, alloc_im_time, calc_mat_Q, RPA_postprocessing_start, get_mat_3c_overl_int_cut, dealloc_im_time
+   PUBLIC :: RPA_postprocessing_nokp, alloc_im_time, calc_mat_Q, RPA_postprocessing_start, get_mat_3c_overl_int_cut, &
+             dealloc_im_time, contract_P_omega_with_mat_L
 
 CONTAINS
 
@@ -1917,4 +1919,48 @@ SUBROUTINE dealloc_im_time(do_mao, do_dbcsr_t, do_ri_sos_laplace_mp2, fm_mo_coef
 
    END SUBROUTINE dealloc_im_time
 
+! **************************************************************************************************
+!> \brief ...
+!> \param mat_P_omega ...
+!> \param mat_L ...
+!> \param mat_work ...
+!> \param eps_filter_im_time ...
+!> \param fm_mat_work ...
+!> \param dimen_RI ...
+!> \param fm_mat_L ...
+!> \param fm_mat_Q ...
+! **************************************************************************************************
+   SUBROUTINE contract_P_omega_with_mat_L(mat_P_omega, mat_L, mat_work, eps_filter_im_time, fm_mat_work, dimen_RI, &
+                                          fm_mat_L, fm_mat_Q)
+
+      TYPE(dbcsr_type), POINTER                          :: mat_P_omega, mat_L, mat_work
+      REAL(KIND=dp), INTENT(IN)                          :: eps_filter_im_time
+      TYPE(cp_fm_type), POINTER                          :: fm_mat_work
+      INTEGER, INTENT(IN)                                :: dimen_RI
+      TYPE(cp_fm_type), POINTER                          :: fm_mat_L, fm_mat_Q
+
+      CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_P_omega_with_mat_L', &
+         routineP = moduleN//':'//routineN
+
+      INTEGER                                            :: handle
+
+      CALL timeset(routineN, handle)
+
+      ! multiplication with RI metric/Coulomb operator
+      CALL dbcsr_multiply("N", "T", 1.0_dp, mat_P_omega, mat_L, &
+                          0.0_dp, mat_work, filter_eps=eps_filter_im_time)
+
+      CALL copy_dbcsr_to_fm(mat_work, fm_mat_work)
+
+      CALL cp_gemm('N', 'N', dimen_RI, dimen_RI, dimen_RI, 1.0_dp, fm_mat_L, fm_mat_work, &
+                   0.0_dp, fm_mat_Q)
+
+      ! Reset mat_work to save memory
+      CALL dbcsr_set(mat_work, 0.0_dp)
+      CALL dbcsr_filter(mat_work, 1.0_dp)
+
+      CALL timestop(handle)
+
+   END SUBROUTINE contract_P_omega_with_mat_L
+
 END MODULE rpa_util