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