ylm_wannier.F90

Go to the documentation of this file.

 !
 ! Copyright (C) 2003-2013 Quantum ESPRESSO and Wannier90 groups
 ! This file is distributed under the terms of the
 ! GNU General Public License. See the file "License"
 ! in the root directory of the present distribution,
 ! or http://www.gnu.org/copyleft/gpl.txt .
 !
 ! routine adapted from the original Pwscf interface PP/pw2wannier90.f90
 !
#include "global.h"
 
module ylm_wannier_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use loct_math_oct_m
 
  implicit none
 
  private
 
  public :: ylm_wannier
 
contains
  subroutine ylm_wannier(ylm, l, mr, rr, xx, nr)
    ! this routine returns in ylm(r) the values at the nr points r(1:3,1:nr)
    ! of the spherical harmonic identified  by indices (l,mr)
    ! in table 3.1 of the wannierf90 specification.
    !
    ! No reference to the particular ylm ordering internal to octopus
    ! is assumed.
    !
    ! If ordering in wannier90 code is changed or extended this should be the
    ! only place to be modified accordingly
    !
 
    ! I/O variables
    !
    integer, intent(in)    :: l, mr, nr
    real(real64),   intent(out)   :: ylm(nr)
    real(real64),   intent(in)    :: rr(nr), xx(3, nr)
    !
    ! local variables
    !
    !     real(real64), EXTERNAL ::  s, pz,px,py, dz2, dxz, dyz, dx2my2, dxy
    !     real(real64), EXTERNAL :: fz3, fxz2, fyz2, fzx2my2, fxyz, fxx2m3y2, fy3x2my2
    real(real64) :: cost, phi
    integer :: ir
    real(real64) :: bs2, bs3, bs6, bs12
    bs2  = m_one / sqrt(m_two)
    bs3  = m_one / sqrt(m_three)
    bs6  = m_one / sqrt(6.0_real64)
    bs12 = m_one / sqrt(12.0_real64)
 
    if (l > 3 .or. l < -5) stop 'ylm_wannier l out of range '
    if (l >= 0) then
      if (mr < 1 .or. mr > 2*l + 1) stop 'ylm_wannier mr out of range'
    else
      if (mr < 1 .or. mr > abs(l) + 1) stop 'ylm_wannier mr out of range'
    end if
 
    do ir = 1, nr
      ! if r=0, direction is undefined => make ylm=0 except for l=0
      if (rr(ir) < r_small) then
        if (l == 0) then
          ylm(ir) = s()
        else
          ylm(ir) = m_zero
        end if
        cycle
      end if
 
      cost =  xx(3, ir) / rr(ir)
      !
      !  beware the atan, it is defined modulo M_PI
      !
      phi = atan2(xx(2, ir), xx(1, ir))
 
      select case (l)
      case (0)   ! s orbital
        ylm(ir) = s()
      case (1)   ! p orbitals
        if (mr == 1) ylm(ir) = pz(cost)
        if (mr == 2) ylm(ir) = px(cost, phi)
        if (mr == 3) ylm(ir) = py(cost, phi)
      case (2)   ! d orbitals
        if (mr == 1) ylm(ir) = dz2(cost)
        if (mr == 2) ylm(ir) = dxz(cost, phi)
        if (mr == 3) ylm(ir) = dyz(cost, phi)
        if (mr == 4) ylm(ir) = dx2my2(cost, phi)
        if (mr == 5) ylm(ir) = dxy(cost, phi)
      case (3)   ! f orbitals
        if (mr == 1) ylm(ir) = fz3(cost)
        if (mr == 2) ylm(ir) = fxz2(cost, phi)
        if (mr == 3) ylm(ir) = fyz2(cost, phi)
        if (mr == 4) ylm(ir) = fzx2my2(cost, phi)
        if (mr == 5) ylm(ir) = fxyz(cost, phi)
        if (mr == 6) ylm(ir) = fxx2m3y2(cost, phi)
        if (mr == 7) ylm(ir) = fy3x2my2(cost, phi)
      case (-1)  !  sp hybrids
        if (mr == 1) ylm(ir) = bs2 * (s() + px(cost, phi))
        if (mr == 2) ylm(ir) = bs2 * (s() - px(cost, phi))
      case (-2)  !  sp2 hybrids
        if (mr == 1) ylm(ir) = bs3 * s() - bs6 * px(cost, phi) + bs2 * py(cost, phi)
        if (mr == 2) ylm(ir) = bs3 * s() - bs6 * px(cost, phi) - bs2 * py(cost, phi)
        if (mr == 3) ylm(ir) = bs3 * s() + m_two * bs6 * px(cost, phi)
      case (-3)  !  sp3 hybrids
        if (mr == 1) ylm(ir) = m_half*(s() + px(cost, phi) + py(cost, phi) + pz(cost))
        if (mr == 2) ylm(ir) = m_half*(s() + px(cost, phi) - py(cost, phi) - pz(cost))
        if (mr == 3) ylm(ir) = m_half*(s() - px(cost, phi) + py(cost, phi) - pz(cost))
        if (mr == 4) ylm(ir) = m_half*(s() - px(cost, phi) - py(cost, phi) + pz(cost))
      case (-4)  !  sp3d hybrids
        if (mr == 1) ylm(ir) = bs3 * s() - bs6 * px(cost, phi) + bs2 * py(cost, phi)
        if (mr == 2) ylm(ir) = bs3 * s() - bs6 * px(cost, phi) - bs2 * py(cost, phi)
        if (mr == 3) ylm(ir) = bs3 * s() + m_two * bs6 * px(cost, phi)
        if (mr == 4) ylm(ir) = bs2 * pz(cost) + bs2 * dz2(cost)
        if (mr == 5) ylm(ir) =-bs2 * pz(cost) + bs2 * dz2(cost)
      case (-5)  ! sp3d2 hybrids
        if (mr == 1) ylm(ir) = bs6 * s() - bs2 * px(cost, phi) -bs12 * dz2(cost) + m_half * dx2my2(cost, phi)
        if (mr == 2) ylm(ir) = bs6 * s() + bs2 * px(cost, phi) -bs12 * dz2(cost) + m_half * dx2my2(cost, phi)
        if (mr == 3) ylm(ir) = bs6 * s() - bs2 * py(cost, phi) -bs12 * dz2(cost) - m_half * dx2my2(cost, phi)
        if (mr == 4) ylm(ir) = bs6 * s() + bs2 * py(cost, phi) -bs12 * dz2(cost) - m_half * dx2my2(cost, phi)
        if (mr == 5) ylm(ir) = bs6 * s() - bs2 * pz(cost) +bs3 * dz2(cost)
        if (mr == 6) ylm(ir) = bs6 * s() + bs2 * pz(cost) +bs3 * dz2(cost)
      end select
 
    end do
  end subroutine ylm_wannier
 
  !======== l = 0 =====================================================================
  function s()
    real(real64) :: s
    s = m_one / sqrt((m_four * m_pi))
  end function s
 
  !======== l = 1 =====================================================================
  function pz(cost)
    real(real64) ::pz, cost
    pz =  sqrt(m_three / (m_four * m_pi)) * cost
  end function pz
 
  function px(cost, phi)
    real(real64) ::px, cost, phi, sint
    sint = sqrt(abs(m_one - cost * cost))
    px =  sqrt(m_three / (m_four * m_pi)) * sint * cos(phi)
  end function px
 
  function py(cost, phi)
    real(real64) ::py, cost, phi, sint
    sint = sqrt(abs(m_one - cost * cost))
    py =  sqrt(m_three / (m_four * m_pi)) * sint * sin(phi)
  end function py
 
  !======== l = 2 =====================================================================
  function dz2(cost)
    real(real64) ::dz2, cost
    dz2 =  sqrt(1.25_real64 / (m_four * m_pi)) * (m_three* cost * cost - m_one)
  end function dz2
 
  function dxz(cost, phi)
    real(real64) ::dxz, cost, phi, sint
    sint = sqrt(abs(m_one - cost * cost))
    dxz =  sqrt(15.0_real64 / (m_four * m_pi)) * sint * cost * cos(phi)
  end function dxz
 
  function dyz(cost, phi)
    real(real64) ::dyz, cost, phi, sint
    sint = sqrt(abs(m_one - cost * cost))
    dyz =  sqrt(15.0_real64 / (m_four * m_pi)) * sint * cost * sin(phi)
  end function dyz
 
  function dx2my2(cost, phi)
    real(real64) ::dx2my2, cost, phi, sint
    sint = sqrt(abs(m_one - cost * cost))
    dx2my2 =  sqrt(3.750_real64 / (m_four * m_pi)) * sint * sint * cos(m_two * phi)
  end function dx2my2
 
  function dxy(cost, phi)
    real(real64) ::dxy, cost, phi, sint
    sint = sqrt(abs(m_one - cost * cost))
    dxy =  sqrt(3.750_real64 / (m_four * m_pi)) * sint * sint * sin(m_two * phi)
  end function dxy
 
  !======== l = 3 =====================================================================
  function fz3(cost)
    real(real64) ::fz3, cost
    fz3 =  0.25_real64 * sqrt(7.0_real64 / m_pi) * (5.0_real64 * cost * cost - m_three) * cost
  end function fz3
 
  function fxz2(cost, phi)
    real(real64) ::fxz2, cost, phi, sint
    sint = sqrt(abs(m_one - cost * cost))
    fxz2 =  0.25_real64 * sqrt(10.5_real64 / m_pi) * (5.0_real64 * cost * cost - m_one) * sint * cos(phi)
  end function fxz2
 
  function fyz2(cost, phi)
    real(real64) ::fyz2, cost, phi, sint
    sint = sqrt(abs(m_one - cost * cost))
    fyz2 =  0.25_real64 * sqrt(10.5_real64 / m_pi) * (5.0_real64 * cost * cost - m_one) * sint * sin(phi)
  end function fyz2
 
  function fzx2my2(cost, phi)
    real(real64) ::fzx2my2, cost, phi, sint
    sint = sqrt(abs(m_one - cost * cost))
    fzx2my2 =  0.25_real64 * sqrt(105.0_real64 / m_pi) * sint * sint * cost * cos(m_two * phi)
  end function fzx2my2
 
  function fxyz(cost, phi)
    real(real64) ::fxyz, cost, phi, sint
    sint = sqrt(abs(m_one - cost * cost))
    fxyz =  0.25_real64*sqrt(105.0_real64 / m_pi) * sint * sint * cost * sin(m_two * phi)
  end function fxyz
 
  function fxx2m3y2(cost, phi)
    real(real64) ::fxx2m3y2, cost, phi, sint
    sint = sqrt(abs(m_one - cost * cost))
    fxx2m3y2 =  0.25_real64 * sqrt(17.5_real64 / m_pi) * sint * sint * sint * cos(m_three * phi)
  end function fxx2m3y2
 
  function fy3x2my2(cost, phi)
    real(real64) ::fy3x2my2, cost, phi, sint
    sint = sqrt(abs(m_one - cost * cost))
    fy3x2my2 =  0.25_real64 * sqrt(17.5_real64 / m_pi) * sint * sint * sint * sin(m_three * phi)
  end function fy3x2my2
 
end module ylm_wannier_oct_m