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