zora.F90

Go to the documentation of this file.

!! Copyright (C) 2023  L. Konecny,  M. Lueders
!!
!! 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 zora_oct_m
 
  use batch_oct_m
  use batch_ops_oct_m
  use debug_oct_m
  use derivatives_oct_m
  use epot_oct_m
  use global_oct_m
  use grid_oct_m
  use hamiltonian_elec_base_oct_m
  use mesh_oct_m
  use messages_oct_m
  use namespace_oct_m
  use profiling_oct_m
  use states_elec_dim_oct_m
  use wfs_elec_oct_m
 
  implicit none
 
  private
  public :: zora_t
 
  integer, parameter, public :: &
    ZORA_NONE                = 0,     &  !< no ZORA
    zora_scalar_relativistic = 1,     &  
    zora_fully_relativistic  = 2         
 
  !
  type :: zora_t
    private
 
    real(real64), allocatable  :: pot(:,:)
    real(real64), allocatable  :: grad_pot(:,:,:)
    real(real64), allocatable  :: soc(:,:,:)
 
    real(real64)   :: so_strength = m_zero
    real(real64)   :: mass        = m_zero
 
    integer :: zora_level  = zora_none
 
  contains
    procedure :: update => zora_update              
    procedure :: dapply_batch => dzora_apply_batch  
    procedure :: zapply_batch => zzora_apply_batch  
    final     :: zora_finalize                      
  end type
 
  interface zora_t
    procedure zora_constructor
  end interface zora_t
 
contains
 
  !
  function zora_constructor(namespace, der, st_d, ep, mass) result(this)
    type(namespace_t),           intent(in)    :: namespace
    type(derivatives_t),         intent(in)    :: der
    type(states_elec_dim_t),     intent(in)    :: st_d
    type(epot_t),                intent(in)    :: ep
    real(real64),                intent(in)    :: mass
    class(zora_t), pointer                     :: this
 
    push_sub(zora_constructor)
    safe_allocate(this)
 
 
    this%so_strength   = ep%so_strength
    this%mass          = mass
 
    select case (ep%reltype)
    case(scalar_relativistic_zora)
      this%zora_level = zora_scalar_relativistic
    case(fully_relativistic_zora)
      this%zora_level = zora_fully_relativistic
    case default
      this%zora_level = zora_none
      pop_sub(zora_constructor)
      return
    end select
 
    if (der%dim /= 3) then
      call messages_not_implemented("ZORA for dimensions /= 3", namespace)
    end if
 
    safe_allocate(this%pot(1:der%mesh%np_part, 1:st_d%nspin))
    safe_allocate(this%grad_pot(1:der%mesh%np, 1:st_d%nspin, 1:der%dim))
 
    if(this%zora_level == zora_fully_relativistic) then
 
      assert(st_d%nspin == 4)
      safe_allocate(this%soc(1:der%mesh%np, 1:st_d%nspin, 1:der%dim))
 
    end if
 
    pop_sub(zora_constructor)
  end function zora_constructor
 
 
  !
  subroutine zora_finalize(this)
    type(zora_t), intent(inout) :: this
 
    push_sub(zora_finalize)
 
    safe_deallocate_a(this%pot)
    safe_deallocate_a(this%grad_pot)
    safe_deallocate_a(this%soc)
 
    pop_sub(zora_finalize)
  end subroutine zora_finalize
 
 
  !
  subroutine zora_update(this, der, potential)
    class(zora_t),         intent(inout) :: this
    class(derivatives_t), intent(in)    :: der             
    real(real64), contiguous,     intent(in)    :: potential(:, :)
 
    real(real64), allocatable :: zora_scpot(:), grad_pot(:,:)
    real(real64)       :: prefactor, two_c2
 
    integer :: ip, idir
    integer :: nspin
 
    push_sub(zora_update)
 
    if (this%zora_level == zora_none) then
      pop_sub(zora_update)
      return
    end if
 
    nspin = size(this%pot, dim=2)
    two_c2 = m_two * p_c**2
 
    ! We need to copy the potential to an array, allocated to np_part for the application of the derivatives.
 
    safe_allocate( zora_scpot(1:der%mesh%np_part) )
    safe_allocate( grad_pot(1:der%mesh%np_part, 1:der%dim) )
 
    !$omp parallel private(ip)
    !$omp do simd schedule(static)
    do ip = 1, der%mesh%np
      ! zora_scpot = V
      zora_scpot(ip) = potential(ip, 1)
      ! zora_pot   = \frac{c^2}{2 m c^2 - V}
      ! included extra factor of 2 because of -1/2 in kinetic energy
      this%pot(ip,1) = two_c2 / (this%mass * two_c2 - zora_scpot(ip) )
    end do
    !$omp end do simd
 
    if (nspin > 1) then
      !$omp do simd schedule(static)
      do ip = 1, der%mesh%np
        this%pot(ip,2) = this%pot(ip,1)
      end do
      !$omp end do simd
    end if
    if (nspin > 2) then
      !$omp do simd schedule(static)
      do ip = 1, der%mesh%np
        this%pot(ip,3) = m_zero
        this%pot(ip,4) = m_zero
      end do
      !$omp end do simd
    end if
    !$omp end parallel
 
    ! zora_grad_pot   = grad( \frac{c^2}{2 m c^2 - V} )
 
    call dderivatives_grad(der, this%pot(:,1), grad_pot)
 
    do idir=1, der%dim
      !$omp parallel private(ip)
      !$omp do simd schedule(static)
      do ip = 1, der%mesh%np
        this%grad_pot(ip, 1, idir)= grad_pot(ip,idir)
      end do
      !$omp end do simd
 
      if (nspin > 1) then
        !$omp do simd schedule(static)
        do ip = 1, der%mesh%np
          this%grad_pot(ip, 2, idir)= grad_pot(ip,idir)
        end do
        !$omp end do simd
      end if
 
      if (nspin > 2) then
        !$omp do simd schedule(static)
        do ip = 1, der%mesh%np
          this%grad_pot(ip, 3, idir)= m_zero
          this%grad_pot(ip, 4, idir)= m_zero
        end do
        !$omp end do simd
      end if
      !$omp end parallel
 
    end do
 
    ! spin-orbit relativistic ZORA contribution
    if (this%zora_level == zora_fully_relativistic) then
      call zora_update_fully_relativistic()
    end if
 
    safe_deallocate_a( zora_scpot )
    safe_deallocate_a( grad_pot )
 
    pop_sub(zora_update)
 
  contains
 
 
    subroutine zora_update_fully_relativistic()
 
      real(real64), allocatable :: grad_v(:,:)
 
      push_sub(zora_update_fully_relativistic)
 
      safe_allocate(grad_v(1:der%mesh%np, 1:der%dim))
 
 
      ! gradient of potenial, only consider scalar potential
      call dderivatives_grad(der, zora_scpot(:), grad_v(:,:))
 
      !$omp parallel private(ip, prefactor)
      !$omp do simd schedule(static)
      do ip = 1, der%mesh%np
 
        ! added extra factor of 2 because of -1/2 in the gradient
        prefactor = this%so_strength *  two_c2 / (this%mass*two_c2 - zora_scpot(ip))**2
 
        grad_v(ip, 1) = grad_v(ip, 1) * prefactor
        grad_v(ip, 2) = grad_v(ip, 2) * prefactor
        grad_v(ip, 3) = grad_v(ip, 3) * prefactor
 
        ! sigma \times zora_grad_v
        ! four-component vector has elements ( up-up, down-down, Re(up-down), Im(up-down) )
        this%soc(ip, 1, 1) = -grad_v(ip, 2)
        this%soc(ip, 2, 1) =  grad_v(ip, 2)
        this%soc(ip, 3, 1) =  m_zero
        this%soc(ip, 4, 1) = -grad_v(ip, 3)
 
        this%soc(ip, 1, 2) =  grad_v(ip, 1)
        this%soc(ip, 2, 2) = -grad_v(ip, 1)
        this%soc(ip, 3, 2) = -grad_v(ip, 3)
        this%soc(ip, 4, 2) =  m_zero
 
        this%soc(ip, 1, 3) =  m_zero
        this%soc(ip, 2, 3) =  m_zero
        this%soc(ip, 3, 3) =  grad_v(ip, 2)
        this%soc(ip, 4, 3) =  grad_v(ip, 1)
 
      end do
      !$omp end do simd
      !$omp end parallel
 
      safe_deallocate_a(grad_v)
 
      pop_sub(zora_update_fully_relativistic)
    end subroutine zora_update_fully_relativistic
 
  end subroutine zora_update
 
#include "undef.F90"
#include "real.F90"
#include "zora_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "zora_inc.F90"
 
 
end module zora_oct_m