Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 : !> \par History
90 : !> 10.2019 created [JGH]
91 : !> \author JGH
92 : ! **************************************************************************************************
93 5370 : SUBROUTINE ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, &
94 : ehartree, exc, h_stress)
95 : TYPE(qs_environment_type), POINTER :: qs_env
96 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: vh_rspace
97 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_rspace, vtau_rspace, vadmm_rspace
98 : REAL(KIND=dp), INTENT(OUT) :: ehartree, exc
99 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT), &
100 : OPTIONAL :: h_stress
101 :
102 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ks_ref_potential'
103 :
104 : INTEGER :: handle, iab, ispin, nspins
105 : REAL(dp) :: eadmm, eovrl, eself
106 : REAL(KIND=dp), DIMENSION(3, 3) :: virial_xc
107 : TYPE(admm_type), POINTER :: admm_env
108 : TYPE(dft_control_type), POINTER :: dft_control
109 : TYPE(mp_para_env_type), POINTER :: para_env
110 : TYPE(pw_c1d_gs_type) :: rho_tot_gspace, v_hartree_gspace
111 : TYPE(pw_c1d_gs_type), POINTER :: rho_core
112 : TYPE(pw_env_type), POINTER :: pw_env
113 : TYPE(pw_grid_type), POINTER :: pw_grid
114 : TYPE(pw_poisson_type), POINTER :: poisson_env
115 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
116 : TYPE(pw_r3d_rs_type) :: v_hartree_rspace
117 1790 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_admm_rspace, v_admm_tau_rspace, &
118 1790 : v_rspace, v_tau_rspace
119 : TYPE(qs_ks_env_type), POINTER :: ks_env
120 : TYPE(qs_rho_type), POINTER :: rho, rho_xc
121 : TYPE(section_vals_type), POINTER :: xc_section
122 : TYPE(virial_type), POINTER :: virial
123 :
124 1790 : CALL timeset(routineN, handle)
125 :
126 : ! get all information on the electronic density
127 1790 : NULLIFY (rho, ks_env)
128 : CALL get_qs_env(qs_env=qs_env, rho=rho, dft_control=dft_control, &
129 1790 : para_env=para_env, ks_env=ks_env, rho_core=rho_core)
130 :
131 1790 : nspins = dft_control%nspins
132 :
133 1790 : NULLIFY (pw_env)
134 1790 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
135 1790 : CPASSERT(ASSOCIATED(pw_env))
136 :
137 1790 : NULLIFY (auxbas_pw_pool, poisson_env)
138 : ! gets the tmp grids
139 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
140 1790 : poisson_env=poisson_env)
141 :
142 : ! Calculate the Hartree potential
143 1790 : CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
144 1790 : CALL auxbas_pw_pool%create_pw(v_hartree_rspace)
145 1790 : CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
146 :
147 : ! Get the total density in g-space [ions + electrons]
148 1790 : CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho)
149 :
150 : CALL pw_poisson_solve(poisson_env, rho_tot_gspace, ehartree, &
151 1790 : v_hartree_gspace, h_stress=h_stress, rho_core=rho_core)
152 1790 : CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
153 1790 : CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
154 :
155 1790 : CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
156 1790 : CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
157 : !
158 1790 : CALL calculate_ecore_self(qs_env, E_self_core=eself)
159 1790 : CALL calculate_ecore_overlap(qs_env, para_env, PRESENT(h_stress), E_overlap_core=eovrl)
160 1790 : ehartree = ehartree + eovrl + eself
161 :
162 : ! v_rspace and v_tau_rspace are generated from the auxbas pool
163 1790 : IF (dft_control%do_admm) THEN
164 390 : CALL get_qs_env(qs_env, admm_env=admm_env)
165 390 : xc_section => admm_env%xc_section_primary
166 : ELSE
167 1400 : xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC")
168 : END IF
169 1790 : NULLIFY (v_rspace, v_tau_rspace)
170 1790 : IF (dft_control%qs_control%gapw_xc) THEN
171 80 : CALL get_qs_env(qs_env=qs_env, rho_xc=rho_xc)
172 : CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_xc, xc_section=xc_section, &
173 80 : vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE.)
174 : ELSE
175 : CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=xc_section, &
176 1710 : vxc_rho=v_rspace, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE.)
177 : END IF
178 :
179 1790 : NULLIFY (v_admm_rspace, v_admm_tau_rspace)
180 1790 : IF (dft_control%do_admm) THEN
181 390 : IF (dft_control%admm_control%aux_exch_func /= do_admm_aux_exch_func_none) THEN
182 : ! For the virial, we have to save the pv_xc component because it will be reset in qs_vxc_create
183 240 : IF (PRESENT(h_stress)) THEN
184 12 : CALL get_qs_env(qs_env, virial=virial)
185 156 : virial_xc = virial%pv_xc
186 : END IF
187 240 : CALL get_admm_env(admm_env, rho_aux_fit=rho)
188 240 : xc_section => admm_env%xc_section_aux
189 : CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=xc_section, &
190 240 : vxc_rho=v_admm_rspace, vxc_tau=v_admm_tau_rspace, exc=eadmm, just_energy=.FALSE.)
191 384 : IF (PRESENT(h_stress)) virial%pv_xc = virial%pv_xc + virial_xc
192 : END IF
193 : END IF
194 :
195 : ! allocate potentials
196 1790 : IF (ASSOCIATED(vh_rspace%pw_grid)) THEN
197 0 : CALL vh_rspace%release()
198 : END IF
199 1790 : IF (ASSOCIATED(vxc_rspace)) THEN
200 0 : DO iab = 1, SIZE(vxc_rspace)
201 0 : CALL vxc_rspace(iab)%release()
202 : END DO
203 : ELSE
204 7374 : ALLOCATE (vxc_rspace(nspins))
205 : END IF
206 1790 : IF (ASSOCIATED(v_tau_rspace)) THEN
207 96 : IF (ASSOCIATED(vtau_rspace)) THEN
208 0 : DO iab = 1, SIZE(vtau_rspace)
209 0 : CALL vtau_rspace(iab)%release()
210 : END DO
211 : ELSE
212 404 : ALLOCATE (vtau_rspace(nspins))
213 : END IF
214 : ELSE
215 1694 : NULLIFY (vtau_rspace)
216 : END IF
217 1790 : IF (ASSOCIATED(v_admm_rspace)) THEN
218 232 : IF (ASSOCIATED(vadmm_rspace)) THEN
219 0 : DO iab = 1, SIZE(vadmm_rspace)
220 0 : CALL vadmm_rspace(iab)%release()
221 : END DO
222 : ELSE
223 956 : ALLOCATE (vadmm_rspace(nspins))
224 : END IF
225 : ELSE
226 1558 : NULLIFY (vadmm_rspace)
227 : END IF
228 :
229 1790 : pw_grid => v_hartree_rspace%pw_grid
230 1790 : CALL vh_rspace%create(pw_grid)
231 3794 : DO ispin = 1, nspins
232 2004 : CALL vxc_rspace(ispin)%create(pw_grid)
233 2004 : IF (ASSOCIATED(vtau_rspace)) THEN
234 116 : CALL vtau_rspace(ispin)%create(pw_grid)
235 : END IF
236 3794 : IF (ASSOCIATED(vadmm_rspace)) THEN
237 260 : CALL vadmm_rspace(ispin)%create(pw_grid)
238 : END IF
239 : END DO
240 : !
241 1790 : CALL pw_transfer(v_hartree_rspace, vh_rspace)
242 1790 : IF (ASSOCIATED(v_rspace)) THEN
243 2960 : DO ispin = 1, nspins
244 1562 : CALL pw_transfer(v_rspace(ispin), vxc_rspace(ispin))
245 1562 : CALL pw_scale(vxc_rspace(ispin), v_rspace(ispin)%pw_grid%dvol)
246 2960 : IF (ASSOCIATED(v_tau_rspace)) THEN
247 116 : CALL pw_transfer(v_tau_rspace(ispin), vtau_rspace(ispin))
248 116 : CALL pw_scale(vtau_rspace(ispin), v_tau_rspace(ispin)%pw_grid%dvol)
249 : END IF
250 : END DO
251 : ELSE
252 834 : DO ispin = 1, nspins
253 834 : CALL pw_zero(vxc_rspace(ispin))
254 : END DO
255 : END IF
256 1790 : IF (ASSOCIATED(v_admm_rspace)) THEN
257 492 : DO ispin = 1, nspins
258 260 : CALL pw_transfer(v_admm_rspace(ispin), vadmm_rspace(ispin))
259 492 : CALL pw_scale(vadmm_rspace(ispin), vadmm_rspace(ispin)%pw_grid%dvol)
260 : END DO
261 : END IF
262 :
263 : ! return pw grids
264 1790 : CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
265 1790 : IF (ASSOCIATED(v_rspace)) THEN
266 2960 : DO ispin = 1, nspins
267 1562 : CALL auxbas_pw_pool%give_back_pw(v_rspace(ispin))
268 2960 : IF (ASSOCIATED(v_tau_rspace)) THEN
269 116 : CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
270 : END IF
271 : END DO
272 1398 : DEALLOCATE (v_rspace)
273 : END IF
274 1790 : IF (ASSOCIATED(v_tau_rspace)) DEALLOCATE (v_tau_rspace)
275 1790 : IF (ASSOCIATED(v_admm_rspace)) THEN
276 492 : DO ispin = 1, nspins
277 492 : CALL auxbas_pw_pool%give_back_pw(v_admm_rspace(ispin))
278 : END DO
279 232 : DEALLOCATE (v_admm_rspace)
280 : END IF
281 :
282 1790 : CALL timestop(handle)
283 :
284 1790 : END SUBROUTINE ks_ref_potential
285 :
286 : ! **************************************************************************************************
287 : !> \brief calculate the Kohn-Sham GAPW reference potentials
288 : !> \param qs_env ...
289 : !> \param local_rho_set ...
290 : !> \param local_rho_set_admm ...
291 : !> \param v_hartree_rspace ...
292 : !> \par History
293 : !> 07.2022 created [JGH]
294 : !> \author JGH
295 : ! **************************************************************************************************
296 1190 : SUBROUTINE ks_ref_potential_atom(qs_env, local_rho_set, local_rho_set_admm, v_hartree_rspace)
297 : TYPE(qs_environment_type), POINTER :: qs_env
298 : TYPE(local_rho_type), POINTER :: local_rho_set, local_rho_set_admm
299 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v_hartree_rspace
300 :
301 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ks_ref_potential_atom'
302 :
303 : INTEGER :: handle, natom, nspins
304 : LOGICAL :: gapw, gapw_xc
305 : REAL(KIND=dp) :: eh1c, exc1, exc1_admm
306 : TYPE(admm_type), POINTER :: admm_env
307 1190 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
308 1190 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_aux, rho_ao_kp
309 : TYPE(dft_control_type), POINTER :: dft_control
310 : TYPE(hartree_local_type), POINTER :: hartree_local
311 : TYPE(mp_para_env_type), POINTER :: para_env
312 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
313 1190 : POINTER :: sab
314 : TYPE(oce_matrix_type), POINTER :: oce
315 : TYPE(pw_env_type), POINTER :: pw_env
316 1190 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
317 : TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit
318 : TYPE(section_vals_type), POINTER :: xc_section
319 :
320 1190 : CALL timeset(routineN, handle)
321 :
322 : ! get all information on the electronic density
323 : CALL get_qs_env(qs_env=qs_env, rho=rho, pw_env=pw_env, &
324 1190 : dft_control=dft_control, para_env=para_env)
325 :
326 1190 : nspins = dft_control%nspins
327 1190 : gapw = dft_control%qs_control%gapw
328 1190 : gapw_xc = dft_control%qs_control%gapw_xc
329 :
330 1190 : IF (gapw .OR. gapw_xc) THEN
331 174 : NULLIFY (hartree_local, local_rho_set, local_rho_set_admm)
332 : CALL get_qs_env(qs_env, &
333 : atomic_kind_set=atomic_kind_set, &
334 174 : qs_kind_set=qs_kind_set)
335 174 : CALL local_rho_set_create(local_rho_set)
336 : CALL allocate_rho_atom_internals(local_rho_set%rho_atom_set, atomic_kind_set, &
337 174 : qs_kind_set, dft_control, para_env)
338 174 : IF (gapw) THEN
339 118 : CALL get_qs_env(qs_env, natom=natom)
340 118 : CALL init_rho0(local_rho_set, qs_env, dft_control%qs_control%gapw_control)
341 118 : CALL rho0_s_grid_create(pw_env, local_rho_set%rho0_mpole)
342 118 : CALL hartree_local_create(hartree_local)
343 118 : CALL init_coulomb_local(hartree_local, natom)
344 : END IF
345 :
346 174 : CALL get_qs_env(qs_env=qs_env, oce=oce, sab_orb=sab)
347 174 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
348 : CALL calculate_rho_atom_coeff(qs_env, rho_ao_kp, local_rho_set%rho_atom_set, &
349 174 : qs_kind_set, oce, sab, para_env)
350 174 : CALL prepare_gapw_den(qs_env, local_rho_set, do_rho0=gapw)
351 :
352 174 : IF (dft_control%do_admm) THEN
353 28 : CALL get_qs_env(qs_env, admm_env=admm_env)
354 28 : xc_section => admm_env%xc_section_primary
355 : ELSE
356 146 : xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC")
357 : END IF
358 : CALL calculate_vxc_atom(qs_env, .FALSE., exc1=exc1, xc_section_external=xc_section, &
359 174 : rho_atom_set_external=local_rho_set%rho_atom_set)
360 :
361 174 : IF (gapw) THEN
362 118 : CALL Vh_1c_gg_integrals(qs_env, eh1c, hartree_local%ecoul_1c, local_rho_set, para_env, .FALSE.)
363 : CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, para_env, calculate_forces=.FALSE., &
364 118 : local_rho_set=local_rho_set)
365 : END IF
366 :
367 174 : IF (dft_control%do_admm) THEN
368 28 : IF (admm_env%do_gapw) THEN
369 28 : CALL local_rho_set_create(local_rho_set_admm)
370 : CALL allocate_rho_atom_internals(local_rho_set_admm%rho_atom_set, atomic_kind_set, &
371 28 : admm_env%admm_gapw_env%admm_kind_set, dft_control, para_env)
372 28 : oce => admm_env%admm_gapw_env%oce
373 28 : sab => admm_env%sab_aux_fit
374 28 : CALL get_admm_env(admm_env, rho_aux_fit=rho_aux_fit)
375 28 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_aux)
376 : CALL calculate_rho_atom_coeff(qs_env, rho_ao_aux, local_rho_set_admm%rho_atom_set, &
377 28 : admm_env%admm_gapw_env%admm_kind_set, oce, sab, para_env)
378 : CALL prepare_gapw_den(qs_env, local_rho_set=local_rho_set_admm, &
379 28 : do_rho0=.FALSE., kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
380 : !compute the potential due to atomic densities
381 28 : xc_section => admm_env%xc_section_aux
382 : CALL calculate_vxc_atom(qs_env, energy_only=.FALSE., exc1=exc1_admm, &
383 : kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
384 : xc_section_external=xc_section, &
385 28 : rho_atom_set_external=local_rho_set_admm%rho_atom_set)
386 : END IF
387 : END IF
388 :
389 : ! clean up
390 174 : CALL hartree_local_release(hartree_local)
391 :
392 : ELSE
393 :
394 1016 : NULLIFY (local_rho_set, local_rho_set_admm)
395 :
396 : END IF
397 :
398 1190 : CALL timestop(handle)
399 :
400 1190 : END SUBROUTINE ks_ref_potential_atom
401 :
402 : ! **************************************************************************************************
403 :
404 : END MODULE qs_ks_reference
|
