(git:419edc0)
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
19 USE cp_output_handling, ONLY: cp_p_file,&
25 USE kinds, ONLY: default_string_length,&
26 dp
27 USE kpoint_types, ONLY: kpoint_release,&
33 USE qs_mo_types, ONLY: get_mo_set,&
35#include "./base/base_uses.f90"
36
37 IMPLICIT NONE
38
39 PRIVATE
40
41 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dos'
42
44
45! **************************************************************************************************
46
47CONTAINS
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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: ehist, hist, occval
72 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 NULLIFY (logger)
77 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 "PRINT%DOS"), cp_p_file)
81 iounit = cp_logger_get_default_io_unit(logger)
82 IF ((.NOT. should_output)) RETURN
83
84 CALL timeset(routinen, handle)
85 iterstep = logger%iter_info%iteration(logger%iter_info%n_rlevel)
86
87 IF (iounit > 0) WRITE (unit=iounit, fmt='(/,(T3,A,T61,I10))') &
88 " Calculate DOS at iteration step ", iterstep
89
90 CALL section_vals_val_get(dft_section, "PRINT%DOS%DELTA_E", r_val=de)
91 CALL section_vals_val_get(dft_section, "PRINT%PDOS%APPEND", l_val=append)
92 CALL section_vals_val_get(dft_section, "PRINT%DOS%NDIGITS", i_val=ndigits)
93 IF (append .AND. iterstep > 1) THEN
94 my_pos = "APPEND"
95 ELSE
96 my_pos = "REWIND"
97 END IF
98 ndigits = min(max(ndigits, 1), 10)
99
100 emin = 1.e10_dp
101 emax = -1.e10_dp
102 nspins = SIZE(mos)
103 nmo(:) = 0
104
105 DO ispin = 1, nspins
106 mo_set => mos(ispin)
107 CALL get_mo_set(mo_set=mo_set, nmo=nmo(ispin), mu=e_fermi(ispin))
108 eigenvalues => mo_set%eigenvalues
109 e1 = minval(eigenvalues(1:nmo(ispin)))
110 e2 = maxval(eigenvalues(1:nmo(ispin)))
111 emin = min(emin, e1)
112 emax = max(emax, e2)
113 END DO
114
115 IF (de > 0.0_dp) THEN
116 nhist = nint((emax - emin)/de) + 1
117 ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
118 hist = 0.0_dp
119 occval = 0.0_dp
120 ehist = 0.0_dp
121 DO ispin = 1, nspins
122 mo_set => mos(ispin)
123 occupation_numbers => mo_set%occupation_numbers
124 eigenvalues => mo_set%eigenvalues
125 DO i = 1, nmo(ispin)
126 eval = eigenvalues(i) - emin
127 iv = nint(eval/de) + 1
128 cpassert((iv > 0) .AND. (iv <= nhist))
129 hist(iv, ispin) = hist(iv, ispin) + 1.0_dp
130 occval(iv, ispin) = occval(iv, ispin) + occupation_numbers(i)
131 END DO
132 hist(:, ispin) = hist(:, ispin)/real(nmo(ispin), kind=dp)
133 END DO
134 DO i = 1, nhist
135 ehist(i, 1:nspins) = emin + (i - 1)*de
136 END DO
137 ELSE
138 nhist = maxval(nmo)
139 ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
140 hist = 0.0_dp
141 occval = 0.0_dp
142 ehist = 0.0_dp
143 DO ispin = 1, nspins
144 mo_set => mos(ispin)
145 occupation_numbers => mo_set%occupation_numbers
146 eigenvalues => mo_set%eigenvalues
147 DO i = 1, nmo(ispin)
148 ehist(i, ispin) = eigenvalues(i)
149 hist(i, ispin) = 1.0_dp
150 occval(i, ispin) = occupation_numbers(i)
151 END DO
152 hist(:, ispin) = hist(:, ispin)/real(nmo(ispin), kind=dp)
153 END DO
154 END IF
155
156 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 file_form="FORMATTED")
160 IF (iw > 0) THEN
161 IF (nspins == 2) THEN
162 WRITE (unit=iw, fmt="(T2,A,I0,A,2F12.6)") &
163 "# DOS at iteration step i = ", iterstep, ", E_Fermi[a.u.] = ", e_fermi(1:2)
164 WRITE (unit=iw, fmt="(T2,A, A)") "# Energy[a.u.] Alpha_Density Occupation", &
165 " Energy[a.u.] Beta_Density Occupation"
166 ! (2(F15.8,2F15.ndigits))
167 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 "# DOS at iteration step i = ", iterstep, ", E_Fermi[a.u.] = ", e_fermi(1)
171 WRITE (unit=iw, fmt="(T2,A)") "# Energy[a.u.] Density Occupation"
172 ! (F15.8,2F15.ndigits)
173 WRITE (unit=fmtstr_data, fmt="(A,I0,A)") "(F15.8,2F15.", ndigits, ")"
174 END IF
175 DO i = 1, nhist
176 IF (nspins == 2) THEN
177 e1 = ehist(i, 1)
178 e2 = ehist(i, 2)
179 ! fmtstr_data == "(2(F15.8,2F15.xx))"
180 WRITE (unit=iw, fmt=fmtstr_data) e1, hist(i, 1), occval(i, 1), &
181 e2, hist(i, 2), occval(i, 2)
182 ELSE
183 eval = ehist(i, 1)
184 ! fmtstr_data == "(F15.8,2F15.xx)"
185 WRITE (unit=iw, fmt=fmtstr_data) eval, hist(i, 1), occval(i, 1)
186 END IF
187 END DO
188 END IF
189 CALL cp_print_key_finished_output(iw, logger, dft_section, "PRINT%DOS")
190 DEALLOCATE (hist, occval, ehist)
191
192 CALL timestop(handle)
193
194 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 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 INTEGER, DIMENSION(:), POINTER :: nkp_grid
218 LOGICAL :: append, explicit, ionode, should_output
219 REAL(kind=dp) :: de, e1, e2, emax, emin, eval, wkp
220 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: ehist, hist, occval
221 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 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 NULLIFY (logger, kpoints)
230 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 "PRINT%DOS"), cp_p_file)
234 iounit = cp_logger_get_default_io_unit(logger)
235 IF ((.NOT. should_output)) RETURN
236
237 CALL timeset(routinen, handle)
238 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 CALL section_vals_val_get(dft_section, "PRINT%DOS%MP_GRID", i_vals=nkp_grid, explicit=explicit)
242
243 IF (explicit) THEN
244 ! make sure is a valid grid
245 DO i = 1, 3
246 IF (nkp_grid(i) < 1) THEN
247 WRITE (unit=iounit, fmt='(T4,A,I3,A,I1)') &
248 "Invalid kpoint grid for DOS ", nkp_grid(i), " in dimension ", i
249 cpabort("")
250 END IF
251 END DO
252 ! calculate orbitals and energies
253 CALL calculate_kp_orbitals(qs_env, kpoints, "MONKHORST-PACK", 0, nkp_grid)
254 ELSE
255 ! use the kpoints from the environment
256 CALL get_qs_env(qs_env, kpoints=kpoints)
257 END IF
258
259 IF (iounit > 0) WRITE (unit=iounit, fmt='(/,(T3,A,T61,I10))') &
260 " Calculate DOS at iteration step ", iterstep
261 IF (iounit > 0) WRITE (unit=iounit, fmt='((T3,A,3I3,A))') &
262 " Using a", kpoints%nkp_grid(:), ' '//trim(kpoints%kp_scheme)//' grid'
263
264 CALL section_vals_val_get(dft_section, "PRINT%DOS%DELTA_E", r_val=de)
265 CALL section_vals_val_get(dft_section, "PRINT%DOS%APPEND", l_val=append)
266 CALL section_vals_val_get(dft_section, "PRINT%DOS%NDIGITS", i_val=ndigits)
267 ! ensure a lower value for the histogram width
268 de = max(de, 0.00001_dp)
269 IF (append .AND. iterstep > 1) THEN
270 my_pos = "APPEND"
271 ELSE
272 my_pos = "REWIND"
273 END IF
274 ndigits = min(max(ndigits, 1), 10)
275
276 CALL get_qs_env(qs_env, dft_control=dft_control)
277 nspins = dft_control%nspins
278 para_env => kpoints%para_env_inter_kp
279
280 emin = 1.e10_dp
281 emax = -1.e10_dp
282 nmo(:) = 0
283 IF (kpoints%nkp /= 0) THEN
284 DO ik = 1, SIZE(kpoints%kp_env)
285 mos => kpoints%kp_env(ik)%kpoint_env%mos
286 cpassert(ASSOCIATED(mos))
287 DO ispin = 1, nspins
288 mo_set => mos(1, ispin)
289 CALL get_mo_set(mo_set=mo_set, nmo=nmo_kp)
290 eigenvalues => mo_set%eigenvalues
291 e1 = minval(eigenvalues(1:nmo_kp))
292 e2 = maxval(eigenvalues(1:nmo_kp))
293 emin = min(emin, e1)
294 emax = max(emax, e2)
295 nmo(ispin) = max(nmo(ispin), nmo_kp)
296 END DO
297 END DO
298 END IF
299 CALL para_env%min(emin)
300 CALL para_env%max(emax)
301 CALL para_env%max(nmo)
302
303 nhist = nint((emax - emin)/de) + 1
304 ALLOCATE (hist(nhist, nspins), occval(nhist, nspins), ehist(nhist, nspins))
305 hist = 0.0_dp
306 occval = 0.0_dp
307 ehist = 0.0_dp
308
309 IF (kpoints%nkp /= 0) THEN
310 DO ik = 1, SIZE(kpoints%kp_env)
311 mos => kpoints%kp_env(ik)%kpoint_env%mos
312 wkp = kpoints%kp_env(ik)%kpoint_env%wkp
313 DO ispin = 1, nspins
314 mo_set => mos(1, ispin)
315 occupation_numbers => mo_set%occupation_numbers
316 eigenvalues => mo_set%eigenvalues
317 DO i = 1, nmo(ispin)
318 eval = eigenvalues(i) - emin
319 iv = nint(eval/de) + 1
320 cpassert((iv > 0) .AND. (iv <= nhist))
321 hist(iv, ispin) = hist(iv, ispin) + wkp
322 occval(iv, ispin) = occval(iv, ispin) + wkp*occupation_numbers(i)
323 END DO
324 END DO
325 END DO
326 END IF
327 CALL para_env%sum(hist)
328 CALL para_env%sum(occval)
329 DO ispin = 1, nspins
330 hist(:, ispin) = hist(:, ispin)/real(nmo(ispin), kind=dp)
331 END DO
332 DO i = 1, nhist
333 ehist(i, 1:nspins) = emin + (i - 1)*de
334 END DO
335
336 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 file_form="FORMATTED")
340 IF (iw > 0) THEN
341 IF (nspins == 2) THEN
342 WRITE (unit=iw, fmt="(T2,A,I0)") "# DOS at iteration step i = ", iterstep
343 WRITE (unit=iw, fmt="(T2,A,A)") "# Energy[a.u.] Alpha_Density Occupation", &
344 " Beta_Density Occupation"
345 ! (F15.8,4F15.ndigits)
346 WRITE (unit=fmtstr_data, fmt="(A,I0,A)") "(F15.8,4F15.", ndigits, ")"
347 ELSE
348 WRITE (unit=iw, fmt="(T2,A,I0)") "# DOS at iteration step i = ", iterstep
349 WRITE (unit=iw, fmt="(T2,A)") "# Energy[a.u.] Density Occupation"
350 ! (F15.8,2F15.ndigits)
351 WRITE (unit=fmtstr_data, fmt="(A,I0,A)") "(F15.8,2F15.", ndigits, ")"
352 END IF
353 DO i = 1, nhist
354 eval = emin + (i - 1)*de
355 IF (nspins == 2) THEN
356 ! fmtstr_data == "(F15.8,4F15.xx)"
357 WRITE (unit=iw, fmt=fmtstr_data) eval, hist(i, 1), occval(i, 1), &
358 hist(i, 2), occval(i, 2)
359 ELSE
360 ! fmtstr_data == "(F15.8,2F15.xx)"
361 WRITE (unit=iw, fmt=fmtstr_data) eval, hist(i, 1), occval(i, 1)
362 END IF
363 END DO
364 END IF
365 CALL cp_print_key_finished_output(iw, logger, dft_section, "PRINT%DOS")
366 DEALLOCATE (hist, occval)
367
368 ! destroy the extra k-point set if it was created
369 IF (explicit) THEN
370 CALL kpoint_release(kpoints)
371 END IF
372
373 CALL timestop(handle)
374
375 END SUBROUTINE calculate_dos_kp
376
377END MODULE qs_dos
378
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.
subroutine, public kpoint_release(kpoint)
Release a kpoint environment, deallocate all data.
Interface to the message passing library MPI.
Calculation of band structures.
subroutine, public calculate_kp_orbitals(qs_env, kpoint, scheme, nadd, mp_grid, kpgeneral, group_size_ext)
diagonalize KS matrices at a set of kpoints
Calculation and writing of density of states.
Definition qs_dos.F:14
subroutine, public calculate_dos(mos, dft_section)
Compute and write density of states.
Definition qs_dos.F:59
subroutine, public calculate_dos_kp(qs_env, dft_section)
Compute and write density of states (kpoints)
Definition qs_dos.F:206
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 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.
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Contains information about kpoints.
stores all the informations relevant to an mpi environment