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