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 570636 : 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_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 142659 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g, rho_m_gspace, rho_struct_g
124 : TYPE(pw_c1d_gs_type), POINTER :: rho_nlcc_g, tmp_g, tmp_g2
125 : TYPE(pw_env_type), POINTER :: pw_env
126 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool, vdw_pw_pool, xc_pw_pool
127 142659 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: my_vxc_rho, my_vxc_tau, rho_m_rspace, &
128 142659 : rho_r, rho_struct_r, tau, tau_struct_r
129 : TYPE(pw_r3d_rs_type), POINTER :: rho_nlcc, tmp_pw, weights
130 : TYPE(virial_type), POINTER :: virial
131 :
132 142659 : CALL timeset(routineN, handle)
133 :
134 142659 : CPASSERT(.NOT. ASSOCIATED(vxc_rho))
135 142659 : CPASSERT(.NOT. ASSOCIATED(vxc_tau))
136 142659 : NULLIFY (dft_control, pw_env, auxbas_pw_pool, xc_pw_pool, vdw_pw_pool, cell, my_vxc_rho, &
137 142659 : tmp_pw, tmp_g, tmp_g2, my_vxc_tau, rho_g, rho_r, tau, rho_m_rspace, &
138 142659 : rho_m_gspace, rho_nlcc, rho_nlcc_g, rho_struct_r, rho_struct_g, tau_struct_r)
139 :
140 142659 : exc = 0.0_dp
141 142659 : my_just_energy = .FALSE.
142 142659 : IF (PRESENT(just_energy)) my_just_energy = just_energy
143 142659 : my_adiabatic_rescale_factor = 1.0_dp
144 142659 : do_adiabatic_rescaling = .FALSE.
145 142659 : IF (PRESENT(adiabatic_rescale_factor)) THEN
146 44 : my_adiabatic_rescale_factor = adiabatic_rescale_factor
147 44 : do_adiabatic_rescaling = .TRUE.
148 : END IF
149 :
150 : CALL get_ks_env(ks_env, &
151 : dft_control=dft_control, &
152 : pw_env=pw_env, &
153 : cell=cell, &
154 : xcint_weights=weights, &
155 : virial=virial, &
156 : rho_nlcc=rho_nlcc, &
157 142659 : rho_nlcc_g=rho_nlcc_g)
158 :
159 : CALL qs_rho_get(rho_struct, &
160 : tau_r_valid=tau_r_valid, &
161 : rho_g_valid=rho_g_valid, &
162 : rho_r=rho_struct_r, &
163 : rho_g=rho_struct_g, &
164 142659 : tau_r=tau_struct_r)
165 :
166 142659 : compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
167 142659 : IF (compute_virial) THEN
168 35958 : virial%pv_xc = 0.0_dp
169 : END IF
170 :
171 : CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", &
172 142659 : i_val=myfun)
173 : CALL section_vals_val_get(xc_section, "VDW_POTENTIAL%POTENTIAL_TYPE", &
174 142659 : i_val=vdw)
175 :
176 142659 : vdW_nl = (vdw == xc_vdw_fun_nonloc)
177 : ! this combination has not been investigated
178 142659 : CPASSERT(.NOT. (do_adiabatic_rescaling .AND. vdW_nl))
179 : ! are the necessary inputs available
180 142659 : IF (.NOT. (PRESENT(dispersion_env) .AND. PRESENT(edisp))) THEN
181 : vdW_nl = .FALSE.
182 : END IF
183 142659 : IF (PRESENT(edisp)) edisp = 0.0_dp
184 :
185 142659 : IF (myfun /= xc_none .OR. vdW_nl) THEN
186 :
187 : ! test if the real space density is available
188 128105 : CPASSERT(ASSOCIATED(rho_struct))
189 128105 : IF (dft_control%nspins /= 1 .AND. dft_control%nspins /= 2) &
190 0 : CPABORT("nspins must be 1 or 2")
191 128105 : mspin = SIZE(rho_struct_r)
192 128105 : IF (dft_control%nspins == 2 .AND. mspin == 1) &
193 0 : CPABORT("Spin count mismatch")
194 :
195 : ! there are some options related to SIC here.
196 : ! Normal DFT computes E(rho_alpha,rho_beta) (or its variant E(2*rho_alpha) for non-LSD)
197 : ! SIC can E(rho_alpha,rho_beta)-b*(E(rho_alpha,rho_beta)-E(rho_beta,rho_beta))
198 : ! or compute E(rho_alpha,rho_beta)-b*E(rho_alpha-rho_beta,0)
199 :
200 : ! my_scaling is the scaling needed of the standard E(rho_alpha,rho_beta) term
201 128105 : my_scaling = 1.0_dp
202 128309 : SELECT CASE (dft_control%sic_method_id)
203 : CASE (sic_none)
204 : ! all fine
205 : CASE (sic_mauri_spz, sic_ad)
206 : ! no idea yet what to do here in that case
207 204 : CPASSERT(.NOT. tau_r_valid)
208 : CASE (sic_mauri_us)
209 92 : my_scaling = 1.0_dp - dft_control%sic_scaling_b
210 : ! no idea yet what to do here in that case
211 92 : CPASSERT(.NOT. tau_r_valid)
212 : CASE (sic_eo)
213 : ! NOTHING TO BE DONE
214 : CASE DEFAULT
215 : ! this case has not yet been treated here
216 128105 : CPABORT("NYI")
217 : END SELECT
218 :
219 128105 : IF (dft_control%sic_scaling_b == 0.0_dp) THEN
220 : sic_scaling_b_zero = .TRUE.
221 : ELSE
222 128005 : sic_scaling_b_zero = .FALSE.
223 : END IF
224 :
225 128105 : IF (PRESENT(pw_env_external)) &
226 0 : pw_env => pw_env_external
227 128105 : CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool, auxbas_pw_pool=auxbas_pw_pool)
228 128105 : uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid, xc_pw_pool%pw_grid)
229 :
230 128105 : IF (.NOT. uf_grid) THEN
231 128091 : rho_r => rho_struct_r
232 :
233 128091 : IF (tau_r_valid) THEN
234 3220 : tau => tau_struct_r
235 : END IF
236 :
237 : ! for gradient corrected functional the density in g space might
238 : ! be useful so if we have it, we pass it in
239 128091 : IF (rho_g_valid) THEN
240 128071 : rho_g => rho_struct_g
241 : END IF
242 : ELSE
243 14 : CPASSERT(rho_g_valid)
244 56 : ALLOCATE (rho_r(mspin))
245 56 : ALLOCATE (rho_g(mspin))
246 28 : DO ispin = 1, mspin
247 14 : CALL xc_pw_pool%create_pw(rho_g(ispin))
248 28 : CALL pw_transfer(rho_struct_g(ispin), rho_g(ispin))
249 : END DO
250 28 : DO ispin = 1, mspin
251 14 : CALL xc_pw_pool%create_pw(rho_r(ispin))
252 28 : CALL pw_transfer(rho_g(ispin), rho_r(ispin))
253 : END DO
254 14 : IF (tau_r_valid) THEN
255 : ! tau with finer grids is not implemented (at least not correctly), which this asserts
256 0 : CPABORT("Tau (MGGA) with finer grids not implemented")
257 : END IF
258 14 : IF (ASSOCIATED(weights)) THEN
259 0 : CPABORT("Accurate integration with finer grids not implemented")
260 : END IF
261 : END IF
262 :
263 : ! add the nlcc densities
264 128105 : IF (ASSOCIATED(rho_nlcc)) THEN
265 566 : factor = 1.0_dp
266 1132 : DO ispin = 1, mspin
267 566 : CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
268 1132 : CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
269 : END DO
270 : END IF
271 :
272 : !
273 : ! here the rho_r, rho_g, tau is what it should be
274 : ! we get back the right my_vxc_rho and my_vxc_tau as required
275 : !
276 128105 : IF (my_just_energy) THEN
277 : exc = xc_exc_calc(rho_r=rho_r, tau=tau, &
278 : rho_g=rho_g, xc_section=xc_section, &
279 10406 : weights=weights, pw_pool=xc_pw_pool)
280 :
281 : ELSE
282 : CALL xc_vxc_pw_create(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_r, &
283 : rho_g=rho_g, tau=tau, exc=exc, &
284 : xc_section=xc_section, &
285 : weights=weights, pw_pool=xc_pw_pool, &
286 : compute_virial=compute_virial, &
287 117699 : virial_xc=virial%pv_xc)
288 : END IF
289 :
290 : ! remove the nlcc densities (keep stuff in original state)
291 128105 : IF (ASSOCIATED(rho_nlcc)) THEN
292 566 : factor = -1.0_dp
293 1132 : DO ispin = 1, mspin
294 566 : CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
295 1132 : CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
296 : END DO
297 : END IF
298 :
299 : ! calclulate non-local vdW functional
300 : ! only if this XC_SECTION has it
301 : ! if yes, we use the dispersion_env from ks_env
302 : ! this is dangerous, as it assumes a special connection xc_section -> qs_env
303 128105 : IF (vdW_nl) THEN
304 388 : CALL get_ks_env(ks_env=ks_env, para_env=para_env)
305 : ! no SIC functionals allowed
306 388 : CPASSERT(dft_control%sic_method_id == sic_none)
307 : !
308 388 : CALL pw_env_get(pw_env, vdw_pw_pool=vdw_pw_pool)
309 388 : IF (my_just_energy) THEN
310 : CALL calculate_dispersion_nonloc(my_vxc_rho, rho_r, rho_g, edisp, dispersion_env, &
311 6 : my_just_energy, vdw_pw_pool, xc_pw_pool, para_env)
312 : ELSE
313 : CALL calculate_dispersion_nonloc(my_vxc_rho, rho_r, rho_g, edisp, dispersion_env, &
314 382 : my_just_energy, vdw_pw_pool, xc_pw_pool, para_env, virial=virial)
315 : END IF
316 : END IF
317 :
318 : !! Apply rescaling to the potential if requested
319 128105 : IF (.NOT. my_just_energy) THEN
320 117699 : IF (do_adiabatic_rescaling) THEN
321 24 : IF (ASSOCIATED(my_vxc_rho)) THEN
322 62 : DO ispin = 1, SIZE(my_vxc_rho)
323 62 : CALL pw_scale(my_vxc_rho(ispin), my_adiabatic_rescale_factor)
324 : END DO
325 : END IF
326 : END IF
327 : END IF
328 :
329 128105 : IF (my_scaling /= 1.0_dp) THEN
330 92 : exc = exc*my_scaling
331 92 : IF (ASSOCIATED(my_vxc_rho)) THEN
332 180 : DO ispin = 1, SIZE(my_vxc_rho)
333 180 : CALL pw_scale(my_vxc_rho(ispin), my_scaling)
334 : END DO
335 : END IF
336 92 : IF (ASSOCIATED(my_vxc_tau)) THEN
337 0 : DO ispin = 1, SIZE(my_vxc_tau)
338 0 : CALL pw_scale(my_vxc_tau(ispin), my_scaling)
339 : END DO
340 : END IF
341 : END IF
342 :
343 : ! we have pw data for the xc, qs_ks requests coeff structure, here we transfer
344 : ! pw -> coeff
345 128105 : IF (ASSOCIATED(my_vxc_rho)) THEN
346 117699 : vxc_rho => my_vxc_rho
347 117699 : NULLIFY (my_vxc_rho)
348 : END IF
349 128105 : IF (ASSOCIATED(my_vxc_tau)) THEN
350 2176 : vxc_tau => my_vxc_tau
351 2176 : NULLIFY (my_vxc_tau)
352 : END IF
353 128105 : IF (uf_grid) THEN
354 28 : DO ispin = 1, SIZE(rho_r)
355 28 : CALL xc_pw_pool%give_back_pw(rho_r(ispin))
356 : END DO
357 14 : DEALLOCATE (rho_r)
358 14 : IF (ASSOCIATED(rho_g)) THEN
359 28 : DO ispin = 1, SIZE(rho_g)
360 28 : CALL xc_pw_pool%give_back_pw(rho_g(ispin))
361 : END DO
362 14 : DEALLOCATE (rho_g)
363 : END IF
364 : END IF
365 :
366 : ! compute again the xc but now for Exc(m,o) and the opposite sign
367 128105 : IF (dft_control%sic_method_id == sic_mauri_spz .AND. .NOT. sic_scaling_b_zero) THEN
368 390 : ALLOCATE (rho_m_rspace(2), rho_m_gspace(2))
369 78 : CALL xc_pw_pool%create_pw(rho_m_gspace(1))
370 78 : CALL xc_pw_pool%create_pw(rho_m_rspace(1))
371 78 : CALL pw_copy(rho_struct_r(1), rho_m_rspace(1))
372 78 : CALL pw_axpy(rho_struct_r(2), rho_m_rspace(1), alpha=-1._dp)
373 78 : CALL pw_copy(rho_struct_g(1), rho_m_gspace(1))
374 78 : CALL pw_axpy(rho_struct_g(2), rho_m_gspace(1), alpha=-1._dp)
375 : ! bit sad, these will be just zero...
376 78 : CALL xc_pw_pool%create_pw(rho_m_gspace(2))
377 78 : CALL xc_pw_pool%create_pw(rho_m_rspace(2))
378 78 : CALL pw_zero(rho_m_rspace(2))
379 78 : CALL pw_zero(rho_m_gspace(2))
380 :
381 78 : IF (my_just_energy) THEN
382 : exc_m = xc_exc_calc(rho_r=rho_m_rspace, tau=tau, &
383 : rho_g=rho_m_gspace, xc_section=xc_section, &
384 24 : weights=weights, pw_pool=xc_pw_pool)
385 : ELSE
386 : ! virial untested
387 54 : CPASSERT(.NOT. compute_virial)
388 : CALL xc_vxc_pw_create(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_m_rspace, &
389 : rho_g=rho_m_gspace, tau=tau, exc=exc_m, &
390 : xc_section=xc_section, &
391 : weights=weights, pw_pool=xc_pw_pool, &
392 : compute_virial=.FALSE., &
393 54 : virial_xc=virial_xc_tmp)
394 : END IF
395 :
396 78 : exc = exc - dft_control%sic_scaling_b*exc_m
397 :
398 : ! and take care of the potential only vxc_rho is taken into account
399 78 : IF (.NOT. my_just_energy) THEN
400 54 : CALL pw_axpy(my_vxc_rho(1), vxc_rho(1), -dft_control%sic_scaling_b)
401 54 : CALL pw_axpy(my_vxc_rho(1), vxc_rho(2), dft_control%sic_scaling_b)
402 54 : CALL my_vxc_rho(1)%release()
403 54 : CALL my_vxc_rho(2)%release()
404 54 : DEALLOCATE (my_vxc_rho)
405 : END IF
406 :
407 234 : DO ispin = 1, 2
408 156 : CALL xc_pw_pool%give_back_pw(rho_m_rspace(ispin))
409 234 : CALL xc_pw_pool%give_back_pw(rho_m_gspace(ispin))
410 : END DO
411 78 : DEALLOCATE (rho_m_rspace)
412 78 : DEALLOCATE (rho_m_gspace)
413 :
414 : END IF
415 :
416 : ! now we have - sum_s N_s * Exc(rho_s/N_s,0)
417 128105 : IF (dft_control%sic_method_id == sic_ad .AND. .NOT. sic_scaling_b_zero) THEN
418 :
419 : ! find out how many elecs we have
420 26 : CALL get_ks_env(ks_env, nelectron_spin=nelec_spin)
421 :
422 130 : ALLOCATE (rho_m_rspace(2), rho_m_gspace(2))
423 78 : DO ispin = 1, 2
424 52 : CALL xc_pw_pool%create_pw(rho_m_gspace(ispin))
425 78 : CALL xc_pw_pool%create_pw(rho_m_rspace(ispin))
426 : END DO
427 :
428 78 : DO ispin = 1, 2
429 52 : IF (nelec_spin(ispin) > 0.0_dp) THEN
430 52 : nelec_s_inv = 1.0_dp/nelec_spin(ispin)
431 : ELSE
432 : ! does it matter if there are no electrons with this spin (H) ?
433 0 : nelec_s_inv = 0.0_dp
434 : END IF
435 52 : CALL pw_copy(rho_struct_r(ispin), rho_m_rspace(1))
436 52 : CALL pw_copy(rho_struct_g(ispin), rho_m_gspace(1))
437 52 : CALL pw_scale(rho_m_rspace(1), nelec_s_inv)
438 52 : CALL pw_scale(rho_m_gspace(1), nelec_s_inv)
439 52 : CALL pw_zero(rho_m_rspace(2))
440 52 : CALL pw_zero(rho_m_gspace(2))
441 :
442 52 : IF (my_just_energy) THEN
443 : exc_m = xc_exc_calc(rho_r=rho_m_rspace, tau=tau, &
444 : rho_g=rho_m_gspace, xc_section=xc_section, &
445 12 : weights=weights, pw_pool=xc_pw_pool)
446 : ELSE
447 : ! virial untested
448 40 : CPASSERT(.NOT. compute_virial)
449 : CALL xc_vxc_pw_create(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_m_rspace, &
450 : rho_g=rho_m_gspace, tau=tau, exc=exc_m, &
451 : xc_section=xc_section, &
452 : weights=weights, pw_pool=xc_pw_pool, &
453 : compute_virial=.FALSE., &
454 40 : virial_xc=virial_xc_tmp)
455 : END IF
456 :
457 52 : exc = exc - dft_control%sic_scaling_b*nelec_spin(ispin)*exc_m
458 :
459 : ! and take care of the potential only vxc_rho is taken into account
460 78 : IF (.NOT. my_just_energy) THEN
461 40 : CALL pw_axpy(my_vxc_rho(1), vxc_rho(ispin), -dft_control%sic_scaling_b)
462 40 : CALL my_vxc_rho(1)%release()
463 40 : CALL my_vxc_rho(2)%release()
464 40 : DEALLOCATE (my_vxc_rho)
465 : END IF
466 : END DO
467 :
468 78 : DO ispin = 1, 2
469 52 : CALL xc_pw_pool%give_back_pw(rho_m_rspace(ispin))
470 78 : CALL xc_pw_pool%give_back_pw(rho_m_gspace(ispin))
471 : END DO
472 26 : DEALLOCATE (rho_m_rspace)
473 26 : DEALLOCATE (rho_m_gspace)
474 :
475 : END IF
476 :
477 : ! compute again the xc but now for Exc(n_down,n_down)
478 128105 : IF (dft_control%sic_method_id == sic_mauri_us .AND. .NOT. sic_scaling_b_zero) THEN
479 276 : ALLOCATE (rho_r(2))
480 92 : rho_r(1) = rho_struct_r(2)
481 92 : rho_r(2) = rho_struct_r(2)
482 92 : IF (rho_g_valid) THEN
483 276 : ALLOCATE (rho_g(2))
484 92 : rho_g(1) = rho_struct_g(2)
485 92 : rho_g(2) = rho_struct_g(2)
486 : END IF
487 :
488 92 : IF (my_just_energy) THEN
489 : exc_m = xc_exc_calc(rho_r=rho_r, tau=tau, &
490 : rho_g=rho_g, xc_section=xc_section, &
491 32 : weights=weights, pw_pool=xc_pw_pool)
492 : ELSE
493 : ! virial untested
494 60 : CPASSERT(.NOT. compute_virial)
495 : CALL xc_vxc_pw_create(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_r, &
496 : rho_g=rho_g, tau=tau, exc=exc_m, &
497 : xc_section=xc_section, &
498 : weights=weights, pw_pool=xc_pw_pool, &
499 : compute_virial=.FALSE., &
500 60 : virial_xc=virial_xc_tmp)
501 : END IF
502 :
503 92 : exc = exc + dft_control%sic_scaling_b*exc_m
504 :
505 : ! and take care of the potential
506 92 : IF (.NOT. my_just_energy) THEN
507 : ! both go to minority spin
508 60 : CALL pw_axpy(my_vxc_rho(1), vxc_rho(2), 2.0_dp*dft_control%sic_scaling_b)
509 60 : CALL my_vxc_rho(1)%release()
510 60 : CALL my_vxc_rho(2)%release()
511 60 : DEALLOCATE (my_vxc_rho)
512 : END IF
513 92 : DEALLOCATE (rho_r, rho_g)
514 :
515 : END IF
516 :
517 : !
518 : ! cleanups
519 : !
520 128105 : IF (uf_grid .AND. (ASSOCIATED(vxc_rho) .OR. ASSOCIATED(vxc_tau))) THEN
521 : BLOCK
522 : TYPE(pw_r3d_rs_type) :: tmp_pw
523 : TYPE(pw_c1d_gs_type) :: tmp_g, tmp_g2
524 14 : CALL xc_pw_pool%create_pw(tmp_g)
525 14 : CALL auxbas_pw_pool%create_pw(tmp_g2)
526 14 : IF (ASSOCIATED(vxc_rho)) THEN
527 28 : DO ispin = 1, SIZE(vxc_rho)
528 14 : CALL auxbas_pw_pool%create_pw(tmp_pw)
529 14 : CALL pw_transfer(vxc_rho(ispin), tmp_g)
530 14 : CALL pw_transfer(tmp_g, tmp_g2)
531 14 : CALL pw_transfer(tmp_g2, tmp_pw)
532 14 : CALL xc_pw_pool%give_back_pw(vxc_rho(ispin))
533 28 : vxc_rho(ispin) = tmp_pw
534 : END DO
535 : END IF
536 14 : IF (ASSOCIATED(vxc_tau)) THEN
537 0 : DO ispin = 1, SIZE(vxc_tau)
538 0 : CALL auxbas_pw_pool%create_pw(tmp_pw)
539 0 : CALL pw_transfer(vxc_tau(ispin), tmp_g)
540 0 : CALL pw_transfer(tmp_g, tmp_g2)
541 0 : CALL pw_transfer(tmp_g2, tmp_pw)
542 0 : CALL xc_pw_pool%give_back_pw(vxc_tau(ispin))
543 0 : vxc_tau(ispin) = tmp_pw
544 : END DO
545 : END IF
546 14 : CALL auxbas_pw_pool%give_back_pw(tmp_g2)
547 28 : CALL xc_pw_pool%give_back_pw(tmp_g)
548 : END BLOCK
549 : END IF
550 128105 : IF (ASSOCIATED(tau) .AND. uf_grid) THEN
551 0 : DO ispin = 1, SIZE(tau)
552 0 : CALL xc_pw_pool%give_back_pw(tau(ispin))
553 : END DO
554 0 : DEALLOCATE (tau)
555 : END IF
556 :
557 : END IF
558 :
559 142659 : CALL timestop(handle)
560 :
561 142659 : END SUBROUTINE qs_vxc_create
562 :
563 : ! **************************************************************************************************
564 : !> \brief calculates the XC density: E_xc(r) - V_xc(r)*rho(r) or E_xc(r)/rho(r)
565 : !> \param ks_env to get all the needed things
566 : !> \param rho_struct density
567 : !> \param xc_section ...
568 : !> \param dispersion_env ...
569 : !> \param xc_ener will contain the xc energy density E_xc(r) - V_xc(r)*rho(r)
570 : !> \param xc_den will contain the xc energy density E_xc(r)/rho(r)
571 : !> \param exc will contain the xc energy density E_xc(r)
572 : !> \param vxc ...
573 : !> \param vtau ...
574 : !> \author JGH
575 : ! **************************************************************************************************
576 490 : SUBROUTINE qs_xc_density(ks_env, rho_struct, xc_section, dispersion_env, &
577 98 : xc_ener, xc_den, exc, vxc, vtau)
578 :
579 : TYPE(qs_ks_env_type), POINTER :: ks_env
580 : TYPE(qs_rho_type), POINTER :: rho_struct
581 : TYPE(section_vals_type), POINTER :: xc_section
582 : TYPE(qs_dispersion_type), OPTIONAL, POINTER :: dispersion_env
583 : TYPE(pw_r3d_rs_type), INTENT(INOUT), OPTIONAL :: xc_ener, xc_den
584 : TYPE(pw_r3d_rs_type), OPTIONAL :: exc
585 : TYPE(pw_r3d_rs_type), DIMENSION(:), OPTIONAL :: vxc, vtau
586 :
587 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_xc_density'
588 :
589 : INTEGER :: handle, ispin, myfun, nspins, vdw
590 : LOGICAL :: rho_g_valid, tau_r_valid, uf_grid, vdW_nl
591 : REAL(KIND=dp) :: edisp, excint, factor, rho_cutoff
592 : REAL(KIND=dp), DIMENSION(3, 3) :: vdum
593 : TYPE(cell_type), POINTER :: cell
594 : TYPE(dft_control_type), POINTER :: dft_control
595 : TYPE(mp_para_env_type), POINTER :: para_env
596 98 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
597 : TYPE(pw_c1d_gs_type), POINTER :: rho_nlcc_g
598 : TYPE(pw_env_type), POINTER :: pw_env
599 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool, vdw_pw_pool, xc_pw_pool
600 : TYPE(pw_r3d_rs_type) :: exc_r
601 98 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau_r, vxc_rho, vxc_tau
602 : TYPE(pw_r3d_rs_type), POINTER :: rho_nlcc, weights
603 :
604 98 : CALL timeset(routineN, handle)
605 :
606 : CALL get_ks_env(ks_env, &
607 : dft_control=dft_control, &
608 : pw_env=pw_env, &
609 : cell=cell, &
610 : xcint_weights=weights, &
611 : rho_nlcc=rho_nlcc, &
612 98 : rho_nlcc_g=rho_nlcc_g)
613 :
614 : CALL qs_rho_get(rho_struct, &
615 : tau_r_valid=tau_r_valid, &
616 : rho_g_valid=rho_g_valid, &
617 : rho_r=rho_r, &
618 : rho_g=rho_g, &
619 98 : tau_r=tau_r)
620 98 : nspins = dft_control%nspins
621 :
622 98 : CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", i_val=myfun)
623 98 : CALL section_vals_val_get(xc_section, "VDW_POTENTIAL%POTENTIAL_TYPE", i_val=vdw)
624 98 : vdW_nl = (vdw == xc_vdw_fun_nonloc)
625 98 : IF (PRESENT(xc_ener)) THEN
626 32 : IF (tau_r_valid) THEN
627 0 : CALL cp_warn(__LOCATION__, "Tau contribution will not be correctly handled")
628 : END IF
629 : END IF
630 98 : IF (vdW_nl) THEN
631 0 : CALL cp_warn(__LOCATION__, "vdW functional contribution will be ignored")
632 : END IF
633 :
634 98 : CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool, auxbas_pw_pool=auxbas_pw_pool)
635 98 : uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid, xc_pw_pool%pw_grid)
636 98 : IF (uf_grid) THEN
637 0 : CALL cp_warn(__LOCATION__, "Fine grid option not possible with local energy")
638 0 : CPABORT("Fine Grid in Local Energy")
639 : END IF
640 :
641 98 : IF (PRESENT(xc_ener)) THEN
642 32 : CALL pw_zero(xc_ener)
643 : END IF
644 98 : IF (PRESENT(xc_den)) THEN
645 66 : CALL pw_zero(xc_den)
646 : END IF
647 98 : IF (PRESENT(exc)) THEN
648 0 : CALL pw_zero(exc)
649 : END IF
650 98 : IF (PRESENT(vxc)) THEN
651 138 : DO ispin = 1, nspins
652 138 : CALL pw_zero(vxc(ispin))
653 : END DO
654 : END IF
655 98 : IF (PRESENT(vtau)) THEN
656 40 : DO ispin = 1, nspins
657 40 : CALL pw_zero(vtau(ispin))
658 : END DO
659 : END IF
660 :
661 98 : IF (myfun /= xc_none) THEN
662 :
663 96 : CPASSERT(ASSOCIATED(rho_struct))
664 96 : CPASSERT(dft_control%sic_method_id == sic_none)
665 :
666 : ! add the nlcc densities
667 96 : IF (ASSOCIATED(rho_nlcc)) THEN
668 0 : factor = 1.0_dp
669 0 : DO ispin = 1, nspins
670 0 : CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
671 0 : CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
672 : END DO
673 : END IF
674 96 : NULLIFY (vxc_rho, vxc_tau)
675 : CALL xc_vxc_pw_create(vxc_rho=vxc_rho, vxc_tau=vxc_tau, rho_r=rho_r, &
676 : rho_g=rho_g, tau=tau_r, exc=excint, &
677 : xc_section=xc_section, &
678 : weights=weights, pw_pool=xc_pw_pool, &
679 : compute_virial=.FALSE., &
680 : virial_xc=vdum, &
681 96 : exc_r=exc_r)
682 : ! calclulate non-local vdW functional
683 : ! only if this XC_SECTION has it
684 : ! if yes, we use the dispersion_env from ks_env
685 : ! this is dangerous, as it assumes a special connection xc_section -> qs_env
686 96 : IF (vdW_nl) THEN
687 0 : CALL get_ks_env(ks_env=ks_env, para_env=para_env)
688 : ! no SIC functionals allowed
689 0 : CPASSERT(dft_control%sic_method_id == sic_none)
690 : !
691 0 : CALL pw_env_get(pw_env, vdw_pw_pool=vdw_pw_pool)
692 : CALL calculate_dispersion_nonloc(vxc_rho, rho_r, rho_g, edisp, dispersion_env, &
693 0 : .FALSE., vdw_pw_pool, xc_pw_pool, para_env)
694 : END IF
695 :
696 : ! remove the nlcc densities (keep stuff in original state)
697 96 : IF (ASSOCIATED(rho_nlcc)) THEN
698 0 : factor = -1.0_dp
699 0 : DO ispin = 1, dft_control%nspins
700 0 : CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
701 0 : CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
702 : END DO
703 : END IF
704 : !
705 96 : IF (PRESENT(xc_den)) THEN
706 64 : CALL pw_copy(exc_r, xc_den)
707 64 : rho_cutoff = 1.E-14_dp
708 64 : CALL calc_xc_density(xc_den, rho_r, rho_cutoff)
709 : END IF
710 96 : IF (PRESENT(xc_ener)) THEN
711 32 : CALL pw_copy(exc_r, xc_ener)
712 64 : DO ispin = 1, nspins
713 64 : CALL pw_multiply(xc_ener, vxc_rho(ispin), rho_r(ispin), alpha=-1.0_dp)
714 : END DO
715 : END IF
716 96 : IF (PRESENT(exc)) THEN
717 0 : CALL pw_copy(exc_r, exc)
718 : END IF
719 96 : IF (PRESENT(vxc)) THEN
720 134 : DO ispin = 1, nspins
721 134 : CALL pw_copy(vxc_rho(ispin), vxc(ispin))
722 : END DO
723 : END IF
724 96 : IF (PRESENT(vtau)) THEN
725 40 : DO ispin = 1, nspins
726 40 : CALL pw_copy(vxc_tau(ispin), vtau(ispin))
727 : END DO
728 : END IF
729 : ! remove arrays
730 96 : IF (ASSOCIATED(vxc_rho)) THEN
731 198 : DO ispin = 1, nspins
732 198 : CALL vxc_rho(ispin)%release()
733 : END DO
734 96 : DEALLOCATE (vxc_rho)
735 : END IF
736 96 : IF (ASSOCIATED(vxc_tau)) THEN
737 40 : DO ispin = 1, nspins
738 40 : CALL vxc_tau(ispin)%release()
739 : END DO
740 20 : DEALLOCATE (vxc_tau)
741 : END IF
742 96 : CALL exc_r%release()
743 : !
744 : END IF
745 :
746 98 : CALL timestop(handle)
747 :
748 98 : END SUBROUTINE qs_xc_density
749 :
750 : END MODULE qs_vxc
|
