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