(git:e68414f)
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-2025 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_create1
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
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 virial=virial, &
155 rho_nlcc=rho_nlcc, &
156 rho_nlcc_g=rho_nlcc_g)
157
158 CALL qs_rho_get(rho_struct, &
159 tau_r_valid=tau_r_valid, &
160 rho_g_valid=rho_g_valid, &
161 rho_r=rho_struct_r, &
162 rho_g=rho_struct_g, &
163 tau_r=tau_struct_r)
164
165 compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
166 IF (compute_virial) THEN
167 virial%pv_xc = 0.0_dp
168 END IF
169
170 CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", &
171 i_val=myfun)
172 CALL section_vals_val_get(xc_section, "VDW_POTENTIAL%POTENTIAL_TYPE", &
173 i_val=vdw)
174
175 vdw_nl = (vdw == xc_vdw_fun_nonloc)
176 ! this combination has not been investigated
177 cpassert(.NOT. (do_adiabatic_rescaling .AND. vdw_nl))
178 ! are the necessary inputs available
179 IF (.NOT. (PRESENT(dispersion_env) .AND. PRESENT(edisp))) THEN
180 vdw_nl = .false.
181 END IF
182 IF (PRESENT(edisp)) edisp = 0.0_dp
183
184 IF (myfun /= xc_none .OR. vdw_nl) THEN
185
186 ! test if the real space density is available
187 cpassert(ASSOCIATED(rho_struct))
188 IF (dft_control%nspins /= 1 .AND. dft_control%nspins /= 2) &
189 cpabort("nspins must be 1 or 2")
190 mspin = SIZE(rho_struct_r)
191 IF (dft_control%nspins == 2 .AND. mspin == 1) &
192 cpabort("Spin count mismatch")
193
194 ! there are some options related to SIC here.
195 ! Normal DFT computes E(rho_alpha,rho_beta) (or its variant E(2*rho_alpha) for non-LSD)
196 ! SIC can E(rho_alpha,rho_beta)-b*(E(rho_alpha,rho_beta)-E(rho_beta,rho_beta))
197 ! or compute E(rho_alpha,rho_beta)-b*E(rho_alpha-rho_beta,0)
198
199 ! my_scaling is the scaling needed of the standard E(rho_alpha,rho_beta) term
200 my_scaling = 1.0_dp
201 SELECT CASE (dft_control%sic_method_id)
202 CASE (sic_none)
203 ! all fine
204 CASE (sic_mauri_spz, sic_ad)
205 ! no idea yet what to do here in that case
206 cpassert(.NOT. tau_r_valid)
207 CASE (sic_mauri_us)
208 my_scaling = 1.0_dp - dft_control%sic_scaling_b
209 ! no idea yet what to do here in that case
210 cpassert(.NOT. tau_r_valid)
211 CASE (sic_eo)
212 ! NOTHING TO BE DONE
213 CASE DEFAULT
214 ! this case has not yet been treated here
215 cpabort("NYI")
216 END SELECT
217
218 IF (dft_control%sic_scaling_b .EQ. 0.0_dp) THEN
219 sic_scaling_b_zero = .true.
220 ELSE
221 sic_scaling_b_zero = .false.
222 END IF
223
224 IF (PRESENT(pw_env_external)) &
225 pw_env => pw_env_external
226 CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool, auxbas_pw_pool=auxbas_pw_pool)
227 uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid, xc_pw_pool%pw_grid)
228
229 IF (.NOT. uf_grid) THEN
230 rho_r => rho_struct_r
231
232 IF (tau_r_valid) THEN
233 tau => tau_struct_r
234 END IF
235
236 ! for gradient corrected functional the density in g space might
237 ! be useful so if we have it, we pass it in
238 IF (rho_g_valid) THEN
239 rho_g => rho_struct_g
240 END IF
241 ELSE
242 cpassert(rho_g_valid)
243 ALLOCATE (rho_r(mspin))
244 ALLOCATE (rho_g(mspin))
245 DO ispin = 1, mspin
246 CALL xc_pw_pool%create_pw(rho_g(ispin))
247 CALL pw_transfer(rho_struct_g(ispin), rho_g(ispin))
248 END DO
249 DO ispin = 1, mspin
250 CALL xc_pw_pool%create_pw(rho_r(ispin))
251 CALL pw_transfer(rho_g(ispin), rho_r(ispin))
252 END DO
253 IF (tau_r_valid) THEN
254 ! tau with finer grids is not implemented (at least not correctly), which this asserts
255 cpabort("tau with finer grids not implemented")
256 END IF
257 END IF
258
259 ! add the nlcc densities
260 IF (ASSOCIATED(rho_nlcc)) THEN
261 factor = 1.0_dp
262 DO ispin = 1, mspin
263 CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
264 CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
265 END DO
266 END IF
267
268 !
269 ! here the rho_r, rho_g, tau is what it should be
270 ! we get back the right my_vxc_rho and my_vxc_tau as required
271 !
272 IF (my_just_energy) THEN
273 exc = xc_exc_calc(rho_r=rho_r, tau=tau, &
274 rho_g=rho_g, xc_section=xc_section, &
275 pw_pool=xc_pw_pool)
276
277 ELSE
278 CALL xc_vxc_pw_create1(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_r, &
279 rho_g=rho_g, tau=tau, exc=exc, &
280 xc_section=xc_section, &
281 pw_pool=xc_pw_pool, &
282 compute_virial=compute_virial, &
283 virial_xc=virial%pv_xc)
284 END IF
285
286 ! remove the nlcc densities (keep stuff in original state)
287 IF (ASSOCIATED(rho_nlcc)) THEN
288 factor = -1.0_dp
289 DO ispin = 1, mspin
290 CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
291 CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
292 END DO
293 END IF
294
295 ! calclulate non-local vdW functional
296 ! only if this XC_SECTION has it
297 ! if yes, we use the dispersion_env from ks_env
298 ! this is dangerous, as it assumes a special connection xc_section -> qs_env
299 IF (vdw_nl) THEN
300 CALL get_ks_env(ks_env=ks_env, para_env=para_env)
301 ! no SIC functionals allowed
302 cpassert(dft_control%sic_method_id == sic_none)
303 !
304 CALL pw_env_get(pw_env, vdw_pw_pool=vdw_pw_pool)
305 IF (my_just_energy) THEN
306 CALL calculate_dispersion_nonloc(my_vxc_rho, rho_r, rho_g, edisp, dispersion_env, &
307 my_just_energy, vdw_pw_pool, xc_pw_pool, para_env)
308 ELSE
309 CALL calculate_dispersion_nonloc(my_vxc_rho, rho_r, rho_g, edisp, dispersion_env, &
310 my_just_energy, vdw_pw_pool, xc_pw_pool, para_env, virial=virial)
311 END IF
312 END IF
313
314 !! Apply rescaling to the potential if requested
315 IF (.NOT. my_just_energy) THEN
316 IF (do_adiabatic_rescaling) THEN
317 IF (ASSOCIATED(my_vxc_rho)) THEN
318 DO ispin = 1, SIZE(my_vxc_rho)
319 CALL pw_scale(my_vxc_rho(ispin), my_adiabatic_rescale_factor)
320 END DO
321 END IF
322 END IF
323 END IF
324
325 IF (my_scaling .NE. 1.0_dp) THEN
326 exc = exc*my_scaling
327 IF (ASSOCIATED(my_vxc_rho)) THEN
328 DO ispin = 1, SIZE(my_vxc_rho)
329 CALL pw_scale(my_vxc_rho(ispin), my_scaling)
330 END DO
331 END IF
332 IF (ASSOCIATED(my_vxc_tau)) THEN
333 DO ispin = 1, SIZE(my_vxc_tau)
334 CALL pw_scale(my_vxc_tau(ispin), my_scaling)
335 END DO
336 END IF
337 END IF
338
339 ! we have pw data for the xc, qs_ks requests coeff structure, here we transfer
340 ! pw -> coeff
341 IF (ASSOCIATED(my_vxc_rho)) THEN
342 vxc_rho => my_vxc_rho
343 NULLIFY (my_vxc_rho)
344 END IF
345 IF (ASSOCIATED(my_vxc_tau)) THEN
346 vxc_tau => my_vxc_tau
347 NULLIFY (my_vxc_tau)
348 END IF
349 IF (uf_grid) THEN
350 DO ispin = 1, SIZE(rho_r)
351 CALL xc_pw_pool%give_back_pw(rho_r(ispin))
352 END DO
353 DEALLOCATE (rho_r)
354 IF (ASSOCIATED(rho_g)) THEN
355 DO ispin = 1, SIZE(rho_g)
356 CALL xc_pw_pool%give_back_pw(rho_g(ispin))
357 END DO
358 DEALLOCATE (rho_g)
359 END IF
360 END IF
361
362 ! compute again the xc but now for Exc(m,o) and the opposite sign
363 IF (dft_control%sic_method_id .EQ. sic_mauri_spz .AND. .NOT. sic_scaling_b_zero) THEN
364 ALLOCATE (rho_m_rspace(2), rho_m_gspace(2))
365 CALL xc_pw_pool%create_pw(rho_m_gspace(1))
366 CALL xc_pw_pool%create_pw(rho_m_rspace(1))
367 CALL pw_copy(rho_struct_r(1), rho_m_rspace(1))
368 CALL pw_axpy(rho_struct_r(2), rho_m_rspace(1), alpha=-1._dp)
369 CALL pw_copy(rho_struct_g(1), rho_m_gspace(1))
370 CALL pw_axpy(rho_struct_g(2), rho_m_gspace(1), alpha=-1._dp)
371 ! bit sad, these will be just zero...
372 CALL xc_pw_pool%create_pw(rho_m_gspace(2))
373 CALL xc_pw_pool%create_pw(rho_m_rspace(2))
374 CALL pw_zero(rho_m_rspace(2))
375 CALL pw_zero(rho_m_gspace(2))
376
377 IF (my_just_energy) THEN
378 exc_m = xc_exc_calc(rho_r=rho_m_rspace, tau=tau, &
379 rho_g=rho_m_gspace, xc_section=xc_section, &
380 pw_pool=xc_pw_pool)
381 ELSE
382 ! virial untested
383 cpassert(.NOT. compute_virial)
384 CALL xc_vxc_pw_create1(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_m_rspace, &
385 rho_g=rho_m_gspace, tau=tau, exc=exc_m, &
386 xc_section=xc_section, &
387 pw_pool=xc_pw_pool, &
388 compute_virial=.false., &
389 virial_xc=virial_xc_tmp)
390 END IF
391
392 exc = exc - dft_control%sic_scaling_b*exc_m
393
394 ! and take care of the potential only vxc_rho is taken into account
395 IF (.NOT. my_just_energy) THEN
396 CALL pw_axpy(my_vxc_rho(1), vxc_rho(1), -dft_control%sic_scaling_b)
397 CALL pw_axpy(my_vxc_rho(1), vxc_rho(2), dft_control%sic_scaling_b)
398 CALL my_vxc_rho(1)%release()
399 CALL my_vxc_rho(2)%release()
400 DEALLOCATE (my_vxc_rho)
401 END IF
402
403 DO ispin = 1, 2
404 CALL xc_pw_pool%give_back_pw(rho_m_rspace(ispin))
405 CALL xc_pw_pool%give_back_pw(rho_m_gspace(ispin))
406 END DO
407 DEALLOCATE (rho_m_rspace)
408 DEALLOCATE (rho_m_gspace)
409
410 END IF
411
412 ! now we have - sum_s N_s * Exc(rho_s/N_s,0)
413 IF (dft_control%sic_method_id .EQ. sic_ad .AND. .NOT. sic_scaling_b_zero) THEN
414
415 ! find out how many elecs we have
416 CALL get_ks_env(ks_env, nelectron_spin=nelec_spin)
417
418 ALLOCATE (rho_m_rspace(2), rho_m_gspace(2))
419 DO ispin = 1, 2
420 CALL xc_pw_pool%create_pw(rho_m_gspace(ispin))
421 CALL xc_pw_pool%create_pw(rho_m_rspace(ispin))
422 END DO
423
424 DO ispin = 1, 2
425 IF (nelec_spin(ispin) .GT. 0.0_dp) THEN
426 nelec_s_inv = 1.0_dp/nelec_spin(ispin)
427 ELSE
428 ! does it matter if there are no electrons with this spin (H) ?
429 nelec_s_inv = 0.0_dp
430 END IF
431 CALL pw_copy(rho_struct_r(ispin), rho_m_rspace(1))
432 CALL pw_copy(rho_struct_g(ispin), rho_m_gspace(1))
433 CALL pw_scale(rho_m_rspace(1), nelec_s_inv)
434 CALL pw_scale(rho_m_gspace(1), nelec_s_inv)
435 CALL pw_zero(rho_m_rspace(2))
436 CALL pw_zero(rho_m_gspace(2))
437
438 IF (my_just_energy) THEN
439 exc_m = xc_exc_calc(rho_r=rho_m_rspace, tau=tau, &
440 rho_g=rho_m_gspace, xc_section=xc_section, &
441 pw_pool=xc_pw_pool)
442 ELSE
443 ! virial untested
444 cpassert(.NOT. compute_virial)
445 CALL xc_vxc_pw_create1(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_m_rspace, &
446 rho_g=rho_m_gspace, tau=tau, exc=exc_m, &
447 xc_section=xc_section, &
448 pw_pool=xc_pw_pool, &
449 compute_virial=.false., &
450 virial_xc=virial_xc_tmp)
451 END IF
452
453 exc = exc - dft_control%sic_scaling_b*nelec_spin(ispin)*exc_m
454
455 ! and take care of the potential only vxc_rho is taken into account
456 IF (.NOT. my_just_energy) THEN
457 CALL pw_axpy(my_vxc_rho(1), vxc_rho(ispin), -dft_control%sic_scaling_b)
458 CALL my_vxc_rho(1)%release()
459 CALL my_vxc_rho(2)%release()
460 DEALLOCATE (my_vxc_rho)
461 END IF
462 END DO
463
464 DO ispin = 1, 2
465 CALL xc_pw_pool%give_back_pw(rho_m_rspace(ispin))
466 CALL xc_pw_pool%give_back_pw(rho_m_gspace(ispin))
467 END DO
468 DEALLOCATE (rho_m_rspace)
469 DEALLOCATE (rho_m_gspace)
470
471 END IF
472
473 ! compute again the xc but now for Exc(n_down,n_down)
474 IF (dft_control%sic_method_id .EQ. sic_mauri_us .AND. .NOT. sic_scaling_b_zero) THEN
475 ALLOCATE (rho_r(2))
476 rho_r(1) = rho_struct_r(2)
477 rho_r(2) = rho_struct_r(2)
478 IF (rho_g_valid) THEN
479 ALLOCATE (rho_g(2))
480 rho_g(1) = rho_struct_g(2)
481 rho_g(2) = rho_struct_g(2)
482 END IF
483
484 IF (my_just_energy) THEN
485 exc_m = xc_exc_calc(rho_r=rho_r, tau=tau, &
486 rho_g=rho_g, xc_section=xc_section, &
487 pw_pool=xc_pw_pool)
488 ELSE
489 ! virial untested
490 cpassert(.NOT. compute_virial)
491 CALL xc_vxc_pw_create1(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_r, &
492 rho_g=rho_g, tau=tau, exc=exc_m, &
493 xc_section=xc_section, &
494 pw_pool=xc_pw_pool, &
495 compute_virial=.false., &
496 virial_xc=virial_xc_tmp)
497 END IF
498
499 exc = exc + dft_control%sic_scaling_b*exc_m
500
501 ! and take care of the potential
502 IF (.NOT. my_just_energy) THEN
503 ! both go to minority spin
504 CALL pw_axpy(my_vxc_rho(1), vxc_rho(2), 2.0_dp*dft_control%sic_scaling_b)
505 CALL my_vxc_rho(1)%release()
506 CALL my_vxc_rho(2)%release()
507 DEALLOCATE (my_vxc_rho)
508 END IF
509 DEALLOCATE (rho_r, rho_g)
510
511 END IF
512
513 !
514 ! cleanups
515 !
516 IF (uf_grid .AND. (ASSOCIATED(vxc_rho) .OR. ASSOCIATED(vxc_tau))) THEN
517 block
518 TYPE(pw_r3d_rs_type) :: tmp_pw
519 TYPE(pw_c1d_gs_type) :: tmp_g, tmp_g2
520 CALL xc_pw_pool%create_pw(tmp_g)
521 CALL auxbas_pw_pool%create_pw(tmp_g2)
522 IF (ASSOCIATED(vxc_rho)) THEN
523 DO ispin = 1, SIZE(vxc_rho)
524 CALL auxbas_pw_pool%create_pw(tmp_pw)
525 CALL pw_transfer(vxc_rho(ispin), tmp_g)
526 CALL pw_transfer(tmp_g, tmp_g2)
527 CALL pw_transfer(tmp_g2, tmp_pw)
528 CALL xc_pw_pool%give_back_pw(vxc_rho(ispin))
529 vxc_rho(ispin) = tmp_pw
530 END DO
531 END IF
532 IF (ASSOCIATED(vxc_tau)) THEN
533 DO ispin = 1, SIZE(vxc_tau)
534 CALL auxbas_pw_pool%create_pw(tmp_pw)
535 CALL pw_transfer(vxc_tau(ispin), tmp_g)
536 CALL pw_transfer(tmp_g, tmp_g2)
537 CALL pw_transfer(tmp_g2, tmp_pw)
538 CALL xc_pw_pool%give_back_pw(vxc_tau(ispin))
539 vxc_tau(ispin) = tmp_pw
540 END DO
541 END IF
542 CALL auxbas_pw_pool%give_back_pw(tmp_g2)
543 CALL xc_pw_pool%give_back_pw(tmp_g)
544 END block
545 END IF
546 IF (ASSOCIATED(tau) .AND. uf_grid) THEN
547 DO ispin = 1, SIZE(tau)
548 CALL xc_pw_pool%give_back_pw(tau(ispin))
549 END DO
550 DEALLOCATE (tau)
551 END IF
552
553 END IF
554
555 CALL timestop(handle)
556
557 END SUBROUTINE qs_vxc_create
558
559! **************************************************************************************************
560!> \brief calculates the XC density: E_xc(r) - V_xc(r)*rho(r) or E_xc(r)/rho(r)
561!> \param ks_env to get all the needed things
562!> \param rho_struct density
563!> \param xc_section ...
564!> \param dispersion_env ...
565!> \param xc_ener will contain the xc energy density E_xc(r) - V_xc(r)*rho(r)
566!> \param xc_den will contain the xc energy density E_xc(r)/rho(r)
567!> \param vxc ...
568!> \param vtau ...
569!> \author JGH
570! **************************************************************************************************
571 SUBROUTINE qs_xc_density(ks_env, rho_struct, xc_section, dispersion_env, &
572 xc_ener, xc_den, vxc, vtau)
573
574 TYPE(qs_ks_env_type), POINTER :: ks_env
575 TYPE(qs_rho_type), POINTER :: rho_struct
576 TYPE(section_vals_type), POINTER :: xc_section
577 TYPE(qs_dispersion_type), OPTIONAL, POINTER :: dispersion_env
578 TYPE(pw_r3d_rs_type), INTENT(INOUT), OPTIONAL :: xc_ener, xc_den
579 TYPE(pw_r3d_rs_type), DIMENSION(:), OPTIONAL :: vxc, vtau
580
581 CHARACTER(len=*), PARAMETER :: routinen = 'qs_xc_density'
582
583 INTEGER :: handle, ispin, myfun, nspins, vdw
584 LOGICAL :: rho_g_valid, tau_r_valid, uf_grid, vdw_nl
585 REAL(kind=dp) :: edisp, exc, factor, rho_cutoff
586 REAL(kind=dp), DIMENSION(3, 3) :: vdum
587 TYPE(cell_type), POINTER :: cell
588 TYPE(dft_control_type), POINTER :: dft_control
589 TYPE(mp_para_env_type), POINTER :: para_env
590 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
591 TYPE(pw_c1d_gs_type), POINTER :: rho_nlcc_g
592 TYPE(pw_env_type), POINTER :: pw_env
593 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool, vdw_pw_pool, xc_pw_pool
594 TYPE(pw_r3d_rs_type) :: exc_r
595 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau_r, vxc_rho, vxc_tau
596 TYPE(pw_r3d_rs_type), POINTER :: rho_nlcc
597
598 CALL timeset(routinen, handle)
599
600 CALL get_ks_env(ks_env, &
601 dft_control=dft_control, &
602 pw_env=pw_env, &
603 cell=cell, &
604 rho_nlcc=rho_nlcc, &
605 rho_nlcc_g=rho_nlcc_g)
606
607 CALL qs_rho_get(rho_struct, &
608 tau_r_valid=tau_r_valid, &
609 rho_g_valid=rho_g_valid, &
610 rho_r=rho_r, &
611 rho_g=rho_g, &
612 tau_r=tau_r)
613 nspins = dft_control%nspins
614
615 CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", i_val=myfun)
616 CALL section_vals_val_get(xc_section, "VDW_POTENTIAL%POTENTIAL_TYPE", i_val=vdw)
617 vdw_nl = (vdw == xc_vdw_fun_nonloc)
618 IF (PRESENT(xc_ener)) THEN
619 IF (tau_r_valid) THEN
620 CALL cp_warn(__location__, "Tau contribution will not be correctly handled")
621 END IF
622 END IF
623 IF (vdw_nl) THEN
624 CALL cp_warn(__location__, "vdW functional contribution will be ignored")
625 END IF
626
627 CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool, auxbas_pw_pool=auxbas_pw_pool)
628 uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid, xc_pw_pool%pw_grid)
629 IF (uf_grid) THEN
630 CALL cp_warn(__location__, "Fine grid option not possible with local energy")
631 cpabort("Fine Grid in Local Energy")
632 END IF
633
634 IF (PRESENT(xc_ener)) THEN
635 CALL pw_zero(xc_ener)
636 END IF
637 IF (PRESENT(xc_den)) THEN
638 CALL pw_zero(xc_den)
639 END IF
640 IF (PRESENT(vxc)) THEN
641 DO ispin = 1, nspins
642 CALL pw_zero(vxc(ispin))
643 END DO
644 END IF
645 IF (PRESENT(vtau)) THEN
646 DO ispin = 1, nspins
647 CALL pw_zero(vtau(ispin))
648 END DO
649 END IF
650
651 IF (myfun /= xc_none) THEN
652
653 cpassert(ASSOCIATED(rho_struct))
654 cpassert(dft_control%sic_method_id == sic_none)
655
656 ! add the nlcc densities
657 IF (ASSOCIATED(rho_nlcc)) THEN
658 factor = 1.0_dp
659 DO ispin = 1, nspins
660 CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
661 CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
662 END DO
663 END IF
664 NULLIFY (vxc_rho, vxc_tau)
665 CALL xc_vxc_pw_create1(vxc_rho=vxc_rho, vxc_tau=vxc_tau, rho_r=rho_r, &
666 rho_g=rho_g, tau=tau_r, exc=exc, &
667 xc_section=xc_section, &
668 pw_pool=xc_pw_pool, &
669 compute_virial=.false., &
670 virial_xc=vdum, &
671 exc_r=exc_r)
672 ! calclulate non-local vdW functional
673 ! only if this XC_SECTION has it
674 ! if yes, we use the dispersion_env from ks_env
675 ! this is dangerous, as it assumes a special connection xc_section -> qs_env
676 IF (vdw_nl) THEN
677 CALL get_ks_env(ks_env=ks_env, para_env=para_env)
678 ! no SIC functionals allowed
679 cpassert(dft_control%sic_method_id == sic_none)
680 !
681 CALL pw_env_get(pw_env, vdw_pw_pool=vdw_pw_pool)
682 CALL calculate_dispersion_nonloc(vxc_rho, rho_r, rho_g, edisp, dispersion_env, &
683 .false., vdw_pw_pool, xc_pw_pool, para_env)
684 END IF
685
686 ! remove the nlcc densities (keep stuff in original state)
687 IF (ASSOCIATED(rho_nlcc)) THEN
688 factor = -1.0_dp
689 DO ispin = 1, dft_control%nspins
690 CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
691 CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
692 END DO
693 END IF
694 !
695 IF (PRESENT(xc_den)) THEN
696 CALL pw_copy(exc_r, xc_den)
697 rho_cutoff = 1.e-14_dp
698 CALL calc_xc_density(xc_den, rho_r, rho_cutoff)
699 END IF
700 IF (PRESENT(xc_ener)) THEN
701 CALL pw_copy(exc_r, xc_ener)
702 DO ispin = 1, nspins
703 CALL pw_multiply(xc_ener, vxc_rho(ispin), rho_r(ispin), alpha=-1.0_dp)
704 END DO
705 END IF
706 IF (PRESENT(vxc)) THEN
707 DO ispin = 1, nspins
708 CALL pw_copy(vxc_rho(ispin), vxc(ispin))
709 END DO
710 END IF
711 IF (PRESENT(vtau)) THEN
712 DO ispin = 1, nspins
713 CALL pw_copy(vxc_tau(ispin), vtau(ispin))
714 END DO
715 END IF
716 ! remove arrays
717 IF (ASSOCIATED(vxc_rho)) THEN
718 DO ispin = 1, nspins
719 CALL vxc_rho(ispin)%release()
720 END DO
721 DEALLOCATE (vxc_rho)
722 END IF
723 IF (ASSOCIATED(vxc_tau)) THEN
724 DO ispin = 1, nspins
725 CALL vxc_tau(ispin)%release()
726 END DO
727 DEALLOCATE (vxc_tau)
728 END IF
729 CALL exc_r%release()
730 !
731 END IF
732
733 CALL timestop(handle)
734
735 END SUBROUTINE qs_xc_density
736
737END 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:570
input_constants::xc_vdw_fun_nonloc
integer, parameter, public xc_vdw_fun_nonloc
Definition input_constants.F:597
input_constants::sic_eo
integer, parameter, public sic_eo
Definition input_constants.F:570
input_constants::sic_mauri_us
integer, parameter, public sic_mauri_us
Definition input_constants.F:570
input_constants::sic_none
integer, parameter, public sic_none
Definition input_constants.F:570
input_constants::xc_none
integer, parameter, public xc_none
Definition input_constants.F:553
input_constants::sic_ad
integer, parameter, public sic_ad
Definition input_constants.F:570
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, 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, 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:336
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, vxc, vtau)
calculates the XC density: E_xc(r) - V_xc(r)*rho(r) or E_xc(r)/rho(r)
Definition qs_vxc.F:573
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, pw_pool)
calculates just the exchange and correlation energy (no vxc)
Definition xc.F:1460
xc::calc_xc_density
subroutine, public calc_xc_density(pot, rho, rho_cutoff)
Definition xc.F:1074
xc::xc_vxc_pw_create1
subroutine, public xc_vxc_pw_create1(vxc_rho, vxc_tau, rho_r, rho_g, tau, exc, xc_section, pw_pool, compute_virial, virial_xc, exc_r)
Exchange and Correlation functional calculations. depending on the selected functional_routine calls ...
Definition xc.F:150
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_control_types::dft_control_type
Definition cp_control_types.F:602
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:763
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