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