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 projector_oct_m
  use accel_oct_m
  use batch_oct_m
  use batch_ops_oct_m
  use boundaries_oct_m
  use debug_oct_m
  use global_oct_m
  use grid_oct_m
  use hgh_projector_oct_m
  use ions_oct_m
  use kb_projector_oct_m
  use kpoints_oct_m
  use mesh_oct_m
  use messages_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use profiling_oct_m
  use ps_oct_m
  use pseudo_oct_m
  use rkb_projector_oct_m
  use pseudopotential_oct_m
  use species_factory_oct_m
  use splines_oct_m
  use states_elec_dim_oct_m
  use submesh_oct_m
  use wfs_elec_oct_m
 
  implicit none
 
  private
  public ::                    &
    projector_t,               &
    projector_is_null,         &
    projector_is,              &
    projector_init,            &
    projector_init_phases,     &
    projector_build,           &
    projector_check_discretization, &
    projector_end,             &
    dproject_psi,              &
    zproject_psi,              &
    dproject_psi_batch,        &
    zproject_psi_batch,        &
    dprojector_matrix_element, &
    zprojector_matrix_element, &
    dprojector_commute_r,      &
    zprojector_commute_r,      &
    dprojector_commute_r_allatoms_alldir, &
    zprojector_commute_r_allatoms_alldir
 
  integer, parameter :: MAX_NPROJECTIONS = 4
  integer, parameter :: MAX_L = 5
 
 
  type projector_t
    private
    integer, public :: type = proj_none
    integer         :: nprojections
    integer, public :: lmax
    integer, public :: lloc
    integer         :: nik
    integer         :: reltype
 
    type(submesh_t), public  :: sphere
 
 
    type(hgh_projector_t), allocatable, public :: hgh_p(:, :)
    type(kb_projector_t),  allocatable, public :: kb_p(:, :)
    type(rkb_projector_t), allocatable, public :: rkb_p(:, :)
    real(real64),          allocatable, public :: uspp_dij(:, :, :)
    real(real64),          allocatable, public :: uspp_qij(:, :, :)
    real(real64),          allocatable, public :: uspp_qfunc(:, :, :)
    integer,                            public :: uspp_qfunc_lmax = -1
    integer,                            public :: uspp_qfunc_nbeta = 0
    complex(real64),       allocatable, public :: phase(:, :, :)
  end type projector_t
 
contains
 
  !---------------------------------------------------------
  logical elemental function projector_is_null(p)
    type(projector_t), intent(in) :: p
 
    projector_is_null = (p%type == proj_none)
  end function projector_is_null
 
  !---------------------------------------------------------
  logical elemental function projector_is(p, type)
    type(projector_t), intent(in) :: p
    integer,           intent(in) :: type
    projector_is = (p%type == type)
  end function projector_is
 
  !---------------------------------------------------------
  subroutine projector_init(p, pseudo, namespace, dim, reltype)
    type(projector_t),       intent(inout) :: p
    type(pseudopotential_t), target, intent(in) :: pseudo
    type(namespace_t),       intent(in)    :: namespace
    integer,                 intent(in)    :: dim
    integer,                 intent(in)    :: reltype
 
    type(ps_t), pointer :: ps
 
    push_sub(projector_init)
 
    ps => pseudo%ps
 
    p%reltype = reltype
    p%lmax = ps%lmax
 
    if (ps%local) then
      p%type = proj_none
      pop_sub(projector_init)
      return
    end if
 
    p%lloc = ps%llocal
 
    p%type = ps%projector_type
 
    ! reltype is epot_oct_m::proj_reltype, which is either NOREL (0) or SPIN_ORBIT (1);
    ! the named constants cannot be used here, as epot_oct_m already uses this module.
    if (p%type == proj_kb .and. reltype == 1) then
      if (ps%relativistic_treatment == proj_j_dependent) then
        p%type = proj_rkb
      else
        call messages_write("Spin-orbit coupling for species '"//trim(pseudo%get_label())//" is not available.")
        call messages_warning(namespace=namespace)
      end if
    end if
 
    select case (p%type)
    case (proj_kb, proj_rkb, proj_uspp)
      p%nprojections = ps%kbc
    case (proj_hgh)
      p%nprojections = 3
    case default
      assert(.false.)
    end select
 
    pop_sub(projector_init)
  end subroutine projector_init
 
  !---------------------------------------------
 
  subroutine projector_init_phases(this, dim, std, bnd, kpoints, vec_pot, vec_pot_var)
    type(projector_t),             intent(inout) :: this
    integer,                       intent(in)    :: dim
    type(states_elec_dim_t),       intent(in)    :: std
    type(boundaries_t),            intent(in)    :: bnd
    type(kpoints_t),               intent(in)    :: kpoints
    real(real64), optional,  allocatable, intent(in)    :: vec_pot(:)
    real(real64), optional,  allocatable, intent(in)    :: vec_pot_var(:, :)
 
    integer :: ns, iq, is, ikpoint
    real(real64)   :: kr, kpoint(dim)
    integer :: nphase, iphase
    real(real64), allocatable :: diff(:,:)
 
    push_sub(projector_init_phases)
 
    ns = this%sphere%np 
    nphase = 1
    if (bnd%spiralBC) nphase = 3
 
    if (.not. allocated(this%phase) .and. ns > 0) then
      safe_allocate(this%phase(1:ns, 1:nphase, std%kpt%start:std%kpt%end))
    end if
 
    ! Construct vectors which translate the submesh point back into the unit cell:
    !   The positions this%sphere%x can lie outside the unit cell, while
    !   this%sphere%mesh%x(this%sphere%map(is), 1:ndim) by construction is the periodic image inside the unit cell.
    !   If a point of the submesh is inside the unit cell, diff(:,is) = 0.
    safe_allocate(diff(1:dim, 1:ns))
    !$omp parallel private(ikpoint, kpoint, iphase, is, kr, iq)
    !$omp do
    do is = 1, ns
      diff(:, is) = this%sphere%rel_x(:,is) + this%sphere%center - this%sphere%mesh%x(:, this%sphere%map(is))
    end do
 
    do iq = std%kpt%start, std%kpt%end
      ikpoint = std%get_kpoint_index(iq)
 
      ! if this fails, it probably means that sb is not compatible with std
      assert(ikpoint <= kpoints_number(kpoints))
 
      kpoint = m_zero
      kpoint(1:dim) = kpoints%get_point(ikpoint)
 
      do iphase = 1, nphase
        !$omp do
        do is = 1, ns
          ! this is only the correction to the global phase, that can
          ! appear if the sphere crossed the boundary of the cell. (diff=0 otherwise)
 
          kr = sum(kpoint(1:dim)*diff(1:dim, is))
 
          if (present(vec_pot)) then
            if (allocated(vec_pot)) kr = kr + sum(vec_pot(1:dim)*diff(1:dim, is))
          end if
 
          if (present(vec_pot_var)) then
            if (allocated(vec_pot_var)) kr = kr + sum(vec_pot_var(1:dim, this%sphere%map(is)) &
              *(this%sphere%rel_x(:, is)+this%sphere%center))
          end if
 
          if (bnd%spiralBC .and. iphase > 1) then
            kr = kr + (2*(iphase-1)-3)*sum(bnd%spiral_q(1:dim)*diff(1:dim, is))
          end if
 
          this%phase(is, iphase, iq) = exp(-m_zi*kr)
        end do
        !$omp end do nowait
      end do
    end do
    !$omp end parallel
 
    safe_deallocate_a(diff)
 
    pop_sub(projector_init_phases)
 
  end subroutine projector_init_phases
 
  !---------------------------------------------------------
  subroutine projector_build(p, ps, so_strength)
    type(projector_t),        intent(inout) :: p
    class(pseudopotential_t), intent(in)    :: ps
    real(real64),             intent(in)    :: so_strength
 
    integer :: ll, mm, ic, ip
 
    push_sub_with_profile(projector_build)
 
    select case (p%type)
 
    case (proj_hgh)
      safe_allocate(p%hgh_p(0:p%lmax, -p%lmax:p%lmax))
      do ll = 0, p%lmax
        if (ll == p%lloc) cycle
        do mm = -ll, ll
          call hgh_projector_init(p%hgh_p(ll, mm), p%sphere, p%reltype, ps, ll, mm, so_strength)
        end do
      end do
 
    case (proj_kb, proj_uspp)
      safe_allocate(p%kb_p(0:p%lmax, -p%lmax:p%lmax))
      if (p%type == proj_uspp) then
        safe_allocate(p%uspp_dij(0:p%lmax, 1:ps%ps%kbc, 1:ps%ps%kbc))
        safe_allocate(p%uspp_qij(0:p%lmax, 1:ps%ps%kbc, 1:ps%ps%kbc))
        p%uspp_dij = ps%ps%h
        p%uspp_qij = ps%ps%qij
        if (allocated(ps%ps%radial_augmentation)) then
          p%uspp_qfunc_lmax = ps%ps%radial_augmentation_lmax
          p%uspp_qfunc_nbeta = ps%ps%radial_augmentation_nprojectors
          safe_allocate(p%uspp_qfunc(1:p%sphere%np, 1:p%uspp_qfunc_nbeta*p%uspp_qfunc_nbeta, 0:p%uspp_qfunc_lmax))
          do ll = 0, p%uspp_qfunc_lmax
            do ic = 1, p%uspp_qfunc_nbeta*p%uspp_qfunc_nbeta
              do ip = 1, p%sphere%np
                if (p%sphere%r(ip) < spline_range_max(ps%ps%radial_augmentation(ic, ll))) then
                  p%uspp_qfunc(ip, ic, ll) = spline_eval(ps%ps%radial_augmentation(ic, ll), p%sphere%r(ip))
                else
                  p%uspp_qfunc(ip, ic, ll) = m_zero
                end if
              end do
            end do
          end do
        end if
      end if
      do ll = 0, p%lmax
        if (ll == p%lloc) cycle
        do mm = -ll, ll
          call kb_projector_init(p%kb_p(ll, mm), p%sphere, ps, ll, mm)
        end do
      end do
 
    case (proj_rkb)
      safe_allocate(p%rkb_p(1:p%lmax, -p%lmax:p%lmax))
      do ll = 1, p%lmax
        if (ll == p%lloc) cycle
        do mm = -ll, ll
          call rkb_projector_init(p%rkb_p(ll, mm), p%sphere, ps, ll, mm, so_strength)
        end do
      end do
      ! for rkb, l = 0 is a normal kb
      if (p%lloc /= 0) then
        safe_allocate(p%kb_p(1, 1))
        call kb_projector_init(p%kb_p(1, 1), p%sphere, ps, 0, 0)
      end if
 
    end select
 
    pop_sub_with_profile(projector_build)
  end subroutine projector_build
 
  !---------------------------------------------------------
  subroutine projector_check_discretization(pseudo, sm, namespace)
    class(pseudopotential_t), target, intent(in) :: pseudo
    type(submesh_t),          intent(in) :: sm
    type(namespace_t),        intent(in) :: namespace
 
    integer :: ll, mm, ic, ir, nc
    real(real64) :: rc, dr, rr, norm_rad, norm_grid, dev
    real(real64), allocatable :: uv(:)
    type(ps_t), pointer :: ps
 
    integer, parameter :: quad_points = 2000
    real(real64), parameter :: norm_dev_warning = 0.05_real64
 
    push_sub(projector_check_discretization)
 
    ps => pseudo%ps
 
    dev = m_zero
    safe_allocate(uv(1:max(sm%np, 1)))
 
    do ll = 0, ps%lmax
      if (ll == ps%llocal) cycle
 
      nc = ps%kbc
      if (ps%projector_type == proj_kb) nc = min(nc, ps%projectors_per_l(ll+1))
 
      do ic = 1, nc
        ! Norm^2 of the projector on the radial grid
        rc = ps%kb(ll, ic)%x_threshold
        if (rc <= m_epsilon) cycle
        dr = rc/quad_points
        norm_rad = m_zero
        do ir = 1, quad_points
          rr = (ir - m_half)*dr
          norm_rad = norm_rad + spline_eval(ps%kb(ll, ic), rr)**2 * rr**2
        end do
        norm_rad = norm_rad*dr
        if (norm_rad <= 1.0e-10_real64) cycle
 
        ! Norm^2 of the same projector discretized on the grid, for all values of m
        do mm = -ll, ll
          call pseudopotential_real_nl_projector(pseudo, sm%np, sm%rel_x, sm%r, ll, mm, ic, uv)
          norm_grid = dsm_nrm2(sm, uv)**2
          dev = max(dev, abs(norm_grid - norm_rad)/norm_rad)
        end do
      end do
    end do
 
    safe_deallocate_a(uv)
 
    write(message(1), '(a,es9.2)') "Debug: projector discretization error for species '" &
      //trim(pseudo%get_label())//"' = ", dev
    call messages_info(1, namespace=namespace, debug_only=.true.)
 
    if (dev > norm_dev_warning) then
      write(message(1), '(a)') "The nonlocal projectors of species '"//trim(pseudo%get_label())//"'"
      write(message(2), '(a)') "are not well resolved on the real-space grid: at least one projector norm"
      write(message(3), '(a,f8.1,a)') "deviates by ", dev*100, "% from its radial-grid value."
      write(message(4), '(a)') "This can produce spurious (ghost) states. Please reduce the grid Spacing."
      call messages_warning(4, namespace=namespace)
    end if
 
    pop_sub(projector_check_discretization)
  end subroutine projector_check_discretization
 
  !---------------------------------------------------------
  subroutine projector_end(p)
    type(projector_t), intent(inout) :: p
 
    integer :: ll, mm
 
    push_sub(projector_end)
 
    call submesh_end(p%sphere)
 
    select case (p%type)
    case (proj_hgh)
      do ll = 0, p%lmax
        if (ll == p%lloc) cycle
        do mm = -ll, ll
          call hgh_projector_end(p%hgh_p(ll, mm))
        end do
      end do
      safe_deallocate_a(p%hgh_p)
 
    case (proj_kb, proj_uspp)
      do ll = 0, p%lmax
        if (ll == p%lloc) cycle
        do mm = -ll, ll
          call kb_projector_end(p%kb_p(ll, mm))
        end do
      end do
      safe_deallocate_a(p%kb_p)
      safe_deallocate_a(p%uspp_dij)
      safe_deallocate_a(p%uspp_qij)
      safe_deallocate_a(p%uspp_qfunc)
 
    case (proj_rkb)
      do ll = 1, p%lmax
        if (ll == p%lloc) cycle
        do mm = -ll, ll
          call rkb_projector_end(p%rkb_p(ll, mm))
        end do
      end do
      safe_deallocate_a(p%rkb_p)
      if (p%lloc /= 0) then
        call kb_projector_end(p%kb_p(1, 1))
        safe_deallocate_a(p%kb_p)
      end if
 
    end select
 
    p%type = proj_none
 
    safe_deallocate_a(p%phase)
 
    pop_sub(projector_end)
  end subroutine projector_end
 
#include "undef.F90"
#include "real.F90"
#include "projector_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "projector_inc.F90"
 
end module projector_oct_m
 
 
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: