From 2c97d88377405e81dc1231e900b5c0047accf21e Mon Sep 17 00:00:00 2001 From: SY Wang Date: Mon, 13 Jul 2026 22:22:13 +0800 Subject: [PATCH] (P)DOS: Restore default behavior of output data (#5564) --- docs/methods/electronic_structure/dos.md | 117 ++++++------- src/input_cp2k_print_dft.F | 159 +++++++++--------- src/qs_dos.F | 175 ++++++++++++++------ src/qs_dos_utils.F | 16 +- src/qs_pdos.F | 199 ++++++++++++++--------- src/qs_scf_post_gpw.F | 21 ++- src/qs_scf_post_tb.F | 15 +- 7 files changed, 412 insertions(+), 290 deletions(-) diff --git a/docs/methods/electronic_structure/dos.md b/docs/methods/electronic_structure/dos.md index bc2adee2e8..08283f4985 100644 --- a/docs/methods/electronic_structure/dos.md +++ b/docs/methods/electronic_structure/dos.md @@ -28,10 +28,10 @@ used to assign the chemical character of the peaks or band edges. ```{important} Most of the DOS/PDOS functionality described below, including the unified -`&DFT%PRINT%DOS` interface and the broadened, directly plottable `.dos` and `.pdos` files, is only -available in CP2K 2026.2 and later versions. In earlier CP2K versions, `.dos` and `.pdos` referred -to state-resolved raw outputs rather than broadened DOS/PDOS curves. The raw output of PDOS is -preserved and now written as `.pdos_raw` through `&DFT%PRINT%DOS%PDOS%RAW` keyword. +`&DFT%PRINT%DOS` interface and the `&CURVE` subsection for broadened, directly plottable +DOS/PDOS curves, is only available in CP2K 2026.2 and later versions. The ordinary `.dos` and +`.pdos` outputs keep the traditional non-broadened data, while `&CURVE` requests additional +broadened curve output. ``` ## Basic CP2K input @@ -46,10 +46,14 @@ additionally write projected DOS files. &DFT &PRINT &DOS - ENERGY_UNIT EV - ENERGY_ZERO AUTO - BROADEN_WIDTH [eV] 0.1 DELTA_E [eV] 0.01 + &CURVE + ENERGY_UNIT EV + ENERGY_ZERO AUTO + &BROADEN + WIDTH [eV] 0.1 + &END BROADEN + &END CURVE &PDOS COMPONENTS F &END PDOS @@ -61,15 +65,18 @@ additionally write projected DOS files. Useful input keywords are documented in the input reference: -- [ENERGY_UNIT](#CP2K_INPUT.FORCE_EVAL.DFT.PRINT.DOS.ENERGY_UNIT): print energies in Hartree or eV. -- [ENERGY_ZERO](#CP2K_INPUT.FORCE_EVAL.DFT.PRINT.DOS.ENERGY_ZERO): choose the reference energy. -- [BROADEN_WIDTH](#CP2K_INPUT.FORCE_EVAL.DFT.PRINT.DOS.BROADEN_WIDTH), - [BROADEN_TYPE](#CP2K_INPUT.FORCE_EVAL.DFT.PRINT.DOS.BROADEN_TYPE), and - [VOIGT_MIXING](#CP2K_INPUT.FORCE_EVAL.DFT.PRINT.DOS.VOIGT_MIXING): choose the line shape. - [DELTA_E](#CP2K_INPUT.FORCE_EVAL.DFT.PRINT.DOS.DELTA_E): choose the energy-grid spacing. - [NLUMO](#CP2K_INPUT.FORCE_EVAL.DFT.PRINT.DOS.NLUMO): requested number of unoccupied states for DOS/PDOS. -- [&PDOS](#CP2K_INPUT.FORCE_EVAL.DFT.PRINT.DOS.LDOS): request projected DOS on each element. +- [&CURVE](#CP2K_INPUT.FORCE_EVAL.DFT.PRINT.DOS.CURVE): request additional broadened, directly + plottable DOS/PDOS curve output. +- [ENERGY_UNIT](#CP2K_INPUT.FORCE_EVAL.DFT.PRINT.DOS.CURVE.ENERGY_UNIT): print curve energies in + Hartree or eV. +- [ENERGY_ZERO](#CP2K_INPUT.FORCE_EVAL.DFT.PRINT.DOS.CURVE.ENERGY_ZERO): choose the reference energy + for curve output. +- [&BROADEN](#CP2K_INPUT.FORCE_EVAL.DFT.PRINT.DOS.CURVE.BROADEN): choose the line shape and + broadening width for curve output. +- [&PDOS](#CP2K_INPUT.FORCE_EVAL.DFT.PRINT.DOS.PDOS): request projected DOS on each element. - [&LDOS](#CP2K_INPUT.FORCE_EVAL.DFT.PRINT.DOS.LDOS) and [&R_LDOS](#CP2K_INPUT.FORCE_EVAL.DFT.PRINT.DOS.R_LDOS): request atom-list or real-space local projected DOS. @@ -78,7 +85,6 @@ For PDOS, there are some additional controls: - [COMPONENTS](#CP2K_INPUT.FORCE_EVAL.DFT.PRINT.DOS.PDOS.COMPONENTS): split angular-momentum channels into individual components. This keyword also applies to `&LDOS` and `&R_LDOS`. -- [RAW](#CP2K_INPUT.FORCE_EVAL.DFT.PRINT.DOS.PDOS.RAW): print raw state-resolved projected weights. ```{note} For diagonalization-based SCF calculations, the unoccupied states used for DOS/PDOS are those @@ -91,8 +97,9 @@ separate post-SCF step. ### Choosing the energy reference -Use [ENERGY_ZERO](#CP2K_INPUT.FORCE_EVAL.DFT.PRINT.DOS.ENERGY_ZERO) to make different calculations -and plots comparable. +Use [ENERGY_ZERO](#CP2K_INPUT.FORCE_EVAL.DFT.PRINT.DOS.CURVE.ENERGY_ZERO) to make different curve +plots comparable. This setting applies to broadened curve output, while ordinary non-curve outputs +are not shifted by this setting. - `FERMI` is normally the appropriate zero for metals and smeared calculations. - `HOCO` places the highest occupied crystal orbital at zero. It is often more convenient for @@ -101,7 +108,7 @@ and plots comparable. - `ABSOLUTE` prints the absolute Kohn-Sham eigenvalues. - `AUTO` selects `FERMI` when smearing or fractional occupations are present, and `HOCO` otherwise. -The output header reports both reference energies and the selected zero, for example: +The curve output header reports both reference energies and the selected zero, for example: ```text # E(Fermi) = 0.223578 a.u. = 6.08391 eV @@ -110,9 +117,10 @@ The output header reports both reference energies and the selected zero, for exa ``` When comparing several structures or charge states, choose the same energy-zero convention for all -plots. If an absolute alignment is needed across different cells or surfaces, additional -electrostatic-potential alignment may be required; changing the DOS zero alone is not a substitute -for such an alignment. +curve plots. In spin-polarized calculations, the two spin channels use the same reference energy so +that alpha and beta curves remain on a common energy axis. If an absolute alignment is needed across +different cells or surfaces, additional electrostatic-potential alignment may be required; changing +the DOS zero alone is not a substitute for such an alignment. ## Broadening, k-points, and gaps @@ -121,26 +129,29 @@ widths, can be affected by the calculation setup. Appropriate choices of the fun k-point sampling, smearing, unoccupied states, and broadening parameters are important for obtaining a reliable interpretation of the electronic structure. -The broadened DOS replaces each discrete eigenvalue by a normalized line shape. -[BROADEN_WIDTH](#CP2K_INPUT.FORCE_EVAL.DFT.PRINT.DOS.BROADEN_WIDTH) is the full width at half -maximum (FWHM). The default value of 0.1 eV is a reasonable choice for most applications. Smaller -values can be useful for resolving narrow features, small gaps, spin splittings, or defect states, -whereas larger values mainly produce smoother curves and may hide such details. +The broadened curve output replaces each discrete eigenvalue by a normalized line shape. The +broadening [WIDTH](#CP2K_INPUT.FORCE_EVAL.DFT.PRINT.DOS.CURVE.BROADEN.WIDTH) is the full width at +half maximum (FWHM). The default value of 0.1 eV is a reasonable choice for most applications. +Smaller values can be useful for resolving narrow features, small gaps, spin splittings, or defect +states, whereas larger values mainly produce smoother curves and may hide such details. Curve output +requires `DELTA_E > 0`; if `DELTA_E <= 0`, CP2K skips the broadened curve output and prints a +warning. DOS and PDOS usually require a denser k-point mesh than calculations aimed only at total energies, forces, or geometry optimization. A too sparse k-point mesh gives a DOS dominated by the discrete sampling of the Brillouin zone, especially near the band edges or the Fermi level. Broadening helps produce a smooth curve, but it should not be used as a substitute for insufficient k-point sampling. -For molecular or very large supercell calculations with only the Gamma point, the DOS is a broadened -representation of discrete levels. It can still be useful for assigning orbital character, but peak -widths are chosen mainly for visualization rather than for Brillouin-zone integration. +For molecular or very large supercell calculations with only the Gamma point, the broadened curve is +a representation of discrete levels. It can still be useful for assigning orbital character, but +peak widths are chosen mainly for visualization rather than for Brillouin-zone integration. ## Interpreting PDOS -With `&PDOS`, CP2K writes projected DOS, usually one file per atomic kind and spin channel. A -`.pdos` file contains a projected `total` column followed by angular-momentum or component-resolved -columns: +With `&PDOS`, CP2K writes projected DOS, usually one file per atomic kind and spin channel. The +ordinary `.pdos` files contain state-resolved projected weights. If `&CURVE` is also requested, CP2K +additionally writes broadened PDOS curve files with a projected `total` column followed by +angular-momentum or component-resolved columns: ```text # Energy[eV] total s p d @@ -153,7 +164,7 @@ into individual components: # Energy[eV] total s py pz px ... ``` -The `total` column in one `.pdos` file is the projected total for that atomic kind, not the total +The `total` column in one PDOS curve file is the projected total for that atomic kind, not the total DOS of the full system. Summing the projected totals over all kinds should reproduce the total DOS up to numerical and output precision, provided that the same states and broadening are used. @@ -178,18 +189,8 @@ contributions: For quantitative orbital populations, integrated PDOS values should be used with care because they depend on the projection scheme, basis set, and chosen energy window. -[`RAW T`](#CP2K_INPUT.FORCE_EVAL.DFT.PRINT.DOS.PDOS.RAW) prints state-resolved projected weights in -addition to the broadened `.pdos` files: - -```text -&DOS - &PDOS - RAW T - &END PDOS -&END DOS -``` - -The `RAW` option is currently not available with k-points calculations. +State-resolved k-point PDOS output is currently not available. K-point PDOS is accumulated with +k-point weights for the broadened curve output. ### Projection used for PDOS @@ -213,38 +214,38 @@ then accumulated with the k-point weights to form the broadened PDOS. ### Using gnuplot -The DOS and PDOS files can be plotted directly with [gnuplot](http://gnuplot.info/), since comment -lines starting with `#` are ignored. For example, if the first column is the energy and the second -column is the total DOS, a simple plot is +The broadened DOS and PDOS curve files can be plotted directly with [gnuplot](http://gnuplot.info/), +since comment lines starting with `#` are ignored. For example, if the first column is the energy +and the second column is the total DOS, a simple plot is ```gnuplot -plot "project-1.dos" using 1:2 with lines title "DOS" +plot "project-curve-1.dos" using 1:2 with lines title "DOS" ``` The plotted energy range can be restricted directly in the `plot` command. For example, to show only states within 5 eV of the chosen energy zero, ```gnuplot -plot [-5:5] "project-1.dos" using 1:2 with lines title "DOS" +plot [-5:5] "project-curve-1.dos" using 1:2 with lines title "DOS" ``` Both the energy and DOS ranges can be limited in the same way: ```gnuplot -plot [-5:5][0:*] "project-1.dos" using 1:2 with lines title "DOS" +plot [-5:5][0:*] "project-curve-1.dos" using 1:2 with lines title "DOS" ``` -For PDOS files, different columns can be plotted together. For example, if the second column is the -kind-resolved total PDOS and the following columns are angular-momentum components, +For PDOS curve files, different columns can be plotted together. For example, if the second column +is the kind-resolved total PDOS and the following columns are angular-momentum components, ```gnuplot -plot [-5:5] "project-k1-1.pdos" using 1:2 with lines title "total", \ +plot [-5:5] "project-k1-curve-1.pdos" using 1:2 with lines title "total", \ "" using 1:3 with lines title "s", \ "" using 1:4 with lines title "p" ``` -When comparing DOS/PDOS with band structures, use the same energy zero, for example `FERMI` for -metals or `HOCO` for many insulating and molecular systems. +When comparing DOS/PDOS curves with band structures, use the same energy zero, for example `FERMI` +for metals or `HOCO` for many insulating and molecular systems. ## Local projected DOS @@ -283,8 +284,8 @@ generic PDOS format should not be interpreted as ordinary ground-state Kohn-Sham A typical XAS_TDP use is to assign intense spectral features to atomic or orbital character, for example to distinguish metal-centered, ligand-centered, or mixed transitions. In this context, -broadening and energy-grid settings remain useful plotting controls, while ordinary DFT-DOS concepts -such as Fermi-level alignment, HOCO alignment, or total DOS are usually not the relevant analysis -target. See [&XAS%PRINT%PDOS](#CP2K_INPUT.FORCE_EVAL.DFT.XAS.PRINT.PDOS) and +`DELTA_E` and `&CURVE` remain useful plotting controls, while ordinary DFT-DOS concepts such as +Fermi-level alignment, HOCO alignment, or total DOS are usually not the relevant analysis target. +See [&XAS%PRINT%PDOS](#CP2K_INPUT.FORCE_EVAL.DFT.XAS.PRINT.PDOS) and [&XAS_TDP%PRINT%PDOS](#CP2K_INPUT.FORCE_EVAL.DFT.XAS_TDP.PRINT.PDOS) for the corresponding input sections. diff --git a/src/input_cp2k_print_dft.F b/src/input_cp2k_print_dft.F index 9e41eec4f3..304ca53990 100644 --- a/src/input_cp2k_print_dft.F +++ b/src/input_cp2k_print_dft.F @@ -2687,9 +2687,9 @@ CONTAINS LOGICAL, INTENT(IN) :: xas_mode TYPE(keyword_type), POINTER :: keyword - TYPE(section_type), POINTER :: subsection + TYPE(section_type), POINTER :: broaden_subsection, pdos_section, subsection - NULLIFY (subsection) + NULLIFY (broaden_subsection, pdos_section, subsection) NULLIFY (keyword) CALL keyword_create(keyword, __LOCATION__, name="APPEND", & @@ -2699,70 +2699,18 @@ CONTAINS CALL section_add_keyword(print_key, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="DELTA_E", & - description="Energy spacing of the broadened DOS/PDOS output grid.", & - usage="DELTA_E 0.0005", type_of_var=real_t, default_r_val=0.001_dp) - CALL section_add_keyword(print_key, keyword) - CALL keyword_release(keyword) - - CALL keyword_create(keyword, __LOCATION__, name="ENERGY_UNIT", & - description="Energy unit used for the printed DOS/PDOS energy axis. "// & - "For broadened spectra, intensities are converted consistently to the selected energy unit.", & - usage="ENERGY_UNIT HARTREE", type_of_var=enum_t, & - enum_c_vals=s2a("HARTREE", "EV"), & - enum_i_vals=[1, 2], & - enum_desc=s2a("Print energies in Hartree (a.u.).", & - "Print energies in electronvolt."), & - default_i_val=1) - CALL section_add_keyword(print_key, keyword) - CALL keyword_release(keyword) - - CALL keyword_create(keyword, __LOCATION__, name="ENERGY_ZERO", & - description="Reference energy used for the printed DOS/PDOS energy axis. "// & - "With AUTO, the Fermi energy is used if smearing is enabled or fractional "// & - "occupations are found; otherwise the highest occupied crystal orbital is used.", & - usage="ENERGY_ZERO AUTO", type_of_var=enum_t, & - enum_c_vals=s2a("AUTO", "ABSOLUTE", "FERMI", "HOCO"), & - enum_i_vals=[1, 2, 3, 4], & - enum_desc=s2a("Choose FERMI for smeared or fractionally occupied systems, otherwise HOCO.", & - "Print absolute orbital energies.", & - "Shift orbital energies by the Fermi energy.", & - "Shift orbital energies by the highest occupied crystal orbital."), & - default_i_val=1) - CALL section_add_keyword(print_key, keyword) - CALL keyword_release(keyword) - - CALL keyword_create(keyword, __LOCATION__, name="BROADEN_TYPE", & - description="Type of broadening function used to produce the DOS/PDOS curve.", & - usage="BROADEN_TYPE GAUSSIAN", type_of_var=enum_t, & - enum_c_vals=s2a("GAUSSIAN", "LORENTZIAN", "PSEUDO_VOIGT"), & - enum_i_vals=[1, 2, 3], & - enum_desc=s2a("Use a Gaussian broadening function.", & - "Use a Lorentzian broadening function.", & - "Use a pseudo-Voigt mixture of Lorentzian and Gaussian functions."), & - default_i_val=1) - CALL section_add_keyword(print_key, keyword) - CALL keyword_release(keyword) - - CALL keyword_create(keyword, __LOCATION__, name="BROADEN_WIDTH", & - description="Full width at half maximum (FWHM) of the DOS/PDOS broadening function.", & - usage="BROADEN_WIDTH [eV] 0.1", type_of_var=real_t, & - default_r_val=cp_unit_to_cp2k(value=0.1_dp, unit_str="eV"), unit_str="eV") - CALL section_add_keyword(print_key, keyword) - CALL keyword_release(keyword) - - CALL keyword_create(keyword, __LOCATION__, name="VOIGT_MIXING", & - description="Lorentzian fraction of the pseudo-Voigt broadening function.", & - usage="VOIGT_MIXING 0.5", type_of_var=real_t, default_r_val=0.5_dp) - CALL section_add_keyword(print_key, keyword) - CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="NDIGITS", & description="Specify the number of digits used to print DOS/PDOS values.", & default_i_val=6) CALL section_add_keyword(print_key, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="DELTA_E", & + description="Energy spacing of the DOS/PDOS output grid.", & + usage="DELTA_E 0.0005", type_of_var=real_t, default_r_val=0.001_dp) + CALL section_add_keyword(print_key, keyword) + CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="NLUMO", & description="Number of unoccupied orbitals to include in the DOS/PDOS (-1=all). "// & "For OT calculations, the requested virtual orbitals are generated after SCF using the "// & @@ -2773,29 +2721,80 @@ CONTAINS CALL keyword_release(keyword) IF (.NOT. xas_mode) THEN - CALL section_create(subsection, __LOCATION__, name="PDOS", & + CALL section_create(pdos_section, __LOCATION__, name="PDOS", & description="Controls the printing of kind-resolved projected DOS.", & - n_keywords=2, n_subsections=0) - CALL keyword_create(keyword, __LOCATION__, name="COMPONENTS", & - description="Print out PDOS distinguishing all angular momentum components.", & - usage="COMPONENTS", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) - CALL section_add_keyword(subsection, keyword) - CALL keyword_release(keyword) - - CALL keyword_create(keyword, __LOCATION__, name="RAW", & - description="Print MO-resolved projected DOS weights to *.pdos_raw.", & - usage="RAW", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) - CALL section_add_keyword(subsection, keyword) - CALL keyword_release(keyword) - CALL section_add_subsection(print_key, subsection) - CALL section_release(subsection) + n_keywords=1, n_subsections=0) ELSE - CALL keyword_create(keyword, __LOCATION__, name="COMPONENTS", & - description="Print out PDOS distinguishing all angular momentum components.", & - usage="COMPONENTS", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) - CALL section_add_keyword(print_key, keyword) - CALL keyword_release(keyword) + pdos_section => print_key END IF + CALL keyword_create(keyword, __LOCATION__, name="COMPONENTS", & + description="Print out PDOS distinguishing all angular momentum components.", & + usage="COMPONENTS", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(pdos_section, keyword) + CALL keyword_release(keyword) + IF (.NOT. xas_mode) THEN + CALL section_add_subsection(print_key, pdos_section) + CALL section_release(pdos_section) + END IF + + CALL section_create(subsection, __LOCATION__, name="CURVE", & + description="Controls the printing of broadened DOS/PDOS curves.", & + n_keywords=2, n_subsections=1, repeats=.FALSE.) + CALL keyword_create(keyword, __LOCATION__, name="ENERGY_UNIT", & + description="Energy unit used for the printed DOS/PDOS energy axis. "// & + "Intensities are converted consistently to the selected energy unit.", & + usage="ENERGY_UNIT HARTREE", type_of_var=enum_t, & + enum_c_vals=s2a("HARTREE", "EV"), & + enum_i_vals=[1, 2], & + enum_desc=s2a("Print energies in Hartree (a.u.).", & + "Print energies in electronvolt."), & + default_i_val=1) + CALL section_add_keyword(subsection, keyword) + CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="ENERGY_ZERO", & + description="Reference energy used for the printed DOS/PDOS energy axis. "// & + "With AUTO, the Fermi energy is used if smearing is enabled or fractional "// & + "occupations are found; otherwise the highest occupied crystal orbital is used.", & + usage="ENERGY_ZERO AUTO", type_of_var=enum_t, & + enum_c_vals=s2a("AUTO", "ABSOLUTE", "FERMI", "HOCO"), & + enum_i_vals=[1, 2, 3, 4], & + enum_desc=s2a("Choose FERMI for smeared or fractionally occupied systems, otherwise HOCO.", & + "Print absolute orbital energies.", & + "Shift orbital energies by the Fermi energy.", & + "Shift orbital energies by the highest occupied crystal orbital."), & + default_i_val=1) + CALL section_add_keyword(subsection, keyword) + CALL keyword_release(keyword) + + CALL section_create(broaden_subsection, __LOCATION__, name="BROADEN", & + description="Controls the line shape used for broadened DOS/PDOS curves.", & + n_keywords=3, n_subsections=0, repeats=.FALSE.) + CALL keyword_create(keyword, __LOCATION__, name="TYPE", & + description="Type of broadening function used to produce the DOS/PDOS curve.", & + usage="TYPE GAUSSIAN", type_of_var=enum_t, & + enum_c_vals=s2a("GAUSSIAN", "LORENTZIAN", "PSEUDO_VOIGT"), & + enum_i_vals=[1, 2, 3], & + enum_desc=s2a("Use a Gaussian broadening function.", & + "Use a Lorentzian broadening function.", & + "Use a pseudo-Voigt mixture of Lorentzian and Gaussian functions."), & + default_i_val=1) + CALL section_add_keyword(broaden_subsection, keyword) + CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="WIDTH", & + description="Full width at half maximum (FWHM) of the DOS/PDOS broadening function.", & + usage="WIDTH [eV] 0.1", type_of_var=real_t, & + default_r_val=cp_unit_to_cp2k(value=0.1_dp, unit_str="eV"), unit_str="eV") + CALL section_add_keyword(broaden_subsection, keyword) + CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="VOIGT_MIXING", & + description="Lorentzian fraction of the pseudo-Voigt broadening function.", & + usage="VOIGT_MIXING 0.5", type_of_var=real_t, default_r_val=0.5_dp) + CALL section_add_keyword(broaden_subsection, keyword) + CALL keyword_release(keyword) + CALL section_add_subsection(subsection, broaden_subsection) + CALL section_release(broaden_subsection) + CALL section_add_subsection(print_key, subsection) + CALL section_release(subsection) CALL keyword_create(keyword, __LOCATION__, name="OUT_EACH_STATE", & variants=["OUT_EACH_MO"], & diff --git a/src/qs_dos.F b/src/qs_dos.F index 651c9e1813..a8a13cdf28 100644 --- a/src/qs_dos.F +++ b/src/qs_dos.F @@ -57,18 +57,19 @@ CONTAINS !> \param dft_section ... !> \param unoccupied_evals ... !> \param smearing_enabled ... +!> \param write_curve_output ... !> \date 26.02.2008 !> \par History: !> \author JGH !> \version 1.0 ! ************************************************************************************************** - SUBROUTINE calculate_dos(mos, dft_section, unoccupied_evals, smearing_enabled) + SUBROUTINE calculate_dos(mos, dft_section, unoccupied_evals, smearing_enabled, write_curve_output) TYPE(mo_set_type), DIMENSION(:), POINTER :: mos TYPE(section_vals_type), POINTER :: dft_section TYPE(cp_1d_r_p_type), DIMENSION(:), OPTIONAL, & POINTER :: unoccupied_evals - LOGICAL, INTENT(IN), OPTIONAL :: smearing_enabled + LOGICAL, INTENT(IN), OPTIONAL :: smearing_enabled, write_curve_output CHARACTER(len=*), PARAMETER :: routineN = 'calculate_dos' @@ -106,18 +107,32 @@ CONTAINS CALL section_vals_val_get(dft_section, "PRINT%DOS%DELTA_E", r_val=de) CALL section_vals_val_get(dft_section, "PRINT%DOS%APPEND", l_val=append) CALL section_vals_val_get(dft_section, "PRINT%DOS%NDIGITS", i_val=ndigits) - CALL section_vals_val_get(dft_section, "PRINT%DOS%ENERGY_UNIT", i_val=energy_unit) - CALL section_vals_val_get(dft_section, "PRINT%DOS%ENERGY_ZERO", i_val=energy_zero) - CALL section_vals_val_get(dft_section, "PRINT%DOS%BROADEN_TYPE", i_val=broaden_type) - CALL section_vals_val_get(dft_section, "PRINT%DOS%BROADEN_WIDTH", r_val=broaden_width) - CALL section_vals_val_get(dft_section, "PRINT%DOS%VOIGT_MIXING", r_val=voigt_mixing) + CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%ENERGY_UNIT", i_val=energy_unit) + CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%ENERGY_ZERO", i_val=energy_zero) + CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%BROADEN%TYPE", i_val=broaden_type) + CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%BROADEN%WIDTH", r_val=broaden_width) + CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%BROADEN%VOIGT_MIXING", r_val=voigt_mixing) IF (append .AND. iterstep > 1) THEN my_pos = "APPEND" ELSE my_pos = "REWIND" END IF ndigits = MIN(MAX(ndigits, 1), 10) - do_broaden = (broaden_width > 0.0_dp) + IF (PRESENT(write_curve_output)) THEN + IF (write_curve_output .AND. de <= 0.0_dp) THEN + CPWARN("Broadened DOS output requires DELTA_E > 0 and will be skipped") + CALL timestop(handle) + RETURN + END IF + IF (write_curve_output .AND. broaden_width <= 0.0_dp) THEN + CPWARN("Broadened DOS output requires a finite WIDTH and will be skipped") + CALL timestop(handle) + RETURN + END IF + END IF + do_broaden = .FALSE. + IF (PRESENT(write_curve_output)) do_broaden = write_curve_output + do_broaden = do_broaden .AND. (broaden_width > 0.0_dp) IF (do_broaden) de = MAX(de, 0.00001_dp) emin = 1.e10_dp @@ -239,21 +254,28 @@ CONTAINS CASE (dos_energy_zero_absolute) energy_ref(:) = 0.0_dp CASE (dos_energy_zero_hoco) - energy_ref(:) = hoco(:) + energy_ref(:) = MAXVAL(hoco(1:nspins)) CASE DEFAULT - energy_ref(:) = e_fermi(:) + energy_ref(:) = MAXVAL(e_fermi(1:nspins)) END SELECT - energy_factor = dos_energy_scale(energy_unit) - density_factor = dos_density_scale(energy_unit) - energy_label = dos_energy_label(energy_unit) + IF (.NOT. do_broaden) energy_ref(:) = 0.0_dp + energy_factor = MERGE(dos_energy_scale(energy_unit), 1.0_dp, do_broaden) + density_factor = MERGE(dos_density_scale(energy_unit), 1.0_dp, do_broaden) + energy_label = dos_energy_label(MERGE(energy_unit, 1, do_broaden)) zero_label = dos_energy_zero_label(resolved_energy_zero) IF (energy_zero == dos_energy_zero_auto) zero_label = "AUTO -> "//TRIM(zero_label) ev_factor = dos_energy_scale(dos_energy_unit_ev) my_act = "WRITE" - iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%DOS", & - extension=".dos", file_position=my_pos, file_action=my_act, & - file_form="FORMATTED") + IF (do_broaden) THEN + iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%DOS", & + extension=".dos", file_position=my_pos, file_action=my_act, & + file_form="FORMATTED", middle_name="curve") + ELSE + iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%DOS", & + extension=".dos", file_position=my_pos, file_action=my_act, & + file_form="FORMATTED") + END IF IF (iw > 0) THEN WRITE (UNIT=iw, FMT="(A,I0)") "# DOS at iteration step i = ", iterstep IF (nspins == 2) THEN @@ -261,33 +283,48 @@ CONTAINS "# E(Fermi) = ", e_fermi(1:2), " a.u. = ", e_fermi(1:2)*ev_factor, " eV" WRITE (UNIT=iw, FMT="(A,2F12.6,A,2F12.6,A)") & "# E(HOCO) = ", hoco(1:2), " a.u. = ", hoco(1:2)*ev_factor, " eV" - WRITE (UNIT=iw, FMT="(A,A)") "# Energy zero: ", TRIM(zero_label) - IF (do_broaden) CALL write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing) - WRITE (UNIT=iw, FMT="(T2,A,A,A)") "# "//TRIM(energy_label)//" Alpha_Density Occupation", & - " "//TRIM(energy_label), " Beta_Density Occupation" - ! (2(F15.8,2F20.ndigits)) - WRITE (UNIT=fmtstr_data, FMT="(A,I0,A)") "(2(F15.8,2F20.", ndigits, "))" + IF (do_broaden) THEN + WRITE (UNIT=iw, FMT="(A,A)") "# Energy zero: ", TRIM(zero_label) + CALL write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing) + END IF + IF (do_broaden .OR. de > 0.0_dp) THEN + WRITE (UNIT=iw, FMT="(A,A,A)") "# "//TRIM(energy_label)//" Alpha_Density Occupation", & + " Beta_Density Occupation" + WRITE (UNIT=fmtstr_data, FMT="(A,I0,A)") "(F15.8,4F20.", ndigits, ")" + ELSE + WRITE (UNIT=iw, FMT="(A,A,A)") "# "//TRIM(energy_label)//" Alpha_Density Occupation", & + " "//TRIM(energy_label), " Beta_Density Occupation" + WRITE (UNIT=fmtstr_data, FMT="(A,I0,A)") "(2(F15.8,2F20.", ndigits, "))" + END IF ELSE WRITE (UNIT=iw, FMT="(A,F12.6,A,F12.6,A)") & "# E(Fermi) = ", e_fermi(1), " a.u. = ", e_fermi(1)*ev_factor, " eV" WRITE (UNIT=iw, FMT="(A,F12.6,A,F12.6,A)") & "# E(HOCO) = ", hoco(1), " a.u. = ", hoco(1)*ev_factor, " eV" - WRITE (UNIT=iw, FMT="(A,A)") "# Energy zero: ", TRIM(zero_label) - IF (do_broaden) CALL write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing) + IF (do_broaden) THEN + WRITE (UNIT=iw, FMT="(A,A)") "# Energy zero: ", TRIM(zero_label) + CALL write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing) + END IF WRITE (UNIT=iw, FMT="(A,A)") "# "//TRIM(energy_label), " Density Occupation" ! (F15.8,2F20.ndigits) WRITE (UNIT=fmtstr_data, FMT="(A,I0,A)") "(F15.8,2F20.", ndigits, ")" END IF DO i = 1, nhist IF (nspins == 2) THEN - e1 = (ehist(i, 1) - energy_ref(1))*energy_factor - e2 = (ehist(i, 2) - energy_ref(2))*energy_factor - ! fmtstr_data == "(2(F15.8,2F20.xx))" IF (do_broaden) THEN - WRITE (UNIT=iw, FMT=fmtstr_data) e1, hist(i, 1)*density_factor, & - occval(i, 1)*density_factor, e2, hist(i, 2)*density_factor, & + eval = (ehist(i, 1) - energy_ref(1))*energy_factor + WRITE (UNIT=iw, FMT=fmtstr_data) eval, hist(i, 1)*density_factor, & + occval(i, 1)*density_factor, hist(i, 2)*density_factor, & occval(i, 2)*density_factor + ELSE IF (de > 0.0_dp) THEN + IF (hist(i, 1) == 0.0_dp .AND. occval(i, 1) == 0.0_dp .AND. & + hist(i, 2) == 0.0_dp .AND. occval(i, 2) == 0.0_dp) CYCLE + eval = ehist(i, 1) + WRITE (UNIT=iw, FMT=fmtstr_data) eval, hist(i, 1), occval(i, 1), & + hist(i, 2), occval(i, 2) ELSE + e1 = ehist(i, 1) + e2 = ehist(i, 2) WRITE (UNIT=iw, FMT=fmtstr_data) e1, hist(i, 1), occval(i, 1), & e2, hist(i, 2), occval(i, 2) END IF @@ -300,6 +337,7 @@ CONTAINS ELSE out_density = hist(i, 1) out_occ = occval(i, 1) + IF (out_density == 0.0_dp .AND. out_occ == 0.0_dp) CYCLE END IF WRITE (UNIT=iw, FMT=fmtstr_data) eval, out_density, out_occ END IF @@ -316,15 +354,17 @@ CONTAINS !> \brief Compute and write density of states (kpoints) !> \param qs_env ... !> \param dft_section ... +!> \param write_curve_output ... !> \date 26.02.2008 !> \par History: !> \author JGH !> \version 1.0 ! ************************************************************************************************** - SUBROUTINE calculate_dos_kp(qs_env, dft_section) + SUBROUTINE calculate_dos_kp(qs_env, dft_section, write_curve_output) TYPE(qs_environment_type), POINTER :: qs_env TYPE(section_vals_type), POINTER :: dft_section + LOGICAL, INTENT(IN), OPTIONAL :: write_curve_output CHARACTER(len=*), PARAMETER :: routineN = 'calculate_dos_kp' @@ -388,20 +428,41 @@ CONTAINS CALL section_vals_val_get(dft_section, "PRINT%DOS%DELTA_E", r_val=de) CALL section_vals_val_get(dft_section, "PRINT%DOS%APPEND", l_val=append) CALL section_vals_val_get(dft_section, "PRINT%DOS%NDIGITS", i_val=ndigits) - CALL section_vals_val_get(dft_section, "PRINT%DOS%ENERGY_UNIT", i_val=energy_unit) - CALL section_vals_val_get(dft_section, "PRINT%DOS%ENERGY_ZERO", i_val=energy_zero) - CALL section_vals_val_get(dft_section, "PRINT%DOS%BROADEN_TYPE", i_val=broaden_type) - CALL section_vals_val_get(dft_section, "PRINT%DOS%BROADEN_WIDTH", r_val=broaden_width) - CALL section_vals_val_get(dft_section, "PRINT%DOS%VOIGT_MIXING", r_val=voigt_mixing) - ! ensure a lower value for the DOS grid width - de = MAX(de, 0.00001_dp) + CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%ENERGY_UNIT", i_val=energy_unit) + CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%ENERGY_ZERO", i_val=energy_zero) + CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%BROADEN%TYPE", i_val=broaden_type) + CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%BROADEN%WIDTH", r_val=broaden_width) + CALL section_vals_val_get(dft_section, "PRINT%DOS%CURVE%BROADEN%VOIGT_MIXING", r_val=voigt_mixing) IF (append .AND. iterstep > 1) THEN my_pos = "APPEND" ELSE my_pos = "REWIND" END IF ndigits = MIN(MAX(ndigits, 1), 10) - do_broaden = (broaden_width > 0.0_dp) + IF (PRESENT(write_curve_output)) THEN + IF (write_curve_output .AND. de <= 0.0_dp) THEN + CPWARN("Broadened k-point DOS output requires DELTA_E > 0 and will be skipped") + CALL timestop(handle) + IF (explicit) CALL kpoint_release(kpoints) + RETURN + END IF + IF (write_curve_output .AND. broaden_width <= 0.0_dp) THEN + CPWARN("Broadened k-point DOS output requires a finite WIDTH and will be skipped") + CALL timestop(handle) + IF (explicit) CALL kpoint_release(kpoints) + RETURN + END IF + ELSE IF (de <= 0.0_dp) THEN + CPWARN("K-point DOS output requires DELTA_E > 0 and will be skipped") + CALL timestop(handle) + IF (explicit) CALL kpoint_release(kpoints) + RETURN + END IF + ! ensure a lower value for the DOS grid width + de = MAX(de, 0.00001_dp) + do_broaden = .FALSE. + IF (PRESENT(write_curve_output)) do_broaden = write_curve_output + do_broaden = do_broaden .AND. (broaden_width > 0.0_dp) CALL get_qs_env(qs_env, dft_control=dft_control) nspins = dft_control%nspins @@ -500,21 +561,28 @@ CONTAINS CASE (dos_energy_zero_absolute) energy_ref(:) = 0.0_dp CASE (dos_energy_zero_hoco) - energy_ref(:) = hoco(:) + energy_ref(:) = MAXVAL(hoco(1:nspins)) CASE DEFAULT - energy_ref(:) = e_fermi(:) + energy_ref(:) = MAXVAL(e_fermi(1:nspins)) END SELECT - energy_factor = dos_energy_scale(energy_unit) - density_factor = dos_density_scale(energy_unit) - energy_label = dos_energy_label(energy_unit) + IF (.NOT. do_broaden) energy_ref(:) = 0.0_dp + energy_factor = MERGE(dos_energy_scale(energy_unit), 1.0_dp, do_broaden) + density_factor = MERGE(dos_density_scale(energy_unit), 1.0_dp, do_broaden) + energy_label = dos_energy_label(MERGE(energy_unit, 1, do_broaden)) zero_label = dos_energy_zero_label(resolved_energy_zero) IF (energy_zero == dos_energy_zero_auto) zero_label = "AUTO -> "//TRIM(zero_label) ev_factor = dos_energy_scale(dos_energy_unit_ev) my_act = "WRITE" - iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%DOS", & - extension=".dos", file_position=my_pos, file_action=my_act, & - file_form="FORMATTED") + IF (do_broaden) THEN + iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%DOS", & + extension=".dos", file_position=my_pos, file_action=my_act, & + file_form="FORMATTED", middle_name="curve") + ELSE + iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%DOS", & + extension=".dos", file_position=my_pos, file_action=my_act, & + file_form="FORMATTED") + END IF IF (iw > 0) THEN WRITE (UNIT=iw, FMT="(A,I0)") "# DOS at iteration step i = ", iterstep IF (nspins == 2) THEN @@ -522,8 +590,10 @@ CONTAINS "# E(Fermi) = ", e_fermi(1:2), " a.u. = ", e_fermi(1:2)*ev_factor, " eV" WRITE (UNIT=iw, FMT="(A,2F12.6,A,2F12.6,A)") & "# E(HOCO) = ", hoco(1:2), " a.u. = ", hoco(1:2)*ev_factor, " eV" - WRITE (UNIT=iw, FMT="(A,A)") "# Energy zero: ", TRIM(zero_label) - IF (do_broaden) CALL write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing) + IF (do_broaden) THEN + WRITE (UNIT=iw, FMT="(A,A)") "# Energy zero: ", TRIM(zero_label) + CALL write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing) + END IF WRITE (UNIT=iw, FMT="(A,A)") "# "//TRIM(energy_label)//" Alpha_Density Occupation", & " Beta_Density Occupation" ! (F15.8,4F20.ndigits) @@ -533,8 +603,10 @@ CONTAINS "# E(Fermi) = ", e_fermi(1), " a.u. = ", e_fermi(1)*ev_factor, " eV" WRITE (UNIT=iw, FMT="(A,F12.6,A,F12.6,A)") & "# E(HOCO) = ", hoco(1), " a.u. = ", hoco(1)*ev_factor, " eV" - WRITE (UNIT=iw, FMT="(A,A)") "# Energy zero: ", TRIM(zero_label) - IF (do_broaden) CALL write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing) + IF (do_broaden) THEN + WRITE (UNIT=iw, FMT="(A,A)") "# Energy zero: ", TRIM(zero_label) + CALL write_broadening_info(iw, broaden_type, broaden_width, voigt_mixing) + END IF WRITE (UNIT=iw, FMT="(A,A)") "# "//TRIM(energy_label), " Density Occupation" ! (F15.8,2F20.ndigits) WRITE (UNIT=fmtstr_data, FMT="(A,I0,A)") "(F15.8,2F20.", ndigits, ")" @@ -547,6 +619,8 @@ CONTAINS WRITE (UNIT=iw, FMT=fmtstr_data) eval, hist(i, 1)*density_factor, & occval(i, 1)*density_factor, hist(i, 2)*density_factor, occval(i, 2)*density_factor ELSE + IF (hist(i, 1) == 0.0_dp .AND. occval(i, 1) == 0.0_dp .AND. & + hist(i, 2) == 0.0_dp .AND. occval(i, 2) == 0.0_dp) CYCLE WRITE (UNIT=iw, FMT=fmtstr_data) eval, hist(i, 1), occval(i, 1), & hist(i, 2), occval(i, 2) END IF @@ -558,6 +632,7 @@ CONTAINS ELSE out_density = hist(i, 1) out_occ = occval(i, 1) + IF (out_density == 0.0_dp .AND. out_occ == 0.0_dp) CYCLE END IF WRITE (UNIT=iw, FMT=fmtstr_data) eval, out_density, out_occ END IF diff --git a/src/qs_dos_utils.F b/src/qs_dos_utils.F index 6a615a1032..556745e277 100644 --- a/src/qs_dos_utils.F +++ b/src/qs_dos_utils.F @@ -12,8 +12,7 @@ MODULE qs_dos_utils USE cp_units, ONLY: cp_unit_from_cp2k USE input_section_types, ONLY: section_vals_get,& section_vals_get_subs_vals,& - section_vals_type,& - section_vals_val_get + section_vals_type USE kinds, ONLY: dp USE mathconstants, ONLY: pi #include "./base/base_uses.f90" @@ -146,28 +145,29 @@ CONTAINS !> \param do_dos_output whether the DOS print key is active !> \param do_projected_dos whether any projected DOS output is requested !> \param do_pdos whether kind-resolved PDOS output is requested -!> \param do_pdos_raw whether raw kind-resolved PDOS output is requested +!> \param do_curve whether broadened curve output is requested ! ************************************************************************************************** - SUBROUTINE get_dos_pdos_flags(dos_section, do_dos_output, do_projected_dos, do_pdos, do_pdos_raw) + SUBROUTINE get_dos_pdos_flags(dos_section, do_dos_output, do_projected_dos, do_pdos, do_curve) TYPE(section_vals_type), POINTER :: dos_section LOGICAL, INTENT(IN) :: do_dos_output - LOGICAL, INTENT(OUT) :: do_projected_dos, do_pdos, do_pdos_raw + LOGICAL, INTENT(OUT) :: do_projected_dos, do_pdos, do_curve INTEGER :: nrep LOGICAL :: has_ldos, has_r_ldos - TYPE(section_vals_type), POINTER :: ldos_section, pdos_section + TYPE(section_vals_type), POINTER :: curve_section, ldos_section, pdos_section do_projected_dos = .FALSE. do_pdos = .FALSE. - do_pdos_raw = .FALSE. + do_curve = .FALSE. has_ldos = .FALSE. has_r_ldos = .FALSE. IF (do_dos_output) THEN pdos_section => section_vals_get_subs_vals(dos_section, "PDOS") CALL section_vals_get(pdos_section, explicit=do_pdos) - IF (do_pdos) CALL section_vals_val_get(pdos_section, "RAW", l_val=do_pdos_raw) + curve_section => section_vals_get_subs_vals(dos_section, "CURVE") + CALL section_vals_get(curve_section, explicit=do_curve) ldos_section => section_vals_get_subs_vals(dos_section, "LDOS") CALL section_vals_get(ldos_section, n_repetition=nrep) has_ldos = (nrep > 0) diff --git a/src/qs_pdos.F b/src/qs_pdos.F index 301f7a1fb0..5a16f5b080 100644 --- a/src/qs_pdos.F +++ b/src/qs_pdos.F @@ -139,7 +139,7 @@ CONTAINS !> \param unoccupied_evals ... !> \param pdos_print_key ... !> \param write_pdos ... -!> \param write_pdos_raw ... +!> \param write_pdos_curve ... !> \date 26.02.2008 !> \par History: !> - Added optional external matrix_shalf to avoid recomputing it (A. Bussy, 09.2019) @@ -151,7 +151,7 @@ CONTAINS ! ************************************************************************************************** SUBROUTINE calculate_projected_dos(mo_set, atomic_kind_set, qs_kind_set, particle_set, qs_env, & dft_section, ispin, xas_mittle, external_matrix_shalf, & - unoccupied_orbs, unoccupied_evals, pdos_print_key, write_pdos, write_pdos_raw) + unoccupied_orbs, unoccupied_evals, pdos_print_key, write_pdos, write_pdos_curve) TYPE(mo_set_type), INTENT(IN) :: mo_set TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set @@ -165,7 +165,7 @@ CONTAINS TYPE(cp_fm_type), INTENT(IN), OPTIONAL, TARGET :: external_matrix_shalf, unoccupied_orbs TYPE(cp_1d_r_p_type), INTENT(IN), OPTIONAL, TARGET :: unoccupied_evals CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: pdos_print_key - LOGICAL, INTENT(IN), OPTIONAL :: write_pdos, write_pdos_raw + LOGICAL, INTENT(IN), OPTIONAL :: write_pdos, write_pdos_curve CHARACTER(len=*), PARAMETER :: routineN = 'calculate_projected_dos' @@ -178,21 +178,20 @@ CONTAINS CHARACTER(LEN=default_string_length), & ALLOCATABLE, DIMENSION(:) :: ldos_index, r_ldos_index INTEGER :: broaden_type, energy_unit, energy_zero, handle, homo, i, iatom, ikind, il, ildos, & - im, imo, in_x, in_y, in_z, ir, irow, iset, isgf, ishell, iso, iterstep, iw, j, jx, jy, & - jz, k, lcomponent, lshell, maxl, maxlgto, my_spin, n_dependent, n_r_ldos, n_rep, nao, & - natom, ncol_global, ndigits, nkind, nldos, nmo, np_tot, npoints, nrow_global, nset, nsgf, & - nvirt, out_each, output_unit, resolved_energy_zero + im, imo, imo_ref, in_x, in_y, in_z, ir, irow, iset, isgf, ishell, iso, ispin_ref, & + iterstep, iw, j, jx, jy, jz, k, lcomponent, lshell, maxl, maxlgto, my_spin, n_dependent, & + n_r_ldos, n_rep, nao, natom, ncol_global, ndigits, nkind, nldos, nmo, nmo_ref, np_tot, & + npoints, nrow_global, nset, nsgf, nvirt, out_each, output_unit, resolved_energy_zero INTEGER, ALLOCATABLE, DIMENSION(:) :: firstrow INTEGER, DIMENSION(:), POINTER :: list, nshell INTEGER, DIMENSION(:, :), POINTER :: bo, l - LOGICAL :: append, calc_matsh, do_broaden, do_ldos, do_r_ldos, do_virt, & - fractional_occupation, ionode, separate_components, should_output, write_raw, & - write_standard + LOGICAL :: append, calc_matsh, do_curve, do_ldos, do_r_ldos, do_virt, fractional_occupation, & + ionode, separate_components, should_output, write_curve, write_pdos_file LOGICAL, DIMENSION(:, :), POINTER :: read_r - REAL(KIND=dp) :: broaden_width, de, dh(3, 3), dvol, e_fermi, energy_factor, energy_ref, & - ev_factor, hoco, r(3), r_vec(3), ratom(3), voigt_mixing - REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, evals_virt, & - occupation_numbers + REAL(KIND=dp) :: broaden_width, de, dh(3, 3), dvol, e_fermi, e_fermi_ref(2), energy_factor, & + energy_ref, ev_factor, hoco, hoco_ref(2), r(3), r_vec(3), ratom(3), voigt_mixing + REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, eval_ref, evals_virt, & + occ_ref, occupation_numbers REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: pdos_array TYPE(cell_type), POINTER :: cell @@ -205,6 +204,7 @@ CONTAINS TYPE(dft_control_type), POINTER :: dft_control TYPE(gto_basis_set_type), POINTER :: orb_basis_set TYPE(ldos_p_type), DIMENSION(:), POINTER :: ldos_p + TYPE(mo_set_type), DIMENSION(:), POINTER :: mos_ref TYPE(mp_para_env_type), POINTER :: para_env TYPE(pw_c1d_gs_type) :: wf_g TYPE(pw_env_type), POINTER :: pw_env @@ -212,18 +212,18 @@ CONTAINS TYPE(pw_pool_type), POINTER :: auxbas_pw_pool TYPE(pw_r3d_rs_type) :: wf_r TYPE(r_ldos_p_type), DIMENSION(:), POINTER :: r_ldos_p - TYPE(section_vals_type), POINTER :: ldos_section + TYPE(section_vals_type), POINTER :: curve_section, ldos_section - NULLIFY (logger) + NULLIFY (logger, mos_ref, eval_ref, occ_ref) logger => cp_get_default_logger() ionode = logger%para_env%is_source() my_print_key = "PRINT%PDOS" IF (PRESENT(pdos_print_key)) my_print_key = TRIM(pdos_print_key) - write_standard = .TRUE. - IF (PRESENT(write_pdos)) write_standard = write_pdos - write_raw = .TRUE. - IF (PRESENT(write_pdos_raw)) write_raw = write_pdos_raw - write_raw = write_raw .AND. write_standard + write_pdos_file = .TRUE. + IF (PRESENT(write_pdos)) write_pdos_file = write_pdos + curve_section => section_vals_get_subs_vals(dft_section, TRIM(my_print_key)//"%CURVE") + CALL section_vals_get(curve_section, explicit=write_curve) + IF (PRESENT(write_pdos_curve)) write_curve = write_pdos_curve should_output = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, & TRIM(my_print_key)), cp_p_file) output_unit = cp_logger_get_default_io_unit(logger) @@ -234,7 +234,7 @@ CONTAINS NULLIFY (context, s_matrix, orb_basis_set, para_env, pdos_array) NULLIFY (eigenvalues, fm_struct_tmp, mo_coeff, vecbuffer, mo_virt) - NULLIFY (ldos_section, list, cell, pw_env, auxbas_pw_pool, evals_virt) + NULLIFY (curve_section, ldos_section, list, cell, pw_env, auxbas_pw_pool, evals_virt) NULLIFY (occupation_numbers, ldos_p, r_ldos_p, dft_control, occupation_numbers) CALL timeset(routineN, handle) @@ -260,15 +260,23 @@ CONTAINS CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%OUT_EACH_MO", i_val=out_each) IF (out_each == -1) out_each = nao + 1 CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%DELTA_E", r_val=de) - CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%BROADEN_TYPE", i_val=broaden_type) - CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%BROADEN_WIDTH", r_val=broaden_width) - CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%VOIGT_MIXING", r_val=voigt_mixing) + CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%CURVE%BROADEN%TYPE", i_val=broaden_type) + CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%CURVE%BROADEN%WIDTH", r_val=broaden_width) + CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%CURVE%BROADEN%VOIGT_MIXING", r_val=voigt_mixing) CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%NDIGITS", i_val=ndigits) - CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%ENERGY_UNIT", i_val=energy_unit) - CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%ENERGY_ZERO", i_val=energy_zero) + CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%CURVE%ENERGY_UNIT", i_val=energy_unit) + CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%CURVE%ENERGY_ZERO", i_val=energy_zero) ndigits = MIN(MAX(ndigits, 1), 10) - do_broaden = (broaden_width > 0.0_dp) - IF (do_broaden) de = MAX(de, 0.00001_dp) + IF (write_curve .AND. de <= 0.0_dp) THEN + CPWARN("Broadened PDOS output requires DELTA_E > 0 and will be skipped") + write_curve = .FALSE. + END IF + IF (write_curve .AND. broaden_width <= 0.0_dp) THEN + CPWARN("Broadened PDOS output requires a finite WIDTH and will be skipped") + write_curve = .FALSE. + END IF + do_curve = write_curve .AND. (broaden_width > 0.0_dp) + IF (do_curve) de = MAX(de, 0.00001_dp) nvirt = 0 NULLIFY (evals_virt) IF (PRESENT(unoccupied_orbs) .AND. PRESENT(unoccupied_evals)) THEN @@ -561,14 +569,39 @@ CONTAINS 1.0e-8_dp) fractional_occupation = .TRUE. END DO IF (hoco < -0.5_dp*HUGE(0.0_dp)) hoco = e_fermi + IF (PRESENT(ispin) .AND. dft_control%nspins == 2) THEN + e_fermi_ref(:) = 0.0_dp + hoco_ref(:) = -HUGE(0.0_dp) + CALL get_qs_env(qs_env=qs_env, mos=mos_ref) + IF (ASSOCIATED(mos_ref)) THEN + DO ispin_ref = 1, dft_control%nspins + CALL get_mo_set(mo_set=mos_ref(ispin_ref), nmo=nmo_ref, mu=e_fermi_ref(ispin_ref)) + eval_ref => mos_ref(ispin_ref)%eigenvalues + occ_ref => mos_ref(ispin_ref)%occupation_numbers + DO imo_ref = 1, nmo_ref + IF (occ_ref(imo_ref) > 1.0e-10_dp) & + hoco_ref(ispin_ref) = MAX(hoco_ref(ispin_ref), eval_ref(imo_ref)) + IF (ABS(occ_ref(imo_ref) - REAL(NINT(occ_ref(imo_ref)), KIND=dp)) > & + 1.0e-8_dp) fractional_occupation = .TRUE. + END DO + IF (hoco_ref(ispin_ref) < -0.5_dp*HUGE(0.0_dp)) hoco_ref(ispin_ref) = e_fermi_ref(ispin_ref) + END DO + ELSE + e_fermi_ref(:) = e_fermi + hoco_ref(:) = hoco + END IF + ELSE + e_fermi_ref(:) = e_fermi + hoco_ref(:) = hoco + END IF resolved_energy_zero = dos_resolve_energy_zero(energy_zero, dft_control%smear, fractional_occupation) SELECT CASE (resolved_energy_zero) CASE (dos_energy_zero_absolute) energy_ref = 0.0_dp CASE (dos_energy_zero_hoco) - energy_ref = hoco + energy_ref = MAXVAL(hoco_ref(1:dft_control%nspins)) CASE DEFAULT - energy_ref = e_fermi + energy_ref = MAXVAL(e_fermi_ref(1:dft_control%nspins)) END SELECT energy_factor = dos_energy_scale(energy_unit) energy_label = dos_energy_label(energy_unit) @@ -732,7 +765,7 @@ CONTAINS my_pos = "REWIND" END IF my_act = "WRITE" - IF (write_standard .OR. write_raw) THEN + IF (write_pdos_file) THEN DO ikind = 1, nkind NULLIFY (orb_basis_set) @@ -755,16 +788,16 @@ CONTAINS my_spin = 1 END IF - IF (write_standard .AND. do_broaden) THEN + IF (write_pdos_file .AND. do_curve) THEN IF (separate_components) THEN - CALL write_broadened_pdos(logger, dft_section, TRIM(my_mittle), my_pos, my_act, iterstep, & + CALL write_broadened_pdos(logger, dft_section, TRIM(my_mittle)//"_curve", my_pos, my_act, iterstep, & e_fermi, hoco, energy_ref, TRIM(zero_label), & "Projected DOS for atomic kind "//TRIM(kind_name), & maxl, .TRUE., pdos_array(1:nsoset(maxl), ikind, :), & eigenvalues, nmo + nvirt, de, broaden_type, broaden_width, & voigt_mixing, ndigits, pdos_print_key=TRIM(my_print_key)) ELSE - CALL write_broadened_pdos(logger, dft_section, TRIM(my_mittle), my_pos, my_act, iterstep, & + CALL write_broadened_pdos(logger, dft_section, TRIM(my_mittle)//"_curve", my_pos, my_act, iterstep, & e_fermi, hoco, energy_ref, TRIM(zero_label), & "Projected DOS for atomic kind "//TRIM(kind_name), & maxl, .FALSE., pdos_array(0:maxl, ikind, :), & @@ -773,9 +806,9 @@ CONTAINS END IF END IF - IF (write_raw) THEN + IF (write_pdos_file) THEN iw = cp_print_key_unit_nr(logger, dft_section, TRIM(my_print_key), & - extension=".pdos_raw", file_position=my_pos, file_action=my_act, & + extension=".pdos", file_position=my_pos, file_action=my_act, & file_form="FORMATTED", middle_name=TRIM(my_mittle)) IF (iw > 0) THEN @@ -795,7 +828,6 @@ CONTAINS "# E(Fermi) = ", e_fermi, " a.u. = ", e_fermi*ev_factor, " eV" WRITE (UNIT=iw, FMT="(A,F12.6,A,F12.6,A)") & "# E(HOCO) = ", hoco, " a.u. = ", hoco*ev_factor, " eV" - WRITE (UNIT=iw, FMT="(A,A)") "# Energy zero: ", TRIM(zero_label) IF (separate_components) THEN ALLOCATE (tmp_str(0:0, 0:maxl, -maxl:maxl)) tmp_str = "" @@ -806,20 +838,20 @@ CONTAINS END DO WRITE (UNIT=iw, FMT=fmtstr2) & - "# MO "//TRIM(energy_label)//" Occupation", "Total", & + "# MO Energy[a.u.] Occupation", "Total", & ((TRIM(tmp_str(0, il, im)), im=-il, il), il=0, maxl) DO imo = 1, nmo + nvirt - WRITE (UNIT=iw, FMT=fmtstr1) imo, (eigenvalues(imo) - energy_ref)*energy_factor, & + WRITE (UNIT=iw, FMT=fmtstr1) imo, eigenvalues(imo), & occupation_numbers(imo), SUM(pdos_array(1:nsoset(maxl), ikind, imo)), & (pdos_array(lshell, ikind, imo), lshell=1, nsoset(maxl)) END DO DEALLOCATE (tmp_str) ELSE WRITE (UNIT=iw, FMT=fmtstr2) & - "# MO "//TRIM(energy_label)//" Occupation", "Total", & + "# MO Energy[a.u.] Occupation", "Total", & (TRIM(l_sym(il)), il=0, maxl) DO imo = 1, nmo + nvirt - WRITE (UNIT=iw, FMT=fmtstr1) imo, (eigenvalues(imo) - energy_ref)*energy_factor, & + WRITE (UNIT=iw, FMT=fmtstr1) imo, eigenvalues(imo), & occupation_numbers(imo), SUM(pdos_array(0:maxl, ikind, imo)), & (pdos_array(lshell, ikind, imo), lshell=0, maxl) END DO @@ -850,9 +882,9 @@ CONTAINS my_spin = 1 END IF - IF (do_broaden) THEN + IF (do_curve) THEN IF (ldos_p(ildos)%ldos%separate_components) THEN - CALL write_broadened_pdos(logger, dft_section, TRIM(my_mittle), my_pos, my_act, iterstep, & + CALL write_broadened_pdos(logger, dft_section, TRIM(my_mittle)//"_curve", my_pos, my_act, iterstep, & e_fermi, hoco, energy_ref, TRIM(zero_label), & "Projected DOS for atom list "//TRIM(ldos_index(ildos)), & ldos_p(ildos)%ldos%maxl, .TRUE., & @@ -860,7 +892,7 @@ CONTAINS eigenvalues, nmo + nvirt, de, broaden_type, broaden_width, & voigt_mixing, ndigits, pdos_print_key=TRIM(my_print_key)) ELSE - CALL write_broadened_pdos(logger, dft_section, TRIM(my_mittle), my_pos, my_act, iterstep, & + CALL write_broadened_pdos(logger, dft_section, TRIM(my_mittle)//"_curve", my_pos, my_act, iterstep, & e_fermi, hoco, energy_ref, TRIM(zero_label), & "Projected DOS for atom list "//TRIM(ldos_index(ildos)), & ldos_p(ildos)%ldos%maxl, .FALSE., & @@ -871,7 +903,7 @@ CONTAINS END IF iw = cp_print_key_unit_nr(logger, dft_section, TRIM(my_print_key), & - extension=".pdos_raw", file_position=my_pos, file_action=my_act, & + extension=".pdos", file_position=my_pos, file_action=my_act, & file_form="FORMATTED", middle_name=TRIM(my_mittle)) IF (iw > 0) THEN @@ -892,7 +924,6 @@ CONTAINS "# E(Fermi) = ", e_fermi, " a.u. = ", e_fermi*ev_factor, " eV" WRITE (UNIT=iw, FMT="(A,F12.6,A,F12.6,A)") & "# E(HOCO) = ", hoco, " a.u. = ", hoco*ev_factor, " eV" - WRITE (UNIT=iw, FMT="(A,A)") "# Energy zero: ", TRIM(zero_label) IF (ldos_p(ildos)%ldos%separate_components) THEN ALLOCATE (tmp_str(0:0, 0:ldos_p(ildos)%ldos%maxl, -ldos_p(ildos)%ldos%maxl:ldos_p(ildos)%ldos%maxl)) tmp_str = "" @@ -903,10 +934,10 @@ CONTAINS END DO WRITE (UNIT=iw, FMT=fmtstr2) & - "# MO "//TRIM(energy_label)//" Occupation", "Total", & + "# MO Energy[a.u.] Occupation", "Total", & ((TRIM(tmp_str(0, il, im)), im=-il, il), il=0, ldos_p(ildos)%ldos%maxl) DO imo = 1, nmo + nvirt - WRITE (UNIT=iw, FMT=fmtstr1) imo, (eigenvalues(imo) - energy_ref)*energy_factor, & + WRITE (UNIT=iw, FMT=fmtstr1) imo, eigenvalues(imo), & occupation_numbers(imo), & SUM(ldos_p(ildos)%ldos%pdos_array(1:nsoset(ldos_p(ildos)%ldos%maxl), imo)), & (ldos_p(ildos)%ldos%pdos_array(lshell, imo), & @@ -915,10 +946,10 @@ CONTAINS DEALLOCATE (tmp_str) ELSE WRITE (UNIT=iw, FMT=fmtstr2) & - "# MO "//TRIM(energy_label)//" Occupation", "Total", & + "# MO Energy[a.u.] Occupation", "Total", & (TRIM(l_sym(il)), il=0, ldos_p(ildos)%ldos%maxl) DO imo = 1, nmo + nvirt - WRITE (UNIT=iw, FMT=fmtstr1) imo, (eigenvalues(imo) - energy_ref)*energy_factor, & + WRITE (UNIT=iw, FMT=fmtstr1) imo, eigenvalues(imo), & occupation_numbers(imo), & SUM(ldos_p(ildos)%ldos%pdos_array(0:ldos_p(ildos)%ldos%maxl, imo)), & (ldos_p(ildos)%ldos%pdos_array(lshell, imo), lshell=0, ldos_p(ildos)%ldos%maxl) @@ -951,7 +982,7 @@ CONTAINS END IF iw = cp_print_key_unit_nr(logger, dft_section, TRIM(my_print_key), & - extension=".pdos_raw", file_position=my_pos, file_action=my_act, & + extension=".pdos", file_position=my_pos, file_action=my_act, & file_form="FORMATTED", middle_name=TRIM(my_mittle)) IF (iw > 0) THEN fmtstr1 = "(I8,2X,2F16.6, (2X,F16.8))" @@ -964,14 +995,13 @@ CONTAINS "# E(Fermi) = ", e_fermi, " a.u. = ", e_fermi*ev_factor, " eV" WRITE (UNIT=iw, FMT="(A,F12.6,A,F12.6,A)") & "# E(HOCO) = ", hoco, " a.u. = ", hoco*ev_factor, " eV" - WRITE (UNIT=iw, FMT="(A,A)") "# Energy zero: ", TRIM(zero_label) - WRITE (UNIT=iw, FMT="(A,A,A)") & - "# MO ", TRIM(energy_label), " Occupation LDOS" + WRITE (UNIT=iw, FMT="(A)") & + "# MO Energy[a.u.] Occupation LDOS" DO imo = 1, nmo + nvirt IF (r_ldos_p(ildos)%ldos%eval_range(1) <= eigenvalues(imo) .AND. & r_ldos_p(ildos)%ldos%eval_range(2) >= eigenvalues(imo)) THEN WRITE (UNIT=iw, FMT="(I8,2X,2F16.6,E20.10,E20.10)") imo, & - (eigenvalues(imo) - energy_ref)*energy_factor, occupation_numbers(imo), & + eigenvalues(imo), occupation_numbers(imo), & r_ldos_p(ildos)%ldos%pdos_array(imo), r_ldos_p(ildos)%ldos%pdos_array(imo)*np_tot END IF END DO @@ -1034,14 +1064,14 @@ CONTAINS !> \param dft_section ... !> \param pdos_print_key ... !> \param write_pdos ... -!> \param write_pdos_raw ... +!> \param write_pdos_curve ... ! ************************************************************************************************** - SUBROUTINE calculate_projected_dos_kp(qs_env, dft_section, pdos_print_key, write_pdos, write_pdos_raw) + SUBROUTINE calculate_projected_dos_kp(qs_env, dft_section, pdos_print_key, write_pdos, write_pdos_curve) TYPE(qs_environment_type), POINTER :: qs_env TYPE(section_vals_type), POINTER :: dft_section CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: pdos_print_key - LOGICAL, INTENT(IN), OPTIONAL :: write_pdos, write_pdos_raw + LOGICAL, INTENT(IN), OPTIONAL :: write_pdos, write_pdos_curve CHARACTER(len=*), PARAMETER :: routineN = 'calculate_projected_dos_kp' @@ -1055,7 +1085,7 @@ CONTAINS INTEGER, ALLOCATABLE, DIMENSION(:) :: ao_comp, ao_kind, ao_l, kind_maxl LOGICAL :: append, fractional_occupation, & separate_components, should_output, & - write_raw, write_standard + write_curve, write_pdos_file REAL(KIND=dp) :: broaden_cutoff, broaden_width, de, e1, & e2, e_fermi(2), emax, emin, & energy_ref(2), hoco(2), voigt_mixing, & @@ -1079,19 +1109,19 @@ CONTAINS TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set TYPE(qs_scf_env_type), POINTER :: scf_env - TYPE(section_vals_type), POINTER :: ldos_section + TYPE(section_vals_type), POINTER :: curve_section, ldos_section NULLIFY (logger, kpoints, dft_control, para_env, atomic_kind_set, qs_kind_set, particle_set) NULLIFY (matrix_s_kp, scf_env) - NULLIFY (kp, mo_set, eigenvalues, fm_struct_tmp, orb_basis_set, ldos_section) + NULLIFY (kp, mo_set, eigenvalues, fm_struct_tmp, orb_basis_set, curve_section, ldos_section) logger => cp_get_default_logger() my_print_key = "PRINT%PDOS" IF (PRESENT(pdos_print_key)) my_print_key = TRIM(pdos_print_key) - write_standard = .TRUE. - IF (PRESENT(write_pdos)) write_standard = write_pdos - write_raw = .FALSE. - IF (PRESENT(write_pdos_raw)) write_raw = write_pdos_raw - write_raw = write_raw .AND. write_standard + write_pdos_file = .TRUE. + IF (PRESENT(write_pdos)) write_pdos_file = write_pdos + curve_section => section_vals_get_subs_vals(dft_section, TRIM(my_print_key)//"%CURVE") + CALL section_vals_get(curve_section, explicit=write_curve) + IF (PRESENT(write_pdos_curve)) write_curve = write_pdos_curve should_output = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, & TRIM(my_print_key)), cp_p_file) output_unit = cp_logger_get_default_io_unit(logger) @@ -1134,12 +1164,17 @@ CONTAINS CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%COMPONENTS", l_val=separate_components) END IF CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%NDIGITS", i_val=ndigits) - CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%BROADEN_TYPE", i_val=broaden_type) - CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%BROADEN_WIDTH", r_val=broaden_width) - CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%VOIGT_MIXING", r_val=voigt_mixing) - CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%ENERGY_ZERO", i_val=energy_zero) + CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%CURVE%BROADEN%TYPE", i_val=broaden_type) + CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%CURVE%BROADEN%WIDTH", r_val=broaden_width) + CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%CURVE%BROADEN%VOIGT_MIXING", r_val=voigt_mixing) + CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%CURVE%ENERGY_ZERO", i_val=energy_zero) ndigits = MIN(MAX(ndigits, 1), 10) - de = MAX(de, 0.00001_dp) + IF (write_curve .AND. de <= 0.0_dp) THEN + CPWARN("Broadened k-point PDOS output requires DELTA_E > 0 and will be skipped") + CALL timestop(handle) + RETURN + END IF + IF (write_curve) de = MAX(de, 0.00001_dp) ldos_section => section_vals_get_subs_vals(dft_section, TRIM(my_print_key)//"%LDOS") CALL section_vals_get(ldos_section, n_repetition=nldos) @@ -1148,11 +1183,15 @@ CONTAINS IF (nldos > 0 .OR. n_r_ldos > 0) THEN CPWARN("LDOS/R_LDOS are not implemented for k-point PDOS and will be ignored") END IF - IF (write_raw) THEN - CPWARN("Raw k-point PDOS output is not implemented and will be ignored") + IF (write_pdos_file) THEN + CPWARN("State-resolved k-point PDOS output is not implemented yet") + END IF + IF (.NOT. write_curve) THEN + CALL timestop(handle) + RETURN END IF IF (broaden_width <= 0.0_dp) THEN - CPWARN("Broadened k-point PDOS output requires a finite BROADEN_WIDTH") + CPWARN("Broadened k-point PDOS output requires a finite WIDTH and will be skipped") CALL timestop(handle) RETURN END IF @@ -1211,9 +1250,9 @@ CONTAINS CASE (dos_energy_zero_absolute) energy_ref(:) = 0.0_dp CASE (dos_energy_zero_hoco) - energy_ref(:) = hoco(:) + energy_ref(:) = MAXVAL(hoco(1:nspins)) CASE DEFAULT - energy_ref(:) = e_fermi(:) + energy_ref(:) = MAXVAL(e_fermi(1:nspins)) END SELECT zero_label = dos_energy_zero_label(resolved_energy_zero) IF (energy_zero == dos_energy_zero_auto) zero_label = "AUTO -> "//TRIM(zero_label) @@ -1302,7 +1341,7 @@ CONTAINS my_act = "WRITE" spin(1) = "ALPHA" spin(2) = "BETA" - IF (write_standard) THEN + IF (write_pdos_file) THEN DO ikind = 1, nkind IF (kind_maxl(ikind) < 0) CYCLE CALL get_atomic_kind(atomic_kind_set(ikind), name=kind_name) @@ -1312,7 +1351,7 @@ CONTAINS ELSE my_mittle = "k"//TRIM(ADJUSTL(cp_to_string(ikind))) END IF - CALL write_broadened_pdos_curve(logger, dft_section, TRIM(my_mittle), my_pos, my_act, & + CALL write_broadened_pdos_curve(logger, dft_section, TRIM(my_mittle)//"_curve", my_pos, my_act, & iterstep, e_fermi(ispin), hoco(ispin), energy_ref(ispin), & TRIM(zero_label), & "K-point projected DOS for atomic kind "//TRIM(kind_name), & @@ -1473,7 +1512,7 @@ CONTAINS ELSE ncomponents = maxl + 1 END IF - CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%ENERGY_UNIT", i_val=energy_unit) + CALL section_vals_val_get(dft_section, TRIM(my_print_key)//"%CURVE%ENERGY_UNIT", i_val=energy_unit) energy_factor = dos_energy_scale(energy_unit) density_factor = dos_density_scale(energy_unit) energy_label = dos_energy_label(energy_unit) diff --git a/src/qs_scf_post_gpw.F b/src/qs_scf_post_gpw.F index ff7e14a629..16cd43613c 100644 --- a/src/qs_scf_post_gpw.F +++ b/src/qs_scf_post_gpw.F @@ -2143,8 +2143,8 @@ CONTAINS INTEGER :: handle, homo, ispin, nlumo_dos, & nlumo_molden, nlumo_required, nlumos, & nmo, output_unit - LOGICAL :: all_equal, defer_molden, do_dos, & - do_kpoints, do_pdos, do_pdos_raw, & + LOGICAL :: all_equal, defer_molden, do_curve, & + do_dos, do_kpoints, do_pdos, & do_projected_dos, explicit REAL(KIND=dp) :: maxocc, s_square, s_square_ideal, & total_abs_spin_dens, total_spin_dens @@ -2266,7 +2266,7 @@ CONTAINS dos_section => section_vals_get_subs_vals(dft_section, "PRINT%DOS") do_dos = BTEST(cp_print_key_should_output(logger%iter_info, dos_section), cp_p_file) - CALL get_dos_pdos_flags(dos_section, do_dos, do_projected_dos, do_pdos, do_pdos_raw) + CALL get_dos_pdos_flags(dos_section, do_dos, do_projected_dos, do_pdos, do_curve) ! For OT calculations, collect the largest request for additional unoccupied ! orbitals among DOS, PDOS, and Molden, and generate them only once. @@ -2339,13 +2339,18 @@ CONTAINS IF (do_dos) THEN IF (do_kpoints) THEN CALL calculate_dos_kp(qs_env, dft_section) + IF (do_curve) CALL calculate_dos_kp(qs_env, dft_section, write_curve_output=.TRUE.) ELSE CALL get_qs_env(qs_env, mos=mos) IF (ASSOCIATED(unoccupied_evals)) THEN CALL calculate_dos(mos, dft_section, unoccupied_evals=unoccupied_evals, & smearing_enabled=dft_control%smear) + IF (do_curve) CALL calculate_dos(mos, dft_section, unoccupied_evals=unoccupied_evals, & + smearing_enabled=dft_control%smear, write_curve_output=.TRUE.) ELSE CALL calculate_dos(mos, dft_section, smearing_enabled=dft_control%smear) + IF (do_curve) CALL calculate_dos(mos, dft_section, smearing_enabled=dft_control%smear, & + write_curve_output=.TRUE.) END IF END IF END IF @@ -2354,7 +2359,7 @@ CONTAINS IF (do_projected_dos) THEN IF (do_kpoints) THEN CALL calculate_projected_dos_kp(qs_env, dft_section, pdos_print_key="PRINT%DOS", & - write_pdos=do_pdos, write_pdos_raw=do_pdos_raw) + write_pdos=do_pdos, write_pdos_curve=do_curve) ELSE CALL get_qs_env(qs_env, & mos=mos, & @@ -2366,11 +2371,11 @@ CONTAINS qs_kind_set, particle_set, qs_env, dft_section, ispin=ispin, & unoccupied_orbs=unoccupied_orbs(ispin), & unoccupied_evals=unoccupied_evals(ispin), & - pdos_print_key="PRINT%DOS", write_pdos=do_pdos, write_pdos_raw=do_pdos_raw) + pdos_print_key="PRINT%DOS", write_pdos=do_pdos, write_pdos_curve=do_curve) ELSE CALL calculate_projected_dos(mos(ispin), atomic_kind_set, & qs_kind_set, particle_set, qs_env, dft_section, ispin=ispin, & - pdos_print_key="PRINT%DOS", write_pdos=do_pdos, write_pdos_raw=do_pdos_raw) + pdos_print_key="PRINT%DOS", write_pdos=do_pdos, write_pdos_curve=do_curve) END IF ELSE IF (ASSOCIATED(unoccupied_orbs)) THEN @@ -2378,11 +2383,11 @@ CONTAINS qs_kind_set, particle_set, qs_env, dft_section, & unoccupied_orbs=unoccupied_orbs(ispin), & unoccupied_evals=unoccupied_evals(ispin), & - pdos_print_key="PRINT%DOS", write_pdos=do_pdos, write_pdos_raw=do_pdos_raw) + pdos_print_key="PRINT%DOS", write_pdos=do_pdos, write_pdos_curve=do_curve) ELSE CALL calculate_projected_dos(mos(ispin), atomic_kind_set, & qs_kind_set, particle_set, qs_env, dft_section, & - pdos_print_key="PRINT%DOS", write_pdos=do_pdos, write_pdos_raw=do_pdos_raw) + pdos_print_key="PRINT%DOS", write_pdos=do_pdos, write_pdos_curve=do_curve) END IF END IF END DO diff --git a/src/qs_scf_post_tb.F b/src/qs_scf_post_tb.F index bacb1090e8..c8a4d4cd64 100644 --- a/src/qs_scf_post_tb.F +++ b/src/qs_scf_post_tb.F @@ -179,8 +179,8 @@ CONTAINS CHARACTER(LEN=default_string_length) :: aname INTEGER :: after, gfn_type, handle, homo, iat, iatom, ikind, img, ispin, iw, nat, natom, & nkind, nlumo_stm, nlumos, nspins, print_level, unit_nr - LOGICAL :: do_cube, do_dos, do_kpoints, do_pdos, do_pdos_raw, do_projected_dos, explicit, & - gfn0, has_homo, omit_headers, print_it, rebuild, vdip + LOGICAL :: do_cube, do_curve, do_dos, do_kpoints, do_pdos, do_projected_dos, explicit, gfn0, & + has_homo, omit_headers, print_it, rebuild, vdip REAL(KIND=dp) :: zeff REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: mcharge, zcharge REAL(KIND=dp), DIMENSION(2, 2) :: homo_lumo @@ -446,13 +446,16 @@ CONTAINS IF (.NOT. no_mos) THEN print_key => section_vals_get_subs_vals(print_section, "DOS") do_dos = BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file) - CALL get_dos_pdos_flags(print_key, do_dos, do_projected_dos, do_pdos, do_pdos_raw) + CALL get_dos_pdos_flags(print_key, do_dos, do_projected_dos, do_pdos, do_curve) IF (do_dos) THEN IF (do_kpoints) THEN CALL calculate_dos_kp(qs_env, dft_section) + IF (do_curve) CALL calculate_dos_kp(qs_env, dft_section, write_curve_output=.TRUE.) ELSE CALL get_qs_env(qs_env, mos=mos) CALL calculate_dos(mos, dft_section, smearing_enabled=dft_control%smear) + IF (do_curve) CALL calculate_dos(mos, dft_section, smearing_enabled=dft_control%smear, & + write_curve_output=.TRUE.) END IF END IF @@ -460,7 +463,7 @@ CONTAINS IF (do_projected_dos) THEN IF (do_kpoints) THEN CALL calculate_projected_dos_kp(qs_env, dft_section, pdos_print_key="PRINT%DOS", & - write_pdos=do_pdos, write_pdos_raw=do_pdos_raw) + write_pdos=do_pdos, write_pdos_curve=do_curve) ELSE CALL get_qs_env(qs_env, mos=mos, matrix_ks=ks_rmpv) DO ispin = 1, dft_control%nspins @@ -480,11 +483,11 @@ CONTAINS IF (dft_control%nspins == 2) THEN CALL calculate_projected_dos(mos(ispin), atomic_kind_set, & qs_kind_set, particle_set, qs_env, dft_section, ispin=ispin, & - pdos_print_key="PRINT%DOS", write_pdos=do_pdos, write_pdos_raw=do_pdos_raw) + pdos_print_key="PRINT%DOS", write_pdos=do_pdos, write_pdos_curve=do_curve) ELSE CALL calculate_projected_dos(mos(ispin), atomic_kind_set, & qs_kind_set, particle_set, qs_env, dft_section, & - pdos_print_key="PRINT%DOS", write_pdos=do_pdos, write_pdos_raw=do_pdos_raw) + pdos_print_key="PRINT%DOS", write_pdos=do_pdos, write_pdos_curve=do_curve) END IF END DO END IF