Line data Source code
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 : ! **************************************************************************************************
12 : MODULE 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 :
81 : CONTAINS
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 29294 : 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 29294 : 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 58588 : REAL(dp), DIMENSION(:, :), POINTER :: rho_nlcc, weight_h, weight_s
121 58588 : REAL(dp), DIMENSION(:, :, :), POINTER :: rho_h, rho_s, tau_h, tau_s, vtau_h, &
122 29294 : vtau_s, vxc_h, vxc_s
123 58588 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: drho_h, drho_s, vxg_h, vxg_s
124 29294 : 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 29294 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: my_kind_set
131 29294 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s, int_hh, int_ss, r_h, r_s
132 29294 : TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r_h_d, r_s_d
133 29294 : 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 29294 : 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 29294 : CALL timeset(routineN, handle)
144 :
145 29294 : NULLIFY (atom_list)
146 29294 : NULLIFY (my_kind_set)
147 29294 : NULLIFY (atomic_kind_set)
148 29294 : NULLIFY (grid_atom)
149 29294 : NULLIFY (harmonics)
150 29294 : NULLIFY (input)
151 29294 : NULLIFY (para_env)
152 29294 : NULLIFY (rho_atom)
153 29294 : NULLIFY (my_rho_atom_set)
154 29294 : NULLIFY (rho_nlcc)
155 :
156 29294 : IF (PRESENT(adiabatic_rescale_factor)) THEN
157 44 : my_adiabatic_rescale_factor = adiabatic_rescale_factor
158 : ELSE
159 29250 : 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 29294 : rho_atom_set=my_rho_atom_set)
169 :
170 29294 : IF (PRESENT(kind_set_external)) my_kind_set => kind_set_external
171 29294 : IF (PRESENT(rho_atom_set_external)) my_rho_atom_set => rho_atom_set_external
172 :
173 29294 : nlcc = has_nlcc(my_kind_set)
174 29294 : accint = dft_control%qs_control%gapw_control%accurate_xcint
175 :
176 29294 : my_xc_section => section_vals_get_subs_vals(input, "DFT%XC")
177 :
178 29294 : IF (PRESENT(xc_section_external)) my_xc_section => xc_section_external
179 :
180 29294 : 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 29294 : i_val=myfun)
183 :
184 29294 : IF (myfun == xc_none) THEN
185 4018 : exc1 = 0.0_dp
186 15686 : my_rho_atom_set(:)%exc_h = 0.0_dp
187 15686 : my_rho_atom_set(:)%exc_s = 0.0_dp
188 : ELSE
189 : CALL section_vals_val_get(my_xc_section, "DENSITY_CUTOFF", &
190 25276 : r_val=density_cut)
191 : CALL section_vals_val_get(my_xc_section, "GRADIENT_CUTOFF", &
192 25276 : r_val=gradient_cut)
193 : CALL section_vals_val_get(my_xc_section, "TAU_CUTOFF", &
194 25276 : r_val=tau_cut)
195 :
196 25276 : lsd = dft_control%lsd
197 25276 : nspins = dft_control%nspins
198 : needs = xc_functionals_get_needs(xc_fun_section, &
199 : lsd=lsd, &
200 25276 : calc_potential=.TRUE.)
201 :
202 25276 : gradient_f = (needs%drho .OR. needs%drho_spin)
203 25276 : tau_f = (needs%tau .OR. needs%tau_spin)
204 :
205 : ! Initialize energy contribution from the one center XC terms to zero
206 25276 : exc1 = 0.0_dp
207 :
208 : ! Nullify some pointers for work-arrays
209 25276 : NULLIFY (rho_h, drho_h, rho_s, drho_s, weight_h, weight_s)
210 25276 : NULLIFY (vxc_h, vxc_s, vxg_h, vxg_s)
211 25276 : NULLIFY (tau_h, tau_s)
212 25276 : NULLIFY (vtau_h, vtau_s)
213 :
214 : ! Here starts the loop over all the atoms
215 :
216 75324 : DO ikind = 1, SIZE(atomic_kind_set)
217 50048 : 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 50048 : harmonics=harmonics, grid_atom=grid_atom)
220 50048 : CALL get_qs_kind(my_kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
221 :
222 50048 : IF (.NOT. paw_atom) CYCLE
223 :
224 46036 : nr = grid_atom%nr
225 46036 : 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 460360 : bounds(1:2, 1:3) = 1
233 46036 : bounds(2, 1) = na
234 46036 : bounds(2, 2) = nr
235 :
236 : ! set integration weights
237 46036 : IF (accint) THEN
238 13282 : weight_h => grid_atom%weight
239 13282 : alpha = dft_control%qs_control%gapw_control%aw(ikind)
240 13282 : IF (ASSOCIATED(grid_atom%gapw_weight_s)) THEN
241 12814 : IF (grid_atom%gapw_weight_alpha /= alpha) DEALLOCATE (grid_atom%gapw_weight_s)
242 : END IF
243 13282 : IF (.NOT. ASSOCIATED(grid_atom%gapw_weight_s)) THEN
244 1872 : ALLOCATE (grid_atom%gapw_weight_s(na, nr))
245 23868 : DO ir = 1, nr
246 23400 : agr = 1.0_dp - EXP(-alpha*grid_atom%rad2(ir))
247 1193868 : grid_atom%gapw_weight_s(:, ir) = grid_atom%weight(:, ir)*agr
248 : END DO
249 468 : grid_atom%gapw_weight_alpha = alpha
250 : END IF
251 13282 : weight_s => grid_atom%gapw_weight_s
252 : ELSE
253 32754 : weight_h => grid_atom%weight
254 32754 : weight_s => grid_atom%weight
255 : END IF
256 :
257 : ! create a place where to put the derivatives
258 46036 : 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 46036 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
262 : CALL xc_rho_set_create(rho_set_s, bounds, rho_cutoff=density_cut, &
263 46036 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
264 :
265 : ! allocate the required 3d arrays where to store rho and drho
266 46036 : CALL xc_rho_set_atom_update(rho_set_h, needs, nspins, bounds)
267 46036 : CALL xc_rho_set_atom_update(rho_set_s, needs, nspins, bounds)
268 :
269 46036 : CALL reallocate(rho_h, 1, na, 1, nr, 1, nspins)
270 46036 : CALL reallocate(rho_s, 1, na, 1, nr, 1, nspins)
271 46036 : CALL reallocate(vxc_h, 1, na, 1, nr, 1, nspins)
272 46036 : CALL reallocate(vxc_s, 1, na, 1, nr, 1, nspins)
273 : !
274 46036 : IF (gradient_f) THEN
275 30016 : CALL reallocate(drho_h, 1, 4, 1, na, 1, nr, 1, nspins)
276 30016 : CALL reallocate(drho_s, 1, 4, 1, na, 1, nr, 1, nspins)
277 30016 : CALL reallocate(vxg_h, 1, 3, 1, na, 1, nr, 1, nspins)
278 30016 : CALL reallocate(vxg_s, 1, 3, 1, na, 1, nr, 1, nspins)
279 : END IF
280 :
281 46036 : IF (tau_f) THEN
282 1218 : CALL create_tau_basis_cache(tau_basis_cache, grid_atom, basis_1c, harmonics)
283 1218 : CALL reallocate(tau_h, 1, na, 1, nr, 1, nspins)
284 1218 : CALL reallocate(tau_s, 1, na, 1, nr, 1, nspins)
285 1218 : CALL reallocate(vtau_h, 1, na, 1, nr, 1, nspins)
286 1218 : 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 46036 : donlcc = .FALSE.
291 46036 : IF (nlcc) THEN
292 458 : NULLIFY (rho_nlcc)
293 458 : rho_nlcc => my_kind_set(ikind)%nlcc_pot
294 458 : IF (ASSOCIATED(rho_nlcc)) donlcc = .TRUE.
295 : END IF
296 :
297 : ! Distribute the atoms of this kind
298 :
299 46036 : num_pe = para_env%num_pe
300 46036 : bo = get_limit(natom, para_env%num_pe, para_env%mepos)
301 :
302 82285 : DO iat = bo(1), bo(2)
303 36249 : iatom = atom_list(iat)
304 :
305 36249 : my_rho_atom_set(iatom)%exc_h = 0.0_dp
306 36249 : my_rho_atom_set(iatom)%exc_s = 0.0_dp
307 :
308 36249 : rho_atom => my_rho_atom_set(iatom)
309 117689623 : rho_h = 0.0_dp
310 117689623 : rho_s = 0.0_dp
311 36249 : IF (gradient_f) THEN
312 22957 : 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 22957 : rho_rad_s_d=r_s_d)
317 335985922 : drho_h = 0.0_dp
318 335985922 : drho_s = 0.0_dp
319 : ELSE
320 13292 : NULLIFY (r_h, r_s)
321 13292 : CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s)
322 13292 : rho_d = 0.0_dp
323 : END IF
324 36249 : IF (tau_f) THEN
325 : !compute tau on the grid all at once
326 849 : CALL calc_tau_atom(tau_h, tau_s, rho_atom, tau_basis_cache, nspins)
327 : ELSE
328 35400 : tau_d = 0.0_dp
329 : END IF
330 :
331 2019239 : 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 1982990 : r_h_d, r_s_d, drho_h, drho_s)
335 2019239 : IF (donlcc) THEN
336 : CALL calc_rho_nlcc(grid_atom, nspins, gradient_f, &
337 8650 : ir, rho_nlcc(:, 1), rho_h, rho_s, rho_nlcc(:, 2), drho_h, drho_s)
338 : END IF
339 : END DO
340 :
341 2019239 : DO ir = 1, nr
342 2019239 : IF (tau_f) THEN
343 43750 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_h, na, ir)
344 43750 : CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_s, na, ir)
345 1939240 : ELSE IF (gradient_f) THEN
346 1121540 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_d, na, ir)
347 1121540 : CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_d, na, ir)
348 : ELSE
349 817700 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, rho_d, tau_d, na, ir)
350 817700 : 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 36249 : 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 36249 : adiabatic_rescale_factor=my_adiabatic_rescale_factor)
361 36249 : rho_atom%exc_h = rho_atom%exc_h + exc_h
362 :
363 : !-------------------!
364 : ! soft atom density !
365 : !-------------------!
366 36249 : 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 36249 : adiabatic_rescale_factor=my_adiabatic_rescale_factor)
370 36249 : rho_atom%exc_s = rho_atom%exc_s + exc_s
371 :
372 : ! Add contributions to the exc energy
373 :
374 36249 : 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 36249 : IF (.NOT. energy_only) THEN
381 34737 : NULLIFY (int_hh, int_ss)
382 34737 : CALL get_rho_atom(rho_atom=rho_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss)
383 34737 : IF (gradient_f) THEN
384 : CALL gaVxcgb_GC(vxc_h, vxc_s, vxg_h, vxg_s, int_hh, int_ss, &
385 21627 : grid_atom, basis_1c, harmonics, nspins)
386 : ELSE
387 : CALL gaVxcgb_noGC(vxc_h, vxc_s, int_hh, int_ss, &
388 13110 : grid_atom, basis_1c, harmonics, nspins)
389 : END IF
390 34737 : IF (tau_f) THEN
391 : CALL dgaVtaudgb(vtau_h, vtau_s, int_hh, int_ss, &
392 849 : tau_basis_cache, nspins)
393 : END IF
394 : END IF ! energy_only
395 82285 : NULLIFY (r_h, r_s, dr_h, dr_s)
396 : END DO ! iat
397 :
398 46036 : IF (tau_f) CALL release_tau_basis_cache(tau_basis_cache)
399 :
400 : ! Release the xc structure used to store the xc derivatives
401 46036 : CALL xc_dset_release(deriv_set)
402 46036 : CALL xc_rho_set_release(rho_set_h)
403 167396 : CALL xc_rho_set_release(rho_set_s)
404 : END DO ! ikind
405 :
406 25276 : CALL para_env%sum(exc1)
407 :
408 25276 : IF (ASSOCIATED(rho_h)) DEALLOCATE (rho_h)
409 25276 : IF (ASSOCIATED(rho_s)) DEALLOCATE (rho_s)
410 25276 : IF (ASSOCIATED(vxc_h)) DEALLOCATE (vxc_h)
411 25276 : IF (ASSOCIATED(vxc_s)) DEALLOCATE (vxc_s)
412 :
413 25276 : IF (gradient_f) THEN
414 16650 : IF (ASSOCIATED(drho_h)) DEALLOCATE (drho_h)
415 16650 : IF (ASSOCIATED(drho_s)) DEALLOCATE (drho_s)
416 16650 : IF (ASSOCIATED(vxg_h)) DEALLOCATE (vxg_h)
417 16650 : IF (ASSOCIATED(vxg_s)) DEALLOCATE (vxg_s)
418 : END IF
419 :
420 25276 : IF (tau_f) THEN
421 670 : IF (ASSOCIATED(tau_h)) DEALLOCATE (tau_h)
422 670 : IF (ASSOCIATED(tau_s)) DEALLOCATE (tau_s)
423 670 : IF (ASSOCIATED(vtau_h)) DEALLOCATE (vtau_h)
424 670 : IF (ASSOCIATED(vtau_s)) DEALLOCATE (vtau_s)
425 : END IF
426 :
427 : END IF !xc_none
428 :
429 29294 : CALL timestop(handle)
430 :
431 1113172 : 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 10 : 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 10 : 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 20 : REAL(dp), DIMENSION(:, :), POINTER :: rho_nlcc, weight_h, weight_s
459 20 : REAL(dp), DIMENSION(:, :, :), POINTER :: rho_h, rho_s, tau_h, tau_s, vtau_h, &
460 10 : vtau_s, vxc_h, vxc_s
461 20 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: drho_h, drho_s, vxg_h, vxg_s
462 10 : 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 10 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: my_kind_set
469 10 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s, int_hh, int_ss, r_h, r_s
470 10 : TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r_h_d, r_s_d
471 10 : 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 10 : 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 10 : CALL timeset(routineN, handle)
482 :
483 10 : NULLIFY (atom_list)
484 10 : NULLIFY (my_kind_set)
485 10 : NULLIFY (atomic_kind_set)
486 10 : NULLIFY (grid_atom)
487 10 : NULLIFY (harmonics)
488 10 : NULLIFY (input)
489 10 : NULLIFY (para_env)
490 10 : NULLIFY (rho_atom)
491 10 : NULLIFY (my_rho_atom_set)
492 10 : 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 10 : rho_atom_set=my_rho_atom_set)
501 :
502 10 : nlcc = has_nlcc(my_kind_set)
503 10 : accint = dft_control%qs_control%gapw_control%accurate_xcint
504 :
505 : my_xc_section => section_vals_get_subs_vals(input, &
506 10 : "PROPERTIES%LINRES%EPR%PRINT%G_TENSOR%XC")
507 10 : 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 10 : i_val=myfun)
510 :
511 10 : IF (myfun == xc_none) THEN
512 0 : exc1 = 0.0_dp
513 0 : my_rho_atom_set(:)%exc_h = 0.0_dp
514 0 : my_rho_atom_set(:)%exc_s = 0.0_dp
515 : ELSE
516 : CALL section_vals_val_get(my_xc_section, "DENSITY_CUTOFF", &
517 10 : r_val=density_cut)
518 : CALL section_vals_val_get(my_xc_section, "GRADIENT_CUTOFF", &
519 10 : r_val=gradient_cut)
520 : CALL section_vals_val_get(my_xc_section, "TAU_CUTOFF", &
521 10 : r_val=tau_cut)
522 :
523 10 : lsd = dft_control%lsd
524 10 : nspins = dft_control%nspins
525 : needs = xc_functionals_get_needs(xc_fun_section, &
526 : lsd=lsd, &
527 10 : calc_potential=.TRUE.)
528 :
529 : ! whatever the xc, if epr_xc, drho_spin is needed
530 10 : needs%drho_spin = .TRUE.
531 :
532 10 : gradient_f = (needs%drho .OR. needs%drho_spin)
533 10 : tau_f = (needs%tau .OR. needs%tau_spin)
534 :
535 : ! Initialize energy contribution from the one center XC terms to zero
536 10 : exc1 = 0.0_dp
537 :
538 : ! Nullify some pointers for work-arrays
539 10 : NULLIFY (rho_h, drho_h, rho_s, drho_s, weight_h, weight_s)
540 10 : NULLIFY (vxc_h, vxc_s, vxg_h, vxg_s)
541 10 : NULLIFY (tau_h, tau_s)
542 10 : NULLIFY (vtau_h, vtau_s)
543 :
544 : ! Here starts the loop over all the atoms
545 :
546 30 : DO ikind = 1, SIZE(atomic_kind_set)
547 20 : 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 20 : harmonics=harmonics, grid_atom=grid_atom)
550 20 : CALL get_qs_kind(my_kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
551 :
552 20 : IF (.NOT. paw_atom) CYCLE
553 :
554 20 : nr = grid_atom%nr
555 20 : 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 200 : bounds(1:2, 1:3) = 1
563 20 : bounds(2, 1) = na
564 20 : bounds(2, 2) = nr
565 :
566 : ! set integration weights
567 20 : IF (accint) THEN
568 0 : weight_h => grid_atom%weight
569 0 : alpha = dft_control%qs_control%gapw_control%aw(ikind)
570 0 : IF (ASSOCIATED(grid_atom%gapw_weight_s)) THEN
571 0 : IF (grid_atom%gapw_weight_alpha /= alpha) DEALLOCATE (grid_atom%gapw_weight_s)
572 : END IF
573 0 : IF (.NOT. ASSOCIATED(grid_atom%gapw_weight_s)) THEN
574 0 : ALLOCATE (grid_atom%gapw_weight_s(na, nr))
575 0 : DO ir = 1, nr
576 0 : agr = 1.0_dp - EXP(-alpha*grid_atom%rad2(ir))
577 0 : grid_atom%gapw_weight_s(:, ir) = grid_atom%weight(:, ir)*agr
578 : END DO
579 0 : grid_atom%gapw_weight_alpha = alpha
580 : END IF
581 0 : weight_s => grid_atom%gapw_weight_s
582 : ELSE
583 20 : weight_h => grid_atom%weight
584 20 : weight_s => grid_atom%weight
585 : END IF
586 :
587 : ! create a place where to put the derivatives
588 20 : 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 20 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
592 : CALL xc_rho_set_create(rho_set_s, bounds, rho_cutoff=density_cut, &
593 20 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
594 :
595 : ! allocate the required 3d arrays where to store rho and drho
596 20 : CALL xc_rho_set_atom_update(rho_set_h, needs, nspins, bounds)
597 20 : CALL xc_rho_set_atom_update(rho_set_s, needs, nspins, bounds)
598 :
599 20 : CALL reallocate(rho_h, 1, na, 1, nr, 1, nspins)
600 20 : CALL reallocate(rho_s, 1, na, 1, nr, 1, nspins)
601 20 : CALL reallocate(vxc_h, 1, na, 1, nr, 1, nspins)
602 20 : CALL reallocate(vxc_s, 1, na, 1, nr, 1, nspins)
603 : !
604 : IF (gradient_f) THEN
605 20 : CALL reallocate(drho_h, 1, 4, 1, na, 1, nr, 1, nspins)
606 20 : CALL reallocate(drho_s, 1, 4, 1, na, 1, nr, 1, nspins)
607 20 : CALL reallocate(vxg_h, 1, 3, 1, na, 1, nr, 1, nspins)
608 20 : CALL reallocate(vxg_s, 1, 3, 1, na, 1, nr, 1, nspins)
609 : END IF
610 :
611 20 : IF (tau_f) THEN
612 0 : CALL create_tau_basis_cache(tau_basis_cache, grid_atom, basis_1c, harmonics)
613 0 : CALL reallocate(tau_h, 1, na, 1, nr, 1, nspins)
614 0 : CALL reallocate(tau_s, 1, na, 1, nr, 1, nspins)
615 0 : CALL reallocate(vtau_h, 1, na, 1, nr, 1, nspins)
616 0 : 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 20 : donlcc = .FALSE.
621 20 : IF (nlcc) THEN
622 0 : NULLIFY (rho_nlcc)
623 0 : rho_nlcc => my_kind_set(ikind)%nlcc_pot
624 0 : IF (ASSOCIATED(rho_nlcc)) donlcc = .TRUE.
625 : END IF
626 :
627 : ! Distribute the atoms of this kind
628 :
629 20 : num_pe = para_env%num_pe
630 20 : bo = get_limit(natom, para_env%num_pe, para_env%mepos)
631 :
632 35 : DO iat = bo(1), bo(2)
633 15 : iatom = atom_list(iat)
634 :
635 15 : my_rho_atom_set(iatom)%exc_h = 0.0_dp
636 15 : my_rho_atom_set(iatom)%exc_s = 0.0_dp
637 :
638 15 : rho_atom => my_rho_atom_set(iatom)
639 76545 : rho_h = 0.0_dp
640 76545 : rho_s = 0.0_dp
641 : IF (gradient_f) THEN
642 15 : 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 15 : rho_rad_s_d=r_s_d)
647 376545 : drho_h = 0.0_dp
648 376545 : 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 15 : IF (tau_f) THEN
655 : !compute tau on the grid all at once
656 0 : CALL calc_tau_atom(tau_h, tau_s, rho_atom, tau_basis_cache, nspins)
657 : ELSE
658 15 : tau_d = 0.0_dp
659 : END IF
660 :
661 765 : 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 750 : r_h_d, r_s_d, drho_h, drho_s)
665 765 : IF (donlcc) THEN
666 : CALL calc_rho_nlcc(grid_atom, nspins, gradient_f, &
667 0 : ir, rho_nlcc(:, 1), rho_h, rho_s, rho_nlcc(:, 2), drho_h, drho_s)
668 : END IF
669 : END DO
670 765 : DO ir = 1, nr
671 765 : IF (tau_f) THEN
672 0 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_h, na, ir)
673 0 : 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 750 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_d, na, ir)
676 750 : 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 15 : 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 15 : lsd, na, nr, exc_h, vxc_h, vxg_h, vtau_h)
689 15 : rho_atom%exc_h = rho_atom%exc_h + exc_h
690 :
691 : !-------------------!
692 : ! soft atom density !
693 : !-------------------!
694 15 : 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 15 : lsd, na, nr, exc_s, vxc_s, vxg_s, vtau_s)
697 15 : rho_atom%exc_s = rho_atom%exc_s + exc_s
698 :
699 45 : DO ispin = 1, nspins
700 135 : DO idir = 1, 3
701 4620 : DO ir = 1, nr
702 229590 : 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 225000 : + 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 229500 : + 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 15 : 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 15 : NULLIFY (int_hh, int_ss)
723 15 : 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 15 : 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 15 : IF (tau_f) THEN
732 : CALL dgaVtaudgb(vtau_h, vtau_s, int_hh, int_ss, &
733 0 : tau_basis_cache, nspins)
734 : END IF
735 35 : NULLIFY (r_h, r_s, dr_h, dr_s)
736 : END DO ! iat
737 :
738 20 : IF (tau_f) CALL release_tau_basis_cache(tau_basis_cache)
739 :
740 : ! Release the xc structure used to store the xc derivatives
741 20 : CALL xc_dset_release(deriv_set)
742 20 : CALL xc_rho_set_release(rho_set_h)
743 70 : CALL xc_rho_set_release(rho_set_s)
744 : END DO ! ikind
745 :
746 10 : CALL para_env%sum(exc1)
747 :
748 10 : IF (ASSOCIATED(rho_h)) DEALLOCATE (rho_h)
749 10 : IF (ASSOCIATED(rho_s)) DEALLOCATE (rho_s)
750 10 : IF (ASSOCIATED(vxc_h)) DEALLOCATE (vxc_h)
751 10 : IF (ASSOCIATED(vxc_s)) DEALLOCATE (vxc_s)
752 :
753 : IF (gradient_f) THEN
754 10 : IF (ASSOCIATED(drho_h)) DEALLOCATE (drho_h)
755 10 : IF (ASSOCIATED(drho_s)) DEALLOCATE (drho_s)
756 10 : IF (ASSOCIATED(vxg_h)) DEALLOCATE (vxg_h)
757 10 : IF (ASSOCIATED(vxg_s)) DEALLOCATE (vxg_s)
758 : END IF
759 :
760 10 : IF (tau_f) THEN
761 0 : IF (ASSOCIATED(tau_h)) DEALLOCATE (tau_h)
762 0 : IF (ASSOCIATED(tau_s)) DEALLOCATE (tau_s)
763 0 : IF (ASSOCIATED(vtau_h)) DEALLOCATE (vtau_h)
764 0 : IF (ASSOCIATED(vtau_s)) DEALLOCATE (vtau_s)
765 : END IF
766 :
767 : END IF !xc_none
768 :
769 10 : CALL timestop(handle)
770 :
771 380 : 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 5860 : 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 5860 : 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 5860 : 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 5860 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: weight_h, weight_s
812 17580 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: rho1_h, rho1_s, rho_h, rho_s, tau1_h, &
813 5860 : tau1_s, tau_h, tau_s
814 11720 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: drho1_h, drho1_s, drho_h, drho_s, vxg_h, &
815 5860 : vxg_s
816 5860 : 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 5860 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: my_kind_set, qs_kind_set
822 5860 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr1_h, dr1_s, dr_h, dr_s, int_hh, &
823 5860 : int_ss, r1_h, r1_s, r_h, r_s
824 5860 : 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 5860 : 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 5860 : CALL timeset(routineN, handle)
836 :
837 5860 : NULLIFY (qs_kind_set)
838 5860 : NULLIFY (rho_h, rho_s, drho_h, drho_s, weight_h, weight_s)
839 5860 : NULLIFY (rho1_h, rho1_s, drho1_h, drho1_s)
840 5860 : NULLIFY (vxc_h, vxc_s, vxg_h, vxg_s)
841 5860 : 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 5860 : atomic_kind_set=atomic_kind_set)
848 :
849 5860 : IF (PRESENT(kind_set_external)) THEN
850 616 : my_kind_set => kind_set_external
851 : ELSE
852 5244 : my_kind_set => qs_kind_set
853 : END IF
854 :
855 5860 : accint = dft_control%qs_control%gapw_control%accurate_xcint
856 :
857 5860 : CALL section_vals_val_get(input, "DFT%LSD", l_val=lsd)
858 : CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", &
859 5860 : r_val=density_cut)
860 : CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", &
861 5860 : r_val=gradient_cut)
862 : CALL section_vals_val_get(xc_section, "TAU_CUTOFF", &
863 5860 : r_val=tau_cut)
864 :
865 5860 : my_do_sf = .FALSE.
866 5860 : IF (PRESENT(do_sf)) my_do_sf = do_sf
867 :
868 : xc_fun_section => section_vals_get_subs_vals(xc_section, &
869 5860 : "XC_FUNCTIONAL")
870 5860 : IF (lsd) THEN
871 134 : nspins = 2
872 : ELSE
873 5726 : nspins = 1
874 : END IF
875 :
876 5860 : scale_rho = .FALSE.
877 5860 : IF (PRESENT(do_tddfpt2) .AND. PRESENT(do_triplet)) THEN
878 2876 : IF (nspins == 1 .AND. do_triplet) THEN
879 310 : lsd = .TRUE.
880 310 : scale_rho = .TRUE.
881 : END IF
882 2984 : ELSEIF (PRESENT(do_triplet)) THEN
883 2686 : IF (nspins == 1 .AND. do_triplet) lsd = .TRUE.
884 : END IF
885 :
886 : needs = xc_functionals_get_needs(xc_fun_section, lsd=lsd, &
887 5860 : calc_potential=.TRUE.)
888 5860 : gradient_functional = needs%drho .OR. needs%drho_spin
889 5860 : tau_f = (needs%tau .OR. needs%tau_spin)
890 5860 : IF (.NOT. tau_f) rtau = 0.0_dp
891 :
892 : ! Here starts the loop over all the atoms
893 18638 : DO ikind = 1, SIZE(atomic_kind_set)
894 :
895 12778 : NULLIFY (atom_list, harmonics, grid_atom)
896 12778 : 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 12778 : harmonics=harmonics, grid_atom=grid_atom)
899 12778 : CALL get_qs_kind(my_kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
900 12778 : IF (.NOT. paw_atom) CYCLE
901 :
902 12070 : nr = grid_atom%nr
903 12070 : na = grid_atom%ng_sphere
904 :
905 : ! set integration weights
906 12070 : IF (accint) THEN
907 4772 : weight_h => grid_atom%weight
908 4772 : alpha = dft_control%qs_control%gapw_control%aw(ikind)
909 4772 : IF (ASSOCIATED(grid_atom%gapw_weight_s)) THEN
910 4768 : IF (grid_atom%gapw_weight_alpha /= alpha) DEALLOCATE (grid_atom%gapw_weight_s)
911 : END IF
912 4772 : IF (.NOT. ASSOCIATED(grid_atom%gapw_weight_s)) THEN
913 16 : ALLOCATE (grid_atom%gapw_weight_s(na, nr))
914 204 : DO ir = 1, nr
915 200 : agr = 1.0_dp - EXP(-alpha*grid_atom%rad2(ir))
916 10204 : grid_atom%gapw_weight_s(:, ir) = grid_atom%weight(:, ir)*agr
917 : END DO
918 4 : grid_atom%gapw_weight_alpha = alpha
919 : END IF
920 4772 : weight_s => grid_atom%gapw_weight_s
921 : ELSE
922 7298 : weight_h => grid_atom%weight
923 7298 : 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 120700 : bounds(1:2, 1:3) = 1
930 12070 : bounds(2, 1) = na
931 12070 : bounds(2, 2) = nr
932 :
933 12070 : CALL xc_dset_create(deriv_set, local_bounds=bounds)
934 : CALL xc_rho_set_create(rho_set_h, bounds, rho_cutoff=density_cut, &
935 12070 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
936 : CALL xc_rho_set_create(rho_set_s, bounds, rho_cutoff=density_cut, &
937 12070 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
938 : CALL xc_rho_set_create(rho1_set_h, bounds, rho_cutoff=density_cut, &
939 12070 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
940 : CALL xc_rho_set_create(rho1_set_s, bounds, rho_cutoff=density_cut, &
941 12070 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
942 :
943 : ! allocate the required 3d arrays where to store rho and drho
944 12070 : IF (nspins == 1 .AND. .NOT. lsd) THEN
945 11292 : CALL xc_rho_set_atom_update(rho_set_h, needs, 1, bounds)
946 11292 : CALL xc_rho_set_atom_update(rho1_set_h, needs, 1, bounds)
947 11292 : CALL xc_rho_set_atom_update(rho_set_s, needs, 1, bounds)
948 11292 : CALL xc_rho_set_atom_update(rho1_set_s, needs, 1, bounds)
949 : ELSE
950 778 : CALL xc_rho_set_atom_update(rho_set_h, needs, 2, bounds)
951 778 : CALL xc_rho_set_atom_update(rho1_set_h, needs, 2, bounds)
952 778 : CALL xc_rho_set_atom_update(rho_set_s, needs, 2, bounds)
953 778 : 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 168980 : rho_s(1:na, 1:nr, 1:nspins), rho1_s(1:na, 1:nr, 1:nspins))
958 :
959 84490 : ALLOCATE (vxc_h(1:na, 1:nr, 1:nspins), vxc_s(1:na, 1:nr, 1:nspins))
960 31404676 : vxc_h = 0.0_dp
961 31404676 : vxc_s = 0.0_dp
962 :
963 12070 : IF (tau_f) THEN
964 0 : 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 0 : tau_s(1:na, 1:nr, 1:nspins), tau1_s(1:na, 1:nr, 1:nspins))
967 0 : ALLOCATE (vtau_h(1:na, 1:nr, 1:nspins), vtau_s(1:na, 1:nr, 1:nspins))
968 : END IF
969 :
970 12070 : 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 120904 : drho_s(1:4, 1:na, 1:nr, 1:nspins), drho1_s(1:4, 1:na, 1:nr, 1:nspins))
973 69088 : 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 3434 : drho_s(1, 1, 1, 1), drho1_s(1, 1, 1, 1))
977 3434 : ALLOCATE (vxg_h(1, 1, 1, 1), vxg_s(1, 1, 1, 1))
978 3434 : rrho = 0.0_dp
979 : END IF
980 88997258 : vxg_h = 0.0_dp
981 88997258 : vxg_s = 0.0_dp
982 :
983 : ! parallelization
984 12070 : local_loop_limit = get_limit(natom, para_env%num_pe, para_env%mepos)
985 :
986 20660 : DO iatom = local_loop_limit(1), local_loop_limit(2) !1,natom
987 8590 : atom = atom_list(iatom)
988 :
989 8590 : rho_atom_set(atom)%exc_h = 0.0_dp
990 8590 : rho_atom_set(atom)%exc_s = 0.0_dp
991 8590 : rho1_atom_set(atom)%exc_h = 0.0_dp
992 8590 : rho1_atom_set(atom)%exc_s = 0.0_dp
993 :
994 8590 : rho_atom => rho_atom_set(atom)
995 8590 : rho1_atom => rho1_atom_set(atom)
996 8590 : NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
997 8590 : NULLIFY (r1_h, r1_s, dr1_h, dr1_s, r1_h_d, r1_s_d)
998 22416574 : rho_h = 0.0_dp
999 22416574 : rho_s = 0.0_dp
1000 22416574 : rho1_h = 0.0_dp
1001 22416574 : rho1_s = 0.0_dp
1002 8590 : 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 6178 : 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 6178 : rho_rad_h_d=r1_h_d, rho_rad_s_d=r1_s_d)
1011 159705102 : drho_h = 0.0_dp; drho_s = 0.0_dp
1012 159705102 : drho1_h = 0.0_dp; drho1_s = 0.0_dp
1013 : ELSE
1014 : CALL get_rho_atom(rho_atom=rho_atom, &
1015 2412 : rho_rad_h=r_h, rho_rad_s=r_s)
1016 : CALL get_rho_atom(rho_atom=rho1_atom, &
1017 2412 : rho_rad_h=r1_h, rho_rad_s=r1_s)
1018 : END IF
1019 :
1020 8590 : rtot = 0.0_dp
1021 :
1022 438090 : 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 429500 : 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 438090 : drho1_h, drho1_s)
1029 : END DO
1030 8590 : IF (tau_f) THEN
1031 0 : CALL calc_tau_atom(tau_h, tau_s, rho_atom, tau_basis_cache, nspins)
1032 0 : CALL calc_tau_atom(tau1_h, tau1_s, rho1_atom, tau_basis_cache, nspins)
1033 : END IF
1034 8590 : IF (scale_rho) THEN
1035 926376 : rho_h = 2.0_dp*rho_h
1036 926376 : rho_s = 2.0_dp*rho_s
1037 363 : IF (gradient_functional) THEN
1038 3426696 : drho_h = 2.0_dp*drho_h
1039 3426696 : drho_s = 2.0_dp*drho_s
1040 : END IF
1041 363 : IF (tau_f) THEN
1042 0 : tau_h = 2.0_dp*tau_h
1043 0 : tau_s = 2.0_dp*tau_s
1044 : END IF
1045 : END IF
1046 :
1047 438090 : DO ir = 1, nr
1048 438090 : IF (tau_f) THEN
1049 0 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_h, na, ir)
1050 0 : CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, drho1_h, tau1_h, na, ir)
1051 0 : CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_s, na, ir)
1052 0 : CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, drho1_s, tau1_s, na, ir)
1053 429500 : ELSE IF (gradient_functional) THEN
1054 308900 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, rtau, na, ir)
1055 308900 : CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, drho1_h, rtau, na, ir)
1056 308900 : CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, rtau, na, ir)
1057 308900 : CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, drho1_s, rtau, na, ir)
1058 : ELSE
1059 120600 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, rrho, rtau, na, ir)
1060 120600 : CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, rrho, rtau, na, ir)
1061 120600 : CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, rrho, rtau, na, ir)
1062 120600 : 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 8590 : 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 8590 : do_sf=my_do_sf)
1076 :
1077 8590 : CALL get_rho_atom(rho_atom=rho1_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss)
1078 8590 : IF (gradient_functional) THEN
1079 : CALL gaVxcgb_GC(vxc_h, vxc_s, vxg_h, vxg_s, int_hh, int_ss, &
1080 6178 : grid_atom, basis_1c, harmonics, nspins)
1081 : ELSE
1082 : CALL gaVxcgb_noGC(vxc_h, vxc_s, int_hh, int_ss, &
1083 2412 : grid_atom, basis_1c, harmonics, nspins)
1084 : END IF
1085 8590 : IF (tau_f) THEN
1086 : CALL dgaVtaudgb(vtau_h, vtau_s, int_hh, int_ss, &
1087 0 : tau_basis_cache, nspins)
1088 : END IF
1089 :
1090 20660 : NULLIFY (r_h, r_s, dr_h, dr_s)
1091 :
1092 : END DO
1093 :
1094 : ! some cleanup
1095 12070 : DEALLOCATE (rho_h, rho_s, rho1_h, rho1_s, vxc_h, vxc_s)
1096 12070 : DEALLOCATE (drho_h, drho_s, vxg_h, vxg_s)
1097 12070 : DEALLOCATE (drho1_h, drho1_s)
1098 12070 : IF (tau_f) THEN
1099 0 : DEALLOCATE (tau_h, tau_s, tau1_h, tau1_s)
1100 0 : DEALLOCATE (vtau_h, vtau_s)
1101 0 : CALL release_tau_basis_cache(tau_basis_cache)
1102 : END IF
1103 :
1104 12070 : CALL xc_dset_release(deriv_set)
1105 12070 : CALL xc_rho_set_release(rho_set_h)
1106 12070 : CALL xc_rho_set_release(rho1_set_h)
1107 12070 : CALL xc_rho_set_release(rho_set_s)
1108 43486 : CALL xc_rho_set_release(rho1_set_s)
1109 : END DO
1110 :
1111 5860 : CALL timestop(handle)
1112 :
1113 433640 : 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 0 : 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 0 : 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 0 : REAL(dp), DIMENSION(:, :), POINTER :: rho_nlcc, weight_h, weight_s
1153 0 : REAL(dp), DIMENSION(:, :, :), POINTER :: rho0_h, rho0_s, rho1_h, rho1_s, rho_h, &
1154 0 : rho_s, tau0_h, tau0_s, tau1_h, tau1_s, &
1155 0 : tau_h, tau_s, vtau_h, vtau_s, vxc_h, &
1156 0 : vxc_s
1157 0 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: drho0_h, drho0_s, drho1_h, drho1_s, &
1158 0 : drho_h, drho_s, vxg_h, vxg_s
1159 : REAL(KIND=dp), DIMENSION(-4:4) :: ak, bl
1160 0 : 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 0 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s, fint_hh, fint_ss, int_hh, &
1167 0 : int_ss, r_h, r_s
1168 0 : 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 0 : 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 0 : CALL timeset(routineN, handle)
1177 :
1178 0 : NULLIFY (vtau_h, vtau_s)
1179 :
1180 0 : ak = 0.0_dp
1181 0 : bl = 0.0_dp
1182 0 : SELECT CASE (accuracy)
1183 : CASE (:4)
1184 0 : nstep = 2
1185 0 : ak(-2:2) = [1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp]/12.0_dp
1186 0 : 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 0 : nstep = 3
1189 0 : 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 0 : 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 0 : 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 0 : 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 0 : 896.0_dp, -112.0_dp, 128.0_dp/9.0_dp, -1.0_dp]/560.0_dp
1197 : END SELECT
1198 0 : oeps1 = 1.0_dp/epsrho
1199 0 : 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 0 : atomic_kind_set=atomic_kind_set)
1205 :
1206 0 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
1207 0 : CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
1208 :
1209 0 : accint = dft_control%qs_control%gapw_control%accurate_xcint
1210 :
1211 0 : IF (myfun == xc_none) THEN
1212 : ! no action needed?
1213 : ELSE
1214 0 : CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", r_val=density_cut)
1215 0 : CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", r_val=gradient_cut)
1216 0 : CALL section_vals_val_get(xc_section, "TAU_CUTOFF", r_val=tau_cut)
1217 :
1218 0 : nlcc = has_nlcc(kind_set)
1219 0 : lsd = dft_control%lsd
1220 0 : nspins = dft_control%nspins
1221 0 : mspins = nspins
1222 0 : IF (is_triplet) THEN
1223 0 : CPASSERT(nspins == 1)
1224 0 : lsd = .TRUE.
1225 0 : mspins = 2
1226 : END IF
1227 0 : needs = xc_functionals_get_needs(xc_fun_section, lsd=lsd, calc_potential=.TRUE.)
1228 0 : gradient_f = (needs%drho .OR. needs%drho_spin)
1229 0 : tau_f = (needs%tau .OR. needs%tau_spin)
1230 :
1231 : ! Here starts the loop over all the atoms
1232 0 : DO ikind = 1, SIZE(atomic_kind_set)
1233 0 : 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 0 : harmonics=harmonics, grid_atom=grid_atom)
1236 0 : CALL get_qs_kind(kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
1237 :
1238 0 : IF (.NOT. paw_atom) CYCLE
1239 :
1240 0 : nr = grid_atom%nr
1241 0 : na = grid_atom%ng_sphere
1242 :
1243 : ! set integration weights
1244 0 : IF (accint) THEN
1245 0 : weight_h => grid_atom%weight
1246 0 : alpha = dft_control%qs_control%gapw_control%aw(ikind)
1247 0 : IF (ASSOCIATED(grid_atom%gapw_weight_s)) THEN
1248 0 : IF (grid_atom%gapw_weight_alpha /= alpha) DEALLOCATE (grid_atom%gapw_weight_s)
1249 : END IF
1250 0 : IF (.NOT. ASSOCIATED(grid_atom%gapw_weight_s)) THEN
1251 0 : ALLOCATE (grid_atom%gapw_weight_s(na, nr))
1252 0 : DO ir = 1, nr
1253 0 : agr = 1.0_dp - EXP(-alpha*grid_atom%rad2(ir))
1254 0 : grid_atom%gapw_weight_s(:, ir) = grid_atom%weight(:, ir)*agr
1255 : END DO
1256 0 : grid_atom%gapw_weight_alpha = alpha
1257 : END IF
1258 0 : weight_s => grid_atom%gapw_weight_s
1259 : ELSE
1260 0 : weight_h => grid_atom%weight
1261 0 : 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 0 : bounds(1:2, 1:3) = 1
1270 0 : bounds(2, 1) = na
1271 0 : bounds(2, 2) = nr
1272 :
1273 : ! create a place where to put the derivatives
1274 0 : 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 0 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
1278 : CALL xc_rho_set_create(rho_set_s, bounds, rho_cutoff=density_cut, &
1279 0 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
1280 :
1281 : ! allocate the required 3d arrays where to store rho and drho
1282 0 : CALL xc_rho_set_atom_update(rho_set_h, needs, mspins, bounds)
1283 0 : 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 0 : rho1_h(na, nr, nspins), rho1_s(na, nr, nspins))
1288 0 : ALLOCATE (vxc_h(na, nr, mspins), vxc_s(na, nr, mspins))
1289 0 : 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 0 : drho1_h(4, na, nr, nspins), drho1_s(4, na, nr, nspins))
1293 0 : ALLOCATE (vxg_h(3, na, nr, mspins), vxg_s(3, na, nr, mspins))
1294 : END IF
1295 0 : IF (tau_f) THEN
1296 0 : 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 0 : tau1_h(na, nr, nspins), tau1_s(na, nr, nspins))
1300 0 : 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 0 : donlcc = .FALSE.
1305 0 : IF (nlcc) THEN
1306 0 : NULLIFY (rho_nlcc)
1307 0 : rho_nlcc => kind_set(ikind)%nlcc_pot
1308 0 : IF (ASSOCIATED(rho_nlcc)) donlcc = .TRUE.
1309 : END IF
1310 :
1311 : ! Distribute the atoms of this kind
1312 0 : num_pe = para_env%num_pe
1313 0 : bo = get_limit(natom, num_pe, para_env%mepos)
1314 :
1315 0 : DO iat = bo(1), bo(2)
1316 0 : iatom = atom_list(iat)
1317 : !
1318 0 : NULLIFY (int_hh, int_ss)
1319 0 : rho0_atom => rho0_atom_set(iatom)
1320 0 : CALL get_rho_atom(rho_atom=rho0_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss)
1321 0 : ALLOCATE (fint_ss(nspins), fint_hh(nspins))
1322 0 : DO ns = 1, nspins
1323 0 : nf = SIZE(int_ss(ns)%r_coef, 1)
1324 0 : ALLOCATE (fint_ss(ns)%r_coef(nf, nf))
1325 0 : nf = SIZE(int_hh(ns)%r_coef, 1)
1326 0 : ALLOCATE (fint_hh(ns)%r_coef(nf, nf))
1327 : END DO
1328 :
1329 : ! RHO0
1330 0 : rho0_h = 0.0_dp
1331 0 : rho0_s = 0.0_dp
1332 0 : rho0_atom => rho0_atom_set(iatom)
1333 0 : IF (gradient_f) THEN
1334 0 : 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 0 : drho_rad_s=dr_s, rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
1337 0 : drho0_h = 0.0_dp
1338 0 : drho0_s = 0.0_dp
1339 : ELSE
1340 0 : NULLIFY (r_h, r_s)
1341 0 : CALL get_rho_atom(rho_atom=rho0_atom, rho_rad_h=r_h, rho_rad_s=r_s)
1342 0 : rho_d = 0.0_dp
1343 : END IF
1344 0 : 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 0 : r_h_d, r_s_d, drho0_h, drho0_s)
1348 0 : IF (donlcc) THEN
1349 : CALL calc_rho_nlcc(grid_atom, nspins, gradient_f, &
1350 0 : ir, rho_nlcc(:, 1), rho0_h, rho0_s, rho_nlcc(:, 2), drho0_h, drho0_s)
1351 : END IF
1352 : END DO
1353 0 : IF (tau_f) THEN
1354 : !compute tau on the grid all at once
1355 0 : CALL calc_tau_atom(tau0_h, tau0_s, rho0_atom, tau_basis_cache, nspins)
1356 : ELSE
1357 0 : tau_d = 0.0_dp
1358 : END IF
1359 : ! RHO1
1360 0 : rho1_h = 0.0_dp
1361 0 : rho1_s = 0.0_dp
1362 0 : rho1_atom => rho1_atom_set(iatom)
1363 0 : IF (gradient_f) THEN
1364 0 : 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 0 : drho_rad_s=dr_s, rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
1367 0 : drho1_h = 0.0_dp
1368 0 : drho1_s = 0.0_dp
1369 : ELSE
1370 0 : NULLIFY (r_h, r_s)
1371 0 : CALL get_rho_atom(rho_atom=rho1_atom, rho_rad_h=r_h, rho_rad_s=r_s)
1372 : END IF
1373 0 : 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 0 : r_h_d, r_s_d, drho1_h, drho1_s)
1377 : END DO
1378 0 : IF (tau_f) THEN
1379 : !compute tau on the grid all at once
1380 0 : CALL calc_tau_atom(tau1_h, tau1_s, rho1_atom, tau_basis_cache, nspins)
1381 : END IF
1382 : ! RHO2
1383 0 : rho2_atom => rho2_atom_set(iatom)
1384 :
1385 0 : DO istep = -nstep, nstep
1386 :
1387 0 : beta = REAL(istep, KIND=dp)*epsrho
1388 :
1389 0 : IF (is_triplet) THEN
1390 0 : rho_h(:, :, 1) = rho0_h(:, :, 1) + beta*rho1_h(:, :, 1)
1391 0 : rho_h(:, :, 2) = rho0_h(:, :, 1)
1392 0 : rho_h = 0.5_dp*rho_h
1393 0 : rho_s(:, :, 1) = rho0_s(:, :, 1) + beta*rho1_s(:, :, 1)
1394 0 : rho_s(:, :, 2) = rho0_s(:, :, 1)
1395 0 : rho_s = 0.5_dp*rho_s
1396 0 : IF (gradient_f) THEN
1397 0 : drho_h(:, :, :, 1) = drho0_h(:, :, :, 1) + beta*drho1_h(:, :, :, 1)
1398 0 : drho_h(:, :, :, 2) = drho0_h(:, :, :, 1)
1399 0 : drho_h = 0.5_dp*drho_h
1400 0 : drho_s(:, :, :, 1) = drho0_s(:, :, :, 1) + beta*drho1_s(:, :, :, 1)
1401 0 : drho_s(:, :, :, 2) = drho0_s(:, :, :, 1)
1402 0 : drho_s = 0.5_dp*drho_s
1403 : END IF
1404 0 : IF (tau_f) THEN
1405 0 : tau_h(:, :, 1) = tau0_h(:, :, 1) + beta*tau1_h(:, :, 1)
1406 0 : tau_h(:, :, 2) = tau0_h(:, :, 1)
1407 0 : tau_h = 0.5_dp*tau0_h
1408 0 : tau_s(:, :, 1) = tau0_s(:, :, 1) + beta*tau1_s(:, :, 1)
1409 0 : tau_s(:, :, 2) = tau0_s(:, :, 1)
1410 0 : tau_s = 0.5_dp*tau0_s
1411 : END IF
1412 : ELSE
1413 0 : rho_h = rho0_h + beta*rho1_h
1414 0 : rho_s = rho0_s + beta*rho1_s
1415 0 : IF (gradient_f) THEN
1416 0 : drho_h = drho0_h + beta*drho1_h
1417 0 : drho_s = drho0_s + beta*drho1_s
1418 : END IF
1419 0 : IF (tau_f) THEN
1420 0 : tau_h = tau0_h + beta*tau1_h
1421 0 : tau_s = tau0_s + beta*tau1_s
1422 : END IF
1423 : END IF
1424 : !
1425 0 : IF (gradient_f) THEN
1426 : drho_h(4, :, :, :) = SQRT( &
1427 : drho_h(1, :, :, :)*drho_h(1, :, :, :) + &
1428 : drho_h(2, :, :, :)*drho_h(2, :, :, :) + &
1429 0 : 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 0 : drho_s(3, :, :, :)*drho_s(3, :, :, :))
1435 : END IF
1436 :
1437 0 : DO ir = 1, nr
1438 0 : IF (tau_f) THEN
1439 0 : CALL fill_rho_set(rho_set_h, lsd, mspins, needs, rho_h, drho_h, tau_h, na, ir)
1440 0 : CALL fill_rho_set(rho_set_s, lsd, mspins, needs, rho_s, drho_s, tau_s, na, ir)
1441 0 : ELSE IF (gradient_f) THEN
1442 0 : CALL fill_rho_set(rho_set_h, lsd, mspins, needs, rho_h, drho_h, tau_d, na, ir)
1443 0 : CALL fill_rho_set(rho_set_s, lsd, mspins, needs, rho_s, drho_s, tau_d, na, ir)
1444 : ELSE
1445 0 : CALL fill_rho_set(rho_set_h, lsd, mspins, needs, rho_h, rho_d, tau_d, na, ir)
1446 0 : 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 0 : 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 0 : lsd, na, nr, exc_h, vxc_h, vxg_h, vtau_h)
1454 0 : IF (is_triplet) THEN
1455 0 : vxc_h(:, :, 1) = vxc_h(:, :, 1) - vxc_h(:, :, 2)
1456 0 : IF (gradient_f) THEN
1457 0 : vxg_h(:, :, :, 1) = vxg_h(:, :, :, 1) - vxg_h(:, :, :, 2)
1458 : END IF
1459 0 : IF (tau_f) THEN
1460 0 : vtau_h(:, :, 1) = vtau_h(:, :, 1) - vtau_h(:, :, 2)
1461 : END IF
1462 : END IF
1463 : ! soft atom density !
1464 0 : 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 0 : lsd, na, nr, exc_s, vxc_s, vxg_s, vtau_s)
1467 0 : IF (is_triplet) THEN
1468 0 : vxc_s(:, :, 1) = vxc_s(:, :, 1) - vxc_s(:, :, 2)
1469 0 : IF (gradient_f) THEN
1470 0 : vxg_s(:, :, :, 1) = vxg_s(:, :, :, 1) - vxg_s(:, :, :, 2)
1471 : END IF
1472 0 : IF (tau_f) THEN
1473 0 : vtau_s(:, :, 1) = vtau_s(:, :, 1) - vtau_s(:, :, 2)
1474 : END IF
1475 : END IF
1476 : ! potentials
1477 0 : DO ns = 1, nspins
1478 0 : fint_hh(ns)%r_coef(:, :) = 0.0_dp
1479 0 : fint_ss(ns)%r_coef(:, :) = 0.0_dp
1480 : END DO
1481 0 : IF (gradient_f) THEN
1482 : CALL gaVxcgb_GC(vxc_h, vxc_s, vxg_h, vxg_s, fint_hh, fint_ss, &
1483 0 : grid_atom, basis_1c, harmonics, nspins)
1484 : ELSE
1485 : CALL gaVxcgb_noGC(vxc_h, vxc_s, fint_hh, fint_ss, &
1486 0 : grid_atom, basis_1c, harmonics, nspins)
1487 : END IF
1488 0 : IF (tau_f) THEN
1489 : CALL dgaVtaudgb(vtau_h, vtau_s, fint_hh, fint_ss, &
1490 0 : tau_basis_cache, nspins)
1491 : END IF
1492 : ! first derivative fxc
1493 0 : NULLIFY (int_hh, int_ss)
1494 0 : CALL get_rho_atom(rho_atom=rho1_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss)
1495 0 : DO ns = 1, nspins
1496 0 : int_ss(ns)%r_coef(:, :) = int_ss(ns)%r_coef(:, :) + oeps1*ak(istep)*fint_ss(ns)%r_coef(:, :)
1497 0 : 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 0 : NULLIFY (int_hh, int_ss)
1501 0 : CALL get_rho_atom(rho_atom=rho2_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss)
1502 0 : DO ns = 1, nspins
1503 0 : int_ss(ns)%r_coef(:, :) = int_ss(ns)%r_coef(:, :) + oeps2*bl(istep)*fint_ss(ns)%r_coef(:, :)
1504 0 : 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 0 : DO ns = 1, nspins
1509 0 : DEALLOCATE (fint_ss(ns)%r_coef)
1510 0 : DEALLOCATE (fint_hh(ns)%r_coef)
1511 : END DO
1512 0 : DEALLOCATE (fint_ss, fint_hh)
1513 :
1514 : END DO ! iat
1515 :
1516 : ! Release the xc structure used to store the xc derivatives
1517 0 : CALL xc_dset_release(deriv_set)
1518 0 : CALL xc_rho_set_release(rho_set_h)
1519 0 : CALL xc_rho_set_release(rho_set_s)
1520 :
1521 0 : DEALLOCATE (rho_h, rho_s, rho0_h, rho0_s, rho1_h, rho1_s)
1522 0 : DEALLOCATE (vxc_h, vxc_s)
1523 0 : IF (gradient_f) THEN
1524 0 : DEALLOCATE (drho_h, drho_s, drho0_h, drho0_s, drho1_h, drho1_s)
1525 0 : DEALLOCATE (vxg_h, vxg_s)
1526 : END IF
1527 0 : IF (tau_f) THEN
1528 0 : DEALLOCATE (tau_h, tau_s, tau0_h, tau0_s, tau1_h, tau1_s)
1529 0 : DEALLOCATE (vtau_h, vtau_s)
1530 0 : CALL release_tau_basis_cache(tau_basis_cache)
1531 : END IF
1532 : END DO ! ikind
1533 :
1534 : END IF !xc_none
1535 :
1536 0 : CALL timestop(handle)
1537 :
1538 0 : 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 114 : 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 114 : 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 114 : 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 228 : REAL(dp), DIMENSION(:, :), POINTER :: rho_nlcc, weight_h, weight_s
1579 228 : REAL(dp), DIMENSION(:, :, :), POINTER :: rho0_h, rho0_s, rho1_h, rho1_s, rho_h, &
1580 228 : rho_s, tau0_h, tau0_s, tau1_h, tau1_s, &
1581 114 : tau_h, tau_s
1582 114 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: drho0_h, drho0_s, drho1_h, drho1_s, &
1583 228 : drho_h, drho_s, vxg_h, vxg_s
1584 : REAL(KIND=dp), DIMENSION(-4:4) :: ak
1585 114 : 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 114 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s, fint_hh, fint_ss, int_hh, &
1592 114 : int_ss, r_h, r_s
1593 114 : 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 114 : 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 114 : CALL timeset(routineN, handle)
1603 :
1604 114 : NULLIFY (vtau_h, vtau_s)
1605 :
1606 114 : ak = 0.0_dp
1607 114 : SELECT CASE (accuracy)
1608 : CASE (:4)
1609 0 : nstep = 2
1610 0 : 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 912 : nstep = 3
1613 912 : 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 0 : 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 114 : 224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp]/280.0_dp
1618 : END SELECT
1619 114 : 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 114 : atomic_kind_set=atomic_kind_set)
1625 :
1626 114 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
1627 114 : CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
1628 :
1629 114 : accint = dft_control%qs_control%gapw_control%accurate_xcint
1630 :
1631 114 : 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 114 : do_triplet=is_triplet, kind_set_external=kind_set)
1637 :
1638 114 : CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", r_val=density_cut)
1639 114 : CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", r_val=gradient_cut)
1640 114 : CALL section_vals_val_get(xc_section, "TAU_CUTOFF", r_val=tau_cut)
1641 :
1642 114 : nlcc = has_nlcc(kind_set)
1643 114 : lsd = dft_control%lsd
1644 114 : nspins = dft_control%nspins
1645 114 : mspins = nspins
1646 114 : IF (is_triplet) THEN
1647 12 : CPASSERT(nspins == 1)
1648 12 : lsd = .TRUE.
1649 12 : mspins = 2
1650 : END IF
1651 114 : needs = xc_functionals_get_needs(xc_fun_section, lsd=lsd, calc_potential=.TRUE.)
1652 114 : gradient_f = (needs%drho .OR. needs%drho_spin)
1653 114 : tau_f = (needs%tau .OR. needs%tau_spin)
1654 :
1655 : ! Here starts the loop over all the atoms
1656 386 : DO ikind = 1, SIZE(atomic_kind_set)
1657 272 : 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 272 : harmonics=harmonics, grid_atom=grid_atom)
1660 272 : CALL get_qs_kind(kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
1661 :
1662 272 : IF (.NOT. paw_atom) CYCLE
1663 :
1664 258 : nr = grid_atom%nr
1665 258 : na = grid_atom%ng_sphere
1666 :
1667 : ! set integration weights
1668 258 : IF (accint) THEN
1669 138 : weight_h => grid_atom%weight
1670 138 : alpha = dft_control%qs_control%gapw_control%aw(ikind)
1671 138 : IF (ASSOCIATED(grid_atom%gapw_weight_s)) THEN
1672 138 : IF (grid_atom%gapw_weight_alpha /= alpha) DEALLOCATE (grid_atom%gapw_weight_s)
1673 : END IF
1674 138 : IF (.NOT. ASSOCIATED(grid_atom%gapw_weight_s)) THEN
1675 0 : ALLOCATE (grid_atom%gapw_weight_s(na, nr))
1676 0 : DO ir = 1, nr
1677 0 : agr = 1.0_dp - EXP(-alpha*grid_atom%rad2(ir))
1678 0 : grid_atom%gapw_weight_s(:, ir) = grid_atom%weight(:, ir)*agr
1679 : END DO
1680 0 : grid_atom%gapw_weight_alpha = alpha
1681 : END IF
1682 138 : weight_s => grid_atom%gapw_weight_s
1683 : ELSE
1684 120 : weight_h => grid_atom%weight
1685 120 : 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 2580 : bounds(1:2, 1:3) = 1
1694 258 : bounds(2, 1) = na
1695 258 : bounds(2, 2) = nr
1696 :
1697 : ! create a place where to put the derivatives
1698 258 : 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 258 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
1702 : CALL xc_rho_set_create(rho_set_s, bounds, rho_cutoff=density_cut, &
1703 258 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
1704 : CALL xc_rho_set_create(rho1_set_h, bounds, rho_cutoff=density_cut, &
1705 258 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
1706 : CALL xc_rho_set_create(rho1_set_s, bounds, rho_cutoff=density_cut, &
1707 258 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
1708 :
1709 : ! allocate the required 3d arrays where to store rho and drho
1710 258 : CALL xc_rho_set_atom_update(rho_set_h, needs, mspins, bounds)
1711 258 : CALL xc_rho_set_atom_update(rho_set_s, needs, mspins, bounds)
1712 258 : CALL xc_rho_set_atom_update(rho1_set_h, needs, mspins, bounds)
1713 258 : 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 5160 : rho1_h(na, nr, nspins), rho1_s(na, nr, nspins))
1718 1806 : ALLOCATE (vxc_h(na, nr, nspins), vxc_s(na, nr, nspins))
1719 258 : 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 3520 : drho1_h(4, na, nr, nspins), drho1_s(4, na, nr, nspins))
1723 1408 : ALLOCATE (vxg_h(3, na, nr, nspins), vxg_s(3, na, nr, nspins))
1724 : END IF
1725 258 : IF (tau_f) THEN
1726 0 : 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 0 : tau1_h(na, nr, nspins), tau1_s(na, nr, nspins))
1730 0 : 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 258 : donlcc = .FALSE.
1735 258 : IF (nlcc) THEN
1736 0 : NULLIFY (rho_nlcc)
1737 0 : rho_nlcc => kind_set(ikind)%nlcc_pot
1738 0 : IF (ASSOCIATED(rho_nlcc)) donlcc = .TRUE.
1739 : END IF
1740 :
1741 : ! Distribute the atoms of this kind
1742 258 : num_pe = para_env%num_pe
1743 258 : bo = get_limit(natom, num_pe, para_env%mepos)
1744 :
1745 441 : DO iat = bo(1), bo(2)
1746 183 : iatom = atom_list(iat)
1747 : !
1748 183 : NULLIFY (int_hh, int_ss)
1749 183 : rho0_atom => rho0_atom_set(iatom)
1750 183 : CALL get_rho_atom(rho_atom=rho0_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss)
1751 1098 : ALLOCATE (fint_ss(nspins), fint_hh(nspins))
1752 366 : DO ns = 1, nspins
1753 183 : nf = SIZE(int_ss(ns)%r_coef, 1)
1754 732 : ALLOCATE (fint_ss(ns)%r_coef(nf, nf))
1755 183 : nf = SIZE(int_hh(ns)%r_coef, 1)
1756 915 : ALLOCATE (fint_hh(ns)%r_coef(nf, nf))
1757 : END DO
1758 :
1759 : ! RHO0
1760 467016 : rho0_h = 0.0_dp
1761 467016 : rho0_s = 0.0_dp
1762 183 : rho0_atom => rho0_atom_set(iatom)
1763 183 : IF (gradient_f) THEN
1764 126 : 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 126 : drho_rad_s=dr_s, rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
1767 1581552 : drho0_h = 0.0_dp
1768 1581552 : drho0_s = 0.0_dp
1769 : ELSE
1770 57 : NULLIFY (r_h, r_s)
1771 57 : CALL get_rho_atom(rho_atom=rho0_atom, rho_rad_h=r_h, rho_rad_s=r_s)
1772 57 : rho_d = 0.0_dp
1773 : END IF
1774 9333 : 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 9150 : r_h_d, r_s_d, drho0_h, drho0_s)
1778 9333 : IF (donlcc) THEN
1779 : CALL calc_rho_nlcc(grid_atom, nspins, gradient_f, &
1780 0 : ir, rho_nlcc(:, 1), rho0_h, rho0_s, rho_nlcc(:, 2), drho0_h, drho0_s)
1781 : END IF
1782 : END DO
1783 183 : IF (tau_f) THEN
1784 : !compute tau on the grid all at once
1785 0 : CALL calc_tau_atom(tau0_h, tau0_s, rho0_atom, tau_basis_cache, nspins)
1786 : ELSE
1787 183 : tau_d = 0.0_dp
1788 : END IF
1789 : ! RHO1
1790 467016 : rho1_h = 0.0_dp
1791 467016 : rho1_s = 0.0_dp
1792 183 : rho1_atom => rho1_atom_set(iatom)
1793 183 : IF (gradient_f) THEN
1794 126 : 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 126 : drho_rad_s=dr_s, rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
1797 1581552 : drho1_h = 0.0_dp
1798 1581552 : drho1_s = 0.0_dp
1799 : ELSE
1800 57 : NULLIFY (r_h, r_s)
1801 57 : CALL get_rho_atom(rho_atom=rho1_atom, rho_rad_h=r_h, rho_rad_s=r_s)
1802 : END IF
1803 9333 : 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 9333 : r_h_d, r_s_d, drho1_h, drho1_s)
1807 : END DO
1808 183 : IF (tau_f) THEN
1809 : !compute tau on the grid all at once
1810 0 : CALL calc_tau_atom(tau1_h, tau1_s, rho1_atom, tau_basis_cache, nspins)
1811 : END IF
1812 :
1813 9333 : DO ir = 1, nr
1814 9333 : IF (tau_f) THEN
1815 0 : CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, drho1_h, tau1_h, na, ir)
1816 0 : CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, drho1_s, tau1_s, na, ir)
1817 9150 : ELSE IF (gradient_f) THEN
1818 6300 : CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, drho1_h, tau_d, na, ir)
1819 6300 : CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, drho1_s, tau_d, na, ir)
1820 : ELSE
1821 2850 : CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, rho_d, tau_d, na, ir)
1822 2850 : 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 183 : rho2_atom => rho2_atom_set(iatom)
1828 :
1829 1464 : DO istep = -nstep, nstep
1830 :
1831 1281 : beta = REAL(istep, KIND=dp)*epsrho
1832 :
1833 6536943 : rho_h = rho0_h + beta*rho1_h
1834 6536943 : rho_s = rho0_s + beta*rho1_s
1835 1281 : IF (gradient_f) THEN
1836 22140846 : drho_h = drho0_h + beta*drho1_h
1837 22140846 : drho_s = drho0_s + beta*drho1_s
1838 : END IF
1839 1281 : IF (tau_f) THEN
1840 0 : tau_h = tau0_h + beta*tau1_h
1841 0 : tau_s = tau0_s + beta*tau1_s
1842 : END IF
1843 : !
1844 1281 : IF (gradient_f) THEN
1845 : drho_h(4, :, :, :) = SQRT( &
1846 : drho_h(1, :, :, :)*drho_h(1, :, :, :) + &
1847 : drho_h(2, :, :, :)*drho_h(2, :, :, :) + &
1848 2250864 : 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 2250864 : drho_s(3, :, :, :)*drho_s(3, :, :, :))
1854 : END IF
1855 :
1856 65331 : DO ir = 1, nr
1857 65331 : IF (tau_f) THEN
1858 0 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_h, na, ir)
1859 0 : CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_s, na, ir)
1860 64050 : ELSE IF (gradient_f) THEN
1861 44100 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_d, na, ir)
1862 44100 : CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_d, na, ir)
1863 : ELSE
1864 19950 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, rho_d, tau_d, na, ir)
1865 19950 : 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 1281 : 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 1281 : do_triplet=is_triplet)
1876 : ! soft atom density !
1877 1281 : 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 1281 : do_triplet=is_triplet)
1883 : ! potentials
1884 2562 : DO ns = 1, nspins
1885 2410793 : fint_hh(ns)%r_coef(:, :) = 0.0_dp
1886 2412074 : fint_ss(ns)%r_coef(:, :) = 0.0_dp
1887 : END DO
1888 1281 : IF (gradient_f) THEN
1889 : CALL gaVxcgb_GC(vxc_h, vxc_s, vxg_h, vxg_s, fint_hh, fint_ss, &
1890 882 : grid_atom, basis_1c, harmonics, nspins)
1891 : ELSE
1892 : CALL gaVxcgb_noGC(vxc_h, vxc_s, fint_hh, fint_ss, &
1893 399 : grid_atom, basis_1c, harmonics, nspins)
1894 : END IF
1895 1281 : IF (tau_f) THEN
1896 : CALL dgaVtaudgb(vtau_h, vtau_s, fint_hh, fint_ss, &
1897 0 : tau_basis_cache, nspins)
1898 : END IF
1899 : ! second derivative gxc
1900 1281 : NULLIFY (int_hh, int_ss)
1901 1281 : CALL get_rho_atom(rho_atom=rho2_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss)
1902 2745 : DO ns = 1, nspins
1903 4820305 : int_ss(ns)%r_coef(:, :) = int_ss(ns)%r_coef(:, :) + oeps1*ak(istep)*fint_ss(ns)%r_coef(:, :)
1904 4821586 : 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 366 : DO ns = 1, nspins
1909 183 : DEALLOCATE (fint_ss(ns)%r_coef)
1910 366 : DEALLOCATE (fint_hh(ns)%r_coef)
1911 : END DO
1912 441 : DEALLOCATE (fint_ss, fint_hh)
1913 :
1914 : END DO ! iat
1915 :
1916 : ! Release the xc structure used to store the xc derivatives
1917 258 : CALL xc_dset_release(deriv_set)
1918 258 : CALL xc_rho_set_release(rho_set_h)
1919 258 : CALL xc_rho_set_release(rho_set_s)
1920 258 : CALL xc_rho_set_release(rho1_set_h)
1921 258 : CALL xc_rho_set_release(rho1_set_s)
1922 :
1923 258 : DEALLOCATE (rho_h, rho_s, rho0_h, rho0_s, rho1_h, rho1_s)
1924 258 : DEALLOCATE (vxc_h, vxc_s)
1925 258 : IF (gradient_f) THEN
1926 176 : DEALLOCATE (drho_h, drho_s, drho0_h, drho0_s, drho1_h, drho1_s)
1927 176 : DEALLOCATE (vxg_h, vxg_s)
1928 : END IF
1929 902 : IF (tau_f) THEN
1930 0 : DEALLOCATE (tau_h, tau_s, tau0_h, tau0_s, tau1_h, tau1_s)
1931 0 : DEALLOCATE (vtau_h, vtau_s)
1932 0 : CALL release_tau_basis_cache(tau_basis_cache)
1933 : END IF
1934 : END DO ! ikind
1935 :
1936 : END IF !xc_none
1937 :
1938 114 : CALL timestop(handle)
1939 :
1940 8436 : 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 2861040 : 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 2861040 : CPASSERT(ASSOCIATED(r_h))
1979 2861040 : CPASSERT(ASSOCIATED(r_s))
1980 2861040 : CPASSERT(ASSOCIATED(rho_h))
1981 2861040 : CPASSERT(ASSOCIATED(rho_s))
1982 2861040 : IF (grad_func) THEN
1983 1796440 : CPASSERT(ASSOCIATED(dr_h))
1984 1796440 : CPASSERT(ASSOCIATED(dr_s))
1985 1796440 : CPASSERT(ASSOCIATED(r_h_d))
1986 1796440 : CPASSERT(ASSOCIATED(r_s_d))
1987 1796440 : CPASSERT(ASSOCIATED(drho_h))
1988 1796440 : CPASSERT(ASSOCIATED(drho_s))
1989 : END IF
1990 :
1991 2861040 : na = grid_atom%ng_sphere
1992 2861040 : rad = grid_atom%rad(ir)
1993 2861040 : urad = grid_atom%oorad2l(ir, 1)
1994 6063820 : DO ispin = 1, nspins
1995 49154440 : DO iso = 1, harmonics%max_iso_not0
1996 2202365920 : DO ia = 1, na
1997 : rho_h(ia, ir, ispin) = rho_h(ia, ir, ispin) + &
1998 2156072520 : r_h(ispin)%r_coef(ir, iso)*harmonics%slm(ia, iso)
1999 : rho_s(ia, ir, ispin) = rho_s(ia, ir, ispin) + &
2000 2199163140 : 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 2861040 : IF (grad_func) THEN
2006 3783570 : DO ispin = 1, nspins
2007 29585300 : DO iso = 1, harmonics%max_iso_not0
2008 1411035320 : 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 1381450020 : 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 1381450020 : 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 1381450020 : 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 1381450020 : 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 1381450020 : 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 1409048190 : harmonics%slm(ia, iso)
2047 :
2048 : END DO ! ia
2049 : END DO ! iso
2050 103218550 : 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 99434980 : 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 101422110 : 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 2861040 : 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 1218 : 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 1218 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: a1, a2, gexp, r1, r2
2083 1218 : REAL(dp), DIMENSION(:, :), POINTER :: slm
2084 1218 : REAL(dp), DIMENSION(:, :, :), POINTER :: dslm_dxyz
2085 :
2086 1218 : NULLIFY (slm, dslm_dxyz)
2087 :
2088 1218 : 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 1218 : zet=tau_cache%zet)
2094 : CALL get_paw_basis_info(basis_1c, o2nindex=tau_cache%o2nindex, &
2095 : n2oindex=tau_cache%n2oindex, &
2096 1218 : nsatbas=tau_cache%nsatbas)
2097 :
2098 1218 : tau_cache%nr = grid_atom%nr
2099 1218 : tau_cache%na = grid_atom%ng_sphere
2100 1218 : slm => harmonics%slm
2101 1218 : dslm_dxyz => harmonics%dslm_dxyz
2102 :
2103 6090 : 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 9744 : r1(tau_cache%nr), r2(tau_cache%nr))
2106 283015104 : tau_cache%grad = 0.0_dp
2107 :
2108 5092 : DO iset = 1, tau_cache%nset
2109 16818 : DO ipgf = 1, tau_cache%npgf(iset)
2110 : starti = (iset - 1)*tau_cache%maxso + &
2111 11726 : (ipgf - 1)*nsoset(tau_cache%lmax(iset))
2112 : gexp(1:tau_cache%nr) = EXP(-tau_cache%zet(ipgf, iset)* &
2113 611026 : grid_atom%rad2(1:tau_cache%nr))
2114 51544 : DO iso = nsoset(tau_cache%lmin(iset) - 1) + 1, nsoset(tau_cache%lmax(iset))
2115 35944 : ip = tau_cache%o2nindex(starti + iso)
2116 35944 : IF (ip == 0) CYCLE
2117 35944 : l = indso(1, iso)
2118 :
2119 1869544 : 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 1869544 : grid_atom%rad2(1:tau_cache%nr)*r1(1:tau_cache%nr)
2122 :
2123 155502 : DO dir = 1, 3
2124 5578056 : a1(1:tau_cache%na) = dslm_dxyz(dir, 1:tau_cache%na, iso)
2125 5578056 : a2(1:tau_cache%na) = harmonics%a(dir, 1:tau_cache%na)*slm(1:tau_cache%na, iso)
2126 5644576 : DO ir = 1, tau_cache%nr
2127 288511032 : DO ia = 1, tau_cache%na
2128 282902400 : igrid = ia + (ir - 1)*tau_cache%na
2129 288403200 : 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 1218 : DEALLOCATE (a1, a2, gexp, r1, r2)
2138 :
2139 1218 : 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 2436 : SUBROUTINE release_tau_basis_cache(tau_cache)
2146 :
2147 : TYPE(tau_basis_cache_type), INTENT(INOUT) :: tau_cache
2148 :
2149 2436 : IF (ALLOCATED(tau_cache%grad)) DEALLOCATE (tau_cache%grad)
2150 2436 : IF (ASSOCIATED(tau_cache%n2oindex)) DEALLOCATE (tau_cache%n2oindex)
2151 2436 : IF (ASSOCIATED(tau_cache%o2nindex)) DEALLOCATE (tau_cache%o2nindex)
2152 2436 : NULLIFY (tau_cache%lmax, tau_cache%lmin, tau_cache%n2oindex, tau_cache%npgf, &
2153 2436 : tau_cache%zet, tau_cache%o2nindex)
2154 2436 : tau_cache%maxso = 0
2155 2436 : tau_cache%na = 0
2156 2436 : tau_cache%nr = 0
2157 2436 : tau_cache%nsatbas = 0
2158 2436 : tau_cache%nset = 0
2159 :
2160 2436 : 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 849 : 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 849 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: work
2185 :
2186 849 : CALL timeset(routineN, handle)
2187 :
2188 849 : CPASSERT(ALLOCATED(tau_cache%grad))
2189 :
2190 : !zeroing tau, assuming it is already allocated
2191 2346956 : tau_h = 0.0_dp
2192 2346956 : tau_s = 0.0_dp
2193 :
2194 849 : nr = tau_cache%nr
2195 849 : na = tau_cache%na
2196 849 : nbas = tau_cache%nsatbas
2197 849 : ngrid = na*nr
2198 3396 : ALLOCATE (work(ngrid, nbas))
2199 :
2200 1706 : DO ispin = 1, nspins
2201 4277 : DO dir = 1, 3
2202 : CALL dgemm('N', 'T', ngrid, nbas, nbas, 0.5_dp, tau_cache%grad(:, :, dir), &
2203 2571 : ngrid, rho_atom%cpc_h(ispin)%r_coef, nbas, 0.0_dp, work, ngrid)
2204 72576 : DO ibas = 1, nbas
2205 3627426 : DO ir = 1, nr
2206 185298555 : DO ia = 1, na
2207 181673700 : igrid = ia + (ir - 1)*na
2208 : tau_h(ia, ir, ispin) = tau_h(ia, ir, ispin) + &
2209 185228550 : 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 2571 : ngrid, rho_atom%cpc_s(ispin)%r_coef, nbas, 0.0_dp, work, ngrid)
2216 73433 : DO ibas = 1, nbas
2217 3627426 : DO ir = 1, nr
2218 185298555 : DO ia = 1, na
2219 181673700 : igrid = ia + (ir - 1)*na
2220 : tau_s(ia, ir, ispin) = tau_s(ia, ir, ispin) + &
2221 185228550 : 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 849 : DEALLOCATE (work)
2229 :
2230 849 : CALL timestop(handle)
2231 :
2232 849 : 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 8650 : SUBROUTINE calc_rho_nlcc(grid_atom, nspins, grad_func, &
2248 8650 : 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 8650 : CPASSERT(ASSOCIATED(rho_h))
2263 8650 : CPASSERT(ASSOCIATED(rho_s))
2264 8650 : IF (grad_func) THEN
2265 8650 : CPASSERT(ASSOCIATED(drho_h))
2266 8650 : CPASSERT(ASSOCIATED(drho_s))
2267 : END IF
2268 :
2269 8650 : na = grid_atom%ng_sphere
2270 8650 : rad = grid_atom%rad(ir)
2271 8650 : urad = grid_atom%oorad2l(ir, 1)
2272 :
2273 8650 : xsp = REAL(nspins, KIND=dp)
2274 8650 : rho = rho_nlcc(ir)/xsp
2275 17300 : DO ispin = 1, nspins
2276 441150 : rho_h(1:na, ir, ispin) = rho_h(1:na, ir, ispin) + rho
2277 449800 : rho_s(1:na, ir, ispin) = rho_s(1:na, ir, ispin) + rho
2278 : END DO ! ispin
2279 :
2280 8650 : IF (grad_func) THEN
2281 8650 : drho = drho_nlcc(ir)/xsp
2282 17300 : DO ispin = 1, nspins
2283 449800 : DO ia = 1, na
2284 432500 : IF (grid_atom%azi(ia) == 0.0_dp) THEN
2285 : dx = 0.0_dp
2286 : dy = 0.0_dp
2287 : ELSE
2288 389250 : dx = grid_atom%sin_pol(ia)*grid_atom%sin_azi(ia)
2289 389250 : dy = grid_atom%sin_pol(ia)*grid_atom%cos_azi(ia)
2290 : END IF
2291 432500 : dz = grid_atom%cos_pol(ia)
2292 : ! components of the gradient of rho1 hard
2293 432500 : drho_h(1, ia, ir, ispin) = drho_h(1, ia, ir, ispin) + drho*dx
2294 432500 : drho_h(2, ia, ir, ispin) = drho_h(2, ia, ir, ispin) + drho*dy
2295 432500 : drho_h(3, ia, ir, ispin) = drho_h(3, ia, ir, ispin) + drho*dz
2296 : ! components of the gradient of rho1 soft
2297 432500 : drho_s(1, ia, ir, ispin) = drho_s(1, ia, ir, ispin) + drho*dx
2298 432500 : drho_s(2, ia, ir, ispin) = drho_s(2, ia, ir, ispin) + drho*dy
2299 432500 : 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 432500 : 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 441150 : 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 8650 : 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 15921 : 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 15921 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
2344 15921 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: g1, g2
2345 15921 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: gg, gVg_h, gVg_s, matso_h, matso_s, vx
2346 15921 : REAL(dp), DIMENSION(:, :), POINTER :: zet
2347 : REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG
2348 :
2349 15921 : CALL timeset(routineN, handle)
2350 :
2351 15921 : 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 15921 : nset=nset, zet=zet)
2356 :
2357 15921 : nr = grid_atom%nr
2358 15921 : na = grid_atom%ng_sphere
2359 15921 : my_CG => harmonics%my_CG
2360 15921 : max_iso_not0 = harmonics%max_iso_not0
2361 15921 : lmax_expansion = indso(1, max_iso_not0)
2362 15921 : max_s_harm = harmonics%max_s_harm
2363 :
2364 111447 : ALLOCATE (g1(nr), g2(nr), gg(nr, 0:2*maxl))
2365 95526 : ALLOCATE (gVg_h(na, 0:2*maxl), gVg_s(na, 0:2*maxl))
2366 : ALLOCATE (matso_h(nsoset(maxl), nsoset(maxl)), &
2367 95526 : matso_s(nsoset(maxl), nsoset(maxl)))
2368 63684 : ALLOCATE (vx(na, nr))
2369 95526 : ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm))
2370 :
2371 958271 : g1 = 0.0_dp
2372 958271 : g2 = 0.0_dp
2373 15921 : m1 = 0
2374 54528 : DO iset1 = 1, nset
2375 38607 : n1 = nsoset(lmax(iset1))
2376 38607 : m2 = 0
2377 160646 : DO iset2 = 1, nset
2378 : CALL get_none0_cg_list(my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
2379 122039 : max_s_harm, lmax_expansion, cg_list, cg_n_list, max_iso_not0_local)
2380 122039 : CPASSERT(max_iso_not0_local <= max_iso_not0)
2381 :
2382 122039 : n2 = nsoset(lmax(iset2))
2383 447916 : DO ipgf1 = 1, npgf(iset1)
2384 325877 : ngau1 = n1*(ipgf1 - 1) + m1
2385 325877 : size1 = nsoset(lmax(iset1)) - nsoset(lmin(iset1) - 1)
2386 325877 : nngau1 = nsoset(lmin(iset1) - 1) + ngau1
2387 :
2388 18743727 : g1(1:nr) = EXP(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))
2389 1514144 : DO ipgf2 = 1, npgf(iset2)
2390 1066228 : ngau2 = n2*(ipgf2 - 1) + m2
2391 :
2392 60868328 : g2(1:nr) = EXP(-zet(ipgf2, iset2)*grid_atom%rad2(1:nr))
2393 1066228 : lmin12 = lmin(iset1) + lmin(iset2)
2394 1066228 : lmax12 = lmax(iset1) + lmax(iset2)
2395 :
2396 : ! reduce expansion local densities
2397 1392105 : IF (lmin12 <= lmax_expansion) THEN
2398 :
2399 262621710 : gg = 0.0_dp
2400 1065283 : IF (lmin12 == 0) THEN
2401 33750591 : gg(1:nr, lmin12) = g1(1:nr)*g2(1:nr)
2402 : ELSE
2403 27069542 : 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 1065283 : IF (lmax12 > lmax_expansion) lmax12 = lmax_expansion
2408 :
2409 1562895 : DO l = lmin12 + 1, lmax12
2410 30944295 : gg(1:nr, l) = grid_atom%rad(1:nr)*gg(:, l - 1)
2411 : END DO
2412 :
2413 2300622 : DO ispin = 1, nspins
2414 1235339 : ld = lmax12 + 1
2415 74505289 : DO ir = 1, nr
2416 3738002789 : 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 1235339 : gg(1:nr, 0:lmax12), nr, 0.0_dp, gVg_h(1:na, 0:lmax12), na)
2420 74505289 : DO ir = 1, nr
2421 3738002789 : 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 1235339 : gg(1:nr, 0:lmax12), nr, 0.0_dp, gVg_s(1:na, 0:lmax12), na)
2425 :
2426 103955651 : matso_h = 0.0_dp
2427 103955651 : matso_s = 0.0_dp
2428 9730568 : DO iso = 1, max_iso_not0_local
2429 26894967 : DO icg = 1, cg_n_list(iso)
2430 17164399 : iso1 = cg_list(1, icg, iso)
2431 17164399 : iso2 = cg_list(2, icg, iso)
2432 17164399 : l = indso(1, iso1) + indso(1, iso2)
2433 :
2434 17164399 : CPASSERT(l <= lmax_expansion)
2435 883879578 : 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 858219950 : 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 875384349 : harmonics%slm(ia, iso)
2444 : END DO
2445 : END DO
2446 : END DO
2447 :
2448 : ! Write in the global matrix
2449 5550889 : DO ic = nsoset(lmin(iset2) - 1) + 1, nsoset(lmax(iset2))
2450 3250267 : iso1 = nsoset(lmin(iset1) - 1) + 1
2451 3250267 : iso2 = ngau2 + ic
2452 : CALL daxpy(size1, 1.0_dp, matso_h(iso1, ic), 1, &
2453 3250267 : int_hh(ispin)%r_coef(nngau1 + 1, iso2), 1)
2454 : CALL daxpy(size1, 1.0_dp, matso_s(iso1, ic), 1, &
2455 4485606 : 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 282685 : m2 = m2 + maxso
2465 : END DO ! iset2
2466 54528 : m1 = m1 + maxso
2467 : END DO ! iset1
2468 :
2469 15921 : DEALLOCATE (g1, g2, gg, matso_h, matso_s, gVg_s, gVg_h, vx)
2470 :
2471 15921 : DEALLOCATE (cg_list, cg_n_list)
2472 :
2473 15921 : CALL timestop(handle)
2474 :
2475 15921 : 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 28702 : 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 28702 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
2510 : REAL(dp) :: urad
2511 28702 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: g1, g2
2512 28702 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: dgg, gg, gVXCg_h, gVXCg_s, matso_h, &
2513 28702 : matso_s
2514 28702 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: gVXGg_h, gVXGg_s
2515 28702 : REAL(dp), DIMENSION(:, :), POINTER :: zet
2516 : REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG
2517 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: my_CG_dxyz
2518 :
2519 28702 : CALL timeset(routineN, handle)
2520 :
2521 28702 : 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 28702 : nset=nset, zet=zet)
2526 :
2527 28702 : nr = grid_atom%nr
2528 28702 : na = grid_atom%ng_sphere
2529 28702 : my_CG => harmonics%my_CG
2530 28702 : my_CG_dxyz => harmonics%my_CG_dxyz
2531 28702 : max_iso_not0 = harmonics%max_iso_not0
2532 28702 : lmax_expansion = indso(1, max_iso_not0)
2533 28702 : max_s_harm = harmonics%max_s_harm
2534 :
2535 258318 : ALLOCATE (g1(nr), g2(nr), gg(nr, 0:2*maxl), dgg(nr, 0:2*maxl))
2536 172212 : ALLOCATE (gVXCg_h(na, 0:2*maxl), gVXCg_s(na, 0:2*maxl))
2537 172212 : 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 258318 : 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 172212 : matso_s(nsoset(maxl), nsoset(maxl)))
2543 :
2544 60426 : DO ispin = 1, nspins
2545 :
2546 1643004 : g1 = 0.0_dp
2547 1643004 : g2 = 0.0_dp
2548 31724 : m1 = 0
2549 141607 : DO iset1 = 1, nset
2550 81181 : n1 = nsoset(lmax(iset1))
2551 81181 : m2 = 0
2552 354926 : DO iset2 = 1, nset
2553 : CALL get_none0_cg_list(my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
2554 273745 : max_s_harm, lmax_expansion, cg_list, cg_n_list, max_iso_not0_local)
2555 273745 : 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 273745 : max_s_harm, lmax_expansion, dcg_list, dcg_n_list, dmax_iso_not0_local)
2558 :
2559 273745 : n2 = nsoset(lmax(iset2))
2560 888315 : DO ipgf1 = 1, npgf(iset1)
2561 614570 : ngau1 = n1*(ipgf1 - 1) + m1
2562 614570 : size1 = nsoset(lmax(iset1)) - nsoset(lmin(iset1) - 1)
2563 614570 : nngau1 = nsoset(lmin(iset1) - 1) + ngau1
2564 :
2565 32064150 : g1(1:nr) = EXP(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))
2566 2432713 : DO ipgf2 = 1, npgf(iset2)
2567 1544398 : ngau2 = n2*(ipgf2 - 1) + m2
2568 :
2569 80790978 : g2(1:nr) = EXP(-zet(ipgf2, iset2)*grid_atom%rad2(1:nr))
2570 1544398 : lmin12 = lmin(iset1) + lmin(iset2)
2571 1544398 : lmax12 = lmax(iset1) + lmax(iset2)
2572 :
2573 : !test reduce expansion local densities
2574 1544398 : IF (lmin12 <= lmax_expansion) THEN
2575 :
2576 338286432 : gg = 0.0_dp
2577 338286432 : dgg = 0.0_dp
2578 :
2579 1543798 : IF (lmin12 == 0) THEN
2580 48750516 : gg(1:nr, lmin12) = g1(1:nr)*g2(1:nr)
2581 : ELSE
2582 32009862 : 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 1543798 : IF (lmax12 > lmax_expansion) lmax12 = lmax_expansion
2587 :
2588 2375378 : DO l = lmin12 + 1, lmax12
2589 43869380 : 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 45413178 : 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 80760378 : gg(1:nr, lmax12)
2596 :
2597 328522708 : gVXCg_h = 0.0_dp
2598 328522708 : gVXCg_s = 0.0_dp
2599 1290216616 : gVXGg_h = 0.0_dp
2600 1290216616 : gVXGg_s = 0.0_dp
2601 :
2602 : ! Cross Term
2603 3919176 : DO l = lmin12, lmax12
2604 122639308 : DO ia = 1, na
2605 6237356910 : 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 6116261400 : 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 6116261400 : vxg_s(3, ia, ir, ispin)*harmonics%a(3, ia))
2619 :
2620 6116261400 : 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 6116261400 : gg(ir, l)*urad
2625 :
2626 : gVXGg_h(2, ia, l) = gVXGg_h(2, ia, l) + &
2627 : vxg_h(2, ia, ir, ispin)* &
2628 6116261400 : gg(ir, l)*urad
2629 :
2630 : gVXGg_h(3, ia, l) = gVXGg_h(3, ia, l) + &
2631 : vxg_h(3, ia, ir, ispin)* &
2632 6116261400 : gg(ir, l)*urad
2633 :
2634 : gVXGg_s(1, ia, l) = gVXGg_s(1, ia, l) + &
2635 : vxg_s(1, ia, ir, ispin)* &
2636 6116261400 : gg(ir, l)*urad
2637 :
2638 : gVXGg_s(2, ia, l) = gVXGg_s(2, ia, l) + &
2639 : vxg_s(2, ia, ir, ispin)* &
2640 6116261400 : gg(ir, l)*urad
2641 :
2642 : gVXGg_s(3, ia, l) = gVXGg_s(3, ia, l) + &
2643 : vxg_s(3, ia, ir, ispin)* &
2644 6234981532 : gg(ir, l)*urad
2645 :
2646 : END DO ! ir
2647 : END DO ! ia
2648 : END DO ! l
2649 :
2650 112813382 : matso_h = 0.0_dp
2651 112813382 : matso_s = 0.0_dp
2652 11129206 : DO iso = 1, max_iso_not0_local
2653 31436911 : DO icg = 1, cg_n_list(iso)
2654 20307705 : iso1 = cg_list(1, icg, iso)
2655 20307705 : iso2 = cg_list(2, icg, iso)
2656 :
2657 20307705 : l = indso(1, iso1) + indso(1, iso2)
2658 :
2659 : !test reduce expansion local densities
2660 20307705 : CPASSERT(l <= lmax_expansion)
2661 1045030023 : DO ia = 1, na
2662 : matso_h(iso1, iso2) = matso_h(iso1, iso2) + &
2663 : gVXCg_h(ia, l)* &
2664 : harmonics%slm(ia, iso)* &
2665 1015136910 : 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 1035444615 : 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 5936116 : DO iso = 1, dmax_iso_not0_local
2679 39057859 : DO icg = 1, dcg_n_list(iso)
2680 33121743 : iso1 = dcg_list(1, icg, iso)
2681 33121743 : iso2 = dcg_list(2, icg, iso)
2682 :
2683 33121743 : l = indso(1, iso1) + indso(1, iso2)
2684 : !test reduce expansion local densities
2685 33121743 : CPASSERT(l <= lmax_expansion)
2686 1692861879 : 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 1655347818 : 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 1688469561 : 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 6090568 : DO ic = nsoset(lmin(iset2) - 1) + 1, nsoset(lmax(iset2))
2710 3931600 : iso1 = nsoset(lmin(iset1) - 1) + 1
2711 3931600 : iso2 = ngau2 + ic
2712 : CALL daxpy(size1, 1.0_dp, matso_h(iso1, ic), 1, &
2713 3931600 : int_hh(ispin)%r_coef(nngau1 + 1, iso2), 1)
2714 : CALL daxpy(size1, 1.0_dp, matso_s(iso1, ic), 1, &
2715 5475998 : int_ss(ispin)%r_coef(nngau1 + 1, iso2), 1)
2716 : END DO
2717 :
2718 : END DO ! ipfg2
2719 : END DO ! ipfg1
2720 902416 : m2 = m2 + maxso
2721 : END DO ! iset2
2722 112905 : m1 = m1 + maxso
2723 : END DO ! iset1
2724 : END DO ! ispin
2725 :
2726 28702 : DEALLOCATE (g1, g2, gg, dgg, matso_h, matso_s, gVXCg_h, gVXCg_s, gVXGg_h, gVXGg_s)
2727 28702 : DEALLOCATE (cg_list, cg_n_list, dcg_list, dcg_n_list)
2728 :
2729 28702 : CALL timestop(handle)
2730 :
2731 28702 : 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 849 : 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 849 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: int_h, int_s, weighted_grad
2758 :
2759 849 : CALL timeset(routineN, handle)
2760 :
2761 849 : CPASSERT(ALLOCATED(tau_cache%grad))
2762 849 : CPASSERT(ASSOCIATED(tau_cache%n2oindex))
2763 :
2764 849 : nr = tau_cache%nr
2765 849 : na = tau_cache%na
2766 849 : nbas = tau_cache%nsatbas
2767 849 : ngrid = na*nr
2768 66894 : max_old_basis = MAXVAL(tau_cache%n2oindex)
2769 7641 : ALLOCATE (int_h(nbas, nbas), int_s(nbas, nbas), weighted_grad(ngrid, nbas))
2770 :
2771 1706 : DO ispin = 1, nspins
2772 857 : CPASSERT(SIZE(int_hh(ispin)%r_coef, 1) >= max_old_basis)
2773 857 : CPASSERT(SIZE(int_hh(ispin)%r_coef, 2) >= max_old_basis)
2774 857 : CPASSERT(SIZE(int_ss(ispin)%r_coef, 1) >= max_old_basis)
2775 857 : CPASSERT(SIZE(int_ss(ispin)%r_coef, 2) >= max_old_basis)
2776 951505 : int_h = 0.0_dp
2777 951505 : int_s = 0.0_dp
2778 3428 : DO dir = 1, 3
2779 72576 : DO ibas = 1, nbas
2780 3627426 : DO ir = 1, nr
2781 185298555 : DO ia = 1, na
2782 181673700 : igrid = ia + (ir - 1)*na
2783 : weighted_grad(igrid, ibas) = vtau_h(ia, ir, ispin)* &
2784 185228550 : 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 2571 : ngrid, weighted_grad, ngrid, 1.0_dp, int_h, nbas)
2790 :
2791 72576 : DO ibas = 1, nbas
2792 3627426 : DO ir = 1, nr
2793 185298555 : DO ia = 1, na
2794 181673700 : igrid = ia + (ir - 1)*na
2795 : weighted_grad(igrid, ibas) = vtau_s(ia, ir, ispin)* &
2796 185228550 : 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 3428 : ngrid, weighted_grad, ngrid, 1.0_dp, int_s, nbas)
2802 : END DO
2803 :
2804 25041 : DO jbas = 1, nbas
2805 23335 : jold = tau_cache%n2oindex(jbas)
2806 951505 : DO ibas = 1, nbas
2807 927313 : iold = tau_cache%n2oindex(ibas)
2808 : int_hh(ispin)%r_coef(iold, jold) = int_hh(ispin)%r_coef(iold, jold) + &
2809 927313 : int_h(ibas, jbas)
2810 : int_ss(ispin)%r_coef(iold, jold) = int_ss(ispin)%r_coef(iold, jold) + &
2811 950648 : int_s(ibas, jbas)
2812 : END DO
2813 : END DO
2814 : END DO
2815 :
2816 849 : DEALLOCATE (int_h, int_s, weighted_grad)
2817 :
2818 849 : CALL timestop(handle)
2819 :
2820 849 : END SUBROUTINE dgaVtaudgb
2821 :
2822 0 : END MODULE qs_vxc_atom
|
