Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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_get_ival,&
31 : section_get_rval,&
32 : section_vals_get_subs_vals,&
33 : section_vals_type
34 : USE kinds, ONLY: dp
35 : USE pw_env_types, ONLY: pw_env_get,&
36 : pw_env_type
37 : USE pw_methods, ONLY: pw_axpy,&
38 : pw_scale,&
39 : pw_zero
40 : USE pw_pool_types, ONLY: pw_pool_type
41 : USE pw_types, ONLY: pw_c1d_gs_type,&
42 : pw_r3d_rs_type
43 : USE qs_ks_types, ONLY: get_ks_env,&
44 : qs_ks_env_type
45 : USE qs_rho_methods, ONLY: qs_rho_copy,&
46 : qs_rho_scale_and_add,&
47 : qs_rho_scale_and_add_b
48 : USE qs_rho_types, ONLY: qs_rho_create,&
49 : qs_rho_get,&
50 : qs_rho_release,&
51 : qs_rho_type
52 : USE qs_vxc, ONLY: qs_vxc_create
53 : USE xc, ONLY: xc_calc_2nd_deriv,&
54 : xc_calc_2nd_deriv_analytical,&
55 : xc_calc_3rd_deriv_analytical,&
56 : xc_prep_2nd_deriv,&
57 : xc_prep_3rd_deriv
58 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
59 : xc_dset_release
60 : USE xc_derivatives, ONLY: xc_functionals_get_needs
61 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
62 : USE xc_rho_set_types, ONLY: xc_rho_set_create,&
63 : xc_rho_set_release,&
64 : xc_rho_set_type,&
65 : xc_rho_set_update
66 : #include "./base/base_uses.f90"
67 :
68 : IMPLICIT NONE
69 :
70 : PRIVATE
71 :
72 : ! *** Public subroutines ***
73 : PUBLIC :: qs_fxc_fdiff, qs_fxc_analytic, qs_fgxc_gdiff, qs_fgxc_analytic, qs_fgxc_create, qs_fgxc_release
74 :
75 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_fxc'
76 :
77 : ! **************************************************************************************************
78 :
79 : CONTAINS
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 19156 : 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 19156 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho0_g, rho1_g
113 19156 : 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 9578 : CALL timeset(routineN, handle)
120 :
121 9578 : CPASSERT(.NOT. ASSOCIATED(v_xc))
122 9578 : CPASSERT(.NOT. ASSOCIATED(v_xc_tau))
123 :
124 9578 : do_sf = .FALSE.
125 9578 : IF (PRESENT(spinflip)) do_sf = spinflip
126 :
127 9578 : CALL qs_rho_get(rho0, rho_r=rho0_r, rho_g=rho0_g, tau_r=tau0_r)
128 9578 : nspins = SIZE(rho0_r)
129 :
130 9578 : lsd = (nspins == 2)
131 : fac = 0._dp
132 9578 : IF (is_triplet .AND. nspins == 1) fac = -1.0_dp
133 :
134 9578 : NULLIFY (rho1_g)
135 95780 : bo = rho1_r(1)%pw_grid%bounds_local
136 9578 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
137 9578 : 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 9578 : 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 9578 : gapw=.FALSE., do_triplet=is_triplet, do_sf=do_sf)
145 9578 : CALL xc_dset_release(deriv_set)
146 9578 : CALL xc_rho_set_release(rho0_set)
147 :
148 9578 : CALL timestop(handle)
149 :
150 210716 : 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 2286 : 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 2286 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_tau_rspace, vxc00
183 : TYPE(qs_rho_type), POINTER :: rhoin
184 :
185 2286 : CALL timeset(routineN, handle)
186 :
187 2286 : CPASSERT(.NOT. ASSOCIATED(fxc_rho))
188 2286 : CPASSERT(.NOT. ASSOCIATED(fxc_tau))
189 2286 : CPASSERT(ASSOCIATED(rho0_struct))
190 2286 : CPASSERT(ASSOCIATED(rho1_struct))
191 :
192 2286 : ak = 0.0_dp
193 2286 : SELECT CASE (accuracy)
194 : CASE (:4)
195 0 : nstep = 2
196 0 : 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 18288 : nstep = 3
199 18288 : 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 0 : 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 2286 : 224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp]/280.0_dp
204 : END SELECT
205 :
206 2286 : CALL get_ks_env(ks_env, dft_control=dft_control, pw_env=pw_env)
207 2286 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
208 :
209 2286 : nspins = dft_control%nspins
210 2286 : exc = 0.0_dp
211 :
212 18288 : DO istep = -nstep, nstep
213 :
214 18288 : IF (ak(istep) /= 0.0_dp) THEN
215 13716 : alpha = 1.0_dp
216 13716 : beta = REAL(istep, KIND=dp)*epsrho
217 : NULLIFY (rhoin)
218 13716 : ALLOCATE (rhoin)
219 13716 : CALL qs_rho_create(rhoin)
220 13716 : NULLIFY (vxc00, v_tau_rspace)
221 13716 : IF (is_triplet) THEN
222 1176 : CPASSERT(nspins == 1)
223 : ! rhoin = (0.5 rho0, 0.5 rho0)
224 1176 : CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, 2)
225 : ! rhoin = (0.5 rho0 + 0.5 rho1, 0.5 rho0)
226 1176 : 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 1176 : vxc_rho=vxc00, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE.)
229 1176 : CALL pw_axpy(vxc00(2), vxc00(1), -1.0_dp)
230 1176 : IF (ASSOCIATED(v_tau_rspace)) CALL pw_axpy(v_tau_rspace(2), v_tau_rspace(1), -1.0_dp)
231 : ELSE
232 12540 : CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
233 12540 : 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 12540 : vxc_rho=vxc00, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE.)
236 : END IF
237 13716 : CALL qs_rho_release(rhoin)
238 13716 : DEALLOCATE (rhoin)
239 13716 : IF (.NOT. ASSOCIATED(fxc_rho)) THEN
240 9444 : ALLOCATE (fxc_rho(nspins))
241 4872 : DO ispin = 1, nspins
242 2586 : CALL auxbas_pw_pool%create_pw(fxc_rho(ispin))
243 4872 : CALL pw_zero(fxc_rho(ispin))
244 : END DO
245 : END IF
246 29232 : DO ispin = 1, nspins
247 29232 : CALL pw_axpy(vxc00(ispin), fxc_rho(ispin), ak(istep))
248 : END DO
249 30408 : DO ispin = 1, SIZE(vxc00)
250 30408 : CALL auxbas_pw_pool%give_back_pw(vxc00(ispin))
251 : END DO
252 13716 : DEALLOCATE (vxc00)
253 13716 : IF (ASSOCIATED(v_tau_rspace)) THEN
254 0 : IF (.NOT. ASSOCIATED(fxc_tau)) THEN
255 0 : ALLOCATE (fxc_tau(nspins))
256 0 : DO ispin = 1, nspins
257 0 : CALL auxbas_pw_pool%create_pw(fxc_tau(ispin))
258 0 : CALL pw_zero(fxc_tau(ispin))
259 : END DO
260 : END IF
261 0 : DO ispin = 1, nspins
262 0 : CALL pw_axpy(v_tau_rspace(ispin), fxc_tau(ispin), ak(istep))
263 : END DO
264 0 : DO ispin = 1, SIZE(v_tau_rspace)
265 0 : CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
266 : END DO
267 0 : DEALLOCATE (v_tau_rspace)
268 : END IF
269 : END IF
270 :
271 : END DO
272 :
273 2286 : oeps1 = 1.0_dp/epsrho
274 4872 : DO ispin = 1, nspins
275 4872 : CALL pw_scale(fxc_rho(ispin), oeps1)
276 : END DO
277 2286 : IF (ASSOCIATED(fxc_tau)) THEN
278 0 : DO ispin = 1, nspins
279 0 : CALL pw_scale(fxc_tau(ispin), oeps1)
280 : END DO
281 : END IF
282 :
283 2286 : CALL timestop(handle)
284 :
285 2286 : 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 0 : 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 0 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho0_g, rho1_g
322 0 : 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 0 : CALL timeset(routineN, handle)
329 :
330 : ! Only rho0 and rho1 should be associated
331 0 : CPASSERT(.NOT. ASSOCIATED(fxc_rho))
332 0 : CPASSERT(.NOT. ASSOCIATED(fxc_tau))
333 0 : CPASSERT(.NOT. ASSOCIATED(gxc_rho))
334 0 : CPASSERT(.NOT. ASSOCIATED(gxc_tau))
335 0 : CPASSERT(ASSOCIATED(rho0_struct))
336 0 : CPASSERT(ASSOCIATED(rho1_struct))
337 :
338 : ! Initialize parameters
339 0 : do_sf = .FALSE.
340 0 : IF (PRESENT(spinflip)) do_sf = spinflip
341 : !
342 : ! Get the values on the gridpoints of the rho0 density
343 0 : CALL qs_rho_get(rho0_struct, rho_r=rho0_r, rho_g=rho0_g, tau_r=tau0_r)
344 0 : nspins = SIZE(rho0_r)
345 0 : lsd = (nspins == 2)
346 : !
347 0 : IF (do_sf) THEN
348 : spindim = 1
349 : ELSE
350 0 : spindim = nspins
351 : END IF
352 : !
353 0 : fac = 0._dp
354 0 : IF (nspins == 1) THEN
355 0 : fac = 1.0_dp
356 : END IF
357 :
358 : ! Read xc functional section and find out what the functional actually needs
359 0 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
360 0 : needs = xc_functionals_get_needs(xc_fun_section, lsd, .TRUE.)
361 :
362 : ! Create fields for the kernel and kernel derivative
363 0 : ALLOCATE (fxc_rho(spindim), gxc_rho(nspins))
364 0 : DO ispin = 1, spindim
365 0 : CALL pw_pool%create_pw(fxc_rho(ispin))
366 0 : CALL pw_zero(fxc_rho(ispin))
367 : END DO
368 0 : DO ispin = 1, nspins
369 0 : CALL pw_pool%create_pw(gxc_rho(ispin))
370 0 : CALL pw_zero(gxc_rho(ispin))
371 : END DO
372 : ! Create fields for mGGA functionals. This implementation is not ready yet!
373 0 : IF (needs%tau .OR. needs%tau_spin) THEN
374 0 : IF (.NOT. ASSOCIATED(tau1_r)) &
375 0 : CPABORT("Tau-dependent functionals requires allocated kinetic energy density grid")
376 0 : ALLOCATE (fxc_tau(spindim), gxc_tau(nspins))
377 0 : DO ispin = 1, spindim
378 0 : CALL pw_pool%create_pw(fxc_tau(ispin))
379 0 : CALL pw_zero(fxc_tau(ispin))
380 : END DO
381 0 : DO ispin = 1, nspins
382 0 : CALL pw_pool%create_pw(gxc_tau(ispin))
383 0 : 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 0 : 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 0 : CALL qs_rho_get(rho1_struct, rho_r=rho1_r, rho_g=rho1_g, tau_r=tau1_r)
402 0 : 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 0 : 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 0 : 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 0 : 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 0 : xc_section, spinflip=do_sf)
435 :
436 0 : CALL xc_dset_release(deriv_set)
437 0 : CALL xc_rho_set_release(rho0_set)
438 0 : CALL xc_rho_set_release(rho1_set)
439 :
440 0 : CALL timestop(handle)
441 :
442 0 : 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 326 : 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 326 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r, tau1_r, v_tau_rspace, vxc00, &
483 326 : vxc00b
484 : TYPE(qs_rho_type), POINTER :: rhoin
485 :
486 326 : CALL timeset(routineN, handle)
487 :
488 326 : CPASSERT(.NOT. ASSOCIATED(fxc_rho))
489 326 : CPASSERT(.NOT. ASSOCIATED(fxc_tau))
490 326 : CPASSERT(.NOT. ASSOCIATED(gxc_rho))
491 326 : CPASSERT(.NOT. ASSOCIATED(gxc_tau))
492 326 : CPASSERT(ASSOCIATED(rho0_struct))
493 326 : CPASSERT(ASSOCIATED(rho1_struct))
494 :
495 326 : do_sf = .FALSE.
496 326 : IF (PRESENT(spinflip)) do_sf = spinflip
497 :
498 326 : ak = 0.0_dp
499 326 : SELECT CASE (accuracy)
500 : CASE (:4)
501 0 : nstep = 2
502 0 : 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 2608 : nstep = 3
505 2608 : 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 0 : 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 326 : 224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp]/280.0_dp
510 : END SELECT
511 :
512 326 : CALL get_ks_env(ks_env, dft_control=dft_control, pw_env=pw_env)
513 326 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
514 :
515 326 : nspins = dft_control%nspins
516 326 : exc = 0.0_dp
517 :
518 326 : IF (do_sf) THEN
519 8 : 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 8 : 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 318 : fxc_rho, fxc_tau)
526 : END IF
527 :
528 2608 : DO istep = -nstep, nstep
529 :
530 2608 : IF (ak(istep) /= 0.0_dp) THEN
531 1956 : alpha = 1.0_dp
532 1956 : beta = REAL(istep, KIND=dp)*epsrho
533 : NULLIFY (rhoin)
534 1956 : ALLOCATE (rhoin)
535 1956 : CALL qs_rho_create(rhoin)
536 1956 : NULLIFY (vxc00, vxc00b, v_tau_rspace)
537 1956 : CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
538 1956 : CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, beta)
539 1956 : 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 48 : vxc00, v_tau_rspace, spinflip=do_sf)
544 : ! variation in beta density
545 48 : CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
546 48 : 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 48 : 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 1908 : fxc_rho=vxc00, fxc_tau=v_tau_rspace)
554 : END IF
555 1956 : CALL qs_rho_release(rhoin)
556 1956 : DEALLOCATE (rhoin)
557 1956 : IF (.NOT. ASSOCIATED(gxc_rho)) THEN
558 1348 : ALLOCATE (gxc_rho(nspins))
559 696 : DO ispin = 1, nspins
560 370 : CALL auxbas_pw_pool%create_pw(gxc_rho(ispin))
561 696 : CALL pw_zero(gxc_rho(ispin))
562 : END DO
563 : END IF
564 1956 : IF (do_sf) THEN
565 48 : CALL pw_axpy(vxc00(1), gxc_rho(1), ak(istep))
566 48 : CALL pw_axpy(vxc00b(1), gxc_rho(2), ak(istep))
567 : ELSE
568 4032 : DO ispin = 1, nspins
569 4032 : CALL pw_axpy(vxc00(ispin), gxc_rho(ispin), ak(istep))
570 : END DO
571 : END IF
572 4128 : DO ispin = 1, SIZE(vxc00)
573 4128 : CALL auxbas_pw_pool%give_back_pw(vxc00(ispin))
574 : END DO
575 1956 : DEALLOCATE (vxc00)
576 1956 : IF (ASSOCIATED(vxc00b)) THEN
577 48 : CALL auxbas_pw_pool%give_back_pw(vxc00b(1))
578 48 : DEALLOCATE (vxc00b)
579 : END IF
580 1956 : IF (ASSOCIATED(v_tau_rspace)) THEN
581 0 : IF (.NOT. ASSOCIATED(gxc_tau)) THEN
582 0 : ALLOCATE (gxc_tau(nspins))
583 0 : DO ispin = 1, nspins
584 0 : CALL auxbas_pw_pool%create_pw(gxc_tau(ispin))
585 0 : CALL pw_zero(gxc_tau(ispin))
586 : END DO
587 : END IF
588 0 : DO ispin = 1, nspins
589 0 : CALL pw_axpy(v_tau_rspace(ispin), gxc_tau(ispin), ak(istep))
590 : END DO
591 0 : DO ispin = 1, SIZE(v_tau_rspace)
592 0 : CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
593 : END DO
594 0 : DEALLOCATE (v_tau_rspace)
595 : END IF
596 : END IF
597 :
598 : END DO
599 :
600 326 : oeps1 = 1.0_dp/epsrho
601 696 : DO ispin = 1, nspins
602 696 : CALL pw_scale(gxc_rho(ispin), oeps1)
603 : END DO
604 326 : IF (ASSOCIATED(gxc_tau)) THEN
605 0 : DO ispin = 1, nspins
606 0 : CALL pw_scale(gxc_tau(ispin), oeps1)
607 : END DO
608 : END IF
609 :
610 326 : CALL timestop(handle)
611 :
612 326 : 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 0 : 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 0 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_tau_rspace, vxc00
647 : TYPE(qs_rho_type), POINTER :: rhoin
648 :
649 0 : CALL timeset(routineN, handle)
650 :
651 0 : CPASSERT(.NOT. ASSOCIATED(fxc_rho))
652 0 : CPASSERT(.NOT. ASSOCIATED(fxc_tau))
653 0 : CPASSERT(.NOT. ASSOCIATED(gxc_rho))
654 0 : CPASSERT(.NOT. ASSOCIATED(gxc_tau))
655 0 : CPASSERT(ASSOCIATED(rho0_struct))
656 0 : CPASSERT(ASSOCIATED(rho1_struct))
657 :
658 0 : ak = 0.0_dp
659 0 : bl = 0.0_dp
660 0 : SELECT CASE (accuracy)
661 : CASE (:4)
662 0 : nstep = 2
663 0 : ak(-2:2) = [1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp]/12.0_dp
664 0 : 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 0 : nstep = 3
667 0 : ak(-3:3) = [-1.0_dp, 9.0_dp, -45.0_dp, 0.0_dp, 45.0_dp, -9.0_dp, 1.0_dp]/60.0_dp
668 0 : bl(-3:3) = [2.0_dp, -27.0_dp, 270.0_dp, -490.0_dp, 270.0_dp, -27.0_dp, 2.0_dp]/180.0_dp
669 : CASE (8:)
670 0 : 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 0 : 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 0 : 896.0_dp, -112.0_dp, 128.0_dp/9.0_dp, -1.0_dp]/560.0_dp
675 : END SELECT
676 :
677 0 : CALL get_ks_env(ks_env, dft_control=dft_control, pw_env=pw_env)
678 0 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
679 :
680 0 : nspins = dft_control%nspins
681 0 : exc = 0.0_dp
682 :
683 0 : DO istep = -nstep, nstep
684 :
685 0 : alpha = 1.0_dp
686 0 : beta = REAL(istep, KIND=dp)*epsrho
687 : NULLIFY (rhoin)
688 0 : ALLOCATE (rhoin)
689 0 : CALL qs_rho_create(rhoin)
690 0 : NULLIFY (vxc00, v_tau_rspace)
691 0 : IF (is_triplet) THEN
692 0 : CPASSERT(nspins == 1)
693 : ! rhoin = (0.5 rho0, 0.5 rho0)
694 0 : CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, 2)
695 : ! rhoin = (0.5 rho0 + 0.5 rho1, 0.5 rho0)
696 0 : 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 0 : vxc_rho=vxc00, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE.)
699 0 : CALL pw_axpy(vxc00(2), vxc00(1), -1.0_dp)
700 : ELSE
701 0 : CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
702 0 : 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 0 : vxc_rho=vxc00, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE.)
705 : END IF
706 0 : CALL qs_rho_release(rhoin)
707 0 : DEALLOCATE (rhoin)
708 0 : IF (.NOT. ASSOCIATED(fxc_rho)) THEN
709 0 : ALLOCATE (fxc_rho(nspins))
710 0 : DO ispin = 1, nspins
711 0 : CALL auxbas_pw_pool%create_pw(fxc_rho(ispin))
712 0 : CALL pw_zero(fxc_rho(ispin))
713 : END DO
714 : END IF
715 0 : IF (.NOT. ASSOCIATED(gxc_rho)) THEN
716 0 : ALLOCATE (gxc_rho(nspins))
717 0 : DO ispin = 1, nspins
718 0 : CALL auxbas_pw_pool%create_pw(gxc_rho(ispin))
719 0 : CALL pw_zero(gxc_rho(ispin))
720 : END DO
721 : END IF
722 0 : CPASSERT(.NOT. ASSOCIATED(v_tau_rspace))
723 0 : DO ispin = 1, nspins
724 0 : IF (ak(istep) /= 0.0_dp) THEN
725 0 : CALL pw_axpy(vxc00(ispin), fxc_rho(ispin), ak(istep))
726 : END IF
727 0 : IF (bl(istep) /= 0.0_dp) THEN
728 0 : CALL pw_axpy(vxc00(ispin), gxc_rho(ispin), bl(istep))
729 : END IF
730 : END DO
731 0 : DO ispin = 1, SIZE(vxc00)
732 0 : CALL auxbas_pw_pool%give_back_pw(vxc00(ispin))
733 : END DO
734 0 : DEALLOCATE (vxc00)
735 :
736 : END DO
737 :
738 0 : oeps1 = 1.0_dp/epsrho
739 0 : oeps2 = 1.0_dp/(epsrho**2)
740 0 : DO ispin = 1, nspins
741 0 : CALL pw_scale(fxc_rho(ispin), oeps1)
742 0 : CALL pw_scale(gxc_rho(ispin), oeps2)
743 : END DO
744 :
745 0 : CALL timestop(handle)
746 :
747 0 : 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 326 : 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 326 : CALL get_ks_env(ks_env, pw_env=pw_env)
767 326 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
768 :
769 326 : IF (ASSOCIATED(fxc_rho)) THEN
770 688 : DO ispin = 1, SIZE(fxc_rho)
771 688 : CALL auxbas_pw_pool%give_back_pw(fxc_rho(ispin))
772 : END DO
773 326 : DEALLOCATE (fxc_rho)
774 : END IF
775 326 : IF (ASSOCIATED(fxc_tau)) THEN
776 0 : DO ispin = 1, SIZE(fxc_tau)
777 0 : CALL auxbas_pw_pool%give_back_pw(fxc_tau(ispin))
778 : END DO
779 0 : DEALLOCATE (fxc_tau)
780 : END IF
781 326 : IF (ASSOCIATED(gxc_rho)) THEN
782 696 : DO ispin = 1, SIZE(gxc_rho)
783 696 : CALL auxbas_pw_pool%give_back_pw(gxc_rho(ispin))
784 : END DO
785 326 : DEALLOCATE (gxc_rho)
786 : END IF
787 326 : IF (ASSOCIATED(gxc_tau)) THEN
788 0 : DO ispin = 1, SIZE(gxc_tau)
789 0 : CALL auxbas_pw_pool%give_back_pw(gxc_tau(ispin))
790 : END DO
791 0 : DEALLOCATE (gxc_tau)
792 : END IF
793 :
794 326 : END SUBROUTINE qs_fgxc_release
795 :
796 : END MODULE qs_fxc
|
