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 Calculate the KS reference potentials
10 : !> \par History
11 : !> 07.2022 created
12 : !> \author JGH
13 : ! **************************************************************************************************
14 : MODULE qs_ks_reference
15 : USE admm_types, ONLY: admm_type,&
16 : get_admm_env
17 : USE atomic_kind_types, ONLY: atomic_kind_type
18 : USE cp_control_types, ONLY: dft_control_type
19 : USE cp_dbcsr_api, ONLY: dbcsr_p_type
20 : USE hartree_local_methods, ONLY: Vh_1c_gg_integrals,&
21 : init_coulomb_local
22 : USE hartree_local_types, ONLY: hartree_local_create,&
23 : hartree_local_release,&
24 : hartree_local_type
25 : USE input_constants, ONLY: do_admm_aux_exch_func_none
26 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
27 : section_vals_type
28 : USE kinds, ONLY: dp
29 : USE message_passing, ONLY: mp_para_env_type
30 : USE pw_env_types, ONLY: pw_env_get,&
31 : pw_env_type
32 : USE pw_grid_types, ONLY: pw_grid_type
33 : USE pw_methods, ONLY: pw_scale,&
34 : pw_transfer,&
35 : pw_zero
36 : USE pw_poisson_methods, ONLY: pw_poisson_solve
37 : USE pw_poisson_types, ONLY: pw_poisson_type
38 : USE pw_pool_types, ONLY: pw_pool_type
39 : USE pw_types, ONLY: pw_c1d_gs_type,&
40 : pw_r3d_rs_type
41 : USE qs_core_energies, ONLY: calculate_ecore_overlap,&
42 : calculate_ecore_self
43 : USE qs_environment_types, ONLY: get_qs_env,&
44 : qs_environment_type
45 : USE qs_gapw_densities, ONLY: prepare_gapw_den
46 : USE qs_kind_types, ONLY: qs_kind_type
47 : USE qs_ks_methods, ONLY: calc_rho_tot_gspace
48 : USE qs_ks_types, ONLY: qs_ks_env_type
49 : USE qs_local_rho_types, ONLY: local_rho_set_create,&
50 : local_rho_type
51 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
52 : USE qs_oce_types, ONLY: oce_matrix_type
53 : USE qs_rho0_ggrid, ONLY: integrate_vhg0_rspace,&
54 : rho0_s_grid_create
55 : USE qs_rho0_methods, ONLY: init_rho0
56 : USE qs_rho_atom_methods, ONLY: allocate_rho_atom_internals,&
57 : calculate_rho_atom_coeff
58 : USE qs_rho_types, ONLY: qs_rho_get,&
59 : qs_rho_type
60 : USE qs_vxc, ONLY: qs_vxc_create
61 : USE qs_vxc_atom, ONLY: calculate_vxc_atom
62 : USE virial_types, ONLY: virial_type
63 : #include "./base/base_uses.f90"
64 :
65 : IMPLICIT NONE
66 :
67 : PRIVATE
68 :
69 : ! *** Global parameters ***
70 :
71 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ks_reference'
72 :
73 : PUBLIC :: ks_ref_potential, ks_ref_potential_atom
74 :
75 : ! **************************************************************************************************
76 :
77 : CONTAINS
78 :
79 : ! **************************************************************************************************
80 : !> \brief calculate the Kohn-Sham reference potential
81 : !> \param qs_env ...
82 : !> \param vh_rspace ...
83 : !> \param vxc_rspace ...
84 : !> \param vtau_rspace ...
85 : !> \param vadmm_rspace ...
86 : !> \param ehartree ...
87 : !> \param exc ...
88 : !> \param h_stress container for the stress tensor of the Hartree term
89 : !> \param vadmm_tau_rspace ADMM auxiliary tau potential
90 : !> \par History
91 : !> 10.2019 created [JGH]
92 : !> \author JGH
93 : ! **************************************************************************************************
94 5892 : SUBROUTINE ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, &
95 : ehartree, exc, h_stress, vadmm_tau_rspace)
96 : TYPE(qs_environment_type), POINTER :: qs_env
97 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: vh_rspace
98 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_rspace, vtau_rspace, vadmm_rspace
99 : REAL(KIND=dp), INTENT(OUT) :: ehartree, exc
100 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), &
101 : OPTIONAL :: h_stress
102 : TYPE(pw_r3d_rs_type), DIMENSION(:), OPTIONAL, &
103 : POINTER :: vadmm_tau_rspace
104 :
105 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ks_ref_potential'
106 :
107 : INTEGER :: handle, iab, ispin, nspins
108 : REAL(dp) :: eadmm, eovrl, eself
109 : REAL(KIND=dp), DIMENSION(3, 3) :: virial_xc
110 : TYPE(admm_type), POINTER :: admm_env
111 : TYPE(dft_control_type), POINTER :: dft_control
112 : TYPE(mp_para_env_type), POINTER :: para_env
113 : TYPE(pw_c1d_gs_type) :: rho_tot_gspace, v_hartree_gspace
114 : TYPE(pw_c1d_gs_type), POINTER :: rho_core
115 : TYPE(pw_env_type), POINTER :: pw_env
116 : TYPE(pw_grid_type), POINTER :: pw_grid
117 : TYPE(pw_poisson_type), POINTER :: poisson_env
118 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
119 : TYPE(pw_r3d_rs_type) :: v_hartree_rspace
120 1964 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_admm_rspace, v_admm_tau_rspace, &
121 1964 : v_rspace, v_tau_rspace
122 : TYPE(qs_ks_env_type), POINTER :: ks_env
123 : TYPE(qs_rho_type), POINTER :: rho, rho_xc
124 : TYPE(section_vals_type), POINTER :: xc_section
125 : TYPE(virial_type), POINTER :: virial
126 :
127 1964 : CALL timeset(routineN, handle)
128 :
129 : ! get all information on the electronic density
130 1964 : NULLIFY (rho, ks_env)
131 : CALL get_qs_env(qs_env=qs_env, rho=rho, dft_control=dft_control, &
132 1964 : para_env=para_env, ks_env=ks_env, rho_core=rho_core)
133 :
134 1964 : nspins = dft_control%nspins
135 :
136 1964 : NULLIFY (pw_env)
137 1964 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
138 1964 : CPASSERT(ASSOCIATED(pw_env))
139 :
140 1964 : NULLIFY (auxbas_pw_pool, poisson_env)
141 : ! gets the tmp grids
142 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
143 1964 : poisson_env=poisson_env)
144 :
145 : ! Calculate the Hartree potential
146 1964 : CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
147 1964 : CALL auxbas_pw_pool%create_pw(v_hartree_rspace)
148 1964 : CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
149 :
150 : ! Get the total density in g-space [ions + electrons]
151 1964 : CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho)
152 :
153 : CALL pw_poisson_solve(poisson_env, rho_tot_gspace, ehartree, &
154 1964 : v_hartree_gspace, h_stress=h_stress, rho_core=rho_core)
155 1964 : CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
156 1964 : CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
157 :
158 1964 : CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
159 1964 : CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
160 : !
161 1964 : CALL calculate_ecore_self(qs_env, E_self_core=eself)
162 1964 : CALL calculate_ecore_overlap(qs_env, para_env, PRESENT(h_stress), E_overlap_core=eovrl)
163 1964 : ehartree = ehartree + eovrl + eself
164 :
165 : ! v_rspace and v_tau_rspace are generated from the auxbas pool
166 1964 : IF (dft_control%do_admm) THEN
167 420 : CALL get_qs_env(qs_env, admm_env=admm_env)
168 420 : xc_section => admm_env%xc_section_primary
169 : ELSE
170 1544 : xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC")
171 : END IF
172 1964 : NULLIFY (v_rspace, v_tau_rspace)
173 1964 : IF (dft_control%qs_control%gapw_xc) THEN
174 92 : CALL get_qs_env(qs_env=qs_env, rho_xc=rho_xc)
175 : CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_xc, xc_section=xc_section, &
176 92 : vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE.)
177 : ELSE
178 : CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=xc_section, &
179 1872 : vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE.)
180 : END IF
181 :
182 1964 : NULLIFY (v_admm_rspace, v_admm_tau_rspace)
183 1964 : IF (dft_control%do_admm) THEN
184 420 : IF (dft_control%admm_control%aux_exch_func /= do_admm_aux_exch_func_none) THEN
185 : ! For the virial, we have to save the pv_xc component because it will be reset in qs_vxc_create
186 248 : IF (PRESENT(h_stress)) THEN
187 12 : CALL get_qs_env(qs_env, virial=virial)
188 156 : virial_xc = virial%pv_xc
189 : END IF
190 248 : CALL get_admm_env(admm_env, rho_aux_fit=rho)
191 248 : xc_section => admm_env%xc_section_aux
192 : CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=xc_section, &
193 248 : vxc_rho=v_admm_rspace, vxc_tau=v_admm_tau_rspace, exc=eadmm, just_energy=.FALSE.)
194 392 : IF (PRESENT(h_stress)) virial%pv_xc = virial%pv_xc + virial_xc
195 : END IF
196 : END IF
197 :
198 : ! allocate potentials
199 1964 : IF (ASSOCIATED(vh_rspace%pw_grid)) THEN
200 0 : CALL vh_rspace%release()
201 : END IF
202 1964 : IF (ASSOCIATED(vxc_rspace)) THEN
203 0 : DO iab = 1, SIZE(vxc_rspace)
204 0 : CALL vxc_rspace(iab)%release()
205 : END DO
206 : ELSE
207 8078 : ALLOCATE (vxc_rspace(nspins))
208 : END IF
209 1964 : IF (ASSOCIATED(v_tau_rspace)) THEN
210 96 : IF (ASSOCIATED(vtau_rspace)) THEN
211 0 : DO iab = 1, SIZE(vtau_rspace)
212 0 : CALL vtau_rspace(iab)%release()
213 : END DO
214 : ELSE
215 404 : ALLOCATE (vtau_rspace(nspins))
216 : END IF
217 : ELSE
218 1868 : NULLIFY (vtau_rspace)
219 : END IF
220 1964 : IF (ASSOCIATED(v_admm_rspace)) THEN
221 240 : IF (ASSOCIATED(vadmm_rspace)) THEN
222 0 : DO iab = 1, SIZE(vadmm_rspace)
223 0 : CALL vadmm_rspace(iab)%release()
224 : END DO
225 : ELSE
226 988 : ALLOCATE (vadmm_rspace(nspins))
227 : END IF
228 : ELSE
229 1724 : NULLIFY (vadmm_rspace)
230 : END IF
231 1964 : IF (PRESENT(vadmm_tau_rspace)) THEN
232 1338 : IF (ASSOCIATED(vadmm_tau_rspace)) THEN
233 0 : DO iab = 1, SIZE(vadmm_tau_rspace)
234 0 : CALL vadmm_tau_rspace(iab)%release()
235 : END DO
236 : END IF
237 1338 : IF (ASSOCIATED(v_admm_tau_rspace)) THEN
238 0 : IF (.NOT. ASSOCIATED(vadmm_tau_rspace)) ALLOCATE (vadmm_tau_rspace(nspins))
239 : ELSE
240 1338 : IF (ASSOCIATED(vadmm_tau_rspace)) DEALLOCATE (vadmm_tau_rspace)
241 1338 : NULLIFY (vadmm_tau_rspace)
242 : END IF
243 : END IF
244 :
245 1964 : pw_grid => v_hartree_rspace%pw_grid
246 1964 : CALL vh_rspace%create(pw_grid)
247 4150 : DO ispin = 1, nspins
248 2186 : CALL vxc_rspace(ispin)%create(pw_grid)
249 2186 : IF (ASSOCIATED(vtau_rspace)) THEN
250 116 : CALL vtau_rspace(ispin)%create(pw_grid)
251 : END IF
252 2186 : IF (ASSOCIATED(vadmm_rspace)) THEN
253 268 : CALL vadmm_rspace(ispin)%create(pw_grid)
254 : END IF
255 4150 : IF (PRESENT(vadmm_tau_rspace)) THEN
256 1452 : IF (ASSOCIATED(vadmm_tau_rspace)) CALL vadmm_tau_rspace(ispin)%create(pw_grid)
257 : END IF
258 : END DO
259 : !
260 1964 : CALL pw_transfer(v_hartree_rspace, vh_rspace)
261 1964 : IF (ASSOCIATED(v_rspace)) THEN
262 3252 : DO ispin = 1, nspins
263 1712 : CALL pw_transfer(v_rspace(ispin), vxc_rspace(ispin))
264 1712 : CALL pw_scale(vxc_rspace(ispin), v_rspace(ispin)%pw_grid%dvol)
265 3252 : IF (ASSOCIATED(v_tau_rspace)) THEN
266 116 : CALL pw_transfer(v_tau_rspace(ispin), vtau_rspace(ispin))
267 116 : CALL pw_scale(vtau_rspace(ispin), v_tau_rspace(ispin)%pw_grid%dvol)
268 : END IF
269 : END DO
270 : ELSE
271 898 : DO ispin = 1, nspins
272 898 : CALL pw_zero(vxc_rspace(ispin))
273 : END DO
274 : END IF
275 1964 : IF (ASSOCIATED(v_admm_rspace)) THEN
276 508 : DO ispin = 1, nspins
277 268 : CALL pw_transfer(v_admm_rspace(ispin), vadmm_rspace(ispin))
278 508 : CALL pw_scale(vadmm_rspace(ispin), vadmm_rspace(ispin)%pw_grid%dvol)
279 : END DO
280 : END IF
281 1964 : IF (PRESENT(vadmm_tau_rspace)) THEN
282 1338 : IF (ASSOCIATED(v_admm_tau_rspace)) THEN
283 0 : DO ispin = 1, nspins
284 0 : CALL pw_transfer(v_admm_tau_rspace(ispin), vadmm_tau_rspace(ispin))
285 0 : CALL pw_scale(vadmm_tau_rspace(ispin), vadmm_tau_rspace(ispin)%pw_grid%dvol)
286 : END DO
287 : END IF
288 : END IF
289 :
290 : ! return pw grids
291 1964 : CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
292 1964 : IF (ASSOCIATED(v_rspace)) THEN
293 3252 : DO ispin = 1, nspins
294 1712 : CALL auxbas_pw_pool%give_back_pw(v_rspace(ispin))
295 3252 : IF (ASSOCIATED(v_tau_rspace)) THEN
296 116 : CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
297 : END IF
298 : END DO
299 1540 : DEALLOCATE (v_rspace)
300 : END IF
301 1964 : IF (ASSOCIATED(v_tau_rspace)) DEALLOCATE (v_tau_rspace)
302 1964 : IF (ASSOCIATED(v_admm_rspace)) THEN
303 508 : DO ispin = 1, nspins
304 508 : CALL auxbas_pw_pool%give_back_pw(v_admm_rspace(ispin))
305 : END DO
306 240 : DEALLOCATE (v_admm_rspace)
307 : END IF
308 1964 : IF (ASSOCIATED(v_admm_tau_rspace)) THEN
309 0 : DO ispin = 1, nspins
310 0 : CALL auxbas_pw_pool%give_back_pw(v_admm_tau_rspace(ispin))
311 : END DO
312 0 : DEALLOCATE (v_admm_tau_rspace)
313 : END IF
314 :
315 1964 : CALL timestop(handle)
316 :
317 1964 : END SUBROUTINE ks_ref_potential
318 :
319 : ! **************************************************************************************************
320 : !> \brief calculate the Kohn-Sham GAPW reference potentials
321 : !> \param qs_env ...
322 : !> \param local_rho_set ...
323 : !> \param local_rho_set_admm ...
324 : !> \param v_hartree_rspace ...
325 : !> \par History
326 : !> 07.2022 created [JGH]
327 : !> \author JGH
328 : ! **************************************************************************************************
329 1338 : SUBROUTINE ks_ref_potential_atom(qs_env, local_rho_set, local_rho_set_admm, v_hartree_rspace)
330 : TYPE(qs_environment_type), POINTER :: qs_env
331 : TYPE(local_rho_type), POINTER :: local_rho_set, local_rho_set_admm
332 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v_hartree_rspace
333 :
334 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ks_ref_potential_atom'
335 :
336 : INTEGER :: handle, natom, nspins
337 : LOGICAL :: gapw, gapw_xc
338 : REAL(KIND=dp) :: eh1c, exc1, exc1_admm
339 : TYPE(admm_type), POINTER :: admm_env
340 1338 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
341 1338 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_aux, rho_ao_kp
342 : TYPE(dft_control_type), POINTER :: dft_control
343 : TYPE(hartree_local_type), POINTER :: hartree_local
344 : TYPE(mp_para_env_type), POINTER :: para_env
345 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
346 1338 : POINTER :: sab
347 : TYPE(oce_matrix_type), POINTER :: oce
348 : TYPE(pw_env_type), POINTER :: pw_env
349 1338 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
350 : TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit
351 : TYPE(section_vals_type), POINTER :: xc_section
352 :
353 1338 : CALL timeset(routineN, handle)
354 :
355 : ! get all information on the electronic density
356 : CALL get_qs_env(qs_env=qs_env, rho=rho, pw_env=pw_env, &
357 1338 : dft_control=dft_control, para_env=para_env)
358 :
359 1338 : nspins = dft_control%nspins
360 1338 : gapw = dft_control%qs_control%gapw
361 1338 : gapw_xc = dft_control%qs_control%gapw_xc
362 :
363 1338 : IF (gapw .OR. gapw_xc) THEN
364 302 : NULLIFY (hartree_local, local_rho_set, local_rho_set_admm)
365 : CALL get_qs_env(qs_env, &
366 : atomic_kind_set=atomic_kind_set, &
367 302 : qs_kind_set=qs_kind_set)
368 302 : CALL local_rho_set_create(local_rho_set)
369 : CALL allocate_rho_atom_internals(local_rho_set%rho_atom_set, atomic_kind_set, &
370 302 : qs_kind_set, dft_control, para_env)
371 302 : IF (gapw) THEN
372 234 : CALL get_qs_env(qs_env, natom=natom)
373 234 : CALL init_rho0(local_rho_set, qs_env, dft_control%qs_control%gapw_control)
374 234 : CALL rho0_s_grid_create(pw_env, local_rho_set%rho0_mpole)
375 234 : CALL hartree_local_create(hartree_local)
376 234 : CALL init_coulomb_local(hartree_local, natom)
377 : END IF
378 :
379 302 : CALL get_qs_env(qs_env=qs_env, oce=oce, sab_orb=sab)
380 302 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
381 : CALL calculate_rho_atom_coeff(qs_env, rho_ao_kp, local_rho_set%rho_atom_set, &
382 302 : qs_kind_set, oce, sab, para_env)
383 302 : CALL prepare_gapw_den(qs_env, local_rho_set, do_rho0=gapw)
384 :
385 302 : IF (dft_control%do_admm) THEN
386 54 : CALL get_qs_env(qs_env, admm_env=admm_env)
387 54 : xc_section => admm_env%xc_section_primary
388 : ELSE
389 248 : xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC")
390 : END IF
391 : CALL calculate_vxc_atom(qs_env, .FALSE., exc1=exc1, xc_section_external=xc_section, &
392 302 : rho_atom_set_external=local_rho_set%rho_atom_set)
393 :
394 302 : IF (gapw) THEN
395 234 : CALL Vh_1c_gg_integrals(qs_env, eh1c, hartree_local%ecoul_1c, local_rho_set, para_env, .FALSE.)
396 : CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, para_env, calculate_forces=.FALSE., &
397 234 : local_rho_set=local_rho_set)
398 : END IF
399 :
400 302 : IF (dft_control%do_admm) THEN
401 54 : IF (admm_env%do_gapw) THEN
402 54 : CALL local_rho_set_create(local_rho_set_admm)
403 : CALL allocate_rho_atom_internals(local_rho_set_admm%rho_atom_set, atomic_kind_set, &
404 54 : admm_env%admm_gapw_env%admm_kind_set, dft_control, para_env)
405 54 : oce => admm_env%admm_gapw_env%oce
406 54 : sab => admm_env%sab_aux_fit
407 54 : CALL get_admm_env(admm_env, rho_aux_fit=rho_aux_fit)
408 54 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_aux)
409 : CALL calculate_rho_atom_coeff(qs_env, rho_ao_aux, local_rho_set_admm%rho_atom_set, &
410 54 : admm_env%admm_gapw_env%admm_kind_set, oce, sab, para_env)
411 : CALL prepare_gapw_den(qs_env, local_rho_set=local_rho_set_admm, &
412 54 : do_rho0=.FALSE., kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
413 : !compute the potential due to atomic densities
414 54 : xc_section => admm_env%xc_section_aux
415 : CALL calculate_vxc_atom(qs_env, energy_only=.FALSE., exc1=exc1_admm, &
416 : kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
417 : xc_section_external=xc_section, &
418 54 : rho_atom_set_external=local_rho_set_admm%rho_atom_set)
419 : END IF
420 : END IF
421 :
422 : ! clean up
423 302 : CALL hartree_local_release(hartree_local)
424 :
425 : ELSE
426 :
427 1036 : NULLIFY (local_rho_set, local_rho_set_admm)
428 :
429 : END IF
430 :
431 1338 : CALL timestop(handle)
432 :
433 1338 : END SUBROUTINE ks_ref_potential_atom
434 :
435 : ! **************************************************************************************************
436 :
437 : END MODULE qs_ks_reference
|
