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