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