output_mxll.F90

Go to the documentation of this file.

!! Copyright (C) 2019 F. Bonafe, R. Jestaedt, H. Appel
!!
!! 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_mxll_oct_m
  use debug_oct_m
  use external_densities_oct_m
  use energy_mxll_oct_m
  use hamiltonian_mxll_oct_m
  use helmholtz_decomposition_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_oct_m
  use io_oct_m
  use io_function_oct_m
  use mesh_oct_m
  use messages_oct_m
  use namespace_oct_m
  use output_oct_m
  use output_low_oct_m
  use parser_oct_m
  use profiling_oct_m
  use space_oct_m
  use states_elec_oct_m
  use states_mxll_oct_m
  use string_oct_m
  use unit_oct_m
  use unit_system_oct_m
 
 
  implicit none
 
  private
  public ::              &
    output_mxll_init,    &
    output_mxll
 
contains
 
  ! ---------------------------------------------------------
  subroutine output_mxll_init(outp, namespace, space)
    type(output_t),            intent(out) :: outp
    type(namespace_t),         intent(in)  :: namespace
    class(space_t),            intent(in)  :: space
    integer :: what_i
    push_sub(output_mxll_init)
 
    !%Variable MaxwellOutput
    !%Type block
    !%Default none
    !%Section Output
    !%Description
    !% Specifies what to print. The output files are written at the end of the run into the output directory for the
    !% Maxwell run.
    !% Time-dependent simulations print only per iteration, including always the last. The frequency of output per iteration
    !% is set by <tt>OutputInterval</tt> and the directory is set by <tt>OutputIterDir</tt>.
    !% Each option must be in a separate row. Optionally individual output formats and output intervals can be defined
    !% for each row or they can be read separately from <tt>OutputFormat</tt> and <tt>MaxwellOutputInterval</tt> variables
    !% in the input file.
    !%
    !% Example:
    !% <br><br><tt>%MaxwellOutput
    !% <br>&nbsp;&nbsp;electric_field
    !% <br>&nbsp;&nbsp;magnetic_field
    !% <br>%<br></tt>
    !% This block supports all the formats of the <tt>Output</tt> block.
    !% See <tt>Output</tt>.
    !%Option electric_field 1
    !% Output of the electric field
    !%Option magnetic_field 2
    !% Output of the magnetic field
    !%Option trans_electric_field 3
    !% Output of the transversal electric field
    !%Option trans_magnetic_field 4
    !% Output of the transversal magnetic field
    !%Option long_electric_field 5
    !% Output of the longitudinal electric field
    !%Option long_magnetic_field 6
    !% Output of the longitudinal magnetic field
    !%Option div_electric_field 7
    !% Output of the divergence of the electric field
    !%Option div_magnetic_field 8
    !% Output of the divergence of the magnetic field
    !%Option poynting_vector 9
    !% Output of the Maxwell Poynting vector
    !%Option maxwell_energy_density 10
    !% Output of the electromagnetic density
    !%Option external_current 11
    !% Output of the external Maxwell current
    !%Option charge_density 12
    !% Output of the charge density calculated by the divergence of the electric field.
    !%Option orbital_angular_momentum 13
    !% Output of the orbital angular momentum
    !%Option vector_potential_mag 14
    !% Output of the vector potential from magnetic field
    !%Option magnetic_field_diff 15
    !% Output of the magnetic field difference
    !%Option total_current_mxll 16
    !% Output of the total current density
    !%End
 
    !%Variable MaxwellOutputInterval
    !%Type integer
    !%Default 50
    !%Section Output
    !%Description
    !% The output requested by variable <tt>MaxwellOutput</tt> is written
    !% to the directory <tt>MaxwellOutputIterDir</tt>
    !% when the iteration number is a multiple of the <tt>MaxwellOutputInterval</tt> variable.
    !% Subdirectories are named Y.X, where Y is <tt>td</tt>, <tt>scf</tt>, or <tt>unocc</tt>, and
    !% X is the iteration number. To use the working directory, specify <tt>"."</tt>
    !% (Output of restart files is instead controlled by <tt>MaxwellRestartWriteInterval</tt>.)
    !% Must be >= 0. If it is 0, then no output is written.
    !% This variable can also be defined inside the <tt>MaxwellOutput</tt> block.
    !% See <tt>MaxwellOutput</tt>.
    !%End
 
    outp%what = .false.
    call io_function_read_what_how_when(namespace, space, outp%what, outp%how, outp%output_interval, &
      'MaxwellOutput', 'OutputFormat', 'MaxwellOutputInterval')
 
    do what_i = lbound(outp%what, 1), ubound(outp%what, 1)
      ! xyz format is not available for Maxwell
      if (bitand(outp%how(what_i), option__outputformat__xyz) /= 0) then
        message(1) = "OutputFormat = xyz is not compatible with Maxwell systems"
        call messages_warning(1)
      end if
    end do
 
 
 
    !%Variable MaxwellOutputIterDir
    !%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>MaxwellOutput</tt>
    !% in the format specified by <tt>OutputHow</tt>.
    !% This information is written while iterating <tt>CalculationMode = maxwell</tt>
    !% according to <tt>OutputInterval</tt>, and has nothing to do with the restart information.
    !%End
    call parse_variable(namespace, 'MaxwellOutputIterDir', "output_iter", outp%iter_dir)
    if (any(outp%what) .and. maxval(outp%output_interval) > 0) then
      call io_mkdir(outp%iter_dir, namespace)
    end if
    call add_last_slash(outp%iter_dir)
 
    !%Variable MaxwellRestartWriteInterval
    !%Type integer
    !%Default 50
    !%Section Execution::IO
    !%Description
    !% Restart data is written when the iteration number is a multiple of the
    !% <tt>MaxwellRestartWriteInterval</tt> variable. (Other output is controlled by <tt>MaxwellOutputInterval</tt>.)
    !%End
    call parse_variable(namespace, 'MaxwellRestartWriteInterval', 50, outp%restart_write_interval)
    if (outp%restart_write_interval <= 0) then
      message(1) = "MaxwellRestartWriteInterval must be > 0."
      call messages_fatal(1, namespace=namespace)
    end if
 
    if (outp%what(option__maxwelloutput__electric_field)) then
      outp%wfs_list = trim("1-3")
    end if
 
    if (outp%what(option__maxwelloutput__magnetic_field)) then
      outp%wfs_list = trim("1-3")
    end if
 
    if (outp%what(option__maxwelloutput__trans_electric_field)) then
      outp%wfs_list = trim("1-3")
    end if
 
    if (outp%what(option__maxwelloutput__trans_magnetic_field)) then
      outp%wfs_list = trim("1-3")
    end if
 
 
    pop_sub(output_mxll_init)
  end subroutine output_mxll_init
 
 
  ! ---------------------------------------------------------
  subroutine output_mxll(outp, namespace, space, gr_mxll, st_mxll, hm_mxll, helmholtz, time, dir, gr_elec, st_elec, hm_elec)
    type(output_t),                     intent(in)    :: outp
    type(namespace_t),                  intent(in)    :: namespace
    class(space_t),                     intent(in)    :: space
    type(states_mxll_t),                intent(inout) :: st_mxll
    type(hamiltonian_mxll_t),           intent(in)    :: hm_mxll
    type(grid_t),                       intent(in)    :: gr_mxll
    type(helmholtz_decomposition_t),    intent(inout) :: helmholtz
    real(real64),                       intent(in)    :: time
    character(len=*),                   intent(in)    :: dir
    type(grid_t),             optional, intent(inout) :: gr_elec
    type(states_elec_t),      optional, intent(inout) :: st_elec
    type(hamiltonian_elec_t), optional, intent(in)    :: hm_elec
 
    push_sub(output_mxll)
 
    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
 
    call output_states_mxll(outp, namespace, space, dir, st_mxll, gr_mxll)
    call output_energy_density_mxll(outp, namespace, space, dir, st_mxll, gr_mxll)
    call output_poynting_vector_orbital_angular_momentum(outp, namespace, space, dir, st_mxll, gr_mxll)
    call output_transverse_rs_state(helmholtz, outp, namespace, space, dir, st_mxll, gr_mxll)
    call output_longitudinal_rs_state(helmholtz, outp, namespace, space, dir, st_mxll, gr_mxll)
    call output_vector_potential_mag(outp, helmholtz, namespace, gr_mxll, space, dir, st_mxll)
    call output_divergence_rs_state(outp, namespace, space, dir, st_mxll, gr_mxll)
    call output_external_current_density(outp, namespace, space, dir, st_mxll, gr_mxll, hm_mxll, time)
    call output_total_current_mxll(outp, namespace, space, dir, st_mxll, gr_mxll, hm_mxll, time, st_mxll%ep, gr_mxll%np)
    call output_charge_density_mxll(outp, namespace, space, dir, st_mxll, gr_mxll, hm_mxll)
 
    if (present(hm_elec) .and. present(gr_elec) .and. present(st_elec)) then
      call output_coupling_potentials(outp, namespace, dir, hm_elec, gr_elec)
      call output_current_density(outp, namespace, dir, st_mxll, gr_mxll, hm_mxll, st_elec, gr_elec, hm_elec, time)
    end if
 
    pop_sub(output_mxll)
  end subroutine output_mxll
 
 
  ! ---------------------------------------------------------
  subroutine output_states_mxll(outp, namespace, space, dir, st, mesh)
    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_mxll_t),      intent(in)    :: st
    class(mesh_t),            intent(in)    :: mesh
 
    integer :: ierr
    character(len=MAX_PATH_LEN) :: fname
    type(unit_t) :: fn_unit
    real(real64), allocatable :: dtmp(:,:)
 
    push_sub(output_states_mxll)
 
    ! Electric field
    if (outp%what(option__maxwelloutput__electric_field)) then
      fn_unit = units_out%energy/units_out%length
      safe_allocate(dtmp(1:mesh%np, 1:space%dim))
      call get_electric_field_state(st%rs_state, mesh, dtmp, st%ep, mesh%np)
      fname = 'e_field'
      call io_function_output_vector(outp%how(option__maxwelloutput__electric_field), dir, fname, namespace, space, &
        mesh, dtmp, fn_unit, ierr)
      safe_deallocate_a(dtmp)
    end if
 
    ! Magnetic field
    if (outp%what(option__maxwelloutput__magnetic_field)) then
      fn_unit = unit_one/units_out%length**2
      safe_allocate(dtmp(1:mesh%np, 1:space%dim))
      call get_magnetic_field_state(st%rs_state, mesh, st%rs_sign, dtmp, st%mu(1:mesh%np), mesh%np)
      fname = 'b_field'
      call io_function_output_vector(outp%how(option__maxwelloutput__magnetic_field), dir, fname, namespace, space, &
        mesh, dtmp, fn_unit, ierr)
      safe_deallocate_a(dtmp)
    end if
 
    pop_sub(output_states_mxll)
  end subroutine output_states_mxll
 
 
  !----------------------------------------------------------
  subroutine output_energy_density_mxll(outp, namespace, space, dir, st, mesh)   !< have to set unit output correctly
    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_mxll_t),  intent(in)    :: st
    class(mesh_t),        intent(in)    :: mesh
 
    integer :: ierr
    real(real64), allocatable :: energy_density(:), e_energy_density(:), b_energy_density(:)
 
    push_sub(output_energy_density_mxll)
 
    ! Maxwell energy density
    if (outp%what(option__maxwelloutput__maxwell_energy_density)) then
      safe_allocate(energy_density(1:mesh%np))
      safe_allocate(e_energy_density(1:mesh%np))
      safe_allocate(b_energy_density(1:mesh%np))
      ! calculate the energy density
      call energy_density_calc(mesh, st, st%rs_state, energy_density, e_energy_density, b_energy_density)
 
      call dio_function_output(outp%how(option__maxwelloutput__maxwell_energy_density), dir, "maxwell_energy_density", &
        namespace, space, mesh, energy_density, units_out%energy/units_out%length**3, ierr)
 
      safe_deallocate_a(energy_density)
      safe_deallocate_a(e_energy_density)
      safe_deallocate_a(b_energy_density)
    end if
 
    pop_sub(output_energy_density_mxll)
  end subroutine output_energy_density_mxll
 
 
  !----------------------------------------------------------
  subroutine output_poynting_vector_orbital_angular_momentum(outp, namespace, space, dir, st, mesh)
    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_mxll_t),       intent(in)    :: st
    class(mesh_t),             intent(in)    :: mesh
 
    integer :: ierr
    character(len=MAX_PATH_LEN) :: fname
    type(unit_t) :: fn_unit
    real(real64), allocatable :: poynting_vector(:,:), oam(:,:)
 
    push_sub(output_poynting_vector_orbital_angular_momentum)
 
    if (outp%what(option__maxwelloutput__poynting_vector).or.outp%what(option__maxwelloutput__orbital_angular_momentum)) then
 
      ! the Poynting vector is needed in both cases, so we compute it here first
      fn_unit = units_out%energy/(units_out%time*units_out%length**2)
      safe_allocate(poynting_vector(1:mesh%np, 1:space%dim))
      call get_poynting_vector(mesh, st, st%rs_state, st%rs_sign, poynting_vector, st%ep, st%mu)
 
      if(outp%what(option__maxwelloutput__poynting_vector)) then
        fname = 'poynting_vector'
        call io_function_output_vector(outp%how(option__maxwelloutput__poynting_vector), dir, fname, namespace, space, &
          mesh, poynting_vector, fn_unit, ierr)
      end if
 
      if(outp%what(option__maxwelloutput__orbital_angular_momentum)) then
        safe_allocate(oam(1:mesh%np, 1:space%dim))
        call get_orbital_angular_momentum(mesh, st, poynting_vector, oam)
        fname = 'orbital_angular_momentum'
        call io_function_output_vector(outp%how(option__maxwelloutput__orbital_angular_momentum), dir, fname, namespace, space, &
          mesh, oam, fn_unit, ierr)
        safe_deallocate_a(oam)
      end if
 
      safe_deallocate_a(poynting_vector)
    end if
 
    pop_sub(output_poynting_vector_orbital_angular_momentum)
  end subroutine output_poynting_vector_orbital_angular_momentum
 
 
  ! add documentation to magnetic field subtraction order
  subroutine output_vector_potential_mag(outp, helmholtz, namespace, gr, space, dir, st)    !< have to set unit output correctly
    type(output_t),                  intent(in)    :: outp
    type(helmholtz_decomposition_t), intent(inout) :: helmholtz
    type(namespace_t),               intent(in)    :: namespace
    type(grid_t),                    intent(in)    :: gr
    class(space_t),                  intent(in)    :: space
    character(len=*),                intent(in)    :: dir
    type(states_mxll_t),             intent(in)    :: st
 
    character(len=MAX_PATH_LEN) :: fname
    integer :: ierr
    type(unit_t) :: fn_unit
    real(real64), allocatable :: dtmp(:,:),vtmp(:,:), mag_fromvec(:,:) , delta(:,:)
 
    push_sub(output_vector_potential_mag)
 
    ! vector potential from magnetic field
    if (outp%what(option__maxwelloutput__vector_potential_mag) .or. &
      outp%what(option__maxwelloutput__magnetic_field_diff)) then
 
      fn_unit = unit_one/units_out%length
      safe_allocate(dtmp(1:gr%np_part, 1:space%dim))
      safe_allocate(vtmp(1:gr%np_part, 1:space%dim))
      call get_magnetic_field_state(st%rs_state, gr, st%rs_sign, dtmp, st%mu(1:gr%np), gr%np)
      ! Get vector potential for magnetic field
      call helmholtz%get_vector_potential(namespace, vtmp, trans_field=dtmp)
 
      if (outp%what(option__maxwelloutput__vector_potential_mag)) then
        message(1) = 'Vector potential is currently missing surface contributions'
        message(2) = 'If in doubt, use magnetic_field_diff output which shows deviation of B field'
        message(3) = 'Large B field deviations mean the calculated vector potential is unreliable'
        call messages_warning(3)
 
        fname = 'vector_potential_mag'
        call io_function_output_vector(outp%how(option__maxwelloutput__vector_potential_mag), dir, fname, namespace, &
          space, gr, vtmp, fn_unit, ierr)
      end if
 
      if (outp%what(option__maxwelloutput__magnetic_field_diff)) then
        fn_unit = unit_one/units_out%length**2
        safe_allocate(mag_fromvec(1:gr%np_part, 1:space%dim))
        safe_allocate(delta(1:gr%np, 1:space%dim))
        ! Get transverse field from vector potential
        call helmholtz%get_trans_field(namespace, mag_fromvec, vector_potential=vtmp)
        ! mag_fromvec is of size np_part, but should be used up to np
        ! Due to how the Helmholtz decomposition works (see its documentation), we have to remove a stencil from dtmp
        call helmholtz%apply_inner_stencil_mask(dtmp)
        delta = dtmp(1:gr%np, 1:space%dim) - mag_fromvec(1:gr%np, 1:space%dim)
        fname = 'magnetic_field_diff'
        call io_function_output_vector(outp%how(option__maxwelloutput__magnetic_field_diff), dir, fname, namespace, &
          space, gr, delta, fn_unit, ierr)
        safe_deallocate_a(mag_fromvec)
        safe_deallocate_a(delta)
      end if
 
      safe_deallocate_a(dtmp)
      safe_deallocate_a(vtmp)
    end if
 
    pop_sub(output_vector_potential_mag)
  end subroutine output_vector_potential_mag
 
 
  !----------------------------------------------------------
  subroutine output_transverse_rs_state(helmholtz, outp, namespace, space, dir, st, gr)    !< have to set unit output correctly
    type(helmholtz_decomposition_t), intent(inout) :: helmholtz
    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_mxll_t),             intent(inout) :: st  ! needs to be inout as st%rs_state is passed to get_trans_field
    type(grid_t),                    intent(in)    :: gr
 
    character(len=MAX_PATH_LEN) :: fname
    integer :: ierr
    type(unit_t) :: fn_unit
    real(real64), allocatable :: dtmp(:,:)
 
    push_sub(output_transverse_rs_state)
 
    if (outp%what(option__maxwelloutput__trans_electric_field) .or. outp%what(option__maxwelloutput__trans_magnetic_field)) then
      call helmholtz%get_trans_field(namespace, st%rs_state_trans, total_field=st%rs_state)
    end if
 
    ! transverse component of the electric field
    if (outp%what(option__maxwelloutput__trans_electric_field)) then
      fn_unit = units_out%energy/units_out%length
      safe_allocate(dtmp(1:gr%np, 1:space%dim))
      call get_electric_field_state(st%rs_state_trans, gr, dtmp, st%ep(1:gr%np), gr%np)
      fname = 'e_field_trans'
      call io_function_output_vector(outp%how(option__maxwelloutput__trans_electric_field), dir, fname, namespace, space, &
        gr, dtmp, fn_unit, ierr)
      safe_deallocate_a(dtmp)
    end if
 
    ! transverse component of the magnetic field
    if (outp%what(option__maxwelloutput__trans_magnetic_field)) then
      fn_unit = unit_one/units_out%length**2
      safe_allocate(dtmp(1:gr%np, 1:space%dim))
      call get_magnetic_field_state(st%rs_state_trans, gr, st%rs_sign, dtmp, st%mu(1:gr%np), gr%np)
      fname = 'b_field_trans'
      call io_function_output_vector(outp%how(option__maxwelloutput__trans_magnetic_field), dir, fname, namespace, space, &
        gr, dtmp, fn_unit, ierr)
      safe_deallocate_a(dtmp)
    end if
 
    pop_sub(output_transverse_rs_state)
  end subroutine output_transverse_rs_state
 
 
  !----------------------------------------------------------
  !This subroutine is left here, will be corrected and called in near future
  subroutine output_longitudinal_rs_state(helmholtz, outp, namespace, space, dir, st, gr)    !< have to set unit output correctly
    type(helmholtz_decomposition_t), intent(inout) :: helmholtz
    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_mxll_t),       intent(inout) :: st   ! needs to be inout as st%rs_state is passed to get_long_field
    type(grid_t),              intent(in)    :: gr
 
    character(len=MAX_PATH_LEN) :: fname
    integer :: ierr
    type(unit_t) :: fn_unit
    real(real64), allocatable :: dtmp(:,:)
 
    push_sub(output_longitudinal_rs_state)
 
    if (outp%what(option__maxwelloutput__long_electric_field) .or. outp%what(option__maxwelloutput__long_magnetic_field)) then
      call helmholtz%get_long_field(namespace, st%rs_state_long, total_field=st%rs_state)
    end if
 
    ! longitudinal component of the electric field
    if (outp%what(option__maxwelloutput__long_electric_field)) then
      fn_unit = units_out%energy/units_out%length
      safe_allocate(dtmp(1:gr%np, 1:space%dim))
      call get_electric_field_state(st%rs_state_long, gr, dtmp, st%ep(1:gr%np), gr%np)
      fname = 'e_field_long'
      call io_function_output_vector(outp%how(option__maxwelloutput__long_electric_field), dir, fname, namespace, space, &
        gr, dtmp, fn_unit, ierr)
      safe_deallocate_a(dtmp)
    end if
 
    ! longitudinal component of the magnetic field
    if (outp%what(option__maxwelloutput__long_magnetic_field)) then
      fn_unit = unit_one/units_out%length**2
      safe_allocate(dtmp(1:gr%np, 1:space%dim))
      call get_magnetic_field_state(st%rs_state_long, gr, st%rs_sign, dtmp, st%mu(1:gr%np), gr%np)
      fname = 'b_field_long'
      call io_function_output_vector(outp%how(option__maxwelloutput__long_magnetic_field), dir, fname, namespace, space, &
        gr, dtmp, fn_unit, ierr)
      safe_deallocate_a(dtmp)
    end if
 
    pop_sub(output_longitudinal_rs_state)
  end subroutine output_longitudinal_rs_state
 
 
  !----------------------------------------------------------
  subroutine output_divergence_rs_state(outp, namespace, space, dir, st, gr)    !< have to set unit output correctly
    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_mxll_t),       intent(in)    :: st
    type(grid_t),              intent(in)    :: gr
 
    integer :: ierr
    type(unit_t) :: fn_unit
    real(real64), allocatable :: dtmp_1(:,:), dtmp_2(:)
 
    push_sub(output_divergence_rs_state)
 
    ! divergence of the electric field
    if (outp%what(option__maxwelloutput__div_electric_field)) then
      fn_unit = units_out%length**(-space%dim)
      safe_allocate(dtmp_1(1:gr%np_part, 1:space%dim))
      safe_allocate(dtmp_2(1:gr%np))
      dtmp_1 = m_zero
      call get_electric_field_state(st%rs_state, gr, dtmp_1(1:gr%np,:), st%ep(1:gr%np), gr%np)
      call get_divergence_field(gr, dtmp_1, dtmp_2, .false.)
      call dio_function_output(outp%how(option__maxwelloutput__div_electric_field), dir, "e_field_div", namespace, space, &
        gr, dtmp_2, fn_unit, ierr)
      safe_deallocate_a(dtmp_1)
      safe_deallocate_a(dtmp_2)
    end if
 
    ! divergence of the magnetic field
    if (outp%what(option__maxwelloutput__div_magnetic_field)) then
      ! unit does not matter, should be zero
      fn_unit = unit_one
      safe_allocate(dtmp_1(1:gr%np_part, 1:space%dim))
      safe_allocate(dtmp_2(1:gr%np))
      dtmp_1 = m_zero
      call get_magnetic_field_state(st%rs_state, gr, st%rs_sign, dtmp_1(1:gr%np,:), st%mu(1:gr%np), gr%np)
      call get_divergence_field(gr, dtmp_1, dtmp_2, .false.)
      call dio_function_output(outp%how(option__maxwelloutput__div_magnetic_field), dir, "b_field_div", namespace, space, &
        gr, dtmp_2, fn_unit, ierr)
      safe_deallocate_a(dtmp_1)
      safe_deallocate_a(dtmp_2)
    end if
 
    pop_sub(output_divergence_rs_state)
  end subroutine output_divergence_rs_state
 
 
  !------------------------------------------------------------
  subroutine output_coupling_potentials(outp, namespace, dir, hm, gr)    !< have to set unit output correctly
    type(output_t),            intent(in)    :: outp
    type(namespace_t),         intent(in)    :: namespace
    character(len=*),          intent(in)    :: dir
    type(hamiltonian_elec_t),  intent(in)    :: hm
    type(grid_t),              intent(in)    :: gr
 
    message(1) = "Maxwell-matter coupling potentials not implemented yet."
    call messages_fatal(1, namespace=namespace)
 
  end subroutine output_coupling_potentials
 
 
  !----------------------------------------------------------
  subroutine output_current_density(outp, namespace, dir, st_mxll, gr_mxll, hm_mxll, st_elec, gr_elec, hm_elec, time)
    type(output_t),            intent(in)    :: outp
    type(namespace_t),         intent(in)    :: namespace
    character(len=*),          intent(in)    :: dir
    type(states_mxll_t),       intent(in)    :: st_mxll
    type(grid_t),              intent(in)    :: gr_mxll
    type(hamiltonian_mxll_t),  intent(in)    :: hm_mxll
    type(states_elec_t),       intent(in)    :: st_elec
    type(grid_t),              intent(in)    :: gr_elec
    type(hamiltonian_elec_t),  intent(in)    :: hm_elec
    real(real64),              intent(in)    :: time
 
    message(1) = "Current density can not yet be calculated in a Maxwell propagation."
    call messages_fatal(1, namespace=namespace)
 
  end subroutine output_current_density
 
 
  !----------------------------------------------------------
  subroutine output_external_current_density(outp, namespace, space, dir, st, mesh, hm, time)
    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_mxll_t),       intent(in)    :: st
    class(mesh_t),             intent(in)    :: mesh
    type(hamiltonian_mxll_t),  intent(in)    :: hm
    real(real64),              intent(in)    :: time
 
    character(len=MAX_PATH_LEN) :: fname
    integer                     :: ierr
    real(real64), allocatable          :: dtmp(:,:)
    type(unit_t)                :: fn_unit
 
    push_sub(output_external_current_density)
 
    ! Maxwell current density
    if (outp%what(option__maxwelloutput__external_current)) then
      if (hm%current_density_ext_flag) then
        fn_unit = (unit_one/units_out%time)/(units_out%length**2)      
        safe_allocate(dtmp(1:mesh%np, 1:space%dim))
        call external_current_calculation(st, space, mesh, time, dtmp)
        fname = 'external_current'
        call io_function_output_vector(outp%how(option__maxwelloutput__external_current), dir, fname, namespace, space, &
          mesh, dtmp, fn_unit, ierr)
        safe_deallocate_a(dtmp)
      end if
    end if
 
    pop_sub(output_external_current_density)
  end subroutine output_external_current_density
 
  subroutine output_total_current_mxll(outp, namespace, space, dir, st, mesh, hm, time, ep_field, np)
    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_mxll_t),       intent(in)    :: st
    class(mesh_t),             intent(in)    :: mesh
    type(hamiltonian_mxll_t),  intent(in)    :: hm
    real(real64),              intent(in)    :: time
    real(real64),   optional, intent(in)    :: ep_field(:)
    integer, optional, intent(in)    :: np
 
    real(real64), allocatable          :: ctmp(:,:)
    character(len=MAX_PATH_LEN) :: fname
    integer                     :: ierr, ip, idim, np_, ff_dim
    type(unit_t)                :: fn_unit
 
    push_sub(output_total_current_mxll)
    safe_allocate(ctmp(1:mesh%np, 1:space%dim))
    if (outp%what(option__maxwelloutput__total_current_mxll)) then
      np_ = optional_default(np, mesh%np)
      ff_dim = space%dim
      if (present(ep_field)) then
        do idim = 1, ff_dim
          do ip = 1, np_
            ctmp(ip, idim) = sqrt(m_two*ep_field(ip)) * real(st%rs_current_density_t1(ip, idim))
          end do
        end do
      else
        do idim = 1, ff_dim
          do ip = 1, np_
            ctmp(ip, idim) = sqrt(m_two*p_ep) * real(st%rs_current_density_t1(ip, idim))
          end do
        end do
      end if
 
      fn_unit = (unit_one/units_out%time)/(units_out%length**2)      
      fname = 'total_current_mxll'
      call io_function_output_vector(outp%how(option__maxwelloutput__total_current_mxll), dir, fname, namespace, space, &
        mesh, ctmp, fn_unit, ierr)
    end if
    safe_deallocate_a(ctmp)
    pop_sub(output_total_current_mxll)
  end subroutine output_total_current_mxll
 
 
  !----------------------------------------------------------
  subroutine output_charge_density_mxll(outp, namespace, space, dir, st, gr, hm)    !< have to set unit output correctly
    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_mxll_t),           intent(in)    :: st
    type(grid_t),                  intent(in)    :: gr
    type(hamiltonian_mxll_t),      intent(in)    :: hm
 
    integer :: ierr
    type(unit_t) :: fn_unit
    real(real64), allocatable :: dtmp_1(:,:), dtmp_2(:)
 
    push_sub(output_charge_density_mxll)
 
    ! charge density calculated by the divergence of the electric field
    if (outp%what(option__maxwelloutput__charge_density)) then
      fn_unit = units_out%length**(-space%dim)
      safe_allocate(dtmp_1(1:gr%np_part,1:space%dim))
      safe_allocate(dtmp_2(1:gr%np))
      call get_electric_field_state(st%rs_state, gr, dtmp_1, st%ep, gr%np)
      call get_divergence_field(gr, dtmp_1, dtmp_2, .true.)
      call dio_function_output(outp%how(option__maxwelloutput__charge_density), dir, "charge_density", namespace, space, &
        gr, dtmp_2(:), fn_unit, ierr)
      safe_deallocate_a(dtmp_1)
      safe_deallocate_a(dtmp_2)
    end if
 
    pop_sub(output_charge_density_mxll)
  end subroutine output_charge_density_mxll
 
end module output_mxll_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: