spectrum.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 spectrum_oct_m
  use batch_oct_m
  use iso_c_binding
  use compressed_sensing_oct_m
  use debug_oct_m
  use fft_oct_m
  use global_oct_m
  use io_oct_m
  use kick_oct_m
  use, intrinsic :: iso_fortran_env
  use lalg_adv_oct_m
  use math_oct_m
  use messages_oct_m
  use minimizer_oct_m
  use namespace_oct_m
  use parser_oct_m
  use pcm_oct_m
  use profiling_oct_m
  use string_oct_m
  use types_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use varinfo_oct_m
 
  implicit none
 
  private
  public ::                        &
    spectrum_t,                    &
    spectrum_init,                 &
    spectrum_cross_section,        &
    spectrum_cross_section_tensor, &
    spectrum_dipole_power,         &
    spectrum_dyn_structure_factor, &
    spectrum_rotatory_strength,    &
    spectrum_hs_from_mult,         &
    spectrum_hs_ar_from_mult,      &
    spectrum_hs_ar_from_acc,       &
    spectrum_hs_from_acc,          &
    spectrum_hsfunction_init,      &
    spectrum_hsfunction_end,       &
    spectrum_hsfunction_min,       &
    spectrum_mult_info,            &
    spectrum_fix_time_limits,      &
    spectrum_count_time_steps,     &
    spectrum_signal_damp,          &
    spectrum_fourier_transform,    &
    spectrum_hs_from_current,      &
    spectrum_nenergy_steps
 
  integer, public, parameter ::    &
    SPECTRUM_DAMP_NONE       = 0,  &
    spectrum_damp_lorentzian = 1,  &
    spectrum_damp_polynomial = 2,  &
    spectrum_damp_gaussian   = 3,  &
    spectrum_damp_sin        = 4
 
  integer, public, parameter ::    &
    SPECTRUM_TRANSFORM_LAPLACE = 1,  &
    spectrum_transform_sin     = 2,  &
    spectrum_transform_cos     = 3
 
  integer, public, parameter :: &
    SPECTRUM_ABSORPTION    = 1, &
    spectrum_energyloss    = 2, &
    spectrum_p_power       = 3, &
    spectrum_rotatory      = 4
 
  integer, public, parameter ::       &
    SPECTRUM_FOURIER            = 1,  &
    spectrum_compressed_sensing = 2
 
  type spectrum_t
    real(real64)   :: start_time
    real(real64)   :: end_time
    real(real64)   :: energy_step
    real(real64)   :: min_energy
    real(real64)   :: max_energy
    integer :: damp
    integer :: transform
    real(real64)   :: damp_factor
    integer :: spectype
    integer :: method
    real(real64)   :: noise
    logical, private :: sigma_diag
  end type spectrum_t
 
  integer :: niter_
  real(real64) :: time_step_, energy_step_
  complex(real64), allocatable :: func_(:),func_ar_(:,:),pos_(:,:),tret_(:), funcw_(:)
  type(fft_t), save :: fft_handler
  integer :: is_, ie_, default
 
contains
 
  ! ---------------------------------------------------------
  subroutine spectrum_init(spectrum, namespace, default_energy_step, default_max_energy)
    type(spectrum_t),           intent(inout) :: spectrum
    type(namespace_t),          intent(in)    :: namespace
    real(real64),     optional, intent(in)    :: default_energy_step
    real(real64),     optional, intent(in)    :: default_max_energy
 
    real(real64) :: fdefault
 
    push_sub(spectrum_init)
 
    call messages_print_with_emphasis(msg="Spectrum Options", namespace=namespace)
 
    !%Variable PropagationSpectrumType
    !%Type integer
    !%Default AbsorptionSpectrum
    !%Section Utilities::oct-propagation_spectrum
    !%Description
    !% Type of spectrum to calculate.
    !%Option AbsorptionSpectrum 1
    !% Photoabsorption spectrum.
    !%Option EnergyLossSpectrum 2
    !% Dynamic structure factor (also known as energy-loss function or spectrum).
    !%Option DipolePower 3
    !% Power spectrum of the dipole moment.
    !%Option RotatoryStrength 4
    !% Rotatory strength spectrum.
    !%End
 
    call parse_variable(namespace, 'PropagationSpectrumType', spectrum_absorption, spectrum%spectype)
 
    if (.not. varinfo_valid_option('PropagationSpectrumType', spectrum%spectype)) then
      call messages_input_error(namespace, 'PropagationSpectrumType')
    end if
    call messages_print_var_option('PropagationSpectrumType', spectrum%spectype, namespace=namespace)
 
    !%Variable SpectrumMethod
    !%Type integer
    !%Default fourier
    !%Section Utilities::oct-propagation_spectrum
    !%Description
    !% Decides which method is used to obtain the spectrum.
    !%Option fourier 1
    !% The standard Fourier transform. Further specified by <tt>PropagationSpectrumTransform</tt>.
    !%Option compressed_sensing 2
    !% (Experimental) Uses the compressed sensing technique.
    !%End
    call parse_variable(namespace, 'SpectrumMethod', spectrum_fourier, spectrum%method)
    if (.not. varinfo_valid_option('SpectrumMethod', spectrum%method)) then
      call messages_input_error(namespace, 'SpectrumMethod')
    end if
    call messages_print_var_option('SpectrumMethod', spectrum%method, namespace=namespace)
 
    if (spectrum%method == spectrum_compressed_sensing) then
      call messages_experimental('compressed sensing', namespace=namespace)
 
      !%Variable SpectrumSignalNoise
      !%Type float
      !%Default 0.0
      !%Section Utilities::oct-propagation_spectrum
      !%Description
      !% For compressed sensing, the signal to process, the
      !% time-dependent dipole in this case, is assumed to have some
      !% noise that is given by this dimensionless quantity.
      !%End
      call parse_variable(namespace, 'SpectrumSignalNoise', 0.0_real64, spectrum%noise)
      call messages_print_var_value('SpectrumSignalNoise', spectrum%noise, namespace=namespace)
    end if
 
 
    !%Variable PropagationSpectrumDampMode
    !%Type integer
    !%Section Utilities::oct-propagation_spectrum
    !%Description
    !% Decides which damping/filtering is to be applied in order to
    !% calculate spectra by calculating a Fourier transform. The
    !% default is polynomial damping, except when <tt>SpectrumMethod = compressed_sensing</tt>.
    !% In that case the default is none.
    !%Option none 0
    !% No filtering at all.
    !%Option exponential 1
    !% Exponential filtering, corresponding to a Lorentzian-shaped spectrum.
    !%Option polynomial 2
    !% Third-order polynomial damping.
    !%Option gaussian 3
    !% Gaussian damping.
    !%End
    default = spectrum_damp_polynomial
    if (spectrum%method == spectrum_compressed_sensing) default = spectrum_damp_none
 
    call parse_variable(namespace, 'PropagationSpectrumDampMode', default, spectrum%damp)
 
    if (.not. varinfo_valid_option('PropagationSpectrumDampMode', spectrum%damp)) then
      call messages_input_error(namespace, 'PropagationSpectrumDampMode')
    end if
    call messages_print_var_option('PropagationSpectrumDampMode', spectrum%damp, namespace=namespace)
 
    if (spectrum%method == spectrum_compressed_sensing .and. spectrum%damp /= spectrum_damp_none) then
      message(1) = 'Using damping with compressed sensing, this is not required'
      message(2) = 'and can introduce noise in the spectra.'
      call messages_warning(2, namespace=namespace)
    end if
 
    !%Variable PropagationSpectrumTransform
    !%Type integer
    !%Default sine
    !%Section Utilities::oct-propagation_spectrum
    !%Description
    !% Decides which transform to perform, if <tt>SpectrumMethod = fourier</tt>.
    !%Option sine 2
    !% Sine transform: <math>\int dt \sin(wt) f(t)</math>. Produces the imaginary part of the polarizability.
    !%Option cosine 3
    !% Cosine transform: <math>\int dt \cos(wt) f(t)</math>. Produces the real part of the polarizability.
    !%Option laplace 1
    !% Real exponential transform: <math>\int dt e^{-wt} f(t)</math>. Produces the real part of the polarizability at imaginary
    !% frequencies, <i>e.g.</i> for Van der Waals <math>C_6</math> coefficients.
    !% This is the only allowed choice for complex scaling.
    !%End
    call parse_variable(namespace, 'PropagationSpectrumTransform', spectrum_transform_sin, spectrum%transform)
    if (.not. varinfo_valid_option('PropagationSpectrumTransform', spectrum%transform)) then
      call messages_input_error(namespace, 'PropagationSpectrumTransform')
    end if
    call messages_print_var_option('PropagationSpectrumTransform', spectrum%transform, namespace=namespace)
 
    !%Variable PropagationSpectrumStartTime
    !%Type float
    !%Default 0.0
    !%Section Utilities::oct-propagation_spectrum
    !%Description
    !% Processing is done for the given function in a time-window that starts at the
    !% value of this variable.
    !%End
    call parse_variable(namespace, 'PropagationSpectrumStartTime',  m_zero, spectrum%start_time, units_inp%time)
    call messages_print_var_value('PropagationSpectrumStartTime', spectrum%start_time, unit = units_out%time, namespace=namespace)
 
    !%Variable PropagationSpectrumEndTime
    !%Type float
    !%Default -1.0 au
    !%Section Utilities::oct-propagation_spectrum
    !%Description
    !% Processing is done for the given function in a time-window that ends at the
    !% value of this variable. If set to a negative value, the maximum value from
    !% the corresponding multipole file will used.
    !%End
    call parse_variable(namespace, 'PropagationSpectrumEndTime', -m_one, spectrum%end_time, units_inp%time)
    call messages_print_var_value('PropagationSpectrumEndTime', spectrum%end_time, unit = units_out%time, namespace=namespace)
 
    !%Variable PropagationSpectrumEnergyStep
    !%Type float
    !%Default 0.01 eV
    !%Section Utilities::oct-propagation_spectrum
    !%Description
    !% Sampling rate for the spectrum. If you supply a number equal or smaller than zero, then
    !% the sampling rate will be <math>2 \pi / T</math>, where <math>T</math> is the total propagation time.
    !%End
    fdefault = 0.01_real64/(m_two*p_ry)
    if (present(default_energy_step)) fdefault = default_energy_step
    call parse_variable(namespace, 'PropagationSpectrumEnergyStep', fdefault, spectrum%energy_step, units_inp%energy)
    call messages_print_var_value('PropagationSpectrumEnergyStep', spectrum%energy_step, unit = units_out%energy, &
      namespace=namespace)
 
    !%Variable PropagationSpectrumMinEnergy
    !%Type float
    !%Default 0
    !%Section Utilities::oct-propagation_spectrum
    !%Description
    !% The Fourier transform is calculated for energies larger than this value.
    !%End
    call parse_variable(namespace, 'PropagationSpectrumMinEnergy', m_zero, spectrum%min_energy, units_inp%energy)
    call messages_print_var_value('PropagationSpectrumMinEnergy', spectrum%min_energy, unit = units_out%energy, &
      namespace=namespace)
 
 
    !%Variable PropagationSpectrumMaxEnergy
    !%Type float
    !%Default 20 eV
    !%Section Utilities::oct-propagation_spectrum
    !%Description
    !% The Fourier transform is calculated for energies smaller than this value.
    !%End
    fdefault = 20.0_real64/(m_two*p_ry)
    if (present(default_max_energy)) fdefault = default_max_energy
    call parse_variable(namespace, 'PropagationSpectrumMaxEnergy', fdefault, spectrum%max_energy, units_inp%energy)
    call messages_print_var_value('PropagationSpectrumMaxEnergy', spectrum%max_energy, unit = units_out%energy, &
      namespace=namespace)
 
    !%Variable PropagationSpectrumDampFactor
    !%Type float
    !%Default -1.0
    !%Section Utilities::oct-propagation_spectrum
    !%Description
    !% If <tt>PropagationSpectrumDampMode = exponential, gaussian</tt>, the damping parameter of the exponential
    !% is fixed through this variable.
    !% Default value ensure that the damping function adquires a 0.0001 value at the end of the propagation time.
    !%End
    call parse_variable(namespace, 'PropagationSpectrumDampFactor', -m_one, spectrum%damp_factor, units_inp%time**(-1))
 
    call messages_print_var_value('PropagationSpectrumDampFactor', spectrum%damp_factor, unit = units_out%time**(-1), &
      namespace=namespace)
 
    !%Variable PropagationSpectrumSigmaDiagonalization
    !%Type logical
    !%Default .false.
    !%Section Utilities::oct-propagation_spectrum
    !%Description
    !% If <tt>PropagationSpectrumSigmaDiagonalization = yes</tt>, the polarizability tensor is diagonalizied.
    !% This variable is only used if the cross_section_tensor is computed.
    !%End
    call parse_variable(namespace, 'PropagationSpectrumSigmaDiagonalization', .false., spectrum%sigma_diag)
    call messages_print_var_value('PropagationSpectrumSigmaDiagonalization', spectrum%sigma_diag, namespace=namespace)
 
    call messages_print_with_emphasis(namespace=namespace)
 
    pop_sub(spectrum_init)
  end subroutine spectrum_init
 
 
  ! ---------------------------------------------------------
  subroutine spectrum_cross_section_tensor(spectrum, namespace, out_file, in_file)
    type(spectrum_t),  intent(inout) :: spectrum
    type(namespace_t), intent(in)    :: namespace
    integer,           intent(in)    :: out_file
    integer,           intent(in)    :: in_file(:)
 
    integer :: nspin, energy_steps, ie, is, equiv_axes, n_files, trash
    real(real64), allocatable :: sigma(:, :, :, :), sigmap(:, :, :, :), sigmau(:, :, :),  &
      sigmav(:, :, :), sigmaw(:, :, :), ip(:, :)
    real(real64) :: dw, dump
    type(kick_t) :: kick
 
    push_sub(spectrum_cross_section_tensor)
 
    n_files = size(in_file)
    equiv_axes = 3 - n_files + 1
 
    call spectrum_cross_section_info(namespace, in_file(1), nspin, kick, energy_steps, dw)
    ! on subsequent calls, do not overwrite energy_steps and dw
    call io_skip_header(in_file(1))
 
    safe_allocate(sigma(1:3, 1:3, 1:energy_steps, 1:nspin))
    safe_allocate(sigmap(1:3, 1:3, 1:energy_steps, 1:nspin))
    safe_allocate(sigmau(1:3,      1:energy_steps, 1:nspin))
    safe_allocate(sigmav(1:3,      1:energy_steps, 1:nspin))
    safe_allocate(sigmaw(1:3,      1:energy_steps, 1:nspin))
    safe_allocate(    ip(1:3, 1:3))
 
    select case (equiv_axes)
 
    case (3)
 
      do ie = 1, energy_steps
        read(in_file(1), *) dump, (sigmau(1:3, ie, is), is = 1, nspin)
      end do
 
      ! The first row of sigma is the vector that we have just read, but properly projected...
      do is = 1, nspin
        do ie = 1, energy_steps
          sigmap(1, 1, ie, is) = sum(sigmau(1:3, ie, is) * kick%pol(1:3, 1))
          sigmap(1, 2, ie, is) = sum(sigmau(1:3, ie, is) * kick%pol(1:3, 2))
          sigmap(1, 3, ie, is) = sum(sigmau(1:3, ie, is) * kick%pol(1:3, 3))
        end do
      end do
 
      ! The diagonal parts are also equal:
      sigmap(2, 2, :, :) = sigmap(1, 1, :, :)
      sigmap(3, 3, :, :) = sigmap(1, 1, :, :)
 
      ! The (2,1) term and (3,1) term are equal by symmetry:
      sigmap(2, 1, :, :) = sigmap(1, 2, :, :)
      sigmap(3, 1, :, :) = sigmap(1, 3, :, :)
 
      ! But for the (2,3) term we need the wprime vector....
      do is = 1, nspin
        do ie = 1, energy_steps
          sigmap(2, 3, ie, is) = sum(sigmau(1:3, ie, is) * kick%wprime(1:3))
          sigmap(3, 2, ie, is) = sigmap(2, 3, ie, is)
        end do
      end do
 
    case (2)
 
      call spectrum_cross_section_info(namespace, in_file(2), ie, kick, trash, dump)
      call io_skip_header(in_file(2))
 
      do ie = 1, energy_steps
        read(in_file(1), *) dump, (sigmau(1:3, ie, is), is = 1, nspin)
        read(in_file(2), *) dump, (sigmaw(1:3, ie, is), is = 1, nspin)
      end do
 
      ! The first row of sigma is the vector that we have just read, but properly projected...
      do is = 1, nspin
        do ie = 1, energy_steps
          sigmap(1, 1, ie, is) = sum(sigmau(1:3, ie, is) * kick%pol(1:3, 1))
          sigmap(1, 2, ie, is) = sum(sigmau(1:3, ie, is) * kick%pol(1:3, 2))
          sigmap(1, 3, ie, is) = sum(sigmau(1:3, ie, is) * kick%pol(1:3, 3))
        end do
      end do
 
      ! The third row of sigma is also the vector that we have just read, but properly projected...
      do is = 1, nspin
        do ie = 1, energy_steps
          sigmap(3, 1, ie, is) = sum(sigmaw(1:3, ie, is) * kick%pol(1:3, 1))
          sigmap(3, 2, ie, is) = sum(sigmaw(1:3, ie, is) * kick%pol(1:3, 2))
          sigmap(3, 3, ie, is) = sum(sigmaw(1:3, ie, is) * kick%pol(1:3, 3))
        end do
      end do
 
      ! The diagonal (2,2) is equal by symmetry to (1,1)
      sigmap(2, 2, :, :) = sigmap(1, 1, :, :)
 
      ! The (2,1) term and (1,2) term are equal; (2,3) and (3,2), also.
      sigmap(2, 1, :, :) = sigmap(1, 2, :, :)
      sigmap(2, 3, :, :) = sigmap(3, 2, :, :)
 
    case default
 
      call spectrum_cross_section_info(namespace, in_file(2), ie, kick, trash, dump)
      call spectrum_cross_section_info(namespace, in_file(3), ie, kick, trash, dump)
      call io_skip_header(in_file(2))
      call io_skip_header(in_file(3))
 
      do ie = 1, energy_steps
        read(in_file(1), *) dump, (sigmau(1:3, ie, is), is = 1, nspin)
        read(in_file(2), *) dump, (sigmav(1:3, ie, is), is = 1, nspin)
        read(in_file(3), *) dump, (sigmaw(1:3, ie, is), is = 1, nspin)
      end do
 
      do is = 1, nspin
        do ie = 1, energy_steps
          sigmap(1, 1, ie, is) = sum(sigmau(1:3, ie, is) * kick%pol(1:3, 1))
          sigmap(1, 2, ie, is) = sum(sigmau(1:3, ie, is) * kick%pol(1:3, 2))
          sigmap(1, 3, ie, is) = sum(sigmau(1:3, ie, is) * kick%pol(1:3, 3))
        end do
      end do
      do is = 1, nspin
        do ie = 1, energy_steps
          sigmap(2, 1, ie, is) = sum(sigmav(1:3, ie, is) * kick%pol(1:3, 1))
          sigmap(2, 2, ie, is) = sum(sigmav(1:3, ie, is) * kick%pol(1:3, 2))
          sigmap(2, 3, ie, is) = sum(sigmav(1:3, ie, is) * kick%pol(1:3, 3))
        end do
      end do
      do is = 1, nspin
        do ie = 1, energy_steps
          sigmap(3, 1, ie, is) = sum(sigmaw(1:3, ie, is) * kick%pol(1:3, 1))
          sigmap(3, 2, ie, is) = sum(sigmaw(1:3, ie, is) * kick%pol(1:3, 2))
          sigmap(3, 3, ie, is) = sum(sigmaw(1:3, ie, is) * kick%pol(1:3, 3))
        end do
      end do
 
    end select
 
    ! And now, perform the necessary transformation.
    ip(1:3, 1:3) = kick%pol(1:3, 1:3)
    call lalg_inverse(3, ip, 'dir')
    do is = 1, nspin
      do ie = 1, energy_steps
        sigma(:, :, ie, is) = matmul(transpose(ip), matmul(sigmap(:, :, ie, is), ip))
      end do
    end do
 
    ! Finally, write down the result
    call spectrum_cross_section_tensor_write(out_file, sigma, nspin, spectrum%energy_step, &
      spectrum%min_energy, energy_steps, kick)
 
    ! Diagonalize sigma tensor
    if (spectrum%sigma_diag) then
      call spectrum_sigma_diagonalize(namespace, sigma, nspin, spectrum%energy_step, spectrum%min_energy, energy_steps, kick)
    end if
 
    safe_deallocate_a(sigma)
    safe_deallocate_a(sigmap)
    safe_deallocate_a(sigmau)
    safe_deallocate_a(sigmav)
    safe_deallocate_a(sigmaw)
    safe_deallocate_a(ip)
 
    pop_sub(spectrum_cross_section_tensor)
  end subroutine spectrum_cross_section_tensor
 
 
  ! ---------------------------------------------------------
  subroutine spectrum_cross_section_tensor_write(out_file, sigma, nspin, energy_step, min_energy, energy_steps, kick)
    integer,                intent(in) :: out_file
    real(real64),           intent(in) :: sigma(:, :, :, :)
    integer,                intent(in) :: nspin
    real(real64),           intent(in) :: energy_step, min_energy
    integer,                intent(in) :: energy_steps
    type(kick_t), optional, intent(in) :: kick
 
    integer :: is, idir, jdir, ie, ii
    real(real64) :: average, anisotropy
    real(real64), allocatable :: pp(:,:), pp2(:,:), ip(:,:)
    logical :: spins_singlet, spins_triplet
    character(len=20) :: header_string
 
    push_sub(spectrum_cross_section_tensor_write)
 
    spins_singlet = .true.
    spins_triplet = .false.
    if (present(kick)) then
      write(out_file, '(a15,i2)')      '# nspin        ', nspin
      call kick_write(kick, out_file)
      select case (kick_get_type(kick))
      case (kick_spin_mode)
        spins_triplet = .true.
        spins_singlet = .false.
      case (kick_spin_density_mode)
        spins_triplet = .true.
      end select
    end if
 
    write(out_file, '(a1, a20)', advance = 'no') '#', str_center("Energy", 20)
    write(out_file, '(a20)', advance = 'no') str_center("(1/3)*Tr[sigma]", 20)
    write(out_file, '(a20)', advance = 'no') str_center("Anisotropy[sigma]", 20)
    if (spins_triplet .and. spins_singlet) then
      write(out_file, '(a20)', advance = 'no') str_center("(1/3)*Tr[sigma-]", 20)
    end if
    do is = 1, nspin
      do idir = 1, 3
        do jdir = 1, 3
          write(header_string,'(a6,i1,a1,i1,a1,i1,a1)') 'sigma(', idir, ',', jdir, ',', is, ')'
          write(out_file, '(a20)', advance = 'no') str_center(trim(header_string), 20)
        end do
      end do
    end do
    write(out_file, '(1x)')
    write(out_file, '(a1,a20)', advance = 'no') '#', str_center('[' // trim(units_abbrev(units_out%energy)) // ']', 20)
    if (spins_triplet .and. spins_singlet) then
      write(out_file, '(a20)', advance = 'no')  str_center('[' // trim(units_abbrev(units_out%length**2)) // ']', 20)
    end if
    do ii = 1, 2 + nspin * 9
      write(out_file, '(a20)', advance = 'no')  str_center('[' // trim(units_abbrev(units_out%length**2)) // ']', 20)
    end do
    write(out_file, '(1x)')
 
    ! The anisotropy (Delta alpha) of a second-rank symmetric tensor alpha (such
    ! as the cross section) is defined as:
    !
    ! (Delta alpha)^2 = (1/3) * (3 * Tr(alpha^2) - (Tr(a))^2)
    !
    ! The reason for this definition is that it is identically equal to:
    !
    ! (Delta alpha)^2 = (1/3) * ( (alpha_1-alpha_2)^2 + (alpha_1-alpha_3)^2 + (alpha_2-alpha_3)^2)
    !
    ! where {alpha_1, alpha_2, alpha_3} are the eigenvalues of alpha. An "isotropic" tensor
    ! is characterized by having three equal eigenvalues, which leads to zero anisotropy. The
    ! more different that the eigenvalues are, the larger the anisotropy is.
 
    safe_allocate(pp(1:3, 1:3))
    if (spins_triplet .and. spins_singlet) then
      safe_allocate(pp2(1:3, 1:3))
    end if
    safe_allocate(ip(1:3, 1:3))
 
    do ie = 1, energy_steps
 
      pp(:, :) = sigma(:, :, ie, 1)
      if (nspin >= 2) then
        if (spins_singlet .and. spins_triplet) then
          pp2(:, :) = pp(:, :) - sigma(:, :, ie, 2)
          pp(:, :)  = pp(:, :) + sigma(:, :, ie, 2)
        elseif (spins_triplet .and. .not. spins_singlet) then
          pp(:, :) = pp(:, :) - sigma(:, :, ie, 2)
        elseif (spins_singlet .and. .not. spins_triplet) then
          pp(:, :) = pp(:, :) + sigma(:, :, ie, 2)
        end if
      end if
 
      average = m_third * (pp(1, 1) + pp(2, 2) + pp(3, 3))
      ip = matmul(pp, pp)
      anisotropy = m_third * (m_three * ( ip(1, 1) + ip(2, 2) + ip(3, 3)) - (m_three * average)**2)
 
      ! Note that the cross-section elements do not have to be transformed to the proper units, since
      ! they have been read from the "cross_section_vector.x", where they are already in the proper units.
      write(out_file,'(3e20.8)', advance = 'no') units_from_atomic(units_out%energy, ((ie - 1) * energy_step + min_energy)), &
        average, sqrt(max(anisotropy, m_zero))
 
      if (spins_singlet .and. spins_triplet) then
        average = m_third * (pp2(1, 1) + pp2(2, 2) + pp2(3, 3))
        write(out_file,'(1e20.8)', advance = 'no') average
      end if
 
      do is = 1, nspin
        write(out_file,'(9e20.8)', advance = 'no') sigma(1:3, 1:3, ie, is)
      end do
      write(out_file, '(1x)')
    end do
 
    safe_deallocate_a(pp)
    if (spins_triplet .and. spins_singlet) then
      safe_deallocate_a(pp2)
    end if
    safe_deallocate_a(ip)
    pop_sub(spectrum_cross_section_tensor_write)
  end subroutine spectrum_cross_section_tensor_write
 
 
  ! ---------------------------------------------------------
  subroutine spectrum_cross_section(spectrum, namespace, in_file, out_file, ref_file)
    type(spectrum_t),  intent(inout) :: spectrum
    type(namespace_t), intent(in)    :: namespace
    integer,           intent(in)    :: in_file
    integer,           intent(in)    :: out_file
    integer, optional, intent(in)    :: ref_file
 
    character(len=20) :: header_string
    integer :: nspin, ref_nspin, lmax, ref_lmax, time_steps, &
      ref_time_steps, istart, iend, ntiter, it, ii, isp, no_e, ie, idir
    real(real64)   :: dt, ref_dt, energy, ewsum, polsum
    type(kick_t) :: kick, ref_kick
    real(real64), allocatable :: dipole(:, :, :), ref_dipole(:, :, :), sigma(:, :, :), sf(:, :)
    type(unit_system_t) :: file_units, ref_file_units
    type(batch_t) :: dipoleb, sigmab
 
    type(pcm_min_t) :: pcm
 
    push_sub(spectrum_cross_section)
 
    ! This function gives us back the unit connected to the "multipoles" file, the header information,
    ! the number of time steps, and the time step.
    call spectrum_mult_info(namespace, in_file, nspin, kick, time_steps, dt, file_units, lmax=lmax)
 
    if (present(ref_file)) then
      call spectrum_mult_info(namespace, ref_file, ref_nspin, ref_kick, &
        ref_time_steps, ref_dt, ref_file_units, lmax = ref_lmax)
      if ((nspin /= ref_nspin)         .or. &
        (time_steps /= ref_time_steps) .or. &
        (.not.(abs(dt-ref_dt)< 1e-10_real64)) .or. &
        (lmax /= ref_lmax)) then
        write(message(1),'(a)') 'The multipoles and reference multipoles files do not match.'
        call messages_fatal(1, namespace=namespace)
      end if
    end if
 
    ! Now we cannot process files that do not contain the dipole, or that contain more than the dipole.
    if (lmax /= 1) then
      message(1) = 'Multipoles file should contain the dipole -- and only the dipole.'
      call messages_fatal(1, namespace=namespace)
    end if
 
    if (kick%function_mode /= kick_function_dipole) then
      message(1) = "Kick function must have been dipole to run this utility."
      call messages_fatal(1, namespace=namespace)
    end if
 
    if (kick%pol_dir < 1) then
      message(1) = "Kick polarization direction is not set. Probably no kick was used."
      call messages_fatal(1, namespace=namespace)
    end if
 
    ! Find out the iteration numbers corresponding to the time limits.
    call spectrum_fix_time_limits(spectrum, time_steps, dt, istart, iend, ntiter)
 
    safe_allocate(dipole(0:time_steps, 1:3, 1:nspin))
    call spectrum_read_dipole(namespace, in_file, dipole)
 
    if (present(ref_file)) then
      safe_allocate(ref_dipole(0:time_steps, 1:3, 1:nspin))
      call spectrum_read_dipole(namespace, ref_file, ref_dipole)
    end if
 
    ! parsing and re-printing to output useful PCM data
    call pcm_min_input_parsing_for_spectrum(pcm, namespace)
 
    ! adding the dipole generated by the PCM polarization charges due solute
    if (pcm%localf) then
      call spectrum_add_pcm_dipole(namespace, dipole, time_steps, nspin)
    end if
 
    ! Now subtract the initial dipole.
    if (present(ref_file)) then
      dipole = dipole - ref_dipole
    else
      do it = 1, time_steps
        dipole(it, :, :) = dipole(it, :, :) - dipole(0, :, :)
      end do
      dipole(0, :, :) = m_zero
    end if
 
    if (spectrum%energy_step <= m_zero) spectrum%energy_step = m_two * m_pi / (dt*time_steps)
 
    ! Get the number of energy steps.
    no_e = spectrum_nenergy_steps(spectrum)
    safe_allocate(sigma(1:no_e, 1:3, 1:nspin))
 
    if (pcm%localf) then
 
      ! for numerical reasons, we cannot add the difference (d(t)-d(t0)) of PCM dipoles here -- although it would look cleaner
 
      ! in the PCM local field case \sigma(\omega) \propto \Im{\alpha(\omega)\epsilon(\omega)} not just \propto \Im{\alpha(\omega)}
      ! since the dielectric function is complex as well, we need to compute both the real and imaginary part of the polarizability
      call spectrum_times_pcm_epsilon(spectrum, pcm, dipole, sigma, nspin, istart, iend, kick%time, dt, no_e)
 
      write(out_file,'(a57)') "Cross-section spectrum contains full local field effects."
 
    else
 
      call batch_init(dipoleb, 3, 1, nspin, dipole)
      call batch_init(sigmab, 3, 1, nspin, sigma)
 
      call spectrum_signal_damp(spectrum%damp, spectrum%damp_factor, istart + 1, iend + 1, kick%time, dt, dipoleb)
      call spectrum_fourier_transform(spectrum%method, spectrum%transform, spectrum%noise, &
        istart + 1, iend + 1, kick%time, dt, dipoleb, spectrum%min_energy, spectrum%max_energy, spectrum%energy_step, sigmab)
 
      call dipoleb%end()
      call sigmab%end()
 
    end if
 
    if (pcm%run_pcm) then
      call spectrum_over_pcm_refraction_index(spectrum, pcm, sigma, nspin, no_e)
    end if
 
 
    safe_deallocate_a(dipole)
    if (present(ref_file)) then
      safe_deallocate_a(ref_dipole)
    end if
 
    safe_allocate(sf(1:no_e, nspin))
 
    if (abs(kick%delta_strength) < 1e-12_real64) kick%delta_strength = m_one
    do ie = 1, no_e
      energy = (ie-1) * spectrum%energy_step + spectrum%min_energy
      do isp = 1, nspin
        sf(ie, isp) = sum(sigma(ie, 1:3, isp)*kick%pol(1:3, kick%pol_dir))
      end do
      sf(ie, 1:nspin) = -sf(ie, 1:nspin) * (energy * m_two) / (m_pi * kick%delta_strength)
      sigma(ie, 1:3, 1:nspin) = -sigma(ie, 1:3, 1:nspin)*(m_four*m_pi*energy/p_c)/kick%delta_strength
    end do
 
    ! The formulae below are only correct in this particular case.
    if (kick_get_type(kick) == kick_density_mode .and. spectrum%transform == spectrum_transform_sin) then
      ewsum = sum(sf(1, 1:nspin))
      polsum = m_zero
 
      do ie = 2, no_e
        energy = (ie-1) * spectrum%energy_step + spectrum%min_energy
        ewsum = ewsum + sum(sf(ie, 1:nspin))
        polsum = polsum + sum(sf(ie, 1:nspin)) / energy**2
      end do
 
      ewsum = ewsum * spectrum%energy_step
      polsum = polsum * spectrum%energy_step
    end if
 
    write(out_file, '(a15,i2)')      '# nspin        ', nspin
    call kick_write(kick, out_file)
    write(out_file, '(a)') '#%'
    write(out_file, '(a,i8)')    '# Number of time steps = ', time_steps
    call spectrum_write_info(spectrum, out_file)
    write(out_file, '(a)') '#%'
    if (kick_get_type(kick) == kick_density_mode .and. spectrum%transform == spectrum_transform_sin) then
      write(out_file, '(a,f16.6)') '# Electronic sum rule       = ', ewsum
      write(out_file, '(a,f16.6,1x,a)') '# Static polarizability (from sum rule) = ', &
        units_from_atomic(units_out%length**3, polsum), trim(units_abbrev(units_out%length))
      write(out_file, '(a)') '#%'
    end if
 
    write(out_file, '(a1,a20)', advance = 'no') '#', str_center("Energy", 20)
    do isp = 1, nspin
      do idir = 1, 3
        write(header_string,'(a6,i1,a8,i1,a1)') 'sigma(', idir, ', nspin=', isp, ')'
        write(out_file, '(a20)', advance = 'no') str_center(trim(header_string), 20)
      end do
    end do
    do isp = 1, nspin
      write(header_string,'(a18,i1,a1)') 'StrengthFunction(', isp, ')'
      write(out_file, '(a20)', advance = 'no') str_center(trim(header_string), 20)
    end do
    write(out_file, '(1x)')
    write(out_file, '(a1,a20)', advance = 'no') '#', str_center('['//trim(units_abbrev(units_out%energy)) // ']', 20)
    do ii = 1, nspin * 3
      write(out_file, '(a20)', advance = 'no') str_center('[' // trim(units_abbrev(units_out%length**2)) // ']', 20)
    end do
    do isp = 1, nspin
      write(out_file, '(a20)', advance = 'no') str_center('[' // trim(units_abbrev(unit_one/units_out%energy)) // ']', 20)
    end do
    write(out_file, '(1x)')
 
    do ie = 1, no_e
      write(out_file,'(e20.8)', advance = 'no') units_from_atomic(units_out%energy, &
        (ie-1) * spectrum%energy_step + spectrum%min_energy)
      do isp = 1, nspin
        write(out_file,'(3e20.8)', advance = 'no') (units_from_atomic(units_out%length**2, sigma(ie, idir, isp)), &
          idir = 1, 3)
      end do
      do isp = 1, nspin
        write(out_file,'(e20.8)', advance = 'no') units_from_atomic(unit_one/units_out%energy, sf(ie, isp))
      end do
      write(out_file, '(1x)')
    end do
 
    safe_deallocate_a(sigma)
    pop_sub(spectrum_cross_section)
  end subroutine spectrum_cross_section
 
  ! ---------------------------------------------------------
 
  subroutine spectrum_read_dipole(namespace, in_file, dipole)
    type(namespace_t), intent(in)    :: namespace
    integer,           intent(in)    :: in_file
    real(real64),      intent(out)   :: dipole(0:, :, :)
 
    integer :: nspin, lmax, time_steps, trash, it, idir, ispin
    real(real64)   :: dt,  dump
    type(kick_t) :: kick
    type(unit_system_t) :: file_units
 
    push_sub(spectrum_read_dipole)
 
    ! This function gives us back the unit connected to the "multipoles" file, the header information,
    ! the number of time steps, and the time step.
    call spectrum_mult_info(namespace, in_file, nspin, kick, time_steps, dt, file_units, lmax = lmax)
 
    ! Read the dipole.
    call io_skip_header(in_file)
 
    do it = 0, time_steps
      dipole(it, :, :) = m_zero
      read(in_file, *) trash, dump, (dump, (dipole(it, idir, ispin), idir = 1, kick%dim), ispin = 1, nspin)
    end do
    dipole(:,:,:) = units_to_atomic(file_units%length, dipole(:,:,:))
 
    pop_sub(spectrum_read_dipole)
 
  end subroutine spectrum_read_dipole
 
  ! ---------------------------------------------------------
 
  subroutine spectrum_add_pcm_dipole(namespace, dipole, time_steps, nspin)
    type(namespace_t), intent(in)    :: namespace
    real(real64),      intent(inout) :: dipole(0:, :, :)
    integer,           intent(in)    :: time_steps
    integer,           intent(in)    :: nspin
 
    type(pcm_t) :: pcm
    real(real64) :: dipole_pcm(1:3)
    integer :: ia, it
 
    ! unit io variables
    integer :: asc_unit_test
    integer :: cavity_unit
    integer :: asc_vs_t_unit, asc_vs_t_unit_check
    integer :: dipole_vs_t_unit_check, dipole_vs_t_unit_check1
    integer :: iocheck
    integer :: aux_int
    real(real64) :: aux_float, aux_float1, aux_vec(1:3)
    character(len=23) :: asc_vs_t_unit_format
    character(len=16) :: asc_vs_t_unit_format_tail
 
    push_sub(spectrum_add_pcm_dipole)
 
    ! reading PCM cavity from standard output file in two steps
 
    ! first step - counting tesserae
    asc_unit_test = io_open(pcm_dir//'ASC_e.dat', namespace, action='read')
    pcm%n_tesserae = 0
    iocheck = 1
    do while(iocheck >= 0)
      read(asc_unit_test,*,iostat=iocheck) aux_vec(1:3), aux_float, aux_int
      if (iocheck >= 0) pcm%n_tesserae = pcm%n_tesserae + 1
    end do
    call io_close(asc_unit_test)
 
    ! intermezzo - allocating PCM tessellation and polarization charges arrays
    safe_allocate(pcm%tess(1:pcm%n_tesserae))
    safe_allocate(pcm%q_e(1:pcm%n_tesserae))
    safe_allocate(pcm%q_e_in(1:pcm%n_tesserae)) ! with auxiliary role
 
    ! second step - reading of PCM tessellation arrays from standard output file
    !               writing the cavity to debug-purpose file
    asc_unit_test = io_open(pcm_dir//'ASC_e.dat', namespace, action='read')
    cavity_unit = io_open(pcm_dir//'cavity_check.xyz', namespace, action='write')
    write(cavity_unit,'(I3)') pcm%n_tesserae
    write(cavity_unit,*)
    do ia = 1, pcm%n_tesserae
      read(asc_unit_test,*) pcm%tess(ia)%point(1:3), aux_float, aux_int
      write(cavity_unit,'(A1,3(1X,F14.8))') 'H', pcm%tess(ia)%point(1:3)
    end do
    call io_close(asc_unit_test)
    call io_close(cavity_unit)
 
    write (asc_vs_t_unit_format_tail,'(I5,A11)') pcm%n_tesserae,'(1X,F14.8))'
    write (asc_vs_t_unit_format,'(A)') '(F14.8,'//trim(adjustl(asc_vs_t_unit_format_tail))
 
    ! Now, summary: * read the time-dependent PCM polarization charges due to solute electrons, pcm%q_e
    !               * compute the real-time dipole generated by pcm%q_e, dipole_pcm
    !               * add it to the real-time molecular dipole
    !               * write the total dipole and its PCM contribution to debug-purpose files
    ! N.B. we assume nuclei fixed in time
 
    ! opening time-dependent PCM charges standard and debug-purpose file
    asc_vs_t_unit = io_open(pcm_dir//'ASC_e_vs_t.dat', namespace, action='read', form='formatted')
    asc_vs_t_unit_check = io_open(pcm_dir//'ASC_e_vs_t_check.dat', namespace, action='write', form='formatted')
 
    ! opening time-dependent PCM and total dipole debug-purpose files
    dipole_vs_t_unit_check = io_open(pcm_dir//'dipole_e_vs_t_check.dat', namespace, action='write', form='formatted')
    dipole_vs_t_unit_check1 = io_open(pcm_dir//'dipole_e_vs_t_check1.dat', namespace, action='write', form='formatted')
 
    ! reading PCM charges for the zeroth-iteration - not used - pcm%q_e_in is only auxiliary here
    read(asc_vs_t_unit,trim(adjustl(asc_vs_t_unit_format))) aux_float1, ( pcm%q_e_in(ia) , ia=1,pcm%n_tesserae)
 
    do it = 1, time_steps
 
      ! reading real-time PCM charges due to electrons per timestep
      read(asc_vs_t_unit,trim(adjustl(asc_vs_t_unit_format))) aux_float, ( pcm%q_e(ia) , ia=1,pcm%n_tesserae)
 
      ! computing real-time PCM dipole per timestep
      call pcm_dipole(dipole_pcm(1:3), -pcm%q_e(1:pcm%n_tesserae), pcm%tess, pcm%n_tesserae)
 
      ! adding PCM dipole to the molecular dipole per timestep
      dipole(it, 1, 1:nspin) = dipole(it, 1, 1:nspin) + dipole_pcm(1)
      dipole(it, 2, 1:nspin) = dipole(it, 2, 1:nspin) + dipole_pcm(2)
      dipole(it, 3, 1:nspin) = dipole(it, 3, 1:nspin) + dipole_pcm(3)
 
      ! since we always have a kick for the optical spectrum in Octopus
      ! the first-iteration dipole should be equal to the zeroth-iteration one
      ! in any case, made them equal by hand
      if (it == 1) then
        dipole(0, 1, 1:nspin) = dipole(1, 1, 1:nspin)
        dipole(0, 2, 1:nspin) = dipole(1, 2, 1:nspin)
        dipole(0, 3, 1:nspin) = dipole(1, 3, 1:nspin)
      end if
 
      ! writing real-time PCM charges and dipole, and the total dipole for debug purposes
      write(asc_vs_t_unit_check,trim(adjustl(asc_vs_t_unit_format))) aux_float, (pcm%q_e(ia), ia=1,pcm%n_tesserae)
      write(dipole_vs_t_unit_check,'(F14.8,3(1X,F14.8))') aux_float, dipole_pcm
      write(dipole_vs_t_unit_check1,'(F14.8,3(1X,F14.8))') aux_float, dipole(it,:,1)
 
    end do
 
    ! closing PCM and debug files
    call io_close(asc_vs_t_unit)
    call io_close(asc_vs_t_unit_check)
    call io_close(dipole_vs_t_unit_check)
    call io_close(dipole_vs_t_unit_check1)
 
    ! deallocating PCM arrays
    safe_deallocate_a(pcm%tess)
    safe_deallocate_a(pcm%q_e)
    safe_deallocate_a(pcm%q_e_in)
 
    pop_sub(spectrum_add_pcm_dipole)
 
  end subroutine spectrum_add_pcm_dipole
 
  ! ---------------------------------------------------------
 
  subroutine spectrum_times_pcm_epsilon(spectrum, pcm, dipole, sigma, nspin, istart, iend, kick_time, dt, no_e)
    type(spectrum_t),   intent(in)    :: spectrum
    type(pcm_min_t)   , intent(in)    :: pcm
    real(real64), allocatable, intent(inout) :: sigma(:, :, :)
    real(real64), allocatable, intent(in)    :: dipole(:, :, :)
    integer,            intent(in)    :: nspin
    real(real64),       intent(in)    :: kick_time
    integer,            intent(in)    :: istart, iend
    real(real64),       intent(in)    :: dt
    integer,            intent(in)    :: no_e
 
    real(real64), allocatable :: sigmap(:, :, :)
    type(batch_t) :: dipoleb, sigmab
 
    integer :: ie
 
    complex(real64), allocatable :: eps(:)
 
    push_sub(spectrum_times_pcm_epsilon)
 
    ! imaginary part of the polarizability
 
    call batch_init(dipoleb, 3, 1, nspin, dipole)
    call batch_init(sigmab, 3, 1, nspin, sigma)
 
    call spectrum_signal_damp(spectrum%damp, spectrum%damp_factor, istart + 1, iend + 1, kick_time, dt, dipoleb)
    call spectrum_fourier_transform(spectrum%method, spectrum_transform_sin, spectrum%noise, &
      istart + 1, iend + 1, kick_time, dt, dipoleb, spectrum%min_energy, spectrum%max_energy, spectrum%energy_step, sigmab)
 
    call dipoleb%end()
    call sigmab%end()
 
    ! real part of the polarizability
 
    safe_allocate(sigmap(1:no_e, 1:3, 1:nspin))
 
    call batch_init(dipoleb, 3, 1, nspin, dipole)
    call batch_init(sigmab, 3, 1, nspin, sigmap)
 
    call spectrum_signal_damp(spectrum%damp, spectrum%damp_factor, istart + 1, iend + 1, kick_time, dt, dipoleb)
    call spectrum_fourier_transform(spectrum%method, spectrum_transform_cos, spectrum%noise, &
      istart + 1, iend + 1, kick_time, dt, dipoleb, spectrum%min_energy, spectrum%max_energy, spectrum%energy_step, sigmab)
 
    call dipoleb%end()
    call sigmab%end()
 
    safe_allocate(eps(1:no_e))
 
    ! multiplying by the dielectric function and taking the imaginary part of the product
 
    do ie = 1, no_e
      call pcm_eps(pcm, eps(ie), (ie-1)*spectrum%energy_step + spectrum%min_energy)
      sigma(ie, 1:3, 1:nspin) = sigma(ie, 1:3, 1:nspin) * real(eps(ie), real64) + sigmap(ie, 1:3, 1:nspin) *aimag(eps(ie))
    end do
 
    safe_deallocate_a(sigmap)
    safe_deallocate_a(eps)
 
    pop_sub(spectrum_times_pcm_epsilon)
 
  end subroutine spectrum_times_pcm_epsilon
 
  ! ---------------------------------------------------------
 
  subroutine spectrum_over_pcm_refraction_index(spectrum, pcm, sigma, nspin, no_e)
    type(spectrum_t),   intent(in)    :: spectrum
    type(pcm_min_t)   , intent(in)    :: pcm
    real(real64), allocatable, intent(inout) :: sigma(:, :, :)
    integer,            intent(in)    :: nspin
    integer,            intent(in)    :: no_e
 
    integer :: ie
 
    complex(real64), allocatable :: eps(:)
 
    push_sub(spectrum_over_pcm_refraction_index)
 
    safe_allocate(eps(1:no_e))
 
    ! dividing by the refraction index - n(\omega)=\sqrt{\frac{|\epsilon(\omega)|+\Re[\epsilon(\omega)]}{2}}
 
    do ie = 1, no_e
      call pcm_eps(pcm, eps(ie), (ie-1)*spectrum%energy_step + spectrum%min_energy)
      sigma(ie, 1:3, 1:nspin) = sigma(ie, 1:3, 1:nspin) / sqrt(0.5_real64 * (abs(eps(ie)) + real(eps(ie), real64)))
    end do
 
    safe_deallocate_a(eps)
 
    pop_sub(spectrum_over_pcm_refraction_index)
 
  end subroutine spectrum_over_pcm_refraction_index
 
  ! ---------------------------------------------------------
  subroutine spectrum_dipole_power(spectrum, namespace, in_file, out_file)
    type(spectrum_t),  intent(inout) :: spectrum
    type(namespace_t), intent(in)    :: namespace
    integer,           intent(in)    :: in_file
    integer,           intent(in)    :: out_file
 
    character(len=20) :: header_string
    integer :: nspin, lmax, time_steps, istart, iend, ntiter, it, ii, isp, no_e, ie, idir
    real(real64)   :: dt
    real(real64), allocatable :: dipole(:, :, :), transform_cos(:, :, :), transform_sin(:, :, :), power(:, :, :)
    type(unit_system_t) :: file_units
    type(batch_t) :: dipoleb, transformb_cos, transformb_sin
    type(kick_t) :: kick
 
    push_sub(spectrum_dipole_power)
 
    ! This function gives us back the unit connected to the "multipoles" file, the header information,
    ! the number of time steps, and the time step.
    call spectrum_mult_info(namespace, in_file, nspin, kick, time_steps, dt, file_units, lmax=lmax)
 
    ! Now we cannot process files that do not contain the dipole, or that contain more than the dipole.
    if (lmax /= 1) then
      message(1) = 'Multipoles file should contain the dipole -- and only the dipole.'
      call messages_fatal(1, namespace=namespace)
    end if
 
    ! Find out the iteration numbers corresponding to the time limits.
    call spectrum_fix_time_limits(spectrum, time_steps, dt, istart, iend, ntiter)
 
    safe_allocate(dipole(0:time_steps, 1:3, 1:nspin))
    call spectrum_read_dipole(namespace, in_file, dipole)
 
    ! Now subtract the initial dipole.
    do it = 1, time_steps
      dipole(it, :, :) = dipole(it, :, :) - dipole(0, :, :)
    end do
    dipole(0, :, :) = m_zero
 
    if (spectrum%energy_step <= m_zero) spectrum%energy_step = m_two * m_pi / (dt*time_steps)
 
    ! Get the number of energy steps.
    no_e = spectrum_nenergy_steps(spectrum)
    safe_allocate(transform_cos(1:no_e, 1:3, 1:nspin))
    safe_allocate(transform_sin(1:no_e, 1:3, 1:nspin))
    safe_allocate(power(1:no_e, 1:3, 1:nspin))
 
 
    call batch_init(dipoleb, 3, 1, nspin, dipole)
    call batch_init(transformb_cos, 3, 1, nspin, transform_cos)
    call batch_init(transformb_sin, 3, 1, nspin, transform_sin)
 
    call spectrum_signal_damp(spectrum%damp, spectrum%damp_factor, istart + 1, iend + 1, spectrum%start_time, dt, dipoleb)
 
    call spectrum_fourier_transform(spectrum%method, spectrum_transform_cos, spectrum%noise, &
      istart + 1, iend + 1, spectrum%start_time, dt, dipoleb, spectrum%min_energy,  &
      spectrum%max_energy, spectrum%energy_step, transformb_cos)
    call spectrum_fourier_transform(spectrum%method, spectrum_transform_sin, spectrum%noise, &
      istart + 1, iend + 1, spectrum%start_time, dt, dipoleb, spectrum%min_energy, &
      spectrum%max_energy, spectrum%energy_step, transformb_sin)
 
    do ie = 1, no_e
      power(ie, :, :) = (transform_sin(ie, :, :)**2 + transform_cos(ie, :, :)**2)
    end do
 
    call dipoleb%end()
    call transformb_cos%end()
    call transformb_sin%end()
 
    safe_deallocate_a(dipole)
    safe_deallocate_a(transform_sin)
    safe_deallocate_a(transform_cos)
 
    write(out_file, '(a15,i2)')      '# nspin        ', nspin
    write(out_file, '(a)') '#%'
    write(out_file, '(a,i8)')    '# Number of time steps = ', time_steps
    call spectrum_write_info(spectrum, out_file)
    write(out_file, '(a)') '#%'
 
    write(out_file, '(a1,a20,1x)', advance = 'no') '#', str_center("Energy", 20)
    do isp = 1, nspin
      do idir = 1, 3
        write(header_string,'(a6,i1,a8,i1,a1)') 'power(', idir, ', nspin=', isp, ')'
        write(out_file, '(a20)', advance = 'no') str_center(trim(header_string), 20)
      end do
    end do
    write(out_file, '(1x)')
    write(out_file, '(a1,a20)', advance = 'no') '#', str_center('['//trim(units_abbrev(units_out%energy)) // ']', 20)
    do ii = 1, nspin * 3
      write(out_file, '(a20)', advance = 'no') str_center('[' // trim(units_abbrev(units_out%length**2)) // ']', 20)
    end do
    write(out_file, '(1x)')
 
    do ie = 1, no_e
      write(out_file,'(e20.8)', advance = 'no') units_from_atomic(units_out%energy, &
        (ie-1) * spectrum%energy_step + spectrum%min_energy)
      do isp = 1, nspin
        write(out_file,'(3e20.8)', advance = 'no') (units_from_atomic(units_out%length**2, power(ie, idir, isp)), &
          idir = 1, 3)
      end do
      write(out_file, '(1x)')
    end do
 
    safe_deallocate_a(power)
 
    pop_sub(spectrum_dipole_power)
  end subroutine spectrum_dipole_power
 
  ! ---------------------------------------------------------
  subroutine spectrum_dyn_structure_factor(spectrum, namespace, in_file_sin, in_file_cos, out_file)
    type(spectrum_t),  intent(inout) :: spectrum
    type(namespace_t), intent(in)    :: namespace
    integer,           intent(in)    :: in_file_sin, in_file_cos
    integer,           intent(in)    :: out_file
 
    character(len=20) :: header_string
    integer :: time_steps, time_steps_sin, time_steps_cos
    integer :: istart, iend, ntiter, it, jj, ii, no_e, ie, trash
    real(real64)   :: dt, dt_sin, dt_cos
    real(real64)   :: dump, dummy1, dummy2, dummy3, dummy4, energy, fsum
    type(kick_t) :: kick
    complex(real64) :: xx
    complex(real64), allocatable :: ftchd(:), chi(:), damp(:)
    type(unit_system_t) :: file_units
    character(len=100) :: line
 
    push_sub(spectrum_dyn_structure_factor)
 
    safe_allocate(kick%qvector(1:3, 1:1))
 
    ! Read information from ftchds.sin file
 
    rewind(in_file_sin)
 
    ! skip two lines
    read(in_file_sin, *)
    read(in_file_sin, *)
    read(in_file_sin, '(15x,i2)')     kick%qkick_mode
    read(in_file_sin, '(10x,3f9.5)')  kick%qvector
    read(in_file_sin, '(15x,f18.12)') kick%delta_strength
 
    ! skip two lines
    read(in_file_sin, *)
    read(in_file_sin, '(a)') line
    call io_skip_header(in_file_sin)
    call io_skip_header(in_file_cos)
 
    ! Figure out the units of the file
    ii = index(line, 'eV')
    if (ii /= 0) then
      call unit_system_get(file_units, units_eva)
    else
      call unit_system_get(file_units, units_atomic)
    end if
 
    ! get time_steps and dt, and make sure that dt is the same in the two files
    call spectrum_count_time_steps(namespace, in_file_sin, time_steps_sin, dt_sin)
    call spectrum_count_time_steps(namespace, in_file_cos, time_steps_cos, dt_cos)
 
    if (.not. is_close(dt_sin, dt_cos)) then
      message(1) = "dt is different in ftchds.cos and ftchds.sin!"
      call messages_fatal(1, namespace=namespace)
    end if
 
    time_steps = min(time_steps_sin, time_steps_cos)
    dt = units_to_atomic(file_units%time, dt_cos) ! units_out is OK
 
    ! Find out the iteration numbers corresponding to the time limits.
    call spectrum_fix_time_limits(spectrum, time_steps, dt, istart, iend, ntiter)
 
    ! Read the f-transformed charge density.
    call io_skip_header(in_file_sin)
    call io_skip_header(in_file_cos)
 
    safe_allocate(ftchd(0:time_steps))
    do it = 0, time_steps
      read(in_file_sin, *) trash, dump, dummy1, dummy2
      read(in_file_cos, *) trash, dump, dummy3, dummy4
      ftchd(it) = cmplx(dummy3-dummy2, dummy4+dummy1, real64)
    end do
 
    ! Now subtract the initial value.
    do it = 1, time_steps
      ftchd(it) = ftchd(it) - ftchd(0)
    end do
    ftchd(0) = m_zero
 
    if (spectrum%energy_step <= m_zero) spectrum%energy_step = m_two * m_pi / (dt*time_steps)
 
    ! Get the number of energy steps.
    no_e = spectrum_nenergy_steps(spectrum)
 
    safe_allocate(chi(1:no_e))
    chi = m_zero
 
    ! Gets the damp function
    safe_allocate(damp(istart:iend))
    do it = istart, iend
      jj = it - istart
      select case (spectrum%damp)
      case (spectrum_damp_none)
        damp(it) = m_one
      case (spectrum_damp_lorentzian)
        damp(it)= exp(-jj * dt * spectrum%damp_factor)
      case (spectrum_damp_polynomial)
        damp(it) = m_one - m_three * (real(jj, real64)  / ntiter)**2 &
          + m_two * (real(jj, real64)  / ntiter)**3
      case (spectrum_damp_gaussian)
        damp(it)= exp(-(jj * dt)**2 * spectrum%damp_factor**2)
      end select
    end do
 
    ! Fourier transformation from time to frequency
    if (abs(kick%delta_strength) < 1.d-12) kick%delta_strength = m_one
    do ie = 1, no_e
      energy = (ie-1) * spectrum%energy_step + spectrum%min_energy
      do it = istart, iend
        jj = it - istart
 
        xx = exp(m_zi * energy * jj * dt)
        chi(ie) = chi(ie) + xx * damp(it) * ftchd(it)
 
      end do
      chi(ie) = chi(ie) * dt / kick%delta_strength / m_pi
    end do
 
    ! Test f-sum rule
    fsum = m_zero
    do ie = 1, no_e
      energy = (ie-1) * spectrum%energy_step + spectrum%min_energy
      fsum = fsum + energy * aimag(chi(ie))
    end do
    fsum = spectrum%energy_step * fsum * 2/sum(kick%qvector(:,1)**2)
 
    write(out_file, '(a)') '#%'
    write(out_file, '(a,i8)')    '# Number of time steps = ', time_steps
    call spectrum_write_info(spectrum, out_file)
    write(out_file, '(a,3f9.5)') '# qvector    : ', kick%qvector
    write(out_file, '(a,f10.4)') '# F-sum rule : ', fsum
    write(out_file, '(a)') '#%'
 
    write(out_file, '(a1,a20)', advance = 'no') '#', str_center("Energy", 20)
    write(header_string,'(a3)') 'chi'
    write(out_file, '(a20)', advance = 'no') str_center(trim(header_string), 20)
    write(out_file, '(1x)')
    write(out_file, '(a1,a20)', advance = 'no') '#', str_center('['//trim(units_abbrev(units_out%energy)) // ']', 20)
    write(out_file, '(a20)', advance = 'no') str_center('[' // trim(units_abbrev(units_out%energy)) // '**(-1)]', 20)
    write(out_file, '(1x)')
 
    do ie = 1, no_e
      write(out_file,'(e20.8)', advance = 'no') units_from_atomic(units_out%energy, &
        (ie-1) * spectrum%energy_step + spectrum%min_energy)
      write(out_file,'(e20.8)', advance = 'no') units_from_atomic(units_out%energy**(-1), aimag(chi(ie)))
      write(out_file, '(1x)')
    end do
 
    safe_deallocate_a(ftchd)
    safe_deallocate_a(chi)
    pop_sub(spectrum_dyn_structure_factor)
 
  end subroutine spectrum_dyn_structure_factor
 
 
  ! ---------------------------------------------------------
  subroutine spectrum_rotatory_strength(spectrum, namespace, in_file, out_file)
    type(spectrum_t),  intent(inout) :: spectrum
    type(namespace_t), intent(in)    :: namespace
    integer,           intent(in)    :: in_file
    integer,           intent(in)    :: out_file
 
    integer :: istart, iend, ntiter, ie, idir, time_steps, no_e, nspin, trash, it
    real(real64) :: dump, dt, energy
    type(kick_t) :: kick
    complex(real64) :: sum1, sum2, sp
    real(real64), allocatable :: angular(:, :), resp(:), imsp(:)
    type(batch_t) :: angularb, respb, imspb
    type(unit_system_t) :: file_units
 
    push_sub(spectrum_rotatory_strength)
 
    call spectrum_mult_info(namespace, in_file, nspin, kick, time_steps, dt, file_units)
    call spectrum_fix_time_limits(spectrum, time_steps, dt, istart, iend, ntiter)
 
    if (kick%dim /= 3) then
      message(1) = "Rotatory strength can only be computed for 3D systems."
      call messages_fatal(1, namespace=namespace)
    end if
 
    ! load angular momentum from file
    safe_allocate(angular(0:time_steps, 1:3))
    call io_skip_header(in_file)
    do ie = 0, time_steps
      read(in_file, *) trash, dump, (angular(ie, idir), idir = 1, 3)
    end do
 
    ! subtract static dipole
    do idir = 1, 3
      angular(:, idir) = angular(:, idir) - angular(0, idir)
    end do
 
    if (spectrum%energy_step <= m_zero) spectrum%energy_step = m_two * m_pi / (dt*time_steps)
 
    no_e = spectrum_nenergy_steps(spectrum)
 
    do it = istart, iend
      angular(it, 1) = sum(angular(it, 1:3)*kick%pol(1:3, kick%pol_dir))
    end do
 
    safe_allocate(resp(1:no_e))
    safe_allocate(imsp(1:no_e))
 
    call batch_init(angularb, angular(:, 1))
    call batch_init(respb, resp)
    call batch_init(imspb, imsp)
 
    call spectrum_signal_damp(spectrum%damp, spectrum%damp_factor, istart + 1, iend + 1, kick%time, dt, angularb)
 
    call spectrum_fourier_transform(spectrum%method, spectrum_transform_cos, spectrum%noise, &
      istart + 1, iend + 1, kick%time, dt, angularb, spectrum%min_energy, spectrum%max_energy, spectrum%energy_step, respb)
    call spectrum_fourier_transform(spectrum%method, spectrum_transform_sin, spectrum%noise, &
      istart + 1, iend + 1, kick%time, dt, angularb, spectrum%min_energy, spectrum%max_energy, spectrum%energy_step, imspb)
 
    call angularb%end()
    call respb%end()
    call imspb%end()
 
    sum1 = m_z0
    sum2 = m_z0
    if (abs(kick%delta_strength) < 1.d-12) kick%delta_strength = m_one
    do ie = 1, no_e
      energy = (ie-1) * spectrum%energy_step + spectrum%min_energy
 
      sp = cmplx(resp(ie), imsp(ie), real64)
 
      sp = sp*m_zi/(m_two*p_c*kick%delta_strength)
 
      sum1 = sum1 + spectrum%energy_step*sp
      sum2 = sum2 + spectrum%energy_step*sp*energy**2
 
      resp(ie) = real(sp, real64)
      imsp(ie) = aimag(sp)
    end do
 
    safe_deallocate_a(angular)
 
    ! print some info
    write(message(1), '(a,i8)')    'Number of time steps = ', ntiter
    write(message(2), '(a,i4)')    'PropagationSpectrumDampMode   = ', spectrum%damp
    write(message(3), '(a,f10.4)') 'PropagationSpectrumDampFactor = ', units_from_atomic(units_out%time**(-1), spectrum%damp_factor)
    write(message(4), '(a,f10.4)') 'PropagationSpectrumStartTime  = ', units_from_atomic(units_out%time, spectrum%start_time)
    write(message(5), '(a,f10.4)') 'PropagationSpectrumEndTime    = ', units_from_atomic(units_out%time, spectrum%end_time)
    write(message(6), '(a,f10.4)') 'PropagationSpectrumMaxEnergy  = ', units_from_atomic(units_inp%energy, spectrum%max_energy)
    write(message(7),'(a,f10.4)')  'PropagationSpectrumEnergyStep = ', units_from_atomic(units_inp%energy, spectrum%energy_step)
    message(8) = ""
    write(message(9), '(a,5e15.6,5e15.6)') 'R(0) sum rule = ', sum1
    write(message(10),'(a,5e15.6,5e15.6)') 'R(2) sum rule = ', sum2
    call messages_info(10, namespace=namespace)
 
 
    ! Output to file
    write(out_file, '(a15,i2)')      '# nspin        ', nspin
    call kick_write(kick, out_file)
    write(out_file, '(a1,a20,a20,a20)') '#', str_center("Energy", 20), str_center("R", 20), str_center("Re[beta]", 20)
    write(out_file, '(a1,a20,a20,a20)') '#', str_center('[' // trim(units_abbrev(units_out%energy)) // ']', 20), &
      str_center('[' // trim(units_abbrev(units_out%length**3)) // ']', 20), &
      str_center('[' // trim(units_abbrev(units_out%length**4)) // ']', 20)
    write(out_file, '(a,5e15.6,5e15.6)') '# R(0) sum rule = ', sum1
    write(out_file, '(a,5e15.6,5e15.6)') '# R(2) sum rule = ', sum2
    do ie = 1, no_e
      write(out_file,'(e20.8,e20.8,e20.8)') units_from_atomic(units_out%energy, (ie-1)*spectrum%energy_step+spectrum%min_energy), &
        units_from_atomic(units_out%length**3, imsp(ie)/m_pi), &
        units_from_atomic(units_out%length**4, resp(ie)*p_c/(m_three*max((ie-1),1)*spectrum%energy_step))
    end do
 
    safe_deallocate_a(resp)
    safe_deallocate_a(imsp)
 
    pop_sub(spectrum_rotatory_strength)
  end subroutine spectrum_rotatory_strength
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine spectrum_hsfunction_init(dt, is, ie, niter, acc)
    real(real64),      intent(in)    :: dt
    integer,           intent(in)    :: is, ie, niter
    complex(real64),   intent(in)    :: acc(:)
 
    integer :: nn(3), j, optimize_parity(3)
    logical :: optimize(3)
 
    push_sub(spectrum_hsfunction_init)
 
    is_ = is
    ie_ = ie
    time_step_ = dt
    niter_ = niter
    energy_step_ = (m_two * m_pi) / (niter * time_step_)
    safe_allocate(func_(0:niter))
    safe_allocate(funcw_(0:niter))
    func_ = m_z0
    func_ = acc
    nn(1:3) = (/ niter, 1, 1 /)
    optimize(1:3) = .false.
    optimize_parity(1:3) = -1
 
    call fft_init(fft_handler, nn(1:3), 1, fft_complex, fftlib_fftw, optimize, optimize_parity)
    call zfft_forward(fft_handler, func_(0:niter-1), funcw_(0:niter-1))
    do j = 0, niter - 1
      funcw_(j) = -abs(funcw_(j))**2 * dt**2
    end do
 
    pop_sub(spectrum_hsfunction_init)
  end subroutine spectrum_hsfunction_init
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine spectrum_hsfunction_end
 
    push_sub(spectrum_hsfunction_end)
 
    call fft_end(fft_handler)
 
    safe_deallocate_a(func_)
    safe_deallocate_a(funcw_)
    pop_sub(spectrum_hsfunction_end)
  end subroutine spectrum_hsfunction_end
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine spectrum_hsfunction_min(namespace, aa, bb, omega_min, func_min)
    type(namespace_t), intent(in)  :: namespace
    real(real64),      intent(in)  :: aa, bb
    real(real64),      intent(out) :: omega_min, func_min
 
    integer :: ierr, ie
    real(real64) :: xx, hsval, minhsval, ww, xa, xb, hxa, hxb
 
    push_sub(spectrum_hsfunction_min)
 
    ! xx should be an initial guess for the minimum. So we do a quick search
    ! that we refine later calling 1dminimize.
    !xx = omega
    !call hsfunction(xx, minhsval)
 
    ierr = 0
 
    ie = int(aa/energy_step_)
    ww = ie * energy_step_
    if (ww < aa) then
      ie = ie + 1
      ww = ie * energy_step_
    end if
    xx = ie * energy_step_
    minhsval = real(funcw_(ie), real64)
    do while(ww <= bb)
      hsval = real(funcw_(ie), real64)
      if (hsval < minhsval) then
        minhsval = hsval
        xx = ww
      end if
      ie = ie + 1
      ww = ie * energy_step_
    end do
 
    ! Around xx, we call some GSL sophisticated search algorithm to find the minimum.
    ! First, we get the value of the function at the extremes of the interval
    xa = max(xx-energy_step_, aa)
    xb = min(xx+energy_step_, bb)
    call hsfunction(xa, hxa)
    call hsfunction(xb, hxb)
 
    if (hxa <= minhsval) then
      xx = xa
      minhsval = hxa
    elseif (hxb <= minhsval) then
      xx = xb
      minhsval = hxb
    else
      call loct_1dminimize(xa, xb, xx, hsfunction, ierr)
    end if
 
    if (ierr /= 0) then
      write(message(1),'(a,f14.6,a)') 'spectrum_hsfunction_min: The maximum at', xx,' was not properly converged.'
      write(message(2),'(a,i12)')      'Error code: ', ierr
      call messages_warning(2, namespace=namespace)
    end if
    call hsfunction(xx, hsval)
    omega_min = xx
    func_min  = hsval
 
    pop_sub(spectrum_hsfunction_min)
  end subroutine spectrum_hsfunction_min
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine hsfunction(omega, power)
    real(real64), intent(in)  :: omega
    real(real64), intent(out) :: power
 
    complex(real64)   :: cc, ez1, ez, zz
    integer :: jj,dir
 
    push_sub(hsfunction)
 
    zz = m_zi * omega * time_step_
    ez  = exp(zz)
 
    if (allocated(func_ar_)) then
      power = m_zero
      do dir = 1, 3
        ez1 = exp((is_ - 1) * zz)
        cc = m_z0
        do jj = is_, ie_
          ! This would be easier, but slower.
          !cc = cc + exp(M_zI * omega * jj * time_step_) * func_(jj)
          ez1 = ez1 * ez
          cc = cc + ez1 * func_ar_(dir,jj) &
            *exp(-m_zi * omega * tret_(jj)) !integrate over the retarded time
        end do
        power = power - abs(cc)**2 * time_step_**2
!        pp(dir) = cc * time_step_
      end do
!      power = - sum(abs(zcross_product(vv_, pp))**2)
    else
      cc = m_z0
      ez1 = exp((is_ - 1) * zz)
      do jj = is_, ie_
        ez1 = ez1 * ez
        cc = cc + ez1 * func_(jj)
      end do
      power = -abs(cc)**2 * time_step_**2
    end if
 
 
    pop_sub(hsfunction)
  end subroutine hsfunction
  ! ---------------------------------------------------------
 
  ! ---------------------------------------------------------
  subroutine spectrum_hsfunction_ar_init(dt, is, ie, niter, acc, pos,tret)
    real(real64),   intent(in) :: dt
    integer, intent(in) :: is, ie, niter
    complex(real64),   intent(in) :: acc(:,:),pos(:,:),tret(:)
 
    push_sub(spectrum_hsfunction_ar_init)
 
    is_ = is
    ie_ = ie
    time_step_ = dt
    safe_allocate(func_ar_(1:3, 0:niter))
    safe_allocate(pos_(1:3, 0:niter))
    safe_allocate(tret_(0:niter))
    func_ar_ = acc
    pos_= pos
    tret_=tret
 
 
    pop_sub(spectrum_hsfunction_ar_init)
  end subroutine spectrum_hsfunction_ar_init
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  ! FIXME: why is this never called?
  subroutine spectrum_hsfunction_ar_end
 
    push_sub(spectrum_hsfunction_ar_end)
 
    safe_deallocate_a(func_ar_)
    safe_deallocate_a(pos_)
    safe_deallocate_a(tret_)
 
    pop_sub(spectrum_hsfunction_ar_end)
  end subroutine spectrum_hsfunction_ar_end
  ! ---------------------------------------------------------
 
  ! ---------------------------------------------------------
  subroutine spectrum_hs_ar_from_acc(spectrum, namespace, out_file, vec, w0)
    type(spectrum_t),  intent(inout) :: spectrum
    type(namespace_t), intent(in)    :: namespace
    character(len=*),  intent(in)    :: out_file
    real(real64),      intent(in)    :: vec(:)
    real(real64),  optional,  intent(in)    :: w0
 
    integer :: istep, trash, iunit, nspin, time_steps, istart, iend, ntiter, lmax, ierr, jj, idir, ispin
    real(real64) :: dt, dump, aa(3)
    complex(real64) :: nn(3)
    type(kick_t) :: kick
    real(real64), allocatable :: dd(:,:)
    complex(real64), allocatable :: acc(:,:),pp(:,:),pos(:,:),tret(:)
    real(real64) :: vv(3)
    type(unit_system_t) :: file_units
 
    push_sub(spectrum_hs_ar_from_acc)
 
    call spectrum_tdfile_info(namespace, 'acceleration', iunit, time_steps, dt)
    call spectrum_fix_time_limits(spectrum, time_steps, dt, istart, iend, ntiter)
 
    ! load dipole from file
    safe_allocate(acc(1:3, 0:time_steps))
    safe_allocate(pp(1:3, 0:time_steps))
    safe_allocate(pos(1:3, 0:time_steps))
    safe_allocate(tret(0:time_steps))
 
    acc = m_zero
    pos = m_zero
    pp = m_zero
    nn = m_zero
    tret = m_zero
 
    call io_skip_header(iunit)
    do istep = 0, time_steps-1
      aa = m_zero
      read(iunit, '(28x,e20.12)', advance = 'no', iostat = ierr) aa(1)
      jj = 2
      do while((ierr == 0) .and. (jj <= 3))
        read(iunit, '(e20.12)', advance = 'no', iostat = ierr) aa(jj)
        jj = jj+1
      end do
 
!      read(iunit, *) trash, dump, aa
 
      acc(:,istep) = units_to_atomic(units_out%acceleration, aa(:))
!      write (*,*) istep, Real(acc(:,istep))
    end do
    close(iunit)
 
 
    ! Try to get the trajectory from multipole file
 
    iunit = io_open('multipoles', namespace, action='read', status='old', die=.false.)
    if (iunit == -1) then
      iunit = io_open('td.general/multipoles', namespace, action='read', status='old')
    end if
    call spectrum_mult_info(namespace, iunit, nspin, kick, time_steps, dt, file_units, lmax=lmax)
    call spectrum_fix_time_limits(spectrum, time_steps, dt, istart, iend, ntiter)
 
    call io_skip_header(iunit)
!    write (*,*)
 
    safe_allocate(dd(1:3, 1:nspin))
    do istep = 0, time_steps-1
      dd = m_zero
      read(iunit, *) trash, dump, (dump, (dd(idir, ispin), idir = 1, kick%dim), ispin = 1, nspin)
      pos(1:3, istep) = -sum(dd(1:3, :),2)
      pos(:,istep) = units_to_atomic(units_out%length, pos(:,istep))
    end do
    safe_deallocate_a(dd)
    pos(:,0) = pos(:,1)
    call io_close(iunit)
 
!    write (*,*)
 
    ! normalize vector and set to global var
!    vv = vec / sqrt(sum(vec(:)**2))
    vv(1:3) = vec(1:3)
 
    pp(:,0) = m_zero
    do istep = 0, time_steps - 1
      nn(:) = vv(:)-pos(:,istep)
      nn(:) = nn(:)/norm2(abs(nn(:)))
      tret(istep) = dot_product(vv(:), real(pos(:,istep), real64))/p_c
      pp(:,istep) = zcross_product(nn, zcross_product(nn, acc(:,istep)))
!       write (*,*) istep, Real(PP(:,istep)),"acc", Real(acc(:,istep))
    end do
 
    call spectrum_hsfunction_ar_init(dt, istart, iend, time_steps, pp, pos, tret)
    call spectrum_hs(spectrum, namespace, out_file, 'a', w0)
    call spectrum_hsfunction_ar_end()
 
    safe_deallocate_a(acc)
    safe_deallocate_a(pp)
    safe_deallocate_a(pos)
    safe_deallocate_a(tret)
 
    pop_sub(spectrum_hs_ar_from_acc)
  end subroutine spectrum_hs_ar_from_acc
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine spectrum_hs_ar_from_mult(spectrum, namespace, out_file, vec, w0)
    type(spectrum_t),  intent(inout) :: spectrum
    type(namespace_t), intent(in)    :: namespace
    character(len=*),  intent(in)    :: out_file
    real(real64),      intent(in)    :: vec(:)
    real(real64),  optional,  intent(in)    :: w0
 
    integer :: istep, trash, iunit, nspin, time_steps, istart, iend, ntiter, lmax, idir, ispin
    real(real64) :: dt, dump
    type(kick_t) :: kick
    real(real64), allocatable :: dd(:,:)
    complex(real64), allocatable :: dipole(:,:), ddipole(:,:), pp(:,:), tret(:)
    complex(real64) :: vv(3)
    type(unit_system_t) :: file_units
 
    push_sub(spectrum_hs_ar_from_mult)
 
 
    iunit = io_open('multipoles', namespace, action='read', status='old', die=.false.)
    if (iunit == -1) then
      iunit = io_open('td.general/multipoles', namespace, action='read', status='old')
    end if
    call spectrum_mult_info(namespace, iunit, nspin, kick, time_steps, dt, file_units, lmax=lmax)
    call spectrum_fix_time_limits(spectrum, time_steps, dt, istart, iend, ntiter)
 
    call io_skip_header(iunit)
 
    ! load dipole from file
    safe_allocate(dipole(1:3, 0:time_steps))
    safe_allocate(ddipole(1:3, 0:time_steps))
    safe_allocate(pp(1:3, 0:time_steps))
    safe_allocate(tret(0:time_steps))
    safe_allocate(dd(1:3, 1:nspin))
 
    dipole = m_z0
    ddipole = m_z0
    pp = m_z0
    tret= m_zero
 
    do istep = 1, time_steps
      dd = m_zero
      read(iunit, *) trash, dump, (dump, (dd(idir, ispin), idir = 1, kick%dim), ispin = 1, nspin)
      dipole(1:3, istep) = -sum(dd(1:3, :),2)
      dipole(:,istep) = units_to_atomic(units_out%length, dipole(:,istep))
    end do
    safe_deallocate_a(dd)
    dipole(:,0) = dipole(:,1)
    call io_close(iunit)
 
    ! we now calculate the acceleration.
    ddipole(:,0) = m_zero
    do istep = 1, time_steps - 1
      ddipole(:,istep) = (dipole(:,istep - 1) + dipole(:,istep + 1) - m_two * dipole(:,istep)) / dt**2
    end do
    call interpolate(dt*(/ -3, -2, -1 /),   &
      ddipole(1,time_steps - 3:time_steps - 1), &
      m_zero, &
      ddipole(1,time_steps))
    call interpolate(dt*(/ -3, -2, -1 /),   &
      ddipole(2,time_steps - 3:time_steps - 1), &
      m_zero, &
      ddipole(2,time_steps))
    call interpolate(dt*(/ -3, -2, -1 /),   &
      ddipole(3,time_steps - 3:time_steps - 1), &
      m_zero, &
      ddipole(3,time_steps))
 
    ! normalize vector and set to global var
    vv(1:3) = vec(1:3) / norm2(vec(1:3))
 
    pp(:,0) = m_zero
    do istep = 1, time_steps - 1
!      write (*,*) istep, istep*dt, Real(ddipole(1,istep)), Real(ddipole(2,istep))
      tret(istep) = dot_product(vv(:), dipole(:,istep))/p_c
      pp(:,istep) = zcross_product(vv, zcross_product(vv, ddipole(:,istep - 1)))
!      PP(istep) = sum(abs(dipole(:,istep))**2)
!      write(*,*) istep, PP(istep)
 
    end do
 
 
    call spectrum_hsfunction_ar_init(dt, istart, iend, time_steps, pp, dipole,tret)
    call spectrum_hs(spectrum, namespace, out_file, 'a', w0)
    call spectrum_hsfunction_ar_end()
 
    safe_deallocate_a(dipole)
    safe_deallocate_a(ddipole)
    safe_deallocate_a(pp)
    safe_deallocate_a(tret)
 
 
    pop_sub(spectrum_hs_ar_from_mult)
  end subroutine spectrum_hs_ar_from_mult
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine spectrum_hs_from_mult(spectrum, namespace, out_file, pol, vec, w0)
    type(spectrum_t),  intent(inout) :: spectrum
    type(namespace_t), intent(in)    :: namespace
    character(len=*),  intent(in)    :: out_file
    character,         intent(in)    :: pol
    real(real64),      intent(in)    :: vec(:)
    real(real64),   optional, intent(in)    :: w0
 
    integer :: istep, trash, iunit, nspin, time_steps, istart, iend, ntiter, lmax, no_e, ie, idir, ispin
    real(real64) :: dt, dump, vv(3)
    type(kick_t) :: kick
    real(real64), allocatable :: dd(:,:)
    real(real64), allocatable :: sps(:), spc(:), racc(:)
    complex(real64), allocatable :: dipole(:), ddipole(:)
    type(batch_t) :: acc_batch, sps_batch, spc_batch
    type(unit_system_t) :: file_units
 
    push_sub(spectrum_hs_from_mult)
 
    iunit = io_open('multipoles', namespace, action='read', status='old', die=.false.)
    if (iunit == -1) then
      iunit = io_open('td.general/multipoles', namespace, action='read', status='old')
    end if
    call spectrum_mult_info(namespace, iunit, nspin, kick, time_steps, dt, file_units, lmax=lmax)
    call spectrum_fix_time_limits(spectrum, time_steps, dt, istart, iend, ntiter)
 
    if (spectrum%energy_step <= m_zero) spectrum%energy_step = m_two * m_pi / (dt*time_steps)
 
    call io_skip_header(iunit)
 
    ! load dipole from file
    safe_allocate(dipole(0:time_steps))
    safe_allocate(ddipole(0:time_steps))
    safe_allocate(dd(1:3, 1:nspin))
 
    vv(1:3) = vec(1:3) / norm2(vec(1:3))
 
    do istep = 1, time_steps
      dd = m_zero
      read(iunit, *) trash, dump, (dump, (dd(idir, ispin), idir = 1, kick%dim), ispin = 1, nspin)
      select case (pol)
      case ('x')
        dipole(istep) = -sum(dd(1, :))
      case ('y')
        dipole(istep) = -sum(dd(2, :))
      case ('z')
        dipole(istep) =  sum(dd(3, :))
      case ('+')
        dipole(istep) = -sum(cmplx(dd(1, :), dd(2, :), real64)) / sqrt(m_two)
      case ('-')
        dipole(istep) = -sum(cmplx(dd(1, :), -dd(2, :), real64)) / sqrt(m_two)
      case ('v')
        dipole(istep) = -sum(vv(1)*dd(1, :) + vv(2)*dd(2, :) + vv(3)*dd(3, :))
      end select
      dipole(istep) = units_to_atomic(units_out%length, dipole(istep))
    end do
    safe_deallocate_a(dd)
    dipole(0) = dipole(1)
    call io_close(iunit)
 
    ! we now calculate the acceleration.
    ddipole(0) = m_zero
    do istep = 1, time_steps - 1
      ddipole(istep) = (dipole(istep - 1) + dipole(istep + 1) - m_two * dipole(istep)) / dt**2
    end do
    call interpolate( dt*(/ -3, -2, -1 /),   &
      ddipole(time_steps - 3:time_steps - 1), &
      m_zero, &
      ddipole(time_steps))
 
    if (present(w0)) then
 
      call spectrum_hsfunction_init(dt, istart, iend, time_steps, ddipole)
      call spectrum_hs(spectrum, namespace, out_file, pol, w0)
      call spectrum_hsfunction_end()
 
    else
 
      safe_allocate(racc(0:time_steps))
      racc = real(ddipole, real64)
 
      no_e = spectrum_nenergy_steps(spectrum)
      safe_allocate(sps(1:no_e))
      safe_allocate(spc(1:no_e))
      sps = m_zero
      spc = m_zero
 
      call batch_init(acc_batch, racc)
      call batch_init(sps_batch, sps)
      call batch_init(spc_batch, spc)
 
      call spectrum_fourier_transform(spectrum%method, spectrum_transform_cos, spectrum%noise, &
        istart + 1, iend + 1, m_zero, dt, acc_batch, spectrum%min_energy, spectrum%max_energy, spectrum%energy_step, spc_batch)
      call spectrum_fourier_transform(spectrum%method, spectrum_transform_sin, spectrum%noise, &
        istart + 1, iend + 1, m_zero, dt, acc_batch, spectrum%min_energy, spectrum%max_energy, spectrum%energy_step, sps_batch)
 
      do ie = 1, no_e
        sps(ie) = (sps(ie)**2 + spc(ie)**2)
      end do
 
      call spectrum_hs_output(spectrum, namespace, out_file, pol, no_e, sps)
 
      call acc_batch%end()
      call sps_batch%end()
      call spc_batch%end()
 
      safe_deallocate_a(racc)
 
    end if
 
    safe_deallocate_a(dipole)
    safe_deallocate_a(ddipole)
 
    pop_sub(spectrum_hs_from_mult)
  end subroutine spectrum_hs_from_mult
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine spectrum_hs_from_acc(spectrum, namespace, out_file, pol, vec, w0)
    type(spectrum_t),  intent(inout) :: spectrum
    type(namespace_t), intent(in)    :: namespace
    character(len=*),  intent(in)    :: out_file
    character,         intent(in)    :: pol
    real(real64),      intent(in)    :: vec(:)
    real(real64),   optional, intent(in)    :: w0
 
    integer :: istep, jj, iunit, time_steps, istart, iend, ntiter, ierr, no_e, ie
    real(real64) :: dt, aa(3), vv(3)
    complex(real64), allocatable :: acc(:)
    real(real64), allocatable :: racc(:), sps(:), spc(:)
    type(batch_t) :: acc_batch, sps_batch, spc_batch
 
    push_sub(spectrum_hs_from_acc)
 
    call spectrum_tdfile_info(namespace, 'acceleration', iunit, time_steps, dt)
    call spectrum_fix_time_limits(spectrum, time_steps, dt, istart, iend, ntiter)
 
    if (spectrum%energy_step <= m_zero) spectrum%energy_step = m_two * m_pi / (dt*time_steps)
 
    ! load dipole from file
    safe_allocate(acc(0:time_steps))
    acc = m_zero
    vv = vec / norm2(vec(:))
    call io_skip_header(iunit)
 
    do istep = 1, time_steps
      aa = m_zero
      read(iunit, '(28x,e20.12)', advance = 'no', iostat = ierr) aa(1)
      jj = 2
      do while((ierr == 0) .and. (jj <= 3))
        read(iunit, '(e20.12)', advance = 'no', iostat = ierr) aa(jj)
        jj = jj + 1
      end do
      select case (pol)
      case ('x')
        acc(istep) = aa(1)
      case ('y')
        acc(istep) = aa(2)
      case ('z')
        acc(istep) = aa(3)
      case ('+')
        acc(istep) = cmplx(aa(1), aa(2), real64) / sqrt(m_two)
      case ('-')
        acc(istep) = cmplx(aa(1), -aa(2), real64) / sqrt(m_two)
      case ('v')
        acc(istep) = vv(1)*aa(1) + vv(2)*aa(2) + vv(3)*aa(3)
      end select
      acc(istep) = units_to_atomic(units_out%acceleration, acc(istep))
    end do
    close(iunit)
 
    if (present(w0)) then
 
      call spectrum_hsfunction_init(dt, istart, iend, time_steps, acc)
      call spectrum_hs(spectrum, namespace, out_file, pol, w0)
      call spectrum_hsfunction_end()
 
    else
 
      safe_allocate(racc(0:time_steps))
      racc = real(acc, real64)
 
      no_e = spectrum_nenergy_steps(spectrum)
      safe_allocate(sps(1:no_e))
      safe_allocate(spc(1:no_e))
      sps = m_zero
      spc = m_zero
 
      call batch_init(acc_batch, racc)
      call batch_init(sps_batch, sps)
      call batch_init(spc_batch, spc)
 
      call spectrum_fourier_transform(spectrum%method, spectrum_transform_cos, spectrum%noise, &
        istart + 1, iend + 1, m_zero, dt, acc_batch, spectrum%min_energy, &
        spectrum%max_energy, spectrum%energy_step, spc_batch)
      call spectrum_fourier_transform(spectrum%method, spectrum_transform_sin, spectrum%noise, &
        istart + 1, iend + 1, m_zero, dt, acc_batch, spectrum%min_energy, &
        spectrum%max_energy, spectrum%energy_step, sps_batch)
 
      do ie = 1, no_e
        sps(ie) = (sps(ie)**2 + spc(ie)**2)
      end do
 
      call spectrum_hs_output(spectrum, namespace, out_file, pol, no_e, sps)
 
      call acc_batch%end()
      call sps_batch%end()
      call spc_batch%end()
 
      safe_deallocate_a(racc)
 
    end if
 
    safe_deallocate_a(acc)
    pop_sub(spectrum_hs_from_acc)
  end subroutine spectrum_hs_from_acc
  ! ---------------------------------------------------------
 
  ! ---------------------------------------------------------
  subroutine spectrum_hs_from_current(spectrum, namespace, out_file, pol, vec, w0)
    type(spectrum_t),  intent(inout) :: spectrum
    type(namespace_t), intent(in)    :: namespace
    character(len=*),  intent(in)    :: out_file
    character,         intent(in)    :: pol
    real(real64),      intent(in)    :: vec(:)
    real(real64),   optional, intent(in)    :: w0
 
    integer :: istep, jj, iunit, time_steps, istart, iend, ntiter, ierr, no_e, ie
    real(real64) :: dt, cc(3), vv(3)
    complex(real64), allocatable :: cur(:)
    real(real64), allocatable :: rcur(:), sps(:), spc(:)
    type(batch_t) :: cur_batch, sps_batch, spc_batch
 
    push_sub(spectrum_hs_from_current)
 
    call spectrum_tdfile_info(namespace, 'total_current', iunit, time_steps, dt)
    call spectrum_fix_time_limits(spectrum, time_steps, dt, istart, iend, ntiter)
 
    if (spectrum%energy_step <= m_zero) spectrum%energy_step = m_two * m_pi / (dt * time_steps)
 
    ! load dipole from file
    safe_allocate(cur(0:time_steps))
    cur = m_zero
    vv = vec / norm2(vec(:))
    call io_skip_header(iunit)
 
    do istep = 1, time_steps
      cc = m_zero
      read(iunit, '(28x,e20.12)', advance = 'no', iostat = ierr) cc(1)
      jj = 2
      do while((ierr == 0) .and. (jj <= 3))
        read(iunit, '(e20.12)', advance = 'no', iostat = ierr) cc(jj)
        jj = jj + 1
      end do
      select case (pol)
      case ('x')
        cur(istep) = cc(1)
      case ('y')
        cur(istep) = cc(2)
      case ('z')
        cur(istep) = cc(3)
      case ('+')
        cur(istep) = cmplx(cc(1), cc(2), real64) / sqrt(m_two)
      case ('-')
        cur(istep) = cmplx(cc(1), -cc(2), real64) / sqrt(m_two)
      case ('v')
        cur(istep) = vv(1)*cc(1) + vv(2)*cc(2) + vv(3)*cc(3)
      end select
      cur(istep) = units_to_atomic(units_out%velocity, cur(istep))
    end do
    close(iunit)
 
    if (present(w0)) then
 
      call spectrum_hsfunction_init(dt, istart, iend, time_steps, cur)
      call spectrum_hs(spectrum, namespace, out_file, pol, w0)
      call spectrum_hsfunction_end()
 
    else
 
      safe_allocate(rcur(0:time_steps))
      rcur = real(cur, real64)
 
      no_e = spectrum_nenergy_steps(spectrum)
      safe_allocate(sps(1:no_e))
      safe_allocate(spc(1:no_e))
      sps = m_zero
      spc = m_zero
 
      call batch_init(cur_batch, rcur)
      call batch_init(sps_batch, sps)
      call batch_init(spc_batch, spc)
 
      call spectrum_fourier_transform(spectrum%method, spectrum_transform_cos, spectrum%noise, &
        istart + 1, iend + 1, m_zero, dt, cur_batch, spectrum%min_energy, spectrum%max_energy, spectrum%energy_step, spc_batch)
      call spectrum_fourier_transform(spectrum%method, spectrum_transform_sin, spectrum%noise, &
        istart + 1, iend + 1, m_zero, dt, cur_batch, spectrum%min_energy, spectrum%max_energy, spectrum%energy_step, sps_batch)
 
      do ie = 1, no_e
        sps(ie) = (sps(ie)**2 + spc(ie)**2) * ((ie-1) * spectrum%energy_step + spectrum%min_energy)**2
      end do
 
      call spectrum_hs_output(spectrum, namespace, out_file, pol, no_e, sps)
 
      call cur_batch%end()
      call sps_batch%end()
      call spc_batch%end()
 
      safe_deallocate_a(rcur)
 
    end if
 
    safe_deallocate_a(cur)
    pop_sub(spectrum_hs_from_current)
  end subroutine spectrum_hs_from_current
 
 
  ! ---------------------------------------------------------
  subroutine spectrum_hs(spectrum, namespace, out_file, pol, w0)
    type(spectrum_t),  intent(inout) :: spectrum
    type(namespace_t), intent(in)    :: namespace
    character(len=*),  intent(in)    :: out_file
    character,         intent(in)    :: pol
    real(real64),  optional,  intent(in)    :: w0
 
    integer :: iunit, no_e, ie
    real(real64)   :: omega, hsval, xx
    real(real64), allocatable :: sp(:)
 
    push_sub(spectrum_hs)
 
    if (present(w0)) then
 
      iunit = io_open(trim(out_file) // "." // trim(pol), namespace, action='write')
      write(iunit, '(a1,a20,a20)') '#', str_center("w", 20), str_center("H(w)", 20)
      write(iunit, '(a1,a20,a20)') '#', &
        str_center('['//trim(units_abbrev(units_out%energy)) // ']', 20), &
        str_center('[('//trim(units_abbrev(units_out%length))//'/' &
        //trim(units_abbrev(units_out%time**2)), 20)
 
      ! output
      omega = w0
      do while(omega <= spectrum%max_energy)
        call spectrum_hsfunction_min(namespace, omega - w0, omega + w0, xx, hsval)
 
        write(iunit, '(1x,2e20.8)') units_from_atomic(units_out%energy, xx), &
          units_from_atomic((units_out%length / units_out%time)**2, -hsval)
 
        ! 2 * w0 because we assume that there are only odd peaks.
        omega = omega + 2 * w0
      end do
      call io_close(iunit)
 
    else
      no_e = spectrum_nenergy_steps(spectrum)
      safe_allocate(sp(1:no_e))
      sp = m_zero
 
      do ie = 1, no_e
        call hsfunction((ie-1) * spectrum%energy_step + spectrum%min_energy, sp(ie))
        sp(ie) = -sp(ie)
      end do
 
      call spectrum_hs_output(spectrum, namespace, out_file, pol, no_e, sp)
 
      safe_deallocate_a(sp)
 
    end if
 
    pop_sub(spectrum_hs)
  end subroutine spectrum_hs
  ! ---------------------------------------------------------
 
 
  subroutine spectrum_hs_output(spectrum, namespace, out_file, pol, no_e, sp)
    type(spectrum_t),  intent(inout) :: spectrum
    type(namespace_t), intent(in)    :: namespace
    character(len=*),  intent(in)    :: out_file
    character,         intent(in)    :: pol
    integer,           intent(in)    :: no_e
    real(real64),      intent(in)    :: sp(:)
 
    integer :: iunit, ie
 
    push_sub(spectrum_hs_output)
 
    ! output
    if (trim(out_file) /= '-') then
      iunit = io_open(trim(out_file) // "." // trim(pol), namespace, action='write')
      write(iunit, '(a1,a20,a20)') '#', str_center("w", 20), str_center("H(w)", 20)
 
      write(iunit, '(a1,a20,a20)') &
        '#', str_center('['//trim(units_abbrev(units_out%energy)) // ']', 20), &
        str_center('[('//trim(units_abbrev(units_out%length))//'/' &
        //trim(units_abbrev(units_out%time**2)), 20)
 
      do ie = 1, no_e
        write(iunit, '(2e15.6)') units_from_atomic(units_out%energy, (ie-1) * spectrum%energy_step + spectrum%min_energy), &
          units_from_atomic((units_out%length / units_out%time)**2, sp(ie))
      end do
 
      call io_close(iunit)
    end if
 
    pop_sub(spectrum_hs_output)
  end subroutine spectrum_hs_output
 
 
  ! ---------------------------------------------------------
  subroutine spectrum_mult_info(namespace, iunit, nspin, kick, time_steps, dt, file_units, lmax)
    type(namespace_t),   intent(in)  :: namespace
    integer,             intent(in)  :: iunit
    integer,             intent(out) :: nspin
    type(kick_t),        intent(out) :: kick
    integer,             intent(out) :: time_steps
    real(real64),        intent(out) :: dt
    type(unit_system_t), intent(out) :: file_units
    integer,   optional, intent(out) :: lmax
 
    integer :: ii
    character(len=100) :: line
 
    push_sub(spectrum_mult_info)
 
    assert(iunit /= -1)
 
    rewind(iunit)
    read(iunit,*)
    read(iunit,*)
    read(iunit, '(15x,i2)') nspin
    if (present(lmax)) then
      read(iunit, '(15x,i2)') lmax
    end if
    call kick_read(kick, iunit, namespace)
    read(iunit, '(a)') line
    read(iunit, '(a)') line
    call io_skip_header(iunit)
 
    ! Figure out the units of the file
    ii = index(line,'eV')
    if (ii /= 0) then
      call unit_system_get(file_units, units_eva)
    else
      call unit_system_get(file_units, units_atomic)
    end if
 
    call spectrum_count_time_steps(namespace, iunit, time_steps, dt)
    dt = units_to_atomic(file_units%time, dt) ! units_out is OK
 
    pop_sub(spectrum_mult_info)
  end subroutine spectrum_mult_info
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine spectrum_count_time_steps(namespace, iunit, time_steps, dt)
    type(namespace_t), intent(in)  :: namespace
    integer,           intent(in)  :: iunit
    integer,           intent(out) :: time_steps
    real(real64),      intent(out) :: dt
 
    real(real64) :: t1, t2, dummy
    integer :: trash
 
    push_sub(count_time_steps)
 
    ! count number of time_steps
    time_steps = 0
    t1 = m_zero
    t2 = m_zero
    do
      read(iunit, *, end=100) trash, dummy
      time_steps = time_steps + 1
      if (time_steps == 1) t1 = dummy
      if (time_steps == 2) t2 = dummy
    end do
100 continue
    dt = (t2 - t1)
    time_steps = time_steps - 1
 
    if (time_steps < 3) then
      message(1) = "Empty file?"
      call messages_fatal(1, namespace=namespace)
    end if
 
    pop_sub(count_time_steps)
  end subroutine spectrum_count_time_steps
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine spectrum_cross_section_info(namespace, iunit, nspin, kick, energy_steps, dw)
    type(namespace_t), intent(in)  :: namespace
    integer,           intent(in)  :: iunit
    integer,           intent(out) :: nspin
    type(kick_t),      intent(out) :: kick
    integer,           intent(out) :: energy_steps
    real(real64),      intent(out) :: dw
 
    real(real64) :: dummy, e1, e2
 
    push_sub(spectrum_cross_section_info)
 
    ! read in number of spin components
    read(iunit, '(15x,i2)') nspin
    call kick_read(kick, iunit, namespace)
    call io_skip_header(iunit)
 
    ! count number of time_steps
    energy_steps = 0
    do
      read(iunit, *, end=100) dummy
      energy_steps = energy_steps + 1
      if (energy_steps == 1) e1 = dummy
      if (energy_steps == 2) e2 = dummy
    end do
100 continue
    dw = units_to_atomic(units_out%energy, e2 - e1)
 
    if (energy_steps < 3) then
      message(1) = "Empty multipole file?"
      call messages_fatal(1, namespace=namespace)
    end if
 
    pop_sub(spectrum_cross_section_info)
  end subroutine spectrum_cross_section_info
 
 
  ! ---------------------------------------------------------
  subroutine spectrum_tdfile_info(namespace, fname, iunit, time_steps, dt)
    type(namespace_t), intent(in)  :: namespace
    character(len=*),  intent(in)  :: fname
    integer,           intent(out) :: iunit, time_steps
    real(real64),      intent(out) :: dt
 
    integer :: trash
    real(real64) :: t1, t2, dummy
    character(len=256) :: filename
 
    push_sub(spectrum_tdfile_info)
 
    ! open files
    filename = trim('td.general/')//trim(fname)
    iunit = io_open(filename, namespace, action='read', status='old', die=.false.)
 
    if (iunit == -1) then
      filename = trim('./')//trim(fname)
      iunit = io_open(filename, namespace, action='read', status='old')
    end if
 
 
    ! read in dipole
    call io_skip_header(iunit)
 
    ! count number of time_steps
    time_steps = 0
    do
      read(iunit, *, end=100) trash, dummy
      time_steps = time_steps + 1
      if (time_steps == 1) t1 = dummy
      if (time_steps == 2) t2 = dummy
    end do
100 continue
    dt = units_to_atomic(units_out%time, t2 - t1) ! units_out is OK
    time_steps = time_steps - 1
 
    if (time_steps < 3) then
      message(1) = "Empty file?"
      call messages_fatal(1, namespace=namespace)
    end if
 
    rewind(iunit)
    pop_sub(spectrum_tdfile_info)
  end subroutine spectrum_tdfile_info
 
 
  ! ---------------------------------------------------------
  subroutine spectrum_fix_time_limits(spectrum, time_steps, dt, istart, iend, ntiter)
    type(spectrum_t), intent(inout) :: spectrum
    integer,          intent(in)    :: time_steps
    real(real64),     intent(in)    :: dt
    integer,          intent(out)   :: istart, iend, ntiter
 
    real(real64) :: ts, te, dummy
 
    push_sub(spectrum_fix_time_limits)
 
    ts = m_zero
    te = time_steps * dt
 
    if (spectrum%start_time < ts) spectrum%start_time = ts
    if (spectrum%start_time > te) spectrum%start_time = te
    if (spectrum%end_time   > te .or. spectrum%end_time <= m_zero) spectrum%end_time = te
    if (spectrum%end_time   < ts) spectrum%end_time   = ts
 
    if (spectrum%end_time < spectrum%start_time) then
      dummy = spectrum%end_time ! swap
      spectrum%end_time = spectrum%start_time
      spectrum%start_time = dummy
    end if
    istart = nint(spectrum%start_time / dt)
    iend = nint(spectrum%end_time / dt)
    ntiter = iend - istart + 1
 
    ! Get default damp factor
    if (spectrum%damp /= spectrum_damp_none .and. spectrum%damp /= spectrum_damp_polynomial &
      .and. is_close(spectrum%damp_factor, -m_one)) then
      select case (spectrum%damp)
      case (spectrum_damp_lorentzian)
        spectrum%damp_factor =  -log(0.0001_real64)/(spectrum%end_time-spectrum%start_time)
      case (spectrum_damp_gaussian)
        spectrum%damp_factor =  sqrt(-log(0.0001_real64)/(spectrum%end_time-spectrum%start_time)**2)
      end select
    end if
 
 
    pop_sub(spectrum_fix_time_limits)
  end subroutine spectrum_fix_time_limits
 
  ! -------------------------------------------------------
 
  subroutine spectrum_signal_damp(damp_type, damp_factor, time_start, time_end, t0, time_step, time_function)
    integer,            intent(in)    :: damp_type
    real(real64),       intent(in)    :: damp_factor
    integer,            intent(in)    :: time_start
    integer,            intent(in)    :: time_end
    real(real64),       intent(in)    :: t0
    real(real64),       intent(in)    :: time_step
    type(batch_t),      intent(inout) :: time_function
 
    integer :: itime, ii
    real(real64)   :: time
    real(real64), allocatable :: weight(:)
 
    push_sub(signal_damp)
 
    assert(time_function%status() == batch_not_packed)
 
    safe_allocate(weight(time_start:time_end))
 
    do itime = time_start, time_end
      time = time_step*(itime-1)
 
      ! Gets the damp function
      select case (damp_type)
      case (spectrum_damp_none)
        weight(itime) = m_one
      case (spectrum_damp_lorentzian)
        if (time < t0) then
          weight(itime) = m_one
        else
          weight(itime) = exp(-(time - t0)*damp_factor)
        end if
      case (spectrum_damp_polynomial)
        if (time < t0) then
          weight(itime) = m_one
        else
          weight(itime) = m_one - m_three*((time - t0) / (time_step * (time_end - 1) - t0))**2 + &
            m_two * ((time - t0) / (time_step * (time_end - 1) - t0))**3
        end if
      case (spectrum_damp_gaussian)
        if (time < t0) then
          weight(itime) = m_one
        else
          weight(itime) = exp(-(time - t0)**2*damp_factor**2)
        end if
      case (spectrum_damp_sin)
        if (time < t0) then
          weight(itime) = m_one
        else
          weight(itime) = sin(-(time - t0)*m_pi/(time_end+t0))
        end if
      end select
    end do
 
    if (time_function%type() == type_cmplx) then
      do ii = 1, time_function%nst_linear
        do itime = time_start, time_end
          time_function%zff_linear(itime, ii) = weight(itime)*time_function%zff_linear(itime, ii)
        end do
      end do
    else
      do ii = 1, time_function%nst_linear
        do itime = time_start, time_end
          time_function%dff_linear(itime, ii) = weight(itime)*time_function%dff_linear(itime, ii)
        end do
      end do
    end if
 
    safe_deallocate_a(weight)
 
    pop_sub(signal_damp)
 
  end subroutine spectrum_signal_damp
  ! -------------------------------------------------------
 
  ! -------------------------------------------------------
  subroutine spectrum_fourier_transform(method, transform, noise, time_start, time_end, t0, time_step, time_function, &
    energy_start, energy_end, energy_step, energy_function)
    integer,                  intent(in)    :: method
    integer,                  intent(in)    :: transform
    real(real64),             intent(in)    :: noise
    integer,                  intent(in)    :: time_start
    integer,                  intent(in)    :: time_end
    real(real64),             intent(in)    :: t0
    real(real64),             intent(in)    :: time_step
    type(batch_t),            intent(in)    :: time_function
    real(real64),             intent(in)    :: energy_start
    real(real64),             intent(in)    :: energy_end
    real(real64),             intent(in)    :: energy_step
    type(batch_t),            intent(inout) :: energy_function
 
    integer :: itime, ienergy, ii, energy_steps
    real(real64)   :: energy, sinz, cosz
    complex(real64) :: ez, eidt
    type(compressed_sensing_t) :: cs
 
    push_sub(fourier_transform)
 
    assert(time_function%nst_linear == energy_function%nst_linear)
    assert(time_function%status() == energy_function%status())
    assert(time_function%status() == batch_not_packed)
    assert(time_function%type() == type_float)
    assert(energy_function%type() == type_float)
 
    energy_steps = nint((energy_end-energy_start) / energy_step) + 1
 
    select case (method)
 
    case (spectrum_fourier)
 
      do ienergy = 1, energy_steps
 
        energy = energy_step*(ienergy - 1) + energy_start
 
        do ii = 1, energy_function%nst_linear
          energy_function%dff_linear(ienergy, ii) = m_zero
        end do
 
        select case (transform)
 
          ! The sine and cosine transforms are computed as the real and imaginary part of the exponential.
          ! One can compute the exponential by successive multiplications, instead of calling the sine or
          ! cosine function at each time step.
        case (spectrum_transform_sin)
 
          eidt = exp(m_zi * energy * time_step)
          ez = exp(m_zi * energy * ((time_start-1)*time_step - t0))
          sinz = aimag(ez)
          do itime = time_start, time_end
            do ii = 1, time_function%nst_linear
              energy_function%dff_linear(ienergy, ii) = &
                energy_function%dff_linear(ienergy, ii) + &
                time_function%dff_linear(itime, ii) * sinz
            end do
            ez = ez * eidt
            sinz = aimag(ez)
          end do
 
        case (spectrum_transform_cos)
 
          eidt = exp(m_zi * energy * time_step)
          ez = exp(m_zi * energy * ( (time_start-1)*time_step - t0))
          cosz = real(ez, real64)
          do itime = time_start, time_end
            do ii = 1, time_function%nst_linear
              energy_function%dff_linear(ienergy, ii) = &
                energy_function%dff_linear(ienergy, ii) + &
                time_function%dff_linear(itime, ii) * cosz
            end do
            ez = ez * eidt
            cosz = real(ez, real64)
          end do
 
        case (spectrum_transform_laplace)
 
          eidt = exp(-energy * time_step)
          ez = exp(-energy * ((time_start - 1) * time_step - t0))
          do itime = time_start, time_end
            do ii = 1, time_function%nst_linear
              energy_function%dff_linear(ienergy, ii) = &
                energy_function%dff_linear(ienergy, ii) + &
                real( time_function%dff_linear(itime, ii) * ez, real64)
            end do
            ez = ez * eidt
          end do
        end select
 
        ! The total sum must be multiplied by time_step in order to get the integral.
        do ii = 1, time_function%nst_linear
          energy_function%dff_linear(ienergy, ii) = &
            energy_function%dff_linear(ienergy, ii) * time_step
        end do
 
 
      end do
 
    case (spectrum_compressed_sensing)
 
      call compressed_sensing_init(cs, transform, &
        time_end - time_start + 1, time_step, time_step*(time_start - 1) - t0, &
        energy_steps, energy_step, energy_start, noise)
 
      do ii = 1, time_function%nst_linear
        call compressed_sensing_spectral_analysis(cs, time_function%dff_linear(:, ii), &
          energy_function%dff_linear(:, ii))
      end do
 
      call compressed_sensing_end(cs)
 
    end select
 
    pop_sub(fourier_transform)
 
  end subroutine spectrum_fourier_transform
 
  ! ---------------------------------------------------------
  subroutine spectrum_sigma_diagonalize(namespace, sigma, nspin, energy_step, min_energy, energy_steps, kick)
    type(namespace_t),      intent(in) :: namespace
    real(real64),           intent(in) :: sigma(:, :, :, :)
    integer,                intent(in) :: nspin
    real(real64),           intent(in) :: energy_step, min_energy
    integer,                intent(in) :: energy_steps
    type(kick_t), optional, intent(in) :: kick
 
    integer :: is, idir, jdir, ie, info, out_file, out_file_t
    real(real64), allocatable :: work(:,:)
    complex(real64), allocatable :: w(:)
    character(len=20) :: header_string
    logical :: spins_singlet, spins_triplet, symmetrize
    real(real64), allocatable :: pp(:,:), pp2(:,:)
 
    push_sub(spectrum_sigma_diagonalize)
 
 
    !%Variable PropagationSpectrumSymmetrizeSigma
    !%Type logical
    !%Default .false.
    !%Section Utilities::oct-propagation_spectrum
    !%Description
    !% The polarizablity tensor has to be real and symmetric. Due to numerical accuracy,
    !% that is not extricly conserved when computing it from different time-propations.
    !% If <tt>PropagationSpectrumSymmetrizeSigma = yes</tt>, the polarizability tensor is
    !% symmetrized before its diagonalizied.
    !% This variable is only used if the cross_section_tensor is computed.
    !%End
    call parse_variable(namespace, 'PropagationSpectrumSymmetrizeSigma', .false., symmetrize)
    call messages_print_var_value('PropagationSpectrumSymmetrizeSigma', symmetrize, namespace=namespace)
 
    spins_singlet = .true.
    spins_triplet = .false.
    if (present(kick)) then
      select case (kick_get_type(kick))
      case (kick_spin_mode)
        spins_triplet = .true.
        spins_singlet = .false.
      case (kick_spin_density_mode)
        spins_triplet = .true.
      end select
    end if
 
    if (spins_singlet .and. spins_triplet) then
      out_file = io_open('cross_section_diagonal-sigma_s', namespace, action='write')
      out_file_t = io_open('cross_section_diagonal-sigma_t', namespace, action='write')
    else
      out_file = io_open('cross_section_diagonal-sigma', namespace, action='write')
    end if
 
    is = 1
    write(out_file, '(a1, a20)', advance = 'no') '#', str_center("Energy", 20)
    do idir = 1, 3
      write(out_file, '(a20)', advance = 'no') str_center("Real part", 20)
      if (.not. symmetrize) write(out_file, '(a20)', advance = 'no') str_center("Imaginary part", 20)
      do jdir = 1, 3
        write(header_string,'(a7,i1,a1,i1,a1,i1,a1)') 'vector(', idir, ',', jdir, ',', is, ')'
        write(out_file, '(a20)', advance = 'no') str_center(trim(header_string), 20)
      end do
    end do
    write(out_file, '(1x)')
    write(out_file, '(a1,a20)', advance = 'no') '#', str_center('[' // trim(units_abbrev(units_out%energy)) // ']', 20)
 
    do idir = 1, 3
      write(out_file, '(a20)', advance = 'no')  str_center('[' // trim(units_abbrev(units_out%length**2)) // ']', 20)
      if (.not. symmetrize) then
        write(out_file, '(a20)', advance = 'no')  str_center('[' // trim(units_abbrev(units_out%length**2)) // ']', 20)
      end if
      do jdir = 1, 3
        write(out_file, '(a20)', advance = 'no')  str_center('[ - ]', 20)
      end do
    end do
    write(out_file, '(1x)')
 
    if (spins_singlet .and. spins_triplet) then
      is = 2
      write(out_file_t, '(a1, a20)', advance = 'no') '#', str_center("Energy", 20)
      do idir = 1, 3
        write(out_file_t, '(a20)', advance = 'no') str_center("Real part", 20)
        if (.not. symmetrize) write(out_file_t, '(a20)', advance = 'no') str_center("Imaginary part", 20)
        do jdir = 1, 3
          write(header_string,'(a7,i1,a1,i1,a1,i1,a1)') 'vector(', idir, ',', jdir, ',', is, ')'
          write(out_file_t, '(a20)', advance = 'no') str_center(trim(header_string), 20)
        end do
      end do
      write(out_file_t, '(1x)')
      write(out_file_t, '(a1,a20)', advance = 'no') '#', str_center('[' // trim(units_abbrev(units_out%energy)) // ']', 20)
 
      do idir = 1, 3
        write(out_file_t, '(a20)', advance = 'no')  str_center('[' // trim(units_abbrev(units_out%length**2)) // ']', 20)
        if (.not. symmetrize) then
          write(out_file_t, '(a20)', advance = 'no')  str_center('[' // trim(units_abbrev(units_out%length**2)) // ']', 20)
        end if
        do jdir = 1, 3
          write(out_file_t, '(a20)', advance = 'no')  str_center('[ - ]', 20)
        end do
      end do
      write(out_file_t, '(1x)')
    end if
 
    safe_allocate(pp(1:3, 1:3))
    if (spins_triplet .and. spins_singlet) then
      safe_allocate(pp2(1:3, 1:3))
    end if
    safe_allocate(w(1:3))
    safe_allocate(work(1:3, 1:3))
    do ie = 1, energy_steps
 
      pp(:, :) = sigma(:, :, ie, 1)
      if (nspin >= 2) then
        if (spins_singlet .and. spins_triplet) then
          pp2(:, :) = pp(:, :) - sigma(:, :, ie, 2)
          pp(:, :)  = pp(:, :) + sigma(:, :, ie, 2)
        elseif (spins_triplet .and. .not. spins_singlet) then
          pp(:, :) = pp(:, :) - sigma(:, :, ie, 2)
        elseif (spins_singlet .and. .not. spins_triplet) then
          pp(:, :) = pp(:, :) + sigma(:, :, ie, 2)
        end if
      end if
 
      if (symmetrize) then
        do idir = 1, 3
          do jdir = idir + 1, 3
            pp(idir, jdir) = (pp(idir, jdir) + pp(jdir, idir)) / m_two
            pp(jdir, idir) = pp(idir, jdir)
          end do
        end do
      end if
 
      work(1:3, 1:3) = pp(1:3, 1:3)
      call lalg_eigensolve_nonh(3, work, w, err_code = info, sort_eigenvectors = .true.)
      ! Note that the cross-section elements do not have to be transformed to the proper units, since
      ! they have been read from the "cross_section_vector.x", where they are already in the proper units.
 
      write(out_file,'(e20.8)', advance = 'no') units_from_atomic(units_out%energy, ((ie-1) * energy_step + min_energy))
      do idir = 3, 1, -1
        if (symmetrize) then
          write(out_file,'(2e20.8)', advance = 'no') real(w(idir), real64)
        else
          write(out_file,'(2e20.8)', advance = 'no') w(idir)
        end if
 
        do jdir = 1, 3
          write(out_file,'(e20.8)', advance = 'no') work(jdir, idir)
        end do
      end do
      write(out_file, '(1x)')
 
      if (spins_singlet .and. spins_triplet) then
        if (symmetrize) then
          do idir = 1, 3
            do jdir = idir + 1, 3
              pp2(idir, jdir) = (pp2(idir, jdir) + pp2(jdir, idir)) / m_two
              pp2(jdir, idir) = pp2(idir, jdir)
            end do
          end do
        end if
        work(1:3, 1:3) = -pp2(1:3, 1:3)
        call lalg_eigensolve_nonh(3, work, w, err_code = info, sort_eigenvectors = .true.)
        ! Note that the cross-section elements do not have to be transformed to the proper units, since
        ! they have been read from the "cross_section_vector.x", where they are already in the proper units.
 
        write(out_file_t,'(e20.8)', advance = 'no') units_from_atomic(units_out%energy, (ie * energy_step + min_energy))
        do idir = 3, 1, -1
          if (symmetrize) then
            write(out_file_t,'(2e20.8)', advance = 'no') real(w(idir), real64)
          else
            write(out_file_t,'(2e20.8)', advance = 'no') w(idir)
          end if
 
          do jdir = 1, 3
            write(out_file_t,'(e20.8)', advance = 'no') work(jdir, idir)
          end do
        end do
        write(out_file_t, '(1x)')
      end if
    end do
 
    call io_close(out_file)
 
    safe_deallocate_a(pp)
    if (spins_triplet .and. spins_singlet) then
      safe_deallocate_a(pp2)
      call io_close(out_file_t)
    end if
    safe_deallocate_a(w)
    safe_deallocate_a(work)
 
    pop_sub(spectrum_sigma_diagonalize)
  end subroutine spectrum_sigma_diagonalize
 
  pure integer function spectrum_nenergy_steps(spectrum) result(no_e)
    type(spectrum_t), intent(in) :: spectrum
 
    no_e = nint((spectrum%max_energy-spectrum%min_energy) / spectrum%energy_step) + 1
  end function spectrum_nenergy_steps
 
  subroutine spectrum_write_info(spectrum, out_file)
    type(spectrum_t), intent(in) :: spectrum
    integer,          intent(in) :: out_file
 
    push_sub(spectrum_write_info)
 
    write(out_file, '(a,i4)')    '# PropagationSpectrumDampMode   = ', spectrum%damp
    write(out_file, '(a,f10.4)') '# PropagationSpectrumDampFactor = ', units_from_atomic(units_out%time**(-1), &
      spectrum%damp_factor)
    write(out_file, '(a,f10.4)') '# PropagationSpectrumStartTime  = ', units_from_atomic(units_out%time, spectrum%start_time)
    write(out_file, '(a,f10.4)') '# PropagationSpectrumEndTime    = ', units_from_atomic(units_out%time, spectrum%end_time)
    write(out_file, '(a,f10.4)') '# PropagationSpectrumMinEnergy  = ', units_from_atomic(units_out%energy, spectrum%min_energy)
    write(out_file, '(a,f10.4)') '# PropagationSpectrumMaxEnergy  = ', units_from_atomic(units_out%energy, spectrum%max_energy)
    write(out_file, '(a,f10.4)') '# PropagationSpectrumEnergyStep = ', units_from_atomic(units_out%energy, spectrum%energy_step)
 
    pop_sub(spectrum_write_info)
  end subroutine spectrum_write_info
 
end module spectrum_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: