Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Auxiliary routines needed for RPA-exchange
10 : !> given blacs_env to another
11 : !> \par History
12 : !> 09.2016 created [Vladimir Rybkin]
13 : !> 03.2019 Renamed [Frederick Stein]
14 : !> 03.2019 Moved Functions from rpa_ri_gpw.F [Frederick Stein]
15 : !> 04.2024 Added open-shell calculations, SOSEX [Frederick Stein]
16 : !> \author Vladimir Rybkin
17 : ! **************************************************************************************************
18 : MODULE rpa_exchange
19 : USE atomic_kind_types, ONLY: atomic_kind_type
20 : USE cell_types, ONLY: cell_type
21 : USE cp_blacs_env, ONLY: cp_blacs_env_type
22 : USE cp_control_types, ONLY: dft_control_type
23 : USE cp_dbcsr_api, ONLY: &
24 : dbcsr_copy, dbcsr_create, dbcsr_get_info, dbcsr_init_p, dbcsr_multiply, dbcsr_p_type, &
25 : dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
26 : USE cp_dbcsr_contrib, ONLY: dbcsr_trace
27 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
28 : USE cp_dbcsr_redistribution, ONLY: cp_array_to_dbcsr_slices
29 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale
30 : USE cp_fm_diag, ONLY: choose_eigv_solver
31 : USE cp_fm_redistribution, ONLY: cp_fm_to_array
32 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
33 : cp_fm_struct_p_type,&
34 : cp_fm_struct_release
35 : USE cp_fm_types, ONLY: cp_fm_create,&
36 : cp_fm_get_info,&
37 : cp_fm_release,&
38 : cp_fm_set_all,&
39 : cp_fm_to_fm,&
40 : cp_fm_to_fm_submat_general,&
41 : cp_fm_type
42 : USE group_dist_types, ONLY: create_group_dist,&
43 : get_group_dist,&
44 : group_dist_d1_type,&
45 : group_dist_proc,&
46 : maxsize,&
47 : release_group_dist
48 : USE hfx_admm_utils, ONLY: tddft_hfx_matrix
49 : USE hfx_types, ONLY: hfx_create,&
50 : hfx_release,&
51 : hfx_type
52 : USE input_constants, ONLY: rpa_exchange_axk,&
53 : rpa_exchange_none,&
54 : rpa_exchange_sosex
55 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
56 : section_vals_type
57 : USE kinds, ONLY: dp,&
58 : int_8
59 : USE local_gemm_api, ONLY: LOCAL_GEMM_PU_GPU
60 : USE mathconstants, ONLY: sqrthalf
61 : USE message_passing, ONLY: mp_para_env_type,&
62 : mp_proc_null
63 : USE mp2_types, ONLY: mp2_type
64 : USE parallel_gemm_api, ONLY: parallel_gemm
65 : USE particle_types, ONLY: particle_type
66 : USE qs_environment_types, ONLY: get_qs_env,&
67 : qs_environment_type
68 : USE qs_kind_types, ONLY: qs_kind_type
69 : USE qs_subsys_types, ONLY: qs_subsys_get,&
70 : qs_subsys_type
71 : USE rpa_util, ONLY: calc_fm_mat_S_rpa,&
72 : remove_scaling_factor_rpa
73 : USE scf_control_types, ONLY: scf_control_type
74 : #include "./base/base_uses.f90"
75 :
76 : IMPLICIT NONE
77 :
78 : PRIVATE
79 :
80 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_exchange'
81 :
82 : PUBLIC :: rpa_exchange_work_type, rpa_exchange_needed_mem
83 :
84 : TYPE rpa_exchange_env_type
85 : PRIVATE
86 : TYPE(qs_environment_type), POINTER :: qs_env => NULL()
87 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_hfx => NULL()
88 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dbcsr_Gamma_munu_P => NULL()
89 : TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: dbcsr_Gamma_inu_P
90 : ! Workaround GCC 8
91 : TYPE(dbcsr_type), DIMENSION(:), POINTER :: mo_coeff_o => NULL()
92 : TYPE(dbcsr_type), DIMENSION(:), POINTER :: mo_coeff_v => NULL()
93 : TYPE(dbcsr_type) :: work_ao
94 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data => NULL()
95 : TYPE(mp_para_env_type), POINTER :: para_env => NULL()
96 : TYPE(section_vals_type), POINTER :: hfx_sections => NULL()
97 : LOGICAL :: my_recalc_hfx_integrals = .FALSE.
98 : REAL(KIND=dp) :: eps_filter = 0.0_dp
99 : TYPE(cp_fm_struct_p_type), DIMENSION(:), ALLOCATABLE :: struct_Gamma
100 : CONTAINS
101 : PROCEDURE, PASS(exchange_env), NON_OVERRIDABLE :: create => hfx_create_subgroup
102 : !PROCEDURE, PASS(exchange_env), NON_OVERRIDABLE :: integrate => integrate_exchange
103 : PROCEDURE, PASS(exchange_env), NON_OVERRIDABLE :: release => hfx_release_subgroup
104 : END TYPE rpa_exchange_env_type
105 :
106 : TYPE dbcsr_matrix_p_set
107 : TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: matrix_set
108 : END TYPE dbcsr_matrix_p_set
109 :
110 : TYPE rpa_exchange_work_type
111 : PRIVATE
112 : INTEGER :: exchange_correction = rpa_exchange_none
113 : TYPE(rpa_exchange_env_type) :: exchange_env
114 : INTEGER, DIMENSION(:), ALLOCATABLE :: homo, virtual, dimen_ia
115 : TYPE(group_dist_d1_type) :: aux_func_dist = group_dist_d1_type()
116 : INTEGER, DIMENSION(:), ALLOCATABLE :: aux2send
117 : INTEGER :: dimen_RI = 0
118 : INTEGER :: block_size = 0
119 : INTEGER :: color_sub = 0
120 : INTEGER :: ngroup = 0
121 : TYPE(cp_fm_type) :: fm_mat_Q_tmp = cp_fm_type()
122 : TYPE(cp_fm_type) :: fm_mat_R_half_gemm = cp_fm_type()
123 : TYPE(cp_fm_type) :: fm_mat_U = cp_fm_type()
124 : TYPE(mp_para_env_type), POINTER :: para_env_sub => NULL()
125 : CONTAINS
126 : PROCEDURE, PUBLIC, PASS(exchange_work), NON_OVERRIDABLE :: create => rpa_exchange_work_create
127 : PROCEDURE, PUBLIC, PASS(exchange_work), NON_OVERRIDABLE :: compute => rpa_exchange_work_compute
128 : PROCEDURE, PUBLIC, PASS(exchange_work), NON_OVERRIDABLE :: release => rpa_exchange_work_release
129 : PROCEDURE, PRIVATE, PASS(exchange_work), NON_OVERRIDABLE :: redistribute_into_subgroups
130 : PROCEDURE, PRIVATE, PASS(exchange_work), NON_OVERRIDABLE :: compute_fm => rpa_exchange_work_compute_fm
131 : PROCEDURE, PRIVATE, PASS(exchange_work), NON_OVERRIDABLE :: compute_hfx => rpa_exchange_work_compute_hfx
132 : END TYPE rpa_exchange_work_type
133 :
134 : CONTAINS
135 :
136 : ! **************************************************************************************************
137 : !> \brief ...
138 : !> \param mp2_env ...
139 : !> \param homo ...
140 : !> \param virtual ...
141 : !> \param dimen_RI ...
142 : !> \param para_env ...
143 : !> \param mem_per_rank ...
144 : !> \param mem_per_repl ...
145 : ! **************************************************************************************************
146 140 : SUBROUTINE rpa_exchange_needed_mem(mp2_env, homo, virtual, dimen_RI, para_env, mem_per_rank, mem_per_repl)
147 : TYPE(mp2_type), INTENT(IN) :: mp2_env
148 : INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
149 : INTEGER, INTENT(IN) :: dimen_RI
150 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
151 : REAL(KIND=dp), INTENT(INOUT) :: mem_per_rank, mem_per_repl
152 :
153 : INTEGER :: block_size
154 :
155 : ! We need the block size and if it is unknown, an upper bound
156 140 : block_size = mp2_env%ri_rpa%exchange_block_size
157 140 : IF (block_size <= 0) block_size = MAX(1, (dimen_RI + para_env%num_pe - 1)/para_env%num_pe)
158 :
159 : ! storage of product matrix (upper bound only as it depends on the square of the potential still unknown block size)
160 310 : mem_per_rank = mem_per_rank + REAL(MAXVAL(homo), KIND=dp)**2*block_size**2*8.0_dp/(1024_dp**2)
161 :
162 : ! work arrays R (2x) and U, copies of Gamma (2x), communication buffer (as expensive as Gamma)
163 : mem_per_repl = mem_per_repl + 3.0_dp*dimen_RI*dimen_RI*8.0_dp/(1024_dp**2) &
164 310 : + 3.0_dp*MAXVAL(homo*virtual)*dimen_RI*8.0_dp/(1024_dp**2)
165 140 : END SUBROUTINE rpa_exchange_needed_mem
166 :
167 : ! **************************************************************************************************
168 : !> \brief ...
169 : !> \param exchange_work ...
170 : !> \param qs_env ...
171 : !> \param para_env_sub ...
172 : !> \param mat_munu ...
173 : !> \param dimen_RI ...
174 : !> \param fm_mat_S ...
175 : !> \param fm_mat_Q ...
176 : !> \param fm_mat_Q_gemm ...
177 : !> \param homo ...
178 : !> \param virtual ...
179 : ! **************************************************************************************************
180 170 : SUBROUTINE rpa_exchange_work_create(exchange_work, qs_env, para_env_sub, mat_munu, dimen_RI, &
181 170 : fm_mat_S, fm_mat_Q, fm_mat_Q_gemm, homo, virtual)
182 : CLASS(rpa_exchange_work_type), INTENT(INOUT) :: exchange_work
183 : TYPE(qs_environment_type), POINTER :: qs_env
184 : TYPE(mp_para_env_type), POINTER, INTENT(IN) :: para_env_sub
185 : TYPE(dbcsr_p_type), INTENT(IN) :: mat_munu
186 : INTEGER, INTENT(IN) :: dimen_RI
187 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mat_S
188 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_Q, fm_mat_Q_gemm
189 : INTEGER, DIMENSION(SIZE(fm_mat_S)), INTENT(IN) :: homo, virtual
190 :
191 : INTEGER :: nspins, aux_global, aux_local, my_process_row, proc, ispin
192 170 : INTEGER, DIMENSION(:), POINTER :: row_indices, aux_distribution_fm
193 : TYPE(cp_blacs_env_type), POINTER :: context
194 :
195 170 : exchange_work%exchange_correction = qs_env%mp2_env%ri_rpa%exchange_correction
196 :
197 170 : IF (exchange_work%exchange_correction == rpa_exchange_none) RETURN
198 :
199 : ASSOCIATE (para_env => fm_mat_S(1)%matrix_struct%para_env)
200 12 : exchange_work%para_env_sub => para_env_sub
201 12 : exchange_work%ngroup = para_env%num_pe/para_env_sub%num_pe
202 12 : exchange_work%color_sub = para_env%mepos/para_env_sub%num_pe
203 : END ASSOCIATE
204 :
205 12 : CALL cp_fm_get_info(fm_mat_S(1), row_indices=row_indices, nrow_locals=aux_distribution_fm, context=context)
206 12 : CALL context%get(my_process_row=my_process_row)
207 :
208 12 : CALL create_group_dist(exchange_work%aux_func_dist, exchange_work%ngroup, dimen_RI)
209 36 : ALLOCATE (exchange_work%aux2send(0:exchange_work%ngroup - 1))
210 36 : exchange_work%aux2send = 0
211 499 : DO aux_local = 1, aux_distribution_fm(my_process_row)
212 487 : aux_global = row_indices(aux_local)
213 487 : proc = group_dist_proc(exchange_work%aux_func_dist, aux_global)
214 499 : exchange_work%aux2send(proc) = exchange_work%aux2send(proc) + 1
215 : END DO
216 :
217 12 : nspins = SIZE(fm_mat_S)
218 :
219 60 : ALLOCATE (exchange_work%homo(nspins), exchange_work%virtual(nspins), exchange_work%dimen_ia(nspins))
220 26 : exchange_work%homo(:) = homo
221 26 : exchange_work%virtual(:) = virtual
222 26 : exchange_work%dimen_ia(:) = homo*virtual
223 12 : exchange_work%dimen_RI = dimen_RI
224 :
225 12 : exchange_work%block_size = qs_env%mp2_env%ri_rpa%exchange_block_size
226 12 : IF (exchange_work%block_size <= 0) exchange_work%block_size = dimen_RI
227 :
228 12 : CALL cp_fm_create(exchange_work%fm_mat_U, fm_mat_Q%matrix_struct, name="fm_mat_U")
229 12 : CALL cp_fm_create(exchange_work%fm_mat_Q_tmp, fm_mat_Q%matrix_struct, name="fm_mat_Q_tmp")
230 12 : CALL cp_fm_create(exchange_work%fm_mat_R_half_gemm, fm_mat_Q_gemm%matrix_struct)
231 :
232 12 : IF (qs_env%mp2_env%ri_rpa%use_hfx_implementation) THEN
233 2 : CALL exchange_work%exchange_env%create(qs_env, mat_munu%matrix, para_env_sub, fm_mat_S)
234 : END IF
235 :
236 12 : IF (ASSOCIATED(qs_env%mp2_env%ri_rpa%mo_coeff_o)) THEN
237 22 : DO ispin = 1, SIZE(qs_env%mp2_env%ri_rpa%mo_coeff_o)
238 22 : CALL dbcsr_release(qs_env%mp2_env%ri_rpa%mo_coeff_o(ispin))
239 : END DO
240 10 : DEALLOCATE (qs_env%mp2_env%ri_rpa%mo_coeff_o)
241 : END IF
242 :
243 12 : IF (ASSOCIATED(qs_env%mp2_env%ri_rpa%mo_coeff_v)) THEN
244 22 : DO ispin = 1, SIZE(qs_env%mp2_env%ri_rpa%mo_coeff_v)
245 22 : CALL dbcsr_release(qs_env%mp2_env%ri_rpa%mo_coeff_v(ispin))
246 : END DO
247 10 : DEALLOCATE (qs_env%mp2_env%ri_rpa%mo_coeff_v)
248 : END IF
249 170 : END SUBROUTINE rpa_exchange_work_create
250 :
251 : ! **************************************************************************************************
252 : !> \brief ... Initializes x_data on a subgroup
253 : !> \param exchange_env ...
254 : !> \param qs_env ...
255 : !> \param mat_munu ...
256 : !> \param para_env_sub ...
257 : !> \param fm_mat_S ...
258 : !> \author Vladimir Rybkin
259 : ! **************************************************************************************************
260 2 : SUBROUTINE hfx_create_subgroup(exchange_env, qs_env, mat_munu, para_env_sub, fm_mat_S)
261 : CLASS(rpa_exchange_env_type), INTENT(INOUT) :: exchange_env
262 : TYPE(dbcsr_type), INTENT(IN) :: mat_munu
263 : TYPE(qs_environment_type), POINTER :: qs_env
264 : TYPE(mp_para_env_type), POINTER, INTENT(IN) :: para_env_sub
265 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mat_S
266 :
267 : CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_create_subgroup'
268 :
269 : INTEGER :: handle, nelectron_total, ispin, &
270 : number_of_aos, nspins, dimen_RI, dimen_ia
271 2 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
272 : TYPE(cell_type), POINTER :: my_cell
273 : TYPE(dft_control_type), POINTER :: dft_control
274 2 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
275 2 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
276 : TYPE(qs_subsys_type), POINTER :: subsys
277 : TYPE(scf_control_type), POINTER :: scf_control
278 : TYPE(section_vals_type), POINTER :: input
279 :
280 2 : CALL timeset(routineN, handle)
281 :
282 2 : exchange_env%mo_coeff_o => qs_env%mp2_env%ri_rpa%mo_coeff_o
283 2 : exchange_env%mo_coeff_v => qs_env%mp2_env%ri_rpa%mo_coeff_v
284 2 : NULLIFY (qs_env%mp2_env%ri_rpa%mo_coeff_o, qs_env%mp2_env%ri_rpa%mo_coeff_v)
285 :
286 2 : nspins = SIZE(exchange_env%mo_coeff_o)
287 :
288 2 : exchange_env%qs_env => qs_env
289 2 : exchange_env%para_env => para_env_sub
290 2 : exchange_env%eps_filter = qs_env%mp2_env%mp2_gpw%eps_filter
291 :
292 2 : NULLIFY (my_cell, atomic_kind_set, particle_set, dft_control, qs_kind_set, scf_control)
293 :
294 : CALL get_qs_env(qs_env, &
295 : subsys=subsys, &
296 : input=input, &
297 : scf_control=scf_control, &
298 2 : nelectron_total=nelectron_total)
299 :
300 : CALL qs_subsys_get(subsys, &
301 : cell=my_cell, &
302 : atomic_kind_set=atomic_kind_set, &
303 : qs_kind_set=qs_kind_set, &
304 2 : particle_set=particle_set)
305 :
306 2 : exchange_env%hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
307 2 : CALL get_qs_env(qs_env, dft_control=dft_control)
308 :
309 : ! Retrieve particle_set and atomic_kind_set
310 : CALL hfx_create(exchange_env%x_data, para_env_sub, exchange_env%hfx_sections, atomic_kind_set, &
311 : qs_kind_set, particle_set, dft_control, my_cell, orb_basis='ORB', &
312 2 : nelectron_total=nelectron_total)
313 :
314 2 : exchange_env%my_recalc_hfx_integrals = .TRUE.
315 :
316 2 : CALL dbcsr_allocate_matrix_set(exchange_env%mat_hfx, nspins)
317 4 : DO ispin = 1, nspins
318 2 : ALLOCATE (exchange_env%mat_hfx(ispin)%matrix)
319 2 : CALL dbcsr_init_p(exchange_env%mat_hfx(ispin)%matrix)
320 : CALL dbcsr_create(exchange_env%mat_hfx(ispin)%matrix, template=mat_munu, &
321 2 : matrix_type=dbcsr_type_no_symmetry)
322 4 : CALL dbcsr_copy(exchange_env%mat_hfx(ispin)%matrix, mat_munu)
323 : END DO
324 :
325 2 : CALL dbcsr_get_info(mat_munu, nfullcols_total=number_of_aos)
326 :
327 : CALL dbcsr_create(exchange_env%work_ao, template=mat_munu, &
328 2 : matrix_type=dbcsr_type_no_symmetry)
329 :
330 8 : ALLOCATE (exchange_env%dbcsr_Gamma_inu_P(nspins))
331 2 : CALL dbcsr_allocate_matrix_set(exchange_env%dbcsr_Gamma_munu_P, nspins)
332 4 : DO ispin = 1, nspins
333 2 : ALLOCATE (exchange_env%dbcsr_Gamma_munu_P(ispin)%matrix)
334 : CALL dbcsr_create(exchange_env%dbcsr_Gamma_munu_P(ispin)%matrix, template=mat_munu, &
335 2 : matrix_type=dbcsr_type_no_symmetry)
336 2 : CALL dbcsr_copy(exchange_env%dbcsr_Gamma_munu_P(ispin)%matrix, mat_munu)
337 2 : CALL dbcsr_set(exchange_env%dbcsr_Gamma_munu_P(ispin)%matrix, 0.0_dp)
338 :
339 2 : CALL dbcsr_create(exchange_env%dbcsr_Gamma_inu_P(ispin), template=exchange_env%mo_coeff_o(ispin))
340 2 : CALL dbcsr_copy(exchange_env%dbcsr_Gamma_inu_P(ispin), exchange_env%mo_coeff_o(ispin))
341 4 : CALL dbcsr_set(exchange_env%dbcsr_Gamma_inu_P(ispin), 0.0_dp)
342 : END DO
343 :
344 8 : ALLOCATE (exchange_env%struct_Gamma(nspins))
345 4 : DO ispin = 1, nspins
346 2 : CALL cp_fm_get_info(fm_mat_S(ispin), nrow_global=dimen_RI, ncol_global=dimen_ia)
347 : CALL cp_fm_struct_create(exchange_env%struct_Gamma(ispin)%struct, template_fmstruct=fm_mat_S(ispin)%matrix_struct, &
348 4 : nrow_global=dimen_ia, ncol_global=dimen_RI)
349 : END DO
350 :
351 2 : CALL timestop(handle)
352 :
353 4 : END SUBROUTINE hfx_create_subgroup
354 :
355 : ! **************************************************************************************************
356 : !> \brief ...
357 : !> \param exchange_work ...
358 : ! **************************************************************************************************
359 170 : SUBROUTINE rpa_exchange_work_release(exchange_work)
360 : CLASS(rpa_exchange_work_type), INTENT(INOUT) :: exchange_work
361 :
362 170 : IF (ALLOCATED(exchange_work%homo)) DEALLOCATE (exchange_work%homo)
363 170 : IF (ALLOCATED(exchange_work%virtual)) DEALLOCATE (exchange_work%virtual)
364 170 : IF (ALLOCATED(exchange_work%dimen_ia)) DEALLOCATE (exchange_work%dimen_ia)
365 170 : NULLIFY (exchange_work%para_env_sub)
366 170 : CALL release_group_dist(exchange_work%aux_func_dist)
367 170 : IF (ALLOCATED(exchange_work%aux2send)) DEALLOCATE (exchange_work%aux2send)
368 170 : CALL cp_fm_release(exchange_work%fm_mat_Q_tmp)
369 170 : CALL cp_fm_release(exchange_work%fm_mat_U)
370 170 : CALL cp_fm_release(exchange_work%fm_mat_R_half_gemm)
371 :
372 170 : CALL exchange_work%exchange_env%release()
373 170 : END SUBROUTINE rpa_exchange_work_release
374 :
375 : ! **************************************************************************************************
376 : !> \brief ...
377 : !> \param exchange_env ...
378 : ! **************************************************************************************************
379 170 : SUBROUTINE hfx_release_subgroup(exchange_env)
380 : CLASS(rpa_exchange_env_type), INTENT(INOUT) :: exchange_env
381 :
382 : INTEGER :: ispin
383 :
384 170 : NULLIFY (exchange_env%para_env, exchange_env%hfx_sections)
385 :
386 170 : IF (ASSOCIATED(exchange_env%x_data)) THEN
387 2 : CALL hfx_release(exchange_env%x_data)
388 2 : NULLIFY (exchange_env%x_data)
389 : END IF
390 :
391 170 : CALL dbcsr_release(exchange_env%work_ao)
392 :
393 170 : IF (ASSOCIATED(exchange_env%dbcsr_Gamma_munu_P)) THEN
394 4 : DO ispin = 1, SIZE(exchange_env%mat_hfx, 1)
395 2 : CALL dbcsr_release(exchange_env%dbcsr_Gamma_munu_P(ispin)%matrix)
396 2 : CALL dbcsr_release(exchange_env%mat_hfx(ispin)%matrix)
397 2 : CALL dbcsr_release(exchange_env%dbcsr_Gamma_inu_P(ispin))
398 2 : CALL dbcsr_release(exchange_env%mo_coeff_o(ispin))
399 2 : CALL dbcsr_release(exchange_env%mo_coeff_v(ispin))
400 2 : DEALLOCATE (exchange_env%mat_hfx(ispin)%matrix)
401 4 : DEALLOCATE (exchange_env%dbcsr_Gamma_munu_P(ispin)%matrix)
402 : END DO
403 2 : DEALLOCATE (exchange_env%mat_hfx, exchange_env%dbcsr_Gamma_munu_P)
404 2 : DEALLOCATE (exchange_env%dbcsr_Gamma_inu_P, exchange_env%mo_coeff_o, exchange_env%mo_coeff_v)
405 2 : NULLIFY (exchange_env%mat_hfx, exchange_env%dbcsr_Gamma_munu_P)
406 : END IF
407 170 : IF (ALLOCATED(exchange_env%struct_Gamma)) THEN
408 4 : DO ispin = 1, SIZE(exchange_env%struct_Gamma)
409 4 : CALL cp_fm_struct_release(exchange_env%struct_Gamma(ispin)%struct)
410 : END DO
411 2 : DEALLOCATE (exchange_env%struct_Gamma)
412 : END IF
413 170 : END SUBROUTINE hfx_release_subgroup
414 :
415 : ! **************************************************************************************************
416 : !> \brief Main driver for RPA-exchange energies
417 : !> \param exchange_work ...
418 : !> \param fm_mat_Q ...
419 : !> \param eig ...
420 : !> \param fm_mat_S ...
421 : !> \param omega ...
422 : !> \param e_exchange_corr exchange energy correction for a quadrature point
423 : !> \param mp2_env ...
424 : !> \author Vladimir Rybkin, 07/2016
425 : ! **************************************************************************************************
426 12 : SUBROUTINE rpa_exchange_work_compute(exchange_work, fm_mat_Q, eig, fm_mat_S, omega, &
427 : e_exchange_corr, mp2_env)
428 : CLASS(rpa_exchange_work_type), INTENT(INOUT) :: exchange_work
429 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_Q
430 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: eig
431 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: fm_mat_S
432 : REAL(KIND=dp), INTENT(IN) :: omega
433 : REAL(KIND=dp), INTENT(INOUT) :: e_exchange_corr
434 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
435 :
436 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rpa_exchange_work_compute'
437 : REAL(KIND=dp), PARAMETER :: thresh = 0.0000001_dp
438 :
439 : INTEGER :: handle, nspins, dimen_RI, iiB
440 12 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenval
441 :
442 12 : IF (exchange_work%exchange_correction == rpa_exchange_none) RETURN
443 :
444 12 : CALL timeset(routineN, handle)
445 :
446 12 : CALL cp_fm_get_info(fm_mat_Q, ncol_global=dimen_RI)
447 :
448 12 : nspins = SIZE(fm_mat_S)
449 :
450 : ! Eigenvalues
451 36 : ALLOCATE (eigenval(dimen_RI))
452 12 : eigenval = 0.0_dp
453 :
454 12 : CALL cp_fm_set_all(matrix=exchange_work%fm_mat_Q_tmp, alpha=0.0_dp)
455 12 : CALL cp_fm_set_all(matrix=exchange_work%fm_mat_U, alpha=0.0_dp)
456 :
457 : ! Copy Q to Q_tmp
458 12 : CALL cp_fm_to_fm(fm_mat_Q, exchange_work%fm_mat_Q_tmp)
459 : ! Diagonalize Q
460 12 : CALL choose_eigv_solver(exchange_work%fm_mat_Q_tmp, exchange_work%fm_mat_U, eigenval)
461 :
462 : ! Calculate diagonal matrix for R_half
463 :
464 : ! Manipulate eigenvalues to get diagonal matrix
465 12 : IF (exchange_work%exchange_correction == rpa_exchange_axk) THEN
466 818 : DO iib = 1, dimen_RI
467 818 : IF (ABS(eigenval(iib)) >= thresh) THEN
468 : eigenval(iib) = &
469 : SQRT((1.0_dp/(eigenval(iib)**2))*LOG(1.0_dp + eigenval(iib)) &
470 718 : - 1.0_dp/(eigenval(iib)*(eigenval(iib) + 1.0_dp)))
471 : ELSE
472 90 : eigenval(iib) = sqrthalf
473 : END IF
474 : END DO
475 2 : ELSE IF (exchange_work%exchange_correction == rpa_exchange_sosex) THEN
476 168 : DO iib = 1, dimen_RI
477 168 : IF (ABS(eigenval(iib)) >= thresh) THEN
478 : eigenval(iib) = &
479 : SQRT(-(1.0_dp/(eigenval(iib)**2))*LOG(1.0_dp + eigenval(iib)) &
480 144 : + 1.0_dp/eigenval(iib))
481 : ELSE
482 22 : eigenval(iib) = sqrthalf
483 : END IF
484 : END DO
485 : ELSE
486 0 : CPABORT("Unknown RPA exchange correction")
487 : END IF
488 :
489 : ! fm_mat_U now contains some sqrt of the required matrix-valued function
490 12 : CALL cp_fm_column_scale(exchange_work%fm_mat_U, eigenval)
491 :
492 : ! Release memory
493 12 : DEALLOCATE (eigenval)
494 :
495 : ! Redistribute fm_mat_U for "rectangular" multiplication: ia*P P*P
496 12 : CALL cp_fm_set_all(matrix=exchange_work%fm_mat_R_half_gemm, alpha=0.0_dp)
497 :
498 : CALL cp_fm_to_fm_submat_general(exchange_work%fm_mat_U, exchange_work%fm_mat_R_half_gemm, dimen_RI, &
499 12 : dimen_RI, 1, 1, 1, 1, exchange_work%fm_mat_U%matrix_struct%context)
500 :
501 12 : IF (mp2_env%ri_rpa%use_hfx_implementation) THEN
502 2 : CALL exchange_work%compute_hfx(fm_mat_S, eig, omega, e_exchange_corr)
503 : ELSE
504 10 : CALL exchange_work%compute_fm(fm_mat_S, eig, omega, e_exchange_corr, mp2_env)
505 : END IF
506 :
507 12 : CALL timestop(handle)
508 :
509 12 : END SUBROUTINE rpa_exchange_work_compute
510 :
511 : ! **************************************************************************************************
512 : !> \brief Main driver for RPA-exchange energies
513 : !> \param exchange_work ...
514 : !> \param fm_mat_S ...
515 : !> \param eig ...
516 : !> \param omega ...
517 : !> \param e_exchange_corr exchange energy correction for a quadrature point
518 : !> \param mp2_env ...
519 : !> \author Frederick Stein, May-June 2024
520 : ! **************************************************************************************************
521 10 : SUBROUTINE rpa_exchange_work_compute_fm(exchange_work, fm_mat_S, eig, omega, &
522 : e_exchange_corr, mp2_env)
523 : CLASS(rpa_exchange_work_type), INTENT(INOUT) :: exchange_work
524 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mat_S
525 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: eig
526 : REAL(KIND=dp), INTENT(IN) :: omega
527 : REAL(KIND=dp), INTENT(INOUT) :: e_exchange_corr
528 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
529 :
530 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rpa_exchange_work_compute_fm'
531 :
532 : INTEGER :: handle, ispin, nspins, P, Q, L_size_Gamma, hom, virt, i, &
533 : send_proc, recv_proc, recv_size, max_aux_size, proc_shift, dimen_ia, &
534 : block_size, P_start, P_end, P_size, Q_start, Q_size, Q_end, handle2, my_aux_size, my_virt
535 10 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), TARGET :: mat_Gamma_3_3D
536 10 : REAL(KIND=dp), POINTER, DIMENSION(:), CONTIGUOUS :: mat_Gamma_3_1D
537 10 : REAL(KIND=dp), POINTER, DIMENSION(:, :), CONTIGUOUS :: mat_Gamma_3_2D
538 10 : REAL(KIND=dp), ALLOCATABLE, TARGET, DIMENSION(:) :: recv_buffer_1D
539 10 : REAL(KIND=dp), POINTER, DIMENSION(:, :), CONTIGUOUS :: recv_buffer_2D
540 10 : REAL(KIND=dp), POINTER, DIMENSION(:, :, :), CONTIGUOUS :: recv_buffer_3D
541 10 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mat_B_iaP
542 10 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), TARGET :: product_matrix_1D
543 10 : REAL(KIND=dp), POINTER, DIMENSION(:, :), CONTIGUOUS :: product_matrix_2D
544 10 : REAL(KIND=dp), POINTER, DIMENSION(:, :, :, :), CONTIGUOUS :: product_matrix_4D
545 : TYPE(cp_fm_type) :: fm_mat_Gamma_3
546 : TYPE(mp_para_env_type), POINTER :: para_env
547 10 : TYPE(group_dist_d1_type) :: virt_dist
548 :
549 10 : CALL timeset(routineN, handle)
550 :
551 10 : nspins = SIZE(fm_mat_S)
552 :
553 10 : CALL get_group_dist(exchange_work%aux_func_dist, exchange_work%color_sub, sizes=my_aux_size)
554 :
555 10 : e_exchange_corr = 0.0_dp
556 10 : max_aux_size = maxsize(exchange_work%aux_func_dist)
557 :
558 : ! local_gemm_ctx has a very large footprint the first time this routine is
559 : ! called.
560 10 : CALL mp2_env%local_gemm_ctx%create(LOCAL_GEMM_PU_GPU)
561 10 : CALL mp2_env%local_gemm_ctx%set_op_threshold_gpu(128*128*128*2)
562 :
563 22 : DO ispin = 1, nspins
564 12 : hom = exchange_work%homo(ispin)
565 12 : virt = exchange_work%virtual(ispin)
566 12 : dimen_ia = hom*virt
567 12 : IF (hom < 1 .OR. virt < 1) CYCLE
568 :
569 12 : CALL cp_fm_get_info(fm_mat_S(ispin), para_env=para_env)
570 :
571 12 : CALL cp_fm_create(fm_mat_Gamma_3, fm_mat_S(ispin)%matrix_struct)
572 12 : CALL cp_fm_set_all(matrix=fm_mat_Gamma_3, alpha=0.0_dp)
573 :
574 : ! Update G with a new value of Omega: in practice, it is G*S
575 :
576 : ! Scale fm_work_iaP
577 : CALL calc_fm_mat_S_rpa(fm_mat_S(ispin), .TRUE., virt, eig(:, ispin), &
578 12 : hom, omega, 0.0_dp)
579 :
580 : ! Calculate Gamma_3: Gamma_3 = G*S*R^(1/2) = G*S*R^(1/2)
581 : CALL parallel_gemm(transa="T", transb="N", m=exchange_work%dimen_RI, n=dimen_ia, k=exchange_work%dimen_RI, alpha=1.0_dp, &
582 : matrix_a=exchange_work%fm_mat_R_half_gemm, matrix_b=fm_mat_S(ispin), beta=0.0_dp, &
583 12 : matrix_c=fm_mat_Gamma_3)
584 :
585 12 : CALL create_group_dist(virt_dist, exchange_work%para_env_sub%num_pe, virt)
586 :
587 : ! Remove extra factor from S after the multiplication (to return to the original matrix)
588 12 : CALL remove_scaling_factor_rpa(fm_mat_S(ispin), virt, eig(:, ispin), hom, omega)
589 :
590 12 : CALL exchange_work%redistribute_into_subgroups(fm_mat_Gamma_3, mat_Gamma_3_3D, ispin, virt_dist)
591 12 : CALL cp_fm_release(fm_mat_Gamma_3)
592 :
593 : ! We need only the pure matrix
594 12 : CALL remove_scaling_factor_rpa(fm_mat_S(ispin), virt, eig(:, ispin), hom, omega)
595 :
596 : ! Reorder matrix from (P, i*a) -> (a, i, P) with P being distributed within subgroups
597 12 : CALL exchange_work%redistribute_into_subgroups(fm_mat_S(ispin), mat_B_iaP, ispin, virt_dist)
598 :
599 : ! Return to the original tensor
600 12 : CALL calc_fm_mat_S_rpa(fm_mat_S(ispin), .TRUE., virt, eig(:, ispin), hom, omega, 0.0_dp)
601 :
602 12 : L_size_Gamma = SIZE(mat_Gamma_3_3D, 3)
603 12 : my_virt = SIZE(mat_Gamma_3_3D, 1)
604 12 : block_size = exchange_work%block_size
605 :
606 12 : mat_Gamma_3_1D(1:INT(my_virt, KIND=int_8)*hom*my_aux_size) => mat_Gamma_3_3D(:, :, 1:my_aux_size)
607 12 : mat_Gamma_3_2D(1:my_virt, 1:hom*my_aux_size) => mat_Gamma_3_1D(1:INT(my_virt, KIND=int_8)*hom*my_aux_size)
608 :
609 0 : ALLOCATE (product_matrix_1D(INT(hom*MIN(block_size, L_size_gamma), KIND=int_8)* &
610 36 : INT(hom*MIN(block_size, max_aux_size), KIND=int_8)))
611 36 : ALLOCATE (recv_buffer_1D(INT(virt, KIND=int_8)*hom*max_aux_size))
612 12 : recv_buffer_2D(1:my_virt, 1:hom*max_aux_size) => recv_buffer_1D(1:INT(virt, KIND=int_8)*hom*max_aux_size)
613 12 : recv_buffer_3D(1:my_virt, 1:hom, 1:max_aux_size) => recv_buffer_1D(1:INT(virt, KIND=int_8)*hom*max_aux_size)
614 36 : DO proc_shift = 0, para_env%num_pe - 1, exchange_work%para_env_sub%num_pe
615 24 : send_proc = MODULO(para_env%mepos + proc_shift, para_env%num_pe)
616 24 : recv_proc = MODULO(para_env%mepos - proc_shift, para_env%num_pe)
617 :
618 24 : CALL get_group_dist(exchange_work%aux_func_dist, recv_proc/exchange_work%para_env_sub%num_pe, sizes=recv_size)
619 :
620 24 : IF (recv_size == 0) recv_proc = mp_proc_null
621 :
622 24 : CALL para_env%sendrecv(mat_B_iaP, send_proc, recv_buffer_3D(:, :, 1:recv_size), recv_proc)
623 :
624 24 : IF (recv_size == 0) CYCLE
625 :
626 1038 : DO P_start = 1, L_size_Gamma, block_size
627 1002 : P_end = MIN(L_size_Gamma, P_start + block_size - 1)
628 1002 : P_size = P_end - P_start + 1
629 43875 : DO Q_start = 1, recv_size, block_size
630 42849 : Q_end = MIN(recv_size, Q_start + block_size - 1)
631 42849 : Q_size = Q_end - Q_start + 1
632 :
633 : ! Reassign product_matrix pointers to enforce contiguity of target array
634 : product_matrix_2D(1:hom*P_size, 1:hom*Q_size) => &
635 42849 : product_matrix_1D(1:INT(hom*P_size, KIND=int_8)*INT(hom*Q_size, KIND=int_8))
636 : product_matrix_4D(1:hom, 1:P_size, 1:hom, 1:Q_size) => &
637 42849 : product_matrix_1D(1:INT(hom*P_size, KIND=int_8)*INT(hom*Q_size, KIND=int_8))
638 :
639 42849 : CALL timeset(routineN//"_gemm", handle2)
640 : CALL mp2_env%local_gemm_ctx%gemm("T", "N", hom*P_size, hom*Q_size, my_virt, 1.0_dp, &
641 : mat_Gamma_3_2D(:, hom*(P_start - 1) + 1:hom*P_end), my_virt, &
642 : recv_buffer_2D(:, hom*(Q_start - 1) + 1:hom*Q_end), my_virt, &
643 42849 : 0.0_dp, product_matrix_2D, hom*P_size)
644 42849 : CALL timestop(handle2)
645 :
646 42849 : CALL timeset(routineN//"_energy", handle2)
647 : !$OMP PARALLEL DO DEFAULT(NONE) SHARED(P_size, Q_size, hom, product_matrix_4D) &
648 42849 : !$OMP COLLAPSE(3) REDUCTION(+: e_exchange_corr) PRIVATE(P, Q, i)
649 : DO P = 1, P_size
650 : DO Q = 1, Q_size
651 : DO i = 1, hom
652 : e_exchange_corr = e_exchange_corr + DOT_PRODUCT(product_matrix_4D(i, P, :, Q), product_matrix_4D(:, P, i, Q))
653 : END DO
654 : END DO
655 : END DO
656 86700 : CALL timestop(handle2)
657 : END DO
658 : END DO
659 : END DO
660 :
661 12 : CALL release_group_dist(virt_dist)
662 12 : IF (ALLOCATED(mat_B_iaP)) DEALLOCATE (mat_B_iaP)
663 12 : IF (ALLOCATED(mat_Gamma_3_3D)) DEALLOCATE (mat_Gamma_3_3D)
664 12 : IF (ALLOCATED(product_matrix_1D)) DEALLOCATE (product_matrix_1D)
665 34 : IF (ALLOCATED(recv_buffer_1D)) DEALLOCATE (recv_buffer_1D)
666 : END DO
667 :
668 10 : CALL mp2_env%local_gemm_ctx%destroy()
669 :
670 10 : IF (nspins == 2) e_exchange_corr = e_exchange_corr*2.0_dp
671 10 : IF (nspins == 1) e_exchange_corr = e_exchange_corr*4.0_dp
672 :
673 10 : CALL timestop(handle)
674 :
675 20 : END SUBROUTINE rpa_exchange_work_compute_fm
676 :
677 : ! **************************************************************************************************
678 : !> \brief Contract RPA-exchange density matrix with HF exchange integrals and evaluate the correction
679 : !> \param exchange_work ...
680 : !> \param fm_mat_S ...
681 : !> \param eig ...
682 : !> \param omega ...
683 : !> \param e_exchange_corr ...
684 : !> \author Vladimir Rybkin, 08/2016
685 : ! **************************************************************************************************
686 2 : SUBROUTINE rpa_exchange_work_compute_hfx(exchange_work, fm_mat_S, eig, omega, e_exchange_corr)
687 : CLASS(rpa_exchange_work_type), INTENT(INOUT) :: exchange_work
688 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: fm_mat_S
689 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: eig
690 : REAL(KIND=dp), INTENT(IN) :: omega
691 : REAL(KIND=dp), INTENT(OUT) :: e_exchange_corr
692 :
693 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rpa_exchange_work_compute_hfx'
694 :
695 : INTEGER :: handle, ispin, my_aux_start, my_aux_end, &
696 : my_aux_size, nspins, L_counter, dimen_ia, hom, virt
697 : REAL(KIND=dp) :: e_exchange_P
698 2 : TYPE(dbcsr_matrix_p_set), DIMENSION(:), ALLOCATABLE :: dbcsr_Gamma_3
699 : TYPE(cp_fm_type) :: fm_mat_Gamma_3
700 : TYPE(mp_para_env_type), POINTER :: para_env
701 :
702 2 : CALL timeset(routineN, handle)
703 :
704 2 : e_exchange_corr = 0.0_dp
705 :
706 2 : nspins = SIZE(fm_mat_S)
707 :
708 2 : CALL get_group_dist(exchange_work%aux_func_dist, exchange_work%color_sub, my_aux_start, my_aux_end, my_aux_size)
709 :
710 8 : ALLOCATE (dbcsr_Gamma_3(nspins))
711 4 : DO ispin = 1, nspins
712 2 : hom = exchange_work%homo(ispin)
713 2 : virt = exchange_work%virtual(ispin)
714 2 : dimen_ia = hom*virt
715 2 : IF (hom < 1 .OR. virt < 1) CYCLE
716 :
717 2 : CALL cp_fm_get_info(fm_mat_S(ispin), para_env=para_env)
718 :
719 2 : CALL cp_fm_create(fm_mat_Gamma_3, exchange_work%exchange_env%struct_Gamma(ispin)%struct)
720 2 : CALL cp_fm_set_all(matrix=fm_mat_Gamma_3, alpha=0.0_dp)
721 :
722 : ! Update G with a new value of Omega: in practice, it is G*S
723 :
724 : ! Scale fm_work_iaP
725 : CALL calc_fm_mat_S_rpa(fm_mat_S(ispin), .TRUE., virt, eig(:, ispin), &
726 2 : hom, omega, 0.0_dp)
727 :
728 : ! Calculate Gamma_3: Gamma_3 = G*S*R^(1/2) = G*S*R^(1/2)
729 : CALL parallel_gemm(transa="T", transb="N", m=dimen_ia, n=exchange_work%dimen_RI, &
730 : k=exchange_work%dimen_RI, alpha=1.0_dp, &
731 : matrix_a=fm_mat_S(ispin), matrix_b=exchange_work%fm_mat_R_half_gemm, beta=0.0_dp, &
732 2 : matrix_c=fm_mat_Gamma_3)
733 :
734 : ! Remove extra factor from S after the multiplication (to return to the original matrix)
735 2 : CALL remove_scaling_factor_rpa(fm_mat_S(ispin), virt, eig(:, ispin), hom, omega)
736 :
737 : ! Copy Gamma_ia_P^3 to dbcsr matrix set
738 4 : BLOCK
739 : INTEGER :: first_row, last_row, nlocal
740 2 : REAL(KIND=dp), ALLOCATABLE :: slices(:, :)
741 2 : TYPE(group_dist_d1_type) :: row_distribution
742 :
743 2 : CALL create_group_dist(row_distribution, exchange_work%para_env_sub%num_pe, dimen_ia)
744 2 : CALL get_group_dist(row_distribution, exchange_work%para_env_sub%mepos, first_row, last_row, nlocal)
745 2 : CALL cp_fm_to_array(slices, first_row, last_row, my_aux_start, my_aux_end, fm_mat_Gamma_3)
746 89 : ALLOCATE (dbcsr_Gamma_3(ispin)%matrix_set(my_aux_size))
747 : CALL cp_array_to_dbcsr_slices(slices, hom, virt, exchange_work%para_env_sub, first_row, last_row, &
748 : exchange_work%exchange_env%mo_coeff_o(ispin), &
749 2 : matrices=dbcsr_Gamma_3(ispin)%matrix_set)
750 10 : CALL release_group_dist(row_distribution)
751 : END BLOCK
752 : END DO
753 :
754 85 : DO L_counter = 1, my_aux_size
755 166 : DO ispin = 1, nspins
756 : ! Do dbcsr multiplication: transform the virtual index
757 : CALL dbcsr_multiply("N", "T", 1.0_dp, exchange_work%exchange_env%mo_coeff_v(ispin), &
758 : dbcsr_Gamma_3(ispin)%matrix_set(L_counter), &
759 : 0.0_dp, exchange_work%exchange_env%dbcsr_Gamma_inu_P(ispin), &
760 83 : filter_eps=exchange_work%exchange_env%eps_filter)
761 :
762 83 : CALL dbcsr_release(dbcsr_Gamma_3(ispin)%matrix_set(L_counter))
763 :
764 : ! Do dbcsr multiplication: transform the occupied index
765 : CALL dbcsr_multiply("N", "T", 0.5_dp, exchange_work%exchange_env%dbcsr_Gamma_inu_P(ispin), &
766 : exchange_work%exchange_env%mo_coeff_o(ispin), &
767 : 0.0_dp, exchange_work%exchange_env%dbcsr_Gamma_munu_P(ispin)%matrix, &
768 83 : filter_eps=exchange_work%exchange_env%eps_filter)
769 : CALL dbcsr_multiply("N", "T", 0.5_dp, exchange_work%exchange_env%mo_coeff_o(ispin), &
770 : exchange_work%exchange_env%dbcsr_Gamma_inu_P(ispin), &
771 : 1.0_dp, exchange_work%exchange_env%dbcsr_Gamma_munu_P(ispin)%matrix, &
772 83 : filter_eps=exchange_work%exchange_env%eps_filter)
773 :
774 166 : CALL dbcsr_set(exchange_work%exchange_env%mat_hfx(ispin)%matrix, 0.0_dp)
775 : END DO
776 :
777 : CALL tddft_hfx_matrix(exchange_work%exchange_env%mat_hfx, exchange_work%exchange_env%dbcsr_Gamma_munu_P, &
778 : exchange_work%exchange_env%qs_env, .FALSE., &
779 : exchange_work%exchange_env%my_recalc_hfx_integrals, &
780 : exchange_work%exchange_env%hfx_sections, exchange_work%exchange_env%x_data, &
781 83 : exchange_work%exchange_env%para_env)
782 :
783 83 : exchange_work%exchange_env%my_recalc_hfx_integrals = .FALSE.
784 168 : DO ispin = 1, nspins
785 : CALL dbcsr_multiply("N", "T", 1.0_dp, exchange_work%exchange_env%mat_hfx(ispin)%matrix, &
786 : exchange_work%exchange_env%dbcsr_Gamma_munu_P(ispin)%matrix, &
787 83 : 0.0_dp, exchange_work%exchange_env%work_ao, filter_eps=exchange_work%exchange_env%eps_filter)
788 83 : CALL dbcsr_trace(exchange_work%exchange_env%work_ao, e_exchange_P)
789 166 : e_exchange_corr = e_exchange_corr - e_exchange_P
790 : END DO
791 : END DO
792 :
793 : IF (nspins == 2) e_exchange_corr = e_exchange_corr
794 2 : IF (nspins == 1) e_exchange_corr = e_exchange_corr*4.0_dp
795 :
796 2 : CALL timestop(handle)
797 :
798 6 : END SUBROUTINE rpa_exchange_work_compute_hfx
799 :
800 : ! **************************************************************************************************
801 : !> \brief ...
802 : !> \param exchange_work ...
803 : !> \param fm_mat ...
804 : !> \param mat ...
805 : !> \param ispin ...
806 : !> \param virt_dist ...
807 : ! **************************************************************************************************
808 24 : SUBROUTINE redistribute_into_subgroups(exchange_work, fm_mat, mat, ispin, virt_dist)
809 : CLASS(rpa_exchange_work_type), INTENT(IN) :: exchange_work
810 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat
811 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
812 : INTENT(OUT) :: mat
813 : INTEGER, INTENT(IN) :: ispin
814 : TYPE(group_dist_d1_type), INTENT(IN) :: virt_dist
815 :
816 : CHARACTER(LEN=*), PARAMETER :: routineN = 'redistribute_into_subgroups'
817 :
818 : INTEGER :: aux_counter, aux_global, aux_local, aux_proc, avirt, dimen_RI, handle, handle2, &
819 : ia_global, ia_local, iocc, max_number_recv, max_number_send, my_aux_end, my_aux_size, &
820 : my_aux_start, my_process_column, my_process_row, my_virt_end, my_virt_size, &
821 : my_virt_start, proc, proc_shift, recv_proc, send_proc, virt_counter, virt_proc, group_size
822 24 : INTEGER, ALLOCATABLE, DIMENSION(:) :: data2send, recv_col_indices, &
823 24 : recv_row_indices, send_aux_indices, send_virt_indices, virt2send
824 : INTEGER, DIMENSION(2) :: recv_shape
825 24 : INTEGER, DIMENSION(:), POINTER :: aux_distribution_fm, col_indices, &
826 24 : ia_distribution_fm, row_indices
827 24 : INTEGER, DIMENSION(:, :), POINTER :: mpi2blacs
828 24 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), TARGET :: recv_buffer, send_buffer
829 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
830 24 : POINTER :: recv_ptr, send_ptr
831 : TYPE(cp_blacs_env_type), POINTER :: context
832 : TYPE(mp_para_env_type), POINTER :: para_env
833 :
834 24 : CALL timeset(routineN, handle)
835 :
836 : CALL cp_fm_get_info(matrix=fm_mat, &
837 : nrow_locals=aux_distribution_fm, &
838 : col_indices=col_indices, &
839 : row_indices=row_indices, &
840 : ncol_locals=ia_distribution_fm, &
841 : context=context, &
842 : nrow_global=dimen_RI, &
843 24 : para_env=para_env)
844 :
845 24 : IF (exchange_work%homo(ispin) <= 0 .OR. exchange_work%virtual(ispin) <= 0) THEN
846 0 : CALL get_group_dist(virt_dist, exchange_work%para_env_sub%mepos, my_virt_start, my_virt_end, my_virt_size)
847 0 : ALLOCATE (mat(exchange_work%homo(ispin), my_virt_size, dimen_RI))
848 0 : CALL timestop(handle)
849 0 : RETURN
850 : END IF
851 :
852 24 : group_size = exchange_work%para_env_sub%num_pe
853 :
854 24 : CALL timeset(routineN//"_prep", handle2)
855 24 : CALL get_group_dist(exchange_work%aux_func_dist, exchange_work%color_sub, my_aux_start, my_aux_end, my_aux_size)
856 24 : CALL get_group_dist(virt_dist, exchange_work%para_env_sub%mepos, my_virt_start, my_virt_end, my_virt_size)
857 24 : CALL context%get(my_process_column=my_process_column, my_process_row=my_process_row, mpi2blacs=mpi2blacs)
858 :
859 : ! Determine the number of columns to send
860 120 : ALLOCATE (send_aux_indices(MAXVAL(exchange_work%aux2send)))
861 72 : ALLOCATE (virt2send(0:group_size - 1))
862 24 : virt2send = 0
863 1924 : DO ia_local = 1, ia_distribution_fm(my_process_column)
864 1900 : ia_global = col_indices(ia_local)
865 1900 : avirt = MOD(ia_global - 1, exchange_work%virtual(ispin)) + 1
866 1900 : proc = group_dist_proc(virt_dist, avirt)
867 1924 : virt2send(proc) = virt2send(proc) + 1
868 : END DO
869 :
870 72 : ALLOCATE (data2send(0:para_env%num_pe - 1))
871 72 : DO aux_proc = 0, exchange_work%ngroup - 1
872 120 : DO virt_proc = 0, group_size - 1
873 96 : data2send(aux_proc*group_size + virt_proc) = exchange_work%aux2send(aux_proc)*virt2send(virt_proc)
874 : END DO
875 : END DO
876 :
877 96 : ALLOCATE (send_virt_indices(MAXVAL(virt2send)))
878 72 : max_number_send = MAXVAL(data2send)
879 :
880 72 : ALLOCATE (send_buffer(INT(max_number_send, KIND=int_8)*exchange_work%homo(ispin)))
881 24 : max_number_recv = max_number_send
882 24 : CALL para_env%max(max_number_recv)
883 72 : ALLOCATE (recv_buffer(max_number_recv))
884 :
885 120 : ALLOCATE (mat(my_virt_size, exchange_work%homo(ispin), my_aux_size))
886 :
887 24 : CALL timestop(handle2)
888 :
889 24 : CALL timeset(routineN//"_own", handle2)
890 : ! Start with own data
891 1026 : DO aux_local = 1, aux_distribution_fm(my_process_row)
892 1002 : aux_global = row_indices(aux_local)
893 1002 : IF (aux_global < my_aux_start .OR. aux_global > my_aux_end) CYCLE
894 40848 : DO ia_local = 1, ia_distribution_fm(my_process_column)
895 40318 : ia_global = fm_mat%matrix_struct%col_indices(ia_local)
896 :
897 40318 : iocc = (ia_global - 1)/exchange_work%virtual(ispin) + 1
898 40318 : avirt = MOD(ia_global - 1, exchange_work%virtual(ispin)) + 1
899 :
900 40318 : IF (my_virt_start > avirt .OR. my_virt_end < avirt) CYCLE
901 :
902 41320 : mat(avirt - my_virt_start + 1, iocc, aux_global - my_aux_start + 1) = fm_mat%local_data(aux_local, ia_local)
903 : END DO
904 : END DO
905 24 : CALL timestop(handle2)
906 :
907 48 : DO proc_shift = 1, para_env%num_pe - 1
908 24 : send_proc = MODULO(para_env%mepos + proc_shift, para_env%num_pe)
909 24 : recv_proc = MODULO(para_env%mepos - proc_shift, para_env%num_pe)
910 :
911 24 : CALL timeset(routineN//"_pack_buffer", handle2)
912 : send_ptr(1:virt2send(MOD(send_proc, group_size)), &
913 : 1:exchange_work%aux2send(send_proc/group_size)) => &
914 : send_buffer(1:INT(virt2send(MOD(send_proc, group_size)), KIND=int_8)* &
915 24 : exchange_work%aux2send(send_proc/group_size))
916 : ! Pack send buffer
917 24 : aux_counter = 0
918 1026 : DO aux_local = 1, aux_distribution_fm(my_process_row)
919 1002 : aux_global = row_indices(aux_local)
920 1002 : proc = group_dist_proc(exchange_work%aux_func_dist, aux_global)
921 1002 : IF (proc /= send_proc/group_size) CYCLE
922 496 : aux_counter = aux_counter + 1
923 496 : virt_counter = 0
924 40016 : DO ia_local = 1, ia_distribution_fm(my_process_column)
925 39520 : ia_global = col_indices(ia_local)
926 39520 : avirt = MOD(ia_global - 1, exchange_work%virtual(ispin)) + 1
927 :
928 39520 : proc = group_dist_proc(virt_dist, avirt)
929 39520 : IF (proc /= MOD(send_proc, group_size)) CYCLE
930 39520 : virt_counter = virt_counter + 1
931 39520 : send_ptr(virt_counter, aux_counter) = fm_mat%local_data(aux_local, ia_local)
932 40016 : send_virt_indices(virt_counter) = ia_global
933 : END DO
934 1026 : send_aux_indices(aux_counter) = aux_global
935 : END DO
936 24 : CALL timestop(handle2)
937 :
938 24 : CALL timeset(routineN//"_ex_size", handle2)
939 24 : recv_shape = [1, 1]
940 72 : CALL para_env%sendrecv(SHAPE(send_ptr), send_proc, recv_shape, recv_proc)
941 24 : CALL timestop(handle2)
942 :
943 72 : IF (SIZE(send_ptr) == 0) send_proc = mp_proc_null
944 72 : IF (PRODUCT(recv_shape) == 0) recv_proc = mp_proc_null
945 :
946 24 : CALL timeset(routineN//"_ex_idx", handle2)
947 120 : ALLOCATE (recv_row_indices(recv_shape(1)), recv_col_indices(recv_shape(2)))
948 24 : CALL para_env%sendrecv(send_virt_indices(1:virt_counter), send_proc, recv_row_indices, recv_proc)
949 24 : CALL para_env%sendrecv(send_aux_indices(1:aux_counter), send_proc, recv_col_indices, recv_proc)
950 24 : CALL timestop(handle2)
951 :
952 : ! Prepare pointer to recv buffer (consider transposition while packing the send buffer)
953 24 : recv_ptr(1:recv_shape(1), 1:MAX(1, recv_shape(2))) => recv_buffer(1:recv_shape(1)*MAX(1, recv_shape(2)))
954 :
955 24 : CALL timeset(routineN//"_sendrecv", handle2)
956 : ! Perform communication
957 24 : CALL para_env%sendrecv(send_ptr, send_proc, recv_ptr, recv_proc)
958 24 : CALL timestop(handle2)
959 :
960 24 : IF (recv_proc == mp_proc_null) THEN
961 0 : DEALLOCATE (recv_row_indices, recv_col_indices)
962 0 : CYCLE
963 : END IF
964 :
965 24 : CALL timeset(routineN//"_unpack", handle2)
966 : ! Unpack receive buffer
967 520 : DO aux_local = 1, SIZE(recv_col_indices)
968 496 : aux_global = recv_col_indices(aux_local)
969 :
970 40040 : DO ia_local = 1, SIZE(recv_row_indices)
971 39520 : ia_global = recv_row_indices(ia_local)
972 :
973 39520 : iocc = (ia_global - 1)/exchange_work%virtual(ispin) + 1
974 39520 : avirt = MOD(ia_global - 1, exchange_work%virtual(ispin)) + 1
975 :
976 40016 : mat(avirt - my_virt_start + 1, iocc, aux_global - my_aux_start + 1) = recv_ptr(ia_local, aux_local)
977 : END DO
978 : END DO
979 24 : CALL timestop(handle2)
980 :
981 24 : IF (ALLOCATED(recv_row_indices)) DEALLOCATE (recv_row_indices)
982 168 : IF (ALLOCATED(recv_col_indices)) DEALLOCATE (recv_col_indices)
983 : END DO
984 :
985 24 : DEALLOCATE (send_aux_indices, send_virt_indices)
986 :
987 24 : CALL timestop(handle)
988 :
989 96 : END SUBROUTINE redistribute_into_subgroups
990 :
991 0 : END MODULE rpa_exchange
|
