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