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_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 :
80 : CONTAINS
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 169425 : 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 169425 : 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 169425 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: rhonl
152 169425 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: drhonl
153 : REAL(KIND=dp), DIMENSION(3, 3) :: virial_xc_tmp
154 : TYPE(cell_type), POINTER :: cell
155 169425 : 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 169425 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
159 : TYPE(pw_c1d_gs_type) :: tmp_g
160 169425 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_direct_g, rho_g, rho_hard_g, &
161 169425 : rho_m_gspace, rho_smooth_g, rho_smooth_model_g, rho_struct_g, tau_direct_g, tau_hard_g, &
162 169425 : 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 169425 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: my_vxc_rho, my_vxc_tau, rho_direct_r, &
169 169425 : rho_hard_r, rho_m_rspace, rho_r, rho_smooth_model_r, rho_smooth_r, rho_struct_r, tau, &
170 169425 : 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 169425 : CALL timeset(routineN, handle)
178 :
179 169425 : CPASSERT(.NOT. ASSOCIATED(vxc_rho))
180 169425 : CPASSERT(.NOT. ASSOCIATED(vxc_tau))
181 169425 : NULLIFY (dft_control, pw_env, auxbas_pw_pool, diagnostics_pw_pool, xc_pw_pool, vdw_pw_pool, &
182 169425 : cell, my_vxc_rho, &
183 169425 : tmp_pw, my_vxc_tau, rho_g, rho_r, tau, rho_m_rspace, &
184 169425 : rho_m_gspace, rho_nlcc, rho_nlcc_g, rho_nlcc_g_use, rho_nlcc_g_xc, &
185 169425 : rho_nlcc_use, rho_nlcc_xc, rho_struct_r, rho_struct_g, tau_struct_g, tau_struct_r, &
186 169425 : weights_use, weights_xc, particle_set, rho_ao_kp, rho_direct_g, rho_direct_r, &
187 169425 : rho_hard_g, rho_hard_r, rho_smooth_g, rho_smooth_model_g, rho_smooth_model_r, &
188 169425 : rho_smooth_r, tau_direct_g, tau_direct_r, tau_hard_g, tau_hard_r, tau_smooth_g, &
189 169425 : tau_smooth_model_g, tau_smooth_model_r, tau_smooth_r, gauxc_section)
190 :
191 169425 : exc = 0.0_dp
192 169425 : my_just_energy = .FALSE.
193 169425 : IF (PRESENT(just_energy)) my_just_energy = just_energy
194 169425 : my_adiabatic_rescale_factor = 1.0_dp
195 169425 : do_adiabatic_rescaling = .FALSE.
196 169425 : IF (PRESENT(adiabatic_rescale_factor)) THEN
197 44 : my_adiabatic_rescale_factor = adiabatic_rescale_factor
198 44 : 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 169425 : rho_nlcc_g=rho_nlcc_g)
211 169425 : rho_nlcc_use => rho_nlcc
212 169425 : rho_nlcc_g_use => rho_nlcc_g
213 169425 : 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 169425 : tau_r=tau_struct_r)
223 169425 : rho_smooth_r => rho_struct_r
224 169425 : rho_smooth_g => rho_struct_g
225 169425 : tau_smooth_r => tau_struct_r
226 169425 : tau_smooth_g => tau_struct_g
227 :
228 169425 : compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
229 169425 : IF (compute_virial) THEN
230 38506 : virial%pv_xc = 0.0_dp
231 : END IF
232 :
233 : CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", &
234 169425 : i_val=myfun)
235 : CALL section_vals_val_get(xc_section, "VDW_POTENTIAL%POTENTIAL_TYPE", &
236 169425 : i_val=vdw)
237 :
238 169425 : vdW_nl = (vdw == xc_vdw_fun_nonloc)
239 : ! this combination has not been investigated
240 169425 : CPASSERT(.NOT. (do_adiabatic_rescaling .AND. vdW_nl))
241 : ! are the necessary inputs available
242 169425 : IF (.NOT. (PRESENT(dispersion_env) .AND. PRESENT(edisp))) THEN
243 : vdW_nl = .FALSE.
244 : END IF
245 169425 : IF (PRESENT(edisp)) edisp = 0.0_dp
246 169425 : native_skala_grid = xc_section_uses_native_skala_evaluator(xc_section)
247 169425 : defer_native_skala_to_atom_composite = .FALSE.
248 169425 : IF (PRESENT(native_skala_defer_to_atom_composite)) THEN
249 292 : defer_native_skala_to_atom_composite = native_skala_defer_to_atom_composite
250 : END IF
251 169425 : 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 169425 : dft_control%qs_control%gapw_xc)
256 169425 : 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 300 : dft_control%qs_control%gapw_xc)
261 : END IF
262 : native_gapw_composite_direct_ao = native_gapw_composite_reference .AND. &
263 169129 : native_skala_gapw_composite_direct_ao(xc_section)
264 169421 : native_grid_diagnostics = .FALSE.
265 4 : IF (native_gapw_composite_reference) THEN
266 4 : gauxc_section => get_gauxc_section(xc_section)
267 4 : CPASSERT(ASSOCIATED(gauxc_section))
268 : CALL section_vals_val_get(gauxc_section, "NATIVE_GRID_DIAGNOSTICS", &
269 4 : l_val=native_grid_diagnostics)
270 : END IF
271 4 : IF (native_gapw_composite_reference) NULLIFY (weights_use)
272 :
273 169425 : IF (myfun /= xc_none .OR. vdW_nl) THEN
274 :
275 : ! test if the real space density is available
276 153953 : CPASSERT(ASSOCIATED(rho_struct))
277 153953 : IF (dft_control%nspins /= 1 .AND. dft_control%nspins /= 2) THEN
278 0 : CPABORT("nspins must be 1 or 2")
279 : END IF
280 153953 : mspin = SIZE(rho_struct_r)
281 153953 : IF (dft_control%nspins == 2 .AND. mspin == 1) THEN
282 0 : 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 153953 : my_scaling = 1.0_dp
292 154047 : 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 94 : CPASSERT(.NOT. tau_r_valid)
298 : CASE (sic_mauri_us)
299 0 : my_scaling = 1.0_dp - dft_control%sic_scaling_b
300 : ! no idea yet what to do here in that case
301 0 : 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 153953 : CPABORT("NYI")
307 : END SELECT
308 :
309 153953 : IF (dft_control%sic_scaling_b == 0.0_dp) THEN
310 : sic_scaling_b_zero = .TRUE.
311 : ELSE
312 153859 : sic_scaling_b_zero = .FALSE.
313 : END IF
314 :
315 153953 : IF (PRESENT(pw_env_external)) THEN
316 0 : pw_env => pw_env_external
317 : END IF
318 153953 : CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool, auxbas_pw_pool=auxbas_pw_pool)
319 :
320 153953 : IF (native_gapw_composite_reference) THEN
321 4 : CPASSERT(tau_r_valid)
322 4 : CPASSERT(PRESENT(qs_env_external))
323 44 : ALLOCATE (rho_hard_r(mspin), rho_hard_g(mspin), tau_hard_r(mspin), tau_hard_g(mspin))
324 8 : DO ispin = 1, mspin
325 4 : CALL auxbas_pw_pool%create_pw(rho_hard_r(ispin))
326 4 : CALL auxbas_pw_pool%create_pw(rho_hard_g(ispin))
327 4 : CALL auxbas_pw_pool%create_pw(tau_hard_r(ispin))
328 8 : CALL auxbas_pw_pool%create_pw(tau_hard_g(ispin))
329 : END DO
330 4 : IF (native_gapw_composite_direct_ao) THEN
331 0 : CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
332 0 : CPASSERT(ASSOCIATED(rho_ao_kp))
333 0 : CPASSERT(SIZE(rho_ao_kp, 2) == 1)
334 0 : DO ispin = 1, mspin
335 : CALL calculate_rho_elec( &
336 : matrix_p=rho_ao_kp(ispin, 1)%matrix, rho=rho_hard_r(ispin), &
337 0 : 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 0 : rho_gspace=tau_hard_g(ispin), ks_env=ks_env, compute_tau=.TRUE.)
341 : END DO
342 : ELSE
343 4 : CALL get_pw_grid_info(rho_hard_g(1)%pw_grid, vol=volume)
344 27652 : 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 4 : rho_source=rho_struct, allow_nonorthorhombic=.TRUE.)
350 4 : rho_composite_hard = composite_hard_integral
351 4 : 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 4 : compute_tau=.TRUE., rho_source=rho_struct, allow_nonorthorhombic=.TRUE.)
357 4 : tau_composite_hard = composite_hard_integral
358 4 : tau_composite_soft = composite_soft_integral
359 4 : IF (para_env%mepos == 0) THEN
360 2 : output_unit = cp_logger_get_default_io_unit()
361 2 : IF (output_unit > 0) THEN
362 : WRITE (UNIT=output_unit, FMT="(/,T2,A,2(1X,ES19.11))") &
363 2 : "SKALA_GPW| Composite rho hard and soft integrals", &
364 4 : rho_composite_hard, rho_composite_soft
365 : WRITE (UNIT=output_unit, FMT="(T2,A,2(1X,ES19.11))") &
366 2 : "SKALA_GPW| Composite tau hard and soft integrals", &
367 4 : tau_composite_hard, tau_composite_soft
368 : END IF
369 : END IF
370 4 : IF (mspin == 1) THEN
371 4 : CALL pw_scale(rho_hard_g(1), 1.0_dp/volume)
372 4 : 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 0 : 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 0 : allow_nonorthorhombic=.TRUE.)
385 0 : CALL pw_scale(rho_hard_g(1), 0.5_dp/volume)
386 0 : CALL pw_scale(rho_hard_g(2), 0.5_dp/volume)
387 0 : CALL auxbas_pw_pool%create_pw(tmp_g)
388 0 : CALL pw_copy(rho_hard_g(1), tmp_g)
389 0 : CALL pw_axpy(rho_hard_g(2), rho_hard_g(1), 1.0_dp)
390 0 : CALL pw_axpy(rho_hard_g(2), tmp_g, -1.0_dp)
391 0 : CALL pw_copy(tmp_g, rho_hard_g(2))
392 0 : CALL pw_scale(tau_hard_g(1), 0.5_dp/volume)
393 0 : CALL pw_scale(tau_hard_g(2), 0.5_dp/volume)
394 0 : CALL pw_copy(tau_hard_g(1), tmp_g)
395 0 : CALL pw_axpy(tau_hard_g(2), tau_hard_g(1), 1.0_dp)
396 0 : CALL pw_axpy(tau_hard_g(2), tmp_g, -1.0_dp)
397 0 : CALL pw_copy(tmp_g, tau_hard_g(2))
398 0 : CALL auxbas_pw_pool%give_back_pw(tmp_g)
399 : END IF
400 : END IF
401 8 : DO ispin = 1, mspin
402 4 : CALL pw_transfer(rho_hard_g(ispin), rho_hard_r(ispin))
403 8 : CALL pw_transfer(tau_hard_g(ispin), tau_hard_r(ispin))
404 : END DO
405 4 : IF (native_grid_diagnostics .AND. .NOT. native_gapw_composite_direct_ao) THEN
406 0 : CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
407 0 : CPASSERT(ASSOCIATED(rho_ao_kp))
408 0 : CPASSERT(SIZE(rho_ao_kp, 2) == 1)
409 : ALLOCATE (rho_direct_r(mspin), rho_direct_g(mspin), &
410 0 : tau_direct_r(mspin), tau_direct_g(mspin))
411 0 : rho_diff_l2 = 0.0_dp
412 0 : rho_diff_max = 0.0_dp
413 0 : rho_ref_l2 = 0.0_dp
414 0 : tau_diff_l2 = 0.0_dp
415 0 : tau_diff_max = 0.0_dp
416 0 : tau_ref_l2 = 0.0_dp
417 0 : direct_rho_integral = 0.0_dp
418 0 : direct_tau_integral = 0.0_dp
419 0 : DO ispin = 1, mspin
420 0 : CALL auxbas_pw_pool%create_pw(rho_direct_r(ispin))
421 0 : CALL auxbas_pw_pool%create_pw(rho_direct_g(ispin))
422 0 : CALL auxbas_pw_pool%create_pw(tau_direct_r(ispin))
423 0 : 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 0 : 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 0 : rho_gspace=tau_direct_g(ispin), ks_env=ks_env, compute_tau=.TRUE.)
430 : direct_rho_integral = direct_rho_integral + &
431 0 : pw_integrate_function(rho_direct_r(ispin))
432 : direct_tau_integral = direct_tau_integral + &
433 0 : pw_integrate_function(tau_direct_r(ispin))
434 0 : DO k = LBOUND(rho_direct_r(ispin)%array, 3), UBOUND(rho_direct_r(ispin)%array, 3)
435 0 : DO j = LBOUND(rho_direct_r(ispin)%array, 2), UBOUND(rho_direct_r(ispin)%array, 2)
436 0 : 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 0 : rho_direct_r(ispin)%array(i, j, k)
439 0 : rho_diff_l2 = rho_diff_l2 + delta*delta
440 0 : rho_diff_max = MAX(rho_diff_max, ABS(delta))
441 0 : 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 0 : tau_direct_r(ispin)%array(i, j, k)
444 0 : tau_diff_l2 = tau_diff_l2 + delta*delta
445 0 : tau_diff_max = MAX(tau_diff_max, ABS(delta))
446 0 : 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 0 : CALL para_env%sum(rho_diff_l2)
452 0 : CALL para_env%sum(rho_ref_l2)
453 0 : CALL para_env%sum(tau_diff_l2)
454 0 : CALL para_env%sum(tau_ref_l2)
455 0 : CALL para_env%max(rho_diff_max)
456 0 : CALL para_env%max(tau_diff_max)
457 0 : rho_diff_l2 = SQRT(rho_diff_l2/MAX(rho_ref_l2, TINY(1.0_dp)))
458 0 : tau_diff_l2 = SQRT(tau_diff_l2/MAX(tau_ref_l2, TINY(1.0_dp)))
459 0 : IF (para_env%mepos == 0) THEN
460 0 : output_unit = cp_logger_get_default_io_unit()
461 0 : IF (output_unit > 0) THEN
462 : WRITE (UNIT=output_unit, FMT="(T2,A,2(1X,ES19.11))") &
463 0 : "SKALA_GPW| Direct AO rho and tau integrals", &
464 0 : direct_rho_integral, direct_tau_integral
465 : WRITE (UNIT=output_unit, FMT="(T2,A,2(1X,ES19.11))") &
466 0 : "SKALA_GPW| Composite/direct relative L2 rho and tau", &
467 0 : rho_diff_l2, tau_diff_l2
468 : WRITE (UNIT=output_unit, FMT="(T2,A,2(1X,ES19.11))") &
469 0 : "SKALA_GPW| Composite/direct Linf rho and tau", &
470 0 : rho_diff_max, tau_diff_max
471 : END IF
472 : END IF
473 0 : DO ispin = 1, mspin
474 0 : CALL auxbas_pw_pool%give_back_pw(rho_direct_r(ispin))
475 0 : CALL auxbas_pw_pool%give_back_pw(rho_direct_g(ispin))
476 0 : CALL auxbas_pw_pool%give_back_pw(tau_direct_r(ispin))
477 0 : CALL auxbas_pw_pool%give_back_pw(tau_direct_g(ispin))
478 : END DO
479 0 : DEALLOCATE (rho_direct_r, rho_direct_g, tau_direct_r, tau_direct_g)
480 : END IF
481 4 : composite_rho_r_integral = 0.0_dp
482 4 : composite_tau_r_integral = 0.0_dp
483 4 : composite_rho_min = HUGE(1.0_dp)
484 4 : composite_rho_max = -HUGE(1.0_dp)
485 4 : composite_tau_min = HUGE(1.0_dp)
486 4 : composite_tau_max = -HUGE(1.0_dp)
487 8 : DO ispin = 1, mspin
488 : composite_rho_r_integral = composite_rho_r_integral + &
489 4 : pw_integrate_function(rho_hard_r(ispin))
490 : composite_tau_r_integral = composite_tau_r_integral + &
491 4 : pw_integrate_function(tau_hard_r(ispin))
492 30052 : composite_rho_min = MIN(composite_rho_min, MINVAL(rho_hard_r(ispin)%array))
493 30052 : composite_rho_max = MAX(composite_rho_max, MAXVAL(rho_hard_r(ispin)%array))
494 30052 : composite_tau_min = MIN(composite_tau_min, MINVAL(tau_hard_r(ispin)%array))
495 30056 : composite_tau_max = MAX(composite_tau_max, MAXVAL(tau_hard_r(ispin)%array))
496 : END DO
497 4 : CALL para_env%min(composite_rho_min)
498 4 : CALL para_env%max(composite_rho_max)
499 4 : CALL para_env%min(composite_tau_min)
500 4 : CALL para_env%max(composite_tau_max)
501 4 : IF (para_env%mepos == 0) THEN
502 2 : output_unit = cp_logger_get_default_io_unit()
503 2 : IF (output_unit > 0) THEN
504 : WRITE (UNIT=output_unit, FMT="(T2,A,1X,ES19.11)") &
505 2 : "SKALA_GPW| Composite real-grid rho integral", composite_rho_r_integral
506 : WRITE (UNIT=output_unit, FMT="(T2,A,1X,ES19.11)") &
507 2 : "SKALA_GPW| Composite real-grid tau integral", composite_tau_r_integral
508 : WRITE (UNIT=output_unit, FMT="(T2,A,2(1X,ES19.11))") &
509 2 : "SKALA_GPW| Composite real-grid rho min and max", &
510 4 : composite_rho_min, composite_rho_max
511 : WRITE (UNIT=output_unit, FMT="(T2,A,2(1X,ES19.11))") &
512 2 : "SKALA_GPW| Composite real-grid tau min and max", &
513 4 : composite_tau_min, composite_tau_max
514 : END IF
515 : END IF
516 4 : rho_struct_r => rho_hard_r
517 4 : rho_struct_g => rho_hard_g
518 4 : tau_struct_r => tau_hard_r
519 4 : tau_struct_g => tau_hard_g
520 4 : rho_g_valid = .TRUE.
521 4 : tau_g_valid = .TRUE.
522 : END IF
523 :
524 153953 : uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid, xc_pw_pool%pw_grid)
525 :
526 153953 : IF (.NOT. uf_grid) THEN
527 152949 : rho_r => rho_struct_r
528 :
529 152949 : IF (tau_r_valid) THEN
530 4598 : 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 152949 : IF (rho_g_valid) THEN
536 152889 : rho_g => rho_struct_g
537 : END IF
538 : ELSE
539 1004 : CPASSERT(rho_g_valid)
540 4016 : ALLOCATE (rho_r(mspin))
541 4016 : ALLOCATE (rho_g(mspin))
542 2008 : DO ispin = 1, mspin
543 1004 : CALL xc_pw_pool%create_pw(rho_g(ispin))
544 2008 : CALL pw_transfer(rho_struct_g(ispin), rho_g(ispin))
545 : END DO
546 2008 : DO ispin = 1, mspin
547 1004 : CALL xc_pw_pool%create_pw(rho_r(ispin))
548 2008 : CALL pw_transfer(rho_g(ispin), rho_r(ispin))
549 : END DO
550 1004 : IF (tau_r_valid) THEN
551 666 : ALLOCATE (tau(mspin))
552 444 : DO ispin = 1, mspin
553 222 : CALL xc_pw_pool%create_pw(tau(ispin))
554 222 : BLOCK
555 : TYPE(pw_c1d_gs_type) :: tau_g_aux, tau_g_xc
556 222 : CALL xc_pw_pool%create_pw(tau_g_xc)
557 222 : IF (tau_g_valid) THEN
558 222 : CALL pw_transfer(tau_struct_g(ispin), tau_g_xc)
559 : ELSE
560 0 : CALL auxbas_pw_pool%create_pw(tau_g_aux)
561 0 : CALL pw_transfer(tau_struct_r(ispin), tau_g_aux)
562 0 : CALL pw_transfer(tau_g_aux, tau_g_xc)
563 0 : CALL auxbas_pw_pool%give_back_pw(tau_g_aux)
564 : END IF
565 222 : CALL pw_transfer(tau_g_xc, tau(ispin))
566 444 : CALL xc_pw_pool%give_back_pw(tau_g_xc)
567 : END BLOCK
568 : END DO
569 : END IF
570 1004 : IF (ASSOCIATED(weights) .AND. .NOT. native_gapw_composite_reference) THEN
571 990 : ALLOCATE (weights_xc)
572 990 : CALL xc_pw_pool%create_pw(weights_xc)
573 : BLOCK
574 : TYPE(pw_c1d_gs_type) :: weights_g_aux, weights_g_xc
575 990 : CALL auxbas_pw_pool%create_pw(weights_g_aux)
576 990 : CALL xc_pw_pool%create_pw(weights_g_xc)
577 990 : CALL pw_transfer(weights, weights_g_aux)
578 990 : CALL pw_transfer(weights_g_aux, weights_g_xc)
579 990 : CALL pw_transfer(weights_g_xc, weights_xc)
580 990 : CALL xc_pw_pool%give_back_pw(weights_g_xc)
581 1980 : CALL auxbas_pw_pool%give_back_pw(weights_g_aux)
582 : END BLOCK
583 990 : weights_use => weights_xc
584 : END IF
585 1004 : IF (ASSOCIATED(rho_nlcc)) THEN
586 56 : CPASSERT(ASSOCIATED(rho_nlcc_g))
587 56 : ALLOCATE (rho_nlcc_g_xc, rho_nlcc_xc)
588 56 : CALL xc_pw_pool%create_pw(rho_nlcc_g_xc)
589 56 : CALL xc_pw_pool%create_pw(rho_nlcc_xc)
590 56 : CALL pw_transfer(rho_nlcc_g, rho_nlcc_g_xc)
591 56 : CALL pw_transfer(rho_nlcc_g_xc, rho_nlcc_xc)
592 56 : rho_nlcc_use => rho_nlcc_xc
593 56 : rho_nlcc_g_use => rho_nlcc_g_xc
594 : END IF
595 : END IF
596 :
597 153953 : IF (native_gapw_composite_reference) THEN
598 4 : composite_rho_r_integral = 0.0_dp
599 4 : composite_tau_r_integral = 0.0_dp
600 8 : DO ispin = 1, mspin
601 : composite_rho_r_integral = composite_rho_r_integral + &
602 4 : pw_integrate_function(rho_r(ispin))
603 : composite_tau_r_integral = composite_tau_r_integral + &
604 8 : pw_integrate_function(tau(ispin))
605 : END DO
606 4 : IF (para_env%mepos == 0) THEN
607 2 : output_unit = cp_logger_get_default_io_unit()
608 2 : IF (output_unit > 0) THEN
609 : WRITE (UNIT=output_unit, FMT="(T2,A,1X,ES19.11)") &
610 2 : "SKALA_GPW| Composite XC-grid rho integral", composite_rho_r_integral
611 : WRITE (UNIT=output_unit, FMT="(T2,A,1X,ES19.11)") &
612 2 : "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 153953 : IF (ASSOCIATED(rho_nlcc_use)) THEN
619 2474 : factor = 1.0_dp
620 4948 : DO ispin = 1, mspin
621 2474 : CALL pw_axpy(rho_nlcc_use, rho_r(ispin), factor)
622 4948 : 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 153953 : 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 1596 : ALLOCATE (my_vxc_rho(mspin), my_vxc_tau(mspin))
635 546 : DO ispin = 1, mspin
636 294 : CALL xc_pw_pool%create_pw(my_vxc_rho(ispin))
637 294 : CALL xc_pw_pool%create_pw(my_vxc_tau(ispin))
638 294 : CALL pw_zero(my_vxc_rho(ispin))
639 546 : CALL pw_zero(my_vxc_tau(ispin))
640 : END DO
641 153701 : 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 88 : just_energy=my_just_energy, atom_force=native_skala_atom_force)
648 153653 : 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 11626 : 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 142027 : 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 153953 : IF (native_gapw_composite_reference .AND. .NOT. my_just_energy) THEN
665 4 : IF (uf_grid) THEN
666 0 : diagnostics_pw_pool => xc_pw_pool
667 : ELSE
668 4 : diagnostics_pw_pool => auxbas_pw_pool
669 : END IF
670 4 : CALL diagnostics_pw_pool%create_pw(tmp_g)
671 8 : DO ispin = 1, mspin
672 4 : CALL pw_transfer(my_vxc_rho(ispin), tmp_g)
673 4 : CALL pw_transfer(tmp_g, my_vxc_rho(ispin))
674 4 : CALL pw_transfer(my_vxc_tau(ispin), tmp_g)
675 8 : CALL pw_transfer(tmp_g, my_vxc_tau(ispin))
676 : END DO
677 4 : CALL diagnostics_pw_pool%give_back_pw(tmp_g)
678 4 : NULLIFY (diagnostics_pw_pool)
679 : END IF
680 :
681 4 : IF (native_gapw_composite_reference .AND. native_grid_diagnostics .AND. &
682 : .NOT. my_just_energy) THEN
683 0 : IF (uf_grid) THEN
684 0 : diagnostics_pw_pool => xc_pw_pool
685 : ELSE
686 0 : diagnostics_pw_pool => auxbas_pw_pool
687 : END IF
688 : ALLOCATE (rho_smooth_model_g(mspin), rho_smooth_model_r(mspin), &
689 0 : tau_smooth_model_g(mspin), tau_smooth_model_r(mspin))
690 0 : DO ispin = 1, mspin
691 0 : CALL diagnostics_pw_pool%create_pw(rho_smooth_model_g(ispin))
692 0 : CALL diagnostics_pw_pool%create_pw(rho_smooth_model_r(ispin))
693 0 : CALL diagnostics_pw_pool%create_pw(tau_smooth_model_g(ispin))
694 0 : CALL diagnostics_pw_pool%create_pw(tau_smooth_model_r(ispin))
695 0 : IF (uf_grid) THEN
696 0 : IF (rho_g_valid) THEN
697 0 : CALL pw_transfer(rho_smooth_g(ispin), rho_smooth_model_g(ispin))
698 : ELSE
699 0 : CALL auxbas_pw_pool%create_pw(tmp_g)
700 0 : CALL pw_transfer(rho_smooth_r(ispin), tmp_g)
701 0 : CALL pw_transfer(tmp_g, rho_smooth_model_g(ispin))
702 0 : CALL auxbas_pw_pool%give_back_pw(tmp_g)
703 : END IF
704 0 : IF (tau_g_valid) THEN
705 0 : CALL pw_transfer(tau_smooth_g(ispin), tau_smooth_model_g(ispin))
706 : ELSE
707 0 : CALL auxbas_pw_pool%create_pw(tmp_g)
708 0 : CALL pw_transfer(tau_smooth_r(ispin), tmp_g)
709 0 : CALL pw_transfer(tmp_g, tau_smooth_model_g(ispin))
710 0 : CALL auxbas_pw_pool%give_back_pw(tmp_g)
711 : END IF
712 : ELSE
713 0 : IF (rho_g_valid) THEN
714 0 : CALL pw_copy(rho_smooth_g(ispin), rho_smooth_model_g(ispin))
715 : ELSE
716 0 : CALL pw_transfer(rho_smooth_r(ispin), rho_smooth_model_g(ispin))
717 : END IF
718 0 : IF (tau_g_valid) THEN
719 0 : CALL pw_copy(tau_smooth_g(ispin), tau_smooth_model_g(ispin))
720 : ELSE
721 0 : CALL pw_transfer(tau_smooth_r(ispin), tau_smooth_model_g(ispin))
722 : END IF
723 : END IF
724 0 : CALL pw_transfer(rho_smooth_model_g(ispin), rho_smooth_model_r(ispin))
725 0 : CALL pw_transfer(tau_smooth_model_g(ispin), tau_smooth_model_r(ispin))
726 0 : IF (ASSOCIATED(rho_nlcc_use)) THEN
727 0 : CALL pw_axpy(rho_nlcc_use, rho_smooth_model_r(ispin), 1.0_dp)
728 0 : 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 0 : diagnostics_pw_pool, particle_set, cell, para_env)
735 0 : DO ispin = 1, mspin
736 0 : CALL diagnostics_pw_pool%give_back_pw(rho_smooth_model_g(ispin))
737 0 : CALL diagnostics_pw_pool%give_back_pw(rho_smooth_model_r(ispin))
738 0 : CALL diagnostics_pw_pool%give_back_pw(tau_smooth_model_g(ispin))
739 0 : CALL diagnostics_pw_pool%give_back_pw(tau_smooth_model_r(ispin))
740 : END DO
741 0 : DEALLOCATE (rho_smooth_model_g, rho_smooth_model_r, &
742 0 : tau_smooth_model_g, tau_smooth_model_r)
743 0 : 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 153953 : IF (vdW_nl) THEN
751 854 : CALL get_ks_env(ks_env=ks_env, para_env=para_env)
752 : ! no SIC functionals allowed
753 854 : CPASSERT(dft_control%sic_method_id == sic_none)
754 : !
755 854 : CALL pw_env_get(pw_env, vdw_pw_pool=vdw_pw_pool)
756 854 : np = vdw_pw_pool%pw_grid%ngpts_local
757 4270 : ALLOCATE (rhonl(np), drhonl(np, 3))
758 854 : CALL prepare_nlvdw_density(rho_r, rhonl, drhonl, vdw_pw_pool, xc_pw_pool)
759 854 : IF (my_just_energy) THEN
760 : CALL calculate_vxc_nlvdw(vxcnl, rhonl, drhonl, edisp, dispersion_env, &
761 6 : my_just_energy, vdw_pw_pool, para_env)
762 : ELSE
763 848 : CALL vdw_pw_pool%create_pw(vxcnl)
764 : CALL calculate_vxc_nlvdw(vxcnl, rhonl, drhonl, edisp, dispersion_env, &
765 848 : my_just_energy, vdw_pw_pool, para_env, virial=virial)
766 848 : CALL return_nlvdw_potential(my_vxc_rho, vxcnl, vdw_pw_pool, xc_pw_pool)
767 848 : CALL vdw_pw_pool%give_back_pw(vxcnl)
768 : END IF
769 854 : DEALLOCATE (rhonl, drhonl)
770 : END IF
771 :
772 : ! remove the nlcc densities (keep stuff in original state)
773 153953 : IF (ASSOCIATED(rho_nlcc_use)) THEN
774 2474 : factor = -1.0_dp
775 4948 : DO ispin = 1, mspin
776 2474 : CALL pw_axpy(rho_nlcc_use, rho_r(ispin), factor)
777 4948 : 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 153953 : IF (.NOT. my_just_energy) THEN
783 142327 : IF (do_adiabatic_rescaling) THEN
784 24 : IF (ASSOCIATED(my_vxc_rho)) THEN
785 62 : DO ispin = 1, SIZE(my_vxc_rho)
786 62 : 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 153953 : IF (my_scaling /= 1.0_dp) THEN
793 0 : exc = exc*my_scaling
794 0 : IF (ASSOCIATED(my_vxc_rho)) THEN
795 0 : DO ispin = 1, SIZE(my_vxc_rho)
796 0 : CALL pw_scale(my_vxc_rho(ispin), my_scaling)
797 : END DO
798 : END IF
799 0 : IF (ASSOCIATED(my_vxc_tau)) THEN
800 0 : DO ispin = 1, SIZE(my_vxc_tau)
801 0 : 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 153953 : IF (ASSOCIATED(my_vxc_rho)) THEN
809 142327 : vxc_rho => my_vxc_rho
810 142327 : NULLIFY (my_vxc_rho)
811 : END IF
812 153953 : IF (ASSOCIATED(my_vxc_tau)) THEN
813 3726 : vxc_tau => my_vxc_tau
814 3726 : NULLIFY (my_vxc_tau)
815 : END IF
816 153953 : IF (uf_grid) THEN
817 2008 : DO ispin = 1, SIZE(rho_r)
818 2008 : CALL xc_pw_pool%give_back_pw(rho_r(ispin))
819 : END DO
820 1004 : DEALLOCATE (rho_r)
821 1004 : IF (ASSOCIATED(rho_g)) THEN
822 2008 : DO ispin = 1, SIZE(rho_g)
823 2008 : CALL xc_pw_pool%give_back_pw(rho_g(ispin))
824 : END DO
825 1004 : 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 153953 : IF (dft_control%sic_method_id == sic_mauri_spz .AND. .NOT. sic_scaling_b_zero) THEN
831 0 : ALLOCATE (rho_m_rspace(2), rho_m_gspace(2))
832 0 : CALL xc_pw_pool%create_pw(rho_m_gspace(1))
833 0 : CALL xc_pw_pool%create_pw(rho_m_rspace(1))
834 0 : CALL pw_copy(rho_struct_r(1), rho_m_rspace(1))
835 0 : CALL pw_axpy(rho_struct_r(2), rho_m_rspace(1), alpha=-1._dp)
836 0 : CALL pw_copy(rho_struct_g(1), rho_m_gspace(1))
837 0 : CALL pw_axpy(rho_struct_g(2), rho_m_gspace(1), alpha=-1._dp)
838 : ! bit sad, these will be just zero...
839 0 : CALL xc_pw_pool%create_pw(rho_m_gspace(2))
840 0 : CALL xc_pw_pool%create_pw(rho_m_rspace(2))
841 0 : CALL pw_zero(rho_m_rspace(2))
842 0 : CALL pw_zero(rho_m_gspace(2))
843 :
844 0 : 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 0 : weights=weights_use, pw_pool=xc_pw_pool)
848 : ELSE
849 : ! virial untested
850 0 : 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 0 : virial_xc=virial_xc_tmp)
857 : END IF
858 :
859 0 : 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 0 : IF (.NOT. my_just_energy) THEN
863 0 : CALL pw_axpy(my_vxc_rho(1), vxc_rho(1), -dft_control%sic_scaling_b)
864 0 : CALL pw_axpy(my_vxc_rho(1), vxc_rho(2), dft_control%sic_scaling_b)
865 0 : CALL my_vxc_rho(1)%release()
866 0 : CALL my_vxc_rho(2)%release()
867 0 : DEALLOCATE (my_vxc_rho)
868 : END IF
869 :
870 0 : DO ispin = 1, 2
871 0 : CALL xc_pw_pool%give_back_pw(rho_m_rspace(ispin))
872 0 : CALL xc_pw_pool%give_back_pw(rho_m_gspace(ispin))
873 : END DO
874 0 : DEALLOCATE (rho_m_rspace)
875 0 : DEALLOCATE (rho_m_gspace)
876 :
877 : END IF
878 :
879 : ! now we have - sum_s N_s * Exc(rho_s/N_s,0)
880 153953 : 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 0 : CALL get_ks_env(ks_env, nelectron_spin=nelec_spin)
884 :
885 0 : ALLOCATE (rho_m_rspace(2), rho_m_gspace(2))
886 0 : DO ispin = 1, 2
887 0 : CALL xc_pw_pool%create_pw(rho_m_gspace(ispin))
888 0 : CALL xc_pw_pool%create_pw(rho_m_rspace(ispin))
889 : END DO
890 :
891 0 : DO ispin = 1, 2
892 0 : IF (nelec_spin(ispin) > 0.0_dp) THEN
893 0 : 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 0 : nelec_s_inv = 0.0_dp
897 : END IF
898 0 : CALL pw_copy(rho_struct_r(ispin), rho_m_rspace(1))
899 0 : CALL pw_copy(rho_struct_g(ispin), rho_m_gspace(1))
900 0 : CALL pw_scale(rho_m_rspace(1), nelec_s_inv)
901 0 : CALL pw_scale(rho_m_gspace(1), nelec_s_inv)
902 0 : CALL pw_zero(rho_m_rspace(2))
903 0 : CALL pw_zero(rho_m_gspace(2))
904 :
905 0 : 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 0 : weights=weights_use, pw_pool=xc_pw_pool)
909 : ELSE
910 : ! virial untested
911 0 : 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 0 : virial_xc=virial_xc_tmp)
918 : END IF
919 :
920 0 : 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 0 : IF (.NOT. my_just_energy) THEN
924 0 : CALL pw_axpy(my_vxc_rho(1), vxc_rho(ispin), -dft_control%sic_scaling_b)
925 0 : CALL my_vxc_rho(1)%release()
926 0 : CALL my_vxc_rho(2)%release()
927 0 : DEALLOCATE (my_vxc_rho)
928 : END IF
929 : END DO
930 :
931 0 : DO ispin = 1, 2
932 0 : CALL xc_pw_pool%give_back_pw(rho_m_rspace(ispin))
933 0 : CALL xc_pw_pool%give_back_pw(rho_m_gspace(ispin))
934 : END DO
935 0 : DEALLOCATE (rho_m_rspace)
936 0 : DEALLOCATE (rho_m_gspace)
937 :
938 : END IF
939 :
940 : ! compute again the xc but now for Exc(n_down,n_down)
941 153953 : IF (dft_control%sic_method_id == sic_mauri_us .AND. .NOT. sic_scaling_b_zero) THEN
942 0 : ALLOCATE (rho_r(2))
943 0 : rho_r(1) = rho_struct_r(2)
944 0 : rho_r(2) = rho_struct_r(2)
945 0 : IF (rho_g_valid) THEN
946 0 : ALLOCATE (rho_g(2))
947 0 : rho_g(1) = rho_struct_g(2)
948 0 : rho_g(2) = rho_struct_g(2)
949 : END IF
950 :
951 0 : 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 0 : weights=weights_use, pw_pool=xc_pw_pool)
955 : ELSE
956 : ! virial untested
957 0 : 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 0 : virial_xc=virial_xc_tmp)
964 : END IF
965 :
966 0 : exc = exc + dft_control%sic_scaling_b*exc_m
967 :
968 : ! and take care of the potential
969 0 : IF (.NOT. my_just_energy) THEN
970 : ! both go to minority spin
971 0 : CALL pw_axpy(my_vxc_rho(1), vxc_rho(2), 2.0_dp*dft_control%sic_scaling_b)
972 0 : CALL my_vxc_rho(1)%release()
973 0 : CALL my_vxc_rho(2)%release()
974 0 : DEALLOCATE (my_vxc_rho)
975 : END IF
976 0 : DEALLOCATE (rho_r, rho_g)
977 :
978 : END IF
979 :
980 : !
981 : ! cleanups
982 : !
983 153953 : 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 1004 : CALL xc_pw_pool%create_pw(tmp_g)
988 1004 : CALL auxbas_pw_pool%create_pw(tmp_g2)
989 1004 : IF (ASSOCIATED(vxc_rho)) THEN
990 2008 : DO ispin = 1, SIZE(vxc_rho)
991 1004 : CALL auxbas_pw_pool%create_pw(tmp_pw)
992 1004 : CALL pw_transfer(vxc_rho(ispin), tmp_g)
993 1004 : CALL pw_transfer(tmp_g, tmp_g2)
994 1004 : CALL pw_transfer(tmp_g2, tmp_pw)
995 1004 : CALL xc_pw_pool%give_back_pw(vxc_rho(ispin))
996 2008 : vxc_rho(ispin) = tmp_pw
997 : END DO
998 : END IF
999 1004 : IF (ASSOCIATED(vxc_tau)) THEN
1000 444 : DO ispin = 1, SIZE(vxc_tau)
1001 222 : CALL auxbas_pw_pool%create_pw(tmp_pw)
1002 222 : CALL pw_transfer(vxc_tau(ispin), tmp_g)
1003 222 : CALL pw_transfer(tmp_g, tmp_g2)
1004 222 : CALL pw_transfer(tmp_g2, tmp_pw)
1005 222 : CALL xc_pw_pool%give_back_pw(vxc_tau(ispin))
1006 444 : vxc_tau(ispin) = tmp_pw
1007 : END DO
1008 : END IF
1009 1004 : CALL auxbas_pw_pool%give_back_pw(tmp_g2)
1010 2008 : CALL xc_pw_pool%give_back_pw(tmp_g)
1011 : END BLOCK
1012 : END IF
1013 153953 : IF (ASSOCIATED(tau) .AND. uf_grid) THEN
1014 444 : DO ispin = 1, SIZE(tau)
1015 444 : CALL xc_pw_pool%give_back_pw(tau(ispin))
1016 : END DO
1017 222 : DEALLOCATE (tau)
1018 : END IF
1019 153953 : IF (ASSOCIATED(weights_xc)) THEN
1020 990 : CALL xc_pw_pool%give_back_pw(weights_xc)
1021 990 : DEALLOCATE (weights_xc)
1022 : END IF
1023 153953 : IF (ASSOCIATED(rho_nlcc_xc)) THEN
1024 56 : CALL xc_pw_pool%give_back_pw(rho_nlcc_xc)
1025 56 : DEALLOCATE (rho_nlcc_xc)
1026 : END IF
1027 153953 : IF (ASSOCIATED(rho_nlcc_g_xc)) THEN
1028 56 : CALL xc_pw_pool%give_back_pw(rho_nlcc_g_xc)
1029 56 : DEALLOCATE (rho_nlcc_g_xc)
1030 : END IF
1031 :
1032 153953 : IF (ASSOCIATED(rho_hard_r)) THEN
1033 8 : DO ispin = 1, SIZE(rho_hard_r)
1034 4 : CALL auxbas_pw_pool%give_back_pw(rho_hard_r(ispin))
1035 4 : CALL auxbas_pw_pool%give_back_pw(rho_hard_g(ispin))
1036 4 : CALL auxbas_pw_pool%give_back_pw(tau_hard_r(ispin))
1037 8 : CALL auxbas_pw_pool%give_back_pw(tau_hard_g(ispin))
1038 : END DO
1039 4 : DEALLOCATE (rho_hard_r, rho_hard_g, tau_hard_r, tau_hard_g)
1040 : END IF
1041 :
1042 : END IF
1043 :
1044 169425 : CALL timestop(handle)
1045 :
1046 338850 : 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 0 : 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 0 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_test_g, tau_test_g, vxc_rho_g, &
1088 0 : vxc_tau_g
1089 0 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_test_r, tau_test, vxc_dummy_rho, &
1090 0 : vxc_dummy_tau
1091 :
1092 0 : NULLIFY (rho_test_g, rho_test_r, tau_test, tau_test_g, vxc_dummy_rho, vxc_dummy_tau, &
1093 0 : vxc_rho_g, vxc_tau_g)
1094 0 : nspin = SIZE(rho_combined_r)
1095 0 : CPASSERT(SIZE(rho_combined_g) == nspin)
1096 0 : CPASSERT(SIZE(tau_combined) == nspin)
1097 0 : CPASSERT(SIZE(rho_smooth_r) == nspin)
1098 0 : CPASSERT(SIZE(rho_smooth_g) == nspin)
1099 0 : CPASSERT(SIZE(tau_smooth) == nspin)
1100 0 : CPASSERT(SIZE(vxc_rho) == nspin)
1101 0 : CPASSERT(SIZE(vxc_tau) == nspin)
1102 :
1103 : ALLOCATE (rho_test_g(nspin), rho_test_r(nspin), tau_test(nspin), tau_test_g(nspin), &
1104 0 : vxc_rho_g(nspin), vxc_tau_g(nspin))
1105 0 : DO ispin = 1, nspin
1106 0 : CALL pw_pool%create_pw(rho_test_g(ispin))
1107 0 : CALL pw_pool%create_pw(rho_test_r(ispin))
1108 0 : CALL pw_pool%create_pw(tau_test(ispin))
1109 0 : CALL pw_pool%create_pw(tau_test_g(ispin))
1110 0 : CALL pw_pool%create_pw(vxc_rho_g(ispin))
1111 0 : CALL pw_pool%create_pw(vxc_tau_g(ispin))
1112 0 : CALL pw_copy(rho_combined_g(ispin), rho_test_g(ispin))
1113 0 : CALL pw_axpy(rho_smooth_g(ispin), rho_test_g(ispin), -1.0_dp)
1114 0 : CALL pw_copy(rho_combined_r(ispin), rho_test_r(ispin))
1115 0 : CALL pw_axpy(rho_smooth_r(ispin), rho_test_r(ispin), -1.0_dp)
1116 0 : CALL pw_copy(tau_combined(ispin), tau_test(ispin))
1117 0 : CALL pw_axpy(tau_smooth(ispin), tau_test(ispin), -1.0_dp)
1118 0 : CALL pw_transfer(tau_test(ispin), tau_test_g(ispin))
1119 0 : CALL pw_transfer(vxc_rho(ispin), vxc_rho_g(ispin))
1120 0 : CALL pw_transfer(vxc_tau(ispin), vxc_tau_g(ispin))
1121 : END DO
1122 0 : rho_contraction = 0.0_dp
1123 0 : rho_g_contraction = 0.0_dp
1124 0 : tau_contraction = 0.0_dp
1125 0 : tau_g_contraction = 0.0_dp
1126 0 : DO ispin = 1, nspin
1127 0 : rho_contraction = rho_contraction + pw_integral_ab(rho_test_r(ispin), vxc_rho(ispin))
1128 0 : rho_g_contraction = rho_g_contraction + pw_integral_ab(rho_test_g(ispin), vxc_rho_g(ispin))
1129 0 : tau_contraction = tau_contraction + pw_integral_ab(tau_test(ispin), vxc_tau(ispin))
1130 0 : tau_g_contraction = tau_g_contraction + pw_integral_ab(tau_test_g(ispin), vxc_tau_g(ispin))
1131 : END DO
1132 :
1133 0 : total_contraction = rho_contraction + tau_contraction
1134 0 : DO istep = 1, SIZE(fd_steps)
1135 0 : fd_step = fd_steps(istep)
1136 0 : DO ispin = 1, nspin
1137 0 : CALL pw_copy(rho_combined_g(ispin), rho_test_g(ispin))
1138 0 : CALL pw_axpy(rho_combined_g(ispin), rho_test_g(ispin), fd_step)
1139 0 : CALL pw_axpy(rho_smooth_g(ispin), rho_test_g(ispin), -fd_step)
1140 0 : CALL pw_copy(rho_combined_r(ispin), rho_test_r(ispin))
1141 0 : CALL pw_axpy(rho_combined_r(ispin), rho_test_r(ispin), fd_step)
1142 0 : CALL pw_axpy(rho_smooth_r(ispin), rho_test_r(ispin), -fd_step)
1143 0 : CALL pw_copy(tau_combined(ispin), tau_test(ispin))
1144 0 : CALL pw_axpy(tau_combined(ispin), tau_test(ispin), fd_step)
1145 0 : 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 0 : just_energy=.TRUE.)
1152 0 : CPASSERT(.NOT. ASSOCIATED(vxc_dummy_rho))
1153 0 : CPASSERT(.NOT. ASSOCIATED(vxc_dummy_tau))
1154 :
1155 0 : DO ispin = 1, nspin
1156 0 : CALL pw_copy(rho_combined_g(ispin), rho_test_g(ispin))
1157 0 : CALL pw_axpy(rho_combined_g(ispin), rho_test_g(ispin), -fd_step)
1158 0 : CALL pw_axpy(rho_smooth_g(ispin), rho_test_g(ispin), fd_step)
1159 0 : CALL pw_copy(rho_combined_r(ispin), rho_test_r(ispin))
1160 0 : CALL pw_axpy(rho_combined_r(ispin), rho_test_r(ispin), -fd_step)
1161 0 : CALL pw_axpy(rho_smooth_r(ispin), rho_test_r(ispin), fd_step)
1162 0 : CALL pw_copy(tau_combined(ispin), tau_test(ispin))
1163 0 : CALL pw_axpy(tau_combined(ispin), tau_test(ispin), -fd_step)
1164 0 : 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 0 : just_energy=.TRUE.)
1171 0 : CPASSERT(.NOT. ASSOCIATED(vxc_dummy_rho))
1172 0 : CPASSERT(.NOT. ASSOCIATED(vxc_dummy_tau))
1173 :
1174 0 : fd_derivative = (exc_plus - exc_minus)/(2.0_dp*fd_step)
1175 0 : IF (para_env%mepos == 0) THEN
1176 0 : output_unit = cp_logger_get_default_io_unit()
1177 0 : IF (output_unit > 0) THEN
1178 : WRITE (output_unit, "(/,T2,A,1X,ES12.4,2(1X,ES20.10))") &
1179 0 : "SKALA_GPW| Composite-direction FD step and energies", &
1180 0 : fd_step, exc_plus, exc_minus
1181 : WRITE (output_unit, "(T2,A,1X,ES20.10)") &
1182 0 : "SKALA_GPW| Composite-direction energy FD", fd_derivative
1183 : WRITE (output_unit, "(T2,A,1X,ES20.10)") &
1184 0 : "SKALA_GPW| Composite-direction FD/VXC difference", &
1185 0 : fd_derivative - total_contraction
1186 : END IF
1187 : END IF
1188 : END DO
1189 0 : IF (para_env%mepos == 0) THEN
1190 0 : output_unit = cp_logger_get_default_io_unit()
1191 0 : IF (output_unit > 0) THEN
1192 : WRITE (output_unit, "(T2,A,1X,ES20.10)") &
1193 0 : "SKALA_GPW| Composite-direction rho VXC contraction", rho_contraction
1194 : WRITE (output_unit, "(T2,A,1X,ES20.10)") &
1195 0 : "SKALA_GPW| Composite-direction rho G-space contraction", rho_g_contraction
1196 : WRITE (output_unit, "(T2,A,1X,ES20.10)") &
1197 0 : "SKALA_GPW| Composite-direction tau VXC contraction", tau_contraction
1198 : WRITE (output_unit, "(T2,A,1X,ES20.10)") &
1199 0 : "SKALA_GPW| Composite-direction tau G-space contraction", tau_g_contraction
1200 : WRITE (output_unit, "(T2,A,1X,ES20.10)") &
1201 0 : "SKALA_GPW| Composite-direction total VXC contraction", total_contraction
1202 : END IF
1203 : END IF
1204 :
1205 0 : DO ispin = 1, nspin
1206 0 : CALL pw_pool%give_back_pw(rho_test_g(ispin))
1207 0 : CALL pw_pool%give_back_pw(rho_test_r(ispin))
1208 0 : CALL pw_pool%give_back_pw(tau_test(ispin))
1209 0 : CALL pw_pool%give_back_pw(tau_test_g(ispin))
1210 0 : CALL pw_pool%give_back_pw(vxc_rho_g(ispin))
1211 0 : CALL pw_pool%give_back_pw(vxc_tau_g(ispin))
1212 : END DO
1213 0 : DEALLOCATE (rho_test_g, rho_test_r, tau_test, tau_test_g, vxc_rho_g, vxc_tau_g)
1214 :
1215 0 : 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 100 : SUBROUTINE qs_xc_density(ks_env, rho_struct, xc_section, dispersion_env, &
1231 100 : 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 100 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: rhonl
1249 100 : 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 100 : 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 100 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, rho_struct_r, tau_r, &
1260 100 : 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 100 : CALL timeset(routineN, handle)
1265 :
1266 100 : NULLIFY (dft_control, pw_env, auxbas_pw_pool, xc_pw_pool, vdw_pw_pool, cell, &
1267 100 : rho_g, rho_struct_g, tau_g, tau_struct_g, rho_nlcc, rho_nlcc_g, &
1268 100 : rho_nlcc_g_use, rho_nlcc_g_xc, rho_nlcc_use, rho_nlcc_xc, rho_r, &
1269 100 : rho_struct_r, tau_r, tau_struct_r, vxc_rho, vxc_tau, weights, &
1270 100 : 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 100 : 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 100 : tau_g=tau_struct_g)
1288 100 : nspins = dft_control%nspins
1289 100 : mspin = SIZE(rho_struct_r)
1290 100 : rho_r => rho_struct_r
1291 100 : rho_g => rho_struct_g
1292 100 : tau_r => tau_struct_r
1293 100 : tau_g => tau_struct_g
1294 100 : rho_nlcc_use => rho_nlcc
1295 100 : rho_nlcc_g_use => rho_nlcc_g
1296 100 : weights_use => weights
1297 :
1298 100 : CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", i_val=myfun)
1299 100 : CALL section_vals_val_get(xc_section, "VDW_POTENTIAL%POTENTIAL_TYPE", i_val=vdw)
1300 100 : vdW_nl = (vdw == xc_vdw_fun_nonloc)
1301 100 : IF (PRESENT(xc_ener)) THEN
1302 34 : IF (tau_r_valid) THEN
1303 0 : CALL cp_warn(__LOCATION__, "Tau contribution will not be correctly handled")
1304 : END IF
1305 : END IF
1306 100 : IF (vdW_nl) THEN
1307 0 : CALL cp_warn(__LOCATION__, "vdW functional contribution will be ignored")
1308 : END IF
1309 :
1310 100 : CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool, auxbas_pw_pool=auxbas_pw_pool)
1311 100 : uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid, xc_pw_pool%pw_grid)
1312 :
1313 100 : IF (PRESENT(xc_ener)) THEN
1314 34 : CALL pw_zero(xc_ener)
1315 : END IF
1316 100 : IF (PRESENT(xc_den)) THEN
1317 66 : CALL pw_zero(xc_den)
1318 : END IF
1319 100 : IF (PRESENT(exc)) THEN
1320 0 : CALL pw_zero(exc)
1321 : END IF
1322 100 : IF (PRESENT(vxc)) THEN
1323 138 : DO ispin = 1, nspins
1324 138 : CALL pw_zero(vxc(ispin))
1325 : END DO
1326 : END IF
1327 100 : IF (PRESENT(vtau)) THEN
1328 40 : DO ispin = 1, nspins
1329 40 : CALL pw_zero(vtau(ispin))
1330 : END DO
1331 : END IF
1332 :
1333 100 : IF (myfun /= xc_none) THEN
1334 :
1335 98 : CPASSERT(ASSOCIATED(rho_struct))
1336 98 : CPASSERT(dft_control%sic_method_id == sic_none)
1337 :
1338 98 : IF (uf_grid) THEN
1339 2 : NULLIFY (rho_r, rho_g, tau_r, tau_g)
1340 2 : IF (rho_g_valid) THEN
1341 2 : CALL create_density_on_pool(xc_pw_pool, rho_struct_g, rho_r, rho_g)
1342 0 : ELSE IF (ASSOCIATED(rho_struct_r)) THEN
1343 0 : CALL create_density_on_pool_from_r(auxbas_pw_pool, xc_pw_pool, rho_struct_r, rho_r, rho_g)
1344 : ELSE
1345 0 : CPABORT("Fine Grid in qs_xc_density requires rho_r or rho_g")
1346 : END IF
1347 2 : IF (tau_r_valid) THEN
1348 0 : IF (tau_g_valid) THEN
1349 0 : CALL create_density_on_pool(xc_pw_pool, tau_struct_g, tau_r, tau_g)
1350 0 : ELSE IF (ASSOCIATED(tau_struct_r)) THEN
1351 0 : CALL create_density_on_pool_from_r(auxbas_pw_pool, xc_pw_pool, tau_struct_r, tau_r, tau_g)
1352 : ELSE
1353 0 : CPABORT("Fine Grid in qs_xc_density requires tau_r or tau_g")
1354 : END IF
1355 : END IF
1356 2 : IF (ASSOCIATED(weights)) THEN
1357 2 : ALLOCATE (weights_xc)
1358 2 : CALL xc_pw_pool%create_pw(weights_xc)
1359 2 : CALL transfer_rspace_between_pools(auxbas_pw_pool, xc_pw_pool, weights, weights_xc)
1360 2 : weights_use => weights_xc
1361 : END IF
1362 2 : IF (ASSOCIATED(rho_nlcc)) THEN
1363 0 : CPASSERT(ASSOCIATED(rho_nlcc_g))
1364 0 : ALLOCATE (rho_nlcc_g_xc, rho_nlcc_xc)
1365 0 : CALL xc_pw_pool%create_pw(rho_nlcc_g_xc)
1366 0 : CALL xc_pw_pool%create_pw(rho_nlcc_xc)
1367 0 : CALL pw_transfer(rho_nlcc_g, rho_nlcc_g_xc)
1368 0 : 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 98 : IF (ASSOCIATED(rho_nlcc_use)) THEN
1376 0 : factor = 1.0_dp
1377 0 : DO ispin = 1, mspin
1378 0 : CALL pw_axpy(rho_nlcc_use, rho_r(ispin), factor)
1379 0 : CALL pw_axpy(rho_nlcc_g_use, rho_g(ispin), factor)
1380 : END DO
1381 : END IF
1382 98 : 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 98 : 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 98 : IF (vdW_nl) THEN
1395 0 : CALL get_ks_env(ks_env=ks_env, para_env=para_env)
1396 : ! no SIC functionals allowed
1397 0 : CPASSERT(dft_control%sic_method_id == sic_none)
1398 : !
1399 0 : CALL pw_env_get(pw_env, vdw_pw_pool=vdw_pw_pool)
1400 0 : np = vdw_pw_pool%pw_grid%ngpts_local
1401 0 : ALLOCATE (rhonl(np), drhonl(np, 3))
1402 0 : CALL prepare_nlvdw_density(rho_r, rhonl, drhonl, vdw_pw_pool, xc_pw_pool)
1403 0 : CALL vdw_pw_pool%create_pw(vxcnl)
1404 : CALL calculate_vxc_nlvdw(vxcnl, rhonl, drhonl, edisp, dispersion_env, &
1405 0 : .FALSE., vdw_pw_pool, para_env)
1406 0 : CALL return_nlvdw_potential(vxc_rho, vxcnl, vdw_pw_pool, xc_pw_pool)
1407 0 : CALL vdw_pw_pool%give_back_pw(vxcnl)
1408 0 : DEALLOCATE (rhonl, drhonl)
1409 : END IF
1410 :
1411 : ! remove the nlcc densities (keep stuff in original state)
1412 98 : IF (ASSOCIATED(rho_nlcc_use)) THEN
1413 0 : factor = -1.0_dp
1414 0 : DO ispin = 1, mspin
1415 0 : CALL pw_axpy(rho_nlcc_use, rho_r(ispin), factor)
1416 0 : CALL pw_axpy(rho_nlcc_g_use, rho_g(ispin), factor)
1417 : END DO
1418 : END IF
1419 : !
1420 98 : IF (PRESENT(xc_den)) THEN
1421 64 : rho_cutoff = 1.E-14_dp
1422 64 : IF (uf_grid) THEN
1423 : BLOCK
1424 : TYPE(pw_r3d_rs_type) :: tmp_pw
1425 0 : CALL xc_pw_pool%create_pw(tmp_pw)
1426 0 : CALL pw_copy(exc_r, tmp_pw)
1427 0 : CALL calc_xc_density(tmp_pw, rho_r, rho_cutoff)
1428 0 : CALL transfer_rspace_between_pools(xc_pw_pool, auxbas_pw_pool, tmp_pw, xc_den)
1429 0 : CALL xc_pw_pool%give_back_pw(tmp_pw)
1430 : END BLOCK
1431 : ELSE
1432 64 : CALL pw_copy(exc_r, xc_den)
1433 64 : CALL calc_xc_density(xc_den, rho_r, rho_cutoff)
1434 : END IF
1435 : END IF
1436 98 : IF (PRESENT(xc_ener)) THEN
1437 34 : IF (uf_grid) THEN
1438 : BLOCK
1439 : TYPE(pw_r3d_rs_type) :: tmp_pw
1440 2 : CALL xc_pw_pool%create_pw(tmp_pw)
1441 2 : CALL pw_copy(exc_r, tmp_pw)
1442 4 : DO ispin = 1, nspins
1443 4 : CALL pw_multiply(tmp_pw, vxc_rho(ispin), rho_r(ispin), alpha=-1.0_dp)
1444 : END DO
1445 2 : CALL transfer_rspace_between_pools(xc_pw_pool, auxbas_pw_pool, tmp_pw, xc_ener)
1446 2 : CALL xc_pw_pool%give_back_pw(tmp_pw)
1447 : END BLOCK
1448 : ELSE
1449 32 : CALL pw_copy(exc_r, xc_ener)
1450 64 : DO ispin = 1, nspins
1451 64 : 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 98 : IF (PRESENT(exc)) THEN
1456 0 : IF (uf_grid) THEN
1457 0 : CALL transfer_rspace_between_pools(xc_pw_pool, auxbas_pw_pool, exc_r, exc)
1458 : ELSE
1459 0 : CALL pw_copy(exc_r, exc)
1460 : END IF
1461 : END IF
1462 98 : IF (PRESENT(vxc)) THEN
1463 134 : DO ispin = 1, nspins
1464 134 : IF (uf_grid) THEN
1465 0 : CALL transfer_rspace_between_pools(xc_pw_pool, auxbas_pw_pool, vxc_rho(ispin), vxc(ispin))
1466 : ELSE
1467 70 : CALL pw_copy(vxc_rho(ispin), vxc(ispin))
1468 : END IF
1469 : END DO
1470 : END IF
1471 98 : IF (PRESENT(vtau) .AND. ASSOCIATED(vxc_tau)) THEN
1472 40 : DO ispin = 1, nspins
1473 40 : IF (uf_grid) THEN
1474 0 : CALL transfer_rspace_between_pools(xc_pw_pool, auxbas_pw_pool, vxc_tau(ispin), vtau(ispin))
1475 : ELSE
1476 20 : CALL pw_copy(vxc_tau(ispin), vtau(ispin))
1477 : END IF
1478 : END DO
1479 : END IF
1480 : ! remove arrays
1481 98 : IF (ASSOCIATED(vxc_rho)) THEN
1482 202 : DO ispin = 1, nspins
1483 202 : CALL vxc_rho(ispin)%release()
1484 : END DO
1485 98 : DEALLOCATE (vxc_rho)
1486 : END IF
1487 98 : IF (ASSOCIATED(vxc_tau)) THEN
1488 40 : DO ispin = 1, nspins
1489 40 : CALL vxc_tau(ispin)%release()
1490 : END DO
1491 20 : DEALLOCATE (vxc_tau)
1492 : END IF
1493 98 : CALL exc_r%release()
1494 98 : IF (uf_grid) THEN
1495 2 : CALL give_back_density_on_pool(xc_pw_pool, rho_r, rho_g)
1496 2 : IF (ASSOCIATED(tau_r)) CALL give_back_density_on_pool(xc_pw_pool, tau_r, tau_g)
1497 2 : IF (ASSOCIATED(weights_xc)) THEN
1498 2 : CALL xc_pw_pool%give_back_pw(weights_xc)
1499 2 : DEALLOCATE (weights_xc)
1500 : END IF
1501 2 : IF (ASSOCIATED(rho_nlcc_xc)) THEN
1502 0 : CALL xc_pw_pool%give_back_pw(rho_nlcc_xc)
1503 0 : DEALLOCATE (rho_nlcc_xc)
1504 : END IF
1505 2 : IF (ASSOCIATED(rho_nlcc_g_xc)) THEN
1506 0 : CALL xc_pw_pool%give_back_pw(rho_nlcc_g_xc)
1507 0 : DEALLOCATE (rho_nlcc_g_xc)
1508 : END IF
1509 : END IF
1510 : !
1511 : END IF
1512 :
1513 100 : CALL timestop(handle)
1514 :
1515 200 : 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 4 : 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 0 : CPASSERT(ASSOCIATED(source_pw_pool))
1531 4 : CPASSERT(ASSOCIATED(target_pw_pool))
1532 :
1533 4 : IF (pw_grid_compare(source_pw_pool%pw_grid, target_pw_pool%pw_grid)) THEN
1534 0 : CALL pw_copy(source, TARGET)
1535 : ELSE
1536 4 : CALL source_pw_pool%create_pw(source_g)
1537 4 : CALL target_pw_pool%create_pw(target_g)
1538 4 : CALL pw_transfer(source, source_g)
1539 4 : CALL pw_transfer(source_g, target_g)
1540 4 : CALL pw_transfer(target_g, TARGET)
1541 4 : CALL target_pw_pool%give_back_pw(target_g)
1542 4 : CALL source_pw_pool%give_back_pw(source_g)
1543 : END IF
1544 :
1545 4 : 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 2 : 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 2 : CPASSERT(ASSOCIATED(pw_pool))
1563 2 : CPASSERT(ASSOCIATED(rho_g_in))
1564 :
1565 2 : nspins = SIZE(rho_g_in)
1566 14 : ALLOCATE (rho_r_out(nspins), rho_g_out(nspins))
1567 4 : DO ispin = 1, nspins
1568 2 : CALL pw_pool%create_pw(rho_g_out(ispin))
1569 2 : CALL pw_pool%create_pw(rho_r_out(ispin))
1570 2 : CALL pw_transfer(rho_g_in(ispin), rho_g_out(ispin))
1571 4 : CALL pw_transfer(rho_g_out(ispin), rho_r_out(ispin))
1572 : END DO
1573 :
1574 2 : 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 0 : 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 0 : CPASSERT(ASSOCIATED(source_pw_pool))
1593 0 : CPASSERT(ASSOCIATED(target_pw_pool))
1594 0 : CPASSERT(ASSOCIATED(rho_r_in))
1595 :
1596 0 : nspins = SIZE(rho_r_in)
1597 0 : ALLOCATE (rho_r_out(nspins), rho_g_out(nspins))
1598 0 : DO ispin = 1, nspins
1599 0 : CALL source_pw_pool%create_pw(rho_g_in)
1600 0 : CALL target_pw_pool%create_pw(rho_g_out(ispin))
1601 0 : CALL target_pw_pool%create_pw(rho_r_out(ispin))
1602 0 : CALL pw_transfer(rho_r_in(ispin), rho_g_in)
1603 0 : CALL pw_transfer(rho_g_in, rho_g_out(ispin))
1604 0 : CALL pw_transfer(rho_g_out(ispin), rho_r_out(ispin))
1605 0 : CALL source_pw_pool%give_back_pw(rho_g_in)
1606 : END DO
1607 :
1608 0 : 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 2 : 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 2 : CPASSERT(ASSOCIATED(pw_pool))
1624 :
1625 2 : IF (ASSOCIATED(rho_r)) THEN
1626 4 : DO ispin = 1, SIZE(rho_r)
1627 4 : CALL pw_pool%give_back_pw(rho_r(ispin))
1628 : END DO
1629 2 : DEALLOCATE (rho_r)
1630 : END IF
1631 2 : IF (ASSOCIATED(rho_g)) THEN
1632 4 : DO ispin = 1, SIZE(rho_g)
1633 4 : CALL pw_pool%give_back_pw(rho_g(ispin))
1634 : END DO
1635 2 : DEALLOCATE (rho_g)
1636 : END IF
1637 :
1638 2 : END SUBROUTINE give_back_density_on_pool
1639 :
1640 : END MODULE qs_vxc
|
