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 Utilities for broadened DOS and PDOS output.
10 : ! **************************************************************************************************
11 : MODULE qs_dos_utils
12 : USE cp_units, ONLY: cp_unit_from_cp2k
13 : USE input_section_types, ONLY: section_vals_get,&
14 : section_vals_get_subs_vals,&
15 : section_vals_type
16 : USE kinds, ONLY: dp
17 : USE mathconstants, ONLY: pi
18 : #include "./base/base_uses.f90"
19 :
20 : IMPLICIT NONE
21 :
22 : PRIVATE
23 :
24 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dos_utils'
25 :
26 : INTEGER, PARAMETER, PUBLIC :: broadening_gaussian = 1, &
27 : broadening_lorentzian = 2, &
28 : broadening_pseudo_voigt = 3
29 : INTEGER, PARAMETER, PUBLIC :: dos_energy_unit_hartree = 1, &
30 : dos_energy_unit_ev = 2
31 : INTEGER, PARAMETER, PUBLIC :: dos_energy_zero_auto = 1, &
32 : dos_energy_zero_absolute = 2, &
33 : dos_energy_zero_fermi = 3, &
34 : dos_energy_zero_hoco = 4
35 :
36 : PUBLIC :: add_broadened_peak, add_broadened_value, broadening_cutoff, broadening_function, &
37 : dos_density_scale, dos_energy_label, dos_energy_scale, dos_energy_zero_label, &
38 : dos_resolve_energy_zero, get_dos_pdos_flags, write_broadening_info
39 :
40 : CONTAINS
41 :
42 : ! **************************************************************************************************
43 : !> \brief Return the conversion factor from internal energy units to the selected DOS energy unit.
44 : !> \param energy_unit ...
45 : !> \return ...
46 : ! **************************************************************************************************
47 326 : FUNCTION dos_energy_scale(energy_unit) RESULT(scale)
48 :
49 : INTEGER, INTENT(IN) :: energy_unit
50 : REAL(KIND=dp) :: scale
51 :
52 446 : SELECT CASE (energy_unit)
53 : CASE (dos_energy_unit_ev)
54 120 : scale = cp_unit_from_cp2k(value=1.0_dp, unit_str="eV")
55 : CASE DEFAULT
56 326 : scale = 1.0_dp
57 : END SELECT
58 :
59 326 : END FUNCTION dos_energy_scale
60 :
61 : ! **************************************************************************************************
62 : !> \brief Return the DOS-density conversion factor for the selected energy unit.
63 : !> \param energy_unit ...
64 : !> \return ...
65 : ! **************************************************************************************************
66 86 : FUNCTION dos_density_scale(energy_unit) RESULT(scale)
67 :
68 : INTEGER, INTENT(IN) :: energy_unit
69 : REAL(KIND=dp) :: scale
70 :
71 86 : scale = 1.0_dp/dos_energy_scale(energy_unit)
72 :
73 86 : END FUNCTION dos_density_scale
74 :
75 : ! **************************************************************************************************
76 : !> \brief Return the energy-column label for DOS-like output.
77 : !> \param energy_unit ...
78 : !> \return ...
79 : ! **************************************************************************************************
80 120 : FUNCTION dos_energy_label(energy_unit) RESULT(label)
81 :
82 : INTEGER, INTENT(IN) :: energy_unit
83 : CHARACTER(LEN=16) :: label
84 :
85 120 : SELECT CASE (energy_unit)
86 : CASE (dos_energy_unit_ev)
87 0 : label = "Energy[eV]"
88 : CASE DEFAULT
89 120 : label = "Energy[a.u.]"
90 : END SELECT
91 :
92 120 : END FUNCTION dos_energy_label
93 :
94 : ! **************************************************************************************************
95 : !> \brief Return the label for the selected DOS energy zero.
96 : !> \param energy_zero ...
97 : !> \return ...
98 : ! **************************************************************************************************
99 120 : FUNCTION dos_energy_zero_label(energy_zero) RESULT(label)
100 :
101 : INTEGER, INTENT(IN) :: energy_zero
102 : CHARACTER(LEN=16) :: label
103 :
104 120 : SELECT CASE (energy_zero)
105 : CASE (dos_energy_zero_absolute)
106 0 : label = "ABSOLUTE"
107 : CASE (dos_energy_zero_hoco)
108 98 : label = "HOCO"
109 : CASE (dos_energy_zero_auto)
110 0 : label = "AUTO"
111 : CASE DEFAULT
112 120 : label = "FERMI"
113 : END SELECT
114 :
115 120 : END FUNCTION dos_energy_zero_label
116 :
117 : ! **************************************************************************************************
118 : !> \brief Resolve AUTO energy-zero selection for DOS-like output.
119 : !> \param energy_zero ...
120 : !> \param smearing_enabled ...
121 : !> \param fractional_occupation ...
122 : !> \return ...
123 : ! **************************************************************************************************
124 120 : FUNCTION dos_resolve_energy_zero(energy_zero, smearing_enabled, fractional_occupation) RESULT(resolved)
125 :
126 : INTEGER, INTENT(IN) :: energy_zero
127 : LOGICAL, INTENT(IN) :: smearing_enabled, fractional_occupation
128 : INTEGER :: resolved
129 :
130 120 : IF (energy_zero == dos_energy_zero_auto) THEN
131 120 : IF (smearing_enabled .OR. fractional_occupation) THEN
132 : resolved = dos_energy_zero_fermi
133 : ELSE
134 98 : resolved = dos_energy_zero_hoco
135 : END IF
136 : ELSE
137 : resolved = energy_zero
138 : END IF
139 :
140 120 : END FUNCTION dos_resolve_energy_zero
141 :
142 : ! **************************************************************************************************
143 : !> \brief Resolve projected-DOS requests from a DOS print section.
144 : !> \param dos_section DOS print section
145 : !> \param do_dos_output whether the DOS print key is active
146 : !> \param do_projected_dos whether any projected DOS output is requested
147 : !> \param do_pdos whether kind-resolved PDOS output is requested
148 : !> \param do_curve whether broadened curve output is requested
149 : ! **************************************************************************************************
150 22953 : SUBROUTINE get_dos_pdos_flags(dos_section, do_dos_output, do_projected_dos, do_pdos, do_curve)
151 :
152 : TYPE(section_vals_type), POINTER :: dos_section
153 : LOGICAL, INTENT(IN) :: do_dos_output
154 : LOGICAL, INTENT(OUT) :: do_projected_dos, do_pdos, do_curve
155 :
156 : INTEGER :: nrep
157 : LOGICAL :: has_ldos, has_r_ldos
158 : TYPE(section_vals_type), POINTER :: curve_section, ldos_section, pdos_section
159 :
160 22953 : do_projected_dos = .FALSE.
161 22953 : do_pdos = .FALSE.
162 22953 : do_curve = .FALSE.
163 22953 : has_ldos = .FALSE.
164 22953 : has_r_ldos = .FALSE.
165 :
166 22953 : IF (do_dos_output) THEN
167 86 : pdos_section => section_vals_get_subs_vals(dos_section, "PDOS")
168 86 : CALL section_vals_get(pdos_section, explicit=do_pdos)
169 86 : curve_section => section_vals_get_subs_vals(dos_section, "CURVE")
170 86 : CALL section_vals_get(curve_section, explicit=do_curve)
171 86 : ldos_section => section_vals_get_subs_vals(dos_section, "LDOS")
172 86 : CALL section_vals_get(ldos_section, n_repetition=nrep)
173 86 : has_ldos = (nrep > 0)
174 86 : ldos_section => section_vals_get_subs_vals(dos_section, "R_LDOS")
175 86 : CALL section_vals_get(ldos_section, n_repetition=nrep)
176 86 : has_r_ldos = (nrep > 0)
177 : END IF
178 :
179 22953 : do_projected_dos = do_pdos .OR. has_ldos .OR. has_r_ldos
180 :
181 22953 : END SUBROUTINE get_dos_pdos_flags
182 :
183 : ! **************************************************************************************************
184 : !> \brief Add a broadened spectral line to a DOS curve.
185 : !> \param dos ...
186 : !> \param occ_dos ...
187 : !> \param emin ...
188 : !> \param de ...
189 : !> \param eig ...
190 : !> \param occ ...
191 : !> \param weight ...
192 : !> \param broaden_type ...
193 : !> \param broaden_width ...
194 : !> \param voigt_mixing ...
195 : ! **************************************************************************************************
196 0 : SUBROUTINE add_broadened_peak(dos, occ_dos, emin, de, eig, occ, weight, broaden_type, &
197 : broaden_width, voigt_mixing)
198 :
199 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: dos, occ_dos
200 : REAL(KIND=dp), INTENT(IN) :: emin, de, eig, occ, weight
201 : INTEGER, INTENT(IN) :: broaden_type
202 : REAL(KIND=dp), INTENT(IN) :: broaden_width, voigt_mixing
203 :
204 : INTEGER :: i, iend, istart, nhist
205 : REAL(KIND=dp) :: eval, line_shape
206 :
207 0 : nhist = SIZE(dos)
208 0 : istart = MAX(1, FLOOR((eig - broadening_cutoff(broaden_type, broaden_width) - emin)/de) + 1)
209 0 : iend = MIN(nhist, CEILING((eig + broadening_cutoff(broaden_type, broaden_width) - emin)/de) + 1)
210 0 : DO i = istart, iend
211 0 : eval = emin + (i - 1)*de
212 0 : line_shape = broadening_function(eval - eig, broaden_type, broaden_width, voigt_mixing)
213 0 : dos(i) = dos(i) + weight*line_shape
214 0 : occ_dos(i) = occ_dos(i) + weight*occ*line_shape
215 : END DO
216 :
217 0 : END SUBROUTINE add_broadened_peak
218 :
219 : ! **************************************************************************************************
220 : !> \brief Add a broadened spectral line with a scalar weight to a curve.
221 : !> \param curve ...
222 : !> \param emin ...
223 : !> \param de ...
224 : !> \param eig ...
225 : !> \param weight ...
226 : !> \param broaden_type ...
227 : !> \param broaden_width ...
228 : !> \param voigt_mixing ...
229 : ! **************************************************************************************************
230 0 : SUBROUTINE add_broadened_value(curve, emin, de, eig, weight, broaden_type, &
231 : broaden_width, voigt_mixing)
232 :
233 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: curve
234 : REAL(KIND=dp), INTENT(IN) :: emin, de, eig, weight
235 : INTEGER, INTENT(IN) :: broaden_type
236 : REAL(KIND=dp), INTENT(IN) :: broaden_width, voigt_mixing
237 :
238 : INTEGER :: i, iend, istart, nhist
239 : REAL(KIND=dp) :: eval, line_shape
240 :
241 0 : nhist = SIZE(curve)
242 0 : istart = MAX(1, FLOOR((eig - broadening_cutoff(broaden_type, broaden_width) - emin)/de) + 1)
243 0 : iend = MIN(nhist, CEILING((eig + broadening_cutoff(broaden_type, broaden_width) - emin)/de) + 1)
244 0 : DO i = istart, iend
245 0 : eval = emin + (i - 1)*de
246 0 : line_shape = broadening_function(eval - eig, broaden_type, broaden_width, voigt_mixing)
247 0 : curve(i) = curve(i) + weight*line_shape
248 : END DO
249 :
250 0 : END SUBROUTINE add_broadened_value
251 :
252 : ! **************************************************************************************************
253 : !> \brief Broadening cutoff used for numerical accumulation.
254 : !> \param broaden_type ...
255 : !> \param broaden_width ...
256 : !> \return ...
257 : ! **************************************************************************************************
258 0 : PURE FUNCTION broadening_cutoff(broaden_type, broaden_width) RESULT(cutoff)
259 :
260 : INTEGER, INTENT(IN) :: broaden_type
261 : REAL(KIND=dp), INTENT(IN) :: broaden_width
262 : REAL(KIND=dp) :: cutoff
263 :
264 0 : IF (broaden_type == broadening_gaussian) THEN
265 0 : cutoff = 8.0_dp*broaden_width
266 : ELSE
267 0 : cutoff = 50.0_dp*broaden_width
268 : END IF
269 :
270 0 : END FUNCTION broadening_cutoff
271 :
272 : ! **************************************************************************************************
273 : !> \brief Normalized broadening function. BROADEN_WIDTH is FWHM.
274 : !> \param delta_e ...
275 : !> \param broaden_type ...
276 : !> \param broaden_width ...
277 : !> \param voigt_mixing ...
278 : !> \return ...
279 : ! **************************************************************************************************
280 0 : PURE FUNCTION broadening_function(delta_e, broaden_type, broaden_width, voigt_mixing) RESULT(value)
281 :
282 : REAL(KIND=dp), INTENT(IN) :: delta_e
283 : INTEGER, INTENT(IN) :: broaden_type
284 : REAL(KIND=dp), INTENT(IN) :: broaden_width, voigt_mixing
285 : REAL(KIND=dp) :: value
286 :
287 : REAL(KIND=dp) :: eta, gamma, gaussian_value, sigma
288 :
289 0 : IF (broaden_width <= 0.0_dp) THEN
290 0 : value = 0.0_dp
291 : RETURN
292 : END IF
293 :
294 0 : sigma = broaden_width/(2.0_dp*SQRT(2.0_dp*LOG(2.0_dp)))
295 0 : gamma = 0.5_dp*broaden_width
296 0 : gaussian_value = EXP(-0.5_dp*(delta_e/sigma)**2)/(sigma*SQRT(2.0_dp*pi))
297 :
298 0 : SELECT CASE (broaden_type)
299 : CASE (broadening_gaussian)
300 0 : value = gaussian_value
301 : CASE (broadening_lorentzian)
302 0 : value = gamma/(pi*(delta_e**2 + gamma**2))
303 : CASE (broadening_pseudo_voigt)
304 0 : eta = MIN(1.0_dp, MAX(0.0_dp, voigt_mixing))
305 0 : value = eta*gamma/(pi*(delta_e**2 + gamma**2)) + (1.0_dp - eta)*gaussian_value
306 : CASE DEFAULT
307 0 : value = gaussian_value
308 : END SELECT
309 :
310 : END FUNCTION broadening_function
311 :
312 : ! **************************************************************************************************
313 : !> \brief Write broadening metadata.
314 : !> \param iw ...
315 : !> \param broaden_type ...
316 : !> \param broaden_width ...
317 : !> \param voigt_mixing ...
318 : ! **************************************************************************************************
319 0 : SUBROUTINE write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing)
320 :
321 : INTEGER, INTENT(IN) :: iw, broaden_type
322 : REAL(KIND=dp), INTENT(IN) :: broaden_width, voigt_mixing
323 :
324 : REAL(KIND=dp) :: broaden_width_ev
325 :
326 0 : broaden_width_ev = cp_unit_from_cp2k(value=broaden_width, unit_str="eV")
327 :
328 0 : SELECT CASE (broaden_type)
329 : CASE (broadening_gaussian)
330 : WRITE (UNIT=iw, FMT="(A,F10.8,A,F8.6,A)") &
331 0 : "# Gaussian broadening, FWHM = ", broaden_width, " a.u. = ", broaden_width_ev, " eV"
332 : CASE (broadening_lorentzian)
333 : WRITE (UNIT=iw, FMT="(A,F10.8,A,F8.6,A)") &
334 0 : "# Lorentzian broadening, FWHM = ", broaden_width, " a.u. = ", broaden_width_ev, " eV"
335 : CASE (broadening_pseudo_voigt)
336 : WRITE (UNIT=iw, FMT="(A,F10.8,A,F8.6,A,F8.4)") &
337 0 : "# Pseudo-Voigt broadening, FWHM = ", broaden_width, " a.u. = ", &
338 0 : broaden_width_ev, " eV, Lorentzian fraction = ", MIN(1.0_dp, MAX(0.0_dp, voigt_mixing))
339 : END SELECT
340 :
341 0 : END SUBROUTINE write_broadening_info
342 :
343 : END MODULE qs_dos_utils
|
