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