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