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