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