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