Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 gamma, ONLY: init_md_ftable
39 : USE hfx_energy_potential, ONLY: coulomb4
40 : USE hfx_pair_list_methods, ONLY: build_pair_list_mp2
41 : USE hfx_screening_methods, ONLY: calc_pair_dist_radii,&
42 : calc_screening_functions
43 : USE hfx_types, ONLY: &
44 : hfx_basis_info_type, hfx_basis_type, hfx_create_neighbor_cells, hfx_pgf_list, &
45 : hfx_pgf_product_list, hfx_potential_type, hfx_screen_coeff_type, hfx_type, log_zero, &
46 : pair_set_list_type
47 : USE input_constants, ONLY: do_potential_TShPSC,&
48 : do_potential_long
49 : USE kinds, ONLY: dp,&
50 : int_8
51 : USE libint_wrapper, ONLY: cp_libint_t
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 6 : 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 :: p_work
453 2908 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: shm_pmax_block, zeta, zetb, zetc, zetd
454 2908 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: sphi_a_ext_set, sphi_b_ext_set, &
455 2908 : sphi_c_ext_set, sphi_d_ext_set
456 2908 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: sphi_a_ext, sphi_b_ext, sphi_c_ext, &
457 2908 : sphi_d_ext
458 : TYPE(cell_type), POINTER :: cell
459 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
460 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_L
461 : TYPE(cp_fm_type) :: fm_matrix_L
462 : TYPE(cp_libint_t) :: private_lib
463 2908 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter
464 2908 : TYPE(hfx_pgf_list), ALLOCATABLE, DIMENSION(:) :: pgf_list_ij, pgf_list_kl
465 : TYPE(hfx_pgf_product_list), ALLOCATABLE, &
466 2908 : DIMENSION(:) :: pgf_product_list
467 : TYPE(hfx_potential_type) :: mp2_potential_parameter
468 : TYPE(hfx_screen_coeff_type), ALLOCATABLE, &
469 2908 : DIMENSION(:, :), TARGET :: radii_pgf_large, screen_coeffs_pgf_large
470 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
471 2908 : POINTER :: screen_coeffs_kind, tmp_R_1, tmp_R_2, &
472 2908 : tmp_screen_pgf1, tmp_screen_pgf2
473 : TYPE(hfx_screen_coeff_type), &
474 2908 : DIMENSION(:, :, :, :), POINTER :: screen_coeffs_set
475 : TYPE(hfx_screen_coeff_type), &
476 2908 : DIMENSION(:, :, :, :, :, :), POINTER :: radii_pgf, screen_coeffs_pgf
477 : TYPE(hfx_type), POINTER :: actual_x_data
478 2908 : TYPE(pair_list_type_mp2) :: list_kl
479 : TYPE(pair_set_list_type), ALLOCATABLE, &
480 2908 : DIMENSION(:) :: set_list_kl
481 2908 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
482 :
483 2908 : CALL timeset(routineN, handle)
484 :
485 2908 : ln_10 = LOG(10.0_dp)
486 :
487 : !! initialize some counters
488 2908 : neris_tmp = 0_int_8
489 2908 : neris_total = 0_int_8
490 2908 : nprim_ints = 0_int_8
491 :
492 : CALL prepare_integral_calc(cell, qs_env, mp2_env, para_env, mp2_potential_parameter, actual_x_data, &
493 : do_periodic, basis_parameter, max_set, particle_set, natom, kind_of, &
494 : nsgf_max, primitive_integrals, ee_work, ee_work2, ee_buffer1, ee_buffer2, &
495 : ee_primitives_tmp, nspins, max_contraction, max_pgf, pgf_list_ij, &
496 : pgf_list_kl, pgf_product_list, nimages, eps_schwarz, log10_eps_schwarz, &
497 : private_lib, p_work, screen_coeffs_set, screen_coeffs_kind, screen_coeffs_pgf, &
498 2908 : radii_pgf, RI_basis_parameter, RI_basis_info)
499 :
500 52344 : ALLOCATE (radii_pgf_large(SIZE(radii_pgf, 1), SIZE(radii_pgf, 2)))
501 52344 : ALLOCATE (screen_coeffs_pgf_large(SIZE(screen_coeffs_pgf, 1), SIZE(screen_coeffs_pgf, 2)))
502 11632 : DO iiB = 1, SIZE(radii_pgf, 1)
503 37804 : DO jjB = 1, SIZE(radii_pgf, 2)
504 87240 : radii_pgf_large(iiB, jjB)%x = 100_dp
505 : END DO
506 : END DO
507 11632 : DO iiB = 1, SIZE(screen_coeffs_pgf, 1)
508 37804 : DO jjB = 1, SIZE(screen_coeffs_pgf, 2)
509 87240 : screen_coeffs_pgf_large(iiB, jjB)%x = 5000_dp
510 : END DO
511 : END DO
512 2908 : tmp_R_1 => radii_pgf_large(:, :)
513 2908 : tmp_R_2 => radii_pgf_large(:, :)
514 2908 : tmp_screen_pgf1 => screen_coeffs_pgf_large(:, :)
515 2908 : tmp_screen_pgf2 => screen_coeffs_pgf_large(:, :)
516 :
517 : ! start computing the L matrix
518 11632 : ALLOCATE (L_full_matrix(RI_dimen, RI_dimen))
519 5349532 : L_full_matrix = 0.0_dp
520 :
521 2908 : counter_L_blocks = 0
522 5816 : DO iatom = 1, natom
523 2908 : jatom = iatom
524 :
525 2908 : ikind = kind_of(iatom)
526 2908 : jkind = kind_of(jatom)
527 11632 : ra = particle_set(iatom)%r(:)
528 11632 : rb = particle_set(jatom)%r(:)
529 :
530 2908 : la_max => RI_basis_parameter(ikind)%lmax
531 2908 : la_min => RI_basis_parameter(ikind)%lmin
532 2908 : npgfa => RI_basis_parameter(ikind)%npgf
533 2908 : nseta = RI_basis_parameter(ikind)%nset
534 2908 : zeta => RI_basis_parameter(ikind)%zet
535 2908 : nsgfa => RI_basis_parameter(ikind)%nsgf
536 2908 : sphi_a_ext => RI_basis_parameter(ikind)%sphi_ext(:, :, :, :)
537 2908 : nsgfl_a => RI_basis_parameter(ikind)%nsgfl
538 2908 : sphi_a_u1 = UBOUND(sphi_a_ext, 1)
539 2908 : sphi_a_u2 = UBOUND(sphi_a_ext, 2)
540 2908 : sphi_a_u3 = UBOUND(sphi_a_ext, 3)
541 :
542 2908 : lb_max => basis_S0(jkind)%lmax
543 2908 : lb_min => basis_S0(jkind)%lmin
544 2908 : npgfb => basis_S0(jkind)%npgf
545 2908 : nsetb = basis_S0(jkind)%nset
546 2908 : zetb => basis_S0(jkind)%zet
547 2908 : nsgfb => basis_S0(jkind)%nsgf
548 2908 : sphi_b_ext => basis_S0(jkind)%sphi_ext(:, :, :, :)
549 2908 : nsgfl_b => basis_S0(jkind)%nsgfl
550 2908 : sphi_b_u1 = UBOUND(sphi_b_ext, 1)
551 2908 : sphi_b_u2 = UBOUND(sphi_b_ext, 2)
552 2908 : sphi_b_u3 = UBOUND(sphi_b_ext, 3)
553 :
554 8724 : DO katom = iatom, natom
555 2908 : latom = katom
556 :
557 2908 : kkind = kind_of(katom)
558 2908 : lkind = kind_of(latom)
559 11632 : rc = particle_set(katom)%r(:)
560 11632 : rd = particle_set(latom)%r(:)
561 :
562 2908 : lc_max => RI_basis_parameter(kkind)%lmax
563 2908 : lc_min => RI_basis_parameter(kkind)%lmin
564 2908 : npgfc => RI_basis_parameter(kkind)%npgf
565 2908 : nsetc = RI_basis_parameter(kkind)%nset
566 2908 : zetc => RI_basis_parameter(kkind)%zet
567 2908 : nsgfc => RI_basis_parameter(kkind)%nsgf
568 2908 : sphi_c_ext => RI_basis_parameter(kkind)%sphi_ext(:, :, :, :)
569 2908 : nsgfl_c => RI_basis_parameter(kkind)%nsgfl
570 2908 : sphi_c_u1 = UBOUND(sphi_c_ext, 1)
571 2908 : sphi_c_u2 = UBOUND(sphi_c_ext, 2)
572 2908 : sphi_c_u3 = UBOUND(sphi_c_ext, 3)
573 :
574 2908 : ld_max => basis_S0(lkind)%lmax
575 2908 : ld_min => basis_S0(lkind)%lmin
576 2908 : npgfd => basis_S0(lkind)%npgf
577 2908 : nsetd = RI_basis_parameter(lkind)%nset
578 2908 : zetd => basis_S0(lkind)%zet
579 2908 : nsgfd => basis_S0(lkind)%nsgf
580 2908 : sphi_d_ext => basis_S0(lkind)%sphi_ext(:, :, :, :)
581 2908 : nsgfl_d => basis_S0(lkind)%nsgfl
582 2908 : sphi_d_u1 = UBOUND(sphi_d_ext, 1)
583 2908 : sphi_d_u2 = UBOUND(sphi_d_ext, 2)
584 2908 : sphi_d_u3 = UBOUND(sphi_d_ext, 3)
585 :
586 2908 : jset = 1
587 2908 : lset = 1
588 45948 : DO iset = 1, nseta
589 :
590 40132 : ncob = npgfb(jset)*ncoset(lb_max(jset))
591 40132 : sphi_a_ext_set => sphi_a_ext(:, :, :, iset)
592 40132 : sphi_b_ext_set => sphi_b_ext(:, :, :, jset)
593 :
594 40132 : L_B_i_start = RI_index_table(iatom, iset)
595 40132 : L_B_i_end = RI_index_table(iatom, iset) + nsgfa(iset) - 1
596 :
597 40132 : kset_start = 1
598 40132 : IF (iatom == katom) kset_start = iset
599 341424 : DO kset = kset_start, nsetc
600 :
601 298384 : counter_L_blocks = counter_L_blocks + 1
602 298384 : IF (MOD(counter_L_blocks, para_env%num_pe) /= para_env%mepos) CYCLE
603 :
604 279006 : sphi_c_ext_set => sphi_c_ext(:, :, :, kset)
605 279006 : sphi_d_ext_set => sphi_d_ext(:, :, :, lset)
606 :
607 279006 : L_B_k_start = RI_index_table(katom, kset)
608 279006 : L_B_k_end = RI_index_table(katom, kset) + nsgfc(kset) - 1
609 :
610 279006 : pmax_entry = 0.0_dp
611 279006 : log10_pmax = pmax_entry
612 279006 : log10_eps_schwarz = log_zero
613 :
614 279006 : max_contraction_val = 1.0000_dp
615 :
616 1116024 : ALLOCATE (L_block(nsgfa(iset), nsgfc(kset)))
617 :
618 : CALL coulomb4(private_lib, ra, rb, rc, rd, npgfa(iset), npgfb(jset), npgfc(kset), npgfd(lset), &
619 : la_min(iset), la_max(iset), lb_min(jset), lb_max(jset), &
620 : lc_min(kset), lc_max(kset), ld_min(lset), ld_max(lset), &
621 : nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
622 : sphi_a_u1, sphi_a_u2, sphi_a_u3, &
623 : sphi_b_u1, sphi_b_u2, sphi_b_u3, &
624 : sphi_c_u1, sphi_c_u2, sphi_c_u3, &
625 : sphi_d_u1, sphi_d_u2, sphi_d_u3, &
626 : zeta(1:npgfa(iset), iset), zetb(1:npgfb(jset), jset), &
627 : zetc(1:npgfc(kset), kset), zetd(1:npgfd(lset), lset), &
628 : primitive_integrals, &
629 : mp2_potential_parameter, &
630 : actual_x_data%neighbor_cells, screen_coeffs_set(1, 1, 1, 1)%x, &
631 : screen_coeffs_set(1, 1, 1, 1)%x, eps_schwarz, &
632 : max_contraction_val, cartesian_estimate, cell, neris_tmp, &
633 : log10_pmax, log10_eps_schwarz, &
634 : tmp_R_1, tmp_R_2, tmp_screen_pgf1, tmp_screen_pgf2, &
635 : pgf_list_ij, pgf_list_kl, pgf_product_list, &
636 : nsgfl_a(:, iset), nsgfl_b(:, jset), &
637 : nsgfl_c(:, kset), nsgfl_d(:, lset), &
638 : sphi_a_ext_set, &
639 : sphi_b_ext_set, &
640 : sphi_c_ext_set, &
641 : sphi_d_ext_set, &
642 : ee_work, ee_work2, ee_buffer1, ee_buffer2, ee_primitives_tmp, &
643 279006 : nimages, do_periodic, p_work)
644 :
645 279006 : primitive_counter = 0
646 558012 : DO llB = 1, nsgfd(lset)
647 1689326 : DO kkB = 1, nsgfc(kset)
648 2541634 : DO jjB = 1, nsgfb(jset)
649 4952738 : DO iiB = 1, nsgfa(iset)
650 2690110 : primitive_counter = primitive_counter + 1
651 3821424 : L_block(iiB, kkB) = primitive_integrals(primitive_counter)
652 : END DO
653 : END DO
654 : END DO
655 : END DO
656 :
657 4100430 : L_full_matrix(L_B_i_start:L_B_i_end, L_B_k_start:L_B_k_end) = L_block
658 3546782 : L_full_matrix(L_B_k_start:L_B_k_end, L_B_i_start:L_B_i_end) = TRANSPOSE(L_block)
659 :
660 338516 : DEALLOCATE (L_block)
661 :
662 : END DO ! kset
663 : END DO ! iset
664 :
665 : END DO !katom
666 : END DO !iatom
667 :
668 : ! now create a fm_matrix for the L matrix
669 2908 : CALL para_env%sum(L_full_matrix)
670 :
671 : ! create a sub blacs_env
672 2908 : NULLIFY (blacs_env)
673 2908 : CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env)
674 :
675 2908 : NULLIFY (fm_struct_L)
676 : CALL cp_fm_struct_create(fm_struct_L, context=blacs_env, nrow_global=RI_dimen, &
677 2908 : ncol_global=RI_dimen, para_env=para_env)
678 2908 : CALL cp_fm_create(fm_matrix_L, fm_struct_L, name="fm_matrix_L")
679 2908 : CALL cp_fm_struct_release(fm_struct_L)
680 2908 : CALL cp_blacs_env_release(blacs_env)
681 :
682 : CALL cp_fm_set_submatrix(fm=fm_matrix_L, new_values=L_full_matrix, start_row=1, start_col=1, &
683 2908 : n_rows=RI_dimen, n_cols=RI_dimen)
684 :
685 : info_chol = 0
686 2908 : CALL cp_fm_cholesky_decompose(matrix=fm_matrix_L, n=RI_dimen, info_out=info_chol)
687 2908 : CPASSERT(info_chol == 0)
688 :
689 : ! triangual invert
690 2908 : CALL cp_fm_triangular_invert(matrix_a=fm_matrix_L, uplo_tr='U')
691 :
692 : ! replicate L matrix to each proc
693 5349532 : L_full_matrix = 0.0_dp
694 2908 : CALL cp_fm_get_submatrix(fm_matrix_L, L_full_matrix, 1, 1, RI_dimen, RI_dimen, .FALSE.)
695 2908 : CALL cp_fm_release(fm_matrix_L)
696 :
697 : ! clean lower part
698 125632 : DO iiB = 1, RI_dimen
699 2676220 : L_full_matrix(iiB + 1:RI_dimen, iiB) = 0.0_dp
700 : END DO
701 :
702 46528 : ALLOCATE (list_kl%elements(natom**2))
703 :
704 8724 : coeffs_kind_max0 = MAXVAL(screen_coeffs_kind(:, :)%x(2))
705 571176 : ALLOCATE (set_list_kl((max_set*natom)**2))
706 :
707 : !! precalculate maximum density matrix elements in blocks
708 8724 : actual_x_data%pmax_block = 0.0_dp
709 2908 : shm_pmax_block => actual_x_data%pmax_block
710 :
711 2908 : shm_atomic_pair_list => actual_x_data%atomic_pair_list
712 :
713 2908 : iatom_start = 1
714 2908 : iatom_end = natom
715 2908 : jatom_start = 1
716 2908 : jatom_end = natom
717 2908 : katom_start = 1
718 2908 : katom_end = natom
719 2908 : latom_start = 1
720 2908 : latom_end = natom
721 :
722 : CALL build_pair_list_mp2(natom, list_kl, set_list_kl, katom_start, katom_end, &
723 : latom_start, latom_end, &
724 : kind_of, basis_parameter, particle_set, &
725 : do_periodic, screen_coeffs_set, screen_coeffs_kind, &
726 : coeffs_kind_max0, log10_eps_schwarz, cell, 0.D+00, &
727 2908 : shm_atomic_pair_list)
728 :
729 2908 : virtual = dimen - occupied
730 :
731 14540 : ALLOCATE (Lai(RI_dimen, virtual, occupied))
732 3654960 : Lai = 0.0_dp
733 :
734 5816 : DO iatom = 1, natom
735 2908 : jatom = iatom
736 :
737 2908 : ikind = kind_of(iatom)
738 2908 : jkind = kind_of(jatom)
739 11632 : ra = particle_set(iatom)%r(:)
740 11632 : rb = particle_set(jatom)%r(:)
741 :
742 2908 : la_max => RI_basis_parameter(ikind)%lmax
743 2908 : la_min => RI_basis_parameter(ikind)%lmin
744 2908 : npgfa => RI_basis_parameter(ikind)%npgf
745 2908 : nseta = RI_basis_parameter(ikind)%nset
746 2908 : zeta => RI_basis_parameter(ikind)%zet
747 2908 : nsgfa => RI_basis_parameter(ikind)%nsgf
748 2908 : sphi_a_ext => RI_basis_parameter(ikind)%sphi_ext(:, :, :, :)
749 2908 : nsgfl_a => RI_basis_parameter(ikind)%nsgfl
750 2908 : sphi_a_u1 = UBOUND(sphi_a_ext, 1)
751 2908 : sphi_a_u2 = UBOUND(sphi_a_ext, 2)
752 2908 : sphi_a_u3 = UBOUND(sphi_a_ext, 3)
753 :
754 2908 : lb_max => basis_S0(jkind)%lmax
755 2908 : lb_min => basis_S0(jkind)%lmin
756 2908 : npgfb => basis_S0(jkind)%npgf
757 2908 : nsetb = basis_S0(jkind)%nset
758 2908 : zetb => basis_S0(jkind)%zet
759 2908 : nsgfb => basis_S0(jkind)%nsgf
760 2908 : sphi_b_ext => basis_S0(jkind)%sphi_ext(:, :, :, :)
761 2908 : nsgfl_b => basis_S0(jkind)%nsgfl
762 2908 : sphi_b_u1 = UBOUND(sphi_b_ext, 1)
763 2908 : sphi_b_u2 = UBOUND(sphi_b_ext, 2)
764 2908 : sphi_b_u3 = UBOUND(sphi_b_ext, 3)
765 :
766 2908 : jset = 1
767 45948 : DO iset = 1, nseta
768 :
769 40132 : counter_L_blocks = counter_L_blocks + 1
770 40132 : IF (MOD(counter_L_blocks, para_env%num_pe) /= para_env%mepos) CYCLE
771 :
772 37518 : ncob = npgfb(jset)*ncoset(lb_max(jset))
773 37518 : sphi_a_ext_set => sphi_a_ext(:, :, :, iset)
774 37518 : sphi_b_ext_set => sphi_b_ext(:, :, :, jset)
775 :
776 37518 : L_B_i_start = RI_index_table(iatom, iset)
777 37518 : L_B_i_end = RI_index_table(iatom, iset) + nsgfa(iset) - 1
778 :
779 187590 : ALLOCATE (BI1(dimen, dimen, nsgfa(iset)))
780 21034572 : BI1 = 0.0_dp
781 :
782 75036 : DO i_list_kl = 1, list_kl%n_element
783 :
784 37518 : katom = list_kl%elements(i_list_kl)%pair(1)
785 37518 : latom = list_kl%elements(i_list_kl)%pair(2)
786 :
787 37518 : i_set_list_kl_start = list_kl%elements(i_list_kl)%set_bounds(1)
788 37518 : i_set_list_kl_stop = list_kl%elements(i_list_kl)%set_bounds(2)
789 37518 : kkind = list_kl%elements(i_list_kl)%kind_pair(1)
790 37518 : lkind = list_kl%elements(i_list_kl)%kind_pair(2)
791 150072 : rc = list_kl%elements(i_list_kl)%r1
792 150072 : rd = list_kl%elements(i_list_kl)%r2
793 :
794 37518 : pmax_atom = 0.0_dp
795 :
796 37518 : lc_max => basis_parameter(kkind)%lmax
797 37518 : lc_min => basis_parameter(kkind)%lmin
798 37518 : npgfc => basis_parameter(kkind)%npgf
799 37518 : zetc => basis_parameter(kkind)%zet
800 37518 : nsgfc => basis_parameter(kkind)%nsgf
801 37518 : sphi_c_ext => basis_parameter(kkind)%sphi_ext(:, :, :, :)
802 37518 : nsgfl_c => basis_parameter(kkind)%nsgfl
803 37518 : sphi_c_u1 = UBOUND(sphi_c_ext, 1)
804 37518 : sphi_c_u2 = UBOUND(sphi_c_ext, 2)
805 37518 : sphi_c_u3 = UBOUND(sphi_c_ext, 3)
806 :
807 37518 : ld_max => basis_parameter(lkind)%lmax
808 37518 : ld_min => basis_parameter(lkind)%lmin
809 37518 : npgfd => basis_parameter(lkind)%npgf
810 37518 : zetd => basis_parameter(lkind)%zet
811 37518 : nsgfd => basis_parameter(lkind)%nsgf
812 37518 : sphi_d_ext => basis_parameter(lkind)%sphi_ext(:, :, :, :)
813 37518 : nsgfl_d => basis_parameter(lkind)%nsgfl
814 37518 : sphi_d_u1 = UBOUND(sphi_d_ext, 1)
815 37518 : sphi_d_u2 = UBOUND(sphi_d_ext, 2)
816 37518 : sphi_d_u3 = UBOUND(sphi_d_ext, 3)
817 :
818 412698 : DO i_set_list_kl = i_set_list_kl_start, i_set_list_kl_stop
819 337662 : kset = set_list_kl(i_set_list_kl)%pair(1)
820 337662 : lset = set_list_kl(i_set_list_kl)%pair(2)
821 :
822 337662 : IF (katom == latom .AND. lset < kset) CYCLE
823 :
824 225108 : orb_k_start = mp2_biel%index_table(katom, kset)
825 225108 : orb_k_end = orb_k_start + nsgfc(kset) - 1
826 225108 : orb_l_start = mp2_biel%index_table(latom, lset)
827 225108 : orb_l_end = orb_l_start + nsgfd(lset) - 1
828 :
829 : !! get max_vals if we screen on initial density
830 225108 : pmax_entry = 0.0_dp
831 :
832 225108 : sphi_c_ext_set => sphi_c_ext(:, :, :, kset)
833 225108 : sphi_d_ext_set => sphi_d_ext(:, :, :, lset)
834 :
835 225108 : log10_pmax = pmax_entry
836 225108 : log10_eps_schwarz = log_zero
837 :
838 225108 : IF (ALLOCATED(MNRS)) DEALLOCATE (MNRS)
839 1125540 : ALLOCATE (MNRS(nsgfd(lset), nsgfc(kset), nsgfa(iset)))
840 :
841 16747380 : MNRS = 0.0_dp
842 :
843 225108 : max_contraction_val = max_contraction(kset, katom)*max_contraction(lset, latom)
844 :
845 : CALL coulomb4(private_lib, ra, rb, rc, rd, npgfa(iset), npgfb(jset), npgfc(kset), npgfd(lset), &
846 : la_min(iset), la_max(iset), lb_min(jset), lb_max(jset), &
847 : lc_min(kset), lc_max(kset), ld_min(lset), ld_max(lset), &
848 : nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
849 : sphi_a_u1, sphi_a_u2, sphi_a_u3, &
850 : sphi_b_u1, sphi_b_u2, sphi_b_u3, &
851 : sphi_c_u1, sphi_c_u2, sphi_c_u3, &
852 : sphi_d_u1, sphi_d_u2, sphi_d_u3, &
853 : zeta(1:npgfa(iset), iset), zetb(1:npgfb(jset), jset), &
854 : zetc(1:npgfc(kset), kset), zetd(1:npgfd(lset), lset), &
855 : primitive_integrals, &
856 : mp2_potential_parameter, &
857 : actual_x_data%neighbor_cells, screen_coeffs_set(kset, kset, kkind, kkind)%x, &
858 : screen_coeffs_set(kset, kset, kkind, kkind)%x, eps_schwarz, &
859 : max_contraction_val, cartesian_estimate, cell, neris_tmp, &
860 : log10_pmax, log10_eps_schwarz, &
861 : tmp_R_1, tmp_R_2, tmp_screen_pgf1, tmp_screen_pgf2, &
862 : pgf_list_ij, pgf_list_kl, pgf_product_list, &
863 : nsgfl_a(:, iset), nsgfl_b(:, jset), &
864 : nsgfl_c(:, kset), nsgfl_d(:, lset), &
865 : sphi_a_ext_set, &
866 : sphi_b_ext_set, &
867 : sphi_c_ext_set, &
868 : sphi_d_ext_set, &
869 : ee_work, ee_work2, ee_buffer1, ee_buffer2, ee_primitives_tmp, &
870 225108 : nimages, do_periodic, p_work)
871 :
872 225108 : nints = nsgfa(iset)*nsgfb(jset)*nsgfc(kset)*nsgfd(lset)
873 225108 : neris_total = neris_total + nints
874 225108 : nprim_ints = nprim_ints + neris_tmp
875 225108 : IF (cartesian_estimate == 0.0_dp) cartesian_estimate = TINY(cartesian_estimate)
876 225108 : estimate_to_store_int = EXPONENT(cartesian_estimate)
877 225108 : estimate_to_store_int = MAX(estimate_to_store_int, -15_int_8)
878 225108 : cartesian_estimate = SET_EXPONENT(1.0_dp, estimate_to_store_int + 1)
879 :
880 225108 : primitive_counter = 0
881 1238094 : DO llB = 1, nsgfd(lset)
882 5477628 : DO kkB = 1, nsgfc(kset)
883 9492054 : DO jjB = 1, nsgfb(jset)
884 21444462 : DO iiB = 1, nsgfa(iset)
885 12965394 : primitive_counter = primitive_counter + 1
886 17204928 : MNRS(llB, kkB, iiB) = primitive_integrals(primitive_counter)
887 : END DO
888 : END DO
889 : END DO
890 : END DO
891 :
892 951054 : DO iiB = 1, nsgfa(iset)
893 16522272 : BI1(orb_l_start:orb_l_end, orb_k_start:orb_k_end, iiB) = MNRS(:, :, iiB)
894 17089410 : BI1(orb_k_start:orb_k_end, orb_l_start:orb_l_end, iiB) = TRANSPOSE(MNRS(:, :, iiB))
895 : END DO
896 :
897 : END DO ! i_set_list_kl
898 : END DO ! i_list_kl
899 :
900 152256 : DO iiB = 1, nsgfa(iset)
901 114738 : BI1(1:virtual, 1:occupied, iiB) = MATMUL(TRANSPOSE(C(1:dimen, occupied + 1:dimen)), &
902 221100126 : MATMUL(BI1(1:dimen, 1:dimen, iiB), C(1:dimen, 1:occupied)))
903 3823872 : Lai(L_B_i_start + iiB - 1, 1:virtual, 1:occupied) = BI1(1:virtual, 1:occupied, iiB)
904 : END DO
905 :
906 43040 : DEALLOCATE (BI1)
907 :
908 : END DO
909 : END DO
910 :
911 2908 : CALL para_env%sum(Lai)
912 :
913 11632 : DO iiB = 1, occupied
914 11632 : IF (MOD(iiB, para_env%num_pe) == para_env%mepos) THEN
915 293843910 : Lai(1:RI_dimen, 1:virtual, iiB) = MATMUL(TRANSPOSE(L_full_matrix), Lai(1:RI_dimen, 1:virtual, iiB))
916 : ELSE
917 237674 : Lai(:, :, iiB) = 0.0_dp
918 : END IF
919 : END DO
920 :
921 2908 : CALL para_env%sum(Lai)
922 :
923 2908 : DEALLOCATE (set_list_kl)
924 :
925 29080 : DO i = 1, max_pgf**2
926 26172 : DEALLOCATE (pgf_list_ij(i)%image_list)
927 29080 : DEALLOCATE (pgf_list_kl(i)%image_list)
928 : END DO
929 :
930 2908 : DEALLOCATE (pgf_list_ij)
931 2908 : DEALLOCATE (pgf_list_kl)
932 2908 : DEALLOCATE (pgf_product_list)
933 :
934 2908 : DEALLOCATE (max_contraction, kind_of)
935 :
936 2908 : DEALLOCATE (ee_work, ee_work2, ee_buffer1, ee_buffer2, ee_primitives_tmp)
937 :
938 2908 : DEALLOCATE (nimages)
939 :
940 2908 : IF (mp2_env%potential_parameter%potential_type == do_potential_TShPSC) THEN
941 0 : init_TShPSC_lmax = -1
942 0 : CALL free_C0()
943 : END IF
944 :
945 2908 : CALL timestop(handle)
946 :
947 14540 : END SUBROUTINE calc_lai_libint
948 :
949 : ! **************************************************************************************************
950 : !> \brief ...
951 : !> \param cell ...
952 : !> \param qs_env ...
953 : !> \param mp2_env ...
954 : !> \param para_env ...
955 : !> \param mp2_potential_parameter ...
956 : !> \param actual_x_data ...
957 : !> \param do_periodic ...
958 : !> \param basis_parameter ...
959 : !> \param max_set ...
960 : !> \param particle_set ...
961 : !> \param natom ...
962 : !> \param kind_of ...
963 : !> \param nsgf_max ...
964 : !> \param primitive_integrals ...
965 : !> \param ee_work ...
966 : !> \param ee_work2 ...
967 : !> \param ee_buffer1 ...
968 : !> \param ee_buffer2 ...
969 : !> \param ee_primitives_tmp ...
970 : !> \param nspins ...
971 : !> \param max_contraction ...
972 : !> \param max_pgf ...
973 : !> \param pgf_list_ij ...
974 : !> \param pgf_list_kl ...
975 : !> \param pgf_product_list ...
976 : !> \param nimages ...
977 : !> \param eps_schwarz ...
978 : !> \param log10_eps_schwarz ...
979 : !> \param private_lib ...
980 : !> \param p_work ...
981 : !> \param screen_coeffs_set ...
982 : !> \param screen_coeffs_kind ...
983 : !> \param screen_coeffs_pgf ...
984 : !> \param radii_pgf ...
985 : !> \param RI_basis_parameter ...
986 : !> \param RI_basis_info ...
987 : ! **************************************************************************************************
988 2951 : SUBROUTINE prepare_integral_calc(cell, qs_env, mp2_env, para_env, mp2_potential_parameter, actual_x_data, &
989 : do_periodic, basis_parameter, max_set, particle_set, natom, kind_of, &
990 : nsgf_max, primitive_integrals, ee_work, ee_work2, ee_buffer1, ee_buffer2, &
991 : ee_primitives_tmp, nspins, max_contraction, max_pgf, pgf_list_ij, &
992 : pgf_list_kl, pgf_product_list, nimages, eps_schwarz, log10_eps_schwarz, &
993 : private_lib, p_work, screen_coeffs_set, screen_coeffs_kind, screen_coeffs_pgf, &
994 : radii_pgf, RI_basis_parameter, RI_basis_info)
995 :
996 : TYPE(cell_type), POINTER :: cell
997 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
998 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
999 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
1000 : TYPE(hfx_potential_type), INTENT(OUT) :: mp2_potential_parameter
1001 : TYPE(hfx_type), POINTER :: actual_x_data
1002 : LOGICAL, INTENT(OUT) :: do_periodic
1003 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter
1004 : INTEGER, INTENT(OUT) :: max_set
1005 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1006 : INTEGER, INTENT(OUT) :: natom
1007 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: kind_of
1008 : INTEGER, INTENT(OUT) :: nsgf_max
1009 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
1010 : INTENT(OUT) :: primitive_integrals, ee_work, ee_work2, &
1011 : ee_buffer1, ee_buffer2, &
1012 : ee_primitives_tmp
1013 : INTEGER, INTENT(OUT) :: nspins
1014 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
1015 : INTENT(OUT) :: max_contraction
1016 : INTEGER, INTENT(OUT) :: max_pgf
1017 : TYPE(hfx_pgf_list), ALLOCATABLE, DIMENSION(:), &
1018 : INTENT(OUT) :: pgf_list_ij, pgf_list_kl
1019 : TYPE(hfx_pgf_product_list), ALLOCATABLE, &
1020 : DIMENSION(:), INTENT(OUT) :: pgf_product_list
1021 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: nimages
1022 : REAL(KIND=dp), INTENT(OUT) :: eps_schwarz, log10_eps_schwarz
1023 : TYPE(cp_libint_t), INTENT(OUT) :: private_lib
1024 : REAL(KIND=dp), DIMENSION(:), POINTER :: p_work
1025 : TYPE(hfx_screen_coeff_type), &
1026 : DIMENSION(:, :, :, :), POINTER :: screen_coeffs_set
1027 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
1028 : POINTER :: screen_coeffs_kind
1029 : TYPE(hfx_screen_coeff_type), &
1030 : DIMENSION(:, :, :, :, :, :), POINTER :: screen_coeffs_pgf, radii_pgf
1031 : TYPE(hfx_basis_type), DIMENSION(:), OPTIONAL, &
1032 : POINTER :: RI_basis_parameter
1033 : TYPE(hfx_basis_info_type), INTENT(IN), OPTIONAL :: RI_basis_info
1034 :
1035 : CHARACTER(LEN=*), PARAMETER :: routineN = 'prepare_integral_calc'
1036 :
1037 : INTEGER :: handle, i_thread, ii, ikind, irep, &
1038 : l_max, n_threads, ncos_max, nkind, &
1039 : nneighbors
1040 2951 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1041 : TYPE(dft_control_type), POINTER :: dft_control
1042 : TYPE(hfx_basis_info_type), POINTER :: basis_info
1043 : TYPE(hfx_type), POINTER :: shm_master_x_data
1044 :
1045 2951 : CALL timeset(routineN, handle)
1046 :
1047 2951 : NULLIFY (dft_control)
1048 :
1049 2951 : irep = 1
1050 :
1051 : CALL get_qs_env(qs_env, &
1052 : atomic_kind_set=atomic_kind_set, &
1053 : cell=cell, &
1054 2951 : dft_control=dft_control)
1055 :
1056 : !! Calculate l_max used in fgamma , because init_md_ftable is definitely not thread safe
1057 2951 : nkind = SIZE(atomic_kind_set, 1)
1058 2951 : l_max = 0
1059 5920 : DO ikind = 1, nkind
1060 15003 : l_max = MAX(l_max, MAXVAL(qs_env%x_data(1, 1)%basis_parameter(ikind)%lmax))
1061 : END DO
1062 2951 : IF (PRESENT(RI_basis_parameter)) THEN
1063 5816 : DO ikind = 1, nkind
1064 45948 : l_max = MAX(l_max, MAXVAL(RI_basis_parameter(ikind)%lmax))
1065 : END DO
1066 : END IF
1067 2951 : l_max = 4*l_max
1068 2951 : CALL init_md_ftable(l_max)
1069 :
1070 2951 : CALL get_potential_parameters(mp2_env, l_max, para_env, mp2_potential_parameter)
1071 :
1072 2951 : n_threads = 1
1073 :
1074 2951 : i_thread = 0
1075 :
1076 2951 : actual_x_data => qs_env%x_data(irep, i_thread + 1)
1077 :
1078 2951 : shm_master_x_data => qs_env%x_data(irep, 1)
1079 :
1080 2951 : do_periodic = actual_x_data%periodic_parameter%do_periodic
1081 :
1082 2951 : IF (do_periodic) THEN
1083 : ! ** Rebuild neighbor lists in case the cell has changed (i.e. NPT MD)
1084 0 : actual_x_data%periodic_parameter%number_of_shells = actual_x_data%periodic_parameter%mode
1085 : CALL hfx_create_neighbor_cells(actual_x_data, actual_x_data%periodic_parameter%number_of_shells_from_input, &
1086 0 : cell, i_thread)
1087 : END IF
1088 :
1089 2951 : basis_parameter => actual_x_data%basis_parameter
1090 2951 : basis_info => actual_x_data%basis_info
1091 :
1092 2951 : max_set = basis_info%max_set
1093 2951 : IF (PRESENT(RI_basis_info)) max_set = MAX(max_set, RI_basis_info%max_set)
1094 :
1095 : CALL get_qs_env(qs_env=qs_env, &
1096 : atomic_kind_set=atomic_kind_set, &
1097 2951 : particle_set=particle_set)
1098 :
1099 2951 : natom = SIZE(particle_set, 1)
1100 :
1101 2951 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
1102 :
1103 : !! precompute maximum nco and allocate scratch
1104 2951 : ncos_max = 0
1105 2951 : nsgf_max = 0
1106 2951 : CALL get_ncos_and_ncsgf(natom, kind_of, basis_parameter, ncos_max, nsgf_max)
1107 2951 : IF (PRESENT(RI_basis_parameter)) THEN
1108 2908 : CALL get_ncos_and_ncsgf(natom, kind_of, RI_basis_parameter, ncos_max, nsgf_max)
1109 : END IF
1110 :
1111 : !! Allocate the arrays for the integrals.
1112 8853 : ALLOCATE (primitive_integrals(nsgf_max**4))
1113 7024366 : primitive_integrals = 0.0_dp
1114 :
1115 8853 : ALLOCATE (ee_work(ncos_max**4))
1116 5902 : ALLOCATE (ee_work2(ncos_max**4))
1117 5902 : ALLOCATE (ee_buffer1(ncos_max**4))
1118 5902 : ALLOCATE (ee_buffer2(ncos_max**4))
1119 5902 : ALLOCATE (ee_primitives_tmp(nsgf_max**4))
1120 :
1121 2951 : nspins = dft_control%nspins
1122 :
1123 2951 : CALL get_max_contraction(max_contraction, max_set, natom, max_pgf, kind_of, basis_parameter)
1124 :
1125 : ! ** Allocate buffers for pgf_lists
1126 2951 : nneighbors = SIZE(actual_x_data%neighbor_cells)
1127 36647 : ALLOCATE (pgf_list_ij(max_pgf**2))
1128 33696 : ALLOCATE (pgf_list_kl(max_pgf**2))
1129 153452 : ALLOCATE (pgf_product_list(nneighbors**3))
1130 8853 : ALLOCATE (nimages(max_pgf**2))
1131 :
1132 30745 : DO ii = 1, max_pgf**2
1133 444704 : ALLOCATE (pgf_list_ij(ii)%image_list(nneighbors))
1134 419861 : ALLOCATE (pgf_list_kl(ii)%image_list(nneighbors))
1135 : END DO
1136 :
1137 : !!! Skipped part
1138 :
1139 : !! Get screening parameter
1140 2951 : eps_schwarz = actual_x_data%screening_parameter%eps_schwarz
1141 2951 : IF (eps_schwarz <= 0.0_dp) THEN
1142 0 : log10_eps_schwarz = log_zero
1143 : ELSE
1144 2951 : log10_eps_schwarz = LOG10(eps_schwarz)
1145 : END IF
1146 :
1147 : !! The number of integrals that fit into the given MAX_MEMORY
1148 :
1149 2951 : private_lib = actual_x_data%lib
1150 :
1151 : !!!!!!!! Missing part on the density matrix
1152 :
1153 : !! Initialize schwarz screening matrices for near field estimates and boxing screening matrices
1154 : !! for far field estimates. The update is only performed if the geomtry of the system changed.
1155 : !! If the system is periodic, then the corresponding routines are called and some variables
1156 : !! are initialized
1157 :
1158 2951 : IF (.NOT. shm_master_x_data%screen_funct_is_initialized) THEN
1159 : CALL calc_pair_dist_radii(qs_env, basis_parameter, &
1160 : shm_master_x_data%pair_dist_radii_pgf, max_set, max_pgf, eps_schwarz, &
1161 0 : n_threads, i_thread)
1162 : CALL calc_screening_functions(qs_env, basis_parameter, private_lib, shm_master_x_data%potential_parameter, &
1163 : shm_master_x_data%screen_funct_coeffs_set, &
1164 : shm_master_x_data%screen_funct_coeffs_kind, &
1165 : shm_master_x_data%screen_funct_coeffs_pgf, &
1166 : shm_master_x_data%pair_dist_radii_pgf, &
1167 0 : max_set, max_pgf, n_threads, i_thread, p_work)
1168 :
1169 0 : shm_master_x_data%screen_funct_is_initialized = .TRUE.
1170 : END IF
1171 :
1172 2951 : screen_coeffs_set => shm_master_x_data%screen_funct_coeffs_set
1173 2951 : screen_coeffs_kind => shm_master_x_data%screen_funct_coeffs_kind
1174 2951 : screen_coeffs_pgf => shm_master_x_data%screen_funct_coeffs_pgf
1175 2951 : radii_pgf => shm_master_x_data%pair_dist_radii_pgf
1176 :
1177 2951 : CALL timestop(handle)
1178 :
1179 8853 : END SUBROUTINE prepare_integral_calc
1180 :
1181 : ! **************************************************************************************************
1182 : !> \brief ...
1183 : !> \param mp2_env ...
1184 : !> \param l_max ...
1185 : !> \param para_env ...
1186 : !> \param mp2_potential_parameter ...
1187 : ! **************************************************************************************************
1188 2951 : SUBROUTINE get_potential_parameters(mp2_env, l_max, para_env, mp2_potential_parameter)
1189 :
1190 : TYPE(mp2_type) :: mp2_env
1191 : INTEGER, INTENT(IN) :: l_max
1192 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
1193 : TYPE(hfx_potential_type), INTENT(OUT) :: mp2_potential_parameter
1194 :
1195 : INTEGER :: unit_id
1196 :
1197 2951 : IF (mp2_env%potential_parameter%potential_type == do_potential_TShPSC) THEN
1198 2 : IF (l_max > init_TShPSC_lmax) THEN
1199 2 : IF (para_env%is_source()) THEN
1200 2 : CALL open_file(unit_number=unit_id, file_name=mp2_env%potential_parameter%filename)
1201 : END IF
1202 2 : CALL init(l_max, unit_id, para_env%mepos, para_env)
1203 2 : IF (para_env%is_source()) THEN
1204 2 : CALL close_file(unit_id)
1205 : END IF
1206 2 : init_TShPSC_lmax = l_max
1207 : END IF
1208 2 : mp2_potential_parameter%cutoff_radius = mp2_env%potential_parameter%cutoff_radius/2.0_dp
1209 2949 : ELSE IF (mp2_env%potential_parameter%potential_type == do_potential_long) THEN
1210 2 : mp2_potential_parameter%omega = mp2_env%potential_parameter%omega
1211 : END IF
1212 :
1213 2951 : mp2_potential_parameter%potential_type = mp2_env%potential_parameter%potential_type
1214 :
1215 2951 : END SUBROUTINE get_potential_parameters
1216 :
1217 : ! **************************************************************************************************
1218 : !> \brief ...
1219 : !> \param natom ...
1220 : !> \param kind_of ...
1221 : !> \param basis_parameter ...
1222 : !> \param ncos_max ...
1223 : !> \param nsgf_max ...
1224 : ! **************************************************************************************************
1225 5859 : SUBROUTINE get_ncos_and_ncsgf(natom, kind_of, basis_parameter, ncos_max, nsgf_max)
1226 :
1227 : INTEGER, INTENT(IN) :: natom
1228 : INTEGER, DIMENSION(:), INTENT(IN) :: kind_of
1229 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter
1230 : INTEGER, INTENT(INOUT) :: ncos_max, nsgf_max
1231 :
1232 : INTEGER :: iatom, ikind, iset, nseta
1233 5859 : INTEGER, DIMENSION(:), POINTER :: la_max, npgfa, nsgfa
1234 :
1235 11762 : DO iatom = 1, natom
1236 5903 : ikind = kind_of(iatom)
1237 5903 : nseta = basis_parameter(ikind)%nset
1238 5903 : npgfa => basis_parameter(ikind)%npgf
1239 5903 : la_max => basis_parameter(ikind)%lmax
1240 5903 : nsgfa => basis_parameter(ikind)%nsgf
1241 61123 : DO iset = 1, nseta
1242 49361 : ncos_max = MAX(ncos_max, ncoset(la_max(iset)))
1243 55264 : nsgf_max = MAX(nsgf_max, nsgfa(iset))
1244 : END DO
1245 : END DO
1246 :
1247 5859 : END SUBROUTINE get_ncos_and_ncsgf
1248 :
1249 : ! **************************************************************************************************
1250 : !> \brief ...
1251 : !> \param max_contraction ...
1252 : !> \param max_set ...
1253 : !> \param natom ...
1254 : !> \param max_pgf ...
1255 : !> \param kind_of ...
1256 : !> \param basis_parameter ...
1257 : ! **************************************************************************************************
1258 2951 : SUBROUTINE get_max_contraction(max_contraction, max_set, natom, max_pgf, kind_of, basis_parameter)
1259 :
1260 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
1261 : INTENT(OUT) :: max_contraction
1262 : INTEGER, INTENT(IN) :: max_set, natom
1263 : INTEGER, INTENT(OUT) :: max_pgf
1264 : INTEGER, DIMENSION(:), INTENT(IN) :: kind_of
1265 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter
1266 :
1267 : INTEGER :: i, jatom, jkind, jset, ncob, nsetb, sgfb
1268 2951 : INTEGER, DIMENSION(:), POINTER :: lb_max, npgfb, nsgfb
1269 2951 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfb
1270 2951 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: sphi_b
1271 :
1272 11804 : ALLOCATE (max_contraction(max_set, natom))
1273 :
1274 46709 : max_contraction = 0.0_dp
1275 2951 : max_pgf = 0
1276 5946 : DO jatom = 1, natom
1277 2995 : jkind = kind_of(jatom)
1278 2995 : lb_max => basis_parameter(jkind)%lmax
1279 2995 : nsetb = basis_parameter(jkind)%nset
1280 2995 : npgfb => basis_parameter(jkind)%npgf
1281 2995 : first_sgfb => basis_parameter(jkind)%first_sgf
1282 2995 : sphi_b => basis_parameter(jkind)%sphi
1283 2995 : nsgfb => basis_parameter(jkind)%nsgf
1284 15175 : DO jset = 1, nsetb
1285 : ! takes the primitive to contracted transformation into account
1286 9229 : ncob = npgfb(jset)*ncoset(lb_max(jset))
1287 9229 : sgfb = first_sgfb(1, jset)
1288 : ! if the primitives are assumed to be all of max_val2, max_val2*p2s_b becomes
1289 : ! the maximum value after multiplication with sphi_b
1290 448135 : max_contraction(jset, jatom) = MAXVAL((/(SUM(ABS(sphi_b(1:ncob, i))), i=sgfb, sgfb + nsgfb(jset) - 1)/))
1291 12224 : max_pgf = MAX(max_pgf, npgfb(jset))
1292 : END DO
1293 : END DO
1294 :
1295 2951 : END SUBROUTINE get_max_contraction
1296 :
1297 114738 : END MODULE mp2_ri_libint
|
