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 lasers_oct_m
  use messages_oct_m
  use mesh_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 wfs_elec_oct_m
  use xc_oct_m
 
  implicit none
 
  private
 
  public ::                    &
    td_magnus,                 &
    td_cfmagnus4
 
  
  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
 
 
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
 
  ! ---------------------------------------------------------
  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 (ions_dyn%is_active() .or. list_has_gauge_field(ext_partners)) then
      message(1) = "The commutator-free Magnus expansion cannot be used with moving ions, cell dynamics, 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%vks_old, 4, time, dt, t1, vhxc1)
    call potential_interpolation_interpolate(tr%vks_old, 4, time, dt, t2, vhxc2)
 
    hm%ks_pot%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%ks_pot%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: