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 Setup Routine for Fxc Potentials
10 : ! https://en.wikipedia.org/wiki/Finite_difference_coefficient
11 : !---------------------------------------------------------------------------------------------------
12 : !Derivative Accuracy 4 3 2 1 0 1 2 3 4
13 : !---------------------------------------------------------------------------------------------------
14 : ! 1 2 -1/2 0 1/2
15 : ! 4 1/12 -2/3 0 2/3 -1/12
16 : ! 6 -1/60 3/20 -3/4 0 3/4 -3/20 1/60
17 : ! 8 1/280 -4/105 1/5 -4/5 0 4/5 -1/5 4/105 -1/280
18 : !---------------------------------------------------------------------------------------------------
19 : ! 2 2 1 -2 1
20 : ! 4 -1/12 4/3 -5/2 4/3 -1/12
21 : ! 6 1/90 -3/20 3/2 -49/18 3/2 -3/20 1/90
22 : ! 8 -1/560 8/315 -1/5 8/5 -205/72 8/5 -1/5 8/315 -1/560
23 : !---------------------------------------------------------------------------------------------------
24 : !> \par History
25 : !> init 17.03.2020
26 : !> complete refactoring 08.2026
27 : !> \author JGH
28 : ! **************************************************************************************************
29 : MODULE qs_fxc
30 :
31 : USE cp_control_types, ONLY: dft_control_type
32 : USE input_section_types, ONLY: section_get_ival,&
33 : section_get_lval,&
34 : section_get_rval,&
35 : section_vals_get_subs_vals,&
36 : section_vals_type
37 : USE kinds, ONLY: dp
38 : USE message_passing, ONLY: mp_para_env_type
39 : USE pw_env_types, ONLY: pw_env_get,&
40 : pw_env_type
41 : USE pw_grids, ONLY: pw_grid_compare
42 : USE pw_methods, ONLY: pw_axpy,&
43 : pw_scale,&
44 : pw_transfer,&
45 : pw_zero
46 : USE pw_pool_types, ONLY: pw_pool_type
47 : USE pw_types, ONLY: pw_c1d_gs_type,&
48 : pw_r3d_rs_type
49 : USE qs_environment_types, ONLY: get_qs_env,&
50 : qs_environment_type
51 : USE qs_fxc_atom, ONLY: fxc_atom_calc
52 : USE qs_kind_types, ONLY: qs_kind_type
53 : USE qs_ks_types, ONLY: qs_ks_env_type
54 : USE qs_rho_atom_types, ONLY: rho_atom_type
55 : USE qs_rho_methods, ONLY: qs_rho_copy,&
56 : qs_rho_scale_and_add,&
57 : qs_rho_transfer
58 : USE qs_rho_types, ONLY: qs_rho_create,&
59 : qs_rho_get,&
60 : qs_rho_release,&
61 : qs_rho_type
62 : USE qs_vxc, ONLY: qs_vxc_create
63 : USE xc, ONLY: xc_calc_2nd_deriv_analytical,&
64 : xc_calc_2nd_deriv_numerical,&
65 : xc_prep_2nd_deriv
66 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
67 : xc_dset_release
68 : USE xc_derivatives, ONLY: xc_functionals_get_needs
69 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
70 : USE xc_rho_set_types, ONLY: xc_rho_set_create,&
71 : xc_rho_set_release,&
72 : xc_rho_set_type,&
73 : xc_rho_set_update
74 : #include "./base/base_uses.f90"
75 :
76 : IMPLICIT NONE
77 :
78 : PRIVATE
79 :
80 : ! *** Public subroutines ***
81 : PUBLIC :: qs_fxc_create, qs_fxc_prep, qs_fxc_apply
82 : PUBLIC :: qs_fxc_fdiff
83 :
84 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_fxc'
85 :
86 : ! **************************************************************************************************
87 :
88 : CONTAINS
89 :
90 : ! **************************************************************************************************
91 : !> \brief ...
92 : !> \param qs_env ...
93 : !> \param rho0_struct ...
94 : !> \param rho1_struct ...
95 : !> \param rho0_atom_set ...
96 : !> \param xc_section ...
97 : !> \param do_onecenter ...
98 : !> \param fxc_rho ...
99 : !> \param fxc_tau ...
100 : !> \param rho1_atom_set ...
101 : !> \param do_scale ...
102 : !> \param is_triplet ...
103 : !> \param spinflip ...
104 : !> \param no_weights ...
105 : !> \param uf_grid_results ...
106 : !> \param pw_env_ext ...
107 : !> \param kind_set_external ...
108 : !> \param para_env_external ...
109 : !> \param compute_virial ...
110 : !> \param virial_xc ...
111 : ! **************************************************************************************************
112 3930 : SUBROUTINE qs_fxc_create(qs_env, rho0_struct, rho1_struct, rho0_atom_set, &
113 : xc_section, do_onecenter, &
114 : fxc_rho, fxc_tau, rho1_atom_set, &
115 : do_scale, is_triplet, spinflip, no_weights, uf_grid_results, &
116 : pw_env_ext, kind_set_external, para_env_external, &
117 : compute_virial, virial_xc)
118 :
119 : TYPE(qs_environment_type), POINTER :: qs_env
120 : TYPE(qs_rho_type), POINTER :: rho0_struct, rho1_struct
121 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho0_atom_set
122 : TYPE(section_vals_type), POINTER :: xc_section
123 : LOGICAL, INTENT(IN) :: do_onecenter
124 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau
125 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho1_atom_set
126 : LOGICAL, INTENT(IN), OPTIONAL :: do_scale, is_triplet, spinflip, &
127 : no_weights, uf_grid_results
128 : TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_ext
129 : TYPE(qs_kind_type), DIMENSION(:), OPTIONAL, &
130 : POINTER :: kind_set_external
131 : TYPE(mp_para_env_type), INTENT(IN), OPTIONAL :: para_env_external
132 : LOGICAL, INTENT(IN), OPTIONAL :: compute_virial
133 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), &
134 : OPTIONAL :: virial_xc
135 :
136 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_fxc_create'
137 :
138 : INTEGER :: handle, ispin, nspins
139 : LOGICAL :: do_virial, do_w, ret_uf, uf_grid
140 : REAL(KIND=dp) :: factor
141 : TYPE(dft_control_type), POINTER :: dft_control
142 3930 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
143 : TYPE(pw_c1d_gs_type), POINTER :: rho_nlcc_g
144 : TYPE(pw_env_type), POINTER :: pw_env
145 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool, xc_pw_pool
146 3930 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho_lo, fxc_rho_uf, fxc_tau_lo, &
147 3930 : fxc_tau_uf, rho_r
148 : TYPE(pw_r3d_rs_type), POINTER :: rho_nlcc, weights, weights_uf
149 : TYPE(qs_rho_type), POINTER :: rho0_uf, rho1_uf
150 :
151 3930 : CALL timeset(routineN, handle)
152 :
153 3930 : do_virial = .FALSE.
154 3930 : IF (PRESENT(compute_virial)) do_virial = compute_virial
155 :
156 3930 : CALL get_qs_env(qs_env, dft_control=dft_control)
157 :
158 3930 : IF (PRESENT(pw_env_ext)) THEN
159 0 : pw_env => pw_env_ext
160 : ELSE
161 3930 : CALL get_qs_env(qs_env, pw_env=pw_env)
162 : END IF
163 3930 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, xc_pw_pool=xc_pw_pool)
164 3930 : uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid, xc_pw_pool%pw_grid)
165 :
166 3930 : nspins = dft_control%nspins
167 3930 : IF (ASSOCIATED(fxc_rho)) THEN
168 0 : CPASSERT(nspins == SIZE(fxc_rho))
169 : END IF
170 3930 : IF (ASSOCIATED(fxc_tau)) THEN
171 0 : CPASSERT(nspins == SIZE(fxc_tau))
172 : END IF
173 :
174 3930 : NULLIFY (rho_nlcc, rho_nlcc_g)
175 3930 : CALL get_qs_env(qs_env, rho_nlcc=rho_nlcc, rho_nlcc_g=rho_nlcc_g)
176 3930 : IF (ASSOCIATED(rho_nlcc)) THEN
177 0 : NULLIFY (rho_r, rho_g)
178 0 : CALL qs_rho_get(rho0_struct, rho_r=rho_r, rho_g=rho_g)
179 0 : factor = 1.0_dp
180 0 : DO ispin = 1, nspins
181 0 : CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
182 0 : CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
183 : END DO
184 : END IF
185 :
186 3930 : do_w = .TRUE.
187 3930 : IF (PRESENT(no_weights)) do_w = .NOT. no_weights
188 62 : IF (do_w) THEN
189 3868 : CALL get_qs_env(qs_env, xcint_weights=weights)
190 : ELSE
191 62 : NULLIFY (weights)
192 : END IF
193 :
194 3930 : NULLIFY (fxc_rho_lo, fxc_tau_lo)
195 3930 : IF (uf_grid) THEN
196 24 : IF (PRESENT(uf_grid_results)) THEN
197 8 : ret_uf = uf_grid_results
198 : ELSE
199 : ret_uf = .FALSE.
200 : END IF
201 24 : NULLIFY (weights_uf)
202 24 : IF (ASSOCIATED(weights)) THEN
203 16 : ALLOCATE (weights_uf)
204 16 : CALL xc_pw_pool%create_pw(weights_uf)
205 : BLOCK
206 : TYPE(pw_c1d_gs_type) :: weights_g, weights_g_uf
207 16 : CALL auxbas_pw_pool%create_pw(weights_g)
208 16 : CALL xc_pw_pool%create_pw(weights_g_uf)
209 16 : CALL pw_transfer(weights, weights_g)
210 16 : CALL pw_transfer(weights_g, weights_g_uf)
211 16 : CALL pw_transfer(weights_g_uf, weights_uf)
212 16 : CALL xc_pw_pool%give_back_pw(weights_g_uf)
213 32 : CALL auxbas_pw_pool%give_back_pw(weights_g)
214 : END BLOCK
215 : END IF
216 : !
217 24 : ALLOCATE (rho0_uf, rho1_uf)
218 24 : CALL qs_rho_create(rho0_uf)
219 24 : CALL qs_rho_create(rho1_uf)
220 24 : CALL qs_rho_transfer(rho0_struct, rho0_uf, auxbas_pw_pool, xc_pw_pool)
221 24 : CALL qs_rho_transfer(rho1_struct, rho1_uf, auxbas_pw_pool, xc_pw_pool)
222 : !
223 24 : NULLIFY (fxc_rho_uf, fxc_tau_uf)
224 : CALL qs_fxc_calculate(rho0_uf, rho1_uf, xc_section, weights_uf, xc_pw_pool, &
225 : fxc_rho_uf, fxc_tau_uf, &
226 : is_triplet=is_triplet, spinflip=spinflip, &
227 24 : compute_virial=do_virial, virial_xc=virial_xc)
228 : !
229 24 : CALL qs_rho_release(rho0_uf)
230 24 : CALL qs_rho_release(rho1_uf)
231 24 : DEALLOCATE (rho0_uf, rho1_uf)
232 24 : IF (ASSOCIATED(weights_uf)) THEN
233 16 : CALL xc_pw_pool%give_back_pw(weights_uf)
234 16 : DEALLOCATE (weights_uf)
235 : END IF
236 24 : IF (ret_uf) THEN
237 8 : fxc_rho_lo => fxc_rho_uf
238 8 : fxc_tau_lo => fxc_tau_uf
239 : ELSE
240 16 : IF (ASSOCIATED(fxc_rho_uf)) THEN
241 64 : ALLOCATE (fxc_rho_lo(nspins))
242 32 : DO ispin = 1, nspins
243 16 : CALL auxbas_pw_pool%create_pw(fxc_rho_lo(ispin))
244 : BLOCK
245 : TYPE(pw_c1d_gs_type) :: fxc_g, fxc_g_uf
246 16 : CALL auxbas_pw_pool%create_pw(fxc_g)
247 16 : CALL xc_pw_pool%create_pw(fxc_g_uf)
248 16 : CALL pw_transfer(fxc_rho_uf(ispin), fxc_g_uf)
249 16 : CALL pw_transfer(fxc_g_uf, fxc_g)
250 16 : CALL pw_transfer(fxc_g, fxc_rho_lo(ispin))
251 16 : CALL xc_pw_pool%give_back_pw(fxc_g_uf)
252 32 : CALL auxbas_pw_pool%give_back_pw(fxc_g)
253 : END BLOCK
254 32 : CALL xc_pw_pool%give_back_pw(fxc_rho_uf(ispin))
255 : END DO
256 16 : DEALLOCATE (fxc_rho_uf)
257 : END IF
258 16 : IF (ASSOCIATED(fxc_tau_uf)) THEN
259 0 : ALLOCATE (fxc_tau_lo(nspins))
260 0 : DO ispin = 1, nspins
261 0 : CALL auxbas_pw_pool%create_pw(fxc_tau_lo(ispin))
262 : BLOCK
263 : TYPE(pw_c1d_gs_type) :: fxc_g, fxc_g_uf
264 0 : CALL auxbas_pw_pool%create_pw(fxc_g)
265 0 : CALL xc_pw_pool%create_pw(fxc_g_uf)
266 0 : CALL pw_transfer(fxc_tau_uf(ispin), fxc_g_uf)
267 0 : CALL pw_transfer(fxc_g_uf, fxc_g)
268 0 : CALL pw_transfer(fxc_g, fxc_tau_lo(ispin))
269 0 : CALL xc_pw_pool%give_back_pw(fxc_g_uf)
270 0 : CALL auxbas_pw_pool%give_back_pw(fxc_g)
271 : END BLOCK
272 0 : CALL xc_pw_pool%give_back_pw(fxc_tau_uf(ispin))
273 : END DO
274 : END IF
275 : END IF
276 :
277 : ELSE
278 : CALL qs_fxc_calculate(rho0_struct, rho1_struct, xc_section, weights, auxbas_pw_pool, &
279 : fxc_rho_lo, fxc_tau_lo, &
280 : is_triplet=is_triplet, spinflip=spinflip, &
281 3906 : compute_virial=do_virial, virial_xc=virial_xc)
282 : END IF
283 :
284 : ! de-apply NLCC density
285 3930 : IF (ASSOCIATED(rho_nlcc)) THEN
286 0 : factor = -1.0_dp
287 0 : DO ispin = 1, nspins
288 0 : CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
289 0 : CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
290 : END DO
291 : END IF
292 :
293 : ! return potentials
294 3930 : IF (ASSOCIATED(fxc_rho)) THEN
295 0 : DO ispin = 1, MIN(SIZE(fxc_rho_lo), SIZE(fxc_rho))
296 0 : CALL pw_transfer(fxc_rho_lo(ispin), fxc_rho(ispin))
297 : END DO
298 0 : DO ispin = 1, SIZE(fxc_rho_lo)
299 0 : CALL auxbas_pw_pool%give_back_pw(fxc_rho_lo(ispin))
300 : END DO
301 0 : DEALLOCATE (fxc_rho_lo)
302 : ELSE
303 3930 : fxc_rho => fxc_rho_lo
304 : END IF
305 3930 : IF (ASSOCIATED(fxc_tau)) THEN
306 0 : IF (ASSOCIATED(fxc_tau_lo)) THEN
307 0 : DO ispin = 1, MIN(SIZE(fxc_tau_lo), SIZE(fxc_tau))
308 0 : CALL pw_transfer(fxc_tau_lo(ispin), fxc_tau(ispin))
309 : END DO
310 0 : DO ispin = 1, SIZE(fxc_tau_lo)
311 0 : CALL auxbas_pw_pool%give_back_pw(fxc_tau_lo(ispin))
312 : END DO
313 0 : DEALLOCATE (fxc_tau_lo)
314 : ELSE
315 0 : DO ispin = 1, nspins
316 0 : CALL pw_zero(fxc_tau(ispin))
317 : END DO
318 : END IF
319 : ELSE
320 3930 : fxc_tau => fxc_tau_lo
321 : END IF
322 :
323 3930 : IF (do_onecenter) THEN
324 : CALL fxc_atom_calc(qs_env, rho0_atom_set, rho1_atom_set, xc_section, &
325 : do_scale=do_scale, do_triplet=is_triplet, do_sf=spinflip, &
326 : para_env_ext=para_env_external, &
327 374 : kind_set_external=kind_set_external)
328 : END IF
329 :
330 3930 : CALL timestop(handle)
331 :
332 3930 : END SUBROUTINE qs_fxc_create
333 :
334 : ! **************************************************************************************************
335 : !> \brief ...
336 : !> \param rho0 ...
337 : !> \param rho1 ...
338 : !> \param xc_section ...
339 : !> \param weights ...
340 : !> \param auxbas_pw_pool ...
341 : !> \param fxc_rho ...
342 : !> \param fxc_tau ...
343 : !> \param is_triplet ...
344 : !> \param spinflip ...
345 : !> \param compute_virial ...
346 : !> \param virial_xc ...
347 : ! **************************************************************************************************
348 3930 : SUBROUTINE qs_fxc_calculate(rho0, rho1, xc_section, weights, auxbas_pw_pool, &
349 : fxc_rho, fxc_tau, is_triplet, spinflip, &
350 : compute_virial, virial_xc)
351 :
352 : TYPE(qs_rho_type), POINTER :: rho0, rho1
353 : TYPE(section_vals_type), POINTER :: xc_section
354 : TYPE(pw_r3d_rs_type), POINTER :: weights
355 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
356 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau
357 : LOGICAL, INTENT(IN), OPTIONAL :: is_triplet, spinflip, compute_virial
358 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), &
359 : OPTIONAL :: virial_xc
360 :
361 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_fxc_calculate'
362 :
363 : INTEGER :: handle, ispin, mspins, nspins
364 : INTEGER, DIMENSION(2, 3) :: bo
365 : LOGICAL :: do_analytic, do_sf, do_triplet, &
366 : do_virial, lsd
367 3930 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: vxg
368 7860 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho0_g, rho1_g
369 11790 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho0_r, rho1_r, tau0_r, tau1_r
370 : TYPE(qs_rho_type), POINTER :: rhot0, rhot1
371 : TYPE(section_vals_type), POINTER :: xc_fun_section
372 : TYPE(xc_derivative_set_type) :: xc_deriv_set
373 : TYPE(xc_rho_cflags_type) :: needs
374 : TYPE(xc_rho_set_type) :: rho0_set, rho1_set
375 :
376 3930 : CALL timeset(routineN, handle)
377 :
378 3930 : do_triplet = .FALSE.
379 3930 : IF (PRESENT(is_triplet)) do_triplet = is_triplet
380 :
381 3930 : do_sf = .FALSE.
382 3930 : IF (PRESENT(spinflip)) do_sf = spinflip
383 :
384 3930 : do_virial = .FALSE.
385 3930 : IF (PRESENT(compute_virial)) do_virial = compute_virial
386 :
387 3930 : do_analytic = section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL")
388 :
389 3930 : CALL qs_rho_get(rho0, rho_r=rho0_r, rho_g=rho0_g, tau_r=tau0_r)
390 3930 : CALL qs_rho_get(rho1, rho_r=rho1_r, tau_r=tau1_r)
391 3930 : NULLIFY (rho1_g)
392 :
393 3930 : mspins = SIZE(rho0_r)
394 3930 : nspins = SIZE(rho0_r)
395 3930 : lsd = (nspins == 2)
396 3930 : IF (nspins == 1 .AND. do_triplet) THEN
397 6 : nspins = 2
398 6 : lsd = .TRUE.
399 3924 : ELSE IF (do_sf) THEN
400 104 : nspins = 1
401 104 : mspins = 1
402 104 : lsd = .TRUE.
403 : END IF
404 :
405 3930 : CPASSERT(.NOT. ASSOCIATED(fxc_rho))
406 3930 : CPASSERT(.NOT. ASSOCIATED(fxc_tau))
407 3930 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
408 3930 : needs = xc_functionals_get_needs(xc_fun_section, lsd, .TRUE.)
409 16148 : ALLOCATE (fxc_rho(mspins))
410 8288 : DO ispin = 1, mspins
411 4358 : CALL auxbas_pw_pool%create_pw(fxc_rho(ispin))
412 8288 : CALL pw_zero(fxc_rho(ispin))
413 : END DO
414 3930 : IF (needs%tau .OR. needs%tau_spin) THEN
415 96 : IF (.NOT. ASSOCIATED(tau1_r)) THEN
416 0 : CPABORT("Tau-dependent functionals requires allocated kinetic energy density grid")
417 : END IF
418 288 : ALLOCATE (fxc_tau(mspins))
419 192 : DO ispin = 1, mspins
420 96 : CALL auxbas_pw_pool%create_pw(fxc_tau(ispin))
421 4026 : CALL pw_zero(fxc_tau(ispin))
422 : END DO
423 : END IF
424 :
425 3930 : IF (mspins == 1 .AND. do_triplet) THEN
426 : ! split the density and response density arrays for triplet calculation
427 6 : ALLOCATE (rhot0)
428 6 : CALL qs_rho_create(rhot0)
429 6 : CALL qs_rho_copy(rho0, rhot0, auxbas_pw_pool, 2, factor=2.0_dp)
430 : !
431 6 : ALLOCATE (rhot1)
432 6 : CALL qs_rho_create(rhot1)
433 6 : CALL qs_rho_copy(rho1, rhot1, auxbas_pw_pool, 2, factor=2.0_dp)
434 : !
435 6 : CALL qs_rho_get(rhot0, rho_r=rho0_r, rho_g=rho0_g, tau_r=tau0_r)
436 6 : CALL qs_rho_get(rhot1, rho_r=rho1_r, tau_r=tau1_r)
437 :
438 : CALL xc_prep_2nd_deriv(xc_deriv_set, rho0_set, rho0_r, auxbas_pw_pool, weights, &
439 6 : xc_section=xc_section, tau_r=tau0_r)
440 60 : bo = rho1_r(1)%pw_grid%bounds_local
441 : ! create the place where to store the argument for the functionals
442 : CALL xc_rho_set_create(rho1_set, bo, &
443 : rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
444 : drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
445 6 : tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
446 :
447 : ! calculate the arguments needed by the functionals
448 : CALL xc_rho_set_update(rho1_set, rho1_r, rho1_g, tau1_r, needs, &
449 : section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
450 : section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
451 6 : auxbas_pw_pool, spinflip=do_sf)
452 : ELSE
453 : CALL xc_prep_2nd_deriv(xc_deriv_set, rho0_set, rho0_r, auxbas_pw_pool, weights, &
454 3924 : xc_section=xc_section, tau_r=tau0_r)
455 39240 : bo = rho1_r(1)%pw_grid%bounds_local
456 : ! create the place where to store the argument for the functionals
457 : CALL xc_rho_set_create(rho1_set, bo, &
458 : rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
459 : drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
460 3924 : tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
461 :
462 : ! calculate the arguments needed by the functionals
463 : CALL xc_rho_set_update(rho1_set, rho1_r, rho1_g, tau1_r, needs, &
464 : section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
465 : section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
466 3924 : auxbas_pw_pool, spinflip=do_sf)
467 : END IF
468 :
469 3930 : IF (mspins == 1 .AND. do_triplet .AND. do_analytic) THEN
470 :
471 : CALL xc_calc_2nd_deriv_analytical(fxc_rho, fxc_tau, xc_deriv_set, rho0_set, &
472 : rho1_set, auxbas_pw_pool, xc_section, &
473 : gapw=.FALSE., vxg=vxg, tddfpt_fac=-1.0_dp, spinflip=do_sf, &
474 6 : compute_virial=compute_virial, virial_xc=virial_xc)
475 :
476 3924 : ELSE IF (do_analytic) THEN
477 :
478 : CALL xc_calc_2nd_deriv_analytical(fxc_rho, fxc_tau, xc_deriv_set, rho0_set, &
479 : rho1_set, auxbas_pw_pool, xc_section, &
480 : gapw=.FALSE., vxg=vxg, spinflip=do_sf, &
481 3856 : compute_virial=compute_virial, virial_xc=virial_xc)
482 :
483 : ELSE
484 :
485 : CALL xc_calc_2nd_deriv_numerical(fxc_rho, fxc_tau, rho0_set, rho1_r, rho1_g, tau1_r, &
486 : auxbas_pw_pool, weights, xc_section, &
487 68 : do_triplet, compute_virial, virial_xc, xc_deriv_set)
488 :
489 : END IF
490 :
491 3930 : IF (mspins == 1 .AND. do_triplet) THEN
492 6 : CALL qs_rho_release(rhot0)
493 6 : DEALLOCATE (rhot0)
494 6 : CALL qs_rho_release(rhot1)
495 6 : DEALLOCATE (rhot1)
496 : END IF
497 :
498 3930 : CALL xc_dset_release(xc_deriv_set)
499 3930 : CALL xc_rho_set_release(rho0_set)
500 3930 : CALL xc_rho_set_release(rho1_set)
501 :
502 3930 : CALL timestop(handle)
503 :
504 168990 : END SUBROUTINE qs_fxc_calculate
505 :
506 : ! **************************************************************************************************
507 : !> \brief ...
508 : !> \param qs_env ...
509 : !> \param rho0_struct ...
510 : !> \param xc_rho_set ...
511 : !> \param xc_deriv_set ...
512 : !> \param xc_section ...
513 : !> \param pw_env_ext ...
514 : !> \param is_triplet ...
515 : ! **************************************************************************************************
516 2820 : SUBROUTINE qs_fxc_prep(qs_env, rho0_struct, xc_rho_set, xc_deriv_set, &
517 : xc_section, pw_env_ext, is_triplet)
518 :
519 : TYPE(qs_environment_type), POINTER :: qs_env
520 : TYPE(qs_rho_type), POINTER :: rho0_struct
521 : TYPE(xc_rho_set_type) :: xc_rho_set
522 : TYPE(xc_derivative_set_type) :: xc_deriv_set
523 : TYPE(section_vals_type), POINTER :: xc_section
524 : TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_ext
525 : LOGICAL, INTENT(IN), OPTIONAL :: is_triplet
526 :
527 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_fxc_prep'
528 :
529 : INTEGER :: handle, ispin, nspins
530 : LOGICAL :: uf_grid
531 : REAL(KIND=dp) :: factor
532 : TYPE(dft_control_type), POINTER :: dft_control
533 2820 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
534 : TYPE(pw_c1d_gs_type), POINTER :: rho_nlcc_g
535 : TYPE(pw_env_type), POINTER :: pw_env
536 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool, xc_pw_pool
537 2820 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
538 : TYPE(pw_r3d_rs_type), POINTER :: rho_nlcc, weights, weights_uf
539 : TYPE(qs_rho_type), POINTER :: rho0_uf
540 :
541 2820 : CALL timeset(routineN, handle)
542 :
543 2820 : CALL get_qs_env(qs_env, dft_control=dft_control)
544 :
545 2820 : IF (PRESENT(pw_env_ext)) THEN
546 2820 : pw_env => pw_env_ext
547 : ELSE
548 0 : CALL get_qs_env(qs_env, pw_env=pw_env)
549 : END IF
550 2820 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, xc_pw_pool=xc_pw_pool)
551 2820 : uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid, xc_pw_pool%pw_grid)
552 :
553 2820 : nspins = dft_control%nspins
554 :
555 2820 : NULLIFY (rho_nlcc, rho_nlcc_g)
556 2820 : CALL get_qs_env(qs_env, rho_nlcc=rho_nlcc, rho_nlcc_g=rho_nlcc_g)
557 2820 : IF (ASSOCIATED(rho_nlcc)) THEN
558 20 : NULLIFY (rho_r, rho_g)
559 20 : CALL qs_rho_get(rho0_struct, rho_r=rho_r, rho_g=rho_g)
560 20 : factor = 1.0_dp
561 40 : DO ispin = 1, nspins
562 20 : CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
563 40 : CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
564 : END DO
565 : END IF
566 :
567 2820 : NULLIFY (weights)
568 2820 : CALL get_qs_env(qs_env, xcint_weights=weights)
569 :
570 2820 : IF (uf_grid) THEN
571 30 : NULLIFY (weights_uf)
572 30 : IF (ASSOCIATED(weights)) THEN
573 30 : ALLOCATE (weights_uf)
574 30 : CALL xc_pw_pool%create_pw(weights_uf)
575 : BLOCK
576 : TYPE(pw_c1d_gs_type) :: weights_g, weights_g_uf
577 30 : CALL auxbas_pw_pool%create_pw(weights_g)
578 30 : CALL xc_pw_pool%create_pw(weights_g_uf)
579 30 : CALL pw_transfer(weights, weights_g)
580 30 : CALL pw_transfer(weights_g, weights_g_uf)
581 30 : CALL pw_transfer(weights_g_uf, weights_uf)
582 30 : CALL xc_pw_pool%give_back_pw(weights_g_uf)
583 60 : CALL auxbas_pw_pool%give_back_pw(weights_g)
584 : END BLOCK
585 : END IF
586 : !
587 30 : ALLOCATE (rho0_uf)
588 30 : CALL qs_rho_create(rho0_uf)
589 30 : CALL qs_rho_transfer(rho0_struct, rho0_uf, auxbas_pw_pool, xc_pw_pool)
590 : !
591 : CALL qs_fxc_deriv(rho0_uf, xc_rho_set, xc_deriv_set, &
592 30 : xc_section, weights_uf, xc_pw_pool, is_triplet)
593 : !
594 30 : CALL qs_rho_release(rho0_uf)
595 30 : DEALLOCATE (rho0_uf)
596 30 : IF (ASSOCIATED(weights_uf)) THEN
597 30 : CALL xc_pw_pool%give_back_pw(weights_uf)
598 30 : DEALLOCATE (weights_uf)
599 : END IF
600 : ELSE
601 : CALL qs_fxc_deriv(rho0_struct, xc_rho_set, xc_deriv_set, &
602 2790 : xc_section, weights, auxbas_pw_pool, is_triplet)
603 : END IF
604 :
605 : ! de-apply NLCC density
606 2820 : IF (ASSOCIATED(rho_nlcc)) THEN
607 20 : factor = -1.0_dp
608 40 : DO ispin = 1, nspins
609 20 : CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
610 40 : CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
611 : END DO
612 : END IF
613 :
614 2820 : CALL timestop(handle)
615 :
616 2820 : END SUBROUTINE qs_fxc_prep
617 :
618 : ! **************************************************************************************************
619 : !> \brief ...
620 : !> \param rho0 ...
621 : !> \param xc_rho_set ...
622 : !> \param xc_deriv_set ...
623 : !> \param xc_section ...
624 : !> \param weights ...
625 : !> \param auxbas_pw_pool ...
626 : !> \param is_triplet ...
627 : ! **************************************************************************************************
628 2820 : SUBROUTINE qs_fxc_deriv(rho0, xc_rho_set, xc_deriv_set, xc_section, weights, auxbas_pw_pool, &
629 : is_triplet)
630 :
631 : TYPE(qs_rho_type), POINTER :: rho0
632 : TYPE(xc_rho_set_type) :: xc_rho_set
633 : TYPE(xc_derivative_set_type) :: xc_deriv_set
634 : TYPE(section_vals_type), POINTER :: xc_section
635 : TYPE(pw_r3d_rs_type), POINTER :: weights
636 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
637 : LOGICAL, INTENT(IN) :: is_triplet
638 :
639 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_fxc_deriv'
640 :
641 : INTEGER :: handle
642 2820 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho0_r, tau0_r
643 : TYPE(qs_rho_type), POINTER :: rhot0
644 :
645 2820 : CALL timeset(routineN, handle)
646 :
647 2820 : NULLIFY (rho0_r, tau0_r)
648 2820 : IF (is_triplet) THEN
649 : ! split the density and response density arrays for triplet calculation
650 110 : ALLOCATE (rhot0)
651 110 : CALL qs_rho_create(rhot0)
652 110 : CALL qs_rho_copy(rho0, rhot0, auxbas_pw_pool, 2, factor=2.0_dp)
653 110 : CALL qs_rho_get(rhot0, rho_r=rho0_r, tau_r=tau0_r)
654 : CALL xc_prep_2nd_deriv(xc_deriv_set, xc_rho_set, rho0_r, auxbas_pw_pool, weights, &
655 110 : xc_section=xc_section, tau_r=tau0_r)
656 110 : CALL qs_rho_release(rhot0)
657 110 : DEALLOCATE (rhot0)
658 : ELSE
659 2710 : CALL qs_rho_get(rho0, rho_r=rho0_r, tau_r=tau0_r)
660 : CALL xc_prep_2nd_deriv(xc_deriv_set, xc_rho_set, rho0_r, auxbas_pw_pool, weights, &
661 2710 : xc_section=xc_section, tau_r=tau0_r)
662 : END IF
663 :
664 2820 : CALL timestop(handle)
665 :
666 2820 : END SUBROUTINE qs_fxc_deriv
667 :
668 : ! **************************************************************************************************
669 : !> \brief ...
670 : !> \param qs_env ...
671 : !> \param xc_deriv_set ...
672 : !> \param xc_rho_set ...
673 : !> \param rho1_struct ...
674 : !> \param rho0_atom_set ...
675 : !> \param xc_section ...
676 : !> \param do_onecenter ...
677 : !> \param fxc_rho ...
678 : !> \param fxc_tau ...
679 : !> \param rho1_atom_set ...
680 : !> \param do_scale ...
681 : !> \param is_triplet ...
682 : !> \param spinflip ...
683 : !> \param pw_env_ext ...
684 : !> \param kind_set_external ...
685 : !> \param para_env_external ...
686 : !> \param compute_virial ...
687 : !> \param virial_xc ...
688 : ! **************************************************************************************************
689 22294 : SUBROUTINE qs_fxc_apply(qs_env, xc_deriv_set, xc_rho_set, rho1_struct, rho0_atom_set, &
690 : xc_section, do_onecenter, fxc_rho, fxc_tau, rho1_atom_set, &
691 : do_scale, is_triplet, spinflip, pw_env_ext, &
692 : kind_set_external, para_env_external, compute_virial, virial_xc)
693 :
694 : TYPE(qs_environment_type), POINTER :: qs_env
695 : TYPE(xc_derivative_set_type) :: xc_deriv_set
696 : TYPE(xc_rho_set_type) :: xc_rho_set
697 : TYPE(qs_rho_type), POINTER :: rho1_struct
698 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho0_atom_set
699 : TYPE(section_vals_type), POINTER :: xc_section
700 : LOGICAL, INTENT(IN) :: do_onecenter
701 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau
702 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho1_atom_set
703 : LOGICAL, INTENT(IN), OPTIONAL :: do_scale, is_triplet, spinflip
704 : TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_ext
705 : TYPE(qs_kind_type), DIMENSION(:), OPTIONAL, &
706 : POINTER :: kind_set_external
707 : TYPE(mp_para_env_type), INTENT(IN), OPTIONAL :: para_env_external
708 : LOGICAL, INTENT(IN), OPTIONAL :: compute_virial
709 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), &
710 : OPTIONAL :: virial_xc
711 :
712 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_fxc_apply'
713 :
714 : INTEGER :: handle, ispin, nspins
715 : LOGICAL :: do_virial, uf_grid
716 : TYPE(dft_control_type), POINTER :: dft_control
717 : TYPE(pw_env_type), POINTER :: pw_env
718 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool, xc_pw_pool
719 22294 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho_lo, fxc_rho_uf, fxc_tau_lo, &
720 22294 : fxc_tau_uf
721 : TYPE(pw_r3d_rs_type), POINTER :: weights, weights_uf
722 : TYPE(qs_rho_type), POINTER :: rho1_uf
723 :
724 22294 : CALL timeset(routineN, handle)
725 :
726 22294 : do_virial = .FALSE.
727 22294 : IF (PRESENT(compute_virial)) do_virial = compute_virial
728 :
729 22294 : CALL get_qs_env(qs_env, dft_control=dft_control)
730 :
731 22294 : IF (PRESENT(pw_env_ext)) THEN
732 9592 : pw_env => pw_env_ext
733 : ELSE
734 12702 : CALL get_qs_env(qs_env, pw_env=pw_env)
735 : END IF
736 22294 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, xc_pw_pool=xc_pw_pool)
737 22294 : uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid, xc_pw_pool%pw_grid)
738 :
739 22294 : nspins = dft_control%nspins
740 22294 : IF (ASSOCIATED(fxc_rho)) THEN
741 9592 : CPASSERT(nspins == SIZE(fxc_rho))
742 : END IF
743 22294 : IF (ASSOCIATED(fxc_tau)) THEN
744 9592 : CPASSERT(nspins == SIZE(fxc_tau))
745 : END IF
746 :
747 22294 : CALL get_qs_env(qs_env, xcint_weights=weights)
748 :
749 22294 : NULLIFY (fxc_rho_lo, fxc_tau_lo)
750 22294 : IF (uf_grid) THEN
751 198 : NULLIFY (weights_uf)
752 198 : IF (ASSOCIATED(weights)) THEN
753 198 : ALLOCATE (weights_uf)
754 198 : CALL xc_pw_pool%create_pw(weights_uf)
755 : BLOCK
756 : TYPE(pw_c1d_gs_type) :: weights_g, weights_g_uf
757 198 : CALL auxbas_pw_pool%create_pw(weights_g)
758 198 : CALL xc_pw_pool%create_pw(weights_g_uf)
759 198 : CALL pw_transfer(weights, weights_g)
760 198 : CALL pw_transfer(weights_g, weights_g_uf)
761 198 : CALL pw_transfer(weights_g_uf, weights_uf)
762 198 : CALL xc_pw_pool%give_back_pw(weights_g_uf)
763 396 : CALL auxbas_pw_pool%give_back_pw(weights_g)
764 : END BLOCK
765 : END IF
766 : !
767 198 : ALLOCATE (rho1_uf)
768 198 : CALL qs_rho_create(rho1_uf)
769 198 : CALL qs_rho_transfer(rho1_struct, rho1_uf, auxbas_pw_pool, xc_pw_pool)
770 : !
771 198 : NULLIFY (fxc_rho_uf, fxc_tau_uf)
772 : CALL qs_fxc_eval(xc_deriv_set, xc_rho_set, rho1_uf, xc_section, &
773 : weights_uf, xc_pw_pool, fxc_rho_uf, fxc_tau_uf, &
774 : is_triplet=is_triplet, spinflip=spinflip, &
775 198 : compute_virial=do_virial, virial_xc=virial_xc)
776 : !
777 198 : CALL qs_rho_release(rho1_uf)
778 198 : DEALLOCATE (rho1_uf)
779 198 : IF (ASSOCIATED(weights_uf)) THEN
780 198 : CALL xc_pw_pool%give_back_pw(weights_uf)
781 198 : DEALLOCATE (weights_uf)
782 : END IF
783 198 : IF (ASSOCIATED(fxc_rho_uf)) THEN
784 792 : ALLOCATE (fxc_rho_lo(nspins))
785 396 : DO ispin = 1, nspins
786 198 : CALL auxbas_pw_pool%create_pw(fxc_rho_lo(ispin))
787 : BLOCK
788 : TYPE(pw_c1d_gs_type) :: fxc_g, fxc_g_uf
789 198 : CALL auxbas_pw_pool%create_pw(fxc_g)
790 198 : CALL xc_pw_pool%create_pw(fxc_g_uf)
791 198 : CALL pw_transfer(fxc_rho_uf(ispin), fxc_g_uf)
792 198 : CALL pw_transfer(fxc_g_uf, fxc_g)
793 198 : CALL pw_transfer(fxc_g, fxc_rho_lo(ispin))
794 198 : CALL xc_pw_pool%give_back_pw(fxc_g_uf)
795 396 : CALL auxbas_pw_pool%give_back_pw(fxc_g)
796 : END BLOCK
797 396 : CALL xc_pw_pool%give_back_pw(fxc_rho_uf(ispin))
798 : END DO
799 198 : DEALLOCATE (fxc_rho_uf)
800 : END IF
801 198 : IF (ASSOCIATED(fxc_tau_uf)) THEN
802 0 : ALLOCATE (fxc_tau_lo(nspins))
803 0 : DO ispin = 1, nspins
804 0 : CALL auxbas_pw_pool%create_pw(fxc_tau_lo(ispin))
805 : BLOCK
806 : TYPE(pw_c1d_gs_type) :: fxc_g, fxc_g_uf
807 0 : CALL auxbas_pw_pool%create_pw(fxc_g)
808 0 : CALL xc_pw_pool%create_pw(fxc_g_uf)
809 0 : CALL pw_transfer(fxc_tau_uf(ispin), fxc_g_uf)
810 0 : CALL pw_transfer(fxc_g_uf, fxc_g)
811 0 : CALL pw_transfer(fxc_g, fxc_tau_lo(ispin))
812 0 : CALL xc_pw_pool%give_back_pw(fxc_g_uf)
813 0 : CALL auxbas_pw_pool%give_back_pw(fxc_g)
814 : END BLOCK
815 0 : CALL xc_pw_pool%give_back_pw(fxc_tau_uf(ispin))
816 : END DO
817 : END IF
818 :
819 : ELSE
820 : CALL qs_fxc_eval(xc_deriv_set, xc_rho_set, rho1_struct, xc_section, &
821 : weights, auxbas_pw_pool, fxc_rho_lo, fxc_tau_lo, &
822 : is_triplet=is_triplet, spinflip=spinflip, &
823 22096 : compute_virial=do_virial, virial_xc=virial_xc)
824 : END IF
825 :
826 : ! return potentials
827 22294 : IF (ASSOCIATED(fxc_rho)) THEN
828 21006 : DO ispin = 1, MIN(SIZE(fxc_rho_lo), SIZE(fxc_rho))
829 21006 : CALL pw_transfer(fxc_rho_lo(ispin), fxc_rho(ispin))
830 : END DO
831 21006 : DO ispin = 1, SIZE(fxc_rho_lo)
832 21006 : CALL auxbas_pw_pool%give_back_pw(fxc_rho_lo(ispin))
833 : END DO
834 9592 : DEALLOCATE (fxc_rho_lo)
835 : ELSE
836 12702 : fxc_rho => fxc_rho_lo
837 : END IF
838 22294 : IF (ASSOCIATED(fxc_tau)) THEN
839 9592 : IF (ASSOCIATED(fxc_tau_lo)) THEN
840 440 : DO ispin = 1, MIN(SIZE(fxc_tau_lo), SIZE(fxc_tau))
841 440 : CALL pw_transfer(fxc_tau_lo(ispin), fxc_tau(ispin))
842 : END DO
843 440 : DO ispin = 1, SIZE(fxc_tau_lo)
844 440 : CALL auxbas_pw_pool%give_back_pw(fxc_tau_lo(ispin))
845 : END DO
846 220 : DEALLOCATE (fxc_tau_lo)
847 : ELSE
848 20876 : DO ispin = 1, nspins
849 20876 : CALL pw_zero(fxc_tau(ispin))
850 : END DO
851 : END IF
852 : ELSE
853 12702 : fxc_tau => fxc_tau_lo
854 : END IF
855 :
856 22294 : IF (do_onecenter) THEN
857 : CALL fxc_atom_calc(qs_env, rho0_atom_set, rho1_atom_set, xc_section, &
858 : do_scale=do_scale, do_triplet=is_triplet, do_sf=spinflip, &
859 : para_env_ext=para_env_external, &
860 5746 : kind_set_external=kind_set_external)
861 : END IF
862 :
863 22294 : CALL timestop(handle)
864 :
865 22294 : END SUBROUTINE qs_fxc_apply
866 :
867 : ! **************************************************************************************************
868 : !> \brief ...
869 : !> \param xc_deriv_set ...
870 : !> \param xc_rho_set ...
871 : !> \param rho1 ...
872 : !> \param xc_section ...
873 : !> \param weights ...
874 : !> \param auxbas_pw_pool ...
875 : !> \param fxc_rho ...
876 : !> \param fxc_tau ...
877 : !> \param is_triplet ...
878 : !> \param spinflip ...
879 : !> \param compute_virial ...
880 : !> \param virial_xc ...
881 : ! **************************************************************************************************
882 22294 : SUBROUTINE qs_fxc_eval(xc_deriv_set, xc_rho_set, rho1, xc_section, weights, auxbas_pw_pool, &
883 : fxc_rho, fxc_tau, is_triplet, spinflip, &
884 : compute_virial, virial_xc)
885 :
886 : TYPE(xc_derivative_set_type) :: xc_deriv_set
887 : TYPE(xc_rho_set_type) :: xc_rho_set
888 : TYPE(qs_rho_type), POINTER :: rho1
889 : TYPE(section_vals_type), POINTER :: xc_section
890 : TYPE(pw_r3d_rs_type), POINTER :: weights
891 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
892 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau
893 : LOGICAL, INTENT(IN), OPTIONAL :: is_triplet, spinflip, compute_virial
894 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), &
895 : OPTIONAL :: virial_xc
896 :
897 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_fxc_eval'
898 :
899 : INTEGER :: handle, ispin, mspins, nspins
900 : INTEGER, DIMENSION(2, 3) :: bo
901 : LOGICAL :: do_analytic, do_sf, do_triplet, &
902 : do_virial, lsd
903 22294 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: vxg
904 22294 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho1_g
905 44588 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r, tau1_r
906 : TYPE(qs_rho_type), POINTER :: rhot1
907 : TYPE(section_vals_type), POINTER :: xc_fun_section
908 : TYPE(xc_rho_cflags_type) :: needs
909 : TYPE(xc_rho_set_type) :: rho1_set
910 :
911 22294 : CALL timeset(routineN, handle)
912 :
913 22294 : do_triplet = .FALSE.
914 22294 : IF (PRESENT(is_triplet)) do_triplet = is_triplet
915 :
916 22294 : do_sf = .FALSE.
917 22294 : IF (PRESENT(spinflip)) do_sf = spinflip
918 :
919 22294 : do_virial = .FALSE.
920 22294 : IF (PRESENT(compute_virial)) do_virial = compute_virial
921 :
922 22294 : do_analytic = section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL")
923 :
924 22294 : CALL qs_rho_get(rho1, rho_r=rho1_r, tau_r=tau1_r)
925 22294 : NULLIFY (rho1_g)
926 :
927 22294 : mspins = SIZE(rho1_r)
928 22294 : nspins = SIZE(rho1_r)
929 22294 : lsd = (nspins == 2)
930 22294 : IF (nspins == 1 .AND. do_triplet) THEN
931 652 : nspins = 2
932 652 : lsd = .TRUE.
933 21642 : ELSE IF (do_sf) THEN
934 310 : nspins = 1
935 310 : mspins = 1
936 310 : lsd = .TRUE.
937 : END IF
938 :
939 22294 : CPASSERT(.NOT. ASSOCIATED(fxc_rho))
940 22294 : CPASSERT(.NOT. ASSOCIATED(fxc_tau))
941 22294 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
942 22294 : needs = xc_functionals_get_needs(xc_fun_section, lsd, .TRUE.)
943 92742 : ALLOCATE (fxc_rho(mspins))
944 48154 : DO ispin = 1, mspins
945 25860 : CALL auxbas_pw_pool%create_pw(fxc_rho(ispin))
946 48154 : CALL pw_zero(fxc_rho(ispin))
947 : END DO
948 22294 : IF (needs%tau .OR. needs%tau_spin) THEN
949 718 : IF (.NOT. ASSOCIATED(tau1_r)) THEN
950 0 : CPABORT("Tau-dependent functionals requires allocated kinetic energy density grid")
951 : END IF
952 2282 : ALLOCATE (fxc_tau(mspins))
953 1564 : DO ispin = 1, mspins
954 846 : CALL auxbas_pw_pool%create_pw(fxc_tau(ispin))
955 23140 : CALL pw_zero(fxc_tau(ispin))
956 : END DO
957 : END IF
958 :
959 22294 : IF (mspins == 1 .AND. do_triplet) THEN
960 : ! split the density and response density arrays for triplet calculation
961 652 : ALLOCATE (rhot1)
962 652 : CALL qs_rho_create(rhot1)
963 652 : CALL qs_rho_copy(rho1, rhot1, auxbas_pw_pool, 2, factor=2.0_dp)
964 : !
965 652 : CALL qs_rho_get(rhot1, rho_r=rho1_r, tau_r=tau1_r)
966 : END IF
967 :
968 222940 : bo = rho1_r(1)%pw_grid%bounds_local
969 : ! create the place where to store the argument for the functionals
970 : CALL xc_rho_set_create(rho1_set, bo, &
971 : rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
972 : drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
973 22294 : tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
974 :
975 : ! calculate the arguments needed by the functionals
976 : CALL xc_rho_set_update(rho1_set, rho1_r, rho1_g, tau1_r, needs, &
977 : section_get_ival(xc_section, "XC_GRID%XC_DERIV"), &
978 : section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO"), &
979 22294 : auxbas_pw_pool, spinflip=do_sf)
980 :
981 22294 : IF (mspins == 1 .AND. do_triplet .AND. do_analytic) THEN
982 :
983 : CALL xc_calc_2nd_deriv_analytical(fxc_rho, fxc_tau, xc_deriv_set, xc_rho_set, &
984 : rho1_set, auxbas_pw_pool, xc_section, &
985 : gapw=.FALSE., vxg=vxg, tddfpt_fac=-1.0_dp, spinflip=do_sf, &
986 632 : compute_virial=compute_virial, virial_xc=virial_xc)
987 :
988 21642 : ELSE IF (do_analytic) THEN
989 :
990 : CALL xc_calc_2nd_deriv_analytical(fxc_rho, fxc_tau, xc_deriv_set, xc_rho_set, &
991 : rho1_set, auxbas_pw_pool, xc_section, &
992 : gapw=.FALSE., vxg=vxg, spinflip=do_sf, &
993 21236 : compute_virial=compute_virial, virial_xc=virial_xc)
994 :
995 : ELSE
996 :
997 : CALL xc_calc_2nd_deriv_numerical(fxc_rho, fxc_tau, xc_rho_set, rho1_r, rho1_g, tau1_r, &
998 : auxbas_pw_pool, weights, xc_section, &
999 426 : do_triplet, compute_virial, virial_xc, xc_deriv_set)
1000 :
1001 : END IF
1002 :
1003 22294 : IF (mspins == 1 .AND. do_triplet) THEN
1004 652 : CALL qs_rho_release(rhot1)
1005 652 : DEALLOCATE (rhot1)
1006 : END IF
1007 22294 : CALL xc_rho_set_release(rho1_set)
1008 :
1009 22294 : CALL timestop(handle)
1010 :
1011 490468 : END SUBROUTINE qs_fxc_eval
1012 :
1013 : ! **************************************************************************************************
1014 : !> \brief ...
1015 : !> \param qs_env ...
1016 : !> \param rho0_struct ...
1017 : !> \param rho1_struct ...
1018 : !> \param xc_section ...
1019 : !> \param accuracy ...
1020 : !> \param fxc_rho ...
1021 : !> \param fxc_tau ...
1022 : !> \param is_triplet ...
1023 : !> \param spinflip ...
1024 : ! **************************************************************************************************
1025 2282 : SUBROUTINE qs_fxc_fdiff(qs_env, rho0_struct, rho1_struct, xc_section, accuracy, &
1026 : fxc_rho, fxc_tau, is_triplet, spinflip)
1027 :
1028 : TYPE(qs_environment_type), POINTER :: qs_env
1029 : TYPE(qs_rho_type), POINTER :: rho0_struct, rho1_struct
1030 : TYPE(section_vals_type), POINTER :: xc_section
1031 : INTEGER, INTENT(IN) :: accuracy
1032 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau
1033 : LOGICAL, INTENT(IN), OPTIONAL :: is_triplet, spinflip
1034 :
1035 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_fxc_fdiff'
1036 : REAL(KIND=dp), PARAMETER :: epsrho = 5.e-4_dp
1037 :
1038 : INTEGER :: handle, ispin, istep, nspins, nstep
1039 : LOGICAL :: do_sf, do_triplet
1040 : REAL(KIND=dp) :: alpha, beta, exc, oeps1
1041 : REAL(KIND=dp), DIMENSION(-4:4) :: ak
1042 : TYPE(dft_control_type), POINTER :: dft_control
1043 : TYPE(pw_env_type), POINTER :: pw_env
1044 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1045 2282 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_tau_rspace, vxc00
1046 : TYPE(qs_ks_env_type), POINTER :: ks_env
1047 : TYPE(qs_rho_type), POINTER :: rhoin
1048 :
1049 2282 : CALL timeset(routineN, handle)
1050 :
1051 2282 : CPASSERT(.NOT. ASSOCIATED(fxc_rho))
1052 2282 : CPASSERT(.NOT. ASSOCIATED(fxc_tau))
1053 2282 : CPASSERT(ASSOCIATED(rho0_struct))
1054 2282 : CPASSERT(ASSOCIATED(rho1_struct))
1055 :
1056 2282 : do_triplet = .FALSE.
1057 2282 : IF (PRESENT(is_triplet)) do_triplet = is_triplet
1058 :
1059 2282 : do_sf = .FALSE.
1060 2282 : IF (PRESENT(spinflip)) do_sf = spinflip
1061 0 : IF (do_sf) THEN
1062 0 : CPABORT("Spin Flip TDDFT only available with analytic 2nd xc derivatives")
1063 : END IF
1064 :
1065 2282 : ak = 0.0_dp
1066 2282 : SELECT CASE (accuracy)
1067 : CASE (:4)
1068 0 : nstep = 2
1069 0 : ak(-2:2) = [1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp]/12.0_dp
1070 : CASE (5:7)
1071 18256 : nstep = 3
1072 18256 : 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
1073 : CASE (8:)
1074 0 : nstep = 4
1075 : ak(-4:4) = [1.0_dp, -32.0_dp/3.0_dp, 56.0_dp, -224.0_dp, 0.0_dp, &
1076 2282 : 224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp]/280.0_dp
1077 : END SELECT
1078 :
1079 2282 : CALL get_qs_env(qs_env, ks_env=ks_env, dft_control=dft_control, pw_env=pw_env)
1080 2282 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1081 :
1082 2282 : nspins = dft_control%nspins
1083 : exc = 0.0_dp
1084 :
1085 18256 : DO istep = -nstep, nstep
1086 :
1087 18256 : IF (ak(istep) /= 0.0_dp) THEN
1088 13692 : alpha = 1.0_dp
1089 13692 : beta = REAL(istep, KIND=dp)*epsrho
1090 : NULLIFY (rhoin)
1091 13692 : ALLOCATE (rhoin)
1092 13692 : CALL qs_rho_create(rhoin)
1093 13692 : NULLIFY (vxc00, v_tau_rspace)
1094 13692 : IF (do_triplet) THEN
1095 1176 : CPASSERT(nspins == 1)
1096 : ! rhoin = (0.5 rho0, 0.5 rho0)
1097 1176 : CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, 2)
1098 : ! rhoin = (0.5 rho0 + 0.5 rho1, 0.5 rho0)
1099 1176 : CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, 0.5_dp*beta)
1100 1176 : CALL qs_vxc_create(ks_env, rhoin, xc_section, vxc00, v_tau_rspace, exc)
1101 1176 : CALL pw_axpy(vxc00(2), vxc00(1), -1.0_dp)
1102 1176 : IF (ASSOCIATED(v_tau_rspace)) CALL pw_axpy(v_tau_rspace(2), v_tau_rspace(1), -1.0_dp)
1103 : ELSE
1104 12516 : CALL qs_rho_copy(rho0_struct, rhoin, auxbas_pw_pool, nspins)
1105 12516 : CALL qs_rho_scale_and_add(rhoin, rho1_struct, alpha, beta)
1106 12516 : CALL qs_vxc_create(ks_env, rhoin, xc_section, vxc00, v_tau_rspace, exc)
1107 : END IF
1108 13692 : CALL qs_rho_release(rhoin)
1109 13692 : DEALLOCATE (rhoin)
1110 13692 : IF (.NOT. ASSOCIATED(fxc_rho)) THEN
1111 9380 : ALLOCATE (fxc_rho(nspins))
1112 4816 : DO ispin = 1, nspins
1113 2534 : CALL auxbas_pw_pool%create_pw(fxc_rho(ispin))
1114 4816 : CALL pw_zero(fxc_rho(ispin))
1115 : END DO
1116 : END IF
1117 28896 : DO ispin = 1, nspins
1118 28896 : CALL pw_axpy(vxc00(ispin), fxc_rho(ispin), ak(istep))
1119 : END DO
1120 30072 : DO ispin = 1, SIZE(vxc00)
1121 30072 : CALL auxbas_pw_pool%give_back_pw(vxc00(ispin))
1122 : END DO
1123 13692 : DEALLOCATE (vxc00)
1124 13692 : IF (ASSOCIATED(v_tau_rspace)) THEN
1125 0 : IF (.NOT. ASSOCIATED(fxc_tau)) THEN
1126 0 : ALLOCATE (fxc_tau(nspins))
1127 0 : DO ispin = 1, nspins
1128 0 : CALL auxbas_pw_pool%create_pw(fxc_tau(ispin))
1129 0 : CALL pw_zero(fxc_tau(ispin))
1130 : END DO
1131 : END IF
1132 0 : DO ispin = 1, nspins
1133 0 : CALL pw_axpy(v_tau_rspace(ispin), fxc_tau(ispin), ak(istep))
1134 : END DO
1135 0 : DO ispin = 1, SIZE(v_tau_rspace)
1136 0 : CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
1137 : END DO
1138 0 : DEALLOCATE (v_tau_rspace)
1139 : END IF
1140 : END IF
1141 :
1142 : END DO
1143 :
1144 2282 : oeps1 = 1.0_dp/epsrho
1145 4816 : DO ispin = 1, nspins
1146 4816 : CALL pw_scale(fxc_rho(ispin), oeps1)
1147 : END DO
1148 2282 : IF (ASSOCIATED(fxc_tau)) THEN
1149 0 : DO ispin = 1, nspins
1150 0 : CALL pw_scale(fxc_tau(ispin), oeps1)
1151 : END DO
1152 : END IF
1153 :
1154 2282 : CALL timestop(handle)
1155 :
1156 2282 : END SUBROUTINE qs_fxc_fdiff
1157 :
1158 : END MODULE qs_fxc
|
