Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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_control_types, ONLY: dft_control_type
16 : USE cp_log_handling, ONLY: cp_get_default_logger,&
17 : cp_logger_get_default_io_unit,&
18 : cp_logger_type
19 : USE cp_output_handling, ONLY: cp_p_file,&
20 : cp_print_key_finished_output,&
21 : cp_print_key_should_output,&
22 : cp_print_key_unit_nr
23 : USE input_section_types, ONLY: section_vals_type,&
24 : section_vals_val_get
25 : USE kinds, ONLY: default_string_length,&
26 : dp
27 : USE kpoint_types, ONLY: kpoint_release,&
28 : kpoint_type
29 : USE message_passing, ONLY: mp_para_env_type
30 : USE qs_band_structure, ONLY: calculate_kp_orbitals
31 : USE qs_environment_types, ONLY: get_qs_env,&
32 : qs_environment_type
33 : USE qs_mo_types, ONLY: get_mo_set,&
34 : mo_set_type
35 : #include "./base/base_uses.f90"
36 :
37 : IMPLICIT NONE
38 :
39 : PRIVATE
40 :
41 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dos'
42 :
43 : PUBLIC :: calculate_dos, calculate_dos_kp
44 :
45 : ! **************************************************************************************************
46 :
47 : CONTAINS
48 :
49 : ! **************************************************************************************************
50 : !> \brief Compute and write density of states
51 : !> \param mos ...
52 : !> \param dft_section ...
53 : !> \date 26.02.2008
54 : !> \par History:
55 : !> \author JGH
56 : !> \version 1.0
57 : ! **************************************************************************************************
58 40 : SUBROUTINE calculate_dos(mos, dft_section)
59 :
60 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
61 : TYPE(section_vals_type), POINTER :: dft_section
62 :
63 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_dos'
64 :
65 : CHARACTER(LEN=20) :: fmtstr_data
66 : CHARACTER(LEN=default_string_length) :: my_act, my_pos
67 : INTEGER :: handle, i, iounit, ispin, iterstep, iv, &
68 : iw, ndigits, nhist, nmo(2), nspins
69 : LOGICAL :: append, ionode, should_output
70 : REAL(KIND=dp) :: de, e1, e2, e_fermi(2), emax, emin, eval
71 40 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ehist, hist, occval
72 40 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, occupation_numbers
73 : TYPE(cp_logger_type), POINTER :: logger
74 : TYPE(mo_set_type), POINTER :: mo_set
75 :
76 40 : NULLIFY (logger)
77 80 : logger => cp_get_default_logger()
78 : ionode = logger%para_env%is_source()
79 : should_output = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, &
80 40 : "PRINT%DOS"), cp_p_file)
81 40 : iounit = cp_logger_get_default_io_unit(logger)
82 40 : IF ((.NOT. should_output)) RETURN
83 :
84 40 : CALL timeset(routineN, handle)
85 40 : iterstep = logger%iter_info%iteration(logger%iter_info%n_rlevel)
86 :
87 40 : IF (iounit > 0) WRITE (UNIT=iounit, FMT='(/,(T3,A,T61,I10))') &
88 20 : " Calculate DOS at iteration step ", iterstep
89 :
90 40 : CALL section_vals_val_get(dft_section, "PRINT%DOS%DELTA_E", r_val=de)
91 40 : CALL section_vals_val_get(dft_section, "PRINT%PDOS%APPEND", l_val=append)
92 40 : CALL section_vals_val_get(dft_section, "PRINT%DOS%NDIGITS", i_val=ndigits)
93 40 : IF (append .AND. iterstep > 1) THEN
94 0 : my_pos = "APPEND"
95 : ELSE
96 40 : my_pos = "REWIND"
97 : END IF
98 40 : ndigits = MIN(MAX(ndigits, 1), 10)
99 :
100 40 : emin = 1.e10_dp
101 40 : emax = -1.e10_dp
102 40 : nspins = SIZE(mos)
103 40 : nmo(:) = 0
104 :
105 80 : DO ispin = 1, nspins
106 40 : mo_set => mos(ispin)
107 40 : CALL get_mo_set(mo_set=mo_set, nmo=nmo(ispin), mu=e_fermi(ispin))
108 40 : eigenvalues => mo_set%eigenvalues
109 468 : e1 = MINVAL(eigenvalues(1:nmo(ispin)))
110 468 : e2 = MAXVAL(eigenvalues(1:nmo(ispin)))
111 40 : emin = MIN(emin, e1)
112 80 : emax = MAX(emax, e2)
113 : END DO
114 :
115 40 : IF (de > 0.0_dp) THEN
116 34 : nhist = NINT((emax - emin)/de) + 1
117 272 : ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
118 47950 : hist = 0.0_dp
119 47950 : occval = 0.0_dp
120 47950 : ehist = 0.0_dp
121 68 : DO ispin = 1, nspins
122 34 : mo_set => mos(ispin)
123 34 : occupation_numbers => mo_set%occupation_numbers
124 34 : eigenvalues => mo_set%eigenvalues
125 284 : DO i = 1, nmo(ispin)
126 250 : eval = eigenvalues(i) - emin
127 250 : iv = NINT(eval/de) + 1
128 250 : CPASSERT((iv > 0) .AND. (iv <= nhist))
129 250 : hist(iv, ispin) = hist(iv, ispin) + 1.0_dp
130 284 : occval(iv, ispin) = occval(iv, ispin) + occupation_numbers(i)
131 : END DO
132 47950 : hist(:, ispin) = hist(:, ispin)/REAL(nmo(ispin), KIND=dp)
133 : END DO
134 47916 : DO i = 1, nhist
135 95798 : ehist(i, 1:nspins) = emin + (i - 1)*de
136 : END DO
137 : ELSE
138 18 : nhist = MAXVAL(nmo)
139 48 : ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
140 150 : hist = 0.0_dp
141 150 : occval = 0.0_dp
142 150 : ehist = 0.0_dp
143 12 : DO ispin = 1, nspins
144 6 : mo_set => mos(ispin)
145 6 : occupation_numbers => mo_set%occupation_numbers
146 6 : eigenvalues => mo_set%eigenvalues
147 144 : DO i = 1, nmo(ispin)
148 138 : ehist(i, ispin) = eigenvalues(i)
149 138 : hist(i, ispin) = 1.0_dp
150 144 : occval(i, ispin) = occupation_numbers(i)
151 : END DO
152 150 : hist(:, ispin) = hist(:, ispin)/REAL(nmo(ispin), KIND=dp)
153 : END DO
154 : END IF
155 :
156 40 : my_act = "WRITE"
157 : iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%DOS", &
158 : extension=".dos", file_position=my_pos, file_action=my_act, &
159 40 : file_form="FORMATTED")
160 40 : IF (iw > 0) THEN
161 20 : IF (nspins == 2) THEN
162 : WRITE (UNIT=iw, FMT="(T2,A,I0,A,2F12.6)") &
163 0 : "# DOS at iteration step i = ", iterstep, ", E_Fermi[a.u.] = ", e_fermi(1:2)
164 0 : WRITE (UNIT=iw, FMT="(T2,A, A)") "# Energy[a.u.] Alpha_Density Occupation", &
165 0 : " Energy[a.u.] Beta_Density Occupation"
166 : ! (2(F15.8,2F15.ndigits))
167 0 : WRITE (UNIT=fmtstr_data, FMT="(A,I0,A)") "(2(F15.8,2F15.", ndigits, "))"
168 : ELSE
169 : WRITE (UNIT=iw, FMT="(T2,A,I0,A,F12.6)") &
170 20 : "# DOS at iteration step i = ", iterstep, ", E_Fermi[a.u.] = ", e_fermi(1)
171 20 : WRITE (UNIT=iw, FMT="(T2,A)") "# Energy[a.u.] Density Occupation"
172 : ! (F15.8,2F15.ndigits)
173 20 : WRITE (UNIT=fmtstr_data, FMT="(A,I0,A)") "(F15.8,2F15.", ndigits, ")"
174 : END IF
175 24030 : DO i = 1, nhist
176 24030 : IF (nspins == 2) THEN
177 0 : e1 = ehist(i, 1)
178 0 : e2 = ehist(i, 2)
179 : ! fmtstr_data == "(2(F15.8,2F15.xx))"
180 0 : WRITE (UNIT=iw, FMT=fmtstr_data) e1, hist(i, 1), occval(i, 1), &
181 0 : e2, hist(i, 2), occval(i, 2)
182 : ELSE
183 24010 : eval = ehist(i, 1)
184 : ! fmtstr_data == "(F15.8,2F15.xx)"
185 24010 : WRITE (UNIT=iw, FMT=fmtstr_data) eval, hist(i, 1), occval(i, 1)
186 : END IF
187 : END DO
188 : END IF
189 40 : CALL cp_print_key_finished_output(iw, logger, dft_section, "PRINT%DOS")
190 40 : DEALLOCATE (hist, occval, ehist)
191 :
192 40 : CALL timestop(handle)
193 :
194 40 : END SUBROUTINE calculate_dos
195 :
196 : ! **************************************************************************************************
197 : !> \brief Compute and write density of states (kpoints)
198 : !> \param qs_env ...
199 : !> \param dft_section ...
200 : !> \date 26.02.2008
201 : !> \par History:
202 : !> \author JGH
203 : !> \version 1.0
204 : ! **************************************************************************************************
205 4 : SUBROUTINE calculate_dos_kp(qs_env, dft_section)
206 :
207 : TYPE(qs_environment_type), POINTER :: qs_env
208 : TYPE(section_vals_type), POINTER :: dft_section
209 :
210 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_dos_kp'
211 :
212 : CHARACTER(LEN=16) :: fmtstr_data
213 : CHARACTER(LEN=default_string_length) :: my_act, my_pos
214 : INTEGER :: handle, i, ik, iounit, ispin, iterstep, &
215 : iv, iw, ndigits, nhist, nmo(2), &
216 : nmo_kp, nspins
217 4 : INTEGER, DIMENSION(:), POINTER :: nkp_grid
218 : LOGICAL :: append, explicit, ionode, should_output
219 : REAL(KIND=dp) :: de, e1, e2, emax, emin, eval, wkp
220 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ehist, hist, occval
221 4 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, occupation_numbers
222 : TYPE(cp_logger_type), POINTER :: logger
223 : TYPE(dft_control_type), POINTER :: dft_control
224 : TYPE(kpoint_type), POINTER :: kpoints
225 4 : TYPE(mo_set_type), DIMENSION(:, :), POINTER :: mos
226 : TYPE(mo_set_type), POINTER :: mo_set
227 : TYPE(mp_para_env_type), POINTER :: para_env
228 :
229 4 : NULLIFY (logger, kpoints)
230 8 : logger => cp_get_default_logger()
231 : ionode = logger%para_env%is_source()
232 : should_output = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, &
233 4 : "PRINT%DOS"), cp_p_file)
234 4 : iounit = cp_logger_get_default_io_unit(logger)
235 4 : IF ((.NOT. should_output)) RETURN
236 :
237 4 : CALL timeset(routineN, handle)
238 4 : iterstep = logger%iter_info%iteration(logger%iter_info%n_rlevel)
239 :
240 : ! check whether the user requested a different MP grid for the DOS
241 4 : CALL section_vals_val_get(dft_section, "PRINT%DOS%MP_GRID", i_vals=nkp_grid, explicit=explicit)
242 :
243 4 : IF (explicit) THEN
244 : ! make sure is a valid grid
245 0 : DO i = 1, 3
246 0 : IF (nkp_grid(i) < 1) THEN
247 : WRITE (UNIT=iounit, FMT='(T4,A,I3,A,I1)') &
248 0 : "Invalid kpoint grid for DOS ", nkp_grid(i), " in dimension ", i
249 0 : CPABORT("")
250 : END IF
251 : END DO
252 : ! calculate orbitals and energies
253 0 : CALL calculate_kp_orbitals(qs_env, kpoints, "MONKHORST-PACK", 0, nkp_grid)
254 : ELSE
255 : ! use the kpoints from the environment
256 4 : CALL get_qs_env(qs_env, kpoints=kpoints)
257 : END IF
258 :
259 4 : IF (iounit > 0) WRITE (UNIT=iounit, FMT='(/,(T3,A,T61,I10))') &
260 2 : " Calculate DOS at iteration step ", iterstep
261 : IF (iounit > 0) WRITE (UNIT=iounit, FMT='((T3,A,3I3,A))') &
262 2 : " Using a", kpoints%nkp_grid(:), ' '//TRIM(kpoints%kp_scheme)//' grid'
263 :
264 4 : CALL section_vals_val_get(dft_section, "PRINT%DOS%DELTA_E", r_val=de)
265 4 : CALL section_vals_val_get(dft_section, "PRINT%DOS%APPEND", l_val=append)
266 4 : CALL section_vals_val_get(dft_section, "PRINT%DOS%NDIGITS", i_val=ndigits)
267 : ! ensure a lower value for the histogram width
268 4 : de = MAX(de, 0.00001_dp)
269 4 : IF (append .AND. iterstep > 1) THEN
270 0 : my_pos = "APPEND"
271 : ELSE
272 4 : my_pos = "REWIND"
273 : END IF
274 4 : ndigits = MIN(MAX(ndigits, 1), 10)
275 :
276 4 : CALL get_qs_env(qs_env, dft_control=dft_control)
277 4 : nspins = dft_control%nspins
278 4 : para_env => kpoints%para_env_inter_kp
279 :
280 4 : emin = 1.e10_dp
281 4 : emax = -1.e10_dp
282 4 : nmo(:) = 0
283 4 : IF (kpoints%nkp /= 0) THEN
284 16 : DO ik = 1, SIZE(kpoints%kp_env)
285 12 : mos => kpoints%kp_env(ik)%kpoint_env%mos
286 12 : CPASSERT(ASSOCIATED(mos))
287 28 : DO ispin = 1, nspins
288 12 : mo_set => mos(1, ispin)
289 12 : CALL get_mo_set(mo_set=mo_set, nmo=nmo_kp)
290 12 : eigenvalues => mo_set%eigenvalues
291 440 : e1 = MINVAL(eigenvalues(1:nmo_kp))
292 440 : e2 = MAXVAL(eigenvalues(1:nmo_kp))
293 12 : emin = MIN(emin, e1)
294 12 : emax = MAX(emax, e2)
295 36 : nmo(ispin) = MAX(nmo(ispin), nmo_kp)
296 : END DO
297 : END DO
298 : END IF
299 4 : CALL para_env%min(emin)
300 4 : CALL para_env%max(emax)
301 4 : CALL para_env%max(nmo)
302 :
303 4 : nhist = NINT((emax - emin)/de) + 1
304 32 : ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
305 9096 : hist = 0.0_dp
306 9096 : occval = 0.0_dp
307 9096 : ehist = 0.0_dp
308 :
309 4 : IF (kpoints%nkp /= 0) THEN
310 16 : DO ik = 1, SIZE(kpoints%kp_env)
311 12 : mos => kpoints%kp_env(ik)%kpoint_env%mos
312 12 : wkp = kpoints%kp_env(ik)%kpoint_env%wkp
313 28 : DO ispin = 1, nspins
314 12 : mo_set => mos(1, ispin)
315 12 : occupation_numbers => mo_set%occupation_numbers
316 12 : eigenvalues => mo_set%eigenvalues
317 440 : DO i = 1, nmo(ispin)
318 416 : eval = eigenvalues(i) - emin
319 416 : iv = NINT(eval/de) + 1
320 416 : CPASSERT((iv > 0) .AND. (iv <= nhist))
321 416 : hist(iv, ispin) = hist(iv, ispin) + wkp
322 428 : occval(iv, ispin) = occval(iv, ispin) + wkp*occupation_numbers(i)
323 : END DO
324 : END DO
325 : END DO
326 : END IF
327 4 : CALL para_env%sum(hist)
328 4 : CALL para_env%sum(occval)
329 8 : DO ispin = 1, nspins
330 9096 : hist(:, ispin) = hist(:, ispin)/REAL(nmo(ispin), KIND=dp)
331 : END DO
332 9092 : DO i = 1, nhist
333 18180 : ehist(i, 1:nspins) = emin + (i - 1)*de
334 : END DO
335 :
336 4 : my_act = "WRITE"
337 : iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%DOS", &
338 : extension=".dos", file_position=my_pos, file_action=my_act, &
339 4 : file_form="FORMATTED")
340 4 : IF (iw > 0) THEN
341 2 : IF (nspins == 2) THEN
342 0 : WRITE (UNIT=iw, FMT="(T2,A,I0)") "# DOS at iteration step i = ", iterstep
343 0 : WRITE (UNIT=iw, FMT="(T2,A,A)") "# Energy[a.u.] Alpha_Density Occupation", &
344 0 : " Beta_Density Occupation"
345 : ! (F15.8,4F15.ndigits)
346 0 : WRITE (UNIT=fmtstr_data, FMT="(A,I0,A)") "(F15.8,4F15.", ndigits, ")"
347 : ELSE
348 2 : WRITE (UNIT=iw, FMT="(T2,A,I0)") "# DOS at iteration step i = ", iterstep
349 2 : WRITE (UNIT=iw, FMT="(T2,A)") "# Energy[a.u.] Density Occupation"
350 : ! (F15.8,2F15.ndigits)
351 2 : WRITE (UNIT=fmtstr_data, FMT="(A,I0,A)") "(F15.8,2F15.", ndigits, ")"
352 : END IF
353 4546 : DO i = 1, nhist
354 4544 : eval = emin + (i - 1)*de
355 4546 : IF (nspins == 2) THEN
356 : ! fmtstr_data == "(F15.8,4F15.xx)"
357 0 : WRITE (UNIT=iw, FMT=fmtstr_data) eval, hist(i, 1), occval(i, 1), &
358 0 : hist(i, 2), occval(i, 2)
359 : ELSE
360 : ! fmtstr_data == "(F15.8,2F15.xx)"
361 4544 : WRITE (UNIT=iw, FMT=fmtstr_data) eval, hist(i, 1), occval(i, 1)
362 : END IF
363 : END DO
364 : END IF
365 4 : CALL cp_print_key_finished_output(iw, logger, dft_section, "PRINT%DOS")
366 4 : DEALLOCATE (hist, occval)
367 :
368 : ! destroy the extra k-point set if it was created
369 4 : IF (explicit) THEN
370 0 : CALL kpoint_release(kpoints)
371 : END IF
372 :
373 4 : CALL timestop(handle)
374 :
375 16 : END SUBROUTINE calculate_dos_kp
376 :
377 : END MODULE qs_dos
378 :
|
