(git:98357aa)
Loading...
Searching...
No Matches
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-2026 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
13MODULE xas_tdp_kernel
14 USE basis_set_types, ONLY: gto_basis_set_p_type
15 USE cp_dbcsr_api, ONLY: &
16 dbcsr_add, dbcsr_complete_redistribute, dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, &
17 dbcsr_distribution_get, dbcsr_distribution_new, dbcsr_distribution_release, &
18 dbcsr_distribution_type, dbcsr_filter, dbcsr_finalize, dbcsr_get_block_p, dbcsr_get_info, &
19 dbcsr_get_num_blocks, dbcsr_get_stored_coordinates, dbcsr_iterator_blocks_left, &
20 dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
21 dbcsr_multiply, dbcsr_p_type, dbcsr_put_block, dbcsr_release, dbcsr_reserve_blocks, &
22 dbcsr_set, dbcsr_transposed, dbcsr_type, dbcsr_type_no_symmetry, dbcsr_type_symmetric
23 USE cp_dbcsr_operations, ONLY: cp_dbcsr_dist2d_to_dist,&
24 dbcsr_deallocate_matrix_set
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
58CONTAINS
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
168 IF (do_sf) ub = 4
169 DO i = lb, ub
170 IF (.NOT. ASSOCIATED(xas_tdp_env%ri_fxc(ri_atom, i)%array)) THEN
171 ALLOCATE (xas_tdp_env%ri_fxc(ri_atom, i)%array(nsgfp, nsgfp))
172 END IF
173 CALL para_env%bcast(xas_tdp_env%ri_fxc(ri_atom, i)%array, source)
174 END DO
175
176 xas_tdp_env%fxc_avail = .true.
177 END IF
178
179 ! Case study on the calculation type
180 IF (do_sg .OR. do_tp) THEN
181 CALL rcs_xc(xc_ker(1)%matrix, xc_ker(2)%matrix, contr1_int, dist, blk_size, &
182 donor_state, xas_tdp_env, xas_tdp_control, qs_env)
183 END IF
184
185 IF (do_sc) THEN
186 CALL sc_os_xc(xc_ker(3)%matrix, contr1_int, dist, blk_size, donor_state, &
187 xas_tdp_env, xas_tdp_control, qs_env)
188 END IF
189
190 IF (do_sf) THEN
191 CALL ondiag_sf_os_xc(xc_ker(4)%matrix, contr1_int, dist, blk_size, donor_state, &
192 xas_tdp_env, xas_tdp_control, qs_env)
193 END IF
194
195 END IF ! do_xc
196
197! Clean-up
198 CALL dbcsr_deallocate_matrix_set(contr1_int)
199
200 CALL timestop(handle)
201
202 END SUBROUTINE kernel_coulomb_xc
203
204! **************************************************************************************************
205!> \brief Create the matrix containing the Coulomb kernel, which is:
206!> (aI_sigma|J_tau b) ~= (aI_sigma|P) * (P|Q) * (Q|J_tau b)
207!> \param coul_ker the Coulomb kernel
208!> \param contr1_int the once contracted RI integrals (aI|P)
209!> \param dist the inherited dbcsr ditribution
210!> \param blk_size the inherited block sizes
211!> \param xas_tdp_env ...
212!> \param xas_tdp_control ...
213!> \param qs_env ...
214! **************************************************************************************************
215 SUBROUTINE coulomb(coul_ker, contr1_int, dist, blk_size, xas_tdp_env, xas_tdp_control, qs_env)
216
217 TYPE(dbcsr_type), INTENT(INOUT) :: coul_ker
218 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: contr1_int
219 TYPE(dbcsr_distribution_type), POINTER :: dist
220 INTEGER, DIMENSION(:), POINTER :: blk_size
221 TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
222 TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
223 TYPE(qs_environment_type), POINTER :: qs_env
224
225 LOGICAL :: quadrants(3)
226 REAL(dp), DIMENSION(:, :), POINTER :: pq
227 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: lhs_int, rhs_int
228 TYPE(dbcsr_type) :: work_mat
229
230 NULLIFY (pq, rhs_int, lhs_int)
231
232 ! Get the inver RI coulomb
233 pq => xas_tdp_env%ri_inv_coul
234
235 ! Create a normal type work matrix
236 CALL dbcsr_create(work_mat, name="WORK", matrix_type=dbcsr_type_no_symmetry, dist=dist, &
237 row_blk_size=blk_size, col_blk_size=blk_size)
238
239 ! Compute the product (aI|P) * (P|Q)^-1 * (Q|Jb) = (aI|Jb)
240 rhs_int => contr1_int ! the incoming contr1_int is not modified
241 ALLOCATE (lhs_int(SIZE(contr1_int)))
242 CALL copy_ri_contr_int(lhs_int, rhs_int) ! RHS containts (Q|JB)^T
243 CALL ri_all_blocks_mm(lhs_int, pq) ! LHS contatins (aI|P)*(P|Q)^-1
244
245 !In the special case of ROKS, same MOs for each spin => put same (aI|Jb) product on the
246 !alpha-alpha, alpha-beta and beta-beta quadrants of the kernel matrix.
247 IF (xas_tdp_control%do_roks) THEN
248 quadrants = [.true., .true., .true.]
249 ELSE
250 quadrants = [.true., .false., .false.]
251 END IF
252 CALL ri_int_product(work_mat, lhs_int, rhs_int, quadrants, qs_env, &
253 eps_filter=xas_tdp_control%eps_filter)
254 CALL dbcsr_finalize(work_mat)
255
256 !Create the symmetric kernel matrix and redistribute work_mat into it
257 CALL dbcsr_create(coul_ker, name="COULOMB KERNEL", matrix_type=dbcsr_type_symmetric, dist=dist, &
258 row_blk_size=blk_size, col_blk_size=blk_size)
259 CALL dbcsr_complete_redistribute(work_mat, coul_ker)
260
261 !clean-up
262 CALL dbcsr_release(work_mat)
263 CALL dbcsr_deallocate_matrix_set(lhs_int)
264
265 END SUBROUTINE coulomb
266
267! **************************************************************************************************
268!> \brief Create the matrix containing the XC kenrel in the spin-conserving open-shell case:
269!> (aI_sigma|fxc|J_tau b) ~= (aI_sigma|P) (P|Q)^-1 (Q|fxc|R) (R|S)^-1 (S|J_tau b)
270!> \param xc_ker the kernel matrix
271!> \param contr1_int_PQ the once contracted RI integrals, with inverse coulomb: (aI_sigma|P) (P|Q)^-1
272!> \param dist inherited dbcsr dist
273!> \param blk_size inherited block sizes
274!> \param donor_state ...
275!> \param xas_tdp_env ...
276!> \param xas_tdp_control ...
277!> \param qs_env ...
278!> note Prior to calling this function, the (Q|fxc|R) integral must be brodcasted to all procs
279! **************************************************************************************************
280 SUBROUTINE sc_os_xc(xc_ker, contr1_int_PQ, dist, blk_size, donor_state, xas_tdp_env, &
281 xas_tdp_control, qs_env)
282
283 TYPE(dbcsr_type), INTENT(INOUT) :: xc_ker
284 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: contr1_int_pq
285 TYPE(dbcsr_distribution_type), POINTER :: dist
286 INTEGER, DIMENSION(:), POINTER :: blk_size
287 TYPE(donor_state_type), POINTER :: donor_state
288 TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
289 TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
290 TYPE(qs_environment_type), POINTER :: qs_env
291
292 INTEGER :: ndo_mo, ri_atom
293 LOGICAL :: quadrants(3)
294 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: lhs_int, rhs_int
295 TYPE(dbcsr_type) :: work_mat
296
297 NULLIFY (lhs_int, rhs_int)
298
299 ! Initialization
300 ndo_mo = donor_state%ndo_mo
301 ri_atom = donor_state%at_index
302 !normal type work matrix such that distribution of all spin quadrants match
303 CALL dbcsr_create(work_mat, name="WORK", matrix_type=dbcsr_type_no_symmetry, dist=dist, &
304 row_blk_size=blk_size, col_blk_size=blk_size)
305
306 rhs_int => contr1_int_pq ! contains [ (aI|P)*(P|Q)^-1 ]^T
307 ALLOCATE (lhs_int(SIZE(contr1_int_pq))) ! will contain (aI|P)*(P|Q)^-1 * (Q|fxc|R)
308
309 ! Case study: UKS or ROKS ?
310 IF (xas_tdp_control%do_uks) THEN
311
312 ! In the case of UKS, donor MOs might be different for different spins. Moreover, the
313 ! fxc itself might change since fxc = fxc_sigma,tau
314 ! => Carfully treat each spin-quadrant separately
315
316 ! alpha-alpha spin quadrant (upper-lefet)
317 quadrants = [.true., .false., .false.]
318
319 ! Copy the alpha part into lhs_int, multiply by the alpha-alpha (Q|fxc|R) and then
320 ! by the alpha part of rhs_int
321 CALL copy_ri_contr_int(lhs_int(1:ndo_mo), rhs_int(1:ndo_mo))
322 CALL ri_all_blocks_mm(lhs_int(1:ndo_mo), xas_tdp_env%ri_fxc(ri_atom, 1)%array)
323 CALL ri_int_product(work_mat, lhs_int(1:ndo_mo), rhs_int(1:ndo_mo), quadrants, qs_env, &
324 eps_filter=xas_tdp_control%eps_filter)
325
326 ! alpha-beta spin quadrant (upper-right)
327 quadrants = [.false., .true., .false.]
328
329 !Copy the alpha part into LHS, multiply by the alpha-beta kernel and the beta part of RHS
330 CALL copy_ri_contr_int(lhs_int(1:ndo_mo), rhs_int(1:ndo_mo))
331 CALL ri_all_blocks_mm(lhs_int(1:ndo_mo), xas_tdp_env%ri_fxc(ri_atom, 2)%array)
332 CALL ri_int_product(work_mat, lhs_int(1:ndo_mo), rhs_int(ndo_mo + 1:2*ndo_mo), &
333 quadrants, qs_env, eps_filter=xas_tdp_control%eps_filter)
334
335 ! beta-beta spin quadrant (lower-right)
336 quadrants = [.false., .false., .true.]
337
338 !Copy the beta part into LHS, multiply by the beta-beta kernel and the beta part of RHS
339 CALL copy_ri_contr_int(lhs_int(ndo_mo + 1:2*ndo_mo), rhs_int(ndo_mo + 1:2*ndo_mo))
340 CALL ri_all_blocks_mm(lhs_int(ndo_mo + 1:2*ndo_mo), xas_tdp_env%ri_fxc(ri_atom, 3)%array)
341 CALL ri_int_product(work_mat, lhs_int(ndo_mo + 1:2*ndo_mo), rhs_int(ndo_mo + 1:2*ndo_mo), &
342 quadrants, qs_env, eps_filter=xas_tdp_control%eps_filter)
343
344 ELSE IF (xas_tdp_control%do_roks) THEN
345
346 ! In the case of ROKS, fxc = fxc_sigma,tau is different for each spin quadrant, but the
347 ! donor MOs remain the same
348
349 ! alpha-alpha kernel in the upper left quadrant
350 quadrants = [.true., .false., .false.]
351
352 !Copy the LHS and multiply by alpha-alpha kernel
353 CALL copy_ri_contr_int(lhs_int, rhs_int)
354 CALL ri_all_blocks_mm(lhs_int, xas_tdp_env%ri_fxc(ri_atom, 1)%array)
355 CALL ri_int_product(work_mat, lhs_int, rhs_int, quadrants, qs_env, &
356 eps_filter=xas_tdp_control%eps_filter)
357
358 ! alpha-beta kernel in the upper-right quadrant
359 quadrants = [.false., .true., .false.]
360
361 !Copy LHS and multiply by the alpha-beta kernel
362 CALL copy_ri_contr_int(lhs_int, rhs_int)
363 CALL ri_all_blocks_mm(lhs_int, xas_tdp_env%ri_fxc(ri_atom, 2)%array)
364 CALL ri_int_product(work_mat, lhs_int, rhs_int, quadrants, qs_env, &
365 eps_filter=xas_tdp_control%eps_filter)
366
367 ! beta-beta kernel in the lower-right quadrant
368 quadrants = [.false., .false., .true.]
369
370 !Copy the LHS and multiply by the beta-beta kernel
371 CALL copy_ri_contr_int(lhs_int, rhs_int)
372 CALL ri_all_blocks_mm(lhs_int, xas_tdp_env%ri_fxc(ri_atom, 3)%array)
373 CALL ri_int_product(work_mat, lhs_int, rhs_int, quadrants, qs_env, &
374 eps_filter=xas_tdp_control%eps_filter)
375
376 END IF
377 CALL dbcsr_finalize(work_mat)
378
379 ! Create a symmetric kernel matrix and redistribute the normal work matrix into it
380 CALL dbcsr_create(xc_ker, name="SC OS XC KERNEL", matrix_type=dbcsr_type_symmetric, dist=dist, &
381 row_blk_size=blk_size, col_blk_size=blk_size)
382 CALL dbcsr_complete_redistribute(work_mat, xc_ker)
383
384 !clean-up
385 CALL dbcsr_deallocate_matrix_set(lhs_int)
386 CALL dbcsr_release(work_mat)
387
388 END SUBROUTINE sc_os_xc
389
390! **************************************************************************************************
391!> \brief Create the matrix containing the on-diagonal spin-flip XC kernel (open-shell), which is:
392!> (a I_sigma|fxc|J_tau b) * delta_sigma,tau, fxc = 1/(rhoa-rhob) * (dE/drhoa - dE/drhob)
393!> 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)
394!> \param xc_ker the kernel matrix
395!> \param contr1_int_PQ the once contracted RI integrals with coulomb product: (aI_sigma|P) (P|Q)^-1
396!> \param dist inherited dbcsr dist
397!> \param blk_size inherited block sizes
398!> \param donor_state ...
399!> \param xas_tdp_env ...
400!> \param xas_tdp_control ...
401!> \param qs_env ...
402!> \note It must be later on multiplied by the spin-swapped Q projector
403!> Prior to calling this function, the (Q|fxc|R) integral must be brodcasted to all procs
404! **************************************************************************************************
405 SUBROUTINE ondiag_sf_os_xc(xc_ker, contr1_int_PQ, dist, blk_size, donor_state, xas_tdp_env, &
406 xas_tdp_control, qs_env)
407
408 TYPE(dbcsr_type), INTENT(INOUT) :: xc_ker
409 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: contr1_int_pq
410 TYPE(dbcsr_distribution_type), POINTER :: dist
411 INTEGER, DIMENSION(:), POINTER :: blk_size
412 TYPE(donor_state_type), POINTER :: donor_state
413 TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
414 TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
415 TYPE(qs_environment_type), POINTER :: qs_env
416
417 INTEGER :: ndo_mo, ri_atom
418 LOGICAL :: quadrants(3)
419 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: lhs_int, rhs_int
420 TYPE(dbcsr_type) :: work_mat
421
422 NULLIFY (lhs_int, rhs_int)
423
424 ! Initialization
425 ndo_mo = donor_state%ndo_mo
426 ri_atom = donor_state%at_index
427 !normal type work matrix such that distribution of all spin quadrants match
428 CALL dbcsr_create(work_mat, name="WORK", matrix_type=dbcsr_type_no_symmetry, dist=dist, &
429 row_blk_size=blk_size, col_blk_size=blk_size)
430
431 !Create a lhs_int, in which the whole (aI_sigma|P) (P|Q)^-1 (Q|fxc|R) will be put
432 !because in spin-flip fxc is spin-independent, can take the product once and for all
433 rhs_int => contr1_int_pq
434 ALLOCATE (lhs_int(SIZE(contr1_int_pq)))
435 CALL copy_ri_contr_int(lhs_int, rhs_int)
436 CALL ri_all_blocks_mm(lhs_int, xas_tdp_env%ri_fxc(ri_atom, 4)%array)
437
438 ! Case study: UKS or ROKS ?
439 IF (xas_tdp_control%do_uks) THEN
440
441 ! In the case of UKS, donor MOs might be different for different spins
442 ! => Carfully treat each spin-quadrant separately
443 ! NO alpha-beta because of the delta_sigma,tau
444
445 ! alpha-alpha spin quadrant (upper-lefet)
446 quadrants = [.true., .false., .false.]
447 CALL ri_int_product(work_mat, lhs_int(1:ndo_mo), rhs_int(1:ndo_mo), quadrants, qs_env, &
448 eps_filter=xas_tdp_control%eps_filter)
449
450 ! beta-beta spin quadrant (lower-right)
451 quadrants = [.false., .false., .true.]
452 CALL ri_int_product(work_mat, lhs_int(ndo_mo + 1:2*ndo_mo), rhs_int(ndo_mo + 1:2*ndo_mo), &
453 quadrants, qs_env, eps_filter=xas_tdp_control%eps_filter)
454
455 ELSE IF (xas_tdp_control%do_roks) THEN
456
457 ! In the case of ROKS, same donor MOs for both spins => can do it all at once
458 ! But NOT the alpha-beta quadrant because of delta_sigma,tau
459
460 quadrants = [.true., .false., .true.]
461 CALL ri_int_product(work_mat, lhs_int, rhs_int, quadrants, qs_env, &
462 eps_filter=xas_tdp_control%eps_filter)
463
464 END IF
465 CALL dbcsr_finalize(work_mat)
466
467 ! Create a symmetric kernel matrix and redistribute the normal work matrix into it
468 CALL dbcsr_create(xc_ker, name="ON-DIAG SF OS XC KERNEL", matrix_type=dbcsr_type_symmetric, &
469 dist=dist, row_blk_size=blk_size, col_blk_size=blk_size)
470 CALL dbcsr_complete_redistribute(work_mat, xc_ker)
471
472 !clean-up
473 CALL dbcsr_deallocate_matrix_set(lhs_int)
474 CALL dbcsr_release(work_mat)
475
476 END SUBROUTINE ondiag_sf_os_xc
477
478! **************************************************************************************************
479!> \brief Create the matrix containing the XC kernel in the restricted closed-shell case, for
480!> singlets: (aI|fxc|Jb) = (aI|P) (P|Q)^-1 (Q|fxc_alp,alp+fxc_alp,bet|R) (R|S)^-1 (S|Jb)
481!> triplets: (aI|fxc|Jb) = (aI|P) (P|Q)^-1 (Q|fxc_alp,alp-fxc_alp,bet|R) (R|S)^-1 (S|Jb)
482!> \param sg_xc_ker the singlet kernel matrix
483!> \param tp_xc_ker the triplet kernel matrix
484!> \param contr1_int_PQ the once contracted RI integrals including inverse 2-center Coulomb prodcut:
485!> (aI|P)*(P|Q)^-1
486!> \param dist inherited dbcsr dist
487!> \param blk_size inherited block sizes
488!> \param donor_state ...
489!> \param xas_tdp_env ...
490!> \param xas_tdp_control ...
491!> \param qs_env ...
492! **************************************************************************************************
493 SUBROUTINE rcs_xc(sg_xc_ker, tp_xc_ker, contr1_int_PQ, dist, blk_size, donor_state, &
494 xas_tdp_env, xas_tdp_control, qs_env)
495
496 TYPE(dbcsr_type), INTENT(INOUT) :: sg_xc_ker, tp_xc_ker
497 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: contr1_int_pq
498 TYPE(dbcsr_distribution_type), POINTER :: dist
499 INTEGER, DIMENSION(:), POINTER :: blk_size
500 TYPE(donor_state_type), POINTER :: donor_state
501 TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
502 TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
503 TYPE(qs_environment_type), POINTER :: qs_env
504
505 INTEGER :: nsgfp, ri_atom
506 LOGICAL :: quadrants(3)
507 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: fxc
508 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: lhs_int, rhs_int
509 TYPE(dbcsr_type) :: work_mat
510
511 NULLIFY (lhs_int, rhs_int)
512
513 ! Initialization
514 ri_atom = donor_state%at_index
515 nsgfp = SIZE(xas_tdp_env%ri_fxc(ri_atom, 1)%array, 1)
516 rhs_int => contr1_int_pq ! RHS contains [ (aI|P)*(P|Q)^-1 ]^T
517 ALLOCATE (lhs_int(SIZE(contr1_int_pq))) ! LHS will contatin (aI|P)*(P|Q)^-1 * (Q|fxc|R)
518
519 ! Work structures
520 ALLOCATE (fxc(nsgfp, nsgfp))
521 CALL dbcsr_create(work_mat, name="WORK", matrix_type=dbcsr_type_no_symmetry, dist=dist, &
522 row_blk_size=blk_size, col_blk_size=blk_size)
523
524 ! Case study: singlet and/or triplet ?
525 IF (xas_tdp_control%do_singlet) THEN
526
527 ! Take the sum of fxc for alpha-alpha and alpha-beta
528 CALL dcopy(nsgfp*nsgfp, xas_tdp_env%ri_fxc(ri_atom, 1)%array, 1, fxc, 1)
529 CALL daxpy(nsgfp*nsgfp, 1.0_dp, xas_tdp_env%ri_fxc(ri_atom, 2)%array, 1, fxc, 1)
530
531 ! Copy the fresh lhs_int = (aI|P) (P|Q)^-1 and multiply by (Q|fxc|R)
532 CALL copy_ri_contr_int(lhs_int, rhs_int)
533 CALL ri_all_blocks_mm(lhs_int, fxc)
534
535 ! Compute the final LHS RHS product => spin-restricted, only upper-left quadrant
536 quadrants = [.true., .false., .false.]
537 CALL ri_int_product(work_mat, lhs_int, rhs_int, quadrants, qs_env, &
538 eps_filter=xas_tdp_control%eps_filter)
539 CALL dbcsr_finalize(work_mat)
540
541 !Create the symmetric kernel matrix and redistribute work_mat into it
542 CALL dbcsr_create(sg_xc_ker, name="XC SINGLET KERNEL", matrix_type=dbcsr_type_symmetric, &
543 dist=dist, row_blk_size=blk_size, col_blk_size=blk_size)
544 CALL dbcsr_complete_redistribute(work_mat, sg_xc_ker)
545
546 END IF
547
548 IF (xas_tdp_control%do_triplet) THEN
549
550 ! Take the difference of fxc for alpha-alpha and alpha-beta
551 CALL dcopy(nsgfp*nsgfp, xas_tdp_env%ri_fxc(ri_atom, 1)%array, 1, fxc, 1)
552 CALL daxpy(nsgfp*nsgfp, -1.0_dp, xas_tdp_env%ri_fxc(ri_atom, 2)%array, 1, fxc, 1)
553
554 ! Copy the fresh lhs_int = (aI|P) (P|Q)^-1 and multiply by (Q|fxc|R)
555 CALL copy_ri_contr_int(lhs_int, rhs_int)
556 CALL ri_all_blocks_mm(lhs_int, fxc)
557
558 ! Compute the final LHS RHS product => spin-restricted, only upper-left quadrant
559 quadrants = [.true., .false., .false.]
560 CALL ri_int_product(work_mat, lhs_int, rhs_int, quadrants, qs_env, &
561 eps_filter=xas_tdp_control%eps_filter)
562 CALL dbcsr_finalize(work_mat)
563
564 !Create the symmetric kernel matrix and redistribute work_mat into it
565 CALL dbcsr_create(tp_xc_ker, name="XC TRIPLET KERNEL", matrix_type=dbcsr_type_symmetric, &
566 dist=dist, row_blk_size=blk_size, col_blk_size=blk_size)
567 CALL dbcsr_complete_redistribute(work_mat, tp_xc_ker)
568
569 END IF
570
571 ! clean-up
572 CALL dbcsr_deallocate_matrix_set(lhs_int)
573 CALL dbcsr_release(work_mat)
574 DEALLOCATE (fxc)
575
576 END SUBROUTINE rcs_xc
577
578! **************************************************************************************************
579!> \brief Computes the exact exchange kernel matrix using RI. Returns an array of 2 matrices,
580!> which are:
581!> 1) the on-diagonal kernel: (ab|I_sigma J_tau) * delta_sigma,tau
582!> 2) the off-diagonal spin-conserving kernel: (aJ_sigma|I_tau b) * delta_sigma,tau
583!> An internal analysis determines which of the above are computed (can range from 0 to 2),
584!> \param ex_ker ...
585!> \param donor_state ...
586!> \param xas_tdp_env ...
587!> \param xas_tdp_control ...
588!> \param qs_env ...
589!> \note In the case of spin-conserving excitation, the kernel must later be multiplied by the
590!> usual Q projector. In the case of spin-flip, one needs to project the excitations coming
591!> from alpha donor MOs on the unoccupied beta MOs. This is done by multiplying by a Q
592!> projector where the alpha-alpha and beta-beta quadrants are swapped
593!> The ex_ker array should be allocated on entry (not the internals)
594! **************************************************************************************************
595 SUBROUTINE kernel_exchange(ex_ker, donor_state, xas_tdp_env, xas_tdp_control, qs_env)
596
597 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ex_ker
598 TYPE(donor_state_type), POINTER :: donor_state
599 TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
600 TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
601 TYPE(qs_environment_type), POINTER :: qs_env
602
603 CHARACTER(len=*), PARAMETER :: routinen = 'kernel_exchange'
604
605 INTEGER :: handle
606 INTEGER, DIMENSION(:), POINTER :: blk_size
607 LOGICAL :: do_off_sc
608 TYPE(dbcsr_distribution_type), POINTER :: dist
609 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: contr1_int
610
611 NULLIFY (contr1_int, dist, blk_size)
612
613 !Don't do anything if no hfx
614 IF (.NOT. xas_tdp_control%do_hfx) RETURN
615
616 CALL timeset(routinen, handle)
617
618 dist => donor_state%dbcsr_dist
619 blk_size => donor_state%blk_size
620
621 !compute the off-diag spin-conserving only if not TDA and anything that is spin-conserving
622 do_off_sc = (.NOT. xas_tdp_control%tamm_dancoff) .AND. &
623 (xas_tdp_control%do_spin_cons .OR. xas_tdp_control%do_singlet .OR. xas_tdp_control%do_triplet)
624
625 ! Need the once contracted integrals (aI|P)
626 CALL contract2_ao_to_domo(contr1_int, "EXCHANGE", donor_state, xas_tdp_env, xas_tdp_control, qs_env)
627
628! The on-diagonal exchange : (ab|P) * (P|Q)^-1 * (Q|I_sigma J_tau) * delta_sigma,tau
629 CALL ondiag_ex(ex_ker(1)%matrix, contr1_int, dist, blk_size, donor_state, xas_tdp_env, &
630 xas_tdp_control, qs_env)
631
632! The off-diag spin-conserving case: (aJ_sigma|P) * (P|Q)^-1 * (Q|I_tau b) * delta_sigma,tau
633 IF (do_off_sc) THEN
634 CALL offdiag_ex_sc(ex_ker(2)%matrix, contr1_int, dist, blk_size, donor_state, &
635 xas_tdp_env, xas_tdp_control, qs_env)
636 END IF
637
638 !clean-up
639 CALL dbcsr_deallocate_matrix_set(contr1_int)
640
641 CALL timestop(handle)
642
643 END SUBROUTINE kernel_exchange
644
645! **************************************************************************************************
646!> \brief Create the matrix containing the on-diagonal exact exchange kernel, which is:
647!> (ab|I_sigma J_tau) * delta_sigma,tau, where a,b are AOs, I_sigma and J_tau are the donor
648!> spin-orbitals. A RI is done: (ab|I_sigma J_tau) = (ab|P) * (P|Q)^-1 * (Q|I_sigma J_tau)
649!> \param ondiag_ex_ker the on-diagonal exchange kernel in dbcsr format
650!> \param contr1_int the already once-contracted RI 3-center integrals (aI_sigma|P)
651!> where each matrix of the array contains the contraction for the donor spin-orbital I_sigma
652!> \param dist the inherited dbcsr distribution
653!> \param blk_size the inherited dbcsr block sizes
654!> \param donor_state ...
655!> \param xas_tdp_env ...
656!> \param xas_tdp_control ...
657!> \param qs_env ...
658!> \note In the presence of a RI metric, we have instead M^-1 * (P|Q) * M^-1
659! **************************************************************************************************
660 SUBROUTINE ondiag_ex(ondiag_ex_ker, contr1_int, dist, blk_size, donor_state, xas_tdp_env, &
661 xas_tdp_control, qs_env)
662
663 TYPE(dbcsr_type), INTENT(INOUT) :: ondiag_ex_ker
664 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: contr1_int
665 TYPE(dbcsr_distribution_type), POINTER :: dist
666 INTEGER, DIMENSION(:), POINTER :: blk_size
667 TYPE(donor_state_type), POINTER :: donor_state
668 TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
669 TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
670 TYPE(qs_environment_type), POINTER :: qs_env
671
672 INTEGER :: group, iblk, iso, jblk, jso, nblk, &
673 ndo_mo, ndo_so, nsgfa, nsgfp, ri_atom, &
674 source
675 INTEGER, DIMENSION(:), POINTER :: col_dist, col_dist_work, row_dist, &
676 row_dist_work
677 INTEGER, DIMENSION(:, :), POINTER :: pgrid
678 LOGICAL :: do_roks, do_uks, found
679 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: coeffs, ri_coeffs
680 REAL(dp), DIMENSION(:, :), POINTER :: aiq, pblock, pq
681 TYPE(dbcsr_distribution_type) :: opt_dbcsr_dist, work_dbcsr_dist
682 TYPE(dbcsr_iterator_type) :: iter
683 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
684 TYPE(dbcsr_type) :: abij, mats_desymm, work_mat
685 TYPE(mp_para_env_type), POINTER :: para_env
686
687 NULLIFY (para_env, matrix_s, pblock, aiq, row_dist, col_dist, row_dist_work, col_dist_work, pgrid)
688
689 ! We want to compute (ab|I_sigma J_tau) = (ab|P) * (P|Q)^-1 * (Q|I_sigma J_tau)
690 ! Already have (cJ_tau|P) stored in contr1_int. Need to further contract the
691 ! AOs with the coeff of the I_alpha spin-orbital.
692
693 ! Initialization
694 ndo_mo = donor_state%ndo_mo
695 ri_atom = donor_state%at_index
696 do_roks = xas_tdp_control%do_roks
697 do_uks = xas_tdp_control%do_uks
698 ndo_so = ndo_mo; IF (do_uks) ndo_so = 2*ndo_mo !if not UKS, same donor MOs for both spins
699 pq => xas_tdp_env%ri_inv_ex
700
701 CALL get_qs_env(qs_env, para_env=para_env, matrix_s=matrix_s, natom=nblk)
702 nsgfp = SIZE(pq, 1)
703 nsgfa = SIZE(donor_state%contract_coeffs, 1)
704 ALLOCATE (coeffs(nsgfp, ndo_so), ri_coeffs(nsgfp, ndo_so))
705
706 ! a and b need to overlap for non-zero (ab|IJ) => same block structure as overlap S
707 ! need compatible distribution_2d with 3c tensor + normal type
708 CALL cp_dbcsr_dist2d_to_dist(xas_tdp_env%opt_dist2d_ex, opt_dbcsr_dist)
709
710 CALL dbcsr_desymmetrize(matrix_s(1)%matrix, mats_desymm)
711
712 CALL dbcsr_create(abij, template=mats_desymm, name="(ab|IJ)", dist=opt_dbcsr_dist)
713 CALL dbcsr_complete_redistribute(mats_desymm, abij)
714
715 CALL dbcsr_release(mats_desymm)
716
717 ! Create a work distribution based on opt_dbcsr_dist, but for full size matrices
718 CALL dbcsr_distribution_get(opt_dbcsr_dist, row_dist=row_dist, col_dist=col_dist, group=group, &
719 pgrid=pgrid)
720
721 ALLOCATE (row_dist_work(ndo_so*nblk))
722 ALLOCATE (col_dist_work(ndo_so*nblk))
723 DO iso = 1, ndo_so
724 row_dist_work((iso - 1)*nblk + 1:iso*nblk) = row_dist(:)
725 col_dist_work((iso - 1)*nblk + 1:iso*nblk) = col_dist(:)
726 END DO
727
728 CALL dbcsr_distribution_new(work_dbcsr_dist, group=group, pgrid=pgrid, row_dist=row_dist_work, &
729 col_dist=col_dist_work)
730
731 CALL dbcsr_create(work_mat, name="WORK", matrix_type=dbcsr_type_no_symmetry, dist=work_dbcsr_dist, &
732 row_blk_size=blk_size, col_blk_size=blk_size)
733
734 ! Loop over donor spin-orbitals. End matrix is symmetric => span only upper half
735 DO iso = 1, ndo_so
736
737 ! take the (aI|Q) block in contr1_int that has a centered on the excited atom
738 CALL dbcsr_get_stored_coordinates(contr1_int(iso)%matrix, ri_atom, ri_atom, source)
739 IF (para_env%mepos == source) THEN
740 CALL dbcsr_get_block_p(contr1_int(iso)%matrix, ri_atom, ri_atom, aiq, found)
741 ELSE
742 ALLOCATE (aiq(nsgfa, nsgfp))
743 END IF
744 CALL para_env%bcast(aiq, source)
745
746 ! get the contraction (Q|IJ) by taking (Q|Ia)*contract_coeffs and put it in coeffs
747 CALL dgemm('T', 'N', nsgfp, ndo_so, nsgfa, 1.0_dp, aiq, nsgfa, donor_state%contract_coeffs, &
748 nsgfa, 0.0_dp, coeffs, nsgfp)
749
750 ! take (P|Q)^-1 * (Q|IJ) and put that in ri_coeffs
751 CALL dgemm('N', 'N', nsgfp, ndo_so, nsgfp, 1.0_dp, pq, nsgfp, coeffs, nsgfp, 0.0_dp, &
752 ri_coeffs, nsgfp)
753
754 IF (.NOT. para_env%mepos == source) DEALLOCATE (aiq)
755
756 DO jso = iso, ndo_so
757
758 ! There is no alpha-beta exchange. In case of UKS, iso,jso span all spin-orbitals
759 ! => CYCLE if iso and jso are indexing MOs with different spin (and we have UKS)
760 IF (do_uks .AND. (iso <= ndo_mo .AND. jso > ndo_mo)) cycle
761
762 ! compute (ab|IJ) = sum_P (ab|P) * (P|Q)^-1 * (Q|IJ)
763 CALL dbcsr_set(abij, 0.0_dp)
764 CALL contract3_ri_to_domos(xas_tdp_env%ri_3c_ex, ri_coeffs(:, jso), abij, ri_atom)
765
766 ! Loop over (ab|IJ) and copy into work. OK because dist are made to match
767 CALL dbcsr_iterator_start(iter, abij)
768 DO WHILE (dbcsr_iterator_blocks_left(iter))
769
770 CALL dbcsr_iterator_next_block(iter, row=iblk, column=jblk)
771 IF (iso == jso .AND. jblk < iblk) cycle
772
773 CALL dbcsr_get_block_p(abij, iblk, jblk, pblock, found)
774
775 IF (found) THEN
776 CALL dbcsr_put_block(work_mat, (iso - 1)*nblk + iblk, (jso - 1)*nblk + jblk, pblock)
777
778 !In case of ROKS, we have (ab|IJ) for alpha-alpha spin, but it is the same for
779 !beta-beta => replicate the blocks (alpha-beta is zero)
780 IF (do_roks) THEN
781 !the beta-beta block
782 CALL dbcsr_put_block(work_mat, (ndo_so + iso - 1)*nblk + iblk, &
783 (ndo_so + jso - 1)*nblk + jblk, pblock)
784 END IF
785 END IF
786
787 END DO !iterator
788 CALL dbcsr_iterator_stop(iter)
789
790 END DO !jso
791 END DO !iso
792
793 CALL dbcsr_finalize(work_mat)
794 CALL dbcsr_create(ondiag_ex_ker, name="ONDIAG EX KERNEL", matrix_type=dbcsr_type_symmetric, &
795 dist=dist, row_blk_size=blk_size, col_blk_size=blk_size)
796 CALL dbcsr_complete_redistribute(work_mat, ondiag_ex_ker)
797
798 !Clean-up
799 CALL dbcsr_release(work_mat)
800 CALL dbcsr_release(abij)
801 CALL dbcsr_distribution_release(opt_dbcsr_dist)
802 CALL dbcsr_distribution_release(work_dbcsr_dist)
803 DEALLOCATE (col_dist_work, row_dist_work)
804
805 END SUBROUTINE ondiag_ex
806
807! **************************************************************************************************
808!> \brief Create the matrix containing the off-diagonal exact exchange kernel in the spin-conserving
809!> case (which also includes excitations from the closed=shell ref state ) This matrix reads:
810!> (aJ_sigma|I_tau b) * delta_sigma,tau , where a, b are AOs and J_sigma, I_tau are the donor
811!> spin-orbital. A RI is done: (aJ_sigma|I_tau b) = (aJ_sigma|P) * (P|Q)^-1 * (Q|I_tau b)
812!> \param offdiag_ex_ker the off-diagonal, spin-conserving exchange kernel in dbcsr format
813!> \param contr1_int the once-contracted RI integrals: (aJ_sigma|P)
814!> \param dist the inherited dbcsr ditribution
815!> \param blk_size the inherited block sizes
816!> \param donor_state ...
817!> \param xas_tdp_env ...
818!> \param xas_tdp_control ...
819!> \param qs_env ...
820! **************************************************************************************************
821 SUBROUTINE offdiag_ex_sc(offdiag_ex_ker, contr1_int, dist, blk_size, donor_state, xas_tdp_env, &
822 xas_tdp_control, qs_env)
823
824 TYPE(dbcsr_type), INTENT(INOUT) :: offdiag_ex_ker
825 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: contr1_int
826 TYPE(dbcsr_distribution_type), POINTER :: dist
827 INTEGER, DIMENSION(:), POINTER :: blk_size
828 TYPE(donor_state_type), POINTER :: donor_state
829 TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
830 TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
831 TYPE(qs_environment_type), POINTER :: qs_env
832
833 INTEGER :: ndo_mo
834 LOGICAL :: do_roks, do_uks, quadrants(3)
835 REAL(dp), DIMENSION(:, :), POINTER :: pq
836 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: lhs_int, rhs_int
837 TYPE(dbcsr_type) :: work_mat
838
839 NULLIFY (pq, lhs_int, rhs_int)
840
841 !Initialization
842 ndo_mo = donor_state%ndo_mo
843 do_roks = xas_tdp_control%do_roks
844 do_uks = xas_tdp_control%do_uks
845 pq => xas_tdp_env%ri_inv_ex
846
847 rhs_int => contr1_int
848 ALLOCATE (lhs_int(SIZE(contr1_int)))
849 CALL copy_ri_contr_int(lhs_int, rhs_int)
850 CALL ri_all_blocks_mm(lhs_int, pq)
851
852 !Given the lhs_int and rhs_int, all we need to do is multiply elements from the former by
853 !the transpose of the later, and put the result in the correct spin quadrants
854
855 !Create a normal type work matrix
856 CALL dbcsr_create(work_mat, name="WORK", matrix_type=dbcsr_type_no_symmetry, dist=dist, &
857 row_blk_size=blk_size, col_blk_size=blk_size)
858
859 !Case study on closed-shell, ROKS or UKS
860 IF (do_roks) THEN
861 !In ROKS, the donor MOs for each spin are the same => copy the product in both the
862 !alpha-alpha and the beta-beta quadrants
863 quadrants = [.true., .false., .true.]
864 CALL ri_int_product(work_mat, lhs_int, rhs_int, quadrants, qs_env, &
865 eps_filter=xas_tdp_control%eps_filter, mo_transpose=.true.)
866
867 ELSE IF (do_uks) THEN
868 !In UKS, the donor MOs are possibly different for each spin => start with the
869 !alpha-alpha product and the perform the beta-beta product separately
870 quadrants = [.true., .false., .false.]
871 CALL ri_int_product(work_mat, lhs_int(1:ndo_mo), rhs_int(1:ndo_mo), quadrants, &
872 qs_env, eps_filter=xas_tdp_control%eps_filter, mo_transpose=.true.)
873
874 quadrants = [.false., .false., .true.]
875 CALL ri_int_product(work_mat, lhs_int(ndo_mo + 1:2*ndo_mo), rhs_int(ndo_mo + 1:2*ndo_mo), &
876 quadrants, qs_env, eps_filter=xas_tdp_control%eps_filter, mo_transpose=.true.)
877 ELSE
878 !In the restricted closed-shell case, only have one spin and a single qudarant
879 quadrants = [.true., .false., .false.]
880 CALL ri_int_product(work_mat, lhs_int, rhs_int, quadrants, qs_env, &
881 eps_filter=xas_tdp_control%eps_filter, mo_transpose=.true.)
882 END IF
883 CALL dbcsr_finalize(work_mat)
884
885 !Create the symmetric kernel matrix and redistribute work_mat into it
886 CALL dbcsr_create(offdiag_ex_ker, name="OFFDIAG EX KERNEL", matrix_type=dbcsr_type_symmetric, &
887 dist=dist, row_blk_size=blk_size, col_blk_size=blk_size)
888 CALL dbcsr_complete_redistribute(work_mat, offdiag_ex_ker)
889
890 !clean-up
891 CALL dbcsr_release(work_mat)
892 CALL dbcsr_deallocate_matrix_set(lhs_int)
893
894 END SUBROUTINE offdiag_ex_sc
895
896! **************************************************************************************************
897!> \brief Reserves the blocks in of a dbcsr matrix as needed for RI 3-center contraction (aI|P)
898!> \param matrices the matrices for which blocks are reserved
899!> \param ri_atom the index of the atom on which RI is done (= all coeffs of I are there, and P too)
900!> \param qs_env ...
901!> \note the end product are normal type matrices that are possibly slightly spraser as matrix_s
902! **************************************************************************************************
903 SUBROUTINE reserve_contraction_blocks(matrices, ri_atom, qs_env)
904
905 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrices
906 INTEGER, INTENT(IN) :: ri_atom
907 TYPE(qs_environment_type), POINTER :: qs_env
908
909 INTEGER :: i, iblk, jblk, max_nblks, nblks
910 INTEGER, ALLOCATABLE, DIMENSION(:) :: reserve_cols, reserve_rows
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)
917
918 ! Initialization
919 CALL get_qs_env(qs_env, matrix_s=matrix_s)
920 CALL dbcsr_get_info(matrices(1)%matrix, distribution=dist)
921
922 ! Need to redistribute matrix_s in the distribution of matrices
923 CALL dbcsr_create(work, template=matrix_s(1)%matrix, dist=dist)
924 CALL dbcsr_complete_redistribute(matrix_s(1)%matrix, work)
925
926 ! Need to desymmetrize matrix as as a template
927 CALL dbcsr_desymmetrize(work, template)
928
929 ! Allocate space for block indicies to reserve.
930 max_nblks = dbcsr_get_num_blocks(template)
931 ALLOCATE (reserve_rows(max_nblks), reserve_cols(max_nblks))
932
933 ! Loop over matrix_s as need a,b to overlap
934 nblks = 0
935 CALL dbcsr_iterator_start(iter, template)
936 DO WHILE (dbcsr_iterator_blocks_left(iter))
937 CALL dbcsr_iterator_next_block(iter, row=iblk, column=jblk)
938 !only have a,b pair if one of them is the ri_atom
939 IF (iblk == ri_atom .OR. jblk == ri_atom) THEN
940 nblks = nblks + 1
941 reserve_rows(nblks) = iblk
942 reserve_cols(nblks) = jblk
943 END IF
944 END DO
945 CALL dbcsr_iterator_stop(iter)
946
947 DO i = 1, SIZE(matrices)
948 CALL dbcsr_reserve_blocks(matrices(i)%matrix, rows=reserve_rows(1:nblks), cols=reserve_cols(1:nblks))
949 END DO
950
951 ! Clean-up
952 CALL dbcsr_release(template)
953 CALL dbcsr_release(work)
954
955 END SUBROUTINE reserve_contraction_blocks
956
957! **************************************************************************************************
958!> \brief Contract the ri 3-center integrals stored in a tensor with repect to the donor MOs coeffs,
959!> for a given excited atom k => (aI|k) = sum_b c_Ib (ab|k)
960!> \param contr_int the contracted integrals as array of dbcsr matrices
961!> \param op_type for which operator type we contract (COULOMB or EXCHANGE)
962!> \param donor_state ...
963!> \param xas_tdp_env ...
964!> \param xas_tdp_control ...
965!> \param qs_env ...
966!> \note In the output matrices, (aI_b|k) is stored at block a,b where I_b is the partial
967!> contraction that only includes coeffs from atom b. Note that the contracted matrix is
968!> not symmetric. To get the fully contracted matrix over b, one need to add the block
969!> columns of (aI_b|k) (get an array of size nao*nsgfp). This step is unnessary in our case
970!> because we assume locality of donor state, and only one column od (aI_b|k) is pouplated
971! **************************************************************************************************
972 SUBROUTINE contract2_ao_to_domo(contr_int, op_type, donor_state, xas_tdp_env, xas_tdp_control, qs_env)
973
974 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: contr_int
975 CHARACTER(len=*), INTENT(IN) :: op_type
976 TYPE(donor_state_type), POINTER :: donor_state
977 TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
978 TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
979 TYPE(qs_environment_type), POINTER :: qs_env
980
981 CHARACTER(len=*), PARAMETER :: routinen = 'contract2_AO_to_doMO'
982
983 INTEGER :: handle, i, imo, ispin, katom, kkind, &
984 natom, ndo_mo, ndo_so, nkind, nspins
985 INTEGER, DIMENSION(:), POINTER :: ri_blk_size, std_blk_size
986 LOGICAL :: do_uks
987 REAL(dp), DIMENSION(:, :), POINTER :: coeffs
988 TYPE(dbcsr_distribution_type) :: opt_dbcsr_dist
989 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrices, matrix_s
990 TYPE(dbcsr_type), POINTER :: ai_p, p_ib, work
991 TYPE(dbt_type), POINTER :: pq_x
992 TYPE(distribution_2d_type), POINTER :: opt_dist2d
993 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: ri_basis
994 TYPE(mp_para_env_type), POINTER :: para_env
995 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
996 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
997
998 NULLIFY (matrix_s, std_blk_size, ri_blk_size, qs_kind_set, ri_basis, pq_x)
999 NULLIFY (ai_p, p_ib, work, matrices, coeffs, opt_dist2d, particle_set)
1000
1001 CALL timeset(routinen, handle)
1002
1003! Initialization
1004 CALL get_qs_env(qs_env, natom=natom, matrix_s=matrix_s, qs_kind_set=qs_kind_set, para_env=para_env)
1005 ndo_mo = donor_state%ndo_mo
1006 kkind = donor_state%kind_index
1007 katom = donor_state%at_index
1008 !by default contract for Coulomb
1009 pq_x => xas_tdp_env%ri_3c_coul
1010 opt_dist2d => xas_tdp_env%opt_dist2d_coul
1011 IF (op_type == "EXCHANGE") THEN
1012 cpassert(ASSOCIATED(xas_tdp_env%ri_3c_ex))
1013 pq_x => xas_tdp_env%ri_3c_ex
1014 opt_dist2d => xas_tdp_env%opt_dist2d_ex
1015 END IF
1016 do_uks = xas_tdp_control%do_uks
1017 nspins = 1; IF (do_uks) nspins = 2
1018 ndo_so = nspins*ndo_mo
1019
1020! contracted integrals block sizes
1021 CALL dbcsr_get_info(matrix_s(1)%matrix, col_blk_size=std_blk_size)
1022 ! getting the block dimensions for the RI basis
1023 CALL get_qs_env(qs_env, particle_set=particle_set, nkind=nkind)
1024 ALLOCATE (ri_basis(nkind), ri_blk_size(natom))
1025 CALL basis_set_list_setup(ri_basis, "RI_XAS", qs_kind_set)
1026 CALL get_particle_set(particle_set, qs_kind_set, nsgf=ri_blk_size, basis=ri_basis)
1027
1028! Create work matrices. Everything that goes into a 3c routine must be compatible with the optimal dist_2d
1029 CALL cp_dbcsr_dist2d_to_dist(opt_dist2d, opt_dbcsr_dist)
1030
1031 ALLOCATE (ai_p, p_ib, work, matrices(2))
1032 CALL dbcsr_create(ai_p, dist=opt_dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, name="(aI|P)", &
1033 row_blk_size=std_blk_size, col_blk_size=ri_blk_size)
1034
1035 CALL dbcsr_create(p_ib, dist=opt_dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, name="(P|Ib)", &
1036 row_blk_size=ri_blk_size, col_blk_size=std_blk_size)
1037
1038 !reserve the blocks (needed for 3c contraction routines)
1039 matrices(1)%matrix => ai_p; matrices(2)%matrix => p_ib
1040 CALL reserve_contraction_blocks(matrices, katom, qs_env)
1041 DEALLOCATE (matrices)
1042
1043 ! Create the contracted integral matrices
1044 ALLOCATE (contr_int(ndo_so))
1045 DO i = 1, ndo_so
1046 ALLOCATE (contr_int(i)%matrix)
1047 CALL dbcsr_create(matrix=contr_int(i)%matrix, template=matrix_s(1)%matrix, &
1048 matrix_type=dbcsr_type_no_symmetry, row_blk_size=std_blk_size, &
1049 col_blk_size=ri_blk_size)
1050 END DO
1051
1052 ! Only take the coeffs for atom on which MOs I,J are localized
1053 coeffs => donor_state%contract_coeffs
1054
1055 DO ispin = 1, nspins
1056
1057! Loop over the donor MOs and contract
1058 DO imo = 1, ndo_mo
1059
1060 ! do the contraction
1061 CALL dbcsr_set(ai_p, 0.0_dp); CALL dbcsr_set(p_ib, 0.0_dp)
1062 CALL contract2_ao_to_domo_low(pq_x, coeffs(:, (ispin - 1)*ndo_mo + imo), ai_p, p_ib, katom)
1063
1064 ! Get the full (aI|P) contracted integrals
1065 CALL dbcsr_transposed(work, p_ib)
1066 CALL dbcsr_add(work, ai_p, 1.0_dp, 1.0_dp)
1067 CALL dbcsr_complete_redistribute(work, contr_int((ispin - 1)*ndo_mo + imo)%matrix)
1068 CALL dbcsr_filter(contr_int((ispin - 1)*ndo_mo + imo)%matrix, 1.0e-16_dp)
1069
1070 CALL dbcsr_release(work)
1071 END DO !imo
1072 END DO !ispin
1073
1074! Clean-up
1075 CALL dbcsr_release(ai_p)
1076 CALL dbcsr_release(p_ib)
1077 CALL dbcsr_distribution_release(opt_dbcsr_dist)
1078 DEALLOCATE (ri_blk_size, ai_p, p_ib, work, ri_basis)
1079
1080 CALL timestop(handle)
1081
1082 END SUBROUTINE contract2_ao_to_domo
1083
1084! **************************************************************************************************
1085!> \brief Contraction of the 3-center integrals (ab|Q) over the RI basis elements Q to get donor MOS
1086!> => (ab|IJ) = sum_X (ab|Q) coeffs_Q
1087!> \param ab_Q the tensor holding the integrals
1088!> \param vec the contraction coefficients
1089!> \param mat_abIJ the matrix holding the (ab|IJ) integrals (blocks must be reserved)
1090!> \param atom_k the atom for which we contract, i.e. we only take RI basis Q centered on atom_k
1091!> \note By construction, distribution of tensor and matrix match, also for OMP threads
1092! **************************************************************************************************
1093 SUBROUTINE contract3_ri_to_domos(ab_Q, vec, mat_abIJ, atom_k)
1094
1095 TYPE(dbt_type) :: ab_q
1096 REAL(dp), DIMENSION(:), INTENT(IN) :: vec
1097 TYPE(dbcsr_type) :: mat_abij
1098 INTEGER, INTENT(IN) :: atom_k
1099
1100 CHARACTER(len=*), PARAMETER :: routinen = 'contract3_RI_to_doMOs'
1101
1102 INTEGER :: handle, i, iatom, ind(3), j, jatom, katom
1103 LOGICAL :: found, t_found
1104 REAL(dp) :: prefac
1105 REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: iabc
1106 REAL(kind=dp), DIMENSION(:, :), POINTER :: pblock
1107 TYPE(dbcsr_type) :: work
1108 TYPE(dbt_iterator_type) :: iter
1109
1110 NULLIFY (pblock)
1111
1112 CALL timeset(routinen, handle)
1113
1114!$OMP PARALLEL DEFAULT(NONE) &
1115!$OMP SHARED(ab_Q,vec,mat_abIJ,atom_k) &
1116!$OMP PRIVATE(iter,ind,iatom,jatom,katom,prefac,iabc,t_found,found,pblock,i,j)
1117 CALL dbt_iterator_start(iter, ab_q)
1118 DO WHILE (dbt_iterator_blocks_left(iter))
1119 CALL dbt_iterator_next_block(iter, ind)
1120
1121 iatom = ind(1)
1122 jatom = ind(2)
1123 katom = ind(3)
1124
1125 IF (.NOT. atom_k == katom) cycle
1126
1127 prefac = 1.0_dp
1128 IF (iatom == jatom) prefac = 0.5_dp
1129
1130 CALL dbt_get_block(ab_q, ind, iabc, t_found)
1131
1132 CALL dbcsr_get_block_p(mat_abij, iatom, jatom, pblock, found)
1133 IF ((.NOT. found) .OR. (.NOT. t_found)) cycle
1134
1135 DO i = 1, SIZE(pblock, 1)
1136 DO j = 1, SIZE(pblock, 2)
1137!$OMP ATOMIC
1138 pblock(i, j) = pblock(i, j) + prefac*dot_product(vec(:), iabc(i, j, :))
1139 END DO
1140 END DO
1141
1142 DEALLOCATE (iabc)
1143 END DO !iter
1144 CALL dbt_iterator_stop(iter)
1145!$OMP END PARALLEL
1146
1147 !matrix only half filled => need to add its transpose
1148 CALL dbcsr_create(work, template=mat_abij)
1149 CALL dbcsr_transposed(work, mat_abij)
1150 CALL dbcsr_add(mat_abij, work, 1.0_dp, 1.0_dp)
1151 CALL dbcsr_release(work)
1152
1153 CALL timestop(handle)
1154
1155 END SUBROUTINE contract3_ri_to_domos
1156
1157! **************************************************************************************************
1158!> \brief Contraction of the 3-center integrals over index 1 and 2, for a given atom_k. The results
1159!> are stored in two matrices, such that (a,b are block indices):
1160!> mat_aIb(ab) = mat_aIb(ab) + sum j_b (i_aj_b|k)*v(j_b) and
1161!> mat_bIa(ba) = mat_bIa(ba) + sum i_a (i_aj_b|k)*v(i_a)
1162!> The block size of the columns of mat_aIb and the rows of mat_bIa are the size of k (RI)
1163!> \param ab_Q the tensor containing the 3-center integrals
1164!> \param vec the contraction coefficients
1165!> \param mat_aIb normal type dbcsr matrix
1166!> \param mat_bIa normal type dbcsr matrix
1167!> \param atom_k the atom for which we contract
1168!> It is assumed that the contraction coefficients for MO I are all on atom_k
1169!> We do the classic thing when we fill half the matrix and add its transposed to get the full
1170!> one, but here, the matrix is not symmetric, hence we explicitely have 2 input matrices
1171!> The distribution of the integrals and the normal dbcsr matrix are compatible out of the box
1172! **************************************************************************************************
1173 SUBROUTINE contract2_ao_to_domo_low(ab_Q, vec, mat_aIb, mat_bIa, atom_k)
1174
1175 TYPE(dbt_type) :: ab_q
1176 REAL(dp), DIMENSION(:), INTENT(IN) :: vec
1177 TYPE(dbcsr_type), INTENT(INOUT) :: mat_aib, mat_bia
1178 INTEGER, INTENT(IN) :: atom_k
1179
1180 CHARACTER(LEN=*), PARAMETER :: routinen = 'contract2_AO_to_doMO_low'
1181
1182 INTEGER :: handle, i, iatom, ind(3), j, jatom, &
1183 katom, s1, s2
1184 INTEGER, DIMENSION(:), POINTER :: atom_blk_size
1185 LOGICAL :: found, t_found
1186 REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: iabc
1187 REAL(dp), DIMENSION(:, :), POINTER :: pblock
1188 TYPE(dbt_iterator_type) :: iter
1189
1190 NULLIFY (atom_blk_size, pblock)
1191
1192 CALL timeset(routinen, handle)
1193
1194 CALL dbcsr_get_info(mat_aib, row_blk_size=atom_blk_size)
1195
1196!$OMP PARALLEL DEFAULT(NONE) &
1197!$OMP SHARED(ab_Q,vec,mat_aIb,mat_bIa,atom_k,atom_blk_size) &
1198!$OMP PRIVATE(iter,ind,iatom,jatom,katom,iabc,t_found,found,s1,s2,j,i,pblock)
1199 CALL dbt_iterator_start(iter, ab_q)
1200 DO WHILE (dbt_iterator_blocks_left(iter))
1201 CALL dbt_iterator_next_block(iter, ind)
1202
1203 iatom = ind(1)
1204 jatom = ind(2)
1205 katom = ind(3)
1206
1207 IF (atom_k /= katom) cycle
1208
1209 CALL dbt_get_block(ab_q, ind, iabc, t_found)
1210 IF (.NOT. t_found) cycle
1211
1212 ! Deal with mat_aIb
1213 IF (jatom == atom_k) THEN
1214 s1 = atom_blk_size(iatom)
1215 s2 = SIZE(iabc, 3)
1216
1217 CALL dbcsr_get_block_p(matrix=mat_aib, row=iatom, col=jatom, block=pblock, found=found)
1218
1219 IF (found) THEN
1220 DO i = 1, s1
1221 DO j = 1, s2
1222!$OMP ATOMIC
1223 pblock(i, j) = pblock(i, j) + dot_product(vec, iabc(i, :, j))
1224 END DO
1225 END DO
1226 END IF
1227 END IF ! jatom == atom_k
1228
1229 ! Deal with mat_bIa, keep block diagonal empty
1230 IF (iatom == jatom) cycle
1231 IF (iatom == atom_k) THEN
1232 s1 = SIZE(iabc, 3)
1233 s2 = atom_blk_size(jatom)
1234
1235 CALL dbcsr_get_block_p(matrix=mat_bia, row=iatom, col=jatom, block=pblock, found=found)
1236
1237 IF (found) THEN
1238 DO i = 1, s1
1239 DO j = 1, s2
1240!$OMP ATOMIC
1241 pblock(i, j) = pblock(i, j) + dot_product(vec, iabc(:, j, i))
1242 END DO
1243 END DO
1244 END IF
1245 END IF !iatom== atom_k
1246
1247 DEALLOCATE (iabc)
1248 END DO !iter
1249 CALL dbt_iterator_stop(iter)
1250!$OMP END PARALLEL
1251
1252 CALL timestop(handle)
1253
1254 END SUBROUTINE contract2_ao_to_domo_low
1255
1256! **************************************************************************************************
1257!> \brief Multiply all the blocks of a contracted RI integral (aI|P) by a matrix of type (P|...|Q)
1258!> \param contr_int the integral array
1259!> \param PQ the smaller matrix to multiply all blocks
1260!> \note It is assumed that all non-zero blocks have the same number of columns. Can pass partial
1261!> arrays, e.g. contr_int(1:3)
1262! **************************************************************************************************
1263 SUBROUTINE ri_all_blocks_mm(contr_int, PQ)
1264
1265 TYPE(dbcsr_p_type), DIMENSION(:) :: contr_int
1266 REAL(dp), DIMENSION(:, :), INTENT(IN) :: pq
1267
1268 INTEGER :: iblk, imo, jblk, ndo_mo, s1, s2
1269 LOGICAL :: found
1270 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: work
1271 REAL(dp), DIMENSION(:, :), POINTER :: pblock
1272 TYPE(dbcsr_iterator_type) :: iter
1273
1274 NULLIFY (pblock)
1275
1276 ndo_mo = SIZE(contr_int)
1277
1278 DO imo = 1, ndo_mo
1279 CALL dbcsr_iterator_start(iter, contr_int(imo)%matrix)
1280 DO WHILE (dbcsr_iterator_blocks_left(iter))
1281
1282 CALL dbcsr_iterator_next_block(iter, row=iblk, column=jblk)
1283 CALL dbcsr_get_block_p(contr_int(imo)%matrix, iblk, jblk, pblock, found)
1284
1285 IF (found) THEN
1286 s1 = SIZE(pblock, 1)
1287 s2 = SIZE(pblock, 2)
1288 ALLOCATE (work(s1, s2))
1289 CALL dgemm('N', 'N', s1, s2, s2, 1.0_dp, pblock, s1, pq, s2, 0.0_dp, work, s1)
1290 CALL dcopy(s1*s2, work, 1, pblock, 1)
1291 DEALLOCATE (work)
1292 END IF
1293
1294 END DO ! dbcsr iterator
1295 CALL dbcsr_iterator_stop(iter)
1296 END DO !imo
1297
1298 END SUBROUTINE ri_all_blocks_mm
1299
1300! **************************************************************************************************
1301!> \brief Copies an (partial) array of contracted RI integrals into anoter one
1302!> \param new_int where the copy is stored
1303!> \param ref_int what is copied
1304!> \note Allocate the matrices of new_int if not done already
1305! **************************************************************************************************
1306 SUBROUTINE copy_ri_contr_int(new_int, ref_int)
1307
1308 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: new_int
1309 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: ref_int
1310
1311 INTEGER :: iso, ndo_so
1312
1313 cpassert(SIZE(new_int) == SIZE(ref_int))
1314 ndo_so = SIZE(ref_int)
1315
1316 DO iso = 1, ndo_so
1317 IF (.NOT. ASSOCIATED(new_int(iso)%matrix)) ALLOCATE (new_int(iso)%matrix)
1318 CALL dbcsr_copy(new_int(iso)%matrix, ref_int(iso)%matrix)
1319 END DO
1320
1321 END SUBROUTINE copy_ri_contr_int
1322
1323! **************************************************************************************************
1324!> \brief Takes the product of contracted integrals and put them in a kernel matrix
1325!> \param kernel the matrix where the products are stored
1326!> \param lhs_int the left-hand side contracted integrals
1327!> \param rhs_int the right-hand side contracted integrals
1328!> \param quadrants on which quadrant(s) on the kernel matrix the product is stored
1329!> \param qs_env ...
1330!> \param eps_filter filter for dbcsr matrix multiplication
1331!> \param mo_transpose whether the MO blocks should be transpose, i.e. (aI|Jb) => (aJ|Ib)
1332!> \note It is assumed that the kerenl matrix is NOT symmetric
1333!> There are three quadrants, corresponding to 1: the upper-left (diagonal), 2: the
1334!> upper-right (off-diagonal) and 3: the lower-right (diagonal).
1335!> Need to finalize the kernel matrix after calling this routine (possibly multiple times)
1336! **************************************************************************************************
1337 SUBROUTINE ri_int_product(kernel, lhs_int, rhs_int, quadrants, qs_env, eps_filter, mo_transpose)
1338
1339 TYPE(dbcsr_type), INTENT(INOUT) :: kernel
1340 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: lhs_int, rhs_int
1341 LOGICAL, DIMENSION(3), INTENT(IN) :: quadrants
1342 TYPE(qs_environment_type), POINTER :: qs_env
1343 REAL(dp), INTENT(IN), OPTIONAL :: eps_filter
1344 LOGICAL, INTENT(IN), OPTIONAL :: mo_transpose
1345
1346 INTEGER :: i, iblk, iso, j, jblk, jso, nblk, ndo_so
1347 LOGICAL :: found, my_mt
1348 REAL(dp), DIMENSION(:, :), POINTER :: pblock
1349 TYPE(dbcsr_iterator_type) :: iter
1350 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1351 TYPE(dbcsr_type) :: prod
1352
1353 NULLIFY (matrix_s, pblock)
1354
1355! Initialization
1356 cpassert(SIZE(lhs_int) == SIZE(rhs_int))
1357 cpassert(any(quadrants))
1358 ndo_so = SIZE(lhs_int)
1359 CALL get_qs_env(qs_env, matrix_s=matrix_s, natom=nblk)
1360 CALL dbcsr_create(prod, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1361 my_mt = .false.
1362 IF (PRESENT(mo_transpose)) my_mt = mo_transpose
1363
1364 ! The kernel matrix is symmetric (even if normal type) => only fill upper half on diagonal
1365 ! quadrants, but the whole thing on upper-right quadrant
1366 DO iso = 1, ndo_so
1367 DO jso = 1, ndo_so
1368
1369 ! If on-diagonal quadrants only, can skip jso < iso
1370 IF (.NOT. quadrants(2) .AND. jso < iso) cycle
1371
1372 i = iso; j = jso
1373 IF (my_mt) THEN
1374 i = jso; j = iso
1375 END IF
1376
1377 ! Take the product lhs*rhs^T
1378 CALL dbcsr_multiply('N', 'T', 1.0_dp, lhs_int(i)%matrix, rhs_int(j)%matrix, &
1379 0.0_dp, prod, filter_eps=eps_filter)
1380
1381 ! Loop over blocks of prod and fill kernel matrix => ok cuz same (but replicated) dist
1382 CALL dbcsr_iterator_start(iter, prod)
1383 DO WHILE (dbcsr_iterator_blocks_left(iter))
1384
1385 CALL dbcsr_iterator_next_block(iter, row=iblk, column=jblk)
1386 IF ((iso == jso .AND. jblk < iblk) .AND. .NOT. quadrants(2)) cycle
1387
1388 CALL dbcsr_get_block_p(prod, iblk, jblk, pblock, found)
1389
1390 IF (found) THEN
1391
1392 ! Case study on quadrant
1393 !upper-left
1394 IF (quadrants(1)) THEN
1395 CALL dbcsr_put_block(kernel, (iso - 1)*nblk + iblk, (jso - 1)*nblk + jblk, pblock)
1396 END IF
1397
1398 !upper-right
1399 IF (quadrants(2)) THEN
1400 CALL dbcsr_put_block(kernel, (iso - 1)*nblk + iblk, (ndo_so + jso - 1)*nblk + jblk, pblock)
1401 END IF
1402
1403 !lower-right
1404 IF (quadrants(3)) THEN
1405 CALL dbcsr_put_block(kernel, (ndo_so + iso - 1)*nblk + iblk, (ndo_so + jso - 1)*nblk + jblk, pblock)
1406 END IF
1407
1408 END IF
1409
1410 END DO ! dbcsr iterator
1411 CALL dbcsr_iterator_stop(iter)
1412
1413 END DO !jso
1414 END DO !iso
1415
1416! Clean-up
1417 CALL dbcsr_release(prod)
1418
1419 END SUBROUTINE ri_int_product
1420
1421END 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:225
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:85
basis_set_types
Definition basis_set_types.F:15
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_transposed
subroutine, public dbcsr_transposed(transposed, normal, shallow_data_copy, transpose_distribution, use_distribution)
...
Definition cp_dbcsr_api.F:1242
cp_dbcsr_api::dbcsr_distribution_release
subroutine, public dbcsr_distribution_release(dist)
...
Definition cp_dbcsr_api.F:643
cp_dbcsr_api::dbcsr_distribution_new
subroutine, public dbcsr_distribution_new(dist, template, group, pgrid, row_dist, col_dist, reuse_arrays)
...
Definition cp_dbcsr_api.F:616
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:944
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1050
cp_dbcsr_api::dbcsr_desymmetrize
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
Definition cp_dbcsr_api.F:525
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:693
cp_dbcsr_api::dbcsr_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1096
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:808
cp_dbcsr_api::dbcsr_reserve_blocks
subroutine, public dbcsr_reserve_blocks(matrix, rows, cols)
...
Definition cp_dbcsr_api.F:1172
cp_dbcsr_api::dbcsr_get_stored_coordinates
subroutine, public dbcsr_get_stored_coordinates(matrix, row, column, processor)
...
Definition cp_dbcsr_api.F:911
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset, transposed)
...
Definition cp_dbcsr_api.F:971
cp_dbcsr_api::dbcsr_filter
subroutine, public dbcsr_filter(matrix, eps)
...
Definition cp_dbcsr_api.F:658
cp_dbcsr_api::dbcsr_finalize
subroutine, public dbcsr_finalize(matrix)
...
Definition cp_dbcsr_api.F:673
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:1012
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1142
cp_dbcsr_api::dbcsr_complete_redistribute
subroutine, public dbcsr_complete_redistribute(matrix, redist)
...
Definition cp_dbcsr_api.F:320
cp_dbcsr_api::dbcsr_get_num_blocks
integer function, public dbcsr_get_num_blocks(matrix)
...
Definition cp_dbcsr_api.F:877
cp_dbcsr_api::dbcsr_put_block
subroutine, public dbcsr_put_block(matrix, row, col, block, summation)
...
Definition cp_dbcsr_api.F:1125
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:254
cp_dbcsr_api::dbcsr_distribution_get
subroutine, public dbcsr_distribution_get(dist, row_dist, col_dist, nrows, ncols, has_threads, group, mynode, numnodes, nprows, npcols, myprow, mypcol, pgrid, subgroups_defined, prow_group, pcol_group)
...
Definition cp_dbcsr_api.F:572
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_dist2d_to_dist
subroutine, public cp_dbcsr_dist2d_to_dist(dist2d, dist)
Creates a DBCSR distribution from a distribution_2d.
Definition cp_dbcsr_operations.F:415
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, ncgf)
Get the components of a particle set.
Definition particle_methods.F:121
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:534
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:973
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:596
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:904
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:1264
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:1455
basis_set_types::gto_basis_set_p_type
Definition basis_set_types.F:94
cp_dbcsr_api::dbcsr_distribution_type
Definition cp_dbcsr_api.F:182
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:188
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
distribution_2d_types::distribution_2d_type
distributes pairs on a 2d grid of processors
Definition distribution_2d_types.F:62
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:750
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
xas_tdp_types::donor_state_type
Type containing informations about a single donor state.
Definition xas_tdp_types.F:273
xas_tdp_types::xas_tdp_control_type
Type containing control information for TDP XAS calculations.
Definition xas_tdp_types.F:102
xas_tdp_types::xas_tdp_env_type
Type containing informations such as inputs and results for TDP XAS calculations.
Definition xas_tdp_types.F:195