(git:a660c7f)
Loading...
Searching...
No Matches
qs_vxc.F
Go to the documentation of this file.
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
10!>
11!>
12!> \par History
13!> refactoring 03-2011 [MI]
14!> \author MI
15! **************************************************************************************************
16MODULE qs_vxc
17
18 USE cell_types, ONLY: cell_type
20 USE input_constants, ONLY: sic_ad,&
21 sic_eo,&
24 sic_none,&
25 xc_none,&
29 USE kinds, ONLY: dp
32 USE pw_env_types, ONLY: pw_env_get,&
34 USE pw_grids, ONLY: pw_grid_compare
35 USE pw_methods, ONLY: pw_axpy,&
36 pw_copy,&
38 pw_scale,&
42 USE pw_types, ONLY: pw_c1d_gs_type,&
46 USE qs_ks_types, ONLY: get_ks_env,&
48 USE qs_rho_types, ONLY: qs_rho_get,&
52 USE virial_types, ONLY: virial_type
53 USE xc, ONLY: calc_xc_density,&
56#include "./base/base_uses.f90"
57
58 IMPLICIT NONE
59
60 PRIVATE
61
62 ! *** Public subroutines ***
64
65 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_vxc'
66
67CONTAINS
68
69! **************************************************************************************************
70!> \brief calculates and allocates the xc potential, already reducing it to
71!> the dependence on rho and the one on tau
72!> \param ks_env to get all the needed things
73!> \param rho_struct density for which v_xc is calculated
74!> \param xc_section ...
75!> \param vxc_rho will contain the v_xc part that depend on rho
76!> (if one of the chosen xc functionals has it it is allocated and you
77!> are responsible for it)
78!> \param vxc_tau will contain the kinetic tau part of v_xc
79!> (if one of the chosen xc functionals has it it is allocated and you
80!> are responsible for it)
81!> \param exc ...
82!> \param just_energy if true calculates just the energy, and does not
83!> allocate v_*_rspace
84!> \param edisp ...
85!> \param dispersion_env ...
86!> \param adiabatic_rescale_factor ...
87!> \param pw_env_external external plane wave environment
88!> \param native_skala_atom_force ...
89!> \par History
90!> - 05.2002 modified to use the mp_allgather function each pe
91!> computes only part of the grid and this is broadcasted to all
92!> instead of summed.
93!> This scales significantly better (e.g. factor 3 on 12 cpus
94!> 32 H2O) [Joost VdV]
95!> - moved to qs_ks_methods [fawzi]
96!> - sic alterations [Joost VandeVondele]
97!> \author Fawzi Mohamed
98! **************************************************************************************************
99 SUBROUTINE qs_vxc_create(ks_env, rho_struct, xc_section, vxc_rho, vxc_tau, exc, &
100 just_energy, edisp, dispersion_env, adiabatic_rescale_factor, &
101 pw_env_external, native_skala_atom_force)
102
103 TYPE(qs_ks_env_type), POINTER :: ks_env
104 TYPE(qs_rho_type), POINTER :: rho_struct
105 TYPE(section_vals_type), POINTER :: xc_section
106 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_rho, vxc_tau
107 REAL(kind=dp), INTENT(out) :: exc
108 LOGICAL, INTENT(in), OPTIONAL :: just_energy
109 REAL(kind=dp), INTENT(out), OPTIONAL :: edisp
110 TYPE(qs_dispersion_type), OPTIONAL, POINTER :: dispersion_env
111 REAL(kind=dp), INTENT(in), OPTIONAL :: adiabatic_rescale_factor
112 TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_external
113 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT), &
114 OPTIONAL :: native_skala_atom_force
115
116 CHARACTER(len=*), PARAMETER :: routinen = 'qs_vxc_create'
117
118 INTEGER :: handle, ispin, mspin, myfun, &
119 nelec_spin(2), vdw
120 LOGICAL :: compute_virial, do_adiabatic_rescaling, my_just_energy, native_skala_grid, &
121 rho_g_valid, sic_scaling_b_zero, tau_g_valid, tau_r_valid, uf_grid, vdw_nl
122 REAL(kind=dp) :: exc_m, factor, &
123 my_adiabatic_rescale_factor, &
124 my_scaling, nelec_s_inv
125 REAL(kind=dp), DIMENSION(3, 3) :: virial_xc_tmp
126 TYPE(cell_type), POINTER :: cell
127 TYPE(dft_control_type), POINTER :: dft_control
128 TYPE(mp_para_env_type), POINTER :: para_env
129 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
130 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g, rho_m_gspace, rho_struct_g, &
131 tau_struct_g
132 TYPE(pw_c1d_gs_type), POINTER :: rho_nlcc_g, rho_nlcc_g_use, &
133 rho_nlcc_g_xc, tmp_g, tmp_g2
134 TYPE(pw_env_type), POINTER :: pw_env
135 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool, vdw_pw_pool, xc_pw_pool
136 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: my_vxc_rho, my_vxc_tau, rho_m_rspace, &
137 rho_r, rho_struct_r, tau, tau_struct_r
138 TYPE(pw_r3d_rs_type), POINTER :: rho_nlcc, rho_nlcc_use, rho_nlcc_xc, &
139 tmp_pw, weights, weights_use, &
140 weights_xc
141 TYPE(virial_type), POINTER :: virial
142
143 CALL timeset(routinen, handle)
144
145 cpassert(.NOT. ASSOCIATED(vxc_rho))
146 cpassert(.NOT. ASSOCIATED(vxc_tau))
147 NULLIFY (dft_control, pw_env, auxbas_pw_pool, xc_pw_pool, vdw_pw_pool, cell, my_vxc_rho, &
148 tmp_pw, tmp_g, tmp_g2, my_vxc_tau, rho_g, rho_r, tau, rho_m_rspace, &
149 rho_m_gspace, rho_nlcc, rho_nlcc_g, rho_nlcc_g_use, rho_nlcc_g_xc, &
150 rho_nlcc_use, rho_nlcc_xc, rho_struct_r, rho_struct_g, tau_struct_g, tau_struct_r, &
151 weights_use, weights_xc, particle_set)
152
153 exc = 0.0_dp
154 my_just_energy = .false.
155 IF (PRESENT(just_energy)) my_just_energy = just_energy
156 my_adiabatic_rescale_factor = 1.0_dp
157 do_adiabatic_rescaling = .false.
158 IF (PRESENT(adiabatic_rescale_factor)) THEN
159 my_adiabatic_rescale_factor = adiabatic_rescale_factor
160 do_adiabatic_rescaling = .true.
161 END IF
162
163 CALL get_ks_env(ks_env, &
164 dft_control=dft_control, &
165 pw_env=pw_env, &
166 cell=cell, &
167 particle_set=particle_set, &
168 xcint_weights=weights, &
169 virial=virial, &
170 rho_nlcc=rho_nlcc, &
171 rho_nlcc_g=rho_nlcc_g)
172 rho_nlcc_use => rho_nlcc
173 rho_nlcc_g_use => rho_nlcc_g
174 weights_use => weights
175
176 CALL qs_rho_get(rho_struct, &
177 tau_r_valid=tau_r_valid, &
178 tau_g_valid=tau_g_valid, &
179 rho_g_valid=rho_g_valid, &
180 rho_r=rho_struct_r, &
181 rho_g=rho_struct_g, &
182 tau_g=tau_struct_g, &
183 tau_r=tau_struct_r)
184
185 compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
186 IF (compute_virial) THEN
187 virial%pv_xc = 0.0_dp
188 END IF
189
190 CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", &
191 i_val=myfun)
192 CALL section_vals_val_get(xc_section, "VDW_POTENTIAL%POTENTIAL_TYPE", &
193 i_val=vdw)
194
195 vdw_nl = (vdw == xc_vdw_fun_nonloc)
196 ! this combination has not been investigated
197 cpassert(.NOT. (do_adiabatic_rescaling .AND. vdw_nl))
198 ! are the necessary inputs available
199 IF (.NOT. (PRESENT(dispersion_env) .AND. PRESENT(edisp))) THEN
200 vdw_nl = .false.
201 END IF
202 IF (PRESENT(edisp)) edisp = 0.0_dp
203 native_skala_grid = xc_section_uses_native_skala_grid(xc_section)
204
205 IF (myfun /= xc_none .OR. vdw_nl) THEN
206
207 ! test if the real space density is available
208 cpassert(ASSOCIATED(rho_struct))
209 IF (dft_control%nspins /= 1 .AND. dft_control%nspins /= 2) THEN
210 cpabort("nspins must be 1 or 2")
211 END IF
212 mspin = SIZE(rho_struct_r)
213 IF (dft_control%nspins == 2 .AND. mspin == 1) THEN
214 cpabort("Spin count mismatch")
215 END IF
216
217 ! there are some options related to SIC here.
218 ! Normal DFT computes E(rho_alpha,rho_beta) (or its variant E(2*rho_alpha) for non-LSD)
219 ! SIC can E(rho_alpha,rho_beta)-b*(E(rho_alpha,rho_beta)-E(rho_beta,rho_beta))
220 ! or compute E(rho_alpha,rho_beta)-b*E(rho_alpha-rho_beta,0)
221
222 ! my_scaling is the scaling needed of the standard E(rho_alpha,rho_beta) term
223 my_scaling = 1.0_dp
224 SELECT CASE (dft_control%sic_method_id)
225 CASE (sic_none)
226 ! all fine
227 CASE (sic_mauri_spz, sic_ad)
228 ! no idea yet what to do here in that case
229 cpassert(.NOT. tau_r_valid)
230 CASE (sic_mauri_us)
231 my_scaling = 1.0_dp - dft_control%sic_scaling_b
232 ! no idea yet what to do here in that case
233 cpassert(.NOT. tau_r_valid)
234 CASE (sic_eo)
235 ! NOTHING TO BE DONE
236 CASE DEFAULT
237 ! this case has not yet been treated here
238 cpabort("NYI")
239 END SELECT
240
241 IF (dft_control%sic_scaling_b == 0.0_dp) THEN
242 sic_scaling_b_zero = .true.
243 ELSE
244 sic_scaling_b_zero = .false.
245 END IF
246
247 IF (PRESENT(pw_env_external)) THEN
248 pw_env => pw_env_external
249 END IF
250 CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool, auxbas_pw_pool=auxbas_pw_pool)
251 uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid, xc_pw_pool%pw_grid)
252
253 IF (.NOT. uf_grid) THEN
254 rho_r => rho_struct_r
255
256 IF (tau_r_valid) THEN
257 tau => tau_struct_r
258 END IF
259
260 ! for gradient corrected functional the density in g space might
261 ! be useful so if we have it, we pass it in
262 IF (rho_g_valid) THEN
263 rho_g => rho_struct_g
264 END IF
265 ELSE
266 cpassert(rho_g_valid)
267 ALLOCATE (rho_r(mspin))
268 ALLOCATE (rho_g(mspin))
269 DO ispin = 1, mspin
270 CALL xc_pw_pool%create_pw(rho_g(ispin))
271 CALL pw_transfer(rho_struct_g(ispin), rho_g(ispin))
272 END DO
273 DO ispin = 1, mspin
274 CALL xc_pw_pool%create_pw(rho_r(ispin))
275 CALL pw_transfer(rho_g(ispin), rho_r(ispin))
276 END DO
277 IF (tau_r_valid) THEN
278 ALLOCATE (tau(mspin))
279 DO ispin = 1, mspin
280 CALL xc_pw_pool%create_pw(tau(ispin))
281 block
282 TYPE(pw_c1d_gs_type) :: tau_g_aux, tau_g_xc
283 CALL xc_pw_pool%create_pw(tau_g_xc)
284 IF (tau_g_valid) THEN
285 CALL pw_transfer(tau_struct_g(ispin), tau_g_xc)
286 ELSE
287 CALL auxbas_pw_pool%create_pw(tau_g_aux)
288 CALL pw_transfer(tau_struct_r(ispin), tau_g_aux)
289 CALL pw_transfer(tau_g_aux, tau_g_xc)
290 CALL auxbas_pw_pool%give_back_pw(tau_g_aux)
291 END IF
292 CALL pw_transfer(tau_g_xc, tau(ispin))
293 CALL xc_pw_pool%give_back_pw(tau_g_xc)
294 END block
295 END DO
296 END IF
297 IF (ASSOCIATED(weights)) THEN
298 ALLOCATE (weights_xc)
299 CALL xc_pw_pool%create_pw(weights_xc)
300 block
301 TYPE(pw_c1d_gs_type) :: weights_g_aux, weights_g_xc
302 CALL auxbas_pw_pool%create_pw(weights_g_aux)
303 CALL xc_pw_pool%create_pw(weights_g_xc)
304 CALL pw_transfer(weights, weights_g_aux)
305 CALL pw_transfer(weights_g_aux, weights_g_xc)
306 CALL pw_transfer(weights_g_xc, weights_xc)
307 CALL xc_pw_pool%give_back_pw(weights_g_xc)
308 CALL auxbas_pw_pool%give_back_pw(weights_g_aux)
309 END block
310 weights_use => weights_xc
311 END IF
312 IF (ASSOCIATED(rho_nlcc)) THEN
313 cpassert(ASSOCIATED(rho_nlcc_g))
314 ALLOCATE (rho_nlcc_g_xc, rho_nlcc_xc)
315 CALL xc_pw_pool%create_pw(rho_nlcc_g_xc)
316 CALL xc_pw_pool%create_pw(rho_nlcc_xc)
317 CALL pw_transfer(rho_nlcc_g, rho_nlcc_g_xc)
318 CALL pw_transfer(rho_nlcc_g_xc, rho_nlcc_xc)
319 rho_nlcc_use => rho_nlcc_xc
320 rho_nlcc_g_use => rho_nlcc_g_xc
321 END IF
322 END IF
323
324 ! add the nlcc densities
325 IF (ASSOCIATED(rho_nlcc_use)) THEN
326 factor = 1.0_dp
327 DO ispin = 1, mspin
328 CALL pw_axpy(rho_nlcc_use, rho_r(ispin), factor)
329 CALL pw_axpy(rho_nlcc_g_use, rho_g(ispin), factor)
330 END DO
331 END IF
332
333 !
334 ! here the rho_r, rho_g, tau is what it should be
335 ! we get back the right my_vxc_rho and my_vxc_tau as required
336 !
337 IF (native_skala_grid) THEN
338 CALL skala_gpw_eval(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, exc=exc, &
339 rho_r=rho_r, rho_g=rho_g, tau=tau, xc_section=xc_section, &
340 weights=weights_use, pw_pool=xc_pw_pool, &
341 particle_set=particle_set, cell=cell, &
342 compute_virial=compute_virial, virial_xc=virial%pv_xc, &
343 just_energy=my_just_energy, atom_force=native_skala_atom_force)
344 ELSE IF (my_just_energy) THEN
345 exc = xc_exc_calc(rho_r=rho_r, tau=tau, &
346 rho_g=rho_g, xc_section=xc_section, &
347 weights=weights_use, pw_pool=xc_pw_pool)
348
349 ELSE
350 CALL xc_vxc_pw_create(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_r, &
351 rho_g=rho_g, tau=tau, exc=exc, &
352 xc_section=xc_section, &
353 weights=weights_use, pw_pool=xc_pw_pool, &
354 compute_virial=compute_virial, &
355 virial_xc=virial%pv_xc)
356 END IF
357
358 ! remove the nlcc densities (keep stuff in original state)
359 IF (ASSOCIATED(rho_nlcc_use)) THEN
360 factor = -1.0_dp
361 DO ispin = 1, mspin
362 CALL pw_axpy(rho_nlcc_use, rho_r(ispin), factor)
363 CALL pw_axpy(rho_nlcc_g_use, rho_g(ispin), factor)
364 END DO
365 END IF
366
367 ! calclulate non-local vdW functional
368 ! only if this XC_SECTION has it
369 ! if yes, we use the dispersion_env from ks_env
370 ! this is dangerous, as it assumes a special connection xc_section -> qs_env
371 IF (vdw_nl) THEN
372 CALL get_ks_env(ks_env=ks_env, para_env=para_env)
373 ! no SIC functionals allowed
374 cpassert(dft_control%sic_method_id == sic_none)
375 !
376 CALL pw_env_get(pw_env, vdw_pw_pool=vdw_pw_pool)
377 IF (my_just_energy) THEN
378 CALL calculate_dispersion_nonloc(my_vxc_rho, rho_r, rho_g, edisp, dispersion_env, &
379 my_just_energy, vdw_pw_pool, xc_pw_pool, para_env)
380 ELSE
381 CALL calculate_dispersion_nonloc(my_vxc_rho, rho_r, rho_g, edisp, dispersion_env, &
382 my_just_energy, vdw_pw_pool, xc_pw_pool, para_env, virial=virial)
383 END IF
384 END IF
385
386 !! Apply rescaling to the potential if requested
387 IF (.NOT. my_just_energy) THEN
388 IF (do_adiabatic_rescaling) THEN
389 IF (ASSOCIATED(my_vxc_rho)) THEN
390 DO ispin = 1, SIZE(my_vxc_rho)
391 CALL pw_scale(my_vxc_rho(ispin), my_adiabatic_rescale_factor)
392 END DO
393 END IF
394 END IF
395 END IF
396
397 IF (my_scaling /= 1.0_dp) THEN
398 exc = exc*my_scaling
399 IF (ASSOCIATED(my_vxc_rho)) THEN
400 DO ispin = 1, SIZE(my_vxc_rho)
401 CALL pw_scale(my_vxc_rho(ispin), my_scaling)
402 END DO
403 END IF
404 IF (ASSOCIATED(my_vxc_tau)) THEN
405 DO ispin = 1, SIZE(my_vxc_tau)
406 CALL pw_scale(my_vxc_tau(ispin), my_scaling)
407 END DO
408 END IF
409 END IF
410
411 ! we have pw data for the xc, qs_ks requests coeff structure, here we transfer
412 ! pw -> coeff
413 IF (ASSOCIATED(my_vxc_rho)) THEN
414 vxc_rho => my_vxc_rho
415 NULLIFY (my_vxc_rho)
416 END IF
417 IF (ASSOCIATED(my_vxc_tau)) THEN
418 vxc_tau => my_vxc_tau
419 NULLIFY (my_vxc_tau)
420 END IF
421 IF (uf_grid) THEN
422 DO ispin = 1, SIZE(rho_r)
423 CALL xc_pw_pool%give_back_pw(rho_r(ispin))
424 END DO
425 DEALLOCATE (rho_r)
426 IF (ASSOCIATED(rho_g)) THEN
427 DO ispin = 1, SIZE(rho_g)
428 CALL xc_pw_pool%give_back_pw(rho_g(ispin))
429 END DO
430 DEALLOCATE (rho_g)
431 END IF
432 END IF
433
434 ! compute again the xc but now for Exc(m,o) and the opposite sign
435 IF (dft_control%sic_method_id == sic_mauri_spz .AND. .NOT. sic_scaling_b_zero) THEN
436 ALLOCATE (rho_m_rspace(2), rho_m_gspace(2))
437 CALL xc_pw_pool%create_pw(rho_m_gspace(1))
438 CALL xc_pw_pool%create_pw(rho_m_rspace(1))
439 CALL pw_copy(rho_struct_r(1), rho_m_rspace(1))
440 CALL pw_axpy(rho_struct_r(2), rho_m_rspace(1), alpha=-1._dp)
441 CALL pw_copy(rho_struct_g(1), rho_m_gspace(1))
442 CALL pw_axpy(rho_struct_g(2), rho_m_gspace(1), alpha=-1._dp)
443 ! bit sad, these will be just zero...
444 CALL xc_pw_pool%create_pw(rho_m_gspace(2))
445 CALL xc_pw_pool%create_pw(rho_m_rspace(2))
446 CALL pw_zero(rho_m_rspace(2))
447 CALL pw_zero(rho_m_gspace(2))
448
449 IF (my_just_energy) THEN
450 exc_m = xc_exc_calc(rho_r=rho_m_rspace, tau=tau, &
451 rho_g=rho_m_gspace, xc_section=xc_section, &
452 weights=weights_use, pw_pool=xc_pw_pool)
453 ELSE
454 ! virial untested
455 cpassert(.NOT. compute_virial)
456 CALL xc_vxc_pw_create(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_m_rspace, &
457 rho_g=rho_m_gspace, tau=tau, exc=exc_m, &
458 xc_section=xc_section, &
459 weights=weights_use, pw_pool=xc_pw_pool, &
460 compute_virial=.false., &
461 virial_xc=virial_xc_tmp)
462 END IF
463
464 exc = exc - dft_control%sic_scaling_b*exc_m
465
466 ! and take care of the potential only vxc_rho is taken into account
467 IF (.NOT. my_just_energy) THEN
468 CALL pw_axpy(my_vxc_rho(1), vxc_rho(1), -dft_control%sic_scaling_b)
469 CALL pw_axpy(my_vxc_rho(1), vxc_rho(2), dft_control%sic_scaling_b)
470 CALL my_vxc_rho(1)%release()
471 CALL my_vxc_rho(2)%release()
472 DEALLOCATE (my_vxc_rho)
473 END IF
474
475 DO ispin = 1, 2
476 CALL xc_pw_pool%give_back_pw(rho_m_rspace(ispin))
477 CALL xc_pw_pool%give_back_pw(rho_m_gspace(ispin))
478 END DO
479 DEALLOCATE (rho_m_rspace)
480 DEALLOCATE (rho_m_gspace)
481
482 END IF
483
484 ! now we have - sum_s N_s * Exc(rho_s/N_s,0)
485 IF (dft_control%sic_method_id == sic_ad .AND. .NOT. sic_scaling_b_zero) THEN
486
487 ! find out how many elecs we have
488 CALL get_ks_env(ks_env, nelectron_spin=nelec_spin)
489
490 ALLOCATE (rho_m_rspace(2), rho_m_gspace(2))
491 DO ispin = 1, 2
492 CALL xc_pw_pool%create_pw(rho_m_gspace(ispin))
493 CALL xc_pw_pool%create_pw(rho_m_rspace(ispin))
494 END DO
495
496 DO ispin = 1, 2
497 IF (nelec_spin(ispin) > 0.0_dp) THEN
498 nelec_s_inv = 1.0_dp/nelec_spin(ispin)
499 ELSE
500 ! does it matter if there are no electrons with this spin (H) ?
501 nelec_s_inv = 0.0_dp
502 END IF
503 CALL pw_copy(rho_struct_r(ispin), rho_m_rspace(1))
504 CALL pw_copy(rho_struct_g(ispin), rho_m_gspace(1))
505 CALL pw_scale(rho_m_rspace(1), nelec_s_inv)
506 CALL pw_scale(rho_m_gspace(1), nelec_s_inv)
507 CALL pw_zero(rho_m_rspace(2))
508 CALL pw_zero(rho_m_gspace(2))
509
510 IF (my_just_energy) THEN
511 exc_m = xc_exc_calc(rho_r=rho_m_rspace, tau=tau, &
512 rho_g=rho_m_gspace, xc_section=xc_section, &
513 weights=weights_use, pw_pool=xc_pw_pool)
514 ELSE
515 ! virial untested
516 cpassert(.NOT. compute_virial)
517 CALL xc_vxc_pw_create(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_m_rspace, &
518 rho_g=rho_m_gspace, tau=tau, exc=exc_m, &
519 xc_section=xc_section, &
520 weights=weights_use, pw_pool=xc_pw_pool, &
521 compute_virial=.false., &
522 virial_xc=virial_xc_tmp)
523 END IF
524
525 exc = exc - dft_control%sic_scaling_b*nelec_spin(ispin)*exc_m
526
527 ! and take care of the potential only vxc_rho is taken into account
528 IF (.NOT. my_just_energy) THEN
529 CALL pw_axpy(my_vxc_rho(1), vxc_rho(ispin), -dft_control%sic_scaling_b)
530 CALL my_vxc_rho(1)%release()
531 CALL my_vxc_rho(2)%release()
532 DEALLOCATE (my_vxc_rho)
533 END IF
534 END DO
535
536 DO ispin = 1, 2
537 CALL xc_pw_pool%give_back_pw(rho_m_rspace(ispin))
538 CALL xc_pw_pool%give_back_pw(rho_m_gspace(ispin))
539 END DO
540 DEALLOCATE (rho_m_rspace)
541 DEALLOCATE (rho_m_gspace)
542
543 END IF
544
545 ! compute again the xc but now for Exc(n_down,n_down)
546 IF (dft_control%sic_method_id == sic_mauri_us .AND. .NOT. sic_scaling_b_zero) THEN
547 ALLOCATE (rho_r(2))
548 rho_r(1) = rho_struct_r(2)
549 rho_r(2) = rho_struct_r(2)
550 IF (rho_g_valid) THEN
551 ALLOCATE (rho_g(2))
552 rho_g(1) = rho_struct_g(2)
553 rho_g(2) = rho_struct_g(2)
554 END IF
555
556 IF (my_just_energy) THEN
557 exc_m = xc_exc_calc(rho_r=rho_r, tau=tau, &
558 rho_g=rho_g, xc_section=xc_section, &
559 weights=weights_use, pw_pool=xc_pw_pool)
560 ELSE
561 ! virial untested
562 cpassert(.NOT. compute_virial)
563 CALL xc_vxc_pw_create(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_r, &
564 rho_g=rho_g, tau=tau, exc=exc_m, &
565 xc_section=xc_section, &
566 weights=weights_use, pw_pool=xc_pw_pool, &
567 compute_virial=.false., &
568 virial_xc=virial_xc_tmp)
569 END IF
570
571 exc = exc + dft_control%sic_scaling_b*exc_m
572
573 ! and take care of the potential
574 IF (.NOT. my_just_energy) THEN
575 ! both go to minority spin
576 CALL pw_axpy(my_vxc_rho(1), vxc_rho(2), 2.0_dp*dft_control%sic_scaling_b)
577 CALL my_vxc_rho(1)%release()
578 CALL my_vxc_rho(2)%release()
579 DEALLOCATE (my_vxc_rho)
580 END IF
581 DEALLOCATE (rho_r, rho_g)
582
583 END IF
584
585 !
586 ! cleanups
587 !
588 IF (uf_grid .AND. (ASSOCIATED(vxc_rho) .OR. ASSOCIATED(vxc_tau))) THEN
589 block
590 TYPE(pw_r3d_rs_type) :: tmp_pw
591 TYPE(pw_c1d_gs_type) :: tmp_g, tmp_g2
592 CALL xc_pw_pool%create_pw(tmp_g)
593 CALL auxbas_pw_pool%create_pw(tmp_g2)
594 IF (ASSOCIATED(vxc_rho)) THEN
595 DO ispin = 1, SIZE(vxc_rho)
596 CALL auxbas_pw_pool%create_pw(tmp_pw)
597 CALL pw_transfer(vxc_rho(ispin), tmp_g)
598 CALL pw_transfer(tmp_g, tmp_g2)
599 CALL pw_transfer(tmp_g2, tmp_pw)
600 CALL xc_pw_pool%give_back_pw(vxc_rho(ispin))
601 vxc_rho(ispin) = tmp_pw
602 END DO
603 END IF
604 IF (ASSOCIATED(vxc_tau)) THEN
605 DO ispin = 1, SIZE(vxc_tau)
606 CALL auxbas_pw_pool%create_pw(tmp_pw)
607 CALL pw_transfer(vxc_tau(ispin), tmp_g)
608 CALL pw_transfer(tmp_g, tmp_g2)
609 CALL pw_transfer(tmp_g2, tmp_pw)
610 CALL xc_pw_pool%give_back_pw(vxc_tau(ispin))
611 vxc_tau(ispin) = tmp_pw
612 END DO
613 END IF
614 CALL auxbas_pw_pool%give_back_pw(tmp_g2)
615 CALL xc_pw_pool%give_back_pw(tmp_g)
616 END block
617 END IF
618 IF (ASSOCIATED(tau) .AND. uf_grid) THEN
619 DO ispin = 1, SIZE(tau)
620 CALL xc_pw_pool%give_back_pw(tau(ispin))
621 END DO
622 DEALLOCATE (tau)
623 END IF
624 IF (ASSOCIATED(weights_xc)) THEN
625 CALL xc_pw_pool%give_back_pw(weights_xc)
626 DEALLOCATE (weights_xc)
627 END IF
628 IF (ASSOCIATED(rho_nlcc_xc)) THEN
629 CALL xc_pw_pool%give_back_pw(rho_nlcc_xc)
630 DEALLOCATE (rho_nlcc_xc)
631 END IF
632 IF (ASSOCIATED(rho_nlcc_g_xc)) THEN
633 CALL xc_pw_pool%give_back_pw(rho_nlcc_g_xc)
634 DEALLOCATE (rho_nlcc_g_xc)
635 END IF
636
637 END IF
638
639 CALL timestop(handle)
640
641 END SUBROUTINE qs_vxc_create
642
643! **************************************************************************************************
644!> \brief calculates the XC density: E_xc(r) - V_xc(r)*rho(r) or E_xc(r)/rho(r)
645!> \param ks_env to get all the needed things
646!> \param rho_struct density
647!> \param xc_section ...
648!> \param dispersion_env ...
649!> \param xc_ener will contain the xc energy density E_xc(r) - V_xc(r)*rho(r)
650!> \param xc_den will contain the xc energy density E_xc(r)/rho(r)
651!> \param exc will contain the xc energy density E_xc(r)
652!> \param vxc ...
653!> \param vtau ...
654!> \author JGH
655! **************************************************************************************************
656 SUBROUTINE qs_xc_density(ks_env, rho_struct, xc_section, dispersion_env, &
657 xc_ener, xc_den, exc, vxc, vtau)
658
659 TYPE(qs_ks_env_type), POINTER :: ks_env
660 TYPE(qs_rho_type), POINTER :: rho_struct
661 TYPE(section_vals_type), POINTER :: xc_section
662 TYPE(qs_dispersion_type), OPTIONAL, POINTER :: dispersion_env
663 TYPE(pw_r3d_rs_type), INTENT(INOUT), OPTIONAL :: xc_ener, xc_den
664 TYPE(pw_r3d_rs_type), OPTIONAL :: exc
665 TYPE(pw_r3d_rs_type), DIMENSION(:), OPTIONAL :: vxc, vtau
666
667 CHARACTER(len=*), PARAMETER :: routinen = 'qs_xc_density'
668
669 INTEGER :: handle, ispin, mspin, myfun, nspins, vdw
670 LOGICAL :: rho_g_valid, tau_g_valid, tau_r_valid, &
671 uf_grid, vdw_nl
672 REAL(kind=dp) :: edisp, excint, factor, rho_cutoff
673 REAL(kind=dp), DIMENSION(3, 3) :: vdum
674 TYPE(cell_type), POINTER :: cell
675 TYPE(dft_control_type), POINTER :: dft_control
676 TYPE(mp_para_env_type), POINTER :: para_env
677 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g, rho_struct_g, tau_g, tau_struct_g
678 TYPE(pw_c1d_gs_type), POINTER :: rho_nlcc_g, rho_nlcc_g_use, rho_nlcc_g_xc
679 TYPE(pw_env_type), POINTER :: pw_env
680 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool, vdw_pw_pool, xc_pw_pool
681 TYPE(pw_r3d_rs_type) :: exc_r
682 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, rho_struct_r, tau_r, &
683 tau_struct_r, vxc_rho, vxc_tau
684 TYPE(pw_r3d_rs_type), POINTER :: rho_nlcc, rho_nlcc_use, rho_nlcc_xc, &
685 weights, weights_use, weights_xc
686
687 CALL timeset(routinen, handle)
688
689 NULLIFY (dft_control, pw_env, auxbas_pw_pool, xc_pw_pool, vdw_pw_pool, cell, &
690 rho_g, rho_struct_g, tau_g, tau_struct_g, rho_nlcc, rho_nlcc_g, &
691 rho_nlcc_g_use, rho_nlcc_g_xc, rho_nlcc_use, rho_nlcc_xc, rho_r, &
692 rho_struct_r, tau_r, tau_struct_r, vxc_rho, vxc_tau, weights, &
693 weights_use, weights_xc)
694
695 CALL get_ks_env(ks_env, &
696 dft_control=dft_control, &
697 pw_env=pw_env, &
698 cell=cell, &
699 xcint_weights=weights, &
700 rho_nlcc=rho_nlcc, &
701 rho_nlcc_g=rho_nlcc_g)
702
703 CALL qs_rho_get(rho_struct, &
704 tau_r_valid=tau_r_valid, &
705 tau_g_valid=tau_g_valid, &
706 rho_g_valid=rho_g_valid, &
707 rho_r=rho_struct_r, &
708 rho_g=rho_struct_g, &
709 tau_r=tau_struct_r, &
710 tau_g=tau_struct_g)
711 nspins = dft_control%nspins
712 mspin = SIZE(rho_struct_r)
713 rho_r => rho_struct_r
714 rho_g => rho_struct_g
715 tau_r => tau_struct_r
716 tau_g => tau_struct_g
717 rho_nlcc_use => rho_nlcc
718 rho_nlcc_g_use => rho_nlcc_g
719 weights_use => weights
720
721 CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", i_val=myfun)
722 CALL section_vals_val_get(xc_section, "VDW_POTENTIAL%POTENTIAL_TYPE", i_val=vdw)
723 vdw_nl = (vdw == xc_vdw_fun_nonloc)
724 IF (PRESENT(xc_ener)) THEN
725 IF (tau_r_valid) THEN
726 CALL cp_warn(__location__, "Tau contribution will not be correctly handled")
727 END IF
728 END IF
729 IF (vdw_nl) THEN
730 CALL cp_warn(__location__, "vdW functional contribution will be ignored")
731 END IF
732
733 CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool, auxbas_pw_pool=auxbas_pw_pool)
734 uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid, xc_pw_pool%pw_grid)
735
736 IF (PRESENT(xc_ener)) THEN
737 CALL pw_zero(xc_ener)
738 END IF
739 IF (PRESENT(xc_den)) THEN
740 CALL pw_zero(xc_den)
741 END IF
742 IF (PRESENT(exc)) THEN
743 CALL pw_zero(exc)
744 END IF
745 IF (PRESENT(vxc)) THEN
746 DO ispin = 1, nspins
747 CALL pw_zero(vxc(ispin))
748 END DO
749 END IF
750 IF (PRESENT(vtau)) THEN
751 DO ispin = 1, nspins
752 CALL pw_zero(vtau(ispin))
753 END DO
754 END IF
755
756 IF (myfun /= xc_none) THEN
757
758 cpassert(ASSOCIATED(rho_struct))
759 cpassert(dft_control%sic_method_id == sic_none)
760
761 IF (uf_grid) THEN
762 NULLIFY (rho_r, rho_g, tau_r, tau_g)
763 IF (rho_g_valid) THEN
764 CALL create_density_on_pool(xc_pw_pool, rho_struct_g, rho_r, rho_g)
765 ELSE IF (ASSOCIATED(rho_struct_r)) THEN
766 CALL create_density_on_pool_from_r(auxbas_pw_pool, xc_pw_pool, rho_struct_r, rho_r, rho_g)
767 ELSE
768 cpabort("Fine Grid in qs_xc_density requires rho_r or rho_g")
769 END IF
770 IF (tau_r_valid) THEN
771 IF (tau_g_valid) THEN
772 CALL create_density_on_pool(xc_pw_pool, tau_struct_g, tau_r, tau_g)
773 ELSE IF (ASSOCIATED(tau_struct_r)) THEN
774 CALL create_density_on_pool_from_r(auxbas_pw_pool, xc_pw_pool, tau_struct_r, tau_r, tau_g)
775 ELSE
776 cpabort("Fine Grid in qs_xc_density requires tau_r or tau_g")
777 END IF
778 END IF
779 IF (ASSOCIATED(weights)) THEN
780 ALLOCATE (weights_xc)
781 CALL xc_pw_pool%create_pw(weights_xc)
782 CALL transfer_rspace_between_pools(auxbas_pw_pool, xc_pw_pool, weights, weights_xc)
783 weights_use => weights_xc
784 END IF
785 IF (ASSOCIATED(rho_nlcc)) THEN
786 cpassert(ASSOCIATED(rho_nlcc_g))
787 ALLOCATE (rho_nlcc_g_xc, rho_nlcc_xc)
788 CALL xc_pw_pool%create_pw(rho_nlcc_g_xc)
789 CALL xc_pw_pool%create_pw(rho_nlcc_xc)
790 CALL pw_transfer(rho_nlcc_g, rho_nlcc_g_xc)
791 CALL pw_transfer(rho_nlcc_g_xc, rho_nlcc_xc)
792 rho_nlcc_use => rho_nlcc_xc
793 rho_nlcc_g_use => rho_nlcc_g_xc
794 END IF
795 END IF
796
797 ! add the nlcc densities
798 IF (ASSOCIATED(rho_nlcc_use)) THEN
799 factor = 1.0_dp
800 DO ispin = 1, mspin
801 CALL pw_axpy(rho_nlcc_use, rho_r(ispin), factor)
802 CALL pw_axpy(rho_nlcc_g_use, rho_g(ispin), factor)
803 END DO
804 END IF
805 NULLIFY (vxc_rho, vxc_tau)
806 CALL xc_vxc_pw_create(vxc_rho=vxc_rho, vxc_tau=vxc_tau, rho_r=rho_r, &
807 rho_g=rho_g, tau=tau_r, exc=excint, &
808 xc_section=xc_section, &
809 weights=weights_use, pw_pool=xc_pw_pool, &
810 compute_virial=.false., &
811 virial_xc=vdum, &
812 exc_r=exc_r)
813 ! calclulate non-local vdW functional
814 ! only if this XC_SECTION has it
815 ! if yes, we use the dispersion_env from ks_env
816 ! this is dangerous, as it assumes a special connection xc_section -> qs_env
817 IF (vdw_nl) THEN
818 CALL get_ks_env(ks_env=ks_env, para_env=para_env)
819 ! no SIC functionals allowed
820 cpassert(dft_control%sic_method_id == sic_none)
821 !
822 CALL pw_env_get(pw_env, vdw_pw_pool=vdw_pw_pool)
823 CALL calculate_dispersion_nonloc(vxc_rho, rho_r, rho_g, edisp, dispersion_env, &
824 .false., vdw_pw_pool, xc_pw_pool, para_env)
825 END IF
826
827 ! remove the nlcc densities (keep stuff in original state)
828 IF (ASSOCIATED(rho_nlcc_use)) THEN
829 factor = -1.0_dp
830 DO ispin = 1, mspin
831 CALL pw_axpy(rho_nlcc_use, rho_r(ispin), factor)
832 CALL pw_axpy(rho_nlcc_g_use, rho_g(ispin), factor)
833 END DO
834 END IF
835 !
836 IF (PRESENT(xc_den)) THEN
837 rho_cutoff = 1.e-14_dp
838 IF (uf_grid) THEN
839 block
840 TYPE(pw_r3d_rs_type) :: tmp_pw
841 CALL xc_pw_pool%create_pw(tmp_pw)
842 CALL pw_copy(exc_r, tmp_pw)
843 CALL calc_xc_density(tmp_pw, rho_r, rho_cutoff)
844 CALL transfer_rspace_between_pools(xc_pw_pool, auxbas_pw_pool, tmp_pw, xc_den)
845 CALL xc_pw_pool%give_back_pw(tmp_pw)
846 END block
847 ELSE
848 CALL pw_copy(exc_r, xc_den)
849 CALL calc_xc_density(xc_den, rho_r, rho_cutoff)
850 END IF
851 END IF
852 IF (PRESENT(xc_ener)) THEN
853 IF (uf_grid) THEN
854 block
855 TYPE(pw_r3d_rs_type) :: tmp_pw
856 CALL xc_pw_pool%create_pw(tmp_pw)
857 CALL pw_copy(exc_r, tmp_pw)
858 DO ispin = 1, nspins
859 CALL pw_multiply(tmp_pw, vxc_rho(ispin), rho_r(ispin), alpha=-1.0_dp)
860 END DO
861 CALL transfer_rspace_between_pools(xc_pw_pool, auxbas_pw_pool, tmp_pw, xc_ener)
862 CALL xc_pw_pool%give_back_pw(tmp_pw)
863 END block
864 ELSE
865 CALL pw_copy(exc_r, xc_ener)
866 DO ispin = 1, nspins
867 CALL pw_multiply(xc_ener, vxc_rho(ispin), rho_r(ispin), alpha=-1.0_dp)
868 END DO
869 END IF
870 END IF
871 IF (PRESENT(exc)) THEN
872 IF (uf_grid) THEN
873 CALL transfer_rspace_between_pools(xc_pw_pool, auxbas_pw_pool, exc_r, exc)
874 ELSE
875 CALL pw_copy(exc_r, exc)
876 END IF
877 END IF
878 IF (PRESENT(vxc)) THEN
879 DO ispin = 1, nspins
880 IF (uf_grid) THEN
881 CALL transfer_rspace_between_pools(xc_pw_pool, auxbas_pw_pool, vxc_rho(ispin), vxc(ispin))
882 ELSE
883 CALL pw_copy(vxc_rho(ispin), vxc(ispin))
884 END IF
885 END DO
886 END IF
887 IF (PRESENT(vtau) .AND. ASSOCIATED(vxc_tau)) THEN
888 DO ispin = 1, nspins
889 IF (uf_grid) THEN
890 CALL transfer_rspace_between_pools(xc_pw_pool, auxbas_pw_pool, vxc_tau(ispin), vtau(ispin))
891 ELSE
892 CALL pw_copy(vxc_tau(ispin), vtau(ispin))
893 END IF
894 END DO
895 END IF
896 ! remove arrays
897 IF (ASSOCIATED(vxc_rho)) THEN
898 DO ispin = 1, nspins
899 CALL vxc_rho(ispin)%release()
900 END DO
901 DEALLOCATE (vxc_rho)
902 END IF
903 IF (ASSOCIATED(vxc_tau)) THEN
904 DO ispin = 1, nspins
905 CALL vxc_tau(ispin)%release()
906 END DO
907 DEALLOCATE (vxc_tau)
908 END IF
909 CALL exc_r%release()
910 IF (uf_grid) THEN
911 CALL give_back_density_on_pool(xc_pw_pool, rho_r, rho_g)
912 IF (ASSOCIATED(tau_r)) CALL give_back_density_on_pool(xc_pw_pool, tau_r, tau_g)
913 IF (ASSOCIATED(weights_xc)) THEN
914 CALL xc_pw_pool%give_back_pw(weights_xc)
915 DEALLOCATE (weights_xc)
916 END IF
917 IF (ASSOCIATED(rho_nlcc_xc)) THEN
918 CALL xc_pw_pool%give_back_pw(rho_nlcc_xc)
919 DEALLOCATE (rho_nlcc_xc)
920 END IF
921 IF (ASSOCIATED(rho_nlcc_g_xc)) THEN
922 CALL xc_pw_pool%give_back_pw(rho_nlcc_g_xc)
923 DEALLOCATE (rho_nlcc_g_xc)
924 END IF
925 END IF
926 !
927 END IF
928
929 CALL timestop(handle)
930
931 END SUBROUTINE qs_xc_density
932
933! **************************************************************************************************
934!> \brief transfers an r-space PW between two pools and writes into an existing target PW
935!> \param source_pw_pool ...
936!> \param target_pw_pool ...
937!> \param source ...
938!> \param TARGET ...
939! **************************************************************************************************
940 SUBROUTINE transfer_rspace_between_pools(source_pw_pool, target_pw_pool, source, TARGET)
941 TYPE(pw_pool_type), POINTER :: source_pw_pool, target_pw_pool
942 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: source, target
943
944 TYPE(pw_c1d_gs_type) :: source_g, target_g
945
946 cpassert(ASSOCIATED(source_pw_pool))
947 cpassert(ASSOCIATED(target_pw_pool))
948
949 IF (pw_grid_compare(source_pw_pool%pw_grid, target_pw_pool%pw_grid)) THEN
950 CALL pw_copy(source, TARGET)
951 ELSE
952 CALL source_pw_pool%create_pw(source_g)
953 CALL target_pw_pool%create_pw(target_g)
954 CALL pw_transfer(source, source_g)
955 CALL pw_transfer(source_g, target_g)
956 CALL pw_transfer(target_g, TARGET)
957 CALL target_pw_pool%give_back_pw(target_g)
958 CALL source_pw_pool%give_back_pw(source_g)
959 END IF
960
961 END SUBROUTINE transfer_rspace_between_pools
962
963! **************************************************************************************************
964!> \brief transfers a g-space density to a given PW pool and creates its r-space representation
965!> \param pw_pool ...
966!> \param rho_g_in ...
967!> \param rho_r_out ...
968!> \param rho_g_out ...
969! **************************************************************************************************
970 SUBROUTINE create_density_on_pool(pw_pool, rho_g_in, rho_r_out, rho_g_out)
971 TYPE(pw_pool_type), POINTER :: pw_pool
972 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_in
973 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_out
974 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_out
975
976 INTEGER :: ispin, nspins
977
978 cpassert(ASSOCIATED(pw_pool))
979 cpassert(ASSOCIATED(rho_g_in))
980
981 nspins = SIZE(rho_g_in)
982 ALLOCATE (rho_r_out(nspins), rho_g_out(nspins))
983 DO ispin = 1, nspins
984 CALL pw_pool%create_pw(rho_g_out(ispin))
985 CALL pw_pool%create_pw(rho_r_out(ispin))
986 CALL pw_transfer(rho_g_in(ispin), rho_g_out(ispin))
987 CALL pw_transfer(rho_g_out(ispin), rho_r_out(ispin))
988 END DO
989
990 END SUBROUTINE create_density_on_pool
991
992! **************************************************************************************************
993!> \brief transfers an r-space density to a given PW pool and creates its g-space representation
994!> \param source_pw_pool ...
995!> \param target_pw_pool ...
996!> \param rho_r_in ...
997!> \param rho_r_out ...
998!> \param rho_g_out ...
999! **************************************************************************************************
1000 SUBROUTINE create_density_on_pool_from_r(source_pw_pool, target_pw_pool, rho_r_in, rho_r_out, rho_g_out)
1001 TYPE(pw_pool_type), POINTER :: source_pw_pool, target_pw_pool
1002 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_in, rho_r_out
1003 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_out
1004
1005 INTEGER :: ispin, nspins
1006 TYPE(pw_c1d_gs_type) :: rho_g_in
1007
1008 cpassert(ASSOCIATED(source_pw_pool))
1009 cpassert(ASSOCIATED(target_pw_pool))
1010 cpassert(ASSOCIATED(rho_r_in))
1011
1012 nspins = SIZE(rho_r_in)
1013 ALLOCATE (rho_r_out(nspins), rho_g_out(nspins))
1014 DO ispin = 1, nspins
1015 CALL source_pw_pool%create_pw(rho_g_in)
1016 CALL target_pw_pool%create_pw(rho_g_out(ispin))
1017 CALL target_pw_pool%create_pw(rho_r_out(ispin))
1018 CALL pw_transfer(rho_r_in(ispin), rho_g_in)
1019 CALL pw_transfer(rho_g_in, rho_g_out(ispin))
1020 CALL pw_transfer(rho_g_out(ispin), rho_r_out(ispin))
1021 CALL source_pw_pool%give_back_pw(rho_g_in)
1022 END DO
1023
1024 END SUBROUTINE create_density_on_pool_from_r
1025
1026! **************************************************************************************************
1027!> \brief returns temporary density arrays to the given PW pool
1028!> \param pw_pool ...
1029!> \param rho_r ...
1030!> \param rho_g ...
1031! **************************************************************************************************
1032 SUBROUTINE give_back_density_on_pool(pw_pool, rho_r, rho_g)
1033 TYPE(pw_pool_type), POINTER :: pw_pool
1034 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
1035 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
1036
1037 INTEGER :: ispin
1038
1039 cpassert(ASSOCIATED(pw_pool))
1040
1041 IF (ASSOCIATED(rho_r)) THEN
1042 DO ispin = 1, SIZE(rho_r)
1043 CALL pw_pool%give_back_pw(rho_r(ispin))
1044 END DO
1045 DEALLOCATE (rho_r)
1046 END IF
1047 IF (ASSOCIATED(rho_g)) THEN
1048 DO ispin = 1, SIZE(rho_g)
1049 CALL pw_pool%give_back_pw(rho_g(ispin))
1050 END DO
1051 DEALLOCATE (rho_g)
1052 END IF
1053
1054 END SUBROUTINE give_back_density_on_pool
1055
1056END MODULE qs_vxc
Handles all functions related to the CELL.
Definition cell_types.F:15
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public sic_mauri_spz
integer, parameter, public xc_vdw_fun_nonloc
integer, parameter, public sic_eo
integer, parameter, public sic_mauri_us
integer, parameter, public sic_none
integer, parameter, public xc_none
integer, parameter, public sic_ad
objects that represent the structure of input sections and the data contained in an input section
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Interface to the message passing library MPI.
Define the data structure for the particle information.
container for various plainwaves related things
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
This module defines the grid data type and some basic operations on it.
Definition pw_grids.F:36
logical function, public pw_grid_compare(grida, gridb)
Check if two pw_grids are equal.
Definition pw_grids.F:148
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Calculation of non local dispersion functionals Some routines adapted from: Copyright (C) 2001-2009 Q...
subroutine, public calculate_dispersion_nonloc(vxc_rho, rho_r, rho_g, edispersion, dispersion_env, energy_only, pw_pool, xc_pw_pool, para_env, virial)
Calculates the non-local vdW functional using the method of Soler For spin polarized cases we use E(a...
Definition of disperson types for DFT calculations.
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, exc_accint, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, rho, rho_xc, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, task_list, task_list_soft, kpoints, do_kpoints, atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, particle_set, energy, force, local_particles, local_molecules, molecule_kind_set, molecule_set, subsys, cp_subsys, virial, results, atprop, nkind, natom, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env, nelectron_total, nelectron_spin)
...
superstucture that hold various representations of the density and keeps track of which ones are vali...
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
subroutine, public qs_vxc_create(ks_env, rho_struct, xc_section, vxc_rho, vxc_tau, exc, just_energy, edisp, dispersion_env, adiabatic_rescale_factor, pw_env_external, native_skala_atom_force)
calculates and allocates the xc potential, already reducing it to the dependence on rho and the one o...
Definition qs_vxc.F:102
subroutine, public qs_xc_density(ks_env, rho_struct, xc_section, dispersion_env, xc_ener, xc_den, exc, vxc, vtau)
calculates the XC density: E_xc(r) - V_xc(r)*rho(r) or E_xc(r)/rho(r)
Definition qs_vxc.F:658
Experimental CP2K-native GPW real-space-grid path for SKALA TorchScript models.
subroutine, public skala_gpw_eval(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau, xc_section, weights, pw_pool, particle_set, cell, compute_virial, virial_xc, just_energy, atom_force)
Evaluate SKALA energy and first derivatives on a CP2K GPW grid.
logical function, public xc_section_uses_native_skala_grid(xc_section)
Return true if the GAUXC subsection requests the CP2K-native GPW grid path.
Exchange and Correlation functional calculations.
Definition xc.F:17
real(kind=dp) function, public xc_exc_calc(rho_r, rho_g, tau, xc_section, weights, pw_pool)
calculates just the exchange and correlation energy (no vxc)
Definition xc.F:791
subroutine, public xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau, xc_section, weights, pw_pool, compute_virial, virial_xc, exc_r)
Exchange and Correlation functional calculations.
Definition xc.F:474
subroutine, public calc_xc_density(pot, rho, rho_cutoff)
Definition xc.F:402
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
stores all the informations relevant to an mpi environment
contained for different pw related things
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
keeps the density in various representations, keeping track of which ones are valid.