Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines to calculate RI-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 : ! **************************************************************************************************
14 : MODULE 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 :
49 : CONTAINS
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 1062 : 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 354 : gd_array, gd_B_virtual, &
74 354 : 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 354 : INTEGER, ALLOCATABLE, DIMENSION(:) :: integ_group_pos2color_sub, my_B_size, &
100 354 : my_B_virtual_end, my_B_virtual_start, num_ij_pairs, sizes_array_orig, virtual
101 354 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: ij_map
102 354 : 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 354 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), TARGET :: buffer_1D
108 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
109 354 : TARGET :: local_ab, local_ba, t_ab
110 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
111 354 : TARGET :: local_i_aL, local_j_aL, Y_i_aP, Y_j_aP
112 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
113 354 : POINTER :: external_ab, external_i_aL
114 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
115 354 : 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 354 : DIMENSION(:) :: B_ia_Q
120 :
121 354 : CALL timeset(routineN, handle)
122 :
123 354 : nspins = SIZE(homo)
124 :
125 1062 : ALLOCATE (virtual(nspins))
126 796 : virtual(:) = nmo - homo(:)
127 :
128 1416 : ALLOCATE (my_B_size(nspins), my_B_virtual_start(nspins), my_B_virtual_end(nspins))
129 796 : DO ispin = 1, nspins
130 : CALL get_group_dist(gd_B_virtual(ispin), para_env_sub%mepos, &
131 796 : my_B_virtual_start(ispin), my_B_virtual_end(ispin), my_B_size(ispin))
132 : END DO
133 :
134 354 : CALL get_group_dist(gd_array, color_sub, my_group_L_start, my_group_L_end, my_group_L_size)
135 :
136 354 : 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 354 : CALL mp2_env%local_gemm_ctx%create(LOCAL_GEMM_PU_GPU)
141 354 : 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 354 : 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 354 : ranges_info_array, comm_exchange, comm_rep, num_integ_group)
158 :
159 : ! We cannot fix the tag because of the recv routine
160 354 : tag = 42
161 :
162 796 : 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 442 : my_group_L_size, ranges_info_array)
167 :
168 1326 : DO ispin = 1, jspin
169 :
170 530 : IF (unit_nr > 0) THEN
171 265 : IF (nspins == 1) THEN
172 133 : WRITE (unit_nr, *) "Start loop run"
173 132 : ELSE IF (ispin == 1 .AND. jspin == 1) THEN
174 44 : WRITE (unit_nr, *) "Start loop run alpha-alpha"
175 88 : ELSE IF (ispin == 1 .AND. jspin == 2) THEN
176 44 : WRITE (unit_nr, *) "Start loop run alpha-beta"
177 44 : ELSE IF (ispin == 2 .AND. jspin == 2) THEN
178 44 : WRITE (unit_nr, *) "Start loop run beta-beta"
179 : END IF
180 265 : CALL m_flush(unit_nr)
181 : END IF
182 :
183 530 : 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 530 : my_beta_beta_case = .FALSE.
190 530 : my_alpha_beta_case = .FALSE.
191 530 : IF (ispin /= jspin) THEN
192 88 : my_alpha_beta_case = .TRUE.
193 442 : ELSE IF (my_open_shell_SS) THEN
194 176 : IF (ispin == 2) my_beta_beta_case = .TRUE.
195 : END IF
196 :
197 530 : amp_fac = mp2_env%scale_S + mp2_env%scale_T
198 530 : 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 530 : 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 530 : 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 530 : block_size, ngroup, ij_map, color_sub, my_ij_pairs, my_open_shell_SS, unit_nr)
218 :
219 1590 : ALLOCATE (num_ij_pairs(0:comm_exchange%num_pe - 1))
220 530 : CALL comm_exchange%allgather(my_ij_pairs, num_ij_pairs)
221 :
222 1170 : 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 530 : local_j_aL, calc_forces, Y_i_aP, Y_j_aP, ispin, jspin)
227 :
228 530 : CALL timeset(routineN//"_RI_loop", handle2)
229 530 : my_Emp2_Cou = 0.0_dp
230 530 : my_Emp2_EX = 0.0_dp
231 530 : t_start = m_walltime()
232 2572 : DO ij_index = 1, max_ij_pairs
233 :
234 : ! Prediction is unreliable if we are in the first step of the loop
235 2042 : IF (unit_nr > 0 .AND. ij_index > 1) THEN
236 742 : decil = ij_index*10/max_ij_pairs
237 742 : IF (decil /= (ij_index - 1)*10/max_ij_pairs) THEN
238 701 : t_new = m_walltime()
239 701 : 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 701 : cp_to_string(decil*10)//". Minutes left: "//cp_to_string(t_new)
242 701 : CALL m_flush(unit_nr)
243 : END IF
244 : END IF
245 :
246 2042 : IF (calc_forces) THEN
247 1617 : Y_i_aP = 0.0_dp
248 1617 : Y_j_aP = 0.0_dp
249 : END IF
250 :
251 2042 : IF (ij_index <= my_ij_pairs) THEN
252 : ! We have work to do
253 1993 : ij_counter = (ij_index - MIN(1, color_sub))*ngroup + color_sub
254 1993 : my_i = ij_map(1, ij_counter)
255 1993 : my_j = ij_map(2, ij_counter)
256 1993 : my_block_size = ij_map(3, ij_counter)
257 :
258 1993 : 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 1993 : BIb_C(ispin)%array(:, :, my_i:my_i + my_block_size - 1))
261 :
262 1993 : 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 1993 : BIb_C(jspin)%array(:, :, my_j:my_j + my_block_size - 1))
265 :
266 : ! collect data from other proc
267 1993 : CALL timeset(routineN//"_comm", handle3)
268 2102 : DO proc_shift = 1, comm_exchange%num_pe - 1
269 109 : proc_send = MODULO(comm_exchange%mepos + proc_shift, comm_exchange%num_pe)
270 109 : proc_receive = MODULO(comm_exchange%mepos - proc_shift, comm_exchange%num_pe)
271 :
272 109 : send_ij_index = num_ij_pairs(proc_send)
273 :
274 109 : CALL get_group_dist(gd_array, proc_receive, sizes=rec_L_size)
275 :
276 2102 : IF (ij_index <= send_ij_index) THEN
277 : ij_counter_send = (ij_index - MIN(1, integ_group_pos2color_sub(proc_send)))*ngroup + &
278 60 : integ_group_pos2color_sub(proc_send)
279 60 : send_i = ij_map(1, ij_counter_send)
280 60 : 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 60 : buffer_1D(1:rec_L_size*my_B_size(ispin)*my_block_size)
285 60 : BI_C_rec = 0.0_dp
286 : CALL comm_exchange%sendrecv(BIb_C(ispin)%array(:, :, send_i:send_i + my_block_size - 1), &
287 60 : 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 60 : 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 60 : buffer_1D(1:INT(rec_L_size, int_8)*my_B_size(jspin)*my_block_size)
295 60 : BI_C_rec = 0.0_dp
296 : CALL comm_exchange%sendrecv(BIb_C(jspin)%array(:, :, send_j:send_j + my_block_size - 1), &
297 60 : 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 60 : 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 49 : buffer_1D(1:INT(rec_L_size, int_8)*my_B_size(ispin)*my_block_size)
308 49 : BI_C_rec = 0.0_dp
309 49 : 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 49 : 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 49 : buffer_1D(1:INT(rec_L_size, int_8)*my_B_size(jspin)*my_block_size)
317 49 : BI_C_rec = 0.0_dp
318 49 : 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 49 : BI_C_rec(:, 1:my_B_size(jspin), 1:my_block_size))
322 :
323 : END IF
324 :
325 : END DO
326 :
327 1993 : CALL timestop(handle3)
328 :
329 : ! loop over the block elements
330 3990 : DO iiB = 1, my_block_size
331 5995 : DO jjB = 1, my_block_size
332 2005 : CALL timeset(routineN//"_expansion", handle3)
333 4002 : 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 2005 : local_ab = 0.0_dp
337 574 : IF ((my_alpha_beta_case) .AND. (calc_forces)) THEN
338 502 : local_ba = 0.0_dp
339 : END IF
340 2005 : 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 2005 : my_B_size(ispin))
345 : ! Additional integrals only for alpha_beta case and forces
346 2005 : IF (my_alpha_beta_case .AND. calc_forces) THEN
347 : local_ba(my_B_virtual_start(jspin):my_B_virtual_end(jspin), :) = &
348 133101 : 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 2279 : DO proc_shift = 1, para_env_sub%num_pe - 1
352 274 : proc_send = MODULO(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
353 274 : 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 274 : rec_B_virtual_end, rec_B_size)
357 :
358 274 : external_i_aL(1:dimen_RI, 1:rec_B_size) => buffer_1D(1:INT(dimen_RI, int_8)*rec_B_size)
359 274 : external_i_aL = 0.0_dp
360 :
361 : CALL para_env_sub%sendrecv(my_local_i_aL, proc_send, &
362 274 : 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 274 : 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 2553 : 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 70 : rec_B_virtual_end, rec_B_size)
374 :
375 70 : external_i_aL(1:dimen_RI, 1:rec_B_size) => buffer_1D(1:INT(dimen_RI, int_8)*rec_B_size)
376 70 : external_i_aL = 0.0_dp
377 :
378 : CALL para_env_sub%sendrecv(my_local_j_aL, proc_send, &
379 70 : 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 70 : 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 2005 : 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 502 : CALL dgemm_counter_stop(dgemm_counter, virtual(ispin), my_B_size(ispin) + my_B_size(jspin), dimen_RI)
390 : ELSE
391 1503 : CALL dgemm_counter_stop(dgemm_counter, virtual(ispin), my_B_size(jspin), dimen_RI)
392 : END IF
393 2005 : CALL timestop(handle3)
394 :
395 : !sample peak memory
396 2005 : CALL m_memory()
397 :
398 2005 : CALL timeset(routineN//"_ener", handle3)
399 : ! calculate coulomb only MP2
400 2005 : sym_fac = 2.0_dp
401 2005 : IF (my_i == my_j) sym_fac = 1.0_dp
402 2005 : IF (my_alpha_beta_case) sym_fac = 0.5_dp
403 30928 : DO b = 1, my_B_size(jspin)
404 28923 : b_global = b + my_B_virtual_start(jspin) - 1
405 585679 : 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 583674 : Eigenval(my_i + iiB - 1, ispin) - Eigenval(my_j + jjB - 1, jspin))
409 : END DO
410 : END DO
411 :
412 2005 : IF (calc_ex) THEN
413 : ! contract integrals with orbital energies for exchange MP2 energy
414 : ! starting with local ...
415 2005 : IF (calc_forces .AND. (.NOT. my_alpha_beta_case)) t_ab = 0.0_dp
416 30328 : DO b = 1, my_B_size(ispin)
417 28323 : b_global = b + my_B_virtual_start(ispin) - 1
418 548367 : DO a = 1, my_B_size(ispin)
419 518039 : 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 518039 : Eigenval(my_i + iiB - 1, ispin) - Eigenval(my_j + jjB - 1, ispin))
423 546362 : 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 302030 : 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 2279 : DO proc_shift = 1, para_env_sub%num_pe - 1
433 274 : proc_send = MODULO(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
434 274 : 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 274 : rec_B_virtual_start, rec_B_virtual_end, rec_B_size)
438 : CALL get_group_dist(gd_B_virtual(ispin), proc_send, &
439 274 : 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 274 : buffer_1D(1:INT(rec_B_size, int_8)*my_B_size(ispin))
443 274 : 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 60536 : external_ab(1:my_B_size(ispin), 1:rec_B_size), proc_receive, tag)
447 :
448 5284 : DO b = 1, my_B_size(ispin)
449 2731 : b_global = b + my_B_virtual_start(ispin) - 1
450 30405 : DO a = 1, rec_B_size
451 27400 : 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 27400 : Eigenval(my_i + iiB - 1, ispin) - Eigenval(my_j + jjB - 1, ispin))
455 30131 : 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 9580 : 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 2005 : CALL timestop(handle3)
466 :
467 4010 : 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 1620 : 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 49 : my_block_size = 1
486 :
487 49 : 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 98 : DO proc_shift = 1, comm_exchange%num_pe - 1
491 49 : proc_send = MODULO(comm_exchange%mepos + proc_shift, comm_exchange%num_pe)
492 49 : proc_receive = MODULO(comm_exchange%mepos - proc_shift, comm_exchange%num_pe)
493 :
494 49 : send_ij_index = num_ij_pairs(proc_send)
495 :
496 98 : 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 49 : integ_group_pos2color_sub(proc_send)
500 49 : send_i = ij_map(1, ij_counter_send)
501 49 : 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 49 : proc_send, tag)
506 : ! occupied j
507 : CALL comm_exchange%send(BIb_C(jspin)%array(:, :, send_j:send_j + my_block_size - 1), &
508 49 : proc_send, tag)
509 : END IF
510 : END DO
511 49 : CALL timestop(handle3)
512 : END IF
513 :
514 : ! redistribute gamma
515 2572 : 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 1617 : 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 1617 : gd_array%sizes, 2, buffer_1D)
526 : END IF
527 :
528 : END DO
529 530 : CALL timestop(handle2)
530 :
531 530 : DEALLOCATE (local_i_aL)
532 530 : DEALLOCATE (local_j_aL)
533 530 : DEALLOCATE (ij_map)
534 530 : DEALLOCATE (num_ij_pairs)
535 530 : DEALLOCATE (local_ab)
536 :
537 530 : IF (calc_forces) THEN
538 384 : DEALLOCATE (Y_i_aP)
539 384 : DEALLOCATE (Y_j_aP)
540 384 : IF (ALLOCATED(t_ab)) THEN
541 306 : DEALLOCATE (t_ab)
542 : END IF
543 384 : 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 384 : integ_group_pos2color_sub, dgemm_counter, buffer_1D)
559 :
560 : END IF
561 :
562 530 : 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 530 : IF (calc_forces .AND. jspin == nspins) THEN
568 1068 : IF (.NOT. ALLOCATED(B_ia_Q)) ALLOCATE (B_ia_Q(nspins))
569 1530 : ALLOCATE (B_ia_Q(ispin)%array(homo(ispin), my_B_size(ispin), my_group_L_size_orig))
570 907278 : B_ia_Q(ispin)%array = 0.0_dp
571 1408 : DO jjB = 1, homo(ispin)
572 17332 : DO iiB = 1, my_B_size(ispin)
573 : B_ia_Q(ispin)%array(jjB, iiB, 1:my_group_L_size_orig) = &
574 724326 : BIb_C(ispin)%array(1:my_group_L_size_orig, iiB, jjB)
575 : END DO
576 : END DO
577 306 : DEALLOCATE (BIb_C(ispin)%array)
578 :
579 : ! sum Gamma and dereplicate
580 1530 : ALLOCATE (BIb_C(ispin)%array(my_B_size(ispin), homo(ispin), my_group_L_size_orig))
581 582 : DO proc_shift = 1, comm_rep%num_pe - 1
582 : ! invert order
583 276 : proc_send = MODULO(comm_rep%mepos - proc_shift, comm_rep%num_pe)
584 276 : proc_receive = MODULO(comm_rep%mepos + proc_shift, comm_rep%num_pe)
585 :
586 276 : start_point = ranges_info_array(3, proc_shift, comm_exchange%mepos)
587 276 : 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 276 : proc_send, BIb_C(ispin)%array, proc_receive, tag)
591 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) &
592 582 : !$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 770352 : BIb_C(ispin)%array(:, :, :) = mp2_env%ri_grad%Gamma_P_ia(ispin)%array(:, :, 1:my_group_L_size_orig)
600 306 : DEALLOCATE (mp2_env%ri_grad%Gamma_P_ia(ispin)%array)
601 306 : CALL MOVE_ALLOC(BIb_C(ispin)%array, mp2_env%ri_grad%Gamma_P_ia(ispin)%array)
602 224 : ELSE IF (jspin == nspins) THEN
603 136 : DEALLOCATE (BIb_C(ispin)%array)
604 : END IF
605 :
606 530 : CALL para_env%sum(my_Emp2_Cou)
607 530 : CALL para_env%sum(my_Emp2_Ex)
608 :
609 1502 : IF (my_open_shell_SS .OR. my_alpha_beta_case) THEN
610 264 : IF (my_alpha_beta_case) THEN
611 88 : Emp2_S = Emp2_S + my_Emp2_Cou
612 88 : Emp2_Cou = Emp2_Cou + my_Emp2_Cou
613 : ELSE
614 176 : my_Emp2_Cou = my_Emp2_Cou*0.25_dp
615 176 : my_Emp2_EX = my_Emp2_EX*0.5_dp
616 176 : Emp2_T = Emp2_T + my_Emp2_Cou + my_Emp2_EX
617 176 : Emp2_Cou = Emp2_Cou + my_Emp2_Cou
618 176 : Emp2_EX = Emp2_EX + my_Emp2_EX
619 : END IF
620 : ELSE
621 266 : Emp2_Cou = Emp2_Cou + my_Emp2_Cou
622 266 : Emp2_EX = Emp2_EX + my_Emp2_EX
623 : END IF
624 : END DO
625 :
626 : END DO
627 :
628 354 : DEALLOCATE (integ_group_pos2color_sub)
629 354 : DEALLOCATE (ranges_info_array)
630 :
631 354 : CALL comm_exchange%free()
632 354 : CALL comm_rep%free()
633 :
634 354 : IF (calc_forces) THEN
635 : ! recover original information (before replication)
636 228 : DEALLOCATE (gd_array%sizes)
637 228 : iiB = SIZE(sizes_array_orig)
638 684 : ALLOCATE (gd_array%sizes(0:iiB - 1))
639 678 : gd_array%sizes(:) = sizes_array_orig
640 228 : DEALLOCATE (sizes_array_orig)
641 :
642 : ! Remove replication from BIb_C and reorder the matrix
643 228 : my_group_L_size = my_group_L_size_orig
644 :
645 : ! B_ia_Q(ispin)%array will be deallocated inside of complete_gamma
646 534 : 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 534 : ispin)
653 : END DO
654 534 : DEALLOCATE (B_ia_Q)
655 :
656 45074 : 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 228 : CALL comm%from_split(para_env, para_env_sub%mepos)
660 534 : DO ispin = 1, nspins
661 : ! P_ab is only replicated over all subgroups
662 306 : CALL comm%sum(mp2_env%ri_grad%P_ab(ispin)%array)
663 : ! P_ij is replicated over all processes
664 534 : CALL para_env%sum(mp2_env%ri_grad%P_ij(ispin)%array)
665 : END DO
666 456 : CALL comm%free()
667 : END BLOCK
668 : END IF
669 :
670 354 : CALL release_group_dist(gd_array)
671 796 : DO ispin = 1, nspins
672 442 : IF (ALLOCATED(BIb_C(ispin)%array)) DEALLOCATE (BIb_C(ispin)%array)
673 796 : 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 354 : IF (calc_forces) DEALLOCATE (mp2_env%ri_grad%PQ_half)
678 354 : IF (calc_forces .AND. .NOT. compare_potential_types(mp2_env%ri_metric, mp2_env%potential_parameter)) THEN
679 8 : DEALLOCATE (mp2_env%ri_grad%operator_half)
680 : END IF
681 :
682 354 : 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 354 : CALL mp2_env%local_gemm_ctx%destroy()
686 354 : CALL timestop(handle)
687 :
688 1062 : 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 4360 : 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 4360 : CALL timeset(routineN, handle)
707 :
708 12100 : DO irep = 1, SIZE(ranges_info_array, 2)
709 7740 : Lstart_pos = ranges_info_array(1, irep)
710 7740 : Lend_pos = ranges_info_array(2, irep)
711 7740 : start_point = ranges_info_array(3, irep)
712 7740 : end_point = ranges_info_array(4, irep)
713 :
714 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) &
715 12100 : !$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 4360 : CALL timestop(handle)
721 :
722 4360 : 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 266 : 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 266 : CALL timeset(routineN, handle)
741 :
742 766 : DO irep = 1, SIZE(ranges_info_array, 2)
743 500 : Lstart_pos = ranges_info_array(1, irep)
744 500 : Lend_pos = ranges_info_array(2, irep)
745 500 : start_point = ranges_info_array(3, irep)
746 500 : end_point = ranges_info_array(4, irep)
747 :
748 : !$OMP PARALLEL WORKSHARE DEFAULT(NONE) &
749 766 : !$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 266 : CALL timestop(handle)
755 :
756 266 : 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 442 : SUBROUTINE replicate_iaK_2intgroup(BIb_C, comm_exchange, comm_rep, homo, sizes_array, my_B_size, &
770 442 : 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 442 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: BIb_C_gather
785 :
786 442 : CALL timeset(routineN, handle)
787 :
788 : ! replication scheme using mpi_allgather
789 : ! get the max L size of the
790 988 : max_L_size = MAXVAL(sizes_array)
791 :
792 2210 : ALLOCATE (BIb_C_copy(max_L_size, my_B_size, homo))
793 442 : BIb_C_copy = 0.0_dp
794 898474 : BIb_C_copy(1:SIZE(BIb_C, 1), 1:my_B_size, 1:homo) = BIb_C
795 :
796 442 : DEALLOCATE (BIb_C)
797 :
798 2652 : ALLOCATE (BIb_C_gather(max_L_size, my_B_size, homo, 0:comm_rep%num_pe - 1))
799 442 : BIb_C_gather = 0.0_dp
800 :
801 442 : CALL comm_rep%allgather(BIb_C_copy, BIb_C_gather)
802 :
803 442 : DEALLOCATE (BIb_C_copy)
804 :
805 2210 : ALLOCATE (BIb_C(my_group_L_size, my_B_size, homo))
806 442 : BIb_C = 0.0_dp
807 :
808 : ! reorder data
809 1204 : DO proc_shift = 0, comm_rep%num_pe - 1
810 762 : proc_receive = MODULO(comm_rep%mepos - proc_shift, comm_rep%num_pe)
811 :
812 762 : start_point = ranges_info_array(3, proc_shift, comm_exchange%mepos)
813 762 : 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 1674771 : BIb_C_gather(1:end_point - start_point + 1, 1:my_B_size, 1:homo, proc_receive)
817 :
818 : END DO
819 :
820 442 : DEALLOCATE (BIb_C_gather)
821 :
822 442 : CALL timestop(handle)
823 :
824 442 : 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 530 : 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 530 : CALL timeset(routineN, handle)
857 :
858 2120 : ALLOCATE (local_ab(virtual(ispin), my_B_size(jspin)))
859 530 : local_ab = 0.0_dp
860 :
861 530 : IF (calc_forces) THEN
862 384 : IF (.NOT. ALLOCATED(mp2_env%ri_grad%P_ij(jspin)%array)) THEN
863 1224 : ALLOCATE (mp2_env%ri_grad%P_ij(jspin)%array(homo(ispin), homo(ispin)))
864 6278 : mp2_env%ri_grad%P_ij(jspin)%array = 0.0_dp
865 : END IF
866 384 : IF (.NOT. ALLOCATED(mp2_env%ri_grad%P_ab(jspin)%array)) THEN
867 1224 : ALLOCATE (mp2_env%ri_grad%P_ab(jspin)%array(my_B_size(jspin), virtual(jspin)))
868 85070 : mp2_env%ri_grad%P_ab(jspin)%array = 0.0_dp
869 : END IF
870 384 : IF (.NOT. ALLOCATED(mp2_env%ri_grad%Gamma_P_ia(jspin)%array)) THEN
871 1530 : ALLOCATE (mp2_env%ri_grad%Gamma_P_ia(jspin)%array(my_B_size(jspin), homo(jspin), my_group_L_size))
872 1475502 : mp2_env%ri_grad%Gamma_P_ia(jspin)%array = 0.0_dp
873 : END IF
874 :
875 384 : IF (ispin == jspin) THEN
876 : ! For non-alpha-beta case we need amplitudes
877 918 : 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 306 : ALLOCATE (local_ba(1, 1))
881 : ELSE
882 : ! We need more integrals
883 312 : ALLOCATE (local_ba(virtual(jspin), my_B_size(ispin)))
884 : END IF
885 : END IF
886 : !
887 :
888 530 : CALL timestop(handle)
889 :
890 530 : 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 530 : 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 530 : CALL timeset(routineN, handle)
921 :
922 2650 : ALLOCATE (local_i_aL(dimen_RI, my_B_size(ispin), block_size))
923 530 : local_i_aL = 0.0_dp
924 2650 : ALLOCATE (local_j_aL(dimen_RI, my_B_size(jspin), block_size))
925 530 : local_j_aL = 0.0_dp
926 :
927 530 : IF (calc_forces) THEN
928 1536 : ALLOCATE (Y_i_aP(my_B_size(ispin), dimen_RI, block_size))
929 384 : 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 1536 : ALLOCATE (Y_j_aP(my_B_size(jspin), dimen_RI, block_size))
933 384 : Y_j_aP = 0.0_dp
934 : END IF
935 : !
936 :
937 530 : CALL timestop(handle)
938 :
939 530 : 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 530 : 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 530 : LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: ij_marker
972 :
973 : ! Calculate the maximum number of ij pairs that have to be computed
974 : ! among groups
975 :
976 530 : CALL timeset(routineN, handle)
977 :
978 530 : IF (.NOT. my_open_shell_ss .AND. .NOT. my_alpha_beta_case) THEN
979 266 : total_ij_pairs = homo*(1 + homo)/2
980 266 : num_IJ_blocks = homo/block_size - 1
981 :
982 266 : first_I_block = 1
983 266 : last_i_block = block_size*(num_IJ_blocks - 1)
984 :
985 266 : first_J_block = block_size + 1
986 266 : last_J_block = block_size*(num_IJ_blocks + 1)
987 :
988 266 : ij_block_counter = 0
989 602 : DO iiB = first_I_block, last_i_block, block_size
990 602 : DO jjB = iiB + block_size, last_J_block, block_size
991 840 : ij_block_counter = ij_block_counter + 1
992 : END DO
993 : END DO
994 :
995 266 : total_ij_block = ij_block_counter
996 266 : num_block_per_group = total_ij_block/ngroup
997 266 : assigned_blocks = num_block_per_group*ngroup
998 :
999 266 : total_ij_pairs_blocks = assigned_blocks + (total_ij_pairs - assigned_blocks*(block_size**2))
1000 :
1001 1064 : ALLOCATE (ij_marker(homo, homo))
1002 4018 : ij_marker = .TRUE.
1003 798 : ALLOCATE (ij_map(3, total_ij_pairs_blocks))
1004 266 : ij_map = 0
1005 266 : ij_counter = 0
1006 266 : my_ij_pairs = 0
1007 602 : DO iiB = first_I_block, last_i_block, block_size
1008 1278 : DO jjB = iiB + block_size, last_J_block, block_size
1009 840 : IF (ij_counter + 1 > assigned_blocks) EXIT
1010 676 : ij_counter = ij_counter + 1
1011 2028 : ij_marker(iiB:iiB + block_size - 1, jjB:jjB + block_size - 1) = .FALSE.
1012 676 : ij_map(1, ij_counter) = iiB
1013 676 : ij_map(2, ij_counter) = jjB
1014 676 : ij_map(3, ij_counter) = block_size
1015 1012 : IF (MOD(ij_counter, ngroup) == color_sub) my_ij_pairs = my_ij_pairs + 1
1016 : END DO
1017 : END DO
1018 1070 : DO iiB = 1, homo
1019 2946 : DO jjB = iiB, homo
1020 2680 : IF (ij_marker(iiB, jjB)) THEN
1021 1200 : ij_counter = ij_counter + 1
1022 1200 : ij_map(1, ij_counter) = iiB
1023 1200 : ij_map(2, ij_counter) = jjB
1024 1200 : ij_map(3, ij_counter) = 1
1025 1200 : IF (MOD(ij_counter, ngroup) == color_sub) my_ij_pairs = my_ij_pairs + 1
1026 : END IF
1027 : END DO
1028 : END DO
1029 266 : DEALLOCATE (ij_marker)
1030 :
1031 264 : 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 176 : total_ij_pairs = homo*(homo - 1)/2
1035 176 : num_IJ_blocks = (homo - 1)/block_size - 1
1036 :
1037 176 : first_I_block = 1
1038 176 : 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 176 : first_J_block = block_size + 2
1042 176 : last_J_block = block_size*(num_IJ_blocks + 1) + 1
1043 :
1044 176 : ij_block_counter = 0
1045 262 : DO iiB = first_I_block, last_i_block, block_size
1046 262 : DO jjB = iiB + block_size + 1, last_J_block, block_size
1047 200 : ij_block_counter = ij_block_counter + 1
1048 : END DO
1049 : END DO
1050 :
1051 176 : total_ij_block = ij_block_counter
1052 176 : num_block_per_group = total_ij_block/ngroup
1053 176 : assigned_blocks = num_block_per_group*ngroup
1054 :
1055 176 : total_ij_pairs_blocks = assigned_blocks + (total_ij_pairs - assigned_blocks*(block_size**2))
1056 :
1057 704 : ALLOCATE (ij_marker(homo, homo))
1058 2984 : ij_marker = .TRUE.
1059 528 : ALLOCATE (ij_map(3, total_ij_pairs_blocks))
1060 176 : ij_map = 0
1061 176 : ij_counter = 0
1062 176 : my_ij_pairs = 0
1063 262 : DO iiB = first_I_block, last_i_block, block_size
1064 458 : DO jjB = iiB + block_size + 1, last_J_block, block_size
1065 200 : IF (ij_counter + 1 > assigned_blocks) EXIT
1066 196 : ij_counter = ij_counter + 1
1067 604 : ij_marker(iiB:iiB + block_size - 1, jjB:jjB + block_size - 1) = .FALSE.
1068 196 : ij_map(1, ij_counter) = iiB
1069 196 : ij_map(2, ij_counter) = jjB
1070 196 : ij_map(3, ij_counter) = block_size
1071 282 : IF (MOD(ij_counter, ngroup) == color_sub) my_ij_pairs = my_ij_pairs + 1
1072 : END DO
1073 : END DO
1074 774 : DO iiB = 1, homo
1075 1580 : DO jjB = iiB + 1, homo
1076 1404 : IF (ij_marker(iiB, jjB)) THEN
1077 598 : ij_counter = ij_counter + 1
1078 598 : ij_map(1, ij_counter) = iiB
1079 598 : ij_map(2, ij_counter) = jjB
1080 598 : ij_map(3, ij_counter) = 1
1081 598 : IF (MOD(ij_counter, ngroup) == color_sub) my_ij_pairs = my_ij_pairs + 1
1082 : END IF
1083 : END DO
1084 : END DO
1085 176 : DEALLOCATE (ij_marker)
1086 :
1087 : ELSE
1088 88 : total_ij_pairs = homo*homo_beta
1089 88 : num_IJ_blocks = homo/block_size
1090 88 : num_IJ_blocks_beta = homo_beta/block_size
1091 :
1092 88 : first_I_block = 1
1093 88 : last_i_block = block_size*(num_IJ_blocks - 1)
1094 :
1095 88 : first_J_block = 1
1096 88 : last_J_block = block_size*(num_IJ_blocks_beta - 1)
1097 :
1098 88 : ij_block_counter = 0
1099 246 : DO iiB = first_I_block, last_i_block, block_size
1100 246 : DO jjB = first_J_block, last_J_block, block_size
1101 284 : ij_block_counter = ij_block_counter + 1
1102 : END DO
1103 : END DO
1104 :
1105 88 : total_ij_block = ij_block_counter
1106 88 : num_block_per_group = total_ij_block/ngroup
1107 88 : assigned_blocks = num_block_per_group*ngroup
1108 :
1109 88 : total_ij_pairs_blocks = assigned_blocks + (total_ij_pairs - assigned_blocks*(block_size**2))
1110 :
1111 352 : ALLOCATE (ij_marker(homo, homo_beta))
1112 1396 : ij_marker = .TRUE.
1113 264 : ALLOCATE (ij_map(3, total_ij_pairs_blocks))
1114 88 : ij_map = 0
1115 88 : ij_counter = 0
1116 88 : my_ij_pairs = 0
1117 246 : DO iiB = first_I_block, last_i_block, block_size
1118 486 : DO jjB = first_J_block, last_J_block, block_size
1119 244 : IF (ij_counter + 1 > assigned_blocks) EXIT
1120 240 : ij_counter = ij_counter + 1
1121 720 : ij_marker(iiB:iiB + block_size - 1, jjB:jjB + block_size - 1) = .FALSE.
1122 240 : ij_map(1, ij_counter) = iiB
1123 240 : ij_map(2, ij_counter) = jjB
1124 240 : ij_map(3, ij_counter) = block_size
1125 398 : IF (MOD(ij_counter, ngroup) == color_sub) my_ij_pairs = my_ij_pairs + 1
1126 : END DO
1127 : END DO
1128 432 : DO iiB = 1, homo
1129 1486 : DO jjB = 1, homo_beta
1130 1398 : IF (ij_marker(iiB, jjB)) THEN
1131 814 : ij_counter = ij_counter + 1
1132 814 : ij_map(1, ij_counter) = iiB
1133 814 : ij_map(2, ij_counter) = jjB
1134 814 : ij_map(3, ij_counter) = 1
1135 814 : IF (MOD(ij_counter, ngroup) == color_sub) my_ij_pairs = my_ij_pairs + 1
1136 : END IF
1137 : END DO
1138 : END DO
1139 88 : DEALLOCATE (ij_marker)
1140 : END IF
1141 :
1142 530 : IF (unit_nr > 0) THEN
1143 265 : IF (block_size == 1) THEN
1144 : WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.1)") &
1145 236 : "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 29 : "RI_INFO| Percentage of ij pairs communicated with block size 1:", &
1149 58 : 100.0_dp*REAL((total_ij_pairs - assigned_blocks*(block_size**2)), KIND=dp)/REAL(total_ij_pairs, KIND=dp)
1150 : END IF
1151 265 : CALL m_flush(unit_nr)
1152 : END IF
1153 :
1154 530 : CALL timestop(handle)
1155 :
1156 530 : 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 354 : 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 354 : INTEGER, ALLOCATABLE, DIMENSION(:) :: new_sizes_array, rep_ends_array, &
1205 : rep_sizes_array, rep_starts_array
1206 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: my_info
1207 :
1208 354 : CALL timeset(routineN, handle)
1209 : !
1210 354 : sub_sub_color = para_env_sub%mepos*num_integ_group + color_sub/integ_group_size
1211 354 : CALL comm_exchange%from_split(para_env, sub_sub_color)
1212 :
1213 : ! create the replication group
1214 354 : sub_sub_color = para_env_sub%mepos*comm_exchange%num_pe + comm_exchange%mepos
1215 354 : 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 1062 : ALLOCATE (rep_ends_array(0:comm_rep%num_pe - 1))
1222 708 : ALLOCATE (rep_starts_array(0:comm_rep%num_pe - 1))
1223 708 : ALLOCATE (rep_sizes_array(0:comm_rep%num_pe - 1))
1224 :
1225 354 : CALL comm_rep%allgather(my_group_L_size, rep_sizes_array)
1226 354 : CALL comm_rep%allgather(my_group_L_start, rep_starts_array)
1227 354 : CALL comm_rep%allgather(my_group_L_end, rep_ends_array)
1228 :
1229 : ! calculate my_new_group_L_size according to sizes_array
1230 354 : my_new_group_L_size = my_group_L_size
1231 :
1232 : ! Info of this process
1233 1062 : ALLOCATE (my_info(4, 0:comm_rep%num_pe - 1))
1234 354 : my_info(1, 0) = my_group_L_start
1235 354 : my_info(2, 0) = my_group_L_end
1236 354 : my_info(3, 0) = 1
1237 354 : my_info(4, 0) = my_group_L_size
1238 :
1239 596 : DO proc_shift = 1, comm_rep%num_pe - 1
1240 242 : proc_receive = MODULO(comm_rep%mepos - proc_shift, comm_rep%num_pe)
1241 :
1242 242 : my_new_group_L_size = my_new_group_L_size + rep_sizes_array(proc_receive)
1243 :
1244 242 : my_info(1, proc_shift) = rep_starts_array(proc_receive)
1245 242 : my_info(2, proc_shift) = rep_ends_array(proc_receive)
1246 242 : my_info(3, proc_shift) = my_info(4, proc_shift - 1) + 1
1247 596 : my_info(4, proc_shift) = my_new_group_L_size
1248 :
1249 : END DO
1250 :
1251 1062 : ALLOCATE (new_sizes_array(0:comm_exchange%num_pe - 1))
1252 1416 : ALLOCATE (ranges_info_array(4, 0:comm_rep%num_pe - 1, 0:comm_exchange%num_pe - 1))
1253 354 : CALL comm_exchange%allgather(my_new_group_L_size, new_sizes_array)
1254 354 : CALL comm_exchange%allgather(my_info, ranges_info_array)
1255 :
1256 354 : DEALLOCATE (rep_sizes_array)
1257 354 : DEALLOCATE (rep_starts_array)
1258 354 : DEALLOCATE (rep_ends_array)
1259 :
1260 1062 : ALLOCATE (integ_group_pos2color_sub(0:comm_exchange%num_pe - 1))
1261 354 : CALL comm_exchange%allgather(color_sub, integ_group_pos2color_sub)
1262 :
1263 354 : IF (calc_forces) THEN
1264 228 : iiB = SIZE(sizes_array)
1265 684 : ALLOCATE (sizes_array_orig(0:iiB - 1))
1266 678 : sizes_array_orig(:) = sizes_array
1267 : END IF
1268 :
1269 354 : my_group_L_size_orig = my_group_L_size
1270 354 : my_group_L_size = my_new_group_L_size
1271 354 : DEALLOCATE (sizes_array)
1272 :
1273 1062 : ALLOCATE (sizes_array(0:integ_group_size - 1))
1274 806 : sizes_array(:) = new_sizes_array
1275 :
1276 354 : DEALLOCATE (new_sizes_array)
1277 : !
1278 354 : CALL timestop(handle)
1279 :
1280 708 : 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 1416 : SUBROUTINE mp2_ri_get_integ_group_size(mp2_env, para_env, para_env_sub, gd_array, gd_B_virtual, &
1299 354 : homo, dimen_RI, unit_nr, &
1300 : integ_group_size, &
1301 : ngroup, num_integ_group, &
1302 354 : 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 354 : CALL timeset(routineN, handle)
1325 :
1326 354 : ngroup = para_env%num_pe/para_env_sub%num_pe
1327 :
1328 354 : calc_group_size = mp2_env%ri_mp2%number_integration_groups <= 0
1329 354 : IF (.NOT. calc_group_size) THEN
1330 10 : IF (MOD(ngroup, mp2_env%ri_mp2%number_integration_groups) /= 0) calc_group_size = .TRUE.
1331 : END IF
1332 :
1333 354 : IF (calc_group_size) THEN
1334 344 : CALL m_memory(mem)
1335 344 : mem_real = (mem + 1024*1024 - 1)/(1024*1024)
1336 344 : CALL para_env%min(mem_real)
1337 344 : mem_real = MIN(mem_real, mp2_env%mp2_memory)
1338 :
1339 344 : mem_base = 0.0_dp
1340 344 : mem_per_blk = 0.0_dp
1341 344 : mem_per_repl = 0.0_dp
1342 344 : 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 770 : maxsize(gd_B_virtual))*8.0_dp/(1024**2)
1347 : ! BIB_C
1348 770 : 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 770 : 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 770 : 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 770 : 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 1196 : mem_base = mem_base + REAL(MAX(dimen_RI, MAXVAL(virtual)), KIND=dp)*MAXVAL(maxsize(gd_B_virtual))*8.0_dp/(1024**2)
1357 :
1358 344 : IF (calc_forces) THEN
1359 : ! Gamma_P_ia
1360 : mem_per_repl = mem_per_repl + SUM(REAL(homo, KIND=dp)*maxsize(gd_array)* &
1361 514 : maxsize(gd_B_virtual))*8.0_dp/(1024**2)
1362 : ! Y_i_aP+Y_j_aP
1363 514 : 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 808 : 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 514 : mem_base = mem_base + SUM(REAL(homo, KIND=dp)*homo)*8.0_dp/(1024**2)
1368 : ! P_ab
1369 514 : 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 808 : 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 770 : block_size = MAX(1, MIN(FLOOR(SQRT(REAL(MINVAL(homo), KIND=dp))), FLOOR(MINVAL(homo)/SQRT(2.0_dp*ngroup))))
1376 344 : IF (mp2_env%ri_mp2%block_size > 0) block_size = mp2_env%ri_mp2%block_size
1377 :
1378 344 : mem_min = mem_base + mem_per_repl + (mem_per_blk + mem_per_repl_blk)*block_size
1379 :
1380 516 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T68,F9.2,A4)') 'RI_INFO| Minimum available memory per MPI process:', &
1381 344 : mem_real, ' MiB'
1382 516 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T68,F9.2,A4)') 'RI_INFO| Minimum required memory per MPI process:', &
1383 344 : 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 344 : 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 1622 : - SUM((REAL(MAXVAL(homo), KIND=dp)/block_size + block_size - 2.0_dp)*homo*virtual)/ngroup
1409 672 : IF (SIZE(homo) == 2) factor = factor - 2.0_dp*PRODUCT(homo)/block_size/ngroup*SUM(homo*virtual)
1410 :
1411 688 : 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 256 : (mem_per_repl + mem_per_repl_blk*block_size)), 1), ngroup)
1415 : ! Convert to an integration group size
1416 256 : 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 256 : DO iiB = MAX(MIN(min_integ_group_size, ngroup), 1), ngroup
1420 : ! check that the ngroup is a multiple of integ_group_size
1421 256 : IF (MOD(ngroup, iiB) == 0) THEN
1422 256 : integ_group_size = iiB
1423 256 : EXIT
1424 : END IF
1425 0 : integ_group_size = integ_group_size + 1
1426 : END DO
1427 : END IF
1428 : ELSE ! We take the user provided group size
1429 10 : integ_group_size = ngroup/mp2_env%ri_mp2%number_integration_groups
1430 : END IF
1431 :
1432 354 : IF (unit_nr > 0) THEN
1433 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
1434 177 : "RI_INFO| Group size for integral replication:", integ_group_size*para_env_sub%num_pe
1435 177 : CALL m_flush(unit_nr)
1436 : END IF
1437 :
1438 354 : num_integ_group = ngroup/integ_group_size
1439 :
1440 354 : CALL timestop(handle)
1441 :
1442 354 : 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 530 : SUBROUTINE mp2_ri_get_block_size(mp2_env, para_env, para_env_sub, gd_array, gd_B_virtual, &
1463 530 : 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 530 : CALL timeset(routineN, handle)
1486 :
1487 530 : ngroup = para_env%num_pe/para_env_sub%num_pe
1488 :
1489 530 : CALL m_memory(mem)
1490 530 : mem_real = (mem + 1024*1024 - 1)/(1024*1024)
1491 530 : CALL para_env%min(mem_real)
1492 :
1493 530 : mem_base = 0.0_dp
1494 530 : mem_per_blk = 0.0_dp
1495 530 : mem_per_repl_blk = 0.0_dp
1496 :
1497 : ! external_ab
1498 1766 : mem_base = mem_base + MAXVAL(maxsize(gd_B_virtual))*MAX(dimen_RI, MAXVAL(virtual))*8.0_dp/(1024**2)
1499 : ! BIB_C_rec
1500 1148 : 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 1148 : 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 1766 : mem_base = mem_base + MAXVAL(maxsize(gd_B_virtual))*MAX(dimen_RI, MAXVAL(virtual))*8.0_dp/(1024**2)
1505 :
1506 530 : IF (calc_forces) THEN
1507 : ! Y_i_aP+Y_j_aP+BIb_C_send
1508 846 : 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 1308 : mem_base = mem_base + MAXVAL(maxsize(gd_B_virtual))*MAX(dimen_RI, MAXVAL(virtual))*8.0_dp/(1024**2)
1511 : END IF
1512 :
1513 530 : best_block_size = 1
1514 :
1515 : ! Here we split the memory half for the communication, half for replication
1516 530 : IF (mp2_env%ri_mp2%block_size > 0) THEN
1517 : best_block_size = mp2_env%ri_mp2%block_size
1518 : ELSE
1519 294 : best_block_size = MAX(FLOOR((mem_real - mem_base)/(mem_per_blk + mem_per_repl_blk*ngroup/num_integ_group)), 1)
1520 :
1521 2988966 : DO
1522 2989260 : IF (SIZE(homo) == 1) THEN
1523 2405280 : IF (.NOT. my_open_shell_ss) THEN
1524 1106908 : num_IJ_blocks = (homo(1)/best_block_size)
1525 1106908 : num_IJ_blocks = (num_IJ_blocks*num_IJ_blocks - num_IJ_blocks)/2
1526 : ELSE
1527 1298372 : num_IJ_blocks = ((homo(1) - 1)/best_block_size)
1528 1298372 : num_IJ_blocks = (num_IJ_blocks*num_IJ_blocks - num_IJ_blocks)/2
1529 : END IF
1530 : ELSE
1531 1751940 : num_ij_blocks = PRODUCT(homo/best_block_size)
1532 : END IF
1533 : ! Enforce at least one large block for each subgroup
1534 2989260 : IF ((num_IJ_blocks >= ngroup .AND. num_IJ_blocks > 0) .OR. best_block_size == 1) THEN
1535 : EXIT
1536 : ELSE
1537 2988966 : best_block_size = best_block_size - 1
1538 : END IF
1539 : END DO
1540 :
1541 294 : IF (SIZE(homo) == 1) THEN
1542 232 : 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 124 : 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 108 : best_block_size = MIN(FLOOR(SQRT(REAL(homo(1), KIND=dp))), best_block_size)
1549 : END IF
1550 : END IF
1551 : END IF
1552 530 : block_size = MAX(1, best_block_size)
1553 :
1554 530 : IF (unit_nr > 0) THEN
1555 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
1556 265 : "RI_INFO| Block size:", block_size
1557 265 : 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 1766 : *MAXVAL(maxsize(gd_B_virtual))
1563 : ! The send buffer has the same size as the recv buffer
1564 530 : IF (calc_forces) buffer_size = buffer_size*2
1565 1590 : ALLOCATE (buffer_1D(buffer_size))
1566 :
1567 530 : CALL timestop(handle)
1568 :
1569 530 : 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 6480 : SUBROUTINE mp2_update_P_gamma(mp2_env, para_env_sub, gd_B_virtual, &
1601 3240 : Eigenval, homo, dimen_RI, iiB, jjB, my_B_size, &
1602 3240 : my_B_virtual_end, my_B_virtual_start, my_i, my_j, virtual, local_ab, &
1603 3240 : t_ab, my_local_i_aL, my_local_j_aL, open_ss, Y_i_aP, Y_j_aP, &
1604 1620 : 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 1620 : POINTER :: external_ab, send_ab
1636 :
1637 1620 : CALL timeset(routineN//"_Pia", handle)
1638 :
1639 1620 : alpha_beta = .NOT. (ispin == jspin)
1640 1620 : IF (open_ss) THEN
1641 : factor = 1.0_dp
1642 : ELSE
1643 1231 : factor = 2.0_dp
1644 : END IF
1645 : ! divide the (ia|jb) integrals by Delta_ij^ab
1646 25159 : DO b = 1, my_B_size(jspin)
1647 23539 : b_global = b + my_B_virtual_start(jspin) - 1
1648 469848 : 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 468228 : Eigenval(my_i + iiB - 1, ispin) - Eigenval(my_j + jjB - 1, jspin))
1652 : END DO
1653 : END DO
1654 1620 : IF (.NOT. (alpha_beta)) THEN
1655 329242 : P_ij_diag = -SUM(local_ab*t_ab)*factor
1656 : ELSE
1657 : ! update diagonal part of P_ij
1658 140606 : P_ij_diag = -SUM(local_ab*local_ab)*mp2_env%scale_S
1659 : ! More integrals needed only for alpha-beta case: local_ba
1660 6987 : DO b = 1, my_B_size(ispin)
1661 6485 : b_global = b + my_B_virtual_start(ispin) - 1
1662 140066 : 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 139564 : 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 1620 : CALL dgemm_counter_start(dgemm_counter)
1673 1620 : 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 1118 : 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 1118 : 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 502 : 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 502 : 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 502 : 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 502 : 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 1620 : 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 823 : 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 823 : mp2_env%ri_grad%P_ij(ispin)%array(my_j + jjB - 1, my_j + jjB - 1) + P_ij_diag
1707 : END IF
1708 1792 : DO proc_shift = 1, para_env_sub%num_pe - 1
1709 172 : proc_send = MODULO(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
1710 172 : proc_receive = MODULO(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
1711 :
1712 172 : CALL get_group_dist(gd_B_virtual(jspin), proc_receive, rec_B_virtual_start, rec_B_virtual_end, rec_B_size)
1713 172 : CALL get_group_dist(gd_B_virtual(jspin), proc_send, send_B_virtual_start, send_B_virtual_end, send_B_size)
1714 :
1715 172 : external_ab(1:virtual(ispin), 1:rec_B_size) => buffer_1D(1:INT(virtual(ispin), int_8)*rec_B_size)
1716 172 : external_ab = 0.0_dp
1717 :
1718 : CALL para_env_sub%sendrecv(local_ab, proc_send, &
1719 172 : external_ab, proc_receive)
1720 :
1721 172 : 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 102 : 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 70 : 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 70 : CALL get_group_dist(gd_B_virtual(ispin), proc_receive, rec_B_virtual_start, rec_B_virtual_end, rec_B_size)
1734 70 : CALL get_group_dist(gd_B_virtual(ispin), proc_send, send_B_virtual_start, send_B_virtual_end, send_B_size)
1735 70 : external_ab(1:virtual(jspin), 1:rec_B_size) => buffer_1D(1:INT(virtual(jspin), int_8)*rec_B_size)
1736 70 : external_ab = 0.0_dp
1737 : CALL para_env_sub%sendrecv(local_ba, proc_send, &
1738 70 : 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 70 : 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 1964 : IF ((my_i /= my_j) .AND. (.NOT. alpha_beta)) THEN
1745 : external_ab(1:my_B_size(ispin), 1:virtual(ispin)) => &
1746 86 : buffer_1D(1:INT(virtual(ispin), int_8)*my_B_size(ispin))
1747 86 : external_ab = 0.0_dp
1748 :
1749 86 : offset = INT(virtual(ispin), int_8)*my_B_size(ispin)
1750 :
1751 86 : send_ab(1:send_B_size, 1:virtual(ispin)) => buffer_1D(offset + 1:offset + INT(send_B_size, int_8)*virtual(ispin))
1752 86 : 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 86 : local_ab, virtual(ispin), 0.0_dp, send_ab, send_B_size)
1757 : CALL para_env_sub%sendrecv(send_ab, proc_send, &
1758 86 : external_ab, proc_receive)
1759 :
1760 18083 : 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 1620 : IF (.NOT. alpha_beta) THEN
1765 1118 : IF (my_i /= my_j) THEN
1766 823 : CALL dgemm_counter_stop(dgemm_counter, 2*my_B_size(ispin), virtual(ispin), virtual(ispin))
1767 : ELSE
1768 295 : CALL dgemm_counter_stop(dgemm_counter, my_B_size(ispin), virtual(ispin), virtual(ispin))
1769 : END IF
1770 : ELSE
1771 1506 : CALL dgemm_counter_stop(dgemm_counter, SUM(my_B_size), virtual(ispin), virtual(jspin))
1772 : END IF
1773 1620 : CALL timestop(handle)
1774 :
1775 : ! Now, Gamma_P_ia (made of Y_ia_P)
1776 :
1777 1620 : CALL timeset(routineN//"_Gamma", handle)
1778 1620 : CALL dgemm_counter_start(dgemm_counter)
1779 1620 : 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 1118 : 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 502 : 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 502 : my_local_i_aL, dimen_RI, 1.0_dp, Y_j_aP, my_B_size(jspin))
1791 : END IF
1792 :
1793 1620 : IF (para_env_sub%num_pe > 1) THEN
1794 172 : external_ab(1:my_B_size(ispin), 1:dimen_RI) => buffer_1D(1:INT(my_B_size(ispin), int_8)*dimen_RI)
1795 172 : external_ab = 0.0_dp
1796 :
1797 172 : offset = INT(my_B_size(ispin), int_8)*dimen_RI
1798 : END IF
1799 : !
1800 1792 : DO proc_shift = 1, para_env_sub%num_pe - 1
1801 172 : proc_send = MODULO(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
1802 172 : proc_receive = MODULO(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
1803 :
1804 172 : CALL get_group_dist(gd_B_virtual(ispin), proc_receive, rec_B_virtual_start, rec_B_virtual_end, rec_B_size)
1805 172 : CALL get_group_dist(gd_B_virtual(ispin), proc_send, send_B_virtual_start, send_B_virtual_end, send_B_size)
1806 :
1807 172 : send_ab(1:send_B_size, 1:dimen_RI) => buffer_1D(offset + 1:offset + INT(dimen_RI, int_8)*send_B_size)
1808 172 : send_ab = 0.0_dp
1809 1792 : 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 102 : my_local_j_aL, dimen_RI, 0.0_dp, send_ab, send_B_size)
1813 102 : CALL para_env_sub%sendrecv(send_ab, proc_send, external_ab, proc_receive)
1814 :
1815 217544 : 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 70 : my_local_j_aL, dimen_RI, 0.0_dp, send_ab, send_B_size)
1822 70 : CALL para_env_sub%sendrecv(send_ab, proc_send, external_ab, proc_receive)
1823 161840 : Y_i_aP(:, :) = Y_i_aP + external_ab
1824 : END IF
1825 : END DO
1826 :
1827 1620 : IF (alpha_beta) THEN
1828 : ! For beta-beta part (in alpha-beta case) we need a new parallel code
1829 502 : IF (para_env_sub%num_pe > 1) THEN
1830 70 : external_ab(1:my_B_size(jspin), 1:dimen_RI) => buffer_1D(1:INT(my_B_size(jspin), int_8)*dimen_RI)
1831 70 : external_ab = 0.0_dp
1832 :
1833 70 : offset = INT(my_B_size(jspin), int_8)*dimen_RI
1834 : END IF
1835 572 : DO proc_shift = 1, para_env_sub%num_pe - 1
1836 70 : proc_send = MODULO(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
1837 70 : proc_receive = MODULO(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
1838 :
1839 70 : CALL get_group_dist(gd_B_virtual(jspin), proc_send, send_B_virtual_start, send_B_virtual_end, send_B_size)
1840 70 : send_ab(1:send_B_size, 1:dimen_RI) => buffer_1D(offset + 1:offset + INT(dimen_RI, int_8)*send_B_size)
1841 70 : 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 70 : my_local_i_aL, dimen_RI, 0.0_dp, send_ab, send_B_size)
1845 70 : CALL para_env_sub%sendrecv(send_ab, proc_send, external_ab, proc_receive)
1846 177812 : 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 502 : CALL dgemm_counter_stop(dgemm_counter, 3*virtual(ispin), dimen_RI, my_B_size(jspin))
1852 : ELSE
1853 1118 : CALL dgemm_counter_stop(dgemm_counter, virtual(ispin), dimen_RI, my_B_size(ispin))
1854 : END IF
1855 :
1856 1620 : IF ((my_i /= my_j) .AND. (.NOT. alpha_beta)) THEN
1857 : ! Alpha-alpha, beta-beta and closed shell
1858 823 : 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 823 : my_local_i_aL, dimen_RI, 1.0_dp, Y_j_aP, my_B_size(ispin))
1862 909 : DO proc_shift = 1, para_env_sub%num_pe - 1
1863 86 : proc_send = MODULO(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
1864 86 : proc_receive = MODULO(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
1865 :
1866 86 : CALL get_group_dist(gd_B_virtual(ispin), proc_receive, rec_B_virtual_start, rec_B_virtual_end, rec_B_size)
1867 :
1868 86 : external_ab(1:dimen_RI, 1:rec_B_size) => buffer_1D(1:INT(dimen_RI, int_8)*rec_B_size)
1869 86 : external_ab = 0.0_dp
1870 :
1871 : CALL para_env_sub%sendrecv(my_local_i_aL, proc_send, &
1872 86 : 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 995 : external_ab, dimen_RI, 1.0_dp, Y_j_aP, my_B_size(ispin))
1878 : END DO
1879 :
1880 823 : CALL dgemm_counter_stop(dgemm_counter, my_B_size(ispin), dimen_RI, virtual(ispin))
1881 : END IF
1882 :
1883 1620 : CALL timestop(handle)
1884 1620 : 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 3234 : 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 3234 : ij_map, ranges_info_array, Y_i_aP, comm_exchange, &
1911 3234 : 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 3234 : 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 3234 : CALL timeset(routineN//"_comm2", handle)
1938 :
1939 3234 : tag = 43
1940 :
1941 3234 : IF (ij_index <= my_ij_pairs) THEN
1942 : ! somethig to send
1943 : ! start with myself
1944 3216 : CALL timeset(routineN//"_comm2_w", handle2)
1945 9230 : DO irep = 0, num_integ_group - 1
1946 6014 : Lstart_pos = ranges_info_array(1, irep, comm_exchange%mepos)
1947 6014 : start_point = ranges_info_array(3, irep, comm_exchange%mepos)
1948 6014 : 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 9230 : !$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 3216 : CALL timestop(handle2)
1964 :
1965 : ! Y_i_aP(my_B_size,dimen_RI,block_size)
1966 :
1967 3314 : DO proc_shift = 1, comm_exchange%num_pe - 1
1968 98 : proc_send = MODULO(comm_exchange%mepos + proc_shift, comm_exchange%num_pe)
1969 98 : proc_receive = MODULO(comm_exchange%mepos - proc_shift, comm_exchange%num_pe)
1970 :
1971 98 : send_L_size = sizes_array(proc_send)
1972 : BI_C_send(1:my_B_size, 1:my_block_size, 1:send_L_size) => &
1973 98 : buffer_1D(1:INT(my_B_size, int_8)*my_block_size*send_L_size)
1974 :
1975 98 : offset = INT(my_B_size, int_8)*my_block_size*send_L_size
1976 :
1977 98 : CALL timeset(routineN//"_comm2_w", handle2)
1978 98 : BI_C_send = 0.0_dp
1979 196 : DO irep = 0, num_integ_group - 1
1980 98 : Lstart_pos = ranges_info_array(1, irep, proc_send)
1981 98 : start_point = ranges_info_array(3, irep, proc_send)
1982 98 : 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 196 : !$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 98 : CALL timestop(handle2)
1996 :
1997 98 : rec_ij_index = num_ij_pairs(proc_receive)
1998 :
1999 3412 : 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 80 : (ij_index - MIN(1, integ_group_pos2color_sub(proc_receive)))*ngroup + integ_group_pos2color_sub(proc_receive)
2003 :
2004 80 : 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 80 : buffer_1D(offset + 1:offset + INT(my_B_size, int_8)*my_block_size*my_group_L_size)
2008 80 : BI_C_rec = 0.0_dp
2009 :
2010 : CALL comm_exchange%sendrecv(BI_C_send, proc_send, &
2011 80 : BI_C_rec, proc_receive, tag)
2012 :
2013 80 : CALL timeset(routineN//"_comm2_w", handle2)
2014 160 : DO irep = 0, num_integ_group - 1
2015 80 : start_point = ranges_info_array(3, irep, comm_exchange%mepos)
2016 80 : 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 160 : !$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 80 : CALL timestop(handle2)
2026 :
2027 : ELSE
2028 : ! we have something to send but nothing to receive
2029 18 : 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 36 : DO proc_shift = 1, comm_exchange%num_pe - 1
2038 18 : proc_send = MODULO(comm_exchange%mepos + proc_shift, comm_exchange%num_pe)
2039 18 : proc_receive = MODULO(comm_exchange%mepos - proc_shift, comm_exchange%num_pe)
2040 18 : rec_ij_index = num_ij_pairs(proc_receive)
2041 :
2042 36 : 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 18 : (ij_index - MIN(1, integ_group_pos2color_sub(proc_receive)))*ngroup + integ_group_pos2color_sub(proc_receive)
2046 :
2047 18 : 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 18 : buffer_1D(1:INT(my_B_size, int_8)*my_block_size*my_group_L_size)
2051 :
2052 18 : BI_C_rec = 0.0_dp
2053 :
2054 18 : CALL comm_exchange%recv(BI_C_rec, proc_receive, tag)
2055 :
2056 18 : CALL timeset(routineN//"_comm2_w", handle2)
2057 36 : DO irep = 0, num_integ_group - 1
2058 18 : start_point = ranges_info_array(3, irep, comm_exchange%mepos)
2059 18 : 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 36 : !$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 18 : CALL timestop(handle2)
2073 :
2074 : END IF
2075 : END DO
2076 :
2077 : END IF
2078 3234 : CALL timestop(handle)
2079 :
2080 3234 : 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 384 : SUBROUTINE quasi_degenerate_P_ij(mp2_env, Eigenval, homo, virtual, open_shell, &
2110 384 : beta_beta, Bib_C, unit_nr, dimen_RI, &
2111 384 : my_B_size, ngroup, my_group_L_size, &
2112 : color_sub, ranges_info_array, comm_exchange, para_env_sub, para_env, &
2113 384 : my_B_virtual_start, my_B_virtual_end, sizes_array, gd_B_virtual, &
2114 384 : 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 384 : INTEGER, ALLOCATABLE, DIMENSION(:) :: num_ijk
2144 384 : 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 384 : TARGET :: local_ab, local_aL_i, local_aL_j, t_ab
2151 384 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: local_aL_k
2152 384 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: BI_C_rec, external_ab, external_aL
2153 384 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: BI_C_rec_3D
2154 :
2155 384 : CALL timeset(routineN//"_ij_sing", handle)
2156 :
2157 384 : tag = 44
2158 384 : tag2 = 45
2159 :
2160 384 : nspins = SIZE(BIb_C)
2161 384 : alpha_beta = (nspins == 2)
2162 :
2163 : ! Set amplitude factor
2164 384 : amp_fac = mp2_env%scale_S + mp2_env%scale_T
2165 384 : IF (open_shell) amp_fac = mp2_env%scale_T
2166 :
2167 794 : ALLOCATE (send_last_k(2, comm_exchange%num_pe - 1))
2168 :
2169 : ! Loop(s) over orbital triplets
2170 846 : DO ispin = 1, nspins
2171 462 : size_B_i = my_B_size(ispin)
2172 462 : IF (ispin == 1 .AND. alpha_beta) THEN
2173 : kspin = 2
2174 : ELSE
2175 384 : kspin = 1
2176 : END IF
2177 462 : 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 618 : SIZE(buffer_1D), my_group_L_size, size_B_k, para_env, virtual(ispin), size_B_i)
2184 :
2185 462 : my_virtual = virtual(ispin)
2186 462 : IF (SIZE(ijk_map, 2) > 0) THEN
2187 98 : max_block_size = ijk_map(4, 1)
2188 : ELSE
2189 : max_block_size = 1
2190 : END IF
2191 :
2192 1848 : ALLOCATE (local_aL_i(dimen_RI, size_B_i))
2193 1386 : ALLOCATE (local_aL_j(dimen_RI, size_B_i))
2194 2310 : ALLOCATE (local_aL_k(dimen_RI, size_B_k, max_block_size))
2195 1848 : ALLOCATE (t_ab(my_virtual, size_B_k))
2196 :
2197 1386 : my_last_k = -1
2198 552 : send_last_k = -1
2199 :
2200 462 : t_start = m_walltime()
2201 618 : DO ijk_index = 1, max_ijk
2202 :
2203 : ! Prediction is unreliable if we are in the first step of the loop
2204 156 : IF (unit_nr > 0 .AND. ijk_index > 1) THEN
2205 18 : decil = ijk_index*10/max_ijk
2206 18 : IF (decil /= (ijk_index - 1)*10/max_ijk) THEN
2207 18 : t_new = m_walltime()
2208 18 : 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 18 : cp_to_string(decil*10)//". Minutes left: "//cp_to_string(t_new)
2211 18 : CALL m_flush(unit_nr)
2212 : END IF
2213 : END IF
2214 :
2215 618 : IF (ijk_index <= my_ijk) THEN
2216 : ! work to be done
2217 154 : ijk_counter = (ijk_index - MIN(1, color_sub))*ngroup + color_sub
2218 154 : my_i = ijk_map(1, ijk_counter)
2219 154 : my_j = ijk_map(2, ijk_counter)
2220 154 : my_k = ijk_map(3, ijk_counter)
2221 154 : block_size = ijk_map(4, ijk_counter)
2222 :
2223 154 : do_recv_i = (ispin /= kspin) .OR. my_i < my_k .OR. my_i > my_k + block_size - 1
2224 154 : do_recv_j = (ispin /= kspin) .OR. my_j < my_k .OR. my_j > my_k + block_size - 1
2225 154 : do_recv_k = my_k /= my_last_k(1) .OR. my_k + block_size - 1 /= my_last_k(2)
2226 154 : my_last_k(1) = my_k
2227 154 : my_last_k(2) = my_k + block_size - 1
2228 :
2229 154 : local_aL_i = 0.0_dp
2230 154 : IF (do_recv_i) THEN
2231 : CALL fill_local_i_aL_2D(local_al_i, ranges_info_array(:, :, comm_exchange%mepos), &
2232 125 : BIb_C(ispin)%array(:, :, my_i))
2233 : END IF
2234 :
2235 154 : local_aL_j = 0.0_dp
2236 154 : IF (do_recv_j) THEN
2237 : CALL fill_local_i_aL_2D(local_al_j, ranges_info_array(:, :, comm_exchange%mepos), &
2238 125 : BIb_C(ispin)%array(:, :, my_j))
2239 : END IF
2240 :
2241 154 : IF (do_recv_k) THEN
2242 146 : 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 146 : BIb_C(kspin)%array(:, :, my_k:my_k + block_size - 1))
2245 : END IF
2246 :
2247 154 : CALL timeset(routineN//"_comm", handle2)
2248 164 : DO proc_shift = 1, comm_exchange%num_pe - 1
2249 10 : proc_send = MODULO(comm_exchange%mepos + proc_shift, comm_exchange%num_pe)
2250 10 : proc_receive = MODULO(comm_exchange%mepos - proc_shift, comm_exchange%num_pe)
2251 :
2252 10 : send_ijk_index = num_ijk(proc_send)
2253 :
2254 10 : rec_L_size = sizes_array(proc_receive)
2255 10 : 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 10 : do_send_i = .FALSE.
2258 10 : do_send_j = .FALSE.
2259 10 : do_send_k = .FALSE.
2260 10 : 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 8 : ngroup + integ_group_pos2color_sub(proc_send)
2264 8 : send_i = ijk_map(1, ijk_counter_send)
2265 8 : send_j = ijk_map(2, ijk_counter_send)
2266 8 : send_k = ijk_map(3, ijk_counter_send)
2267 :
2268 8 : do_send_i = (ispin /= kspin) .OR. send_i < send_k .OR. send_i > send_k + block_size - 1
2269 8 : do_send_j = (ispin /= kspin) .OR. send_j < send_k .OR. send_j > send_k + block_size - 1
2270 8 : do_send_k = send_k /= send_last_k(1, proc_shift) .OR. send_k + block_size - 1 /= send_last_k(2, proc_shift)
2271 8 : send_last_k(1, proc_shift) = send_k
2272 8 : send_last_k(2, proc_shift) = send_k + block_size - 1
2273 : END IF
2274 :
2275 : ! occupied i
2276 722 : BI_C_rec = 0.0_dp
2277 10 : IF (do_send_i) THEN
2278 6 : IF (do_recv_i) THEN
2279 : CALL comm_exchange%sendrecv(BIb_C(ispin)%array(:, :, send_i), proc_send, &
2280 288 : BI_C_rec, proc_receive, tag)
2281 : ELSE
2282 2 : CALL comm_exchange%send(BIb_C(ispin)%array(:, :, send_i), proc_send, tag)
2283 : END IF
2284 4 : ELSE IF (do_recv_i) THEN
2285 580 : CALL comm_exchange%recv(BI_C_rec, proc_receive, tag)
2286 : END IF
2287 10 : IF (do_recv_i) THEN
2288 8 : 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 722 : BI_C_rec = 0.0_dp
2293 10 : IF (do_send_j) THEN
2294 8 : IF (do_recv_j) THEN
2295 : CALL comm_exchange%sendrecv(BIb_C(ispin)%array(:, :, send_j), proc_send, &
2296 576 : BI_C_rec, proc_receive, tag)
2297 : ELSE
2298 0 : CALL comm_exchange%send(BIb_C(ispin)%array(:, :, send_j), proc_send, tag)
2299 : END IF
2300 2 : ELSE IF (do_recv_j) THEN
2301 0 : CALL comm_exchange%recv(BI_C_rec, proc_receive, tag)
2302 : END IF
2303 8 : IF (do_recv_j) THEN
2304 8 : 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 10 : buffer_1D(1:INT(rec_L_size, KIND=int_8)*size_B_k*block_size)
2310 10 : IF (do_send_k) THEN
2311 8 : IF (do_recv_k) THEN
2312 : CALL comm_exchange%sendrecv(BIb_C(kspin)%array(:, :, send_k:send_k + block_size - 1), proc_send, &
2313 726 : BI_C_rec_3D, proc_receive, tag)
2314 : ELSE
2315 0 : CALL comm_exchange%send(BI_C_rec, proc_receive, tag)
2316 : END IF
2317 2 : ELSE IF (do_recv_k) THEN
2318 294 : CALL comm_exchange%recv(BI_C_rec_3D, proc_receive, tag)
2319 : END IF
2320 164 : IF (do_recv_k) THEN
2321 10 : 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 20868 : IF (.NOT. do_recv_i) local_aL_i(:, :) = local_aL_k(:, :, my_i - my_k + 1)
2326 20868 : IF (.NOT. do_recv_j) local_aL_j(:, :) = local_aL_k(:, :, my_j - my_k + 1)
2327 154 : CALL timestop(handle2)
2328 :
2329 : ! expand integrals
2330 376 : DO kkB = 1, block_size
2331 222 : CALL timeset(routineN//"_exp_ik", handle2)
2332 222 : CALL dgemm_counter_start(dgemm_counter)
2333 666 : ALLOCATE (local_ab(my_virtual, size_B_k))
2334 222 : 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 222 : 0.0_dp, local_ab(my_B_virtual_start(ispin):my_B_virtual_end(ispin), 1:size_B_k), size_B_i)
2338 222 : DO proc_shift = 1, para_env_sub%num_pe - 1
2339 0 : proc_send = MODULO(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
2340 0 : proc_receive = MODULO(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
2341 :
2342 0 : CALL get_group_dist(gd_B_virtual(ispin), proc_receive, rec_B_virtual_start, rec_B_virtual_end, rec_B_size)
2343 :
2344 0 : 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 0 : 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 222 : 0.0_dp, local_ab(rec_B_virtual_start:rec_B_virtual_end, 1:size_B_k), rec_B_size)
2352 : END DO
2353 222 : CALL dgemm_counter_stop(dgemm_counter, my_virtual, size_B_k, dimen_RI)
2354 222 : CALL timestop(handle2)
2355 :
2356 : ! Amplitudes
2357 222 : CALL timeset(routineN//"_tab", handle2)
2358 222 : t_ab = 0.0_dp
2359 : ! Alpha-alpha, beta-beta and closed shell
2360 222 : IF (.NOT. alpha_beta) THEN
2361 1156 : DO b = 1, size_B_k
2362 1038 : b_global = b + my_B_virtual_start(1) - 1
2363 17250 : DO a = 1, my_B_size(1)
2364 16094 : 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 17132 : - Eigenval(homo(1) + a_global, 1) - Eigenval(homo(1) + b_global, 1))
2368 : END DO
2369 : END DO
2370 : ELSE
2371 1030 : DO b = 1, size_B_k
2372 926 : b_global = b + my_B_virtual_start(kspin) - 1
2373 15388 : DO a = 1, my_B_size(ispin)
2374 14358 : 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 15284 : - Eigenval(homo(ispin) + a_global, ispin) - Eigenval(homo(kspin) + b_global, kspin))
2378 : END DO
2379 : END DO
2380 : END IF
2381 :
2382 222 : IF (.NOT. alpha_beta) THEN
2383 118 : DO proc_shift = 1, para_env_sub%num_pe - 1
2384 0 : proc_send = MODULO(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
2385 0 : proc_receive = MODULO(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
2386 0 : CALL get_group_dist(gd_B_virtual(1), proc_receive, rec_B_virtual_start, rec_B_virtual_end, rec_B_size)
2387 0 : CALL get_group_dist(gd_B_virtual(1), proc_send, send_B_virtual_start, send_B_virtual_end, send_B_size)
2388 :
2389 0 : 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 0 : external_ab(1:size_B_i, 1:rec_B_size), proc_receive, tag)
2392 :
2393 118 : DO b = 1, my_B_size(1)
2394 0 : b_global = b + my_B_virtual_start(1) - 1
2395 0 : DO a = 1, rec_B_size
2396 0 : 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 0 : - 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 222 : CALL timestop(handle2)
2405 :
2406 : ! Expand the second set of integrals
2407 222 : CALL timeset(routineN//"_exp_jk", handle2)
2408 222 : local_ab = 0.0_dp
2409 222 : 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 222 : 0.0_dp, local_ab(my_B_virtual_start(ispin):my_B_virtual_end(ispin), 1:size_B_k), size_B_i)
2413 222 : DO proc_shift = 1, para_env_sub%num_pe - 1
2414 0 : proc_send = MODULO(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
2415 0 : proc_receive = MODULO(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
2416 :
2417 0 : CALL get_group_dist(gd_B_virtual(ispin), proc_receive, rec_B_virtual_start, rec_B_virtual_end, rec_B_size)
2418 :
2419 0 : 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 0 : 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 222 : 0.0_dp, local_ab(rec_B_virtual_start:rec_B_virtual_end, 1:size_B_k), rec_B_size)
2426 : END DO
2427 222 : CALL dgemm_counter_stop(dgemm_counter, my_virtual, size_B_k, dimen_RI)
2428 222 : CALL timestop(handle2)
2429 :
2430 222 : CALL timeset(routineN//"_Pij", handle2)
2431 2186 : DO b = 1, size_B_k
2432 1964 : b_global = b + my_B_virtual_start(kspin) - 1
2433 32638 : DO a = 1, my_B_size(ispin)
2434 30452 : 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 32416 : - Eigenval(homo(ispin) + a_global, ispin) - Eigenval(homo(kspin) + b_global, kspin))
2438 : END DO
2439 : END DO
2440 : !
2441 32638 : P_ij_elem = SUM(local_ab*t_ab)
2442 222 : DEALLOCATE (local_ab)
2443 222 : IF ((.NOT. open_shell) .AND. (.NOT. alpha_beta)) THEN
2444 4 : P_ij_elem = P_ij_elem*2.0_dp
2445 : END IF
2446 222 : IF (beta_beta) THEN
2447 34 : 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 34 : 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 188 : 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 188 : 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 1264 : CALL timestop(handle2)
2454 : END DO
2455 : ELSE
2456 2 : CALL timeset(routineN//"_comm", handle2)
2457 : ! no work to be done, possible messeges to be exchanged
2458 4 : DO proc_shift = 1, comm_exchange%num_pe - 1
2459 2 : proc_send = MODULO(comm_exchange%mepos + proc_shift, comm_exchange%num_pe)
2460 2 : proc_receive = MODULO(comm_exchange%mepos - proc_shift, comm_exchange%num_pe)
2461 :
2462 2 : send_ijk_index = num_ijk(proc_send)
2463 :
2464 4 : 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 2 : integ_group_pos2color_sub(proc_send)
2468 2 : send_i = ijk_map(1, ijk_counter_send)
2469 2 : send_j = ijk_map(2, ijk_counter_send)
2470 2 : send_k = ijk_map(3, ijk_counter_send)
2471 2 : block_size = ijk_map(4, ijk_counter_send)
2472 :
2473 2 : do_send_i = (ispin /= kspin) .OR. send_i < send_k .OR. send_i > send_k + block_size - 1
2474 2 : do_send_j = (ispin /= kspin) .OR. send_j < send_k .OR. send_j > send_k + block_size - 1
2475 : ! occupied i
2476 2 : IF (do_send_i) THEN
2477 2 : CALL comm_exchange%send(BIb_C(ispin)%array(:, :, send_i), proc_send, tag)
2478 : END IF
2479 : ! occupied j
2480 2 : IF (do_send_j) THEN
2481 0 : CALL comm_exchange%send(BIb_C(ispin)%array(:, :, send_j), proc_send, tag)
2482 : END IF
2483 : ! occupied k
2484 2 : 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 2 : CALL timestop(handle2)
2489 : END IF
2490 : END DO ! ijk_index loop
2491 462 : DEALLOCATE (local_aL_i)
2492 462 : DEALLOCATE (local_aL_j)
2493 462 : DEALLOCATE (local_aL_k)
2494 462 : DEALLOCATE (t_ab)
2495 846 : DEALLOCATE (ijk_map)
2496 : END DO ! over number of loops (ispin)
2497 384 : CALL timestop(handle)
2498 :
2499 768 : 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 462 : 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 462 : LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: ijk_marker
2547 :
2548 1386 : ALLOCATE (num_ijk(0:comm_exchange%num_pe - 1))
2549 :
2550 462 : num_sing_ij = 0
2551 2092 : DO iiB = 1, homo
2552 : ! diagonal elements already updated
2553 4368 : DO jjB = iiB + 1, homo
2554 3906 : IF (ABS(Eigenval(jjB) - Eigenval(iiB)) < mp2_env%ri_grad%eps_canonical) THEN
2555 114 : num_sing_ij = num_sing_ij + 1
2556 : END IF
2557 : END DO
2558 : END DO
2559 :
2560 462 : IF (unit_nr > 0) THEN
2561 231 : IF (do_print_alpha) THEN
2562 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
2563 153 : "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 78 : "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 462 : max_block_size = buffer_size/(my_group_L_size*B_size_k)
2572 :
2573 : ! Second limit: memory
2574 462 : CALL m_memory(mem)
2575 : ! Convert to number of doubles
2576 462 : mem = mem/8
2577 : ! Remove local_ab (2x) and local_aL_i (2x)
2578 462 : mem = mem - 2_int_8*(virtual*B_size_k + B_size_i*my_group_L_size)
2579 462 : 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 462 : 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 462 : block_size = 1
2586 462 : min_communication_volume = 3*homo_beta*num_sing_ij
2587 462 : communication_steps = 3*homo_beta*num_sing_ij
2588 1174 : DO iiB = max_block_size, 2, -1
2589 712 : max_num_k_blocks = homo_beta/iiB*num_sing_ij
2590 712 : num_k_blocks = max_num_k_blocks - MOD(max_num_k_blocks, ngroup)
2591 712 : communication_volume = num_k_blocks*(2 + iiB) + 3*(homo_beta*num_sing_ij - iiB*num_k_blocks)
2592 712 : my_steps = num_k_blocks + homo_beta*num_sing_ij - iiB*num_k_blocks
2593 1174 : IF (communication_volume < min_communication_volume) THEN
2594 52 : block_size = iiB
2595 52 : min_communication_volume = communication_volume
2596 52 : communication_steps = my_steps
2597 660 : ELSE IF (communication_volume == min_communication_volume .AND. my_steps < communication_steps) THEN
2598 58 : block_size = iiB
2599 58 : communication_steps = my_steps
2600 : END IF
2601 : END DO
2602 :
2603 462 : IF (unit_nr > 0) THEN
2604 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
2605 231 : "MO_INFO| Block size:", block_size
2606 231 : CALL m_flush(unit_nr)
2607 : END IF
2608 :
2609 : ! Calculate number of large blocks
2610 462 : max_num_k_blocks = homo_beta/block_size*num_sing_ij
2611 462 : num_k_blocks = max_num_k_blocks - MOD(max_num_k_blocks, ngroup)
2612 :
2613 462 : total_ijk = num_k_blocks + homo_beta*num_sing_ij - num_k_blocks*block_size
2614 1022 : ALLOCATE (ijk_map(4, total_ijk))
2615 462 : ijk_map = 0
2616 1484 : ALLOCATE (ijk_marker(homo_beta, num_sing_ij))
2617 1020 : ijk_marker = .TRUE.
2618 :
2619 462 : my_ijk = 0
2620 462 : ijk_counter = 0
2621 462 : ij_counter = 0
2622 2092 : DO iiB = 1, homo
2623 : ! diagonal elements already updated
2624 4368 : DO jjB = iiB + 1, homo
2625 2276 : IF (ABS(Eigenval(jjB) - Eigenval(iiB)) >= mp2_env%ri_grad%eps_canonical) CYCLE
2626 114 : ij_counter = ij_counter + 1
2627 1880 : DO kkB = 1, homo_beta - MOD(homo_beta, block_size), block_size
2628 172 : IF (ijk_counter + 1 > num_k_blocks) EXIT
2629 136 : ijk_counter = ijk_counter + 1
2630 408 : ijk_marker(kkB:kkB + block_size - 1, ij_counter) = .FALSE.
2631 136 : ijk_map(1, ijk_counter) = iiB
2632 136 : ijk_map(2, ijk_counter) = jjB
2633 136 : ijk_map(3, ijk_counter) = kkB
2634 136 : ijk_map(4, ijk_counter) = block_size
2635 2412 : 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 2092 : DO iiB = 1, homo
2641 : ! diagonal elements already updated
2642 4368 : DO jjB = iiB + 1, homo
2643 2276 : IF (ABS(Eigenval(jjB) - Eigenval(iiB)) >= mp2_env%ri_grad%eps_canonical) CYCLE
2644 114 : ij_counter = ij_counter + 1
2645 2188 : DO kkB = 1, homo_beta
2646 2720 : IF (ijk_marker(kkB, ij_counter)) THEN
2647 172 : ijk_counter = ijk_counter + 1
2648 172 : ijk_map(1, ijk_counter) = iiB
2649 172 : ijk_map(2, ijk_counter) = jjB
2650 172 : ijk_map(3, ijk_counter) = kkB
2651 172 : ijk_map(4, ijk_counter) = 1
2652 172 : 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 462 : DEALLOCATE (ijk_marker)
2659 :
2660 462 : CALL comm_exchange%allgather(my_ijk, num_ijk)
2661 954 : max_ijk = MAXVAL(num_ijk)
2662 :
2663 462 : END SUBROUTINE Find_quasi_degenerate_ij
2664 :
2665 : END MODULE mp2_ri_gpw
|
