2003-01-20 05:58:50 +08:00
|
|
|
!
|
2009-08-01 21:41:56 +08:00
|
|
|
! Copyright (C) 2001-2007 Quantum ESPRESSO group
|
2003-01-20 05:58:50 +08:00
|
|
|
! 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 .
|
|
|
|
|
!
|
2004-12-13 21:11:22 +08:00
|
|
|
!
|
|
|
|
|
!--------------------------------------------------------------------
|
|
|
|
|
subroutine sph_bes (msh, r, q, l, jl)
|
|
|
|
|
!--------------------------------------------------------------------
|
2021-02-25 00:44:56 +08:00
|
|
|
!! Spherical Bessel function.
|
2003-01-20 05:58:50 +08:00
|
|
|
!
|
2021-04-03 23:51:18 +08:00
|
|
|
USE upf_kinds, only: DP
|
|
|
|
|
USE upf_const, ONLY : eps14
|
2003-01-20 05:58:50 +08:00
|
|
|
!
|
2004-12-13 21:11:22 +08:00
|
|
|
implicit none
|
2003-01-20 05:58:50 +08:00
|
|
|
!
|
2021-02-25 00:44:56 +08:00
|
|
|
integer :: msh
|
|
|
|
|
!! number of grid points points
|
|
|
|
|
integer :: l
|
|
|
|
|
!! angular momentum (-1 <= l <= 6)
|
|
|
|
|
real(DP) :: r (msh)
|
|
|
|
|
!! radial grid
|
|
|
|
|
real(DP) :: q
|
|
|
|
|
!! q
|
|
|
|
|
real(DP) :: jl (msh)
|
|
|
|
|
!! Output: Spherical Bessel function \(j_l(q*r(i))\)
|
2003-01-20 05:58:50 +08:00
|
|
|
!
|
2006-07-26 21:22:14 +08:00
|
|
|
! xseries = convergence radius of the series for small x of j_l(x)
|
2007-05-03 16:38:01 +08:00
|
|
|
real(DP) :: x, xl, xseries = 0.05_dp
|
2006-07-26 21:22:14 +08:00
|
|
|
integer :: ir, ir0
|
|
|
|
|
integer, external:: semifact
|
2003-01-20 05:58:50 +08:00
|
|
|
!
|
2004-11-17 20:05:53 +08:00
|
|
|
#if defined (__MASS)
|
2005-08-28 22:09:42 +08:00
|
|
|
real(DP) :: qr(msh), sin_qr(msh), cos_qr(msh)
|
2003-03-08 01:22:36 +08:00
|
|
|
#endif
|
2006-07-28 21:00:51 +08:00
|
|
|
|
|
|
|
|
! case q=0
|
|
|
|
|
|
2005-05-27 21:49:19 +08:00
|
|
|
if (abs (q) < eps14) then
|
2004-12-13 21:11:22 +08:00
|
|
|
if (l == -1) then
|
|
|
|
|
call errore ('sph_bes', 'j_{-1}(0) ?!?', 1)
|
|
|
|
|
elseif (l == 0) then
|
|
|
|
|
jl(:) = 1.d0
|
|
|
|
|
else
|
|
|
|
|
jl(:) = 0.d0
|
|
|
|
|
endif
|
2006-07-28 01:06:26 +08:00
|
|
|
return
|
|
|
|
|
end if
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2006-07-28 21:00:51 +08:00
|
|
|
! case l=-1
|
2006-07-28 18:38:23 +08:00
|
|
|
|
2006-07-28 01:06:26 +08:00
|
|
|
if (l == - 1) then
|
2006-08-10 00:01:03 +08:00
|
|
|
if (abs (q * r (1) ) < eps14) call errore ('sph_bes', 'j_{-1}(0) ?!?',1)
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2004-11-17 20:05:53 +08:00
|
|
|
#if defined (__MASS)
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2006-08-10 00:01:03 +08:00
|
|
|
qr = q * r
|
|
|
|
|
call vcos( cos_qr, qr, msh)
|
|
|
|
|
jl = cos_qr / qr
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2003-03-08 01:22:36 +08:00
|
|
|
#else
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2006-08-10 00:01:03 +08:00
|
|
|
jl (:) = cos (q * r (:) ) / (q * r (:) )
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2003-03-08 01:22:36 +08:00
|
|
|
#endif
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2006-08-10 00:01:03 +08:00
|
|
|
return
|
|
|
|
|
|
2006-07-28 01:06:26 +08:00
|
|
|
end if
|
|
|
|
|
|
2006-07-28 21:00:51 +08:00
|
|
|
! series expansion for small values of the argument
|
2006-08-10 00:01:03 +08:00
|
|
|
! ir0 is the first grid point for which q*r(ir0) > xseries
|
|
|
|
|
! notice that for small q it may happen that q*r(msh) < xseries !
|
2006-07-28 21:00:51 +08:00
|
|
|
|
2006-08-10 00:01:03 +08:00
|
|
|
ir0 = msh+1
|
2006-07-28 21:00:51 +08:00
|
|
|
do ir = 1, msh
|
|
|
|
|
if ( abs (q * r (ir) ) > xseries ) then
|
|
|
|
|
ir0 = ir
|
|
|
|
|
exit
|
|
|
|
|
end if
|
|
|
|
|
end do
|
|
|
|
|
|
2006-07-28 01:06:26 +08:00
|
|
|
do ir = 1, ir0 - 1
|
|
|
|
|
x = q * r (ir)
|
2006-11-20 18:28:59 +08:00
|
|
|
if ( l == 0 ) then
|
|
|
|
|
xl = 1.0_dp
|
|
|
|
|
else
|
|
|
|
|
xl = x**l
|
|
|
|
|
end if
|
|
|
|
|
jl (ir) = xl/semifact(2*l+1) * &
|
2006-07-26 21:22:14 +08:00
|
|
|
( 1.0_dp - x**2/1.0_dp/2.0_dp/(2.0_dp*l+3) * &
|
|
|
|
|
( 1.0_dp - x**2/2.0_dp/2.0_dp/(2.0_dp*l+5) * &
|
|
|
|
|
( 1.0_dp - x**2/3.0_dp/2.0_dp/(2.0_dp*l+7) * &
|
|
|
|
|
( 1.0_dp - x**2/4.0_dp/2.0_dp/(2.0_dp*l+9) ) ) ) )
|
2006-07-28 01:06:26 +08:00
|
|
|
end do
|
2006-08-10 00:01:03 +08:00
|
|
|
|
|
|
|
|
! the following shouldn't be needed but do you trust compilers
|
|
|
|
|
! to do the right thing in this special case ? I don't - PG
|
|
|
|
|
|
|
|
|
|
if ( ir0 > msh ) return
|
|
|
|
|
|
2006-07-28 01:06:26 +08:00
|
|
|
if (l == 0) then
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2004-11-17 20:05:53 +08:00
|
|
|
#if defined (__MASS)
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2006-07-28 01:06:26 +08:00
|
|
|
qr = q * r
|
|
|
|
|
call vsin( sin_qr, qr, msh)
|
|
|
|
|
jl (ir0:) = sin_qr(ir0:) / (q * r (ir0:) )
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2003-03-08 01:22:36 +08:00
|
|
|
#else
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2006-07-28 01:06:26 +08:00
|
|
|
jl (ir0:) = sin (q * r (ir0:) ) / (q * r (ir0:) )
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2003-03-08 01:22:36 +08:00
|
|
|
#endif
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2006-07-28 01:06:26 +08:00
|
|
|
elseif (l == 1) then
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2004-11-17 20:05:53 +08:00
|
|
|
#if defined (__MASS)
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2006-07-28 01:06:26 +08:00
|
|
|
qr = q * r
|
|
|
|
|
call vcos( cos_qr, qr, msh)
|
|
|
|
|
call vsin( sin_qr, qr, msh)
|
|
|
|
|
jl (ir0:) = ( sin_qr(ir0:) / (q * r (ir0:) ) - &
|
|
|
|
|
cos_qr(ir0:) ) / (q * r (ir0:) )
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2003-03-08 01:22:36 +08:00
|
|
|
#else
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2006-07-28 01:06:26 +08:00
|
|
|
jl (ir0:) = (sin (q * r (ir0:) ) / (q * r (ir0:) ) - &
|
|
|
|
|
cos (q * r (ir0:) ) ) / (q * r (ir0:) )
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2003-03-08 01:22:36 +08:00
|
|
|
#endif
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2006-07-28 01:06:26 +08:00
|
|
|
elseif (l == 2) then
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2004-11-17 20:05:53 +08:00
|
|
|
#if defined (__MASS)
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2006-07-28 01:06:26 +08:00
|
|
|
qr = q * r
|
|
|
|
|
call vcos( cos_qr, qr, msh)
|
|
|
|
|
call vsin( sin_qr, qr, msh)
|
|
|
|
|
jl (ir0:) = ( (3.d0 / (q*r(ir0:)) - (q*r(ir0:)) ) * sin_qr(ir0: ) - &
|
|
|
|
|
3.d0 * cos_qr(ir0:) ) / (q*r(ir0:))**2
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2003-03-08 01:22:36 +08:00
|
|
|
#else
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2006-07-28 01:06:26 +08:00
|
|
|
jl (ir0:) = ( (3.d0 / (q*r(ir0:)) - (q*r(ir0:)) ) * sin (q*r(ir0:)) - &
|
|
|
|
|
3.d0 * cos (q*r(ir0:)) ) / (q*r(ir0:))**2
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2003-03-08 01:22:36 +08:00
|
|
|
#endif
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2006-07-28 01:06:26 +08:00
|
|
|
elseif (l == 3) then
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2004-11-17 20:05:53 +08:00
|
|
|
#if defined (__MASS)
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2006-07-28 01:06:26 +08:00
|
|
|
qr = q * r
|
|
|
|
|
call vcos( cos_qr, qr, msh)
|
|
|
|
|
call vsin( sin_qr, qr, msh)
|
|
|
|
|
jl (ir0:) = (sin_qr (ir0:) * &
|
|
|
|
|
(15.d0 / (q*r(ir0:)) - 6.d0 * (q*r(ir0:)) ) + &
|
|
|
|
|
cos_qr (ir0:) * ( (q*r(ir0:))**2 - 15.d0) ) / &
|
|
|
|
|
(q*r(ir0:))**3
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2003-03-08 01:22:36 +08:00
|
|
|
#else
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2006-07-28 01:06:26 +08:00
|
|
|
jl (ir0:) = (sin (q*r(ir0:)) * &
|
|
|
|
|
(15.d0 / (q*r(ir0:)) - 6.d0 * (q*r(ir0:)) ) + &
|
|
|
|
|
cos (q*r(ir0:)) * ( (q*r(ir0:))**2 - 15.d0) ) / &
|
|
|
|
|
(q*r(ir0:)) **3
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2003-03-08 01:22:36 +08:00
|
|
|
#endif
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2006-07-28 01:06:26 +08:00
|
|
|
elseif (l == 4) then
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2004-11-17 20:05:53 +08:00
|
|
|
#if defined (__MASS)
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2006-07-28 01:06:26 +08:00
|
|
|
qr = q * r
|
|
|
|
|
call vcos( cos_qr, qr, msh)
|
|
|
|
|
call vsin( sin_qr, qr, msh)
|
|
|
|
|
jl (ir0:) = (sin_qr (ir0:) * &
|
|
|
|
|
(105.d0 - 45.d0 * (q*r(ir0:))**2 + (q*r(ir0:))**4) + &
|
|
|
|
|
cos_qr (ir0:) * &
|
|
|
|
|
(10.d0 * (q*r(ir0:))**3 - 105.d0 * (q*r(ir0:))) ) / &
|
|
|
|
|
(q*r(ir0:))**5
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2003-03-08 01:22:36 +08:00
|
|
|
#else
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2006-07-28 01:06:26 +08:00
|
|
|
jl (ir0:) = (sin (q*r(ir0:)) * &
|
|
|
|
|
(105.d0 - 45.d0 * (q*r(ir0:))**2 + (q*r(ir0:))**4) + &
|
|
|
|
|
cos (q*r(ir0:)) * &
|
|
|
|
|
(10.d0 * (q*r(ir0:))**3 - 105.d0 * (q*r(ir0:))) ) / &
|
|
|
|
|
(q*r(ir0:))**5
|
2003-03-08 01:22:36 +08:00
|
|
|
#endif
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2006-07-28 01:06:26 +08:00
|
|
|
elseif (l == 5) then
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2004-11-17 20:05:53 +08:00
|
|
|
#if defined (__MASS)
|
2006-07-28 01:06:26 +08:00
|
|
|
qr = q * r
|
|
|
|
|
call vcos( cos_qr, qr, msh)
|
|
|
|
|
call vsin( sin_qr, qr, msh)
|
|
|
|
|
jl (ir0:) = (-cos_qr(ir0:) - &
|
|
|
|
|
(945.d0*cos_qr(ir0:)) / (q*r(ir0:)) ** 4 + &
|
|
|
|
|
(105.d0*cos_qr(ir0:)) / (q*r(ir0:)) ** 2 + &
|
|
|
|
|
(945.d0*sin_qr(ir0:)) / (q*r(ir0:)) ** 5 - &
|
|
|
|
|
(420.d0*sin_qr(ir0:)) / (q*r(ir0:)) ** 3 + &
|
|
|
|
|
( 15.d0*sin_qr(ir0:)) / (q*r(ir0:)) ) / (q*r(ir0:))
|
2003-03-08 01:22:36 +08:00
|
|
|
#else
|
2006-07-28 01:06:26 +08:00
|
|
|
jl (ir0:) = (-cos(q*r(ir0:)) - &
|
|
|
|
|
(945.d0*cos(q*r(ir0:))) / (q*r(ir0:)) ** 4 + &
|
|
|
|
|
(105.d0*cos(q*r(ir0:))) / (q*r(ir0:)) ** 2 + &
|
|
|
|
|
(945.d0*sin(q*r(ir0:))) / (q*r(ir0:)) ** 5 - &
|
|
|
|
|
(420.d0*sin(q*r(ir0:))) / (q*r(ir0:)) ** 3 + &
|
|
|
|
|
( 15.d0*sin(q*r(ir0:))) / (q*r(ir0:)) ) / (q*r(ir0:))
|
2003-03-08 01:22:36 +08:00
|
|
|
#endif
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2006-07-28 01:06:26 +08:00
|
|
|
elseif (l == 6) then
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2004-11-17 20:05:53 +08:00
|
|
|
#if defined (__MASS)
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2006-07-28 01:06:26 +08:00
|
|
|
qr = q * r
|
|
|
|
|
call vcos( cos_qr, qr, msh)
|
|
|
|
|
call vsin( sin_qr, qr, msh)
|
|
|
|
|
jl (ir0:) = ((-10395.d0*cos_qr(ir0:)) / (q*r(ir0:))**5 + &
|
|
|
|
|
( 1260.d0*cos_qr(ir0:)) / (q*r(ir0:))**3 - &
|
|
|
|
|
( 21.d0*cos_qr(ir0:)) / (q*r(ir0:)) - &
|
|
|
|
|
sin_qr(ir0:) + &
|
|
|
|
|
( 10395.d0*sin_qr(ir0:)) / (q*r(ir0:))**6 - &
|
|
|
|
|
( 4725.d0*sin_qr(ir0:)) / (q*r(ir0:))**4 + &
|
|
|
|
|
( 210.d0*sin_qr(ir0:)) / (q*r(ir0:))**2 ) / (q*r(ir0:))
|
2003-03-08 01:22:36 +08:00
|
|
|
#else
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2006-07-28 01:06:26 +08:00
|
|
|
jl (ir0:) = ((-10395.d0*cos(q*r(ir0:))) / (q*r(ir0:))**5 + &
|
|
|
|
|
( 1260.d0*cos(q*r(ir0:))) / (q*r(ir0:))**3 - &
|
|
|
|
|
( 21.d0*cos(q*r(ir0:))) / (q*r(ir0:)) - &
|
|
|
|
|
sin(q*r(ir0:)) + &
|
|
|
|
|
( 10395.d0*sin(q*r(ir0:))) / (q*r(ir0:))**6 - &
|
|
|
|
|
( 4725.d0*sin(q*r(ir0:))) / (q*r(ir0:))**4 + &
|
|
|
|
|
( 210.d0*sin(q*r(ir0:))) / (q*r(ir0:))**2 ) / (q*r(ir0:))
|
2003-03-08 01:22:36 +08:00
|
|
|
#endif
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2006-07-28 01:06:26 +08:00
|
|
|
else
|
2004-12-13 21:11:22 +08:00
|
|
|
|
2006-07-28 01:06:26 +08:00
|
|
|
call errore ('sph_bes', 'not implemented', abs(l))
|
2004-12-13 21:11:22 +08:00
|
|
|
|
|
|
|
|
endif
|
2003-01-20 05:58:50 +08:00
|
|
|
!
|
2004-12-13 21:11:22 +08:00
|
|
|
return
|
|
|
|
|
end subroutine sph_bes
|
2007-10-30 16:18:46 +08:00
|
|
|
|
|
|
|
|
integer function semifact(n)
|
|
|
|
|
! semifact(n) = n!!
|
|
|
|
|
implicit none
|
|
|
|
|
integer :: n, i
|
|
|
|
|
|
|
|
|
|
semifact = 1
|
|
|
|
|
do i = n, 1, -2
|
|
|
|
|
semifact = i*semifact
|
|
|
|
|
end do
|
|
|
|
|
return
|
|
|
|
|
end function semifact
|
2021-04-05 01:01:26 +08:00
|
|
|
!
|
|
|
|
|
SUBROUTINE sph_dbes1 ( nr, r, xg, l, jl, djl )
|
|
|
|
|
!
|
|
|
|
|
!! Calculates \(x*dj_l(x)/dx\) using the recursion formula:
|
|
|
|
|
!! $$ dj_l(x)/dx = l/x j_l(x) - j_{l+1}(x) $$
|
|
|
|
|
!! for \(l=0\), and for \(l>0\):
|
|
|
|
|
!! $$ dj_l(x)/dx = j_{l-1}(x) - (l+1)/x j_l(x) $$
|
|
|
|
|
!! Requires \(j_l(r)\) in input.
|
|
|
|
|
!! Used only in CP. Note that upflib uses numerical differentiation.
|
|
|
|
|
!
|
|
|
|
|
USE upf_kinds, only: DP
|
|
|
|
|
USE upf_const, ONLY : eps8
|
|
|
|
|
!
|
|
|
|
|
IMPLICIT NONE
|
|
|
|
|
INTEGER, INTENT(IN) :: l, nr
|
|
|
|
|
REAL (DP), INTENT(IN) :: xg, jl(nr), r(nr)
|
|
|
|
|
REAL (DP), INTENT(OUT):: djl(nr)
|
|
|
|
|
!
|
|
|
|
|
if ( xg < eps8 ) then
|
|
|
|
|
!
|
|
|
|
|
! special case q=0
|
|
|
|
|
! note that x*dj_l(x)/dx = 0 for x = 0
|
|
|
|
|
!
|
|
|
|
|
djl(:) = 0.0d0
|
|
|
|
|
else
|
|
|
|
|
!
|
|
|
|
|
if ( l > 0 ) then
|
|
|
|
|
call sph_bes ( nr, r, xg, l-1, djl )
|
|
|
|
|
djl(:) = djl(:) * (xg * r(:) ) - (l+1) * jl(:)
|
|
|
|
|
else if ( l == 0 ) then
|
|
|
|
|
call sph_bes ( nr, r, xg, l+1, djl )
|
|
|
|
|
djl(:) = - djl(:) * (xg * r(:) )
|
|
|
|
|
else
|
|
|
|
|
call errore('sph_dbes','l < 0 not implemented', abs(l) )
|
|
|
|
|
end if
|
|
|
|
|
end if
|
|
|
|
|
!
|
|
|
|
|
end SUBROUTINE sph_dbes1
|
