(git:9ea9339)
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
19 USE cp_control_types, ONLY: dft_control_type
20 USE input_constants, ONLY: sic_ad,&
21 sic_eo,&
22 sic_mauri_spz,&
23 sic_mauri_us,&
24 sic_none,&
25 xc_none,&
26 xc_vdw_fun_nonloc
27 USE input_section_types, ONLY: section_vals_type,&
28 section_vals_val_get
29 USE kinds, ONLY: dp
30 USE message_passing, ONLY: mp_para_env_type
31 USE pw_env_types, ONLY: pw_env_get,&
32 pw_env_type
33 USE pw_grids, ONLY: pw_grid_compare
34 USE pw_methods, ONLY: pw_axpy,&
35 pw_copy,&
36 pw_multiply,&
37 pw_scale,&
38 pw_transfer,&
39 pw_zero
40 USE pw_pool_types, ONLY: pw_pool_type
41 USE pw_types, ONLY: pw_c1d_gs_type,&
42 pw_r3d_rs_type
43 USE qs_dispersion_nonloc, ONLY: calculate_dispersion_nonloc
44 USE qs_dispersion_types, ONLY: qs_dispersion_type
45 USE qs_ks_types, ONLY: get_ks_env,&
46 qs_ks_env_type
47 USE qs_rho_types, ONLY: qs_rho_get,&
48 qs_rho_type
49 USE virial_types, ONLY: virial_type
50 USE xc, ONLY: calc_xc_density,&
51 xc_exc_calc,&
52 xc_vxc_pw_create
53#include "./base/base_uses.f90"
54
55 IMPLICIT NONE
56
57 PRIVATE
58
59 ! *** Public subroutines ***
60 PUBLIC :: qs_vxc_create, qs_xc_density
61
62 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_vxc'
63
64CONTAINS
65
66! **************************************************************************************************
67!> \brief calculates and allocates the xc potential, already reducing it to
68!> the dependence on rho and the one on tau
69!> \param ks_env to get all the needed things
70!> \param rho_struct density for which v_xc is calculated
71!> \param xc_section ...
72!> \param vxc_rho will contain the v_xc part that depend on rho
73!> (if one of the chosen xc functionals has it it is allocated and you
74!> are responsible for it)
75!> \param vxc_tau will contain the kinetic tau part of v_xc
76!> (if one of the chosen xc functionals has it it is allocated and you
77!> are responsible for it)
78!> \param exc ...
79!> \param just_energy if true calculates just the energy, and does not
80!> allocate v_*_rspace
81!> \param edisp ...
82!> \param dispersion_env ...
83!> \param adiabatic_rescale_factor ...
84!> \param pw_env_external external plane wave environment
85!> \par History
86!> - 05.2002 modified to use the mp_allgather function each pe
87!> computes only part of the grid and this is broadcasted to all
88!> instead of summed.
89!> This scales significantly better (e.g. factor 3 on 12 cpus
90!> 32 H2O) [Joost VdV]
91!> - moved to qs_ks_methods [fawzi]
92!> - sic alterations [Joost VandeVondele]
93!> \author Fawzi Mohamed
94! **************************************************************************************************
95 SUBROUTINE qs_vxc_create(ks_env, rho_struct, xc_section, vxc_rho, vxc_tau, exc, &
96 just_energy, edisp, dispersion_env, adiabatic_rescale_factor, &
97 pw_env_external)
98
99 TYPE(qs_ks_env_type), POINTER :: ks_env
100 TYPE(qs_rho_type), POINTER :: rho_struct
101 TYPE(section_vals_type), POINTER :: xc_section
102 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_rho, vxc_tau
103 REAL(kind=dp), INTENT(out) :: exc
104 LOGICAL, INTENT(in), OPTIONAL :: just_energy
105 REAL(kind=dp), INTENT(out), OPTIONAL :: edisp
106 TYPE(qs_dispersion_type), OPTIONAL, POINTER :: dispersion_env
107 REAL(kind=dp), INTENT(in), OPTIONAL :: adiabatic_rescale_factor
108 TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_external
109
110 CHARACTER(len=*), PARAMETER :: routinen = 'qs_vxc_create'
111
112 INTEGER :: handle, ispin, mspin, myfun, &
113 nelec_spin(2), vdw
114 LOGICAL :: compute_virial, do_adiabatic_rescaling, my_just_energy, rho_g_valid, &
115 sic_scaling_b_zero, tau_r_valid, uf_grid, vdw_nl
116 REAL(kind=dp) :: exc_m, factor, &
117 my_adiabatic_rescale_factor, &
118 my_scaling, nelec_s_inv
119 REAL(kind=dp), DIMENSION(3, 3) :: virial_xc_tmp
120 TYPE(cell_type), POINTER :: cell
121 TYPE(dft_control_type), POINTER :: dft_control
122 TYPE(mp_para_env_type), POINTER :: para_env
123 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g, rho_m_gspace, rho_struct_g
124 TYPE(pw_c1d_gs_type), POINTER :: rho_nlcc_g, tmp_g, tmp_g2
125 TYPE(pw_env_type), POINTER :: pw_env
126 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool, vdw_pw_pool, xc_pw_pool
127 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: my_vxc_rho, my_vxc_tau, rho_m_rspace, &
128 rho_r, rho_struct_r, tau, tau_struct_r
129 TYPE(pw_r3d_rs_type), POINTER :: rho_nlcc, tmp_pw, weights
130 TYPE(virial_type), POINTER :: virial
131
132 CALL timeset(routinen, handle)
133
134 cpassert(.NOT. ASSOCIATED(vxc_rho))
135 cpassert(.NOT. ASSOCIATED(vxc_tau))
136 NULLIFY (dft_control, pw_env, auxbas_pw_pool, xc_pw_pool, vdw_pw_pool, cell, my_vxc_rho, &
137 tmp_pw, tmp_g, tmp_g2, my_vxc_tau, rho_g, rho_r, tau, rho_m_rspace, &
138 rho_m_gspace, rho_nlcc, rho_nlcc_g, rho_struct_r, rho_struct_g, tau_struct_r)
139
140 exc = 0.0_dp
141 my_just_energy = .false.
142 IF (PRESENT(just_energy)) my_just_energy = just_energy
143 my_adiabatic_rescale_factor = 1.0_dp
144 do_adiabatic_rescaling = .false.
145 IF (PRESENT(adiabatic_rescale_factor)) THEN
146 my_adiabatic_rescale_factor = adiabatic_rescale_factor
147 do_adiabatic_rescaling = .true.
148 END IF
149
150 CALL get_ks_env(ks_env, &
151 dft_control=dft_control, &
152 pw_env=pw_env, &
153 cell=cell, &
154 xcint_weights=weights, &
155 virial=virial, &
156 rho_nlcc=rho_nlcc, &
157 rho_nlcc_g=rho_nlcc_g)
158
159 CALL qs_rho_get(rho_struct, &
160 tau_r_valid=tau_r_valid, &
161 rho_g_valid=rho_g_valid, &
162 rho_r=rho_struct_r, &
163 rho_g=rho_struct_g, &
164 tau_r=tau_struct_r)
165
166 compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
167 IF (compute_virial) THEN
168 virial%pv_xc = 0.0_dp
169 END IF
170
171 CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", &
172 i_val=myfun)
173 CALL section_vals_val_get(xc_section, "VDW_POTENTIAL%POTENTIAL_TYPE", &
174 i_val=vdw)
175
176 vdw_nl = (vdw == xc_vdw_fun_nonloc)
177 ! this combination has not been investigated
178 cpassert(.NOT. (do_adiabatic_rescaling .AND. vdw_nl))
179 ! are the necessary inputs available
180 IF (.NOT. (PRESENT(dispersion_env) .AND. PRESENT(edisp))) THEN
181 vdw_nl = .false.
182 END IF
183 IF (PRESENT(edisp)) edisp = 0.0_dp
184
185 IF (myfun /= xc_none .OR. vdw_nl) THEN
186
187 ! test if the real space density is available
188 cpassert(ASSOCIATED(rho_struct))
189 IF (dft_control%nspins /= 1 .AND. dft_control%nspins /= 2) &
190 cpabort("nspins must be 1 or 2")
191 mspin = SIZE(rho_struct_r)
192 IF (dft_control%nspins == 2 .AND. mspin == 1) &
193 cpabort("Spin count mismatch")
194
195 ! there are some options related to SIC here.
196 ! Normal DFT computes E(rho_alpha,rho_beta) (or its variant E(2*rho_alpha) for non-LSD)
197 ! SIC can E(rho_alpha,rho_beta)-b*(E(rho_alpha,rho_beta)-E(rho_beta,rho_beta))
198 ! or compute E(rho_alpha,rho_beta)-b*E(rho_alpha-rho_beta,0)
199
200 ! my_scaling is the scaling needed of the standard E(rho_alpha,rho_beta) term
201 my_scaling = 1.0_dp
202 SELECT CASE (dft_control%sic_method_id)
203 CASE (sic_none)
204 ! all fine
205 CASE (sic_mauri_spz, sic_ad)
206 ! no idea yet what to do here in that case
207 cpassert(.NOT. tau_r_valid)
208 CASE (sic_mauri_us)
209 my_scaling = 1.0_dp - dft_control%sic_scaling_b
210 ! no idea yet what to do here in that case
211 cpassert(.NOT. tau_r_valid)
212 CASE (sic_eo)
213 ! NOTHING TO BE DONE
214 CASE DEFAULT
215 ! this case has not yet been treated here
216 cpabort("NYI")
217 END SELECT
218
219 IF (dft_control%sic_scaling_b == 0.0_dp) THEN
220 sic_scaling_b_zero = .true.
221 ELSE
222 sic_scaling_b_zero = .false.
223 END IF
224
225 IF (PRESENT(pw_env_external)) &
226 pw_env => pw_env_external
227 CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool, auxbas_pw_pool=auxbas_pw_pool)
228 uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid, xc_pw_pool%pw_grid)
229
230 IF (.NOT. uf_grid) THEN
231 rho_r => rho_struct_r
232
233 IF (tau_r_valid) THEN
234 tau => tau_struct_r
235 END IF
236
237 ! for gradient corrected functional the density in g space might
238 ! be useful so if we have it, we pass it in
239 IF (rho_g_valid) THEN
240 rho_g => rho_struct_g
241 END IF
242 ELSE
243 cpassert(rho_g_valid)
244 ALLOCATE (rho_r(mspin))
245 ALLOCATE (rho_g(mspin))
246 DO ispin = 1, mspin
247 CALL xc_pw_pool%create_pw(rho_g(ispin))
248 CALL pw_transfer(rho_struct_g(ispin), rho_g(ispin))
249 END DO
250 DO ispin = 1, mspin
251 CALL xc_pw_pool%create_pw(rho_r(ispin))
252 CALL pw_transfer(rho_g(ispin), rho_r(ispin))
253 END DO
254 IF (tau_r_valid) THEN
255 ! tau with finer grids is not implemented (at least not correctly), which this asserts
256 cpabort("Tau (MGGA) with finer grids not implemented")
257 END IF
258 IF (ASSOCIATED(weights)) THEN
259 cpabort("Accurate integration with finer grids not implemented")
260 END IF
261 END IF
262
263 ! add the nlcc densities
264 IF (ASSOCIATED(rho_nlcc)) THEN
265 factor = 1.0_dp
266 DO ispin = 1, mspin
267 CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
268 CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
269 END DO
270 END IF
271
272 !
273 ! here the rho_r, rho_g, tau is what it should be
274 ! we get back the right my_vxc_rho and my_vxc_tau as required
275 !
276 IF (my_just_energy) THEN
277 exc = xc_exc_calc(rho_r=rho_r, tau=tau, &
278 rho_g=rho_g, xc_section=xc_section, &
279 weights=weights, pw_pool=xc_pw_pool)
280
281 ELSE
282 CALL xc_vxc_pw_create(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_r, &
283 rho_g=rho_g, tau=tau, exc=exc, &
284 xc_section=xc_section, &
285 weights=weights, pw_pool=xc_pw_pool, &
286 compute_virial=compute_virial, &
287 virial_xc=virial%pv_xc)
288 END IF
289
290 ! remove the nlcc densities (keep stuff in original state)
291 IF (ASSOCIATED(rho_nlcc)) THEN
292 factor = -1.0_dp
293 DO ispin = 1, mspin
294 CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
295 CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
296 END DO
297 END IF
298
299 ! calclulate non-local vdW functional
300 ! only if this XC_SECTION has it
301 ! if yes, we use the dispersion_env from ks_env
302 ! this is dangerous, as it assumes a special connection xc_section -> qs_env
303 IF (vdw_nl) THEN
304 CALL get_ks_env(ks_env=ks_env, para_env=para_env)
305 ! no SIC functionals allowed
306 cpassert(dft_control%sic_method_id == sic_none)
307 !
308 CALL pw_env_get(pw_env, vdw_pw_pool=vdw_pw_pool)
309 IF (my_just_energy) THEN
310 CALL calculate_dispersion_nonloc(my_vxc_rho, rho_r, rho_g, edisp, dispersion_env, &
311 my_just_energy, vdw_pw_pool, xc_pw_pool, para_env)
312 ELSE
313 CALL calculate_dispersion_nonloc(my_vxc_rho, rho_r, rho_g, edisp, dispersion_env, &
314 my_just_energy, vdw_pw_pool, xc_pw_pool, para_env, virial=virial)
315 END IF
316 END IF
317
318 !! Apply rescaling to the potential if requested
319 IF (.NOT. my_just_energy) THEN
320 IF (do_adiabatic_rescaling) THEN
321 IF (ASSOCIATED(my_vxc_rho)) THEN
322 DO ispin = 1, SIZE(my_vxc_rho)
323 CALL pw_scale(my_vxc_rho(ispin), my_adiabatic_rescale_factor)
324 END DO
325 END IF
326 END IF
327 END IF
328
329 IF (my_scaling /= 1.0_dp) THEN
330 exc = exc*my_scaling
331 IF (ASSOCIATED(my_vxc_rho)) THEN
332 DO ispin = 1, SIZE(my_vxc_rho)
333 CALL pw_scale(my_vxc_rho(ispin), my_scaling)
334 END DO
335 END IF
336 IF (ASSOCIATED(my_vxc_tau)) THEN
337 DO ispin = 1, SIZE(my_vxc_tau)
338 CALL pw_scale(my_vxc_tau(ispin), my_scaling)
339 END DO
340 END IF
341 END IF
342
343 ! we have pw data for the xc, qs_ks requests coeff structure, here we transfer
344 ! pw -> coeff
345 IF (ASSOCIATED(my_vxc_rho)) THEN
346 vxc_rho => my_vxc_rho
347 NULLIFY (my_vxc_rho)
348 END IF
349 IF (ASSOCIATED(my_vxc_tau)) THEN
350 vxc_tau => my_vxc_tau
351 NULLIFY (my_vxc_tau)
352 END IF
353 IF (uf_grid) THEN
354 DO ispin = 1, SIZE(rho_r)
355 CALL xc_pw_pool%give_back_pw(rho_r(ispin))
356 END DO
357 DEALLOCATE (rho_r)
358 IF (ASSOCIATED(rho_g)) THEN
359 DO ispin = 1, SIZE(rho_g)
360 CALL xc_pw_pool%give_back_pw(rho_g(ispin))
361 END DO
362 DEALLOCATE (rho_g)
363 END IF
364 END IF
365
366 ! compute again the xc but now for Exc(m,o) and the opposite sign
367 IF (dft_control%sic_method_id == sic_mauri_spz .AND. .NOT. sic_scaling_b_zero) THEN
368 ALLOCATE (rho_m_rspace(2), rho_m_gspace(2))
369 CALL xc_pw_pool%create_pw(rho_m_gspace(1))
370 CALL xc_pw_pool%create_pw(rho_m_rspace(1))
371 CALL pw_copy(rho_struct_r(1), rho_m_rspace(1))
372 CALL pw_axpy(rho_struct_r(2), rho_m_rspace(1), alpha=-1._dp)
373 CALL pw_copy(rho_struct_g(1), rho_m_gspace(1))
374 CALL pw_axpy(rho_struct_g(2), rho_m_gspace(1), alpha=-1._dp)
375 ! bit sad, these will be just zero...
376 CALL xc_pw_pool%create_pw(rho_m_gspace(2))
377 CALL xc_pw_pool%create_pw(rho_m_rspace(2))
378 CALL pw_zero(rho_m_rspace(2))
379 CALL pw_zero(rho_m_gspace(2))
380
381 IF (my_just_energy) THEN
382 exc_m = xc_exc_calc(rho_r=rho_m_rspace, tau=tau, &
383 rho_g=rho_m_gspace, xc_section=xc_section, &
384 weights=weights, pw_pool=xc_pw_pool)
385 ELSE
386 ! virial untested
387 cpassert(.NOT. compute_virial)
388 CALL xc_vxc_pw_create(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_m_rspace, &
389 rho_g=rho_m_gspace, tau=tau, exc=exc_m, &
390 xc_section=xc_section, &
391 weights=weights, pw_pool=xc_pw_pool, &
392 compute_virial=.false., &
393 virial_xc=virial_xc_tmp)
394 END IF
395
396 exc = exc - dft_control%sic_scaling_b*exc_m
397
398 ! and take care of the potential only vxc_rho is taken into account
399 IF (.NOT. my_just_energy) THEN
400 CALL pw_axpy(my_vxc_rho(1), vxc_rho(1), -dft_control%sic_scaling_b)
401 CALL pw_axpy(my_vxc_rho(1), vxc_rho(2), dft_control%sic_scaling_b)
402 CALL my_vxc_rho(1)%release()
403 CALL my_vxc_rho(2)%release()
404 DEALLOCATE (my_vxc_rho)
405 END IF
406
407 DO ispin = 1, 2
408 CALL xc_pw_pool%give_back_pw(rho_m_rspace(ispin))
409 CALL xc_pw_pool%give_back_pw(rho_m_gspace(ispin))
410 END DO
411 DEALLOCATE (rho_m_rspace)
412 DEALLOCATE (rho_m_gspace)
413
414 END IF
415
416 ! now we have - sum_s N_s * Exc(rho_s/N_s,0)
417 IF (dft_control%sic_method_id == sic_ad .AND. .NOT. sic_scaling_b_zero) THEN
418
419 ! find out how many elecs we have
420 CALL get_ks_env(ks_env, nelectron_spin=nelec_spin)
421
422 ALLOCATE (rho_m_rspace(2), rho_m_gspace(2))
423 DO ispin = 1, 2
424 CALL xc_pw_pool%create_pw(rho_m_gspace(ispin))
425 CALL xc_pw_pool%create_pw(rho_m_rspace(ispin))
426 END DO
427
428 DO ispin = 1, 2
429 IF (nelec_spin(ispin) > 0.0_dp) THEN
430 nelec_s_inv = 1.0_dp/nelec_spin(ispin)
431 ELSE
432 ! does it matter if there are no electrons with this spin (H) ?
433 nelec_s_inv = 0.0_dp
434 END IF
435 CALL pw_copy(rho_struct_r(ispin), rho_m_rspace(1))
436 CALL pw_copy(rho_struct_g(ispin), rho_m_gspace(1))
437 CALL pw_scale(rho_m_rspace(1), nelec_s_inv)
438 CALL pw_scale(rho_m_gspace(1), nelec_s_inv)
439 CALL pw_zero(rho_m_rspace(2))
440 CALL pw_zero(rho_m_gspace(2))
441
442 IF (my_just_energy) THEN
443 exc_m = xc_exc_calc(rho_r=rho_m_rspace, tau=tau, &
444 rho_g=rho_m_gspace, xc_section=xc_section, &
445 weights=weights, pw_pool=xc_pw_pool)
446 ELSE
447 ! virial untested
448 cpassert(.NOT. compute_virial)
449 CALL xc_vxc_pw_create(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_m_rspace, &
450 rho_g=rho_m_gspace, tau=tau, exc=exc_m, &
451 xc_section=xc_section, &
452 weights=weights, pw_pool=xc_pw_pool, &
453 compute_virial=.false., &
454 virial_xc=virial_xc_tmp)
455 END IF
456
457 exc = exc - dft_control%sic_scaling_b*nelec_spin(ispin)*exc_m
458
459 ! and take care of the potential only vxc_rho is taken into account
460 IF (.NOT. my_just_energy) THEN
461 CALL pw_axpy(my_vxc_rho(1), vxc_rho(ispin), -dft_control%sic_scaling_b)
462 CALL my_vxc_rho(1)%release()
463 CALL my_vxc_rho(2)%release()
464 DEALLOCATE (my_vxc_rho)
465 END IF
466 END DO
467
468 DO ispin = 1, 2
469 CALL xc_pw_pool%give_back_pw(rho_m_rspace(ispin))
470 CALL xc_pw_pool%give_back_pw(rho_m_gspace(ispin))
471 END DO
472 DEALLOCATE (rho_m_rspace)
473 DEALLOCATE (rho_m_gspace)
474
475 END IF
476
477 ! compute again the xc but now for Exc(n_down,n_down)
478 IF (dft_control%sic_method_id == sic_mauri_us .AND. .NOT. sic_scaling_b_zero) THEN
479 ALLOCATE (rho_r(2))
480 rho_r(1) = rho_struct_r(2)
481 rho_r(2) = rho_struct_r(2)
482 IF (rho_g_valid) THEN
483 ALLOCATE (rho_g(2))
484 rho_g(1) = rho_struct_g(2)
485 rho_g(2) = rho_struct_g(2)
486 END IF
487
488 IF (my_just_energy) THEN
489 exc_m = xc_exc_calc(rho_r=rho_r, tau=tau, &
490 rho_g=rho_g, xc_section=xc_section, &
491 weights=weights, pw_pool=xc_pw_pool)
492 ELSE
493 ! virial untested
494 cpassert(.NOT. compute_virial)
495 CALL xc_vxc_pw_create(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_r, &
496 rho_g=rho_g, tau=tau, exc=exc_m, &
497 xc_section=xc_section, &
498 weights=weights, pw_pool=xc_pw_pool, &
499 compute_virial=.false., &
500 virial_xc=virial_xc_tmp)
501 END IF
502
503 exc = exc + dft_control%sic_scaling_b*exc_m
504
505 ! and take care of the potential
506 IF (.NOT. my_just_energy) THEN
507 ! both go to minority spin
508 CALL pw_axpy(my_vxc_rho(1), vxc_rho(2), 2.0_dp*dft_control%sic_scaling_b)
509 CALL my_vxc_rho(1)%release()
510 CALL my_vxc_rho(2)%release()
511 DEALLOCATE (my_vxc_rho)
512 END IF
513 DEALLOCATE (rho_r, rho_g)
514
515 END IF
516
517 !
518 ! cleanups
519 !
520 IF (uf_grid .AND. (ASSOCIATED(vxc_rho) .OR. ASSOCIATED(vxc_tau))) THEN
521 block
522 TYPE(pw_r3d_rs_type) :: tmp_pw
523 TYPE(pw_c1d_gs_type) :: tmp_g, tmp_g2
524 CALL xc_pw_pool%create_pw(tmp_g)
525 CALL auxbas_pw_pool%create_pw(tmp_g2)
526 IF (ASSOCIATED(vxc_rho)) THEN
527 DO ispin = 1, SIZE(vxc_rho)
528 CALL auxbas_pw_pool%create_pw(tmp_pw)
529 CALL pw_transfer(vxc_rho(ispin), tmp_g)
530 CALL pw_transfer(tmp_g, tmp_g2)
531 CALL pw_transfer(tmp_g2, tmp_pw)
532 CALL xc_pw_pool%give_back_pw(vxc_rho(ispin))
533 vxc_rho(ispin) = tmp_pw
534 END DO
535 END IF
536 IF (ASSOCIATED(vxc_tau)) THEN
537 DO ispin = 1, SIZE(vxc_tau)
538 CALL auxbas_pw_pool%create_pw(tmp_pw)
539 CALL pw_transfer(vxc_tau(ispin), tmp_g)
540 CALL pw_transfer(tmp_g, tmp_g2)
541 CALL pw_transfer(tmp_g2, tmp_pw)
542 CALL xc_pw_pool%give_back_pw(vxc_tau(ispin))
543 vxc_tau(ispin) = tmp_pw
544 END DO
545 END IF
546 CALL auxbas_pw_pool%give_back_pw(tmp_g2)
547 CALL xc_pw_pool%give_back_pw(tmp_g)
548 END block
549 END IF
550 IF (ASSOCIATED(tau) .AND. uf_grid) THEN
551 DO ispin = 1, SIZE(tau)
552 CALL xc_pw_pool%give_back_pw(tau(ispin))
553 END DO
554 DEALLOCATE (tau)
555 END IF
556
557 END IF
558
559 CALL timestop(handle)
560
561 END SUBROUTINE qs_vxc_create
562
563! **************************************************************************************************
564!> \brief calculates the XC density: E_xc(r) - V_xc(r)*rho(r) or E_xc(r)/rho(r)
565!> \param ks_env to get all the needed things
566!> \param rho_struct density
567!> \param xc_section ...
568!> \param dispersion_env ...
569!> \param xc_ener will contain the xc energy density E_xc(r) - V_xc(r)*rho(r)
570!> \param xc_den will contain the xc energy density E_xc(r)/rho(r)
571!> \param exc will contain the xc energy density E_xc(r)
572!> \param vxc ...
573!> \param vtau ...
574!> \author JGH
575! **************************************************************************************************
576 SUBROUTINE qs_xc_density(ks_env, rho_struct, xc_section, dispersion_env, &
577 xc_ener, xc_den, exc, vxc, vtau)
578
579 TYPE(qs_ks_env_type), POINTER :: ks_env
580 TYPE(qs_rho_type), POINTER :: rho_struct
581 TYPE(section_vals_type), POINTER :: xc_section
582 TYPE(qs_dispersion_type), OPTIONAL, POINTER :: dispersion_env
583 TYPE(pw_r3d_rs_type), INTENT(INOUT), OPTIONAL :: xc_ener, xc_den
584 TYPE(pw_r3d_rs_type), OPTIONAL :: exc
585 TYPE(pw_r3d_rs_type), DIMENSION(:), OPTIONAL :: vxc, vtau
586
587 CHARACTER(len=*), PARAMETER :: routinen = 'qs_xc_density'
588
589 INTEGER :: handle, ispin, myfun, nspins, vdw
590 LOGICAL :: rho_g_valid, tau_r_valid, uf_grid, vdw_nl
591 REAL(kind=dp) :: edisp, excint, factor, rho_cutoff
592 REAL(kind=dp), DIMENSION(3, 3) :: vdum
593 TYPE(cell_type), POINTER :: cell
594 TYPE(dft_control_type), POINTER :: dft_control
595 TYPE(mp_para_env_type), POINTER :: para_env
596 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
597 TYPE(pw_c1d_gs_type), POINTER :: rho_nlcc_g
598 TYPE(pw_env_type), POINTER :: pw_env
599 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool, vdw_pw_pool, xc_pw_pool
600 TYPE(pw_r3d_rs_type) :: exc_r
601 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau_r, vxc_rho, vxc_tau
602 TYPE(pw_r3d_rs_type), POINTER :: rho_nlcc, weights
603
604 CALL timeset(routinen, handle)
605
606 CALL get_ks_env(ks_env, &
607 dft_control=dft_control, &
608 pw_env=pw_env, &
609 cell=cell, &
610 xcint_weights=weights, &
611 rho_nlcc=rho_nlcc, &
612 rho_nlcc_g=rho_nlcc_g)
613
614 CALL qs_rho_get(rho_struct, &
615 tau_r_valid=tau_r_valid, &
616 rho_g_valid=rho_g_valid, &
617 rho_r=rho_r, &
618 rho_g=rho_g, &
619 tau_r=tau_r)
620 nspins = dft_control%nspins
621
622 CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", i_val=myfun)
623 CALL section_vals_val_get(xc_section, "VDW_POTENTIAL%POTENTIAL_TYPE", i_val=vdw)
624 vdw_nl = (vdw == xc_vdw_fun_nonloc)
625 IF (PRESENT(xc_ener)) THEN
626 IF (tau_r_valid) THEN
627 CALL cp_warn(__location__, "Tau contribution will not be correctly handled")
628 END IF
629 END IF
630 IF (vdw_nl) THEN
631 CALL cp_warn(__location__, "vdW functional contribution will be ignored")
632 END IF
633
634 CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool, auxbas_pw_pool=auxbas_pw_pool)
635 uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid, xc_pw_pool%pw_grid)
636 IF (uf_grid) THEN
637 CALL cp_warn(__location__, "Fine grid option not possible with local energy")
638 cpabort("Fine Grid in Local Energy")
639 END IF
640
641 IF (PRESENT(xc_ener)) THEN
642 CALL pw_zero(xc_ener)
643 END IF
644 IF (PRESENT(xc_den)) THEN
645 CALL pw_zero(xc_den)
646 END IF
647 IF (PRESENT(exc)) THEN
648 CALL pw_zero(exc)
649 END IF
650 IF (PRESENT(vxc)) THEN
651 DO ispin = 1, nspins
652 CALL pw_zero(vxc(ispin))
653 END DO
654 END IF
655 IF (PRESENT(vtau)) THEN
656 DO ispin = 1, nspins
657 CALL pw_zero(vtau(ispin))
658 END DO
659 END IF
660
661 IF (myfun /= xc_none) THEN
662
663 cpassert(ASSOCIATED(rho_struct))
664 cpassert(dft_control%sic_method_id == sic_none)
665
666 ! add the nlcc densities
667 IF (ASSOCIATED(rho_nlcc)) THEN
668 factor = 1.0_dp
669 DO ispin = 1, nspins
670 CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
671 CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
672 END DO
673 END IF
674 NULLIFY (vxc_rho, vxc_tau)
675 CALL xc_vxc_pw_create(vxc_rho=vxc_rho, vxc_tau=vxc_tau, rho_r=rho_r, &
676 rho_g=rho_g, tau=tau_r, exc=excint, &
677 xc_section=xc_section, &
678 weights=weights, pw_pool=xc_pw_pool, &
679 compute_virial=.false., &
680 virial_xc=vdum, &
681 exc_r=exc_r)
682 ! calclulate non-local vdW functional
683 ! only if this XC_SECTION has it
684 ! if yes, we use the dispersion_env from ks_env
685 ! this is dangerous, as it assumes a special connection xc_section -> qs_env
686 IF (vdw_nl) THEN
687 CALL get_ks_env(ks_env=ks_env, para_env=para_env)
688 ! no SIC functionals allowed
689 cpassert(dft_control%sic_method_id == sic_none)
690 !
691 CALL pw_env_get(pw_env, vdw_pw_pool=vdw_pw_pool)
692 CALL calculate_dispersion_nonloc(vxc_rho, rho_r, rho_g, edisp, dispersion_env, &
693 .false., vdw_pw_pool, xc_pw_pool, para_env)
694 END IF
695
696 ! remove the nlcc densities (keep stuff in original state)
697 IF (ASSOCIATED(rho_nlcc)) THEN
698 factor = -1.0_dp
699 DO ispin = 1, dft_control%nspins
700 CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
701 CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
702 END DO
703 END IF
704 !
705 IF (PRESENT(xc_den)) THEN
706 CALL pw_copy(exc_r, xc_den)
707 rho_cutoff = 1.e-14_dp
708 CALL calc_xc_density(xc_den, rho_r, rho_cutoff)
709 END IF
710 IF (PRESENT(xc_ener)) THEN
711 CALL pw_copy(exc_r, xc_ener)
712 DO ispin = 1, nspins
713 CALL pw_multiply(xc_ener, vxc_rho(ispin), rho_r(ispin), alpha=-1.0_dp)
714 END DO
715 END IF
716 IF (PRESENT(exc)) THEN
717 CALL pw_copy(exc_r, exc)
718 END IF
719 IF (PRESENT(vxc)) THEN
720 DO ispin = 1, nspins
721 CALL pw_copy(vxc_rho(ispin), vxc(ispin))
722 END DO
723 END IF
724 IF (PRESENT(vtau)) THEN
725 DO ispin = 1, nspins
726 CALL pw_copy(vxc_tau(ispin), vtau(ispin))
727 END DO
728 END IF
729 ! remove arrays
730 IF (ASSOCIATED(vxc_rho)) THEN
731 DO ispin = 1, nspins
732 CALL vxc_rho(ispin)%release()
733 END DO
734 DEALLOCATE (vxc_rho)
735 END IF
736 IF (ASSOCIATED(vxc_tau)) THEN
737 DO ispin = 1, nspins
738 CALL vxc_tau(ispin)%release()
739 END DO
740 DEALLOCATE (vxc_tau)
741 END IF
742 CALL exc_r%release()
743 !
744 END IF
745
746 CALL timestop(handle)
747
748 END SUBROUTINE qs_xc_density
749
750END MODULE qs_vxc
pw_methods::pw_axpy
Definition pw_methods.F:165
pw_methods::pw_copy
Definition pw_methods.F:146
pw_methods::pw_multiply
Definition pw_methods.F:184
pw_methods::pw_scale
Definition pw_methods.F:94
pw_methods::pw_transfer
Definition pw_methods.F:304
pw_methods::pw_zero
Definition pw_methods.F:83
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::sic_mauri_spz
integer, parameter, public sic_mauri_spz
Definition input_constants.F:572
input_constants::xc_vdw_fun_nonloc
integer, parameter, public xc_vdw_fun_nonloc
Definition input_constants.F:599
input_constants::sic_eo
integer, parameter, public sic_eo
Definition input_constants.F:572
input_constants::sic_mauri_us
integer, parameter, public sic_mauri_us
Definition input_constants.F:572
input_constants::sic_none
integer, parameter, public sic_none
Definition input_constants.F:572
input_constants::xc_none
integer, parameter, public xc_none
Definition input_constants.F:555
input_constants::sic_ad
integer, parameter, public sic_ad
Definition input_constants.F:572
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_val_get
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
Definition input_section_types.F:1047
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_env_types::pw_env_get
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
Definition pw_env_types.F:113
pw_grids
This module defines the grid data type and some basic operations on it.
Definition pw_grids.F:36
pw_grids::pw_grid_compare
logical function, public pw_grid_compare(grida, gridb)
Check if two pw_grids are equal.
Definition pw_grids.F:148
pw_methods
Definition pw_methods.F:25
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_types
Definition pw_types.F:24
qs_dispersion_nonloc
Calculation of non local dispersion functionals Some routines adapted from: Copyright (C) 2001-2009 Q...
Definition qs_dispersion_nonloc.F:19
qs_dispersion_nonloc::calculate_dispersion_nonloc
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 qs_dispersion_nonloc.F:166
qs_dispersion_types
Definition of disperson types for DFT calculations.
Definition qs_dispersion_types.F:12
qs_ks_types
Definition qs_ks_types.F:15
qs_ks_types::get_ks_env
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)
...
Definition qs_ks_types.F:341
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
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...
Definition qs_rho_types.F:229
qs_vxc
Definition qs_vxc.F:16
qs_vxc::qs_vxc_create
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)
calculates and allocates the xc potential, already reducing it to the dependence on rho and the one o...
Definition qs_vxc.F:98
qs_vxc::qs_xc_density
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:578
virial_types
Definition virial_types.F:13
xc
Exchange and Correlation functional calculations.
Definition xc.F:17
xc::xc_exc_calc
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:787
xc::xc_vxc_pw_create
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:470
xc::calc_xc_density
subroutine, public calc_xc_density(pot, rho, rho_cutoff)
Definition xc.F:398
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_control_types::dft_control_type
Definition cp_control_types.F:631
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:743
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_dispersion_types::qs_dispersion_type
Definition qs_dispersion_types.F:34
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:137
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
virial_types::virial_type
Definition virial_types.F:26