(git:ccc2433)
xas_tdp_kernel.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 All the kernel specific subroutines for XAS TDP calculations
10 !> \author A. Bussy (03.2019)
11 ! **************************************************************************************************
12 
13 MODULE xas_tdp_kernel
14  USE basis_set_types, ONLY: gto_basis_set_p_type
15  USE cp_dbcsr_operations, ONLY: cp_dbcsr_dist2d_to_dist,&
16  dbcsr_deallocate_matrix_set
17  USE dbcsr_api, ONLY: &
18  dbcsr_add, dbcsr_complete_redistribute, dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, &
19  dbcsr_distribution_get, dbcsr_distribution_new, dbcsr_distribution_release, &
20  dbcsr_distribution_type, dbcsr_filter, dbcsr_finalize, dbcsr_get_block_p, dbcsr_get_info, &
21  dbcsr_get_stored_coordinates, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
22  dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, &
23  dbcsr_p_type, dbcsr_put_block, dbcsr_release, dbcsr_reserve_block2d, dbcsr_set, &
24  dbcsr_transposed, dbcsr_type, dbcsr_type_no_symmetry, dbcsr_type_symmetric
25  USE dbt_api, ONLY: dbt_get_block,&
26  dbt_iterator_blocks_left,&
27  dbt_iterator_next_block,&
28  dbt_iterator_start,&
29  dbt_iterator_stop,&
30  dbt_iterator_type,&
31  dbt_type
32  USE distribution_2d_types, ONLY: distribution_2d_type
33  USE kinds, ONLY: dp
34  USE message_passing, ONLY: mp_para_env_type
35  USE particle_methods, ONLY: get_particle_set
36  USE particle_types, ONLY: particle_type
37  USE qs_environment_types, ONLY: get_qs_env,&
38  qs_environment_type
39  USE qs_integral_utils, ONLY: basis_set_list_setup
40  USE qs_kind_types, ONLY: qs_kind_type
41  USE util, ONLY: get_limit
42  USE xas_tdp_types, ONLY: donor_state_type,&
43  get_proc_batch_sizes,&
44  xas_tdp_control_type,&
45  xas_tdp_env_type
46 
47 !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
48 #include "./base/base_uses.f90"
49 
50  IMPLICIT NONE
51  PRIVATE
52 
53  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xas_tdp_kernel'
54 
55  PUBLIC :: kernel_coulomb_xc, kernel_exchange, contract2_ao_to_domo, &
56  reserve_contraction_blocks, ri_all_blocks_mm
57 
58 CONTAINS
59 
60 ! **************************************************************************************************
61 !> \brief Computes, if asked for it, the Coulomb and XC kernel matrices, in the usuall matrix format
62 !> \param coul_ker pointer the the Coulomb kernel matrix (can be void pointer)
63 !> \param xc_ker array of pointer to the different xc kernels (5 of them):
64 !> 1) the restricted closed-shell singlet kernel
65 !> 2) the restricted closed-shell triplet kernel
66 !> 3) the spin-conserving open-shell xc kernel
67 !> 4) the on-diagonal spin-flip open-shell xc kernel
68 !> \param donor_state ...
69 !> \param xas_tdp_env ...
70 !> \param xas_tdp_control ...
71 !> \param qs_env ...
72 !> \note Coulomb and xc kernel are put together in the same routine because they use the same RI
73 !> Coulomb: (aI|Jb) = (aI|P) (P|Q)^-1 (Q|Jb)
74 !> XC : (aI|fxc|Jb) = (aI|P) (P|Q)^-1 (Q|fxc|R) (R|S)^-1 (S|Jb)
75 !> In the above formula, a,b label the sgfs
76 !> The routine analyses the xas_tdp_control to know which kernel must be computed and how
77 !> (open-shell, singlet, triplet, ROKS, LSD, etc...)
78 !> On entry, the pointers should be allocated
79 ! **************************************************************************************************
80  SUBROUTINE kernel_coulomb_xc(coul_ker, xc_ker, donor_state, xas_tdp_env, xas_tdp_control, qs_env)
81 
82  TYPE(dbcsr_type), INTENT(INOUT) :: coul_ker
83  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: xc_ker
84  TYPE(donor_state_type), POINTER :: donor_state
85  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
86  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
87  TYPE(qs_environment_type), POINTER :: qs_env
88 
89  CHARACTER(len=*), PARAMETER :: routinen = 'kernel_coulomb_xc'
90 
91  INTEGER :: batch_size, bo(2), handle, i, ibatch, &
92  iex, lb, natom, nbatch, ndo_mo, &
93  ndo_so, nex_atom, nsgfp, ri_atom, &
94  source, ub
95  INTEGER, DIMENSION(:), POINTER :: blk_size
96  LOGICAL :: do_coulomb, do_sc, do_sf, do_sg, do_tp, &
97  do_xc, found
98  REAL(dp), DIMENSION(:, :), POINTER :: pq
99  TYPE(dbcsr_distribution_type), POINTER :: dist
100  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: contr1_int
101  TYPE(mp_para_env_type), POINTER :: para_env
102 
103  NULLIFY (contr1_int, pq, para_env, dist, blk_size)
104 
105 ! Initialization
106  ndo_mo = donor_state%ndo_mo
107  do_xc = xas_tdp_control%do_xc
108  do_sg = xas_tdp_control%do_singlet
109  do_tp = xas_tdp_control%do_triplet
110  do_sc = xas_tdp_control%do_spin_cons
111  do_sf = xas_tdp_control%do_spin_flip
112  ndo_so = ndo_mo; IF (xas_tdp_control%do_uks) ndo_so = 2*ndo_mo
113  ri_atom = donor_state%at_index
114  CALL get_qs_env(qs_env, natom=natom, para_env=para_env)
115  do_coulomb = xas_tdp_control%do_coulomb
116  dist => donor_state%dbcsr_dist
117  blk_size => donor_state%blk_size
118 
119 ! If no Coulomb nor xc, simply exit
120  IF ((.NOT. do_coulomb) .AND. (.NOT. do_xc)) RETURN
121 
122  CALL timeset(routinen, handle)
123 
124 ! Contract the RI 3-center integrals once to get (aI|P)
125  CALL contract2_ao_to_domo(contr1_int, "COULOMB", donor_state, xas_tdp_env, xas_tdp_control, qs_env)
126 
127 ! Deal with the Coulomb case
128  IF (do_coulomb) CALL coulomb(coul_ker, contr1_int, dist, blk_size, xas_tdp_env, &
129  xas_tdp_control, qs_env)
130 
131 ! Deal with the XC case
132  IF (do_xc) THEN
133 
134  ! In the end, we compute: (aI|fxc|Jb) = (aI|P) (P|Q)^-1 (Q|fxc|R) (R|S)^-1 (S|Jb)
135  ! where fxc can take different spin contributions.
136 
137  ! Precompute the product (aI|P) * (P|Q)^-1 and store it in contr1_int
138  pq => xas_tdp_env%ri_inv_coul
139  CALL ri_all_blocks_mm(contr1_int, pq)
140 
141  ! If not already done (e.g. when multpile donor states for a given excited atom), broadcast
142  ! the RI matrix (Q|fxc|R) on all procs
143  IF (.NOT. xas_tdp_env%fxc_avail) THEN
144  ! Find on which processor the integrals (Q|fxc|R) for this atom are stored
145  nsgfp = SIZE(pq, 1)
146  CALL get_qs_env(qs_env, para_env=para_env)
147  found = .false.
148  nex_atom = SIZE(xas_tdp_env%ex_atom_indices)
149  CALL get_proc_batch_sizes(batch_size, nbatch, nex_atom, para_env%num_pe)
150 
151  DO ibatch = 0, nbatch - 1
152 
153  bo = get_limit(nex_atom, nbatch, ibatch)
154  DO iex = bo(1), bo(2)
155 
156  IF (xas_tdp_env%ex_atom_indices(iex) == ri_atom) THEN
157  source = ibatch*batch_size
158  found = .true. !but simply take the first
159  EXIT
160  END IF
161  END DO !iex
162  IF (found) EXIT
163  END DO !ip
164 
165  ! Broadcast the integrals to all procs (deleted after all donor states for this atoms are treated)
166  lb = 1; IF (do_sf .AND. .NOT. do_sc) lb = 4
167  ub = 2; IF (do_sc) ub = 3; IF (do_sf) ub = 4
168  DO i = lb, ub
169  IF (.NOT. ASSOCIATED(xas_tdp_env%ri_fxc(ri_atom, i)%array)) THEN
170  ALLOCATE (xas_tdp_env%ri_fxc(ri_atom, i)%array(nsgfp, nsgfp))
171  END IF
172  CALL para_env%bcast(xas_tdp_env%ri_fxc(ri_atom, i)%array, source)
173  END DO
174 
175  xas_tdp_env%fxc_avail = .true.
176  END IF
177 
178  ! Case study on the calculation type
179  IF (do_sg .OR. do_tp) THEN
180  CALL rcs_xc(xc_ker(1)%matrix, xc_ker(2)%matrix, contr1_int, dist, blk_size, &
181  donor_state, xas_tdp_env, xas_tdp_control, qs_env)
182  END IF
183 
184  IF (do_sc) THEN
185  CALL sc_os_xc(xc_ker(3)%matrix, contr1_int, dist, blk_size, donor_state, &
186  xas_tdp_env, xas_tdp_control, qs_env)
187  END IF
188 
189  IF (do_sf) THEN
190  CALL ondiag_sf_os_xc(xc_ker(4)%matrix, contr1_int, dist, blk_size, donor_state, &
191  xas_tdp_env, xas_tdp_control, qs_env)
192  END IF
193 
194  END IF ! do_xc
195 
196 ! Clean-up
197  CALL dbcsr_deallocate_matrix_set(contr1_int)
198 
199  CALL timestop(handle)
200 
201  END SUBROUTINE kernel_coulomb_xc
202 
203 ! **************************************************************************************************
204 !> \brief Create the matrix containing the Coulomb kernel, which is:
205 !> (aI_sigma|J_tau b) ~= (aI_sigma|P) * (P|Q) * (Q|J_tau b)
206 !> \param coul_ker the Coulomb kernel
207 !> \param contr1_int the once contracted RI integrals (aI|P)
208 !> \param dist the inherited dbcsr ditribution
209 !> \param blk_size the inherited block sizes
210 !> \param xas_tdp_env ...
211 !> \param xas_tdp_control ...
212 !> \param qs_env ...
213 ! **************************************************************************************************
214  SUBROUTINE coulomb(coul_ker, contr1_int, dist, blk_size, xas_tdp_env, xas_tdp_control, qs_env)
215 
216  TYPE(dbcsr_type), INTENT(INOUT) :: coul_ker
217  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: contr1_int
218  TYPE(dbcsr_distribution_type), POINTER :: dist
219  INTEGER, DIMENSION(:), POINTER :: blk_size
220  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
221  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
222  TYPE(qs_environment_type), POINTER :: qs_env
223 
224  LOGICAL :: quadrants(3)
225  REAL(dp), DIMENSION(:, :), POINTER :: pq
226  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: lhs_int, rhs_int
227  TYPE(dbcsr_type) :: work_mat
228 
229  NULLIFY (pq, rhs_int, lhs_int)
230 
231  ! Get the inver RI coulomb
232  pq => xas_tdp_env%ri_inv_coul
233 
234  ! Create a normal type work matrix
235  CALL dbcsr_create(work_mat, name="WORK", matrix_type=dbcsr_type_no_symmetry, dist=dist, &
236  row_blk_size=blk_size, col_blk_size=blk_size)
237 
238  ! Compute the product (aI|P) * (P|Q)^-1 * (Q|Jb) = (aI|Jb)
239  rhs_int => contr1_int ! the incoming contr1_int is not modified
240  ALLOCATE (lhs_int(SIZE(contr1_int)))
241  CALL copy_ri_contr_int(lhs_int, rhs_int) ! RHS containts (Q|JB)^T
242  CALL ri_all_blocks_mm(lhs_int, pq) ! LHS contatins (aI|P)*(P|Q)^-1
243 
244  !In the special case of ROKS, same MOs for each spin => put same (aI|Jb) product on the
245  !alpha-alpha, alpha-beta and beta-beta quadrants of the kernel matrix.
246  IF (xas_tdp_control%do_roks) THEN
247  quadrants = [.true., .true., .true.]
248  ELSE
249  quadrants = [.true., .false., .false.]
250  END IF
251  CALL ri_int_product(work_mat, lhs_int, rhs_int, quadrants, qs_env, &
252  eps_filter=xas_tdp_control%eps_filter)
253  CALL dbcsr_finalize(work_mat)
254 
255  !Create the symmetric kernel matrix and redistribute work_mat into it
256  CALL dbcsr_create(coul_ker, name="COULOMB KERNEL", matrix_type=dbcsr_type_symmetric, dist=dist, &
257  row_blk_size=blk_size, col_blk_size=blk_size)
258  CALL dbcsr_complete_redistribute(work_mat, coul_ker)
259 
260  !clean-up
261  CALL dbcsr_release(work_mat)
262  CALL dbcsr_deallocate_matrix_set(lhs_int)
263 
264  END SUBROUTINE coulomb
265 
266 ! **************************************************************************************************
267 !> \brief Create the matrix containing the XC kenrel in the spin-conserving open-shell case:
268 !> (aI_sigma|fxc|J_tau b) ~= (aI_sigma|P) (P|Q)^-1 (Q|fxc|R) (R|S)^-1 (S|J_tau b)
269 !> \param xc_ker the kernel matrix
270 !> \param contr1_int_PQ the once contracted RI integrals, with inverse coulomb: (aI_sigma|P) (P|Q)^-1
271 !> \param dist inherited dbcsr dist
272 !> \param blk_size inherited block sizes
273 !> \param donor_state ...
274 !> \param xas_tdp_env ...
275 !> \param xas_tdp_control ...
276 !> \param qs_env ...
277 !> note Prior to calling this function, the (Q|fxc|R) integral must be brodcasted to all procs
278 ! **************************************************************************************************
279  SUBROUTINE sc_os_xc(xc_ker, contr1_int_PQ, dist, blk_size, donor_state, xas_tdp_env, &
280  xas_tdp_control, qs_env)
281 
282  TYPE(dbcsr_type), INTENT(INOUT) :: xc_ker
283  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: contr1_int_pq
284  TYPE(dbcsr_distribution_type), POINTER :: dist
285  INTEGER, DIMENSION(:), POINTER :: blk_size
286  TYPE(donor_state_type), POINTER :: donor_state
287  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
288  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
289  TYPE(qs_environment_type), POINTER :: qs_env
290 
291  INTEGER :: ndo_mo, ri_atom
292  LOGICAL :: quadrants(3)
293  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: lhs_int, rhs_int
294  TYPE(dbcsr_type) :: work_mat
295 
296  NULLIFY (lhs_int, rhs_int)
297 
298  ! Initialization
299  ndo_mo = donor_state%ndo_mo
300  ri_atom = donor_state%at_index
301  !normal type work matrix such that distribution of all spin quadrants match
302  CALL dbcsr_create(work_mat, name="WORK", matrix_type=dbcsr_type_no_symmetry, dist=dist, &
303  row_blk_size=blk_size, col_blk_size=blk_size)
304 
305  rhs_int => contr1_int_pq ! contains [ (aI|P)*(P|Q)^-1 ]^T
306  ALLOCATE (lhs_int(SIZE(contr1_int_pq))) ! will contain (aI|P)*(P|Q)^-1 * (Q|fxc|R)
307 
308  ! Case study: UKS or ROKS ?
309  IF (xas_tdp_control%do_uks) THEN
310 
311  ! In the case of UKS, donor MOs might be different for different spins. Moreover, the
312  ! fxc itself might change since fxc = fxc_sigma,tau
313  ! => Carfully treat each spin-quadrant separately
314 
315  ! alpha-alpha spin quadrant (upper-lefet)
316  quadrants = [.true., .false., .false.]
317 
318  ! Copy the alpha part into lhs_int, multiply by the alpha-alpha (Q|fxc|R) and then
319  ! by the alpha part of rhs_int
320  CALL copy_ri_contr_int(lhs_int(1:ndo_mo), rhs_int(1:ndo_mo))
321  CALL ri_all_blocks_mm(lhs_int(1:ndo_mo), xas_tdp_env%ri_fxc(ri_atom, 1)%array)
322  CALL ri_int_product(work_mat, lhs_int(1:ndo_mo), rhs_int(1:ndo_mo), quadrants, qs_env, &
323  eps_filter=xas_tdp_control%eps_filter)
324 
325  ! alpha-beta spin quadrant (upper-right)
326  quadrants = [.false., .true., .false.]
327 
328  !Copy the alpha part into LHS, multiply by the alpha-beta kernel and the beta part of RHS
329  CALL copy_ri_contr_int(lhs_int(1:ndo_mo), rhs_int(1:ndo_mo))
330  CALL ri_all_blocks_mm(lhs_int(1:ndo_mo), xas_tdp_env%ri_fxc(ri_atom, 2)%array)
331  CALL ri_int_product(work_mat, lhs_int(1:ndo_mo), rhs_int(ndo_mo + 1:2*ndo_mo), &
332  quadrants, qs_env, eps_filter=xas_tdp_control%eps_filter)
333 
334  ! beta-beta spin quadrant (lower-right)
335  quadrants = [.false., .false., .true.]
336 
337  !Copy the beta part into LHS, multiply by the beta-beta kernel and the beta part of RHS
338  CALL copy_ri_contr_int(lhs_int(ndo_mo + 1:2*ndo_mo), rhs_int(ndo_mo + 1:2*ndo_mo))
339  CALL ri_all_blocks_mm(lhs_int(ndo_mo + 1:2*ndo_mo), xas_tdp_env%ri_fxc(ri_atom, 3)%array)
340  CALL ri_int_product(work_mat, lhs_int(ndo_mo + 1:2*ndo_mo), rhs_int(ndo_mo + 1:2*ndo_mo), &
341  quadrants, qs_env, eps_filter=xas_tdp_control%eps_filter)
342 
343  ELSE IF (xas_tdp_control%do_roks) THEN
344 
345  ! In the case of ROKS, fxc = fxc_sigma,tau is different for each spin quadrant, but the
346  ! donor MOs remain the same
347 
348  ! alpha-alpha kernel in the upper left quadrant
349  quadrants = [.true., .false., .false.]
350 
351  !Copy the LHS and multiply by alpha-alpha kernel
352  CALL copy_ri_contr_int(lhs_int, rhs_int)
353  CALL ri_all_blocks_mm(lhs_int, xas_tdp_env%ri_fxc(ri_atom, 1)%array)
354  CALL ri_int_product(work_mat, lhs_int, rhs_int, quadrants, qs_env, &
355  eps_filter=xas_tdp_control%eps_filter)
356 
357  ! alpha-beta kernel in the upper-right quadrant
358  quadrants = [.false., .true., .false.]
359 
360  !Copy LHS and multiply by the alpha-beta kernel
361  CALL copy_ri_contr_int(lhs_int, rhs_int)
362  CALL ri_all_blocks_mm(lhs_int, xas_tdp_env%ri_fxc(ri_atom, 2)%array)
363  CALL ri_int_product(work_mat, lhs_int, rhs_int, quadrants, qs_env, &
364  eps_filter=xas_tdp_control%eps_filter)
365 
366  ! beta-beta kernel in the lower-right quadrant
367  quadrants = [.false., .false., .true.]
368 
369  !Copy the LHS and multiply by the beta-beta kernel
370  CALL copy_ri_contr_int(lhs_int, rhs_int)
371  CALL ri_all_blocks_mm(lhs_int, xas_tdp_env%ri_fxc(ri_atom, 3)%array)
372  CALL ri_int_product(work_mat, lhs_int, rhs_int, quadrants, qs_env, &
373  eps_filter=xas_tdp_control%eps_filter)
374 
375  END IF
376  CALL dbcsr_finalize(work_mat)
377 
378  ! Create a symmetric kernel matrix and redistribute the normal work matrix into it
379  CALL dbcsr_create(xc_ker, name="SC OS XC KERNEL", matrix_type=dbcsr_type_symmetric, dist=dist, &
380  row_blk_size=blk_size, col_blk_size=blk_size)
381  CALL dbcsr_complete_redistribute(work_mat, xc_ker)
382 
383  !clean-up
384  CALL dbcsr_deallocate_matrix_set(lhs_int)
385  CALL dbcsr_release(work_mat)
386 
387  END SUBROUTINE sc_os_xc
388 
389 ! **************************************************************************************************
390 !> \brief Create the matrix containing the on-diagonal spin-flip XC kernel (open-shell), which is:
391 !> (a I_sigma|fxc|J_tau b) * delta_sigma,tau, fxc = 1/(rhoa-rhob) * (dE/drhoa - dE/drhob)
392 !> with RI: (a I_sigma|fxc|J_tau b) ~= (aI_sigma|P) (P|Q)^-1 (Q|fxc|R) (R|S)^-1 (S|J_tau b)
393 !> \param xc_ker the kernel matrix
394 !> \param contr1_int_PQ the once contracted RI integrals with coulomb product: (aI_sigma|P) (P|Q)^-1
395 !> \param dist inherited dbcsr dist
396 !> \param blk_size inherited block sizes
397 !> \param donor_state ...
398 !> \param xas_tdp_env ...
399 !> \param xas_tdp_control ...
400 !> \param qs_env ...
401 !> \note It must be later on multiplied by the spin-swapped Q projector
402 !> Prior to calling this function, the (Q|fxc|R) integral must be brodcasted to all procs
403 ! **************************************************************************************************
404  SUBROUTINE ondiag_sf_os_xc(xc_ker, contr1_int_PQ, dist, blk_size, donor_state, xas_tdp_env, &
405  xas_tdp_control, qs_env)
406 
407  TYPE(dbcsr_type), INTENT(INOUT) :: xc_ker
408  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: contr1_int_pq
409  TYPE(dbcsr_distribution_type), POINTER :: dist
410  INTEGER, DIMENSION(:), POINTER :: blk_size
411  TYPE(donor_state_type), POINTER :: donor_state
412  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
413  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
414  TYPE(qs_environment_type), POINTER :: qs_env
415 
416  INTEGER :: ndo_mo, ri_atom
417  LOGICAL :: quadrants(3)
418  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: lhs_int, rhs_int
419  TYPE(dbcsr_type) :: work_mat
420 
421  NULLIFY (lhs_int, rhs_int)
422 
423  ! Initialization
424  ndo_mo = donor_state%ndo_mo
425  ri_atom = donor_state%at_index
426  !normal type work matrix such that distribution of all spin quadrants match
427  CALL dbcsr_create(work_mat, name="WORK", matrix_type=dbcsr_type_no_symmetry, dist=dist, &
428  row_blk_size=blk_size, col_blk_size=blk_size)
429 
430  !Create a lhs_int, in which the whole (aI_sigma|P) (P|Q)^-1 (Q|fxc|R) will be put
431  !because in spin-flip fxc is spin-independent, can take the product once and for all
432  rhs_int => contr1_int_pq
433  ALLOCATE (lhs_int(SIZE(contr1_int_pq)))
434  CALL copy_ri_contr_int(lhs_int, rhs_int)
435  CALL ri_all_blocks_mm(lhs_int, xas_tdp_env%ri_fxc(ri_atom, 4)%array)
436 
437  ! Case study: UKS or ROKS ?
438  IF (xas_tdp_control%do_uks) THEN
439 
440  ! In the case of UKS, donor MOs might be different for different spins
441  ! => Carfully treat each spin-quadrant separately
442  ! NO alpha-beta because of the delta_sigma,tau
443 
444  ! alpha-alpha spin quadrant (upper-lefet)
445  quadrants = [.true., .false., .false.]
446  CALL ri_int_product(work_mat, lhs_int(1:ndo_mo), rhs_int(1:ndo_mo), quadrants, qs_env, &
447  eps_filter=xas_tdp_control%eps_filter)
448 
449  ! beta-beta spin quadrant (lower-right)
450  quadrants = [.false., .false., .true.]
451  CALL ri_int_product(work_mat, lhs_int(ndo_mo + 1:2*ndo_mo), rhs_int(ndo_mo + 1:2*ndo_mo), &
452  quadrants, qs_env, eps_filter=xas_tdp_control%eps_filter)
453 
454  ELSE IF (xas_tdp_control%do_roks) THEN
455 
456  ! In the case of ROKS, same donor MOs for both spins => can do it all at once
457  ! But NOT the alpha-beta quadrant because of delta_sigma,tau
458 
459  quadrants = [.true., .false., .true.]
460  CALL ri_int_product(work_mat, lhs_int, rhs_int, quadrants, qs_env, &
461  eps_filter=xas_tdp_control%eps_filter)
462 
463  END IF
464  CALL dbcsr_finalize(work_mat)
465 
466  ! Create a symmetric kernel matrix and redistribute the normal work matrix into it
467  CALL dbcsr_create(xc_ker, name="ON-DIAG SF OS XC KERNEL", matrix_type=dbcsr_type_symmetric, &
468  dist=dist, row_blk_size=blk_size, col_blk_size=blk_size)
469  CALL dbcsr_complete_redistribute(work_mat, xc_ker)
470 
471  !clean-up
472  CALL dbcsr_deallocate_matrix_set(lhs_int)
473  CALL dbcsr_release(work_mat)
474 
475  END SUBROUTINE ondiag_sf_os_xc
476 
477 ! **************************************************************************************************
478 !> \brief Create the matrix containing the XC kernel in the restricted closed-shell case, for
479 !> singlets: (aI|fxc|Jb) = (aI|P) (P|Q)^-1 (Q|fxc_alp,alp+fxc_alp,bet|R) (R|S)^-1 (S|Jb)
480 !> triplets: (aI|fxc|Jb) = (aI|P) (P|Q)^-1 (Q|fxc_alp,alp-fxc_alp,bet|R) (R|S)^-1 (S|Jb)
481 !> \param sg_xc_ker the singlet kernel matrix
482 !> \param tp_xc_ker the triplet kernel matrix
483 !> \param contr1_int_PQ the once contracted RI integrals including inverse 2-center Coulomb prodcut:
484 !> (aI|P)*(P|Q)^-1
485 !> \param dist inherited dbcsr dist
486 !> \param blk_size inherited block sizes
487 !> \param donor_state ...
488 !> \param xas_tdp_env ...
489 !> \param xas_tdp_control ...
490 !> \param qs_env ...
491 ! **************************************************************************************************
492  SUBROUTINE rcs_xc(sg_xc_ker, tp_xc_ker, contr1_int_PQ, dist, blk_size, donor_state, &
493  xas_tdp_env, xas_tdp_control, qs_env)
494 
495  TYPE(dbcsr_type), INTENT(INOUT) :: sg_xc_ker, tp_xc_ker
496  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: contr1_int_pq
497  TYPE(dbcsr_distribution_type), POINTER :: dist
498  INTEGER, DIMENSION(:), POINTER :: blk_size
499  TYPE(donor_state_type), POINTER :: donor_state
500  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
501  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
502  TYPE(qs_environment_type), POINTER :: qs_env
503 
504  INTEGER :: nsgfp, ri_atom
505  LOGICAL :: quadrants(3)
506  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: fxc
507  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: lhs_int, rhs_int
508  TYPE(dbcsr_type) :: work_mat
509 
510  NULLIFY (lhs_int, rhs_int)
511 
512  ! Initialization
513  ri_atom = donor_state%at_index
514  nsgfp = SIZE(xas_tdp_env%ri_fxc(ri_atom, 1)%array, 1)
515  rhs_int => contr1_int_pq ! RHS contains [ (aI|P)*(P|Q)^-1 ]^T
516  ALLOCATE (lhs_int(SIZE(contr1_int_pq))) ! LHS will contatin (aI|P)*(P|Q)^-1 * (Q|fxc|R)
517 
518  ! Work structures
519  ALLOCATE (fxc(nsgfp, nsgfp))
520  CALL dbcsr_create(work_mat, name="WORK", matrix_type=dbcsr_type_no_symmetry, dist=dist, &
521  row_blk_size=blk_size, col_blk_size=blk_size)
522 
523  ! Case study: singlet and/or triplet ?
524  IF (xas_tdp_control%do_singlet) THEN
525 
526  ! Take the sum of fxc for alpha-alpha and alpha-beta
527  CALL dcopy(nsgfp*nsgfp, xas_tdp_env%ri_fxc(ri_atom, 1)%array, 1, fxc, 1)
528  CALL daxpy(nsgfp*nsgfp, 1.0_dp, xas_tdp_env%ri_fxc(ri_atom, 2)%array, 1, fxc, 1)
529 
530  ! Copy the fresh lhs_int = (aI|P) (P|Q)^-1 and multiply by (Q|fxc|R)
531  CALL copy_ri_contr_int(lhs_int, rhs_int)
532  CALL ri_all_blocks_mm(lhs_int, fxc)
533 
534  ! Compute the final LHS RHS product => spin-restricted, only upper-left quadrant
535  quadrants = [.true., .false., .false.]
536  CALL ri_int_product(work_mat, lhs_int, rhs_int, quadrants, qs_env, &
537  eps_filter=xas_tdp_control%eps_filter)
538  CALL dbcsr_finalize(work_mat)
539 
540  !Create the symmetric kernel matrix and redistribute work_mat into it
541  CALL dbcsr_create(sg_xc_ker, name="XC SINGLET KERNEL", matrix_type=dbcsr_type_symmetric, &
542  dist=dist, row_blk_size=blk_size, col_blk_size=blk_size)
543  CALL dbcsr_complete_redistribute(work_mat, sg_xc_ker)
544 
545  END IF
546 
547  IF (xas_tdp_control%do_triplet) THEN
548 
549  ! Take the difference of fxc for alpha-alpha and alpha-beta
550  CALL dcopy(nsgfp*nsgfp, xas_tdp_env%ri_fxc(ri_atom, 1)%array, 1, fxc, 1)
551  CALL daxpy(nsgfp*nsgfp, -1.0_dp, xas_tdp_env%ri_fxc(ri_atom, 2)%array, 1, fxc, 1)
552 
553  ! Copy the fresh lhs_int = (aI|P) (P|Q)^-1 and multiply by (Q|fxc|R)
554  CALL copy_ri_contr_int(lhs_int, rhs_int)
555  CALL ri_all_blocks_mm(lhs_int, fxc)
556 
557  ! Compute the final LHS RHS product => spin-restricted, only upper-left quadrant
558  quadrants = [.true., .false., .false.]
559  CALL ri_int_product(work_mat, lhs_int, rhs_int, quadrants, qs_env, &
560  eps_filter=xas_tdp_control%eps_filter)
561  CALL dbcsr_finalize(work_mat)
562 
563  !Create the symmetric kernel matrix and redistribute work_mat into it
564  CALL dbcsr_create(tp_xc_ker, name="XC TRIPLET KERNEL", matrix_type=dbcsr_type_symmetric, &
565  dist=dist, row_blk_size=blk_size, col_blk_size=blk_size)
566  CALL dbcsr_complete_redistribute(work_mat, tp_xc_ker)
567 
568  END IF
569 
570  ! clean-up
571  CALL dbcsr_deallocate_matrix_set(lhs_int)
572  CALL dbcsr_release(work_mat)
573  DEALLOCATE (fxc)
574 
575  END SUBROUTINE rcs_xc
576 
577 ! **************************************************************************************************
578 !> \brief Computes the exact exchange kernel matrix using RI. Returns an array of 2 matrices,
579 !> which are:
580 !> 1) the on-diagonal kernel: (ab|I_sigma J_tau) * delta_sigma,tau
581 !> 2) the off-diagonal spin-conserving kernel: (aJ_sigma|I_tau b) * delta_sigma,tau
582 !> An internal analysis determines which of the above are computed (can range from 0 to 2),
583 !> \param ex_ker ...
584 !> \param donor_state ...
585 !> \param xas_tdp_env ...
586 !> \param xas_tdp_control ...
587 !> \param qs_env ...
588 !> \note In the case of spin-conserving excitation, the kernel must later be multiplied by the
589 !> usual Q projector. In the case of spin-flip, one needs to project the excitations coming
590 !> from alpha donor MOs on the unoccupied beta MOs. This is done by multiplying by a Q
591 !> projector where the alpha-alpha and beta-beta quadrants are swapped
592 !> The ex_ker array should be allocated on entry (not the internals)
593 ! **************************************************************************************************
594  SUBROUTINE kernel_exchange(ex_ker, donor_state, xas_tdp_env, xas_tdp_control, qs_env)
595 
596  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ex_ker
597  TYPE(donor_state_type), POINTER :: donor_state
598  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
599  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
600  TYPE(qs_environment_type), POINTER :: qs_env
601 
602  CHARACTER(len=*), PARAMETER :: routinen = 'kernel_exchange'
603 
604  INTEGER :: handle
605  INTEGER, DIMENSION(:), POINTER :: blk_size
606  LOGICAL :: do_off_sc
607  TYPE(dbcsr_distribution_type), POINTER :: dist
608  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: contr1_int
609 
610  NULLIFY (contr1_int, dist, blk_size)
611 
612  !Don't do anything if no hfx
613  IF (.NOT. xas_tdp_control%do_hfx) RETURN
614 
615  CALL timeset(routinen, handle)
616 
617  dist => donor_state%dbcsr_dist
618  blk_size => donor_state%blk_size
619 
620  !compute the off-diag spin-conserving only if not TDA and anything that is spin-conserving
621  do_off_sc = (.NOT. xas_tdp_control%tamm_dancoff) .AND. &
622  (xas_tdp_control%do_spin_cons .OR. xas_tdp_control%do_singlet .OR. xas_tdp_control%do_triplet)
623 
624  ! Need the once contracted integrals (aI|P)
625  CALL contract2_ao_to_domo(contr1_int, "EXCHANGE", donor_state, xas_tdp_env, xas_tdp_control, qs_env)
626 
627 ! The on-diagonal exchange : (ab|P) * (P|Q)^-1 * (Q|I_sigma J_tau) * delta_sigma,tau
628  CALL ondiag_ex(ex_ker(1)%matrix, contr1_int, dist, blk_size, donor_state, xas_tdp_env, &
629  xas_tdp_control, qs_env)
630 
631 ! The off-diag spin-conserving case: (aJ_sigma|P) * (P|Q)^-1 * (Q|I_tau b) * delta_sigma,tau
632  IF (do_off_sc) THEN
633  CALL offdiag_ex_sc(ex_ker(2)%matrix, contr1_int, dist, blk_size, donor_state, &
634  xas_tdp_env, xas_tdp_control, qs_env)
635  END IF
636 
637  !clean-up
638  CALL dbcsr_deallocate_matrix_set(contr1_int)
639 
640  CALL timestop(handle)
641 
642  END SUBROUTINE kernel_exchange
643 
644 ! **************************************************************************************************
645 !> \brief Create the matrix containing the on-diagonal exact exchange kernel, which is:
646 !> (ab|I_sigma J_tau) * delta_sigma,tau, where a,b are AOs, I_sigma and J_tau are the donor
647 !> spin-orbitals. A RI is done: (ab|I_sigma J_tau) = (ab|P) * (P|Q)^-1 * (Q|I_sigma J_tau)
648 !> \param ondiag_ex_ker the on-diagonal exchange kernel in dbcsr format
649 !> \param contr1_int the already once-contracted RI 3-center integrals (aI_sigma|P)
650 !> where each matrix of the array contains the contraction for the donor spin-orbital I_sigma
651 !> \param dist the inherited dbcsr distribution
652 !> \param blk_size the inherited dbcsr block sizes
653 !> \param donor_state ...
654 !> \param xas_tdp_env ...
655 !> \param xas_tdp_control ...
656 !> \param qs_env ...
657 !> \note In the presence of a RI metric, we have instead M^-1 * (P|Q) * M^-1
658 ! **************************************************************************************************
659  SUBROUTINE ondiag_ex(ondiag_ex_ker, contr1_int, dist, blk_size, donor_state, xas_tdp_env, &
660  xas_tdp_control, qs_env)
661 
662  TYPE(dbcsr_type), INTENT(INOUT) :: ondiag_ex_ker
663  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: contr1_int
664  TYPE(dbcsr_distribution_type), POINTER :: dist
665  INTEGER, DIMENSION(:), POINTER :: blk_size
666  TYPE(donor_state_type), POINTER :: donor_state
667  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
668  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
669  TYPE(qs_environment_type), POINTER :: qs_env
670 
671  INTEGER :: blk, group, iblk, iso, jblk, jso, nblk, &
672  ndo_mo, ndo_so, nsgfa, nsgfp, ri_atom, &
673  source
674  INTEGER, DIMENSION(:), POINTER :: col_dist, col_dist_work, row_dist, &
675  row_dist_work
676  INTEGER, DIMENSION(:, :), POINTER :: pgrid
677  LOGICAL :: do_roks, do_uks, found
678  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: coeffs, ri_coeffs
679  REAL(dp), DIMENSION(:, :), POINTER :: aiq, pblock, pq
680  TYPE(dbcsr_distribution_type) :: opt_dbcsr_dist, work_dbcsr_dist
681  TYPE(dbcsr_iterator_type) :: iter
682  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
683  TYPE(dbcsr_type) :: abij, mats_desymm, work_mat
684  TYPE(mp_para_env_type), POINTER :: para_env
685 
686  NULLIFY (para_env, matrix_s, pblock, aiq, row_dist, col_dist, row_dist_work, col_dist_work, pgrid)
687 
688  ! We want to compute (ab|I_sigma J_tau) = (ab|P) * (P|Q)^-1 * (Q|I_sigma J_tau)
689  ! Already have (cJ_tau|P) stored in contr1_int. Need to further contract the
690  ! AOs with the coeff of the I_alpha spin-orbital.
691 
692  ! Initialization
693  ndo_mo = donor_state%ndo_mo
694  ri_atom = donor_state%at_index
695  do_roks = xas_tdp_control%do_roks
696  do_uks = xas_tdp_control%do_uks
697  ndo_so = ndo_mo; IF (do_uks) ndo_so = 2*ndo_mo !if not UKS, same donor MOs for both spins
698  pq => xas_tdp_env%ri_inv_ex
699 
700  CALL get_qs_env(qs_env, para_env=para_env, matrix_s=matrix_s, natom=nblk)
701  nsgfp = SIZE(pq, 1)
702  nsgfa = SIZE(donor_state%contract_coeffs, 1)
703  ALLOCATE (coeffs(nsgfp, ndo_so), ri_coeffs(nsgfp, ndo_so))
704 
705  ! a and b need to overlap for non-zero (ab|IJ) => same block structure as overlap S
706  ! need compatible distribution_2d with 3c tensor + normal type
707  CALL cp_dbcsr_dist2d_to_dist(xas_tdp_env%opt_dist2d_ex, opt_dbcsr_dist)
708 
709  CALL dbcsr_desymmetrize(matrix_s(1)%matrix, mats_desymm)
710 
711  CALL dbcsr_create(abij, template=mats_desymm, name="(ab|IJ)", dist=opt_dbcsr_dist)
712  CALL dbcsr_complete_redistribute(mats_desymm, abij)
713 
714  CALL dbcsr_release(mats_desymm)
715 
716  ! Create a work distribution based on opt_dbcsr_dist, but for full size matrices
717  CALL dbcsr_distribution_get(opt_dbcsr_dist, row_dist=row_dist, col_dist=col_dist, group=group, &
718  pgrid=pgrid)
719 
720  ALLOCATE (row_dist_work(ndo_so*nblk))
721  ALLOCATE (col_dist_work(ndo_so*nblk))
722  DO iso = 1, ndo_so
723  row_dist_work((iso - 1)*nblk + 1:iso*nblk) = row_dist(:)
724  col_dist_work((iso - 1)*nblk + 1:iso*nblk) = col_dist(:)
725  END DO
726 
727  CALL dbcsr_distribution_new(work_dbcsr_dist, group=group, pgrid=pgrid, row_dist=row_dist_work, &
728  col_dist=col_dist_work)
729 
730  CALL dbcsr_create(work_mat, name="WORK", matrix_type=dbcsr_type_no_symmetry, dist=work_dbcsr_dist, &
731  row_blk_size=blk_size, col_blk_size=blk_size)
732 
733  ! Loop over donor spin-orbitals. End matrix is symmetric => span only upper half
734  DO iso = 1, ndo_so
735 
736  ! take the (aI|Q) block in contr1_int that has a centered on the excited atom
737  CALL dbcsr_get_stored_coordinates(contr1_int(iso)%matrix, ri_atom, ri_atom, source)
738  IF (para_env%mepos == source) THEN
739  CALL dbcsr_get_block_p(contr1_int(iso)%matrix, ri_atom, ri_atom, aiq, found)
740  ELSE
741  ALLOCATE (aiq(nsgfa, nsgfp))
742  END IF
743  CALL para_env%bcast(aiq, source)
744 
745  ! get the contraction (Q|IJ) by taking (Q|Ia)*contract_coeffs and put it in coeffs
746  CALL dgemm('T', 'N', nsgfp, ndo_so, nsgfa, 1.0_dp, aiq, nsgfa, donor_state%contract_coeffs, &
747  nsgfa, 0.0_dp, coeffs, nsgfp)
748 
749  ! take (P|Q)^-1 * (Q|IJ) and put that in ri_coeffs
750  CALL dgemm('N', 'N', nsgfp, ndo_so, nsgfp, 1.0_dp, pq, nsgfp, coeffs, nsgfp, 0.0_dp, &
751  ri_coeffs, nsgfp)
752 
753  IF (.NOT. para_env%mepos == source) DEALLOCATE (aiq)
754 
755  DO jso = iso, ndo_so
756 
757  ! There is no alpha-beta exchange. In case of UKS, iso,jso span all spin-orbitals
758  ! => CYCLE if iso and jso are indexing MOs with different spin (and we have UKS)
759  IF (do_uks .AND. (iso <= ndo_mo .AND. jso > ndo_mo)) cycle
760 
761  ! compute (ab|IJ) = sum_P (ab|P) * (P|Q)^-1 * (Q|IJ)
762  CALL dbcsr_set(abij, 0.0_dp)
763  CALL contract3_ri_to_domos(xas_tdp_env%ri_3c_ex, ri_coeffs(:, jso), abij, ri_atom)
764 
765  ! Loop over (ab|IJ) and copy into work. OK because dist are made to match
766  CALL dbcsr_iterator_start(iter, abij)
767  DO WHILE (dbcsr_iterator_blocks_left(iter))
768 
769  CALL dbcsr_iterator_next_block(iter, row=iblk, column=jblk, blk=blk)
770  IF (iso == jso .AND. jblk < iblk) cycle
771 
772  CALL dbcsr_get_block_p(abij, iblk, jblk, pblock, found)
773 
774  IF (found) THEN
775  CALL dbcsr_put_block(work_mat, (iso - 1)*nblk + iblk, (jso - 1)*nblk + jblk, pblock)
776 
777  !In case of ROKS, we have (ab|IJ) for alpha-alpha spin, but it is the same for
778  !beta-beta => replicate the blocks (alpha-beta is zero)
779  IF (do_roks) THEN
780  !the beta-beta block
781  CALL dbcsr_put_block(work_mat, (ndo_so + iso - 1)*nblk + iblk, &
782  (ndo_so + jso - 1)*nblk + jblk, pblock)
783  END IF
784  END IF
785 
786  END DO !iterator
787  CALL dbcsr_iterator_stop(iter)
788 
789  END DO !jso
790  END DO !iso
791 
792  CALL dbcsr_finalize(work_mat)
793  CALL dbcsr_create(ondiag_ex_ker, name="ONDIAG EX KERNEL", matrix_type=dbcsr_type_symmetric, &
794  dist=dist, row_blk_size=blk_size, col_blk_size=blk_size)
795  CALL dbcsr_complete_redistribute(work_mat, ondiag_ex_ker)
796 
797  !Clean-up
798  CALL dbcsr_release(work_mat)
799  CALL dbcsr_release(abij)
800  CALL dbcsr_distribution_release(opt_dbcsr_dist)
801  CALL dbcsr_distribution_release(work_dbcsr_dist)
802  DEALLOCATE (col_dist_work, row_dist_work)
803 
804  END SUBROUTINE ondiag_ex
805 
806 ! **************************************************************************************************
807 !> \brief Create the matrix containing the off-diagonal exact exchange kernel in the spin-conserving
808 !> case (which also includes excitations from the closed=shell ref state ) This matrix reads:
809 !> (aJ_sigma|I_tau b) * delta_sigma,tau , where a, b are AOs and J_sigma, I_tau are the donor
810 !> spin-orbital. A RI is done: (aJ_sigma|I_tau b) = (aJ_sigma|P) * (P|Q)^-1 * (Q|I_tau b)
811 !> \param offdiag_ex_ker the off-diagonal, spin-conserving exchange kernel in dbcsr format
812 !> \param contr1_int the once-contracted RI integrals: (aJ_sigma|P)
813 !> \param dist the inherited dbcsr ditribution
814 !> \param blk_size the inherited block sizes
815 !> \param donor_state ...
816 !> \param xas_tdp_env ...
817 !> \param xas_tdp_control ...
818 !> \param qs_env ...
819 ! **************************************************************************************************
820  SUBROUTINE offdiag_ex_sc(offdiag_ex_ker, contr1_int, dist, blk_size, donor_state, xas_tdp_env, &
821  xas_tdp_control, qs_env)
822 
823  TYPE(dbcsr_type), INTENT(INOUT) :: offdiag_ex_ker
824  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: contr1_int
825  TYPE(dbcsr_distribution_type), POINTER :: dist
826  INTEGER, DIMENSION(:), POINTER :: blk_size
827  TYPE(donor_state_type), POINTER :: donor_state
828  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
829  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
830  TYPE(qs_environment_type), POINTER :: qs_env
831 
832  INTEGER :: ndo_mo
833  LOGICAL :: do_roks, do_uks, quadrants(3)
834  REAL(dp), DIMENSION(:, :), POINTER :: pq
835  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: lhs_int, rhs_int
836  TYPE(dbcsr_type) :: work_mat
837 
838  NULLIFY (pq, lhs_int, rhs_int)
839 
840  !Initialization
841  ndo_mo = donor_state%ndo_mo
842  do_roks = xas_tdp_control%do_roks
843  do_uks = xas_tdp_control%do_uks
844  pq => xas_tdp_env%ri_inv_ex
845 
846  rhs_int => contr1_int
847  ALLOCATE (lhs_int(SIZE(contr1_int)))
848  CALL copy_ri_contr_int(lhs_int, rhs_int)
849  CALL ri_all_blocks_mm(lhs_int, pq)
850 
851  !Given the lhs_int and rhs_int, all we need to do is multiply elements from the former by
852  !the transpose of the later, and put the result in the correct spin quadrants
853 
854  !Create a normal type work matrix
855  CALL dbcsr_create(work_mat, name="WORK", matrix_type=dbcsr_type_no_symmetry, dist=dist, &
856  row_blk_size=blk_size, col_blk_size=blk_size)
857 
858  !Case study on closed-shell, ROKS or UKS
859  IF (do_roks) THEN
860  !In ROKS, the donor MOs for each spin are the same => copy the product in both the
861  !alpha-alpha and the beta-beta quadrants
862  quadrants = [.true., .false., .true.]
863  CALL ri_int_product(work_mat, lhs_int, rhs_int, quadrants, qs_env, &
864  eps_filter=xas_tdp_control%eps_filter, mo_transpose=.true.)
865 
866  ELSE IF (do_uks) THEN
867  !In UKS, the donor MOs are possibly different for each spin => start with the
868  !alpha-alpha product and the perform the beta-beta product separately
869  quadrants = [.true., .false., .false.]
870  CALL ri_int_product(work_mat, lhs_int(1:ndo_mo), rhs_int(1:ndo_mo), quadrants, &
871  qs_env, eps_filter=xas_tdp_control%eps_filter, mo_transpose=.true.)
872 
873  quadrants = [.false., .false., .true.]
874  CALL ri_int_product(work_mat, lhs_int(ndo_mo + 1:2*ndo_mo), rhs_int(ndo_mo + 1:2*ndo_mo), &
875  quadrants, qs_env, eps_filter=xas_tdp_control%eps_filter, mo_transpose=.true.)
876  ELSE
877  !In the restricted closed-shell case, only have one spin and a single qudarant
878  quadrants = [.true., .false., .false.]
879  CALL ri_int_product(work_mat, lhs_int, rhs_int, quadrants, qs_env, &
880  eps_filter=xas_tdp_control%eps_filter, mo_transpose=.true.)
881  END IF
882  CALL dbcsr_finalize(work_mat)
883 
884  !Create the symmetric kernel matrix and redistribute work_mat into it
885  CALL dbcsr_create(offdiag_ex_ker, name="OFFDIAG EX KERNEL", matrix_type=dbcsr_type_symmetric, &
886  dist=dist, row_blk_size=blk_size, col_blk_size=blk_size)
887  CALL dbcsr_complete_redistribute(work_mat, offdiag_ex_ker)
888 
889  !clean-up
890  CALL dbcsr_release(work_mat)
891  CALL dbcsr_deallocate_matrix_set(lhs_int)
892 
893  END SUBROUTINE offdiag_ex_sc
894 
895 ! **************************************************************************************************
896 !> \brief Reserves the blocks in of a dbcsr matrix as needed for RI 3-center contraction (aI|P)
897 !> \param matrices the matrices for which blocks are reserved
898 !> \param ri_atom the index of the atom on which RI is done (= all coeffs of I are there, and P too)
899 !> \param qs_env ...
900 !> \note the end product are normal type matrices that are possibly slightly spraser as matrix_s
901 ! **************************************************************************************************
902  SUBROUTINE reserve_contraction_blocks(matrices, ri_atom, qs_env)
903 
904  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrices
905  INTEGER, INTENT(IN) :: ri_atom
906  TYPE(qs_environment_type), POINTER :: qs_env
907 
908  INTEGER :: blk, i, iblk, jblk, n
909  LOGICAL :: found
910  REAL(dp), DIMENSION(:, :), POINTER :: pblock_m, pblock_s
911  TYPE(dbcsr_distribution_type) :: dist
912  TYPE(dbcsr_iterator_type) :: iter
913  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
914  TYPE(dbcsr_type) :: template, work
915 
916  NULLIFY (matrix_s, pblock_s, pblock_m)
917 
918 ! Initialization
919  CALL get_qs_env(qs_env, matrix_s=matrix_s)
920  n = SIZE(matrices)
921  CALL dbcsr_get_info(matrices(1)%matrix, distribution=dist)
922 
923 ! Need to redistribute matrix_s in the distribution of matrices
924  CALL dbcsr_create(work, template=matrix_s(1)%matrix, dist=dist)
925  CALL dbcsr_complete_redistribute(matrix_s(1)%matrix, work)
926 
927 ! Need to desymmetrize matrix as as a template
928  CALL dbcsr_desymmetrize(work, template)
929 
930 ! Loop over matrix_s as need a,b to overlap
931  CALL dbcsr_iterator_start(iter, template)
932  DO WHILE (dbcsr_iterator_blocks_left(iter))
933 
934  CALL dbcsr_iterator_next_block(iter, row=iblk, column=jblk, blk=blk)
935  !only have a,b pair if one of them is the ri_atom
936  IF (iblk .NE. ri_atom .AND. jblk .NE. ri_atom) cycle
937 
938  CALL dbcsr_get_block_p(template, iblk, jblk, pblock_s, found)
939 
940  IF (found) THEN
941  DO i = 1, n
942  NULLIFY (pblock_m)
943  CALL dbcsr_reserve_block2d(matrices(i)%matrix, iblk, jblk, pblock_m)
944  pblock_m = 0.0_dp
945  END DO
946  END IF
947 
948  END DO !dbcsr iter
949  CALL dbcsr_iterator_stop(iter)
950  DO i = 1, n
951  CALL dbcsr_finalize(matrices(i)%matrix)
952  END DO
953 
954 ! Clean-up
955  CALL dbcsr_release(template)
956  CALL dbcsr_release(work)
957 
958  END SUBROUTINE reserve_contraction_blocks
959 
960 ! **************************************************************************************************
961 !> \brief Contract the ri 3-center integrals stored in a tensor with repect to the donor MOs coeffs,
962 !> for a given excited atom k => (aI|k) = sum_b c_Ib (ab|k)
963 !> \param contr_int the contracted integrals as array of dbcsr matrices
964 !> \param op_type for which operator type we contract (COULOMB or EXCHANGE)
965 !> \param donor_state ...
966 !> \param xas_tdp_env ...
967 !> \param xas_tdp_control ...
968 !> \param qs_env ...
969 !> \note In the output matrices, (aI_b|k) is stored at block a,b where I_b is the partial
970 !> contraction that only includes coeffs from atom b. Note that the contracted matrix is
971 !> not symmetric. To get the fully contracted matrix over b, one need to add the block
972 !> columns of (aI_b|k) (get an array of size nao*nsgfp). This step is unnessary in our case
973 !> because we assume locality of donor state, and only one column od (aI_b|k) is pouplated
974 ! **************************************************************************************************
975  SUBROUTINE contract2_ao_to_domo(contr_int, op_type, donor_state, xas_tdp_env, xas_tdp_control, qs_env)
976 
977  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: contr_int
978  CHARACTER(len=*), INTENT(IN) :: op_type
979  TYPE(donor_state_type), POINTER :: donor_state
980  TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
981  TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
982  TYPE(qs_environment_type), POINTER :: qs_env
983 
984  CHARACTER(len=*), PARAMETER :: routinen = 'contract2_AO_to_doMO'
985 
986  INTEGER :: handle, i, imo, ispin, katom, kkind, &
987  natom, ndo_mo, ndo_so, nkind, nspins
988  INTEGER, DIMENSION(:), POINTER :: ri_blk_size, std_blk_size
989  LOGICAL :: do_uks
990  REAL(dp), DIMENSION(:, :), POINTER :: coeffs
991  TYPE(dbcsr_distribution_type) :: opt_dbcsr_dist
992  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrices, matrix_s
993  TYPE(dbcsr_type), POINTER :: ai_p, p_ib, work
994  TYPE(dbt_type), POINTER :: pq_x
995  TYPE(distribution_2d_type), POINTER :: opt_dist2d
996  TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: ri_basis
997  TYPE(mp_para_env_type), POINTER :: para_env
998  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
999  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1000 
1001  NULLIFY (matrix_s, std_blk_size, ri_blk_size, qs_kind_set, ri_basis, pq_x)
1002  NULLIFY (ai_p, p_ib, work, matrices, coeffs, opt_dist2d, particle_set)
1003 
1004  CALL timeset(routinen, handle)
1005 
1006 ! Initialization
1007  CALL get_qs_env(qs_env, natom=natom, matrix_s=matrix_s, qs_kind_set=qs_kind_set, para_env=para_env)
1008  ndo_mo = donor_state%ndo_mo
1009  kkind = donor_state%kind_index
1010  katom = donor_state%at_index
1011  !by default contract for Coulomb
1012  pq_x => xas_tdp_env%ri_3c_coul
1013  opt_dist2d => xas_tdp_env%opt_dist2d_coul
1014  IF (op_type == "EXCHANGE") THEN
1015  cpassert(ASSOCIATED(xas_tdp_env%ri_3c_ex))
1016  pq_x => xas_tdp_env%ri_3c_ex
1017  opt_dist2d => xas_tdp_env%opt_dist2d_ex
1018  END IF
1019  do_uks = xas_tdp_control%do_uks
1020  nspins = 1; IF (do_uks) nspins = 2
1021  ndo_so = nspins*ndo_mo
1022 
1023 ! contracted integrals block sizes
1024  CALL dbcsr_get_info(matrix_s(1)%matrix, col_blk_size=std_blk_size)
1025  ! getting the block dimensions for the RI basis
1026  CALL get_qs_env(qs_env, particle_set=particle_set, nkind=nkind)
1027  ALLOCATE (ri_basis(nkind), ri_blk_size(natom))
1028  CALL basis_set_list_setup(ri_basis, "RI_XAS", qs_kind_set)
1029  CALL get_particle_set(particle_set, qs_kind_set, nsgf=ri_blk_size, basis=ri_basis)
1030 
1031 ! Create work matrices. Everything that goes into a 3c routine must be compatible with the optimal dist_2d
1032  CALL cp_dbcsr_dist2d_to_dist(opt_dist2d, opt_dbcsr_dist)
1033 
1034  ALLOCATE (ai_p, p_ib, work, matrices(2))
1035  CALL dbcsr_create(ai_p, dist=opt_dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, name="(aI|P)", &
1036  row_blk_size=std_blk_size, col_blk_size=ri_blk_size)
1037 
1038  CALL dbcsr_create(p_ib, dist=opt_dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, name="(P|Ib)", &
1039  row_blk_size=ri_blk_size, col_blk_size=std_blk_size)
1040 
1041  !reserve the blocks (needed for 3c contraction routines)
1042  matrices(1)%matrix => ai_p; matrices(2)%matrix => p_ib
1043  CALL reserve_contraction_blocks(matrices, katom, qs_env)
1044  DEALLOCATE (matrices)
1045 
1046  ! Create the contracted integral matrices
1047  ALLOCATE (contr_int(ndo_so))
1048  DO i = 1, ndo_so
1049  ALLOCATE (contr_int(i)%matrix)
1050  CALL dbcsr_create(matrix=contr_int(i)%matrix, template=matrix_s(1)%matrix, &
1051  matrix_type=dbcsr_type_no_symmetry, row_blk_size=std_blk_size, &
1052  col_blk_size=ri_blk_size)
1053  END DO
1054 
1055  ! Only take the coeffs for atom on which MOs I,J are localized
1056  coeffs => donor_state%contract_coeffs
1057 
1058  DO ispin = 1, nspins
1059 
1060 ! Loop over the donor MOs and contract
1061  DO imo = 1, ndo_mo
1062 
1063  ! do the contraction
1064  CALL dbcsr_set(ai_p, 0.0_dp); CALL dbcsr_set(p_ib, 0.0_dp)
1065  CALL contract2_ao_to_domo_low(pq_x, coeffs(:, (ispin - 1)*ndo_mo + imo), ai_p, p_ib, katom)
1066 
1067  ! Get the full (aI|P) contracted integrals
1068  CALL dbcsr_transposed(work, p_ib)
1069  CALL dbcsr_add(work, ai_p, 1.0_dp, 1.0_dp)
1070  CALL dbcsr_complete_redistribute(work, contr_int((ispin - 1)*ndo_mo + imo)%matrix)
1071  CALL dbcsr_filter(contr_int((ispin - 1)*ndo_mo + imo)%matrix, 1.0e-16_dp)
1072 
1073  CALL dbcsr_release(work)
1074  END DO !imo
1075  END DO !ispin
1076 
1077 ! Clean-up
1078  CALL dbcsr_release(ai_p)
1079  CALL dbcsr_release(p_ib)
1080  CALL dbcsr_distribution_release(opt_dbcsr_dist)
1081  DEALLOCATE (ri_blk_size, ai_p, p_ib, work, ri_basis)
1082 
1083  CALL timestop(handle)
1084 
1085  END SUBROUTINE contract2_ao_to_domo
1086 
1087 ! **************************************************************************************************
1088 !> \brief Contraction of the 3-center integrals (ab|Q) over the RI basis elements Q to get donor MOS
1089 !> => (ab|IJ) = sum_X (ab|Q) coeffs_Q
1090 !> \param ab_Q the tensor holding the integrals
1091 !> \param vec the contraction coefficients
1092 !> \param mat_abIJ the matrix holding the (ab|IJ) integrals (blocks must be reserved)
1093 !> \param atom_k the atom for which we contract, i.e. we only take RI basis Q centered on atom_k
1094 !> \note By construction, distribution of tensor and matrix match, also for OMP threads
1095 ! **************************************************************************************************
1096  SUBROUTINE contract3_ri_to_domos(ab_Q, vec, mat_abIJ, atom_k)
1097 
1098  TYPE(dbt_type) :: ab_q
1099  REAL(dp), DIMENSION(:), INTENT(IN) :: vec
1100  TYPE(dbcsr_type) :: mat_abij
1101  INTEGER, INTENT(IN) :: atom_k
1102 
1103  CHARACTER(len=*), PARAMETER :: routinen = 'contract3_RI_to_doMOs'
1104 
1105  INTEGER :: handle, i, iatom, ind(3), j, jatom, katom
1106  LOGICAL :: found, t_found
1107  REAL(dp) :: prefac
1108  REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: iabc
1109  REAL(kind=dp), DIMENSION(:, :), POINTER :: pblock
1110  TYPE(dbcsr_type) :: work
1111  TYPE(dbt_iterator_type) :: iter
1112 
1113  NULLIFY (pblock)
1114 
1115  CALL timeset(routinen, handle)
1116 
1117 !$OMP PARALLEL DEFAULT(NONE) &
1118 !$OMP SHARED(ab_Q,vec,mat_abIJ,atom_k) &
1119 !$OMP PRIVATE(iter,ind,iatom,jatom,katom,prefac,iabc,t_found,found,pblock,i,j)
1120  CALL dbt_iterator_start(iter, ab_q)
1121  DO WHILE (dbt_iterator_blocks_left(iter))
1122  CALL dbt_iterator_next_block(iter, ind)
1123 
1124  iatom = ind(1)
1125  jatom = ind(2)
1126  katom = ind(3)
1127 
1128  IF (.NOT. atom_k == katom) cycle
1129 
1130  prefac = 1.0_dp
1131  IF (iatom == jatom) prefac = 0.5_dp
1132 
1133  CALL dbt_get_block(ab_q, ind, iabc, t_found)
1134 
1135  CALL dbcsr_get_block_p(mat_abij, iatom, jatom, pblock, found)
1136  IF ((.NOT. found) .OR. (.NOT. t_found)) cycle
1137 
1138  DO i = 1, SIZE(pblock, 1)
1139  DO j = 1, SIZE(pblock, 2)
1140 !$OMP ATOMIC
1141  pblock(i, j) = pblock(i, j) + prefac*dot_product(vec(:), iabc(i, j, :))
1142  END DO
1143  END DO
1144 
1145  DEALLOCATE (iabc)
1146  END DO !iter
1147  CALL dbt_iterator_stop(iter)
1148 !$OMP END PARALLEL
1149 
1150  !matrix only half filled => need to add its transpose
1151  CALL dbcsr_create(work, template=mat_abij)
1152  CALL dbcsr_transposed(work, mat_abij)
1153  CALL dbcsr_add(mat_abij, work, 1.0_dp, 1.0_dp)
1154  CALL dbcsr_release(work)
1155 
1156  CALL timestop(handle)
1157 
1158  END SUBROUTINE contract3_ri_to_domos
1159 
1160 ! **************************************************************************************************
1161 !> \brief Contraction of the 3-center integrals over index 1 and 2, for a given atom_k. The results
1162 !> are stored in two matrices, such that (a,b are block indices):
1163 !> mat_aIb(ab) = mat_aIb(ab) + sum j_b (i_aj_b|k)*v(j_b) and
1164 !> mat_bIa(ba) = mat_bIa(ba) + sum i_a (i_aj_b|k)*v(i_a)
1165 !> The block size of the columns of mat_aIb and the rows of mat_bIa are the size of k (RI)
1166 !> \param ab_Q the tensor containing the 3-center integrals
1167 !> \param vec the contraction coefficients
1168 !> \param mat_aIb normal type dbcsr matrix
1169 !> \param mat_bIa normal type dbcsr matrix
1170 !> \param atom_k the atom for which we contract
1171 !> It is assumed that the contraction coefficients for MO I are all on atom_k
1172 !> We do the classic thing when we fill half the matrix and add its transposed to get the full
1173 !> one, but here, the matrix is not symmetric, hence we explicitely have 2 input matrices
1174 !> The distribution of the integrals and the normal dbcsr matrix are compatible out of the box
1175 ! **************************************************************************************************
1176  SUBROUTINE contract2_ao_to_domo_low(ab_Q, vec, mat_aIb, mat_bIa, atom_k)
1177 
1178  TYPE(dbt_type) :: ab_q
1179  REAL(dp), DIMENSION(:), INTENT(IN) :: vec
1180  TYPE(dbcsr_type), INTENT(INOUT) :: mat_aib, mat_bia
1181  INTEGER, INTENT(IN) :: atom_k
1182 
1183  CHARACTER(LEN=*), PARAMETER :: routinen = 'contract2_AO_to_doMO_low'
1184 
1185  INTEGER :: handle, i, iatom, ind(3), j, jatom, &
1186  katom, s1, s2
1187  INTEGER, DIMENSION(:), POINTER :: atom_blk_size
1188  LOGICAL :: found, t_found
1189  REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: iabc
1190  REAL(dp), DIMENSION(:, :), POINTER :: pblock
1191  TYPE(dbt_iterator_type) :: iter
1192 
1193  NULLIFY (atom_blk_size, pblock)
1194 
1195  CALL timeset(routinen, handle)
1196 
1197  CALL dbcsr_get_info(mat_aib, row_blk_size=atom_blk_size)
1198 
1199 !$OMP PARALLEL DEFAULT(NONE) &
1200 !$OMP SHARED(ab_Q,vec,mat_aIb,mat_bIa,atom_k,atom_blk_size) &
1201 !$OMP PRIVATE(iter,ind,iatom,jatom,katom,iabc,t_found,found,s1,s2,j,i,pblock)
1202  CALL dbt_iterator_start(iter, ab_q)
1203  DO WHILE (dbt_iterator_blocks_left(iter))
1204  CALL dbt_iterator_next_block(iter, ind)
1205 
1206  iatom = ind(1)
1207  jatom = ind(2)
1208  katom = ind(3)
1209 
1210  IF (atom_k .NE. katom) cycle
1211 
1212  CALL dbt_get_block(ab_q, ind, iabc, t_found)
1213  IF (.NOT. t_found) cycle
1214 
1215  ! Deal with mat_aIb
1216  IF (jatom == atom_k) THEN
1217  s1 = atom_blk_size(iatom)
1218  s2 = SIZE(iabc, 3)
1219 
1220  CALL dbcsr_get_block_p(matrix=mat_aib, row=iatom, col=jatom, block=pblock, found=found)
1221 
1222  IF (found) THEN
1223  DO i = 1, s1
1224  DO j = 1, s2
1225 !$OMP ATOMIC
1226  pblock(i, j) = pblock(i, j) + dot_product(vec, iabc(i, :, j))
1227  END DO
1228  END DO
1229  END IF
1230  END IF ! jatom == atom_k
1231 
1232  ! Deal with mat_bIa, keep block diagonal empty
1233  IF (iatom == jatom) cycle
1234  IF (iatom == atom_k) THEN
1235  s1 = SIZE(iabc, 3)
1236  s2 = atom_blk_size(jatom)
1237 
1238  CALL dbcsr_get_block_p(matrix=mat_bia, row=iatom, col=jatom, block=pblock, found=found)
1239 
1240  IF (found) THEN
1241  DO i = 1, s1
1242  DO j = 1, s2
1243 !$OMP ATOMIC
1244  pblock(i, j) = pblock(i, j) + dot_product(vec, iabc(:, j, i))
1245  END DO
1246  END DO
1247  END IF
1248  END IF !iatom== atom_k
1249 
1250  DEALLOCATE (iabc)
1251  END DO !iter
1252  CALL dbt_iterator_stop(iter)
1253 !$OMP END PARALLEL
1254 
1255  CALL timestop(handle)
1256 
1257  END SUBROUTINE contract2_ao_to_domo_low
1258 
1259 ! **************************************************************************************************
1260 !> \brief Multiply all the blocks of a contracted RI integral (aI|P) by a matrix of type (P|...|Q)
1261 !> \param contr_int the integral array
1262 !> \param PQ the smaller matrix to multiply all blocks
1263 !> \note It is assumed that all non-zero blocks have the same number of columns. Can pass partial
1264 !> arrays, e.g. contr_int(1:3)
1265 ! **************************************************************************************************
1266  SUBROUTINE ri_all_blocks_mm(contr_int, PQ)
1267 
1268  TYPE(dbcsr_p_type), DIMENSION(:) :: contr_int
1269  REAL(dp), DIMENSION(:, :), INTENT(IN) :: pq
1270 
1271  INTEGER :: blk, iblk, imo, jblk, ndo_mo, s1, s2
1272  LOGICAL :: found
1273  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: work
1274  REAL(dp), DIMENSION(:, :), POINTER :: pblock
1275  TYPE(dbcsr_iterator_type) :: iter
1276 
1277  NULLIFY (pblock)
1278 
1279  ndo_mo = SIZE(contr_int)
1280 
1281  DO imo = 1, ndo_mo
1282  CALL dbcsr_iterator_start(iter, contr_int(imo)%matrix)
1283  DO WHILE (dbcsr_iterator_blocks_left(iter))
1284 
1285  CALL dbcsr_iterator_next_block(iter, row=iblk, column=jblk, blk=blk)
1286  CALL dbcsr_get_block_p(contr_int(imo)%matrix, iblk, jblk, pblock, found)
1287 
1288  IF (found) THEN
1289  s1 = SIZE(pblock, 1)
1290  s2 = SIZE(pblock, 2)
1291  ALLOCATE (work(s1, s2))
1292  CALL dgemm('N', 'N', s1, s2, s2, 1.0_dp, pblock, s1, pq, s2, 0.0_dp, work, s1)
1293  CALL dcopy(s1*s2, work, 1, pblock, 1)
1294  DEALLOCATE (work)
1295  END IF
1296 
1297  END DO ! dbcsr iterator
1298  CALL dbcsr_iterator_stop(iter)
1299  END DO !imo
1300 
1301  END SUBROUTINE ri_all_blocks_mm
1302 
1303 ! **************************************************************************************************
1304 !> \brief Copies an (partial) array of contracted RI integrals into anoter one
1305 !> \param new_int where the copy is stored
1306 !> \param ref_int what is copied
1307 !> \note Allocate the matrices of new_int if not done already
1308 ! **************************************************************************************************
1309  SUBROUTINE copy_ri_contr_int(new_int, ref_int)
1310 
1311  TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: new_int
1312  TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: ref_int
1313 
1314  INTEGER :: iso, ndo_so
1315 
1316  cpassert(SIZE(new_int) == SIZE(ref_int))
1317  ndo_so = SIZE(ref_int)
1318 
1319  DO iso = 1, ndo_so
1320  IF (.NOT. ASSOCIATED(new_int(iso)%matrix)) ALLOCATE (new_int(iso)%matrix)
1321  CALL dbcsr_copy(new_int(iso)%matrix, ref_int(iso)%matrix)
1322  END DO
1323 
1324  END SUBROUTINE copy_ri_contr_int
1325 
1326 ! **************************************************************************************************
1327 !> \brief Takes the product of contracted integrals and put them in a kernel matrix
1328 !> \param kernel the matrix where the products are stored
1329 !> \param lhs_int the left-hand side contracted integrals
1330 !> \param rhs_int the right-hand side contracted integrals
1331 !> \param quadrants on which quadrant(s) on the kernel matrix the product is stored
1332 !> \param qs_env ...
1333 !> \param eps_filter filter for dbcsr matrix multiplication
1334 !> \param mo_transpose whether the MO blocks should be transpose, i.e. (aI|Jb) => (aJ|Ib)
1335 !> \note It is assumed that the kerenl matrix is NOT symmetric
1336 !> There are three quadrants, corresponding to 1: the upper-left (diagonal), 2: the
1337 !> upper-right (off-diagonal) and 3: the lower-right (diagonal).
1338 !> Need to finalize the kernel matrix after calling this routine (possibly multiple times)
1339 ! **************************************************************************************************
1340  SUBROUTINE ri_int_product(kernel, lhs_int, rhs_int, quadrants, qs_env, eps_filter, mo_transpose)
1341 
1342  TYPE(dbcsr_type), INTENT(INOUT) :: kernel
1343  TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: lhs_int, rhs_int
1344  LOGICAL, DIMENSION(3), INTENT(IN) :: quadrants
1345  TYPE(qs_environment_type), POINTER :: qs_env
1346  REAL(dp), INTENT(IN), OPTIONAL :: eps_filter
1347  LOGICAL, INTENT(IN), OPTIONAL :: mo_transpose
1348 
1349  INTEGER :: blk, i, iblk, iso, j, jblk, jso, nblk, &
1350  ndo_so
1351  LOGICAL :: found, my_mt
1352  REAL(dp), DIMENSION(:, :), POINTER :: pblock
1353  TYPE(dbcsr_iterator_type) :: iter
1354  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1355  TYPE(dbcsr_type) :: prod
1356 
1357  NULLIFY (matrix_s, pblock)
1358 
1359 ! Initialization
1360  cpassert(SIZE(lhs_int) == SIZE(rhs_int))
1361  cpassert(any(quadrants))
1362  ndo_so = SIZE(lhs_int)
1363  CALL get_qs_env(qs_env, matrix_s=matrix_s, natom=nblk)
1364  CALL dbcsr_create(prod, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1365  my_mt = .false.
1366  IF (PRESENT(mo_transpose)) my_mt = mo_transpose
1367 
1368  ! The kernel matrix is symmetric (even if normal type) => only fill upper half on diagonal
1369  ! quadrants, but the whole thing on upper-right quadrant
1370  DO iso = 1, ndo_so
1371  DO jso = 1, ndo_so
1372 
1373  ! If on-diagonal quadrants only, can skip jso < iso
1374  IF (.NOT. quadrants(2) .AND. jso < iso) cycle
1375 
1376  i = iso; j = jso;
1377  IF (my_mt) THEN
1378  i = jso; j = iso;
1379  END IF
1380 
1381  ! Take the product lhs*rhs^T
1382  CALL dbcsr_multiply('N', 'T', 1.0_dp, lhs_int(i)%matrix, rhs_int(j)%matrix, &
1383  0.0_dp, prod, filter_eps=eps_filter)
1384 
1385  ! Loop over blocks of prod and fill kernel matrix => ok cuz same (but replicated) dist
1386  CALL dbcsr_iterator_start(iter, prod)
1387  DO WHILE (dbcsr_iterator_blocks_left(iter))
1388 
1389  CALL dbcsr_iterator_next_block(iter, row=iblk, column=jblk, blk=blk)
1390  IF ((iso == jso .AND. jblk < iblk) .AND. .NOT. quadrants(2)) cycle
1391 
1392  CALL dbcsr_get_block_p(prod, iblk, jblk, pblock, found)
1393 
1394  IF (found) THEN
1395 
1396  ! Case study on quadrant
1397  !upper-left
1398  IF (quadrants(1)) THEN
1399  CALL dbcsr_put_block(kernel, (iso - 1)*nblk + iblk, (jso - 1)*nblk + jblk, pblock)
1400  END IF
1401 
1402  !upper-right
1403  IF (quadrants(2)) THEN
1404  CALL dbcsr_put_block(kernel, (iso - 1)*nblk + iblk, (ndo_so + jso - 1)*nblk + jblk, pblock)
1405  END IF
1406 
1407  !lower-right
1408  IF (quadrants(3)) THEN
1409  CALL dbcsr_put_block(kernel, (ndo_so + iso - 1)*nblk + iblk, (ndo_so + jso - 1)*nblk + jblk, pblock)
1410  END IF
1411 
1412  END IF
1413 
1414  END DO ! dbcsr iterator
1415  CALL dbcsr_iterator_stop(iter)
1416 
1417  END DO !jso
1418  END DO !iso
1419 
1420 ! Clean-up
1421  CALL dbcsr_release(prod)
1422 
1423  END SUBROUTINE ri_int_product
1424 
1425 END MODULE xas_tdp_kernel
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
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition: cp_dbcsr_operations.F:89
basis_set_types
Definition: basis_set_types.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
dbt_api
This is the start of a dbt_api, all publically needed functions are exported here....
Definition: dbt_api.F:17
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
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
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
util
All kind of helpful little routines.
Definition: util.F:14
util::get_limit
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Definition: util.F:333
xas_tdp_kernel
All the kernel specific subroutines for XAS TDP calculations.
Definition: xas_tdp_kernel.F:13
xas_tdp_kernel::kernel_coulomb_xc
subroutine, public kernel_coulomb_xc(coul_ker, xc_ker, donor_state, xas_tdp_env, xas_tdp_control, qs_env)
Computes, if asked for it, the Coulomb and XC kernel matrices, in the usuall matrix format.
Definition: xas_tdp_kernel.F:81
xas_tdp_kernel::contract2_ao_to_domo
subroutine, public contract2_ao_to_domo(contr_int, op_type, donor_state, xas_tdp_env, xas_tdp_control, qs_env)
Contract the ri 3-center integrals stored in a tensor with repect to the donor MOs coeffs,...
Definition: xas_tdp_kernel.F:976
xas_tdp_kernel::kernel_exchange
subroutine, public kernel_exchange(ex_ker, donor_state, xas_tdp_env, xas_tdp_control, qs_env)
Computes the exact exchange kernel matrix using RI. Returns an array of 2 matrices,...
Definition: xas_tdp_kernel.F:595
xas_tdp_kernel::reserve_contraction_blocks
subroutine, public reserve_contraction_blocks(matrices, ri_atom, qs_env)
Reserves the blocks in of a dbcsr matrix as needed for RI 3-center contraction (aI|P)
Definition: xas_tdp_kernel.F:903
xas_tdp_kernel::ri_all_blocks_mm
subroutine, public ri_all_blocks_mm(contr_int, PQ)
Multiply all the blocks of a contracted RI integral (aI|P) by a matrix of type (P|....
Definition: xas_tdp_kernel.F:1267
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
xas_tdp_types::get_proc_batch_sizes
subroutine, public get_proc_batch_sizes(batch_size, nbatch, nex_atom, nprocs)
Uses heuristics to determine a good batching of the processros for fxc integration.
Definition: xas_tdp_types.F:1427