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