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