(git:ae513ba)
Loading...
Searching...
No Matches
qs_fxc.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Setup Routine for Fxc Potentials
10! https://en.wikipedia.org/wiki/Finite_difference_coefficient
11!---------------------------------------------------------------------------------------------------
12!Derivative Accuracy 4 3 2 1 0 1 2 3 4
13!---------------------------------------------------------------------------------------------------
14! 1 2 -1/2 0 1/2
15! 4 1/12 -2/3 0 2/3 -1/12
16! 6 -1/60 3/20 -3/4 0 3/4 -3/20 1/60
17! 8 1/280 -4/105 1/5 -4/5 0 4/5 -1/5 4/105 -1/280
18!---------------------------------------------------------------------------------------------------
19! 2 2 1 -2 1
20! 4 -1/12 4/3 -5/2 4/3 -1/12
21! 6 1/90 -3/20 3/2 -49/18 3/2 -3/20 1/90
22! 8 -1/560 8/315 -1/5 8/5 -205/72 8/5 -1/5 8/315 -1/560
23!---------------------------------------------------------------------------------------------------
24!> \par History
25!> init 17.03.2020
26!> complete refactoring 08.2026
27!> \author JGH
28! **************************************************************************************************
29MODULE qs_fxc
30
31 USE cp_control_types, ONLY: dft_control_type
32 USE input_section_types, ONLY: section_get_ival,&
33 section_get_lval,&
34 section_get_rval,&
35 section_vals_get_subs_vals,&
36 section_vals_type
37 USE kinds, ONLY: dp
38 USE pw_env_types, ONLY: pw_env_get,&
39 pw_env_type
40 USE pw_grids, ONLY: pw_grid_compare
41 USE pw_methods, ONLY: pw_axpy,&
42 pw_scale,&
43 pw_transfer,&
44 pw_zero
45 USE pw_pool_types, ONLY: pw_pool_type
46 USE pw_types, ONLY: pw_c1d_gs_type,&
47 pw_r3d_rs_type
48 USE qs_environment_types, ONLY: get_qs_env,&
49 qs_environment_type
50 USE qs_ks_types, ONLY: qs_ks_env_type
51 USE qs_rho_methods, ONLY: qs_rho_copy,&
52 qs_rho_scale_and_add,&
53 qs_rho_transfer
54 USE qs_rho_types, ONLY: qs_rho_create,&
55 qs_rho_get,&
56 qs_rho_release,&
57 qs_rho_type
58 USE qs_vxc, ONLY: qs_vxc_create
59 USE xc, ONLY: xc_calc_2nd_deriv_analytical,&
60 xc_calc_2nd_deriv_numerical,&
61 xc_prep_2nd_deriv
62 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
63 xc_dset_release
64 USE xc_derivatives, ONLY: xc_functionals_get_needs
65 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
66 USE xc_rho_set_types, ONLY: xc_rho_set_create,&
67 xc_rho_set_release,&
68 xc_rho_set_type,&
69 xc_rho_set_update
70#include "./base/base_uses.f90"
71
72 IMPLICIT NONE
73
74 PRIVATE
75
76 ! *** Public subroutines ***
77 PUBLIC :: qs_fxc_create, qs_fxc_prep, qs_fxc_apply
78 PUBLIC :: qs_fxc_fdiff
79
80 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_fxc'
81
82! **************************************************************************************************
83
84CONTAINS
85
86! **************************************************************************************************
87!> \brief ...
88!> \param qs_env ...
89!> \param rho0_struct ...
90!> \param rho1_struct ...
91!> \param xc_section ...
92!> \param fxc_rho ...
93!> \param fxc_tau ...
94!> \param is_triplet ...
95!> \param spinflip ...
96!> \param no_weights ...
97!> \param uf_grid_results ...
98!> \param pw_env_ext ...
99!> \param compute_virial ...
100!> \param virial_xc ...
101! **************************************************************************************************
102 SUBROUTINE qs_fxc_create(qs_env, rho0_struct, rho1_struct, xc_section, fxc_rho, fxc_tau, &
103 is_triplet, spinflip, no_weights, uf_grid_results, &
104 pw_env_ext, compute_virial, virial_xc)
105
106 TYPE(qs_environment_type), POINTER :: qs_env
107 TYPE(qs_rho_type), POINTER :: rho0_struct, rho1_struct
108 TYPE(section_vals_type), POINTER :: xc_section
109 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau
110 LOGICAL, INTENT(IN), OPTIONAL :: is_triplet, spinflip, no_weights, &
111 uf_grid_results
112 TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_ext
113 LOGICAL, INTENT(IN), OPTIONAL :: compute_virial
114 REAL(kind=dp), DIMENSION(3, 3), INTENT(INOUT), &
115 OPTIONAL :: virial_xc
116
117 CHARACTER(len=*), PARAMETER :: routinen = 'qs_fxc_create'
118
119 INTEGER :: handle, ispin, nspins
120 LOGICAL :: do_virial, do_w, ret_uf, uf_grid
121 REAL(kind=dp) :: factor
122 TYPE(dft_control_type), POINTER :: dft_control
123 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
124 TYPE(pw_c1d_gs_type), POINTER :: rho_nlcc_g
125 TYPE(pw_env_type), POINTER :: pw_env
126 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool, xc_pw_pool
127 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho_lo, fxc_rho_uf, fxc_tau_lo, &
128 fxc_tau_uf, rho_r
129 TYPE(pw_r3d_rs_type), POINTER :: rho_nlcc, weights, weights_uf
130 TYPE(qs_rho_type), POINTER :: rho0_uf, rho1_uf
131
132 CALL timeset(routinen, handle)
133
134 do_virial = .false.
135 IF (PRESENT(compute_virial)) do_virial = compute_virial
136
137 CALL get_qs_env(qs_env, dft_control=dft_control)
138
139 IF (PRESENT(pw_env_ext)) THEN
140 pw_env => pw_env_ext
141 ELSE
142 CALL get_qs_env(qs_env, pw_env=pw_env)
143 END IF
144 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, xc_pw_pool=xc_pw_pool)
145 uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid, xc_pw_pool%pw_grid)
146
147 nspins = dft_control%nspins
148 IF (ASSOCIATED(fxc_rho)) THEN
149 cpassert(nspins == SIZE(fxc_rho))
150 END IF
151 IF (ASSOCIATED(fxc_tau)) THEN
152 cpassert(nspins == SIZE(fxc_tau))
153 END IF
154
155 NULLIFY (rho_nlcc, rho_nlcc_g)
156 CALL get_qs_env(qs_env, rho_nlcc=rho_nlcc, rho_nlcc_g=rho_nlcc_g)
157 IF (ASSOCIATED(rho_nlcc)) THEN
158 NULLIFY (rho_r, rho_g)
159 CALL qs_rho_get(rho0_struct, rho_r=rho_r, rho_g=rho_g)
160 factor = 1.0_dp
161 DO ispin = 1, nspins
162 CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
163 CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
164 END DO
165 END IF
166
167 do_w = .true.
168 IF (PRESENT(no_weights)) do_w = .NOT. no_weights
169 IF (do_w) THEN
170 CALL get_qs_env(qs_env, xcint_weights=weights)
171 ELSE
172 NULLIFY (weights)
173 END IF
174
175 NULLIFY (fxc_rho_lo, fxc_tau_lo)
176 IF (uf_grid) THEN
177 IF (PRESENT(uf_grid_results)) THEN
178 ret_uf = uf_grid_results
179 ELSE
180 ret_uf = .false.
181 END IF
182 NULLIFY (weights_uf)
183 IF (ASSOCIATED(weights)) THEN
184 ALLOCATE (weights_uf)
185 CALL xc_pw_pool%create_pw(weights_uf)
186 block
187 TYPE(pw_c1d_gs_type) :: weights_g, weights_g_uf
188 CALL auxbas_pw_pool%create_pw(weights_g)
189 CALL xc_pw_pool%create_pw(weights_g_uf)
190 CALL pw_transfer(weights, weights_g)
191 CALL pw_transfer(weights_g, weights_g_uf)
192 CALL pw_transfer(weights_g_uf, weights_uf)
193 CALL xc_pw_pool%give_back_pw(weights_g_uf)
194 CALL auxbas_pw_pool%give_back_pw(weights_g)
195 END block
196 END IF
197 !
198 ALLOCATE (rho0_uf, rho1_uf)
199 CALL qs_rho_create(rho0_uf)
200 CALL qs_rho_create(rho1_uf)
201 CALL qs_rho_transfer(rho0_struct, rho0_uf, auxbas_pw_pool, xc_pw_pool)
202 CALL qs_rho_transfer(rho1_struct, rho1_uf, auxbas_pw_pool, xc_pw_pool)
203 !
204 NULLIFY (fxc_rho_uf, fxc_tau_uf)
205 CALL qs_fxc_calculate(rho0_uf, rho1_uf, xc_section, weights_uf, xc_pw_pool, &
206 fxc_rho_uf, fxc_tau_uf, &
207 is_triplet=is_triplet, spinflip=spinflip, &
208 compute_virial=do_virial, virial_xc=virial_xc)
209 !
210 CALL qs_rho_release(rho0_uf)
211 CALL qs_rho_release(rho1_uf)
212 DEALLOCATE (rho0_uf, rho1_uf)
213 IF (ASSOCIATED(weights_uf)) THEN
214 CALL xc_pw_pool%give_back_pw(weights_uf)
215 DEALLOCATE (weights_uf)
216 END IF
217 IF (ret_uf) THEN
218 fxc_rho_lo => fxc_rho_uf
219 fxc_tau_lo => fxc_tau_uf
220 ELSE
221 IF (ASSOCIATED(fxc_rho_uf)) THEN
222 ALLOCATE (fxc_rho_lo(nspins))
223 DO ispin = 1, nspins
224 CALL auxbas_pw_pool%create_pw(fxc_rho_lo(ispin))
225 block
226 TYPE(pw_c1d_gs_type) :: fxc_g, fxc_g_uf
227 CALL auxbas_pw_pool%create_pw(fxc_g)
228 CALL xc_pw_pool%create_pw(fxc_g_uf)
229 CALL pw_transfer(fxc_rho_uf(ispin), fxc_g_uf)
230 CALL pw_transfer(fxc_g_uf, fxc_g)
231 CALL pw_transfer(fxc_g, fxc_rho_lo(ispin))
232 CALL xc_pw_pool%give_back_pw(fxc_g_uf)
233 CALL auxbas_pw_pool%give_back_pw(fxc_g)
234 END block
235 CALL xc_pw_pool%give_back_pw(fxc_rho_uf(ispin))
236 END DO
237 DEALLOCATE (fxc_rho_uf)
238 END IF
239 IF (ASSOCIATED(fxc_tau_uf)) THEN
240 ALLOCATE (fxc_tau_lo(nspins))
241 DO ispin = 1, nspins
242 CALL auxbas_pw_pool%create_pw(fxc_tau_lo(ispin))
243 block
244 TYPE(pw_c1d_gs_type) :: fxc_g, fxc_g_uf
245 CALL auxbas_pw_pool%create_pw(fxc_g)
246 CALL xc_pw_pool%create_pw(fxc_g_uf)
247 CALL pw_transfer(fxc_tau_uf(ispin), fxc_g_uf)
248 CALL pw_transfer(fxc_g_uf, fxc_g)
249 CALL pw_transfer(fxc_g, fxc_tau_lo(ispin))
250 CALL xc_pw_pool%give_back_pw(fxc_g_uf)
251 CALL auxbas_pw_pool%give_back_pw(fxc_g)
252 END block
253 CALL xc_pw_pool%give_back_pw(fxc_tau_uf(ispin))
254 END DO
255 END IF
256 END IF
257
258 ELSE
259 CALL qs_fxc_calculate(rho0_struct, rho1_struct, xc_section, weights, auxbas_pw_pool, &
260 fxc_rho_lo, fxc_tau_lo, &
261 is_triplet=is_triplet, spinflip=spinflip, &
262 compute_virial=do_virial, virial_xc=virial_xc)
263 END IF
264
265 ! de-apply NLCC density
266 IF (ASSOCIATED(rho_nlcc)) THEN
267 factor = -1.0_dp
268 DO ispin = 1, nspins
269 CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
270 CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
271 END DO
272 END IF
273
274 ! return potentials
275 IF (ASSOCIATED(fxc_rho)) THEN
276 DO ispin = 1, min(SIZE(fxc_rho_lo), SIZE(fxc_rho))
277 CALL pw_transfer(fxc_rho_lo(ispin), fxc_rho(ispin))
278 END DO
279 DO ispin = 1, SIZE(fxc_rho_lo)
280 CALL auxbas_pw_pool%give_back_pw(fxc_rho_lo(ispin))
281 END DO
282 DEALLOCATE (fxc_rho_lo)
283 ELSE
284 fxc_rho => fxc_rho_lo
285 END IF
286 IF (ASSOCIATED(fxc_tau)) THEN
287 IF (ASSOCIATED(fxc_tau_lo)) THEN
288 DO ispin = 1, min(SIZE(fxc_tau_lo), SIZE(fxc_tau))
289 CALL pw_transfer(fxc_tau_lo(ispin), fxc_tau(ispin))
290 END DO
291 DO ispin = 1, SIZE(fxc_tau_lo)
292 CALL auxbas_pw_pool%give_back_pw(fxc_tau_lo(ispin))
293 END DO
294 DEALLOCATE (fxc_tau_lo)
295 ELSE
296 DO ispin = 1, nspins
297 CALL pw_zero(fxc_tau(ispin))
298 END DO
299 END IF
300 ELSE
301 fxc_tau => fxc_tau_lo
302 END IF
303
304 CALL timestop(handle)
305
306 END SUBROUTINE qs_fxc_create
307
308! **************************************************************************************************
309!> \brief ...
310!> \param rho0 ...
311!> \param rho1 ...
312!> \param xc_section ...
313!> \param weights ...
314!> \param auxbas_pw_pool ...
315!> \param fxc_rho ...
316!> \param fxc_tau ...
317!> \param is_triplet ...
318!> \param spinflip ...
319!> \param compute_virial ...
320!> \param virial_xc ...
321! **************************************************************************************************
322 SUBROUTINE qs_fxc_calculate(rho0, rho1, xc_section, weights, auxbas_pw_pool, &
323 fxc_rho, fxc_tau, is_triplet, spinflip, &
324 compute_virial, virial_xc)
325
326 TYPE(qs_rho_type), POINTER :: rho0, rho1
327 TYPE(section_vals_type), POINTER :: xc_section
328 TYPE(pw_r3d_rs_type), POINTER :: weights
329 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
330 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau
331 LOGICAL, INTENT(IN), OPTIONAL :: is_triplet, spinflip, compute_virial
332 REAL(kind=dp), DIMENSION(3, 3), INTENT(INOUT), &
333 OPTIONAL :: virial_xc
334
335 CHARACTER(len=*), PARAMETER :: routinen = 'qs_fxc_calculate'
336
337 INTEGER :: handle, ispin, mspins, nspins
338 INTEGER, DIMENSION(2, 3) :: bo
339 LOGICAL :: do_analytic, do_sf, do_triplet, &
340 do_virial, lsd
341 REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: vxg
342 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho0_g, rho1_g
343 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho0_r, rho1_r, tau0_r, tau1_r
344 TYPE(qs_rho_type), POINTER :: rhot0, rhot1
345 TYPE(section_vals_type), POINTER :: xc_fun_section
346 TYPE(xc_derivative_set_type) :: xc_deriv_set
347 TYPE(xc_rho_cflags_type) :: needs
348 TYPE(xc_rho_set_type) :: rho0_set, rho1_set
349
350 CALL timeset(routinen, handle)
351
352 do_triplet = .false.
353 IF (PRESENT(is_triplet)) do_triplet = is_triplet
354
355 do_sf = .false.
356 IF (PRESENT(spinflip)) do_sf = spinflip
357
358 do_virial = .false.
359 IF (PRESENT(compute_virial)) do_virial = compute_virial
360
361 do_analytic = section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL")
362
363 CALL qs_rho_get(rho0, rho_r=rho0_r, rho_g=rho0_g, tau_r=tau0_r)
364 CALL qs_rho_get(rho1, rho_r=rho1_r, tau_r=tau1_r)
365 NULLIFY (rho1_g)
366
367 mspins = SIZE(rho0_r)
368 nspins = SIZE(rho0_r)
369 lsd = (nspins == 2)
370 IF (nspins == 1 .AND. do_triplet) THEN
371 nspins = 2
372 lsd = .true.
373 ELSE IF (do_sf) THEN
374 nspins = 1
375 mspins = 1
376 lsd = .true.
377 END IF
378
379 cpassert(.NOT. ASSOCIATED(fxc_rho))
380 cpassert(.NOT. ASSOCIATED(fxc_tau))
381 xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
382 needs = xc_functionals_get_needs(xc_fun_section, lsd, .true.)
383 ALLOCATE (fxc_rho(mspins))
384 DO ispin = 1, mspins
385 CALL auxbas_pw_pool%create_pw(fxc_rho(ispin))
386 CALL pw_zero(fxc_rho(ispin))
387 END DO
388 IF (needs%tau .OR. needs%tau_spin) THEN
389 IF (.NOT. ASSOCIATED(tau1_r)) THEN
390 cpabort("Tau-dependent functionals requires allocated kinetic energy density grid")
391 END IF
392 ALLOCATE (fxc_tau(mspins))
393 DO ispin = 1, mspins
394 CALL auxbas_pw_pool%create_pw(fxc_tau(ispin))
395 CALL pw_zero(fxc_tau(ispin))
396 END DO
397 END IF
398
399 IF (mspins == 1 .AND. do_triplet) THEN
400 ! split the density and response density arrays for triplet calculation
401 ALLOCATE (rhot0)
402 CALL qs_rho_create(rhot0)
403 CALL qs_rho_copy(rho0, rhot0, auxbas_pw_pool, 2, factor=2.0_dp)
404 !
405 ALLOCATE (rhot1)
406 CALL qs_rho_create(rhot1)
407 CALL qs_rho_copy(rho1, rhot1, auxbas_pw_pool, 2, factor=2.0_dp)
408 !
409 CALL qs_rho_get(rhot0, rho_r=rho0_r, rho_g=rho0_g, tau_r=tau0_r)
410 CALL qs_rho_get(rhot1, rho_r=rho1_r, tau_r=tau1_r)
411
412 CALL xc_prep_2nd_deriv(xc_deriv_set, rho0_set, rho0_r, auxbas_pw_pool, weights, &
413 xc_section=xc_section, tau_r=tau0_r)
414 bo = rho1_r(1)%pw_grid%bounds_local
415 ! create the place where to store the argument for the functionals
416 CALL xc_rho_set_create(rho1_set, bo, &
417 rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
418 drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
419 tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
420
421 ! calculate the arguments needed by the functionals
422 CALL xc_rho_set_update(rho1_set, rho1_r, rho1_g, tau1_r, needs, &
423 section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
424 section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
425 auxbas_pw_pool, spinflip=do_sf)
426 ELSE
427 CALL xc_prep_2nd_deriv(xc_deriv_set, rho0_set, rho0_r, auxbas_pw_pool, weights, &
428 xc_section=xc_section, tau_r=tau0_r)
429 bo = rho1_r(1)%pw_grid%bounds_local
430 ! create the place where to store the argument for the functionals
431 CALL xc_rho_set_create(rho1_set, bo, &
432 rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
433 drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
434 tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
435
436 ! calculate the arguments needed by the functionals
437 CALL xc_rho_set_update(rho1_set, rho1_r, rho1_g, tau1_r, needs, &
438 section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
439 section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
440 auxbas_pw_pool, spinflip=do_sf)
441 END IF
442
443 IF (mspins == 1 .AND. do_triplet .AND. do_analytic) THEN
444
445 CALL xc_calc_2nd_deriv_analytical(fxc_rho, fxc_tau, xc_deriv_set, rho0_set, &
446 rho1_set, auxbas_pw_pool, xc_section, &
447 gapw=.false., vxg=vxg, tddfpt_fac=-1.0_dp, spinflip=do_sf, &
448 compute_virial=compute_virial, virial_xc=virial_xc)
449
450 ELSE IF (do_analytic) THEN
451
452 CALL xc_calc_2nd_deriv_analytical(fxc_rho, fxc_tau, xc_deriv_set, rho0_set, &
453 rho1_set, auxbas_pw_pool, xc_section, &
454 gapw=.false., vxg=vxg, spinflip=do_sf, &
455 compute_virial=compute_virial, virial_xc=virial_xc)
456
457 ELSE
458
459 CALL xc_calc_2nd_deriv_numerical(fxc_rho, fxc_tau, rho0_set, rho1_r, rho1_g, tau1_r, &
460 auxbas_pw_pool, weights, xc_section, &
461 do_triplet, compute_virial, virial_xc, xc_deriv_set)
462
463 END IF
464
465 IF (mspins == 1 .AND. do_triplet) THEN
466 CALL qs_rho_release(rhot0)
467 DEALLOCATE (rhot0)
468 CALL qs_rho_release(rhot1)
469 DEALLOCATE (rhot1)
470 END IF
471
472 CALL xc_dset_release(xc_deriv_set)
473 CALL xc_rho_set_release(rho0_set)
474 CALL xc_rho_set_release(rho1_set)
475
476 CALL timestop(handle)
477
478 END SUBROUTINE qs_fxc_calculate
479
480! **************************************************************************************************
481!> \brief ...
482!> \param qs_env ...
483!> \param rho0_struct ...
484!> \param xc_rho_set ...
485!> \param xc_deriv_set ...
486!> \param xc_section ...
487!> \param pw_env_ext ...
488!> \param is_triplet ...
489! **************************************************************************************************
490 SUBROUTINE qs_fxc_prep(qs_env, rho0_struct, xc_rho_set, xc_deriv_set, &
491 xc_section, pw_env_ext, is_triplet)
492
493 TYPE(qs_environment_type), POINTER :: qs_env
494 TYPE(qs_rho_type), POINTER :: rho0_struct
495 TYPE(xc_rho_set_type) :: xc_rho_set
496 TYPE(xc_derivative_set_type) :: xc_deriv_set
497 TYPE(section_vals_type), POINTER :: xc_section
498 TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_ext
499 LOGICAL, INTENT(IN), OPTIONAL :: is_triplet
500
501 CHARACTER(len=*), PARAMETER :: routinen = 'qs_fxc_prep'
502
503 INTEGER :: handle, ispin, nspins
504 LOGICAL :: uf_grid
505 REAL(kind=dp) :: factor
506 TYPE(dft_control_type), POINTER :: dft_control
507 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
508 TYPE(pw_c1d_gs_type), POINTER :: rho_nlcc_g
509 TYPE(pw_env_type), POINTER :: pw_env
510 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool, xc_pw_pool
511 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
512 TYPE(pw_r3d_rs_type), POINTER :: rho_nlcc, weights, weights_uf
513 TYPE(qs_rho_type), POINTER :: rho0_uf
514
515 CALL timeset(routinen, handle)
516
517 CALL get_qs_env(qs_env, dft_control=dft_control)
518
519 IF (PRESENT(pw_env_ext)) THEN
520 pw_env => pw_env_ext
521 ELSE
522 CALL get_qs_env(qs_env, pw_env=pw_env)
523 END IF
524 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, xc_pw_pool=xc_pw_pool)
525 uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid, xc_pw_pool%pw_grid)
526
527 nspins = dft_control%nspins
528
529 NULLIFY (rho_nlcc, rho_nlcc_g)
530 CALL get_qs_env(qs_env, rho_nlcc=rho_nlcc, rho_nlcc_g=rho_nlcc_g)
531 IF (ASSOCIATED(rho_nlcc)) THEN
532 NULLIFY (rho_r, rho_g)
533 CALL qs_rho_get(rho0_struct, rho_r=rho_r, rho_g=rho_g)
534 factor = 1.0_dp
535 DO ispin = 1, nspins
536 CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
537 CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
538 END DO
539 END IF
540
541 NULLIFY (weights)
542 CALL get_qs_env(qs_env, xcint_weights=weights)
543
544 IF (uf_grid) THEN
545 NULLIFY (weights_uf)
546 IF (ASSOCIATED(weights)) THEN
547 ALLOCATE (weights_uf)
548 CALL xc_pw_pool%create_pw(weights_uf)
549 block
550 TYPE(pw_c1d_gs_type) :: weights_g, weights_g_uf
551 CALL auxbas_pw_pool%create_pw(weights_g)
552 CALL xc_pw_pool%create_pw(weights_g_uf)
553 CALL pw_transfer(weights, weights_g)
554 CALL pw_transfer(weights_g, weights_g_uf)
555 CALL pw_transfer(weights_g_uf, weights_uf)
556 CALL xc_pw_pool%give_back_pw(weights_g_uf)
557 CALL auxbas_pw_pool%give_back_pw(weights_g)
558 END block
559 END IF
560 !
561 ALLOCATE (rho0_uf)
562 CALL qs_rho_create(rho0_uf)
563 CALL qs_rho_transfer(rho0_struct, rho0_uf, auxbas_pw_pool, xc_pw_pool)
564 !
565 CALL qs_fxc_deriv(rho0_uf, xc_rho_set, xc_deriv_set, &
566 xc_section, weights_uf, xc_pw_pool, is_triplet)
567 !
568 CALL qs_rho_release(rho0_uf)
569 DEALLOCATE (rho0_uf)
570 IF (ASSOCIATED(weights_uf)) THEN
571 CALL xc_pw_pool%give_back_pw(weights_uf)
572 DEALLOCATE (weights_uf)
573 END IF
574 ELSE
575 CALL qs_fxc_deriv(rho0_struct, xc_rho_set, xc_deriv_set, &
576 xc_section, weights, auxbas_pw_pool, is_triplet)
577 END IF
578
579 ! de-apply NLCC density
580 IF (ASSOCIATED(rho_nlcc)) THEN
581 factor = -1.0_dp
582 DO ispin = 1, nspins
583 CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
584 CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
585 END DO
586 END IF
587
588 CALL timestop(handle)
589
590 END SUBROUTINE qs_fxc_prep
591
592! **************************************************************************************************
593!> \brief ...
594!> \param rho0 ...
595!> \param xc_rho_set ...
596!> \param xc_deriv_set ...
597!> \param xc_section ...
598!> \param weights ...
599!> \param auxbas_pw_pool ...
600!> \param is_triplet ...
601! **************************************************************************************************
602 SUBROUTINE qs_fxc_deriv(rho0, xc_rho_set, xc_deriv_set, xc_section, weights, auxbas_pw_pool, &
603 is_triplet)
604
605 TYPE(qs_rho_type), POINTER :: rho0
606 TYPE(xc_rho_set_type) :: xc_rho_set
607 TYPE(xc_derivative_set_type) :: xc_deriv_set
608 TYPE(section_vals_type), POINTER :: xc_section
609 TYPE(pw_r3d_rs_type), POINTER :: weights
610 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
611 LOGICAL, INTENT(IN) :: is_triplet
612
613 CHARACTER(len=*), PARAMETER :: routinen = 'qs_fxc_deriv'
614
615 INTEGER :: handle
616 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho0_r, tau0_r
617 TYPE(qs_rho_type), POINTER :: rhot0
618
619 CALL timeset(routinen, handle)
620
621 NULLIFY (rho0_r, tau0_r)
622 IF (is_triplet) THEN
623 ! split the density and response density arrays for triplet calculation
624 ALLOCATE (rhot0)
625 CALL qs_rho_create(rhot0)
626 CALL qs_rho_copy(rho0, rhot0, auxbas_pw_pool, 2, factor=2.0_dp)
627 CALL qs_rho_get(rhot0, rho_r=rho0_r, tau_r=tau0_r)
628 CALL xc_prep_2nd_deriv(xc_deriv_set, xc_rho_set, rho0_r, auxbas_pw_pool, weights, &
629 xc_section=xc_section, tau_r=tau0_r)
630 CALL qs_rho_release(rhot0)
631 DEALLOCATE (rhot0)
632 ELSE
633 CALL qs_rho_get(rho0, rho_r=rho0_r, tau_r=tau0_r)
634 CALL xc_prep_2nd_deriv(xc_deriv_set, xc_rho_set, rho0_r, auxbas_pw_pool, weights, &
635 xc_section=xc_section, tau_r=tau0_r)
636 END IF
637
638 CALL timestop(handle)
639
640 END SUBROUTINE qs_fxc_deriv
641
642! **************************************************************************************************
643!> \brief ...
644!> \param qs_env ...
645!> \param xc_deriv_set ...
646!> \param xc_rho_set ...
647!> \param rho1_struct ...
648!> \param xc_section ...
649!> \param fxc_rho ...
650!> \param fxc_tau ...
651!> \param is_triplet ...
652!> \param spinflip ...
653!> \param pw_env_ext ...
654!> \param compute_virial ...
655!> \param virial_xc ...
656! **************************************************************************************************
657 SUBROUTINE qs_fxc_apply(qs_env, xc_deriv_set, xc_rho_set, rho1_struct, xc_section, fxc_rho, fxc_tau, &
658 is_triplet, spinflip, pw_env_ext, compute_virial, virial_xc)
659
660 TYPE(qs_environment_type), POINTER :: qs_env
661 TYPE(xc_derivative_set_type) :: xc_deriv_set
662 TYPE(xc_rho_set_type) :: xc_rho_set
663 TYPE(qs_rho_type), POINTER :: rho1_struct
664 TYPE(section_vals_type), POINTER :: xc_section
665 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau
666 LOGICAL, INTENT(IN), OPTIONAL :: is_triplet, spinflip
667 TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_ext
668 LOGICAL, INTENT(IN), OPTIONAL :: compute_virial
669 REAL(kind=dp), DIMENSION(3, 3), INTENT(INOUT), &
670 OPTIONAL :: virial_xc
671
672 CHARACTER(len=*), PARAMETER :: routinen = 'qs_fxc_apply'
673
674 INTEGER :: handle, ispin, nspins
675 LOGICAL :: do_virial, uf_grid
676 TYPE(dft_control_type), POINTER :: dft_control
677 TYPE(pw_env_type), POINTER :: pw_env
678 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool, xc_pw_pool
679 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho_lo, fxc_rho_uf, fxc_tau_lo, &
680 fxc_tau_uf
681 TYPE(pw_r3d_rs_type), POINTER :: weights, weights_uf
682 TYPE(qs_rho_type), POINTER :: rho1_uf
683
684 CALL timeset(routinen, handle)
685
686 do_virial = .false.
687 IF (PRESENT(compute_virial)) do_virial = compute_virial
688
689 CALL get_qs_env(qs_env, dft_control=dft_control)
690
691 IF (PRESENT(pw_env_ext)) THEN
692 pw_env => pw_env_ext
693 ELSE
694 CALL get_qs_env(qs_env, pw_env=pw_env)
695 END IF
696 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, xc_pw_pool=xc_pw_pool)
697 uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid, xc_pw_pool%pw_grid)
698
699 nspins = dft_control%nspins
700 IF (ASSOCIATED(fxc_rho)) THEN
701 cpassert(nspins == SIZE(fxc_rho))
702 END IF
703 IF (ASSOCIATED(fxc_tau)) THEN
704 cpassert(nspins == SIZE(fxc_tau))
705 END IF
706
707 CALL get_qs_env(qs_env, xcint_weights=weights)
708
709 NULLIFY (fxc_rho_lo, fxc_tau_lo)
710 IF (uf_grid) THEN
711 NULLIFY (weights_uf)
712 IF (ASSOCIATED(weights)) THEN
713 ALLOCATE (weights_uf)
714 CALL xc_pw_pool%create_pw(weights_uf)
715 block
716 TYPE(pw_c1d_gs_type) :: weights_g, weights_g_uf
717 CALL auxbas_pw_pool%create_pw(weights_g)
718 CALL xc_pw_pool%create_pw(weights_g_uf)
719 CALL pw_transfer(weights, weights_g)
720 CALL pw_transfer(weights_g, weights_g_uf)
721 CALL pw_transfer(weights_g_uf, weights_uf)
722 CALL xc_pw_pool%give_back_pw(weights_g_uf)
723 CALL auxbas_pw_pool%give_back_pw(weights_g)
724 END block
725 END IF
726 !
727 ALLOCATE (rho1_uf)
728 CALL qs_rho_create(rho1_uf)
729 CALL qs_rho_transfer(rho1_struct, rho1_uf, auxbas_pw_pool, xc_pw_pool)
730 !
731 NULLIFY (fxc_rho_uf, fxc_tau_uf)
732 CALL qs_fxc_eval(xc_deriv_set, xc_rho_set, rho1_uf, xc_section, &
733 weights_uf, xc_pw_pool, fxc_rho_uf, fxc_tau_uf, &
734 is_triplet=is_triplet, spinflip=spinflip, &
735 compute_virial=do_virial, virial_xc=virial_xc)
736 !
737 CALL qs_rho_release(rho1_uf)
738 DEALLOCATE (rho1_uf)
739 IF (ASSOCIATED(weights_uf)) THEN
740 CALL xc_pw_pool%give_back_pw(weights_uf)
741 DEALLOCATE (weights_uf)
742 END IF
743 IF (ASSOCIATED(fxc_rho_uf)) THEN
744 ALLOCATE (fxc_rho_lo(nspins))
745 DO ispin = 1, nspins
746 CALL auxbas_pw_pool%create_pw(fxc_rho_lo(ispin))
747 block
748 TYPE(pw_c1d_gs_type) :: fxc_g, fxc_g_uf
749 CALL auxbas_pw_pool%create_pw(fxc_g)
750 CALL xc_pw_pool%create_pw(fxc_g_uf)
751 CALL pw_transfer(fxc_rho_uf(ispin), fxc_g_uf)
752 CALL pw_transfer(fxc_g_uf, fxc_g)
753 CALL pw_transfer(fxc_g, fxc_rho_lo(ispin))
754 CALL xc_pw_pool%give_back_pw(fxc_g_uf)
755 CALL auxbas_pw_pool%give_back_pw(fxc_g)
756 END block
757 CALL xc_pw_pool%give_back_pw(fxc_rho_uf(ispin))
758 END DO
759 DEALLOCATE (fxc_rho_uf)
760 END IF
761 IF (ASSOCIATED(fxc_tau_uf)) THEN
762 ALLOCATE (fxc_tau_lo(nspins))
763 DO ispin = 1, nspins
764 CALL auxbas_pw_pool%create_pw(fxc_tau_lo(ispin))
765 block
766 TYPE(pw_c1d_gs_type) :: fxc_g, fxc_g_uf
767 CALL auxbas_pw_pool%create_pw(fxc_g)
768 CALL xc_pw_pool%create_pw(fxc_g_uf)
769 CALL pw_transfer(fxc_tau_uf(ispin), fxc_g_uf)
770 CALL pw_transfer(fxc_g_uf, fxc_g)
771 CALL pw_transfer(fxc_g, fxc_tau_lo(ispin))
772 CALL xc_pw_pool%give_back_pw(fxc_g_uf)
773 CALL auxbas_pw_pool%give_back_pw(fxc_g)
774 END block
775 CALL xc_pw_pool%give_back_pw(fxc_tau_uf(ispin))
776 END DO
777 END IF
778
779 ELSE
780 CALL qs_fxc_eval(xc_deriv_set, xc_rho_set, rho1_struct, xc_section, &
781 weights, auxbas_pw_pool, fxc_rho_lo, fxc_tau_lo, &
782 is_triplet=is_triplet, spinflip=spinflip, &
783 compute_virial=do_virial, virial_xc=virial_xc)
784 END IF
785
786 ! return potentials
787 IF (ASSOCIATED(fxc_rho)) THEN
788 DO ispin = 1, min(SIZE(fxc_rho_lo), SIZE(fxc_rho))
789 CALL pw_transfer(fxc_rho_lo(ispin), fxc_rho(ispin))
790 END DO
791 DO ispin = 1, SIZE(fxc_rho_lo)
792 CALL auxbas_pw_pool%give_back_pw(fxc_rho_lo(ispin))
793 END DO
794 DEALLOCATE (fxc_rho_lo)
795 ELSE
796 fxc_rho => fxc_rho_lo
797 END IF
798 IF (ASSOCIATED(fxc_tau)) THEN
799 IF (ASSOCIATED(fxc_tau_lo)) THEN
800 DO ispin = 1, min(SIZE(fxc_tau_lo), SIZE(fxc_tau))
801 CALL pw_transfer(fxc_tau_lo(ispin), fxc_tau(ispin))
802 END DO
803 DO ispin = 1, SIZE(fxc_tau_lo)
804 CALL auxbas_pw_pool%give_back_pw(fxc_tau_lo(ispin))
805 END DO
806 DEALLOCATE (fxc_tau_lo)
807 ELSE
808 DO ispin = 1, nspins
809 CALL pw_zero(fxc_tau(ispin))
810 END DO
811 END IF
812 ELSE
813 fxc_tau => fxc_tau_lo
814 END IF
815
816 CALL timestop(handle)
817
818 END SUBROUTINE qs_fxc_apply
819
820! **************************************************************************************************
821!> \brief ...
822!> \param xc_deriv_set ...
823!> \param xc_rho_set ...
824!> \param rho1 ...
825!> \param xc_section ...
826!> \param weights ...
827!> \param auxbas_pw_pool ...
828!> \param fxc_rho ...
829!> \param fxc_tau ...
830!> \param is_triplet ...
831!> \param spinflip ...
832!> \param compute_virial ...
833!> \param virial_xc ...
834! **************************************************************************************************
835 SUBROUTINE qs_fxc_eval(xc_deriv_set, xc_rho_set, rho1, xc_section, weights, auxbas_pw_pool, &
836 fxc_rho, fxc_tau, is_triplet, spinflip, &
837 compute_virial, virial_xc)
838
839 TYPE(xc_derivative_set_type) :: xc_deriv_set
840 TYPE(xc_rho_set_type) :: xc_rho_set
841 TYPE(qs_rho_type), POINTER :: rho1
842 TYPE(section_vals_type), POINTER :: xc_section
843 TYPE(pw_r3d_rs_type), POINTER :: weights
844 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
845 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau
846 LOGICAL, INTENT(IN), OPTIONAL :: is_triplet, spinflip, compute_virial
847 REAL(kind=dp), DIMENSION(3, 3), INTENT(INOUT), &
848 OPTIONAL :: virial_xc
849
850 CHARACTER(len=*), PARAMETER :: routinen = 'qs_fxc_eval'
851
852 INTEGER :: handle, ispin, mspins, nspins
853 INTEGER, DIMENSION(2, 3) :: bo
854 LOGICAL :: do_analytic, do_sf, do_triplet, &
855 do_virial, lsd
856 REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: vxg
857 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho1_g
858 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r, tau1_r
859 TYPE(qs_rho_type), POINTER :: rhot1
860 TYPE(section_vals_type), POINTER :: xc_fun_section
861 TYPE(xc_rho_cflags_type) :: needs
862 TYPE(xc_rho_set_type) :: rho1_set
863
864 CALL timeset(routinen, handle)
865
866 do_triplet = .false.
867 IF (PRESENT(is_triplet)) do_triplet = is_triplet
868
869 do_sf = .false.
870 IF (PRESENT(spinflip)) do_sf = spinflip
871
872 do_virial = .false.
873 IF (PRESENT(compute_virial)) do_virial = compute_virial
874
875 do_analytic = section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL")
876
877 CALL qs_rho_get(rho1, rho_r=rho1_r, tau_r=tau1_r)
878 NULLIFY (rho1_g)
879
880 mspins = SIZE(rho1_r)
881 nspins = SIZE(rho1_r)
882 lsd = (nspins == 2)
883 IF (nspins == 1 .AND. do_triplet) THEN
884 nspins = 2
885 lsd = .true.
886 ELSE IF (do_sf) THEN
887 nspins = 1
888 mspins = 1
889 lsd = .true.
890 END IF
891
892 cpassert(.NOT. ASSOCIATED(fxc_rho))
893 cpassert(.NOT. ASSOCIATED(fxc_tau))
894 xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
895 needs = xc_functionals_get_needs(xc_fun_section, lsd, .true.)
896 ALLOCATE (fxc_rho(mspins))
897 DO ispin = 1, mspins
898 CALL auxbas_pw_pool%create_pw(fxc_rho(ispin))
899 CALL pw_zero(fxc_rho(ispin))
900 END DO
901 IF (needs%tau .OR. needs%tau_spin) THEN
902 IF (.NOT. ASSOCIATED(tau1_r)) THEN
903 cpabort("Tau-dependent functionals requires allocated kinetic energy density grid")
904 END IF
905 ALLOCATE (fxc_tau(mspins))
906 DO ispin = 1, mspins
907 CALL auxbas_pw_pool%create_pw(fxc_tau(ispin))
908 CALL pw_zero(fxc_tau(ispin))
909 END DO
910 END IF
911
912 IF (mspins == 1 .AND. do_triplet) THEN
913 ! split the density and response density arrays for triplet calculation
914 ALLOCATE (rhot1)
915 CALL qs_rho_create(rhot1)
916 CALL qs_rho_copy(rho1, rhot1, auxbas_pw_pool, 2, factor=2.0_dp)
917 !
918 CALL qs_rho_get(rhot1, rho_r=rho1_r, tau_r=tau1_r)
919 END IF
920
921 bo = rho1_r(1)%pw_grid%bounds_local
922 ! create the place where to store the argument for the functionals
923 CALL xc_rho_set_create(rho1_set, bo, &
924 rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
925 drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
926 tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
927
928 ! calculate the arguments needed by the functionals
929 CALL xc_rho_set_update(rho1_set, rho1_r, rho1_g, tau1_r, needs, &
930 section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
931 section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
932 auxbas_pw_pool, spinflip=do_sf)
933
934 IF (mspins == 1 .AND. do_triplet .AND. do_analytic) THEN
935
936 CALL xc_calc_2nd_deriv_analytical(fxc_rho, fxc_tau, xc_deriv_set, xc_rho_set, &
937 rho1_set, auxbas_pw_pool, xc_section, &
938 gapw=.false., vxg=vxg, tddfpt_fac=-1.0_dp, spinflip=do_sf, &
939 compute_virial=compute_virial, virial_xc=virial_xc)
940
941 ELSE IF (do_analytic) THEN
942
943 CALL xc_calc_2nd_deriv_analytical(fxc_rho, fxc_tau, xc_deriv_set, xc_rho_set, &
944 rho1_set, auxbas_pw_pool, xc_section, &
945 gapw=.false., vxg=vxg, spinflip=do_sf, &
946 compute_virial=compute_virial, virial_xc=virial_xc)
947
948 ELSE
949
950 CALL xc_calc_2nd_deriv_numerical(fxc_rho, fxc_tau, xc_rho_set, rho1_r, rho1_g, tau1_r, &
951 auxbas_pw_pool, weights, xc_section, &
952 do_triplet, compute_virial, virial_xc, xc_deriv_set)
953
954 END IF
955
956 IF (mspins == 1 .AND. do_triplet) THEN
957 CALL qs_rho_release(rhot1)
958 DEALLOCATE (rhot1)
959 END IF
960 CALL xc_rho_set_release(rho1_set)
961
962 CALL timestop(handle)
963
964 END SUBROUTINE qs_fxc_eval
965
966! **************************************************************************************************
967!> \brief ...
968!> \param qs_env ...
969!> \param rho0_struct ...
970!> \param rho1_struct ...
971!> \param xc_section ...
972!> \param accuracy ...
973!> \param fxc_rho ...
974!> \param fxc_tau ...
975!> \param is_triplet ...
976!> \param spinflip ...
977! **************************************************************************************************
978 SUBROUTINE qs_fxc_fdiff(qs_env, rho0_struct, rho1_struct, xc_section, accuracy, &
979 fxc_rho, fxc_tau, is_triplet, spinflip)
980
981 TYPE(qs_environment_type), POINTER :: qs_env
982 TYPE(qs_rho_type), POINTER :: rho0_struct, rho1_struct
983 TYPE(section_vals_type), POINTER :: xc_section
984 INTEGER, INTENT(IN) :: accuracy
985 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau
986 LOGICAL, INTENT(IN), OPTIONAL :: is_triplet, spinflip
987
988 CHARACTER(len=*), PARAMETER :: routinen = 'qs_fxc_fdiff'
989 REAL(kind=dp), PARAMETER :: epsrho = 5.e-4_dp
990
991 INTEGER :: handle, ispin, istep, nspins, nstep
992 LOGICAL :: do_sf, do_triplet
993 REAL(kind=dp) :: alpha, beta, exc, oeps1
994 REAL(kind=dp), DIMENSION(-4:4) :: ak
995 TYPE(dft_control_type), POINTER :: dft_control
996 TYPE(pw_env_type), POINTER :: pw_env
997 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
998 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_tau_rspace, vxc00
999 TYPE(qs_ks_env_type), POINTER :: ks_env
1000 TYPE(qs_rho_type), POINTER :: rhoin
1001
1002 CALL timeset(routinen, handle)
1003
1004 cpassert(.NOT. ASSOCIATED(fxc_rho))
1005 cpassert(.NOT. ASSOCIATED(fxc_tau))
1006 cpassert(ASSOCIATED(rho0_struct))
1007 cpassert(ASSOCIATED(rho1_struct))
1008
1009 do_triplet = .false.
1010 IF (PRESENT(is_triplet)) do_triplet = is_triplet
1011
1012 do_sf = .false.
1013 IF (PRESENT(spinflip)) do_sf = spinflip
1014 IF (do_sf) THEN
1015 cpabort("Spin Flip TDDFT only available with analytic 2nd xc derivatives")
1016 END IF
1017
1018 ak = 0.0_dp
1019 SELECT CASE (accuracy)
1020 CASE (:4)
1021 nstep = 2
1022 ak(-2:2) = [1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp]/12.0_dp
1023 CASE (5:7)
1024 nstep = 3
1025 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
1026 CASE (8:)
1027 nstep = 4
1028 ak(-4:4) = [1.0_dp, -32.0_dp/3.0_dp, 56.0_dp, -224.0_dp, 0.0_dp, &
1029 224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp]/280.0_dp
1030 END SELECT
1031
1032 CALL get_qs_env(qs_env, ks_env=ks_env, dft_control=dft_control, pw_env=pw_env)
1033 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1034
1035 nspins = dft_control%nspins
1036 exc = 0.0_dp
1037
1038 DO istep = -nstep, nstep
1039
1040 IF (ak(istep) /= 0.0_dp) THEN
1041 alpha = 1.0_dp
1042 beta = real(istep, kind=dp)*epsrho
1043 NULLIFY (rhoin)
1044 ALLOCATE (rhoin)
1045 CALL qs_rho_create(rhoin)
1046 NULLIFY (vxc00, v_tau_rspace)
1047 IF (do_triplet) THEN
1048 cpassert(nspins == 1)
1049 ! rhoin = (0.5 rho0, 0.5 rho0)
1050 CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, 2)
1051 ! rhoin = (0.5 rho0 + 0.5 rho1, 0.5 rho0)
1052 CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, 0.5_dp*beta)
1053 CALL qs_vxc_create(ks_env, rhoin, xc_section, vxc00, v_tau_rspace, exc)
1054 CALL pw_axpy(vxc00(2), vxc00(1), -1.0_dp)
1055 IF (ASSOCIATED(v_tau_rspace)) CALL pw_axpy(v_tau_rspace(2), v_tau_rspace(1), -1.0_dp)
1056 ELSE
1057 CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
1058 CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, beta)
1059 CALL qs_vxc_create(ks_env, rhoin, xc_section, vxc00, v_tau_rspace, exc)
1060 END IF
1061 CALL qs_rho_release(rhoin)
1062 DEALLOCATE (rhoin)
1063 IF (.NOT. ASSOCIATED(fxc_rho)) THEN
1064 ALLOCATE (fxc_rho(nspins))
1065 DO ispin = 1, nspins
1066 CALL auxbas_pw_pool%create_pw(fxc_rho(ispin))
1067 CALL pw_zero(fxc_rho(ispin))
1068 END DO
1069 END IF
1070 DO ispin = 1, nspins
1071 CALL pw_axpy(vxc00(ispin), fxc_rho(ispin), ak(istep))
1072 END DO
1073 DO ispin = 1, SIZE(vxc00)
1074 CALL auxbas_pw_pool%give_back_pw(vxc00(ispin))
1075 END DO
1076 DEALLOCATE (vxc00)
1077 IF (ASSOCIATED(v_tau_rspace)) THEN
1078 IF (.NOT. ASSOCIATED(fxc_tau)) THEN
1079 ALLOCATE (fxc_tau(nspins))
1080 DO ispin = 1, nspins
1081 CALL auxbas_pw_pool%create_pw(fxc_tau(ispin))
1082 CALL pw_zero(fxc_tau(ispin))
1083 END DO
1084 END IF
1085 DO ispin = 1, nspins
1086 CALL pw_axpy(v_tau_rspace(ispin), fxc_tau(ispin), ak(istep))
1087 END DO
1088 DO ispin = 1, SIZE(v_tau_rspace)
1089 CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
1090 END DO
1091 DEALLOCATE (v_tau_rspace)
1092 END IF
1093 END IF
1094
1095 END DO
1096
1097 oeps1 = 1.0_dp/epsrho
1098 DO ispin = 1, nspins
1099 CALL pw_scale(fxc_rho(ispin), oeps1)
1100 END DO
1101 IF (ASSOCIATED(fxc_tau)) THEN
1102 DO ispin = 1, nspins
1103 CALL pw_scale(fxc_tau(ispin), oeps1)
1104 END DO
1105 END IF
1106
1107 CALL timestop(handle)
1108
1109 END SUBROUTINE qs_fxc_fdiff
1110
1111END MODULE qs_fxc
pw_methods::pw_axpy
Definition pw_methods.F:165
pw_methods::pw_scale
Definition pw_methods.F:94
pw_methods::pw_transfer
Definition pw_methods.F:304
pw_methods::pw_zero
Definition pw_methods.F:83
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_get_rval
real(kind=dp) function, public section_get_rval(section_vals, keyword_name)
...
Definition input_section_types.F:953
input_section_types::section_get_ival
integer function, public section_get_ival(section_vals, keyword_name)
...
Definition input_section_types.F:985
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
input_section_types::section_get_lval
logical function, public section_get_lval(section_vals, keyword_name)
...
Definition input_section_types.F:1017
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition pw_env_types.F:113
pw_grids
This module defines the grid data type and some basic operations on it.
Definition pw_grids.F:36
pw_grids::pw_grid_compare
logical function, public pw_grid_compare(grida, gridb)
Check if two pw_grids are equal.
Definition pw_grids.F:148
pw_methods
Definition pw_methods.F:25
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_types
Definition pw_types.F:24
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:534
qs_fxc
Setup Routine for Fxc Potentials.
Definition qs_fxc.F:29
qs_fxc::qs_fxc_create
subroutine, public qs_fxc_create(qs_env, rho0_struct, rho1_struct, xc_section, fxc_rho, fxc_tau, is_triplet, spinflip, no_weights, uf_grid_results, pw_env_ext, compute_virial, virial_xc)
...
Definition qs_fxc.F:105
qs_fxc::qs_fxc_prep
subroutine, public qs_fxc_prep(qs_env, rho0_struct, xc_rho_set, xc_deriv_set, xc_section, pw_env_ext, is_triplet)
...
Definition qs_fxc.F:492
qs_fxc::qs_fxc_fdiff
subroutine, public qs_fxc_fdiff(qs_env, rho0_struct, rho1_struct, xc_section, accuracy, fxc_rho, fxc_tau, is_triplet, spinflip)
...
Definition qs_fxc.F:980
qs_fxc::qs_fxc_apply
subroutine, public qs_fxc_apply(qs_env, xc_deriv_set, xc_rho_set, rho1_struct, xc_section, fxc_rho, fxc_tau, is_triplet, spinflip, pw_env_ext, compute_virial, virial_xc)
...
Definition qs_fxc.F:659
qs_ks_types
Definition qs_ks_types.F:15
qs_rho_methods
methods of the rho structure (defined in qs_rho_types)
Definition qs_rho_methods.F:15
qs_rho_methods::qs_rho_copy
subroutine, public qs_rho_copy(rho_input, rho_output, auxbas_pw_pool, mspin, factor)
Allocate a density structure and fill it with data from an input structure SIZE(rho_input) == mspin =...
Definition qs_rho_methods.F:779
qs_rho_methods::qs_rho_scale_and_add
subroutine, public qs_rho_scale_and_add(rhoa, rhob, alpha, beta)
rhoa = alpha*rhoa+beta*rhob
Definition qs_rho_methods.F:1445
qs_rho_methods::qs_rho_transfer
subroutine, public qs_rho_transfer(rho_input, rho_output, in_pw_pool, out_pw_pool)
Allocate a density structure and fill it with data from an input structure Transfer all data to input...
Definition qs_rho_methods.F:1078
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
qs_rho_types::qs_rho_create
subroutine, public qs_rho_create(rho)
Allocates a new instance of rho.
Definition qs_rho_types.F:101
qs_rho_types::qs_rho_release
subroutine, public qs_rho_release(rho_struct)
releases a rho_struct by decreasing the reference count by one and deallocating if it reaches 0 (to b...
Definition qs_rho_types.F:115
qs_vxc
Definition qs_vxc.F:16
qs_vxc::qs_vxc_create
subroutine, public qs_vxc_create(ks_env, rho_struct, xc_section, vxc_rho, vxc_tau, exc, just_energy, edisp, dispersion_env, adiabatic_rescale_factor, pw_env_external, native_skala_atom_force, qs_env_external, native_gapw_composite_override, native_skala_defer_to_atom_composite)
calculates and allocates the xc potential, already reducing it to the dependence on rho and the one o...
Definition qs_vxc.F:118
xc_derivative_set_types
represent a group ofunctional derivatives
Definition xc_derivative_set_types.F:14
xc_derivative_set_types::xc_dset_release
subroutine, public xc_dset_release(derivative_set)
releases a derivative set
Definition xc_derivative_set_types.F:143
xc_derivatives
Definition xc_derivatives.F:9
xc_derivatives::xc_functionals_get_needs
type(xc_rho_cflags_type) function, public xc_functionals_get_needs(functionals, lsd, calc_potential)
...
Definition xc_derivatives.F:770
xc_rho_cflags_types
contains the structure
Definition xc_rho_cflags_types.F:14
xc_rho_set_types
contains the structure
Definition xc_rho_set_types.F:14
xc_rho_set_types::xc_rho_set_create
subroutine, public xc_rho_set_create(rho_set, local_bounds, rho_cutoff, drho_cutoff, tau_cutoff)
allocates and does (minimal) initialization of a rho_set
Definition xc_rho_set_types.F:111
xc_rho_set_types::xc_rho_set_release
subroutine, public xc_rho_set_release(rho_set, pw_pool)
releases the given rho_set
Definition xc_rho_set_types.F:129
xc_rho_set_types::xc_rho_set_update
subroutine, public xc_rho_set_update(rho_set, rho_r, rho_g, tau, needs, xc_deriv_method_id, xc_rho_smooth_id, pw_pool, spinflip)
updates the given rho set with the density given by rho_r (and rho_g). The rho set will contain the c...
Definition xc_rho_set_types.F:691
xc
Exchange and Correlation functional calculations.
Definition xc.F:17
xc::xc_prep_2nd_deriv
subroutine, public xc_prep_2nd_deriv(deriv_set, rho_set, rho_r, pw_pool, weights, xc_section, tau_r)
Prepare objects for the calculation of the 2nd derivatives of the density functional....
Definition xc.F:5539
xc::xc_calc_2nd_deriv_analytical
subroutine, public xc_calc_2nd_deriv_analytical(v_xc, v_xc_tau, deriv_set, rho_set, rho1_set, pw_pool, xc_section, gapw, vxg, tddfpt_fac, compute_virial, virial_xc, spinflip)
Calculates the second derivative of E_xc at rho in the direction rho1 (if you see the second derivati...
Definition xc.F:2056
xc::xc_calc_2nd_deriv_numerical
subroutine, public xc_calc_2nd_deriv_numerical(v_xc, v_tau, rho_set, rho1_r, rho1_g, tau1_r, pw_pool, weights, xc_section, do_triplet, calc_virial, virial_xc, deriv_set)
calculates 2nd derivative numerically
Definition xc.F:1062
cp_control_types::dft_control_type
Definition cp_control_types.F:757
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:137
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
xc_derivative_set_types::xc_derivative_set_type
A derivative set contains the different derivatives of a xc-functional in form of a linked list.
Definition xc_derivative_set_types.F:48
xc_rho_cflags_types::xc_rho_cflags_type
contains a flag for each component of xc_rho_set, so that you can use it to tell which components you...
Definition xc_rho_cflags_types.F:48
xc_rho_set_types::xc_rho_set_type
represent a density, with all the representation and data needed to perform a functional evaluation
Definition xc_rho_set_types.F:78