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