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