(git:d4bfb39)
Loading...
Searching...
No Matches
qs_vxc_atom.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief routines that build the integrals of the Vxc potential calculated
10!> for the atomic density in the basis set of spherical primitives
11! **************************************************************************************************
12MODULE qs_vxc_atom
13 USE atomic_kind_types, ONLY: atomic_kind_type,&
14 get_atomic_kind
15 USE basis_set_types, ONLY: get_gto_basis_set,&
16 gto_basis_set_type
17 USE cp_control_types, ONLY: dft_control_type
18 USE input_constants, ONLY: xc_none
19 USE input_section_types, ONLY: section_vals_get_subs_vals,&
20 section_vals_type,&
21 section_vals_val_get
22 USE kinds, ONLY: dp
23 USE memory_utilities, ONLY: reallocate
24 USE message_passing, ONLY: mp_para_env_type
25 USE orbital_pointers, ONLY: indso,&
26 nsoset
27 USE particle_types, ONLY: particle_type
28 USE paw_basis_types, ONLY: get_paw_basis_info
29 USE qs_environment_types, ONLY: get_qs_env,&
30 qs_environment_type
31 USE qs_force_types, ONLY: qs_force_type
32 USE qs_grid_atom, ONLY: grid_atom_type
33 USE qs_harmonics_atom, ONLY: get_none0_cg_list,&
34 harmonics_atom_type
35 USE qs_kind_types, ONLY: get_qs_kind,&
36 has_nlcc,&
37 qs_kind_type
38 USE qs_linres_types, ONLY: nablavks_atom_type
39 USE qs_rho_atom_types, ONLY: get_rho_atom,&
40 rho_atom_coeff,&
41 rho_atom_type
42 USE skala_gpw_functional, ONLY: native_skala_gapw_density_partition,&
43 skala_gapw_atom_vxc_of_r,&
44 skala_gapw_density_partition_hard_minus_soft,&
45 skala_gapw_density_partition_hard_only,&
46 skala_gapw_density_partition_none,&
47 skala_gapw_density_partition_soft_only,&
48 xc_section_uses_gauxc_model
49 USE util, ONLY: get_limit
50 USE virial_types, ONLY: virial_type
51 USE xc_atom, ONLY: fill_rho_set,&
52 vxc_of_r_epr,&
53 vxc_of_r_new,&
54 xc_2nd_deriv_of_r,&
55 xc_rho_set_atom_update
56 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
57 xc_dset_create,&
58 xc_dset_release,&
59 xc_dset_zero_all
60 USE xc_derivatives, ONLY: xc_functionals_get_needs
61 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
62 USE xc_rho_set_types, ONLY: xc_rho_set_create,&
63 xc_rho_set_release,&
64 xc_rho_set_type
65#include "./base/base_uses.f90"
66
67 IMPLICIT NONE
68
69 PRIVATE
70
71 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_vxc_atom'
72
73 TYPE tau_basis_cache_type
74 INTEGER :: maxso = 0, na = 0, nr = 0, nsatbas = 0, &
75 nset = 0
76 INTEGER, DIMENSION(:), POINTER :: lmax => null(), lmin => null(), &
77 n2oindex => null(), npgf => null(), &
78 o2nindex => null()
79 REAL(dp), DIMENSION(:, :), POINTER :: zet => null()
80 REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: grad
81 END TYPE tau_basis_cache_type
82
83 PUBLIC :: calculate_vxc_atom, &
84 calculate_vxc_atom_epr, &
85 calculate_xc_2nd_deriv_atom, &
86 calc_rho_angular, &
87 calculate_gfxc_atom, &
88 gfxc_atom_diff, &
89 gavxcgb_nogc
90
91CONTAINS
92
93! **************************************************************************************************
94!> \brief ...
95!> \param qs_env ...
96!> \param energy_only ...
97!> \param exc1 the on-body ex energy contribution
98!> \param adiabatic_rescale_factor ...
99!> \param kind_set_external provides a non-default kind_set to use
100!> \param rho_atom_set_external provides a non-default atomic density set to use
101!> \param xc_section_external provides an external non-default XC
102!> \param calculate_forces ...
103! **************************************************************************************************
104 SUBROUTINE calculate_vxc_atom(qs_env, energy_only, exc1, &
105 adiabatic_rescale_factor, kind_set_external, &
106 rho_atom_set_external, xc_section_external, calculate_forces)
107
108 TYPE(qs_environment_type), POINTER :: qs_env
109 LOGICAL, INTENT(IN) :: energy_only
110 REAL(dp), INTENT(INOUT) :: exc1
111 REAL(dp), INTENT(IN), OPTIONAL :: adiabatic_rescale_factor
112 TYPE(qs_kind_type), DIMENSION(:), OPTIONAL, &
113 POINTER :: kind_set_external
114 TYPE(rho_atom_type), DIMENSION(:), OPTIONAL, &
115 POINTER :: rho_atom_set_external
116 TYPE(section_vals_type), OPTIONAL, POINTER :: xc_section_external
117 LOGICAL, INTENT(IN), OPTIONAL :: calculate_forces
118
119 CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_vxc_atom'
120
121 INTEGER :: bo(2), gapw_density_partition, handle, &
122 iat, iatom, idir, ikind, ir, jdir, &
123 myfun, na, natom, nr, nspins, num_pe
124 INTEGER, DIMENSION(2, 3) :: bounds
125 INTEGER, DIMENSION(:), POINTER :: atom_list
126 LOGICAL :: accint, donlcc, evaluate_hard, evaluate_soft, gradient_f, lsd, &
127 my_calculate_forces, nlcc, paw_atom, skala_atom_grid, tau_f, use_virial
128 REAL(dp) :: agr, alpha, density_cut, exc_h, exc_s, &
129 gradient_cut, &
130 my_adiabatic_rescale_factor, tau_cut
131 REAL(dp), DIMENSION(1, 1, 1) :: tau_d
132 REAL(dp), DIMENSION(1, 1, 1, 1) :: rho_d
133 REAL(dp), DIMENSION(3) :: skala_atom_force_h, skala_atom_force_s
134 REAL(dp), DIMENSION(3, 3) :: skala_atom_virial, skala_atom_virial_h, &
135 skala_atom_virial_s
136 REAL(dp), DIMENSION(:, :), POINTER :: rho_nlcc, weight_h, weight_s
137 REAL(dp), DIMENSION(:, :, :), POINTER :: rho_h, rho_s, tau_h, tau_s, vtau_h, &
138 vtau_s, vxc_h, vxc_s
139 REAL(dp), DIMENSION(:, :, :, :), POINTER :: drho_h, drho_s, vxg_h, vxg_s
140 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
141 TYPE(dft_control_type), POINTER :: dft_control
142 TYPE(grid_atom_type), POINTER :: grid_atom
143 TYPE(gto_basis_set_type), POINTER :: basis_1c
144 TYPE(harmonics_atom_type), POINTER :: harmonics
145 TYPE(mp_para_env_type), POINTER :: para_env
146 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
147 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
148 TYPE(qs_kind_type), DIMENSION(:), POINTER :: my_kind_set
149 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s, int_hh, int_ss, r_h, r_s
150 TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r_h_d, r_s_d
151 TYPE(rho_atom_type), DIMENSION(:), POINTER :: my_rho_atom_set
152 TYPE(rho_atom_type), POINTER :: rho_atom
153 TYPE(section_vals_type), POINTER :: input, my_xc_section, xc_fun_section
154 TYPE(tau_basis_cache_type) :: tau_basis_cache
155 TYPE(virial_type), POINTER :: virial
156 TYPE(xc_derivative_set_type) :: deriv_set
157 TYPE(xc_rho_cflags_type) :: needs
158 TYPE(xc_rho_set_type) :: rho_set_h, rho_set_s
159
160! -------------------------------------------------------------------------
161
162 CALL timeset(routinen, handle)
163
164 NULLIFY (atom_list)
165 NULLIFY (my_kind_set)
166 NULLIFY (atomic_kind_set)
167 NULLIFY (grid_atom)
168 NULLIFY (force)
169 NULLIFY (harmonics)
170 NULLIFY (input)
171 NULLIFY (para_env)
172 NULLIFY (particle_set)
173 NULLIFY (rho_atom)
174 NULLIFY (my_rho_atom_set)
175 NULLIFY (rho_nlcc)
176 NULLIFY (virial)
177 my_calculate_forces = .false.
178 IF (PRESENT(calculate_forces)) my_calculate_forces = calculate_forces
179
180 IF (PRESENT(adiabatic_rescale_factor)) THEN
181 my_adiabatic_rescale_factor = adiabatic_rescale_factor
182 ELSE
183 my_adiabatic_rescale_factor = 1.0_dp
184 END IF
185
186 CALL get_qs_env(qs_env=qs_env, &
187 dft_control=dft_control, &
188 para_env=para_env, &
189 atomic_kind_set=atomic_kind_set, &
190 qs_kind_set=my_kind_set, &
191 input=input, &
192 particle_set=particle_set, &
193 virial=virial, &
194 rho_atom_set=my_rho_atom_set, &
195 force=force)
196
197 IF (PRESENT(kind_set_external)) my_kind_set => kind_set_external
198 IF (PRESENT(rho_atom_set_external)) my_rho_atom_set => rho_atom_set_external
199
200 nlcc = has_nlcc(my_kind_set)
201 accint = dft_control%qs_control%gapw_control%accurate_xcint
202
203 my_xc_section => section_vals_get_subs_vals(input, "DFT%XC")
204
205 IF (PRESENT(xc_section_external)) my_xc_section => xc_section_external
206
207 xc_fun_section => section_vals_get_subs_vals(my_xc_section, "XC_FUNCTIONAL")
208 CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", &
209 i_val=myfun)
210 skala_atom_grid = xc_section_uses_gauxc_model(my_xc_section)
211 gapw_density_partition = skala_gapw_density_partition_hard_minus_soft
212 IF (skala_atom_grid) THEN
213 gapw_density_partition = native_skala_gapw_density_partition(my_xc_section)
214 END IF
215 use_virial = ASSOCIATED(virial)
216 IF (use_virial) use_virial = my_calculate_forces .AND. &
217 virial%pv_calculate .AND. (.NOT. virial%pv_numer)
218
219 IF (myfun == xc_none) THEN
220 exc1 = 0.0_dp
221 my_rho_atom_set(:)%exc_h = 0.0_dp
222 my_rho_atom_set(:)%exc_s = 0.0_dp
223 ELSE
224 CALL section_vals_val_get(my_xc_section, "DENSITY_CUTOFF", &
225 r_val=density_cut)
226 CALL section_vals_val_get(my_xc_section, "GRADIENT_CUTOFF", &
227 r_val=gradient_cut)
228 CALL section_vals_val_get(my_xc_section, "TAU_CUTOFF", &
229 r_val=tau_cut)
230
231 lsd = dft_control%lsd
232 nspins = dft_control%nspins
233 needs = xc_functionals_get_needs(xc_fun_section, &
234 lsd=lsd, &
235 calc_potential=.true.)
236
237 gradient_f = (needs%drho .OR. needs%drho_spin) .OR. skala_atom_grid
238 tau_f = (needs%tau .OR. needs%tau_spin) .OR. skala_atom_grid
239
240 ! Initialize energy contribution from the one center XC terms to zero
241 exc1 = 0.0_dp
242
243 ! Nullify some pointers for work-arrays
244 NULLIFY (rho_h, drho_h, rho_s, drho_s, weight_h, weight_s)
245 NULLIFY (vxc_h, vxc_s, vxg_h, vxg_s)
246 NULLIFY (tau_h, tau_s)
247 NULLIFY (vtau_h, vtau_s)
248
249 ! Here starts the loop over all the atoms
250
251 DO ikind = 1, SIZE(atomic_kind_set)
252 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
253 CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom, &
254 harmonics=harmonics, grid_atom=grid_atom)
255 CALL get_qs_kind(my_kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
256
257 IF (.NOT. paw_atom) cycle
258
259 nr = grid_atom%nr
260 na = grid_atom%ng_sphere
261
262 ! Prepare the structures needed to calculate and store the xc derivatives
263
264 ! Array dimension: here anly one dimensional arrays are used,
265 ! i.e. only the first column of deriv_data is read.
266 ! The other to dimensions are set to size equal 1
267 bounds(1:2, 1:3) = 1
268 bounds(2, 1) = na
269 bounds(2, 2) = nr
270
271 ! set integration weights
272 IF (accint) THEN
273 weight_h => grid_atom%weight
274 alpha = dft_control%qs_control%gapw_control%aw(ikind)
275 IF (ASSOCIATED(grid_atom%gapw_weight_s)) THEN
276 IF (grid_atom%gapw_weight_alpha /= alpha) DEALLOCATE (grid_atom%gapw_weight_s)
277 END IF
278 IF (.NOT. ASSOCIATED(grid_atom%gapw_weight_s)) THEN
279 ALLOCATE (grid_atom%gapw_weight_s(na, nr))
280 DO ir = 1, nr
281 agr = 1.0_dp - exp(-alpha*grid_atom%rad2(ir))
282 grid_atom%gapw_weight_s(:, ir) = grid_atom%weight(:, ir)*agr
283 END DO
284 grid_atom%gapw_weight_alpha = alpha
285 END IF
286 weight_s => grid_atom%gapw_weight_s
287 ELSE
288 weight_h => grid_atom%weight
289 weight_s => grid_atom%weight
290 END IF
291
292 ! create a place where to put the derivatives
293 CALL xc_dset_create(deriv_set, local_bounds=bounds)
294 ! create the place where to store the argument for the functionals
295 CALL xc_rho_set_create(rho_set_h, bounds, rho_cutoff=density_cut, &
296 drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
297 CALL xc_rho_set_create(rho_set_s, bounds, rho_cutoff=density_cut, &
298 drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
299
300 ! allocate the required 3d arrays where to store rho and drho
301 CALL xc_rho_set_atom_update(rho_set_h, needs, nspins, bounds)
302 CALL xc_rho_set_atom_update(rho_set_s, needs, nspins, bounds)
303
304 CALL reallocate(rho_h, 1, na, 1, nr, 1, nspins)
305 CALL reallocate(rho_s, 1, na, 1, nr, 1, nspins)
306 CALL reallocate(vxc_h, 1, na, 1, nr, 1, nspins)
307 CALL reallocate(vxc_s, 1, na, 1, nr, 1, nspins)
308 !
309 IF (gradient_f) THEN
310 CALL reallocate(drho_h, 1, 4, 1, na, 1, nr, 1, nspins)
311 CALL reallocate(drho_s, 1, 4, 1, na, 1, nr, 1, nspins)
312 CALL reallocate(vxg_h, 1, 3, 1, na, 1, nr, 1, nspins)
313 CALL reallocate(vxg_s, 1, 3, 1, na, 1, nr, 1, nspins)
314 END IF
315
316 IF (tau_f) THEN
317 CALL create_tau_basis_cache(tau_basis_cache, grid_atom, basis_1c, harmonics)
318 CALL reallocate(tau_h, 1, na, 1, nr, 1, nspins)
319 CALL reallocate(tau_s, 1, na, 1, nr, 1, nspins)
320 CALL reallocate(vtau_h, 1, na, 1, nr, 1, nspins)
321 CALL reallocate(vtau_s, 1, na, 1, nr, 1, nspins)
322 END IF
323
324 ! NLCC: prepare rho and drho of the core charge for this KIND
325 donlcc = .false.
326 IF (nlcc) THEN
327 NULLIFY (rho_nlcc)
328 rho_nlcc => my_kind_set(ikind)%nlcc_pot
329 IF (ASSOCIATED(rho_nlcc)) donlcc = .true.
330 END IF
331
332 ! Distribute the atoms of this kind
333
334 num_pe = para_env%num_pe
335 bo = get_limit(natom, para_env%num_pe, para_env%mepos)
336
337 DO iat = bo(1), bo(2)
338 iatom = atom_list(iat)
339
340 my_rho_atom_set(iatom)%exc_h = 0.0_dp
341 my_rho_atom_set(iatom)%exc_s = 0.0_dp
342
343 rho_atom => my_rho_atom_set(iatom)
344 rho_h = 0.0_dp
345 rho_s = 0.0_dp
346 IF (gradient_f) THEN
347 NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
348 CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, &
349 rho_rad_s=r_s, drho_rad_h=dr_h, &
350 drho_rad_s=dr_s, rho_rad_h_d=r_h_d, &
351 rho_rad_s_d=r_s_d)
352 drho_h = 0.0_dp
353 drho_s = 0.0_dp
354 ELSE
355 NULLIFY (r_h, r_s)
356 CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s)
357 rho_d = 0.0_dp
358 END IF
359 IF (tau_f) THEN
360 !compute tau on the grid all at once
361 CALL calc_tau_atom(tau_h, tau_s, rho_atom, tau_basis_cache, nspins)
362 ELSE
363 tau_d = 0.0_dp
364 END IF
365
366 DO ir = 1, nr
367 CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_f, &
368 ir, r_h, r_s, rho_h, rho_s, dr_h, dr_s, &
369 r_h_d, r_s_d, drho_h, drho_s)
370 IF (donlcc) THEN
371 CALL calc_rho_nlcc(grid_atom, nspins, gradient_f, &
372 ir, rho_nlcc(:, 1), rho_h, rho_s, rho_nlcc(:, 2), drho_h, drho_s)
373 END IF
374 END DO
375
376 DO ir = 1, nr
377 IF (tau_f) THEN
378 CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_h, na, ir)
379 CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_s, na, ir)
380 ELSE IF (gradient_f) THEN
381 CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_d, na, ir)
382 CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_d, na, ir)
383 ELSE
384 CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, rho_d, tau_d, na, ir)
385 CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, rho_d, tau_d, na, ir)
386 END IF
387 END DO
388
389 evaluate_hard = .true.
390 evaluate_soft = .true.
391 skala_atom_force_h = 0.0_dp
392 skala_atom_force_s = 0.0_dp
393 skala_atom_virial_h = 0.0_dp
394 skala_atom_virial_s = 0.0_dp
395 IF (skala_atom_grid) THEN
396 SELECT CASE (gapw_density_partition)
397 CASE (skala_gapw_density_partition_hard_minus_soft)
398 CONTINUE
399 CASE (skala_gapw_density_partition_hard_only)
400 evaluate_soft = .false.
401 CASE (skala_gapw_density_partition_soft_only)
402 evaluate_hard = .false.
403 CASE (skala_gapw_density_partition_none)
404 evaluate_hard = .false.
405 evaluate_soft = .false.
406 CASE DEFAULT
407 CALL cp_abort(__location__, &
408 "Unknown GAUXC%NATIVE_GRID_GAPW_DENSITY_PARTITION value.")
409 END SELECT
410 END IF
411
412 !-------------------!
413 ! hard atom density !
414 !-------------------!
415 CALL xc_dset_zero_all(deriv_set)
416 IF (.NOT. evaluate_hard) THEN
417 exc_h = 0.0_dp
418 IF (.NOT. energy_only) THEN
419 vxc_h = 0.0_dp
420 IF (ASSOCIATED(vxg_h)) vxg_h = 0.0_dp
421 IF (ASSOCIATED(vtau_h)) vtau_h = 0.0_dp
422 END IF
423 ELSE IF (skala_atom_grid) THEN
424 CALL skala_gapw_atom_vxc_of_r( &
425 my_xc_section, grid_atom, para_env, particle_set(iatom)%r, &
426 rho_h, drho_h, tau_h, weight_h, lsd, nspins, na, nr, &
427 exc_h, vxc_h, vxg_h, vtau_h, energy_only=energy_only, &
428 atom_force=skala_atom_force_h, atom_virial=skala_atom_virial_h)
429 ELSE
430 CALL vxc_of_r_new(xc_fun_section, rho_set_h, deriv_set, 1, needs, weight_h, &
431 lsd, na, nr, exc_h, vxc_h, vxg_h, vtau_h, energy_only=energy_only, &
432 adiabatic_rescale_factor=my_adiabatic_rescale_factor)
433 END IF
434 rho_atom%exc_h = rho_atom%exc_h + exc_h
435
436 !-------------------!
437 ! soft atom density !
438 !-------------------!
439 CALL xc_dset_zero_all(deriv_set)
440 IF (.NOT. evaluate_soft) THEN
441 exc_s = 0.0_dp
442 IF (.NOT. energy_only) THEN
443 vxc_s = 0.0_dp
444 IF (ASSOCIATED(vxg_s)) vxg_s = 0.0_dp
445 IF (ASSOCIATED(vtau_s)) vtau_s = 0.0_dp
446 END IF
447 ELSE IF (skala_atom_grid) THEN
448 CALL skala_gapw_atom_vxc_of_r( &
449 my_xc_section, grid_atom, para_env, particle_set(iatom)%r, &
450 rho_s, drho_s, tau_s, weight_s, lsd, nspins, na, nr, &
451 exc_s, vxc_s, vxg_s, vtau_s, energy_only=energy_only, &
452 atom_force=skala_atom_force_s, atom_virial=skala_atom_virial_s)
453 ELSE
454 CALL vxc_of_r_new(xc_fun_section, rho_set_s, deriv_set, 1, needs, weight_s, &
455 lsd, na, nr, exc_s, vxc_s, vxg_s, vtau_s, energy_only=energy_only, &
456 adiabatic_rescale_factor=my_adiabatic_rescale_factor)
457 END IF
458 rho_atom%exc_s = rho_atom%exc_s + exc_s
459
460 ! Add contributions to the exc energy
461
462 exc1 = exc1 + rho_atom%exc_h - rho_atom%exc_s
463 IF (skala_atom_grid .AND. my_calculate_forces .AND. ASSOCIATED(force)) THEN
464 force(ikind)%rho_elec(:, iat) = force(ikind)%rho_elec(:, iat) + &
465 skala_atom_force_h - skala_atom_force_s
466 END IF
467 IF (skala_atom_grid .AND. use_virial) THEN
468 skala_atom_virial = skala_atom_virial_h - skala_atom_virial_s
469 DO idir = 1, 3
470 DO jdir = 1, 3
471 virial%pv_gapw(idir, jdir) = virial%pv_gapw(idir, jdir) + &
472 skala_atom_virial(idir, jdir)
473 virial%pv_virial(idir, jdir) = virial%pv_virial(idir, jdir) + &
474 skala_atom_virial(idir, jdir)
475 END DO
476 END DO
477 END IF
478
479 ! Integration to get the matrix elements relative to the vxc_atom
480 ! here the products with the primitives is done: gaVxcgb
481 ! internal transformation to get the integral in cartesian Gaussians
482
483 IF (.NOT. energy_only) THEN
484 NULLIFY (int_hh, int_ss)
485 CALL get_rho_atom(rho_atom=rho_atom, ga_vlocal_gb_h=int_hh, ga_vlocal_gb_s=int_ss)
486 IF (gradient_f) THEN
487 CALL gavxcgb_gc(vxc_h, vxc_s, vxg_h, vxg_s, int_hh, int_ss, &
488 grid_atom, basis_1c, harmonics, nspins)
489 ELSE
490 CALL gavxcgb_nogc(vxc_h, vxc_s, int_hh, int_ss, &
491 grid_atom, basis_1c, harmonics, nspins)
492 END IF
493 IF (tau_f) THEN
494 CALL dgavtaudgb(vtau_h, vtau_s, int_hh, int_ss, &
495 tau_basis_cache, nspins)
496 END IF
497 END IF ! energy_only
498 NULLIFY (r_h, r_s, dr_h, dr_s)
499 END DO ! iat
500
501 IF (tau_f) CALL release_tau_basis_cache(tau_basis_cache)
502
503 ! Release the xc structure used to store the xc derivatives
504 CALL xc_dset_release(deriv_set)
505 CALL xc_rho_set_release(rho_set_h)
506 CALL xc_rho_set_release(rho_set_s)
507 END DO ! ikind
508
509 CALL para_env%sum(exc1)
510
511 IF (ASSOCIATED(rho_h)) DEALLOCATE (rho_h)
512 IF (ASSOCIATED(rho_s)) DEALLOCATE (rho_s)
513 IF (ASSOCIATED(vxc_h)) DEALLOCATE (vxc_h)
514 IF (ASSOCIATED(vxc_s)) DEALLOCATE (vxc_s)
515
516 IF (gradient_f) THEN
517 IF (ASSOCIATED(drho_h)) DEALLOCATE (drho_h)
518 IF (ASSOCIATED(drho_s)) DEALLOCATE (drho_s)
519 IF (ASSOCIATED(vxg_h)) DEALLOCATE (vxg_h)
520 IF (ASSOCIATED(vxg_s)) DEALLOCATE (vxg_s)
521 END IF
522
523 IF (tau_f) THEN
524 IF (ASSOCIATED(tau_h)) DEALLOCATE (tau_h)
525 IF (ASSOCIATED(tau_s)) DEALLOCATE (tau_s)
526 IF (ASSOCIATED(vtau_h)) DEALLOCATE (vtau_h)
527 IF (ASSOCIATED(vtau_s)) DEALLOCATE (vtau_s)
528 END IF
529
530 END IF !xc_none
531
532 CALL timestop(handle)
533
534 END SUBROUTINE calculate_vxc_atom
535
536! **************************************************************************************************
537!> \brief ...
538!> \param qs_env ...
539!> \param exc1 the on-body ex energy contribution
540!> \param gradient_atom_set ...
541! **************************************************************************************************
542 SUBROUTINE calculate_vxc_atom_epr(qs_env, exc1, gradient_atom_set)
543
544 TYPE(qs_environment_type), POINTER :: qs_env
545 REAL(dp), INTENT(INOUT) :: exc1
546 TYPE(nablavks_atom_type), DIMENSION(:), POINTER :: gradient_atom_set
547
548 CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_vxc_atom_epr'
549
550 INTEGER :: bo(2), handle, ia, iat, iatom, idir, &
551 ikind, ir, ispin, myfun, na, natom, &
552 nr, nspins, num_pe
553 INTEGER, DIMENSION(2, 3) :: bounds
554 INTEGER, DIMENSION(:), POINTER :: atom_list
555 LOGICAL :: accint, donlcc, gradient_f, lsd, nlcc, &
556 paw_atom, tau_f
557 REAL(dp) :: agr, alpha, density_cut, exc_h, exc_s, &
558 gradient_cut, tau_cut
559 REAL(dp), DIMENSION(1, 1, 1) :: tau_d
560 REAL(dp), DIMENSION(1, 1, 1, 1) :: rho_d
561 REAL(dp), DIMENSION(:, :), POINTER :: rho_nlcc, weight_h, weight_s
562 REAL(dp), DIMENSION(:, :, :), POINTER :: rho_h, rho_s, tau_h, tau_s, vtau_h, &
563 vtau_s, vxc_h, vxc_s
564 REAL(dp), DIMENSION(:, :, :, :), POINTER :: drho_h, drho_s, vxg_h, vxg_s
565 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
566 TYPE(dft_control_type), POINTER :: dft_control
567 TYPE(grid_atom_type), POINTER :: grid_atom
568 TYPE(gto_basis_set_type), POINTER :: basis_1c
569 TYPE(harmonics_atom_type), POINTER :: harmonics
570 TYPE(mp_para_env_type), POINTER :: para_env
571 TYPE(qs_kind_type), DIMENSION(:), POINTER :: my_kind_set
572 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s, int_hh, int_ss, r_h, r_s
573 TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r_h_d, r_s_d
574 TYPE(rho_atom_type), DIMENSION(:), POINTER :: my_rho_atom_set
575 TYPE(rho_atom_type), POINTER :: rho_atom
576 TYPE(section_vals_type), POINTER :: input, my_xc_section, xc_fun_section
577 TYPE(tau_basis_cache_type) :: tau_basis_cache
578 TYPE(xc_derivative_set_type) :: deriv_set
579 TYPE(xc_rho_cflags_type) :: needs
580 TYPE(xc_rho_set_type) :: rho_set_h, rho_set_s
581
582! -------------------------------------------------------------------------
583
584 CALL timeset(routinen, handle)
585
586 NULLIFY (atom_list)
587 NULLIFY (my_kind_set)
588 NULLIFY (atomic_kind_set)
589 NULLIFY (grid_atom)
590 NULLIFY (harmonics)
591 NULLIFY (input)
592 NULLIFY (para_env)
593 NULLIFY (rho_atom)
594 NULLIFY (my_rho_atom_set)
595 NULLIFY (rho_nlcc)
596
597 CALL get_qs_env(qs_env=qs_env, &
598 dft_control=dft_control, &
599 para_env=para_env, &
600 atomic_kind_set=atomic_kind_set, &
601 qs_kind_set=my_kind_set, &
602 input=input, &
603 rho_atom_set=my_rho_atom_set)
604
605 nlcc = has_nlcc(my_kind_set)
606 accint = dft_control%qs_control%gapw_control%accurate_xcint
607
608 my_xc_section => section_vals_get_subs_vals(input, &
609 "PROPERTIES%LINRES%EPR%PRINT%G_TENSOR%XC")
610 xc_fun_section => section_vals_get_subs_vals(my_xc_section, "XC_FUNCTIONAL")
611 CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", &
612 i_val=myfun)
613
614 IF (myfun == xc_none) THEN
615 exc1 = 0.0_dp
616 my_rho_atom_set(:)%exc_h = 0.0_dp
617 my_rho_atom_set(:)%exc_s = 0.0_dp
618 ELSE
619 CALL section_vals_val_get(my_xc_section, "DENSITY_CUTOFF", &
620 r_val=density_cut)
621 CALL section_vals_val_get(my_xc_section, "GRADIENT_CUTOFF", &
622 r_val=gradient_cut)
623 CALL section_vals_val_get(my_xc_section, "TAU_CUTOFF", &
624 r_val=tau_cut)
625
626 lsd = dft_control%lsd
627 nspins = dft_control%nspins
628 needs = xc_functionals_get_needs(xc_fun_section, &
629 lsd=lsd, &
630 calc_potential=.true.)
631
632 ! whatever the xc, if epr_xc, drho_spin is needed
633 needs%drho_spin = .true.
634
635 gradient_f = (needs%drho .OR. needs%drho_spin)
636 tau_f = (needs%tau .OR. needs%tau_spin)
637
638 ! Initialize energy contribution from the one center XC terms to zero
639 exc1 = 0.0_dp
640
641 ! Nullify some pointers for work-arrays
642 NULLIFY (rho_h, drho_h, rho_s, drho_s, weight_h, weight_s)
643 NULLIFY (vxc_h, vxc_s, vxg_h, vxg_s)
644 NULLIFY (tau_h, tau_s)
645 NULLIFY (vtau_h, vtau_s)
646
647 ! Here starts the loop over all the atoms
648
649 DO ikind = 1, SIZE(atomic_kind_set)
650 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
651 CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom, &
652 harmonics=harmonics, grid_atom=grid_atom)
653 CALL get_qs_kind(my_kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
654
655 IF (.NOT. paw_atom) cycle
656
657 nr = grid_atom%nr
658 na = grid_atom%ng_sphere
659
660 ! Prepare the structures needed to calculate and store the xc derivatives
661
662 ! Array dimension: here anly one dimensional arrays are used,
663 ! i.e. only the first column of deriv_data is read.
664 ! The other to dimensions are set to size equal 1
665 bounds(1:2, 1:3) = 1
666 bounds(2, 1) = na
667 bounds(2, 2) = nr
668
669 ! set integration weights
670 IF (accint) THEN
671 weight_h => grid_atom%weight
672 alpha = dft_control%qs_control%gapw_control%aw(ikind)
673 IF (ASSOCIATED(grid_atom%gapw_weight_s)) THEN
674 IF (grid_atom%gapw_weight_alpha /= alpha) DEALLOCATE (grid_atom%gapw_weight_s)
675 END IF
676 IF (.NOT. ASSOCIATED(grid_atom%gapw_weight_s)) THEN
677 ALLOCATE (grid_atom%gapw_weight_s(na, nr))
678 DO ir = 1, nr
679 agr = 1.0_dp - exp(-alpha*grid_atom%rad2(ir))
680 grid_atom%gapw_weight_s(:, ir) = grid_atom%weight(:, ir)*agr
681 END DO
682 grid_atom%gapw_weight_alpha = alpha
683 END IF
684 weight_s => grid_atom%gapw_weight_s
685 ELSE
686 weight_h => grid_atom%weight
687 weight_s => grid_atom%weight
688 END IF
689
690 ! create a place where to put the derivatives
691 CALL xc_dset_create(deriv_set, local_bounds=bounds)
692 ! create the place where to store the argument for the functionals
693 CALL xc_rho_set_create(rho_set_h, bounds, rho_cutoff=density_cut, &
694 drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
695 CALL xc_rho_set_create(rho_set_s, bounds, rho_cutoff=density_cut, &
696 drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
697
698 ! allocate the required 3d arrays where to store rho and drho
699 CALL xc_rho_set_atom_update(rho_set_h, needs, nspins, bounds)
700 CALL xc_rho_set_atom_update(rho_set_s, needs, nspins, bounds)
701
702 CALL reallocate(rho_h, 1, na, 1, nr, 1, nspins)
703 CALL reallocate(rho_s, 1, na, 1, nr, 1, nspins)
704 CALL reallocate(vxc_h, 1, na, 1, nr, 1, nspins)
705 CALL reallocate(vxc_s, 1, na, 1, nr, 1, nspins)
706 !
707 IF (gradient_f) THEN
708 CALL reallocate(drho_h, 1, 4, 1, na, 1, nr, 1, nspins)
709 CALL reallocate(drho_s, 1, 4, 1, na, 1, nr, 1, nspins)
710 CALL reallocate(vxg_h, 1, 3, 1, na, 1, nr, 1, nspins)
711 CALL reallocate(vxg_s, 1, 3, 1, na, 1, nr, 1, nspins)
712 END IF
713
714 IF (tau_f) THEN
715 CALL create_tau_basis_cache(tau_basis_cache, grid_atom, basis_1c, harmonics)
716 CALL reallocate(tau_h, 1, na, 1, nr, 1, nspins)
717 CALL reallocate(tau_s, 1, na, 1, nr, 1, nspins)
718 CALL reallocate(vtau_h, 1, na, 1, nr, 1, nspins)
719 CALL reallocate(vtau_s, 1, na, 1, nr, 1, nspins)
720 END IF
721
722 ! NLCC: prepare rho and drho of the core charge for this KIND
723 donlcc = .false.
724 IF (nlcc) THEN
725 NULLIFY (rho_nlcc)
726 rho_nlcc => my_kind_set(ikind)%nlcc_pot
727 IF (ASSOCIATED(rho_nlcc)) donlcc = .true.
728 END IF
729
730 ! Distribute the atoms of this kind
731
732 num_pe = para_env%num_pe
733 bo = get_limit(natom, para_env%num_pe, para_env%mepos)
734
735 DO iat = bo(1), bo(2)
736 iatom = atom_list(iat)
737
738 my_rho_atom_set(iatom)%exc_h = 0.0_dp
739 my_rho_atom_set(iatom)%exc_s = 0.0_dp
740
741 rho_atom => my_rho_atom_set(iatom)
742 rho_h = 0.0_dp
743 rho_s = 0.0_dp
744 IF (gradient_f) THEN
745 NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
746 CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, &
747 rho_rad_s=r_s, drho_rad_h=dr_h, &
748 drho_rad_s=dr_s, rho_rad_h_d=r_h_d, &
749 rho_rad_s_d=r_s_d)
750 drho_h = 0.0_dp
751 drho_s = 0.0_dp
752 ELSE
753 NULLIFY (r_h, r_s)
754 CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s)
755 rho_d = 0.0_dp
756 END IF
757 IF (tau_f) THEN
758 !compute tau on the grid all at once
759 CALL calc_tau_atom(tau_h, tau_s, rho_atom, tau_basis_cache, nspins)
760 ELSE
761 tau_d = 0.0_dp
762 END IF
763
764 DO ir = 1, nr
765 CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_f, &
766 ir, r_h, r_s, rho_h, rho_s, dr_h, dr_s, &
767 r_h_d, r_s_d, drho_h, drho_s)
768 IF (donlcc) THEN
769 CALL calc_rho_nlcc(grid_atom, nspins, gradient_f, &
770 ir, rho_nlcc(:, 1), rho_h, rho_s, rho_nlcc(:, 2), drho_h, drho_s)
771 END IF
772 END DO
773 DO ir = 1, nr
774 IF (tau_f) THEN
775 CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_h, na, ir)
776 CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_s, na, ir)
777 ELSE IF (gradient_f) THEN
778 CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_d, na, ir)
779 CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_d, na, ir)
780 ELSE
781 CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, rho_d, tau_d, na, ir)
782 CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, rho_d, tau_d, na, ir)
783 END IF
784 END DO
785
786 !-------------------!
787 ! hard atom density !
788 !-------------------!
789 CALL xc_dset_zero_all(deriv_set)
790 CALL vxc_of_r_epr(xc_fun_section, rho_set_h, deriv_set, needs, weight_h, &
791 lsd, na, nr, exc_h, vxc_h, vxg_h, vtau_h)
792 rho_atom%exc_h = rho_atom%exc_h + exc_h
793
794 !-------------------!
795 ! soft atom density !
796 !-------------------!
797 CALL xc_dset_zero_all(deriv_set)
798 CALL vxc_of_r_epr(xc_fun_section, rho_set_s, deriv_set, needs, weight_s, &
799 lsd, na, nr, exc_s, vxc_s, vxg_s, vtau_s)
800 rho_atom%exc_s = rho_atom%exc_s + exc_s
801
802 DO ispin = 1, nspins
803 DO idir = 1, 3
804 DO ir = 1, nr
805 DO ia = 1, na
806 gradient_atom_set(iatom)%nablavks_vec_rad_h(idir, ispin)%r_coef(ir, ia) = &
807 gradient_atom_set(iatom)%nablavks_vec_rad_h(idir, ispin)%r_coef(ir, ia) &
808 + vxg_h(idir, ia, ir, ispin)
809 gradient_atom_set(iatom)%nablavks_vec_rad_s(idir, ispin)%r_coef(ir, ia) = &
810 gradient_atom_set(iatom)%nablavks_vec_rad_s(idir, ispin)%r_coef(ir, ia) &
811 + vxg_s(idir, ia, ir, ispin)
812 END DO ! ia
813 END DO ! ir
814 END DO ! idir
815 END DO ! ispin
816
817 ! Add contributions to the exc energy
818
819 exc1 = exc1 + rho_atom%exc_h - rho_atom%exc_s
820
821 ! Integration to get the matrix elements relative to the vxc_atom
822 ! here the products with the primitives is done: gaVxcgb
823 ! internal transformation to get the integral in cartesian Gaussians
824
825 NULLIFY (int_hh, int_ss)
826 CALL get_rho_atom(rho_atom=rho_atom, ga_vlocal_gb_h=int_hh, ga_vlocal_gb_s=int_ss)
827 IF (gradient_f) THEN
828 CALL gavxcgb_gc(vxc_h, vxc_s, vxg_h, vxg_s, int_hh, int_ss, &
829 grid_atom, basis_1c, harmonics, nspins)
830 ELSE
831 CALL gavxcgb_nogc(vxc_h, vxc_s, int_hh, int_ss, &
832 grid_atom, basis_1c, harmonics, nspins)
833 END IF
834 IF (tau_f) THEN
835 CALL dgavtaudgb(vtau_h, vtau_s, int_hh, int_ss, &
836 tau_basis_cache, nspins)
837 END IF
838 NULLIFY (r_h, r_s, dr_h, dr_s)
839 END DO ! iat
840
841 IF (tau_f) CALL release_tau_basis_cache(tau_basis_cache)
842
843 ! Release the xc structure used to store the xc derivatives
844 CALL xc_dset_release(deriv_set)
845 CALL xc_rho_set_release(rho_set_h)
846 CALL xc_rho_set_release(rho_set_s)
847 END DO ! ikind
848
849 CALL para_env%sum(exc1)
850
851 IF (ASSOCIATED(rho_h)) DEALLOCATE (rho_h)
852 IF (ASSOCIATED(rho_s)) DEALLOCATE (rho_s)
853 IF (ASSOCIATED(vxc_h)) DEALLOCATE (vxc_h)
854 IF (ASSOCIATED(vxc_s)) DEALLOCATE (vxc_s)
855
856 IF (gradient_f) THEN
857 IF (ASSOCIATED(drho_h)) DEALLOCATE (drho_h)
858 IF (ASSOCIATED(drho_s)) DEALLOCATE (drho_s)
859 IF (ASSOCIATED(vxg_h)) DEALLOCATE (vxg_h)
860 IF (ASSOCIATED(vxg_s)) DEALLOCATE (vxg_s)
861 END IF
862
863 IF (tau_f) THEN
864 IF (ASSOCIATED(tau_h)) DEALLOCATE (tau_h)
865 IF (ASSOCIATED(tau_s)) DEALLOCATE (tau_s)
866 IF (ASSOCIATED(vtau_h)) DEALLOCATE (vtau_h)
867 IF (ASSOCIATED(vtau_s)) DEALLOCATE (vtau_s)
868 END IF
869
870 END IF !xc_none
871
872 CALL timestop(handle)
873
874 END SUBROUTINE calculate_vxc_atom_epr
875
876! **************************************************************************************************
877!> \brief ...
878!> \param rho_atom_set ...
879!> \param rho1_atom_set ...
880!> \param qs_env ...
881!> \param xc_section ...
882!> \param para_env ...
883!> \param do_tddfpt2 New implementation of TDDFT.
884!> \param do_triplet ...
885!> \param do_sf ...
886!> \param kind_set_external ...
887! **************************************************************************************************
888 SUBROUTINE calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, xc_section, para_env, &
889 do_tddfpt2, do_triplet, do_sf, kind_set_external)
890
891 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set, rho1_atom_set
892 TYPE(qs_environment_type), POINTER :: qs_env
893 TYPE(section_vals_type), POINTER :: xc_section
894 TYPE(mp_para_env_type), INTENT(IN) :: para_env
895 LOGICAL, INTENT(IN), OPTIONAL :: do_tddfpt2, do_triplet, do_sf
896 TYPE(qs_kind_type), DIMENSION(:), OPTIONAL, &
897 POINTER :: kind_set_external
898
899 CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_xc_2nd_deriv_atom'
900
901 INTEGER :: atom, handle, iatom, ikind, ir, na, &
902 natom, nr, nspins
903 INTEGER, DIMENSION(2) :: local_loop_limit
904 INTEGER, DIMENSION(2, 3) :: bounds
905 INTEGER, DIMENSION(:), POINTER :: atom_list
906 LOGICAL :: accint, gradient_functional, lsd, &
907 my_do_sf, paw_atom, scale_rho, tau_f
908 REAL(kind=dp) :: agr, alpha, density_cut, gradient_cut, &
909 rtot, tau_cut
910 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
911 POINTER :: vtau_h, vtau_s, vxc_h, vxc_s
912 REAL(kind=dp), DIMENSION(1, 1, 1) :: rtau
913 REAL(kind=dp), DIMENSION(1, 1, 1, 1) :: rrho
914 REAL(kind=dp), DIMENSION(:, :), POINTER :: weight_h, weight_s
915 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: rho1_h, rho1_s, rho_h, rho_s, tau1_h, &
916 tau1_s, tau_h, tau_s
917 REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: drho1_h, drho1_s, drho_h, drho_s, vxg_h, &
918 vxg_s
919 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
920 TYPE(dft_control_type), POINTER :: dft_control
921 TYPE(grid_atom_type), POINTER :: grid_atom
922 TYPE(gto_basis_set_type), POINTER :: basis_1c
923 TYPE(harmonics_atom_type), POINTER :: harmonics
924 TYPE(qs_kind_type), DIMENSION(:), POINTER :: my_kind_set, qs_kind_set
925 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr1_h, dr1_s, dr_h, dr_s, int_hh, &
926 int_ss, r1_h, r1_s, r_h, r_s
927 TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r1_h_d, r1_s_d, r_h_d, r_s_d
928 TYPE(rho_atom_type), POINTER :: rho1_atom, rho_atom
929 TYPE(section_vals_type), POINTER :: input, xc_fun_section
930 TYPE(tau_basis_cache_type) :: tau_basis_cache
931 TYPE(xc_derivative_set_type) :: deriv_set
932 TYPE(xc_rho_cflags_type) :: needs
933 TYPE(xc_rho_set_type) :: rho1_set_h, rho1_set_s, rho_set_h, &
934 rho_set_s
935
936! -------------------------------------------------------------------------
937
938 CALL timeset(routinen, handle)
939
940 NULLIFY (qs_kind_set)
941 NULLIFY (rho_h, rho_s, drho_h, drho_s, weight_h, weight_s)
942 NULLIFY (rho1_h, rho1_s, drho1_h, drho1_s)
943 NULLIFY (vxc_h, vxc_s, vxg_h, vxg_s)
944 NULLIFY (tau_h, tau_s, tau1_h, tau1_s, vtau_h, vtau_s)
945
946 CALL get_qs_env(qs_env=qs_env, &
947 input=input, &
948 dft_control=dft_control, &
949 qs_kind_set=qs_kind_set, &
950 atomic_kind_set=atomic_kind_set)
951
952 IF (PRESENT(kind_set_external)) THEN
953 my_kind_set => kind_set_external
954 ELSE
955 my_kind_set => qs_kind_set
956 END IF
957
958 accint = dft_control%qs_control%gapw_control%accurate_xcint
959
960 CALL section_vals_val_get(input, "DFT%LSD", l_val=lsd)
961 CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", &
962 r_val=density_cut)
963 CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", &
964 r_val=gradient_cut)
965 CALL section_vals_val_get(xc_section, "TAU_CUTOFF", &
966 r_val=tau_cut)
967
968 my_do_sf = .false.
969 IF (PRESENT(do_sf)) my_do_sf = do_sf
970
971 xc_fun_section => section_vals_get_subs_vals(xc_section, &
972 "XC_FUNCTIONAL")
973 IF (lsd) THEN
974 nspins = 2
975 ELSE
976 nspins = 1
977 END IF
978
979 scale_rho = .false.
980 IF (PRESENT(do_tddfpt2) .AND. PRESENT(do_triplet)) THEN
981 IF (nspins == 1 .AND. do_triplet) THEN
982 lsd = .true.
983 scale_rho = .true.
984 END IF
985 ELSEIF (PRESENT(do_triplet)) THEN
986 IF (nspins == 1 .AND. do_triplet) lsd = .true.
987 END IF
988
989 needs = xc_functionals_get_needs(xc_fun_section, lsd=lsd, &
990 calc_potential=.true.)
991 gradient_functional = needs%drho .OR. needs%drho_spin
992 tau_f = (needs%tau .OR. needs%tau_spin)
993 IF (.NOT. tau_f) rtau = 0.0_dp
994
995 ! Here starts the loop over all the atoms
996 DO ikind = 1, SIZE(atomic_kind_set)
997
998 NULLIFY (atom_list, harmonics, grid_atom)
999 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
1000 CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom, &
1001 harmonics=harmonics, grid_atom=grid_atom)
1002 CALL get_qs_kind(my_kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
1003 IF (.NOT. paw_atom) cycle
1004
1005 nr = grid_atom%nr
1006 na = grid_atom%ng_sphere
1007
1008 ! set integration weights
1009 IF (accint) THEN
1010 weight_h => grid_atom%weight
1011 alpha = dft_control%qs_control%gapw_control%aw(ikind)
1012 IF (ASSOCIATED(grid_atom%gapw_weight_s)) THEN
1013 IF (grid_atom%gapw_weight_alpha /= alpha) DEALLOCATE (grid_atom%gapw_weight_s)
1014 END IF
1015 IF (.NOT. ASSOCIATED(grid_atom%gapw_weight_s)) THEN
1016 ALLOCATE (grid_atom%gapw_weight_s(na, nr))
1017 DO ir = 1, nr
1018 agr = 1.0_dp - exp(-alpha*grid_atom%rad2(ir))
1019 grid_atom%gapw_weight_s(:, ir) = grid_atom%weight(:, ir)*agr
1020 END DO
1021 grid_atom%gapw_weight_alpha = alpha
1022 END IF
1023 weight_s => grid_atom%gapw_weight_s
1024 ELSE
1025 weight_h => grid_atom%weight
1026 weight_s => grid_atom%weight
1027 END IF
1028
1029 ! Array dimension: here anly one dimensional arrays are used,
1030 ! i.e. only the first column of deriv_data is read.
1031 ! The other to dimensions are set to size equal 1.
1032 bounds(1:2, 1:3) = 1
1033 bounds(2, 1) = na
1034 bounds(2, 2) = nr
1035
1036 CALL xc_dset_create(deriv_set, local_bounds=bounds)
1037 CALL xc_rho_set_create(rho_set_h, bounds, rho_cutoff=density_cut, &
1038 drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
1039 CALL xc_rho_set_create(rho_set_s, bounds, rho_cutoff=density_cut, &
1040 drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
1041 CALL xc_rho_set_create(rho1_set_h, bounds, rho_cutoff=density_cut, &
1042 drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
1043 CALL xc_rho_set_create(rho1_set_s, bounds, rho_cutoff=density_cut, &
1044 drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
1045
1046 ! allocate the required 3d arrays where to store rho and drho
1047 IF (nspins == 1 .AND. .NOT. lsd) THEN
1048 CALL xc_rho_set_atom_update(rho_set_h, needs, 1, bounds)
1049 CALL xc_rho_set_atom_update(rho1_set_h, needs, 1, bounds)
1050 CALL xc_rho_set_atom_update(rho_set_s, needs, 1, bounds)
1051 CALL xc_rho_set_atom_update(rho1_set_s, needs, 1, bounds)
1052 ELSE
1053 CALL xc_rho_set_atom_update(rho_set_h, needs, 2, bounds)
1054 CALL xc_rho_set_atom_update(rho1_set_h, needs, 2, bounds)
1055 CALL xc_rho_set_atom_update(rho_set_s, needs, 2, bounds)
1056 CALL xc_rho_set_atom_update(rho1_set_s, needs, 2, bounds)
1057 END IF
1058
1059 ALLOCATE (rho_h(1:na, 1:nr, 1:nspins), rho1_h(1:na, 1:nr, 1:nspins), &
1060 rho_s(1:na, 1:nr, 1:nspins), rho1_s(1:na, 1:nr, 1:nspins))
1061
1062 ALLOCATE (vxc_h(1:na, 1:nr, 1:nspins), vxc_s(1:na, 1:nr, 1:nspins))
1063 vxc_h = 0.0_dp
1064 vxc_s = 0.0_dp
1065
1066 IF (tau_f) THEN
1067 CALL create_tau_basis_cache(tau_basis_cache, grid_atom, basis_1c, harmonics)
1068 ALLOCATE (tau_h(1:na, 1:nr, 1:nspins), tau1_h(1:na, 1:nr, 1:nspins), &
1069 tau_s(1:na, 1:nr, 1:nspins), tau1_s(1:na, 1:nr, 1:nspins))
1070 ALLOCATE (vtau_h(1:na, 1:nr, 1:nspins), vtau_s(1:na, 1:nr, 1:nspins))
1071 END IF
1072
1073 IF (gradient_functional) THEN
1074 ALLOCATE (drho_h(1:4, 1:na, 1:nr, 1:nspins), drho1_h(1:4, 1:na, 1:nr, 1:nspins), &
1075 drho_s(1:4, 1:na, 1:nr, 1:nspins), drho1_s(1:4, 1:na, 1:nr, 1:nspins))
1076 ALLOCATE (vxg_h(1:3, 1:na, 1:nr, 1:nspins), vxg_s(1:3, 1:na, 1:nr, 1:nspins))
1077 ELSE
1078 ALLOCATE (drho_h(1, 1, 1, 1), drho1_h(1, 1, 1, 1), &
1079 drho_s(1, 1, 1, 1), drho1_s(1, 1, 1, 1))
1080 ALLOCATE (vxg_h(1, 1, 1, 1), vxg_s(1, 1, 1, 1))
1081 rrho = 0.0_dp
1082 END IF
1083 vxg_h = 0.0_dp
1084 vxg_s = 0.0_dp
1085
1086 ! parallelization
1087 local_loop_limit = get_limit(natom, para_env%num_pe, para_env%mepos)
1088
1089 DO iatom = local_loop_limit(1), local_loop_limit(2) !1,natom
1090 atom = atom_list(iatom)
1091
1092 rho_atom_set(atom)%exc_h = 0.0_dp
1093 rho_atom_set(atom)%exc_s = 0.0_dp
1094 rho1_atom_set(atom)%exc_h = 0.0_dp
1095 rho1_atom_set(atom)%exc_s = 0.0_dp
1096
1097 rho_atom => rho_atom_set(atom)
1098 rho1_atom => rho1_atom_set(atom)
1099 NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
1100 NULLIFY (r1_h, r1_s, dr1_h, dr1_s, r1_h_d, r1_s_d)
1101 rho_h = 0.0_dp
1102 rho_s = 0.0_dp
1103 rho1_h = 0.0_dp
1104 rho1_s = 0.0_dp
1105 IF (gradient_functional) THEN
1106 CALL get_rho_atom(rho_atom=rho_atom, &
1107 rho_rad_h=r_h, rho_rad_s=r_s, &
1108 drho_rad_h=dr_h, drho_rad_s=dr_s, &
1109 rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
1110 CALL get_rho_atom(rho_atom=rho1_atom, &
1111 rho_rad_h=r1_h, rho_rad_s=r1_s, &
1112 drho_rad_h=dr1_h, drho_rad_s=dr1_s, &
1113 rho_rad_h_d=r1_h_d, rho_rad_s_d=r1_s_d)
1114 drho_h = 0.0_dp; drho_s = 0.0_dp
1115 drho1_h = 0.0_dp; drho1_s = 0.0_dp
1116 ELSE
1117 CALL get_rho_atom(rho_atom=rho_atom, &
1118 rho_rad_h=r_h, rho_rad_s=r_s)
1119 CALL get_rho_atom(rho_atom=rho1_atom, &
1120 rho_rad_h=r1_h, rho_rad_s=r1_s)
1121 END IF
1122
1123 rtot = 0.0_dp
1124
1125 DO ir = 1, nr
1126 CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_functional, &
1127 ir, r_h, r_s, rho_h, rho_s, dr_h, dr_s, r_h_d, r_s_d, &
1128 drho_h, drho_s)
1129 CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_functional, &
1130 ir, r1_h, r1_s, rho1_h, rho1_s, dr1_h, dr1_s, r1_h_d, r1_s_d, &
1131 drho1_h, drho1_s)
1132 END DO
1133 IF (tau_f) THEN
1134 CALL calc_tau_atom(tau_h, tau_s, rho_atom, tau_basis_cache, nspins)
1135 CALL calc_tau_atom(tau1_h, tau1_s, rho1_atom, tau_basis_cache, nspins)
1136 END IF
1137 IF (scale_rho) THEN
1138 rho_h = 2.0_dp*rho_h
1139 rho_s = 2.0_dp*rho_s
1140 IF (gradient_functional) THEN
1141 drho_h = 2.0_dp*drho_h
1142 drho_s = 2.0_dp*drho_s
1143 END IF
1144 IF (tau_f) THEN
1145 tau_h = 2.0_dp*tau_h
1146 tau_s = 2.0_dp*tau_s
1147 END IF
1148 END IF
1149
1150 DO ir = 1, nr
1151 IF (tau_f) THEN
1152 CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_h, na, ir)
1153 CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, drho1_h, tau1_h, na, ir)
1154 CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_s, na, ir)
1155 CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, drho1_s, tau1_s, na, ir)
1156 ELSE IF (gradient_functional) THEN
1157 CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, rtau, na, ir)
1158 CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, drho1_h, rtau, na, ir)
1159 CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, rtau, na, ir)
1160 CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, drho1_s, rtau, na, ir)
1161 ELSE
1162 CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, rrho, rtau, na, ir)
1163 CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, rrho, rtau, na, ir)
1164 CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, rrho, rtau, na, ir)
1165 CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, rrho, rtau, na, ir)
1166 END IF
1167 END DO
1168
1169 CALL xc_2nd_deriv_of_r(xc_section=xc_section, &
1170 rho_set=rho_set_h, rho1_set=rho1_set_h, &
1171 deriv_set=deriv_set, &
1172 w=weight_h, vxc=vxc_h, vxg=vxg_h, vtau=vtau_h, do_triplet=do_triplet, &
1173 do_sf=my_do_sf)
1174 CALL xc_2nd_deriv_of_r(xc_section=xc_section, &
1175 rho_set=rho_set_s, rho1_set=rho1_set_s, &
1176 deriv_set=deriv_set, &
1177 w=weight_s, vxc=vxc_s, vxg=vxg_s, vtau=vtau_s, do_triplet=do_triplet, &
1178 do_sf=my_do_sf)
1179
1180 CALL get_rho_atom(rho_atom=rho1_atom, ga_vlocal_gb_h=int_hh, ga_vlocal_gb_s=int_ss)
1181 IF (gradient_functional) THEN
1182 CALL gavxcgb_gc(vxc_h, vxc_s, vxg_h, vxg_s, int_hh, int_ss, &
1183 grid_atom, basis_1c, harmonics, nspins)
1184 ELSE
1185 CALL gavxcgb_nogc(vxc_h, vxc_s, int_hh, int_ss, &
1186 grid_atom, basis_1c, harmonics, nspins)
1187 END IF
1188 IF (tau_f) THEN
1189 CALL dgavtaudgb(vtau_h, vtau_s, int_hh, int_ss, &
1190 tau_basis_cache, nspins)
1191 END IF
1192
1193 NULLIFY (r_h, r_s, dr_h, dr_s)
1194
1195 END DO
1196
1197 ! some cleanup
1198 DEALLOCATE (rho_h, rho_s, rho1_h, rho1_s, vxc_h, vxc_s)
1199 DEALLOCATE (drho_h, drho_s, vxg_h, vxg_s)
1200 DEALLOCATE (drho1_h, drho1_s)
1201 IF (tau_f) THEN
1202 DEALLOCATE (tau_h, tau_s, tau1_h, tau1_s)
1203 DEALLOCATE (vtau_h, vtau_s)
1204 CALL release_tau_basis_cache(tau_basis_cache)
1205 END IF
1206
1207 CALL xc_dset_release(deriv_set)
1208 CALL xc_rho_set_release(rho_set_h)
1209 CALL xc_rho_set_release(rho1_set_h)
1210 CALL xc_rho_set_release(rho_set_s)
1211 CALL xc_rho_set_release(rho1_set_s)
1212 END DO
1213
1214 CALL timestop(handle)
1215
1216 END SUBROUTINE calculate_xc_2nd_deriv_atom
1217
1218! **************************************************************************************************
1219!> \brief ...
1220!> \param qs_env ...
1221!> \param rho0_atom_set ...
1222!> \param rho1_atom_set ...
1223!> \param rho2_atom_set ...
1224!> \param kind_set ...
1225!> \param xc_section ...
1226!> \param is_triplet ...
1227!> \param accuracy ...
1228! **************************************************************************************************
1229 SUBROUTINE calculate_gfxc_atom(qs_env, rho0_atom_set, rho1_atom_set, rho2_atom_set, &
1230 kind_set, xc_section, is_triplet, accuracy)
1231
1232 TYPE(qs_environment_type), POINTER :: qs_env
1233 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho0_atom_set, rho1_atom_set, &
1234 rho2_atom_set
1235 TYPE(qs_kind_type), DIMENSION(:), POINTER :: kind_set
1236 TYPE(section_vals_type), OPTIONAL, POINTER :: xc_section
1237 LOGICAL, INTENT(IN) :: is_triplet
1238 INTEGER, INTENT(IN) :: accuracy
1239
1240 CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_gfxc_atom'
1241 REAL(kind=dp), PARAMETER :: epsrho = 5.e-4_dp
1242
1243 INTEGER :: bo(2), handle, iat, iatom, ikind, ir, &
1244 istep, mspins, myfun, na, natom, nf, &
1245 nr, ns, nspins, nstep, num_pe
1246 INTEGER, DIMENSION(2, 3) :: bounds
1247 INTEGER, DIMENSION(:), POINTER :: atom_list
1248 LOGICAL :: accint, donlcc, gradient_f, lsd, nlcc, &
1249 paw_atom, tau_f
1250 REAL(dp) :: agr, alpha, beta, density_cut, exc_h, &
1251 exc_s, gradient_cut, oeps1, oeps2, &
1252 tau_cut
1253 REAL(dp), DIMENSION(1, 1, 1) :: tau_d
1254 REAL(dp), DIMENSION(1, 1, 1, 1) :: rho_d
1255 REAL(dp), DIMENSION(:, :), POINTER :: rho_nlcc, weight_h, weight_s
1256 REAL(dp), DIMENSION(:, :, :), POINTER :: rho0_h, rho0_s, rho1_h, rho1_s, rho_h, &
1257 rho_s, tau0_h, tau0_s, tau1_h, tau1_s, &
1258 tau_h, tau_s, vtau_h, vtau_s, vxc_h, &
1259 vxc_s
1260 REAL(dp), DIMENSION(:, :, :, :), POINTER :: drho0_h, drho0_s, drho1_h, drho1_s, &
1261 drho_h, drho_s, vxg_h, vxg_s
1262 REAL(kind=dp), DIMENSION(-4:4) :: ak, bl
1263 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1264 TYPE(dft_control_type), POINTER :: dft_control
1265 TYPE(grid_atom_type), POINTER :: grid_atom
1266 TYPE(gto_basis_set_type), POINTER :: basis_1c
1267 TYPE(harmonics_atom_type), POINTER :: harmonics
1268 TYPE(mp_para_env_type), POINTER :: para_env
1269 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s, fint_hh, fint_ss, int_hh, &
1270 int_ss, r_h, r_s
1271 TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r_h_d, r_s_d
1272 TYPE(rho_atom_type), POINTER :: rho0_atom, rho1_atom, rho2_atom
1273 TYPE(section_vals_type), POINTER :: xc_fun_section
1274 TYPE(tau_basis_cache_type) :: tau_basis_cache
1275 TYPE(xc_derivative_set_type) :: deriv_set
1276 TYPE(xc_rho_cflags_type) :: needs
1277 TYPE(xc_rho_set_type) :: rho_set_h, rho_set_s
1278
1279 CALL timeset(routinen, handle)
1280
1281 NULLIFY (vtau_h, vtau_s)
1282
1283 ak = 0.0_dp
1284 bl = 0.0_dp
1285 SELECT CASE (accuracy)
1286 CASE (:4)
1287 nstep = 2
1288 ak(-2:2) = [1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp]/12.0_dp
1289 bl(-2:2) = [-1.0_dp, 16.0_dp, -30.0_dp, 16.0_dp, -1.0_dp]/12.0_dp
1290 CASE (5:7)
1291 nstep = 3
1292 ak(-3:3) = [-1.0_dp, 9.0_dp, -45.0_dp, 0.0_dp, 45.0_dp, -9.0_dp, 1.0_dp]/60.0_dp
1293 bl(-3:3) = [2.0_dp, -27.0_dp, 270.0_dp, -490.0_dp, 270.0_dp, -27.0_dp, 2.0_dp]/180.0_dp
1294 CASE (8:)
1295 nstep = 4
1296 ak(-4:4) = [1.0_dp, -32.0_dp/3.0_dp, 56.0_dp, -224.0_dp, 0.0_dp, &
1297 224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp]/280.0_dp
1298 bl(-4:4) = [-1.0_dp, 128.0_dp/9.0_dp, -112.0_dp, 896.0_dp, -14350.0_dp/9.0_dp, &
1299 896.0_dp, -112.0_dp, 128.0_dp/9.0_dp, -1.0_dp]/560.0_dp
1300 END SELECT
1301 oeps1 = 1.0_dp/epsrho
1302 oeps2 = 1.0_dp/(epsrho**2)
1303
1304 CALL get_qs_env(qs_env=qs_env, &
1305 dft_control=dft_control, &
1306 para_env=para_env, &
1307 atomic_kind_set=atomic_kind_set)
1308
1309 xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
1310 CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
1311
1312 accint = dft_control%qs_control%gapw_control%accurate_xcint
1313
1314 IF (myfun == xc_none) THEN
1315 ! no action needed?
1316 ELSE
1317 CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", r_val=density_cut)
1318 CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", r_val=gradient_cut)
1319 CALL section_vals_val_get(xc_section, "TAU_CUTOFF", r_val=tau_cut)
1320
1321 nlcc = has_nlcc(kind_set)
1322 lsd = dft_control%lsd
1323 nspins = dft_control%nspins
1324 mspins = nspins
1325 IF (is_triplet) THEN
1326 cpassert(nspins == 1)
1327 lsd = .true.
1328 mspins = 2
1329 END IF
1330 needs = xc_functionals_get_needs(xc_fun_section, lsd=lsd, calc_potential=.true.)
1331 gradient_f = (needs%drho .OR. needs%drho_spin)
1332 tau_f = (needs%tau .OR. needs%tau_spin)
1333
1334 ! Here starts the loop over all the atoms
1335 DO ikind = 1, SIZE(atomic_kind_set)
1336 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
1337 CALL get_qs_kind(kind_set(ikind), paw_atom=paw_atom, &
1338 harmonics=harmonics, grid_atom=grid_atom)
1339 CALL get_qs_kind(kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
1340
1341 IF (.NOT. paw_atom) cycle
1342
1343 nr = grid_atom%nr
1344 na = grid_atom%ng_sphere
1345
1346 ! set integration weights
1347 IF (accint) THEN
1348 weight_h => grid_atom%weight
1349 alpha = dft_control%qs_control%gapw_control%aw(ikind)
1350 IF (ASSOCIATED(grid_atom%gapw_weight_s)) THEN
1351 IF (grid_atom%gapw_weight_alpha /= alpha) DEALLOCATE (grid_atom%gapw_weight_s)
1352 END IF
1353 IF (.NOT. ASSOCIATED(grid_atom%gapw_weight_s)) THEN
1354 ALLOCATE (grid_atom%gapw_weight_s(na, nr))
1355 DO ir = 1, nr
1356 agr = 1.0_dp - exp(-alpha*grid_atom%rad2(ir))
1357 grid_atom%gapw_weight_s(:, ir) = grid_atom%weight(:, ir)*agr
1358 END DO
1359 grid_atom%gapw_weight_alpha = alpha
1360 END IF
1361 weight_s => grid_atom%gapw_weight_s
1362 ELSE
1363 weight_h => grid_atom%weight
1364 weight_s => grid_atom%weight
1365 END IF
1366
1367 ! Prepare the structures needed to calculate and store the xc derivatives
1368
1369 ! Array dimension: here anly one dimensional arrays are used,
1370 ! i.e. only the first column of deriv_data is read.
1371 ! The other to dimensions are set to size equal 1
1372 bounds(1:2, 1:3) = 1
1373 bounds(2, 1) = na
1374 bounds(2, 2) = nr
1375
1376 ! create a place where to put the derivatives
1377 CALL xc_dset_create(deriv_set, local_bounds=bounds)
1378 ! create the place where to store the argument for the functionals
1379 CALL xc_rho_set_create(rho_set_h, bounds, rho_cutoff=density_cut, &
1380 drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
1381 CALL xc_rho_set_create(rho_set_s, bounds, rho_cutoff=density_cut, &
1382 drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
1383
1384 ! allocate the required 3d arrays where to store rho and drho
1385 CALL xc_rho_set_atom_update(rho_set_h, needs, mspins, bounds)
1386 CALL xc_rho_set_atom_update(rho_set_s, needs, mspins, bounds)
1387
1388 ALLOCATE (rho_h(na, nr, mspins), rho_s(na, nr, mspins), &
1389 rho0_h(na, nr, nspins), rho0_s(na, nr, nspins), &
1390 rho1_h(na, nr, nspins), rho1_s(na, nr, nspins))
1391 ALLOCATE (vxc_h(na, nr, mspins), vxc_s(na, nr, mspins))
1392 IF (gradient_f) THEN
1393 ALLOCATE (drho_h(4, na, nr, mspins), drho_s(4, na, nr, mspins), &
1394 drho0_h(4, na, nr, nspins), drho0_s(4, na, nr, nspins), &
1395 drho1_h(4, na, nr, nspins), drho1_s(4, na, nr, nspins))
1396 ALLOCATE (vxg_h(3, na, nr, mspins), vxg_s(3, na, nr, mspins))
1397 END IF
1398 IF (tau_f) THEN
1399 CALL create_tau_basis_cache(tau_basis_cache, grid_atom, basis_1c, harmonics)
1400 ALLOCATE (tau_h(na, nr, mspins), tau_s(na, nr, mspins), &
1401 tau0_h(na, nr, nspins), tau0_s(na, nr, nspins), &
1402 tau1_h(na, nr, nspins), tau1_s(na, nr, nspins))
1403 ALLOCATE (vtau_h(na, nr, mspins), vtau_s(na, nr, mspins))
1404 END IF
1405 !
1406 ! NLCC: prepare rho and drho of the core charge for this KIND
1407 donlcc = .false.
1408 IF (nlcc) THEN
1409 NULLIFY (rho_nlcc)
1410 rho_nlcc => kind_set(ikind)%nlcc_pot
1411 IF (ASSOCIATED(rho_nlcc)) donlcc = .true.
1412 END IF
1413
1414 ! Distribute the atoms of this kind
1415 num_pe = para_env%num_pe
1416 bo = get_limit(natom, num_pe, para_env%mepos)
1417
1418 DO iat = bo(1), bo(2)
1419 iatom = atom_list(iat)
1420 !
1421 NULLIFY (int_hh, int_ss)
1422 rho0_atom => rho0_atom_set(iatom)
1423 CALL get_rho_atom(rho_atom=rho0_atom, ga_vlocal_gb_h=int_hh, ga_vlocal_gb_s=int_ss)
1424 ALLOCATE (fint_ss(nspins), fint_hh(nspins))
1425 DO ns = 1, nspins
1426 nf = SIZE(int_ss(ns)%r_coef, 1)
1427 ALLOCATE (fint_ss(ns)%r_coef(nf, nf))
1428 nf = SIZE(int_hh(ns)%r_coef, 1)
1429 ALLOCATE (fint_hh(ns)%r_coef(nf, nf))
1430 END DO
1431
1432 ! RHO0
1433 rho0_h = 0.0_dp
1434 rho0_s = 0.0_dp
1435 rho0_atom => rho0_atom_set(iatom)
1436 IF (gradient_f) THEN
1437 NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
1438 CALL get_rho_atom(rho_atom=rho0_atom, rho_rad_h=r_h, rho_rad_s=r_s, drho_rad_h=dr_h, &
1439 drho_rad_s=dr_s, rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
1440 drho0_h = 0.0_dp
1441 drho0_s = 0.0_dp
1442 ELSE
1443 NULLIFY (r_h, r_s)
1444 CALL get_rho_atom(rho_atom=rho0_atom, rho_rad_h=r_h, rho_rad_s=r_s)
1445 rho_d = 0.0_dp
1446 END IF
1447 DO ir = 1, nr
1448 CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_f, &
1449 ir, r_h, r_s, rho0_h, rho0_s, dr_h, dr_s, &
1450 r_h_d, r_s_d, drho0_h, drho0_s)
1451 IF (donlcc) THEN
1452 CALL calc_rho_nlcc(grid_atom, nspins, gradient_f, &
1453 ir, rho_nlcc(:, 1), rho0_h, rho0_s, rho_nlcc(:, 2), drho0_h, drho0_s)
1454 END IF
1455 END DO
1456 IF (tau_f) THEN
1457 !compute tau on the grid all at once
1458 CALL calc_tau_atom(tau0_h, tau0_s, rho0_atom, tau_basis_cache, nspins)
1459 ELSE
1460 tau_d = 0.0_dp
1461 END IF
1462 ! RHO1
1463 rho1_h = 0.0_dp
1464 rho1_s = 0.0_dp
1465 rho1_atom => rho1_atom_set(iatom)
1466 IF (gradient_f) THEN
1467 NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
1468 CALL get_rho_atom(rho_atom=rho1_atom, rho_rad_h=r_h, rho_rad_s=r_s, drho_rad_h=dr_h, &
1469 drho_rad_s=dr_s, rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
1470 drho1_h = 0.0_dp
1471 drho1_s = 0.0_dp
1472 ELSE
1473 NULLIFY (r_h, r_s)
1474 CALL get_rho_atom(rho_atom=rho1_atom, rho_rad_h=r_h, rho_rad_s=r_s)
1475 END IF
1476 DO ir = 1, nr
1477 CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_f, &
1478 ir, r_h, r_s, rho1_h, rho1_s, dr_h, dr_s, &
1479 r_h_d, r_s_d, drho1_h, drho1_s)
1480 END DO
1481 IF (tau_f) THEN
1482 !compute tau on the grid all at once
1483 CALL calc_tau_atom(tau1_h, tau1_s, rho1_atom, tau_basis_cache, nspins)
1484 END IF
1485 ! RHO2
1486 rho2_atom => rho2_atom_set(iatom)
1487
1488 DO istep = -nstep, nstep
1489
1490 beta = real(istep, kind=dp)*epsrho
1491
1492 IF (is_triplet) THEN
1493 rho_h(:, :, 1) = rho0_h(:, :, 1) + beta*rho1_h(:, :, 1)
1494 rho_h(:, :, 2) = rho0_h(:, :, 1)
1495 rho_h = 0.5_dp*rho_h
1496 rho_s(:, :, 1) = rho0_s(:, :, 1) + beta*rho1_s(:, :, 1)
1497 rho_s(:, :, 2) = rho0_s(:, :, 1)
1498 rho_s = 0.5_dp*rho_s
1499 IF (gradient_f) THEN
1500 drho_h(:, :, :, 1) = drho0_h(:, :, :, 1) + beta*drho1_h(:, :, :, 1)
1501 drho_h(:, :, :, 2) = drho0_h(:, :, :, 1)
1502 drho_h = 0.5_dp*drho_h
1503 drho_s(:, :, :, 1) = drho0_s(:, :, :, 1) + beta*drho1_s(:, :, :, 1)
1504 drho_s(:, :, :, 2) = drho0_s(:, :, :, 1)
1505 drho_s = 0.5_dp*drho_s
1506 END IF
1507 IF (tau_f) THEN
1508 tau_h(:, :, 1) = tau0_h(:, :, 1) + beta*tau1_h(:, :, 1)
1509 tau_h(:, :, 2) = tau0_h(:, :, 1)
1510 tau_h = 0.5_dp*tau0_h
1511 tau_s(:, :, 1) = tau0_s(:, :, 1) + beta*tau1_s(:, :, 1)
1512 tau_s(:, :, 2) = tau0_s(:, :, 1)
1513 tau_s = 0.5_dp*tau0_s
1514 END IF
1515 ELSE
1516 rho_h = rho0_h + beta*rho1_h
1517 rho_s = rho0_s + beta*rho1_s
1518 IF (gradient_f) THEN
1519 drho_h = drho0_h + beta*drho1_h
1520 drho_s = drho0_s + beta*drho1_s
1521 END IF
1522 IF (tau_f) THEN
1523 tau_h = tau0_h + beta*tau1_h
1524 tau_s = tau0_s + beta*tau1_s
1525 END IF
1526 END IF
1527 !
1528 IF (gradient_f) THEN
1529 drho_h(4, :, :, :) = sqrt( &
1530 drho_h(1, :, :, :)*drho_h(1, :, :, :) + &
1531 drho_h(2, :, :, :)*drho_h(2, :, :, :) + &
1532 drho_h(3, :, :, :)*drho_h(3, :, :, :))
1533
1534 drho_s(4, :, :, :) = sqrt( &
1535 drho_s(1, :, :, :)*drho_s(1, :, :, :) + &
1536 drho_s(2, :, :, :)*drho_s(2, :, :, :) + &
1537 drho_s(3, :, :, :)*drho_s(3, :, :, :))
1538 END IF
1539
1540 DO ir = 1, nr
1541 IF (tau_f) THEN
1542 CALL fill_rho_set(rho_set_h, lsd, mspins, needs, rho_h, drho_h, tau_h, na, ir)
1543 CALL fill_rho_set(rho_set_s, lsd, mspins, needs, rho_s, drho_s, tau_s, na, ir)
1544 ELSE IF (gradient_f) THEN
1545 CALL fill_rho_set(rho_set_h, lsd, mspins, needs, rho_h, drho_h, tau_d, na, ir)
1546 CALL fill_rho_set(rho_set_s, lsd, mspins, needs, rho_s, drho_s, tau_d, na, ir)
1547 ELSE
1548 CALL fill_rho_set(rho_set_h, lsd, mspins, needs, rho_h, rho_d, tau_d, na, ir)
1549 CALL fill_rho_set(rho_set_s, lsd, mspins, needs, rho_s, rho_d, tau_d, na, ir)
1550 END IF
1551 END DO
1552
1553 ! hard atom density !
1554 CALL xc_dset_zero_all(deriv_set)
1555 CALL vxc_of_r_new(xc_fun_section, rho_set_h, deriv_set, 1, needs, weight_h, &
1556 lsd, na, nr, exc_h, vxc_h, vxg_h, vtau_h)
1557 IF (is_triplet) THEN
1558 vxc_h(:, :, 1) = vxc_h(:, :, 1) - vxc_h(:, :, 2)
1559 IF (gradient_f) THEN
1560 vxg_h(:, :, :, 1) = vxg_h(:, :, :, 1) - vxg_h(:, :, :, 2)
1561 END IF
1562 IF (tau_f) THEN
1563 vtau_h(:, :, 1) = vtau_h(:, :, 1) - vtau_h(:, :, 2)
1564 END IF
1565 END IF
1566 ! soft atom density !
1567 CALL xc_dset_zero_all(deriv_set)
1568 CALL vxc_of_r_new(xc_fun_section, rho_set_s, deriv_set, 1, needs, weight_s, &
1569 lsd, na, nr, exc_s, vxc_s, vxg_s, vtau_s)
1570 IF (is_triplet) THEN
1571 vxc_s(:, :, 1) = vxc_s(:, :, 1) - vxc_s(:, :, 2)
1572 IF (gradient_f) THEN
1573 vxg_s(:, :, :, 1) = vxg_s(:, :, :, 1) - vxg_s(:, :, :, 2)
1574 END IF
1575 IF (tau_f) THEN
1576 vtau_s(:, :, 1) = vtau_s(:, :, 1) - vtau_s(:, :, 2)
1577 END IF
1578 END IF
1579 ! potentials
1580 DO ns = 1, nspins
1581 fint_hh(ns)%r_coef(:, :) = 0.0_dp
1582 fint_ss(ns)%r_coef(:, :) = 0.0_dp
1583 END DO
1584 IF (gradient_f) THEN
1585 CALL gavxcgb_gc(vxc_h, vxc_s, vxg_h, vxg_s, fint_hh, fint_ss, &
1586 grid_atom, basis_1c, harmonics, nspins)
1587 ELSE
1588 CALL gavxcgb_nogc(vxc_h, vxc_s, fint_hh, fint_ss, &
1589 grid_atom, basis_1c, harmonics, nspins)
1590 END IF
1591 IF (tau_f) THEN
1592 CALL dgavtaudgb(vtau_h, vtau_s, fint_hh, fint_ss, &
1593 tau_basis_cache, nspins)
1594 END IF
1595 ! first derivative fxc
1596 NULLIFY (int_hh, int_ss)
1597 CALL get_rho_atom(rho_atom=rho1_atom, ga_vlocal_gb_h=int_hh, ga_vlocal_gb_s=int_ss)
1598 DO ns = 1, nspins
1599 int_ss(ns)%r_coef(:, :) = int_ss(ns)%r_coef(:, :) + oeps1*ak(istep)*fint_ss(ns)%r_coef(:, :)
1600 int_hh(ns)%r_coef(:, :) = int_hh(ns)%r_coef(:, :) + oeps1*ak(istep)*fint_hh(ns)%r_coef(:, :)
1601 END DO
1602 ! second derivative gxc
1603 NULLIFY (int_hh, int_ss)
1604 CALL get_rho_atom(rho_atom=rho2_atom, ga_vlocal_gb_h=int_hh, ga_vlocal_gb_s=int_ss)
1605 DO ns = 1, nspins
1606 int_ss(ns)%r_coef(:, :) = int_ss(ns)%r_coef(:, :) + oeps2*bl(istep)*fint_ss(ns)%r_coef(:, :)
1607 int_hh(ns)%r_coef(:, :) = int_hh(ns)%r_coef(:, :) + oeps2*bl(istep)*fint_hh(ns)%r_coef(:, :)
1608 END DO
1609 END DO
1610 !
1611 DO ns = 1, nspins
1612 DEALLOCATE (fint_ss(ns)%r_coef)
1613 DEALLOCATE (fint_hh(ns)%r_coef)
1614 END DO
1615 DEALLOCATE (fint_ss, fint_hh)
1616
1617 END DO ! iat
1618
1619 ! Release the xc structure used to store the xc derivatives
1620 CALL xc_dset_release(deriv_set)
1621 CALL xc_rho_set_release(rho_set_h)
1622 CALL xc_rho_set_release(rho_set_s)
1623
1624 DEALLOCATE (rho_h, rho_s, rho0_h, rho0_s, rho1_h, rho1_s)
1625 DEALLOCATE (vxc_h, vxc_s)
1626 IF (gradient_f) THEN
1627 DEALLOCATE (drho_h, drho_s, drho0_h, drho0_s, drho1_h, drho1_s)
1628 DEALLOCATE (vxg_h, vxg_s)
1629 END IF
1630 IF (tau_f) THEN
1631 DEALLOCATE (tau_h, tau_s, tau0_h, tau0_s, tau1_h, tau1_s)
1632 DEALLOCATE (vtau_h, vtau_s)
1633 CALL release_tau_basis_cache(tau_basis_cache)
1634 END IF
1635 END DO ! ikind
1636
1637 END IF !xc_none
1638
1639 CALL timestop(handle)
1640
1641 END SUBROUTINE calculate_gfxc_atom
1642
1643! **************************************************************************************************
1644!> \brief ...
1645!> \param qs_env ...
1646!> \param rho0_atom_set ...
1647!> \param rho1_atom_set ...
1648!> \param rho2_atom_set ...
1649!> \param kind_set ...
1650!> \param xc_section ...
1651!> \param is_triplet ...
1652!> \param accuracy ...
1653!> \param epsrho ...
1654! **************************************************************************************************
1655 SUBROUTINE gfxc_atom_diff(qs_env, rho0_atom_set, rho1_atom_set, rho2_atom_set, &
1656 kind_set, xc_section, is_triplet, accuracy, epsrho)
1657
1658 TYPE(qs_environment_type), POINTER :: qs_env
1659 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho0_atom_set, rho1_atom_set, &
1660 rho2_atom_set
1661 TYPE(qs_kind_type), DIMENSION(:), POINTER :: kind_set
1662 TYPE(section_vals_type), OPTIONAL, POINTER :: xc_section
1663 LOGICAL, INTENT(IN) :: is_triplet
1664 INTEGER, INTENT(IN) :: accuracy
1665 REAL(kind=dp), INTENT(IN) :: epsrho
1666
1667 CHARACTER(LEN=*), PARAMETER :: routinen = 'gfxc_atom_diff'
1668
1669 INTEGER :: bo(2), handle, iat, iatom, ikind, ir, &
1670 istep, mspins, myfun, na, natom, nf, &
1671 nr, ns, nspins, nstep, num_pe
1672 INTEGER, DIMENSION(2, 3) :: bounds
1673 INTEGER, DIMENSION(:), POINTER :: atom_list
1674 LOGICAL :: accint, donlcc, gradient_f, lsd, nlcc, &
1675 paw_atom, tau_f
1676 REAL(dp) :: agr, alpha, beta, density_cut, &
1677 gradient_cut, oeps1, tau_cut
1678 REAL(dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: vtau_h, vtau_s, vxc_h, vxc_s
1679 REAL(dp), DIMENSION(1, 1, 1) :: tau_d
1680 REAL(dp), DIMENSION(1, 1, 1, 1) :: rho_d
1681 REAL(dp), DIMENSION(:, :), POINTER :: rho_nlcc, weight_h, weight_s
1682 REAL(dp), DIMENSION(:, :, :), POINTER :: rho0_h, rho0_s, rho1_h, rho1_s, rho_h, &
1683 rho_s, tau0_h, tau0_s, tau1_h, tau1_s, &
1684 tau_h, tau_s
1685 REAL(dp), DIMENSION(:, :, :, :), POINTER :: drho0_h, drho0_s, drho1_h, drho1_s, &
1686 drho_h, drho_s, vxg_h, vxg_s
1687 REAL(kind=dp), DIMENSION(-4:4) :: ak
1688 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1689 TYPE(dft_control_type), POINTER :: dft_control
1690 TYPE(grid_atom_type), POINTER :: grid_atom
1691 TYPE(gto_basis_set_type), POINTER :: basis_1c
1692 TYPE(harmonics_atom_type), POINTER :: harmonics
1693 TYPE(mp_para_env_type), POINTER :: para_env
1694 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s, fint_hh, fint_ss, int_hh, &
1695 int_ss, r_h, r_s
1696 TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r_h_d, r_s_d
1697 TYPE(rho_atom_type), POINTER :: rho0_atom, rho1_atom, rho2_atom
1698 TYPE(section_vals_type), POINTER :: xc_fun_section
1699 TYPE(tau_basis_cache_type) :: tau_basis_cache
1700 TYPE(xc_derivative_set_type) :: deriv_set
1701 TYPE(xc_rho_cflags_type) :: needs
1702 TYPE(xc_rho_set_type) :: rho1_set_h, rho1_set_s, rho_set_h, &
1703 rho_set_s
1704
1705 CALL timeset(routinen, handle)
1706
1707 NULLIFY (vtau_h, vtau_s)
1708
1709 ak = 0.0_dp
1710 SELECT CASE (accuracy)
1711 CASE (:4)
1712 nstep = 2
1713 ak(-2:2) = [1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp]/12.0_dp
1714 CASE (5:7)
1715 nstep = 3
1716 ak(-3:3) = [-1.0_dp, 9.0_dp, -45.0_dp, 0.0_dp, 45.0_dp, -9.0_dp, 1.0_dp]/60.0_dp
1717 CASE (8:)
1718 nstep = 4
1719 ak(-4:4) = [1.0_dp, -32.0_dp/3.0_dp, 56.0_dp, -224.0_dp, 0.0_dp, &
1720 224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp]/280.0_dp
1721 END SELECT
1722 oeps1 = 1.0_dp/epsrho
1723
1724 CALL get_qs_env(qs_env=qs_env, &
1725 dft_control=dft_control, &
1726 para_env=para_env, &
1727 atomic_kind_set=atomic_kind_set)
1728
1729 xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
1730 CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
1731
1732 accint = dft_control%qs_control%gapw_control%accurate_xcint
1733
1734 IF (myfun == xc_none) THEN
1735 ! no action needed?
1736 ELSE
1737 ! calculate fxc
1738 CALL calculate_xc_2nd_deriv_atom(rho0_atom_set, rho1_atom_set, qs_env, xc_section, para_env, &
1739 do_triplet=is_triplet, kind_set_external=kind_set)
1740
1741 CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", r_val=density_cut)
1742 CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", r_val=gradient_cut)
1743 CALL section_vals_val_get(xc_section, "TAU_CUTOFF", r_val=tau_cut)
1744
1745 nlcc = has_nlcc(kind_set)
1746 lsd = dft_control%lsd
1747 nspins = dft_control%nspins
1748 mspins = nspins
1749 IF (is_triplet) THEN
1750 cpassert(nspins == 1)
1751 lsd = .true.
1752 mspins = 2
1753 END IF
1754 needs = xc_functionals_get_needs(xc_fun_section, lsd=lsd, calc_potential=.true.)
1755 gradient_f = (needs%drho .OR. needs%drho_spin)
1756 tau_f = (needs%tau .OR. needs%tau_spin)
1757
1758 ! Here starts the loop over all the atoms
1759 DO ikind = 1, SIZE(atomic_kind_set)
1760 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
1761 CALL get_qs_kind(kind_set(ikind), paw_atom=paw_atom, &
1762 harmonics=harmonics, grid_atom=grid_atom)
1763 CALL get_qs_kind(kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
1764
1765 IF (.NOT. paw_atom) cycle
1766
1767 nr = grid_atom%nr
1768 na = grid_atom%ng_sphere
1769
1770 ! set integration weights
1771 IF (accint) THEN
1772 weight_h => grid_atom%weight
1773 alpha = dft_control%qs_control%gapw_control%aw(ikind)
1774 IF (ASSOCIATED(grid_atom%gapw_weight_s)) THEN
1775 IF (grid_atom%gapw_weight_alpha /= alpha) DEALLOCATE (grid_atom%gapw_weight_s)
1776 END IF
1777 IF (.NOT. ASSOCIATED(grid_atom%gapw_weight_s)) THEN
1778 ALLOCATE (grid_atom%gapw_weight_s(na, nr))
1779 DO ir = 1, nr
1780 agr = 1.0_dp - exp(-alpha*grid_atom%rad2(ir))
1781 grid_atom%gapw_weight_s(:, ir) = grid_atom%weight(:, ir)*agr
1782 END DO
1783 grid_atom%gapw_weight_alpha = alpha
1784 END IF
1785 weight_s => grid_atom%gapw_weight_s
1786 ELSE
1787 weight_h => grid_atom%weight
1788 weight_s => grid_atom%weight
1789 END IF
1790
1791 ! Prepare the structures needed to calculate and store the xc derivatives
1792
1793 ! Array dimension: here anly one dimensional arrays are used,
1794 ! i.e. only the first column of deriv_data is read.
1795 ! The other to dimensions are set to size equal 1
1796 bounds(1:2, 1:3) = 1
1797 bounds(2, 1) = na
1798 bounds(2, 2) = nr
1799
1800 ! create a place where to put the derivatives
1801 CALL xc_dset_create(deriv_set, local_bounds=bounds)
1802 ! create the place where to store the argument for the functionals
1803 CALL xc_rho_set_create(rho_set_h, bounds, rho_cutoff=density_cut, &
1804 drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
1805 CALL xc_rho_set_create(rho_set_s, bounds, rho_cutoff=density_cut, &
1806 drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
1807 CALL xc_rho_set_create(rho1_set_h, bounds, rho_cutoff=density_cut, &
1808 drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
1809 CALL xc_rho_set_create(rho1_set_s, bounds, rho_cutoff=density_cut, &
1810 drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
1811
1812 ! allocate the required 3d arrays where to store rho and drho
1813 CALL xc_rho_set_atom_update(rho_set_h, needs, mspins, bounds)
1814 CALL xc_rho_set_atom_update(rho_set_s, needs, mspins, bounds)
1815 CALL xc_rho_set_atom_update(rho1_set_h, needs, mspins, bounds)
1816 CALL xc_rho_set_atom_update(rho1_set_s, needs, mspins, bounds)
1817
1818 ALLOCATE (rho_h(na, nr, nspins), rho_s(na, nr, nspins), &
1819 rho0_h(na, nr, nspins), rho0_s(na, nr, nspins), &
1820 rho1_h(na, nr, nspins), rho1_s(na, nr, nspins))
1821 ALLOCATE (vxc_h(na, nr, nspins), vxc_s(na, nr, nspins))
1822 IF (gradient_f) THEN
1823 ALLOCATE (drho_h(4, na, nr, nspins), drho_s(4, na, nr, nspins), &
1824 drho0_h(4, na, nr, nspins), drho0_s(4, na, nr, nspins), &
1825 drho1_h(4, na, nr, nspins), drho1_s(4, na, nr, nspins))
1826 ALLOCATE (vxg_h(3, na, nr, nspins), vxg_s(3, na, nr, nspins))
1827 END IF
1828 IF (tau_f) THEN
1829 CALL create_tau_basis_cache(tau_basis_cache, grid_atom, basis_1c, harmonics)
1830 ALLOCATE (tau_h(na, nr, nspins), tau_s(na, nr, nspins), &
1831 tau0_h(na, nr, nspins), tau0_s(na, nr, nspins), &
1832 tau1_h(na, nr, nspins), tau1_s(na, nr, nspins))
1833 ALLOCATE (vtau_h(na, nr, nspins), vtau_s(na, nr, nspins))
1834 END IF
1835 !
1836 ! NLCC: prepare rho and drho of the core charge for this KIND
1837 donlcc = .false.
1838 IF (nlcc) THEN
1839 NULLIFY (rho_nlcc)
1840 rho_nlcc => kind_set(ikind)%nlcc_pot
1841 IF (ASSOCIATED(rho_nlcc)) donlcc = .true.
1842 END IF
1843
1844 ! Distribute the atoms of this kind
1845 num_pe = para_env%num_pe
1846 bo = get_limit(natom, num_pe, para_env%mepos)
1847
1848 DO iat = bo(1), bo(2)
1849 iatom = atom_list(iat)
1850 !
1851 NULLIFY (int_hh, int_ss)
1852 rho0_atom => rho0_atom_set(iatom)
1853 CALL get_rho_atom(rho_atom=rho0_atom, ga_vlocal_gb_h=int_hh, ga_vlocal_gb_s=int_ss)
1854 ALLOCATE (fint_ss(nspins), fint_hh(nspins))
1855 DO ns = 1, nspins
1856 nf = SIZE(int_ss(ns)%r_coef, 1)
1857 ALLOCATE (fint_ss(ns)%r_coef(nf, nf))
1858 nf = SIZE(int_hh(ns)%r_coef, 1)
1859 ALLOCATE (fint_hh(ns)%r_coef(nf, nf))
1860 END DO
1861
1862 ! RHO0
1863 rho0_h = 0.0_dp
1864 rho0_s = 0.0_dp
1865 rho0_atom => rho0_atom_set(iatom)
1866 IF (gradient_f) THEN
1867 NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
1868 CALL get_rho_atom(rho_atom=rho0_atom, rho_rad_h=r_h, rho_rad_s=r_s, drho_rad_h=dr_h, &
1869 drho_rad_s=dr_s, rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
1870 drho0_h = 0.0_dp
1871 drho0_s = 0.0_dp
1872 ELSE
1873 NULLIFY (r_h, r_s)
1874 CALL get_rho_atom(rho_atom=rho0_atom, rho_rad_h=r_h, rho_rad_s=r_s)
1875 rho_d = 0.0_dp
1876 END IF
1877 DO ir = 1, nr
1878 CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_f, &
1879 ir, r_h, r_s, rho0_h, rho0_s, dr_h, dr_s, &
1880 r_h_d, r_s_d, drho0_h, drho0_s)
1881 IF (donlcc) THEN
1882 CALL calc_rho_nlcc(grid_atom, nspins, gradient_f, &
1883 ir, rho_nlcc(:, 1), rho0_h, rho0_s, rho_nlcc(:, 2), drho0_h, drho0_s)
1884 END IF
1885 END DO
1886 IF (tau_f) THEN
1887 !compute tau on the grid all at once
1888 CALL calc_tau_atom(tau0_h, tau0_s, rho0_atom, tau_basis_cache, nspins)
1889 ELSE
1890 tau_d = 0.0_dp
1891 END IF
1892 ! RHO1
1893 rho1_h = 0.0_dp
1894 rho1_s = 0.0_dp
1895 rho1_atom => rho1_atom_set(iatom)
1896 IF (gradient_f) THEN
1897 NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
1898 CALL get_rho_atom(rho_atom=rho1_atom, rho_rad_h=r_h, rho_rad_s=r_s, drho_rad_h=dr_h, &
1899 drho_rad_s=dr_s, rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
1900 drho1_h = 0.0_dp
1901 drho1_s = 0.0_dp
1902 ELSE
1903 NULLIFY (r_h, r_s)
1904 CALL get_rho_atom(rho_atom=rho1_atom, rho_rad_h=r_h, rho_rad_s=r_s)
1905 END IF
1906 DO ir = 1, nr
1907 CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_f, &
1908 ir, r_h, r_s, rho1_h, rho1_s, dr_h, dr_s, &
1909 r_h_d, r_s_d, drho1_h, drho1_s)
1910 END DO
1911 IF (tau_f) THEN
1912 !compute tau on the grid all at once
1913 CALL calc_tau_atom(tau1_h, tau1_s, rho1_atom, tau_basis_cache, nspins)
1914 END IF
1915
1916 DO ir = 1, nr
1917 IF (tau_f) THEN
1918 CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, drho1_h, tau1_h, na, ir)
1919 CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, drho1_s, tau1_s, na, ir)
1920 ELSE IF (gradient_f) THEN
1921 CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, drho1_h, tau_d, na, ir)
1922 CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, drho1_s, tau_d, na, ir)
1923 ELSE
1924 CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, rho_d, tau_d, na, ir)
1925 CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, rho_d, tau_d, na, ir)
1926 END IF
1927 END DO
1928
1929 ! RHO2
1930 rho2_atom => rho2_atom_set(iatom)
1931
1932 DO istep = -nstep, nstep
1933
1934 beta = real(istep, kind=dp)*epsrho
1935
1936 rho_h = rho0_h + beta*rho1_h
1937 rho_s = rho0_s + beta*rho1_s
1938 IF (gradient_f) THEN
1939 drho_h = drho0_h + beta*drho1_h
1940 drho_s = drho0_s + beta*drho1_s
1941 END IF
1942 IF (tau_f) THEN
1943 tau_h = tau0_h + beta*tau1_h
1944 tau_s = tau0_s + beta*tau1_s
1945 END IF
1946 !
1947 IF (gradient_f) THEN
1948 drho_h(4, :, :, :) = sqrt( &
1949 drho_h(1, :, :, :)*drho_h(1, :, :, :) + &
1950 drho_h(2, :, :, :)*drho_h(2, :, :, :) + &
1951 drho_h(3, :, :, :)*drho_h(3, :, :, :))
1952
1953 drho_s(4, :, :, :) = sqrt( &
1954 drho_s(1, :, :, :)*drho_s(1, :, :, :) + &
1955 drho_s(2, :, :, :)*drho_s(2, :, :, :) + &
1956 drho_s(3, :, :, :)*drho_s(3, :, :, :))
1957 END IF
1958
1959 DO ir = 1, nr
1960 IF (tau_f) THEN
1961 CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_h, na, ir)
1962 CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_s, na, ir)
1963 ELSE IF (gradient_f) THEN
1964 CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_d, na, ir)
1965 CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_d, na, ir)
1966 ELSE
1967 CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, rho_d, tau_d, na, ir)
1968 CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, rho_d, tau_d, na, ir)
1969 END IF
1970 END DO
1971
1972 ! hard atom density !
1973 CALL xc_dset_zero_all(deriv_set)
1974 CALL xc_2nd_deriv_of_r(xc_section=xc_section, &
1975 rho_set=rho_set_h, rho1_set=rho1_set_h, &
1976 deriv_set=deriv_set, &
1977 w=weight_h, vxc=vxc_h, vxg=vxg_h, vtau=vtau_h, &
1978 do_triplet=is_triplet)
1979 ! soft atom density !
1980 CALL xc_dset_zero_all(deriv_set)
1981 CALL xc_2nd_deriv_of_r(xc_section=xc_section, &
1982 rho_set=rho_set_s, rho1_set=rho1_set_s, &
1983 deriv_set=deriv_set, &
1984 w=weight_s, vxc=vxc_s, vxg=vxg_s, vtau=vtau_s, &
1985 do_triplet=is_triplet)
1986 ! potentials
1987 DO ns = 1, nspins
1988 fint_hh(ns)%r_coef(:, :) = 0.0_dp
1989 fint_ss(ns)%r_coef(:, :) = 0.0_dp
1990 END DO
1991 IF (gradient_f) THEN
1992 CALL gavxcgb_gc(vxc_h, vxc_s, vxg_h, vxg_s, fint_hh, fint_ss, &
1993 grid_atom, basis_1c, harmonics, nspins)
1994 ELSE
1995 CALL gavxcgb_nogc(vxc_h, vxc_s, fint_hh, fint_ss, &
1996 grid_atom, basis_1c, harmonics, nspins)
1997 END IF
1998 IF (tau_f) THEN
1999 CALL dgavtaudgb(vtau_h, vtau_s, fint_hh, fint_ss, &
2000 tau_basis_cache, nspins)
2001 END IF
2002 ! second derivative gxc
2003 NULLIFY (int_hh, int_ss)
2004 CALL get_rho_atom(rho_atom=rho2_atom, ga_vlocal_gb_h=int_hh, ga_vlocal_gb_s=int_ss)
2005 DO ns = 1, nspins
2006 int_ss(ns)%r_coef(:, :) = int_ss(ns)%r_coef(:, :) + oeps1*ak(istep)*fint_ss(ns)%r_coef(:, :)
2007 int_hh(ns)%r_coef(:, :) = int_hh(ns)%r_coef(:, :) + oeps1*ak(istep)*fint_hh(ns)%r_coef(:, :)
2008 END DO
2009 END DO
2010 !
2011 DO ns = 1, nspins
2012 DEALLOCATE (fint_ss(ns)%r_coef)
2013 DEALLOCATE (fint_hh(ns)%r_coef)
2014 END DO
2015 DEALLOCATE (fint_ss, fint_hh)
2016
2017 END DO ! iat
2018
2019 ! Release the xc structure used to store the xc derivatives
2020 CALL xc_dset_release(deriv_set)
2021 CALL xc_rho_set_release(rho_set_h)
2022 CALL xc_rho_set_release(rho_set_s)
2023 CALL xc_rho_set_release(rho1_set_h)
2024 CALL xc_rho_set_release(rho1_set_s)
2025
2026 DEALLOCATE (rho_h, rho_s, rho0_h, rho0_s, rho1_h, rho1_s)
2027 DEALLOCATE (vxc_h, vxc_s)
2028 IF (gradient_f) THEN
2029 DEALLOCATE (drho_h, drho_s, drho0_h, drho0_s, drho1_h, drho1_s)
2030 DEALLOCATE (vxg_h, vxg_s)
2031 END IF
2032 IF (tau_f) THEN
2033 DEALLOCATE (tau_h, tau_s, tau0_h, tau0_s, tau1_h, tau1_s)
2034 DEALLOCATE (vtau_h, vtau_s)
2035 CALL release_tau_basis_cache(tau_basis_cache)
2036 END IF
2037 END DO ! ikind
2038
2039 END IF !xc_none
2040
2041 CALL timestop(handle)
2042
2043 END SUBROUTINE gfxc_atom_diff
2044
2045! **************************************************************************************************
2046!> \brief ...
2047!> \param grid_atom ...
2048!> \param harmonics ...
2049!> \param nspins ...
2050!> \param grad_func ...
2051!> \param ir ...
2052!> \param r_h ...
2053!> \param r_s ...
2054!> \param rho_h ...
2055!> \param rho_s ...
2056!> \param dr_h ...
2057!> \param dr_s ...
2058!> \param r_h_d ...
2059!> \param r_s_d ...
2060!> \param drho_h ...
2061!> \param drho_s ...
2062! **************************************************************************************************
2063 SUBROUTINE calc_rho_angular(grid_atom, harmonics, nspins, grad_func, &
2064 ir, r_h, r_s, rho_h, rho_s, &
2065 dr_h, dr_s, r_h_d, r_s_d, drho_h, drho_s)
2066
2067 TYPE(grid_atom_type), POINTER :: grid_atom
2068 TYPE(harmonics_atom_type), POINTER :: harmonics
2069 INTEGER, INTENT(IN) :: nspins
2070 LOGICAL, INTENT(IN) :: grad_func
2071 INTEGER, INTENT(IN) :: ir
2072 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: r_h, r_s
2073 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: rho_h, rho_s
2074 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s
2075 TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r_h_d, r_s_d
2076 REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: drho_h, drho_s
2077
2078 INTEGER :: ia, iso, ispin, na
2079 REAL(kind=dp) :: rad, urad
2080
2081 cpassert(ASSOCIATED(r_h))
2082 cpassert(ASSOCIATED(r_s))
2083 cpassert(ASSOCIATED(rho_h))
2084 cpassert(ASSOCIATED(rho_s))
2085 IF (grad_func) THEN
2086 cpassert(ASSOCIATED(dr_h))
2087 cpassert(ASSOCIATED(dr_s))
2088 cpassert(ASSOCIATED(r_h_d))
2089 cpassert(ASSOCIATED(r_s_d))
2090 cpassert(ASSOCIATED(drho_h))
2091 cpassert(ASSOCIATED(drho_s))
2092 END IF
2093
2094 na = grid_atom%ng_sphere
2095 rad = grid_atom%rad(ir)
2096 urad = grid_atom%oorad2l(ir, 1)
2097 DO ispin = 1, nspins
2098 DO iso = 1, harmonics%max_iso_not0
2099 DO ia = 1, na
2100 rho_h(ia, ir, ispin) = rho_h(ia, ir, ispin) + &
2101 r_h(ispin)%r_coef(ir, iso)*harmonics%slm(ia, iso)
2102 rho_s(ia, ir, ispin) = rho_s(ia, ir, ispin) + &
2103 r_s(ispin)%r_coef(ir, iso)*harmonics%slm(ia, iso)
2104 END DO ! ia
2105 END DO ! iso
2106 END DO ! ispin
2107
2108 IF (grad_func) THEN
2109 DO ispin = 1, nspins
2110 DO iso = 1, harmonics%max_iso_not0
2111 DO ia = 1, na
2112
2113 ! components of the gradient of rho1 hard
2114 drho_h(1, ia, ir, ispin) = drho_h(1, ia, ir, ispin) + &
2115 dr_h(ispin)%r_coef(ir, iso)* &
2116 harmonics%a(1, ia)*harmonics%slm(ia, iso) + &
2117 r_h_d(1, ispin)%r_coef(ir, iso)* &
2118 harmonics%slm(ia, iso)
2119
2120 drho_h(2, ia, ir, ispin) = drho_h(2, ia, ir, ispin) + &
2121 dr_h(ispin)%r_coef(ir, iso)* &
2122 harmonics%a(2, ia)*harmonics%slm(ia, iso) + &
2123 r_h_d(2, ispin)%r_coef(ir, iso)* &
2124 harmonics%slm(ia, iso)
2125
2126 drho_h(3, ia, ir, ispin) = drho_h(3, ia, ir, ispin) + &
2127 dr_h(ispin)%r_coef(ir, iso)* &
2128 harmonics%a(3, ia)*harmonics%slm(ia, iso) + &
2129 r_h_d(3, ispin)%r_coef(ir, iso)* &
2130 harmonics%slm(ia, iso)
2131
2132 ! components of the gradient of rho1 soft
2133 drho_s(1, ia, ir, ispin) = drho_s(1, ia, ir, ispin) + &
2134 dr_s(ispin)%r_coef(ir, iso)* &
2135 harmonics%a(1, ia)*harmonics%slm(ia, iso) + &
2136 r_s_d(1, ispin)%r_coef(ir, iso)* &
2137 harmonics%slm(ia, iso)
2138
2139 drho_s(2, ia, ir, ispin) = drho_s(2, ia, ir, ispin) + &
2140 dr_s(ispin)%r_coef(ir, iso)* &
2141 harmonics%a(2, ia)*harmonics%slm(ia, iso) + &
2142 r_s_d(2, ispin)%r_coef(ir, iso)* &
2143 harmonics%slm(ia, iso)
2144
2145 drho_s(3, ia, ir, ispin) = drho_s(3, ia, ir, ispin) + &
2146 dr_s(ispin)%r_coef(ir, iso)* &
2147 harmonics%a(3, ia)*harmonics%slm(ia, iso) + &
2148 r_s_d(3, ispin)%r_coef(ir, iso)* &
2149 harmonics%slm(ia, iso)
2150
2151 END DO ! ia
2152 END DO ! iso
2153 DO ia = 1, na
2154 drho_h(4, ia, ir, ispin) = sqrt( &
2155 drho_h(1, ia, ir, ispin)*drho_h(1, ia, ir, ispin) + &
2156 drho_h(2, ia, ir, ispin)*drho_h(2, ia, ir, ispin) + &
2157 drho_h(3, ia, ir, ispin)*drho_h(3, ia, ir, ispin))
2158
2159 drho_s(4, ia, ir, ispin) = sqrt( &
2160 drho_s(1, ia, ir, ispin)*drho_s(1, ia, ir, ispin) + &
2161 drho_s(2, ia, ir, ispin)*drho_s(2, ia, ir, ispin) + &
2162 drho_s(3, ia, ir, ispin)*drho_s(3, ia, ir, ispin))
2163 END DO ! ia
2164 END DO ! ispin
2165 END IF
2166
2167 END SUBROUTINE calc_rho_angular
2168
2169! **************************************************************************************************
2170!> \brief Precompute radial and angular factors for GAPW meta-GGA tau contractions
2171!> \param tau_cache precomputed compact one-center gradient basis
2172!> \param grid_atom atom-centered integration grid
2173!> \param basis_1c GAPW one-center basis
2174!> \param harmonics spherical harmonics on the atom-centered grid
2175! **************************************************************************************************
2176 SUBROUTINE create_tau_basis_cache(tau_cache, grid_atom, basis_1c, harmonics)
2177
2178 TYPE(tau_basis_cache_type), INTENT(INOUT) :: tau_cache
2179 TYPE(grid_atom_type), POINTER :: grid_atom
2180 TYPE(gto_basis_set_type), POINTER :: basis_1c
2181 TYPE(harmonics_atom_type), POINTER :: harmonics
2182
2183 INTEGER :: dir, ia, igrid, ip, ipgf, ir, iset, iso, &
2184 l, starti
2185 REAL(dp), ALLOCATABLE, DIMENSION(:) :: a1, a2, gexp, r1, r2
2186 REAL(dp), DIMENSION(:, :), POINTER :: slm
2187 REAL(dp), DIMENSION(:, :, :), POINTER :: dslm_dxyz
2188
2189 NULLIFY (slm, dslm_dxyz)
2190
2191 CALL release_tau_basis_cache(tau_cache)
2192
2193 CALL get_gto_basis_set(gto_basis_set=basis_1c, lmax=tau_cache%lmax, &
2194 lmin=tau_cache%lmin, maxso=tau_cache%maxso, &
2195 npgf=tau_cache%npgf, nset=tau_cache%nset, &
2196 zet=tau_cache%zet)
2197 CALL get_paw_basis_info(basis_1c, o2nindex=tau_cache%o2nindex, &
2198 n2oindex=tau_cache%n2oindex, &
2199 nsatbas=tau_cache%nsatbas)
2200
2201 tau_cache%nr = grid_atom%nr
2202 tau_cache%na = grid_atom%ng_sphere
2203 slm => harmonics%slm
2204 dslm_dxyz => harmonics%dslm_dxyz
2205
2206 ALLOCATE (tau_cache%grad(tau_cache%na*tau_cache%nr, tau_cache%nsatbas, 3))
2207 ALLOCATE (a1(tau_cache%na), a2(tau_cache%na), gexp(tau_cache%nr), &
2208 r1(tau_cache%nr), r2(tau_cache%nr))
2209 tau_cache%grad = 0.0_dp
2210
2211 DO iset = 1, tau_cache%nset
2212 DO ipgf = 1, tau_cache%npgf(iset)
2213 starti = (iset - 1)*tau_cache%maxso + &
2214 (ipgf - 1)*nsoset(tau_cache%lmax(iset))
2215 gexp(1:tau_cache%nr) = exp(-tau_cache%zet(ipgf, iset)* &
2216 grid_atom%rad2(1:tau_cache%nr))
2217 DO iso = nsoset(tau_cache%lmin(iset) - 1) + 1, nsoset(tau_cache%lmax(iset))
2218 ip = tau_cache%o2nindex(starti + iso)
2219 IF (ip == 0) cycle
2220 l = indso(1, iso)
2221
2222 r1(1:tau_cache%nr) = grid_atom%rad(1:tau_cache%nr)**(l - 1)*gexp(1:tau_cache%nr)
2223 r2(1:tau_cache%nr) = -2.0_dp*tau_cache%zet(ipgf, iset)* &
2224 grid_atom%rad2(1:tau_cache%nr)*r1(1:tau_cache%nr)
2225
2226 DO dir = 1, 3
2227 a1(1:tau_cache%na) = dslm_dxyz(dir, 1:tau_cache%na, iso)
2228 a2(1:tau_cache%na) = harmonics%a(dir, 1:tau_cache%na)*slm(1:tau_cache%na, iso)
2229 DO ir = 1, tau_cache%nr
2230 DO ia = 1, tau_cache%na
2231 igrid = ia + (ir - 1)*tau_cache%na
2232 tau_cache%grad(igrid, ip, dir) = r1(ir)*a1(ia) + r2(ir)*a2(ia)
2233 END DO
2234 END DO
2235 END DO
2236 END DO
2237 END DO
2238 END DO
2239
2240 DEALLOCATE (a1, a2, gexp, r1, r2)
2241
2242 END SUBROUTINE create_tau_basis_cache
2243
2244! **************************************************************************************************
2245!> \brief Release precomputed GAPW meta-GGA tau factors
2246!> \param tau_cache precomputed compact one-center gradient basis
2247! **************************************************************************************************
2248 SUBROUTINE release_tau_basis_cache(tau_cache)
2249
2250 TYPE(tau_basis_cache_type), INTENT(INOUT) :: tau_cache
2251
2252 IF (ALLOCATED(tau_cache%grad)) DEALLOCATE (tau_cache%grad)
2253 IF (ASSOCIATED(tau_cache%n2oindex)) DEALLOCATE (tau_cache%n2oindex)
2254 IF (ASSOCIATED(tau_cache%o2nindex)) DEALLOCATE (tau_cache%o2nindex)
2255 NULLIFY (tau_cache%lmax, tau_cache%lmin, tau_cache%n2oindex, tau_cache%npgf, &
2256 tau_cache%zet, tau_cache%o2nindex)
2257 tau_cache%maxso = 0
2258 tau_cache%na = 0
2259 tau_cache%nr = 0
2260 tau_cache%nsatbas = 0
2261 tau_cache%nset = 0
2262
2263 END SUBROUTINE release_tau_basis_cache
2264
2265! **************************************************************************************************
2266!> \brief Computes tau hard and soft on the atomic grids for meta-GGA calculations
2267!> \param tau_h the hard part of tau
2268!> \param tau_s the soft part of tau
2269!> \param rho_atom atom-centered density matrices
2270!> \param tau_cache precomputed compact one-center gradient basis
2271!> \param nspins number of spin channels
2272!> \note This is a rewrite to correct a meta-GGA GAPW bug. This is more brute force than the original,
2273!> which was done along in qs_rho_atom_methods.F, but makes sure that no corner is cut in
2274!> terms of accuracy (A. Bussy)
2275! **************************************************************************************************
2276 SUBROUTINE calc_tau_atom(tau_h, tau_s, rho_atom, tau_cache, nspins)
2277
2278 REAL(dp), DIMENSION(:, :, :), INTENT(INOUT) :: tau_h, tau_s
2279 TYPE(rho_atom_type), POINTER :: rho_atom
2280 TYPE(tau_basis_cache_type), INTENT(IN) :: tau_cache
2281 INTEGER, INTENT(IN) :: nspins
2282
2283 CHARACTER(len=*), PARAMETER :: routinen = 'calc_tau_atom'
2284
2285 INTEGER :: dir, handle, ia, ibas, igrid, ir, ispin, &
2286 na, nbas, ngrid, nr
2287 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: work
2288
2289 CALL timeset(routinen, handle)
2290
2291 cpassert(ALLOCATED(tau_cache%grad))
2292
2293 !zeroing tau, assuming it is already allocated
2294 tau_h = 0.0_dp
2295 tau_s = 0.0_dp
2296
2297 nr = tau_cache%nr
2298 na = tau_cache%na
2299 nbas = tau_cache%nsatbas
2300 ngrid = na*nr
2301 ALLOCATE (work(ngrid, nbas))
2302
2303 DO ispin = 1, nspins
2304 DO dir = 1, 3
2305 CALL dgemm('N', 'T', ngrid, nbas, nbas, 0.5_dp, tau_cache%grad(:, :, dir), &
2306 ngrid, rho_atom%cpc_h(ispin)%r_coef, nbas, 0.0_dp, work, ngrid)
2307 DO ibas = 1, nbas
2308 DO ir = 1, nr
2309 DO ia = 1, na
2310 igrid = ia + (ir - 1)*na
2311 tau_h(ia, ir, ispin) = tau_h(ia, ir, ispin) + &
2312 tau_cache%grad(igrid, ibas, dir)*work(igrid, ibas)
2313 END DO
2314 END DO
2315 END DO
2316
2317 CALL dgemm('N', 'T', ngrid, nbas, nbas, 0.5_dp, tau_cache%grad(:, :, dir), &
2318 ngrid, rho_atom%cpc_s(ispin)%r_coef, nbas, 0.0_dp, work, ngrid)
2319 DO ibas = 1, nbas
2320 DO ir = 1, nr
2321 DO ia = 1, na
2322 igrid = ia + (ir - 1)*na
2323 tau_s(ia, ir, ispin) = tau_s(ia, ir, ispin) + &
2324 tau_cache%grad(igrid, ibas, dir)*work(igrid, ibas)
2325 END DO
2326 END DO
2327 END DO
2328 END DO
2329 END DO
2330
2331 DEALLOCATE (work)
2332
2333 CALL timestop(handle)
2334
2335 END SUBROUTINE calc_tau_atom
2336
2337! **************************************************************************************************
2338!> \brief ...
2339!> \param grid_atom ...
2340!> \param nspins ...
2341!> \param grad_func ...
2342!> \param ir ...
2343!> \param rho_nlcc ...
2344!> \param rho_h ...
2345!> \param rho_s ...
2346!> \param drho_nlcc ...
2347!> \param drho_h ...
2348!> \param drho_s ...
2349! **************************************************************************************************
2350 SUBROUTINE calc_rho_nlcc(grid_atom, nspins, grad_func, &
2351 ir, rho_nlcc, rho_h, rho_s, drho_nlcc, drho_h, drho_s)
2352
2353 TYPE(grid_atom_type), POINTER :: grid_atom
2354 INTEGER, INTENT(IN) :: nspins
2355 LOGICAL, INTENT(IN) :: grad_func
2356 INTEGER, INTENT(IN) :: ir
2357 REAL(kind=dp), DIMENSION(:) :: rho_nlcc
2358 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: rho_h, rho_s
2359 REAL(kind=dp), DIMENSION(:) :: drho_nlcc
2360 REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: drho_h, drho_s
2361
2362 INTEGER :: ia, ispin, na
2363 REAL(kind=dp) :: drho, dx, dy, dz, rad, rho, urad, xsp
2364
2365 cpassert(ASSOCIATED(rho_h))
2366 cpassert(ASSOCIATED(rho_s))
2367 IF (grad_func) THEN
2368 cpassert(ASSOCIATED(drho_h))
2369 cpassert(ASSOCIATED(drho_s))
2370 END IF
2371
2372 na = grid_atom%ng_sphere
2373 rad = grid_atom%rad(ir)
2374 urad = grid_atom%oorad2l(ir, 1)
2375
2376 xsp = real(nspins, kind=dp)
2377 rho = rho_nlcc(ir)/xsp
2378 DO ispin = 1, nspins
2379 rho_h(1:na, ir, ispin) = rho_h(1:na, ir, ispin) + rho
2380 rho_s(1:na, ir, ispin) = rho_s(1:na, ir, ispin) + rho
2381 END DO ! ispin
2382
2383 IF (grad_func) THEN
2384 drho = drho_nlcc(ir)/xsp
2385 DO ispin = 1, nspins
2386 DO ia = 1, na
2387 IF (grid_atom%azi(ia) == 0.0_dp) THEN
2388 dx = 0.0_dp
2389 dy = 0.0_dp
2390 ELSE
2391 dx = grid_atom%sin_pol(ia)*grid_atom%sin_azi(ia)
2392 dy = grid_atom%sin_pol(ia)*grid_atom%cos_azi(ia)
2393 END IF
2394 dz = grid_atom%cos_pol(ia)
2395 ! components of the gradient of rho1 hard
2396 drho_h(1, ia, ir, ispin) = drho_h(1, ia, ir, ispin) + drho*dx
2397 drho_h(2, ia, ir, ispin) = drho_h(2, ia, ir, ispin) + drho*dy
2398 drho_h(3, ia, ir, ispin) = drho_h(3, ia, ir, ispin) + drho*dz
2399 ! components of the gradient of rho1 soft
2400 drho_s(1, ia, ir, ispin) = drho_s(1, ia, ir, ispin) + drho*dx
2401 drho_s(2, ia, ir, ispin) = drho_s(2, ia, ir, ispin) + drho*dy
2402 drho_s(3, ia, ir, ispin) = drho_s(3, ia, ir, ispin) + drho*dz
2403 ! norm of gradient
2404 drho_h(4, ia, ir, ispin) = sqrt( &
2405 drho_h(1, ia, ir, ispin)*drho_h(1, ia, ir, ispin) + &
2406 drho_h(2, ia, ir, ispin)*drho_h(2, ia, ir, ispin) + &
2407 drho_h(3, ia, ir, ispin)*drho_h(3, ia, ir, ispin))
2408
2409 drho_s(4, ia, ir, ispin) = sqrt( &
2410 drho_s(1, ia, ir, ispin)*drho_s(1, ia, ir, ispin) + &
2411 drho_s(2, ia, ir, ispin)*drho_s(2, ia, ir, ispin) + &
2412 drho_s(3, ia, ir, ispin)*drho_s(3, ia, ir, ispin))
2413 END DO ! ia
2414 END DO ! ispin
2415 END IF
2416
2417 END SUBROUTINE calc_rho_nlcc
2418
2419! **************************************************************************************************
2420!> \brief ...
2421!> \param vxc_h ...
2422!> \param vxc_s ...
2423!> \param int_hh ...
2424!> \param int_ss ...
2425!> \param grid_atom ...
2426!> \param basis_1c ...
2427!> \param harmonics ...
2428!> \param nspins ...
2429! **************************************************************************************************
2430 SUBROUTINE gavxcgb_nogc(vxc_h, vxc_s, int_hh, int_ss, grid_atom, basis_1c, harmonics, nspins)
2431
2432 REAL(dp), DIMENSION(:, :, :), POINTER :: vxc_h, vxc_s
2433 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_hh, int_ss
2434 TYPE(grid_atom_type), POINTER :: grid_atom
2435 TYPE(gto_basis_set_type), POINTER :: basis_1c
2436 TYPE(harmonics_atom_type), POINTER :: harmonics
2437 INTEGER, INTENT(IN) :: nspins
2438
2439 CHARACTER(len=*), PARAMETER :: routinen = 'gaVxcgb_noGC'
2440
2441 INTEGER :: handle, ia, ic, icg, ipgf1, ipgf2, ir, iset1, iset2, iso, iso1, iso2, ispin, l, &
2442 ld, lmax12, lmax_expansion, lmin12, m1, m2, max_iso_not0, max_iso_not0_local, max_s_harm, &
2443 maxl, maxso, n1, n2, na, ngau1, ngau2, nngau1, nr, nset, size1
2444 INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list
2445 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list
2446 INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
2447 REAL(dp), ALLOCATABLE, DIMENSION(:) :: g1, g2
2448 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: gg, gvg_h, gvg_s, matso_h, matso_s, vx
2449 REAL(dp), DIMENSION(:, :), POINTER :: zet
2450 REAL(dp), DIMENSION(:, :, :), POINTER :: my_cg
2451
2452 CALL timeset(routinen, handle)
2453
2454 NULLIFY (lmin, lmax, npgf, zet, my_cg)
2455
2456 CALL get_gto_basis_set(gto_basis_set=basis_1c, lmax=lmax, lmin=lmin, &
2457 maxso=maxso, maxl=maxl, npgf=npgf, &
2458 nset=nset, zet=zet)
2459
2460 nr = grid_atom%nr
2461 na = grid_atom%ng_sphere
2462 my_cg => harmonics%my_CG
2463 max_iso_not0 = harmonics%max_iso_not0
2464 lmax_expansion = indso(1, max_iso_not0)
2465 max_s_harm = harmonics%max_s_harm
2466
2467 ALLOCATE (g1(nr), g2(nr), gg(nr, 0:2*maxl))
2468 ALLOCATE (gvg_h(na, 0:2*maxl), gvg_s(na, 0:2*maxl))
2469 ALLOCATE (matso_h(nsoset(maxl), nsoset(maxl)), &
2470 matso_s(nsoset(maxl), nsoset(maxl)))
2471 ALLOCATE (vx(na, nr))
2472 ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm))
2473
2474 g1 = 0.0_dp
2475 g2 = 0.0_dp
2476 m1 = 0
2477 DO iset1 = 1, nset
2478 n1 = nsoset(lmax(iset1))
2479 m2 = 0
2480 DO iset2 = 1, nset
2481 CALL get_none0_cg_list(my_cg, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
2482 max_s_harm, lmax_expansion, cg_list, cg_n_list, max_iso_not0_local)
2483 cpassert(max_iso_not0_local <= max_iso_not0)
2484
2485 n2 = nsoset(lmax(iset2))
2486 DO ipgf1 = 1, npgf(iset1)
2487 ngau1 = n1*(ipgf1 - 1) + m1
2488 size1 = nsoset(lmax(iset1)) - nsoset(lmin(iset1) - 1)
2489 nngau1 = nsoset(lmin(iset1) - 1) + ngau1
2490
2491 g1(1:nr) = exp(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))
2492 DO ipgf2 = 1, npgf(iset2)
2493 ngau2 = n2*(ipgf2 - 1) + m2
2494
2495 g2(1:nr) = exp(-zet(ipgf2, iset2)*grid_atom%rad2(1:nr))
2496 lmin12 = lmin(iset1) + lmin(iset2)
2497 lmax12 = lmax(iset1) + lmax(iset2)
2498
2499 ! reduce expansion local densities
2500 IF (lmin12 <= lmax_expansion) THEN
2501
2502 gg = 0.0_dp
2503 IF (lmin12 == 0) THEN
2504 gg(1:nr, lmin12) = g1(1:nr)*g2(1:nr)
2505 ELSE
2506 gg(1:nr, lmin12) = grid_atom%rad2l(1:nr, lmin12)*g1(1:nr)*g2(1:nr)
2507 END IF
2508
2509 ! limit the expansion of the local densities to a max L
2510 IF (lmax12 > lmax_expansion) lmax12 = lmax_expansion
2511
2512 DO l = lmin12 + 1, lmax12
2513 gg(1:nr, l) = grid_atom%rad(1:nr)*gg(:, l - 1)
2514 END DO
2515
2516 DO ispin = 1, nspins
2517 ld = lmax12 + 1
2518 DO ir = 1, nr
2519 vx(1:na, ir) = vxc_h(1:na, ir, ispin)
2520 END DO
2521 CALL dgemm('N', 'N', na, ld, nr, 1.0_dp, vx(1:na, 1:nr), na, &
2522 gg(1:nr, 0:lmax12), nr, 0.0_dp, gvg_h(1:na, 0:lmax12), na)
2523 DO ir = 1, nr
2524 vx(1:na, ir) = vxc_s(1:na, ir, ispin)
2525 END DO
2526 CALL dgemm('N', 'N', na, ld, nr, 1.0_dp, vx(1:na, 1:nr), na, &
2527 gg(1:nr, 0:lmax12), nr, 0.0_dp, gvg_s(1:na, 0:lmax12), na)
2528
2529 matso_h = 0.0_dp
2530 matso_s = 0.0_dp
2531 DO iso = 1, max_iso_not0_local
2532 DO icg = 1, cg_n_list(iso)
2533 iso1 = cg_list(1, icg, iso)
2534 iso2 = cg_list(2, icg, iso)
2535 l = indso(1, iso1) + indso(1, iso2)
2536
2537 cpassert(l <= lmax_expansion)
2538 DO ia = 1, na
2539 matso_h(iso1, iso2) = matso_h(iso1, iso2) + &
2540 gvg_h(ia, l)* &
2541 my_cg(iso1, iso2, iso)* &
2542 harmonics%slm(ia, iso)
2543 matso_s(iso1, iso2) = matso_s(iso1, iso2) + &
2544 gvg_s(ia, l)* &
2545 my_cg(iso1, iso2, iso)* &
2546 harmonics%slm(ia, iso)
2547 END DO
2548 END DO
2549 END DO
2550
2551 ! Write in the global matrix
2552 DO ic = nsoset(lmin(iset2) - 1) + 1, nsoset(lmax(iset2))
2553 iso1 = nsoset(lmin(iset1) - 1) + 1
2554 iso2 = ngau2 + ic
2555 CALL daxpy(size1, 1.0_dp, matso_h(iso1, ic), 1, &
2556 int_hh(ispin)%r_coef(nngau1 + 1, iso2), 1)
2557 CALL daxpy(size1, 1.0_dp, matso_s(iso1, ic), 1, &
2558 int_ss(ispin)%r_coef(nngau1 + 1, iso2), 1)
2559 END DO
2560
2561 END DO ! ispin
2562
2563 END IF ! lmax_expansion
2564
2565 END DO ! ipfg2
2566 END DO ! ipfg1
2567 m2 = m2 + maxso
2568 END DO ! iset2
2569 m1 = m1 + maxso
2570 END DO ! iset1
2571
2572 DEALLOCATE (g1, g2, gg, matso_h, matso_s, gvg_s, gvg_h, vx)
2573
2574 DEALLOCATE (cg_list, cg_n_list)
2575
2576 CALL timestop(handle)
2577
2578 END SUBROUTINE gavxcgb_nogc
2579
2580! **************************************************************************************************
2581!> \brief ...
2582!> \param vxc_h ...
2583!> \param vxc_s ...
2584!> \param vxg_h ...
2585!> \param vxg_s ...
2586!> \param int_hh ...
2587!> \param int_ss ...
2588!> \param grid_atom ...
2589!> \param basis_1c ...
2590!> \param harmonics ...
2591!> \param nspins ...
2592! **************************************************************************************************
2593 SUBROUTINE gavxcgb_gc(vxc_h, vxc_s, vxg_h, vxg_s, int_hh, int_ss, &
2594 grid_atom, basis_1c, harmonics, nspins)
2595
2596 REAL(dp), DIMENSION(:, :, :), POINTER :: vxc_h, vxc_s
2597 REAL(dp), DIMENSION(:, :, :, :), POINTER :: vxg_h, vxg_s
2598 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_hh, int_ss
2599 TYPE(grid_atom_type), POINTER :: grid_atom
2600 TYPE(gto_basis_set_type), POINTER :: basis_1c
2601 TYPE(harmonics_atom_type), POINTER :: harmonics
2602 INTEGER, INTENT(IN) :: nspins
2603
2604 CHARACTER(len=*), PARAMETER :: routinen = 'gaVxcgb_GC'
2605
2606 INTEGER :: dmax_iso_not0_local, handle, ia, ic, icg, ipgf1, ipgf2, ir, iset1, iset2, iso, &
2607 iso1, iso2, ispin, l, lmax12, lmax_expansion, lmin12, m1, m2, max_iso_not0, &
2608 max_iso_not0_local, max_s_harm, maxl, maxso, n1, n2, na, ngau1, ngau2, nngau1, nr, nset, &
2609 size1
2610 INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list, dcg_n_list
2611 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list, dcg_list
2612 INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
2613 REAL(dp) :: urad
2614 REAL(dp), ALLOCATABLE, DIMENSION(:) :: g1, g2
2615 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: dgg, gg, gvxcg_h, gvxcg_s, matso_h, &
2616 matso_s
2617 REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: gvxgg_h, gvxgg_s
2618 REAL(dp), DIMENSION(:, :), POINTER :: zet
2619 REAL(dp), DIMENSION(:, :, :), POINTER :: my_cg
2620 REAL(dp), DIMENSION(:, :, :, :), POINTER :: my_cg_dxyz
2621
2622 CALL timeset(routinen, handle)
2623
2624 NULLIFY (lmin, lmax, npgf, zet, my_cg, my_cg_dxyz)
2625
2626 CALL get_gto_basis_set(gto_basis_set=basis_1c, lmax=lmax, lmin=lmin, &
2627 maxso=maxso, maxl=maxl, npgf=npgf, &
2628 nset=nset, zet=zet)
2629
2630 nr = grid_atom%nr
2631 na = grid_atom%ng_sphere
2632 my_cg => harmonics%my_CG
2633 my_cg_dxyz => harmonics%my_CG_dxyz
2634 max_iso_not0 = harmonics%max_iso_not0
2635 lmax_expansion = indso(1, max_iso_not0)
2636 max_s_harm = harmonics%max_s_harm
2637
2638 ALLOCATE (g1(nr), g2(nr), gg(nr, 0:2*maxl), dgg(nr, 0:2*maxl))
2639 ALLOCATE (gvxcg_h(na, 0:2*maxl), gvxcg_s(na, 0:2*maxl))
2640 ALLOCATE (gvxgg_h(3, na, 0:2*maxl), gvxgg_s(3, na, 0:2*maxl))
2641 ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm), &
2642 dcg_list(2, nsoset(maxl)**2, max_s_harm), dcg_n_list(max_s_harm))
2643
2644 ALLOCATE (matso_h(nsoset(maxl), nsoset(maxl)), &
2645 matso_s(nsoset(maxl), nsoset(maxl)))
2646
2647 DO ispin = 1, nspins
2648
2649 g1 = 0.0_dp
2650 g2 = 0.0_dp
2651 m1 = 0
2652 DO iset1 = 1, nset
2653 n1 = nsoset(lmax(iset1))
2654 m2 = 0
2655 DO iset2 = 1, nset
2656 CALL get_none0_cg_list(my_cg, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
2657 max_s_harm, lmax_expansion, cg_list, cg_n_list, max_iso_not0_local)
2658 cpassert(max_iso_not0_local <= max_iso_not0)
2659 CALL get_none0_cg_list(my_cg_dxyz, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
2660 max_s_harm, lmax_expansion, dcg_list, dcg_n_list, dmax_iso_not0_local)
2661
2662 n2 = nsoset(lmax(iset2))
2663 DO ipgf1 = 1, npgf(iset1)
2664 ngau1 = n1*(ipgf1 - 1) + m1
2665 size1 = nsoset(lmax(iset1)) - nsoset(lmin(iset1) - 1)
2666 nngau1 = nsoset(lmin(iset1) - 1) + ngau1
2667
2668 g1(1:nr) = exp(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))
2669 DO ipgf2 = 1, npgf(iset2)
2670 ngau2 = n2*(ipgf2 - 1) + m2
2671
2672 g2(1:nr) = exp(-zet(ipgf2, iset2)*grid_atom%rad2(1:nr))
2673 lmin12 = lmin(iset1) + lmin(iset2)
2674 lmax12 = lmax(iset1) + lmax(iset2)
2675
2676 !test reduce expansion local densities
2677 IF (lmin12 <= lmax_expansion) THEN
2678
2679 gg = 0.0_dp
2680 dgg = 0.0_dp
2681
2682 IF (lmin12 == 0) THEN
2683 gg(1:nr, lmin12) = g1(1:nr)*g2(1:nr)
2684 ELSE
2685 gg(1:nr, lmin12) = grid_atom%rad2l(1:nr, lmin12)*g1(1:nr)*g2(1:nr)
2686 END IF
2687
2688 !test reduce expansion local densities
2689 IF (lmax12 > lmax_expansion) lmax12 = lmax_expansion
2690
2691 DO l = lmin12 + 1, lmax12
2692 gg(1:nr, l) = grid_atom%rad(1:nr)*gg(:, l - 1)
2693 dgg(1:nr, l - 1) = dgg(1:nr, l - 1) - 2.0_dp*(zet(ipgf1, iset1) + &
2694 zet(ipgf2, iset2))*gg(1:nr, l)
2695 END DO
2696 dgg(1:nr, lmax12) = dgg(1:nr, lmax12) - 2.0_dp*(zet(ipgf1, iset1) + &
2697 zet(ipgf2, iset2))*grid_atom%rad(1:nr)* &
2698 gg(1:nr, lmax12)
2699
2700 gvxcg_h = 0.0_dp
2701 gvxcg_s = 0.0_dp
2702 gvxgg_h = 0.0_dp
2703 gvxgg_s = 0.0_dp
2704
2705 ! Cross Term
2706 DO l = lmin12, lmax12
2707 DO ia = 1, na
2708 DO ir = 1, nr
2709 gvxcg_h(ia, l) = gvxcg_h(ia, l) + &
2710 gg(ir, l)*vxc_h(ia, ir, ispin) + &
2711 dgg(ir, l)* &
2712 (vxg_h(1, ia, ir, ispin)*harmonics%a(1, ia) + &
2713 vxg_h(2, ia, ir, ispin)*harmonics%a(2, ia) + &
2714 vxg_h(3, ia, ir, ispin)*harmonics%a(3, ia))
2715
2716 gvxcg_s(ia, l) = gvxcg_s(ia, l) + &
2717 gg(ir, l)*vxc_s(ia, ir, ispin) + &
2718 dgg(ir, l)* &
2719 (vxg_s(1, ia, ir, ispin)*harmonics%a(1, ia) + &
2720 vxg_s(2, ia, ir, ispin)*harmonics%a(2, ia) + &
2721 vxg_s(3, ia, ir, ispin)*harmonics%a(3, ia))
2722
2723 urad = grid_atom%oorad2l(ir, 1)
2724
2725 gvxgg_h(1, ia, l) = gvxgg_h(1, ia, l) + &
2726 vxg_h(1, ia, ir, ispin)* &
2727 gg(ir, l)*urad
2728
2729 gvxgg_h(2, ia, l) = gvxgg_h(2, ia, l) + &
2730 vxg_h(2, ia, ir, ispin)* &
2731 gg(ir, l)*urad
2732
2733 gvxgg_h(3, ia, l) = gvxgg_h(3, ia, l) + &
2734 vxg_h(3, ia, ir, ispin)* &
2735 gg(ir, l)*urad
2736
2737 gvxgg_s(1, ia, l) = gvxgg_s(1, ia, l) + &
2738 vxg_s(1, ia, ir, ispin)* &
2739 gg(ir, l)*urad
2740
2741 gvxgg_s(2, ia, l) = gvxgg_s(2, ia, l) + &
2742 vxg_s(2, ia, ir, ispin)* &
2743 gg(ir, l)*urad
2744
2745 gvxgg_s(3, ia, l) = gvxgg_s(3, ia, l) + &
2746 vxg_s(3, ia, ir, ispin)* &
2747 gg(ir, l)*urad
2748
2749 END DO ! ir
2750 END DO ! ia
2751 END DO ! l
2752
2753 matso_h = 0.0_dp
2754 matso_s = 0.0_dp
2755 DO iso = 1, max_iso_not0_local
2756 DO icg = 1, cg_n_list(iso)
2757 iso1 = cg_list(1, icg, iso)
2758 iso2 = cg_list(2, icg, iso)
2759
2760 l = indso(1, iso1) + indso(1, iso2)
2761
2762 !test reduce expansion local densities
2763 cpassert(l <= lmax_expansion)
2764 DO ia = 1, na
2765 matso_h(iso1, iso2) = matso_h(iso1, iso2) + &
2766 gvxcg_h(ia, l)* &
2767 harmonics%slm(ia, iso)* &
2768 my_cg(iso1, iso2, iso)
2769 matso_s(iso1, iso2) = matso_s(iso1, iso2) + &
2770 gvxcg_s(ia, l)* &
2771 harmonics%slm(ia, iso)* &
2772 my_cg(iso1, iso2, iso)
2773 END DO ! ia
2774
2775 !test reduce expansion local densities
2776
2777 END DO
2778
2779 END DO ! iso
2780
2781 DO iso = 1, dmax_iso_not0_local
2782 DO icg = 1, dcg_n_list(iso)
2783 iso1 = dcg_list(1, icg, iso)
2784 iso2 = dcg_list(2, icg, iso)
2785
2786 l = indso(1, iso1) + indso(1, iso2)
2787 !test reduce expansion local densities
2788 cpassert(l <= lmax_expansion)
2789 DO ia = 1, na
2790 matso_h(iso1, iso2) = matso_h(iso1, iso2) + &
2791 (gvxgg_h(1, ia, l)*my_cg_dxyz(1, iso1, iso2, iso) + &
2792 gvxgg_h(2, ia, l)*my_cg_dxyz(2, iso1, iso2, iso) + &
2793 gvxgg_h(3, ia, l)*my_cg_dxyz(3, iso1, iso2, iso))* &
2794 harmonics%slm(ia, iso)
2795
2796 matso_s(iso1, iso2) = matso_s(iso1, iso2) + &
2797 (gvxgg_s(1, ia, l)*my_cg_dxyz(1, iso1, iso2, iso) + &
2798 gvxgg_s(2, ia, l)*my_cg_dxyz(2, iso1, iso2, iso) + &
2799 gvxgg_s(3, ia, l)*my_cg_dxyz(3, iso1, iso2, iso))* &
2800 harmonics%slm(ia, iso)
2801
2802 END DO ! ia
2803
2804 !test reduce expansion local densities
2805
2806 END DO ! icg
2807 END DO ! iso
2808 !test reduce expansion local densities
2809 END IF ! lmax_expansion
2810
2811 ! Write in the global matrix
2812 DO ic = nsoset(lmin(iset2) - 1) + 1, nsoset(lmax(iset2))
2813 iso1 = nsoset(lmin(iset1) - 1) + 1
2814 iso2 = ngau2 + ic
2815 CALL daxpy(size1, 1.0_dp, matso_h(iso1, ic), 1, &
2816 int_hh(ispin)%r_coef(nngau1 + 1, iso2), 1)
2817 CALL daxpy(size1, 1.0_dp, matso_s(iso1, ic), 1, &
2818 int_ss(ispin)%r_coef(nngau1 + 1, iso2), 1)
2819 END DO
2820
2821 END DO ! ipfg2
2822 END DO ! ipfg1
2823 m2 = m2 + maxso
2824 END DO ! iset2
2825 m1 = m1 + maxso
2826 END DO ! iset1
2827 END DO ! ispin
2828
2829 DEALLOCATE (g1, g2, gg, dgg, matso_h, matso_s, gvxcg_h, gvxcg_s, gvxgg_h, gvxgg_s)
2830 DEALLOCATE (cg_list, cg_n_list, dcg_list, dcg_n_list)
2831
2832 CALL timestop(handle)
2833
2834 END SUBROUTINE gavxcgb_gc
2835
2836! **************************************************************************************************
2837!> \brief Integrates 0.5 * grad_ga .dot. (V_tau * grad_gb) on the atomic grid for meta-GGA
2838!> \param vtau_h the hard tau potential
2839!> \param vtau_s the soft tau potential
2840!> \param int_hh hard one-center matrix contribution
2841!> \param int_ss soft one-center matrix contribution
2842!> \param tau_cache precomputed compact one-center gradient basis
2843!> \param nspins number of spin channels
2844!> \note This is a rewrite to correct meta-GGA GAPW bug. This is more brute force than the original
2845!> but makes sure that no corner is cut in terms of accuracy (A. Bussy)
2846! **************************************************************************************************
2847 SUBROUTINE dgavtaudgb(vtau_h, vtau_s, int_hh, int_ss, &
2848 tau_cache, nspins)
2849
2850 REAL(dp), DIMENSION(:, :, :), POINTER :: vtau_h, vtau_s
2851 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_hh, int_ss
2852 TYPE(tau_basis_cache_type), INTENT(IN) :: tau_cache
2853 INTEGER, INTENT(IN) :: nspins
2854
2855 CHARACTER(len=*), PARAMETER :: routinen = 'dgaVtaudgb'
2856
2857 INTEGER :: dir, handle, ia, ibas, igrid, iold, ir, &
2858 ispin, jbas, jold, max_old_basis, na, &
2859 nbas, ngrid, nr
2860 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: int_h, int_s, weighted_grad
2861
2862 CALL timeset(routinen, handle)
2863
2864 cpassert(ALLOCATED(tau_cache%grad))
2865 cpassert(ASSOCIATED(tau_cache%n2oindex))
2866
2867 nr = tau_cache%nr
2868 na = tau_cache%na
2869 nbas = tau_cache%nsatbas
2870 ngrid = na*nr
2871 max_old_basis = maxval(tau_cache%n2oindex)
2872 ALLOCATE (int_h(nbas, nbas), int_s(nbas, nbas), weighted_grad(ngrid, nbas))
2873
2874 DO ispin = 1, nspins
2875 cpassert(SIZE(int_hh(ispin)%r_coef, 1) >= max_old_basis)
2876 cpassert(SIZE(int_hh(ispin)%r_coef, 2) >= max_old_basis)
2877 cpassert(SIZE(int_ss(ispin)%r_coef, 1) >= max_old_basis)
2878 cpassert(SIZE(int_ss(ispin)%r_coef, 2) >= max_old_basis)
2879 int_h = 0.0_dp
2880 int_s = 0.0_dp
2881 DO dir = 1, 3
2882 DO ibas = 1, nbas
2883 DO ir = 1, nr
2884 DO ia = 1, na
2885 igrid = ia + (ir - 1)*na
2886 weighted_grad(igrid, ibas) = vtau_h(ia, ir, ispin)* &
2887 tau_cache%grad(igrid, ibas, dir)
2888 END DO
2889 END DO
2890 END DO
2891 CALL dgemm('T', 'N', nbas, nbas, ngrid, 0.5_dp, tau_cache%grad(:, :, dir), &
2892 ngrid, weighted_grad, ngrid, 1.0_dp, int_h, nbas)
2893
2894 DO ibas = 1, nbas
2895 DO ir = 1, nr
2896 DO ia = 1, na
2897 igrid = ia + (ir - 1)*na
2898 weighted_grad(igrid, ibas) = vtau_s(ia, ir, ispin)* &
2899 tau_cache%grad(igrid, ibas, dir)
2900 END DO
2901 END DO
2902 END DO
2903 CALL dgemm('T', 'N', nbas, nbas, ngrid, 0.5_dp, tau_cache%grad(:, :, dir), &
2904 ngrid, weighted_grad, ngrid, 1.0_dp, int_s, nbas)
2905 END DO
2906
2907 DO jbas = 1, nbas
2908 jold = tau_cache%n2oindex(jbas)
2909 DO ibas = 1, nbas
2910 iold = tau_cache%n2oindex(ibas)
2911 int_hh(ispin)%r_coef(iold, jold) = int_hh(ispin)%r_coef(iold, jold) + &
2912 int_h(ibas, jbas)
2913 int_ss(ispin)%r_coef(iold, jold) = int_ss(ispin)%r_coef(iold, jold) + &
2914 int_s(ibas, jbas)
2915 END DO
2916 END DO
2917 END DO
2918
2919 DEALLOCATE (int_h, int_s, weighted_grad)
2920
2921 CALL timestop(handle)
2922
2923 END SUBROUTINE dgavtaudgb
2924
2925END MODULE qs_vxc_atom
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:225
memory_utilities::reallocate
Definition memory_utilities.F:27
qs_harmonics_atom::get_none0_cg_list
Definition qs_harmonics_atom.F:54
atom
Definition atom.F:9
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
basis_set_types
Definition basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum, ccon)
...
Definition basis_set_types.F:643
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::xc_none
integer, parameter, public xc_none
Definition input_constants.F:572
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:731
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1047
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
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::nsoset
integer, dimension(:), allocatable, public nsoset
Definition orbital_pointers.F:42
orbital_pointers::indso
integer, dimension(:, :), allocatable, public indso
Definition orbital_pointers.F:43
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
paw_basis_types
Definition paw_basis_types.F:13
paw_basis_types::get_paw_basis_info
subroutine, public get_paw_basis_info(basis_1c, o2nindex, n2oindex, nsatbas)
Return some info on the PAW basis derived from a GTO basis set.
Definition paw_basis_types.F:40
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:530
qs_force_types
Definition qs_force_types.F:13
qs_grid_atom
Definition qs_grid_atom.F:8
qs_harmonics_atom
Definition qs_harmonics_atom.F:8
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::has_nlcc
logical function, public has_nlcc(qs_kind_set)
finds if a given qs run needs to use nlcc
Definition qs_kind_types.F:3973
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:468
qs_linres_types
Type definitiona for linear response calculations.
Definition qs_linres_types.F:12
qs_rho_atom_types
Definition qs_rho_atom_types.F:8
qs_rho_atom_types::get_rho_atom
subroutine, public get_rho_atom(rho_atom, cpc_h, cpc_s, rho_rad_h, rho_rad_s, drho_rad_h, drho_rad_s, vrho_rad_h, vrho_rad_s, rho_rad_h_d, rho_rad_s_d, ga_vlocal_gb_h, ga_vlocal_gb_s, int_scr_h, int_scr_s)
...
Definition qs_rho_atom_types.F:195
qs_vxc_atom
routines that build the integrals of the Vxc potential calculated for the atomic density in the basis...
Definition qs_vxc_atom.F:12
qs_vxc_atom::calculate_xc_2nd_deriv_atom
subroutine, public calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, xc_section, para_env, do_tddfpt2, do_triplet, do_sf, kind_set_external)
...
Definition qs_vxc_atom.F:890
qs_vxc_atom::calculate_gfxc_atom
subroutine, public calculate_gfxc_atom(qs_env, rho0_atom_set, rho1_atom_set, rho2_atom_set, kind_set, xc_section, is_triplet, accuracy)
...
Definition qs_vxc_atom.F:1231
qs_vxc_atom::calculate_vxc_atom
subroutine, public calculate_vxc_atom(qs_env, energy_only, exc1, adiabatic_rescale_factor, kind_set_external, rho_atom_set_external, xc_section_external, calculate_forces)
...
Definition qs_vxc_atom.F:107
qs_vxc_atom::calc_rho_angular
subroutine, public calc_rho_angular(grid_atom, harmonics, nspins, grad_func, ir, r_h, r_s, rho_h, rho_s, dr_h, dr_s, r_h_d, r_s_d, drho_h, drho_s)
...
Definition qs_vxc_atom.F:2066
qs_vxc_atom::calculate_vxc_atom_epr
subroutine, public calculate_vxc_atom_epr(qs_env, exc1, gradient_atom_set)
...
Definition qs_vxc_atom.F:543
qs_vxc_atom::gavxcgb_nogc
subroutine, public gavxcgb_nogc(vxc_h, vxc_s, int_hh, int_ss, grid_atom, basis_1c, harmonics, nspins)
...
Definition qs_vxc_atom.F:2431
qs_vxc_atom::gfxc_atom_diff
subroutine, public gfxc_atom_diff(qs_env, rho0_atom_set, rho1_atom_set, rho2_atom_set, kind_set, xc_section, is_triplet, accuracy, epsrho)
...
Definition qs_vxc_atom.F:1657
skala_gpw_functional
Experimental CP2K-native GPW real-space-grid path for SKALA TorchScript models.
Definition skala_gpw_functional.F:11
skala_gpw_functional::skala_gapw_density_partition_soft_only
integer, parameter, public skala_gapw_density_partition_soft_only
Definition skala_gpw_functional.F:72
skala_gpw_functional::skala_gapw_density_partition_none
integer, parameter, public skala_gapw_density_partition_none
Definition skala_gpw_functional.F:72
skala_gpw_functional::xc_section_uses_gauxc_model
logical function, public xc_section_uses_gauxc_model(xc_section)
Return true if the GAUXC subsection requests a Skala-style model.
Definition skala_gpw_functional.F:117
skala_gpw_functional::skala_gapw_density_partition_hard_minus_soft
integer, parameter, public skala_gapw_density_partition_hard_minus_soft
Definition skala_gpw_functional.F:72
skala_gpw_functional::native_skala_gapw_density_partition
integer function, public native_skala_gapw_density_partition(xc_section)
Return the hard/soft GAPW one-center density partition for native SKALA.
Definition skala_gpw_functional.F:140
skala_gpw_functional::skala_gapw_atom_vxc_of_r
subroutine, public skala_gapw_atom_vxc_of_r(xc_section, grid_atom, group, atom_coord, rho, drho, tau, weights, lsd, nspins, na, nr, exc, vxc, vxg, vtau, energy_only, atom_force, atom_virial)
Evaluate SKALA on a GAPW one-center atomic grid.
Definition skala_gpw_functional.F:649
skala_gpw_functional::skala_gapw_density_partition_hard_only
integer, parameter, public skala_gapw_density_partition_hard_only
Definition skala_gpw_functional.F:72
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
virial_types
Definition virial_types.F:13
xc_atom
Definition xc_atom.F:9
xc_atom::vxc_of_r_epr
subroutine, public vxc_of_r_epr(xc_fun_section, rho_set, deriv_set, needs, w, lsd, na, nr, exc, vxc, vxg, vtau)
Specific EPR version of vxc_of_r_new.
Definition xc_atom.F:297
xc_atom::vxc_of_r_new
subroutine, public vxc_of_r_new(xc_fun_section, rho_set, deriv_set, deriv_order, needs, w, lsd, na, nr, exc, vxc, vxg, vtau, energy_only, adiabatic_rescale_factor)
...
Definition xc_atom.F:64
xc_atom::xc_rho_set_atom_update
subroutine, public xc_rho_set_atom_update(rho_set, needs, nspins, bo)
...
Definition xc_atom.F:507
xc_atom::xc_2nd_deriv_of_r
subroutine, public xc_2nd_deriv_of_r(rho_set, rho1_set, xc_section, deriv_set, w, vxc, vxg, vtau, do_triplet, do_sf)
...
Definition xc_atom.F:403
xc_atom::fill_rho_set
subroutine, public fill_rho_set(rho_set, lsd, nspins, needs, rho, drho, tau, na, ir)
...
Definition xc_atom.F:667
xc_derivative_set_types
represent a group ofunctional derivatives
Definition xc_derivative_set_types.F:14
xc_derivative_set_types::xc_dset_zero_all
subroutine, public xc_dset_zero_all(deriv_set)
...
Definition xc_derivative_set_types.F:165
xc_derivative_set_types::xc_dset_release
subroutine, public xc_dset_release(derivative_set)
releases a derivative set
Definition xc_derivative_set_types.F:142
xc_derivative_set_types::xc_dset_create
subroutine, public xc_dset_create(derivative_set, pw_pool, local_bounds)
creates a derivative set object
Definition xc_derivative_set_types.F:111
xc_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:610
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
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cp_control_types::dft_control_type
Definition cp_control_types.F:751
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:743
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:220
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_grid_atom::grid_atom_type
Definition qs_grid_atom.F:44
qs_harmonics_atom::harmonics_atom_type
Definition qs_harmonics_atom.F:32
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:177
qs_linres_types::nablavks_atom_type
Definition qs_linres_types.F:213
qs_rho_atom_types::rho_atom_coeff
Definition qs_rho_atom_types.F:19
qs_rho_atom_types::rho_atom_type
Definition qs_rho_atom_types.F:23
virial_types::virial_type
Definition virial_types.F:26
xc_derivative_set_types::xc_derivative_set_type
A derivative set contains the different derivatives of a xc-functional in form of a linked list.
Definition xc_derivative_set_types.F:48
xc_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