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 Fxc kernel calculated
10 : !> for the atomic density in the basis set of spherical primitives
11 : ! **************************************************************************************************
12 : MODULE qs_fgxc_atom
13 : USE atomic_kind_types, ONLY: atomic_kind_type,&
14 : get_atomic_kind
15 : USE basis_set_types, ONLY: gto_basis_set_type
16 : USE cp_control_types, ONLY: dft_control_type
17 : USE input_constants, ONLY: xc_none
18 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
19 : section_vals_type,&
20 : section_vals_val_get
21 : USE kinds, ONLY: dp
22 : USE message_passing, ONLY: mp_para_env_type
23 : USE qs_environment_types, ONLY: get_qs_env,&
24 : qs_environment_type
25 : USE qs_fxc_atom, ONLY: fxc_atom_calc
26 : USE qs_grid_atom, ONLY: grid_atom_type
27 : USE qs_harmonics_atom, ONLY: harmonics_atom_type
28 : USE qs_kind_types, ONLY: get_qs_kind,&
29 : has_nlcc,&
30 : qs_kind_type
31 : USE qs_rho_atom_types, ONLY: get_rho_atom,&
32 : rho_atom_coeff,&
33 : rho_atom_type
34 : USE qs_vxc_atom_utils, ONLY: &
35 : calc_rho_angular, calc_rho_nlcc, calc_tau_atom, create_tau_basis_cache, dgaVtaudgb, &
36 : gaVxcgb_GC, gaVxcgb_noGC, release_tau_basis_cache, tau_basis_cache_type
37 : USE util, ONLY: get_limit
38 : USE xc_atom, ONLY: fill_rho_set,&
39 : vxc_of_r_new,&
40 : xc_2nd_deriv_of_r,&
41 : xc_rho_set_atom_update
42 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
43 : xc_dset_create,&
44 : xc_dset_release,&
45 : xc_dset_zero_all
46 : USE xc_derivatives, ONLY: xc_functionals_get_needs
47 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
48 : USE xc_rho_set_types, ONLY: xc_rho_set_create,&
49 : xc_rho_set_release,&
50 : xc_rho_set_type
51 : #include "./base/base_uses.f90"
52 :
53 : IMPLICIT NONE
54 :
55 : PRIVATE
56 :
57 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_fgxc_atom'
58 :
59 : PUBLIC :: fgxc_atom_calc, &
60 : fgxc_atom_diff
61 :
62 : ! **************************************************************************************************
63 :
64 : CONTAINS
65 :
66 : ! **************************************************************************************************
67 : !> \brief ...
68 : !> \param qs_env ...
69 : !> \param rho0_atom_set ...
70 : !> \param rho1_atom_set ...
71 : !> \param rho2_atom_set ...
72 : !> \param kind_set ...
73 : !> \param xc_section ...
74 : !> \param is_triplet ...
75 : !> \param accuracy ...
76 : ! **************************************************************************************************
77 0 : SUBROUTINE fgxc_atom_calc(qs_env, rho0_atom_set, rho1_atom_set, rho2_atom_set, &
78 : kind_set, xc_section, is_triplet, accuracy)
79 :
80 : TYPE(qs_environment_type), POINTER :: qs_env
81 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho0_atom_set, rho1_atom_set, &
82 : rho2_atom_set
83 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: kind_set
84 : TYPE(section_vals_type), OPTIONAL, POINTER :: xc_section
85 : LOGICAL, INTENT(IN) :: is_triplet
86 : INTEGER, INTENT(IN) :: accuracy
87 :
88 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fgxc_atom_calc'
89 : REAL(KIND=dp), PARAMETER :: epsrho = 5.e-4_dp
90 :
91 : INTEGER :: bo(2), handle, iat, iatom, ikind, ir, &
92 : istep, mspins, myfun, na, natom, nf, &
93 : nr, ns, nspins, nstep, num_pe
94 : INTEGER, DIMENSION(2, 3) :: bounds
95 0 : INTEGER, DIMENSION(:), POINTER :: atom_list
96 : LOGICAL :: accint, donlcc, gradient_f, lsd, nlcc, &
97 : paw_atom, tau_f
98 : REAL(dp) :: agr, alpha, beta, density_cut, exc_h, &
99 : exc_s, gradient_cut, oeps1, oeps2, &
100 : tau_cut
101 : REAL(dp), DIMENSION(1, 1, 1) :: tau_d
102 : REAL(dp), DIMENSION(1, 1, 1, 1) :: rho_d
103 0 : REAL(dp), DIMENSION(:, :), POINTER :: rho_nlcc, weight_h, weight_s
104 0 : REAL(dp), DIMENSION(:, :, :), POINTER :: rho0_h, rho0_s, rho1_h, rho1_s, rho_h, &
105 0 : rho_s, tau0_h, tau0_s, tau1_h, tau1_s, &
106 0 : tau_h, tau_s, vtau_h, vtau_s, vxc_h, &
107 0 : vxc_s
108 0 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: drho0_h, drho0_s, drho1_h, drho1_s, &
109 0 : drho_h, drho_s, vxg_h, vxg_s
110 : REAL(KIND=dp), DIMENSION(-4:4) :: ak, bl
111 0 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
112 : TYPE(dft_control_type), POINTER :: dft_control
113 : TYPE(grid_atom_type), POINTER :: grid_atom
114 : TYPE(gto_basis_set_type), POINTER :: basis_1c
115 : TYPE(harmonics_atom_type), POINTER :: harmonics
116 : TYPE(mp_para_env_type), POINTER :: para_env
117 0 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s, fint_hh, fint_ss, int_hh, &
118 0 : int_ss, r_h, r_s
119 0 : TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r_h_d, r_s_d
120 : TYPE(rho_atom_type), POINTER :: rho0_atom, rho1_atom, rho2_atom
121 : TYPE(section_vals_type), POINTER :: xc_fun_section
122 0 : TYPE(tau_basis_cache_type) :: tau_basis_cache
123 : TYPE(xc_derivative_set_type) :: deriv_set
124 : TYPE(xc_rho_cflags_type) :: needs
125 : TYPE(xc_rho_set_type) :: rho_set_h, rho_set_s
126 :
127 0 : CALL timeset(routineN, handle)
128 :
129 0 : NULLIFY (vtau_h, vtau_s)
130 :
131 0 : ak = 0.0_dp
132 0 : bl = 0.0_dp
133 0 : SELECT CASE (accuracy)
134 : CASE (:4)
135 0 : nstep = 2
136 0 : ak(-2:2) = [1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp]/12.0_dp
137 0 : bl(-2:2) = [-1.0_dp, 16.0_dp, -30.0_dp, 16.0_dp, -1.0_dp]/12.0_dp
138 : CASE (5:7)
139 0 : nstep = 3
140 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
141 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
142 : CASE (8:)
143 0 : nstep = 4
144 : ak(-4:4) = [1.0_dp, -32.0_dp/3.0_dp, 56.0_dp, -224.0_dp, 0.0_dp, &
145 0 : 224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp]/280.0_dp
146 : bl(-4:4) = [-1.0_dp, 128.0_dp/9.0_dp, -112.0_dp, 896.0_dp, -14350.0_dp/9.0_dp, &
147 0 : 896.0_dp, -112.0_dp, 128.0_dp/9.0_dp, -1.0_dp]/560.0_dp
148 : END SELECT
149 0 : oeps1 = 1.0_dp/epsrho
150 0 : oeps2 = 1.0_dp/(epsrho**2)
151 :
152 : CALL get_qs_env(qs_env=qs_env, &
153 : dft_control=dft_control, &
154 : para_env=para_env, &
155 0 : atomic_kind_set=atomic_kind_set)
156 :
157 0 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
158 0 : CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
159 :
160 0 : accint = dft_control%qs_control%gapw_control%accurate_xcint
161 :
162 0 : IF (myfun == xc_none) THEN
163 : ! no action needed?
164 : ELSE
165 0 : CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", r_val=density_cut)
166 0 : CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", r_val=gradient_cut)
167 0 : CALL section_vals_val_get(xc_section, "TAU_CUTOFF", r_val=tau_cut)
168 :
169 0 : nlcc = has_nlcc(kind_set)
170 0 : lsd = dft_control%lsd
171 0 : nspins = dft_control%nspins
172 0 : mspins = nspins
173 0 : IF (is_triplet) THEN
174 0 : CPASSERT(nspins == 1)
175 0 : lsd = .TRUE.
176 0 : mspins = 2
177 : END IF
178 0 : needs = xc_functionals_get_needs(xc_fun_section, lsd=lsd, calc_potential=.TRUE.)
179 0 : gradient_f = (needs%drho .OR. needs%drho_spin)
180 0 : tau_f = (needs%tau .OR. needs%tau_spin)
181 :
182 : ! Here starts the loop over all the atoms
183 0 : DO ikind = 1, SIZE(atomic_kind_set)
184 0 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
185 : CALL get_qs_kind(kind_set(ikind), paw_atom=paw_atom, &
186 0 : harmonics=harmonics, grid_atom=grid_atom)
187 0 : CALL get_qs_kind(kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
188 :
189 0 : IF (.NOT. paw_atom) CYCLE
190 :
191 0 : nr = grid_atom%nr
192 0 : na = grid_atom%ng_sphere
193 :
194 : ! set integration weights
195 0 : IF (accint) THEN
196 0 : weight_h => grid_atom%weight
197 0 : alpha = dft_control%qs_control%gapw_control%aw(ikind)
198 0 : IF (ASSOCIATED(grid_atom%gapw_weight_s)) THEN
199 0 : IF (grid_atom%gapw_weight_alpha /= alpha) DEALLOCATE (grid_atom%gapw_weight_s)
200 : END IF
201 0 : IF (.NOT. ASSOCIATED(grid_atom%gapw_weight_s)) THEN
202 0 : ALLOCATE (grid_atom%gapw_weight_s(na, nr))
203 0 : DO ir = 1, nr
204 0 : agr = 1.0_dp - EXP(-alpha*grid_atom%rad2(ir))
205 0 : grid_atom%gapw_weight_s(:, ir) = grid_atom%weight(:, ir)*agr
206 : END DO
207 0 : grid_atom%gapw_weight_alpha = alpha
208 : END IF
209 0 : weight_s => grid_atom%gapw_weight_s
210 : ELSE
211 0 : weight_h => grid_atom%weight
212 0 : weight_s => grid_atom%weight
213 : END IF
214 :
215 : ! Prepare the structures needed to calculate and store the xc derivatives
216 :
217 : ! Array dimension: here only one dimensional arrays are used,
218 : ! i.e. only the first column of deriv_data is read.
219 : ! The other to dimensions are set to size equal 1
220 0 : bounds(1:2, 1:3) = 1
221 0 : bounds(2, 1) = na
222 0 : bounds(2, 2) = nr
223 :
224 : ! create a place where to put the derivatives
225 0 : CALL xc_dset_create(deriv_set, local_bounds=bounds)
226 : ! create the place where to store the argument for the functionals
227 : CALL xc_rho_set_create(rho_set_h, bounds, rho_cutoff=density_cut, &
228 0 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
229 : CALL xc_rho_set_create(rho_set_s, bounds, rho_cutoff=density_cut, &
230 0 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
231 :
232 : ! allocate the required 3d arrays where to store rho and drho
233 0 : CALL xc_rho_set_atom_update(rho_set_h, needs, mspins, bounds)
234 0 : CALL xc_rho_set_atom_update(rho_set_s, needs, mspins, bounds)
235 :
236 : ALLOCATE (rho_h(na, nr, mspins), rho_s(na, nr, mspins), &
237 : rho0_h(na, nr, nspins), rho0_s(na, nr, nspins), &
238 0 : rho1_h(na, nr, nspins), rho1_s(na, nr, nspins))
239 0 : ALLOCATE (vxc_h(na, nr, mspins), vxc_s(na, nr, mspins))
240 0 : IF (gradient_f) THEN
241 : ALLOCATE (drho_h(4, na, nr, mspins), drho_s(4, na, nr, mspins), &
242 : drho0_h(4, na, nr, nspins), drho0_s(4, na, nr, nspins), &
243 0 : drho1_h(4, na, nr, nspins), drho1_s(4, na, nr, nspins))
244 0 : ALLOCATE (vxg_h(3, na, nr, mspins), vxg_s(3, na, nr, mspins))
245 : END IF
246 0 : IF (tau_f) THEN
247 0 : CALL create_tau_basis_cache(tau_basis_cache, grid_atom, basis_1c, harmonics)
248 : ALLOCATE (tau_h(na, nr, mspins), tau_s(na, nr, mspins), &
249 : tau0_h(na, nr, nspins), tau0_s(na, nr, nspins), &
250 0 : tau1_h(na, nr, nspins), tau1_s(na, nr, nspins))
251 0 : ALLOCATE (vtau_h(na, nr, mspins), vtau_s(na, nr, mspins))
252 : END IF
253 : !
254 : ! NLCC: prepare rho and drho of the core charge for this KIND
255 0 : donlcc = .FALSE.
256 0 : IF (nlcc) THEN
257 0 : NULLIFY (rho_nlcc)
258 0 : rho_nlcc => kind_set(ikind)%nlcc_pot
259 0 : IF (ASSOCIATED(rho_nlcc)) donlcc = .TRUE.
260 : END IF
261 :
262 : ! Distribute the atoms of this kind
263 0 : num_pe = para_env%num_pe
264 0 : bo = get_limit(natom, num_pe, para_env%mepos)
265 :
266 0 : DO iat = bo(1), bo(2)
267 0 : iatom = atom_list(iat)
268 : !
269 0 : NULLIFY (int_hh, int_ss)
270 0 : rho0_atom => rho0_atom_set(iatom)
271 0 : CALL get_rho_atom(rho_atom=rho0_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss)
272 0 : ALLOCATE (fint_ss(nspins), fint_hh(nspins))
273 0 : DO ns = 1, nspins
274 0 : nf = SIZE(int_ss(ns)%r_coef, 1)
275 0 : ALLOCATE (fint_ss(ns)%r_coef(nf, nf))
276 0 : nf = SIZE(int_hh(ns)%r_coef, 1)
277 0 : ALLOCATE (fint_hh(ns)%r_coef(nf, nf))
278 : END DO
279 :
280 : ! RHO0
281 0 : rho0_h = 0.0_dp
282 0 : rho0_s = 0.0_dp
283 0 : rho0_atom => rho0_atom_set(iatom)
284 0 : IF (gradient_f) THEN
285 0 : NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
286 : CALL get_rho_atom(rho_atom=rho0_atom, rho_rad_h=r_h, rho_rad_s=r_s, drho_rad_h=dr_h, &
287 0 : drho_rad_s=dr_s, rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
288 0 : drho0_h = 0.0_dp
289 0 : drho0_s = 0.0_dp
290 : ELSE
291 0 : NULLIFY (r_h, r_s)
292 0 : CALL get_rho_atom(rho_atom=rho0_atom, rho_rad_h=r_h, rho_rad_s=r_s)
293 0 : rho_d = 0.0_dp
294 : END IF
295 0 : DO ir = 1, nr
296 : CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_f, &
297 : ir, r_h, r_s, rho0_h, rho0_s, dr_h, dr_s, &
298 0 : r_h_d, r_s_d, drho0_h, drho0_s)
299 0 : IF (donlcc) THEN
300 : CALL calc_rho_nlcc(grid_atom, nspins, gradient_f, &
301 0 : ir, rho_nlcc(:, 1), rho0_h, rho0_s, rho_nlcc(:, 2), drho0_h, drho0_s)
302 : END IF
303 : END DO
304 0 : IF (tau_f) THEN
305 : !compute tau on the grid all at once
306 0 : CALL calc_tau_atom(tau0_h, tau0_s, rho0_atom, tau_basis_cache, nspins)
307 : ELSE
308 0 : tau_d = 0.0_dp
309 : END IF
310 : ! RHO1
311 0 : rho1_h = 0.0_dp
312 0 : rho1_s = 0.0_dp
313 0 : rho1_atom => rho1_atom_set(iatom)
314 0 : IF (gradient_f) THEN
315 0 : NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
316 : CALL get_rho_atom(rho_atom=rho1_atom, rho_rad_h=r_h, rho_rad_s=r_s, drho_rad_h=dr_h, &
317 0 : drho_rad_s=dr_s, rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
318 0 : drho1_h = 0.0_dp
319 0 : drho1_s = 0.0_dp
320 : ELSE
321 0 : NULLIFY (r_h, r_s)
322 0 : CALL get_rho_atom(rho_atom=rho1_atom, rho_rad_h=r_h, rho_rad_s=r_s)
323 : END IF
324 0 : DO ir = 1, nr
325 : CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_f, &
326 : ir, r_h, r_s, rho1_h, rho1_s, dr_h, dr_s, &
327 0 : r_h_d, r_s_d, drho1_h, drho1_s)
328 : END DO
329 0 : IF (tau_f) THEN
330 : !compute tau on the grid all at once
331 0 : CALL calc_tau_atom(tau1_h, tau1_s, rho1_atom, tau_basis_cache, nspins)
332 : END IF
333 : ! RHO2
334 0 : rho2_atom => rho2_atom_set(iatom)
335 :
336 0 : DO istep = -nstep, nstep
337 :
338 0 : beta = REAL(istep, KIND=dp)*epsrho
339 :
340 0 : IF (is_triplet) THEN
341 0 : rho_h(:, :, 1) = rho0_h(:, :, 1) + beta*rho1_h(:, :, 1)
342 0 : rho_h(:, :, 2) = rho0_h(:, :, 1)
343 0 : rho_h = 0.5_dp*rho_h
344 0 : rho_s(:, :, 1) = rho0_s(:, :, 1) + beta*rho1_s(:, :, 1)
345 0 : rho_s(:, :, 2) = rho0_s(:, :, 1)
346 0 : rho_s = 0.5_dp*rho_s
347 0 : IF (gradient_f) THEN
348 0 : drho_h(:, :, :, 1) = drho0_h(:, :, :, 1) + beta*drho1_h(:, :, :, 1)
349 0 : drho_h(:, :, :, 2) = drho0_h(:, :, :, 1)
350 0 : drho_h = 0.5_dp*drho_h
351 0 : drho_s(:, :, :, 1) = drho0_s(:, :, :, 1) + beta*drho1_s(:, :, :, 1)
352 0 : drho_s(:, :, :, 2) = drho0_s(:, :, :, 1)
353 0 : drho_s = 0.5_dp*drho_s
354 : END IF
355 0 : IF (tau_f) THEN
356 0 : tau_h(:, :, 1) = tau0_h(:, :, 1) + beta*tau1_h(:, :, 1)
357 0 : tau_h(:, :, 2) = tau0_h(:, :, 1)
358 0 : tau_h = 0.5_dp*tau0_h
359 0 : tau_s(:, :, 1) = tau0_s(:, :, 1) + beta*tau1_s(:, :, 1)
360 0 : tau_s(:, :, 2) = tau0_s(:, :, 1)
361 0 : tau_s = 0.5_dp*tau0_s
362 : END IF
363 : ELSE
364 0 : rho_h = rho0_h + beta*rho1_h
365 0 : rho_s = rho0_s + beta*rho1_s
366 0 : IF (gradient_f) THEN
367 0 : drho_h = drho0_h + beta*drho1_h
368 0 : drho_s = drho0_s + beta*drho1_s
369 : END IF
370 0 : IF (tau_f) THEN
371 0 : tau_h = tau0_h + beta*tau1_h
372 0 : tau_s = tau0_s + beta*tau1_s
373 : END IF
374 : END IF
375 : !
376 0 : IF (gradient_f) THEN
377 0 : drho_h(4, :, :, :) = NORM2(drho_h(1:3, :, :, :), 1)
378 0 : drho_s(4, :, :, :) = NORM2(drho_s(1:3, :, :, :), 1)
379 : END IF
380 :
381 0 : DO ir = 1, nr
382 0 : IF (tau_f) THEN
383 0 : CALL fill_rho_set(rho_set_h, lsd, mspins, needs, rho_h, drho_h, tau_h, na, ir)
384 0 : CALL fill_rho_set(rho_set_s, lsd, mspins, needs, rho_s, drho_s, tau_s, na, ir)
385 0 : ELSE IF (gradient_f) THEN
386 0 : CALL fill_rho_set(rho_set_h, lsd, mspins, needs, rho_h, drho_h, tau_d, na, ir)
387 0 : CALL fill_rho_set(rho_set_s, lsd, mspins, needs, rho_s, drho_s, tau_d, na, ir)
388 : ELSE
389 0 : CALL fill_rho_set(rho_set_h, lsd, mspins, needs, rho_h, rho_d, tau_d, na, ir)
390 0 : CALL fill_rho_set(rho_set_s, lsd, mspins, needs, rho_s, rho_d, tau_d, na, ir)
391 : END IF
392 : END DO
393 :
394 : ! hard atom density !
395 0 : CALL xc_dset_zero_all(deriv_set)
396 : CALL vxc_of_r_new(xc_fun_section, rho_set_h, deriv_set, 1, needs, weight_h, &
397 0 : lsd, na, nr, exc_h, vxc_h, vxg_h, vtau_h)
398 0 : IF (is_triplet) THEN
399 0 : vxc_h(:, :, 1) = vxc_h(:, :, 1) - vxc_h(:, :, 2)
400 0 : IF (gradient_f) THEN
401 0 : vxg_h(:, :, :, 1) = vxg_h(:, :, :, 1) - vxg_h(:, :, :, 2)
402 : END IF
403 0 : IF (tau_f) THEN
404 0 : vtau_h(:, :, 1) = vtau_h(:, :, 1) - vtau_h(:, :, 2)
405 : END IF
406 : END IF
407 : ! soft atom density !
408 0 : CALL xc_dset_zero_all(deriv_set)
409 : CALL vxc_of_r_new(xc_fun_section, rho_set_s, deriv_set, 1, needs, weight_s, &
410 0 : lsd, na, nr, exc_s, vxc_s, vxg_s, vtau_s)
411 0 : IF (is_triplet) THEN
412 0 : vxc_s(:, :, 1) = vxc_s(:, :, 1) - vxc_s(:, :, 2)
413 0 : IF (gradient_f) THEN
414 0 : vxg_s(:, :, :, 1) = vxg_s(:, :, :, 1) - vxg_s(:, :, :, 2)
415 : END IF
416 0 : IF (tau_f) THEN
417 0 : vtau_s(:, :, 1) = vtau_s(:, :, 1) - vtau_s(:, :, 2)
418 : END IF
419 : END IF
420 : ! potentials
421 0 : DO ns = 1, nspins
422 0 : fint_hh(ns)%r_coef(:, :) = 0.0_dp
423 0 : fint_ss(ns)%r_coef(:, :) = 0.0_dp
424 : END DO
425 0 : IF (gradient_f) THEN
426 : CALL gaVxcgb_GC(vxc_h, vxc_s, vxg_h, vxg_s, fint_hh, fint_ss, &
427 0 : grid_atom, basis_1c, harmonics, nspins)
428 : ELSE
429 : CALL gaVxcgb_noGC(vxc_h, vxc_s, fint_hh, fint_ss, &
430 0 : grid_atom, basis_1c, harmonics, nspins)
431 : END IF
432 0 : IF (tau_f) THEN
433 : CALL dgaVtaudgb(vtau_h, vtau_s, fint_hh, fint_ss, &
434 0 : tau_basis_cache, nspins)
435 : END IF
436 : ! first derivative fxc
437 0 : NULLIFY (int_hh, int_ss)
438 0 : CALL get_rho_atom(rho_atom=rho1_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss)
439 0 : DO ns = 1, nspins
440 0 : int_ss(ns)%r_coef(:, :) = int_ss(ns)%r_coef(:, :) + oeps1*ak(istep)*fint_ss(ns)%r_coef(:, :)
441 0 : int_hh(ns)%r_coef(:, :) = int_hh(ns)%r_coef(:, :) + oeps1*ak(istep)*fint_hh(ns)%r_coef(:, :)
442 : END DO
443 : ! second derivative gxc
444 0 : NULLIFY (int_hh, int_ss)
445 0 : CALL get_rho_atom(rho_atom=rho2_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss)
446 0 : DO ns = 1, nspins
447 0 : int_ss(ns)%r_coef(:, :) = int_ss(ns)%r_coef(:, :) + oeps2*bl(istep)*fint_ss(ns)%r_coef(:, :)
448 0 : int_hh(ns)%r_coef(:, :) = int_hh(ns)%r_coef(:, :) + oeps2*bl(istep)*fint_hh(ns)%r_coef(:, :)
449 : END DO
450 : END DO
451 : !
452 0 : DO ns = 1, nspins
453 0 : DEALLOCATE (fint_ss(ns)%r_coef)
454 0 : DEALLOCATE (fint_hh(ns)%r_coef)
455 : END DO
456 0 : DEALLOCATE (fint_ss, fint_hh)
457 :
458 : END DO ! iat
459 :
460 : ! Release the xc structure used to store the xc derivatives
461 0 : CALL xc_dset_release(deriv_set)
462 0 : CALL xc_rho_set_release(rho_set_h)
463 0 : CALL xc_rho_set_release(rho_set_s)
464 :
465 0 : DEALLOCATE (rho_h, rho_s, rho0_h, rho0_s, rho1_h, rho1_s)
466 0 : DEALLOCATE (vxc_h, vxc_s)
467 0 : IF (gradient_f) THEN
468 0 : DEALLOCATE (drho_h, drho_s, drho0_h, drho0_s, drho1_h, drho1_s)
469 0 : DEALLOCATE (vxg_h, vxg_s)
470 : END IF
471 0 : IF (tau_f) THEN
472 0 : DEALLOCATE (tau_h, tau_s, tau0_h, tau0_s, tau1_h, tau1_s)
473 0 : DEALLOCATE (vtau_h, vtau_s)
474 0 : CALL release_tau_basis_cache(tau_basis_cache)
475 : END IF
476 : END DO ! ikind
477 :
478 : END IF !xc_none
479 :
480 0 : CALL timestop(handle)
481 :
482 0 : END SUBROUTINE fgxc_atom_calc
483 :
484 : ! **************************************************************************************************
485 : !> \brief ...
486 : !> \param qs_env ...
487 : !> \param rho0_atom_set ...
488 : !> \param rho1_atom_set ...
489 : !> \param rho2_atom_set ...
490 : !> \param kind_set ...
491 : !> \param xc_section ...
492 : !> \param is_triplet ...
493 : !> \param accuracy ...
494 : !> \param epsrho ...
495 : ! **************************************************************************************************
496 114 : SUBROUTINE fgxc_atom_diff(qs_env, rho0_atom_set, rho1_atom_set, rho2_atom_set, &
497 : kind_set, xc_section, is_triplet, accuracy, epsrho)
498 :
499 : TYPE(qs_environment_type), POINTER :: qs_env
500 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho0_atom_set, rho1_atom_set, &
501 : rho2_atom_set
502 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: kind_set
503 : TYPE(section_vals_type), OPTIONAL, POINTER :: xc_section
504 : LOGICAL, INTENT(IN) :: is_triplet
505 : INTEGER, INTENT(IN) :: accuracy
506 : REAL(KIND=dp), INTENT(IN) :: epsrho
507 :
508 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fgxc_atom_diff'
509 :
510 : INTEGER :: bo(2), handle, iat, iatom, ikind, ir, &
511 : istep, mspins, myfun, na, natom, nf, &
512 : nr, ns, nspins, nstep, num_pe
513 : INTEGER, DIMENSION(2, 3) :: bounds
514 114 : INTEGER, DIMENSION(:), POINTER :: atom_list
515 : LOGICAL :: accint, donlcc, gradient_f, lsd, nlcc, &
516 : paw_atom, tau_f
517 : REAL(dp) :: agr, alpha, beta, density_cut, &
518 : gradient_cut, oeps1, tau_cut
519 114 : REAL(dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: vtau_h, vtau_s, vxc_h, vxc_s
520 : REAL(dp), DIMENSION(1, 1, 1) :: tau_d
521 : REAL(dp), DIMENSION(1, 1, 1, 1) :: rho_d
522 228 : REAL(dp), DIMENSION(:, :), POINTER :: rho_nlcc, weight_h, weight_s
523 228 : REAL(dp), DIMENSION(:, :, :), POINTER :: rho0_h, rho0_s, rho1_h, rho1_s, rho_h, &
524 228 : rho_s, tau0_h, tau0_s, tau1_h, tau1_s, &
525 114 : tau_h, tau_s
526 114 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: drho0_h, drho0_s, drho1_h, drho1_s, &
527 228 : drho_h, drho_s, vxg_h, vxg_s
528 : REAL(KIND=dp), DIMENSION(-4:4) :: ak
529 114 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
530 : TYPE(dft_control_type), POINTER :: dft_control
531 : TYPE(grid_atom_type), POINTER :: grid_atom
532 : TYPE(gto_basis_set_type), POINTER :: basis_1c
533 : TYPE(harmonics_atom_type), POINTER :: harmonics
534 : TYPE(mp_para_env_type), POINTER :: para_env
535 114 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s, fint_hh, fint_ss, int_hh, &
536 114 : int_ss, r_h, r_s
537 114 : TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r_h_d, r_s_d
538 : TYPE(rho_atom_type), POINTER :: rho0_atom, rho1_atom, rho2_atom
539 : TYPE(section_vals_type), POINTER :: xc_fun_section
540 114 : TYPE(tau_basis_cache_type) :: tau_basis_cache
541 : TYPE(xc_derivative_set_type) :: deriv_set
542 : TYPE(xc_rho_cflags_type) :: needs
543 : TYPE(xc_rho_set_type) :: rho1_set_h, rho1_set_s, rho_set_h, &
544 : rho_set_s
545 :
546 114 : CALL timeset(routineN, handle)
547 :
548 114 : NULLIFY (vtau_h, vtau_s)
549 :
550 114 : ak = 0.0_dp
551 114 : SELECT CASE (accuracy)
552 : CASE (:4)
553 0 : nstep = 2
554 0 : ak(-2:2) = [1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp]/12.0_dp
555 : CASE (5:7)
556 912 : nstep = 3
557 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
558 : CASE (8:)
559 0 : nstep = 4
560 : ak(-4:4) = [1.0_dp, -32.0_dp/3.0_dp, 56.0_dp, -224.0_dp, 0.0_dp, &
561 114 : 224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp]/280.0_dp
562 : END SELECT
563 114 : oeps1 = 1.0_dp/epsrho
564 :
565 : CALL get_qs_env(qs_env=qs_env, &
566 : dft_control=dft_control, &
567 : para_env=para_env, &
568 114 : atomic_kind_set=atomic_kind_set)
569 :
570 114 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
571 114 : CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
572 :
573 114 : accint = dft_control%qs_control%gapw_control%accurate_xcint
574 :
575 114 : IF (myfun == xc_none) THEN
576 : ! no action needed?
577 : ELSE
578 : ! calculate fxc
579 : CALL fxc_atom_calc(qs_env, rho0_atom_set, rho1_atom_set, xc_section, para_env, &
580 114 : do_triplet=is_triplet, kind_set_external=kind_set)
581 :
582 114 : CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", r_val=density_cut)
583 114 : CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", r_val=gradient_cut)
584 114 : CALL section_vals_val_get(xc_section, "TAU_CUTOFF", r_val=tau_cut)
585 :
586 114 : nlcc = has_nlcc(kind_set)
587 114 : lsd = dft_control%lsd
588 114 : nspins = dft_control%nspins
589 114 : mspins = nspins
590 114 : IF (is_triplet) THEN
591 12 : CPASSERT(nspins == 1)
592 12 : lsd = .TRUE.
593 12 : mspins = 2
594 : END IF
595 114 : needs = xc_functionals_get_needs(xc_fun_section, lsd=lsd, calc_potential=.TRUE.)
596 114 : gradient_f = (needs%drho .OR. needs%drho_spin)
597 114 : tau_f = (needs%tau .OR. needs%tau_spin)
598 :
599 : ! Here starts the loop over all the atoms
600 386 : DO ikind = 1, SIZE(atomic_kind_set)
601 272 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
602 : CALL get_qs_kind(kind_set(ikind), paw_atom=paw_atom, &
603 272 : harmonics=harmonics, grid_atom=grid_atom)
604 272 : CALL get_qs_kind(kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
605 :
606 272 : IF (.NOT. paw_atom) CYCLE
607 :
608 258 : nr = grid_atom%nr
609 258 : na = grid_atom%ng_sphere
610 :
611 : ! set integration weights
612 258 : IF (accint) THEN
613 138 : weight_h => grid_atom%weight
614 138 : alpha = dft_control%qs_control%gapw_control%aw(ikind)
615 138 : IF (ASSOCIATED(grid_atom%gapw_weight_s)) THEN
616 138 : IF (grid_atom%gapw_weight_alpha /= alpha) DEALLOCATE (grid_atom%gapw_weight_s)
617 : END IF
618 138 : IF (.NOT. ASSOCIATED(grid_atom%gapw_weight_s)) THEN
619 0 : ALLOCATE (grid_atom%gapw_weight_s(na, nr))
620 0 : DO ir = 1, nr
621 0 : agr = 1.0_dp - EXP(-alpha*grid_atom%rad2(ir))
622 0 : grid_atom%gapw_weight_s(:, ir) = grid_atom%weight(:, ir)*agr
623 : END DO
624 0 : grid_atom%gapw_weight_alpha = alpha
625 : END IF
626 138 : weight_s => grid_atom%gapw_weight_s
627 : ELSE
628 120 : weight_h => grid_atom%weight
629 120 : weight_s => grid_atom%weight
630 : END IF
631 :
632 : ! Prepare the structures needed to calculate and store the xc derivatives
633 :
634 : ! Array dimension: here anly one dimensional arrays are used,
635 : ! i.e. only the first column of deriv_data is read.
636 : ! The other two dimensions are set to size equal 1
637 2580 : bounds(1:2, 1:3) = 1
638 258 : bounds(2, 1) = na
639 258 : bounds(2, 2) = nr
640 :
641 : ! create a place where to put the derivatives
642 258 : CALL xc_dset_create(deriv_set, local_bounds=bounds)
643 : ! create the place where to store the argument for the functionals
644 : CALL xc_rho_set_create(rho_set_h, bounds, rho_cutoff=density_cut, &
645 258 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
646 : CALL xc_rho_set_create(rho_set_s, bounds, rho_cutoff=density_cut, &
647 258 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
648 : CALL xc_rho_set_create(rho1_set_h, bounds, rho_cutoff=density_cut, &
649 258 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
650 : CALL xc_rho_set_create(rho1_set_s, bounds, rho_cutoff=density_cut, &
651 258 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
652 :
653 : ! allocate the required 3d arrays where to store rho and drho
654 258 : CALL xc_rho_set_atom_update(rho_set_h, needs, mspins, bounds)
655 258 : CALL xc_rho_set_atom_update(rho_set_s, needs, mspins, bounds)
656 258 : CALL xc_rho_set_atom_update(rho1_set_h, needs, mspins, bounds)
657 258 : CALL xc_rho_set_atom_update(rho1_set_s, needs, mspins, bounds)
658 :
659 : ALLOCATE (rho_h(na, nr, nspins), rho_s(na, nr, nspins), &
660 : rho0_h(na, nr, nspins), rho0_s(na, nr, nspins), &
661 5160 : rho1_h(na, nr, nspins), rho1_s(na, nr, nspins))
662 1806 : ALLOCATE (vxc_h(na, nr, nspins), vxc_s(na, nr, nspins))
663 258 : IF (gradient_f) THEN
664 : ALLOCATE (drho_h(4, na, nr, nspins), drho_s(4, na, nr, nspins), &
665 : drho0_h(4, na, nr, nspins), drho0_s(4, na, nr, nspins), &
666 3520 : drho1_h(4, na, nr, nspins), drho1_s(4, na, nr, nspins))
667 1408 : ALLOCATE (vxg_h(3, na, nr, nspins), vxg_s(3, na, nr, nspins))
668 : END IF
669 258 : IF (tau_f) THEN
670 0 : CALL create_tau_basis_cache(tau_basis_cache, grid_atom, basis_1c, harmonics)
671 : ALLOCATE (tau_h(na, nr, nspins), tau_s(na, nr, nspins), &
672 : tau0_h(na, nr, nspins), tau0_s(na, nr, nspins), &
673 0 : tau1_h(na, nr, nspins), tau1_s(na, nr, nspins))
674 0 : ALLOCATE (vtau_h(na, nr, nspins), vtau_s(na, nr, nspins))
675 : END IF
676 : !
677 : ! NLCC: prepare rho and drho of the core charge for this KIND
678 258 : donlcc = .FALSE.
679 258 : IF (nlcc) THEN
680 0 : NULLIFY (rho_nlcc)
681 0 : rho_nlcc => kind_set(ikind)%nlcc_pot
682 0 : IF (ASSOCIATED(rho_nlcc)) donlcc = .TRUE.
683 : END IF
684 :
685 : ! Distribute the atoms of this kind
686 258 : num_pe = para_env%num_pe
687 258 : bo = get_limit(natom, num_pe, para_env%mepos)
688 :
689 441 : DO iat = bo(1), bo(2)
690 183 : iatom = atom_list(iat)
691 : !
692 183 : NULLIFY (int_hh, int_ss)
693 183 : rho0_atom => rho0_atom_set(iatom)
694 183 : CALL get_rho_atom(rho_atom=rho0_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss)
695 1098 : ALLOCATE (fint_ss(nspins), fint_hh(nspins))
696 366 : DO ns = 1, nspins
697 183 : nf = SIZE(int_ss(ns)%r_coef, 1)
698 732 : ALLOCATE (fint_ss(ns)%r_coef(nf, nf))
699 183 : nf = SIZE(int_hh(ns)%r_coef, 1)
700 915 : ALLOCATE (fint_hh(ns)%r_coef(nf, nf))
701 : END DO
702 :
703 : ! RHO0
704 467016 : rho0_h = 0.0_dp
705 467016 : rho0_s = 0.0_dp
706 183 : rho0_atom => rho0_atom_set(iatom)
707 183 : IF (gradient_f) THEN
708 126 : NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
709 : CALL get_rho_atom(rho_atom=rho0_atom, rho_rad_h=r_h, rho_rad_s=r_s, drho_rad_h=dr_h, &
710 126 : drho_rad_s=dr_s, rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
711 1581552 : drho0_h = 0.0_dp
712 1581552 : drho0_s = 0.0_dp
713 : ELSE
714 57 : NULLIFY (r_h, r_s)
715 57 : CALL get_rho_atom(rho_atom=rho0_atom, rho_rad_h=r_h, rho_rad_s=r_s)
716 57 : rho_d = 0.0_dp
717 : END IF
718 9333 : DO ir = 1, nr
719 : CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_f, &
720 : ir, r_h, r_s, rho0_h, rho0_s, dr_h, dr_s, &
721 9150 : r_h_d, r_s_d, drho0_h, drho0_s)
722 9333 : IF (donlcc) THEN
723 : CALL calc_rho_nlcc(grid_atom, nspins, gradient_f, &
724 0 : ir, rho_nlcc(:, 1), rho0_h, rho0_s, rho_nlcc(:, 2), drho0_h, drho0_s)
725 : END IF
726 : END DO
727 183 : IF (tau_f) THEN
728 : !compute tau on the grid all at once
729 0 : CALL calc_tau_atom(tau0_h, tau0_s, rho0_atom, tau_basis_cache, nspins)
730 : ELSE
731 183 : tau_d = 0.0_dp
732 : END IF
733 : ! RHO1
734 467016 : rho1_h = 0.0_dp
735 467016 : rho1_s = 0.0_dp
736 183 : rho1_atom => rho1_atom_set(iatom)
737 183 : IF (gradient_f) THEN
738 126 : NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
739 : CALL get_rho_atom(rho_atom=rho1_atom, rho_rad_h=r_h, rho_rad_s=r_s, drho_rad_h=dr_h, &
740 126 : drho_rad_s=dr_s, rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
741 1581552 : drho1_h = 0.0_dp
742 1581552 : drho1_s = 0.0_dp
743 : ELSE
744 57 : NULLIFY (r_h, r_s)
745 57 : CALL get_rho_atom(rho_atom=rho1_atom, rho_rad_h=r_h, rho_rad_s=r_s)
746 : END IF
747 9333 : DO ir = 1, nr
748 : CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_f, &
749 : ir, r_h, r_s, rho1_h, rho1_s, dr_h, dr_s, &
750 9333 : r_h_d, r_s_d, drho1_h, drho1_s)
751 : END DO
752 183 : IF (tau_f) THEN
753 : !compute tau on the grid all at once
754 0 : CALL calc_tau_atom(tau1_h, tau1_s, rho1_atom, tau_basis_cache, nspins)
755 : END IF
756 :
757 9333 : DO ir = 1, nr
758 9333 : IF (tau_f) THEN
759 0 : CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, drho1_h, tau1_h, na, ir)
760 0 : CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, drho1_s, tau1_s, na, ir)
761 9150 : ELSE IF (gradient_f) THEN
762 6300 : CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, drho1_h, tau_d, na, ir)
763 6300 : CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, drho1_s, tau_d, na, ir)
764 : ELSE
765 2850 : CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, rho_d, tau_d, na, ir)
766 2850 : CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, rho_d, tau_d, na, ir)
767 : END IF
768 : END DO
769 :
770 : ! RHO2
771 183 : rho2_atom => rho2_atom_set(iatom)
772 :
773 1464 : DO istep = -nstep, nstep
774 :
775 1281 : beta = REAL(istep, KIND=dp)*epsrho
776 :
777 6536943 : rho_h = rho0_h + beta*rho1_h
778 6536943 : rho_s = rho0_s + beta*rho1_s
779 1281 : IF (gradient_f) THEN
780 22140846 : drho_h = drho0_h + beta*drho1_h
781 22140846 : drho_s = drho0_s + beta*drho1_s
782 : END IF
783 1281 : IF (tau_f) THEN
784 0 : tau_h = tau0_h + beta*tau1_h
785 0 : tau_s = tau0_s + beta*tau1_s
786 : END IF
787 : !
788 1281 : IF (gradient_f) THEN
789 2250864 : drho_h(4, :, :, :) = NORM2(drho_h(1:3, :, :, :), 1)
790 2250864 : drho_s(4, :, :, :) = NORM2(drho_s(1:3, :, :, :), 1)
791 : END IF
792 :
793 65331 : DO ir = 1, nr
794 65331 : IF (tau_f) THEN
795 0 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_h, na, ir)
796 0 : CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_s, na, ir)
797 64050 : ELSE IF (gradient_f) THEN
798 44100 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_d, na, ir)
799 44100 : CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_d, na, ir)
800 : ELSE
801 19950 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, rho_d, tau_d, na, ir)
802 19950 : CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, rho_d, tau_d, na, ir)
803 : END IF
804 : END DO
805 :
806 : ! hard atom density !
807 1281 : CALL xc_dset_zero_all(deriv_set)
808 : CALL xc_2nd_deriv_of_r(xc_section=xc_section, &
809 : rho_set=rho_set_h, rho1_set=rho1_set_h, &
810 : deriv_set=deriv_set, &
811 : w=weight_h, vxc=vxc_h, vxg=vxg_h, vtau=vtau_h, &
812 1281 : do_triplet=is_triplet)
813 : ! soft atom density !
814 1281 : CALL xc_dset_zero_all(deriv_set)
815 : CALL xc_2nd_deriv_of_r(xc_section=xc_section, &
816 : rho_set=rho_set_s, rho1_set=rho1_set_s, &
817 : deriv_set=deriv_set, &
818 : w=weight_s, vxc=vxc_s, vxg=vxg_s, vtau=vtau_s, &
819 1281 : do_triplet=is_triplet)
820 : ! potentials
821 2562 : DO ns = 1, nspins
822 2410793 : fint_hh(ns)%r_coef(:, :) = 0.0_dp
823 2412074 : fint_ss(ns)%r_coef(:, :) = 0.0_dp
824 : END DO
825 1281 : IF (gradient_f) THEN
826 : CALL gaVxcgb_GC(vxc_h, vxc_s, vxg_h, vxg_s, fint_hh, fint_ss, &
827 882 : grid_atom, basis_1c, harmonics, nspins)
828 : ELSE
829 : CALL gaVxcgb_noGC(vxc_h, vxc_s, fint_hh, fint_ss, &
830 399 : grid_atom, basis_1c, harmonics, nspins)
831 : END IF
832 1281 : IF (tau_f) THEN
833 : CALL dgaVtaudgb(vtau_h, vtau_s, fint_hh, fint_ss, &
834 0 : tau_basis_cache, nspins)
835 : END IF
836 : ! second derivative gxc
837 1281 : NULLIFY (int_hh, int_ss)
838 1281 : CALL get_rho_atom(rho_atom=rho2_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss)
839 2745 : DO ns = 1, nspins
840 4820305 : int_ss(ns)%r_coef(:, :) = int_ss(ns)%r_coef(:, :) + oeps1*ak(istep)*fint_ss(ns)%r_coef(:, :)
841 4821586 : int_hh(ns)%r_coef(:, :) = int_hh(ns)%r_coef(:, :) + oeps1*ak(istep)*fint_hh(ns)%r_coef(:, :)
842 : END DO
843 : END DO
844 : !
845 366 : DO ns = 1, nspins
846 183 : DEALLOCATE (fint_ss(ns)%r_coef)
847 366 : DEALLOCATE (fint_hh(ns)%r_coef)
848 : END DO
849 441 : DEALLOCATE (fint_ss, fint_hh)
850 :
851 : END DO ! iat
852 :
853 : ! Release the xc structure used to store the xc derivatives
854 258 : CALL xc_dset_release(deriv_set)
855 258 : CALL xc_rho_set_release(rho_set_h)
856 258 : CALL xc_rho_set_release(rho_set_s)
857 258 : CALL xc_rho_set_release(rho1_set_h)
858 258 : CALL xc_rho_set_release(rho1_set_s)
859 :
860 258 : DEALLOCATE (rho_h, rho_s, rho0_h, rho0_s, rho1_h, rho1_s)
861 258 : DEALLOCATE (vxc_h, vxc_s)
862 258 : IF (gradient_f) THEN
863 176 : DEALLOCATE (drho_h, drho_s, drho0_h, drho0_s, drho1_h, drho1_s)
864 176 : DEALLOCATE (vxg_h, vxg_s)
865 : END IF
866 902 : IF (tau_f) THEN
867 0 : DEALLOCATE (tau_h, tau_s, tau0_h, tau0_s, tau1_h, tau1_s)
868 0 : DEALLOCATE (vtau_h, vtau_s)
869 0 : CALL release_tau_basis_cache(tau_basis_cache)
870 : END IF
871 : END DO ! ikind
872 :
873 : END IF !xc_none
874 :
875 114 : CALL timestop(handle)
876 :
877 8436 : END SUBROUTINE fgxc_atom_diff
878 :
879 : END MODULE qs_fgxc_atom
|
