(git:1f9fd2c)
Loading...
Searching...
No Matches
xas_tdp_integrals.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief 3-center integrals machinery for the XAS_TDP method
10!> \author A. Bussy (03.2020)
11! **************************************************************************************************
12
13MODULE 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
114CONTAINS
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 SUBROUTINE create_pqx_tensor(pq_X, ab_nl, ac_nl, matrix_dist, blk_size_1, blk_size_2, &
132 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: idx1, idx2, idx3
147 INTEGER, DIMENSION(3) :: pdims
148 INTEGER, DIMENSION(:), POINTER :: col_dist, row_dist
149 INTEGER, DIMENSION(:, :), POINTER :: mat_pgrid
150 LOGICAL :: my_sort_bc, symmetric
151 REAL(dp), DIMENSION(3) :: rab, rac, rbc
152 TYPE(dbt_distribution_type) :: t_dist
153 TYPE(dbt_pgrid_type) :: t_pgrid
154 TYPE(mp_comm_type) :: group
155 TYPE(neighbor_list_iterator_p_type), &
156 DIMENSION(:), POINTER :: ab_iter, ac_iter
157
158 NULLIFY (ab_iter, ac_iter, col_dist, row_dist, mat_pgrid)
159
160 CALL timeset(routinen, handle)
161
162 my_sort_bc = .false.
163 IF (PRESENT(only_bc_same_center)) my_sort_bc = only_bc_same_center
164
165 CALL get_neighbor_list_set_p(ab_nl, symmetric=symmetric)
166 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 row_dist=row_dist, col_dist=col_dist)
171 CALL group%set_handle(group_handle)
172
173 !create the corresponding tensor proc grid and dist
174 pdims(1) = SIZE(mat_pgrid, 1); pdims(2) = SIZE(mat_pgrid, 2); pdims(3) = 1
175 CALL dbt_pgrid_create(group, pdims, t_pgrid)
176
177 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 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 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 CALL neighbor_list_iterator_create(ab_iter, ab_nl)
187 CALL neighbor_list_iterator_create(ac_iter, ac_nl, search=.true.)
188 nblk = 0
189 DO WHILE (neighbor_list_iterate(ab_iter) == 0)
190 CALL get_iterator_info(ab_iter, ikind=ikind, iatom=iatom, nkind=nkind, r=rab)
191
192 DO kkind = 1, nkind
193 CALL nl_set_sub_iterator(ac_iter, ikind, kkind, iatom)
194
195 DO WHILE (nl_sub_iterate(ac_iter) == 0)
196
197 IF (my_sort_bc) THEN
198 !we check for rbc or rac because of symmetry in ab_nl
199 CALL get_iterator_info(ac_iter, r=rac)
200 rbc(:) = rac(:) - rab(:)
201 IF (.NOT. (all(abs(rbc) <= 1.0e-8_dp) .OR. all(abs(rac) <= 1.0e-8_dp))) cycle
202
203 END IF
204
205 nblk = nblk + 1
206 END DO !ac_iter
207 END DO !kkind
208 END DO !ab_iter
209 CALL neighbor_list_iterator_release(ab_iter)
210 CALL neighbor_list_iterator_release(ac_iter)
211
212 ALLOCATE (idx1(nblk), idx2(nblk), idx3(nblk))
213
214 !actually reserve the blocks
215 CALL neighbor_list_iterator_create(ab_iter, ab_nl)
216 CALL neighbor_list_iterator_create(ac_iter, ac_nl, search=.true.)
217 nblk = 0
218 DO WHILE (neighbor_list_iterate(ab_iter) == 0)
219 CALL get_iterator_info(ab_iter, ikind=ikind, iatom=iatom, jatom=jatom, nkind=nkind, r=rab)
220
221 DO kkind = 1, nkind
222 CALL nl_set_sub_iterator(ac_iter, ikind, kkind, iatom)
223
224 DO WHILE (nl_sub_iterate(ac_iter) == 0)
225 CALL get_iterator_info(ac_iter, jatom=katom, r=rac)
226
227 IF (my_sort_bc) THEN
228 !we check for rbc or rac because of symmetry in ab_nl
229 CALL get_iterator_info(ac_iter, r=rac)
230 rbc(:) = rac(:) - rab(:)
231 IF (.NOT. (all(abs(rbc) <= 1.0e-8_dp) .OR. all(abs(rac) <= 1.0e-8_dp))) cycle
232
233 END IF
234
235 nblk = nblk + 1
236
237 idx1(nblk) = iatom
238 idx2(nblk) = jatom
239 idx3(nblk) = katom
240
241 END DO !ac_iter
242 END DO !kkind
243 END DO !ab_iter
244 CALL neighbor_list_iterator_release(ab_iter)
245 CALL neighbor_list_iterator_release(ac_iter)
246
247!TODO: Parallelize creation of block list.
248!$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 CALL dbt_finalize(pq_x)
255
256 !clean-up
257 CALL dbt_distribution_destroy(t_dist)
258 CALL dbt_pgrid_destroy(t_pgrid)
259
260 CALL timestop(handle)
261
262 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 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 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, lc_max, &
299 lc_min, npgfa, npgfb, npgfc, nsgfa, &
300 nsgfb, nsgfc
301 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 REAL(dp), ALLOCATABLE, DIMENSION(:) :: ccp_buffer, cpp_buffer
307 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: max_contr, max_contra, max_contrb, &
308 max_contrc
309 REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: iabc, sabc, work
310 REAL(dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b, set_radius_c
311 REAL(dp), DIMENSION(:, :), POINTER :: rpgf_a, rpgf_b, rpgf_c, zeta, zetb, zetc
312 TYPE(cp_2d_r_p_type), DIMENSION(:, :), POINTER :: spb, spc, tspa
313 TYPE(cp_libint_t) :: lib
314 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 NULLIFY (basis_set, basis_set_list, para_env, la_max, la_min)
322 NULLIFY (lb_max, lb_min, lc_max, lc_min, npgfa, npgfb, npgfc, nsgfa, nsgfb, nsgfc)
323 NULLIFY (first_sgfa, first_sgfb, first_sgfc, set_radius_a, set_radius_b, set_radius_c)
324 NULLIFY (rpgf_a, rpgf_b, rpgf_c, zeta, zetb, zetc)
325 NULLIFY (basis_set_a, basis_set_b, basis_set_c, tspa, spb, spc)
326
327 CALL timeset(routinen, handle)
328
329 !Need the max l for each basis for libint (and overall max #of sets for screening)
330 nbasis = SIZE(basis_set_list_a)
331 max_nset = 0
332 maxli = 0
333 DO ibasis = 1, nbasis
334 CALL get_gto_basis_set(gto_basis_set=basis_set_list_a(ibasis)%gto_basis_set, &
335 maxl=imax, nset=iset, nsgf_set=nsgfa)
336 maxli = max(maxli, imax)
337 max_nset = max(max_nset, iset)
338 END DO
339 maxlj = 0
340 DO ibasis = 1, nbasis
341 CALL get_gto_basis_set(gto_basis_set=basis_set_list_b(ibasis)%gto_basis_set, &
342 maxl=imax, nset=iset, nsgf_set=nsgfb, npgf=npgfb)
343 maxlj = max(maxlj, imax)
344 max_nset = max(max_nset, iset)
345 END DO
346 maxlk = 0
347 DO ibasis = 1, nbasis
348 CALL get_gto_basis_set(gto_basis_set=basis_set_list_c(ibasis)%gto_basis_set, &
349 maxl=imax, nset=iset, npgf=npgfc)
350 maxlk = max(maxlk, imax)
351 max_nset = max(max_nset, iset)
352 END DO
353 m_max = maxli + maxlj + maxlk
354
355 !Screening
356 do_screen = .false.
357 IF (PRESENT(eps_screen)) THEN
358 do_screen = .true.
359 my_eps_screen = eps_screen
360 END IF
361 screen_radius = 0.0_dp
362 IF (potential_parameter%potential_type == do_potential_truncated .OR. &
363 potential_parameter%potential_type == do_potential_short) THEN
364
365 screen_radius = potential_parameter%cutoff_radius*cutoff_screen_factor
366 ELSE IF (potential_parameter%potential_type == do_potential_coulomb) THEN
367
368 screen_radius = 1000000.0_dp
369 END IF
370
371 !get maximum contraction values for abc_contract screening
372 IF (do_screen) THEN
373
374 !Allocate max_contraction arrays such that we have a specific value for each set/kind
375 ALLOCATE (max_contr(max_nset, nbasis), max_contra(max_nset, nbasis), &
376 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 DO ilist = 1, 3
380
381 IF (ilist == 1) basis_set_list => basis_set_list_a
382 IF (ilist == 2) basis_set_list => basis_set_list_b
383 IF (ilist == 3) basis_set_list => basis_set_list_c
384
385 max_contr = 0.0_dp
386
387 DO ibasis = 1, nbasis
388 basis_set => basis_set_list(ibasis)%gto_basis_set
389
390 DO iset = 1, basis_set%nset
391
392 ncoa = basis_set%npgf(iset)*ncoset(basis_set%lmax(iset))
393 sgfa = basis_set%first_sgf(1, iset)
394 egfa = sgfa + basis_set%nsgf_set(iset) - 1
395
396 max_contr(iset, ibasis) = &
397 maxval([(sum(abs(basis_set%sphi(1:ncoa, i))), i=sgfa, egfa)])
398
399 END DO !iset
400 END DO !ibasis
401
402 IF (ilist == 1) max_contra(:, :) = max_contr(:, :)
403 IF (ilist == 2) max_contrb(:, :) = max_contr(:, :)
404 IF (ilist == 3) max_contrc(:, :) = max_contr(:, :)
405 END DO !ilist
406 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 ALLOCATE (tspa(max_nset, nbasis), spb(max_nset, nbasis), spc(max_nset, nbasis))
412 DO ibasis = 1, nbasis
413 DO iset = 1, max_nset
414 NULLIFY (tspa(iset, ibasis)%array)
415 NULLIFY (spb(iset, ibasis)%array)
416 NULLIFY (spc(iset, ibasis)%array)
417 END DO
418 END DO
419
420 DO ilist = 1, 3
421
422 DO ibasis = 1, nbasis
423 IF (ilist == 1) basis_set => basis_set_list_a(ibasis)%gto_basis_set
424 IF (ilist == 2) basis_set => basis_set_list_b(ibasis)%gto_basis_set
425 IF (ilist == 3) basis_set => basis_set_list_c(ibasis)%gto_basis_set
426
427 DO iset = 1, basis_set%nset
428
429 ncoa = basis_set%npgf(iset)*ncoset(basis_set%lmax(iset))
430 sgfa = basis_set%first_sgf(1, iset)
431 egfa = sgfa + basis_set%nsgf_set(iset) - 1
432
433 IF (ilist == 1) THEN
434 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 basis_set%nsgf_set(iset))
441#else
442 tspa(iset, ibasis)%array(:, :) = transpose(basis_set%sphi(1:ncoa, sgfa:egfa))
443#endif
444 ELSE IF (ilist == 2) THEN
445 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 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 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 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 my_sort_bc = .false.
471 IF (PRESENT(only_bc_same_center)) my_sort_bc = only_bc_same_center
472
473 !Init the truncated Coulomb operator
474 CALL get_qs_env(qs_env, para_env=para_env)
475 IF (potential_parameter%potential_type == do_potential_truncated) THEN
476
477 !open the file only if necessary
478 IF (m_max > get_lmax_init()) THEN
479 IF (para_env%mepos == 0) THEN
480 CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
481 END IF
482 CALL init(m_max, unit_id, para_env%mepos, para_env)
483 IF (para_env%mepos == 0) THEN
484 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 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!$ nthread = omp_get_max_threads()
497
498 ALLOCATE (o3c)
499 CALL init_o3c_container(o3c, 1, basis_set_list_a, basis_set_list_b, basis_set_list_c, &
500 ab_nl, ac_nl, only_bc_same_center=my_sort_bc)
501 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!$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 CALL o3c_iterator_release(o3c_iterator)
650 CALL release_o3c_container(o3c)
651 DEALLOCATE (o3c)
652
653 DO iset = 1, max_nset
654 DO ibasis = 1, nbasis
655 IF (ASSOCIATED(tspa(iset, ibasis)%array)) DEALLOCATE (tspa(iset, ibasis)%array)
656 IF (ASSOCIATED(spb(iset, ibasis)%array)) DEALLOCATE (spb(iset, ibasis)%array)
657 IF (ASSOCIATED(spc(iset, ibasis)%array)) DEALLOCATE (spc(iset, ibasis)%array)
658 END DO
659 END DO
660 DEALLOCATE (tspa, spb, spc)
661
662 CALL timestop(handle)
663
664 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 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 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius
692 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 TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
697 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
698 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
699
700 NULLIFY (atomic_kind_set, distribution_1d, distribution_2d, molecule_set, particle_set, cell)
701
702! Initialization
703 CALL get_qs_env(qs_env, nkind=nkind)
704 CALL section_vals_val_get(qs_env%input, "DFT%SUBCELLS", r_val=subcells)
705
706 my_restrictb = .false.
707 IF (PRESENT(excited_atoms)) THEN
708 my_restrictb = .true.
709 END IF
710
711 ALLOCATE (a_present(nkind), b_present(nkind))
712 a_present = .false.
713 b_present = .false.
714 ALLOCATE (a_radius(nkind), b_radius(nkind))
715 a_radius = 0.0_dp
716 b_radius = 0.0_dp
717
718! Set up the radii
719 DO ikind = 1, nkind
720 IF (ASSOCIATED(basis_a(ikind)%gto_basis_set)) THEN
721 a_present(ikind) = .true.
722 CALL get_gto_basis_set(basis_a(ikind)%gto_basis_set, kind_radius=a_radius(ikind))
723 END IF
724
725 IF (ASSOCIATED(basis_b(ikind)%gto_basis_set)) THEN
726 b_present(ikind) = .true.
727 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 ALLOCATE (pair_radius(nkind, nkind))
732 pair_radius = 0.0_dp
733 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 particle_set=particle_set, molecule_set=molecule_set)
739
740 !use an external distribution_2d if required
741 IF (PRESENT(ext_dist2d)) distribution_2d => ext_dist2d
742
743 ALLOCATE (atom2d(nkind))
744 CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
745 molecule_set, .false., particle_set)
746
747 IF (my_restrictb) THEN
748
749 CALL build_neighbor_lists(ab_list, particle_set, atom2d, cell, pair_radius, subcells, &
750 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 nlname="XAS_TDP_ovlp_nl")
756
757 END IF
758! Clean-up
759 CALL atom2d_cleanup(atom2d)
760
761 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 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 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius
794 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 TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
799 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
800 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
801
802 NULLIFY (atomic_kind_set, distribution_1d, distribution_2d, molecule_set, particle_set, cell)
803
804! Initialization
805 CALL get_qs_env(qs_env, nkind=nkind)
806 CALL section_vals_val_get(qs_env%input, "DFT%SUBCELLS", r_val=subcells)
807 sort_atoms = .false.
808 IF ((PRESENT(excited_atoms))) sort_atoms = .true.
809
810 ALLOCATE (a_present(nkind), c_present(nkind))
811 a_present = .false.
812 c_present = .false.
813 ALLOCATE (a_radius(nkind), c_radius(nkind))
814 a_radius = 0.0_dp
815 c_radius = 0.0_dp
816
817! Set up the radii, depending on the operator type
818 IF (op_type == do_potential_id) THEN
819
820 !overlap => use the kind radius for both a and c
821 DO ikind = 1, nkind
822 !orbital basis set
823 IF (ASSOCIATED(basis_a(ikind)%gto_basis_set)) THEN
824 a_present(ikind) = .true.
825 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 IF (ASSOCIATED(basis_c(ikind)%gto_basis_set)) THEN
829 c_present(ikind) = .true.
830 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 ELSE IF (op_type == do_potential_coulomb) THEN
835
836 !Coulomb operator, virtually infinite range => set c_radius to arbitrarily large number
837 DO ikind = 1, nkind
838 IF (ASSOCIATED(basis_c(ikind)%gto_basis_set)) THEN
839 c_present(ikind) = .true.
840 c_radius(ikind) = 1000000.0_dp
841 END IF
842 IF (ASSOCIATED(basis_a(ikind)%gto_basis_set)) a_present(ikind) = .true.
843 END DO !ikind
844
845 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 DO ikind = 1, nkind
849 IF (ASSOCIATED(basis_a(ikind)%gto_basis_set)) THEN
850 a_present(ikind) = .true.
851 CALL get_gto_basis_set(basis_a(ikind)%gto_basis_set, kind_radius=a_radius(ikind))
852 END IF
853 IF (ASSOCIATED(basis_c(ikind)%gto_basis_set)) THEN
854 c_present(ikind) = .true.
855 CALL get_gto_basis_set(basis_c(ikind)%gto_basis_set, kind_radius=c_radius(ikind))
856 c_radius(ikind) = c_radius(ikind) + x_range
857 END IF
858 END DO !ikind
859
860 ELSE
861 cpabort("Operator not known")
862 END IF
863
864 ALLOCATE (pair_radius(nkind, nkind))
865 pair_radius = 0.0_dp
866 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 particle_set=particle_set, molecule_set=molecule_set)
872
873 !use an external distribution_2d if required
874 IF (PRESENT(ext_dist2d)) distribution_2d => ext_dist2d
875
876 ALLOCATE (atom2d(nkind))
877 CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
878 molecule_set, .false., particle_set)
879
880 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 nlname="XAS_TDP_3c_nl")
884 ELSE
885 CALL build_neighbor_lists(ac_list, particle_set, atom2d, cell, pair_radius, subcells, &
886 operator_type="ABC", nlname="XAS_TDP_3c_nl")
887 END IF
888
889! Clean-up
890 CALL atom2d_cleanup(atom2d)
891
892 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: nparticle_local_col, nparticle_local_row
923 INTEGER, DIMENSION(:, :), POINTER :: col_dist, row_dist
924 LOGICAL :: my_sort_bc
925 REAL(dp) :: cost, rab(3), rac(3), rbc(3)
926 REAL(dp), ALLOCATABLE, DIMENSION(:) :: col_cost, col_proc_cost, row_cost, &
927 row_proc_cost
928 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 DIMENSION(:), POINTER :: ab_iter, ac_iter
933 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
934 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
935
936 NULLIFY (para_env, col_dist, row_dist, blacs_env, qs_kind_set, particle_set)
937 NULLIFY (local_particle_col, local_particle_row, ab_iter, ac_iter)
938
939 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 my_sort_bc = .false.
946 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 qs_kind_set=qs_kind_set, particle_set=particle_set, nkind=nkind)
950
951 myprow = blacs_env%mepos(1) + 1
952 mypcol = blacs_env%mepos(2) + 1
953 nprow = blacs_env%num_pe(1)
954 npcol = blacs_env%num_pe(2)
955
956 ALLOCATE (col_cost(natom), row_cost(natom))
957 col_cost = 0.0_dp; row_cost = 0.0_dp
958
959 ALLOCATE (row_dist(natom, 2), col_dist(natom, 2))
960 row_dist = 1; col_dist = 1
961
962 CALL neighbor_list_iterator_create(ab_iter, ab_list)
963 CALL neighbor_list_iterator_create(ac_iter, ac_list, search=.true.)
964 DO WHILE (neighbor_list_iterate(ab_iter) == 0)
965 CALL get_iterator_info(ab_iter, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rab)
966
967 DO kkind = 1, nkind
968 CALL nl_set_sub_iterator(ac_iter, ikind, kkind, iatom)
969
970 DO WHILE (nl_sub_iterate(ac_iter) == 0)
971
972 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 CALL get_iterator_info(ac_iter, r=rac)
975 rbc(:) = rac(:) - rab(:)
976 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 nsgfa = basis_set_a(ikind)%gto_basis_set%nsgf
982 nsgfb = basis_set_b(jkind)%gto_basis_set%nsgf
983 nsgfc = basis_set_c(kkind)%gto_basis_set%nsgf
984
985 cost = real(nsgfa*nsgfb*nsgfc, dp)
986
987 row_cost(iatom) = row_cost(iatom) + cost
988 col_cost(jatom) = col_cost(jatom) + cost
989
990 END DO !ac_iter
991 END DO !kkind
992 END DO !ab_iter
993 CALL neighbor_list_iterator_release(ab_iter)
994 CALL neighbor_list_iterator_release(ac_iter)
995
996 CALL para_env%sum(row_cost)
997 CALL para_env%sum(col_cost)
998
999 !Distribute the cost as evenly as possible
1000 ALLOCATE (col_proc_cost(npcol), row_proc_cost(nprow))
1001 col_proc_cost = 0.0_dp; row_proc_cost = 0.0_dp
1002 DO i = 1, natom
1003 iatom = maxloc(row_cost, 1)
1004 ip = minloc(row_proc_cost, 1)
1005 row_proc_cost(ip) = row_proc_cost(ip) + row_cost(iatom)
1006 row_dist(iatom, 1) = ip
1007 row_cost(iatom) = 0.0_dp
1008
1009 iatom = maxloc(col_cost, 1)
1010 ip = minloc(col_proc_cost, 1)
1011 col_proc_cost(ip) = col_proc_cost(ip) + col_cost(iatom)
1012 col_dist(iatom, 1) = ip
1013 col_cost(iatom) = 0.0_dp
1014 END DO
1015
1016 !the usual stuff
1017 ALLOCATE (local_particle_col(nkind), local_particle_row(nkind))
1018 ALLOCATE (nparticle_local_row(nkind), nparticle_local_col(nkind))
1019 nparticle_local_row = 0; nparticle_local_col = 0
1020
1021 DO iatom = 1, natom
1022 ikind = particle_set(iatom)%atomic_kind%kind_number
1023
1024 IF (row_dist(iatom, 1) == myprow) nparticle_local_row(ikind) = nparticle_local_row(ikind) + 1
1025 IF (col_dist(iatom, 1) == mypcol) nparticle_local_col(ikind) = nparticle_local_col(ikind) + 1
1026 END DO
1027
1028 DO ikind = 1, nkind
1029 n = nparticle_local_row(ikind)
1030 ALLOCATE (local_particle_row(ikind)%array(n))
1031
1032 n = nparticle_local_col(ikind)
1033 ALLOCATE (local_particle_col(ikind)%array(n))
1034 END DO
1035
1036 nparticle_local_row = 0; nparticle_local_col = 0
1037 DO iatom = 1, natom
1038 ikind = particle_set(iatom)%atomic_kind%kind_number
1039
1040 IF (row_dist(iatom, 1) == myprow) THEN
1041 nparticle_local_row(ikind) = nparticle_local_row(ikind) + 1
1042 local_particle_row(ikind)%array(nparticle_local_row(ikind)) = iatom
1043 END IF
1044 IF (col_dist(iatom, 1) == mypcol) THEN
1045 nparticle_local_col(ikind) = nparticle_local_col(ikind) + 1
1046 local_particle_col(ikind)%array(nparticle_local_col(ikind)) = iatom
1047 END IF
1048 END DO
1049
1050 !Finally create the dist_2d
1051 row_dist(:, 1) = row_dist(:, 1) - 1
1052 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 local_cols_ptr=local_particle_col, blacs_env=blacs_env)
1056
1057 CALL timestop(handle)
1058
1059 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 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 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 POINTER :: ab_list, ac_list
1087 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1088 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1089
1090 NULLIFY (basis_set_ri, basis_set_orb, ac_list, ab_list, qs_kind_set, para_env, particle_set)
1091
1092 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 natom=natom, particle_set=particle_set)
1097
1098! Build the basis set lists
1099 ALLOCATE (basis_set_ri(nkind))
1100 ALLOCATE (basis_set_orb(nkind))
1101 CALL basis_set_list_setup(basis_set_ri, "RI_XAS", qs_kind_set)
1102 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 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 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 basis_set_orb, basis_set_ri, qs_env)
1112 CALL release_neighbor_list_sets(ab_list)
1113 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 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 ext_dist2d=xas_tdp_env%opt_dist2d_ex)
1122
1123! Allocate, init and compute the integrals.
1124 ALLOCATE (blk_size_orb(natom), blk_size_ri(natom))
1125 CALL cp_dbcsr_dist2d_to_dist(xas_tdp_env%opt_dist2d_ex, opt_dbcsr_dist)
1126 CALL get_particle_set(particle_set, qs_kind_set, nsgf=blk_size_orb, basis=basis_set_orb)
1127 CALL get_particle_set(particle_set, qs_kind_set, nsgf=blk_size_ri, basis=basis_set_ri)
1128
1129 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 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 eps_screen=xas_tdp_control%eps_screen)
1135
1136! Clean-up
1137 CALL release_neighbor_list_sets(ab_list)
1138 CALL release_neighbor_list_sets(ac_list)
1139 CALL dbcsr_distribution_release(opt_dbcsr_dist)
1140 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 CALL para_env%sync()
1145
1146 CALL timestop(handle)
1147
1148 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 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 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 POINTER :: ab_list, ac_list
1173 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1174 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1175
1176 NULLIFY (basis_set_ri, basis_set_orb, ac_list, ab_list, qs_kind_set, para_env, particle_set)
1177
1178 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 natom=natom, particle_set=particle_set)
1183
1184! Build the basis set lists
1185 ALLOCATE (basis_set_ri(nkind))
1186 ALLOCATE (basis_set_orb(nkind))
1187 CALL basis_set_list_setup(basis_set_ri, "RI_XAS", qs_kind_set)
1188 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 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 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 basis_set_orb, basis_set_ri, qs_env, only_bc_same_center=.true.)
1197 CALL release_neighbor_list_sets(ab_list)
1198 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 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 ext_dist2d=xas_tdp_env%opt_dist2d_coul)
1212
1213! Allocate, init and compute the integrals
1214 ALLOCATE (blk_size_orb(natom), blk_size_ri(natom))
1215 CALL cp_dbcsr_dist2d_to_dist(xas_tdp_env%opt_dist2d_coul, opt_dbcsr_dist)
1216 CALL get_particle_set(particle_set, qs_kind_set, nsgf=blk_size_orb, basis=basis_set_orb)
1217 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 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 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 basis_set_ri, pot, qs_env, only_bc_same_center=.true.)
1225
1226! Clean-up
1227 CALL release_neighbor_list_sets(ab_list)
1228 CALL release_neighbor_list_sets(ac_list)
1229 CALL dbcsr_distribution_release(opt_dbcsr_dist)
1230 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 CALL para_env%sync()
1235
1236 CALL timestop(handle)
1237
1238 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 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 INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf_set, nsgf_set
1263 INTEGER, DIMENSION(:, :), POINTER :: first_sgf
1264 REAL(dp) :: r(3)
1265 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: metric, pq, work
1266 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 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1271
1272 NULLIFY (ri_basis, qs_kind_set, para_env, lmin, lmax, npgf_set, zet, rpgf, first_sgf)
1273 NULLIFY (sphi, nsgf_set)
1274
1275! Initialization
1276 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, para_env=para_env)
1277 IF (ASSOCIATED(xas_tdp_env%ri_inv_ex)) THEN
1278 DEALLOCATE (xas_tdp_env%ri_inv_ex)
1279 END IF
1280
1281! Get the RI basis of interest and its quantum numbers
1282 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 nsgf_set=nsgf_set, sphi=sphi, nset=nset)
1286 ALLOCATE (metric(nsgf, nsgf))
1287 metric = 0.0_dp
1288
1289 r = 0.0_dp
1290
1291 !init the libint 2-center object
1292 CALL cp_libint_init_2eri(lib, maxl)
1293 CALL cp_libint_set_contrdepth(lib, 1)
1294
1295 !make sure the truncted coulomb is initialized
1296 IF (xas_tdp_control%ri_m_potential%potential_type == do_potential_truncated) THEN
1297
1298 IF (2*maxl + 1 > get_lmax_init()) THEN
1299 IF (para_env%mepos == 0) THEN
1300 CALL open_file(unit_number=unit_id, file_name=xas_tdp_control%ri_m_potential%filename)
1301 END IF
1302 CALL init(2*maxl + 1, unit_id, para_env%mepos, para_env)
1303 IF (para_env%mepos == 0) THEN
1304 CALL close_file(unit_id)
1305 END IF
1306 END IF
1307 END IF
1308
1309! Compute the RI metric
1310 DO pset = 1, nset
1311 ncop = npgf_set(pset)*ncoset(lmax(pset))
1312 sgfp = first_sgf(1, pset)
1313 egfp = sgfp + nsgf_set(pset) - 1
1314
1315 DO qset = 1, nset
1316 ncoq = npgf_set(qset)*ncoset(lmax(qset))
1317 sgfq = first_sgf(1, qset)
1318 egfq = sgfq + nsgf_set(qset) - 1
1319
1320 ALLOCATE (work(ncop, ncoq))
1321 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 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 ncop, ncoq, nsgf_set(pset), nsgf_set(qset))
1329
1330 DEALLOCATE (work)
1331 END DO !qset
1332 END DO !pset
1333
1334! Inverting (to M^-1)
1335 CALL invmat_symm(metric)
1336
1337 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 ALLOCATE (xas_tdp_env%ri_inv_ex(nsgf, nsgf))
1341 xas_tdp_env%ri_inv_ex(:, :) = metric(:, :)
1342 CALL cp_libint_cleanup_2eri(lib)
1343 RETURN
1344
1345 END IF
1346
1347 !make sure the truncted coulomb is initialized
1348 IF (xas_tdp_control%x_potential%potential_type == do_potential_truncated) THEN
1349
1350 IF (2*maxl + 1 > get_lmax_init()) THEN
1351 IF (para_env%mepos == 0) THEN
1352 CALL open_file(unit_number=unit_id, file_name=xas_tdp_control%x_potential%filename)
1353 END IF
1354 CALL init(2*maxl + 1, unit_id, para_env%mepos, para_env)
1355 IF (para_env%mepos == 0) THEN
1356 CALL close_file(unit_id)
1357 END IF
1358 END IF
1359 END IF
1360
1361! Compute the proper exchange 2-center
1362 ALLOCATE (pq(nsgf, nsgf))
1363 pq = 0.0_dp
1364
1365 DO pset = 1, nset
1366 ncop = npgf_set(pset)*ncoset(lmax(pset))
1367 sgfp = first_sgf(1, pset)
1368 egfp = sgfp + nsgf_set(pset) - 1
1369
1370 DO qset = 1, nset
1371 ncoq = npgf_set(qset)*ncoset(lmax(qset))
1372 sgfq = first_sgf(1, qset)
1373 egfq = sgfq + nsgf_set(qset) - 1
1374
1375 ALLOCATE (work(ncop, ncoq))
1376 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 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 ncop, ncoq, nsgf_set(pset), nsgf_set(qset))
1384
1385 DEALLOCATE (work)
1386 END DO !qset
1387 END DO !pset
1388
1389! Compute and store M^-1 (P|Q) M^-1
1390 ALLOCATE (xas_tdp_env%ri_inv_ex(nsgf, nsgf))
1391 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 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 0.0_dp, pq, nsgf)
1397 xas_tdp_env%ri_inv_ex(:, :) = pq(:, :)
1398
1399 CALL cp_libint_cleanup_2eri(lib)
1400
1401 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 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 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1422
1423 NULLIFY (ri_basis, qs_kind_set)
1424
1425! Initialization
1426 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
1427 IF (ASSOCIATED(xas_tdp_env%ri_inv_coul)) THEN
1428 DEALLOCATE (xas_tdp_env%ri_inv_coul)
1429 END IF
1430
1431! Get the RI basis of interest and its quantum numbers
1432 CALL get_qs_kind(qs_kind_set(ex_kind), basis_set=ri_basis, basis_type="RI_XAS")
1433 CALL get_gto_basis_set(ri_basis, nsgf=nsgf)
1434 ALLOCATE (xas_tdp_env%ri_inv_coul(nsgf, nsgf))
1435 xas_tdp_env%ri_inv_coul = 0.0_dp
1436
1437 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 calculate_forces=.false.)
1441 cpassert(ASSOCIATED(xas_tdp_control))
1442 ELSE
1443 CALL periodic_ri_coulomb2(xas_tdp_env%ri_inv_coul, ri_basis, qs_env)
1444 END IF
1445
1446! Inverting
1447 CALL invmat_symm(xas_tdp_env%ri_inv_coul)
1448
1449 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 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 INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nsgf
1466 INTEGER, DIMENSION(:, :), POINTER :: first_sgf
1467 REAL(dp) :: r(3)
1468 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: hpq
1469 REAL(dp), DIMENSION(:, :), POINTER :: sphi, zet
1470 TYPE(cell_type), POINTER :: cell
1471 TYPE(cp_eri_mme_param) :: mme_param
1472 TYPE(mp_para_env_type), POINTER :: para_env
1473 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1474
1475 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 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 do_error_est=.false.)
1487 mme_param%do_calib = .true.
1488
1489 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 nsgf_set=nsgf, sphi=sphi, first_sgf=first_sgf, maxco=maxco)
1493
1494 r = 0.0_dp
1495 ALLOCATE (hpq(nset*maxco, nset*maxco))
1496 hpq = 0.0_dp
1497
1498 DO pset = 1, nset
1499 ncop = npgf(pset)*ncoset(lmax(pset))
1500 sgfp = first_sgf(1, pset)
1501
1502 DO qset = 1, nset
1503 ncoq = npgf(qset)*ncoset(lmax(qset))
1504 sgfq = first_sgf(1, qset)
1505
1506 DO ppgf = 1, npgf(pset)
1507 op = (pset - 1)*maxco + (ppgf - 1)*ncoset(lmax(pset))
1508 DO qpgf = 1, npgf(qset)
1509 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 op, oq)
1514
1515 END DO !qpgf
1516 END DO ! ppgf
1517
1518 !contraction into sgfs
1519 op = (pset - 1)*maxco + 1
1520 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 ncop, ncoq, nsgf(pset), nsgf(qset))
1525
1526 END DO !qset
1527 END DO !pset
1528
1529 !celan-up
1530 CALL eri_mme_release(mme_param%par)
1531
1532 END SUBROUTINE periodic_ri_coulomb2
1533
1534END MODULE xas_tdp_integrals
imax
static int imax(int x, int y)
Returns the larger of two given integers (missing from the C standard)
Definition dbm_internal.h:14
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:225
idx2
#define idx2(a, i, j)
Definition grid_dgemm_tensor_local.h:258
idx3
#define idx3(a, i, j, k)
Definition grid_dgemm_tensor_local.h:257
cp_eri_mme_interface::cp_eri_mme_set_params
Definition cp_eri_mme_interface.F:74
qs_neighbor_list_types::nl_sub_iterate
Definition qs_neighbor_list_types.F:96
ai_contraction_sphi
Contraction of integrals over primitive Cartesian Gaussians based on the contraction matrix sphi whic...
Definition ai_contraction_sphi.F:15
ai_contraction_sphi::abc_contract_xsmm
subroutine, public abc_contract_xsmm(abcint, sabc, sphi_a, sphi_b, sphi_c, ncoa, ncob, ncoc, nsgfa, nsgfb, nsgfc, cpp_buffer, ccp_buffer, prefac, pstfac)
3-center contraction routine from primitive cartesian Gaussians to spherical Gaussian functions using...
Definition ai_contraction_sphi.F:250
ai_contraction_sphi::ab_contract
subroutine, public ab_contract(abint, sab, sphi_a, sphi_b, ncoa, ncob, nsgfa, nsgfb)
contract overlap integrals (a,b) and transfer to spherical Gaussians
Definition ai_contraction_sphi.F:47
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
basis_set_types
Definition basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum, ccon)
...
Definition basis_set_types.F:643
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
constants_operator
constants for the different operators of the 2c-integrals
Definition constants_operator.F:14
constants_operator::operator_coulomb
integer, parameter, public operator_coulomb
Definition constants_operator.F:22
cp_array_utils
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Definition cp_array_utils.F:18
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_distribution_release
subroutine, public dbcsr_distribution_release(dist)
...
Definition cp_dbcsr_api.F:642
cp_dbcsr_api::dbcsr_distribution_get
subroutine, public dbcsr_distribution_get(dist, row_dist, col_dist, nrows, ncols, has_threads, group, mynode, numnodes, nprows, npcols, myprow, mypcol, pgrid, subgroups_defined, prow_group, pcol_group)
...
Definition cp_dbcsr_api.F:571
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_dist2d_to_dist
subroutine, public cp_dbcsr_dist2d_to_dist(dist2d, dist)
Creates a DBCSR distribution from a distribution_2d.
Definition cp_dbcsr_operations.F:415
cp_eri_mme_interface
Interface to Minimax-Ewald method for periodic ERI's to be used in CP2K.
Definition cp_eri_mme_interface.F:12
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition cp_files.F:311
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition cp_files.F:122
dbt_api
This is the start of a dbt_api, all publically needed functions are exported here....
Definition dbt_api.F:17
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition distribution_1d_types.F:24
distribution_2d_types
stores a mapping of 2D info (e.g. matrix) on a 2D processor distribution (i.e. blacs grid) where cpus...
Definition distribution_2d_types.F:15
distribution_2d_types::distribution_2d_create
subroutine, public distribution_2d_create(distribution_2d, blacs_env, local_rows_ptr, n_local_rows, local_cols_ptr, row_distribution_ptr, col_distribution_ptr, n_local_cols, n_row_distribution, n_col_distribution)
initializes the distribution_2d
Definition distribution_2d_types.F:102
eri_mme_integrate
Minimax-Ewald (MME) method for calculating 2-center and 3-center electron repulsion integrals (ERI) o...
Definition eri_mme_integrate.F:23
eri_mme_integrate::eri_mme_2c_integrate
subroutine, public eri_mme_2c_integrate(param, la_min, la_max, lb_min, lb_max, zeta, zetb, rab, hab, o1, o2, g_count, r_count, normalize, potential, pot_par)
Low-level integration routine for 2-center ERIs.
Definition eri_mme_integrate.F:72
eri_mme_types
Types and initialization / release routines for Minimax-Ewald method for electron repulsion integrals...
Definition eri_mme_types.F:16
eri_mme_types::eri_mme_release
subroutine, public eri_mme_release(param)
...
Definition eri_mme_types.F:471
eri_mme_types::eri_mme_init
subroutine, public eri_mme_init(param, n_minimax, cutoff, do_calib_cutoff, do_error_est, cutoff_min, cutoff_max, cutoff_eps, cutoff_delta, sum_precision, debug, debug_delta, debug_nsum, unit_nr, print_calib)
...
Definition eri_mme_types.F:116
gamma
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition gamma.F:15
gamma::init_md_ftable
subroutine, public init_md_ftable(nmax)
Initialize a table of F_n(t) values in the range 0 <= t <= 12 with a stepsize of 0....
Definition gamma.F:540
generic_os_integrals
Calculation of contracted, spherical Gaussian integrals using the (OS) integral scheme....
Definition generic_os_integrals.F:18
generic_os_integrals::int_operators_r12_ab_os
subroutine, public int_operators_r12_ab_os(r12_operator, vab, dvab, rab, fba, fbb, omega, r_cutoff, calculate_forces)
Calcululates the two-center integrals of the type (a|O(r12)|b) using the OS scheme.
Definition generic_os_integrals.F:78
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::do_potential_truncated
integer, parameter, public do_potential_truncated
Definition input_constants.F:834
input_constants::do_potential_id
integer, parameter, public do_potential_id
Definition input_constants.F:834
input_constants::do_potential_coulomb
integer, parameter, public do_potential_coulomb
Definition input_constants.F:834
input_constants::do_potential_short
integer, parameter, public do_potential_short
Definition input_constants.F:834
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1047
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
libint_2c_3c
2- and 3-center electron repulsion integral routines based on libint2 Currently available operators: ...
Definition libint_2c_3c.F:14
libint_2c_3c::eri_2center
subroutine, public eri_2center(int_ab, la_min, la_max, npgfa, zeta, rpgfa, ra, lb_min, lb_max, npgfb, zetb, rpgfb, rb, dab, lib, potential_parameter)
Computes the 2-center electron repulsion integrals (a|b) for a given set of cartesian gaussian orbita...
Definition libint_2c_3c.F:862
libint_2c_3c::eri_3center
subroutine, public eri_3center(int_abc, la_min, la_max, npgfa, zeta, rpgfa, ra, lb_min, lb_max, npgfb, zetb, rpgfb, rb, lc_min, lc_max, npgfc, zetc, rpgfc, rc, dab, dac, dbc, lib, potential_parameter, int_abc_ext)
Computes the 3-center electron repulsion integrals (ab|c) for a given set of cartesian gaussian orbit...
Definition libint_2c_3c.F:115
libint_2c_3c::cutoff_screen_factor
real(kind=dp), parameter, public cutoff_screen_factor
Definition libint_2c_3c.F:49
libint_wrapper
Interface to the Libint-Library or a c++ wrapper.
Definition libint_wrapper.F:15
libint_wrapper::cp_libint_init_2eri
subroutine, public cp_libint_init_2eri(lib, max_am)
Definition libint_wrapper.F:1175
libint_wrapper::cp_libint_init_3eri
subroutine, public cp_libint_init_3eri(lib, max_am)
Definition libint_wrapper.F:1139
libint_wrapper::cp_libint_cleanup_3eri
subroutine, public cp_libint_cleanup_3eri(lib)
Definition libint_wrapper.F:1207
libint_wrapper::cp_libint_set_contrdepth
subroutine, public cp_libint_set_contrdepth(lib, contrdepth)
Definition libint_wrapper.F:1109
libint_wrapper::cp_libint_cleanup_2eri
subroutine, public cp_libint_cleanup_2eri(lib)
Definition libint_wrapper.F:1237
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::invmat_symm
subroutine, public invmat_symm(a, potrf, uplo)
returns inverse of real symmetric, positive definite matrix
Definition mathlib.F:587
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
molecule_types
Define the data structure for the molecule information.
Definition molecule_types.F:16
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
particle_methods
Define methods related to particle_type.
Definition particle_methods.F:14
particle_methods::get_particle_set
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis, ncgf)
Get the components of a particle set.
Definition particle_methods.F:121
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:530
qs_integral_utils
Some utility functions for the calculation of integrals.
Definition qs_integral_utils.F:14
qs_integral_utils::basis_set_list_setup
subroutine, public basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
Set up an easy accessible list of the basis sets for all kinds.
Definition qs_integral_utils.F:149
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:468
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_neighbor_list_types::release_neighbor_list_sets
subroutine, public release_neighbor_list_sets(nlists)
releases an array of neighbor_list_sets
Definition qs_neighbor_list_types.F:1095
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition qs_neighbor_list_types.F:165
qs_neighbor_list_types::nl_set_sub_iterator
subroutine, public nl_set_sub_iterator(iterator_set, ikind, jkind, iatom, mepos)
...
Definition qs_neighbor_list_types.F:288
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition qs_neighbor_list_types.F:251
qs_neighbor_list_types::get_neighbor_list_set_p
subroutine, public get_neighbor_list_set_p(neighbor_list_sets, nlist, symmetric)
Return the components of the first neighbor list set.
Definition qs_neighbor_list_types.F:995
qs_neighbor_list_types::neighbor_list_iterate
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
Definition qs_neighbor_list_types.F:351
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition qs_neighbor_list_types.F:549
qs_neighbor_lists
Generate the atomic neighbor lists.
Definition qs_neighbor_lists.F:19
qs_neighbor_lists::atom2d_cleanup
subroutine, public atom2d_cleanup(atom2d)
free the internals of atom2d
Definition qs_neighbor_lists.F:132
qs_neighbor_lists::pair_radius_setup
subroutine, public pair_radius_setup(present_a, present_b, radius_a, radius_b, pair_radius, prmin)
...
Definition qs_neighbor_lists.F:1665
qs_neighbor_lists::build_neighbor_lists
subroutine, public build_neighbor_lists(ab_list, particle_set, atom, cell, pair_radius, subcells, mic, symmetric, molecular, subset_of_mol, current_subset, operator_type, nlname, atomb_to_keep)
Build simple pair neighbor lists.
Definition qs_neighbor_lists.F:1067
qs_neighbor_lists::atom2d_build
subroutine, public atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, molecule_set, molecule_only, particle_set)
Build some distribution structure of atoms, refactored from build_qs_neighbor_lists.
Definition qs_neighbor_lists.F:174
qs_o3c_types
3-center overlap type integrals containers
Definition qs_o3c_types.F:12
qs_o3c_types::get_o3c_iterator_info
subroutine, public get_o3c_iterator_info(o3c_iterator, mepos, iatom, jatom, katom, ikind, jkind, kkind, rij, rik, cellj, cellk, integral, tvec, force_i, force_j, force_k)
...
Definition qs_o3c_types.F:418
qs_o3c_types::o3c_iterator_create
subroutine, public o3c_iterator_create(o3c, o3c_iterator, nthread)
...
Definition qs_o3c_types.F:357
qs_o3c_types::release_o3c_container
subroutine, public release_o3c_container(o3c_container)
...
Definition qs_o3c_types.F:225
qs_o3c_types::o3c_iterator_release
subroutine, public o3c_iterator_release(o3c_iterator)
...
Definition qs_o3c_types.F:384
qs_o3c_types::o3c_iterate
integer function, public o3c_iterate(o3c_iterator, mepos)
...
Definition qs_o3c_types.F:473
qs_o3c_types::init_o3c_container
subroutine, public init_o3c_container(o3c, nspin, basis_set_list_a, basis_set_list_b, basis_set_list_c, sab_nl, sac_nl, only_bc_same_center)
...
Definition qs_o3c_types.F:117
t_c_g0
This module computes the basic integrals for the truncated coulomb operator.
Definition t_c_g0.F:58
t_c_g0::init
subroutine, public init(nder, iunit, mepos, group)
...
Definition t_c_g0.F:1361
t_c_g0::get_lmax_init
integer function, public get_lmax_init()
Returns the value of nderiv_init so that one can check if opening the potential file is worhtwhile.
Definition t_c_g0.F:1468
xas_tdp_integrals
3-center integrals machinery for the XAS_TDP method
Definition xas_tdp_integrals.F:13
xas_tdp_integrals::build_xas_tdp_3c_nl
subroutine, public build_xas_tdp_3c_nl(ac_list, basis_a, basis_c, op_type, qs_env, excited_atoms, x_range, ext_dist2d)
Builds a neighbor lists set taylored for 3-center integral within XAS TDP, such that only excited ato...
Definition xas_tdp_integrals.F:778
xas_tdp_integrals::compute_ri_coulomb2_int
subroutine, public compute_ri_coulomb2_int(ex_kind, xas_tdp_env, xas_tdp_control, qs_env)
Computes the two-center Coulomb integral needed for the RI in kernel calculation. Stores the integral...
Definition xas_tdp_integrals.F:1413
xas_tdp_integrals::build_xas_tdp_ovlp_nl
subroutine, public build_xas_tdp_ovlp_nl(ab_list, basis_a, basis_b, qs_env, excited_atoms, ext_dist2d)
Builds a neighbor lists set for overlaping 2-center S_ab, where b is restricted on a a given list of ...
Definition xas_tdp_integrals.F:678
xas_tdp_integrals::fill_pqx_tensor
subroutine, public fill_pqx_tensor(pq_x, ab_nl, ac_nl, basis_set_list_a, basis_set_list_b, basis_set_list_c, potential_parameter, qs_env, only_bc_same_center, eps_screen)
Fills the given 3 dimensional (pq|X) tensor with integrals.
Definition xas_tdp_integrals.F:282
xas_tdp_integrals::compute_ri_3c_coulomb
subroutine, public compute_ri_3c_coulomb(xas_tdp_env, qs_env)
Computes the RI Coulomb 3-center integrals (ab|c), where c is from the RI_XAS basis and centered on t...
Definition xas_tdp_integrals.F:1159
xas_tdp_integrals::compute_ri_exchange2_int
subroutine, public compute_ri_exchange2_int(ex_kind, xas_tdp_env, xas_tdp_control, qs_env)
Computes the two-center Exchange integral needed for the RI in kernel calculation....
Definition xas_tdp_integrals.F:1254
xas_tdp_integrals::create_pqx_tensor
subroutine, public create_pqx_tensor(pq_x, ab_nl, ac_nl, matrix_dist, blk_size_1, blk_size_2, blk_size_3, only_bc_same_center)
Prepares a tensor to hold 3-center integrals (pq|X), where p,q are distributed according to the given...
Definition xas_tdp_integrals.F:133
xas_tdp_integrals::get_opt_3c_dist2d
subroutine, public get_opt_3c_dist2d(opt_3c_dist2d, ab_list, ac_list, basis_set_a, basis_set_b, basis_set_c, qs_env, only_bc_same_center)
Returns an optimized distribution_2d for the given neighbor lists based on an evaluation of the cost ...
Definition xas_tdp_integrals.F:908
xas_tdp_integrals::compute_ri_3c_exchange
subroutine, public compute_ri_3c_exchange(ex_atoms, xas_tdp_env, xas_tdp_control, qs_env)
Computes the RI exchange 3-center integrals (ab|c), where c is from the RI_XAS basis and centered on ...
Definition xas_tdp_integrals.F:1072
xas_tdp_types
Define XAS TDP control type and associated create, release, etc subroutines, as well as XAS TDP envir...
Definition xas_tdp_types.F:13
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_p_type
Definition basis_set_types.F:94
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_array_utils::cp_1d_i_p_type
represent a pointer to a 1d array
Definition cp_array_utils.F:221
cp_array_utils::cp_2d_r_p_type
represent a pointer to a 2d array
Definition cp_array_utils.F:109
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_dbcsr_api::dbcsr_distribution_type
Definition cp_dbcsr_api.F:182
cp_eri_mme_interface::cp_eri_mme_param
Definition cp_eri_mme_interface.F:79
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
distribution_2d_types::distribution_2d_type
distributes pairs on a 2d grid of processors
Definition distribution_2d_types.F:62
libint_2c_3c::libint_potential_type
Definition libint_2c_3c.F:67
libint_wrapper::cp_libint_t
Definition libint_wrapper.F:71
message_passing::mp_comm_type
Definition message_passing.F:155
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:743
molecule_types::molecule_type
Definition molecule_types.F:159
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:220
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:177
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_neighbor_lists::local_atoms_type
Definition qs_neighbor_lists.F:105
qs_o3c_types::o3c_container_type
Definition qs_o3c_types.F:56
qs_o3c_types::o3c_iterator_type
Definition qs_o3c_types.F:75
xas_tdp_types::xas_tdp_control_type
Type containing control information for TDP XAS calculations.
Definition xas_tdp_types.F:102
xas_tdp_types::xas_tdp_env_type
Type containing informations such as inputs and results for TDP XAS calculations.
Definition xas_tdp_types.F:195