Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief This module deals with all the integrals done on local atomic grids in xas_tdp. This is
10 : !> mostly used to compute the xc kernel matrix elements wrt two RI basis elements (centered
11 : !> on the same excited atom) <P|fxc(r)|Q>, where the kernel fxc is purely a function of the
12 : !> ground state density and r. This is also used to compute the SOC matrix elements in the
13 : !> orbital basis
14 : ! **************************************************************************************************
15 : MODULE xas_tdp_atom
16 : USE ai_contraction_sphi, ONLY: ab_contract
17 : USE atom_operators, ONLY: calculate_model_potential
18 : USE basis_set_types, ONLY: get_gto_basis_set,&
19 : gto_basis_set_p_type,&
20 : gto_basis_set_type
21 : USE cell_types, ONLY: cell_type,&
22 : pbc
23 : USE cp_array_utils, ONLY: cp_1d_i_p_type,&
24 : cp_1d_r_p_type,&
25 : cp_2d_r_p_type,&
26 : cp_3d_r_p_type
27 : USE cp_blacs_env, ONLY: cp_blacs_env_type
28 : USE cp_control_types, ONLY: dft_control_type,&
29 : qs_control_type
30 : USE cp_dbcsr_api, ONLY: &
31 : dbcsr_copy, dbcsr_create, dbcsr_distribution_get, dbcsr_distribution_new, &
32 : dbcsr_distribution_release, dbcsr_distribution_type, dbcsr_filter, dbcsr_finalize, &
33 : dbcsr_get_block_p, dbcsr_get_stored_coordinates, dbcsr_iterator_blocks_left, &
34 : dbcsr_iterator_next_block, dbcsr_iterator_readonly_start, dbcsr_iterator_start, &
35 : dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, dbcsr_put_block, dbcsr_release, &
36 : dbcsr_replicate_all, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, &
37 : dbcsr_type_no_symmetry, dbcsr_type_symmetric
38 : USE cp_dbcsr_cholesky, ONLY: cp_dbcsr_cholesky_decompose,&
39 : cp_dbcsr_cholesky_invert
40 : USE cp_dbcsr_operations, ONLY: dbcsr_deallocate_matrix_set
41 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
42 : USE dbt_api, ONLY: dbt_destroy,&
43 : dbt_get_block,&
44 : dbt_iterator_blocks_left,&
45 : dbt_iterator_next_block,&
46 : dbt_iterator_start,&
47 : dbt_iterator_stop,&
48 : dbt_iterator_type,&
49 : dbt_type
50 : USE input_constants, ONLY: do_potential_id
51 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
52 : section_vals_type
53 : USE kinds, ONLY: default_string_length,&
54 : dp
55 : USE lebedev, ONLY: deallocate_lebedev_grids,&
56 : get_number_of_lebedev_grid,&
57 : init_lebedev_grids,&
58 : lebedev_grid
59 : USE libint_2c_3c, ONLY: libint_potential_type
60 : USE mathlib, ONLY: get_diag,&
61 : invmat_symm
62 : USE memory_utilities, ONLY: reallocate
63 : USE message_passing, ONLY: mp_comm_type,&
64 : mp_para_env_type,&
65 : mp_request_type,&
66 : mp_waitall
67 : USE orbital_pointers, ONLY: indco,&
68 : indso,&
69 : nco,&
70 : ncoset,&
71 : nso,&
72 : nsoset
73 : USE orbital_transformation_matrices, ONLY: orbtramat
74 : USE particle_methods, ONLY: get_particle_set
75 : USE particle_types, ONLY: particle_type
76 : USE physcon, ONLY: c_light_au
77 : USE qs_environment_types, ONLY: get_qs_env,&
78 : qs_environment_type
79 : USE qs_grid_atom, ONLY: allocate_grid_atom,&
80 : create_grid_atom,&
81 : grid_atom_type
82 : USE qs_harmonics_atom, ONLY: allocate_harmonics_atom,&
83 : create_harmonics_atom,&
84 : get_maxl_CG,&
85 : get_none0_cg_list,&
86 : harmonics_atom_type
87 : USE qs_integral_utils, ONLY: basis_set_list_setup
88 : USE qs_kind_types, ONLY: get_qs_kind,&
89 : get_qs_kind_set,&
90 : qs_kind_type
91 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type,&
92 : release_neighbor_list_sets
93 : USE qs_overlap, ONLY: build_overlap_matrix_simple
94 : USE qs_rho_types, ONLY: qs_rho_get,&
95 : qs_rho_type
96 : USE qs_tddfpt2_soc_types, ONLY: soc_atom_env_type
97 : USE spherical_harmonics, ONLY: clebsch_gordon,&
98 : clebsch_gordon_deallocate,&
99 : clebsch_gordon_init
100 : USE util, ONLY: get_limit,&
101 : sort_unique
102 : USE xas_tdp_integrals, ONLY: build_xas_tdp_3c_nl,&
103 : build_xas_tdp_ovlp_nl,&
104 : create_pqX_tensor,&
105 : fill_pqx_tensor
106 : USE xas_tdp_types, ONLY: batch_info_type,&
107 : get_proc_batch_sizes,&
108 : release_batch_info,&
109 : xas_atom_env_type,&
110 : xas_tdp_control_type,&
111 : xas_tdp_env_type
112 : USE xc, ONLY: divide_by_norm_drho
113 : USE xc_atom, ONLY: xc_rho_set_atom_update
114 : USE xc_derivative_desc, ONLY: deriv_norm_drho,&
115 : deriv_norm_drhoa,&
116 : deriv_norm_drhob,&
117 : deriv_rhoa,&
118 : deriv_rhob
119 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
120 : xc_dset_create,&
121 : xc_dset_get_derivative,&
122 : xc_dset_release
123 : USE xc_derivative_types, ONLY: xc_derivative_get,&
124 : xc_derivative_p_type,&
125 : xc_derivative_type
126 : USE xc_derivatives, ONLY: xc_functionals_eval,&
127 : xc_functionals_get_needs
128 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
129 : USE xc_rho_set_types, ONLY: xc_rho_set_create,&
130 : xc_rho_set_release,&
131 : xc_rho_set_type
132 :
133 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
134 : #include "./base/base_uses.f90"
135 :
136 : IMPLICIT NONE
137 :
138 : PRIVATE
139 :
140 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xas_tdp_atom'
141 :
142 : PUBLIC :: init_xas_atom_env, &
143 : integrate_fxc_atoms, &
144 : integrate_soc_atoms, &
145 : calculate_density_coeffs, &
146 : compute_sphi_so, &
147 : truncate_radial_grid, &
148 : init_xas_atom_grid_harmo
149 :
150 : CONTAINS
151 :
152 : ! **************************************************************************************************
153 : !> \brief Initializes a xas_atom_env type given the qs_enxas_atom_env, qs_envv
154 : !> \param xas_atom_env the xas_atom_env to initialize
155 : !> \param xas_tdp_env ...
156 : !> \param xas_tdp_control ...
157 : !> \param qs_env ...
158 : !> \param ltddfpt ...
159 : ! **************************************************************************************************
160 62 : SUBROUTINE init_xas_atom_env(xas_atom_env, xas_tdp_env, xas_tdp_control, qs_env, ltddfpt)
161 :
162 : TYPE(xas_atom_env_type), POINTER :: xas_atom_env
163 : TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
164 : TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
165 : TYPE(qs_environment_type), POINTER :: qs_env
166 : LOGICAL, OPTIONAL :: ltddfpt
167 :
168 : CHARACTER(len=*), PARAMETER :: routineN = 'init_xas_atom_env'
169 :
170 : INTEGER :: handle, ikind, natom, nex_atoms, &
171 : nex_kinds, nkind, nspins
172 : LOGICAL :: do_xc
173 62 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
174 62 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
175 :
176 62 : NULLIFY (qs_kind_set, particle_set)
177 :
178 62 : CALL timeset(routineN, handle)
179 :
180 : ! Initializing the type
181 62 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, natom=natom, particle_set=particle_set)
182 :
183 62 : nkind = SIZE(qs_kind_set)
184 62 : nex_kinds = xas_tdp_env%nex_kinds
185 62 : nex_atoms = xas_tdp_env%nex_atoms
186 62 : do_xc = xas_tdp_control%do_xc
187 62 : IF (PRESENT(ltddfpt)) THEN
188 0 : IF (ltddfpt) do_xc = .FALSE.
189 : END IF
190 62 : nspins = 1; IF (xas_tdp_control%do_uks .OR. xas_tdp_control%do_roks) nspins = 2
191 :
192 62 : xas_atom_env%nspins = nspins
193 62 : xas_atom_env%ri_radius = xas_tdp_control%ri_radius
194 284 : ALLOCATE (xas_atom_env%grid_atom_set(nkind))
195 284 : ALLOCATE (xas_atom_env%harmonics_atom_set(nkind))
196 284 : ALLOCATE (xas_atom_env%ri_sphi_so(nkind))
197 284 : ALLOCATE (xas_atom_env%orb_sphi_so(nkind))
198 62 : IF (do_xc) THEN
199 240 : ALLOCATE (xas_atom_env%gr(nkind))
200 240 : ALLOCATE (xas_atom_env%ga(nkind))
201 240 : ALLOCATE (xas_atom_env%dgr1(nkind))
202 240 : ALLOCATE (xas_atom_env%dgr2(nkind))
203 240 : ALLOCATE (xas_atom_env%dga1(nkind))
204 240 : ALLOCATE (xas_atom_env%dga2(nkind))
205 : END IF
206 :
207 62 : xas_atom_env%excited_atoms => xas_tdp_env%ex_atom_indices
208 62 : xas_atom_env%excited_kinds => xas_tdp_env%ex_kind_indices
209 :
210 : ! Allocate and initialize the atomic grids and harmonics
211 62 : CALL init_xas_atom_grid_harmo(xas_atom_env, xas_tdp_control%grid_info, do_xc, qs_env)
212 :
213 : ! Compute the contraction coefficients for spherical orbitals
214 160 : DO ikind = 1, nkind
215 98 : NULLIFY (xas_atom_env%orb_sphi_so(ikind)%array, xas_atom_env%ri_sphi_so(ikind)%array)
216 98 : CALL compute_sphi_so(ikind, "ORB", xas_atom_env%orb_sphi_so(ikind)%array, qs_env)
217 198 : IF (ANY(xas_atom_env%excited_kinds == ikind)) THEN
218 68 : CALL compute_sphi_so(ikind, "RI_XAS", xas_atom_env%ri_sphi_so(ikind)%array, qs_env)
219 : END IF
220 : END DO !ikind
221 :
222 : ! Compute the coefficients to expand the density in the RI_XAS basis, if requested
223 62 : IF (do_xc .AND. (.NOT. xas_tdp_control%xps_only)) THEN
224 52 : CALL calculate_density_coeffs(xas_atom_env=xas_atom_env, qs_env=qs_env)
225 : END IF
226 :
227 62 : CALL timestop(handle)
228 :
229 62 : END SUBROUTINE init_xas_atom_env
230 :
231 : ! **************************************************************************************************
232 : !> \brief Initializes the atomic grids and harmonics for the RI atomic calculations
233 : !> \param xas_atom_env ...
234 : !> \param grid_info ...
235 : !> \param do_xc Whether the xc kernel will ne computed on the atomic grids. If not, the harmonics
236 : !> are built for the orbital basis for all kinds.
237 : !> \param qs_env ...
238 : !> \note Largely inspired by init_rho_atom subroutine
239 : ! **************************************************************************************************
240 62 : SUBROUTINE init_xas_atom_grid_harmo(xas_atom_env, grid_info, do_xc, qs_env)
241 :
242 : TYPE(xas_atom_env_type), POINTER :: xas_atom_env
243 : CHARACTER(len=default_string_length), &
244 : DIMENSION(:, :), POINTER :: grid_info
245 : LOGICAL, INTENT(IN) :: do_xc
246 : TYPE(qs_environment_type), POINTER :: qs_env
247 :
248 : CHARACTER(LEN=default_string_length) :: kind_name
249 : INTEGER :: igrid, ikind, il, iso, iso1, iso2, l1, l1l2, l2, la, lc1, lc2, lcleb, ll, llmax, &
250 : lp, m1, m2, max_s_harm, max_s_set, maxl, maxlgto, maxs, mm, mp, na, nr, quadrature, stat
251 : REAL(dp) :: kind_radius
252 62 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rga
253 62 : REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG
254 : TYPE(dft_control_type), POINTER :: dft_control
255 : TYPE(grid_atom_type), POINTER :: grid_atom
256 : TYPE(gto_basis_set_type), POINTER :: tmp_basis
257 : TYPE(harmonics_atom_type), POINTER :: harmonics
258 : TYPE(qs_control_type), POINTER :: qs_control
259 62 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
260 :
261 62 : NULLIFY (my_CG, qs_kind_set, tmp_basis, grid_atom, harmonics, qs_control, dft_control)
262 :
263 : ! Initialization of some integer for the CG coeff generation
264 62 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, dft_control=dft_control)
265 62 : IF (do_xc) THEN
266 52 : CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto, basis_type="RI_XAS")
267 : ELSE
268 10 : CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto, basis_type="ORB")
269 : END IF
270 62 : qs_control => dft_control%qs_control
271 :
272 : !maximum expansion
273 62 : llmax = 2*maxlgto
274 62 : max_s_harm = nsoset(llmax)
275 62 : max_s_set = nsoset(maxlgto)
276 62 : lcleb = llmax
277 :
278 : ! Allocate and compute the CG coeffs (copied from init_rho_atom)
279 62 : CALL clebsch_gordon_init(lcleb)
280 62 : CALL reallocate(my_CG, 1, max_s_set, 1, max_s_set, 1, max_s_harm)
281 :
282 186 : ALLOCATE (rga(lcleb, 2))
283 300 : DO lc1 = 0, maxlgto
284 1282 : DO iso1 = nsoset(lc1 - 1) + 1, nsoset(lc1)
285 982 : l1 = indso(1, iso1)
286 982 : m1 = indso(2, iso1)
287 5598 : DO lc2 = 0, maxlgto
288 26382 : DO iso2 = nsoset(lc2 - 1) + 1, nsoset(lc2)
289 21022 : l2 = indso(1, iso2)
290 21022 : m2 = indso(2, iso2)
291 21022 : CALL clebsch_gordon(l1, m1, l2, m2, rga)
292 21022 : IF (l1 + l2 > llmax) THEN
293 : l1l2 = llmax
294 : ELSE
295 : l1l2 = l1 + l2
296 : END IF
297 21022 : mp = m1 + m2
298 21022 : mm = m1 - m2
299 21022 : IF (m1*m2 < 0 .OR. (m1*m2 == 0 .AND. (m1 < 0 .OR. m2 < 0))) THEN
300 10020 : mp = -ABS(mp)
301 10020 : mm = -ABS(mm)
302 : ELSE
303 11002 : mp = ABS(mp)
304 11002 : mm = ABS(mm)
305 : END IF
306 101890 : DO lp = MOD(l1 + l2, 2), l1l2, 2
307 76490 : il = lp/2 + 1
308 76490 : IF (ABS(mp) <= lp) THEN
309 52134 : IF (mp >= 0) THEN
310 29474 : iso = nsoset(lp - 1) + lp + 1 + mp
311 : ELSE
312 22660 : iso = nsoset(lp - 1) + lp + 1 - ABS(mp)
313 : END IF
314 52134 : my_CG(iso1, iso2, iso) = rga(il, 1)
315 : END IF
316 97512 : IF (mp /= mm .AND. ABS(mm) <= lp) THEN
317 34660 : IF (mm >= 0) THEN
318 22352 : iso = nsoset(lp - 1) + lp + 1 + mm
319 : ELSE
320 12308 : iso = nsoset(lp - 1) + lp + 1 - ABS(mm)
321 : END IF
322 34660 : my_CG(iso1, iso2, iso) = rga(il, 2)
323 : END IF
324 : END DO
325 : END DO ! iso2
326 : END DO ! lc2
327 : END DO ! iso1
328 : END DO ! lc1
329 62 : DEALLOCATE (rga)
330 62 : CALL clebsch_gordon_deallocate()
331 :
332 : ! Create the Lebedev grids and compute the spherical harmonics
333 62 : CALL init_lebedev_grids()
334 62 : quadrature = qs_control%gapw_control%quadrature
335 :
336 160 : DO ikind = 1, SIZE(xas_atom_env%grid_atom_set)
337 :
338 : ! Allocate the grid and the harmonics for this kind
339 98 : NULLIFY (xas_atom_env%grid_atom_set(ikind)%grid_atom)
340 98 : NULLIFY (xas_atom_env%harmonics_atom_set(ikind)%harmonics_atom)
341 98 : CALL allocate_grid_atom(xas_atom_env%grid_atom_set(ikind)%grid_atom)
342 98 : CALL allocate_harmonics_atom(xas_atom_env%harmonics_atom_set(ikind)%harmonics_atom)
343 :
344 98 : NULLIFY (grid_atom, harmonics)
345 98 : grid_atom => xas_atom_env%grid_atom_set(ikind)%grid_atom
346 98 : harmonics => xas_atom_env%harmonics_atom_set(ikind)%harmonics_atom
347 :
348 : ! Initialize some integers
349 98 : CALL get_qs_kind(qs_kind_set(ikind), ngrid_rad=nr, ngrid_ang=na, name=kind_name)
350 :
351 : !take the grid dimension given as input, if none, take the GAPW ones above
352 296 : IF (ANY(grid_info == kind_name)) THEN
353 68 : DO igrid = 1, SIZE(grid_info, 1)
354 68 : IF (grid_info(igrid, 1) == kind_name) THEN
355 :
356 : !hack to convert string into integer
357 62 : READ (grid_info(igrid, 2), *, iostat=stat) na
358 62 : IF (stat /= 0) CPABORT("The 'na' value for the GRID keyword must be an integer")
359 62 : READ (grid_info(igrid, 3), *, iostat=stat) nr
360 62 : IF (stat /= 0) CPABORT("The 'nr' value for the GRID keyword must be an integer")
361 : EXIT
362 : END IF
363 : END DO
364 68 : ELSE IF (do_xc .AND. ANY(xas_atom_env%excited_kinds == ikind)) THEN
365 : !need good integration grids for the xc kernel, but taking the default GAPW grid
366 0 : CPWARN("The default (and possibly small) GAPW grid is being used for one excited KIND")
367 : END IF
368 :
369 98 : ll = get_number_of_lebedev_grid(n=na)
370 98 : na = lebedev_grid(ll)%n
371 98 : la = lebedev_grid(ll)%l
372 98 : grid_atom%ng_sphere = na
373 98 : grid_atom%nr = nr
374 :
375 : ! If this is an excited kind, create the harmonics with the RI_XAS basis, otherwise the ORB
376 136 : IF (ANY(xas_atom_env%excited_kinds == ikind) .AND. do_xc) THEN
377 58 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=tmp_basis, basis_type="RI_XAS")
378 : ELSE
379 40 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=tmp_basis, basis_type="ORB")
380 : END IF
381 98 : CALL get_gto_basis_set(gto_basis_set=tmp_basis, maxl=maxl, kind_radius=kind_radius)
382 :
383 98 : CALL create_grid_atom(grid_atom, nr, na, llmax, ll, quadrature)
384 98 : CALL truncate_radial_grid(grid_atom, kind_radius)
385 :
386 98 : maxs = nsoset(maxl)
387 : CALL create_harmonics_atom(harmonics, &
388 : my_CG, na, llmax, maxs, max_s_harm, ll, grid_atom%wa, &
389 98 : grid_atom%azi, grid_atom%pol)
390 356 : CALL get_maxl_CG(harmonics, tmp_basis, llmax, max_s_harm)
391 : END DO
392 :
393 62 : CALL deallocate_lebedev_grids()
394 62 : DEALLOCATE (my_CG)
395 :
396 62 : END SUBROUTINE init_xas_atom_grid_harmo
397 :
398 : ! **************************************************************************************************
399 : !> \brief Reduces the radial extension of an atomic grid such that it only covers a given radius
400 : !> \param grid_atom the atomic grid from which we truncate the radial part
401 : !> \param max_radius the maximal radial extension of the resulting grid
402 : !> \note Since the RI density used for <P|fxc|Q> is only of quality close to the atom, and the
403 : !> integrand only non-zero within the radius of the gaussian P,Q. One can reduce the grid
404 : !> extansion to the largest radius of the RI basis set
405 : ! **************************************************************************************************
406 106 : SUBROUTINE truncate_radial_grid(grid_atom, max_radius)
407 :
408 : TYPE(grid_atom_type), POINTER :: grid_atom
409 : REAL(dp), INTENT(IN) :: max_radius
410 :
411 : INTEGER :: first_ir, ir, llmax_p1, na, new_nr, nr
412 :
413 106 : nr = grid_atom%nr
414 106 : na = grid_atom%ng_sphere
415 106 : llmax_p1 = SIZE(grid_atom%rad2l, 2) - 1
416 :
417 : ! Find the index corresponding to the limiting radius (small ir => large radius)
418 2778 : DO ir = 1, nr
419 2778 : IF (grid_atom%rad(ir) < max_radius) THEN
420 : first_ir = ir
421 : EXIT
422 : END IF
423 : END DO
424 106 : new_nr = nr - first_ir + 1
425 :
426 : ! Reallcoate everything that depends on r
427 106 : grid_atom%nr = new_nr
428 :
429 21148 : grid_atom%rad(1:new_nr) = grid_atom%rad(first_ir:nr)
430 21148 : grid_atom%rad2(1:new_nr) = grid_atom%rad2(first_ir:nr)
431 21148 : grid_atom%wr(1:new_nr) = grid_atom%wr(first_ir:nr)
432 3512444 : grid_atom%weight(:, 1:new_nr) = grid_atom%weight(:, first_ir:nr)
433 171548 : grid_atom%rad2l(1:new_nr, :) = grid_atom%rad2l(first_ir:nr, :)
434 171548 : grid_atom%oorad2l(1:new_nr, :) = grid_atom%oorad2l(first_ir:nr, :)
435 :
436 106 : CALL reallocate(grid_atom%rad, 1, new_nr)
437 106 : CALL reallocate(grid_atom%rad2, 1, new_nr)
438 106 : CALL reallocate(grid_atom%wr, 1, new_nr)
439 106 : CALL reallocate(grid_atom%weight, 1, na, 1, new_nr)
440 106 : CALL reallocate(grid_atom%rad2l, 1, new_nr, 0, llmax_p1)
441 106 : CALL reallocate(grid_atom%oorad2l, 1, new_nr, 0, llmax_p1)
442 :
443 106 : END SUBROUTINE truncate_radial_grid
444 :
445 : ! **************************************************************************************************
446 : !> \brief Computes the contraction coefficients to go from spherical orbitals to sgf for a given
447 : !> atomic kind and a given basis type.
448 : !> \param ikind the kind for which we compute the coefficients
449 : !> \param basis_type the type of basis for which we compute
450 : !> \param sphi_so where the new contraction coefficients are stored (not yet allocated)
451 : !> \param qs_env ...
452 : ! **************************************************************************************************
453 174 : SUBROUTINE compute_sphi_so(ikind, basis_type, sphi_so, qs_env)
454 :
455 : INTEGER, INTENT(IN) :: ikind
456 : CHARACTER(len=*), INTENT(IN) :: basis_type
457 : REAL(dp), DIMENSION(:, :), POINTER :: sphi_so
458 : TYPE(qs_environment_type), POINTER :: qs_env
459 :
460 : INTEGER :: ico, ipgf, iset, iso, l, lx, ly, lz, &
461 : maxso, nset, sgfi, start_c, start_s
462 174 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nsgf_set
463 174 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf
464 : REAL(dp) :: factor
465 174 : REAL(dp), DIMENSION(:, :), POINTER :: sphi
466 : TYPE(gto_basis_set_type), POINTER :: basis
467 174 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
468 :
469 174 : NULLIFY (basis, lmax, lmin, npgf, nsgf_set, qs_kind_set, first_sgf, sphi)
470 :
471 174 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
472 174 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis, basis_type=basis_type)
473 : CALL get_gto_basis_set(basis, lmax=lmax, nset=nset, npgf=npgf, maxso=maxso, lmin=lmin, &
474 174 : nsgf_set=nsgf_set, sphi=sphi, first_sgf=first_sgf)
475 :
476 1356 : ALLOCATE (sphi_so(maxso, SUM(nsgf_set)))
477 625934 : sphi_so = 0.0_dp
478 :
479 834 : DO iset = 1, nset
480 660 : sgfi = first_sgf(1, iset)
481 2948 : DO ipgf = 1, npgf(iset)
482 2114 : start_s = (ipgf - 1)*nsoset(lmax(iset))
483 2114 : start_c = (ipgf - 1)*ncoset(lmax(iset))
484 6452 : DO l = lmin(iset), lmax(iset)
485 15674 : DO iso = 1, nso(l)
486 61618 : DO ico = 1, nco(l)
487 48058 : lx = indco(1, ico + ncoset(l - 1))
488 48058 : ly = indco(2, ico + ncoset(l - 1))
489 48058 : lz = indco(3, ico + ncoset(l - 1))
490 : !MK factor = orbtramat(l)%s2c(iso, ico) &
491 : !MK *SQRT(4.0_dp*pi/dfac(2*l + 1)*dfac(2*lx - 1)*dfac(2*ly - 1)*dfac(2*lz - 1))
492 48058 : factor = orbtramat(l)%slm_inv(iso, ico)
493 : sphi_so(start_s + nsoset(l - 1) + iso, sgfi:sgfi + nsgf_set(iset) - 1) = &
494 : sphi_so(start_s + nsoset(l - 1) + iso, sgfi:sgfi + nsgf_set(iset) - 1) + &
495 5222690 : factor*sphi(start_c + ncoset(l - 1) + ico, sgfi:sgfi + nsgf_set(iset) - 1)
496 : END DO ! ico
497 : END DO ! iso
498 : END DO ! l
499 : END DO ! ipgf
500 : END DO ! iset
501 :
502 348 : END SUBROUTINE compute_sphi_so
503 :
504 : ! **************************************************************************************************
505 : !> \brief Find the neighbors of a given set of atoms based on the non-zero blocks of a provided
506 : !> overlap matrix. Optionally returns an array containing the indices of all involved atoms
507 : !> (the given subset plus all their neighbors, without repetition) AND/OR an array of arrays
508 : !> providing the indices of the neighbors of each input atom.
509 : !> \param base_atoms the set of atoms for which we search neighbors
510 : !> \param mat_s the overlap matrix used to find neighbors
511 : !> \param radius the cutoff radius after which atoms are not considered neighbors
512 : !> \param qs_env ...
513 : !> \param all_neighbors the array uniquely contatining all indices of all atoms involved
514 : !> \param neighbor_set array of arrays containing the neighbors of all given atoms
515 : ! **************************************************************************************************
516 52 : SUBROUTINE find_neighbors(base_atoms, mat_s, radius, qs_env, all_neighbors, neighbor_set)
517 :
518 : INTEGER, DIMENSION(:), INTENT(INOUT) :: base_atoms
519 : TYPE(dbcsr_type), INTENT(IN) :: mat_s
520 : REAL(dp) :: radius
521 : TYPE(qs_environment_type), POINTER :: qs_env
522 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: all_neighbors
523 : TYPE(cp_1d_i_p_type), DIMENSION(:), OPTIONAL, &
524 : POINTER :: neighbor_set
525 :
526 : INTEGER :: i, iat, ibase, iblk, jblk, mepos, natom, &
527 : nb, nbase
528 52 : INTEGER, ALLOCATABLE, DIMENSION(:) :: blk_to_base, inb, who_is_there
529 52 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: n_neighbors
530 52 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: is_base_atom
531 : REAL(dp) :: dist2, rad2, ri(3), rij(3), rj(3)
532 : TYPE(cell_type), POINTER :: cell
533 52 : TYPE(cp_1d_i_p_type), DIMENSION(:), POINTER :: my_neighbor_set
534 : TYPE(dbcsr_iterator_type) :: iter
535 : TYPE(mp_para_env_type), POINTER :: para_env
536 52 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
537 :
538 52 : NULLIFY (particle_set, para_env, my_neighbor_set, cell)
539 :
540 : ! Initialization
541 52 : CALL get_qs_env(qs_env, para_env=para_env, natom=natom, particle_set=particle_set, cell=cell)
542 52 : mepos = para_env%mepos
543 52 : nbase = SIZE(base_atoms)
544 : !work with the neighbor_set structure, see at the end if we keep it
545 218 : ALLOCATE (my_neighbor_set(nbase))
546 52 : rad2 = radius**2
547 :
548 208 : ALLOCATE (blk_to_base(natom), is_base_atom(natom))
549 924 : blk_to_base = 0; is_base_atom = .FALSE.
550 114 : DO ibase = 1, nbase
551 62 : blk_to_base(base_atoms(ibase)) = ibase
552 114 : is_base_atom(base_atoms(ibase)) = .TRUE.
553 : END DO
554 :
555 : ! First loop over S => count the number of neighbors
556 208 : ALLOCATE (n_neighbors(nbase, 0:para_env%num_pe - 1))
557 280 : n_neighbors = 0
558 :
559 52 : CALL dbcsr_iterator_readonly_start(iter, mat_s)
560 3874 : DO WHILE (dbcsr_iterator_blocks_left(iter))
561 :
562 3822 : CALL dbcsr_iterator_next_block(iter, row=iblk, column=jblk)
563 :
564 : !avoid self-neighbors
565 3822 : IF (iblk == jblk) CYCLE
566 :
567 : !test distance
568 3617 : ri = pbc(particle_set(iblk)%r, cell)
569 3617 : rj = pbc(particle_set(jblk)%r, cell)
570 3617 : rij = pbc(ri, rj, cell)
571 14468 : dist2 = SUM(rij**2)
572 3617 : IF (dist2 > rad2) CYCLE
573 :
574 26 : IF (is_base_atom(iblk)) THEN
575 14 : ibase = blk_to_base(iblk)
576 14 : n_neighbors(ibase, mepos) = n_neighbors(ibase, mepos) + 1
577 : END IF
578 78 : IF (is_base_atom(jblk)) THEN
579 3 : ibase = blk_to_base(jblk)
580 3 : n_neighbors(ibase, mepos) = n_neighbors(ibase, mepos) + 1
581 : END IF
582 :
583 : END DO !iter
584 52 : CALL dbcsr_iterator_stop(iter)
585 52 : CALL para_env%sum(n_neighbors)
586 :
587 : ! Allocate the neighbor_set arrays at the correct length
588 114 : DO ibase = 1, nbase
589 260 : ALLOCATE (my_neighbor_set(ibase)%array(SUM(n_neighbors(ibase, :))))
590 148 : my_neighbor_set(ibase)%array = 0
591 : END DO
592 :
593 : ! Loop a second time over S, this time fill the neighbors details
594 52 : CALL dbcsr_iterator_readonly_start(iter, mat_s)
595 156 : ALLOCATE (inb(nbase))
596 114 : inb = 1
597 3874 : DO WHILE (dbcsr_iterator_blocks_left(iter))
598 :
599 3822 : CALL dbcsr_iterator_next_block(iter, row=iblk, column=jblk)
600 3822 : IF (iblk == jblk) CYCLE
601 :
602 : !test distance
603 3617 : ri = pbc(particle_set(iblk)%r, cell)
604 3617 : rj = pbc(particle_set(jblk)%r, cell)
605 3617 : rij = pbc(ri, rj, cell)
606 14468 : dist2 = SUM(rij**2)
607 3617 : IF (dist2 > rad2) CYCLE
608 :
609 26 : IF (is_base_atom(iblk)) THEN
610 14 : ibase = blk_to_base(iblk)
611 23 : my_neighbor_set(ibase)%array(SUM(n_neighbors(ibase, 0:mepos - 1)) + inb(ibase)) = jblk
612 14 : inb(ibase) = inb(ibase) + 1
613 : END IF
614 78 : IF (is_base_atom(jblk)) THEN
615 3 : ibase = blk_to_base(jblk)
616 4 : my_neighbor_set(ibase)%array(SUM(n_neighbors(ibase, 0:mepos - 1)) + inb(ibase)) = iblk
617 3 : inb(ibase) = inb(ibase) + 1
618 : END IF
619 :
620 : END DO !iter
621 52 : CALL dbcsr_iterator_stop(iter)
622 :
623 : ! Make sure the info is shared among the procs
624 114 : DO ibase = 1, nbase
625 182 : CALL para_env%sum(my_neighbor_set(ibase)%array)
626 : END DO
627 :
628 : ! Gather all indices if asked for it
629 52 : IF (PRESENT(all_neighbors)) THEN
630 156 : ALLOCATE (who_is_there(natom))
631 462 : who_is_there = 0
632 :
633 114 : DO ibase = 1, nbase
634 62 : who_is_there(base_atoms(ibase)) = 1
635 148 : DO nb = 1, SIZE(my_neighbor_set(ibase)%array)
636 96 : who_is_there(my_neighbor_set(ibase)%array(nb)) = 1
637 : END DO
638 : END DO
639 :
640 104 : ALLOCATE (all_neighbors(natom))
641 52 : i = 0
642 462 : DO iat = 1, natom
643 462 : IF (who_is_there(iat) == 1) THEN
644 88 : i = i + 1
645 88 : all_neighbors(i) = iat
646 : END IF
647 : END DO
648 52 : CALL reallocate(all_neighbors, 1, i)
649 : END IF
650 :
651 : ! If not asked for the neighbor set, deallocate it
652 52 : IF (PRESENT(neighbor_set)) THEN
653 52 : neighbor_set => my_neighbor_set
654 : ELSE
655 0 : DO ibase = 1, nbase
656 0 : DEALLOCATE (my_neighbor_set(ibase)%array)
657 : END DO
658 0 : DEALLOCATE (my_neighbor_set)
659 : END IF
660 :
661 208 : END SUBROUTINE find_neighbors
662 :
663 : ! **************************************************************************************************
664 : !> \brief Returns the RI inverse overlap for a subset of the RI_XAS matrix spaning a given
665 : !> excited atom and its neighbors.
666 : !> \param ri_sinv the inverse overlap as a dbcsr matrix
667 : !> \param whole_s the whole RI overlap matrix
668 : !> \param neighbors the indeces of the excited atom and their neighbors
669 : !> \param idx_to_nb array telling where any atom can be found in neighbors (if there at all)
670 : !> \param basis_set_ri the RI basis set list for all kinds
671 : !> \param qs_env ...
672 : !> \note It is assumed that the neighbors are sorted, the output matrix is assumed to be small and
673 : !> is replicated on all processors
674 : ! **************************************************************************************************
675 62 : SUBROUTINE get_exat_ri_sinv(ri_sinv, whole_s, neighbors, idx_to_nb, basis_set_ri, qs_env)
676 :
677 : TYPE(dbcsr_type) :: ri_sinv, whole_s
678 : INTEGER, DIMENSION(:), INTENT(IN) :: neighbors, idx_to_nb
679 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_ri
680 : TYPE(qs_environment_type), POINTER :: qs_env
681 :
682 : CHARACTER(len=*), PARAMETER :: routineN = 'get_exat_ri_sinv'
683 :
684 : INTEGER :: blk, dest, group_handle, handle, iat, &
685 : ikind, inb, ir, is, jat, jnb, natom, &
686 : nnb, npcols, nprows, source, tag
687 62 : INTEGER, DIMENSION(:), POINTER :: col_dist, nsgf, row_dist
688 62 : INTEGER, DIMENSION(:, :), POINTER :: pgrid
689 : LOGICAL :: found_risinv, found_whole
690 62 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: is_neighbor
691 : REAL(dp) :: ri(3), rij(3), rj(3)
692 62 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: radius
693 62 : REAL(dp), DIMENSION(:, :), POINTER :: block_risinv, block_whole
694 : TYPE(cell_type), POINTER :: cell
695 62 : TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:) :: recv_buff, send_buff
696 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
697 : TYPE(dbcsr_distribution_type) :: sinv_dist
698 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
699 : TYPE(dbcsr_iterator_type) :: iter
700 : TYPE(mp_comm_type) :: group
701 : TYPE(mp_para_env_type), POINTER :: para_env
702 62 : TYPE(mp_request_type), ALLOCATABLE, DIMENSION(:) :: recv_req, send_req
703 62 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
704 :
705 62 : NULLIFY (pgrid, dbcsr_dist, row_dist, col_dist, nsgf, particle_set, block_whole, block_risinv)
706 62 : NULLIFY (cell, para_env, blacs_env)
707 :
708 62 : CALL timeset(routineN, handle)
709 :
710 : CALL get_qs_env(qs_env, dbcsr_dist=dbcsr_dist, particle_set=particle_set, natom=natom, &
711 62 : para_env=para_env, blacs_env=blacs_env, cell=cell)
712 62 : nnb = SIZE(neighbors)
713 434 : ALLOCATE (nsgf(nnb), is_neighbor(natom), radius(nnb))
714 686 : is_neighbor = .FALSE.
715 158 : DO inb = 1, nnb
716 96 : iat = neighbors(inb)
717 96 : ikind = particle_set(iat)%atomic_kind%kind_number
718 96 : CALL get_gto_basis_set(basis_set_ri(ikind)%gto_basis_set, nsgf=nsgf(inb), kind_radius=radius(inb))
719 158 : is_neighbor(iat) = .TRUE.
720 : END DO
721 :
722 : !Create the ri_sinv matrix based on some arbitrary dbcsr_dist
723 62 : CALL dbcsr_distribution_get(dbcsr_dist, group=group_handle, pgrid=pgrid, nprows=nprows, npcols=npcols)
724 62 : CALL group%set_handle(group_handle)
725 :
726 186 : ALLOCATE (row_dist(nnb), col_dist(nnb))
727 158 : DO inb = 1, nnb
728 96 : row_dist(inb) = MODULO(nprows - inb, nprows)
729 158 : col_dist(inb) = MODULO(npcols - inb, npcols)
730 : END DO
731 :
732 : CALL dbcsr_distribution_new(sinv_dist, group=group_handle, pgrid=pgrid, row_dist=row_dist, &
733 62 : col_dist=col_dist)
734 :
735 : CALL dbcsr_create(matrix=ri_sinv, name="RI_SINV", matrix_type=dbcsr_type_symmetric, &
736 62 : dist=sinv_dist, row_blk_size=nsgf, col_blk_size=nsgf)
737 : !reserving the blocks in the correct pattern
738 158 : DO inb = 1, nnb
739 96 : ri = pbc(particle_set(neighbors(inb))%r, cell)
740 320 : DO jnb = inb, nnb
741 :
742 : !do the atom overlap ?
743 162 : rj = pbc(particle_set(neighbors(jnb))%r, cell)
744 162 : rij = pbc(ri, rj, cell)
745 648 : IF (SUM(rij**2) > (radius(inb) + radius(jnb))**2) CYCLE
746 :
747 162 : CALL dbcsr_get_stored_coordinates(ri_sinv, inb, jnb, blk)
748 258 : IF (para_env%mepos == blk) THEN
749 324 : ALLOCATE (block_risinv(nsgf(inb), nsgf(jnb)))
750 276734 : block_risinv = 0.0_dp
751 81 : CALL dbcsr_put_block(ri_sinv, inb, jnb, block_risinv)
752 81 : DEALLOCATE (block_risinv)
753 : END IF
754 : END DO
755 : END DO
756 62 : CALL dbcsr_finalize(ri_sinv)
757 :
758 62 : CALL dbcsr_distribution_release(sinv_dist)
759 62 : DEALLOCATE (row_dist, col_dist)
760 :
761 : !prepare the send and recv buffers we will need for change of dist between the two matrices
762 : !worst case scenario: all neighbors are on same procs => need to send nnb**2 messages
763 766 : ALLOCATE (send_buff(nnb**2), recv_buff(nnb**2))
764 828 : ALLOCATE (send_req(nnb**2), recv_req(nnb**2))
765 62 : is = 0; ir = 0
766 :
767 : !Loop over the whole RI overlap matrix and pick the blocks we need
768 62 : CALL dbcsr_iterator_start(iter, whole_s)
769 7468 : DO WHILE (dbcsr_iterator_blocks_left(iter))
770 :
771 7406 : CALL dbcsr_iterator_next_block(iter, row=iat, column=jat)
772 7406 : CALL dbcsr_get_block_p(whole_s, iat, jat, block_whole, found_whole)
773 :
774 : !only interested in neighbors
775 7406 : IF (.NOT. found_whole) CYCLE
776 7406 : IF (.NOT. is_neighbor(iat)) CYCLE
777 240 : IF (.NOT. is_neighbor(jat)) CYCLE
778 :
779 81 : inb = idx_to_nb(iat)
780 81 : jnb = idx_to_nb(jat)
781 :
782 : !If blocks are on the same proc for both matrices, simply copy
783 81 : CALL dbcsr_get_block_p(ri_sinv, inb, jnb, block_risinv, found_risinv)
784 143 : IF (found_risinv) THEN
785 11 : CALL dcopy(nsgf(inb)*nsgf(jnb), block_whole, 1, block_risinv, 1)
786 : ELSE
787 :
788 : !send the block with unique tag to the proc where inb,jnb is in ri_sinv
789 70 : CALL dbcsr_get_stored_coordinates(ri_sinv, inb, jnb, dest)
790 70 : is = is + 1
791 70 : send_buff(is)%array => block_whole
792 70 : tag = natom*inb + jnb
793 70 : CALL group%isend(msgin=send_buff(is)%array, dest=dest, request=send_req(is), tag=tag)
794 :
795 : END IF
796 :
797 : END DO !dbcsr iter
798 62 : CALL dbcsr_iterator_stop(iter)
799 :
800 : !Loop over ri_sinv and receive all those blocks
801 62 : CALL dbcsr_iterator_start(iter, ri_sinv)
802 143 : DO WHILE (dbcsr_iterator_blocks_left(iter))
803 :
804 81 : CALL dbcsr_iterator_next_block(iter, row=inb, column=jnb)
805 81 : CALL dbcsr_get_block_p(ri_sinv, inb, jnb, block_risinv, found_risinv)
806 :
807 81 : IF (.NOT. found_risinv) CYCLE
808 :
809 81 : iat = neighbors(inb)
810 81 : jat = neighbors(jnb)
811 :
812 : !If blocks are on the same proc on both matrices do nothing
813 81 : CALL dbcsr_get_stored_coordinates(whole_s, iat, jat, source)
814 81 : IF (para_env%mepos == source) CYCLE
815 :
816 70 : tag = natom*inb + jnb
817 70 : ir = ir + 1
818 70 : recv_buff(ir)%array => block_risinv
819 : CALL group%irecv(msgout=recv_buff(ir)%array, source=source, request=recv_req(ir), &
820 143 : tag=tag)
821 :
822 : END DO
823 62 : CALL dbcsr_iterator_stop(iter)
824 :
825 : !make sure that all comm is over before proceeding
826 62 : CALL mp_waitall(send_req(1:is))
827 62 : CALL mp_waitall(recv_req(1:ir))
828 :
829 : !Invert. 2 cases: with or without neighbors. If no neighbors, easier to invert on one proc and
830 : !avoid the whole fm to dbcsr to fm that is quite expensive
831 62 : IF (nnb == 1) THEN
832 :
833 50 : CALL dbcsr_get_block_p(ri_sinv, 1, 1, block_risinv, found_risinv)
834 50 : IF (found_risinv) THEN
835 25 : CALL invmat_symm(block_risinv)
836 : END IF
837 :
838 : ELSE
839 12 : CALL cp_dbcsr_cholesky_decompose(ri_sinv, para_env=para_env, blacs_env=blacs_env)
840 12 : CALL cp_dbcsr_cholesky_invert(ri_sinv, para_env=para_env, blacs_env=blacs_env, uplo_to_full=.TRUE.)
841 12 : CALL dbcsr_filter(ri_sinv, 1.E-10_dp) !make sure ri_sinv is sparse coming out of fm routines
842 : END IF
843 62 : CALL dbcsr_replicate_all(ri_sinv)
844 :
845 : !clean-up
846 62 : DEALLOCATE (nsgf)
847 :
848 62 : CALL timestop(handle)
849 :
850 372 : END SUBROUTINE get_exat_ri_sinv
851 :
852 : ! **************************************************************************************************
853 : !> \brief Compute the coefficients to project the density on a partial RI_XAS basis
854 : !> \param xas_atom_env ...
855 : !> \param qs_env ...
856 : !> \note The density is n = sum_ab P_ab*phi_a*phi_b, the RI basis covers the products of orbital sgfs
857 : !> => n = sum_ab sum_cd P_ab (phi_a phi_b xi_c) S_cd^-1 xi_d
858 : !> = sum_d coeff_d xi_d , where xi are the RI basis func.
859 : !> In this case, with the partial RI projection, the RI basis is restricted to an excited atom
860 : !> and its neighbors at a time. Leads to smaller overlap matrix to invert and less 3-center
861 : !> overlap to compute. The procedure is repeated for each excited atom
862 : ! **************************************************************************************************
863 52 : SUBROUTINE calculate_density_coeffs(xas_atom_env, qs_env)
864 :
865 : TYPE(xas_atom_env_type), POINTER :: xas_atom_env
866 : TYPE(qs_environment_type), POINTER :: qs_env
867 :
868 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_density_coeffs'
869 :
870 : INTEGER :: exat, handle, i, iat, iatom, iex, inb, &
871 : ind(3), ispin, jatom, jnb, katom, &
872 : natom, nex, nkind, nnb, nspins, &
873 : output_unit, ri_at
874 52 : INTEGER, ALLOCATABLE, DIMENSION(:) :: blk_size_orb, blk_size_ri, idx_to_nb, &
875 52 : neighbors
876 52 : INTEGER, DIMENSION(:), POINTER :: all_ri_atoms
877 : LOGICAL :: pmat_found, pmat_foundt, sinv_found, &
878 : sinv_foundt, tensor_found, unique
879 : REAL(dp) :: factor, prefac
880 52 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: work2
881 52 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: work1
882 52 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: t_block
883 104 : REAL(dp), DIMENSION(:, :), POINTER :: pmat_block, pmat_blockt, sinv_block, &
884 52 : sinv_blockt
885 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
886 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
887 104 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: overlap, rho_ao
888 : TYPE(dbcsr_type) :: ri_sinv
889 : TYPE(dbcsr_type), POINTER :: ri_mats
890 : TYPE(dbt_iterator_type) :: iter
891 468 : TYPE(dbt_type) :: pqX
892 52 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_orb, basis_set_ri
893 : TYPE(libint_potential_type) :: pot
894 : TYPE(mp_para_env_type), POINTER :: para_env
895 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
896 104 : POINTER :: ab_list, ac_list, sab_ri
897 52 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
898 52 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
899 : TYPE(qs_rho_type), POINTER :: rho
900 :
901 52 : NULLIFY (qs_kind_set, basis_set_ri, basis_set_orb, ac_list, rho, rho_ao, sab_ri, ri_mats)
902 52 : NULLIFY (particle_set, para_env, all_ri_atoms, overlap, pmat_blockt)
903 52 : NULLIFY (blacs_env, pmat_block, ab_list, dbcsr_dist, sinv_block, sinv_blockt)
904 :
905 : !Idea: We don't do a full RI here as it would be too expensive in many ways (inversion of a
906 : ! large matrix, many 3-center overlaps, density coefficients for all atoms, etc...)
907 : ! Instead, we go excited atom by excited atom and only do a RI expansion on its basis
908 : ! and that of its closest neighbors (defined through RI_RADIUS), such that we only have
909 : ! very small matrices to invert and only a few (abc) overlp integrals with c on the
910 : ! excited atom its neighbors. This is physically sound since we only need the density
911 : ! well defined on the excited atom grid as we do (aI|P)*(P|Q)^-1*(Q|fxc|R)*(R|S)^-1*(S|Jb)
912 :
913 52 : CALL timeset(routineN, handle)
914 :
915 : CALL get_qs_env(qs_env, nkind=nkind, qs_kind_set=qs_kind_set, rho=rho, &
916 : natom=natom, particle_set=particle_set, dbcsr_dist=dbcsr_dist, &
917 52 : para_env=para_env, blacs_env=blacs_env)
918 52 : nspins = xas_atom_env%nspins
919 52 : nex = SIZE(xas_atom_env%excited_atoms)
920 52 : CALL qs_rho_get(rho, rho_ao=rho_ao)
921 :
922 : ! Create the needed neighbor list and basis set lists.
923 240 : ALLOCATE (basis_set_ri(nkind))
924 188 : ALLOCATE (basis_set_orb(nkind))
925 52 : CALL basis_set_list_setup(basis_set_ri, "RI_XAS", qs_kind_set)
926 52 : CALL basis_set_list_setup(basis_set_orb, "ORB", qs_kind_set)
927 :
928 : ! Compute the RI overlap matrix on the whole system
929 52 : CALL build_xas_tdp_ovlp_nl(sab_ri, basis_set_ri, basis_set_ri, qs_env)
930 52 : CALL build_overlap_matrix_simple(qs_env%ks_env, overlap, basis_set_ri, basis_set_ri, sab_ri)
931 52 : ri_mats => overlap(1)%matrix
932 52 : CALL release_neighbor_list_sets(sab_ri)
933 :
934 : ! Get the neighbors of the excited atoms (= all the atoms where density coeffs are needed)
935 : CALL find_neighbors(xas_atom_env%excited_atoms, ri_mats, xas_atom_env%ri_radius, &
936 52 : qs_env, all_neighbors=all_ri_atoms, neighbor_set=xas_atom_env%exat_neighbors)
937 :
938 : !keep in mind that double occupation is included in rho_ao in case of closed-shell
939 52 : factor = 0.5_dp; IF (nspins == 2) factor = 1.0_dp
940 :
941 : ! Allocate space for the projected density coefficients. On all ri_atoms.
942 : ! Note: the sub-region where we project the density changes from excited atom to excited atom
943 : ! => need different sets of RI coeffs
944 156 : ALLOCATE (blk_size_ri(natom))
945 52 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=blk_size_ri, basis=basis_set_ri)
946 1038 : ALLOCATE (xas_atom_env%ri_dcoeff(natom, nspins, nex))
947 114 : DO iex = 1, nex
948 182 : DO ispin = 1, nspins
949 778 : DO iat = 1, natom
950 648 : NULLIFY (xas_atom_env%ri_dcoeff(iat, ispin, iex)%array)
951 : IF ((.NOT. ANY(xas_atom_env%exat_neighbors(iex)%array == iat)) &
952 750 : .AND. (.NOT. xas_atom_env%excited_atoms(iex) == iat)) CYCLE
953 360 : ALLOCATE (xas_atom_env%ri_dcoeff(iat, ispin, iex)%array(blk_size_ri(iat)))
954 6956 : xas_atom_env%ri_dcoeff(iat, ispin, iex)%array = 0.0_dp
955 : END DO
956 : END DO
957 : END DO
958 :
959 52 : output_unit = cp_logger_get_default_io_unit()
960 52 : IF (output_unit > 0) THEN
961 : WRITE (output_unit, FMT="(/,T7,A,/,T7,A)") &
962 26 : "Excited atom, natoms in RI_REGION:", &
963 52 : "---------------------------------"
964 : END IF
965 :
966 : !We go atom by atom, first computing the integrals themselves that we put into a tensor, then we do
967 : !the contraction with the density. We do that in the original dist, which is optimized for overlap
968 :
969 156 : ALLOCATE (blk_size_orb(natom))
970 52 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=blk_size_orb, basis=basis_set_orb)
971 :
972 114 : DO iex = 1, nex
973 :
974 : !get neighbors of current atom
975 62 : exat = xas_atom_env%excited_atoms(iex)
976 62 : nnb = 1 + SIZE(xas_atom_env%exat_neighbors(iex)%array)
977 186 : ALLOCATE (neighbors(nnb))
978 62 : neighbors(1) = exat
979 96 : neighbors(2:nnb) = xas_atom_env%exat_neighbors(iex)%array(:)
980 62 : CALL sort_unique(neighbors, unique)
981 :
982 : !link the atoms to their position in neighbors
983 186 : ALLOCATE (idx_to_nb(natom))
984 686 : idx_to_nb = 0
985 158 : DO inb = 1, nnb
986 158 : idx_to_nb(neighbors(inb)) = inb
987 : END DO
988 :
989 62 : IF (output_unit > 0) THEN
990 : WRITE (output_unit, FMT="(T7,I12,I21)") &
991 31 : exat, nnb
992 : END IF
993 :
994 : !Get the neighbor lists for the overlap integrals (abc), centers c on the current
995 : !excited atom and its neighbors defined by RI_RADIUS
996 62 : CALL build_xas_tdp_ovlp_nl(ab_list, basis_set_orb, basis_set_orb, qs_env)
997 : CALL build_xas_tdp_3c_nl(ac_list, basis_set_orb, basis_set_ri, do_potential_id, &
998 62 : qs_env, excited_atoms=neighbors)
999 :
1000 : !Compute the 3-center overlap integrals
1001 62 : pot%potential_type = do_potential_id
1002 :
1003 : CALL create_pqX_tensor(pqX, ab_list, ac_list, dbcsr_dist, blk_size_orb, blk_size_orb, &
1004 62 : blk_size_ri)
1005 : CALL fill_pqX_tensor(pqX, ab_list, ac_list, basis_set_orb, basis_set_orb, basis_set_ri, &
1006 62 : pot, qs_env)
1007 :
1008 : !Compute the RI inverse overlap matrix on the reduced RI basis that spans the excited
1009 : !atom and its neighbors, ri_sinv is replicated over all procs
1010 62 : CALL get_exat_ri_sinv(ri_sinv, ri_mats, neighbors, idx_to_nb, basis_set_ri, qs_env)
1011 :
1012 : !Do the actual contraction: coeff_y = sum_pq sum_x P_pq (phi_p phi_q xi_x) S_xy^-1
1013 :
1014 : !$OMP PARALLEL DEFAULT(NONE) &
1015 : !$OMP SHARED(pqX,rho_ao,ri_sinv,xas_atom_env) &
1016 : !$OMP SHARED(blk_size_ri,idx_to_nb,nspins,nnb,neighbors,iex,factor) &
1017 : !$OMP PRIVATE(iter,ind,t_block,tensor_found,iatom,jatom,katom,inb,prefac,ispin) &
1018 : !$OMP PRIVATE(pmat_block,pmat_found,pmat_blockt,pmat_foundt,work1,work2,jnb,ri_at) &
1019 62 : !$OMP PRIVATE(sinv_block,sinv_found,sinv_blockt,sinv_foundt)
1020 : CALL dbt_iterator_start(iter, pqX)
1021 : DO WHILE (dbt_iterator_blocks_left(iter))
1022 : CALL dbt_iterator_next_block(iter, ind)
1023 : CALL dbt_get_block(pqX, ind, t_block, tensor_found)
1024 :
1025 : iatom = ind(1)
1026 : jatom = ind(2)
1027 : katom = ind(3)
1028 : inb = idx_to_nb(katom)
1029 :
1030 : !non-diagonal elements need to be counted twice
1031 : prefac = 2.0_dp
1032 : IF (iatom == jatom) prefac = 1.0_dp
1033 :
1034 : DO ispin = 1, nspins
1035 :
1036 : !rho_ao is symmetric, block can be in either location
1037 : CALL dbcsr_get_block_p(rho_ao(ispin)%matrix, iatom, jatom, pmat_block, pmat_found)
1038 : CALL dbcsr_get_block_p(rho_ao(ispin)%matrix, jatom, iatom, pmat_blockt, pmat_foundt)
1039 : IF ((.NOT. pmat_found) .AND. (.NOT. pmat_foundt)) CYCLE
1040 :
1041 : ALLOCATE (work1(blk_size_ri(katom), 1))
1042 : work1 = 0.0_dp
1043 :
1044 : !first contraction with the density matrix
1045 : IF (pmat_found) THEN
1046 : DO i = 1, blk_size_ri(katom)
1047 : work1(i, 1) = prefac*SUM(pmat_block(:, :)*t_block(:, :, i))
1048 : END DO
1049 : ELSE
1050 : DO i = 1, blk_size_ri(katom)
1051 : work1(i, 1) = prefac*SUM(TRANSPOSE(pmat_blockt(:, :))*t_block(:, :, i))
1052 : END DO
1053 : END IF
1054 :
1055 : !loop over neighbors
1056 : DO jnb = 1, nnb
1057 :
1058 : ri_at = neighbors(jnb)
1059 :
1060 : !ri_sinv is a symmetric matrix => actual block is one of the two
1061 : CALL dbcsr_get_block_p(ri_sinv, inb, jnb, sinv_block, sinv_found)
1062 : CALL dbcsr_get_block_p(ri_sinv, jnb, inb, sinv_blockt, sinv_foundt)
1063 : IF ((.NOT. sinv_found) .AND. (.NOT. sinv_foundt)) CYCLE
1064 :
1065 : !second contraction with the inverse RI overlap
1066 : ALLOCATE (work2(SIZE(xas_atom_env%ri_dcoeff(ri_at, ispin, iex)%array)))
1067 : work2 = 0.0_dp
1068 :
1069 : IF (sinv_found) THEN
1070 : DO i = 1, blk_size_ri(katom)
1071 : work2(:) = work2(:) + factor*work1(i, 1)*sinv_block(i, :)
1072 : END DO
1073 : ELSE
1074 : DO i = 1, blk_size_ri(katom)
1075 : work2(:) = work2(:) + factor*work1(i, 1)*sinv_blockt(:, i)
1076 : END DO
1077 : END IF
1078 : DO i = 1, SIZE(work2)
1079 : !$OMP ATOMIC
1080 : xas_atom_env%ri_dcoeff(ri_at, ispin, iex)%array(i) = &
1081 : xas_atom_env%ri_dcoeff(ri_at, ispin, iex)%array(i) + work2(i)
1082 : END DO
1083 :
1084 : DEALLOCATE (work2)
1085 : END DO !jnb
1086 :
1087 : DEALLOCATE (work1)
1088 : END DO
1089 :
1090 : DEALLOCATE (t_block)
1091 : END DO !iter
1092 : CALL dbt_iterator_stop(iter)
1093 : !$OMP END PARALLEL
1094 :
1095 : !clean-up
1096 62 : CALL dbcsr_release(ri_sinv)
1097 62 : CALL dbt_destroy(pqX)
1098 62 : CALL release_neighbor_list_sets(ab_list)
1099 62 : CALL release_neighbor_list_sets(ac_list)
1100 176 : DEALLOCATE (neighbors, idx_to_nb)
1101 :
1102 : END DO !iex
1103 :
1104 : !making sure all procs have the same info
1105 114 : DO iex = 1, nex
1106 182 : DO ispin = 1, nspins
1107 778 : DO iat = 1, natom
1108 : IF ((.NOT. ANY(xas_atom_env%exat_neighbors(iex)%array == iat)) &
1109 750 : .AND. (.NOT. xas_atom_env%excited_atoms(iex) == iat)) CYCLE
1110 14252 : CALL para_env%sum(xas_atom_env%ri_dcoeff(iat, ispin, iex)%array)
1111 : END DO !iat
1112 : END DO !ispin
1113 : END DO !iex
1114 :
1115 : ! clean-up
1116 52 : CALL dbcsr_deallocate_matrix_set(overlap)
1117 52 : DEALLOCATE (basis_set_ri, basis_set_orb, all_ri_atoms)
1118 :
1119 52 : CALL timestop(handle)
1120 :
1121 156 : END SUBROUTINE calculate_density_coeffs
1122 :
1123 : ! **************************************************************************************************
1124 : !> \brief Evaluates the density on a given atomic grid
1125 : !> \param rho_set where the densities are stored
1126 : !> \param ri_dcoeff the arrays containing the RI density coefficients of this atom, for each spin
1127 : !> \param atom_kind the kind of the atom in question
1128 : !> \param do_gga whether the gradient of the density should also be put on the grid
1129 : !> \param batch_info how the so are distributed
1130 : !> \param xas_atom_env ...
1131 : !> \param qs_env ...
1132 : !> \note The density is expressed as n = sum_d coeff_d*xi_d. Knowing the coordinate of each grid
1133 : !> grid point, one can simply evaluate xi_d(r)
1134 : ! **************************************************************************************************
1135 52 : SUBROUTINE put_density_on_atomic_grid(rho_set, ri_dcoeff, atom_kind, do_gga, batch_info, &
1136 : xas_atom_env, qs_env)
1137 :
1138 : TYPE(xc_rho_set_type), INTENT(INOUT) :: rho_set
1139 : TYPE(cp_1d_r_p_type), DIMENSION(:) :: ri_dcoeff
1140 : INTEGER, INTENT(IN) :: atom_kind
1141 : LOGICAL, INTENT(IN) :: do_gga
1142 : TYPE(batch_info_type) :: batch_info
1143 : TYPE(xas_atom_env_type), POINTER :: xas_atom_env
1144 : TYPE(qs_environment_type), POINTER :: qs_env
1145 :
1146 : CHARACTER(len=*), PARAMETER :: routineN = 'put_density_on_atomic_grid'
1147 :
1148 : INTEGER :: dir, handle, ipgf, iset, iso, iso_proc, &
1149 : ispin, maxso, n, na, nr, nset, nsgfi, &
1150 : nsoi, nspins, sgfi, starti
1151 52 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nsgf_set
1152 52 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf
1153 52 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: so
1154 52 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: dso
1155 52 : REAL(dp), DIMENSION(:, :), POINTER :: dgr1, dgr2, ga, gr, ri_sphi_so, zet
1156 52 : REAL(dp), DIMENSION(:, :, :), POINTER :: dga1, dga2, rhoa, rhob
1157 52 : TYPE(cp_1d_r_p_type), DIMENSION(:), POINTER :: ri_dcoeff_so
1158 : TYPE(grid_atom_type), POINTER :: grid_atom
1159 : TYPE(gto_basis_set_type), POINTER :: ri_basis
1160 52 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1161 :
1162 52 : NULLIFY (grid_atom, ri_basis, qs_kind_set, lmax, npgf, zet, nsgf_set, ri_sphi_so)
1163 52 : NULLIFY (lmin, first_sgf, rhoa, rhob, ga, gr, dgr1, dgr2, dga1, dga2, ri_dcoeff_so)
1164 :
1165 52 : CALL timeset(routineN, handle)
1166 :
1167 : ! Strategy: it makes sense to evaluate the spherical orbital on the grid (because of symmetry)
1168 : ! From there, one can directly contract into sgf using ri_sphi_so and then take the weight
1169 : ! The spherical orbital were precomputed and split in a purely radial and a purely
1170 : ! angular part. The full values on each grid point are obtain through gemm
1171 :
1172 : ! Generalities
1173 52 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
1174 52 : CALL get_qs_kind(qs_kind_set(atom_kind), basis_set=ri_basis, basis_type="RI_XAS")
1175 : CALL get_gto_basis_set(ri_basis, lmax=lmax, npgf=npgf, zet=zet, nset=nset, nsgf_set=nsgf_set, &
1176 52 : first_sgf=first_sgf, lmin=lmin, maxso=maxso)
1177 :
1178 : ! Get the grid and the info we need from it
1179 52 : grid_atom => xas_atom_env%grid_atom_set(atom_kind)%grid_atom
1180 52 : na = grid_atom%ng_sphere
1181 52 : nr = grid_atom%nr
1182 52 : n = na*nr
1183 52 : nspins = xas_atom_env%nspins
1184 52 : ri_sphi_so => xas_atom_env%ri_sphi_so(atom_kind)%array
1185 :
1186 : ! Point to the rho_set densities
1187 52 : rhoa => rho_set%rhoa
1188 52 : rhob => rho_set%rhob
1189 2983732 : rhoa = 0.0_dp; rhob = 0.0_dp;
1190 52 : IF (do_gga) THEN
1191 112 : DO dir = 1, 3
1192 2709744 : rho_set%drhoa(dir)%array = 0.0_dp
1193 2709772 : rho_set%drhob(dir)%array = 0.0_dp
1194 : END DO
1195 : END IF
1196 :
1197 : ! Point to the precomputed SO
1198 52 : ga => xas_atom_env%ga(atom_kind)%array
1199 52 : gr => xas_atom_env%gr(atom_kind)%array
1200 52 : IF (do_gga) THEN
1201 28 : dga1 => xas_atom_env%dga1(atom_kind)%array
1202 28 : dga2 => xas_atom_env%dga2(atom_kind)%array
1203 28 : dgr1 => xas_atom_env%dgr1(atom_kind)%array
1204 28 : dgr2 => xas_atom_env%dgr2(atom_kind)%array
1205 : ELSE
1206 24 : dga1 => xas_atom_env%dga1(atom_kind)%array
1207 24 : dga2 => xas_atom_env%dga2(atom_kind)%array
1208 24 : dgr1 => xas_atom_env%dgr1(atom_kind)%array
1209 24 : dgr2 => xas_atom_env%dgr2(atom_kind)%array
1210 : END IF
1211 :
1212 : ! Need to express the ri_dcoeffs in terms of so (and not sgf)
1213 214 : ALLOCATE (ri_dcoeff_so(nspins))
1214 110 : DO ispin = 1, nspins
1215 174 : ALLOCATE (ri_dcoeff_so(ispin)%array(nset*maxso))
1216 12046 : ri_dcoeff_so(ispin)%array = 0.0_dp
1217 :
1218 : !for a given so, loop over sgf and sum
1219 302 : DO iset = 1, nset
1220 192 : sgfi = first_sgf(1, iset)
1221 192 : nsoi = npgf(iset)*nsoset(lmax(iset))
1222 192 : nsgfi = nsgf_set(iset)
1223 :
1224 : CALL dgemv('N', nsoi, nsgfi, 1.0_dp, ri_sphi_so(1:nsoi, sgfi:sgfi + nsgfi - 1), nsoi, &
1225 : ri_dcoeff(ispin)%array(sgfi:sgfi + nsgfi - 1), 1, 0.0_dp, &
1226 250 : ri_dcoeff_so(ispin)%array((iset - 1)*maxso + 1:(iset - 1)*maxso + nsoi), 1)
1227 :
1228 : END DO
1229 : END DO
1230 :
1231 : !allocate space to store the spherical orbitals on the grid
1232 208 : ALLOCATE (so(na, nr))
1233 164 : IF (do_gga) ALLOCATE (dso(na, nr, 3))
1234 :
1235 : ! Loop over the spherical orbitals on this proc
1236 6034 : DO iso_proc = 1, batch_info%nso_proc(atom_kind)
1237 5982 : iset = batch_info%so_proc_info(atom_kind)%array(1, iso_proc)
1238 5982 : ipgf = batch_info%so_proc_info(atom_kind)%array(2, iso_proc)
1239 5982 : iso = batch_info%so_proc_info(atom_kind)%array(3, iso_proc)
1240 5982 : IF (iso < 0) CYCLE
1241 :
1242 3024 : starti = (iset - 1)*maxso + (ipgf - 1)*nsoset(lmax(iset))
1243 :
1244 : !the spherical orbital on the grid
1245 : CALL dgemm('N', 'T', na, nr, 1, 1.0_dp, ga(:, iso_proc:iso_proc), na, &
1246 3024 : gr(:, iso_proc:iso_proc), nr, 0.0_dp, so(:, :), na)
1247 :
1248 : !the gradient on the grid
1249 3024 : IF (do_gga) THEN
1250 :
1251 7564 : DO dir = 1, 3
1252 : CALL dgemm('N', 'T', na, nr, 1, 1.0_dp, dga1(:, iso_proc:iso_proc, dir), na, &
1253 5673 : dgr1(:, iso_proc:iso_proc), nr, 0.0_dp, dso(:, :, dir), na)
1254 : CALL dgemm('N', 'T', na, nr, 1, 1.0_dp, dga2(:, iso_proc:iso_proc, dir), na, &
1255 7564 : dgr2(:, iso_proc:iso_proc), nr, 1.0_dp, dso(:, :, dir), na)
1256 : END DO
1257 : END IF
1258 :
1259 : !put the so on the grid with the approriate coefficients and sum
1260 3024 : CALL daxpy(n, ri_dcoeff_so(1)%array(starti + iso), so, 1, rhoa(:, :, 1), 1)
1261 :
1262 3024 : IF (nspins == 2) THEN
1263 378 : CALL daxpy(n, ri_dcoeff_so(2)%array(starti + iso), so, 1, rhob(:, :, 1), 1)
1264 : END IF
1265 :
1266 3076 : IF (do_gga) THEN
1267 :
1268 : !put the gradient of the so on the grid with correspond RI coeff
1269 7564 : DO dir = 1, 3
1270 : CALL daxpy(n, ri_dcoeff_so(1)%array(starti + iso), dso(:, :, dir), &
1271 5673 : 1, rho_set%drhoa(dir)%array(:, :, 1), 1)
1272 :
1273 7564 : IF (nspins == 2) THEN
1274 : CALL daxpy(n, ri_dcoeff_so(2)%array(starti + iso), dso(:, :, dir), &
1275 1134 : 1, rho_set%drhob(dir)%array(:, :, 1), 1)
1276 : END IF
1277 : END DO !dir
1278 : END IF !do_gga
1279 :
1280 : END DO
1281 :
1282 : ! Treat spin restricted case (=> copy alpha into beta)
1283 52 : IF (nspins == 1) THEN
1284 46 : CALL dcopy(n, rhoa(:, :, 1), 1, rhob(:, :, 1), 1)
1285 :
1286 46 : IF (do_gga) THEN
1287 88 : DO dir = 1, 3
1288 88 : CALL dcopy(n, rho_set%drhoa(dir)%array(:, :, 1), 1, rho_set%drhob(dir)%array(:, :, 1), 1)
1289 : END DO
1290 : END IF
1291 : END IF
1292 :
1293 : ! Note: sum over procs is done outside
1294 :
1295 : ! clean-up
1296 110 : DO ispin = 1, nspins
1297 110 : DEALLOCATE (ri_dcoeff_so(ispin)%array)
1298 : END DO
1299 52 : DEALLOCATE (ri_dcoeff_so)
1300 :
1301 52 : CALL timestop(handle)
1302 :
1303 156 : END SUBROUTINE put_density_on_atomic_grid
1304 :
1305 : ! **************************************************************************************************
1306 : !> \brief Adds the density of a given source atom with source kind (with ri_dcoeff) on the atomic
1307 : !> grid belonging to another target atom of target kind. The evaluations of the basis
1308 : !> function first requires the evaluation of the x,y,z coordinates on each grid point of
1309 : !> target atom wrt to the position of source atom
1310 : !> \param rho_set where the densities are stored
1311 : !> \param ri_dcoeff the arrays containing the RI density coefficient of source_iat, for each spin
1312 : !> \param source_iat the index of the source atom
1313 : !> \param source_ikind the kind of the source atom
1314 : !> \param target_iat the index of the target atom
1315 : !> \param target_ikind the kind of the target atom
1316 : !> \param sr starting r index for the local grid
1317 : !> \param er ending r index for the local grid
1318 : !> \param do_gga whether the gradient of the density is needed
1319 : !> \param xas_atom_env ...
1320 : !> \param qs_env ...
1321 : ! **************************************************************************************************
1322 28 : SUBROUTINE put_density_on_other_grid(rho_set, ri_dcoeff, source_iat, source_ikind, target_iat, &
1323 : target_ikind, sr, er, do_gga, xas_atom_env, qs_env)
1324 :
1325 : TYPE(xc_rho_set_type), INTENT(INOUT) :: rho_set
1326 : TYPE(cp_1d_r_p_type), DIMENSION(:) :: ri_dcoeff
1327 : INTEGER, INTENT(IN) :: source_iat, source_ikind, target_iat, &
1328 : target_ikind, sr, er
1329 : LOGICAL, INTENT(IN) :: do_gga
1330 : TYPE(xas_atom_env_type), POINTER :: xas_atom_env
1331 : TYPE(qs_environment_type), POINTER :: qs_env
1332 :
1333 : CHARACTER(len=*), PARAMETER :: routineN = 'put_density_on_other_grid'
1334 :
1335 : INTEGER :: dir, handle, ia, ico, ipgf, ir, iset, &
1336 : isgf, lx, ly, lz, n, na, nr, nset, &
1337 : nspins, sgfi, start
1338 28 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nsgf_set
1339 28 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf
1340 : REAL(dp) :: rmom
1341 28 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: sgf
1342 28 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: co, dsgf, pos
1343 28 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dco
1344 : REAL(dp), DIMENSION(3) :: rs, rst, rt
1345 28 : REAL(dp), DIMENSION(:, :), POINTER :: ri_sphi, zet
1346 28 : REAL(dp), DIMENSION(:, :, :), POINTER :: rhoa, rhob
1347 : TYPE(cell_type), POINTER :: cell
1348 : TYPE(grid_atom_type), POINTER :: grid_atom
1349 : TYPE(gto_basis_set_type), POINTER :: ri_basis
1350 : TYPE(harmonics_atom_type), POINTER :: harmonics
1351 28 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1352 28 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1353 :
1354 28 : NULLIFY (qs_kind_set, ri_basis, lmax, npgf, nsgf_set, lmin, zet, first_sgf, grid_atom)
1355 28 : NULLIFY (harmonics, rhoa, rhob, particle_set, cell, ri_sphi)
1356 :
1357 : !Same logic as the put_density_on_own_grid routine. Loop over orbitals, put them on the grid,
1358 : !contract into sgf and daxpy with coeff. Notable difference: use cartesian orbitals instead of
1359 : !spherical, since the center of the grid is not the origin and thus, spherical symmetry can't
1360 : !be exploited so well
1361 :
1362 28 : CALL timeset(routineN, handle)
1363 :
1364 : !Generalities
1365 28 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, cell=cell)
1366 : !want basis of the source atom
1367 28 : CALL get_qs_kind(qs_kind_set(source_ikind), basis_set=ri_basis, basis_type="RI_XAS")
1368 : CALL get_gto_basis_set(ri_basis, lmax=lmax, npgf=npgf, zet=zet, nset=nset, nsgf_set=nsgf_set, &
1369 28 : first_sgf=first_sgf, lmin=lmin, sphi=ri_sphi)
1370 :
1371 : ! Want the grid and harmonics of the target atom
1372 28 : grid_atom => xas_atom_env%grid_atom_set(target_ikind)%grid_atom
1373 28 : harmonics => xas_atom_env%harmonics_atom_set(target_ikind)%harmonics_atom
1374 28 : na = grid_atom%ng_sphere
1375 28 : nr = er - sr + 1
1376 28 : n = na*nr
1377 28 : nspins = xas_atom_env%nspins
1378 :
1379 : ! Point to the rho_set densities
1380 28 : rhoa => rho_set%rhoa
1381 28 : rhob => rho_set%rhob
1382 :
1383 : ! Need the source-target position vector
1384 28 : rs = pbc(particle_set(source_iat)%r, cell)
1385 28 : rt = pbc(particle_set(target_iat)%r, cell)
1386 28 : rst = pbc(rs, rt, cell)
1387 :
1388 : ! Precompute the positions on the target grid
1389 140 : ALLOCATE (pos(na, sr:er, 4))
1390 : !$OMP PARALLEL DO COLLAPSE(2) SCHEDULE(STATIC) DEFAULT(NONE), &
1391 : !$OMP SHARED(na,sr,er,pos,harmonics,grid_atom,rst), &
1392 28 : !$OMP PRIVATE(ia,ir)
1393 : DO ir = sr, er
1394 : DO ia = 1, na
1395 : pos(ia, ir, 1:3) = harmonics%a(:, ia)*grid_atom%rad(ir) + rst
1396 : pos(ia, ir, 4) = pos(ia, ir, 1)**2 + pos(ia, ir, 2)**2 + pos(ia, ir, 3)**2
1397 : END DO
1398 : END DO
1399 : !$OMP END PARALLEL DO
1400 :
1401 : ! Loop over the cartesian gaussian functions and evaluate them
1402 84 : DO iset = 1, nset
1403 :
1404 : !allocate space to store the cartesian orbtial on the grid
1405 280 : ALLOCATE (co(na, sr:er, npgf(iset)*ncoset(lmax(iset))))
1406 336 : IF (do_gga) ALLOCATE (dco(na, sr:er, 3, npgf(iset)*ncoset(lmax(iset))))
1407 :
1408 : !$OMP PARALLEL DEFAULT(NONE), &
1409 : !$OMP SHARED(co,npgf,ncoset,lmax,lmin,indco,pos,zet,iset,na,sr,er,do_gga,dco), &
1410 56 : !$OMP PRIVATE(ipgf,start,ico,lx,ly,lz,ia,ir,rmom)
1411 :
1412 : !$OMP DO COLLAPSE(2) SCHEDULE(STATIC)
1413 : DO ir = sr, er
1414 : DO ia = 1, na
1415 : co(ia, ir, :) = 0.0_dp
1416 : IF (do_gga) THEN
1417 : dco(ia, ir, :, :) = 0.0_dp
1418 : END IF
1419 : END DO
1420 : END DO
1421 : !$OMP END DO NOWAIT
1422 :
1423 : DO ipgf = 1, npgf(iset)
1424 : start = (ipgf - 1)*ncoset(lmax(iset))
1425 :
1426 : !loop over the cartesian orbitals
1427 : DO ico = ncoset(lmin(iset) - 1) + 1, ncoset(lmax(iset))
1428 : lx = indco(1, ico)
1429 : ly = indco(2, ico)
1430 : lz = indco(3, ico)
1431 :
1432 : ! compute g = x**lx * y**ly * z**lz * exp(-zet * r**2)
1433 : !$OMP DO COLLAPSE(2) SCHEDULE(STATIC)
1434 : DO ir = sr, er
1435 : DO ia = 1, na
1436 : rmom = EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
1437 : IF (lx /= 0) rmom = rmom*pos(ia, ir, 1)**lx
1438 : IF (ly /= 0) rmom = rmom*pos(ia, ir, 2)**ly
1439 : IF (lz /= 0) rmom = rmom*pos(ia, ir, 3)**lz
1440 : co(ia, ir, start + ico) = rmom
1441 : !co(ia, ir, start + ico) = pos(ia, ir, 1)**lx*pos(ia, ir, 2)**ly*pos(ia, ir, 3)**lz &
1442 : ! *EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
1443 : END DO
1444 : END DO
1445 : !$OMP END DO NOWAIT
1446 :
1447 : IF (do_gga) THEN
1448 : !the gradient: dg_x = lx*x**(lx-1) * y**ly * z**lz * exp(-zet * r**2)
1449 : ! -2*zet* x**(lx+1) * y**ly * z**lz * exp(-zet * r**2)
1450 : ! = (lx*x**(lx-1) - 2*zet*x**(lx+1)) * y**ly * z**lz * exp(-zet * r**2)
1451 :
1452 : !x direction, special case if lx == 0
1453 : IF (lx == 0) THEN
1454 : !$OMP DO COLLAPSE(2) SCHEDULE(STATIC)
1455 : DO ir = sr, er
1456 : DO ia = 1, na
1457 : rmom = -2.0_dp*pos(ia, ir, 1)*zet(ipgf, iset)*EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
1458 : IF (ly /= 0) rmom = rmom*pos(ia, ir, 2)**ly
1459 : IF (lz /= 0) rmom = rmom*pos(ia, ir, 3)**lz
1460 : dco(ia, ir, 1, start + ico) = rmom
1461 : ! dco(ia, ir, 1, start + ico) = -2.0_dp*pos(ia, ir, 1)*zet(ipgf, iset) &
1462 : ! *pos(ia, ir, 2)**ly*pos(ia, ir, 3)**lz &
1463 : ! *EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
1464 : END DO
1465 : END DO
1466 : !$OMP END DO NOWAIT
1467 : ELSE
1468 : !$OMP DO COLLAPSE(2) SCHEDULE(STATIC)
1469 : DO ir = sr, er
1470 : DO ia = 1, na
1471 : IF (lx /= 1) THEN
1472 : rmom = (lx*pos(ia, ir, 1)**(lx - 1) - 2.0_dp*pos(ia, ir, 1)**(lx + 1)* &
1473 : zet(ipgf, iset))*EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
1474 : ELSE
1475 : rmom = (1.0_dp - 2.0_dp*pos(ia, ir, 1)**2*zet(ipgf, iset))* &
1476 : EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
1477 : END IF
1478 : IF (ly /= 0) rmom = rmom*pos(ia, ir, 2)**ly
1479 : IF (lz /= 0) rmom = rmom*pos(ia, ir, 3)**lz
1480 : dco(ia, ir, 1, start + ico) = rmom
1481 : ! dco(ia, ir, 1, start + ico) = (lx*pos(ia, ir, 1)**(lx - 1) &
1482 : ! - 2.0_dp*pos(ia, ir, 1)**(lx + 1)*zet(ipgf, iset)) &
1483 : ! *pos(ia, ir, 2)**ly*pos(ia, ir, 3)**lz &
1484 : ! *EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
1485 : END DO
1486 : END DO
1487 : !$OMP END DO NOWAIT
1488 : END IF !lx == 0
1489 :
1490 : !y direction, special case if ly == 0
1491 : IF (ly == 0) THEN
1492 : !$OMP DO COLLAPSE(2) SCHEDULE(STATIC)
1493 : DO ir = sr, er
1494 : DO ia = 1, na
1495 : rmom = -2.0_dp*pos(ia, ir, 2)*zet(ipgf, iset)*EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
1496 : IF (lx /= 0) rmom = rmom*pos(ia, ir, 1)**lx
1497 : IF (lz /= 0) rmom = rmom*pos(ia, ir, 3)**lz
1498 : dco(ia, ir, 2, start + ico) = rmom
1499 : ! dco(ia, ir, 2, start + ico) = -2.0_dp*pos(ia, ir, 2)*zet(ipgf, iset) &
1500 : ! *pos(ia, ir, 1)**lx*pos(ia, ir, 3)**lz &
1501 : ! *EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
1502 : END DO
1503 : END DO
1504 : !$OMP END DO NOWAIT
1505 : ELSE
1506 : !$OMP DO COLLAPSE(2) SCHEDULE(STATIC)
1507 : DO ir = sr, er
1508 : DO ia = 1, na
1509 : IF (ly /= 1) THEN
1510 : rmom = (ly*pos(ia, ir, 2)**(ly - 1) - 2.0_dp*pos(ia, ir, 2)**(ly + 1)*zet(ipgf, iset)) &
1511 : *EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
1512 : ELSE
1513 : rmom = (1.0_dp - 2.0_dp*pos(ia, ir, 2)**2*zet(ipgf, iset)) &
1514 : *EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
1515 : END IF
1516 : IF (lx /= 0) rmom = rmom*pos(ia, ir, 1)**lx
1517 : IF (lz /= 0) rmom = rmom*pos(ia, ir, 3)**lz
1518 : dco(ia, ir, 2, start + ico) = rmom
1519 : ! dco(ia, ir, 2, start + ico) = (ly*pos(ia, ir, 2)**(ly - 1) &
1520 : ! - 2.0_dp*pos(ia, ir, 2)**(ly + 1)*zet(ipgf, iset)) &
1521 : ! *pos(ia, ir, 1)**lx*pos(ia, ir, 3)**lz &
1522 : ! *EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
1523 : END DO
1524 : END DO
1525 : !$OMP END DO NOWAIT
1526 : END IF !ly == 0
1527 :
1528 : !z direction, special case if lz == 0
1529 : IF (lz == 0) THEN
1530 : !$OMP DO COLLAPSE(2) SCHEDULE(STATIC)
1531 : DO ir = sr, er
1532 : DO ia = 1, na
1533 : rmom = -2.0_dp*pos(ia, ir, 3)*zet(ipgf, iset)*EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
1534 : IF (lx /= 0) rmom = rmom*pos(ia, ir, 1)**lx
1535 : IF (ly /= 0) rmom = rmom*pos(ia, ir, 2)**ly
1536 : dco(ia, ir, 3, start + ico) = rmom
1537 : ! dco(ia, ir, 3, start + ico) = -2.0_dp*pos(ia, ir, 3)*zet(ipgf, iset) &
1538 : ! *pos(ia, ir, 1)**lx*pos(ia, ir, 2)**ly &
1539 : ! *EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
1540 : END DO
1541 : END DO
1542 : !$OMP END DO NOWAIT
1543 : ELSE
1544 : !$OMP DO COLLAPSE(2) SCHEDULE(STATIC)
1545 : DO ir = sr, er
1546 : DO ia = 1, na
1547 : IF (lz /= 1) THEN
1548 : rmom = (lz*pos(ia, ir, 3)**(lz - 1) - 2.0_dp*pos(ia, ir, 3)**(lz + 1)* &
1549 : zet(ipgf, iset))*EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
1550 : ELSE
1551 : rmom = (1.0_dp - 2.0_dp*pos(ia, ir, 3)**2*zet(ipgf, iset))* &
1552 : EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
1553 : END IF
1554 : IF (lx /= 0) rmom = rmom*pos(ia, ir, 1)**lx
1555 : IF (ly /= 0) rmom = rmom*pos(ia, ir, 2)**ly
1556 : dco(ia, ir, 3, start + ico) = rmom
1557 : ! dco(ia, ir, 3, start + ico) = (lz*pos(ia, ir, 3)**(lz - 1) &
1558 : ! - 2.0_dp*pos(ia, ir, 3)**(lz + 1)*zet(ipgf, iset)) &
1559 : ! *pos(ia, ir, 1)**lx*pos(ia, ir, 2)**ly &
1560 : ! *EXP(-zet(ipgf, iset)*pos(ia, ir, 4))
1561 : END DO
1562 : END DO
1563 : !$OMP END DO NOWAIT
1564 : END IF !lz == 0
1565 :
1566 : END IF !gga
1567 :
1568 : END DO !ico
1569 : END DO !ipgf
1570 :
1571 : !$OMP END PARALLEL
1572 :
1573 : !contract the co into sgf
1574 224 : ALLOCATE (sgf(na, sr:er))
1575 280 : IF (do_gga) ALLOCATE (dsgf(na, sr:er, 3))
1576 56 : sgfi = first_sgf(1, iset) - 1
1577 :
1578 706 : DO isgf = 1, nsgf_set(iset)
1579 25199993 : sgf = 0.0_dp
1580 75601279 : IF (do_gga) dsgf = 0.0_dp
1581 :
1582 4048 : DO ipgf = 1, npgf(iset)
1583 3398 : start = (ipgf - 1)*ncoset(lmax(iset))
1584 27330 : DO ico = ncoset(lmin(iset) - 1) + 1, ncoset(lmax(iset))
1585 26680 : CALL daxpy(n, ri_sphi(start + ico, sgfi + isgf), co(:, sr:er, start + ico), 1, sgf(:, sr:er), 1)
1586 : END DO !ico
1587 : END DO !ipgf
1588 :
1589 : !add the density to the grid
1590 650 : CALL daxpy(n, ri_dcoeff(1)%array(sgfi + isgf), sgf(:, sr:er), 1, rhoa(:, sr:er, 1), 1)
1591 :
1592 650 : IF (nspins == 2) THEN
1593 414 : CALL daxpy(n, ri_dcoeff(2)%array(sgfi + isgf), sgf(:, sr:er), 1, rhob(:, sr:er, 1), 1)
1594 : END IF
1595 :
1596 : !deal with the gradient
1597 706 : IF (do_gga) THEN
1598 :
1599 4048 : DO ipgf = 1, npgf(iset)
1600 3398 : start = (ipgf - 1)*ncoset(lmax(iset))
1601 27330 : DO ico = ncoset(lmin(iset) - 1) + 1, ncoset(lmax(iset))
1602 96526 : DO dir = 1, 3
1603 : CALL daxpy(n, ri_sphi(start + ico, sgfi + isgf), dco(:, sr:er, dir, start + ico), &
1604 93128 : 1, dsgf(:, sr:er, dir), 1)
1605 : END DO
1606 : END DO !ico
1607 : END DO !ipgf
1608 :
1609 2600 : DO dir = 1, 3
1610 : CALL daxpy(n, ri_dcoeff(1)%array(sgfi + isgf), dsgf(:, sr:er, dir), 1, &
1611 1950 : rho_set%drhoa(dir)%array(:, sr:er, 1), 1)
1612 :
1613 2600 : IF (nspins == 2) THEN
1614 : CALL daxpy(n, ri_dcoeff(2)%array(sgfi + isgf), dsgf(:, sr:er, dir), 1, &
1615 1242 : rho_set%drhob(dir)%array(:, sr:er, 1), 1)
1616 : END IF
1617 : END DO
1618 : END IF !do_gga
1619 :
1620 : END DO !isgf
1621 :
1622 56 : DEALLOCATE (co, sgf)
1623 84 : IF (do_gga) DEALLOCATE (dco, dsgf)
1624 : END DO !iset
1625 :
1626 : !Treat spin-restricted case (copy alpha into beta)
1627 28 : IF (nspins == 1) THEN
1628 10 : CALL dcopy(n, rhoa(:, sr:er, 1), 1, rhob(:, sr:er, 1), 1)
1629 :
1630 10 : IF (do_gga) THEN
1631 40 : DO dir = 1, 3
1632 40 : CALL dcopy(n, rho_set%drhoa(dir)%array(:, sr:er, 1), 1, rho_set%drhob(dir)%array(:, sr:er, 1), 1)
1633 : END DO
1634 : END IF
1635 : END IF
1636 :
1637 28 : CALL timestop(handle)
1638 :
1639 84 : END SUBROUTINE put_density_on_other_grid
1640 :
1641 : ! **************************************************************************************************
1642 : !> \brief Computes the norm of the density gradient on the atomic grid
1643 : !> \param rho_set ...
1644 : !> \param atom_kind ...
1645 : !> \param xas_atom_env ...
1646 : !> \note GGA is assumed
1647 : ! **************************************************************************************************
1648 28 : SUBROUTINE compute_norm_drho(rho_set, atom_kind, xas_atom_env)
1649 :
1650 : TYPE(xc_rho_set_type), INTENT(INOUT) :: rho_set
1651 : INTEGER, INTENT(IN) :: atom_kind
1652 : TYPE(xas_atom_env_type), POINTER :: xas_atom_env
1653 :
1654 : INTEGER :: dir, ia, ir, n, na, nr, nspins
1655 :
1656 28 : na = xas_atom_env%grid_atom_set(atom_kind)%grid_atom%ng_sphere
1657 28 : nr = xas_atom_env%grid_atom_set(atom_kind)%grid_atom%nr
1658 28 : n = na*nr
1659 28 : nspins = xas_atom_env%nspins
1660 :
1661 903248 : rho_set%norm_drhoa = 0.0_dp
1662 903248 : rho_set%norm_drhob = 0.0_dp
1663 903248 : rho_set%norm_drho = 0.0_dp
1664 :
1665 112 : DO dir = 1, 3
1666 12856 : DO ir = 1, nr
1667 2709660 : DO ia = 1, na
1668 : rho_set%norm_drhoa(ia, ir, 1) = rho_set%norm_drhoa(ia, ir, 1) &
1669 2709576 : + rho_set%drhoa(dir)%array(ia, ir, 1)**2
1670 : END DO !ia
1671 : END DO !ir
1672 : END DO !dir
1673 903248 : rho_set%norm_drhoa = SQRT(rho_set%norm_drhoa)
1674 :
1675 28 : IF (nspins == 1) THEN
1676 : !spin-restricted
1677 22 : CALL dcopy(n, rho_set%norm_drhoa(:, :, 1), 1, rho_set%norm_drhob(:, :, 1), 1)
1678 : ELSE
1679 24 : DO dir = 1, 3
1680 4398 : DO ir = 1, nr
1681 1325340 : DO ia = 1, na
1682 : rho_set%norm_drhob(ia, ir, 1) = rho_set%norm_drhob(ia, ir, 1) &
1683 1325322 : + rho_set%drhob(dir)%array(ia, ir, 1)**2
1684 : END DO
1685 : END DO
1686 : END DO
1687 441786 : rho_set%norm_drhob = SQRT(rho_set%norm_drhob)
1688 : END IF
1689 :
1690 112 : DO dir = 1, 3
1691 12856 : DO ir = 1, nr
1692 2709660 : DO ia = 1, na
1693 : rho_set%norm_drho(ia, ir, 1) = rho_set%norm_drho(ia, ir, 1) + &
1694 : (rho_set%drhoa(dir)%array(ia, ir, 1) + &
1695 2709576 : rho_set%drhob(dir)%array(ia, ir, 1))**2
1696 : END DO
1697 : END DO
1698 : END DO
1699 903248 : rho_set%norm_drho = SQRT(rho_set%norm_drho)
1700 :
1701 28 : END SUBROUTINE compute_norm_drho
1702 :
1703 : ! **************************************************************************************************
1704 : !> \brief Precomputes the spherical orbitals of the RI basis on the excited atom grids
1705 : !> \param do_gga whether the gradient needs to be computed for GGA or not
1706 : !> \param batch_info the parallelization info to complete with so distribution info
1707 : !> \param xas_atom_env ...
1708 : !> \param qs_env ...
1709 : !> \note the functions are split in a purely angular part of size na and a purely radial part of
1710 : !> size nr. The full function on the grid can simply be obtained with dgemm and we save space
1711 : ! **************************************************************************************************
1712 52 : SUBROUTINE precompute_so_dso(do_gga, batch_info, xas_atom_env, qs_env)
1713 :
1714 : LOGICAL, INTENT(IN) :: do_gga
1715 : TYPE(batch_info_type) :: batch_info
1716 : TYPE(xas_atom_env_type), POINTER :: xas_atom_env
1717 : TYPE(qs_environment_type), POINTER :: qs_env
1718 :
1719 : CHARACTER(len=*), PARAMETER :: routineN = 'precompute_so_dso'
1720 :
1721 : INTEGER :: bo(2), dir, handle, ikind, ipgf, iset, &
1722 : iso, iso_proc, l, maxso, n, na, nkind, &
1723 : nr, nset, nso_proc, nsotot, starti
1724 52 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nsgf_set
1725 52 : INTEGER, DIMENSION(:, :), POINTER :: so_proc_info
1726 52 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: rexp
1727 52 : REAL(dp), DIMENSION(:, :), POINTER :: dgr1, dgr2, ga, gr, slm, zet
1728 52 : REAL(dp), DIMENSION(:, :, :), POINTER :: dga1, dga2, dslm_dxyz
1729 : TYPE(grid_atom_type), POINTER :: grid_atom
1730 : TYPE(gto_basis_set_type), POINTER :: ri_basis
1731 : TYPE(harmonics_atom_type), POINTER :: harmonics
1732 : TYPE(mp_para_env_type), POINTER :: para_env
1733 52 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1734 :
1735 52 : NULLIFY (qs_kind_set, harmonics, grid_atom, slm, dslm_dxyz, ri_basis, lmax, lmin, npgf)
1736 52 : NULLIFY (nsgf_set, zet, para_env, so_proc_info)
1737 :
1738 52 : CALL timeset(routineN, handle)
1739 :
1740 52 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, para_env=para_env)
1741 52 : nkind = SIZE(qs_kind_set)
1742 :
1743 240 : ALLOCATE (batch_info%so_proc_info(nkind))
1744 156 : ALLOCATE (batch_info%nso_proc(nkind))
1745 156 : ALLOCATE (batch_info%so_bo(2, nkind))
1746 :
1747 136 : DO ikind = 1, nkind
1748 :
1749 84 : NULLIFY (xas_atom_env%ga(ikind)%array)
1750 84 : NULLIFY (xas_atom_env%gr(ikind)%array)
1751 84 : NULLIFY (xas_atom_env%dga1(ikind)%array)
1752 84 : NULLIFY (xas_atom_env%dga2(ikind)%array)
1753 84 : NULLIFY (xas_atom_env%dgr1(ikind)%array)
1754 84 : NULLIFY (xas_atom_env%dgr2(ikind)%array)
1755 :
1756 84 : NULLIFY (batch_info%so_proc_info(ikind)%array)
1757 :
1758 118 : IF (.NOT. ANY(xas_atom_env%excited_kinds == ikind)) CYCLE
1759 :
1760 : !grid info
1761 58 : harmonics => xas_atom_env%harmonics_atom_set(ikind)%harmonics_atom
1762 58 : grid_atom => xas_atom_env%grid_atom_set(ikind)%grid_atom
1763 :
1764 58 : na = grid_atom%ng_sphere
1765 58 : nr = grid_atom%nr
1766 58 : n = na*nr
1767 :
1768 58 : slm => harmonics%slm
1769 58 : dslm_dxyz => harmonics%dslm_dxyz
1770 :
1771 : !basis info
1772 58 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis, basis_type="RI_XAS")
1773 : CALL get_gto_basis_set(ri_basis, lmax=lmax, npgf=npgf, zet=zet, nset=nset, &
1774 58 : nsgf_set=nsgf_set, lmin=lmin, maxso=maxso)
1775 58 : nsotot = maxso*nset
1776 :
1777 : !we split all so among the processors of the batch
1778 58 : bo = get_limit(nsotot, batch_info%batch_size, batch_info%ipe)
1779 58 : nso_proc = bo(2) - bo(1) + 1
1780 174 : batch_info%so_bo(:, ikind) = bo
1781 58 : batch_info%nso_proc(ikind) = nso_proc
1782 :
1783 : !store info about the so's set, pgf and index
1784 174 : ALLOCATE (batch_info%so_proc_info(ikind)%array(3, nso_proc))
1785 58 : so_proc_info => batch_info%so_proc_info(ikind)%array
1786 26290 : so_proc_info = -1 !default is -1 => set so value to zero
1787 250 : DO iset = 1, nset
1788 1188 : DO ipgf = 1, npgf(iset)
1789 938 : starti = (iset - 1)*maxso + (ipgf - 1)*nsoset(lmax(iset))
1790 6820 : DO iso = nsoset(lmin(iset) - 1) + 1, nsoset(lmax(iset))
1791 :
1792 : !only consider so that are on this proc
1793 5690 : IF (starti + iso < bo(1) .OR. starti + iso > bo(2)) CYCLE
1794 3462 : iso_proc = starti + iso - bo(1) + 1
1795 3462 : so_proc_info(1, iso_proc) = iset
1796 3462 : so_proc_info(2, iso_proc) = ipgf
1797 6628 : so_proc_info(3, iso_proc) = iso
1798 :
1799 : END DO
1800 : END DO
1801 : END DO
1802 :
1803 : !Put the gaussians and their gradient as purely angular or radial arrays
1804 232 : ALLOCATE (xas_atom_env%ga(ikind)%array(na, nso_proc))
1805 232 : ALLOCATE (xas_atom_env%gr(ikind)%array(nr, nso_proc))
1806 2235148 : xas_atom_env%ga(ikind)%array = 0.0_dp; xas_atom_env%gr(ikind)%array = 0.0_dp
1807 58 : IF (do_gga) THEN
1808 170 : ALLOCATE (xas_atom_env%dga1(ikind)%array(na, nso_proc, 3))
1809 102 : ALLOCATE (xas_atom_env%dgr1(ikind)%array(nr, nso_proc))
1810 102 : ALLOCATE (xas_atom_env%dga2(ikind)%array(na, nso_proc, 3))
1811 102 : ALLOCATE (xas_atom_env%dgr2(ikind)%array(nr, nso_proc))
1812 2586058 : xas_atom_env%dga1(ikind)%array = 0.0_dp; xas_atom_env%dgr1(ikind)%array = 0.0_dp
1813 2586058 : xas_atom_env%dga2(ikind)%array = 0.0_dp; xas_atom_env%dgr2(ikind)%array = 0.0_dp
1814 : END IF
1815 :
1816 58 : ga => xas_atom_env%ga(ikind)%array
1817 58 : gr => xas_atom_env%gr(ikind)%array
1818 58 : dga1 => xas_atom_env%dga1(ikind)%array
1819 58 : dga2 => xas_atom_env%dga2(ikind)%array
1820 58 : dgr1 => xas_atom_env%dgr1(ikind)%array
1821 58 : dgr2 => xas_atom_env%dgr2(ikind)%array
1822 :
1823 174 : ALLOCATE (rexp(nr))
1824 :
1825 6616 : DO iso_proc = 1, nso_proc
1826 6558 : iset = so_proc_info(1, iso_proc)
1827 6558 : ipgf = so_proc_info(2, iso_proc)
1828 6558 : iso = so_proc_info(3, iso_proc)
1829 6558 : IF (iso < 0) CYCLE
1830 :
1831 3462 : l = indso(1, iso)
1832 :
1833 : !The gaussian is g = r^l * Ylm * exp(-a*r^2)
1834 :
1835 : !radial part of the gaussian
1836 494499 : rexp(1:nr) = EXP(-zet(ipgf, iset)*grid_atom%rad2(1:nr))
1837 985536 : gr(1:nr, iso_proc) = grid_atom%rad(1:nr)**l*rexp(1:nr)
1838 :
1839 : !angular part of the gaussian
1840 1232334 : ga(1:na, iso_proc) = slm(1:na, iso)
1841 :
1842 : !For the gradient, devide in 2 parts: dg/dx = d/dx(r^l * Ylm) * exp(-a*r^2)
1843 : ! + r^l * Ylm * d/dx(exp(-a*r^2))
1844 : !Note: we make this choice of separation because of cartesian coordinates, where
1845 : ! g = x^lx * y^ly * z^lz * exp(-a*r^2) and r^(l-1)*dslm_dxyz = d/dx(r^l * Ylm)
1846 :
1847 3520 : IF (do_gga) THEN
1848 : !radial part of the gradient => same in all three direction
1849 679285 : dgr1(1:nr, iso_proc) = grid_atom%rad(1:nr)**(l - 1)*rexp(1:nr)
1850 679285 : dgr2(1:nr, iso_proc) = -2.0_dp*zet(ipgf, iset)*grid_atom%rad(1:nr)**(l + 1)*rexp(1:nr)
1851 :
1852 : !angular part of the gradient
1853 9316 : DO dir = 1, 3
1854 2469399 : dga1(1:na, iso_proc, dir) = dslm_dxyz(dir, 1:na, iso)
1855 2471728 : dga2(1:na, iso_proc, dir) = harmonics%a(dir, 1:na)*slm(1:na, iso)
1856 : END DO
1857 : END IF
1858 :
1859 : END DO !iso_proc
1860 :
1861 194 : DEALLOCATE (rexp)
1862 : END DO !ikind
1863 :
1864 52 : CALL timestop(handle)
1865 :
1866 104 : END SUBROUTINE precompute_so_dso
1867 :
1868 : ! **************************************************************************************************
1869 : !> \brief Integrate the xc kernel as a function of r on the atomic grids for the RI_XAS basis
1870 : !> \param int_fxc the global array containing the (P|fxc|Q) integrals, for all spin configurations
1871 : !> \param xas_atom_env ...
1872 : !> \param xas_tdp_control ...
1873 : !> \param qs_env ...
1874 : !> \note Note that if closed-shell, alpha-alpha term and beta-beta terms are the same
1875 : !> Store the (P|fxc|Q) integrals on the processor they were computed on
1876 : !> int_fxc(1)%matrix is alpha-alpha, 2: alpha-beta, 3: beta-beta
1877 : ! **************************************************************************************************
1878 52 : SUBROUTINE integrate_fxc_atoms(int_fxc, xas_atom_env, xas_tdp_control, qs_env)
1879 :
1880 : TYPE(cp_2d_r_p_type), DIMENSION(:, :), POINTER :: int_fxc
1881 : TYPE(xas_atom_env_type), POINTER :: xas_atom_env
1882 : TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
1883 : TYPE(qs_environment_type), POINTER :: qs_env
1884 :
1885 : CHARACTER(len=*), PARAMETER :: routineN = 'integrate_fxc_atoms'
1886 :
1887 : INTEGER :: batch_size, dir, er, ex_bo(2), handle, i, iatom, ibatch, iex, ikind, inb, ipe, &
1888 : mepos, na, natom, nb, nb_bo(2), nbatch, nbk, nex_atom, nr, num_pe, sr
1889 : INTEGER, DIMENSION(2, 3) :: bounds
1890 52 : INTEGER, DIMENSION(:), POINTER :: exat_neighbors
1891 : LOGICAL :: do_gga, do_sc, do_sf
1892 52 : TYPE(batch_info_type) :: batch_info
1893 52 : TYPE(cp_1d_r_p_type), DIMENSION(:, :), POINTER :: ri_dcoeff
1894 : TYPE(dft_control_type), POINTER :: dft_control
1895 : TYPE(mp_para_env_type), POINTER :: para_env
1896 52 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1897 52 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1898 : TYPE(section_vals_type), POINTER :: input, xc_functionals
1899 : TYPE(xc_derivative_set_type) :: deriv_set
1900 : TYPE(xc_rho_cflags_type) :: needs
1901 : TYPE(xc_rho_set_type) :: rho_set
1902 :
1903 52 : NULLIFY (particle_set, qs_kind_set, dft_control, para_env, exat_neighbors)
1904 52 : NULLIFY (input, xc_functionals)
1905 :
1906 52 : CALL timeset(routineN, handle)
1907 :
1908 : ! Initialize
1909 : CALL get_qs_env(qs_env, particle_set=particle_set, qs_kind_set=qs_kind_set, natom=natom, &
1910 52 : dft_control=dft_control, input=input, para_env=para_env)
1911 2004 : ALLOCATE (int_fxc(natom, 4))
1912 462 : DO iatom = 1, natom
1913 2102 : DO i = 1, 4
1914 2050 : NULLIFY (int_fxc(iatom, i)%array)
1915 : END DO
1916 : END DO
1917 52 : nex_atom = SIZE(xas_atom_env%excited_atoms)
1918 : !spin conserving in the general sense here
1919 52 : do_sc = xas_tdp_control%do_spin_cons .OR. xas_tdp_control%do_singlet .OR. xas_tdp_control%do_triplet
1920 52 : do_sf = xas_tdp_control%do_spin_flip
1921 :
1922 : ! Get some info on the functionals
1923 52 : IF (qs_env%do_rixs) THEN
1924 14 : xc_functionals => section_vals_get_subs_vals(input, "PROPERTIES%RIXS%XAS_TDP%KERNEL%XC_FUNCTIONAL")
1925 : ELSE
1926 38 : xc_functionals => section_vals_get_subs_vals(input, "DFT%XAS_TDP%KERNEL%XC_FUNCTIONAL")
1927 : END IF
1928 : ! ask for lsd in any case
1929 52 : needs = xc_functionals_get_needs(xc_functionals, lsd=.TRUE., calc_potential=.TRUE.)
1930 52 : do_gga = needs%drho_spin !because either LDA or GGA, and the former does not need gradient
1931 :
1932 : ! Distribute the excited atoms over batches of processors
1933 : ! Then, the spherical orbital of the RI basis are distributed among the procs of the batch, making
1934 : ! the GGA integration very efficient
1935 52 : num_pe = para_env%num_pe
1936 52 : mepos = para_env%mepos
1937 :
1938 : !create a batch_info_type
1939 52 : CALL get_proc_batch_sizes(batch_size, nbatch, nex_atom, num_pe)
1940 :
1941 : !the batch index
1942 52 : ibatch = mepos/batch_size
1943 : !the proc index within the batch
1944 52 : ipe = MODULO(mepos, batch_size)
1945 :
1946 52 : batch_info%batch_size = batch_size
1947 52 : batch_info%nbatch = nbatch
1948 52 : batch_info%ibatch = ibatch
1949 52 : batch_info%ipe = ipe
1950 :
1951 : !create a subcommunicator for this batch
1952 52 : CALL batch_info%para_env%from_split(para_env, ibatch)
1953 :
1954 : ! Precompute the spherical orbital of the RI basis (and maybe their gradient) on the grids of the
1955 : ! excited atoms. Needed for the GGA integration and to actually put the density on the grid
1956 52 : CALL precompute_so_dso(do_gga, batch_info, xas_atom_env, qs_env)
1957 :
1958 : !distribute the excted atoms over the batches
1959 52 : ex_bo = get_limit(nex_atom, nbatch, ibatch)
1960 :
1961 : ! Looping over the excited atoms
1962 104 : DO iex = ex_bo(1), ex_bo(2)
1963 :
1964 52 : iatom = xas_atom_env%excited_atoms(iex)
1965 52 : ikind = particle_set(iatom)%atomic_kind%kind_number
1966 52 : exat_neighbors => xas_atom_env%exat_neighbors(iex)%array
1967 52 : ri_dcoeff => xas_atom_env%ri_dcoeff(:, :, iex)
1968 :
1969 : ! General grid/basis info
1970 52 : na = xas_atom_env%grid_atom_set(ikind)%grid_atom%ng_sphere
1971 52 : nr = xas_atom_env%grid_atom_set(ikind)%grid_atom%nr
1972 :
1973 : ! Creating a xc_rho_set to store the density and dset for the kernel
1974 520 : bounds(1:2, 1:3) = 1
1975 52 : bounds(2, 1) = na
1976 52 : bounds(2, 2) = nr
1977 :
1978 : CALL xc_rho_set_create(rho_set=rho_set, local_bounds=bounds, &
1979 : rho_cutoff=dft_control%qs_control%eps_rho_rspace, &
1980 52 : drho_cutoff=dft_control%qs_control%eps_rho_rspace)
1981 52 : CALL xc_dset_create(deriv_set, local_bounds=bounds)
1982 :
1983 : ! allocate internals of the rho_set
1984 52 : CALL xc_rho_set_atom_update(rho_set, needs, nspins=2, bo=bounds)
1985 :
1986 : ! Put the density, and possibly its gradient, on the grid (for this atom)
1987 : CALL put_density_on_atomic_grid(rho_set, ri_dcoeff(iatom, :), ikind, &
1988 52 : do_gga, batch_info, xas_atom_env, qs_env)
1989 :
1990 : ! Take the neighboring atom contributions to the density (and gradient)
1991 : ! distribute the grid among the procs (for best load balance)
1992 52 : nb_bo = get_limit(nr, batch_size, ipe)
1993 52 : sr = nb_bo(1); er = nb_bo(2)
1994 80 : DO inb = 1, SIZE(exat_neighbors)
1995 :
1996 28 : nb = exat_neighbors(inb)
1997 28 : nbk = particle_set(nb)%atomic_kind%kind_number
1998 : CALL put_density_on_other_grid(rho_set, ri_dcoeff(nb, :), nb, nbk, iatom, &
1999 80 : ikind, sr, er, do_gga, xas_atom_env, qs_env)
2000 :
2001 : END DO
2002 :
2003 : ! make sure contributions from different procs are summed up
2004 52 : CALL batch_info%para_env%sum(rho_set%rhoa)
2005 52 : CALL batch_info%para_env%sum(rho_set%rhob)
2006 52 : IF (do_gga) THEN
2007 112 : DO dir = 1, 3
2008 84 : CALL batch_info%para_env%sum(rho_set%drhoa(dir)%array)
2009 112 : CALL batch_info%para_env%sum(rho_set%drhob(dir)%array)
2010 : END DO
2011 : END IF
2012 :
2013 : ! In case of GGA, also need the norm of the density gradient
2014 52 : IF (do_gga) CALL compute_norm_drho(rho_set, ikind, xas_atom_env)
2015 :
2016 : ! Compute the required derivatives
2017 : CALL xc_functionals_eval(xc_functionals, lsd=.TRUE., rho_set=rho_set, deriv_set=deriv_set, &
2018 52 : deriv_order=2)
2019 :
2020 : !spin-conserving (LDA part)
2021 52 : IF (do_sc) THEN
2022 52 : CALL integrate_sc_fxc(int_fxc, iatom, ikind, deriv_set, xas_atom_env, qs_env)
2023 : END IF
2024 :
2025 : !spin-flip (LDA part)
2026 52 : IF (do_sf) THEN
2027 0 : CALL integrate_sf_fxc(int_fxc, iatom, ikind, rho_set, deriv_set, xas_atom_env, qs_env)
2028 : END IF
2029 :
2030 : !Gradient correction (note: spin-flip only keeps the lda part, aka ALDA0)
2031 52 : IF (do_gga .AND. do_sc) THEN
2032 : CALL integrate_gga_fxc(int_fxc, iatom, ikind, batch_info, rho_set, deriv_set, &
2033 28 : xas_atom_env, qs_env)
2034 : END IF
2035 :
2036 : ! Clean-up
2037 52 : CALL xc_dset_release(deriv_set)
2038 156 : CALL xc_rho_set_release(rho_set)
2039 : END DO !iex
2040 :
2041 52 : CALL release_batch_info(batch_info)
2042 :
2043 : !Not necessary to sync, but makes sure that any load inbalance is reported here
2044 52 : CALL para_env%sync()
2045 :
2046 52 : CALL timestop(handle)
2047 :
2048 1248 : END SUBROUTINE integrate_fxc_atoms
2049 :
2050 : ! **************************************************************************************************
2051 : !> \brief Integrate the gradient correction part of the xc kernel on the atomic grid
2052 : !> \param int_fxc the array containing the (P|fxc|Q) integrals
2053 : !> \param iatom the index of the current excited atom
2054 : !> \param ikind the index of the current excited kind
2055 : !> \param batch_info how the so are distributed over the processor batch
2056 : !> \param rho_set the variable contatinind the density and its gradient
2057 : !> \param deriv_set the functional derivatives
2058 : !> \param xas_atom_env ...
2059 : !> \param qs_env ...
2060 : !> \note Ignored in case of pure LDA, added on top of the LDA kernel in case of GGA
2061 : ! **************************************************************************************************
2062 28 : SUBROUTINE integrate_gga_fxc(int_fxc, iatom, ikind, batch_info, rho_set, deriv_set, &
2063 : xas_atom_env, qs_env)
2064 :
2065 : TYPE(cp_2d_r_p_type), DIMENSION(:, :), POINTER :: int_fxc
2066 : INTEGER, INTENT(IN) :: iatom, ikind
2067 : TYPE(batch_info_type) :: batch_info
2068 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
2069 : TYPE(xc_derivative_set_type), INTENT(INOUT) :: deriv_set
2070 : TYPE(xas_atom_env_type), POINTER :: xas_atom_env
2071 : TYPE(qs_environment_type), POINTER :: qs_env
2072 :
2073 : CHARACTER(len=*), PARAMETER :: routineN = 'integrate_gga_fxc'
2074 :
2075 : INTEGER :: bo(2), dir, handle, i, ia, ir, jpgf, &
2076 : jset, jso, l, maxso, na, nr, nset, &
2077 : nsgf, nsoi, nsotot, startj, ub
2078 28 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
2079 28 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: rexp
2080 28 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: int_sgf, res, so, work
2081 28 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: dso
2082 28 : REAL(dp), DIMENSION(:, :), POINTER :: dgr1, dgr2, ga, gr, ri_sphi_so, weight, &
2083 28 : zet
2084 28 : REAL(dp), DIMENSION(:, :, :), POINTER :: dga1, dga2
2085 28 : TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:) :: int_so, vxc
2086 : TYPE(cp_3d_r_p_type), ALLOCATABLE, DIMENSION(:) :: vxg
2087 : TYPE(grid_atom_type), POINTER :: grid_atom
2088 : TYPE(gto_basis_set_type), POINTER :: ri_basis
2089 : TYPE(harmonics_atom_type), POINTER :: harmonics
2090 : TYPE(mp_para_env_type), POINTER :: para_env
2091 28 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2092 :
2093 28 : NULLIFY (grid_atom, ri_basis, qs_kind_set, ga, gr, dgr1, dgr2, lmax, lmin, npgf)
2094 28 : NULLIFY (weight, ri_sphi_so, dga1, dga2, para_env, harmonics, zet)
2095 :
2096 : !Strategy: we need to compute <phi_i|fxc|phij>, most of existing application of the 2nd
2097 : ! functional derivative involve the response density, and the expression of the
2098 : ! integral int (fxc*n^1) is well known. We substitute the spherical orbital phi_j
2099 : ! in place of n^1 in the formula and thus perform the first integration. Then
2100 : ! we obtain something in the form int (fxc*phi_j) = Vxc - div (Vxg) that we can
2101 : ! put on the grid and treat like a potential. The second integration is done by
2102 : ! using the divergence theorem and numerical integration:
2103 : ! <phi_i|fxc|phi_j> = int phi_i*(Vxc - div(Vxg)) = int phi_i*Vxc + grad(phi_i).Vxg
2104 : ! Note the sign change and the dot product.
2105 :
2106 28 : CALL timeset(routineN, handle)
2107 :
2108 : !If closed shell, only compute f_aa and f_ab (ub = 2)
2109 28 : ub = 2
2110 28 : IF (xas_atom_env%nspins == 2) ub = 3
2111 :
2112 : !Get the necessary grid info
2113 28 : harmonics => xas_atom_env%harmonics_atom_set(ikind)%harmonics_atom
2114 28 : grid_atom => xas_atom_env%grid_atom_set(ikind)%grid_atom
2115 28 : na = grid_atom%ng_sphere
2116 28 : nr = grid_atom%nr
2117 28 : weight => grid_atom%weight
2118 :
2119 : !get the ri_basis info
2120 28 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, para_env=para_env)
2121 28 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis, basis_type="RI_XAS")
2122 :
2123 : CALL get_gto_basis_set(gto_basis_set=ri_basis, lmax=lmax, lmin=lmin, nsgf=nsgf, &
2124 28 : maxso=maxso, npgf=npgf, nset=nset, zet=zet)
2125 28 : nsotot = nset*maxso
2126 :
2127 : !Point to the precomputed so
2128 28 : ga => xas_atom_env%ga(ikind)%array
2129 28 : gr => xas_atom_env%gr(ikind)%array
2130 28 : dgr1 => xas_atom_env%dgr1(ikind)%array
2131 28 : dgr2 => xas_atom_env%dgr2(ikind)%array
2132 28 : dga1 => xas_atom_env%dga1(ikind)%array
2133 28 : dga2 => xas_atom_env%dga2(ikind)%array
2134 :
2135 : !Before integration, wanna pre-divide all relevant derivastives by the nrom of the gradient
2136 28 : CALL divide_by_norm_drho(deriv_set, rho_set, lsd=.TRUE.)
2137 :
2138 : !Wanna integrate <phi_i|fxc|phi_j>, start looping over phi_j and do the first integration, then
2139 : !collect vxc and vxg and loop over phi_i for the second integration
2140 : !Note: we do not use the CG coefficients because they are only useful when there is a product
2141 : ! of Gaussians, which is not really the case here
2142 : !Note: the spherical orbitals for phi_i are distributed among the prcos of the current batch
2143 :
2144 112 : ALLOCATE (so(na, nr))
2145 140 : ALLOCATE (dso(na, nr, 3))
2146 84 : ALLOCATE (rexp(nr))
2147 :
2148 118 : ALLOCATE (vxc(ub))
2149 90 : ALLOCATE (vxg(ub))
2150 90 : ALLOCATE (int_so(ub))
2151 90 : DO i = 1, ub
2152 186 : ALLOCATE (vxc(i)%array(na, nr))
2153 186 : ALLOCATE (vxg(i)%array(na, nr, 3))
2154 248 : ALLOCATE (int_so(i)%array(nsotot, nsotot))
2155 11424552 : vxc(i)%array = 0.0_dp; vxg(i)%array = 0.0_dp; int_so(i)%array = 0.0_dp
2156 : END DO
2157 :
2158 84 : DO jset = 1, nset
2159 496 : DO jpgf = 1, npgf(jset)
2160 412 : startj = (jset - 1)*maxso + (jpgf - 1)*nsoset(lmax(jset))
2161 3600 : DO jso = nsoset(lmin(jset) - 1) + 1, nsoset(lmax(jset))
2162 3132 : l = indso(1, jso)
2163 :
2164 : !put the so phi_j and its gradient on the grid
2165 : !more efficient to recompute it rather than mp_bcast each chunk
2166 :
2167 494634 : rexp(1:nr) = EXP(-zet(jpgf, jset)*grid_atom%rad2(1:nr))
2168 : !$OMP PARALLEL DO COLLAPSE(2) DEFAULT(NONE), &
2169 : !$OMP SHARED(nr,na,so,dso,grid_atom,l,rexp,harmonics,jso,zet,jset,jpgf), &
2170 3132 : !$OMP PRIVATE(ir,ia,dir)
2171 : DO ir = 1, nr
2172 : DO ia = 1, na
2173 :
2174 : !so
2175 : so(ia, ir) = grid_atom%rad(ir)**l*rexp(ir)*harmonics%slm(ia, jso)
2176 :
2177 : !dso
2178 : dso(ia, ir, :) = 0.0_dp
2179 : DO dir = 1, 3
2180 : dso(ia, ir, dir) = dso(ia, ir, dir) &
2181 : + grid_atom%rad(ir)**(l - 1)*rexp(ir)*harmonics%dslm_dxyz(dir, ia, jso) &
2182 : - 2.0_dp*zet(jpgf, jset)*grid_atom%rad(ir)**(l + 1)*rexp(ir) &
2183 : *harmonics%a(dir, ia)*harmonics%slm(ia, jso)
2184 : END DO
2185 : END DO
2186 : END DO
2187 : !$OMP END PARALLEL DO
2188 :
2189 : !Perform the first integration (analytically)
2190 3132 : CALL get_vxc_vxg(vxc, vxg, so, dso, na, nr, rho_set, deriv_set, weight)
2191 :
2192 : !For a given phi_j, compute the second integration with all phi_i at once
2193 : !=> allows for efficient gemm to take place, especially since so are distributed
2194 3132 : nsoi = batch_info%nso_proc(ikind)
2195 9396 : bo = batch_info%so_bo(:, ikind)
2196 21924 : ALLOCATE (res(nsoi, nsoi), work(na, nsoi))
2197 120696000 : res = 0.0_dp; work = 0.0_dp
2198 :
2199 10152 : DO i = 1, ub
2200 :
2201 : !integrate so*Vxc and store in the int_so
2202 : CALL dgemm('N', 'N', na, nsoi, nr, 1.0_dp, vxc(i)%array(:, :), na, &
2203 7020 : gr(:, 1:nsoi), nr, 0.0_dp, work, na)
2204 : CALL dgemm('T', 'N', nsoi, nsoi, na, 1.0_dp, work, na, &
2205 7020 : ga(:, 1:nsoi), na, 0.0_dp, res, nsoi)
2206 802416 : int_so(i)%array(bo(1):bo(2), startj + jso) = get_diag(res)
2207 :
2208 31212 : DO dir = 1, 3
2209 :
2210 : ! integrate and sum up Vxg*dso
2211 : CALL dgemm('N', 'N', na, nsoi, nr, 1.0_dp, vxg(i)%array(:, :, dir), na, &
2212 21060 : dgr1(:, 1:nsoi), nr, 0.0_dp, work, na)
2213 : CALL dgemm('T', 'N', nsoi, nsoi, na, 1.0_dp, work, na, &
2214 21060 : dga1(:, 1:nsoi, dir), na, 0.0_dp, res, nsoi)
2215 21060 : CALL daxpy(nsoi, 1.0_dp, get_diag(res), 1, int_so(i)%array(bo(1):bo(2), startj + jso), 1)
2216 :
2217 : CALL dgemm('N', 'N', na, nsoi, nr, 1.0_dp, vxg(i)%array(:, :, dir), na, &
2218 21060 : dgr2(:, 1:nsoi), nr, 0.0_dp, work, na)
2219 : CALL dgemm('T', 'N', nsoi, nsoi, na, 1.0_dp, work, na, &
2220 21060 : dga2(:, 1:nsoi, dir), na, 0.0_dp, res, nsoi)
2221 28080 : CALL daxpy(nsoi, 1.0_dp, get_diag(res), 1, int_so(i)%array(bo(1):bo(2), startj + jso), 1)
2222 :
2223 : END DO
2224 :
2225 : END DO !i
2226 3544 : DEALLOCATE (res, work)
2227 :
2228 : END DO !jso
2229 : END DO !jpgf
2230 : END DO !jset
2231 :
2232 : !Contract into sgf and add to already computed LDA part of int_fxc
2233 28 : ri_sphi_so => xas_atom_env%ri_sphi_so(ikind)%array
2234 112 : ALLOCATE (int_sgf(nsgf, nsgf))
2235 90 : DO i = 1, ub
2236 4863350 : CALL batch_info%para_env%sum(int_so(i)%array)
2237 62 : CALL contract_so2sgf(int_sgf, int_so(i)%array, ri_basis, ri_sphi_so)
2238 90 : CALL daxpy(nsgf*nsgf, 1.0_dp, int_sgf, 1, int_fxc(iatom, i)%array, 1)
2239 : END DO
2240 :
2241 : !Clean-up
2242 90 : DO i = 1, ub
2243 62 : DEALLOCATE (vxc(i)%array)
2244 62 : DEALLOCATE (vxg(i)%array)
2245 90 : DEALLOCATE (int_so(i)%array)
2246 : END DO
2247 28 : DEALLOCATE (vxc, vxg, int_so)
2248 :
2249 28 : CALL timestop(handle)
2250 :
2251 84 : END SUBROUTINE integrate_gga_fxc
2252 :
2253 : ! **************************************************************************************************
2254 : !> \brief Computes the first integration of the GGA part of <phi_i|fxc|phi_j>, i.e. int fxc*phi_j.
2255 : !> The result is of the form Vxc - div(Vxg). Up to 3 results are returned, correspoinding to
2256 : !> f_aa, f_ab and (if open-shell) f_bb
2257 : !> \param vxc ...
2258 : !> \param vxg ...
2259 : !> \param so the spherical orbital on the grid
2260 : !> \param dso the derivative of the spherical orbital on the grid
2261 : !> \param na ...
2262 : !> \param nr ...
2263 : !> \param rho_set ...
2264 : !> \param deriv_set ...
2265 : !> \param weight the grid weight
2266 : !> \note This method is extremely similar to xc_calc_2nd_deriv of xc.F, but because it is a special
2267 : !> case that can be further optimized and because the interface of the original routine does
2268 : !> not fit this code, it has been re-written (no pw, no rho1_set but just the so, etc...)
2269 : ! **************************************************************************************************
2270 3132 : SUBROUTINE get_vxc_vxg(vxc, vxg, so, dso, na, nr, rho_set, deriv_set, weight)
2271 :
2272 : TYPE(cp_2d_r_p_type), DIMENSION(:) :: vxc
2273 : TYPE(cp_3d_r_p_type), DIMENSION(:) :: vxg
2274 : REAL(dp), DIMENSION(:, :) :: so
2275 : REAL(dp), DIMENSION(:, :, :) :: dso
2276 : INTEGER, INTENT(IN) :: na, nr
2277 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
2278 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
2279 : REAL(dp), DIMENSION(:, :), POINTER :: weight
2280 :
2281 : CHARACTER(len=*), PARAMETER :: routineN = 'get_vxc_vxg'
2282 :
2283 : INTEGER :: dir, handle, i, ia, ir, ub
2284 3132 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: dot_proda, dot_prodb, tmp
2285 3132 : REAL(dp), DIMENSION(:, :, :), POINTER :: d1e, d2e, norm_drhoa, norm_drhob
2286 : TYPE(xc_derivative_type), POINTER :: deriv
2287 :
2288 3132 : NULLIFY (norm_drhoa, norm_drhob, d2e, d1e, deriv)
2289 :
2290 3132 : CALL timeset(routineN, handle)
2291 :
2292 : !Note: this routines follows the order of the terms in equaiton (A.7) of Thomas Chassaing
2293 : ! thesis, except for the pure LDA terms that are dropped. The n^(1)_a and n^(1)_b
2294 : ! response densities are replaced by the spherical orbital.
2295 : ! The usual spin ordering is used: aa => 1, ab => 2 , bb => 3
2296 :
2297 : !point to the relevant components of rho_set
2298 3132 : ub = SIZE(vxc)
2299 3132 : norm_drhoa => rho_set%norm_drhoa
2300 3132 : norm_drhob => rho_set%norm_drhob
2301 :
2302 : !Some init
2303 10152 : DO i = 1, ub
2304 274657620 : vxc(i)%array = 0.0_dp
2305 823983012 : vxg(i)%array = 0.0_dp
2306 : END DO
2307 :
2308 25056 : ALLOCATE (tmp(na, nr), dot_proda(na, nr), dot_prodb(na, nr))
2309 218993340 : dot_proda = 0.0_dp; dot_prodb = 0.0_dp
2310 :
2311 : !Strategy: most terms are either multiplied by drhoa or drhob => group those first and then
2312 : ! multiply. Also most terms are multiplied by the dot product grad_n . grad_so, so
2313 : ! precompute it as well
2314 :
2315 : !$OMP PARALLEL DEFAULT(NONE), &
2316 : !$OMP SHARED(dot_proda,dot_prodb,tmp,vxc,vxg,deriv_set,rho_set,na,nr,norm_drhoa,norm_drhob,dso, &
2317 : !$OMP so,weight,ub), &
2318 3132 : !$OMP PRIVATE(ia,ir,dir,deriv,d1e,d2e)
2319 :
2320 : !Precompute the very common dot products grad_na . grad_so and grad_nb . grad_so
2321 : DO dir = 1, 3
2322 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2323 : DO ir = 1, nr
2324 : DO ia = 1, na
2325 : dot_proda(ia, ir) = dot_proda(ia, ir) + rho_set%drhoa(dir)%array(ia, ir, 1)*dso(ia, ir, dir)
2326 : dot_prodb(ia, ir) = dot_prodb(ia, ir) + rho_set%drhob(dir)%array(ia, ir, 1)*dso(ia, ir, dir)
2327 : END DO !ia
2328 : END DO !ir
2329 : !$OMP END DO NOWAIT
2330 : END DO !dir
2331 :
2332 : !Deal with f_aa
2333 :
2334 : !Vxc, first term
2335 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_norm_drhoa])
2336 : IF (ASSOCIATED(deriv)) THEN
2337 : CALL xc_derivative_get(deriv, deriv_data=d2e)
2338 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2339 : DO ir = 1, nr
2340 : DO ia = 1, na
2341 : vxc(1)%array(ia, ir) = d2e(ia, ir, 1)*dot_proda(ia, ir)
2342 : END DO !ia
2343 : END DO !ir
2344 : !$OMP END DO NOWAIT
2345 : END IF
2346 :
2347 : !Vxc, second term
2348 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_norm_drho])
2349 :
2350 : IF (ASSOCIATED(deriv)) THEN
2351 : CALL xc_derivative_get(deriv, deriv_data=d2e)
2352 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2353 : DO ir = 1, nr
2354 : DO ia = 1, na
2355 : vxc(1)%array(ia, ir) = vxc(1)%array(ia, ir) + d2e(ia, ir, 1)*dot_prodb(ia, ir)
2356 : END DO !ia
2357 : END DO !ir
2358 : !$OMP END DO NOWAIT
2359 : END IF
2360 :
2361 : !Vxc, take the grid weight into acocunt
2362 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2363 : DO ir = 1, nr
2364 : DO ia = 1, na
2365 : vxc(1)%array(ia, ir) = vxc(1)%array(ia, ir)*weight(ia, ir)
2366 : END DO !ia
2367 : END DO !ir
2368 : !$OMP END DO NOWAIT
2369 :
2370 : !Vxg, first term (to be multiplied by drhoa)
2371 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_norm_drhoa])
2372 : IF (ASSOCIATED(deriv)) THEN
2373 : CALL xc_derivative_get(deriv, deriv_data=d2e)
2374 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2375 : DO ir = 1, nr
2376 : DO ia = 1, na
2377 : tmp(ia, ir) = d2e(ia, ir, 1)*so(ia, ir)
2378 : END DO !ia
2379 : END DO !ir
2380 : !$OMP END DO NOWAIT
2381 : END IF
2382 :
2383 : !Vxg, second term (to be multiplied by drhoa)
2384 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_norm_drho])
2385 : IF (ASSOCIATED(deriv)) THEN
2386 : CALL xc_derivative_get(deriv, deriv_data=d2e)
2387 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2388 : DO ir = 1, nr
2389 : DO ia = 1, na
2390 : tmp(ia, ir) = tmp(ia, ir) + d2e(ia, ir, 1)*dot_prodb(ia, ir)
2391 : END DO !ia
2392 : END DO !ir
2393 : !$OMP END DO NOWAIT
2394 : END IF
2395 :
2396 : !Vxg, third term (to be multiplied by drhoa)
2397 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_norm_drhoa])
2398 : IF (ASSOCIATED(deriv)) THEN
2399 : CALL xc_derivative_get(deriv, deriv_data=d2e)
2400 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2401 : DO ir = 1, nr
2402 : DO ia = 1, na
2403 : tmp(ia, ir) = tmp(ia, ir) + d2e(ia, ir, 1)*dot_proda(ia, ir)
2404 : END DO !ia
2405 : END DO !ir
2406 : !$OMP END DO NOWAIT
2407 : END IF
2408 :
2409 : !Vxg, fourth term (to be multiplied by drhoa)
2410 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa])
2411 : IF (ASSOCIATED(deriv)) THEN
2412 : CALL xc_derivative_get(deriv, deriv_data=d1e)
2413 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2414 : DO ir = 1, nr
2415 : DO ia = 1, na
2416 : tmp(ia, ir) = tmp(ia, ir) - d1e(ia, ir, 1)*dot_proda(ia, ir) &
2417 : /MAX(norm_drhoa(ia, ir, 1), rho_set%drho_cutoff)**2
2418 : END DO !ia
2419 : END DO !ir
2420 : !$OMP END DO NOWAIT
2421 : END IF
2422 :
2423 : !put tmp*drhoa in Vxg (so that we can reuse it for drhob terms)
2424 : DO dir = 1, 3
2425 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2426 : DO ir = 1, nr
2427 : DO ia = 1, na
2428 : vxg(1)%array(ia, ir, dir) = tmp(ia, ir)*rho_set%drhoa(dir)%array(ia, ir, 1)
2429 : END DO !ia
2430 : END DO !ir
2431 : !$OMP END DO NOWAIT
2432 : END DO !dir
2433 :
2434 : !Vxg, fifth term (to be multiplied by drhob)
2435 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_norm_drho])
2436 : IF (ASSOCIATED(deriv)) THEN
2437 : CALL xc_derivative_get(deriv, deriv_data=d2e)
2438 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2439 : DO ir = 1, nr
2440 : DO ia = 1, na
2441 : tmp(ia, ir) = d2e(ia, ir, 1)*so(ia, ir)
2442 : END DO !ia
2443 : END DO !ir
2444 : !$OMP END DO NOWAIT
2445 : END IF
2446 :
2447 : !Vxg, sixth term (to be multiplied by drhob)
2448 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_norm_drho])
2449 : IF (ASSOCIATED(deriv)) THEN
2450 : CALL xc_derivative_get(deriv, deriv_data=d2e)
2451 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2452 : DO ir = 1, nr
2453 : DO ia = 1, na
2454 : tmp(ia, ir) = tmp(ia, ir) + d2e(ia, ir, 1)*dot_proda(ia, ir)
2455 : END DO !ia
2456 : END DO !ir
2457 : !$OMP END DO NOWAIT
2458 : END IF
2459 :
2460 : !Vxg, seventh term (to be multiplied by drhob)
2461 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drho])
2462 : IF (ASSOCIATED(deriv)) THEN
2463 : CALL xc_derivative_get(deriv, deriv_data=d2e)
2464 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2465 : DO ir = 1, nr
2466 : DO ia = 1, na
2467 : tmp(ia, ir) = tmp(ia, ir) + d2e(ia, ir, 1)*dot_prodb(ia, ir)
2468 : END DO !ia
2469 : END DO !ir
2470 : !$OMP END DO NOWAIT
2471 : END IF
2472 :
2473 : !put tmp*drhob in Vxg
2474 : DO dir = 1, 3
2475 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2476 : DO ir = 1, nr
2477 : DO ia = 1, na
2478 : vxg(1)%array(ia, ir, dir) = vxg(1)%array(ia, ir, dir) + tmp(ia, ir)*rho_set%drhob(dir)%array(ia, ir, 1)
2479 : END DO !ia
2480 : END DO !ir
2481 : !$OMP END DO NOWAIT
2482 : END DO !dir
2483 :
2484 : !Vxg, last term
2485 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa])
2486 : IF (ASSOCIATED(deriv)) THEN
2487 : CALL xc_derivative_get(deriv, deriv_data=d1e)
2488 : DO dir = 1, 3
2489 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2490 : DO ir = 1, nr
2491 : DO ia = 1, na
2492 : vxg(1)%array(ia, ir, dir) = vxg(1)%array(ia, ir, dir) + d1e(ia, ir, 1)*dso(ia, ir, dir)
2493 : END DO !ia
2494 : END DO !ir
2495 : !$OMP END DO NOWAIT
2496 : END DO !dir
2497 : END IF
2498 :
2499 : !Vxg, take the grid weight into account
2500 : DO dir = 1, 3
2501 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2502 : DO ir = 1, nr
2503 : DO ia = 1, na
2504 : vxg(1)%array(ia, ir, dir) = vxg(1)%array(ia, ir, dir)*weight(ia, ir)
2505 : END DO !ia
2506 : END DO !ir
2507 : !$OMP END DO NOWAIT
2508 : END DO !dir
2509 :
2510 : !Deal with fab
2511 :
2512 : !Vxc, first term
2513 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_norm_drhob])
2514 : IF (ASSOCIATED(deriv)) THEN
2515 : CALL xc_derivative_get(deriv, deriv_data=d2e)
2516 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2517 : DO ir = 1, nr
2518 : DO ia = 1, na
2519 : vxc(2)%array(ia, ir) = d2e(ia, ir, 1)*dot_prodb(ia, ir)
2520 : END DO !ia
2521 : END DO !ir
2522 : !$OMP END DO NOWAIT
2523 : END IF
2524 :
2525 : !Vxc, second term
2526 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_norm_drho])
2527 : IF (ASSOCIATED(deriv)) THEN
2528 : CALL xc_derivative_get(deriv, deriv_data=d2e)
2529 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2530 : DO ir = 1, nr
2531 : DO ia = 1, na
2532 : vxc(2)%array(ia, ir) = vxc(2)%array(ia, ir) + d2e(ia, ir, 1)*dot_proda(ia, ir)
2533 : END DO !ia
2534 : END DO !ir
2535 : !$OMP END DO NOWAIT
2536 : END IF
2537 :
2538 : !Vxc, take the grid weight into acocunt
2539 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2540 : DO ir = 1, nr
2541 : DO ia = 1, na
2542 : vxc(2)%array(ia, ir) = vxc(2)%array(ia, ir)*weight(ia, ir)
2543 : END DO !ia
2544 : END DO !ir
2545 : !$OMP END DO NOWAIT
2546 :
2547 : !Vxg, first term (to be multiplied by drhoa)
2548 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_norm_drhoa])
2549 : IF (ASSOCIATED(deriv)) THEN
2550 : CALL xc_derivative_get(deriv, deriv_data=d2e)
2551 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2552 : DO ir = 1, nr
2553 : DO ia = 1, na
2554 : tmp(ia, ir) = d2e(ia, ir, 1)*so(ia, ir)
2555 : END DO !ia
2556 : END DO !ir
2557 : !$OMP END DO NOWAIT
2558 : END IF
2559 :
2560 : !Vxg, second term (to be multiplied by drhoa)
2561 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_norm_drho])
2562 : IF (ASSOCIATED(deriv)) THEN
2563 : CALL xc_derivative_get(deriv, deriv_data=d2e)
2564 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2565 : DO ir = 1, nr
2566 : DO ia = 1, na
2567 : tmp(ia, ir) = tmp(ia, ir) + d2e(ia, ir, 1)*dot_proda(ia, ir)
2568 : END DO !ia
2569 : END DO !ir
2570 : !$OMP END DO NOWAIT
2571 : END IF
2572 :
2573 : !Vxg, third term (to be multiplied by drhoa)
2574 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_norm_drhob])
2575 : IF (ASSOCIATED(deriv)) THEN
2576 : CALL xc_derivative_get(deriv, deriv_data=d2e)
2577 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2578 : DO ir = 1, nr
2579 : DO ia = 1, na
2580 : tmp(ia, ir) = tmp(ia, ir) + d2e(ia, ir, 1)*dot_prodb(ia, ir)
2581 : END DO !ia
2582 : END DO !ir
2583 : !$OMP END DO NOWAIT
2584 : END IF
2585 :
2586 : !put tmp*drhoa in Vxg
2587 : DO dir = 1, 3
2588 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2589 : DO ir = 1, nr
2590 : DO ia = 1, na
2591 : vxg(2)%array(ia, ir, dir) = tmp(ia, ir)*rho_set%drhoa(dir)%array(ia, ir, 1)
2592 : END DO !ia
2593 : END DO !ir
2594 : !$OMP END DO NOWAIT
2595 : END DO !dir
2596 :
2597 : !Vxg, fourth term (to be multiplied by drhob)
2598 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_norm_drho])
2599 : IF (ASSOCIATED(deriv)) THEN
2600 : CALL xc_derivative_get(deriv, deriv_data=d2e)
2601 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2602 : DO ir = 1, nr
2603 : DO ia = 1, na
2604 : tmp(ia, ir) = d2e(ia, ir, 1)*so(ia, ir)
2605 : END DO !ia
2606 : END DO !ir
2607 : !$OMP END DO NOWAIT
2608 : END IF
2609 :
2610 : !Vxg, fifth term (to be multiplied by drhob)
2611 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drhob])
2612 : IF (ASSOCIATED(deriv)) THEN
2613 : CALL xc_derivative_get(deriv, deriv_data=d2e)
2614 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2615 : DO ir = 1, nr
2616 : DO ia = 1, na
2617 : tmp(ia, ir) = tmp(ia, ir) + d2e(ia, ir, 1)*dot_prodb(ia, ir)
2618 : END DO !ia
2619 : END DO !ir
2620 : !$OMP END DO NOWAIT
2621 : END IF
2622 :
2623 : !Vxg, sixth term (to be multiplied by drhob)
2624 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drho])
2625 : IF (ASSOCIATED(deriv)) THEN
2626 : CALL xc_derivative_get(deriv, deriv_data=d2e)
2627 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2628 : DO ir = 1, nr
2629 : DO ia = 1, na
2630 : tmp(ia, ir) = tmp(ia, ir) + d2e(ia, ir, 1)*dot_proda(ia, ir)
2631 : END DO !ia
2632 : END DO !ir
2633 : !$OMP END DO NOWAIT
2634 : END IF
2635 :
2636 : !put tmp*drhob in Vxg
2637 : DO dir = 1, 3
2638 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2639 : DO ir = 1, nr
2640 : DO ia = 1, na
2641 : vxg(2)%array(ia, ir, dir) = vxg(2)%array(ia, ir, dir) + tmp(ia, ir)*rho_set%drhob(dir)%array(ia, ir, 1)
2642 : END DO
2643 : END DO
2644 : !$OMP END DO NOWAIT
2645 : END DO
2646 :
2647 : !Vxg, last term
2648 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho])
2649 : IF (ASSOCIATED(deriv)) THEN
2650 : CALL xc_derivative_get(deriv, deriv_data=d1e)
2651 : DO dir = 1, 3
2652 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2653 : DO ir = 1, nr
2654 : DO ia = 1, na
2655 : vxg(2)%array(ia, ir, dir) = vxg(2)%array(ia, ir, dir) + d1e(ia, ir, 1)*dso(ia, ir, dir)
2656 : END DO !ia
2657 : END DO !ir
2658 : !$OMP END DO NOWAIT
2659 : END DO !dir
2660 : END IF
2661 :
2662 : !Vxg, take the grid weight into account
2663 : DO dir = 1, 3
2664 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2665 : DO ir = 1, nr
2666 : DO ia = 1, na
2667 : vxg(2)%array(ia, ir, dir) = vxg(2)%array(ia, ir, dir)*weight(ia, ir)
2668 : END DO !ia
2669 : END DO !ir
2670 : !$OMP END DO NOWAIT
2671 : END DO !dir
2672 :
2673 : !Deal with f_bb, if so required
2674 : IF (ub == 3) THEN
2675 :
2676 : !Vxc, first term
2677 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_norm_drhob])
2678 : IF (ASSOCIATED(deriv)) THEN
2679 : CALL xc_derivative_get(deriv, deriv_data=d2e)
2680 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2681 : DO ir = 1, nr
2682 : DO ia = 1, na
2683 : vxc(3)%array(ia, ir) = d2e(ia, ir, 1)*dot_prodb(ia, ir)
2684 : END DO !ia
2685 : END DO !ir
2686 : !$OMP END DO NOWAIT
2687 : END IF
2688 :
2689 : !Vxc, second term
2690 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_norm_drho])
2691 : IF (ASSOCIATED(deriv)) THEN
2692 : CALL xc_derivative_get(deriv, deriv_data=d2e)
2693 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2694 : DO ir = 1, nr
2695 : DO ia = 1, na
2696 : vxc(3)%array(ia, ir) = vxc(3)%array(ia, ir) + d2e(ia, ir, 1)*dot_proda(ia, ir)
2697 : END DO !i
2698 : END DO !ir
2699 : !$OMP END DO NOWAIT
2700 : END IF
2701 :
2702 : !Vxc, take the grid weight into acocunt
2703 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2704 : DO ir = 1, nr
2705 : DO ia = 1, na
2706 : vxc(3)%array(ia, ir) = vxc(3)%array(ia, ir)*weight(ia, ir)
2707 : END DO !ia
2708 : END DO !ir
2709 : !$OMP END DO NOWAIT
2710 :
2711 : !Vxg, first term (to be multiplied by drhob)
2712 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_norm_drhob])
2713 : IF (ASSOCIATED(deriv)) THEN
2714 : CALL xc_derivative_get(deriv, deriv_data=d2e)
2715 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2716 : DO ir = 1, nr
2717 : DO ia = 1, na
2718 : tmp(ia, ir) = d2e(ia, ir, 1)*so(ia, ir)
2719 : END DO !ia
2720 : END DO !ir
2721 : !$OMP END DO NOWAIT
2722 : END IF
2723 :
2724 : !Vxg, second term (to be multiplied by drhob)
2725 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_norm_drho])
2726 : IF (ASSOCIATED(deriv)) THEN
2727 : CALL xc_derivative_get(deriv, deriv_data=d2e)
2728 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2729 : DO ir = 1, nr
2730 : DO ia = 1, na
2731 : tmp(ia, ir) = tmp(ia, ir) + d2e(ia, ir, 1)*dot_proda(ia, ir)
2732 : END DO !ia
2733 : END DO !ir
2734 : !$OMP END DO NOWAIT
2735 : END IF
2736 :
2737 : !Vxg, third term (to be multiplied by drhob)
2738 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_norm_drhob])
2739 : IF (ASSOCIATED(deriv)) THEN
2740 : CALL xc_derivative_get(deriv, deriv_data=d2e)
2741 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2742 : DO ir = 1, nr
2743 : DO ia = 1, na
2744 : tmp(ia, ir) = tmp(ia, ir) + d2e(ia, ir, 1)*dot_prodb(ia, ir)
2745 : END DO !ia
2746 : END DO !ir
2747 : !$OMP END DO NOWAIT
2748 : END IF
2749 :
2750 : !Vxg, fourth term (to be multiplied by drhob)
2751 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob])
2752 : IF (ASSOCIATED(deriv)) THEN
2753 : CALL xc_derivative_get(deriv, deriv_data=d1e)
2754 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2755 : DO ir = 1, nr
2756 : DO ia = 1, na
2757 : tmp(ia, ir) = tmp(ia, ir) - d1e(ia, ir, 1)*dot_prodb(ia, ir) &
2758 : /MAX(norm_drhob(ia, ir, 1), rho_set%drho_cutoff)**2
2759 : END DO !ia
2760 : END DO !ir
2761 : !$OMP END DO NOWAIT
2762 : END IF
2763 :
2764 : !put tmp*drhob in Vxg (so that we can reuse it for drhoa terms)
2765 : DO dir = 1, 3
2766 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2767 : DO ir = 1, nr
2768 : DO ia = 1, na
2769 : vxg(3)%array(ia, ir, dir) = tmp(ia, ir)*rho_set%drhob(dir)%array(ia, ir, 1)
2770 : END DO !ia
2771 : END DO !ir
2772 : !$OMP END DO NOWAIT
2773 : END DO !dir
2774 :
2775 : !Vxg, fifth term (to be multiplied by drhoa)
2776 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_norm_drho])
2777 : IF (ASSOCIATED(deriv)) THEN
2778 : CALL xc_derivative_get(deriv, deriv_data=d2e)
2779 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2780 : DO ir = 1, nr
2781 : DO ia = 1, na
2782 : tmp(ia, ir) = d2e(ia, ir, 1)*so(ia, ir)
2783 : END DO !ia
2784 : END DO !ir
2785 : !$OMP END DO NOWAIT
2786 : END IF
2787 :
2788 : !Vxg, sixth term (to be multiplied by drhoa)
2789 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_norm_drho])
2790 : IF (ASSOCIATED(deriv)) THEN
2791 : CALL xc_derivative_get(deriv, deriv_data=d2e)
2792 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2793 : DO ir = 1, nr
2794 : DO ia = 1, na
2795 : tmp(ia, ir) = tmp(ia, ir) + d2e(ia, ir, 1)*dot_prodb(ia, ir)
2796 : END DO !ia
2797 : END DO !ir
2798 : !$OMP END DO NOWAIT
2799 : END IF
2800 :
2801 : !Vxg, seventh term (to be multiplied by drhoa)
2802 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_norm_drho])
2803 : IF (ASSOCIATED(deriv)) THEN
2804 : CALL xc_derivative_get(deriv, deriv_data=d2e)
2805 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2806 : DO ir = 1, nr
2807 : DO ia = 1, na
2808 : tmp(ia, ir) = tmp(ia, ir) + d2e(ia, ir, 1)*dot_proda(ia, ir)
2809 : END DO !ia
2810 : END DO !ir
2811 : !$OMP END DO NOWAIT
2812 : END IF
2813 :
2814 : !put tmp*drhoa in Vxg
2815 : DO dir = 1, 3
2816 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2817 : DO ir = 1, nr
2818 : DO ia = 1, na
2819 : vxg(3)%array(ia, ir, dir) = vxg(3)%array(ia, ir, dir) + &
2820 : tmp(ia, ir)*rho_set%drhoa(dir)%array(ia, ir, 1)
2821 : END DO !ia
2822 : END DO !ir
2823 : !$OMP END DO NOWAIT
2824 : END DO !dir
2825 :
2826 : !Vxg, last term
2827 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob])
2828 : IF (ASSOCIATED(deriv)) THEN
2829 : CALL xc_derivative_get(deriv, deriv_data=d1e)
2830 : DO dir = 1, 3
2831 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2832 : DO ir = 1, nr
2833 : DO ia = 1, na
2834 : vxg(3)%array(ia, ir, dir) = vxg(3)%array(ia, ir, dir) + d1e(ia, ir, 1)*dso(ia, ir, dir)
2835 : END DO !ia
2836 : END DO !ir
2837 : !$OMP END DO NOWAIT
2838 : END DO !dir
2839 : END IF
2840 :
2841 : !Vxg, take the grid weight into account
2842 : DO dir = 1, 3
2843 : !$OMP DO SCHEDULE(STATIC) COLLAPSE(2)
2844 : DO ir = 1, nr
2845 : DO ia = 1, na
2846 : vxg(3)%array(ia, ir, dir) = vxg(3)%array(ia, ir, dir)*weight(ia, ir)
2847 : END DO !ia
2848 : END DO !ir
2849 : !$OMP END DO NOWAIT
2850 : END DO !dir
2851 :
2852 : END IF !f_bb
2853 :
2854 : !$OMP END PARALLEL
2855 :
2856 3132 : CALL timestop(handle)
2857 :
2858 6264 : END SUBROUTINE get_vxc_vxg
2859 :
2860 : ! **************************************************************************************************
2861 : !> \brief Integrate the fxc kernel in the spin-conserving case, be it closed- or open-shell
2862 : !> \param int_fxc the array containing the (P|fxc|Q) integrals
2863 : !> \param iatom the index of the current excited atom
2864 : !> \param ikind the index of the current excited kind
2865 : !> \param deriv_set the set of functional derivatives
2866 : !> \param xas_atom_env ...
2867 : !> \param qs_env ...
2868 : ! **************************************************************************************************
2869 52 : SUBROUTINE integrate_sc_fxc(int_fxc, iatom, ikind, deriv_set, xas_atom_env, qs_env)
2870 :
2871 : TYPE(cp_2d_r_p_type), DIMENSION(:, :), POINTER :: int_fxc
2872 : INTEGER, INTENT(IN) :: iatom, ikind
2873 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
2874 : TYPE(xas_atom_env_type), POINTER :: xas_atom_env
2875 : TYPE(qs_environment_type), POINTER :: qs_env
2876 :
2877 : INTEGER :: i, maxso, na, nr, nset, nsotot, nspins, &
2878 : ri_nsgf
2879 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: fxc, int_so
2880 : REAL(dp), DIMENSION(:, :), POINTER :: ri_sphi_so
2881 52 : TYPE(cp_3d_r_p_type), ALLOCATABLE, DIMENSION(:) :: d2e
2882 : TYPE(dft_control_type), POINTER :: dft_control
2883 : TYPE(grid_atom_type), POINTER :: grid_atom
2884 : TYPE(gto_basis_set_type), POINTER :: ri_basis
2885 52 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2886 208 : TYPE(xc_derivative_p_type), DIMENSION(3) :: derivs
2887 :
2888 52 : NULLIFY (grid_atom, ri_basis, ri_sphi_so, qs_kind_set, dft_control)
2889 :
2890 : ! Initialization
2891 52 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, dft_control=dft_control)
2892 52 : grid_atom => xas_atom_env%grid_atom_set(ikind)%grid_atom
2893 52 : na = grid_atom%ng_sphere
2894 52 : nr = grid_atom%nr
2895 52 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis, basis_type="RI_XAS")
2896 52 : CALL get_gto_basis_set(ri_basis, nset=nset, maxso=maxso, nsgf=ri_nsgf)
2897 52 : nsotot = nset*maxso
2898 52 : ri_sphi_so => xas_atom_env%ri_sphi_so(ikind)%array
2899 52 : nspins = dft_control%nspins
2900 :
2901 : ! Get the second derivatives
2902 260 : ALLOCATE (d2e(3))
2903 52 : derivs(1)%deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa])
2904 52 : derivs(2)%deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob])
2905 52 : derivs(3)%deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob])
2906 208 : DO i = 1, 3
2907 208 : IF (ASSOCIATED(derivs(i)%deriv)) THEN
2908 156 : CALL xc_derivative_get(derivs(i)%deriv, deriv_data=d2e(i)%array)
2909 : END IF
2910 : END DO
2911 :
2912 : ! Allocate some work arrays
2913 208 : ALLOCATE (fxc(na, nr))
2914 208 : ALLOCATE (int_so(nsotot, nsotot))
2915 :
2916 : ! Integrate for all three derivatives, taking the grid weight into account
2917 : ! If closed shell, do not need to integrate beta-beta as it is the same as alpha-alpha
2918 162 : DO i = 1, nspins + 1
2919 7818122 : int_so = 0.0_dp
2920 110 : IF (ASSOCIATED(derivs(i)%deriv)) THEN
2921 3425408 : fxc(:, :) = d2e(i)%array(:, :, 1)*grid_atom%weight(:, :)
2922 110 : CALL integrate_so_prod(int_so, fxc, ikind, xas_atom_env, qs_env)
2923 : END IF
2924 :
2925 : !contract into sgf. Array allocated on current processor only
2926 440 : ALLOCATE (int_fxc(iatom, i)%array(ri_nsgf, ri_nsgf))
2927 815322 : int_fxc(iatom, i)%array = 0.0_dp
2928 162 : CALL contract_so2sgf(int_fxc(iatom, i)%array, int_so, ri_basis, ri_sphi_so)
2929 : END DO
2930 :
2931 156 : END SUBROUTINE integrate_sc_fxc
2932 :
2933 : ! **************************************************************************************************
2934 : !> \brief Integrate the fxc kernel in the spin-flip case (open-shell assumed). The spin-flip LDA
2935 : !> kernel reads: fxc = 1/(rhoa - rhob) * (dE/drhoa - dE/drhob)
2936 : !> \param int_fxc the array containing the (P|fxc|Q) integrals
2937 : !> \param iatom the index of the current excited atom
2938 : !> \param ikind the index of the current excited kind
2939 : !> \param rho_set the density in the atomic grid
2940 : !> \param deriv_set the set of functional derivatives
2941 : !> \param xas_atom_env ...
2942 : !> \param qs_env ...
2943 : ! **************************************************************************************************
2944 0 : SUBROUTINE integrate_sf_fxc(int_fxc, iatom, ikind, rho_set, deriv_set, xas_atom_env, qs_env)
2945 :
2946 : TYPE(cp_2d_r_p_type), DIMENSION(:, :), POINTER :: int_fxc
2947 : INTEGER, INTENT(IN) :: iatom, ikind
2948 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
2949 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
2950 : TYPE(xas_atom_env_type), POINTER :: xas_atom_env
2951 : TYPE(qs_environment_type), POINTER :: qs_env
2952 :
2953 : INTEGER :: ia, ir, maxso, na, nr, nset, nsotot, &
2954 : ri_nsgf
2955 0 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: fxc, int_so
2956 : REAL(dp), DIMENSION(:, :), POINTER :: ri_sphi_so
2957 0 : REAL(dp), DIMENSION(:, :, :), POINTER :: rhoa, rhob
2958 0 : TYPE(cp_3d_r_p_type), ALLOCATABLE, DIMENSION(:) :: d1e, d2e
2959 : TYPE(dft_control_type), POINTER :: dft_control
2960 : TYPE(grid_atom_type), POINTER :: grid_atom
2961 : TYPE(gto_basis_set_type), POINTER :: ri_basis
2962 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2963 : TYPE(xc_derivative_type), POINTER :: deriv
2964 :
2965 0 : NULLIFY (grid_atom, deriv, ri_basis, ri_sphi_so, qs_kind_set, rhoa, rhob, dft_control)
2966 :
2967 : ! Initialization
2968 0 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, dft_control=dft_control)
2969 0 : grid_atom => xas_atom_env%grid_atom_set(ikind)%grid_atom
2970 0 : na = grid_atom%ng_sphere
2971 0 : nr = grid_atom%nr
2972 0 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis, basis_type="RI_XAS")
2973 0 : CALL get_gto_basis_set(ri_basis, nset=nset, maxso=maxso, nsgf=ri_nsgf)
2974 0 : nsotot = nset*maxso
2975 0 : ri_sphi_so => xas_atom_env%ri_sphi_so(ikind)%array
2976 0 : rhoa => rho_set%rhoa
2977 0 : rhob => rho_set%rhob
2978 :
2979 0 : ALLOCATE (d1e(2))
2980 0 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa])
2981 0 : IF (ASSOCIATED(deriv)) THEN
2982 0 : CALL xc_derivative_get(deriv, deriv_data=d1e(1)%array)
2983 : END IF
2984 0 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob])
2985 0 : IF (ASSOCIATED(deriv)) THEN
2986 0 : CALL xc_derivative_get(deriv, deriv_data=d1e(2)%array)
2987 : END IF
2988 :
2989 : ! In case rhoa -> rhob, take the limit, which involves the (already computed) second derivatives
2990 0 : ALLOCATE (d2e(3))
2991 0 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa])
2992 0 : IF (ASSOCIATED(deriv)) THEN
2993 0 : CALL xc_derivative_get(deriv, deriv_data=d2e(1)%array)
2994 : END IF
2995 0 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhob])
2996 0 : IF (ASSOCIATED(deriv)) THEN
2997 0 : CALL xc_derivative_get(deriv, deriv_data=d2e(2)%array)
2998 : END IF
2999 0 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob])
3000 0 : IF (ASSOCIATED(deriv)) THEN
3001 0 : CALL xc_derivative_get(deriv, deriv_data=d2e(3)%array)
3002 : END IF
3003 :
3004 : !Compute the kernel on the grid. Already take weight into acocunt there
3005 0 : ALLOCATE (fxc(na, nr))
3006 : !$OMP PARALLEL DO COLLAPSE(2) SCHEDULE(STATIC) DEFAULT(NONE), &
3007 : !$OMP SHARED(grid_atom,fxc,d1e,d2e,dft_control,na,nr,rhoa,rhob), &
3008 0 : !$OMP PRIVATE(ia,ir)
3009 : DO ir = 1, nr
3010 : DO ia = 1, na
3011 :
3012 : !Need to be careful not to divide by zero. Assume that if rhoa == rhob, then
3013 : !take the limit fxc = 0.5* (f_aa + f_bb - 2f_ab)
3014 : IF (ABS(rhoa(ia, ir, 1) - rhob(ia, ir, 1)) > dft_control%qs_control%eps_rho_rspace) THEN
3015 : fxc(ia, ir) = grid_atom%weight(ia, ir)/(rhoa(ia, ir, 1) - rhob(ia, ir, 1)) &
3016 : *(d1e(1)%array(ia, ir, 1) - d1e(2)%array(ia, ir, 1))
3017 : ELSE
3018 : fxc(ia, ir) = 0.5_dp*grid_atom%weight(ia, ir)* &
3019 : (d2e(1)%array(ia, ir, 1) + d2e(3)%array(ia, ir, 1) - 2._dp*d2e(2)%array(ia, ir, 1))
3020 : END IF
3021 :
3022 : END DO
3023 : END DO
3024 : !$OMP END PARALLEL DO
3025 :
3026 : ! Integrate wrt to so
3027 0 : ALLOCATE (int_so(nsotot, nsotot))
3028 0 : int_so = 0.0_dp
3029 0 : CALL integrate_so_prod(int_so, fxc, ikind, xas_atom_env, qs_env)
3030 :
3031 : ! Contract into sgf. Array located on current processor only
3032 0 : ALLOCATE (int_fxc(iatom, 4)%array(ri_nsgf, ri_nsgf))
3033 0 : int_fxc(iatom, 4)%array = 0.0_dp
3034 0 : CALL contract_so2sgf(int_fxc(iatom, 4)%array, int_so, ri_basis, ri_sphi_so)
3035 :
3036 0 : END SUBROUTINE integrate_sf_fxc
3037 :
3038 : ! **************************************************************************************************
3039 : !> \brief Contract spherical orbitals to spherical Gaussians (so to sgf)
3040 : !> \param int_sgf the array with the sgf integrals
3041 : !> \param int_so the array with the so integrals (to contract)
3042 : !> \param basis the corresponding gto basis set
3043 : !> \param sphi_so the contraction coefficients for the s:
3044 : ! **************************************************************************************************
3045 304 : SUBROUTINE contract_so2sgf(int_sgf, int_so, basis, sphi_so)
3046 :
3047 : REAL(dp), DIMENSION(:, :) :: int_sgf, int_so
3048 : TYPE(gto_basis_set_type), POINTER :: basis
3049 : REAL(dp), DIMENSION(:, :) :: sphi_so
3050 :
3051 : INTEGER :: iset, jset, maxso, nset, nsoi, nsoj, &
3052 : sgfi, sgfj, starti, startj
3053 304 : INTEGER, DIMENSION(:), POINTER :: lmax, npgf, nsgf_set
3054 304 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf
3055 :
3056 304 : NULLIFY (nsgf_set, npgf, lmax, first_sgf)
3057 :
3058 : CALL get_gto_basis_set(basis, nset=nset, maxso=maxso, nsgf_set=nsgf_set, first_sgf=first_sgf, &
3059 304 : npgf=npgf, lmax=lmax)
3060 :
3061 1400 : DO iset = 1, nset
3062 1096 : starti = (iset - 1)*maxso + 1
3063 1096 : nsoi = npgf(iset)*nsoset(lmax(iset))
3064 1096 : sgfi = first_sgf(1, iset)
3065 :
3066 9484 : DO jset = 1, nset
3067 8084 : startj = (jset - 1)*maxso + 1
3068 8084 : nsoj = npgf(jset)*nsoset(lmax(jset))
3069 8084 : sgfj = first_sgf(1, jset)
3070 :
3071 : CALL ab_contract(int_sgf(sgfi:sgfi + nsgf_set(iset) - 1, sgfj:sgfj + nsgf_set(jset) - 1), &
3072 : int_so(starti:starti + nsoi - 1, startj:startj + nsoj - 1), &
3073 : sphi_so(:, sgfi:), sphi_so(:, sgfj:), nsoi, nsoj, &
3074 9180 : nsgf_set(iset), nsgf_set(jset))
3075 : END DO !jset
3076 : END DO !iset
3077 :
3078 304 : END SUBROUTINE contract_so2sgf
3079 :
3080 : ! **************************************************************************************************
3081 : !> \brief Integrate the product of spherical gaussian orbitals with the xc kernel on the atomic grid
3082 : !> \param intso the integral in terms of spherical orbitals
3083 : !> \param fxc the xc kernel at each grid point
3084 : !> \param ikind the kind of the atom we integrate for
3085 : !> \param xas_atom_env ...
3086 : !> \param qs_env ...
3087 : !> \note Largely copied from gaVxcgb_noGC. Rewritten here because we need our own atomic grid,
3088 : !> harmonics, basis set and we do not need the soft vxc. Could have tweaked the original, but
3089 : !> it would have been messy. Also we do not need rho_atom (too big and fancy for us)
3090 : !> We also go over the whole range of angular momentum l
3091 : ! **************************************************************************************************
3092 110 : SUBROUTINE integrate_so_prod(intso, fxc, ikind, xas_atom_env, qs_env)
3093 :
3094 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: intso
3095 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: fxc
3096 : INTEGER, INTENT(IN) :: ikind
3097 : TYPE(xas_atom_env_type), POINTER :: xas_atom_env
3098 : TYPE(qs_environment_type), POINTER :: qs_env
3099 :
3100 : CHARACTER(len=*), PARAMETER :: routineN = 'integrate_so_prod'
3101 :
3102 : INTEGER :: handle, ia, ic, icg, ipgf1, ipgf2, iset1, iset2, iso, iso1, iso2, l, ld, lmax12, &
3103 : lmin12, m1, m2, max_iso_not0, max_iso_not0_local, max_s_harm, maxl, maxso, n1, n2, na, &
3104 : ngau1, ngau2, nngau1, nr, nset, size1
3105 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list
3106 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list
3107 110 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
3108 110 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: g1, g2
3109 110 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: gfxcg, gg, matso
3110 110 : REAL(dp), DIMENSION(:, :), POINTER :: zet
3111 : REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG
3112 : TYPE(grid_atom_type), POINTER :: grid_atom
3113 : TYPE(gto_basis_set_type), POINTER :: ri_basis
3114 : TYPE(harmonics_atom_type), POINTER :: harmonics
3115 110 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3116 :
3117 110 : CALL timeset(routineN, handle)
3118 :
3119 110 : NULLIFY (grid_atom, harmonics, ri_basis, qs_kind_set, lmax, lmin, npgf, zet, my_CG)
3120 :
3121 : ! Initialization
3122 110 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
3123 110 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis, basis_type="RI_XAS")
3124 110 : grid_atom => xas_atom_env%grid_atom_set(ikind)%grid_atom
3125 110 : harmonics => xas_atom_env%harmonics_atom_set(ikind)%harmonics_atom
3126 :
3127 : CALL get_gto_basis_set(ri_basis, lmax=lmax, lmin=lmin, maxso=maxso, maxl=maxl, npgf=npgf, &
3128 110 : nset=nset, zet=zet)
3129 :
3130 110 : na = grid_atom%ng_sphere
3131 110 : nr = grid_atom%nr
3132 110 : my_CG => harmonics%my_CG
3133 110 : max_iso_not0 = harmonics%max_iso_not0
3134 110 : max_s_harm = harmonics%max_s_harm
3135 110 : CPASSERT(2*maxl <= indso(1, max_iso_not0))
3136 :
3137 770 : ALLOCATE (g1(nr), g2(nr), gg(nr, 0:2*maxl))
3138 440 : ALLOCATE (gfxcg(na, 0:2*maxl))
3139 440 : ALLOCATE (matso(nsoset(maxl), nsoset(maxl)))
3140 660 : ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm))
3141 :
3142 16300 : g1 = 0.0_dp
3143 16300 : g2 = 0.0_dp
3144 110 : m1 = 0
3145 : ! Loop over the product of so
3146 482 : DO iset1 = 1, nset
3147 372 : n1 = nsoset(lmax(iset1))
3148 372 : m2 = 0
3149 4308 : DO iset2 = 1, nset
3150 : CALL get_none0_cg_list(my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
3151 : max_s_harm, lmax(iset1) + lmax(iset2), cg_list, cg_n_list, &
3152 3936 : max_iso_not0_local)
3153 3936 : CPASSERT(max_iso_not0_local <= max_iso_not0)
3154 :
3155 3936 : n2 = nsoset(lmax(iset2))
3156 10616 : DO ipgf1 = 1, npgf(iset1)
3157 6680 : ngau1 = n1*(ipgf1 - 1) + m1
3158 6680 : size1 = nsoset(lmax(iset1)) - nsoset(lmin(iset1) - 1)
3159 6680 : nngau1 = nsoset(lmin(iset1) - 1) + ngau1
3160 :
3161 1122464 : g1(:) = EXP(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))
3162 39088 : DO ipgf2 = 1, npgf(iset2)
3163 28472 : ngau2 = n2*(ipgf2 - 1) + m2
3164 :
3165 4033592 : g2(:) = EXP(-zet(ipgf2, iset2)*grid_atom%rad2(1:nr))
3166 28472 : lmin12 = lmin(iset1) + lmin(iset2)
3167 28472 : lmax12 = lmax(iset1) + lmax(iset2)
3168 :
3169 : !get the gaussian product
3170 30574088 : gg = 0.0_dp
3171 28472 : IF (lmin12 == 0) THEN
3172 2254084 : gg(:, lmin12) = g1(:)*g2(:)
3173 : ELSE
3174 1779508 : gg(:, lmin12) = grid_atom%rad2l(1:nr, lmin12)*g1(:)*g2(:)
3175 : END IF
3176 :
3177 87864 : DO l = lmin12 + 1, lmax12
3178 8932664 : gg(:, l) = grid_atom%rad(1:nr)*gg(:, l - 1)
3179 : END DO
3180 :
3181 28472 : ld = lmax12 + 1
3182 : CALL dgemm('N', 'N', na, ld, nr, 1.0_dp, fxc(1:na, 1:nr), na, gg(:, 0:lmax12), &
3183 28472 : nr, 0.0_dp, gfxcg(1:na, 0:lmax12), na)
3184 :
3185 : !integrate
3186 11735352 : matso = 0.0_dp
3187 613708 : DO iso = 1, max_iso_not0_local
3188 3400726 : DO icg = 1, cg_n_list(iso)
3189 2787018 : iso1 = cg_list(1, icg, iso)
3190 2787018 : iso2 = cg_list(2, icg, iso)
3191 2787018 : l = indso(1, iso1) + indso(1, iso2)
3192 :
3193 604173530 : DO ia = 1, na
3194 : matso(iso1, iso2) = matso(iso1, iso2) + gfxcg(ia, l)* &
3195 603588294 : my_CG(iso1, iso2, iso)*harmonics%slm(ia, iso)
3196 : END DO !ia
3197 : END DO !icg
3198 : END DO !iso
3199 :
3200 : !write in integral matrix
3201 211304 : DO ic = nsoset(lmin(iset2) - 1) + 1, nsoset(lmax(iset2))
3202 176152 : iso1 = nsoset(lmin(iset1) - 1) + 1
3203 176152 : iso2 = ngau2 + ic
3204 204624 : CALL daxpy(size1, 1.0_dp, matso(iso1:, ic), 1, intso(nngau1 + 1:, iso2), 1)
3205 : END DO !ic
3206 :
3207 : END DO !ipgf2
3208 : END DO ! ipgf1
3209 4308 : m2 = m2 + maxso
3210 : END DO !iset2
3211 482 : m1 = m1 + maxso
3212 : END DO !iset1
3213 :
3214 110 : CALL timestop(handle)
3215 :
3216 330 : END SUBROUTINE integrate_so_prod
3217 :
3218 : ! **************************************************************************************************
3219 : !> \brief This routine computes the integral of a potential V wrt the derivative of the spherical
3220 : !> orbitals, that is <df/dx|V|dg/dy> on the atomic grid.
3221 : !> \param intso the integral in terms of the spherical orbitals (well, their derivative)
3222 : !> \param V the potential (put on the grid and weighted) to integrate
3223 : !> \param ikind the atomic kind for which we integrate
3224 : !> \param qs_env ...
3225 : !> \param soc_atom_env ...
3226 : ! **************************************************************************************************
3227 44 : SUBROUTINE integrate_so_dxdy_prod(intso, V, ikind, qs_env, soc_atom_env)
3228 :
3229 : REAL(dp), DIMENSION(:, :, :), INTENT(INOUT) :: intso
3230 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: V
3231 : INTEGER, INTENT(IN) :: ikind
3232 : TYPE(qs_environment_type), POINTER :: qs_env
3233 : TYPE(soc_atom_env_type), OPTIONAL, POINTER :: soc_atom_env
3234 :
3235 : INTEGER :: i, ipgf, iset, iso, j, jpgf, jset, jso, &
3236 : k, l, maxso, na, nr, nset, starti, &
3237 : startj
3238 44 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
3239 44 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: fga, fgr, r1, r2, work
3240 44 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: a1, a2
3241 44 : REAL(dp), DIMENSION(:, :), POINTER :: slm, zet
3242 44 : REAL(dp), DIMENSION(:, :, :), POINTER :: dslm_dxyz
3243 : TYPE(grid_atom_type), POINTER :: grid_atom
3244 : TYPE(gto_basis_set_type), POINTER :: basis
3245 : TYPE(harmonics_atom_type), POINTER :: harmonics
3246 44 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3247 :
3248 44 : NULLIFY (grid_atom, harmonics, basis, qs_kind_set, dslm_dxyz, slm, lmin, lmax, npgf, zet)
3249 :
3250 44 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
3251 44 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis, basis_type="ORB")
3252 44 : IF (PRESENT(soc_atom_env)) THEN
3253 8 : grid_atom => soc_atom_env%grid_atom_set(ikind)%grid_atom
3254 8 : harmonics => soc_atom_env%harmonics_atom_set(ikind)%harmonics_atom
3255 : ELSE
3256 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis, basis_type="ORB", harmonics=harmonics, &
3257 36 : grid_atom=grid_atom)
3258 : END IF
3259 :
3260 44 : na = grid_atom%ng_sphere
3261 44 : nr = grid_atom%nr
3262 :
3263 44 : slm => harmonics%slm
3264 44 : dslm_dxyz => harmonics%dslm_dxyz
3265 :
3266 : ! Getting what we need from the orbital basis
3267 : CALL get_gto_basis_set(gto_basis_set=basis, lmax=lmax, lmin=lmin, &
3268 44 : maxso=maxso, npgf=npgf, nset=nset, zet=zet)
3269 :
3270 : ! Separate the functions into purely r and purely angular parts, compute them all
3271 : ! and use matrix mutliplication for the integral. We use f for x derivative and g for y
3272 :
3273 : ! Separating the functions. Note that the radial part is the same for x and y derivatives
3274 308 : ALLOCATE (a1(na, nset*maxso, 3), a2(na, nset*maxso, 3))
3275 264 : ALLOCATE (r1(nr, nset*maxso), r2(nr, nset*maxso))
3276 929324 : a1 = 0.0_dp; a2 = 0.0_dp
3277 299588 : r1 = 0.0_dp; r2 = 0.0_dp
3278 :
3279 244 : DO iset = 1, nset
3280 722 : DO ipgf = 1, npgf(iset)
3281 478 : starti = (iset - 1)*maxso + (ipgf - 1)*nsoset(lmax(iset))
3282 1930 : DO iso = nsoset(lmin(iset) - 1) + 1, nsoset(lmax(iset))
3283 1252 : l = indso(1, iso)
3284 :
3285 : ! The x derivative of the spherical orbital, divided in angular and radial parts
3286 : ! Two of each are needed because d/dx(r^l Y_lm) * exp(-al*r^2) + r^l Y_lm * ! d/dx(exp-al*r^2)
3287 :
3288 : ! the purely radial part of d/dx(r^l Y_lm) * exp(-al*r^2) (same for y)
3289 61182 : r1(1:nr, starti + iso) = grid_atom%rad(1:nr)**(l - 1)*EXP(-zet(ipgf, iset)*grid_atom%rad2(1:nr))
3290 :
3291 : ! the purely radial part of r^l Y_lm * d/dx(exp-al*r^2) (same for y)
3292 : r2(1:nr, starti + iso) = -2.0_dp*zet(ipgf, iset)*grid_atom%rad(1:nr)**(l + 1) &
3293 61182 : *EXP(-zet(ipgf, iset)*grid_atom%rad2(1:nr))
3294 :
3295 5486 : DO i = 1, 3
3296 : ! the purely angular part of d/dx(r^l Y_lm) * exp(-al*r^2)
3297 191556 : a1(1:na, starti + iso, i) = dslm_dxyz(i, 1:na, iso)
3298 :
3299 : ! the purely angular part of r^l Y_lm * d/dx(exp-al*r^2)
3300 192808 : a2(1:na, starti + iso, i) = harmonics%a(i, 1:na)*slm(1:na, iso)
3301 : END DO
3302 :
3303 : END DO !iso
3304 : END DO !ipgf
3305 : END DO !iset
3306 :
3307 : ! Do the integration in terms of so using matrix products
3308 1308380 : intso = 0.0_dp
3309 132 : ALLOCATE (fga(na, 1))
3310 132 : ALLOCATE (fgr(nr, 1))
3311 88 : ALLOCATE (work(na, 1))
3312 6714 : fga = 0.0_dp; fgr = 0.0_dp; work = 0.0_dp
3313 :
3314 244 : DO iset = 1, nset
3315 1544 : DO jset = 1, nset
3316 4470 : DO ipgf = 1, npgf(iset)
3317 2970 : starti = (iset - 1)*maxso + (ipgf - 1)*nsoset(lmax(iset))
3318 11472 : DO jpgf = 1, npgf(jset)
3319 7202 : startj = (jset - 1)*maxso + (jpgf - 1)*nsoset(lmax(jset))
3320 :
3321 31778 : DO i = 1, 3
3322 21606 : j = MOD(i, 3) + 1
3323 21606 : k = MOD(i + 1, 3) + 1
3324 :
3325 84584 : DO iso = nsoset(lmin(iset) - 1) + 1, nsoset(lmax(iset))
3326 234174 : DO jso = nsoset(lmin(jset) - 1) + 1, nsoset(lmax(jset))
3327 :
3328 : !Two component per function => 4 terms in total
3329 :
3330 : ! take r1*a1(j) * V * r1*a1(k)
3331 7624818 : fgr(1:nr, 1) = r1(1:nr, starti + iso)*r1(1:nr, startj + jso)
3332 7996392 : fga(1:na, 1) = a1(1:na, starti + iso, j)*a1(1:na, startj + jso, k)
3333 :
3334 156792 : CALL dgemm('N', 'N', na, 1, nr, 1.0_dp, V, na, fgr, nr, 0.0_dp, work, na)
3335 : CALL dgemm('T', 'N', 1, 1, na, 1.0_dp, work, na, fga, na, 0.0_dp, &
3336 156792 : intso(starti + iso:, startj + jso, i), 1)
3337 :
3338 : ! add r1*a1(j) * V * r2*a2(k)
3339 7624818 : fgr(1:nr, 1) = r1(1:nr, starti + iso)*r2(1:nr, startj + jso)
3340 7996392 : fga(1:na, 1) = a1(1:na, starti + iso, j)*a2(1:na, startj + jso, k)
3341 :
3342 156792 : CALL dgemm('N', 'N', na, 1, nr, 1.0_dp, V, na, fgr, nr, 0.0_dp, work, na)
3343 : CALL dgemm('T', 'N', 1, 1, na, 1.0_dp, work, na, fga, na, 1.0_dp, &
3344 156792 : intso(starti + iso:, startj + jso, i), 1)
3345 :
3346 : ! add r2*a2(j) * V * r1*a1(k)
3347 7624818 : fgr(1:nr, 1) = r2(1:nr, starti + iso)*r1(1:nr, startj + jso)
3348 7996392 : fga(1:na, 1) = a2(1:na, starti + iso, j)*a1(1:na, startj + jso, k)
3349 :
3350 156792 : CALL dgemm('N', 'N', na, 1, nr, 1.0_dp, V, na, fgr, nr, 0.0_dp, work, na)
3351 : CALL dgemm('T', 'N', 1, 1, na, 1.0_dp, work, na, fga, na, 1.0_dp, &
3352 156792 : intso(starti + iso:, startj + jso, i), 1)
3353 :
3354 : ! add the last term: r2*a2(j) * V * r2*a2(k)
3355 7624818 : fgr(1:nr, 1) = r2(1:nr, starti + iso)*r2(1:nr, startj + jso)
3356 7996392 : fga(1:na, 1) = a2(1:na, starti + iso, j)*a2(1:na, startj + jso, k)
3357 :
3358 156792 : CALL dgemm('N', 'N', na, 1, nr, 1.0_dp, V, na, fgr, nr, 0.0_dp, work, na)
3359 : CALL dgemm('T', 'N', 1, 1, na, 1.0_dp, work, na, fga, na, 1.0_dp, &
3360 212568 : intso(starti + iso:, startj + jso, i), 1)
3361 :
3362 : END DO !jso
3363 : END DO !iso
3364 :
3365 : END DO !i
3366 : END DO !jpgf
3367 : END DO !ipgf
3368 : END DO !jset
3369 : END DO !iset
3370 :
3371 176 : DO i = 1, 3
3372 2616584 : intso(:, :, i) = intso(:, :, i) - TRANSPOSE(intso(:, :, i))
3373 : END DO
3374 :
3375 132 : END SUBROUTINE integrate_so_dxdy_prod
3376 :
3377 : ! **************************************************************************************************
3378 : !> \brief Computes the SOC matrix elements with respect to the ORB basis set for each atomic kind
3379 : !> and put them as the block diagonal of dbcsr_matrix
3380 : !> \param matrix_soc the matrix where the SOC is stored
3381 : !> \param xas_atom_env ...
3382 : !> \param qs_env ...
3383 : !> \param soc_atom_env ...
3384 : !> \note We compute: <da_dx|V\(4c^2-2V)|db_dy> - <da_dy|V\(4c^2-2V)|db_dx>, where V is a model
3385 : !> potential on the atomic grid. What we get is purely imaginary
3386 : ! **************************************************************************************************
3387 30 : SUBROUTINE integrate_soc_atoms(matrix_soc, xas_atom_env, qs_env, soc_atom_env)
3388 :
3389 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_soc
3390 : TYPE(xas_atom_env_type), OPTIONAL, POINTER :: xas_atom_env
3391 : TYPE(qs_environment_type), POINTER :: qs_env
3392 : TYPE(soc_atom_env_type), OPTIONAL, POINTER :: soc_atom_env
3393 :
3394 : CHARACTER(len=*), PARAMETER :: routineN = 'integrate_soc_atoms'
3395 :
3396 : INTEGER :: handle, i, iat, ikind, ir, jat, maxso, &
3397 : na, nkind, nr, nset, nsgf
3398 : LOGICAL :: all_potential_present
3399 : REAL(dp) :: zeff
3400 30 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: Vr
3401 30 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: V
3402 30 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: intso
3403 30 : REAL(dp), DIMENSION(:, :), POINTER :: sphi_so
3404 30 : REAL(dp), DIMENSION(:, :, :), POINTER :: intsgf
3405 30 : TYPE(cp_3d_r_p_type), DIMENSION(:), POINTER :: int_soc
3406 : TYPE(dbcsr_iterator_type) :: iter
3407 30 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
3408 : TYPE(grid_atom_type), POINTER :: grid
3409 : TYPE(gto_basis_set_type), POINTER :: basis
3410 30 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3411 30 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3412 :
3413 : !!DEBUG
3414 30 : CALL timeset(routineN, handle)
3415 :
3416 30 : NULLIFY (int_soc, basis, qs_kind_set, sphi_so, matrix_s, grid)
3417 30 : NULLIFY (particle_set)
3418 :
3419 : ! Initialization
3420 : CALL get_qs_env(qs_env, nkind=nkind, qs_kind_set=qs_kind_set, matrix_s=matrix_s, &
3421 30 : particle_set=particle_set)
3422 :
3423 : ! all_potential_present
3424 30 : CALL get_qs_kind_set(qs_kind_set, all_potential_present=all_potential_present)
3425 :
3426 : ! Loop over the kinds to compute the integrals
3427 134 : ALLOCATE (int_soc(nkind))
3428 74 : DO ikind = 1, nkind
3429 44 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis, basis_type="ORB", zeff=zeff)
3430 44 : IF (PRESENT(soc_atom_env)) THEN
3431 8 : grid => soc_atom_env%grid_atom_set(ikind)%grid_atom
3432 : ELSE
3433 36 : CALL get_qs_kind(qs_kind_set(ikind), grid_atom=grid)
3434 : END IF
3435 44 : CALL get_gto_basis_set(basis, nset=nset, maxso=maxso)
3436 220 : ALLOCATE (intso(nset*maxso, nset*maxso, 3))
3437 :
3438 : ! compute the model potential on the grid
3439 44 : nr = grid%nr
3440 44 : na = grid%ng_sphere
3441 :
3442 132 : ALLOCATE (Vr(nr))
3443 44 : CALL calculate_model_potential(Vr, grid, zeff)
3444 :
3445 : ! Compute V/(4c^2-2V) and weight it
3446 176 : ALLOCATE (V(na, nr))
3447 109082 : V = 0.0_dp
3448 2182 : DO ir = 1, nr
3449 : CALL daxpy(na, Vr(ir)/(4.0_dp*c_light_au**2 - 2.0_dp*Vr(ir)), grid%weight(:, ir), 1, &
3450 2182 : V(:, ir), 1)
3451 : END DO
3452 44 : DEALLOCATE (Vr)
3453 :
3454 : ! compute the integral <da_dx|...|db_dy> in terms of so
3455 44 : IF (PRESENT(xas_atom_env)) THEN
3456 36 : CALL integrate_so_dxdy_prod(intso, V, ikind, qs_env)
3457 : ELSE
3458 8 : CALL integrate_so_dxdy_prod(intso, V, ikind, qs_env, soc_atom_env)
3459 : END IF
3460 44 : DEALLOCATE (V)
3461 :
3462 : ! contract in terms of sgf
3463 44 : CALL get_gto_basis_set(basis, nsgf=nsgf)
3464 220 : ALLOCATE (int_soc(ikind)%array(nsgf, nsgf, 3))
3465 44 : intsgf => int_soc(ikind)%array
3466 44 : IF (PRESENT(xas_atom_env)) THEN
3467 36 : sphi_so => xas_atom_env%orb_sphi_so(ikind)%array
3468 : ELSE
3469 8 : sphi_so => soc_atom_env%orb_sphi_so(ikind)%array
3470 : END IF
3471 35888 : intsgf = 0.0_dp
3472 :
3473 176 : DO i = 1, 3
3474 176 : CALL contract_so2sgf(intsgf(:, :, i), intso(:, :, i), basis, sphi_so)
3475 : END DO
3476 :
3477 206 : DEALLOCATE (intso)
3478 : END DO !ikind
3479 :
3480 : ! Build the matrix_soc based on the matrix_s (but anti-symmetric)
3481 30 : IF ((PRESENT(xas_atom_env)) .OR. all_potential_present) THEN
3482 112 : DO i = 1, 3
3483 : CALL dbcsr_create(matrix_soc(i)%matrix, name="SOC MATRIX", template=matrix_s(1)%matrix, &
3484 112 : matrix_type=dbcsr_type_antisymmetric)
3485 : END DO
3486 : ! Iterate over its diagonal blocks and fill=it
3487 28 : CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
3488 158 : DO WHILE (dbcsr_iterator_blocks_left(iter))
3489 :
3490 130 : CALL dbcsr_iterator_next_block(iter, row=iat, column=jat)
3491 130 : IF (.NOT. iat == jat) CYCLE
3492 46 : ikind = particle_set(iat)%atomic_kind%kind_number
3493 :
3494 212 : DO i = 1, 3
3495 268 : CALL dbcsr_put_block(matrix_soc(i)%matrix, iat, iat, int_soc(ikind)%array(:, :, i))
3496 : END DO
3497 :
3498 : END DO !iat
3499 28 : CALL dbcsr_iterator_stop(iter)
3500 : ELSE ! pseudopotentials here
3501 8 : DO i = 1, 3
3502 : CALL dbcsr_create(matrix_soc(i)%matrix, name="SOC MATRIX", template=matrix_s(1)%matrix, &
3503 6 : matrix_type=dbcsr_type_no_symmetry)
3504 6 : CALL dbcsr_set(matrix_soc(i)%matrix, 0.0_dp)
3505 8 : CALL dbcsr_copy(matrix_soc(i)%matrix, soc_atom_env%soc_pp(i, 1)%matrix)
3506 : END DO
3507 : END IF
3508 :
3509 120 : DO i = 1, 3
3510 120 : CALL dbcsr_finalize(matrix_soc(i)%matrix)
3511 : END DO
3512 :
3513 : ! Clean-up
3514 74 : DO ikind = 1, nkind
3515 74 : DEALLOCATE (int_soc(ikind)%array)
3516 : END DO
3517 30 : DEALLOCATE (int_soc)
3518 :
3519 30 : CALL timestop(handle)
3520 :
3521 90 : END SUBROUTINE integrate_soc_atoms
3522 :
3523 : END MODULE xas_tdp_atom
|
