states_elec_io.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!!
!! 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 states_elec_io_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
  use io_oct_m
  use ions_oct_m
  use, intrinsic :: iso_fortran_env
  use kpoints_oct_m
  use math_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use orbitalset_oct_m
  use orbitalset_utils_oct_m
  use parser_oct_m
  use phase_oct_m
  use profiling_oct_m
  use smear_oct_m
  use sort_oct_m
  use space_oct_m
  use species_oct_m
  use states_abst_oct_m
  use states_elec_oct_m
  use states_elec_dim_oct_m
  use submesh_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use utils_oct_m
  use types_oct_m
  use wfs_elec_oct_m
 
  implicit none
 
  private
 
  public ::                           &
    states_elec_write_eigenvalues,         &
    states_elec_write_tpa,                 &
    states_elec_write_bandstructure,       &
    states_elec_write_gaps
 
contains
 
  ! ---------------------------------------------------------
  !
  subroutine states_elec_write_eigenvalues(nst, st, space, kpoints, error, st_start, compact, iunit, namespace)
    integer,                     intent(in) :: nst
    type(states_elec_t),         intent(in) :: st
    class(space_t),              intent(in) :: space
    type(kpoints_t),             intent(in) :: kpoints
    real(real64),      optional, intent(in) :: error(:,:)
    integer,           optional, intent(in) :: st_start
    logical,           optional, intent(in) :: compact
    integer,           optional, intent(in) :: iunit
    type(namespace_t), optional, intent(in) :: namespace
 
    integer :: ik, ikk, ist, ns, is, idir, st_start_, iflat, iqn, homo_index, not_printed
    logical :: print_eigenval
    real(real64) :: kpoint(space%dim), max_error
    character(len=120) :: tmp_str(max(space%dim, 3)), cspin
 
    real(real64), allocatable :: flat_eigenval(:)
    integer, allocatable :: flat_indices(:, :)
    integer, parameter :: print_range = 8
 
    push_sub(states_elec_write_eigenvalues)
 
    st_start_ = 1
    if (present(st_start)) st_start_ = st_start
    assert(nst <= st%nst)
 
 
    if (.not. mpi_grp_is_root(mpi_world)) then
      pop_sub(states_elec_write_eigenvalues)
      return
    end if
 
    if (.not. optional_default(compact, .false.)) then
 
      ns = 1
      if (st%d%nspin == 2) ns = 2
 
      message(1) = 'Eigenvalues [' // trim(units_abbrev(units_out%energy)) // ']'
      call messages_info(1, iunit=iunit, namespace=namespace)
 
      if (st%d%ispin == spinors) then
        write(message(1), '(a4,1x,a5,1x,a12,1x,a12,2x,a4,4x,a4,4x,a4)')   &
          '#st',' Spin',' Eigenvalue', 'Occupation ', '<Sx>', '<Sy>', '<Sz>'
      else
        write(message(1), '(a4,1x,a5,1x,a12,4x,a12)')       &
          '#st',' Spin',' Eigenvalue', 'Occupation'
      end if
      if (present(error)) then
        write(message(1),'(a,a10)') trim(message(1)), ' Error'
      end if
      call messages_info(1, iunit=iunit, namespace=namespace)
 
      do ik = 1, st%nik, ns
        if (space%is_periodic()) then
          ikk = st%d%get_kpoint_index(ik)
          kpoint(1:space%dim) = kpoints%get_point(ikk, absolute_coordinates = .false.)
          write(message(1), '(a,i8,a)') '#k =', ikk, ', k = ('
          do idir = 1, space%dim
            write(tmp_str(1), '(f10.6)') kpoint(idir)
            message(1) = trim(message(1))//trim(tmp_str(1))
            if (idir < space%dim) message(1) = trim(message(1))//','
          end do
          message(1) = trim(message(1))//')'
          call messages_info(1, iunit=iunit, namespace=namespace)
        end if
 
        do ist = st_start_, nst
          do is = 0, ns-1
            if (is == 0) cspin = 'up'
            if (is == 1) cspin = 'dn'
            if (st%d%ispin == unpolarized .or. st%d%ispin == spinors) cspin = '--'
 
            write(tmp_str(1), '(i4,3x,a2)') ist, trim(cspin)
            if (st%d%ispin == spinors) then
              write(tmp_str(2), '(1x,f12.6,5x,f5.2,3x,3f8.4)') &
                units_from_atomic(units_out%energy, st%eigenval(ist, ik)), st%occ(ist, ik), st%spin(1:3, ist, ik)
              if (present(error)) write(tmp_str(3), '(a3,es8.1,a1)')'  (', error(ist, ik), ')'
            else
              write(tmp_str(2), '(1x,f12.6,3x,f12.6)') &
                units_from_atomic(units_out%energy, st%eigenval(ist, ik+is)), st%occ(ist, ik+is)
              if (present(error)) write(tmp_str(3), '(a7,es8.1,a1)')'      (', error(ist, ik+is), ')'
            end if
            if (present(error)) then
              message(1) = trim(tmp_str(1))//trim(tmp_str(2))//trim(tmp_str(3))
            else
              message(1) = trim(tmp_str(1))//trim(tmp_str(2))
            end if
            call messages_info(1, iunit=iunit, namespace=namespace)
          end do
        end do
      end do
 
    else
 
      call messages_info(1, iunit=iunit, namespace=namespace)
 
      safe_allocate(flat_eigenval(1:st%nik*nst))
      safe_allocate(flat_indices(1:2, 1:st%nik*nst))
 
      iflat = 1
      do iqn = 1, st%nik
        do ist = 1, nst
 
          flat_eigenval(iflat) = st%eigenval(ist, iqn)
          flat_indices(1:2, iflat) = (/iqn, ist/)
 
          iflat = iflat + 1
        end do
      end do
 
      call sort(flat_eigenval, flat_indices(:, :))
 
      homo_index = st%nik*nst
      do iflat = 1, st%nik*nst
        iqn = flat_indices(1, iflat)
        ist = flat_indices(2, iflat)
        if (abs(st%occ(ist, iqn)) < 0.1_real64) then
          homo_index = iflat - 1
          exit
        end if
      end do
 
      tmp_str(1) = '#  State'
 
      if (space%is_periodic()) tmp_str(1) = trim(tmp_str(1))//'  KPoint'
 
      if (st%d%ispin == spin_polarized) tmp_str(1) = trim(tmp_str(1))//'  Spin'
 
      tmp_str(1) = trim(tmp_str(1))//'  Eigenvalue ['// trim(units_abbrev(units_out%energy)) // ']'
 
      tmp_str(1) = trim(tmp_str(1))//'  Occupation'
 
      if (st%d%ispin == spinors) then
        tmp_str(1) = trim(tmp_str(1))//'      <Sx>     <Sy>     <Sz>'
      end if
 
      if (present(error)) tmp_str(1) = trim(tmp_str(1))//'    Error'
 
      call messages_write(tmp_str(1))
      call messages_info(iunit=iunit, namespace=namespace)
 
      not_printed = 0
      max_error = 0.0_real64
      do iflat = 1, st%nik*nst
        iqn = flat_indices(1, iflat)
        ist = flat_indices(2, iflat)
        ik = st%d%get_kpoint_index(iqn)
        is = st%d%get_spin_index(iqn)
 
        print_eigenval = iflat <= print_range
        print_eigenval = print_eigenval .or. st%nik*nst - iflat < print_range
        print_eigenval = print_eigenval .or. abs(iflat - homo_index) <= print_range
 
        if (print_eigenval) then
 
          if (not_printed > 0) then
            call messages_write('')
            call messages_new_line()
            call messages_write('  [output of ')
            call messages_write(not_printed)
            call messages_write(' eigenvalues skipped')
            if (present(error)) then
              call messages_write(': maximum error =')
              call messages_write(max_error, fmt = '(es7.1)', align_left = .true.)
            end if
            call messages_write(']')
            call messages_new_line()
            call messages_write('')
            call messages_info(iunit=iunit, namespace=namespace)
 
            not_printed = 0
            max_error = 0.0_real64
 
          end if
 
          write(tmp_str(1), '(i7)') ist
 
          if (space%is_periodic()) then
            write(tmp_str(1), '(2a,i7)') trim(tmp_str(1)), ' ', ik
          end if
 
          if (st%d%ispin == spin_polarized) then
            if (is == 1) cspin = '   up'
            if (is == 2) cspin = '   dn'
            write(tmp_str(1), '(2a,a5)') trim(tmp_str(1)), ' ', cspin
          end if
 
 
          if (len(units_abbrev(units_out%energy)) == 1) then
            write(tmp_str(1), '(2a,f14.6)') trim(tmp_str(1)), ' ', units_from_atomic(units_out%energy, st%eigenval(ist, iqn))
          else
            write(tmp_str(1), '(2a,f15.6)') trim(tmp_str(1)), ' ', units_from_atomic(units_out%energy, st%eigenval(ist, iqn))
          end if
 
          write(tmp_str(1), '(2a,f11.6)') trim(tmp_str(1)), ' ', st%occ(ist, iqn)
 
          if (st%d%ispin == spinors) then
            write(tmp_str(1), '(2a,3f9.4)') trim(tmp_str(1)), ' ', st%spin(1:3, ist, iqn)
          end if
 
          if (present(error)) then
            write(tmp_str(1), '(2a,es8.1,a)') trim(tmp_str(1)), '   (', error(ist, iqn), ')'
          end if
 
          call messages_write(tmp_str(1))
          call messages_info(iunit=iunit, namespace=namespace)
 
        else
 
          not_printed = not_printed + 1
 
          if (present(error)) then
            max_error = max(max_error, error(ist, iqn))
          end if
 
        end if
 
      end do
 
      if (nst*st%nik > 1) call print_dos()
 
      safe_deallocate_a(flat_indices)
      safe_deallocate_a(flat_eigenval)
 
    end if
 
    call smear_write_info(st%smear, namespace, iunit)
 
    pop_sub(states_elec_write_eigenvalues)
 
  contains
 
    subroutine print_dos()
 
      integer, parameter :: ndiv = 70, height = 10
      integer :: histogram(1:ndiv), iev, ien, iline, ife
      real(real64) :: dhistogram(1:ndiv)
      character(len=ndiv) :: line
      real(real64) :: emin, emax, de, maxhist
 
      push_sub(states_elec_write_eigenvalues.print_dos)
 
      emin = flat_eigenval(1)
      emax = flat_eigenval(st%nik*nst)
      de = (emax - emin)/(ndiv - m_one)
 
      if (de < m_epsilon) then
        pop_sub(states_elec_write_eigenvalues.print_dos)
        return
      end if
 
      ife = nint((st%smear%e_fermi - emin)/de) + 1
 
      histogram = 0
      dhistogram = m_zero
      do iev = 1, st%nik*nst
        ien = nint((flat_eigenval(iev) - emin)/de) + 1
        assert(ien >= 1)
        assert(ien <= ndiv)
        dhistogram(ien) = dhistogram(ien) + st%kweights(flat_indices(1, iev))*kpoints%full%npoints
      end do
 
      !normalize
      maxhist = maxval(dhistogram)
      do ien = 1, ndiv
        histogram(ien) = nint(dhistogram(ien)*height/maxhist)
      end do
 
      call messages_new_line()
      call messages_write('Density of states:')
      call messages_new_line()
      call messages_info(iunit=iunit, namespace=namespace)
 
      !print histogram
      do iline = height, 1, -1
        do ien = 1, ndiv
          if (histogram(ien) >= iline) then
            call messages_write('%')
          else
            call messages_write('-')
          end if
        end do
        call messages_info(iunit=iunit, namespace=namespace)
      end do
 
      line(1:ndiv) = ' '
      line(ife:ife) = '^'
      call messages_write(line)
      call messages_new_line()
      call messages_info(iunit=iunit, namespace=namespace)
 
      pop_sub(states_elec_write_eigenvalues.print_dos)
    end subroutine print_dos
 
  end subroutine states_elec_write_eigenvalues
 
  ! ---------------------------------------------------------
  !
  subroutine states_elec_write_gaps(iunit, st, space)
    integer,             intent(in) :: iunit
    type(states_elec_t), intent(in) :: st
    class(space_t),      intent(in) :: space
 
    integer :: ik, ist
 
    real(real64) :: homo, lumo, egdir, egindir
    integer :: homok, lumok, egdirk
 
    push_sub(states_elec_write_gaps)
 
    if (.not. mpi_grp_is_root(mpi_world) .or. .not. space%is_periodic() &
      .or. (st%smear%method  /=  smear_semiconductor .and. st%smear%method  /=  smear_fixed_occ)) then
      pop_sub(states_elec_write_gaps)
      return
    end if
 
    !We do not want to compute the bandgap in case there is no extra states
    if (ceiling(st%qtot/st%smear%el_per_state) == st%nst) then
      pop_sub(states_elec_write_gaps)
      return
    end if
 
    homo = -m_huge
    homok = -1
    lumo =  m_huge
    lumok = -1
    egdir = m_huge
    egdirk = -1
 
    !TODO: This definition depends on the type of smearing
    do ik = 1, st%nik
      if (abs(st%kweights(ik)) < m_epsilon) cycle
      do ist = 1,st%nst-1
        if (st%occ(ist, ik) > m_epsilon .and. st%eigenval(ist, ik) > homo) then
          homo = st%eigenval(ist, ik)
          homok = ik
        end if
 
        if (st%occ(ist+1, ik) <= m_epsilon .and. st%eigenval(ist+1, ik) < lumo) then
          lumo = st%eigenval(ist+1, ik)
          lumok = ik
        end if
 
        if (st%occ(ist, ik) > m_epsilon .and. st%occ(ist+1, ik) <= m_epsilon &
          .and. (st%eigenval(ist+1, ik) - st%eigenval(ist, ik)) < egdir) then
          egdir = (st%eigenval(ist+1, ik) - st%eigenval(ist, ik))
          egdirk = ik
        end if
      end do
    end do
 
    egindir = lumo-homo
 
    if (lumok == -1 .or. egdir <= m_epsilon) then
      write(message(1),'(a)') 'The system seems to have no gap.'
      call messages_info(1, iunit)
    else
      write(message(1),'(a,i5,a,f7.4,a)') 'Direct gap at ik=', egdirk, ' of ', &
        units_from_atomic(units_out%energy, egdir),  ' ' // trim(units_abbrev(units_out%energy))
      write(message(2),'(a,i5,a,i5,a,f7.4,a)') 'Indirect gap between ik=', homok, ' and ik=', lumok, &
        ' of ', units_from_atomic(units_out%energy, egindir),  ' ' // trim(units_abbrev(units_out%energy))
      call messages_info(2, iunit)
    end if
 
 
    pop_sub(states_elec_write_gaps)
  end subroutine states_elec_write_gaps
 
 
  ! ---------------------------------------------------------
  !
  subroutine states_elec_write_tpa(dir, namespace, space, mesh, st)
    character(len=*),    intent(in) :: dir
    type(namespace_t),   intent(in) :: namespace
    class(space_t),      intent(in) :: space
    class(mesh_t),       intent(in) :: mesh
    type(states_elec_t), intent(in) :: st
 
    type(block_t) :: blk
    integer       :: ncols, icoord, ist, ik, tpa_initialst, tpa_initialk
    integer       :: iunit
 
    real(real64), allocatable  :: ff(:)
    complex(real64), allocatable  :: cff(:)
    real(real64), allocatable  :: osc(:)
    real(real64)        :: transition_energy, osc_strength, dsf
 
    logical             :: use_qvector
    real(real64), allocatable  :: qvector(:), psi_initial(:, :), psi_ist(:, :)
 
    push_sub(states_elec_write_tpa)
 
    use_qvector = .false.
 
    ! find the orbital with half-occupation
    tpa_initialst = -1
    do ist = 1, st%nst
      do ik = 1, st%nik
        if (abs(st%occ(ist,ik)-m_half)  <  m_min_occ) then
          tpa_initialst = ist
          tpa_initialk  = ik
        end if
      end do
    end do
 
    ! make sure that half-occupancy was found
    if (tpa_initialst == -1) then
      if (mpi_grp_is_root(mpi_world)) then
 
        call messages_write('No orbital with half-occupancy found. TPA output is not written.')
        call messages_warning(namespace=namespace)
 
      end if
      pop_sub(states_elec_write_tpa)
      return
    end if
 
    !%Variable MomentumTransfer
    !%Type block
    !%Section Output
    !%Description
    !% Momentum-transfer vector <math>\vec{q}</math> to be used when calculating matrix elements
    !% <math>\left< f \left| e^{i \vec{q} \cdot \vec{r}} \right| i \right></math>.
    !% This enables the calculation of the dynamical structure factor,
    !% which is closely related to generalized oscillator strengths.
    !% If the vector is not given, but TPA output is requested (<tt>Output = TPA</tt>),
    !% only the oscillator strengths are written in the output file.
    !% For example, to use <math>\vec{q}</math> = (0.1, 0.2, 0.3), set
    !%
    !% <tt>%MomentumTransfer
    !% <br>&nbsp;&nbsp; 0.1 | 0.2 | 0.3
    !% <br>%</tt>
    !%End
    if (parse_block(namespace, 'MomentumTransfer', blk) == 0) then
 
      ! check if input makes sense
      ncols = parse_block_cols(blk, 0)
 
      if (ncols /= space%dim) then ! wrong size
 
        if (mpi_grp_is_root(mpi_world)) then
          call messages_write('Inconsistent size of momentum-transfer vector. It will not be used in the TPA calculation.')
          call messages_warning(namespace=namespace)
        end if
 
      else ! correct size
 
        use_qvector = .true.
        safe_allocate(qvector(1:space%dim))
 
        do icoord = 1, space%dim    !for x,y,z
          call parse_block_float(blk, 0, icoord-1, qvector(icoord))
          qvector(icoord) = units_to_atomic(unit_one / units_inp%length, qvector(icoord))
        end do
 
      end if
 
    end if
 
    ! calculate the matrix elements
 
    safe_allocate(ff(1:mesh%np))
    if (use_qvector) then
      safe_allocate(cff(1:mesh%np))
    end if
    safe_allocate(osc(1:space%dim))
 
    ! root writes output to file
 
    if (mpi_grp_is_root(mpi_world)) then
 
      iunit = io_open(trim(dir)//'/'//trim('tpa_xas'), namespace, action='write')
 
      ! header
      if (use_qvector) then
        write (message(1),'(a1,a30,3(es14.5,1x),a1)') '#', ' momentum-transfer vector : (', &
          (units_from_atomic(unit_one / units_out%length, qvector(icoord)), icoord=1, space%dim),')'
        select case (space%dim)
        case (1)
          write(message(2), '(a1,4(a15,1x))') '#', 'E' , '<x>', '<f>', 'S(q,omega)'
        case (2)
          write(message(2), '(a1,5(a15,1x))') '#', 'E' , '<x>', '<y>', '<f>', 'S(q,omega)'
        case (3)
          write(message(2), '(a1,6(a15,1x))') '#', 'E' , '<x>', '<y>', '<z>', '<f>', 'S(q,omega)'
        end select
        call messages_info(2, iunit, namespace=namespace)
      else
        select case (space%dim)
        case (1)
          write(message(1), '(a1,3(a15,1x))') '#', 'E' , '<x>', '<f>'
        case (2)
          write(message(1), '(a1,4(a15,1x))') '#', 'E' , '<x>', '<y>', '<f>'
        case (3)
          write(message(1), '(a1,5(a15,1x))') '#', 'E' , '<x>', '<y>', '<z>', '<f>'
        end select
        call messages_info(1, iunit, namespace=namespace)
      end if
 
    end if
 
    safe_allocate(psi_initial(1:mesh%np, 1:st%d%dim))
    safe_allocate(psi_ist(1:mesh%np, 1:st%d%dim))
 
    call states_elec_get_state(st, mesh, tpa_initialst, tpa_initialk, psi_initial)
 
    ! loop through every state
    do ist = 1, st%nst
 
      call states_elec_get_state(st, mesh, ist, tpa_initialk, psi_ist)
 
      ! final states are the unoccupied ones
      if (abs(st%occ(ist,tpa_initialk))  <  m_min_occ) then
 
        osc_strength = m_zero
        transition_energy = st%eigenval(ist, tpa_initialk) - st%eigenval(tpa_initialst, tpa_initialk)
 
        ! dipole matrix elements <f|x|i> etc. -> oscillator strengths
        do icoord = 1, space%dim    ! for x,y,z
 
          ff(1:mesh%np) = psi_initial(1:mesh%np, 1)*mesh%x(1:mesh%np, icoord)* &
            psi_ist(1:mesh%np, 1)
          osc(icoord)  = dmf_integrate(mesh, ff)
          osc_strength = osc_strength + m_two/real(space%dim, real64) *transition_energy*abs(osc(icoord))**2
 
        end do
 
        ! matrix elements <f|exp(iq.r)|i> -> dynamic structure factor
        if (use_qvector) then
 
          cff(1:mesh%np) = psi_initial(1:mesh%np, 1)*psi_ist(1:mesh%np, 1)
 
          do icoord = 1, space%dim    ! for x,y,z
            cff(1:mesh%np) = cff(1:mesh%np)*exp(m_zi*mesh%x(1:mesh%np, icoord)*qvector(icoord))
          end do
 
          dsf = abs(zmf_integrate(mesh, cff))**2
        end if
 
        ! write oscillator strengths (+ dynamic structure factor if qvector is given) into file
        if (mpi_grp_is_root(mpi_world)) then
 
          if (use_qvector) then
            write(message(1), '(1x,6(es15.8,1x))') units_from_atomic(units_out%energy, transition_energy), osc(:), osc_strength, &
              units_from_atomic(unit_one/units_out%energy, dsf)
          else
            write(message(1), '(1x,6(es15.8,1x))') units_from_atomic(units_out%energy, transition_energy), osc(:), osc_strength
          end if
 
          call messages_info(1, iunit, namespace=namespace)
 
        end if
 
      end if
 
    end do
 
    ! finally close the file
    if (mpi_grp_is_root(mpi_world)) then
      call io_close(iunit)
    end if
 
    safe_deallocate_a(psi_initial)
    safe_deallocate_a(psi_ist)
 
    safe_deallocate_a(ff)
    if (use_qvector) then
      safe_deallocate_a(cff)
    end if
    safe_deallocate_a(osc)
    if (use_qvector) then
      safe_deallocate_a(qvector)
    end if
 
    pop_sub(states_elec_write_tpa)
 
  end subroutine states_elec_write_tpa
 
 
  ! ---------------------------------------------------------
  !
  subroutine states_elec_write_bandstructure(dir, namespace, nst, st, ions, mesh, kpoints, &
    phase, vec_pot, vec_pot_var)
    character(len=*),             intent(in)  :: dir
    type(namespace_t),            intent(in)  :: namespace
    integer,                      intent(in)  :: nst
    type(states_elec_t),  target, intent(in)  :: st
    type(ions_t),         target, intent(in)  :: ions
    class(mesh_t),                intent(in)  :: mesh
    type(kpoints_t),              intent(in)  :: kpoints
    type(phase_t),                intent(in)  :: phase
    real(real64), optional, allocatable, intent(in)  :: vec_pot(:)
    !                                                                     !! dimension (1:mesh\%box\%dim)
    real(real64), optional, allocatable, intent(in)  :: vec_pot_var(:, :)
    !                                                                     !! dimensions (1:mesh\%box\%dim, 1:os\%sphere\%ns)
 
    integer :: idir, ist, ik, ns, is,npath, ib, ind
    integer, allocatable :: iunit(:)
    real(real64)   :: red_kpoint(mesh%box%dim)
    character(len=80) :: filename
 
    logical :: projection
    integer :: ii, ll, mm, nn, work, norb, work2
    integer :: ia, iorb, idim, maxnorb
    real(real64), allocatable :: ddot(:,:,:)
    complex(real64), allocatable :: zdot(:,:,:)
    real(real64), allocatable :: weight(:,:,:,:,:)
    type(orbitalset_t) :: os
    type(wfs_elec_t), pointer :: epsib
 
 
    push_sub(states_elec_write_bandstructure)
 
    ! shortcuts
    ns = 1
    if (st%d%nspin == 2) ns = 2
 
    !%Variable BandStructureComputeProjections
    !%Type logical
    !%Default false
    !%Section Output
    !%Description
    !% Determines if projections of wavefunctions on the atomic orbitals
    !% are computed or not for obtaining the orbital resolved band-structure.
    !%End
    call parse_variable(namespace, 'BandStructureComputeProjections', .false., projection)
 
 
    if (mpi_grp_is_root(mpi_world)) then
      safe_allocate(iunit(0:ns-1))
 
      do is = 0, ns-1
        if (ns > 1) then
          write(filename, '(a,i1.1,a)') 'bandstructure-sp', is+1
        else
          write(filename, '(a)') 'bandstructure'
        end if
        iunit(is) = io_open(trim(dir)//'/'//trim(filename), namespace, action='write')
 
        ! write header
        write(iunit(is),'(a)',advance='no') '# coord. '
        do idir = 1, mesh%box%dim
          write(iunit(is),'(3a)',advance='no') 'k', index2axis(idir), ' '
        end do
        if (.not. projection) then
          write(iunit(is),'(a,i6,3a)') '(red. coord.), bands:', nst, ' [', trim(units_abbrev(units_out%energy)), ']'
        else
          write(iunit(is),'(a,i6,3a)',advance='no') '(red. coord.), bands:', nst, ' [', trim(units_abbrev(units_out%energy)), '] '
          do ia = 1, ions%natoms
            work = orbitalset_utils_count(ions%atom(ia)%species)
            do norb = 1, work
              work2 = orbitalset_utils_count(ions%atom(ia)%species, norb)
              write(iunit(is),'(a, i3.3,a,i1.1,a)',advance='no') 'w(at=',ia,',os=',norb,') '
            end do
          end do
          write(iunit(is),'(a)') ''
        end if
      end do
    end if
 
    npath = kpoints%nkpt_in_path()*ns
 
 
    !We need to compute the projections of each wavefunctions on the localized basis
    if (projection) then
 
      maxnorb = 0
      do ia = 1, ions%natoms
        maxnorb = max(maxnorb, orbitalset_utils_count(ions%atom(ia)%species))
      end do
 
      safe_allocate(weight(1:st%nik,1:st%nst, 1:(2*max_l+1), 1:maxnorb, 1:ions%natoms))
      weight(1:st%nik,1:st%nst, 1:(2*max_l+1), 1:maxnorb, 1:ions%natoms) = m_zero
 
      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)
 
          !We count the orbitals
          work2 = 0
          do iorb = 1, ions%atom(ia)%species%get_niwfs()
            call ions%atom(ia)%species%get_iwf_ilm(iorb, 1, ii, ll, mm)
            call ions%atom(ia)%species%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(ions%atom(ia)%species, mesh, iorb, 1, &
                option__aotruncation__ao_full, 0.01_real64)
              work2 = work2 + 1
            end if
          end do
          os%norbs = work2
          os%ndim = 1
          os%use_submesh = .false.
          os%spec => ions%atom(ia)%species
 
          do iorb = 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, iorb, os%radius, os%ndim, use_mesh=.not.os%use_submesh, &
                normalize = .true.)
            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, iorb, os%radius, os%ndim, &
                use_mesh = .not. phase%is_allocated() .and. .not. os%use_submesh, &
                normalize = .true.)
            end if
          end do !iorb
 
          if (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 = max(pad_pow2(os%sphere%np), 1)
              if(.not. os%use_submesh) os%ldorbs = 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*os%norbs)
              end do
            end if
 
            call orbitalset_update_phase(os, mesh%box%dim, st%d%kpt, kpoints, (st%d%ispin==spin_polarized), &
              vec_pot, vec_pot_var)
          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
 
          do ik = st%d%kpt%start, st%d%kpt%end
            if (ik < st%nik - npath + 1) cycle ! We only want points inside the k-point path
 
            do ib = st%group%block_start, st%group%block_end
 
              if (phase%is_allocated()) then
                safe_allocate(epsib)
                call st%group%psib(ib, ik)%copy_to(epsib)
                call 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(ist, ind, iorb, norb, ia) = weight(ik, ind, iorb, norb, ia) + 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(ik, ind, iorb, norb, ia) = weight(ik, ind, iorb, norb, ia) + abs(zdot(idim, iorb, ist))**2
                    end do
                  end do
                end do
              end if
 
              if (phase%is_allocated()) then
                call epsib%end(copy=.false.)
                safe_deallocate_p(epsib)
              end if
 
            end do
          end do
 
          safe_deallocate_a(ddot)
          safe_deallocate_a(zdot)
 
          call orbitalset_end(os)
        end do !norb
 
        if (st%parallel_in_states .or. st%d%kpt%parallel) then
          call comm_allreduce(st%st_kpt_mpi_grp, weight(:,:,:,:,ia))
        end if
      end do !ia
 
    end if  !projection
 
    if (mpi_grp_is_root(mpi_world)) then
      ! output bands
      do ik = st%nik-npath+1, st%nik, ns
        do is = 0, ns - 1
          red_kpoint(1:mesh%box%dim) = kpoints%get_point(st%d%get_kpoint_index(ik + is), absolute_coordinates=.false.)
          write(iunit(is),'(1x)',advance='no')
          if (st%nik > npath) then
            write(iunit(is),'(f14.8)',advance='no') kpoints_get_path_coord(kpoints, &
              st%d%get_kpoint_index(ik + is)-st%d%get_kpoint_index(st%nik -npath))
          else
            write(iunit(is),'(f14.8)',advance='no') kpoints_get_path_coord(kpoints, &
              st%d%get_kpoint_index(ik + is))
          end if
          do idir = 1, mesh%box%dim
            write(iunit(is),'(f14.8)',advance='no') red_kpoint(idir)
          end do
          do ist = 1, nst
            write(iunit(is),'(1x,f14.8)',advance='no') units_from_atomic(units_out%energy, st%eigenval(ist, ik + is))
          end do
          if (projection) then
            do ia = 1, ions%natoms
              work = orbitalset_utils_count(ions%atom(ia)%species)
              do norb = 1, work
                work2 = orbitalset_utils_count(ions%atom(ia)%species, norb)
                do iorb = 1, work2
                  do ist = 1, nst
                    write(iunit(is),'(es15.8)',advance='no') weight(ik+is,ist,iorb,norb,ia)
                  end do
                end do
              end do
            end do
          end if
        end do
 
        do is = 0, ns-1
          write(iunit(is), '(a)')
        end do
      end do
 
      do is = 0, ns-1
        call io_close(iunit(is))
      end do
    end if
 
    if (projection) then
      safe_deallocate_a(weight)
    end if
 
    safe_deallocate_a(iunit)
 
 
    pop_sub(states_elec_write_bandstructure)
 
  end subroutine states_elec_write_bandstructure
 
end module states_elec_io_oct_m
 
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: