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