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