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 the 3 and 2 center ERI's needed in the RI
10 : !> approximation using libint
11 : !> \par History
12 : !> 08.2013 created [Mauro Del Ben]
13 : !> \author Mauro Del Ben
14 : ! **************************************************************************************************
15 : MODULE mp2_ri_libint
16 : USE atomic_kind_types, ONLY: atomic_kind_type,&
17 : get_atomic_kind_set
18 : USE basis_set_types, ONLY: get_gto_basis_set,&
19 : gto_basis_set_type,&
20 : init_aux_basis_set
21 : USE cell_types, ONLY: cell_type
22 : USE cp_blacs_env, ONLY: cp_blacs_env_create,&
23 : cp_blacs_env_release,&
24 : cp_blacs_env_type
25 : USE cp_control_types, ONLY: dft_control_type
26 : USE cp_files, ONLY: close_file,&
27 : open_file
28 : USE cp_fm_basic_linalg, ONLY: cp_fm_triangular_invert
29 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose
30 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
31 : cp_fm_struct_release,&
32 : cp_fm_struct_type
33 : USE cp_fm_types, ONLY: cp_fm_create,&
34 : cp_fm_get_submatrix,&
35 : cp_fm_release,&
36 : cp_fm_set_submatrix,&
37 : cp_fm_type
38 : USE four_center_integral_interface, ONLY: four_center_integral_context_type
39 : USE four_center_integral_methods, ONLY: coulomb4
40 : USE gamma, ONLY: init_md_ftable
41 : USE hfx_pair_list_methods, ONLY: build_pair_list_mp2
42 : USE hfx_screening_methods, ONLY: calc_pair_dist_radii,&
43 : calc_screening_functions
44 : USE hfx_types, ONLY: &
45 : hfx_basis_info_type, hfx_basis_type, hfx_create_neighbor_cells, hfx_pgf_list, &
46 : hfx_pgf_product_list, hfx_potential_type, hfx_screen_coeff_type, hfx_type, log_zero, &
47 : pair_set_list_type
48 : USE input_constants, ONLY: do_potential_TShPSC,&
49 : do_potential_long
50 : USE kinds, ONLY: dp,&
51 : int_8
52 : USE message_passing, ONLY: mp_para_env_type
53 : USE mp2_types, ONLY: init_TShPSC_lmax,&
54 : mp2_biel_type,&
55 : mp2_type,&
56 : pair_list_type_mp2
57 : USE orbital_pointers, ONLY: nco,&
58 : ncoset,&
59 : nso
60 : USE particle_types, ONLY: particle_type
61 : USE qs_environment_types, ONLY: get_qs_env,&
62 : qs_environment_type
63 : USE qs_kind_types, ONLY: get_qs_kind,&
64 : get_qs_kind_set,&
65 : qs_kind_type
66 : USE t_sh_p_s_c, ONLY: free_C0,&
67 : init
68 : #include "./base/base_uses.f90"
69 :
70 : IMPLICIT NONE
71 : PRIVATE
72 :
73 : PUBLIC :: libint_ri_mp2, read_RI_basis_set, release_RI_basis_set, prepare_integral_calc
74 :
75 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mp2_ri_libint'
76 :
77 : !***
78 :
79 : CONTAINS
80 :
81 : ! **************************************************************************************************
82 : !> \brief ...
83 : !> \param dimen ...
84 : !> \param RI_dimen ...
85 : !> \param occupied ...
86 : !> \param natom ...
87 : !> \param mp2_biel ...
88 : !> \param mp2_env ...
89 : !> \param C ...
90 : !> \param kind_of ...
91 : !> \param RI_basis_parameter ...
92 : !> \param RI_basis_info ...
93 : !> \param basis_S0 ...
94 : !> \param RI_index_table ...
95 : !> \param qs_env ...
96 : !> \param para_env ...
97 : !> \param Lai ...
98 : ! **************************************************************************************************
99 5816 : SUBROUTINE libint_ri_mp2(dimen, RI_dimen, occupied, natom, mp2_biel, mp2_env, C, &
100 2908 : kind_of, &
101 : RI_basis_parameter, RI_basis_info, basis_S0, RI_index_table, &
102 : qs_env, para_env, &
103 : Lai)
104 : INTEGER :: dimen, RI_dimen, occupied, natom
105 : TYPE(mp2_biel_type) :: mp2_biel
106 : TYPE(mp2_type) :: mp2_env
107 : REAL(KIND=dp), DIMENSION(dimen, dimen) :: C
108 : INTEGER, DIMENSION(:) :: kind_of
109 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: RI_basis_parameter
110 : TYPE(hfx_basis_info_type) :: RI_basis_info
111 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_S0
112 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: RI_index_table
113 : TYPE(qs_environment_type), POINTER :: qs_env
114 : TYPE(mp_para_env_type), POINTER :: para_env
115 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Lai
116 :
117 : CHARACTER(LEN=*), PARAMETER :: routineN = 'libint_ri_mp2'
118 :
119 : INTEGER :: handle, nkind
120 :
121 2908 : CALL timeset(routineN, handle)
122 :
123 : ! Get the RI basis set and store in to a nice form
124 2908 : IF (.NOT. (ASSOCIATED(RI_basis_parameter))) THEN
125 0 : IF (ALLOCATED(RI_index_table)) DEALLOCATE (RI_index_table)
126 0 : IF (ASSOCIATED(basis_S0)) DEALLOCATE (basis_S0)
127 : CALL read_RI_basis_set(qs_env, RI_basis_parameter, RI_basis_info, &
128 : natom, nkind, kind_of, RI_index_table, RI_dimen, &
129 0 : basis_S0)
130 : END IF
131 :
132 : CALL calc_lai_libint(mp2_env, qs_env, para_env, &
133 : mp2_biel, dimen, C, occupied, &
134 : RI_basis_parameter, RI_basis_info, RI_index_table, RI_dimen, basis_S0, &
135 2908 : Lai)
136 :
137 2908 : CALL timestop(handle)
138 :
139 2908 : END SUBROUTINE libint_ri_mp2
140 :
141 : ! **************************************************************************************************
142 : !> \brief Read the auxiliary basis set for RI approxiamtion
143 : !> \param qs_env ...
144 : !> \param RI_basis_parameter ...
145 : !> \param RI_basis_info ...
146 : !> \param natom ...
147 : !> \param nkind ...
148 : !> \param kind_of ...
149 : !> \param RI_index_table ...
150 : !> \param RI_dimen ...
151 : !> \param basis_S0 ...
152 : !> \par History
153 : !> 08.2013 created [Mauro Del Ben]
154 : !> \author Mauro Del Ben
155 : ! **************************************************************************************************
156 6 : SUBROUTINE read_RI_basis_set(qs_env, RI_basis_parameter, RI_basis_info, &
157 6 : natom, nkind, kind_of, RI_index_table, RI_dimen, &
158 : basis_S0)
159 : TYPE(qs_environment_type), POINTER :: qs_env
160 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: RI_basis_parameter
161 : TYPE(hfx_basis_info_type) :: RI_basis_info
162 : INTEGER :: natom, nkind
163 : INTEGER, DIMENSION(:) :: kind_of
164 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: RI_index_table
165 : INTEGER :: RI_dimen
166 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_S0
167 :
168 : INTEGER :: co_counter, counter, i, iatom, ikind, ipgf, iset, j, k, la, max_am_kind, &
169 : max_coeff, max_nsgfl, max_pgf, max_pgf_kind, max_set, nl_count, nset, nseta, offset_a, &
170 : offset_a1, s_offset_nl_a, sgfa, so_counter
171 6 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nshell
172 6 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, nl_a
173 6 : REAL(dp), DIMENSION(:, :), POINTER :: sphi_a
174 : TYPE(gto_basis_set_type), POINTER :: orb_basis_a
175 6 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
176 : TYPE(qs_kind_type), POINTER :: atom_kind
177 :
178 6 : NULLIFY (RI_basis_parameter)
179 :
180 6 : NULLIFY (qs_kind_set)
181 6 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
182 :
183 6 : nkind = SIZE(qs_kind_set, 1)
184 24 : ALLOCATE (RI_basis_parameter(nkind))
185 18 : ALLOCATE (basis_S0(nkind))
186 12 : max_set = 0
187 12 : DO ikind = 1, nkind
188 6 : NULLIFY (atom_kind)
189 6 : atom_kind => qs_kind_set(ikind)
190 : ! here we reset the initial RI basis such that we can
191 : ! work with non-normalized auxiliary basis functions
192 : CALL get_qs_kind(qs_kind=atom_kind, &
193 6 : basis_set=orb_basis_a, basis_type="RI_AUX")
194 6 : IF (.NOT. (ASSOCIATED(orb_basis_a))) THEN
195 0 : CPABORT("Initial RI auxiliary basis not specified.")
196 : END IF
197 264 : orb_basis_a%gcc = 1.0_dp
198 6 : orb_basis_a%norm_type = 1
199 6 : CALL init_aux_basis_set(orb_basis_a)
200 :
201 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
202 : maxsgf=RI_basis_info%max_sgf, &
203 : maxnset=RI_basis_info%max_set, &
204 : maxlgto=RI_basis_info%max_am, &
205 6 : basis_type="RI_AUX")
206 : CALL get_gto_basis_set(gto_basis_set=orb_basis_a, &
207 : lmax=RI_basis_parameter(ikind)%lmax, &
208 : lmin=RI_basis_parameter(ikind)%lmin, &
209 : npgf=RI_basis_parameter(ikind)%npgf, &
210 : nset=RI_basis_parameter(ikind)%nset, &
211 : zet=RI_basis_parameter(ikind)%zet, &
212 : nsgf_set=RI_basis_parameter(ikind)%nsgf, &
213 : first_sgf=RI_basis_parameter(ikind)%first_sgf, &
214 : sphi=RI_basis_parameter(ikind)%sphi, &
215 : nsgf=RI_basis_parameter(ikind)%nsgf_total, &
216 : l=RI_basis_parameter(ikind)%nl, &
217 : nshell=RI_basis_parameter(ikind)%nshell, &
218 : set_radius=RI_basis_parameter(ikind)%set_radius, &
219 : pgf_radius=RI_basis_parameter(ikind)%pgf_radius, &
220 6 : kind_radius=RI_basis_parameter(ikind)%kind_radius)
221 :
222 6 : max_set = MAX(max_set, RI_basis_parameter(ikind)%nset)
223 :
224 6 : basis_S0(ikind)%kind_radius = RI_basis_parameter(ikind)%kind_radius
225 6 : basis_S0(ikind)%nset = 1
226 6 : basis_S0(ikind)%nsgf_total = 1
227 6 : ALLOCATE (basis_S0(ikind)%set_radius(1))
228 6 : basis_S0(ikind)%set_radius(1) = RI_basis_parameter(ikind)%kind_radius
229 6 : ALLOCATE (basis_S0(ikind)%lmax(1))
230 6 : basis_S0(ikind)%lmax(1) = 0
231 6 : ALLOCATE (basis_S0(ikind)%lmin(1))
232 6 : basis_S0(ikind)%lmin(1) = 0
233 6 : ALLOCATE (basis_S0(ikind)%npgf(1))
234 6 : basis_S0(ikind)%npgf(1) = 1
235 6 : ALLOCATE (basis_S0(ikind)%nsgf(1))
236 6 : basis_S0(ikind)%nsgf(1) = 1
237 6 : ALLOCATE (basis_S0(ikind)%nshell(1))
238 6 : basis_S0(ikind)%nshell(1) = 1
239 6 : ALLOCATE (basis_S0(ikind)%pgf_radius(1, 1))
240 6 : basis_S0(ikind)%pgf_radius(1, 1) = RI_basis_parameter(ikind)%kind_radius
241 6 : ALLOCATE (basis_S0(ikind)%sphi(1, 1))
242 6 : basis_S0(ikind)%sphi(1, 1) = 1.0_dp
243 6 : ALLOCATE (basis_S0(ikind)%zet(1, 1))
244 6 : basis_S0(ikind)%zet(1, 1) = 0.0_dp
245 6 : ALLOCATE (basis_S0(ikind)%first_sgf(1, 1))
246 6 : basis_S0(ikind)%first_sgf(1, 1) = 1
247 6 : ALLOCATE (basis_S0(ikind)%nl(1, 1))
248 6 : basis_S0(ikind)%nl(1, 1) = 0
249 :
250 6 : ALLOCATE (basis_S0(ikind)%nsgfl(0:0, 1))
251 18 : basis_S0(ikind)%nsgfl = 1
252 6 : ALLOCATE (basis_S0(ikind)%sphi_ext(1, 0:0, 1, 1))
253 12 : basis_S0(ikind)%sphi_ext(1, 0, 1, 1) = 1.0_dp
254 :
255 : END DO
256 6 : RI_basis_info%max_set = max_set
257 :
258 12 : DO ikind = 1, nkind
259 24 : ALLOCATE (RI_basis_parameter(ikind)%nsgfl(0:RI_basis_info%max_am, max_set))
260 436 : RI_basis_parameter(ikind)%nsgfl = 0
261 6 : nset = RI_basis_parameter(ikind)%nset
262 6 : nshell => RI_basis_parameter(ikind)%nshell
263 98 : DO iset = 1, nset
264 436 : DO i = 0, RI_basis_info%max_am
265 344 : nl_count = 0
266 688 : DO j = 1, nshell(iset)
267 688 : IF (RI_basis_parameter(ikind)%nl(j, iset) == i) nl_count = nl_count + 1
268 : END DO
269 430 : RI_basis_parameter(ikind)%nsgfl(i, iset) = nl_count
270 : END DO
271 : END DO
272 : END DO
273 :
274 : max_nsgfl = 0
275 : max_pgf = 0
276 12 : DO ikind = 1, nkind
277 6 : max_coeff = 0
278 6 : max_am_kind = 0
279 6 : max_pgf_kind = 0
280 6 : npgfa => RI_basis_parameter(ikind)%npgf
281 6 : nseta = RI_basis_parameter(ikind)%nset
282 6 : nl_a => RI_basis_parameter(ikind)%nsgfl
283 6 : la_max => RI_basis_parameter(ikind)%lmax
284 6 : la_min => RI_basis_parameter(ikind)%lmin
285 92 : DO iset = 1, nseta
286 86 : max_pgf_kind = MAX(max_pgf_kind, npgfa(iset))
287 : max_pgf = MAX(max_pgf, npgfa(iset))
288 178 : DO la = la_min(iset), la_max(iset)
289 86 : max_nsgfl = MAX(max_nsgfl, nl_a(la, iset))
290 86 : max_coeff = MAX(max_coeff, nso(la)*nl_a(la, iset)*nco(la))
291 172 : max_am_kind = MAX(max_am_kind, la)
292 : END DO
293 : END DO
294 36 : ALLOCATE (RI_basis_parameter(ikind)%sphi_ext(max_coeff, 0:max_am_kind, max_pgf_kind, nseta))
295 24608 : RI_basis_parameter(ikind)%sphi_ext = 0.0_dp
296 : END DO
297 :
298 12 : DO ikind = 1, nkind
299 6 : sphi_a => RI_basis_parameter(ikind)%sphi
300 6 : nseta = RI_basis_parameter(ikind)%nset
301 6 : la_max => RI_basis_parameter(ikind)%lmax
302 6 : la_min => RI_basis_parameter(ikind)%lmin
303 6 : npgfa => RI_basis_parameter(ikind)%npgf
304 6 : first_sgfa => RI_basis_parameter(ikind)%first_sgf
305 6 : nl_a => RI_basis_parameter(ikind)%nsgfl
306 98 : DO iset = 1, nseta
307 86 : sgfa = first_sgfa(1, iset)
308 178 : DO ipgf = 1, npgfa(iset)
309 86 : offset_a1 = (ipgf - 1)*ncoset(la_max(iset))
310 86 : s_offset_nl_a = 0
311 258 : DO la = la_min(iset), la_max(iset)
312 86 : offset_a = offset_a1 + ncoset(la - 1)
313 : co_counter = 0
314 86 : co_counter = co_counter + 1
315 86 : so_counter = 0
316 352 : DO k = sgfa + s_offset_nl_a, sgfa + s_offset_nl_a + nso(la)*nl_a(la, iset) - 1
317 1778 : DO i = offset_a + 1, offset_a + nco(la)
318 1426 : so_counter = so_counter + 1
319 1692 : RI_basis_parameter(ikind)%sphi_ext(so_counter, la, ipgf, iset) = sphi_a(i, k)
320 : END DO
321 : END DO
322 172 : s_offset_nl_a = s_offset_nl_a + nso(la)*(nl_a(la, iset))
323 : END DO
324 : END DO
325 : END DO
326 : END DO
327 :
328 24 : ALLOCATE (RI_index_table(natom, max_set))
329 178 : RI_index_table = -HUGE(0)
330 6 : counter = 0
331 6 : RI_dimen = 0
332 12 : DO iatom = 1, natom
333 6 : ikind = kind_of(iatom)
334 6 : nset = RI_basis_parameter(ikind)%nset
335 98 : DO iset = 1, nset
336 86 : RI_index_table(iatom, iset) = counter + 1
337 86 : counter = counter + RI_basis_parameter(ikind)%nsgf(iset)
338 92 : RI_dimen = RI_dimen + RI_basis_parameter(ikind)%nsgf(iset)
339 : END DO
340 : END DO
341 :
342 6 : END SUBROUTINE read_RI_basis_set
343 :
344 : ! **************************************************************************************************
345 : !> \brief Release the auxiliary basis set for RI approxiamtion (to be used
346 : !> only in the case of basis optimization)
347 : !> \param RI_basis_parameter ...
348 : !> \param basis_S0 ...
349 : !> \par History
350 : !> 08.2013 created [Mauro Del Ben]
351 : !> \author Mauro Del Ben
352 : ! **************************************************************************************************
353 6 : SUBROUTINE release_RI_basis_set(RI_basis_parameter, basis_S0)
354 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: RI_basis_parameter, basis_S0
355 :
356 : INTEGER :: i
357 :
358 : ! RI basis
359 :
360 12 : DO i = 1, SIZE(RI_basis_parameter)
361 6 : DEALLOCATE (RI_basis_parameter(i)%nsgfl)
362 12 : DEALLOCATE (RI_basis_parameter(i)%sphi_ext)
363 : END DO
364 6 : DEALLOCATE (RI_basis_parameter)
365 :
366 : ! S0 basis
367 12 : DO i = 1, SIZE(basis_S0)
368 6 : DEALLOCATE (basis_S0(i)%set_radius)
369 6 : DEALLOCATE (basis_S0(i)%lmax)
370 6 : DEALLOCATE (basis_S0(i)%lmin)
371 6 : DEALLOCATE (basis_S0(i)%npgf)
372 6 : DEALLOCATE (basis_S0(i)%nsgf)
373 6 : DEALLOCATE (basis_S0(i)%nshell)
374 6 : DEALLOCATE (basis_S0(i)%pgf_radius)
375 6 : DEALLOCATE (basis_S0(i)%sphi)
376 6 : DEALLOCATE (basis_S0(i)%zet)
377 6 : DEALLOCATE (basis_S0(i)%first_sgf)
378 6 : DEALLOCATE (basis_S0(i)%nl)
379 6 : DEALLOCATE (basis_S0(i)%nsgfl)
380 12 : DEALLOCATE (basis_S0(i)%sphi_ext)
381 : END DO
382 6 : DEALLOCATE (basis_S0)
383 :
384 6 : END SUBROUTINE release_RI_basis_set
385 :
386 : ! **************************************************************************************************
387 : !> \brief ...
388 : !> \param mp2_env ...
389 : !> \param qs_env ...
390 : !> \param para_env ...
391 : !> \param mp2_biel ...
392 : !> \param dimen ...
393 : !> \param C ...
394 : !> \param occupied ...
395 : !> \param RI_basis_parameter ...
396 : !> \param RI_basis_info ...
397 : !> \param RI_index_table ...
398 : !> \param RI_dimen ...
399 : !> \param basis_S0 ...
400 : !> \param Lai ...
401 : !> \par History
402 : !> 08.2013 created [Mauro Del Ben]
403 : !> \author Mauro Del Ben
404 : ! **************************************************************************************************
405 2908 : SUBROUTINE calc_lai_libint(mp2_env, qs_env, para_env, &
406 2908 : mp2_biel, dimen, C, occupied, &
407 : RI_basis_parameter, RI_basis_info, RI_index_table, RI_dimen, basis_S0, &
408 : Lai)
409 :
410 : TYPE(mp2_type) :: mp2_env
411 : TYPE(qs_environment_type), POINTER :: qs_env
412 : TYPE(mp_para_env_type), POINTER :: para_env
413 : TYPE(mp2_biel_type) :: mp2_biel
414 : INTEGER :: dimen
415 : REAL(KIND=dp), DIMENSION(dimen, dimen) :: C
416 : INTEGER :: occupied
417 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: RI_basis_parameter
418 : TYPE(hfx_basis_info_type) :: RI_basis_info
419 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: RI_index_table
420 : INTEGER :: RI_dimen
421 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_S0
422 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Lai
423 :
424 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_lai_libint'
425 :
426 : INTEGER :: counter_L_blocks, handle, i, i_list_kl, i_set_list_kl, i_set_list_kl_start, &
427 : i_set_list_kl_stop, iatom, iatom_end, iatom_start, iiB, ikind, info_chol, iset, jatom, &
428 : jatom_end, jatom_start, jjB, jkind, jset, katom, katom_end, katom_start, kkB, kkind, &
429 : kset, kset_start, L_B_i_end, L_B_i_start, L_B_k_end, L_B_k_start, latom, latom_end, &
430 : latom_start, lkind, llB, lset, max_pgf, max_set, natom, ncob, nints, nseta, nsetb, nsetc, &
431 : nsetd, nsgf_max, nspins, orb_k_end, orb_k_start, orb_l_end, orb_l_start, &
432 : primitive_counter, sphi_a_u1, sphi_a_u2, sphi_a_u3, sphi_b_u1, sphi_b_u2, sphi_b_u3
433 : INTEGER :: sphi_c_u1, sphi_c_u2, sphi_c_u3, sphi_d_u1, sphi_d_u2, sphi_d_u3, virtual
434 : INTEGER(int_8) :: estimate_to_store_int, neris_tmp, &
435 : neris_total, nprim_ints
436 2908 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of, nimages
437 2908 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, lc_max, &
438 2908 : lc_min, ld_max, ld_min, npgfa, npgfb, &
439 2908 : npgfc, npgfd, nsgfa, nsgfb, nsgfc, &
440 2908 : nsgfd
441 2908 : INTEGER, DIMENSION(:, :), POINTER :: nsgfl_a, nsgfl_b, nsgfl_c, nsgfl_d
442 : LOGICAL :: do_periodic
443 2908 : LOGICAL, DIMENSION(:, :), POINTER :: shm_atomic_pair_list
444 : REAL(KIND=dp) :: cartesian_estimate, coeffs_kind_max0, eps_schwarz, ln_10, &
445 : log10_eps_schwarz, log10_pmax, max_contraction_val, pmax_atom, pmax_entry, ra(3), rb(3), &
446 : rc(3), rd(3)
447 2908 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ee_buffer1, ee_buffer2, &
448 2908 : ee_primitives_tmp, ee_work, ee_work2, &
449 2908 : primitive_integrals
450 2908 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: L_block, L_full_matrix, max_contraction
451 2908 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: BI1, MNRS
452 2908 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: shm_pmax_block, zeta, zetb, zetc, zetd
453 2908 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: sphi_a_ext_set, sphi_b_ext_set, &
454 2908 : sphi_c_ext_set, sphi_d_ext_set
455 2908 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: sphi_a_ext, sphi_b_ext, sphi_c_ext, &
456 2908 : sphi_d_ext
457 : TYPE(cell_type), POINTER :: cell
458 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
459 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_L
460 : TYPE(cp_fm_type) :: fm_matrix_L
461 : TYPE(four_center_integral_context_type) :: private_context
462 2908 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter
463 2908 : TYPE(hfx_pgf_list), ALLOCATABLE, DIMENSION(:) :: pgf_list_ij, pgf_list_kl
464 : TYPE(hfx_pgf_product_list), ALLOCATABLE, &
465 2908 : DIMENSION(:) :: pgf_product_list
466 : TYPE(hfx_potential_type) :: mp2_potential_parameter
467 : TYPE(hfx_screen_coeff_type), ALLOCATABLE, &
468 2908 : DIMENSION(:, :), TARGET :: radii_pgf_large, screen_coeffs_pgf_large
469 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
470 2908 : POINTER :: screen_coeffs_kind, tmp_R_1, tmp_R_2, &
471 2908 : tmp_screen_pgf1, tmp_screen_pgf2
472 : TYPE(hfx_screen_coeff_type), &
473 2908 : DIMENSION(:, :, :, :), POINTER :: screen_coeffs_set
474 : TYPE(hfx_screen_coeff_type), &
475 2908 : DIMENSION(:, :, :, :, :, :), POINTER :: radii_pgf, screen_coeffs_pgf
476 : TYPE(hfx_type), POINTER :: actual_x_data
477 2908 : TYPE(pair_list_type_mp2) :: list_kl
478 : TYPE(pair_set_list_type), ALLOCATABLE, &
479 2908 : DIMENSION(:) :: set_list_kl
480 2908 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
481 :
482 2908 : CALL timeset(routineN, handle)
483 :
484 2908 : ln_10 = LOG(10.0_dp)
485 :
486 : !! initialize some counters
487 2908 : neris_tmp = 0_int_8
488 2908 : neris_total = 0_int_8
489 2908 : nprim_ints = 0_int_8
490 :
491 : CALL prepare_integral_calc(cell, qs_env, mp2_env, para_env, mp2_potential_parameter, actual_x_data, &
492 : do_periodic, basis_parameter, max_set, particle_set, natom, kind_of, &
493 : nsgf_max, primitive_integrals, ee_work, ee_work2, ee_buffer1, ee_buffer2, &
494 : ee_primitives_tmp, nspins, max_contraction, max_pgf, pgf_list_ij, &
495 : pgf_list_kl, pgf_product_list, nimages, eps_schwarz, log10_eps_schwarz, &
496 : private_context, screen_coeffs_set, screen_coeffs_kind, screen_coeffs_pgf, &
497 2908 : radii_pgf, RI_basis_parameter, RI_basis_info)
498 :
499 52344 : ALLOCATE (radii_pgf_large(SIZE(radii_pgf, 1), SIZE(radii_pgf, 2)))
500 52344 : ALLOCATE (screen_coeffs_pgf_large(SIZE(screen_coeffs_pgf, 1), SIZE(screen_coeffs_pgf, 2)))
501 11632 : DO iiB = 1, SIZE(radii_pgf, 1)
502 37804 : DO jjB = 1, SIZE(radii_pgf, 2)
503 87240 : radii_pgf_large(iiB, jjB)%x = 100_dp
504 : END DO
505 : END DO
506 11632 : DO iiB = 1, SIZE(screen_coeffs_pgf, 1)
507 37804 : DO jjB = 1, SIZE(screen_coeffs_pgf, 2)
508 87240 : screen_coeffs_pgf_large(iiB, jjB)%x = 5000_dp
509 : END DO
510 : END DO
511 2908 : tmp_R_1 => radii_pgf_large(:, :)
512 2908 : tmp_R_2 => radii_pgf_large(:, :)
513 2908 : tmp_screen_pgf1 => screen_coeffs_pgf_large(:, :)
514 2908 : tmp_screen_pgf2 => screen_coeffs_pgf_large(:, :)
515 :
516 : ! start computing the L matrix
517 11632 : ALLOCATE (L_full_matrix(RI_dimen, RI_dimen))
518 2908 : L_full_matrix = 0.0_dp
519 :
520 2908 : counter_L_blocks = 0
521 5816 : DO iatom = 1, natom
522 2908 : jatom = iatom
523 :
524 2908 : ikind = kind_of(iatom)
525 2908 : jkind = kind_of(jatom)
526 11632 : ra = particle_set(iatom)%r(:)
527 11632 : rb = particle_set(jatom)%r(:)
528 :
529 2908 : la_max => RI_basis_parameter(ikind)%lmax
530 2908 : la_min => RI_basis_parameter(ikind)%lmin
531 2908 : npgfa => RI_basis_parameter(ikind)%npgf
532 2908 : nseta = RI_basis_parameter(ikind)%nset
533 2908 : zeta => RI_basis_parameter(ikind)%zet
534 2908 : nsgfa => RI_basis_parameter(ikind)%nsgf
535 2908 : sphi_a_ext => RI_basis_parameter(ikind)%sphi_ext(:, :, :, :)
536 2908 : nsgfl_a => RI_basis_parameter(ikind)%nsgfl
537 2908 : sphi_a_u1 = UBOUND(sphi_a_ext, 1)
538 2908 : sphi_a_u2 = UBOUND(sphi_a_ext, 2)
539 2908 : sphi_a_u3 = UBOUND(sphi_a_ext, 3)
540 :
541 2908 : lb_max => basis_S0(jkind)%lmax
542 2908 : lb_min => basis_S0(jkind)%lmin
543 2908 : npgfb => basis_S0(jkind)%npgf
544 2908 : nsetb = basis_S0(jkind)%nset
545 2908 : zetb => basis_S0(jkind)%zet
546 2908 : nsgfb => basis_S0(jkind)%nsgf
547 2908 : sphi_b_ext => basis_S0(jkind)%sphi_ext(:, :, :, :)
548 2908 : nsgfl_b => basis_S0(jkind)%nsgfl
549 2908 : sphi_b_u1 = UBOUND(sphi_b_ext, 1)
550 2908 : sphi_b_u2 = UBOUND(sphi_b_ext, 2)
551 2908 : sphi_b_u3 = UBOUND(sphi_b_ext, 3)
552 :
553 8724 : DO katom = iatom, natom
554 2908 : latom = katom
555 :
556 2908 : kkind = kind_of(katom)
557 2908 : lkind = kind_of(latom)
558 11632 : rc = particle_set(katom)%r(:)
559 11632 : rd = particle_set(latom)%r(:)
560 :
561 2908 : lc_max => RI_basis_parameter(kkind)%lmax
562 2908 : lc_min => RI_basis_parameter(kkind)%lmin
563 2908 : npgfc => RI_basis_parameter(kkind)%npgf
564 2908 : nsetc = RI_basis_parameter(kkind)%nset
565 2908 : zetc => RI_basis_parameter(kkind)%zet
566 2908 : nsgfc => RI_basis_parameter(kkind)%nsgf
567 2908 : sphi_c_ext => RI_basis_parameter(kkind)%sphi_ext(:, :, :, :)
568 2908 : nsgfl_c => RI_basis_parameter(kkind)%nsgfl
569 2908 : sphi_c_u1 = UBOUND(sphi_c_ext, 1)
570 2908 : sphi_c_u2 = UBOUND(sphi_c_ext, 2)
571 2908 : sphi_c_u3 = UBOUND(sphi_c_ext, 3)
572 :
573 2908 : ld_max => basis_S0(lkind)%lmax
574 2908 : ld_min => basis_S0(lkind)%lmin
575 2908 : npgfd => basis_S0(lkind)%npgf
576 2908 : nsetd = RI_basis_parameter(lkind)%nset
577 2908 : zetd => basis_S0(lkind)%zet
578 2908 : nsgfd => basis_S0(lkind)%nsgf
579 2908 : sphi_d_ext => basis_S0(lkind)%sphi_ext(:, :, :, :)
580 2908 : nsgfl_d => basis_S0(lkind)%nsgfl
581 2908 : sphi_d_u1 = UBOUND(sphi_d_ext, 1)
582 2908 : sphi_d_u2 = UBOUND(sphi_d_ext, 2)
583 2908 : sphi_d_u3 = UBOUND(sphi_d_ext, 3)
584 :
585 2908 : jset = 1
586 2908 : lset = 1
587 45948 : DO iset = 1, nseta
588 :
589 40132 : ncob = npgfb(jset)*ncoset(lb_max(jset))
590 40132 : sphi_a_ext_set => sphi_a_ext(:, :, :, iset)
591 40132 : sphi_b_ext_set => sphi_b_ext(:, :, :, jset)
592 :
593 40132 : L_B_i_start = RI_index_table(iatom, iset)
594 40132 : L_B_i_end = RI_index_table(iatom, iset) + nsgfa(iset) - 1
595 :
596 40132 : kset_start = 1
597 40132 : IF (iatom == katom) kset_start = iset
598 341424 : DO kset = kset_start, nsetc
599 :
600 298384 : counter_L_blocks = counter_L_blocks + 1
601 298384 : IF (MOD(counter_L_blocks, para_env%num_pe) /= para_env%mepos) CYCLE
602 :
603 279006 : sphi_c_ext_set => sphi_c_ext(:, :, :, kset)
604 279006 : sphi_d_ext_set => sphi_d_ext(:, :, :, lset)
605 :
606 279006 : L_B_k_start = RI_index_table(katom, kset)
607 279006 : L_B_k_end = RI_index_table(katom, kset) + nsgfc(kset) - 1
608 :
609 279006 : pmax_entry = 0.0_dp
610 279006 : log10_pmax = pmax_entry
611 279006 : log10_eps_schwarz = log_zero
612 :
613 279006 : max_contraction_val = 1.0000_dp
614 :
615 1116024 : ALLOCATE (L_block(nsgfa(iset), nsgfc(kset)))
616 :
617 : CALL coulomb4(private_context, ra, rb, rc, rd, npgfa(iset), npgfb(jset), npgfc(kset), npgfd(lset), &
618 : la_min(iset), la_max(iset), lb_min(jset), lb_max(jset), &
619 : lc_min(kset), lc_max(kset), ld_min(lset), ld_max(lset), &
620 : nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
621 : sphi_a_u1, sphi_a_u2, sphi_a_u3, &
622 : sphi_b_u1, sphi_b_u2, sphi_b_u3, &
623 : sphi_c_u1, sphi_c_u2, sphi_c_u3, &
624 : sphi_d_u1, sphi_d_u2, sphi_d_u3, &
625 : zeta(1:npgfa(iset), iset), zetb(1:npgfb(jset), jset), &
626 : zetc(1:npgfc(kset), kset), zetd(1:npgfd(lset), lset), &
627 : primitive_integrals, &
628 : mp2_potential_parameter, &
629 : actual_x_data%neighbor_cells, screen_coeffs_set(1, 1, 1, 1)%x, &
630 : screen_coeffs_set(1, 1, 1, 1)%x, eps_schwarz, &
631 : max_contraction_val, cartesian_estimate, cell, neris_tmp, &
632 : log10_pmax, log10_eps_schwarz, &
633 : tmp_R_1, tmp_R_2, tmp_screen_pgf1, tmp_screen_pgf2, &
634 : pgf_list_ij, pgf_list_kl, pgf_product_list, &
635 : nsgfl_a(:, iset), nsgfl_b(:, jset), &
636 : nsgfl_c(:, kset), nsgfl_d(:, lset), &
637 : sphi_a_ext_set, &
638 : sphi_b_ext_set, &
639 : sphi_c_ext_set, &
640 : sphi_d_ext_set, &
641 : ee_work, ee_work2, ee_buffer1, ee_buffer2, ee_primitives_tmp, &
642 279006 : nimages, do_periodic)
643 :
644 279006 : primitive_counter = 0
645 558012 : DO llB = 1, nsgfd(lset)
646 1689326 : DO kkB = 1, nsgfc(kset)
647 2541634 : DO jjB = 1, nsgfb(jset)
648 4952738 : DO iiB = 1, nsgfa(iset)
649 2690110 : primitive_counter = primitive_counter + 1
650 3821424 : L_block(iiB, kkB) = primitive_integrals(primitive_counter)
651 : END DO
652 : END DO
653 : END DO
654 : END DO
655 :
656 4100430 : L_full_matrix(L_B_i_start:L_B_i_end, L_B_k_start:L_B_k_end) = L_block
657 3546782 : L_full_matrix(L_B_k_start:L_B_k_end, L_B_i_start:L_B_i_end) = TRANSPOSE(L_block)
658 :
659 338516 : DEALLOCATE (L_block)
660 :
661 : END DO ! kset
662 : END DO ! iset
663 :
664 : END DO !katom
665 : END DO !iatom
666 :
667 : ! now create a fm_matrix for the L matrix
668 2908 : CALL para_env%sum(L_full_matrix)
669 :
670 : ! create a sub blacs_env
671 2908 : NULLIFY (blacs_env)
672 2908 : CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env)
673 :
674 2908 : NULLIFY (fm_struct_L)
675 : CALL cp_fm_struct_create(fm_struct_L, context=blacs_env, nrow_global=RI_dimen, &
676 2908 : ncol_global=RI_dimen, para_env=para_env)
677 2908 : CALL cp_fm_create(fm_matrix_L, fm_struct_L, name="fm_matrix_L")
678 2908 : CALL cp_fm_struct_release(fm_struct_L)
679 2908 : CALL cp_blacs_env_release(blacs_env)
680 :
681 : CALL cp_fm_set_submatrix(fm=fm_matrix_L, new_values=L_full_matrix, start_row=1, start_col=1, &
682 2908 : n_rows=RI_dimen, n_cols=RI_dimen)
683 :
684 : info_chol = 0
685 2908 : CALL cp_fm_cholesky_decompose(matrix=fm_matrix_L, n=RI_dimen, info_out=info_chol)
686 2908 : CPASSERT(info_chol == 0)
687 :
688 : ! triangual invert
689 2908 : CALL cp_fm_triangular_invert(matrix_a=fm_matrix_L, uplo_tr='U')
690 :
691 : ! replicate L matrix to each proc
692 2908 : L_full_matrix = 0.0_dp
693 2908 : CALL cp_fm_get_submatrix(fm_matrix_L, L_full_matrix, 1, 1, RI_dimen, RI_dimen, .FALSE.)
694 2908 : CALL cp_fm_release(fm_matrix_L)
695 :
696 : ! clean lower part
697 125632 : DO iiB = 1, RI_dimen
698 2676220 : L_full_matrix(iiB + 1:RI_dimen, iiB) = 0.0_dp
699 : END DO
700 :
701 46528 : ALLOCATE (list_kl%elements(natom**2))
702 :
703 8724 : coeffs_kind_max0 = MAXVAL(screen_coeffs_kind(:, :)%x(2))
704 571176 : ALLOCATE (set_list_kl((max_set*natom)**2))
705 :
706 : !! precalculate maximum density matrix elements in blocks
707 8724 : actual_x_data%pmax_block = 0.0_dp
708 2908 : shm_pmax_block => actual_x_data%pmax_block
709 :
710 2908 : shm_atomic_pair_list => actual_x_data%atomic_pair_list
711 :
712 2908 : iatom_start = 1
713 2908 : iatom_end = natom
714 2908 : jatom_start = 1
715 2908 : jatom_end = natom
716 2908 : katom_start = 1
717 2908 : katom_end = natom
718 2908 : latom_start = 1
719 2908 : latom_end = natom
720 :
721 : CALL build_pair_list_mp2(natom, list_kl, set_list_kl, katom_start, katom_end, &
722 : latom_start, latom_end, &
723 : kind_of, basis_parameter, particle_set, &
724 : do_periodic, screen_coeffs_set, screen_coeffs_kind, &
725 : coeffs_kind_max0, log10_eps_schwarz, cell, 0.D+00, &
726 2908 : shm_atomic_pair_list)
727 :
728 2908 : virtual = dimen - occupied
729 :
730 14540 : ALLOCATE (Lai(RI_dimen, virtual, occupied))
731 2908 : Lai = 0.0_dp
732 :
733 5816 : DO iatom = 1, natom
734 2908 : jatom = iatom
735 :
736 2908 : ikind = kind_of(iatom)
737 2908 : jkind = kind_of(jatom)
738 11632 : ra = particle_set(iatom)%r(:)
739 11632 : rb = particle_set(jatom)%r(:)
740 :
741 2908 : la_max => RI_basis_parameter(ikind)%lmax
742 2908 : la_min => RI_basis_parameter(ikind)%lmin
743 2908 : npgfa => RI_basis_parameter(ikind)%npgf
744 2908 : nseta = RI_basis_parameter(ikind)%nset
745 2908 : zeta => RI_basis_parameter(ikind)%zet
746 2908 : nsgfa => RI_basis_parameter(ikind)%nsgf
747 2908 : sphi_a_ext => RI_basis_parameter(ikind)%sphi_ext(:, :, :, :)
748 2908 : nsgfl_a => RI_basis_parameter(ikind)%nsgfl
749 2908 : sphi_a_u1 = UBOUND(sphi_a_ext, 1)
750 2908 : sphi_a_u2 = UBOUND(sphi_a_ext, 2)
751 2908 : sphi_a_u3 = UBOUND(sphi_a_ext, 3)
752 :
753 2908 : lb_max => basis_S0(jkind)%lmax
754 2908 : lb_min => basis_S0(jkind)%lmin
755 2908 : npgfb => basis_S0(jkind)%npgf
756 2908 : nsetb = basis_S0(jkind)%nset
757 2908 : zetb => basis_S0(jkind)%zet
758 2908 : nsgfb => basis_S0(jkind)%nsgf
759 2908 : sphi_b_ext => basis_S0(jkind)%sphi_ext(:, :, :, :)
760 2908 : nsgfl_b => basis_S0(jkind)%nsgfl
761 2908 : sphi_b_u1 = UBOUND(sphi_b_ext, 1)
762 2908 : sphi_b_u2 = UBOUND(sphi_b_ext, 2)
763 2908 : sphi_b_u3 = UBOUND(sphi_b_ext, 3)
764 :
765 2908 : jset = 1
766 45948 : DO iset = 1, nseta
767 :
768 40132 : counter_L_blocks = counter_L_blocks + 1
769 40132 : IF (MOD(counter_L_blocks, para_env%num_pe) /= para_env%mepos) CYCLE
770 :
771 37518 : ncob = npgfb(jset)*ncoset(lb_max(jset))
772 37518 : sphi_a_ext_set => sphi_a_ext(:, :, :, iset)
773 37518 : sphi_b_ext_set => sphi_b_ext(:, :, :, jset)
774 :
775 37518 : L_B_i_start = RI_index_table(iatom, iset)
776 37518 : L_B_i_end = RI_index_table(iatom, iset) + nsgfa(iset) - 1
777 :
778 187590 : ALLOCATE (BI1(dimen, dimen, nsgfa(iset)))
779 37518 : BI1 = 0.0_dp
780 :
781 75036 : DO i_list_kl = 1, list_kl%n_element
782 :
783 37518 : katom = list_kl%elements(i_list_kl)%pair(1)
784 37518 : latom = list_kl%elements(i_list_kl)%pair(2)
785 :
786 37518 : i_set_list_kl_start = list_kl%elements(i_list_kl)%set_bounds(1)
787 37518 : i_set_list_kl_stop = list_kl%elements(i_list_kl)%set_bounds(2)
788 37518 : kkind = list_kl%elements(i_list_kl)%kind_pair(1)
789 37518 : lkind = list_kl%elements(i_list_kl)%kind_pair(2)
790 150072 : rc = list_kl%elements(i_list_kl)%r1
791 150072 : rd = list_kl%elements(i_list_kl)%r2
792 :
793 37518 : pmax_atom = 0.0_dp
794 :
795 37518 : lc_max => basis_parameter(kkind)%lmax
796 37518 : lc_min => basis_parameter(kkind)%lmin
797 37518 : npgfc => basis_parameter(kkind)%npgf
798 37518 : zetc => basis_parameter(kkind)%zet
799 37518 : nsgfc => basis_parameter(kkind)%nsgf
800 37518 : sphi_c_ext => basis_parameter(kkind)%sphi_ext(:, :, :, :)
801 37518 : nsgfl_c => basis_parameter(kkind)%nsgfl
802 37518 : sphi_c_u1 = UBOUND(sphi_c_ext, 1)
803 37518 : sphi_c_u2 = UBOUND(sphi_c_ext, 2)
804 37518 : sphi_c_u3 = UBOUND(sphi_c_ext, 3)
805 :
806 37518 : ld_max => basis_parameter(lkind)%lmax
807 37518 : ld_min => basis_parameter(lkind)%lmin
808 37518 : npgfd => basis_parameter(lkind)%npgf
809 37518 : zetd => basis_parameter(lkind)%zet
810 37518 : nsgfd => basis_parameter(lkind)%nsgf
811 37518 : sphi_d_ext => basis_parameter(lkind)%sphi_ext(:, :, :, :)
812 37518 : nsgfl_d => basis_parameter(lkind)%nsgfl
813 37518 : sphi_d_u1 = UBOUND(sphi_d_ext, 1)
814 37518 : sphi_d_u2 = UBOUND(sphi_d_ext, 2)
815 37518 : sphi_d_u3 = UBOUND(sphi_d_ext, 3)
816 :
817 412698 : DO i_set_list_kl = i_set_list_kl_start, i_set_list_kl_stop
818 337662 : kset = set_list_kl(i_set_list_kl)%pair(1)
819 337662 : lset = set_list_kl(i_set_list_kl)%pair(2)
820 :
821 337662 : IF (katom == latom .AND. lset < kset) CYCLE
822 :
823 225108 : orb_k_start = mp2_biel%index_table(katom, kset)
824 225108 : orb_k_end = orb_k_start + nsgfc(kset) - 1
825 225108 : orb_l_start = mp2_biel%index_table(latom, lset)
826 225108 : orb_l_end = orb_l_start + nsgfd(lset) - 1
827 :
828 : !! get max_vals if we screen on initial density
829 225108 : pmax_entry = 0.0_dp
830 :
831 225108 : sphi_c_ext_set => sphi_c_ext(:, :, :, kset)
832 225108 : sphi_d_ext_set => sphi_d_ext(:, :, :, lset)
833 :
834 225108 : log10_pmax = pmax_entry
835 225108 : log10_eps_schwarz = log_zero
836 :
837 225108 : IF (ALLOCATED(MNRS)) DEALLOCATE (MNRS)
838 1125540 : ALLOCATE (MNRS(nsgfd(lset), nsgfc(kset), nsgfa(iset)))
839 :
840 225108 : MNRS = 0.0_dp
841 :
842 225108 : max_contraction_val = max_contraction(kset, katom)*max_contraction(lset, latom)
843 :
844 : CALL coulomb4(private_context, ra, rb, rc, rd, npgfa(iset), npgfb(jset), npgfc(kset), npgfd(lset), &
845 : la_min(iset), la_max(iset), lb_min(jset), lb_max(jset), &
846 : lc_min(kset), lc_max(kset), ld_min(lset), ld_max(lset), &
847 : nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
848 : sphi_a_u1, sphi_a_u2, sphi_a_u3, &
849 : sphi_b_u1, sphi_b_u2, sphi_b_u3, &
850 : sphi_c_u1, sphi_c_u2, sphi_c_u3, &
851 : sphi_d_u1, sphi_d_u2, sphi_d_u3, &
852 : zeta(1:npgfa(iset), iset), zetb(1:npgfb(jset), jset), &
853 : zetc(1:npgfc(kset), kset), zetd(1:npgfd(lset), lset), &
854 : primitive_integrals, &
855 : mp2_potential_parameter, &
856 : actual_x_data%neighbor_cells, screen_coeffs_set(kset, kset, kkind, kkind)%x, &
857 : screen_coeffs_set(kset, kset, kkind, kkind)%x, eps_schwarz, &
858 : max_contraction_val, cartesian_estimate, cell, neris_tmp, &
859 : log10_pmax, log10_eps_schwarz, &
860 : tmp_R_1, tmp_R_2, tmp_screen_pgf1, tmp_screen_pgf2, &
861 : pgf_list_ij, pgf_list_kl, pgf_product_list, &
862 : nsgfl_a(:, iset), nsgfl_b(:, jset), &
863 : nsgfl_c(:, kset), nsgfl_d(:, lset), &
864 : sphi_a_ext_set, &
865 : sphi_b_ext_set, &
866 : sphi_c_ext_set, &
867 : sphi_d_ext_set, &
868 : ee_work, ee_work2, ee_buffer1, ee_buffer2, ee_primitives_tmp, &
869 225108 : nimages, do_periodic)
870 :
871 225108 : nints = nsgfa(iset)*nsgfb(jset)*nsgfc(kset)*nsgfd(lset)
872 225108 : neris_total = neris_total + nints
873 225108 : nprim_ints = nprim_ints + neris_tmp
874 225108 : IF (cartesian_estimate == 0.0_dp) cartesian_estimate = TINY(cartesian_estimate)
875 225108 : estimate_to_store_int = EXPONENT(cartesian_estimate)
876 225108 : estimate_to_store_int = MAX(estimate_to_store_int, -15_int_8)
877 225108 : cartesian_estimate = SET_EXPONENT(1.0_dp, estimate_to_store_int + 1)
878 :
879 225108 : primitive_counter = 0
880 1238094 : DO llB = 1, nsgfd(lset)
881 5477628 : DO kkB = 1, nsgfc(kset)
882 9492054 : DO jjB = 1, nsgfb(jset)
883 21444462 : DO iiB = 1, nsgfa(iset)
884 12965394 : primitive_counter = primitive_counter + 1
885 17204928 : MNRS(llB, kkB, iiB) = primitive_integrals(primitive_counter)
886 : END DO
887 : END DO
888 : END DO
889 : END DO
890 :
891 951054 : DO iiB = 1, nsgfa(iset)
892 16522272 : BI1(orb_l_start:orb_l_end, orb_k_start:orb_k_end, iiB) = MNRS(:, :, iiB)
893 17089410 : BI1(orb_k_start:orb_k_end, orb_l_start:orb_l_end, iiB) = TRANSPOSE(MNRS(:, :, iiB))
894 : END DO
895 :
896 : END DO ! i_set_list_kl
897 : END DO ! i_list_kl
898 :
899 152256 : DO iiB = 1, nsgfa(iset)
900 114738 : BI1(1:virtual, 1:occupied, iiB) = MATMUL(TRANSPOSE(C(1:dimen, occupied + 1:dimen)), &
901 216281130 : MATMUL(BI1(1:dimen, 1:dimen, iiB), C(1:dimen, 1:occupied)))
902 3823872 : Lai(L_B_i_start + iiB - 1, 1:virtual, 1:occupied) = BI1(1:virtual, 1:occupied, iiB)
903 : END DO
904 :
905 43040 : DEALLOCATE (BI1)
906 :
907 : END DO
908 : END DO
909 :
910 2908 : CALL para_env%sum(Lai)
911 :
912 11632 : DO iiB = 1, occupied
913 11632 : IF (MOD(iiB, para_env%num_pe) == para_env%mepos) THEN
914 293843910 : Lai(1:RI_dimen, 1:virtual, iiB) = MATMUL(TRANSPOSE(L_full_matrix), Lai(1:RI_dimen, 1:virtual, iiB))
915 : ELSE
916 237674 : Lai(:, :, iiB) = 0.0_dp
917 : END IF
918 : END DO
919 :
920 2908 : CALL para_env%sum(Lai)
921 :
922 2908 : DEALLOCATE (set_list_kl)
923 :
924 29080 : DO i = 1, max_pgf**2
925 26172 : DEALLOCATE (pgf_list_ij(i)%image_list)
926 29080 : DEALLOCATE (pgf_list_kl(i)%image_list)
927 : END DO
928 :
929 2908 : DEALLOCATE (pgf_list_ij)
930 2908 : DEALLOCATE (pgf_list_kl)
931 2908 : DEALLOCATE (pgf_product_list)
932 :
933 2908 : DEALLOCATE (max_contraction, kind_of)
934 :
935 2908 : DEALLOCATE (ee_work, ee_work2, ee_buffer1, ee_buffer2, ee_primitives_tmp)
936 :
937 2908 : DEALLOCATE (nimages)
938 :
939 2908 : IF (mp2_env%potential_parameter%potential_type == do_potential_TShPSC) THEN
940 0 : init_TShPSC_lmax = -1
941 0 : CALL free_C0()
942 : END IF
943 :
944 2908 : CALL timestop(handle)
945 :
946 11632 : END SUBROUTINE calc_lai_libint
947 :
948 : ! **************************************************************************************************
949 : !> \brief ...
950 : !> \param cell ...
951 : !> \param qs_env ...
952 : !> \param mp2_env ...
953 : !> \param para_env ...
954 : !> \param mp2_potential_parameter ...
955 : !> \param actual_x_data ...
956 : !> \param do_periodic ...
957 : !> \param basis_parameter ...
958 : !> \param max_set ...
959 : !> \param particle_set ...
960 : !> \param natom ...
961 : !> \param kind_of ...
962 : !> \param nsgf_max ...
963 : !> \param primitive_integrals ...
964 : !> \param ee_work ...
965 : !> \param ee_work2 ...
966 : !> \param ee_buffer1 ...
967 : !> \param ee_buffer2 ...
968 : !> \param ee_primitives_tmp ...
969 : !> \param nspins ...
970 : !> \param max_contraction ...
971 : !> \param max_pgf ...
972 : !> \param pgf_list_ij ...
973 : !> \param pgf_list_kl ...
974 : !> \param pgf_product_list ...
975 : !> \param nimages ...
976 : !> \param eps_schwarz ...
977 : !> \param log10_eps_schwarz ...
978 : !> \param private_context borrowed per-thread four-center integral context
979 : !> \param screen_coeffs_set ...
980 : !> \param screen_coeffs_kind ...
981 : !> \param screen_coeffs_pgf ...
982 : !> \param radii_pgf ...
983 : !> \param RI_basis_parameter ...
984 : !> \param RI_basis_info ...
985 : ! **************************************************************************************************
986 2951 : SUBROUTINE prepare_integral_calc(cell, qs_env, mp2_env, para_env, mp2_potential_parameter, actual_x_data, &
987 : do_periodic, basis_parameter, max_set, particle_set, natom, kind_of, &
988 : nsgf_max, primitive_integrals, ee_work, ee_work2, ee_buffer1, ee_buffer2, &
989 : ee_primitives_tmp, nspins, max_contraction, max_pgf, pgf_list_ij, &
990 : pgf_list_kl, pgf_product_list, nimages, eps_schwarz, log10_eps_schwarz, &
991 : private_context, screen_coeffs_set, screen_coeffs_kind, screen_coeffs_pgf, &
992 : radii_pgf, RI_basis_parameter, RI_basis_info)
993 :
994 : TYPE(cell_type), POINTER :: cell
995 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
996 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
997 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
998 : TYPE(hfx_potential_type), INTENT(OUT) :: mp2_potential_parameter
999 : TYPE(hfx_type), POINTER :: actual_x_data
1000 : LOGICAL, INTENT(OUT) :: do_periodic
1001 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter
1002 : INTEGER, INTENT(OUT) :: max_set
1003 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1004 : INTEGER, INTENT(OUT) :: natom
1005 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: kind_of
1006 : INTEGER, INTENT(OUT) :: nsgf_max
1007 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
1008 : INTENT(OUT) :: primitive_integrals, ee_work, ee_work2, &
1009 : ee_buffer1, ee_buffer2, &
1010 : ee_primitives_tmp
1011 : INTEGER, INTENT(OUT) :: nspins
1012 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
1013 : INTENT(OUT) :: max_contraction
1014 : INTEGER, INTENT(OUT) :: max_pgf
1015 : TYPE(hfx_pgf_list), ALLOCATABLE, DIMENSION(:), &
1016 : INTENT(OUT) :: pgf_list_ij, pgf_list_kl
1017 : TYPE(hfx_pgf_product_list), ALLOCATABLE, &
1018 : DIMENSION(:), INTENT(OUT) :: pgf_product_list
1019 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: nimages
1020 : REAL(KIND=dp), INTENT(OUT) :: eps_schwarz, log10_eps_schwarz
1021 : TYPE(four_center_integral_context_type), &
1022 : INTENT(OUT) :: private_context
1023 : TYPE(hfx_screen_coeff_type), &
1024 : DIMENSION(:, :, :, :), POINTER :: screen_coeffs_set
1025 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
1026 : POINTER :: screen_coeffs_kind
1027 : TYPE(hfx_screen_coeff_type), &
1028 : DIMENSION(:, :, :, :, :, :), POINTER :: screen_coeffs_pgf, radii_pgf
1029 : TYPE(hfx_basis_type), DIMENSION(:), OPTIONAL, &
1030 : POINTER :: RI_basis_parameter
1031 : TYPE(hfx_basis_info_type), INTENT(IN), OPTIONAL :: RI_basis_info
1032 :
1033 : CHARACTER(LEN=*), PARAMETER :: routineN = 'prepare_integral_calc'
1034 :
1035 : INTEGER :: handle, i_thread, ii, ikind, irep, &
1036 : l_max, n_threads, ncos_max, nkind, &
1037 : nneighbors
1038 2951 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1039 : TYPE(dft_control_type), POINTER :: dft_control
1040 : TYPE(hfx_basis_info_type), POINTER :: basis_info
1041 : TYPE(hfx_type), POINTER :: shm_master_x_data
1042 :
1043 2951 : CALL timeset(routineN, handle)
1044 :
1045 2951 : NULLIFY (dft_control)
1046 :
1047 2951 : irep = 1
1048 :
1049 : CALL get_qs_env(qs_env, &
1050 : atomic_kind_set=atomic_kind_set, &
1051 : cell=cell, &
1052 2951 : dft_control=dft_control)
1053 :
1054 : !! Calculate l_max used in fgamma , because init_md_ftable is definitely not thread safe
1055 2951 : nkind = SIZE(atomic_kind_set, 1)
1056 2951 : l_max = 0
1057 5920 : DO ikind = 1, nkind
1058 15003 : l_max = MAX(l_max, MAXVAL(qs_env%x_data(1, 1)%basis_parameter(ikind)%lmax))
1059 : END DO
1060 2951 : IF (PRESENT(RI_basis_parameter)) THEN
1061 5816 : DO ikind = 1, nkind
1062 45948 : l_max = MAX(l_max, MAXVAL(RI_basis_parameter(ikind)%lmax))
1063 : END DO
1064 : END IF
1065 2951 : l_max = 4*l_max
1066 2951 : CALL init_md_ftable(l_max)
1067 :
1068 2951 : CALL get_potential_parameters(mp2_env, l_max, para_env, mp2_potential_parameter)
1069 :
1070 2951 : n_threads = 1
1071 :
1072 2951 : i_thread = 0
1073 :
1074 2951 : actual_x_data => qs_env%x_data(irep, i_thread + 1)
1075 :
1076 2951 : shm_master_x_data => qs_env%x_data(irep, 1)
1077 :
1078 2951 : do_periodic = actual_x_data%periodic_parameter%do_periodic
1079 :
1080 2951 : IF (do_periodic) THEN
1081 : ! ** Rebuild neighbor lists in case the cell has changed (i.e. NPT MD)
1082 0 : actual_x_data%periodic_parameter%number_of_shells = actual_x_data%periodic_parameter%mode
1083 : CALL hfx_create_neighbor_cells(actual_x_data, actual_x_data%periodic_parameter%number_of_shells_from_input, &
1084 0 : cell, i_thread)
1085 : END IF
1086 :
1087 2951 : basis_parameter => actual_x_data%basis_parameter
1088 2951 : basis_info => actual_x_data%basis_info
1089 :
1090 2951 : max_set = basis_info%max_set
1091 2951 : IF (PRESENT(RI_basis_info)) max_set = MAX(max_set, RI_basis_info%max_set)
1092 :
1093 : CALL get_qs_env(qs_env=qs_env, &
1094 : atomic_kind_set=atomic_kind_set, &
1095 2951 : particle_set=particle_set)
1096 :
1097 2951 : natom = SIZE(particle_set, 1)
1098 :
1099 2951 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
1100 :
1101 : !! precompute maximum nco and allocate scratch
1102 2951 : ncos_max = 0
1103 2951 : nsgf_max = 0
1104 2951 : CALL get_ncos_and_ncsgf(natom, kind_of, basis_parameter, ncos_max, nsgf_max)
1105 2951 : IF (PRESENT(RI_basis_parameter)) THEN
1106 2908 : CALL get_ncos_and_ncsgf(natom, kind_of, RI_basis_parameter, ncos_max, nsgf_max)
1107 : END IF
1108 :
1109 : !! Allocate the arrays for the integrals.
1110 8853 : ALLOCATE (primitive_integrals(nsgf_max**4))
1111 2951 : primitive_integrals = 0.0_dp
1112 :
1113 8853 : ALLOCATE (ee_work(ncos_max**4))
1114 5902 : ALLOCATE (ee_work2(ncos_max**4))
1115 5902 : ALLOCATE (ee_buffer1(ncos_max**4))
1116 5902 : ALLOCATE (ee_buffer2(ncos_max**4))
1117 5902 : ALLOCATE (ee_primitives_tmp(nsgf_max**4))
1118 :
1119 2951 : nspins = dft_control%nspins
1120 :
1121 2951 : CALL get_max_contraction(max_contraction, max_set, natom, max_pgf, kind_of, basis_parameter)
1122 :
1123 : ! ** Allocate buffers for pgf_lists
1124 2951 : nneighbors = SIZE(actual_x_data%neighbor_cells)
1125 36647 : ALLOCATE (pgf_list_ij(max_pgf**2))
1126 33696 : ALLOCATE (pgf_list_kl(max_pgf**2))
1127 159354 : ALLOCATE (pgf_product_list(nneighbors**3))
1128 8853 : ALLOCATE (nimages(max_pgf**2))
1129 :
1130 30745 : DO ii = 1, max_pgf**2
1131 444704 : ALLOCATE (pgf_list_ij(ii)%image_list(nneighbors))
1132 419861 : ALLOCATE (pgf_list_kl(ii)%image_list(nneighbors))
1133 : END DO
1134 :
1135 : !!! Skipped part
1136 :
1137 : !! Get screening parameter
1138 2951 : eps_schwarz = actual_x_data%screening_parameter%eps_schwarz
1139 2951 : IF (eps_schwarz <= 0.0_dp) THEN
1140 0 : log10_eps_schwarz = log_zero
1141 : ELSE
1142 2951 : log10_eps_schwarz = LOG10(eps_schwarz)
1143 : END IF
1144 :
1145 : !! The number of integrals that fit into the given MAX_MEMORY
1146 :
1147 2951 : private_context = actual_x_data%eri4_context
1148 :
1149 : !!!!!!!! Missing part on the density matrix
1150 :
1151 : !! Initialize schwarz screening matrices for near field estimates and boxing screening matrices
1152 : !! for far field estimates. The update is only performed if the geomtry of the system changed.
1153 : !! If the system is periodic, then the corresponding routines are called and some variables
1154 : !! are initialized
1155 :
1156 2951 : IF (.NOT. shm_master_x_data%screen_funct_is_initialized) THEN
1157 : CALL calc_pair_dist_radii(qs_env, basis_parameter, &
1158 : shm_master_x_data%pair_dist_radii_pgf, max_set, max_pgf, eps_schwarz, &
1159 0 : n_threads, i_thread)
1160 : CALL calc_screening_functions(qs_env, basis_parameter, private_context, shm_master_x_data%potential_parameter, &
1161 : shm_master_x_data%screen_funct_coeffs_set, &
1162 : shm_master_x_data%screen_funct_coeffs_kind, &
1163 : shm_master_x_data%screen_funct_coeffs_pgf, &
1164 : shm_master_x_data%pair_dist_radii_pgf, &
1165 0 : max_set, max_pgf, n_threads, i_thread)
1166 :
1167 0 : shm_master_x_data%screen_funct_is_initialized = .TRUE.
1168 : END IF
1169 :
1170 2951 : screen_coeffs_set => shm_master_x_data%screen_funct_coeffs_set
1171 2951 : screen_coeffs_kind => shm_master_x_data%screen_funct_coeffs_kind
1172 2951 : screen_coeffs_pgf => shm_master_x_data%screen_funct_coeffs_pgf
1173 2951 : radii_pgf => shm_master_x_data%pair_dist_radii_pgf
1174 :
1175 2951 : CALL timestop(handle)
1176 :
1177 5902 : END SUBROUTINE prepare_integral_calc
1178 :
1179 : ! **************************************************************************************************
1180 : !> \brief ...
1181 : !> \param mp2_env ...
1182 : !> \param l_max ...
1183 : !> \param para_env ...
1184 : !> \param mp2_potential_parameter ...
1185 : ! **************************************************************************************************
1186 2951 : SUBROUTINE get_potential_parameters(mp2_env, l_max, para_env, mp2_potential_parameter)
1187 :
1188 : TYPE(mp2_type) :: mp2_env
1189 : INTEGER, INTENT(IN) :: l_max
1190 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
1191 : TYPE(hfx_potential_type), INTENT(OUT) :: mp2_potential_parameter
1192 :
1193 : INTEGER :: unit_id
1194 :
1195 2951 : IF (mp2_env%potential_parameter%potential_type == do_potential_TShPSC) THEN
1196 2 : IF (l_max > init_TShPSC_lmax) THEN
1197 2 : IF (para_env%is_source()) THEN
1198 2 : CALL open_file(unit_number=unit_id, file_name=mp2_env%potential_parameter%filename)
1199 : END IF
1200 2 : CALL init(l_max, unit_id, para_env%mepos, para_env)
1201 2 : IF (para_env%is_source()) THEN
1202 2 : CALL close_file(unit_id)
1203 : END IF
1204 2 : init_TShPSC_lmax = l_max
1205 : END IF
1206 2 : mp2_potential_parameter%cutoff_radius = mp2_env%potential_parameter%cutoff_radius/2.0_dp
1207 2949 : ELSE IF (mp2_env%potential_parameter%potential_type == do_potential_long) THEN
1208 2 : mp2_potential_parameter%omega = mp2_env%potential_parameter%omega
1209 : END IF
1210 :
1211 2951 : mp2_potential_parameter%potential_type = mp2_env%potential_parameter%potential_type
1212 :
1213 2951 : END SUBROUTINE get_potential_parameters
1214 :
1215 : ! **************************************************************************************************
1216 : !> \brief ...
1217 : !> \param natom ...
1218 : !> \param kind_of ...
1219 : !> \param basis_parameter ...
1220 : !> \param ncos_max ...
1221 : !> \param nsgf_max ...
1222 : ! **************************************************************************************************
1223 5859 : SUBROUTINE get_ncos_and_ncsgf(natom, kind_of, basis_parameter, ncos_max, nsgf_max)
1224 :
1225 : INTEGER, INTENT(IN) :: natom
1226 : INTEGER, DIMENSION(:), INTENT(IN) :: kind_of
1227 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter
1228 : INTEGER, INTENT(INOUT) :: ncos_max, nsgf_max
1229 :
1230 : INTEGER :: iatom, ikind, iset, nseta
1231 5859 : INTEGER, DIMENSION(:), POINTER :: la_max, npgfa, nsgfa
1232 :
1233 11762 : DO iatom = 1, natom
1234 5903 : ikind = kind_of(iatom)
1235 5903 : nseta = basis_parameter(ikind)%nset
1236 5903 : npgfa => basis_parameter(ikind)%npgf
1237 5903 : la_max => basis_parameter(ikind)%lmax
1238 5903 : nsgfa => basis_parameter(ikind)%nsgf
1239 61123 : DO iset = 1, nseta
1240 49361 : ncos_max = MAX(ncos_max, ncoset(la_max(iset)))
1241 55264 : nsgf_max = MAX(nsgf_max, nsgfa(iset))
1242 : END DO
1243 : END DO
1244 :
1245 5859 : END SUBROUTINE get_ncos_and_ncsgf
1246 :
1247 : ! **************************************************************************************************
1248 : !> \brief ...
1249 : !> \param max_contraction ...
1250 : !> \param max_set ...
1251 : !> \param natom ...
1252 : !> \param max_pgf ...
1253 : !> \param kind_of ...
1254 : !> \param basis_parameter ...
1255 : ! **************************************************************************************************
1256 2951 : SUBROUTINE get_max_contraction(max_contraction, max_set, natom, max_pgf, kind_of, basis_parameter)
1257 :
1258 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
1259 : INTENT(OUT) :: max_contraction
1260 : INTEGER, INTENT(IN) :: max_set, natom
1261 : INTEGER, INTENT(OUT) :: max_pgf
1262 : INTEGER, DIMENSION(:), INTENT(IN) :: kind_of
1263 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter
1264 :
1265 : INTEGER :: i, jatom, jkind, jset, ncob, nsetb, sgfb
1266 2951 : INTEGER, DIMENSION(:), POINTER :: lb_max, npgfb, nsgfb
1267 2951 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfb
1268 2951 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: sphi_b
1269 :
1270 11804 : ALLOCATE (max_contraction(max_set, natom))
1271 :
1272 2951 : max_contraction = 0.0_dp
1273 2951 : max_pgf = 0
1274 5946 : DO jatom = 1, natom
1275 2995 : jkind = kind_of(jatom)
1276 2995 : lb_max => basis_parameter(jkind)%lmax
1277 2995 : nsetb = basis_parameter(jkind)%nset
1278 2995 : npgfb => basis_parameter(jkind)%npgf
1279 2995 : first_sgfb => basis_parameter(jkind)%first_sgf
1280 2995 : sphi_b => basis_parameter(jkind)%sphi
1281 2995 : nsgfb => basis_parameter(jkind)%nsgf
1282 15175 : DO jset = 1, nsetb
1283 : ! takes the primitive to contracted transformation into account
1284 9229 : ncob = npgfb(jset)*ncoset(lb_max(jset))
1285 9229 : sgfb = first_sgfb(1, jset)
1286 : ! if the primitives are assumed to be all of max_val2, max_val2*p2s_b becomes
1287 : ! the maximum value after multiplication with sphi_b
1288 438906 : max_contraction(jset, jatom) = MAXVAL([(SUM(ABS(sphi_b(1:ncob, i))), i=sgfb, sgfb + nsgfb(jset) - 1)])
1289 12224 : max_pgf = MAX(max_pgf, npgfb(jset))
1290 : END DO
1291 : END DO
1292 :
1293 2951 : END SUBROUTINE get_max_contraction
1294 :
1295 114738 : END MODULE mp2_ri_libint
|
