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 batch_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_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 messages_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 v_ks_oct_m
  use xc_oct_m
 
  implicit none
 
  private
 
  public ::                    &
    td_magnus,                 &
    td_cfmagnus4
 
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(:, :, :), pot(:)
    type(lasers_t), pointer :: lasers
 
    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))
 
    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%vksold, 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
 
    tr%vmagnus(:, :, 2)  = m_half*(vaux(:, :, 1) + vaux(:, :, 2))
    tr%vmagnus(:, :, 1) = (sqrt(m_three)/12.0_real64)*dt*(vaux(:, :, 2) - vaux(:, :, 1))
 
    call propagation_ops_elec_fuse_density_exp_apply(tr%te, namespace, st, gr, hm, dt, vmagnus = tr%vmagnus)
 
    safe_deallocate_a(vaux)
    pop_sub(propagator_dt.td_magnus)
  end subroutine td_magnus
 
 
  ! ---------------------------------------------------------
  subroutine td_cfmagnus4(ks, namespace, space, hm, gr, st, tr, time, dt, ions_dyn, ions, ext_partners, 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(in)    :: 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
    integer,                          intent(in)    :: iter
 
    real(real64) :: alpha1, alpha2, c1, c2, t1, t2
    real(real64), allocatable :: vhxc1(:, :), vhxc2(:, :)
 
    if (ion_dynamics_ions_move(ions_dyn) .or. list_has_gauge_field(ext_partners)) then
      message(1) = "The commutator-free Magnus expansion cannot be used with moving ions or gauge fields"
      call messages_fatal(1, namespace=namespace)
    end if
 
    push_sub(propagator_dt.td_cfmagnus4)
 
    if (iter < 4) then
      call td_explicit_runge_kutta4(ks, namespace, space, hm, gr, st, time, dt, ions_dyn, ions, ext_partners)
      pop_sub(propagator_dt.td_cfmagnus4)
      return
    end if
 
 
    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
 
    safe_allocate(vhxc1(1:gr%np, 1:st%d%nspin))
    safe_allocate(vhxc2(1:gr%np, 1:st%d%nspin))
 
    call potential_interpolation_interpolate(tr%vksold, 4, time, dt, t1, vhxc1)
    call potential_interpolation_interpolate(tr%vksold, 4, time, dt, t2, vhxc2)
 
    hm%vhxc = m_two * (alpha2 * vhxc1 + alpha1 * vhxc2)
    call hamiltonian_elec_update_with_ext_pot(hm, gr, space, ext_partners, (/ t1, t2 /), (/ m_two * alpha2, m_two * alpha1 /))
    ! propagate by dt/2
    call propagation_ops_elec_exp_apply(tr%te, namespace, st, gr, hm, m_half*dt)
 
    hm%vhxc = m_two * (alpha1 * vhxc1 + alpha2 * vhxc2)
    call hamiltonian_elec_update_with_ext_pot(hm, gr, space, ext_partners, (/ t1, t2 /), (/ m_two * alpha1, m_two * alpha2 /))
    ! propagate by dt/2
    !TODO: fuse this with density calc
    call propagation_ops_elec_exp_apply(tr%te, namespace, st, gr, hm, m_half * dt)
 
    call density_calc(st, gr, st%rho)
 
    safe_deallocate_a(vhxc1)
    safe_deallocate_a(vhxc2)
    pop_sub(propagator_dt.td_cfmagnus4)
  end subroutine td_cfmagnus4
 
end module propagator_magnus_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: