(git:0de0cc2)
qs_vxc.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 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 ! **************************************************************************************************
16 MODULE 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,&
23  sic_mauri_us,&
24  sic_none,&
25  xc_none,&
27  USE input_section_types, ONLY: section_vals_type,&
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
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,&
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 
64 CONTAINS
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 
737 END 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.
container for various plainwaves related things
Definition: pw_env_types.F:14
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
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:147
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
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, 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_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, 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:330
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
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
Definition: qs_vxc.F:16
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
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
Exchange and Correlation functional calculations.
Definition: xc.F:17
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:1449
subroutine, public calc_xc_density(pot, rho, rho_cutoff)
Definition: xc.F:1062
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