rkb_projector.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 rkb_projector_oct_m
  use debug_oct_m
  use global_oct_m
  use mesh_oct_m
  use messages_oct_m
  use mpi_oct_m
  use submesh_oct_m
  use profiling_oct_m
  use ps_oct_m
  use pseudopotential_oct_m
  use species_oct_m
 
  implicit none
 
  private
  public :: &
    rkb_projector_t,    &
    rkb_projector_init, &
    rkb_project,        &
    rkb_project_bra,    &
    rkb_project_ket,    &
    rkb_projector_end
 
  type rkb_projector_t
    private
    integer            :: n_s
    integer, public    :: n_c
    complex(real64), allocatable :: bra(:, :)
    complex(real64), public, allocatable :: ket(:, :, :, :)
    real(real64), public, allocatable :: f(:, :, :)
  end type rkb_projector_t
 
 
contains
 
  ! ---------------------------------------------------------
  subroutine rkb_projector_init(rkb_p, sm, spec, l, lm, so_strength)
    type(rkb_projector_t),    intent(inout) :: rkb_p
    type(submesh_t),          intent(in)    :: sm
    class(pseudopotential_t), intent(in)    :: spec
    integer,                  intent(in)    :: l, lm
    real(real64),             intent(in)    :: so_strength
 
    integer :: i
 
    push_sub(rkb_projector_init)
 
    rkb_p%n_s = sm%np
    rkb_p%n_c = spec%ps%projectors_per_l(l+1)
    assert(mod(rkb_p%n_c, 2) == 0)
    assert(spec%ps%relativistic_treatment == proj_j_dependent)
 
    !Allocate memory
    safe_allocate(rkb_p%bra(1:rkb_p%n_s, 1:rkb_p%n_c))
    safe_allocate(rkb_p%ket(1:rkb_p%n_s, 1:rkb_p%n_c, 1:2, 1:2))
    safe_allocate(rkb_p%f(1:rkb_p%n_c, 1:2, 1:2))
 
    ! Build projectors
    ! He only need P^{m,m}, P^{m+1,m}, and P^{m-1,m}, for l+1/2 and for l-1/2.
    ! If we use Vanderbilt-type multi-projector per l, we repeat this for each pair of (l+1/2,l-1/2) projectors.
    ! For instead ONCVPSP uses 2 projector per l, so rkb_p%n_c is 4
    do i = 1, rkb_p%n_c
      call pseudopotential_nl_projector(spec, rkb_p%n_s, sm%rel_x, sm%r, l, lm, i, rkb_p%ket(:, i, 1, 1))
      rkb_p%bra(:, i) = conjg(rkb_p%ket(:, i, 1, 1))
      rkb_p%ket(:, i, 2, 2) = rkb_p%ket(:, i, 1, 1)
 
      if (lm /= l) then
        call pseudopotential_nl_projector(spec, rkb_p%n_s, sm%rel_x, sm%r, l, lm+1, i, rkb_p%ket(:, i, 2, 1))
      else
        rkb_p%ket(:, i, 2, 1) = m_z0
      end if
      if (lm /= -l) then
        call pseudopotential_nl_projector(spec, rkb_p%n_s, sm%rel_x, sm%r, l, lm-1, i, rkb_p%ket(:, i, 1, 2))
      else
        rkb_p%ket(:, i, 1, 2) = m_z0
      end if
    end do
 
    ! The l- and m-dependent prefactors are included in the KB energies
    ! See Eq. B.7 in the thesis of M. Oliveira
    do i = 0, rkb_p%n_c/2-1
      rkb_p%f(i*2+1, 1, 1) = real(l + so_strength*lm + 1, real64)                    ! S_{\up\up}
      rkb_p%f(i*2+1, 2, 1) = so_strength*sqrt(real((l + lm + 1)*(l - lm), real64))   ! S_{\down\up}
      rkb_p%f(i*2+1, 1, 2) = so_strength*sqrt(real((l - lm + 1)*(l + lm), real64))   ! S_{\up\down}
      rkb_p%f(i*2+1, 2, 2) = real(l - so_strength*lm + 1, real64)                    ! S_{\down\down}
      rkb_p%f(i*2+2, 1, 1) = real(l - so_strength*lm, real64)                        ! S_{\up\up}
      rkb_p%f(i*2+2, 2, 1) = -so_strength*sqrt(real((l + lm + 1)*(l - lm), real64))  ! S_{\down\up}
      rkb_p%f(i*2+2, 1, 2) = -so_strength*sqrt(real((l - lm + 1)*(l + lm), real64))  ! S_{\up\down}
      rkb_p%f(i*2+2, 2, 2) = real(l + so_strength*lm, real64)                        ! S_{\down\down}
    end do
    rkb_p%f = rkb_p%f/real(2*l + 1, real64)
 
    ! Multiply by the KB energies
    do i = 1, rkb_p%n_c
      rkb_p%f(i, :, :) = rkb_p%f(i, :, :) * spec%ps%h(l, i, i)
    end do
 
    pop_sub(rkb_projector_init)
  end subroutine rkb_projector_init
 
  ! ---------------------------------------------------------
  subroutine rkb_projector_end(rkb_p)
    type(rkb_projector_t), intent(inout) :: rkb_p
 
    push_sub(rkb_projector_end)
 
    safe_deallocate_a(rkb_p%bra)
    safe_deallocate_a(rkb_p%ket)
    safe_deallocate_a(rkb_p%f)
 
    pop_sub(rkb_projector_end)
  end subroutine rkb_projector_end
 
  ! ---------------------------------------------------------
  subroutine rkb_project(mesh, sm, rkb_p, psi, ppsi)
    type(mesh_t),          intent(in)    :: mesh
    type(submesh_t),       intent(in)    :: sm
    type(rkb_projector_t), intent(in)    :: rkb_p
    complex(real64),       intent(in)    :: psi(:, :)
    complex(real64),       intent(inout) :: ppsi(:, :)
 
    complex(real64) :: uvpsi(1:2, 1:rkb_p%n_c)
 
    push_sub(rkb_project)
 
    call rkb_project_bra(mesh, sm, rkb_p, psi, uvpsi)
    call mesh%allreduce(uvpsi)
 
    call rkb_project_ket(rkb_p, uvpsi, ppsi)
 
    pop_sub(rkb_project)
  end subroutine rkb_project
 
  ! ---------------------------------------------------------
  subroutine rkb_project_bra(mesh, sm, rkb_p, psi, uvpsi)
    type(mesh_t),          intent(in)  :: mesh
    type(submesh_t),       intent(in)  :: sm
    type(rkb_projector_t), intent(in)  :: rkb_p
    complex(real64),       intent(in)  :: psi(:, :)
    complex(real64),       intent(out) :: uvpsi(:,:)
 
    integer :: idim, n_s, is, ic
 
    complex(real64), allocatable :: bra(:, :)
 
#ifndef HAVE_OPENMP
    push_sub(rkb_project_bra)
#endif
 
    call profiling_in("RKB_PROJECT_BRA")
 
    uvpsi(1:2, 1:rkb_p%n_c) = m_zero
 
    n_s = rkb_p%n_s
 
    if (mesh%use_curvilinear) then
      safe_allocate(bra(1:n_s, 1:rkb_p%n_c))
      do ic = 1, rkb_p%n_c
        bra(1:n_s, ic) = rkb_p%bra(1:n_s, ic)*mesh%vol_pp(sm%map(1:n_s))
        do idim = 1, 2
          do is = 1, n_s
            uvpsi(idim, ic) = uvpsi(idim, ic) + psi(is, idim)*bra(is, ic)
          end do
        end do
      end do
      safe_deallocate_a(bra)
    else
      do ic = 1, rkb_p%n_c
        do idim = 1, 2
          do is = 1, n_s
            uvpsi(idim, ic) = uvpsi(idim, ic) + psi(is, idim)*rkb_p%bra(is, ic)
          end do
        end do
      end do
    end if
 
    uvpsi(1:2, 1:rkb_p%n_c) = uvpsi(1:2, 1:rkb_p%n_c)*mesh%volume_element
 
    safe_deallocate_a(bra)
 
    call profiling_out("RKB_PROJECT_BRA")
 
#ifndef HAVE_OPENMP
    pop_sub(rkb_project_bra)
#endif
  end subroutine rkb_project_bra
 
  ! ---------------------------------------------------------
  subroutine rkb_project_ket(rkb_p, uvpsi, psi)
    type(rkb_projector_t), intent(in)    :: rkb_p
    complex(real64),       intent(in)    :: uvpsi(:, :)
    complex(real64),       intent(inout) :: psi(:, :)
 
    integer :: idim, jdim, n_s, is, ic
    complex(real64) :: aa
    complex(real64) :: weight(2,2,rkb_p%n_c)
 
#ifndef HAVE_OPENMP
    push_sub(rkb_project_ket)
#endif
 
    call profiling_in("RKB_PROJECT_KET")
 
    n_s = rkb_p%n_s
 
    ! Weight the projectors
    do jdim = 1, 2
      do ic = 1, rkb_p%n_c
        do idim = 1, 2
          weight(jdim, idim, ic) = rkb_p%f(ic, jdim, idim)*uvpsi(idim, ic)
        end do
      end do
    end do
 
    do ic = 1, rkb_p%n_c
      do jdim = 1, 2
        do is = 1, n_s
          aa = m_z0
          do idim = 1, 2
            aa = aa + weight(jdim, idim, ic)*rkb_p%ket(is, ic, jdim, idim)
          end do
          psi(is, jdim) = psi(is, jdim) + aa
        end do
      end do
    end do
 
    call profiling_out("RKB_PROJECT_KET")
 
#ifndef HAVE_OPENMP
    pop_sub(rkb_project_ket)
#endif
  end subroutine rkb_project_ket
 
end module rkb_projector_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: