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