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