(git:58e3e09)
qs_dos.F
Go to the documentation of this file.
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
18  cp_logger_type
19  USE cp_output_handling, ONLY: cp_p_file,&
23  USE input_section_types, ONLY: section_vals_type,&
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 
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  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 
349 END MODULE qs_dos
350 
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
various routines to log and control the output. The idea is that decisions about where to log should ...
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...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
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)
...
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,...
integer, parameter, public cp_p_file
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...
objects that represent the structure of input sections and the data contained in an input section
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
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public default_string_length
Definition: kinds.F:57
Types and basic routines needed for a kpoint calculation.
Definition: kpoint_types.F:15
Interface to the message passing library MPI.
Calculation and writing of density of states.
Definition: qs_dos.F:14
subroutine, public calculate_dos_kp(kpoints, qs_env, dft_section)
Compute and write density of states (kpoints)
Definition: qs_dos.F:205
subroutine, public calculate_dos(mos, dft_section)
Compute and write density of states.
Definition: qs_dos.F:57
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_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, 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, rhs)
Get the QUICKSTEP environment.
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
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