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