output_me.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!! Copyright (C) 2023 F. Troisi
!!
!! 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 output_me_oct_m
  use debug_oct_m
  use elec_matrix_elements_oct_m
  use electron_space_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_oct_m
  use io_function_oct_m
  use io_oct_m
  use ions_oct_m
  use kpoints_oct_m
  use mpi_oct_m
  use mpi_lib_oct_m
  use messages_oct_m
  use namespace_oct_m
  use output_low_oct_m
  use parser_oct_m
  use profiling_oct_m
  use singularity_oct_m
  use space_oct_m
  use states_abst_oct_m
  use states_elec_oct_m
  use states_elec_dim_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use utils_oct_m
  use varinfo_oct_m
  use xc_oct_m
 
  implicit none
 
  private
  public ::           &
    output_me_init,   &
    output_me
 
contains
 
  ! ---------------------------------------------------------
  subroutine output_me_init(this, namespace, space, st, gr, nst)
    type(output_me_t),   intent(out) :: this
    type(namespace_t),   intent(in)  :: namespace
    class(space_t),      intent(in)  :: space
    type(states_elec_t), intent(in)  :: st
    type(grid_t),        intent(in)  :: gr
    integer,             intent(in)  :: nst
 
    integer(int64) :: how(0:MAX_OUTPUT_TYPES)
    integer :: output_interval(0:MAX_OUTPUT_TYPES)
 
    push_sub(output_me_init)
 
    !%Variable OutputMatrixElements
    !%Type block
    !%Default none
    !%Section Output
    !%Description
    !% Specifies what matrix elements to print.
    !% Enabled only if <tt>Output</tt> block includes <tt>matrix_elements</tt>.
    !% The output files go into the <tt>static</tt> directory, except when
    !% running a time-dependent simulation, when the directory <tt>td.XXXXXXX</tt> is used.
    !%
    !% Example:
    !% <br><br><tt>%OutputMatrixElements
    !% <br>&nbsp;&nbsp;momentum
    !% <br>&nbsp;&nbsp;ks_multipoles
    !% <br>%<br></tt>
    !%
    !% It is possible to specify only compute the matrix elements for some of the states
    !% using the variables <tt>OutptMEStart</tt> and <tt>OutputMEEnd</tt>.
    !%Option momentum 1
    !% Momentum. Filename: <tt>ks_me_momentum</tt>.
    !% This does not include the nonlocal part of the pseudopotential, nor any nonlocal operator.
    !% These are included in the velocity operator.
    !% In case of a vector potential, this returns the matrix element of the canonical momentum
    !% <math>\hat{\boldsymbol{\Pi}} = \hat{\mathbf{p}}-\mathbf{A}</math>.
    !%Option ang_momentum 2
    !% Dimensionless angular momentum <math>\vec{r} \times \vec{k}</math>. Filename: <tt>ks_me_angular_momentum</tt>.
    !%Option one_body 3
    !% <math>\left< i \left| \hat{T} + V_{ext} \right| j \right></math>. Not available with states parallelization.
    !%Option two_body 4
    !% <math>\left< ij \left| \frac{1}{\left|\vec{r}_1-\vec{r}_2\right|} \right| kl \right></math>.
    !% Not available with states parallelization.
    !% Not available with states parallelization. For periodic system, this is not available for k-point parallelization neither.
    !%Option two_body_exc_k 5
    !% <math>\left< n1-k1, n2-k2 \left| \frac{1}{\left|\vec{r}_1-\vec{r}_2\right|} \right| n2-k1 n1-k2 \right></math>.
    !% Not available with states parallelization. For periodic system, this is not available for k-point parallelization neither.
    !%Option ks_multipoles 6
    !% See <tt>OutputMEMultipoles</tt>. Not available with states parallelization.
    !%Option dipole 7
    !% Prints the dipole matrix elements. Not available with states parallelization.
    !% For periodic systems, the intraband terms (dipole matrix elements between degenerated states)
    !% are set to zero, and only the absolute value of the dipole matrix element is printed.
    !% Not yet supported for spinors.
    !%Option velocity 8
    !% Velocity matrix elements, defined as the matrix elements of <math>\hat{v} = \hat{p}/m -i/\hbar[\hat{V}_{NL},\vec{r}]</math>.
    !%End
 
    this%what = .false.
    call io_function_read_what_how_when(namespace, space, this%what, how, output_interval, &
      'OutputMatrixElements')
 
    if (st%parallel_in_states) then
      if (this%what(option__outputmatrixelements__two_body)) then
        call messages_not_implemented("OutputMatrixElements=two_body is not implemented in states parallelization.", &
          namespace=namespace)
      end if
      if (this%what(option__outputmatrixelements__dipole)) then
        call messages_not_implemented("OutputMatrixElements=dipole is not implemented in states parallelization.", &
          namespace=namespace)
      end if
    end if
 
    if (space%dim /= 3) this%what(option__outputmatrixelements__ang_momentum) = .false.
 
    if (this%what(option__outputmatrixelements__ks_multipoles)) then
      !%Variable OutputMEMultipoles
      !%Type integer
      !%Default 1
      !%Section Output
      !%Description
      !% This variable decides which multipole moments are printed out for
      !% <tt>OutputMatrixElements = ks_multipoles</tt>:
      !%
      !% In 3D, if, for example, <tt>OutputMEMultipoles = 1</tt>, then the program will print three
      !% files, <tt>ks_me_multipoles.x</tt> (<tt>x</tt>=1,2,3), containing
      !% respectively the (1,-1), (1,0) and (1,1) multipole matrix elements
      !% between Kohn-Sham states.
      !%
      !% In 2D, this variable is ignored: it will always print two files,
      !% <tt>ks_me_multipoles.i</tt> (<tt>i</tt>=1,2), containing the <math>x</math> and
      !% <math>y</math> dipole matrix elements.
      !%
      !% In 1D, if, for example, <tt>OutputMEMultipoles = 2</tt>, the program will print two files, containing the
      !% <math>x</math> and <math>x^2</math> matrix elements between Kohn-Sham states.
      !%End
      call parse_variable(namespace, 'OutputMEMultipoles', 1, this%ks_multipoles)
    end if
 
    !%Variable OutputMEStart
    !%Type integer
    !%Default 1
    !%Section Output
    !%Description
    !% Specifies the state/band index for starting to compute the matrix element.
    !% So far, this is only used for dipole matrix elements.
    !%End
    call parse_variable(namespace, 'OutputMEStart', 1, this%st_start)
    assert(this%st_start > 0 .and. this%st_start <= nst)
 
    !%Variable OutputMEEnd
    !%Type integer
    !%Default 1
    !%Section Output
    !%Description
    !% Specifies the highest state/band index used to compute the matrix element.
    !% So far, this is only used for dipole matrix elements.
    !%End
    call parse_variable(namespace, 'OutputMEEnd', nst, this%st_end)
    if (this%st_end <= 0 .or. this%st_end > nst .or. this%st_start > this%st_end) then
      message(1) = "OutputMEEnd must be positive, smaller or equal to the number of states, and larger than OutputMEStart"
      call messages_fatal(1, namespace=namespace)
    end if
    this%nst = this%st_end - this%st_start +1
 
    pop_sub(output_me_init)
  end subroutine output_me_init
 
 
  ! ---------------------------------------------------------
  subroutine output_me(this, namespace, space, dir, st, gr, ions, hm)
    type(output_me_t),        intent(in)    :: this
    type(namespace_t),        intent(in)    :: namespace
    class(space_t),           intent(in)    :: space
    character(len=*),         intent(in)    :: dir
    type(states_elec_t),      intent(inout) :: st
    type(grid_t),             intent(in)    :: gr
    type(ions_t),             intent(in)    :: ions
    type(hamiltonian_elec_t), intent(inout) :: hm
 
    integer :: id, ll, iunit
    character(len=256) :: fname
    real(real64), allocatable :: doneint(:), dtwoint(:)
    complex(real64), allocatable :: zoneint(:), ztwoint(:)
    integer, allocatable :: iindex(:,:), jindex(:,:), kindex(:,:), lindex(:,:)
    type(singularity_t) :: singul
 
    push_sub(output_me)
 
    if (this%what(option__outputmatrixelements__momentum)) then
      write(fname,'(2a)') trim(dir), '/ks_me_momentum'
      call output_me_out_momentum(gr, st, space, fname, namespace, hm, ions)
    end if
 
    if (this%what(option__outputmatrixelements__velocity)) then
      write(fname,'(2a)') trim(dir), '/ks_me_velocity'
      call output_me_out_momentum(gr, st, space, fname, namespace, hm, ions, non_local=.true.)
    end if
 
    if (this%what(option__outputmatrixelements__ang_momentum)) then
      write(fname,'(2a)') trim(dir), '/ks_me_angular_momentum'
      call output_me_out_ang_momentum(gr, st, space, fname, namespace, hm%kpoints)
    end if
 
    if (this%what(option__outputmatrixelements__ks_multipoles)) then
      call output_me_out_ks_multipoles(gr, st, space, dir, this, namespace)
    end if
 
    if (this%what(option__outputmatrixelements__dipole)) then
      assert(.not. st%parallel_in_states)
      call output_me_out_dipole(gr, st, space, dir, namespace, hm, ions, this%st_start, this%st_end)
    end if
 
    if (this%what(option__outputmatrixelements__one_body)) then
      message(1) = "Computing one-body matrix elements"
      call messages_info(1, namespace=namespace)
 
      if (st%parallel_in_states) call messages_not_implemented("OutputMatrixElements=one_body with states parallelization", &
        namespace=namespace)
      if (st%d%kpt%parallel) call messages_not_implemented("OutputMatrixElements=one_body with k-points parallelization", &
        namespace=namespace)
      if (family_is_mgga_with_exc(hm%xc)) then
        call messages_not_implemented("OutputMatrixElements=one_body with MGGA", namespace=namespace)
      end if
      ! how to do this properly? states_elec_matrix
      iunit = io_open(trim(dir)//'/output_me_one_body', namespace, action='write')
 
      id = st%nst*(st%nst+1)/2
 
      safe_allocate(iindex(1:id,1))
      safe_allocate(jindex(1:id,1))
 
      if (states_are_real(st)) then
        safe_allocate(doneint(1:id))
        call one_body_me(gr, st, namespace, hm, iindex(:,1), jindex(:,1), doneint)
        if (gr%mpi_grp%is_root()) then
          do ll = 1, id
            write(iunit, *) iindex(ll,1), jindex(ll,1), doneint(ll)
          end do
        end if
        safe_deallocate_a(doneint)
      else
        safe_allocate(zoneint(1:id))
        call one_body_me(gr, st, namespace, hm, iindex(:,1), jindex(:,1), zoneint)
        if (gr%mpi_grp%is_root()) then
          do ll = 1, id
            write(iunit, *) iindex(ll,1), jindex(ll,1), zoneint(ll)
          end do
        end if
        safe_deallocate_a(zoneint)
      end if
 
      safe_deallocate_a(iindex)
      safe_deallocate_a(jindex)
      call io_close(iunit)
    end if
 
    if (this%what(option__outputmatrixelements__two_body) .or. this%what(option__outputmatrixelements__two_body_exc_k)) then
      message(1) = "Computing two-body matrix elements"
      call messages_info(1, namespace=namespace)
 
      assert(.not. st%parallel_in_states)
      ! two_body_me needs the exchange Poisson solver in double precision
      if (hm%exxop%compute_single_precision) then
        call messages_not_implemented("OutputMatrixElements=two_body with ExchangeSinglePrecision", namespace=namespace)
      end if
      if (this%what(option__outputmatrixelements__two_body)) then
        if (st%parallel_in_states)  then
          call messages_not_implemented("OutputMatrixElements=two_body with states parallelization", namespace=namespace)
        end if
        if (st%d%kpt%parallel) then
          call messages_not_implemented("OutputMatrixElements=two_body with k-points parallelization", namespace=namespace)
        end if
        ! how to do this properly? states_matrix
        iunit = io_open(trim(dir)//'/output_me_two_body', namespace, action='write')
        if (gr%mpi_grp%is_root()) then
          write(iunit, '(a)') '# (n1, k1)  (n2, k2)  (n3, k3)  (n4, k4)  (n1-k1, n2-k2|n3-k3, n4-k4)'
        end if
      else
        if (st%parallel_in_states) then
          call messages_not_implemented("OutputMatrixElements=two_body_exc_k with states parallelization", namespace=namespace)
        end if
        if (st%d%kpt%parallel) then
          call messages_not_implemented("OutputMatrixElements=two_body_exc_k with k-points parallelization", namespace=namespace)
        end if
        ! how to do this properly? states_matrix
        iunit = io_open(trim(dir)//'/output_me_two_body_density', namespace, action='write')
        if (gr%mpi_grp%is_root()) then
          write(iunit, '(a)') '#(n1, k1)  (n2, k2)  (n1-k1, n1-k2|n2-k2, n2-k1)'
        end if
      end if
 
      if (this%what(option__outputmatrixelements__two_body)) then
        if (states_are_real(st)) then
          id = st%nik*this%nst*(st%nik*this%nst+1)*(st%nik**2*this%nst**2+st%nik*this%nst+2)/8
        else
          id = (st%nik*this%nst)**4
        end if
      else
        id = (st%nik*this%nst)**2
      end if
 
      if (states_are_complex(st)) then
        call singularity_init(singul, namespace, space, st, hm%kpoints)
      end if
 
      safe_allocate(iindex(1:2, 1:id))
      safe_allocate(jindex(1:2, 1:id))
      safe_allocate(kindex(1:2, 1:id))
      safe_allocate(lindex(1:2, 1:id))
 
      if (states_are_real(st)) then
        safe_allocate(dtwoint(1:id))
        call two_body_me(gr, st, space, namespace, hm%kpoints, hm%exxop%psolver, this%st_start, this%st_end, &
          iindex, jindex, kindex, lindex, dtwoint)
        if (gr%mpi_grp%is_root()) then
          do ll = 1, id
            write(iunit, '(4(i4,i5,1x),es22.12)') iindex(1:2,ll), jindex(1:2,ll), kindex(1:2,ll), lindex(1:2,ll), dtwoint(ll)
          end do
        end if
        safe_deallocate_a(dtwoint)
      else
        safe_allocate(ztwoint(1:id))
        if (hm%phase%is_allocated()) then
          !We cannot pass the phase array like that if kpt%start is not 1.
          assert(.not. st%d%kpt%parallel)
          call two_body_me(gr, st, space, namespace, hm%kpoints, hm%exxop%psolver, this%st_start, this%st_end, iindex, jindex, &
            kindex, lindex, ztwoint, phase = hm%phase, singularity = singul, &
            exc_k = (this%what(option__outputmatrixelements__two_body_exc_k)))
        else
          call two_body_me(gr, st, space, namespace, hm%kpoints, hm%exxop%psolver, this%st_start, this%st_end, &
            iindex, jindex, kindex, lindex, ztwoint, exc_k = (this%what(option__outputmatrixelements__two_body_exc_k)))
        end if
 
        if (gr%mpi_grp%is_root()) then
          if (this%what(option__outputmatrixelements__two_body)) then
            do ll = 1, id
              write(iunit, '(4(i4,i5,1x),2es22.12)') iindex(1:2,ll), jindex(1:2,ll), kindex(1:2,ll), lindex(1:2,ll), ztwoint(ll)
            end do
          else
            do ll = 1, id
              write(iunit, '(2(i4,i5),2es22.12)') iindex(1:2,ll), kindex(1:2,ll), ztwoint(ll)
            end do
          end if
        end if
        safe_deallocate_a(ztwoint)
      end if
 
      safe_deallocate_a(iindex)
      safe_deallocate_a(jindex)
      safe_deallocate_a(kindex)
      safe_deallocate_a(lindex)
      call io_close(iunit)
 
      if (states_are_complex(st)) then
        call singularity_end(singul)
      end if
    end if
 
    pop_sub(output_me)
  end subroutine output_me
 
 
  ! ---------------------------------------------------------
  subroutine output_me_out_momentum(gr, st, space, fname, namespace, hm, ions, non_local)
    type(grid_t),                 intent(in) :: gr
    type(states_elec_t),          intent(in) :: st
    type(space_t),                intent(in) :: space
    character(len=*),             intent(in) :: fname
    type(namespace_t),            intent(in) :: namespace
    type(hamiltonian_elec_t),     intent(in) :: hm
    type(ions_t),                 intent(in) :: ions
    logical, optional,            intent(in) :: non_local
 
    integer            :: ik, ist, is, ns, iunit, idir, dim
    character(len=80)  :: cspin, str_tmp
    real(real64)       :: kpoint(space%dim)
    real(real64), allocatable :: momentum(:,:,:)
 
    push_sub(output_me_out_momentum)
 
    dim = space%dim
 
    safe_allocate(momentum(1:dim, 1:st%nst, 1:st%nik))
 
    ! Compute matrix elements
    call elec_momentum_me(namespace, gr, st, space, hm, ions, momentum, non_local)
 
    iunit = io_open(fname, namespace, action='write')
 
    ns = 1
    if (st%d%nspin == 2) ns = 2
 
    if (optional_default(non_local, .false.)) then
      write(message(1),'(a)') 'Velocity of the KS states [a.u.]:'
    else
      write(message(1),'(a)') 'Momentum of the KS states [a.u.]:'
    end if
    call messages_info(1, iunit)
    if (st%nik > ns) then
      message(1) = 'k-points [' // trim(units_abbrev(unit_one/units_out%length)) // ']'
      call messages_info(1, iunit)
    end if
 
    do ik = 1, st%nik, ns
      kpoint(:) = hm%kpoints%get_point(st%d%get_kpoint_index(ik))
 
      if (st%nik > ns) then
        write(message(1), '(a,i4, a)') '#k =', ik, ', k = ('
        do idir = 1, dim
          write(str_tmp, '(f12.6)') units_from_atomic(unit_one/units_out%length, kpoint(idir))
          message(1) = trim(message(1)) // trim(str_tmp)
          if (idir == dim) then
            message(1) = trim(message(1)) // ")"
          else
            message(1) = trim(message(1)) // ","
          end if
        end do
        call messages_info(1, iunit)
      end if
 
      write(message(1), '(a4,1x,a5)') '#st',' Spin'
      do idir = 1, dim
        if (optional_default(non_local, .false.)) then
          write(str_tmp, '(a,a1,a)') '        <v', index2axis(idir), '>'
        else
          write(str_tmp, '(a,a1,a)') '        <p', index2axis(idir), '>'
        end if
        message(1) = trim(message(1)) // trim(str_tmp)
      end do
      write(str_tmp, '(4x,a12,1x)') 'Occupation '
      message(1) = trim(message(1)) // trim(str_tmp)
      call messages_info(1, iunit)
 
      do ist = 1, st%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(message(1), '(i4,3x,a2,1x)') ist, trim(cspin)
          do idir = 1, dim
            write(str_tmp, '(f12.6)') momentum(idir, ist, ik+is)
            message(1) = trim(message(1)) // trim(str_tmp)
          end do
          write(str_tmp, '(3x,f12.6)') st%occ(ist, ik+is)
          message(1) = trim(message(1)) // trim(str_tmp)
          call messages_info(1, iunit)
 
        end do
      end do
 
      write(message(1),'(a)') ''
      call messages_info(1, iunit)
 
    end do
 
    safe_deallocate_a(momentum)
    call io_close(iunit)
 
    pop_sub(output_me_out_momentum)
  end subroutine output_me_out_momentum
 
 
  ! ---------------------------------------------------------
  subroutine output_me_out_ang_momentum(gr, st, space, fname, namespace, kpoints)
    type(grid_t),                 intent(in) :: gr
    type(states_elec_t),          intent(in) :: st
    type(space_t),                intent(in) :: space
    character(len=*),             intent(in) :: fname
    type(namespace_t),            intent(in) :: namespace
    type(kpoints_t),              intent(in) :: kpoints
 
    integer            :: iunit, ik, ist, is, ns, idir, k_start, k_end, k_n, dim, st_start, st_end, nst
    character(len=80)  :: tmp_str(space%dim), cspin
    real(real64)       :: angular(3), lsquare, kpoint(space%dim)
    real(real64), allocatable :: ang(:, :, :), ang2(:, :)
#if defined(HAVE_MPI)
    integer            :: tmp
#endif
    real(real64), allocatable :: lang(:, :)
 
    push_sub(output_me_out_ang_momentum)
 
    ns = 1
    if (st%d%nspin == 2) ns = 2
    assert(space%dim == 3)
    ! Define useful quantities
    dim = space%dim ! Space dimension
    st_start = st%st_start ! Starting state
    st_end = st%st_end ! Last state
    nst = st%nst ! Number of states
 
    iunit = io_open(fname, namespace, action='write')
 
    if (st%system_grp%is_root()) then
      write(iunit,'(a)') 'Warning: When non-local pseudopotentials are used '
      write(iunit,'(a)') '         the numbers below may be meaningless.    '
      write(iunit,'(a)') '                                                  '
      write(iunit,'(a)') 'Angular Momentum of the KS states [dimensionless]:'
      ! r x k is dimensionless. we do not include hbar.
      if (st%nik > ns) then
        message(1) = 'k-points [' // trim(units_abbrev(unit_one/units_out%length)) // ']'
        call messages_info(1, iunit)
      end if
    end if
 
    safe_allocate(ang(1:nst, 1:st%nik, 1:3))
    safe_allocate(ang2(1:nst, 1:st%nik))
 
    ! Compute matrix elements
    call elec_angular_momentum_me(gr, st, space, ang, ang2)
 
    k_start = st%d%kpt%start
    k_end = st%d%kpt%end
    do idir = 1, 3
      angular(idir) = states_elec_eigenvalues_sum(st, ang(st_start:st_end, k_start:k_end, idir))
    end do
    lsquare = states_elec_eigenvalues_sum(st, ang2(st_start:st_end, k_start:k_end))
 
    if (st%d%kpt%parallel) then
      k_n = st%d%kpt%nlocal()
 
      assert(.not. st%parallel_in_states)
 
      ! note: could use lmpi_gen_allgatherv here?
      safe_allocate(lang(1:st%lnst, 1:k_n))
      do idir = 1, 3
        lang(1:st%lnst, 1:k_n) = ang(st_start:st_end, k_start:k_end, idir)
        call st%d%kpt%mpi_grp%allgatherv(lang, nst*k_n, mpi_double_precision, ang(:, :, idir), &
          st%d%kpt%num(:)*nst, (st%d%kpt%range(1, :) - 1)*nst, mpi_double_precision)
      end do
      lang(1:st%lnst, 1:k_n) = ang2(st_start:st_end, k_start:k_end)
      call st%d%kpt%mpi_grp%allgatherv(lang, nst*k_n, mpi_double_precision, ang2, &
        st%d%kpt%num(:)*nst, (st%d%kpt%range(1, :) - 1)*nst, mpi_double_precision)
      safe_deallocate_a(lang)
    end if
 
    if (st%parallel_in_states) then
      safe_allocate(lang(1:st%lnst, 1))
    end if
 
    do ik = 1, st%nik, ns
      if (st%nik > ns) then
 
        kpoint(:) = kpoints%get_point(st%d%get_kpoint_index(ik))
 
        write(message(1), '(a,i4, a)') '#k =', ik, ', k = ('
        do idir = 1, dim
          write(tmp_str(1), '(f12.6, a)') units_from_atomic(unit_one/units_out%length, kpoint(idir)), ','
          message(1) = trim(message(1)) // trim(tmp_str(1))
          if (idir == dim) then
            message(1) = trim(message(1)) // ")"
          else
            message(1) = trim(message(1)) // ","
          end if
        end do
        call messages_info(1, iunit)
      end if
 
      ! Exchange ang and ang2.
      if (st%parallel_in_states) then
        assert(.not. st%d%kpt%parallel)
#if defined(HAVE_MPI)
        do is = 1, ns
          do idir = 1, 3
            lang(1:st%lnst, 1) = ang(st_start:st_end, ik+is-1, idir)
            call lmpi_gen_allgatherv(st%lnst, lang(:, 1), tmp, ang(:, ik+is-1, idir), st%mpi_grp)
          end do
          lang(1:st%lnst, 1) = ang2(st_start:st_end, ik+is-1)
          call lmpi_gen_allgatherv(st%lnst, lang(:, 1), tmp, ang2(:, ik+is-1), st%mpi_grp)
        end do
#endif
      end if
 
      write(message(1), '(a4,1x,a5,4a12,4x,a12,1x)')       &
        '#st',' Spin','        <Lx>', '        <Ly>', '        <Lz>', '        <L2>', 'Occupation '
      call messages_info(1, iunit)
 
      if (st%system_grp%is_root()) then
        do ist = 1, 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)
            write(tmp_str(2), '(1x,4f12.6,3x,f12.6)') &
              (ang(ist, ik+is, idir), idir = 1, 3), ang2(ist, ik+is), st%occ(ist, ik+is)
            message(1) = trim(tmp_str(1))//trim(tmp_str(2))
            call messages_info(1, iunit)
          end do
        end do
      end if
      write(message(1),'(a)') ''
      call messages_info(1, iunit)
 
    end do
 
    write(message(1),'(a)') 'Total Angular Momentum L [dimensionless]'
    write(message(2),'(10x,4f12.6)') angular(1:3), lsquare
    call messages_info(2, iunit)
 
    call io_close(iunit)
 
    if (st%parallel_in_states) then
      safe_deallocate_a(lang)
    end if
 
    safe_deallocate_a(ang)
    safe_deallocate_a(ang2)
 
    pop_sub(output_me_out_ang_momentum)
  end subroutine output_me_out_ang_momentum
 
  subroutine output_me_out_dipole(gr, st, space, dir, namespace, hm, ions, st_start, st_end)
    type(grid_t),                 intent(in) :: gr
    type(states_elec_t),          intent(in) :: st
    type(space_t),                intent(in) :: space
    character(len=*),             intent(in) :: dir
    type(namespace_t),            intent(in) :: namespace
    type(hamiltonian_elec_t),     intent(in) :: hm
    type(ions_t),                 intent(in) :: ions
    integer,                      intent(in) :: st_start, st_end
 
    integer :: ik, idir, ist, jst, iunit
    character(len=256) :: fname
    real(real64), allocatable :: ddip_elements(:,:,:)
    complex(real64), allocatable :: zdip_elements(:,:,:)
 
    push_sub(output_me_out_dipole)
 
    ! The content of each file should be clear from the header of each file.
    do ik = st%d%kpt%start, st%d%kpt%end
      write(fname,'(i8)') ik
      write(fname,'(a)') trim(dir)//'/ks_me_dipole.k'//trim(adjustl(fname))//'_'
 
      if (states_are_real(st)) then
        ! Allocate dipole array
        safe_allocate(ddip_elements(1:space%dim, st_start:st_end, st_start:st_end))
        ! Compute dipole elements
        call dipole_me(gr, st, namespace, hm, ions, ik, st_start, st_end, ddip_elements)
        ! Output elements
        do idir = 1, space%dim
          iunit = io_open(trim(fname)//index2axis(idir), namespace, action = 'write')
          write(iunit, '(a)') '# Dipole matrix elements file: |<Phi_i | r | Phi_j>|'
          write(iunit, '(a,i8)')      '# ik =', ik
          write(iunit, '(a)')    '# Units = ['//trim(units_abbrev(units_out%length))//']'
 
          do ist = st_start, st_end
            do jst = st_start, st_end
              write(iunit, fmt='(f20.12)', advance = 'no') units_from_atomic(units_out%length, abs(ddip_elements(idir, ist, jst)))
              write(iunit, fmt='(a)', advance = 'no') '  '
            end do
            write(iunit, '(a)') '' ! New line
          end do
          call io_close(iunit)
        end do
        safe_deallocate_a(ddip_elements)
      else
        ! Allocate dipole array
        safe_allocate(zdip_elements(1:space%dim, st_start:st_end, st_start:st_end))
        ! Compute dipole elements
        call dipole_me(gr, st, namespace, hm, ions, ik, st_start, st_end, zdip_elements)
        ! Output elements
        do idir = 1, space%dim
          iunit = io_open(trim(fname)//index2axis(idir), namespace, action = 'write')
          write(iunit, '(a)') '# Dipole matrix elements file: |<Phi_i | r | Phi_j>|'
          write(iunit, '(a,i8)')      '# ik =', ik
          write(iunit, '(a)')    '# Units = ['//trim(units_abbrev(units_out%length))//']'
 
          do ist = st_start, st_end
            do jst = st_start, st_end
              write(iunit, fmt='(f20.12)', advance = 'no') units_from_atomic(units_out%length, abs(zdip_elements(idir, ist, jst)))
              write(iunit, fmt='(a)', advance = 'no') '  '
            end do
            write(iunit, '(a)') '' ! New line
          end do
          call io_close(iunit)
        end do
        safe_deallocate_a(zdip_elements)
      end if
    end do
 
    pop_sub(output_me_out_dipole)
  end subroutine output_me_out_dipole
 
  subroutine output_me_out_ks_multipoles(gr, st, space, dir, this, namespace)
    type(grid_t),                 intent(in) :: gr
    type(states_elec_t),          intent(in) :: st
    type(space_t),                intent(in) :: space
    character(len=*),             intent(in) :: dir
    type(output_me_t),            intent(in) :: this
    type(namespace_t),            intent(in) :: namespace
 
    integer :: id, ll, mm, ik, iunit
    character(len=256) :: fname
    real(real64), allocatable :: delements(:,:)
    complex(real64), allocatable :: zelements(:,:)
 
    push_sub(output_me_out_ks_multipoles)
 
    ! The content of each file should be clear from the header of each file.
    id = 1
    mm = 0
    do ik = 1, st%nik
      select case (space%dim)
      case (3)
        do ll = 1, this%ks_multipoles
          do mm = -ll, ll
            ! Define file name
            write(fname,'(i4)') id
            write(fname,'(a)') trim(dir) // '/ks_me_multipoles.' // trim(adjustl(fname))
            iunit = io_open(fname, namespace, action = 'write')
            ! Print header
            call print_ks_multipole_header(iunit, ll, mm, ik)
            ! Compute and print elements
            if (states_are_real(st)) then
              safe_allocate(delements(1:st%nst, 1:st%nst))
              call ks_multipoles_3d(gr, st, space, ll, mm, ik, delements)
              call dprint_ks_multipole(iunit, st%nst, delements, ll)
              safe_deallocate_a(delements)
            else
              safe_allocate(zelements(1:st%nst, 1:st%nst))
              call ks_multipoles_3d(gr, st, space, ll, mm, ik, zelements)
              call zprint_ks_multipole(iunit, st%nst, zelements, ll)
              safe_deallocate_a(zelements)
            end if
            ! Close file
            call io_close(iunit)
            id = id + 1
          end do
        end do
      case (2)
        do ll = 1, 2
          mm = 0
          ! Define file name
          write(fname,'(i4)') id
          write(fname,'(a)') trim(dir) // '/ks_me_multipoles.' // trim(adjustl(fname))
          iunit = io_open(fname, namespace, action = 'write')
          ! Print header
          call print_ks_multipole_header(iunit, ll, mm, ik)
          ! Compute and print elements
          if (states_are_real(st)) then
            safe_allocate(delements(1:st%nst, 1:st%nst))
            call ks_multipoles_2d(gr, st, ll, ik, delements)
            call dprint_ks_multipole(iunit, st%nst, delements)
            safe_deallocate_a(delements)
          else
            safe_allocate(zelements(1:st%nst, 1:st%nst))
            call ks_multipoles_2d(gr, st, ll, ik, zelements)
            call zprint_ks_multipole(iunit, st%nst, zelements)
            safe_deallocate_a(zelements)
          end if
          ! Close file
          call io_close(iunit)
          id = id + 1
        end do
      case (1)
        do ll = 1, this%ks_multipoles
          mm = 0
          ! Define file name
          write(fname,'(i4)') id
          write(fname,'(a)') trim(dir) // '/ks_me_multipoles.' // trim(adjustl(fname))
          iunit = io_open(fname, namespace, action = 'write')
          ! Print header
          call print_ks_multipole_header(iunit, ll, mm, ik)
          ! Compute and print elements
          if (states_are_real(st)) then
            safe_allocate(delements(1:st%nst, 1:st%nst))
            call ks_multipoles_1d(gr, st, ll, ik, delements)
            call dprint_ks_multipole(iunit, st%nst, delements, ll)
            safe_deallocate_a(delements)
          else
            safe_allocate(zelements(1:st%nst, 1:st%nst))
            call ks_multipoles_1d(gr, st, ll, ik, zelements)
            call zprint_ks_multipole(iunit, st%nst, zelements, ll)
            safe_deallocate_a(zelements)
          end if
          ! Close file
          call io_close(iunit)
          id = id + 1
        end do
      end select
    end do
 
    pop_sub(output_me_out_ks_multipoles)
 
  contains
 
    subroutine print_ks_multipole_header(iunit, ll, mm, ik)
      integer, intent(in) :: iunit, ll, mm, ik
 
      write(iunit, fmt = '(a)') '# Multipole matrix elements file: <Phi_i | r**l * Y_{lm}(theta,phi) | Phi_j>'
      write(iunit, fmt = '(a,i2,a,i2)') '# l =', ll, '; m =', mm
      write(iunit, fmt = '(a,i8)')      '# ik =', ik
      if (ll>1) then
        write(iunit, fmt = '(a,i1)') '# Units = ['//trim(units_abbrev(units_out%length))//']^', ll
      else
        write(iunit, fmt = '(a)')    '# Units = ['//trim(units_abbrev(units_out%length))//']'
      end if
    end subroutine print_ks_multipole_header
 
    subroutine dprint_ks_multipole(iunit, nst, elements, ll)
      integer, intent(in) :: iunit, nst
      real(real64), contiguous, intent(in) :: elements(:,:)
      integer, optional, intent(in) :: ll
 
      integer :: ist, jst
      real(real64)   :: element
 
      do ist = 1, nst
        do jst = 1, nst
          element = units_from_atomic(units_out%length**optional_default(ll, 1), elements(ist, jst))
          write(iunit, fmt='(f20.12, a)', advance = 'no') element, '   '
        end do
        write(iunit, '(a)') ''
      end do
    end subroutine dprint_ks_multipole
 
    subroutine zprint_ks_multipole(iunit, nst, elements, ll)
      integer, intent(in) :: iunit, nst
      complex(real64), contiguous, intent(in) :: elements(:,:)
      integer, optional, intent(in) :: ll
 
      integer :: ist, jst
      complex(real64)   :: element
 
      do ist = 1, nst
        do jst = 1, nst
          element = units_from_atomic(units_out%length**optional_default(ll, 1), elements(ist, jst))
          write(iunit, fmt='(f20.12,a1,f20.12,a)', advance = 'no') real(element),',',aimag(element), '   '
        end do
        write(iunit, '(a)') ''
      end do
    end subroutine zprint_ks_multipole
 
  end subroutine output_me_out_ks_multipoles
 
 
#include "undef.F90"
#include "real.F90"
 
#include "undef.F90"
#include "complex.F90"
 
end module output_me_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: