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 19772 : 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 19772 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho0_g, rho1_g
113 19772 : 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 9886 : CALL timeset(routineN, handle)
120 :
121 9886 : CPASSERT(.NOT. ASSOCIATED(v_xc))
122 9886 : CPASSERT(.NOT. ASSOCIATED(v_xc_tau))
123 :
124 9886 : do_sf = .FALSE.
125 9886 : IF (PRESENT(spinflip)) do_sf = spinflip
126 :
127 9886 : CALL qs_rho_get(rho0, rho_r=rho0_r, rho_g=rho0_g, tau_r=tau0_r)
128 9886 : nspins = SIZE(rho0_r)
129 :
130 9886 : lsd = (nspins == 2)
131 : fac = 0._dp
132 9886 : IF (is_triplet .AND. nspins == 1) fac = -1.0_dp
133 :
134 9886 : NULLIFY (rho1_g)
135 98860 : bo = rho1_r(1)%pw_grid%bounds_local
136 9886 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
137 9886 : 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 9886 : 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 9886 : gapw=.FALSE., do_triplet=is_triplet, do_sf=do_sf)
145 9886 : CALL xc_dset_release(deriv_set)
146 9886 : CALL xc_rho_set_release(rho0_set)
147 :
148 9886 : CALL timestop(handle)
149 :
150 217492 : 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 2378 : 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 2378 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_tau_rspace, vxc00
183 : TYPE(qs_rho_type), POINTER :: rhoin
184 :
185 2378 : CALL timeset(routineN, handle)
186 :
187 2378 : CPASSERT(.NOT. ASSOCIATED(fxc_rho))
188 2378 : CPASSERT(.NOT. ASSOCIATED(fxc_tau))
189 2378 : CPASSERT(ASSOCIATED(rho0_struct))
190 2378 : CPASSERT(ASSOCIATED(rho1_struct))
191 :
192 2378 : ak = 0.0_dp
193 2378 : 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 19024 : nstep = 3
199 19024 : 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 2378 : 224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp]/280.0_dp
204 : END SELECT
205 :
206 2378 : CALL get_ks_env(ks_env, dft_control=dft_control, pw_env=pw_env)
207 2378 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
208 :
209 2378 : nspins = dft_control%nspins
210 : exc = 0.0_dp
211 :
212 19024 : DO istep = -nstep, nstep
213 :
214 19024 : IF (ak(istep) /= 0.0_dp) THEN
215 14268 : alpha = 1.0_dp
216 14268 : beta = REAL(istep, KIND=dp)*epsrho
217 : NULLIFY (rhoin)
218 14268 : ALLOCATE (rhoin)
219 14268 : CALL qs_rho_create(rhoin)
220 14268 : NULLIFY (vxc00, v_tau_rspace)
221 14268 : 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 13092 : CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
233 13092 : 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 13092 : vxc_rho=vxc00, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE.)
236 : END IF
237 14268 : CALL qs_rho_release(rhoin)
238 14268 : DEALLOCATE (rhoin)
239 14268 : IF (.NOT. ASSOCIATED(fxc_rho)) THEN
240 9818 : ALLOCATE (fxc_rho(nspins))
241 5062 : DO ispin = 1, nspins
242 2684 : CALL auxbas_pw_pool%create_pw(fxc_rho(ispin))
243 5062 : CALL pw_zero(fxc_rho(ispin))
244 : END DO
245 : END IF
246 30372 : DO ispin = 1, nspins
247 30372 : CALL pw_axpy(vxc00(ispin), fxc_rho(ispin), ak(istep))
248 : END DO
249 31548 : DO ispin = 1, SIZE(vxc00)
250 31548 : CALL auxbas_pw_pool%give_back_pw(vxc00(ispin))
251 : END DO
252 14268 : DEALLOCATE (vxc00)
253 14268 : 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 2378 : oeps1 = 1.0_dp/epsrho
274 5062 : DO ispin = 1, nspins
275 5062 : CALL pw_scale(fxc_rho(ispin), oeps1)
276 : END DO
277 2378 : 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 2378 : CALL timestop(handle)
284 :
285 2378 : 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)) THEN
375 0 : CPABORT("Tau-dependent functionals requires allocated kinetic energy density grid")
376 : END IF
377 0 : ALLOCATE (fxc_tau(spindim), gxc_tau(nspins))
378 0 : DO ispin = 1, spindim
379 0 : CALL pw_pool%create_pw(fxc_tau(ispin))
380 0 : CALL pw_zero(fxc_tau(ispin))
381 : END DO
382 0 : DO ispin = 1, nspins
383 0 : CALL pw_pool%create_pw(gxc_tau(ispin))
384 0 : CALL pw_zero(gxc_tau(ispin))
385 : END DO
386 : END IF
387 :
388 : ! Build rho0_set
389 : ! calculate the arguments needed by the functionals
390 : ! Needs
391 : ! deriv_set xc_derivative_set_type just declared
392 : ! rho0_set xc_rho_set_type just declared
393 : ! rho0_r pw_type calculated by qs_rho_get
394 : ! pw_pool given by the calling subroutine
395 : ! xc_section given by the calling subroutine
396 : ! tau0_r pw_type calculated by qs_rho_get
397 : CALL xc_prep_3rd_deriv(deriv_set, rho0_set, rho0_r, pw_pool, weights, &
398 0 : xc_section, tau_r=tau0_r, do_sf=do_sf)
399 :
400 : ! Build rho1_set
401 : ! Get the values on the gridpoints of the rho1 density
402 0 : CALL qs_rho_get(rho1_struct, rho_r=rho1_r, rho_g=rho1_g, tau_r=tau1_r)
403 0 : bo = rho1_r(1)%pw_grid%bounds_local
404 : ! create the place where to store the argument for the functionals
405 : ! Needs
406 : ! rho1_set xc_rho_set_type just declared
407 : ! bo 2x3 integer matrix should have bounds_local or rho1_r
408 : CALL xc_rho_set_create(rho1_set, bo, &
409 : rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
410 : drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
411 0 : tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
412 : ! calculate the arguments needed by the functionals
413 : ! Needs
414 : ! rho1_set object created by xc_rho_set_create
415 : ! rho1_r,rho1_g,tau1_r pw_type values of rho1 in real space grid
416 : ! needs xc_rho_cflags_type defined through xc_functionals_get_needs
417 : ! pw_pool Given by the calling subroutine
418 : CALL xc_rho_set_update(rho1_set, rho1_r, rho1_g, tau1_r, needs, &
419 : section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
420 : section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
421 0 : pw_pool, spinflip=do_sf)
422 :
423 : ! Calculate exchange correlation kernel
424 : ! Needs
425 : ! fxc_rho, fxc_tau pw_type not associated
426 : ! deriv_set created and defined by xc_prep_3rd_deriv
427 : ! rho0_set xc_rho_set_type build by xc_prep_3rd_deriv
428 : ! rho1_set xc_rho_set_type build by xc_rho_set_create/update
429 : ! pw_pool needs to be given by the calling subroutine
430 : ! xc_section needs to be given by the calling subroutine
431 : CALL xc_calc_2nd_deriv_analytical(fxc_rho, fxc_tau, deriv_set, rho0_set, rho1_set, pw_pool, &
432 0 : xc_section, .FALSE., spinflip=do_sf, tddfpt_fac=fac)
433 : ! Calculate exchange correlation kernel derivative
434 : CALL xc_calc_3rd_deriv_analytical(gxc_rho, gxc_tau, deriv_set, rho0_set, rho1_set, pw_pool, &
435 0 : xc_section, spinflip=do_sf)
436 :
437 0 : CALL xc_dset_release(deriv_set)
438 0 : CALL xc_rho_set_release(rho0_set)
439 0 : CALL xc_rho_set_release(rho1_set)
440 :
441 0 : CALL timestop(handle)
442 :
443 0 : END SUBROUTINE qs_fgxc_analytic
444 :
445 : ! **************************************************************************************************
446 : !> \brief ...
447 : !> \param ks_env ...
448 : !> \param rho0_struct ...
449 : !> \param rho1_struct ...
450 : !> \param xc_section ...
451 : !> \param accuracy ...
452 : !> \param epsrho ...
453 : !> \param is_triplet ...
454 : !> \param weights ...
455 : !> \param fxc_rho ...
456 : !> \param fxc_tau ...
457 : !> \param gxc_rho ...
458 : !> \param gxc_tau ...
459 : !> \param spinflip ...
460 : ! **************************************************************************************************
461 334 : SUBROUTINE qs_fgxc_gdiff(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, epsrho, &
462 : is_triplet, weights, fxc_rho, fxc_tau, gxc_rho, gxc_tau, spinflip)
463 :
464 : TYPE(qs_ks_env_type), POINTER :: ks_env
465 : TYPE(qs_rho_type), POINTER :: rho0_struct, rho1_struct
466 : TYPE(section_vals_type), POINTER :: xc_section
467 : INTEGER, INTENT(IN) :: accuracy
468 : REAL(KIND=dp), INTENT(IN) :: epsrho
469 : LOGICAL, INTENT(IN) :: is_triplet
470 : TYPE(pw_r3d_rs_type), POINTER :: weights
471 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau, gxc_rho, gxc_tau
472 : LOGICAL, OPTIONAL :: spinflip
473 :
474 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_fgxc_gdiff'
475 :
476 : INTEGER :: handle, ispin, istep, nspins, nstep
477 : LOGICAL :: do_sf
478 : REAL(KIND=dp) :: alpha, beta, exc, oeps1
479 : REAL(KIND=dp), DIMENSION(-4:4) :: ak
480 : TYPE(dft_control_type), POINTER :: dft_control
481 : TYPE(pw_env_type), POINTER :: pw_env
482 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
483 334 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r, tau1_r, v_tau_rspace, vxc00, &
484 334 : vxc00b
485 : TYPE(qs_rho_type), POINTER :: rhoin
486 :
487 334 : CALL timeset(routineN, handle)
488 :
489 334 : CPASSERT(.NOT. ASSOCIATED(fxc_rho))
490 334 : CPASSERT(.NOT. ASSOCIATED(fxc_tau))
491 334 : CPASSERT(.NOT. ASSOCIATED(gxc_rho))
492 334 : CPASSERT(.NOT. ASSOCIATED(gxc_tau))
493 334 : CPASSERT(ASSOCIATED(rho0_struct))
494 334 : CPASSERT(ASSOCIATED(rho1_struct))
495 :
496 334 : do_sf = .FALSE.
497 334 : IF (PRESENT(spinflip)) do_sf = spinflip
498 :
499 334 : ak = 0.0_dp
500 334 : SELECT CASE (accuracy)
501 : CASE (:4)
502 0 : nstep = 2
503 0 : ak(-2:2) = [1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp]/12.0_dp
504 : CASE (5:7)
505 2672 : nstep = 3
506 2672 : 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
507 : CASE (8:)
508 0 : nstep = 4
509 : ak(-4:4) = [1.0_dp, -32.0_dp/3.0_dp, 56.0_dp, -224.0_dp, 0.0_dp, &
510 334 : 224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp]/280.0_dp
511 : END SELECT
512 :
513 334 : CALL get_ks_env(ks_env, dft_control=dft_control, pw_env=pw_env)
514 334 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
515 :
516 334 : nspins = dft_control%nspins
517 334 : exc = 0.0_dp
518 :
519 334 : IF (do_sf) THEN
520 8 : CALL qs_rho_get(rho1_struct, rho_r=rho1_r, tau_r=tau1_r)
521 : CALL qs_fxc_analytic(rho0_struct, rho1_r, tau1_r, xc_section, &
522 : weights, auxbas_pw_pool, is_triplet, &
523 8 : fxc_rho, fxc_tau, spinflip=do_sf)
524 : ELSE
525 : CALL qs_fxc_fdiff(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, is_triplet, &
526 326 : fxc_rho, fxc_tau)
527 : END IF
528 :
529 2672 : DO istep = -nstep, nstep
530 :
531 2672 : IF (ak(istep) /= 0.0_dp) THEN
532 2004 : alpha = 1.0_dp
533 2004 : beta = REAL(istep, KIND=dp)*epsrho
534 : NULLIFY (rhoin)
535 2004 : ALLOCATE (rhoin)
536 2004 : CALL qs_rho_create(rhoin)
537 2004 : NULLIFY (vxc00, vxc00b, v_tau_rspace)
538 2004 : CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
539 2004 : CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, beta)
540 2004 : IF (do_sf) THEN
541 : ! variation in alpha density
542 : CALL qs_fxc_analytic(rhoin, rho1_r, tau1_r, &
543 : xc_section, weights, auxbas_pw_pool, is_triplet, &
544 48 : vxc00, v_tau_rspace, spinflip=do_sf)
545 : ! variation in beta density
546 48 : CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
547 48 : CALL qs_rho_scale_and_add_b(rhoin, rho1_struct, alpha, beta)
548 : CALL qs_fxc_analytic(rhoin, rho1_r, tau1_r, &
549 : xc_section, weights, auxbas_pw_pool, is_triplet, &
550 48 : vxc00b, v_tau_rspace, spinflip=do_sf)
551 : ELSE
552 : CALL qs_fxc_fdiff(ks_env=ks_env, rho0_struct=rhoin, rho1_struct=rho1_struct, &
553 : xc_section=xc_section, accuracy=accuracy, is_triplet=is_triplet, &
554 1956 : fxc_rho=vxc00, fxc_tau=v_tau_rspace)
555 : END IF
556 2004 : CALL qs_rho_release(rhoin)
557 2004 : DEALLOCATE (rhoin)
558 2004 : IF (.NOT. ASSOCIATED(gxc_rho)) THEN
559 1380 : ALLOCATE (gxc_rho(nspins))
560 712 : DO ispin = 1, nspins
561 378 : CALL auxbas_pw_pool%create_pw(gxc_rho(ispin))
562 712 : CALL pw_zero(gxc_rho(ispin))
563 : END DO
564 : END IF
565 2004 : IF (do_sf) THEN
566 48 : CALL pw_axpy(vxc00(1), gxc_rho(1), ak(istep))
567 48 : CALL pw_axpy(vxc00b(1), gxc_rho(2), ak(istep))
568 : ELSE
569 4128 : DO ispin = 1, nspins
570 4128 : CALL pw_axpy(vxc00(ispin), gxc_rho(ispin), ak(istep))
571 : END DO
572 : END IF
573 4224 : DO ispin = 1, SIZE(vxc00)
574 4224 : CALL auxbas_pw_pool%give_back_pw(vxc00(ispin))
575 : END DO
576 2004 : DEALLOCATE (vxc00)
577 2004 : IF (ASSOCIATED(vxc00b)) THEN
578 48 : CALL auxbas_pw_pool%give_back_pw(vxc00b(1))
579 48 : DEALLOCATE (vxc00b)
580 : END IF
581 2004 : IF (ASSOCIATED(v_tau_rspace)) THEN
582 0 : IF (.NOT. ASSOCIATED(gxc_tau)) THEN
583 0 : ALLOCATE (gxc_tau(nspins))
584 0 : DO ispin = 1, nspins
585 0 : CALL auxbas_pw_pool%create_pw(gxc_tau(ispin))
586 0 : CALL pw_zero(gxc_tau(ispin))
587 : END DO
588 : END IF
589 0 : DO ispin = 1, nspins
590 0 : CALL pw_axpy(v_tau_rspace(ispin), gxc_tau(ispin), ak(istep))
591 : END DO
592 0 : DO ispin = 1, SIZE(v_tau_rspace)
593 0 : CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
594 : END DO
595 0 : DEALLOCATE (v_tau_rspace)
596 : END IF
597 : END IF
598 :
599 : END DO
600 :
601 334 : oeps1 = 1.0_dp/epsrho
602 712 : DO ispin = 1, nspins
603 712 : CALL pw_scale(gxc_rho(ispin), oeps1)
604 : END DO
605 334 : IF (ASSOCIATED(gxc_tau)) THEN
606 0 : DO ispin = 1, nspins
607 0 : CALL pw_scale(gxc_tau(ispin), oeps1)
608 : END DO
609 : END IF
610 :
611 334 : CALL timestop(handle)
612 :
613 334 : END SUBROUTINE qs_fgxc_gdiff
614 :
615 : ! **************************************************************************************************
616 : !> \brief ...
617 : !> \param ks_env ...
618 : !> \param rho0_struct ...
619 : !> \param rho1_struct ...
620 : !> \param xc_section ...
621 : !> \param accuracy ...
622 : !> \param is_triplet ...
623 : !> \param fxc_rho ...
624 : !> \param fxc_tau ...
625 : !> \param gxc_rho ...
626 : !> \param gxc_tau ...
627 : ! **************************************************************************************************
628 0 : SUBROUTINE qs_fgxc_create(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, is_triplet, &
629 : fxc_rho, fxc_tau, gxc_rho, gxc_tau)
630 :
631 : TYPE(qs_ks_env_type), POINTER :: ks_env
632 : TYPE(qs_rho_type), POINTER :: rho0_struct, rho1_struct
633 : TYPE(section_vals_type), POINTER :: xc_section
634 : INTEGER, INTENT(IN) :: accuracy
635 : LOGICAL, INTENT(IN) :: is_triplet
636 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau, gxc_rho, gxc_tau
637 :
638 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_fgxc_create'
639 : REAL(KIND=dp), PARAMETER :: epsrho = 5.e-4_dp
640 :
641 : INTEGER :: handle, ispin, istep, nspins, nstep
642 : REAL(KIND=dp) :: alpha, beta, exc, oeps1, oeps2
643 : REAL(KIND=dp), DIMENSION(-4:4) :: ak, bl
644 : TYPE(dft_control_type), POINTER :: dft_control
645 : TYPE(pw_env_type), POINTER :: pw_env
646 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
647 0 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_tau_rspace, vxc00
648 : TYPE(qs_rho_type), POINTER :: rhoin
649 :
650 0 : CALL timeset(routineN, handle)
651 :
652 0 : CPASSERT(.NOT. ASSOCIATED(fxc_rho))
653 0 : CPASSERT(.NOT. ASSOCIATED(fxc_tau))
654 0 : CPASSERT(.NOT. ASSOCIATED(gxc_rho))
655 0 : CPASSERT(.NOT. ASSOCIATED(gxc_tau))
656 0 : CPASSERT(ASSOCIATED(rho0_struct))
657 0 : CPASSERT(ASSOCIATED(rho1_struct))
658 :
659 0 : ak = 0.0_dp
660 0 : bl = 0.0_dp
661 0 : SELECT CASE (accuracy)
662 : CASE (:4)
663 0 : nstep = 2
664 0 : ak(-2:2) = [1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp]/12.0_dp
665 0 : bl(-2:2) = [-1.0_dp, 16.0_dp, -30.0_dp, 16.0_dp, -1.0_dp]/12.0_dp
666 : CASE (5:7)
667 0 : nstep = 3
668 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
669 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
670 : CASE (8:)
671 0 : nstep = 4
672 : ak(-4:4) = [1.0_dp, -32.0_dp/3.0_dp, 56.0_dp, -224.0_dp, 0.0_dp, &
673 0 : 224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp]/280.0_dp
674 : bl(-4:4) = [-1.0_dp, 128.0_dp/9.0_dp, -112.0_dp, 896.0_dp, -14350.0_dp/9.0_dp, &
675 0 : 896.0_dp, -112.0_dp, 128.0_dp/9.0_dp, -1.0_dp]/560.0_dp
676 : END SELECT
677 :
678 0 : CALL get_ks_env(ks_env, dft_control=dft_control, pw_env=pw_env)
679 0 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
680 :
681 0 : nspins = dft_control%nspins
682 : exc = 0.0_dp
683 :
684 0 : DO istep = -nstep, nstep
685 :
686 0 : alpha = 1.0_dp
687 0 : beta = REAL(istep, KIND=dp)*epsrho
688 : NULLIFY (rhoin)
689 0 : ALLOCATE (rhoin)
690 0 : CALL qs_rho_create(rhoin)
691 0 : NULLIFY (vxc00, v_tau_rspace)
692 0 : IF (is_triplet) THEN
693 0 : CPASSERT(nspins == 1)
694 : ! rhoin = (0.5 rho0, 0.5 rho0)
695 0 : CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, 2)
696 : ! rhoin = (0.5 rho0 + 0.5 rho1, 0.5 rho0)
697 0 : CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, 0.5_dp*beta)
698 : CALL qs_vxc_create(ks_env=ks_env, rho_struct=rhoin, xc_section=xc_section, &
699 0 : vxc_rho=vxc00, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE.)
700 0 : CALL pw_axpy(vxc00(2), vxc00(1), -1.0_dp)
701 : ELSE
702 0 : CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
703 0 : CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, beta)
704 : CALL qs_vxc_create(ks_env=ks_env, rho_struct=rhoin, xc_section=xc_section, &
705 0 : vxc_rho=vxc00, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE.)
706 : END IF
707 0 : CALL qs_rho_release(rhoin)
708 0 : DEALLOCATE (rhoin)
709 0 : IF (.NOT. ASSOCIATED(fxc_rho)) THEN
710 0 : ALLOCATE (fxc_rho(nspins))
711 0 : DO ispin = 1, nspins
712 0 : CALL auxbas_pw_pool%create_pw(fxc_rho(ispin))
713 0 : CALL pw_zero(fxc_rho(ispin))
714 : END DO
715 : END IF
716 0 : IF (.NOT. ASSOCIATED(gxc_rho)) THEN
717 0 : ALLOCATE (gxc_rho(nspins))
718 0 : DO ispin = 1, nspins
719 0 : CALL auxbas_pw_pool%create_pw(gxc_rho(ispin))
720 0 : CALL pw_zero(gxc_rho(ispin))
721 : END DO
722 : END IF
723 0 : CPASSERT(.NOT. ASSOCIATED(v_tau_rspace))
724 0 : DO ispin = 1, nspins
725 0 : IF (ak(istep) /= 0.0_dp) THEN
726 0 : CALL pw_axpy(vxc00(ispin), fxc_rho(ispin), ak(istep))
727 : END IF
728 0 : IF (bl(istep) /= 0.0_dp) THEN
729 0 : CALL pw_axpy(vxc00(ispin), gxc_rho(ispin), bl(istep))
730 : END IF
731 : END DO
732 0 : DO ispin = 1, SIZE(vxc00)
733 0 : CALL auxbas_pw_pool%give_back_pw(vxc00(ispin))
734 : END DO
735 0 : DEALLOCATE (vxc00)
736 :
737 : END DO
738 :
739 0 : oeps1 = 1.0_dp/epsrho
740 0 : oeps2 = 1.0_dp/(epsrho**2)
741 0 : DO ispin = 1, nspins
742 0 : CALL pw_scale(fxc_rho(ispin), oeps1)
743 0 : CALL pw_scale(gxc_rho(ispin), oeps2)
744 : END DO
745 :
746 0 : CALL timestop(handle)
747 :
748 0 : END SUBROUTINE qs_fgxc_create
749 :
750 : ! **************************************************************************************************
751 : !> \brief ...
752 : !> \param ks_env ...
753 : !> \param fxc_rho ...
754 : !> \param fxc_tau ...
755 : !> \param gxc_rho ...
756 : !> \param gxc_tau ...
757 : ! **************************************************************************************************
758 334 : SUBROUTINE qs_fgxc_release(ks_env, fxc_rho, fxc_tau, gxc_rho, gxc_tau)
759 :
760 : TYPE(qs_ks_env_type), POINTER :: ks_env
761 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau, gxc_rho, gxc_tau
762 :
763 : INTEGER :: ispin
764 : TYPE(pw_env_type), POINTER :: pw_env
765 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
766 :
767 334 : CALL get_ks_env(ks_env, pw_env=pw_env)
768 334 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
769 :
770 334 : IF (ASSOCIATED(fxc_rho)) THEN
771 704 : DO ispin = 1, SIZE(fxc_rho)
772 704 : CALL auxbas_pw_pool%give_back_pw(fxc_rho(ispin))
773 : END DO
774 334 : DEALLOCATE (fxc_rho)
775 : END IF
776 334 : IF (ASSOCIATED(fxc_tau)) THEN
777 0 : DO ispin = 1, SIZE(fxc_tau)
778 0 : CALL auxbas_pw_pool%give_back_pw(fxc_tau(ispin))
779 : END DO
780 0 : DEALLOCATE (fxc_tau)
781 : END IF
782 334 : IF (ASSOCIATED(gxc_rho)) THEN
783 712 : DO ispin = 1, SIZE(gxc_rho)
784 712 : CALL auxbas_pw_pool%give_back_pw(gxc_rho(ispin))
785 : END DO
786 334 : DEALLOCATE (gxc_rho)
787 : END IF
788 334 : IF (ASSOCIATED(gxc_tau)) THEN
789 0 : DO ispin = 1, SIZE(gxc_tau)
790 0 : CALL auxbas_pw_pool%give_back_pw(gxc_tau(ispin))
791 : END DO
792 0 : DEALLOCATE (gxc_tau)
793 : END IF
794 :
795 334 : END SUBROUTINE qs_fgxc_release
796 :
797 : END MODULE qs_fxc
|
