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