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