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 : !> \param write_curve_output ...
61 : !> \date 26.02.2008
62 : !> \par History:
63 : !> \author JGH
64 : !> \version 1.0
65 : ! **************************************************************************************************
66 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 66 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ehist, hist, occval
89 66 : 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 66 : NULLIFY (logger)
94 132 : 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 66 : "PRINT%DOS"), cp_p_file)
98 66 : iounit = cp_logger_get_default_io_unit(logger)
99 66 : IF ((.NOT. should_output)) RETURN
100 :
101 66 : CALL timeset(routineN, handle)
102 66 : iterstep = logger%iter_info%iteration(logger%iter_info%n_rlevel)
103 :
104 66 : IF (iounit > 0) WRITE (UNIT=iounit, FMT='(/,(T3,A,T61,I10))') &
105 33 : " Calculate DOS at iteration step ", iterstep
106 :
107 66 : CALL section_vals_val_get(dft_section, "PRINT%DOS%DELTA_E", r_val=de)
108 66 : CALL section_vals_val_get(dft_section, "PRINT%DOS%APPEND", l_val=append)
109 66 : CALL section_vals_val_get(dft_section, "PRINT%DOS%NDIGITS", i_val=ndigits)
110 66 : CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%ENERGY_UNIT", i_val=energy_unit)
111 66 : CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%ENERGY_ZERO", i_val=energy_zero)
112 66 : CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%BROADEN%TYPE", i_val=broaden_type)
113 66 : CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%BROADEN%WIDTH", r_val=broaden_width)
114 66 : CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%BROADEN%VOIGT_MIXING", r_val=voigt_mixing)
115 66 : IF (append .AND. iterstep > 1) THEN
116 10 : my_pos = "APPEND"
117 : ELSE
118 56 : my_pos = "REWIND"
119 : END IF
120 66 : ndigits = MIN(MAX(ndigits, 1), 10)
121 66 : IF (PRESENT(write_curve_output)) THEN
122 0 : IF (write_curve_output .AND. de <= 0.0_dp) THEN
123 0 : CPWARN("Broadened DOS output requires DELTA_E > 0 and will be skipped")
124 0 : CALL timestop(handle)
125 0 : RETURN
126 : END IF
127 0 : IF (write_curve_output .AND. broaden_width <= 0.0_dp) THEN
128 0 : CPWARN("Broadened DOS output requires a finite WIDTH and will be skipped")
129 0 : CALL timestop(handle)
130 0 : RETURN
131 : END IF
132 : END IF
133 : do_broaden = .FALSE.
134 : IF (PRESENT(write_curve_output)) do_broaden = write_curve_output
135 66 : do_broaden = do_broaden .AND. (broaden_width > 0.0_dp)
136 0 : IF (do_broaden) de = MAX(de, 0.00001_dp)
137 :
138 66 : emin = 1.e10_dp
139 66 : emax = -1.e10_dp
140 66 : nspins = SIZE(mos)
141 66 : nmo(:) = 0
142 66 : nvirt(:) = 0
143 66 : nstates(:) = 0
144 198 : hoco(:) = -HUGE(0.0_dp)
145 66 : fractional_occupation = .FALSE.
146 66 : smear_on = .FALSE.
147 66 : IF (PRESENT(smearing_enabled)) smear_on = smearing_enabled
148 :
149 132 : DO ispin = 1, nspins
150 66 : mo_set => mos(ispin)
151 66 : CALL get_mo_set(mo_set=mo_set, nmo=nmo(ispin), mu=e_fermi(ispin))
152 66 : eigenvalues => mo_set%eigenvalues
153 66 : occupation_numbers => mo_set%occupation_numbers
154 642 : DO i = 1, nmo(ispin)
155 576 : IF (occupation_numbers(i) > 1.0e-10_dp) hoco(ispin) = MAX(hoco(ispin), eigenvalues(i))
156 576 : IF (ABS(occupation_numbers(i) - REAL(NINT(occupation_numbers(i)), KIND=dp)) > &
157 74 : 1.0e-8_dp) fractional_occupation = .TRUE.
158 : END DO
159 66 : IF (hoco(ispin) < -0.5_dp*HUGE(0.0_dp)) hoco(ispin) = e_fermi(ispin)
160 66 : IF (PRESENT(unoccupied_evals)) THEN
161 2 : IF (ASSOCIATED(unoccupied_evals(ispin)%array)) nvirt(ispin) = SIZE(unoccupied_evals(ispin)%array)
162 : END IF
163 66 : nstates(ispin) = nmo(ispin) + nvirt(ispin)
164 642 : e1 = MINVAL(eigenvalues(1:nmo(ispin)))
165 642 : e2 = MAXVAL(eigenvalues(1:nmo(ispin)))
166 66 : IF (nvirt(ispin) > 0) THEN
167 12 : e1 = MIN(e1, MINVAL(unoccupied_evals(ispin)%array(1:nvirt(ispin))))
168 12 : e2 = MAX(e2, MAXVAL(unoccupied_evals(ispin)%array(1:nvirt(ispin))))
169 : END IF
170 66 : emin = MIN(emin, e1)
171 132 : emax = MAX(emax, e2)
172 : END DO
173 :
174 66 : IF (do_broaden) THEN
175 0 : broaden_cutoff = broadening_cutoff(broaden_type, broaden_width)
176 0 : emin = emin - broaden_cutoff
177 0 : emax = emax + broaden_cutoff
178 0 : nhist = NINT((emax - emin)/de) + 1
179 0 : ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
180 0 : hist = 0.0_dp
181 0 : occval = 0.0_dp
182 0 : ehist = 0.0_dp
183 0 : DO ispin = 1, nspins
184 0 : mo_set => mos(ispin)
185 0 : occupation_numbers => mo_set%occupation_numbers
186 0 : eigenvalues => mo_set%eigenvalues
187 0 : 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 0 : voigt_mixing)
191 : END DO
192 0 : 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 0 : broaden_type, broaden_width, voigt_mixing)
196 : END DO
197 : END DO
198 0 : DO i = 1, nhist
199 0 : ehist(i, 1:nspins) = emin + (i - 1)*de
200 : END DO
201 66 : ELSE IF (de > 0.0_dp) THEN
202 66 : nhist = NINT((emax - emin)/de) + 1
203 528 : ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
204 66 : hist = 0.0_dp
205 66 : occval = 0.0_dp
206 66 : ehist = 0.0_dp
207 132 : DO ispin = 1, nspins
208 66 : mo_set => mos(ispin)
209 66 : occupation_numbers => mo_set%occupation_numbers
210 66 : eigenvalues => mo_set%eigenvalues
211 642 : DO i = 1, nmo(ispin)
212 576 : eval = eigenvalues(i) - emin
213 576 : iv = NINT(eval/de) + 1
214 576 : CPASSERT((iv > 0) .AND. (iv <= nhist))
215 576 : hist(iv, ispin) = hist(iv, ispin) + 1.0_dp
216 642 : occval(iv, ispin) = occval(iv, ispin) + occupation_numbers(i)
217 : END DO
218 76 : DO i = 1, nvirt(ispin)
219 10 : eval = unoccupied_evals(ispin)%array(i) - emin
220 10 : iv = NINT(eval/de) + 1
221 10 : CPASSERT((iv > 0) .AND. (iv <= nhist))
222 76 : hist(iv, ispin) = hist(iv, ispin) + 1.0_dp
223 : END DO
224 73352 : hist(:, ispin) = hist(:, ispin)/REAL(nstates(ispin), KIND=dp)
225 : END DO
226 73286 : DO i = 1, nhist
227 146506 : ehist(i, 1:nspins) = emin + (i - 1)*de
228 : END DO
229 : ELSE
230 0 : nhist = MAXVAL(nstates)
231 0 : ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
232 0 : hist = 0.0_dp
233 0 : occval = 0.0_dp
234 0 : ehist = 0.0_dp
235 0 : DO ispin = 1, nspins
236 0 : mo_set => mos(ispin)
237 0 : occupation_numbers => mo_set%occupation_numbers
238 0 : eigenvalues => mo_set%eigenvalues
239 0 : DO i = 1, nmo(ispin)
240 0 : ehist(i, ispin) = eigenvalues(i)
241 0 : hist(i, ispin) = 1.0_dp
242 0 : occval(i, ispin) = occupation_numbers(i)
243 : END DO
244 0 : DO i = 1, nvirt(ispin)
245 0 : ehist(nmo(ispin) + i, ispin) = unoccupied_evals(ispin)%array(i)
246 0 : hist(nmo(ispin) + i, ispin) = 1.0_dp
247 : END DO
248 0 : hist(:, ispin) = hist(:, ispin)/REAL(nstates(ispin), KIND=dp)
249 : END DO
250 : END IF
251 :
252 66 : resolved_energy_zero = dos_resolve_energy_zero(energy_zero, smear_on, fractional_occupation)
253 0 : SELECT CASE (resolved_energy_zero)
254 : CASE (dos_energy_zero_absolute)
255 0 : energy_ref(:) = 0.0_dp
256 : CASE (dos_energy_zero_hoco)
257 248 : energy_ref(:) = MAXVAL(hoco(1:nspins))
258 : CASE DEFAULT
259 78 : energy_ref(:) = MAXVAL(e_fermi(1:nspins))
260 : END SELECT
261 66 : IF (.NOT. do_broaden) energy_ref(:) = 0.0_dp
262 66 : energy_factor = MERGE(dos_energy_scale(energy_unit), 1.0_dp, do_broaden)
263 66 : density_factor = MERGE(dos_density_scale(energy_unit), 1.0_dp, do_broaden)
264 66 : energy_label = dos_energy_label(MERGE(energy_unit, 1, do_broaden))
265 66 : zero_label = dos_energy_zero_label(resolved_energy_zero)
266 66 : IF (energy_zero == dos_energy_zero_auto) zero_label = "AUTO -> "//TRIM(zero_label)
267 66 : ev_factor = dos_energy_scale(dos_energy_unit_ev)
268 :
269 66 : my_act = "WRITE"
270 66 : 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 0 : 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 66 : file_form="FORMATTED")
278 : END IF
279 66 : IF (iw > 0) THEN
280 33 : WRITE (UNIT=iw, FMT="(A,I0)") "# DOS at iteration step i = ", iterstep
281 33 : IF (nspins == 2) THEN
282 : WRITE (UNIT=iw, FMT="(A,2F12.6,A,2F12.6,A)") &
283 0 : "# 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 0 : "# E(HOCO) = ", hoco(1:2), " a.u. = ", hoco(1:2)*ev_factor, " eV"
286 0 : IF (do_broaden) THEN
287 0 : WRITE (UNIT=iw, FMT="(A,A)") "# Energy zero: ", TRIM(zero_label)
288 0 : CALL write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing)
289 : END IF
290 0 : IF (do_broaden .OR. de > 0.0_dp) THEN
291 0 : WRITE (UNIT=iw, FMT="(A,A,A)") "# "//TRIM(energy_label)//" Alpha_Density Occupation", &
292 0 : " Beta_Density Occupation"
293 0 : WRITE (UNIT=fmtstr_data, FMT="(A,I0,A)") "(F15.8,4F20.", ndigits, ")"
294 : ELSE
295 0 : WRITE (UNIT=iw, FMT="(A,A,A)") "# "//TRIM(energy_label)//" Alpha_Density Occupation", &
296 0 : " "//TRIM(energy_label), " Beta_Density Occupation"
297 0 : 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 33 : "# 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 33 : "# E(HOCO) = ", hoco(1), " a.u. = ", hoco(1)*ev_factor, " eV"
304 33 : IF (do_broaden) THEN
305 0 : WRITE (UNIT=iw, FMT="(A,A)") "# Energy zero: ", TRIM(zero_label)
306 0 : CALL write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing)
307 : END IF
308 33 : WRITE (UNIT=iw, FMT="(A,A)") "# "//TRIM(energy_label), " Density Occupation"
309 : ! (F15.8,2F20.ndigits)
310 33 : WRITE (UNIT=fmtstr_data, FMT="(A,I0,A)") "(F15.8,2F20.", ndigits, ")"
311 : END IF
312 36643 : DO i = 1, nhist
313 36643 : IF (nspins == 2) THEN
314 0 : IF (do_broaden) THEN
315 0 : eval = (ehist(i, 1) - energy_ref(1))*energy_factor
316 0 : WRITE (UNIT=iw, FMT=fmtstr_data) eval, hist(i, 1)*density_factor, &
317 0 : occval(i, 1)*density_factor, hist(i, 2)*density_factor, &
318 0 : occval(i, 2)*density_factor
319 0 : ELSE IF (de > 0.0_dp) THEN
320 : IF (hist(i, 1) == 0.0_dp .AND. occval(i, 1) == 0.0_dp .AND. &
321 0 : hist(i, 2) == 0.0_dp .AND. occval(i, 2) == 0.0_dp) CYCLE
322 0 : eval = ehist(i, 1)
323 0 : WRITE (UNIT=iw, FMT=fmtstr_data) eval, hist(i, 1), occval(i, 1), &
324 0 : hist(i, 2), occval(i, 2)
325 : ELSE
326 0 : e1 = ehist(i, 1)
327 0 : e2 = ehist(i, 2)
328 0 : WRITE (UNIT=iw, FMT=fmtstr_data) e1, hist(i, 1), occval(i, 1), &
329 0 : e2, hist(i, 2), occval(i, 2)
330 : END IF
331 : ELSE
332 36610 : eval = (ehist(i, 1) - energy_ref(1))*energy_factor
333 : ! fmtstr_data == "(F15.8,2F20.xx)"
334 36610 : IF (do_broaden) THEN
335 0 : out_density = hist(i, 1)*density_factor
336 0 : out_occ = occval(i, 1)*density_factor
337 : ELSE
338 36610 : out_density = hist(i, 1)
339 36610 : out_occ = occval(i, 1)
340 36610 : IF (out_density == 0.0_dp .AND. out_occ == 0.0_dp) CYCLE
341 : END IF
342 271 : WRITE (UNIT=iw, FMT=fmtstr_data) eval, out_density, out_occ
343 : END IF
344 : END DO
345 : END IF
346 66 : CALL cp_print_key_finished_output(iw, logger, dft_section, "PRINT%DOS")
347 66 : DEALLOCATE (hist, occval, ehist)
348 :
349 66 : CALL timestop(handle)
350 :
351 66 : 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 20 : 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 20 : 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 20 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ehist, hist, occval
385 20 : 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 20 : 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 20 : NULLIFY (logger, kpoints)
394 40 : 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 20 : "PRINT%DOS"), cp_p_file)
398 20 : iounit = cp_logger_get_default_io_unit(logger)
399 20 : IF ((.NOT. should_output)) RETURN
400 :
401 20 : CALL timeset(routineN, handle)
402 20 : 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 20 : CALL section_vals_val_get(dft_section, "PRINT%DOS%MP_GRID", i_vals=nkp_grid, explicit=explicit)
406 :
407 20 : IF (explicit) THEN
408 : ! make sure is a valid grid
409 0 : DO i = 1, 3
410 0 : IF (nkp_grid(i) < 1) THEN
411 : WRITE (UNIT=err, FMT='(T4,A,I3,A,I1)') &
412 0 : "Invalid kpoint grid for DOS ", nkp_grid(i), " in dimension ", i
413 0 : CPABORT(TRIM(err))
414 : END IF
415 : END DO
416 : ! calculate orbitals and energies
417 0 : CALL calculate_kp_orbitals(qs_env, kpoints, "MONKHORST-PACK", 0, nkp_grid)
418 : ELSE
419 : ! use the kpoints from the environment
420 20 : CALL get_qs_env(qs_env, kpoints=kpoints)
421 : END IF
422 :
423 20 : IF (iounit > 0) WRITE (UNIT=iounit, FMT='(/,(T3,A,T61,I10))') &
424 10 : " Calculate DOS at iteration step ", iterstep
425 : IF (iounit > 0) WRITE (UNIT=iounit, FMT='((T3,A,3I3,A))') &
426 40 : " Using a", kpoints%nkp_grid(:), ' '//TRIM(kpoints%kp_scheme)//' grid'
427 :
428 20 : CALL section_vals_val_get(dft_section, "PRINT%DOS%DELTA_E", r_val=de)
429 20 : CALL section_vals_val_get(dft_section, "PRINT%DOS%APPEND", l_val=append)
430 20 : CALL section_vals_val_get(dft_section, "PRINT%DOS%NDIGITS", i_val=ndigits)
431 20 : CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%ENERGY_UNIT", i_val=energy_unit)
432 20 : CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%ENERGY_ZERO", i_val=energy_zero)
433 20 : CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%BROADEN%TYPE", i_val=broaden_type)
434 20 : CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%BROADEN%WIDTH", r_val=broaden_width)
435 20 : CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%BROADEN%VOIGT_MIXING", r_val=voigt_mixing)
436 20 : IF (append .AND. iterstep > 1) THEN
437 0 : my_pos = "APPEND"
438 : ELSE
439 20 : my_pos = "REWIND"
440 : END IF
441 20 : ndigits = MIN(MAX(ndigits, 1), 10)
442 20 : IF (PRESENT(write_curve_output)) THEN
443 0 : IF (write_curve_output .AND. de <= 0.0_dp) THEN
444 0 : CPWARN("Broadened k-point DOS output requires DELTA_E > 0 and will be skipped")
445 0 : CALL timestop(handle)
446 0 : IF (explicit) CALL kpoint_release(kpoints)
447 0 : RETURN
448 : END IF
449 0 : IF (write_curve_output .AND. broaden_width <= 0.0_dp) THEN
450 0 : CPWARN("Broadened k-point DOS output requires a finite WIDTH and will be skipped")
451 0 : CALL timestop(handle)
452 0 : IF (explicit) CALL kpoint_release(kpoints)
453 0 : RETURN
454 : END IF
455 20 : ELSE IF (de <= 0.0_dp) THEN
456 0 : CPWARN("K-point DOS output requires DELTA_E > 0 and will be skipped")
457 0 : CALL timestop(handle)
458 0 : IF (explicit) CALL kpoint_release(kpoints)
459 0 : RETURN
460 : END IF
461 : ! ensure a lower value for the DOS grid width
462 20 : de = MAX(de, 0.00001_dp)
463 20 : do_broaden = .FALSE.
464 20 : IF (PRESENT(write_curve_output)) do_broaden = write_curve_output
465 0 : do_broaden = do_broaden .AND. (broaden_width > 0.0_dp)
466 :
467 20 : CALL get_qs_env(qs_env, dft_control=dft_control)
468 20 : nspins = dft_control%nspins
469 20 : para_env => kpoints%para_env_inter_kp
470 :
471 20 : emin = 1.e10_dp
472 20 : emax = -1.e10_dp
473 20 : nmo(:) = 0
474 20 : e_fermi(:) = 0.0_dp
475 60 : hoco(:) = -HUGE(0.0_dp)
476 20 : fractional_occupation = .FALSE.
477 20 : IF (kpoints%nkp /= 0) THEN
478 132 : DO ik = 1, SIZE(kpoints%kp_env)
479 112 : mos => kpoints%kp_env(ik)%kpoint_env%mos
480 112 : CPASSERT(ASSOCIATED(mos))
481 244 : DO ispin = 1, nspins
482 112 : mo_set => mos(1, ispin)
483 112 : CALL get_mo_set(mo_set=mo_set, nmo=nmo_kp, mu=e_fermi(ispin))
484 112 : eigenvalues => mo_set%eigenvalues
485 112 : occupation_numbers => mo_set%occupation_numbers
486 2334 : DO i = 1, nmo_kp
487 2222 : IF (occupation_numbers(i) > 1.0e-10_dp) hoco(ispin) = MAX(hoco(ispin), eigenvalues(i))
488 2222 : IF (ABS(occupation_numbers(i) - REAL(NINT(occupation_numbers(i)), KIND=dp)) > &
489 314 : 1.0e-8_dp) fractional_occupation = .TRUE.
490 : END DO
491 2334 : e1 = MINVAL(eigenvalues(1:nmo_kp))
492 2334 : e2 = MAXVAL(eigenvalues(1:nmo_kp))
493 112 : emin = MIN(emin, e1)
494 112 : emax = MAX(emax, e2)
495 336 : nmo(ispin) = MAX(nmo(ispin), nmo_kp)
496 : END DO
497 : END DO
498 : END IF
499 20 : CALL para_env%min(emin)
500 20 : CALL para_env%max(emax)
501 20 : CALL para_env%max(nmo)
502 20 : CALL para_env%max(e_fermi)
503 20 : CALL para_env%max(hoco)
504 20 : fractional_occupation_int = MERGE(1, 0, fractional_occupation)
505 20 : CALL para_env%max(fractional_occupation_int)
506 20 : fractional_occupation = (fractional_occupation_int /= 0)
507 40 : DO ispin = 1, nspins
508 40 : IF (hoco(ispin) < -0.5_dp*HUGE(0.0_dp)) hoco(ispin) = e_fermi(ispin)
509 : END DO
510 :
511 20 : IF (do_broaden) THEN
512 0 : broaden_cutoff = broadening_cutoff(broaden_type, broaden_width)
513 0 : emin = emin - broaden_cutoff
514 0 : emax = emax + broaden_cutoff
515 : END IF
516 20 : nhist = NINT((emax - emin)/de) + 1
517 160 : ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
518 20 : hist = 0.0_dp
519 20 : occval = 0.0_dp
520 20 : ehist = 0.0_dp
521 :
522 20 : IF (kpoints%nkp /= 0) THEN
523 132 : DO ik = 1, SIZE(kpoints%kp_env)
524 112 : mos => kpoints%kp_env(ik)%kpoint_env%mos
525 112 : wkp = kpoints%kp_env(ik)%kpoint_env%wkp
526 244 : DO ispin = 1, nspins
527 112 : mo_set => mos(1, ispin)
528 112 : occupation_numbers => mo_set%occupation_numbers
529 112 : eigenvalues => mo_set%eigenvalues
530 224 : IF (do_broaden) THEN
531 0 : 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 0 : voigt_mixing)
535 : END DO
536 : ELSE
537 2334 : DO i = 1, nmo(ispin)
538 2222 : eval = eigenvalues(i) - emin
539 2222 : iv = NINT(eval/de) + 1
540 2222 : CPASSERT((iv > 0) .AND. (iv <= nhist))
541 2222 : hist(iv, ispin) = hist(iv, ispin) + wkp
542 2334 : 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 20 : CALL para_env%sum(hist)
549 20 : CALL para_env%sum(occval)
550 20 : IF (.NOT. do_broaden) THEN
551 40 : DO ispin = 1, nspins
552 137324 : hist(:, ispin) = hist(:, ispin)/REAL(nmo(ispin), KIND=dp)
553 : END DO
554 : END IF
555 137304 : DO i = 1, nhist
556 274588 : ehist(i, 1:nspins) = emin + (i - 1)*de
557 : END DO
558 :
559 20 : resolved_energy_zero = dos_resolve_energy_zero(energy_zero, dft_control%smear, fractional_occupation)
560 0 : SELECT CASE (resolved_energy_zero)
561 : CASE (dos_energy_zero_absolute)
562 0 : energy_ref(:) = 0.0_dp
563 : CASE (dos_energy_zero_hoco)
564 16 : energy_ref(:) = MAXVAL(hoco(1:nspins))
565 : CASE DEFAULT
566 68 : energy_ref(:) = MAXVAL(e_fermi(1:nspins))
567 : END SELECT
568 20 : IF (.NOT. do_broaden) energy_ref(:) = 0.0_dp
569 20 : energy_factor = MERGE(dos_energy_scale(energy_unit), 1.0_dp, do_broaden)
570 20 : density_factor = MERGE(dos_density_scale(energy_unit), 1.0_dp, do_broaden)
571 20 : energy_label = dos_energy_label(MERGE(energy_unit, 1, do_broaden))
572 20 : zero_label = dos_energy_zero_label(resolved_energy_zero)
573 20 : IF (energy_zero == dos_energy_zero_auto) zero_label = "AUTO -> "//TRIM(zero_label)
574 20 : ev_factor = dos_energy_scale(dos_energy_unit_ev)
575 :
576 20 : my_act = "WRITE"
577 20 : 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 0 : 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 20 : file_form="FORMATTED")
585 : END IF
586 20 : IF (iw > 0) THEN
587 10 : WRITE (UNIT=iw, FMT="(A,I0)") "# DOS at iteration step i = ", iterstep
588 10 : IF (nspins == 2) THEN
589 : WRITE (UNIT=iw, FMT="(A,2F12.6,A,2F12.6,A)") &
590 0 : "# 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 0 : "# E(HOCO) = ", hoco(1:2), " a.u. = ", hoco(1:2)*ev_factor, " eV"
593 0 : IF (do_broaden) THEN
594 0 : WRITE (UNIT=iw, FMT="(A,A)") "# Energy zero: ", TRIM(zero_label)
595 0 : CALL write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing)
596 : END IF
597 0 : WRITE (UNIT=iw, FMT="(A,A)") "# "//TRIM(energy_label)//" Alpha_Density Occupation", &
598 0 : " Beta_Density Occupation"
599 : ! (F15.8,4F20.ndigits)
600 0 : 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 10 : "# 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 10 : "# E(HOCO) = ", hoco(1), " a.u. = ", hoco(1)*ev_factor, " eV"
606 10 : IF (do_broaden) THEN
607 0 : WRITE (UNIT=iw, FMT="(A,A)") "# Energy zero: ", TRIM(zero_label)
608 0 : CALL write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing)
609 : END IF
610 10 : WRITE (UNIT=iw, FMT="(A,A)") "# "//TRIM(energy_label), " Density Occupation"
611 : ! (F15.8,2F20.ndigits)
612 10 : WRITE (UNIT=fmtstr_data, FMT="(A,I0,A)") "(F15.8,2F20.", ndigits, ")"
613 : END IF
614 68652 : DO i = 1, nhist
615 68642 : eval = (ehist(i, 1) - energy_ref(1))*energy_factor
616 68652 : IF (nspins == 2) THEN
617 : ! fmtstr_data == "(F15.8,4F20.xx)"
618 0 : IF (do_broaden) THEN
619 0 : WRITE (UNIT=iw, FMT=fmtstr_data) eval, hist(i, 1)*density_factor, &
620 0 : 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 0 : hist(i, 2) == 0.0_dp .AND. occval(i, 2) == 0.0_dp) CYCLE
624 0 : WRITE (UNIT=iw, FMT=fmtstr_data) eval, hist(i, 1), occval(i, 1), &
625 0 : hist(i, 2), occval(i, 2)
626 : END IF
627 : ELSE
628 : ! fmtstr_data == "(F15.8,2F20.xx)"
629 68642 : IF (do_broaden) THEN
630 0 : out_density = hist(i, 1)*density_factor
631 0 : out_occ = occval(i, 1)*density_factor
632 : ELSE
633 68642 : out_density = hist(i, 1)
634 68642 : out_occ = occval(i, 1)
635 68642 : IF (out_density == 0.0_dp .AND. out_occ == 0.0_dp) CYCLE
636 : END IF
637 813 : WRITE (UNIT=iw, FMT=fmtstr_data) eval, out_density, out_occ
638 : END IF
639 : END DO
640 : END IF
641 20 : CALL cp_print_key_finished_output(iw, logger, dft_section, "PRINT%DOS")
642 20 : DEALLOCATE (hist, occval, ehist)
643 :
644 : ! destroy the extra k-point set if it was created
645 20 : IF (explicit) THEN
646 0 : CALL kpoint_release(kpoints)
647 : END IF
648 :
649 20 : CALL timestop(handle)
650 :
651 80 : END SUBROUTINE calculate_dos_kp
652 :
653 : END MODULE qs_dos
|
