(git:039bc03)
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 https://en.wikipedia.org/wiki/Finite_difference_coefficient
10!---------------------------------------------------------------------------------------------------
11!Derivative Accuracy 4 3 2 1 0 1 2 3 4
12!---------------------------------------------------------------------------------------------------
13! 1 2 -1/2 0 1/2
14! 4 1/12 -2/3 0 2/3 -1/12
15! 6 -1/60 3/20 -3/4 0 3/4 -3/20 1/60
16! 8 1/280 -4/105 1/5 -4/5 0 4/5 -1/5 4/105 -1/280
17!---------------------------------------------------------------------------------------------------
18! 2 2 1 -2 1
19! 4 -1/12 4/3 -5/2 4/3 -1/12
20! 6 1/90 -3/20 3/2 -49/18 3/2 -3/20 1/90
21! 8 -1/560 8/315 -1/5 8/5 -205/72 8/5 -1/5 8/315 -1/560
22!---------------------------------------------------------------------------------------------------
23!> \par History
24!> init 17.03.2020
25!> \author JGH
26! **************************************************************************************************
27MODULE qs_fxc
28
34 USE kinds, ONLY: dp
35 USE pw_env_types, ONLY: pw_env_get,&
37 USE pw_methods, ONLY: pw_axpy,&
38 pw_scale,&
41 USE pw_types, ONLY: pw_c1d_gs_type,&
43 USE qs_ks_types, ONLY: get_ks_env,&
45 USE qs_rho_methods, ONLY: qs_rho_copy,&
48 USE qs_rho_types, ONLY: qs_rho_create,&
52 USE qs_vxc, ONLY: qs_vxc_create
53 USE xc, ONLY: xc_calc_2nd_deriv,&
66#include "./base/base_uses.f90"
67
68 IMPLICIT NONE
69
70 PRIVATE
71
72 ! *** Public subroutines ***
74
75 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_fxc'
76
77! **************************************************************************************************
78
79CONTAINS
80
81! **************************************************************************************************
82!> \brief ...
83!> \param rho0 ...
84!> \param rho1_r ...
85!> \param tau1_r ...
86!> \param xc_section ...
87!> \param weights ...
88!> \param auxbas_pw_pool ...
89!> \param is_triplet ...
90!> \param v_xc ...
91!> \param v_xc_tau ...
92!> \param spinflip ...
93! **************************************************************************************************
94 SUBROUTINE qs_fxc_analytic(rho0, rho1_r, tau1_r, xc_section, weights, auxbas_pw_pool, &
95 is_triplet, v_xc, v_xc_tau, spinflip)
96
97 TYPE(qs_rho_type), POINTER :: rho0
98 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r, tau1_r
99 TYPE(section_vals_type), POINTER :: xc_section
100 TYPE(pw_r3d_rs_type), POINTER :: weights
101 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
102 LOGICAL, INTENT(IN) :: is_triplet
103 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_xc, v_xc_tau
104 LOGICAL, OPTIONAL :: spinflip
105
106 CHARACTER(len=*), PARAMETER :: routinen = 'qs_fxc_analytic'
107
108 INTEGER :: handle, nspins
109 INTEGER, DIMENSION(2, 3) :: bo
110 LOGICAL :: do_sf, lsd
111 REAL(kind=dp) :: fac
112 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho0_g, rho1_g
113 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho0_r, tau0_r
114 TYPE(section_vals_type), POINTER :: xc_fun_section
115 TYPE(xc_derivative_set_type) :: deriv_set
116 TYPE(xc_rho_cflags_type) :: needs
117 TYPE(xc_rho_set_type) :: rho0_set
118
119 CALL timeset(routinen, handle)
120
121 cpassert(.NOT. ASSOCIATED(v_xc))
122 cpassert(.NOT. ASSOCIATED(v_xc_tau))
123
124 do_sf = .false.
125 IF (PRESENT(spinflip)) do_sf = spinflip
126
127 CALL qs_rho_get(rho0, rho_r=rho0_r, rho_g=rho0_g, tau_r=tau0_r)
128 nspins = SIZE(rho0_r)
129
130 lsd = (nspins == 2)
131 fac = 0._dp
132 IF (is_triplet .AND. nspins == 1) fac = -1.0_dp
133
134 NULLIFY (rho1_g)
135 bo = rho1_r(1)%pw_grid%bounds_local
136 xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
137 needs = xc_functionals_get_needs(xc_fun_section, lsd, .true.)
138 ! calculate the arguments needed by the functionals and the values of the functional on the grid
139 CALL xc_prep_2nd_deriv(deriv_set, rho0_set, rho0_r, auxbas_pw_pool, weights, &
140 xc_section=xc_section, tau_r=tau0_r)
141 ! Folds the density rho1 with the functional
142 CALL xc_calc_2nd_deriv(v_xc, v_xc_tau, deriv_set, rho0_set, rho1_r, rho1_g, tau1_r, &
143 auxbas_pw_pool, weights, xc_section=xc_section, &
144 gapw=.false., do_triplet=is_triplet, do_sf=do_sf)
145 CALL xc_dset_release(deriv_set)
146 CALL xc_rho_set_release(rho0_set)
147
148 CALL timestop(handle)
149
150 END SUBROUTINE qs_fxc_analytic
151
152! **************************************************************************************************
153!> \brief ...
154!> \param ks_env ...
155!> \param rho0_struct ...
156!> \param rho1_struct ...
157!> \param xc_section ...
158!> \param accuracy ...
159!> \param is_triplet ...
160!> \param fxc_rho ...
161!> \param fxc_tau ...
162! **************************************************************************************************
163 SUBROUTINE qs_fxc_fdiff(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, &
164 is_triplet, fxc_rho, fxc_tau)
165
166 TYPE(qs_ks_env_type), POINTER :: ks_env
167 TYPE(qs_rho_type), POINTER :: rho0_struct, rho1_struct
168 TYPE(section_vals_type), POINTER :: xc_section
169 INTEGER, INTENT(IN) :: accuracy
170 LOGICAL, INTENT(IN) :: is_triplet
171 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau
172
173 CHARACTER(len=*), PARAMETER :: routinen = 'qs_fxc_fdiff'
174 REAL(kind=dp), PARAMETER :: epsrho = 5.e-4_dp
175
176 INTEGER :: handle, ispin, istep, nspins, nstep
177 REAL(kind=dp) :: alpha, beta, exc, oeps1
178 REAL(kind=dp), DIMENSION(-4:4) :: ak
179 TYPE(dft_control_type), POINTER :: dft_control
180 TYPE(pw_env_type), POINTER :: pw_env
181 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
182 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_tau_rspace, vxc00
183 TYPE(qs_rho_type), POINTER :: rhoin
184
185 CALL timeset(routinen, handle)
186
187 cpassert(.NOT. ASSOCIATED(fxc_rho))
188 cpassert(.NOT. ASSOCIATED(fxc_tau))
189 cpassert(ASSOCIATED(rho0_struct))
190 cpassert(ASSOCIATED(rho1_struct))
191
192 ak = 0.0_dp
193 SELECT CASE (accuracy)
194 CASE (:4)
195 nstep = 2
196 ak(-2:2) = [1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp]/12.0_dp
197 CASE (5:7)
198 nstep = 3
199 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
200 CASE (8:)
201 nstep = 4
202 ak(-4:4) = [1.0_dp, -32.0_dp/3.0_dp, 56.0_dp, -224.0_dp, 0.0_dp, &
203 224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp]/280.0_dp
204 END SELECT
205
206 CALL get_ks_env(ks_env, dft_control=dft_control, pw_env=pw_env)
207 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
208
209 nspins = dft_control%nspins
210 exc = 0.0_dp
211
212 DO istep = -nstep, nstep
213
214 IF (ak(istep) /= 0.0_dp) THEN
215 alpha = 1.0_dp
216 beta = real(istep, kind=dp)*epsrho
217 NULLIFY (rhoin)
218 ALLOCATE (rhoin)
219 CALL qs_rho_create(rhoin)
220 NULLIFY (vxc00, v_tau_rspace)
221 IF (is_triplet) THEN
222 cpassert(nspins == 1)
223 ! rhoin = (0.5 rho0, 0.5 rho0)
224 CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, 2)
225 ! rhoin = (0.5 rho0 + 0.5 rho1, 0.5 rho0)
226 CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, 0.5_dp*beta)
227 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rhoin, xc_section=xc_section, &
228 vxc_rho=vxc00, vxc_tau=v_tau_rspace, exc=exc, just_energy=.false.)
229 CALL pw_axpy(vxc00(2), vxc00(1), -1.0_dp)
230 IF (ASSOCIATED(v_tau_rspace)) CALL pw_axpy(v_tau_rspace(2), v_tau_rspace(1), -1.0_dp)
231 ELSE
232 CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
233 CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, beta)
234 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rhoin, xc_section=xc_section, &
235 vxc_rho=vxc00, vxc_tau=v_tau_rspace, exc=exc, just_energy=.false.)
236 END IF
237 CALL qs_rho_release(rhoin)
238 DEALLOCATE (rhoin)
239 IF (.NOT. ASSOCIATED(fxc_rho)) THEN
240 ALLOCATE (fxc_rho(nspins))
241 DO ispin = 1, nspins
242 CALL auxbas_pw_pool%create_pw(fxc_rho(ispin))
243 CALL pw_zero(fxc_rho(ispin))
244 END DO
245 END IF
246 DO ispin = 1, nspins
247 CALL pw_axpy(vxc00(ispin), fxc_rho(ispin), ak(istep))
248 END DO
249 DO ispin = 1, SIZE(vxc00)
250 CALL auxbas_pw_pool%give_back_pw(vxc00(ispin))
251 END DO
252 DEALLOCATE (vxc00)
253 IF (ASSOCIATED(v_tau_rspace)) THEN
254 IF (.NOT. ASSOCIATED(fxc_tau)) THEN
255 ALLOCATE (fxc_tau(nspins))
256 DO ispin = 1, nspins
257 CALL auxbas_pw_pool%create_pw(fxc_tau(ispin))
258 CALL pw_zero(fxc_tau(ispin))
259 END DO
260 END IF
261 DO ispin = 1, nspins
262 CALL pw_axpy(v_tau_rspace(ispin), fxc_tau(ispin), ak(istep))
263 END DO
264 DO ispin = 1, SIZE(v_tau_rspace)
265 CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
266 END DO
267 DEALLOCATE (v_tau_rspace)
268 END IF
269 END IF
270
271 END DO
272
273 oeps1 = 1.0_dp/epsrho
274 DO ispin = 1, nspins
275 CALL pw_scale(fxc_rho(ispin), oeps1)
276 END DO
277 IF (ASSOCIATED(fxc_tau)) THEN
278 DO ispin = 1, nspins
279 CALL pw_scale(fxc_tau(ispin), oeps1)
280 END DO
281 END IF
282
283 CALL timestop(handle)
284
285 END SUBROUTINE qs_fxc_fdiff
286
287! **************************************************************************************************
288!> \brief Calculates the values at the grid points in real space (r_i), of the second and third
289!> functional derivatives of the exchange-correlation energy functional.
290!> fxc_rho(r_i) = fxc[n](r_i)*n^(1)(r_i) ! Second functional derivative
291!> gxc_rho(r_i) = n^(1)(r_i)*gxc[n](r_i)*n^(1)(r_i) ! Third functional derivative
292!> \param rho0_struct Ground state density, n(r).
293!> \param rho1_struct Density used to fold the functional derivatives, n^(1)(r).
294!> \param xc_section ...
295!> \param weights ...
296!> \param pw_pool ...
297!> \param fxc_rho Second functional derivative with respect to the density, n(r).
298!> \param fxc_tau mGGA contribution to the second functional derivative with respect to the density.
299!> \param gxc_rho Third functional derivative with respect to the density, n(r).
300!> \param gxc_tau mGGA contribution to the third functional derivative with respect to the density.
301!> \param spinflip Flag used to activate the spin-flip noncollinear kernel and kernel derivatives.
302!> \par History
303!> * 07.2024 Created [LHS]
304! **************************************************************************************************
305 SUBROUTINE qs_fgxc_analytic(rho0_struct, rho1_struct, xc_section, weights, pw_pool, &
306 fxc_rho, fxc_tau, gxc_rho, gxc_tau, spinflip)
307
308 TYPE(qs_rho_type), POINTER :: rho0_struct, rho1_struct
309 TYPE(section_vals_type), POINTER :: xc_section
310 TYPE(pw_r3d_rs_type), POINTER :: weights
311 TYPE(pw_pool_type), POINTER :: pw_pool
312 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau, gxc_rho, gxc_tau
313 LOGICAL, INTENT(IN), OPTIONAL :: spinflip
314
315 CHARACTER(len=*), PARAMETER :: routinen = 'qs_fgxc_analytic'
316
317 INTEGER :: handle, ispin, nspins, spindim
318 INTEGER, DIMENSION(2, 3) :: bo
319 LOGICAL :: do_sf, lsd
320 REAL(kind=dp) :: fac
321 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho0_g, rho1_g
322 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho0_r, rho1_r, tau0_r, tau1_r
323 TYPE(section_vals_type), POINTER :: xc_fun_section
324 TYPE(xc_derivative_set_type) :: deriv_set
325 TYPE(xc_rho_cflags_type) :: needs
326 TYPE(xc_rho_set_type) :: rho0_set, rho1_set
327
328 CALL timeset(routinen, handle)
329
330 ! Only rho0 and rho1 should be associated
331 cpassert(.NOT. ASSOCIATED(fxc_rho))
332 cpassert(.NOT. ASSOCIATED(fxc_tau))
333 cpassert(.NOT. ASSOCIATED(gxc_rho))
334 cpassert(.NOT. ASSOCIATED(gxc_tau))
335 cpassert(ASSOCIATED(rho0_struct))
336 cpassert(ASSOCIATED(rho1_struct))
337
338 ! Initialize parameters
339 do_sf = .false.
340 IF (PRESENT(spinflip)) do_sf = spinflip
341 !
342 ! Get the values on the gridpoints of the rho0 density
343 CALL qs_rho_get(rho0_struct, rho_r=rho0_r, rho_g=rho0_g, tau_r=tau0_r)
344 nspins = SIZE(rho0_r)
345 lsd = (nspins == 2)
346 !
347 IF (do_sf) THEN
348 spindim = 1
349 ELSE
350 spindim = nspins
351 END IF
352 !
353 fac = 0._dp
354 IF (nspins == 1) THEN
355 fac = 1.0_dp
356 END IF
357
358 ! Read xc functional section and find out what the functional actually needs
359 xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
360 needs = xc_functionals_get_needs(xc_fun_section, lsd, .true.)
361
362 ! Create fields for the kernel and kernel derivative
363 ALLOCATE (fxc_rho(spindim), gxc_rho(nspins))
364 DO ispin = 1, spindim
365 CALL pw_pool%create_pw(fxc_rho(ispin))
366 CALL pw_zero(fxc_rho(ispin))
367 END DO
368 DO ispin = 1, nspins
369 CALL pw_pool%create_pw(gxc_rho(ispin))
370 CALL pw_zero(gxc_rho(ispin))
371 END DO
372 ! Create fields for mGGA functionals. This implementation is not ready yet!
373 IF (needs%tau .OR. needs%tau_spin) THEN
374 IF (.NOT. ASSOCIATED(tau1_r)) &
375 cpabort("Tau-dependent functionals requires allocated kinetic energy density grid")
376 ALLOCATE (fxc_tau(spindim), gxc_tau(nspins))
377 DO ispin = 1, spindim
378 CALL pw_pool%create_pw(fxc_tau(ispin))
379 CALL pw_zero(fxc_tau(ispin))
380 END DO
381 DO ispin = 1, nspins
382 CALL pw_pool%create_pw(gxc_tau(ispin))
383 CALL pw_zero(gxc_tau(ispin))
384 END DO
385 END IF
386
387 ! Build rho0_set
388 ! calculate the arguments needed by the functionals
389 ! Needs
390 ! deriv_set xc_derivative_set_type just declared
391 ! rho0_set xc_rho_set_type just declared
392 ! rho0_r pw_type calculated by qs_rho_get
393 ! pw_pool given by the calling subroutine
394 ! xc_section given by the calling subroutine
395 ! tau0_r pw_type calculated by qs_rho_get
396 CALL xc_prep_3rd_deriv(deriv_set, rho0_set, rho0_r, pw_pool, weights, &
397 xc_section, tau_r=tau0_r, do_sf=do_sf)
398
399 ! Build rho1_set
400 ! Get the values on the gridpoints of the rho1 density
401 CALL qs_rho_get(rho1_struct, rho_r=rho1_r, rho_g=rho1_g, tau_r=tau1_r)
402 bo = rho1_r(1)%pw_grid%bounds_local
403 ! create the place where to store the argument for the functionals
404 ! Needs
405 ! rho1_set xc_rho_set_type just declared
406 ! bo 2x3 integer matrix should have bounds_local or rho1_r
407 CALL xc_rho_set_create(rho1_set, bo, &
408 rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
409 drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
410 tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
411 ! calculate the arguments needed by the functionals
412 ! Needs
413 ! rho1_set object created by xc_rho_set_create
414 ! rho1_r,rho1_g,tau1_r pw_type values of rho1 in real space grid
415 ! needs xc_rho_cflags_type defined through xc_functionals_get_needs
416 ! pw_pool Given by the calling subroutine
417 CALL xc_rho_set_update(rho1_set, rho1_r, rho1_g, tau1_r, needs, &
418 section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
419 section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
420 pw_pool, spinflip=do_sf)
421
422 ! Calculate exchange correlation kernel
423 ! Needs
424 ! fxc_rho, fxc_tau pw_type not associated
425 ! deriv_set created and defined by xc_prep_3rd_deriv
426 ! rho0_set xc_rho_set_type build by xc_prep_3rd_deriv
427 ! rho1_set xc_rho_set_type build by xc_rho_set_create/update
428 ! pw_pool needs to be given by the calling subroutine
429 ! xc_section needs to be given by the calling subroutine
430 CALL xc_calc_2nd_deriv_analytical(fxc_rho, fxc_tau, deriv_set, rho0_set, rho1_set, pw_pool, &
431 xc_section, .false., spinflip=do_sf, tddfpt_fac=fac)
432 ! Calculate exchange correlation kernel derivative
433 CALL xc_calc_3rd_deriv_analytical(gxc_rho, gxc_tau, deriv_set, rho0_set, rho1_set, pw_pool, &
434 xc_section, spinflip=do_sf)
435
436 CALL xc_dset_release(deriv_set)
437 CALL xc_rho_set_release(rho0_set)
438 CALL xc_rho_set_release(rho1_set)
439
440 CALL timestop(handle)
441
442 END SUBROUTINE qs_fgxc_analytic
443
444! **************************************************************************************************
445!> \brief ...
446!> \param ks_env ...
447!> \param rho0_struct ...
448!> \param rho1_struct ...
449!> \param xc_section ...
450!> \param accuracy ...
451!> \param epsrho ...
452!> \param is_triplet ...
453!> \param weights ...
454!> \param fxc_rho ...
455!> \param fxc_tau ...
456!> \param gxc_rho ...
457!> \param gxc_tau ...
458!> \param spinflip ...
459! **************************************************************************************************
460 SUBROUTINE qs_fgxc_gdiff(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, epsrho, &
461 is_triplet, weights, fxc_rho, fxc_tau, gxc_rho, gxc_tau, spinflip)
462
463 TYPE(qs_ks_env_type), POINTER :: ks_env
464 TYPE(qs_rho_type), POINTER :: rho0_struct, rho1_struct
465 TYPE(section_vals_type), POINTER :: xc_section
466 INTEGER, INTENT(IN) :: accuracy
467 REAL(kind=dp), INTENT(IN) :: epsrho
468 LOGICAL, INTENT(IN) :: is_triplet
469 TYPE(pw_r3d_rs_type), POINTER :: weights
470 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau, gxc_rho, gxc_tau
471 LOGICAL, OPTIONAL :: spinflip
472
473 CHARACTER(len=*), PARAMETER :: routinen = 'qs_fgxc_gdiff'
474
475 INTEGER :: handle, ispin, istep, nspins, nstep
476 LOGICAL :: do_sf
477 REAL(kind=dp) :: alpha, beta, exc, oeps1
478 REAL(kind=dp), DIMENSION(-4:4) :: ak
479 TYPE(dft_control_type), POINTER :: dft_control
480 TYPE(pw_env_type), POINTER :: pw_env
481 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
482 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r, tau1_r, v_tau_rspace, vxc00, &
483 vxc00b
484 TYPE(qs_rho_type), POINTER :: rhoin
485
486 CALL timeset(routinen, handle)
487
488 cpassert(.NOT. ASSOCIATED(fxc_rho))
489 cpassert(.NOT. ASSOCIATED(fxc_tau))
490 cpassert(.NOT. ASSOCIATED(gxc_rho))
491 cpassert(.NOT. ASSOCIATED(gxc_tau))
492 cpassert(ASSOCIATED(rho0_struct))
493 cpassert(ASSOCIATED(rho1_struct))
494
495 do_sf = .false.
496 IF (PRESENT(spinflip)) do_sf = spinflip
497
498 ak = 0.0_dp
499 SELECT CASE (accuracy)
500 CASE (:4)
501 nstep = 2
502 ak(-2:2) = [1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp]/12.0_dp
503 CASE (5:7)
504 nstep = 3
505 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
506 CASE (8:)
507 nstep = 4
508 ak(-4:4) = [1.0_dp, -32.0_dp/3.0_dp, 56.0_dp, -224.0_dp, 0.0_dp, &
509 224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp]/280.0_dp
510 END SELECT
511
512 CALL get_ks_env(ks_env, dft_control=dft_control, pw_env=pw_env)
513 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
514
515 nspins = dft_control%nspins
516 exc = 0.0_dp
517
518 IF (do_sf) THEN
519 CALL qs_rho_get(rho1_struct, rho_r=rho1_r, tau_r=tau1_r)
520 CALL qs_fxc_analytic(rho0_struct, rho1_r, tau1_r, xc_section, &
521 weights, auxbas_pw_pool, is_triplet, &
522 fxc_rho, fxc_tau, spinflip=do_sf)
523 ELSE
524 CALL qs_fxc_fdiff(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, is_triplet, &
525 fxc_rho, fxc_tau)
526 END IF
527
528 DO istep = -nstep, nstep
529
530 IF (ak(istep) /= 0.0_dp) THEN
531 alpha = 1.0_dp
532 beta = real(istep, kind=dp)*epsrho
533 NULLIFY (rhoin)
534 ALLOCATE (rhoin)
535 CALL qs_rho_create(rhoin)
536 NULLIFY (vxc00, vxc00b, v_tau_rspace)
537 CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
538 CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, beta)
539 IF (do_sf) THEN
540 ! variation in alpha density
541 CALL qs_fxc_analytic(rhoin, rho1_r, tau1_r, &
542 xc_section, weights, auxbas_pw_pool, is_triplet, &
543 vxc00, v_tau_rspace, spinflip=do_sf)
544 ! variation in beta density
545 CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
546 CALL qs_rho_scale_and_add_b(rhoin, rho1_struct, alpha, beta)
547 CALL qs_fxc_analytic(rhoin, rho1_r, tau1_r, &
548 xc_section, weights, auxbas_pw_pool, is_triplet, &
549 vxc00b, v_tau_rspace, spinflip=do_sf)
550 ELSE
551 CALL qs_fxc_fdiff(ks_env=ks_env, rho0_struct=rhoin, rho1_struct=rho1_struct, &
552 xc_section=xc_section, accuracy=accuracy, is_triplet=is_triplet, &
553 fxc_rho=vxc00, fxc_tau=v_tau_rspace)
554 END IF
555 CALL qs_rho_release(rhoin)
556 DEALLOCATE (rhoin)
557 IF (.NOT. ASSOCIATED(gxc_rho)) THEN
558 ALLOCATE (gxc_rho(nspins))
559 DO ispin = 1, nspins
560 CALL auxbas_pw_pool%create_pw(gxc_rho(ispin))
561 CALL pw_zero(gxc_rho(ispin))
562 END DO
563 END IF
564 IF (do_sf) THEN
565 CALL pw_axpy(vxc00(1), gxc_rho(1), ak(istep))
566 CALL pw_axpy(vxc00b(1), gxc_rho(2), ak(istep))
567 ELSE
568 DO ispin = 1, nspins
569 CALL pw_axpy(vxc00(ispin), gxc_rho(ispin), ak(istep))
570 END DO
571 END IF
572 DO ispin = 1, SIZE(vxc00)
573 CALL auxbas_pw_pool%give_back_pw(vxc00(ispin))
574 END DO
575 DEALLOCATE (vxc00)
576 IF (ASSOCIATED(vxc00b)) THEN
577 CALL auxbas_pw_pool%give_back_pw(vxc00b(1))
578 DEALLOCATE (vxc00b)
579 END IF
580 IF (ASSOCIATED(v_tau_rspace)) THEN
581 IF (.NOT. ASSOCIATED(gxc_tau)) THEN
582 ALLOCATE (gxc_tau(nspins))
583 DO ispin = 1, nspins
584 CALL auxbas_pw_pool%create_pw(gxc_tau(ispin))
585 CALL pw_zero(gxc_tau(ispin))
586 END DO
587 END IF
588 DO ispin = 1, nspins
589 CALL pw_axpy(v_tau_rspace(ispin), gxc_tau(ispin), ak(istep))
590 END DO
591 DO ispin = 1, SIZE(v_tau_rspace)
592 CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
593 END DO
594 DEALLOCATE (v_tau_rspace)
595 END IF
596 END IF
597
598 END DO
599
600 oeps1 = 1.0_dp/epsrho
601 DO ispin = 1, nspins
602 CALL pw_scale(gxc_rho(ispin), oeps1)
603 END DO
604 IF (ASSOCIATED(gxc_tau)) THEN
605 DO ispin = 1, nspins
606 CALL pw_scale(gxc_tau(ispin), oeps1)
607 END DO
608 END IF
609
610 CALL timestop(handle)
611
612 END SUBROUTINE qs_fgxc_gdiff
613
614! **************************************************************************************************
615!> \brief ...
616!> \param ks_env ...
617!> \param rho0_struct ...
618!> \param rho1_struct ...
619!> \param xc_section ...
620!> \param accuracy ...
621!> \param is_triplet ...
622!> \param fxc_rho ...
623!> \param fxc_tau ...
624!> \param gxc_rho ...
625!> \param gxc_tau ...
626! **************************************************************************************************
627 SUBROUTINE qs_fgxc_create(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, is_triplet, &
628 fxc_rho, fxc_tau, gxc_rho, gxc_tau)
629
630 TYPE(qs_ks_env_type), POINTER :: ks_env
631 TYPE(qs_rho_type), POINTER :: rho0_struct, rho1_struct
632 TYPE(section_vals_type), POINTER :: xc_section
633 INTEGER, INTENT(IN) :: accuracy
634 LOGICAL, INTENT(IN) :: is_triplet
635 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau, gxc_rho, gxc_tau
636
637 CHARACTER(len=*), PARAMETER :: routinen = 'qs_fgxc_create'
638 REAL(kind=dp), PARAMETER :: epsrho = 5.e-4_dp
639
640 INTEGER :: handle, ispin, istep, nspins, nstep
641 REAL(kind=dp) :: alpha, beta, exc, oeps1, oeps2
642 REAL(kind=dp), DIMENSION(-4:4) :: ak, bl
643 TYPE(dft_control_type), POINTER :: dft_control
644 TYPE(pw_env_type), POINTER :: pw_env
645 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
646 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_tau_rspace, vxc00
647 TYPE(qs_rho_type), POINTER :: rhoin
648
649 CALL timeset(routinen, handle)
650
651 cpassert(.NOT. ASSOCIATED(fxc_rho))
652 cpassert(.NOT. ASSOCIATED(fxc_tau))
653 cpassert(.NOT. ASSOCIATED(gxc_rho))
654 cpassert(.NOT. ASSOCIATED(gxc_tau))
655 cpassert(ASSOCIATED(rho0_struct))
656 cpassert(ASSOCIATED(rho1_struct))
657
658 ak = 0.0_dp
659 bl = 0.0_dp
660 SELECT CASE (accuracy)
661 CASE (:4)
662 nstep = 2
663 ak(-2:2) = [1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp]/12.0_dp
664 bl(-2:2) = [-1.0_dp, 16.0_dp, -30.0_dp, 16.0_dp, -1.0_dp]/12.0_dp
665 CASE (5:7)
666 nstep = 3
667 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
668 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
669 CASE (8:)
670 nstep = 4
671 ak(-4:4) = [1.0_dp, -32.0_dp/3.0_dp, 56.0_dp, -224.0_dp, 0.0_dp, &
672 224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp]/280.0_dp
673 bl(-4:4) = [-1.0_dp, 128.0_dp/9.0_dp, -112.0_dp, 896.0_dp, -14350.0_dp/9.0_dp, &
674 896.0_dp, -112.0_dp, 128.0_dp/9.0_dp, -1.0_dp]/560.0_dp
675 END SELECT
676
677 CALL get_ks_env(ks_env, dft_control=dft_control, pw_env=pw_env)
678 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
679
680 nspins = dft_control%nspins
681 exc = 0.0_dp
682
683 DO istep = -nstep, nstep
684
685 alpha = 1.0_dp
686 beta = real(istep, kind=dp)*epsrho
687 NULLIFY (rhoin)
688 ALLOCATE (rhoin)
689 CALL qs_rho_create(rhoin)
690 NULLIFY (vxc00, v_tau_rspace)
691 IF (is_triplet) THEN
692 cpassert(nspins == 1)
693 ! rhoin = (0.5 rho0, 0.5 rho0)
694 CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, 2)
695 ! rhoin = (0.5 rho0 + 0.5 rho1, 0.5 rho0)
696 CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, 0.5_dp*beta)
697 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rhoin, xc_section=xc_section, &
698 vxc_rho=vxc00, vxc_tau=v_tau_rspace, exc=exc, just_energy=.false.)
699 CALL pw_axpy(vxc00(2), vxc00(1), -1.0_dp)
700 ELSE
701 CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
702 CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, beta)
703 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rhoin, xc_section=xc_section, &
704 vxc_rho=vxc00, vxc_tau=v_tau_rspace, exc=exc, just_energy=.false.)
705 END IF
706 CALL qs_rho_release(rhoin)
707 DEALLOCATE (rhoin)
708 IF (.NOT. ASSOCIATED(fxc_rho)) THEN
709 ALLOCATE (fxc_rho(nspins))
710 DO ispin = 1, nspins
711 CALL auxbas_pw_pool%create_pw(fxc_rho(ispin))
712 CALL pw_zero(fxc_rho(ispin))
713 END DO
714 END IF
715 IF (.NOT. ASSOCIATED(gxc_rho)) THEN
716 ALLOCATE (gxc_rho(nspins))
717 DO ispin = 1, nspins
718 CALL auxbas_pw_pool%create_pw(gxc_rho(ispin))
719 CALL pw_zero(gxc_rho(ispin))
720 END DO
721 END IF
722 cpassert(.NOT. ASSOCIATED(v_tau_rspace))
723 DO ispin = 1, nspins
724 IF (ak(istep) /= 0.0_dp) THEN
725 CALL pw_axpy(vxc00(ispin), fxc_rho(ispin), ak(istep))
726 END IF
727 IF (bl(istep) /= 0.0_dp) THEN
728 CALL pw_axpy(vxc00(ispin), gxc_rho(ispin), bl(istep))
729 END IF
730 END DO
731 DO ispin = 1, SIZE(vxc00)
732 CALL auxbas_pw_pool%give_back_pw(vxc00(ispin))
733 END DO
734 DEALLOCATE (vxc00)
735
736 END DO
737
738 oeps1 = 1.0_dp/epsrho
739 oeps2 = 1.0_dp/(epsrho**2)
740 DO ispin = 1, nspins
741 CALL pw_scale(fxc_rho(ispin), oeps1)
742 CALL pw_scale(gxc_rho(ispin), oeps2)
743 END DO
744
745 CALL timestop(handle)
746
747 END SUBROUTINE qs_fgxc_create
748
749! **************************************************************************************************
750!> \brief ...
751!> \param ks_env ...
752!> \param fxc_rho ...
753!> \param fxc_tau ...
754!> \param gxc_rho ...
755!> \param gxc_tau ...
756! **************************************************************************************************
757 SUBROUTINE qs_fgxc_release(ks_env, fxc_rho, fxc_tau, gxc_rho, gxc_tau)
758
759 TYPE(qs_ks_env_type), POINTER :: ks_env
760 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau, gxc_rho, gxc_tau
761
762 INTEGER :: ispin
763 TYPE(pw_env_type), POINTER :: pw_env
764 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
765
766 CALL get_ks_env(ks_env, pw_env=pw_env)
767 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
768
769 IF (ASSOCIATED(fxc_rho)) THEN
770 DO ispin = 1, SIZE(fxc_rho)
771 CALL auxbas_pw_pool%give_back_pw(fxc_rho(ispin))
772 END DO
773 DEALLOCATE (fxc_rho)
774 END IF
775 IF (ASSOCIATED(fxc_tau)) THEN
776 DO ispin = 1, SIZE(fxc_tau)
777 CALL auxbas_pw_pool%give_back_pw(fxc_tau(ispin))
778 END DO
779 DEALLOCATE (fxc_tau)
780 END IF
781 IF (ASSOCIATED(gxc_rho)) THEN
782 DO ispin = 1, SIZE(gxc_rho)
783 CALL auxbas_pw_pool%give_back_pw(gxc_rho(ispin))
784 END DO
785 DEALLOCATE (gxc_rho)
786 END IF
787 IF (ASSOCIATED(gxc_tau)) THEN
788 DO ispin = 1, SIZE(gxc_tau)
789 CALL auxbas_pw_pool%give_back_pw(gxc_tau(ispin))
790 END DO
791 DEALLOCATE (gxc_tau)
792 END IF
793
794 END SUBROUTINE qs_fgxc_release
795
796END MODULE qs_fxc
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition grid_common.h:51
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
objects that represent the structure of input sections and the data contained in an input section
real(kind=dp) function, public section_get_rval(section_vals, keyword_name)
...
integer function, public section_get_ival(section_vals, keyword_name)
...
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
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
container for various plainwaves related things
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
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
https://en.wikipedia.org/wiki/Finite_difference_coefficient
Definition qs_fxc.F:27
subroutine, public qs_fgxc_create(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, is_triplet, fxc_rho, fxc_tau, gxc_rho, gxc_tau)
...
Definition qs_fxc.F:629
subroutine, public qs_fgxc_gdiff(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, epsrho, is_triplet, weights, fxc_rho, fxc_tau, gxc_rho, gxc_tau, spinflip)
...
Definition qs_fxc.F:462
subroutine, public qs_fxc_analytic(rho0, rho1_r, tau1_r, xc_section, weights, auxbas_pw_pool, is_triplet, v_xc, v_xc_tau, spinflip)
...
Definition qs_fxc.F:96
subroutine, public qs_fgxc_release(ks_env, fxc_rho, fxc_tau, gxc_rho, gxc_tau)
...
Definition qs_fxc.F:758
subroutine, public qs_fxc_fdiff(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, is_triplet, fxc_rho, fxc_tau)
...
Definition qs_fxc.F:165
subroutine, public qs_fgxc_analytic(rho0_struct, rho1_struct, xc_section, weights, pw_pool, fxc_rho, fxc_tau, gxc_rho, gxc_tau, spinflip)
Calculates the values at the grid points in real space (r_i), of the second and third functional deri...
Definition qs_fxc.F:307
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, exc_accint, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, rho, rho_xc, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, task_list, task_list_soft, kpoints, do_kpoints, atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, particle_set, energy, force, local_particles, local_molecules, molecule_kind_set, molecule_set, subsys, cp_subsys, virial, results, atprop, nkind, natom, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env, nelectron_total, nelectron_spin)
...
methods of the rho structure (defined in qs_rho_types)
subroutine, public qs_rho_scale_and_add_b(rhoa, rhob, alpha, beta)
rhoa(2) = alpha*rhoa(2)+beta*rhob(1)
subroutine, public qs_rho_copy(rho_input, rho_output, auxbas_pw_pool, mspin)
Allocate a density structure and fill it with data from an input structure SIZE(rho_input) == mspin =...
subroutine, public qs_rho_scale_and_add(rhoa, rhob, alpha, beta)
rhoa = alpha*rhoa+beta*rhob
superstucture that hold various representations of the density and keeps track of which ones are vali...
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...
subroutine, public qs_rho_create(rho)
Allocates a new instance of rho.
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...
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)
calculates and allocates the xc potential, already reducing it to the dependence on rho and the one o...
Definition qs_vxc.F:98
represent a group ofunctional derivatives
subroutine, public xc_dset_release(derivative_set)
releases a derivative set
type(xc_rho_cflags_type) function, public xc_functionals_get_needs(functionals, lsd, calc_potential)
...
contains the structure
contains the structure
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
subroutine, public xc_rho_set_release(rho_set, pw_pool)
releases the given rho_set
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...
Exchange and Correlation functional calculations.
Definition xc.F:17
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:5533
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:2051
subroutine, public xc_calc_2nd_deriv(v_xc, v_xc_tau, deriv_set, rho_set, rho1_r, rho1_g, tau1_r, pw_pool, weights, xc_section, gapw, vxg, do_excitations, do_sf, do_triplet, compute_virial, virial_xc)
Caller routine to calculate the second order potential in the direction of rho1_r.
Definition xc.F:924
subroutine, public xc_calc_3rd_deriv_analytical(v_xc, v_xc_tau, deriv_set, rho_set, rho1_set, pw_pool, xc_section, spinflip)
Calculates the third functional derivative of the exchange-correlation functional,...
Definition xc.F:4729
subroutine, public xc_prep_3rd_deriv(deriv_set, rho_set, rho_r, pw_pool, weights, xc_section, tau_r, do_sf)
Prepare deriv_set for the calculation of the 3rd derivatives of the density functional....
Definition xc.F:5593
contained for different pw related things
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
keeps the density in various representations, keeping track of which ones are valid.
A derivative set contains the different derivatives of a xc-functional in form of a linked list.
contains a flag for each component of xc_rho_set, so that you can use it to tell which components you...
represent a density, with all the representation and data needed to perform a functional evaluation