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