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