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 Utility subroutine for qs energy calculation
10 : !> \par History
11 : !> none
12 : !> \author MK (29.10.2002)
13 : ! **************************************************************************************************
14 : MODULE qs_energy_utils
15 : USE atomic_kind_types, ONLY: atomic_kind_type
16 : USE atprop_types, ONLY: atprop_array_add,&
17 : atprop_array_init,&
18 : atprop_type
19 : USE cp_control_types, ONLY: dft_control_type
20 : USE cp_control_utils, ONLY: read_ddapc_section
21 : USE dbcsr_api, ONLY: dbcsr_copy,&
22 : dbcsr_create,&
23 : dbcsr_p_type,&
24 : dbcsr_release,&
25 : dbcsr_set
26 : USE et_coupling, ONLY: calc_et_coupling
27 : USE et_coupling_proj, ONLY: calc_et_coupling_proj
28 : USE hartree_local_methods, ONLY: Vh_1c_gg_integrals
29 : USE hartree_local_types, ONLY: ecoul_1center_type
30 : USE input_section_types, ONLY: section_vals_get,&
31 : section_vals_get_subs_vals,&
32 : section_vals_type
33 : USE kinds, ONLY: dp
34 : USE message_passing, ONLY: mp_para_env_type
35 : USE mulliken, ONLY: atom_trace
36 : USE post_scf_bandstructure_methods, ONLY: post_scf_bandstructure
37 : USE pw_env_types, ONLY: pw_env_get,&
38 : pw_env_type
39 : USE pw_methods, ONLY: pw_axpy,&
40 : pw_scale
41 : USE pw_pool_types, ONLY: pw_pool_type
42 : USE pw_types, ONLY: pw_r3d_rs_type
43 : USE qs_dispersion_types, ONLY: qs_dispersion_type
44 : USE qs_energy_types, ONLY: qs_energy_type
45 : USE qs_environment_types, ONLY: get_qs_env,&
46 : qs_environment_type
47 : USE qs_integrate_potential, ONLY: integrate_v_core_rspace,&
48 : integrate_v_rspace
49 : USE qs_kind_types, ONLY: qs_kind_type
50 : USE qs_ks_atom, ONLY: update_ks_atom
51 : USE qs_ks_methods, ONLY: qs_ks_update_qs_env
52 : USE qs_ks_types, ONLY: qs_ks_env_type
53 : USE qs_linres_module, ONLY: linres_calculation_low
54 : USE qs_local_rho_types, ONLY: local_rho_type
55 : USE qs_rho0_ggrid, ONLY: integrate_vhg0_rspace
56 : USE qs_rho_atom_types, ONLY: rho_atom_type,&
57 : zero_rho_atom_integrals
58 : USE qs_rho_types, ONLY: qs_rho_get,&
59 : qs_rho_type
60 : USE qs_scf, ONLY: scf
61 : USE qs_tddfpt2_methods, ONLY: tddfpt
62 : USE qs_vxc, ONLY: qs_xc_density
63 : USE qs_vxc_atom, ONLY: calculate_vxc_atom
64 : USE tip_scan_methods, ONLY: tip_scanning
65 : USE xas_methods, ONLY: xas
66 : USE xas_tdp_methods, ONLY: xas_tdp
67 : USE xc_derivatives, ONLY: xc_functionals_get_needs
68 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
69 : #include "./base/base_uses.f90"
70 :
71 : IMPLICIT NONE
72 :
73 : PRIVATE
74 :
75 : ! *** Global parameters ***
76 :
77 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_energy_utils'
78 :
79 : PUBLIC :: qs_energies_properties
80 :
81 : CONTAINS
82 :
83 : ! **************************************************************************************************
84 : !> \brief Refactoring of qs_energies_scf. Moves computation of properties
85 : !> into separate subroutine
86 : !> \param qs_env ...
87 : !> \param calc_forces ...
88 : !> \par History
89 : !> 05.2013 created [Florian Schiffmann]
90 : ! **************************************************************************************************
91 :
92 75372 : SUBROUTINE qs_energies_properties(qs_env, calc_forces)
93 : TYPE(qs_environment_type), POINTER :: qs_env
94 : LOGICAL, INTENT(IN) :: calc_forces
95 :
96 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_energies_properties'
97 :
98 : INTEGER :: handle, natom
99 : LOGICAL :: do_et, do_et_proj, &
100 : do_post_scf_bandstructure, do_tip_scan
101 : REAL(KIND=dp) :: ekts
102 : TYPE(atprop_type), POINTER :: atprop
103 : TYPE(dft_control_type), POINTER :: dft_control
104 : TYPE(mp_para_env_type), POINTER :: para_env
105 : TYPE(pw_env_type), POINTER :: pw_env
106 : TYPE(pw_r3d_rs_type), POINTER :: v_hartree_rspace
107 : TYPE(qs_energy_type), POINTER :: energy
108 : TYPE(section_vals_type), POINTER :: input, post_scf_bands_section, &
109 : proj_section, rest_b_section, &
110 : tip_scan_section
111 :
112 18843 : NULLIFY (atprop, energy, pw_env)
113 18843 : CALL timeset(routineN, handle)
114 :
115 : ! atomic energies using Mulliken partition
116 : CALL get_qs_env(qs_env, &
117 : dft_control=dft_control, &
118 : input=input, &
119 : atprop=atprop, &
120 : energy=energy, &
121 : v_hartree_rspace=v_hartree_rspace, &
122 : para_env=para_env, &
123 18843 : pw_env=pw_env)
124 18843 : IF (atprop%energy) THEN
125 140 : CALL qs_energies_mulliken(qs_env)
126 140 : CALL get_qs_env(qs_env, natom=natom)
127 : IF (.NOT. dft_control%qs_control%semi_empirical .AND. &
128 140 : .NOT. dft_control%qs_control%xtb .AND. &
129 : .NOT. dft_control%qs_control%dftb) THEN
130 : ! Nuclear charge correction
131 36 : CALL integrate_v_core_rspace(v_hartree_rspace, qs_env)
132 36 : IF (.NOT. ASSOCIATED(atprop%ateb)) THEN
133 8 : CALL atprop_array_init(atprop%ateb, natom)
134 : END IF
135 : ! Kohn-Sham Functional corrections
136 36 : CALL ks_xc_correction(qs_env)
137 : END IF
138 140 : CALL atprop_array_add(atprop%atener, atprop%ateb)
139 140 : CALL atprop_array_add(atprop%atener, atprop%ateself)
140 140 : CALL atprop_array_add(atprop%atener, atprop%atexc)
141 140 : CALL atprop_array_add(atprop%atener, atprop%atecoul)
142 140 : CALL atprop_array_add(atprop%atener, atprop%atevdw)
143 140 : CALL atprop_array_add(atprop%atener, atprop%ategcp)
144 140 : CALL atprop_array_add(atprop%atener, atprop%atecc)
145 140 : CALL atprop_array_add(atprop%atener, atprop%ate1c)
146 : ! entropic energy
147 140 : ekts = energy%kts/REAL(natom, KIND=dp)/REAL(para_env%num_pe, KIND=dp)
148 7504 : atprop%atener(:) = atprop%atener(:) + ekts
149 : END IF
150 :
151 : ! ET coupling - projection-operator approach
152 18843 : NULLIFY (proj_section)
153 : proj_section => &
154 18843 : section_vals_get_subs_vals(input, "PROPERTIES%ET_COUPLING%PROJECTION")
155 18843 : CALL section_vals_get(proj_section, explicit=do_et_proj)
156 18843 : IF (do_et_proj) THEN
157 10 : CALL calc_et_coupling_proj(qs_env)
158 : END IF
159 :
160 : ! ********** Calculate the electron transfer coupling elements********
161 18843 : do_et = .FALSE.
162 18843 : do_et = dft_control%qs_control%et_coupling_calc
163 18843 : IF (do_et) THEN
164 0 : qs_env%et_coupling%energy = energy%total
165 0 : qs_env%et_coupling%keep_matrix = .TRUE.
166 0 : qs_env%et_coupling%first_run = .TRUE.
167 0 : CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., just_energy=.TRUE.)
168 0 : qs_env%et_coupling%first_run = .FALSE.
169 0 : IF (dft_control%qs_control%ddapc_restraint) THEN
170 0 : rest_b_section => section_vals_get_subs_vals(input, "PROPERTIES%ET_COUPLING%DDAPC_RESTRAINT_B")
171 : CALL read_ddapc_section(qs_control=dft_control%qs_control, &
172 0 : ddapc_restraint_section=rest_b_section)
173 : END IF
174 0 : CALL scf(qs_env=qs_env)
175 0 : qs_env%et_coupling%keep_matrix = .TRUE.
176 :
177 0 : CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., just_energy=.TRUE.)
178 0 : CALL calc_et_coupling(qs_env)
179 : END IF
180 :
181 : !Properties
182 18843 : IF (dft_control%do_xas_calculation) THEN
183 42 : CALL xas(qs_env, dft_control)
184 : END IF
185 :
186 18843 : IF (dft_control%do_xas_tdp_calculation) THEN
187 50 : CALL xas_tdp(qs_env)
188 : END IF
189 :
190 : ! Compute Linear Response properties as post-scf
191 18843 : IF (.NOT. qs_env%linres_run) THEN
192 18343 : CALL linres_calculation_low(qs_env)
193 : END IF
194 :
195 18843 : IF (dft_control%tddfpt2_control%enabled) THEN
196 1054 : CALL tddfpt(qs_env, calc_forces)
197 : END IF
198 :
199 : ! post-SCF bandstructure calculation from higher level methods
200 18843 : NULLIFY (post_scf_bands_section)
201 18843 : post_scf_bands_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%BANDSTRUCTURE")
202 18843 : CALL section_vals_get(post_scf_bands_section, explicit=do_post_scf_bandstructure)
203 18843 : IF (do_post_scf_bandstructure) THEN
204 18 : CALL post_scf_bandstructure(qs_env, post_scf_bands_section)
205 : END IF
206 :
207 : ! tip scan
208 18843 : NULLIFY (tip_scan_section)
209 18843 : tip_scan_section => section_vals_get_subs_vals(input, "PROPERTIES%TIP_SCAN")
210 18843 : CALL section_vals_get(tip_scan_section, explicit=do_tip_scan)
211 18843 : IF (do_tip_scan) THEN
212 0 : CALL tip_scanning(qs_env, tip_scan_section)
213 : END IF
214 :
215 18843 : CALL timestop(handle)
216 :
217 18843 : END SUBROUTINE qs_energies_properties
218 :
219 : ! **************************************************************************************************
220 : !> \brief Use a simple Mulliken-like energy decomposition
221 : !> \param qs_env ...
222 : !> \date 07.2011
223 : !> \author JHU
224 : !> \version 1.0
225 : ! **************************************************************************************************
226 140 : SUBROUTINE qs_energies_mulliken(qs_env)
227 :
228 : TYPE(qs_environment_type), POINTER :: qs_env
229 :
230 : INTEGER :: ispin
231 : TYPE(atprop_type), POINTER :: atprop
232 140 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_h, matrix_ks, rho_ao
233 : TYPE(qs_rho_type), POINTER :: rho
234 :
235 140 : NULLIFY (atprop, matrix_h, matrix_ks, rho, rho_ao)
236 : CALL get_qs_env(qs_env=qs_env, matrix_ks=matrix_ks, matrix_h=matrix_h, &
237 140 : rho=rho, atprop=atprop)
238 140 : CALL qs_rho_get(rho, rho_ao=rho_ao)
239 :
240 140 : IF (atprop%energy) THEN
241 : ! E = 0.5*Tr(H*P+F*P)
242 7504 : atprop%atener = 0._dp
243 284 : DO ispin = 1, SIZE(rho_ao)
244 : CALL atom_trace(matrix_h(1)%matrix, rho_ao(ispin)%matrix, &
245 144 : 0.5_dp, atprop%atener)
246 : CALL atom_trace(matrix_ks(ispin)%matrix, rho_ao(ispin)%matrix, &
247 284 : 0.5_dp, atprop%atener)
248 : END DO
249 : END IF
250 :
251 140 : END SUBROUTINE qs_energies_mulliken
252 :
253 : ! **************************************************************************************************
254 : !> \brief ...
255 : !> \param qs_env ...
256 : ! **************************************************************************************************
257 36 : SUBROUTINE ks_xc_correction(qs_env)
258 : TYPE(qs_environment_type), POINTER :: qs_env
259 :
260 : CHARACTER(len=*), PARAMETER :: routineN = 'ks_xc_correction'
261 :
262 : INTEGER :: handle, iatom, ispin, natom, nspins
263 : LOGICAL :: gapw, gapw_xc
264 : REAL(KIND=dp) :: eh1, exc1
265 36 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
266 : TYPE(atprop_type), POINTER :: atprop
267 36 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, rho_ao, xcmat
268 36 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p
269 : TYPE(dft_control_type), POINTER :: dft_control
270 36 : TYPE(ecoul_1center_type), DIMENSION(:), POINTER :: ecoul_1c
271 : TYPE(local_rho_type), POINTER :: local_rho_set
272 : TYPE(mp_para_env_type), POINTER :: para_env
273 : TYPE(pw_env_type), POINTER :: pw_env
274 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
275 : TYPE(pw_r3d_rs_type) :: xc_den
276 36 : TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: vtau, vxc
277 : TYPE(pw_r3d_rs_type), POINTER :: v_hartree_rspace
278 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
279 36 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
280 : TYPE(qs_ks_env_type), POINTER :: ks_env
281 : TYPE(qs_rho_type), POINTER :: rho_struct
282 36 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
283 : TYPE(section_vals_type), POINTER :: xc_fun_section, xc_section
284 : TYPE(xc_rho_cflags_type) :: needs
285 :
286 36 : CALL timeset(routineN, handle)
287 :
288 36 : CALL get_qs_env(qs_env, ks_env=ks_env, dft_control=dft_control, pw_env=pw_env, atprop=atprop)
289 :
290 36 : IF (atprop%energy) THEN
291 :
292 36 : nspins = dft_control%nspins
293 36 : xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC")
294 36 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
295 36 : needs = xc_functionals_get_needs(xc_fun_section, (nspins == 2), .TRUE.)
296 36 : gapw = dft_control%qs_control%gapw
297 36 : gapw_xc = dft_control%qs_control%gapw_xc
298 :
299 : ! Nuclear charge correction
300 36 : CALL get_qs_env(qs_env, v_hartree_rspace=v_hartree_rspace)
301 36 : IF (gapw .OR. gapw_xc) THEN
302 : CALL get_qs_env(qs_env=qs_env, local_rho_set=local_rho_set, &
303 : rho_atom_set=rho_atom_set, ecoul_1c=ecoul_1c, &
304 12 : natom=natom, para_env=para_env)
305 12 : CALL zero_rho_atom_integrals(rho_atom_set)
306 12 : CALL calculate_vxc_atom(qs_env, .FALSE., exc1)
307 12 : IF (gapw) THEN
308 8 : CALL Vh_1c_gg_integrals(qs_env, eh1, ecoul_1c, local_rho_set, para_env, tddft=.FALSE.)
309 8 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
310 : CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, para_env, calculate_forces=.FALSE., &
311 8 : local_rho_set=local_rho_set, atener=atprop%ateb)
312 : END IF
313 : END IF
314 :
315 36 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
316 36 : CALL auxbas_pw_pool%create_pw(xc_den)
317 148 : ALLOCATE (vxc(nspins))
318 76 : DO ispin = 1, nspins
319 76 : CALL auxbas_pw_pool%create_pw(vxc(ispin))
320 : END DO
321 36 : IF (needs%tau .OR. needs%tau_spin) THEN
322 48 : ALLOCATE (vtau(nspins))
323 24 : DO ispin = 1, nspins
324 48 : CALL auxbas_pw_pool%create_pw(vtau(ispin))
325 : END DO
326 : END IF
327 :
328 36 : IF (gapw_xc) THEN
329 4 : CALL get_qs_env(qs_env, rho_xc=rho_struct, dispersion_env=dispersion_env)
330 : ELSE
331 32 : CALL get_qs_env(qs_env, rho=rho_struct, dispersion_env=dispersion_env)
332 : END IF
333 36 : IF (needs%tau .OR. needs%tau_spin) THEN
334 : CALL qs_xc_density(ks_env, rho_struct, xc_section, dispersion_env=dispersion_env, &
335 12 : xc_den=xc_den, vxc=vxc, vtau=vtau)
336 : ELSE
337 : CALL qs_xc_density(ks_env, rho_struct, xc_section, dispersion_env=dispersion_env, &
338 24 : xc_den=xc_den, vxc=vxc)
339 : END IF
340 36 : CALL get_qs_env(qs_env, rho=rho_struct)
341 36 : CALL qs_rho_get(rho_struct, rho_ao=rho_ao)
342 36 : CALL get_qs_env(qs_env, natom=natom, matrix_s=matrix_s)
343 36 : CALL atprop_array_init(atprop%atexc, natom)
344 148 : ALLOCATE (xcmat(nspins))
345 76 : DO ispin = 1, nspins
346 40 : ALLOCATE (xcmat(ispin)%matrix)
347 40 : CALL dbcsr_create(xcmat(ispin)%matrix, template=matrix_s(1)%matrix)
348 40 : CALL dbcsr_copy(xcmat(ispin)%matrix, matrix_s(1)%matrix)
349 40 : CALL dbcsr_set(xcmat(ispin)%matrix, 0.0_dp)
350 40 : CALL pw_scale(vxc(ispin), -0.5_dp)
351 40 : CALL pw_axpy(xc_den, vxc(ispin))
352 40 : CALL pw_scale(vxc(ispin), vxc(ispin)%pw_grid%dvol)
353 : CALL integrate_v_rspace(qs_env=qs_env, v_rspace=vxc(ispin), hmat=xcmat(ispin), &
354 40 : calculate_forces=.FALSE., gapw=(gapw .OR. gapw_xc))
355 76 : IF (needs%tau .OR. needs%tau_spin) THEN
356 12 : CALL pw_scale(vtau(ispin), -0.5_dp*vtau(ispin)%pw_grid%dvol)
357 : CALL integrate_v_rspace(qs_env=qs_env, v_rspace=vtau(ispin), &
358 : hmat=xcmat(ispin), calculate_forces=.FALSE., &
359 12 : gapw=(gapw .OR. gapw_xc), compute_tau=.TRUE.)
360 : END IF
361 : END DO
362 36 : IF (gapw .OR. gapw_xc) THEN
363 : ! remove one-center potential matrix part
364 12 : CALL qs_rho_get(rho_struct, rho_ao_kp=matrix_p)
365 12 : CALL update_ks_atom(qs_env, xcmat, matrix_p, forces=.FALSE., kscale=-0.5_dp)
366 12 : CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom_set)
367 12 : CALL atprop_array_init(atprop%ate1c, natom)
368 48 : atprop%ate1c = 0.0_dp
369 48 : DO iatom = 1, natom
370 : atprop%ate1c(iatom) = atprop%ate1c(iatom) + &
371 48 : rho_atom_set(iatom)%exc_h - rho_atom_set(iatom)%exc_s
372 : END DO
373 12 : IF (gapw) THEN
374 8 : CALL get_qs_env(qs_env=qs_env, ecoul_1c=ecoul_1c)
375 32 : DO iatom = 1, natom
376 : atprop%ate1c(iatom) = atprop%ate1c(iatom) + &
377 : ecoul_1c(iatom)%ecoul_1_h - ecoul_1c(iatom)%ecoul_1_s + &
378 32 : ecoul_1c(iatom)%ecoul_1_z - ecoul_1c(iatom)%ecoul_1_0
379 : END DO
380 : END IF
381 : END IF
382 76 : DO ispin = 1, nspins
383 40 : CALL atom_trace(xcmat(ispin)%matrix, rho_ao(ispin)%matrix, 1.0_dp, atprop%atexc)
384 40 : CALL dbcsr_release(xcmat(ispin)%matrix)
385 76 : DEALLOCATE (xcmat(ispin)%matrix)
386 : END DO
387 36 : DEALLOCATE (xcmat)
388 :
389 36 : CALL auxbas_pw_pool%give_back_pw(xc_den)
390 76 : DO ispin = 1, nspins
391 76 : CALL auxbas_pw_pool%give_back_pw(vxc(ispin))
392 : END DO
393 36 : IF (needs%tau .OR. needs%tau_spin) THEN
394 24 : DO ispin = 1, nspins
395 48 : CALL auxbas_pw_pool%give_back_pw(vtau(ispin))
396 : END DO
397 : END IF
398 :
399 : END IF
400 :
401 36 : CALL timestop(handle)
402 :
403 72 : END SUBROUTINE ks_xc_correction
404 :
405 : END MODULE qs_energy_utils
|
