propagator_magnus.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_magnus_oct_m
  use absorbing_boundaries_oct_m
  use batch_oct_m
  use batch_ops_oct_m
  use debug_oct_m
  use density_oct_m
  use electron_space_oct_m
  use exponential_oct_m
  use ext_partner_list_oct_m
  use gauge_field_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_abst_oct_m
  use hamiltonian_elec_oct_m
  use hamiltonian_elec_base_oct_m
  use interaction_partner_oct_m
  use ion_dynamics_oct_m
  use ions_oct_m
  use, intrinsic :: iso_fortran_env
  use ks_potential_oct_m
  use lda_u_io_oct_m
  use lda_u_oct_m
  use lasers_oct_m
  use messages_oct_m
  use mesh_oct_m
  use multicomm_oct_m
  use namespace_oct_m
  use parser_oct_m
  use potential_interpolation_oct_m
  use profiling_oct_m
  use propagation_ops_elec_oct_m
  use propagator_base_oct_m
  use propagator_rk_oct_m
  use space_oct_m
  use states_elec_oct_m
  use states_elec_group_oct_m
  use v_ks_oct_m
  use wfs_elec_oct_m
  use xc_oct_m
 
  implicit none
 
  private
 
  public ::                    &
    td_magnus,                 &
    td_cfmagnus4,              &
    td_cfmagnus4_sc,           &
    td_explicit_magnus3
 
  
  type, extends(hamiltonian_operator_t) :: magnus4_operator_t
    private
    real(real64), allocatable :: vmagnus(:, :, :)
  contains
    procedure :: apply => magnus4_operator_apply
  end type magnus4_operator_t
 
  interface magnus4_operator_t
    procedure magnus4_operator_constructor
  end interface magnus4_operator_t
 
  
  type, extends(hamiltonian_operator_t) :: magnus3_linear_operator_t
    private
    type(hamiltonian_elec_t), pointer :: hm1 => null()
    type(hamiltonian_elec_t), pointer :: hm2 => null()
    real(real64) :: c1 = m_zero
    real(real64) :: c2 = m_zero
  contains
    procedure :: apply => magnus3_linear_operator_apply
  end type magnus3_linear_operator_t
 
  interface magnus3_linear_operator_t
    procedure magnus3_linear_operator_constructor
  end interface magnus3_linear_operator_t
 
  
  type, extends(hamiltonian_operator_t) :: magnus3_operator_t
    private
    type(hamiltonian_elec_t), pointer :: hm1 => null()
    type(hamiltonian_elec_t), pointer :: hm2 => null()
    type(hamiltonian_elec_t), pointer :: hm3 => null()
    type(hamiltonian_elec_t), pointer :: hm4 => null()
    real(real64) :: dt = m_zero
  contains
    procedure :: apply => magnus3_operator_apply
  end type magnus3_operator_t
 
  interface magnus3_operator_t
    procedure magnus3_operator_constructor
  end interface magnus3_operator_t
 
 
contains
 
  ! ---------------------------------------------------------
  subroutine td_magnus(hm, ext_partners, gr, st, tr, namespace, time, dt)
    type(hamiltonian_elec_t),         intent(inout) :: hm
    type(partner_list_t),             intent(in)    :: ext_partners
    type(grid_t),                     intent(in)    :: gr
    type(states_elec_t),              intent(inout) :: st
    type(propagator_base_t),          intent(inout) :: tr
    type(namespace_t),                intent(in)    :: namespace
    real(real64),                     intent(in)    :: time
    real(real64),                     intent(in)    :: dt
 
    integer :: j, is, i
    real(real64) :: atime(2)
    real(real64), allocatable :: vaux(:, :, :), vmagnus(:, :, :), pot(:)
    type(lasers_t), pointer :: lasers
    class(operator_t), pointer :: op
 
    push_sub(propagator_dt.td_magnus)
 
    assert(.not. family_is_mgga_with_exc(hm%xc))
 
    safe_allocate(vaux(1:gr%np, 1:st%d%nspin, 1:2))
    safe_allocate(vmagnus(1:gr%np, 1:st%d%nspin, 1:2))
 
    atime(1) = (m_half-sqrt(m_three)/6.0_real64)*dt
    atime(2) = (m_half+sqrt(m_three)/6.0_real64)*dt
 
    if (hm%theory_level /= independent_particles) then
      do j = 1, 2
        call potential_interpolation_interpolate(tr%vks_old, 3, time, dt, time - dt + atime(j), vaux(:, :, j))
      end do
    else
      vaux = m_zero
    end if
 
    do j = 1, 2
      ! WARNING: This should be carefully tested, and extended to allow for velocity-gauge laser fields.
      lasers => list_get_lasers(ext_partners)
      if(associated(lasers)) then
        do i = 1, lasers%no_lasers
          select case (laser_kind(lasers%lasers(i)))
          case (e_field_electric)
            safe_allocate(pot(1:gr%np))
            pot = m_zero
            call laser_potential(lasers%lasers(i), gr, pot, time - dt + atime(j))
            do is = 1, st%d%nspin
              vaux(:, is, j) = vaux(:, is, j) + pot(:)
            end do
            safe_deallocate_a(pot)
          case (e_field_magnetic, e_field_vector_potential)
            write(message(1),'(a)') 'The Magnus propagator cannot be used with magnetic fields, or'
            write(message(2),'(a)') 'with an electric field described in the velocity gauge.'
            call messages_fatal(2, namespace=namespace)
          end select
        end do
      end if
    end do
 
    vmagnus(:, :, 2)  = m_half*(vaux(:, :, 1) + vaux(:, :, 2))
    vmagnus(:, :, 1) = (sqrt(m_three)/12.0_real64)*dt*(vaux(:, :, 2) - vaux(:, :, 1))
 
    op => magnus4_operator_t(namespace, gr, hm, vmagnus)
    call propagation_ops_elec_fuse_density_exp_apply(tr%te, namespace, st, gr, hm, dt, op=op)
    safe_deallocate_p(op)
 
    safe_deallocate_a(vaux)
    safe_deallocate_a(vmagnus)
    pop_sub(propagator_dt.td_magnus)
  end subroutine td_magnus
 
  function magnus4_operator_constructor(namespace, mesh, hm, vmagnus) result(this)
    type(namespace_t), target,         intent(in) :: namespace
    class(mesh_t), target,             intent(in) :: mesh
    class(hamiltonian_abst_t), target, intent(in) :: hm
    real(real64),                      intent(in) :: vmagnus(:, :, :)
    type(magnus4_operator_t), pointer :: this
 
    push_sub(magnus4_operator_constructor)
 
    allocate(this)
    call this%init(namespace, mesh, hm)
 
    select type(hm)
    class is (hamiltonian_elec_t)
      assert(size(vmagnus, 1) >= mesh%np)
      assert(size(vmagnus, 2) == hm%d%nspin)
      assert(size(vmagnus, 3) == 2)
    class default
      write(message(1), '(a)') 'Magnus operators only implemented for electrons at the moment'
      call messages_fatal(1, namespace=namespace)
    end select
    safe_allocate_source(this%vmagnus, vmagnus)
 
    pop_sub(magnus4_operator_constructor)
  end function magnus4_operator_constructor
 
  subroutine magnus4_operator_apply(this, psib, hpsib)
    class(magnus4_operator_t), intent(in) :: this
    class(batch_t), intent(inout) :: psib
    class(batch_t), intent(inout) :: hpsib
 
    integer :: ispin
    type(wfs_elec_t) :: auxpsib, aux2psib
 
    push_sub(magnus4_operator_apply)
 
    select type(hm => this%hm)
    class is (hamiltonian_elec_t)
      select type (psib)
      class is (wfs_elec_t)
        select type (hpsib)
        class is (wfs_elec_t)
          ! We will assume, for the moment, no spinors.
          if (hm%d%dim /= 1) call messages_not_implemented("Magnus with spinors", namespace=this%namespace)
 
          assert(psib%nst == hpsib%nst)
 
          ispin = hm%d%get_spin_index(psib%ik)
 
          call hpsib%copy_to(auxpsib, copy_data=.false.)
          call hpsib%copy_to(aux2psib, copy_data=.false.)
 
          ! Compute (T + Vnl)|psi> and store it
          call zhamiltonian_elec_apply_batch(hm, this%namespace, this%mesh, psib, auxpsib, &
            terms = term_kinetic + term_non_local_potential)
 
          ! H|psi>  =  (T + Vnl)|psi> + Vpsl|psi> + Vmagnus(t2)|psi> + Vborders|psi>
          call auxpsib%copy_data_to(this%mesh%np, hpsib)
          call zhamiltonian_elec_external(hm, this%mesh, psib, hpsib)
          if (hm%abs_boundaries%abtype == imaginary_absorbing) then
            call batch_mul(this%mesh%np, hm%abs_boundaries%mf(1:this%mesh%np), psib, aux2psib)
            call batch_axpy(this%mesh%np, m_zi, aux2psib, hpsib)
          end if
          call batch_mul(this%mesh%np, this%vmagnus(1:this%mesh%np, ispin, 2), psib, aux2psib)
          call batch_axpy(this%mesh%np, m_one, aux2psib, hpsib)
 
          ! Add first term of the commutator:  - i Vmagnus(t1) (T + Vnl) |psi>
          call batch_mul(this%mesh%np, this%vmagnus(1:this%mesh%np, ispin, 1), auxpsib, aux2psib)
          call batch_axpy(this%mesh%np, -m_zi, aux2psib, hpsib)
 
          ! Add second term of commutator:  i (T + Vnl) Vmagnus(t1) |psi>
          call batch_mul(this%mesh%np, this%vmagnus(1:this%mesh%np, ispin, 1), psib, auxpsib)
          call zhamiltonian_elec_apply_batch(hm, this%namespace, this%mesh, auxpsib, aux2psib, &
            terms = term_kinetic + term_non_local_potential)
          call batch_axpy(this%mesh%np, m_zi, aux2psib, hpsib)
 
          call auxpsib%end()
          call aux2psib%end()
        class default
          message(1) = "Internal error: imcompatible batch_t in magnus4_operator_apply for argument hpsib."
          call messages_fatal(1, namespace=this%namespace)
        end select
      class default
        message(1) = "Internal error: imcompatible batch_t in magnus4_operator_apply for argument psib."
        call messages_fatal(1, namespace=this%namespace)
      end select
    end select
 
    pop_sub(magnus4_operator_apply)
  end subroutine magnus4_operator_apply
 
  ! ---------------------------------------------------------
  function magnus3_linear_operator_constructor(namespace, mesh, hm1, hm2, c1, c2) result(this)
    type(namespace_t), target,        intent(in) :: namespace
    class(mesh_t), target,            intent(in) :: mesh
    type(hamiltonian_elec_t), target, intent(in) :: hm1
    type(hamiltonian_elec_t), target, intent(in) :: hm2
    real(real64),                     intent(in) :: c1
    real(real64),                     intent(in) :: c2
    type(magnus3_linear_operator_t), pointer :: this
 
    push_sub(magnus3_linear_operator_constructor)
 
    allocate(this)
    call this%init(namespace, mesh, hm1)
    this%hm1 => hm1
    this%hm2 => hm2
    this%c1 = c1
    this%c2 = c2
 
    pop_sub(magnus3_linear_operator_constructor)
  end function magnus3_linear_operator_constructor
 
  ! ---------------------------------------------------------
  subroutine magnus3_apply_stage_in_ref_phase(namespace, mesh, hm_ref, hm_stage, psib, hpsib)
    type(namespace_t),         intent(in)    :: namespace
    class(mesh_t),             intent(in)    :: mesh
    type(hamiltonian_elec_t),  intent(in)    :: hm_ref
    type(hamiltonian_elec_t),  intent(in)    :: hm_stage
    type(wfs_elec_t),          intent(inout) :: psib
    type(wfs_elec_t),          intent(inout) :: hpsib
 
    push_sub(magnus3_apply_stage_in_ref_phase)
 
    if (hm_ref%phase%is_allocated()) call hm_ref%phase%unset_phase_corr(mesh, psib)
    if (hm_stage%phase%is_allocated()) call hm_stage%phase%set_phase_corr(mesh, psib)
 
    call zhamiltonian_elec_apply_batch(hm_stage, namespace, mesh, psib, hpsib)
 
    if (hm_stage%phase%is_allocated()) then
      call hm_stage%phase%unset_phase_corr(mesh, hpsib)
      call hm_stage%phase%unset_phase_corr(mesh, psib)
    end if
    if (hm_ref%phase%is_allocated()) then
      call hm_ref%phase%set_phase_corr(mesh, hpsib)
      call hm_ref%phase%set_phase_corr(mesh, psib)
    end if
 
    pop_sub(magnus3_apply_stage_in_ref_phase)
  end subroutine magnus3_apply_stage_in_ref_phase
 
  ! ---------------------------------------------------------
  subroutine magnus3_linear_operator_apply(this, psib, hpsib)
    class(magnus3_linear_operator_t), intent(in) :: this
    class(batch_t), intent(inout) :: psib
    class(batch_t), intent(inout) :: hpsib
 
    type(wfs_elec_t) :: auxpsib
 
    push_sub(magnus3_linear_operator_apply)
 
    select type (hm => this%hm)
    class is (hamiltonian_elec_t)
      select type (psib)
      class is (wfs_elec_t)
        select type (hpsib)
        class is (wfs_elec_t)
          if (hm%d%dim /= 1) call messages_not_implemented("Explicit Magnus 3 with spinors", namespace=this%namespace)
 
          assert(psib%nst == hpsib%nst)
 
          call hpsib%copy_to(auxpsib, copy_data=.false.)
 
          call magnus3_apply_stage_in_ref_phase(this%namespace, this%mesh, this%hm1, this%hm1, psib, hpsib)
          call batch_scal(this%mesh%np, this%c1, hpsib)
 
          call magnus3_apply_stage_in_ref_phase(this%namespace, this%mesh, this%hm1, this%hm2, psib, auxpsib)
          call batch_axpy(this%mesh%np, this%c2, auxpsib, hpsib)
 
          call auxpsib%end()
        class default
          message(1) = "Internal error: imcompatible batch_t in magnus3_linear_operator_apply for argument hpsib."
          call messages_fatal(1, namespace=this%namespace)
        end select
      class default
        message(1) = "Internal error: imcompatible batch_t in magnus3_linear_operator_apply for argument psib."
        call messages_fatal(1, namespace=this%namespace)
      end select
    end select
 
    pop_sub(magnus3_linear_operator_apply)
  end subroutine magnus3_linear_operator_apply
 
  ! ---------------------------------------------------------
  function magnus3_operator_constructor(namespace, mesh, hm1, hm2, hm3, hm4, dt) result(this)
    type(namespace_t), target,        intent(in) :: namespace
    class(mesh_t), target,            intent(in) :: mesh
    type(hamiltonian_elec_t), target, intent(in) :: hm1
    type(hamiltonian_elec_t), target, intent(in) :: hm2
    type(hamiltonian_elec_t), target, intent(in) :: hm3
    type(hamiltonian_elec_t), target, intent(in) :: hm4
    real(real64),                     intent(in) :: dt
    type(magnus3_operator_t), pointer :: this
 
    push_sub(magnus3_operator_constructor)
 
    allocate(this)
    call this%init(namespace, mesh, hm1)
    this%hm1 => hm1
    this%hm2 => hm2
    this%hm3 => hm3
    this%hm4 => hm4
    this%dt = dt
 
    pop_sub(magnus3_operator_constructor)
  end function magnus3_operator_constructor
 
  ! ---------------------------------------------------------
  subroutine magnus3_operator_apply(this, psib, hpsib)
    class(magnus3_operator_t), intent(in) :: this
    class(batch_t), intent(inout) :: psib
    class(batch_t), intent(inout) :: hpsib
 
    type(wfs_elec_t) :: k1psib, k2psib, k3psib, k4psib
    type(wfs_elec_t) :: work1psib, work2psib, sum12psib, commpsib
 
    push_sub(magnus3_operator_apply)
 
    select type (hm => this%hm)
    class is (hamiltonian_elec_t)
      select type (psib)
      class is (wfs_elec_t)
        select type (hpsib)
        class is (wfs_elec_t)
          if (hm%d%dim /= 1) call messages_not_implemented("Explicit Magnus 3 with spinors", namespace=this%namespace)
 
          assert(psib%nst == hpsib%nst)
 
          call hpsib%copy_to(k1psib, copy_data=.false.)
          call hpsib%copy_to(k2psib, copy_data=.false.)
          call hpsib%copy_to(k3psib, copy_data=.false.)
          call hpsib%copy_to(k4psib, copy_data=.false.)
          call hpsib%copy_to(work1psib, copy_data=.false.)
          call hpsib%copy_to(work2psib, copy_data=.false.)
          call hpsib%copy_to(sum12psib, copy_data=.false.)
          call hpsib%copy_to(commpsib, copy_data=.false.)
 
          ! Precompute stage derivatives once:
          !   k1 = H1 psi, k2 = H2 psi, k3 = H3 psi, k4 = H4 psi.
          ! no phase treatment needed for first stage
          call zhamiltonian_elec_apply_batch(this%hm1, this%namespace, this%mesh, psib, k1psib)
          call magnus3_apply_stage_in_ref_phase(this%namespace, this%mesh, this%hm1, this%hm2, psib, k2psib)
          call magnus3_apply_stage_in_ref_phase(this%namespace, this%mesh, this%hm1, this%hm3, psib, k3psib)
          call magnus3_apply_stage_in_ref_phase(this%namespace, this%mesh, this%hm1, this%hm4, psib, k4psib)
 
          ! Weighted sum in eq. 22: (k1 + 4 k3 + k4)/6.
          call k1psib%copy_data_to(this%mesh%np, hpsib)
          call batch_scal(this%mesh%np, m_one/6.0_real64, hpsib)
          call batch_axpy(this%mesh%np, m_four/6.0_real64, k3psib, hpsib)
          call batch_axpy(this%mesh%np, m_one/6.0_real64, k4psib, hpsib)
 
          ! Commutator term from eq. 22:
          !   -1/3 [u3, k3] - 1/12 [u4, k4]
          ! where u3 = (k1 + k2)/4 and u4 = k2.
 
          ! [H1 + H2, H3] psi
          call magnus3_apply_stage_in_ref_phase(this%namespace, this%mesh, this%hm1, this%hm1, k3psib, work1psib)
          call magnus3_apply_stage_in_ref_phase(this%namespace, this%mesh, this%hm1, this%hm2, k3psib, work2psib)
          call batch_axpy(this%mesh%np, m_one, work2psib, work1psib)
 
          call k1psib%copy_data_to(this%mesh%np, sum12psib)
          call batch_axpy(this%mesh%np, m_one, k2psib, sum12psib)
          call magnus3_apply_stage_in_ref_phase(this%namespace, this%mesh, this%hm1, this%hm3, sum12psib, commpsib)
          call batch_axpy(this%mesh%np, -m_one, commpsib, work1psib)
 
          ! [H2, H4] psi
          call magnus3_apply_stage_in_ref_phase(this%namespace, this%mesh, this%hm1, this%hm2, k4psib, work2psib)
          call magnus3_apply_stage_in_ref_phase(this%namespace, this%mesh, this%hm1, this%hm4, k2psib, commpsib)
          call batch_axpy(this%mesh%np, -m_one, commpsib, work2psib)
 
          call batch_axpy(this%mesh%np, m_one, work2psib, work1psib)
          call batch_axpy(this%mesh%np, (m_zi*this%dt)/12.0_real64, work1psib, hpsib)
 
          call k1psib%end()
          call k2psib%end()
          call k3psib%end()
          call k4psib%end()
          call work1psib%end()
          call work2psib%end()
          call sum12psib%end()
          call commpsib%end()
        class default
          message(1) = "Internal error: imcompatible batch_t in magnus3_operator_apply for argument hpsib."
          call messages_fatal(1, namespace=this%namespace)
        end select
      class default
        message(1) = "Internal error: imcompatible batch_t in magnus3_operator_apply for argument psib."
        call messages_fatal(1, namespace=this%namespace)
      end select
    end select
 
    pop_sub(magnus3_operator_apply)
  end subroutine magnus3_operator_apply
 
  ! ---------------------------------------------------------
  subroutine td_explicit_magnus3(ks, namespace, space, hm, gr, st, tr, time, dt, ions_dyn, ions, ext_partners, mc)
    type(v_ks_t),                     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
    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(gauge_field_t), pointer :: gfield
    type(states_elec_t) :: st0, st2, st3, st4
    type(hamiltonian_elec_t), target :: hm1, hm2, hm3, hm4
    class(operator_t), pointer :: op
 
    push_sub(propagator_dt.td_explicit_magnus3)
 
    ! Stage u1 / k1 (eq. 22):
    ! Keep Y_n and snapshot H1 = H(t_n, Y_n).
    call states_elec_copy(st0, st)
    call hamiltonian_elec_copy(hm1, hm)
 
    ! Stage u2 = 1/2 k1 -> k2 (eq. 22):
    ! Predict Y at t_n + dt/2 using H1, then build hm2 from that Y at the midpoint
    call states_elec_copy(st2, st0)
    call propagation_ops_elec_fuse_density_exp_apply(tr%te, namespace, st2, gr, hm1, m_half*dt)
 
    call propagation_ops_elec_move_ions(tr%propagation_ops_elec, gr, hm, st, namespace, space, ions_dyn, ions, &
      ext_partners, mc, time - m_half*dt, m_half*dt, save_pos = .true.)
 
    gfield => list_get_gauge_field(ext_partners)
    if (associated(gfield)) then
      call propagation_ops_elec_propagate_gauge_field(tr%propagation_ops_elec, gfield, m_half*dt, time, save_gf=.true.)
    end if
 
    if (hm%theory_level /= independent_particles) then
      call v_ks_calc(ks, namespace, space, hm, st2, ions, ext_partners, &
        calc_current = .false., calc_energy = .false., calc_eigenval = .false.)
    end if
    call propagation_ops_elec_update_hamiltonian(namespace, space, st2, gr, hm, ext_partners, time - m_half*dt)
    call hamiltonian_elec_copy(hm2, hm)
 
    ! Stage u3 = 1/4 (k1 + k2) -> k3 (eq. 22):
    ! Propagate with 1/4 H1 + 1/4 H2 and build hm3 at the midpoint from st3
    call states_elec_copy(st3, st0)
    op => magnus3_linear_operator_t(namespace, gr, hm1, hm2, m_one/4.0_real64, m_one/4.0_real64)
    call propagation_ops_elec_fuse_density_exp_apply(tr%te, namespace, st3, gr, hm1, dt, op=op)
    safe_deallocate_p(op)
 
    if (hm%theory_level /= independent_particles) then
      call v_ks_calc(ks, namespace, space, hm, st3, ions, ext_partners, &
        calc_current = .false., calc_energy = .false., calc_eigenval = .false.)
    end if
    call propagation_ops_elec_update_hamiltonian(namespace, space, st3, gr, hm, ext_partners, time - m_half*dt)
    call hamiltonian_elec_copy(hm3, hm)
 
    ! Stage u4 = k2 -> k4 (eq. 22):
    ! Propagate Y_n with H2 to t_n + dt, then rebuild H4 at the end of the step.
    call states_elec_copy(st4, st0)
    call propagation_ops_elec_fuse_density_exp_apply(tr%te, namespace, st4, gr, hm2, dt)
 
    call propagation_ops_elec_move_ions(tr%propagation_ops_elec, gr, hm, st, namespace, space, ions_dyn, ions, &
      ext_partners, mc, time, m_half*dt)
    if (associated(gfield)) then
      call propagation_ops_elec_propagate_gauge_field(tr%propagation_ops_elec, gfield, m_half*dt, time)
    end if
 
    if (hm%theory_level /= independent_particles) then
      call v_ks_calc(ks, namespace, space, hm, st4, ions, ext_partners, &
        calc_current = .false., calc_energy = .false., calc_eigenval = .false.)
    end if
    call propagation_ops_elec_update_hamiltonian(namespace, space, st4, gr, hm, ext_partners, time)
    call hamiltonian_elec_copy(hm4, hm)
 
    ! Final stage v3 (eq. 22):
    ! Combine k1, k3, k4 and the commutator correction
    !   -1/3 [u3, k3] - 1/12 [u4, k4].
    op => magnus3_operator_t(namespace, gr, hm1, hm2, hm3, hm4, dt)
    call propagation_ops_elec_fuse_density_exp_apply(tr%te, namespace, st, gr, hm1, dt, op=op)
    safe_deallocate_p(op)
 
    ! Restore ions and gauge field to their saved values at t_n.
    call propagation_ops_elec_restore_ions(tr%propagation_ops_elec, ions_dyn, ions)
    if (associated(gfield)) then
      call propagation_ops_elec_restore_gauge_field(tr%propagation_ops_elec, namespace, space, hm, gr, ext_partners)
    end if
 
    call states_elec_end(st0)
    call states_elec_end(st2)
    call states_elec_end(st3)
    call states_elec_end(st4)
    call hamiltonian_elec_end(hm1)
    call hamiltonian_elec_end(hm2)
    call hamiltonian_elec_end(hm3)
    call hamiltonian_elec_end(hm4)
 
    pop_sub(propagator_dt.td_explicit_magnus3)
  end subroutine td_explicit_magnus3
 
  ! ---------------------------------------------------------
  subroutine td_cfmagnus4(ks, namespace, space, hm, gr, st, tr, time, dt, ions_dyn, ions, ext_partners, mc, iter)
    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
    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
    integer,                          intent(in)    :: iter
 
    real(real64) :: alpha1, alpha2, c1, c2, t1, t2, dt1, dt2
    type(hamiltonian_elec_t), target :: hm1, hm2
    class(operator_t), pointer :: op
    type(gauge_field_t), pointer :: gfield
 
    push_sub(propagator_dt.td_cfmagnus4)
 
    alpha1 = (m_three - m_two * sqrt(m_three))/12.0_real64
    alpha2 = (m_three + m_two * sqrt(m_three))/12.0_real64
    c1 = m_half - sqrt(m_three)/6.0_real64
    c2 = m_half + sqrt(m_three)/6.0_real64
    t1 = time - dt + c1*dt
    t2 = time - dt + c2*dt
 
    gfield => list_get_gauge_field(ext_partners)
    dt1 = c1*dt
    dt2 = (c2 - c1)*dt
 
    call propagation_ops_elec_move_ions(tr%propagation_ops_elec, gr, hm, st, namespace, space, ions_dyn, ions, &
      ext_partners, mc, t1, dt1, save_pos=.true.)
    if (associated(gfield)) then
      call propagation_ops_elec_propagate_gauge_field(tr%propagation_ops_elec, gfield, dt1, time, save_gf=.true.)
    end if
 
    call hm%ks_pot%interpolate_potentials(tr%vks_old, 4, time, dt, t1)
    call propagation_ops_elec_update_hamiltonian(namespace, space, st, gr, hm, ext_partners, t1)
    call hamiltonian_elec_copy(hm1, hm)
 
    call propagation_ops_elec_move_ions(tr%propagation_ops_elec, gr, hm, st, namespace, space, ions_dyn, ions, &
      ext_partners, mc, t2, dt2)
    if (associated(gfield)) then
      call propagation_ops_elec_propagate_gauge_field(tr%propagation_ops_elec, gfield, dt2, time)
    end if
 
    call hm%ks_pot%interpolate_potentials(tr%vks_old, 4, time, dt, t2)
    call propagation_ops_elec_update_hamiltonian(namespace, space, st, gr, hm, ext_partners, t2)
    call hamiltonian_elec_copy(hm2, hm)
 
    call propagation_ops_elec_restore_ions(tr%propagation_ops_elec, ions_dyn, ions)
    if (associated(gfield)) then
      call propagation_ops_elec_restore_gauge_field(tr%propagation_ops_elec, namespace, space, hm, gr, ext_partners)
    end if
 
    op => magnus3_linear_operator_t(namespace, gr, hm1, hm2, m_two*alpha2, m_two*alpha1)
    call propagation_ops_elec_exp_apply(tr%te, namespace, st, gr, hm1, m_half*dt, op=op)
    safe_deallocate_p(op)
 
    op => magnus3_linear_operator_t(namespace, gr, hm1, hm2, m_two*alpha1, m_two*alpha2)
    call propagation_ops_elec_fuse_density_exp_apply(tr%te, namespace, st, gr, hm1, m_half*dt, op=op)
    safe_deallocate_p(op)
 
    call hamiltonian_elec_end(hm1)
    call hamiltonian_elec_end(hm2)
 
    pop_sub(propagator_dt.td_cfmagnus4)
  end subroutine td_cfmagnus4
 
  ! ---------------------------------------------------------
  subroutine td_cfmagnus4_sc(ks, namespace, space, hm, gr, st, tr, time, dt, ions_dyn, ions, ext_partners, mc, &
    iter, sctol, scsteps)
    type(v_ks_t),                     intent(inout) :: ks
    type(namespace_t),                intent(in)    :: namespace
    type(electron_space_t),           intent(in)    :: space
    type(hamiltonian_elec_t),         intent(inout) :: hm
    type(grid_t),                     intent(inout) :: gr
    type(states_elec_t),              intent(inout) :: st
    type(propagator_base_t),          intent(inout) :: tr
    real(real64),                     intent(in)    :: time
    real(real64),                     intent(in)    :: dt
    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
    integer,                          intent(in)    :: iter
    real(real64),                     intent(in)    :: sctol
    integer,                optional, intent(out)   :: scsteps
 
    real(real64) :: diff
    integer :: iter_sc
    type(states_elec_t) :: st0
    integer, parameter :: niter = 10
    type(xc_copied_potentials_t) :: vhxc_pred
 
    push_sub(propagator_dt.td_cfmagnus4_sc)
 
    assert(hm%theory_level /= independent_particles)
 
    call messages_new_line()
    call messages_write('        Self-consistency iteration:')
    call messages_info(namespace=namespace)
 
    ! store copy of initial states
    call states_elec_copy(st0, st)
 
    do iter_sc = 1, niter
      ! Keep a copy of the current end-of-step prediction (history 0)
      ! to evaluate SC convergence after the propagation.
      call hm%ks_pot%get_interpolated_potentials(tr%vks_old, 0, storage=vhxc_pred)
 
      call td_cfmagnus4(ks, namespace, space, hm, gr, st, tr, time, dt, ions_dyn, ions, ext_partners, mc, iter)
 
      call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
        time = time, calc_current = .false., calc_energy = .false., calc_eigenval = .false.)
      call lda_u_update_occ_matrices(hm%lda_u, namespace, gr, st, hm%phase, hm%energy)
 
      diff = hm%ks_pot%check_convergence(vhxc_pred, gr, st%rho, st%qtot)
 
      ! Feed the updated end-of-step potential back into history 0 so
      ! the next SC iteration interpolates from the latest prediction.
      call hm%ks_pot%set_interpolated_potentials(tr%vks_old, 0)
 
      call messages_write('          step ')
      call messages_write(iter_sc)
      call messages_write(', residue = ')
      call messages_write(abs(diff), fmt = '(1x,es9.2)')
      call messages_info(namespace=namespace)
 
      if (diff <= sctol) exit
      if (iter_sc /= niter) then
        ! restore initial states for next iteration
        call states_elec_group_copy(st0%d, st0%group, st%group)
      end if
    end do
 
    if (hm%lda_u_level /= dft_u_none) then
      call lda_u_write_u(hm%lda_u, namespace=namespace)
      call lda_u_write_v(hm%lda_u, namespace=namespace)
    end if
 
    call messages_info(namespace=namespace)
 
    if (present(scsteps)) scsteps = iter_sc
 
    call states_elec_end(st0)
 
    pop_sub(propagator_dt.td_cfmagnus4_sc)
  end subroutine td_cfmagnus4_sc
 
end module propagator_magnus_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: