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