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 3-center integrals machinery for the XAS_TDP method
10 : !> \author A. Bussy (03.2020)
11 : ! **************************************************************************************************
12 :
13 : MODULE xas_tdp_integrals
14 : USE OMP_LIB, ONLY: omp_get_max_threads, &
15 : omp_get_num_threads, &
16 : omp_get_thread_num
17 : USE ai_contraction_sphi, ONLY: ab_contract, &
18 : abc_contract_xsmm
19 : USE atomic_kind_types, ONLY: atomic_kind_type
20 : USE basis_set_types, ONLY: get_gto_basis_set, &
21 : gto_basis_set_p_type, &
22 : gto_basis_set_type
23 : USE cell_types, ONLY: cell_type
24 : USE constants_operator, ONLY: operator_coulomb
25 : USE cp_array_utils, ONLY: cp_1d_i_p_type, &
26 : cp_2d_r_p_type
27 : USE cp_blacs_env, ONLY: cp_blacs_env_type
28 : USE cp_dbcsr_api, ONLY: dbcsr_distribution_get, &
29 : dbcsr_distribution_release, &
30 : dbcsr_distribution_type
31 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_dist2d_to_dist
32 : USE cp_eri_mme_interface, ONLY: cp_eri_mme_param, &
33 : cp_eri_mme_set_params
34 : USE cp_files, ONLY: close_file, &
35 : open_file
36 : USE dbt_api, ONLY: &
37 : dbt_create, dbt_distribution_destroy, dbt_distribution_new, dbt_distribution_type, &
38 : dbt_finalize, dbt_pgrid_create, dbt_pgrid_destroy, dbt_pgrid_type, dbt_put_block, &
39 : dbt_reserve_blocks, dbt_type
40 : USE distribution_1d_types, ONLY: distribution_1d_type
41 : USE distribution_2d_types, ONLY: distribution_2d_create, &
42 : distribution_2d_type
43 : USE eri_mme_integrate, ONLY: eri_mme_2c_integrate
44 : USE eri_mme_types, ONLY: eri_mme_init, &
45 : eri_mme_release
46 : USE gamma, ONLY: init_md_ftable
47 : USE generic_os_integrals, ONLY: int_operators_r12_ab_os
48 : USE input_constants, ONLY: do_potential_coulomb, &
49 : do_potential_id, &
50 : do_potential_short, &
51 : do_potential_truncated
52 : USE input_section_types, ONLY: section_vals_val_get
53 : USE kinds, ONLY: dp
54 : USE libint_2c_3c, ONLY: cutoff_screen_factor, &
55 : eri_2center, &
56 : eri_3center, &
57 : libint_potential_type
58 : USE libint_wrapper, ONLY: cp_libint_cleanup_2eri, &
59 : cp_libint_cleanup_3eri, &
60 : cp_libint_init_2eri, &
61 : cp_libint_init_3eri, &
62 : cp_libint_set_contrdepth, &
63 : cp_libint_t
64 : USE mathlib, ONLY: invmat_symm
65 : USE message_passing, ONLY: mp_comm_type, &
66 : mp_para_env_type
67 : USE molecule_types, ONLY: molecule_type
68 : USE orbital_pointers, ONLY: ncoset
69 : USE particle_methods, ONLY: get_particle_set
70 : USE particle_types, ONLY: particle_type
71 : USE qs_environment_types, ONLY: get_qs_env, &
72 : qs_environment_type
73 : USE qs_integral_utils, ONLY: basis_set_list_setup
74 : USE qs_kind_types, ONLY: get_qs_kind, &
75 : qs_kind_type
76 : USE qs_neighbor_list_types, ONLY: &
77 : get_iterator_info, get_neighbor_list_set_p, neighbor_list_iterate, &
78 : neighbor_list_iterator_create, neighbor_list_iterator_p_type, &
79 : neighbor_list_iterator_release, neighbor_list_set_p_type, nl_set_sub_iterator, &
80 : nl_sub_iterate, release_neighbor_list_sets
81 : USE qs_neighbor_lists, ONLY: atom2d_build, &
82 : atom2d_cleanup, &
83 : build_neighbor_lists, &
84 : local_atoms_type, &
85 : pair_radius_setup
86 : USE qs_o3c_types, ONLY: get_o3c_iterator_info, &
87 : init_o3c_container, &
88 : o3c_container_type, &
89 : o3c_iterate, &
90 : o3c_iterator_create, &
91 : o3c_iterator_release, &
92 : o3c_iterator_type, &
93 : release_o3c_container
94 : USE t_c_g0, ONLY: get_lmax_init, &
95 : init
96 : USE xas_tdp_types, ONLY: xas_tdp_control_type, &
97 : xas_tdp_env_type
98 : #if defined(__LIBXS)
99 : USE ISO_C_BINDING, ONLY: C_LOC
100 : USE libxs, ONLY: libxs_matcopy, &
101 : libxs_otrans
102 : #endif
103 : #include "./base/base_uses.f90"
104 :
105 : IMPLICIT NONE
106 : PRIVATE
107 :
108 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xas_tdp_integrals'
109 :
110 : PUBLIC :: create_pqX_tensor, fill_pqX_tensor, compute_ri_3c_coulomb, compute_ri_3c_exchange, &
111 : build_xas_tdp_3c_nl, build_xas_tdp_ovlp_nl, get_opt_3c_dist2d, &
112 : compute_ri_coulomb2_int, compute_ri_exchange2_int
113 :
114 : CONTAINS
115 :
116 : ! **************************************************************************************************
117 : !> \brief Prepares a tensor to hold 3-center integrals (pq|X), where p,q are distributed according
118 : !> to the given 2d dbcsr distribution of the given . The third dimension of the tensor is
119 : !> iteslf not distributed (i.e. the t_pgrid's third dimension has size 1). The blocks are
120 : !> reserved according to the neighbor lists
121 : !> \param pq_X the tensor to store the integrals
122 : !> \param ab_nl the 1st and 2nd center neighbor list
123 : !> \param ac_nl the 1st and 3rd center neighbor list
124 : !> \param matrix_dist ...
125 : !> \param blk_size_1 the block size in the first dimension
126 : !> \param blk_size_2 the block size in the second dimension
127 : !> \param blk_size_3 the block size in the third dimension
128 : !> \param only_bc_same_center only keep block if, for the corresponding integral (ab|c), b and c
129 : !> share the same center, i.e. r_bc = 0.0
130 : ! **************************************************************************************************
131 2024 : SUBROUTINE create_pqX_tensor(pq_X, ab_nl, ac_nl, matrix_dist, blk_size_1, blk_size_2, &
132 184 : blk_size_3, only_bc_same_center)
133 :
134 : TYPE(dbt_type), INTENT(OUT) :: pq_X
135 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
136 : POINTER :: ab_nl, ac_nl
137 : TYPE(dbcsr_distribution_type), INTENT(IN) :: matrix_dist
138 : INTEGER, DIMENSION(:), INTENT(IN) :: blk_size_1, blk_size_2, blk_size_3
139 : LOGICAL, INTENT(IN), OPTIONAL :: only_bc_same_center
140 :
141 : CHARACTER(len=*), PARAMETER :: routineN = 'create_pqX_tensor'
142 :
143 : INTEGER :: A, b, group_handle, handle, i, iatom, &
144 : ikind, jatom, katom, kkind, nblk, &
145 : nblk_3, nblk_per_thread, nkind
146 184 : INTEGER, ALLOCATABLE, DIMENSION(:) :: idx1, idx2, idx3
147 : INTEGER, DIMENSION(3) :: pdims
148 184 : INTEGER, DIMENSION(:), POINTER :: col_dist, row_dist
149 184 : INTEGER, DIMENSION(:, :), POINTER :: mat_pgrid
150 : LOGICAL :: my_sort_bc, symmetric
151 : REAL(dp), DIMENSION(3) :: rab, rac, rbc
152 1656 : TYPE(dbt_distribution_type) :: t_dist
153 552 : TYPE(dbt_pgrid_type) :: t_pgrid
154 : TYPE(mp_comm_type) :: group
155 : TYPE(neighbor_list_iterator_p_type), &
156 184 : DIMENSION(:), POINTER :: ab_iter, ac_iter
157 :
158 184 : NULLIFY (ab_iter, ac_iter, col_dist, row_dist, mat_pgrid)
159 :
160 184 : CALL timeset(routineN, handle)
161 :
162 184 : my_sort_bc = .FALSE.
163 184 : IF (PRESENT(only_bc_same_center)) my_sort_bc = only_bc_same_center
164 :
165 184 : CALL get_neighbor_list_set_p(ab_nl, symmetric=symmetric)
166 184 : CPASSERT(symmetric)
167 :
168 : !get 2D distribution info from matrix_dist
169 : CALL dbcsr_distribution_get(matrix_dist, pgrid=mat_pgrid, group=group_handle, &
170 184 : row_dist=row_dist, col_dist=col_dist)
171 184 : CALL group%set_handle(group_handle)
172 :
173 : !create the corresponding tensor proc grid and dist
174 184 : pdims(1) = SIZE(mat_pgrid, 1); pdims(2) = SIZE(mat_pgrid, 2); pdims(3) = 1
175 184 : CALL dbt_pgrid_create(group, pdims, t_pgrid)
176 :
177 184 : nblk_3 = SIZE(blk_size_3)
178 : CALL dbt_distribution_new(t_dist, t_pgrid, nd_dist_1=row_dist, nd_dist_2=col_dist, &
179 2964 : nd_dist_3=[(0, i=1, nblk_3)])
180 :
181 : !create the tensor itself.
182 : CALL dbt_create(pq_X, name="(pq|X)", dist=t_dist, map1_2d=[1, 2], map2_2d=[3], &
183 184 : blk_size_1=blk_size_1, blk_size_2=blk_size_2, blk_size_3=blk_size_3)
184 :
185 : !count the blocks to reserve !note: dbcsr takes care of only keeping unique indices
186 184 : CALL neighbor_list_iterator_create(ab_iter, ab_nl)
187 184 : CALL neighbor_list_iterator_create(ac_iter, ac_nl, search=.TRUE.)
188 184 : nblk = 0
189 17538 : DO WHILE (neighbor_list_iterate(ab_iter) == 0)
190 17354 : CALL get_iterator_info(ab_iter, ikind=ikind, iatom=iatom, nkind=nkind, r=rab)
191 :
192 51004 : DO kkind = 1, nkind
193 33466 : CALL nl_set_sub_iterator(ac_iter, ikind, kkind, iatom)
194 :
195 72362 : DO WHILE (nl_sub_iterate(ac_iter) == 0)
196 :
197 21542 : IF (my_sort_bc) THEN
198 : !we check for rbc or rac because of symmetry in ab_nl
199 724 : CALL get_iterator_info(ac_iter, r=rac)
200 2896 : rbc(:) = rac(:) - rab(:)
201 3196 : IF (.NOT. (ALL(ABS(rbc) <= 1.0E-8_dp) .OR. ALL(ABS(rac) <= 1.0E-8_dp))) CYCLE
202 :
203 : END IF
204 :
205 21542 : nblk = nblk + 1
206 : END DO !ac_iter
207 : END DO !kkind
208 : END DO !ab_iter
209 184 : CALL neighbor_list_iterator_release(ab_iter)
210 184 : CALL neighbor_list_iterator_release(ac_iter)
211 :
212 833 : ALLOCATE (idx1(nblk), idx2(nblk), idx3(nblk))
213 :
214 : !actually reserve the blocks
215 184 : CALL neighbor_list_iterator_create(ab_iter, ab_nl)
216 184 : CALL neighbor_list_iterator_create(ac_iter, ac_nl, search=.TRUE.)
217 184 : nblk = 0
218 17538 : DO WHILE (neighbor_list_iterate(ab_iter) == 0)
219 17354 : CALL get_iterator_info(ab_iter, ikind=ikind, iatom=iatom, jatom=jatom, nkind=nkind, r=rab)
220 :
221 51004 : DO kkind = 1, nkind
222 33466 : CALL nl_set_sub_iterator(ac_iter, ikind, kkind, iatom)
223 :
224 72362 : DO WHILE (nl_sub_iterate(ac_iter) == 0)
225 21542 : CALL get_iterator_info(ac_iter, jatom=katom, r=rac)
226 :
227 21542 : IF (my_sort_bc) THEN
228 : !we check for rbc or rac because of symmetry in ab_nl
229 724 : CALL get_iterator_info(ac_iter, r=rac)
230 2896 : rbc(:) = rac(:) - rab(:)
231 3196 : IF (.NOT. (ALL(ABS(rbc) <= 1.0E-8_dp) .OR. ALL(ABS(rac) <= 1.0E-8_dp))) CYCLE
232 :
233 : END IF
234 :
235 21248 : nblk = nblk + 1
236 :
237 21248 : idx1(nblk) = iatom
238 21248 : idx2(nblk) = jatom
239 21542 : idx3(nblk) = katom
240 :
241 : END DO !ac_iter
242 : END DO !kkind
243 : END DO !ab_iter
244 184 : CALL neighbor_list_iterator_release(ab_iter)
245 184 : CALL neighbor_list_iterator_release(ac_iter)
246 :
247 : !TODO: Parallelize creation of block list.
248 184 : !$OMP PARALLEL DEFAULT(NONE) SHARED(pq_X, nblk, idx1, idx2, idx3) PRIVATE(nblk_per_thread,A,b)
249 : nblk_per_thread = nblk/omp_get_num_threads() + 1
250 : a = omp_get_thread_num()*nblk_per_thread + 1
251 : b = MIN(a + nblk_per_thread, nblk)
252 : CALL dbt_reserve_blocks(pq_X, idx1(a:b), idx2(a:b), idx3(a:b))
253 : !$OMP END PARALLEL
254 184 : CALL dbt_finalize(pq_X)
255 :
256 : !clean-up
257 184 : CALL dbt_distribution_destroy(t_dist)
258 184 : CALL dbt_pgrid_destroy(t_pgrid)
259 :
260 184 : CALL timestop(handle)
261 :
262 920 : END SUBROUTINE create_pqX_tensor
263 :
264 : ! **************************************************************************************************
265 : !> \brief Fills the given 3 dimensional (pq|X) tensor with integrals
266 : !> \param pq_X the tensor to fill
267 : !> \param ab_nl the neighbor list for the first 2 centers
268 : !> \param ac_nl the neighbor list for the first and third centers
269 : !> \param basis_set_list_a basis sets for first center
270 : !> \param basis_set_list_b basis sets for second center
271 : !> \param basis_set_list_c basis sets for third center
272 : !> \param potential_parameter the operator for the integrals
273 : !> \param qs_env ...
274 : !> \param only_bc_same_center same as in create_pqX_tensor
275 : !> \param eps_screen threshold for possible screening
276 : !> \note The following indices are happily mixed within this routine: First center i,a,p
277 : !> Second center: j,b,q Third center: k,c,X
278 : ! **************************************************************************************************
279 184 : SUBROUTINE fill_pqX_tensor(pq_X, ab_nl, ac_nl, basis_set_list_a, basis_set_list_b, &
280 : basis_set_list_c, potential_parameter, qs_env, &
281 : only_bc_same_center, eps_screen)
282 :
283 : TYPE(dbt_type) :: pq_X
284 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
285 : POINTER :: ab_nl, ac_nl
286 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b, &
287 : basis_set_list_c
288 : TYPE(libint_potential_type) :: potential_parameter
289 : TYPE(qs_environment_type), POINTER :: qs_env
290 : LOGICAL, INTENT(IN), OPTIONAL :: only_bc_same_center
291 : REAL(dp), INTENT(IN), OPTIONAL :: eps_screen
292 :
293 : CHARACTER(len=*), PARAMETER :: routineN = 'fill_pqX_tensor'
294 :
295 : INTEGER :: egfa, egfb, egfc, handle, i, iatom, ibasis, ikind, ilist, imax, iset, jatom, &
296 : jkind, jset, katom, kkind, kset, m_max, max_nset, maxli, maxlj, maxlk, mepos, nbasis, &
297 : ncoa, ncob, ncoc, ni, nj, nk, nseta, nsetb, nsetc, nthread, sgfa, sgfb, sgfc, unit_id
298 184 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, lc_max, &
299 184 : lc_min, npgfa, npgfb, npgfc, nsgfa, &
300 184 : nsgfb, nsgfc
301 184 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb, first_sgfc
302 : LOGICAL :: do_screen, my_sort_bc
303 : REAL(dp) :: dij, dik, djk, my_eps_screen, ri(3), &
304 : rij(3), rik(3), rj(3), rjk(3), rk(3), &
305 : sabc_ext, screen_radius
306 184 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: ccp_buffer, cpp_buffer
307 184 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: max_contr, max_contra, max_contrb, &
308 184 : max_contrc
309 184 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: iabc, sabc, work
310 184 : REAL(dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b, set_radius_c
311 184 : REAL(dp), DIMENSION(:, :), POINTER :: rpgf_a, rpgf_b, rpgf_c, zeta, zetb, zetc
312 184 : TYPE(cp_2d_r_p_type), DIMENSION(:, :), POINTER :: spb, spc, tspa
313 : TYPE(cp_libint_t) :: lib
314 184 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
315 : TYPE(gto_basis_set_type), POINTER :: basis_set, basis_set_a, basis_set_b, &
316 : basis_set_c
317 : TYPE(mp_para_env_type), POINTER :: para_env
318 : TYPE(o3c_container_type), POINTER :: o3c
319 : TYPE(o3c_iterator_type) :: o3c_iterator
320 :
321 184 : NULLIFY (basis_set, basis_set_list, para_env, la_max, la_min)
322 184 : NULLIFY (lb_max, lb_min, lc_max, lc_min, npgfa, npgfb, npgfc, nsgfa, nsgfb, nsgfc)
323 184 : NULLIFY (first_sgfa, first_sgfb, first_sgfc, set_radius_a, set_radius_b, set_radius_c)
324 184 : NULLIFY (rpgf_a, rpgf_b, rpgf_c, zeta, zetb, zetc)
325 184 : NULLIFY (basis_set_a, basis_set_b, basis_set_c, tspa, spb, spc)
326 :
327 184 : CALL timeset(routineN, handle)
328 :
329 : !Need the max l for each basis for libint (and overall max #of sets for screening)
330 184 : nbasis = SIZE(basis_set_list_a)
331 184 : max_nset = 0
332 184 : maxli = 0
333 492 : DO ibasis = 1, nbasis
334 : CALL get_gto_basis_set(gto_basis_set=basis_set_list_a(ibasis)%gto_basis_set, &
335 308 : maxl=imax, nset=iset, nsgf_set=nsgfa)
336 308 : maxli = MAX(maxli, imax)
337 492 : max_nset = MAX(max_nset, iset)
338 : END DO
339 : maxlj = 0
340 492 : DO ibasis = 1, nbasis
341 : CALL get_gto_basis_set(gto_basis_set=basis_set_list_b(ibasis)%gto_basis_set, &
342 308 : maxl=imax, nset=iset, nsgf_set=nsgfb, npgf=npgfb)
343 308 : maxlj = MAX(maxlj, imax)
344 492 : max_nset = MAX(max_nset, iset)
345 : END DO
346 : maxlk = 0
347 492 : DO ibasis = 1, nbasis
348 : CALL get_gto_basis_set(gto_basis_set=basis_set_list_c(ibasis)%gto_basis_set, &
349 308 : maxl=imax, nset=iset, npgf=npgfc)
350 308 : maxlk = MAX(maxlk, imax)
351 492 : max_nset = MAX(max_nset, iset)
352 : END DO
353 184 : m_max = maxli + maxlj + maxlk
354 :
355 : !Screening
356 184 : do_screen = .FALSE.
357 184 : IF (PRESENT(eps_screen)) THEN
358 54 : do_screen = .TRUE.
359 54 : my_eps_screen = eps_screen
360 : END IF
361 184 : screen_radius = 0.0_dp
362 184 : IF (potential_parameter%potential_type == do_potential_truncated .OR. &
363 : potential_parameter%potential_type == do_potential_short) THEN
364 :
365 12 : screen_radius = potential_parameter%cutoff_radius*cutoff_screen_factor
366 172 : ELSE IF (potential_parameter%potential_type == do_potential_coulomb) THEN
367 :
368 102 : screen_radius = 1000000.0_dp
369 : END IF
370 :
371 : !get maximum contraction values for abc_contract screening
372 184 : IF (do_screen) THEN
373 :
374 : !Allocate max_contraction arrays such that we have a specific value for each set/kind
375 0 : ALLOCATE (max_contr(max_nset, nbasis), max_contra(max_nset, nbasis), &
376 540 : max_contrb(max_nset, nbasis), max_contrc(max_nset, nbasis))
377 :
378 : !Not the most elegent, but better than copying 3 times the same
379 216 : DO ilist = 1, 3
380 :
381 162 : IF (ilist == 1) basis_set_list => basis_set_list_a
382 162 : IF (ilist == 2) basis_set_list => basis_set_list_b
383 162 : IF (ilist == 3) basis_set_list => basis_set_list_c
384 :
385 162 : max_contr = 0.0_dp
386 :
387 438 : DO ibasis = 1, nbasis
388 276 : basis_set => basis_set_list(ibasis)%gto_basis_set
389 :
390 1326 : DO iset = 1, basis_set%nset
391 :
392 888 : ncoa = basis_set%npgf(iset)*ncoset(basis_set%lmax(iset))
393 888 : sgfa = basis_set%first_sgf(1, iset)
394 888 : egfa = sgfa + basis_set%nsgf_set(iset) - 1
395 :
396 : max_contr(iset, ibasis) = &
397 561714 : MAXVAL([(SUM(ABS(basis_set%sphi(1:ncoa, i))), i=sgfa, egfa)])
398 :
399 : END DO !iset
400 : END DO !ibasis
401 :
402 700 : IF (ilist == 1) max_contra(:, :) = max_contr(:, :)
403 700 : IF (ilist == 2) max_contrb(:, :) = max_contr(:, :)
404 754 : IF (ilist == 3) max_contrc(:, :) = max_contr(:, :)
405 : END DO !ilist
406 54 : DEALLOCATE (max_contr)
407 : END IF !do_screen
408 :
409 : !To minimize memory ops in contraction, we need to pre-allocate buffers, pre-tranpose sphi_a
410 : !and also trim sphi in general to have contiguous arrays
411 6872 : ALLOCATE (tspa(max_nset, nbasis), spb(max_nset, nbasis), spc(max_nset, nbasis))
412 492 : DO ibasis = 1, nbasis
413 1984 : DO iset = 1, max_nset
414 1492 : NULLIFY (tspa(iset, ibasis)%array)
415 1492 : NULLIFY (spb(iset, ibasis)%array)
416 1800 : NULLIFY (spc(iset, ibasis)%array)
417 : END DO
418 : END DO
419 :
420 736 : DO ilist = 1, 3
421 :
422 1660 : DO ibasis = 1, nbasis
423 924 : IF (ilist == 1) basis_set => basis_set_list_a(ibasis)%gto_basis_set
424 924 : IF (ilist == 2) basis_set => basis_set_list_b(ibasis)%gto_basis_set
425 924 : IF (ilist == 3) basis_set => basis_set_list_c(ibasis)%gto_basis_set
426 :
427 4536 : DO iset = 1, basis_set%nset
428 :
429 3060 : ncoa = basis_set%npgf(iset)*ncoset(basis_set%lmax(iset))
430 3060 : sgfa = basis_set%first_sgf(1, iset)
431 3060 : egfa = sgfa + basis_set%nsgf_set(iset) - 1
432 :
433 3984 : IF (ilist == 1) THEN
434 4584 : ALLOCATE (tspa(iset, ibasis)%array(basis_set%nsgf_set(iset), ncoa))
435 : #if defined(__LIBXS)
436 : CALL libxs_otrans(C_LOC(tspa(iset, ibasis)%array), &
437 : C_LOC(basis_set%sphi(1, sgfa)), &
438 : 8, ncoa, egfa - sgfa + 1, &
439 : SIZE(basis_set%sphi, 1), &
440 1146 : basis_set%nsgf_set(iset))
441 : #else
442 : tspa(iset, ibasis)%array(:, :) = TRANSPOSE(basis_set%sphi(1:ncoa, sgfa:egfa))
443 : #endif
444 1914 : ELSE IF (ilist == 2) THEN
445 4584 : ALLOCATE (spb(iset, ibasis)%array(ncoa, basis_set%nsgf_set(iset)))
446 : #if defined(__LIBXS)
447 : CALL libxs_matcopy(C_LOC(spb(iset, ibasis)%array), &
448 : C_LOC(basis_set%sphi(1, sgfa)), &
449 : 8, ncoa, egfa - sgfa + 1, &
450 1146 : SIZE(basis_set%sphi, 1), ncoa)
451 : #else
452 : spb(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoa, sgfa:egfa)
453 : #endif
454 : ELSE
455 3072 : ALLOCATE (spc(iset, ibasis)%array(ncoa, basis_set%nsgf_set(iset)))
456 : #if defined(__LIBXS)
457 : CALL libxs_matcopy(C_LOC(spc(iset, ibasis)%array), &
458 : C_LOC(basis_set%sphi(1, sgfa)), &
459 : 8, ncoa, egfa - sgfa + 1, &
460 768 : SIZE(basis_set%sphi, 1), ncoa)
461 : #else
462 : spc(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoa, sgfa:egfa)
463 : #endif
464 : END IF
465 :
466 : END DO !iset
467 : END DO !ibasis
468 : END DO !ilist
469 :
470 184 : my_sort_bc = .FALSE.
471 184 : IF (PRESENT(only_bc_same_center)) my_sort_bc = only_bc_same_center
472 :
473 : !Init the truncated Coulomb operator
474 184 : CALL get_qs_env(qs_env, para_env=para_env)
475 184 : IF (potential_parameter%potential_type == do_potential_truncated) THEN
476 :
477 : !open the file only if necessary
478 6 : IF (m_max > get_lmax_init()) THEN
479 0 : IF (para_env%mepos == 0) THEN
480 0 : CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
481 : END IF
482 0 : CALL init(m_max, unit_id, para_env%mepos, para_env)
483 0 : IF (para_env%mepos == 0) THEN
484 0 : CALL close_file(unit_id)
485 : END IF
486 : END IF
487 : END IF
488 :
489 : !Inint the initial gamma function before the OMP region as it is not thread safe
490 184 : CALL init_md_ftable(nmax=m_max)
491 :
492 : !Strategy: we use the o3c iterator because it is OMP parallelized and also takes the
493 : ! only_bc_same_center argument. Only the dbcsr_put_block is critical
494 :
495 : nthread = 1
496 184 : !$ nthread = omp_get_max_threads()
497 :
498 184 : ALLOCATE (o3c)
499 : CALL init_o3c_container(o3c, 1, basis_set_list_a, basis_set_list_b, basis_set_list_c, &
500 184 : ab_nl, ac_nl, only_bc_same_center=my_sort_bc)
501 184 : CALL o3c_iterator_create(o3c, o3c_iterator, nthread=nthread)
502 :
503 : !$OMP PARALLEL DEFAULT(NONE) &
504 : !$OMP SHARED (pq_X,do_screen,max_nset,basis_set_list_a,max_contra,max_contrb,max_contrc,&
505 : !$OMP basis_set_list_b, basis_set_list_c,ncoset,screen_radius,potential_parameter,&
506 : !$OMP my_eps_screen,maxli,maxlj,maxlk,my_sort_bc,nthread,o3c,o3c_iterator,tspa,spb,spc) &
507 : !$OMP PRIVATE (lib,i,mepos,work,iset,ncoa,sgfa,egfa,nseta,iatom,ikind,jatom,jkind,katom,kkind,&
508 : !$OMP rij,rik,rjk,basis_set_a,nsetb,la_max,la_min,lb_max,lb_min,lc_max,lc_min,npgfa,&
509 : !$OMP npgfb,npgfc,nsgfa,nsgfb,nsgfc,ri,rk,first_sgfa,first_sgfb,first_sgfc,nsetc,&
510 : !$OMP set_radius_a,set_radius_b,set_radius_c,rj,rpgf_a,rpgf_b,rpgf_c,zeta,zetb,&
511 : !$OMP zetc,basis_set_b,basis_set_c,dij,dik,djk,ni,nj,nk,iabc,sabc,jset,kset,&
512 184 : !$OMP ncob,ncoc,sgfb,sgfc,egfb,egfc,sabc_ext,cpp_buffer,ccp_buffer)
513 :
514 : mepos = 0
515 : !$ mepos = omp_get_thread_num()
516 :
517 : !each thread need its own libint object (internals may change at different rates)
518 : CALL cp_libint_init_3eri(lib, MAX(maxli, maxlj, maxlk))
519 : CALL cp_libint_set_contrdepth(lib, 1)
520 :
521 : DO WHILE (o3c_iterate(o3c_iterator, mepos=mepos) == 0)
522 : CALL get_o3c_iterator_info(o3c_iterator, mepos=mepos, ikind=ikind, jkind=jkind, kkind=kkind, &
523 : iatom=iatom, jatom=jatom, katom=katom, rij=rij, rik=rik)
524 :
525 : !get first center basis info
526 : basis_set_a => basis_set_list_a(ikind)%gto_basis_set
527 : first_sgfa => basis_set_a%first_sgf
528 : la_max => basis_set_a%lmax
529 : la_min => basis_set_a%lmin
530 : npgfa => basis_set_a%npgf
531 : nseta = basis_set_a%nset
532 : nsgfa => basis_set_a%nsgf_set
533 : zeta => basis_set_a%zet
534 : rpgf_a => basis_set_a%pgf_radius
535 : set_radius_a => basis_set_a%set_radius
536 : ni = SUM(nsgfa)
537 : !second center basis info
538 : basis_set_b => basis_set_list_b(jkind)%gto_basis_set
539 : first_sgfb => basis_set_b%first_sgf
540 : lb_max => basis_set_b%lmax
541 : lb_min => basis_set_b%lmin
542 : npgfb => basis_set_b%npgf
543 : nsetb = basis_set_b%nset
544 : nsgfb => basis_set_b%nsgf_set
545 : zetb => basis_set_b%zet
546 : rpgf_b => basis_set_b%pgf_radius
547 : set_radius_b => basis_set_b%set_radius
548 : nj = SUM(nsgfb)
549 : !third center basis info
550 : basis_set_c => basis_set_list_c(kkind)%gto_basis_set
551 : first_sgfc => basis_set_c%first_sgf
552 : lc_max => basis_set_c%lmax
553 : lc_min => basis_set_c%lmin
554 : npgfc => basis_set_c%npgf
555 : nsetc = basis_set_c%nset
556 : nsgfc => basis_set_c%nsgf_set
557 : zetc => basis_set_c%zet
558 : rpgf_c => basis_set_c%pgf_radius
559 : set_radius_c => basis_set_c%set_radius
560 : nk = SUM(nsgfc)
561 :
562 : !position and distances, only relative pos matter for libint
563 : rjk = rik - rij
564 : ri = 0.0_dp
565 : rj = rij ! ri + rij
566 : rk = rik ! ri + rik
567 :
568 : djk = NORM2(rjk)
569 : dij = NORM2(rij)
570 : dik = NORM2(rik)
571 :
572 : !sgf integrals
573 : ALLOCATE (iabc(ni, nj, nk))
574 : iabc(:, :, :) = 0.0_dp
575 :
576 : DO iset = 1, nseta
577 : ncoa = npgfa(iset)*ncoset(la_max(iset))
578 : sgfa = first_sgfa(1, iset)
579 : egfa = sgfa + nsgfa(iset) - 1
580 :
581 : DO jset = 1, nsetb
582 : ncob = npgfb(jset)*ncoset(lb_max(jset))
583 : sgfb = first_sgfb(1, jset)
584 : egfb = sgfb + nsgfb(jset) - 1
585 :
586 : !screening (overlap)
587 : IF (set_radius_a(iset) + set_radius_b(jset) < dij) CYCLE
588 :
589 : DO kset = 1, nsetc
590 : ncoc = npgfc(kset)*ncoset(lc_max(kset))
591 : sgfc = first_sgfc(1, kset)
592 : egfc = sgfc + nsgfc(kset) - 1
593 :
594 : !screening (potential)
595 : IF (set_radius_a(iset) + set_radius_c(kset) + screen_radius < dik) CYCLE
596 : IF (set_radius_b(jset) + set_radius_c(kset) + screen_radius < djk) CYCLE
597 :
598 : !pgf integrals
599 : ALLOCATE (sabc(ncoa, ncob, ncoc))
600 : sabc(:, :, :) = 0.0_dp
601 :
602 : IF (do_screen) THEN
603 : CALL eri_3center(sabc, la_min(iset), la_max(iset), npgfa(iset), zeta(:, iset), &
604 : rpgf_a(:, iset), ri, lb_min(jset), lb_max(jset), npgfb(jset), &
605 : zetb(:, jset), rpgf_b(:, jset), rj, lc_min(kset), lc_max(kset), &
606 : npgfc(kset), zetc(:, kset), rpgf_c(:, kset), rk, dij, dik, &
607 : djk, lib, potential_parameter, int_abc_ext=sabc_ext)
608 : IF (my_eps_screen > sabc_ext*(max_contra(iset, ikind)* &
609 : max_contrb(jset, jkind)* &
610 : max_contrc(kset, kkind))) THEN
611 : DEALLOCATE (sabc)
612 : CYCLE
613 : END IF
614 : ELSE
615 : CALL eri_3center(sabc, la_min(iset), la_max(iset), npgfa(iset), zeta(:, iset), &
616 : rpgf_a(:, iset), ri, lb_min(jset), lb_max(jset), npgfb(jset), &
617 : zetb(:, jset), rpgf_b(:, jset), rj, lc_min(kset), lc_max(kset), &
618 : npgfc(kset), zetc(:, kset), rpgf_c(:, kset), rk, dij, dik, &
619 : djk, lib, potential_parameter)
620 : END IF
621 :
622 : ALLOCATE (work(nsgfa(iset), nsgfb(jset), nsgfc(kset)))
623 :
624 : CALL abc_contract_xsmm(work, sabc, tspa(iset, ikind)%array, spb(jset, jkind)%array, &
625 : spc(kset, kkind)%array, ncoa, ncob, ncoc, nsgfa(iset), &
626 : nsgfb(jset), nsgfc(kset), cpp_buffer, ccp_buffer)
627 :
628 : iabc(sgfa:egfa, sgfb:egfb, sgfc:egfc) = work(:, :, :)
629 : DEALLOCATE (sabc, work)
630 :
631 : END DO !kset
632 : END DO !jset
633 : END DO !iset
634 :
635 : !Add the integral to the proper tensor block
636 : !$OMP CRITICAL(omp_fill_pqX_tensor)
637 : CALL dbt_put_block(pq_X, [iatom, jatom, katom], SHAPE(iabc), iabc, summation=.TRUE.)
638 : !$OMP END CRITICAL(omp_fill_pqX_tensor)
639 :
640 : DEALLOCATE (iabc)
641 : END DO !o3c_iterator
642 :
643 : IF (ALLOCATED(ccp_buffer)) DEALLOCATE (ccp_buffer)
644 : IF (ALLOCATED(cpp_buffer)) DEALLOCATE (cpp_buffer)
645 :
646 : CALL cp_libint_cleanup_3eri(lib)
647 :
648 : !$OMP END PARALLEL
649 184 : CALL o3c_iterator_release(o3c_iterator)
650 184 : CALL release_o3c_container(o3c)
651 184 : DEALLOCATE (o3c)
652 :
653 1122 : DO iset = 1, max_nset
654 2614 : DO ibasis = 1, nbasis
655 1492 : IF (ASSOCIATED(tspa(iset, ibasis)%array)) DEALLOCATE (tspa(iset, ibasis)%array)
656 1492 : IF (ASSOCIATED(spb(iset, ibasis)%array)) DEALLOCATE (spb(iset, ibasis)%array)
657 2430 : IF (ASSOCIATED(spc(iset, ibasis)%array)) DEALLOCATE (spc(iset, ibasis)%array)
658 : END DO
659 : END DO
660 184 : DEALLOCATE (tspa, spb, spc)
661 :
662 184 : CALL timestop(handle)
663 :
664 368 : END SUBROUTINE fill_pqX_tensor
665 :
666 : ! **************************************************************************************************
667 : !> \brief Builds a neighbor lists set for overlaping 2-center S_ab, where b is restricted on a
668 : !> a given list of atoms. Used for Coulomb RI where (aI|P) = sum_b C_bI (ab|P), where
669 : !> contraction coeff C_bI is assumed to be non-zero only on excited atoms
670 : !> \param ab_list the neighbor list
671 : !> \param basis_a basis set list for atom a
672 : !> \param basis_b basis set list for atom b
673 : !> \param qs_env ...
674 : !> \param excited_atoms the indices of the excited atoms on which b is centered
675 : !> \param ext_dist2d use an external distribution2d
676 : ! **************************************************************************************************
677 356 : SUBROUTINE build_xas_tdp_ovlp_nl(ab_list, basis_a, basis_b, qs_env, excited_atoms, ext_dist2d)
678 :
679 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
680 : POINTER :: ab_list
681 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_a, basis_b
682 : TYPE(qs_environment_type), POINTER :: qs_env
683 : INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: excited_atoms
684 : TYPE(distribution_2d_type), OPTIONAL, POINTER :: ext_dist2d
685 :
686 : INTEGER :: ikind, nkind
687 : LOGICAL :: my_restrictb
688 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: a_present, b_present
689 : REAL(dp) :: subcells
690 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: a_radius, b_radius
691 356 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius
692 356 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
693 : TYPE(cell_type), POINTER :: cell
694 : TYPE(distribution_1d_type), POINTER :: distribution_1d
695 : TYPE(distribution_2d_type), POINTER :: distribution_2d
696 356 : TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
697 356 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
698 356 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
699 :
700 356 : NULLIFY (atomic_kind_set, distribution_1d, distribution_2d, molecule_set, particle_set, cell)
701 :
702 : ! Initialization
703 356 : CALL get_qs_env(qs_env, nkind=nkind)
704 356 : CALL section_vals_val_get(qs_env%input, "DFT%SUBCELLS", r_val=subcells)
705 :
706 356 : my_restrictb = .FALSE.
707 356 : IF (PRESENT(excited_atoms)) THEN
708 124 : my_restrictb = .TRUE.
709 : END IF
710 :
711 1424 : ALLOCATE (a_present(nkind), b_present(nkind))
712 356 : a_present = .FALSE.
713 356 : b_present = .FALSE.
714 1424 : ALLOCATE (a_radius(nkind), b_radius(nkind))
715 356 : a_radius = 0.0_dp
716 356 : b_radius = 0.0_dp
717 :
718 : ! Set up the radii
719 946 : DO ikind = 1, nkind
720 590 : IF (ASSOCIATED(basis_a(ikind)%gto_basis_set)) THEN
721 590 : a_present(ikind) = .TRUE.
722 590 : CALL get_gto_basis_set(basis_a(ikind)%gto_basis_set, kind_radius=a_radius(ikind))
723 : END IF
724 :
725 946 : IF (ASSOCIATED(basis_b(ikind)%gto_basis_set)) THEN
726 590 : b_present(ikind) = .TRUE.
727 590 : CALL get_gto_basis_set(basis_b(ikind)%gto_basis_set, kind_radius=b_radius(ikind))
728 : END IF
729 : END DO !ikind
730 :
731 1424 : ALLOCATE (pair_radius(nkind, nkind))
732 356 : pair_radius = 0.0_dp
733 356 : CALL pair_radius_setup(a_present, b_present, a_radius, b_radius, pair_radius)
734 :
735 : ! Set up the nl
736 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, cell=cell, &
737 : distribution_2d=distribution_2d, local_particles=distribution_1d, &
738 356 : particle_set=particle_set, molecule_set=molecule_set)
739 :
740 : !use an external distribution_2d if required
741 356 : IF (PRESENT(ext_dist2d)) distribution_2d => ext_dist2d
742 :
743 1658 : ALLOCATE (atom2d(nkind))
744 : CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
745 356 : molecule_set, .FALSE., particle_set)
746 :
747 356 : IF (my_restrictb) THEN
748 :
749 : CALL build_neighbor_lists(ab_list, particle_set, atom2d, cell, pair_radius, subcells, &
750 124 : atomb_to_keep=excited_atoms, nlname="XAS_TDP_ovlp_nl")
751 :
752 : ELSE
753 :
754 : CALL build_neighbor_lists(ab_list, particle_set, atom2d, cell, pair_radius, subcells, &
755 232 : nlname="XAS_TDP_ovlp_nl")
756 :
757 : END IF
758 : ! Clean-up
759 356 : CALL atom2d_cleanup(atom2d)
760 :
761 712 : END SUBROUTINE build_xas_tdp_ovlp_nl
762 :
763 : ! **************************************************************************************************
764 : !> \brief Builds a neighbor lists set taylored for 3-center integral within XAS TDP, such that only
765 : !> excited atoms are taken into account for the list_c
766 : !> \param ac_list the neighbor list ready for 3-center integrals
767 : !> \param basis_a basis set list for atom a
768 : !> \param basis_c basis set list for atom c
769 : !> \param op_type to indicate whther the list should be built with overlap, Coulomb or else in mind
770 : !> \param qs_env ...
771 : !> \param excited_atoms the indices of the excited atoms to consider (if not given, all atoms are taken)
772 : !> \param x_range in case some truncated/screened operator is used, gives its range
773 : !> \param ext_dist2d external distribution_2d to be used
774 : !> \note Based on setup_neighbor_list with added features
775 : ! **************************************************************************************************
776 300 : SUBROUTINE build_xas_tdp_3c_nl(ac_list, basis_a, basis_c, op_type, qs_env, excited_atoms, &
777 : x_range, ext_dist2d)
778 :
779 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
780 : POINTER :: ac_list
781 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_a, basis_c
782 : INTEGER, INTENT(IN) :: op_type
783 : TYPE(qs_environment_type), POINTER :: qs_env
784 : INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: excited_atoms
785 : REAL(dp), INTENT(IN), OPTIONAL :: x_range
786 : TYPE(distribution_2d_type), OPTIONAL, POINTER :: ext_dist2d
787 :
788 : INTEGER :: ikind, nkind
789 : LOGICAL :: sort_atoms
790 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: a_present, c_present
791 : REAL(dp) :: subcells
792 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: a_radius, c_radius
793 300 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius
794 300 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
795 : TYPE(cell_type), POINTER :: cell
796 : TYPE(distribution_1d_type), POINTER :: distribution_1d
797 : TYPE(distribution_2d_type), POINTER :: distribution_2d
798 300 : TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
799 300 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
800 300 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
801 :
802 300 : NULLIFY (atomic_kind_set, distribution_1d, distribution_2d, molecule_set, particle_set, cell)
803 :
804 : ! Initialization
805 300 : CALL get_qs_env(qs_env, nkind=nkind)
806 300 : CALL section_vals_val_get(qs_env%input, "DFT%SUBCELLS", r_val=subcells)
807 300 : sort_atoms = .FALSE.
808 300 : IF ((PRESENT(excited_atoms))) sort_atoms = .TRUE.
809 :
810 1200 : ALLOCATE (a_present(nkind), c_present(nkind))
811 300 : a_present = .FALSE.
812 300 : c_present = .FALSE.
813 1200 : ALLOCATE (a_radius(nkind), c_radius(nkind))
814 300 : a_radius = 0.0_dp
815 300 : c_radius = 0.0_dp
816 :
817 : ! Set up the radii, depending on the operator type
818 300 : IF (op_type == do_potential_id) THEN
819 :
820 : !overlap => use the kind radius for both a and c
821 514 : DO ikind = 1, nkind
822 : !orbital basis set
823 318 : IF (ASSOCIATED(basis_a(ikind)%gto_basis_set)) THEN
824 318 : a_present(ikind) = .TRUE.
825 318 : CALL get_gto_basis_set(basis_a(ikind)%gto_basis_set, kind_radius=a_radius(ikind))
826 : END IF
827 : !RI_XAS basis set
828 514 : IF (ASSOCIATED(basis_c(ikind)%gto_basis_set)) THEN
829 318 : c_present(ikind) = .TRUE.
830 318 : CALL get_gto_basis_set(basis_c(ikind)%gto_basis_set, kind_radius=c_radius(ikind))
831 : END IF
832 : END DO !ikind
833 :
834 104 : ELSE IF (op_type == do_potential_coulomb) THEN
835 :
836 : !Coulomb operator, virtually infinite range => set c_radius to arbitrarily large number
837 236 : DO ikind = 1, nkind
838 156 : IF (ASSOCIATED(basis_c(ikind)%gto_basis_set)) THEN
839 156 : c_present(ikind) = .TRUE.
840 156 : c_radius(ikind) = 1000000.0_dp
841 : END IF
842 236 : IF (ASSOCIATED(basis_a(ikind)%gto_basis_set)) a_present(ikind) = .TRUE.
843 : END DO !ikind
844 :
845 24 : ELSE IF (op_type == do_potential_truncated .OR. op_type == do_potential_short) THEN
846 :
847 : !Truncated coulomb/short range: set c_radius to x_range + the kind_radii
848 48 : DO ikind = 1, nkind
849 24 : IF (ASSOCIATED(basis_a(ikind)%gto_basis_set)) THEN
850 24 : a_present(ikind) = .TRUE.
851 24 : CALL get_gto_basis_set(basis_a(ikind)%gto_basis_set, kind_radius=a_radius(ikind))
852 : END IF
853 48 : IF (ASSOCIATED(basis_c(ikind)%gto_basis_set)) THEN
854 24 : c_present(ikind) = .TRUE.
855 24 : CALL get_gto_basis_set(basis_c(ikind)%gto_basis_set, kind_radius=c_radius(ikind))
856 24 : c_radius(ikind) = c_radius(ikind) + x_range
857 : END IF
858 : END DO !ikind
859 :
860 : ELSE
861 0 : CPABORT("Operator not known")
862 : END IF
863 :
864 1200 : ALLOCATE (pair_radius(nkind, nkind))
865 300 : pair_radius = 0.0_dp
866 300 : CALL pair_radius_setup(a_present, c_present, a_radius, c_radius, pair_radius)
867 :
868 : ! Actually setup the list
869 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, cell=cell, &
870 : distribution_2d=distribution_2d, local_particles=distribution_1d, &
871 300 : particle_set=particle_set, molecule_set=molecule_set)
872 :
873 : !use an external distribution_2d if required
874 300 : IF (PRESENT(ext_dist2d)) distribution_2d => ext_dist2d
875 :
876 1398 : ALLOCATE (atom2d(nkind))
877 : CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
878 300 : molecule_set, .FALSE., particle_set)
879 :
880 300 : IF (sort_atoms) THEN
881 : CALL build_neighbor_lists(ac_list, particle_set, atom2d, cell, pair_radius, subcells, &
882 : operator_type="ABC", atomb_to_keep=excited_atoms, &
883 300 : nlname="XAS_TDP_3c_nl")
884 : ELSE
885 : CALL build_neighbor_lists(ac_list, particle_set, atom2d, cell, pair_radius, subcells, &
886 0 : operator_type="ABC", nlname="XAS_TDP_3c_nl")
887 : END IF
888 :
889 : ! Clean-up
890 300 : CALL atom2d_cleanup(atom2d)
891 :
892 600 : END SUBROUTINE build_xas_tdp_3c_nl
893 :
894 : ! **************************************************************************************************
895 : !> \brief Returns an optimized distribution_2d for the given neighbor lists based on an evaluation
896 : !> of the cost of the corresponding 3-center integrals
897 : !> \param opt_3c_dist2d the optimized distribution_2d
898 : !> \param ab_list ...
899 : !> \param ac_list ...
900 : !> \param basis_set_a ...
901 : !> \param basis_set_b ...
902 : !> \param basis_set_c ...
903 : !> \param qs_env ...
904 : !> \param only_bc_same_center ...
905 : ! **************************************************************************************************
906 116 : SUBROUTINE get_opt_3c_dist2d(opt_3c_dist2d, ab_list, ac_list, basis_set_a, basis_set_b, &
907 : basis_set_c, qs_env, only_bc_same_center)
908 :
909 : TYPE(distribution_2d_type), POINTER :: opt_3c_dist2d
910 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
911 : POINTER :: ab_list, ac_list
912 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_a, basis_set_b, basis_set_c
913 : TYPE(qs_environment_type), POINTER :: qs_env
914 : LOGICAL, INTENT(IN), OPTIONAL :: only_bc_same_center
915 :
916 : CHARACTER(len=*), PARAMETER :: routineN = 'get_opt_3c_dist2d'
917 :
918 : INTEGER :: handle, i, iatom, ikind, ip, jatom, &
919 : jkind, kkind, mypcol, myprow, n, &
920 : natom, nkind, npcol, nprow, nsgfa, &
921 : nsgfb, nsgfc
922 116 : INTEGER, ALLOCATABLE, DIMENSION(:) :: nparticle_local_col, nparticle_local_row
923 116 : INTEGER, DIMENSION(:, :), POINTER :: col_dist, row_dist
924 : LOGICAL :: my_sort_bc
925 : REAL(dp) :: cost, rab(3), rac(3), rbc(3)
926 116 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: col_cost, col_proc_cost, row_cost, &
927 116 : row_proc_cost
928 116 : TYPE(cp_1d_i_p_type), DIMENSION(:), POINTER :: local_particle_col, local_particle_row
929 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
930 : TYPE(mp_para_env_type), POINTER :: para_env
931 : TYPE(neighbor_list_iterator_p_type), &
932 116 : DIMENSION(:), POINTER :: ab_iter, ac_iter
933 116 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
934 116 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
935 :
936 116 : NULLIFY (para_env, col_dist, row_dist, blacs_env, qs_kind_set, particle_set)
937 116 : NULLIFY (local_particle_col, local_particle_row, ab_iter, ac_iter)
938 :
939 116 : CALL timeset(routineN, handle)
940 :
941 : !Idea: create a,b and a,c nl_iterators in the original dist, then loop over them and compute the
942 : ! cost of each ab pairs and project on the col/row. Then distribute the atom col/row to
943 : ! the proc col/row in order to spread out the cost as much as possible
944 :
945 116 : my_sort_bc = .FALSE.
946 116 : IF (PRESENT(only_bc_same_center)) my_sort_bc = only_bc_same_center
947 :
948 : CALL get_qs_env(qs_env, natom=natom, para_env=para_env, blacs_env=blacs_env, &
949 116 : qs_kind_set=qs_kind_set, particle_set=particle_set, nkind=nkind)
950 :
951 116 : myprow = blacs_env%mepos(1) + 1
952 116 : mypcol = blacs_env%mepos(2) + 1
953 116 : nprow = blacs_env%num_pe(1)
954 116 : npcol = blacs_env%num_pe(2)
955 :
956 464 : ALLOCATE (col_cost(natom), row_cost(natom))
957 116 : col_cost = 0.0_dp; row_cost = 0.0_dp
958 :
959 464 : ALLOCATE (row_dist(natom, 2), col_dist(natom, 2))
960 3212 : row_dist = 1; col_dist = 1
961 :
962 116 : CALL neighbor_list_iterator_create(ab_iter, ab_list)
963 116 : CALL neighbor_list_iterator_create(ac_iter, ac_list, search=.TRUE.)
964 1256 : DO WHILE (neighbor_list_iterate(ab_iter) == 0)
965 1140 : CALL get_iterator_info(ab_iter, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rab)
966 :
967 2846 : DO kkind = 1, nkind
968 1590 : CALL nl_set_sub_iterator(ac_iter, ikind, kkind, iatom)
969 :
970 7135 : DO WHILE (nl_sub_iterate(ac_iter) == 0)
971 :
972 4405 : IF (my_sort_bc) THEN
973 : !only take a,b,c if b,c (or a,c because of symmetry) share the same center
974 724 : CALL get_iterator_info(ac_iter, r=rac)
975 2896 : rbc(:) = rac(:) - rab(:)
976 3196 : IF (.NOT. (ALL(ABS(rbc) <= 1.0E-8_dp) .OR. ALL(ABS(rac) <= 1.0E-8_dp))) CYCLE
977 :
978 : END IF
979 :
980 : !Use the size of integral as measure as contraciton cost seems to dominate
981 4111 : nsgfa = basis_set_a(ikind)%gto_basis_set%nsgf
982 4111 : nsgfb = basis_set_b(jkind)%gto_basis_set%nsgf
983 4111 : nsgfc = basis_set_c(kkind)%gto_basis_set%nsgf
984 :
985 4111 : cost = REAL(nsgfa*nsgfb*nsgfc, dp)
986 :
987 4111 : row_cost(iatom) = row_cost(iatom) + cost
988 4405 : col_cost(jatom) = col_cost(jatom) + cost
989 :
990 : END DO !ac_iter
991 : END DO !kkind
992 : END DO !ab_iter
993 116 : CALL neighbor_list_iterator_release(ab_iter)
994 116 : CALL neighbor_list_iterator_release(ac_iter)
995 :
996 116 : CALL para_env%sum(row_cost)
997 116 : CALL para_env%sum(col_cost)
998 :
999 : !Distribute the cost as evenly as possible
1000 580 : ALLOCATE (col_proc_cost(npcol), row_proc_cost(nprow))
1001 116 : col_proc_cost = 0.0_dp; row_proc_cost = 0.0_dp
1002 774 : DO i = 1, natom
1003 21708 : iatom = MAXLOC(row_cost, 1)
1004 2632 : ip = MINLOC(row_proc_cost, 1)
1005 658 : row_proc_cost(ip) = row_proc_cost(ip) + row_cost(iatom)
1006 658 : row_dist(iatom, 1) = ip
1007 658 : row_cost(iatom) = 0.0_dp
1008 :
1009 21708 : iatom = MAXLOC(col_cost, 1)
1010 1974 : ip = MINLOC(col_proc_cost, 1)
1011 658 : col_proc_cost(ip) = col_proc_cost(ip) + col_cost(iatom)
1012 658 : col_dist(iatom, 1) = ip
1013 774 : col_cost(iatom) = 0.0_dp
1014 : END DO
1015 :
1016 : !the usual stuff
1017 844 : ALLOCATE (local_particle_col(nkind), local_particle_row(nkind))
1018 464 : ALLOCATE (nparticle_local_row(nkind), nparticle_local_col(nkind))
1019 116 : nparticle_local_row = 0; nparticle_local_col = 0
1020 :
1021 774 : DO iatom = 1, natom
1022 658 : ikind = particle_set(iatom)%atomic_kind%kind_number
1023 :
1024 658 : IF (row_dist(iatom, 1) == myprow) nparticle_local_row(ikind) = nparticle_local_row(ikind) + 1
1025 774 : IF (col_dist(iatom, 1) == mypcol) nparticle_local_col(ikind) = nparticle_local_col(ikind) + 1
1026 : END DO
1027 :
1028 306 : DO ikind = 1, nkind
1029 190 : n = nparticle_local_row(ikind)
1030 505 : ALLOCATE (local_particle_row(ikind)%array(n))
1031 :
1032 190 : n = nparticle_local_col(ikind)
1033 686 : ALLOCATE (local_particle_col(ikind)%array(n))
1034 : END DO
1035 :
1036 116 : nparticle_local_row = 0; nparticle_local_col = 0
1037 774 : DO iatom = 1, natom
1038 658 : ikind = particle_set(iatom)%atomic_kind%kind_number
1039 :
1040 658 : IF (row_dist(iatom, 1) == myprow) THEN
1041 329 : nparticle_local_row(ikind) = nparticle_local_row(ikind) + 1
1042 329 : local_particle_row(ikind)%array(nparticle_local_row(ikind)) = iatom
1043 : END IF
1044 774 : IF (col_dist(iatom, 1) == mypcol) THEN
1045 658 : nparticle_local_col(ikind) = nparticle_local_col(ikind) + 1
1046 658 : local_particle_col(ikind)%array(nparticle_local_col(ikind)) = iatom
1047 : END IF
1048 : END DO
1049 :
1050 : !Finally create the dist_2d
1051 774 : row_dist(:, 1) = row_dist(:, 1) - 1
1052 774 : col_dist(:, 1) = col_dist(:, 1) - 1
1053 : CALL distribution_2d_create(opt_3c_dist2d, row_distribution_ptr=row_dist, &
1054 : col_distribution_ptr=col_dist, local_rows_ptr=local_particle_row, &
1055 116 : local_cols_ptr=local_particle_col, blacs_env=blacs_env)
1056 :
1057 116 : CALL timestop(handle)
1058 :
1059 232 : END SUBROUTINE get_opt_3c_dist2d
1060 :
1061 : ! **************************************************************************************************
1062 : !> \brief Computes the RI exchange 3-center integrals (ab|c), where c is from the RI_XAS basis and
1063 : !> centered on excited atoms and kind. The operator used is that of the RI metric
1064 : !> \param ex_atoms excited atoms on which the third center is located
1065 : !> \param xas_tdp_env ...
1066 : !> \param xas_tdp_control ...
1067 : !> \param qs_env ...
1068 : !> \note This routine is called once for each excited atom. Because there are many different a,b
1069 : !> pairs involved, load balance is ok. This allows memory saving
1070 : ! **************************************************************************************************
1071 54 : SUBROUTINE compute_ri_3c_exchange(ex_atoms, xas_tdp_env, xas_tdp_control, qs_env)
1072 :
1073 : INTEGER, DIMENSION(:), INTENT(IN) :: ex_atoms
1074 : TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
1075 : TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
1076 : TYPE(qs_environment_type), POINTER :: qs_env
1077 :
1078 : CHARACTER(len=*), PARAMETER :: routineN = 'compute_ri_3c_exchange'
1079 :
1080 : INTEGER :: handle, natom, nkind
1081 : INTEGER, ALLOCATABLE, DIMENSION(:) :: blk_size_orb, blk_size_ri
1082 : TYPE(dbcsr_distribution_type) :: opt_dbcsr_dist
1083 54 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_orb, basis_set_ri
1084 : TYPE(mp_para_env_type), POINTER :: para_env
1085 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1086 54 : POINTER :: ab_list, ac_list
1087 54 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1088 54 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1089 :
1090 54 : NULLIFY (basis_set_ri, basis_set_orb, ac_list, ab_list, qs_kind_set, para_env, particle_set)
1091 :
1092 54 : CALL timeset(routineN, handle)
1093 :
1094 : ! Take what we need from the qs_env
1095 : CALL get_qs_env(qs_env, nkind=nkind, qs_kind_set=qs_kind_set, para_env=para_env, &
1096 54 : natom=natom, particle_set=particle_set)
1097 :
1098 : ! Build the basis set lists
1099 254 : ALLOCATE (basis_set_ri(nkind))
1100 200 : ALLOCATE (basis_set_orb(nkind))
1101 54 : CALL basis_set_list_setup(basis_set_ri, "RI_XAS", qs_kind_set)
1102 54 : CALL basis_set_list_setup(basis_set_orb, "ORB", qs_kind_set)
1103 :
1104 : ! Get the optimized distribution 2d for theses integrals (and store it in xas_tdp_env)
1105 54 : CALL build_xas_tdp_ovlp_nl(ab_list, basis_set_orb, basis_set_orb, qs_env)
1106 : CALL build_xas_tdp_3c_nl(ac_list, basis_set_orb, basis_set_ri, &
1107 : xas_tdp_control%ri_m_potential%potential_type, qs_env, &
1108 54 : excited_atoms=ex_atoms, x_range=xas_tdp_control%ri_m_potential%cutoff_radius)
1109 :
1110 : CALL get_opt_3c_dist2d(xas_tdp_env%opt_dist2d_ex, ab_list, ac_list, basis_set_orb, &
1111 54 : basis_set_orb, basis_set_ri, qs_env)
1112 54 : CALL release_neighbor_list_sets(ab_list)
1113 54 : CALL release_neighbor_list_sets(ac_list)
1114 :
1115 : ! Build the ab and ac centers neighbor lists based on the optimized distribution
1116 : CALL build_xas_tdp_ovlp_nl(ab_list, basis_set_orb, basis_set_orb, qs_env, &
1117 54 : ext_dist2d=xas_tdp_env%opt_dist2d_ex)
1118 : CALL build_xas_tdp_3c_nl(ac_list, basis_set_orb, basis_set_ri, &
1119 : xas_tdp_control%ri_m_potential%potential_type, qs_env, &
1120 : excited_atoms=ex_atoms, x_range=xas_tdp_control%ri_m_potential%cutoff_radius, &
1121 54 : ext_dist2d=xas_tdp_env%opt_dist2d_ex)
1122 :
1123 : ! Allocate, init and compute the integrals.
1124 216 : ALLOCATE (blk_size_orb(natom), blk_size_ri(natom))
1125 54 : CALL cp_dbcsr_dist2d_to_dist(xas_tdp_env%opt_dist2d_ex, opt_dbcsr_dist)
1126 54 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=blk_size_orb, basis=basis_set_orb)
1127 54 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=blk_size_ri, basis=basis_set_ri)
1128 :
1129 540 : ALLOCATE (xas_tdp_env%ri_3c_ex)
1130 : CALL create_pqX_tensor(xas_tdp_env%ri_3c_ex, ab_list, ac_list, opt_dbcsr_dist, blk_size_orb, &
1131 54 : blk_size_orb, blk_size_ri)
1132 : CALL fill_pqX_tensor(xas_tdp_env%ri_3c_ex, ab_list, ac_list, basis_set_orb, basis_set_orb, &
1133 : basis_set_ri, xas_tdp_control%ri_m_potential, qs_env, &
1134 54 : eps_screen=xas_tdp_control%eps_screen)
1135 :
1136 : ! Clean-up
1137 54 : CALL release_neighbor_list_sets(ab_list)
1138 54 : CALL release_neighbor_list_sets(ac_list)
1139 54 : CALL dbcsr_distribution_release(opt_dbcsr_dist)
1140 54 : DEALLOCATE (basis_set_ri, basis_set_orb)
1141 :
1142 : !not strictly necessary but avoid having any load unbalance here being reported in the
1143 : !timings for other routines
1144 54 : CALL para_env%sync()
1145 :
1146 54 : CALL timestop(handle)
1147 :
1148 162 : END SUBROUTINE compute_ri_3c_exchange
1149 :
1150 : ! **************************************************************************************************
1151 : !> \brief Computes the RI Coulomb 3-center integrals (ab|c), where c is from the RI_XAS basis and
1152 : !> centered on the excited atoms of xas_tdp_env
1153 : !> \param xas_tdp_env ...
1154 : !> \param qs_env ...
1155 : !> \note The ri_3c_coul tensor of xas_tdp_env is defined and allocated here. Only computed once
1156 : !> for the whole system (for optimized load balance). Ok because not too much memory needed
1157 : ! **************************************************************************************************
1158 62 : SUBROUTINE compute_ri_3c_coulomb(xas_tdp_env, qs_env)
1159 :
1160 : TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
1161 : TYPE(qs_environment_type), POINTER :: qs_env
1162 :
1163 : CHARACTER(len=*), PARAMETER :: routineN = 'compute_ri_3c_coulomb'
1164 :
1165 : INTEGER :: handle, natom, nkind
1166 : INTEGER, ALLOCATABLE, DIMENSION(:) :: blk_size_orb, blk_size_ri
1167 : TYPE(dbcsr_distribution_type) :: opt_dbcsr_dist
1168 62 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_orb, basis_set_ri
1169 : TYPE(libint_potential_type) :: pot
1170 : TYPE(mp_para_env_type), POINTER :: para_env
1171 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1172 62 : POINTER :: ab_list, ac_list
1173 62 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1174 62 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1175 :
1176 62 : NULLIFY (basis_set_ri, basis_set_orb, ac_list, ab_list, qs_kind_set, para_env, particle_set)
1177 :
1178 62 : CALL timeset(routineN, handle)
1179 :
1180 : ! Take what we need from the qs_env
1181 : CALL get_qs_env(qs_env, nkind=nkind, qs_kind_set=qs_kind_set, para_env=para_env, &
1182 62 : natom=natom, particle_set=particle_set)
1183 :
1184 : ! Build the basis set lists
1185 284 : ALLOCATE (basis_set_ri(nkind))
1186 222 : ALLOCATE (basis_set_orb(nkind))
1187 62 : CALL basis_set_list_setup(basis_set_ri, "RI_XAS", qs_kind_set)
1188 62 : CALL basis_set_list_setup(basis_set_orb, "ORB", qs_kind_set)
1189 :
1190 : ! Get the optimized distribution 2d for these integrals (and store it in xas_tdp_env)
1191 : CALL build_xas_tdp_ovlp_nl(ab_list, basis_set_orb, basis_set_orb, qs_env, &
1192 62 : excited_atoms=xas_tdp_env%ex_atom_indices)
1193 : CALL build_xas_tdp_3c_nl(ac_list, basis_set_orb, basis_set_ri, do_potential_id, &
1194 62 : qs_env, excited_atoms=xas_tdp_env%ex_atom_indices)
1195 : CALL get_opt_3c_dist2d(xas_tdp_env%opt_dist2d_coul, ab_list, ac_list, basis_set_orb, &
1196 62 : basis_set_orb, basis_set_ri, qs_env, only_bc_same_center=.TRUE.)
1197 62 : CALL release_neighbor_list_sets(ab_list)
1198 62 : CALL release_neighbor_list_sets(ac_list)
1199 :
1200 : ! Build a neighbor list for the ab centers. Assume (aI|c) = sum_b c_bI (ab|c), with c_bI only
1201 : ! non-zero for b centered on the same atom as c => build overlap nl, but only keeping b if centered
1202 : ! on an excited atom
1203 : CALL build_xas_tdp_ovlp_nl(ab_list, basis_set_orb, basis_set_orb, qs_env, &
1204 : excited_atoms=xas_tdp_env%ex_atom_indices, &
1205 62 : ext_dist2d=xas_tdp_env%opt_dist2d_coul)
1206 :
1207 : ! Build a neighbor list for the ac centers. Since we later contract as (aI|c) and we assume I is
1208 : ! very localized on the same atom as c, we take a,c as neighbors if they overlap
1209 : CALL build_xas_tdp_3c_nl(ac_list, basis_set_orb, basis_set_ri, do_potential_id, &
1210 : qs_env, excited_atoms=xas_tdp_env%ex_atom_indices, &
1211 62 : ext_dist2d=xas_tdp_env%opt_dist2d_coul)
1212 :
1213 : ! Allocate, init and compute the integrals
1214 248 : ALLOCATE (blk_size_orb(natom), blk_size_ri(natom))
1215 62 : CALL cp_dbcsr_dist2d_to_dist(xas_tdp_env%opt_dist2d_coul, opt_dbcsr_dist)
1216 62 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=blk_size_orb, basis=basis_set_orb)
1217 62 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=blk_size_ri, basis=basis_set_ri)
1218 : pot%potential_type = do_potential_coulomb
1219 :
1220 620 : ALLOCATE (xas_tdp_env%ri_3c_coul)
1221 : CALL create_pqX_tensor(xas_tdp_env%ri_3c_coul, ab_list, ac_list, opt_dbcsr_dist, blk_size_orb, &
1222 62 : blk_size_orb, blk_size_ri, only_bc_same_center=.TRUE.)
1223 : CALL fill_pqX_tensor(xas_tdp_env%ri_3c_coul, ab_list, ac_list, basis_set_orb, basis_set_orb, &
1224 62 : basis_set_ri, pot, qs_env, only_bc_same_center=.TRUE.)
1225 :
1226 : ! Clean-up
1227 62 : CALL release_neighbor_list_sets(ab_list)
1228 62 : CALL release_neighbor_list_sets(ac_list)
1229 62 : CALL dbcsr_distribution_release(opt_dbcsr_dist)
1230 62 : DEALLOCATE (basis_set_ri, basis_set_orb)
1231 :
1232 : !not strictly necessary but avoid having any load unbalance here being reported in the
1233 : !timings for other routines
1234 62 : CALL para_env%sync()
1235 :
1236 62 : CALL timestop(handle)
1237 :
1238 186 : END SUBROUTINE compute_ri_3c_coulomb
1239 :
1240 : ! **************************************************************************************************
1241 : !> \brief Computes the two-center Exchange integral needed for the RI in kernel calculation. Stores
1242 : !> the integrals in the xas_tdp_env as global (small) arrays. Does that for a given excited
1243 : !> kind. The quantity stored is M^-1 (P|Q) M^-1, where M is the RI metric. If the metric is
1244 : !> the same as the exchange potential, then we end up with the V-approximation (P|Q)^-1
1245 : !> By default (if no metric), the ri_m_potential is a copy of the x_potential
1246 : !> \param ex_kind ...
1247 : !> \param xas_tdp_env ...
1248 : !> \param xas_tdp_control ...
1249 : !> \param qs_env ...
1250 : !> \note Computes all these integrals in non-PBCs as we assume that the range is short enough that
1251 : !> atoms do not exchange with their periodic images
1252 : ! **************************************************************************************************
1253 54 : SUBROUTINE compute_ri_exchange2_int(ex_kind, xas_tdp_env, xas_tdp_control, qs_env)
1254 :
1255 : INTEGER, INTENT(IN) :: ex_kind
1256 : TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
1257 : TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
1258 : TYPE(qs_environment_type), POINTER :: qs_env
1259 :
1260 : INTEGER :: egfp, egfq, maxl, ncop, ncoq, nset, &
1261 : nsgf, pset, qset, sgfp, sgfq, unit_id
1262 54 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf_set, nsgf_set
1263 54 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf
1264 : REAL(dp) :: r(3)
1265 54 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: metric, pq, work
1266 54 : REAL(dp), DIMENSION(:, :), POINTER :: rpgf, sphi, zet
1267 : TYPE(cp_libint_t) :: lib
1268 : TYPE(gto_basis_set_type), POINTER :: ri_basis
1269 : TYPE(mp_para_env_type), POINTER :: para_env
1270 54 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1271 :
1272 54 : NULLIFY (ri_basis, qs_kind_set, para_env, lmin, lmax, npgf_set, zet, rpgf, first_sgf)
1273 54 : NULLIFY (sphi, nsgf_set)
1274 :
1275 : ! Initialization
1276 54 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, para_env=para_env)
1277 54 : IF (ASSOCIATED(xas_tdp_env%ri_inv_ex)) THEN
1278 8 : DEALLOCATE (xas_tdp_env%ri_inv_ex)
1279 : END IF
1280 :
1281 : ! Get the RI basis of interest and its quantum numbers
1282 54 : CALL get_qs_kind(qs_kind_set(ex_kind), basis_set=ri_basis, basis_type="RI_XAS")
1283 : CALL get_gto_basis_set(ri_basis, nsgf=nsgf, maxl=maxl, npgf=npgf_set, lmin=lmin, &
1284 : lmax=lmax, zet=zet, pgf_radius=rpgf, first_sgf=first_sgf, &
1285 54 : nsgf_set=nsgf_set, sphi=sphi, nset=nset)
1286 216 : ALLOCATE (metric(nsgf, nsgf))
1287 54 : metric = 0.0_dp
1288 :
1289 54 : r = 0.0_dp
1290 :
1291 : !init the libint 2-center object
1292 54 : CALL cp_libint_init_2eri(lib, maxl)
1293 54 : CALL cp_libint_set_contrdepth(lib, 1)
1294 :
1295 : !make sure the truncted coulomb is initialized
1296 54 : IF (xas_tdp_control%ri_m_potential%potential_type == do_potential_truncated) THEN
1297 :
1298 6 : IF (2*maxl + 1 > get_lmax_init()) THEN
1299 6 : IF (para_env%mepos == 0) THEN
1300 3 : CALL open_file(unit_number=unit_id, file_name=xas_tdp_control%ri_m_potential%filename)
1301 : END IF
1302 6 : CALL init(2*maxl + 1, unit_id, para_env%mepos, para_env)
1303 6 : IF (para_env%mepos == 0) THEN
1304 3 : CALL close_file(unit_id)
1305 : END IF
1306 : END IF
1307 : END IF
1308 :
1309 : ! Compute the RI metric
1310 162 : DO pset = 1, nset
1311 108 : ncop = npgf_set(pset)*ncoset(lmax(pset))
1312 108 : sgfp = first_sgf(1, pset)
1313 108 : egfp = sgfp + nsgf_set(pset) - 1
1314 :
1315 378 : DO qset = 1, nset
1316 216 : ncoq = npgf_set(qset)*ncoset(lmax(qset))
1317 216 : sgfq = first_sgf(1, qset)
1318 216 : egfq = sgfq + nsgf_set(qset) - 1
1319 :
1320 864 : ALLOCATE (work(ncop, ncoq))
1321 216 : work = 0.0_dp
1322 :
1323 : CALL eri_2center(work, lmin(pset), lmax(pset), npgf_set(pset), zet(:, pset), rpgf(:, pset), &
1324 : r, lmin(qset), lmax(qset), npgf_set(qset), zet(:, qset), rpgf(:, qset), &
1325 216 : r, 0.0_dp, lib, xas_tdp_control%ri_m_potential)
1326 :
1327 : CALL ab_contract(metric(sgfp:egfp, sgfq:egfq), work, sphi(:, sgfp:), sphi(:, sgfq:), &
1328 216 : ncop, ncoq, nsgf_set(pset), nsgf_set(qset))
1329 :
1330 324 : DEALLOCATE (work)
1331 : END DO !qset
1332 : END DO !pset
1333 :
1334 : ! Inverting (to M^-1)
1335 54 : CALL invmat_symm(metric)
1336 :
1337 54 : IF (.NOT. xas_tdp_control%do_ri_metric) THEN
1338 :
1339 : !If no metric, then x_pot = ri_m_pot and (P|Q)^-1 = M^-1 (V-approximation)
1340 144 : ALLOCATE (xas_tdp_env%ri_inv_ex(nsgf, nsgf))
1341 445396 : xas_tdp_env%ri_inv_ex(:, :) = metric(:, :)
1342 48 : CALL cp_libint_cleanup_2eri(lib)
1343 48 : RETURN
1344 :
1345 : END IF
1346 :
1347 : !make sure the truncted coulomb is initialized
1348 6 : IF (xas_tdp_control%x_potential%potential_type == do_potential_truncated) THEN
1349 :
1350 2 : IF (2*maxl + 1 > get_lmax_init()) THEN
1351 0 : IF (para_env%mepos == 0) THEN
1352 0 : CALL open_file(unit_number=unit_id, file_name=xas_tdp_control%x_potential%filename)
1353 : END IF
1354 0 : CALL init(2*maxl + 1, unit_id, para_env%mepos, para_env)
1355 0 : IF (para_env%mepos == 0) THEN
1356 0 : CALL close_file(unit_id)
1357 : END IF
1358 : END IF
1359 : END IF
1360 :
1361 : ! Compute the proper exchange 2-center
1362 18 : ALLOCATE (pq(nsgf, nsgf))
1363 6 : pq = 0.0_dp
1364 :
1365 18 : DO pset = 1, nset
1366 12 : ncop = npgf_set(pset)*ncoset(lmax(pset))
1367 12 : sgfp = first_sgf(1, pset)
1368 12 : egfp = sgfp + nsgf_set(pset) - 1
1369 :
1370 42 : DO qset = 1, nset
1371 24 : ncoq = npgf_set(qset)*ncoset(lmax(qset))
1372 24 : sgfq = first_sgf(1, qset)
1373 24 : egfq = sgfq + nsgf_set(qset) - 1
1374 :
1375 96 : ALLOCATE (work(ncop, ncoq))
1376 24 : work = 0.0_dp
1377 :
1378 : CALL eri_2center(work, lmin(pset), lmax(pset), npgf_set(pset), zet(:, pset), rpgf(:, pset), &
1379 : r, lmin(qset), lmax(qset), npgf_set(qset), zet(:, qset), rpgf(:, qset), &
1380 24 : r, 0.0_dp, lib, xas_tdp_control%x_potential)
1381 :
1382 : CALL ab_contract(pq(sgfp:egfp, sgfq:egfq), work, sphi(:, sgfp:), sphi(:, sgfq:), &
1383 24 : ncop, ncoq, nsgf_set(pset), nsgf_set(qset))
1384 :
1385 36 : DEALLOCATE (work)
1386 : END DO !qset
1387 : END DO !pset
1388 :
1389 : ! Compute and store M^-1 (P|Q) M^-1
1390 18 : ALLOCATE (xas_tdp_env%ri_inv_ex(nsgf, nsgf))
1391 22698 : xas_tdp_env%ri_inv_ex = 0.0_dp
1392 :
1393 : CALL dgemm('N', 'N', nsgf, nsgf, nsgf, 1.0_dp, metric, nsgf, pq, nsgf, &
1394 6 : 0.0_dp, xas_tdp_env%ri_inv_ex, nsgf)
1395 : CALL dgemm('N', 'N', nsgf, nsgf, nsgf, 1.0_dp, xas_tdp_env%ri_inv_ex, nsgf, metric, nsgf, &
1396 6 : 0.0_dp, pq, nsgf)
1397 22698 : xas_tdp_env%ri_inv_ex(:, :) = pq(:, :)
1398 :
1399 6 : CALL cp_libint_cleanup_2eri(lib)
1400 :
1401 162 : END SUBROUTINE compute_ri_exchange2_int
1402 :
1403 : ! **************************************************************************************************
1404 : !> \brief Computes the two-center Coulomb integral needed for the RI in kernel calculation. Stores
1405 : !> the integrals (P|Q)^-1 in the xas_tdp_env as global (small) arrays. Does that for a given
1406 : !> excited kind
1407 : !> \param ex_kind ...
1408 : !> \param xas_tdp_env ...
1409 : !> \param xas_tdp_control ...
1410 : !> \param qs_env ...
1411 : ! **************************************************************************************************
1412 74 : SUBROUTINE compute_ri_coulomb2_int(ex_kind, xas_tdp_env, xas_tdp_control, qs_env)
1413 :
1414 : INTEGER, INTENT(IN) :: ex_kind
1415 : TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
1416 : TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
1417 : TYPE(qs_environment_type), POINTER :: qs_env
1418 :
1419 : INTEGER :: nsgf
1420 : TYPE(gto_basis_set_type), POINTER :: ri_basis
1421 74 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1422 :
1423 74 : NULLIFY (ri_basis, qs_kind_set)
1424 :
1425 : ! Initialization
1426 74 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
1427 74 : IF (ASSOCIATED(xas_tdp_env%ri_inv_coul)) THEN
1428 8 : DEALLOCATE (xas_tdp_env%ri_inv_coul)
1429 : END IF
1430 :
1431 : ! Get the RI basis of interest and its quantum numbers
1432 74 : CALL get_qs_kind(qs_kind_set(ex_kind), basis_set=ri_basis, basis_type="RI_XAS")
1433 74 : CALL get_gto_basis_set(ri_basis, nsgf=nsgf)
1434 296 : ALLOCATE (xas_tdp_env%ri_inv_coul(nsgf, nsgf))
1435 660186 : xas_tdp_env%ri_inv_coul = 0.0_dp
1436 :
1437 74 : IF (.NOT. xas_tdp_control%is_periodic) THEN
1438 : CALL int_operators_r12_ab_os(r12_operator=operator_coulomb, vab=xas_tdp_env%ri_inv_coul, &
1439 : rab=[0.0_dp, 0.0_dp, 0.0_dp], fba=ri_basis, fbb=ri_basis, &
1440 58 : calculate_forces=.FALSE.)
1441 58 : CPASSERT(ASSOCIATED(xas_tdp_control))
1442 : ELSE
1443 16 : CALL periodic_ri_coulomb2(xas_tdp_env%ri_inv_coul, ri_basis, qs_env)
1444 : END IF
1445 :
1446 : ! Inverting
1447 74 : CALL invmat_symm(xas_tdp_env%ri_inv_coul)
1448 :
1449 148 : END SUBROUTINE compute_ri_coulomb2_int
1450 :
1451 : ! **************************************************************************************************
1452 : !> \brief Computes the two-center inverse coulomb integral in the case of PBCs
1453 : !> \param ri_coul2 ...
1454 : !> \param ri_basis ...
1455 : !> \param qs_env ...
1456 : ! **************************************************************************************************
1457 16 : SUBROUTINE periodic_ri_coulomb2(ri_coul2, ri_basis, qs_env)
1458 :
1459 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: ri_coul2
1460 : TYPE(gto_basis_set_type), POINTER :: ri_basis
1461 : TYPE(qs_environment_type), POINTER :: qs_env
1462 :
1463 : INTEGER :: maxco, ncop, ncoq, nset, op, oq, ppgf, &
1464 : pset, qpgf, qset, sgfp, sgfq
1465 32 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nsgf
1466 16 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf
1467 : REAL(dp) :: r(3)
1468 16 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: hpq
1469 16 : REAL(dp), DIMENSION(:, :), POINTER :: sphi, zet
1470 : TYPE(cell_type), POINTER :: cell
1471 416 : TYPE(cp_eri_mme_param) :: mme_param
1472 : TYPE(mp_para_env_type), POINTER :: para_env
1473 16 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1474 :
1475 16 : NULLIFY (cell, qs_kind_set, lmin, lmax, nsgf, npgf, zet, sphi, first_sgf)
1476 :
1477 : ! Use eri_mme for this. Don't want to annoy the user with a full input section just for this
1478 : ! tiny bit => initialize our own eri_mme section with the defaults
1479 :
1480 16 : CALL get_qs_env(qs_env, cell=cell, qs_kind_set=qs_kind_set, para_env=para_env)
1481 :
1482 : CALL eri_mme_init(mme_param%par, n_minimax=20, cutoff=300.0_dp, do_calib_cutoff=.TRUE., &
1483 : cutoff_min=10.0_dp, cutoff_max=10000.0_dp, cutoff_eps=0.01_dp, &
1484 : cutoff_delta=0.9_dp, sum_precision=1.0E-12_dp, debug=.FALSE., &
1485 : debug_delta=1.0E-6_dp, debug_nsum=1000000, unit_nr=0, print_calib=.FALSE., &
1486 16 : do_error_est=.FALSE.)
1487 16 : mme_param%do_calib = .TRUE.
1488 :
1489 16 : CALL cp_eri_mme_set_params(mme_param, cell, qs_kind_set, basis_type_1="RI_XAS", para_env=para_env)
1490 :
1491 : CALL get_gto_basis_set(ri_basis, lmax=lmax, npgf=npgf, zet=zet, lmin=lmin, nset=nset, &
1492 16 : nsgf_set=nsgf, sphi=sphi, first_sgf=first_sgf, maxco=maxco)
1493 :
1494 16 : r = 0.0_dp
1495 64 : ALLOCATE (hpq(nset*maxco, nset*maxco))
1496 16 : hpq = 0.0_dp
1497 :
1498 48 : DO pset = 1, nset
1499 32 : ncop = npgf(pset)*ncoset(lmax(pset))
1500 32 : sgfp = first_sgf(1, pset)
1501 :
1502 112 : DO qset = 1, nset
1503 64 : ncoq = npgf(qset)*ncoset(lmax(qset))
1504 64 : sgfq = first_sgf(1, qset)
1505 :
1506 608 : DO ppgf = 1, npgf(pset)
1507 544 : op = (pset - 1)*maxco + (ppgf - 1)*ncoset(lmax(pset))
1508 5232 : DO qpgf = 1, npgf(qset)
1509 4624 : oq = (qset - 1)*maxco + (qpgf - 1)*ncoset(lmax(qset))
1510 :
1511 : CALL eri_mme_2c_integrate(mme_param%par, lmin(pset), lmax(pset), lmin(qset), &
1512 : lmax(qset), zet(ppgf, pset), zet(qpgf, qset), r, hpq, &
1513 5168 : op, oq)
1514 :
1515 : END DO !qpgf
1516 : END DO ! ppgf
1517 :
1518 : !contraction into sgfs
1519 64 : op = (pset - 1)*maxco + 1
1520 64 : oq = (qset - 1)*maxco + 1
1521 :
1522 : CALL ab_contract(ri_coul2(sgfp:sgfp + nsgf(pset) - 1, sgfq:sgfq + nsgf(qset) - 1), &
1523 : hpq(op:op + ncop - 1, oq:oq + ncoq - 1), sphi(:, sgfp:), sphi(:, sgfq:), &
1524 96 : ncop, ncoq, nsgf(pset), nsgf(qset))
1525 :
1526 : END DO !qset
1527 : END DO !pset
1528 :
1529 : !celan-up
1530 16 : CALL eri_mme_release(mme_param%par)
1531 :
1532 64 : END SUBROUTINE periodic_ri_coulomb2
1533 :
1534 : END MODULE xas_tdp_integrals
|
