propagator_etrs.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_etrs_oct_m
  use accel_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 grid_oct_m
  use gauge_field_oct_m
  use global_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 ks_potential_oct_m
  use lda_u_oct_m
  use lda_u_io_oct_m
  use math_oct_m
  use messages_oct_m
  use mesh_function_oct_m
  use multicomm_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use propagator_base_oct_m
  use space_oct_m
  use states_elec_calc_oct_m
  use states_elec_dim_oct_m
  use states_elec_group_oct_m
  use states_elec_oct_m
  use types_oct_m
  use v_ks_oct_m
  use wfs_elec_oct_m
  use propagation_ops_elec_oct_m
  use xc_oct_m
 
  implicit none
 
  private
 
  public ::                       &
    td_etrs,                      &
    td_etrs_sc,                   &
    td_aetrs,                     &
    td_aetrs_sc,                  &
    td_caetrs
 
contains
 
  subroutine td_etrs(ks, namespace, space, hm, ext_partners, gr, st, tr, time, dt, &
    ions_dyn, ions, 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),         intent(inout) :: hm
    type(partner_list_t),             intent(in)    :: ext_partners
    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(multicomm_t),                intent(inout) :: mc
 
    push_sub(td_etrs)
 
    if (hm%theory_level == independent_particles) then
      call td_etrs_independent_particles(namespace, space, hm, ext_partners, gr, st, tr, time, dt, &
        ions_dyn, ions, mc)
    else
      call td_etrs_predictor(ks, namespace, space, hm, ext_partners, gr, st, tr, time, dt, &
        ions_dyn, ions, mc)
      call propagation_ops_elec_fuse_density_exp_apply(tr%te, namespace, st, gr, hm, m_half*dt)
    endif
 
    pop_sub(td_etrs)
 
  end subroutine td_etrs
 
 
  subroutine td_etrs_independent_particles(namespace, space, hm, ext_partners, gr, st, tr, time, dt, &
    ions_dyn, ions, mc)
    type(namespace_t),                intent(in)    :: namespace
    type(electron_space_t),           intent(in)    :: space
    type(hamiltonian_elec_t),         intent(inout) :: hm
    type(partner_list_t),             intent(in)    :: ext_partners
    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(multicomm_t),                intent(inout) :: mc
 
    type(gauge_field_t), pointer :: gfield
 
    push_sub(td_etrs_independent_particles)
 
    assert(hm%theory_level == independent_particles)
 
    call propagation_ops_elec_exp_apply(tr%te, namespace, st, gr, hm, m_half*dt)
 
    ! first move the ions to time t
    call propagation_ops_elec_move_ions(tr%propagation_ops_elec, gr, hm, st, namespace, space, ions_dyn, ions, &
      ext_partners, mc, time, dt)
 
    gfield => list_get_gauge_field(ext_partners)
    if(associated(gfield)) then
      call propagation_ops_elec_propagate_gauge_field(tr%propagation_ops_elec, gfield, dt, time)
    end if
 
    call propagation_ops_elec_update_hamiltonian(namespace, space, st, gr, hm, ext_partners, time)
 
    call propagation_ops_elec_fuse_density_exp_apply(tr%te, namespace, st, gr, hm, m_half*dt)
 
    pop_sub(td_etrs_independent_particles)
 
  end subroutine td_etrs_independent_particles
 
 
  subroutine td_etrs_predictor(ks, namespace, space, hm, ext_partners, gr, st, tr, time, dt, &
    ions_dyn, ions, 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),         intent(inout) :: hm
    type(partner_list_t),             intent(in)    :: ext_partners
    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(multicomm_t),                intent(inout) :: mc
 
    type(gauge_field_t), pointer :: gfield
 
    push_sub(td_etrs_predictor)
 
    if (v_ks_is_orbital_dependent(ks) .or. hm%lda_u_level /= dft_u_none) then
      ! Full-step states, allocated once and reused every step.
      ! With ACE, v_ks_calc leaves hm%exxop%st pointing at them.
      if (.not. allocated(tr%st)) then
        safe_allocate(tr%st)
        call states_elec_copy(tr%st, st)
      end if
 
      ! st contains \psi_{n+1/2} = exp(-i dt/2  H_n) \psi_n
      ! tr%st contains \tilde{\psi}_{n+1} = exp(-i dt H_n) psi_n
      ! Predictor density from states propagated by the full time step, \tilde{\psi}_{n+1}
      call propagation_ops_elec_fuse_density_exp_apply(tr%te, namespace, st, gr, hm, m_half*dt, dt, st_td=tr%st)
 
      ! Compute potentials, V_{n+1}, from predictor density and states at dt
      call v_ks_calc(ks, namespace, space, hm, tr%st, ions, ext_partners, &
        calc_current = .false., calc_energy = .false., calc_eigenval = .false.)
    else
      ! st contains \psi_{n+1/2} = exp(-i dt/2  H_n) \psi_n
      ! Predictor density from states propagated by the full time step, \tilde{\psi}_{n+1}
      call propagation_ops_elec_fuse_density_exp_apply(tr%te, namespace, st, gr, hm, m_half*dt, dt)
 
      ! Compute potentials, V_{n+1}, from predictor density at dt
      call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
        calc_current = .false., calc_energy = .false., calc_eigenval = .false.)
    endif
 
    ! Advance ions and gauge field to t_{n+1}. Neither reads the Hamiltonian. The ions use
    ! forces from the end of the previous step, and epot_generate rebuilds the ionic terms
    call propagation_ops_elec_move_ions(tr%propagation_ops_elec, gr, hm, st, namespace, space, ions_dyn, ions, &
      ext_partners, mc, time, dt)
 
    gfield => list_get_gauge_field(ext_partners)
    if(associated(gfield)) then
      call propagation_ops_elec_propagate_gauge_field(tr%propagation_ops_elec, gfield, dt, time)
    end if
 
    ! Assemble H_{n+1} from V_{n+1}, the moved ions and the new gauge field
    if (v_ks_is_orbital_dependent(ks) .or. hm%lda_u_level /= dft_u_none) then
      ! Maxwell and LDA+U require predictor density and states, respectively
      call propagation_ops_elec_update_hamiltonian(namespace, space, tr%st, gr, hm, ext_partners, time)
    else
      call propagation_ops_elec_update_hamiltonian(namespace, space, st, gr, hm, ext_partners, time)
    endif
 
    pop_sub(td_etrs_predictor)
 
  end subroutine td_etrs_predictor
 
 
  subroutine td_etrs_sc(ks, namespace, space, hm, ext_partners, gr, st, tr, time, dt, &
    ions_dyn, ions, mc, 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(partner_list_t),             intent(in)    :: ext_partners
    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(multicomm_t),                intent(inout) :: mc
    real(real64),                     intent(in)    :: sctol
    integer,                optional, intent(out)   :: scsteps
 
    integer, parameter :: niter = 10
    integer      :: ik, ib, iter
    logical      :: orbital_dependent
    real(real64) :: diff
    class(wfs_elec_t), allocatable :: psi2(:, :)
    type(xc_copied_potentials_t) :: vhxc_t2
 
    push_sub(td_etrs_sc)
 
    assert(hm%theory_level /= independent_particles)
 
    call messages_new_line()
    call messages_write('        Self-consistency iteration:')
    call messages_info(namespace=namespace)
 
    call td_etrs_predictor(ks, namespace, space, hm, ext_partners, gr, st, tr, time, dt, &
      ions_dyn, ions, mc)
 
    ! Store the state at half iteration
    safe_allocate_type_array(wfs_elec_t, psi2, (st%group%block_start:st%group%block_end, st%d%kpt%start:st%d%kpt%end))
    do ik = st%d%kpt%start, st%d%kpt%end
      do ib = st%group%block_start, st%group%block_end
        call st%group%psib(ib, ik)%copy_to(psi2(ib, ik), copy_data=.true.)
      end do
    end do
 
    orbital_dependent = v_ks_is_orbital_dependent(ks) .or. hm%lda_u_level /= dft_u_none
 
    do iter = 1, niter
      if (.not. orbital_dependent) call hm%ks_pot%store_potentials(vhxc_t2)
 
      ! Predict the states and density at the full step, n+1
      call propagation_ops_elec_fuse_density_exp_apply(tr%te, namespace, st, gr, hm, m_half * dt)
 
      if (orbital_dependent) then
        ! Change of the full-step states since the previous iterate; tr%st holds the predictor at iter 1
        diff = m_two * states_elec_calc_orbital_difference(gr, st, st%group%psib, tr%st%group%psib) / abs(dt)
        write(message(1), '(a,i0,a,es9.2)') '          step ', iter, ', residue = ', abs(diff)
        call messages_info(1, namespace=namespace)
        if (abs(diff) <= sctol) exit
        ! Keep this iterate as the reference for the next comparison
        do ik = st%d%kpt%start, st%d%kpt%end
          do ib = st%group%block_start, st%group%block_end
            call st%group%psib(ib, ik)%copy_data_to(gr%np, tr%st%group%psib(ib, ik))
          end do
        end do
      endif
 
      ! Rebuild potentials and H at n+1
      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)
 
      if (.not. orbital_dependent) then
        ! Check how much the potential changed
        diff = hm%ks_pot%check_convergence(vhxc_t2, gr, st%rho, st%qtot)
        write(message(1), '(a,i0,a,es9.2)') '          step ', iter, ', residue = ', abs(diff)
        call messages_info(1, namespace=namespace)
        if (abs(diff) <= sctol) exit
      endif
 
      if (iter /= niter) then
        ! we are not converged, restore the states at half iteration
        do ik = st%d%kpt%start, st%d%kpt%end
          do ib = st%group%block_start, st%group%block_end
            call psi2(ib, ik)%copy_data_to(gr%np, st%group%psib(ib, ik))
          end do
        end do
      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
 
    ! print an empty line
    call messages_info(namespace=namespace)
 
    if (present(scsteps)) scsteps = iter
 
    do ik = st%d%kpt%start, st%d%kpt%end
      do ib = st%group%block_start, st%group%block_end
        call psi2(ib, ik)%end()
      end do
    end do
 
    safe_deallocate_a(psi2)
 
    pop_sub(td_etrs_sc)
  end subroutine td_etrs_sc
 
  subroutine td_aetrs(namespace, space, hm, gr, st, tr, time, dt, ions_dyn, ions, ext_partners, mc)
    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
 
    type(gauge_field_t), pointer :: gfield
 
    push_sub(td_aetrs)
 
    ! propagate half of the time step with H(time - dt)
    call propagation_ops_elec_exp_apply(tr%te, namespace, st, gr, hm, m_half*dt)
 
    !Get the potentials from the interpolation
    call propagation_ops_elec_interpolate_get(hm, tr%vks_old)
 
    ! move the ions to time t
    call propagation_ops_elec_move_ions(tr%propagation_ops_elec, gr, hm, st, namespace, space, ions_dyn, &
      ions, ext_partners, mc, time, dt)
 
    !Propagate gauge field
    gfield => list_get_gauge_field(ext_partners)
    if(associated(gfield)) then
      call propagation_ops_elec_propagate_gauge_field(tr%propagation_ops_elec, gfield, dt, time)
    end if
 
    !Update Hamiltonian
    call propagation_ops_elec_update_hamiltonian(namespace, space, st, gr, hm, ext_partners, time)
 
    !Do the time propagation for the second half of the time step
    call propagation_ops_elec_fuse_density_exp_apply(tr%te, namespace, st, gr, hm, m_half*dt)
 
    pop_sub(td_aetrs)
  end subroutine td_aetrs
 
  ! ---------------------------------------------------------
  subroutine td_aetrs_sc(ks, namespace, space, hm, ext_partners, gr, st, tr, time, dt, &
    ions_dyn, ions, mc, 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(partner_list_t),             intent(in)    :: ext_partners
    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(multicomm_t),                intent(inout) :: mc
    real(real64),                     intent(in)    :: sctol
    integer,                optional, intent(out)   :: scsteps
 
    real(real64) :: diff
    integer :: iter, ib, ik
    type(states_elec_t) :: st0
    integer, parameter :: niter = 10
    logical :: orbital_dependent
    type(gauge_field_t), pointer :: gfield
    type(xc_copied_potentials_t) :: vhxc_t2
 
    push_sub(td_aetrs_sc)
 
    assert(hm%theory_level /= independent_particles)
 
    call messages_new_line()
    call messages_write('        Self-consistency iteration:')
    call messages_info(namespace=namespace)
 
    ! First half-step with H(t - dt)
    call propagation_ops_elec_exp_apply(tr%te, namespace, st, gr, hm, m_half*dt)
 
    ! Initial guess for H(t) from extrapolated potentials
    call propagation_ops_elec_interpolate_get(hm, tr%vks_old)
 
    ! Move ions and gauge to t only once before SC loop
    call propagation_ops_elec_move_ions(tr%propagation_ops_elec, gr, hm, st, namespace, space, ions_dyn, ions, &
      ext_partners, mc, time, dt)
 
    gfield => list_get_gauge_field(ext_partners)
    if (associated(gfield)) then
      call propagation_ops_elec_propagate_gauge_field(tr%propagation_ops_elec, gfield, dt, time)
    end if
 
    call propagation_ops_elec_update_hamiltonian(namespace, space, st, gr, hm, ext_partners, time)
 
    orbital_dependent = v_ks_is_orbital_dependent(ks) .or. hm%lda_u_level /= dft_u_none
 
    ! Store the half-step state to restart each SC iteration
    call states_elec_copy(st0, st)
 
    ! Reference states for the orbital residue, constructed on the first
    ! self-consistent step and reused for the rest of the run
    if (orbital_dependent .and. .not. allocated(tr%st)) then
      safe_allocate(tr%st)
      call states_elec_copy(tr%st, st)
    end if
 
    diff = m_huge
 
    do iter = 1, niter
      if (.not. orbital_dependent) call hm%ks_pot%store_potentials(vhxc_t2)
 
      call propagation_ops_elec_fuse_density_exp_apply(tr%te, namespace, st, gr, hm, m_half*dt)
 
      if (orbital_dependent) then
        call messages_write('          step ')
        call messages_write(iter)
        ! Nothing to evaluate for the first step, as there is no extrapolation for wave functions
        if (iter > 1) then
          ! Change of the full-step states since the previous iteration
          diff = m_two * states_elec_calc_orbital_difference(gr, st, st%group%psib, tr%st%group%psib) / abs(dt)
          call messages_write(', residue = ')
          call messages_write(abs(diff), fmt = '(1x,es9.2)')
        endif
        call messages_info(namespace=namespace)
        if (abs(diff) <= sctol) exit
        ! Keep this iterate as the reference for the next comparison
        do ik = st%d%kpt%start, st%d%kpt%end
          do ib = st%group%block_start, st%group%block_end
            call st%group%psib(ib, ik)%copy_data_to(gr%np, tr%st%group%psib(ib, ik))
          end do
        end do
      endif
 
      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)
 
      if (.not. orbital_dependent) then
        diff = hm%ks_pot%check_convergence(vhxc_t2, gr, st%rho, st%qtot)
        write(message(1), '(a,i0,a,es9.2)') '          step ', iter, ', residue = ', abs(diff)
        call messages_info(1, namespace=namespace)
        if (diff <= sctol) exit
      endif
 
      if (iter /= 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
 
    call states_elec_end(st0)
 
    pop_sub(td_aetrs_sc)
  end subroutine td_aetrs_sc
 
  ! ---------------------------------------------------------
  subroutine td_caetrs(ks, namespace, space, hm, ext_partners, gr, st, tr, time, dt, &
    ions_dyn, ions, 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),         intent(inout) :: hm
    type(partner_list_t),             intent(in)    :: ext_partners
    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(multicomm_t),                intent(inout) :: mc
 
    integer :: ik, ispin, ip, ib
    real(real64) :: vv
    complex(real64), allocatable :: phase(:)
    type(density_calc_t)  :: dens_calc
    integer(int64) :: pnp
    integer(int64), dimension(3) :: gsizes, bsizes
    type(accel_mem_t) :: phase_buff
    type(gauge_field_t), pointer :: gfield
    type(xc_copied_potentials_t) :: vold
 
    push_sub(td_caetrs)
 
    ! Get the interpolated KS potentials into vold
    call hm%ks_pot%get_interpolated_potentials(tr%vks_old, 2, storage=vold)
    ! And set it to the Hamiltonian
    call hm%ks_pot%restore_potentials(vold)
 
    call propagation_ops_elec_update_hamiltonian(namespace, space, st, gr, hm, ext_partners, time - dt)
 
    call v_ks_calc_start(ks, namespace, space, hm, st, ions, ions%latt, ext_partners, &
      time = time - dt, calc_energy = .false.)
 
    ! propagate half of the time step with H(time - dt)
    call propagation_ops_elec_exp_apply(tr%te, namespace, st, gr, hm, m_half*dt)
 
    call v_ks_calc_finish(ks, hm, namespace, space, ions%latt, st, ext_partners)
 
    call hm%ks_pot%set_interpolated_potentials(tr%vks_old, 1)
 
    call hm%ks_pot%perform_interpolation(tr%vks_old, (/time - dt, time - m_two*dt, time - m_three*dt/), time)
 
    ! Replace vold by 0.5(vhxc+vold)
    call hm%ks_pot%mix_potentials(vold, dt)
 
    ! copy vold to a cl buffer
    if (accel_is_enabled() .and. hm%apply_packed()) then
      if (family_is_mgga_with_exc(hm%xc)) then
        call messages_not_implemented('CAETRS propagator with accel and MGGA with energy functionals', namespace=namespace)
      end if
      pnp = accel_padded_size(gr%np)
      call accel_create_buffer(phase_buff, accel_mem_read_only, type_float, pnp*st%d%nspin)
      call vold%copy_vhxc_to_buffer(int(gr%np, int64), st%d%nspin, pnp, phase_buff)
    end if
 
    !Get the potentials from the interpolator
    call propagation_ops_elec_interpolate_get(hm, tr%vks_old)
 
    ! move the ions to time t
    call propagation_ops_elec_move_ions(tr%propagation_ops_elec, gr, hm, st, namespace, space, ions_dyn, &
      ions, ext_partners, mc, time, dt)
 
    gfield => list_get_gauge_field(ext_partners)
    if(associated(gfield)) then
      call propagation_ops_elec_propagate_gauge_field(tr%propagation_ops_elec, gfield, dt, time)
    end if
 
    call propagation_ops_elec_update_hamiltonian(namespace, space, st, gr, hm, ext_partners, time)
 
    call density_calc_init(dens_calc, st, gr, st%rho)
 
    ! propagate the other half with H(t)
    do ik = st%d%kpt%start, st%d%kpt%end
      ispin = st%d%get_spin_index(ik)
 
      do ib = st%group%block_start, st%group%block_end
        if (hm%apply_packed()) then
          call st%group%psib(ib, ik)%do_pack()
          if (hamiltonian_elec_inh_term(hm)) call hm%inh_st%group%psib(ib, ik)%do_pack()
        end if
 
        call profiling_in("CAETRS_PHASE")
        select case (st%group%psib(ib, ik)%status())
        case (batch_not_packed, batch_packed)
          safe_allocate(phase(1:gr%np))
          do ip = 1, gr%np
            vv = vold%vhxc(ip, ispin)
            phase(ip) = cmplx(cos(vv), -sin(vv), real64)
          end do
          call batch_mul_mf(gr%np, phase, st%group%psib(ib, ik), st%group%psib(ib, ik))
          safe_deallocate_a(phase)
        case (batch_device_packed)
          call accel_set_kernel_arg(kernel_phase, 0, pnp*(ispin - 1))
          call accel_set_kernel_arg(kernel_phase, 1, gr%np)
          call accel_set_kernel_arg(kernel_phase, 2, phase_buff)
          call accel_set_kernel_arg(kernel_phase, 3, st%group%psib(ib, ik)%ff_device)
          call accel_set_kernel_arg(kernel_phase, 4, log2(st%group%psib(ib, ik)%pack_size(1)))
 
          ! Compute the grid (extend to another dimensions if the size of the problem is too big)
          call accel_grid_size_extend_dim(pnp, st%group%psib(ib, ik)%pack_size(1), gsizes, bsizes)
 
          call accel_kernel_run(kernel_phase, gsizes, bsizes)
 
        end select
        call profiling_out("CAETRS_PHASE")
 
        call hm%phase%set_phase_corr(gr, st%group%psib(ib, ik))
        if (hamiltonian_elec_inh_term(hm)) then
          call tr%te%apply_batch(namespace, gr, hm, st%group%psib(ib, ik), m_half*dt, &
            inh_psib = hm%inh_st%group%psib(ib, ik))
        else
          call tr%te%apply_batch(namespace, gr, hm, st%group%psib(ib, ik), m_half*dt)
        end if
        call hm%phase%unset_phase_corr(gr, st%group%psib(ib, ik))
 
        call density_calc_accumulate(dens_calc, st%group%psib(ib, ik))
 
        if (hm%apply_packed()) then
          call st%group%psib(ib, ik)%do_unpack()
          if (hamiltonian_elec_inh_term(hm)) call hm%inh_st%group%psib(ib, ik)%do_unpack()
        end if
      end do
    end do
 
    call density_calc_end(dens_calc)
 
    if (accel_is_enabled() .and. hm%apply_packed()) then
      call accel_free_buffer(phase_buff)
    end if
 
    pop_sub(td_caetrs)
  end subroutine td_caetrs
 
end module propagator_etrs_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: