(git:98357aa)
Loading...
Searching...
No Matches
mp2_ri_gpw.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-GPW-MP2 energy using pw
10!> \par History
11!> 06.2012 created [Mauro Del Ben]
12!> 03.2019 Refactored from mp2_ri_gpw [Frederick Stein]
13! **************************************************************************************************
14MODULE mp2_ri_gpw
15 USE cp_log_handling, ONLY: cp_to_string
16 USE dgemm_counter_types, ONLY: dgemm_counter_init,&
17 dgemm_counter_start,&
18 dgemm_counter_stop,&
19 dgemm_counter_type,&
20 dgemm_counter_write
21 USE group_dist_types, ONLY: get_group_dist,&
22 group_dist_d1_type,&
23 maxsize,&
24 release_group_dist
25 USE kinds, ONLY: dp,&
26 int_8
27 USE libint_2c_3c, ONLY: compare_potential_types
28 USE local_gemm_api, ONLY: local_gemm_pu_gpu
29 USE machine, ONLY: m_flush,&
30 m_memory,&
31 m_walltime
32 USE message_passing, ONLY: mp_comm_type,&
33 mp_para_env_type
34 USE mp2_ri_grad_util, ONLY: complete_gamma
35 USE mp2_types, ONLY: mp2_type,&
36 three_dim_real_array
37
38!$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
39#include "./base/base_uses.f90"
40
41 IMPLICIT NONE
42
43 PRIVATE
44
45 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mp2_ri_gpw'
46
47 PUBLIC :: mp2_ri_gpw_compute_en
48
49CONTAINS
50
51! **************************************************************************************************
52!> \brief ...
53!> \param Emp2_Cou ...
54!> \param Emp2_EX ...
55!> \param Emp2_S ...
56!> \param Emp2_T ...
57!> \param BIb_C ...
58!> \param mp2_env ...
59!> \param para_env ...
60!> \param para_env_sub ...
61!> \param color_sub ...
62!> \param gd_array ...
63!> \param gd_B_virtual ...
64!> \param Eigenval ...
65!> \param nmo ...
66!> \param homo ...
67!> \param dimen_RI ...
68!> \param unit_nr ...
69!> \param calc_forces ...
70!> \param calc_ex ...
71! **************************************************************************************************
72 SUBROUTINE mp2_ri_gpw_compute_en(Emp2_Cou, Emp2_EX, Emp2_S, Emp2_T, BIb_C, mp2_env, para_env, para_env_sub, color_sub, &
73 gd_array, gd_B_virtual, &
74 Eigenval, nmo, homo, dimen_RI, unit_nr, calc_forces, calc_ex)
75 REAL(kind=dp), INTENT(INOUT) :: emp2_cou, emp2_ex, emp2_s, emp2_t
76 TYPE(three_dim_real_array), DIMENSION(:), &
77 INTENT(INOUT) :: bib_c
78 TYPE(mp2_type) :: mp2_env
79 TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env, para_env_sub
80 INTEGER, INTENT(IN) :: color_sub
81 TYPE(group_dist_d1_type), INTENT(INOUT) :: gd_array
82 INTEGER, DIMENSION(:), INTENT(IN) :: homo
83 INTEGER, INTENT(IN) :: nmo
84 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: eigenval
85 TYPE(group_dist_d1_type), DIMENSION(SIZE(homo)), &
86 INTENT(INOUT) :: gd_b_virtual
87 INTEGER, INTENT(IN) :: dimen_ri, unit_nr
88 LOGICAL, INTENT(IN) :: calc_forces, calc_ex
89
90 CHARACTER(LEN=*), PARAMETER :: routinen = 'mp2_ri_gpw_compute_en'
91
92 INTEGER :: a, a_global, b, b_global, block_size, decil, end_point, handle, handle2, handle3, &
93 iib, ij_counter, ij_counter_send, ij_index, integ_group_size, ispin, jjb, jspin, &
94 max_ij_pairs, my_block_size, my_group_l_end, my_group_l_size, my_group_l_size_orig, &
95 my_group_l_start, my_i, my_ij_pairs, my_j, my_new_group_l_size, ngroup, nspins, &
96 num_integ_group, proc_receive, proc_send, proc_shift, rec_b_size, rec_b_virtual_end, &
97 rec_b_virtual_start, rec_l_size, send_b_size, send_b_virtual_end, send_b_virtual_start, &
98 send_i, send_ij_index, send_j, start_point, tag, total_ij_pairs
99 INTEGER, ALLOCATABLE, DIMENSION(:) :: integ_group_pos2color_sub, my_b_size, &
100 my_b_virtual_end, my_b_virtual_start, num_ij_pairs, sizes_array_orig, virtual
101 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: ij_map
102 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: ranges_info_array
103 LOGICAL :: my_alpha_beta_case, my_beta_beta_case, &
104 my_open_shell_ss
105 REAL(kind=dp) :: amp_fac, my_emp2_cou, my_emp2_ex, &
106 sym_fac, t_new, t_start
107 REAL(kind=dp), ALLOCATABLE, DIMENSION(:), TARGET :: buffer_1d
108 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
109 TARGET :: local_ab, local_ba, t_ab
110 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
111 TARGET :: local_i_al, local_j_al, y_i_ap, y_j_ap
112 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
113 POINTER :: external_ab, external_i_al
114 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
115 POINTER :: bi_c_rec
116 TYPE(dgemm_counter_type) :: dgemm_counter
117 TYPE(mp_comm_type) :: comm_exchange, comm_rep
118 TYPE(three_dim_real_array), ALLOCATABLE, &
119 DIMENSION(:) :: b_ia_q
120
121 CALL timeset(routinen, handle)
122
123 nspins = SIZE(homo)
124
125 ALLOCATE (virtual(nspins))
126 virtual(:) = nmo - homo(:)
127
128 ALLOCATE (my_b_size(nspins), my_b_virtual_start(nspins), my_b_virtual_end(nspins))
129 DO ispin = 1, nspins
130 CALL get_group_dist(gd_b_virtual(ispin), para_env_sub%mepos, &
131 my_b_virtual_start(ispin), my_b_virtual_end(ispin), my_b_size(ispin))
132 END DO
133
134 CALL get_group_dist(gd_array, color_sub, my_group_l_start, my_group_l_end, my_group_l_size)
135
136 CALL dgemm_counter_init(dgemm_counter, unit_nr, mp2_env%ri_mp2%print_dgemm_info)
137
138 ! local_gemm_ctx has a very footprint the first time this routine is
139 ! called.
140 CALL mp2_env%local_gemm_ctx%create(local_gemm_pu_gpu)
141 CALL mp2_env%local_gemm_ctx%set_op_threshold_gpu(128*128*128*2)
142
143 CALL mp2_ri_get_integ_group_size( &
144 mp2_env, para_env, para_env_sub, gd_array, gd_b_virtual, &
145 homo, dimen_ri, unit_nr, &
146 integ_group_size, ngroup, &
147 num_integ_group, virtual, calc_forces)
148
149 ! now create a group that contains all the proc that have the same virtual starting point
150 ! in the integ group
151 CALL mp2_ri_create_group( &
152 para_env, para_env_sub, color_sub, &
153 gd_array%sizes, calc_forces, &
154 integ_group_size, my_group_l_end, &
155 my_group_l_size, my_group_l_size_orig, my_group_l_start, my_new_group_l_size, &
156 integ_group_pos2color_sub, sizes_array_orig, &
157 ranges_info_array, comm_exchange, comm_rep, num_integ_group)
158
159 ! We cannot fix the tag because of the recv routine
160 tag = 42
161
162 DO jspin = 1, nspins
163
164 CALL replicate_iak_2intgroup(bib_c(jspin)%array, comm_exchange, comm_rep, &
165 homo(jspin), gd_array%sizes, my_b_size(jspin), &
166 my_group_l_size, ranges_info_array)
167
168 DO ispin = 1, jspin
169
170 IF (unit_nr > 0) THEN
171 IF (nspins == 1) THEN
172 WRITE (unit_nr, *) "Start loop run"
173 ELSE IF (ispin == 1 .AND. jspin == 1) THEN
174 WRITE (unit_nr, *) "Start loop run alpha-alpha"
175 ELSE IF (ispin == 1 .AND. jspin == 2) THEN
176 WRITE (unit_nr, *) "Start loop run alpha-beta"
177 ELSE IF (ispin == 2 .AND. jspin == 2) THEN
178 WRITE (unit_nr, *) "Start loop run beta-beta"
179 END IF
180 CALL m_flush(unit_nr)
181 END IF
182
183 my_open_shell_ss = (nspins == 2) .AND. (ispin == jspin)
184
185 ! t_ab = amp_fac*(:,a|:,b)-(:,b|:,a)
186 ! If we calculate the gradient we need to distinguish
187 ! between alpha-alpha and beta-beta cases for UMP2
188
189 my_beta_beta_case = .false.
190 my_alpha_beta_case = .false.
191 IF (ispin /= jspin) THEN
192 my_alpha_beta_case = .true.
193 ELSE IF (my_open_shell_ss) THEN
194 IF (ispin == 2) my_beta_beta_case = .true.
195 END IF
196
197 amp_fac = mp2_env%scale_S + mp2_env%scale_T
198 IF (my_alpha_beta_case .OR. my_open_shell_ss) amp_fac = mp2_env%scale_T
199
200 CALL mp2_ri_allocate_no_blk(local_ab, t_ab, mp2_env, homo, virtual, my_b_size, &
201 my_group_l_size, calc_forces, ispin, jspin, local_ba)
202
203 CALL mp2_ri_get_block_size( &
204 mp2_env, para_env, para_env_sub, gd_array, gd_b_virtual(ispin:jspin), &
205 homo(ispin:jspin), virtual(ispin:jspin), dimen_ri, unit_nr, block_size, &
206 ngroup, num_integ_group, my_open_shell_ss, calc_forces, buffer_1d)
207
208 ! *****************************************************************
209 ! ********** REPLICATION-BLOCKED COMMUNICATION SCHEME ***********
210 ! *****************************************************************
211 ! introduce block size, the number of occupied orbitals has to be a
212 ! multiple of the block size
213
214 ! Calculate the maximum number of ij pairs that have to be computed
215 ! among groups
216 CALL mp2_ri_communication(my_alpha_beta_case, total_ij_pairs, homo(ispin), homo(jspin), &
217 block_size, ngroup, ij_map, color_sub, my_ij_pairs, my_open_shell_ss, unit_nr)
218
219 ALLOCATE (num_ij_pairs(0:comm_exchange%num_pe - 1))
220 CALL comm_exchange%allgather(my_ij_pairs, num_ij_pairs)
221
222 max_ij_pairs = maxval(num_ij_pairs)
223
224 ! start real stuff
225 CALL mp2_ri_allocate_blk(dimen_ri, my_b_size, block_size, local_i_al, &
226 local_j_al, calc_forces, y_i_ap, y_j_ap, ispin, jspin)
227
228 CALL timeset(routinen//"_RI_loop", handle2)
229 my_emp2_cou = 0.0_dp
230 my_emp2_ex = 0.0_dp
231 t_start = m_walltime()
232 DO ij_index = 1, max_ij_pairs
233
234 ! Prediction is unreliable if we are in the first step of the loop
235 IF (unit_nr > 0 .AND. ij_index > 1) THEN
236 decil = ij_index*10/max_ij_pairs
237 IF (decil /= (ij_index - 1)*10/max_ij_pairs) THEN
238 t_new = m_walltime()
239 t_new = (t_new - t_start)/60.0_dp*(max_ij_pairs - ij_index + 1)/(ij_index - 1)
240 WRITE (unit_nr, fmt="(T3,A)") "Percentage of finished loop: "// &
241 cp_to_string(decil*10)//". Minutes left: "//cp_to_string(t_new)
242 CALL m_flush(unit_nr)
243 END IF
244 END IF
245
246 IF (calc_forces) THEN
247 y_i_ap = 0.0_dp
248 y_j_ap = 0.0_dp
249 END IF
250
251 IF (ij_index <= my_ij_pairs) THEN
252 ! We have work to do
253 ij_counter = (ij_index - min(1, color_sub))*ngroup + color_sub
254 my_i = ij_map(1, ij_counter)
255 my_j = ij_map(2, ij_counter)
256 my_block_size = ij_map(3, ij_counter)
257
258 local_i_al = 0.0_dp
259 CALL fill_local_i_al(local_i_al(:, :, 1:my_block_size), ranges_info_array(:, :, comm_exchange%mepos), &
260 bib_c(ispin)%array(:, :, my_i:my_i + my_block_size - 1))
261
262 local_j_al = 0.0_dp
263 CALL fill_local_i_al(local_j_al(:, :, 1:my_block_size), ranges_info_array(:, :, comm_exchange%mepos), &
264 bib_c(jspin)%array(:, :, my_j:my_j + my_block_size - 1))
265
266 ! collect data from other proc
267 CALL timeset(routinen//"_comm", handle3)
268 DO proc_shift = 1, comm_exchange%num_pe - 1
269 proc_send = modulo(comm_exchange%mepos + proc_shift, comm_exchange%num_pe)
270 proc_receive = modulo(comm_exchange%mepos - proc_shift, comm_exchange%num_pe)
271
272 send_ij_index = num_ij_pairs(proc_send)
273
274 CALL get_group_dist(gd_array, proc_receive, sizes=rec_l_size)
275
276 IF (ij_index <= send_ij_index) THEN
277 ij_counter_send = (ij_index - min(1, integ_group_pos2color_sub(proc_send)))*ngroup + &
278 integ_group_pos2color_sub(proc_send)
279 send_i = ij_map(1, ij_counter_send)
280 send_j = ij_map(2, ij_counter_send)
281
282 ! occupied i
283 bi_c_rec(1:rec_l_size, 1:my_b_size(ispin), 1:my_block_size) => &
284 buffer_1d(1:rec_l_size*my_b_size(ispin)*my_block_size)
285 bi_c_rec = 0.0_dp
286 CALL comm_exchange%sendrecv(bib_c(ispin)%array(:, :, send_i:send_i + my_block_size - 1), &
287 proc_send, bi_c_rec, proc_receive, tag)
288
289 CALL fill_local_i_al(local_i_al(:, :, 1:my_block_size), ranges_info_array(:, :, proc_receive), &
290 bi_c_rec(:, 1:my_b_size(ispin), :))
291
292 ! occupied j
293 bi_c_rec(1:rec_l_size, 1:my_b_size(jspin), 1:my_block_size) => &
294 buffer_1d(1:int(rec_l_size, int_8)*my_b_size(jspin)*my_block_size)
295 bi_c_rec = 0.0_dp
296 CALL comm_exchange%sendrecv(bib_c(jspin)%array(:, :, send_j:send_j + my_block_size - 1), &
297 proc_send, bi_c_rec, proc_receive, tag)
298
299 CALL fill_local_i_al(local_j_al(:, :, 1:my_block_size), ranges_info_array(:, :, proc_receive), &
300 bi_c_rec(:, 1:my_b_size(jspin), :))
301
302 ELSE
303 ! we send nothing while we know that we have to receive something
304
305 ! occupied i
306 bi_c_rec(1:rec_l_size, 1:my_b_size(ispin), 1:my_block_size) => &
307 buffer_1d(1:int(rec_l_size, int_8)*my_b_size(ispin)*my_block_size)
308 bi_c_rec = 0.0_dp
309 CALL comm_exchange%recv(bi_c_rec, proc_receive, tag)
310
311 CALL fill_local_i_al(local_i_al(:, :, 1:my_block_size), ranges_info_array(:, :, proc_receive), &
312 bi_c_rec(:, 1:my_b_size(ispin), 1:my_block_size))
313
314 ! occupied j
315 bi_c_rec(1:rec_l_size, 1:my_b_size(jspin), 1:my_block_size) => &
316 buffer_1d(1:int(rec_l_size, int_8)*my_b_size(jspin)*my_block_size)
317 bi_c_rec = 0.0_dp
318 CALL comm_exchange%recv(bi_c_rec, proc_receive, tag)
319
320 CALL fill_local_i_al(local_j_al(:, :, 1:my_block_size), ranges_info_array(:, :, proc_receive), &
321 bi_c_rec(:, 1:my_b_size(jspin), 1:my_block_size))
322
323 END IF
324
325 END DO
326
327 CALL timestop(handle3)
328
329 ! loop over the block elements
330 DO iib = 1, my_block_size
331 DO jjb = 1, my_block_size
332 CALL timeset(routinen//"_expansion", handle3)
333 associate(my_local_i_al => local_i_al(:, :, iib), my_local_j_al => local_j_al(:, :, jjb))
334
335 ! calculate the integrals (ia|jb) strating from my local data ...
336 local_ab = 0.0_dp
337 IF ((my_alpha_beta_case) .AND. (calc_forces)) THEN
338 local_ba = 0.0_dp
339 END IF
340 CALL dgemm_counter_start(dgemm_counter)
341 CALL mp2_env%local_gemm_ctx%gemm('T', 'N', my_b_size(ispin), my_b_size(jspin), dimen_ri, 1.0_dp, &
342 my_local_i_al, dimen_ri, my_local_j_al, dimen_ri, &
343 0.0_dp, local_ab(my_b_virtual_start(ispin):my_b_virtual_end(ispin), :), &
344 my_b_size(ispin))
345 ! Additional integrals only for alpha_beta case and forces
346 IF (my_alpha_beta_case .AND. calc_forces) THEN
347 local_ba(my_b_virtual_start(jspin):my_b_virtual_end(jspin), :) = &
348 transpose(local_ab(my_b_virtual_start(ispin):my_b_virtual_end(ispin), :))
349 END IF
350 ! ... and from the other of my subgroup
351 DO proc_shift = 1, para_env_sub%num_pe - 1
352 proc_send = modulo(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
353 proc_receive = modulo(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
354
355 CALL get_group_dist(gd_b_virtual(ispin), proc_receive, rec_b_virtual_start, &
356 rec_b_virtual_end, rec_b_size)
357
358 external_i_al(1:dimen_ri, 1:rec_b_size) => buffer_1d(1:int(dimen_ri, int_8)*rec_b_size)
359 external_i_al = 0.0_dp
360
361 CALL para_env_sub%sendrecv(my_local_i_al, proc_send, &
362 external_i_al, proc_receive, tag)
363
364 CALL mp2_env%local_gemm_ctx%gemm( &
365 'T', 'N', rec_b_size, my_b_size(jspin), dimen_ri, 1.0_dp, &
366 external_i_al, dimen_ri, my_local_j_al, dimen_ri, &
367 0.0_dp, local_ab(rec_b_virtual_start:rec_b_virtual_end, 1:my_b_size(jspin)), rec_b_size)
368
369 ! Additional integrals only for alpha_beta case and forces
370 IF (my_alpha_beta_case .AND. calc_forces) THEN
371
372 CALL get_group_dist(gd_b_virtual(jspin), proc_receive, rec_b_virtual_start, &
373 rec_b_virtual_end, rec_b_size)
374
375 external_i_al(1:dimen_ri, 1:rec_b_size) => buffer_1d(1:int(dimen_ri, int_8)*rec_b_size)
376 external_i_al = 0.0_dp
377
378 CALL para_env_sub%sendrecv(my_local_j_al, proc_send, &
379 external_i_al, proc_receive, tag)
380
381 CALL mp2_env%local_gemm_ctx%gemm('T', 'N', rec_b_size, my_b_size(ispin), dimen_ri, 1.0_dp, &
382 external_i_al, dimen_ri, my_local_i_al, dimen_ri, &
383 0.0_dp, local_ba(rec_b_virtual_start:rec_b_virtual_end, 1:my_b_size(ispin)), rec_b_size)
384 END IF
385 END DO
386 IF (my_alpha_beta_case .AND. calc_forces) THEN
387 ! Is just an approximation, but the call does not allow it,
388 ! it ought to be (virtual_i*B_size_j+virtual_j*B_size_i)*dimen_RI
389 CALL dgemm_counter_stop(dgemm_counter, virtual(ispin), my_b_size(ispin) + my_b_size(jspin), dimen_ri)
390 ELSE
391 CALL dgemm_counter_stop(dgemm_counter, virtual(ispin), my_b_size(jspin), dimen_ri)
392 END IF
393 CALL timestop(handle3)
394
395 !sample peak memory
396 CALL m_memory()
397
398 CALL timeset(routinen//"_ener", handle3)
399 ! calculate coulomb only MP2
400 sym_fac = 2.0_dp
401 IF (my_i == my_j) sym_fac = 1.0_dp
402 IF (my_alpha_beta_case) sym_fac = 0.5_dp
403 DO b = 1, my_b_size(jspin)
404 b_global = b + my_b_virtual_start(jspin) - 1
405 DO a = 1, virtual(ispin)
406 my_emp2_cou = my_emp2_cou - sym_fac*2.0_dp*local_ab(a, b)**2/ &
407 (eigenval(homo(ispin) + a, ispin) + eigenval(homo(jspin) + b_global, jspin) - &
408 eigenval(my_i + iib - 1, ispin) - eigenval(my_j + jjb - 1, jspin))
409 END DO
410 END DO
411
412 IF (calc_ex) THEN
413 ! contract integrals with orbital energies for exchange MP2 energy
414 ! starting with local ...
415 IF (calc_forces .AND. (.NOT. my_alpha_beta_case)) t_ab = 0.0_dp
416 DO b = 1, my_b_size(ispin)
417 b_global = b + my_b_virtual_start(ispin) - 1
418 DO a = 1, my_b_size(ispin)
419 a_global = a + my_b_virtual_start(ispin) - 1
420 my_emp2_ex = my_emp2_ex + sym_fac*local_ab(a_global, b)*local_ab(b_global, a)/ &
421 (eigenval(homo(ispin) + a_global, ispin) + eigenval(homo(ispin) + b_global, ispin) - &
422 eigenval(my_i + iib - 1, ispin) - eigenval(my_j + jjb - 1, ispin))
423 IF (calc_forces .AND. (.NOT. my_alpha_beta_case)) THEN
424 t_ab(a_global, b) = -(amp_fac*local_ab(a_global, b) - mp2_env%scale_T*local_ab(b_global, a))/ &
425 (eigenval(homo(ispin) + a_global, ispin) + &
426 eigenval(homo(ispin) + b_global, ispin) - &
427 eigenval(my_i + iib - 1, ispin) - eigenval(my_j + jjb - 1, ispin))
428 END IF
429 END DO
430 END DO
431 ! ... and then with external data
432 DO proc_shift = 1, para_env_sub%num_pe - 1
433 proc_send = modulo(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
434 proc_receive = modulo(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
435
436 CALL get_group_dist(gd_b_virtual(ispin), proc_receive, &
437 rec_b_virtual_start, rec_b_virtual_end, rec_b_size)
438 CALL get_group_dist(gd_b_virtual(ispin), proc_send, &
439 send_b_virtual_start, send_b_virtual_end, send_b_size)
440
441 external_ab(1:my_b_size(ispin), 1:rec_b_size) => &
442 buffer_1d(1:int(rec_b_size, int_8)*my_b_size(ispin))
443 external_ab = 0.0_dp
444
445 CALL para_env_sub%sendrecv(local_ab(send_b_virtual_start:send_b_virtual_end, 1:my_b_size(ispin)), proc_send, &
446 external_ab(1:my_b_size(ispin), 1:rec_b_size), proc_receive, tag)
447
448 DO b = 1, my_b_size(ispin)
449 b_global = b + my_b_virtual_start(ispin) - 1
450 DO a = 1, rec_b_size
451 a_global = a + rec_b_virtual_start - 1
452 my_emp2_ex = my_emp2_ex + sym_fac*local_ab(a_global, b)*external_ab(b, a)/ &
453 (eigenval(homo(ispin) + a_global, ispin) + eigenval(homo(ispin) + b_global, ispin) - &
454 eigenval(my_i + iib - 1, ispin) - eigenval(my_j + jjb - 1, ispin))
455 IF (calc_forces .AND. (.NOT. my_alpha_beta_case)) THEN
456 t_ab(a_global, b) = -(amp_fac*local_ab(a_global, b) - mp2_env%scale_T*external_ab(b, a))/ &
457 (eigenval(homo(ispin) + a_global, ispin) + &
458 eigenval(homo(ispin) + b_global, ispin) - &
459 eigenval(my_i + iib - 1, ispin) - eigenval(my_j + jjb - 1, ispin))
460 END IF
461 END DO
462 END DO
463 END DO
464 END IF
465 CALL timestop(handle3)
466
467 IF (calc_forces) THEN
468 ! update P_ab, Gamma_P_ia
469 CALL mp2_update_p_gamma(mp2_env, para_env_sub, gd_b_virtual, &
470 eigenval, homo, dimen_ri, iib, jjb, my_b_size, &
471 my_b_virtual_end, my_b_virtual_start, my_i, my_j, virtual, &
472 local_ab, t_ab, my_local_i_al, my_local_j_al, &
473 my_open_shell_ss, y_i_ap(:, :, iib), y_j_ap(:, :, jjb), local_ba, &
474 ispin, jspin, dgemm_counter, buffer_1d)
475
476 END IF
477
478 END associate
479
480 END DO ! jjB
481 END DO ! iiB
482
483 ELSE
484 ! We need it later in case of gradients
485 my_block_size = 1
486
487 CALL timeset(routinen//"_comm", handle3)
488 ! No work to do and we know that we have to receive nothing, but send something
489 ! send data to other proc
490 DO proc_shift = 1, comm_exchange%num_pe - 1
491 proc_send = modulo(comm_exchange%mepos + proc_shift, comm_exchange%num_pe)
492 proc_receive = modulo(comm_exchange%mepos - proc_shift, comm_exchange%num_pe)
493
494 send_ij_index = num_ij_pairs(proc_send)
495
496 IF (ij_index <= send_ij_index) THEN
497 ! something to send
498 ij_counter_send = (ij_index - min(1, integ_group_pos2color_sub(proc_send)))*ngroup + &
499 integ_group_pos2color_sub(proc_send)
500 send_i = ij_map(1, ij_counter_send)
501 send_j = ij_map(2, ij_counter_send)
502
503 ! occupied i
504 CALL comm_exchange%send(bib_c(ispin)%array(:, :, send_i:send_i + my_block_size - 1), &
505 proc_send, tag)
506 ! occupied j
507 CALL comm_exchange%send(bib_c(jspin)%array(:, :, send_j:send_j + my_block_size - 1), &
508 proc_send, tag)
509 END IF
510 END DO
511 CALL timestop(handle3)
512 END IF
513
514 ! redistribute gamma
515 IF (calc_forces) THEN
516 CALL mp2_redistribute_gamma(mp2_env%ri_grad%Gamma_P_ia(ispin)%array, ij_index, my_b_size(ispin), &
517 my_block_size, my_group_l_size, my_i, my_ij_pairs, ngroup, &
518 num_integ_group, integ_group_pos2color_sub, num_ij_pairs, &
519 ij_map, ranges_info_array, y_i_ap(:, :, 1:my_block_size), comm_exchange, &
520 gd_array%sizes, 1, buffer_1d)
521 CALL mp2_redistribute_gamma(mp2_env%ri_grad%Gamma_P_ia(jspin)%array, ij_index, my_b_size(jspin), &
522 my_block_size, my_group_l_size, my_j, my_ij_pairs, ngroup, &
523 num_integ_group, integ_group_pos2color_sub, num_ij_pairs, &
524 ij_map, ranges_info_array, y_j_ap(:, :, 1:my_block_size), comm_exchange, &
525 gd_array%sizes, 2, buffer_1d)
526 END IF
527
528 END DO
529 CALL timestop(handle2)
530
531 DEALLOCATE (local_i_al)
532 DEALLOCATE (local_j_al)
533 DEALLOCATE (ij_map)
534 DEALLOCATE (num_ij_pairs)
535 DEALLOCATE (local_ab)
536
537 IF (calc_forces) THEN
538 DEALLOCATE (y_i_ap)
539 DEALLOCATE (y_j_ap)
540 IF (ALLOCATED(t_ab)) THEN
541 DEALLOCATE (t_ab)
542 END IF
543 DEALLOCATE (local_ba)
544
545 ! here we check if there are almost degenerate ij
546 ! pairs and we update P_ij with these contribution.
547 ! If all pairs are degenerate with each other this step will scale O(N^6),
548 ! if the number of degenerate pairs scales linearly with the system size
549 ! this step will scale O(N^5).
550 ! Start counting the number of almost degenerate ij pairs according
551 ! to eps_canonical
552 CALL quasi_degenerate_p_ij( &
553 mp2_env, eigenval(:, ispin:jspin), homo(ispin:jspin), virtual(ispin:jspin), my_open_shell_ss, &
554 my_beta_beta_case, bib_c(ispin:jspin), unit_nr, dimen_ri, &
555 my_b_size(ispin:jspin), ngroup, my_group_l_size, &
556 color_sub, ranges_info_array, comm_exchange, para_env_sub, para_env, &
557 my_b_virtual_start(ispin:jspin), my_b_virtual_end(ispin:jspin), gd_array%sizes, gd_b_virtual(ispin:jspin), &
558 integ_group_pos2color_sub, dgemm_counter, buffer_1d)
559
560 END IF
561
562 DEALLOCATE (buffer_1d)
563
564 ! Dereplicate BIb_C and Gamma_P_ia to save memory
565 ! These matrices will not be needed in that fashion anymore
566 ! B_ia_Q will needed later
567 IF (calc_forces .AND. jspin == nspins) THEN
568 IF (.NOT. ALLOCATED(b_ia_q)) ALLOCATE (b_ia_q(nspins))
569 ALLOCATE (b_ia_q(ispin)%array(homo(ispin), my_b_size(ispin), my_group_l_size_orig))
570 b_ia_q(ispin)%array = 0.0_dp
571 DO jjb = 1, homo(ispin)
572 DO iib = 1, my_b_size(ispin)
573 b_ia_q(ispin)%array(jjb, iib, 1:my_group_l_size_orig) = &
574 bib_c(ispin)%array(1:my_group_l_size_orig, iib, jjb)
575 END DO
576 END DO
577 DEALLOCATE (bib_c(ispin)%array)
578
579 ! sum Gamma and dereplicate
580 ALLOCATE (bib_c(ispin)%array(my_b_size(ispin), homo(ispin), my_group_l_size_orig))
581 DO proc_shift = 1, comm_rep%num_pe - 1
582 ! invert order
583 proc_send = modulo(comm_rep%mepos - proc_shift, comm_rep%num_pe)
584 proc_receive = modulo(comm_rep%mepos + proc_shift, comm_rep%num_pe)
585
586 start_point = ranges_info_array(3, proc_shift, comm_exchange%mepos)
587 end_point = ranges_info_array(4, proc_shift, comm_exchange%mepos)
588
589 CALL comm_rep%sendrecv(mp2_env%ri_grad%Gamma_P_ia(ispin)%array(:, :, start_point:end_point), &
590 proc_send, bib_c(ispin)%array, proc_receive, tag)
591!$OMP PARALLEL WORKSHARE DEFAULT(NONE) &
592!$OMP SHARED(mp2_env,BIb_C,ispin,homo,my_B_size,my_group_L_size_orig)
593 mp2_env%ri_grad%Gamma_P_ia(ispin)%array(:, :, 1:my_group_l_size_orig) = &
594 mp2_env%ri_grad%Gamma_P_ia(ispin)%array(:, :, 1:my_group_l_size_orig) &
595 + bib_c(ispin)%array(:, :, :)
596!$OMP END PARALLEL WORKSHARE
597 END DO
598
599 bib_c(ispin)%array(:, :, :) = mp2_env%ri_grad%Gamma_P_ia(ispin)%array(:, :, 1:my_group_l_size_orig)
600 DEALLOCATE (mp2_env%ri_grad%Gamma_P_ia(ispin)%array)
601 CALL move_alloc(bib_c(ispin)%array, mp2_env%ri_grad%Gamma_P_ia(ispin)%array)
602 ELSE IF (jspin == nspins) THEN
603 DEALLOCATE (bib_c(ispin)%array)
604 END IF
605
606 CALL para_env%sum(my_emp2_cou)
607 CALL para_env%sum(my_emp2_ex)
608
609 IF (my_open_shell_ss .OR. my_alpha_beta_case) THEN
610 IF (my_alpha_beta_case) THEN
611 emp2_s = emp2_s + my_emp2_cou
612 emp2_cou = emp2_cou + my_emp2_cou
613 ELSE
614 my_emp2_cou = my_emp2_cou*0.25_dp
615 my_emp2_ex = my_emp2_ex*0.5_dp
616 emp2_t = emp2_t + my_emp2_cou + my_emp2_ex
617 emp2_cou = emp2_cou + my_emp2_cou
618 emp2_ex = emp2_ex + my_emp2_ex
619 END IF
620 ELSE
621 emp2_cou = emp2_cou + my_emp2_cou
622 emp2_ex = emp2_ex + my_emp2_ex
623 END IF
624 END DO
625
626 END DO
627
628 DEALLOCATE (integ_group_pos2color_sub)
629 DEALLOCATE (ranges_info_array)
630
631 CALL comm_exchange%free()
632 CALL comm_rep%free()
633
634 IF (calc_forces) THEN
635 ! recover original information (before replication)
636 DEALLOCATE (gd_array%sizes)
637 iib = SIZE(sizes_array_orig)
638 ALLOCATE (gd_array%sizes(0:iib - 1))
639 gd_array%sizes(:) = sizes_array_orig
640 DEALLOCATE (sizes_array_orig)
641
642 ! Remove replication from BIb_C and reorder the matrix
643 my_group_l_size = my_group_l_size_orig
644
645 ! B_ia_Q(ispin)%array will be deallocated inside of complete_gamma
646 DO ispin = 1, nspins
647 CALL complete_gamma(mp2_env, b_ia_q(ispin)%array, dimen_ri, homo(ispin), &
648 virtual(ispin), para_env, para_env_sub, ngroup, &
649 my_group_l_size, my_group_l_start, my_group_l_end, &
650 my_b_size(ispin), my_b_virtual_start(ispin), &
651 gd_array, gd_b_virtual(ispin), &
652 ispin)
653 END DO
654 DEALLOCATE (b_ia_q)
655
656 IF (nspins == 1) mp2_env%ri_grad%P_ab(1)%array(:, :) = mp2_env%ri_grad%P_ab(1)%array(:, :)*2.0_dp
657 block
658 TYPE(mp_comm_type) :: comm
659 CALL comm%from_split(para_env, para_env_sub%mepos)
660 DO ispin = 1, nspins
661 ! P_ab is only replicated over all subgroups
662 CALL comm%sum(mp2_env%ri_grad%P_ab(ispin)%array)
663 ! P_ij is replicated over all processes
664 CALL para_env%sum(mp2_env%ri_grad%P_ij(ispin)%array)
665 END DO
666 CALL comm%free()
667 END block
668 END IF
669
670 CALL release_group_dist(gd_array)
671 DO ispin = 1, nspins
672 IF (ALLOCATED(bib_c(ispin)%array)) DEALLOCATE (bib_c(ispin)%array)
673 CALL release_group_dist(gd_b_virtual(ispin))
674 END DO
675
676 ! We do not need this matrix later, so deallocate it here to safe memory
677 IF (calc_forces) DEALLOCATE (mp2_env%ri_grad%PQ_half)
678 IF (calc_forces .AND. .NOT. compare_potential_types(mp2_env%ri_metric, mp2_env%potential_parameter)) THEN
679 DEALLOCATE (mp2_env%ri_grad%operator_half)
680 END IF
681
682 CALL dgemm_counter_write(dgemm_counter, para_env)
683
684 ! release memory allocated by local_gemm when run on GPU. local_gemm_ctx is null on cpu only runs
685 CALL mp2_env%local_gemm_ctx%destroy()
686 CALL timestop(handle)
687
688 END SUBROUTINE mp2_ri_gpw_compute_en
689
690! **************************************************************************************************
691!> \brief ...
692!> \param local_i_aL ...
693!> \param ranges_info_array ...
694!> \param BIb_C_rec ...
695! **************************************************************************************************
696 SUBROUTINE fill_local_i_al(local_i_aL, ranges_info_array, BIb_C_rec)
697 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: local_i_al
698 INTEGER, DIMENSION(:, :), INTENT(IN) :: ranges_info_array
699 REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: bib_c_rec
700
701 CHARACTER(LEN=*), PARAMETER :: routinen = 'fill_local_i_aL'
702
703 INTEGER :: end_point, handle, irep, lend_pos, &
704 lstart_pos, start_point
705
706 CALL timeset(routinen, handle)
707
708 DO irep = 1, SIZE(ranges_info_array, 2)
709 lstart_pos = ranges_info_array(1, irep)
710 lend_pos = ranges_info_array(2, irep)
711 start_point = ranges_info_array(3, irep)
712 end_point = ranges_info_array(4, irep)
713
714!$OMP PARALLEL WORKSHARE DEFAULT(NONE) &
715!$OMP SHARED(BIb_C_rec,local_i_aL,Lstart_pos,Lend_pos,start_point,end_point)
716 local_i_al(lstart_pos:lend_pos, :, :) = bib_c_rec(start_point:end_point, :, :)
717!$OMP END PARALLEL WORKSHARE
718 END DO
719
720 CALL timestop(handle)
721
722 END SUBROUTINE fill_local_i_al
723
724! **************************************************************************************************
725!> \brief ...
726!> \param local_i_aL ...
727!> \param ranges_info_array ...
728!> \param BIb_C_rec ...
729! **************************************************************************************************
730 SUBROUTINE fill_local_i_al_2d(local_i_aL, ranges_info_array, BIb_C_rec)
731 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: local_i_al
732 INTEGER, DIMENSION(:, :), INTENT(IN) :: ranges_info_array
733 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: bib_c_rec
734
735 CHARACTER(LEN=*), PARAMETER :: routinen = 'fill_local_i_aL_2D'
736
737 INTEGER :: end_point, handle, irep, lend_pos, &
738 lstart_pos, start_point
739
740 CALL timeset(routinen, handle)
741
742 DO irep = 1, SIZE(ranges_info_array, 2)
743 lstart_pos = ranges_info_array(1, irep)
744 lend_pos = ranges_info_array(2, irep)
745 start_point = ranges_info_array(3, irep)
746 end_point = ranges_info_array(4, irep)
747
748!$OMP PARALLEL WORKSHARE DEFAULT(NONE) &
749!$OMP SHARED(BIb_C_rec,local_i_aL,Lstart_pos,Lend_pos,start_point,end_point)
750 local_i_al(lstart_pos:lend_pos, :) = bib_c_rec(start_point:end_point, :)
751!$OMP END PARALLEL WORKSHARE
752 END DO
753
754 CALL timestop(handle)
755
756 END SUBROUTINE fill_local_i_al_2d
757
758! **************************************************************************************************
759!> \brief ...
760!> \param BIb_C ...
761!> \param comm_exchange ...
762!> \param comm_rep ...
763!> \param homo ...
764!> \param sizes_array ...
765!> \param my_B_size ...
766!> \param my_group_L_size ...
767!> \param ranges_info_array ...
768! **************************************************************************************************
769 SUBROUTINE replicate_iak_2intgroup(BIb_C, comm_exchange, comm_rep, homo, sizes_array, my_B_size, &
770 my_group_L_size, ranges_info_array)
771 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
772 INTENT(INOUT) :: bib_c
773 TYPE(mp_comm_type), INTENT(IN) :: comm_exchange, comm_rep
774 INTEGER, INTENT(IN) :: homo
775 INTEGER, DIMENSION(:), INTENT(IN) :: sizes_array
776 INTEGER, INTENT(IN) :: my_b_size, my_group_l_size
777 INTEGER, DIMENSION(:, 0:, 0:), INTENT(IN) :: ranges_info_array
778
779 CHARACTER(LEN=*), PARAMETER :: routinen = 'replicate_iaK_2intgroup'
780
781 INTEGER :: end_point, handle, max_l_size, &
782 proc_receive, proc_shift, start_point
783 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: bib_c_copy
784 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: bib_c_gather
785
786 CALL timeset(routinen, handle)
787
788 ! replication scheme using mpi_allgather
789 ! get the max L size of the
790 max_l_size = maxval(sizes_array)
791
792 ALLOCATE (bib_c_copy(max_l_size, my_b_size, homo))
793 bib_c_copy = 0.0_dp
794 bib_c_copy(1:SIZE(bib_c, 1), 1:my_b_size, 1:homo) = bib_c
795
796 DEALLOCATE (bib_c)
797
798 ALLOCATE (bib_c_gather(max_l_size, my_b_size, homo, 0:comm_rep%num_pe - 1))
799 bib_c_gather = 0.0_dp
800
801 CALL comm_rep%allgather(bib_c_copy, bib_c_gather)
802
803 DEALLOCATE (bib_c_copy)
804
805 ALLOCATE (bib_c(my_group_l_size, my_b_size, homo))
806 bib_c = 0.0_dp
807
808 ! reorder data
809 DO proc_shift = 0, comm_rep%num_pe - 1
810 proc_receive = modulo(comm_rep%mepos - proc_shift, comm_rep%num_pe)
811
812 start_point = ranges_info_array(3, proc_shift, comm_exchange%mepos)
813 end_point = ranges_info_array(4, proc_shift, comm_exchange%mepos)
814
815 bib_c(start_point:end_point, 1:my_b_size, 1:homo) = &
816 bib_c_gather(1:end_point - start_point + 1, 1:my_b_size, 1:homo, proc_receive)
817
818 END DO
819
820 DEALLOCATE (bib_c_gather)
821
822 CALL timestop(handle)
823
824 END SUBROUTINE replicate_iak_2intgroup
825
826! **************************************************************************************************
827!> \brief ...
828!> \param local_ab ...
829!> \param t_ab ...
830!> \param mp2_env ...
831!> \param homo ...
832!> \param virtual ...
833!> \param my_B_size ...
834!> \param my_group_L_size ...
835!> \param calc_forces ...
836!> \param ispin ...
837!> \param jspin ...
838!> \param local_ba ...
839! **************************************************************************************************
840 SUBROUTINE mp2_ri_allocate_no_blk(local_ab, t_ab, mp2_env, homo, virtual, my_B_size, &
841 my_group_L_size, calc_forces, ispin, jspin, local_ba)
842 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
843 INTENT(OUT) :: local_ab, t_ab
844 TYPE(mp2_type) :: mp2_env
845 INTEGER, INTENT(IN) :: homo(2), virtual(2), my_b_size(2), &
846 my_group_l_size
847 LOGICAL, INTENT(IN) :: calc_forces
848 INTEGER, INTENT(IN) :: ispin, jspin
849 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
850 INTENT(OUT) :: local_ba
851
852 CHARACTER(LEN=*), PARAMETER :: routinen = 'mp2_ri_allocate_no_blk'
853
854 INTEGER :: handle
855
856 CALL timeset(routinen, handle)
857
858 ALLOCATE (local_ab(virtual(ispin), my_b_size(jspin)))
859 local_ab = 0.0_dp
860
861 IF (calc_forces) THEN
862 IF (.NOT. ALLOCATED(mp2_env%ri_grad%P_ij(jspin)%array)) THEN
863 ALLOCATE (mp2_env%ri_grad%P_ij(jspin)%array(homo(ispin), homo(ispin)))
864 mp2_env%ri_grad%P_ij(jspin)%array = 0.0_dp
865 END IF
866 IF (.NOT. ALLOCATED(mp2_env%ri_grad%P_ab(jspin)%array)) THEN
867 ALLOCATE (mp2_env%ri_grad%P_ab(jspin)%array(my_b_size(jspin), virtual(jspin)))
868 mp2_env%ri_grad%P_ab(jspin)%array = 0.0_dp
869 END IF
870 IF (.NOT. ALLOCATED(mp2_env%ri_grad%Gamma_P_ia(jspin)%array)) THEN
871 ALLOCATE (mp2_env%ri_grad%Gamma_P_ia(jspin)%array(my_b_size(jspin), homo(jspin), my_group_l_size))
872 mp2_env%ri_grad%Gamma_P_ia(jspin)%array = 0.0_dp
873 END IF
874
875 IF (ispin == jspin) THEN
876 ! For non-alpha-beta case we need amplitudes
877 ALLOCATE (t_ab(virtual(ispin), my_b_size(jspin)))
878
879 ! That is just a dummy. In that way, we can pass it as array to other routines w/o requirement for allocatable array
880 ALLOCATE (local_ba(1, 1))
881 ELSE
882 ! We need more integrals
883 ALLOCATE (local_ba(virtual(jspin), my_b_size(ispin)))
884 END IF
885 END IF
886 !
887
888 CALL timestop(handle)
889
890 END SUBROUTINE mp2_ri_allocate_no_blk
891
892! **************************************************************************************************
893!> \brief ...
894!> \param dimen_RI ...
895!> \param my_B_size ...
896!> \param block_size ...
897!> \param local_i_aL ...
898!> \param local_j_aL ...
899!> \param calc_forces ...
900!> \param Y_i_aP ...
901!> \param Y_j_aP ...
902!> \param ispin ...
903!> \param jspin ...
904! **************************************************************************************************
905 SUBROUTINE mp2_ri_allocate_blk(dimen_RI, my_B_size, block_size, &
906 local_i_aL, local_j_aL, calc_forces, &
907 Y_i_aP, Y_j_aP, ispin, jspin)
908 INTEGER, INTENT(IN) :: dimen_ri, my_b_size(2), block_size
909 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
910 INTENT(OUT) :: local_i_al, local_j_al
911 LOGICAL, INTENT(IN) :: calc_forces
912 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
913 INTENT(OUT) :: y_i_ap, y_j_ap
914 INTEGER, INTENT(IN) :: ispin, jspin
915
916 CHARACTER(LEN=*), PARAMETER :: routinen = 'mp2_ri_allocate_blk'
917
918 INTEGER :: handle
919
920 CALL timeset(routinen, handle)
921
922 ALLOCATE (local_i_al(dimen_ri, my_b_size(ispin), block_size))
923 local_i_al = 0.0_dp
924 ALLOCATE (local_j_al(dimen_ri, my_b_size(jspin), block_size))
925 local_j_al = 0.0_dp
926
927 IF (calc_forces) THEN
928 ALLOCATE (y_i_ap(my_b_size(ispin), dimen_ri, block_size))
929 y_i_ap = 0.0_dp
930 ! For closed-shell, alpha-alpha and beta-beta my_B_size_beta=my_b_size
931 ! Not for alpha-beta case: Y_j_aP_beta is sent and received as Y_j_aP
932 ALLOCATE (y_j_ap(my_b_size(jspin), dimen_ri, block_size))
933 y_j_ap = 0.0_dp
934 END IF
935 !
936
937 CALL timestop(handle)
938
939 END SUBROUTINE mp2_ri_allocate_blk
940
941! **************************************************************************************************
942!> \brief ...
943!> \param my_alpha_beta_case ...
944!> \param total_ij_pairs ...
945!> \param homo ...
946!> \param homo_beta ...
947!> \param block_size ...
948!> \param ngroup ...
949!> \param ij_map ...
950!> \param color_sub ...
951!> \param my_ij_pairs ...
952!> \param my_open_shell_SS ...
953!> \param unit_nr ...
954! **************************************************************************************************
955 SUBROUTINE mp2_ri_communication(my_alpha_beta_case, total_ij_pairs, homo, homo_beta, &
956 block_size, ngroup, ij_map, color_sub, my_ij_pairs, my_open_shell_SS, unit_nr)
957 LOGICAL, INTENT(IN) :: my_alpha_beta_case
958 INTEGER, INTENT(OUT) :: total_ij_pairs
959 INTEGER, INTENT(IN) :: homo, homo_beta, block_size, ngroup
960 INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(OUT) :: ij_map
961 INTEGER, INTENT(IN) :: color_sub
962 INTEGER, INTENT(OUT) :: my_ij_pairs
963 LOGICAL, INTENT(IN) :: my_open_shell_ss
964 INTEGER, INTENT(IN) :: unit_nr
965
966 CHARACTER(LEN=*), PARAMETER :: routinen = 'mp2_ri_communication'
967
968 INTEGER :: assigned_blocks, first_i_block, first_j_block, handle, iib, ij_block_counter, &
969 ij_counter, jjb, last_i_block, last_j_block, num_block_per_group, num_ij_blocks, &
970 num_ij_blocks_beta, total_ij_block, total_ij_pairs_blocks
971 LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: ij_marker
972
973! Calculate the maximum number of ij pairs that have to be computed
974! among groups
975
976 CALL timeset(routinen, handle)
977
978 IF (.NOT. my_open_shell_ss .AND. .NOT. my_alpha_beta_case) THEN
979 total_ij_pairs = homo*(1 + homo)/2
980 num_ij_blocks = homo/block_size - 1
981
982 first_i_block = 1
983 last_i_block = block_size*(num_ij_blocks - 1)
984
985 first_j_block = block_size + 1
986 last_j_block = block_size*(num_ij_blocks + 1)
987
988 ij_block_counter = 0
989 DO iib = first_i_block, last_i_block, block_size
990 DO jjb = iib + block_size, last_j_block, block_size
991 ij_block_counter = ij_block_counter + 1
992 END DO
993 END DO
994
995 total_ij_block = ij_block_counter
996 num_block_per_group = total_ij_block/ngroup
997 assigned_blocks = num_block_per_group*ngroup
998
999 total_ij_pairs_blocks = assigned_blocks + (total_ij_pairs - assigned_blocks*(block_size**2))
1000
1001 ALLOCATE (ij_marker(homo, homo))
1002 ij_marker = .true.
1003 ALLOCATE (ij_map(3, total_ij_pairs_blocks))
1004 ij_map = 0
1005 ij_counter = 0
1006 my_ij_pairs = 0
1007 DO iib = first_i_block, last_i_block, block_size
1008 DO jjb = iib + block_size, last_j_block, block_size
1009 IF (ij_counter + 1 > assigned_blocks) EXIT
1010 ij_counter = ij_counter + 1
1011 ij_marker(iib:iib + block_size - 1, jjb:jjb + block_size - 1) = .false.
1012 ij_map(1, ij_counter) = iib
1013 ij_map(2, ij_counter) = jjb
1014 ij_map(3, ij_counter) = block_size
1015 IF (mod(ij_counter, ngroup) == color_sub) my_ij_pairs = my_ij_pairs + 1
1016 END DO
1017 END DO
1018 DO iib = 1, homo
1019 DO jjb = iib, homo
1020 IF (ij_marker(iib, jjb)) THEN
1021 ij_counter = ij_counter + 1
1022 ij_map(1, ij_counter) = iib
1023 ij_map(2, ij_counter) = jjb
1024 ij_map(3, ij_counter) = 1
1025 IF (mod(ij_counter, ngroup) == color_sub) my_ij_pairs = my_ij_pairs + 1
1026 END IF
1027 END DO
1028 END DO
1029 DEALLOCATE (ij_marker)
1030
1031 ELSE IF (.NOT. my_alpha_beta_case) THEN
1032 ! THese are the cases alpha/alpha and beta/beta
1033 ! We do not have to consider the diagonal elements
1034 total_ij_pairs = homo*(homo - 1)/2
1035 num_ij_blocks = (homo - 1)/block_size - 1
1036
1037 first_i_block = 1
1038 last_i_block = block_size*(num_ij_blocks - 1)
1039
1040 ! We shift the blocks to prevent the calculation of the diagonal elements which always give zero
1041 first_j_block = block_size + 2
1042 last_j_block = block_size*(num_ij_blocks + 1) + 1
1043
1044 ij_block_counter = 0
1045 DO iib = first_i_block, last_i_block, block_size
1046 DO jjb = iib + block_size + 1, last_j_block, block_size
1047 ij_block_counter = ij_block_counter + 1
1048 END DO
1049 END DO
1050
1051 total_ij_block = ij_block_counter
1052 num_block_per_group = total_ij_block/ngroup
1053 assigned_blocks = num_block_per_group*ngroup
1054
1055 total_ij_pairs_blocks = assigned_blocks + (total_ij_pairs - assigned_blocks*(block_size**2))
1056
1057 ALLOCATE (ij_marker(homo, homo))
1058 ij_marker = .true.
1059 ALLOCATE (ij_map(3, total_ij_pairs_blocks))
1060 ij_map = 0
1061 ij_counter = 0
1062 my_ij_pairs = 0
1063 DO iib = first_i_block, last_i_block, block_size
1064 DO jjb = iib + block_size + 1, last_j_block, block_size
1065 IF (ij_counter + 1 > assigned_blocks) EXIT
1066 ij_counter = ij_counter + 1
1067 ij_marker(iib:iib + block_size - 1, jjb:jjb + block_size - 1) = .false.
1068 ij_map(1, ij_counter) = iib
1069 ij_map(2, ij_counter) = jjb
1070 ij_map(3, ij_counter) = block_size
1071 IF (mod(ij_counter, ngroup) == color_sub) my_ij_pairs = my_ij_pairs + 1
1072 END DO
1073 END DO
1074 DO iib = 1, homo
1075 DO jjb = iib + 1, homo
1076 IF (ij_marker(iib, jjb)) THEN
1077 ij_counter = ij_counter + 1
1078 ij_map(1, ij_counter) = iib
1079 ij_map(2, ij_counter) = jjb
1080 ij_map(3, ij_counter) = 1
1081 IF (mod(ij_counter, ngroup) == color_sub) my_ij_pairs = my_ij_pairs + 1
1082 END IF
1083 END DO
1084 END DO
1085 DEALLOCATE (ij_marker)
1086
1087 ELSE
1088 total_ij_pairs = homo*homo_beta
1089 num_ij_blocks = homo/block_size
1090 num_ij_blocks_beta = homo_beta/block_size
1091
1092 first_i_block = 1
1093 last_i_block = block_size*(num_ij_blocks - 1)
1094
1095 first_j_block = 1
1096 last_j_block = block_size*(num_ij_blocks_beta - 1)
1097
1098 ij_block_counter = 0
1099 DO iib = first_i_block, last_i_block, block_size
1100 DO jjb = first_j_block, last_j_block, block_size
1101 ij_block_counter = ij_block_counter + 1
1102 END DO
1103 END DO
1104
1105 total_ij_block = ij_block_counter
1106 num_block_per_group = total_ij_block/ngroup
1107 assigned_blocks = num_block_per_group*ngroup
1108
1109 total_ij_pairs_blocks = assigned_blocks + (total_ij_pairs - assigned_blocks*(block_size**2))
1110
1111 ALLOCATE (ij_marker(homo, homo_beta))
1112 ij_marker = .true.
1113 ALLOCATE (ij_map(3, total_ij_pairs_blocks))
1114 ij_map = 0
1115 ij_counter = 0
1116 my_ij_pairs = 0
1117 DO iib = first_i_block, last_i_block, block_size
1118 DO jjb = first_j_block, last_j_block, block_size
1119 IF (ij_counter + 1 > assigned_blocks) EXIT
1120 ij_counter = ij_counter + 1
1121 ij_marker(iib:iib + block_size - 1, jjb:jjb + block_size - 1) = .false.
1122 ij_map(1, ij_counter) = iib
1123 ij_map(2, ij_counter) = jjb
1124 ij_map(3, ij_counter) = block_size
1125 IF (mod(ij_counter, ngroup) == color_sub) my_ij_pairs = my_ij_pairs + 1
1126 END DO
1127 END DO
1128 DO iib = 1, homo
1129 DO jjb = 1, homo_beta
1130 IF (ij_marker(iib, jjb)) THEN
1131 ij_counter = ij_counter + 1
1132 ij_map(1, ij_counter) = iib
1133 ij_map(2, ij_counter) = jjb
1134 ij_map(3, ij_counter) = 1
1135 IF (mod(ij_counter, ngroup) == color_sub) my_ij_pairs = my_ij_pairs + 1
1136 END IF
1137 END DO
1138 END DO
1139 DEALLOCATE (ij_marker)
1140 END IF
1141
1142 IF (unit_nr > 0) THEN
1143 IF (block_size == 1) THEN
1144 WRITE (unit=unit_nr, fmt="(T3,A,T66,F15.1)") &
1145 "RI_INFO| Percentage of ij pairs communicated with block size 1:", 100.0_dp
1146 ELSE
1147 WRITE (unit=unit_nr, fmt="(T3,A,T66,F15.1)") &
1148 "RI_INFO| Percentage of ij pairs communicated with block size 1:", &
1149 100.0_dp*real((total_ij_pairs - assigned_blocks*(block_size**2)), kind=dp)/real(total_ij_pairs, kind=dp)
1150 END IF
1151 CALL m_flush(unit_nr)
1152 END IF
1153
1154 CALL timestop(handle)
1155
1156 END SUBROUTINE mp2_ri_communication
1157
1158! **************************************************************************************************
1159!> \brief ...
1160!> \param para_env ...
1161!> \param para_env_sub ...
1162!> \param color_sub ...
1163!> \param sizes_array ...
1164!> \param calc_forces ...
1165!> \param integ_group_size ...
1166!> \param my_group_L_end ...
1167!> \param my_group_L_size ...
1168!> \param my_group_L_size_orig ...
1169!> \param my_group_L_start ...
1170!> \param my_new_group_L_size ...
1171!> \param integ_group_pos2color_sub ...
1172!> \param sizes_array_orig ...
1173!> \param ranges_info_array ...
1174!> \param comm_exchange ...
1175!> \param comm_rep ...
1176!> \param num_integ_group ...
1177! **************************************************************************************************
1178 SUBROUTINE mp2_ri_create_group(para_env, para_env_sub, color_sub, &
1179 sizes_array, calc_forces, &
1180 integ_group_size, my_group_L_end, &
1181 my_group_L_size, my_group_L_size_orig, my_group_L_start, my_new_group_L_size, &
1182 integ_group_pos2color_sub, &
1183 sizes_array_orig, ranges_info_array, comm_exchange, comm_rep, num_integ_group)
1184 TYPE(mp_para_env_type), INTENT(IN) :: para_env, para_env_sub
1185 INTEGER, INTENT(IN) :: color_sub
1186 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(INOUT) :: sizes_array
1187 LOGICAL, INTENT(IN) :: calc_forces
1188 INTEGER, INTENT(IN) :: integ_group_size, my_group_l_end
1189 INTEGER, INTENT(INOUT) :: my_group_l_size
1190 INTEGER, INTENT(OUT) :: my_group_l_size_orig
1191 INTEGER, INTENT(IN) :: my_group_l_start
1192 INTEGER, INTENT(INOUT) :: my_new_group_l_size
1193 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: integ_group_pos2color_sub, &
1194 sizes_array_orig
1195 INTEGER, ALLOCATABLE, DIMENSION(:, :, :), &
1196 INTENT(OUT) :: ranges_info_array
1197 TYPE(mp_comm_type), INTENT(OUT) :: comm_exchange, comm_rep
1198 INTEGER, INTENT(IN) :: num_integ_group
1199
1200 CHARACTER(LEN=*), PARAMETER :: routinen = 'mp2_ri_create_group'
1201
1202 INTEGER :: handle, iib, proc_receive, proc_shift, &
1203 sub_sub_color
1204 INTEGER, ALLOCATABLE, DIMENSION(:) :: new_sizes_array, rep_ends_array, &
1205 rep_sizes_array, rep_starts_array
1206 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: my_info
1207
1208 CALL timeset(routinen, handle)
1209 !
1210 sub_sub_color = para_env_sub%mepos*num_integ_group + color_sub/integ_group_size
1211 CALL comm_exchange%from_split(para_env, sub_sub_color)
1212
1213 ! create the replication group
1214 sub_sub_color = para_env_sub%mepos*comm_exchange%num_pe + comm_exchange%mepos
1215 CALL comm_rep%from_split(para_env, sub_sub_color)
1216
1217 ! create the new limits for K according to the size
1218 ! of the integral group
1219
1220 ! info array for replication
1221 ALLOCATE (rep_ends_array(0:comm_rep%num_pe - 1))
1222 ALLOCATE (rep_starts_array(0:comm_rep%num_pe - 1))
1223 ALLOCATE (rep_sizes_array(0:comm_rep%num_pe - 1))
1224
1225 CALL comm_rep%allgather(my_group_l_size, rep_sizes_array)
1226 CALL comm_rep%allgather(my_group_l_start, rep_starts_array)
1227 CALL comm_rep%allgather(my_group_l_end, rep_ends_array)
1228
1229 ! calculate my_new_group_L_size according to sizes_array
1230 my_new_group_l_size = my_group_l_size
1231
1232 ! Info of this process
1233 ALLOCATE (my_info(4, 0:comm_rep%num_pe - 1))
1234 my_info(1, 0) = my_group_l_start
1235 my_info(2, 0) = my_group_l_end
1236 my_info(3, 0) = 1
1237 my_info(4, 0) = my_group_l_size
1238
1239 DO proc_shift = 1, comm_rep%num_pe - 1
1240 proc_receive = modulo(comm_rep%mepos - proc_shift, comm_rep%num_pe)
1241
1242 my_new_group_l_size = my_new_group_l_size + rep_sizes_array(proc_receive)
1243
1244 my_info(1, proc_shift) = rep_starts_array(proc_receive)
1245 my_info(2, proc_shift) = rep_ends_array(proc_receive)
1246 my_info(3, proc_shift) = my_info(4, proc_shift - 1) + 1
1247 my_info(4, proc_shift) = my_new_group_l_size
1248
1249 END DO
1250
1251 ALLOCATE (new_sizes_array(0:comm_exchange%num_pe - 1))
1252 ALLOCATE (ranges_info_array(4, 0:comm_rep%num_pe - 1, 0:comm_exchange%num_pe - 1))
1253 CALL comm_exchange%allgather(my_new_group_l_size, new_sizes_array)
1254 CALL comm_exchange%allgather(my_info, ranges_info_array)
1255
1256 DEALLOCATE (rep_sizes_array)
1257 DEALLOCATE (rep_starts_array)
1258 DEALLOCATE (rep_ends_array)
1259
1260 ALLOCATE (integ_group_pos2color_sub(0:comm_exchange%num_pe - 1))
1261 CALL comm_exchange%allgather(color_sub, integ_group_pos2color_sub)
1262
1263 IF (calc_forces) THEN
1264 iib = SIZE(sizes_array)
1265 ALLOCATE (sizes_array_orig(0:iib - 1))
1266 sizes_array_orig(:) = sizes_array
1267 END IF
1268
1269 my_group_l_size_orig = my_group_l_size
1270 my_group_l_size = my_new_group_l_size
1271 DEALLOCATE (sizes_array)
1272
1273 ALLOCATE (sizes_array(0:integ_group_size - 1))
1274 sizes_array(:) = new_sizes_array
1275
1276 DEALLOCATE (new_sizes_array)
1277 !
1278 CALL timestop(handle)
1279
1280 END SUBROUTINE mp2_ri_create_group
1281
1282! **************************************************************************************************
1283!> \brief ...
1284!> \param mp2_env ...
1285!> \param para_env ...
1286!> \param para_env_sub ...
1287!> \param gd_array ...
1288!> \param gd_B_virtual ...
1289!> \param homo ...
1290!> \param dimen_RI ...
1291!> \param unit_nr ...
1292!> \param integ_group_size ...
1293!> \param ngroup ...
1294!> \param num_integ_group ...
1295!> \param virtual ...
1296!> \param calc_forces ...
1297! **************************************************************************************************
1298 SUBROUTINE mp2_ri_get_integ_group_size(mp2_env, para_env, para_env_sub, gd_array, gd_B_virtual, &
1299 homo, dimen_RI, unit_nr, &
1300 integ_group_size, &
1301 ngroup, num_integ_group, &
1302 virtual, calc_forces)
1303 TYPE(mp2_type) :: mp2_env
1304 TYPE(mp_para_env_type), INTENT(IN) :: para_env, para_env_sub
1305 TYPE(group_dist_d1_type), INTENT(IN) :: gd_array
1306 TYPE(group_dist_d1_type), DIMENSION(:), INTENT(IN) :: gd_b_virtual
1307 INTEGER, DIMENSION(:), INTENT(IN) :: homo
1308 INTEGER, INTENT(IN) :: dimen_ri, unit_nr
1309 INTEGER, INTENT(OUT) :: integ_group_size, ngroup, num_integ_group
1310 INTEGER, DIMENSION(:), INTENT(IN) :: virtual
1311 LOGICAL, INTENT(IN) :: calc_forces
1312
1313 CHARACTER(LEN=*), PARAMETER :: routinen = 'mp2_ri_get_integ_group_size'
1314
1315 INTEGER :: block_size, handle, iib, &
1316 max_repl_group_size, &
1317 min_integ_group_size
1318 INTEGER(KIND=int_8) :: mem
1319 LOGICAL :: calc_group_size
1320 REAL(kind=dp) :: factor, mem_base, mem_min, mem_per_blk, &
1321 mem_per_repl, mem_per_repl_blk, &
1322 mem_real
1323
1324 CALL timeset(routinen, handle)
1325
1326 ngroup = para_env%num_pe/para_env_sub%num_pe
1327
1328 calc_group_size = mp2_env%ri_mp2%number_integration_groups <= 0
1329 IF (.NOT. calc_group_size) THEN
1330 IF (mod(ngroup, mp2_env%ri_mp2%number_integration_groups) /= 0) calc_group_size = .true.
1331 END IF
1332
1333 IF (calc_group_size) THEN
1334 CALL m_memory(mem)
1335 mem_real = (mem + 1024*1024 - 1)/(1024*1024)
1336 CALL para_env%min(mem_real)
1337 mem_real = min(mem_real, mp2_env%mp2_memory)
1338
1339 mem_base = 0.0_dp
1340 mem_per_blk = 0.0_dp
1341 mem_per_repl = 0.0_dp
1342 mem_per_repl_blk = 0.0_dp
1343
1344 ! BIB_C_copy
1345 mem_per_repl = mem_per_repl + maxval(max(real(homo, kind=dp)*maxsize(gd_array), real(dimen_ri, kind=dp))* &
1346 maxsize(gd_b_virtual))*8.0_dp/(1024**2)
1347 ! BIB_C
1348 mem_per_repl = mem_per_repl + sum(real(homo, kind=dp)*maxsize(gd_b_virtual))*maxsize(gd_array)*8.0_dp/(1024**2)
1349 ! BIB_C_rec
1350 mem_per_repl_blk = mem_per_repl_blk + real(maxval(maxsize(gd_b_virtual)), kind=dp)*maxsize(gd_array)*8.0_dp/(1024**2)
1351 ! local_i_aL+local_j_aL
1352 mem_per_blk = mem_per_blk + 2.0_dp*maxval(maxsize(gd_b_virtual))*real(dimen_ri, kind=dp)*8.0_dp/(1024**2)
1353 ! local_ab
1354 mem_base = mem_base + maxval(real(virtual, kind=dp)*maxsize(gd_b_virtual))*8.0_dp/(1024**2)
1355 ! external_ab/external_i_aL
1356 mem_base = mem_base + real(max(dimen_ri, maxval(virtual)), kind=dp)*maxval(maxsize(gd_b_virtual))*8.0_dp/(1024**2)
1357
1358 IF (calc_forces) THEN
1359 ! Gamma_P_ia
1360 mem_per_repl = mem_per_repl + sum(real(homo, kind=dp)*maxsize(gd_array)* &
1361 maxsize(gd_b_virtual))*8.0_dp/(1024**2)
1362 ! Y_i_aP+Y_j_aP
1363 mem_per_blk = mem_per_blk + 2.0_dp*maxval(maxsize(gd_b_virtual))*dimen_ri*8.0_dp/(1024**2)
1364 ! local_ba/t_ab
1365 mem_base = mem_base + real(maxval(maxsize(gd_b_virtual)), kind=dp)*max(dimen_ri, maxval(virtual))*8.0_dp/(1024**2)
1366 ! P_ij
1367 mem_base = mem_base + sum(real(homo, kind=dp)*homo)*8.0_dp/(1024**2)
1368 ! P_ab
1369 mem_base = mem_base + sum(real(virtual, kind=dp)*maxsize(gd_b_virtual))*8.0_dp/(1024**2)
1370 ! send_ab/send_i_aL
1371 mem_base = mem_base + real(max(dimen_ri, maxval(virtual)), kind=dp)*maxval(maxsize(gd_b_virtual))*8.0_dp/(1024**2)
1372 END IF
1373
1374 ! This a first guess based on the assumption of optimal block sizes
1375 block_size = max(1, min(floor(sqrt(real(minval(homo), kind=dp))), floor(minval(homo)/sqrt(2.0_dp*ngroup))))
1376 IF (mp2_env%ri_mp2%block_size > 0) block_size = mp2_env%ri_mp2%block_size
1377
1378 mem_min = mem_base + mem_per_repl + (mem_per_blk + mem_per_repl_blk)*block_size
1379
1380 IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T68,F9.2,A4)') 'RI_INFO| Minimum available memory per MPI process:', &
1381 mem_real, ' MiB'
1382 IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T68,F9.2,A4)') 'RI_INFO| Minimum required memory per MPI process:', &
1383 mem_min, ' MiB'
1384
1385 ! We use the following communication model
1386 ! Comm(replication)+Comm(collection of data for ij pair)+Comm(contraction)
1387 ! One can show that the costs of the contraction step are independent of the block size and the replication group size
1388 ! With gradients, the other two steps are carried out twice (Y_i_aP -> Gamma_i_aP, and dereplication)
1389 ! NL ... number of RI basis functions
1390 ! NR ... replication group size
1391 ! NG ... number of sub groups
1392 ! NB ... Block size
1393 ! o ... number of occupied orbitals
1394 ! Then, we have the communication costs (in multiples of the original BIb_C matrix)
1395 ! (NR/NG)+(1-(NR/NG))*(o/NB+NB-2)/NG = (NR/NG)*(1-(o/NB+NB-2)/NG)+(o/NB+NB-2)/NG
1396 ! and with gradients
1397 ! 2*(NR/NG)+2*(1-(NR/NG))*(o/NB+NB-2)/NG = (NR/NG)*(1-(o/NB+NB-2)/NG)+(o/NB+NB-2)/NG
1398 ! We are looking for the minimum of the communication volume,
1399 ! thus, if the prefactor of (NR/NG) is smaller than zero, use the largest possible replication group size.
1400 ! If the factor is larger than zero, set the replication group size to 1.
1401 ! (For small systems and a large number of subgroups)
1402 ! Replication group size = 1 implies that the integration group size equals the number of subgroups
1403
1404 integ_group_size = ngroup
1405
1406 ! Multiply everything by homo*virtual to consider differences between spin channels in case of open-shell calculations
1407 factor = real(sum(homo*virtual), kind=dp) &
1408 - sum((real(maxval(homo), kind=dp)/block_size + block_size - 2.0_dp)*homo*virtual)/ngroup
1409 IF (SIZE(homo) == 2) factor = factor - 2.0_dp*product(homo)/block_size/ngroup*sum(homo*virtual)
1410
1411 IF (factor <= 0.0_dp) THEN
1412 ! Remove the fixed memory and divide by the memory per replication group size
1413 max_repl_group_size = min(max(floor((mem_real - mem_base - mem_per_blk*block_size)/ &
1414 (mem_per_repl + mem_per_repl_blk*block_size)), 1), ngroup)
1415 ! Convert to an integration group size
1416 min_integ_group_size = ngroup/max_repl_group_size
1417
1418 ! Ensure that the integration group size is a divisor of the number of sub groups
1419 DO iib = max(min(min_integ_group_size, ngroup), 1), ngroup
1420 ! check that the ngroup is a multiple of integ_group_size
1421 IF (mod(ngroup, iib) == 0) THEN
1422 integ_group_size = iib
1423 EXIT
1424 END IF
1425 integ_group_size = integ_group_size + 1
1426 END DO
1427 END IF
1428 ELSE ! We take the user provided group size
1429 integ_group_size = ngroup/mp2_env%ri_mp2%number_integration_groups
1430 END IF
1431
1432 IF (unit_nr > 0) THEN
1433 WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
1434 "RI_INFO| Group size for integral replication:", integ_group_size*para_env_sub%num_pe
1435 CALL m_flush(unit_nr)
1436 END IF
1437
1438 num_integ_group = ngroup/integ_group_size
1439
1440 CALL timestop(handle)
1441
1442 END SUBROUTINE mp2_ri_get_integ_group_size
1443
1444! **************************************************************************************************
1445!> \brief ...
1446!> \param mp2_env ...
1447!> \param para_env ...
1448!> \param para_env_sub ...
1449!> \param gd_array ...
1450!> \param gd_B_virtual ...
1451!> \param homo ...
1452!> \param virtual ...
1453!> \param dimen_RI ...
1454!> \param unit_nr ...
1455!> \param block_size ...
1456!> \param ngroup ...
1457!> \param num_integ_group ...
1458!> \param my_open_shell_ss ...
1459!> \param calc_forces ...
1460!> \param buffer_1D ...
1461! **************************************************************************************************
1462 SUBROUTINE mp2_ri_get_block_size(mp2_env, para_env, para_env_sub, gd_array, gd_B_virtual, &
1463 homo, virtual, dimen_RI, unit_nr, &
1464 block_size, ngroup, num_integ_group, &
1465 my_open_shell_ss, calc_forces, buffer_1D)
1466 TYPE(mp2_type) :: mp2_env
1467 TYPE(mp_para_env_type), INTENT(IN) :: para_env, para_env_sub
1468 TYPE(group_dist_d1_type), INTENT(IN) :: gd_array
1469 TYPE(group_dist_d1_type), DIMENSION(:), INTENT(IN) :: gd_b_virtual
1470 INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
1471 INTEGER, INTENT(IN) :: dimen_ri, unit_nr
1472 INTEGER, INTENT(OUT) :: block_size, ngroup
1473 INTEGER, INTENT(IN) :: num_integ_group
1474 LOGICAL, INTENT(IN) :: my_open_shell_ss, calc_forces
1475 REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
1476 INTENT(OUT) :: buffer_1d
1477
1478 CHARACTER(LEN=*), PARAMETER :: routinen = 'mp2_ri_get_block_size'
1479
1480 INTEGER :: best_block_size, handle, num_ij_blocks
1481 INTEGER(KIND=int_8) :: buffer_size, mem
1482 REAL(kind=dp) :: mem_base, mem_per_blk, mem_per_repl_blk, &
1483 mem_real
1484
1485 CALL timeset(routinen, handle)
1486
1487 ngroup = para_env%num_pe/para_env_sub%num_pe
1488
1489 CALL m_memory(mem)
1490 mem_real = (mem + 1024*1024 - 1)/(1024*1024)
1491 CALL para_env%min(mem_real)
1492
1493 mem_base = 0.0_dp
1494 mem_per_blk = 0.0_dp
1495 mem_per_repl_blk = 0.0_dp
1496
1497 ! external_ab
1498 mem_base = mem_base + maxval(maxsize(gd_b_virtual))*max(dimen_ri, maxval(virtual))*8.0_dp/(1024**2)
1499 ! BIB_C_rec
1500 mem_per_repl_blk = mem_per_repl_blk + real(maxval(maxsize(gd_b_virtual)), kind=dp)*maxsize(gd_array)*8.0_dp/(1024**2)
1501 ! local_i_aL+local_j_aL
1502 mem_per_blk = mem_per_blk + 2.0_dp*maxval(maxsize(gd_b_virtual))*real(dimen_ri, kind=dp)*8.0_dp/(1024**2)
1503 ! Copy to keep arrays contiguous
1504 mem_base = mem_base + maxval(maxsize(gd_b_virtual))*max(dimen_ri, maxval(virtual))*8.0_dp/(1024**2)
1505
1506 IF (calc_forces) THEN
1507 ! Y_i_aP+Y_j_aP+BIb_C_send
1508 mem_per_blk = mem_per_blk + 3.0_dp*maxval(maxsize(gd_b_virtual))*dimen_ri*8.0_dp/(1024**2)
1509 ! send_ab
1510 mem_base = mem_base + maxval(maxsize(gd_b_virtual))*max(dimen_ri, maxval(virtual))*8.0_dp/(1024**2)
1511 END IF
1512
1513 best_block_size = 1
1514
1515 ! Here we split the memory half for the communication, half for replication
1516 IF (mp2_env%ri_mp2%block_size > 0) THEN
1517 best_block_size = mp2_env%ri_mp2%block_size
1518 ELSE
1519 best_block_size = max(floor((mem_real - mem_base)/(mem_per_blk + mem_per_repl_blk*ngroup/num_integ_group)), 1)
1520
1521 DO
1522 IF (SIZE(homo) == 1) THEN
1523 IF (.NOT. my_open_shell_ss) THEN
1524 num_ij_blocks = (homo(1)/best_block_size)
1525 num_ij_blocks = (num_ij_blocks*num_ij_blocks - num_ij_blocks)/2
1526 ELSE
1527 num_ij_blocks = ((homo(1) - 1)/best_block_size)
1528 num_ij_blocks = (num_ij_blocks*num_ij_blocks - num_ij_blocks)/2
1529 END IF
1530 ELSE
1531 num_ij_blocks = product(homo/best_block_size)
1532 END IF
1533 ! Enforce at least one large block for each subgroup
1534 IF ((num_ij_blocks >= ngroup .AND. num_ij_blocks > 0) .OR. best_block_size == 1) THEN
1535 EXIT
1536 ELSE
1537 best_block_size = best_block_size - 1
1538 END IF
1539 END DO
1540
1541 IF (SIZE(homo) == 1) THEN
1542 IF (my_open_shell_ss) THEN
1543 ! check that best_block_size is not bigger than sqrt(homo-1)
1544 ! Diagonal elements do not have to be considered
1545 best_block_size = min(floor(sqrt(real(homo(1) - 1, kind=dp))), best_block_size)
1546 ELSE
1547 ! check that best_block_size is not bigger than sqrt(homo)
1548 best_block_size = min(floor(sqrt(real(homo(1), kind=dp))), best_block_size)
1549 END IF
1550 END IF
1551 END IF
1552 block_size = max(1, best_block_size)
1553
1554 IF (unit_nr > 0) THEN
1555 WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
1556 "RI_INFO| Block size:", block_size
1557 CALL m_flush(unit_nr)
1558 END IF
1559
1560 ! Determine recv buffer size (BI_C_recv, external_i_aL, external_ab)
1561 buffer_size = max(int(maxsize(gd_array), kind=int_8)*block_size, int(max(dimen_ri, maxval(virtual)), kind=int_8)) &
1562 *maxval(maxsize(gd_b_virtual))
1563 ! The send buffer has the same size as the recv buffer
1564 IF (calc_forces) buffer_size = buffer_size*2
1565 ALLOCATE (buffer_1d(buffer_size))
1566
1567 CALL timestop(handle)
1568
1569 END SUBROUTINE mp2_ri_get_block_size
1570
1571! **************************************************************************************************
1572!> \brief ...
1573!> \param mp2_env ...
1574!> \param para_env_sub ...
1575!> \param gd_B_virtual ...
1576!> \param Eigenval ...
1577!> \param homo ...
1578!> \param dimen_RI ...
1579!> \param iiB ...
1580!> \param jjB ...
1581!> \param my_B_size ...
1582!> \param my_B_virtual_end ...
1583!> \param my_B_virtual_start ...
1584!> \param my_i ...
1585!> \param my_j ...
1586!> \param virtual ...
1587!> \param local_ab ...
1588!> \param t_ab ...
1589!> \param my_local_i_aL ...
1590!> \param my_local_j_aL ...
1591!> \param open_ss ...
1592!> \param Y_i_aP ...
1593!> \param Y_j_aP ...
1594!> \param local_ba ...
1595!> \param ispin ...
1596!> \param jspin ...
1597!> \param dgemm_counter ...
1598!> \param buffer_1D ...
1599! **************************************************************************************************
1600 SUBROUTINE mp2_update_p_gamma(mp2_env, para_env_sub, gd_B_virtual, &
1601 Eigenval, homo, dimen_RI, iiB, jjB, my_B_size, &
1602 my_B_virtual_end, my_B_virtual_start, my_i, my_j, virtual, local_ab, &
1603 t_ab, my_local_i_aL, my_local_j_aL, open_ss, Y_i_aP, Y_j_aP, &
1604 local_ba, ispin, jspin, dgemm_counter, buffer_1D)
1605 TYPE(mp2_type) :: mp2_env
1606 TYPE(mp_para_env_type), INTENT(IN) :: para_env_sub
1607 TYPE(group_dist_d1_type), DIMENSION(:), INTENT(IN) :: gd_b_virtual
1608 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: eigenval
1609 INTEGER, DIMENSION(:), INTENT(IN) :: homo
1610 INTEGER, INTENT(IN) :: dimen_ri, iib, jjb
1611 INTEGER, DIMENSION(:), INTENT(IN) :: my_b_size, my_b_virtual_end, &
1612 my_b_virtual_start
1613 INTEGER, INTENT(IN) :: my_i, my_j
1614 INTEGER, DIMENSION(:), INTENT(IN) :: virtual
1615 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
1616 INTENT(INOUT), TARGET :: local_ab
1617 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
1618 INTENT(IN), TARGET :: t_ab, my_local_i_al, my_local_j_al
1619 LOGICAL, INTENT(IN) :: open_ss
1620 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
1621 INTENT(INOUT), TARGET :: y_i_ap, y_j_ap, local_ba
1622 INTEGER, INTENT(IN) :: ispin, jspin
1623 TYPE(dgemm_counter_type), INTENT(INOUT) :: dgemm_counter
1624 REAL(kind=dp), CONTIGUOUS, DIMENSION(:), TARGET :: buffer_1d
1625
1626 CHARACTER(LEN=*), PARAMETER :: routinen = 'mp2_update_P_gamma'
1627
1628 INTEGER :: a, b, b_global, handle, proc_receive, proc_send, proc_shift, rec_b_size, &
1629 rec_b_virtual_end, rec_b_virtual_start, send_b_size, send_b_virtual_end, &
1630 send_b_virtual_start
1631 INTEGER(KIND=int_8) :: offset
1632 LOGICAL :: alpha_beta
1633 REAL(kind=dp) :: factor, p_ij_diag
1634 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
1635 POINTER :: external_ab, send_ab
1636
1637 CALL timeset(routinen//"_Pia", handle)
1638
1639 alpha_beta = .NOT. (ispin == jspin)
1640 IF (open_ss) THEN
1641 factor = 1.0_dp
1642 ELSE
1643 factor = 2.0_dp
1644 END IF
1645 ! divide the (ia|jb) integrals by Delta_ij^ab
1646 DO b = 1, my_b_size(jspin)
1647 b_global = b + my_b_virtual_start(jspin) - 1
1648 DO a = 1, virtual(ispin)
1649 local_ab(a, b) = -local_ab(a, b)/ &
1650 (eigenval(homo(ispin) + a, ispin) + eigenval(homo(jspin) + b_global, jspin) - &
1651 eigenval(my_i + iib - 1, ispin) - eigenval(my_j + jjb - 1, jspin))
1652 END DO
1653 END DO
1654 IF (.NOT. (alpha_beta)) THEN
1655 p_ij_diag = -sum(local_ab*t_ab)*factor
1656 ELSE
1657 ! update diagonal part of P_ij
1658 p_ij_diag = -sum(local_ab*local_ab)*mp2_env%scale_S
1659 ! More integrals needed only for alpha-beta case: local_ba
1660 DO b = 1, my_b_size(ispin)
1661 b_global = b + my_b_virtual_start(ispin) - 1
1662 DO a = 1, virtual(jspin)
1663 local_ba(a, b) = -local_ba(a, b)/ &
1664 (eigenval(homo(jspin) + a, jspin) + eigenval(homo(ispin) + b_global, ispin) - &
1665 eigenval(my_i + iib - 1, ispin) - eigenval(my_j + jjb - 1, jspin))
1666 END DO
1667 END DO
1668 END IF
1669
1670 ! P_ab and add diagonal part of P_ij
1671
1672 CALL dgemm_counter_start(dgemm_counter)
1673 IF (.NOT. (alpha_beta)) THEN
1674 CALL mp2_env%local_gemm_ctx%gemm('T', 'N', my_b_size(ispin), my_b_size(ispin), virtual(ispin), 1.0_dp, &
1675 t_ab, virtual(ispin), local_ab, virtual(ispin), &
1676 1.0_dp, mp2_env%ri_grad%P_ab(ispin)%array(:, &
1677 my_b_virtual_start(ispin):my_b_virtual_end(ispin)), my_b_size(ispin))
1678 mp2_env%ri_grad%P_ij(ispin)%array(my_i + iib - 1, my_i + iib - 1) = &
1679 mp2_env%ri_grad%P_ij(ispin)%array(my_i + iib - 1, my_i + iib - 1) + p_ij_diag
1680 ELSE
1681 CALL mp2_env%local_gemm_ctx%gemm('T', 'N', my_b_size(ispin), my_b_size(ispin), virtual(jspin), mp2_env%scale_S, &
1682 local_ba, virtual(jspin), local_ba, virtual(jspin), 1.0_dp, &
1683 mp2_env%ri_grad%P_ab(ispin)%array(:, my_b_virtual_start(ispin):my_b_virtual_end(ispin)), my_b_size(ispin))
1684
1685 mp2_env%ri_grad%P_ij(ispin)%array(my_i + iib - 1, my_i + iib - 1) = &
1686 mp2_env%ri_grad%P_ij(ispin)%array(my_i + iib - 1, my_i + iib - 1) + p_ij_diag
1687
1688 CALL mp2_env%local_gemm_ctx%gemm('T', 'N', my_b_size(jspin), my_b_size(jspin), virtual(ispin), mp2_env%scale_S, &
1689 local_ab, virtual(ispin), local_ab, virtual(ispin), 1.0_dp, &
1690 mp2_env%ri_grad%P_ab(jspin)%array(:, my_b_virtual_start(jspin):my_b_virtual_end(jspin)), my_b_size(jspin))
1691
1692 mp2_env%ri_grad%P_ij(jspin)%array(my_j + jjb - 1, my_j + jjb - 1) = &
1693 mp2_env%ri_grad%P_ij(jspin)%array(my_j + jjb - 1, my_j + jjb - 1) + p_ij_diag
1694 END IF
1695 ! The summation is over unique pairs. In alpha-beta case, all pairs are unique: subroutine is called for
1696 ! both i^alpha,j^beta and i^beta,j^alpha. Formally, my_i can be equal to my_j, but they are different
1697 ! due to spin in alpha-beta case.
1698 IF ((my_i /= my_j) .AND. (.NOT. alpha_beta)) THEN
1699
1700 CALL mp2_env%local_gemm_ctx%gemm('N', 'T', my_b_size(ispin), virtual(ispin), my_b_size(ispin), 1.0_dp, &
1701 t_ab(my_b_virtual_start(ispin):my_b_virtual_end(ispin), :), my_b_size(ispin), &
1702 local_ab, virtual(ispin), &
1703 1.0_dp, mp2_env%ri_grad%P_ab(ispin)%array, my_b_size(ispin))
1704
1705 mp2_env%ri_grad%P_ij(ispin)%array(my_j + jjb - 1, my_j + jjb - 1) = &
1706 mp2_env%ri_grad%P_ij(ispin)%array(my_j + jjb - 1, my_j + jjb - 1) + p_ij_diag
1707 END IF
1708 DO proc_shift = 1, para_env_sub%num_pe - 1
1709 proc_send = modulo(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
1710 proc_receive = modulo(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
1711
1712 CALL get_group_dist(gd_b_virtual(jspin), proc_receive, rec_b_virtual_start, rec_b_virtual_end, rec_b_size)
1713 CALL get_group_dist(gd_b_virtual(jspin), proc_send, send_b_virtual_start, send_b_virtual_end, send_b_size)
1714
1715 external_ab(1:virtual(ispin), 1:rec_b_size) => buffer_1d(1:int(virtual(ispin), int_8)*rec_b_size)
1716 external_ab = 0.0_dp
1717
1718 CALL para_env_sub%sendrecv(local_ab, proc_send, &
1719 external_ab, proc_receive)
1720
1721 IF (.NOT. (alpha_beta)) THEN
1722 CALL mp2_env%local_gemm_ctx%gemm('T', 'N', my_b_size(ispin), rec_b_size, virtual(ispin), 1.0_dp, &
1723 t_ab, virtual(ispin), external_ab, virtual(ispin), &
1724 1.0_dp, mp2_env%ri_grad%P_ab(ispin)%array(:, rec_b_virtual_start:rec_b_virtual_end), my_b_size(ispin))
1725 ELSE
1726 CALL mp2_env%local_gemm_ctx%gemm('T', 'N', my_b_size(jspin), rec_b_size, virtual(ispin), mp2_env%scale_S, &
1727 local_ab, virtual(ispin), external_ab, virtual(ispin), &
1728 1.0_dp, mp2_env%ri_grad%P_ab(jspin)%array(:, rec_b_virtual_start:rec_b_virtual_end), &
1729 my_b_size(jspin))
1730
1731 ! For alpha-beta part of alpha-density we need a new parallel code
1732 ! And new external_ab (of a different size)
1733 CALL get_group_dist(gd_b_virtual(ispin), proc_receive, rec_b_virtual_start, rec_b_virtual_end, rec_b_size)
1734 CALL get_group_dist(gd_b_virtual(ispin), proc_send, send_b_virtual_start, send_b_virtual_end, send_b_size)
1735 external_ab(1:virtual(jspin), 1:rec_b_size) => buffer_1d(1:int(virtual(jspin), int_8)*rec_b_size)
1736 external_ab = 0.0_dp
1737 CALL para_env_sub%sendrecv(local_ba, proc_send, &
1738 external_ab, proc_receive)
1739 CALL mp2_env%local_gemm_ctx%gemm('T', 'N', my_b_size(ispin), rec_b_size, virtual(jspin), mp2_env%scale_S, &
1740 local_ba, virtual(jspin), external_ab, virtual(jspin), &
1741 1.0_dp, mp2_env%ri_grad%P_ab(ispin)%array(:, rec_b_virtual_start:rec_b_virtual_end), my_b_size(ispin))
1742 END IF
1743
1744 IF ((my_i /= my_j) .AND. (.NOT. alpha_beta)) THEN
1745 external_ab(1:my_b_size(ispin), 1:virtual(ispin)) => &
1746 buffer_1d(1:int(virtual(ispin), int_8)*my_b_size(ispin))
1747 external_ab = 0.0_dp
1748
1749 offset = int(virtual(ispin), int_8)*my_b_size(ispin)
1750
1751 send_ab(1:send_b_size, 1:virtual(ispin)) => buffer_1d(offset + 1:offset + int(send_b_size, int_8)*virtual(ispin))
1752 send_ab = 0.0_dp
1753
1754 CALL mp2_env%local_gemm_ctx%gemm('N', 'T', send_b_size, virtual(ispin), my_b_size(ispin), 1.0_dp, &
1755 t_ab(send_b_virtual_start:send_b_virtual_end, :), send_b_size, &
1756 local_ab, virtual(ispin), 0.0_dp, send_ab, send_b_size)
1757 CALL para_env_sub%sendrecv(send_ab, proc_send, &
1758 external_ab, proc_receive)
1759
1760 mp2_env%ri_grad%P_ab(ispin)%array(:, :) = mp2_env%ri_grad%P_ab(ispin)%array + external_ab
1761 END IF
1762
1763 END DO
1764 IF (.NOT. alpha_beta) THEN
1765 IF (my_i /= my_j) THEN
1766 CALL dgemm_counter_stop(dgemm_counter, 2*my_b_size(ispin), virtual(ispin), virtual(ispin))
1767 ELSE
1768 CALL dgemm_counter_stop(dgemm_counter, my_b_size(ispin), virtual(ispin), virtual(ispin))
1769 END IF
1770 ELSE
1771 CALL dgemm_counter_stop(dgemm_counter, sum(my_b_size), virtual(ispin), virtual(jspin))
1772 END IF
1773 CALL timestop(handle)
1774
1775 ! Now, Gamma_P_ia (made of Y_ia_P)
1776
1777 CALL timeset(routinen//"_Gamma", handle)
1778 CALL dgemm_counter_start(dgemm_counter)
1779 IF (.NOT. alpha_beta) THEN
1780 ! Alpha-alpha, beta-beta and closed shell
1781 CALL mp2_env%local_gemm_ctx%gemm('N', 'T', my_b_size(ispin), dimen_ri, my_b_size(ispin), 1.0_dp, &
1782 t_ab(my_b_virtual_start(ispin):my_b_virtual_end(ispin), :), my_b_size(ispin), &
1783 my_local_j_al, dimen_ri, 1.0_dp, y_i_ap, my_b_size(ispin))
1784 ELSE ! Alpha-beta
1785 CALL mp2_env%local_gemm_ctx%gemm('N', 'T', my_b_size(ispin), dimen_ri, my_b_size(jspin), mp2_env%scale_S, &
1786 local_ab(my_b_virtual_start(ispin):my_b_virtual_end(ispin), :), my_b_size(ispin), &
1787 my_local_j_al, dimen_ri, 1.0_dp, y_i_ap, my_b_size(ispin))
1788 CALL mp2_env%local_gemm_ctx%gemm('T', 'T', my_b_size(jspin), dimen_ri, my_b_size(ispin), mp2_env%scale_S, &
1789 local_ab(my_b_virtual_start(ispin):my_b_virtual_end(ispin), :), my_b_size(ispin), &
1790 my_local_i_al, dimen_ri, 1.0_dp, y_j_ap, my_b_size(jspin))
1791 END IF
1792
1793 IF (para_env_sub%num_pe > 1) THEN
1794 external_ab(1:my_b_size(ispin), 1:dimen_ri) => buffer_1d(1:int(my_b_size(ispin), int_8)*dimen_ri)
1795 external_ab = 0.0_dp
1796
1797 offset = int(my_b_size(ispin), int_8)*dimen_ri
1798 END IF
1799 !
1800 DO proc_shift = 1, para_env_sub%num_pe - 1
1801 proc_send = modulo(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
1802 proc_receive = modulo(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
1803
1804 CALL get_group_dist(gd_b_virtual(ispin), proc_receive, rec_b_virtual_start, rec_b_virtual_end, rec_b_size)
1805 CALL get_group_dist(gd_b_virtual(ispin), proc_send, send_b_virtual_start, send_b_virtual_end, send_b_size)
1806
1807 send_ab(1:send_b_size, 1:dimen_ri) => buffer_1d(offset + 1:offset + int(dimen_ri, int_8)*send_b_size)
1808 send_ab = 0.0_dp
1809 IF (.NOT. alpha_beta) THEN
1810 CALL mp2_env%local_gemm_ctx%gemm('N', 'T', send_b_size, dimen_ri, my_b_size(ispin), 1.0_dp, &
1811 t_ab(send_b_virtual_start:send_b_virtual_end, :), send_b_size, &
1812 my_local_j_al, dimen_ri, 0.0_dp, send_ab, send_b_size)
1813 CALL para_env_sub%sendrecv(send_ab, proc_send, external_ab, proc_receive)
1814
1815 y_i_ap(:, :) = y_i_ap + external_ab
1816
1817 ELSE ! Alpha-beta case
1818 ! Alpha-alpha part
1819 CALL mp2_env%local_gemm_ctx%gemm('N', 'T', send_b_size, dimen_ri, my_b_size(jspin), mp2_env%scale_S, &
1820 local_ab(send_b_virtual_start:send_b_virtual_end, :), send_b_size, &
1821 my_local_j_al, dimen_ri, 0.0_dp, send_ab, send_b_size)
1822 CALL para_env_sub%sendrecv(send_ab, proc_send, external_ab, proc_receive)
1823 y_i_ap(:, :) = y_i_ap + external_ab
1824 END IF
1825 END DO
1826
1827 IF (alpha_beta) THEN
1828 ! For beta-beta part (in alpha-beta case) we need a new parallel code
1829 IF (para_env_sub%num_pe > 1) THEN
1830 external_ab(1:my_b_size(jspin), 1:dimen_ri) => buffer_1d(1:int(my_b_size(jspin), int_8)*dimen_ri)
1831 external_ab = 0.0_dp
1832
1833 offset = int(my_b_size(jspin), int_8)*dimen_ri
1834 END IF
1835 DO proc_shift = 1, para_env_sub%num_pe - 1
1836 proc_send = modulo(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
1837 proc_receive = modulo(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
1838
1839 CALL get_group_dist(gd_b_virtual(jspin), proc_send, send_b_virtual_start, send_b_virtual_end, send_b_size)
1840 send_ab(1:send_b_size, 1:dimen_ri) => buffer_1d(offset + 1:offset + int(dimen_ri, int_8)*send_b_size)
1841 send_ab = 0.0_dp
1842 CALL mp2_env%local_gemm_ctx%gemm('N', 'T', send_b_size, dimen_ri, my_b_size(ispin), mp2_env%scale_S, &
1843 local_ba(send_b_virtual_start:send_b_virtual_end, :), send_b_size, &
1844 my_local_i_al, dimen_ri, 0.0_dp, send_ab, send_b_size)
1845 CALL para_env_sub%sendrecv(send_ab, proc_send, external_ab, proc_receive)
1846 y_j_ap(:, :) = y_j_ap + external_ab
1847
1848 END DO
1849
1850 ! Here, we just use approximate bounds. For large systems virtual(ispin) is approx virtual(jspin), same for B_size
1851 CALL dgemm_counter_stop(dgemm_counter, 3*virtual(ispin), dimen_ri, my_b_size(jspin))
1852 ELSE
1853 CALL dgemm_counter_stop(dgemm_counter, virtual(ispin), dimen_ri, my_b_size(ispin))
1854 END IF
1855
1856 IF ((my_i /= my_j) .AND. (.NOT. alpha_beta)) THEN
1857 ! Alpha-alpha, beta-beta and closed shell
1858 CALL dgemm_counter_start(dgemm_counter)
1859 CALL mp2_env%local_gemm_ctx%gemm('T', 'T', my_b_size(ispin), dimen_ri, my_b_size(ispin), 1.0_dp, &
1860 t_ab(my_b_virtual_start(ispin):my_b_virtual_end(ispin), :), my_b_size(ispin), &
1861 my_local_i_al, dimen_ri, 1.0_dp, y_j_ap, my_b_size(ispin))
1862 DO proc_shift = 1, para_env_sub%num_pe - 1
1863 proc_send = modulo(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
1864 proc_receive = modulo(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
1865
1866 CALL get_group_dist(gd_b_virtual(ispin), proc_receive, rec_b_virtual_start, rec_b_virtual_end, rec_b_size)
1867
1868 external_ab(1:dimen_ri, 1:rec_b_size) => buffer_1d(1:int(dimen_ri, int_8)*rec_b_size)
1869 external_ab = 0.0_dp
1870
1871 CALL para_env_sub%sendrecv(my_local_i_al, proc_send, &
1872 external_ab, proc_receive)
1873
1874 ! Alpha-alpha, beta-beta and closed shell
1875 CALL mp2_env%local_gemm_ctx%gemm('T', 'T', my_b_size(ispin), dimen_ri, rec_b_size, 1.0_dp, &
1876 t_ab(rec_b_virtual_start:rec_b_virtual_end, :), rec_b_size, &
1877 external_ab, dimen_ri, 1.0_dp, y_j_ap, my_b_size(ispin))
1878 END DO
1879
1880 CALL dgemm_counter_stop(dgemm_counter, my_b_size(ispin), dimen_ri, virtual(ispin))
1881 END IF
1882
1883 CALL timestop(handle)
1884 END SUBROUTINE mp2_update_p_gamma
1885
1886! **************************************************************************************************
1887!> \brief ...
1888!> \param Gamma_P_ia ...
1889!> \param ij_index ...
1890!> \param my_B_size ...
1891!> \param my_block_size ...
1892!> \param my_group_L_size ...
1893!> \param my_i ...
1894!> \param my_ij_pairs ...
1895!> \param ngroup ...
1896!> \param num_integ_group ...
1897!> \param integ_group_pos2color_sub ...
1898!> \param num_ij_pairs ...
1899!> \param ij_map ...
1900!> \param ranges_info_array ...
1901!> \param Y_i_aP ...
1902!> \param comm_exchange ...
1903!> \param sizes_array ...
1904!> \param spin ...
1905!> \param buffer_1D ...
1906! **************************************************************************************************
1907 SUBROUTINE mp2_redistribute_gamma(Gamma_P_ia, ij_index, my_B_size, &
1908 my_block_size, my_group_L_size, my_i, my_ij_pairs, ngroup, &
1909 num_integ_group, integ_group_pos2color_sub, num_ij_pairs, &
1910 ij_map, ranges_info_array, Y_i_aP, comm_exchange, &
1911 sizes_array, spin, buffer_1D)
1912
1913 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: gamma_p_ia
1914 INTEGER, INTENT(IN) :: ij_index, my_b_size, my_block_size, &
1915 my_group_l_size, my_i, my_ij_pairs, &
1916 ngroup, num_integ_group
1917 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(IN) :: integ_group_pos2color_sub, num_ij_pairs
1918 INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(IN) :: ij_map
1919 INTEGER, ALLOCATABLE, DIMENSION(:, :, :), &
1920 INTENT(IN) :: ranges_info_array
1921 REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: y_i_ap
1922 TYPE(mp_comm_type), INTENT(IN) :: comm_exchange
1923 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(IN) :: sizes_array
1924 INTEGER, INTENT(IN) :: spin
1925 REAL(kind=dp), CONTIGUOUS, DIMENSION(:), TARGET :: buffer_1d
1926
1927 CHARACTER(LEN=*), PARAMETER :: routinen = 'mp2_redistribute_gamma'
1928
1929 INTEGER :: end_point, handle, handle2, iib, ij_counter_rec, irep, kkk, lll, lstart_pos, &
1930 proc_receive, proc_send, proc_shift, rec_i, rec_ij_index, send_l_size, start_point, tag
1931 INTEGER(KIND=int_8) :: offset
1932 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
1933 POINTER :: bi_c_rec, bi_c_send
1934
1935! In alpha-beta case Y_i_aP_beta is sent as Y_j_aP
1936
1937 CALL timeset(routinen//"_comm2", handle)
1938
1939 tag = 43
1940
1941 IF (ij_index <= my_ij_pairs) THEN
1942 ! somethig to send
1943 ! start with myself
1944 CALL timeset(routinen//"_comm2_w", handle2)
1945 DO irep = 0, num_integ_group - 1
1946 lstart_pos = ranges_info_array(1, irep, comm_exchange%mepos)
1947 start_point = ranges_info_array(3, irep, comm_exchange%mepos)
1948 end_point = ranges_info_array(4, irep, comm_exchange%mepos)
1949!$OMP PARALLEL DO DEFAULT(NONE) &
1950!$OMP PRIVATE(kkk,lll,iiB) &
1951!$OMP SHARED(start_point,end_point,Lstart_pos,my_block_size,&
1952!$OMP Gamma_P_ia,my_i,my_B_size,Y_i_aP)
1953 DO kkk = start_point, end_point
1954 lll = kkk - start_point + lstart_pos
1955 DO iib = 1, my_block_size
1956 gamma_p_ia(1:my_b_size, my_i + iib - 1, kkk) = &
1957 gamma_p_ia(1:my_b_size, my_i + iib - 1, kkk) + &
1958 y_i_ap(1:my_b_size, lll, iib)
1959 END DO
1960 END DO
1961!$OMP END PARALLEL DO
1962 END DO
1963 CALL timestop(handle2)
1964
1965 ! Y_i_aP(my_B_size,dimen_RI,block_size)
1966
1967 DO proc_shift = 1, comm_exchange%num_pe - 1
1968 proc_send = modulo(comm_exchange%mepos + proc_shift, comm_exchange%num_pe)
1969 proc_receive = modulo(comm_exchange%mepos - proc_shift, comm_exchange%num_pe)
1970
1971 send_l_size = sizes_array(proc_send)
1972 bi_c_send(1:my_b_size, 1:my_block_size, 1:send_l_size) => &
1973 buffer_1d(1:int(my_b_size, int_8)*my_block_size*send_l_size)
1974
1975 offset = int(my_b_size, int_8)*my_block_size*send_l_size
1976
1977 CALL timeset(routinen//"_comm2_w", handle2)
1978 bi_c_send = 0.0_dp
1979 DO irep = 0, num_integ_group - 1
1980 lstart_pos = ranges_info_array(1, irep, proc_send)
1981 start_point = ranges_info_array(3, irep, proc_send)
1982 end_point = ranges_info_array(4, irep, proc_send)
1983!$OMP PARALLEL DO DEFAULT(NONE) &
1984!$OMP PRIVATE(kkk,lll,iiB) &
1985!$OMP SHARED(start_point,end_point,Lstart_pos,my_block_size,&
1986!$OMP BI_C_send,my_B_size,Y_i_aP)
1987 DO kkk = start_point, end_point
1988 lll = kkk - start_point + lstart_pos
1989 DO iib = 1, my_block_size
1990 bi_c_send(1:my_b_size, iib, kkk) = y_i_ap(1:my_b_size, lll, iib)
1991 END DO
1992 END DO
1993!$OMP END PARALLEL DO
1994 END DO
1995 CALL timestop(handle2)
1996
1997 rec_ij_index = num_ij_pairs(proc_receive)
1998
1999 IF (ij_index <= rec_ij_index) THEN
2000 ! we know that proc_receive has something to send for us, let's see what
2001 ij_counter_rec = &
2002 (ij_index - min(1, integ_group_pos2color_sub(proc_receive)))*ngroup + integ_group_pos2color_sub(proc_receive)
2003
2004 rec_i = ij_map(spin, ij_counter_rec)
2005
2006 bi_c_rec(1:my_b_size, 1:my_block_size, 1:my_group_l_size) => &
2007 buffer_1d(offset + 1:offset + int(my_b_size, int_8)*my_block_size*my_group_l_size)
2008 bi_c_rec = 0.0_dp
2009
2010 CALL comm_exchange%sendrecv(bi_c_send, proc_send, &
2011 bi_c_rec, proc_receive, tag)
2012
2013 CALL timeset(routinen//"_comm2_w", handle2)
2014 DO irep = 0, num_integ_group - 1
2015 start_point = ranges_info_array(3, irep, comm_exchange%mepos)
2016 end_point = ranges_info_array(4, irep, comm_exchange%mepos)
2017!$OMP PARALLEL WORKSHARE DEFAULT(NONE) &
2018!$OMP SHARED(start_point,end_point,my_block_size,&
2019!$OMP Gamma_P_ia,rec_i,iiB,my_B_size,BI_C_rec)
2020 gamma_p_ia(:, rec_i:rec_i + my_block_size - 1, start_point:end_point) = &
2021 gamma_p_ia(:, rec_i:rec_i + my_block_size - 1, start_point:end_point) + &
2022 bi_c_rec(1:my_b_size, :, start_point:end_point)
2023!$OMP END PARALLEL WORKSHARE
2024 END DO
2025 CALL timestop(handle2)
2026
2027 ELSE
2028 ! we have something to send but nothing to receive
2029 CALL comm_exchange%send(bi_c_send, proc_send, tag)
2030
2031 END IF
2032
2033 END DO
2034
2035 ELSE
2036 ! noting to send check if we have to receive
2037 DO proc_shift = 1, comm_exchange%num_pe - 1
2038 proc_send = modulo(comm_exchange%mepos + proc_shift, comm_exchange%num_pe)
2039 proc_receive = modulo(comm_exchange%mepos - proc_shift, comm_exchange%num_pe)
2040 rec_ij_index = num_ij_pairs(proc_receive)
2041
2042 IF (ij_index <= rec_ij_index) THEN
2043 ! we know that proc_receive has something to send for us, let's see what
2044 ij_counter_rec = &
2045 (ij_index - min(1, integ_group_pos2color_sub(proc_receive)))*ngroup + integ_group_pos2color_sub(proc_receive)
2046
2047 rec_i = ij_map(spin, ij_counter_rec)
2048
2049 bi_c_rec(1:my_b_size, 1:my_block_size, 1:my_group_l_size) => &
2050 buffer_1d(1:int(my_b_size, int_8)*my_block_size*my_group_l_size)
2051
2052 bi_c_rec = 0.0_dp
2053
2054 CALL comm_exchange%recv(bi_c_rec, proc_receive, tag)
2055
2056 CALL timeset(routinen//"_comm2_w", handle2)
2057 DO irep = 0, num_integ_group - 1
2058 start_point = ranges_info_array(3, irep, comm_exchange%mepos)
2059 end_point = ranges_info_array(4, irep, comm_exchange%mepos)
2060#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
2061!$OMP PARALLEL WORKSHARE DEFAULT(NONE) &
2062!$OMP SHARED(start_point,end_point,my_block_size,&
2063!$OMP Gamma_P_ia,rec_i,my_B_size,BI_C_rec)
2064#endif
2065 gamma_p_ia(:, rec_i:rec_i + my_block_size - 1, start_point:end_point) = &
2066 gamma_p_ia(:, rec_i:rec_i + my_block_size - 1, start_point:end_point) + &
2067 bi_c_rec(1:my_b_size, :, start_point:end_point)
2068#if !defined(__INTEL_LLVM_COMPILER) || (20250000 <= __INTEL_LLVM_COMPILER)
2069!$OMP END PARALLEL WORKSHARE
2070#endif
2071 END DO
2072 CALL timestop(handle2)
2073
2074 END IF
2075 END DO
2076
2077 END IF
2078 CALL timestop(handle)
2079
2080 END SUBROUTINE mp2_redistribute_gamma
2081
2082! **************************************************************************************************
2083!> \brief ...
2084!> \param mp2_env ...
2085!> \param Eigenval ...
2086!> \param homo ...
2087!> \param virtual ...
2088!> \param open_shell ...
2089!> \param beta_beta ...
2090!> \param Bib_C ...
2091!> \param unit_nr ...
2092!> \param dimen_RI ...
2093!> \param my_B_size ...
2094!> \param ngroup ...
2095!> \param my_group_L_size ...
2096!> \param color_sub ...
2097!> \param ranges_info_array ...
2098!> \param comm_exchange ...
2099!> \param para_env_sub ...
2100!> \param para_env ...
2101!> \param my_B_virtual_start ...
2102!> \param my_B_virtual_end ...
2103!> \param sizes_array ...
2104!> \param gd_B_virtual ...
2105!> \param integ_group_pos2color_sub ...
2106!> \param dgemm_counter ...
2107!> \param buffer_1D ...
2108! **************************************************************************************************
2109 SUBROUTINE quasi_degenerate_p_ij(mp2_env, Eigenval, homo, virtual, open_shell, &
2110 beta_beta, Bib_C, unit_nr, dimen_RI, &
2111 my_B_size, ngroup, my_group_L_size, &
2112 color_sub, ranges_info_array, comm_exchange, para_env_sub, para_env, &
2113 my_B_virtual_start, my_B_virtual_end, sizes_array, gd_B_virtual, &
2114 integ_group_pos2color_sub, dgemm_counter, buffer_1D)
2115 TYPE(mp2_type) :: mp2_env
2116 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: eigenval
2117 INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
2118 LOGICAL, INTENT(IN) :: open_shell, beta_beta
2119 TYPE(three_dim_real_array), DIMENSION(:), &
2120 INTENT(IN) :: bib_c
2121 INTEGER, INTENT(IN) :: unit_nr, dimen_ri
2122 INTEGER, DIMENSION(:), INTENT(IN) :: my_b_size
2123 INTEGER, INTENT(IN) :: ngroup, my_group_l_size, color_sub
2124 INTEGER, ALLOCATABLE, DIMENSION(:, :, :), &
2125 INTENT(IN) :: ranges_info_array
2126 TYPE(mp_comm_type), INTENT(IN) :: comm_exchange
2127 TYPE(mp_para_env_type), INTENT(IN) :: para_env_sub, para_env
2128 INTEGER, DIMENSION(:), INTENT(IN) :: my_b_virtual_start, my_b_virtual_end
2129 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(IN) :: sizes_array
2130 TYPE(group_dist_d1_type), DIMENSION(:), INTENT(IN) :: gd_b_virtual
2131 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(IN) :: integ_group_pos2color_sub
2132 TYPE(dgemm_counter_type), INTENT(INOUT) :: dgemm_counter
2133 REAL(kind=dp), CONTIGUOUS, DIMENSION(:), TARGET :: buffer_1d
2134
2135 CHARACTER(LEN=*), PARAMETER :: routinen = 'quasi_degenerate_P_ij'
2136
2137 INTEGER :: a, a_global, b, b_global, block_size, decil, handle, handle2, ijk_counter, &
2138 ijk_counter_send, ijk_index, ispin, kkb, kspin, max_block_size, max_ijk, my_i, my_ijk, &
2139 my_j, my_k, my_last_k(2), my_virtual, nspins, proc_receive, proc_send, proc_shift, &
2140 rec_b_size, rec_b_virtual_end, rec_b_virtual_start, rec_l_size, send_b_size, &
2141 send_b_virtual_end, send_b_virtual_start, send_i, send_ijk_index, send_j, send_k, &
2142 size_b_i, size_b_k, tag, tag2
2143 INTEGER, ALLOCATABLE, DIMENSION(:) :: num_ijk
2144 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: ijk_map, send_last_k
2145 LOGICAL :: alpha_beta, do_recv_i, do_recv_j, &
2146 do_recv_k, do_send_i, do_send_j, &
2147 do_send_k
2148 REAL(kind=dp) :: amp_fac, p_ij_elem, t_new, t_start
2149 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
2150 TARGET :: local_ab, local_al_i, local_al_j, t_ab
2151 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: local_al_k
2152 REAL(kind=dp), DIMENSION(:, :), POINTER :: bi_c_rec, external_ab, external_al
2153 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: bi_c_rec_3d
2154
2155 CALL timeset(routinen//"_ij_sing", handle)
2156
2157 tag = 44
2158 tag2 = 45
2159
2160 nspins = SIZE(bib_c)
2161 alpha_beta = (nspins == 2)
2162
2163 ! Set amplitude factor
2164 amp_fac = mp2_env%scale_S + mp2_env%scale_T
2165 IF (open_shell) amp_fac = mp2_env%scale_T
2166
2167 ALLOCATE (send_last_k(2, comm_exchange%num_pe - 1))
2168
2169 ! Loop(s) over orbital triplets
2170 DO ispin = 1, nspins
2171 size_b_i = my_b_size(ispin)
2172 IF (ispin == 1 .AND. alpha_beta) THEN
2173 kspin = 2
2174 ELSE
2175 kspin = 1
2176 END IF
2177 size_b_k = my_b_size(kspin)
2178
2179 ! Find the number of quasi-degenerate orbitals and orbital triplets
2180
2181 CALL find_quasi_degenerate_ij(my_ijk, homo(ispin), homo(kspin), eigenval(:, ispin), mp2_env, ijk_map, unit_nr, ngroup, &
2182 .NOT. beta_beta .AND. ispin /= 2, comm_exchange, num_ijk, max_ijk, color_sub, &
2183 SIZE(buffer_1d), my_group_l_size, size_b_k, para_env, virtual(ispin), size_b_i)
2184
2185 my_virtual = virtual(ispin)
2186 IF (SIZE(ijk_map, 2) > 0) THEN
2187 max_block_size = ijk_map(4, 1)
2188 ELSE
2189 max_block_size = 1
2190 END IF
2191
2192 ALLOCATE (local_al_i(dimen_ri, size_b_i))
2193 ALLOCATE (local_al_j(dimen_ri, size_b_i))
2194 ALLOCATE (local_al_k(dimen_ri, size_b_k, max_block_size))
2195 ALLOCATE (t_ab(my_virtual, size_b_k))
2196
2197 my_last_k = -1
2198 send_last_k = -1
2199
2200 t_start = m_walltime()
2201 DO ijk_index = 1, max_ijk
2202
2203 ! Prediction is unreliable if we are in the first step of the loop
2204 IF (unit_nr > 0 .AND. ijk_index > 1) THEN
2205 decil = ijk_index*10/max_ijk
2206 IF (decil /= (ijk_index - 1)*10/max_ijk) THEN
2207 t_new = m_walltime()
2208 t_new = (t_new - t_start)/60.0_dp*(max_ijk - ijk_index + 1)/(ijk_index - 1)
2209 WRITE (unit_nr, fmt="(T3,A)") "Percentage of finished loop: "// &
2210 cp_to_string(decil*10)//". Minutes left: "//cp_to_string(t_new)
2211 CALL m_flush(unit_nr)
2212 END IF
2213 END IF
2214
2215 IF (ijk_index <= my_ijk) THEN
2216 ! work to be done
2217 ijk_counter = (ijk_index - min(1, color_sub))*ngroup + color_sub
2218 my_i = ijk_map(1, ijk_counter)
2219 my_j = ijk_map(2, ijk_counter)
2220 my_k = ijk_map(3, ijk_counter)
2221 block_size = ijk_map(4, ijk_counter)
2222
2223 do_recv_i = (ispin /= kspin) .OR. my_i < my_k .OR. my_i > my_k + block_size - 1
2224 do_recv_j = (ispin /= kspin) .OR. my_j < my_k .OR. my_j > my_k + block_size - 1
2225 do_recv_k = my_k /= my_last_k(1) .OR. my_k + block_size - 1 /= my_last_k(2)
2226 my_last_k(1) = my_k
2227 my_last_k(2) = my_k + block_size - 1
2228
2229 local_al_i = 0.0_dp
2230 IF (do_recv_i) THEN
2231 CALL fill_local_i_al_2d(local_al_i, ranges_info_array(:, :, comm_exchange%mepos), &
2232 bib_c(ispin)%array(:, :, my_i))
2233 END IF
2234
2235 local_al_j = 0.0_dp
2236 IF (do_recv_j) THEN
2237 CALL fill_local_i_al_2d(local_al_j, ranges_info_array(:, :, comm_exchange%mepos), &
2238 bib_c(ispin)%array(:, :, my_j))
2239 END IF
2240
2241 IF (do_recv_k) THEN
2242 local_al_k = 0.0_dp
2243 CALL fill_local_i_al(local_al_k(:, :, 1:block_size), ranges_info_array(:, :, comm_exchange%mepos), &
2244 bib_c(kspin)%array(:, :, my_k:my_k + block_size - 1))
2245 END IF
2246
2247 CALL timeset(routinen//"_comm", handle2)
2248 DO proc_shift = 1, comm_exchange%num_pe - 1
2249 proc_send = modulo(comm_exchange%mepos + proc_shift, comm_exchange%num_pe)
2250 proc_receive = modulo(comm_exchange%mepos - proc_shift, comm_exchange%num_pe)
2251
2252 send_ijk_index = num_ijk(proc_send)
2253
2254 rec_l_size = sizes_array(proc_receive)
2255 bi_c_rec(1:rec_l_size, 1:size_b_i) => buffer_1d(1:int(rec_l_size, kind=int_8)*size_b_i)
2256
2257 do_send_i = .false.
2258 do_send_j = .false.
2259 do_send_k = .false.
2260 IF (ijk_index <= send_ijk_index) THEN
2261 ! something to send
2262 ijk_counter_send = (ijk_index - min(1, integ_group_pos2color_sub(proc_send)))* &
2263 ngroup + integ_group_pos2color_sub(proc_send)
2264 send_i = ijk_map(1, ijk_counter_send)
2265 send_j = ijk_map(2, ijk_counter_send)
2266 send_k = ijk_map(3, ijk_counter_send)
2267
2268 do_send_i = (ispin /= kspin) .OR. send_i < send_k .OR. send_i > send_k + block_size - 1
2269 do_send_j = (ispin /= kspin) .OR. send_j < send_k .OR. send_j > send_k + block_size - 1
2270 do_send_k = send_k /= send_last_k(1, proc_shift) .OR. send_k + block_size - 1 /= send_last_k(2, proc_shift)
2271 send_last_k(1, proc_shift) = send_k
2272 send_last_k(2, proc_shift) = send_k + block_size - 1
2273 END IF
2274
2275 ! occupied i
2276 bi_c_rec = 0.0_dp
2277 IF (do_send_i) THEN
2278 IF (do_recv_i) THEN
2279 CALL comm_exchange%sendrecv(bib_c(ispin)%array(:, :, send_i), proc_send, &
2280 bi_c_rec, proc_receive, tag)
2281 ELSE
2282 CALL comm_exchange%send(bib_c(ispin)%array(:, :, send_i), proc_send, tag)
2283 END IF
2284 ELSE IF (do_recv_i) THEN
2285 CALL comm_exchange%recv(bi_c_rec, proc_receive, tag)
2286 END IF
2287 IF (do_recv_i) THEN
2288 CALL fill_local_i_al_2d(local_al_i, ranges_info_array(:, :, proc_receive), bi_c_rec)
2289 END IF
2290
2291 ! occupied j
2292 bi_c_rec = 0.0_dp
2293 IF (do_send_j) THEN
2294 IF (do_recv_j) THEN
2295 CALL comm_exchange%sendrecv(bib_c(ispin)%array(:, :, send_j), proc_send, &
2296 bi_c_rec, proc_receive, tag)
2297 ELSE
2298 CALL comm_exchange%send(bib_c(ispin)%array(:, :, send_j), proc_send, tag)
2299 END IF
2300 ELSE IF (do_recv_j) THEN
2301 CALL comm_exchange%recv(bi_c_rec, proc_receive, tag)
2302 END IF
2303 IF (do_recv_j) THEN
2304 CALL fill_local_i_al_2d(local_al_j, ranges_info_array(:, :, proc_receive), bi_c_rec)
2305 END IF
2306
2307 ! occupied k
2308 bi_c_rec_3d(1:rec_l_size, 1:size_b_k, 1:block_size) => &
2309 buffer_1d(1:int(rec_l_size, kind=int_8)*size_b_k*block_size)
2310 IF (do_send_k) THEN
2311 IF (do_recv_k) THEN
2312 CALL comm_exchange%sendrecv(bib_c(kspin)%array(:, :, send_k:send_k + block_size - 1), proc_send, &
2313 bi_c_rec_3d, proc_receive, tag)
2314 ELSE
2315 CALL comm_exchange%send(bi_c_rec, proc_receive, tag)
2316 END IF
2317 ELSE IF (do_recv_k) THEN
2318 CALL comm_exchange%recv(bi_c_rec_3d, proc_receive, tag)
2319 END IF
2320 IF (do_recv_k) THEN
2321 CALL fill_local_i_al(local_al_k(:, :, 1:block_size), ranges_info_array(:, :, proc_receive), bi_c_rec_3d)
2322 END IF
2323 END DO
2324
2325 IF (.NOT. do_recv_i) local_al_i(:, :) = local_al_k(:, :, my_i - my_k + 1)
2326 IF (.NOT. do_recv_j) local_al_j(:, :) = local_al_k(:, :, my_j - my_k + 1)
2327 CALL timestop(handle2)
2328
2329 ! expand integrals
2330 DO kkb = 1, block_size
2331 CALL timeset(routinen//"_exp_ik", handle2)
2332 CALL dgemm_counter_start(dgemm_counter)
2333 ALLOCATE (local_ab(my_virtual, size_b_k))
2334 local_ab = 0.0_dp
2335 CALL mp2_env%local_gemm_ctx%gemm('T', 'N', size_b_i, size_b_k, dimen_ri, 1.0_dp, &
2336 local_al_i, dimen_ri, local_al_k(:, :, kkb), dimen_ri, &
2337 0.0_dp, local_ab(my_b_virtual_start(ispin):my_b_virtual_end(ispin), 1:size_b_k), size_b_i)
2338 DO proc_shift = 1, para_env_sub%num_pe - 1
2339 proc_send = modulo(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
2340 proc_receive = modulo(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
2341
2342 CALL get_group_dist(gd_b_virtual(ispin), proc_receive, rec_b_virtual_start, rec_b_virtual_end, rec_b_size)
2343
2344 external_al(1:dimen_ri, 1:rec_b_size) => buffer_1d(1:int(dimen_ri, kind=int_8)*rec_b_size)
2345
2346 CALL comm_exchange%sendrecv(local_al_i, proc_send, &
2347 external_al, proc_receive, tag)
2348
2349 CALL mp2_env%local_gemm_ctx%gemm('T', 'N', rec_b_size, size_b_k, dimen_ri, 1.0_dp, &
2350 external_al, dimen_ri, local_al_k(:, :, kkb), dimen_ri, &
2351 0.0_dp, local_ab(rec_b_virtual_start:rec_b_virtual_end, 1:size_b_k), rec_b_size)
2352 END DO
2353 CALL dgemm_counter_stop(dgemm_counter, my_virtual, size_b_k, dimen_ri)
2354 CALL timestop(handle2)
2355
2356 ! Amplitudes
2357 CALL timeset(routinen//"_tab", handle2)
2358 t_ab = 0.0_dp
2359 ! Alpha-alpha, beta-beta and closed shell
2360 IF (.NOT. alpha_beta) THEN
2361 DO b = 1, size_b_k
2362 b_global = b + my_b_virtual_start(1) - 1
2363 DO a = 1, my_b_size(1)
2364 a_global = a + my_b_virtual_start(1) - 1
2365 t_ab(a_global, b) = (amp_fac*local_ab(a_global, b) - mp2_env%scale_T*local_ab(b_global, a))/ &
2366 (eigenval(my_i, 1) + eigenval(my_k + kkb - 1, 1) &
2367 - eigenval(homo(1) + a_global, 1) - eigenval(homo(1) + b_global, 1))
2368 END DO
2369 END DO
2370 ELSE
2371 DO b = 1, size_b_k
2372 b_global = b + my_b_virtual_start(kspin) - 1
2373 DO a = 1, my_b_size(ispin)
2374 a_global = a + my_b_virtual_start(ispin) - 1
2375 t_ab(a_global, b) = mp2_env%scale_S*local_ab(a_global, b)/ &
2376 (eigenval(my_i, ispin) + eigenval(my_k + kkb - 1, kspin) &
2377 - eigenval(homo(ispin) + a_global, ispin) - eigenval(homo(kspin) + b_global, kspin))
2378 END DO
2379 END DO
2380 END IF
2381
2382 IF (.NOT. alpha_beta) THEN
2383 DO proc_shift = 1, para_env_sub%num_pe - 1
2384 proc_send = modulo(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
2385 proc_receive = modulo(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
2386 CALL get_group_dist(gd_b_virtual(1), proc_receive, rec_b_virtual_start, rec_b_virtual_end, rec_b_size)
2387 CALL get_group_dist(gd_b_virtual(1), proc_send, send_b_virtual_start, send_b_virtual_end, send_b_size)
2388
2389 external_ab(1:size_b_i, 1:rec_b_size) => buffer_1d(1:int(size_b_i, kind=int_8)*rec_b_size)
2390 CALL para_env_sub%sendrecv(local_ab(send_b_virtual_start:send_b_virtual_end, 1:size_b_k), proc_send, &
2391 external_ab(1:size_b_i, 1:rec_b_size), proc_receive, tag)
2392
2393 DO b = 1, my_b_size(1)
2394 b_global = b + my_b_virtual_start(1) - 1
2395 DO a = 1, rec_b_size
2396 a_global = a + rec_b_virtual_start - 1
2397 t_ab(a_global, b) = (amp_fac*local_ab(a_global, b) - mp2_env%scale_T*external_ab(b, a))/ &
2398 (eigenval(my_i, 1) + eigenval(my_k + kkb - 1, 1) &
2399 - eigenval(homo(1) + a_global, 1) - eigenval(homo(1) + b_global, 1))
2400 END DO
2401 END DO
2402 END DO
2403 END IF
2404 CALL timestop(handle2)
2405
2406 ! Expand the second set of integrals
2407 CALL timeset(routinen//"_exp_jk", handle2)
2408 local_ab = 0.0_dp
2409 CALL dgemm_counter_start(dgemm_counter)
2410 CALL mp2_env%local_gemm_ctx%gemm('T', 'N', size_b_i, size_b_k, dimen_ri, 1.0_dp, &
2411 local_al_j, dimen_ri, local_al_k(:, :, kkb), dimen_ri, &
2412 0.0_dp, local_ab(my_b_virtual_start(ispin):my_b_virtual_end(ispin), 1:size_b_k), size_b_i)
2413 DO proc_shift = 1, para_env_sub%num_pe - 1
2414 proc_send = modulo(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
2415 proc_receive = modulo(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
2416
2417 CALL get_group_dist(gd_b_virtual(ispin), proc_receive, rec_b_virtual_start, rec_b_virtual_end, rec_b_size)
2418
2419 external_al(1:dimen_ri, 1:rec_b_size) => buffer_1d(1:int(dimen_ri, kind=int_8)*rec_b_size)
2420
2421 CALL comm_exchange%sendrecv(local_al_j, proc_send, &
2422 external_al, proc_receive, tag)
2423 CALL mp2_env%local_gemm_ctx%gemm('T', 'N', rec_b_size, size_b_k, dimen_ri, 1.0_dp, &
2424 external_al, dimen_ri, local_al_k(:, :, kkb), dimen_ri, &
2425 0.0_dp, local_ab(rec_b_virtual_start:rec_b_virtual_end, 1:size_b_k), rec_b_size)
2426 END DO
2427 CALL dgemm_counter_stop(dgemm_counter, my_virtual, size_b_k, dimen_ri)
2428 CALL timestop(handle2)
2429
2430 CALL timeset(routinen//"_Pij", handle2)
2431 DO b = 1, size_b_k
2432 b_global = b + my_b_virtual_start(kspin) - 1
2433 DO a = 1, my_b_size(ispin)
2434 a_global = a + my_b_virtual_start(ispin) - 1
2435 local_ab(a_global, b) = &
2436 local_ab(a_global, b)/(eigenval(my_j, ispin) + eigenval(my_k + kkb - 1, kspin) &
2437 - eigenval(homo(ispin) + a_global, ispin) - eigenval(homo(kspin) + b_global, kspin))
2438 END DO
2439 END DO
2440 !
2441 p_ij_elem = sum(local_ab*t_ab)
2442 DEALLOCATE (local_ab)
2443 IF ((.NOT. open_shell) .AND. (.NOT. alpha_beta)) THEN
2444 p_ij_elem = p_ij_elem*2.0_dp
2445 END IF
2446 IF (beta_beta) THEN
2447 mp2_env%ri_grad%P_ij(2)%array(my_i, my_j) = mp2_env%ri_grad%P_ij(2)%array(my_i, my_j) - p_ij_elem
2448 mp2_env%ri_grad%P_ij(2)%array(my_j, my_i) = mp2_env%ri_grad%P_ij(2)%array(my_j, my_i) - p_ij_elem
2449 ELSE
2450 mp2_env%ri_grad%P_ij(ispin)%array(my_i, my_j) = mp2_env%ri_grad%P_ij(ispin)%array(my_i, my_j) - p_ij_elem
2451 mp2_env%ri_grad%P_ij(ispin)%array(my_j, my_i) = mp2_env%ri_grad%P_ij(ispin)%array(my_j, my_i) - p_ij_elem
2452 END IF
2453 CALL timestop(handle2)
2454 END DO
2455 ELSE
2456 CALL timeset(routinen//"_comm", handle2)
2457 ! no work to be done, possible messeges to be exchanged
2458 DO proc_shift = 1, comm_exchange%num_pe - 1
2459 proc_send = modulo(comm_exchange%mepos + proc_shift, comm_exchange%num_pe)
2460 proc_receive = modulo(comm_exchange%mepos - proc_shift, comm_exchange%num_pe)
2461
2462 send_ijk_index = num_ijk(proc_send)
2463
2464 IF (ijk_index <= send_ijk_index) THEN
2465 ! somethig to send
2466 ijk_counter_send = (ijk_index - min(1, integ_group_pos2color_sub(proc_send)))*ngroup + &
2467 integ_group_pos2color_sub(proc_send)
2468 send_i = ijk_map(1, ijk_counter_send)
2469 send_j = ijk_map(2, ijk_counter_send)
2470 send_k = ijk_map(3, ijk_counter_send)
2471 block_size = ijk_map(4, ijk_counter_send)
2472
2473 do_send_i = (ispin /= kspin) .OR. send_i < send_k .OR. send_i > send_k + block_size - 1
2474 do_send_j = (ispin /= kspin) .OR. send_j < send_k .OR. send_j > send_k + block_size - 1
2475 ! occupied i
2476 IF (do_send_i) THEN
2477 CALL comm_exchange%send(bib_c(ispin)%array(:, :, send_i), proc_send, tag)
2478 END IF
2479 ! occupied j
2480 IF (do_send_j) THEN
2481 CALL comm_exchange%send(bib_c(ispin)%array(:, :, send_j), proc_send, tag)
2482 END IF
2483 ! occupied k
2484 CALL comm_exchange%send(bib_c(kspin)%array(:, :, send_k:send_k + block_size - 1), proc_send, tag)
2485 END IF
2486
2487 END DO ! proc loop
2488 CALL timestop(handle2)
2489 END IF
2490 END DO ! ijk_index loop
2491 DEALLOCATE (local_al_i)
2492 DEALLOCATE (local_al_j)
2493 DEALLOCATE (local_al_k)
2494 DEALLOCATE (t_ab)
2495 DEALLOCATE (ijk_map)
2496 END DO ! over number of loops (ispin)
2497 CALL timestop(handle)
2498
2499 END SUBROUTINE quasi_degenerate_p_ij
2500
2501! **************************************************************************************************
2502!> \brief ...
2503!> \param my_ijk ...
2504!> \param homo ...
2505!> \param homo_beta ...
2506!> \param Eigenval ...
2507!> \param mp2_env ...
2508!> \param ijk_map ...
2509!> \param unit_nr ...
2510!> \param ngroup ...
2511!> \param do_print_alpha ...
2512!> \param comm_exchange ...
2513!> \param num_ijk ...
2514!> \param max_ijk ...
2515!> \param color_sub ...
2516!> \param buffer_size ...
2517!> \param my_group_L_size ...
2518!> \param B_size_k ...
2519!> \param para_env ...
2520!> \param virtual ...
2521!> \param B_size_i ...
2522! **************************************************************************************************
2523 SUBROUTINE find_quasi_degenerate_ij(my_ijk, homo, homo_beta, Eigenval, mp2_env, ijk_map, unit_nr, ngroup, &
2524 do_print_alpha, comm_exchange, num_ijk, max_ijk, color_sub, &
2525 buffer_size, my_group_L_size, B_size_k, para_env, virtual, B_size_i)
2526
2527 INTEGER, INTENT(OUT) :: my_ijk
2528 INTEGER, INTENT(IN) :: homo, homo_beta
2529 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: eigenval
2530 TYPE(mp2_type), INTENT(IN) :: mp2_env
2531 INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(OUT) :: ijk_map
2532 INTEGER, INTENT(IN) :: unit_nr, ngroup
2533 LOGICAL, INTENT(IN) :: do_print_alpha
2534 TYPE(mp_comm_type), INTENT(IN) :: comm_exchange
2535 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: num_ijk
2536 INTEGER, INTENT(OUT) :: max_ijk
2537 INTEGER, INTENT(IN) :: color_sub, buffer_size, my_group_l_size, &
2538 b_size_k
2539 TYPE(mp_para_env_type), INTENT(IN) :: para_env
2540 INTEGER, INTENT(IN) :: virtual, b_size_i
2541
2542 INTEGER :: block_size, communication_steps, communication_volume, iib, ij_counter, &
2543 ijk_counter, jjb, kkb, max_block_size, max_num_k_blocks, min_communication_volume, &
2544 my_steps, num_k_blocks, num_sing_ij, total_ijk
2545 INTEGER(KIND=int_8) :: mem
2546 LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: ijk_marker
2547
2548 ALLOCATE (num_ijk(0:comm_exchange%num_pe - 1))
2549
2550 num_sing_ij = 0
2551 DO iib = 1, homo
2552 ! diagonal elements already updated
2553 DO jjb = iib + 1, homo
2554 IF (abs(eigenval(jjb) - eigenval(iib)) < mp2_env%ri_grad%eps_canonical) THEN
2555 num_sing_ij = num_sing_ij + 1
2556 END IF
2557 END DO
2558 END DO
2559
2560 IF (unit_nr > 0) THEN
2561 IF (do_print_alpha) THEN
2562 WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
2563 "MO_INFO| Number of ij pairs below EPS_CANONICAL:", num_sing_ij
2564 ELSE
2565 WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
2566 "MO_INFO| Number of ij pairs (spin beta) below EPS_CANONICAL:", num_sing_ij
2567 END IF
2568 END IF
2569
2570 ! Determine the block size, first guess: use available buffer
2571 max_block_size = buffer_size/(my_group_l_size*b_size_k)
2572
2573 ! Second limit: memory
2574 CALL m_memory(mem)
2575 ! Convert to number of doubles
2576 mem = mem/8
2577 ! Remove local_ab (2x) and local_aL_i (2x)
2578 mem = mem - 2_int_8*(virtual*b_size_k + b_size_i*my_group_l_size)
2579 max_block_size = min(max_block_size, max(1, int(mem/(my_group_l_size*b_size_k), kind(max_block_size))))
2580
2581 ! Exchange the limit
2582 CALL para_env%min(max_block_size)
2583
2584 ! Find now the block size which minimizes the communication volume and then the number of communication steps
2585 block_size = 1
2586 min_communication_volume = 3*homo_beta*num_sing_ij
2587 communication_steps = 3*homo_beta*num_sing_ij
2588 DO iib = max_block_size, 2, -1
2589 max_num_k_blocks = homo_beta/iib*num_sing_ij
2590 num_k_blocks = max_num_k_blocks - mod(max_num_k_blocks, ngroup)
2591 communication_volume = num_k_blocks*(2 + iib) + 3*(homo_beta*num_sing_ij - iib*num_k_blocks)
2592 my_steps = num_k_blocks + homo_beta*num_sing_ij - iib*num_k_blocks
2593 IF (communication_volume < min_communication_volume) THEN
2594 block_size = iib
2595 min_communication_volume = communication_volume
2596 communication_steps = my_steps
2597 ELSE IF (communication_volume == min_communication_volume .AND. my_steps < communication_steps) THEN
2598 block_size = iib
2599 communication_steps = my_steps
2600 END IF
2601 END DO
2602
2603 IF (unit_nr > 0) THEN
2604 WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
2605 "MO_INFO| Block size:", block_size
2606 CALL m_flush(unit_nr)
2607 END IF
2608
2609 ! Calculate number of large blocks
2610 max_num_k_blocks = homo_beta/block_size*num_sing_ij
2611 num_k_blocks = max_num_k_blocks - mod(max_num_k_blocks, ngroup)
2612
2613 total_ijk = num_k_blocks + homo_beta*num_sing_ij - num_k_blocks*block_size
2614 ALLOCATE (ijk_map(4, total_ijk))
2615 ijk_map = 0
2616 ALLOCATE (ijk_marker(homo_beta, num_sing_ij))
2617 ijk_marker = .true.
2618
2619 my_ijk = 0
2620 ijk_counter = 0
2621 ij_counter = 0
2622 DO iib = 1, homo
2623 ! diagonal elements already updated
2624 DO jjb = iib + 1, homo
2625 IF (abs(eigenval(jjb) - eigenval(iib)) >= mp2_env%ri_grad%eps_canonical) cycle
2626 ij_counter = ij_counter + 1
2627 DO kkb = 1, homo_beta - mod(homo_beta, block_size), block_size
2628 IF (ijk_counter + 1 > num_k_blocks) EXIT
2629 ijk_counter = ijk_counter + 1
2630 ijk_marker(kkb:kkb + block_size - 1, ij_counter) = .false.
2631 ijk_map(1, ijk_counter) = iib
2632 ijk_map(2, ijk_counter) = jjb
2633 ijk_map(3, ijk_counter) = kkb
2634 ijk_map(4, ijk_counter) = block_size
2635 IF (mod(ijk_counter, ngroup) == color_sub) my_ijk = my_ijk + 1
2636 END DO
2637 END DO
2638 END DO
2639 ij_counter = 0
2640 DO iib = 1, homo
2641 ! diagonal elements already updated
2642 DO jjb = iib + 1, homo
2643 IF (abs(eigenval(jjb) - eigenval(iib)) >= mp2_env%ri_grad%eps_canonical) cycle
2644 ij_counter = ij_counter + 1
2645 DO kkb = 1, homo_beta
2646 IF (ijk_marker(kkb, ij_counter)) THEN
2647 ijk_counter = ijk_counter + 1
2648 ijk_map(1, ijk_counter) = iib
2649 ijk_map(2, ijk_counter) = jjb
2650 ijk_map(3, ijk_counter) = kkb
2651 ijk_map(4, ijk_counter) = 1
2652 IF (mod(ijk_counter, ngroup) == color_sub) my_ijk = my_ijk + 1
2653 END IF
2654 END DO
2655 END DO
2656 END DO
2657
2658 DEALLOCATE (ijk_marker)
2659
2660 CALL comm_exchange%allgather(my_ijk, num_ijk)
2661 max_ijk = maxval(num_ijk)
2662
2663 END SUBROUTINE find_quasi_degenerate_ij
2664
2665END MODULE mp2_ri_gpw
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_log_handling::cp_to_string
Definition cp_log_handling.F:90
group_dist_types::get_group_dist
Definition group_dist_types.F:47
group_dist_types::release_group_dist
Definition group_dist_types.F:53
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
dgemm_counter_types
Counters to determine the performance of parallel DGEMMs.
Definition dgemm_counter_types.F:14
dgemm_counter_types::dgemm_counter_init
elemental subroutine, public dgemm_counter_init(dgemm_counter, unit_nr, print_info)
Initialize a dgemm_counter.
Definition dgemm_counter_types.F:49
dgemm_counter_types::dgemm_counter_write
subroutine, public dgemm_counter_write(dgemm_counter, para_env)
calculate and print flop rates
Definition dgemm_counter_types.F:103
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
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
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:141
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
mp2_ri_gpw
Routines to calculate RI-GPW-MP2 energy using pw.
Definition mp2_ri_gpw.F:14
mp2_ri_gpw::mp2_ri_gpw_compute_en
subroutine, public mp2_ri_gpw_compute_en(emp2_cou, emp2_ex, emp2_s, emp2_t, bib_c, mp2_env, para_env, para_env_sub, color_sub, gd_array, gd_b_virtual, eigenval, nmo, homo, dimen_ri, unit_nr, calc_forces, calc_ex)
...
Definition mp2_ri_gpw.F:75
mp2_ri_grad_util
Routines for calculating RI-MP2 gradients.
Definition mp2_ri_grad_util.F:13
mp2_ri_grad_util::complete_gamma
subroutine, public complete_gamma(mp2_env, b_ia_q, dimen_ri, homo, virtual, para_env, para_env_sub, ngroup, my_group_l_size, my_group_l_start, my_group_l_end, my_b_size, my_b_virtual_start, gd_array, gd_b_virtual, kspin)
complete the calculation of the Gamma matrices
Definition mp2_ri_grad_util.F:87
mp2_types
Types needed for MP2 calculations.
Definition mp2_types.F:14
dgemm_counter_types::dgemm_counter_type
Definition dgemm_counter_types.F:32
group_dist_types::group_dist_d1_type
Definition group_dist_types.F:35
message_passing::mp_comm_type
Definition message_passing.F:161
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:750
mp2_types::mp2_type
Definition mp2_types.F:345
mp2_types::three_dim_real_array
Definition mp2_types.F:87