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