hamiltonian_elec_base.F90

Go to the documentation of this file.

!! Copyright (C) 2009 X. Andrade
!!
!! 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 hamiltonian_elec_base_oct_m
  use accel_oct_m
  use batch_oct_m
  use batch_ops_oct_m
  use blas_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use electron_space_oct_m
  use epot_oct_m
  use global_oct_m
  use hardware_oct_m
  use math_oct_m
  use mesh_oct_m
  use messages_oct_m
  use mpi_oct_m
  use nl_operator_oct_m
  use profiling_oct_m
  use ps_oct_m
  use space_oct_m
  use states_elec_oct_m
  use states_elec_dim_oct_m
  use types_oct_m
  use wfs_elec_oct_m
 
  implicit none
 
  private
 
  public ::                                         &
    hamiltonian_elec_base_t,                        &
    hamiltonian_elec_base_accel_rebuild,            &
    dhamiltonian_elec_base_local_sub,               &
    zhamiltonian_elec_base_local_sub
 
  type hamiltonian_elec_base_t
    private
    integer                                       :: nspin
    real(real64)                                  :: mass
    real(real64)                                  :: rashba_coupling
    type(nl_operator_t),      pointer,     public :: kinetic
    real(real64),             allocatable, public :: potential(:, :)
    real(real64),             allocatable, public :: Impotential(:, :)
    real(real64),             allocatable, public :: uniform_magnetic_field(:)
    real(real64),             allocatable, public :: magnetic_field(:, :)
    real(real64),             allocatable, public :: zeeman_pot(:, :)
    real(real64),             allocatable, public :: uniform_vector_potential(:)
    real(real64),             allocatable, public :: vector_potential(:, :)
    type(accel_mem_t)                             :: potential_accel
    type(accel_mem_t)                             :: impotential_accel
  contains
 
    procedure :: init => hamiltonian_elec_base_init
    
    procedure :: end => hamiltonian_elec_base_end
    
    procedure :: clear => hamiltonian_elec_base_clear
    
    procedure :: update_magnetic_terms => hamiltonian_elec_base_update_magnetic_terms
    
    procedure :: allocate_field => hamiltonian_elec_base_allocate
    
    procedure :: accel_copy_pot => hamiltonian_elec_base_accel_copy_pot
    
    procedure :: has_magnetic => hamiltonian_elec_base_has_magnetic
    
    procedure :: has_zeeman => hamiltonian_elec_base_has_zeeman
    
    procedure :: has_vector_potential => hamiltonian_elec_base_has_vector_potential
    
    procedure :: calc_rashba => hamiltonian_elec_base_rashba
    
    procedure :: dcalc_magnetic => dhamiltonian_elec_base_magnetic
    
    procedure :: zcalc_magnetic => zhamiltonian_elec_base_magnetic
    
    procedure :: dcalc_local => dhamiltonian_elec_base_local
    
    procedure :: zcalc_local => zhamiltonian_elec_base_local
    
  end type hamiltonian_elec_base_t
 
  integer, parameter, public ::          &
    TERM_ALL                 = huge(1),  &
    term_kinetic             =   1,      &
    term_local_potential     =   2,      &
    term_non_local_potential =   4,      &
    term_others              =   8,      &
    term_local_external      =  16,      &
    term_mgga                =  32,      &
    term_dft_u               =  64,      &
    term_rdmft_occ           = 128
 
  integer, parameter, public ::            &
    FIELD_POTENTIAL                = 1,    &
    field_vector_potential         = 2,    &
    field_uniform_vector_potential = 4,    &
    field_uniform_magnetic_field   = 8,    &
    field_magnetic_field           = 16
 
contains
 
  ! ---------------------------------------------------------
  subroutine hamiltonian_elec_base_init(this, nspin, mass, rashba_coupling)
    class(hamiltonian_elec_base_t), intent(inout) :: this
    integer,                  intent(in)    :: nspin
    real(real64),             intent(in)    :: mass
    real(real64),             intent(in)    :: rashba_coupling
 
    push_sub(hamiltonian_elec_base_init)
 
    this%nspin = nspin
    this%mass  = mass
    this%rashba_coupling = rashba_coupling
 
    pop_sub(hamiltonian_elec_base_init)
  end subroutine hamiltonian_elec_base_init
 
  ! ---------------------------------------------------------
  subroutine hamiltonian_elec_base_end(this)
    class(hamiltonian_elec_base_t), intent(inout) :: this
 
    push_sub(hamiltonian_elec_base_end)
 
    if (allocated(this%potential) .and. accel_is_enabled()) then
      call accel_free_buffer(this%potential_accel)
      if (allocated(this%Impotential)) then
        call accel_free_buffer(this%impotential_accel)
      end if
    end if
 
    safe_deallocate_a(this%potential)
    safe_deallocate_a(this%Impotential)
    safe_deallocate_a(this%vector_potential)
    safe_deallocate_a(this%uniform_vector_potential)
    safe_deallocate_a(this%uniform_magnetic_field)
    safe_deallocate_a(this%magnetic_field)
    safe_deallocate_a(this%zeeman_pot)
 
    pop_sub(hamiltonian_elec_base_end)
  end subroutine hamiltonian_elec_base_end
 
  ! ----------------------------------------------------------
  !
  !
  subroutine hamiltonian_elec_base_clear(this, np)
    class(hamiltonian_elec_base_t), intent(inout) :: this
    integer,                       intent(in)    :: np
 
    integer :: ip, ispin
 
    push_sub(hamiltonian_elec_clear)
 
    if (allocated(this%potential)) then
      do ispin = 1, this%nspin
        !$omp parallel do schedule(static)
        do ip = 1, np
          this%potential(ip, ispin) = m_zero
        end do
      end do
    end if
 
    if (allocated(this%Impotential)) then
      do ispin = 1, this%nspin
        !$omp parallel do schedule(static)
        do ip = 1, np
          this%Impotential(ip, ispin) = m_zero
        end do
      end do
    end if
 
    if (allocated(this%uniform_vector_potential)) then
      this%uniform_vector_potential = m_zero
    end if
 
    if (allocated(this%vector_potential)) then
      this%vector_potential = m_zero
    end if
 
    if (allocated(this%uniform_magnetic_field)) then
      this%uniform_magnetic_field = m_zero
    end if
 
    if (allocated(this%magnetic_field)) then
      this%magnetic_field = m_zero
    end if
 
    if (allocated(this%zeeman_pot)) then
      this%zeeman_pot = m_zero
    end if
 
    pop_sub(hamiltonian_elec_clear)
  end subroutine hamiltonian_elec_base_clear
 
 
  ! ---------------------------------------------------------------
  !
  subroutine hamiltonian_elec_base_allocate(this, mesh, field, complex_potential)
    class(hamiltonian_elec_base_t), intent(inout) :: this
    class(mesh_t),            intent(in)    :: mesh
    integer,                  intent(in)    :: field
    logical,                  intent(in)    :: complex_potential
 
    integer :: ip, ispin
 
    push_sub(hamiltonian_elec_base_allocate)
 
    if (bitand(field_potential, field) /= 0) then
      if (.not. allocated(this%potential)) then
        safe_allocate(this%potential(1:mesh%np, 1:this%nspin))
        do ispin = 1, this%nspin
          !$omp parallel do schedule(static)
          do ip = 1, mesh%np
            this%potential(ip, ispin) = m_zero
          end do
        end do
        if (complex_potential) then
          safe_allocate(this%Impotential(1:mesh%np, 1:this%nspin))
          do ispin = 1, this%nspin
            !$omp parallel do schedule(static)
            do ip = 1, mesh%np
              this%Impotential(ip, ispin) = m_zero
            end do
          end do
        end if
        if (accel_is_enabled()) then
          call accel_create_buffer(this%potential_accel, accel_mem_read_only, type_float, accel_padded_size(mesh%np)*this%nspin)
          if (complex_potential) then
            call accel_create_buffer(this%impotential_accel, accel_mem_read_only, type_float, &
              accel_padded_size(mesh%np)*this%nspin)
          end if
        end if
      end if
    end if
 
    if (bitand(field_uniform_vector_potential, field) /= 0) then
      if (.not. allocated(this%uniform_vector_potential)) then
        safe_allocate(this%uniform_vector_potential(1:mesh%box%dim))
        this%uniform_vector_potential = m_zero
      end if
    end if
 
    if (bitand(field_vector_potential, field) /= 0) then
      if (.not. allocated(this%vector_potential)) then
        safe_allocate(this%vector_potential(1:mesh%box%dim, 1:mesh%np))
        this%vector_potential = m_zero
      end if
    end if
 
    if (bitand(field_uniform_magnetic_field, field) /= 0) then
      if (.not. allocated(this%uniform_magnetic_field)) then
        safe_allocate(this%uniform_magnetic_field(1:max(mesh%box%dim, 3)))
        this%uniform_magnetic_field = m_zero
      end if
    end if
 
    if ((bitand(field_magnetic_field, field) /= 0) .or. &
      bitand(field_uniform_magnetic_field, field) /= 0) then
      if (.not. allocated(this%magnetic_field)) then
        safe_allocate(this%magnetic_field(1:mesh%np, 1:max(mesh%box%dim, 3)))
        this%magnetic_field = m_zero
        safe_allocate(this%zeeman_pot(1:mesh%np, 1:this%nspin))
        this%zeeman_pot = m_zero
      end if
    end if
 
    pop_sub(hamiltonian_elec_base_allocate)
  end subroutine hamiltonian_elec_base_allocate
 
  ! ----------------------------------------------------------
  !
  !
  subroutine hamiltonian_elec_base_update_magnetic_terms(this, mesh, gyromagnetic_ratio, ispin)
    class(hamiltonian_elec_base_t), intent(inout) :: this
    class(mesh_t),                  intent(in)    :: mesh
    real(real64),                   intent(in)    :: gyromagnetic_ratio
    integer,                        intent(in)    :: ispin
 
    real(real64) :: cc, b_perp2
    integer :: idir, ip
    logical, save :: warning_shown = .false.
 
    push_sub(hamiltonian_elec_base_update_magnetic_terms)
 
    if (allocated(this%uniform_vector_potential) .and. allocated(this%vector_potential)) then
      ! copy the uniform vector potential onto the non-uniform one
      do idir = 1, mesh%box%dim
        !$omp parallel do schedule(static)
        do ip = 1, mesh%np
          this%vector_potential(idir, ip) = &
            this%vector_potential(idir, ip) + this%uniform_vector_potential(idir)
        end do
      end do
      safe_deallocate_a(this%uniform_vector_potential)
    end if
 
 
    if (.not. allocated(this%magnetic_field)) then
      pop_sub(hamiltonian_elec_base_update_magnetic_terms)
      return
    end if
 
    ! add the Zeeman term
    ! It reads in the Pauli equation as  $\frac{|e|\hbar}{2m} \sigma\cdot\vec{B}$
    ! The gyromagnetic ratio is twice the magnetic moment of the electron in units of Bohr magneton
    ! This is why we multiply by $\frac{|g_e|}{2}$.
    ! The factor 1/c comes from the factor 1/c in front of the vector potential, which leads to
    ! the Zeeman term starting from the Dirac Hamiltonian.
    cc = m_half/p_c*gyromagnetic_ratio*m_half
 
    if (allocated(this%uniform_magnetic_field) ) then
      !$omp parallel private(ip)
      do idir = 1, 3
        !$omp do
        do ip = 1, mesh%np
          this%magnetic_field(ip, idir) = this%magnetic_field(ip, idir) + this%uniform_magnetic_field(idir)
        end do
        !$omp end do nowait
      end do
      !$omp end parallel
    end if
 
    select case (ispin)
    case (spin_polarized)
      ! The collinear spin axis is z, as for spinors: only B_z enters, with its sign
      b_perp2 = m_zero
      !$omp parallel do reduction(+:b_perp2)
      do ip = 1, mesh%np
        this%zeeman_pot(ip, 1) =   cc*this%magnetic_field(ip, 3)
        this%zeeman_pot(ip, 2) = - cc*this%magnetic_field(ip, 3)
        b_perp2 = b_perp2 + this%magnetic_field(ip, 1)**2 + this%magnetic_field(ip, 2)**2
      end do
      call mesh%allreduce(b_perp2)
      if (b_perp2 > m_epsilon .and. .not. warning_shown) then
        message(1) = "The magnetic field has x or y components, which do not enter the Zeeman term"
        message(2) = "for spin-polarized calculations (spin axis z). Use spinors to include them."
        call messages_warning(2)
        warning_shown = .true.
      end if
 
    case (spinors)
      !$omp parallel do
      do ip = 1, mesh%np
        this%zeeman_pot(ip, 1) =  cc*this%magnetic_field(ip, 3)
        this%zeeman_pot(ip, 2) = - cc*this%magnetic_field(ip, 3)
        this%zeeman_pot(ip, 3) =   cc*this%magnetic_field(ip, 1)
        this%zeeman_pot(ip, 4) = - cc*this%magnetic_field(ip, 2)
      end do
    end select
    safe_deallocate_a(this%uniform_magnetic_field)
 
    pop_sub(hamiltonian_elec_base_update_magnetic_terms)
  end subroutine hamiltonian_elec_base_update_magnetic_terms
 
 
  !--------------------------------------------------------
  !
  subroutine hamiltonian_elec_base_accel_copy_pot(this, mesh)
    class(hamiltonian_elec_base_t), intent(inout) :: this
    class(mesh_t),            intent(in)    :: mesh
 
    integer :: offset, ispin
 
    push_sub(hamiltonian_elec_base_accel_copy_pot)
 
    if (allocated(this%potential) .and. accel_is_enabled()) then
      offset = 0
      do ispin = 1, this%nspin
        call accel_write_buffer(this%potential_accel, mesh%np, this%potential(:, ispin), offset = offset)
        if(allocated(this%Impotential)) then
          call accel_write_buffer(this%impotential_accel, mesh%np, this%Impotential(:, ispin), offset = offset)
        end if
        offset = offset + accel_padded_size(mesh%np)
      end do
    end if
 
    pop_sub(hamiltonian_elec_base_accel_copy_pot)
  end subroutine hamiltonian_elec_base_accel_copy_pot
 
  ! ---------------------------------------------------------
  subroutine hamiltonian_elec_base_accel_rebuild(this, mesh)
    class(hamiltonian_elec_base_t), intent(inout) :: this
    class(mesh_t),                  intent(in)    :: mesh
 
    logical :: complex_potential
 
    push_sub(hamiltonian_elec_base_accel_rebuild)
 
    if (.not. accel_is_enabled()) then
      pop_sub(hamiltonian_elec_base_accel_rebuild)
      return
    end if
 
    call accel_detach_buffer(this%potential_accel)
    call accel_detach_buffer(this%impotential_accel)
 
    if (.not. allocated(this%potential)) then
      pop_sub(hamiltonian_elec_base_accel_rebuild)
      return
    end if
 
    complex_potential = allocated(this%Impotential)
 
    call accel_create_buffer(this%potential_accel, accel_mem_read_only, type_float, &
      accel_padded_size(mesh%np)*this%nspin)
    if (complex_potential) then
      call accel_create_buffer(this%impotential_accel, accel_mem_read_only, type_float, &
        accel_padded_size(mesh%np)*this%nspin)
    end if
 
    call this%accel_copy_pot(mesh)
 
    pop_sub(hamiltonian_elec_base_accel_rebuild)
  end subroutine hamiltonian_elec_base_accel_rebuild
 
  ! ----------------------------------------------------------------------------------
  !
  logical pure function hamiltonian_elec_base_has_magnetic(this) result(has_magnetic)
    class(hamiltonian_elec_base_t), intent(in) :: this
 
    has_magnetic = allocated(this%vector_potential) &
      .or. allocated(this%uniform_magnetic_field) &
      .or. allocated(this%magnetic_field)
 
  end function hamiltonian_elec_base_has_magnetic
 
  ! ----------------------------------------------------------------------------------
  !
  logical pure function hamiltonian_elec_base_has_zeeman(this) result(has_zeeman)
    class(hamiltonian_elec_base_t), intent(in) :: this
 
    has_zeeman = allocated(this%zeeman_pot)
 
  end function hamiltonian_elec_base_has_zeeman
 
  ! ----------------------------------------------------------------------------------
  !
  logical pure function hamiltonian_elec_base_has_vector_potential(this) result(has_vector_potential)
    class(hamiltonian_elec_base_t), intent(in) :: this
 
    has_vector_potential = allocated(this%vector_potential) &
      .or. allocated(this%uniform_vector_potential)
 
  end function hamiltonian_elec_base_has_vector_potential
 
  ! ---------------------------------------------------------------------------------------
  subroutine hamiltonian_elec_base_rashba(this, mesh, der, std, psib, vpsib)
    class(hamiltonian_elec_base_t),  intent(in)    :: this
    class(mesh_t),                  intent(in)    :: mesh
    type(derivatives_t),            intent(in)    :: der
    type(states_elec_dim_t),        intent(in)    :: std
    type(wfs_elec_t), target,       intent(in)    :: psib
    type(wfs_elec_t), target,       intent(inout) :: vpsib
 
    integer :: ist, idim, ip
    complex(real64), allocatable :: psi(:, :), vpsi(:, :), grad(:, :, :)
 
    push_sub(hamiltonian_elec_base_rashba)
 
    if (abs(this%rashba_coupling) < m_epsilon) then
      pop_sub(hamiltonian_elec_base_rashba)
      return
    end if
    assert(std%ispin == spinors)
    assert(mesh%box%dim == 2)
    assert(psib%type() == type_cmplx)
    assert(vpsib%type() == type_cmplx)
 
    safe_allocate(psi(1:mesh%np_part, 1:std%dim))
    safe_allocate(vpsi(1:mesh%np, 1:std%dim))
    safe_allocate(grad(1:mesh%np, 1:mesh%box%dim, 1:std%dim))
 
    do ist = 1, psib%nst
      call batch_get_state(psib, ist, mesh%np_part, psi)
      call batch_get_state(vpsib, ist, mesh%np, vpsi)
 
      do idim = 1, std%dim
        call zderivatives_grad(der, psi(:, idim), grad(:, :, idim), ghost_update = .false., set_bc = .false.)
      end do
 
      if (allocated(this%vector_potential)) then
        do ip = 1, mesh%np
          vpsi(ip, 1) = vpsi(ip, 1) - &
            (this%rashba_coupling) * cmplx(this%vector_potential(2, ip), this%vector_potential(1, ip), real64) * psi(ip, 2)
          vpsi(ip, 2) = vpsi(ip, 2) - &
            (this%rashba_coupling) * cmplx(this%vector_potential(2, ip), -this%vector_potential(1, ip), real64) * psi(ip, 1)
        end do
      end if
 
      do ip = 1, mesh%np
        vpsi(ip, 1) = vpsi(ip, 1) - &
          this%rashba_coupling*( grad(ip, 1, 2) - m_zi*grad(ip, 2, 2))
        vpsi(ip, 2) = vpsi(ip, 2) + &
          this%rashba_coupling*( grad(ip, 1, 1) + m_zi*grad(ip, 2, 1))
      end do
 
      call batch_set_state(vpsib, ist, mesh%np, vpsi)
    end do
 
    safe_deallocate_a(grad)
    safe_deallocate_a(vpsi)
    safe_deallocate_a(psi)
 
    pop_sub(hamiltonian_elec_base_rashba)
  end subroutine hamiltonian_elec_base_rashba
 
#include "undef.F90"
#include "real.F90"
#include "hamiltonian_elec_base_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "hamiltonian_elec_base_inc.F90"
 
end module hamiltonian_elec_base_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: