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 : section_vals_val_get
17 : USE kinds, ONLY: dp
18 : USE mathconstants, ONLY: pi
19 : #include "./base/base_uses.f90"
20 :
21 : IMPLICIT NONE
22 :
23 : PRIVATE
24 :
25 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dos_utils'
26 :
27 : INTEGER, PARAMETER, PUBLIC :: broadening_gaussian = 1, &
28 : broadening_lorentzian = 2, &
29 : broadening_pseudo_voigt = 3
30 : INTEGER, PARAMETER, PUBLIC :: dos_energy_unit_hartree = 1, &
31 : dos_energy_unit_ev = 2
32 : INTEGER, PARAMETER, PUBLIC :: dos_energy_zero_auto = 1, &
33 : dos_energy_zero_absolute = 2, &
34 : dos_energy_zero_fermi = 3, &
35 : dos_energy_zero_hoco = 4
36 :
37 : PUBLIC :: add_broadened_peak, add_broadened_value, broadening_cutoff, broadening_function, &
38 : dos_density_scale, dos_energy_label, dos_energy_scale, dos_energy_zero_label, &
39 : dos_resolve_energy_zero, get_dos_pdos_flags, write_broadening_info
40 :
41 : CONTAINS
42 :
43 : ! **************************************************************************************************
44 : !> \brief Return the conversion factor from internal energy units to the selected DOS energy unit.
45 : !> \param energy_unit ...
46 : !> \return ...
47 : ! **************************************************************************************************
48 632 : FUNCTION dos_energy_scale(energy_unit) RESULT(scale)
49 :
50 : INTEGER, INTENT(IN) :: energy_unit
51 : REAL(KIND=dp) :: scale
52 :
53 854 : SELECT CASE (energy_unit)
54 : CASE (dos_energy_unit_ev)
55 222 : scale = cp_unit_from_cp2k(value=1.0_dp, unit_str="eV")
56 : CASE DEFAULT
57 632 : scale = 1.0_dp
58 : END SELECT
59 :
60 632 : END FUNCTION dos_energy_scale
61 :
62 : ! **************************************************************************************************
63 : !> \brief Return the DOS-density conversion factor for the selected energy unit.
64 : !> \param energy_unit ...
65 : !> \return ...
66 : ! **************************************************************************************************
67 188 : FUNCTION dos_density_scale(energy_unit) RESULT(scale)
68 :
69 : INTEGER, INTENT(IN) :: energy_unit
70 : REAL(KIND=dp) :: scale
71 :
72 188 : scale = 1.0_dp/dos_energy_scale(energy_unit)
73 :
74 188 : END FUNCTION dos_density_scale
75 :
76 : ! **************************************************************************************************
77 : !> \brief Return the energy-column label for DOS-like output.
78 : !> \param energy_unit ...
79 : !> \return ...
80 : ! **************************************************************************************************
81 222 : FUNCTION dos_energy_label(energy_unit) RESULT(label)
82 :
83 : INTEGER, INTENT(IN) :: energy_unit
84 : CHARACTER(LEN=16) :: label
85 :
86 222 : SELECT CASE (energy_unit)
87 : CASE (dos_energy_unit_ev)
88 0 : label = "Energy[eV]"
89 : CASE DEFAULT
90 222 : label = "Energy[a.u.]"
91 : END SELECT
92 :
93 222 : END FUNCTION dos_energy_label
94 :
95 : ! **************************************************************************************************
96 : !> \brief Return the label for the selected DOS energy zero.
97 : !> \param energy_zero ...
98 : !> \return ...
99 : ! **************************************************************************************************
100 136 : FUNCTION dos_energy_zero_label(energy_zero) RESULT(label)
101 :
102 : INTEGER, INTENT(IN) :: energy_zero
103 : CHARACTER(LEN=16) :: label
104 :
105 136 : SELECT CASE (energy_zero)
106 : CASE (dos_energy_zero_absolute)
107 0 : label = "ABSOLUTE"
108 : CASE (dos_energy_zero_hoco)
109 100 : label = "HOCO"
110 : CASE (dos_energy_zero_auto)
111 0 : label = "AUTO"
112 : CASE DEFAULT
113 136 : label = "FERMI"
114 : END SELECT
115 :
116 136 : END FUNCTION dos_energy_zero_label
117 :
118 : ! **************************************************************************************************
119 : !> \brief Resolve AUTO energy-zero selection for DOS-like output.
120 : !> \param energy_zero ...
121 : !> \param smearing_enabled ...
122 : !> \param fractional_occupation ...
123 : !> \return ...
124 : ! **************************************************************************************************
125 136 : FUNCTION dos_resolve_energy_zero(energy_zero, smearing_enabled, fractional_occupation) RESULT(resolved)
126 :
127 : INTEGER, INTENT(IN) :: energy_zero
128 : LOGICAL, INTENT(IN) :: smearing_enabled, fractional_occupation
129 : INTEGER :: resolved
130 :
131 136 : IF (energy_zero == dos_energy_zero_auto) THEN
132 136 : IF (smearing_enabled .OR. fractional_occupation) THEN
133 : resolved = dos_energy_zero_fermi
134 : ELSE
135 100 : resolved = dos_energy_zero_hoco
136 : END IF
137 : ELSE
138 : resolved = energy_zero
139 : END IF
140 :
141 136 : END FUNCTION dos_resolve_energy_zero
142 :
143 : ! **************************************************************************************************
144 : !> \brief Resolve projected-DOS requests from a DOS print section.
145 : !> \param dos_section DOS print section
146 : !> \param do_dos_output whether the DOS print key is active
147 : !> \param do_projected_dos whether any projected DOS output is requested
148 : !> \param do_pdos whether kind-resolved PDOS output is requested
149 : !> \param do_pdos_raw whether raw kind-resolved PDOS output is requested
150 : ! **************************************************************************************************
151 22185 : SUBROUTINE get_dos_pdos_flags(dos_section, do_dos_output, do_projected_dos, do_pdos, do_pdos_raw)
152 :
153 : TYPE(section_vals_type), POINTER :: dos_section
154 : LOGICAL, INTENT(IN) :: do_dos_output
155 : LOGICAL, INTENT(OUT) :: do_projected_dos, do_pdos, do_pdos_raw
156 :
157 : INTEGER :: nrep
158 : LOGICAL :: has_ldos, has_r_ldos
159 : TYPE(section_vals_type), POINTER :: ldos_section, pdos_section
160 :
161 22185 : do_projected_dos = .FALSE.
162 22185 : do_pdos = .FALSE.
163 22185 : do_pdos_raw = .FALSE.
164 22185 : has_ldos = .FALSE.
165 22185 : has_r_ldos = .FALSE.
166 :
167 22185 : IF (do_dos_output) THEN
168 86 : pdos_section => section_vals_get_subs_vals(dos_section, "PDOS")
169 86 : CALL section_vals_get(pdos_section, explicit=do_pdos)
170 86 : IF (do_pdos) CALL section_vals_val_get(pdos_section, "RAW", l_val=do_pdos_raw)
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 22185 : do_projected_dos = do_pdos .OR. has_ldos .OR. has_r_ldos
180 :
181 22185 : 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 2808 : 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 2808 : nhist = SIZE(dos)
208 2808 : istart = MAX(1, FLOOR((eig - broadening_cutoff(broaden_type, broaden_width) - emin)/de) + 1)
209 2808 : iend = MIN(nhist, CEILING((eig + broadening_cutoff(broaden_type, broaden_width) - emin)/de) + 1)
210 72841 : DO i = istart, iend
211 70033 : eval = emin + (i - 1)*de
212 70033 : line_shape = broadening_function(eval - eig, broaden_type, broaden_width, voigt_mixing)
213 70033 : dos(i) = dos(i) + weight*line_shape
214 72841 : occ_dos(i) = occ_dos(i) + weight*occ*line_shape
215 : END DO
216 :
217 2808 : 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 5460 : 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 5460 : nhist = SIZE(curve)
242 5460 : istart = MAX(1, FLOOR((eig - broadening_cutoff(broaden_type, broaden_width) - emin)/de) + 1)
243 5460 : iend = MIN(nhist, CEILING((eig + broadening_cutoff(broaden_type, broaden_width) - emin)/de) + 1)
244 57428 : DO i = istart, iend
245 51968 : eval = emin + (i - 1)*de
246 51968 : line_shape = broadening_function(eval - eig, broaden_type, broaden_width, voigt_mixing)
247 57428 : curve(i) = curve(i) + weight*line_shape
248 : END DO
249 :
250 5460 : 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 8456 : 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 8456 : IF (broaden_type == broadening_gaussian) THEN
265 8456 : cutoff = 8.0_dp*broaden_width
266 : ELSE
267 0 : cutoff = 50.0_dp*broaden_width
268 : END IF
269 :
270 8456 : 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 152009 : 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 152009 : IF (broaden_width <= 0.0_dp) THEN
290 152009 : value = 0.0_dp
291 : RETURN
292 : END IF
293 :
294 152009 : sigma = broaden_width/(2.0_dp*SQRT(2.0_dp*LOG(2.0_dp)))
295 152009 : gamma = 0.5_dp*broaden_width
296 152009 : gaussian_value = EXP(-0.5_dp*(delta_e/sigma)**2)/(sigma*SQRT(2.0_dp*pi))
297 :
298 152009 : 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 152009 : 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 94 : 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 94 : broaden_width_ev = cp_unit_from_cp2k(value=broaden_width, unit_str="eV")
327 :
328 188 : SELECT CASE (broaden_type)
329 : CASE (broadening_gaussian)
330 : WRITE (UNIT=iw, FMT="(A,F10.8,A,F8.6,A)") &
331 94 : "# 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 94 : broaden_width_ev, " eV, Lorentzian fraction = ", MIN(1.0_dp, MAX(0.0_dp, voigt_mixing))
339 : END SELECT
340 :
341 94 : END SUBROUTINE write_broadening_info
342 :
343 : END MODULE qs_dos_utils
|
