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