(git:fdbe441)
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!> \param write_curve_output ...
61!> \date 26.02.2008
62!> \par History:
63!> \author JGH
64!> \version 1.0
65! **************************************************************************************************
66 SUBROUTINE calculate_dos(mos, dft_section, unoccupied_evals, smearing_enabled, write_curve_output)
67
68 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
69 TYPE(section_vals_type), POINTER :: dft_section
70 TYPE(cp_1d_r_p_type), DIMENSION(:), OPTIONAL, &
71 POINTER :: unoccupied_evals
72 LOGICAL, INTENT(IN), OPTIONAL :: smearing_enabled, write_curve_output
73
74 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_dos'
75
76 CHARACTER(LEN=16) :: energy_label
77 CHARACTER(LEN=20) :: fmtstr_data
78 CHARACTER(LEN=32) :: zero_label
79 CHARACTER(LEN=default_string_length) :: my_act, my_pos
80 INTEGER :: broaden_type, energy_unit, energy_zero, handle, i, iounit, ispin, iterstep, iv, &
81 iw, ndigits, nhist, nmo(2), nspins, nstates(2), nvirt(2), resolved_energy_zero
82 LOGICAL :: append, do_broaden, &
83 fractional_occupation, ionode, &
84 should_output, smear_on
85 REAL(kind=dp) :: broaden_cutoff, broaden_width, de, density_factor, e1, e2, e_fermi(2), &
86 emax, emin, energy_factor, energy_ref(2), ev_factor, eval, hoco(2), out_density, out_occ, &
87 voigt_mixing
88 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: ehist, hist, occval
89 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues, occupation_numbers
90 TYPE(cp_logger_type), POINTER :: logger
91 TYPE(mo_set_type), POINTER :: mo_set
92
93 NULLIFY (logger)
94 logger => cp_get_default_logger()
95 ionode = logger%para_env%is_source()
96 should_output = btest(cp_print_key_should_output(logger%iter_info, dft_section, &
97 "PRINT%DOS"), cp_p_file)
98 iounit = cp_logger_get_default_io_unit(logger)
99 IF ((.NOT. should_output)) RETURN
100
101 CALL timeset(routinen, handle)
102 iterstep = logger%iter_info%iteration(logger%iter_info%n_rlevel)
103
104 IF (iounit > 0) WRITE (unit=iounit, fmt='(/,(T3,A,T61,I10))') &
105 " Calculate DOS at iteration step ", iterstep
106
107 CALL section_vals_val_get(dft_section, "PRINT%DOS%DELTA_E", r_val=de)
108 CALL section_vals_val_get(dft_section, "PRINT%DOS%APPEND", l_val=append)
109 CALL section_vals_val_get(dft_section, "PRINT%DOS%NDIGITS", i_val=ndigits)
110 CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%ENERGY_UNIT", i_val=energy_unit)
111 CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%ENERGY_ZERO", i_val=energy_zero)
112 CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%BROADEN%TYPE", i_val=broaden_type)
113 CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%BROADEN%WIDTH", r_val=broaden_width)
114 CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%BROADEN%VOIGT_MIXING", r_val=voigt_mixing)
115 IF (append .AND. iterstep > 1) THEN
116 my_pos = "APPEND"
117 ELSE
118 my_pos = "REWIND"
119 END IF
120 ndigits = min(max(ndigits, 1), 10)
121 IF (PRESENT(write_curve_output)) THEN
122 IF (write_curve_output .AND. de <= 0.0_dp) THEN
123 cpwarn("Broadened DOS output requires DELTA_E > 0 and will be skipped")
124 CALL timestop(handle)
125 RETURN
126 END IF
127 IF (write_curve_output .AND. broaden_width <= 0.0_dp) THEN
128 cpwarn("Broadened DOS output requires a finite WIDTH and will be skipped")
129 CALL timestop(handle)
130 RETURN
131 END IF
132 END IF
133 do_broaden = .false.
134 IF (PRESENT(write_curve_output)) do_broaden = write_curve_output
135 do_broaden = do_broaden .AND. (broaden_width > 0.0_dp)
136 IF (do_broaden) de = max(de, 0.00001_dp)
137
138 emin = 1.e10_dp
139 emax = -1.e10_dp
140 nspins = SIZE(mos)
141 nmo(:) = 0
142 nvirt(:) = 0
143 nstates(:) = 0
144 hoco(:) = -huge(0.0_dp)
145 fractional_occupation = .false.
146 smear_on = .false.
147 IF (PRESENT(smearing_enabled)) smear_on = smearing_enabled
148
149 DO ispin = 1, nspins
150 mo_set => mos(ispin)
151 CALL get_mo_set(mo_set=mo_set, nmo=nmo(ispin), mu=e_fermi(ispin))
152 eigenvalues => mo_set%eigenvalues
153 occupation_numbers => mo_set%occupation_numbers
154 DO i = 1, nmo(ispin)
155 IF (occupation_numbers(i) > 1.0e-10_dp) hoco(ispin) = max(hoco(ispin), eigenvalues(i))
156 IF (abs(occupation_numbers(i) - real(nint(occupation_numbers(i)), kind=dp)) > &
157 1.0e-8_dp) fractional_occupation = .true.
158 END DO
159 IF (hoco(ispin) < -0.5_dp*huge(0.0_dp)) hoco(ispin) = e_fermi(ispin)
160 IF (PRESENT(unoccupied_evals)) THEN
161 IF (ASSOCIATED(unoccupied_evals(ispin)%array)) nvirt(ispin) = SIZE(unoccupied_evals(ispin)%array)
162 END IF
163 nstates(ispin) = nmo(ispin) + nvirt(ispin)
164 e1 = minval(eigenvalues(1:nmo(ispin)))
165 e2 = maxval(eigenvalues(1:nmo(ispin)))
166 IF (nvirt(ispin) > 0) THEN
167 e1 = min(e1, minval(unoccupied_evals(ispin)%array(1:nvirt(ispin))))
168 e2 = max(e2, maxval(unoccupied_evals(ispin)%array(1:nvirt(ispin))))
169 END IF
170 emin = min(emin, e1)
171 emax = max(emax, e2)
172 END DO
173
174 IF (do_broaden) THEN
175 broaden_cutoff = broadening_cutoff(broaden_type, broaden_width)
176 emin = emin - broaden_cutoff
177 emax = emax + broaden_cutoff
178 nhist = nint((emax - emin)/de) + 1
179 ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
180 hist = 0.0_dp
181 occval = 0.0_dp
182 ehist = 0.0_dp
183 DO ispin = 1, nspins
184 mo_set => mos(ispin)
185 occupation_numbers => mo_set%occupation_numbers
186 eigenvalues => mo_set%eigenvalues
187 DO i = 1, nmo(ispin)
188 CALL add_broadened_peak(hist(:, ispin), occval(:, ispin), emin, de, eigenvalues(i), &
189 occupation_numbers(i), 1.0_dp, broaden_type, broaden_width, &
190 voigt_mixing)
191 END DO
192 DO i = 1, nvirt(ispin)
193 CALL add_broadened_peak(hist(:, ispin), occval(:, ispin), emin, de, &
194 unoccupied_evals(ispin)%array(i), 0.0_dp, 1.0_dp, &
195 broaden_type, broaden_width, voigt_mixing)
196 END DO
197 END DO
198 DO i = 1, nhist
199 ehist(i, 1:nspins) = emin + (i - 1)*de
200 END DO
201 ELSE IF (de > 0.0_dp) THEN
202 nhist = nint((emax - emin)/de) + 1
203 ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
204 hist = 0.0_dp
205 occval = 0.0_dp
206 ehist = 0.0_dp
207 DO ispin = 1, nspins
208 mo_set => mos(ispin)
209 occupation_numbers => mo_set%occupation_numbers
210 eigenvalues => mo_set%eigenvalues
211 DO i = 1, nmo(ispin)
212 eval = eigenvalues(i) - emin
213 iv = nint(eval/de) + 1
214 cpassert((iv > 0) .AND. (iv <= nhist))
215 hist(iv, ispin) = hist(iv, ispin) + 1.0_dp
216 occval(iv, ispin) = occval(iv, ispin) + occupation_numbers(i)
217 END DO
218 DO i = 1, nvirt(ispin)
219 eval = unoccupied_evals(ispin)%array(i) - emin
220 iv = nint(eval/de) + 1
221 cpassert((iv > 0) .AND. (iv <= nhist))
222 hist(iv, ispin) = hist(iv, ispin) + 1.0_dp
223 END DO
224 hist(:, ispin) = hist(:, ispin)/real(nstates(ispin), kind=dp)
225 END DO
226 DO i = 1, nhist
227 ehist(i, 1:nspins) = emin + (i - 1)*de
228 END DO
229 ELSE
230 nhist = maxval(nstates)
231 ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
232 hist = 0.0_dp
233 occval = 0.0_dp
234 ehist = 0.0_dp
235 DO ispin = 1, nspins
236 mo_set => mos(ispin)
237 occupation_numbers => mo_set%occupation_numbers
238 eigenvalues => mo_set%eigenvalues
239 DO i = 1, nmo(ispin)
240 ehist(i, ispin) = eigenvalues(i)
241 hist(i, ispin) = 1.0_dp
242 occval(i, ispin) = occupation_numbers(i)
243 END DO
244 DO i = 1, nvirt(ispin)
245 ehist(nmo(ispin) + i, ispin) = unoccupied_evals(ispin)%array(i)
246 hist(nmo(ispin) + i, ispin) = 1.0_dp
247 END DO
248 hist(:, ispin) = hist(:, ispin)/real(nstates(ispin), kind=dp)
249 END DO
250 END IF
251
252 resolved_energy_zero = dos_resolve_energy_zero(energy_zero, smear_on, fractional_occupation)
253 SELECT CASE (resolved_energy_zero)
254 CASE (dos_energy_zero_absolute)
255 energy_ref(:) = 0.0_dp
256 CASE (dos_energy_zero_hoco)
257 energy_ref(:) = maxval(hoco(1:nspins))
258 CASE DEFAULT
259 energy_ref(:) = maxval(e_fermi(1:nspins))
260 END SELECT
261 IF (.NOT. do_broaden) energy_ref(:) = 0.0_dp
262 energy_factor = merge(dos_energy_scale(energy_unit), 1.0_dp, do_broaden)
263 density_factor = merge(dos_density_scale(energy_unit), 1.0_dp, do_broaden)
264 energy_label = dos_energy_label(merge(energy_unit, 1, do_broaden))
265 zero_label = dos_energy_zero_label(resolved_energy_zero)
266 IF (energy_zero == dos_energy_zero_auto) zero_label = "AUTO -> "//trim(zero_label)
267 ev_factor = dos_energy_scale(dos_energy_unit_ev)
268
269 my_act = "WRITE"
270 IF (do_broaden) THEN
271 iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%DOS", &
272 extension=".dos", file_position=my_pos, file_action=my_act, &
273 file_form="FORMATTED", middle_name="curve")
274 ELSE
275 iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%DOS", &
276 extension=".dos", file_position=my_pos, file_action=my_act, &
277 file_form="FORMATTED")
278 END IF
279 IF (iw > 0) THEN
280 WRITE (unit=iw, fmt="(A,I0)") "# DOS at iteration step i = ", iterstep
281 IF (nspins == 2) THEN
282 WRITE (unit=iw, fmt="(A,2F12.6,A,2F12.6,A)") &
283 "# E(Fermi) = ", e_fermi(1:2), " a.u. = ", e_fermi(1:2)*ev_factor, " eV"
284 WRITE (unit=iw, fmt="(A,2F12.6,A,2F12.6,A)") &
285 "# E(HOCO) = ", hoco(1:2), " a.u. = ", hoco(1:2)*ev_factor, " eV"
286 IF (do_broaden) THEN
287 WRITE (unit=iw, fmt="(A,A)") "# Energy zero: ", trim(zero_label)
288 CALL write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing)
289 END IF
290 IF (do_broaden .OR. de > 0.0_dp) THEN
291 WRITE (unit=iw, fmt="(A,A,A)") "# "//trim(energy_label)//" Alpha_Density Occupation", &
292 " Beta_Density Occupation"
293 WRITE (unit=fmtstr_data, fmt="(A,I0,A)") "(F15.8,4F20.", ndigits, ")"
294 ELSE
295 WRITE (unit=iw, fmt="(A,A,A)") "# "//trim(energy_label)//" Alpha_Density Occupation", &
296 " "//trim(energy_label), " Beta_Density Occupation"
297 WRITE (unit=fmtstr_data, fmt="(A,I0,A)") "(2(F15.8,2F20.", ndigits, "))"
298 END IF
299 ELSE
300 WRITE (unit=iw, fmt="(A,F12.6,A,F12.6,A)") &
301 "# E(Fermi) = ", e_fermi(1), " a.u. = ", e_fermi(1)*ev_factor, " eV"
302 WRITE (unit=iw, fmt="(A,F12.6,A,F12.6,A)") &
303 "# E(HOCO) = ", hoco(1), " a.u. = ", hoco(1)*ev_factor, " eV"
304 IF (do_broaden) THEN
305 WRITE (unit=iw, fmt="(A,A)") "# Energy zero: ", trim(zero_label)
306 CALL write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing)
307 END IF
308 WRITE (unit=iw, fmt="(A,A)") "# "//trim(energy_label), " Density Occupation"
309 ! (F15.8,2F20.ndigits)
310 WRITE (unit=fmtstr_data, fmt="(A,I0,A)") "(F15.8,2F20.", ndigits, ")"
311 END IF
312 DO i = 1, nhist
313 IF (nspins == 2) THEN
314 IF (do_broaden) THEN
315 eval = (ehist(i, 1) - energy_ref(1))*energy_factor
316 WRITE (unit=iw, fmt=fmtstr_data) eval, hist(i, 1)*density_factor, &
317 occval(i, 1)*density_factor, hist(i, 2)*density_factor, &
318 occval(i, 2)*density_factor
319 ELSE IF (de > 0.0_dp) THEN
320 IF (hist(i, 1) == 0.0_dp .AND. occval(i, 1) == 0.0_dp .AND. &
321 hist(i, 2) == 0.0_dp .AND. occval(i, 2) == 0.0_dp) cycle
322 eval = ehist(i, 1)
323 WRITE (unit=iw, fmt=fmtstr_data) eval, hist(i, 1), occval(i, 1), &
324 hist(i, 2), occval(i, 2)
325 ELSE
326 e1 = ehist(i, 1)
327 e2 = ehist(i, 2)
328 WRITE (unit=iw, fmt=fmtstr_data) e1, hist(i, 1), occval(i, 1), &
329 e2, hist(i, 2), occval(i, 2)
330 END IF
331 ELSE
332 eval = (ehist(i, 1) - energy_ref(1))*energy_factor
333 ! fmtstr_data == "(F15.8,2F20.xx)"
334 IF (do_broaden) THEN
335 out_density = hist(i, 1)*density_factor
336 out_occ = occval(i, 1)*density_factor
337 ELSE
338 out_density = hist(i, 1)
339 out_occ = occval(i, 1)
340 IF (out_density == 0.0_dp .AND. out_occ == 0.0_dp) cycle
341 END IF
342 WRITE (unit=iw, fmt=fmtstr_data) eval, out_density, out_occ
343 END IF
344 END DO
345 END IF
346 CALL cp_print_key_finished_output(iw, logger, dft_section, "PRINT%DOS")
347 DEALLOCATE (hist, occval, ehist)
348
349 CALL timestop(handle)
350
351 END SUBROUTINE calculate_dos
352
353! **************************************************************************************************
354!> \brief Compute and write density of states (kpoints)
355!> \param qs_env ...
356!> \param dft_section ...
357!> \param write_curve_output ...
358!> \date 26.02.2008
359!> \par History:
360!> \author JGH
361!> \version 1.0
362! **************************************************************************************************
363 SUBROUTINE calculate_dos_kp(qs_env, dft_section, write_curve_output)
364
365 TYPE(qs_environment_type), POINTER :: qs_env
366 TYPE(section_vals_type), POINTER :: dft_section
367 LOGICAL, INTENT(IN), OPTIONAL :: write_curve_output
368
369 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_dos_kp'
370
371 CHARACTER(LEN=16) :: energy_label, fmtstr_data
372 CHARACTER(LEN=32) :: zero_label
373 CHARACTER(LEN=default_string_length) :: err, my_act, my_pos
374 INTEGER :: broaden_type, energy_unit, energy_zero, fractional_occupation_int, handle, i, ik, &
375 iounit, ispin, iterstep, iv, iw, ndigits, nhist, nmo(2), nmo_kp, nspins, &
376 resolved_energy_zero
377 INTEGER, DIMENSION(:), POINTER :: nkp_grid
378 LOGICAL :: append, do_broaden, explicit, &
379 fractional_occupation, ionode, &
380 should_output
381 REAL(kind=dp) :: broaden_cutoff, broaden_width, de, density_factor, e1, e2, e_fermi(2), &
382 emax, emin, energy_factor, energy_ref(2), ev_factor, eval, hoco(2), out_density, out_occ, &
383 voigt_mixing, wkp
384 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: ehist, hist, occval
385 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues, occupation_numbers
386 TYPE(cp_logger_type), POINTER :: logger
387 TYPE(dft_control_type), POINTER :: dft_control
388 TYPE(kpoint_type), POINTER :: kpoints
389 TYPE(mo_set_type), DIMENSION(:, :), POINTER :: mos
390 TYPE(mo_set_type), POINTER :: mo_set
391 TYPE(mp_para_env_type), POINTER :: para_env
392
393 NULLIFY (logger, kpoints)
394 logger => cp_get_default_logger()
395 ionode = logger%para_env%is_source()
396 should_output = btest(cp_print_key_should_output(logger%iter_info, dft_section, &
397 "PRINT%DOS"), cp_p_file)
398 iounit = cp_logger_get_default_io_unit(logger)
399 IF ((.NOT. should_output)) RETURN
400
401 CALL timeset(routinen, handle)
402 iterstep = logger%iter_info%iteration(logger%iter_info%n_rlevel)
403
404 ! check whether the user requested a different MP grid for the DOS
405 CALL section_vals_val_get(dft_section, "PRINT%DOS%MP_GRID", i_vals=nkp_grid, explicit=explicit)
406
407 IF (explicit) THEN
408 ! make sure is a valid grid
409 DO i = 1, 3
410 IF (nkp_grid(i) < 1) THEN
411 WRITE (unit=err, fmt='(T4,A,I3,A,I1)') &
412 "Invalid kpoint grid for DOS ", nkp_grid(i), " in dimension ", i
413 cpabort(trim(err))
414 END IF
415 END DO
416 ! calculate orbitals and energies
417 CALL calculate_kp_orbitals(qs_env, kpoints, "MONKHORST-PACK", 0, nkp_grid)
418 ELSE
419 ! use the kpoints from the environment
420 CALL get_qs_env(qs_env, kpoints=kpoints)
421 END IF
422
423 IF (iounit > 0) WRITE (unit=iounit, fmt='(/,(T3,A,T61,I10))') &
424 " Calculate DOS at iteration step ", iterstep
425 IF (iounit > 0) WRITE (unit=iounit, fmt='((T3,A,3I3,A))') &
426 " Using a", kpoints%nkp_grid(:), ' '//trim(kpoints%kp_scheme)//' grid'
427
428 CALL section_vals_val_get(dft_section, "PRINT%DOS%DELTA_E", r_val=de)
429 CALL section_vals_val_get(dft_section, "PRINT%DOS%APPEND", l_val=append)
430 CALL section_vals_val_get(dft_section, "PRINT%DOS%NDIGITS", i_val=ndigits)
431 CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%ENERGY_UNIT", i_val=energy_unit)
432 CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%ENERGY_ZERO", i_val=energy_zero)
433 CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%BROADEN%TYPE", i_val=broaden_type)
434 CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%BROADEN%WIDTH", r_val=broaden_width)
435 CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%BROADEN%VOIGT_MIXING", r_val=voigt_mixing)
436 IF (append .AND. iterstep > 1) THEN
437 my_pos = "APPEND"
438 ELSE
439 my_pos = "REWIND"
440 END IF
441 ndigits = min(max(ndigits, 1), 10)
442 IF (PRESENT(write_curve_output)) THEN
443 IF (write_curve_output .AND. de <= 0.0_dp) THEN
444 cpwarn("Broadened k-point DOS output requires DELTA_E > 0 and will be skipped")
445 CALL timestop(handle)
446 IF (explicit) CALL kpoint_release(kpoints)
447 RETURN
448 END IF
449 IF (write_curve_output .AND. broaden_width <= 0.0_dp) THEN
450 cpwarn("Broadened k-point DOS output requires a finite WIDTH and will be skipped")
451 CALL timestop(handle)
452 IF (explicit) CALL kpoint_release(kpoints)
453 RETURN
454 END IF
455 ELSE IF (de <= 0.0_dp) THEN
456 cpwarn("K-point DOS output requires DELTA_E > 0 and will be skipped")
457 CALL timestop(handle)
458 IF (explicit) CALL kpoint_release(kpoints)
459 RETURN
460 END IF
461 ! ensure a lower value for the DOS grid width
462 de = max(de, 0.00001_dp)
463 do_broaden = .false.
464 IF (PRESENT(write_curve_output)) do_broaden = write_curve_output
465 do_broaden = do_broaden .AND. (broaden_width > 0.0_dp)
466
467 CALL get_qs_env(qs_env, dft_control=dft_control)
468 nspins = dft_control%nspins
469 para_env => kpoints%para_env_inter_kp
470
471 emin = 1.e10_dp
472 emax = -1.e10_dp
473 nmo(:) = 0
474 e_fermi(:) = 0.0_dp
475 hoco(:) = -huge(0.0_dp)
476 fractional_occupation = .false.
477 IF (kpoints%nkp /= 0) THEN
478 DO ik = 1, SIZE(kpoints%kp_env)
479 mos => kpoints%kp_env(ik)%kpoint_env%mos
480 cpassert(ASSOCIATED(mos))
481 DO ispin = 1, nspins
482 mo_set => mos(1, ispin)
483 CALL get_mo_set(mo_set=mo_set, nmo=nmo_kp, mu=e_fermi(ispin))
484 eigenvalues => mo_set%eigenvalues
485 occupation_numbers => mo_set%occupation_numbers
486 DO i = 1, nmo_kp
487 IF (occupation_numbers(i) > 1.0e-10_dp) hoco(ispin) = max(hoco(ispin), eigenvalues(i))
488 IF (abs(occupation_numbers(i) - real(nint(occupation_numbers(i)), kind=dp)) > &
489 1.0e-8_dp) fractional_occupation = .true.
490 END DO
491 e1 = minval(eigenvalues(1:nmo_kp))
492 e2 = maxval(eigenvalues(1:nmo_kp))
493 emin = min(emin, e1)
494 emax = max(emax, e2)
495 nmo(ispin) = max(nmo(ispin), nmo_kp)
496 END DO
497 END DO
498 END IF
499 CALL para_env%min(emin)
500 CALL para_env%max(emax)
501 CALL para_env%max(nmo)
502 CALL para_env%max(e_fermi)
503 CALL para_env%max(hoco)
504 fractional_occupation_int = merge(1, 0, fractional_occupation)
505 CALL para_env%max(fractional_occupation_int)
506 fractional_occupation = (fractional_occupation_int /= 0)
507 DO ispin = 1, nspins
508 IF (hoco(ispin) < -0.5_dp*huge(0.0_dp)) hoco(ispin) = e_fermi(ispin)
509 END DO
510
511 IF (do_broaden) THEN
512 broaden_cutoff = broadening_cutoff(broaden_type, broaden_width)
513 emin = emin - broaden_cutoff
514 emax = emax + broaden_cutoff
515 END IF
516 nhist = nint((emax - emin)/de) + 1
517 ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
518 hist = 0.0_dp
519 occval = 0.0_dp
520 ehist = 0.0_dp
521
522 IF (kpoints%nkp /= 0) THEN
523 DO ik = 1, SIZE(kpoints%kp_env)
524 mos => kpoints%kp_env(ik)%kpoint_env%mos
525 wkp = kpoints%kp_env(ik)%kpoint_env%wkp
526 DO ispin = 1, nspins
527 mo_set => mos(1, ispin)
528 occupation_numbers => mo_set%occupation_numbers
529 eigenvalues => mo_set%eigenvalues
530 IF (do_broaden) THEN
531 DO i = 1, nmo(ispin)
532 CALL add_broadened_peak(hist(:, ispin), occval(:, ispin), emin, de, eigenvalues(i), &
533 occupation_numbers(i), wkp, broaden_type, broaden_width, &
534 voigt_mixing)
535 END DO
536 ELSE
537 DO i = 1, nmo(ispin)
538 eval = eigenvalues(i) - emin
539 iv = nint(eval/de) + 1
540 cpassert((iv > 0) .AND. (iv <= nhist))
541 hist(iv, ispin) = hist(iv, ispin) + wkp
542 occval(iv, ispin) = occval(iv, ispin) + wkp*occupation_numbers(i)
543 END DO
544 END IF
545 END DO
546 END DO
547 END IF
548 CALL para_env%sum(hist)
549 CALL para_env%sum(occval)
550 IF (.NOT. do_broaden) THEN
551 DO ispin = 1, nspins
552 hist(:, ispin) = hist(:, ispin)/real(nmo(ispin), kind=dp)
553 END DO
554 END IF
555 DO i = 1, nhist
556 ehist(i, 1:nspins) = emin + (i - 1)*de
557 END DO
558
559 resolved_energy_zero = dos_resolve_energy_zero(energy_zero, dft_control%smear, fractional_occupation)
560 SELECT CASE (resolved_energy_zero)
561 CASE (dos_energy_zero_absolute)
562 energy_ref(:) = 0.0_dp
563 CASE (dos_energy_zero_hoco)
564 energy_ref(:) = maxval(hoco(1:nspins))
565 CASE DEFAULT
566 energy_ref(:) = maxval(e_fermi(1:nspins))
567 END SELECT
568 IF (.NOT. do_broaden) energy_ref(:) = 0.0_dp
569 energy_factor = merge(dos_energy_scale(energy_unit), 1.0_dp, do_broaden)
570 density_factor = merge(dos_density_scale(energy_unit), 1.0_dp, do_broaden)
571 energy_label = dos_energy_label(merge(energy_unit, 1, do_broaden))
572 zero_label = dos_energy_zero_label(resolved_energy_zero)
573 IF (energy_zero == dos_energy_zero_auto) zero_label = "AUTO -> "//trim(zero_label)
574 ev_factor = dos_energy_scale(dos_energy_unit_ev)
575
576 my_act = "WRITE"
577 IF (do_broaden) THEN
578 iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%DOS", &
579 extension=".dos", file_position=my_pos, file_action=my_act, &
580 file_form="FORMATTED", middle_name="curve")
581 ELSE
582 iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%DOS", &
583 extension=".dos", file_position=my_pos, file_action=my_act, &
584 file_form="FORMATTED")
585 END IF
586 IF (iw > 0) THEN
587 WRITE (unit=iw, fmt="(A,I0)") "# DOS at iteration step i = ", iterstep
588 IF (nspins == 2) THEN
589 WRITE (unit=iw, fmt="(A,2F12.6,A,2F12.6,A)") &
590 "# E(Fermi) = ", e_fermi(1:2), " a.u. = ", e_fermi(1:2)*ev_factor, " eV"
591 WRITE (unit=iw, fmt="(A,2F12.6,A,2F12.6,A)") &
592 "# E(HOCO) = ", hoco(1:2), " a.u. = ", hoco(1:2)*ev_factor, " eV"
593 IF (do_broaden) THEN
594 WRITE (unit=iw, fmt="(A,A)") "# Energy zero: ", trim(zero_label)
595 CALL write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing)
596 END IF
597 WRITE (unit=iw, fmt="(A,A)") "# "//trim(energy_label)//" Alpha_Density Occupation", &
598 " Beta_Density Occupation"
599 ! (F15.8,4F20.ndigits)
600 WRITE (unit=fmtstr_data, fmt="(A,I0,A)") "(F15.8,4F20.", ndigits, ")"
601 ELSE
602 WRITE (unit=iw, fmt="(A,F12.6,A,F12.6,A)") &
603 "# E(Fermi) = ", e_fermi(1), " a.u. = ", e_fermi(1)*ev_factor, " eV"
604 WRITE (unit=iw, fmt="(A,F12.6,A,F12.6,A)") &
605 "# E(HOCO) = ", hoco(1), " a.u. = ", hoco(1)*ev_factor, " eV"
606 IF (do_broaden) THEN
607 WRITE (unit=iw, fmt="(A,A)") "# Energy zero: ", trim(zero_label)
608 CALL write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing)
609 END IF
610 WRITE (unit=iw, fmt="(A,A)") "# "//trim(energy_label), " Density Occupation"
611 ! (F15.8,2F20.ndigits)
612 WRITE (unit=fmtstr_data, fmt="(A,I0,A)") "(F15.8,2F20.", ndigits, ")"
613 END IF
614 DO i = 1, nhist
615 eval = (ehist(i, 1) - energy_ref(1))*energy_factor
616 IF (nspins == 2) THEN
617 ! fmtstr_data == "(F15.8,4F20.xx)"
618 IF (do_broaden) THEN
619 WRITE (unit=iw, fmt=fmtstr_data) eval, hist(i, 1)*density_factor, &
620 occval(i, 1)*density_factor, hist(i, 2)*density_factor, occval(i, 2)*density_factor
621 ELSE
622 IF (hist(i, 1) == 0.0_dp .AND. occval(i, 1) == 0.0_dp .AND. &
623 hist(i, 2) == 0.0_dp .AND. occval(i, 2) == 0.0_dp) cycle
624 WRITE (unit=iw, fmt=fmtstr_data) eval, hist(i, 1), occval(i, 1), &
625 hist(i, 2), occval(i, 2)
626 END IF
627 ELSE
628 ! fmtstr_data == "(F15.8,2F20.xx)"
629 IF (do_broaden) THEN
630 out_density = hist(i, 1)*density_factor
631 out_occ = occval(i, 1)*density_factor
632 ELSE
633 out_density = hist(i, 1)
634 out_occ = occval(i, 1)
635 IF (out_density == 0.0_dp .AND. out_occ == 0.0_dp) cycle
636 END IF
637 WRITE (unit=iw, fmt=fmtstr_data) eval, out_density, out_occ
638 END IF
639 END DO
640 END IF
641 CALL cp_print_key_finished_output(iw, logger, dft_section, "PRINT%DOS")
642 DEALLOCATE (hist, occval, ehist)
643
644 ! destroy the extra k-point set if it was created
645 IF (explicit) THEN
646 CALL kpoint_release(kpoints)
647 END IF
648
649 CALL timestop(handle)
650
651 END SUBROUTINE calculate_dos_kp
652
653END 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:523
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:1073
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:1052
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:288
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:31
qs_dos_utils::dos_energy_zero_auto
integer, parameter, public dos_energy_zero_auto
Definition qs_dos_utils.F:31
qs_dos_utils::dos_energy_unit_ev
integer, parameter, public dos_energy_unit_ev
Definition qs_dos_utils.F:29
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:81
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:125
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:48
qs_dos_utils::dos_energy_zero_absolute
integer, parameter, public dos_energy_zero_absolute
Definition qs_dos_utils.F:31
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:67
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:100
qs_dos
Calculation and writing of density of states.
Definition qs_dos.F:14
qs_dos::calculate_dos_kp
subroutine, public calculate_dos_kp(qs_env, dft_section, write_curve_output)
Compute and write density of states (kpoints)
Definition qs_dos.F:364
qs_dos::calculate_dos
subroutine, public calculate_dos(mos, dft_section, unoccupied_evals, smearing_enabled, write_curve_output)
Compute and write density of states.
Definition qs_dos.F:67
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, matrix_vhxc, 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:534
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:757
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:175
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:750
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
qs_mo_types::mo_set_type
Definition qs_mo_types.F:47