td_write.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 td_write_oct_m
  use blas_oct_m
  use comm_oct_m
  use current_oct_m
  use debug_oct_m
  use density_oct_m
  use derivatives_oct_m
  use elec_matrix_elements_oct_m
  use electron_space_oct_m
  use excited_states_oct_m
  use ext_partner_list_oct_m
  use gauge_field_oct_m
  use global_oct_m
  use grid_oct_m
  use output_oct_m
  use output_low_oct_m
  use output_mxll_oct_m
  use output_modelmb_oct_m
  use hamiltonian_elec_oct_m
  use hamiltonian_mxll_oct_m
  use helmholtz_decomposition_m
  use interaction_partner_oct_m
  use io_function_oct_m
  use io_oct_m
  use ion_dynamics_oct_m
  use ions_oct_m
  use kick_oct_m
  use, intrinsic :: iso_fortran_env
  use kpoints_oct_m
  use lasers_oct_m
  use lalg_basic_oct_m
  use lalg_adv_oct_m
  use lda_u_oct_m
  use loct_math_oct_m
  use magnetic_oct_m
  use math_oct_m
  use mesh_batch_oct_m
  use mesh_function_oct_m
  use mesh_oct_m
  use messages_oct_m
  use modelmb_exchange_syms_oct_m
  use mpi_oct_m
  use mpi_lib_oct_m
  use multicomm_oct_m
  use mxll_elec_coupling_oct_m
  use namespace_oct_m
  use parser_oct_m
  use partial_charges_oct_m
  use perturbation_oct_m
  use perturbation_magnetic_oct_m
  use profiling_oct_m
  use restart_oct_m
  use space_oct_m
  use states_elec_oct_m
  use states_elec_calc_oct_m
  use states_elec_dim_oct_m
  use states_elec_restart_oct_m
  use states_mxll_oct_m
  use states_mxll_restart_oct_m
  use td_calc_oct_m
  use td_write_low_oct_m
  use td_write_classical_oct_m
  use types_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 write_iter_oct_m
 
  implicit none
 
  private
  public ::         &
    td_write_t,     &
    td_write_init,  &
    td_write_end,   &
    td_write_iter,  &
    td_write_data,  &
    td_write_kick,  &
    td_write_output, &
    td_write_mxll_end, &
    td_write_mxll_init, &
    td_write_mxll_iter, &
    td_write_mxll_free_data
 
  integer, parameter, public ::   &
    OUT_MULTIPOLES         =  1, &
    out_angular            =  2, &
    out_spin               =  3, &
    out_populations        =  4, &
    out_coords             =  5, &
    out_acc                =  6, &
    out_laser              =  7, &
    out_energy             =  8, &
    out_proj               =  9, &
    out_magnets            = 10, &
    out_gauge_field        = 11, &
    out_temperature        = 12, &
    out_ftchd              = 13, &
    out_vel                = 14, &
    out_eigs               = 15, &
    out_ion_ch             = 16, &
    out_total_current      = 17, &
    out_partial_charges    = 18, &
    out_kp_proj            = 19, &
    out_floquet            = 20, &
    out_n_ex               = 21, &
    out_separate_coords    = 22, &
    out_separate_velocity  = 23, &
    out_separate_forces    = 24, &
    out_total_heat_current = 25, &
    out_tot_m              = 26, &
    out_q                  = 27, &
    out_mxll_field         = 28, &
    out_norm_ks            = 29, &
  ! This constant should always equal the last integer above
  ! as it defines the total numer of output options
    out_max                = 29
 
 
  ! Entries must be consistent with the ordering above
  character(len=100) :: td_file_name(OUT_MAX) = [character(100) :: &
    "NULL", &          ! Do not include OUT_MULTIPOLES (extension is a function of state index)
    "angular", &
    "spin", &
    "populations", &
    "coordinates", &
    "acceleration", &
    "NULL", &          ! Do not include OUT_LASER
    "energy", &
    "projections", &
    "magnetic_moments", &
    "gauge_field", &
    "temperature", &
    "NULL", &          ! Do not include OUT_FTCHD
    "velocity", &
    "eigenvalues", &
    "ion_ch", &
    "total_current", &
    "partial_charges", &
    "projections", &
    "floquet_bands", &
    "n_ex", &
    "onlyCoordinates", &
    "onlyVelocities", &
    "onlyForces", &
    "total_heat_current", &
    "total_magnetization", &
    "photons_q", &
    "maxwell_dipole_field", &
    "norm_wavefunctions"]
 
  integer, parameter ::       &
    OUT_DFTU_EFFECTIVE_U = 1, &
    out_dftu_max         = 1
 
  
  integer, parameter ::   &
    OUT_MAXWELL_TOTAL_E_FIELD =  1, &
    out_maxwell_total_b_field =  4, &
    out_maxwell_long_e_field  =  7, &
    out_maxwell_long_b_field  = 10, &
    out_maxwell_trans_e_field = 13, &
    out_maxwell_trans_b_field = 16, &
    out_maxwell_energy        = 19, &
    out_e_field_surface_x     = 20, &
    out_e_field_surface_y     = 21, &
    out_e_field_surface_z     = 22, &
    out_b_field_surface_x     = 23, &
    out_b_field_surface_y     = 24, &
    out_b_field_surface_z     = 25
 
  integer, parameter, public ::   &
    OUT_MAXWELL_MAX             = 25
 
  integer, parameter ::      &
    MAXWELL_TOTAL_FIELD = 1, &
    maxwell_long_field  = 2, &
    maxwell_trans_field = 3
 
  integer, parameter ::  &
    maxwell_e_field = 1, &
    maxwell_b_field = 2
 
  type td_write_prop_t
    type(c_ptr)              :: handle
    type(c_ptr), allocatable :: mult_handles(:)
    type(mpi_grp_t)          :: mpi_grp
    integer                  :: hand_start
    integer                  :: hand_end
    logical                  :: write = .false.
    logical                  :: resolve_states = .false.    
  end type td_write_prop_t
 
  type td_write_t
    private
    type(td_write_prop_t), public :: out(out_max)
    type(td_write_prop_t)         :: out_dftu(out_dftu_max)
 
    integer        :: lmax
    real(real64)   :: lmm_r
    type(states_elec_t) :: gs_st
    !                                calculate the projections(s) onto it.
    integer        :: n_excited_states
    type(excited_states_t), allocatable :: excited_st(:)
    integer :: compute_interval
    type(elec_matrix_elements_t) :: elec_me
  end type td_write_t
 
contains
 
 
  subroutine td_write_kick(outp, namespace, space, mesh, kick, ions, iter)
    type(output_t),         intent(in) :: outp
    type(namespace_t),      intent(in) :: namespace
    type(electron_space_t), intent(in) :: space
    class(mesh_t),          intent(in) :: mesh
    type(kick_t),           intent(in) :: kick
    type(ions_t),           intent(in) :: ions
    integer,                intent(in) :: iter
 
    complex(real64), allocatable :: kick_function(:)
    character(len=256) :: filename
    integer :: err
 
    push_sub(td_write_kick)
 
    write(filename, '(a,i7.7)') "td.", iter  ! name of directory
    if (outp%what(option__output__delta_perturbation)) then
      safe_allocate(kick_function(1:mesh%np))
      call kick_function_get(space, mesh, kick, kick_function, 1)
      call zio_function_output(outp%how(option__output__delta_perturbation), filename, "kick_function", namespace, &
        space, mesh, kick_function(:), units_out%energy, err, pos=ions%pos, atoms=ions%atom)
      safe_deallocate_a(kick_function)
    end if
 
    pop_sub(td_write_kick)
  end subroutine td_write_kick
 
 
  subroutine td_write_init(writ, namespace, space, outp, gr, st, hm, ions, ext_partners, ks, ions_move, &
    with_gauge_field, kick, iter, max_iter, dt, mc)
    type(td_write_t), target, intent(out)   :: writ
    type(namespace_t),        intent(in)    :: namespace
    type(electron_space_t),   intent(in)    :: space
    type(output_t),           intent(inout) :: outp
    type(grid_t),             intent(in)    :: gr
    type(states_elec_t),      intent(inout) :: st
    type(hamiltonian_elec_t), intent(inout) :: hm
    type(ions_t),             intent(in)    :: ions
    type(partner_list_t),     intent(in)    :: ext_partners
    type(v_ks_t),             intent(inout) :: ks
    logical,                  intent(in)    :: ions_move
    logical,                  intent(in)    :: with_gauge_field
    type(kick_t),             intent(in)    :: kick
    integer,                  intent(in)    :: iter
    integer,                  intent(in)    :: max_iter
    real(real64),             intent(in)    :: dt
    type(multicomm_t),        intent(in)    :: mc
 
    real(real64) :: rmin
    integer :: ierr, first, ii, ist, jj, flags, iout, default, ifile
    logical :: output_options(max_output_types)
    integer :: output_interval(0:max_output_types)
    integer(int64) :: how(0:max_output_types)
    type(block_t) :: blk
    character(len=MAX_PATH_LEN) :: filename
    type(restart_t) :: restart_gs
    logical :: resolve_states
    logical, allocatable :: skip(:)
 
    push_sub(td_write_init)
 
    !%Variable TDOutput
    !%Type block
    !%Default multipoles + energy (+ others depending on other options)
    !%Section Time-Dependent::TD Output
    !%Description
    !% Defines what should be output during the time-dependent
    !% simulation. Many of the options can increase the computational
    !% cost of the simulation, so only use the ones that you need. In
    !% most cases the default value is enough, as it is adapted to the
    !% details of the TD run. If the ions are allowed to be moved, additionally
    !% the geometry and the temperature are output. If a laser is
    !% included it will output by default.
    !%
    !% Note: the output files generated by this option are updated
    !% every <tt>RestartWriteInterval</tt> steps.
    !%
    !% Example:
    !% <br><br><tt>%TDOutput
    !% <br>&nbsp;&nbsp;multipoles
    !% <br>&nbsp;&nbsp;energy
    !% <br>%<br></tt>
    !%
    !%Option multipoles 1
    !% Outputs the (electric) multipole moments of the density to the file <tt>td.general/multipoles</tt>.
    !% This is required to, <i>e.g.</i>, calculate optical absorption spectra of finite systems. The
    !% maximum value of <math>l</math> can be set with the variable <tt>TDMultipoleLmax</tt>.
    !%Option angular 2
    !% Outputs the orbital angular momentum of the system to <tt>td.general/angular</tt>.
    !% This can be used to calculate circular dicrhoism. Not implemented for periodic systems.
    !% In case of a nonlocal pseudopotential, the gauge used can be specified using
    !% the variable <tt>MagneticGaugeCorrection</tt>.
    !%Option spin 3
    !% (Experimental) Outputs the expectation value of the spin, which can be used to calculate magnetic
    !% circular dichroism.
    !%Option populations 4
    !% (Experimental) Outputs the projection of the time-dependent
    !% Kohn-Sham Slater determinant onto the ground state (or
    !% approximations to the excited states) to the file
    !% <tt>td.general/populations</tt>. Note that the calculation of
    !% populations is expensive in memory and computer time, so it
    !% should only be used if it is really needed. See <tt>TDExcitedStatesToProject</tt>.
    !%Option geometry 5
    !% If set (and if the atoms are allowed to move), outputs the coordinates, velocities,
    !% and forces of the atoms to the the file <tt>td.general/coordinates</tt>. On by default if <tt>MoveIons = yes</tt>.
    !%Option dipole_acceleration 6
    !% When set, outputs the acceleration of the electronic dipole, calculated from the Ehrenfest theorem,
    !% in the file <tt>td.general/acceleration</tt>. This file can then be
    !% processed by the utility <tt>oct-harmonic-spectrum</tt> in order to obtain the harmonic spectrum.
    !%Option laser 7
    !% If set, outputs the laser field to the file <tt>td.general/laser</tt>.
    !% On by default if <tt>TDExternalFields</tt> is set.
    !%Option energy 8
    !% If set, <tt>octopus</tt> outputs the different components of the energy
    !% to the file <tt>td.general/energy</tt>. Will be zero except for every <tt>TDEnergyUpdateIter</tt> iterations.
    !%Option td_occup 9
    !% (Experimental) If set, outputs the projections of the
    !% time-dependent Kohn-Sham wavefunctions onto the static
    !% (zero-time) wavefunctions to the file
    !% <tt>td.general/projections.XXX</tt>. Only use this option if
    !% you really need it, as it might be computationally expensive. See <tt>TDProjStateStart</tt>.
    !% The output interval of this quantity is controled by the variable <tt>TDOutputComputeInterval</tt>
    !% In case of states parallelization, all the ground-state states are stored by each task.
    !%Option local_mag_moments 10
    !% If set, outputs the local magnetic moments, integrated in sphere centered around each atom.
    !% The radius of the sphere can be set with <tt>LocalMagneticMomentsSphereRadius</tt>.
    !%Option gauge_field 11
    !% If set, outputs the vector gauge field corresponding to a spatially uniform (but time-dependent)
    !% external electrical potential. This is only useful in a time-dependent periodic run.
    !% On by default if <tt>GaugeVectorField</tt> is set.
    !%Option temperature 12
    !% If set, the ionic temperature at each step is printed. On by default if <tt>MoveIons = yes</tt>.
    !%Option ftchd 13
    !% Write Fourier transform of the electron density to the file <tt>ftchd.X</tt>,
    !% where X depends on the kick (e.g. with sin-shaped perturbation X=sin).
    !% This is needed for calculating the dynamic structure factor.
    !% In the case that the kick mode is qbessel, the written quantity is integral over
    !% density, multiplied by spherical Bessel function times real spherical harmonic.
    !% On by default if <tt>TDMomentumTransfer</tt> is set.
    !%Option dipole_velocity 14
    !% When set, outputs the dipole velocity, calculated from the Ehrenfest theorem,
    !% in the file <tt>td.general/velocity</tt>. This file can then be
    !% processed by the utility <tt>oct-harmonic-spectrum</tt> in order to obtain the harmonic spectrum.
    !%Option eigenvalues 15
    !% Write the KS eigenvalues.
    !%Option ionization_channels 16
    !% Write the multiple-ionization channels using the KS orbital densities as proposed in
    !% C. Ullrich, Journal of Molecular Structure: THEOCHEM 501, 315 (2000).
    !%Option total_current 17
    !% Output the total current (average of the current density over the cell).
    !%Option partial_charges 18
    !% Bader and Hirshfeld partial charges. The output file is called 'td.general/partial_charges'.
    !%Option td_kpoint_occup 19
    !% Project propagated Kohn-Sham states to the states at t=0 given in the directory
    !% restart_proj (see %RestartOptions). This is an alternative to the option
    !% td_occup, with a formating more suitable for k-points and works only in
    !% k- and/or state parallelization
    !%Option td_floquet 20
    !% Compute non-interacting Floquet bandstructure according to further options:
    !% TDFloquetFrequency, TDFloquetSample, TDFloquetDimension.
    !% This is done only once per td-run at t=0.
    !% works only in k- and/or state parallelization
    !%Option n_excited_el 21
    !% Output the number of excited electrons, based on the projections
    !% of the time evolved wave-functions on the ground-state wave-functions.
    !% The output interval of this quantity is controled by the variable <tt>TDOutputComputeInterval</tt>
    !% Note that if one sets RecalculateGSDuringEvolution=yes,
    !% the code will recompute the GS states
    !% and use them for the computing the number of excited electrons.
    !% This is useful if ions are moving or if one wants to get the number of excited electrons in the basis
    !% of the instantaneous eigenvalues of the Hamiltonian (Houston states).
    !%Option coordinates_sep 22
    !% Writes geometries in a separate file.
    !%Option velocities_sep 23
    !% Writes velocities in a separate file.
    !%Option forces_sep 24
    !% Writes forces in a separate file.
    !%Option total_heat_current 25
    !% Output the total heat current (average of the heat current density over the cell).
    !%Option total_magnetization 26
    !% Writes the total magnetization, where the total magnetization is calculated at the momentum
    !% defined by <tt>TDMomentumTransfer</tt>.
    !% This is used to extract the magnon frequency in case of a magnon kick.
    !%Option photons_q 27
    !% Writes photons_q in a separate file.
    !%Option maxwell_field 28
    !% Writes total electric field (if coupling is in length_geuge) or vector potential
    !% (if coupling is in velocity_gauge) coming from the interaction with Maxwell systems
    !% (only if Maxwell-TDDFT coupling is in dipole).
    !%Option norm_ks_orbitals 29
    !% Writes the norm of each Kohn-Sham orbital.
    !% The data is ordered per row as:
    !%  Iteration time (state 1 kpoint 1) (state2 kpoint1) ... (state-Nstates kpoint1) (state1 kpoint2) ... (state-Nstates kpoint-nkpt)
    !% noting that the kpoint index will also include the spin index for spin-polarised calculations.
    !%End
 
    !defaults
    output_options = .false.
    output_options(out_multipoles) = .true.
    output_options(out_energy) = .true.
    if (ions_move) then
      output_options(out_coords) = .true.
      output_options(out_temperature) = .true.
    end if
    if (with_gauge_field) output_options(out_gauge_field) = .true.
    if (list_has_lasers(ext_partners)) output_options(out_laser) = .true.
    if (kick%qkick_mode /= qkickmode_none) output_options(out_ftchd) = .true.
 
    call io_function_read_what_how_when(namespace, space, output_options, how, output_interval, &
      'TDOutput')
 
    ! Define which files to write
    do iout = 1, out_max
      writ%out(iout)%write = output_options(iout)
    end do
 
    ! experimental stuff
    if (writ%out(out_spin)%write) call messages_experimental('TDOutput = spin', namespace=namespace)
    if (writ%out(out_populations)%write) call messages_experimental('TDOutput = populations', namespace=namespace)
    if (writ%out(out_proj)%write) call messages_experimental('TDOutput = td_occup', namespace=namespace)
    if (writ%out(out_ion_ch)%write) call messages_experimental('TDOutput = ionization_channels', namespace=namespace)
    if (writ%out(out_partial_charges)%write) call messages_experimental('TDOutput = partial_charges', namespace=namespace)
    if (writ%out(out_kp_proj)%write) call messages_experimental('TDOutput = td_kpoint_occup', namespace=namespace)
    if (writ%out(out_floquet)%write) call messages_experimental('TDOutput = td_floquet', namespace=namespace)
    if (writ%out(out_n_ex)%write) call messages_experimental('TDOutput = n_excited_el', namespace=namespace)
    if (writ%out(out_q)%write) call messages_experimental('TDOutput = photons_q', namespace=namespace)
    if (writ%out(out_mxll_field)%write) call messages_experimental('TDOutput = maxwell_field', namespace=namespace)
 
    if (space%is_periodic() .and. writ%out(out_angular)%write) then
      call messages_not_implemented("TDOutput angular for periodic systems", namespace=namespace)
    end if
 
    !See comment in zstates_elec_mpdotp
    if (space%is_periodic() .and. writ%out(out_populations)%write) then
      call messages_not_implemented("TDOutput populations for periodic systems", namespace=namespace)
    end if
 
    if (writ%out(out_kp_proj)%write.or.writ%out(out_floquet)%write) then
      ! make sure this is not domain distributed
      if (gr%np /= gr%np_global) then
        message(1) = "TDOutput option td_kpoint_occup and td_floquet do not work with domain parallelization"
        call messages_fatal(1, namespace=namespace)
      end if
    end if
 
    if ((writ%out(out_separate_forces)%write .or. writ%out(out_coords)%write) .and. space%periodic_dim == 1) then
      call messages_input_error(namespace, 'TDOutput', &
        'Forces for systems periodic in 1D are not currently implemented and options that output the forces are not allowed.')
    end if
 
    ! NTD: The implementation of the option should be really redone properly
    if (writ%out(out_kp_proj)%write .and. hm%kpoints%nik_skip == 0) then
      message(1) = "TDOutput option td_kpoint_occup only work with zero-weight k-points at the moment."
      call messages_fatal(1, namespace=namespace)
    end if
 
    !%Variable TDOutputResolveStates
    !%Type logical
    !%Default No
    !%Section Time-Dependent::TD Output
    !%Description
    !% Defines whether the output should be resolved by states.
    !%
    !% So far only TDOutput = multipoles is supported.
    !%
    !%End
    call parse_variable(namespace, 'TDOutputResolveStates', .false., resolve_states)
    if (.not. writ%out(out_multipoles)%write) then
      write(message(1),'(a)') "TDOutputResolveStates works only for TDOutput = multipoles."
      call messages_warning(1, namespace=namespace)
    end if
 
    !%Variable TDMultipoleLmax
    !%Type integer
    !%Default 1
    !%Section Time-Dependent::TD Output
    !%Description
    !% Maximum electric multipole of the density output to the file <tt>td.general/multipoles</tt>
    !% during a time-dependent simulation. Must be non-negative.
    !%End
    call parse_variable(namespace, 'TDMultipoleLmax', 1, writ%lmax)
    if (writ%lmax < 0) then
      write(message(1), '(a,i6,a)') "Input: '", writ%lmax, "' is not a valid TDMultipoleLmax."
      message(2) = '(Must be TDMultipoleLmax >= 0 )'
      call messages_fatal(2, namespace=namespace)
    end if
    call messages_obsolete_variable(namespace, 'TDDipoleLmax', 'TDMultipoleLmax')
 
    ! Compatibility test
    if ((writ%out(out_acc)%write) .and. ions_move) then
      message(1) = 'If harmonic spectrum is to be calculated, atoms should not be allowed to move.'
      call messages_fatal(1, namespace=namespace)
    end if
 
    if ((writ%out(out_q)%write) .and. .not.(ks%has_photons)) then
      message(1) = 'If q(t) is to be calculated, you need to allow for photon modes.'
      call messages_fatal(1, namespace=namespace)
    end if
 
    if ((writ%out(out_mxll_field)%write) .and. .not. (hm%mxll%coupling_mode == velocity_gauge_dipole &
      .or. hm%mxll%add_electric_dip)) then
      message(1) = 'If the dipolar field has to be written, MaxwellCouplingMode has to be'
      message(2) = '"lenght_gauge_dipole" or "velocity_gauge_dipole" and at least one Maxwell system'
      message(3) = 'must be present.'
      call messages_fatal(3, namespace=namespace)
    end if
 
    rmin = ions%min_distance()
 
    ! This variable is documented in scf/scf.F90
    call parse_variable(namespace, 'LocalMagneticMomentsSphereRadius', min(m_half*rmin, lmm_r_single_atom), writ%lmm_r, &
      unit=units_inp%length)
 
    if (writ%out(out_proj)%write .or. writ%out(out_populations)%write &
      .or.writ%out(out_kp_proj)%write .or. writ%out(out_n_ex)%write) then
 
      if (st%parallel_in_states .and. writ%out(out_populations)%write) then
        message(1) = "Option TDOutput = populations is not implemented for parallel in states."
        call messages_fatal(1, namespace=namespace)
      end if
 
      if (writ%out(out_proj)%write .or. writ%out(out_populations)%write) then
        call states_elec_copy(writ%gs_st, st, exclude_wfns = .true., exclude_eigenval = .true.)
      else
        ! we want the same layout of gs_st as st
        call states_elec_copy(writ%gs_st, st, special=.true.)
      end if
 
      ! clean up all the stuff we have to reallocate
      safe_deallocate_a(writ%gs_st%node)
 
      call restart_init(restart_gs, namespace, restart_proj, restart_type_load, mc, ierr, mesh=gr)
 
      if (writ%out(out_proj)%write .or. writ%out(out_populations)%write) then
        if (ierr == 0) then
          call states_elec_look(restart_gs, ii, jj, writ%gs_st%nst, ierr)
        end if
        writ%gs_st%st_end = writ%gs_st%nst
        if (ierr /= 0) then
          message(1) = "Unable to read states information."
          call messages_fatal(1, namespace=namespace)
        end if
 
        writ%gs_st%st_start = 1
        ! do this only when not calculating populations, since all states are needed then
        if (.not. writ%out(out_populations)%write) then
          ! We will store the ground-state Kohn-Sham system for all processors.
          !%Variable TDProjStateStart
          !%Type integer
          !%Default 1
          !%Section Time-Dependent::TD Output
          !%Description
          !% To be used with <tt>TDOutput = td_occup</tt>. Not available if <tt>TDOutput = populations</tt>.
          !% Only output projections to states above <tt>TDProjStateStart</tt>. Usually one is only interested
          !% in particle-hole projections around the HOMO, so there is no need to calculate (and store)
          !% the projections of all TD states onto all static states. This sets a lower limit. The upper limit
          !% is set by the number of states in the propagation and the number of unoccupied states
          !% available.
          !%End
          call parse_variable(namespace, 'TDProjStateStart', 1, writ%gs_st%st_start)
 
          if (st%parallel_in_states .and. writ%out(out_proj)%write .and. writ%gs_st%st_start > 1) then
            call messages_not_implemented("TDOutput = td_occup for parallel in states with TDProjStateStart > 1", &
              namespace=namespace)
          end if
        end if
 
        writ%gs_st%lnst = writ%gs_st%st_end - writ%gs_st%st_start + 1
 
        call states_elec_deallocate_wfns(writ%gs_st)
 
        writ%gs_st%parallel_in_states = .false.
 
        ! allocate memory
        safe_allocate(writ%gs_st%occ(1:writ%gs_st%nst, 1:writ%gs_st%nik))
        safe_allocate(writ%gs_st%eigenval(1:writ%gs_st%nst, 1:writ%gs_st%nik))
 
        !We want all the task to have all the states
        !States can be distibuted for the states we propagate.
        safe_allocate(writ%gs_st%node(1:writ%gs_st%nst))
        writ%gs_st%node(:)  = 0
 
        writ%gs_st%eigenval = huge(writ%gs_st%eigenval)
        writ%gs_st%occ      = m_zero
        if (writ%gs_st%d%ispin == spinors) then
          safe_deallocate_a(writ%gs_st%spin)
          safe_allocate(writ%gs_st%spin(1:3, 1:writ%gs_st%nst, 1:writ%gs_st%nik))
        end if
 
        safe_allocate(skip(1:writ%gs_st%nst))
        skip = .false.
        skip(1:writ%gs_st%st_start-1) = .true.
 
        call states_elec_allocate_wfns(writ%gs_st, gr, type_cmplx, packed=.true., skip=skip)
 
        safe_deallocate_a(skip)
      end if
 
      call states_elec_load(restart_gs, namespace, space, writ%gs_st, gr, hm%kpoints, ierr, label = ': gs for TDOutput')
 
      if (ierr /= 0 .and. ierr /= (writ%gs_st%st_end-writ%gs_st%st_start+1)*writ%gs_st%nik &
        *writ%gs_st%d%dim*writ%gs_st%mpi_grp%size) then
        message(1) = "Unable to read wavefunctions for TDOutput."
        call messages_fatal(1, namespace=namespace)
      end if
      call restart_end(restart_gs)
    end if
 
    ! Build the excited states...
    if (writ%out(out_populations)%write) then
      !%Variable TDExcitedStatesToProject
      !%Type block
      !%Section Time-Dependent::TD Output
      !%Description
      !% <b>[WARNING: This is a *very* experimental feature]</b>
      !% To be used with <tt>TDOutput = populations</tt>.
      !% The population of the excited states
      !% (as defined by <Phi_I|Phi(t)> where |Phi(t)> is the many-body time-dependent state at
      !% time <i>t</i>, and |Phi_I> is the excited state of interest) can be approximated -- it is not clear
      !% how well -- by substituting for those real many-body states the time-dependent Kohn-Sham
      !% determinant and some modification of the Kohn-Sham ground-state determinant (<i>e.g.</i>,
      !% a simple HOMO-LUMO substitution, or the Casida ansatz for excited states in linear-response
      !% theory. If you set <tt>TDOutput</tt> to contain <tt>populations</tt>, you may ask for these approximated
      !% populations for a number of excited states, which will be described in the files specified
      !% in this block: each line should be the name of a file that contains one excited state.
      !%
      !% This file structure is the one written by the casida run mode, in the files in the directory <tt>*_excitations</tt>.
      !% The file describes the "promotions" from occupied
      !% to unoccupied levels that change the initial Slater determinant
      !% structure specified in ground_state. These promotions are a set
      !% of electron-hole pairs. Each line in the file, after an optional header, has four
      !% columns:
      !%
      !% <i>i  a  <math>\sigma</math> weight</i>
      !%
      !% where <i>i</i> should be an occupied state, <i>a</i> an unoccupied one, and <math>\sigma</math>
      !% the spin of the corresponding orbital. This pair is then associated with a
      !% creation-annihilation pair <math>a^{\dagger}_{a,\sigma} a_{i,\sigma}</math>, so that the
      !% excited state will be a linear combination in the form:
      !%
      !% <math>\left|{\rm ExcitedState}\right> =
      !% \sum weight(i,a,\sigma) a^{\dagger}_{a,\sigma} a_{i,\sigma} \left|{\rm GroundState}\right></math>
      !%
      !% where <i>weight</i> is the number in the fourth column.
      !% These weights should be normalized to one; otherwise the routine
      !% will normalize them, and write a warning.
      !%End
      if (parse_block(namespace, 'TDExcitedStatesToProject', blk) == 0) then
        writ%n_excited_states = parse_block_n(blk)
        safe_allocate(writ%excited_st(1:writ%n_excited_states))
        do ist = 1, writ%n_excited_states
          call parse_block_string(blk, ist-1, 0, filename)
          call excited_states_init(writ%excited_st(ist), writ%gs_st, trim(filename), namespace)
        end do
      else
        writ%n_excited_states = 0
      end if
    end if
 
    !%Variable TDOutputComputeInterval
    !%Type integer
    !%Default 50
    !%Section Time-Dependent::TD Output
    !%Description
    !% The TD output requested are computed
    !% when the iteration number is a multiple of the <tt>TDOutputComputeInterval</tt> variable.
    !% Must be >= 0. If it is 0, then no output is written.
    !% Implemented only for projections and number of excited electrons for the moment.
    !%End
    call parse_variable(namespace, 'TDOutputComputeInterval', 50, writ%compute_interval)
    if (writ%compute_interval < 0) then
      message(1) = "TDOutputComputeInterval must be >= 0."
      call messages_fatal(1, namespace=namespace)
    end if
 
    ! ----------------------------------------------
    ! Initialize write file handlers
    ! ----------------------------------------------
 
    ! Step
    if (iter == 0) then
      first = 0
    else
      first = iter + 1
    end if
 
    ! Root dir for TD files
    call io_mkdir('td.general', namespace)
 
    ! ----------------------------------------------
    ! Initialize write that are MPI agnostic
    ! ----------------------------------------------
 
    ! Default
    writ%out(:)%mpi_grp = mpi_world
 
    if (mpi_grp_is_root(mpi_world)) then
 
      do ifile = 1, out_max
        ! Exceptions that are handled later
        if (any([out_multipoles, out_laser, out_ftchd] == ifile)) cycle
 
        if (writ%out(ifile)%write) then
          call write_iter_init(writ%out(ifile)%handle, &
            first, &
            units_from_atomic(units_out%time, dt),  &
            trim(io_workpath("td.general/"//trim(td_file_name(ifile)), namespace)))
        end if
 
      enddo
 
      ! Exceptions
 
      ! Multipoles, when MPI is not state-parallel
      if (writ%out(out_multipoles)%write .and. .not. resolve_states) then
        call write_iter_init(writ%out(out_multipoles)%handle, &
          first, units_from_atomic(units_out%time, dt), &
          trim(io_workpath("td.general/multipoles", namespace)))
      end if
 
      if (writ%out(out_ftchd)%write) then
        select case (kick%qkick_mode)
        case (qkickmode_sin)
          write(filename, '(a)') 'td.general/ftchd.sin'
        case (qkickmode_cos)
          write(filename, '(a)') 'td.general/ftchd.cos'
        case (qkickmode_bessel)
          write(filename, '(a, SP, I0.3, a, I0.3)') 'td.general/ftchd.l', kick%qbessel_l, '_m', kick%qbessel_m
        case default
          write(filename, '(a)') 'td.general/ftchd'
        end select
        call write_iter_init(writ%out(out_ftchd)%handle, &
          first, units_from_atomic(units_out%time, dt), trim(io_workpath(filename, namespace)))
      end if
 
      if (writ%out(out_laser)%write) then
        if(associated(list_get_lasers(ext_partners))) then
          ! The laser file is written for the full propagation in one go, so that
          ! the user can check that the laser is correct and as intended before letting
          ! the code run for a possibly large period of time. This is done even after
          ! a restart, so that it takes into account any changes to max_iter.
          call io_rm("td.general/laser", namespace=namespace)
          call write_iter_init(writ%out(out_laser)%handle, 0, &
            units_from_atomic(units_out%time, dt),  &
            trim(io_workpath("td.general/laser", namespace)))
          do ii = 0, max_iter
            call td_write_laser(writ%out(out_laser)%handle, space, list_get_lasers(ext_partners), dt, ii)
          end do
          call write_iter_end(writ%out(out_laser)%handle)
        end if
      end if
 
    end if  ! mpi_grp_is_root(mpi_world)
 
    ! ----------------------------------------------
    ! Initialize write that are MPI group-specific
    ! ----------------------------------------------
 
    if (writ%out(out_multipoles)%write .and. resolve_states) then
      !resolve state contribution on multipoles
      writ%out(out_multipoles)%hand_start       = st%st_start
      writ%out(out_multipoles)%hand_end         = st%st_end
      writ%out(out_multipoles)%resolve_states   = .true.
      writ%out(out_multipoles)%mpi_grp          = gr%mpi_grp
 
      safe_allocate(writ%out(out_multipoles)%mult_handles(st%st_start:st%st_end))
 
      if (mpi_grp_is_root(writ%out(out_multipoles)%mpi_grp)) then
        do ist = st%st_start, st%st_end
          write(filename, '(a,i4.4)') 'td.general/multipoles-ist', ist
          call write_iter_init(writ%out(out_multipoles)%mult_handles(ist), &
            first, units_from_atomic(units_out%time, dt), &
            trim(io_workpath(filename, namespace)))
        end do
      end if
    end if
 
    ! Misc operations
 
    if (writ%out(out_total_current)%write .or. writ%out(out_total_heat_current)%write) then
      !TODO: we should only compute the current here, not v_ks
      call v_ks_calculate_current(ks, .true.)
      call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
        calc_eigenval=.false., time = iter*dt, calc_energy = .false.)
    end if
 
    if (writ%out(out_n_ex)%write .and. writ%compute_interval > 0) then
      call io_mkdir(outp%iter_dir, namespace)
    end if
 
    if (all(outp%how == 0) .and. writ%out(out_n_ex)%write) then
      call io_function_read_what_how_when(namespace, space, outp%what, outp%how, outp%output_interval)
    end if
 
    ! --------------------------------------------------------
    ! All steps of the above routine, but specific to DFT+U
    ! --------------------------------------------------------
 
    !%Variable TDOutputDFTU
    !%Type flag
    !%Default none
    !%Section Time-Dependent::TD Output
    !%Description
    !% Defines what should be output during the time-dependent
    !% simulation, related to DFT+U.
    !%
    !% Note: the output files generated by this option are updated
    !% every <tt>RestartWriteInterval</tt> steps.
    !%Option effective_u 1
    !% Writes the effective U for each orbital set as a function of time.
    !%End
    default = 0
    if (hm%lda_u_level == dft_u_acbn0) default = default + 2**(out_dftu_effective_u - 1)
    call parse_variable(namespace, 'TDOutputDFTU', default, flags)
 
    if (.not. varinfo_valid_option('TDOutputDFTU', flags, is_flag = .true.)) then
      call messages_input_error(namespace, 'TDOutputDFTU')
    end if
 
    do iout = 1, out_dftu_max
      writ%out_dftu(iout)%write = (iand(flags, 2**(iout - 1)) /= 0)
    end do
 
    if (mpi_grp_is_root(mpi_world)) then
      if (writ%out_dftu(out_dftu_effective_u)%write) then
        call write_iter_init(writ%out_dftu(out_dftu_effective_u)%handle, &
          first, units_from_atomic(units_out%time, dt),  &
          trim(io_workpath("td.general/effectiveU", namespace)))
      end if
    end if
 
    ! Init matrix elements calculator
    call writ%elec_me%init(gr, space, st)
 
    pop_sub(td_write_init)
  end subroutine td_write_init
 
 
  ! ---------------------------------------------------------
  subroutine td_write_end(writ)
    type(td_write_t), intent(inout) :: writ
 
    integer :: ist, iout
 
    push_sub(td_write_end)
 
    do iout = 1, out_max
      if (iout == out_laser) cycle
      if (writ%out(iout)%write) then
        if (mpi_grp_is_root(writ%out(iout)%mpi_grp)) then
          if (writ%out(iout)%resolve_states) then
            do ist = writ%out(iout)%hand_start, writ%out(iout)%hand_end
              call write_iter_end(writ%out(iout)%mult_handles(ist))
            end do
            safe_deallocate_a(writ%out(iout)%mult_handles)
          else
            call write_iter_end(writ%out(iout)%handle)
          end if
        end if
      end if
    end do
 
    if (mpi_grp_is_root(mpi_world)) then
      do iout = 1, out_dftu_max
        if (writ%out_dftu(iout)%write) call write_iter_end(writ%out_dftu(iout)%handle)
      end do
 
    end if
 
    if (writ%out(out_populations)%write) then
      do ist = 1, writ%n_excited_states
        call excited_states_kill(writ%excited_st(ist))
      end do
      writ%n_excited_states = 0
    end if
 
    if (writ%out(out_proj)%write .or. writ%out(out_populations)%write &
      .or. writ%out(out_n_ex)%write .or. writ%out(out_kp_proj)%write) then
      call states_elec_end(writ%gs_st)
    end if
 
    pop_sub(td_write_end)
  end subroutine td_write_end
 
 
  ! ---------------------------------------------------------
  subroutine td_write_iter(writ, namespace, space, outp, gr, st, hm, ions, ext_partners, kick, ks, dt, iter, mc, recalculate_gs)
    type(td_write_t),         intent(inout) :: writ
    type(namespace_t),        intent(in)    :: namespace
    class(space_t),           intent(in)    :: space
    type(output_t),           intent(in)    :: outp
    type(grid_t),             intent(in)    :: gr
    type(states_elec_t),      intent(inout) :: st
    type(hamiltonian_elec_t), intent(inout) :: hm
    type(ions_t),             intent(inout) :: ions
    type(partner_list_t),     intent(in)    :: ext_partners
    type(kick_t),             intent(in)    :: kick
    type(v_ks_t),             intent(in)    :: ks
    real(real64),             intent(in)    :: dt
    integer,                  intent(in)    :: iter
    type(multicomm_t),        intent(in)    :: mc
    logical,                  intent(in)    :: recalculate_gs
 
    type(gauge_field_t), pointer :: gfield
    integer :: ierr
    type(restart_t) :: restart_gs
 
    push_sub_with_profile(td_write_iter)
 
    if (writ%out(out_multipoles)%write) then
      call td_write_multipole(writ%out(out_multipoles), space, gr, ions, st, writ%lmax, kick, iter)
    end if
 
    if (writ%out(out_ftchd)%write) then
      call td_write_ftchd(writ%out(out_ftchd)%handle, space, gr, st, kick, iter)
    end if
 
    if (writ%out(out_angular)%write) then
      call td_write_angular(writ%out(out_angular)%handle, namespace, space, gr, ions, hm, st, kick, iter)
    end if
 
    if (writ%out(out_spin)%write) then
      call td_write_spin(writ%out(out_spin)%handle, gr, st, iter)
    end if
 
    if (writ%out(out_magnets)%write) then
      call td_write_local_magnetic_moments(writ%out(out_magnets)%handle, gr, st, ions, writ%lmm_r, iter)
    end if
 
    if (writ%out(out_tot_m)%write) then
      call td_write_tot_mag(writ%out(out_tot_m)%handle, gr, st, kick, iter)
    end if
 
    if (writ%out(out_proj)%write .and. mod(iter, writ%compute_interval) == 0) then
      if (mpi_grp_is_root(mpi_world)) call write_iter_set(writ%out(out_proj)%handle, iter)
      call td_write_proj(writ%out(out_proj)%handle, space, gr, ions, st, writ%gs_st, kick, iter)
    end if
 
    if (writ%out(out_floquet)%write) then
      call td_write_floquet(namespace, space, hm, ext_partners, gr, st, iter)
    end if
 
    if (writ%out(out_kp_proj)%write .and. mod(iter, writ%compute_interval) == 0) then
      call td_write_proj_kp(gr, hm%kpoints, st, writ%gs_st, namespace, iter)
    end if
 
    if (writ%out(out_coords)%write) then
      call td_write_coordinates(writ%out(out_coords)%handle, ions%natoms, ions%space, &
        ions%pos, ions%vel, ions%tot_force, iter)
    end if
 
    if (writ%out(out_separate_coords)%write) then
      call td_write_sep_coordinates(writ%out(out_separate_coords)%handle, ions%natoms, ions%space, &
        ions%pos, ions%vel, ions%tot_force, iter, 1)
    end if
 
    if (writ%out(out_separate_velocity)%write) then
      call td_write_sep_coordinates(writ%out(out_separate_velocity)%handle, ions%natoms, ions%space, &
        ions%pos, ions%vel, ions%tot_force, iter, 2)
    end if
 
    if (writ%out(out_separate_forces)%write) then
      call td_write_sep_coordinates(writ%out(out_separate_forces)%handle, ions%natoms, ions%space, &
        ions%pos, ions%vel, ions%tot_force, iter, 3)
    end if
 
    if (writ%out(out_temperature)%write) then
      call td_write_temperature(writ%out(out_temperature)%handle, ions, iter)
    end if
 
    if (writ%out(out_populations)%write) then
      call td_write_populations(writ%out(out_populations)%handle, namespace, space, gr, st, writ, dt, iter)
    end if
 
    if (writ%out(out_acc)%write) then
      call td_write_acc(writ%out(out_acc)%handle, namespace, space, gr, ions, st, hm, ext_partners, dt, iter)
    end if
 
    if (writ%out(out_vel)%write) then
      call td_write_vel(writ%out(out_vel)%handle, writ%elec_me, hm%kpoints, iter)
    end if
 
    ! td_write_laser no longer called here, because the whole laser is printed
    ! out at the beginning.
 
    if (writ%out(out_energy)%write) then
      call td_write_energy(writ%out(out_energy)%handle, hm, iter, ions%kinetic_energy)
    end if
 
    if (writ%out(out_gauge_field)%write) then
      gfield => list_get_gauge_field(ext_partners)
      if(associated(gfield)) then
        call gauge_field_output_write(gfield, writ%out(out_gauge_field)%handle, iter)
      end if
    end if
 
    if (writ%out(out_eigs)%write) then
      call td_write_eigs(writ%out(out_eigs)%handle, st, iter)
    end if
 
    if (writ%out(out_ion_ch)%write) then
      call td_write_ionch(writ%out(out_ion_ch)%handle, gr, st, iter)
    end if
 
    if (writ%out(out_total_current)%write) then
      call td_write_total_current(writ%out(out_total_current)%handle, space, gr, st, iter)
    end if
 
    if (writ%out(out_total_heat_current)%write) then
      call td_write_total_heat_current(writ%out(out_total_heat_current)%handle, space, hm, gr, st, iter)
    end if
 
    if (writ%out(out_partial_charges)%write) then
      call td_write_partial_charges(writ%out(out_partial_charges)%handle, gr, st, &
        ions, iter)
    end if
 
    if (writ%out(out_n_ex)%write .and. mod(iter, writ%compute_interval) == 0) then
      if (mpi_grp_is_root(mpi_world)) call write_iter_set(writ%out(out_n_ex)%handle, iter)
      if (recalculate_gs) then ! Load recomputed GS states
        call restart_init(restart_gs, namespace, restart_proj, restart_type_load, mc, ierr, mesh=gr)
        call states_elec_load(restart_gs, namespace, space, writ%gs_st, gr, hm%kpoints, &
          ierr, label = ': Houston states for TDOutput')
        call restart_end(restart_gs)
      end if
 
      call td_write_n_ex(writ%out(out_n_ex)%handle, outp, namespace, gr, hm%kpoints, st, writ%gs_st, iter)
    end if
 
    if (writ%out(out_norm_ks)%write) then
      call td_write_norm_ks_orbitals(writ%out(out_norm_ks)%handle, gr, hm%kpoints, st, iter)
    end if
 
    !DFT+U outputs
    if (writ%out_dftu(out_dftu_effective_u)%write) then
      call td_write_effective_u(writ%out_dftu(out_dftu_effective_u)%handle, hm%lda_u, iter)
    end if
 
    if (writ%out(out_q)%write .and. ks%has_photons) then
      call td_write_q(writ%out(out_q)%handle, space, ks, iter)
    end if
 
    if (writ%out(out_mxll_field)%write .and. hm%mxll%calc_field_dip) then
      call td_write_mxll_field(writ%out(out_mxll_field)%handle, space, hm, dt, iter)
    end if
 
    pop_sub_with_profile(td_write_iter)
  end subroutine td_write_iter
 
 
  ! ---------------------------------------------------------
  subroutine td_write_data(writ)
    type(td_write_t),     intent(inout) :: writ
 
    integer :: iout, ii
 
    push_sub(td_write_data)
    call profiling_in("TD_WRITE_DATA")
 
    do iout = 1, out_max
      if (iout == out_laser) cycle
      if (writ%out(iout)%write) then
        if (mpi_grp_is_root(writ%out(iout)%mpi_grp)) then
          if (writ%out(iout)%resolve_states) then
            do ii = writ%out(iout)%hand_start, writ%out(iout)%hand_end
              call write_iter_flush(writ%out(iout)%mult_handles(ii))
            end do
          else
            call write_iter_flush(writ%out(iout)%handle)
          end if
        end if
      end if
    end do
 
    if (mpi_grp_is_root(mpi_world)) then
      do iout = 1, out_dftu_max
        if (writ%out_dftu(iout)%write) call write_iter_flush(writ%out_dftu(iout)%handle)
      end do
    end if
 
    call profiling_out("TD_WRITE_DATA")
    pop_sub(td_write_data)
  end subroutine td_write_data
 
  ! ---------------------------------------------------------
  subroutine td_write_output(namespace, space, gr, st, hm, ks, outp, ions, ext_partners, iter, dt)
    type(namespace_t),        intent(in)    :: namespace
    type(electron_space_t),   intent(in)    :: space
    type(grid_t),             intent(in)    :: gr
    type(states_elec_t),      intent(inout) :: st
    type(hamiltonian_elec_t), intent(inout) :: hm
    type(v_ks_t),             intent(inout) :: ks
    type(output_t),           intent(in)    :: outp
    type(ions_t),             intent(in)    :: ions
    type(partner_list_t),     intent(in)    :: ext_partners
    integer,                  intent(in)    :: iter
    real(real64), optional,          intent(in)    :: dt
 
    character(len=256) :: filename
 
    push_sub(td_write_output)
    call profiling_in("TD_WRITE_OUTPUT")
 
    ! TODO this now overwrites wf inside st. If this is not wanted need to add an optional overwrite=no flag
    if (st%modelmbparticles%nparticle > 0) then
      call modelmb_sym_all_states(space, gr, st)
    end if
 
    ! now write down the rest
    write(filename, '(a,a,i7.7)') trim(outp%iter_dir),"td.", iter  ! name of directory
 
    call output_all(outp, namespace, space, filename, gr, ions, iter, st, hm, ks)
 
    call output_modelmb(outp, namespace, space, filename, gr, ions, iter, st)
 
    if (present(dt)) then
      call output_scalar_pot(outp, namespace, space, filename, gr, ions, ext_partners, iter*dt)
    else
      if (iter == 0) call output_scalar_pot(outp, namespace, space, filename, gr, ions, ext_partners)
    end if
 
    call profiling_out("TD_WRITE_OUTPUT")
    pop_sub(td_write_output)
  end subroutine td_write_output
 
  ! ---------------------------------------------------------
  subroutine td_write_spin(out_spin, mesh, st, iter)
    type(c_ptr),         intent(inout) :: out_spin
    class(mesh_t),       intent(in)    :: mesh
    type(states_elec_t), intent(in)    :: st
    integer,             intent(in)    :: iter
 
    character(len=130) :: aux
    real(real64) :: spin(3)
 
    push_sub(td_write_spin)
 
    ! The expectation value of the spin operator is half the total magnetic moment
    ! This has to be calculated by all nodes
    call magnetic_moment(mesh, st, st%rho, spin)
    spin = m_half*spin
 
    if (mpi_grp_is_root(mpi_world)) then ! only first node outputs
 
      if (iter == 0) then
        call td_write_print_header_init(out_spin)
 
        !second line -> columns name
        call write_iter_header_start(out_spin)
        if (st%d%ispin == spinors) then
          write(aux, '(a2,18x)') 'Sx'
          call write_iter_header(out_spin, aux)
          write(aux, '(a2,18x)') 'Sy'
          call write_iter_header(out_spin, aux)
        end if
        write(aux, '(a2,18x)') 'Sz'
        call write_iter_header(out_spin, aux)
        call write_iter_nl(out_spin)
 
        call td_write_print_header_end(out_spin)
      end if
 
      call write_iter_start(out_spin)
      select case (st%d%ispin)
      case (spin_polarized)
        call write_iter_double(out_spin, spin(3), 1)
      case (spinors)
        call write_iter_double(out_spin, spin(1:3), 3)
      end select
      call write_iter_nl(out_spin)
 
    end if
 
    pop_sub(td_write_spin)
  end subroutine td_write_spin
 
 
  ! ---------------------------------------------------------
  subroutine td_write_local_magnetic_moments(out_magnets, gr, st, ions, lmm_r, iter)
    type(c_ptr),              intent(inout) :: out_magnets
    type(grid_t),             intent(in)    :: gr
    type(states_elec_t),      intent(in)    :: st
    type(ions_t),             intent(in)    :: ions
    real(real64),             intent(in)    :: lmm_r
    integer,                  intent(in)    :: iter
 
    integer :: ia
    character(len=50) :: aux
    real(real64), allocatable :: lmm(:,:)
 
    push_sub(td_write_local_magnetic_moments)
 
    !get the atoms` magnetization. This has to be calculated by all nodes
    safe_allocate(lmm(1:3, 1:ions%natoms))
    call magnetic_local_moments(gr, st, ions, gr%der%boundaries, st%rho, lmm_r, lmm)
 
    if (mpi_grp_is_root(mpi_world)) then ! only first node outputs
 
      if (iter == 0) then
        call td_write_print_header_init(out_magnets)
 
        !second line -> columns name
        call write_iter_header_start(out_magnets)
        do ia = 1, ions%natoms
          if (st%d%ispin == spinors) then
            write(aux, '(a2,i2.2,16x)') 'mx', ia
            call write_iter_header(out_magnets, aux)
            write(aux, '(a2,i2.2,16x)') 'my', ia
            call write_iter_header(out_magnets, aux)
          end if
          write(aux, '(a2,i2.2,16x)') 'mz', ia
          call write_iter_header(out_magnets, aux)
        end do
        call write_iter_nl(out_magnets)
 
        call td_write_print_header_end(out_magnets)
      end if
 
      call write_iter_start(out_magnets)
      do ia = 1, ions%natoms
        select case (st%d%ispin)
        case (spin_polarized)
          call write_iter_double(out_magnets, lmm(3, ia), 1)
        case (spinors)
          call write_iter_double(out_magnets, lmm(1:3, ia), 3)
        end select
      end do
      call write_iter_nl(out_magnets)
      safe_deallocate_a(lmm)
    end if
 
    pop_sub(td_write_local_magnetic_moments)
  end subroutine td_write_local_magnetic_moments
 
  ! ---------------------------------------------------------
  subroutine td_write_tot_mag(out_magnets, mesh, st, kick, iter)
    type(c_ptr),              intent(inout) :: out_magnets
    class(mesh_t),            intent(in)    :: mesh
    type(states_elec_t),      intent(in)    :: st
    type(kick_t),             intent(in)    :: kick
    integer,                  intent(in)    :: iter
 
    complex(real64), allocatable :: tm(:,:)
    integer :: ii, iq
 
    push_sub(td_write_tot_mag)
 
    safe_allocate(tm(1:6,1:kick%nqvec))
 
    do iq = 1, kick%nqvec
      call magnetic_total_magnetization(mesh, st, kick%qvector(:,iq), tm(1:6,iq))
    end do
 
    if (mpi_grp_is_root(mpi_world)) then ! only first node outputs
 
      if (iter == 0) then
        call td_write_print_header_init(out_magnets)
        call kick_write(kick, out = out_magnets)
 
        !second line -> columns name
        call write_iter_header_start(out_magnets)
        call write_iter_header(out_magnets, 'Re[m_x(q)]')
        call write_iter_header(out_magnets, 'Im[m_x(q)]')
        call write_iter_header(out_magnets, 'Re[m_y(q)]')
        call write_iter_header(out_magnets, 'Im[m_y(q)]')
        call write_iter_header(out_magnets, 'Re[m_z(q)]')
        call write_iter_header(out_magnets, 'Im[m_z(q)]')
        call write_iter_header(out_magnets, 'Re[m_x(-q)]')
        call write_iter_header(out_magnets, 'Im[m_x(-q)]')
        call write_iter_header(out_magnets, 'Re[m_y(-q)]')
        call write_iter_header(out_magnets, 'Im[m_y(-q)]')
        call write_iter_header(out_magnets, 'Re[m_z(-q)]')
        call write_iter_header(out_magnets, 'Im[m_z(-q)]')
        call write_iter_nl(out_magnets)
 
        call td_write_print_header_end(out_magnets)
      end if
 
      call write_iter_start(out_magnets)
      do iq = 1, kick%nqvec
        do ii = 1, 6
          call write_iter_double(out_magnets, real(tm(ii, iq), real64), 1)
          call write_iter_double(out_magnets, aimag(tm(ii, iq)), 1)
        end do
      end do
      call write_iter_nl(out_magnets)
    end if
 
    safe_deallocate_a(tm)
 
    pop_sub(td_write_tot_mag)
  end subroutine td_write_tot_mag
 
 
  ! ---------------------------------------------------------
  subroutine td_write_angular(out_angular, namespace, space, gr, ions, hm, st, kick, iter)
    type(c_ptr),              intent(inout) :: out_angular
    type(namespace_t),        intent(in)    :: namespace
    class(space_t),           intent(in)    :: space
    type(grid_t),             intent(in)    :: gr
    type(ions_t),             intent(inout) :: ions
    type(hamiltonian_elec_t), intent(inout) :: hm
    type(states_elec_t),      intent(inout) :: st
    type(kick_t),             intent(in)    :: kick
    integer,                  intent(in)    :: iter
 
    integer :: idir
    character(len=130) :: aux
    real(real64) :: angular(3)
    class(perturbation_magnetic_t), pointer  :: angular_momentum
 
    push_sub(td_write_angular)
 
    angular_momentum => perturbation_magnetic_t(namespace, ions)
    do idir = 1, 3
      call angular_momentum%setup_dir(idir)
      !we have to multiply by 2, because the perturbation returns L/2
      angular(idir) = &
        m_two*real(angular_momentum%zstates_elec_expectation_value(namespace, space, gr, hm, st), real64)
    end do
    safe_deallocate_p(angular_momentum)
 
    if (mpi_grp_is_root(mpi_world)) then ! Only first node outputs
 
      if (iter == 0) then
        call td_write_print_header_init(out_angular)
 
        write(aux, '(a15,i2)')      '# nspin        ', st%d%nspin
        call write_iter_string(out_angular, aux)
        call write_iter_nl(out_angular)
 
        call kick_write(kick, out = out_angular)
 
        !second line -> columns name
        call write_iter_header_start(out_angular)
        write(aux, '(a4,18x)') '<Lx>'
        call write_iter_header(out_angular, aux)
        write(aux, '(a4,18x)') '<Ly>'
        call write_iter_header(out_angular, aux)
        write(aux, '(a4,18x)') '<Lz>'
        call write_iter_header(out_angular, aux)
        call write_iter_nl(out_angular)
 
        !third line -> should hold the units.
        call write_iter_string(out_angular, '#[Iter n.]')
        call write_iter_header(out_angular, '[' // trim(units_abbrev(units_out%time)) // ']')
        call write_iter_header(out_angular, '[hbar]')
        call write_iter_header(out_angular, '[hbar]')
        call write_iter_header(out_angular, '[hbar]')
        call write_iter_nl(out_angular)
 
        call td_write_print_header_end(out_angular)
      end if
 
      call write_iter_start(out_angular)
      call write_iter_double(out_angular, angular(1:3), 3)
      call write_iter_nl(out_angular)
 
    end if
 
    pop_sub(td_write_angular)
  end subroutine td_write_angular
 
 
  subroutine td_write_multipole(out_multip, space, gr, ions, st, lmax, kick, iter)
    type(td_write_prop_t),    intent(inout) :: out_multip
    class(space_t),           intent(in)    :: space
    type(grid_t),             intent(in)    :: gr
    type(ions_t),             intent(in)    :: ions
    type(states_elec_t),      intent(in)    :: st
    integer,                  intent(in)    :: lmax
    type(kick_t),             intent(in)    :: kick
    integer,                  intent(in)    :: iter
 
    integer :: ist
    real(real64), allocatable :: rho(:,:)
 
    push_sub(td_write_multipole)
 
    if (out_multip%resolve_states) then
      safe_allocate(rho(1:gr%np_part, 1:st%d%nspin))
      rho(:,:)     = m_zero
 
      do ist = st%st_start, st%st_end
        call density_calc(st, gr, rho, istin = ist)
        call td_write_multipole_r(out_multip%mult_handles(ist), space, gr, ions, st, lmax, kick, rho, iter, &
          mpi_grp = out_multip%mpi_grp)
      end do
 
      safe_deallocate_a(rho)
 
    else
      if (allocated(st%frozen_rho)) then
        safe_allocate(rho(1:gr%np, 1:st%d%nspin))
        call lalg_copy(gr%np, st%d%nspin, st%rho, rho)
        call lalg_axpy(gr%np, st%d%nspin, m_one, st%frozen_rho, rho)
 
        call td_write_multipole_r(out_multip%handle, space, gr, ions, st, lmax, kick, rho, iter)
 
        safe_deallocate_a(rho)
      else
        call td_write_multipole_r(out_multip%handle, space, gr, ions, st, lmax, kick, st%rho, iter)
      end if
 
    end if
 
    pop_sub(td_write_multipole)
  end subroutine td_write_multipole
 
 
  subroutine td_write_multipole_r(out_multip, space, mesh, ions, st, lmax, kick, rho, iter, mpi_grp)
    type(c_ptr),         intent(inout) :: out_multip
    class(space_t),      intent(in)    :: space
    class(mesh_t),       intent(in)    :: mesh
    type(ions_t),        intent(in)    :: ions
    type(states_elec_t), intent(in)    :: st
    integer,             intent(in)    :: lmax
    type(kick_t),        intent(in)    :: kick
    real(real64),        intent(in)    :: rho(:,:)
    integer,             intent(in)    :: iter
    type(mpi_grp_t), optional, intent(in) :: mpi_grp
 
 
    integer :: is, idir, ll, mm, add_lm
    character(len=120) :: aux
    real(real64) :: ionic_dipole(ions%space%dim)
    real(real64), allocatable :: multipole(:,:)
    type(mpi_grp_t)    :: mpi_grp_
 
    push_sub(td_write_multipole_r)
 
    ! We cannot output multipoles beyond the dipole for higher dimensions
    assert(.not. (lmax > 1 .and. space%dim > 3))
 
    mpi_grp_ = mpi_world
    if (present(mpi_grp)) mpi_grp_ = mpi_grp
 
    if (mpi_grp_is_root(mpi_grp_).and.iter == 0) then
      call td_write_print_header_init(out_multip)
 
      write(aux, '(a15,i2)')      '# nspin        ', st%d%nspin
      call write_iter_string(out_multip, aux)
      call write_iter_nl(out_multip)
 
      write(aux, '(a15,i2)')      '# lmax         ', lmax
      call write_iter_string(out_multip, aux)
      call write_iter_nl(out_multip)
 
      call kick_write(kick, out = out_multip)
 
      call write_iter_header_start(out_multip)
 
      do is = 1, st%d%nspin
        write(aux,'(a18,i1,a1)') 'Electronic charge(', is,')'
        call write_iter_header(out_multip, aux)
        if (lmax > 0) then
          do idir = 1, space%dim
            write(aux, '(4a1,i1,a1)') '<', index2axis(idir), '>', '(', is,')'
            call write_iter_header(out_multip, aux)
          end do
        end if
        do ll = 2, lmax
          do mm = -ll, ll
            write(aux, '(a2,i2,a4,i2,a2,i1,a1)') 'l=', ll, ', m=', mm, ' (', is,')'
            call write_iter_header(out_multip, aux)
          end do
        end do
      end do
      call write_iter_nl(out_multip)
 
      ! units
      call write_iter_string(out_multip, '#[Iter n.]')
      call write_iter_header(out_multip, '[' // trim(units_abbrev(units_out%time)) // ']')
 
      do is = 1, st%d%nspin
        call write_iter_header(out_multip, 'Electrons')
        if (lmax > 0) then
          do idir = 1, space%dim
            call write_iter_header(out_multip, '[' // trim(units_abbrev(units_out%length)) // ']')
          end do
        end if
        do ll = 2, lmax
          do mm = -ll, ll
            write(aux, '(a,a2,i1)') trim(units_abbrev(units_out%length)), "**", ll
            call write_iter_header(out_multip, '[' // trim(aux) // ']')
          end do
        end do
      end do
      call write_iter_nl(out_multip)
 
      call td_write_print_header_end(out_multip)
    end if
 
    if (space%dim > 3 .and. lmax == 1) then
      ! For higher dimensions we can only have up to the dipole
      safe_allocate(multipole(1:space%dim+1, 1:st%d%nspin))
    else
      safe_allocate(multipole(1:(lmax + 1)**2, 1:st%d%nspin))
    end if
    multipole(:,:) = m_zero
 
    do is = 1, st%d%nspin
      call dmf_multipoles(mesh, rho(:,is), lmax, multipole(:,is))
    end do
 
    if (lmax > 0) then
      ionic_dipole = ions%dipole()
      do is = 1, st%d%nspin
        multipole(2:space%dim+1, is) = -ionic_dipole(1:space%dim)/st%d%nspin - multipole(2:space%dim+1, is)
      end do
    end if
 
    if (mpi_grp_is_root(mpi_grp_)) then
      call write_iter_start(out_multip)
      do is = 1, st%d%nspin
        call write_iter_double(out_multip, units_from_atomic(units_out%length**0, multipole(1, is)), 1)
        if (lmax > 0) then
          do idir = 1, space%dim
            call write_iter_double(out_multip, units_from_atomic(units_out%length, multipole(1+idir, is)), 1)
          end do
        end if
        add_lm = space%dim + 2
        do ll = 2, lmax
          do mm = -ll, ll
            call write_iter_double(out_multip, units_from_atomic(units_out%length**ll, multipole(add_lm, is)), 1)
            add_lm = add_lm + 1
          end do
        end do
      end do
      call write_iter_nl(out_multip)
    end if
 
    safe_deallocate_a(multipole)
    pop_sub(td_write_multipole_r)
  end subroutine td_write_multipole_r
 
  ! ---------------------------------------------------------
  subroutine td_write_ftchd(out_ftchd, space, mesh, st, kick, iter)
    type(c_ptr),         intent(inout) :: out_ftchd
    class(space_t),      intent(in)    :: space
    class(mesh_t),       intent(in)    :: mesh
    type(states_elec_t), intent(in)    :: st
    type(kick_t),        intent(in)    :: kick
    integer,             intent(in)    :: iter
 
    integer :: is, ip, idir
    character(len=120) :: aux, aux2
    real(real64)   :: ftchd_bessel
    complex(real64)   :: ftchd
    real(real64)   :: ylm
    real(real64), allocatable :: integrand_bessel(:)
    complex(real64), allocatable :: integrand(:)
 
    push_sub(td_write_ftchd)
 
    if (mpi_grp_is_root(mpi_world).and.iter == 0) then
      call td_write_print_header_init(out_ftchd)
 
      write(aux,'(a15, i2)') '# qkickmode    ', kick%qkick_mode
      call write_iter_string(out_ftchd, aux)
      call write_iter_nl(out_ftchd)
 
      if (kick%qkick_mode == qkickmode_bessel) then
        write(aux,'(a15, i0.3, 1x, i0.3)') '# ll, mm       ', kick%qbessel_l, kick%qbessel_m
        call write_iter_string(out_ftchd, aux)
        call write_iter_nl(out_ftchd)
      end if
 
      if (kick%qkick_mode == qkickmode_bessel) then
        write(aux, '(a15, f9.6)') '# qlength      ', kick%qlength
      else ! sin or cos
        write(aux, '(a15)')       '# qvector      '
        do idir = 1, space%dim
          write(aux2, '(f9.5)') kick%qvector(idir,1)
          aux = trim(aux) // trim(aux2)
        end do
      end if
      call write_iter_string(out_ftchd, aux)
      call write_iter_nl(out_ftchd)
 
      write(aux, '(a15,f18.12)')  '# kick strength', kick%delta_strength
      call write_iter_string(out_ftchd, aux)
      call write_iter_nl(out_ftchd)
 
      call write_iter_header_start(out_ftchd)
      if (kick%qkick_mode == qkickmode_bessel) then
        write(aux,'(a17)') 'int(j_l*Y_lm*rho)'
      else
        write(aux,'(a12)') 'Real, Imag'
      end if
      call write_iter_header(out_ftchd, aux)
      call write_iter_nl(out_ftchd)
 
      ! units
      call write_iter_string(out_ftchd, '#[Iter n.]')
      call write_iter_header(out_ftchd, '[' // trim(units_abbrev(units_out%time)) // ']')
      call write_iter_nl(out_ftchd)
      call td_write_print_header_end(out_ftchd)
 
    end if
 
    ftchd = m_zero
 
    ! If kick mode is exp, sin, or cos, apply the normal Fourier transform
    if (kick%qkick_mode /= qkickmode_bessel) then
      safe_allocate(integrand(1:mesh%np))
      integrand = m_zero
      do is = 1, st%d%nspin
        do ip = 1, mesh%np
          integrand(ip) = integrand(ip) + st%rho(ip, is) * exp(-m_zi*sum(mesh%x(ip, 1:space%dim)*kick%qvector(1:space%dim, 1)))
        end do
      end do
      ftchd = zmf_integrate(mesh, integrand)
      safe_deallocate_a(integrand)
    else
      ftchd_bessel = m_zero
      safe_allocate(integrand_bessel(1:mesh%np))
      integrand_bessel = m_zero
      do is = 1, st%d%nspin
        do ip = 1, mesh%np
          call ylmr_real(mesh%x(ip, 1:3), kick%qbessel_l, kick%qbessel_m, ylm)
          integrand_bessel(ip) = integrand_bessel(ip) + st%rho(ip, is) * &
            loct_sph_bessel(kick%qbessel_l, kick%qlength*norm2(mesh%x(ip, :)))*ylm
        end do
      end do
      ftchd_bessel = dmf_integrate(mesh, integrand_bessel)
      safe_deallocate_a(integrand_bessel)
    end if
 
    if (mpi_grp_is_root(mpi_world)) then
      call write_iter_start(out_ftchd)
      if (kick%qkick_mode == qkickmode_bessel) then
        call write_iter_double(out_ftchd, ftchd_bessel, 1)
      else ! exp, sin, cos
        call write_iter_double(out_ftchd, real(ftchd), 1)
        call write_iter_double(out_ftchd, aimag(ftchd), 1)
      end if
      call write_iter_nl(out_ftchd)
    end if
 
    pop_sub(td_write_ftchd)
  end subroutine td_write_ftchd
 
  ! ---------------------------------------------------------
  subroutine td_write_temperature(out_temperature, ions, iter)
    type(c_ptr),       intent(inout) :: out_temperature
    type(ions_t),      intent(in) :: ions
    integer,           intent(in) :: iter
 
    if (.not. mpi_grp_is_root(mpi_world)) return ! only first node outputs
 
    push_sub(td_write_temperature)
 
    if (iter == 0) then
      call td_write_print_header_init(out_temperature)
 
      ! first line: column names
      call write_iter_header_start(out_temperature)
      call write_iter_header(out_temperature, 'Temperature')
      call write_iter_nl(out_temperature)
 
      ! second line: units
      call write_iter_string(out_temperature, '#[Iter n.]')
      call write_iter_header(out_temperature, '[' // trim(units_abbrev(units_out%time)) // ']')
      call write_iter_string(out_temperature, '        [K]')
      call write_iter_nl(out_temperature)
 
      call td_write_print_header_end(out_temperature)
    end if
 
    call write_iter_start(out_temperature)
 
    call write_iter_double(out_temperature, units_from_atomic(unit_kelvin, ion_dynamics_temperature(ions)), 1)
 
    call write_iter_nl(out_temperature)
 
    pop_sub(td_write_temperature)
  end subroutine td_write_temperature
 
 
  ! ---------------------------------------------------------
  subroutine td_write_populations(out_populations, namespace, space, mesh, st, writ, dt, iter)
    type(c_ptr),            intent(inout) :: out_populations
    type(namespace_t),      intent(in)    :: namespace
    class(space_t),         intent(in)    :: space
    class(mesh_t),          intent(in)    :: mesh
    type(states_elec_t),    intent(inout) :: st
    type(td_write_t),       intent(in)    :: writ
    real(real64),           intent(in)    :: dt
    integer,                intent(in)    :: iter
 
    integer :: ist
    character(len=6) :: excited_name
    complex(real64) :: gsp
    complex(real64), allocatable :: excited_state_p(:)
    complex(real64), allocatable :: dotprodmatrix(:, :, :)
 
 
    push_sub(td_write_populations)
 
    safe_allocate(dotprodmatrix(1:writ%gs_st%nst, 1:st%nst, 1:st%nik))
    call zstates_elec_matrix(writ%gs_st, st, mesh, dotprodmatrix)
 
 
    !See comment in zstates_elec_mpdotp
    assert(.not. space%is_periodic())
 
    ! all processors calculate the projection
    gsp = zstates_elec_mpdotp(namespace, mesh, writ%gs_st, st, dotprodmatrix)
 
    if (writ%n_excited_states > 0) then
      safe_allocate(excited_state_p(1:writ%n_excited_states))
      do ist = 1, writ%n_excited_states
        excited_state_p(ist) = zstates_elec_mpdotp(namespace, mesh, writ%excited_st(ist), st, dotprodmatrix)
      end do
    end if
 
    if (mpi_grp_is_root(mpi_world)) then
      if (iter == 0) then
        call td_write_print_header_init(out_populations)
 
        ! first line -> column names
        call write_iter_header_start(out_populations)
        call write_iter_header(out_populations, 'Re<Phi_gs|Phi(t)>')
        call write_iter_header(out_populations, 'Im<Phi_gs|Phi(t)>')
        do ist = 1, writ%n_excited_states
          write(excited_name,'(a2,i3,a1)') 'P(', ist,')'
          call write_iter_header(out_populations, 'Re<'//excited_name//'|Phi(t)>')
          call write_iter_header(out_populations, 'Im<'//excited_name//'|Phi(t)>')
        end do
        call write_iter_nl(out_populations)
 
        ! second line -> units
        call write_iter_string(out_populations, '#[Iter n.]')
        call write_iter_header(out_populations, '[' // trim(units_abbrev(units_out%time)) // ']')
        call write_iter_nl(out_populations)
 
        call td_write_print_header_end(out_populations)
      end if
 
      ! cannot call write_iter_start, for the step is not 1
      call write_iter_int(out_populations, iter, 1)
      call write_iter_double(out_populations, units_from_atomic(units_out%time, iter*dt),  1)
      call write_iter_double(out_populations, real(gsp),  1)
      call write_iter_double(out_populations, aimag(gsp), 1)
      do ist = 1, writ%n_excited_states
        call write_iter_double(out_populations, real(excited_state_p(ist)),  1)
        call write_iter_double(out_populations, aimag(excited_state_p(ist)), 1)
      end do
      call write_iter_nl(out_populations)
    end if
 
    if (writ%n_excited_states > 0) then
      safe_deallocate_a(excited_state_p)
    end if
    safe_deallocate_a(dotprodmatrix)
    pop_sub(td_write_populations)
  end subroutine td_write_populations
 
 
  ! ---------------------------------------------------------
  subroutine td_write_acc(out_acc, namespace, space, gr, ions, st, hm, ext_partners, dt, iter)
    type(c_ptr),              intent(inout) :: out_acc
    type(namespace_t),        intent(in)    :: namespace
    class(space_t),           intent(in)    :: space
    type(grid_t),             intent(in)    :: gr
    type(ions_t),             intent(inout) :: ions
    type(states_elec_t),      intent(inout) :: st
    type(hamiltonian_elec_t), intent(inout) :: hm
    type(partner_list_t),     intent(in)    :: ext_partners
    real(real64),             intent(in)    :: dt
    integer,                  intent(in)    :: iter
 
    integer :: idim
    character(len=7) :: aux
    real(real64) :: acc(space%dim)
 
    push_sub(td_write_acc)
 
    if (iter == 0 .and. mpi_grp_is_root(mpi_world)) then
      call td_write_print_header_init(out_acc)
 
      ! first line -> column names
      call write_iter_header_start(out_acc)
      do idim = 1, space%dim
        write(aux, '(a4,i1,a1)') 'Acc(', idim, ')'
        call write_iter_header(out_acc, aux)
      end do
      call write_iter_nl(out_acc)
 
      ! second line: units
      call write_iter_string(out_acc, '#[Iter n.]')
      call write_iter_header(out_acc, '[' // trim(units_abbrev(units_out%time)) // ']')
      do idim = 1, space%dim
        call write_iter_header(out_acc, '[' // trim(units_abbrev(units_out%acceleration)) // ']')
      end do
      call write_iter_nl(out_acc)
      call td_write_print_header_end(out_acc)
    end if
 
    call td_calc_tacc(namespace, space, gr, ions, ext_partners, st, hm, acc, dt*iter)
 
    if (mpi_grp_is_root(mpi_world)) then
      call write_iter_start(out_acc)
      acc = units_from_atomic(units_out%acceleration, acc)
      call write_iter_double(out_acc, acc, space%dim)
      call write_iter_nl(out_acc)
    end if
 
    pop_sub(td_write_acc)
  end subroutine td_write_acc
 
  ! ---------------------------------------------------------
  subroutine td_write_vel(out_vel, elec_me, kpoints, iter)
    type(c_ptr),                  intent(inout) :: out_vel
    type(elec_matrix_elements_t), intent(in)    :: elec_me
    type(kpoints_t),              intent(in)    :: kpoints
    integer,                      intent(in)    :: iter
 
    integer :: idim
    character(len=7) :: aux
    real(real64) :: vel(elec_me%space%dim)
 
    push_sub(td_write_vel)
 
    if (iter == 0 .and. mpi_grp_is_root(mpi_world)) then
      call td_write_print_header_init(out_vel)
 
      ! first line -> column names
      call write_iter_header_start(out_vel)
      do idim = 1, elec_me%space%dim
        write(aux, '(a4,i1,a1)') 'Vel(', idim, ')'
        call write_iter_header(out_vel, aux)
      end do
      call write_iter_nl(out_vel)
 
      ! second line: units
      call write_iter_string(out_vel, '#[Iter n.]')
      call write_iter_header(out_vel, '[' // trim(units_abbrev(units_out%time)) // ']')
      do idim = 1, elec_me%space%dim
        call write_iter_header(out_vel, '[' // trim(units_abbrev(units_out%velocity)) // ']')
      end do
      call write_iter_nl(out_vel)
      call td_write_print_header_end(out_vel)
    end if
 
    call td_calc_tvel(elec_me, kpoints, vel)
 
    if (mpi_grp_is_root(mpi_world)) then
      call write_iter_start(out_vel)
      vel = units_from_atomic(units_out%velocity, vel)
      call write_iter_double(out_vel, vel, elec_me%space%dim)
      call write_iter_nl(out_vel)
    end if
 
    pop_sub(td_write_vel)
  end subroutine td_write_vel
 
 
  ! ---------------------------------------------------------
  subroutine td_write_laser(out_laser, space, lasers, dt, iter)
    type(c_ptr),              intent(inout) :: out_laser
    class(space_t),           intent(in)    :: space
    type(lasers_t),           intent(inout)    :: lasers
    real(real64),             intent(in)    :: dt
    integer,                  intent(in)    :: iter
 
    integer :: il, idir
    real(real64) :: field(space%dim)
    real(real64) :: ndfield(space%dim)
    character(len=80) :: aux
 
    if (.not. mpi_grp_is_root(mpi_world)) return ! only first node outputs
 
    ! no PUSH SUB, called too often
 
    if (iter == 0) then
      call td_write_print_header_init(out_laser)
 
      ! first line
      write(aux, '(a7,e20.12,3a)') '# dt = ', units_from_atomic(units_out%time, dt), &
        " [", trim(units_abbrev(units_out%time)), "]"
      call write_iter_string(out_laser, aux)
      call write_iter_nl(out_laser)
 
      call write_iter_header_start(out_laser)
      do il = 1, lasers%no_lasers
        select case (laser_kind(lasers%lasers(il)))
        case (e_field_electric)
          do idir = 1, space%dim
            write(aux, '(a,i1,a)') 'E(', idir, ')'
            call write_iter_header(out_laser, aux)
          end do
        case (e_field_magnetic)
          do idir = 1, space%dim
            write(aux, '(a,i1,a)') 'B(', idir, ')'
            call write_iter_header(out_laser, aux)
          end do
        case (e_field_vector_potential)
          do idir = 1, space%dim
            write(aux, '(a,i1,a)') 'A(', idir, ')'
            call write_iter_header(out_laser, aux)
          end do
        case (e_field_scalar_potential)
          write(aux, '(a,i1,a)') 'e(t)'
          call write_iter_header(out_laser, aux)
        end select
      end do
 
      if (lasers_with_nondipole_field(lasers)) then
        do idir = 1, space%dim
          write(aux, '(a,i1,a)') 'A^M(', idir, ')'
          call write_iter_header(out_laser, aux)
        end do
      end if
      call write_iter_nl(out_laser)
 
      call write_iter_string(out_laser, '#[Iter n.]')
      call write_iter_header(out_laser, '[' // trim(units_abbrev(units_out%time)) // ']')
 
      ! Note that we do not print out units of E, B, or A, but rather units of e*E, e*B, e*A.
      ! (force, force, and energy, respectively). The reason is that the units of E, B or A
      ! are ugly.
      do il = 1, lasers%no_lasers
        select case (laser_kind(lasers%lasers(il)))
        case (e_field_electric, e_field_magnetic)
          aux = '[' // trim(units_abbrev(units_out%force)) // ']'
          do idir = 1, space%dim
            call write_iter_header(out_laser, aux)
          end do
        case (e_field_vector_potential)
          aux = '[' // trim(units_abbrev(units_out%energy)) // ']'
          do idir = 1, space%dim
            call write_iter_header(out_laser, aux)
          end do
        case (e_field_scalar_potential)
          aux = '[adim]'
          call write_iter_header(out_laser, aux)
        end select
      end do
 
      if (lasers_with_nondipole_field(lasers)) then
        aux = '[' // trim(units_abbrev(units_out%energy)) // ']'
        do idir = 1, space%dim
          call write_iter_header(out_laser, aux)
        end do
      end if
 
 
      call write_iter_nl(out_laser)
 
      call td_write_print_header_end(out_laser)
    end if
 
    call write_iter_start(out_laser)
 
    do il = 1, lasers%no_lasers
      field = m_zero
      call laser_field(lasers%lasers(il), field(1:space%dim), iter*dt)
      select case (laser_kind(lasers%lasers(il)))
      case (e_field_electric, e_field_magnetic)
        field = units_from_atomic(units_out%force, field)
        call write_iter_double(out_laser, field, space%dim)
      case (e_field_vector_potential)
        field = units_from_atomic(units_out%energy, field)
        call write_iter_double(out_laser, field, space%dim)
      case (e_field_scalar_potential)
        call write_iter_double(out_laser, field(1), 1)
      end select
    end do
 
    if (lasers_with_nondipole_field(lasers)) then
      call lasers_nondipole_laser_field_step(lasers, ndfield, iter*dt)
      call lasers_set_nondipole_parameters(lasers, ndfield, iter*dt)
      call write_iter_double(out_laser, ndfield, space%dim)
    end if
    call write_iter_nl(out_laser)
 
  end subroutine td_write_laser
 
 
  ! ---------------------------------------------------------
  subroutine td_write_energy(out_energy, hm, iter, ke)
    type(c_ptr),         intent(inout)   :: out_energy
    type(hamiltonian_elec_t), intent(in) :: hm
    integer,             intent(in)      :: iter
    real(real64), intent(in)                    :: ke
 
    integer :: ii
 
    integer :: n_columns
 
    if (.not. mpi_grp_is_root(mpi_world)) return ! only first node outputs
 
    push_sub(td_write_energy)
 
    n_columns = 9
 
    if (iter == 0) then
      call td_write_print_header_init(out_energy)
 
      ! first line -> column names
      call write_iter_header_start(out_energy)
      call write_iter_header(out_energy, 'Total')
      call write_iter_header(out_energy, 'Kinetic (ions)')
      call write_iter_header(out_energy, 'Ion-Ion')
      call write_iter_header(out_energy, 'Electronic')
      call write_iter_header(out_energy, 'Eigenvalues')
      call write_iter_header(out_energy, 'Hartree')
      call write_iter_header(out_energy, 'Int[n v_xc]')
      call write_iter_header(out_energy, 'Exchange')
      call write_iter_header(out_energy, 'Correlation')
 
      if (hm%pcm%run_pcm) then
        call write_iter_header(out_energy, 'E_M-solvent')
        n_columns = n_columns + 1
      end if
 
      if (hm%lda_u_level /= dft_u_none) then
        call write_iter_header(out_energy, 'Hubbard')
        n_columns = n_columns + 1
      end if
 
      call write_iter_nl(out_energy)
 
      ! units
 
      call write_iter_string(out_energy, '#[Iter n.]')
      call write_iter_header(out_energy, '[' // trim(units_abbrev(units_out%time)) // ']')
 
      do ii = 1, n_columns
        call write_iter_header(out_energy, '[' // trim(units_abbrev(units_out%energy)) // ']')
      end do
      call write_iter_nl(out_energy)
 
 
      call td_write_print_header_end(out_energy)
    end if
 
    call write_iter_start(out_energy)
    call write_iter_double(out_energy, units_from_atomic(units_out%energy, hm%energy%total+ke), 1)
    call write_iter_double(out_energy, units_from_atomic(units_out%energy, ke), 1)
    call write_iter_double(out_energy, units_from_atomic(units_out%energy, hm%ep%eii), 1)
    call write_iter_double(out_energy, units_from_atomic(units_out%energy, hm%energy%total-hm%ep%eii), 1)
    call write_iter_double(out_energy, units_from_atomic(units_out%energy, hm%energy%eigenvalues), 1)
    call write_iter_double(out_energy, units_from_atomic(units_out%energy, hm%energy%hartree), 1)
    call write_iter_double(out_energy, units_from_atomic(units_out%energy, hm%energy%intnvxc), 1)
    call write_iter_double(out_energy, units_from_atomic(units_out%energy, hm%energy%exchange), 1)
    call write_iter_double(out_energy, units_from_atomic(units_out%energy, hm%energy%correlation), 1)
 
    !adding the molecule-solvent electrostatic interaction
    if (hm%pcm%run_pcm) call write_iter_double(out_energy, &
      units_from_atomic(units_out%energy, hm%energy%int_ee_pcm + hm%energy%int_en_pcm + &
      hm%energy%int_nn_pcm + hm%energy%int_ne_pcm), 1)
 
    if (hm%lda_u_level /= dft_u_none) then
      call write_iter_double(out_energy, units_from_atomic(units_out%energy, hm%energy%dft_u), 1)
    end if
 
    call write_iter_nl(out_energy)
 
    pop_sub(td_write_energy)
  end subroutine td_write_energy
 
  ! ---------------------------------------------------------
  subroutine td_write_eigs(out_eigs, st, iter)
    type(c_ptr),         intent(inout) :: out_eigs
    type(states_elec_t), intent(in)    :: st
    integer,             intent(in)    :: iter
 
    integer             :: ii, is
    character(len=68)   :: buf
 
    push_sub(td_write_eigs)
 
    if (.not. mpi_grp_is_root(mpi_world)) then
      pop_sub(td_write_eigs)
      return ! only first node outputs
    end if
 
 
    if (iter == 0) then
      call td_write_print_header_init(out_eigs)
 
      write(buf, '(a15,i2)')      '# nst          ', st%nst
      call write_iter_string(out_eigs, buf)
      call write_iter_nl(out_eigs)
 
      write(buf, '(a15,i2)')      '# nspin        ', st%d%nspin
      call write_iter_string(out_eigs, buf)
      call write_iter_nl(out_eigs)
 
      ! first line -> column names
      call write_iter_header_start(out_eigs)
      do is = 1, st%d%kpt%nglobal
        do ii = 1, st%nst
          write(buf, '(a,i4)') 'Eigenvalue ',ii
          call write_iter_header(out_eigs, buf)
        end do
      end do
      call write_iter_nl(out_eigs)
 
      ! second line: units
      call write_iter_string(out_eigs, '#[Iter n.]')
      call write_iter_header(out_eigs, '[' // trim(units_abbrev(units_out%time)) // ']')
      do is = 1, st%d%kpt%nglobal
        do ii = 1, st%nst
          call write_iter_header(out_eigs, '[' // trim(units_abbrev(units_out%energy)) // ']')
        end do
      end do
      call write_iter_nl(out_eigs)
      call td_write_print_header_end(out_eigs)
    end if
 
    call write_iter_start(out_eigs)
    do is = 1, st%d%kpt%nglobal
      do ii =1 , st%nst
        call write_iter_double(out_eigs, units_from_atomic(units_out%energy, st%eigenval(ii,is)), 1)
      end do
    end do
    call write_iter_nl(out_eigs)
 
    pop_sub(td_write_eigs)
  end subroutine td_write_eigs
 
  ! ---------------------------------------------------------
  subroutine td_write_ionch(out_ionch, mesh, st, iter)
    type(c_ptr),         intent(inout) :: out_ionch
    class(mesh_t),       intent(in)    :: mesh
    type(states_elec_t), intent(in)    :: st
    integer,             intent(in)    :: iter
 
    integer             :: ii, ist, Nch, ik, idim
    character(len=68)   :: buf
    real(real64), allocatable  :: ch(:), occ(:)
    real(real64), allocatable :: occbuf(:)
 
    push_sub(td_write_ionch)
 
 
    nch =  st%nst * st%d%kpt%nglobal * st%d%dim
    safe_allocate(ch(0: nch))
    safe_allocate(occ(0: nch))
 
    occ(:) = m_zero
    ii = 1
    do ik = 1, st%nik
      do ist = 1, st%nst
        do idim = 1, st%d%dim
          if (st%st_start <= ist .and. ist <= st%st_end .and. &
            st%d%kpt%start <= ik .and. ik <= st%d%kpt%end) then
            occ(ii) = st%occ(ist, ik)
          end if
          ii = ii+1
        end do
      end do
    end do
 
 
    if (st%parallel_in_states) then
      safe_allocate(occbuf(0: nch))
      occbuf(:) = m_zero
      call st%mpi_grp%allreduce(occ(0), occbuf(0), nch+1, mpi_double_precision, mpi_sum)
      occ(:) = occbuf(:)
      safe_deallocate_a(occbuf)
    end if
 
    !Calculate the channels
    call td_calc_ionch(mesh, st, ch, nch)
 
 
    if (.not. mpi_grp_is_root(mpi_world)) then
      safe_deallocate_a(ch)
      pop_sub(td_write_ionch)
      return ! only first node outputs
    end if
 
 
    if (iter == 0) then
      call td_write_print_header_init(out_ionch)
 
      ! first line -> column names
      call write_iter_header_start(out_ionch)
 
      do ii = 0, nch
        if (occ(ii)>m_zero .or. ii == 0) then
          write(buf, '(a,f4.1,a)') 'Pion(',occ(ii)*ii,'+, t)'
          call write_iter_header(out_ionch, buf)
        end if
      end do
      call write_iter_nl(out_ionch)
 
      ! second line: units
      call write_iter_string(out_ionch, '#[Iter n.]')
      call write_iter_header(out_ionch, '[' // trim(units_abbrev(units_out%time)) // ']')
      do ii = 0, nch
        if (occ(ii)>m_zero .or. ii == 0) then
          call write_iter_header(out_ionch, '[' // trim(units_abbrev(unit_one)) // ']')
        end if
      end do
      call write_iter_nl(out_ionch)
      call td_write_print_header_end(out_ionch)
    end if
 
    call write_iter_start(out_ionch)
    do ii =0 , nch
      if (occ(ii)>m_zero .or. ii == 0) then
        call write_iter_double(out_ionch, units_from_atomic(unit_one, ch(ii)), 1)
      end if
    end do
    call write_iter_nl(out_ionch)
 
    safe_deallocate_a(ch)
    safe_deallocate_a(occ)
 
    pop_sub(td_write_ionch)
  end subroutine td_write_ionch
 
  ! ---------------------------------------------------------
  subroutine td_write_proj(out_proj, space, mesh, ions, st, gs_st, kick, iter)
    type(c_ptr),         intent(inout) :: out_proj
    class(space_t),      intent(in)    :: space
    class(mesh_t),       intent(in)    :: mesh
    type(ions_t),        intent(in)    :: ions
    type(states_elec_t), intent(inout) :: st
    type(states_elec_t), intent(in)    :: gs_st
    type(kick_t),        intent(in)    :: kick
    integer,             intent(in)    :: iter
 
    complex(real64), allocatable :: projections(:,:,:)
    character(len=80) :: aux
    integer :: ik, ist, uist, idir
 
    push_sub(td_write_proj)
 
    if (iter == 0) then
      if (mpi_grp_is_root(mpi_world)) then
        call td_write_print_header_init(out_proj)
 
        write(aux, '(a15,i2)')      '# nspin        ', st%d%nspin
        call write_iter_string(out_proj, aux)
        call write_iter_nl(out_proj)
 
        call kick_write(kick, out = out_proj)
 
        call write_iter_string(out_proj, "#%")
        call write_iter_nl(out_proj)
 
        write(aux, '(a,i8)') "# nik  ", st%nik
        call write_iter_string(out_proj, aux)
        call write_iter_nl(out_proj)
 
        write(aux, '(a,2i8)') "#  st  ", gs_st%st_start, st%nst
        call write_iter_string(out_proj, aux)
        call write_iter_nl(out_proj)
 
        write(aux, '(a,2i8)') "# ust  ", gs_st%st_start, gs_st%st_end
        call write_iter_string(out_proj, aux)
        call write_iter_nl(out_proj)
 
        do ik = 1, st%nik
          call write_iter_string(out_proj, "# w(ik)*occ(ist,ik)  ")
          do ist = gs_st%st_start, st%nst
            call write_iter_double(out_proj, st%kweights(ik)*st%occ(ist, ik), 1)
          end do
          call write_iter_nl(out_proj)
        end do
 
        call write_iter_header_start(out_proj)
        do ik = 1, st%nik
          do ist = gs_st%st_start, st%nst
            do uist = gs_st%st_start, gs_st%st_end
              write(aux, '(i4,a,i4)') ist, ' -> ', uist
              call write_iter_header(out_proj, 'Re {'//trim(aux)//'}')
              call write_iter_header(out_proj, 'Im {'//trim(aux)//'}')
            end do
          end do
        end do
        call write_iter_nl(out_proj)
 
      end if
 
      !The dipole matrix elements cannot be computed like that for solids
      if (.not. space%is_periodic()) then
 
        safe_allocate(projections(1:st%nst, gs_st%st_start:gs_st%st_end, 1:st%nik))
        do idir = 1, space%dim
          projections = m_z0
 
          call dipole_matrix_elements(idir)
 
          if (mpi_grp_is_root(mpi_world)) then
            write(aux, '(a,i1,a)') "<i|x_", idir, "|a>"
            call write_iter_string(out_proj, "# ------")
            call write_iter_header(out_proj, aux)
            do ik = 1, st%nik
              do ist = gs_st%st_start, st%st_end
                do uist = gs_st%st_start, gs_st%st_end
                  call write_iter_double(out_proj, real(projections(ist, uist, ik), real64), 1)
                  call write_iter_double(out_proj, aimag(projections(ist, uist, ik)), 1)
                end do
              end do
            end do
            call write_iter_nl(out_proj)
 
          end if
        end do
        safe_deallocate_a(projections)
 
      end if
 
      if (mpi_grp_is_root(mpi_world)) then
        call td_write_print_header_end(out_proj)
      end if
 
    end if
 
    safe_allocate(projections(1:st%nst, gs_st%st_start:gs_st%st_end, 1:st%nik))
    projections(:,:,:) = m_z0
    call calc_projections(mesh, st, gs_st, projections)
 
    if (mpi_grp_is_root(mpi_world)) then
      call write_iter_start(out_proj)
      do ik = 1, st%nik
        do ist = gs_st%st_start, st%nst
          do uist = gs_st%st_start, gs_st%st_end
            call write_iter_double(out_proj, real(projections(ist, uist, ik), real64), 1)
            call write_iter_double(out_proj, aimag(projections(ist, uist, ik)), 1)
          end do
        end do
      end do
      call write_iter_nl(out_proj)
    end if
 
    safe_deallocate_a(projections)
    pop_sub(td_write_proj)
 
  contains
    ! ---------------------------------------------------------
    subroutine dipole_matrix_elements(dir)
      integer, intent(in) :: dir
 
      integer :: uist, ist, ik, idim
      real(real64)   :: n_dip(space%dim)
      complex(real64), allocatable :: xpsi(:,:)
      complex(real64), allocatable :: psi(:, :), gspsi(:, :)
 
      push_sub(td_write_proj.dipole_matrix_elements)
 
      safe_allocate(psi(1:mesh%np, 1:st%d%dim))
      safe_allocate(gspsi(1:mesh%np, 1:st%d%dim))
      safe_allocate(xpsi(1:mesh%np, 1:st%d%dim))
 
      do ik = st%d%kpt%start, st%d%kpt%end
        do ist = st%st_start, st%st_end
          call states_elec_get_state(st, mesh, ist, ik, psi)
          do uist = gs_st%st_start, gs_st%st_end
            call states_elec_get_state(gs_st, mesh, uist, ik, gspsi)
 
            do idim = 1, st%d%dim
              xpsi(1:mesh%np, idim) = mesh%x(1:mesh%np, dir)*gspsi(1:mesh%np, idim)
            end do
            projections(ist, uist, ik) = -zmf_dotp(mesh, st%d%dim, psi, xpsi, reduce = .false.)
 
          end do
        end do
      end do
 
      safe_deallocate_a(xpsi)
      safe_deallocate_a(gspsi)
      safe_deallocate_a(psi)
 
      call comm_allreduce(st%dom_st_kpt_mpi_grp, projections)
 
      ! n_dip is not defined for more than space%dim
      n_dip = ions%dipole()
      do ik = 1, st%nik
        do ist = gs_st%st_start, st%nst
          do uist = gs_st%st_start, gs_st%st_end
            projections(ist, uist, ik) = projections(ist, uist, ik) - n_dip(dir)
          end do
        end do
      end do
 
 
      pop_sub(td_write_proj.dipole_matrix_elements)
    end subroutine dipole_matrix_elements
 
  end subroutine td_write_proj
 
  ! ---------------------------------------------------------
  ! ---------------------------------------------------------
  subroutine td_write_n_ex(out_nex, outp, namespace, mesh, kpoints, st, gs_st, iter)
    type(c_ptr),         intent(inout) :: out_nex
    type(output_t),      intent(in)    :: outp
    type(namespace_t),   intent(in)    :: namespace
    class(mesh_t),       intent(in)    :: mesh
    type(kpoints_t),     intent(in)    :: kpoints
    type(states_elec_t), intent(inout) :: st
    type(states_elec_t), intent(inout) :: gs_st
    integer,             intent(in)    :: iter
 
    complex(real64), allocatable :: projections(:,:)
    character(len=80) :: aux, dir
    integer :: ik, ikpt, ist, uist, err
    real(real64) :: Nex, weight
    integer :: gs_nst
    real(real64), allocatable :: Nex_kpt(:)
 
 
    push_sub(td_write_n_ex)
 
    if (iter == 0) then
      if (mpi_grp_is_root(mpi_world)) then
        call td_write_print_header_init(out_nex)
 
        write(aux, '(a15,i2)')      '# nspin        ', st%d%nspin
        call write_iter_string(out_nex, aux)
        call write_iter_nl(out_nex)
 
        call write_iter_string(out_nex, "#%")
        call write_iter_nl(out_nex)
 
        write(aux, '(a,i8)') "# nik  ", st%nik
        call write_iter_string(out_nex, aux)
        call write_iter_nl(out_nex)
 
        write(aux, '(a,2i8)') "#  st  ", gs_st%st_start, st%nst
        call write_iter_string(out_nex, aux)
        call write_iter_nl(out_nex)
 
        write(aux, '(a,2i8)') "# ust  ", gs_st%st_start, gs_st%st_end
        call write_iter_string(out_nex, aux)
        call write_iter_nl(out_nex)
 
        call write_iter_header_start(out_nex)
        call write_iter_header(out_nex, '#  iter t Nex(t)')
        call write_iter_nl(out_nex)
 
      end if
 
      if (mpi_grp_is_root(mpi_world)) then
        call td_write_print_header_end(out_nex)
      end if
 
    end if
 
    !We only need the occupied GS states
    if (.not. st%parallel_in_states) then
      do ik = st%d%kpt%start, st%d%kpt%end
        do ist = 1, gs_st%nst
          if (gs_st%occ(ist, ik) > m_min_occ) gs_nst = ist
        end do
      end do
    else
      gs_nst = gs_st%nst
    end if
 
    safe_allocate(projections(1:gs_nst, 1:st%nst))
 
    safe_allocate(nex_kpt(1:st%nik))
    nex_kpt = m_zero
    do ik = st%d%kpt%start, st%d%kpt%end
      ikpt = st%d%get_kpoint_index(ik)
      call zstates_elec_calc_projections(st, gs_st, namespace, mesh, ik, projections, gs_nst)
      do ist = 1, gs_nst
        weight = st%kweights(ik) * gs_st%occ(ist, ik)/ st%smear%el_per_state
        do uist = st%st_start, st%st_end
          nex_kpt(ikpt) = nex_kpt(ikpt) - weight * st%occ(uist, ik) * abs(projections(ist, uist))**2
        end do
      end do
      if (st%d%ispin == spin_polarized) then
        nex_kpt(ikpt) = nex_kpt(ikpt) + sum(st%occ(st%st_start:st%st_end, ik))*m_half*st%kweights(ik)
      else
        nex_kpt(ikpt) = nex_kpt(ikpt) + sum(st%occ(st%st_start:st%st_end, ik))*st%kweights(ik)
      end if
    end do
 
    if (st%parallel_in_states .or. st%d%kpt%parallel) then
      call comm_allreduce(st%st_kpt_mpi_grp, nex_kpt)
    end if
 
    nex = sum(nex_kpt)
 
    if (mpi_grp_is_root(mpi_world)) then
      call write_iter_start(out_nex)
      call write_iter_double(out_nex, nex, 1)
      call write_iter_nl(out_nex)
 
      ! now write down the k-resolved part
      write(dir, '(a,a,i7.7)') trim(outp%iter_dir),"td.", iter  ! name of directory
 
      call io_function_output_global_bz(outp%how(option__output__current_kpt) &
        + outp%how(option__output__density_kpt), dir, "n_excited_el_kpt", namespace, &
        kpoints, nex_kpt, unit_one, err)
    end if
 
    safe_deallocate_a(projections)
    safe_deallocate_a(nex_kpt)
 
    pop_sub(td_write_n_ex)
  end subroutine td_write_n_ex
 
  ! ---------------------------------------------------------
  ! ---------------------------------------------------------
  subroutine calc_projections(mesh, st, gs_st, projections)
    class(mesh_t),       intent(in)    :: mesh
    type(states_elec_t), intent(inout) :: st
    type(states_elec_t), intent(in)    :: gs_st
    complex(real64),     intent(inout) :: projections(1:st%nst, gs_st%st_start:gs_st%nst, 1:st%nik)
 
    integer :: uist, ist, ik
    complex(real64), allocatable :: psi(:, :), gspsi(:, :)
    push_sub(calc_projections)
 
    safe_allocate(psi(1:mesh%np, 1:st%d%dim))
    safe_allocate(gspsi(1:mesh%np, 1:st%d%dim))
 
    projections(:,:,:) = m_zero
 
    do ik = st%d%kpt%start, st%d%kpt%end
      do ist = st%st_start, st%st_end
        call states_elec_get_state(st, mesh, ist, ik, psi)
        do uist = gs_st%st_start, gs_st%nst
          call states_elec_get_state(gs_st, mesh, uist, ik, gspsi)
          projections(ist, uist, ik) = zmf_dotp(mesh, st%d%dim, psi, gspsi, reduce = .false.)
        end do
      end do
    end do
 
    safe_deallocate_a(psi)
    safe_deallocate_a(gspsi)
 
    call comm_allreduce(st%dom_st_kpt_mpi_grp, projections)
 
    pop_sub(calc_projections)
  end subroutine calc_projections
 
 
  subroutine td_write_proj_kp(mesh, kpoints, st, gs_st, namespace, iter)
    class(mesh_t),       intent(in)    :: mesh
    type(kpoints_t),     intent(in)    :: kpoints
    type(states_elec_t), intent(in)    :: st
    type(states_elec_t), intent(inout) :: gs_st
    type(namespace_t),   intent(in)    :: namespace
    integer,             intent(in)    :: iter
 
    complex(real64), allocatable :: proj(:,:), psi(:,:,:), gs_psi(:,:,:), temp_state(:,:)
    character(len=80) :: filename1, filename2
    integer :: ik,ist, jst, file, idim, nk_proj
 
    push_sub(td_write_proj_kp)
 
    write(filename1,'(I10)') iter
    filename1 = 'td.general/projections_iter_'//trim(adjustl(filename1))
    file = 9845623
 
    safe_allocate(proj(1:gs_st%nst, 1:gs_st%nst))
    safe_allocate(psi(1:gs_st%nst,1:gs_st%d%dim,1:mesh%np))
    safe_allocate(gs_psi(1:gs_st%nst,1:gs_st%d%dim,1:mesh%np))
    safe_allocate(temp_state(1:mesh%np,1:gs_st%d%dim))
 
    ! Project only k-points that have a zero weight.
    ! Why? It is unlikely that one is interested in the projections
    ! of the Monkhorst-Pack kpoints, but instead we assume that
    ! the user has specified a k-path with zero weights
    nk_proj = kpoints%nik_skip
 
    do ik = kpoints%reduced%npoints-nk_proj+1, kpoints%reduced%npoints
      ! reset arrays
      psi(1:gs_st%nst, 1:gs_st%d%dim, 1:mesh%np)= m_zero
      gs_psi(1:gs_st%nst, 1:gs_st%d%dim, 1:mesh%np)= m_zero
      ! open file for writing
      if (mpi_world%rank == 0) then
        write(filename2,'(I10)') ik
        filename2 = trim(adjustl(filename1))//'_ik_'//trim(adjustl(filename2))
        file = io_open(filename2, namespace, action='write')
      end if
      ! get all states at ik that are locally stored (ground state and td-states)
      do ist=gs_st%st_start,gs_st%st_end
        if (state_kpt_is_local(gs_st, ist, ik)) then
          call states_elec_get_state(st, mesh, ist, ik,temp_state)
          do idim = 1,gs_st%d%dim
            psi(ist,idim,1:mesh%np) =  temp_state(1:mesh%np,idim)
          end do
          call states_elec_get_state(gs_st, mesh, ist, ik, temp_state)
          do idim = 1,gs_st%d%dim
            gs_psi(ist,idim,1:mesh%np) =  temp_state(1:mesh%np,idim)
          end do
        end if
      end do
      ! collect states at ik from all processes in one array
      call comm_allreduce(mpi_world, psi)
      call comm_allreduce(mpi_world, gs_psi)
 
      ! compute the overlaps as a matrix product
      assert(mesh%np_global*gs_st%d%dim < huge(0_int32))
      proj(1:gs_st%nst, 1:gs_st%nst) = m_zero
      call zgemm('n',                         &
        'c',                                  &
        gs_st%nst,                            &
        gs_st%nst,                            &
        i8_to_i4(mesh%np_global*gs_st%d%dim), &
        cmplx(mesh%volume_element, m_zero, real64) , &
        psi(1, 1, 1),                         &
        ubound(psi, dim = 1),                 &
        gs_psi(1, 1, 1),                      &
        ubound(gs_psi, dim = 1),              &
        m_z0,                                 &
        proj(1, 1),                           &
        ubound(proj, dim = 1))
 
      ! write to file
      if (mpi_world%rank == 0) then
        do ist = 1, gs_st%nst
          do jst = 1, gs_st%nst
            write(file,'(I3,1x,I3,1x,e13.6,1x,e13.6,2x)') ist, jst, proj(ist,jst)
          end do
        end do
        call io_close(file)
      end if
 
    end do! ik
 
    safe_deallocate_a(proj)
    safe_deallocate_a(psi)
    safe_deallocate_a(gs_psi)
    safe_deallocate_a(temp_state)
 
    pop_sub(td_write_proj_kp)
  end subroutine td_write_proj_kp
 
  !---------------------------------------
  subroutine td_write_floquet(namespace, space, hm, ext_partners, mesh, st, iter)
    type(namespace_t),        intent(in)    :: namespace
    class(space_t),           intent(in)    :: space
    type(hamiltonian_elec_t), intent(inout) :: hm
    type(partner_list_t),     intent(in)    :: ext_partners
    class(mesh_t),            intent(in)    :: mesh
    type(states_elec_t),      intent(inout) :: st
    integer,                  intent(in)    :: iter
 
    complex(real64), allocatable :: hmss(:,:), psi(:,:,:), hpsi(:,:,:), temp_state1(:,:)
    complex(real64), allocatable :: HFloquet(:,:,:), HFloq_eff(:,:), temp(:,:)
    real(real64), allocatable :: eigenval(:), bands(:,:)
    character(len=80) :: filename
    integer :: it, nT, ik, ist, in, im, file, idim, nik, ik_count
    integer :: Forder, Fdim, m0, n0, n1, nst, ii, jj, lim_nst
    logical :: downfolding
    type(states_elec_t) :: hm_st
 
    real(real64) :: dt, Tcycle, omega
 
    push_sub(td_write_floquet)
 
    downfolding = .false.
 
    ! this does not depend on propagation, so we do it only once
    if (.not. iter == 0) then
      pop_sub(td_write_floquet)
      return
    end if
 
    nst = st%nst
 
    !for now no domain distributionallowed
    assert(mesh%np == mesh%np_global)
 
    ! this is used to initialize the hpsi (more effiecient ways?)
    call states_elec_copy(hm_st, st)
 
    !%Variable TDFloquetFrequency
    !%Type float
    !%Default 0
    !%Section Time-Dependent::TD Output
    !%Description
    !% Frequency for the Floquet analysis, this should be the carrier frequency or integer multiples of it.
    !% Other options will work, but likely be nonsense.
    !%
    !%End
    call parse_variable(namespace, 'TDFloquetFrequency', m_zero, omega, units_inp%energy)
    call messages_print_var_value('Frequency used for Floquet analysis', omega, namespace=namespace)
    if (abs(omega) <= m_epsilon) then
      message(1) = "Please give a non-zero value for TDFloquetFrequency"
      call messages_fatal(1, namespace=namespace)
    end if
 
    ! get time of one cycle
    tcycle = m_two * m_pi / omega
 
    !%Variable TDFloquetSample
    !%Type integer
    !%Default 20
    !%Section Time-Dependent::TD Output
    !%Description
    !% Number of points on which one Floquet cycle is sampled in the time-integral of the Floquet analysis.
    !%
    !%End
    call parse_variable(namespace, 'TDFloquetSample',20 ,nt)
    call messages_print_var_value('Number of Floquet time-sampling points', nt, namespace=namespace)
    dt = tcycle/real(nt, real64)
 
    !%Variable TDFloquetDimension
    !%Type integer
    !%Default -1
    !%Section Time-Dependent::TD Output
    !%Description
    !% Order of Floquet Hamiltonian. If negative number is given, downfolding is performed.
    !%End
    call parse_variable(namespace, 'TDFloquetDimension',-1,forder)
    if (forder .ge. 0) then
      call messages_print_var_value('Order of multiphoton Floquet-Hamiltonian', forder, namespace=namespace)
      !Dimension of multiphoton Floquet-Hamiltonian
      fdim = 2 * forder + 1
    else
      message(1) = 'Floquet-Hamiltonian is downfolded'
      call messages_info(1, namespace=namespace)
      downfolding = .true.
      forder = 1
      fdim = 3
    end if
 
    dt = tcycle/real(nt, real64)
 
    ! we are only interested for k-point with zero weight
    nik = hm%kpoints%nik_skip
 
    safe_allocate(hmss(1:nst,1:nst))
    safe_allocate( psi(1:nst,1:st%d%dim,1:mesh%np))
    safe_allocate(hpsi(1:nst,1:st%d%dim,1:mesh%np))
    safe_allocate(temp_state1(1:mesh%np,1:st%d%dim))
 
    ! multiphoton Floquet Hamiltonian, layout:
    !     (H_{-n,-m} ...  H_{-n,0} ...  H_{-n,m})
    !     (    .      .      .      .      .    )
    ! H = (H_{0,-m}  ...  H_{0,0}  ...  H_{0,m} )
    !     (    .      .      .      .      .    )
    !     (H_{n,-m}  ...  H_{n,0}  ...  H_{n,m} )
    safe_allocate(hfloquet(1:nik,1:nst*fdim, 1:nst*fdim))
    hfloquet(1:nik,1:nst*fdim, 1:nst*fdim) = m_zero
 
    ! perform time-integral over one cycle
    do it = 1, nt
      ! get non-interacting Hamiltonian at time (offset by one cycle to allow for ramp)
      call hm%update(mesh, namespace, space, ext_partners, time=tcycle+it*dt)
      ! get hpsi
      call zhamiltonian_elec_apply_all(hm, namespace, mesh, st, hm_st)
 
      ! project Hamiltonian into grounstates for zero weight k-points
      ik_count = 0
 
      do ik = hm%kpoints%reduced%npoints-nik+1, hm%kpoints%reduced%npoints
        ik_count = ik_count + 1
 
        psi(1:nst, 1:st%d%dim, 1:mesh%np)= m_zero
        hpsi(1:nst, 1:st%d%dim, 1:mesh%np)= m_zero
 
        do ist = st%st_start, st%st_end
          if (state_kpt_is_local(st, ist, ik)) then
            call states_elec_get_state(st, mesh, ist, ik,temp_state1)
            do idim = 1, st%d%dim
              psi(ist, idim, 1:mesh%np) = temp_state1(1:mesh%np,idim)
            end do
            call states_elec_get_state(hm_st, mesh, ist, ik,temp_state1)
            do idim = 1, st%d%dim
              hpsi(ist, idim,1:mesh%np) = temp_state1(1:mesh%np,idim)
            end do
          end if
        end do
        call comm_allreduce(mpi_world, psi)
        call comm_allreduce(mpi_world, hpsi)
        assert(mesh%np_global*st%d%dim < huge(0_int32))
        hmss(1:nst,1:nst) = m_zero
        call zgemm( 'n',                        &
          'c',                                  &
          nst,                                  &
          nst,                                  &
          i8_to_i4(mesh%np_global*st%d%dim),    &
          cmplx(mesh%volume_element, m_zero, real64) , &
          hpsi(1, 1, 1),                        &
          ubound(hpsi, dim = 1),                &
          psi(1, 1, 1),                         &
          ubound(psi, dim = 1),                 &
          m_z0,                                 &
          hmss(1, 1),                           &
          ubound(hmss, dim = 1))
 
        hmss(1:nst,1:nst) = conjg(hmss(1:nst,1:nst))
 
        ! accumulate the Floqeut integrals
        do in = -forder, forder
          do im = -forder, forder
            ii = (in+forder) * nst
            jj = (im+forder) * nst
            hfloquet(ik_count, ii+1:ii+nst, jj+1:jj+nst) =  &
              hfloquet(ik_count, ii+1:ii+nst, jj+1:jj+nst) + hmss(1:nst, 1:nst) * exp(-(in-im)*m_zi*omega*it*dt)
            ! diagonal term
            if (in == im) then
              do ist = 1, nst
                hfloquet(ik_count, ii+ist, ii+ist) = hfloquet(ik_count, ii+ist, ii+ist) + in*omega
              end do
            end if
          end do
        end do
      end do !ik
 
    end do ! it
 
    hfloquet(:,:,:) = m_one/nt*hfloquet(:,:,:)
 
    ! diagonalize Floquet Hamiltonian
    if (downfolding) then
      ! here perform downfolding
      safe_allocate(hfloq_eff(1:nst,1:nst))
      safe_allocate(eigenval(1:nst))
      safe_allocate(bands(1:nik,1:nst))
 
      hfloq_eff(1:nst,1:nst) = m_zero
      do ik = 1, nik
        ! the HFloquet blocks are copied directly out of the super matrix
        m0 = nst ! the m=0 start position
        n0 = nst ! the n=0 start postion
        n1 = 2*nst ! the n=+1 start postion
        hfloq_eff(1:nst, 1:nst) = hfloquet(ik, n0+1:n0+nst, m0+1:m0+nst) + &
          m_one/omega*(matmul(hfloquet(ik, 1:nst, m0+1:m0+nst), hfloquet(ik, n1+1:n1+nst, m0+1:m0+nst))- &
          matmul(hfloquet(ik, n1+1:n1+nst, m0+1:m0+nst), hfloquet(ik, 1:nst, m0+1:m0+nst)))
 
        call lalg_eigensolve(nst, hfloq_eff, eigenval)
        bands(ik,1:nst) = eigenval(1:nst)
      end do
      safe_deallocate_a(hfloq_eff)
    else
      ! the full Floquet
      safe_allocate(eigenval(1:nst*fdim))
      safe_allocate(bands(1:nik,1:nst*fdim))
      safe_allocate(temp(1:nst*fdim, 1:nst*fdim))
 
      do ik = 1, nik
        temp(1:nst*fdim, 1:nst*fdim) = hfloquet(ik, 1:nst*fdim, 1:nst*fdim)
        call lalg_eigensolve(nst*fdim, temp, eigenval)
        bands(ik, 1:nst*fdim) = eigenval(1:nst*fdim)
      end do
    end if
 
    !write bandstructure to file
    if (downfolding) then
      lim_nst = nst
      filename = "downfolded_floquet_bands"
    else
      lim_nst = nst*fdim
      filename = "floquet_bands"
    end if
    ! write bands (full or downfolded)
    if (mpi_world%rank == 0) then
      file = 987254
      file = io_open(filename, namespace, action = 'write')
      do ik = 1,nik
        do ist = 1,lim_nst
          write(file,'(e13.6, 1x)',advance='no') bands(ik,ist)
        end do
        write(file,'(1x)')
      end do
      call io_close(file)
    end if
 
    if (.not. downfolding) then
      ! for the full Floquet case compute also the trivially shifted
      ! Floquet bands for reference (i.e. setting H_{nm}=0 for n!=m)
      bands(1:nik, 1:nst*fdim) = m_zero
      do ik = 1, nik
        temp(1:nst*fdim,1:nst*fdim) = m_zero
        do jj = 0, fdim - 1
          ii = jj * nst
          temp(ii+1:ii+nst, ii+1:ii+nst) = hfloquet(ik, ii+1:ii+nst, ii+1:ii+nst)
        end do
        call lalg_eigensolve(nst*fdim, temp, eigenval)
        bands(ik, 1:nst*fdim) = eigenval(1:nst*fdim)
      end do
 
      if (mpi_world%rank == 0) then
        filename = 'trivial_floquet_bands'
        file = io_open(filename, namespace, action = 'write')
        do ik=1,nik
          do ist = 1,lim_nst
            write(file,'(e13.6, 1x)', advance='no') bands(ik,ist)
          end do
          write(file,'(1x)')
        end do
        call io_close(file)
      end if
    end if
 
    ! reset time in Hamiltonian
    call hm%update(mesh, namespace, space, ext_partners, time=m_zero)
 
    safe_deallocate_a(hmss)
    safe_deallocate_a(psi)
    safe_deallocate_a(hpsi)
    safe_deallocate_a(temp_state1)
    safe_deallocate_a(hfloquet)
    safe_deallocate_a(eigenval)
    safe_deallocate_a(bands)
    safe_deallocate_a(temp)
 
    pop_sub(td_write_floquet)
 
  end subroutine td_write_floquet
 
  ! ---------------------------------------------------------
  subroutine td_write_total_current(out_total_current, space, mesh, st, iter)
    type(c_ptr),         intent(inout) :: out_total_current
    class(space_t),      intent(in)    :: space
    class(mesh_t),       intent(in)    :: mesh
    type(states_elec_t), intent(in)    :: st
    integer,             intent(in)    :: iter
 
    integer :: idir, ispin
    character(len=50) :: aux
    real(real64) :: total_current(space%dim), abs_current(space%dim)
 
    push_sub(td_write_total_current)
 
    if (mpi_grp_is_root(mpi_world) .and. iter == 0) then
      call td_write_print_header_init(out_total_current)
 
      ! first line: column names
      call write_iter_header_start(out_total_current)
 
      do idir = 1, space%dim
        write(aux, '(a2,i1,a1)') 'I(', idir, ')'
        call write_iter_header(out_total_current, aux)
      end do
 
      do idir = 1, space%dim
        write(aux, '(a10,i1,a1)') 'IntAbs(j)(', idir, ')'
        call write_iter_header(out_total_current, aux)
      end do
 
      do ispin = 1, st%d%nspin
        do idir = 1, space%dim
          write(aux, '(a4,i1,a1,i1,a1)') 'I-sp', ispin, '(', idir, ')'
          call write_iter_header(out_total_current, aux)
        end do
      end do
 
      call write_iter_nl(out_total_current)
 
      call td_write_print_header_end(out_total_current)
    end if
 
    assert(allocated(st%current))
 
    if (mpi_grp_is_root(mpi_world)) then
      call write_iter_start(out_total_current)
    end if
 
    total_current = 0.0_real64
    do idir = 1, space%dim
      do ispin = 1, st%d%spin_channels
        total_current(idir) =  total_current(idir) + dmf_integrate(mesh, st%current(:, idir, ispin), reduce = .false.)
      end do
      total_current(idir) = units_from_atomic(units_out%length/units_out%time, total_current(idir))
    end do
    if (mesh%parallel_in_domains) then
      call mesh%allreduce(total_current, dim = space%dim)
    end if
 
    abs_current = 0.0_real64
    do idir = 1, space%dim
      do ispin = 1, st%d%spin_channels
        abs_current(idir) =  abs_current(idir) + dmf_integrate(mesh, abs(st%current(:, idir, ispin)), reduce = .false.)
      end do
      abs_current(idir) = units_from_atomic(units_out%length/units_out%time, abs_current(idir))
    end do
    if (mesh%parallel_in_domains) then
      call mesh%allreduce(abs_current, dim = space%dim)
    end if
 
    if (mpi_grp_is_root(mpi_world)) then
      call write_iter_double(out_total_current, total_current, space%dim)
      call write_iter_double(out_total_current, abs_current, space%dim)
    end if
 
    do ispin = 1, st%d%nspin
      total_current = 0.0_real64
      do idir = 1, space%dim
        total_current(idir) = units_from_atomic(units_out%length/units_out%time, &
          dmf_integrate(mesh, st%current(:, idir, ispin), reduce = .false.))
      end do
      if (mesh%parallel_in_domains) then
        call mesh%allreduce(total_current, dim = space%dim)
      end if
      if (mpi_grp_is_root(mpi_world)) then
        call write_iter_double(out_total_current, total_current, space%dim)
      end if
    end do
 
    if (mpi_grp_is_root(mpi_world)) then
      call write_iter_nl(out_total_current)
    end if
 
    pop_sub(td_write_total_current)
  end subroutine td_write_total_current
 
  ! ---------------------------------------------------------
 
  subroutine td_write_total_heat_current(write_obj, space, hm, gr, st, iter)
    type(c_ptr),              intent(inout) :: write_obj
    class(space_t),           intent(in)    :: space
    type(hamiltonian_elec_t), intent(inout) :: hm
    type(grid_t),             intent(in)    :: gr
    type(states_elec_t),      intent(in)    :: st
    integer,                  intent(in)    :: iter
 
    integer :: idir, ispin
    character(len=50) :: aux
    real(real64), allocatable :: heat_current(:, :, :)
    real(real64) :: total_current(space%dim)
 
    push_sub(td_write_total_current)
 
    if (mpi_grp_is_root(mpi_world) .and. iter == 0) then
      call td_write_print_header_init(write_obj)
 
      ! first line: column names
      call write_iter_header_start(write_obj)
 
      do idir = 1, space%dim
        write(aux, '(a2,i1,a1)') 'Jh(', idir, ')'
        call write_iter_header(write_obj, aux)
      end do
 
      call write_iter_nl(write_obj)
 
      call td_write_print_header_end(write_obj)
    end if
 
    safe_allocate(heat_current(1:gr%np, 1:space%dim, 1:st%d%nspin))
 
    call current_heat_calculate(space, gr%der, hm, st, heat_current)
 
    if (mpi_grp_is_root(mpi_world)) call write_iter_start(write_obj)
 
    total_current = 0.0_real64
    do idir = 1, space%dim
      do ispin = 1, st%d%spin_channels
        total_current(idir) =  total_current(idir) + dmf_integrate(gr, heat_current(:, idir, ispin))
      end do
      total_current(idir) = units_from_atomic(units_out%energy*units_out%length/units_out%time, total_current(idir))
    end do
 
    safe_deallocate_a(heat_current)
 
    if (mpi_grp_is_root(mpi_world)) call write_iter_double(write_obj, total_current, space%dim)
 
    if (mpi_grp_is_root(mpi_world)) call write_iter_nl(write_obj)
 
    pop_sub(td_write_total_current)
  end subroutine td_write_total_heat_current
 
 
  ! ---------------------------------------------------------
  subroutine td_write_partial_charges(out_partial_charges, mesh, st, ions, iter)
    type(c_ptr),             intent(inout) :: out_partial_charges
    class(mesh_t),           intent(in)    :: mesh
    type(states_elec_t),     intent(in)    :: st
    type(ions_t),            intent(in)    :: ions
    integer,                 intent(in)    :: iter
 
    integer :: idir
    character(len=50) :: aux
    real(real64), allocatable :: hirshfeld_charges(:)
 
    push_sub(td_write_partial_charges)
 
    safe_allocate(hirshfeld_charges(1:ions%natoms))
 
    call partial_charges_calculate(mesh, st, ions, hirshfeld_charges)
 
    if (mpi_grp_is_root(mpi_world)) then
 
      if (iter == 0) then
 
        call td_write_print_header_init(out_partial_charges)
 
        ! first line: column names
        call write_iter_header_start(out_partial_charges)
 
        do idir = 1, ions%natoms
          write(aux, '(a13,i3,a1)') 'hirshfeld(atom=', idir, ')'
          call write_iter_header(out_partial_charges, aux)
        end do
 
        call write_iter_nl(out_partial_charges)
 
        call td_write_print_header_end(out_partial_charges)
      end if
 
      call write_iter_start(out_partial_charges)
 
      call write_iter_double(out_partial_charges, hirshfeld_charges, ions%natoms)
 
      call write_iter_nl(out_partial_charges)
    end if
 
    safe_deallocate_a(hirshfeld_charges)
 
    pop_sub(td_write_partial_charges)
  end subroutine td_write_partial_charges
 
  ! ---------------------------------------------------------
  subroutine td_write_q(out_q, space, ks, iter)
    type(c_ptr),             intent(inout) :: out_q
    class(space_t),          intent(in)    :: space
    type(v_ks_t),            intent(in)    :: ks
    integer,                 intent(in)    :: iter
 
    integer :: ii
    character(len=50) :: aux
 
    push_sub(td_write_q)
 
    if (mpi_grp_is_root(mpi_world)) then
      if (iter == 0) then
        call td_write_print_header_init(out_q)
        call write_iter_header_start(out_q)
        do ii = 1, ks%pt%nmodes
          write(aux, '(a1,i3,a3)') 'q', ii, '(t)'
          call write_iter_header(out_q, aux)
        end do
        do ii = 1, ks%pt%nmodes
          write(aux, '(a1,i3,a3)') 'p', ii, '(t)'
          call write_iter_header(out_q, aux)
        end do
        do ii = 1, space%dim
          write(aux, '(a3,i3,a3)') 'f_pt', ii, '(t)'
          call write_iter_header(out_q, aux)
        end do
        call write_iter_nl(out_q)
        call td_write_print_header_end(out_q)
      end if
 
      call write_iter_start(out_q)
      call write_iter_double(out_q, ks%pt_mx%pt_q, ks%pt%nmodes)
      call write_iter_double(out_q, ks%pt_mx%pt_p, ks%pt%nmodes)
      call write_iter_double(out_q, ks%pt_mx%fmf, space%dim)
      call write_iter_nl(out_q)
    end if
 
    pop_sub(td_write_q)
  end subroutine td_write_q
 
 
  ! ---------------------------------------------------------
  subroutine td_write_mxll_field(out_mxll, space, hm, dt, iter)
    type(c_ptr),             intent(inout) :: out_mxll
    class(space_t),          intent(in)    :: space
    type(hamiltonian_elec_t),intent(in)    :: hm
    real(real64),            intent(in)    :: dt
    integer,                 intent(in)    :: iter
 
    integer :: idir
    real(real64) :: field(space%dim)
    character(len=80) :: aux
    character(len=1) :: field_char
 
    push_sub(td_write_mxll_field)
 
    if (.not. mpi_grp_is_root(mpi_world)) then
      pop_sub(td_write_mxll_field)
      return ! only first node outputs
    end if
 
    if (iter == 0) then
      call td_write_print_header_init(out_mxll)
 
      write(aux, '(a7,e20.12,3a)') '# dt = ', units_from_atomic(units_out%time, dt), &
        " [", trim(units_abbrev(units_out%time)), "]"
      call write_iter_string(out_mxll, aux)
      call write_iter_nl(out_mxll)
 
      call write_iter_header_start(out_mxll)
      select case (hm%mxll%coupling_mode)
      case (length_gauge_dipole, multipolar_expansion)
        if (hm%mxll%add_electric_dip) field_char = 'E'
        if (hm%mxll%add_magnetic_dip) field_char = 'B'
        do idir = 1, space%dim
          write(aux, '(a,i1,a)') field_char // '(', idir, ')'
          call write_iter_header(out_mxll, aux)
        end do
      case (velocity_gauge_dipole)
        do idir = 1, space%dim
          write(aux, '(a,i1,a)') 'A(', idir, ')'
          call write_iter_header(out_mxll, aux)
        end do
      end select
      call write_iter_nl(out_mxll)
 
      call write_iter_string(out_mxll, '#[Iter n.]')
      call write_iter_header(out_mxll, '[' // trim(units_abbrev(units_out%time)) // ']')
 
      ! Note that we do not print out units of E, A or B, but rather units of e*E, e*A or e*B
      ! (force, energy and magnetic flux density, respectively).
      select case (hm%mxll%coupling_mode)
      case (length_gauge_dipole, multipolar_expansion)
        if (hm%mxll%add_electric_dip) aux = '[' // trim(units_abbrev(units_out%force)) // ']'
        if (hm%mxll%add_magnetic_dip) aux = '[' // trim(units_abbrev(unit_one/units_out%length**2)) // ']'
        do idir = 1, space%dim
          call write_iter_header(out_mxll, aux)
        end do
      case (velocity_gauge_dipole)
        aux = '[' // trim(units_abbrev(units_out%energy)) // ']'
        do idir = 1, space%dim
          call write_iter_header(out_mxll, aux)
        end do
      end select
      call write_iter_nl(out_mxll)
      call td_write_print_header_end(out_mxll)
    end if
 
    call write_iter_start(out_mxll)
 
    field = m_zero
    select case (hm%mxll%coupling_mode)
    case (length_gauge_dipole, multipolar_expansion)
      if (hm%mxll%add_electric_dip) field = units_from_atomic(units_out%force, hm%mxll%e_field_dip)
      if (hm%mxll%add_magnetic_dip) field = units_from_atomic(unit_one/units_out%length**2, hm%mxll%b_field_dip)
      call write_iter_double(out_mxll, field, space%dim)
    case (velocity_gauge_dipole)
      field = units_from_atomic(units_out%energy, hm%mxll%vec_pot_dip)
      call write_iter_double(out_mxll, field, space%dim)
    end select
    call write_iter_nl(out_mxll)
 
    pop_sub(td_write_mxll_field)
  end subroutine td_write_mxll_field
 
 
  ! ---------------------------------------------------------
  subroutine td_write_effective_u(out_coords, lda_u, iter)
    type(c_ptr),       intent(inout) :: out_coords
    type(lda_u_t),     intent(in)    :: lda_u
    integer,           intent(in)    :: iter
 
    integer :: ios, inn
    character(len=50) :: aux
 
    if (.not. mpi_grp_is_root(mpi_world)) return ! only first node outputs
 
    push_sub(td_write_effective_u)
 
    if (iter == 0) then
      call td_write_print_header_init(out_coords)
 
      ! first line: column names
      call write_iter_header_start(out_coords)
 
      do ios = 1, lda_u%norbsets
        write(aux, '(a2,i3,a1)') 'Ueff(', ios, ')'
        call write_iter_header(out_coords, aux)
      end do
 
      do ios = 1, lda_u%norbsets
        write(aux, '(a2,i3,a1)') 'U(', ios, ')'
        call write_iter_header(out_coords, aux)
      end do
 
      do ios = 1, lda_u%norbsets
        write(aux, '(a2,i3,a1)') 'J(', ios, ')'
        call write_iter_header(out_coords, aux)
      end do
 
      if (lda_u%intersite) then
        do ios = 1, lda_u%norbsets
          do inn = 1, lda_u%orbsets(ios)%nneighbors
            write(aux, '(a2,i3,a1,i3,a1)') 'V(', ios,'-', inn, ')'
            call write_iter_header(out_coords, aux)
          end do
        end do
      end if
 
 
      call write_iter_nl(out_coords)
 
      ! second line: units
      call write_iter_string(out_coords, '#[Iter n.]')
      call write_iter_header(out_coords, '[' // trim(units_abbrev(units_out%time)) // ']')
      call write_iter_string(out_coords, &
        'Effective U '   // trim(units_abbrev(units_out%energy))   //   &
        ', U in '// trim(units_abbrev(units_out%energy)) //   &
        ', J in '    // trim(units_abbrev(units_out%energy)))
      call write_iter_nl(out_coords)
 
      call td_write_print_header_end(out_coords)
    end if
 
    call write_iter_start(out_coords)
 
    do ios = 1, lda_u%norbsets
      call write_iter_double(out_coords, units_from_atomic(units_out%energy,  &
        lda_u%orbsets(ios)%Ueff), 1)
    end do
 
    do ios = 1, lda_u%norbsets
      call write_iter_double(out_coords, units_from_atomic(units_out%energy,  &
        lda_u%orbsets(ios)%Ubar), 1)
    end do
 
    do ios = 1, lda_u%norbsets
      call write_iter_double(out_coords, units_from_atomic(units_out%energy,  &
        lda_u%orbsets(ios)%Jbar), 1)
    end do
 
    if (lda_u%intersite) then
      do ios = 1, lda_u%norbsets
        do inn = 1, lda_u%orbsets(ios)%nneighbors
          call write_iter_double(out_coords, units_from_atomic(units_out%energy,  &
            lda_u%orbsets(ios)%V_ij(inn,0)), 1)
        end do
      end do
    end if
 
    call write_iter_nl(out_coords)
 
    pop_sub(td_write_effective_u)
  end subroutine td_write_effective_u
 
  subroutine td_write_norm_ks_orbitals(file_handle, grid, kpoints, st, iter)
    type(c_ptr),              intent(inout) :: file_handle
    type(grid_t),             intent(in)    :: grid
    type(kpoints_t),          intent(in)    :: kpoints
    type(states_elec_t),      intent(in)    :: st
    integer,                  intent(in)    :: iter
 
    integer            :: ik_ispin, ist
    character(len=7)   :: nkpt_str, nst_str
    character(len=7)   :: ik_str, ist_str
    real(real64), allocatable :: norm_ks(:, :)
    real(real64)       :: n_electrons
 
    push_sub(td_write_norm_ks_orbitals)
 
    safe_allocate(norm_ks(1:st%nst, 1:st%nik))
    call states_elec_calc_norms(grid, kpoints, st, norm_ks)
 
    if (mpi_grp_is_root(mpi_world)) then
      ! Header
      if (iter == 0) then
        call td_write_print_header_init(file_handle)
 
        ! Line 1
        write(nkpt_str, '(I7)') st%nik
        write(nst_str, '(I7)') st%nst
        call write_iter_string(file_handle, '# Dimensions. (nstates, nkpt * nspin):')
        call write_iter_string(file_handle,  trim(adjustl(nst_str)) // ' ' // trim(adjustl(nkpt_str)))
        call write_iter_nl(file_handle)
 
        ! Line 2
        call write_iter_string(file_handle,  '# Norm ordering: (istate, ikpoint_spin)')
        call write_iter_nl(file_handle)
 
        ! Line 3
        call write_iter_header_start(file_handle)
        call write_iter_header(file_handle, 'N_electrons')
        do ik_ispin = 1, st%nik
          do ist = 1, st%nst
            write(ik_str, '(I7)') ik_ispin
            write(ist_str, '(I7)') ist
            call write_iter_header(file_handle, &
              'Norm (' // trim(ist_str) // ',' // trim(ik_str) // ')')
          enddo
        enddo
        call write_iter_nl(file_handle)
        call td_write_print_header_end(file_handle)
      end if
 
      n_electrons = sum(st%occ * norm_ks**2)
 
      ! Output norms for time step `iter`
      call write_iter_start(file_handle)
      call write_iter_double(file_handle, n_electrons, 1)
      do ik_ispin = 1, st%nik
        call write_iter_double(file_handle, norm_ks(:, ik_ispin), size(norm_ks, 1))
      enddo
      call write_iter_nl(file_handle)
 
    end if
 
    safe_deallocate_a(norm_ks)
 
    pop_sub(td_write_norm_ks_orbitals)
 
  end subroutine td_write_norm_ks_orbitals
 
  ! ---------------------------------------------------------
  subroutine td_write_mxll_init(writ, namespace, iter, dt)
    type(td_write_t),         intent(out)   :: writ
    type(namespace_t),        intent(in)    :: namespace
    integer,                  intent(in)    :: iter
    real(real64),             intent(in)    :: dt
 
    integer :: default, flags, iout, first
 
    push_sub(td_write_mxll_init)
 
    !%Variable MaxwellTDOutput
    !%Type flag
    !%Default maxwell_energy
    !%Section Maxwell::Output
    !%Description
    !% Defines what should be output during the time-dependent
    !% Maxwell simulation. Many of the options can increase the computational
    !% cost of the simulation, so only use the ones that you need. In
    !% most cases the default value is enough, as it is adapted to the
    !% details of the TD run.
    !% WARNING: the calculation of the longitudinal or transverse E and B fields
    !% can be very expensive, so please consider using the MaxwellOutput block
    !% to calculate and output these quantities at certain timesteps.
    !%Option maxwell_total_e_field 1
    !% Output of the total (longitudinal plus transverse) electric field at
    !% the points specified in the MaxwellFieldsCoordinate block
    !%Option maxwell_total_b_field 8
    !% Output of the total (longitudinal plus transverse) magnetic field at
    !% the points specified in the MaxwellFieldsCoordinate block
    !%Option maxwell_longitudinal_e_field 64
    !% Output of the longitudinal electric field at the points
    !% specified in the MaxwellFieldsCoordinate block  (can slow down the run)
    !%Option maxwell_longitudinal_b_field 512
    !% Output of the longitudinal magnetic field at the points
    !% specified in the MaxwellFieldsCoordinate block  (can slow down the run)
    !%Option maxwell_transverse_e_field 4096
    !% Output of the transverse electric field at the points
    !% specified in the MaxwellFieldsCoordinate block  (can slow down the run)
    !%Option maxwell_transverse_b_field 32768
    !% Output of the transverse magnetic field at the points
    !% specified in the MaxwellFieldsCoordinate block  (can slow down the run)
    !%Option maxwell_energy 262144
    !% Output of the electromagnetic field energy into the folder <tt>td.general/maxwell</tt>.
    !% WARNING: the transverse and longitudinal energies will be correct only if you request
    !% the longitudinal and transverse E or B fields as output. Otherwise they will be set to
    !% zero.
    !%Option e_field_surface_x 524288
    !% Output of the E field sliced along the plane x=0 for each field component
    !%Option e_field_surface_y 1048576
    !% Output of the E field sliced along the plane y=0 for each field component
    !%Option e_field_surface_z 2097152
    !% Output of the E field sliced along the plane z=0 for each field component
    !%Option b_field_surface_x 4194304
    !% Output of the B field sliced along the plane x=0 for each field component
    !%Option b_field_surface_y 8388608
    !% Output of the B field sliced along the plane y=0 for each field component
    !%Option b_field_surface_z 16777216
    !% Output of the B field sliced along the plane z=0 for each field component
    !%End
 
    default = 2**(out_maxwell_energy - 1)
    call parse_variable(namespace, 'MaxwellTDOutput', default, flags)
 
    if (.not. varinfo_valid_option('MaxwellTDOutput', flags, is_flag = .true.)) then
      call messages_input_error(namespace, 'MaxwellTDOutput')
    end if
 
    do iout = out_maxwell_total_e_field, out_maxwell_trans_b_field, 3
      writ%out(iout)%write     = (iand(flags, 2**(iout     - 1)) /= 0)
      if (writ%out(iout)%write) then
        writ%out(iout + 1)%write = .true.
        writ%out(iout + 2)%write = .true.
      end if
    end do
 
    do iout = out_maxwell_energy, out_maxwell_max
      writ%out(iout)%write = (iand(flags, 2**(iout - 1)) /= 0)
    end do
 
    if (iter == 0) then
      first = 0
    else
      first = iter + 1
    end if
 
    call io_mkdir('td.general', namespace)
 
    ! total E field
    if (writ%out(out_maxwell_total_e_field)%write) then
      call write_iter_init(writ%out(out_maxwell_total_e_field)%handle, first, &
        units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/total_e_field_x", namespace)))
      call write_iter_init(writ%out(out_maxwell_total_e_field + 1)%handle, first, &
        units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/total_e_field_y", namespace)))
      call write_iter_init(writ%out(out_maxwell_total_e_field + 2)%handle, first, &
        units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/total_e_field_z", namespace)))
    end if
 
    ! total B field
    if (writ%out(out_maxwell_total_b_field)%write) then
      call write_iter_init(writ%out(out_maxwell_total_b_field)%handle, first, &
        units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/total_b_field_x", namespace)))
      call write_iter_init(writ%out(out_maxwell_total_b_field + 1)%handle, first, &
        units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/total_b_field_y", namespace)))
      call write_iter_init(writ%out(out_maxwell_total_b_field + 2)%handle, first, &
        units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/total_b_field_z", namespace)))
    end if
 
    ! longitudinal E field
    if (writ%out(out_maxwell_long_e_field)%write) then
      call write_iter_init(writ%out(out_maxwell_long_e_field)%handle, first, &
        units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/longitudinal_e_field_x", namespace)))
      call write_iter_init(writ%out(out_maxwell_long_e_field + 1)%handle, first, &
        units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/longitudinal_e_field_y", namespace)))
      call write_iter_init(writ%out(out_maxwell_long_e_field + 2)%handle, first, &
        units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/longitudinal_e_field_z", namespace)))
    end if
 
    ! longitudinal B field
    if (writ%out(out_maxwell_long_b_field)%write) then
      call write_iter_init(writ%out(out_maxwell_long_b_field)%handle, first, &
        units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/longitudinal_b_field_x", namespace)))
      call write_iter_init(writ%out(out_maxwell_long_b_field + 1)%handle, first, &
        units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/longitudinal_b_field_y", namespace)))
      call write_iter_init(writ%out(out_maxwell_long_b_field + 2)%handle, first, &
        units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/longitudinal_b_field_z", namespace)))
    end if
 
    ! transverse E field
    if (writ%out(out_maxwell_trans_e_field)%write) then
      call write_iter_init(writ%out(out_maxwell_trans_e_field)%handle, first, &
        units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/transverse_e_field_x", namespace)))
      call write_iter_init(writ%out(out_maxwell_trans_e_field + 1)%handle, first, &
        units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/transverse_e_field_y", namespace)))
      call write_iter_init(writ%out(out_maxwell_trans_e_field + 2)%handle, first, &
        units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/transverse_e_field_z", namespace)))
    end if
 
    ! transverse B field
    if (writ%out(out_maxwell_trans_b_field)%write) then
      call write_iter_init(writ%out(out_maxwell_trans_b_field)%handle, first, &
        units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/transverse_b_field_x", namespace)))
      call write_iter_init(writ%out(out_maxwell_trans_b_field + 1)%handle, first, &
        units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/transverse_b_field_y", namespace)))
      call write_iter_init(writ%out(out_maxwell_trans_b_field + 2)%handle, first, &
        units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/transverse_b_field_z", namespace)))
    end if
 
    ! Maxwell energy
    if (writ%out(out_maxwell_energy)%write) then
      call write_iter_init(writ%out(out_maxwell_energy)%handle, first, &
        units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/maxwell_energy", namespace)))
    end if
 
    if (writ%out(out_e_field_surface_x)%write) then
      call write_iter_init(writ%out(out_e_field_surface_x)%handle, first, &
        units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/electric_field_surface-x", namespace)))
    end if
 
    if (writ%out(out_e_field_surface_y)%write) then
      call write_iter_init(writ%out(out_e_field_surface_y)%handle, first, &
        units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/electric_field_surface-y", namespace)))
    end if
 
    if (writ%out(out_e_field_surface_z)%write) then
      call write_iter_init(writ%out(out_e_field_surface_z)%handle, first, &
        units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/electric_field_surface-z", namespace)))
    end if
 
    if (writ%out(out_b_field_surface_x)%write) then
      call write_iter_init(writ%out(out_b_field_surface_x)%handle, first, &
        units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/magnetic_field_surface-x", namespace)))
    end if
 
    if (writ%out(out_b_field_surface_y)%write) then
      call write_iter_init(writ%out(out_b_field_surface_y)%handle, first, &
        units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/magnetic_field_surface-y", namespace)))
    end if
 
    if (writ%out(out_b_field_surface_z)%write) then
      call write_iter_init(writ%out(out_b_field_surface_z)%handle, first, &
        units_from_atomic(units_out%time, dt), trim(io_workpath("td.general/magnetic_field_surface-z", namespace)))
    end if
 
    pop_sub(td_write_mxll_init)
  end subroutine td_write_mxll_init
 
 
  ! ---------------------------------------------------------
  subroutine td_write_mxll_end(writ)
    type(td_write_t), intent(inout) :: writ
 
    integer :: iout
 
    push_sub(td_write_mxll_end)
 
    do iout = 1, out_maxwell_max
      if (writ%out(iout)%write)  call write_iter_end(writ%out(iout)%handle)
    end do
 
    pop_sub(td_write_mxll_end)
  end subroutine td_write_mxll_end
 
 
  ! ---------------------------------------------------------
  subroutine td_write_mxll_iter(writ, space, gr, st, hm, helmholtz, dt, iter, namespace)
    type(td_write_t),                intent(inout) :: writ
    class(space_t),                  intent(in)    :: space
    type(grid_t),                    intent(inout) :: gr
    type(states_mxll_t),             intent(inout) :: st
    type(hamiltonian_mxll_t),        intent(inout) :: hm
    type(helmholtz_decomposition_t), intent(inout) :: helmholtz
    real(real64),                    intent(in)    :: dt
    integer,                         intent(in)    :: iter
    type(namespace_t),               intent(in)    :: namespace
 
 
    push_sub(td_write_mxll_iter)
 
    call profiling_in("TD_WRITE_ITER_MAXWELL")
 
    if (writ%out(out_maxwell_trans_e_field)%write .or. writ%out(out_maxwell_trans_b_field)%write) then
      call helmholtz%get_trans_field(namespace, st%rs_state_trans, total_field=st%rs_state)
      call get_rs_state_at_point(st%selected_points_rs_state_trans(:,:), st%rs_state_trans, &
        st%selected_points_coordinate(:,:), st, gr)
      !TODO: calculate transverse energy
    else
      hm%energy%energy_trans = m_zero
    end if
 
    if (writ%out(out_maxwell_long_e_field)%write .or. writ%out(out_maxwell_long_b_field)%write) then
      call helmholtz%get_long_field(namespace, st%rs_state_long, total_field=st%rs_state)
      call get_rs_state_at_point(st%selected_points_rs_state_long(:,:), st%rs_state_long, &
        st%selected_points_coordinate(:,:), st, gr)
      !TODO: calculate longitudinal energy
    else
      hm%energy%energy_long = m_zero
    end if
 
    ! total E field
    if (writ%out(out_maxwell_total_e_field)%write) then
      call td_write_fields(writ%out(out_maxwell_total_e_field)%handle, space, st, iter, dt, &
        maxwell_e_field, maxwell_total_field, 1)
      call td_write_fields(writ%out(out_maxwell_total_e_field + 1)%handle, space, st, iter, dt, &
        maxwell_e_field, maxwell_total_field, 2)
      call td_write_fields(writ%out(out_maxwell_total_e_field + 2)%handle, space, st, iter, dt, &
        maxwell_e_field, maxwell_total_field, 3)
    end if
 
    ! total B field
    if (writ%out(out_maxwell_total_b_field)%write) then
      call td_write_fields(writ%out(out_maxwell_total_b_field)%handle, space, st, iter, dt, &
        maxwell_b_field, maxwell_total_field, 1)
      call td_write_fields(writ%out(out_maxwell_total_b_field + 1)%handle, space, st, iter, dt, &
        maxwell_b_field, maxwell_total_field, 2)
      call td_write_fields(writ%out(out_maxwell_total_b_field + 2)%handle, space, st, iter, dt, &
        maxwell_b_field, maxwell_total_field, 3)
    end if
 
    ! Longitudinal E field
    if (writ%out(out_maxwell_long_e_field)%write) then
      call td_write_fields(writ%out(out_maxwell_long_e_field)%handle, space, st, iter, dt, &
        maxwell_e_field, maxwell_long_field, 1)
      call td_write_fields(writ%out(out_maxwell_long_e_field + 1)%handle, space, st, iter, dt, &
        maxwell_e_field, maxwell_long_field, 2)
      call td_write_fields(writ%out(out_maxwell_long_e_field + 2)%handle, space, st, iter, dt, &
        maxwell_e_field, maxwell_long_field, 3)
    end if
 
    ! Longitudinal B field
    if (writ%out(out_maxwell_long_b_field)%write) then
      call td_write_fields(writ%out(out_maxwell_long_b_field)%handle, space, st, iter, dt, &
        maxwell_b_field, maxwell_long_field, 1)
      call td_write_fields(writ%out(out_maxwell_long_b_field + 1)%handle, space, st, iter, dt, &
        maxwell_b_field, maxwell_long_field, 2)
      call td_write_fields(writ%out(out_maxwell_long_b_field + 2)%handle, space, st, iter, dt, &
        maxwell_b_field, maxwell_long_field, 3)
    end if
 
    ! Transverse E field
    if (writ%out(out_maxwell_trans_e_field)%write) then
      call td_write_fields(writ%out(out_maxwell_trans_e_field)%handle, space, st, iter, dt, &
        maxwell_e_field, maxwell_trans_field, 1)
      call td_write_fields(writ%out(out_maxwell_trans_e_field + 1)%handle, space, st, iter, dt, &
        maxwell_e_field, maxwell_trans_field, 2)
      call td_write_fields(writ%out(out_maxwell_trans_e_field + 2)%handle, space, st, iter, dt, &
        maxwell_e_field, maxwell_trans_field, 3)
    end if
 
    ! Transverse B field
    if (writ%out(out_maxwell_trans_b_field)%write) then
      call td_write_fields(writ%out(out_maxwell_trans_b_field)%handle, space, st, iter, dt, &
        maxwell_b_field, maxwell_trans_field, 1)
      call td_write_fields(writ%out(out_maxwell_trans_b_field + 1)%handle, space, st, iter, dt, &
        maxwell_b_field, maxwell_trans_field, 2)
      call td_write_fields(writ%out(out_maxwell_trans_b_field + 2)%handle, space, st, iter, dt, &
        maxwell_b_field, maxwell_trans_field, 3)
    end if
 
    ! Maxwell energy
    if (writ%out(out_maxwell_energy)%write) then
      call td_write_maxwell_energy(writ%out(out_maxwell_energy)%handle, hm, iter)
    end if
 
    if (writ%out(out_e_field_surface_x)%write) then
      call td_write_electric_field_box_surface(writ%out(out_e_field_surface_x)%handle, st, 1, iter)
    end if
 
    if (writ%out(out_e_field_surface_y)%write) then
      call td_write_electric_field_box_surface(writ%out(out_e_field_surface_y)%handle, st, 2, iter)
    end if
 
    if (writ%out(out_e_field_surface_z)%write) then
      call td_write_electric_field_box_surface(writ%out(out_e_field_surface_z)%handle, st, 3, iter)
    end if
 
    if (writ%out(out_b_field_surface_x)%write) then
      call td_write_magnetic_field_box_surface(writ%out(out_b_field_surface_x)%handle, st, 1, iter)
    end if
 
    if (writ%out(out_b_field_surface_y)%write) then
      call td_write_magnetic_field_box_surface(writ%out(out_b_field_surface_y)%handle, st, 2, iter)
    end if
 
    if (writ%out(out_b_field_surface_z)%write) then
      call td_write_magnetic_field_box_surface(writ%out(out_b_field_surface_z)%handle, st, 3, iter)
    end if
 
    call profiling_out("TD_WRITE_ITER_MAXWELL")
 
    pop_sub(td_write_mxll_iter)
  end subroutine td_write_mxll_iter
 
 
  ! ---------------------------------------------------------
  subroutine td_write_maxwell_energy(out_maxwell_energy, hm, iter)
    type(c_ptr),                   intent(inout) :: out_maxwell_energy
    type(hamiltonian_mxll_t),      intent(in)    :: hm
    integer,                       intent(in)    :: iter
 
    integer :: ii
 
    integer :: n_columns
 
    if (.not. mpi_grp_is_root(mpi_world)) return ! only first node outputs
 
    push_sub(td_write_maxwell_energy)
 
    n_columns = 8
 
    if (iter == 0) then
      call td_write_print_header_init(out_maxwell_energy)
 
      ! first line -> column names
      call write_iter_header_start(out_maxwell_energy)
      call write_iter_header(out_maxwell_energy, 'Total')
      call write_iter_header(out_maxwell_energy, 'E**2')
      call write_iter_header(out_maxwell_energy, 'B**2')
      call write_iter_header(out_maxwell_energy, 'Total+Boundaries')
      call write_iter_header(out_maxwell_energy, 'Boundaries')
      call write_iter_header(out_maxwell_energy, 'Transversal')
      call write_iter_header(out_maxwell_energy, 'Longitudinal')
      call write_iter_header(out_maxwell_energy, 'Incident Waves')
 
      call write_iter_nl(out_maxwell_energy)
 
      ! units
 
      call write_iter_string(out_maxwell_energy, '#[Iter n.]')
      call write_iter_header(out_maxwell_energy, '[' // trim(units_abbrev(units_out%time)) // ']')
 
      do ii = 1, n_columns
        call write_iter_header(out_maxwell_energy, '[' // trim(units_abbrev(units_out%energy)) // ']')
      end do
      call write_iter_nl(out_maxwell_energy)
 
      call td_write_print_header_end(out_maxwell_energy)
    end if
 
    call write_iter_start(out_maxwell_energy)
    call write_iter_double(out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%energy), 1)
    call write_iter_double(out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%e_energy), 1)
    call write_iter_double(out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%b_energy), 1)
    call write_iter_double(out_maxwell_energy, units_from_atomic(units_out%energy, &
      hm%energy%energy+hm%energy%boundaries), 1)
    call write_iter_double(out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%boundaries), 1)
    call write_iter_double(out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%energy_trans), 1)
    call write_iter_double(out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%energy_long), 1)
    call write_iter_double(out_maxwell_energy, units_from_atomic(units_out%energy, hm%energy%energy_plane_waves), 1)
    call write_iter_nl(out_maxwell_energy)
 
    pop_sub(td_write_maxwell_energy)
  end subroutine td_write_maxwell_energy
 
 
  ! ---------------------------------------------------------
  subroutine td_write_electric_field_box_surface(out_field_surf, st, dim, iter)
    type(c_ptr),                   intent(inout) :: out_field_surf
    type(states_mxll_t),           intent(in)    :: st
    integer,                       intent(in)    :: dim
    integer,                       intent(in)    :: iter
 
    integer :: ii
 
    integer :: n_columns
 
    if (.not. mpi_grp_is_root(mpi_world)) return ! only first node outputs
 
    push_sub(td_write_electric_field_box_surface)
 
    n_columns = 12
 
    if (iter == 0) then
      call td_write_print_header_init(out_field_surf)
 
      ! first line -> column names
      call write_iter_header_start(out_field_surf)
      call write_iter_header(out_field_surf, '- x direction')
      call write_iter_header(out_field_surf, '+ x direction')
      call write_iter_header(out_field_surf, '- y direction')
      call write_iter_header(out_field_surf, '+ y direction')
      call write_iter_header(out_field_surf, '- z direction')
      call write_iter_header(out_field_surf, '+ z direction')
      call write_iter_header(out_field_surf, '- x dir. p. w.')
      call write_iter_header(out_field_surf, '+ x dir. p. w.')
      call write_iter_header(out_field_surf, '- y dir. p. w.')
      call write_iter_header(out_field_surf, '+ y dir. p. w.')
      call write_iter_header(out_field_surf, '- z dir. p. w.')
      call write_iter_header(out_field_surf, '+ z dir. p. w.')
 
      call write_iter_nl(out_field_surf)
 
      ! units
      call write_iter_string(out_field_surf, '#[Iter n.]')
      call write_iter_header(out_field_surf, '[' // trim(units_abbrev(units_out%time)) // ']')
 
      do ii = 1, n_columns
        call write_iter_header(out_field_surf, '[' // trim(units_abbrev(units_out%energy/units_out%length)) // ']')
      end do
      call write_iter_nl(out_field_surf)
 
      call td_write_print_header_end(out_field_surf)
    end if
 
    call write_iter_start(out_field_surf)
    call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
      st%electric_field_box_surface(1,1,dim)), 1)
    call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
      st%electric_field_box_surface(2,1,dim)), 1)
    call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
      st%electric_field_box_surface(1,2,dim)), 1)
    call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
      st%electric_field_box_surface(2,2,dim)), 1)
    call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
      st%electric_field_box_surface(1,3,dim)), 1)
    call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
      st%electric_field_box_surface(2,3,dim)), 1)
    call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
      st%electric_field_box_surface_plane_waves(1,1,dim)), 1)
    call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
      st%electric_field_box_surface_plane_waves(2,1,dim)), 1)
    call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
      st%electric_field_box_surface_plane_waves(1,2,dim)), 1)
    call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
      st%electric_field_box_surface_plane_waves(2,2,dim)), 1)
    call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
      st%electric_field_box_surface_plane_waves(1,3,dim)), 1)
    call write_iter_double(out_field_surf, units_from_atomic(units_out%energy/units_out%length, &
      st%electric_field_box_surface_plane_waves(2,3,dim)), 1)
    call write_iter_nl(out_field_surf)
 
    pop_sub(td_write_electric_field_box_surface)
  end subroutine td_write_electric_field_box_surface
 
 
  ! ---------------------------------------------------------
  subroutine td_write_magnetic_field_box_surface(out_field_surf, st, dim, iter)
    type(c_ptr),                   intent(inout) :: out_field_surf
    type(states_mxll_t),           intent(in)    :: st
    integer,                       intent(in)    :: dim
    integer,                       intent(in)    :: iter
 
    integer :: ii
 
    integer :: n_columns
 
    if (.not. mpi_grp_is_root(mpi_world)) return ! only first node outputs
 
    push_sub(td_write_magnetic_field_box_surface)
 
    n_columns = 12
 
    if (iter == 0) then
      call td_write_print_header_init(out_field_surf)
 
      ! first line -> column names
      call write_iter_header_start(out_field_surf)
      call write_iter_header(out_field_surf, '- x direction')
      call write_iter_header(out_field_surf, '+ x direction')
      call write_iter_header(out_field_surf, '- y direction')
      call write_iter_header(out_field_surf, '+ y direction')
      call write_iter_header(out_field_surf, '- z direction')
      call write_iter_header(out_field_surf, '+ z direction')
      call write_iter_header(out_field_surf, '- x dir. p. w.')
      call write_iter_header(out_field_surf, '+ x dir. p. w.')
      call write_iter_header(out_field_surf, '- y dir. p. w.')
      call write_iter_header(out_field_surf, '+ y dir. p. w.')
      call write_iter_header(out_field_surf, '- z dir. p. w.')
      call write_iter_header(out_field_surf, '+ z dir. p. w.')
 
      call write_iter_nl(out_field_surf)
 
      ! units
      call write_iter_string(out_field_surf, '#[Iter n.]')
      call write_iter_header(out_field_surf, '[' // trim(units_abbrev(units_out%time)) // ']')
 
      do ii = 1, n_columns
        call write_iter_header(out_field_surf, '[' // trim(units_abbrev(unit_one/units_out%length**2)) // ']')
      end do
      call write_iter_nl(out_field_surf)
 
      call td_write_print_header_end(out_field_surf)
    end if
 
    call write_iter_start(out_field_surf)
    call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
      st%magnetic_field_box_surface(1,1,dim)), 1)
    call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
      st%magnetic_field_box_surface(2,1,dim)), 1)
    call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
      st%magnetic_field_box_surface(1,2,dim)), 1)
    call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
      st%magnetic_field_box_surface(2,2,dim)), 1)
    call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
      st%magnetic_field_box_surface(1,3,dim)), 1)
    call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
      st%magnetic_field_box_surface(2,3,dim)), 1)
    call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
      st%magnetic_field_box_surface_plane_waves(1,1,dim)), 1)
    call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
      st%magnetic_field_box_surface_plane_waves(2,1,dim)), 1)
    call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
      st%magnetic_field_box_surface_plane_waves(1,2,dim)), 1)
    call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
      st%magnetic_field_box_surface_plane_waves(2,2,dim)), 1)
    call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
      st%magnetic_field_box_surface_plane_waves(1,3,dim)), 1)
    call write_iter_double(out_field_surf, units_from_atomic(unit_one/units_out%length**2, &
      st%magnetic_field_box_surface_plane_waves(2,3,dim)), 1)
    call write_iter_nl(out_field_surf)
 
    pop_sub(td_write_magnetic_field_box_surface)
  end subroutine td_write_magnetic_field_box_surface
 
  ! ---------------------------------------------------------
  subroutine td_write_fields(out_fields, space, st, iter, dt, e_or_b_field, field_type, idir)
    type(c_ptr),         intent(inout) :: out_fields
    class(space_t),      intent(in)    :: space
    type(states_mxll_t), intent(in)    :: st
    integer,             intent(in)    :: iter
    real(real64),        intent(in)    :: dt
    integer,             intent(in)    :: e_or_b_field
    integer,             intent(in)    :: field_type
    integer,             intent(in)    :: idir
 
 
    integer :: id
    real(real64) :: field(space%dim), selected_field
    character(len=80) :: aux
 
    if (.not. mpi_grp_is_root(mpi_world)) return ! only first node outputs
 
    push_sub(td_write_fields)
 
    if (iter == 0) then
      call td_write_print_header_init(out_fields)
 
      ! first line
      write(aux, '(a7,e20.12,3a)') '# dt = ', units_from_atomic(units_out%time, dt), &
        " [", trim(units_abbrev(units_out%time)), "]"
      call write_iter_string(out_fields, aux)
      call write_iter_nl(out_fields)
 
      call write_iter_header_start(out_fields)
 
      do id = 1, st%selected_points_number
        select case (e_or_b_field)
        case (maxwell_e_field)
          write(aux, '(a,i1,a)') 'E(', id, ')'
        case (maxwell_b_field)
          write(aux, '(a,i1,a)') 'B(', id, ')'
        end select
        call write_iter_header(out_fields, aux)
      end do
 
      call write_iter_nl(out_fields)
      call write_iter_string(out_fields, '#[Iter n.]')
      call write_iter_header(out_fields, '          [' // trim(units_abbrev(units_out%time)) // ']')
 
      ! Note that we do not print out units of E, B, or A, but rather units of e*E, e*B, e*A.
      ! (force, force, and energy, respectively). The reason is that the units of E, B or A
      ! are ugly.
      aux = '          [' // trim(units_abbrev(units_out%force)) // ']'
      do id = 1, st%selected_points_number
        call write_iter_header(out_fields, aux)
      end do
      call write_iter_nl(out_fields)
      call td_write_print_header_end(out_fields)
    end if
 
    call write_iter_start(out_fields)
 
    do id = 1, st%selected_points_number
      select case (e_or_b_field)
      case (maxwell_e_field)
        ! Output of electric field at selected point
        select case (field_type)
        case (maxwell_total_field)
          call get_electric_field_vector(st%selected_points_rs_state(:,id), field(1:st%dim))
        case (maxwell_long_field)
          call get_electric_field_vector(st%selected_points_rs_state_long(:,id), field(1:st%dim))
        case (maxwell_trans_field)
          call get_electric_field_vector(st%selected_points_rs_state_trans(:,id), field(1:st%dim))
        end select
        selected_field = units_from_atomic(units_out%energy/units_out%length, field(idir))
      case (maxwell_b_field)
        ! Output of magnetic field at selected point
        select case (field_type)
        case (maxwell_total_field)
          call get_magnetic_field_vector(st%selected_points_rs_state(:,id), st%rs_sign, field(1:st%dim))
        case (maxwell_long_field)
          call get_magnetic_field_vector(st%selected_points_rs_state_long(:,id), st%rs_sign, field(1:st%dim))
        case (maxwell_trans_field)
          call get_magnetic_field_vector(st%selected_points_rs_state_trans(:,id), st%rs_sign, field(1:st%dim))
        end select
        selected_field = units_from_atomic(unit_one/units_out%length**2, field(idir))
      end select
      call write_iter_double(out_fields, selected_field, 1)
    end do
 
    call write_iter_nl(out_fields)
 
    pop_sub(td_write_fields)
  end subroutine td_write_fields
 
 
  !----------------------------------------------------------
  subroutine td_write_mxll_free_data(writ, namespace, space, gr, st, hm, helmholtz, outp, iter, time)
    type(td_write_t),                intent(inout) :: writ
    type(namespace_t),               intent(in)    :: namespace
    class(space_t),                  intent(in)    :: space
    type(grid_t),                    intent(inout) :: gr
    type(states_mxll_t),             intent(inout) :: st
    type(hamiltonian_mxll_t),        intent(inout) :: hm
    type(helmholtz_decomposition_t), intent(inout) :: helmholtz
    type(output_t),                  intent(in)    :: outp
    integer,                         intent(in)    :: iter
    real(real64),                    intent(in)    :: time
 
    character(len=256) :: filename
 
    push_sub(td_write_maxwell_free_data)
    call profiling_in("TD_WRITE_MAXWELL_DATA")
 
    ! now write down the rest
    write(filename, '(a,a,i7.7)') trim(outp%iter_dir),"td.", iter  ! name of directory
 
    call output_mxll(outp, namespace, space, gr, st, hm, helmholtz, time, filename)
 
    call profiling_out("TD_WRITE_MAXWELL_DATA")
    pop_sub(td_write_maxwell_free_data)
  end subroutine td_write_mxll_free_data
 
end module td_write_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: