propagator_elec.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 propagator_elec_oct_m
  use debug_oct_m
  use electron_space_oct_m
  use energy_calc_oct_m
  use exponential_oct_m
  use ext_partner_list_oct_m
  use forces_oct_m
  use gauge_field_oct_m
  use grid_oct_m
  use global_oct_m
  use hamiltonian_elec_oct_m
  use interaction_partner_oct_m
  use ion_dynamics_oct_m
  use ions_oct_m
  use, intrinsic :: iso_fortran_env
  use lasers_oct_m
  use lda_u_oct_m
  use ks_potential_oct_m
  use parser_oct_m
  use mesh_function_oct_m
  use messages_oct_m
  use multicomm_oct_m
  use namespace_oct_m
  use opt_control_state_oct_m
  use output_low_oct_m
  use potential_interpolation_oct_m
  use profiling_oct_m
  use propagation_ops_elec_oct_m
  use propagator_base_oct_m
  use propagator_cn_oct_m
  use propagator_etrs_oct_m
  use propagator_expmid_oct_m
  use propagator_magnus_oct_m
  use propagator_qoct_oct_m
  use propagator_rk_oct_m
  use propagator_verlet_oct_m
  use scf_oct_m
  use sparskit_oct_m
  use space_oct_m
  use states_elec_oct_m
  use stress_oct_m
  use td_write_oct_m
  use v_ks_oct_m
  use varinfo_oct_m
  use xc_oct_m
 
  implicit none
 
  private
  public ::                         &
    propagator_elec_init,                &
    propagator_elec_end,                 &
    propagator_elec_copy,                &
    propagator_elec_dt,                  &
    propagator_elec_set_scf_prop,        &
    propagator_elec_remove_scf_prop,     &
    propagator_elec_ions_are_propagated, &
    propagator_elec_dt_bo
 
contains
 
  ! ---------------------------------------------------------
  subroutine propagator_elec_copy(tro, tri)
    type(propagator_base_t), intent(inout) :: tro
    type(propagator_base_t), intent(in)    :: tri
 
    push_sub(propagator_elec_copy)
 
    tro%method = tri%method
 
    select case (tro%method)
    case (prop_crank_nicolson_sparskit)
      tro%tdsk_size = tri%tdsk_size
      call sparskit_solver_copy(tro%tdsk, tri%tdsk)
 
    case (prop_runge_kutta4)
      tro%tdsk_size = tri%tdsk_size
      call sparskit_solver_copy(tro%tdsk, tri%tdsk)
 
    case (prop_runge_kutta2)
      tro%tdsk_size = tri%tdsk_size
      call sparskit_solver_copy(tro%tdsk, tri%tdsk)
 
    end select
 
    call potential_interpolation_copy(tro%vks_old, tri%vks_old)
 
    call exponential_copy(tro%te, tri%te)
    tro%scf_propagation_steps = tri%scf_propagation_steps
 
    tro%scf_threshold = tri%scf_threshold
    pop_sub(propagator_elec_copy)
  end subroutine propagator_elec_copy
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine propagator_elec_init(gr, namespace, st, tr, ks_pot, have_fields, family_is_mgga_with_exc, relax_cell)
    type(grid_t),            intent(in)    :: gr
    type(namespace_t),       intent(in)    :: namespace
    type(states_elec_t),     intent(in)    :: st
    type(propagator_base_t), intent(inout) :: tr
    type(ks_potential_t),    intent(in)    :: ks_pot
    logical,                 intent(in)    :: have_fields
    logical,                 intent(in)    :: family_is_mgga_with_exc
    logical,                 intent(in)    :: relax_cell
 
    push_sub(propagator_elec_init)
 
    !%Variable TDPropagator
    !%Type integer
    !%Default etrs
    !%Section Time-Dependent::Propagation
    !%Description
    !% This variable determines which algorithm will be used to approximate
    !% the evolution operator <math>U(t+\delta t, t)</math>. That is, given
    !% <math>\psi(\tau)</math> and <math>H(\tau)</math> for <math>\tau \le t</math>,
    !% calculate <math>t+\delta t</math>. Note that in general the Hamiltonian
    !% is not known at times in the interior of the interval <math>[t,t+\delta t]</math>.
    !% This is due to the self-consistent nature of the time-dependent Kohn-Sham problem:
    !% the Hamiltonian at a given time <math>\tau</math> is built from the
    !% "solution" wavefunctions at that time.
    !%
    !% Some methods, however, do require the knowledge of the Hamiltonian at some
    !% point of the interval <math>[t,t+\delta t]</math>. This problem is solved by making
    !% use of extrapolation: given a number <math>l</math> of time steps previous to time
    !% <math>t</math>, this information is used to build the Hamiltonian at arbitrary times
    !% within <math>[t,t+\delta t]</math>. To be fully precise, one should then proceed
    !% <i>self-consistently</i>: the obtained Hamiltonian at time <math>t+\delta t</math>
    !% may then be used to interpolate the Hamiltonian, and repeat the evolution
    !% algorithm with this new information. Whenever iterating the procedure does
    !% not change the solution wavefunctions, the cycle is stopped. In practice,
    !% in <tt>Octopus</tt> we perform a second-order extrapolation without a
    !% self-consistency check, except for the first two iterations, where obviously
    !% the extrapolation is not reliable.
    !%
    !% The proliferation of methods is certainly excessive. The reason for it is that
    !% the propagation algorithm is currently a topic of active development. We
    !% hope that in the future the optimal schemes are clearly identified. In the
    !% mean time, if you do not feel like testing, use the default choices and
    !% make sure the time step is small enough.
    !%
    !% More details about the propagation methods could be found in
    !% A. Castro, M.A.L. Marques, and A. Rubio, J. Chem. Phys 121 3425-3433 (2004)
    !% and
    !% A. Pueyo et al., J. Chem. Theory Comput. 14, 6, 3040 (2018)
    !%
    !%Option etrs 2
    !% The idea is to make use of time-reversal symmetry from the beginning:
    !%
    !% <math>
    !%   \exp \left(-i\delta t H_{n} / 2 \right)\psi_n = \exp \left(i\delta t H_{n+1} / 2 \right)\psi_{n+1},
    !% </math>
    !%
    !% and then invert to obtain:
    !%
    !% <math>
    !%   \psi_{n+1} = \exp \left(-i\delta t H_{n+1} / 2 \right) \exp \left(-i\delta t H_{n} / 2 \right)\psi_{n}.
    !% </math>
    !%
    !% But we need to know <math>H_{n+1}</math>, which can only be known exactly through the solution
    !% <math>\psi_{n+1}</math>. What we do is to estimate it by performing a single exponential:
    !% <math>\psi\^{\*}\_{n+1}=\exp \left( -i\delta t H_{n} \right) \psi_n</math>, and then
    !% <math>H_{n+1} = H[\psi\^{\*}_{n+1}]</math>. Thus no extrapolation is performed in this case.
    !%Option aetrs 3
    !% Approximated Enforced Time-Reversal Symmetry (AETRS).
    !% A modification of previous method to make it faster.
    !% It is based on extrapolation of the time-dependent potentials. It is faster
    !% by about 40%.
    !% The only difference is the procedure to estimate <math>H_{n+1}</math>: in this case
    !% it is extrapolated via a second-order polynomial by making use of the
    !% Hamiltonian at time <math>t-2\delta t</math>, <math>t-\delta t</math> and <math>t</math>.
    !%Option caetrs 12
    !% (experimental) Corrected Approximated Enforced Time-Reversal
    !% Symmetry (AETRS), this is the previous propagator but including
    !% a correction step to the exponential.
    !%Option exp_mid 4
    !% Exponential Midpoint Rule (EM).
    !% This is maybe the simplest method, but it is very well grounded theoretically:
    !% it is unitary (if the exponential is performed correctly) and preserves
    !% time-reversal symmetry (if the self-consistency problem is dealt with correctly).
    !% It is defined as:
    !% <math>
    !%   U_{\rm EM}(t+\delta t, t) = \exp \left( -i\delta t H_{t+\delta t/2}\right)\,.
    !% </math>
    !%Option crank_nicholson 5
    !%Option crank_nicolson 5
    !% Classical Crank-Nicolson propagator.
    !% <math>
    !%  (1 + i\delta t H_{n+1/2} / 2) \psi_{n+1} = (1 - i\delta t H_{n+1/2} / 2) \psi_{n}
    !% </math>
    !%Option crank_nicholson_sparskit 6
    !%Option crank_nicolson_sparskit 6
    !% Classical Crank-Nicolson propagator. Requires the SPARSKIT library.
    !% <math>
    !%  (1 + i\delta t H_{n+1/2} / 2) \psi_{n+1} = (1 - i\delta t H_{n+1/2} / 2) \psi_{n}
    !% </math>
    !%Option magnus 7
    !% Magnus Expansion (M4).
    !% This is the most sophisticated approach. It is a fourth-order scheme (a feature
    !% which it shares with the ST scheme; the other schemes are in principle second-order).
    !% It is tailored for making use of very large time steps, or equivalently,
    !% dealing with problem with very high-frequency time-dependence.
    !% It is still in a experimental state; we are not yet sure of when it is
    !% advantageous.
    !%Option qoct_tddft_propagator 10
    !% WARNING: EXPERIMENTAL
    !%Option runge_kutta4 13
    !% WARNING: EXPERIMENTAL. Implicit Gauss-Legendre 4th order Runge-Kutta.
    !%Option runge_kutta2 14
    !% WARNING: EXPERIMENTAL. Implicit 2nd order Runge-Kutta (trapezoidal rule).
    !% Similar, but not identical, to Crank-Nicolson method.
    !%Option expl_runge_kutta4 15
    !% WARNING: EXPERIMENTAL. Explicit RK4 method.
    !%Option cfmagnus4 16
    !% Commutator-free magnus expansion. This method has been developed for linear
    !% systems and involves products of exponentials to avoid commutators. To apply it
    !% to our non-linear equations, we extrapolate the potential from previous timesteps
    !% to the times needed.The fourth-ordermethod implemented here follows equation 54
    !% in Gómez Pueyo, A., Marques, M. A. L., Rubio, A. & Castro, A. Propagators for the
    !% Time-Dependent Kohn–Sham Equations: Multistep, Runge–Kutta, Exponential Runge–Kutta,
    !% and Commutator Free Magnus Methods. J. Chem. Theory Comput. 14, 3040–3052 (2018)
    !% which corresponds to equation 43 in Blanes, S. & Moan, P. C. Fourth- and sixth-order
    !% commutator-free Magnus integrators for linear and non-linear dynamical systems.
    !% Applied Numerical Mathematics 56, 1519–1537 (2006):
    !% <math>
    !% \psi(t+\Delta t) = \exp(-i \Delta t ( \alpha_1 H(t_1) + \alpha_2 H(t_2) ) )
    !%                    \exp(-i \Delta t ( \alpha_2 H(t_1) + \alpha_1 H(t_2) ) ) \psi(t),
    !% </math>
    !% with $\alpha_{1,2} = \frac{3\mp 2\sqrt{3}}{12}$ and $t_i = t_0 + c_i \Delta t$,
    !% $c_{1,2} = \frac{1}{2} \mp \frac{\sqrt{3}}{6}$.
    !% Still in experimental state.
    !%Option exp_mid_predictor 17
    !% Exponential Midpoint Rule (EM).
    !% This is maybe the simplest method, but it is very well grounded theoretically:
    !% it is unitary (if the exponential is performed correctly) and preserves
    !% time-reversal symmetry (if the self-consistency problem is dealt with correctly).
    !% It is defined as:
    !% <math>
    !%   U_{\rm EM}(t+\delta t, t) = \exp \left( -i\delta t H_{t+\delta t/2}\right)\,.
    !% </math>
    !% As opposed to the exp_mid propagator, $H_{t+\delta t/2}$ is not obtained by
    !% extrapolation from previous timesteps, but by predicting by applying the exponential with
    !% $\delta t/2$ and computing the corresponding density and potential.
    !% It is equivalent to an explicit Magnus integration of second order, see equation 20 in
    !% Casas, F. & Iserles, A. Explicit Magnus expansions for nonlinear equations. J. Phys. A: Math. Gen. 39,
    !% 5445 (2006) [https://dx.doi.org/10.1088/0305-4470/39/19/S07]
    !%Option explicit_magnus3 18
    !% Third-order explicit Magnus expansion.
    !% This method follows equations 22 in
    !% Casas, F. & Iserles, A. Explicit Magnus expansions for nonlinear equations.
    !% J. Phys. A: Math. Gen. 39, 5445 (2006) [https://dx.doi.org/10.1088/0305-4470/39/19/S07]
    !% Applying eq. 22 to our notation for the Schrödinger equation yields:
    !% <math>
    !% \psi_0 = \psi(t) \\
    !% H_1 = H(t, \psi_0) \\
    !% \psi_2 = \exp(-i \frac{\Delta t}{2} H_1) \psi_0 \\
    !% H_2 = H(t+\frac{\Delta t}{2}, \psi_2) \\
    !% \psi_3 = \exp(-i \frac{\Delta t}{4} (H_1 + H_2)) \psi_0 \\
    !% H_3 = H(t+\frac{\Delta t}{2}, \psi_3) \\
    !% \psi_4 = \exp(-i \Delta t H_2) \psi_0 \\
    !% H_4 = H(t+\Delta t, \psi_4) \\
    !% \psi(t+\Delta t) = \exp(-i \Delta t (  \frac{1}{6}(H_1 + 4 H_3 + H_4)
    !%                  + \frac{i \Delta t}{12} ( [H_1 + H_2, H_3] + [H_2, H_4]))) \psi_0
    !% </math>
    !% Still in experimental state.
    !%End
    call messages_obsolete_variable(namespace, 'TDEvolutionMethod', 'TDPropagator')
 
    call parse_variable(namespace, 'TDPropagator', prop_etrs, tr%method)
    if (.not. varinfo_valid_option('TDPropagator', tr%method)) call messages_input_error(namespace, 'TDPropagator')
 
    select case (tr%method)
    case (prop_etrs)
    case (prop_aetrs)
    case (prop_caetrs)
      call messages_experimental("CAETRS propagator", namespace=namespace)
    case (prop_exponential_midpoint)
    case (prop_exponential_midpoint_predictor)
    case (prop_crank_nicolson)
      ! set up pointer for zmf_dotu_aux, zmf_nrm2_aux
      call mesh_init_mesh_aux(gr)
    case (prop_runge_kutta4)
      ! set up pointer for zmf_dotu_aux, zmf_nrm2_aux
      call mesh_init_mesh_aux(gr)
      sp_distdot_mode = 3
      tr%tdsk_size = 2 * st%d%dim * gr%np * (st%st_end - st%st_start + 1) * (st%d%kpt%end - st%d%kpt%start + 1)
      call sparskit_solver_init(namespace, tr%tdsk_size, tr%tdsk, .true.)
 
#ifndef HAVE_SPARSKIT
      message(1) = 'Octopus was not compiled with support for the SPARSKIT library. This'
      message(2) = 'library is required if the "runge_kutta4" propagator is selected.'
      message(3) = 'Try using a different propagation scheme or recompile with SPARSKIT support.'
      call messages_fatal(3, namespace=namespace)
#endif
 
      call messages_experimental("Runge-Kutta 4 propagator", namespace=namespace)
    case (prop_runge_kutta2)
 
      ! set up pointer for zmf_dotu_aux, zmf_nrm2_aux
      call mesh_init_mesh_aux(gr)
      sp_distdot_mode = 2
      tr%tdsk_size = st%d%dim * gr%np * (st%st_end - st%st_start + 1) * (st%d%kpt%end - st%d%kpt%start + 1)
      call sparskit_solver_init(namespace, tr%tdsk_size, tr%tdsk, .true.)
 
#ifndef HAVE_SPARSKIT
      message(1) = 'Octopus was not compiled with support for the SPARSKIT library. This'
      message(2) = 'library is required if the "runge_kutta2" propagator is selected.'
      message(3) = 'Try using a different propagation scheme or recompile with SPARSKIT support.'
      call messages_fatal(3, namespace=namespace)
#endif
 
      call messages_experimental("Runge-Kutta 2 propagator", namespace=namespace)
    case (prop_crank_nicolson_sparskit)
      ! set up pointer for zmf_dotu_aux
      call mesh_init_mesh_aux(gr)
      sp_distdot_mode = 1
      tr%tdsk_size = st%d%dim*gr%np
      call sparskit_solver_init(namespace, st%d%dim*gr%np, tr%tdsk, .true.)
 
#ifndef HAVE_SPARSKIT
      message(1) = 'Octopus was not compiled with support for the SPARSKIT library. This'
      message(2) = 'library is required if the "crank_nicolson_sparskit" propagator is selected.'
      message(3) = 'Try using a different propagation scheme or recompile with SPARSKIT support.'
      call messages_fatal(3, namespace=namespace)
#endif
    case (prop_magnus)
      call messages_experimental("Magnus propagator", namespace=namespace)
      if (family_is_mgga_with_exc) then
        message(1) = "Magnus propagator with MGGA"
        call messages_fatal(1, namespace=namespace)
      end if
    case (prop_qoct_tddft_propagator)
      call messages_experimental("QOCT+TDDFT propagator", namespace=namespace)
    case (prop_explicit_runge_kutta4)
      call messages_experimental("explicit Runge-Kutta 4 propagator", namespace=namespace)
    case (prop_cfmagnus4)
      call messages_experimental("Commutator-Free Magnus propagator", namespace=namespace)
    case (prop_explicit_magnus3)
      call messages_experimental("Explicit Magnus 3 propagator", namespace=namespace)
    case default
      call messages_input_error(namespace, 'TDPropagator')
    end select
    call messages_print_var_option('TDPropagator', tr%method, namespace=namespace)
 
    if (have_fields) then
      if (tr%method /= prop_etrs .and. &
        tr%method /= prop_aetrs .and. &
        tr%method /= prop_exponential_midpoint .and. &
        tr%method /= prop_exponential_midpoint_predictor .and. &
        tr%method /= prop_cfmagnus4 .and. &
        tr%method /= prop_explicit_magnus3 .and. &
        tr%method /= prop_qoct_tddft_propagator .and. &
        tr%method /= prop_crank_nicolson .and. &
        tr%method /= prop_runge_kutta4 .and. &
        tr%method /= prop_explicit_runge_kutta4 .and. &
        tr%method /= prop_runge_kutta2 .and. &
        tr%method /= prop_crank_nicolson_sparskit) then
        message(1) = "To move the ions or put in a gauge field, use one of the following propagators:"
        message(2) = "etrs, aetrs, cfmagnus4, crank-nicolson, rk2, rk4, or exp_mid."
        call messages_fatal(2, namespace=namespace)
      end if
    end if
 
    if (relax_cell) then
      if (tr%method /= prop_etrs .and. &
        tr%method /= prop_aetrs .and. &
        tr%method /= prop_exponential_midpoint .and. &
        tr%method /= prop_exponential_midpoint_predictor .and. &
        tr%method /= prop_cfmagnus4 .and. &
        tr%method /= prop_explicit_magnus3 .and. &
        tr%method /= prop_crank_nicolson .and. &
        tr%method /= prop_crank_nicolson_sparskit) then
        message(1) = "To relax the cell, use the one of the following propagators:"
        message(2) = "etrs, aetrs, cfmagnus4, crank-nicolson, or exp_mid."
        call messages_fatal(2, namespace=namespace)
      end if
    end if
 
    select case (tr%method)
    case (prop_cfmagnus4)
      call ks_pot%init_interpolation(tr%vks_old, order = 4)
    case default
      call ks_pot%init_interpolation(tr%vks_old)
    end select
 
    call exponential_init(tr%te, namespace) ! initialize propagator
 
    call messages_obsolete_variable(namespace, 'TDSelfConsistentSteps', 'TDStepsWithSelfConsistency')
 
    !%Variable TDStepsWithSelfConsistency
    !%Type integer
    !%Default 0
    !%Section Time-Dependent::Propagation
    !%Description
    !% Since the KS propagator is non-linear, each propagation step
    !% should be performed self-consistently.  In practice, for most
    !% purposes this is not necessary, except perhaps in the first
    !% iterations. This variable holds the number of propagation steps
    !% for which the propagation is done self-consistently.
    !%
    !% The special value <tt>all_steps</tt> forces self-consistency to
    !% be imposed on all propagation steps. A value of 0 means that
    !% self-consistency will not be imposed.  The default is 0.
    !%Option all_steps -1
    !% Self-consistency is imposed for all propagation steps.
    !%End
 
    call parse_variable(namespace, 'TDStepsWithSelfConsistency', 0, tr%scf_propagation_steps)
 
    if (tr%scf_propagation_steps == -1) tr%scf_propagation_steps = huge(tr%scf_propagation_steps)
    if (tr%scf_propagation_steps < 0) call messages_input_error(namespace, 'TDStepsWithSelfConsistency', 'Cannot be negative')
 
    if (tr%scf_propagation_steps /= 0) then
      call messages_experimental('TDStepsWithSelfConsistency', namespace=namespace)
 
      if (tr%method /= prop_etrs .and. tr%method /= prop_aetrs .and. tr%method /= prop_cfmagnus4) then
        call messages_write('TDStepsWithSelfConsistency only works with the ETRS/AETRS/CFM4 propagators')
        call messages_fatal(namespace=namespace)
      end if
    end if
 
    !%Variable TDSCFThreshold
    !%Type float
    !%Default 1.0e-6
    !%Section Time-Dependent::Propagation
    !%Description
    !% Since the KS propagator is non-linear, each propagation step
    !% should be performed self-consistently.  In practice, for most
    !% purposes this is not necessary, except perhaps in the first
    !% iterations. This variable holds the number of propagation steps
    !% for which the propagation is done self-consistently.
    !%
    !% The self consistency has to be measured against some accuracy
    !% threshold. This variable controls the value of that threshold.
    !%End
    call parse_variable(namespace, 'TDSCFThreshold', 1.0e-6_real64, tr%scf_threshold)
 
    pop_sub(propagator_elec_init)
  end subroutine propagator_elec_init
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine propagator_elec_set_scf_prop(tr, threshold)
    type(propagator_base_t), intent(inout) :: tr
    real(real64), optional,         intent(in)    :: threshold
 
    push_sub(propagator_elec_set_scf_prop)
 
    tr%scf_propagation_steps = huge(1)
    if (present(threshold)) then
      tr%scf_threshold = threshold
    end if
 
    pop_sub(propagator_elec_set_scf_prop)
  end subroutine propagator_elec_set_scf_prop
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine propagator_elec_remove_scf_prop(tr)
    type(propagator_base_t), intent(inout) :: tr
 
    push_sub(propagator_elec_remove_scf_prop)
 
    tr%scf_propagation_steps = -1
 
    pop_sub(propagator_elec_remove_scf_prop)
  end subroutine propagator_elec_remove_scf_prop
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  subroutine propagator_elec_end(tr)
    type(propagator_base_t), intent(inout) :: tr
 
    push_sub(propagator_elec_end)
 
    call potential_interpolation_end(tr%vks_old)
 
    select case (tr%method)
    case (prop_runge_kutta4, prop_runge_kutta2)
      call propagator_rk_end()
      call sparskit_solver_end(tr%tdsk)
 
    case (prop_crank_nicolson_sparskit)
      call sparskit_solver_end(tr%tdsk)
 
    end select
 
    if (allocated(tr%st)) then
      call states_elec_end(tr%st)
      safe_deallocate_a(tr%st)
    end if
 
    pop_sub(propagator_elec_end)
  end subroutine propagator_elec_end
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
  ! ---------------------------------------------------------
  subroutine propagator_elec_dt(ks, namespace, space, hm, gr, st, tr, time, dt, nt, &
    ions_dyn, ions, ext_partners, mc, outp, write_handler, scsteps, update_energy, qcchi)
    type(v_ks_t),                        target, intent(inout) :: ks
    type(namespace_t),                           intent(in)    :: namespace
    type(electron_space_t),                      intent(in)    :: space
    type(hamiltonian_elec_t),            target, intent(inout) :: hm
    type(grid_t),                        target, intent(inout) :: gr
    type(states_elec_t),                 target, intent(inout) :: st
    type(propagator_base_t),             target, intent(inout) :: tr
    real(real64),                                intent(in)    :: time
    real(real64),                                intent(in)    :: dt
    integer,                                     intent(in)    :: nt
    type(ion_dynamics_t),                        intent(inout) :: ions_dyn
    type(ions_t),                                intent(inout) :: ions
    type(partner_list_t),                        intent(in)    :: ext_partners
    type(multicomm_t),                           intent(inout) :: mc
    type(output_t),                              intent(in)    :: outp
    type(td_write_t),                            intent(in)    :: write_handler
    integer,                   optional,         intent(out)   :: scsteps
    logical,                   optional,         intent(in)    :: update_energy
    type(opt_control_state_t), optional, target, intent(inout) :: qcchi
    type(gauge_field_t), pointer :: gfield
    type(lasers_t),      pointer :: lasers
    logical :: generate, update_energy_
    real(real64)   :: am(space%dim)
 
    call profiling_in("TD_PROPAGATOR")
    push_sub(propagator_elec_dt)
 
    update_energy_ = optional_default(update_energy, .true.)
 
    call hm%ks_pot%interpolation_new(tr%vks_old, time, dt)
 
    if (present(scsteps)) scsteps = 1
 
    select case (tr%method)
    case (prop_etrs)
      if (self_consistent_step()) then
        call td_etrs_sc(ks, namespace, space, hm, ext_partners, gr, st, tr, time, dt, &
          ions_dyn, ions, mc, tr%scf_threshold, scsteps)
      else
        call td_etrs(ks, namespace, space, hm, ext_partners, gr, st, tr, time, dt, &
          ions_dyn, ions, mc)
      end if
    case (prop_aetrs)
      if (self_consistent_step()) then
        call td_aetrs_sc(ks, namespace, space, hm, ext_partners, gr, st, tr, time, dt, &
          ions_dyn, ions, mc, tr%scf_threshold, scsteps)
      else
        call td_aetrs(namespace, space, hm, gr, st, tr, time, dt, ions_dyn, ions, ext_partners, mc)
      end if
    case (prop_caetrs)
      call td_caetrs(ks, namespace, space, hm, ext_partners, gr, st, tr, time, dt, ions_dyn, ions, mc)
    case (prop_exponential_midpoint)
      call exponential_midpoint(hm, namespace, space, gr, st, tr, time, dt, ions_dyn, ions, ext_partners, mc)
    case (prop_exponential_midpoint_predictor)
      call exponential_midpoint_predictor(hm, ks, namespace, space, gr, st, tr, time, dt, ions_dyn, ions, &
        ext_partners, mc)
    case (prop_crank_nicolson)
      call td_crank_nicolson(hm, namespace, space, gr, st, tr, time, dt, ions_dyn, ions, ext_partners, mc, .false.)
    case (prop_runge_kutta4)
      call td_runge_kutta4(ks, namespace, space, hm, gr, st, tr, time, dt, ions_dyn, ions, ext_partners)
    case (prop_runge_kutta2)
      call td_runge_kutta2(ks, namespace, space, hm, gr, st, tr, time, dt, ions_dyn, ions, ext_partners)
    case (prop_crank_nicolson_sparskit)
      call td_crank_nicolson(hm, namespace, space, gr, st, tr, time, dt, ions_dyn, ions, ext_partners, mc, .true.)
    case (prop_magnus)
      call td_magnus(hm, ext_partners, gr, st, tr, namespace, time, dt)
    case (prop_qoct_tddft_propagator)
      call td_qoct_tddft_propagator(hm, namespace, space, gr, st, tr, time, dt, ions_dyn, ions, ext_partners, mc)
    case (prop_explicit_runge_kutta4)
      if (present(qcchi)) then
        call td_explicit_runge_kutta4(ks, namespace, space, hm, gr, st, time, dt, ions_dyn, ions, ext_partners, qcchi)
      else
        call td_explicit_runge_kutta4(ks, namespace, space, hm, gr, st, time, dt, ions_dyn, ions, ext_partners)
      end if
    case (prop_cfmagnus4)
      if (self_consistent_step()) then
        call td_cfmagnus4_sc(ks, namespace, space, hm, gr, st, tr, time, dt, ions_dyn, ions, ext_partners, mc, &
          nt, tr%scf_threshold, scsteps)
      else
        call td_cfmagnus4(ks, namespace, space, hm, gr, st, tr, time, dt, ions_dyn, ions, ext_partners, mc, nt)
      end if
    case (prop_explicit_magnus3)
      call td_explicit_magnus3(ks, namespace, space, hm, gr, st, tr, time, dt, ions_dyn, ions, ext_partners, mc)
    end select
 
    generate = .false.
    if (ions_dyn%is_active()) then
      if (.not. propagator_elec_ions_are_propagated(tr)) then
        call propagation_ops_elec_propagate_ions_and_cell(gr, hm, st, namespace, space, &
          ions_dyn, ions, mc, abs(nt*dt), ions_dyn%ionic_scale*dt)
        generate = .true.
      end if
    end if
 
    gfield => list_get_gauge_field(ext_partners)
    if(associated(gfield)) then
      if (gauge_field_is_propagated(gfield) .and. .not. propagator_elec_ions_are_propagated(tr)) then
        call gauge_field_do_algorithmic_operation(gfield, op_verlet_compute_acc, dt, time)
        call propagation_ops_elec_update_hamiltonian(namespace, space, st, gr, hm, ext_partners, time)
      end if
    end if
 
    if (generate .or. ions%has_time_dependent_species()) then
      call hamiltonian_elec_epot_generate(hm, namespace,  space, gr, ions, ext_partners, st, time = abs(nt*dt))
    end if
 
    call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
      calc_eigenval = update_energy_, time = abs(nt*dt), calc_energy = update_energy_)
 
    ! compute time integral of the paramagnetic current
    call ks%v_ks_photons%add_mean_field(ks%gr, space, hm%kpoints, st, abs(nt*dt), dt)
 
    if (update_energy_) call energy_calc_total(namespace, space, hm, gr, st, ext_partners, iunit = -1)
 
    ! Recalculate forces, update velocities...
    if ((ions_dyn%ions_move() .and. tr%method .ne. prop_explicit_runge_kutta4) .or. &
      outp%what(option__output__forces) .or. write_handler%out(out_separate_forces)%write) then
      call forces_calculate(gr, namespace, ions, hm, ext_partners, st, ks, t = abs(nt*dt), dt = dt)
    end if
 
    if (ions_dyn%cell_relax() .or. outp%what(option__output__stress)) then
      call stress_calculate(namespace, gr, hm, st, ions, ks, ext_partners)
      if(ions_dyn%cell_relax()) then
        call ions_dyn%update_stress(ions%space, st%stress_tensors%total, ions%latt%rlattice, ions%latt%rcell_volume)
      end if
    end if
 
    if( ions_dyn%is_active()) then
      call ion_dynamics_propagate_vel(ions_dyn, ions, atoms_moved = generate)
      call ions%update_kinetic_energy()
    end if
 
    if(associated(gfield)) then
      if(gauge_field_is_propagated(gfield)) then
        call gauge_field_get_force(gfield, gr, st%d%spin_channels, st%current, ks%xc%lrc)
        call gauge_field_do_algorithmic_operation(gfield, op_verlet_compute_vel, dt, time)
      end if
    end if
 
    !We update the occupation matrices
    call lda_u_update_occ_matrices(hm%lda_u, namespace, gr, st, hm%phase, hm%energy)
 
    !Updates Nondipole Vector potential
    lasers => list_get_lasers(ext_partners)
    if(associated(lasers)) then
      if (lasers_with_nondipole_field(lasers)) then
        call lasers_nondipole_laser_field_step(lasers, am, time)
        call lasers_set_nondipole_parameters(lasers,am,time)
      end if
    end if
    pop_sub(propagator_elec_dt)
    call profiling_out("TD_PROPAGATOR")
 
  contains
 
    ! ---------------------------------------------------------
    logical pure function self_consistent_step() result(scs)
      scs = (hm%theory_level /= independent_particles) .and. (time <= tr%scf_propagation_steps*abs(dt) + m_epsilon)
    end function self_consistent_step
    ! ---------------------------------------------------------
 
  end subroutine propagator_elec_dt
  ! ---------------------------------------------------------
 
 
 
 
  ! ---------------------------------------------------------
  logical pure function propagator_elec_ions_are_propagated(tr) result(propagated)
    type(propagator_base_t), intent(in) :: tr
 
    select case (tr%method)
    case (prop_etrs, prop_aetrs, prop_caetrs, prop_explicit_runge_kutta4)
      propagated = .true.
    case default
      propagated = .false.
    end select
 
  end function propagator_elec_ions_are_propagated
  ! ---------------------------------------------------------
 
 
  ! ---------------------------------------------------------
 
  subroutine propagator_elec_dt_bo(scf, namespace, space, gr, ks, st, hm, ions, ext_partners, mc, iter, dt, ions_dyn, scsteps)
    type(scf_t),              intent(inout) :: scf
    type(namespace_t),        intent(in)    :: namespace
    type(electron_space_t),   intent(in)    :: space
    type(grid_t),             intent(inout) :: gr
    type(v_ks_t),             intent(inout) :: ks
    type(states_elec_t),      intent(inout) :: st
    type(hamiltonian_elec_t), intent(inout) :: hm
    type(ions_t),             intent(inout) :: ions
    type(partner_list_t),     intent(in)    :: ext_partners
    type(multicomm_t),        intent(inout) :: mc
    integer,                  intent(in)    :: iter
    real(real64),             intent(in)    :: dt
    type(ion_dynamics_t),     intent(inout) :: ions_dyn
    integer,                  intent(inout) :: scsteps
 
    type(gauge_field_t), pointer :: gfield
 
    push_sub(propagator_elec_dt_bo)
 
    ! move the hamiltonian to time t
    call propagation_ops_elec_propagate_ions_and_cell(gr, hm, st, namespace, space, &
      ions_dyn, ions, mc, iter*dt, dt)
 
    call hamiltonian_elec_epot_generate(hm, namespace, space, gr, ions, ext_partners, st, time = iter*dt)
    ! now calculate the eigenfunctions
    call scf_run(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, &
      verbosity = verb_compact, iters_done = scsteps)
 
    ! Stress tensor calculation and uptade
    if (ions_dyn%cell_relax()) then
      if (.not. scf%calc_stress) then
        call stress_calculate(namespace, gr, hm, st, ions, ks, ext_partners)
      end if
      call ions_dyn%update_stress(ions%space, st%stress_tensors%total, ions%latt%rlattice, ions%latt%rcell_volume)
    end if
 
    gfield => list_get_gauge_field(ext_partners)
    if(associated(gfield)) then
      if (gauge_field_is_propagated(gfield)) then
        call gauge_field_do_algorithmic_operation(gfield, op_verlet_compute_acc, dt, iter*dt)
      end if
    end if
 
    !TODO: we should update the occupation matrices here
    if (hm%lda_u_level /= dft_u_none) then
      call messages_not_implemented("DFT+U with propagator_elec_dt_bo", namespace=namespace)
    end if
 
    call hamiltonian_elec_epot_generate(hm, namespace, space, gr, ions, ext_partners, st, time = iter*dt)
 
    ! update Hamiltonian and eigenvalues (fermi is *not* called)
    call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
      calc_eigenval = .true., time = iter*dt, calc_energy = .true.)
 
    ! Get the energies.
    call energy_calc_total(namespace, space, hm, gr, st, ext_partners, iunit = -1)
 
    call ion_dynamics_propagate_vel(ions_dyn, ions)
    call hamiltonian_elec_epot_generate(hm, namespace, space, gr, ions, ext_partners, st, time = iter*dt)
    call ions%update_kinetic_energy()
 
    if(associated(gfield)) then
      if (gauge_field_is_propagated(gfield)) then
        call gauge_field_do_algorithmic_operation(gfield, op_verlet_compute_vel, dt, iter*dt)
      end if
    end if
 
    pop_sub(propagator_elec_dt_bo)
  end subroutine propagator_elec_dt_bo
 
end module propagator_elec_oct_m
 
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: