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 Routines to calculate RI-RPA and SOS-MP2 gradients
10 : !> \par History
11 : !> 10.2021 created [Frederick Stein]
12 : ! **************************************************************************************************
13 : MODULE rpa_grad
14 : USE cp_array_utils, ONLY: cp_1d_r_cp_type,&
15 : cp_3d_r_cp_type
16 : USE cp_blacs_env, ONLY: cp_blacs_env_type
17 : USE cp_fm_basic_linalg, ONLY: cp_fm_geadd,&
18 : cp_fm_scale_and_add,&
19 : cp_fm_uplo_to_full
20 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_invert
21 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
22 : cp_fm_struct_get,&
23 : cp_fm_struct_release,&
24 : cp_fm_struct_type
25 : USE cp_fm_types, ONLY: cp_fm_create,&
26 : cp_fm_get_info,&
27 : cp_fm_release,&
28 : cp_fm_set_all,&
29 : cp_fm_to_fm,&
30 : cp_fm_to_fm_submat_general,&
31 : cp_fm_type
32 : USE dgemm_counter_types, ONLY: dgemm_counter_start,&
33 : dgemm_counter_stop,&
34 : dgemm_counter_type
35 : USE group_dist_types, ONLY: create_group_dist,&
36 : get_group_dist,&
37 : group_dist_d1_type,&
38 : group_dist_proc,&
39 : maxsize,&
40 : release_group_dist
41 : USE kahan_sum, ONLY: accurate_dot_product,&
42 : accurate_dot_product_2
43 : USE kinds, ONLY: dp,&
44 : int_8
45 : USE libint_2c_3c, ONLY: compare_potential_types
46 : USE local_gemm_api, ONLY: LOCAL_GEMM_PU_GPU,&
47 : local_gemm_ctxt_type
48 : USE machine, ONLY: m_flush,&
49 : m_memory
50 : USE mathconstants, ONLY: pi
51 : USE message_passing, ONLY: mp_comm_type,&
52 : mp_para_env_type,&
53 : mp_request_null,&
54 : mp_request_type,&
55 : mp_waitall,&
56 : mp_waitany
57 : USE mp2_laplace, ONLY: calc_fm_mat_s_laplace
58 : USE mp2_ri_grad_util, ONLY: array2fm,&
59 : create_dbcsr_gamma,&
60 : fm2array,&
61 : prepare_redistribution
62 : USE mp2_types, ONLY: mp2_type,&
63 : one_dim_int_array,&
64 : two_dim_int_array,&
65 : two_dim_real_array
66 : USE parallel_gemm_api, ONLY: parallel_gemm
67 : USE qs_environment_types, ONLY: get_qs_env,&
68 : qs_environment_type
69 : USE rpa_util, ONLY: calc_fm_mat_S_rpa,&
70 : remove_scaling_factor_rpa
71 : USE util, ONLY: get_limit
72 : #include "./base/base_uses.f90"
73 :
74 : IMPLICIT NONE
75 :
76 : PRIVATE
77 :
78 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_grad'
79 :
80 : PUBLIC :: rpa_grad_needed_mem, rpa_grad_type, rpa_grad_create, rpa_grad_finalize, rpa_grad_matrix_operations, rpa_grad_copy_Q
81 :
82 : TYPE sos_mp2_grad_work_type
83 : PRIVATE
84 : INTEGER, DIMENSION(:, :), ALLOCATABLE :: pair_list
85 : TYPE(one_dim_int_array), DIMENSION(:), ALLOCATABLE :: index2send, index2recv
86 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: P
87 : END TYPE sos_mp2_grad_work_type
88 :
89 : TYPE rpa_grad_work_type
90 : TYPE(cp_fm_type) :: fm_mat_Q_copy = cp_fm_type()
91 : TYPE(one_dim_int_array), DIMENSION(:, :), ALLOCATABLE :: index2send
92 : TYPE(two_dim_int_array), DIMENSION(:, :), ALLOCATABLE :: index2recv
93 : TYPE(group_dist_d1_type), DIMENSION(:), ALLOCATABLE :: gd_homo, gd_virtual
94 : INTEGER, DIMENSION(2) :: grid = -1, mepos = -1
95 : TYPE(two_dim_real_array), DIMENSION(:), ALLOCATABLE :: P_ij, P_ab
96 : END TYPE rpa_grad_work_type
97 :
98 : TYPE rpa_grad_type
99 : PRIVATE
100 : TYPE(cp_fm_type) :: fm_Gamma_PQ = cp_fm_type()
101 : TYPE(cp_fm_type), DIMENSION(:), ALLOCATABLE :: fm_Y
102 : TYPE(sos_mp2_grad_work_type), ALLOCATABLE, DIMENSION(:) :: sos_mp2_work_occ, sos_mp2_work_virt
103 : TYPE(rpa_grad_work_type) :: rpa_work
104 : END TYPE rpa_grad_type
105 :
106 : INTEGER, PARAMETER :: spla_threshold = 128*128*128*2
107 : INTEGER, PARAMETER :: blksize_threshold = 4
108 :
109 : CONTAINS
110 :
111 : ! **************************************************************************************************
112 : !> \brief Calculates the necessary minimum memory for the Gradient code ion MiB
113 : !> \param homo ...
114 : !> \param virtual ...
115 : !> \param dimen_RI ...
116 : !> \param mem_per_rank ...
117 : !> \param mem_per_repl ...
118 : !> \param do_ri_sos_laplace_mp2 ...
119 : !> \return ...
120 : ! **************************************************************************************************
121 44 : PURE SUBROUTINE rpa_grad_needed_mem(homo, virtual, dimen_RI, mem_per_rank, mem_per_repl, do_ri_sos_laplace_mp2)
122 : INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
123 : INTEGER, INTENT(IN) :: dimen_RI
124 : REAL(KIND=dp), INTENT(INOUT) :: mem_per_rank, mem_per_repl
125 : LOGICAL, INTENT(IN) :: do_ri_sos_laplace_mp2
126 :
127 : REAL(KIND=dp) :: mem_iaK, mem_KL, mem_pab, mem_pij
128 :
129 96 : mem_iaK = SUM(REAL(virtual, KIND=dp)*homo)*dimen_RI
130 96 : mem_pij = SUM(REAL(homo, KIND=dp)**2)
131 96 : mem_pab = SUM(REAL(virtual, KIND=dp)**2)
132 44 : mem_KL = REAL(dimen_RI, KIND=dp)*dimen_RI
133 :
134 : ! Required matrices iaK
135 : ! Ytot_iaP = sum_tau Y_iaP(tau)
136 : ! Y_iaP(tau) = S_iaP(tau)*Q_PQ(tau) (work array)
137 : ! Required matrices density matrices
138 : ! Pij (local)
139 : ! Pab (local)
140 : ! Additionally with SOS-MP2
141 : ! Send and receive buffers for degenerate orbital pairs (rough estimate: everything)
142 : ! Additionally with RPA
143 : ! copy of work matrix
144 : ! receive buffer for calculation of density matrix
145 : ! copy of matrix Q
146 44 : mem_per_rank = mem_per_rank + (mem_pij + mem_pab)*8.0_dp/(1024**2)
147 44 : mem_per_repl = mem_per_repl + (mem_iaK + 2.0_dp*mem_iaK/SIZE(homo) + mem_KL)*8.0_dp/(1024**2)
148 44 : IF (.NOT. do_ri_sos_laplace_mp2) THEN
149 24 : mem_per_repl = mem_per_rank + (mem_iaK/SIZE(homo) + mem_KL)*8.0_dp/(1024**2)
150 : END IF
151 :
152 44 : END SUBROUTINE rpa_grad_needed_mem
153 :
154 : ! **************************************************************************************************
155 : !> \brief Creates the arrays of a rpa_grad_type
156 : !> \param rpa_grad ...
157 : !> \param fm_mat_Q ...
158 : !> \param fm_mat_S ...
159 : !> \param homo ...
160 : !> \param virtual ...
161 : !> \param mp2_env ...
162 : !> \param Eigenval ...
163 : !> \param unit_nr ...
164 : !> \param do_ri_sos_laplace_mp2 ...
165 : ! **************************************************************************************************
166 308 : SUBROUTINE rpa_grad_create(rpa_grad, fm_mat_Q, fm_mat_S, &
167 44 : homo, virtual, mp2_env, Eigenval, unit_nr, do_ri_sos_laplace_mp2)
168 : TYPE(rpa_grad_type), INTENT(OUT) :: rpa_grad
169 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_Q
170 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mat_S
171 : INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
172 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
173 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: Eigenval
174 : INTEGER, INTENT(IN) :: unit_nr
175 : LOGICAL, INTENT(IN) :: do_ri_sos_laplace_mp2
176 :
177 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rpa_grad_create'
178 :
179 : INTEGER :: handle, ispin, nrow_local, nspins
180 :
181 44 : CALL timeset(routineN, handle)
182 :
183 44 : CALL cp_fm_create(rpa_grad%fm_Gamma_PQ, matrix_struct=fm_mat_Q%matrix_struct)
184 44 : CALL cp_fm_set_all(rpa_grad%fm_Gamma_PQ, 0.0_dp)
185 :
186 44 : nspins = SIZE(fm_mat_S)
187 :
188 184 : ALLOCATE (rpa_grad%fm_Y(nspins))
189 96 : DO ispin = 1, nspins
190 96 : CALL cp_fm_create(rpa_grad%fm_Y(ispin), fm_mat_S(ispin)%matrix_struct, set_zero=.TRUE.)
191 : END DO
192 :
193 44 : IF (do_ri_sos_laplace_mp2) THEN
194 : CALL sos_mp2_work_type_create(rpa_grad%sos_mp2_work_occ, rpa_grad%sos_mp2_work_virt, &
195 20 : unit_nr, Eigenval, homo, virtual, mp2_env%ri_grad%eps_canonical, fm_mat_S)
196 : ELSE
197 24 : CALL rpa_work_type_create(rpa_grad%rpa_work, fm_mat_Q, fm_mat_S, homo, virtual)
198 : END IF
199 :
200 : ! Set blocksize
201 44 : CALL cp_fm_struct_get(fm_mat_S(1)%matrix_struct, nrow_local=nrow_local)
202 44 : IF (mp2_env%ri_grad%dot_blksize < 1) mp2_env%ri_grad%dot_blksize = nrow_local
203 44 : mp2_env%ri_grad%dot_blksize = MIN(mp2_env%ri_grad%dot_blksize, nrow_local)
204 44 : IF (unit_nr > 0) THEN
205 22 : WRITE (unit_nr, '(T3,A,T75,I6)') 'GRAD_INFO| Block size for the contraction:', mp2_env%ri_grad%dot_blksize
206 22 : CALL m_flush(unit_nr)
207 : END IF
208 44 : CALL fm_mat_S(1)%matrix_struct%para_env%sync()
209 :
210 44 : CALL timestop(handle)
211 :
212 88 : END SUBROUTINE rpa_grad_create
213 :
214 : ! **************************************************************************************************
215 : !> \brief ...
216 : !> \param sos_mp2_work_occ ...
217 : !> \param sos_mp2_work_virt ...
218 : !> \param unit_nr ...
219 : !> \param Eigenval ...
220 : !> \param homo ...
221 : !> \param virtual ...
222 : !> \param eps_degenerate ...
223 : !> \param fm_mat_S ...
224 : ! **************************************************************************************************
225 20 : SUBROUTINE sos_mp2_work_type_create(sos_mp2_work_occ, sos_mp2_work_virt, unit_nr, &
226 20 : Eigenval, homo, virtual, eps_degenerate, fm_mat_S)
227 : TYPE(sos_mp2_grad_work_type), ALLOCATABLE, &
228 : DIMENSION(:), INTENT(OUT) :: sos_mp2_work_occ, sos_mp2_work_virt
229 : INTEGER, INTENT(IN) :: unit_nr
230 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: Eigenval
231 : INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
232 : REAL(KIND=dp), INTENT(IN) :: eps_degenerate
233 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mat_S
234 :
235 : CHARACTER(LEN=*), PARAMETER :: routineN = 'sos_mp2_work_type_create'
236 :
237 : INTEGER :: handle, ispin, nspins
238 :
239 20 : CALL timeset(routineN, handle)
240 :
241 20 : nspins = SIZE(fm_mat_S)
242 128 : ALLOCATE (sos_mp2_work_occ(nspins), sos_mp2_work_virt(nspins))
243 44 : DO ispin = 1, nspins
244 :
245 : CALL create_list_nearly_degen_pairs(Eigenval(1:homo(ispin), ispin), &
246 24 : eps_degenerate, sos_mp2_work_occ(ispin)%pair_list)
247 24 : IF (unit_nr > 0) WRITE (unit_nr, "(T3,A,T75,i6)") &
248 12 : "MO_INFO| Number of ij pairs below EPS_CANONICAL:", SIZE(sos_mp2_work_occ(ispin)%pair_list, 2)
249 72 : ALLOCATE (sos_mp2_work_occ(ispin)%P(homo(ispin) + SIZE(sos_mp2_work_occ(ispin)%pair_list, 2)))
250 164 : sos_mp2_work_occ(ispin)%P = 0.0_dp
251 24 : CALL prepare_comm_Pij(sos_mp2_work_occ(ispin), virtual(ispin), fm_mat_S(ispin))
252 :
253 : CALL create_list_nearly_degen_pairs(Eigenval(homo(ispin) + 1:, ispin), &
254 24 : eps_degenerate, sos_mp2_work_virt(ispin)%pair_list)
255 24 : IF (unit_nr > 0) WRITE (unit_nr, "(T3,A,T75,i6)") &
256 12 : "MO_INFO| Number of ab pairs below EPS_CANONICAL:", SIZE(sos_mp2_work_virt(ispin)%pair_list, 2)
257 72 : ALLOCATE (sos_mp2_work_virt(ispin)%P(virtual(ispin) + SIZE(sos_mp2_work_virt(ispin)%pair_list, 2)))
258 1136 : sos_mp2_work_virt(ispin)%P = 0.0_dp
259 44 : CALL prepare_comm_Pab(sos_mp2_work_virt(ispin), virtual(ispin), fm_mat_S(ispin))
260 : END DO
261 :
262 20 : CALL timestop(handle)
263 :
264 20 : END SUBROUTINE sos_mp2_work_type_create
265 :
266 : ! **************************************************************************************************
267 : !> \brief ...
268 : !> \param rpa_work ...
269 : !> \param fm_mat_Q ...
270 : !> \param fm_mat_S ...
271 : !> \param homo ...
272 : !> \param virtual ...
273 : ! **************************************************************************************************
274 144 : SUBROUTINE rpa_work_type_create(rpa_work, fm_mat_Q, fm_mat_S, homo, virtual)
275 : TYPE(rpa_grad_work_type), INTENT(OUT) :: rpa_work
276 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_Q
277 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mat_S
278 : INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
279 :
280 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rpa_work_type_create'
281 :
282 : INTEGER :: avirt, col_global, col_local, handle, iocc, ispin, my_a, my_a_end, my_a_size, &
283 : my_a_start, my_i, my_i_end, my_i_size, my_i_start, my_pcol, ncol_local, nspins, &
284 : num_pe_col, proc_homo, proc_homo_send, proc_recv, proc_send, proc_virtual, &
285 : proc_virtual_send
286 24 : INTEGER, ALLOCATABLE, DIMENSION(:) :: data2recv, data2send
287 24 : INTEGER, DIMENSION(:), POINTER :: col_indices
288 :
289 24 : CALL timeset(routineN, handle)
290 :
291 24 : CALL cp_fm_create(rpa_work%fm_mat_Q_copy, matrix_struct=fm_mat_Q%matrix_struct)
292 :
293 24 : CALL fm_mat_S(1)%matrix_struct%context%get(number_of_process_columns=num_pe_col, my_process_column=my_pcol)
294 :
295 24 : nspins = SIZE(fm_mat_S)
296 :
297 0 : ALLOCATE (rpa_work%index2send(0:num_pe_col - 1, nspins), &
298 0 : rpa_work%index2recv(0:num_pe_col - 1, nspins), &
299 0 : rpa_work%gd_homo(nspins), rpa_work%gd_virtual(nspins), &
300 : data2send(0:num_pe_col - 1), data2recv(0:num_pe_col - 1), &
301 608 : rpa_work%P_ij(nspins), rpa_work%P_ab(nspins))
302 :
303 : ! Determine new process grid
304 24 : proc_homo = MAX(1, CEILING(SQRT(REAL(num_pe_col, KIND=dp))))
305 24 : DO WHILE (MOD(num_pe_col, proc_homo) /= 0)
306 0 : proc_homo = proc_homo - 1
307 : END DO
308 24 : proc_virtual = num_pe_col/proc_homo
309 :
310 24 : rpa_work%grid(1) = proc_virtual
311 24 : rpa_work%grid(2) = proc_homo
312 :
313 24 : rpa_work%mepos(1) = MOD(my_pcol, proc_virtual)
314 24 : rpa_work%mepos(2) = my_pcol/proc_virtual
315 :
316 52 : DO ispin = 1, nspins
317 :
318 : ! Determine distributions of the orbitals
319 28 : CALL create_group_dist(rpa_work%gd_homo(ispin), proc_homo, homo(ispin))
320 28 : CALL create_group_dist(rpa_work%gd_virtual(ispin), proc_virtual, virtual(ispin))
321 :
322 28 : CALL cp_fm_struct_get(fm_mat_S(ispin)%matrix_struct, ncol_local=ncol_local, col_indices=col_indices)
323 :
324 28 : data2send = 0
325 : ! Count the amount of data2send to each process
326 2200 : DO col_local = 1, ncol_local
327 2172 : col_global = col_indices(col_local)
328 :
329 2172 : iocc = (col_global - 1)/virtual(ispin) + 1
330 2172 : avirt = col_global - (iocc - 1)*virtual(ispin)
331 :
332 2172 : proc_homo_send = group_dist_proc(rpa_work%gd_homo(ispin), iocc)
333 2172 : proc_virtual_send = group_dist_proc(rpa_work%gd_virtual(ispin), avirt)
334 :
335 2172 : proc_send = proc_homo_send*proc_virtual + proc_virtual_send
336 :
337 2200 : data2send(proc_send) = data2send(proc_send) + 1
338 : END DO
339 :
340 56 : DO proc_send = 0, num_pe_col - 1
341 112 : ALLOCATE (rpa_work%index2send(proc_send, ispin)%array(data2send(proc_send)))
342 : END DO
343 :
344 : ! Prepare the indices
345 28 : data2send = 0
346 2200 : DO col_local = 1, ncol_local
347 2172 : col_global = col_indices(col_local)
348 :
349 2172 : iocc = (col_global - 1)/virtual(ispin) + 1
350 2172 : avirt = col_global - (iocc - 1)*virtual(ispin)
351 :
352 2172 : proc_homo_send = group_dist_proc(rpa_work%gd_homo(ispin), iocc)
353 2172 : proc_virtual_send = group_dist_proc(rpa_work%gd_virtual(ispin), avirt)
354 :
355 2172 : proc_send = proc_homo_send*proc_virtual + proc_virtual_send
356 :
357 2172 : data2send(proc_send) = data2send(proc_send) + 1
358 :
359 2200 : rpa_work%index2send(proc_send, ispin)%array(data2send(proc_send)) = col_local
360 : END DO
361 :
362 : ! Count the amount of data2recv from each process
363 28 : CALL get_group_dist(rpa_work%gd_homo(ispin), my_pcol/proc_virtual, my_i_start, my_i_end, my_i_size)
364 28 : CALL get_group_dist(rpa_work%gd_virtual(ispin), MOD(my_pcol, proc_virtual), my_a_start, my_a_end, my_a_size)
365 :
366 28 : data2recv = 0
367 136 : DO my_i = my_i_start, my_i_end
368 2308 : DO my_a = my_a_start, my_a_end
369 2172 : proc_recv = fm_mat_S(ispin)%matrix_struct%g2p_col((my_i - 1)*virtual(ispin) + my_a)
370 2280 : data2recv(proc_recv) = data2recv(proc_recv) + 1
371 : END DO
372 : END DO
373 :
374 56 : DO proc_recv = 0, num_pe_col - 1
375 112 : ALLOCATE (rpa_work%index2recv(proc_recv, ispin)%array(2, data2recv(proc_recv)))
376 : END DO
377 :
378 28 : data2recv = 0
379 136 : DO my_i = my_i_start, my_i_end
380 2308 : DO my_a = my_a_start, my_a_end
381 2172 : proc_recv = fm_mat_S(ispin)%matrix_struct%g2p_col((my_i - 1)*virtual(ispin) + my_a)
382 2172 : data2recv(proc_recv) = data2recv(proc_recv) + 1
383 :
384 2172 : rpa_work%index2recv(proc_recv, ispin)%array(2, data2recv(proc_recv)) = my_i - my_i_start + 1
385 2280 : rpa_work%index2recv(proc_recv, ispin)%array(1, data2recv(proc_recv)) = my_a - my_a_start + 1
386 : END DO
387 : END DO
388 :
389 0 : ALLOCATE (rpa_work%P_ij(ispin)%array(my_i_size, homo(ispin)), &
390 196 : rpa_work%P_ab(ispin)%array(my_a_size, virtual(ispin)))
391 556 : rpa_work%P_ij(ispin)%array(:, :) = 0.0_dp
392 12384 : rpa_work%P_ab(ispin)%array(:, :) = 0.0_dp
393 :
394 : END DO
395 :
396 24 : DEALLOCATE (data2send, data2recv)
397 :
398 24 : CALL timestop(handle)
399 :
400 48 : END SUBROUTINE rpa_work_type_create
401 :
402 : ! **************************************************************************************************
403 : !> \brief ...
404 : !> \param Eigenval ...
405 : !> \param eps_degen ...
406 : !> \param pair_list ...
407 : ! **************************************************************************************************
408 48 : SUBROUTINE create_list_nearly_degen_pairs(Eigenval, eps_degen, pair_list)
409 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval
410 : REAL(KIND=dp), INTENT(IN) :: eps_degen
411 : INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(OUT) :: pair_list
412 :
413 : INTEGER :: my_i, my_j, num_orbitals, num_pairs, &
414 : pair_counter
415 :
416 48 : num_orbitals = SIZE(Eigenval)
417 :
418 : ! Determine number of nearly degenerate orbital pairs
419 : ! Trivial cases: diagonal elements
420 48 : num_pairs = 0
421 640 : DO my_i = 1, num_orbitals
422 11576 : DO my_j = 1, num_orbitals
423 10936 : IF (my_i == my_j) CYCLE
424 10936 : IF (ABS(Eigenval(my_i) - Eigenval(my_j)) < eps_degen) num_pairs = num_pairs + 1
425 : END DO
426 : END DO
427 104 : ALLOCATE (pair_list(2, num_pairs))
428 :
429 : ! Print the required pairs
430 48 : pair_counter = 1
431 640 : DO my_i = 1, num_orbitals
432 11576 : DO my_j = 1, num_orbitals
433 10936 : IF (my_i == my_j) CYCLE
434 10936 : IF (ABS(Eigenval(my_i) - Eigenval(my_j)) < eps_degen) THEN
435 660 : pair_list(1, pair_counter) = my_i
436 660 : pair_list(2, pair_counter) = my_j
437 660 : pair_counter = pair_counter + 1
438 : END IF
439 : END DO
440 : END DO
441 :
442 48 : END SUBROUTINE create_list_nearly_degen_pairs
443 :
444 : ! **************************************************************************************************
445 : !> \brief ...
446 : !> \param sos_mp2_work ...
447 : !> \param virtual ...
448 : !> \param fm_mat_S ...
449 : ! **************************************************************************************************
450 24 : SUBROUTINE prepare_comm_Pij(sos_mp2_work, virtual, fm_mat_S)
451 : TYPE(sos_mp2_grad_work_type), INTENT(INOUT) :: sos_mp2_work
452 : INTEGER, INTENT(IN) :: virtual
453 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S
454 :
455 : CHARACTER(LEN=*), PARAMETER :: routineN = 'prepare_comm_Pij'
456 :
457 : INTEGER :: avirt, col_global, col_local, counter, handle, ij_counter, iocc, my_i, my_j, &
458 : my_pcol, my_prow, ncol_local, nrow_local, num_ij_pairs, num_pe_col, pcol, pcol_recv, &
459 : pcol_send, proc_shift, tag
460 : INTEGER, ALLOCATABLE, DIMENSION(:) :: data2recv, data2send
461 24 : INTEGER, DIMENSION(:), POINTER :: col_indices, ncol_locals
462 24 : INTEGER, DIMENSION(:, :), POINTER :: blacs2mpi
463 : TYPE(cp_blacs_env_type), POINTER :: context
464 : TYPE(mp_comm_type) :: comm_exchange
465 : TYPE(mp_para_env_type), POINTER :: para_env
466 :
467 24 : CALL timeset(routineN, handle)
468 :
469 24 : tag = 44
470 :
471 24 : CALL fm_mat_S%matrix_struct%context%get(number_of_process_columns=num_pe_col)
472 0 : ALLOCATE (sos_mp2_work%index2send(0:num_pe_col - 1), &
473 144 : sos_mp2_work%index2recv(0:num_pe_col - 1))
474 :
475 72 : ALLOCATE (data2send(0:num_pe_col - 1))
476 48 : ALLOCATE (data2recv(0:num_pe_col - 1))
477 :
478 : CALL cp_fm_struct_get(fm_mat_S%matrix_struct, para_env=para_env, ncol_locals=ncol_locals, &
479 : ncol_local=ncol_local, col_indices=col_indices, &
480 24 : context=context, nrow_local=nrow_local)
481 : CALL context%get(my_process_row=my_prow, my_process_column=my_pcol, &
482 24 : blacs2mpi=blacs2mpi)
483 :
484 24 : num_ij_pairs = SIZE(sos_mp2_work%pair_list, 2)
485 :
486 24 : IF (num_ij_pairs > 0) THEN
487 :
488 4 : CALL comm_exchange%from_split(para_env, my_prow)
489 :
490 4 : data2send = 0
491 4 : data2recv = 0
492 :
493 8 : DO proc_shift = 0, num_pe_col - 1
494 4 : pcol_send = MODULO(my_pcol + proc_shift, num_pe_col)
495 :
496 4 : counter = 0
497 308 : DO col_local = 1, ncol_local
498 304 : col_global = col_indices(col_local)
499 :
500 304 : iocc = MAX(1, col_global - 1)/virtual + 1
501 304 : avirt = col_global - (iocc - 1)*virtual
502 :
503 764 : DO ij_counter = 1, num_ij_pairs
504 :
505 760 : my_i = sos_mp2_work%pair_list(1, ij_counter)
506 760 : my_j = sos_mp2_work%pair_list(2, ij_counter)
507 :
508 760 : IF (iocc /= my_j) CYCLE
509 304 : pcol = fm_mat_S%matrix_struct%g2p_col((my_i - 1)*virtual + avirt)
510 304 : IF (pcol /= pcol_send) CYCLE
511 :
512 304 : counter = counter + 1
513 :
514 760 : EXIT
515 :
516 : END DO
517 : END DO
518 8 : data2send(pcol_send) = counter
519 : END DO
520 :
521 4 : CALL comm_exchange%alltoall(data2send, data2recv, 1)
522 :
523 8 : DO proc_shift = 0, num_pe_col - 1
524 4 : pcol_send = MODULO(my_pcol + proc_shift, num_pe_col)
525 4 : pcol_recv = MODULO(my_pcol - proc_shift, num_pe_col)
526 :
527 : ! Collect indices and exchange
528 12 : ALLOCATE (sos_mp2_work%index2send(pcol_send)%array(data2send(pcol_send)))
529 :
530 4 : counter = 0
531 308 : DO col_local = 1, ncol_local
532 304 : col_global = col_indices(col_local)
533 :
534 304 : iocc = MAX(1, col_global - 1)/virtual + 1
535 304 : avirt = col_global - (iocc - 1)*virtual
536 :
537 764 : DO ij_counter = 1, num_ij_pairs
538 :
539 760 : my_i = sos_mp2_work%pair_list(1, ij_counter)
540 760 : my_j = sos_mp2_work%pair_list(2, ij_counter)
541 :
542 760 : IF (iocc /= my_j) CYCLE
543 304 : pcol = fm_mat_S%matrix_struct%g2p_col((my_i - 1)*virtual + avirt)
544 304 : IF (pcol /= pcol_send) CYCLE
545 :
546 304 : counter = counter + 1
547 :
548 304 : sos_mp2_work%index2send(pcol_send)%array(counter) = col_global
549 :
550 760 : EXIT
551 :
552 : END DO
553 : END DO
554 :
555 12 : ALLOCATE (sos_mp2_work%index2recv(pcol_recv)%array(data2recv(pcol_recv)))
556 : !
557 : CALL para_env%sendrecv(sos_mp2_work%index2send(pcol_send)%array, blacs2mpi(my_prow, pcol_send), &
558 4 : sos_mp2_work%index2recv(pcol_recv)%array, blacs2mpi(my_prow, pcol_recv), tag)
559 :
560 : ! Convert to global coordinates to local coordinates as we always work with them
561 312 : DO counter = 1, data2send(pcol_send)
562 : sos_mp2_work%index2send(pcol_send)%array(counter) = &
563 308 : fm_mat_S%matrix_struct%g2l_col(sos_mp2_work%index2send(pcol_send)%array(counter))
564 : END DO
565 : END DO
566 :
567 4 : CALL comm_exchange%free()
568 : END IF
569 :
570 24 : DEALLOCATE (data2send, data2recv)
571 :
572 24 : CALL timestop(handle)
573 :
574 48 : END SUBROUTINE prepare_comm_Pij
575 :
576 : ! **************************************************************************************************
577 : !> \brief ...
578 : !> \param sos_mp2_work ...
579 : !> \param virtual ...
580 : !> \param fm_mat_S ...
581 : ! **************************************************************************************************
582 24 : SUBROUTINE prepare_comm_Pab(sos_mp2_work, virtual, fm_mat_S)
583 : TYPE(sos_mp2_grad_work_type), INTENT(INOUT) :: sos_mp2_work
584 : INTEGER, INTENT(IN) :: virtual
585 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S
586 :
587 : CHARACTER(LEN=*), PARAMETER :: routineN = 'prepare_comm_Pab'
588 :
589 : INTEGER :: ab_counter, avirt, col_global, col_local, counter, handle, iocc, my_a, my_b, &
590 : my_pcol, my_prow, ncol_local, nrow_local, num_ab_pairs, num_pe_col, pcol, pcol_recv, &
591 : pcol_send, proc_shift, tag
592 24 : INTEGER, ALLOCATABLE, DIMENSION(:) :: data2recv, data2send
593 24 : INTEGER, DIMENSION(:), POINTER :: col_indices, ncol_locals
594 24 : INTEGER, DIMENSION(:, :), POINTER :: blacs2mpi
595 : TYPE(cp_blacs_env_type), POINTER :: context
596 : TYPE(mp_comm_type) :: comm_exchange
597 : TYPE(mp_para_env_type), POINTER :: para_env
598 :
599 24 : CALL timeset(routineN, handle)
600 :
601 24 : tag = 44
602 :
603 24 : CALL fm_mat_S%matrix_struct%context%get(number_of_process_columns=num_pe_col)
604 0 : ALLOCATE (sos_mp2_work%index2send(0:num_pe_col - 1), &
605 144 : sos_mp2_work%index2recv(0:num_pe_col - 1))
606 :
607 24 : num_ab_pairs = SIZE(sos_mp2_work%pair_list, 2)
608 24 : IF (num_ab_pairs > 0) THEN
609 :
610 : CALL cp_fm_struct_get(fm_mat_S%matrix_struct, para_env=para_env, ncol_locals=ncol_locals, &
611 : ncol_local=ncol_local, col_indices=col_indices, &
612 4 : context=context, nrow_local=nrow_local)
613 : CALL context%get(my_process_row=my_prow, my_process_column=my_pcol, &
614 4 : blacs2mpi=blacs2mpi)
615 :
616 4 : CALL comm_exchange%from_split(para_env, my_prow)
617 :
618 12 : ALLOCATE (data2send(0:num_pe_col - 1))
619 8 : ALLOCATE (data2recv(0:num_pe_col - 1))
620 :
621 4 : data2send = 0
622 4 : data2recv = 0
623 8 : DO proc_shift = 0, num_pe_col - 1
624 4 : pcol_send = MODULO(my_pcol + proc_shift, num_pe_col)
625 4 : pcol_recv = MODULO(my_pcol - proc_shift, num_pe_col)
626 :
627 4 : counter = 0
628 308 : DO col_local = 1, ncol_local
629 304 : col_global = col_indices(col_local)
630 :
631 304 : iocc = MAX(1, col_global - 1)/virtual + 1
632 304 : avirt = col_global - (iocc - 1)*virtual
633 :
634 15476 : DO ab_counter = 1, num_ab_pairs
635 :
636 15472 : my_a = sos_mp2_work%pair_list(1, ab_counter)
637 15472 : my_b = sos_mp2_work%pair_list(2, ab_counter)
638 :
639 15472 : IF (avirt /= my_b) CYCLE
640 304 : pcol = fm_mat_S%matrix_struct%g2p_col((iocc - 1)*virtual + my_a)
641 304 : IF (pcol /= pcol_send) CYCLE
642 :
643 304 : counter = counter + 1
644 :
645 15472 : EXIT
646 :
647 : END DO
648 : END DO
649 8 : data2send(pcol_send) = counter
650 : END DO
651 :
652 4 : CALL comm_exchange%alltoall(data2send, data2recv, 1)
653 :
654 8 : DO proc_shift = 0, num_pe_col - 1
655 4 : pcol_send = MODULO(my_pcol + proc_shift, num_pe_col)
656 4 : pcol_recv = MODULO(my_pcol - proc_shift, num_pe_col)
657 :
658 : ! Collect indices and exchange
659 12 : ALLOCATE (sos_mp2_work%index2send(pcol_send)%array(data2send(pcol_send)))
660 :
661 4 : counter = 0
662 308 : DO col_local = 1, ncol_local
663 304 : col_global = col_indices(col_local)
664 :
665 304 : iocc = MAX(1, col_global - 1)/virtual + 1
666 304 : avirt = col_global - (iocc - 1)*virtual
667 :
668 15476 : DO ab_counter = 1, num_ab_pairs
669 :
670 15472 : my_a = sos_mp2_work%pair_list(1, ab_counter)
671 15472 : my_b = sos_mp2_work%pair_list(2, ab_counter)
672 :
673 15472 : IF (avirt /= my_b) CYCLE
674 304 : pcol = fm_mat_S%matrix_struct%g2p_col((iocc - 1)*virtual + my_a)
675 304 : IF (pcol /= pcol_send) CYCLE
676 :
677 304 : counter = counter + 1
678 :
679 304 : sos_mp2_work%index2send(pcol_send)%array(counter) = col_global
680 :
681 15472 : EXIT
682 :
683 : END DO
684 : END DO
685 :
686 12 : ALLOCATE (sos_mp2_work%index2recv(pcol_recv)%array(data2recv(pcol_recv)))
687 : !
688 : CALL para_env%sendrecv(sos_mp2_work%index2send(pcol_send)%array, blacs2mpi(my_prow, pcol_send), &
689 4 : sos_mp2_work%index2recv(pcol_recv)%array, blacs2mpi(my_prow, pcol_recv), tag)
690 :
691 : ! Convert to global coordinates to local coordinates as we always work with them
692 312 : DO counter = 1, data2send(pcol_send)
693 : sos_mp2_work%index2send(pcol_send)%array(counter) = &
694 308 : fm_mat_S%matrix_struct%g2l_col(sos_mp2_work%index2send(pcol_send)%array(counter))
695 : END DO
696 : END DO
697 :
698 4 : CALL comm_exchange%free()
699 8 : DEALLOCATE (data2send, data2recv)
700 :
701 : END IF
702 :
703 24 : CALL timestop(handle)
704 :
705 48 : END SUBROUTINE prepare_comm_Pab
706 :
707 : ! **************************************************************************************************
708 : !> \brief ...
709 : !> \param fm_mat_Q ...
710 : !> \param rpa_grad ...
711 : ! **************************************************************************************************
712 56 : SUBROUTINE rpa_grad_copy_Q(fm_mat_Q, rpa_grad)
713 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_Q
714 : TYPE(rpa_grad_type), INTENT(INOUT) :: rpa_grad
715 :
716 56 : CALL cp_fm_to_fm(fm_mat_Q, rpa_grad%rpa_work%fm_mat_Q_copy)
717 :
718 56 : END SUBROUTINE rpa_grad_copy_Q
719 :
720 : ! **************************************************************************************************
721 : !> \brief ...
722 : !> \param mp2_env ...
723 : !> \param rpa_grad ...
724 : !> \param do_ri_sos_laplace_mp2 ...
725 : !> \param fm_mat_Q ...
726 : !> \param fm_mat_Q_gemm ...
727 : !> \param dgemm_counter ...
728 : !> \param fm_mat_S ...
729 : !> \param omega ...
730 : !> \param homo ...
731 : !> \param virtual ...
732 : !> \param Eigenval ...
733 : !> \param weight ...
734 : !> \param unit_nr ...
735 : ! **************************************************************************************************
736 106 : SUBROUTINE rpa_grad_matrix_operations(mp2_env, rpa_grad, do_ri_sos_laplace_mp2, fm_mat_Q, fm_mat_Q_gemm, &
737 106 : dgemm_counter, fm_mat_S, omega, homo, virtual, Eigenval, weight, unit_nr)
738 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
739 : TYPE(rpa_grad_type), INTENT(INOUT) :: rpa_grad
740 : LOGICAL, INTENT(IN) :: do_ri_sos_laplace_mp2
741 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mat_Q, fm_mat_Q_gemm
742 : TYPE(dgemm_counter_type), INTENT(INOUT) :: dgemm_counter
743 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mat_S
744 : REAL(KIND=dp), INTENT(IN) :: omega
745 : INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
746 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: Eigenval
747 : REAL(KIND=dp), INTENT(IN) :: weight
748 : INTEGER, INTENT(IN) :: unit_nr
749 :
750 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rpa_grad_matrix_operations'
751 :
752 : INTEGER :: col_global, col_local, dimen_ia, &
753 : dimen_RI, handle, handle2, ispin, &
754 : jspin, ncol_local, nrow_local, nspins, &
755 : row_local
756 106 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
757 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
758 106 : TARGET :: mat_S_3D, mat_work_iaP_3D
759 : TYPE(cp_fm_type) :: fm_work_iaP, fm_work_PQ
760 :
761 106 : CALL timeset(routineN, handle)
762 :
763 106 : nspins = SIZE(fm_mat_Q)
764 :
765 : CALL cp_fm_get_info(fm_mat_Q(1), nrow_global=dimen_RI, nrow_local=nrow_local, ncol_local=ncol_local, &
766 106 : col_indices=col_indices, row_indices=row_indices)
767 :
768 106 : IF (.NOT. do_ri_sos_laplace_mp2) THEN
769 56 : CALL cp_fm_create(fm_work_PQ, fm_mat_Q(1)%matrix_struct)
770 :
771 : ! calculate [1+Q(iw')]^-1
772 56 : CALL cp_fm_cholesky_invert(fm_mat_Q(1))
773 : ! symmetrize the result, fm_work_PQ is only a work matrix
774 56 : CALL cp_fm_uplo_to_full(fm_mat_Q(1), fm_work_PQ)
775 :
776 56 : CALL cp_fm_release(fm_work_PQ)
777 :
778 4816 : DO col_local = 1, ncol_local
779 4760 : col_global = col_indices(col_local)
780 191944 : DO row_local = 1, nrow_local
781 191888 : IF (col_global == row_indices(row_local)) THEN
782 4096 : fm_mat_Q(1)%local_data(row_local, col_local) = fm_mat_Q(1)%local_data(row_local, col_local) - 1.0_dp
783 4096 : EXIT
784 : END IF
785 : END DO
786 : END DO
787 :
788 56 : CALL timeset(routineN//"_PQ", handle2)
789 56 : CALL dgemm_counter_start(dgemm_counter)
790 : CALL parallel_gemm(transa="N", transb="N", m=dimen_RI, n=dimen_RI, k=dimen_RI, alpha=weight, &
791 : matrix_a=rpa_grad%rpa_work%fm_mat_Q_copy, matrix_b=fm_mat_Q(1), beta=1.0_dp, &
792 56 : matrix_c=rpa_grad%fm_Gamma_PQ)
793 56 : CALL dgemm_counter_stop(dgemm_counter, dimen_RI, dimen_RI, dimen_RI)
794 56 : CALL timestop(handle2)
795 :
796 : CALL cp_fm_to_fm_submat_general(fm_mat_Q(1), fm_mat_Q_gemm(1), dimen_RI, dimen_RI, 1, 1, 1, 1, &
797 56 : fm_mat_Q_gemm(1)%matrix_struct%context)
798 : END IF
799 :
800 228 : DO ispin = 1, nspins
801 122 : IF (do_ri_sos_laplace_mp2) THEN
802 : ! The spin of the other Q matrix is always the other spin
803 58 : jspin = nspins - ispin + 1
804 : ELSE
805 : ! or the first matrix in the case of RPA
806 : jspin = 1
807 : END IF
808 :
809 122 : IF (do_ri_sos_laplace_mp2) THEN
810 58 : CALL timeset(routineN//"_PQ", handle2)
811 58 : CALL dgemm_counter_start(dgemm_counter)
812 : CALL parallel_gemm(transa="N", transb="N", m=dimen_RI, n=dimen_RI, k=dimen_RI, alpha=weight, &
813 : matrix_a=fm_mat_Q(ispin), matrix_b=fm_mat_Q(jspin), beta=1.0_dp, &
814 58 : matrix_c=rpa_grad%fm_Gamma_PQ)
815 58 : CALL dgemm_counter_stop(dgemm_counter, dimen_RI, dimen_RI, dimen_RI)
816 58 : CALL timestop(handle2)
817 :
818 : CALL cp_fm_to_fm_submat_general(fm_mat_Q(jspin), fm_mat_Q_gemm(jspin), dimen_RI, dimen_RI, 1, 1, 1, 1, &
819 58 : fm_mat_Q_gemm(jspin)%matrix_struct%context)
820 : ELSE
821 : CALL calc_fm_mat_S_rpa(fm_mat_S(ispin), .TRUE., virtual(ispin), Eigenval(:, ispin), &
822 64 : homo(ispin), omega, 0.0_dp)
823 : END IF
824 :
825 122 : CALL timeset(routineN//"_contr_S", handle2)
826 122 : CALL cp_fm_create(fm_work_iaP, rpa_grad%fm_Y(ispin)%matrix_struct)
827 :
828 122 : CALL cp_fm_get_info(fm_mat_S(ispin), ncol_global=dimen_ia)
829 :
830 122 : CALL dgemm_counter_start(dgemm_counter)
831 : CALL parallel_gemm(transa="N", transb="N", m=dimen_RI, n=dimen_ia, k=dimen_RI, alpha=1.0_dp, &
832 : matrix_a=fm_mat_Q_gemm(jspin), matrix_b=fm_mat_S(ispin), beta=0.0_dp, &
833 122 : matrix_c=fm_work_iaP)
834 122 : CALL dgemm_counter_stop(dgemm_counter, dimen_ia, dimen_RI, dimen_RI)
835 122 : CALL timestop(handle2)
836 :
837 350 : IF (do_ri_sos_laplace_mp2) THEN
838 : CALL calc_P_sos_mp2(homo(ispin), fm_mat_S(ispin), fm_work_iaP, &
839 : rpa_grad%sos_mp2_work_occ(ispin), rpa_grad%sos_mp2_work_virt(ispin), &
840 58 : omega, weight, virtual(ispin), Eigenval(:, ispin), mp2_env%ri_grad%dot_blksize)
841 :
842 58 : CALL calc_fm_mat_S_laplace(fm_work_iaP, homo(ispin), virtual(ispin), Eigenval(:, ispin), omega)
843 :
844 58 : CALL cp_fm_scale_and_add(1.0_dp, rpa_grad%fm_Y(ispin), -weight, fm_work_iaP)
845 :
846 58 : CALL cp_fm_release(fm_work_iaP)
847 : ELSE
848 : ! To save memory, we add it now
849 64 : CALL cp_fm_scale_and_add(1.0_dp, rpa_grad%fm_Y(ispin), -weight, fm_work_iaP)
850 :
851 : ! Redistribute both matrices and deallocate fm_work_iaP
852 : CALL redistribute_fm_mat_S(rpa_grad%rpa_work%index2send(:, ispin), rpa_grad%rpa_work%index2recv(:, ispin), &
853 : fm_work_iaP, mat_work_iaP_3D, &
854 : rpa_grad%rpa_work%gd_homo(ispin), rpa_grad%rpa_work%gd_virtual(ispin), &
855 64 : rpa_grad%rpa_work%mepos)
856 64 : CALL cp_fm_release(fm_work_iaP)
857 :
858 : CALL redistribute_fm_mat_S(rpa_grad%rpa_work%index2send(:, ispin), rpa_grad%rpa_work%index2recv(:, ispin), &
859 : fm_mat_S(ispin), mat_S_3D, &
860 : rpa_grad%rpa_work%gd_homo(ispin), rpa_grad%rpa_work%gd_virtual(ispin), &
861 64 : rpa_grad%rpa_work%mepos)
862 :
863 : ! Now collect the density matrix
864 : CALL calc_P_rpa(mat_S_3D, mat_work_iaP_3D, rpa_grad%rpa_work%gd_homo(ispin), rpa_grad%rpa_work%gd_virtual(ispin), &
865 : rpa_grad%rpa_work%grid, rpa_grad%rpa_work%mepos, &
866 : fm_mat_S(ispin)%matrix_struct, &
867 : rpa_grad%rpa_work%P_ij(ispin)%array, rpa_grad%rpa_work%P_ab(ispin)%array, &
868 64 : weight, omega, Eigenval(:, ispin), homo(ispin), unit_nr, mp2_env)
869 :
870 64 : DEALLOCATE (mat_work_iaP_3D, mat_S_3D)
871 :
872 64 : CALL remove_scaling_factor_rpa(fm_mat_S(ispin), virtual(ispin), Eigenval(:, ispin), homo(ispin), omega)
873 :
874 : END IF
875 :
876 : END DO
877 :
878 106 : CALL timestop(handle)
879 :
880 212 : END SUBROUTINE rpa_grad_matrix_operations
881 :
882 : ! **************************************************************************************************
883 : !> \brief ...
884 : !> \param homo ...
885 : !> \param fm_mat_S ...
886 : !> \param fm_work_iaP ...
887 : !> \param sos_mp2_work_occ ...
888 : !> \param sos_mp2_work_virt ...
889 : !> \param omega ...
890 : !> \param weight ...
891 : !> \param virtual ...
892 : !> \param Eigenval ...
893 : !> \param dot_blksize ...
894 : ! **************************************************************************************************
895 290 : SUBROUTINE calc_P_sos_mp2(homo, fm_mat_S, fm_work_iaP, sos_mp2_work_occ, sos_mp2_work_virt, &
896 58 : omega, weight, virtual, Eigenval, dot_blksize)
897 : INTEGER, INTENT(IN) :: homo
898 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S, fm_work_iaP
899 : TYPE(sos_mp2_grad_work_type), INTENT(INOUT) :: sos_mp2_work_occ, sos_mp2_work_virt
900 : REAL(KIND=dp), INTENT(IN) :: omega, weight
901 : INTEGER, INTENT(IN) :: virtual
902 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval
903 : INTEGER, INTENT(IN) :: dot_blksize
904 :
905 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_P_sos_mp2'
906 :
907 : INTEGER :: avirt, col_global, col_local, handle, &
908 : handle2, iocc, my_a, my_i, ncol_local, &
909 : nrow_local, num_ab_pairs, num_ij_pairs
910 58 : INTEGER, DIMENSION(:), POINTER :: col_indices
911 : REAL(KIND=dp) :: ddot, trace
912 :
913 58 : CALL timeset(routineN, handle)
914 :
915 58 : CALL cp_fm_get_info(fm_mat_S, col_indices=col_indices, ncol_local=ncol_local, nrow_local=nrow_local)
916 :
917 58 : CALL timeset(routineN//"_Pij_diag", handle2)
918 282 : DO my_i = 1, homo
919 : ! Collect the contributions of the matrix elements
920 :
921 224 : trace = 0.0_dp
922 :
923 17864 : DO col_local = 1, ncol_local
924 17640 : col_global = col_indices(col_local)
925 :
926 17640 : iocc = MAX(1, col_global - 1)/virtual + 1
927 17640 : avirt = col_global - (iocc - 1)*virtual
928 :
929 17640 : IF (iocc == my_i) trace = trace + &
930 4784 : ddot(nrow_local, fm_mat_S%local_data(:, col_local), 1, fm_work_iaP%local_data(:, col_local), 1)
931 : END DO
932 :
933 282 : sos_mp2_work_occ%P(my_i) = sos_mp2_work_occ%P(my_i) - trace*omega*weight
934 :
935 : END DO
936 58 : CALL timestop(handle2)
937 :
938 58 : CALL timeset(routineN//"_Pab_diag", handle2)
939 1248 : DO my_a = 1, virtual
940 : ! Collect the contributions of the matrix elements
941 :
942 1190 : trace = 0.0_dp
943 :
944 95270 : DO col_local = 1, ncol_local
945 94080 : col_global = col_indices(col_local)
946 :
947 94080 : iocc = MAX(1, col_global - 1)/virtual + 1
948 94080 : avirt = col_global - (iocc - 1)*virtual
949 :
950 94080 : IF (avirt == my_a) trace = trace + &
951 5750 : ddot(nrow_local, fm_mat_S%local_data(:, col_local), 1, fm_work_iaP%local_data(:, col_local), 1)
952 : END DO
953 :
954 1248 : sos_mp2_work_virt%P(my_a) = sos_mp2_work_virt%P(my_a) + trace*omega*weight
955 :
956 : END DO
957 58 : CALL timestop(handle2)
958 :
959 : ! Loop over list and carry out operations
960 58 : num_ij_pairs = SIZE(sos_mp2_work_occ%pair_list, 2)
961 58 : num_ab_pairs = SIZE(sos_mp2_work_virt%pair_list, 2)
962 58 : IF (num_ij_pairs > 0) THEN
963 : CALL calc_Pij_degen(fm_work_iaP, fm_mat_S, sos_mp2_work_occ%pair_list, &
964 : virtual, sos_mp2_work_occ%P(homo + 1:), Eigenval(:homo), omega, weight, &
965 8 : sos_mp2_work_occ%index2send, sos_mp2_work_occ%index2recv, dot_blksize)
966 : END IF
967 58 : IF (num_ab_pairs > 0) THEN
968 : CALL calc_Pab_degen(fm_work_iaP, fm_mat_S, sos_mp2_work_virt%pair_list, &
969 : virtual, sos_mp2_work_virt%P(virtual + 1:), Eigenval(homo + 1:), omega, weight, &
970 8 : sos_mp2_work_virt%index2send, sos_mp2_work_virt%index2recv, dot_blksize)
971 : END IF
972 :
973 58 : CALL timestop(handle)
974 :
975 58 : END SUBROUTINE calc_P_sos_mp2
976 :
977 : ! **************************************************************************************************
978 : !> \brief ...
979 : !> \param mat_S_1D ...
980 : !> \param mat_work_iaP_3D ...
981 : !> \param gd_homo ...
982 : !> \param gd_virtual ...
983 : !> \param grid ...
984 : !> \param mepos ...
985 : !> \param fm_struct_S ...
986 : !> \param P_ij ...
987 : !> \param P_ab ...
988 : !> \param weight ...
989 : !> \param omega ...
990 : !> \param Eigenval ...
991 : !> \param homo ...
992 : !> \param unit_nr ...
993 : !> \param mp2_env ...
994 : ! **************************************************************************************************
995 64 : SUBROUTINE calc_P_rpa(mat_S_1D, mat_work_iaP_3D, gd_homo, gd_virtual, grid, mepos, &
996 64 : fm_struct_S, P_ij, P_ab, weight, omega, Eigenval, homo, unit_nr, mp2_env)
997 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT), TARGET :: mat_S_1D
998 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: mat_work_iaP_3D
999 : TYPE(group_dist_d1_type), INTENT(IN) :: gd_homo, gd_virtual
1000 : INTEGER, DIMENSION(2), INTENT(IN) :: grid, mepos
1001 : TYPE(cp_fm_struct_type), INTENT(IN), POINTER :: fm_struct_S
1002 : REAL(KIND=dp), DIMENSION(:, :) :: P_ij, P_ab
1003 : REAL(KIND=dp), INTENT(IN) :: weight, omega
1004 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval
1005 : INTEGER, INTENT(IN) :: homo, unit_nr
1006 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
1007 :
1008 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_P_rpa'
1009 :
1010 : INTEGER :: completed, handle, handle2, my_a_end, my_a_size, my_a_start, my_i_end, my_i_size, &
1011 : my_i_start, my_P_size, my_prow, number_of_parallel_channels, proc_a_recv, proc_a_send, &
1012 : proc_i_recv, proc_i_send, proc_recv, proc_send, proc_shift, recv_a_end, recv_a_size, &
1013 : recv_a_start, recv_i_end, recv_i_size, recv_i_start, tag
1014 : INTEGER(KIND=int_8) :: mem, number_of_elements_per_blk
1015 64 : INTEGER, ALLOCATABLE, DIMENSION(:) :: procs_recv
1016 64 : INTEGER, DIMENSION(:, :), POINTER :: blacs2mpi
1017 : REAL(KIND=dp) :: mem_per_block, mem_real
1018 64 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), TARGET :: buffer_compens_1D
1019 64 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: mat_S_3D
1020 64 : TYPE(cp_1d_r_cp_type), ALLOCATABLE, DIMENSION(:) :: buffer_1D
1021 64 : TYPE(cp_3d_r_cp_type), ALLOCATABLE, DIMENSION(:) :: buffer_3D
1022 : TYPE(mp_para_env_type), POINTER :: para_env
1023 64 : TYPE(mp_request_type), ALLOCATABLE, DIMENSION(:) :: recv_requests, send_requests
1024 :
1025 64 : CALL timeset(routineN, handle)
1026 :
1027 : ! We allocate it at every step to reduce potential memory conflicts with COSMA
1028 64 : IF (mp2_env%ri_grad%dot_blksize >= blksize_threshold) THEN
1029 48 : CALL mp2_env%local_gemm_ctx%create(LOCAL_GEMM_PU_GPU)
1030 48 : CALL mp2_env%local_gemm_ctx%set_op_threshold_gpu(spla_threshold)
1031 : END IF
1032 :
1033 64 : tag = 47
1034 :
1035 64 : my_P_size = SIZE(mat_work_iaP_3D, 1)
1036 :
1037 64 : CALL cp_fm_struct_get(fm_struct_S, para_env=para_env)
1038 64 : CALL fm_struct_S%context%get(my_process_row=my_prow, blacs2mpi=blacs2mpi, para_env=para_env)
1039 :
1040 64 : CALL get_group_dist(gd_virtual, mepos(1), my_a_start, my_a_end, my_a_size)
1041 64 : CALL get_group_dist(gd_homo, mepos(2), my_i_start, my_i_end, my_i_size)
1042 :
1043 : ! We have to remap the indices because mp_sendrecv requires a 3D array (because of mat_work_iaP_3D)
1044 : ! and dgemm requires 2D arrays
1045 : ! Fortran 2008 does allow pointer remapping independently of the ranks but GCC 7 does not properly support it
1046 64 : mat_S_3D(1:my_P_size, 1:my_a_size, 1:my_i_size) => mat_S_1D(1:INT(my_P_size, int_8)*my_a_size*my_i_size)
1047 :
1048 : number_of_elements_per_blk = MAX(INT(maxsize(gd_homo), KIND=int_8)*my_a_size, &
1049 64 : INT(maxsize(gd_virtual), KIND=int_8)*my_i_size)*my_P_size
1050 :
1051 : ! Determine the available memory and estimate the number of possible parallel communication channels
1052 64 : CALL m_memory(mem)
1053 64 : mem_real = REAL(mem, KIND=dp)
1054 64 : mem_per_block = REAL(number_of_elements_per_blk, KIND=dp)*8.0_dp
1055 192 : number_of_parallel_channels = MAX(1, MIN(MAXVAL(grid) - 1, FLOOR(mem_real/mem_per_block)))
1056 64 : CALL para_env%min(number_of_parallel_channels)
1057 64 : IF (mp2_env%ri_grad%max_parallel_comm > 0) THEN
1058 64 : number_of_parallel_channels = MIN(number_of_parallel_channels, mp2_env%ri_grad%max_parallel_comm)
1059 : END IF
1060 :
1061 64 : IF (unit_nr > 0) THEN
1062 32 : WRITE (unit_nr, '(T3,A,T75,I6)') 'GRAD_INFO| Number of parallel communication channels:', number_of_parallel_channels
1063 32 : CALL m_flush(unit_nr)
1064 : END IF
1065 64 : CALL para_env%sync()
1066 :
1067 256 : ALLOCATE (buffer_1D(number_of_parallel_channels))
1068 128 : DO proc_shift = 1, number_of_parallel_channels
1069 256 : ALLOCATE (buffer_1D(proc_shift)%array(number_of_elements_per_blk))
1070 : END DO
1071 :
1072 256 : ALLOCATE (buffer_3D(number_of_parallel_channels))
1073 :
1074 : ! Allocate buffers for vector version of kahan summation
1075 64 : IF (mp2_env%ri_grad%dot_blksize >= blksize_threshold) THEN
1076 144 : ALLOCATE (buffer_compens_1D(2*MAX(my_a_size*maxsize(gd_virtual), my_i_size*maxsize(gd_homo))))
1077 : END IF
1078 :
1079 64 : IF (number_of_parallel_channels > 1) THEN
1080 0 : ALLOCATE (procs_recv(number_of_parallel_channels))
1081 0 : ALLOCATE (recv_requests(number_of_parallel_channels))
1082 0 : ALLOCATE (send_requests(MAXVAL(grid) - 1))
1083 : END IF
1084 :
1085 64 : IF (number_of_parallel_channels > 1 .AND. grid(1) > 1) THEN
1086 0 : CALL timeset(routineN//"_comm_a", handle2)
1087 0 : recv_requests(:) = mp_request_null
1088 0 : procs_recv(:) = -1
1089 0 : DO proc_shift = 1, MIN(grid(1) - 1, number_of_parallel_channels)
1090 0 : proc_a_recv = MODULO(mepos(1) - proc_shift, grid(1))
1091 0 : proc_recv = mepos(2)*grid(1) + proc_a_recv
1092 :
1093 0 : CALL get_group_dist(gd_virtual, proc_a_recv, recv_a_start, recv_a_end, recv_a_size)
1094 :
1095 : buffer_3D(proc_shift)%array(1:my_P_size, 1:recv_a_size, 1:my_i_size) => &
1096 0 : buffer_1D(proc_shift)%array(1:INT(my_P_size, KIND=int_8)*recv_a_size*my_i_size)
1097 :
1098 : CALL para_env%irecv(buffer_3D(proc_shift)%array, blacs2mpi(my_prow, proc_recv), &
1099 0 : recv_requests(proc_shift), tag)
1100 :
1101 0 : procs_recv(proc_shift) = proc_a_recv
1102 : END DO
1103 :
1104 0 : send_requests(:) = mp_request_null
1105 0 : DO proc_shift = 1, grid(1) - 1
1106 0 : proc_a_send = MODULO(mepos(1) + proc_shift, grid(1))
1107 0 : proc_send = mepos(2)*grid(1) + proc_a_send
1108 :
1109 : CALL para_env%isend(mat_work_iaP_3D, blacs2mpi(my_prow, proc_send), &
1110 0 : send_requests(proc_shift), tag)
1111 : END DO
1112 0 : CALL timestop(handle2)
1113 : END IF
1114 :
1115 : CALL calc_P_rpa_a(P_ab(:, my_a_start:my_a_end), &
1116 : mat_S_3D, mat_work_iaP_3D, &
1117 : mp2_env%ri_grad%dot_blksize, buffer_compens_1D, mp2_env%local_gemm_ctx, &
1118 : Eigenval(homo + my_a_start:homo + my_a_end), Eigenval(my_i_start:my_i_end), &
1119 64 : Eigenval(homo + my_a_start:homo + my_a_end), omega, weight)
1120 :
1121 64 : DO proc_shift = 1, grid(1) - 1
1122 0 : CALL timeset(routineN//"_comm_a", handle2)
1123 0 : IF (number_of_parallel_channels > 1) THEN
1124 0 : CALL mp_waitany(recv_requests, completed)
1125 :
1126 0 : CALL get_group_dist(gd_virtual, procs_recv(completed), recv_a_start, recv_a_end, recv_a_size)
1127 : ELSE
1128 0 : proc_a_send = MODULO(mepos(1) + proc_shift, grid(1))
1129 0 : proc_a_recv = MODULO(mepos(1) - proc_shift, grid(1))
1130 :
1131 0 : proc_send = mepos(2)*grid(1) + proc_a_send
1132 0 : proc_recv = mepos(2)*grid(1) + proc_a_recv
1133 :
1134 0 : CALL get_group_dist(gd_virtual, proc_a_recv, recv_a_start, recv_a_end, recv_a_size)
1135 :
1136 : buffer_3D(1)%array(1:my_P_size, 1:recv_a_size, 1:my_i_size) => &
1137 0 : buffer_1D(1)%array(1:INT(my_P_size, KIND=int_8)*recv_a_size*my_i_size)
1138 :
1139 : CALL para_env%sendrecv(mat_work_iaP_3D, blacs2mpi(my_prow, proc_send), &
1140 0 : buffer_3D(1)%array, blacs2mpi(my_prow, proc_recv), tag)
1141 0 : completed = 1
1142 : END IF
1143 0 : CALL timestop(handle2)
1144 :
1145 : CALL calc_P_rpa_a(P_ab(:, recv_a_start:recv_a_end), &
1146 : mat_S_3D, buffer_3D(completed)%array, &
1147 : mp2_env%ri_grad%dot_blksize, buffer_compens_1D, mp2_env%local_gemm_ctx, &
1148 : Eigenval(homo + my_a_start:homo + my_a_end), Eigenval(my_i_start:my_i_end), &
1149 0 : Eigenval(homo + recv_a_start:homo + recv_a_end), omega, weight)
1150 :
1151 64 : IF (number_of_parallel_channels > 1 .AND. number_of_parallel_channels + proc_shift < grid(1)) THEN
1152 0 : proc_a_recv = MODULO(mepos(1) - proc_shift - number_of_parallel_channels, grid(1))
1153 0 : proc_recv = mepos(2)*grid(1) + proc_a_recv
1154 :
1155 0 : CALL get_group_dist(gd_virtual, proc_a_recv, recv_a_start, recv_a_end, recv_a_size)
1156 :
1157 : buffer_3D(completed)%array(1:my_P_size, 1:recv_a_size, 1:my_i_size) => &
1158 0 : buffer_1D(completed)%array(1:INT(my_P_size, KIND=int_8)*recv_a_size*my_i_size)
1159 :
1160 : CALL para_env%irecv(buffer_3D(completed)%array, blacs2mpi(my_prow, proc_recv), &
1161 0 : recv_requests(completed), tag)
1162 :
1163 0 : procs_recv(completed) = proc_a_recv
1164 : END IF
1165 : END DO
1166 :
1167 64 : IF (number_of_parallel_channels > 1 .AND. grid(1) > 1) THEN
1168 0 : CALL mp_waitall(send_requests)
1169 : END IF
1170 :
1171 64 : IF (number_of_parallel_channels > 1 .AND. grid(2) > 1) THEN
1172 0 : recv_requests(:) = mp_request_null
1173 0 : procs_recv(:) = -1
1174 0 : DO proc_shift = 1, MIN(grid(2) - 1, number_of_parallel_channels)
1175 0 : proc_i_recv = MODULO(mepos(2) - proc_shift, grid(2))
1176 0 : proc_recv = proc_i_recv*grid(1) + mepos(1)
1177 :
1178 0 : CALL get_group_dist(gd_homo, proc_i_recv, recv_i_start, recv_i_end, recv_i_size)
1179 :
1180 : buffer_3D(proc_shift)%array(1:my_P_size, 1:my_a_size, 1:recv_i_size) => &
1181 0 : buffer_1D(proc_shift)%array(1:INT(my_P_size, KIND=int_8)*my_a_size*recv_i_size)
1182 :
1183 : CALL para_env%irecv(buffer_3D(proc_shift)%array, blacs2mpi(my_prow, proc_recv), &
1184 0 : recv_requests(proc_shift), tag)
1185 :
1186 0 : procs_recv(proc_shift) = proc_i_recv
1187 : END DO
1188 :
1189 0 : send_requests(:) = mp_request_null
1190 0 : DO proc_shift = 1, grid(2) - 1
1191 0 : proc_i_send = MODULO(mepos(2) + proc_shift, grid(2))
1192 0 : proc_send = proc_i_send*grid(1) + mepos(1)
1193 :
1194 : CALL para_env%isend(mat_work_iaP_3D, blacs2mpi(my_prow, proc_send), &
1195 0 : send_requests(proc_shift), tag)
1196 : END DO
1197 : END IF
1198 :
1199 : CALL calc_P_rpa_i(P_ij(:, my_i_start:my_i_end), &
1200 : mat_S_3D, mat_work_iaP_3D, &
1201 : mp2_env%ri_grad%dot_blksize, buffer_compens_1D, mp2_env%local_gemm_ctx, &
1202 : Eigenval(homo + my_a_start:homo + my_a_end), Eigenval(my_i_start:my_i_end), &
1203 64 : Eigenval(my_i_start:my_i_end), omega, weight)
1204 :
1205 64 : DO proc_shift = 1, grid(2) - 1
1206 0 : CALL timeset(routineN//"_comm_i", handle2)
1207 0 : IF (number_of_parallel_channels > 1) THEN
1208 0 : CALL mp_waitany(recv_requests, completed)
1209 :
1210 0 : CALL get_group_dist(gd_homo, procs_recv(completed), recv_i_start, recv_i_end, recv_i_size)
1211 : ELSE
1212 0 : proc_i_send = MODULO(mepos(2) + proc_shift, grid(2))
1213 0 : proc_i_recv = MODULO(mepos(2) - proc_shift, grid(2))
1214 :
1215 0 : proc_send = proc_i_send*grid(1) + mepos(1)
1216 0 : proc_recv = proc_i_recv*grid(1) + mepos(1)
1217 :
1218 0 : CALL get_group_dist(gd_homo, proc_i_recv, recv_i_start, recv_i_end, recv_i_size)
1219 :
1220 : buffer_3D(1)%array(1:my_P_size, 1:my_a_size, 1:recv_i_size) => &
1221 0 : buffer_1D(1)%array(1:INT(my_P_size, KIND=int_8)*my_a_size*recv_i_size)
1222 :
1223 : CALL para_env%sendrecv(mat_work_iaP_3D, blacs2mpi(my_prow, proc_send), &
1224 0 : buffer_3D(1)%array, blacs2mpi(my_prow, proc_recv), tag)
1225 0 : completed = 1
1226 : END IF
1227 0 : CALL timestop(handle2)
1228 :
1229 : CALL calc_P_rpa_i(P_ij(:, recv_i_start:recv_i_end), &
1230 : mat_S_3D, buffer_3D(completed)%array, &
1231 : mp2_env%ri_grad%dot_blksize, buffer_compens_1D, mp2_env%local_gemm_ctx, &
1232 : Eigenval(homo + my_a_start:homo + my_a_end), Eigenval(my_i_start:my_i_end), &
1233 0 : Eigenval(recv_i_start:recv_i_end), omega, weight)
1234 :
1235 64 : IF (number_of_parallel_channels > 1 .AND. number_of_parallel_channels + proc_shift < grid(2)) THEN
1236 0 : proc_i_recv = MODULO(mepos(2) - proc_shift - number_of_parallel_channels, grid(2))
1237 0 : proc_recv = proc_i_recv*grid(1) + mepos(1)
1238 :
1239 0 : CALL get_group_dist(gd_homo, proc_i_recv, recv_i_start, recv_a_end, recv_i_size)
1240 :
1241 : buffer_3D(completed)%array(1:my_P_size, 1:my_a_size, 1:recv_i_size) => &
1242 0 : buffer_1D(completed)%array(1:INT(my_P_size, KIND=int_8)*my_a_size*recv_i_size)
1243 :
1244 : CALL para_env%irecv(buffer_3D(completed)%array, blacs2mpi(my_prow, proc_recv), &
1245 0 : recv_requests(completed), tag)
1246 :
1247 0 : procs_recv(completed) = proc_i_recv
1248 : END IF
1249 : END DO
1250 :
1251 64 : IF (number_of_parallel_channels > 1 .AND. grid(2) > 1) THEN
1252 0 : CALL mp_waitall(send_requests)
1253 : END IF
1254 :
1255 64 : IF (number_of_parallel_channels > 1) THEN
1256 0 : DEALLOCATE (procs_recv)
1257 0 : DEALLOCATE (recv_requests)
1258 0 : DEALLOCATE (send_requests)
1259 : END IF
1260 :
1261 64 : IF (mp2_env%ri_grad%dot_blksize >= blksize_threshold) THEN
1262 : ! release memory allocated by local_gemm when run on GPU. local_gemm_ctx is null on cpu only runs
1263 48 : CALL mp2_env%local_gemm_ctx%destroy()
1264 48 : DEALLOCATE (buffer_compens_1D)
1265 : END IF
1266 :
1267 128 : DO proc_shift = 1, number_of_parallel_channels
1268 64 : NULLIFY (buffer_3D(proc_shift)%array)
1269 128 : DEALLOCATE (buffer_1D(proc_shift)%array)
1270 : END DO
1271 64 : DEALLOCATE (buffer_3D, buffer_1D)
1272 :
1273 64 : CALL timestop(handle)
1274 :
1275 192 : END SUBROUTINE calc_P_rpa
1276 :
1277 : ! **************************************************************************************************
1278 : !> \brief ...
1279 : !> \param P_ab ...
1280 : !> \param mat_S ...
1281 : !> \param mat_work ...
1282 : !> \param dot_blksize ...
1283 : !> \param buffer_1D ...
1284 : !> \param local_gemm_ctx ...
1285 : !> \param my_eval_virt ...
1286 : !> \param my_eval_occ ...
1287 : !> \param recv_eval_virt ...
1288 : !> \param omega ...
1289 : !> \param weight ...
1290 : ! **************************************************************************************************
1291 64 : SUBROUTINE calc_P_rpa_a(P_ab, mat_S, mat_work, dot_blksize, buffer_1D, local_gemm_ctx, &
1292 64 : my_eval_virt, my_eval_occ, recv_eval_virt, omega, weight)
1293 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: P_ab
1294 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: mat_S, mat_work
1295 : INTEGER, INTENT(IN) :: dot_blksize
1296 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
1297 : INTENT(INOUT), TARGET :: buffer_1D
1298 : TYPE(local_gemm_ctxt_type), INTENT(INOUT) :: local_gemm_ctx
1299 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: my_eval_virt, my_eval_occ, recv_eval_virt
1300 : REAL(KIND=dp), INTENT(IN) :: omega, weight
1301 :
1302 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_P_rpa_a'
1303 :
1304 : INTEGER :: handle, my_a, my_a_size, my_i, &
1305 : my_i_size, my_P_size, P_end, P_start, &
1306 : recv_a_size, stripesize
1307 64 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: buffer_compens, buffer_unscaled
1308 :
1309 64 : CALL timeset(routineN, handle)
1310 :
1311 64 : my_i_size = SIZE(mat_S, 3)
1312 64 : recv_a_size = SIZE(mat_work, 2)
1313 64 : my_a_size = SIZE(mat_S, 2)
1314 64 : my_P_size = SIZE(mat_S, 1)
1315 :
1316 64 : IF (dot_blksize >= blksize_threshold) THEN
1317 48 : buffer_compens(1:my_a_size, 1:recv_a_size) => buffer_1D(1:my_a_size*recv_a_size)
1318 22208 : buffer_compens = 0.0_dp
1319 48 : buffer_unscaled(1:my_a_size, 1:recv_a_size) => buffer_1D(my_a_size*recv_a_size + 1:2*my_a_size*recv_a_size)
1320 :
1321 : ! This loop imitates the actual tensor contraction
1322 232 : DO my_i = 1, my_i_size
1323 416 : DO P_start = 1, my_P_size, dot_blksize
1324 184 : stripesize = MIN(dot_blksize, my_P_size - P_start + 1)
1325 184 : P_end = P_start + stripesize - 1
1326 :
1327 : CALL local_gemm_ctx%gemm("T", "N", my_a_size, recv_a_size, stripesize, &
1328 : -weight, mat_S(P_start:P_end, :, my_i), stripesize, &
1329 : mat_work(P_start:P_end, :, my_i), stripesize, &
1330 184 : 0.0_dp, buffer_unscaled, my_a_size)
1331 :
1332 : CALL scale_buffer_and_add_compens_virt(buffer_unscaled, buffer_compens, omega, &
1333 184 : my_eval_virt, recv_eval_virt, my_eval_occ(my_i))
1334 :
1335 368 : CALL kahan_step(buffer_compens, P_ab)
1336 : END DO
1337 : END DO
1338 : ELSE
1339 : BLOCK
1340 : INTEGER :: recv_a
1341 : REAL(KIND=dp) :: tmp, e_i, e_a, e_b, omega2, my_compens, my_p, s
1342 16 : omega2 = -omega**2
1343 : !$OMP PARALLEL DO COLLAPSE(2) DEFAULT(NONE)&
1344 : !$OMP SHARED(my_a_size,recv_a_size,my_i_size,mat_S,my_eval_virt,recv_eval_virt,my_eval_occ,omega2,&
1345 : !$OMP P_ab,weight,mat_work)&
1346 16 : !$OMP PRIVATE(tmp,my_a,recv_a,my_i,e_a,e_b,e_i,my_compens,my_p,s)
1347 : DO my_a = 1, my_a_size
1348 : DO recv_a = 1, recv_a_size
1349 : e_a = my_eval_virt(my_a)
1350 : e_b = recv_eval_virt(recv_a)
1351 : my_p = P_ab(my_a, recv_a)
1352 : my_compens = 0.0_dp
1353 : DO my_i = 1, my_i_size
1354 : e_i = -my_eval_occ(my_i)
1355 : tmp = -weight*accurate_dot_product(mat_S(:, my_a, my_i), mat_work(:, recv_a, my_i)) &
1356 : *(1.0_dp + omega2/((e_a + e_i)*(e_b + e_i))) - my_compens
1357 : s = my_p + tmp
1358 : my_compens = (s - my_p) - tmp
1359 : my_p = s
1360 : END DO
1361 : P_ab(my_a, recv_a) = my_p
1362 : END DO
1363 : END DO
1364 : END BLOCK
1365 : END IF
1366 :
1367 64 : CALL timestop(handle)
1368 :
1369 64 : END SUBROUTINE calc_P_rpa_a
1370 :
1371 : ! **************************************************************************************************
1372 : !> \brief ...
1373 : !> \param P_ij ...
1374 : !> \param mat_S ...
1375 : !> \param mat_work ...
1376 : !> \param dot_blksize ...
1377 : !> \param buffer_1D ...
1378 : !> \param local_gemm_ctx ...
1379 : !> \param my_eval_virt ...
1380 : !> \param my_eval_occ ...
1381 : !> \param recv_eval_occ ...
1382 : !> \param omega ...
1383 : !> \param weight ...
1384 : ! **************************************************************************************************
1385 64 : SUBROUTINE calc_P_rpa_i(P_ij, mat_S, mat_work, dot_blksize, buffer_1D, local_gemm_ctx, &
1386 64 : my_eval_virt, my_eval_occ, recv_eval_occ, omega, weight)
1387 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: P_ij
1388 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: mat_S, mat_work
1389 : INTEGER, INTENT(IN) :: dot_blksize
1390 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
1391 : INTENT(INOUT), TARGET :: buffer_1D
1392 : TYPE(local_gemm_ctxt_type), INTENT(INOUT) :: local_gemm_ctx
1393 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: my_eval_virt, my_eval_occ, recv_eval_occ
1394 : REAL(KIND=dp), INTENT(IN) :: omega, weight
1395 :
1396 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_P_rpa_i'
1397 :
1398 : INTEGER :: handle, my_a, my_a_size, my_i, &
1399 : my_i_size, my_P_size, P_end, P_start, &
1400 : recv_i_size, stripesize
1401 64 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: buffer_compens, buffer_unscaled
1402 :
1403 64 : CALL timeset(routineN, handle)
1404 :
1405 64 : my_i_size = SIZE(mat_S, 3)
1406 64 : recv_i_size = SIZE(mat_work, 3)
1407 64 : my_a_size = SIZE(mat_S, 2)
1408 64 : my_P_size = SIZE(mat_S, 1)
1409 :
1410 64 : IF (dot_blksize >= blksize_threshold) THEN
1411 48 : buffer_compens(1:my_i_size, 1:recv_i_size) => buffer_1D(1:my_i_size*recv_i_size)
1412 944 : buffer_compens = 0.0_dp
1413 48 : buffer_unscaled(1:my_i_size, 1:recv_i_size) => buffer_1D(my_i_size*recv_i_size + 1:2*my_i_size*recv_i_size)
1414 :
1415 : ! This loop imitates the actual tensor contraction
1416 1048 : DO my_a = 1, my_a_size
1417 2048 : DO P_start = 1, my_P_size, dot_blksize
1418 1000 : stripesize = MIN(dot_blksize, my_P_size - P_start + 1)
1419 1000 : P_end = P_start + stripesize - 1
1420 :
1421 : CALL local_gemm_ctx%gemm("T", "N", my_i_size, recv_i_size, stripesize, &
1422 : weight, mat_S(P_start:P_end, my_a, :), stripesize, &
1423 : mat_work(P_start:P_end, my_a, :), stripesize, &
1424 1000 : 0.0_dp, buffer_unscaled, my_i_size)
1425 :
1426 : CALL scale_buffer_and_add_compens_occ(buffer_unscaled, buffer_compens, omega, &
1427 1000 : my_eval_occ, recv_eval_occ, my_eval_virt(my_a))
1428 :
1429 2000 : CALL kahan_step(buffer_compens, P_ij)
1430 : END DO
1431 : END DO
1432 : ELSE
1433 : BLOCK
1434 : REAL(KIND=dp) :: tmp, e_i, e_a, e_j, omega2, my_compens, my_p, s
1435 : INTEGER :: recv_i
1436 16 : omega2 = -omega**2
1437 : !$OMP PARALLEL DO COLLAPSE(2) DEFAULT(NONE)&
1438 : !$OMP SHARED(my_a_size,recv_i_size,my_i_size,mat_S,my_eval_occ,my_eval_virt,omega2,&
1439 : !$OMP recv_eval_occ,P_ij,weight,mat_work)&
1440 16 : !$OMP PRIVATE(tmp,my_a,recv_i,my_i,e_i,e_j,e_a,my_compens,my_p,s)
1441 : DO my_i = 1, my_i_size
1442 : DO recv_i = 1, recv_i_size
1443 : e_i = my_eval_occ(my_i)
1444 : e_j = recv_eval_occ(recv_i)
1445 : my_p = P_ij(my_i, recv_i)
1446 : my_compens = 0.0_dp
1447 : DO my_a = 1, my_a_size
1448 : e_a = my_eval_virt(my_a)
1449 : tmp = weight*accurate_dot_product(mat_S(:, my_a, my_i), mat_work(:, my_a, recv_i)) &
1450 : *(1.0_dp + omega2/((e_a - e_i)*(e_a - e_j))) - my_compens
1451 : s = my_p + tmp
1452 : my_compens = (s - my_p) - tmp
1453 : my_p = s
1454 : END DO
1455 : P_ij(my_i, recv_i) = my_p
1456 : END DO
1457 : END DO
1458 : END BLOCK
1459 : END IF
1460 :
1461 64 : CALL timestop(handle)
1462 :
1463 64 : END SUBROUTINE calc_P_rpa_i
1464 :
1465 : ! **************************************************************************************************
1466 : !> \brief ...
1467 : !> \param compens ...
1468 : !> \param P ...
1469 : ! **************************************************************************************************
1470 1184 : SUBROUTINE kahan_step(compens, P)
1471 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: compens, P
1472 :
1473 : CHARACTER(LEN=*), PARAMETER :: routineN = 'kahan_step'
1474 :
1475 : INTEGER :: handle, i, j
1476 : REAL(KIND=dp) :: my_compens, my_p, s
1477 :
1478 1184 : CALL timeset(routineN, handle)
1479 :
1480 1184 : !$OMP PARALLEL DO DEFAULT(NONE) SHARED(P,compens) PRIVATE(i,my_p,my_compens,s, j) COLLAPSE(2)
1481 : DO j = 1, SIZE(compens, 2)
1482 : DO i = 1, SIZE(compens, 1)
1483 : my_p = P(i, j)
1484 : my_compens = compens(i, j)
1485 : s = my_p + my_compens
1486 : compens(i, j) = (s - my_p) - my_compens
1487 : P(i, j) = s
1488 : END DO
1489 : END DO
1490 : !$OMP END PARALLEL DO
1491 :
1492 1184 : CALL timestop(handle)
1493 :
1494 1184 : END SUBROUTINE kahan_step
1495 :
1496 : ! **************************************************************************************************
1497 : !> \brief ...
1498 : !> \param buffer_unscaled ...
1499 : !> \param buffer_compens ...
1500 : !> \param omega ...
1501 : !> \param my_eval_virt ...
1502 : !> \param recv_eval_virt ...
1503 : !> \param my_eval_occ ...
1504 : ! **************************************************************************************************
1505 184 : SUBROUTINE scale_buffer_and_add_compens_virt(buffer_unscaled, buffer_compens, omega, &
1506 184 : my_eval_virt, recv_eval_virt, my_eval_occ)
1507 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: buffer_unscaled
1508 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: buffer_compens
1509 : REAL(KIND=dp), INTENT(IN) :: omega
1510 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: my_eval_virt, recv_eval_virt
1511 : REAL(KIND=dp), INTENT(IN) :: my_eval_occ
1512 :
1513 : CHARACTER(LEN=*), PARAMETER :: routineN = 'scale_buffer_and_add_compens_virt'
1514 :
1515 : INTEGER :: handle, my_a, my_b
1516 :
1517 184 : CALL timeset(routineN, handle)
1518 :
1519 : !$OMP PARALLEL DO DEFAULT(NONE) SHARED(buffer_unscaled,buffer_compens,omega,&
1520 184 : !$OMP my_eval_virt,recv_eval_virt,my_eval_occ) PRIVATE(my_a,my_b)
1521 : DO my_b = 1, SIZE(buffer_compens, 2)
1522 : DO my_a = 1, SIZE(buffer_compens, 1)
1523 : buffer_compens(my_a, my_b) = buffer_unscaled(my_a, my_b) &
1524 : *(1.0_dp - omega**2/((my_eval_virt(my_a) - my_eval_occ)*(recv_eval_virt(my_b) - my_eval_occ))) &
1525 : - buffer_compens(my_a, my_b)
1526 : END DO
1527 : END DO
1528 : !$OMP END PARALLEL DO
1529 :
1530 184 : CALL timestop(handle)
1531 :
1532 184 : END SUBROUTINE scale_buffer_and_add_compens_virt
1533 :
1534 : ! **************************************************************************************************
1535 : !> \brief ...
1536 : !> \param buffer_unscaled ...
1537 : !> \param buffer_compens ...
1538 : !> \param omega ...
1539 : !> \param my_eval_occ ...
1540 : !> \param recv_eval_occ ...
1541 : !> \param my_eval_virt ...
1542 : ! **************************************************************************************************
1543 1000 : SUBROUTINE scale_buffer_and_add_compens_occ(buffer_unscaled, buffer_compens, omega, &
1544 1000 : my_eval_occ, recv_eval_occ, my_eval_virt)
1545 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: buffer_unscaled
1546 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: buffer_compens
1547 : REAL(KIND=dp), INTENT(IN) :: omega
1548 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: my_eval_occ, recv_eval_occ
1549 : REAL(KIND=dp), INTENT(IN) :: my_eval_virt
1550 :
1551 : CHARACTER(LEN=*), PARAMETER :: routineN = 'scale_buffer_and_add_compens_occ'
1552 :
1553 : INTEGER :: handle, my_i, my_j
1554 :
1555 1000 : CALL timeset(routineN, handle)
1556 :
1557 : !$OMP PARALLEL DO DEFAULT(NONE) SHARED(buffer_compens,buffer_unscaled,omega,&
1558 1000 : !$OMP my_eval_virt,my_eval_occ,recv_eval_occ) PRIVATE(my_i,my_j)
1559 : DO my_j = 1, SIZE(buffer_compens, 2)
1560 : DO my_i = 1, SIZE(buffer_compens, 1)
1561 : buffer_compens(my_i, my_j) = buffer_unscaled(my_i, my_j) &
1562 : *(1.0_dp - omega**2/((my_eval_virt - my_eval_occ(my_i))*(my_eval_virt - recv_eval_occ(my_j)))) &
1563 : - buffer_compens(my_i, my_j)
1564 : END DO
1565 : END DO
1566 : !$OMP END PARALLEL DO
1567 :
1568 1000 : CALL timestop(handle)
1569 :
1570 1000 : END SUBROUTINE scale_buffer_and_add_compens_occ
1571 :
1572 : ! **************************************************************************************************
1573 : !> \brief ...
1574 : !> \param x ...
1575 : !> \return ...
1576 : ! **************************************************************************************************
1577 1320 : ELEMENTAL FUNCTION sinh_over_x(x) RESULT(res)
1578 : REAL(KIND=dp), INTENT(IN) :: x
1579 : REAL(KIND=dp) :: res
1580 :
1581 : ! Calculate sinh(x)/x
1582 : ! Split the intervall to prevent numerical instabilities
1583 1320 : IF (ABS(x) > 3.0e-4_dp) THEN
1584 1318 : res = SINH(x)/x
1585 : ELSE
1586 2 : res = 1.0_dp + x**2/6.0_dp
1587 : END IF
1588 :
1589 1320 : END FUNCTION sinh_over_x
1590 :
1591 : ! **************************************************************************************************
1592 : !> \brief ...
1593 : !> \param fm_work_iaP ...
1594 : !> \param fm_mat_S ...
1595 : !> \param pair_list ...
1596 : !> \param virtual ...
1597 : !> \param P_ij ...
1598 : !> \param Eigenval ...
1599 : !> \param omega ...
1600 : !> \param weight ...
1601 : !> \param index2send ...
1602 : !> \param index2recv ...
1603 : !> \param dot_blksize ...
1604 : ! **************************************************************************************************
1605 8 : SUBROUTINE calc_Pij_degen(fm_work_iaP, fm_mat_S, pair_list, virtual, P_ij, Eigenval, &
1606 8 : omega, weight, index2send, index2recv, dot_blksize)
1607 : TYPE(cp_fm_type), INTENT(IN) :: fm_work_iaP, fm_mat_S
1608 : INTEGER, DIMENSION(:, :), INTENT(IN) :: pair_list
1609 : INTEGER, INTENT(IN) :: virtual
1610 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: P_ij
1611 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval
1612 : REAL(KIND=dp), INTENT(IN) :: omega, weight
1613 : TYPE(one_dim_int_array), DIMENSION(0:), INTENT(IN) :: index2send, index2recv
1614 : INTEGER, INTENT(IN) :: dot_blksize
1615 :
1616 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_Pij_degen'
1617 :
1618 : INTEGER :: avirt, col_global, col_local, counter, handle, handle2, ij_counter, iocc, &
1619 : my_col_local, my_i, my_j, my_pcol, my_prow, ncol_local, nrow_local, num_ij_pairs, &
1620 : num_pe_col, pcol, pcol_recv, pcol_send, proc_shift, recv_size, send_size, &
1621 : size_recv_buffer, size_send_buffer, tag
1622 8 : INTEGER, DIMENSION(:), POINTER :: col_indices, ncol_locals
1623 8 : INTEGER, DIMENSION(:, :), POINTER :: blacs2mpi
1624 : REAL(KIND=dp) :: trace
1625 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: buffer_recv, buffer_send
1626 : TYPE(cp_blacs_env_type), POINTER :: context
1627 : TYPE(mp_para_env_type), POINTER :: para_env
1628 :
1629 8 : CALL timeset(routineN, handle)
1630 :
1631 : CALL cp_fm_struct_get(fm_work_iaP%matrix_struct, para_env=para_env, ncol_locals=ncol_locals, &
1632 : ncol_local=ncol_local, col_indices=col_indices, &
1633 8 : context=context, nrow_local=nrow_local)
1634 : CALL context%get(my_process_row=my_prow, my_process_column=my_pcol, &
1635 8 : number_of_process_columns=num_pe_col, blacs2mpi=blacs2mpi)
1636 :
1637 8 : num_ij_pairs = SIZE(pair_list, 2)
1638 :
1639 8 : tag = 42
1640 :
1641 104 : DO ij_counter = 1, num_ij_pairs
1642 :
1643 96 : my_i = pair_list(1, ij_counter)
1644 96 : my_j = pair_list(2, ij_counter)
1645 :
1646 96 : trace = 0.0_dp
1647 :
1648 7392 : DO col_local = 1, ncol_local
1649 7296 : col_global = col_indices(col_local)
1650 :
1651 7296 : iocc = MAX(1, col_global - 1)/virtual + 1
1652 7296 : avirt = col_global - (iocc - 1)*virtual
1653 :
1654 7296 : IF (iocc /= my_j) CYCLE
1655 1824 : pcol = fm_work_iaP%matrix_struct%g2p_col((my_i - 1)*virtual + avirt)
1656 1824 : IF (pcol /= my_pcol) CYCLE
1657 :
1658 1824 : my_col_local = fm_work_iaP%matrix_struct%g2l_col((my_i - 1)*virtual + avirt)
1659 :
1660 : trace = trace + accurate_dot_product_2(fm_mat_S%local_data(:, my_col_local), fm_work_iaP%local_data(:, col_local), &
1661 7392 : dot_blksize)
1662 : END DO
1663 :
1664 104 : P_ij(ij_counter) = P_ij(ij_counter) - trace*sinh_over_x(0.5_dp*(Eigenval(my_i) - Eigenval(my_j))*omega)*omega*weight
1665 :
1666 : END DO
1667 :
1668 8 : IF (num_pe_col > 1) THEN
1669 : size_send_buffer = 0
1670 : size_recv_buffer = 0
1671 0 : DO proc_shift = 1, num_pe_col - 1
1672 0 : pcol_send = MODULO(my_pcol + proc_shift, num_pe_col)
1673 0 : pcol_recv = MODULO(my_pcol - proc_shift, num_pe_col)
1674 :
1675 0 : IF (ALLOCATED(index2send(pcol_send)%array)) THEN
1676 0 : size_send_buffer = MAX(size_send_buffer, SIZE(index2send(pcol_send)%array))
1677 : END IF
1678 :
1679 0 : IF (ALLOCATED(index2recv(pcol_recv)%array)) THEN
1680 0 : size_recv_buffer = MAX(size_recv_buffer, SIZE(index2recv(pcol_recv)%array))
1681 : END IF
1682 : END DO
1683 :
1684 0 : ALLOCATE (buffer_send(nrow_local, size_send_buffer), buffer_recv(nrow_local, size_recv_buffer))
1685 :
1686 0 : DO proc_shift = 1, num_pe_col - 1
1687 0 : pcol_send = MODULO(my_pcol + proc_shift, num_pe_col)
1688 0 : pcol_recv = MODULO(my_pcol - proc_shift, num_pe_col)
1689 :
1690 : ! Collect data and exchange
1691 0 : send_size = 0
1692 0 : IF (ALLOCATED(index2send(pcol_send)%array)) send_size = SIZE(index2send(pcol_send)%array)
1693 :
1694 0 : DO counter = 1, send_size
1695 0 : buffer_send(:, counter) = fm_work_iaP%local_data(:, index2send(pcol_send)%array(counter))
1696 : END DO
1697 :
1698 0 : recv_size = 0
1699 0 : IF (ALLOCATED(index2recv(pcol_recv)%array)) recv_size = SIZE(index2recv(pcol_recv)%array)
1700 0 : IF (recv_size > 0) THEN
1701 0 : CALL timeset(routineN//"_send", handle2)
1702 0 : IF (send_size > 0) THEN
1703 : CALL para_env%sendrecv(buffer_send(:, :send_size), blacs2mpi(my_prow, pcol_send), &
1704 0 : buffer_recv(:, :recv_size), blacs2mpi(my_prow, pcol_recv), tag)
1705 : ELSE
1706 0 : CALL para_env%recv(buffer_recv(:, :recv_size), blacs2mpi(my_prow, pcol_recv), tag)
1707 : END IF
1708 0 : CALL timestop(handle2)
1709 :
1710 0 : DO ij_counter = 1, num_ij_pairs
1711 : ! Collect the contributions of the matrix elements
1712 :
1713 0 : my_i = pair_list(1, ij_counter)
1714 0 : my_j = pair_list(2, ij_counter)
1715 :
1716 0 : trace = 0.0_dp
1717 :
1718 0 : DO col_local = 1, recv_size
1719 0 : col_global = index2recv(pcol_recv)%array(col_local)
1720 :
1721 0 : iocc = MAX(1, col_global - 1)/virtual + 1
1722 0 : IF (iocc /= my_j) CYCLE
1723 0 : avirt = col_global - (iocc - 1)*virtual
1724 0 : pcol = fm_work_iaP%matrix_struct%g2p_col((my_i - 1)*virtual + avirt)
1725 0 : IF (pcol /= my_pcol) CYCLE
1726 :
1727 0 : my_col_local = fm_work_iaP%matrix_struct%g2l_col((my_i - 1)*virtual + avirt)
1728 :
1729 : trace = trace + accurate_dot_product_2(fm_mat_S%local_data(:, my_col_local), buffer_recv(:, col_local), &
1730 0 : dot_blksize)
1731 : END DO
1732 :
1733 : P_ij(ij_counter) = P_ij(ij_counter) &
1734 0 : - trace*sinh_over_x(0.5_dp*(Eigenval(my_i) - Eigenval(my_j))*omega)*omega*weight
1735 : END DO
1736 0 : ELSE IF (send_size > 0) THEN
1737 0 : CALL timeset(routineN//"_send", handle2)
1738 0 : CALL para_env%send(buffer_send(:, :send_size), blacs2mpi(my_prow, pcol_send), tag)
1739 0 : CALL timestop(handle2)
1740 : END IF
1741 : END DO
1742 0 : IF (ALLOCATED(buffer_send)) DEALLOCATE (buffer_send)
1743 0 : IF (ALLOCATED(buffer_recv)) DEALLOCATE (buffer_recv)
1744 : END IF
1745 :
1746 8 : CALL timestop(handle)
1747 :
1748 16 : END SUBROUTINE calc_Pij_degen
1749 :
1750 : ! **************************************************************************************************
1751 : !> \brief ...
1752 : !> \param fm_work_iaP ...
1753 : !> \param fm_mat_S ...
1754 : !> \param pair_list ...
1755 : !> \param virtual ...
1756 : !> \param P_ab ...
1757 : !> \param Eigenval ...
1758 : !> \param omega ...
1759 : !> \param weight ...
1760 : !> \param index2send ...
1761 : !> \param index2recv ...
1762 : !> \param dot_blksize ...
1763 : ! **************************************************************************************************
1764 8 : SUBROUTINE calc_Pab_degen(fm_work_iaP, fm_mat_S, pair_list, virtual, P_ab, Eigenval, &
1765 8 : omega, weight, index2send, index2recv, dot_blksize)
1766 : TYPE(cp_fm_type), INTENT(IN) :: fm_work_iaP, fm_mat_S
1767 : INTEGER, DIMENSION(:, :), INTENT(IN) :: pair_list
1768 : INTEGER, INTENT(IN) :: virtual
1769 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: P_ab
1770 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval
1771 : REAL(KIND=dp), INTENT(IN) :: omega, weight
1772 : TYPE(one_dim_int_array), DIMENSION(0:), INTENT(IN) :: index2send, index2recv
1773 : INTEGER, INTENT(IN) :: dot_blksize
1774 :
1775 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_Pab_degen'
1776 :
1777 : INTEGER :: ab_counter, avirt, col_global, col_local, counter, handle, handle2, iocc, my_a, &
1778 : my_b, my_col_local, my_pcol, my_prow, ncol_local, nrow_local, num_ab_pairs, num_pe_col, &
1779 : pcol, pcol_recv, pcol_send, proc_shift, recv_size, send_size, size_recv_buffer, &
1780 : size_send_buffer, tag
1781 8 : INTEGER, DIMENSION(:), POINTER :: col_indices, ncol_locals
1782 8 : INTEGER, DIMENSION(:, :), POINTER :: blacs2mpi
1783 : REAL(KIND=dp) :: trace
1784 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: buffer_recv, buffer_send
1785 : TYPE(cp_blacs_env_type), POINTER :: context
1786 : TYPE(mp_para_env_type), POINTER :: para_env
1787 :
1788 8 : CALL timeset(routineN, handle)
1789 :
1790 : CALL cp_fm_struct_get(fm_work_iaP%matrix_struct, para_env=para_env, ncol_locals=ncol_locals, &
1791 : ncol_local=ncol_local, col_indices=col_indices, &
1792 8 : context=context, nrow_local=nrow_local)
1793 : CALL context%get(my_process_row=my_prow, my_process_column=my_pcol, &
1794 8 : number_of_process_columns=num_pe_col, blacs2mpi=blacs2mpi)
1795 :
1796 8 : num_ab_pairs = SIZE(pair_list, 2)
1797 :
1798 8 : tag = 43
1799 :
1800 1232 : DO ab_counter = 1, num_ab_pairs
1801 :
1802 1224 : my_a = pair_list(1, ab_counter)
1803 1224 : my_b = pair_list(2, ab_counter)
1804 :
1805 1224 : trace = 0.0_dp
1806 :
1807 94248 : DO col_local = 1, ncol_local
1808 93024 : col_global = col_indices(col_local)
1809 :
1810 93024 : iocc = MAX(1, col_global - 1)/virtual + 1
1811 93024 : avirt = col_global - (iocc - 1)*virtual
1812 :
1813 93024 : IF (avirt /= my_b) CYCLE
1814 4896 : pcol = fm_work_iaP%matrix_struct%g2p_col((iocc - 1)*virtual + my_a)
1815 4896 : IF (pcol /= my_pcol) CYCLE
1816 4896 : my_col_local = fm_work_iaP%matrix_struct%g2l_col((iocc - 1)*virtual + my_a)
1817 :
1818 : trace = trace + accurate_dot_product_2(fm_mat_S%local_data(:, my_col_local), fm_work_iaP%local_data(:, col_local), &
1819 94248 : dot_blksize)
1820 :
1821 : END DO
1822 :
1823 : P_ab(ab_counter) = P_ab(ab_counter) &
1824 1232 : + trace*sinh_over_x(0.5_dp*(Eigenval(my_a) - Eigenval(my_b))*omega)*omega*weight
1825 :
1826 : END DO
1827 :
1828 8 : IF (num_pe_col > 1) THEN
1829 : size_send_buffer = 0
1830 : size_recv_buffer = 0
1831 0 : DO proc_shift = 1, num_pe_col - 1
1832 0 : pcol_send = MODULO(my_pcol + proc_shift, num_pe_col)
1833 0 : pcol_recv = MODULO(my_pcol - proc_shift, num_pe_col)
1834 :
1835 0 : IF (ALLOCATED(index2send(pcol_send)%array)) THEN
1836 0 : size_send_buffer = MAX(size_send_buffer, SIZE(index2send(pcol_send)%array))
1837 : END IF
1838 :
1839 0 : IF (ALLOCATED(index2recv(pcol_recv)%array)) THEN
1840 0 : size_recv_buffer = MAX(size_recv_buffer, SIZE(index2recv(pcol_recv)%array))
1841 : END IF
1842 : END DO
1843 :
1844 0 : ALLOCATE (buffer_send(nrow_local, size_send_buffer), buffer_recv(nrow_local, size_recv_buffer))
1845 :
1846 0 : DO proc_shift = 1, num_pe_col - 1
1847 0 : pcol_send = MODULO(my_pcol + proc_shift, num_pe_col)
1848 0 : pcol_recv = MODULO(my_pcol - proc_shift, num_pe_col)
1849 :
1850 : ! Collect data and exchange
1851 0 : send_size = 0
1852 0 : IF (ALLOCATED(index2send(pcol_send)%array)) send_size = SIZE(index2send(pcol_send)%array)
1853 :
1854 0 : DO counter = 1, send_size
1855 0 : buffer_send(:, counter) = fm_work_iaP%local_data(:, index2send(pcol_send)%array(counter))
1856 : END DO
1857 :
1858 0 : recv_size = 0
1859 0 : IF (ALLOCATED(index2recv(pcol_recv)%array)) recv_size = SIZE(index2recv(pcol_recv)%array)
1860 0 : IF (recv_size > 0) THEN
1861 0 : CALL timeset(routineN//"_send", handle2)
1862 0 : IF (send_size > 0) THEN
1863 : CALL para_env%sendrecv(buffer_send(:, :send_size), blacs2mpi(my_prow, pcol_send), &
1864 0 : buffer_recv(:, :recv_size), blacs2mpi(my_prow, pcol_recv), tag)
1865 : ELSE
1866 0 : CALL para_env%recv(buffer_recv(:, :recv_size), blacs2mpi(my_prow, pcol_recv), tag)
1867 : END IF
1868 0 : CALL timestop(handle2)
1869 :
1870 0 : DO ab_counter = 1, num_ab_pairs
1871 : ! Collect the contributions of the matrix elements
1872 :
1873 0 : my_a = pair_list(1, ab_counter)
1874 0 : my_b = pair_list(2, ab_counter)
1875 :
1876 0 : trace = 0.0_dp
1877 :
1878 0 : DO col_local = 1, SIZE(index2recv(pcol_recv)%array)
1879 0 : col_global = index2recv(pcol_recv)%array(col_local)
1880 :
1881 0 : iocc = MAX(1, col_global - 1)/virtual + 1
1882 0 : avirt = col_global - (iocc - 1)*virtual
1883 0 : IF (avirt /= my_b) CYCLE
1884 0 : pcol = fm_work_iaP%matrix_struct%g2p_col((iocc - 1)*virtual + my_a)
1885 0 : IF (pcol /= my_pcol) CYCLE
1886 :
1887 0 : my_col_local = fm_work_iaP%matrix_struct%g2l_col((iocc - 1)*virtual + my_a)
1888 :
1889 : trace = trace + accurate_dot_product_2(fm_mat_S%local_data(:, my_col_local), buffer_recv(:, col_local), &
1890 0 : dot_blksize)
1891 : END DO
1892 :
1893 : P_ab(ab_counter) = P_ab(ab_counter) &
1894 0 : + trace*sinh_over_x(0.5_dp*(Eigenval(my_a) - Eigenval(my_b))*omega)*omega*weight
1895 :
1896 : END DO
1897 0 : ELSE IF (send_size > 0) THEN
1898 0 : CALL timeset(routineN//"_send", handle2)
1899 0 : CALL para_env%send(buffer_send(:, :send_size), blacs2mpi(my_prow, pcol_send), tag)
1900 0 : CALL timestop(handle2)
1901 : END IF
1902 : END DO
1903 0 : IF (ALLOCATED(buffer_send)) DEALLOCATE (buffer_send)
1904 0 : IF (ALLOCATED(buffer_recv)) DEALLOCATE (buffer_recv)
1905 : END IF
1906 :
1907 8 : CALL timestop(handle)
1908 :
1909 16 : END SUBROUTINE calc_Pab_degen
1910 :
1911 : ! **************************************************************************************************
1912 : !> \brief ...
1913 : !> \param index2send ...
1914 : !> \param index2recv ...
1915 : !> \param fm_mat_S ...
1916 : !> \param mat_S_3D ...
1917 : !> \param gd_homo ...
1918 : !> \param gd_virtual ...
1919 : !> \param mepos ...
1920 : ! **************************************************************************************************
1921 128 : SUBROUTINE redistribute_fm_mat_S(index2send, index2recv, fm_mat_S, mat_S_3D, gd_homo, gd_virtual, mepos)
1922 : TYPE(one_dim_int_array), DIMENSION(0:), INTENT(IN) :: index2send
1923 : TYPE(two_dim_int_array), DIMENSION(0:), INTENT(IN) :: index2recv
1924 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S
1925 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
1926 : INTENT(OUT) :: mat_S_3D
1927 : TYPE(group_dist_d1_type), INTENT(IN) :: gd_homo, gd_virtual
1928 : INTEGER, DIMENSION(2), INTENT(IN) :: mepos
1929 :
1930 : CHARACTER(LEN=*), PARAMETER :: routineN = 'redistribute_fm_mat_S'
1931 :
1932 : INTEGER :: col_local, handle, my_a, my_homo, my_i, my_pcol, my_prow, my_virtual, nrow_local, &
1933 : num_pe_col, proc_recv, proc_send, proc_shift, recv_size, send_size, size_recv_buffer, &
1934 : size_send_buffer, tag
1935 128 : INTEGER, DIMENSION(:, :), POINTER :: blacs2mpi
1936 128 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: buffer_recv, buffer_send
1937 : TYPE(mp_para_env_type), POINTER :: para_env
1938 :
1939 128 : CALL timeset(routineN, handle)
1940 :
1941 128 : tag = 46
1942 :
1943 : CALL fm_mat_S%matrix_struct%context%get(my_process_row=my_prow, my_process_column=my_pcol, &
1944 128 : number_of_process_columns=num_pe_col, blacs2mpi=blacs2mpi)
1945 :
1946 128 : CALL cp_fm_struct_get(fm_mat_S%matrix_struct, nrow_local=nrow_local, para_env=para_env)
1947 :
1948 128 : CALL get_group_dist(gd_homo, mepos(2), sizes=my_homo)
1949 128 : CALL get_group_dist(gd_virtual, mepos(1), sizes=my_virtual)
1950 :
1951 640 : ALLOCATE (mat_S_3D(nrow_local, my_virtual, my_homo))
1952 :
1953 128 : IF (ALLOCATED(index2send(my_pcol)%array)) THEN
1954 10032 : DO col_local = 1, SIZE(index2send(my_pcol)%array)
1955 9904 : my_a = index2recv(my_pcol)%array(1, col_local)
1956 9904 : my_i = index2recv(my_pcol)%array(2, col_local)
1957 769440 : mat_S_3D(:, my_a, my_i) = fm_mat_S%local_data(:, index2send(my_pcol)%array(col_local))
1958 : END DO
1959 : END IF
1960 :
1961 128 : IF (num_pe_col > 1) THEN
1962 : size_send_buffer = 0
1963 : size_recv_buffer = 0
1964 0 : DO proc_shift = 1, num_pe_col - 1
1965 0 : proc_send = MODULO(my_pcol + proc_shift, num_pe_col)
1966 0 : proc_recv = MODULO(my_pcol - proc_shift, num_pe_col)
1967 :
1968 0 : send_size = 0
1969 0 : IF (ALLOCATED(index2send(proc_send)%array)) send_size = SIZE(index2send(proc_send)%array)
1970 0 : size_send_buffer = MAX(size_send_buffer, send_size)
1971 :
1972 0 : recv_size = 0
1973 0 : IF (ALLOCATED(index2recv(proc_recv)%array)) recv_size = SIZE(index2recv(proc_recv)%array)
1974 0 : size_recv_buffer = MAX(size_recv_buffer, recv_size)
1975 :
1976 : END DO
1977 :
1978 0 : ALLOCATE (buffer_send(nrow_local, size_send_buffer), buffer_recv(nrow_local, size_recv_buffer))
1979 :
1980 0 : DO proc_shift = 1, num_pe_col - 1
1981 0 : proc_send = MODULO(my_pcol + proc_shift, num_pe_col)
1982 0 : proc_recv = MODULO(my_pcol - proc_shift, num_pe_col)
1983 :
1984 0 : send_size = 0
1985 0 : IF (ALLOCATED(index2send(proc_send)%array)) send_size = SIZE(index2send(proc_send)%array)
1986 0 : DO col_local = 1, send_size
1987 0 : buffer_send(:, col_local) = fm_mat_S%local_data(:, index2send(proc_send)%array(col_local))
1988 : END DO
1989 :
1990 0 : recv_size = 0
1991 0 : IF (ALLOCATED(index2recv(proc_recv)%array)) recv_size = SIZE(index2recv(proc_recv)%array, 2)
1992 0 : IF (recv_size > 0) THEN
1993 0 : IF (send_size > 0) THEN
1994 : CALL para_env%sendrecv(buffer_send(:, :send_size), blacs2mpi(my_prow, proc_send), &
1995 0 : buffer_recv(:, :recv_size), blacs2mpi(my_prow, proc_recv), tag)
1996 : ELSE
1997 0 : CALL para_env%recv(buffer_recv(:, :recv_size), blacs2mpi(my_prow, proc_recv), tag)
1998 : END IF
1999 :
2000 0 : DO col_local = 1, recv_size
2001 0 : my_a = index2recv(proc_recv)%array(1, col_local)
2002 0 : my_i = index2recv(proc_recv)%array(2, col_local)
2003 0 : mat_S_3D(:, my_a, my_i) = buffer_recv(:, col_local)
2004 : END DO
2005 0 : ELSE IF (send_size > 0) THEN
2006 0 : CALL para_env%send(buffer_send(:, :send_size), blacs2mpi(my_prow, proc_send), tag)
2007 : END IF
2008 :
2009 : END DO
2010 :
2011 0 : IF (ALLOCATED(buffer_send)) DEALLOCATE (buffer_send)
2012 0 : IF (ALLOCATED(buffer_recv)) DEALLOCATE (buffer_recv)
2013 : END IF
2014 :
2015 128 : CALL timestop(handle)
2016 :
2017 384 : END SUBROUTINE redistribute_fm_mat_S
2018 :
2019 : ! **************************************************************************************************
2020 : !> \brief ...
2021 : !> \param rpa_grad ...
2022 : !> \param mp2_env ...
2023 : !> \param para_env_sub ...
2024 : !> \param para_env ...
2025 : !> \param qs_env ...
2026 : !> \param gd_array ...
2027 : !> \param color_sub ...
2028 : !> \param do_ri_sos_laplace_mp2 ...
2029 : !> \param homo ...
2030 : !> \param virtual ...
2031 : ! **************************************************************************************************
2032 44 : SUBROUTINE rpa_grad_finalize(rpa_grad, mp2_env, para_env_sub, para_env, qs_env, gd_array, &
2033 44 : color_sub, do_ri_sos_laplace_mp2, homo, virtual)
2034 : TYPE(rpa_grad_type), INTENT(INOUT) :: rpa_grad
2035 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
2036 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env_sub, para_env
2037 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
2038 : TYPE(group_dist_d1_type) :: gd_array
2039 : INTEGER, INTENT(IN) :: color_sub
2040 : LOGICAL, INTENT(IN) :: do_ri_sos_laplace_mp2
2041 : INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
2042 :
2043 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rpa_grad_finalize'
2044 :
2045 : INTEGER :: dimen_ia, dimen_RI, handle, iiB, ispin, my_group_L_end, my_group_L_size, &
2046 : my_group_L_start, my_ia_end, my_ia_size, my_ia_start, my_P_end, my_P_size, my_P_start, &
2047 : ngroup, nspins, pos_group, pos_sub, proc
2048 44 : INTEGER, ALLOCATABLE, DIMENSION(:) :: pos_info
2049 44 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: group_grid_2_mepos, mepos_2_grid_group
2050 : REAL(KIND=dp) :: my_scale
2051 44 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Gamma_2D
2052 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
2053 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
2054 : TYPE(cp_fm_type) :: fm_G_P_ia, fm_PQ, fm_PQ_2, fm_PQ_half, &
2055 : fm_work_PQ, fm_work_PQ_2, fm_Y, &
2056 : operator_half
2057 44 : TYPE(group_dist_d1_type) :: gd_array_new, gd_ia, gd_P, gd_P_new
2058 :
2059 44 : CALL timeset(routineN, handle)
2060 :
2061 : ! Release unnecessary matrices to save memory for next steps
2062 :
2063 44 : nspins = SIZE(rpa_grad%fm_Y)
2064 :
2065 : ! Scaling factor is required to scale the density matrices and the Gamma matrices later
2066 44 : IF (do_ri_sos_laplace_mp2) THEN
2067 20 : my_scale = mp2_env%scale_s
2068 : ELSE
2069 24 : my_scale = -mp2_env%ri_rpa%scale_rpa/(2.0_dp*pi)
2070 24 : IF (mp2_env%ri_rpa%minimax_quad) my_scale = my_scale/2.0_dp
2071 : END IF
2072 :
2073 44 : IF (do_ri_sos_laplace_mp2) THEN
2074 : CALL sos_mp2_grad_finalize(rpa_grad%sos_mp2_work_occ, rpa_grad%sos_mp2_work_virt, &
2075 20 : para_env, para_env_sub, homo, virtual, mp2_env)
2076 : ELSE
2077 : CALL rpa_grad_work_finalize(rpa_grad%rpa_work, mp2_env, homo, &
2078 24 : virtual, para_env, para_env_sub)
2079 : END IF
2080 :
2081 44 : CALL get_qs_env(qs_env, blacs_env=blacs_env)
2082 :
2083 44 : CALL cp_fm_get_info(rpa_grad%fm_Gamma_PQ, ncol_global=dimen_RI)
2084 :
2085 44 : NULLIFY (fm_struct)
2086 : CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=dimen_RI, &
2087 44 : ncol_global=dimen_RI, para_env=para_env)
2088 44 : CALL cp_fm_create(fm_PQ, fm_struct)
2089 44 : CALL cp_fm_create(fm_work_PQ, fm_struct)
2090 44 : IF (.NOT. compare_potential_types(mp2_env%ri_metric, mp2_env%potential_parameter)) THEN
2091 4 : CALL cp_fm_create(fm_PQ_2, fm_struct)
2092 : END IF
2093 44 : CALL cp_fm_struct_release(fm_struct)
2094 44 : CALL cp_fm_set_all(fm_PQ, 0.0_dp)
2095 :
2096 : ! We still have to left- and right multiply it with PQhalf
2097 44 : CALL dereplicate_and_sum_fm(rpa_grad%fm_Gamma_PQ, fm_PQ)
2098 :
2099 44 : ngroup = para_env%num_pe/para_env_sub%num_pe
2100 :
2101 : CALL prepare_redistribution(para_env, para_env_sub, ngroup, &
2102 44 : group_grid_2_mepos, mepos_2_grid_group, pos_info=pos_info)
2103 :
2104 : ! Create fm_PQ_half
2105 44 : CALL create_group_dist(gd_P, para_env_sub%num_pe, dimen_RI)
2106 44 : CALL get_group_dist(gd_P, para_env_sub%mepos, my_P_start, my_P_end, my_P_size)
2107 :
2108 44 : CALL get_group_dist(gd_array, color_sub, my_group_L_start, my_group_L_end, my_group_L_size)
2109 :
2110 44 : CALL create_group_dist(gd_P_new, para_env%num_pe)
2111 44 : CALL create_group_dist(gd_array_new, para_env%num_pe)
2112 :
2113 132 : DO proc = 0, para_env%num_pe - 1
2114 : ! calculate position of the group
2115 88 : pos_group = proc/para_env_sub%num_pe
2116 : ! calculate position in the subgroup
2117 88 : pos_sub = pos_info(proc)
2118 : ! 1 -> rows, 2 -> cols
2119 88 : CALL get_group_dist(gd_array, pos_group, gd_array_new, proc)
2120 132 : CALL get_group_dist(gd_P, pos_sub, gd_P_new, proc)
2121 : END DO
2122 :
2123 44 : DEALLOCATE (pos_info)
2124 44 : CALL release_group_dist(gd_P)
2125 :
2126 : CALL array2fm(mp2_env%ri_grad%PQ_half, fm_PQ%matrix_struct, &
2127 : my_P_start, my_P_end, &
2128 : my_group_L_start, my_group_L_end, &
2129 : gd_P_new, gd_array_new, &
2130 : group_grid_2_mepos, para_env_sub%num_pe, ngroup, &
2131 44 : fm_PQ_half)
2132 :
2133 44 : IF (.NOT. compare_potential_types(mp2_env%ri_metric, mp2_env%potential_parameter)) THEN
2134 : CALL array2fm(mp2_env%ri_grad%operator_half, fm_PQ%matrix_struct, my_P_start, my_P_end, &
2135 : my_group_L_start, my_group_L_end, &
2136 : gd_P_new, gd_array_new, &
2137 : group_grid_2_mepos, para_env_sub%num_pe, ngroup, &
2138 4 : operator_half)
2139 : END IF
2140 :
2141 : ! deallocate the info array
2142 44 : CALL release_group_dist(gd_P_new)
2143 44 : CALL release_group_dist(gd_array_new)
2144 :
2145 44 : IF (compare_potential_types(mp2_env%ri_metric, mp2_env%potential_parameter)) THEN
2146 : ! Finish Gamma_PQ
2147 : CALL parallel_gemm(transa="N", transb="T", m=dimen_RI, n=dimen_RI, k=dimen_RI, alpha=1.0_dp, &
2148 : matrix_a=fm_PQ, matrix_b=fm_PQ_half, beta=0.0_dp, &
2149 40 : matrix_c=fm_work_PQ)
2150 :
2151 : CALL parallel_gemm(transa="N", transb="N", m=dimen_RI, n=dimen_RI, k=dimen_RI, alpha=-my_scale, &
2152 : matrix_a=fm_PQ_half, matrix_b=fm_work_PQ, beta=0.0_dp, &
2153 40 : matrix_c=fm_PQ)
2154 :
2155 40 : CALL cp_fm_release(fm_work_PQ)
2156 : ELSE
2157 : CALL parallel_gemm(transa="N", transb="T", m=dimen_RI, n=dimen_RI, k=dimen_RI, alpha=1.0_dp, &
2158 : matrix_a=fm_PQ, matrix_b=operator_half, beta=0.0_dp, &
2159 4 : matrix_c=fm_work_PQ)
2160 :
2161 : CALL parallel_gemm(transa="N", transb="N", m=dimen_RI, n=dimen_RI, k=dimen_RI, alpha=my_scale, &
2162 : matrix_a=operator_half, matrix_b=fm_work_PQ, beta=0.0_dp, &
2163 4 : matrix_c=fm_PQ)
2164 4 : CALL cp_fm_release(operator_half)
2165 :
2166 4 : CALL cp_fm_create(fm_work_PQ_2, fm_PQ%matrix_struct, name="fm_Gamma_PQ_2")
2167 : CALL parallel_gemm(transa="N", transb="N", m=dimen_RI, n=dimen_RI, k=dimen_RI, alpha=-my_scale, &
2168 : matrix_a=fm_PQ_half, matrix_b=fm_work_PQ, beta=0.0_dp, &
2169 4 : matrix_c=fm_work_PQ_2)
2170 4 : CALL cp_fm_to_fm(fm_work_PQ_2, fm_PQ_2)
2171 4 : CALL cp_fm_geadd(1.0_dp, "T", fm_work_PQ_2, 1.0_dp, fm_PQ_2)
2172 4 : CALL cp_fm_release(fm_work_PQ_2)
2173 4 : CALL cp_fm_release(fm_work_PQ)
2174 : END IF
2175 :
2176 176 : ALLOCATE (mp2_env%ri_grad%Gamma_PQ(my_P_size, my_group_L_size))
2177 : CALL fm2array(mp2_env%ri_grad%Gamma_PQ, &
2178 : my_P_start, my_P_end, &
2179 : my_group_L_start, my_group_L_end, &
2180 : group_grid_2_mepos, mepos_2_grid_group, &
2181 : para_env_sub%num_pe, ngroup, &
2182 44 : fm_PQ)
2183 :
2184 44 : IF (.NOT. compare_potential_types(mp2_env%ri_metric, mp2_env%potential_parameter)) THEN
2185 12 : ALLOCATE (mp2_env%ri_grad%Gamma_PQ_2(my_P_size, my_group_L_size))
2186 : CALL fm2array(mp2_env%ri_grad%Gamma_PQ_2, my_P_start, my_P_end, &
2187 : my_group_L_start, my_group_L_end, &
2188 : group_grid_2_mepos, mepos_2_grid_group, &
2189 : para_env_sub%num_pe, ngroup, &
2190 4 : fm_PQ_2)
2191 : END IF
2192 :
2193 : ! Now, Gamma_Pia
2194 2830 : ALLOCATE (mp2_env%ri_grad%G_P_ia(my_group_L_size, nspins))
2195 96 : DO ispin = 1, nspins
2196 2698 : DO iiB = 1, my_group_L_size
2197 2654 : NULLIFY (mp2_env%ri_grad%G_P_ia(iiB, ispin)%matrix)
2198 : END DO
2199 : END DO
2200 :
2201 : ! Redistribute the Y matrix
2202 96 : DO ispin = 1, nspins
2203 : ! Collect all data of columns for the own sub group locally
2204 52 : CALL cp_fm_get_info(rpa_grad%fm_Y(ispin), ncol_global=dimen_ia)
2205 :
2206 52 : CALL get_qs_env(qs_env, blacs_env=blacs_env)
2207 :
2208 52 : NULLIFY (fm_struct)
2209 52 : CALL cp_fm_struct_create(fm_struct, template_fmstruct=fm_PQ_half%matrix_struct, nrow_global=dimen_ia)
2210 52 : CALL cp_fm_create(fm_Y, fm_struct)
2211 52 : CALL cp_fm_struct_release(fm_struct)
2212 52 : CALL cp_fm_set_all(fm_Y, 0.0_dp)
2213 :
2214 52 : CALL dereplicate_and_sum_fm(rpa_grad%fm_Y(ispin), fm_Y)
2215 :
2216 52 : CALL cp_fm_create(fm_G_P_ia, fm_Y%matrix_struct)
2217 52 : CALL cp_fm_set_all(fm_G_P_ia, 0.0_dp)
2218 :
2219 : CALL parallel_gemm(transa="N", transb="T", m=dimen_ia, n=dimen_RI, k=dimen_RI, alpha=my_scale, &
2220 : matrix_a=fm_Y, matrix_b=fm_PQ_half, beta=0.0_dp, &
2221 52 : matrix_c=fm_G_P_ia)
2222 :
2223 52 : CALL cp_fm_release(fm_Y)
2224 :
2225 52 : CALL create_group_dist(gd_ia, para_env_sub%num_pe, dimen_ia)
2226 52 : CALL get_group_dist(gd_ia, para_env_sub%mepos, my_ia_start, my_ia_end, my_ia_size)
2227 :
2228 : CALL fm2array(Gamma_2D, my_ia_start, my_ia_end, &
2229 : my_group_L_start, my_group_L_end, &
2230 : group_grid_2_mepos, mepos_2_grid_group, &
2231 : para_env_sub%num_pe, ngroup, &
2232 52 : fm_G_P_ia)
2233 :
2234 : ! create the Gamma_ia_P in DBCSR style
2235 : CALL create_dbcsr_gamma(Gamma_2D, homo(ispin), virtual(ispin), dimen_ia, para_env_sub, &
2236 : my_ia_start, my_ia_end, my_group_L_size, gd_ia, &
2237 52 : mp2_env%ri_grad%G_P_ia(:, ispin), mp2_env%ri_grad%mo_coeff_o(ispin)%matrix)
2238 :
2239 304 : CALL release_group_dist(gd_ia)
2240 :
2241 : END DO
2242 44 : DEALLOCATE (rpa_grad%fm_Y)
2243 44 : CALL cp_fm_release(fm_PQ_half)
2244 :
2245 44 : CALL timestop(handle)
2246 :
2247 264 : END SUBROUTINE rpa_grad_finalize
2248 :
2249 : ! **************************************************************************************************
2250 : !> \brief ...
2251 : !> \param sos_mp2_work_occ ...
2252 : !> \param sos_mp2_work_virt ...
2253 : !> \param para_env ...
2254 : !> \param para_env_sub ...
2255 : !> \param homo ...
2256 : !> \param virtual ...
2257 : !> \param mp2_env ...
2258 : ! **************************************************************************************************
2259 20 : SUBROUTINE sos_mp2_grad_finalize(sos_mp2_work_occ, sos_mp2_work_virt, para_env, para_env_sub, homo, virtual, mp2_env)
2260 : TYPE(sos_mp2_grad_work_type), ALLOCATABLE, &
2261 : DIMENSION(:), INTENT(INOUT) :: sos_mp2_work_occ, sos_mp2_work_virt
2262 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env, para_env_sub
2263 : INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
2264 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
2265 :
2266 : CHARACTER(LEN=*), PARAMETER :: routineN = 'sos_mp2_grad_finalize'
2267 :
2268 : INTEGER :: ab_counter, handle, ij_counter, ispin, &
2269 : itmp(2), my_a, my_b, my_B_end, &
2270 : my_B_size, my_B_start, my_i, my_j, &
2271 : nspins, pcol
2272 : REAL(KIND=dp) :: my_scale
2273 :
2274 20 : CALL timeset(routineN, handle)
2275 :
2276 20 : nspins = SIZE(sos_mp2_work_occ)
2277 20 : my_scale = mp2_env%scale_s
2278 :
2279 44 : DO ispin = 1, nspins
2280 48 : DO pcol = 0, SIZE(sos_mp2_work_occ(ispin)%index2send, 1) - 1
2281 24 : IF (ALLOCATED(sos_mp2_work_occ(ispin)%index2send(pcol)%array)) THEN
2282 4 : DEALLOCATE (sos_mp2_work_occ(ispin)%index2send(pcol)%array)
2283 : END IF
2284 24 : IF (ALLOCATED(sos_mp2_work_occ(ispin)%index2send(pcol)%array)) THEN
2285 0 : DEALLOCATE (sos_mp2_work_occ(ispin)%index2send(pcol)%array)
2286 : END IF
2287 24 : IF (ALLOCATED(sos_mp2_work_virt(ispin)%index2recv(pcol)%array)) THEN
2288 4 : DEALLOCATE (sos_mp2_work_virt(ispin)%index2recv(pcol)%array)
2289 : END IF
2290 48 : IF (ALLOCATED(sos_mp2_work_virt(ispin)%index2recv(pcol)%array)) THEN
2291 0 : DEALLOCATE (sos_mp2_work_virt(ispin)%index2recv(pcol)%array)
2292 : END IF
2293 : END DO
2294 0 : DEALLOCATE (sos_mp2_work_occ(ispin)%index2send, &
2295 0 : sos_mp2_work_occ(ispin)%index2recv, &
2296 0 : sos_mp2_work_virt(ispin)%index2send, &
2297 140 : sos_mp2_work_virt(ispin)%index2recv)
2298 : END DO
2299 :
2300 : ! Sum P_ij and P_ab and redistribute them
2301 44 : DO ispin = 1, nspins
2302 24 : CALL para_env%sum(sos_mp2_work_occ(ispin)%P)
2303 :
2304 96 : ALLOCATE (mp2_env%ri_grad%P_ij(ispin)%array(homo(ispin), homo(ispin)))
2305 472 : mp2_env%ri_grad%P_ij(ispin)%array = 0.0_dp
2306 116 : DO my_i = 1, homo(ispin)
2307 116 : mp2_env%ri_grad%P_ij(ispin)%array(my_i, my_i) = my_scale*sos_mp2_work_occ(ispin)%P(my_i)
2308 : END DO
2309 72 : DO ij_counter = 1, SIZE(sos_mp2_work_occ(ispin)%pair_list, 2)
2310 48 : my_i = sos_mp2_work_occ(ispin)%pair_list(1, ij_counter)
2311 48 : my_j = sos_mp2_work_occ(ispin)%pair_list(2, ij_counter)
2312 :
2313 72 : mp2_env%ri_grad%P_ij(ispin)%array(my_i, my_j) = my_scale*sos_mp2_work_occ(ispin)%P(homo(ispin) + ij_counter)
2314 : END DO
2315 24 : DEALLOCATE (sos_mp2_work_occ(ispin)%P, sos_mp2_work_occ(ispin)%pair_list)
2316 :
2317 : ! Symmetrize P_ij
2318 : mp2_env%ri_grad%P_ij(ispin)%array(:, :) = 0.5_dp*(mp2_env%ri_grad%P_ij(ispin)%array + &
2319 920 : TRANSPOSE(mp2_env%ri_grad%P_ij(ispin)%array))
2320 :
2321 : ! The first index of P_ab has to be distributed within the subgroups,
2322 : ! so sum it up first and add the required elements later
2323 24 : CALL para_env%sum(sos_mp2_work_virt(ispin)%P)
2324 :
2325 24 : itmp = get_limit(virtual(ispin), para_env_sub%num_pe, para_env_sub%mepos)
2326 24 : my_B_size = itmp(2) - itmp(1) + 1
2327 24 : my_B_start = itmp(1)
2328 24 : my_B_end = itmp(2)
2329 :
2330 96 : ALLOCATE (mp2_env%ri_grad%P_ab(ispin)%array(my_B_size, virtual(ispin)))
2331 10382 : mp2_env%ri_grad%P_ab(ispin)%array = 0.0_dp
2332 486 : DO my_a = itmp(1), itmp(2)
2333 486 : mp2_env%ri_grad%P_ab(ispin)%array(my_a - itmp(1) + 1, my_a) = my_scale*sos_mp2_work_virt(ispin)%P(my_a)
2334 : END DO
2335 636 : DO ab_counter = 1, SIZE(sos_mp2_work_virt(ispin)%pair_list, 2)
2336 612 : my_a = sos_mp2_work_virt(ispin)%pair_list(1, ab_counter)
2337 612 : my_b = sos_mp2_work_virt(ispin)%pair_list(2, ab_counter)
2338 :
2339 612 : IF (my_a >= itmp(1) .AND. my_a <= itmp(2)) mp2_env%ri_grad%P_ab(ispin)%array(my_a - itmp(1) + 1, my_b) = &
2340 636 : my_scale*sos_mp2_work_virt(ispin)%P(virtual(ispin) + ab_counter)
2341 : END DO
2342 :
2343 24 : DEALLOCATE (sos_mp2_work_virt(ispin)%P, sos_mp2_work_virt(ispin)%pair_list)
2344 :
2345 : ! Symmetrize P_ab
2346 44 : IF (para_env_sub%num_pe > 1) THEN
2347 12 : BLOCK
2348 : INTEGER :: send_a_start, send_a_end, send_a_size, &
2349 : recv_a_start, recv_a_end, recv_a_size, proc_shift, proc_send, proc_recv
2350 4 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE, TARGET :: buffer_send_1D
2351 4 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: buffer_send
2352 4 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: buffer_recv
2353 4 : TYPE(group_dist_d1_type) :: gd_virtual_sub
2354 :
2355 4 : CALL create_group_dist(gd_virtual_sub, para_env_sub%num_pe, virtual(ispin))
2356 :
2357 : mp2_env%ri_grad%P_ab(ispin)%array(:, my_B_start:my_B_end) = &
2358 : 0.5_dp*(mp2_env%ri_grad%P_ab(ispin)%array(:, my_B_start:my_B_end) &
2359 804 : + TRANSPOSE(mp2_env%ri_grad%P_ab(ispin)%array(:, my_B_start:my_B_end)))
2360 :
2361 12 : ALLOCATE (buffer_send_1D(my_B_size*maxsize(gd_virtual_sub)))
2362 16 : ALLOCATE (buffer_recv(my_B_size, maxsize(gd_virtual_sub)))
2363 :
2364 8 : DO proc_shift = 1, para_env_sub%num_pe - 1
2365 :
2366 4 : proc_send = MODULO(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
2367 4 : proc_recv = MODULO(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
2368 :
2369 4 : CALL get_group_dist(gd_virtual_sub, proc_send, send_a_start, send_a_end, send_a_size)
2370 4 : CALL get_group_dist(gd_virtual_sub, proc_recv, recv_a_start, recv_a_end, recv_a_size)
2371 :
2372 4 : buffer_send(1:send_a_size, 1:my_B_size) => buffer_send_1D(1:my_B_size*send_a_size)
2373 :
2374 402 : buffer_send(:send_a_size, :) = TRANSPOSE(mp2_env%ri_grad%P_ab(ispin)%array(:, send_a_start:send_a_end))
2375 : CALL para_env_sub%sendrecv(buffer_send(:send_a_size, :), proc_send, &
2376 402 : buffer_recv(:, :recv_a_size), proc_recv)
2377 :
2378 : mp2_env%ri_grad%P_ab(ispin)%array(:, recv_a_start:recv_a_end) = &
2379 410 : 0.5_dp*(mp2_env%ri_grad%P_ab(ispin)%array(:, recv_a_start:recv_a_end) + buffer_recv(:, 1:recv_a_size))
2380 :
2381 : END DO
2382 :
2383 4 : DEALLOCATE (buffer_send_1D, buffer_recv)
2384 :
2385 16 : CALL release_group_dist(gd_virtual_sub)
2386 : END BLOCK
2387 : ELSE
2388 : mp2_env%ri_grad%P_ab(ispin)%array(:, :) = 0.5_dp*(mp2_env%ri_grad%P_ab(ispin)%array + &
2389 19140 : TRANSPOSE(mp2_env%ri_grad%P_ab(ispin)%array))
2390 : END IF
2391 :
2392 : END DO
2393 68 : DEALLOCATE (sos_mp2_work_occ, sos_mp2_work_virt)
2394 20 : IF (nspins == 1) THEN
2395 336 : mp2_env%ri_grad%P_ij(1)%array(:, :) = 2.0_dp*mp2_env%ri_grad%P_ij(1)%array
2396 5374 : mp2_env%ri_grad%P_ab(1)%array(:, :) = 2.0_dp*mp2_env%ri_grad%P_ab(1)%array
2397 : END IF
2398 :
2399 20 : CALL timestop(handle)
2400 :
2401 20 : END SUBROUTINE sos_mp2_grad_finalize
2402 :
2403 : ! **************************************************************************************************
2404 : !> \brief ...
2405 : !> \param rpa_work ...
2406 : !> \param mp2_env ...
2407 : !> \param homo ...
2408 : !> \param virtual ...
2409 : !> \param para_env ...
2410 : !> \param para_env_sub ...
2411 : ! **************************************************************************************************
2412 24 : SUBROUTINE rpa_grad_work_finalize(rpa_work, mp2_env, homo, virtual, para_env, para_env_sub)
2413 : TYPE(rpa_grad_work_type), INTENT(INOUT) :: rpa_work
2414 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
2415 : INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
2416 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env, para_env_sub
2417 :
2418 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rpa_grad_work_finalize'
2419 :
2420 : INTEGER :: handle, ispin, itmp(2), my_a_end, my_a_size, my_a_start, my_B_end, my_B_size, &
2421 : my_B_start, my_i_end, my_i_size, my_i_start, nspins, proc, proc_recv, proc_send, &
2422 : proc_shift, recv_a_end, recv_a_size, recv_a_start, recv_end, recv_start, send_a_end, &
2423 : send_a_size, send_a_start, send_end, send_start, size_recv_buffer, size_send_buffer
2424 : REAL(KIND=dp) :: my_scale
2425 24 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: buffer_recv, buffer_send
2426 24 : TYPE(group_dist_d1_type) :: gd_a_sub, gd_virtual_sub
2427 :
2428 24 : CALL timeset(routineN, handle)
2429 :
2430 24 : nspins = SIZE(homo)
2431 24 : my_scale = mp2_env%ri_rpa%scale_rpa/(2.0_dp*pi)
2432 24 : IF (mp2_env%ri_rpa%minimax_quad) my_scale = my_scale/2.0_dp
2433 :
2434 24 : CALL cp_fm_release(rpa_work%fm_mat_Q_copy)
2435 :
2436 52 : DO ispin = 1, nspins
2437 80 : DO proc = 0, SIZE(rpa_work%index2send, 1) - 1
2438 28 : IF (ALLOCATED(rpa_work%index2send(proc, ispin)%array)) DEALLOCATE (rpa_work%index2send(proc, ispin)%array)
2439 56 : IF (ALLOCATED(rpa_work%index2recv(proc, ispin)%array)) DEALLOCATE (rpa_work%index2recv(proc, ispin)%array)
2440 : END DO
2441 : END DO
2442 80 : DEALLOCATE (rpa_work%index2send, rpa_work%index2recv)
2443 :
2444 52 : DO ispin = 1, nspins
2445 28 : CALL get_group_dist(rpa_work%gd_homo(ispin), rpa_work%mepos(2), my_i_start, my_i_end, my_i_size)
2446 28 : CALL release_group_dist(rpa_work%gd_homo(ispin))
2447 :
2448 112 : ALLOCATE (mp2_env%ri_grad%P_ij(ispin)%array(homo(ispin), homo(ispin)))
2449 556 : mp2_env%ri_grad%P_ij(ispin)%array = 0.0_dp
2450 556 : mp2_env%ri_grad%P_ij(ispin)%array(my_i_start:my_i_end, :) = my_scale*rpa_work%P_ij(ispin)%array
2451 28 : DEALLOCATE (rpa_work%P_ij(ispin)%array)
2452 28 : CALL para_env%sum(mp2_env%ri_grad%P_ij(ispin)%array)
2453 :
2454 : ! Symmetrize P_ij
2455 : mp2_env%ri_grad%P_ij(ispin)%array(:, :) = 0.5_dp*(mp2_env%ri_grad%P_ij(ispin)%array + &
2456 1084 : TRANSPOSE(mp2_env%ri_grad%P_ij(ispin)%array))
2457 :
2458 28 : itmp = get_limit(virtual(ispin), para_env_sub%num_pe, para_env_sub%mepos)
2459 28 : my_B_start = itmp(1)
2460 28 : my_B_end = itmp(2)
2461 28 : my_B_size = my_B_end - my_B_start + 1
2462 :
2463 112 : ALLOCATE (mp2_env%ri_grad%P_ab(ispin)%array(my_B_size, virtual(ispin)))
2464 11610 : mp2_env%ri_grad%P_ab(ispin)%array = 0.0_dp
2465 :
2466 28 : CALL get_group_dist(rpa_work%gd_virtual(ispin), rpa_work%mepos(1), my_a_start, my_a_end, my_a_size)
2467 28 : CALL release_group_dist(rpa_work%gd_virtual(ispin))
2468 : ! This group dist contains the info which parts of Pab a process currently owns
2469 28 : CALL create_group_dist(gd_a_sub, my_a_start, my_a_end, my_a_size, para_env_sub)
2470 : ! This group dist contains the info which parts of Pab a process is supposed to own later
2471 28 : CALL create_group_dist(gd_virtual_sub, para_env_sub%num_pe, virtual(ispin))
2472 :
2473 : ! Calculate local indices of the common range of own matrix and send process
2474 28 : send_start = MAX(1, my_B_start - my_a_start + 1)
2475 28 : send_end = MIN(my_a_size, my_B_end - my_a_start + 1)
2476 :
2477 : ! Same for recv process but with reverse positions
2478 28 : recv_start = MAX(1, my_a_start - my_B_start + 1)
2479 28 : recv_end = MIN(my_B_size, my_a_end - my_B_start + 1)
2480 :
2481 : mp2_env%ri_grad%P_ab(ispin)%array(recv_start:recv_end, :) = &
2482 11610 : my_scale*rpa_work%P_ab(ispin)%array(send_start:send_end, :)
2483 :
2484 28 : IF (para_env_sub%num_pe > 1) THEN
2485 : size_send_buffer = 0
2486 : size_recv_buffer = 0
2487 8 : DO proc_shift = 1, para_env_sub%num_pe - 1
2488 4 : proc_send = MODULO(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
2489 4 : proc_recv = MODULO(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
2490 :
2491 4 : CALL get_group_dist(gd_virtual_sub, proc_send, send_a_start, send_a_end)
2492 4 : CALL get_group_dist(gd_a_sub, proc_recv, recv_a_start, recv_a_end)
2493 :
2494 : ! Calculate local indices of the common range of own matrix and send process
2495 4 : send_start = MAX(1, send_a_start - my_a_start + 1)
2496 4 : send_end = MIN(my_a_size, send_a_end - my_a_start + 1)
2497 :
2498 4 : size_send_buffer = MAX(size_send_buffer, MAX(send_end - send_start + 1, 0))
2499 :
2500 : ! Same for recv process but with reverse positions
2501 4 : recv_start = MAX(1, recv_a_start - my_B_start + 1)
2502 4 : recv_end = MIN(my_B_size, recv_a_end - my_B_start + 1)
2503 :
2504 8 : size_recv_buffer = MAX(size_recv_buffer, MAX(recv_end - recv_start + 1, 0))
2505 : END DO
2506 28 : ALLOCATE (buffer_send(size_send_buffer, virtual(ispin)), buffer_recv(size_recv_buffer, virtual(ispin)))
2507 :
2508 8 : DO proc_shift = 1, para_env_sub%num_pe - 1
2509 4 : proc_send = MODULO(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
2510 4 : proc_recv = MODULO(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
2511 :
2512 4 : CALL get_group_dist(gd_virtual_sub, proc_send, send_a_start, send_a_end)
2513 4 : CALL get_group_dist(gd_a_sub, proc_recv, recv_a_start, recv_a_end)
2514 :
2515 : ! Calculate local indices of the common range of own matrix and send process
2516 4 : send_start = MAX(1, send_a_start - my_a_start + 1)
2517 4 : send_end = MIN(my_a_size, send_a_end - my_a_start + 1)
2518 802 : buffer_send(1:MAX(send_end - send_start + 1, 0), :) = rpa_work%P_ab(ispin)%array(send_start:send_end, :)
2519 :
2520 : ! Same for recv process but with reverse positions
2521 4 : recv_start = MAX(1, recv_a_start - my_B_start + 1)
2522 4 : recv_end = MIN(my_B_size, recv_a_end - my_B_start + 1)
2523 :
2524 : CALL para_env_sub%sendrecv(buffer_send(1:MAX(send_end - send_start + 1, 0), :), proc_send, &
2525 1600 : buffer_recv(1:MAX(recv_end - recv_start + 1, 0), :), proc_recv)
2526 :
2527 : mp2_env%ri_grad%P_ab(ispin)%array(recv_start:recv_end, :) = &
2528 : mp2_env%ri_grad%P_ab(ispin)%array(recv_start:recv_end, :) + &
2529 810 : my_scale*buffer_recv(1:MAX(recv_end - recv_start + 1, 0), :)
2530 :
2531 : END DO
2532 :
2533 4 : IF (ALLOCATED(buffer_send)) DEALLOCATE (buffer_send)
2534 4 : IF (ALLOCATED(buffer_recv)) DEALLOCATE (buffer_recv)
2535 : END IF
2536 28 : DEALLOCATE (rpa_work%P_ab(ispin)%array)
2537 :
2538 28 : CALL release_group_dist(gd_a_sub)
2539 :
2540 : BLOCK
2541 : TYPE(mp_comm_type) :: comm_exchange
2542 28 : CALL comm_exchange%from_split(para_env, para_env_sub%mepos)
2543 28 : CALL comm_exchange%sum(mp2_env%ri_grad%P_ab(ispin)%array)
2544 56 : CALL comm_exchange%free()
2545 : END BLOCK
2546 :
2547 : ! Symmetrize P_ab
2548 28 : IF (para_env_sub%num_pe > 1) THEN
2549 : BLOCK
2550 4 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE, TARGET :: buffer_send_1D
2551 4 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: buffer_send
2552 4 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: buffer_recv
2553 :
2554 : mp2_env%ri_grad%P_ab(ispin)%array(:, my_B_start:my_B_end) = &
2555 : 0.5_dp*(mp2_env%ri_grad%P_ab(ispin)%array(:, my_B_start:my_B_end) &
2556 804 : + TRANSPOSE(mp2_env%ri_grad%P_ab(ispin)%array(:, my_B_start:my_B_end)))
2557 :
2558 12 : ALLOCATE (buffer_send_1D(my_B_size*maxsize(gd_virtual_sub)))
2559 16 : ALLOCATE (buffer_recv(my_B_size, maxsize(gd_virtual_sub)))
2560 :
2561 8 : DO proc_shift = 1, para_env_sub%num_pe - 1
2562 :
2563 4 : proc_send = MODULO(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
2564 4 : proc_recv = MODULO(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
2565 :
2566 4 : CALL get_group_dist(gd_virtual_sub, proc_send, send_a_start, send_a_end, send_a_size)
2567 4 : CALL get_group_dist(gd_virtual_sub, proc_recv, recv_a_start, recv_a_end, recv_a_size)
2568 :
2569 4 : buffer_send(1:send_a_size, 1:my_B_size) => buffer_send_1D(1:my_B_size*send_a_size)
2570 :
2571 402 : buffer_send(:send_a_size, :) = TRANSPOSE(mp2_env%ri_grad%P_ab(ispin)%array(:, send_a_start:send_a_end))
2572 : CALL para_env_sub%sendrecv(buffer_send(:send_a_size, :), proc_send, &
2573 402 : buffer_recv(:, :recv_a_size), proc_recv)
2574 :
2575 : mp2_env%ri_grad%P_ab(ispin)%array(:, recv_a_start:recv_a_end) = &
2576 410 : 0.5_dp*(mp2_env%ri_grad%P_ab(ispin)%array(:, recv_a_start:recv_a_end) + buffer_recv(:, 1:recv_a_size))
2577 :
2578 : END DO
2579 :
2580 8 : DEALLOCATE (buffer_send_1D, buffer_recv)
2581 : END BLOCK
2582 : ELSE
2583 : mp2_env%ri_grad%P_ab(ispin)%array(:, :) = 0.5_dp*(mp2_env%ri_grad%P_ab(ispin)%array + &
2584 21592 : TRANSPOSE(mp2_env%ri_grad%P_ab(ispin)%array))
2585 : END IF
2586 :
2587 108 : CALL release_group_dist(gd_virtual_sub)
2588 :
2589 : END DO
2590 136 : DEALLOCATE (rpa_work%gd_homo, rpa_work%gd_virtual, rpa_work%P_ij, rpa_work%P_ab)
2591 24 : IF (nspins == 1) THEN
2592 420 : mp2_env%ri_grad%P_ij(1)%array(:, :) = 2.0_dp*mp2_env%ri_grad%P_ij(1)%array
2593 6602 : mp2_env%ri_grad%P_ab(1)%array(:, :) = 2.0_dp*mp2_env%ri_grad%P_ab(1)%array
2594 : END IF
2595 :
2596 24 : CALL timestop(handle)
2597 24 : END SUBROUTINE rpa_grad_work_finalize
2598 :
2599 : ! **************************************************************************************************
2600 : !> \brief Dereplicate data from fm_sub and collect in fm_global, overlapping data will be added
2601 : !> \param fm_sub replicated matrix, all subgroups have the same size, will be release on output
2602 : !> \param fm_global global matrix, on output it will contain the sum of the replicated matrices redistributed
2603 : ! **************************************************************************************************
2604 96 : SUBROUTINE dereplicate_and_sum_fm(fm_sub, fm_global)
2605 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_sub, fm_global
2606 :
2607 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dereplicate_and_sum_fm'
2608 :
2609 : INTEGER :: col_local, elements2recv_col, elements2recv_row, elements2send_col, &
2610 : elements2send_row, handle, handle2, mypcol_global, myprow_global, ncol_local_global, &
2611 : ncol_local_sub, npcol_global, npcol_sub, nprow_global, nprow_sub, nrow_local_global, &
2612 : nrow_local_sub, pcol_recv, pcol_send, proc_recv, proc_send, proc_send_global, proc_shift, &
2613 : prow_recv, prow_send, row_local, tag
2614 : INTEGER(int_8) :: size_recv_buffer, size_send_buffer
2615 96 : INTEGER, ALLOCATABLE, DIMENSION(:) :: data2recv_col, data2recv_row, &
2616 96 : data2send_col, data2send_row, &
2617 96 : subgroup2mepos
2618 96 : INTEGER, DIMENSION(:), POINTER :: col_indices_global, col_indices_sub, &
2619 96 : row_indices_global, row_indices_sub
2620 96 : INTEGER, DIMENSION(:, :), POINTER :: blacs2mpi_global, blacs2mpi_sub, &
2621 96 : mpi2blacs_global, mpi2blacs_sub
2622 96 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), TARGET :: recv_buffer_1D, send_buffer_1D
2623 96 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: recv_buffer, send_buffer
2624 : TYPE(mp_para_env_type), POINTER :: para_env, para_env_sub
2625 96 : TYPE(one_dim_int_array), ALLOCATABLE, DIMENSION(:) :: index2recv_col, index2recv_row, &
2626 96 : index2send_col, index2send_row
2627 :
2628 96 : CALL timeset(routineN, handle)
2629 :
2630 96 : tag = 1
2631 :
2632 96 : nprow_sub = fm_sub%matrix_struct%context%num_pe(1)
2633 96 : npcol_sub = fm_sub%matrix_struct%context%num_pe(2)
2634 :
2635 96 : myprow_global = fm_global%matrix_struct%context%mepos(1)
2636 96 : mypcol_global = fm_global%matrix_struct%context%mepos(2)
2637 96 : nprow_global = fm_global%matrix_struct%context%num_pe(1)
2638 96 : npcol_global = fm_global%matrix_struct%context%num_pe(2)
2639 :
2640 : CALL cp_fm_get_info(fm_sub, col_indices=col_indices_sub, row_indices=row_indices_sub, &
2641 96 : nrow_local=nrow_local_sub, ncol_local=ncol_local_sub)
2642 96 : CALL cp_fm_struct_get(fm_sub%matrix_struct, para_env=para_env_sub)
2643 : CALL cp_fm_struct_get(fm_global%matrix_struct, para_env=para_env, &
2644 : col_indices=col_indices_global, row_indices=row_indices_global, &
2645 96 : nrow_local=nrow_local_global, ncol_local=ncol_local_global)
2646 96 : CALL fm_sub%matrix_struct%context%get(blacs2mpi=blacs2mpi_sub, mpi2blacs=mpi2blacs_sub)
2647 96 : CALL fm_global%matrix_struct%context%get(blacs2mpi=blacs2mpi_global, mpi2blacs=mpi2blacs_global)
2648 :
2649 96 : IF (para_env%num_pe /= para_env_sub%num_pe) THEN
2650 : BLOCK
2651 : TYPE(mp_comm_type) :: comm_exchange
2652 80 : comm_exchange = fm_sub%matrix_struct%context%interconnect(para_env)
2653 80 : CALL comm_exchange%sum(fm_sub%local_data)
2654 160 : CALL comm_exchange%free()
2655 : END BLOCK
2656 : END IF
2657 :
2658 288 : ALLOCATE (subgroup2mepos(0:para_env_sub%num_pe - 1))
2659 96 : CALL para_env_sub%allgather(para_env%mepos, subgroup2mepos)
2660 :
2661 96 : CALL timeset(routineN//"_data2", handle2)
2662 : ! Create a map how much data has to be sent to what process coordinate, interchange rows and columns to transpose the matrices
2663 96 : CALL get_elements2send_col(data2send_col, fm_global%matrix_struct, row_indices_sub, index2send_col)
2664 96 : CALL get_elements2send_row(data2send_row, fm_global%matrix_struct, col_indices_sub, index2send_row)
2665 :
2666 : ! Create a map how much data has to be sent to what process coordinate, interchange rows and columns to transpose the matrices
2667 : ! Do the reverse for the recieve processes
2668 96 : CALL get_elements2send_col(data2recv_col, fm_sub%matrix_struct, row_indices_global, index2recv_col)
2669 96 : CALL get_elements2send_row(data2recv_row, fm_sub%matrix_struct, col_indices_global, index2recv_row)
2670 96 : CALL timestop(handle2)
2671 :
2672 96 : CALL timeset(routineN//"_local", handle2)
2673 : ! Loop over local data and transpose
2674 96 : prow_send = mpi2blacs_global(1, para_env%mepos)
2675 96 : pcol_send = mpi2blacs_global(2, para_env%mepos)
2676 96 : prow_recv = mpi2blacs_sub(1, para_env_sub%mepos)
2677 96 : pcol_recv = mpi2blacs_sub(2, para_env_sub%mepos)
2678 96 : elements2recv_col = data2recv_col(pcol_recv)
2679 96 : elements2recv_row = data2recv_row(prow_recv)
2680 :
2681 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(row_local,col_local) &
2682 : !$OMP SHARED(elements2recv_col,elements2recv_row,recv_buffer,fm_global,&
2683 : !$OMP index2recv_col,index2recv_row,pcol_recv,prow_recv, &
2684 96 : !$OMP fm_sub,index2send_col,index2send_row,pcol_send,prow_send)
2685 : DO col_local = 1, elements2recv_col
2686 : DO row_local = 1, elements2recv_row
2687 : fm_global%local_data(index2recv_col(pcol_recv)%array(col_local), &
2688 : index2recv_row(prow_recv)%array(row_local)) &
2689 : = fm_sub%local_data(index2send_col(pcol_send)%array(row_local), &
2690 : index2send_row(prow_send)%array(col_local))
2691 : END DO
2692 : END DO
2693 : !$OMP END PARALLEL DO
2694 96 : CALL timestop(handle2)
2695 :
2696 96 : IF (para_env_sub%num_pe > 1) THEN
2697 : size_send_buffer = 0_int_8
2698 : size_recv_buffer = 0_int_8
2699 : ! Loop over all processes in para_env_sub
2700 32 : DO proc_shift = 1, para_env_sub%num_pe - 1
2701 16 : proc_send = MODULO(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
2702 16 : proc_recv = MODULO(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
2703 :
2704 16 : proc_send_global = subgroup2mepos(proc_send)
2705 16 : prow_send = mpi2blacs_global(1, proc_send_global)
2706 16 : pcol_send = mpi2blacs_global(2, proc_send_global)
2707 16 : elements2send_col = data2send_col(pcol_send)
2708 16 : elements2send_row = data2send_row(prow_send)
2709 :
2710 16 : size_send_buffer = MAX(size_send_buffer, INT(elements2send_col, int_8)*elements2send_row)
2711 :
2712 16 : prow_recv = mpi2blacs_sub(1, proc_recv)
2713 16 : pcol_recv = mpi2blacs_sub(2, proc_recv)
2714 16 : elements2recv_col = data2recv_col(pcol_recv)
2715 16 : elements2recv_row = data2recv_row(prow_recv)
2716 :
2717 32 : size_recv_buffer = MAX(size_recv_buffer, INT(elements2recv_col, int_8)*elements2recv_row)
2718 : END DO
2719 80 : ALLOCATE (send_buffer_1D(size_send_buffer), recv_buffer_1D(size_recv_buffer))
2720 :
2721 : ! Loop over all processes in para_env_sub
2722 32 : DO proc_shift = 1, para_env_sub%num_pe - 1
2723 16 : proc_send = MODULO(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
2724 16 : proc_recv = MODULO(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
2725 :
2726 16 : proc_send_global = subgroup2mepos(proc_send)
2727 16 : prow_send = mpi2blacs_global(1, proc_send_global)
2728 16 : pcol_send = mpi2blacs_global(2, proc_send_global)
2729 16 : elements2send_col = data2send_col(pcol_send)
2730 16 : elements2send_row = data2send_row(prow_send)
2731 :
2732 16 : CALL timeset(routineN//"_pack", handle2)
2733 : ! Loop over local data and pack the buffer
2734 : ! Transpose the matrix already
2735 16 : send_buffer(1:elements2send_row, 1:elements2send_col) => send_buffer_1D(1:INT(elements2send_row, int_8)*elements2send_col)
2736 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(row_local,col_local) &
2737 : !$OMP SHARED(elements2send_col,elements2send_row,send_buffer,fm_sub,&
2738 16 : !$OMP index2send_col,index2send_row,pcol_send,prow_send)
2739 : DO row_local = 1, elements2send_col
2740 : DO col_local = 1, elements2send_row
2741 : send_buffer(col_local, row_local) = &
2742 : fm_sub%local_data(index2send_col(pcol_send)%array(row_local), &
2743 : index2send_row(prow_send)%array(col_local))
2744 : END DO
2745 : END DO
2746 : !$OMP END PARALLEL DO
2747 16 : CALL timestop(handle2)
2748 :
2749 16 : prow_recv = mpi2blacs_sub(1, proc_recv)
2750 16 : pcol_recv = mpi2blacs_sub(2, proc_recv)
2751 16 : elements2recv_col = data2recv_col(pcol_recv)
2752 16 : elements2recv_row = data2recv_row(prow_recv)
2753 :
2754 : ! Send data
2755 16 : recv_buffer(1:elements2recv_col, 1:elements2recv_row) => recv_buffer_1D(1:INT(elements2recv_row, int_8)*elements2recv_col)
2756 80 : IF (SIZE(recv_buffer) > 0_int_8) THEN
2757 48 : IF (SIZE(send_buffer) > 0_int_8) THEN
2758 54048 : CALL para_env_sub%sendrecv(send_buffer, proc_send, recv_buffer, proc_recv, tag)
2759 : ELSE
2760 0 : CALL para_env_sub%recv(recv_buffer, proc_recv, tag)
2761 : END IF
2762 :
2763 16 : CALL timeset(routineN//"_unpack", handle2)
2764 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(row_local,col_local) &
2765 : !$OMP SHARED(elements2recv_col,elements2recv_row,recv_buffer,fm_global,&
2766 16 : !$OMP index2recv_col,index2recv_row,pcol_recv,prow_recv)
2767 : DO col_local = 1, elements2recv_col
2768 : DO row_local = 1, elements2recv_row
2769 : fm_global%local_data(index2recv_col(pcol_recv)%array(col_local), &
2770 : index2recv_row(prow_recv)%array(row_local)) &
2771 : = recv_buffer(col_local, row_local)
2772 : END DO
2773 : END DO
2774 : !$OMP END PARALLEL DO
2775 16 : CALL timestop(handle2)
2776 0 : ELSE IF (SIZE(send_buffer) > 0_int_8) THEN
2777 0 : CALL para_env_sub%send(send_buffer, proc_send, tag)
2778 : END IF
2779 : END DO
2780 : END IF
2781 :
2782 96 : DEALLOCATE (data2send_col, data2send_row, data2recv_col, data2recv_row)
2783 192 : DO proc_shift = 0, npcol_global - 1
2784 192 : DEALLOCATE (index2send_col(proc_shift)%array)
2785 : END DO
2786 192 : DO proc_shift = 0, npcol_sub - 1
2787 192 : DEALLOCATE (index2recv_col(proc_shift)%array)
2788 : END DO
2789 288 : DO proc_shift = 0, nprow_global - 1
2790 288 : DEALLOCATE (index2send_row(proc_shift)%array)
2791 : END DO
2792 208 : DO proc_shift = 0, nprow_sub - 1
2793 208 : DEALLOCATE (index2recv_row(proc_shift)%array)
2794 : END DO
2795 592 : DEALLOCATE (index2send_col, index2recv_col, index2send_row, index2recv_row)
2796 :
2797 96 : CALL cp_fm_release(fm_sub)
2798 :
2799 96 : CALL timestop(handle)
2800 :
2801 480 : END SUBROUTINE dereplicate_and_sum_fm
2802 :
2803 : ! **************************************************************************************************
2804 : !> \brief ...
2805 : !> \param data2send ...
2806 : !> \param struct_global ...
2807 : !> \param indices_sub ...
2808 : !> \param index2send ...
2809 : ! **************************************************************************************************
2810 192 : SUBROUTINE get_elements2send_col(data2send, struct_global, indices_sub, index2send)
2811 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: data2send
2812 : TYPE(cp_fm_struct_type), INTENT(INOUT) :: struct_global
2813 : INTEGER, DIMENSION(:), INTENT(IN) :: indices_sub
2814 : TYPE(one_dim_int_array), ALLOCATABLE, &
2815 : DIMENSION(:), INTENT(OUT) :: index2send
2816 :
2817 : INTEGER :: i_global, i_local, np_global, proc
2818 :
2819 192 : CALL struct_global%context%get(number_of_process_columns=np_global)
2820 :
2821 576 : ALLOCATE (data2send(0:np_global - 1))
2822 192 : data2send = 0
2823 11750 : DO i_local = 1, SIZE(indices_sub)
2824 11558 : i_global = indices_sub(i_local)
2825 11558 : proc = struct_global%g2p_col(i_global)
2826 11750 : data2send(proc) = data2send(proc) + 1
2827 : END DO
2828 :
2829 768 : ALLOCATE (index2send(0:np_global - 1))
2830 384 : DO proc = 0, np_global - 1
2831 576 : ALLOCATE (index2send(proc)%array(data2send(proc)))
2832 : ! We want to crash if there is an error
2833 11942 : index2send(proc)%array = -1
2834 : END DO
2835 :
2836 192 : data2send = 0
2837 11750 : DO i_local = 1, SIZE(indices_sub)
2838 11558 : i_global = indices_sub(i_local)
2839 11558 : proc = struct_global%g2p_col(i_global)
2840 11558 : data2send(proc) = data2send(proc) + 1
2841 11750 : index2send(proc)%array(data2send(proc)) = i_local
2842 : END DO
2843 :
2844 192 : END SUBROUTINE get_elements2send_col
2845 :
2846 : ! **************************************************************************************************
2847 : !> \brief ...
2848 : !> \param data2send ...
2849 : !> \param struct_global ...
2850 : !> \param indices_sub ...
2851 : !> \param index2send ...
2852 : ! **************************************************************************************************
2853 192 : SUBROUTINE get_elements2send_row(data2send, struct_global, indices_sub, index2send)
2854 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: data2send
2855 : TYPE(cp_fm_struct_type), INTENT(INOUT) :: struct_global
2856 : INTEGER, DIMENSION(:), INTENT(IN) :: indices_sub
2857 : TYPE(one_dim_int_array), ALLOCATABLE, &
2858 : DIMENSION(:), INTENT(OUT) :: index2send
2859 :
2860 : INTEGER :: i_global, i_local, np_global, proc
2861 :
2862 192 : CALL struct_global%context%get(number_of_process_rows=np_global)
2863 :
2864 576 : ALLOCATE (data2send(0:np_global - 1))
2865 192 : data2send = 0
2866 16332 : DO i_local = 1, SIZE(indices_sub)
2867 16140 : i_global = indices_sub(i_local)
2868 16140 : proc = struct_global%g2p_row(i_global)
2869 16332 : data2send(proc) = data2send(proc) + 1
2870 : END DO
2871 :
2872 880 : ALLOCATE (index2send(0:np_global - 1))
2873 496 : DO proc = 0, np_global - 1
2874 912 : ALLOCATE (index2send(proc)%array(data2send(proc)))
2875 : ! We want to crash if there is an error
2876 16636 : index2send(proc)%array = -1
2877 : END DO
2878 :
2879 192 : data2send = 0
2880 16332 : DO i_local = 1, SIZE(indices_sub)
2881 16140 : i_global = indices_sub(i_local)
2882 16140 : proc = struct_global%g2p_row(i_global)
2883 16140 : data2send(proc) = data2send(proc) + 1
2884 16332 : index2send(proc)%array(data2send(proc)) = i_local
2885 : END DO
2886 :
2887 192 : END SUBROUTINE get_elements2send_row
2888 :
2889 0 : END MODULE rpa_grad
|
