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 input_constants, ONLY: sic_ad,&
21 : sic_eo,&
22 : sic_mauri_spz,&
23 : sic_mauri_us,&
24 : sic_none,&
25 : xc_none,&
26 : xc_vdw_fun_nonloc
27 : USE input_section_types, ONLY: section_vals_type,&
28 : section_vals_val_get
29 : USE kinds, ONLY: dp
30 : USE message_passing, ONLY: mp_para_env_type
31 : USE particle_types, ONLY: particle_type
32 : USE pw_env_types, ONLY: pw_env_get,&
33 : pw_env_type
34 : USE pw_grids, ONLY: pw_grid_compare
35 : USE pw_methods, ONLY: pw_axpy,&
36 : pw_copy,&
37 : pw_multiply,&
38 : pw_scale,&
39 : pw_transfer,&
40 : pw_zero
41 : USE pw_pool_types, ONLY: pw_pool_type
42 : USE pw_types, ONLY: pw_c1d_gs_type,&
43 : pw_r3d_rs_type
44 : USE qs_dispersion_nonloc, ONLY: calculate_dispersion_nonloc
45 : USE qs_dispersion_types, ONLY: qs_dispersion_type
46 : USE qs_ks_types, ONLY: get_ks_env,&
47 : qs_ks_env_type
48 : USE qs_rho_types, ONLY: qs_rho_get,&
49 : qs_rho_type
50 : USE skala_gpw_functional, ONLY: skala_gpw_eval,&
51 : xc_section_uses_native_skala_grid
52 : USE virial_types, ONLY: virial_type
53 : USE xc, ONLY: calc_xc_density,&
54 : xc_exc_calc,&
55 : xc_vxc_pw_create
56 : #include "./base/base_uses.f90"
57 :
58 : IMPLICIT NONE
59 :
60 : PRIVATE
61 :
62 : ! *** Public subroutines ***
63 : PUBLIC :: qs_vxc_create, qs_xc_density
64 :
65 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_vxc'
66 :
67 : CONTAINS
68 :
69 : ! **************************************************************************************************
70 : !> \brief calculates and allocates the xc potential, already reducing it to
71 : !> the dependence on rho and the one on tau
72 : !> \param ks_env to get all the needed things
73 : !> \param rho_struct density for which v_xc is calculated
74 : !> \param xc_section ...
75 : !> \param vxc_rho will contain the v_xc part that depend on rho
76 : !> (if one of the chosen xc functionals has it it is allocated and you
77 : !> are responsible for it)
78 : !> \param vxc_tau will contain the kinetic tau part of v_xc
79 : !> (if one of the chosen xc functionals has it it is allocated and you
80 : !> are responsible for it)
81 : !> \param exc ...
82 : !> \param just_energy if true calculates just the energy, and does not
83 : !> allocate v_*_rspace
84 : !> \param edisp ...
85 : !> \param dispersion_env ...
86 : !> \param adiabatic_rescale_factor ...
87 : !> \param pw_env_external external plane wave environment
88 : !> \param native_skala_atom_force ...
89 : !> \par History
90 : !> - 05.2002 modified to use the mp_allgather function each pe
91 : !> computes only part of the grid and this is broadcasted to all
92 : !> instead of summed.
93 : !> This scales significantly better (e.g. factor 3 on 12 cpus
94 : !> 32 H2O) [Joost VdV]
95 : !> - moved to qs_ks_methods [fawzi]
96 : !> - sic alterations [Joost VandeVondele]
97 : !> \author Fawzi Mohamed
98 : ! **************************************************************************************************
99 760365 : SUBROUTINE qs_vxc_create(ks_env, rho_struct, xc_section, vxc_rho, vxc_tau, exc, &
100 : just_energy, edisp, dispersion_env, adiabatic_rescale_factor, &
101 152073 : pw_env_external, native_skala_atom_force)
102 :
103 : TYPE(qs_ks_env_type), POINTER :: ks_env
104 : TYPE(qs_rho_type), POINTER :: rho_struct
105 : TYPE(section_vals_type), POINTER :: xc_section
106 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_rho, vxc_tau
107 : REAL(KIND=dp), INTENT(out) :: exc
108 : LOGICAL, INTENT(in), OPTIONAL :: just_energy
109 : REAL(KIND=dp), INTENT(out), OPTIONAL :: edisp
110 : TYPE(qs_dispersion_type), OPTIONAL, POINTER :: dispersion_env
111 : REAL(KIND=dp), INTENT(in), OPTIONAL :: adiabatic_rescale_factor
112 : TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_external
113 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT), &
114 : OPTIONAL :: native_skala_atom_force
115 :
116 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_vxc_create'
117 :
118 : INTEGER :: handle, ispin, mspin, myfun, &
119 : nelec_spin(2), vdw
120 : LOGICAL :: compute_virial, do_adiabatic_rescaling, my_just_energy, native_skala_grid, &
121 : rho_g_valid, sic_scaling_b_zero, tau_g_valid, tau_r_valid, uf_grid, vdW_nl
122 : REAL(KIND=dp) :: exc_m, factor, &
123 : my_adiabatic_rescale_factor, &
124 : my_scaling, nelec_s_inv
125 : REAL(KIND=dp), DIMENSION(3, 3) :: virial_xc_tmp
126 : TYPE(cell_type), POINTER :: cell
127 : TYPE(dft_control_type), POINTER :: dft_control
128 : TYPE(mp_para_env_type), POINTER :: para_env
129 152073 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
130 152073 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g, rho_m_gspace, rho_struct_g, &
131 152073 : tau_struct_g
132 : TYPE(pw_c1d_gs_type), POINTER :: rho_nlcc_g, rho_nlcc_g_use, &
133 : rho_nlcc_g_xc, tmp_g, tmp_g2
134 : TYPE(pw_env_type), POINTER :: pw_env
135 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool, vdw_pw_pool, xc_pw_pool
136 152073 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: my_vxc_rho, my_vxc_tau, rho_m_rspace, &
137 152073 : rho_r, rho_struct_r, tau, tau_struct_r
138 : TYPE(pw_r3d_rs_type), POINTER :: rho_nlcc, rho_nlcc_use, rho_nlcc_xc, &
139 : tmp_pw, weights, weights_use, &
140 : weights_xc
141 : TYPE(virial_type), POINTER :: virial
142 :
143 152073 : CALL timeset(routineN, handle)
144 :
145 152073 : CPASSERT(.NOT. ASSOCIATED(vxc_rho))
146 152073 : CPASSERT(.NOT. ASSOCIATED(vxc_tau))
147 152073 : NULLIFY (dft_control, pw_env, auxbas_pw_pool, xc_pw_pool, vdw_pw_pool, cell, my_vxc_rho, &
148 152073 : tmp_pw, tmp_g, tmp_g2, my_vxc_tau, rho_g, rho_r, tau, rho_m_rspace, &
149 152073 : rho_m_gspace, rho_nlcc, rho_nlcc_g, rho_nlcc_g_use, rho_nlcc_g_xc, &
150 152073 : rho_nlcc_use, rho_nlcc_xc, rho_struct_r, rho_struct_g, tau_struct_g, tau_struct_r, &
151 152073 : weights_use, weights_xc, particle_set)
152 :
153 152073 : exc = 0.0_dp
154 152073 : my_just_energy = .FALSE.
155 152073 : IF (PRESENT(just_energy)) my_just_energy = just_energy
156 152073 : my_adiabatic_rescale_factor = 1.0_dp
157 152073 : do_adiabatic_rescaling = .FALSE.
158 152073 : IF (PRESENT(adiabatic_rescale_factor)) THEN
159 44 : my_adiabatic_rescale_factor = adiabatic_rescale_factor
160 44 : do_adiabatic_rescaling = .TRUE.
161 : END IF
162 :
163 : CALL get_ks_env(ks_env, &
164 : dft_control=dft_control, &
165 : pw_env=pw_env, &
166 : cell=cell, &
167 : particle_set=particle_set, &
168 : xcint_weights=weights, &
169 : virial=virial, &
170 : rho_nlcc=rho_nlcc, &
171 152073 : rho_nlcc_g=rho_nlcc_g)
172 152073 : rho_nlcc_use => rho_nlcc
173 152073 : rho_nlcc_g_use => rho_nlcc_g
174 152073 : weights_use => weights
175 :
176 : CALL qs_rho_get(rho_struct, &
177 : tau_r_valid=tau_r_valid, &
178 : tau_g_valid=tau_g_valid, &
179 : rho_g_valid=rho_g_valid, &
180 : rho_r=rho_struct_r, &
181 : rho_g=rho_struct_g, &
182 : tau_g=tau_struct_g, &
183 152073 : tau_r=tau_struct_r)
184 :
185 152073 : compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
186 152073 : IF (compute_virial) THEN
187 36426 : virial%pv_xc = 0.0_dp
188 : END IF
189 :
190 : CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", &
191 152073 : i_val=myfun)
192 : CALL section_vals_val_get(xc_section, "VDW_POTENTIAL%POTENTIAL_TYPE", &
193 152073 : i_val=vdw)
194 :
195 152073 : vdW_nl = (vdw == xc_vdw_fun_nonloc)
196 : ! this combination has not been investigated
197 152073 : CPASSERT(.NOT. (do_adiabatic_rescaling .AND. vdW_nl))
198 : ! are the necessary inputs available
199 152073 : IF (.NOT. (PRESENT(dispersion_env) .AND. PRESENT(edisp))) THEN
200 : vdW_nl = .FALSE.
201 : END IF
202 152073 : IF (PRESENT(edisp)) edisp = 0.0_dp
203 152073 : native_skala_grid = xc_section_uses_native_skala_grid(xc_section)
204 152073 : IF (native_skala_grid .AND. ASSOCIATED(rho_nlcc_use)) THEN
205 : CALL cp_abort(__LOCATION__, &
206 : "Native SKALA grid evaluation with NLCC pseudopotentials is not implemented. "// &
207 0 : "The frozen core density would need a SKALA-consistent feature definition.")
208 : END IF
209 :
210 152073 : IF (myfun /= xc_none .OR. vdW_nl) THEN
211 :
212 : ! test if the real space density is available
213 137431 : CPASSERT(ASSOCIATED(rho_struct))
214 137431 : IF (dft_control%nspins /= 1 .AND. dft_control%nspins /= 2) &
215 0 : CPABORT("nspins must be 1 or 2")
216 137431 : mspin = SIZE(rho_struct_r)
217 137431 : IF (dft_control%nspins == 2 .AND. mspin == 1) &
218 0 : CPABORT("Spin count mismatch")
219 :
220 : ! there are some options related to SIC here.
221 : ! Normal DFT computes E(rho_alpha,rho_beta) (or its variant E(2*rho_alpha) for non-LSD)
222 : ! SIC can E(rho_alpha,rho_beta)-b*(E(rho_alpha,rho_beta)-E(rho_beta,rho_beta))
223 : ! or compute E(rho_alpha,rho_beta)-b*E(rho_alpha-rho_beta,0)
224 :
225 : ! my_scaling is the scaling needed of the standard E(rho_alpha,rho_beta) term
226 137431 : my_scaling = 1.0_dp
227 137635 : SELECT CASE (dft_control%sic_method_id)
228 : CASE (sic_none)
229 : ! all fine
230 : CASE (sic_mauri_spz, sic_ad)
231 : ! no idea yet what to do here in that case
232 204 : CPASSERT(.NOT. tau_r_valid)
233 : CASE (sic_mauri_us)
234 92 : my_scaling = 1.0_dp - dft_control%sic_scaling_b
235 : ! no idea yet what to do here in that case
236 92 : CPASSERT(.NOT. tau_r_valid)
237 : CASE (sic_eo)
238 : ! NOTHING TO BE DONE
239 : CASE DEFAULT
240 : ! this case has not yet been treated here
241 137431 : CPABORT("NYI")
242 : END SELECT
243 :
244 137431 : IF (dft_control%sic_scaling_b == 0.0_dp) THEN
245 : sic_scaling_b_zero = .TRUE.
246 : ELSE
247 137331 : sic_scaling_b_zero = .FALSE.
248 : END IF
249 :
250 137431 : IF (PRESENT(pw_env_external)) &
251 0 : pw_env => pw_env_external
252 137431 : CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool, auxbas_pw_pool=auxbas_pw_pool)
253 137431 : uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid, xc_pw_pool%pw_grid)
254 :
255 137431 : IF (.NOT. uf_grid) THEN
256 136413 : rho_r => rho_struct_r
257 :
258 136413 : IF (tau_r_valid) THEN
259 3644 : tau => tau_struct_r
260 : END IF
261 :
262 : ! for gradient corrected functional the density in g space might
263 : ! be useful so if we have it, we pass it in
264 136413 : IF (rho_g_valid) THEN
265 136313 : rho_g => rho_struct_g
266 : END IF
267 : ELSE
268 1018 : CPASSERT(rho_g_valid)
269 4072 : ALLOCATE (rho_r(mspin))
270 4072 : ALLOCATE (rho_g(mspin))
271 2036 : DO ispin = 1, mspin
272 1018 : CALL xc_pw_pool%create_pw(rho_g(ispin))
273 2036 : CALL pw_transfer(rho_struct_g(ispin), rho_g(ispin))
274 : END DO
275 2036 : DO ispin = 1, mspin
276 1018 : CALL xc_pw_pool%create_pw(rho_r(ispin))
277 2036 : CALL pw_transfer(rho_g(ispin), rho_r(ispin))
278 : END DO
279 1018 : IF (tau_r_valid) THEN
280 750 : ALLOCATE (tau(mspin))
281 500 : DO ispin = 1, mspin
282 250 : CALL xc_pw_pool%create_pw(tau(ispin))
283 250 : BLOCK
284 : TYPE(pw_c1d_gs_type) :: tau_g_aux, tau_g_xc
285 250 : CALL xc_pw_pool%create_pw(tau_g_xc)
286 250 : IF (tau_g_valid) THEN
287 250 : CALL pw_transfer(tau_struct_g(ispin), tau_g_xc)
288 : ELSE
289 0 : CALL auxbas_pw_pool%create_pw(tau_g_aux)
290 0 : CALL pw_transfer(tau_struct_r(ispin), tau_g_aux)
291 0 : CALL pw_transfer(tau_g_aux, tau_g_xc)
292 0 : CALL auxbas_pw_pool%give_back_pw(tau_g_aux)
293 : END IF
294 250 : CALL pw_transfer(tau_g_xc, tau(ispin))
295 500 : CALL xc_pw_pool%give_back_pw(tau_g_xc)
296 : END BLOCK
297 : END DO
298 : END IF
299 1018 : IF (ASSOCIATED(weights)) THEN
300 1004 : ALLOCATE (weights_xc)
301 1004 : CALL xc_pw_pool%create_pw(weights_xc)
302 : BLOCK
303 : TYPE(pw_c1d_gs_type) :: weights_g_aux, weights_g_xc
304 1004 : CALL auxbas_pw_pool%create_pw(weights_g_aux)
305 1004 : CALL xc_pw_pool%create_pw(weights_g_xc)
306 1004 : CALL pw_transfer(weights, weights_g_aux)
307 1004 : CALL pw_transfer(weights_g_aux, weights_g_xc)
308 1004 : CALL pw_transfer(weights_g_xc, weights_xc)
309 1004 : CALL xc_pw_pool%give_back_pw(weights_g_xc)
310 2008 : CALL auxbas_pw_pool%give_back_pw(weights_g_aux)
311 : END BLOCK
312 1004 : weights_use => weights_xc
313 : END IF
314 1018 : IF (ASSOCIATED(rho_nlcc)) THEN
315 28 : CPASSERT(ASSOCIATED(rho_nlcc_g))
316 28 : ALLOCATE (rho_nlcc_g_xc, rho_nlcc_xc)
317 28 : CALL xc_pw_pool%create_pw(rho_nlcc_g_xc)
318 28 : CALL xc_pw_pool%create_pw(rho_nlcc_xc)
319 28 : CALL pw_transfer(rho_nlcc_g, rho_nlcc_g_xc)
320 28 : CALL pw_transfer(rho_nlcc_g_xc, rho_nlcc_xc)
321 : rho_nlcc_use => rho_nlcc_xc
322 : rho_nlcc_g_use => rho_nlcc_g_xc
323 : END IF
324 : END IF
325 :
326 : ! add the nlcc densities
327 137403 : IF (ASSOCIATED(rho_nlcc_use)) THEN
328 594 : factor = 1.0_dp
329 1188 : DO ispin = 1, mspin
330 594 : CALL pw_axpy(rho_nlcc_use, rho_r(ispin), factor)
331 1188 : CALL pw_axpy(rho_nlcc_g_use, rho_g(ispin), factor)
332 : END DO
333 : END IF
334 :
335 : !
336 : ! here the rho_r, rho_g, tau is what it should be
337 : ! we get back the right my_vxc_rho and my_vxc_tau as required
338 : !
339 137431 : IF (native_skala_grid) THEN
340 : CALL skala_gpw_eval(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, exc=exc, &
341 : rho_r=rho_r, rho_g=rho_g, tau=tau, xc_section=xc_section, &
342 : weights=weights_use, pw_pool=xc_pw_pool, &
343 : particle_set=particle_set, cell=cell, &
344 : compute_virial=compute_virial, virial_xc=virial%pv_xc, &
345 234 : just_energy=my_just_energy, atom_force=native_skala_atom_force)
346 137311 : ELSE IF (my_just_energy) THEN
347 : exc = xc_exc_calc(rho_r=rho_r, tau=tau, &
348 : rho_g=rho_g, xc_section=xc_section, &
349 10458 : weights=weights_use, pw_pool=xc_pw_pool)
350 :
351 : ELSE
352 : CALL xc_vxc_pw_create(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_r, &
353 : rho_g=rho_g, tau=tau, exc=exc, &
354 : xc_section=xc_section, &
355 : weights=weights_use, pw_pool=xc_pw_pool, &
356 : compute_virial=compute_virial, &
357 126853 : virial_xc=virial%pv_xc)
358 : END IF
359 :
360 : ! remove the nlcc densities (keep stuff in original state)
361 137431 : IF (ASSOCIATED(rho_nlcc_use)) THEN
362 594 : factor = -1.0_dp
363 1188 : DO ispin = 1, mspin
364 594 : CALL pw_axpy(rho_nlcc_use, rho_r(ispin), factor)
365 1188 : CALL pw_axpy(rho_nlcc_g_use, rho_g(ispin), factor)
366 : END DO
367 : END IF
368 :
369 : ! calclulate non-local vdW functional
370 : ! only if this XC_SECTION has it
371 : ! if yes, we use the dispersion_env from ks_env
372 : ! this is dangerous, as it assumes a special connection xc_section -> qs_env
373 137431 : IF (vdW_nl) THEN
374 422 : CALL get_ks_env(ks_env=ks_env, para_env=para_env)
375 : ! no SIC functionals allowed
376 422 : CPASSERT(dft_control%sic_method_id == sic_none)
377 : !
378 422 : CALL pw_env_get(pw_env, vdw_pw_pool=vdw_pw_pool)
379 422 : IF (my_just_energy) THEN
380 : CALL calculate_dispersion_nonloc(my_vxc_rho, rho_r, rho_g, edisp, dispersion_env, &
381 6 : my_just_energy, vdw_pw_pool, xc_pw_pool, para_env)
382 : ELSE
383 : CALL calculate_dispersion_nonloc(my_vxc_rho, rho_r, rho_g, edisp, dispersion_env, &
384 416 : my_just_energy, vdw_pw_pool, xc_pw_pool, para_env, virial=virial)
385 : END IF
386 : END IF
387 :
388 : !! Apply rescaling to the potential if requested
389 137431 : IF (.NOT. my_just_energy) THEN
390 126973 : IF (do_adiabatic_rescaling) THEN
391 24 : IF (ASSOCIATED(my_vxc_rho)) THEN
392 62 : DO ispin = 1, SIZE(my_vxc_rho)
393 62 : CALL pw_scale(my_vxc_rho(ispin), my_adiabatic_rescale_factor)
394 : END DO
395 : END IF
396 : END IF
397 : END IF
398 :
399 137431 : IF (my_scaling /= 1.0_dp) THEN
400 92 : exc = exc*my_scaling
401 92 : IF (ASSOCIATED(my_vxc_rho)) THEN
402 180 : DO ispin = 1, SIZE(my_vxc_rho)
403 180 : CALL pw_scale(my_vxc_rho(ispin), my_scaling)
404 : END DO
405 : END IF
406 92 : IF (ASSOCIATED(my_vxc_tau)) THEN
407 0 : DO ispin = 1, SIZE(my_vxc_tau)
408 0 : CALL pw_scale(my_vxc_tau(ispin), my_scaling)
409 : END DO
410 : END IF
411 : END IF
412 :
413 : ! we have pw data for the xc, qs_ks requests coeff structure, here we transfer
414 : ! pw -> coeff
415 137431 : IF (ASSOCIATED(my_vxc_rho)) THEN
416 126973 : vxc_rho => my_vxc_rho
417 126973 : NULLIFY (my_vxc_rho)
418 : END IF
419 137431 : IF (ASSOCIATED(my_vxc_tau)) THEN
420 2816 : vxc_tau => my_vxc_tau
421 2816 : NULLIFY (my_vxc_tau)
422 : END IF
423 137431 : IF (uf_grid) THEN
424 2036 : DO ispin = 1, SIZE(rho_r)
425 2036 : CALL xc_pw_pool%give_back_pw(rho_r(ispin))
426 : END DO
427 1018 : DEALLOCATE (rho_r)
428 1018 : IF (ASSOCIATED(rho_g)) THEN
429 2036 : DO ispin = 1, SIZE(rho_g)
430 2036 : CALL xc_pw_pool%give_back_pw(rho_g(ispin))
431 : END DO
432 1018 : DEALLOCATE (rho_g)
433 : END IF
434 : END IF
435 :
436 : ! compute again the xc but now for Exc(m,o) and the opposite sign
437 137431 : IF (dft_control%sic_method_id == sic_mauri_spz .AND. .NOT. sic_scaling_b_zero) THEN
438 390 : ALLOCATE (rho_m_rspace(2), rho_m_gspace(2))
439 78 : CALL xc_pw_pool%create_pw(rho_m_gspace(1))
440 78 : CALL xc_pw_pool%create_pw(rho_m_rspace(1))
441 78 : CALL pw_copy(rho_struct_r(1), rho_m_rspace(1))
442 78 : CALL pw_axpy(rho_struct_r(2), rho_m_rspace(1), alpha=-1._dp)
443 78 : CALL pw_copy(rho_struct_g(1), rho_m_gspace(1))
444 78 : CALL pw_axpy(rho_struct_g(2), rho_m_gspace(1), alpha=-1._dp)
445 : ! bit sad, these will be just zero...
446 78 : CALL xc_pw_pool%create_pw(rho_m_gspace(2))
447 78 : CALL xc_pw_pool%create_pw(rho_m_rspace(2))
448 78 : CALL pw_zero(rho_m_rspace(2))
449 78 : CALL pw_zero(rho_m_gspace(2))
450 :
451 78 : IF (my_just_energy) THEN
452 : exc_m = xc_exc_calc(rho_r=rho_m_rspace, tau=tau, &
453 : rho_g=rho_m_gspace, xc_section=xc_section, &
454 24 : weights=weights_use, pw_pool=xc_pw_pool)
455 : ELSE
456 : ! virial untested
457 54 : CPASSERT(.NOT. compute_virial)
458 : CALL xc_vxc_pw_create(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_m_rspace, &
459 : rho_g=rho_m_gspace, tau=tau, exc=exc_m, &
460 : xc_section=xc_section, &
461 : weights=weights_use, pw_pool=xc_pw_pool, &
462 : compute_virial=.FALSE., &
463 54 : virial_xc=virial_xc_tmp)
464 : END IF
465 :
466 78 : exc = exc - dft_control%sic_scaling_b*exc_m
467 :
468 : ! and take care of the potential only vxc_rho is taken into account
469 78 : IF (.NOT. my_just_energy) THEN
470 54 : CALL pw_axpy(my_vxc_rho(1), vxc_rho(1), -dft_control%sic_scaling_b)
471 54 : CALL pw_axpy(my_vxc_rho(1), vxc_rho(2), dft_control%sic_scaling_b)
472 54 : CALL my_vxc_rho(1)%release()
473 54 : CALL my_vxc_rho(2)%release()
474 54 : DEALLOCATE (my_vxc_rho)
475 : END IF
476 :
477 234 : DO ispin = 1, 2
478 156 : CALL xc_pw_pool%give_back_pw(rho_m_rspace(ispin))
479 234 : CALL xc_pw_pool%give_back_pw(rho_m_gspace(ispin))
480 : END DO
481 78 : DEALLOCATE (rho_m_rspace)
482 78 : DEALLOCATE (rho_m_gspace)
483 :
484 : END IF
485 :
486 : ! now we have - sum_s N_s * Exc(rho_s/N_s,0)
487 137431 : IF (dft_control%sic_method_id == sic_ad .AND. .NOT. sic_scaling_b_zero) THEN
488 :
489 : ! find out how many elecs we have
490 26 : CALL get_ks_env(ks_env, nelectron_spin=nelec_spin)
491 :
492 130 : ALLOCATE (rho_m_rspace(2), rho_m_gspace(2))
493 78 : DO ispin = 1, 2
494 52 : CALL xc_pw_pool%create_pw(rho_m_gspace(ispin))
495 78 : CALL xc_pw_pool%create_pw(rho_m_rspace(ispin))
496 : END DO
497 :
498 78 : DO ispin = 1, 2
499 52 : IF (nelec_spin(ispin) > 0.0_dp) THEN
500 52 : nelec_s_inv = 1.0_dp/nelec_spin(ispin)
501 : ELSE
502 : ! does it matter if there are no electrons with this spin (H) ?
503 0 : nelec_s_inv = 0.0_dp
504 : END IF
505 52 : CALL pw_copy(rho_struct_r(ispin), rho_m_rspace(1))
506 52 : CALL pw_copy(rho_struct_g(ispin), rho_m_gspace(1))
507 52 : CALL pw_scale(rho_m_rspace(1), nelec_s_inv)
508 52 : CALL pw_scale(rho_m_gspace(1), nelec_s_inv)
509 52 : CALL pw_zero(rho_m_rspace(2))
510 52 : CALL pw_zero(rho_m_gspace(2))
511 :
512 52 : IF (my_just_energy) THEN
513 : exc_m = xc_exc_calc(rho_r=rho_m_rspace, tau=tau, &
514 : rho_g=rho_m_gspace, xc_section=xc_section, &
515 12 : weights=weights_use, pw_pool=xc_pw_pool)
516 : ELSE
517 : ! virial untested
518 40 : CPASSERT(.NOT. compute_virial)
519 : CALL xc_vxc_pw_create(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_m_rspace, &
520 : rho_g=rho_m_gspace, tau=tau, exc=exc_m, &
521 : xc_section=xc_section, &
522 : weights=weights_use, pw_pool=xc_pw_pool, &
523 : compute_virial=.FALSE., &
524 40 : virial_xc=virial_xc_tmp)
525 : END IF
526 :
527 52 : exc = exc - dft_control%sic_scaling_b*nelec_spin(ispin)*exc_m
528 :
529 : ! and take care of the potential only vxc_rho is taken into account
530 78 : IF (.NOT. my_just_energy) THEN
531 40 : CALL pw_axpy(my_vxc_rho(1), vxc_rho(ispin), -dft_control%sic_scaling_b)
532 40 : CALL my_vxc_rho(1)%release()
533 40 : CALL my_vxc_rho(2)%release()
534 40 : DEALLOCATE (my_vxc_rho)
535 : END IF
536 : END DO
537 :
538 78 : DO ispin = 1, 2
539 52 : CALL xc_pw_pool%give_back_pw(rho_m_rspace(ispin))
540 78 : CALL xc_pw_pool%give_back_pw(rho_m_gspace(ispin))
541 : END DO
542 26 : DEALLOCATE (rho_m_rspace)
543 26 : DEALLOCATE (rho_m_gspace)
544 :
545 : END IF
546 :
547 : ! compute again the xc but now for Exc(n_down,n_down)
548 137431 : IF (dft_control%sic_method_id == sic_mauri_us .AND. .NOT. sic_scaling_b_zero) THEN
549 276 : ALLOCATE (rho_r(2))
550 92 : rho_r(1) = rho_struct_r(2)
551 92 : rho_r(2) = rho_struct_r(2)
552 92 : IF (rho_g_valid) THEN
553 276 : ALLOCATE (rho_g(2))
554 92 : rho_g(1) = rho_struct_g(2)
555 92 : rho_g(2) = rho_struct_g(2)
556 : END IF
557 :
558 92 : IF (my_just_energy) THEN
559 : exc_m = xc_exc_calc(rho_r=rho_r, tau=tau, &
560 : rho_g=rho_g, xc_section=xc_section, &
561 32 : weights=weights_use, pw_pool=xc_pw_pool)
562 : ELSE
563 : ! virial untested
564 60 : CPASSERT(.NOT. compute_virial)
565 : CALL xc_vxc_pw_create(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_r, &
566 : rho_g=rho_g, tau=tau, exc=exc_m, &
567 : xc_section=xc_section, &
568 : weights=weights_use, pw_pool=xc_pw_pool, &
569 : compute_virial=.FALSE., &
570 60 : virial_xc=virial_xc_tmp)
571 : END IF
572 :
573 92 : exc = exc + dft_control%sic_scaling_b*exc_m
574 :
575 : ! and take care of the potential
576 92 : IF (.NOT. my_just_energy) THEN
577 : ! both go to minority spin
578 60 : CALL pw_axpy(my_vxc_rho(1), vxc_rho(2), 2.0_dp*dft_control%sic_scaling_b)
579 60 : CALL my_vxc_rho(1)%release()
580 60 : CALL my_vxc_rho(2)%release()
581 60 : DEALLOCATE (my_vxc_rho)
582 : END IF
583 92 : DEALLOCATE (rho_r, rho_g)
584 :
585 : END IF
586 :
587 : !
588 : ! cleanups
589 : !
590 137431 : IF (uf_grid .AND. (ASSOCIATED(vxc_rho) .OR. ASSOCIATED(vxc_tau))) THEN
591 : BLOCK
592 : TYPE(pw_r3d_rs_type) :: tmp_pw
593 : TYPE(pw_c1d_gs_type) :: tmp_g, tmp_g2
594 1018 : CALL xc_pw_pool%create_pw(tmp_g)
595 1018 : CALL auxbas_pw_pool%create_pw(tmp_g2)
596 1018 : IF (ASSOCIATED(vxc_rho)) THEN
597 2036 : DO ispin = 1, SIZE(vxc_rho)
598 1018 : CALL auxbas_pw_pool%create_pw(tmp_pw)
599 1018 : CALL pw_transfer(vxc_rho(ispin), tmp_g)
600 1018 : CALL pw_transfer(tmp_g, tmp_g2)
601 1018 : CALL pw_transfer(tmp_g2, tmp_pw)
602 1018 : CALL xc_pw_pool%give_back_pw(vxc_rho(ispin))
603 2036 : vxc_rho(ispin) = tmp_pw
604 : END DO
605 : END IF
606 1018 : IF (ASSOCIATED(vxc_tau)) THEN
607 500 : DO ispin = 1, SIZE(vxc_tau)
608 250 : CALL auxbas_pw_pool%create_pw(tmp_pw)
609 250 : CALL pw_transfer(vxc_tau(ispin), tmp_g)
610 250 : CALL pw_transfer(tmp_g, tmp_g2)
611 250 : CALL pw_transfer(tmp_g2, tmp_pw)
612 250 : CALL xc_pw_pool%give_back_pw(vxc_tau(ispin))
613 500 : vxc_tau(ispin) = tmp_pw
614 : END DO
615 : END IF
616 1018 : CALL auxbas_pw_pool%give_back_pw(tmp_g2)
617 2036 : CALL xc_pw_pool%give_back_pw(tmp_g)
618 : END BLOCK
619 : END IF
620 137431 : IF (ASSOCIATED(tau) .AND. uf_grid) THEN
621 500 : DO ispin = 1, SIZE(tau)
622 500 : CALL xc_pw_pool%give_back_pw(tau(ispin))
623 : END DO
624 250 : DEALLOCATE (tau)
625 : END IF
626 137431 : IF (ASSOCIATED(weights_xc)) THEN
627 1004 : CALL xc_pw_pool%give_back_pw(weights_xc)
628 1004 : DEALLOCATE (weights_xc)
629 : END IF
630 137431 : IF (ASSOCIATED(rho_nlcc_xc)) THEN
631 28 : CALL xc_pw_pool%give_back_pw(rho_nlcc_xc)
632 28 : DEALLOCATE (rho_nlcc_xc)
633 : END IF
634 137431 : IF (ASSOCIATED(rho_nlcc_g_xc)) THEN
635 28 : CALL xc_pw_pool%give_back_pw(rho_nlcc_g_xc)
636 28 : DEALLOCATE (rho_nlcc_g_xc)
637 : END IF
638 :
639 : END IF
640 :
641 152073 : CALL timestop(handle)
642 :
643 152073 : END SUBROUTINE qs_vxc_create
644 :
645 : ! **************************************************************************************************
646 : !> \brief calculates the XC density: E_xc(r) - V_xc(r)*rho(r) or E_xc(r)/rho(r)
647 : !> \param ks_env to get all the needed things
648 : !> \param rho_struct density
649 : !> \param xc_section ...
650 : !> \param dispersion_env ...
651 : !> \param xc_ener will contain the xc energy density E_xc(r) - V_xc(r)*rho(r)
652 : !> \param xc_den will contain the xc energy density E_xc(r)/rho(r)
653 : !> \param exc will contain the xc energy density E_xc(r)
654 : !> \param vxc ...
655 : !> \param vtau ...
656 : !> \author JGH
657 : ! **************************************************************************************************
658 500 : SUBROUTINE qs_xc_density(ks_env, rho_struct, xc_section, dispersion_env, &
659 100 : xc_ener, xc_den, exc, vxc, vtau)
660 :
661 : TYPE(qs_ks_env_type), POINTER :: ks_env
662 : TYPE(qs_rho_type), POINTER :: rho_struct
663 : TYPE(section_vals_type), POINTER :: xc_section
664 : TYPE(qs_dispersion_type), OPTIONAL, POINTER :: dispersion_env
665 : TYPE(pw_r3d_rs_type), INTENT(INOUT), OPTIONAL :: xc_ener, xc_den
666 : TYPE(pw_r3d_rs_type), OPTIONAL :: exc
667 : TYPE(pw_r3d_rs_type), DIMENSION(:), OPTIONAL :: vxc, vtau
668 :
669 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_xc_density'
670 :
671 : INTEGER :: handle, ispin, mspin, myfun, nspins, vdw
672 : LOGICAL :: rho_g_valid, tau_g_valid, tau_r_valid, &
673 : uf_grid, vdW_nl
674 : REAL(KIND=dp) :: edisp, excint, factor, rho_cutoff
675 : REAL(KIND=dp), DIMENSION(3, 3) :: vdum
676 : TYPE(cell_type), POINTER :: cell
677 : TYPE(dft_control_type), POINTER :: dft_control
678 : TYPE(mp_para_env_type), POINTER :: para_env
679 100 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g, rho_struct_g, tau_g, tau_struct_g
680 : TYPE(pw_c1d_gs_type), POINTER :: rho_nlcc_g, rho_nlcc_g_use, rho_nlcc_g_xc
681 : TYPE(pw_env_type), POINTER :: pw_env
682 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool, vdw_pw_pool, xc_pw_pool
683 : TYPE(pw_r3d_rs_type) :: exc_r
684 100 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, rho_struct_r, tau_r, &
685 100 : tau_struct_r, vxc_rho, vxc_tau
686 : TYPE(pw_r3d_rs_type), POINTER :: rho_nlcc, rho_nlcc_use, rho_nlcc_xc, &
687 : weights, weights_use, weights_xc
688 :
689 100 : CALL timeset(routineN, handle)
690 :
691 100 : NULLIFY (dft_control, pw_env, auxbas_pw_pool, xc_pw_pool, vdw_pw_pool, cell, &
692 100 : rho_g, rho_struct_g, tau_g, tau_struct_g, rho_nlcc, rho_nlcc_g, &
693 100 : rho_nlcc_g_use, rho_nlcc_g_xc, rho_nlcc_use, rho_nlcc_xc, rho_r, &
694 100 : rho_struct_r, tau_r, tau_struct_r, vxc_rho, vxc_tau, weights, &
695 100 : weights_use, weights_xc)
696 :
697 : CALL get_ks_env(ks_env, &
698 : dft_control=dft_control, &
699 : pw_env=pw_env, &
700 : cell=cell, &
701 : xcint_weights=weights, &
702 : rho_nlcc=rho_nlcc, &
703 100 : rho_nlcc_g=rho_nlcc_g)
704 :
705 : CALL qs_rho_get(rho_struct, &
706 : tau_r_valid=tau_r_valid, &
707 : tau_g_valid=tau_g_valid, &
708 : rho_g_valid=rho_g_valid, &
709 : rho_r=rho_struct_r, &
710 : rho_g=rho_struct_g, &
711 : tau_r=tau_struct_r, &
712 100 : tau_g=tau_struct_g)
713 100 : nspins = dft_control%nspins
714 100 : mspin = SIZE(rho_struct_r)
715 100 : rho_r => rho_struct_r
716 100 : rho_g => rho_struct_g
717 100 : tau_r => tau_struct_r
718 100 : tau_g => tau_struct_g
719 100 : rho_nlcc_use => rho_nlcc
720 100 : rho_nlcc_g_use => rho_nlcc_g
721 100 : weights_use => weights
722 :
723 100 : CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", i_val=myfun)
724 100 : CALL section_vals_val_get(xc_section, "VDW_POTENTIAL%POTENTIAL_TYPE", i_val=vdw)
725 100 : vdW_nl = (vdw == xc_vdw_fun_nonloc)
726 100 : IF (PRESENT(xc_ener)) THEN
727 34 : IF (tau_r_valid) THEN
728 0 : CALL cp_warn(__LOCATION__, "Tau contribution will not be correctly handled")
729 : END IF
730 : END IF
731 100 : IF (vdW_nl) THEN
732 0 : CALL cp_warn(__LOCATION__, "vdW functional contribution will be ignored")
733 : END IF
734 :
735 100 : CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool, auxbas_pw_pool=auxbas_pw_pool)
736 100 : uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid, xc_pw_pool%pw_grid)
737 :
738 100 : IF (PRESENT(xc_ener)) THEN
739 34 : CALL pw_zero(xc_ener)
740 : END IF
741 100 : IF (PRESENT(xc_den)) THEN
742 66 : CALL pw_zero(xc_den)
743 : END IF
744 100 : IF (PRESENT(exc)) THEN
745 0 : CALL pw_zero(exc)
746 : END IF
747 100 : IF (PRESENT(vxc)) THEN
748 138 : DO ispin = 1, nspins
749 138 : CALL pw_zero(vxc(ispin))
750 : END DO
751 : END IF
752 100 : IF (PRESENT(vtau)) THEN
753 40 : DO ispin = 1, nspins
754 40 : CALL pw_zero(vtau(ispin))
755 : END DO
756 : END IF
757 :
758 100 : IF (myfun /= xc_none) THEN
759 :
760 98 : CPASSERT(ASSOCIATED(rho_struct))
761 98 : CPASSERT(dft_control%sic_method_id == sic_none)
762 :
763 98 : IF (uf_grid) THEN
764 2 : NULLIFY (rho_r, rho_g, tau_r, tau_g)
765 2 : IF (rho_g_valid) THEN
766 2 : CALL create_density_on_pool(xc_pw_pool, rho_struct_g, rho_r, rho_g)
767 0 : ELSEIF (ASSOCIATED(rho_struct_r)) THEN
768 0 : CALL create_density_on_pool_from_r(auxbas_pw_pool, xc_pw_pool, rho_struct_r, rho_r, rho_g)
769 : ELSE
770 0 : CPABORT("Fine Grid in qs_xc_density requires rho_r or rho_g")
771 : END IF
772 2 : IF (tau_r_valid) THEN
773 0 : IF (tau_g_valid) THEN
774 0 : CALL create_density_on_pool(xc_pw_pool, tau_struct_g, tau_r, tau_g)
775 0 : ELSEIF (ASSOCIATED(tau_struct_r)) THEN
776 0 : CALL create_density_on_pool_from_r(auxbas_pw_pool, xc_pw_pool, tau_struct_r, tau_r, tau_g)
777 : ELSE
778 0 : CPABORT("Fine Grid in qs_xc_density requires tau_r or tau_g")
779 : END IF
780 : END IF
781 2 : IF (ASSOCIATED(weights)) THEN
782 2 : ALLOCATE (weights_xc)
783 2 : CALL xc_pw_pool%create_pw(weights_xc)
784 2 : CALL transfer_rspace_between_pools(auxbas_pw_pool, xc_pw_pool, weights, weights_xc)
785 2 : weights_use => weights_xc
786 : END IF
787 2 : IF (ASSOCIATED(rho_nlcc)) THEN
788 0 : CPASSERT(ASSOCIATED(rho_nlcc_g))
789 0 : ALLOCATE (rho_nlcc_g_xc, rho_nlcc_xc)
790 0 : CALL xc_pw_pool%create_pw(rho_nlcc_g_xc)
791 0 : CALL xc_pw_pool%create_pw(rho_nlcc_xc)
792 0 : CALL pw_transfer(rho_nlcc_g, rho_nlcc_g_xc)
793 0 : CALL pw_transfer(rho_nlcc_g_xc, rho_nlcc_xc)
794 : rho_nlcc_use => rho_nlcc_xc
795 : rho_nlcc_g_use => rho_nlcc_g_xc
796 : END IF
797 : END IF
798 :
799 : ! add the nlcc densities
800 98 : IF (ASSOCIATED(rho_nlcc_use)) THEN
801 0 : factor = 1.0_dp
802 0 : DO ispin = 1, mspin
803 0 : CALL pw_axpy(rho_nlcc_use, rho_r(ispin), factor)
804 0 : CALL pw_axpy(rho_nlcc_g_use, rho_g(ispin), factor)
805 : END DO
806 : END IF
807 98 : NULLIFY (vxc_rho, vxc_tau)
808 : CALL xc_vxc_pw_create(vxc_rho=vxc_rho, vxc_tau=vxc_tau, rho_r=rho_r, &
809 : rho_g=rho_g, tau=tau_r, exc=excint, &
810 : xc_section=xc_section, &
811 : weights=weights_use, pw_pool=xc_pw_pool, &
812 : compute_virial=.FALSE., &
813 : virial_xc=vdum, &
814 98 : exc_r=exc_r)
815 : ! calclulate non-local vdW functional
816 : ! only if this XC_SECTION has it
817 : ! if yes, we use the dispersion_env from ks_env
818 : ! this is dangerous, as it assumes a special connection xc_section -> qs_env
819 98 : IF (vdW_nl) THEN
820 0 : CALL get_ks_env(ks_env=ks_env, para_env=para_env)
821 : ! no SIC functionals allowed
822 0 : CPASSERT(dft_control%sic_method_id == sic_none)
823 : !
824 0 : CALL pw_env_get(pw_env, vdw_pw_pool=vdw_pw_pool)
825 : CALL calculate_dispersion_nonloc(vxc_rho, rho_r, rho_g, edisp, dispersion_env, &
826 0 : .FALSE., vdw_pw_pool, xc_pw_pool, para_env)
827 : END IF
828 :
829 : ! remove the nlcc densities (keep stuff in original state)
830 98 : IF (ASSOCIATED(rho_nlcc_use)) THEN
831 0 : factor = -1.0_dp
832 0 : DO ispin = 1, mspin
833 0 : CALL pw_axpy(rho_nlcc_use, rho_r(ispin), factor)
834 0 : CALL pw_axpy(rho_nlcc_g_use, rho_g(ispin), factor)
835 : END DO
836 : END IF
837 : !
838 98 : IF (PRESENT(xc_den)) THEN
839 64 : rho_cutoff = 1.E-14_dp
840 64 : IF (uf_grid) THEN
841 : BLOCK
842 : TYPE(pw_r3d_rs_type) :: tmp_pw
843 0 : CALL xc_pw_pool%create_pw(tmp_pw)
844 0 : CALL pw_copy(exc_r, tmp_pw)
845 0 : CALL calc_xc_density(tmp_pw, rho_r, rho_cutoff)
846 0 : CALL transfer_rspace_between_pools(xc_pw_pool, auxbas_pw_pool, tmp_pw, xc_den)
847 0 : CALL xc_pw_pool%give_back_pw(tmp_pw)
848 : END BLOCK
849 : ELSE
850 64 : CALL pw_copy(exc_r, xc_den)
851 64 : CALL calc_xc_density(xc_den, rho_r, rho_cutoff)
852 : END IF
853 : END IF
854 98 : IF (PRESENT(xc_ener)) THEN
855 34 : IF (uf_grid) THEN
856 : BLOCK
857 : TYPE(pw_r3d_rs_type) :: tmp_pw
858 2 : CALL xc_pw_pool%create_pw(tmp_pw)
859 2 : CALL pw_copy(exc_r, tmp_pw)
860 4 : DO ispin = 1, nspins
861 4 : CALL pw_multiply(tmp_pw, vxc_rho(ispin), rho_r(ispin), alpha=-1.0_dp)
862 : END DO
863 2 : CALL transfer_rspace_between_pools(xc_pw_pool, auxbas_pw_pool, tmp_pw, xc_ener)
864 2 : CALL xc_pw_pool%give_back_pw(tmp_pw)
865 : END BLOCK
866 : ELSE
867 32 : CALL pw_copy(exc_r, xc_ener)
868 64 : DO ispin = 1, nspins
869 64 : CALL pw_multiply(xc_ener, vxc_rho(ispin), rho_r(ispin), alpha=-1.0_dp)
870 : END DO
871 : END IF
872 : END IF
873 98 : IF (PRESENT(exc)) THEN
874 0 : IF (uf_grid) THEN
875 0 : CALL transfer_rspace_between_pools(xc_pw_pool, auxbas_pw_pool, exc_r, exc)
876 : ELSE
877 0 : CALL pw_copy(exc_r, exc)
878 : END IF
879 : END IF
880 98 : IF (PRESENT(vxc)) THEN
881 134 : DO ispin = 1, nspins
882 134 : IF (uf_grid) THEN
883 0 : CALL transfer_rspace_between_pools(xc_pw_pool, auxbas_pw_pool, vxc_rho(ispin), vxc(ispin))
884 : ELSE
885 70 : CALL pw_copy(vxc_rho(ispin), vxc(ispin))
886 : END IF
887 : END DO
888 : END IF
889 98 : IF (PRESENT(vtau) .AND. ASSOCIATED(vxc_tau)) THEN
890 40 : DO ispin = 1, nspins
891 40 : IF (uf_grid) THEN
892 0 : CALL transfer_rspace_between_pools(xc_pw_pool, auxbas_pw_pool, vxc_tau(ispin), vtau(ispin))
893 : ELSE
894 20 : CALL pw_copy(vxc_tau(ispin), vtau(ispin))
895 : END IF
896 : END DO
897 : END IF
898 : ! remove arrays
899 98 : IF (ASSOCIATED(vxc_rho)) THEN
900 202 : DO ispin = 1, nspins
901 202 : CALL vxc_rho(ispin)%release()
902 : END DO
903 98 : DEALLOCATE (vxc_rho)
904 : END IF
905 98 : IF (ASSOCIATED(vxc_tau)) THEN
906 40 : DO ispin = 1, nspins
907 40 : CALL vxc_tau(ispin)%release()
908 : END DO
909 20 : DEALLOCATE (vxc_tau)
910 : END IF
911 98 : CALL exc_r%release()
912 98 : IF (uf_grid) THEN
913 2 : CALL give_back_density_on_pool(xc_pw_pool, rho_r, rho_g)
914 2 : IF (ASSOCIATED(tau_r)) CALL give_back_density_on_pool(xc_pw_pool, tau_r, tau_g)
915 2 : IF (ASSOCIATED(weights_xc)) THEN
916 2 : CALL xc_pw_pool%give_back_pw(weights_xc)
917 2 : DEALLOCATE (weights_xc)
918 : END IF
919 2 : IF (ASSOCIATED(rho_nlcc_xc)) THEN
920 0 : CALL xc_pw_pool%give_back_pw(rho_nlcc_xc)
921 0 : DEALLOCATE (rho_nlcc_xc)
922 : END IF
923 2 : IF (ASSOCIATED(rho_nlcc_g_xc)) THEN
924 0 : CALL xc_pw_pool%give_back_pw(rho_nlcc_g_xc)
925 0 : DEALLOCATE (rho_nlcc_g_xc)
926 : END IF
927 : END IF
928 : !
929 : END IF
930 :
931 100 : CALL timestop(handle)
932 :
933 100 : END SUBROUTINE qs_xc_density
934 :
935 : ! **************************************************************************************************
936 : !> \brief transfers an r-space PW between two pools and writes into an existing target PW
937 : !> \param source_pw_pool ...
938 : !> \param target_pw_pool ...
939 : !> \param source ...
940 : !> \param TARGET ...
941 : ! **************************************************************************************************
942 4 : SUBROUTINE transfer_rspace_between_pools(source_pw_pool, target_pw_pool, source, TARGET)
943 : TYPE(pw_pool_type), POINTER :: source_pw_pool, target_pw_pool
944 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: source, TARGET
945 :
946 : TYPE(pw_c1d_gs_type) :: source_g, target_g
947 :
948 0 : CPASSERT(ASSOCIATED(source_pw_pool))
949 4 : CPASSERT(ASSOCIATED(target_pw_pool))
950 :
951 4 : IF (pw_grid_compare(source_pw_pool%pw_grid, target_pw_pool%pw_grid)) THEN
952 0 : CALL pw_copy(source, TARGET)
953 : ELSE
954 4 : CALL source_pw_pool%create_pw(source_g)
955 4 : CALL target_pw_pool%create_pw(target_g)
956 4 : CALL pw_transfer(source, source_g)
957 4 : CALL pw_transfer(source_g, target_g)
958 4 : CALL pw_transfer(target_g, TARGET)
959 4 : CALL target_pw_pool%give_back_pw(target_g)
960 4 : CALL source_pw_pool%give_back_pw(source_g)
961 : END IF
962 :
963 4 : END SUBROUTINE transfer_rspace_between_pools
964 :
965 : ! **************************************************************************************************
966 : !> \brief transfers a g-space density to a given PW pool and creates its r-space representation
967 : !> \param pw_pool ...
968 : !> \param rho_g_in ...
969 : !> \param rho_r_out ...
970 : !> \param rho_g_out ...
971 : ! **************************************************************************************************
972 2 : SUBROUTINE create_density_on_pool(pw_pool, rho_g_in, rho_r_out, rho_g_out)
973 : TYPE(pw_pool_type), POINTER :: pw_pool
974 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_in
975 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_out
976 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_out
977 :
978 : INTEGER :: ispin, nspins
979 :
980 2 : CPASSERT(ASSOCIATED(pw_pool))
981 2 : CPASSERT(ASSOCIATED(rho_g_in))
982 :
983 2 : nspins = SIZE(rho_g_in)
984 14 : ALLOCATE (rho_r_out(nspins), rho_g_out(nspins))
985 4 : DO ispin = 1, nspins
986 2 : CALL pw_pool%create_pw(rho_g_out(ispin))
987 2 : CALL pw_pool%create_pw(rho_r_out(ispin))
988 2 : CALL pw_transfer(rho_g_in(ispin), rho_g_out(ispin))
989 4 : CALL pw_transfer(rho_g_out(ispin), rho_r_out(ispin))
990 : END DO
991 :
992 2 : END SUBROUTINE create_density_on_pool
993 :
994 : ! **************************************************************************************************
995 : !> \brief transfers an r-space density to a given PW pool and creates its g-space representation
996 : !> \param source_pw_pool ...
997 : !> \param target_pw_pool ...
998 : !> \param rho_r_in ...
999 : !> \param rho_r_out ...
1000 : !> \param rho_g_out ...
1001 : ! **************************************************************************************************
1002 0 : SUBROUTINE create_density_on_pool_from_r(source_pw_pool, target_pw_pool, rho_r_in, rho_r_out, rho_g_out)
1003 : TYPE(pw_pool_type), POINTER :: source_pw_pool, target_pw_pool
1004 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_in, rho_r_out
1005 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_out
1006 :
1007 : INTEGER :: ispin, nspins
1008 : TYPE(pw_c1d_gs_type) :: rho_g_in
1009 :
1010 0 : CPASSERT(ASSOCIATED(source_pw_pool))
1011 0 : CPASSERT(ASSOCIATED(target_pw_pool))
1012 0 : CPASSERT(ASSOCIATED(rho_r_in))
1013 :
1014 0 : nspins = SIZE(rho_r_in)
1015 0 : ALLOCATE (rho_r_out(nspins), rho_g_out(nspins))
1016 0 : DO ispin = 1, nspins
1017 0 : CALL source_pw_pool%create_pw(rho_g_in)
1018 0 : CALL target_pw_pool%create_pw(rho_g_out(ispin))
1019 0 : CALL target_pw_pool%create_pw(rho_r_out(ispin))
1020 0 : CALL pw_transfer(rho_r_in(ispin), rho_g_in)
1021 0 : CALL pw_transfer(rho_g_in, rho_g_out(ispin))
1022 0 : CALL pw_transfer(rho_g_out(ispin), rho_r_out(ispin))
1023 0 : CALL source_pw_pool%give_back_pw(rho_g_in)
1024 : END DO
1025 :
1026 0 : END SUBROUTINE create_density_on_pool_from_r
1027 :
1028 : ! **************************************************************************************************
1029 : !> \brief returns temporary density arrays to the given PW pool
1030 : !> \param pw_pool ...
1031 : !> \param rho_r ...
1032 : !> \param rho_g ...
1033 : ! **************************************************************************************************
1034 2 : SUBROUTINE give_back_density_on_pool(pw_pool, rho_r, rho_g)
1035 : TYPE(pw_pool_type), POINTER :: pw_pool
1036 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
1037 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
1038 :
1039 : INTEGER :: ispin
1040 :
1041 2 : CPASSERT(ASSOCIATED(pw_pool))
1042 :
1043 2 : IF (ASSOCIATED(rho_r)) THEN
1044 4 : DO ispin = 1, SIZE(rho_r)
1045 4 : CALL pw_pool%give_back_pw(rho_r(ispin))
1046 : END DO
1047 2 : DEALLOCATE (rho_r)
1048 : END IF
1049 2 : IF (ASSOCIATED(rho_g)) THEN
1050 4 : DO ispin = 1, SIZE(rho_g)
1051 4 : CALL pw_pool%give_back_pw(rho_g(ispin))
1052 : END DO
1053 2 : DEALLOCATE (rho_g)
1054 : END IF
1055 :
1056 2 : END SUBROUTINE give_back_density_on_pool
1057 :
1058 : END MODULE qs_vxc
|
