td.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module td_oct_m
  use absorbing_boundaries_oct_m
  use boundaries_oct_m
  use calc_mode_par_oct_m
  use current_oct_m
  use classical_particle_oct_m
  use debug_oct_m
  use density_oct_m
  use energy_calc_oct_m
  use electrons_ground_state_oct_m
  use electron_space_oct_m
  use epot_oct_m
  use ext_partner_list_oct_m
  use forces_oct_m
  use gauge_field_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_oct_m
  use interaction_partner_oct_m
  use io_oct_m
  use ion_dynamics_oct_m
  use ions_oct_m
  use kick_oct_m
  use, intrinsic :: iso_fortran_env
  use lasers_oct_m
  use lda_u_oct_m
  use lda_u_io_oct_m
  use linked_list_oct_m
  use loct_oct_m
  use maxwell_boundary_op_oct_m
  use mesh_oct_m
  use messages_oct_m
  use mpi_oct_m
  use multicomm_oct_m
  use namespace_oct_m
  use output_oct_m
  use output_low_oct_m
  use parser_oct_m
  use pes_oct_m
  use photon_mode_mf_oct_m
  use photon_mode_oct_m
  use poisson_oct_m
  use potential_interpolation_oct_m
  use profiling_oct_m
  use propagator_oct_m
  use propagator_elec_oct_m
  use propagator_base_oct_m
  use restart_oct_m
  use scf_oct_m
  use scissor_oct_m
  use space_oct_m
  use states_abst_oct_m
  use states_elec_oct_m
  use states_elec_restart_oct_m
  use stress_oct_m
  use td_write_oct_m
  use types_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use v_ks_oct_m
  use varinfo_oct_m
  use walltimer_oct_m
  use xc_oct_m
 
  implicit none
 
  private
  public ::                    &
    td_t,                      &
    td_run,                    &
    td_run_init,               &
    td_init,                   &
    td_init_run,               &
    td_end,                    &
    td_end_run,                &
    td_write_iter,             &
    td_check_point,            &
    td_dump,                   &
    td_allocate_wavefunctions, &
    td_init_gaugefield,        &
    td_load_restart_from_gs,   &
    td_load_restart_from_td,   &
    td_init_with_wavefunctions,&
    td_get_from_scratch,       &
    td_set_from_scratch
 
  integer, parameter, public :: &
    EHRENFEST = 1,   &
    bo        = 2
 
  type td_t
    private
    type(propagator_base_t), public :: tr
    type(scf_t),             public :: scf
    type(ion_dynamics_t),    public :: ions_dyn
    real(real64),            public :: dt
    integer,                 public :: max_iter
    integer,                 public :: iter
    logical,                 public :: recalculate_gs
 
    type(pes_t),             public :: pesv
 
    integer,                 public :: dynamics
    integer,                 public :: energy_update_iter
    real(real64)                    :: scissor
 
    logical                         :: freeze_occ
    logical                         :: freeze_u
    integer                         :: freeze_orbitals
 
    logical                         :: from_scratch = .false.
 
    type(td_write_t), public        :: write_handler
    type(restart_t)                 :: restart_load
    type(restart_t)                 :: restart_dump
  end type td_t
 
 
contains
 
  subroutine td_run_init()
 
    push_sub(td_run_init)
 
    call calc_mode_par%set_parallelization(p_strategy_states, default = .true.)
 
    pop_sub(td_run_init)
  end subroutine td_run_init
 
  ! ---------------------------------------------------------
 
  subroutine td_init(td, namespace, space, gr, ions, st, ks, hm, ext_partners, outp)
    type(td_t),               intent(inout) :: td
    type(namespace_t),        intent(in)    :: namespace
    class(space_t),           intent(in)    :: space
    type(grid_t),             intent(in)    :: gr
    type(ions_t),             intent(inout) :: ions
    type(states_elec_t),      intent(in)    :: st
    type(v_ks_t),             intent(in)    :: ks
    type(hamiltonian_elec_t), intent(in)    :: hm
    type(partner_list_t),     intent(in)    :: ext_partners
    type(output_t),           intent(in)    :: outp
 
    integer :: default
    real(real64)   :: propagation_time
    type(lasers_t), pointer :: lasers
    logical :: symmetrize
 
    push_sub(td_init)
 
    if (hm%pcm%run_pcm) call messages_experimental("PCM for CalculationMode = td", namespace=namespace)
 
    symmetrize = hm%kpoints%use_symmetries .or. st%symmetrize_density
    call ion_dynamics_init(td%ions_dyn, namespace, ions, symmetrize, gr%symmetrizer)
 
    if (td%ions_dyn%ions_move()) then
      if (hm%kpoints%use_symmetries) then
        message(1) = "KPoints symmetries cannot be used with moving ions."
        message(2) = "Please set KPointsSymmetries = no."
        call messages_fatal(2, namespace=namespace)
      end if
      if (st%symmetrize_density) then
        message(1) = "Symmetrization of the density cannot be used with moving ions."
        message(2) = "Please set SymmetrizeDensity = no."
        call messages_fatal(2, namespace=namespace)
      end if
    end if
 
    td%iter = 0
 
    !%Variable TDTimeStep
    !%Type float
    !%Section Time-Dependent::Propagation
    !%Description
    !% The time-step for the time propagation. For most propagators you
    !% want to use the largest value that is possible without the
    !% evolution becoming unstable.
    !%
    !% While prior versions of Octopus used to have a default time step, however,
    !% this now needs to be systematically defined.
    !%End
 
    call parse_variable(namespace, 'TDTimeStep', -m_one, td%dt, unit = units_inp%time)
 
    if (td%dt <= m_zero) then
      write(message(1),'(a)') 'A positive value for TDTimeStep must be defined in the input file.'
      call messages_fatal(1, namespace=namespace)
    end if
 
    call messages_print_var_value('TDTimeStep', td%dt, unit = units_out%time, namespace=namespace)
 
 
    if (parse_is_defined(namespace, 'TDMaxSteps') .and. parse_is_defined(namespace, 'TDPropagationTime')) then
      call messages_write('You cannot set TDMaxSteps and TDPropagationTime at the same time')
      call messages_fatal(namespace=namespace)
    end if
 
    !%Variable TDPropagationTime
    !%Type float
    !%Section Time-Dependent::Propagation
    !%Description
    !% The length of the time propagation. You cannot set this variable
    !% at the same time as <tt>TDMaxSteps</tt>. By default this variable will
    !% not be used.
    !%
    !% The units for this variable are <math>\hbar</math>/Hartree (or <math>\hbar</math>/eV if you
    !% selected <tt>ev_angstrom</tt> as input units). The approximate conversions to
    !% femtoseconds are 1 fs = 41.34 <math>\hbar</math>/Hartree = 1.52 <math>\hbar</math>/eV.
    !%End
    call parse_variable(namespace, 'TDPropagationTime', -1.0_real64, propagation_time, unit = units_inp%time)
 
    call messages_obsolete_variable(namespace, 'TDMaximumIter', 'TDMaxSteps')
 
    !%Variable TDMaxSteps
    !%Type integer
    !%Default 1500
    !%Section Time-Dependent::Propagation
    !%Description
    !% Number of time-propagation steps that will be performed. You
    !% cannot use this variable together with <tt>TDPropagationTime</tt>.
    !%End
    default = 1500
    if (propagation_time > m_zero) default = nint(propagation_time/td%dt)
    call parse_variable(namespace, 'TDMaxSteps', default, td%max_iter)
 
    if (propagation_time <= m_zero) propagation_time = td%dt*td%max_iter
 
    call messages_print_var_value('TDPropagationTime', propagation_time, unit = units_out%time, namespace=namespace)
    call messages_print_var_value('TDMaxSteps', td%max_iter, namespace=namespace)
 
    if (td%max_iter < 1) then
      write(message(1), '(a,i6,a)') "Input: '", td%max_iter, "' is not a valid value for TDMaxSteps."
      message(2) = '(TDMaxSteps <= 1)'
      call messages_fatal(2, namespace=namespace)
    end if
 
    td%iter = 0
 
    td%dt = td%dt
 
    lasers => list_get_lasers(ext_partners)
 
    ! now the photoelectron stuff
    call pes_init(td%pesv, namespace, space, gr, gr%box, st, outp%restart_write_interval, hm%kpoints, &
      hm%abs_boundaries, ext_partners, td%max_iter, td%dt)
 
    !%Variable TDDynamics
    !%Type integer
    !%Default ehrenfest
    !%Section Time-Dependent::Propagation
    !%Description
    !% Type of dynamics to follow during a time propagation.
    !% For BO, you must set <tt>MoveIons = yes</tt>.
    !%Option ehrenfest 1
    !% Ehrenfest dynamics.
    !%Option bo 2
    !% Born-Oppenheimer (Experimental).
    !%End
 
    call parse_variable(namespace, 'TDDynamics', ehrenfest, td%dynamics)
    if (.not. varinfo_valid_option('TDDynamics', td%dynamics)) call messages_input_error(namespace, 'TDDynamics')
    call messages_print_var_option('TDDynamics', td%dynamics, namespace=namespace)
    if (td%dynamics .ne. ehrenfest) then
      if (.not. td%ions_dyn%is_active()) then
        message(1) = "TDDynamics=bo can only be used if MoveIons=yes or CellDynamics=yes"
        call messages_fatal(1, namespace=namespace)
      end if
    end if
 
    !%Variable RecalculateGSDuringEvolution
    !%Type logical
    !%Default no
    !%Section Time-Dependent::Propagation
    !%Description
    !% In order to calculate some information about the system along the
    !% evolution (e.g. projection onto the ground-state KS determinant,
    !% projection of the TDKS spin-orbitals onto the ground-state KS
    !% spin-orbitals), the ground-state KS orbitals are needed. If the
    !% ionic potential changes -- that is, the ions move -- one may want
    !% to recalculate the ground state. You may do this by setting this
    !% variable.
    !%
    !% The recalculation is not done every time step, but only every
    !% <tt>RestartWriteInterval</tt> time steps.
    !%End
    call parse_variable(namespace, 'RecalculateGSDuringEvolution', .false., td%recalculate_gs)
    if (hm%lda_u_level /= dft_u_none .and. td%recalculate_gs) then
      call messages_not_implemented("DFT+U with RecalculateGSDuringEvolution=yes", namespace=namespace)
    end if
 
    !%Variable TDScissor
    !%Type float
    !%Default 0.0
    !%Section Time-Dependent
    !%Description
    !% (experimental) If set, a scissor operator will be applied in the
    !% Hamiltonian, shifting the excitation energies by the amount
    !% specified. By default, it is not applied.
    !%End
    call parse_variable(namespace, 'TDScissor', m_zero, td%scissor)
    call messages_print_var_value('TDScissor', td%scissor, namespace=namespace)
 
    call propagator_elec_init(gr, namespace, st, td%tr, hm%ks_pot, td%ions_dyn%is_active() .and.&
      list_has_gauge_field(ext_partners), family_is_mgga_with_exc(ks%xc), td%ions_dyn%cell_relax())
 
    if (associated(lasers) .and. mpi_grp_is_root(mpi_world)) then
      call messages_print_with_emphasis(msg="Time-dependent external fields", namespace=namespace)
      call laser_write_info(lasers%lasers, dt=td%dt, max_iter=td%max_iter, namespace=namespace)
      call messages_print_with_emphasis(namespace=namespace)
    end if
 
    !%Variable TDEnergyUpdateIter
    !%Type integer
    !%Section Time-Dependent::Propagation
    !%Description
    !% This variable controls after how many iterations Octopus
    !% updates the total energy during a time-propagation run. For
    !% iterations where the energy is not updated, the last calculated
    !% value is reported. If you set this variable to 1, the energy
    !% will be calculated in each step.
    !%End
 
    default = 10
    call parse_variable(namespace, 'TDEnergyUpdateIter', default, td%energy_update_iter)
 
    if (gr%der%boundaries%spiralBC .and. hm%ep%reltype == spin_orbit) then
      message(1) = "Generalized Bloch theorem cannot be used with spin-orbit coupling."
      call messages_fatal(1, namespace=namespace)
    end if
 
    if (gr%der%boundaries%spiralBC) then
      if (any(abs(hm%kick%easy_axis(1:2)) > m_epsilon)) then
        message(1) = "Generalized Bloch theorem cannot be used for an easy axis not along the z direction."
        call messages_fatal(1, namespace=namespace)
      end if
    end if
 
    !%Variable TDFreezeOrbitals
    !%Type integer
    !%Default 0
    !%Section Time-Dependent
    !%Description
    !% (Experimental) You have the possibility of "freezing" a number of orbitals during a time-propagation.
    !% The Hartree and exchange-correlation potential due to these orbitals (which
    !% will be the lowest-energy ones) will be added during the propagation, but the orbitals
    !% will not be propagated.
    !%Option sae -1
    !% Single-active-electron approximation. This option is only valid for time-dependent
    !% calculations (<tt>CalculationMode = td</tt>). Also, the nuclei should not move.
    !% The idea is that all orbitals except the last one are frozen. The orbitals are to
    !% be read from a previous ground-state calculation. The active orbital is then treated
    !% as independent (whether it contains one electron or two) -- although it will
    !% feel the Hartree and exchange-correlation potentials from the ground-state electronic
    !% configuration.
    !%
    !% It is almost equivalent to setting <tt>TDFreezeOrbitals = N-1</tt>, where <tt>N</tt> is the number
    !% of orbitals, but not completely.
    !%End
    call parse_variable(namespace, 'TDFreezeOrbitals', 0, td%freeze_orbitals)
 
    if (td%freeze_orbitals /= 0) then
      call messages_experimental('TDFreezeOrbitals', namespace=namespace)
 
      if (hm%lda_u_level /= dft_u_none) then
        call messages_not_implemented('TDFreezeOrbitals with DFT+U', namespace=namespace)
      end if
    end if
 
 
    pop_sub(td_init)
    nullify(lasers)
  end subroutine td_init
 
  ! ---------------------------------------------------------
  subroutine td_init_run(td, namespace, mc, gr, ions, st, ks, hm, ext_partners, outp, space, from_scratch)
    type(td_t),               intent(inout) :: td
    type(namespace_t),        intent(in)    :: namespace
    type(multicomm_t),        intent(inout) :: mc
    type(grid_t),             intent(inout) :: gr
    type(ions_t),             intent(inout) :: ions
    type(states_elec_t),      intent(inout) :: st
    type(v_ks_t),             intent(inout) :: ks
    type(hamiltonian_elec_t), intent(inout) :: hm
    type(partner_list_t),     intent(in)    :: ext_partners
    type(output_t),           intent(inout) :: outp
    type(electron_space_t),   intent(in)    :: space
    logical,                  intent(inout) :: from_scratch
    push_sub(td_init_run)
 
    ! NOTE: please do not change code in this function, but only in functions
    ! called from here because the logic of this function is replicated in the
    ! multisystem framework in different places
 
    call td_allocate_wavefunctions(td, namespace, mc, gr, ions, st, hm, space)
    call td_init_gaugefield(td, namespace, gr, st, ks, hm, ext_partners, space)
 
    td%from_scratch = from_scratch
 
    if (.not. td%from_scratch) then
      call td_load_restart_from_td(td, namespace, space, mc, gr, ext_partners, st, ks, hm, td%from_scratch)
      if (td%from_scratch) then
        message(1) = "Unable to read time-dependent restart information: Starting from scratch"
        call messages_warning(1, namespace=namespace)
      end if
    end if
 
    if (td%iter >= td%max_iter) then
      message(1) = "All requested iterations have already been done. Use FromScratch = yes if you want to redo them."
      call messages_info(1, namespace=namespace)
      call states_elec_deallocate_wfns(st)
      td%iter = td%iter + 1
      if ((td%ions_dyn%is_active()) .and. td%recalculate_gs) call restart_end(td%restart_load)
      pop_sub(td_init_run)
      return
    end if
 
    if (td%from_scratch) then
      call td_load_restart_from_gs(td, namespace, space, mc, gr, ext_partners, st, ks, hm)
    end if
 
    call td_init_with_wavefunctions(td, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, outp, td%from_scratch)
 
    pop_sub(td_init_run)
  end subroutine td_init_run
 
  ! ---------------------------------------------------------
  subroutine td_allocate_wavefunctions(td, namespace, mc, gr, ions, st, hm, space)
    type(td_t),               intent(inout) :: td
    type(namespace_t),        intent(in)    :: namespace
    type(multicomm_t),        intent(inout) :: mc
    type(grid_t),             intent(inout) :: gr
    type(ions_t),             intent(inout) :: ions
    type(states_elec_t),      intent(inout) :: st
    type(hamiltonian_elec_t), intent(inout) :: hm
    class(space_t),           intent(in)    :: space
 
    push_sub(td_allocate_wavefunctions)
 
    ! Allocate wavefunctions during time-propagation
    if (td%dynamics == ehrenfest) then
      !Note: this is not really clean to do this
      if (hm%lda_u_level /= dft_u_none .and. states_are_real(st)) then
        call lda_u_end(hm%lda_u)
        !complex wfs are required for Ehrenfest
        call states_elec_allocate_wfns(st, gr, type_cmplx, packed=.true.)
        call lda_u_init(hm%lda_u, namespace, space, hm%lda_u_level, gr, ions, st, mc, &
          hm%kpoints, hm%phase%is_allocated())
      else
        !complex wfs are required for Ehrenfest
        call states_elec_allocate_wfns(st, gr, type_cmplx, packed=.true.)
      end if
    else
      call states_elec_allocate_wfns(st, gr, packed=.true.)
      call scf_init(td%scf, namespace, gr, ions, st, mc, hm, space)
    end if
 
    pop_sub(td_allocate_wavefunctions)
  end subroutine td_allocate_wavefunctions
 
  ! ---------------------------------------------------------
  subroutine td_init_gaugefield(td, namespace, gr, st, ks, hm, ext_partners, space)
    type(td_t),               intent(inout) :: td
    type(namespace_t),        intent(in)    :: namespace
    type(grid_t),             intent(inout) :: gr
    type(states_elec_t),      intent(inout) :: st
    type(v_ks_t),             intent(inout) :: ks
    type(hamiltonian_elec_t), intent(inout) :: hm
    type(partner_list_t),     intent(in)    :: ext_partners
    class(space_t),           intent(in)    :: space
 
    type(gauge_field_t), pointer :: gfield
 
    push_sub(td_init_gaugefield)
 
    gfield => list_get_gauge_field(ext_partners)
    if(associated(gfield)) then
      if (gauge_field_is_used(gfield)) then
        !if the gauge field is applied, we need to tell v_ks to calculate the current
        call v_ks_calculate_current(ks, .true.)
 
        ! initialize the vector field and update the hamiltonian
        call gauge_field_init_vec_pot(gfield, st%qtot)
        call hm%update(gr, namespace, space, ext_partners, time = td%dt*td%iter)
 
      end if
    end if
 
    pop_sub(td_init_gaugefield)
  end subroutine td_init_gaugefield
 
  ! ---------------------------------------------------------
  subroutine td_end(td)
    type(td_t),               intent(inout) :: td
 
    push_sub(td_end)
 
    call pes_end(td%pesv)
    call propagator_elec_end(td%tr)  ! clean the evolution method
    call ion_dynamics_end(td%ions_dyn)
 
    if (td%dynamics == bo) call scf_end(td%scf)
 
    pop_sub(td_end)
  end subroutine td_end
 
  ! ---------------------------------------------------------
  subroutine td_end_run(td, st, hm)
    type(td_t),               intent(inout) :: td
    type(states_elec_t),      intent(inout) :: st
    type(hamiltonian_elec_t), intent(inout) :: hm
 
    push_sub(td_end_run)
 
    if (st%pack_states .and. hm%apply_packed()) call st%unpack()
 
    call restart_end(td%restart_dump)
    call td_write_end(td%write_handler)
 
    ! free memory
    call states_elec_deallocate_wfns(st)
    if ((td%ions_dyn%is_active()).and. td%recalculate_gs) call restart_end(td%restart_load)
 
    pop_sub(td_end_run)
  end subroutine td_end_run
 
  ! ---------------------------------------------------------
  subroutine td_run(td, namespace, mc, gr, ions, st, ks, hm, ext_partners, outp, space, from_scratch)
    type(td_t),               intent(inout) :: td
    type(namespace_t),        intent(in)    :: namespace
    type(multicomm_t),        intent(inout) :: mc
    type(grid_t),             intent(inout) :: gr
    type(ions_t),             intent(inout) :: ions
    type(states_elec_t),      intent(inout) :: st
    type(v_ks_t),             intent(inout) :: ks
    type(hamiltonian_elec_t), intent(inout) :: hm
    type(partner_list_t),     intent(in)    :: ext_partners
    type(output_t),           intent(inout) :: outp
    type(electron_space_t),   intent(in)    :: space
    logical,                  intent(inout) :: from_scratch
 
    logical                      :: stopping
    integer                      :: iter, scsteps
    real(real64)                 :: etime
 
    push_sub(td_run)
 
    etime = loct_clock()
    ! This is the time-propagation loop. It starts at t=0 and finishes
    ! at td%max_iter*dt. The index i runs from 1 to td%max_iter, and
    ! step "iter" means propagation from (iter-1)*dt to iter*dt.
    propagation: do iter = td%iter, td%max_iter
 
      stopping = clean_stop(mc%master_comm) .or. walltimer_alarm(mc%master_comm)
 
      call profiling_in("TIME_STEP")
 
      if (iter > 1) then
        if (((iter-1)*td%dt <= hm%kick%time) .and. (iter*td%dt > hm%kick%time)) then
          if (.not. hm%pcm%localf) then
            call kick_apply(space, gr, st, td%ions_dyn, ions, hm%kick, hm%psolver, hm%kpoints)
          else
            call kick_apply(space, gr, st, td%ions_dyn, ions, hm%kick, hm%psolver, hm%kpoints, pcm = hm%pcm)
          end if
          call td_write_kick(outp, namespace, space, gr, hm%kick, ions, iter)
          !We activate the sprial BC only after the kick,
          !to be sure that the first iteration corresponds to the ground state
          if (gr%der%boundaries%spiralBC) gr%der%boundaries%spiral = .true.
        end if
      end if
 
      ! time iterate the system, one time step.
      select case (td%dynamics)
      case (ehrenfest)
        call propagator_elec_dt(ks, namespace, space, hm, gr, st, td%tr, iter*td%dt, td%dt, iter, td%ions_dyn, &
          ions, ext_partners, mc, outp, td%write_handler, scsteps = scsteps, &
          update_energy = (mod(iter, td%energy_update_iter) == 0) .or. (iter == td%max_iter))
      case (bo)
        call propagator_elec_dt_bo(td%scf, namespace, space, gr, ks, st, hm, ions, ext_partners, mc, outp, iter, &
          td%dt, td%ions_dyn, scsteps)
      end select
 
      !Apply mask absorbing boundaries
      if (hm%abs_boundaries%abtype == mask_absorbing) then
        if (states_are_real(st)) then
          call dvmask(gr, hm, st)
        else
          call zvmask(gr, hm, st)
        end if
      end if
 
      !Photoelectron stuff
      if (td%pesv%calc_spm .or. td%pesv%calc_mask .or. td%pesv%calc_flux) then
        call pes_calc(td%pesv, namespace, space, gr, st, td%dt, iter, gr%der, hm%kpoints, ext_partners, stopping)
      end if
 
      call td_write_iter(td%write_handler, namespace, space, outp, gr, st, hm, ions, ext_partners, &
        hm%kick, ks, td%dt, iter, mc, td%recalculate_gs)
 
      ! write down data
      call td_check_point(td, namespace, mc, gr, ions, st, ks, hm, ext_partners, outp, space, &
        iter, scsteps, etime, stopping, from_scratch)
 
      ! check if debug mode should be enabled or disabled on the fly
      call io_debug_on_the_fly(namespace)
 
      call profiling_out("TIME_STEP")
      if (stopping) exit
 
    end do propagation
 
    pop_sub(td_run)
  end subroutine td_run
 
  subroutine td_print_header(namespace)
    type(namespace_t), intent(in)    :: namespace
 
    push_sub(td_print_header)
 
    write(message(1), '(a7,1x,a14,a14,a10,a17)') 'Iter ', 'Time ', 'Energy ', 'SC Steps', 'Elapsed Time '
 
    call messages_info(1, namespace=namespace)
    call messages_print_with_emphasis(namespace=namespace)
 
    pop_sub(td_print_header)
  end subroutine td_print_header
 
  ! ---------------------------------------------------------
  subroutine td_check_point(td, namespace, mc, gr, ions, st, ks, hm, ext_partners, outp, space, &
    iter, scsteps, etime, stopping, from_scratch)
    type(td_t),               intent(inout) :: td
    type(namespace_t),        intent(in)    :: namespace
    type(multicomm_t),        intent(in)    :: mc
    type(grid_t),             intent(inout) :: gr
    type(ions_t),             intent(inout) :: ions
    type(states_elec_t),      intent(inout) :: st
    type(v_ks_t),             intent(inout) :: ks
    type(hamiltonian_elec_t), intent(inout) :: hm
    type(partner_list_t),     intent(in)    :: ext_partners
    type(output_t),           intent(in)    :: outp
    type(electron_space_t),   intent(in)    :: space
    integer,                  intent(in)    :: iter
    integer,                  intent(in)    :: scsteps
    real(real64),             intent(inout) :: etime
    logical,                  intent(in)    :: stopping
    logical,                  intent(inout) :: from_scratch
 
    integer :: ierr
 
    push_sub(td_check_point)
 
    call td_print_message(td, namespace, ions, hm, iter, scsteps, etime)
 
    if (outp%anything_now(iter)) then ! output
      call td_write_output(namespace, space, gr, st, hm, ks, outp, ions, ext_partners, iter, td%dt)
    end if
 
    if (mod(iter, outp%restart_write_interval) == 0 .or. iter == td%max_iter .or. stopping) then ! restart
      !if (iter == td%max_iter) outp%iter = ii - 1
      call td_write_data(td%write_handler)
      call td_dump(td, namespace, space, gr, st, hm, ks, ext_partners, iter, ierr)
      if (ierr /= 0) then
        message(1) = "Unable to write time-dependent restart information."
        call messages_warning(1, namespace=namespace)
      end if
 
      call pes_output(td%pesv, namespace, space, gr, st, iter, outp, td%dt, ions)
 
      if ((td%ions_dyn%is_active()) .and. td%recalculate_gs) then
        call messages_print_with_emphasis(msg='Recalculating the ground state.', namespace=namespace)
        from_scratch = .false.
        call states_elec_deallocate_wfns(st)
        call electrons_ground_state_run(namespace, mc, gr, ions, ext_partners, st, ks, hm, outp, space, from_scratch)
        call states_elec_allocate_wfns(st, gr, packed=.true.)
        call td_load(td%restart_load, namespace, space, gr, st, hm, ext_partners, td, ks, ierr)
        if (ierr /= 0) then
          message(1) = "Unable to load TD states."
          call messages_fatal(1, namespace=namespace)
        end if
        call density_calc(st, gr, st%rho)
        call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners,  &
          calc_eigenval=.true., time = iter*td%dt, calc_energy=.true.)
        call forces_calculate(gr, namespace, ions, hm, ext_partners, st, ks, t = iter*td%dt, dt = td%dt)
        assert(.not. td%ions_dyn%cell_relax())
        call messages_print_with_emphasis(msg="Time-dependent simulation proceeds", namespace=namespace)
        call td_print_header(namespace)
      end if
    end if
 
    pop_sub(td_check_point)
  end subroutine td_check_point
 
  ! ---------------------------------------------------------
  subroutine td_print_message(td, namespace, ions, hm, iter, scsteps, etime)
    type(td_t),               intent(inout) :: td
    type(namespace_t),        intent(in)    :: namespace
    type(ions_t),             intent(inout) :: ions
    type(hamiltonian_elec_t), intent(inout) :: hm
    integer,                  intent(in)    :: iter
    integer,                  intent(in)    :: scsteps
    real(real64),             intent(inout) :: etime
 
    push_sub(td_print_message)
 
    write(message(1), '(i7,1x,2f14.6,i10,f14.3)') iter, units_from_atomic(units_out%time, iter*td%dt), &
      units_from_atomic(units_out%energy, hm%energy%total + ions%kinetic_energy), scsteps, &
      loct_clock() - etime
    call messages_info(1, namespace=namespace)
    call td_update_elapsed_time(etime)
 
    pop_sub(td_print_message)
  end subroutine td_print_message
 
  ! ---------------------------------------------------------
  subroutine td_update_elapsed_time(etime)
    real(real64), intent(inout) :: etime
 
    push_sub(td_update_elapsed_time)
 
    etime = loct_clock()
 
    pop_sub(td_update_elapsed_time)
  end subroutine td_update_elapsed_time
 
  ! ---------------------------------------------------------
  subroutine td_init_with_wavefunctions(td, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, outp, from_scratch)
    type(td_t),                  intent(inout) :: td
    type(namespace_t),           intent(in)    :: namespace
    type(electron_space_t),      intent(in)    :: space
    type(multicomm_t),           intent(in)    :: mc
    type(grid_t),                intent(inout) :: gr
    type(ions_t),                intent(inout) :: ions
    type(partner_list_t),        intent(in)    :: ext_partners
    type(states_elec_t), target, intent(inout) :: st
    type(v_ks_t),                intent(inout) :: ks
    type(hamiltonian_elec_t),    intent(inout) :: hm
    type(output_t),              intent(inout) :: outp
    logical,                     intent(in)    :: from_scratch
 
    integer :: ierr
    real(real64) :: x
    real(real64) :: ndinitial(space%dim)
    logical :: freeze_hxc, freeze_occ, freeze_u
    type(restart_t) :: restart, restart_frozen
    type(gauge_field_t), pointer :: gfield
    type(lasers_t),      pointer :: lasers
    push_sub(td_init_with_wavefunctions)
 
    !We activate the sprial BC only after the kick,
    !to be sure that the first iteration corresponds to the ground state
    if (gr%der%boundaries%spiralBC) then
      if ((td%iter-1)*td%dt > hm%kick%time) then
        gr%der%boundaries%spiral = .true.
      end if
      hm%vnl%spin => st%spin
      hm%phase%spin => st%spin
      !We fill st%spin. In case of restart, we read it in td_load
      if (from_scratch) call states_elec_fermi(st, namespace, gr)
    end if
 
    if (from_scratch) then
      ! Initialize the occupation matrices and U for DFT+U
      ! This must be called before parsing TDFreezeOccupations and TDFreezeU
      ! in order that the code does properly the initialization.
      call lda_u_update_occ_matrices(hm%lda_u, namespace, gr, st, hm%hm_base, hm%phase, hm%energy)
    end if
 
    if (td%freeze_orbitals > 0) then
      if (from_scratch) then
        ! In this case, we first freeze the orbitals, then calculate the Hxc potential.
        call states_elec_freeze_orbitals(st, namespace, space, gr, mc, hm%kpoints, &
          td%freeze_orbitals, family_is_mgga(ks%xc_family))
      else
        call restart_init(restart, namespace, restart_td, restart_type_load, mc, ierr, mesh=gr)
        if (ierr == 0) then
          call td_load_frozen(namespace, restart, space, gr, st, hm, ierr)
        end if
        if (ierr /= 0) then
          td%iter = 0
          message(1) = "Unable to read frozen restart information."
          call messages_fatal(1, namespace=namespace)
        end if
        call restart_end(restart)
      end if
      write(message(1),'(a,i4,a,i4,a)') 'Info: The lowest', td%freeze_orbitals, &
        ' orbitals have been frozen.', st%nst, ' will be propagated.'
      call messages_info(1, namespace=namespace)
      call states_elec_freeze_adjust_qtot(st)
      call density_calc(st, gr, st%rho)
      call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_eigenval=.true., time = td%iter*td%dt)
    else if (td%freeze_orbitals < 0) then
      ! This means SAE approximation. We calculate the Hxc first, then freeze all
      ! orbitals minus one.
      write(message(1),'(a)') 'Info: The single-active-electron approximation will be used.'
      call messages_info(1, namespace=namespace)
      call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_eigenval=.true., time = td%iter*td%dt)
      if (from_scratch) then
        call states_elec_freeze_orbitals(st, namespace, space, gr, mc, hm%kpoints, st%nst-1, family_is_mgga(ks%xc_family))
      else
        call messages_not_implemented("TDFreezeOrbials < 0 with FromScratch=no", namespace=namespace)
      end if
      call states_elec_freeze_adjust_qtot(st)
      call v_ks_freeze_hxc(ks)
      call density_calc(st, gr, st%rho)
    else
      ! Normal run.
      call density_calc(st, gr, st%rho)
      call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_eigenval=.true., time = td%iter*td%dt)
    end if
 
    !%Variable TDFreezeHXC
    !%Type logical
    !%Default no
    !%Section Time-Dependent
    !%Description
    !% The electrons are evolved as independent particles feeling the Hartree and
    !% exchange-correlation potentials from the ground-state electronic configuration.
    !%End
    call parse_variable(namespace, 'TDFreezeHXC', .false., freeze_hxc)
    if (freeze_hxc) then
      write(message(1),'(a)') 'Info: Freezing Hartree and exchange-correlation potentials.'
      call messages_info(1, namespace=namespace)
 
      if (.not. from_scratch) then
 
        call restart_init(restart_frozen, namespace, restart_gs, restart_type_load, mc, ierr, mesh=gr, exact=.true.)
        call states_elec_load(restart_frozen, namespace, space, st, gr, hm%kpoints, ierr, label = ": gs")
        call states_elec_transform(st, namespace, space, restart_frozen, gr, hm%kpoints)
        call restart_end(restart_frozen)
 
        call density_calc(st, gr, st%rho)
        call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_eigenval=.true., time = td%iter*td%dt)
 
        call restart_init(restart_frozen, namespace, restart_td, restart_type_load, mc, ierr, mesh=gr)
        call states_elec_load(restart_frozen, namespace, space, st, gr, hm%kpoints, ierr, iter=td%iter, label = ": td")
        call restart_end(restart_frozen)
        call hm%ks_pot%run_zero_iter(td%tr%vks_old)
 
      end if
 
      call v_ks_freeze_hxc(ks)
 
    end if
 
    x = minval(st%eigenval(st%st_start, :))
    if (st%parallel_in_states) then
      call st%mpi_grp%bcast(x, 1, mpi_double_precision, 0)
    end if
    call hm%update_span(gr%spacing(1:space%dim), x, namespace)
    ! initialize Fermi energy
    call states_elec_fermi(st, namespace, gr, compute_spin = .not. gr%der%boundaries%spiralBC)
    call energy_calc_total(namespace, space, hm, gr, st, ext_partners)
 
    !%Variable TDFreezeDFTUOccupations
    !%Type logical
    !%Default no
    !%Section Time-Dependent
    !%Description
    !% The occupation matrices than enters in the DFT+U potential
    !% are not evolved during the time evolution.
    !%End
    call parse_variable(namespace, 'TDFreezeDFTUOccupations', .false., freeze_occ)
    if (freeze_occ) then
      write(message(1),'(a)') 'Info: Freezing DFT+U occupation matrices that enters in the DFT+U potential.'
      call messages_info(1, namespace=namespace)
      call lda_u_freeze_occ(hm%lda_u)
 
      !In this case we should reload GS wavefunctions
      if (hm%lda_u_level /= dft_u_none .and. .not. from_scratch) then
        call restart_init(restart_frozen, namespace, restart_gs, restart_type_load, mc, ierr, mesh=gr)
        call lda_u_load(restart_frozen, hm%lda_u, st, hm%energy%dft_u, ierr, occ_only = .true.)
        call restart_end(restart_frozen)
      end if
    end if
 
    !%Variable TDFreezeU
    !%Type logical
    !%Default no
    !%Section Time-Dependent
    !%Description
    !% The effective U of DFT+U is not evolved during the time evolution.
    !%End
    call parse_variable(namespace, 'TDFreezeU', .false., freeze_u)
    if (freeze_u) then
      write(message(1),'(a)') 'Info: Freezing the effective U of DFT+U.'
      call messages_info(1, namespace=namespace)
      call lda_u_freeze_u(hm%lda_u)
 
      !In this case we should reload GS wavefunctions
      if (hm%lda_u_level == dft_u_acbn0 .and. .not. from_scratch) then
        call restart_init(restart_frozen, namespace, restart_gs, restart_type_load, mc, ierr, mesh=gr)
        call lda_u_load(restart_frozen, hm%lda_u, st, hm%energy%dft_u, ierr, u_only = .true.)
        call restart_end(restart_frozen)
        write(message(1),'(a)') 'Loaded GS effective U of DFT+U'
        call messages_info(1, namespace=namespace)
        call lda_u_write_u(hm%lda_u, namespace=namespace)
        call lda_u_write_v(hm%lda_u, namespace=namespace)
      end if
    end if
 
    ! This needs to be called before the calculation of the forces,
    ! as we need to test of we output the forces or not
    call td_write_init(td%write_handler, namespace, space, outp, gr, st, hm, ions, ext_partners, &
      ks, td%ions_dyn%is_active(), &
      list_has_gauge_field(ext_partners), hm%kick, td%iter, td%max_iter, td%dt, mc)
 
    ! Resets the nondipole integration after laser-file has been written.
    lasers => list_get_lasers(ext_partners)
    if(associated(lasers)) then
      if (lasers_with_nondipole_field(lasers)) then
        ndinitial(1:space%dim)=m_zero
        call lasers_set_nondipole_parameters(lasers,ndinitial,m_zero)
      end if
    end if
    nullify(lasers)
 
    call td_init_ions_and_forces(td, namespace, space, gr, ions, ext_partners, st, ks, hm, outp)
 
    if (td%scissor > m_epsilon) then
      call scissor_init(hm%scissor, namespace, space, st, gr, hm%d, hm%kpoints, hm%phase, td%scissor, mc)
    end if
 
    if (td%iter == 0) call td_run_zero_iter(td, namespace, space, gr, ions, st, ks, hm, ext_partners, outp, mc)
 
    gfield => list_get_gauge_field(ext_partners)
    if(associated(gfield)) then
      if (gauge_field_is_propagated(gfield)) then
        if(ks%xc%kernel_lrc_alpha > m_epsilon) then
          call gauge_field_get_force(gfield, gr, st%d%spin_channels, st%current, ks%xc%kernel_lrc_alpha)
          call messages_experimental('TD-LRC kernel')
        else
          call gauge_field_get_force(gfield, gr, st%d%spin_channels, st%current)
        endif
      endif
    end if
 
    !call td_check_trotter(td, sys, h)
    td%iter = td%iter + 1
 
    call restart_init(td%restart_dump, namespace, restart_td, restart_type_dump, mc, ierr, mesh=gr)
    if (td%ions_dyn%is_active() .and. td%recalculate_gs) then
      ! We will also use the TD restart directory as temporary storage during the time propagation
      call restart_init(td%restart_load, namespace, restart_td, restart_type_load, mc, ierr, mesh=gr)
    end if
 
    call messages_print_with_emphasis(msg="Time-Dependent Simulation", namespace=namespace)
    call td_print_header(namespace)
 
    if (td%pesv%calc_spm .or. td%pesv%calc_mask .and. from_scratch) then
      call pes_init_write(td%pesv,gr,st, namespace)
    end if
 
    if (st%pack_states .and. hm%apply_packed()) call st%pack()
 
    pop_sub(td_init_with_wavefunctions)
  end subroutine td_init_with_wavefunctions
 
  ! ---------------------------------------------------------
  subroutine td_init_ions_and_forces(td, namespace, space, gr, ions, ext_partners, st, ks, hm, outp)
    type(td_t),                  intent(inout) :: td
    type(namespace_t),           intent(in)    :: namespace
    type(electron_space_t),      intent(in)    :: space
    type(grid_t),                intent(inout) :: gr
    type(ions_t),                intent(inout) :: ions
    type(partner_list_t),        intent(in)    :: ext_partners
    type(states_elec_t), target, intent(inout) :: st
    type(v_ks_t),                intent(inout) :: ks
    type(hamiltonian_elec_t),    intent(inout) :: hm
    type(output_t),              intent(inout) :: outp
 
    push_sub(td_init_ions_and_forces)
 
    ! Calculate initial forces and kinetic energy
    if (td%ions_dyn%ions_move()) then
      if (td%iter > 0) then
        call td_read_coordinates(td, namespace, ions)
        if (ion_dynamics_drive_ions(td%ions_dyn)) then
          call ion_dynamics_propagate(td%ions_dyn, ions, td%iter*td%dt, td%dt, namespace)
        end if
        call hamiltonian_elec_epot_generate(hm, namespace, space, gr, ions, ext_partners, st, time = td%iter*td%dt)
        ! recompute potential because the ions have moved
        call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_eigenval=.true., time = td%iter*td%dt)
      end if
 
      call forces_calculate(gr, namespace, ions, hm, ext_partners, st, ks, t = td%iter*td%dt, dt = td%dt)
 
      call ions%update_kinetic_energy()
    else
      if (outp%what(option__output__forces) .or. td%write_handler%out(out_separate_forces)%write) then
        call forces_calculate(gr, namespace, ions, hm, ext_partners, st, ks, t = td%iter*td%dt, dt = td%dt)
      end if
    end if
 
    if (outp%what(option__output__stress) .or. td%ions_dyn%cell_relax()) then
      call stress_calculate(namespace, gr, hm, st, ions, ks, ext_partners)
      if (td%ions_dyn%cell_relax()) then
        call td%ions_dyn%update_stress(ions%space, st%stress_tensors%total, ions%latt%rlattice, ions%latt%rcell_volume)
      end if
    end if
 
    pop_sub(td_init_ions_and_forces)
  end subroutine td_init_ions_and_forces
 
  ! ---------------------------------------------------------
  subroutine td_load_restart_from_td(td, namespace, space, mc, gr, ext_partners, st, ks, hm, from_scratch)
    type(td_t),                  intent(inout) :: td
    type(namespace_t),           intent(in)    :: namespace
    class(space_t),              intent(in)    :: space
    type(multicomm_t),           intent(in)    :: mc
    type(grid_t),                intent(inout) :: gr
    type(partner_list_t),        intent(in)    :: ext_partners
    type(states_elec_t), target, intent(inout) :: st
    type(v_ks_t),                intent(inout) :: ks
    type(hamiltonian_elec_t),    intent(inout) :: hm
    logical,                     intent(inout) :: from_scratch
 
    integer :: ierr
    type(restart_t) :: restart
 
    push_sub(td_load_restart_from_td)
 
    !We redistribute the states before the restarting
    if (td%freeze_orbitals > 0) then
      call states_elec_freeze_redistribute_states(st, namespace, gr, mc, td%freeze_orbitals)
    end if
 
    call restart_init(restart, namespace, restart_td, restart_type_load, mc, ierr, mesh=gr)
    if (ierr == 0) then
      call td_load(restart, namespace, space, gr, st, hm, ext_partners, td, ks, ierr)
    end if
    call restart_end(restart)
    if (ierr /= 0) then
      from_scratch = .true.
      td%iter = 0
    end if
 
    pop_sub(td_load_restart_from_td)
  end subroutine td_load_restart_from_td
 
  ! ---------------------------------------------------------
  subroutine td_load_restart_from_gs(td, namespace, space, mc, gr, ext_partners, st, ks, hm)
    type(td_t),                  intent(inout) :: td
    type(namespace_t),           intent(in)    :: namespace
    class(space_t),              intent(in)    :: space
    type(multicomm_t),           intent(in)    :: mc
    type(grid_t),                intent(inout) :: gr
    type(partner_list_t),        intent(in)    :: ext_partners
    type(states_elec_t), target, intent(inout) :: st
    type(v_ks_t),                intent(inout) :: ks
    type(hamiltonian_elec_t),    intent(inout) :: hm
 
    integer :: ierr
    type(restart_t) :: restart
 
    push_sub(td_load_restart_from_gs)
 
    call restart_init(restart, namespace, restart_gs, restart_type_load, mc, ierr, mesh=gr, exact=.true.)
 
    if (.not. st%only_userdef_istates) then
      if (ierr == 0) then
        call states_elec_load(restart, namespace, space, st, gr, hm%kpoints, ierr, label = ": gs")
      end if
      if (ierr /= 0) then
        message(1) = 'Unable to read ground-state wavefunctions.'
        call messages_fatal(1, namespace=namespace)
      end if
    end if
 
    ! check if we should deploy user-defined wavefunctions.
    ! according to the settings in the input file the routine
    ! overwrites orbitals that were read from restart/gs
    if (parse_is_defined(namespace, 'UserDefinedStates')) then
      call states_elec_read_user_def_orbitals(gr, namespace, space, st)
    end if
 
    call states_elec_transform(st, namespace, space, restart, gr, hm%kpoints)
    call restart_end(restart)
 
    pop_sub(td_load_restart_from_gs)
  end subroutine td_load_restart_from_gs
 
  ! ---------------------------------------------------------
  subroutine td_run_zero_iter(td, namespace, space, gr, ions, st, ks, hm, ext_partners, outp, mc)
    type(td_t),               intent(inout) :: td
    type(namespace_t),        intent(in)    :: namespace
    type(electron_space_t),   intent(in)    :: space
    type(grid_t),             intent(inout) :: gr
    type(ions_t),             intent(inout) :: ions
    type(states_elec_t),      intent(inout) :: st
    type(v_ks_t),             intent(inout) :: ks
    type(hamiltonian_elec_t), intent(inout) :: hm
    type(partner_list_t),     intent(in)    :: ext_partners
    type(output_t),           intent(in)    :: outp
    type(multicomm_t),        intent(in)    :: mc
 
    push_sub(td_run_zero_iter)
 
    call td_write_iter(td%write_handler, namespace, space, outp, gr, st, hm, ions, ext_partners, &
      hm%kick, ks, td%dt, 0, mc, td%recalculate_gs)
 
    ! I apply the delta electric field *after* td_write_iter, otherwise the
    ! dipole matrix elements in write_proj are wrong
    if (abs(hm%kick%time)  <=  m_epsilon) then
      if (.not. hm%pcm%localf) then
        call kick_apply(space, gr, st, td%ions_dyn, ions, hm%kick, hm%psolver, hm%kpoints)
      else
        call kick_apply(space, gr, st, td%ions_dyn, ions, hm%kick, hm%psolver, hm%kpoints, pcm = hm%pcm)
      end if
      call td_write_kick(outp, namespace, space, gr, hm%kick, ions, 0)
 
      !We activate the sprial BC only after the kick
      if (gr%der%boundaries%spiralBC) then
        gr%der%boundaries%spiral = .true.
      end if
    end if
    call hm%ks_pot%run_zero_iter(td%tr%vks_old)
 
    if (any(outp%output_interval > 0)) then
      call td_write_data(td%write_handler)
      call td_write_output(namespace, space, gr, st, hm, ks, outp, ions, ext_partners, 0)
    end if
 
    pop_sub(td_run_zero_iter)
  end subroutine td_run_zero_iter
 
 
  ! ---------------------------------------------------------
  subroutine td_read_coordinates(td, namespace, ions)
    type(td_t),               intent(in)    :: td
    type(namespace_t),        intent(in)    :: namespace
    type(ions_t),             intent(inout) :: ions
 
    integer :: iatom, iter, iunit
 
    push_sub(td_read_coordinates)
 
    iunit = io_open('td.general/coordinates', namespace, action='read', status='old', die=.false.)
    if (iunit == -1) then
      message(1) = "Could not open file '"//trim(io_workpath('td.general/coordinates', namespace))//"'."
      message(2) = "Starting simulation from initial geometry."
      call messages_warning(2, namespace=namespace)
      pop_sub(td_read_coordinates)
      return
    end if
 
    call io_skip_header(iunit)
    do iter = 0, td%iter - 1
      read(iunit, *) ! skip previous iterations... sorry, but no portable seek in Fortran
    end do
    read(iunit, '(32x)', advance='no') ! skip the time index.
 
    do iatom = 1, ions%natoms
      read(iunit, '(3es24.16)', advance='no') ions%pos(:, iatom)
      ions%pos(:, iatom) = units_to_atomic(units_out%length, ions%pos(:, iatom))
    end do
    do iatom = 1, ions%natoms
      read(iunit, '(3es24.16)', advance='no') ions%vel(:, iatom)
      ions%vel(:, iatom) = units_to_atomic(units_out%velocity, ions%vel(:, iatom))
    end do
    do iatom = 1, ions%natoms
      read(iunit, '(3es24.16)', advance='no') ions%tot_force(:, iatom)
      ions%tot_force(:, iatom) = units_to_atomic(units_out%force, ions%tot_force(:, iatom))
    end do
 
    call io_close(iunit)
 
    pop_sub(td_read_coordinates)
  end subroutine td_read_coordinates
 
  ! ---------------------------------------------------------
  subroutine td_dump(td, namespace, space, gr, st, hm, ks, ext_partners, iter, ierr)
    type(td_t),               intent(in)  :: td
    type(namespace_t),        intent(in)  :: namespace
    class(space_t),           intent(in)  :: space
    type(grid_t),             intent(in)  :: gr
    type(states_elec_t),      intent(in)  :: st
    type(hamiltonian_elec_t), intent(in)  :: hm
    type(v_ks_t),             intent(in)  :: ks
    type(partner_list_t),     intent(in)  :: ext_partners
    integer,                  intent(in)  :: iter
    integer,                  intent(out) :: ierr
 
    type(gauge_field_t), pointer :: gfield
    integer :: err, err2
 
    push_sub(td_dump)
 
    ierr = 0
 
    if (restart_skip(td%restart_dump)) then
      pop_sub(td_dump)
      return
    end if
 
    message(1) = "Debug: Writing td restart."
    call messages_info(1, namespace=namespace, debug_only=.true.)
 
    ! first write resume file
    call states_elec_dump(td%restart_dump, space, st, gr, hm%kpoints, err, iter=iter)
    if (err /= 0) ierr = ierr + 1
 
    call states_elec_dump_rho(td%restart_dump, space, st, gr, ierr, iter=iter)
    if (err /= 0) ierr = ierr + 1
 
    if (hm%lda_u_level /= dft_u_none) then
      call lda_u_dump(td%restart_dump, namespace, hm%lda_u, st, gr, ierr)
      if (err /= 0) ierr = ierr + 1
    end if
 
    call potential_interpolation_dump(td%tr%vks_old, space, td%restart_dump, gr, st%d%nspin, err2)
    if (err2 /= 0) ierr = ierr + 2
 
    call pes_dump(td%pesv, namespace, td%restart_dump, st, gr, err)
    if (err /= 0) ierr = ierr + 4
 
    ! Gauge field restart
    gfield => list_get_gauge_field(ext_partners)
    if(associated(gfield)) then
      call gauge_field_dump(td%restart_dump, gfield, ierr)
    end if
 
    if (gr%der%boundaries%spiralBC) then
      call states_elec_dump_spin(td%restart_dump, st, err)
      if (err /= 0) ierr = ierr + 8
    end if
 
    if (ks%has_photons) then
      call mf_photons_dump(td%restart_dump, ks%pt_mx, gr, td%dt, ks%pt, err)
      if (err /= 0) ierr = ierr + 16
    end if
 
    if (ks%xc_photon /= 0) then
      ! photon-free mean field
      call ks%xc_photons%mf_dump(td%restart_dump, err)
      if (err /= 0) ierr = ierr + 32
    end if
 
    if (allocated(st%frozen_rho)) then
      call states_elec_dump_frozen(td%restart_dump, space, st, gr, ierr)
    end if
    if (err /= 0) ierr = ierr + 64
 
    if (td%ions_dyn%ions_move() .or. td%ions_dyn%cell_relax()) then
      call ion_dynamics_dump(td%ions_dyn, td%restart_dump, err)
    end if
    if (err /= 0) ierr = ierr + 128
 
    message(1) = "Debug: Writing td restart done."
    call messages_info(1, namespace=namespace, debug_only=.true.)
 
    pop_sub(td_dump)
  end subroutine td_dump
 
  ! ---------------------------------------------------------
  subroutine td_load(restart, namespace, space, gr, st, hm, ext_partners, td, ks, ierr)
    type(restart_t),          intent(in)    :: restart
    type(namespace_t),        intent(in)    :: namespace
    class(space_t),           intent(in)    :: space
    type(grid_t),             intent(in)    :: gr
    type(states_elec_t),      intent(inout) :: st
    type(hamiltonian_elec_t), intent(inout) :: hm
    type(partner_list_t),     intent(in)    :: ext_partners
    type(td_t),               intent(inout) :: td
    type(v_ks_t),             intent(inout) :: ks
    integer,                  intent(out)   :: ierr
 
    integer :: err, err2
    type(gauge_field_t), pointer :: gfield
    push_sub(td_load)
 
    ierr = 0
 
    if (restart_skip(restart)) then
      ierr = -1
      pop_sub(td_load)
      return
    end if
 
    message(1) = "Debug: Reading td restart."
    call messages_info(1, namespace=namespace, debug_only=.true.)
 
    ! Read states
    call states_elec_load(restart, namespace, space, st, gr, hm%kpoints, err, iter=td%iter, label = ": td")
    if (err /= 0) then
      ierr = ierr + 1
    end if
 
    ! read potential from previous interactions
    call potential_interpolation_load(td%tr%vks_old, namespace, space, restart, gr, st%d%nspin, err2)
    if (err2 /= 0) ierr = ierr + 2
 
    if (hm%lda_u_level /= dft_u_none) then
      call lda_u_load(restart, hm%lda_u, st, hm%energy%dft_u, ierr)
      if (err /= 0) ierr = ierr + 1
    end if
 
 
    ! read PES restart
    if (td%pesv%calc_spm .or. td%pesv%calc_mask .or. td%pesv%calc_flux) then
      call pes_load(td%pesv, namespace, restart, st, err)
      if (err /= 0) ierr = ierr + 4
    end if
 
    ! Gauge field restart
    gfield => list_get_gauge_field(ext_partners)
    if (associated(gfield)) then
      call gauge_field_load(restart, gfield, err)
      if (err /= 0) then
        ierr = ierr + 8
      else
        call hm%update(gr, namespace, space, ext_partners, time = td%dt*td%iter)
      end if
    end if
 
    ! add photon restart
    if (ks%has_photons) then
      call mf_photons_load(restart, ks%pt_mx, gr, err)
    end if
    if (err /= 0) ierr = ierr + 16
 
    if (ks%xc_photon /= 0) then
      call ks%xc_photons%mf_load(restart, space, err)
    end if
    if (err /= 0) ierr = ierr + 32
 
    if (gr%der%boundaries%spiralBC) then
      call states_elec_load_spin(restart, st, err)
      !To ensure back compatibility, if the file is not present, we use the
      !current states to get the spins
      if (err /= 0) call states_elec_fermi(st, namespace, gr)
    end if
 
    if (td%ions_dyn%is_active()) then
      call ion_dynamics_load(td%ions_dyn, restart, err)
    end if
    if (err /= 0) ierr = ierr + 64
 
    message(1) = "Debug: Reading td restart done."
    call messages_info(1, namespace=namespace, debug_only=.true.)
 
    pop_sub(td_load)
  end subroutine td_load
 
  ! ---------------------------------------------------------
  subroutine td_load_frozen(namespace, restart, space, mesh, st, hm, ierr)
    type(namespace_t),        intent(in)    :: namespace
    type(restart_t),          intent(in)    :: restart
    class(space_t),           intent(in)    :: space
    class(mesh_t),            intent(in)    :: mesh
    type(states_elec_t),      intent(inout) :: st
    type(hamiltonian_elec_t), intent(inout) :: hm
    integer,                  intent(out)   :: ierr
 
    push_sub(td_load_frozen)
 
    ierr = 0
 
    if (restart_skip(restart)) then
      ierr = -1
      pop_sub(td_load_frozen)
      return
    end if
 
    message(1) = "Debug: Reading td frozen restart."
    call messages_info(1, namespace=namespace, debug_only=.true.)
 
    safe_allocate(st%frozen_rho(1:mesh%np, 1:st%d%nspin))
    if (family_is_mgga(hm%xc%family)) then
      safe_allocate(st%frozen_tau(1:mesh%np, 1:st%d%nspin))
      safe_allocate(st%frozen_gdens(1:mesh%np, 1:space%dim, 1:st%d%nspin))
      safe_allocate(st%frozen_ldens(1:mesh%np, 1:st%d%nspin))
    end if
 
    call states_elec_load_frozen(restart, space, st, mesh, ierr)
 
    message(1) = "Debug: Reading td frozen restart done."
    call messages_info(1, namespace=namespace, debug_only=.true.)
 
    pop_sub(td_load_frozen)
  end subroutine td_load_frozen
 
  ! ---------------------------------------------------------
  logical function td_get_from_scratch(td)
    type(td_t), intent(in) :: td
 
    push_sub(td_get_from_scratch)
 
    td_get_from_scratch = td%from_scratch
 
    pop_sub(td_get_from_scratch)
  end function td_get_from_scratch
 
  ! ---------------------------------------------------------
  subroutine td_set_from_scratch(td, from_scratch)
    type(td_t), intent(inout) :: td
    logical,    intent(in)    :: from_scratch
 
    push_sub(td_set_from_scratch)
 
    td%from_scratch = from_scratch
 
    pop_sub(td_set_from_scratch)
  end subroutine td_set_from_scratch
end module td_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: