Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 auxbas_pw_pool ...
88 : !> \param is_triplet ...
89 : !> \param v_xc ...
90 : !> \param v_xc_tau ...
91 : !> \param spinflip ...
92 : ! **************************************************************************************************
93 17248 : SUBROUTINE qs_fxc_analytic(rho0, rho1_r, tau1_r, xc_section, auxbas_pw_pool, is_triplet, v_xc, v_xc_tau, spinflip)
94 :
95 : TYPE(qs_rho_type), POINTER :: rho0
96 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r, tau1_r
97 : TYPE(section_vals_type), POINTER :: xc_section
98 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
99 : LOGICAL, INTENT(IN) :: is_triplet
100 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_xc, v_xc_tau
101 : LOGICAL, OPTIONAL :: spinflip
102 :
103 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_fxc_analytic'
104 :
105 : INTEGER :: handle, nspins
106 : INTEGER, DIMENSION(2, 3) :: bo
107 : LOGICAL :: do_sf, lsd
108 : REAL(KIND=dp) :: fac
109 17248 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho0_g, rho1_g
110 17248 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho0_r, tau0_r
111 : TYPE(section_vals_type), POINTER :: xc_fun_section
112 : TYPE(xc_derivative_set_type) :: deriv_set
113 : TYPE(xc_rho_cflags_type) :: needs
114 : TYPE(xc_rho_set_type) :: rho0_set
115 :
116 8624 : CALL timeset(routineN, handle)
117 :
118 8624 : CPASSERT(.NOT. ASSOCIATED(v_xc))
119 8624 : CPASSERT(.NOT. ASSOCIATED(v_xc_tau))
120 :
121 8624 : do_sf = .FALSE.
122 8624 : IF (PRESENT(spinflip)) do_sf = spinflip
123 :
124 8624 : CALL qs_rho_get(rho0, rho_r=rho0_r, rho_g=rho0_g, tau_r=tau0_r)
125 8624 : nspins = SIZE(rho0_r)
126 :
127 8624 : lsd = (nspins == 2)
128 : fac = 0._dp
129 8624 : IF (is_triplet .AND. nspins == 1) fac = -1.0_dp
130 :
131 8624 : NULLIFY (rho1_g)
132 86240 : bo = rho1_r(1)%pw_grid%bounds_local
133 8624 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
134 8624 : needs = xc_functionals_get_needs(xc_fun_section, lsd, .TRUE.)
135 : ! calculate the arguments needed by the functionals and the values of the functional on the grid
136 8624 : CALL xc_prep_2nd_deriv(deriv_set, rho0_set, rho0_r, auxbas_pw_pool, xc_section=xc_section, tau_r=tau0_r)
137 : ! Folds the density rho1 with the functional
138 : CALL xc_calc_2nd_deriv(v_xc, v_xc_tau, deriv_set, rho0_set, rho1_r, rho1_g, tau1_r, &
139 : auxbas_pw_pool, xc_section=xc_section, gapw=.FALSE., do_triplet=is_triplet, &
140 8624 : do_sf=do_sf)
141 8624 : CALL xc_dset_release(deriv_set)
142 8624 : CALL xc_rho_set_release(rho0_set)
143 :
144 8624 : CALL timestop(handle)
145 :
146 189728 : END SUBROUTINE qs_fxc_analytic
147 :
148 : ! **************************************************************************************************
149 : !> \brief ...
150 : !> \param ks_env ...
151 : !> \param rho0_struct ...
152 : !> \param rho1_struct ...
153 : !> \param xc_section ...
154 : !> \param accuracy ...
155 : !> \param is_triplet ...
156 : !> \param fxc_rho ...
157 : !> \param fxc_tau ...
158 : ! **************************************************************************************************
159 1908 : SUBROUTINE qs_fxc_fdiff(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, is_triplet, &
160 : fxc_rho, fxc_tau)
161 :
162 : TYPE(qs_ks_env_type), POINTER :: ks_env
163 : TYPE(qs_rho_type), POINTER :: rho0_struct, rho1_struct
164 : TYPE(section_vals_type), POINTER :: xc_section
165 : INTEGER, INTENT(IN) :: accuracy
166 : LOGICAL, INTENT(IN) :: is_triplet
167 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau
168 :
169 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_fxc_fdiff'
170 : REAL(KIND=dp), PARAMETER :: epsrho = 5.e-4_dp
171 :
172 : INTEGER :: handle, ispin, istep, nspins, nstep
173 : REAL(KIND=dp) :: alpha, beta, exc, oeps1
174 : REAL(KIND=dp), DIMENSION(-4:4) :: ak
175 : TYPE(dft_control_type), POINTER :: dft_control
176 : TYPE(pw_env_type), POINTER :: pw_env
177 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
178 1908 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_tau_rspace, vxc00
179 : TYPE(qs_rho_type), POINTER :: rhoin
180 :
181 1908 : CALL timeset(routineN, handle)
182 :
183 1908 : CPASSERT(.NOT. ASSOCIATED(fxc_rho))
184 1908 : CPASSERT(.NOT. ASSOCIATED(fxc_tau))
185 1908 : CPASSERT(ASSOCIATED(rho0_struct))
186 1908 : CPASSERT(ASSOCIATED(rho1_struct))
187 :
188 1908 : ak = 0.0_dp
189 1908 : SELECT CASE (accuracy)
190 : CASE (:4)
191 0 : nstep = 2
192 0 : ak(-2:2) = [1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp]/12.0_dp
193 : CASE (5:7)
194 15264 : nstep = 3
195 15264 : 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
196 : CASE (8:)
197 0 : nstep = 4
198 : ak(-4:4) = [1.0_dp, -32.0_dp/3.0_dp, 56.0_dp, -224.0_dp, 0.0_dp, &
199 1908 : 224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp]/280.0_dp
200 : END SELECT
201 :
202 1908 : CALL get_ks_env(ks_env, dft_control=dft_control, pw_env=pw_env)
203 1908 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
204 :
205 1908 : nspins = dft_control%nspins
206 1908 : exc = 0.0_dp
207 :
208 15264 : DO istep = -nstep, nstep
209 :
210 15264 : IF (ak(istep) /= 0.0_dp) THEN
211 11448 : alpha = 1.0_dp
212 11448 : beta = REAL(istep, KIND=dp)*epsrho
213 : NULLIFY (rhoin)
214 11448 : ALLOCATE (rhoin)
215 11448 : CALL qs_rho_create(rhoin)
216 11448 : NULLIFY (vxc00, v_tau_rspace)
217 11448 : IF (is_triplet) THEN
218 924 : CPASSERT(nspins == 1)
219 : ! rhoin = (0.5 rho0, 0.5 rho0)
220 924 : CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, 2)
221 : ! rhoin = (0.5 rho0 + 0.5 rho1, 0.5 rho0)
222 924 : CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, 0.5_dp*beta)
223 : CALL qs_vxc_create(ks_env=ks_env, rho_struct=rhoin, xc_section=xc_section, &
224 924 : vxc_rho=vxc00, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE.)
225 924 : CALL pw_axpy(vxc00(2), vxc00(1), -1.0_dp)
226 924 : IF (ASSOCIATED(v_tau_rspace)) CALL pw_axpy(v_tau_rspace(2), v_tau_rspace(1), -1.0_dp)
227 : ELSE
228 10524 : CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
229 10524 : CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, beta)
230 : CALL qs_vxc_create(ks_env=ks_env, rho_struct=rhoin, xc_section=xc_section, &
231 10524 : vxc_rho=vxc00, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE.)
232 : END IF
233 11448 : CALL qs_rho_release(rhoin)
234 11448 : DEALLOCATE (rhoin)
235 11448 : IF (.NOT. ASSOCIATED(fxc_rho)) THEN
236 7932 : ALLOCATE (fxc_rho(nspins))
237 4116 : DO ispin = 1, nspins
238 2208 : CALL auxbas_pw_pool%create_pw(fxc_rho(ispin))
239 4116 : CALL pw_zero(fxc_rho(ispin))
240 : END DO
241 : END IF
242 24696 : DO ispin = 1, nspins
243 24696 : CALL pw_axpy(vxc00(ispin), fxc_rho(ispin), ak(istep))
244 : END DO
245 25620 : DO ispin = 1, SIZE(vxc00)
246 25620 : CALL auxbas_pw_pool%give_back_pw(vxc00(ispin))
247 : END DO
248 11448 : DEALLOCATE (vxc00)
249 11448 : IF (ASSOCIATED(v_tau_rspace)) THEN
250 0 : IF (.NOT. ASSOCIATED(fxc_tau)) THEN
251 0 : ALLOCATE (fxc_tau(nspins))
252 0 : DO ispin = 1, nspins
253 0 : CALL auxbas_pw_pool%create_pw(fxc_tau(ispin))
254 0 : CALL pw_zero(fxc_tau(ispin))
255 : END DO
256 : END IF
257 0 : DO ispin = 1, nspins
258 0 : CALL pw_axpy(v_tau_rspace(ispin), fxc_tau(ispin), ak(istep))
259 : END DO
260 0 : DO ispin = 1, SIZE(v_tau_rspace)
261 0 : CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
262 : END DO
263 0 : DEALLOCATE (v_tau_rspace)
264 : END IF
265 : END IF
266 :
267 : END DO
268 :
269 1908 : oeps1 = 1.0_dp/epsrho
270 4116 : DO ispin = 1, nspins
271 4116 : CALL pw_scale(fxc_rho(ispin), oeps1)
272 : END DO
273 1908 : IF (ASSOCIATED(fxc_tau)) THEN
274 0 : DO ispin = 1, nspins
275 0 : CALL pw_scale(fxc_tau(ispin), oeps1)
276 : END DO
277 : END IF
278 :
279 1908 : CALL timestop(handle)
280 :
281 1908 : END SUBROUTINE qs_fxc_fdiff
282 :
283 : ! **************************************************************************************************
284 : !> \brief Calculates the values at the grid points in real space (r_i), of the second and third
285 : !> functional derivatives of the exchange-correlation energy functional.
286 : !> fxc_rho(r_i) = fxc[n](r_i)*n^(1)(r_i) ! Second functional derivative
287 : !> gxc_rho(r_i) = n^(1)(r_i)*gxc[n](r_i)*n^(1)(r_i) ! Third functional derivative
288 : !> \param rho0_struct Ground state density, n(r).
289 : !> \param rho1_struct Density used to fold the functional derivatives, n^(1)(r).
290 : !> \param xc_section ...
291 : !> \param pw_pool ...
292 : !> \param fxc_rho Second functional derivative with respect to the density, n(r).
293 : !> \param fxc_tau mGGA contribution to the second functional derivative with respect to the density.
294 : !> \param gxc_rho Third functional derivative with respect to the density, n(r).
295 : !> \param gxc_tau mGGA contribution to the third functional derivative with respect to the density.
296 : !> \param spinflip Flag used to activate the spin-flip noncollinear kernel and kernel derivatives.
297 : !> \par History
298 : !> * 07.2024 Created [LHS]
299 : ! **************************************************************************************************
300 0 : SUBROUTINE qs_fgxc_analytic(rho0_struct, rho1_struct, xc_section, pw_pool, &
301 : fxc_rho, fxc_tau, gxc_rho, gxc_tau, spinflip)
302 :
303 : TYPE(qs_rho_type), POINTER :: rho0_struct, rho1_struct
304 : TYPE(section_vals_type), POINTER :: xc_section
305 : TYPE(pw_pool_type), POINTER :: pw_pool
306 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau, gxc_rho, gxc_tau
307 : LOGICAL, INTENT(IN), OPTIONAL :: spinflip
308 :
309 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_fgxc_analytic'
310 :
311 : INTEGER :: handle, ispin, nspins, spindim
312 : INTEGER, DIMENSION(2, 3) :: bo
313 : LOGICAL :: do_sf, lsd
314 : REAL(KIND=dp) :: fac
315 0 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho0_g, rho1_g
316 0 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho0_r, rho1_r, tau0_r, tau1_r
317 : TYPE(section_vals_type), POINTER :: xc_fun_section
318 : TYPE(xc_derivative_set_type) :: deriv_set
319 : TYPE(xc_rho_cflags_type) :: needs
320 : TYPE(xc_rho_set_type) :: rho0_set, rho1_set
321 :
322 0 : CALL timeset(routineN, handle)
323 :
324 : ! Only rho0 and rho1 should be associated
325 0 : CPASSERT(.NOT. ASSOCIATED(fxc_rho))
326 0 : CPASSERT(.NOT. ASSOCIATED(fxc_tau))
327 0 : CPASSERT(.NOT. ASSOCIATED(gxc_rho))
328 0 : CPASSERT(.NOT. ASSOCIATED(gxc_tau))
329 0 : CPASSERT(ASSOCIATED(rho0_struct))
330 0 : CPASSERT(ASSOCIATED(rho1_struct))
331 :
332 : ! Initialize parameters
333 0 : do_sf = .FALSE.
334 0 : IF (PRESENT(spinflip)) do_sf = spinflip
335 : !
336 : ! Get the values on the gridpoints of the rho0 density
337 0 : CALL qs_rho_get(rho0_struct, rho_r=rho0_r, rho_g=rho0_g, tau_r=tau0_r)
338 0 : nspins = SIZE(rho0_r)
339 0 : lsd = (nspins == 2)
340 : !
341 0 : IF (do_sf) THEN
342 : spindim = 1
343 : ELSE
344 0 : spindim = nspins
345 : END IF
346 : !
347 0 : fac = 0._dp
348 0 : IF (nspins == 1) THEN
349 0 : fac = 1.0_dp
350 : END IF
351 :
352 : ! Read xc functional section and find out what the functional actually needs
353 0 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
354 0 : needs = xc_functionals_get_needs(xc_fun_section, lsd, .TRUE.)
355 :
356 : ! Create fields for the kernel and kernel derivative
357 0 : ALLOCATE (fxc_rho(spindim), gxc_rho(nspins))
358 0 : DO ispin = 1, spindim
359 0 : CALL pw_pool%create_pw(fxc_rho(ispin))
360 0 : CALL pw_zero(fxc_rho(ispin))
361 : END DO
362 0 : DO ispin = 1, nspins
363 0 : CALL pw_pool%create_pw(gxc_rho(ispin))
364 0 : CALL pw_zero(gxc_rho(ispin))
365 : END DO
366 : ! Create fields for mGGA functionals. This implementation is not ready yet!
367 0 : IF (needs%tau .OR. needs%tau_spin) THEN
368 0 : IF (.NOT. ASSOCIATED(tau1_r)) &
369 0 : CPABORT("Tau-dependent functionals requires allocated kinetic energy density grid")
370 0 : ALLOCATE (fxc_tau(spindim), gxc_tau(nspins))
371 0 : DO ispin = 1, spindim
372 0 : CALL pw_pool%create_pw(fxc_tau(ispin))
373 0 : CALL pw_zero(fxc_tau(ispin))
374 : END DO
375 0 : DO ispin = 1, nspins
376 0 : CALL pw_pool%create_pw(gxc_tau(ispin))
377 0 : CALL pw_zero(gxc_tau(ispin))
378 : END DO
379 : END IF
380 :
381 : ! Build rho0_set
382 : ! calculate the arguments needed by the functionals
383 : ! Needs
384 : ! deriv_set xc_derivative_set_type just declared
385 : ! rho0_set xc_rho_set_type just declared
386 : ! rho0_r pw_type calculated by qs_rho_get
387 : ! pw_pool given by the calling subroutine
388 : ! xc_section given by the calling subroutine
389 : ! tau0_r pw_type calculated by qs_rho_get
390 : CALL xc_prep_3rd_deriv(deriv_set, rho0_set, rho0_r, pw_pool, xc_section, tau_r=tau0_r, &
391 0 : do_sf=do_sf)
392 :
393 : ! Build rho1_set
394 : ! Get the values on the gridpoints of the rho1 density
395 0 : CALL qs_rho_get(rho1_struct, rho_r=rho1_r, rho_g=rho1_g, tau_r=tau1_r)
396 0 : bo = rho1_r(1)%pw_grid%bounds_local
397 : ! create the place where to store the argument for the functionals
398 : ! Needs
399 : ! rho1_set xc_rho_set_type just declared
400 : ! bo 2x3 integer matrix should have bounds_local or rho1_r
401 : CALL xc_rho_set_create(rho1_set, bo, &
402 : rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
403 : drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
404 0 : tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
405 : ! calculate the arguments needed by the functionals
406 : ! Needs
407 : ! rho1_set object created by xc_rho_set_create
408 : ! rho1_r,rho1_g,tau1_r pw_type values of rho1 in real space grid
409 : ! needs xc_rho_cflags_type defined through xc_functionals_get_needs
410 : ! pw_pool Given by the calling subroutine
411 : CALL xc_rho_set_update(rho1_set, rho1_r, rho1_g, tau1_r, needs, &
412 : section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
413 : section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
414 0 : pw_pool, spinflip=do_sf)
415 :
416 : ! Calculate exchange correlation kernel
417 : ! Needs
418 : ! fxc_rho, fxc_tau pw_type not associated
419 : ! deriv_set created and defined by xc_prep_3rd_deriv
420 : ! rho0_set xc_rho_set_type build by xc_prep_3rd_deriv
421 : ! rho1_set xc_rho_set_type build by xc_rho_set_create/update
422 : ! pw_pool needs to be given by the calling subroutine
423 : ! xc_section needs to be given by the calling subroutine
424 : CALL xc_calc_2nd_deriv_analytical(fxc_rho, fxc_tau, deriv_set, rho0_set, rho1_set, pw_pool, &
425 0 : xc_section, .FALSE., spinflip=do_sf, tddfpt_fac=fac)
426 : ! Calculate exchange correlation kernel derivative
427 : CALL xc_calc_3rd_deriv_analytical(gxc_rho, gxc_tau, deriv_set, rho0_set, rho1_set, pw_pool, &
428 0 : xc_section, spinflip=do_sf)
429 :
430 0 : CALL xc_dset_release(deriv_set)
431 0 : CALL xc_rho_set_release(rho0_set)
432 0 : CALL xc_rho_set_release(rho1_set)
433 :
434 0 : CALL timestop(handle)
435 :
436 0 : END SUBROUTINE qs_fgxc_analytic
437 :
438 : ! **************************************************************************************************
439 : !> \brief ...
440 : !> \param ks_env ...
441 : !> \param rho0_struct ...
442 : !> \param rho1_struct ...
443 : !> \param xc_section ...
444 : !> \param accuracy ...
445 : !> \param epsrho ...
446 : !> \param is_triplet ...
447 : !> \param fxc_rho ...
448 : !> \param fxc_tau ...
449 : !> \param gxc_rho ...
450 : !> \param gxc_tau ...
451 : !> \param spinflip ...
452 : ! **************************************************************************************************
453 272 : SUBROUTINE qs_fgxc_gdiff(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, epsrho, &
454 : is_triplet, fxc_rho, fxc_tau, gxc_rho, gxc_tau, spinflip)
455 :
456 : TYPE(qs_ks_env_type), POINTER :: ks_env
457 : TYPE(qs_rho_type), POINTER :: rho0_struct, rho1_struct
458 : TYPE(section_vals_type), POINTER :: xc_section
459 : INTEGER, INTENT(IN) :: accuracy
460 : REAL(KIND=dp), INTENT(IN) :: epsrho
461 : LOGICAL, INTENT(IN) :: is_triplet
462 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau, gxc_rho, gxc_tau
463 : LOGICAL, OPTIONAL :: spinflip
464 :
465 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_fgxc_gdiff'
466 :
467 : INTEGER :: handle, ispin, istep, nspins, nstep
468 : LOGICAL :: do_sf
469 : REAL(KIND=dp) :: alpha, beta, exc, oeps1
470 : REAL(KIND=dp), DIMENSION(-4:4) :: ak
471 : TYPE(dft_control_type), POINTER :: dft_control
472 : TYPE(pw_env_type), POINTER :: pw_env
473 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
474 272 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r, tau1_r, v_tau_rspace, vxc00, &
475 272 : vxc00b
476 : TYPE(qs_rho_type), POINTER :: rhoin
477 :
478 272 : CALL timeset(routineN, handle)
479 :
480 272 : CPASSERT(.NOT. ASSOCIATED(fxc_rho))
481 272 : CPASSERT(.NOT. ASSOCIATED(fxc_tau))
482 272 : CPASSERT(.NOT. ASSOCIATED(gxc_rho))
483 272 : CPASSERT(.NOT. ASSOCIATED(gxc_tau))
484 272 : CPASSERT(ASSOCIATED(rho0_struct))
485 272 : CPASSERT(ASSOCIATED(rho1_struct))
486 :
487 272 : do_sf = .FALSE.
488 272 : IF (PRESENT(spinflip)) do_sf = spinflip
489 :
490 272 : ak = 0.0_dp
491 272 : SELECT CASE (accuracy)
492 : CASE (:4)
493 0 : nstep = 2
494 0 : ak(-2:2) = [1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp]/12.0_dp
495 : CASE (5:7)
496 2176 : nstep = 3
497 2176 : 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
498 : CASE (8:)
499 0 : nstep = 4
500 : ak(-4:4) = [1.0_dp, -32.0_dp/3.0_dp, 56.0_dp, -224.0_dp, 0.0_dp, &
501 272 : 224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp]/280.0_dp
502 : END SELECT
503 :
504 272 : CALL get_ks_env(ks_env, dft_control=dft_control, pw_env=pw_env)
505 272 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
506 :
507 272 : nspins = dft_control%nspins
508 272 : exc = 0.0_dp
509 :
510 272 : IF (do_sf) THEN
511 8 : CALL qs_rho_get(rho1_struct, rho_r=rho1_r, tau_r=tau1_r)
512 : CALL qs_fxc_analytic(rho0_struct, rho1_r, tau1_r, xc_section, &
513 8 : auxbas_pw_pool, is_triplet, fxc_rho, fxc_tau, spinflip=do_sf)
514 : ELSE
515 : CALL qs_fxc_fdiff(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, is_triplet, &
516 264 : fxc_rho, fxc_tau)
517 : END IF
518 :
519 2176 : DO istep = -nstep, nstep
520 :
521 2176 : IF (ak(istep) /= 0.0_dp) THEN
522 1632 : alpha = 1.0_dp
523 1632 : beta = REAL(istep, KIND=dp)*epsrho
524 : NULLIFY (rhoin)
525 1632 : ALLOCATE (rhoin)
526 1632 : CALL qs_rho_create(rhoin)
527 1632 : NULLIFY (vxc00, vxc00b, v_tau_rspace)
528 1632 : CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
529 1632 : CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, beta)
530 1632 : IF (do_sf) THEN
531 : ! variation in alpha density
532 : CALL qs_fxc_analytic(rhoin, rho1_r, tau1_r, &
533 : xc_section, auxbas_pw_pool, is_triplet, &
534 48 : vxc00, v_tau_rspace, spinflip=do_sf)
535 : ! variation in beta density
536 48 : CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
537 48 : CALL qs_rho_scale_and_add_b(rhoin, rho1_struct, alpha, beta)
538 : CALL qs_fxc_analytic(rhoin, rho1_r, tau1_r, &
539 : xc_section, auxbas_pw_pool, is_triplet, &
540 48 : vxc00b, v_tau_rspace, spinflip=do_sf)
541 : ELSE
542 : CALL qs_fxc_fdiff(ks_env=ks_env, rho0_struct=rhoin, rho1_struct=rho1_struct, &
543 : xc_section=xc_section, accuracy=accuracy, is_triplet=is_triplet, &
544 1584 : fxc_rho=vxc00, fxc_tau=v_tau_rspace)
545 : END IF
546 1632 : CALL qs_rho_release(rhoin)
547 1632 : DEALLOCATE (rhoin)
548 1632 : IF (.NOT. ASSOCIATED(gxc_rho)) THEN
549 1132 : ALLOCATE (gxc_rho(nspins))
550 588 : DO ispin = 1, nspins
551 316 : CALL auxbas_pw_pool%create_pw(gxc_rho(ispin))
552 588 : CALL pw_zero(gxc_rho(ispin))
553 : END DO
554 : END IF
555 1632 : IF (do_sf) THEN
556 48 : CALL pw_axpy(vxc00(1), gxc_rho(1), ak(istep))
557 48 : CALL pw_axpy(vxc00b(1), gxc_rho(2), ak(istep))
558 : ELSE
559 3384 : DO ispin = 1, nspins
560 3384 : CALL pw_axpy(vxc00(ispin), gxc_rho(ispin), ak(istep))
561 : END DO
562 : END IF
563 3480 : DO ispin = 1, SIZE(vxc00)
564 3480 : CALL auxbas_pw_pool%give_back_pw(vxc00(ispin))
565 : END DO
566 1632 : DEALLOCATE (vxc00)
567 1632 : IF (ASSOCIATED(vxc00b)) THEN
568 48 : CALL auxbas_pw_pool%give_back_pw(vxc00b(1))
569 48 : DEALLOCATE (vxc00b)
570 : END IF
571 1632 : IF (ASSOCIATED(v_tau_rspace)) THEN
572 0 : IF (.NOT. ASSOCIATED(gxc_tau)) THEN
573 0 : ALLOCATE (gxc_tau(nspins))
574 0 : DO ispin = 1, nspins
575 0 : CALL auxbas_pw_pool%create_pw(gxc_tau(ispin))
576 0 : CALL pw_zero(gxc_tau(ispin))
577 : END DO
578 : END IF
579 0 : DO ispin = 1, nspins
580 0 : CALL pw_axpy(v_tau_rspace(ispin), gxc_tau(ispin), ak(istep))
581 : END DO
582 0 : DO ispin = 1, SIZE(v_tau_rspace)
583 0 : CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
584 : END DO
585 0 : DEALLOCATE (v_tau_rspace)
586 : END IF
587 : END IF
588 :
589 : END DO
590 :
591 272 : oeps1 = 1.0_dp/epsrho
592 588 : DO ispin = 1, nspins
593 588 : CALL pw_scale(gxc_rho(ispin), oeps1)
594 : END DO
595 272 : IF (ASSOCIATED(gxc_tau)) THEN
596 0 : DO ispin = 1, nspins
597 0 : CALL pw_scale(gxc_tau(ispin), oeps1)
598 : END DO
599 : END IF
600 :
601 272 : CALL timestop(handle)
602 :
603 272 : END SUBROUTINE qs_fgxc_gdiff
604 :
605 : ! **************************************************************************************************
606 : !> \brief ...
607 : !> \param ks_env ...
608 : !> \param rho0_struct ...
609 : !> \param rho1_struct ...
610 : !> \param xc_section ...
611 : !> \param accuracy ...
612 : !> \param is_triplet ...
613 : !> \param fxc_rho ...
614 : !> \param fxc_tau ...
615 : !> \param gxc_rho ...
616 : !> \param gxc_tau ...
617 : ! **************************************************************************************************
618 0 : SUBROUTINE qs_fgxc_create(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, is_triplet, &
619 : fxc_rho, fxc_tau, gxc_rho, gxc_tau)
620 :
621 : TYPE(qs_ks_env_type), POINTER :: ks_env
622 : TYPE(qs_rho_type), POINTER :: rho0_struct, rho1_struct
623 : TYPE(section_vals_type), POINTER :: xc_section
624 : INTEGER, INTENT(IN) :: accuracy
625 : LOGICAL, INTENT(IN) :: is_triplet
626 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau, gxc_rho, gxc_tau
627 :
628 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_fgxc_create'
629 : REAL(KIND=dp), PARAMETER :: epsrho = 5.e-4_dp
630 :
631 : INTEGER :: handle, ispin, istep, nspins, nstep
632 : REAL(KIND=dp) :: alpha, beta, exc, oeps1, oeps2
633 : REAL(KIND=dp), DIMENSION(-4:4) :: ak, bl
634 : TYPE(dft_control_type), POINTER :: dft_control
635 : TYPE(pw_env_type), POINTER :: pw_env
636 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
637 0 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_tau_rspace, vxc00
638 : TYPE(qs_rho_type), POINTER :: rhoin
639 :
640 0 : CALL timeset(routineN, handle)
641 :
642 0 : CPASSERT(.NOT. ASSOCIATED(fxc_rho))
643 0 : CPASSERT(.NOT. ASSOCIATED(fxc_tau))
644 0 : CPASSERT(.NOT. ASSOCIATED(gxc_rho))
645 0 : CPASSERT(.NOT. ASSOCIATED(gxc_tau))
646 0 : CPASSERT(ASSOCIATED(rho0_struct))
647 0 : CPASSERT(ASSOCIATED(rho1_struct))
648 :
649 0 : ak = 0.0_dp
650 0 : bl = 0.0_dp
651 0 : SELECT CASE (accuracy)
652 : CASE (:4)
653 0 : nstep = 2
654 0 : ak(-2:2) = [1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp]/12.0_dp
655 0 : bl(-2:2) = [-1.0_dp, 16.0_dp, -30.0_dp, 16.0_dp, -1.0_dp]/12.0_dp
656 : CASE (5:7)
657 0 : nstep = 3
658 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
659 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
660 : CASE (8:)
661 0 : nstep = 4
662 : ak(-4:4) = [1.0_dp, -32.0_dp/3.0_dp, 56.0_dp, -224.0_dp, 0.0_dp, &
663 0 : 224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp]/280.0_dp
664 : bl(-4:4) = [-1.0_dp, 128.0_dp/9.0_dp, -112.0_dp, 896.0_dp, -14350.0_dp/9.0_dp, &
665 0 : 896.0_dp, -112.0_dp, 128.0_dp/9.0_dp, -1.0_dp]/560.0_dp
666 : END SELECT
667 :
668 0 : CALL get_ks_env(ks_env, dft_control=dft_control, pw_env=pw_env)
669 0 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
670 :
671 0 : nspins = dft_control%nspins
672 0 : exc = 0.0_dp
673 :
674 0 : DO istep = -nstep, nstep
675 :
676 0 : alpha = 1.0_dp
677 0 : beta = REAL(istep, KIND=dp)*epsrho
678 : NULLIFY (rhoin)
679 0 : ALLOCATE (rhoin)
680 0 : CALL qs_rho_create(rhoin)
681 0 : NULLIFY (vxc00, v_tau_rspace)
682 0 : IF (is_triplet) THEN
683 0 : CPASSERT(nspins == 1)
684 : ! rhoin = (0.5 rho0, 0.5 rho0)
685 0 : CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, 2)
686 : ! rhoin = (0.5 rho0 + 0.5 rho1, 0.5 rho0)
687 0 : CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, 0.5_dp*beta)
688 : CALL qs_vxc_create(ks_env=ks_env, rho_struct=rhoin, xc_section=xc_section, &
689 0 : vxc_rho=vxc00, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE.)
690 0 : CALL pw_axpy(vxc00(2), vxc00(1), -1.0_dp)
691 : ELSE
692 0 : CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
693 0 : CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, beta)
694 : CALL qs_vxc_create(ks_env=ks_env, rho_struct=rhoin, xc_section=xc_section, &
695 0 : vxc_rho=vxc00, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE.)
696 : END IF
697 0 : CALL qs_rho_release(rhoin)
698 0 : DEALLOCATE (rhoin)
699 0 : IF (.NOT. ASSOCIATED(fxc_rho)) THEN
700 0 : ALLOCATE (fxc_rho(nspins))
701 0 : DO ispin = 1, nspins
702 0 : CALL auxbas_pw_pool%create_pw(fxc_rho(ispin))
703 0 : CALL pw_zero(fxc_rho(ispin))
704 : END DO
705 : END IF
706 0 : IF (.NOT. ASSOCIATED(gxc_rho)) THEN
707 0 : ALLOCATE (gxc_rho(nspins))
708 0 : DO ispin = 1, nspins
709 0 : CALL auxbas_pw_pool%create_pw(gxc_rho(ispin))
710 0 : CALL pw_zero(gxc_rho(ispin))
711 : END DO
712 : END IF
713 0 : CPASSERT(.NOT. ASSOCIATED(v_tau_rspace))
714 0 : DO ispin = 1, nspins
715 0 : IF (ak(istep) /= 0.0_dp) THEN
716 0 : CALL pw_axpy(vxc00(ispin), fxc_rho(ispin), ak(istep))
717 : END IF
718 0 : IF (bl(istep) /= 0.0_dp) THEN
719 0 : CALL pw_axpy(vxc00(ispin), gxc_rho(ispin), bl(istep))
720 : END IF
721 : END DO
722 0 : DO ispin = 1, SIZE(vxc00)
723 0 : CALL auxbas_pw_pool%give_back_pw(vxc00(ispin))
724 : END DO
725 0 : DEALLOCATE (vxc00)
726 :
727 : END DO
728 :
729 0 : oeps1 = 1.0_dp/epsrho
730 0 : oeps2 = 1.0_dp/(epsrho**2)
731 0 : DO ispin = 1, nspins
732 0 : CALL pw_scale(fxc_rho(ispin), oeps1)
733 0 : CALL pw_scale(gxc_rho(ispin), oeps2)
734 : END DO
735 :
736 0 : CALL timestop(handle)
737 :
738 0 : END SUBROUTINE qs_fgxc_create
739 :
740 : ! **************************************************************************************************
741 : !> \brief ...
742 : !> \param ks_env ...
743 : !> \param fxc_rho ...
744 : !> \param fxc_tau ...
745 : !> \param gxc_rho ...
746 : !> \param gxc_tau ...
747 : ! **************************************************************************************************
748 272 : SUBROUTINE qs_fgxc_release(ks_env, fxc_rho, fxc_tau, gxc_rho, gxc_tau)
749 :
750 : TYPE(qs_ks_env_type), POINTER :: ks_env
751 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau, gxc_rho, gxc_tau
752 :
753 : INTEGER :: ispin
754 : TYPE(pw_env_type), POINTER :: pw_env
755 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
756 :
757 272 : CALL get_ks_env(ks_env, pw_env=pw_env)
758 272 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
759 :
760 272 : IF (ASSOCIATED(fxc_rho)) THEN
761 580 : DO ispin = 1, SIZE(fxc_rho)
762 580 : CALL auxbas_pw_pool%give_back_pw(fxc_rho(ispin))
763 : END DO
764 272 : DEALLOCATE (fxc_rho)
765 : END IF
766 272 : IF (ASSOCIATED(fxc_tau)) THEN
767 0 : DO ispin = 1, SIZE(fxc_tau)
768 0 : CALL auxbas_pw_pool%give_back_pw(fxc_tau(ispin))
769 : END DO
770 0 : DEALLOCATE (fxc_tau)
771 : END IF
772 272 : IF (ASSOCIATED(gxc_rho)) THEN
773 588 : DO ispin = 1, SIZE(gxc_rho)
774 588 : CALL auxbas_pw_pool%give_back_pw(gxc_rho(ispin))
775 : END DO
776 272 : DEALLOCATE (gxc_rho)
777 : END IF
778 272 : IF (ASSOCIATED(gxc_tau)) THEN
779 0 : DO ispin = 1, SIZE(gxc_tau)
780 0 : CALL auxbas_pw_pool%give_back_pw(gxc_tau(ispin))
781 : END DO
782 0 : DEALLOCATE (gxc_tau)
783 : END IF
784 :
785 272 : END SUBROUTINE qs_fgxc_release
786 :
787 : END MODULE qs_fxc
|
