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 MP2 energy
10 : !> \par History
11 : !> 06.2011 created [Mauro Del Ben]
12 : !> \author Mauro Del Ben
13 : ! **************************************************************************************************
14 : MODULE mp2_direct_method
15 : USE cell_types, ONLY: cell_type
16 : USE cp_log_handling, ONLY: cp_get_default_logger,&
17 : cp_logger_type
18 : USE four_center_integral_interface, ONLY: four_center_integral_context_type
19 : USE four_center_integral_methods, ONLY: coulomb4
20 : USE hfx_load_balance_methods, ONLY: cost_model,&
21 : p1_energy,&
22 : p2_energy,&
23 : p3_energy
24 : USE hfx_pair_list_methods, ONLY: build_pair_list_mp2
25 : USE hfx_types, ONLY: hfx_basis_type,&
26 : hfx_pgf_list,&
27 : hfx_pgf_product_list,&
28 : hfx_potential_type,&
29 : hfx_screen_coeff_type,&
30 : hfx_type,&
31 : pair_set_list_type
32 : USE input_constants, ONLY: do_potential_TShPSC
33 : USE kinds, ONLY: dp,&
34 : int_8
35 : USE machine, ONLY: m_flush
36 : USE mathconstants, ONLY: zero
37 : USE message_passing, ONLY: mp_para_env_release,&
38 : mp_para_env_type
39 : USE mp2_ri_libint, ONLY: prepare_integral_calc
40 : USE mp2_types, ONLY: init_TShPSC_lmax,&
41 : mp2_biel_type,&
42 : mp2_type,&
43 : pair_list_type_mp2
44 : USE orbital_pointers, ONLY: ncoset
45 : USE particle_types, ONLY: particle_type
46 : USE qs_environment_types, ONLY: qs_environment_type
47 : USE t_sh_p_s_c, ONLY: free_C0
48 : #include "./base/base_uses.f90"
49 :
50 : IMPLICIT NONE
51 : PRIVATE
52 :
53 : PUBLIC :: mp2_direct_energy, mp2_canonical_direct_single_batch
54 :
55 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mp2_direct_method'
56 :
57 : !***
58 :
59 : CONTAINS
60 :
61 : ! **************************************************************************************************
62 : !> \brief ...
63 : !> \param dimen ...
64 : !> \param occ_i ...
65 : !> \param occ_j ...
66 : !> \param mp2_biel ...
67 : !> \param mp2_env ...
68 : !> \param C_i ...
69 : !> \param Auto_i ...
70 : !> \param Emp2 ...
71 : !> \param Emp2_Cou ...
72 : !> \param Emp2_ex ...
73 : !> \param qs_env ...
74 : !> \param para_env ...
75 : !> \param unit_nr ...
76 : !> \param C_j ...
77 : !> \param Auto_j ...
78 : ! **************************************************************************************************
79 26 : SUBROUTINE mp2_direct_energy(dimen, occ_i, occ_j, mp2_biel, mp2_env, C_i, Auto_i, Emp2, Emp2_Cou, Emp2_ex, &
80 26 : qs_env, para_env, unit_nr, C_j, Auto_j)
81 : INTEGER :: dimen, occ_i, occ_j
82 : TYPE(mp2_biel_type) :: mp2_biel
83 : TYPE(mp2_type) :: mp2_env
84 : REAL(KIND=dp), DIMENSION(dimen, dimen) :: C_i
85 : REAL(KIND=dp), DIMENSION(dimen) :: Auto_i
86 : REAL(KIND=dp) :: Emp2, Emp2_Cou, Emp2_ex
87 : TYPE(qs_environment_type), POINTER :: qs_env
88 : TYPE(mp_para_env_type), POINTER :: para_env
89 : INTEGER :: unit_nr
90 : REAL(KIND=dp), DIMENSION(dimen, dimen), OPTIONAL :: C_j
91 : REAL(KIND=dp), DIMENSION(dimen), OPTIONAL :: Auto_j
92 :
93 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mp2_direct_energy'
94 : REAL(KIND=dp), PARAMETER :: zero = 0.D+00
95 :
96 : INTEGER :: batch_number, color_sub, counter, elements_ij_proc, group_counter, handle, i, &
97 : i_batch, i_batch_start, i_group_counter, j, j_batch_start, j_group_counter, last_batch, &
98 : max_batch_number, max_batch_size, max_set, minimum_memory_needed, my_batch_size, &
99 : my_I_batch_size, my_I_occupied_end, my_I_occupied_start, my_J_batch_size, &
100 : my_J_occupied_end, my_J_occupied_start, natom, Ni_occupied, Nj_occupied, number_groups, &
101 : number_i_subset, number_j_subset, one, sqrt_number_groups, total_I_size_batch_group, &
102 : total_J_size_batch_group, virt_i, virt_j
103 26 : INTEGER, ALLOCATABLE, DIMENSION(:) :: batch_sizes, batch_sizes_tmp, &
104 26 : vector_batch_I_size_group, &
105 26 : vector_batch_J_size_group
106 26 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: ij_list_proc, ij_list_proc_temp, &
107 26 : ij_matrix
108 : LOGICAL :: alpha_beta_case
109 : TYPE(mp_para_env_type), POINTER :: para_env_sub
110 :
111 26 : CALL timeset(routineN, handle)
112 :
113 26 : alpha_beta_case = .FALSE.
114 26 : IF (PRESENT(C_j) .AND. PRESENT(Auto_j)) alpha_beta_case = .TRUE.
115 :
116 26 : IF (unit_nr > 0 .AND. mp2_env%potential_parameter%potential_type == do_potential_TShPSC) THEN
117 1 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T64,F12.6,A5)') 'Truncated MP2 method, Rt=', &
118 2 : mp2_env%potential_parameter%cutoff_radius, ' Bohr'
119 : END IF
120 :
121 : ! create the local para env
122 : ! each para_env_sub corresponds to a group that is going to compute
123 : ! all the integrals. To each group a batch I is assigned and the
124 : ! communication takes place only inside the group
125 26 : number_groups = para_env%num_pe/mp2_env%mp2_num_proc
126 26 : IF (number_groups*mp2_env%mp2_num_proc /= para_env%num_pe) THEN
127 0 : CPABORT(" The number of processors needs to be a multiple of the processors per group. ")
128 : END IF
129 26 : IF (number_groups > occ_i*occ_j) THEN
130 2 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A)') 'Number of groups greater then the number of IJ pairs!'
131 2 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A)') 'Consider using more processors per group for improved efficiency'
132 : END IF
133 :
134 26 : color_sub = para_env%mepos/mp2_env%mp2_num_proc
135 26 : ALLOCATE (para_env_sub)
136 26 : CALL para_env_sub%from_split(para_env, color_sub)
137 :
138 : ! calculate the maximal size of the batch, according to the maximum RS size
139 26 : max_set = SIZE(mp2_biel%index_table, 2)
140 26 : minimum_memory_needed = (8*(max_set**4))/1024**2
141 26 : IF (minimum_memory_needed > mp2_env%mp2_memory) THEN
142 0 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T66,F12.6,A3)') 'Memory required below the minimum, new memory:', &
143 0 : minimum_memory_needed, 'MiB'
144 0 : mp2_env%mp2_memory = minimum_memory_needed
145 : END IF
146 :
147 : ! Distribute the batches over the groups in
148 : ! a rectangular fashion, bigger size for J index
149 : ! the sizes of the I batches should be as small as possible
150 26 : sqrt_number_groups = INT(SQRT(REAL(number_groups, KIND=dp)))
151 26 : DO i = 1, number_groups
152 26 : IF (MOD(number_groups, i) == 0) THEN
153 26 : IF (sqrt_number_groups/i <= 1) THEN
154 : number_j_subset = i
155 : EXIT
156 : END IF
157 : END IF
158 : END DO
159 26 : number_i_subset = number_groups/number_j_subset
160 :
161 26 : IF (number_i_subset < number_j_subset) THEN
162 0 : number_i_subset = number_j_subset
163 0 : number_j_subset = number_groups/number_i_subset
164 : END IF
165 :
166 : ! Distribute the I index and the J index over groups
167 26 : total_I_size_batch_group = occ_i/number_i_subset
168 26 : IF (total_I_size_batch_group < 1) total_I_size_batch_group = 1
169 78 : ALLOCATE (vector_batch_I_size_group(0:number_i_subset - 1))
170 :
171 26 : vector_batch_I_size_group = 0
172 78 : DO i = 0, number_i_subset - 1
173 78 : vector_batch_I_size_group(i) = total_I_size_batch_group
174 : END DO
175 78 : IF (SUM(vector_batch_I_size_group) /= occ_i) THEN
176 54 : one = 1
177 54 : IF (SUM(vector_batch_I_size_group) > occ_i) one = -1
178 18 : i = -1
179 : DO
180 18 : i = i + 1
181 18 : vector_batch_I_size_group(i) = vector_batch_I_size_group(i) + one
182 54 : IF (SUM(vector_batch_I_size_group) == occ_i) EXIT
183 18 : IF (i == number_i_subset - 1) i = -1
184 : END DO
185 : END IF
186 :
187 26 : total_J_size_batch_group = occ_j/number_j_subset
188 26 : IF (total_J_size_batch_group < 1) total_J_size_batch_group = 1
189 78 : ALLOCATE (vector_batch_J_size_group(0:number_j_subset - 1))
190 :
191 26 : vector_batch_J_size_group = 0
192 52 : DO i = 0, number_J_subset - 1
193 52 : vector_batch_J_size_group(i) = total_J_size_batch_group
194 : END DO
195 52 : IF (SUM(vector_batch_J_size_group) /= occ_j) THEN
196 0 : one = 1
197 0 : IF (SUM(vector_batch_J_size_group) > occ_j) one = -1
198 0 : i = -1
199 : DO
200 0 : i = i + 1
201 0 : vector_batch_J_size_group(i) = vector_batch_J_size_group(i) + one
202 0 : IF (SUM(vector_batch_J_size_group) == occ_j) EXIT
203 0 : IF (i == number_J_subset - 1) i = -1
204 : END DO
205 : END IF
206 :
207 : ! now the starting and ending I and J occupied orbitals are assigned to each group
208 26 : group_counter = 0
209 26 : i_group_counter = 0
210 26 : my_I_occupied_start = 1
211 39 : DO i = 0, number_i_subset - 1
212 : my_J_occupied_start = 1
213 : j_group_counter = 0
214 52 : DO j = 0, number_j_subset - 1
215 39 : group_counter = group_counter + 1
216 39 : IF (color_sub == group_counter - 1) EXIT
217 13 : my_J_occupied_start = my_J_occupied_start + vector_batch_J_size_group(j)
218 52 : j_group_counter = j_group_counter + 1
219 : END DO
220 39 : IF (color_sub == group_counter - 1) EXIT
221 13 : my_I_occupied_start = my_I_occupied_start + vector_batch_I_size_group(i)
222 39 : i_group_counter = i_group_counter + 1
223 : END DO
224 26 : my_I_occupied_end = my_I_occupied_start + vector_batch_I_size_group(i_group_counter) - 1
225 26 : my_I_batch_size = vector_batch_I_size_group(i_group_counter)
226 26 : my_J_occupied_end = my_J_occupied_start + vector_batch_J_size_group(j_group_counter) - 1
227 26 : my_J_batch_size = vector_batch_J_size_group(j_group_counter)
228 :
229 26 : DEALLOCATE (vector_batch_I_size_group)
230 26 : DEALLOCATE (vector_batch_J_size_group)
231 :
232 : max_batch_size = MIN( &
233 : MAX(1, &
234 : INT(mp2_env%mp2_memory*INT(1024, KIND=int_8)**2/ &
235 : (8*(2*dimen - occ_i)*INT(dimen, KIND=int_8)*my_J_batch_size/para_env_sub%num_pe))) &
236 26 : , my_I_batch_size)
237 26 : IF (max_batch_size < 1) THEN
238 1 : max_batch_size = INT((8*(occ_i + 1)*INT(dimen, KIND=int_8)**2/para_env%num_pe)/1024**2)
239 1 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T72,I6,A3)') 'More memory required, at least:', max_batch_size, 'MiB'
240 1 : max_batch_size = 1
241 : END IF
242 :
243 : ! create the size of the batches inside the group
244 26 : my_batch_size = my_I_batch_size
245 77 : ALLOCATE (batch_sizes(my_batch_size))
246 :
247 85 : batch_sizes = -HUGE(0)
248 : batch_number = 0
249 51 : DO i = 1, my_batch_size
250 45 : IF (i*max_batch_size > my_batch_size) EXIT
251 25 : batch_number = batch_number + 1
252 51 : batch_sizes(i) = max_batch_size
253 : END DO
254 26 : last_batch = my_batch_size - max_batch_size*batch_number
255 26 : IF (last_batch > 0) THEN
256 0 : batch_number = batch_number + 1
257 0 : batch_sizes(batch_number) = last_batch
258 : END IF
259 :
260 77 : ALLOCATE (batch_sizes_tmp(batch_number))
261 51 : batch_sizes_tmp(1:batch_number) = batch_sizes(1:batch_number)
262 26 : DEALLOCATE (batch_sizes)
263 51 : ALLOCATE (batch_sizes(batch_number))
264 51 : batch_sizes(:) = batch_sizes_tmp
265 26 : DEALLOCATE (batch_sizes_tmp)
266 :
267 51 : max_batch_size = MAXVAL(batch_sizes)
268 26 : CALL para_env%max(max_batch_size)
269 26 : max_batch_number = batch_number
270 26 : CALL para_env%max(max_batch_number)
271 26 : IF (unit_nr > 0) THEN
272 13 : WRITE (unit_nr, '(T3,A,T76,I5)') 'Maximum used batch size: ', max_batch_size
273 13 : WRITE (unit_nr, '(T3,A,T76,I5)') 'Number of integral recomputations: ', max_batch_number
274 13 : CALL m_flush(unit_nr)
275 : END IF
276 :
277 : ! Batches sizes exceed the occupied orbitals allocated for group
278 51 : CPASSERT(SUM(batch_sizes) <= my_batch_size)
279 :
280 26 : virt_i = dimen - occ_i
281 26 : virt_j = dimen - occ_j
282 26 : natom = SIZE(mp2_biel%index_table, 1)
283 :
284 26 : CALL para_env%sync()
285 26 : Emp2 = zero
286 26 : Emp2_Cou = zero
287 26 : Emp2_ex = zero
288 26 : i_batch_start = my_I_occupied_start - 1
289 26 : j_batch_start = my_J_occupied_start - 1
290 26 : Nj_occupied = my_J_batch_size
291 51 : DO i_batch = 1, batch_number
292 :
293 25 : Ni_occupied = batch_sizes(i_batch)
294 :
295 25 : counter = -1
296 100 : ALLOCATE (ij_matrix(Ni_occupied, Nj_occupied))
297 :
298 25 : ij_matrix = 0
299 84 : DO i = 1, Ni_occupied
300 402 : DO j = 1, Nj_occupied
301 318 : counter = counter + 1
302 377 : IF (MOD(counter, para_env_sub%num_pe) == para_env_sub%mepos) THEN
303 318 : ij_matrix(i, j) = ij_matrix(i, j) + 1
304 : END IF
305 : END DO
306 : END DO
307 :
308 75 : ALLOCATE (ij_list_proc_temp(Ni_occupied*occ_j, 2))
309 :
310 25 : elements_ij_proc = 0
311 84 : DO i = 1, Ni_occupied
312 402 : DO j = 1, Nj_occupied
313 318 : IF (ij_matrix(i, j) == 0) CYCLE
314 318 : elements_ij_proc = elements_ij_proc + 1
315 318 : ij_list_proc_temp(elements_ij_proc, 1) = i
316 377 : ij_list_proc_temp(elements_ij_proc, 2) = j
317 : END DO
318 : END DO
319 25 : DEALLOCATE (ij_matrix)
320 :
321 75 : ALLOCATE (ij_list_proc(elements_ij_proc, 2))
322 343 : DO i = 1, elements_ij_proc
323 318 : ij_list_proc(i, 1) = ij_list_proc_temp(i, 1)
324 343 : ij_list_proc(i, 2) = ij_list_proc_temp(i, 2)
325 : END DO
326 25 : DEALLOCATE (ij_list_proc_temp)
327 :
328 25 : IF (.NOT. alpha_beta_case) THEN
329 : CALL mp2_canonical_direct_single_batch(Emp2, Emp2_Cou, Emp2_ex, mp2_env, qs_env, para_env_sub, &
330 : mp2_biel, dimen, C_i, Auto_i, i_batch_start, Ni_occupied, occ_i, &
331 21 : elements_ij_proc, ij_list_proc, Nj_occupied, j_batch_start)
332 : ELSE
333 : CALL mp2_canonical_direct_single_batch(Emp2, Emp2_Cou, Emp2_ex, mp2_env, qs_env, para_env_sub, &
334 : mp2_biel, dimen, C_i, Auto_i, i_batch_start, Ni_occupied, occ_i, &
335 : elements_ij_proc, ij_list_proc, Nj_occupied, j_batch_start, &
336 4 : occ_j, C_j, Auto_j)
337 : END IF
338 :
339 25 : i_batch_start = i_batch_start + Ni_occupied
340 :
341 51 : DEALLOCATE (ij_list_proc)
342 :
343 : END DO
344 :
345 26 : CALL para_env%sum(Emp2_Cou)
346 26 : CALL para_env%sum(Emp2_Ex)
347 26 : CALL para_env%sum(Emp2)
348 :
349 26 : CALL mp_para_env_release(para_env_sub)
350 :
351 26 : CALL timestop(handle)
352 :
353 52 : END SUBROUTINE mp2_direct_energy
354 :
355 : ! **************************************************************************************************
356 : !> \brief ...
357 : !> \param Emp2 ...
358 : !> \param Emp2_Cou ...
359 : !> \param Emp2_ex ...
360 : !> \param mp2_env ...
361 : !> \param qs_env ...
362 : !> \param para_env ...
363 : !> \param mp2_biel ...
364 : !> \param dimen ...
365 : !> \param C ...
366 : !> \param Auto ...
367 : !> \param i_batch_start ...
368 : !> \param Ni_occupied ...
369 : !> \param occupied ...
370 : !> \param elements_ij_proc ...
371 : !> \param ij_list_proc ...
372 : !> \param Nj_occupied ...
373 : !> \param j_batch_start ...
374 : !> \param occupied_beta ...
375 : !> \param C_beta ...
376 : !> \param Auto_beta ...
377 : !> \param Integ_MP2 ...
378 : !> \par History
379 : !> 06.2011 created [Mauro Del Ben]
380 : !> \author Mauro Del Ben
381 : ! **************************************************************************************************
382 43 : SUBROUTINE mp2_canonical_direct_single_batch(Emp2, Emp2_Cou, Emp2_ex, mp2_env, qs_env, para_env, &
383 86 : mp2_biel, dimen, C, Auto, i_batch_start, Ni_occupied, &
384 43 : occupied, elements_ij_proc, ij_list_proc, Nj_occupied, j_batch_start, &
385 10 : occupied_beta, C_beta, Auto_beta, Integ_MP2)
386 :
387 : REAL(KIND=dp), INTENT(INOUT) :: Emp2, Emp2_Cou, Emp2_ex
388 : TYPE(mp2_type) :: mp2_env
389 : TYPE(qs_environment_type), POINTER :: qs_env
390 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
391 : TYPE(mp2_biel_type), INTENT(IN) :: mp2_biel
392 : INTEGER, INTENT(IN) :: dimen
393 : REAL(KIND=dp), DIMENSION(dimen, dimen), INTENT(IN) :: C
394 : REAL(KIND=dp), DIMENSION(dimen), INTENT(IN) :: Auto
395 : INTEGER, INTENT(IN) :: i_batch_start, Ni_occupied, occupied, &
396 : elements_ij_proc
397 : INTEGER, DIMENSION(elements_ij_proc, 2), &
398 : INTENT(IN) :: ij_list_proc
399 : INTEGER, INTENT(IN) :: Nj_occupied, j_batch_start
400 : INTEGER, INTENT(IN), OPTIONAL :: occupied_beta
401 : REAL(KIND=dp), DIMENSION(dimen, dimen), &
402 : INTENT(IN), OPTIONAL :: C_beta
403 : REAL(KIND=dp), DIMENSION(dimen), INTENT(IN), &
404 : OPTIONAL :: Auto_beta
405 : REAL(KIND=dp), ALLOCATABLE, &
406 : DIMENSION(:, :, :, :), INTENT(OUT), OPTIONAL :: Integ_MP2
407 :
408 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mp2_canonical_direct_single_batch'
409 :
410 : INTEGER :: case_index, counter_proc, elements_ij_proc_rec, elements_kl_proc, global_counter, &
411 : handle, i, i_list_ij, i_list_kl, i_set_list_ij, i_set_list_ij_start, i_set_list_ij_stop, &
412 : i_set_list_kl, i_set_list_kl_start, i_set_list_kl_stop, i_start, iatom, iatom_end, &
413 : iatom_start, iiB, ij_elem_max, ikind, index_ij_rec, index_ij_send, index_kl, &
414 : index_proc_ij, index_proc_shift, iset, jatom, jatom_end, jatom_start, jjB, jkind, jset, &
415 : katom, katom_end, katom_start, kkB, kkind, kset, latom, latom_end, latom_start, lkind, &
416 : llB, lset, max_num_call_sec_transf, max_pgf, max_set, multiple
417 : INTEGER :: my_num_call_sec_transf, natom, ncob, nints, nseta, nsetb, nsgf_max, nspins, &
418 : primitive_counter, proc_num, proc_receive, proc_send, R_offset_rec, Rsize_rec, &
419 : S_offset_rec, same_size_kl_index, sphi_a_u1, sphi_a_u2, sphi_a_u3, sphi_b_u1, sphi_b_u2, &
420 : sphi_b_u3, sphi_c_u1, sphi_c_u2, sphi_c_u3, sphi_d_u1, sphi_d_u2, sphi_d_u3, Ssize_rec, &
421 : step_size, total_num_RS_task
422 : INTEGER(int_8) :: estimate_to_store_int, neris_tmp, &
423 : neris_total, nprim_ints
424 43 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of, nimages, proc_num_task, &
425 43 : same_size_kl_elements_counter
426 43 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: kl_list_proc, task_counter_RS, &
427 43 : task_counter_RS_temp
428 : INTEGER, DIMENSION(4) :: RS_counter_temp
429 : INTEGER, DIMENSION(5) :: size_parameter_rec, size_parameter_send
430 43 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, lc_max, &
431 43 : lc_min, ld_max, ld_min, npgfa, npgfb, &
432 43 : npgfc, npgfd, nsgfa, nsgfb, nsgfc, &
433 43 : nsgfd
434 43 : INTEGER, DIMENSION(:, :), POINTER :: nsgfl_a, nsgfl_b, nsgfl_c, nsgfl_d
435 : LOGICAL :: alpha_beta_case, case_send_receive, &
436 : copy_integrals, do_periodic
437 43 : LOGICAL, DIMENSION(:, :), POINTER :: shm_atomic_pair_list
438 : REAL(KIND=dp) :: cartesian_estimate, coeffs_kind_max0, cost_tmp, eps_schwarz, ln_10, &
439 : log10_eps_schwarz, log10_pmax, max_contraction_val, max_val1, max_val2, max_val2_set, &
440 : pmax_atom, pmax_entry, ra(3), rab2, rb(3), rc(3), rcd2, rd(3), screen_kind_ij, &
441 : screen_kind_kl
442 43 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cost_RS, cost_RS_temp, ee_buffer1, &
443 43 : ee_buffer2, ee_primitives_tmp, &
444 43 : ee_work, ee_work2, primitive_integrals
445 43 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: BIb_RS_mat_rec, C_beta_T, max_contraction
446 43 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: BIb, BIb_RS_mat_rec_big, zero_mat_big
447 43 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: BI1, MNRS
448 43 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: shm_pmax_block, zeta, zetb, zetc, zetd
449 43 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: sphi_a_ext_set, sphi_b_ext_set, &
450 43 : sphi_c_ext_set, sphi_d_ext_set
451 43 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: sphi_a_ext, sphi_b_ext, sphi_c_ext, &
452 43 : sphi_d_ext
453 86 : REAL(KIND=dp), DIMENSION(dimen, 2) :: zero_mat
454 86 : REAL(KIND=dp), DIMENSION(dimen, dimen) :: C_T
455 : TYPE(cell_type), POINTER :: cell
456 : TYPE(cp_logger_type), POINTER :: logger
457 : TYPE(four_center_integral_context_type) :: private_context
458 43 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter
459 43 : TYPE(hfx_pgf_list), ALLOCATABLE, DIMENSION(:) :: pgf_list_ij, pgf_list_kl
460 : TYPE(hfx_pgf_product_list), ALLOCATABLE, &
461 43 : DIMENSION(:) :: pgf_product_list
462 : TYPE(hfx_potential_type) :: mp2_potential_parameter
463 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
464 43 : POINTER :: screen_coeffs_kind, tmp_R_1, tmp_R_2, &
465 43 : tmp_screen_pgf1, tmp_screen_pgf2
466 : TYPE(hfx_screen_coeff_type), &
467 43 : DIMENSION(:, :, :, :), POINTER :: screen_coeffs_set
468 : TYPE(hfx_screen_coeff_type), &
469 43 : DIMENSION(:, :, :, :, :, :), POINTER :: radii_pgf, screen_coeffs_pgf
470 : TYPE(hfx_type), POINTER :: actual_x_data
471 43 : TYPE(pair_list_type_mp2) :: list_ij, list_kl
472 : TYPE(pair_set_list_type), ALLOCATABLE, &
473 43 : DIMENSION(:) :: set_list_ij, set_list_kl
474 43 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
475 :
476 43 : CALL timeset(routineN, handle)
477 :
478 : ! The Integ_MP2 will contain the (ia|jb) integrals, necessary for example
479 : ! for the RI-MP2 basis optimization. In this case the number of ij batches
480 : ! has to be equal to 1 (all integrals over molecular orbitals are computed
481 : ! in a single step).
482 43 : copy_integrals = .FALSE.
483 43 : IF (PRESENT(Integ_MP2)) copy_integrals = .TRUE.
484 :
485 43 : alpha_beta_case = .FALSE.
486 :
487 : CALL prepare_integral_calc(cell, qs_env, mp2_env, para_env, mp2_potential_parameter, actual_x_data, &
488 : do_periodic, basis_parameter, max_set, particle_set, natom, kind_of, &
489 : nsgf_max, primitive_integrals, ee_work, ee_work2, ee_buffer1, ee_buffer2, &
490 : ee_primitives_tmp, nspins, max_contraction, max_pgf, pgf_list_ij, &
491 : pgf_list_kl, pgf_product_list, nimages, eps_schwarz, log10_eps_schwarz, &
492 : private_context, screen_coeffs_set, screen_coeffs_kind, screen_coeffs_pgf, &
493 43 : radii_pgf)
494 :
495 43 : ln_10 = LOG(10.0_dp)
496 :
497 43 : neris_tmp = 0_int_8
498 43 : neris_total = 0_int_8
499 43 : nprim_ints = 0_int_8
500 :
501 : !!!!!!!!!
502 880 : ALLOCATE (list_ij%elements(natom**2))
503 837 : ALLOCATE (list_kl%elements(natom**2))
504 : !!!!!!!!!
505 :
506 201 : coeffs_kind_max0 = MAXVAL(screen_coeffs_kind(:, :)%x(2))
507 15190 : ALLOCATE (set_list_ij((max_set*natom)**2))
508 15147 : ALLOCATE (set_list_kl((max_set*natom)**2))
509 :
510 : !! precalculate maximum density matrix elements in blocks
511 365 : actual_x_data%pmax_block = 0.0_dp
512 43 : shm_pmax_block => actual_x_data%pmax_block
513 :
514 43 : shm_atomic_pair_list => actual_x_data%atomic_pair_list
515 :
516 43 : iatom_start = 1
517 43 : iatom_end = natom
518 43 : jatom_start = 1
519 43 : jatom_end = natom
520 43 : katom_start = 1
521 43 : katom_end = natom
522 43 : latom_start = 1
523 43 : latom_end = natom
524 :
525 : CALL build_pair_list_mp2(natom, list_ij, set_list_ij, iatom_start, iatom_end, &
526 : jatom_start, jatom_end, &
527 : kind_of, basis_parameter, particle_set, &
528 : do_periodic, screen_coeffs_set, screen_coeffs_kind, &
529 : coeffs_kind_max0, log10_eps_schwarz, cell, 0.D+00, &
530 43 : shm_atomic_pair_list)
531 :
532 : CALL build_pair_list_mp2(natom, list_kl, set_list_kl, katom_start, katom_end, &
533 : latom_start, latom_end, &
534 : kind_of, basis_parameter, particle_set, &
535 : do_periodic, screen_coeffs_set, screen_coeffs_kind, &
536 : coeffs_kind_max0, log10_eps_schwarz, cell, 0.D+00, &
537 43 : shm_atomic_pair_list)
538 :
539 215 : ALLOCATE (BIb(dimen, dimen, elements_ij_proc))
540 43 : BIb = 0.0D+00
541 64195 : C_T = TRANSPOSE(C)
542 :
543 43 : IF (PRESENT(occupied_beta) .AND. PRESENT(C_beta) .AND. PRESENT(Auto_beta)) THEN
544 40 : ALLOCATE (C_beta_T(dimen, dimen))
545 9322 : C_beta_T(:, :) = TRANSPOSE(C_beta)
546 : alpha_beta_case = .TRUE.
547 : END IF
548 :
549 43 : ij_elem_max = elements_ij_proc
550 43 : CALL para_env%max(ij_elem_max)
551 :
552 : ! calculate the minimum multiple of num_pe >= to Ni_occupied*occupied, in such a way
553 : ! that the i, j loop is performed exactly the same number of time for each procewssor
554 43 : multiple = 0
555 : DO
556 402 : multiple = multiple + para_env%num_pe
557 402 : IF (multiple >= Ni_occupied*Nj_occupied) EXIT
558 : END DO
559 :
560 : ! proc_num_task contains the number of times the second occupied
561 : ! orbital transformation is called for each processor, needs for balancing
562 : ! the point to point send
563 129 : ALLOCATE (proc_num_task(0:para_env%num_pe - 1))
564 :
565 43 : proc_num_task = 0
566 :
567 43 : counter_proc = 0
568 204 : DO i_list_ij = 1, list_ij%n_element
569 161 : iatom = list_ij%elements(i_list_ij)%pair(1)
570 161 : jatom = list_ij%elements(i_list_ij)%pair(2)
571 161 : i_set_list_ij_start = list_ij%elements(i_list_ij)%set_bounds(1)
572 161 : i_set_list_ij_stop = list_ij%elements(i_list_ij)%set_bounds(2)
573 161 : ikind = list_ij%elements(i_list_ij)%kind_pair(1)
574 161 : jkind = list_ij%elements(i_list_ij)%kind_pair(2)
575 :
576 161 : nsgfb => basis_parameter(jkind)%nsgf
577 161 : nsgfa => basis_parameter(ikind)%nsgf
578 :
579 6381 : DO i_set_list_ij = i_set_list_ij_start, i_set_list_ij_stop
580 6177 : iset = set_list_ij(i_set_list_ij)%pair(1)
581 6177 : jset = set_list_ij(i_set_list_ij)%pair(2)
582 6177 : IF (iatom == jatom .AND. jset < iset) CYCLE
583 :
584 4707 : counter_proc = counter_proc + 1
585 4707 : proc_num = MOD(counter_proc, para_env%num_pe)
586 :
587 6338 : proc_num_task(proc_num) = proc_num_task(proc_num) + 1
588 :
589 : END DO
590 : END DO
591 : ! calculate the exact maximum number of calls to the second occupied
592 : ! orbital transformation
593 : ! max_num_call_sec_transf=MAXVAL(proc_num_task)
594 :
595 : ! distribute the RS pair over all processor
596 129 : ALLOCATE (kl_list_proc(proc_num_task(para_env%mepos), 3))
597 :
598 43 : kl_list_proc = 0
599 :
600 43 : counter_proc = 0
601 43 : elements_kl_proc = 0
602 204 : DO i_list_ij = 1, list_ij%n_element
603 161 : iatom = list_ij%elements(i_list_ij)%pair(1)
604 161 : jatom = list_ij%elements(i_list_ij)%pair(2)
605 161 : i_set_list_ij_start = list_ij%elements(i_list_ij)%set_bounds(1)
606 161 : i_set_list_ij_stop = list_ij%elements(i_list_ij)%set_bounds(2)
607 161 : ikind = list_ij%elements(i_list_ij)%kind_pair(1)
608 161 : jkind = list_ij%elements(i_list_ij)%kind_pair(2)
609 :
610 161 : nsgfb => basis_parameter(jkind)%nsgf
611 161 : nsgfa => basis_parameter(ikind)%nsgf
612 :
613 6381 : DO i_set_list_ij = i_set_list_ij_start, i_set_list_ij_stop
614 6177 : iset = set_list_ij(i_set_list_ij)%pair(1)
615 6177 : jset = set_list_ij(i_set_list_ij)%pair(2)
616 6177 : IF (iatom == jatom .AND. jset < iset) CYCLE
617 :
618 4707 : counter_proc = counter_proc + 1
619 4707 : proc_num = MOD(counter_proc, para_env%num_pe)
620 :
621 4868 : IF (proc_num == para_env%mepos) THEN
622 4653 : elements_kl_proc = elements_kl_proc + 1
623 4653 : kl_list_proc(elements_kl_proc, 1) = i_list_ij
624 4653 : kl_list_proc(elements_kl_proc, 2) = i_set_list_ij
625 4653 : kl_list_proc(elements_kl_proc, 3) = counter_proc
626 : END IF
627 :
628 : END DO
629 : END DO
630 :
631 104 : total_num_RS_task = SUM(proc_num_task)
632 129 : ALLOCATE (task_counter_RS(total_num_RS_task, 4))
633 :
634 129 : ALLOCATE (cost_RS(total_num_RS_task))
635 :
636 43 : task_counter_RS = 0
637 43 : cost_RS = 0.0_dp
638 :
639 129 : DO case_index = 1, 2
640 :
641 : my_num_call_sec_transf = 0
642 9392 : DO index_kl = 1, elements_kl_proc
643 :
644 9306 : i_list_ij = kl_list_proc(index_kl, 1)
645 9306 : i_set_list_ij = kl_list_proc(index_kl, 2)
646 :
647 9306 : iatom = list_ij%elements(i_list_ij)%pair(1)
648 9306 : jatom = list_ij%elements(i_list_ij)%pair(2)
649 9306 : i_set_list_ij_start = list_ij%elements(i_list_ij)%set_bounds(1)
650 9306 : i_set_list_ij_stop = list_ij%elements(i_list_ij)%set_bounds(2)
651 9306 : ikind = list_ij%elements(i_list_ij)%kind_pair(1)
652 9306 : jkind = list_ij%elements(i_list_ij)%kind_pair(2)
653 37224 : ra = list_ij%elements(i_list_ij)%r1
654 37224 : rb = list_ij%elements(i_list_ij)%r2
655 9306 : rab2 = list_ij%elements(i_list_ij)%dist2
656 :
657 9306 : la_max => basis_parameter(ikind)%lmax
658 9306 : la_min => basis_parameter(ikind)%lmin
659 9306 : npgfa => basis_parameter(ikind)%npgf
660 9306 : nseta = basis_parameter(ikind)%nset
661 9306 : zeta => basis_parameter(ikind)%zet
662 9306 : nsgfa => basis_parameter(ikind)%nsgf
663 9306 : sphi_a_ext => basis_parameter(ikind)%sphi_ext(:, :, :, :)
664 9306 : nsgfl_a => basis_parameter(ikind)%nsgfl
665 9306 : sphi_a_u1 = UBOUND(sphi_a_ext, 1)
666 9306 : sphi_a_u2 = UBOUND(sphi_a_ext, 2)
667 9306 : sphi_a_u3 = UBOUND(sphi_a_ext, 3)
668 :
669 9306 : lb_max => basis_parameter(jkind)%lmax
670 9306 : lb_min => basis_parameter(jkind)%lmin
671 9306 : npgfb => basis_parameter(jkind)%npgf
672 9306 : nsetb = basis_parameter(jkind)%nset
673 9306 : zetb => basis_parameter(jkind)%zet
674 9306 : nsgfb => basis_parameter(jkind)%nsgf
675 9306 : sphi_b_ext => basis_parameter(jkind)%sphi_ext(:, :, :, :)
676 9306 : nsgfl_b => basis_parameter(jkind)%nsgfl
677 9306 : sphi_b_u1 = UBOUND(sphi_b_ext, 1)
678 9306 : sphi_b_u2 = UBOUND(sphi_b_ext, 2)
679 9306 : sphi_b_u3 = UBOUND(sphi_b_ext, 3)
680 :
681 9306 : iset = set_list_ij(i_set_list_ij)%pair(1)
682 9306 : jset = set_list_ij(i_set_list_ij)%pair(2)
683 :
684 9306 : ncob = npgfb(jset)*ncoset(lb_max(jset))
685 : max_val1 = screen_coeffs_set(jset, iset, jkind, ikind)%x(1)*rab2 + &
686 9306 : screen_coeffs_set(jset, iset, jkind, ikind)%x(2)
687 :
688 9306 : sphi_a_ext_set => sphi_a_ext(:, :, :, iset)
689 9306 : sphi_b_ext_set => sphi_b_ext(:, :, :, jset)
690 :
691 9306 : IF (case_index == 1) THEN
692 4653 : global_counter = kl_list_proc(index_kl, 3)
693 4653 : task_counter_RS(global_counter, 1) = i_list_ij
694 4653 : task_counter_RS(global_counter, 2) = i_set_list_ij
695 4653 : task_counter_RS(global_counter, 3) = nsgfb(jset)*nsgfa(iset)
696 : END IF
697 :
698 9306 : IF (ALLOCATED(BI1)) DEALLOCATE (BI1)
699 55836 : ALLOCATE (BI1(dimen, Ni_occupied, nsgfb(jset), nsgfa(iset)))
700 :
701 9306 : BI1 = 0.D+00
702 :
703 71604 : DO i_list_kl = 1, list_kl%n_element
704 :
705 62298 : katom = list_kl%elements(i_list_kl)%pair(1)
706 62298 : latom = list_kl%elements(i_list_kl)%pair(2)
707 :
708 62298 : i_set_list_kl_start = list_kl%elements(i_list_kl)%set_bounds(1)
709 62298 : i_set_list_kl_stop = list_kl%elements(i_list_kl)%set_bounds(2)
710 62298 : kkind = list_kl%elements(i_list_kl)%kind_pair(1)
711 62298 : lkind = list_kl%elements(i_list_kl)%kind_pair(2)
712 249192 : rc = list_kl%elements(i_list_kl)%r1
713 249192 : rd = list_kl%elements(i_list_kl)%r2
714 62298 : rcd2 = list_kl%elements(i_list_kl)%dist2
715 :
716 62298 : pmax_atom = 0.0_dp
717 :
718 : screen_kind_ij = screen_coeffs_kind(jkind, ikind)%x(1)*rab2 + &
719 62298 : screen_coeffs_kind(jkind, ikind)%x(2)
720 : screen_kind_kl = screen_coeffs_kind(lkind, kkind)%x(1)*rcd2 + &
721 62298 : screen_coeffs_kind(lkind, kkind)%x(2)
722 :
723 : !!!!! Change the loop order
724 62298 : IF (max_val1 + screen_kind_kl + pmax_atom < log10_eps_schwarz) CYCLE
725 : !!!!!
726 62274 : IF (screen_kind_ij + screen_kind_kl + pmax_atom < log10_eps_schwarz) CYCLE
727 :
728 62274 : lc_max => basis_parameter(kkind)%lmax
729 62274 : lc_min => basis_parameter(kkind)%lmin
730 62274 : npgfc => basis_parameter(kkind)%npgf
731 62274 : zetc => basis_parameter(kkind)%zet
732 62274 : nsgfc => basis_parameter(kkind)%nsgf
733 62274 : sphi_c_ext => basis_parameter(kkind)%sphi_ext(:, :, :, :)
734 62274 : nsgfl_c => basis_parameter(kkind)%nsgfl
735 62274 : sphi_c_u1 = UBOUND(sphi_c_ext, 1)
736 62274 : sphi_c_u2 = UBOUND(sphi_c_ext, 2)
737 62274 : sphi_c_u3 = UBOUND(sphi_c_ext, 3)
738 :
739 62274 : ld_max => basis_parameter(lkind)%lmax
740 62274 : ld_min => basis_parameter(lkind)%lmin
741 62274 : npgfd => basis_parameter(lkind)%npgf
742 62274 : zetd => basis_parameter(lkind)%zet
743 62274 : nsgfd => basis_parameter(lkind)%nsgf
744 62274 : sphi_d_ext => basis_parameter(lkind)%sphi_ext(:, :, :, :)
745 62274 : nsgfl_d => basis_parameter(lkind)%nsgfl
746 62274 : sphi_d_u1 = UBOUND(sphi_d_ext, 1)
747 62274 : sphi_d_u2 = UBOUND(sphi_d_ext, 2)
748 62274 : sphi_d_u3 = UBOUND(sphi_d_ext, 3)
749 :
750 2761710 : DO i_set_list_kl = i_set_list_kl_start, i_set_list_kl_stop
751 2690130 : kset = set_list_kl(i_set_list_kl)%pair(1)
752 2690130 : lset = set_list_kl(i_set_list_kl)%pair(2)
753 :
754 2690130 : IF (katom == latom .AND. lset < kset) CYCLE
755 :
756 : max_val2_set = (screen_coeffs_set(lset, kset, lkind, kkind)%x(1)*rcd2 + &
757 2076726 : screen_coeffs_set(lset, kset, lkind, kkind)%x(2))
758 2076726 : max_val2 = max_val1 + max_val2_set
759 :
760 : !! Near field screening
761 2076726 : IF (max_val2 + pmax_atom < log10_eps_schwarz) CYCLE
762 2023634 : sphi_c_ext_set => sphi_c_ext(:, :, :, kset)
763 2023634 : sphi_d_ext_set => sphi_d_ext(:, :, :, lset)
764 : !! get max_vals if we screen on initial density
765 2023634 : pmax_entry = 0.0_dp
766 :
767 2023634 : log10_pmax = pmax_entry
768 2023634 : max_val2 = max_val2 + log10_pmax
769 : IF (max_val2 < log10_eps_schwarz) CYCLE
770 2023634 : pmax_entry = EXP(log10_pmax*ln_10)
771 :
772 2085932 : IF (case_index == 2) THEN
773 1011817 : IF (ALLOCATED(MNRS)) DEALLOCATE (MNRS)
774 6070902 : ALLOCATE (MNRS(nsgfd(lset), nsgfc(kset), nsgfb(jset), nsgfa(iset)))
775 :
776 1011817 : MNRS = 0.D+00
777 :
778 : max_contraction_val = max_contraction(iset, iatom)* &
779 : max_contraction(jset, jatom)* &
780 : max_contraction(kset, katom)* &
781 1011817 : max_contraction(lset, latom)*pmax_entry
782 1011817 : tmp_R_1 => radii_pgf(:, :, jset, iset, jkind, ikind)
783 1011817 : tmp_R_2 => radii_pgf(:, :, lset, kset, lkind, kkind)
784 1011817 : tmp_screen_pgf1 => screen_coeffs_pgf(:, :, jset, iset, jkind, ikind)
785 1011817 : tmp_screen_pgf2 => screen_coeffs_pgf(:, :, lset, kset, lkind, kkind)
786 :
787 : CALL coulomb4(private_context, ra, rb, rc, rd, npgfa(iset), npgfb(jset), npgfc(kset), npgfd(lset), &
788 : la_min(iset), la_max(iset), lb_min(jset), lb_max(jset), &
789 : lc_min(kset), lc_max(kset), ld_min(lset), ld_max(lset), &
790 : nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
791 : sphi_a_u1, sphi_a_u2, sphi_a_u3, &
792 : sphi_b_u1, sphi_b_u2, sphi_b_u3, &
793 : sphi_c_u1, sphi_c_u2, sphi_c_u3, &
794 : sphi_d_u1, sphi_d_u2, sphi_d_u3, &
795 : zeta(1:npgfa(iset), iset), zetb(1:npgfb(jset), jset), &
796 : zetc(1:npgfc(kset), kset), zetd(1:npgfd(lset), lset), &
797 : primitive_integrals, &
798 : mp2_potential_parameter, &
799 : actual_x_data%neighbor_cells, screen_coeffs_set(jset, iset, jkind, ikind)%x, &
800 : screen_coeffs_set(lset, kset, lkind, kkind)%x, eps_schwarz, &
801 : max_contraction_val, cartesian_estimate, cell, neris_tmp, &
802 : log10_pmax, log10_eps_schwarz, &
803 : tmp_R_1, tmp_R_2, tmp_screen_pgf1, tmp_screen_pgf2, &
804 : pgf_list_ij, pgf_list_kl, pgf_product_list, &
805 : nsgfl_a(:, iset), nsgfl_b(:, jset), &
806 : nsgfl_c(:, kset), nsgfl_d(:, lset), &
807 : sphi_a_ext_set, &
808 : sphi_b_ext_set, &
809 : sphi_c_ext_set, &
810 : sphi_d_ext_set, &
811 : ee_work, ee_work2, ee_buffer1, ee_buffer2, ee_primitives_tmp, &
812 1011817 : nimages, do_periodic)
813 :
814 1011817 : nints = nsgfa(iset)*nsgfb(jset)*nsgfc(kset)*nsgfd(lset)
815 1011817 : neris_total = neris_total + nints
816 1011817 : nprim_ints = nprim_ints + neris_tmp
817 1011817 : IF (cartesian_estimate == 0.0_dp) cartesian_estimate = TINY(cartesian_estimate)
818 1011817 : estimate_to_store_int = EXPONENT(cartesian_estimate)
819 1011817 : estimate_to_store_int = MAX(estimate_to_store_int, -15_int_8)
820 1011817 : cartesian_estimate = SET_EXPONENT(1.0_dp, estimate_to_store_int + 1)
821 :
822 1011817 : IF (cartesian_estimate < eps_schwarz) CYCLE
823 :
824 : primitive_counter = 0
825 3502120 : DO llB = 1, nsgfd(lset)
826 10172593 : DO kkB = 1, nsgfc(kset)
827 28163930 : DO jjB = 1, nsgfb(jset)
828 80468051 : DO iiB = 1, nsgfa(iset)
829 54796116 : primitive_counter = primitive_counter + 1
830 73797578 : MNRS(llB, kkB, jjB, iiB) = primitive_integrals(primitive_counter)
831 : END DO
832 : END DO
833 : END DO
834 : END DO
835 :
836 : CALL transform_occupied_orbitals_first(dimen, iatom, jatom, katom, latom, &
837 : iset, jset, kset, lset, &
838 : nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
839 : i_batch_start, Ni_occupied, &
840 1010125 : MNRS, C_T, mp2_biel, BI1)
841 : ELSE
842 1011817 : task_counter_RS(global_counter, 4) = task_counter_RS(global_counter, 4) + 1
843 :
844 1011817 : cost_tmp = 0.0_dp
845 : cost_tmp = cost_model(nsgfd(lset), nsgfc(kset), nsgfb(jset), nsgfa(iset), &
846 : npgfd(lset), npgfc(kset), npgfb(jset), npgfa(iset), &
847 : max_val2/log10_eps_schwarz, &
848 1011817 : p1_energy, p2_energy, p3_energy)
849 1011817 : cost_RS(global_counter) = cost_RS(global_counter) + cost_tmp
850 : END IF
851 :
852 : END DO ! i_set_list_kl
853 : END DO ! i_list_kl
854 :
855 9392 : IF (case_index == 2) THEN
856 4653 : my_num_call_sec_transf = my_num_call_sec_transf + 1
857 4653 : IF (.NOT. alpha_beta_case) THEN
858 3993 : IF (.NOT. mp2_env%direct_canonical%big_send) THEN
859 : CALL transform_occupied_orbitals_second(dimen, iatom, jatom, iset, jset, &
860 : nsgfa(iset), nsgfb(jset), Ni_occupied, Nj_occupied, j_batch_start, &
861 : BI1, C_T, mp2_biel, para_env, elements_ij_proc, &
862 0 : multiple, BIb)
863 : ELSE
864 : CALL transform_occupied_orbitals_second_big( &
865 : dimen, iatom, jatom, iset, jset, &
866 : nsgfa(iset), nsgfb(jset), Ni_occupied, Nj_occupied, j_batch_start, &
867 3993 : ij_elem_max, BI1, C_T, mp2_biel, para_env, elements_ij_proc, BIb)
868 : END IF
869 : ELSE
870 660 : IF (.NOT. mp2_env%direct_canonical%big_send) THEN
871 : CALL transform_occupied_orbitals_second(dimen, iatom, jatom, iset, jset, &
872 : nsgfa(iset), nsgfb(jset), Ni_occupied, Nj_occupied, j_batch_start, &
873 : BI1, C_beta_T, mp2_biel, para_env, elements_ij_proc, &
874 0 : multiple, BIb)
875 : ELSE
876 : CALL transform_occupied_orbitals_second_big( &
877 : dimen, iatom, jatom, iset, jset, &
878 : nsgfa(iset), nsgfb(jset), Ni_occupied, Nj_occupied, j_batch_start, &
879 660 : ij_elem_max, BI1, C_beta_T, mp2_biel, para_env, elements_ij_proc, BIb)
880 : END IF
881 : END IF
882 : END IF
883 :
884 : END DO !i_list_ij
885 :
886 129 : IF (case_index == 1) THEN
887 43 : CALL para_env%sum(task_counter_RS)
888 43 : CALL para_env%sum(cost_RS)
889 86 : ALLOCATE (task_counter_RS_temp(total_num_RS_task, 4))
890 :
891 86 : ALLOCATE (cost_RS_temp(total_num_RS_task))
892 :
893 43 : step_size = 1
894 129 : ALLOCATE (same_size_kl_elements_counter((nsgf_max**2 + 1)/step_size + 1))
895 :
896 43 : same_size_kl_elements_counter = 0
897 :
898 43 : same_size_kl_index = 0
899 43 : global_counter = 0
900 43 : DO iiB = nsgf_max**2 + 1, 0, -step_size
901 149802 : DO jjB = 1, total_num_RS_task
902 149802 : IF (task_counter_RS(jjB, 3) > iiB - step_size .AND. task_counter_RS(jjB, 3) <= iiB) THEN
903 4707 : global_counter = global_counter + 1
904 23535 : task_counter_RS_temp(global_counter, 1:4) = task_counter_RS(jjB, 1:4)
905 4707 : cost_RS_temp(global_counter) = cost_RS(jjB)
906 : END IF
907 : END DO
908 1329 : same_size_kl_index = same_size_kl_index + 1
909 1329 : same_size_kl_elements_counter(same_size_kl_index) = global_counter
910 : END DO
911 :
912 43 : DEALLOCATE (task_counter_RS)
913 43 : DEALLOCATE (cost_RS)
914 :
915 43 : i_start = 1
916 1372 : DO same_size_kl_index = 1, SIZE(same_size_kl_elements_counter)
917 6036 : DO iiB = i_start, same_size_kl_elements_counter(same_size_kl_index)
918 100828 : DO jjB = iiB + 1, same_size_kl_elements_counter(same_size_kl_index)
919 :
920 99499 : IF (cost_RS_temp(jjB) >= cost_RS_temp(iiB)) THEN
921 214310 : RS_counter_temp = task_counter_RS_temp(iiB, 1:4)
922 214310 : task_counter_RS_temp(iiB, 1:4) = task_counter_RS_temp(jjB, 1:4)
923 214310 : task_counter_RS_temp(jjB, 1:4) = RS_counter_temp
924 :
925 42862 : cost_tmp = cost_RS_temp(iiB)
926 42862 : cost_RS_temp(iiB) = cost_RS_temp(jjB)
927 42862 : cost_RS_temp(jjB) = cost_tmp
928 : END IF
929 : END DO
930 : END DO
931 1372 : i_start = same_size_kl_elements_counter(same_size_kl_index) + 1
932 : END DO
933 :
934 43 : proc_num_task = 0
935 4750 : DO counter_proc = 1, total_num_RS_task
936 4707 : proc_num = MOD(counter_proc, para_env%num_pe)
937 4750 : proc_num_task(proc_num) = proc_num_task(proc_num) + 1
938 : END DO
939 :
940 104 : max_num_call_sec_transf = MAXVAL(proc_num_task)
941 :
942 43 : DEALLOCATE (kl_list_proc)
943 129 : ALLOCATE (kl_list_proc(proc_num_task(para_env%mepos), 2))
944 :
945 43 : kl_list_proc = 0
946 :
947 43 : elements_kl_proc = 0
948 4750 : DO counter_proc = 1, total_num_RS_task
949 4707 : proc_num = MOD(counter_proc, para_env%num_pe)
950 4750 : IF (proc_num == para_env%mepos) THEN
951 4653 : elements_kl_proc = elements_kl_proc + 1
952 4653 : kl_list_proc(elements_kl_proc, 1) = task_counter_RS_temp(counter_proc, 1)
953 4653 : kl_list_proc(elements_kl_proc, 2) = task_counter_RS_temp(counter_proc, 2)
954 : END IF
955 : END DO
956 :
957 43 : DEALLOCATE (task_counter_RS_temp)
958 43 : DEALLOCATE (cost_RS_temp)
959 : END IF
960 : END DO ! case_index
961 :
962 43 : size_parameter_send(1) = 1
963 43 : size_parameter_send(2) = 1
964 43 : size_parameter_send(3) = 0
965 43 : size_parameter_send(4) = 0
966 43 : size_parameter_send(5) = elements_ij_proc
967 :
968 43 : IF (mp2_env%direct_canonical%big_send) THEN
969 172 : ALLOCATE (zero_mat_big(dimen, 2, ij_elem_max))
970 :
971 : END IF
972 :
973 43 : DO iiB = my_num_call_sec_transf + 1, max_num_call_sec_transf
974 43 : DO index_proc_shift = 0, para_env%num_pe - 1
975 :
976 0 : proc_send = MODULO(para_env%mepos + index_proc_shift, para_env%num_pe)
977 0 : proc_receive = MODULO(para_env%mepos - index_proc_shift, para_env%num_pe)
978 :
979 0 : case_send_receive = (proc_send /= para_env%mepos)
980 :
981 0 : IF (case_send_receive) THEN
982 : ! the processor starts to send (and receive?)
983 :
984 0 : CALL para_env%sendrecv(size_parameter_send, proc_send, size_parameter_rec, proc_receive)
985 :
986 0 : Rsize_rec = size_parameter_rec(1)
987 0 : Ssize_rec = size_parameter_rec(2)
988 0 : R_offset_rec = size_parameter_rec(3)
989 0 : S_offset_rec = size_parameter_rec(4)
990 0 : elements_ij_proc_rec = size_parameter_rec(5)
991 0 : IF (.NOT. mp2_env%direct_canonical%big_send) THEN
992 0 : ALLOCATE (BIb_RS_mat_rec(dimen, Rsize_rec + Ssize_rec))
993 : ELSE
994 0 : ALLOCATE (BIb_RS_mat_rec_big(dimen, Rsize_rec + Ssize_rec, ij_elem_max))
995 : END IF
996 : ELSE
997 0 : elements_ij_proc_rec = elements_ij_proc
998 : END IF
999 :
1000 0 : IF (.NOT. mp2_env%direct_canonical%big_send) THEN
1001 0 : index_ij_send = 0
1002 0 : index_ij_rec = 0
1003 0 : DO index_proc_ij = proc_send + 1, multiple, para_env%num_pe
1004 :
1005 0 : zero_mat = 0.D+00
1006 :
1007 0 : IF (case_send_receive) THEN
1008 :
1009 0 : CALL para_env%sendrecv(zero_mat, proc_send, BIb_RS_mat_rec, proc_receive)
1010 :
1011 0 : index_ij_rec = index_ij_rec + 1
1012 0 : IF (index_ij_rec <= elements_ij_proc .AND. elements_ij_proc > 0) THEN
1013 :
1014 : BIb(1:dimen, R_offset_rec + 1:R_offset_rec + Rsize_rec, index_ij_rec) = &
1015 : BIb(1:dimen, R_offset_rec + 1:R_offset_rec + Rsize_rec, index_ij_rec) + &
1016 0 : BIb_RS_mat_rec(1:dimen, 1:Rsize_rec)
1017 :
1018 : BIb(1:dimen, S_offset_rec + 1:S_offset_rec + Ssize_rec, index_ij_rec) = &
1019 : BIb(1:dimen, S_offset_rec + 1:S_offset_rec + Ssize_rec, index_ij_rec) + &
1020 0 : BIb_RS_mat_rec(1:dimen, Rsize_rec + 1:Rsize_rec + Ssize_rec)
1021 :
1022 : END IF
1023 : END IF
1024 :
1025 : END DO
1026 : ELSE
1027 0 : zero_mat_big = 0.D+00
1028 :
1029 0 : IF (case_send_receive) THEN
1030 :
1031 0 : CALL para_env%sendrecv(zero_mat_big, proc_send, BIb_RS_mat_rec_big, proc_receive)
1032 :
1033 : BIb(1:dimen, R_offset_rec + 1:R_offset_rec + Rsize_rec, 1:elements_ij_proc) = &
1034 : BIb(1:dimen, R_offset_rec + 1:R_offset_rec + Rsize_rec, 1:elements_ij_proc) + &
1035 0 : BIb_RS_mat_rec_big(1:dimen, 1:Rsize_rec, 1:elements_ij_proc)
1036 :
1037 : BIb(1:dimen, S_offset_rec + 1:S_offset_rec + Ssize_rec, 1:elements_ij_proc) = &
1038 : BIb(1:dimen, S_offset_rec + 1:S_offset_rec + Ssize_rec, 1:elements_ij_proc) + &
1039 0 : BIb_RS_mat_rec_big(1:dimen, Rsize_rec + 1:Rsize_rec + Ssize_rec, 1:elements_ij_proc)
1040 :
1041 : END IF
1042 : END IF
1043 :
1044 0 : IF (case_send_receive) THEN
1045 0 : IF (.NOT. mp2_env%direct_canonical%big_send) THEN
1046 0 : DEALLOCATE (BIb_RS_mat_rec)
1047 : ELSE
1048 0 : DEALLOCATE (BIb_RS_mat_rec_big)
1049 : END IF
1050 : END IF
1051 :
1052 : END DO
1053 : END DO
1054 :
1055 43 : IF (mp2_env%direct_canonical%big_send) THEN
1056 43 : DEALLOCATE (zero_mat_big)
1057 : END IF
1058 :
1059 43 : logger => cp_get_default_logger()
1060 :
1061 43 : DEALLOCATE (primitive_integrals)
1062 :
1063 43 : IF (.NOT. alpha_beta_case) THEN
1064 : CALL transform_virtual_orbitals_and_accumulate(dimen, occupied, dimen - occupied, i_batch_start, &
1065 : j_batch_start, BIb, C, Auto, elements_ij_proc, ij_list_proc, &
1066 33 : nspins, Emp2, Emp2_Cou, Emp2_ex)
1067 : ELSE
1068 : CALL transform_virtual_orbitals_and_accumulate_ABcase( &
1069 : dimen, occupied, occupied_beta, dimen - occupied, dimen - occupied_beta, &
1070 : i_batch_start, j_batch_start, &
1071 : BIb, C, C_beta, Auto, Auto_beta, &
1072 10 : elements_ij_proc, ij_list_proc, Emp2, Emp2_Cou)
1073 10 : DEALLOCATE (C_beta_T)
1074 : END IF
1075 :
1076 43 : IF (copy_integrals) THEN
1077 18 : IF (.NOT. alpha_beta_case) THEN
1078 72 : ALLOCATE (Integ_MP2(dimen - occupied, dimen - occupied, occupied, occupied))
1079 12 : Integ_MP2 = 0.0_dp
1080 72 : DO i = 1, elements_ij_proc
1081 60 : iiB = ij_list_proc(i, 1)
1082 60 : jjB = ij_list_proc(i, 2)
1083 5976 : Integ_MP2(:, :, iiB + i_batch_start, jjB + j_batch_start) = BIb(1:dimen - occupied, 1:dimen - occupied, i)
1084 : END DO
1085 : ELSE
1086 36 : ALLOCATE (Integ_MP2(dimen - occupied, dimen - occupied_beta, occupied, occupied_beta))
1087 6 : Integ_MP2 = 0.0_dp
1088 30 : DO i = 1, elements_ij_proc
1089 24 : iiB = ij_list_proc(i, 1)
1090 24 : jjB = ij_list_proc(i, 2)
1091 2670 : Integ_MP2(:, :, iiB + i_batch_start, jjB + j_batch_start) = BIb(1:dimen - occupied, 1:dimen - occupied_beta, i)
1092 : END DO
1093 : END IF
1094 : END IF
1095 43 : DEALLOCATE (BIb)
1096 :
1097 43 : DEALLOCATE (set_list_ij, set_list_kl)
1098 :
1099 1665 : DO i = 1, max_pgf**2
1100 1622 : DEALLOCATE (pgf_list_ij(i)%image_list)
1101 1665 : DEALLOCATE (pgf_list_kl(i)%image_list)
1102 : END DO
1103 :
1104 43 : DEALLOCATE (pgf_list_ij)
1105 43 : DEALLOCATE (pgf_list_kl)
1106 43 : DEALLOCATE (pgf_product_list)
1107 :
1108 43 : DEALLOCATE (max_contraction, kind_of)
1109 :
1110 43 : DEALLOCATE (ee_work, ee_work2, ee_buffer1, ee_buffer2, ee_primitives_tmp)
1111 :
1112 43 : DEALLOCATE (nimages)
1113 :
1114 43 : IF (mp2_env%potential_parameter%potential_type == do_potential_TShPSC) THEN
1115 2 : init_TShPSC_lmax = -1
1116 2 : CALL free_C0()
1117 : END IF
1118 :
1119 43 : CALL timestop(handle)
1120 :
1121 96 : END SUBROUTINE mp2_canonical_direct_single_batch
1122 :
1123 : ! **************************************************************************************************
1124 : !> \brief ...
1125 : !> \param dimen ...
1126 : !> \param latom ...
1127 : !> \param katom ...
1128 : !> \param jatom ...
1129 : !> \param iatom ...
1130 : !> \param lset ...
1131 : !> \param kset ...
1132 : !> \param jset ...
1133 : !> \param iset ...
1134 : !> \param Ssize ...
1135 : !> \param Rsize ...
1136 : !> \param Nsize ...
1137 : !> \param Msize ...
1138 : !> \param i_batch_start ...
1139 : !> \param Ni_occupied ...
1140 : !> \param MNRS ...
1141 : !> \param C_T ...
1142 : !> \param mp2_biel ...
1143 : !> \param BI1 ...
1144 : ! **************************************************************************************************
1145 1010125 : SUBROUTINE transform_occupied_orbitals_first(dimen, latom, katom, jatom, iatom, &
1146 : lset, kset, jset, iset, &
1147 : Ssize, Rsize, Nsize, Msize, &
1148 : i_batch_start, Ni_occupied, &
1149 1010125 : MNRS, C_T, mp2_biel, BI1)
1150 :
1151 : INTEGER, INTENT(IN) :: dimen, latom, katom, jatom, iatom, lset, &
1152 : kset, jset, iset, Ssize, Rsize, Nsize, &
1153 : Msize, i_batch_start, Ni_occupied
1154 : REAL(KIND=dp), &
1155 : DIMENSION(Msize, Nsize, Rsize, Ssize), &
1156 : INTENT(IN) :: MNRS
1157 : REAL(KIND=dp), DIMENSION(dimen, dimen), INTENT(IN) :: C_T
1158 : TYPE(mp2_biel_type), INTENT(IN) :: mp2_biel
1159 : REAL(KIND=dp), &
1160 : DIMENSION(dimen, Ni_occupied, Rsize, Ssize), &
1161 : INTENT(INOUT) :: BI1
1162 :
1163 : INTEGER :: i, i_global, m, M_global, M_offset, &
1164 : M_start, n, N_global, N_offset, r, &
1165 : R_offset, R_start, s, S_offset
1166 : REAL(KIND=dp) :: MNRS_element
1167 :
1168 1010125 : N_offset = mp2_biel%index_table(jatom, jset) - 1
1169 1010125 : M_offset = mp2_biel%index_table(iatom, iset) - 1
1170 1010125 : S_offset = mp2_biel%index_table(latom, lset) - 1
1171 1010125 : R_offset = mp2_biel%index_table(katom, kset) - 1
1172 :
1173 3549662 : DO S = 1, Ssize
1174 2539537 : R_start = 1
1175 2539537 : IF (katom == latom .AND. kset == lset) R_start = S
1176 9925984 : DO R = R_start, Rsize
1177 :
1178 : ! fast i don't know why
1179 25959912 : DO N = 1, Nsize
1180 17044053 : N_global = N + N_offset
1181 17044053 : M_start = 1
1182 17044053 : IF (iatom == jatom .AND. iset == jset) THEN
1183 1686383 : M = N
1184 1686383 : M_global = M + M_offset
1185 1686383 : MNRS_element = MNRS(M, N, R, S)
1186 7429233 : DO i = 1, Ni_occupied
1187 5742850 : i_global = i + i_batch_start
1188 7429233 : BI1(N_global, i, R, S) = BI1(N_global, i, R, S) + C_T(i_global, M_global)*MNRS_element
1189 : END DO
1190 1686383 : M_start = N + 1
1191 : END IF
1192 :
1193 71498951 : DO M = M_start, Msize
1194 48078576 : M_global = M + M_offset
1195 48078576 : MNRS_element = MNRS(M, N, R, S)
1196 247861449 : DO i = 1, Ni_occupied
1197 182738820 : i_global = i + i_batch_start
1198 182738820 : BI1(N_global, i, R, S) = BI1(N_global, i, R, S) + C_T(i_global, M_global)*MNRS_element
1199 230817396 : BI1(M_global, i, R, S) = BI1(M_global, i, R, S) + C_T(i_global, N_global)*MNRS_element
1200 : END DO
1201 : END DO
1202 : END DO
1203 :
1204 : END DO
1205 : END DO
1206 :
1207 1010125 : END SUBROUTINE transform_occupied_orbitals_first
1208 :
1209 : ! **************************************************************************************************
1210 : !> \brief ...
1211 : !> \param dimen ...
1212 : !> \param latom ...
1213 : !> \param katom ...
1214 : !> \param lset ...
1215 : !> \param kset ...
1216 : !> \param Ssize ...
1217 : !> \param Rsize ...
1218 : !> \param Ni_occupied ...
1219 : !> \param Nj_occupied ...
1220 : !> \param j_batch_start ...
1221 : !> \param BI1 ...
1222 : !> \param C_T ...
1223 : !> \param mp2_biel ...
1224 : !> \param para_env ...
1225 : !> \param elements_ij_proc ...
1226 : !> \param multiple ...
1227 : !> \param BIb ...
1228 : ! **************************************************************************************************
1229 0 : SUBROUTINE transform_occupied_orbitals_second(dimen, latom, katom, lset, kset, &
1230 : Ssize, Rsize, Ni_occupied, Nj_occupied, j_batch_start, &
1231 0 : BI1, C_T, mp2_biel, para_env, &
1232 : elements_ij_proc, &
1233 0 : multiple, BIb)
1234 :
1235 : INTEGER, INTENT(IN) :: dimen, latom, katom, lset, kset, Ssize, &
1236 : Rsize, Ni_occupied, Nj_occupied, &
1237 : j_batch_start
1238 : REAL(KIND=dp), &
1239 : DIMENSION(dimen, Ni_occupied, Rsize, Ssize), &
1240 : INTENT(IN) :: BI1
1241 : REAL(KIND=dp), DIMENSION(dimen, dimen), INTENT(IN) :: C_T
1242 : TYPE(mp2_biel_type), INTENT(IN) :: mp2_biel
1243 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
1244 : INTEGER, INTENT(IN) :: elements_ij_proc, multiple
1245 : REAL(KIND=dp), &
1246 : DIMENSION(dimen, dimen, elements_ij_proc), &
1247 : INTENT(INOUT) :: BIb
1248 :
1249 : CHARACTER(LEN=*), PARAMETER :: routineN = 'transform_occupied_orbitals_second'
1250 :
1251 : INTEGER :: elements_ij_proc_rec, handle, i, index_ij_rec, index_ij_send, index_proc_ij, &
1252 : index_proc_shift, j, n, proc_receive, proc_send, r, R_global, R_offset, R_offset_rec, &
1253 : R_start, Rsize_rec, s, S_global, S_offset, S_offset_rec, Ssize_rec
1254 : INTEGER, DIMENSION(5) :: size_parameter_rec, size_parameter_send
1255 : LOGICAL :: case_send_receive
1256 : REAL(KIND=dp) :: C_T_R, C_T_S
1257 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: BIb_RS_mat_rec
1258 0 : REAL(KIND=dp), DIMENSION(dimen, Rsize+Ssize) :: BIb_RS_mat
1259 :
1260 0 : CALL timeset(routineN, handle)
1261 :
1262 0 : S_offset = mp2_biel%index_table(latom, lset) - 1
1263 0 : R_offset = mp2_biel%index_table(katom, kset) - 1
1264 :
1265 0 : size_parameter_send(1) = Rsize
1266 0 : size_parameter_send(2) = Ssize
1267 0 : size_parameter_send(3) = R_offset
1268 0 : size_parameter_send(4) = S_offset
1269 0 : size_parameter_send(5) = elements_ij_proc
1270 :
1271 0 : DO index_proc_shift = 0, para_env%num_pe - 1
1272 :
1273 0 : proc_send = MODULO(para_env%mepos + index_proc_shift, para_env%num_pe)
1274 0 : proc_receive = MODULO(para_env%mepos - index_proc_shift, para_env%num_pe)
1275 :
1276 0 : case_send_receive = (proc_send /= para_env%mepos)
1277 :
1278 0 : IF (case_send_receive) THEN
1279 : ! the processor starts to send (and receive?)
1280 :
1281 0 : CALL para_env%sendrecv(size_parameter_send, proc_send, size_parameter_rec, proc_receive)
1282 :
1283 0 : Rsize_rec = size_parameter_rec(1)
1284 0 : Ssize_rec = size_parameter_rec(2)
1285 0 : R_offset_rec = size_parameter_rec(3)
1286 0 : S_offset_rec = size_parameter_rec(4)
1287 0 : elements_ij_proc_rec = size_parameter_rec(5)
1288 0 : ALLOCATE (BIb_RS_mat_rec(dimen, Rsize_rec + Ssize_rec))
1289 :
1290 : ELSE
1291 0 : elements_ij_proc_rec = elements_ij_proc
1292 : END IF
1293 :
1294 0 : index_ij_send = 0
1295 0 : index_ij_rec = 0
1296 0 : DO index_proc_ij = proc_send + 1, multiple, para_env%num_pe
1297 :
1298 0 : BIb_RS_mat = zero
1299 0 : IF (index_proc_ij <= Ni_occupied*Nj_occupied) THEN
1300 :
1301 0 : index_ij_send = index_ij_send + 1
1302 :
1303 0 : i = (index_proc_ij - 1)/Nj_occupied + 1
1304 0 : j = index_proc_ij - (i - 1)*Nj_occupied + j_batch_start
1305 :
1306 0 : DO S = 1, Ssize
1307 0 : S_global = S + S_offset
1308 0 : R_start = 1
1309 0 : IF (katom == latom .AND. kset == lset) R_start = S
1310 0 : DO R = R_start, Rsize
1311 0 : R_global = R + R_offset
1312 :
1313 0 : IF (R_global /= S_global) THEN
1314 0 : C_T_R = C_T(j, R_global)
1315 0 : C_T_S = C_T(j, S_global)
1316 0 : DO N = 1, dimen
1317 0 : BIb_RS_mat(N, R) = BIb_RS_mat(N, R) + C_T_S*BI1(N, i, R, S)
1318 : END DO
1319 0 : DO N = 1, dimen
1320 0 : BIb_RS_mat(N, Rsize + S) = BIb_RS_mat(N, Rsize + S) + C_T_R*BI1(N, i, R, S)
1321 : END DO
1322 : ELSE
1323 0 : C_T_S = C_T(j, S_global)
1324 0 : DO N = 1, dimen
1325 0 : BIb_RS_mat(N, R) = BIb_RS_mat(N, R) + C_T_S*BI1(N, i, R, S)
1326 : END DO
1327 : END IF
1328 :
1329 : END DO
1330 : END DO
1331 :
1332 : END IF
1333 :
1334 0 : IF (case_send_receive) THEN
1335 :
1336 0 : CALL para_env%sendrecv(BIb_RS_mat, proc_send, BIb_RS_mat_rec, proc_receive)
1337 :
1338 0 : index_ij_rec = index_ij_rec + 1
1339 0 : IF (index_ij_rec <= elements_ij_proc .AND. elements_ij_proc > 0) THEN
1340 :
1341 : BIb(1:dimen, R_offset_rec + 1:R_offset_rec + Rsize_rec, index_ij_rec) = &
1342 : BIb(1:dimen, R_offset_rec + 1:R_offset_rec + Rsize_rec, index_ij_rec) + &
1343 0 : BIb_RS_mat_rec(1:dimen, 1:Rsize_rec)
1344 :
1345 : BIb(1:dimen, S_offset_rec + 1:S_offset_rec + Ssize_rec, index_ij_rec) = &
1346 : BIb(1:dimen, S_offset_rec + 1:S_offset_rec + Ssize_rec, index_ij_rec) + &
1347 0 : BIb_RS_mat_rec(1:dimen, Rsize_rec + 1:Rsize_rec + Ssize_rec)
1348 :
1349 : END IF
1350 : ELSE
1351 : ! the processor is the sender and receiver itself
1352 0 : IF (index_ij_send <= elements_ij_proc .AND. elements_ij_proc > 0) THEN
1353 :
1354 : BIb(1:dimen, R_offset + 1:R_offset + Rsize, index_ij_send) = &
1355 0 : BIb(1:dimen, R_offset + 1:R_offset + Rsize, index_ij_send) + BIb_RS_mat(1:dimen, 1:Rsize)
1356 :
1357 : BIb(1:dimen, S_offset + 1:S_offset + Ssize, index_ij_send) = &
1358 0 : BIb(1:dimen, S_offset + 1:S_offset + Ssize, index_ij_send) + BIb_RS_mat(1:dimen, Rsize + 1:Rsize + Ssize)
1359 :
1360 : END IF
1361 : END IF
1362 :
1363 : END DO ! loop over the ij of the processor
1364 :
1365 0 : IF (case_send_receive) THEN
1366 0 : DEALLOCATE (BIb_RS_mat_rec)
1367 : END IF
1368 :
1369 : END DO ! loop over the processor starting from itself
1370 :
1371 0 : CALL timestop(handle)
1372 :
1373 0 : END SUBROUTINE transform_occupied_orbitals_second
1374 :
1375 : ! **************************************************************************************************
1376 : !> \brief ...
1377 : !> \param dimen ...
1378 : !> \param latom ...
1379 : !> \param katom ...
1380 : !> \param lset ...
1381 : !> \param kset ...
1382 : !> \param Ssize ...
1383 : !> \param Rsize ...
1384 : !> \param Ni_occupied ...
1385 : !> \param Nj_occupied ...
1386 : !> \param j_batch_start ...
1387 : !> \param ij_elem_max ...
1388 : !> \param BI1 ...
1389 : !> \param C_T ...
1390 : !> \param mp2_biel ...
1391 : !> \param para_env ...
1392 : !> \param elements_ij_proc ...
1393 : !> \param BIb ...
1394 : ! **************************************************************************************************
1395 4653 : SUBROUTINE transform_occupied_orbitals_second_big(dimen, latom, katom, lset, kset, &
1396 : Ssize, Rsize, Ni_occupied, Nj_occupied, j_batch_start, &
1397 4653 : ij_elem_max, BI1, C_T, mp2_biel, para_env, &
1398 4653 : elements_ij_proc, BIb)
1399 :
1400 : INTEGER, INTENT(IN) :: dimen, latom, katom, lset, kset, Ssize, &
1401 : Rsize, Ni_occupied, Nj_occupied, &
1402 : j_batch_start, ij_elem_max
1403 : REAL(KIND=dp), &
1404 : DIMENSION(dimen, Ni_occupied, Rsize, Ssize), &
1405 : INTENT(IN) :: BI1
1406 : REAL(KIND=dp), DIMENSION(dimen, dimen), INTENT(IN) :: C_T
1407 : TYPE(mp2_biel_type), INTENT(IN) :: mp2_biel
1408 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
1409 : INTEGER, INTENT(IN) :: elements_ij_proc
1410 : REAL(KIND=dp), &
1411 : DIMENSION(dimen, dimen, elements_ij_proc), &
1412 : INTENT(INOUT) :: BIb
1413 :
1414 : CHARACTER(LEN=*), PARAMETER :: routineN = 'transform_occupied_orbitals_second_big'
1415 :
1416 : INTEGER :: elements_ij_proc_rec, handle, i, index_ij_rec, index_ij_send, index_proc_ij, &
1417 : index_proc_shift, j, n, proc_receive, proc_send, r, R_global, R_offset, R_offset_rec, &
1418 : R_start, Rsize_rec, s, S_global, S_offset, S_offset_rec, Ssize_rec
1419 : INTEGER, DIMENSION(5) :: size_parameter_rec, size_parameter_send
1420 : LOGICAL :: case_send_receive
1421 : REAL(KIND=dp) :: C_T_R, C_T_S
1422 4653 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: BIb_RS_mat_rec
1423 : REAL(KIND=dp), &
1424 9306 : DIMENSION(dimen, Rsize+Ssize, ij_elem_max) :: BIb_RS_mat
1425 :
1426 4653 : CALL timeset(routineN, handle)
1427 :
1428 4653 : S_offset = mp2_biel%index_table(latom, lset) - 1
1429 4653 : R_offset = mp2_biel%index_table(katom, kset) - 1
1430 :
1431 4653 : size_parameter_send(1) = Rsize
1432 4653 : size_parameter_send(2) = Ssize
1433 4653 : size_parameter_send(3) = R_offset
1434 4653 : size_parameter_send(4) = S_offset
1435 4653 : size_parameter_send(5) = elements_ij_proc
1436 :
1437 9360 : DO index_proc_shift = 0, para_env%num_pe - 1
1438 :
1439 4707 : proc_send = MODULO(para_env%mepos + index_proc_shift, para_env%num_pe)
1440 4707 : proc_receive = MODULO(para_env%mepos - index_proc_shift, para_env%num_pe)
1441 :
1442 4707 : case_send_receive = (proc_send /= para_env%mepos)
1443 :
1444 4707 : IF (case_send_receive) THEN
1445 : ! the processor starts to send (and receive?)
1446 :
1447 54 : CALL para_env%sendrecv(size_parameter_send, proc_send, size_parameter_rec, proc_receive)
1448 :
1449 54 : Rsize_rec = size_parameter_rec(1)
1450 54 : Ssize_rec = size_parameter_rec(2)
1451 54 : R_offset_rec = size_parameter_rec(3)
1452 54 : S_offset_rec = size_parameter_rec(4)
1453 54 : elements_ij_proc_rec = size_parameter_rec(5)
1454 270 : ALLOCATE (BIb_RS_mat_rec(dimen, Rsize_rec + Ssize_rec, ij_elem_max))
1455 : ELSE
1456 4707 : elements_ij_proc_rec = elements_ij_proc
1457 : END IF
1458 :
1459 4707 : index_ij_send = 0
1460 4707 : index_ij_rec = 0
1461 30728964 : BIb_RS_mat = zero
1462 :
1463 81312 : DO index_proc_ij = proc_send + 1, Ni_occupied*Nj_occupied, para_env%num_pe
1464 :
1465 76605 : index_ij_send = index_ij_send + 1
1466 :
1467 76605 : i = (index_proc_ij - 1)/Nj_occupied + 1
1468 76605 : j = index_proc_ij - (i - 1)*Nj_occupied + j_batch_start
1469 :
1470 290238 : DO S = 1, Ssize
1471 208926 : S_global = S + S_offset
1472 208926 : R_start = 1
1473 208926 : IF (katom == latom .AND. kset == lset) R_start = S
1474 907160 : DO R = R_start, Rsize
1475 621629 : R_global = R + R_offset
1476 :
1477 830555 : IF (R_global /= S_global) THEN
1478 601864 : C_T_R = C_T(j, R_global)
1479 601864 : C_T_S = C_T(j, S_global)
1480 45769383 : DO N = 1, dimen
1481 45769383 : BIb_RS_mat(N, R, index_ij_send) = BIb_RS_mat(N, R, index_ij_send) + C_T_S*BI1(N, i, R, S)
1482 : END DO
1483 45769383 : DO N = 1, dimen
1484 45769383 : BIb_RS_mat(N, Rsize + S, index_ij_send) = BIb_RS_mat(N, Rsize + S, index_ij_send) + C_T_R*BI1(N, i, R, S)
1485 : END DO
1486 : ELSE
1487 19765 : C_T_S = C_T(j, S_global)
1488 1295658 : DO N = 1, dimen
1489 1295658 : BIb_RS_mat(N, R, index_ij_send) = BIb_RS_mat(N, R, index_ij_send) + C_T_S*BI1(N, i, R, S)
1490 : END DO
1491 : END IF
1492 :
1493 : END DO
1494 : END DO
1495 :
1496 : END DO
1497 :
1498 9360 : IF (case_send_receive) THEN
1499 :
1500 54 : CALL para_env%sendrecv(BIb_RS_mat, proc_send, BIb_RS_mat_rec, proc_receive)
1501 :
1502 : BIb(1:dimen, R_offset_rec + 1:R_offset_rec + Rsize_rec, 1:elements_ij_proc) = &
1503 : BIb(1:dimen, R_offset_rec + 1:R_offset_rec + Rsize_rec, 1:elements_ij_proc) + &
1504 16182 : BIb_RS_mat_rec(1:dimen, 1:Rsize_rec, 1:elements_ij_proc)
1505 :
1506 : BIb(1:dimen, S_offset_rec + 1:S_offset_rec + Ssize_rec, 1:elements_ij_proc) = &
1507 : BIb(1:dimen, S_offset_rec + 1:S_offset_rec + Ssize_rec, 1:elements_ij_proc) + &
1508 15006 : BIb_RS_mat_rec(1:dimen, Rsize_rec + 1:Rsize_rec + Ssize_rec, 1:elements_ij_proc)
1509 :
1510 54 : DEALLOCATE (BIb_RS_mat_rec)
1511 : ELSE
1512 : ! the processor is the sender and receiver itself
1513 : BIb(1:dimen, R_offset + 1:R_offset + Rsize, 1:elements_ij_proc) = &
1514 : BIb(1:dimen, R_offset + 1:R_offset + Rsize, 1:elements_ij_proc) + &
1515 16293273 : BIb_RS_mat(1:dimen, 1:Rsize, 1:elements_ij_proc)
1516 :
1517 : BIb(1:dimen, S_offset + 1:S_offset + Ssize, 1:elements_ij_proc) = &
1518 : BIb(1:dimen, S_offset + 1:S_offset + Ssize, 1:elements_ij_proc) + &
1519 14485815 : BIb_RS_mat(1:dimen, Rsize + 1:Rsize + Ssize, 1:elements_ij_proc)
1520 :
1521 : END IF
1522 :
1523 : END DO ! loop over the processor starting from itself
1524 :
1525 4653 : CALL timestop(handle)
1526 :
1527 4653 : END SUBROUTINE transform_occupied_orbitals_second_big
1528 :
1529 : ! **************************************************************************************************
1530 : !> \brief ...
1531 : !> \param dimen ...
1532 : !> \param occupied ...
1533 : !> \param virtual ...
1534 : !> \param i_batch_start ...
1535 : !> \param j_batch_start ...
1536 : !> \param BIb ...
1537 : !> \param C ...
1538 : !> \param Auto ...
1539 : !> \param elements_ij_proc ...
1540 : !> \param ij_list_proc ...
1541 : !> \param nspins ...
1542 : !> \param Emp2 ...
1543 : !> \param Emp2_Cou ...
1544 : !> \param Emp2_ex ...
1545 : ! **************************************************************************************************
1546 33 : SUBROUTINE transform_virtual_orbitals_and_accumulate(dimen, occupied, virtual, i_batch_start, &
1547 33 : j_batch_start, BIb, C, Auto, elements_ij_proc, &
1548 33 : ij_list_proc, nspins, Emp2, Emp2_Cou, Emp2_ex)
1549 :
1550 : INTEGER, INTENT(IN) :: dimen, occupied, virtual, i_batch_start, &
1551 : j_batch_start
1552 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
1553 : INTENT(INOUT) :: BIb
1554 : REAL(KIND=dp), DIMENSION(dimen, dimen), INTENT(IN) :: C
1555 : REAL(KIND=dp), DIMENSION(dimen), INTENT(IN) :: Auto
1556 : INTEGER, INTENT(IN) :: elements_ij_proc
1557 : INTEGER, DIMENSION(elements_ij_proc, 2), &
1558 : INTENT(IN) :: ij_list_proc
1559 : INTEGER, INTENT(IN) :: nspins
1560 : REAL(KIND=dp), INTENT(INOUT) :: Emp2, Emp2_Cou, Emp2_ex
1561 :
1562 : CHARACTER(LEN=*), PARAMETER :: routineN = 'transform_virtual_orbitals_and_accumulate'
1563 : REAL(KIND=dp), PARAMETER :: zero = 0.0_dp
1564 :
1565 : INTEGER :: a, a_global, b, b_global, handle, i, &
1566 : i_global, index_ij, j, j_global
1567 : REAL(KIND=dp) :: iajb, ibja, parz, two
1568 33 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: BIa
1569 :
1570 33 : CALL timeset(routineN, handle)
1571 :
1572 132 : ALLOCATE (BIa(dimen, virtual))
1573 :
1574 33 : BIa = zero
1575 389 : DO index_ij = 1, elements_ij_proc
1576 :
1577 : CALL DGEMM('T', 'N', dimen, virtual, dimen, 1.0_dp, Bib(1, 1, index_ij), &
1578 356 : dimen, C(1, occupied + 1), dimen, 0.0_dp, Bia(1, 1), dimen)
1579 1100036 : Bib(1:dimen, 1:virtual, index_ij) = Bia(1:dimen, 1:virtual)
1580 :
1581 : END DO
1582 :
1583 33 : DEALLOCATE (BIa)
1584 132 : ALLOCATE (BIa(virtual, virtual))
1585 :
1586 33 : BIa = zero
1587 389 : DO index_ij = 1, elements_ij_proc
1588 :
1589 : CALL DGEMM('T', 'N', virtual, virtual, dimen, 1.0_dp, Bib(1, 1, index_ij), dimen, C(1, occupied + 1), dimen, 0.0_dp, &
1590 356 : BIa(1, 1), virtual)
1591 980141 : BIb(1:virtual, 1:virtual, index_ij) = BIa(1:virtual, 1:virtual)
1592 :
1593 : END DO
1594 :
1595 33 : two = 2.0_dp/nspins
1596 389 : DO index_ij = 1, elements_ij_proc
1597 356 : i = ij_list_proc(index_ij, 1)
1598 356 : j = ij_list_proc(index_ij, 2)
1599 356 : i_global = i + i_batch_start
1600 356 : j_global = j + j_batch_start
1601 16600 : DO a = 1, virtual
1602 16211 : a_global = a + occupied
1603 980108 : DO b = 1, virtual
1604 963541 : b_global = b + occupied
1605 963541 : iajb = BIb(a, b, index_ij)
1606 963541 : ibja = BIb(b, a, index_ij)
1607 963541 : parz = iajb/(Auto(i_global) + Auto(j_global) - Auto(a_global) - Auto(b_global))
1608 : ! parz=parz*(two*iajb-ibja) !Full
1609 : ! parz=parz*(iajb) !Coulomb
1610 : ! parz=parz*(ibja) !Coulomb
1611 : ! Emp2=Emp2+parz/nspins
1612 963541 : Emp2_Cou = Emp2_Cou + parz*two*(iajb)/nspins
1613 963541 : Emp2_ex = Emp2_ex - parz*(ibja)/nspins
1614 979752 : Emp2 = Emp2 + parz*(two*iajb - ibja)/nspins
1615 : END DO
1616 : END DO
1617 : END DO
1618 :
1619 33 : DEALLOCATE (BIa)
1620 :
1621 33 : CALL timestop(handle)
1622 :
1623 33 : END SUBROUTINE transform_virtual_orbitals_and_accumulate
1624 :
1625 : ! **************************************************************************************************
1626 : !> \brief ...
1627 : !> \param dimen ...
1628 : !> \param occ_i ...
1629 : !> \param occ_j ...
1630 : !> \param virt_i ...
1631 : !> \param virt_j ...
1632 : !> \param i_batch_start ...
1633 : !> \param j_batch_start ...
1634 : !> \param BIb ...
1635 : !> \param C_i ...
1636 : !> \param C_j ...
1637 : !> \param Auto_i ...
1638 : !> \param Auto_j ...
1639 : !> \param elements_ij_proc ...
1640 : !> \param ij_list_proc ...
1641 : !> \param Emp2 ...
1642 : !> \param Emp2_Cou ...
1643 : ! **************************************************************************************************
1644 10 : SUBROUTINE transform_virtual_orbitals_and_accumulate_ABcase(dimen, occ_i, occ_j, virt_i, virt_j, i_batch_start, &
1645 10 : j_batch_start, BIb, C_i, C_j, Auto_i, Auto_j, elements_ij_proc, &
1646 10 : ij_list_proc, Emp2, Emp2_Cou)
1647 :
1648 : INTEGER, INTENT(IN) :: dimen, occ_i, occ_j, virt_i, virt_j, &
1649 : i_batch_start, j_batch_start
1650 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
1651 : INTENT(INOUT) :: BIb
1652 : REAL(KIND=dp), DIMENSION(dimen, dimen), INTENT(IN) :: C_i, C_j
1653 : REAL(KIND=dp), DIMENSION(dimen), INTENT(IN) :: Auto_i, Auto_j
1654 : INTEGER, INTENT(IN) :: elements_ij_proc
1655 : INTEGER, DIMENSION(elements_ij_proc, 2), &
1656 : INTENT(IN) :: ij_list_proc
1657 : REAL(KIND=dp), INTENT(INOUT) :: Emp2, Emp2_Cou
1658 :
1659 : CHARACTER(LEN=*), PARAMETER :: routineN = 'transform_virtual_orbitals_and_accumulate_ABcase'
1660 : REAL(KIND=dp), PARAMETER :: two = 2.D+00, zero = 0.D+00
1661 :
1662 : INTEGER :: a, a_global, b, b_global, handle, i, &
1663 : i_global, index_ij, j, j_global, n, s
1664 : REAL(KIND=dp) :: iajb, parz
1665 10 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: BIa
1666 :
1667 10 : CALL timeset(routineN, handle)
1668 :
1669 40 : ALLOCATE (BIa(dimen, virt_i))
1670 :
1671 56 : DO index_ij = 1, elements_ij_proc
1672 :
1673 1236 : DO a = 1, virt_i
1674 1190 : a_global = a + occ_i
1675 51930 : DO S = 1, dimen
1676 : parz = zero
1677 2449752 : DO N = 1, dimen
1678 2449752 : parz = parz + C_i(N, a_global)*BIb(N, S, index_ij)
1679 : END DO
1680 51884 : BIa(S, a) = parz
1681 : END DO
1682 : END DO
1683 :
1684 51940 : BIb(1:dimen, 1:virt_i, index_ij) = BIa
1685 :
1686 : END DO
1687 :
1688 10 : DEALLOCATE (BIa)
1689 40 : ALLOCATE (BIa(virt_i, virt_j))
1690 :
1691 56 : DO index_ij = 1, elements_ij_proc
1692 :
1693 1236 : DO a = 1, virt_i
1694 47824 : DO b = 1, virt_j
1695 46588 : b_global = b + occ_j
1696 46588 : parz = zero
1697 2257144 : DO S = 1, dimen
1698 2257144 : parz = parz + C_j(S, b_global)*BIb(S, a, index_ij)
1699 : END DO
1700 47778 : BIa(a, b) = parz
1701 : END DO
1702 : END DO
1703 :
1704 47904 : BIb(1:virt_i, 1:virt_j, index_ij) = BIa
1705 :
1706 : END DO
1707 :
1708 56 : DO index_ij = 1, elements_ij_proc
1709 46 : i = ij_list_proc(index_ij, 1)
1710 46 : j = ij_list_proc(index_ij, 2)
1711 46 : i_global = i + i_batch_start
1712 46 : j_global = j + j_batch_start
1713 1246 : DO a = 1, virt_i
1714 1190 : a_global = a + occ_i
1715 47824 : DO b = 1, virt_j
1716 46588 : b_global = b + occ_j
1717 46588 : iajb = BIb(a, b, index_ij)
1718 46588 : parz = iajb*iajb/(Auto_i(i_global) + Auto_j(j_global) - Auto_i(a_global) - Auto_j(b_global))
1719 46588 : Emp2_Cou = Emp2_Cou + parz/two
1720 47778 : Emp2 = Emp2 + parz/two
1721 : END DO
1722 : END DO
1723 : END DO
1724 :
1725 10 : DEALLOCATE (BIa)
1726 :
1727 10 : CALL timestop(handle)
1728 :
1729 10 : END SUBROUTINE transform_virtual_orbitals_and_accumulate_ABcase
1730 :
1731 : END MODULE mp2_direct_method
|
