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