mxll_elec_coupling.F90

Go to the documentation of this file.

!! Copyright (C) 2023 F. Bonafé
!!
!! 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 mxll_elec_coupling_oct_m
  use batch_ops_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_base_oct_m
  use, intrinsic :: iso_fortran_env
  use math_oct_m
  use mesh_oct_m
  use messages_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use space_oct_m
  use states_elec_dim_oct_m
  use wfs_elec_oct_m
 
  implicit none
 
  private
  public ::                          &
    mxll_coupling_t,                 &
    mxll_coupling_init,              &
    mxll_coupling_calc,              &
    mxll_coupling_end,               &
    mxll_coupling_copy,              &
    set_electric_quadrupole_pot,     &
    magnetic_dipole_coupling
 
  type :: mxll_coupling_t
    integer :: coupling_mode
    integer :: coupling_terms
    integer :: dipole_field
    real(real64), allocatable :: e_field(:,:)
    real(real64), allocatable :: b_field(:,:)
    real(real64), allocatable :: vec_pot(:,:)
    real(real64), allocatable :: e_field_dip(:)
    real(real64), allocatable :: vec_pot_dip(:)
    real(real64), allocatable :: b_field_dip(:)
    real(real64), allocatable :: e_quadrupole_pot(:)
    logical            :: calc_field_dip = .false.   
    logical            :: add_electric_dip = .false.   
    logical            :: add_electric_quad = .false.   
    logical            :: add_magnetic_dip = .false.   
    logical            :: add_zeeman = .false.   
    logical            :: add_non_uniform_vec_pot = .false. 
    real(real64) :: center_of_mass(1:3)
    integer :: center_of_mass_ip
    integer :: center_of_mass_rankmin
    logical :: test_equad = .false.
    type(derivatives_t), pointer, private :: der
    real(real64), private :: mass
  end type mxll_coupling_t
 
  ! See explanation in the MaxwellCouplingMode variable description.
  integer, public, parameter :: &
    NO_MAXWELL_COUPLING      = 0,  &
    length_gauge_dipole      = 1,  &
    multipolar_expansion     = 2,  &
    velocity_gauge_dipole    = 3,  &
    full_minimal_coupling    = 4
 
  integer, public, parameter :: &
    DIPOLE_AVERAGE    = 0,  &
    dipole_at_com     = 1
 
 
contains
 
  ! ------------------------------------------------------------------------
  subroutine mxll_coupling_init(this, d, gr, namespace, mass)
    type(mxll_coupling_t),  intent(inout)  :: this
    type(states_elec_dim_t),intent(in)     :: d
    type(grid_t),   target, intent(in)     :: gr
    type(namespace_t),      intent(in)     :: namespace
    real(real64),           intent(in)     :: mass
 
    integer :: terms_default, fmc_default, multipolar_terms, pauli_terms
 
    push_sub(mxll_coupling_init)
 
    this%der => gr%der
    this%mass = mass
 
    !%Variable MaxwellCouplingMode
    !%Type integer
    !%Default none
    !%Section Hamiltonian
    !%Description
    !% This variable selects the level of light-matter coupling with electromagnetic fields from the Maxwell solver.
    !% Each option defines a coupling Hamiltonian <math> H_{int}</math>.
    !%Option none 0
    !% No coupling.
    !%Option length_gauge_dipole 1
    !% Dipole level in length gauge with transverse electric field E(t):
    !% <math> H_{int} = -r.E </math>
    !% The option MaxwellDipoleField should be set to define if <math>E</math> is calculated as an average at evaluated
    !% at the center of mass (by default, dipole coupling with the average electric field will be used).
    !%Option multipolar_expansion 2
    !% The option MultipoleExpansionTerms selects the terms to be included (electric dipole, electric quadrupole
    !% and/or magnetic dipole).
    !%Option velocity_gauge_dipole 3
    !% Coupling at dipole level in velocity gauge with transverse vector potential A(t).
    !% <math> H_{int} = 1/2m (2 A(t).p + A^2(t)) </math>
    !% The option MaxwellDipoleField should be set.
    !%Option full_minimal_coupling 4
    !% Full vector potential A(r,t) coupled to the momentum.
    !% <math> H_{int} = 1/2m (2 A(r,t).p + A^2(r,t)) </math>
    !%End
    call parse_variable(namespace, 'MaxwellCouplingMode', no_maxwell_coupling, &
      this%coupling_mode)
    this%add_electric_dip = (this%coupling_mode == length_gauge_dipole)
 
    if (this%coupling_mode == multipolar_expansion) then
 
      !%Variable MultipolarExpansionTerms
      !%Type flag
      !%Default electric_dipole + electric_quadrupole + magnetic_dipole
      !%Section Hamiltonian
      !%Description
      !% Terms to be included in the multipolar expansion.
      !% For this type of coupling to make sense, the E field has to be calculated at the center of mass (not averaged).
      !%Option electric_dipole bit(1)
      !% Adds electric dipole term in length gauge, uses transverse electric field E(t):
      !% <math> H_{int} = -e (r.E) </math>
      !%Option electric_quadrupole bit(2)
      !% Adds electric quadrupole term:
      !% <math> H_{int} = (1/2) e (r . Q . r) </math>
      !% where Q is the outer product of gradient and electric field:
      !% <math> Q_{ij} = partial_i E_j |_{r=r_0} </math>
      !%Option magnetic_dipole bit(3)
      !% Adds magnetic dipole term:
      !% <math> H_{int} = -i (e hbar / 2m) sum_i (B(r_0) . (r x grad)) </math>
      !%End
      terms_default = &
        option__multipolarexpansionterms__electric_dipole + &
        option__multipolarexpansionterms__electric_quadrupole + &
        option__multipolarexpansionterms__magnetic_dipole
 
      call parse_variable(namespace, 'MultipolarExpansionTerms', terms_default, multipolar_terms)
 
      if (bitand(multipolar_terms, option__multipolarexpansionterms__electric_dipole) /= 0) then
        this%add_electric_dip = .true.
      end if
      if (bitand(multipolar_terms, option__multipolarexpansionterms__electric_quadrupole) /=0) then
        this%add_electric_quad = .true.
      end if
      if (bitand(multipolar_terms, option__multipolarexpansionterms__magnetic_dipole) /=0) then
        this%add_magnetic_dip = .true.
      end if
 
    end if
 
    if (this%coupling_mode == full_minimal_coupling) then
 
      !%Variable PauliHamiltonianTerms
      !%Type flag
      !%Default non_uniform_vector_potential + zeeman_term
      !%Section Hamiltonian
      !%Description
      !% Terms to be included in the Pauli Hamiltonnian.
      !%Option non_uniform_vector_potential bit(1)
      !% Adds non-uniform vector potential to the canonical momentum.
      !%Option zeeman_term bit(2)
      !% Adds the non-uniform Zeeman potential.
      !%End
 
      if (d%nspin == 2) then
        fmc_default = option__paulihamiltonianterms__non_uniform_vector_potential + &
          option__paulihamiltonianterms__zeeman_term
      else
        fmc_default = option__paulihamiltonianterms__non_uniform_vector_potential
      end if
 
      call parse_variable(namespace, 'PauliHamiltonianTerms', fmc_default, pauli_terms)
      if (bitand(pauli_terms, option__paulihamiltonianterms__non_uniform_vector_potential) /= 0) then
        this%add_non_uniform_vec_pot = .true.
      end if
      if (bitand(pauli_terms, option__paulihamiltonianterms__zeeman_term) /=0) then
        this%add_zeeman = .true.
      end if
 
    end if
 
    !%Variable MaxwellDipoleField
    !%Type integer
    !%Default average
    !%Section Hamiltonian
    !%Description
    !% This variable selects the method to get the E field vector at the position of the system
    !% if Maxwell-matter coupling at dipole level within length gauge is done.
    !%Option average 0
    !% Transverse E field is averaged in the volume occupied by the matter system.
    !%Option center_of_mass 1
    !% Tranverse E field is evaluated at the center of mass of the matter system (at a single point).
    !%End
    call parse_variable(namespace, 'MaxwellDipoleField', dipole_average, &
      this%dipole_field)
 
    if (this%coupling_mode /= no_maxwell_coupling) then
      safe_allocate(this%e_field(1:gr%np, 1:gr%box%dim))
      safe_allocate(this%b_field(1:gr%np, 1:gr%box%dim))
      safe_allocate(this%vec_pot(1:gr%np, 1:gr%box%dim))
      this%e_field = m_zero
      this%b_field = m_zero
      this%vec_pot = m_zero
      safe_allocate(this%e_field_dip(1:gr%box%dim))
      this%e_field_dip = m_zero
      safe_allocate(this%b_field_dip(1:gr%box%dim))
      this%b_field_dip = m_zero
      safe_allocate(this%vec_pot_dip(1:gr%box%dim))
      this%vec_pot_dip = m_zero
    end if
 
    if (this%coupling_mode == multipolar_expansion .and. this%add_electric_quad) then
      safe_allocate(this%e_quadrupole_pot(1:gr%np))
    end if
 
    call mxll_quadrupole_test_init(this, gr, namespace)
 
    pop_sub(mxll_coupling_init)
 
  end subroutine mxll_coupling_init
 
  ! ------------------------------------------------------------------------
  subroutine mxll_quadrupole_test_init(this, gr, namespace)
    type(mxll_coupling_t),  intent(inout)  :: this
    type(grid_t),           intent(in)     :: gr
    type(namespace_t),      intent(in)     :: namespace
 
    integer :: ip, rankmin
    real(real64) :: dmin
 
    push_sub(mxll_quadrupole_test_init)
 
    !%Variable MaxwellTestQuadrupole
    !%Type logical
    !%Default no
    !%Section Maxwell
    !%Description
    !% Override Maxwell field with linear E field for testing.
    !%End
    call parse_variable(namespace, 'MaxwellTestQuadrupole', .false., this%test_equad)
 
    if (this%test_equad) then
      safe_allocate(this%e_field(1:gr%np, 1:gr%box%dim))
      this%e_field = m_zero
      do ip = 1, gr%np
        this%e_field(ip,1) = 0.02_real64 * gr%x(1, ip)
      end do
      this%calc_field_dip = .true.
      this%center_of_mass(1:3) = m_zero
      this%center_of_mass_ip = mesh_nearest_point(gr, this%center_of_mass, dmin, rankmin)
      this%center_of_mass_rankmin = rankmin
      safe_allocate(this%e_quadrupole_pot(1:gr%np))
      this%e_quadrupole_pot = m_zero
    end if
 
    pop_sub(mxll_quadrupole_test_init)
 
  end subroutine mxll_quadrupole_test_init
 
  ! ------------------------------------------------------------------------
  subroutine mxll_coupling_calc(this, hm_base, mesh, d, space)
    type(mxll_coupling_t),         intent(inout) :: this
    type(hamiltonian_elec_base_t), intent(inout) :: hm_base
    type(mesh_t),                  intent(in)    :: mesh
    type(states_elec_dim_t),       intent(in)    :: d
    class(space_t),                intent(in)    :: space
 
    integer :: ip, ispin, idir
 
    push_sub(mxll_coupling_calc)
 
    ! coupling to Maxwell field
    select case (this%coupling_mode)
 
    case (length_gauge_dipole, multipolar_expansion)
 
      if (this%add_electric_dip) then
        do ispin = 1, d%spin_channels
          do ip = 1, mesh%np
            ! The -1 sign is missing here. Check epot.F90 for the explanation.
            hm_base%potential(ip, ispin) = hm_base%potential(ip, ispin) + &
              sum(this%e_field_dip(1:space%dim) * (mesh%x(1:space%dim, ip) - this%center_of_mass(1:space%dim)))
          end do
        end do
      end if
 
      if (this%add_electric_quad .and. this%calc_field_dip) then
        ! the calc_field_dip condition prevents calling the routine before initializing the interactions
        call set_electric_quadrupole_pot(this, mesh)
        do ispin = 1, d%spin_channels
          do ip = 1, mesh%np
            hm_base%potential(ip, ispin) = hm_base%potential(ip, ispin) + this%e_quadrupole_pot(ip)
          end do
        end do
      end if
 
      ! magnetic_dipole term is added in hamiltonian_elec_inc.F90
 
    case (velocity_gauge_dipole)
      call hm_base%allocate_field(mesh, field_uniform_vector_potential, .false.)
      ! The minus sign is here is for the wrong convention of Octopus for the phase
      hm_base%uniform_vector_potential(1:space%dim) = hm_base%uniform_vector_potential(1:space%dim) - &
        this%vec_pot_dip(1:space%dim)
 
    case (full_minimal_coupling)
      call hm_base%allocate_field(mesh, field_vector_potential, .false.)
      if (this%add_non_uniform_vec_pot) then
        do idir = 1, mesh%box%dim
          hm_base%vector_potential(idir, :) = hm_base%vector_potential(idir, :) + this%vec_pot(:, idir)
        end do
      end if
      if (this%add_zeeman) then
        call hm_base%allocate_field(mesh, field_magnetic_field, .false.)
        hm_base%magnetic_field(:,:) = hm_base%magnetic_field(:,:) + this%b_field(:,:)
      end if
 
    end select
 
    pop_sub(mxll_coupling_calc)
 
  end subroutine mxll_coupling_calc
 
  ! ------------------------------------------------------------------------
  subroutine mxll_coupling_end(this)
    type(mxll_coupling_t), intent(inout) :: this
 
    push_sub(mxll_coupling_end)
 
    safe_deallocate_a(this%e_field)
    safe_deallocate_a(this%b_field)
    safe_deallocate_a(this%vec_pot)
    safe_deallocate_a(this%e_field_dip)
    safe_deallocate_a(this%b_field_dip)
    safe_deallocate_a(this%vec_pot_dip)
    safe_deallocate_a(this%e_quadrupole_pot)
 
    pop_sub(mxll_coupling_end)
  end subroutine mxll_coupling_end
 
  ! ------------------------------------------------------------------------
  subroutine mxll_coupling_copy(this_out, this_in, der_target)
    type(mxll_coupling_t),             intent(out) :: this_out
    type(mxll_coupling_t),             intent(in)  :: this_in
    type(derivatives_t), target,       intent(in)  :: der_target
 
    push_sub(mxll_coupling_copy)
 
    this_out%coupling_mode = this_in%coupling_mode
    this_out%coupling_terms = this_in%coupling_terms
    this_out%dipole_field = this_in%dipole_field
    this_out%calc_field_dip = this_in%calc_field_dip
    this_out%add_electric_dip = this_in%add_electric_dip
    this_out%add_electric_quad = this_in%add_electric_quad
    this_out%add_magnetic_dip = this_in%add_magnetic_dip
    this_out%add_zeeman = this_in%add_zeeman
    this_out%add_non_uniform_vec_pot = this_in%add_non_uniform_vec_pot
    this_out%center_of_mass = this_in%center_of_mass
    this_out%center_of_mass_ip = this_in%center_of_mass_ip
    this_out%center_of_mass_rankmin = this_in%center_of_mass_rankmin
    this_out%test_equad = this_in%test_equad
 
    this_out%der => der_target
    this_out%mass = this_in%mass
 
    if (allocated(this_in%e_field)) then
      safe_allocate(this_out%e_field(1:size(this_in%e_field, dim=1), 1:size(this_in%e_field, dim=2)))
      this_out%e_field = this_in%e_field
    end if
 
    if (allocated(this_in%b_field)) then
      safe_allocate(this_out%b_field(1:size(this_in%b_field, dim=1), 1:size(this_in%b_field, dim=2)))
      this_out%b_field = this_in%b_field
    end if
 
    if (allocated(this_in%vec_pot)) then
      safe_allocate(this_out%vec_pot(1:size(this_in%vec_pot, dim=1), 1:size(this_in%vec_pot, dim=2)))
      this_out%vec_pot = this_in%vec_pot
    end if
 
    if (allocated(this_in%e_field_dip)) then
      safe_allocate(this_out%e_field_dip(1:size(this_in%e_field_dip)))
      this_out%e_field_dip = this_in%e_field_dip
    end if
 
    if (allocated(this_in%vec_pot_dip)) then
      safe_allocate(this_out%vec_pot_dip(1:size(this_in%vec_pot_dip)))
      this_out%vec_pot_dip = this_in%vec_pot_dip
    end if
 
    if (allocated(this_in%b_field_dip)) then
      safe_allocate(this_out%b_field_dip(1:size(this_in%b_field_dip)))
      this_out%b_field_dip = this_in%b_field_dip
    end if
 
    if (allocated(this_in%e_quadrupole_pot)) then
      safe_allocate(this_out%e_quadrupole_pot(1:size(this_in%e_quadrupole_pot)))
      this_out%e_quadrupole_pot = this_in%e_quadrupole_pot
    end if
 
    pop_sub(mxll_coupling_copy)
  end subroutine mxll_coupling_copy
 
 
  ! ------------------------------------------------------------------------
  subroutine set_electric_quadrupole_pot(this, mesh)
    type(mxll_coupling_t),    intent(inout)  :: this
    type(mesh_t),             intent(in)     :: mesh
 
    real(real64), allocatable   :: e_field_quadrupole_tensor(:,:), tmp_partial(:,:), this_e_field(:,:)
    real(real64)  :: r(3), tensor_dot_rr(3)
    integer :: ip, i, dims
 
    push_sub(set_electric_quadrupole_pot)
 
    call profiling_in('SET_ELECTRIC_QUADRUPOLE')
 
    safe_allocate(e_field_quadrupole_tensor(1:mesh%box%dim, 1:mesh%box%dim))
    safe_allocate(this_e_field(1:mesh%np_part, 1:mesh%box%dim))
    safe_allocate(tmp_partial(1:mesh%np, 1:mesh%box%dim))
 
    this_e_field(1:mesh%np,:) = this%e_field(1:mesh%np,:)
 
    e_field_quadrupole_tensor = m_zero
    dims = this%der%dim
    do i = 1, dims
      call dderivatives_grad(this%der, this_e_field(:, i), tmp_partial)
      if (mesh%mpi_grp%rank == this%center_of_mass_rankmin) then
        e_field_quadrupole_tensor(i, 1:dims) = tmp_partial(this%center_of_mass_ip, 1:dims)
      end if
    end do
    call mesh%allreduce(e_field_quadrupole_tensor)
 
    do ip = 1, mesh%np
      r(:) = mesh%x(:, ip) - this%center_of_mass(:)
      tensor_dot_rr = matmul(e_field_quadrupole_tensor, r)
      ! This term is +1/2 r.Q.r and not -1/2 r.Q.r (as e=+1 in Octopus)
      this%e_quadrupole_pot(ip) = m_half * dot_product(r, tensor_dot_rr)
    end do
 
    safe_deallocate_a(e_field_quadrupole_tensor)
    safe_deallocate_a(tmp_partial)
    safe_deallocate_a(this_e_field)
 
    call profiling_out('SET_ELECTRIC_QUADRUPOLE')
    pop_sub(set_electric_quadrupole_pot)
 
  end subroutine set_electric_quadrupole_pot
 
 
  ! ------------------------------------------------------------------------
  subroutine magnetic_dipole_coupling(this, psib, hpsib)
    type(mxll_coupling_t),  intent(in)    :: this
    type(wfs_elec_t),       intent(inout) :: psib
    type(wfs_elec_t),       intent(inout) :: hpsib
 
    integer             :: ip, idir
    real(real64)        :: rtmp(3)
    real(real64), allocatable  :: bxr(:,:)
    type(wfs_elec_t) :: gradb(this%der%dim)
    type(wfs_elec_t) :: bxr_dot_gradb
 
    push_sub(magnetic_dipole_coupling)
 
    safe_allocate(bxr(1:this%der%mesh%np,1:3))
    call profiling_in('MAGNETIC_DIPOLE')
 
    do ip = 1, this%der%mesh%np
      rtmp(:) = this%der%mesh%x(:, ip) - this%center_of_mass
      bxr(ip,:) = dcross_product(this%b_field_dip, rtmp)
    end do
 
    call zderivatives_batch_grad(this%der, psib, gradb)
 
    call hpsib%copy_to(bxr_dot_gradb)
    do idir = 1, this%der%dim
      call batch_mul_mf(this%der%mesh%np, bxr(:,idir), gradb(idir), bxr_dot_gradb)
      call batch_axpy(this%der%mesh%np, m_zi*m_half/this%mass, bxr_dot_gradb, hpsib)
    end do
 
    call bxr_dot_gradb%end()
    do idir = 1, this%der%dim
      call gradb(idir)%end()
    end do
    safe_deallocate_a(bxr)
 
    call profiling_out('MAGNETIC_DIPOLE')
    pop_sub(magnetic_dipole_coupling)
  end subroutine magnetic_dipole_coupling
 
end module mxll_elec_coupling_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: