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 : !> complete refactoring 08.2026
26 : !> \author JGH
27 : ! **************************************************************************************************
28 : MODULE qs_fgxc
29 :
30 : USE cp_control_types, ONLY: dft_control_type
31 : USE input_section_types, ONLY: section_get_ival,&
32 : section_get_lval,&
33 : section_get_rval,&
34 : section_vals_get_subs_vals,&
35 : section_vals_type
36 : USE kinds, ONLY: dp
37 : USE pw_env_types, ONLY: pw_env_get,&
38 : pw_env_type
39 : USE pw_grids, ONLY: pw_grid_compare
40 : USE pw_methods, ONLY: pw_axpy,&
41 : pw_scale,&
42 : pw_transfer,&
43 : pw_zero
44 : USE pw_pool_types, ONLY: pw_pool_type
45 : USE pw_types, ONLY: pw_c1d_gs_type,&
46 : pw_r3d_rs_type
47 : USE qs_environment_types, ONLY: get_qs_env,&
48 : qs_environment_type
49 : USE qs_fgxc_atom, ONLY: fgxc_atom_calc,&
50 : fgxc_atom_diff
51 : USE qs_fxc, ONLY: qs_fxc_create,&
52 : qs_fxc_fdiff
53 : USE qs_kind_types, ONLY: qs_kind_type
54 : USE qs_ks_types, ONLY: qs_ks_env_type
55 : USE qs_rho_atom_types, ONLY: rho_atom_type
56 : USE qs_rho_methods, ONLY: qs_rho_copy,&
57 : qs_rho_scale_and_add,&
58 : qs_rho_scale_and_add_b,&
59 : qs_rho_transfer
60 : USE qs_rho_types, ONLY: qs_rho_create,&
61 : qs_rho_get,&
62 : qs_rho_release,&
63 : qs_rho_type
64 : USE qs_vxc, ONLY: qs_vxc_create
65 : USE xc, ONLY: xc_calc_2nd_deriv_analytical,&
66 : xc_calc_3rd_deriv_analytical,&
67 : xc_prep_3rd_deriv
68 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
69 : xc_dset_release
70 : USE xc_derivatives, ONLY: xc_functionals_get_needs
71 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
72 : USE xc_rho_set_types, ONLY: xc_rho_set_create,&
73 : xc_rho_set_release,&
74 : xc_rho_set_type,&
75 : xc_rho_set_update
76 : #include "./base/base_uses.f90"
77 :
78 : IMPLICIT NONE
79 :
80 : PRIVATE
81 :
82 : ! *** Public subroutines ***
83 : PUBLIC :: qs_fgxc_create, qs_fgxc_release
84 :
85 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_fgxc'
86 :
87 : ! **************************************************************************************************
88 :
89 : CONTAINS
90 :
91 : ! **************************************************************************************************
92 : !> \brief ...
93 : !> \param qs_env ...
94 : !> \param rho0_struct ...
95 : !> \param rho1_struct ...
96 : !> \param rho0_atom_set ...
97 : !> \param xc_section ...
98 : !> \param do_onecenter ...
99 : !> \param is_triplet ...
100 : !> \param spinflip ...
101 : !> \param fxc_rho ...
102 : !> \param fxc_tau ...
103 : !> \param gxc_rho ...
104 : !> \param gxc_tau ...
105 : !> \param rhof_atom_set ...
106 : !> \param rhog_atom_set ...
107 : !> \param qs_kind_set ...
108 : ! **************************************************************************************************
109 334 : SUBROUTINE qs_fgxc_create(qs_env, rho0_struct, rho1_struct, rho0_atom_set, &
110 : xc_section, do_onecenter, is_triplet, spinflip, &
111 : fxc_rho, fxc_tau, gxc_rho, gxc_tau, &
112 : rhof_atom_set, rhog_atom_set, qs_kind_set)
113 :
114 : TYPE(qs_environment_type), POINTER :: qs_env
115 : TYPE(qs_rho_type), POINTER :: rho0_struct, rho1_struct
116 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho0_atom_set
117 : TYPE(section_vals_type), POINTER :: xc_section
118 : LOGICAL, INTENT(IN) :: do_onecenter, is_triplet, spinflip
119 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau, gxc_rho, gxc_tau
120 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rhof_atom_set, rhog_atom_set
121 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
122 :
123 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_fgxc_create'
124 :
125 : INTEGER :: handle, ispin, nspins, nsteps, order
126 : LOGICAL :: analytic_2nd_deriv, analytic_3rd_deriv, &
127 : has_rho, has_tau, uf_grid
128 : REAL(KIND=dp) :: eps_delta, factor
129 : TYPE(dft_control_type), POINTER :: dft_control
130 334 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
131 : TYPE(pw_c1d_gs_type), POINTER :: rho_nlcc_g
132 : TYPE(pw_env_type), POINTER :: pw_env
133 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool, uf_pw_pool, xc_pw_pool
134 334 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho_uf, fxc_tau_uf, gxc_rho_uf, &
135 334 : gxc_tau_uf, rho_r
136 : TYPE(pw_r3d_rs_type), POINTER :: rho_nlcc, weights, weights_uf
137 : TYPE(qs_rho_type), POINTER :: rho0_uf, rho1_uf
138 :
139 334 : CALL timeset(routineN, handle)
140 :
141 : ! Only rho0 and rho1 should be associated
142 334 : CPASSERT(.NOT. ASSOCIATED(fxc_rho))
143 334 : CPASSERT(.NOT. ASSOCIATED(fxc_tau))
144 334 : CPASSERT(.NOT. ASSOCIATED(gxc_rho))
145 334 : CPASSERT(.NOT. ASSOCIATED(gxc_tau))
146 334 : CPASSERT(ASSOCIATED(rho0_struct))
147 334 : CPASSERT(ASSOCIATED(rho1_struct))
148 :
149 334 : nsteps = section_get_ival(xc_section, "NSTEPS")
150 334 : order = 2*nsteps
151 334 : eps_delta = section_get_rval(xc_section, "STEP_SIZE")
152 334 : analytic_2nd_deriv = section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL")
153 334 : analytic_3rd_deriv = section_get_lval(xc_section, "3RD_DERIV_ANALYTICAL")
154 334 : IF (analytic_3rd_deriv .AND. .NOT. analytic_2nd_deriv) THEN
155 : CALL cp_warn(__LOCATION__, "XC Analytic 3rd derivatives cannot be combined with "// &
156 0 : "numeric 2nd derivatives.")
157 0 : CPABORT("XC derivative options")
158 : END IF
159 :
160 334 : CALL get_qs_env(qs_env, dft_control=dft_control)
161 334 : nspins = dft_control%nspins
162 :
163 334 : CALL get_qs_env(qs_env, pw_env=pw_env)
164 334 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, xc_pw_pool=xc_pw_pool)
165 334 : uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid, xc_pw_pool%pw_grid)
166 :
167 334 : NULLIFY (fxc_rho_uf, fxc_tau_uf, gxc_rho_uf, gxc_tau_uf)
168 334 : IF (analytic_2nd_deriv .AND. analytic_3rd_deriv) THEN
169 :
170 : ! NLCC
171 0 : NULLIFY (rho_nlcc, rho_nlcc_g)
172 0 : CALL get_qs_env(qs_env, rho_nlcc=rho_nlcc, rho_nlcc_g=rho_nlcc_g)
173 0 : IF (ASSOCIATED(rho_nlcc)) THEN
174 0 : NULLIFY (rho_r, rho_g)
175 0 : CALL qs_rho_get(rho0_struct, rho_r=rho_r, rho_g=rho_g)
176 0 : factor = 1.0_dp
177 0 : DO ispin = 1, nspins
178 0 : CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
179 0 : CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
180 : END DO
181 : END IF
182 :
183 : ! Integration weights
184 0 : NULLIFY (weights)
185 0 : CALL get_qs_env(qs_env, xcint_weights=weights)
186 :
187 0 : NULLIFY (weights_uf)
188 0 : IF (uf_grid) THEN
189 0 : IF (ASSOCIATED(weights)) THEN
190 0 : ALLOCATE (weights_uf)
191 0 : CALL xc_pw_pool%create_pw(weights_uf)
192 : BLOCK
193 : TYPE(pw_c1d_gs_type) :: weights_g, weights_g_uf
194 0 : CALL auxbas_pw_pool%create_pw(weights_g)
195 0 : CALL xc_pw_pool%create_pw(weights_g_uf)
196 0 : CALL pw_transfer(weights, weights_g)
197 0 : CALL pw_transfer(weights_g, weights_g_uf)
198 0 : CALL pw_transfer(weights_g_uf, weights_uf)
199 0 : CALL xc_pw_pool%give_back_pw(weights_g_uf)
200 0 : CALL auxbas_pw_pool%give_back_pw(weights_g)
201 : END BLOCK
202 : END IF
203 : !
204 0 : ALLOCATE (rho0_uf, rho1_uf)
205 0 : CALL qs_rho_create(rho0_uf)
206 0 : CALL qs_rho_create(rho1_uf)
207 0 : CALL qs_rho_transfer(rho0_struct, rho0_uf, auxbas_pw_pool, xc_pw_pool)
208 0 : CALL qs_rho_transfer(rho1_struct, rho1_uf, auxbas_pw_pool, xc_pw_pool)
209 0 : uf_pw_pool => xc_pw_pool
210 : ELSE
211 0 : weights_uf => weights
212 0 : rho0_uf => rho0_struct
213 0 : rho1_uf => rho1_struct
214 0 : uf_pw_pool => auxbas_pw_pool
215 : END IF
216 :
217 : ! only spinflip is activated !?
218 0 : IF (.NOT. spinflip) THEN
219 0 : CPABORT("Analytic 3rd EXC derivatives not available")
220 : END IF
221 : CALL qs_fgxc_analytic(rho0_uf, rho1_uf, xc_section, uf_pw_pool, weights_uf, &
222 : is_triplet, spinflip, &
223 0 : fxc_rho_uf, fxc_tau_uf, gxc_rho_uf, gxc_tau_uf)
224 :
225 : ! check for results
226 0 : has_rho = .FALSE.
227 0 : IF (ASSOCIATED(fxc_rho_uf) .AND. ASSOCIATED(gxc_rho_uf)) has_rho = .TRUE.
228 0 : has_tau = .FALSE.
229 0 : IF (ASSOCIATED(fxc_tau_uf) .AND. ASSOCIATED(gxc_tau_uf)) has_tau = .TRUE.
230 : ! transfer to standard grid
231 0 : IF (uf_grid) THEN
232 : ! rho
233 0 : IF (has_rho) THEN
234 0 : ALLOCATE (fxc_rho(nspins), gxc_rho(nspins))
235 0 : DO ispin = 1, nspins
236 0 : CALL auxbas_pw_pool%create_pw(fxc_rho(ispin))
237 0 : CALL auxbas_pw_pool%create_pw(gxc_rho(ispin))
238 : BLOCK
239 : TYPE(pw_c1d_gs_type) :: xc_g, xc_g_uf
240 0 : CALL auxbas_pw_pool%create_pw(xc_g)
241 0 : CALL xc_pw_pool%create_pw(xc_g_uf)
242 0 : CALL pw_transfer(fxc_rho_uf(ispin), xc_g_uf)
243 0 : CALL pw_transfer(xc_g_uf, xc_g)
244 0 : CALL pw_transfer(xc_g, fxc_rho(ispin))
245 0 : CALL pw_transfer(gxc_rho_uf(ispin), xc_g_uf)
246 0 : CALL pw_transfer(xc_g_uf, xc_g)
247 0 : CALL pw_transfer(xc_g, gxc_rho(ispin))
248 0 : CALL xc_pw_pool%give_back_pw(xc_g_uf)
249 0 : CALL auxbas_pw_pool%give_back_pw(xc_g)
250 : END BLOCK
251 0 : CALL xc_pw_pool%give_back_pw(fxc_rho_uf(ispin))
252 0 : CALL xc_pw_pool%give_back_pw(gxc_rho_uf(ispin))
253 : END DO
254 0 : DEALLOCATE (fxc_rho_uf, gxc_rho_uf)
255 : ELSE
256 0 : NULLIFY (fxc_rho, gxc_rho)
257 : END IF
258 : ! tau
259 0 : IF (has_tau) THEN
260 0 : ALLOCATE (fxc_tau(nspins), gxc_tau(nspins))
261 0 : DO ispin = 1, nspins
262 0 : CALL auxbas_pw_pool%create_pw(fxc_tau(ispin))
263 0 : CALL auxbas_pw_pool%create_pw(gxc_tau(ispin))
264 : BLOCK
265 : TYPE(pw_c1d_gs_type) :: xc_g, xc_g_uf
266 0 : CALL auxbas_pw_pool%create_pw(xc_g)
267 0 : CALL xc_pw_pool%create_pw(xc_g_uf)
268 0 : CALL pw_transfer(fxc_tau_uf(ispin), xc_g_uf)
269 0 : CALL pw_transfer(xc_g_uf, xc_g)
270 0 : CALL pw_transfer(xc_g, fxc_tau(ispin))
271 0 : CALL pw_transfer(gxc_tau_uf(ispin), xc_g_uf)
272 0 : CALL pw_transfer(xc_g_uf, xc_g)
273 0 : CALL pw_transfer(xc_g, gxc_tau(ispin))
274 0 : CALL xc_pw_pool%give_back_pw(xc_g_uf)
275 0 : CALL auxbas_pw_pool%give_back_pw(xc_g)
276 : END BLOCK
277 0 : CALL xc_pw_pool%give_back_pw(fxc_tau_uf(ispin))
278 0 : CALL xc_pw_pool%give_back_pw(gxc_tau_uf(ispin))
279 : END DO
280 0 : DEALLOCATE (fxc_tau_uf, gxc_tau_uf)
281 : ELSE
282 0 : NULLIFY (fxc_tau, gxc_tau)
283 : END IF
284 : ELSE
285 0 : fxc_rho => fxc_rho_uf
286 0 : fxc_tau => fxc_tau_uf
287 0 : gxc_rho => gxc_rho_uf
288 0 : gxc_tau => gxc_tau_uf
289 : END IF
290 :
291 : ! de-apply NLCC density
292 0 : IF (ASSOCIATED(rho_nlcc)) THEN
293 0 : factor = -1.0_dp
294 0 : DO ispin = 1, nspins
295 0 : CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
296 0 : CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
297 : END DO
298 : END IF
299 :
300 334 : ELSE IF (analytic_2nd_deriv .AND. .NOT. analytic_3rd_deriv) THEN
301 :
302 334 : NULLIFY (weights)
303 334 : CALL get_qs_env(qs_env, xcint_weights=weights)
304 : CALL qs_fgxc_gdiff(qs_env, rho0_struct, rho1_struct, xc_section, &
305 : order, eps_delta, is_triplet, spinflip, &
306 334 : fxc_rho, fxc_tau, gxc_rho, gxc_tau)
307 :
308 0 : ELSE IF (.NOT. (analytic_2nd_deriv .OR. analytic_3rd_deriv)) THEN
309 :
310 : CALL qs_fgxc_fgdiff(qs_env, rho0_struct, rho1_struct, xc_section, &
311 : order, is_triplet, spinflip, &
312 0 : fxc_rho, fxc_tau, gxc_rho, gxc_tau)
313 :
314 : ELSE
315 0 : CPABORT("EXC derivative option not available")
316 : END IF
317 :
318 334 : IF (do_onecenter) THEN
319 114 : IF (analytic_2nd_deriv .AND. analytic_3rd_deriv) THEN
320 0 : CPABORT("Analytic 3rd EXC derivatives not available")
321 114 : ELSE IF (analytic_2nd_deriv) THEN
322 : CALL fgxc_atom_diff(qs_env, rho0_atom_set, rhof_atom_set, rhog_atom_set, &
323 114 : qs_kind_set, xc_section, is_triplet, order, eps_delta)
324 0 : ELSE IF (analytic_3rd_deriv) THEN
325 : CALL cp_warn(__LOCATION__, "XC Analytic 3rd derivatives cannot be combined with "// &
326 0 : "numeric 2nd derivatives.")
327 0 : CPABORT("XC derivative options")
328 : ELSE
329 : CALL fgxc_atom_calc(qs_env, rho0_atom_set, rhof_atom_set, rhog_atom_set, &
330 0 : qs_kind_set, xc_section, is_triplet, order)
331 : END IF
332 : END IF
333 :
334 334 : CALL timestop(handle)
335 :
336 334 : END SUBROUTINE qs_fgxc_create
337 :
338 : ! **************************************************************************************************
339 : !> \brief Calculates the values at the grid points in real space (r_i), of the second and third
340 : !> functional derivatives of the exchange-correlation energy functional.
341 : !> fxc_rho(r_i) = fxc[n](r_i)*n^(1)(r_i) ! Second functional derivative
342 : !> gxc_rho(r_i) = n^(1)(r_i)*gxc[n](r_i)*n^(1)(r_i) ! Third functional derivative
343 : !> \param rho0_struct Ground state density, n(r).
344 : !> \param rho1_struct Density used to fold the functional derivatives, n^(1)(r).
345 : !> \param xc_section ...
346 : !> \param pw_pool ...
347 : !> \param weights ...
348 : !> \param is_triplet Flag used to activate triplet kernel (not available)
349 : !> \param spinflip Flag used to activate the spin-flip noncollinear kernel and kernel derivatives.
350 : !> \param fxc_rho Second functional derivative with respect to the density, n(r).
351 : !> \param fxc_tau mGGA contribution to the second functional derivative with respect to the density.
352 : !> \param gxc_rho Third functional derivative with respect to the density, n(r).
353 : !> \param gxc_tau mGGA contribution to the third functional derivative with respect to the density.
354 : !> \par History
355 : !> * 07.2024 Created [LHS]
356 : ! **************************************************************************************************
357 0 : SUBROUTINE qs_fgxc_analytic(rho0_struct, rho1_struct, xc_section, pw_pool, weights, &
358 : is_triplet, spinflip, fxc_rho, fxc_tau, gxc_rho, gxc_tau)
359 :
360 : TYPE(qs_rho_type), POINTER :: rho0_struct, rho1_struct
361 : TYPE(section_vals_type), POINTER :: xc_section
362 : TYPE(pw_pool_type), POINTER :: pw_pool
363 : TYPE(pw_r3d_rs_type), POINTER :: weights
364 : LOGICAL, INTENT(IN) :: is_triplet, spinflip
365 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau, gxc_rho, gxc_tau
366 :
367 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_fgxc_analytic'
368 :
369 : INTEGER :: handle, ispin, nspins, spindim
370 : INTEGER, DIMENSION(2, 3) :: bo
371 : LOGICAL :: do_sf, lsd
372 : REAL(KIND=dp) :: fac
373 0 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho0_g, rho1_g
374 0 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho0_r, rho1_r, tau0_r, tau1_r
375 : TYPE(section_vals_type), POINTER :: xc_fun_section
376 : TYPE(xc_derivative_set_type) :: deriv_set
377 : TYPE(xc_rho_cflags_type) :: needs
378 : TYPE(xc_rho_set_type) :: rho0_set, rho1_set
379 :
380 0 : CALL timeset(routineN, handle)
381 :
382 : ! Only rho0 and rho1 should be associated
383 0 : CPASSERT(.NOT. ASSOCIATED(fxc_rho))
384 0 : CPASSERT(.NOT. ASSOCIATED(fxc_tau))
385 0 : CPASSERT(.NOT. ASSOCIATED(gxc_rho))
386 0 : CPASSERT(.NOT. ASSOCIATED(gxc_tau))
387 0 : CPASSERT(ASSOCIATED(rho0_struct))
388 0 : CPASSERT(ASSOCIATED(rho1_struct))
389 :
390 : ! only tested for spinflip
391 0 : CPASSERT(.NOT. is_triplet)
392 0 : CPASSERT(spinflip)
393 :
394 : ! Initialize parameters
395 0 : do_sf = spinflip
396 : !
397 : ! Get the values on the gridpoints of the rho0 density
398 0 : CALL qs_rho_get(rho0_struct, rho_r=rho0_r, rho_g=rho0_g, tau_r=tau0_r)
399 0 : nspins = SIZE(rho0_r)
400 0 : lsd = (nspins == 2)
401 : !
402 0 : IF (do_sf) THEN
403 : spindim = 1
404 : ELSE
405 0 : spindim = nspins
406 : END IF
407 : !
408 0 : fac = 0._dp
409 0 : IF (nspins == 1) THEN
410 0 : fac = 1.0_dp
411 : END IF
412 :
413 : ! Read xc functional section and find out what the functional actually needs
414 0 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
415 0 : needs = xc_functionals_get_needs(xc_fun_section, lsd, .TRUE.)
416 :
417 : ! Create fields for the kernel and kernel derivative
418 0 : ALLOCATE (fxc_rho(spindim), gxc_rho(nspins))
419 0 : DO ispin = 1, spindim
420 0 : CALL pw_pool%create_pw(fxc_rho(ispin))
421 0 : CALL pw_zero(fxc_rho(ispin))
422 : END DO
423 0 : DO ispin = 1, nspins
424 0 : CALL pw_pool%create_pw(gxc_rho(ispin))
425 0 : CALL pw_zero(gxc_rho(ispin))
426 : END DO
427 : ! Create fields for mGGA functionals. This implementation is not ready yet!
428 0 : IF (needs%tau .OR. needs%tau_spin) THEN
429 0 : IF (.NOT. ASSOCIATED(tau1_r)) THEN
430 0 : CPABORT("Tau-dependent functionals requires allocated kinetic energy density grid")
431 : END IF
432 0 : ALLOCATE (fxc_tau(spindim), gxc_tau(nspins))
433 0 : DO ispin = 1, spindim
434 0 : CALL pw_pool%create_pw(fxc_tau(ispin))
435 0 : CALL pw_zero(fxc_tau(ispin))
436 : END DO
437 0 : DO ispin = 1, nspins
438 0 : CALL pw_pool%create_pw(gxc_tau(ispin))
439 0 : CALL pw_zero(gxc_tau(ispin))
440 : END DO
441 : END IF
442 :
443 : ! Build rho0_set
444 : ! calculate the arguments needed by the functionals
445 : ! Needs
446 : ! deriv_set xc_derivative_set_type just declared
447 : ! rho0_set xc_rho_set_type just declared
448 : ! rho0_r pw_type calculated by qs_rho_get
449 : ! pw_pool given by the calling subroutine
450 : ! xc_section given by the calling subroutine
451 : ! tau0_r pw_type calculated by qs_rho_get
452 : CALL xc_prep_3rd_deriv(deriv_set, rho0_set, rho0_r, pw_pool, weights, &
453 0 : xc_section, tau_r=tau0_r, do_sf=do_sf)
454 :
455 : ! Build rho1_set
456 : ! Get the values on the gridpoints of the rho1 density
457 0 : CALL qs_rho_get(rho1_struct, rho_r=rho1_r, rho_g=rho1_g, tau_r=tau1_r)
458 0 : bo = rho1_r(1)%pw_grid%bounds_local
459 : ! create the place where to store the argument for the functionals
460 : ! Needs
461 : ! rho1_set xc_rho_set_type just declared
462 : ! bo 2x3 integer matrix should have bounds_local or rho1_r
463 : CALL xc_rho_set_create(rho1_set, bo, &
464 : rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
465 : drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
466 0 : tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
467 : ! calculate the arguments needed by the functionals
468 : ! Needs
469 : ! rho1_set object created by xc_rho_set_create
470 : ! rho1_r,rho1_g,tau1_r pw_type values of rho1 in real space grid
471 : ! needs xc_rho_cflags_type defined through xc_functionals_get_needs
472 : ! pw_pool Given by the calling subroutine
473 : CALL xc_rho_set_update(rho1_set, rho1_r, rho1_g, tau1_r, needs, &
474 : section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
475 : section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
476 0 : pw_pool, spinflip=do_sf)
477 :
478 : ! Calculate exchange correlation kernel
479 : ! Needs
480 : ! fxc_rho, fxc_tau pw_type not associated
481 : ! deriv_set created and defined by xc_prep_3rd_deriv
482 : ! rho0_set xc_rho_set_type build by xc_prep_3rd_deriv
483 : ! rho1_set xc_rho_set_type build by xc_rho_set_create/update
484 : ! pw_pool needs to be given by the calling subroutine
485 : ! xc_section needs to be given by the calling subroutine
486 : CALL xc_calc_2nd_deriv_analytical(fxc_rho, fxc_tau, deriv_set, rho0_set, rho1_set, pw_pool, &
487 0 : xc_section, .FALSE., spinflip=do_sf, tddfpt_fac=fac)
488 : ! Calculate exchange correlation kernel derivative
489 : CALL xc_calc_3rd_deriv_analytical(gxc_rho, gxc_tau, deriv_set, rho0_set, rho1_set, pw_pool, &
490 0 : xc_section, spinflip=do_sf)
491 :
492 0 : CALL xc_dset_release(deriv_set)
493 0 : CALL xc_rho_set_release(rho0_set)
494 0 : CALL xc_rho_set_release(rho1_set)
495 :
496 0 : CALL timestop(handle)
497 :
498 0 : END SUBROUTINE qs_fgxc_analytic
499 :
500 : ! **************************************************************************************************
501 : !> \brief ...
502 : !> \param qs_env ...
503 : !> \param rho0_struct ...
504 : !> \param rho1_struct ...
505 : !> \param xc_section ...
506 : !> \param accuracy ...
507 : !> \param epsrho ...
508 : !> \param is_triplet ...
509 : !> \param spinflip ...
510 : !> \param fxc_rho ...
511 : !> \param fxc_tau ...
512 : !> \param gxc_rho ...
513 : !> \param gxc_tau ...
514 : ! **************************************************************************************************
515 334 : SUBROUTINE qs_fgxc_gdiff(qs_env, rho0_struct, rho1_struct, xc_section, accuracy, epsrho, &
516 : is_triplet, spinflip, fxc_rho, fxc_tau, gxc_rho, gxc_tau)
517 :
518 : TYPE(qs_environment_type), POINTER :: qs_env
519 : TYPE(qs_rho_type), POINTER :: rho0_struct, rho1_struct
520 : TYPE(section_vals_type), POINTER :: xc_section
521 : INTEGER, INTENT(IN) :: accuracy
522 : REAL(KIND=dp), INTENT(IN) :: epsrho
523 : LOGICAL, INTENT(IN) :: is_triplet, spinflip
524 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau, gxc_rho, gxc_tau
525 :
526 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_fgxc_gdiff'
527 :
528 : INTEGER :: handle, ispin, istep, nspins, nstep
529 : LOGICAL :: do_sf
530 : REAL(KIND=dp) :: alpha, beta, exc, oeps1
531 : REAL(KIND=dp), DIMENSION(-4:4) :: ak
532 : TYPE(dft_control_type), POINTER :: dft_control
533 : TYPE(pw_env_type), POINTER :: pw_env
534 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
535 334 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_tau_rspace, vxc00, vxc00b
536 : TYPE(qs_rho_type), POINTER :: rhoin
537 334 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho0_atom_set, rho1_atom_set
538 :
539 334 : CALL timeset(routineN, handle)
540 :
541 334 : CPASSERT(.NOT. ASSOCIATED(fxc_rho))
542 334 : CPASSERT(.NOT. ASSOCIATED(fxc_tau))
543 334 : CPASSERT(.NOT. ASSOCIATED(gxc_rho))
544 334 : CPASSERT(.NOT. ASSOCIATED(gxc_tau))
545 334 : CPASSERT(ASSOCIATED(rho0_struct))
546 334 : CPASSERT(ASSOCIATED(rho1_struct))
547 :
548 334 : do_sf = spinflip
549 :
550 334 : ak = 0.0_dp
551 334 : SELECT CASE (accuracy)
552 : CASE (:4)
553 0 : nstep = 2
554 0 : ak(-2:2) = [1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp]/12.0_dp
555 : CASE (5:7)
556 2672 : nstep = 3
557 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
558 : CASE (8:)
559 0 : nstep = 4
560 : ak(-4:4) = [1.0_dp, -32.0_dp/3.0_dp, 56.0_dp, -224.0_dp, 0.0_dp, &
561 334 : 224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp]/280.0_dp
562 : END SELECT
563 :
564 334 : CALL get_qs_env(qs_env, dft_control=dft_control, pw_env=pw_env)
565 334 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
566 :
567 334 : nspins = dft_control%nspins
568 334 : exc = 0.0_dp
569 :
570 334 : IF (do_sf) THEN
571 : CALL qs_fxc_create(qs_env, rho0_struct, rho1_struct, rho0_atom_set, xc_section, &
572 : .FALSE., fxc_rho, fxc_tau, rho1_atom_set, &
573 8 : is_triplet=is_triplet, spinflip=do_sf)
574 : ELSE
575 : CALL qs_fxc_fdiff(qs_env, rho0_struct, rho1_struct, xc_section, &
576 326 : accuracy, fxc_rho, fxc_tau, is_triplet)
577 : END IF
578 :
579 2672 : DO istep = -nstep, nstep
580 :
581 2672 : IF (ak(istep) /= 0.0_dp) THEN
582 2004 : alpha = 1.0_dp
583 2004 : beta = REAL(istep, KIND=dp)*epsrho
584 : NULLIFY (rhoin)
585 2004 : ALLOCATE (rhoin)
586 2004 : CALL qs_rho_create(rhoin)
587 2004 : NULLIFY (vxc00, vxc00b, v_tau_rspace)
588 2004 : CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
589 2004 : CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, beta)
590 2004 : IF (do_sf) THEN
591 : ! variation in alpha density
592 : CALL qs_fxc_create(qs_env, rhoin, rho1_struct, rho0_atom_set, xc_section, &
593 : .FALSE., vxc00, v_tau_rspace, rho1_atom_set, &
594 48 : is_triplet=is_triplet, spinflip=do_sf)
595 : ! variation in beta density
596 48 : CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
597 48 : CALL qs_rho_scale_and_add_b(rhoin, rho1_struct, alpha, beta)
598 : CALL qs_fxc_create(qs_env, rhoin, rho1_struct, rho0_atom_set, xc_section, &
599 : .FALSE., vxc00b, v_tau_rspace, rho1_atom_set, &
600 48 : is_triplet=is_triplet, spinflip=do_sf)
601 : ELSE
602 : CALL qs_fxc_fdiff(qs_env, rhoin, rho1_struct, xc_section, accuracy, &
603 1956 : fxc_rho=vxc00, fxc_tau=v_tau_rspace, is_triplet=is_triplet)
604 : END IF
605 2004 : CALL qs_rho_release(rhoin)
606 2004 : DEALLOCATE (rhoin)
607 2004 : IF (.NOT. ASSOCIATED(gxc_rho)) THEN
608 1380 : ALLOCATE (gxc_rho(nspins))
609 712 : DO ispin = 1, nspins
610 378 : CALL auxbas_pw_pool%create_pw(gxc_rho(ispin))
611 712 : CALL pw_zero(gxc_rho(ispin))
612 : END DO
613 : END IF
614 2004 : IF (do_sf) THEN
615 48 : CALL pw_axpy(vxc00(1), gxc_rho(1), ak(istep))
616 48 : CALL pw_axpy(vxc00b(1), gxc_rho(2), ak(istep))
617 : ELSE
618 4128 : DO ispin = 1, nspins
619 4128 : CALL pw_axpy(vxc00(ispin), gxc_rho(ispin), ak(istep))
620 : END DO
621 : END IF
622 4224 : DO ispin = 1, SIZE(vxc00)
623 4224 : CALL auxbas_pw_pool%give_back_pw(vxc00(ispin))
624 : END DO
625 2004 : DEALLOCATE (vxc00)
626 2004 : IF (ASSOCIATED(vxc00b)) THEN
627 48 : CALL auxbas_pw_pool%give_back_pw(vxc00b(1))
628 48 : DEALLOCATE (vxc00b)
629 : END IF
630 2004 : IF (ASSOCIATED(v_tau_rspace)) THEN
631 0 : IF (.NOT. ASSOCIATED(gxc_tau)) THEN
632 0 : ALLOCATE (gxc_tau(nspins))
633 0 : DO ispin = 1, nspins
634 0 : CALL auxbas_pw_pool%create_pw(gxc_tau(ispin))
635 0 : CALL pw_zero(gxc_tau(ispin))
636 : END DO
637 : END IF
638 0 : DO ispin = 1, nspins
639 0 : CALL pw_axpy(v_tau_rspace(ispin), gxc_tau(ispin), ak(istep))
640 : END DO
641 0 : DO ispin = 1, SIZE(v_tau_rspace)
642 0 : CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
643 : END DO
644 0 : DEALLOCATE (v_tau_rspace)
645 : END IF
646 : END IF
647 :
648 : END DO
649 :
650 334 : oeps1 = 1.0_dp/epsrho
651 712 : DO ispin = 1, nspins
652 712 : CALL pw_scale(gxc_rho(ispin), oeps1)
653 : END DO
654 334 : IF (ASSOCIATED(gxc_tau)) THEN
655 0 : DO ispin = 1, nspins
656 0 : CALL pw_scale(gxc_tau(ispin), oeps1)
657 : END DO
658 : END IF
659 :
660 334 : CALL timestop(handle)
661 :
662 334 : END SUBROUTINE qs_fgxc_gdiff
663 :
664 : ! **************************************************************************************************
665 : !> \brief ...
666 : !> \param qs_env ...
667 : !> \param rho0_struct ...
668 : !> \param rho1_struct ...
669 : !> \param xc_section ...
670 : !> \param accuracy ...
671 : !> \param is_triplet ...
672 : !> \param spinflip ...
673 : !> \param fxc_rho ...
674 : !> \param fxc_tau ...
675 : !> \param gxc_rho ...
676 : !> \param gxc_tau ...
677 : ! **************************************************************************************************
678 0 : SUBROUTINE qs_fgxc_fgdiff(qs_env, rho0_struct, rho1_struct, xc_section, &
679 : accuracy, is_triplet, spinflip, &
680 : fxc_rho, fxc_tau, gxc_rho, gxc_tau)
681 :
682 : TYPE(qs_environment_type), POINTER :: qs_env
683 : TYPE(qs_rho_type), POINTER :: rho0_struct, rho1_struct
684 : TYPE(section_vals_type), POINTER :: xc_section
685 : INTEGER, INTENT(IN) :: accuracy
686 : LOGICAL, INTENT(IN) :: is_triplet, spinflip
687 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau, gxc_rho, gxc_tau
688 :
689 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_fgxc_fgdiff'
690 : REAL(KIND=dp), PARAMETER :: epsrho = 5.e-4_dp
691 :
692 : INTEGER :: handle, ispin, istep, nspins, nstep
693 : REAL(KIND=dp) :: alpha, beta, exc, oeps1, oeps2
694 : REAL(KIND=dp), DIMENSION(-4:4) :: ak, bl
695 : TYPE(dft_control_type), POINTER :: dft_control
696 : TYPE(pw_env_type), POINTER :: pw_env
697 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
698 0 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_tau_rspace, vxc00
699 : TYPE(qs_ks_env_type), POINTER :: ks_env
700 : TYPE(qs_rho_type), POINTER :: rhoin
701 :
702 0 : CALL timeset(routineN, handle)
703 :
704 0 : CPASSERT(.NOT. ASSOCIATED(fxc_rho))
705 0 : CPASSERT(.NOT. ASSOCIATED(fxc_tau))
706 0 : CPASSERT(.NOT. ASSOCIATED(gxc_rho))
707 0 : CPASSERT(.NOT. ASSOCIATED(gxc_tau))
708 0 : CPASSERT(ASSOCIATED(rho0_struct))
709 0 : CPASSERT(ASSOCIATED(rho1_struct))
710 :
711 : ! no spinflip tested
712 0 : CPASSERT(.NOT. spinflip)
713 :
714 0 : ak = 0.0_dp
715 0 : bl = 0.0_dp
716 0 : SELECT CASE (accuracy)
717 : CASE (:4)
718 0 : nstep = 2
719 0 : ak(-2:2) = [1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp]/12.0_dp
720 0 : bl(-2:2) = [-1.0_dp, 16.0_dp, -30.0_dp, 16.0_dp, -1.0_dp]/12.0_dp
721 : CASE (5:7)
722 0 : nstep = 3
723 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
724 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
725 : CASE (8:)
726 0 : nstep = 4
727 : ak(-4:4) = [1.0_dp, -32.0_dp/3.0_dp, 56.0_dp, -224.0_dp, 0.0_dp, &
728 0 : 224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp]/280.0_dp
729 : bl(-4:4) = [-1.0_dp, 128.0_dp/9.0_dp, -112.0_dp, 896.0_dp, -14350.0_dp/9.0_dp, &
730 0 : 896.0_dp, -112.0_dp, 128.0_dp/9.0_dp, -1.0_dp]/560.0_dp
731 : END SELECT
732 :
733 0 : CALL get_qs_env(qs_env, ks_env=ks_env, dft_control=dft_control, pw_env=pw_env)
734 0 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
735 :
736 0 : nspins = dft_control%nspins
737 : exc = 0.0_dp
738 :
739 0 : DO istep = -nstep, nstep
740 :
741 0 : alpha = 1.0_dp
742 0 : beta = REAL(istep, KIND=dp)*epsrho
743 : NULLIFY (rhoin)
744 0 : ALLOCATE (rhoin)
745 0 : CALL qs_rho_create(rhoin)
746 0 : NULLIFY (vxc00, v_tau_rspace)
747 0 : IF (is_triplet) THEN
748 0 : CPASSERT(nspins == 1)
749 : ! rhoin = (0.5 rho0, 0.5 rho0)
750 0 : CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, 2)
751 : ! rhoin = (0.5 rho0 + 0.5 rho1, 0.5 rho0)
752 0 : CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, 0.5_dp*beta)
753 0 : CALL qs_vxc_create(ks_env, rhoin, xc_section, vxc00, v_tau_rspace, exc)
754 0 : CALL pw_axpy(vxc00(2), vxc00(1), -1.0_dp)
755 : ELSE
756 0 : CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
757 0 : CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, beta)
758 0 : CALL qs_vxc_create(ks_env, rhoin, xc_section, vxc00, v_tau_rspace, exc)
759 : END IF
760 0 : CALL qs_rho_release(rhoin)
761 0 : DEALLOCATE (rhoin)
762 0 : IF (.NOT. ASSOCIATED(fxc_rho)) THEN
763 0 : ALLOCATE (fxc_rho(nspins))
764 0 : DO ispin = 1, nspins
765 0 : CALL auxbas_pw_pool%create_pw(fxc_rho(ispin))
766 0 : CALL pw_zero(fxc_rho(ispin))
767 : END DO
768 : END IF
769 0 : IF (.NOT. ASSOCIATED(gxc_rho)) THEN
770 0 : ALLOCATE (gxc_rho(nspins))
771 0 : DO ispin = 1, nspins
772 0 : CALL auxbas_pw_pool%create_pw(gxc_rho(ispin))
773 0 : CALL pw_zero(gxc_rho(ispin))
774 : END DO
775 : END IF
776 0 : CPASSERT(.NOT. ASSOCIATED(v_tau_rspace))
777 0 : DO ispin = 1, nspins
778 0 : IF (ak(istep) /= 0.0_dp) THEN
779 0 : CALL pw_axpy(vxc00(ispin), fxc_rho(ispin), ak(istep))
780 : END IF
781 0 : IF (bl(istep) /= 0.0_dp) THEN
782 0 : CALL pw_axpy(vxc00(ispin), gxc_rho(ispin), bl(istep))
783 : END IF
784 : END DO
785 0 : DO ispin = 1, SIZE(vxc00)
786 0 : CALL auxbas_pw_pool%give_back_pw(vxc00(ispin))
787 : END DO
788 0 : DEALLOCATE (vxc00)
789 :
790 : END DO
791 :
792 0 : oeps1 = 1.0_dp/epsrho
793 0 : oeps2 = 1.0_dp/(epsrho**2)
794 0 : DO ispin = 1, nspins
795 0 : CALL pw_scale(fxc_rho(ispin), oeps1)
796 0 : CALL pw_scale(gxc_rho(ispin), oeps2)
797 : END DO
798 :
799 0 : CALL timestop(handle)
800 :
801 0 : END SUBROUTINE qs_fgxc_fgdiff
802 :
803 : ! **************************************************************************************************
804 : !> \brief ...
805 : !> \param fxc_rho ...
806 : !> \param fxc_tau ...
807 : !> \param gxc_rho ...
808 : !> \param gxc_tau ...
809 : !> \param pw_pool ...
810 : ! **************************************************************************************************
811 334 : SUBROUTINE qs_fgxc_release(fxc_rho, fxc_tau, gxc_rho, gxc_tau, pw_pool)
812 :
813 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau, gxc_rho, gxc_tau
814 : TYPE(pw_pool_type), POINTER :: pw_pool
815 :
816 : INTEGER :: ispin
817 :
818 334 : IF (ASSOCIATED(fxc_rho)) THEN
819 704 : DO ispin = 1, SIZE(fxc_rho)
820 704 : CALL pw_pool%give_back_pw(fxc_rho(ispin))
821 : END DO
822 334 : DEALLOCATE (fxc_rho)
823 : END IF
824 334 : IF (ASSOCIATED(fxc_tau)) THEN
825 0 : DO ispin = 1, SIZE(fxc_tau)
826 0 : CALL pw_pool%give_back_pw(fxc_tau(ispin))
827 : END DO
828 0 : DEALLOCATE (fxc_tau)
829 : END IF
830 334 : IF (ASSOCIATED(gxc_rho)) THEN
831 712 : DO ispin = 1, SIZE(gxc_rho)
832 712 : CALL pw_pool%give_back_pw(gxc_rho(ispin))
833 : END DO
834 334 : DEALLOCATE (gxc_rho)
835 : END IF
836 334 : IF (ASSOCIATED(gxc_tau)) THEN
837 0 : DO ispin = 1, SIZE(gxc_tau)
838 0 : CALL pw_pool%give_back_pw(gxc_tau(ispin))
839 : END DO
840 0 : DEALLOCATE (gxc_tau)
841 : END IF
842 :
843 334 : END SUBROUTINE qs_fgxc_release
844 :
845 : END MODULE qs_fgxc
|
