dos.F90

Go to the documentation of this file.

!! Copyright (C) 2017 N. Tancogne-Dejean
!! Copyright (C) 2025 N. Tancogne-Dejean
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module dos_oct_m
  use accel_oct_m
  use atomic_orbital_oct_m
  use box_oct_m
  use comm_oct_m
  use debug_oct_m
  use electron_space_oct_m
  use global_oct_m
  use hamiltonian_elec_oct_m
  use io_oct_m
  use io_function_oct_m
  use ions_oct_m
  use, intrinsic :: iso_fortran_env
  use kpoints_oct_m
  use lalg_basic_oct_m
  use math_oct_m
  use mesh_oct_m
  use messages_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use orbitalset_oct_m
  use orbitalset_utils_oct_m
  use output_low_oct_m
  use parser_oct_m
  use profiling_oct_m
  use simplex_oct_m, only: simplex_dos
  use smear_oct_m, only: smear_t, smear_delta_function, &
    smear_fermi_dirac, &
    smear_cold, &
    smear_methfessel_paxton, &
    smear_spline, &
    smear_gaussian, &
    smear_lorentzian
  use species_oct_m
  use states_abst_oct_m
  use states_elec_oct_m
  use submesh_oct_m
  use types_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use wfs_elec_oct_m
 
  implicit none
 
  private
  public ::                    &
    dos_t,                     &
    dos_init,                  &
    dos_write_dos,             &
    dos_write_jdos,            &
    dos_write_ldos
 
  type dos_t
    private
    real(real64)   :: emin
    real(real64)   :: emax
    integer :: epoints
    real(real64)   :: de
 
    real(real64)   :: gamma
 
    logical :: computepdos
 
    integer :: ldos_nenergies = -1
    real(real64), allocatable :: ldos_energies(:)
 
    integer(int64) :: method
    integer :: smear_func
  contains
    final :: dos_end
  end type dos_t
 
contains
 
  subroutine dos_init(this, namespace, st, kpoints)
    type(dos_t),         intent(out)   :: this
    type(namespace_t),   intent(in)    :: namespace
    type(states_elec_t), intent(in)    :: st
    type(kpoints_t),     intent(in)    :: kpoints
 
    real(real64) :: evalmin, evalmax, eextend
    integer :: npath, ie
    type(block_t) :: blk
 
    push_sub(dos_init)
 
    !The range of the dos is only calculated for physical points,
    !without the one from a k-point path
    npath = kpoints%nkpt_in_path()
    if (st%nik > npath) then
      evalmin = minval(st%eigenval(1:st%nst, 1:(st%nik-npath)))
      evalmax = maxval(st%eigenval(1:st%nst, 1:(st%nik-npath)))
    else !In case we only have a path, e.g., a bandstructure calculation
      evalmin = minval(st%eigenval(1:st%nst, 1:st%nik))
      evalmax = maxval(st%eigenval(1:st%nst, 1:st%nik))
    end if
    ! we extend the energy mesh by this amount
    eextend  = (evalmax - evalmin) / m_four
 
    !! The numbers 21 and 22 for tetrahedra options are reused in src/states/smear.F90
    !%Variable DOSMethod
    !%Type integer
    !%Default smear
    !%Section Output
    !%Description
    !% Selects the method that is used to calculate the DOS.
    !% Currently not used for PDOS and LDOS.
    !%Option smear 0
    !% Smearing with the smearing function selected by <tt>DOSSmearingFunction</tt>.
    !%Option tetrahedra 21
    !% Linear tetrahedron method as in P. E. Bloechl, et al., <i>Phys. Rev. B</i> <b>49</b>, 16223 (1994).
    !% Requires a regular Monkhorst-Pack generated by Octopus.
    !%Option tetrahedra_opt 22
    !% Improved tetrahedron method as in M. Kawamura, et al., <i>Phys. Rev. B</i> <b>89</b>, 094515 (2014).
    !% Requires a regular Monkhorst-Pack generated by Octopus.
    !%End
    call parse_variable(namespace, 'DOSMethod', option__dosmethod__smear, this%method)
 
    !%Variable DOSSmearingFunction
    !%Type integer
    !%Default lorentzian_smearing
    !%Section Output
    !%Description
    !% This is the function used to smear the energy levels in the DOS calculation.
    !%Option fermi_dirac 2
    !% Simple Fermi-Dirac distribution. In this case, <tt>DOSGamma</tt> has
    !% the meaning of an electronic temperature. DN Mermin, <i>Phys. Rev.</i> <b>137</b>, A1441 (1965).
    !%Option cold_smearing 3
    !% N Marzari, D Vanderbilt, A De Vita, and MC Payne, <i>Phys. Rev. Lett.</i> <b>82</b>, 3296 (1999).
    !%Option methfessel_paxton 4
    !% M Methfessel and AT Paxton, <i>Phys. Rev. B</i> <b>40</b>, 3616 (1989).
    !% The expansion order of the polynomial is fixed to 1.
    !% Occupations may be negative.
    !%Option spline_smearing 5
    !% Nearly identical to Gaussian smearing.
    !% JM Holender, MJ Gillan, MC Payne, and AD Simpson, <i>Phys. Rev. B</i> <b>52</b>, 967 (1995).
    !%Option gaussian_smearing 7
    !% Gaussian smearing.
    !%Option lorentzian_smearing 8
    !% Lorentzian smearing.
    !%End
    call parse_variable(namespace, 'DOSSmearingFunction', smear_lorentzian, this%smear_func)
 
    !%Variable DOSEnergyMin
    !%Type float
    !%Section Output
    !%Description
    !% Lower bound for the energy mesh of the DOS.
    !% The default is the lowest eigenvalue, minus a quarter of the total range of eigenvalues.
    !% This is ignored for the joint density of states, and the minimal energy is always set to zero.
    !%End
    call parse_variable(namespace, 'DOSEnergyMin', evalmin - eextend, this%emin, units_inp%energy)
 
    !%Variable DOSEnergyMax
    !%Type float
    !%Section Output
    !%Description
    !% Upper bound for the energy mesh of the DOS.
    !% The default is the highest eigenvalue, plus a quarter of the total range of eigenvalues.
    !%End
    call parse_variable(namespace, 'DOSEnergyMax', evalmax + eextend, this%emax, units_inp%energy)
 
    !%Variable DOSEnergyPoints
    !%Type integer
    !%Default 500
    !%Section Output
    !%Description
    !% Determines how many energy points <tt>Octopus</tt> should use for
    !% the DOS energy grid.
    !%End
    call parse_variable(namespace, 'DOSEnergyPoints', 500, this%epoints)
 
    !%Variable DOSGamma
    !%Type float
    !%Default 0.008 Ha
    !%Section Output
    !%Description
    !% Determines the width of the Lorentzian which is used for the DOS sum.
    !%End
    call parse_variable(namespace, 'DOSGamma', 0.008_real64, this%gamma)
 
    !%Variable DOSComputePDOS
    !%Type logical
    !%Default false
    !%Section Output
    !%Description
    !% Determines if projected dos are computed or not.
    !% At the moment, the PDOS is computed from the bare pseudo-atomic orbitals, directly taken from
    !% the pseudopotentials. The orbitals are not orthonormalized, in order to preserve their
    !% atomic orbitals character. As a consequence, the sum of the different PDOS does not integrate
    !% to the total DOS.
    !%
    !% The radii of the orbitals are controled by the threshold defined by <tt>AOThreshold</tt>,
    !% and the fact that they are normalized or not by <tt>AONormalize</tt>.
    !%End
    call parse_variable(namespace, 'DOSComputePDOS', .false., this%computepdos)
 
    ! spacing for energy mesh
    this%de = (this%emax - this%emin) / (this%epoints - 1)
 
    !%Variable LDOSEnergies
    !%Type block
    !%Section Output
    !%Description
    !% Specifies the energies at which the LDOS is computed.
    !%End
    if (parse_block(global_namespace, 'LDOSEnergies', blk) == 0) then
      ! There is one high symmetry k-point per line
      this%ldos_nenergies = parse_block_cols(blk, 0)
 
      safe_allocate(this%ldos_energies(1:this%ldos_nenergies))
      do ie = 1, this%ldos_nenergies
        call parse_block_float(blk, 0, ie-1, this%ldos_energies(ie))
      end do
      call parse_block_end(blk)
    else
      this%ldos_nenergies = -1
    end if
 
    pop_sub(dos_init)
  end subroutine dos_init
 
  subroutine dos_end(this)
    type(dos_t),  intent(inout) :: this
 
    push_sub(dos_end)
 
    safe_deallocate_a(this%ldos_energies)
    this%ldos_nenergies = -1
 
    pop_sub(dos_end)
  end subroutine
 
  ! ---------------------------------------------------------
  subroutine dos_write_dos(this, dir, st, box, ions, mesh, hm, namespace)
    type(dos_t),                 intent(in) :: this
    character(len=*),            intent(in) :: dir
    type(states_elec_t), target, intent(in) :: st
    class(box_t),                intent(in) :: box
    type(ions_t),        target, intent(in) :: ions
    class(mesh_t),               intent(in) :: mesh
    type(hamiltonian_elec_t),    intent(in) :: hm
    type(namespace_t),           intent(in) :: namespace
 
    integer :: ie, ik, ist, is, ns, maxdos, ib, ind
    integer, allocatable :: iunit(:)
    real(real64)   :: energy
    real(real64), allocatable :: tdos(:)
    real(real64), allocatable :: dos(:,:,:)
    character(len=64)  :: filename,format_str
    logical :: normalize
 
    integer :: ii, ll, mm, nn, work, norb, work2
    integer :: ia, iorb, idim
    real(real64)   :: threshold
    real(real64), allocatable :: ddot(:,:,:)
    complex(real64), allocatable :: zdot(:,:,:)
    real(real64), allocatable :: weight(:,:,:)
    type(orbitalset_t) :: os
    type(wfs_elec_t), pointer :: epsib
 
    type(smear_t) :: smear
    real(real64) :: e_simplex(20)
    real(real64) :: dos_simplex
 
    push_sub(dos_write_dos)
 
    ! shortcuts
    ns = 1
    if (st%d%nspin == 2) ns = 2
 
    ! set up smearing parameters
    smear%method = this%smear_func
    smear%MP_n = 1
 
    if (st%system_grp%is_root()) then
      ! space for state-dependent DOS
      safe_allocate(dos(1:this%epoints, 1:st%nst, 0:ns-1))
      safe_allocate(iunit(0:ns-1))
 
      ! compute band/spin-resolved density of states
      do ist = 1, st%nst
 
        do is = 0, ns-1
          if (ns > 1) then
            write(filename, '(a,i5.5,a,i1.1,a)') 'dos-', ist, '-', is+1,'.dat'
          else
            write(filename, '(a,i5.5,a)') 'dos-', ist, '.dat'
          end if
          iunit(is) = io_open(trim(dir)//'/'//trim(filename), namespace, action='write')
          ! write header
          write(iunit(is), '(3a)') '# energy [', trim(units_abbrev(units_out%energy)), '], band-resolved DOS'
        end do
 
        do ie = 1, this%epoints
          energy = this%emin + (ie - 1) * this%de
          dos(ie, ist, :) = m_zero
          select case (this%method)
          case (option__dosmethod__smear)
            ! sum up Lorentzians
            do ik = 1, st%nik, ns
              do is = 0, ns-1
                dos(ie, ist, is) = dos(ie, ist, is) + st%kweights(ik+is) / this%gamma * &
                  smear_delta_function(smear, (energy - st%eigenval(ist, ik+is)) / this%gamma)
              end do
            end do
          case (option__dosmethod__tetrahedra, option__dosmethod__tetrahedra_opt)
            ! sum up tetrahedra
            assert(associated(hm%kpoints%simplex))
            call profiling_in("DOS_WRITE_TETRAHEDRA")
            !$omp parallel do private(ii, is, ll, ik, e_simplex, dos_simplex) reduction(+:dos)
            do ii = 1, hm%kpoints%simplex%n_simplices
              do is = 0, ns - 1
                do ll = 1, hm%kpoints%simplex%sdim
                  ik = hm%kpoints%simplex%simplices(ii, ll)
                  ik = ns * (ik - 1) + 1
                  e_simplex(ll) = st%eigenval(ist, ik+is)
                end do
                call simplex_dos(hm%kpoints%simplex%rdim, e_simplex(1:hm%kpoints%simplex%sdim), energy, dos_simplex)
                dos(ie, ist, is) = dos(ie, ist, is) + dos_simplex / hm%kpoints%simplex%n_simplices
              end do
            end do
            !$omp end parallel do
            call profiling_out("DOS_WRITE_TETRAHEDRA")
          case default
            assert(.false.)
          end select
          do is = 0, ns-1
            write(message(1), '(2es25.16E3)') units_from_atomic(units_out%energy, energy), &
              units_from_atomic(unit_one / units_out%energy, dos(ie, ist, is))
            call messages_info(1, iunit(is))
          end do
        end do
 
        do is = 0, ns-1
          call io_close(iunit(is))
        end do
      end do
 
      safe_allocate(tdos(1))
 
      ! for spin-polarized calculations also output spin-resolved tDOS
      if (st%d%nspin == spin_polarized) then
        do is = 0, ns-1
          write(filename, '(a,i1.1,a)') 'total-dos-', is+1,'.dat'
          iunit(is) = io_open(trim(dir)//'/'//trim(filename), namespace, action='write')
          ! write header
          write(iunit(is), '(3a)') '# energy [', trim(units_abbrev(units_out%energy)), '], total DOS (spin-resolved)'
 
          do ie = 1, this%epoints
            energy = this%emin + (ie - 1) * this%de
            tdos(1) = m_zero
            do ist = 1, st%nst
              tdos(1) = tdos(1) + dos(ie, ist, is)
            end do
            write(message(1), '(2es25.16E3)') units_from_atomic(units_out%energy, energy), &
              units_from_atomic(unit_one / units_out%energy, tdos(1))
            call messages_info(1, iunit(is))
          end do
 
          call io_close(iunit(is))
        end do
      end if
 
 
      iunit(0) = io_open(trim(dir)//'/'//'total-dos.dat', namespace, action='write')
      write(iunit(0), '(3a)') '# energy [', trim(units_abbrev(units_out%energy)), '], total DOS'
 
      ! compute total density of states
      do ie = 1, this%epoints
        energy = this%emin + (ie - 1) * this%de
        tdos(1) = m_zero
        do ist = 1, st%nst
          do is = 0, ns-1
            tdos(1) = tdos(1) + dos(ie, ist, is)
          end do
        end do
        write(message(1), '(2es25.16E3)') units_from_atomic(units_out%energy, energy), &
          units_from_atomic(unit_one / units_out%energy, tdos(1))
        call messages_info(1, iunit(0))
      end do
 
      call io_close(iunit(0))
 
      safe_deallocate_a(tdos)
 
 
      ! write Fermi file
      iunit(0) = io_open(trim(dir)//'/'//'total-dos-efermi.dat', namespace, action='write')
      write(message(1), '(3a)') '# Fermi energy [', trim(units_abbrev(units_out%energy)), &
        '] in a format compatible with total-dos.dat'
 
      ! this is the maximum that tdos can reach
      maxdos = st%smear%el_per_state * st%nst
 
      write(message(2), '(2es25.16E3)') units_from_atomic(units_out%energy, st%smear%e_fermi), m_zero
      write(message(3), '(es25.16E3,i7)') units_from_atomic(units_out%energy, st%smear%e_fermi), maxdos
 
      call messages_info(3, iunit(0))
      call io_close(iunit(0))
 
    end if
 
    if (this%computepdos) then
 
      !These variables are defined in basis_set/orbitalbasis.F90
      call parse_variable(namespace, 'AOThreshold', 0.01_real64, threshold)
      call parse_variable(namespace, 'AONormalize', .true., normalize)
 
      do ia = 1, ions%natoms
        !We first count how many orbital set we have
        work = orbitalset_utils_count(ions%atom(ia)%species)
 
        !We loop over the orbital sets of the atom ia
        do norb = 1, work
          call orbitalset_init(os)
          os%spec => ions%atom(ia)%species
 
          !We count the orbitals
          work2 = 0
          do iorb = 1, os%spec%get_niwfs()
            call os%spec%get_iwf_ilm(iorb, 1, ii, ll, mm)
            call os%spec%get_iwf_n(iorb, 1, nn)
            if (ii == norb) then
              os%ll = ll
              os%nn = nn
              os%ii = ii
              os%radius = atomic_orbital_get_radius(os%spec, mesh, iorb, 1, &
                option__aotruncation__ao_full, threshold)
              work2 = work2 + 1
            end if
          end do
          os%norbs = work2
          os%ndim = 1
          os%use_submesh = .false.
          os%allocated_on_mesh = .true.
          os%spec => ions%atom(ia)%species
 
          do work = 1, os%norbs
            ! We obtain the orbital
            if (states_are_real(st)) then
              call dget_atomic_orbital(namespace, ions%space, ions%latt, ions%pos(:,ia), &
                ions%atom(ia)%species, mesh, os%sphere, os%ii, os%ll, os%jj, &
                os, work, os%radius, os%ndim, use_mesh=.not.os%use_submesh, &
                normalize = normalize)
            else
              call zget_atomic_orbital(namespace, ions%space, ions%latt, ions%pos(:,ia), &
                ions%atom(ia)%species, mesh, os%sphere, os%ii, os%ll, os%jj, &
                os, work, os%radius, os%ndim, &
                use_mesh = .not. hm%phase%is_allocated() .and. .not. os%use_submesh, &
                normalize = normalize)
            end if
          end do
 
          if (hm%phase%is_allocated()) then
            ! In case of complex wavefunction, we allocate the array for the phase correction
            safe_allocate(os%phase(1:os%sphere%np, st%d%kpt%start:st%d%kpt%end))
            os%phase(:,:) = m_zero
 
            if (.not. os%use_submesh) then
              safe_allocate(os%eorb_mesh(1:mesh%np, 1:os%norbs, 1:os%ndim, st%d%kpt%start:st%d%kpt%end))
              os%eorb_mesh(:,:,:,:) = m_zero
            else
              safe_allocate(os%eorb_submesh(1:os%sphere%np, 1:os%ndim, 1:os%norbs, st%d%kpt%start:st%d%kpt%end))
              os%eorb_submesh(:,:,:,:) = m_zero
            end if
 
            if (accel_is_enabled() .and. st%d%dim == 1) then
              os%ldorbs_eorb = max(pad_pow2(os%sphere%np), 1)
              if(.not. os%use_submesh) os%ldorbs_eorb = max(pad_pow2(os%sphere%mesh%np), 1)
 
              safe_allocate(os%buff_eorb(st%d%kpt%start:st%d%kpt%end))
              do ik= st%d%kpt%start, st%d%kpt%end
                call accel_create_buffer(os%buff_eorb(ik), accel_mem_read_only, type_cmplx, os%ldorbs_eorb*os%norbs)
              end do
            end if
 
            call orbitalset_update_phase(os, box%dim, st%d%kpt, hm%kpoints, (st%d%ispin==spin_polarized), &
              vec_pot = hm%hm_base%uniform_vector_potential, &
              vec_pot_var = hm%hm_base%vector_potential)
          else
            if (states_are_real(st)) then
              call dorbitalset_transfer_to_device(os, st%d%kpt, .not. os%use_submesh)
            else
              call zorbitalset_transfer_to_device(os, st%d%kpt, .not. os%use_submesh)
            end if
          end if
 
          if (st%system_grp%is_root()) then
            if (os%nn /= 0) then
              write(filename,'(a, i4.4, a1, a, i1.1, a1,a)') 'pdos-at', ia, '-', trim(os%spec%get_label()), &
                os%nn, l_notation(os%ll), '.dat'
            else
              write(filename,'(a,  i4.4, a1, a, a1,a)') 'pdos-at', ia, '-', trim(os%spec%get_label()), &
                l_notation(os%ll), '.dat'
            end if
 
            iunit(0) = io_open(trim(dir)//'/'//trim(filename), namespace, action='write')
            ! write header
            write(iunit(0), '(3a)') '# energy [', trim(units_abbrev(units_out%energy)), &
              '], projected DOS (total and orbital resolved)'
          end if
 
          if (states_are_real(st)) then
            safe_allocate(ddot(1:st%d%dim, 1:os%norbs, 1:st%block_size))
          else
            safe_allocate(zdot(1:st%d%dim, 1:os%norbs, 1:st%block_size))
          end if
 
          safe_allocate(weight(1:os%norbs,1:st%nik,1:st%nst))
          weight(1:os%norbs,1:st%nik,1:st%nst) = m_zero
 
          do ik = st%d%kpt%start, st%d%kpt%end
            do ib = st%group%block_start, st%group%block_end
 
              if (hm%phase%is_allocated()) then
                safe_allocate(epsib)
                call st%group%psib(ib, ik)%copy_to(epsib)
                call hm%phase%apply_to(mesh, mesh%np, .false., epsib, src = st%group%psib(ib, ik))
              else
                epsib => st%group%psib(ib, ik)
              end if
 
              if (states_are_real(st)) then
                call dorbitalset_get_coeff_batch(os, st%d%dim, epsib, ddot(:, :, :))
                do ist = 1, st%group%psib(ib, ik)%nst
                  ind = st%group%psib(ib, ik)%ist(ist)
                  do iorb = 1, os%norbs
                    do idim = 1, st%d%dim
                      weight(iorb, ik, ind) = weight(iorb, ik, ind) + st%kweights(ik) * abs(ddot(idim, iorb, ist))**2
                    end do
                  end do
                end do
              else
                call zorbitalset_get_coeff_batch(os, st%d%dim, epsib, zdot(:, :, :))
                do ist = 1, st%group%psib(ib, ik)%nst
                  ind = st%group%psib(ib, ik)%ist(ist)
                  do iorb = 1, os%norbs
                    do idim = 1, st%d%dim
                      weight(iorb, ik, ind) = weight(iorb, ik, ind) + st%kweights(ik) * abs(zdot(idim, iorb, ist))**2
                    end do
                  end do
                end do
              end if
 
              if (hm%phase%is_allocated()) then
                call epsib%end(copy=.false.)
                safe_deallocate_p(epsib)
              end if
 
            end do
          end do
 
          if (st%parallel_in_states .or. st%d%kpt%parallel) then
            call comm_allreduce(st%st_kpt_mpi_grp, weight)
          end if
 
          safe_deallocate_a(ddot)
          safe_deallocate_a(zdot)
 
          if (st%system_grp%is_root()) then
            write(format_str,'(a,i5,a)') '(', os%norbs+2, 'es25.16E3)'
            safe_allocate(tdos(1:os%norbs))
            do ie = 1, this%epoints
              energy = this%emin + (ie - 1) * this%de
              do iorb = 1, os%norbs
                tdos(iorb) = m_zero
                do ist = 1, st%nst
                  do ik = 1, st%nik
                    tdos(iorb) = tdos(iorb) + weight(iorb,ik,ist) / this%gamma * &
                      smear_delta_function(smear, (energy - st%eigenval(ist, ik))/this%gamma)
                  end do
                end do
              end do
 
              write(iunit(0), trim(format_str)) units_from_atomic(units_out%energy, energy), &
                units_from_atomic(unit_one / units_out%energy, sum(tdos)), &
                (units_from_atomic(unit_one / units_out%energy, tdos(iorb)), iorb=1,os%norbs)
            end do
            safe_deallocate_a(tdos)
            call io_close(iunit(0))
          end if
 
          call orbitalset_end(os)
          safe_deallocate_a(weight)
 
        end do
 
      end do
 
    end if
 
    safe_deallocate_a(iunit)
    safe_deallocate_a(dos)
 
    pop_sub(dos_write_dos)
  end subroutine dos_write_dos
 
  ! ---------------------------------------------------------
  subroutine dos_write_jdos(this, dir, st, ions, hm, namespace)
    type(dos_t),               intent(in) :: this
    character(len=*),         intent(in) :: dir
    type(states_elec_t),      intent(in) :: st
    type(ions_t),     target, intent(in) :: ions
    type(hamiltonian_elec_t), intent(in) :: hm
    type(namespace_t),        intent(in) :: namespace
 
    integer :: ie, ik, val, cond, is, ns
    integer, allocatable :: iunit(:)
    real(real64) :: energy
    real(real64) :: tjdos(1)
    real(real64), allocatable :: jdos(:,:)
    character(len=64)  :: filename
 
    type(smear_t) :: smear
 
    push_sub(dos_write_jdos)
 
    ! shortcuts
    ns = 1
    if (st%d%nspin == 2) ns = 2
 
    ! set up smearing parameters
    smear%method = this%smear_func
    smear%MP_n = 1
 
    if (st%system_grp%is_root()) then
      ! space for state-dependent DOS
      safe_allocate(jdos(1:this%epoints, 0:ns-1))
      safe_allocate(iunit(0:ns-1))
      jdos = m_zero
 
      ! compute band/spin-resolved density of states
      do val = 1, st%nst
        do cond = val, st%nst
          do ik = 1, st%nik, ns
            do is = 0, ns-1
              if(st%occ(val, ik+is) < m_epsilon) cycle
              if(st%occ(cond, ik+is) > m_epsilon) cycle
              do ie = 1, this%epoints
                energy = (ie - 1) * this%de
                ! sum up Lorentzians
                jdos(ie, is) = jdos(ie, is) + st%kweights(ik+is) / this%gamma * &
                  smear_delta_function(smear, (energy - (st%eigenval(cond, ik+is)-st%eigenval(val, ik+is)))/this%gamma)
              end do
            end do
          end do
        end do
      end do
 
      ! for spin-polarized calculations also output spin-resolved tDOS
      if (st%d%nspin > 1) then
        do is = 0, ns-1
          write(filename, '(a,i1.1,a)') 'total-jdos-', is+1,'.dat'
          iunit(is) = io_open(trim(dir)//'/'//trim(filename), namespace, action='write')
          ! write header
          write(iunit(is), '(3a)') '# energy [', trim(units_abbrev(units_out%energy)), '], total JDOS (spin-resolved)'
 
          do ie = 1, this%epoints
            energy = (ie - 1) * this%de
            write(message(1), '(2es25.16E3)') units_from_atomic(units_out%energy, energy), &
              units_from_atomic(unit_one / units_out%energy, jdos(ie, is))
            call messages_info(1, iunit(is))
          end do
 
          call io_close(iunit(is))
        end do
      end if
 
 
      iunit(0) = io_open(trim(dir)//'/'//'total-jdos.dat', namespace, action='write')
      write(iunit(0), '(3a)') '# energy [', trim(units_abbrev(units_out%energy)), '], total JDOS'
 
      ! compute total joint density of states
      do ie = 1, this%epoints
        energy = (ie - 1) * this%de
        tjdos(1) = m_zero
        do is = 0, ns-1
          tjdos(1) = tjdos(1) + jdos(ie, is)
        end do
        write(message(1), '(2es25.16E3)') units_from_atomic(units_out%energy, energy), &
          units_from_atomic(unit_one / units_out%energy, tjdos(1))
        call messages_info(1, iunit(0))
      end do
 
      call io_close(iunit(0))
    end if
 
    safe_deallocate_a(iunit)
    safe_deallocate_a(jdos)
 
    pop_sub(dos_write_jdos)
  end subroutine dos_write_jdos
 
 
  ! ---------------------------------------------------------
  subroutine dos_write_ldos(this, dir, st, ions, mesh, how, namespace)
    type(dos_t),               intent(in) :: this
    character(len=*),         intent(in) :: dir
    type(states_elec_t),      intent(in) :: st
    type(ions_t),     target, intent(in) :: ions
    class(mesh_t),            intent(in) :: mesh
    integer(int64),           intent(in) :: how
    type(namespace_t),        intent(in) :: namespace
 
    integer :: ie, ik, ist, is, ns, ip, ierr
    character(len=MAX_PATH_LEN)  :: fname, name
    real(real64) :: weight
    real(real64), allocatable :: ldos(:,:,:), dpsi(:,:), abs_psi2(:)
    complex(real64), allocatable :: zpsi(:,:)
    type(unit_t) :: fn_unit
 
    type(smear_t) :: smear
 
    push_sub(dos_write_ldos)
 
    if (this%ldos_nenergies < 1) then
      message(1) = "LDOSEnergies must be defined for Output=ldos"
      call messages_fatal(1, namespace=namespace)
    end if
 
    ! shortcuts
    ns = 1
    if (st%d%nspin == 2) ns = 2
 
    ! set up smearing parameters
    smear%method = this%smear_func
    smear%MP_n = 1
 
    fn_unit = units_out%length**(-ions%space%dim) / units_out%energy
 
    ! space for state-dependent DOS
    safe_allocate(ldos(1:mesh%np, 1:this%ldos_nenergies, 1:ns))
    ldos = m_zero
 
    safe_allocate(abs_psi2(1:mesh%np))
    if (states_are_real(st)) then
      safe_allocate(dpsi(1:mesh%np, 1:st%d%dim))
    else
      safe_allocate(zpsi(1:mesh%np, 1:st%d%dim))
    end if
 
    ! compute band/spin-resolved density of states
    do ik = st%d%kpt%start, st%d%kpt%end
      is = st%d%get_spin_index(ik)
      do ist = st%st_start, st%st_end
 
        ! TODO: This will be replaced by a poin-wise batch multiplication
        if (states_are_real(st)) then
          call states_elec_get_state(st, mesh, ist, ik, dpsi)
          do ip = 1, mesh%np
            abs_psi2(ip) = dpsi(ip, 1)**2
          end do
          if (st%d%dim > 1) then
            abs_psi2(ip) =  abs_psi2(ip) + dpsi(ip, 2)**2
          end if
        else
          call states_elec_get_state(st, mesh, ist, ik, zpsi)
          do ip = 1, mesh%np
            abs_psi2(ip) = real(conjg(zpsi(ip, 1)) * zpsi(ip, 1), real64)
          end do
          if (st%d%dim > 1) then
            abs_psi2(ip) =  abs_psi2(ip) + real(conjg(zpsi(ip, 2)) * zpsi(ip, 2), real64)
          end if
        end if
 
        ! sum up Lorentzians
        do ie = 1, this%ldos_nenergies
          weight = st%kweights(ik) / this%gamma * &
            smear_delta_function(smear, (this%ldos_energies(ie) - st%eigenval(ist, ik))/this%gamma)
          call lalg_axpy(mesh%np, weight, abs_psi2, ldos(:, ie, is))
        end do
      end do
    end do !ist
 
    safe_deallocate_a(dpsi)
    safe_deallocate_a(zpsi)
    safe_deallocate_a(abs_psi2)
 
    if (st%parallel_in_states .or. st%d%kpt%parallel) then
      call comm_allreduce(st%st_kpt_mpi_grp, ldos)
    end if
 
    do is = 1, ns
      do ie = 1, this%ldos_nenergies
        write(name, '(a,i5.5)') 'ldos_en-', ie
        fname = get_filename_with_spin(name, st%d%nspin, is)
 
        call dio_function_output(how, dir, fname, namespace, ions%space, mesh, &
          ldos(:, ie, is), fn_unit, ierr, pos=ions%pos, atoms=ions%atom, grp = st%dom_st_kpt_mpi_grp)
      end do
    end do
 
    safe_deallocate_a(ldos)
 
    pop_sub(dos_write_ldos)
  end subroutine dos_write_ldos
 
 
end module dos_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: