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