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