comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-21 06:25:15 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

(P)DOS: Restore default behavior of output data (#5564)

Browse source
This commit is contained in:
SY Wang 2026-07-13 22:22:13 +08:00 • committed by GitHub
parent d2d9c0e403
commit 2c97d88377
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
7 changed files with 412 additions and 290 deletions
Download patch file Download diff file Expand all files Collapse all files

117
docs/methods/electronic_structure/dos.md
View file

@ -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.

159
src/input_cp2k_print_dft.F
View file

@ -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"], &

175
src/qs_dos.F
View file

@ -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

16
src/qs_dos_utils.F
View file

@ -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)

199
src/qs_pdos.F
View file

@ -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)

21
src/qs_scf_post_gpw.F
View file

@ -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

15
src/qs_scf_post_tb.F
View file

@ -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