(git:c641773)
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 cp_dbcsr_api, ONLY: dbcsr_p_type
21 USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
22 USE gapw_gspace_reconstruction, ONLY: calculate_rhotot_elec_gspace
23 USE input_constants, ONLY: sic_ad,&
24 sic_eo,&
25 sic_mauri_spz,&
26 sic_mauri_us,&
27 sic_none,&
28 xc_none,&
29 xc_vdw_fun_nonloc
30 USE input_section_types, ONLY: section_vals_type,&
31 section_vals_val_get
32 USE kinds, ONLY: dp
33 USE message_passing, ONLY: mp_para_env_type
34 USE particle_types, ONLY: particle_type
35 USE pw_env_types, ONLY: pw_env_get,&
36 pw_env_type
37 USE pw_grids, ONLY: get_pw_grid_info,&
38 pw_grid_compare
39 USE pw_methods, ONLY: pw_axpy,&
40 pw_copy,&
41 pw_integral_ab,&
42 pw_integrate_function,&
43 pw_multiply,&
44 pw_scale,&
45 pw_transfer,&
46 pw_zero
47 USE pw_pool_types, ONLY: pw_pool_type
48 USE pw_types, ONLY: pw_c1d_gs_type,&
49 pw_r3d_rs_type
50 USE qs_collocate_density, ONLY: calculate_rho_elec
51 USE qs_dispersion_nonloc, ONLY: calculate_vxc_nlvdw,&
52 prepare_nlvdw_density,&
53 return_nlvdw_potential
54 USE qs_dispersion_types, ONLY: qs_dispersion_type
55 USE qs_environment_types, ONLY: qs_environment_type
56 USE qs_ks_types, ONLY: get_ks_env,&
57 qs_ks_env_type
58 USE qs_rho_types, ONLY: qs_rho_get,&
59 qs_rho_type
60 USE skala_gpw_functional, ONLY: get_gauxc_section,&
61 native_skala_gapw_composite_direct_ao,&
62 native_skala_gapw_composite_reference,&
63 skala_gpw_eval,&
64 xc_section_uses_native_skala_evaluator
65 USE virial_types, ONLY: virial_type
66 USE xc, ONLY: calc_xc_density,&
67 xc_exc_calc,&
68 xc_vxc_pw_create
69#include "./base/base_uses.f90"
70
71 IMPLICIT NONE
72
73 PRIVATE
74
75 ! *** Public subroutines ***
76 PUBLIC :: qs_vxc_create, qs_xc_density
77
78 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_vxc'
79
80CONTAINS
81
82! **************************************************************************************************
83!> \brief calculates and allocates the xc potential, already reducing it to
84!> the dependence on rho and the one on tau
85!> \param ks_env to get all the needed things
86!> \param rho_struct density for which v_xc is calculated
87!> \param xc_section ...
88!> \param vxc_rho will contain the v_xc part that depend on rho
89!> (if one of the chosen xc functionals has it it is allocated and you
90!> are responsible for it)
91!> \param vxc_tau will contain the kinetic tau part of v_xc
92!> (if one of the chosen xc functionals has it it is allocated and you
93!> are responsible for it)
94!> \param exc ...
95!> \param just_energy if true calculates just the energy, and does not
96!> allocate v_*_rspace
97!> \param edisp ...
98!> \param dispersion_env ...
99!> \param adiabatic_rescale_factor ...
100!> \param pw_env_external external plane wave environment
101!> \param native_skala_atom_force ...
102!> \param qs_env_external parent Quickstep environment for GAPW reconstruction
103!> \param native_gapw_composite_override ...
104!> \param native_skala_defer_to_atom_composite prepare VXC storage while deferring the
105!> native SKALA evaluation to the atom-centered GAPW composite path
106!> \par History
107!> - 05.2002 modified to use the mp_allgather function each pe
108!> computes only part of the grid and this is broadcasted to all
109!> instead of summed.
110!> This scales significantly better (e.g. factor 3 on 12 cpus
111!> 32 H2O) [Joost VdV]
112!> - moved to qs_ks_methods [fawzi]
113!> - sic alterations [Joost VandeVondele]
114!> \author Fawzi Mohamed
115! **************************************************************************************************
116 SUBROUTINE qs_vxc_create(ks_env, rho_struct, xc_section, vxc_rho, vxc_tau, exc, &
117 just_energy, edisp, dispersion_env, adiabatic_rescale_factor, &
118 pw_env_external, native_skala_atom_force, qs_env_external, &
119 native_gapw_composite_override, native_skala_defer_to_atom_composite)
120
121 TYPE(qs_ks_env_type), POINTER :: ks_env
122 TYPE(qs_rho_type), POINTER :: rho_struct
123 TYPE(section_vals_type), POINTER :: xc_section
124 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_rho, vxc_tau
125 REAL(kind=dp), INTENT(out) :: exc
126 LOGICAL, INTENT(in), OPTIONAL :: just_energy
127 REAL(kind=dp), INTENT(out), OPTIONAL :: edisp
128 TYPE(qs_dispersion_type), OPTIONAL, POINTER :: dispersion_env
129 REAL(kind=dp), INTENT(in), OPTIONAL :: adiabatic_rescale_factor
130 TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_external
131 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT), &
132 OPTIONAL :: native_skala_atom_force
133 TYPE(qs_environment_type), OPTIONAL, POINTER :: qs_env_external
134 LOGICAL, INTENT(in), OPTIONAL :: native_gapw_composite_override, &
135 native_skala_defer_to_atom_composite
136
137 CHARACTER(len=*), PARAMETER :: routinen = 'qs_vxc_create'
138
139 INTEGER :: handle, i, ispin, j, k, mspin, myfun, &
140 nelec_spin(2), np, output_unit, vdw
141 LOGICAL :: compute_virial, defer_native_skala_to_atom_composite, do_adiabatic_rescaling, &
142 my_just_energy, native_gapw_composite_direct_ao, native_gapw_composite_reference, &
143 native_grid_diagnostics, native_skala_grid, rho_g_valid, sic_scaling_b_zero, tau_g_valid, &
144 tau_r_valid, uf_grid, vdw_nl
145 REAL(kind=dp) :: composite_hard_integral, composite_rho_max, composite_rho_min, &
146 composite_rho_r_integral, composite_soft_integral, composite_tau_max, composite_tau_min, &
147 composite_tau_r_integral, delta, direct_rho_integral, direct_tau_integral, exc_m, factor, &
148 my_adiabatic_rescale_factor, my_scaling, nelec_s_inv, q_max, rho_composite_hard, &
149 rho_composite_soft, rho_diff_l2, rho_diff_max, rho_ref_l2, tau_composite_hard, &
150 tau_composite_soft, tau_diff_l2, tau_diff_max, tau_ref_l2, volume
151 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: rhonl
152 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: drhonl
153 REAL(kind=dp), DIMENSION(3, 3) :: virial_xc_tmp
154 TYPE(cell_type), POINTER :: cell
155 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
156 TYPE(dft_control_type), POINTER :: dft_control
157 TYPE(mp_para_env_type), POINTER :: para_env
158 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
159 TYPE(pw_c1d_gs_type) :: tmp_g
160 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_direct_g, rho_g, rho_hard_g, &
161 rho_m_gspace, rho_smooth_g, rho_smooth_model_g, rho_struct_g, tau_direct_g, tau_hard_g, &
162 tau_smooth_g, tau_smooth_model_g, tau_struct_g
163 TYPE(pw_c1d_gs_type), POINTER :: rho_nlcc_g, rho_nlcc_g_use, rho_nlcc_g_xc
164 TYPE(pw_env_type), POINTER :: pw_env
165 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool, diagnostics_pw_pool, &
166 vdw_pw_pool, xc_pw_pool
167 TYPE(pw_r3d_rs_type) :: vxcnl
168 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: my_vxc_rho, my_vxc_tau, rho_direct_r, &
169 rho_hard_r, rho_m_rspace, rho_r, rho_smooth_model_r, rho_smooth_r, rho_struct_r, tau, &
170 tau_direct_r, tau_hard_r, tau_smooth_model_r, tau_smooth_r, tau_struct_r
171 TYPE(pw_r3d_rs_type), POINTER :: rho_nlcc, rho_nlcc_use, rho_nlcc_xc, &
172 tmp_pw, weights, weights_use, &
173 weights_xc
174 TYPE(section_vals_type), POINTER :: gauxc_section
175 TYPE(virial_type), POINTER :: virial
176
177 CALL timeset(routinen, handle)
178
179 cpassert(.NOT. ASSOCIATED(vxc_rho))
180 cpassert(.NOT. ASSOCIATED(vxc_tau))
181 NULLIFY (dft_control, pw_env, auxbas_pw_pool, diagnostics_pw_pool, xc_pw_pool, vdw_pw_pool, &
182 cell, my_vxc_rho, &
183 tmp_pw, my_vxc_tau, rho_g, rho_r, tau, rho_m_rspace, &
184 rho_m_gspace, rho_nlcc, rho_nlcc_g, rho_nlcc_g_use, rho_nlcc_g_xc, &
185 rho_nlcc_use, rho_nlcc_xc, rho_struct_r, rho_struct_g, tau_struct_g, tau_struct_r, &
186 weights_use, weights_xc, particle_set, rho_ao_kp, rho_direct_g, rho_direct_r, &
187 rho_hard_g, rho_hard_r, rho_smooth_g, rho_smooth_model_g, rho_smooth_model_r, &
188 rho_smooth_r, tau_direct_g, tau_direct_r, tau_hard_g, tau_hard_r, tau_smooth_g, &
189 tau_smooth_model_g, tau_smooth_model_r, tau_smooth_r, gauxc_section)
190
191 exc = 0.0_dp
192 my_just_energy = .false.
193 IF (PRESENT(just_energy)) my_just_energy = just_energy
194 my_adiabatic_rescale_factor = 1.0_dp
195 do_adiabatic_rescaling = .false.
196 IF (PRESENT(adiabatic_rescale_factor)) THEN
197 my_adiabatic_rescale_factor = adiabatic_rescale_factor
198 do_adiabatic_rescaling = .true.
199 END IF
200
201 CALL get_ks_env(ks_env, &
202 dft_control=dft_control, &
203 pw_env=pw_env, &
204 cell=cell, &
205 para_env=para_env, &
206 particle_set=particle_set, &
207 xcint_weights=weights, &
208 virial=virial, &
209 rho_nlcc=rho_nlcc, &
210 rho_nlcc_g=rho_nlcc_g)
211 rho_nlcc_use => rho_nlcc
212 rho_nlcc_g_use => rho_nlcc_g
213 weights_use => weights
214
215 CALL qs_rho_get(rho_struct, &
216 tau_r_valid=tau_r_valid, &
217 tau_g_valid=tau_g_valid, &
218 rho_g_valid=rho_g_valid, &
219 rho_r=rho_struct_r, &
220 rho_g=rho_struct_g, &
221 tau_g=tau_struct_g, &
222 tau_r=tau_struct_r)
223 rho_smooth_r => rho_struct_r
224 rho_smooth_g => rho_struct_g
225 tau_smooth_r => tau_struct_r
226 tau_smooth_g => tau_struct_g
227
228 compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
229 IF (compute_virial) THEN
230 virial%pv_xc = 0.0_dp
231 END IF
232
233 CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", &
234 i_val=myfun)
235 CALL section_vals_val_get(xc_section, "VDW_POTENTIAL%POTENTIAL_TYPE", &
236 i_val=vdw)
237
238 vdw_nl = (vdw == xc_vdw_fun_nonloc)
239 ! this combination has not been investigated
240 cpassert(.NOT. (do_adiabatic_rescaling .AND. vdw_nl))
241 ! are the necessary inputs available
242 IF (.NOT. (PRESENT(dispersion_env) .AND. PRESENT(edisp))) THEN
243 vdw_nl = .false.
244 END IF
245 IF (PRESENT(edisp)) edisp = 0.0_dp
246 native_skala_grid = xc_section_uses_native_skala_evaluator(xc_section)
247 defer_native_skala_to_atom_composite = .false.
248 IF (PRESENT(native_skala_defer_to_atom_composite)) THEN
249 defer_native_skala_to_atom_composite = native_skala_defer_to_atom_composite
250 END IF
251 cpassert(.NOT. defer_native_skala_to_atom_composite .OR. native_skala_grid)
252 native_gapw_composite_reference = native_skala_grid .AND. &
253 native_skala_gapw_composite_reference(xc_section) .AND. &
254 (dft_control%qs_control%gapw .OR. &
255 dft_control%qs_control%gapw_xc)
256 IF (PRESENT(native_gapw_composite_override)) THEN
257 native_gapw_composite_reference = native_skala_grid .AND. &
258 native_gapw_composite_override .AND. &
259 (dft_control%qs_control%gapw .OR. &
260 dft_control%qs_control%gapw_xc)
261 END IF
262 native_gapw_composite_direct_ao = native_gapw_composite_reference .AND. &
263 native_skala_gapw_composite_direct_ao(xc_section)
264 native_grid_diagnostics = .false.
265 IF (native_gapw_composite_reference) THEN
266 gauxc_section => get_gauxc_section(xc_section)
267 cpassert(ASSOCIATED(gauxc_section))
268 CALL section_vals_val_get(gauxc_section, "NATIVE_GRID_DIAGNOSTICS", &
269 l_val=native_grid_diagnostics)
270 END IF
271 IF (native_gapw_composite_reference) NULLIFY (weights_use)
272
273 IF (myfun /= xc_none .OR. vdw_nl) THEN
274
275 ! test if the real space density is available
276 cpassert(ASSOCIATED(rho_struct))
277 IF (dft_control%nspins /= 1 .AND. dft_control%nspins /= 2) THEN
278 cpabort("nspins must be 1 or 2")
279 END IF
280 mspin = SIZE(rho_struct_r)
281 IF (dft_control%nspins == 2 .AND. mspin == 1) THEN
282 cpabort("Spin count mismatch")
283 END IF
284
285 ! there are some options related to SIC here.
286 ! Normal DFT computes E(rho_alpha,rho_beta) (or its variant E(2*rho_alpha) for non-LSD)
287 ! SIC can E(rho_alpha,rho_beta)-b*(E(rho_alpha,rho_beta)-E(rho_beta,rho_beta))
288 ! or compute E(rho_alpha,rho_beta)-b*E(rho_alpha-rho_beta,0)
289
290 ! my_scaling is the scaling needed of the standard E(rho_alpha,rho_beta) term
291 my_scaling = 1.0_dp
292 SELECT CASE (dft_control%sic_method_id)
293 CASE (sic_none)
294 ! all fine
295 CASE (sic_mauri_spz, sic_ad)
296 ! no idea yet what to do here in that case
297 cpassert(.NOT. tau_r_valid)
298 CASE (sic_mauri_us)
299 my_scaling = 1.0_dp - dft_control%sic_scaling_b
300 ! no idea yet what to do here in that case
301 cpassert(.NOT. tau_r_valid)
302 CASE (sic_eo)
303 ! NOTHING TO BE DONE
304 CASE DEFAULT
305 ! this case has not yet been treated here
306 cpabort("NYI")
307 END SELECT
308
309 IF (dft_control%sic_scaling_b == 0.0_dp) THEN
310 sic_scaling_b_zero = .true.
311 ELSE
312 sic_scaling_b_zero = .false.
313 END IF
314
315 IF (PRESENT(pw_env_external)) THEN
316 pw_env => pw_env_external
317 END IF
318 CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool, auxbas_pw_pool=auxbas_pw_pool)
319
320 IF (native_gapw_composite_reference) THEN
321 cpassert(tau_r_valid)
322 cpassert(PRESENT(qs_env_external))
323 ALLOCATE (rho_hard_r(mspin), rho_hard_g(mspin), tau_hard_r(mspin), tau_hard_g(mspin))
324 DO ispin = 1, mspin
325 CALL auxbas_pw_pool%create_pw(rho_hard_r(ispin))
326 CALL auxbas_pw_pool%create_pw(rho_hard_g(ispin))
327 CALL auxbas_pw_pool%create_pw(tau_hard_r(ispin))
328 CALL auxbas_pw_pool%create_pw(tau_hard_g(ispin))
329 END DO
330 IF (native_gapw_composite_direct_ao) THEN
331 CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
332 cpassert(ASSOCIATED(rho_ao_kp))
333 cpassert(SIZE(rho_ao_kp, 2) == 1)
334 DO ispin = 1, mspin
335 CALL calculate_rho_elec( &
336 matrix_p=rho_ao_kp(ispin, 1)%matrix, rho=rho_hard_r(ispin), &
337 rho_gspace=rho_hard_g(ispin), ks_env=ks_env)
338 CALL calculate_rho_elec( &
339 matrix_p=rho_ao_kp(ispin, 1)%matrix, rho=tau_hard_r(ispin), &
340 rho_gspace=tau_hard_g(ispin), ks_env=ks_env, compute_tau=.true.)
341 END DO
342 ELSE
343 CALL get_pw_grid_info(rho_hard_g(1)%pw_grid, vol=volume)
344 q_max = sqrt(maxval(rho_hard_g(1)%pw_grid%gsq))
345 CALL calculate_rhotot_elec_gspace( &
346 qs_env=qs_env_external, auxbas_pw_pool=auxbas_pw_pool, &
347 rhotot_elec_gspace=rho_hard_g(1), q_max=q_max, &
348 rho_hard=composite_hard_integral, rho_soft=composite_soft_integral, &
349 rho_source=rho_struct, allow_nonorthorhombic=.true.)
350 rho_composite_hard = composite_hard_integral
351 rho_composite_soft = composite_soft_integral
352 CALL calculate_rhotot_elec_gspace( &
353 qs_env=qs_env_external, auxbas_pw_pool=auxbas_pw_pool, &
354 rhotot_elec_gspace=tau_hard_g(1), q_max=q_max, &
355 rho_hard=composite_hard_integral, rho_soft=composite_soft_integral, &
356 compute_tau=.true., rho_source=rho_struct, allow_nonorthorhombic=.true.)
357 tau_composite_hard = composite_hard_integral
358 tau_composite_soft = composite_soft_integral
359 IF (para_env%mepos == 0) THEN
360 output_unit = cp_logger_get_default_io_unit()
361 IF (output_unit > 0) THEN
362 WRITE (unit=output_unit, fmt="(/,T2,A,2(1X,ES19.11))") &
363 "SKALA_GPW| Composite rho hard and soft integrals", &
364 rho_composite_hard, rho_composite_soft
365 WRITE (unit=output_unit, fmt="(T2,A,2(1X,ES19.11))") &
366 "SKALA_GPW| Composite tau hard and soft integrals", &
367 tau_composite_hard, tau_composite_soft
368 END IF
369 END IF
370 IF (mspin == 1) THEN
371 CALL pw_scale(rho_hard_g(1), 1.0_dp/volume)
372 CALL pw_scale(tau_hard_g(1), 1.0_dp/volume)
373 ELSE
374 CALL calculate_rhotot_elec_gspace( &
375 qs_env=qs_env_external, auxbas_pw_pool=auxbas_pw_pool, &
376 rhotot_elec_gspace=rho_hard_g(2), q_max=q_max, &
377 rho_hard=composite_hard_integral, rho_soft=composite_soft_integral, &
378 fsign=-1.0_dp, rho_source=rho_struct, allow_nonorthorhombic=.true.)
379 CALL calculate_rhotot_elec_gspace( &
380 qs_env=qs_env_external, auxbas_pw_pool=auxbas_pw_pool, &
381 rhotot_elec_gspace=tau_hard_g(2), q_max=q_max, &
382 rho_hard=composite_hard_integral, rho_soft=composite_soft_integral, &
383 fsign=-1.0_dp, compute_tau=.true., rho_source=rho_struct, &
384 allow_nonorthorhombic=.true.)
385 CALL pw_scale(rho_hard_g(1), 0.5_dp/volume)
386 CALL pw_scale(rho_hard_g(2), 0.5_dp/volume)
387 CALL auxbas_pw_pool%create_pw(tmp_g)
388 CALL pw_copy(rho_hard_g(1), tmp_g)
389 CALL pw_axpy(rho_hard_g(2), rho_hard_g(1), 1.0_dp)
390 CALL pw_axpy(rho_hard_g(2), tmp_g, -1.0_dp)
391 CALL pw_copy(tmp_g, rho_hard_g(2))
392 CALL pw_scale(tau_hard_g(1), 0.5_dp/volume)
393 CALL pw_scale(tau_hard_g(2), 0.5_dp/volume)
394 CALL pw_copy(tau_hard_g(1), tmp_g)
395 CALL pw_axpy(tau_hard_g(2), tau_hard_g(1), 1.0_dp)
396 CALL pw_axpy(tau_hard_g(2), tmp_g, -1.0_dp)
397 CALL pw_copy(tmp_g, tau_hard_g(2))
398 CALL auxbas_pw_pool%give_back_pw(tmp_g)
399 END IF
400 END IF
401 DO ispin = 1, mspin
402 CALL pw_transfer(rho_hard_g(ispin), rho_hard_r(ispin))
403 CALL pw_transfer(tau_hard_g(ispin), tau_hard_r(ispin))
404 END DO
405 IF (native_grid_diagnostics .AND. .NOT. native_gapw_composite_direct_ao) THEN
406 CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
407 cpassert(ASSOCIATED(rho_ao_kp))
408 cpassert(SIZE(rho_ao_kp, 2) == 1)
409 ALLOCATE (rho_direct_r(mspin), rho_direct_g(mspin), &
410 tau_direct_r(mspin), tau_direct_g(mspin))
411 rho_diff_l2 = 0.0_dp
412 rho_diff_max = 0.0_dp
413 rho_ref_l2 = 0.0_dp
414 tau_diff_l2 = 0.0_dp
415 tau_diff_max = 0.0_dp
416 tau_ref_l2 = 0.0_dp
417 direct_rho_integral = 0.0_dp
418 direct_tau_integral = 0.0_dp
419 DO ispin = 1, mspin
420 CALL auxbas_pw_pool%create_pw(rho_direct_r(ispin))
421 CALL auxbas_pw_pool%create_pw(rho_direct_g(ispin))
422 CALL auxbas_pw_pool%create_pw(tau_direct_r(ispin))
423 CALL auxbas_pw_pool%create_pw(tau_direct_g(ispin))
424 CALL calculate_rho_elec( &
425 matrix_p=rho_ao_kp(ispin, 1)%matrix, rho=rho_direct_r(ispin), &
426 rho_gspace=rho_direct_g(ispin), ks_env=ks_env)
427 CALL calculate_rho_elec( &
428 matrix_p=rho_ao_kp(ispin, 1)%matrix, rho=tau_direct_r(ispin), &
429 rho_gspace=tau_direct_g(ispin), ks_env=ks_env, compute_tau=.true.)
430 direct_rho_integral = direct_rho_integral + &
431 pw_integrate_function(rho_direct_r(ispin))
432 direct_tau_integral = direct_tau_integral + &
433 pw_integrate_function(tau_direct_r(ispin))
434 DO k = lbound(rho_direct_r(ispin)%array, 3), ubound(rho_direct_r(ispin)%array, 3)
435 DO j = lbound(rho_direct_r(ispin)%array, 2), ubound(rho_direct_r(ispin)%array, 2)
436 DO i = lbound(rho_direct_r(ispin)%array, 1), ubound(rho_direct_r(ispin)%array, 1)
437 delta = rho_hard_r(ispin)%array(i, j, k) - &
438 rho_direct_r(ispin)%array(i, j, k)
439 rho_diff_l2 = rho_diff_l2 + delta*delta
440 rho_diff_max = max(rho_diff_max, abs(delta))
441 rho_ref_l2 = rho_ref_l2 + rho_direct_r(ispin)%array(i, j, k)**2
442 delta = tau_hard_r(ispin)%array(i, j, k) - &
443 tau_direct_r(ispin)%array(i, j, k)
444 tau_diff_l2 = tau_diff_l2 + delta*delta
445 tau_diff_max = max(tau_diff_max, abs(delta))
446 tau_ref_l2 = tau_ref_l2 + tau_direct_r(ispin)%array(i, j, k)**2
447 END DO
448 END DO
449 END DO
450 END DO
451 CALL para_env%sum(rho_diff_l2)
452 CALL para_env%sum(rho_ref_l2)
453 CALL para_env%sum(tau_diff_l2)
454 CALL para_env%sum(tau_ref_l2)
455 CALL para_env%max(rho_diff_max)
456 CALL para_env%max(tau_diff_max)
457 rho_diff_l2 = sqrt(rho_diff_l2/max(rho_ref_l2, tiny(1.0_dp)))
458 tau_diff_l2 = sqrt(tau_diff_l2/max(tau_ref_l2, tiny(1.0_dp)))
459 IF (para_env%mepos == 0) THEN
460 output_unit = cp_logger_get_default_io_unit()
461 IF (output_unit > 0) THEN
462 WRITE (unit=output_unit, fmt="(T2,A,2(1X,ES19.11))") &
463 "SKALA_GPW| Direct AO rho and tau integrals", &
464 direct_rho_integral, direct_tau_integral
465 WRITE (unit=output_unit, fmt="(T2,A,2(1X,ES19.11))") &
466 "SKALA_GPW| Composite/direct relative L2 rho and tau", &
467 rho_diff_l2, tau_diff_l2
468 WRITE (unit=output_unit, fmt="(T2,A,2(1X,ES19.11))") &
469 "SKALA_GPW| Composite/direct Linf rho and tau", &
470 rho_diff_max, tau_diff_max
471 END IF
472 END IF
473 DO ispin = 1, mspin
474 CALL auxbas_pw_pool%give_back_pw(rho_direct_r(ispin))
475 CALL auxbas_pw_pool%give_back_pw(rho_direct_g(ispin))
476 CALL auxbas_pw_pool%give_back_pw(tau_direct_r(ispin))
477 CALL auxbas_pw_pool%give_back_pw(tau_direct_g(ispin))
478 END DO
479 DEALLOCATE (rho_direct_r, rho_direct_g, tau_direct_r, tau_direct_g)
480 END IF
481 composite_rho_r_integral = 0.0_dp
482 composite_tau_r_integral = 0.0_dp
483 composite_rho_min = huge(1.0_dp)
484 composite_rho_max = -huge(1.0_dp)
485 composite_tau_min = huge(1.0_dp)
486 composite_tau_max = -huge(1.0_dp)
487 DO ispin = 1, mspin
488 composite_rho_r_integral = composite_rho_r_integral + &
489 pw_integrate_function(rho_hard_r(ispin))
490 composite_tau_r_integral = composite_tau_r_integral + &
491 pw_integrate_function(tau_hard_r(ispin))
492 composite_rho_min = min(composite_rho_min, minval(rho_hard_r(ispin)%array))
493 composite_rho_max = max(composite_rho_max, maxval(rho_hard_r(ispin)%array))
494 composite_tau_min = min(composite_tau_min, minval(tau_hard_r(ispin)%array))
495 composite_tau_max = max(composite_tau_max, maxval(tau_hard_r(ispin)%array))
496 END DO
497 CALL para_env%min(composite_rho_min)
498 CALL para_env%max(composite_rho_max)
499 CALL para_env%min(composite_tau_min)
500 CALL para_env%max(composite_tau_max)
501 IF (para_env%mepos == 0) THEN
502 output_unit = cp_logger_get_default_io_unit()
503 IF (output_unit > 0) THEN
504 WRITE (unit=output_unit, fmt="(T2,A,1X,ES19.11)") &
505 "SKALA_GPW| Composite real-grid rho integral", composite_rho_r_integral
506 WRITE (unit=output_unit, fmt="(T2,A,1X,ES19.11)") &
507 "SKALA_GPW| Composite real-grid tau integral", composite_tau_r_integral
508 WRITE (unit=output_unit, fmt="(T2,A,2(1X,ES19.11))") &
509 "SKALA_GPW| Composite real-grid rho min and max", &
510 composite_rho_min, composite_rho_max
511 WRITE (unit=output_unit, fmt="(T2,A,2(1X,ES19.11))") &
512 "SKALA_GPW| Composite real-grid tau min and max", &
513 composite_tau_min, composite_tau_max
514 END IF
515 END IF
516 rho_struct_r => rho_hard_r
517 rho_struct_g => rho_hard_g
518 tau_struct_r => tau_hard_r
519 tau_struct_g => tau_hard_g
520 rho_g_valid = .true.
521 tau_g_valid = .true.
522 END IF
523
524 uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid, xc_pw_pool%pw_grid)
525
526 IF (.NOT. uf_grid) THEN
527 rho_r => rho_struct_r
528
529 IF (tau_r_valid) THEN
530 tau => tau_struct_r
531 END IF
532
533 ! for gradient corrected functional the density in g space might
534 ! be useful so if we have it, we pass it in
535 IF (rho_g_valid) THEN
536 rho_g => rho_struct_g
537 END IF
538 ELSE
539 cpassert(rho_g_valid)
540 ALLOCATE (rho_r(mspin))
541 ALLOCATE (rho_g(mspin))
542 DO ispin = 1, mspin
543 CALL xc_pw_pool%create_pw(rho_g(ispin))
544 CALL pw_transfer(rho_struct_g(ispin), rho_g(ispin))
545 END DO
546 DO ispin = 1, mspin
547 CALL xc_pw_pool%create_pw(rho_r(ispin))
548 CALL pw_transfer(rho_g(ispin), rho_r(ispin))
549 END DO
550 IF (tau_r_valid) THEN
551 ALLOCATE (tau(mspin))
552 DO ispin = 1, mspin
553 CALL xc_pw_pool%create_pw(tau(ispin))
554 block
555 TYPE(pw_c1d_gs_type) :: tau_g_aux, tau_g_xc
556 CALL xc_pw_pool%create_pw(tau_g_xc)
557 IF (tau_g_valid) THEN
558 CALL pw_transfer(tau_struct_g(ispin), tau_g_xc)
559 ELSE
560 CALL auxbas_pw_pool%create_pw(tau_g_aux)
561 CALL pw_transfer(tau_struct_r(ispin), tau_g_aux)
562 CALL pw_transfer(tau_g_aux, tau_g_xc)
563 CALL auxbas_pw_pool%give_back_pw(tau_g_aux)
564 END IF
565 CALL pw_transfer(tau_g_xc, tau(ispin))
566 CALL xc_pw_pool%give_back_pw(tau_g_xc)
567 END block
568 END DO
569 END IF
570 IF (ASSOCIATED(weights) .AND. .NOT. native_gapw_composite_reference) THEN
571 ALLOCATE (weights_xc)
572 CALL xc_pw_pool%create_pw(weights_xc)
573 block
574 TYPE(pw_c1d_gs_type) :: weights_g_aux, weights_g_xc
575 CALL auxbas_pw_pool%create_pw(weights_g_aux)
576 CALL xc_pw_pool%create_pw(weights_g_xc)
577 CALL pw_transfer(weights, weights_g_aux)
578 CALL pw_transfer(weights_g_aux, weights_g_xc)
579 CALL pw_transfer(weights_g_xc, weights_xc)
580 CALL xc_pw_pool%give_back_pw(weights_g_xc)
581 CALL auxbas_pw_pool%give_back_pw(weights_g_aux)
582 END block
583 weights_use => weights_xc
584 END IF
585 IF (ASSOCIATED(rho_nlcc)) THEN
586 cpassert(ASSOCIATED(rho_nlcc_g))
587 ALLOCATE (rho_nlcc_g_xc, rho_nlcc_xc)
588 CALL xc_pw_pool%create_pw(rho_nlcc_g_xc)
589 CALL xc_pw_pool%create_pw(rho_nlcc_xc)
590 CALL pw_transfer(rho_nlcc_g, rho_nlcc_g_xc)
591 CALL pw_transfer(rho_nlcc_g_xc, rho_nlcc_xc)
592 rho_nlcc_use => rho_nlcc_xc
593 rho_nlcc_g_use => rho_nlcc_g_xc
594 END IF
595 END IF
596
597 IF (native_gapw_composite_reference) THEN
598 composite_rho_r_integral = 0.0_dp
599 composite_tau_r_integral = 0.0_dp
600 DO ispin = 1, mspin
601 composite_rho_r_integral = composite_rho_r_integral + &
602 pw_integrate_function(rho_r(ispin))
603 composite_tau_r_integral = composite_tau_r_integral + &
604 pw_integrate_function(tau(ispin))
605 END DO
606 IF (para_env%mepos == 0) THEN
607 output_unit = cp_logger_get_default_io_unit()
608 IF (output_unit > 0) THEN
609 WRITE (unit=output_unit, fmt="(T2,A,1X,ES19.11)") &
610 "SKALA_GPW| Composite XC-grid rho integral", composite_rho_r_integral
611 WRITE (unit=output_unit, fmt="(T2,A,1X,ES19.11)") &
612 "SKALA_GPW| Composite XC-grid tau integral", composite_tau_r_integral
613 END IF
614 END IF
615 END IF
616
617 ! NLCC augments rho and grad(rho); tau remains valence-only.
618 IF (ASSOCIATED(rho_nlcc_use)) THEN
619 factor = 1.0_dp
620 DO ispin = 1, mspin
621 CALL pw_axpy(rho_nlcc_use, rho_r(ispin), factor)
622 CALL pw_axpy(rho_nlcc_g_use, rho_g(ispin), factor)
623 END DO
624 END IF
625
626 !
627 ! here the rho_r, rho_g, tau is what it should be
628 ! we get back the right my_vxc_rho and my_vxc_tau as required
629 !
630 IF (defer_native_skala_to_atom_composite) THEN
631 ! The atom-centered composite path supplies the complete SKALA energy and
632 ! VXC below qs_vxc_create. Keep only storage for that VXC (and any additive
633 ! nonlocal-dispersion potential) instead of evaluating the smooth grid here.
634 ALLOCATE (my_vxc_rho(mspin), my_vxc_tau(mspin))
635 DO ispin = 1, mspin
636 CALL xc_pw_pool%create_pw(my_vxc_rho(ispin))
637 CALL xc_pw_pool%create_pw(my_vxc_tau(ispin))
638 CALL pw_zero(my_vxc_rho(ispin))
639 CALL pw_zero(my_vxc_tau(ispin))
640 END DO
641 ELSE IF (native_skala_grid) THEN
642 CALL skala_gpw_eval(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, exc=exc, &
643 rho_r=rho_r, rho_g=rho_g, tau=tau, xc_section=xc_section, &
644 weights=weights_use, pw_pool=xc_pw_pool, &
645 particle_set=particle_set, cell=cell, &
646 compute_virial=compute_virial, virial_xc=virial%pv_xc, &
647 just_energy=my_just_energy, atom_force=native_skala_atom_force)
648 ELSE IF (my_just_energy) THEN
649 exc = xc_exc_calc(rho_r=rho_r, tau=tau, &
650 rho_g=rho_g, xc_section=xc_section, &
651 weights=weights_use, pw_pool=xc_pw_pool)
652
653 ELSE
654 CALL xc_vxc_pw_create(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_r, &
655 rho_g=rho_g, tau=tau, exc=exc, &
656 xc_section=xc_section, &
657 weights=weights_use, pw_pool=xc_pw_pool, &
658 compute_virial=compute_virial, &
659 virial_xc=virial%pv_xc)
660 END IF
661
662 ! The composite density is assembled through an R-to-G-to-R projection. Apply the
663 ! discrete adjoint of that projection before integrating VXC with the smooth AO density.
664 IF (native_gapw_composite_reference .AND. .NOT. my_just_energy) THEN
665 IF (uf_grid) THEN
666 diagnostics_pw_pool => xc_pw_pool
667 ELSE
668 diagnostics_pw_pool => auxbas_pw_pool
669 END IF
670 CALL diagnostics_pw_pool%create_pw(tmp_g)
671 DO ispin = 1, mspin
672 CALL pw_transfer(my_vxc_rho(ispin), tmp_g)
673 CALL pw_transfer(tmp_g, my_vxc_rho(ispin))
674 CALL pw_transfer(my_vxc_tau(ispin), tmp_g)
675 CALL pw_transfer(tmp_g, my_vxc_tau(ispin))
676 END DO
677 CALL diagnostics_pw_pool%give_back_pw(tmp_g)
678 NULLIFY (diagnostics_pw_pool)
679 END IF
680
681 IF (native_gapw_composite_reference .AND. native_grid_diagnostics .AND. &
682 .NOT. my_just_energy) THEN
683 IF (uf_grid) THEN
684 diagnostics_pw_pool => xc_pw_pool
685 ELSE
686 diagnostics_pw_pool => auxbas_pw_pool
687 END IF
688 ALLOCATE (rho_smooth_model_g(mspin), rho_smooth_model_r(mspin), &
689 tau_smooth_model_g(mspin), tau_smooth_model_r(mspin))
690 DO ispin = 1, mspin
691 CALL diagnostics_pw_pool%create_pw(rho_smooth_model_g(ispin))
692 CALL diagnostics_pw_pool%create_pw(rho_smooth_model_r(ispin))
693 CALL diagnostics_pw_pool%create_pw(tau_smooth_model_g(ispin))
694 CALL diagnostics_pw_pool%create_pw(tau_smooth_model_r(ispin))
695 IF (uf_grid) THEN
696 IF (rho_g_valid) THEN
697 CALL pw_transfer(rho_smooth_g(ispin), rho_smooth_model_g(ispin))
698 ELSE
699 CALL auxbas_pw_pool%create_pw(tmp_g)
700 CALL pw_transfer(rho_smooth_r(ispin), tmp_g)
701 CALL pw_transfer(tmp_g, rho_smooth_model_g(ispin))
702 CALL auxbas_pw_pool%give_back_pw(tmp_g)
703 END IF
704 IF (tau_g_valid) THEN
705 CALL pw_transfer(tau_smooth_g(ispin), tau_smooth_model_g(ispin))
706 ELSE
707 CALL auxbas_pw_pool%create_pw(tmp_g)
708 CALL pw_transfer(tau_smooth_r(ispin), tmp_g)
709 CALL pw_transfer(tmp_g, tau_smooth_model_g(ispin))
710 CALL auxbas_pw_pool%give_back_pw(tmp_g)
711 END IF
712 ELSE
713 IF (rho_g_valid) THEN
714 CALL pw_copy(rho_smooth_g(ispin), rho_smooth_model_g(ispin))
715 ELSE
716 CALL pw_transfer(rho_smooth_r(ispin), rho_smooth_model_g(ispin))
717 END IF
718 IF (tau_g_valid) THEN
719 CALL pw_copy(tau_smooth_g(ispin), tau_smooth_model_g(ispin))
720 ELSE
721 CALL pw_transfer(tau_smooth_r(ispin), tau_smooth_model_g(ispin))
722 END IF
723 END IF
724 CALL pw_transfer(rho_smooth_model_g(ispin), rho_smooth_model_r(ispin))
725 CALL pw_transfer(tau_smooth_model_g(ispin), tau_smooth_model_r(ispin))
726 IF (ASSOCIATED(rho_nlcc_use)) THEN
727 CALL pw_axpy(rho_nlcc_use, rho_smooth_model_r(ispin), 1.0_dp)
728 CALL pw_axpy(rho_nlcc_g_use, rho_smooth_model_g(ispin), 1.0_dp)
729 END IF
730 END DO
731 CALL diagnose_gapw_composite_direction( &
732 rho_r, rho_g, tau, rho_smooth_model_r, rho_smooth_model_g, &
733 tau_smooth_model_r, my_vxc_rho, my_vxc_tau, xc_section, weights_use, &
734 diagnostics_pw_pool, particle_set, cell, para_env)
735 DO ispin = 1, mspin
736 CALL diagnostics_pw_pool%give_back_pw(rho_smooth_model_g(ispin))
737 CALL diagnostics_pw_pool%give_back_pw(rho_smooth_model_r(ispin))
738 CALL diagnostics_pw_pool%give_back_pw(tau_smooth_model_g(ispin))
739 CALL diagnostics_pw_pool%give_back_pw(tau_smooth_model_r(ispin))
740 END DO
741 DEALLOCATE (rho_smooth_model_g, rho_smooth_model_r, &
742 tau_smooth_model_g, tau_smooth_model_r)
743 NULLIFY (diagnostics_pw_pool)
744 END IF
745
746 ! calclulate non-local vdW functional
747 ! only if this XC_SECTION has it
748 ! if yes, we use the dispersion_env from ks_env
749 ! this is dangerous, as it assumes a special connection xc_section -> qs_env
750 IF (vdw_nl) THEN
751 CALL get_ks_env(ks_env=ks_env, para_env=para_env)
752 ! no SIC functionals allowed
753 cpassert(dft_control%sic_method_id == sic_none)
754 !
755 CALL pw_env_get(pw_env, vdw_pw_pool=vdw_pw_pool)
756 np = vdw_pw_pool%pw_grid%ngpts_local
757 ALLOCATE (rhonl(np), drhonl(np, 3))
758 CALL prepare_nlvdw_density(rho_r, rhonl, drhonl, vdw_pw_pool, xc_pw_pool)
759 IF (my_just_energy) THEN
760 CALL calculate_vxc_nlvdw(vxcnl, rhonl, drhonl, edisp, dispersion_env, &
761 my_just_energy, vdw_pw_pool, para_env)
762 ELSE
763 CALL vdw_pw_pool%create_pw(vxcnl)
764 CALL calculate_vxc_nlvdw(vxcnl, rhonl, drhonl, edisp, dispersion_env, &
765 my_just_energy, vdw_pw_pool, para_env, virial=virial)
766 CALL return_nlvdw_potential(my_vxc_rho, vxcnl, vdw_pw_pool, xc_pw_pool)
767 CALL vdw_pw_pool%give_back_pw(vxcnl)
768 END IF
769 DEALLOCATE (rhonl, drhonl)
770 END IF
771
772 ! remove the nlcc densities (keep stuff in original state)
773 IF (ASSOCIATED(rho_nlcc_use)) THEN
774 factor = -1.0_dp
775 DO ispin = 1, mspin
776 CALL pw_axpy(rho_nlcc_use, rho_r(ispin), factor)
777 CALL pw_axpy(rho_nlcc_g_use, rho_g(ispin), factor)
778 END DO
779 END IF
780
781 !! Apply rescaling to the potential if requested
782 IF (.NOT. my_just_energy) THEN
783 IF (do_adiabatic_rescaling) THEN
784 IF (ASSOCIATED(my_vxc_rho)) THEN
785 DO ispin = 1, SIZE(my_vxc_rho)
786 CALL pw_scale(my_vxc_rho(ispin), my_adiabatic_rescale_factor)
787 END DO
788 END IF
789 END IF
790 END IF
791
792 IF (my_scaling /= 1.0_dp) THEN
793 exc = exc*my_scaling
794 IF (ASSOCIATED(my_vxc_rho)) THEN
795 DO ispin = 1, SIZE(my_vxc_rho)
796 CALL pw_scale(my_vxc_rho(ispin), my_scaling)
797 END DO
798 END IF
799 IF (ASSOCIATED(my_vxc_tau)) THEN
800 DO ispin = 1, SIZE(my_vxc_tau)
801 CALL pw_scale(my_vxc_tau(ispin), my_scaling)
802 END DO
803 END IF
804 END IF
805
806 ! we have pw data for the xc, qs_ks requests coeff structure, here we transfer
807 ! pw -> coeff
808 IF (ASSOCIATED(my_vxc_rho)) THEN
809 vxc_rho => my_vxc_rho
810 NULLIFY (my_vxc_rho)
811 END IF
812 IF (ASSOCIATED(my_vxc_tau)) THEN
813 vxc_tau => my_vxc_tau
814 NULLIFY (my_vxc_tau)
815 END IF
816 IF (uf_grid) THEN
817 DO ispin = 1, SIZE(rho_r)
818 CALL xc_pw_pool%give_back_pw(rho_r(ispin))
819 END DO
820 DEALLOCATE (rho_r)
821 IF (ASSOCIATED(rho_g)) THEN
822 DO ispin = 1, SIZE(rho_g)
823 CALL xc_pw_pool%give_back_pw(rho_g(ispin))
824 END DO
825 DEALLOCATE (rho_g)
826 END IF
827 END IF
828
829 ! compute again the xc but now for Exc(m,o) and the opposite sign
830 IF (dft_control%sic_method_id == sic_mauri_spz .AND. .NOT. sic_scaling_b_zero) THEN
831 ALLOCATE (rho_m_rspace(2), rho_m_gspace(2))
832 CALL xc_pw_pool%create_pw(rho_m_gspace(1))
833 CALL xc_pw_pool%create_pw(rho_m_rspace(1))
834 CALL pw_copy(rho_struct_r(1), rho_m_rspace(1))
835 CALL pw_axpy(rho_struct_r(2), rho_m_rspace(1), alpha=-1._dp)
836 CALL pw_copy(rho_struct_g(1), rho_m_gspace(1))
837 CALL pw_axpy(rho_struct_g(2), rho_m_gspace(1), alpha=-1._dp)
838 ! bit sad, these will be just zero...
839 CALL xc_pw_pool%create_pw(rho_m_gspace(2))
840 CALL xc_pw_pool%create_pw(rho_m_rspace(2))
841 CALL pw_zero(rho_m_rspace(2))
842 CALL pw_zero(rho_m_gspace(2))
843
844 IF (my_just_energy) THEN
845 exc_m = xc_exc_calc(rho_r=rho_m_rspace, tau=tau, &
846 rho_g=rho_m_gspace, xc_section=xc_section, &
847 weights=weights_use, pw_pool=xc_pw_pool)
848 ELSE
849 ! virial untested
850 cpassert(.NOT. compute_virial)
851 CALL xc_vxc_pw_create(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_m_rspace, &
852 rho_g=rho_m_gspace, tau=tau, exc=exc_m, &
853 xc_section=xc_section, &
854 weights=weights_use, pw_pool=xc_pw_pool, &
855 compute_virial=.false., &
856 virial_xc=virial_xc_tmp)
857 END IF
858
859 exc = exc - dft_control%sic_scaling_b*exc_m
860
861 ! and take care of the potential only vxc_rho is taken into account
862 IF (.NOT. my_just_energy) THEN
863 CALL pw_axpy(my_vxc_rho(1), vxc_rho(1), -dft_control%sic_scaling_b)
864 CALL pw_axpy(my_vxc_rho(1), vxc_rho(2), dft_control%sic_scaling_b)
865 CALL my_vxc_rho(1)%release()
866 CALL my_vxc_rho(2)%release()
867 DEALLOCATE (my_vxc_rho)
868 END IF
869
870 DO ispin = 1, 2
871 CALL xc_pw_pool%give_back_pw(rho_m_rspace(ispin))
872 CALL xc_pw_pool%give_back_pw(rho_m_gspace(ispin))
873 END DO
874 DEALLOCATE (rho_m_rspace)
875 DEALLOCATE (rho_m_gspace)
876
877 END IF
878
879 ! now we have - sum_s N_s * Exc(rho_s/N_s,0)
880 IF (dft_control%sic_method_id == sic_ad .AND. .NOT. sic_scaling_b_zero) THEN
881
882 ! find out how many elecs we have
883 CALL get_ks_env(ks_env, nelectron_spin=nelec_spin)
884
885 ALLOCATE (rho_m_rspace(2), rho_m_gspace(2))
886 DO ispin = 1, 2
887 CALL xc_pw_pool%create_pw(rho_m_gspace(ispin))
888 CALL xc_pw_pool%create_pw(rho_m_rspace(ispin))
889 END DO
890
891 DO ispin = 1, 2
892 IF (nelec_spin(ispin) > 0.0_dp) THEN
893 nelec_s_inv = 1.0_dp/nelec_spin(ispin)
894 ELSE
895 ! does it matter if there are no electrons with this spin (H) ?
896 nelec_s_inv = 0.0_dp
897 END IF
898 CALL pw_copy(rho_struct_r(ispin), rho_m_rspace(1))
899 CALL pw_copy(rho_struct_g(ispin), rho_m_gspace(1))
900 CALL pw_scale(rho_m_rspace(1), nelec_s_inv)
901 CALL pw_scale(rho_m_gspace(1), nelec_s_inv)
902 CALL pw_zero(rho_m_rspace(2))
903 CALL pw_zero(rho_m_gspace(2))
904
905 IF (my_just_energy) THEN
906 exc_m = xc_exc_calc(rho_r=rho_m_rspace, tau=tau, &
907 rho_g=rho_m_gspace, xc_section=xc_section, &
908 weights=weights_use, pw_pool=xc_pw_pool)
909 ELSE
910 ! virial untested
911 cpassert(.NOT. compute_virial)
912 CALL xc_vxc_pw_create(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_m_rspace, &
913 rho_g=rho_m_gspace, tau=tau, exc=exc_m, &
914 xc_section=xc_section, &
915 weights=weights_use, pw_pool=xc_pw_pool, &
916 compute_virial=.false., &
917 virial_xc=virial_xc_tmp)
918 END IF
919
920 exc = exc - dft_control%sic_scaling_b*nelec_spin(ispin)*exc_m
921
922 ! and take care of the potential only vxc_rho is taken into account
923 IF (.NOT. my_just_energy) THEN
924 CALL pw_axpy(my_vxc_rho(1), vxc_rho(ispin), -dft_control%sic_scaling_b)
925 CALL my_vxc_rho(1)%release()
926 CALL my_vxc_rho(2)%release()
927 DEALLOCATE (my_vxc_rho)
928 END IF
929 END DO
930
931 DO ispin = 1, 2
932 CALL xc_pw_pool%give_back_pw(rho_m_rspace(ispin))
933 CALL xc_pw_pool%give_back_pw(rho_m_gspace(ispin))
934 END DO
935 DEALLOCATE (rho_m_rspace)
936 DEALLOCATE (rho_m_gspace)
937
938 END IF
939
940 ! compute again the xc but now for Exc(n_down,n_down)
941 IF (dft_control%sic_method_id == sic_mauri_us .AND. .NOT. sic_scaling_b_zero) THEN
942 ALLOCATE (rho_r(2))
943 rho_r(1) = rho_struct_r(2)
944 rho_r(2) = rho_struct_r(2)
945 IF (rho_g_valid) THEN
946 ALLOCATE (rho_g(2))
947 rho_g(1) = rho_struct_g(2)
948 rho_g(2) = rho_struct_g(2)
949 END IF
950
951 IF (my_just_energy) THEN
952 exc_m = xc_exc_calc(rho_r=rho_r, tau=tau, &
953 rho_g=rho_g, xc_section=xc_section, &
954 weights=weights_use, pw_pool=xc_pw_pool)
955 ELSE
956 ! virial untested
957 cpassert(.NOT. compute_virial)
958 CALL xc_vxc_pw_create(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_r, &
959 rho_g=rho_g, tau=tau, exc=exc_m, &
960 xc_section=xc_section, &
961 weights=weights_use, pw_pool=xc_pw_pool, &
962 compute_virial=.false., &
963 virial_xc=virial_xc_tmp)
964 END IF
965
966 exc = exc + dft_control%sic_scaling_b*exc_m
967
968 ! and take care of the potential
969 IF (.NOT. my_just_energy) THEN
970 ! both go to minority spin
971 CALL pw_axpy(my_vxc_rho(1), vxc_rho(2), 2.0_dp*dft_control%sic_scaling_b)
972 CALL my_vxc_rho(1)%release()
973 CALL my_vxc_rho(2)%release()
974 DEALLOCATE (my_vxc_rho)
975 END IF
976 DEALLOCATE (rho_r, rho_g)
977
978 END IF
979
980 !
981 ! cleanups
982 !
983 IF (uf_grid .AND. (ASSOCIATED(vxc_rho) .OR. ASSOCIATED(vxc_tau))) THEN
984 block
985 TYPE(pw_r3d_rs_type) :: tmp_pw
986 TYPE(pw_c1d_gs_type) :: tmp_g, tmp_g2
987 CALL xc_pw_pool%create_pw(tmp_g)
988 CALL auxbas_pw_pool%create_pw(tmp_g2)
989 IF (ASSOCIATED(vxc_rho)) THEN
990 DO ispin = 1, SIZE(vxc_rho)
991 CALL auxbas_pw_pool%create_pw(tmp_pw)
992 CALL pw_transfer(vxc_rho(ispin), tmp_g)
993 CALL pw_transfer(tmp_g, tmp_g2)
994 CALL pw_transfer(tmp_g2, tmp_pw)
995 CALL xc_pw_pool%give_back_pw(vxc_rho(ispin))
996 vxc_rho(ispin) = tmp_pw
997 END DO
998 END IF
999 IF (ASSOCIATED(vxc_tau)) THEN
1000 DO ispin = 1, SIZE(vxc_tau)
1001 CALL auxbas_pw_pool%create_pw(tmp_pw)
1002 CALL pw_transfer(vxc_tau(ispin), tmp_g)
1003 CALL pw_transfer(tmp_g, tmp_g2)
1004 CALL pw_transfer(tmp_g2, tmp_pw)
1005 CALL xc_pw_pool%give_back_pw(vxc_tau(ispin))
1006 vxc_tau(ispin) = tmp_pw
1007 END DO
1008 END IF
1009 CALL auxbas_pw_pool%give_back_pw(tmp_g2)
1010 CALL xc_pw_pool%give_back_pw(tmp_g)
1011 END block
1012 END IF
1013 IF (ASSOCIATED(tau) .AND. uf_grid) THEN
1014 DO ispin = 1, SIZE(tau)
1015 CALL xc_pw_pool%give_back_pw(tau(ispin))
1016 END DO
1017 DEALLOCATE (tau)
1018 END IF
1019 IF (ASSOCIATED(weights_xc)) THEN
1020 CALL xc_pw_pool%give_back_pw(weights_xc)
1021 DEALLOCATE (weights_xc)
1022 END IF
1023 IF (ASSOCIATED(rho_nlcc_xc)) THEN
1024 CALL xc_pw_pool%give_back_pw(rho_nlcc_xc)
1025 DEALLOCATE (rho_nlcc_xc)
1026 END IF
1027 IF (ASSOCIATED(rho_nlcc_g_xc)) THEN
1028 CALL xc_pw_pool%give_back_pw(rho_nlcc_g_xc)
1029 DEALLOCATE (rho_nlcc_g_xc)
1030 END IF
1031
1032 IF (ASSOCIATED(rho_hard_r)) THEN
1033 DO ispin = 1, SIZE(rho_hard_r)
1034 CALL auxbas_pw_pool%give_back_pw(rho_hard_r(ispin))
1035 CALL auxbas_pw_pool%give_back_pw(rho_hard_g(ispin))
1036 CALL auxbas_pw_pool%give_back_pw(tau_hard_r(ispin))
1037 CALL auxbas_pw_pool%give_back_pw(tau_hard_g(ispin))
1038 END DO
1039 DEALLOCATE (rho_hard_r, rho_hard_g, tau_hard_r, tau_hard_g)
1040 END IF
1041
1042 END IF
1043
1044 CALL timestop(handle)
1045
1046 END SUBROUTINE qs_vxc_create
1047
1048! **************************************************************************************************
1049!> \brief Check the SKALA energy derivative along the reconstructed GAPW one-center correction.
1050!> \param rho_combined_r combined smooth plus hard-minus-soft density
1051!> \param rho_combined_g reciprocal representation of rho_combined_r
1052!> \param tau_combined combined kinetic-energy density
1053!> \param rho_smooth_r smooth density on the model grid
1054!> \param rho_smooth_g reciprocal representation of rho_smooth_r
1055!> \param tau_smooth smooth kinetic-energy density on the model grid
1056!> \param vxc_rho derivative of the SKALA energy with respect to rho
1057!> \param vxc_tau derivative of the SKALA energy with respect to tau
1058!> \param xc_section ...
1059!> \param weights ...
1060!> \param pw_pool ...
1061!> \param particle_set ...
1062!> \param cell ...
1063!> \param para_env ...
1064! **************************************************************************************************
1065 SUBROUTINE diagnose_gapw_composite_direction( &
1066 rho_combined_r, rho_combined_g, tau_combined, rho_smooth_r, rho_smooth_g, tau_smooth, &
1067 vxc_rho, vxc_tau, xc_section, weights, pw_pool, particle_set, cell, para_env)
1068
1069 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_combined_r
1070 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_combined_g
1071 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: tau_combined, rho_smooth_r
1072 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_smooth_g
1073 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: tau_smooth, vxc_rho, vxc_tau
1074 TYPE(section_vals_type), POINTER :: xc_section
1075 TYPE(pw_r3d_rs_type), POINTER :: weights
1076 TYPE(pw_pool_type), POINTER :: pw_pool
1077 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1078 TYPE(cell_type), POINTER :: cell
1079 TYPE(mp_para_env_type), POINTER :: para_env
1080
1081 REAL(kind=dp), DIMENSION(3), PARAMETER :: fd_steps = [1.0e-2_dp, 3.0e-3_dp, 1.0e-3_dp]
1082
1083 INTEGER :: ispin, istep, nspin, output_unit
1084 REAL(kind=dp) :: exc_minus, exc_plus, fd_derivative, fd_step, rho_contraction, &
1085 rho_g_contraction, tau_contraction, tau_g_contraction, total_contraction
1086 REAL(kind=dp), DIMENSION(3, 3) :: virial_dummy
1087 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_test_g, tau_test_g, vxc_rho_g, &
1088 vxc_tau_g
1089 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_test_r, tau_test, vxc_dummy_rho, &
1090 vxc_dummy_tau
1091
1092 NULLIFY (rho_test_g, rho_test_r, tau_test, tau_test_g, vxc_dummy_rho, vxc_dummy_tau, &
1093 vxc_rho_g, vxc_tau_g)
1094 nspin = SIZE(rho_combined_r)
1095 cpassert(SIZE(rho_combined_g) == nspin)
1096 cpassert(SIZE(tau_combined) == nspin)
1097 cpassert(SIZE(rho_smooth_r) == nspin)
1098 cpassert(SIZE(rho_smooth_g) == nspin)
1099 cpassert(SIZE(tau_smooth) == nspin)
1100 cpassert(SIZE(vxc_rho) == nspin)
1101 cpassert(SIZE(vxc_tau) == nspin)
1102
1103 ALLOCATE (rho_test_g(nspin), rho_test_r(nspin), tau_test(nspin), tau_test_g(nspin), &
1104 vxc_rho_g(nspin), vxc_tau_g(nspin))
1105 DO ispin = 1, nspin
1106 CALL pw_pool%create_pw(rho_test_g(ispin))
1107 CALL pw_pool%create_pw(rho_test_r(ispin))
1108 CALL pw_pool%create_pw(tau_test(ispin))
1109 CALL pw_pool%create_pw(tau_test_g(ispin))
1110 CALL pw_pool%create_pw(vxc_rho_g(ispin))
1111 CALL pw_pool%create_pw(vxc_tau_g(ispin))
1112 CALL pw_copy(rho_combined_g(ispin), rho_test_g(ispin))
1113 CALL pw_axpy(rho_smooth_g(ispin), rho_test_g(ispin), -1.0_dp)
1114 CALL pw_copy(rho_combined_r(ispin), rho_test_r(ispin))
1115 CALL pw_axpy(rho_smooth_r(ispin), rho_test_r(ispin), -1.0_dp)
1116 CALL pw_copy(tau_combined(ispin), tau_test(ispin))
1117 CALL pw_axpy(tau_smooth(ispin), tau_test(ispin), -1.0_dp)
1118 CALL pw_transfer(tau_test(ispin), tau_test_g(ispin))
1119 CALL pw_transfer(vxc_rho(ispin), vxc_rho_g(ispin))
1120 CALL pw_transfer(vxc_tau(ispin), vxc_tau_g(ispin))
1121 END DO
1122 rho_contraction = 0.0_dp
1123 rho_g_contraction = 0.0_dp
1124 tau_contraction = 0.0_dp
1125 tau_g_contraction = 0.0_dp
1126 DO ispin = 1, nspin
1127 rho_contraction = rho_contraction + pw_integral_ab(rho_test_r(ispin), vxc_rho(ispin))
1128 rho_g_contraction = rho_g_contraction + pw_integral_ab(rho_test_g(ispin), vxc_rho_g(ispin))
1129 tau_contraction = tau_contraction + pw_integral_ab(tau_test(ispin), vxc_tau(ispin))
1130 tau_g_contraction = tau_g_contraction + pw_integral_ab(tau_test_g(ispin), vxc_tau_g(ispin))
1131 END DO
1132
1133 total_contraction = rho_contraction + tau_contraction
1134 DO istep = 1, SIZE(fd_steps)
1135 fd_step = fd_steps(istep)
1136 DO ispin = 1, nspin
1137 CALL pw_copy(rho_combined_g(ispin), rho_test_g(ispin))
1138 CALL pw_axpy(rho_combined_g(ispin), rho_test_g(ispin), fd_step)
1139 CALL pw_axpy(rho_smooth_g(ispin), rho_test_g(ispin), -fd_step)
1140 CALL pw_copy(rho_combined_r(ispin), rho_test_r(ispin))
1141 CALL pw_axpy(rho_combined_r(ispin), rho_test_r(ispin), fd_step)
1142 CALL pw_axpy(rho_smooth_r(ispin), rho_test_r(ispin), -fd_step)
1143 CALL pw_copy(tau_combined(ispin), tau_test(ispin))
1144 CALL pw_axpy(tau_combined(ispin), tau_test(ispin), fd_step)
1145 CALL pw_axpy(tau_smooth(ispin), tau_test(ispin), -fd_step)
1146 END DO
1147 CALL skala_gpw_eval( &
1148 vxc_rho=vxc_dummy_rho, vxc_tau=vxc_dummy_tau, exc=exc_plus, rho_r=rho_test_r, &
1149 rho_g=rho_test_g, tau=tau_test, xc_section=xc_section, weights=weights, pw_pool=pw_pool, &
1150 particle_set=particle_set, cell=cell, compute_virial=.false., virial_xc=virial_dummy, &
1151 just_energy=.true.)
1152 cpassert(.NOT. ASSOCIATED(vxc_dummy_rho))
1153 cpassert(.NOT. ASSOCIATED(vxc_dummy_tau))
1154
1155 DO ispin = 1, nspin
1156 CALL pw_copy(rho_combined_g(ispin), rho_test_g(ispin))
1157 CALL pw_axpy(rho_combined_g(ispin), rho_test_g(ispin), -fd_step)
1158 CALL pw_axpy(rho_smooth_g(ispin), rho_test_g(ispin), fd_step)
1159 CALL pw_copy(rho_combined_r(ispin), rho_test_r(ispin))
1160 CALL pw_axpy(rho_combined_r(ispin), rho_test_r(ispin), -fd_step)
1161 CALL pw_axpy(rho_smooth_r(ispin), rho_test_r(ispin), fd_step)
1162 CALL pw_copy(tau_combined(ispin), tau_test(ispin))
1163 CALL pw_axpy(tau_combined(ispin), tau_test(ispin), -fd_step)
1164 CALL pw_axpy(tau_smooth(ispin), tau_test(ispin), fd_step)
1165 END DO
1166 CALL skala_gpw_eval( &
1167 vxc_rho=vxc_dummy_rho, vxc_tau=vxc_dummy_tau, exc=exc_minus, rho_r=rho_test_r, &
1168 rho_g=rho_test_g, tau=tau_test, xc_section=xc_section, weights=weights, pw_pool=pw_pool, &
1169 particle_set=particle_set, cell=cell, compute_virial=.false., virial_xc=virial_dummy, &
1170 just_energy=.true.)
1171 cpassert(.NOT. ASSOCIATED(vxc_dummy_rho))
1172 cpassert(.NOT. ASSOCIATED(vxc_dummy_tau))
1173
1174 fd_derivative = (exc_plus - exc_minus)/(2.0_dp*fd_step)
1175 IF (para_env%mepos == 0) THEN
1176 output_unit = cp_logger_get_default_io_unit()
1177 IF (output_unit > 0) THEN
1178 WRITE (output_unit, "(/,T2,A,1X,ES12.4,2(1X,ES20.10))") &
1179 "SKALA_GPW| Composite-direction FD step and energies", &
1180 fd_step, exc_plus, exc_minus
1181 WRITE (output_unit, "(T2,A,1X,ES20.10)") &
1182 "SKALA_GPW| Composite-direction energy FD", fd_derivative
1183 WRITE (output_unit, "(T2,A,1X,ES20.10)") &
1184 "SKALA_GPW| Composite-direction FD/VXC difference", &
1185 fd_derivative - total_contraction
1186 END IF
1187 END IF
1188 END DO
1189 IF (para_env%mepos == 0) THEN
1190 output_unit = cp_logger_get_default_io_unit()
1191 IF (output_unit > 0) THEN
1192 WRITE (output_unit, "(T2,A,1X,ES20.10)") &
1193 "SKALA_GPW| Composite-direction rho VXC contraction", rho_contraction
1194 WRITE (output_unit, "(T2,A,1X,ES20.10)") &
1195 "SKALA_GPW| Composite-direction rho G-space contraction", rho_g_contraction
1196 WRITE (output_unit, "(T2,A,1X,ES20.10)") &
1197 "SKALA_GPW| Composite-direction tau VXC contraction", tau_contraction
1198 WRITE (output_unit, "(T2,A,1X,ES20.10)") &
1199 "SKALA_GPW| Composite-direction tau G-space contraction", tau_g_contraction
1200 WRITE (output_unit, "(T2,A,1X,ES20.10)") &
1201 "SKALA_GPW| Composite-direction total VXC contraction", total_contraction
1202 END IF
1203 END IF
1204
1205 DO ispin = 1, nspin
1206 CALL pw_pool%give_back_pw(rho_test_g(ispin))
1207 CALL pw_pool%give_back_pw(rho_test_r(ispin))
1208 CALL pw_pool%give_back_pw(tau_test(ispin))
1209 CALL pw_pool%give_back_pw(tau_test_g(ispin))
1210 CALL pw_pool%give_back_pw(vxc_rho_g(ispin))
1211 CALL pw_pool%give_back_pw(vxc_tau_g(ispin))
1212 END DO
1213 DEALLOCATE (rho_test_g, rho_test_r, tau_test, tau_test_g, vxc_rho_g, vxc_tau_g)
1214
1215 END SUBROUTINE diagnose_gapw_composite_direction
1216
1217! **************************************************************************************************
1218!> \brief calculates the XC density: E_xc(r) - V_xc(r)*rho(r) or E_xc(r)/rho(r)
1219!> \param ks_env to get all the needed things
1220!> \param rho_struct density
1221!> \param xc_section ...
1222!> \param dispersion_env ...
1223!> \param xc_ener will contain the xc energy density E_xc(r) - V_xc(r)*rho(r)
1224!> \param xc_den will contain the xc energy density E_xc(r)/rho(r)
1225!> \param exc will contain the xc energy density E_xc(r)
1226!> \param vxc ...
1227!> \param vtau ...
1228!> \author JGH
1229! **************************************************************************************************
1230 SUBROUTINE qs_xc_density(ks_env, rho_struct, xc_section, dispersion_env, &
1231 xc_ener, xc_den, exc, vxc, vtau)
1232
1233 TYPE(qs_ks_env_type), POINTER :: ks_env
1234 TYPE(qs_rho_type), POINTER :: rho_struct
1235 TYPE(section_vals_type), POINTER :: xc_section
1236 TYPE(qs_dispersion_type), OPTIONAL, POINTER :: dispersion_env
1237 TYPE(pw_r3d_rs_type), INTENT(INOUT), OPTIONAL :: xc_ener, xc_den
1238 TYPE(pw_r3d_rs_type), OPTIONAL :: exc
1239 TYPE(pw_r3d_rs_type), DIMENSION(:), OPTIONAL :: vxc, vtau
1240
1241 CHARACTER(len=*), PARAMETER :: routinen = 'qs_xc_density'
1242
1243 INTEGER :: handle, ispin, mspin, myfun, np, nspins, &
1244 vdw
1245 LOGICAL :: rho_g_valid, tau_g_valid, tau_r_valid, &
1246 uf_grid, vdw_nl
1247 REAL(kind=dp) :: edisp, excint, factor, rho_cutoff
1248 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: rhonl
1249 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: drhonl
1250 REAL(kind=dp), DIMENSION(3, 3) :: vdum
1251 TYPE(cell_type), POINTER :: cell
1252 TYPE(dft_control_type), POINTER :: dft_control
1253 TYPE(mp_para_env_type), POINTER :: para_env
1254 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g, rho_struct_g, tau_g, tau_struct_g
1255 TYPE(pw_c1d_gs_type), POINTER :: rho_nlcc_g, rho_nlcc_g_use, rho_nlcc_g_xc
1256 TYPE(pw_env_type), POINTER :: pw_env
1257 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool, vdw_pw_pool, xc_pw_pool
1258 TYPE(pw_r3d_rs_type) :: exc_r, vxcnl
1259 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, rho_struct_r, tau_r, &
1260 tau_struct_r, vxc_rho, vxc_tau
1261 TYPE(pw_r3d_rs_type), POINTER :: rho_nlcc, rho_nlcc_use, rho_nlcc_xc, &
1262 weights, weights_use, weights_xc
1263
1264 CALL timeset(routinen, handle)
1265
1266 NULLIFY (dft_control, pw_env, auxbas_pw_pool, xc_pw_pool, vdw_pw_pool, cell, &
1267 rho_g, rho_struct_g, tau_g, tau_struct_g, rho_nlcc, rho_nlcc_g, &
1268 rho_nlcc_g_use, rho_nlcc_g_xc, rho_nlcc_use, rho_nlcc_xc, rho_r, &
1269 rho_struct_r, tau_r, tau_struct_r, vxc_rho, vxc_tau, weights, &
1270 weights_use, weights_xc)
1271
1272 CALL get_ks_env(ks_env, &
1273 dft_control=dft_control, &
1274 pw_env=pw_env, &
1275 cell=cell, &
1276 xcint_weights=weights, &
1277 rho_nlcc=rho_nlcc, &
1278 rho_nlcc_g=rho_nlcc_g)
1279
1280 CALL qs_rho_get(rho_struct, &
1281 tau_r_valid=tau_r_valid, &
1282 tau_g_valid=tau_g_valid, &
1283 rho_g_valid=rho_g_valid, &
1284 rho_r=rho_struct_r, &
1285 rho_g=rho_struct_g, &
1286 tau_r=tau_struct_r, &
1287 tau_g=tau_struct_g)
1288 nspins = dft_control%nspins
1289 mspin = SIZE(rho_struct_r)
1290 rho_r => rho_struct_r
1291 rho_g => rho_struct_g
1292 tau_r => tau_struct_r
1293 tau_g => tau_struct_g
1294 rho_nlcc_use => rho_nlcc
1295 rho_nlcc_g_use => rho_nlcc_g
1296 weights_use => weights
1297
1298 CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", i_val=myfun)
1299 CALL section_vals_val_get(xc_section, "VDW_POTENTIAL%POTENTIAL_TYPE", i_val=vdw)
1300 vdw_nl = (vdw == xc_vdw_fun_nonloc)
1301 IF (PRESENT(xc_ener)) THEN
1302 IF (tau_r_valid) THEN
1303 CALL cp_warn(__location__, "Tau contribution will not be correctly handled")
1304 END IF
1305 END IF
1306 IF (vdw_nl) THEN
1307 CALL cp_warn(__location__, "vdW functional contribution will be ignored")
1308 END IF
1309
1310 CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool, auxbas_pw_pool=auxbas_pw_pool)
1311 uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid, xc_pw_pool%pw_grid)
1312
1313 IF (PRESENT(xc_ener)) THEN
1314 CALL pw_zero(xc_ener)
1315 END IF
1316 IF (PRESENT(xc_den)) THEN
1317 CALL pw_zero(xc_den)
1318 END IF
1319 IF (PRESENT(exc)) THEN
1320 CALL pw_zero(exc)
1321 END IF
1322 IF (PRESENT(vxc)) THEN
1323 DO ispin = 1, nspins
1324 CALL pw_zero(vxc(ispin))
1325 END DO
1326 END IF
1327 IF (PRESENT(vtau)) THEN
1328 DO ispin = 1, nspins
1329 CALL pw_zero(vtau(ispin))
1330 END DO
1331 END IF
1332
1333 IF (myfun /= xc_none) THEN
1334
1335 cpassert(ASSOCIATED(rho_struct))
1336 cpassert(dft_control%sic_method_id == sic_none)
1337
1338 IF (uf_grid) THEN
1339 NULLIFY (rho_r, rho_g, tau_r, tau_g)
1340 IF (rho_g_valid) THEN
1341 CALL create_density_on_pool(xc_pw_pool, rho_struct_g, rho_r, rho_g)
1342 ELSE IF (ASSOCIATED(rho_struct_r)) THEN
1343 CALL create_density_on_pool_from_r(auxbas_pw_pool, xc_pw_pool, rho_struct_r, rho_r, rho_g)
1344 ELSE
1345 cpabort("Fine Grid in qs_xc_density requires rho_r or rho_g")
1346 END IF
1347 IF (tau_r_valid) THEN
1348 IF (tau_g_valid) THEN
1349 CALL create_density_on_pool(xc_pw_pool, tau_struct_g, tau_r, tau_g)
1350 ELSE IF (ASSOCIATED(tau_struct_r)) THEN
1351 CALL create_density_on_pool_from_r(auxbas_pw_pool, xc_pw_pool, tau_struct_r, tau_r, tau_g)
1352 ELSE
1353 cpabort("Fine Grid in qs_xc_density requires tau_r or tau_g")
1354 END IF
1355 END IF
1356 IF (ASSOCIATED(weights)) THEN
1357 ALLOCATE (weights_xc)
1358 CALL xc_pw_pool%create_pw(weights_xc)
1359 CALL transfer_rspace_between_pools(auxbas_pw_pool, xc_pw_pool, weights, weights_xc)
1360 weights_use => weights_xc
1361 END IF
1362 IF (ASSOCIATED(rho_nlcc)) THEN
1363 cpassert(ASSOCIATED(rho_nlcc_g))
1364 ALLOCATE (rho_nlcc_g_xc, rho_nlcc_xc)
1365 CALL xc_pw_pool%create_pw(rho_nlcc_g_xc)
1366 CALL xc_pw_pool%create_pw(rho_nlcc_xc)
1367 CALL pw_transfer(rho_nlcc_g, rho_nlcc_g_xc)
1368 CALL pw_transfer(rho_nlcc_g_xc, rho_nlcc_xc)
1369 rho_nlcc_use => rho_nlcc_xc
1370 rho_nlcc_g_use => rho_nlcc_g_xc
1371 END IF
1372 END IF
1373
1374 ! add the nlcc densities
1375 IF (ASSOCIATED(rho_nlcc_use)) THEN
1376 factor = 1.0_dp
1377 DO ispin = 1, mspin
1378 CALL pw_axpy(rho_nlcc_use, rho_r(ispin), factor)
1379 CALL pw_axpy(rho_nlcc_g_use, rho_g(ispin), factor)
1380 END DO
1381 END IF
1382 NULLIFY (vxc_rho, vxc_tau)
1383 CALL xc_vxc_pw_create(vxc_rho=vxc_rho, vxc_tau=vxc_tau, rho_r=rho_r, &
1384 rho_g=rho_g, tau=tau_r, exc=excint, &
1385 xc_section=xc_section, &
1386 weights=weights_use, pw_pool=xc_pw_pool, &
1387 compute_virial=.false., &
1388 virial_xc=vdum, &
1389 exc_r=exc_r)
1390 ! calclulate non-local vdW functional
1391 ! only if this XC_SECTION has it
1392 ! if yes, we use the dispersion_env from ks_env
1393 ! this is dangerous, as it assumes a special connection xc_section -> qs_env
1394 IF (vdw_nl) THEN
1395 CALL get_ks_env(ks_env=ks_env, para_env=para_env)
1396 ! no SIC functionals allowed
1397 cpassert(dft_control%sic_method_id == sic_none)
1398 !
1399 CALL pw_env_get(pw_env, vdw_pw_pool=vdw_pw_pool)
1400 np = vdw_pw_pool%pw_grid%ngpts_local
1401 ALLOCATE (rhonl(np), drhonl(np, 3))
1402 CALL prepare_nlvdw_density(rho_r, rhonl, drhonl, vdw_pw_pool, xc_pw_pool)
1403 CALL vdw_pw_pool%create_pw(vxcnl)
1404 CALL calculate_vxc_nlvdw(vxcnl, rhonl, drhonl, edisp, dispersion_env, &
1405 .false., vdw_pw_pool, para_env)
1406 CALL return_nlvdw_potential(vxc_rho, vxcnl, vdw_pw_pool, xc_pw_pool)
1407 CALL vdw_pw_pool%give_back_pw(vxcnl)
1408 DEALLOCATE (rhonl, drhonl)
1409 END IF
1410
1411 ! remove the nlcc densities (keep stuff in original state)
1412 IF (ASSOCIATED(rho_nlcc_use)) THEN
1413 factor = -1.0_dp
1414 DO ispin = 1, mspin
1415 CALL pw_axpy(rho_nlcc_use, rho_r(ispin), factor)
1416 CALL pw_axpy(rho_nlcc_g_use, rho_g(ispin), factor)
1417 END DO
1418 END IF
1419 !
1420 IF (PRESENT(xc_den)) THEN
1421 rho_cutoff = 1.e-14_dp
1422 IF (uf_grid) THEN
1423 block
1424 TYPE(pw_r3d_rs_type) :: tmp_pw
1425 CALL xc_pw_pool%create_pw(tmp_pw)
1426 CALL pw_copy(exc_r, tmp_pw)
1427 CALL calc_xc_density(tmp_pw, rho_r, rho_cutoff)
1428 CALL transfer_rspace_between_pools(xc_pw_pool, auxbas_pw_pool, tmp_pw, xc_den)
1429 CALL xc_pw_pool%give_back_pw(tmp_pw)
1430 END block
1431 ELSE
1432 CALL pw_copy(exc_r, xc_den)
1433 CALL calc_xc_density(xc_den, rho_r, rho_cutoff)
1434 END IF
1435 END IF
1436 IF (PRESENT(xc_ener)) THEN
1437 IF (uf_grid) THEN
1438 block
1439 TYPE(pw_r3d_rs_type) :: tmp_pw
1440 CALL xc_pw_pool%create_pw(tmp_pw)
1441 CALL pw_copy(exc_r, tmp_pw)
1442 DO ispin = 1, nspins
1443 CALL pw_multiply(tmp_pw, vxc_rho(ispin), rho_r(ispin), alpha=-1.0_dp)
1444 END DO
1445 CALL transfer_rspace_between_pools(xc_pw_pool, auxbas_pw_pool, tmp_pw, xc_ener)
1446 CALL xc_pw_pool%give_back_pw(tmp_pw)
1447 END block
1448 ELSE
1449 CALL pw_copy(exc_r, xc_ener)
1450 DO ispin = 1, nspins
1451 CALL pw_multiply(xc_ener, vxc_rho(ispin), rho_r(ispin), alpha=-1.0_dp)
1452 END DO
1453 END IF
1454 END IF
1455 IF (PRESENT(exc)) THEN
1456 IF (uf_grid) THEN
1457 CALL transfer_rspace_between_pools(xc_pw_pool, auxbas_pw_pool, exc_r, exc)
1458 ELSE
1459 CALL pw_copy(exc_r, exc)
1460 END IF
1461 END IF
1462 IF (PRESENT(vxc)) THEN
1463 DO ispin = 1, nspins
1464 IF (uf_grid) THEN
1465 CALL transfer_rspace_between_pools(xc_pw_pool, auxbas_pw_pool, vxc_rho(ispin), vxc(ispin))
1466 ELSE
1467 CALL pw_copy(vxc_rho(ispin), vxc(ispin))
1468 END IF
1469 END DO
1470 END IF
1471 IF (PRESENT(vtau) .AND. ASSOCIATED(vxc_tau)) THEN
1472 DO ispin = 1, nspins
1473 IF (uf_grid) THEN
1474 CALL transfer_rspace_between_pools(xc_pw_pool, auxbas_pw_pool, vxc_tau(ispin), vtau(ispin))
1475 ELSE
1476 CALL pw_copy(vxc_tau(ispin), vtau(ispin))
1477 END IF
1478 END DO
1479 END IF
1480 ! remove arrays
1481 IF (ASSOCIATED(vxc_rho)) THEN
1482 DO ispin = 1, nspins
1483 CALL vxc_rho(ispin)%release()
1484 END DO
1485 DEALLOCATE (vxc_rho)
1486 END IF
1487 IF (ASSOCIATED(vxc_tau)) THEN
1488 DO ispin = 1, nspins
1489 CALL vxc_tau(ispin)%release()
1490 END DO
1491 DEALLOCATE (vxc_tau)
1492 END IF
1493 CALL exc_r%release()
1494 IF (uf_grid) THEN
1495 CALL give_back_density_on_pool(xc_pw_pool, rho_r, rho_g)
1496 IF (ASSOCIATED(tau_r)) CALL give_back_density_on_pool(xc_pw_pool, tau_r, tau_g)
1497 IF (ASSOCIATED(weights_xc)) THEN
1498 CALL xc_pw_pool%give_back_pw(weights_xc)
1499 DEALLOCATE (weights_xc)
1500 END IF
1501 IF (ASSOCIATED(rho_nlcc_xc)) THEN
1502 CALL xc_pw_pool%give_back_pw(rho_nlcc_xc)
1503 DEALLOCATE (rho_nlcc_xc)
1504 END IF
1505 IF (ASSOCIATED(rho_nlcc_g_xc)) THEN
1506 CALL xc_pw_pool%give_back_pw(rho_nlcc_g_xc)
1507 DEALLOCATE (rho_nlcc_g_xc)
1508 END IF
1509 END IF
1510 !
1511 END IF
1512
1513 CALL timestop(handle)
1514
1515 END SUBROUTINE qs_xc_density
1516
1517! **************************************************************************************************
1518!> \brief transfers an r-space PW between two pools and writes into an existing target PW
1519!> \param source_pw_pool ...
1520!> \param target_pw_pool ...
1521!> \param source ...
1522!> \param TARGET ...
1523! **************************************************************************************************
1524 SUBROUTINE transfer_rspace_between_pools(source_pw_pool, target_pw_pool, source, TARGET)
1525 TYPE(pw_pool_type), POINTER :: source_pw_pool, target_pw_pool
1526 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: source, target
1527
1528 TYPE(pw_c1d_gs_type) :: source_g, target_g
1529
1530 cpassert(ASSOCIATED(source_pw_pool))
1531 cpassert(ASSOCIATED(target_pw_pool))
1532
1533 IF (pw_grid_compare(source_pw_pool%pw_grid, target_pw_pool%pw_grid)) THEN
1534 CALL pw_copy(source, TARGET)
1535 ELSE
1536 CALL source_pw_pool%create_pw(source_g)
1537 CALL target_pw_pool%create_pw(target_g)
1538 CALL pw_transfer(source, source_g)
1539 CALL pw_transfer(source_g, target_g)
1540 CALL pw_transfer(target_g, TARGET)
1541 CALL target_pw_pool%give_back_pw(target_g)
1542 CALL source_pw_pool%give_back_pw(source_g)
1543 END IF
1544
1545 END SUBROUTINE transfer_rspace_between_pools
1546
1547! **************************************************************************************************
1548!> \brief transfers a g-space density to a given PW pool and creates its r-space representation
1549!> \param pw_pool ...
1550!> \param rho_g_in ...
1551!> \param rho_r_out ...
1552!> \param rho_g_out ...
1553! **************************************************************************************************
1554 SUBROUTINE create_density_on_pool(pw_pool, rho_g_in, rho_r_out, rho_g_out)
1555 TYPE(pw_pool_type), POINTER :: pw_pool
1556 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_in
1557 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_out
1558 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_out
1559
1560 INTEGER :: ispin, nspins
1561
1562 cpassert(ASSOCIATED(pw_pool))
1563 cpassert(ASSOCIATED(rho_g_in))
1564
1565 nspins = SIZE(rho_g_in)
1566 ALLOCATE (rho_r_out(nspins), rho_g_out(nspins))
1567 DO ispin = 1, nspins
1568 CALL pw_pool%create_pw(rho_g_out(ispin))
1569 CALL pw_pool%create_pw(rho_r_out(ispin))
1570 CALL pw_transfer(rho_g_in(ispin), rho_g_out(ispin))
1571 CALL pw_transfer(rho_g_out(ispin), rho_r_out(ispin))
1572 END DO
1573
1574 END SUBROUTINE create_density_on_pool
1575
1576! **************************************************************************************************
1577!> \brief transfers an r-space density to a given PW pool and creates its g-space representation
1578!> \param source_pw_pool ...
1579!> \param target_pw_pool ...
1580!> \param rho_r_in ...
1581!> \param rho_r_out ...
1582!> \param rho_g_out ...
1583! **************************************************************************************************
1584 SUBROUTINE create_density_on_pool_from_r(source_pw_pool, target_pw_pool, rho_r_in, rho_r_out, rho_g_out)
1585 TYPE(pw_pool_type), POINTER :: source_pw_pool, target_pw_pool
1586 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_in, rho_r_out
1587 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_out
1588
1589 INTEGER :: ispin, nspins
1590 TYPE(pw_c1d_gs_type) :: rho_g_in
1591
1592 cpassert(ASSOCIATED(source_pw_pool))
1593 cpassert(ASSOCIATED(target_pw_pool))
1594 cpassert(ASSOCIATED(rho_r_in))
1595
1596 nspins = SIZE(rho_r_in)
1597 ALLOCATE (rho_r_out(nspins), rho_g_out(nspins))
1598 DO ispin = 1, nspins
1599 CALL source_pw_pool%create_pw(rho_g_in)
1600 CALL target_pw_pool%create_pw(rho_g_out(ispin))
1601 CALL target_pw_pool%create_pw(rho_r_out(ispin))
1602 CALL pw_transfer(rho_r_in(ispin), rho_g_in)
1603 CALL pw_transfer(rho_g_in, rho_g_out(ispin))
1604 CALL pw_transfer(rho_g_out(ispin), rho_r_out(ispin))
1605 CALL source_pw_pool%give_back_pw(rho_g_in)
1606 END DO
1607
1608 END SUBROUTINE create_density_on_pool_from_r
1609
1610! **************************************************************************************************
1611!> \brief returns temporary density arrays to the given PW pool
1612!> \param pw_pool ...
1613!> \param rho_r ...
1614!> \param rho_g ...
1615! **************************************************************************************************
1616 SUBROUTINE give_back_density_on_pool(pw_pool, rho_r, rho_g)
1617 TYPE(pw_pool_type), POINTER :: pw_pool
1618 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
1619 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
1620
1621 INTEGER :: ispin
1622
1623 cpassert(ASSOCIATED(pw_pool))
1624
1625 IF (ASSOCIATED(rho_r)) THEN
1626 DO ispin = 1, SIZE(rho_r)
1627 CALL pw_pool%give_back_pw(rho_r(ispin))
1628 END DO
1629 DEALLOCATE (rho_r)
1630 END IF
1631 IF (ASSOCIATED(rho_g)) THEN
1632 DO ispin = 1, SIZE(rho_g)
1633 CALL pw_pool%give_back_pw(rho_g(ispin))
1634 END DO
1635 DEALLOCATE (rho_g)
1636 END IF
1637
1638 END SUBROUTINE give_back_density_on_pool
1639
1640END MODULE qs_vxc
pw_methods::pw_axpy
Definition pw_methods.F:165
pw_methods::pw_copy
Definition pw_methods.F:146
pw_methods::pw_integral_ab
Definition pw_methods.F:222
pw_methods::pw_integrate_function
Definition pw_methods.F:116
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
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:523
gapw_gspace_reconstruction
GAPW reciprocal-space reconstruction and its discrete adjoint.
Definition gapw_gspace_reconstruction.F:12
gapw_gspace_reconstruction::calculate_rhotot_elec_gspace
subroutine, public calculate_rhotot_elec_gspace(qs_env, auxbas_pw_pool, rhotot_elec_gspace, q_max, rho_hard, rho_soft, fsign, compute_tau, rho_source, allow_nonorthorhombic)
The total electronic density in reciprocal space (g-space) is calculated.
Definition gapw_gspace_reconstruction.F:88
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::sic_mauri_spz
integer, parameter, public sic_mauri_spz
Definition input_constants.F:622
input_constants::xc_vdw_fun_nonloc
integer, parameter, public xc_vdw_fun_nonloc
Definition input_constants.F:649
input_constants::sic_eo
integer, parameter, public sic_eo
Definition input_constants.F:622
input_constants::sic_mauri_us
integer, parameter, public sic_mauri_us
Definition input_constants.F:622
input_constants::sic_none
integer, parameter, public sic_none
Definition input_constants.F:622
input_constants::xc_none
integer, parameter, public xc_none
Definition input_constants.F:605
input_constants::sic_ad
integer, parameter, public sic_ad
Definition input_constants.F:622
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:1052
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:27
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
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_grids::get_pw_grid_info
subroutine, public get_pw_grid_info(pw_grid, id_nr, mode, vol, dvol, npts, ngpts, ngpts_cut, dr, cutoff, orthorhombic, gvectors, gsquare)
Access to information stored in the pw_grid_type.
Definition pw_grids.F:185
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_collocate_density
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
Definition qs_collocate_density.F:38
qs_collocate_density::calculate_rho_elec
subroutine, public calculate_rho_elec(matrix_p, matrix_p_kp, rho, rho_gspace, total_rho, ks_env, soft_valid, compute_tau, compute_grad, basis_type, der_type, idir, task_list_external, pw_env_external)
computes the density corresponding to a given density matrix on the grid
Definition qs_collocate_density.F:1712
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::prepare_nlvdw_density
subroutine, public prepare_nlvdw_density(rho_r, rho, drho, pw_pool, xc_pw_pool)
Transform the input density to a 1d grid used by the nl-vdW routines.
Definition qs_dispersion_nonloc.F:535
qs_dispersion_nonloc::calculate_vxc_nlvdw
subroutine, public calculate_vxc_nlvdw(vxc_r, rho, drho, edispersion, dispersion_env, energy_only, 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:164
qs_dispersion_nonloc::return_nlvdw_potential
subroutine, public return_nlvdw_potential(vxc_rho, vxc_r, pw_pool, xc_pw_pool)
Take the nl-vdW potential and return it on the XC grid and for all spin channels.
Definition qs_dispersion_nonloc.F:683
qs_dispersion_types
Definition of disperson types for DFT calculations.
Definition qs_dispersion_types.F:12
qs_environment_types
Definition qs_environment_types.F:14
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_vhxc, 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:343
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, native_skala_atom_force, qs_env_external, native_gapw_composite_override, native_skala_defer_to_atom_composite)
calculates and allocates the xc potential, already reducing it to the dependence on rho and the one o...
Definition qs_vxc.F:120
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:1232
skala_gpw_functional
Experimental CP2K-native GPW real-space-grid path for SKALA TorchScript models.
Definition skala_gpw_functional.F:11
skala_gpw_functional::skala_gpw_eval
subroutine, public skala_gpw_eval(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau, xc_section, weights, pw_pool, particle_set, cell, compute_virial, virial_xc, just_energy, atom_force)
Evaluate SKALA energy and first derivatives on a CP2K GPW grid.
Definition skala_gpw_functional.F:429
skala_gpw_functional::native_skala_gapw_composite_direct_ao
logical function, public native_skala_gapw_composite_direct_ao(xc_section)
Return true if the GAPW composite reference uses direct full-ORB collocation.
Definition skala_gpw_functional.F:260
skala_gpw_functional::xc_section_uses_native_skala_evaluator
logical function, public xc_section_uses_native_skala_evaluator(xc_section)
Return true when SKALA must be evaluated by the CP2K-native grid machinery.
Definition skala_gpw_functional.F:166
skala_gpw_functional::get_gauxc_section
type(section_vals_type) function, pointer, public get_gauxc_section(xc_section)
Return the first GAUXC functional subsection, if present.
Definition skala_gpw_functional.F:3061
skala_gpw_functional::native_skala_gapw_composite_reference
logical function, public native_skala_gapw_composite_reference(xc_section)
Return true if native SKALA should use the full GAPW ORB density on one common grid.
Definition skala_gpw_functional.F:181
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:795
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:478
xc::calc_xc_density
subroutine, public calc_xc_density(pot, rho, rho_cutoff)
Definition xc.F:406
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:786
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
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:774
particle_types::particle_type
Definition particle_types.F:35
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_environment_types::qs_environment_type
Definition qs_environment_types.F:223
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