(git:ed6f26b)
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-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! **************************************************************************************************
14MODULE 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
45CONTAINS
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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: ehist, hist, occval
70 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 NULLIFY (logger)
75 logger => cp_get_default_logger()
76 ionode = logger%para_env%is_source()
77 should_output = btest(cp_print_key_should_output(logger%iter_info, dft_section, &
78 "PRINT%DOS"), cp_p_file)
79 iounit = cp_logger_get_default_io_unit(logger)
80 IF ((.NOT. should_output)) RETURN
81
82 CALL timeset(routinen, handle)
83 iterstep = logger%iter_info%iteration(logger%iter_info%n_rlevel)
84
85 IF (iounit > 0) WRITE (unit=iounit, fmt='(/,(T3,A,T61,I10))') &
86 " Calculate DOS at iteration step ", iterstep
87
88 CALL section_vals_val_get(dft_section, "PRINT%DOS%DELTA_E", r_val=de)
89 CALL section_vals_val_get(dft_section, "PRINT%PDOS%APPEND", l_val=append)
90 CALL section_vals_val_get(dft_section, "PRINT%DOS%NDIGITS", i_val=ndigits)
91 IF (append .AND. iterstep > 1) THEN
92 my_pos = "APPEND"
93 ELSE
94 my_pos = "REWIND"
95 END IF
96 ndigits = min(max(ndigits, 1), 10)
97
98 emin = 1.e10_dp
99 emax = -1.e10_dp
100 nspins = SIZE(mos)
101 nmo(:) = 0
102
103 DO ispin = 1, nspins
104 mo_set => mos(ispin)
105 CALL get_mo_set(mo_set=mo_set, nmo=nmo(ispin), mu=e_fermi(ispin))
106 eigenvalues => mo_set%eigenvalues
107 e1 = minval(eigenvalues(1:nmo(ispin)))
108 e2 = maxval(eigenvalues(1:nmo(ispin)))
109 emin = min(emin, e1)
110 emax = max(emax, e2)
111 END DO
112
113 IF (de > 0.0_dp) THEN
114 nhist = nint((emax - emin)/de) + 1
115 ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
116 hist = 0.0_dp
117 occval = 0.0_dp
118 ehist = 0.0_dp
119 DO ispin = 1, nspins
120 mo_set => mos(ispin)
121 occupation_numbers => mo_set%occupation_numbers
122 eigenvalues => mo_set%eigenvalues
123 DO i = 1, nmo(ispin)
124 eval = eigenvalues(i) - emin
125 iv = nint(eval/de) + 1
126 cpassert((iv > 0) .AND. (iv <= nhist))
127 hist(iv, ispin) = hist(iv, ispin) + 1.0_dp
128 occval(iv, ispin) = occval(iv, ispin) + occupation_numbers(i)
129 END DO
130 hist(:, ispin) = hist(:, ispin)/real(nmo(ispin), kind=dp)
131 END DO
132 DO i = 1, nhist
133 ehist(i, 1:nspins) = emin + (i - 1)*de
134 END DO
135 ELSE
136 nhist = maxval(nmo)
137 ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
138 hist = 0.0_dp
139 occval = 0.0_dp
140 ehist = 0.0_dp
141 DO ispin = 1, nspins
142 mo_set => mos(ispin)
143 occupation_numbers => mo_set%occupation_numbers
144 eigenvalues => mo_set%eigenvalues
145 DO i = 1, nmo(ispin)
146 ehist(i, ispin) = eigenvalues(i)
147 hist(i, ispin) = 1.0_dp
148 occval(i, ispin) = occupation_numbers(i)
149 END DO
150 hist(:, ispin) = hist(:, ispin)/real(nmo(ispin), kind=dp)
151 END DO
152 END IF
153
154 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 file_form="FORMATTED")
158 IF (iw > 0) THEN
159 IF (nspins == 2) THEN
160 WRITE (unit=iw, fmt="(T2,A,I0,A,2F12.6)") &
161 "# DOS at iteration step i = ", iterstep, ", E_Fermi[a.u.] = ", e_fermi(1:2)
162 WRITE (unit=iw, fmt="(T2,A, A)") "# Energy[a.u.] Alpha_Density Occupation", &
163 " Energy[a.u.] Beta_Density Occupation"
164 ! (2(F15.8,2F15.ndigits))
165 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 "# DOS at iteration step i = ", iterstep, ", E_Fermi[a.u.] = ", e_fermi(1)
169 WRITE (unit=iw, fmt="(T2,A)") "# Energy[a.u.] Density Occupation"
170 ! (F15.8,2F15.ndigits)
171 WRITE (unit=fmtstr_data, fmt="(A,I0,A)") "(F15.8,2F15.", ndigits, ")"
172 END IF
173 DO i = 1, nhist
174 IF (nspins == 2) THEN
175 e1 = ehist(i, 1)
176 e2 = ehist(i, 2)
177 ! fmtstr_data == "(2(F15.8,2F15.xx))"
178 WRITE (unit=iw, fmt=fmtstr_data) e1, hist(i, 1), occval(i, 1), &
179 e2, hist(i, 2), occval(i, 2)
180 ELSE
181 eval = ehist(i, 1)
182 ! fmtstr_data == "(F15.8,2F15.xx)"
183 WRITE (unit=iw, fmt=fmtstr_data) eval, hist(i, 1), occval(i, 1)
184 END IF
185 END DO
186 END IF
187 CALL cp_print_key_finished_output(iw, logger, dft_section, "PRINT%DOS")
188 DEALLOCATE (hist, occval, ehist)
189
190 CALL timestop(handle)
191
192 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: ehist, hist, occval
220 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues, occupation_numbers
221 TYPE(cp_logger_type), POINTER :: logger
222 TYPE(dft_control_type), POINTER :: dft_control
223 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 NULLIFY (logger)
228 logger => cp_get_default_logger()
229 ionode = logger%para_env%is_source()
230 should_output = btest(cp_print_key_should_output(logger%iter_info, dft_section, &
231 "PRINT%DOS"), cp_p_file)
232 iounit = cp_logger_get_default_io_unit(logger)
233 IF ((.NOT. should_output)) RETURN
234
235 CALL timeset(routinen, handle)
236 iterstep = logger%iter_info%iteration(logger%iter_info%n_rlevel)
237
238 IF (iounit > 0) WRITE (unit=iounit, fmt='(/,(T3,A,T61,I10))') &
239 " Calculate DOS at iteration step ", iterstep
240
241 CALL section_vals_val_get(dft_section, "PRINT%DOS%DELTA_E", r_val=de)
242 CALL section_vals_val_get(dft_section, "PRINT%DOS%APPEND", l_val=append)
243 CALL section_vals_val_get(dft_section, "PRINT%DOS%NDIGITS", i_val=ndigits)
244 ! ensure a lower value for the histogram width
245 de = max(de, 0.00001_dp)
246 IF (append .AND. iterstep > 1) THEN
247 my_pos = "APPEND"
248 ELSE
249 my_pos = "REWIND"
250 END IF
251 ndigits = min(max(ndigits, 1), 10)
252
253 CALL get_qs_env(qs_env, dft_control=dft_control)
254 nspins = dft_control%nspins
255 para_env => kpoints%para_env_inter_kp
256
257 emin = 1.e10_dp
258 emax = -1.e10_dp
259 nmo(:) = 0
260 IF (kpoints%nkp /= 0) THEN
261 DO ik = 1, SIZE(kpoints%kp_env)
262 mos => kpoints%kp_env(ik)%kpoint_env%mos
263 cpassert(ASSOCIATED(mos))
264 DO ispin = 1, nspins
265 mo_set => mos(1, ispin)
266 CALL get_mo_set(mo_set=mo_set, nmo=nmo_kp)
267 eigenvalues => mo_set%eigenvalues
268 e1 = minval(eigenvalues(1:nmo_kp))
269 e2 = maxval(eigenvalues(1:nmo_kp))
270 emin = min(emin, e1)
271 emax = max(emax, e2)
272 nmo(ispin) = max(nmo(ispin), nmo_kp)
273 END DO
274 END DO
275 END IF
276 CALL para_env%min(emin)
277 CALL para_env%max(emax)
278 CALL para_env%max(nmo)
279
280 nhist = nint((emax - emin)/de) + 1
281 ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
282 hist = 0.0_dp
283 occval = 0.0_dp
284 ehist = 0.0_dp
285
286 IF (kpoints%nkp /= 0) THEN
287 DO ik = 1, SIZE(kpoints%kp_env)
288 mos => kpoints%kp_env(ik)%kpoint_env%mos
289 wkp = kpoints%kp_env(ik)%kpoint_env%wkp
290 DO ispin = 1, nspins
291 mo_set => mos(1, ispin)
292 occupation_numbers => mo_set%occupation_numbers
293 eigenvalues => mo_set%eigenvalues
294 DO i = 1, nmo(ispin)
295 eval = eigenvalues(i) - emin
296 iv = nint(eval/de) + 1
297 cpassert((iv > 0) .AND. (iv <= nhist))
298 hist(iv, ispin) = hist(iv, ispin) + wkp
299 occval(iv, ispin) = occval(iv, ispin) + wkp*occupation_numbers(i)
300 END DO
301 END DO
302 END DO
303 END IF
304 CALL para_env%sum(hist)
305 CALL para_env%sum(occval)
306 DO ispin = 1, nspins
307 hist(:, ispin) = hist(:, ispin)/real(nmo(ispin), kind=dp)
308 END DO
309 DO i = 1, nhist
310 ehist(i, 1:nspins) = emin + (i - 1)*de
311 END DO
312
313 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 file_form="FORMATTED")
317 IF (iw > 0) THEN
318 IF (nspins == 2) THEN
319 WRITE (unit=iw, fmt="(T2,A,I0)") "# DOS at iteration step i = ", iterstep
320 WRITE (unit=iw, fmt="(T2,A,A)") "# Energy[a.u.] Alpha_Density Occupation", &
321 " Beta_Density Occupation"
322 ! (F15.8,4F15.ndigits)
323 WRITE (unit=fmtstr_data, fmt="(A,I0,A)") "(F15.8,4F15.", ndigits, ")"
324 ELSE
325 WRITE (unit=iw, fmt="(T2,A,I0)") "# DOS at iteration step i = ", iterstep
326 WRITE (unit=iw, fmt="(T2,A)") "# Energy[a.u.] Density Occupation"
327 ! (F15.8,2F15.ndigits)
328 WRITE (unit=fmtstr_data, fmt="(A,I0,A)") "(F15.8,2F15.", ndigits, ")"
329 END IF
330 DO i = 1, nhist
331 eval = emin + (i - 1)*de
332 IF (nspins == 2) THEN
333 ! fmtstr_data == "(F15.8,4F15.xx)"
334 WRITE (unit=iw, fmt=fmtstr_data) eval, hist(i, 1), occval(i, 1), &
335 hist(i, 2), occval(i, 2)
336 ELSE
337 ! fmtstr_data == "(F15.8,2F15.xx)"
338 WRITE (unit=iw, fmt=fmtstr_data) eval, hist(i, 1), occval(i, 1)
339 END IF
340 END DO
341 END IF
342 CALL cp_print_key_finished_output(iw, logger, dft_section, "PRINT%DOS")
343 DEALLOCATE (hist, occval)
344
345 CALL timestop(handle)
346
347 END SUBROUTINE calculate_dos_kp
348
349END MODULE qs_dos
350
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:515
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:1063
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:1047
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
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
qs_dos
Calculation and writing of density of states.
Definition qs_dos.F:14
qs_dos::calculate_dos_kp
subroutine, public calculate_dos_kp(kpoints, qs_env, dft_section)
Compute and write density of states (kpoints)
Definition qs_dos.F:205
qs_dos::calculate_dos
subroutine, public calculate_dos(mos, dft_section)
Compute and write density of states.
Definition qs_dos.F:57
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, 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, 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, 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, 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, ecoul_1c, rho0_s_rs, rho0_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)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
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:397
cp_control_types::dft_control_type
Definition cp_control_types.F:570
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:150
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:217
qs_mo_types::mo_set_type
Definition qs_mo_types.F:46