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 grid_oct_m, only: grid_t, dgrid_symmetrize_scalar_field, dgrid_symmetrize_single
  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_t, 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_pdos,            &
    dos_write_jdos,            &
    dos_write_ldos
 
  type dos_t
    private
    real(real64)   :: emin
    real(real64)   :: emax
    integer :: epoints
    real(real64)   :: de
 
    real(real64)   :: gamma
 
    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
 
    !%Variable DOSMethod
    !%Type integer
    !%Default smear
    !%Section Output
    !%Description
    !% Selects the method that is used to calculate the DOS.
    !%Option smear 0
    !% Smearing with the smearing function selected by <tt>DOSSmearingFunction</tt>.
    !%Option tetrahedra
    !% 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
    !% 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
    !% Gaussian smearing.
    !%Option lorentzian_smearing
    !% 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)
 
    ! DOSComputePDOS = yes -> Output = pdos
    call messages_obsolete_variable(namespace, 'DOSComputePDOS' , 'Output')
 
    ! 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, nvertices
    integer, allocatable :: iunit(:)
    real(real64)   :: energy
    real(real64), allocatable :: tdos(:)
    real(real64), allocatable :: dos(:,:,:)
    character(len=64)  :: filename
 
    integer :: ii, ll
 
    type(smear_t) :: smear
    type(simplex_t), pointer :: simplex
    real(real64) :: e_simplex(20)
    real(real64), allocatable :: energies(:), dos_simplex_batch(:,:), dos_thread(:,:)
 
    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))
      dos(:, :, :) = m_zero
      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
 
        select case (this%method)
        case (option__dosmethod__smear)
          do ie = 1, this%epoints
            energy = this%emin + (ie - 1) * this%de
            ! 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
          end do
 
        case (option__dosmethod__tetrahedra, option__dosmethod__tetrahedra_opt)
          assert(associated(hm%kpoints%reduced%simplex))
          simplex => hm%kpoints%reduced%simplex
          nvertices = simplex%rdim + 1
          assert(nvertices <= 4)
          call profiling_in("DOS_WRITE_TETRAHEDRA")
          safe_allocate(energies(1:this%epoints))
          do ie = 1, this%epoints
            energies(ie) = this%emin + (ie - 1) * this%de
          end do
          !$omp parallel private(ii, is, ll, ik, ie, e_simplex, dos_simplex_batch, dos_thread)
          safe_allocate(dos_simplex_batch(1:nvertices, 1:this%epoints))
          safe_allocate(dos_thread(1:this%epoints, 0:ns-1))
          dos_thread = m_zero
          !$omp do schedule(static)
          do ii = 1, simplex%n_simplices
            do is = 0, ns - 1
              do ll = 1, simplex%sdim
                ik = simplex%simplices(ii, ll)
                ik = ns * (ik - 1) + 1
                e_simplex(ll) = st%eigenval(ist, ik+is)
              end do
 
              call simplex_dos(simplex%rdim, e_simplex(1:simplex%sdim), energies, dos_simplex_batch)
 
              do ie = 1, this%epoints
                dos_thread(ie, is) = dos_thread(ie, is) + sum(dos_simplex_batch(1:nvertices, ie)) / simplex%n_points
              end do
            end do
          end do
          !$omp end do
          !$omp critical
          dos(:, ist, :) = dos(:, ist, :) + dos_thread
          !$omp end critical
          safe_deallocate_a(dos_thread)
          safe_deallocate_a(dos_simplex_batch)
          !$omp end parallel
          safe_deallocate_a(energies)
          call profiling_out("DOS_WRITE_TETRAHEDRA")
 
        case default
          assert(.false.)
        end select
 
        do ie = 1, this%epoints
          energy = this%emin + (ie - 1) * this%de
          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
 
 
    safe_deallocate_a(iunit)
    safe_deallocate_a(dos)
 
    pop_sub(dos_write_dos)
  end subroutine dos_write_dos
 
  ! ---------------------------------------------------------
  subroutine dos_write_pdos(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, ib, ind, nvertices
    integer :: iunit
    real(real64)   :: energy, threshold
    character(len=64)  :: filename, format_str
    logical :: normalize
 
    integer :: ii, ll, mm, nn, work, norb, work2
    integer :: ia, iorb, idim
    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
    type(simplex_t), pointer :: simplex
    real(real64) :: e_simplex(20), a_simplex(4)
    real(real64), allocatable :: energies(:), dos_simplex_batch(:,:), pdos(:,:), pdos_thread(:,:)
 
    push_sub(dos_write_pdos)
 
    ! shortcuts
    ns = 1
    if (st%d%nspin == 2) ns = 2
 
    ! set up smearing parameters
    smear%method = this%smear_func
    smear%MP_n = 1
 
    !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 = io_open(trim(dir)//'/'//trim(filename), namespace, action='write')
          ! write header
          write(iunit, '(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(pdos(1:os%norbs, 1:this%epoints))
          pdos = m_zero
          select case (this%method)
          case (option__dosmethod__smear)
            do iorb = 1, os%norbs
              do ist = 1, st%nst
                do ik = 1, st%nik
                  do ie = 1, this%epoints
                    energy = this%emin + (ie - 1) * this%de
                    pdos(iorb, ie) = pdos(iorb, ie) + weight(iorb,ik,ist) / this%gamma * &
                      smear_delta_function(smear, (energy - st%eigenval(ist, ik))/this%gamma)
                  end do
                end do
              end do
            end do
          case (option__dosmethod__tetrahedra, option__dosmethod__tetrahedra_opt)
            assert(associated(hm%kpoints%reduced%simplex))
            simplex => hm%kpoints%reduced%simplex
            nvertices = simplex%rdim + 1
            assert(nvertices <= 4)
            safe_allocate(energies(1:this%epoints))
            do ie = 1, this%epoints
              energies(ie) = this%emin + (ie - 1) * this%de
            end do
            call profiling_in("DOS_WRITE_PDOS_TETRAHEDRA")
            !$omp parallel private(ist, ii, is, ll, ik, iorb, ie, e_simplex, a_simplex, dos_simplex_batch, pdos_thread)
            safe_allocate(dos_simplex_batch(1:nvertices, 1:this%epoints))
            safe_allocate(pdos_thread(1:os%norbs, 1:this%epoints))
            pdos_thread = m_zero
            !$omp do collapse(2) schedule(static)
            do ist = 1, st%nst
              do ii = 1, simplex%n_simplices
                do is = 0, ns - 1
                  do ll = 1, simplex%sdim
                    ik = simplex%simplices(ii, ll)
                    ik = ns * (ik - 1) + 1
                    e_simplex(ll) = st%eigenval(ist, ik+is)
                  end do
                  call simplex_dos(simplex%rdim, e_simplex(1:simplex%sdim), energies, dos_simplex_batch)
                  do iorb = 1, os%norbs
                    do ll = 1, nvertices
                      ik = simplex%simplices(ii, ll)
                      ik = ns * (ik - 1) + 1
                      if (st%kweights(ik+is) > m_epsilon) then
                        a_simplex(ll) = weight(iorb, ik+is, ist) / st%kweights(ik+is)
                      else
                        a_simplex(ll) = m_zero
                      end if
                    end do
                    do ie = 1, this%epoints
                      pdos_thread(iorb, ie) = pdos_thread(iorb, ie) + &
                        sum(a_simplex(1:nvertices) * dos_simplex_batch(1:nvertices, ie)) / simplex%n_points
                    end do
                  end do
                end do
              end do
            end do
            !$omp end do
            !$omp critical
            pdos(:,:) = pdos(:,:) + pdos_thread(:,:)
            !$omp end critical
            safe_deallocate_a(pdos_thread)
            safe_deallocate_a(dos_simplex_batch)
            !$omp end parallel
            call profiling_out("DOS_WRITE_PDOS_TETRAHEDRA")
            safe_deallocate_a(energies)
          case default
            assert(.false.)
          end select
          do ie = 1, this%epoints
            energy = this%emin + (ie - 1) * this%de
            write(iunit, trim(format_str)) units_from_atomic(units_out%energy, energy), &
              units_from_atomic(unit_one / units_out%energy, sum(pdos(:,ie))), &
              (units_from_atomic(unit_one / units_out%energy, pdos(iorb,ie)), iorb=1,os%norbs)
          end do
          safe_deallocate_a(pdos)
          call io_close(iunit)
        end if
 
        call orbitalset_end(os)
        safe_deallocate_a(weight)
 
      end do
 
    end do
 
    pop_sub(dos_write_pdos)
  end subroutine dos_write_pdos
 
  ! ---------------------------------------------------------
  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, ll, ii, nvertices
    integer, allocatable :: iunit(:)
    real(real64) :: energy
    real(real64) :: tjdos(1), e_simplex(20), occ_simplex(4)
    real(real64), allocatable :: jdos(:,:), energies(:), dos_simplex(:,:)
    character(len=64)  :: filename
 
    type(smear_t) :: smear
    type(simplex_t), pointer :: simplex
 
    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
 
      select case (this%method)
      case (option__dosmethod__smear)
        ! 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
      case (option__dosmethod__tetrahedra, option__dosmethod__tetrahedra_opt)
        assert(associated(hm%kpoints%reduced%simplex))
        simplex => hm%kpoints%reduced%simplex
        nvertices = simplex%rdim + 1
        assert(nvertices <= 4)
        call profiling_in("DOS_WRITE_JDOS_TETRAHEDRA")
        safe_allocate(energies(1:this%epoints))
        do ie = 1, this%epoints
          energies(ie) = (ie - 1) * this%de
        end do
        !$omp parallel reduction(+:jdos) &
        !$omp          private(val, cond, ie, ii, is, ll, ik, e_simplex, occ_simplex, dos_simplex)
        safe_allocate(dos_simplex(1:nvertices, 1:this%epoints))
        !$omp do collapse(3) schedule(dynamic)
        do val = 1, st%nst
          do cond = 1, st%nst
            do ii = 1, simplex%n_simplices
              if (cond < val) cycle
              do is = 0, ns - 1
                do ll = 1, simplex%sdim
                  ik = simplex%simplices(ii, ll)
                  ik = ns * (ik - 1) + 1
                  e_simplex(ll) = st%eigenval(cond, ik+is) - st%eigenval(val, ik+is)
                end do
 
                do ll = 1, nvertices
                  ik = simplex%simplices(ii, ll)
                  ik = ns * (ik - 1) + 1
                  if (st%occ(val, ik+is) > m_epsilon .and. st%occ(cond, ik+is) < m_epsilon) then
                    occ_simplex(ll) = m_one
                  else
                    occ_simplex(ll) = m_zero
                  end if
                end do
 
                ! use bitwise comparison for exact zero because direct comparison is forbidden by -Wcompare-reals
                if (all(transfer(occ_simplex(1:nvertices), [0_int64]) == 0_int64)) cycle
 
                call simplex_dos(simplex%rdim, e_simplex(1:simplex%sdim), energies, dos_simplex)
                do ie = 1, this%epoints
                  jdos(ie, is) = jdos(ie, is) + sum(occ_simplex(1:nvertices) * dos_simplex(1:nvertices, ie)) / &
                    simplex%n_points
                end do
              end do
            end do
          end do
        end do
        !$omp end do
        safe_deallocate_a(dos_simplex)
        !$omp end parallel
        safe_deallocate_a(energies)
        call profiling_out("DOS_WRITE_JDOS_TETRAHEDRA")
      case default
        assert(.false.)
      end select
 
      ! 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, gr, hm, 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
    type(grid_t),             intent(in) :: gr
    type(hamiltonian_elec_t), intent(in) :: hm
    integer(int64),           intent(in) :: how
    type(namespace_t),        intent(in) :: namespace
 
    integer :: ie, ik, ist, is, ns, ip, ifull, ierr, ii, ll, nvertices, ikpoint
    integer :: iop, iiop, nops, nops_max
    character(len=MAX_PATH_LEN)  :: fname, name
    logical :: has_local_corner
    real(real64) :: weight, e_simplex(20)
    real(real64), allocatable :: ldos(:,:,:), dpsi(:,:), abs_psi2(:), abs_psi2_symm(:), ldos_weights(:,:,:)
    real(real64), allocatable :: dos_simplex(:,:)
    complex(real64), allocatable :: zpsi(:,:)
    type(unit_t) :: fn_unit
 
    type(smear_t) :: smear
    type(simplex_t), pointer :: simplex
 
    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:gr%np, 1:this%ldos_nenergies, 1:ns))
    ldos = m_zero
 
    safe_allocate(abs_psi2(1:gr%np))
    if (states_are_real(st)) then
      safe_allocate(dpsi(1:gr%np, 1:st%d%dim))
    else
      safe_allocate(zpsi(1:gr%np, 1:st%d%dim))
    end if
 
    select case (this%method)
    case (option__dosmethod__smear)
      nops_max = 1
      safe_allocate(ldos_weights(1:this%ldos_nenergies, st%d%kpt%start:st%d%kpt%end, 1:nops_max))
    case (option__dosmethod__tetrahedra, option__dosmethod__tetrahedra_opt)
      assert(associated(hm%kpoints%full%simplex))
      simplex => hm%kpoints%full%simplex
      nvertices = simplex%rdim + 1
      assert(nvertices <= 4)
      nops_max = 1
      do ifull = 1, hm%kpoints%full%npoints
        nops_max = max(nops_max, hm%kpoints%get_full_symmetry_op_index(ifull))
      end do
      safe_allocate(ldos_weights(1:this%ldos_nenergies, st%d%kpt%start:st%d%kpt%end, 1:nops_max))
      safe_allocate(abs_psi2_symm(1:gr%np))
    case default
      assert(.false.)
    end select
 
    do ist = st%st_start, st%st_end
      ldos_weights(:, :, :) = m_zero
 
      select case (this%method)
      case (option__dosmethod__smear)
        do ik = st%d%kpt%start, st%d%kpt%end
          do ie = 1, this%ldos_nenergies
            ldos_weights(ie, ik, 1) = st%kweights(ik) / this%gamma * &
              smear_delta_function(smear, (this%ldos_energies(ie) - st%eigenval(ist, ik))/this%gamma)
          end do
        end do
      case (option__dosmethod__tetrahedra, option__dosmethod__tetrahedra_opt)
        call profiling_in('DOS_WRITE_LDOS_TETRAHEDRA')
        !$omp parallel reduction(+:ldos_weights) &
        !$omp          private(ii, is, ll, ik, ie, e_simplex, dos_simplex, ifull, iiop, has_local_corner)
        safe_allocate(dos_simplex(1:nvertices, 1:this%ldos_nenergies))
        !$omp do schedule(dynamic)
        do ii = 1, simplex%n_simplices
          do is = 0, ns - 1
            has_local_corner = .false.
            do ll = 1, simplex%sdim
              ik = hm%kpoints%get_equiv(simplex%simplices(ii, ll))
              ik = ns * (ik - 1) + 1 + is
              e_simplex(ll) = st%eigenval(ist, ik)
              if (ik >= st%d%kpt%start .and. ik <= st%d%kpt%end) has_local_corner = .true.
            end do
            if (.not. has_local_corner) cycle
 
            call simplex_dos(simplex%rdim, e_simplex(1:simplex%sdim), this%ldos_energies, &
              dos_simplex(1:nvertices, 1:this%ldos_nenergies))
            do ie = 1, this%ldos_nenergies
              do ll = 1, nvertices
                ifull = simplex%simplices(ii, ll)
                iiop = hm%kpoints%get_full_symmetry_op_index(ifull)
                ik = hm%kpoints%get_equiv(ifull)
                ik = ns * (ik - 1) + 1 + is
                if (ik >= st%d%kpt%start .and. ik <= st%d%kpt%end) then
                  ldos_weights(ie, ik, iiop) = ldos_weights(ie, ik, iiop) + dos_simplex(ll, ie) / simplex%n_points
                end if
              end do
            end do
          end do
        end do
        !$omp end do
        safe_deallocate_a(dos_simplex)
        !$omp end parallel
        call profiling_out('DOS_WRITE_LDOS_TETRAHEDRA')
      case default
        assert(.false.)
      end select
 
      do ik = st%d%kpt%start, st%d%kpt%end
        is = st%d%get_spin_index(ik)
 
        if (states_are_real(st)) then
          call states_elec_get_state(st, gr, ist, ik, dpsi)
          do ip = 1, gr%np
            abs_psi2(ip) = dpsi(ip, 1)**2
          end do
          if (st%d%dim > 1) then
            do ip = 1, gr%np
              abs_psi2(ip) =  abs_psi2(ip) + dpsi(ip, 2)**2
            end do
          end if
        else
          call states_elec_get_state(st, gr, ist, ik, zpsi)
          do ip = 1, gr%np
            abs_psi2(ip) = real(conjg(zpsi(ip, 1)) * zpsi(ip, 1), real64)
          end do
          if (st%d%dim > 1) then
            do ip = 1, gr%np
              abs_psi2(ip) =  abs_psi2(ip) + real(conjg(zpsi(ip, 2)) * zpsi(ip, 2), real64)
            end do
          end if
        end if
 
        select case (this%method)
        case (option__dosmethod__smear)
          do ie = 1, this%ldos_nenergies
            weight = ldos_weights(ie, ik, 1)
            call lalg_axpy(gr%np, weight, abs_psi2, ldos(:, ie, is))
          end do
        case (option__dosmethod__tetrahedra, option__dosmethod__tetrahedra_opt)
          ikpoint = st%d%get_kpoint_index(ik)
          nops = kpoints_get_num_symmetry_ops(hm%kpoints, ikpoint)
          do iiop = 1, nops
            iop = abs(kpoints_get_symmetry_ops(hm%kpoints, ikpoint, iiop))
            call dgrid_symmetrize_single(gr, iop, abs_psi2, abs_psi2_symm, suppress_warning = .true.)
            do ie = 1, this%ldos_nenergies
              weight = ldos_weights(ie, ik, iiop)
              call lalg_axpy(gr%np, weight, abs_psi2_symm, ldos(:, ie, is))
            end do
          end do
        case default
          assert(.false.)
        end select
      end do
    end do
 
    safe_deallocate_a(ldos_weights)
    safe_deallocate_a(dpsi)
    safe_deallocate_a(zpsi)
    safe_deallocate_a(abs_psi2)
    safe_deallocate_a(abs_psi2_symm)
 
    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)
 
        if (hm%kpoints%use_symmetries .and. &
          .not. (this%method == option__dosmethod__tetrahedra .or. this%method == option__dosmethod__tetrahedra_opt)) then
          call dgrid_symmetrize_scalar_field(gr, ldos(:, ie, is), suppress_warning = .true.)
        end if
        call dio_function_output(how, dir, fname, namespace, ions%space, gr, &
          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: