(git:b76ce4e)
Loading...
Searching...
No Matches
qs_dos.F
Go to the documentation of this file.
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 Calculation and writing of density of states
10!> \par History
11!> -
12!> \author JGH
13! **************************************************************************************************
14MODULE qs_dos
15 USE cp_array_utils, ONLY: cp_1d_r_p_type
16 USE cp_control_types, ONLY: dft_control_type
17 USE cp_log_handling, ONLY: cp_get_default_logger,&
18 cp_logger_get_default_io_unit,&
19 cp_logger_type
20 USE cp_output_handling, ONLY: cp_p_file,&
21 cp_print_key_finished_output,&
22 cp_print_key_should_output,&
23 cp_print_key_unit_nr
24 USE input_section_types, ONLY: section_vals_type,&
25 section_vals_val_get
26 USE kinds, ONLY: default_string_length,&
27 dp
28 USE kpoint_types, ONLY: kpoint_release,&
29 kpoint_type
30 USE message_passing, ONLY: mp_para_env_type
31 USE qs_band_structure, ONLY: calculate_kp_orbitals
32 USE qs_dos_utils, ONLY: &
33 add_broadened_peak, broadening_cutoff, dos_density_scale, dos_energy_label, &
34 dos_energy_scale, dos_energy_unit_ev, dos_energy_zero_absolute, dos_energy_zero_auto, &
35 dos_energy_zero_hoco, dos_energy_zero_label, dos_resolve_energy_zero, write_broadening_info
36 USE qs_environment_types, ONLY: get_qs_env,&
37 qs_environment_type
38 USE qs_mo_types, ONLY: get_mo_set,&
39 mo_set_type
40#include "./base/base_uses.f90"
41
42 IMPLICIT NONE
43
44 PRIVATE
45
46 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dos'
47
48 PUBLIC :: calculate_dos, calculate_dos_kp
49
50! **************************************************************************************************
51
52CONTAINS
53
54! **************************************************************************************************
55!> \brief Compute and write density of states
56!> \param mos ...
57!> \param dft_section ...
58!> \param unoccupied_evals ...
59!> \param smearing_enabled ...
60!> \date 26.02.2008
61!> \par History:
62!> \author JGH
63!> \version 1.0
64! **************************************************************************************************
65 SUBROUTINE calculate_dos(mos, dft_section, unoccupied_evals, smearing_enabled)
66
67 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
68 TYPE(section_vals_type), POINTER :: dft_section
69 TYPE(cp_1d_r_p_type), DIMENSION(:), OPTIONAL, &
70 POINTER :: unoccupied_evals
71 LOGICAL, INTENT(IN), OPTIONAL :: smearing_enabled
72
73 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_dos'
74
75 CHARACTER(LEN=16) :: energy_label
76 CHARACTER(LEN=20) :: fmtstr_data
77 CHARACTER(LEN=32) :: zero_label
78 CHARACTER(LEN=default_string_length) :: my_act, my_pos
79 INTEGER :: broaden_type, energy_unit, energy_zero, handle, i, iounit, ispin, iterstep, iv, &
80 iw, ndigits, nhist, nmo(2), nspins, nstates(2), nvirt(2), resolved_energy_zero
81 LOGICAL :: append, do_broaden, &
82 fractional_occupation, ionode, &
83 should_output, smear_on
84 REAL(kind=dp) :: broaden_cutoff, broaden_width, de, density_factor, e1, e2, e_fermi(2), &
85 emax, emin, energy_factor, energy_ref(2), ev_factor, eval, hoco(2), out_density, out_occ, &
86 voigt_mixing
87 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: ehist, hist, occval
88 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues, occupation_numbers
89 TYPE(cp_logger_type), POINTER :: logger
90 TYPE(mo_set_type), POINTER :: mo_set
91
92 NULLIFY (logger)
93 logger => cp_get_default_logger()
94 ionode = logger%para_env%is_source()
95 should_output = btest(cp_print_key_should_output(logger%iter_info, dft_section, &
96 "PRINT%DOS"), cp_p_file)
97 iounit = cp_logger_get_default_io_unit(logger)
98 IF ((.NOT. should_output)) RETURN
99
100 CALL timeset(routinen, handle)
101 iterstep = logger%iter_info%iteration(logger%iter_info%n_rlevel)
102
103 IF (iounit > 0) WRITE (unit=iounit, fmt='(/,(T3,A,T61,I10))') &
104 " Calculate DOS at iteration step ", iterstep
105
106 CALL section_vals_val_get(dft_section, "PRINT%DOS%DELTA_E", r_val=de)
107 CALL section_vals_val_get(dft_section, "PRINT%DOS%APPEND", l_val=append)
108 CALL section_vals_val_get(dft_section, "PRINT%DOS%NDIGITS", i_val=ndigits)
109 CALL section_vals_val_get(dft_section, "PRINT%DOS%ENERGY_UNIT", i_val=energy_unit)
110 CALL section_vals_val_get(dft_section, "PRINT%DOS%ENERGY_ZERO", i_val=energy_zero)
111 CALL section_vals_val_get(dft_section, "PRINT%DOS%BROADEN_TYPE", i_val=broaden_type)
112 CALL section_vals_val_get(dft_section, "PRINT%DOS%BROADEN_WIDTH", r_val=broaden_width)
113 CALL section_vals_val_get(dft_section, "PRINT%DOS%VOIGT_MIXING", r_val=voigt_mixing)
114 IF (append .AND. iterstep > 1) THEN
115 my_pos = "APPEND"
116 ELSE
117 my_pos = "REWIND"
118 END IF
119 ndigits = min(max(ndigits, 1), 10)
120 do_broaden = (broaden_width > 0.0_dp)
121 IF (do_broaden) de = max(de, 0.00001_dp)
122
123 emin = 1.e10_dp
124 emax = -1.e10_dp
125 nspins = SIZE(mos)
126 nmo(:) = 0
127 nvirt(:) = 0
128 nstates(:) = 0
129 hoco(:) = -huge(0.0_dp)
130 fractional_occupation = .false.
131 smear_on = .false.
132 IF (PRESENT(smearing_enabled)) smear_on = smearing_enabled
133
134 DO ispin = 1, nspins
135 mo_set => mos(ispin)
136 CALL get_mo_set(mo_set=mo_set, nmo=nmo(ispin), mu=e_fermi(ispin))
137 eigenvalues => mo_set%eigenvalues
138 occupation_numbers => mo_set%occupation_numbers
139 DO i = 1, nmo(ispin)
140 IF (occupation_numbers(i) > 1.0e-10_dp) hoco(ispin) = max(hoco(ispin), eigenvalues(i))
141 IF (abs(occupation_numbers(i) - real(nint(occupation_numbers(i)), kind=dp)) > &
142 1.0e-8_dp) fractional_occupation = .true.
143 END DO
144 IF (hoco(ispin) < -0.5_dp*huge(0.0_dp)) hoco(ispin) = e_fermi(ispin)
145 IF (PRESENT(unoccupied_evals)) THEN
146 IF (ASSOCIATED(unoccupied_evals(ispin)%array)) nvirt(ispin) = SIZE(unoccupied_evals(ispin)%array)
147 END IF
148 nstates(ispin) = nmo(ispin) + nvirt(ispin)
149 e1 = minval(eigenvalues(1:nmo(ispin)))
150 e2 = maxval(eigenvalues(1:nmo(ispin)))
151 IF (nvirt(ispin) > 0) THEN
152 e1 = min(e1, minval(unoccupied_evals(ispin)%array(1:nvirt(ispin))))
153 e2 = max(e2, maxval(unoccupied_evals(ispin)%array(1:nvirt(ispin))))
154 END IF
155 emin = min(emin, e1)
156 emax = max(emax, e2)
157 END DO
158
159 IF (do_broaden) THEN
160 broaden_cutoff = broadening_cutoff(broaden_type, broaden_width)
161 emin = emin - broaden_cutoff
162 emax = emax + broaden_cutoff
163 nhist = nint((emax - emin)/de) + 1
164 ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
165 hist = 0.0_dp
166 occval = 0.0_dp
167 ehist = 0.0_dp
168 DO ispin = 1, nspins
169 mo_set => mos(ispin)
170 occupation_numbers => mo_set%occupation_numbers
171 eigenvalues => mo_set%eigenvalues
172 DO i = 1, nmo(ispin)
173 CALL add_broadened_peak(hist(:, ispin), occval(:, ispin), emin, de, eigenvalues(i), &
174 occupation_numbers(i), 1.0_dp, broaden_type, broaden_width, &
175 voigt_mixing)
176 END DO
177 DO i = 1, nvirt(ispin)
178 CALL add_broadened_peak(hist(:, ispin), occval(:, ispin), emin, de, &
179 unoccupied_evals(ispin)%array(i), 0.0_dp, 1.0_dp, &
180 broaden_type, broaden_width, voigt_mixing)
181 END DO
182 END DO
183 DO i = 1, nhist
184 ehist(i, 1:nspins) = emin + (i - 1)*de
185 END DO
186 ELSE IF (de > 0.0_dp) THEN
187 nhist = nint((emax - emin)/de) + 1
188 ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
189 hist = 0.0_dp
190 occval = 0.0_dp
191 ehist = 0.0_dp
192 DO ispin = 1, nspins
193 mo_set => mos(ispin)
194 occupation_numbers => mo_set%occupation_numbers
195 eigenvalues => mo_set%eigenvalues
196 DO i = 1, nmo(ispin)
197 eval = eigenvalues(i) - emin
198 iv = nint(eval/de) + 1
199 cpassert((iv > 0) .AND. (iv <= nhist))
200 hist(iv, ispin) = hist(iv, ispin) + 1.0_dp
201 occval(iv, ispin) = occval(iv, ispin) + occupation_numbers(i)
202 END DO
203 DO i = 1, nvirt(ispin)
204 eval = unoccupied_evals(ispin)%array(i) - emin
205 iv = nint(eval/de) + 1
206 cpassert((iv > 0) .AND. (iv <= nhist))
207 hist(iv, ispin) = hist(iv, ispin) + 1.0_dp
208 END DO
209 hist(:, ispin) = hist(:, ispin)/real(nstates(ispin), kind=dp)
210 END DO
211 DO i = 1, nhist
212 ehist(i, 1:nspins) = emin + (i - 1)*de
213 END DO
214 ELSE
215 nhist = maxval(nstates)
216 ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
217 hist = 0.0_dp
218 occval = 0.0_dp
219 ehist = 0.0_dp
220 DO ispin = 1, nspins
221 mo_set => mos(ispin)
222 occupation_numbers => mo_set%occupation_numbers
223 eigenvalues => mo_set%eigenvalues
224 DO i = 1, nmo(ispin)
225 ehist(i, ispin) = eigenvalues(i)
226 hist(i, ispin) = 1.0_dp
227 occval(i, ispin) = occupation_numbers(i)
228 END DO
229 DO i = 1, nvirt(ispin)
230 ehist(nmo(ispin) + i, ispin) = unoccupied_evals(ispin)%array(i)
231 hist(nmo(ispin) + i, ispin) = 1.0_dp
232 END DO
233 hist(:, ispin) = hist(:, ispin)/real(nstates(ispin), kind=dp)
234 END DO
235 END IF
236
237 resolved_energy_zero = dos_resolve_energy_zero(energy_zero, smear_on, fractional_occupation)
238 SELECT CASE (resolved_energy_zero)
239 CASE (dos_energy_zero_absolute)
240 energy_ref(:) = 0.0_dp
241 CASE (dos_energy_zero_hoco)
242 energy_ref(:) = hoco(:)
243 CASE DEFAULT
244 energy_ref(:) = e_fermi(:)
245 END SELECT
246 energy_factor = dos_energy_scale(energy_unit)
247 density_factor = dos_density_scale(energy_unit)
248 energy_label = dos_energy_label(energy_unit)
249 zero_label = dos_energy_zero_label(resolved_energy_zero)
250 IF (energy_zero == dos_energy_zero_auto) zero_label = "AUTO -> "//trim(zero_label)
251 ev_factor = dos_energy_scale(dos_energy_unit_ev)
252
253 my_act = "WRITE"
254 iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%DOS", &
255 extension=".dos", file_position=my_pos, file_action=my_act, &
256 file_form="FORMATTED")
257 IF (iw > 0) THEN
258 WRITE (unit=iw, fmt="(A,I0)") "# DOS at iteration step i = ", iterstep
259 IF (nspins == 2) THEN
260 WRITE (unit=iw, fmt="(A,2F12.6,A,2F12.6,A)") &
261 "# E(Fermi) = ", e_fermi(1:2), " a.u. = ", e_fermi(1:2)*ev_factor, " eV"
262 WRITE (unit=iw, fmt="(A,2F12.6,A,2F12.6,A)") &
263 "# E(HOCO) = ", hoco(1:2), " a.u. = ", hoco(1:2)*ev_factor, " eV"
264 WRITE (unit=iw, fmt="(A,A)") "# Energy zero: ", trim(zero_label)
265 IF (do_broaden) CALL write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing)
266 WRITE (unit=iw, fmt="(T2,A,A,A)") "# "//trim(energy_label)//" Alpha_Density Occupation", &
267 " "//trim(energy_label), " Beta_Density Occupation"
268 ! (2(F15.8,2F20.ndigits))
269 WRITE (unit=fmtstr_data, fmt="(A,I0,A)") "(2(F15.8,2F20.", ndigits, "))"
270 ELSE
271 WRITE (unit=iw, fmt="(A,F12.6,A,F12.6,A)") &
272 "# E(Fermi) = ", e_fermi(1), " a.u. = ", e_fermi(1)*ev_factor, " eV"
273 WRITE (unit=iw, fmt="(A,F12.6,A,F12.6,A)") &
274 "# E(HOCO) = ", hoco(1), " a.u. = ", hoco(1)*ev_factor, " eV"
275 WRITE (unit=iw, fmt="(A,A)") "# Energy zero: ", trim(zero_label)
276 IF (do_broaden) CALL write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing)
277 WRITE (unit=iw, fmt="(A,A)") "# "//trim(energy_label), " Density Occupation"
278 ! (F15.8,2F20.ndigits)
279 WRITE (unit=fmtstr_data, fmt="(A,I0,A)") "(F15.8,2F20.", ndigits, ")"
280 END IF
281 DO i = 1, nhist
282 IF (nspins == 2) THEN
283 e1 = (ehist(i, 1) - energy_ref(1))*energy_factor
284 e2 = (ehist(i, 2) - energy_ref(2))*energy_factor
285 ! fmtstr_data == "(2(F15.8,2F20.xx))"
286 IF (do_broaden) THEN
287 WRITE (unit=iw, fmt=fmtstr_data) e1, hist(i, 1)*density_factor, &
288 occval(i, 1)*density_factor, e2, hist(i, 2)*density_factor, &
289 occval(i, 2)*density_factor
290 ELSE
291 WRITE (unit=iw, fmt=fmtstr_data) e1, hist(i, 1), occval(i, 1), &
292 e2, hist(i, 2), occval(i, 2)
293 END IF
294 ELSE
295 eval = (ehist(i, 1) - energy_ref(1))*energy_factor
296 ! fmtstr_data == "(F15.8,2F20.xx)"
297 IF (do_broaden) THEN
298 out_density = hist(i, 1)*density_factor
299 out_occ = occval(i, 1)*density_factor
300 ELSE
301 out_density = hist(i, 1)
302 out_occ = occval(i, 1)
303 END IF
304 WRITE (unit=iw, fmt=fmtstr_data) eval, out_density, out_occ
305 END IF
306 END DO
307 END IF
308 CALL cp_print_key_finished_output(iw, logger, dft_section, "PRINT%DOS")
309 DEALLOCATE (hist, occval, ehist)
310
311 CALL timestop(handle)
312
313 END SUBROUTINE calculate_dos
314
315! **************************************************************************************************
316!> \brief Compute and write density of states (kpoints)
317!> \param qs_env ...
318!> \param dft_section ...
319!> \date 26.02.2008
320!> \par History:
321!> \author JGH
322!> \version 1.0
323! **************************************************************************************************
324 SUBROUTINE calculate_dos_kp(qs_env, dft_section)
325
326 TYPE(qs_environment_type), POINTER :: qs_env
327 TYPE(section_vals_type), POINTER :: dft_section
328
329 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_dos_kp'
330
331 CHARACTER(LEN=16) :: energy_label, fmtstr_data
332 CHARACTER(LEN=32) :: zero_label
333 CHARACTER(LEN=default_string_length) :: my_act, my_pos
334 INTEGER :: broaden_type, energy_unit, energy_zero, fractional_occupation_int, handle, i, ik, &
335 iounit, ispin, iterstep, iv, iw, ndigits, nhist, nmo(2), nmo_kp, nspins, &
336 resolved_energy_zero
337 INTEGER, DIMENSION(:), POINTER :: nkp_grid
338 LOGICAL :: append, do_broaden, explicit, &
339 fractional_occupation, ionode, &
340 should_output
341 REAL(kind=dp) :: broaden_cutoff, broaden_width, de, density_factor, e1, e2, e_fermi(2), &
342 emax, emin, energy_factor, energy_ref(2), ev_factor, eval, hoco(2), out_density, out_occ, &
343 voigt_mixing, wkp
344 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: ehist, hist, occval
345 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues, occupation_numbers
346 TYPE(cp_logger_type), POINTER :: logger
347 TYPE(dft_control_type), POINTER :: dft_control
348 TYPE(kpoint_type), POINTER :: kpoints
349 TYPE(mo_set_type), DIMENSION(:, :), POINTER :: mos
350 TYPE(mo_set_type), POINTER :: mo_set
351 TYPE(mp_para_env_type), POINTER :: para_env
352
353 NULLIFY (logger, kpoints)
354 logger => cp_get_default_logger()
355 ionode = logger%para_env%is_source()
356 should_output = btest(cp_print_key_should_output(logger%iter_info, dft_section, &
357 "PRINT%DOS"), cp_p_file)
358 iounit = cp_logger_get_default_io_unit(logger)
359 IF ((.NOT. should_output)) RETURN
360
361 CALL timeset(routinen, handle)
362 iterstep = logger%iter_info%iteration(logger%iter_info%n_rlevel)
363
364 ! check whether the user requested a different MP grid for the DOS
365 CALL section_vals_val_get(dft_section, "PRINT%DOS%MP_GRID", i_vals=nkp_grid, explicit=explicit)
366
367 IF (explicit) THEN
368 ! make sure is a valid grid
369 DO i = 1, 3
370 IF (nkp_grid(i) < 1) THEN
371 WRITE (unit=iounit, fmt='(T4,A,I3,A,I1)') &
372 "Invalid kpoint grid for DOS ", nkp_grid(i), " in dimension ", i
373 cpabort("")
374 END IF
375 END DO
376 ! calculate orbitals and energies
377 CALL calculate_kp_orbitals(qs_env, kpoints, "MONKHORST-PACK", 0, nkp_grid)
378 ELSE
379 ! use the kpoints from the environment
380 CALL get_qs_env(qs_env, kpoints=kpoints)
381 END IF
382
383 IF (iounit > 0) WRITE (unit=iounit, fmt='(/,(T3,A,T61,I10))') &
384 " Calculate DOS at iteration step ", iterstep
385 IF (iounit > 0) WRITE (unit=iounit, fmt='((T3,A,3I3,A))') &
386 " Using a", kpoints%nkp_grid(:), ' '//trim(kpoints%kp_scheme)//' grid'
387
388 CALL section_vals_val_get(dft_section, "PRINT%DOS%DELTA_E", r_val=de)
389 CALL section_vals_val_get(dft_section, "PRINT%DOS%APPEND", l_val=append)
390 CALL section_vals_val_get(dft_section, "PRINT%DOS%NDIGITS", i_val=ndigits)
391 CALL section_vals_val_get(dft_section, "PRINT%DOS%ENERGY_UNIT", i_val=energy_unit)
392 CALL section_vals_val_get(dft_section, "PRINT%DOS%ENERGY_ZERO", i_val=energy_zero)
393 CALL section_vals_val_get(dft_section, "PRINT%DOS%BROADEN_TYPE", i_val=broaden_type)
394 CALL section_vals_val_get(dft_section, "PRINT%DOS%BROADEN_WIDTH", r_val=broaden_width)
395 CALL section_vals_val_get(dft_section, "PRINT%DOS%VOIGT_MIXING", r_val=voigt_mixing)
396 ! ensure a lower value for the DOS grid width
397 de = max(de, 0.00001_dp)
398 IF (append .AND. iterstep > 1) THEN
399 my_pos = "APPEND"
400 ELSE
401 my_pos = "REWIND"
402 END IF
403 ndigits = min(max(ndigits, 1), 10)
404 do_broaden = (broaden_width > 0.0_dp)
405
406 CALL get_qs_env(qs_env, dft_control=dft_control)
407 nspins = dft_control%nspins
408 para_env => kpoints%para_env_inter_kp
409
410 emin = 1.e10_dp
411 emax = -1.e10_dp
412 nmo(:) = 0
413 e_fermi(:) = 0.0_dp
414 hoco(:) = -huge(0.0_dp)
415 fractional_occupation = .false.
416 IF (kpoints%nkp /= 0) THEN
417 DO ik = 1, SIZE(kpoints%kp_env)
418 mos => kpoints%kp_env(ik)%kpoint_env%mos
419 cpassert(ASSOCIATED(mos))
420 DO ispin = 1, nspins
421 mo_set => mos(1, ispin)
422 CALL get_mo_set(mo_set=mo_set, nmo=nmo_kp, mu=e_fermi(ispin))
423 eigenvalues => mo_set%eigenvalues
424 occupation_numbers => mo_set%occupation_numbers
425 DO i = 1, nmo_kp
426 IF (occupation_numbers(i) > 1.0e-10_dp) hoco(ispin) = max(hoco(ispin), eigenvalues(i))
427 IF (abs(occupation_numbers(i) - real(nint(occupation_numbers(i)), kind=dp)) > &
428 1.0e-8_dp) fractional_occupation = .true.
429 END DO
430 e1 = minval(eigenvalues(1:nmo_kp))
431 e2 = maxval(eigenvalues(1:nmo_kp))
432 emin = min(emin, e1)
433 emax = max(emax, e2)
434 nmo(ispin) = max(nmo(ispin), nmo_kp)
435 END DO
436 END DO
437 END IF
438 CALL para_env%min(emin)
439 CALL para_env%max(emax)
440 CALL para_env%max(nmo)
441 CALL para_env%max(e_fermi)
442 CALL para_env%max(hoco)
443 fractional_occupation_int = merge(1, 0, fractional_occupation)
444 CALL para_env%max(fractional_occupation_int)
445 fractional_occupation = (fractional_occupation_int /= 0)
446 DO ispin = 1, nspins
447 IF (hoco(ispin) < -0.5_dp*huge(0.0_dp)) hoco(ispin) = e_fermi(ispin)
448 END DO
449
450 IF (do_broaden) THEN
451 broaden_cutoff = broadening_cutoff(broaden_type, broaden_width)
452 emin = emin - broaden_cutoff
453 emax = emax + broaden_cutoff
454 END IF
455 nhist = nint((emax - emin)/de) + 1
456 ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
457 hist = 0.0_dp
458 occval = 0.0_dp
459 ehist = 0.0_dp
460
461 IF (kpoints%nkp /= 0) THEN
462 DO ik = 1, SIZE(kpoints%kp_env)
463 mos => kpoints%kp_env(ik)%kpoint_env%mos
464 wkp = kpoints%kp_env(ik)%kpoint_env%wkp
465 DO ispin = 1, nspins
466 mo_set => mos(1, ispin)
467 occupation_numbers => mo_set%occupation_numbers
468 eigenvalues => mo_set%eigenvalues
469 IF (do_broaden) THEN
470 DO i = 1, nmo(ispin)
471 CALL add_broadened_peak(hist(:, ispin), occval(:, ispin), emin, de, eigenvalues(i), &
472 occupation_numbers(i), wkp, broaden_type, broaden_width, &
473 voigt_mixing)
474 END DO
475 ELSE
476 DO i = 1, nmo(ispin)
477 eval = eigenvalues(i) - emin
478 iv = nint(eval/de) + 1
479 cpassert((iv > 0) .AND. (iv <= nhist))
480 hist(iv, ispin) = hist(iv, ispin) + wkp
481 occval(iv, ispin) = occval(iv, ispin) + wkp*occupation_numbers(i)
482 END DO
483 END IF
484 END DO
485 END DO
486 END IF
487 CALL para_env%sum(hist)
488 CALL para_env%sum(occval)
489 IF (.NOT. do_broaden) THEN
490 DO ispin = 1, nspins
491 hist(:, ispin) = hist(:, ispin)/real(nmo(ispin), kind=dp)
492 END DO
493 END IF
494 DO i = 1, nhist
495 ehist(i, 1:nspins) = emin + (i - 1)*de
496 END DO
497
498 resolved_energy_zero = dos_resolve_energy_zero(energy_zero, dft_control%smear, fractional_occupation)
499 SELECT CASE (resolved_energy_zero)
500 CASE (dos_energy_zero_absolute)
501 energy_ref(:) = 0.0_dp
502 CASE (dos_energy_zero_hoco)
503 energy_ref(:) = hoco(:)
504 CASE DEFAULT
505 energy_ref(:) = e_fermi(:)
506 END SELECT
507 energy_factor = dos_energy_scale(energy_unit)
508 density_factor = dos_density_scale(energy_unit)
509 energy_label = dos_energy_label(energy_unit)
510 zero_label = dos_energy_zero_label(resolved_energy_zero)
511 IF (energy_zero == dos_energy_zero_auto) zero_label = "AUTO -> "//trim(zero_label)
512 ev_factor = dos_energy_scale(dos_energy_unit_ev)
513
514 my_act = "WRITE"
515 iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%DOS", &
516 extension=".dos", file_position=my_pos, file_action=my_act, &
517 file_form="FORMATTED")
518 IF (iw > 0) THEN
519 WRITE (unit=iw, fmt="(A,I0)") "# DOS at iteration step i = ", iterstep
520 IF (nspins == 2) THEN
521 WRITE (unit=iw, fmt="(A,2F12.6,A,2F12.6,A)") &
522 "# E(Fermi) = ", e_fermi(1:2), " a.u. = ", e_fermi(1:2)*ev_factor, " eV"
523 WRITE (unit=iw, fmt="(A,2F12.6,A,2F12.6,A)") &
524 "# E(HOCO) = ", hoco(1:2), " a.u. = ", hoco(1:2)*ev_factor, " eV"
525 WRITE (unit=iw, fmt="(A,A)") "# Energy zero: ", trim(zero_label)
526 IF (do_broaden) CALL write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing)
527 WRITE (unit=iw, fmt="(A,A)") "# "//trim(energy_label)//" Alpha_Density Occupation", &
528 " Beta_Density Occupation"
529 ! (F15.8,4F20.ndigits)
530 WRITE (unit=fmtstr_data, fmt="(A,I0,A)") "(F15.8,4F20.", ndigits, ")"
531 ELSE
532 WRITE (unit=iw, fmt="(A,F12.6,A,F12.6,A)") &
533 "# E(Fermi) = ", e_fermi(1), " a.u. = ", e_fermi(1)*ev_factor, " eV"
534 WRITE (unit=iw, fmt="(A,F12.6,A,F12.6,A)") &
535 "# E(HOCO) = ", hoco(1), " a.u. = ", hoco(1)*ev_factor, " eV"
536 WRITE (unit=iw, fmt="(A,A)") "# Energy zero: ", trim(zero_label)
537 IF (do_broaden) CALL write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing)
538 WRITE (unit=iw, fmt="(A,A)") "# "//trim(energy_label), " Density Occupation"
539 ! (F15.8,2F20.ndigits)
540 WRITE (unit=fmtstr_data, fmt="(A,I0,A)") "(F15.8,2F20.", ndigits, ")"
541 END IF
542 DO i = 1, nhist
543 eval = (ehist(i, 1) - energy_ref(1))*energy_factor
544 IF (nspins == 2) THEN
545 ! fmtstr_data == "(F15.8,4F20.xx)"
546 IF (do_broaden) THEN
547 WRITE (unit=iw, fmt=fmtstr_data) eval, hist(i, 1)*density_factor, &
548 occval(i, 1)*density_factor, hist(i, 2)*density_factor, occval(i, 2)*density_factor
549 ELSE
550 WRITE (unit=iw, fmt=fmtstr_data) eval, hist(i, 1), occval(i, 1), &
551 hist(i, 2), occval(i, 2)
552 END IF
553 ELSE
554 ! fmtstr_data == "(F15.8,2F20.xx)"
555 IF (do_broaden) THEN
556 out_density = hist(i, 1)*density_factor
557 out_occ = occval(i, 1)*density_factor
558 ELSE
559 out_density = hist(i, 1)
560 out_occ = occval(i, 1)
561 END IF
562 WRITE (unit=iw, fmt=fmtstr_data) eval, out_density, out_occ
563 END IF
564 END DO
565 END IF
566 CALL cp_print_key_finished_output(iw, logger, dft_section, "PRINT%DOS")
567 DEALLOCATE (hist, occval, ehist)
568
569 ! destroy the extra k-point set if it was created
570 IF (explicit) THEN
571 CALL kpoint_release(kpoints)
572 END IF
573
574 CALL timestop(handle)
575
576 END SUBROUTINE calculate_dos_kp
577
578END MODULE qs_dos
cp_array_utils
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Definition cp_array_utils.F:18
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:515
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1067
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:345
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1047
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
kpoint_types::kpoint_release
subroutine, public kpoint_release(kpoint)
Release a kpoint environment, deallocate all data.
Definition kpoint_types.F:284
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
qs_band_structure
Calculation of band structures.
Definition qs_band_structure.F:14
qs_band_structure::calculate_kp_orbitals
subroutine, public calculate_kp_orbitals(qs_env, kpoint, scheme, nadd, mp_grid, kpgeneral, group_size_ext, kp_shift, gamma_centered)
diagonalize KS matrices at a set of kpoints
Definition qs_band_structure.F:309
qs_dos_utils
Utilities for broadened DOS and PDOS output.
Definition qs_dos_utils.F:11
qs_dos_utils::add_broadened_peak
subroutine, public add_broadened_peak(dos, occ_dos, emin, de, eig, occ, weight, broaden_type, broaden_width, voigt_mixing)
Add a broadened spectral line to a DOS curve.
Definition qs_dos_utils.F:198
qs_dos_utils::write_broadening_info
subroutine, public write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing)
Write broadening metadata.
Definition qs_dos_utils.F:320
qs_dos_utils::dos_energy_zero_hoco
integer, parameter, public dos_energy_zero_hoco
Definition qs_dos_utils.F:32
qs_dos_utils::dos_energy_zero_auto
integer, parameter, public dos_energy_zero_auto
Definition qs_dos_utils.F:32
qs_dos_utils::dos_energy_unit_ev
integer, parameter, public dos_energy_unit_ev
Definition qs_dos_utils.F:30
qs_dos_utils::dos_energy_label
character(len=16) function, public dos_energy_label(energy_unit)
Return the energy-column label for DOS-like output.
Definition qs_dos_utils.F:82
qs_dos_utils::dos_resolve_energy_zero
integer function, public dos_resolve_energy_zero(energy_zero, smearing_enabled, fractional_occupation)
Resolve AUTO energy-zero selection for DOS-like output.
Definition qs_dos_utils.F:126
qs_dos_utils::dos_energy_scale
real(kind=dp) function, public dos_energy_scale(energy_unit)
Return the conversion factor from internal energy units to the selected DOS energy unit.
Definition qs_dos_utils.F:49
qs_dos_utils::dos_energy_zero_absolute
integer, parameter, public dos_energy_zero_absolute
Definition qs_dos_utils.F:32
qs_dos_utils::broadening_cutoff
pure real(kind=dp) function, public broadening_cutoff(broaden_type, broaden_width)
Broadening cutoff used for numerical accumulation.
Definition qs_dos_utils.F:259
qs_dos_utils::dos_density_scale
real(kind=dp) function, public dos_density_scale(energy_unit)
Return the DOS-density conversion factor for the selected energy unit.
Definition qs_dos_utils.F:68
qs_dos_utils::dos_energy_zero_label
character(len=16) function, public dos_energy_zero_label(energy_zero)
Return the label for the selected DOS energy zero.
Definition qs_dos_utils.F:101
qs_dos
Calculation and writing of density of states.
Definition qs_dos.F:14
qs_dos::calculate_dos
subroutine, public calculate_dos(mos, dft_section, unoccupied_evals, smearing_enabled)
Compute and write density of states.
Definition qs_dos.F:66
qs_dos::calculate_dos_kp
subroutine, public calculate_dos_kp(qs_env, dft_section)
Compute and write density of states (kpoints)
Definition qs_dos.F:325
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:530
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition qs_mo_types.F:401
cp_array_utils::cp_1d_r_p_type
represent a pointer to a 1d array
Definition cp_array_utils.F:99
cp_control_types::dft_control_type
Definition cp_control_types.F:744
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:171
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:743
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:220
qs_mo_types::mo_set_type
Definition qs_mo_types.F:47