output.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 output_oct_m
  use accel_oct_m
  use basins_oct_m
  use box_oct_m
  use comm_oct_m
  use cube_oct_m
  use cube_function_oct_m
  use current_oct_m
  use debug_oct_m
  use density_oct_m
  use derivatives_oct_m
  use dos_oct_m
  use electron_space_oct_m
  use elf_oct_m
#if defined(HAVE_ETSF_IO)
  use etsf_io
  use etsf_io_tools
#endif
  use exchange_operator_oct_m
  use ext_partner_list_oct_m
  use fft_oct_m
  use fourier_shell_oct_m
  use fourier_space_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_oct_m
  use interaction_partner_oct_m
  use io_oct_m
  use io_function_oct_m
  use ions_oct_m
  use, intrinsic :: iso_fortran_env
  use kpoints_oct_m
  use ks_potential_oct_m
  use lasers_oct_m
  use lda_u_oct_m
  use lda_u_io_oct_m
  use linear_response_oct_m
  use loct_oct_m
  use math_oct_m
  use magnetic_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 output_low_oct_m
  use output_me_oct_m
  use output_berkeleygw_oct_m
  use parser_oct_m
  use poisson_oct_m
  use profiling_oct_m
  use smear_oct_m
  use space_oct_m
  use species_oct_m
  use states_abst_oct_m
  use states_elec_oct_m
  use states_elec_io_oct_m
  use stress_oct_m
  use string_oct_m
  use submesh_oct_m
  use symm_op_oct_m
  use symmetries_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use utils_oct_m
  use varinfo_oct_m
  use v_ks_oct_m
  use xc_oct_m
  use xc_oep_oct_m
  use xc_oep_photon_oct_m
  use xc_f03_lib_m
  use xc_vxc_oct_m
 
  implicit none
 
  private
  public ::              &
    output_bgw_t,        &
    output_init,         &
    output_states,       &
    output_hamiltonian,  &
    output_all,          &
    output_current_flow, &
    doutput_lr,          &
    zoutput_lr,          &
    output_scalar_pot,   &
    output_needs_current, &
    output_need_exchange
 
contains
 
  subroutine output_init(outp, namespace, space, st, gr, nst, ks)
    type(output_t),            intent(out)   :: outp
    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
    type(v_ks_t),              intent(inout) :: ks
 
    type(block_t) :: blk
    real(real64) :: norm
    character(len=80) :: nst_string, default
#ifndef HAVE_ETSF_IO
    integer :: iout
#endif
 
    push_sub(output_init)
    outp%what = .false.
 
    call io_function_read_what_how_when(namespace, space, outp%what, outp%how, outp%output_interval)
 
    if (outp%what(option__output__wfs_fourier)) then
      if (accel_is_enabled()) then
        message(1) = "Wave functions in Fourier space not supported on GPUs."
        call messages_fatal(1, namespace=namespace)
      end if
      call messages_experimental("Wave-functions in Fourier space", namespace=namespace)
    end if
 
    ! cannot calculate the ELF in 1D
    if (outp%what(option__output__elf) .or. outp%what(option__output__elf_basins)) then
      if (space%dim /= 2 .and. space%dim /= 3) then
        outp%what(option__output__elf) = .false.
        outp%what(option__output__elf_basins) = .false.
        write(message(1), '(a)') 'Cannot calculate ELF except in 2D and 3D.'
        call messages_warning(1, namespace=namespace)
      end if
    end if
 
 
    if (outp%what(option__output__mmb_wfs)) then
      call messages_experimental("Model many-body wfs", namespace=namespace)
    end if
 
    if (outp%what(option__output__xc_torque)) then
      if (st%d%ispin /= spinors) then
        write(message(1), '(a)') 'The output xc_torque can only be computed for spinors.'
        call messages_fatal(1, namespace=namespace)
      end if
      if (space%dim /= 3) then
        write(message(1), '(a)') 'The output xc_torque can only be computed in the 3D case.'
        call messages_fatal(1, namespace=namespace)
      end if
    end if
    if (outp%what(option__output__mmb_den)) then
      call messages_experimental("Model many-body density matrix", namespace=namespace)
      ! NOTES:
      !   could be made into block to be able to specify which dimensions to trace
      !   in principle all combinations are interesting, but this means we need to
      !   be able to output density matrices for multiple particles or multiple
      !   dimensions. The current 1D 1-particle case is simple.
    end if
 
    if (outp%what(option__output__energy_density)) call messages_experimental("'Output = energy_density'", namespace=namespace)
    if (outp%what(option__output__heat_current)) call messages_experimental("'Output = heat_current'", namespace=namespace)
 
    if (outp%what(option__output__wfs) .or. outp%what(option__output__wfs_sqmod)) then
 
      !%Variable OutputWfsNumber
      !%Type string
      !%Default all states
      !%Section Output
      !%Description
      !% Which wavefunctions to print, in list form: <i>i.e.</i>, "1-5" to print the first
      !% five states, "2,3" to print the second and the third state, etc.
      !% If more states are specified than available, extra ones will be ignored.
      !%End
 
      write(nst_string,'(i6)') nst
      write(default,'(a,a)') "1-", trim(adjustl(nst_string))
      call parse_variable(namespace, 'OutputWfsNumber', default, outp%wfs_list)
    end if
 
    if (parse_block(namespace, 'CurrentThroughPlane', blk) == 0) then
      if (.not. outp%what(option__output__j_flow)) then
        outp%what(option__output__j_flow) = .true.
        call parse_variable(namespace, 'OutputInterval', 50, outp%output_interval(option__output__j_flow))
      end if
 
      !%Variable CurrentThroughPlane
      !%Type block
      !%Section Output
      !%Description
      !% The code can calculate current
      !% traversing a user-defined portion of a plane, as specified by this block.
      !% A small plain-text file <tt>current-flow</tt> will be written containing this information.
      !% Only available for 1D, 2D, or 3D.
      !% In the format below, <tt>origin</tt> is a point in the plane.
      !% <tt>u</tt> and <tt>v</tt> are the (dimensionless) vectors defining the plane;
      !% they will be normalized. <tt>spacing</tt> is the fineness of the mesh
      !% on the plane. Integers <tt>nu</tt> and <tt>mu</tt> are the length and
      !% width of the portion of the plane, in units of <tt>spacing</tt>.
      !% Thus, the grid points included in the plane are
      !% <tt>x_ij = origin + i*spacing*u + j*spacing*v</tt>,
      !% for <tt>nu <= i <= mu</tt> and <tt>nv <= j <= mv</tt>.
      !% Analogously, in the 2D case, the current flow is calculated through a line;
      !% in the 1D case, the current flow is calculated through a point. Note that the spacing
      !% can differ from the one used in the main calculation; an interpolation will be performed.
      !%
      !% Example (3D):
      !%
      !% <tt>%CurrentThroughPlane
      !% <br>&nbsp;&nbsp; 0.0 | 0.0 | 0.0  # origin
      !% <br>&nbsp;&nbsp; 0.0 | 1.0 | 0.0  # u
      !% <br>&nbsp;&nbsp; 0.0 | 0.0 | 1.0  # v
      !% <br>&nbsp;&nbsp; 0.2              # spacing
      !% <br>&nbsp;&nbsp; 0 | 50           # nu | mu
      !% <br>&nbsp;&nbsp; -50 | 50         # nv | mv
      !% <br>%</tt>
      !%
      !% Example (2D):
      !%
      !% <tt>%CurrentThroughPlane
      !% <br>&nbsp;&nbsp; 0.0 | 0.0        # origin
      !% <br>&nbsp;&nbsp; 1.0 | 0.0        # u
      !% <br>&nbsp;&nbsp; 0.2              # spacing
      !% <br>&nbsp;&nbsp; 0 | 50           # nu | mu
      !% <br>%</tt>
      !%
      !% Example (1D):
      !%
      !% <tt>%CurrentThroughPlane
      !% <br>&nbsp;&nbsp; 0.0              # origin
      !% <br>%</tt>
      !%
      !%End
 
      select case (space%dim)
      case (3)
 
        call parse_block_float(blk, 0, 0, outp%plane%origin(1), units_inp%length)
        call parse_block_float(blk, 0, 1, outp%plane%origin(2), units_inp%length)
        call parse_block_float(blk, 0, 2, outp%plane%origin(3), units_inp%length)
        call parse_block_float(blk, 1, 0, outp%plane%u(1))
        call parse_block_float(blk, 1, 1, outp%plane%u(2))
        call parse_block_float(blk, 1, 2, outp%plane%u(3))
        call parse_block_float(blk, 2, 0, outp%plane%v(1))
        call parse_block_float(blk, 2, 1, outp%plane%v(2))
        call parse_block_float(blk, 2, 2, outp%plane%v(3))
        call parse_block_float(blk, 3, 0, outp%plane%spacing, units_inp%length)
        call parse_block_integer(blk, 4, 0, outp%plane%nu)
        call parse_block_integer(blk, 4, 1, outp%plane%mu)
        call parse_block_integer(blk, 5, 0, outp%plane%nv)
        call parse_block_integer(blk, 5, 1, outp%plane%mv)
 
        norm = norm2(outp%plane%u(1:3))
        if (norm < m_epsilon) then
          write(message(1), '(a)') 'u-vector for CurrentThroughPlane cannot have norm zero.'
          call messages_fatal(1, namespace=namespace)
        end if
        outp%plane%u(1:3) = outp%plane%u(1:3) / norm
 
        norm = norm2(outp%plane%v(1:3))
        if (norm < m_epsilon) then
          write(message(1), '(a)') 'v-vector for CurrentThroughPlane cannot have norm zero.'
          call messages_fatal(1, namespace=namespace)
        end if
        outp%plane%v(1:3) = outp%plane%v(1:3) / norm
 
        outp%plane%n(1:3) = dcross_product(outp%plane%u(1:3), outp%plane%v(1:3))
 
      case (2)
 
        call parse_block_float(blk, 0, 0, outp%line%origin(1), units_inp%length)
        call parse_block_float(blk, 0, 1, outp%line%origin(2), units_inp%length)
        call parse_block_float(blk, 1, 0, outp%line%u(1))
        call parse_block_float(blk, 1, 1, outp%line%u(2))
        call parse_block_float(blk, 2, 0, outp%line%spacing, units_inp%length)
        call parse_block_integer(blk, 3, 0, outp%line%nu)
        call parse_block_integer(blk, 3, 1, outp%line%mu)
 
        norm = norm2(outp%line%u(1:2))
        if (norm < m_epsilon) then
          write(message(1), '(a)') 'u-vector for CurrentThroughPlane cannot have norm zero.'
          call messages_fatal(1, namespace=namespace)
        end if
        outp%line%u(1:2) = outp%line%u(1:2) / norm
 
        outp%line%n(1) = -outp%line%u(2)
        outp%line%n(2) =  outp%line%u(1)
 
      case (1)
 
        call parse_block_float(blk, 0, 0, outp%line%origin(1), units_inp%length)
 
      case default
 
        call messages_not_implemented("CurrentThroughPlane for 4D or higher", namespace=namespace)
 
      end select
      call parse_block_end(blk)
    end if
 
    if (outp%what(option__output__matrix_elements)) then
      call output_me_init(outp%me, namespace, space, st, gr, nst)
    else
      outp%me%what = .false.
    end if
 
    if (outp%what(option__output__berkeleygw)) then
      if (accel_is_enabled()) then
        message(1) = "BerkeleyGW is not compatible with GPUs."
        call messages_fatal(1, namespace=namespace)
      end if
      call output_berkeleygw_init(nst, namespace, outp%bgw, space%periodic_dim)
    end if
 
    ! required for output_hamiltonian()
    if (outp%what(option__output__potential_gradient) .and. .not. outp%what(option__output__potential)) then
      outp%what(option__output__potential) = .true.
      outp%output_interval(option__output__potential) = outp%output_interval(option__output__potential_gradient)
    end if
 
 
    !%Variable OutputDuringSCF
    !%Type logical
    !%Default no
    !%Section Output
    !%Description
    !% During <tt>gs</tt> and <tt>unocc</tt> runs, if this variable is set to yes,
    !% output will be written after every <tt>OutputInterval</tt> iterations.
    !%End
    call parse_variable(namespace, 'OutputDuringSCF', .false., outp%duringscf)
 
    if (parse_is_defined(namespace, 'RestartWriteInterval')) then
      write(message(1), '(a)') 'Input variable RestartWriteInterval is obsolete.'
      write(message(2), '(a)') 'Restart files are now written in periods of the wallclock time'
      write(message(3), '(a)') 'given by RestartWallTimePeriod, so you can simply delete this variable.'
      call messages_fatal(3, only_root_writes=.true., namespace=namespace)
    end if
 
    !%Variable OutputIterDir
    !%Default "output_iter"
    !%Type string
    !%Section Output
    !%Description
    !% The name of the directory where <tt>Octopus</tt> stores information
    !% such as the density, forces, etc. requested by variable <tt>Output</tt>
    !% in the format specified by <tt>OutputFormat</tt>.
    !% This information is written while iterating <tt>CalculationMode = gs</tt>, <tt>unocc</tt>, or <tt>td</tt>,
    !% according to <tt>OutputInterval</tt>, and has nothing to do with the restart information.
    !%End
    call parse_variable(namespace, 'OutputIterDir', "output_iter", outp%iter_dir)
    if (any(outp%what) .and. any(outp%output_interval > 0)) then
      call io_mkdir(outp%iter_dir, namespace)
    end if
    call add_last_slash(outp%iter_dir)
 
    ! At this point, we don`t know whether the states will be real or complex.
    ! We therefore pass .false. to states_are_real, and need to check for real states later.
 
    if (output_needs_current(outp, .false.)) then
      call v_ks_calculate_current(ks, .true.)
    else
      call v_ks_calculate_current(ks, .false.)
    end if
 
    if (outp%what(option__output__current_dia)) then
      message(1) = "The diamagnetic current will be calculated only if CalculateDiamagneticCurrent = yes."
      call messages_warning(1, namespace=namespace)
    end if
 
#ifndef HAVE_ETSF_IO
    do iout = lbound(outp%how, 1), ubound(outp%how, 1)
      if (bitand(outp%how(iout), option__outputformat__etsf) /= 0) then
        message(1) = "ETSF output can only be used if Octopus is compiled with ETSF-IO"
        call messages_fatal(1, namespace=namespace)
      end if
    end do
#endif
 
    pop_sub(output_init)
  end subroutine output_init
 
  ! ---------------------------------------------------------
  subroutine output_all(outp, namespace, space, dir, gr, ions, iter, st, hm, ks)
    type(output_t),           intent(in)    :: outp
    type(namespace_t),        intent(in)    :: namespace
    class(space_t),           intent(in)    :: space
    character(len=*),         intent(in)    :: dir
    type(grid_t),             intent(in)    :: gr
    type(ions_t),             intent(in)    :: ions
    integer,                  intent(in)    :: iter
    type(states_elec_t),      intent(inout) :: st
    type(hamiltonian_elec_t), intent(inout) :: hm
    type(v_ks_t),             intent(inout) :: ks
 
    integer :: idir, ierr, iout, iunit
    character(len=MAX_PATH_LEN) :: fname
 
    push_sub(output_all)
    call profiling_in("OUTPUT_ALL")
 
    if (any(outp%what)) then
      message(1) = "Info: Writing output to " // trim(dir)
      call messages_info(1, namespace=namespace)
      call io_mkdir(dir, namespace)
    end if
 
    if (outp%what_now(option__output__mesh_r, iter)) then
      do idir = 1, space%dim
        write(fname, '(a,a)') 'mesh_r-', index2axis(idir)
        call dio_function_output(outp%how(option__output__mesh_r), dir, fname, namespace, space, &
          gr, gr%x_t(:,idir), units_out%length, ierr, pos=ions%pos, atoms=ions%atom)
      end do
    end if
 
    call output_states(outp, namespace, space, dir, st, gr, ions, hm, iter)
    call output_hamiltonian(outp, namespace, space, dir, hm, st, gr%der, ions, gr, iter, st%st_kpt_mpi_grp)
 
    ! We can only test against the theory level here, as it is not set in the call of output_init().
    if (outp%what_now(option__output__el_pressure, iter)) then
      if(ks%theory_level /= kohn_sham_dft) then
        call messages_not_implemented("el_pressure for TheoryLevel different from kohn_sham", namespace=namespace)
      end if
      if (st%d%spin_channels > 1) then
        call messages_not_implemented("el_pressure for spin-polarized or spinors", namespace=namespace)
      end if
    end if
 
    !hm not initialized yet when calling output_init()
    if (outp%what(option__output__kanamoriu) .and. hm%lda_u_level /= dft_u_acbn0) then
      message(1) = "kanamoriU output can only be computed for DFTULevel = dft_u_acbn0"
      call messages_fatal(1, namespace=namespace)
    end if
 
 
    call output_localization_funct(outp, namespace, space, dir, st, hm, gr, ions, iter)
 
    if (outp%what_now(option__output__j_flow, iter)) then
      call output_current_flow(outp, namespace, space, dir, gr, st, hm%kpoints)
    end if
 
    if (outp%what_now(option__output__geometry, iter)) then
      if (bitand(outp%how(option__output__geometry), option__outputformat__xcrysden) /= 0) then
        call write_xsf_geometry_file(dir, "geometry", space, ions%latt, ions%pos, ions%atom, gr, namespace)
      end if
      if (bitand(outp%how(option__output__geometry), option__outputformat__xyz) /= 0) then
        call ions%write_xyz(trim(dir)//'/geometry')
        if (ions%space%is_periodic()) then
          call ions%write_crystal(dir)
        end if
      end if
      if (bitand(outp%how(option__output__geometry), option__outputformat__vtk) /= 0) then
        call ions%write_vtk_geometry(trim(dir)//'/geometry')
      end if
      if (bitand(outp%how(option__output__geometry), option__outputformat__poscar) /= 0) then
        call ions%write_poscar(trim(dir)//'/POSCAR')
      end if
    end if
 
    if (outp%what_now(option__output__forces, iter)) then
      if (bitand(outp%how(option__output__forces), option__outputformat__bild) /= 0) then
        call ions%write_bild_forces_file(dir, "forces")
      else
        call write_xsf_geometry_file(dir, "forces", space, ions%latt, ions%pos, ions%atom, &
          gr, namespace, total_forces = ions%tot_force)
      end if
    end if
 
    if (outp%what_now(option__output__matrix_elements, iter)) then
      call output_me(outp%me, namespace, space, dir, st, gr, ions, hm)
    end if
 
    do iout = lbound(outp%how, 1), ubound(outp%how, 1)
      if (bitand(outp%how(iout), option__outputformat__etsf) /= 0) then
        call output_etsf(outp, namespace, space, dir, st, gr, hm%kpoints, ions, iter)
        exit
      end if
    end do
 
    if (outp%what_now(option__output__berkeleygw, iter)) then
      call output_berkeleygw(outp%bgw, namespace, space, dir, st, gr, ks, hm, ions)
    end if
 
    if (outp%what_now(option__output__energy_density, iter)) then
      call output_energy_density(outp, namespace, space, dir, hm, ks, st, ions, gr)
    end if
 
    if (outp%what_now(option__output__stress, iter)) then
      call io_mkdir(dir, namespace)
      iunit = io_open(trim(dir)//'/stress', namespace, action='write')
      call output_stress(iunit, space%periodic_dim, st%stress_tensors)
      call io_close(iunit)
    end if
 
    if (hm%lda_u_level /= dft_u_none) then
      if (outp%what_now(option__output__occ_matrices, iter))&
        call lda_u_write_occupation_matrices(dir, hm%lda_u, st, namespace)
 
      if (outp%what_now(option__output__effectiveu, iter))&
        call lda_u_write_effectiveu(dir, hm%lda_u, st, namespace)
 
      if (outp%what_now(option__output__magnetization, iter))&
        call lda_u_write_magnetization(dir, hm%lda_u, ions, gr, st, namespace)
 
      if (outp%what_now(option__output__local_orbitals, iter))&
        call output_dftu_orbitals(outp, dir, namespace, space, hm%lda_u, st, gr, ions, hm%phase%is_allocated())
 
      if (outp%what_now(option__output__kanamoriu, iter))&
        call lda_u_write_kanamoriu(dir, st, hm%lda_u, namespace)
 
      if (ks%oep_photon%level == oep_level_full) then
        if (outp%what_now(option__output__photon_correlator, iter)) then
          write(fname, '(a)') 'photon_correlator'
          call dio_function_output(outp%how(option__output__photon_correlator), dir, trim(fname), namespace, space, &
            gr, ks%oep_photon%pt%correlator(:,1), units_out%length, ierr, pos=ions%pos, atoms=ions%atom)
        end if
      end if
    end if
 
    ! We can only test against the theory level here, as it is not set in the call of output_init().
    if (outp%what_now(option__output__xc_torque, iter)) then
      if (ks%theory_level /= kohn_sham_dft .and. ks%theory_level /= generalized_kohn_sham_dft) then
        write(message(1), '(a)') 'The output xc_torque can only be computed when there is a xc potential.'
        call messages_fatal(1, namespace=namespace)
      end if
    end if
 
    call output_xc_torque(outp, namespace, dir, gr, hm, st, ions, ions%space)
 
    call profiling_out("OUTPUT_ALL")
    pop_sub(output_all)
  end subroutine output_all
 
 
  ! ---------------------------------------------------------
  subroutine output_localization_funct(outp, namespace, space, dir, st, hm, gr, ions, iter)
    type(output_t),           intent(in)    :: outp
    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(hamiltonian_elec_t), intent(in)    :: hm
    type(grid_t),             intent(in)    :: gr
    type(ions_t),             intent(in)    :: ions
    integer,                  intent(in)    :: iter
 
    real(real64), allocatable :: f_loc(:,:)
    character(len=MAX_PATH_LEN) :: fname
    integer :: is, ierr, imax
    type(mpi_grp_t) :: mpi_grp
 
    push_sub(output_localization_funct)
 
    mpi_grp = st%dom_st_kpt_mpi_grp
 
    ! if SPIN_POLARIZED, the ELF contains one extra channel: the total ELF
    imax = st%d%nspin
    if (st%d%ispin == spin_polarized) imax = 3
 
    safe_allocate(f_loc(1:gr%np, 1:imax))
 
    ! First the ELF in real space
    if (outp%what_now(option__output__elf, iter) .or. outp%what_now(option__output__elf_basins, iter)) then
      assert(space%dim /= 1)
 
      call elf_calc(space, st, gr, hm%kpoints, f_loc)
 
      ! output ELF in real space
      if (outp%what_now(option__output__elf, iter)) then
        write(fname, '(a)') 'elf_rs'
        call dio_function_output(outp%how(option__output__elf), dir, trim(fname), namespace, space, gr, &
          f_loc(:,imax), unit_one, ierr, pos=ions%pos, atoms=ions%atom, grp = mpi_grp)
        ! this quantity is dimensionless
 
        if (st%d%ispin /= unpolarized) then
          do is = 1, 2
            write(fname, '(a,i1)') 'elf_rs-sp', is
            call dio_function_output(outp%how(option__output__elf), dir, trim(fname), namespace, space, gr, &
              f_loc(:, is), unit_one, ierr, pos=ions%pos, atoms=ions%atom, grp = mpi_grp)
            ! this quantity is dimensionless
          end do
        end if
      end if
 
      if (outp%what_now(option__output__elf_basins, iter)) then
        call out_basins(f_loc(:,1), "elf_rs_basins", outp%how(option__output__elf_basins))
      end if
    end if
 
    ! Now Bader analysis
    if (outp%what_now(option__output__bader, iter)) then
      do is = 1, st%d%nspin
        call dderivatives_lapl(gr%der, st%rho(:,is), f_loc(:,is))
 
        fname = get_filename_with_spin('bader', st%d%nspin, is)
 
        call dio_function_output(outp%how(option__output__bader), dir, trim(fname), namespace, space, gr, &
          f_loc(:,is), units_out%length**(-2 - space%dim), ierr, &
          pos=ions%pos, atoms=ions%atom, grp = mpi_grp)
 
        fname = get_filename_with_spin('bader_basins', st%d%nspin, is)
        call out_basins(f_loc(:, is), fname, outp%how(option__output__bader))
      end do
    end if
 
    ! Now the pressure
    if (outp%what_now(option__output__el_pressure, iter)) then
      call calc_electronic_pressure(st, hm, gr, f_loc(:,1))
      call dio_function_output(outp%how(option__output__el_pressure), dir, "el_pressure", namespace, space, gr, &
        f_loc(:,1), unit_one, ierr, pos=ions%pos, atoms=ions%atom, grp = mpi_grp)
      ! this quantity is dimensionless
    end if
 
    safe_deallocate_a(f_loc)
 
    pop_sub(output_localization_funct)
 
  contains
    ! ---------------------------------------------------------
    subroutine out_basins(ff, filename, output_how)
      real(real64),     intent(in)    :: ff(:)
      character(len=*), intent(in)    :: filename
      integer(int64),      intent(in)    :: output_how
 
      character(len=MAX_PATH_LEN) :: fname
      type(basins_t)     :: basins
      integer            :: iunit
 
      push_sub(output_localization_funct.out_basins)
 
      call basins_init(basins, namespace, gr)
      call basins_analyze(basins, namespace, gr, ff(:), st%rho, 0.01_real64)
 
      call dio_function_output(output_how, dir, trim(filename), namespace, space, gr, &
        real(basins%map, real64) , unit_one, ierr, pos=ions%pos, atoms=ions%atom, grp = mpi_grp)
      ! this quantity is dimensionless
 
      write(fname,'(4a)') trim(dir), '/', trim(filename), '.info'
      iunit = io_open(trim(fname), namespace, action = 'write')
      call basins_write(basins, gr, iunit)
      call io_close(iunit)
 
      call basins_end(basins)
 
      pop_sub(output_localization_funct.out_basins)
    end subroutine out_basins
 
  end subroutine output_localization_funct
 
 
  ! ---------------------------------------------------------
  subroutine calc_electronic_pressure(st, hm, gr, pressure)
    type(states_elec_t),      intent(inout) :: st
    type(hamiltonian_elec_t), intent(in)    :: hm
    type(grid_t),             intent(in)    :: gr
    real(real64),             intent(out)   :: pressure(:)
 
    real(real64), allocatable :: rho(:,:), lrho(:), tau(:,:)
    real(real64)   :: p_tf, dens
    integer :: is, ii
 
    push_sub(calc_electronic_pressure)
 
    safe_allocate( rho(1:gr%np_part, 1:st%d%nspin))
    safe_allocate(lrho(1:gr%np))
    safe_allocate( tau(1:gr%np, 1:st%d%nspin))
 
    rho = m_zero
    call density_calc(st, gr, rho)
    call states_elec_calc_quantities(gr, st, hm%kpoints, .false., kinetic_energy_density = tau)
 
    pressure = m_zero
    do is = 1, st%d%spin_channels
      lrho = m_zero
      call dderivatives_lapl(gr%der, rho(:, is), lrho)
 
      pressure(:) = pressure(:) + &
        tau(:, is)/m_three - lrho(:)/m_four
    end do
 
    do ii = 1, gr%np
      dens = sum(rho(ii,1:st%d%spin_channels))
 
      p_tf = m_two/m_five*(m_three*m_pi**2)**(m_two/m_three)* &
        dens**(m_five/m_three)
 
      ! add XC pressure
      ! FIXME: Not correct for spinors and spin-polarized here
      pressure(ii) = pressure(ii) + (dens*hm%ks_pot%vxc(ii,1) - hm%energy%exchange - hm%energy%correlation)
 
      pressure(ii) = pressure(ii)/p_tf
      pressure(ii) = m_half*(m_one + pressure(ii)/sqrt(m_one + pressure(ii)**2))
    end do
 
    pop_sub(calc_electronic_pressure)
  end subroutine calc_electronic_pressure
 
 
  ! ---------------------------------------------------------
  subroutine output_energy_density(outp, namespace, space, dir, hm, ks, st, ions, gr)
    type(output_t),            intent(in) :: outp
    type(namespace_t),         intent(in) :: namespace
    class(space_t),            intent(in) :: space
    character(len=*),          intent(in) :: dir
    type(hamiltonian_elec_t),  intent(in) :: hm
    type(v_ks_t),              intent(inout) :: ks
    type(states_elec_t),       intent(in) :: st
    type(ions_t),              intent(in) :: ions
    type(grid_t),              intent(in) :: gr
 
    integer :: is, ierr, ip
    character(len=MAX_PATH_LEN) :: fname
    type(unit_t) :: fn_unit
    real(real64), allocatable :: energy_density(:, :)
    real(real64), allocatable :: ex_density(:)
    real(real64), allocatable :: ec_density(:)
 
    push_sub(output_energy_density)
 
    fn_unit = units_out%energy*units_out%length**(-space%dim)
    safe_allocate(energy_density(1:gr%np, 1:st%d%nspin))
 
    ! the kinetic energy density
    call states_elec_calc_quantities(gr, st, hm%kpoints, .true., kinetic_energy_density = energy_density)
 
    ! the external potential energy density
    do is = 1, st%d%nspin
      do ip = 1, gr%np
        energy_density(ip, is) = energy_density(ip, is) + st%rho(ip, is)*hm%ep%vpsl(ip)
      end do
    end do
 
    ! the hartree energy density
    do is = 1, st%d%nspin
      do ip = 1, gr%np
        energy_density(ip, is) = energy_density(ip, is) + m_half*st%rho(ip, is)*hm%ks_pot%vhartree(ip)
      end do
    end do
 
    ! the XC energy density
    safe_allocate(ex_density(1:gr%np))
    safe_allocate(ec_density(1:gr%np))
 
    call xc_get_vxc(gr, ks%xc, st, hm%kpoints, hm%psolver, namespace, space, st%rho, st%d%ispin, &
      hm%ions%latt%rcell_volume, ex_density = ex_density, ec_density = ec_density)
    do is = 1, st%d%nspin
      do ip = 1, gr%np
        energy_density(ip, is) = energy_density(ip, is) + ex_density(ip) + ec_density(ip)
      end do
    end do
 
    safe_deallocate_a(ex_density)
    safe_deallocate_a(ec_density)
 
    do is = 1, st%d%spin_channels
      fname = get_filename_with_spin('energy_density', st%d%nspin, is)
      call dio_function_output(outp%how(option__output__energy_density), dir, trim(fname), namespace, space, gr, &
        energy_density(:, is), unit_one, ierr, pos=ions%pos, atoms=ions%atom, grp = st%dom_st_kpt_mpi_grp)
    end do
    safe_deallocate_a(energy_density)
 
    pop_sub(output_energy_density)
  end subroutine output_energy_density
 
  !--------------------------------------------------------------
 
  logical function output_need_exchange(outp) result(need_exx)
    type(output_t),         intent(in)    :: outp
 
    need_exx =(outp%what(option__output__berkeleygw) &
      .or. outp%me%what(option__outputmatrixelements__two_body) &
      .or. outp%me%what(option__outputmatrixelements__two_body_exc_k))
  end function output_need_exchange
 
 
  ! ---------------------------------------------------------
  subroutine output_dftu_orbitals(outp, dir, namespace, space, this, st, mesh, ions, has_phase)
    type(output_t),      intent(in) :: outp
    character(len=*),    intent(in) :: dir
    type(namespace_t),   intent(in) :: namespace
    class(space_t),      intent(in) :: space
    type(lda_u_t),       intent(in) :: this
    type(states_elec_t), intent(in) :: st
    class(mesh_t),       intent(in) :: mesh
    type(ions_t),        intent(in) :: ions
    logical,             intent(in) :: has_phase
 
    integer :: ios, im, ik, idim, ierr
    complex(real64), allocatable :: tmp(:)
    real(real64), allocatable :: dtmp(:)
    type(orbitalset_t), pointer :: os
    type(unit_t) :: fn_unit
    character(len=MAX_PATH_LEN) :: fname
 
    push_sub(output_dftu_orbitals)
 
    fn_unit = sqrt(units_out%length**(-space%dim))
 
    if (this%basis%use_submesh) then
      if (states_are_real(st)) then
        safe_allocate(dtmp(1:mesh%np))
      else
        safe_allocate(tmp(1:mesh%np))
      end if
    end if
 
    do ios = 1, this%norbsets
      os => this%orbsets(ios)
      do ik = st%d%kpt%start, st%d%kpt%end
        do im = 1, this%orbsets(ios)%norbs
          do idim = 1, min(os%ndim, st%d%dim)
            if (st%nik > 1) then
              if (min(os%ndim, st%d%dim) > 1) then
                write(fname, '(a,i1,a,i3.3,a,i8.8,a,i1)') 'orb', im, '-os', ios, '-k', ik, '-sp', idim
              else
                write(fname, '(a,i1,a,i3.3,a,i8.8)') 'orb', im, '-os', ios, '-k', ik
              end if
            else
              if (min(os%ndim, st%d%dim) > 1) then
                write(fname, '(a,i1,a,i3.3,a,i1)') 'orb', im, '-os', ios, '-sp', idim
              else
                write(fname, '(a,i1,a,i3.3)') 'orb', im, '-os', ios
              end if
            end if
            if (has_phase) then
              if (.not. this%basis%use_submesh) then
                call zio_function_output(outp%how(option__output__local_orbitals), dir, fname, namespace, space, &
                  mesh, os%eorb_mesh(1:mesh%np,im,idim,ik), fn_unit, ierr, pos=ions%pos, atoms=ions%atom)
              else
                tmp = m_z0
                call submesh_add_to_mesh(os%sphere, os%eorb_submesh(1:os%sphere%np,idim,im,ik), tmp)
                call zio_function_output(outp%how(option__output__local_orbitals), dir, fname, namespace, space, &
                  mesh, tmp, fn_unit, ierr, pos=ions%pos, atoms=ions%atom)
              end if
            else
              if (.not. this%basis%use_submesh) then
                if (states_are_real(st)) then
                  call dio_function_output(outp%how(option__output__local_orbitals), dir, fname, namespace, space, mesh, &
                    os%dorb(1:mesh%np,idim,im), fn_unit, ierr, pos=ions%pos, atoms=ions%atom)
                else
                  call zio_function_output(outp%how(option__output__local_orbitals), dir, fname, namespace, space, mesh, &
                    os%zorb(1:mesh%np,idim,im), fn_unit, ierr, pos=ions%pos, atoms=ions%atom)
                end if
              else
                if (states_are_real(st)) then
                  dtmp = m_zero
                  call submesh_add_to_mesh(os%sphere, os%dorb(1:os%sphere%np,idim,im), dtmp)
                  call dio_function_output(outp%how(option__output__local_orbitals), dir, fname, namespace, space, &
                    mesh, dtmp, fn_unit, ierr, pos=ions%pos, atoms=ions%atom)
                else
                  tmp = m_z0
                  call submesh_add_to_mesh(os%sphere, os%zorb(1:os%sphere%np,idim,im), tmp)
                  call zio_function_output(outp%how(option__output__local_orbitals), dir, fname, namespace, space, &
                    mesh, tmp, fn_unit, ierr, pos=ions%pos, atoms=ions%atom)
                end if
              end if
            end if
          end do
        end do
      end do
    end do
 
    safe_deallocate_a(tmp)
    safe_deallocate_a(dtmp)
 
    pop_sub(output_dftu_orbitals)
  end subroutine output_dftu_orbitals
 
  ! ---------------------------------------------------------
  logical function output_needs_current(outp, states_are_real)
    type(output_t),      intent(in) :: outp
    logical,             intent(in) :: states_are_real
 
    output_needs_current = .false.
 
    if (outp%what(option__output__current) &
      .or. outp%what(option__output__current_dia) &
      .or. outp%what(option__output__heat_current) &
      .or. outp%what(option__output__current_kpt)) then
      if (.not. states_are_real) then
        output_needs_current = .true.
      else
        message(1) = 'No current density output for real states since it is identically zero.'
        call messages_warning(1)
      end if
    end if
 
 
  end function
 
#include "output_etsf_inc.F90"
 
#include "output_states_inc.F90"
 
#include "output_h_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "output_linear_response_inc.F90"
 
#include "undef.F90"
#include "real.F90"
#include "output_linear_response_inc.F90"
 
end module output_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: