(git:1a29073)
qs_fxc.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 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 ! **************************************************************************************************
27 MODULE qs_fxc
28 
29  USE cp_control_types, ONLY: dft_control_type
30  USE input_section_types, ONLY: section_vals_get_subs_vals,&
31  section_vals_type
32  USE kinds, ONLY: dp
33  USE pw_env_types, ONLY: pw_env_get,&
34  pw_env_type
35  USE pw_methods, ONLY: pw_axpy,&
36  pw_scale,&
37  pw_zero
38  USE pw_pool_types, ONLY: pw_pool_type
39  USE pw_types, ONLY: pw_c1d_gs_type,&
40  pw_r3d_rs_type
41  USE qs_ks_types, ONLY: get_ks_env,&
42  qs_ks_env_type
43  USE qs_rho_methods, ONLY: qs_rho_copy,&
44  qs_rho_scale_and_add
45  USE qs_rho_types, ONLY: qs_rho_create,&
46  qs_rho_get,&
47  qs_rho_release,&
48  qs_rho_type
49  USE qs_vxc, ONLY: qs_vxc_create
50  USE xc, ONLY: xc_calc_2nd_deriv,&
51  xc_prep_2nd_deriv
52  USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
53  xc_dset_release
54  USE xc_derivatives, ONLY: xc_functionals_get_needs
55  USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
56  USE xc_rho_set_types, ONLY: xc_rho_set_release,&
57  xc_rho_set_type
58 #include "./base/base_uses.f90"
59 
60  IMPLICIT NONE
61 
62  PRIVATE
63 
64  ! *** Public subroutines ***
65  PUBLIC :: qs_fxc_fdiff, qs_fxc_analytic, qs_fgxc_gdiff, qs_fgxc_create, qs_fgxc_release
66 
67  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_fxc'
68 
69 ! **************************************************************************************************
70 
71 CONTAINS
72 
73 ! **************************************************************************************************
74 !> \brief ...
75 !> \param rho0 ...
76 !> \param rho1_r ...
77 !> \param tau1_r ...
78 !> \param xc_section ...
79 !> \param auxbas_pw_pool ...
80 !> \param is_triplet ...
81 !> \param v_xc ...
82 !> \param v_xc_tau ...
83 ! **************************************************************************************************
84  SUBROUTINE qs_fxc_analytic(rho0, rho1_r, tau1_r, xc_section, auxbas_pw_pool, is_triplet, v_xc, v_xc_tau)
85 
86  TYPE(qs_rho_type), POINTER :: rho0
87  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r, tau1_r
88  TYPE(section_vals_type), POINTER :: xc_section
89  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
90  LOGICAL, INTENT(IN) :: is_triplet
91  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_xc, v_xc_tau
92 
93  CHARACTER(len=*), PARAMETER :: routinen = 'qs_fxc_analytic'
94 
95  INTEGER :: handle, nspins
96  INTEGER, DIMENSION(2, 3) :: bo
97  LOGICAL :: lsd
98  REAL(kind=dp) :: fac
99  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho0_g, rho1_g
100  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho0_r, tau0_r
101  TYPE(section_vals_type), POINTER :: xc_fun_section
102  TYPE(xc_derivative_set_type) :: deriv_set
103  TYPE(xc_rho_cflags_type) :: needs
104  TYPE(xc_rho_set_type) :: rho0_set
105 
106  CALL timeset(routinen, handle)
107 
108  cpassert(.NOT. ASSOCIATED(v_xc))
109  cpassert(.NOT. ASSOCIATED(v_xc_tau))
110 
111  CALL qs_rho_get(rho0, rho_r=rho0_r, rho_g=rho0_g, tau_r=tau0_r)
112  nspins = SIZE(rho0_r)
113 
114  lsd = (nspins == 2)
115  fac = 0._dp
116  IF (is_triplet .AND. nspins == 1) fac = -1.0_dp
117 
118  NULLIFY (rho1_g)
119  bo = rho1_r(1)%pw_grid%bounds_local
120  xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
121  needs = xc_functionals_get_needs(xc_fun_section, lsd, .true.)
122  ! calculate the arguments needed by the functionals
123  CALL xc_prep_2nd_deriv(deriv_set, rho0_set, rho0_r, auxbas_pw_pool, xc_section=xc_section, tau_r=tau0_r)
124  CALL xc_calc_2nd_deriv(v_xc, v_xc_tau, deriv_set, rho0_set, rho1_r, rho1_g, tau1_r, &
125  auxbas_pw_pool, xc_section=xc_section, gapw=.false., do_triplet=is_triplet)
126  CALL xc_dset_release(deriv_set)
127  CALL xc_rho_set_release(rho0_set)
128 
129  CALL timestop(handle)
130 
131  END SUBROUTINE qs_fxc_analytic
132 
133 ! **************************************************************************************************
134 !> \brief ...
135 !> \param ks_env ...
136 !> \param rho0_struct ...
137 !> \param rho1_struct ...
138 !> \param xc_section ...
139 !> \param accuracy ...
140 !> \param is_triplet ...
141 !> \param fxc_rho ...
142 !> \param fxc_tau ...
143 ! **************************************************************************************************
144  SUBROUTINE qs_fxc_fdiff(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, is_triplet, &
145  fxc_rho, fxc_tau)
146 
147  TYPE(qs_ks_env_type), POINTER :: ks_env
148  TYPE(qs_rho_type), POINTER :: rho0_struct, rho1_struct
149  TYPE(section_vals_type), POINTER :: xc_section
150  INTEGER, INTENT(IN) :: accuracy
151  LOGICAL, INTENT(IN) :: is_triplet
152  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau
153 
154  CHARACTER(len=*), PARAMETER :: routinen = 'qs_fxc_fdiff'
155  REAL(kind=dp), PARAMETER :: epsrho = 5.e-4_dp
156 
157  INTEGER :: handle, ispin, istep, nspins, nstep
158  REAL(kind=dp) :: alpha, beta, exc, oeps1
159  REAL(kind=dp), DIMENSION(-4:4) :: ak
160  TYPE(dft_control_type), POINTER :: dft_control
161  TYPE(pw_env_type), POINTER :: pw_env
162  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
163  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_tau_rspace, vxc00
164  TYPE(qs_rho_type), POINTER :: rhoin
165 
166  CALL timeset(routinen, handle)
167 
168  cpassert(.NOT. ASSOCIATED(fxc_rho))
169  cpassert(.NOT. ASSOCIATED(fxc_tau))
170  cpassert(ASSOCIATED(rho0_struct))
171  cpassert(ASSOCIATED(rho1_struct))
172 
173  ak = 0.0_dp
174  SELECT CASE (accuracy)
175  CASE (:4)
176  nstep = 2
177  ak(-2:2) = (/1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp/)/12.0_dp
178  CASE (5:7)
179  nstep = 3
180  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
181  CASE (8:)
182  nstep = 4
183  ak(-4:4) = (/1.0_dp, -32.0_dp/3.0_dp, 56.0_dp, -224.0_dp, 0.0_dp, &
184  224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp/)/280.0_dp
185  END SELECT
186 
187  CALL get_ks_env(ks_env, dft_control=dft_control, pw_env=pw_env)
188  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
189 
190  nspins = dft_control%nspins
191  exc = 0.0_dp
192 
193  DO istep = -nstep, nstep
194 
195  IF (ak(istep) /= 0.0_dp) THEN
196  alpha = 1.0_dp
197  beta = real(istep, kind=dp)*epsrho
198  NULLIFY (rhoin)
199  ALLOCATE (rhoin)
200  CALL qs_rho_create(rhoin)
201  NULLIFY (vxc00, v_tau_rspace)
202  IF (is_triplet) THEN
203  cpassert(nspins == 1)
204  ! rhoin = (0.5 rho0, 0.5 rho0)
205  CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, 2)
206  ! rhoin = (0.5 rho0 + 0.5 rho1, 0.5 rho0)
207  CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, 0.5_dp*beta)
208  CALL qs_vxc_create(ks_env=ks_env, rho_struct=rhoin, xc_section=xc_section, &
209  vxc_rho=vxc00, vxc_tau=v_tau_rspace, exc=exc, just_energy=.false.)
210  CALL pw_axpy(vxc00(2), vxc00(1), -1.0_dp)
211  IF (ASSOCIATED(v_tau_rspace)) CALL pw_axpy(v_tau_rspace(2), v_tau_rspace(1), -1.0_dp)
212  ELSE
213  CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
214  CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, beta)
215  CALL qs_vxc_create(ks_env=ks_env, rho_struct=rhoin, xc_section=xc_section, &
216  vxc_rho=vxc00, vxc_tau=v_tau_rspace, exc=exc, just_energy=.false.)
217  END IF
218  CALL qs_rho_release(rhoin)
219  DEALLOCATE (rhoin)
220  IF (.NOT. ASSOCIATED(fxc_rho)) THEN
221  ALLOCATE (fxc_rho(nspins))
222  DO ispin = 1, nspins
223  CALL auxbas_pw_pool%create_pw(fxc_rho(ispin))
224  CALL pw_zero(fxc_rho(ispin))
225  END DO
226  END IF
227  DO ispin = 1, nspins
228  CALL pw_axpy(vxc00(ispin), fxc_rho(ispin), ak(istep))
229  END DO
230  DO ispin = 1, SIZE(vxc00)
231  CALL auxbas_pw_pool%give_back_pw(vxc00(ispin))
232  END DO
233  DEALLOCATE (vxc00)
234  IF (ASSOCIATED(v_tau_rspace)) THEN
235  IF (.NOT. ASSOCIATED(fxc_tau)) THEN
236  ALLOCATE (fxc_tau(nspins))
237  DO ispin = 1, nspins
238  CALL auxbas_pw_pool%create_pw(fxc_tau(ispin))
239  CALL pw_zero(fxc_tau(ispin))
240  END DO
241  END IF
242  DO ispin = 1, nspins
243  CALL pw_axpy(v_tau_rspace(ispin), fxc_tau(ispin), ak(istep))
244  END DO
245  DO ispin = 1, SIZE(v_tau_rspace)
246  CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
247  END DO
248  DEALLOCATE (v_tau_rspace)
249  END IF
250  END IF
251 
252  END DO
253 
254  oeps1 = 1.0_dp/epsrho
255  DO ispin = 1, nspins
256  CALL pw_scale(fxc_rho(ispin), oeps1)
257  END DO
258  IF (ASSOCIATED(fxc_tau)) THEN
259  DO ispin = 1, nspins
260  CALL pw_scale(fxc_tau(ispin), oeps1)
261  END DO
262  END IF
263 
264  CALL timestop(handle)
265 
266  END SUBROUTINE qs_fxc_fdiff
267 
268 ! **************************************************************************************************
269 !> \brief ...
270 !> \param ks_env ...
271 !> \param rho0_struct ...
272 !> \param rho1_struct ...
273 !> \param xc_section ...
274 !> \param accuracy ...
275 !> \param epsrho ...
276 !> \param is_triplet ...
277 !> \param fxc_rho ...
278 !> \param fxc_tau ...
279 !> \param gxc_rho ...
280 !> \param gxc_tau ...
281 ! **************************************************************************************************
282  SUBROUTINE qs_fgxc_gdiff(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, epsrho, &
283  is_triplet, fxc_rho, fxc_tau, gxc_rho, gxc_tau)
284 
285  TYPE(qs_ks_env_type), POINTER :: ks_env
286  TYPE(qs_rho_type), POINTER :: rho0_struct, rho1_struct
287  TYPE(section_vals_type), POINTER :: xc_section
288  INTEGER, INTENT(IN) :: accuracy
289  REAL(kind=dp), INTENT(IN) :: epsrho
290  LOGICAL, INTENT(IN) :: is_triplet
291  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau, gxc_rho, gxc_tau
292 
293  CHARACTER(len=*), PARAMETER :: routinen = 'qs_fgxc_gdiff'
294 
295  INTEGER :: handle, ispin, istep, nspins, nstep
296  REAL(kind=dp) :: alpha, beta, exc, oeps1
297  REAL(kind=dp), DIMENSION(-4:4) :: ak
298  TYPE(dft_control_type), POINTER :: dft_control
299  TYPE(pw_env_type), POINTER :: pw_env
300  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
301  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_tau_rspace, vxc00
302  TYPE(qs_rho_type), POINTER :: rhoin
303 
304  CALL timeset(routinen, handle)
305 
306  cpassert(.NOT. ASSOCIATED(fxc_rho))
307  cpassert(.NOT. ASSOCIATED(fxc_tau))
308  cpassert(.NOT. ASSOCIATED(gxc_rho))
309  cpassert(.NOT. ASSOCIATED(gxc_tau))
310  cpassert(ASSOCIATED(rho0_struct))
311  cpassert(ASSOCIATED(rho1_struct))
312 
313  ak = 0.0_dp
314  SELECT CASE (accuracy)
315  CASE (:4)
316  nstep = 2
317  ak(-2:2) = (/1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp/)/12.0_dp
318  CASE (5:7)
319  nstep = 3
320  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
321  CASE (8:)
322  nstep = 4
323  ak(-4:4) = (/1.0_dp, -32.0_dp/3.0_dp, 56.0_dp, -224.0_dp, 0.0_dp, &
324  224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp/)/280.0_dp
325  END SELECT
326 
327  CALL get_ks_env(ks_env, dft_control=dft_control, pw_env=pw_env)
328  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
329 
330  nspins = dft_control%nspins
331  exc = 0.0_dp
332 
333  CALL qs_fxc_fdiff(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, is_triplet, &
334  fxc_rho, fxc_tau)
335 
336  DO istep = -nstep, nstep
337 
338  IF (ak(istep) /= 0.0_dp) THEN
339  alpha = 1.0_dp
340  beta = real(istep, kind=dp)*epsrho
341  NULLIFY (rhoin)
342  ALLOCATE (rhoin)
343  CALL qs_rho_create(rhoin)
344  NULLIFY (vxc00, v_tau_rspace)
345  CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
346  CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, beta)
347  CALL qs_fxc_fdiff(ks_env=ks_env, rho0_struct=rhoin, rho1_struct=rho1_struct, &
348  xc_section=xc_section, accuracy=accuracy, is_triplet=is_triplet, &
349  fxc_rho=vxc00, fxc_tau=v_tau_rspace)
350  CALL qs_rho_release(rhoin)
351  DEALLOCATE (rhoin)
352  IF (.NOT. ASSOCIATED(gxc_rho)) THEN
353  ALLOCATE (gxc_rho(nspins))
354  DO ispin = 1, nspins
355  CALL auxbas_pw_pool%create_pw(gxc_rho(ispin))
356  CALL pw_zero(gxc_rho(ispin))
357  END DO
358  END IF
359  DO ispin = 1, nspins
360  CALL pw_axpy(vxc00(ispin), gxc_rho(ispin), ak(istep))
361  END DO
362  DO ispin = 1, SIZE(vxc00)
363  CALL auxbas_pw_pool%give_back_pw(vxc00(ispin))
364  END DO
365  DEALLOCATE (vxc00)
366  IF (ASSOCIATED(v_tau_rspace)) THEN
367  IF (.NOT. ASSOCIATED(gxc_tau)) THEN
368  ALLOCATE (gxc_tau(nspins))
369  DO ispin = 1, nspins
370  CALL auxbas_pw_pool%create_pw(gxc_tau(ispin))
371  CALL pw_zero(gxc_tau(ispin))
372  END DO
373  END IF
374  DO ispin = 1, nspins
375  CALL pw_axpy(v_tau_rspace(ispin), gxc_tau(ispin), ak(istep))
376  END DO
377  DO ispin = 1, SIZE(v_tau_rspace)
378  CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
379  END DO
380  DEALLOCATE (v_tau_rspace)
381  END IF
382  END IF
383 
384  END DO
385 
386  oeps1 = 1.0_dp/epsrho
387  DO ispin = 1, nspins
388  CALL pw_scale(gxc_rho(ispin), oeps1)
389  END DO
390  IF (ASSOCIATED(gxc_tau)) THEN
391  DO ispin = 1, nspins
392  CALL pw_scale(gxc_tau(ispin), oeps1)
393  END DO
394  END IF
395 
396  CALL timestop(handle)
397 
398  END SUBROUTINE qs_fgxc_gdiff
399 
400 ! **************************************************************************************************
401 !> \brief ...
402 !> \param ks_env ...
403 !> \param rho0_struct ...
404 !> \param rho1_struct ...
405 !> \param xc_section ...
406 !> \param accuracy ...
407 !> \param is_triplet ...
408 !> \param fxc_rho ...
409 !> \param fxc_tau ...
410 !> \param gxc_rho ...
411 !> \param gxc_tau ...
412 ! **************************************************************************************************
413  SUBROUTINE qs_fgxc_create(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, is_triplet, &
414  fxc_rho, fxc_tau, gxc_rho, gxc_tau)
415 
416  TYPE(qs_ks_env_type), POINTER :: ks_env
417  TYPE(qs_rho_type), POINTER :: rho0_struct, rho1_struct
418  TYPE(section_vals_type), POINTER :: xc_section
419  INTEGER, INTENT(IN) :: accuracy
420  LOGICAL, INTENT(IN) :: is_triplet
421  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau, gxc_rho, gxc_tau
422 
423  CHARACTER(len=*), PARAMETER :: routinen = 'qs_fgxc_create'
424  REAL(kind=dp), PARAMETER :: epsrho = 5.e-4_dp
425 
426  INTEGER :: handle, ispin, istep, nspins, nstep
427  REAL(kind=dp) :: alpha, beta, exc, oeps1, oeps2
428  REAL(kind=dp), DIMENSION(-4:4) :: ak, bl
429  TYPE(dft_control_type), POINTER :: dft_control
430  TYPE(pw_env_type), POINTER :: pw_env
431  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
432  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_tau_rspace, vxc00
433  TYPE(qs_rho_type), POINTER :: rhoin
434 
435  CALL timeset(routinen, handle)
436 
437  cpassert(.NOT. ASSOCIATED(fxc_rho))
438  cpassert(.NOT. ASSOCIATED(fxc_tau))
439  cpassert(.NOT. ASSOCIATED(gxc_rho))
440  cpassert(.NOT. ASSOCIATED(gxc_tau))
441  cpassert(ASSOCIATED(rho0_struct))
442  cpassert(ASSOCIATED(rho1_struct))
443 
444  ak = 0.0_dp
445  bl = 0.0_dp
446  SELECT CASE (accuracy)
447  CASE (:4)
448  nstep = 2
449  ak(-2:2) = (/1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp/)/12.0_dp
450  bl(-2:2) = (/-1.0_dp, 16.0_dp, -30.0_dp, 16.0_dp, -1.0_dp/)/12.0_dp
451  CASE (5:7)
452  nstep = 3
453  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
454  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
455  CASE (8:)
456  nstep = 4
457  ak(-4:4) = (/1.0_dp, -32.0_dp/3.0_dp, 56.0_dp, -224.0_dp, 0.0_dp, &
458  224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp/)/280.0_dp
459  bl(-4:4) = (/-1.0_dp, 128.0_dp/9.0_dp, -112.0_dp, 896.0_dp, -14350.0_dp/9.0_dp, &
460  896.0_dp, -112.0_dp, 128.0_dp/9.0_dp, -1.0_dp/)/560.0_dp
461  END SELECT
462 
463  CALL get_ks_env(ks_env, dft_control=dft_control, pw_env=pw_env)
464  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
465 
466  nspins = dft_control%nspins
467  exc = 0.0_dp
468 
469  DO istep = -nstep, nstep
470 
471  alpha = 1.0_dp
472  beta = real(istep, kind=dp)*epsrho
473  NULLIFY (rhoin)
474  ALLOCATE (rhoin)
475  CALL qs_rho_create(rhoin)
476  NULLIFY (vxc00, v_tau_rspace)
477  IF (is_triplet) THEN
478  cpassert(nspins == 1)
479  ! rhoin = (0.5 rho0, 0.5 rho0)
480  CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, 2)
481  ! rhoin = (0.5 rho0 + 0.5 rho1, 0.5 rho0)
482  CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, 0.5_dp*beta)
483  CALL qs_vxc_create(ks_env=ks_env, rho_struct=rhoin, xc_section=xc_section, &
484  vxc_rho=vxc00, vxc_tau=v_tau_rspace, exc=exc, just_energy=.false.)
485  CALL pw_axpy(vxc00(2), vxc00(1), -1.0_dp)
486  ELSE
487  CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
488  CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, beta)
489  CALL qs_vxc_create(ks_env=ks_env, rho_struct=rhoin, xc_section=xc_section, &
490  vxc_rho=vxc00, vxc_tau=v_tau_rspace, exc=exc, just_energy=.false.)
491  END IF
492  CALL qs_rho_release(rhoin)
493  DEALLOCATE (rhoin)
494  IF (.NOT. ASSOCIATED(fxc_rho)) THEN
495  ALLOCATE (fxc_rho(nspins))
496  DO ispin = 1, nspins
497  CALL auxbas_pw_pool%create_pw(fxc_rho(ispin))
498  CALL pw_zero(fxc_rho(ispin))
499  END DO
500  END IF
501  IF (.NOT. ASSOCIATED(gxc_rho)) THEN
502  ALLOCATE (gxc_rho(nspins))
503  DO ispin = 1, nspins
504  CALL auxbas_pw_pool%create_pw(gxc_rho(ispin))
505  CALL pw_zero(gxc_rho(ispin))
506  END DO
507  END IF
508  cpassert(.NOT. ASSOCIATED(v_tau_rspace))
509  DO ispin = 1, nspins
510  IF (ak(istep) /= 0.0_dp) THEN
511  CALL pw_axpy(vxc00(ispin), fxc_rho(ispin), ak(istep))
512  END IF
513  IF (bl(istep) /= 0.0_dp) THEN
514  CALL pw_axpy(vxc00(ispin), gxc_rho(ispin), bl(istep))
515  END IF
516  END DO
517  DO ispin = 1, SIZE(vxc00)
518  CALL auxbas_pw_pool%give_back_pw(vxc00(ispin))
519  END DO
520  DEALLOCATE (vxc00)
521 
522  END DO
523 
524  oeps1 = 1.0_dp/epsrho
525  oeps2 = 1.0_dp/(epsrho**2)
526  DO ispin = 1, nspins
527  CALL pw_scale(fxc_rho(ispin), oeps1)
528  CALL pw_scale(gxc_rho(ispin), oeps2)
529  END DO
530 
531  CALL timestop(handle)
532 
533  END SUBROUTINE qs_fgxc_create
534 
535 ! **************************************************************************************************
536 !> \brief ...
537 !> \param ks_env ...
538 !> \param fxc_rho ...
539 !> \param fxc_tau ...
540 !> \param gxc_rho ...
541 !> \param gxc_tau ...
542 ! **************************************************************************************************
543  SUBROUTINE qs_fgxc_release(ks_env, fxc_rho, fxc_tau, gxc_rho, gxc_tau)
544 
545  TYPE(qs_ks_env_type), POINTER :: ks_env
546  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau, gxc_rho, gxc_tau
547 
548  INTEGER :: ispin
549  TYPE(pw_env_type), POINTER :: pw_env
550  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
551 
552  CALL get_ks_env(ks_env, pw_env=pw_env)
553  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
554 
555  IF (ASSOCIATED(fxc_rho)) THEN
556  DO ispin = 1, SIZE(fxc_rho)
557  CALL auxbas_pw_pool%give_back_pw(fxc_rho(ispin))
558  END DO
559  DEALLOCATE (fxc_rho)
560  END IF
561  IF (ASSOCIATED(fxc_tau)) THEN
562  DO ispin = 1, SIZE(fxc_tau)
563  CALL auxbas_pw_pool%give_back_pw(fxc_tau(ispin))
564  END DO
565  DEALLOCATE (fxc_tau)
566  END IF
567  IF (ASSOCIATED(gxc_rho)) THEN
568  DO ispin = 1, SIZE(gxc_rho)
569  CALL auxbas_pw_pool%give_back_pw(gxc_rho(ispin))
570  END DO
571  DEALLOCATE (gxc_rho)
572  END IF
573  IF (ASSOCIATED(gxc_tau)) THEN
574  DO ispin = 1, SIZE(gxc_tau)
575  CALL auxbas_pw_pool%give_back_pw(gxc_tau(ispin))
576  END DO
577  DEALLOCATE (gxc_tau)
578  END IF
579 
580  END SUBROUTINE qs_fgxc_release
581 
582 END MODULE qs_fxc
fac
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition: grid_common.h:48
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_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:724
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_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_fxc
https://en.wikipedia.org/wiki/Finite_difference_coefficient
Definition: qs_fxc.F:27
qs_fxc::qs_fgxc_gdiff
subroutine, public qs_fgxc_gdiff(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, epsrho, is_triplet, fxc_rho, fxc_tau, gxc_rho, gxc_tau)
...
Definition: qs_fxc.F:284
qs_fxc::qs_fgxc_create
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:415
qs_fxc::qs_fxc_analytic
subroutine, public qs_fxc_analytic(rho0, rho1_r, tau1_r, xc_section, auxbas_pw_pool, is_triplet, v_xc, v_xc_tau)
...
Definition: qs_fxc.F:85
qs_fxc::qs_fgxc_release
subroutine, public qs_fgxc_release(ks_env, fxc_rho, fxc_tau, gxc_rho, gxc_tau)
...
Definition: qs_fxc.F:544
qs_fxc::qs_fxc_fdiff
subroutine, public qs_fxc_fdiff(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, is_triplet, fxc_rho, fxc_tau)
...
Definition: qs_fxc.F:146
qs_ks_types
Definition: qs_ks_types.F:15
qs_ks_types::get_ks_env
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, 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, 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_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, 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)
...
Definition: qs_ks_types.F:330
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)
Allocate a density structure and fill it with data from an input structure SIZE(rho_input) == mspin =...
Definition: qs_rho_methods.F:766
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:1060
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:99
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:113
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)
calculates and allocates the xc potential, already reducing it to the dependence on rho and the one o...
Definition: qs_vxc.F:98
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:144
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:600
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_release
subroutine, public xc_rho_set_release(rho_set, pw_pool)
releases the given rho_set
Definition: xc_rho_set_types.F:129
xc
Exchange and Correlation functional calculations.
Definition: xc.F:17
xc::xc_calc_2nd_deriv
subroutine, public xc_calc_2nd_deriv(v_xc, v_xc_tau, deriv_set, rho_set, rho1_r, rho1_g, tau1_r, pw_pool, xc_section, gapw, vxg, lsd_singlets, do_excitations, do_triplet, do_tddft, compute_virial, virial_xc)
Caller routine to calculate the second order potential in the direction of rho1_r.
Definition: xc.F:1523
xc::xc_prep_2nd_deriv
subroutine, public xc_prep_2nd_deriv(deriv_set, rho_set, rho_r, pw_pool, xc_section, tau_r)
Prepare objects for the calculation of the 2nd derivatives of the density functional....
Definition: xc.F:5364