kdotp.F90

Go to the documentation of this file.

!!! Copyright (C) 2008-2010 David Strubbe
!!
!! 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 kdotp_oct_m
  use debug_oct_m
  use global_oct_m
  use grid_oct_m
  use output_oct_m
  use hamiltonian_elec_oct_m
  use io_oct_m
  use, intrinsic :: iso_fortran_env
  use kdotp_calc_oct_m
  use kpoints_oct_m
  use ks_potential_oct_m
  use lalg_adv_oct_m
  use linear_response_oct_m
  use mesh_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use mpi_oct_m
  use multisystem_basic_oct_m
  use namespace_oct_m
  use parser_oct_m
  use perturbation_oct_m
  use perturbation_kdotp_oct_m
  use perturbation_none_oct_m
  use profiling_oct_m
  use restart_oct_m
  use smear_oct_m
  use states_abst_oct_m
  use states_elec_oct_m
  use states_elec_dim_oct_m
  use states_elec_restart_oct_m
  use sternheimer_oct_m
  use electrons_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use utils_oct_m
  use v_ks_oct_m
 
  implicit none
 
  private
 
  public :: &
    kdotp_lr_run,       &
    int2str
 
  type kdotp_t
    private
    class(perturbation_t), pointer :: pert => null()
    class(perturbation_t), pointer :: pert2 => null()
 
    real(real64), allocatable :: eff_mass_inv(:,:,:,:)
    real(real64), allocatable :: velocity(:,:,:)
 
    type(lr_t), allocatable :: lr(:,:)
 
    type(lr_t), allocatable :: lr2(:,:,:)
 
    logical :: converged
    integer :: occ_solution_method
    real(real64)   :: degen_thres
    real(real64)   :: eta
  end type kdotp_t
 
contains
 
  ! ---------------------------------------------------------
  subroutine kdotp_lr_run(system, from_scratch)
    class(*),        intent(inout) :: system
    logical,         intent(in)    :: from_scratch
 
    push_sub(kdotp_lr_run)
 
    select type (system)
    class is (multisystem_basic_t)
      message(1) = "CalculationMode = kdotp not implemented for multi-system calculations"
      call messages_fatal(1, namespace=system%namespace)
    type is (electrons_t)
      call kdotp_lr_run_legacy(system, from_scratch)
    end select
 
    pop_sub(kdotp_lr_run)
  end subroutine kdotp_lr_run
 
  ! ---------------------------------------------------------
  subroutine kdotp_lr_run_legacy(sys, fromScratch)
    type(electrons_t),   intent(inout) :: sys
    logical,             intent(in)    :: fromScratch
 
    type(kdotp_t)           :: kdotp_vars
    type(sternheimer_t)     :: sh, sh2
    logical                 :: calc_eff_mass, calc_2nd_order, complex_response
 
    integer              :: idir, idir2, ierr, pdim, ispin
    character(len=100)   :: str_tmp
    real(real64)         :: errornorm
    type(restart_t)      :: restart_load, restart_dump
 
    class(perturbation_t), pointer :: pert2  ! for the second direction in second-order kdotp
 
    push_sub(kdotp_lr_run_legacy)
 
    call messages_experimental("k.p perturbation and calculation of effective masses", namespace=sys%namespace)
 
    if (sys%hm%pcm%run_pcm) then
      call messages_not_implemented("PCM for CalculationMode /= gs or td", namespace=sys%namespace)
    end if
 
    !TODO: This test belongs to the pert.F90 file, in the case of the velocity operator not been defined
    ! from the Hamiltonian
    ! In this case, there are other terms missing (MGGA, DFT+U for instance).
    if (sys%hm%theory_level == hartree_fock) then
      call messages_not_implemented('Commutator of Fock operator', namespace=sys%namespace)
    end if
    if (sys%hm%theory_level == generalized_kohn_sham_dft) then
      call messages_not_implemented('k.p with generalized Kohn-Sham DFT', namespace=sys%namespace)
    end if
 
 
    if (sys%kpoints%use_symmetries) then
      call messages_experimental("KPoints symmetries with CalculationMode = kdotp", namespace=sys%namespace)
    end if
 
    pdim = sys%space%periodic_dim
 
    if (.not. sys%space%is_periodic()) then
      message(1) = "k.p perturbation cannot be used for a finite system."
      call messages_fatal(1, namespace=sys%namespace)
    end if
 
    kdotp_vars%pert => perturbation_kdotp_t(sys%namespace, sys%ions)
    safe_allocate(kdotp_vars%lr(1, 1:pdim))
 
    call parse_input()
 
    if (calc_2nd_order) then
      kdotp_vars%pert2 => perturbation_none_t(sys%namespace)
      pert2 => perturbation_kdotp_t(sys%namespace, sys%ions)
      call kdotp_vars%pert2%setup_dir(1) ! direction is irrelevant
      safe_allocate(kdotp_vars%lr2(1, 1:pdim, 1:pdim))
    end if
 
    !Read ground-state wavefunctions
    complex_response = (abs(kdotp_vars%eta) > m_epsilon) .or. states_are_complex(sys%st)
    call restart_init(restart_load, sys%namespace, restart_gs, restart_type_load, sys%mc, ierr, mesh=sys%gr, exact=.true.)
    if (ierr == 0) then
      call states_elec_look_and_load(restart_load, sys%namespace, sys%space, sys%st, sys%gr, sys%kpoints, &
        is_complex = complex_response)
      call restart_end(restart_load)
    else
      message(1) = "A previous gs calculation is required."
      call messages_fatal(1, namespace=sys%namespace)
    end if
 
    ! Use of ForceComplex will make this true after states_elec_look_and_load even if it was not before.
    ! Otherwise, this line is a tautology.
    complex_response = states_are_complex(sys%st)
 
    ! Start restart. Note: we are going to use the same directory to read and write.
    ! Therefore, restart_dump must be initialized first to make sure the directory
    ! exists when we initialize restart_load.
    call restart_init(restart_dump, sys%namespace, restart_kdotp, restart_type_dump, sys%mc, ierr, mesh=sys%gr)
    ! no problem if this fails
    call restart_init(restart_load, sys%namespace, restart_kdotp, restart_type_load, sys%mc, ierr, mesh=sys%gr)
 
    ! setup Hamiltonian
    message(1) = 'Info: Setting up Hamiltonian for linear response.'
    call messages_info(1, namespace=sys%namespace)
    call v_ks_h_setup(sys%namespace, sys%space, sys%gr, sys%ions, sys%ext_partners, sys%st, sys%ks, sys%hm)
 
    if (states_are_real(sys%st)) then
      message(1) = 'Info: Using real wavefunctions.'
    else
      message(1) = 'Info: Using complex wavefunctions.'
    end if
    call messages_info(1, namespace=sys%namespace)
 
    ! Band velocities
    message(1) = 'Calculating band velocities.'
    call messages_info(1, namespace=sys%namespace)
 
    ! Compute band velocities
    safe_allocate(kdotp_vars%velocity(1:pdim, 1:sys%st%nst, 1:sys%st%nik))
    call calc_band_velocity(sys%namespace, sys%space, sys%gr, sys%st, sys%hm, kdotp_vars%pert, &
      kdotp_vars%velocity(:,:,:))
 
    ! Output of the band velocities
    if (mpi_grp_is_root(mpi_world)) then
      call io_mkdir(kdotp_dir, sys%namespace) ! data output
      call kdotp_write_band_velocity(sys%st, pdim, kdotp_vars%velocity(:,:,:), sys%namespace)
    end if
    safe_deallocate_a(kdotp_vars%velocity)
 
 
    call sternheimer_obsolete_variables(sys%namespace, 'KdotP_', 'KdotP')
    call sternheimer_init(sh, sys%namespace, sys%space, sys%gr, sys%st, sys%hm, sys%ks%xc, sys%mc, complex_response, &
      set_ham_var = 0, set_occ_response = (kdotp_vars%occ_solution_method == 0), &
      set_last_occ_response = (kdotp_vars%occ_solution_method == 0), occ_response_by_sternheimer = .true.)
    ! ham_var_set = 0 results in HamiltonianVariation = V_ext_only
    if (calc_2nd_order) then
      call sternheimer_init(sh2, sys%namespace, sys%space, sys%gr, sys%st, sys%hm, sys%ks%xc, sys%mc, complex_response, &
        set_ham_var = 0, set_occ_response = .false., set_last_occ_response = .false.)
    end if
 
    do idir = 1, pdim
      call lr_init(kdotp_vars%lr(1, idir))
      call lr_allocate(kdotp_vars%lr(1, idir), sys%st, sys%gr)
 
      if (calc_2nd_order) then
        do idir2 = idir, pdim
          call lr_init(kdotp_vars%lr2(1, idir, idir2))
          call lr_allocate(kdotp_vars%lr2(1, idir, idir2), sys%st, sys%gr)
        end do
      end if
 
      ! load wavefunctions
      if (.not. fromscratch) then
        str_tmp = kdotp_wfs_tag(idir)
        call restart_open_dir(restart_load, wfs_tag_sigma(sys%namespace, str_tmp, 1), ierr)
        if (ierr == 0) call states_elec_load(restart_load, sys%namespace, sys%space, sys%st, sys%gr, sys%kpoints, &
          ierr, lr=kdotp_vars%lr(1, idir))
        call restart_close_dir(restart_load)
 
        if (ierr /= 0) then
          message(1) = "Unable to read response wavefunctions from '"//trim(wfs_tag_sigma(sys%namespace, str_tmp, 1))//"'."
          call messages_warning(1, namespace=sys%namespace)
        end if
 
        if (calc_2nd_order) then
          do idir2 = idir, pdim
            str_tmp = kdotp_wfs_tag(idir, idir2)
            call restart_open_dir(restart_load, wfs_tag_sigma(sys%namespace, str_tmp, 1), ierr)
            if (ierr == 0) then
              call states_elec_load(restart_load, sys%namespace, sys%space, sys%st, sys%gr, sys%kpoints, &
                ierr, lr=kdotp_vars%lr2(1, idir, idir2))
            end if
            call restart_close_dir(restart_load)
 
            if (ierr /= 0) then
              message(1) = "Unable to read response wavefunctions from '"//trim(wfs_tag_sigma(sys%namespace, str_tmp, 1))//"'."
              call messages_warning(1, namespace=sys%namespace)
            end if
          end do
        end if
      end if
    end do
 
    call info()
    message(1) = "Info: Calculating k.p linear response of ground-state wavefunctions."
    call messages_info(1, namespace=sys%namespace)
    kdotp_vars%converged = .true.
 
    ! solve the Sternheimer equation
    do idir = 1, pdim
      write(message(1), '(3a)') 'Info: Calculating response for the ', index2axis(idir), &
        '-direction.'
      call messages_info(1, namespace=sys%namespace)
      call kdotp_vars%pert%setup_dir(idir)
 
      if (states_are_real(sys%st)) then
        call dsternheimer_solve(sh, sys%namespace, sys%space, sys%gr, sys%kpoints, sys%st, sys%hm, sys%mc, &
          kdotp_vars%lr(1:1, idir), 1, m_zero, kdotp_vars%pert, restart_dump, "", kdotp_wfs_tag(idir), &
          have_restart_rho = .false.)
        if (kdotp_vars%occ_solution_method == 1) then
          call dkdotp_add_occ(sys%namespace, sys%space, sys%gr, sys%st, sys%hm, kdotp_vars%pert, &
            kdotp_vars%lr(1, idir), kdotp_vars%degen_thres)
        end if
      else
        call zsternheimer_solve(sh, sys%namespace, sys%space, sys%gr, sys%kpoints, sys%st, sys%hm, sys%mc, &
          kdotp_vars%lr(1:1, idir), 1, m_zi * kdotp_vars%eta, kdotp_vars%pert, restart_dump, "", &
          kdotp_wfs_tag(idir), have_restart_rho = .false.)
        if (kdotp_vars%occ_solution_method == 1) then
          call zkdotp_add_occ(sys%namespace, sys%space, sys%gr, sys%st, sys%hm, kdotp_vars%pert, &
            kdotp_vars%lr(1, idir), kdotp_vars%degen_thres)
        end if
      end if
 
      kdotp_vars%converged = kdotp_vars%converged .and. sternheimer_has_converged(sh)
 
      errornorm = m_zero
      if (states_are_real(sys%st)) then
        call doutput_lr(sys%outp, sys%namespace, sys%space, kdotp_dir, sys%st, sys%gr, kdotp_vars%lr(1, idir), idir, 1, &
          sys%ions, units_out%force)
 
        do ispin = 1, sys%st%d%nspin
          errornorm = hypot(errornorm, dmf_nrm2(sys%gr, kdotp_vars%lr(1, idir)%ddl_rho(:, ispin)))
        end do
      else
        call zoutput_lr(sys%outp, sys%namespace, sys%space, kdotp_dir, sys%st, sys%gr, kdotp_vars%lr(1, idir), idir, 1, &
          sys%ions, units_out%force)
 
        do ispin = 1, sys%st%d%nspin
          errornorm = hypot(errornorm, zmf_nrm2(sys%gr, kdotp_vars%lr(1, idir)%zdl_rho(:, ispin)))
        end do
      end if
 
      write(message(1),'(a,g13.6)') "Norm of relative density variation = ", errornorm / sys%st%qtot
      call messages_info(1, namespace=sys%namespace)
 
      if (calc_2nd_order) then
        ! by equality of mixed partial derivatives, kdotp_vars%lr2(idir, idir2) = kdotp_vars%lr2(idir2, idir)
        do idir2 = idir, pdim
          write(message(1), '(3a)') 'Info: Calculating second-order response in the ', index2axis(idir2), &
            '-direction.'
          call messages_info(1, namespace=sys%namespace)
 
          call pert2%setup_dir(idir2)
 
          if (states_are_real(sys%st)) then
            call dsternheimer_solve_order2(sh, sh, sh2, sys%namespace, sys%space, sys%gr, sys%kpoints, sys%st, sys%hm, &
              sys%mc, kdotp_vars%lr(1:1, idir), kdotp_vars%lr(1:1, idir2), &
              1, m_zero, m_zero, kdotp_vars%pert, pert2, &
              kdotp_vars%lr2(1:1, idir, idir2), kdotp_vars%pert2, restart_dump, "", kdotp_wfs_tag(idir, idir2), &
              have_restart_rho = .false., have_exact_freq = .true.)
          else
            call zsternheimer_solve_order2(sh, sh, sh2, sys%namespace, sys%space, sys%gr, sys%kpoints, sys%st, sys%hm, &
              sys%mc, kdotp_vars%lr(1:1, idir), kdotp_vars%lr(1:1, idir2), &
              1, m_zi * kdotp_vars%eta, m_zi * kdotp_vars%eta, kdotp_vars%pert, pert2, &
              kdotp_vars%lr2(1:1, idir, idir2), kdotp_vars%pert2, restart_dump, "", kdotp_wfs_tag(idir, idir2), &
              have_restart_rho = .false., have_exact_freq = .true.)
          end if
 
        end do
        message(1) = ""
        call messages_info(1, namespace=sys%namespace)
      end if
    end do ! idir
 
    ! calculate effective masses
    if (calc_eff_mass) then
      message(1) = "Info: Calculating effective masses."
      call messages_info(1, namespace=sys%namespace)
 
      safe_allocate(kdotp_vars%eff_mass_inv(1:pdim, 1:pdim, 1:sys%st%nst, 1:sys%st%nik))
      kdotp_vars%eff_mass_inv(:,:,:,:) = m_zero
 
 
      if (states_are_real(sys%st)) then
        call dcalc_eff_mass_inv(sys%namespace, sys%space, sys%gr, sys%st, sys%hm, kdotp_vars%lr, &
          kdotp_vars%pert, kdotp_vars%eff_mass_inv, kdotp_vars%degen_thres)
      else
        call zcalc_eff_mass_inv(sys%namespace, sys%space, sys%gr, sys%st, sys%hm, kdotp_vars%lr, &
          kdotp_vars%pert, kdotp_vars%eff_mass_inv, kdotp_vars%degen_thres)
      end if
 
      call kdotp_write_degeneracies(sys%st, kdotp_vars%degen_thres, sys%namespace)
      call kdotp_write_eff_mass(sys%st, sys%kpoints, kdotp_vars, sys%namespace, sys%space%periodic_dim)
 
      safe_deallocate_a(kdotp_vars%eff_mass_inv)
    end if
 
    ! clean up some things
    do idir = 1, pdim
      call lr_dealloc(kdotp_vars%lr(1, idir))
 
      if (calc_2nd_order) then
        do idir2 = idir, pdim
          call lr_dealloc(kdotp_vars%lr2(1, idir, idir2))
        end do
      end if
    end do
 
    call restart_end(restart_load)
    call restart_end(restart_dump)
    call sternheimer_end(sh)
    deallocate(kdotp_vars%pert)
    safe_deallocate_a(kdotp_vars%lr)
 
    if (calc_2nd_order) then
      call sternheimer_end(sh2)
      safe_deallocate_p(pert2)
      deallocate(kdotp_vars%pert2)
      safe_deallocate_a(kdotp_vars%lr2)
    end if
 
    call states_elec_deallocate_wfns(sys%st)
 
    pop_sub(kdotp_lr_run_legacy)
 
  contains
 
    ! ---------------------------------------------------------
 
    subroutine parse_input()
 
      push_sub(kdotp_lr_run_legacy.parse_input)
 
      !%Variable KdotPOccupiedSolutionMethod
      !%Type integer
      !%Default sternheimer_eqn
      !%Section Linear Response::KdotP
      !%Description
      !% Method of calculating the contribution of the projection of the
      !% linear-response wavefunctions in the occupied subspace.
      !%Option sternheimer_eqn 0
      !% The Sternheimer equation is solved including the occupied subspace,
      !% to get the full linear-response wavefunctions.
      !%Option sum_over_states 1
      !% The Sternheimer equation is solved only in the unoccupied subspace,
      !% and a sum-over-states perturbation-theory expression is used to
      !% evaluate the contributions in the occupied subspace.
      !%End
 
      call messages_obsolete_variable(sys%namespace, 'KdotP_OccupiedSolutionMethod', 'KdotPOccupiedSolutionMethod')
 
      call parse_variable(sys%namespace, 'KdotPOccupiedSolutionMethod', 0, kdotp_vars%occ_solution_method)
      if (kdotp_vars%occ_solution_method == 1 .and. .not. smear_is_semiconducting(sys%st%smear)) then
        call messages_not_implemented('KdotPOccupiedSolutionMethod = sum_over_states for non-semiconducting smearing', &
          namespace=sys%namespace)
      end if
 
      !%Variable DegeneracyThreshold
      !%Type float
      !%Default 1e-5
      !%Section States
      !%Description
      !% States with energy <math>E_i</math> and <math>E_j</math> will be considered degenerate
      !% if <math> \left| E_i - E_j \right| < </math><tt>DegeneracyThreshold</tt>.
      !%End
      call parse_variable(sys%namespace, 'DegeneracyThreshold', units_from_atomic(units_inp%energy, 1e-5_real64), &
        kdotp_vars%degen_thres)
      kdotp_vars%degen_thres = units_to_atomic(units_inp%energy, kdotp_vars%degen_thres)
 
      !%Variable KdotPEta
      !%Type float
      !%Default 0.0
      !%Section Linear Response::KdotP
      !%Description
      !% Imaginary frequency added to Sternheimer equation which may improve convergence.
      !% Not recommended.
      !%End
      call messages_obsolete_variable(sys%namespace, 'KdotP_Eta', 'KdotPEta')
      call parse_variable(sys%namespace, 'KdotPEta', m_zero, kdotp_vars%eta)
      kdotp_vars%eta = units_to_atomic(units_inp%energy, kdotp_vars%eta)
 
      !%Variable KdotPCalculateEffectiveMasses
      !%Type logical
      !%Default true
      !%Section Linear Response::KdotP
      !%Description
      !% If true, uses <tt>kdotp</tt> perturbations of ground-state wavefunctions
      !% to calculate effective masses. It is not correct for degenerate states.
      !%End
      call messages_obsolete_variable(sys%namespace, 'KdotP_CalculateEffectiveMasses', 'KdotPCalculateEffectiveMasses')
      call parse_variable(sys%namespace, 'KdotPCalculateEffectiveMasses', .true., calc_eff_mass)
 
      !%Variable KdotPCalcSecondOrder
      !%Type logical
      !%Default false
      !%Section Linear Response::KdotP
      !%Description
      !% If true, calculates second-order response of wavefunctions as well as first-order response.
      !% Note that the second derivative of the Hamiltonian is NOT included in this calculation.
      !% This is needed for a subsequent run in <tt>CalculationMode = em_resp</tt> with <tt>EMHyperpol</tt>.
      !%End
      call parse_variable(sys%namespace, 'KdotPCalcSecondOrder', .false., calc_2nd_order)
 
      pop_sub(kdotp_lr_run_legacy.parse_input)
 
    end subroutine parse_input
 
    ! ---------------------------------------------------------
    subroutine info()
 
      push_sub(kdotp_lr_run_legacy.info)
 
      call kdotp_vars%pert%info()
 
      call messages_print_with_emphasis(msg='k.p perturbation theory', namespace=sys%namespace)
 
      if (kdotp_vars%occ_solution_method == 0) then
        message(1) = 'Occupied solution method: Sternheimer equation.'
      else
        message(1) = 'Occupied solution method: sum over states.'
      end if
 
      call messages_info(1, namespace=sys%namespace)
 
      call messages_print_with_emphasis(namespace=sys%namespace)
 
      pop_sub(kdotp_lr_run_legacy.info)
 
    end subroutine info
 
  end subroutine kdotp_lr_run_legacy
 
  ! ---------------------------------------------------------
  subroutine kdotp_write_band_velocity(st, periodic_dim, velocity, namespace)
    type(states_elec_t), intent(inout) :: st
    integer,             intent(in)    :: periodic_dim
    real(real64),        intent(in)    :: velocity(:,:,:)
    type(namespace_t),   intent(inout) :: namespace
 
    character(len=80) :: filename, tmp
    integer :: iunit, ik, ist, ik2, ispin, idir
 
    if (.not. mpi_grp_is_root(mpi_world)) return ! only first node outputs
 
    push_sub(kdotp_write_band_velocity)
 
    write(filename, '(a)') kdotp_dir//'velocity'
    iunit = io_open(trim(filename), namespace, action='write')
    write(iunit,'(a)') '# Band velocities'
 
    do ik = 1, st%nik
      ispin = st%d%get_spin_index(ik)
      ik2 = st%d%get_kpoint_index(ik)
      tmp = int2str(ik2)
 
      write(iunit,'(a,i1,a,a)') '# spin = ', ispin, ', k-point = ', trim(tmp)
 
      write(iunit,'(a)',advance='no') '# state    energy       '
      do idir = 1, periodic_dim
        write(iunit,'(3a)',advance='no') 'vg(', trim(index2axis(idir)), ')       '
      end do
      write(iunit,'(a)')
 
      write(iunit,'(3a)',advance='no')       '#           [', trim(units_abbrev(units_out%energy)), ']     '
      do idir = 1, periodic_dim
        write(iunit,'(3a)',advance='no') '[', trim(units_abbrev(units_out%velocity)), '] '
      end do
      write(iunit,'(a)')
 
      do ist = 1, st%nst
        write(iunit,'(i5,f12.5,3f12.5)') ist, units_from_atomic(units_out%energy, st%eigenval(ist, ik)), &
          velocity(1:periodic_dim, ist, ik)
      end do
    end do
 
    call io_close(iunit)
    pop_sub(kdotp_write_band_velocity)
  end subroutine kdotp_write_band_velocity
 
  ! ---------------------------------------------------------
  subroutine kdotp_write_eff_mass(st, kpoints, kdotp_vars, namespace, periodic_dim)
    type(states_elec_t),  intent(inout) :: st
    type(kpoints_t),      intent(in)    :: kpoints
    type(kdotp_t),        intent(inout) :: kdotp_vars
    type(namespace_t),    intent(in)    :: namespace
    integer,              intent(in)    :: periodic_dim
 
    character(len=80) :: filename, tmp
    integer :: iunit, ik, ist, ik2, ispin
 
    if (.not. mpi_grp_is_root(mpi_world)) return ! only first node outputs
 
    push_sub(kdotp_write_eff_mass)
 
    do ik = 1, st%nik
      ispin = st%d%get_spin_index(ik)
      ik2 = st%d%get_kpoint_index(ik)
 
      tmp = int2str(ik2)
      write(filename, '(3a, i1)') kdotp_dir//'kpoint_', trim(tmp), '_', ispin
      iunit = io_open(trim(filename), namespace, action='write')
 
      write(iunit,'(a, i10)')    '# spin    index = ', ispin
      write(iunit,'(a, i10)')    '# k-point index = ', ik2
      write(iunit,'(a, 99f12.8)') '# k-point coordinates = ', kpoints%get_point(ik2)
      if (.not. kdotp_vars%converged) write(iunit, '(a)') "# WARNING: not converged"
 
      write(iunit,'(a)')
      write(iunit,'(a)') '# Inverse effective-mass tensors'
      do ist = 1, st%nst
        write(iunit,'(a)')
        tmp = int2str(ist)
        write(iunit,'(a, a, a, f12.8, a, a)') 'State #', trim(tmp), ', Energy = ', &
          units_from_atomic(units_out%energy, st%eigenval(ist, ik)), ' ', units_abbrev(units_out%energy)
        call output_tensor(kdotp_vars%eff_mass_inv(:, :, ist, ik), periodic_dim, unit_one, iunit=iunit)
      end do
 
      write(iunit,'(a)')
      write(iunit,'(a)') '# Effective-mass tensors'
      do ist = 1, st%nst
        write(iunit,'(a)')
        tmp = int2str(ist)
        write(iunit,'(a, a, a, f12.8, a, a)') 'State #', trim(tmp), ', Energy = ', &
          units_from_atomic(units_out%energy, st%eigenval(ist, ik)), ' ', units_abbrev(units_out%energy)
        call lalg_inverse(periodic_dim, kdotp_vars%eff_mass_inv(:, :, ist, ik), 'dir')
        call output_tensor(kdotp_vars%eff_mass_inv(:, :, ist, ik), periodic_dim, unit_one, iunit=iunit)
      end do
 
      call io_close(iunit)
    end do
 
    pop_sub(kdotp_write_eff_mass)
  end subroutine kdotp_write_eff_mass
 
  ! ---------------------------------------------------------
  subroutine kdotp_write_degeneracies(st, threshold, namespace)
    type(states_elec_t), intent(inout) :: st
    real(real64),        intent(in)    :: threshold
    type(namespace_t),   intent(in)    :: namespace
 
    character(len=80) :: tmp
    integer :: ik, ist, ist2, ik2, ispin
 
    if (.not. mpi_grp_is_root(mpi_world)) return ! only first node outputs
 
    push_sub(kdotp_write_degeneracies)
 
    call messages_print_with_emphasis(msg='Degenerate subspaces', namespace=namespace)
 
    do ik = 1, st%nik
      ispin = st%d%get_spin_index(ik)
      ik2 = st%d%get_kpoint_index(ik)
 
      tmp = int2str(ik2)
      write(message(1), '(3a, i1)') 'k-point ', trim(tmp), ', spin ', ispin
      call messages_info(1, namespace=namespace)
 
      ist = 1
      do while (ist <= st%nst)
        ! test for degeneracies
        write(message(1),'(a)') '===='
        tmp = int2str(ist)
        write(message(2),'(a, a, a, f12.8, a, a)') 'State #', trim(tmp), ', Energy = ', &
          units_from_atomic(units_out%energy, st%eigenval(ist, ik)), ' ', units_abbrev(units_out%energy)
        call messages_info(2, namespace=namespace)
 
        ist2 = ist + 1
        do while (ist2 <= st%nst .and. &
          abs(st%eigenval(min(ist2, st%nst), ik) - st%eigenval(ist, ik)) < threshold)
          tmp = int2str(ist2)
          write(message(1),'(a, a, a, f12.8, a, a)') 'State #', trim(tmp), ', Energy = ', &
            units_from_atomic(units_out%energy, st%eigenval(ist2, ik)), ' ', units_abbrev(units_out%energy)
          call messages_info(1, namespace=namespace)
          ist2 = ist2 + 1
        end do
 
        ist = ist2
      end do
 
      write(message(1),'()')
      call messages_info(1, namespace=namespace)
    end do
 
    message(1) = "Velocities and effective masses are not correct within degenerate subspaces."
    call messages_warning(1, namespace=namespace)
 
    pop_sub(kdotp_write_degeneracies)
  end subroutine kdotp_write_degeneracies
 
  ! ---------------------------------------------------------
  character(len=12) pure function int2str(ii) result(str)
    integer, intent(in) :: ii
 
    write(str, '(i11)') ii
    str = trim(adjustl(str))
 
  end function int2str
 
end module kdotp_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: