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 Calculation and writing of density of states
10 : !> \par History
11 : !> -
12 : !> \author JGH
13 : ! **************************************************************************************************
14 : MODULE 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 :
52 : CONTAINS
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 66 : 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 66 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ehist, hist, occval
88 66 : 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 66 : NULLIFY (logger)
93 132 : 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 66 : "PRINT%DOS"), cp_p_file)
97 66 : iounit = cp_logger_get_default_io_unit(logger)
98 66 : IF ((.NOT. should_output)) RETURN
99 :
100 66 : CALL timeset(routineN, handle)
101 66 : iterstep = logger%iter_info%iteration(logger%iter_info%n_rlevel)
102 :
103 66 : IF (iounit > 0) WRITE (UNIT=iounit, FMT='(/,(T3,A,T61,I10))') &
104 33 : " Calculate DOS at iteration step ", iterstep
105 :
106 66 : CALL section_vals_val_get(dft_section, "PRINT%DOS%DELTA_E", r_val=de)
107 66 : CALL section_vals_val_get(dft_section, "PRINT%DOS%APPEND", l_val=append)
108 66 : CALL section_vals_val_get(dft_section, "PRINT%DOS%NDIGITS", i_val=ndigits)
109 66 : CALL section_vals_val_get(dft_section, "PRINT%DOS%ENERGY_UNIT", i_val=energy_unit)
110 66 : CALL section_vals_val_get(dft_section, "PRINT%DOS%ENERGY_ZERO", i_val=energy_zero)
111 66 : CALL section_vals_val_get(dft_section, "PRINT%DOS%BROADEN_TYPE", i_val=broaden_type)
112 66 : CALL section_vals_val_get(dft_section, "PRINT%DOS%BROADEN_WIDTH", r_val=broaden_width)
113 66 : CALL section_vals_val_get(dft_section, "PRINT%DOS%VOIGT_MIXING", r_val=voigt_mixing)
114 66 : IF (append .AND. iterstep > 1) THEN
115 10 : my_pos = "APPEND"
116 : ELSE
117 56 : my_pos = "REWIND"
118 : END IF
119 66 : ndigits = MIN(MAX(ndigits, 1), 10)
120 66 : do_broaden = (broaden_width > 0.0_dp)
121 66 : IF (do_broaden) de = MAX(de, 0.00001_dp)
122 :
123 66 : emin = 1.e10_dp
124 66 : emax = -1.e10_dp
125 66 : nspins = SIZE(mos)
126 66 : nmo(:) = 0
127 66 : nvirt(:) = 0
128 66 : nstates(:) = 0
129 198 : hoco(:) = -HUGE(0.0_dp)
130 66 : fractional_occupation = .FALSE.
131 66 : smear_on = .FALSE.
132 66 : IF (PRESENT(smearing_enabled)) smear_on = smearing_enabled
133 :
134 132 : DO ispin = 1, nspins
135 66 : mo_set => mos(ispin)
136 66 : CALL get_mo_set(mo_set=mo_set, nmo=nmo(ispin), mu=e_fermi(ispin))
137 66 : eigenvalues => mo_set%eigenvalues
138 66 : occupation_numbers => mo_set%occupation_numbers
139 642 : DO i = 1, nmo(ispin)
140 576 : IF (occupation_numbers(i) > 1.0e-10_dp) hoco(ispin) = MAX(hoco(ispin), eigenvalues(i))
141 576 : IF (ABS(occupation_numbers(i) - REAL(NINT(occupation_numbers(i)), KIND=dp)) > &
142 74 : 1.0e-8_dp) fractional_occupation = .TRUE.
143 : END DO
144 66 : IF (hoco(ispin) < -0.5_dp*HUGE(0.0_dp)) hoco(ispin) = e_fermi(ispin)
145 66 : IF (PRESENT(unoccupied_evals)) THEN
146 2 : IF (ASSOCIATED(unoccupied_evals(ispin)%array)) nvirt(ispin) = SIZE(unoccupied_evals(ispin)%array)
147 : END IF
148 66 : nstates(ispin) = nmo(ispin) + nvirt(ispin)
149 642 : e1 = MINVAL(eigenvalues(1:nmo(ispin)))
150 642 : e2 = MAXVAL(eigenvalues(1:nmo(ispin)))
151 66 : IF (nvirt(ispin) > 0) THEN
152 12 : e1 = MIN(e1, MINVAL(unoccupied_evals(ispin)%array(1:nvirt(ispin))))
153 12 : e2 = MAX(e2, MAXVAL(unoccupied_evals(ispin)%array(1:nvirt(ispin))))
154 : END IF
155 66 : emin = MIN(emin, e1)
156 132 : emax = MAX(emax, e2)
157 : END DO
158 :
159 66 : IF (do_broaden) THEN
160 66 : broaden_cutoff = broadening_cutoff(broaden_type, broaden_width)
161 66 : emin = emin - broaden_cutoff
162 66 : emax = emax + broaden_cutoff
163 66 : nhist = NINT((emax - emin)/de) + 1
164 528 : ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
165 66 : hist = 0.0_dp
166 66 : occval = 0.0_dp
167 66 : ehist = 0.0_dp
168 132 : DO ispin = 1, nspins
169 66 : mo_set => mos(ispin)
170 66 : occupation_numbers => mo_set%occupation_numbers
171 66 : eigenvalues => mo_set%eigenvalues
172 642 : 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 642 : voigt_mixing)
176 : END DO
177 142 : 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 76 : broaden_type, broaden_width, voigt_mixing)
181 : END DO
182 : END DO
183 76854 : DO i = 1, nhist
184 153642 : ehist(i, 1:nspins) = emin + (i - 1)*de
185 : END DO
186 0 : ELSE IF (de > 0.0_dp) THEN
187 0 : nhist = NINT((emax - emin)/de) + 1
188 0 : ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
189 0 : hist = 0.0_dp
190 0 : occval = 0.0_dp
191 0 : ehist = 0.0_dp
192 0 : DO ispin = 1, nspins
193 0 : mo_set => mos(ispin)
194 0 : occupation_numbers => mo_set%occupation_numbers
195 0 : eigenvalues => mo_set%eigenvalues
196 0 : DO i = 1, nmo(ispin)
197 0 : eval = eigenvalues(i) - emin
198 0 : iv = NINT(eval/de) + 1
199 0 : CPASSERT((iv > 0) .AND. (iv <= nhist))
200 0 : hist(iv, ispin) = hist(iv, ispin) + 1.0_dp
201 0 : occval(iv, ispin) = occval(iv, ispin) + occupation_numbers(i)
202 : END DO
203 0 : DO i = 1, nvirt(ispin)
204 0 : eval = unoccupied_evals(ispin)%array(i) - emin
205 0 : iv = NINT(eval/de) + 1
206 0 : CPASSERT((iv > 0) .AND. (iv <= nhist))
207 0 : hist(iv, ispin) = hist(iv, ispin) + 1.0_dp
208 : END DO
209 0 : hist(:, ispin) = hist(:, ispin)/REAL(nstates(ispin), KIND=dp)
210 : END DO
211 0 : DO i = 1, nhist
212 0 : ehist(i, 1:nspins) = emin + (i - 1)*de
213 : END DO
214 : ELSE
215 0 : nhist = MAXVAL(nstates)
216 0 : ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
217 0 : hist = 0.0_dp
218 0 : occval = 0.0_dp
219 0 : ehist = 0.0_dp
220 0 : DO ispin = 1, nspins
221 0 : mo_set => mos(ispin)
222 0 : occupation_numbers => mo_set%occupation_numbers
223 0 : eigenvalues => mo_set%eigenvalues
224 0 : DO i = 1, nmo(ispin)
225 0 : ehist(i, ispin) = eigenvalues(i)
226 0 : hist(i, ispin) = 1.0_dp
227 0 : occval(i, ispin) = occupation_numbers(i)
228 : END DO
229 0 : DO i = 1, nvirt(ispin)
230 0 : ehist(nmo(ispin) + i, ispin) = unoccupied_evals(ispin)%array(i)
231 0 : hist(nmo(ispin) + i, ispin) = 1.0_dp
232 : END DO
233 0 : hist(:, ispin) = hist(:, ispin)/REAL(nstates(ispin), KIND=dp)
234 : END DO
235 : END IF
236 :
237 66 : resolved_energy_zero = dos_resolve_energy_zero(energy_zero, smear_on, fractional_occupation)
238 0 : SELECT CASE (resolved_energy_zero)
239 : CASE (dos_energy_zero_absolute)
240 0 : energy_ref(:) = 0.0_dp
241 : CASE (dos_energy_zero_hoco)
242 62 : energy_ref(:) = hoco(:)
243 : CASE DEFAULT
244 66 : energy_ref(:) = e_fermi(:)
245 : END SELECT
246 66 : energy_factor = dos_energy_scale(energy_unit)
247 66 : density_factor = dos_density_scale(energy_unit)
248 66 : energy_label = dos_energy_label(energy_unit)
249 66 : zero_label = dos_energy_zero_label(resolved_energy_zero)
250 66 : IF (energy_zero == dos_energy_zero_auto) zero_label = "AUTO -> "//TRIM(zero_label)
251 66 : ev_factor = dos_energy_scale(dos_energy_unit_ev)
252 :
253 66 : 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 66 : file_form="FORMATTED")
257 66 : IF (iw > 0) THEN
258 33 : WRITE (UNIT=iw, FMT="(A,I0)") "# DOS at iteration step i = ", iterstep
259 33 : IF (nspins == 2) THEN
260 : WRITE (UNIT=iw, FMT="(A,2F12.6,A,2F12.6,A)") &
261 0 : "# 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 0 : "# E(HOCO) = ", hoco(1:2), " a.u. = ", hoco(1:2)*ev_factor, " eV"
264 0 : WRITE (UNIT=iw, FMT="(A,A)") "# Energy zero: ", TRIM(zero_label)
265 0 : IF (do_broaden) CALL write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing)
266 0 : WRITE (UNIT=iw, FMT="(T2,A,A,A)") "# "//TRIM(energy_label)//" Alpha_Density Occupation", &
267 0 : " "//TRIM(energy_label), " Beta_Density Occupation"
268 : ! (2(F15.8,2F20.ndigits))
269 0 : 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 33 : "# 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 33 : "# E(HOCO) = ", hoco(1), " a.u. = ", hoco(1)*ev_factor, " eV"
275 33 : WRITE (UNIT=iw, FMT="(A,A)") "# Energy zero: ", TRIM(zero_label)
276 33 : IF (do_broaden) CALL write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing)
277 33 : WRITE (UNIT=iw, FMT="(A,A)") "# "//TRIM(energy_label), " Density Occupation"
278 : ! (F15.8,2F20.ndigits)
279 33 : WRITE (UNIT=fmtstr_data, FMT="(A,I0,A)") "(F15.8,2F20.", ndigits, ")"
280 : END IF
281 38427 : DO i = 1, nhist
282 38427 : IF (nspins == 2) THEN
283 0 : e1 = (ehist(i, 1) - energy_ref(1))*energy_factor
284 0 : e2 = (ehist(i, 2) - energy_ref(2))*energy_factor
285 : ! fmtstr_data == "(2(F15.8,2F20.xx))"
286 0 : IF (do_broaden) THEN
287 0 : WRITE (UNIT=iw, FMT=fmtstr_data) e1, hist(i, 1)*density_factor, &
288 0 : occval(i, 1)*density_factor, e2, hist(i, 2)*density_factor, &
289 0 : occval(i, 2)*density_factor
290 : ELSE
291 0 : WRITE (UNIT=iw, FMT=fmtstr_data) e1, hist(i, 1), occval(i, 1), &
292 0 : e2, hist(i, 2), occval(i, 2)
293 : END IF
294 : ELSE
295 38394 : eval = (ehist(i, 1) - energy_ref(1))*energy_factor
296 : ! fmtstr_data == "(F15.8,2F20.xx)"
297 38394 : IF (do_broaden) THEN
298 38394 : out_density = hist(i, 1)*density_factor
299 38394 : out_occ = occval(i, 1)*density_factor
300 : ELSE
301 0 : out_density = hist(i, 1)
302 0 : out_occ = occval(i, 1)
303 : END IF
304 38394 : WRITE (UNIT=iw, FMT=fmtstr_data) eval, out_density, out_occ
305 : END IF
306 : END DO
307 : END IF
308 66 : CALL cp_print_key_finished_output(iw, logger, dft_section, "PRINT%DOS")
309 66 : DEALLOCATE (hist, occval, ehist)
310 :
311 66 : CALL timestop(handle)
312 :
313 66 : 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 20 : 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 20 : 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 20 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ehist, hist, occval
345 20 : 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 20 : 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 20 : NULLIFY (logger, kpoints)
354 40 : 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 20 : "PRINT%DOS"), cp_p_file)
358 20 : iounit = cp_logger_get_default_io_unit(logger)
359 20 : IF ((.NOT. should_output)) RETURN
360 :
361 20 : CALL timeset(routineN, handle)
362 20 : 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 20 : CALL section_vals_val_get(dft_section, "PRINT%DOS%MP_GRID", i_vals=nkp_grid, explicit=explicit)
366 :
367 20 : IF (explicit) THEN
368 : ! make sure is a valid grid
369 0 : DO i = 1, 3
370 0 : IF (nkp_grid(i) < 1) THEN
371 : WRITE (UNIT=iounit, FMT='(T4,A,I3,A,I1)') &
372 0 : "Invalid kpoint grid for DOS ", nkp_grid(i), " in dimension ", i
373 0 : CPABORT("")
374 : END IF
375 : END DO
376 : ! calculate orbitals and energies
377 0 : CALL calculate_kp_orbitals(qs_env, kpoints, "MONKHORST-PACK", 0, nkp_grid)
378 : ELSE
379 : ! use the kpoints from the environment
380 20 : CALL get_qs_env(qs_env, kpoints=kpoints)
381 : END IF
382 :
383 20 : IF (iounit > 0) WRITE (UNIT=iounit, FMT='(/,(T3,A,T61,I10))') &
384 10 : " Calculate DOS at iteration step ", iterstep
385 : IF (iounit > 0) WRITE (UNIT=iounit, FMT='((T3,A,3I3,A))') &
386 40 : " Using a", kpoints%nkp_grid(:), ' '//TRIM(kpoints%kp_scheme)//' grid'
387 :
388 20 : CALL section_vals_val_get(dft_section, "PRINT%DOS%DELTA_E", r_val=de)
389 20 : CALL section_vals_val_get(dft_section, "PRINT%DOS%APPEND", l_val=append)
390 20 : CALL section_vals_val_get(dft_section, "PRINT%DOS%NDIGITS", i_val=ndigits)
391 20 : CALL section_vals_val_get(dft_section, "PRINT%DOS%ENERGY_UNIT", i_val=energy_unit)
392 20 : CALL section_vals_val_get(dft_section, "PRINT%DOS%ENERGY_ZERO", i_val=energy_zero)
393 20 : CALL section_vals_val_get(dft_section, "PRINT%DOS%BROADEN_TYPE", i_val=broaden_type)
394 20 : CALL section_vals_val_get(dft_section, "PRINT%DOS%BROADEN_WIDTH", r_val=broaden_width)
395 20 : 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 20 : de = MAX(de, 0.00001_dp)
398 20 : IF (append .AND. iterstep > 1) THEN
399 0 : my_pos = "APPEND"
400 : ELSE
401 20 : my_pos = "REWIND"
402 : END IF
403 20 : ndigits = MIN(MAX(ndigits, 1), 10)
404 20 : do_broaden = (broaden_width > 0.0_dp)
405 :
406 20 : CALL get_qs_env(qs_env, dft_control=dft_control)
407 20 : nspins = dft_control%nspins
408 20 : para_env => kpoints%para_env_inter_kp
409 :
410 20 : emin = 1.e10_dp
411 20 : emax = -1.e10_dp
412 20 : nmo(:) = 0
413 20 : e_fermi(:) = 0.0_dp
414 60 : hoco(:) = -HUGE(0.0_dp)
415 20 : fractional_occupation = .FALSE.
416 20 : IF (kpoints%nkp /= 0) THEN
417 132 : DO ik = 1, SIZE(kpoints%kp_env)
418 112 : mos => kpoints%kp_env(ik)%kpoint_env%mos
419 112 : CPASSERT(ASSOCIATED(mos))
420 244 : DO ispin = 1, nspins
421 112 : mo_set => mos(1, ispin)
422 112 : CALL get_mo_set(mo_set=mo_set, nmo=nmo_kp, mu=e_fermi(ispin))
423 112 : eigenvalues => mo_set%eigenvalues
424 112 : occupation_numbers => mo_set%occupation_numbers
425 2334 : DO i = 1, nmo_kp
426 2222 : IF (occupation_numbers(i) > 1.0e-10_dp) hoco(ispin) = MAX(hoco(ispin), eigenvalues(i))
427 2222 : IF (ABS(occupation_numbers(i) - REAL(NINT(occupation_numbers(i)), KIND=dp)) > &
428 314 : 1.0e-8_dp) fractional_occupation = .TRUE.
429 : END DO
430 2334 : e1 = MINVAL(eigenvalues(1:nmo_kp))
431 2334 : e2 = MAXVAL(eigenvalues(1:nmo_kp))
432 112 : emin = MIN(emin, e1)
433 112 : emax = MAX(emax, e2)
434 336 : nmo(ispin) = MAX(nmo(ispin), nmo_kp)
435 : END DO
436 : END DO
437 : END IF
438 20 : CALL para_env%min(emin)
439 20 : CALL para_env%max(emax)
440 20 : CALL para_env%max(nmo)
441 20 : CALL para_env%max(e_fermi)
442 20 : CALL para_env%max(hoco)
443 20 : fractional_occupation_int = MERGE(1, 0, fractional_occupation)
444 20 : CALL para_env%max(fractional_occupation_int)
445 20 : fractional_occupation = (fractional_occupation_int /= 0)
446 40 : DO ispin = 1, nspins
447 40 : IF (hoco(ispin) < -0.5_dp*HUGE(0.0_dp)) hoco(ispin) = e_fermi(ispin)
448 : END DO
449 :
450 20 : IF (do_broaden) THEN
451 20 : broaden_cutoff = broadening_cutoff(broaden_type, broaden_width)
452 20 : emin = emin - broaden_cutoff
453 20 : emax = emax + broaden_cutoff
454 : END IF
455 20 : nhist = NINT((emax - emin)/de) + 1
456 160 : ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
457 20 : hist = 0.0_dp
458 20 : occval = 0.0_dp
459 20 : ehist = 0.0_dp
460 :
461 20 : IF (kpoints%nkp /= 0) THEN
462 132 : DO ik = 1, SIZE(kpoints%kp_env)
463 112 : mos => kpoints%kp_env(ik)%kpoint_env%mos
464 112 : wkp = kpoints%kp_env(ik)%kpoint_env%wkp
465 244 : DO ispin = 1, nspins
466 112 : mo_set => mos(1, ispin)
467 112 : occupation_numbers => mo_set%occupation_numbers
468 112 : eigenvalues => mo_set%eigenvalues
469 224 : IF (do_broaden) THEN
470 2334 : 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 2334 : voigt_mixing)
474 : END DO
475 : ELSE
476 0 : DO i = 1, nmo(ispin)
477 0 : eval = eigenvalues(i) - emin
478 0 : iv = NINT(eval/de) + 1
479 0 : CPASSERT((iv > 0) .AND. (iv <= nhist))
480 0 : hist(iv, ispin) = hist(iv, ispin) + wkp
481 0 : 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 20 : CALL para_env%sum(hist)
488 20 : CALL para_env%sum(occval)
489 20 : IF (.NOT. do_broaden) THEN
490 0 : DO ispin = 1, nspins
491 0 : hist(:, ispin) = hist(:, ispin)/REAL(nmo(ispin), KIND=dp)
492 : END DO
493 : END IF
494 137738 : DO i = 1, nhist
495 275456 : ehist(i, 1:nspins) = emin + (i - 1)*de
496 : END DO
497 :
498 20 : resolved_energy_zero = dos_resolve_energy_zero(energy_zero, dft_control%smear, fractional_occupation)
499 0 : SELECT CASE (resolved_energy_zero)
500 : CASE (dos_energy_zero_absolute)
501 0 : energy_ref(:) = 0.0_dp
502 : CASE (dos_energy_zero_hoco)
503 4 : energy_ref(:) = hoco(:)
504 : CASE DEFAULT
505 20 : energy_ref(:) = e_fermi(:)
506 : END SELECT
507 20 : energy_factor = dos_energy_scale(energy_unit)
508 20 : density_factor = dos_density_scale(energy_unit)
509 20 : energy_label = dos_energy_label(energy_unit)
510 20 : zero_label = dos_energy_zero_label(resolved_energy_zero)
511 20 : IF (energy_zero == dos_energy_zero_auto) zero_label = "AUTO -> "//TRIM(zero_label)
512 20 : ev_factor = dos_energy_scale(dos_energy_unit_ev)
513 :
514 20 : 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 20 : file_form="FORMATTED")
518 20 : IF (iw > 0) THEN
519 10 : WRITE (UNIT=iw, FMT="(A,I0)") "# DOS at iteration step i = ", iterstep
520 10 : IF (nspins == 2) THEN
521 : WRITE (UNIT=iw, FMT="(A,2F12.6,A,2F12.6,A)") &
522 0 : "# 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 0 : "# E(HOCO) = ", hoco(1:2), " a.u. = ", hoco(1:2)*ev_factor, " eV"
525 0 : WRITE (UNIT=iw, FMT="(A,A)") "# Energy zero: ", TRIM(zero_label)
526 0 : IF (do_broaden) CALL write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing)
527 0 : WRITE (UNIT=iw, FMT="(A,A)") "# "//TRIM(energy_label)//" Alpha_Density Occupation", &
528 0 : " Beta_Density Occupation"
529 : ! (F15.8,4F20.ndigits)
530 0 : 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 10 : "# 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 10 : "# E(HOCO) = ", hoco(1), " a.u. = ", hoco(1)*ev_factor, " eV"
536 10 : WRITE (UNIT=iw, FMT="(A,A)") "# Energy zero: ", TRIM(zero_label)
537 10 : IF (do_broaden) CALL write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing)
538 10 : WRITE (UNIT=iw, FMT="(A,A)") "# "//TRIM(energy_label), " Density Occupation"
539 : ! (F15.8,2F20.ndigits)
540 10 : WRITE (UNIT=fmtstr_data, FMT="(A,I0,A)") "(F15.8,2F20.", ndigits, ")"
541 : END IF
542 68869 : DO i = 1, nhist
543 68859 : eval = (ehist(i, 1) - energy_ref(1))*energy_factor
544 68869 : IF (nspins == 2) THEN
545 : ! fmtstr_data == "(F15.8,4F20.xx)"
546 0 : IF (do_broaden) THEN
547 0 : WRITE (UNIT=iw, FMT=fmtstr_data) eval, hist(i, 1)*density_factor, &
548 0 : occval(i, 1)*density_factor, hist(i, 2)*density_factor, occval(i, 2)*density_factor
549 : ELSE
550 0 : WRITE (UNIT=iw, FMT=fmtstr_data) eval, hist(i, 1), occval(i, 1), &
551 0 : hist(i, 2), occval(i, 2)
552 : END IF
553 : ELSE
554 : ! fmtstr_data == "(F15.8,2F20.xx)"
555 68859 : IF (do_broaden) THEN
556 68859 : out_density = hist(i, 1)*density_factor
557 68859 : out_occ = occval(i, 1)*density_factor
558 : ELSE
559 0 : out_density = hist(i, 1)
560 0 : out_occ = occval(i, 1)
561 : END IF
562 68859 : WRITE (UNIT=iw, FMT=fmtstr_data) eval, out_density, out_occ
563 : END IF
564 : END DO
565 : END IF
566 20 : CALL cp_print_key_finished_output(iw, logger, dft_section, "PRINT%DOS")
567 20 : DEALLOCATE (hist, occval, ehist)
568 :
569 : ! destroy the extra k-point set if it was created
570 20 : IF (explicit) THEN
571 0 : CALL kpoint_release(kpoints)
572 : END IF
573 :
574 20 : CALL timestop(handle)
575 :
576 80 : END SUBROUTINE calculate_dos_kp
577 :
578 : END MODULE qs_dos
|
