Commit 0766e4de authored by Miriam Hinzen's avatar Miriam Hinzen

Merge branch 'develop' of iffgit.fz-juelich.de:fleur/fleur into develop

parents 0f570b2b 590e6a79
......@@ -56,7 +56,7 @@ kpoints/tetcon.f kpoints/kvecon.f init/boxdim.f global/radsra.f math/intgr.F glo
)
set(inpgen_F90 ${inpgen_F90} global/constants.f90 io/xsf_io.f90
eigen/vec_for_lo.f90 eigen/orthoglo.F90 juDFT/usage_data.F90 math/ylm4.f90
eigen/vec_for_lo.f90 eigen/orthoglo.F90 juDFT/usage_data.F90 math/ylm4.F90
global/enpara.f90 global/chkmt.f90 inpgen/inpgen.f90 inpgen/set_inp.f90 inpgen/inpgen_help.f90 io/rw_inp.f90 juDFT/juDFT.F90 global/find_enpara.f90
inpgen/closure.f90 inpgen/inpgen_arguments.F90
juDFT/info.F90 juDFT/stop.F90 juDFT/args.F90 juDFT/time.F90 juDFT/init.F90 juDFT/sysinfo.F90 io/w_inpXML.f90 kpoints/julia.f90 global/utility.F90
......
......@@ -507,7 +507,7 @@ MODULE m_corespec_eval
orpref(1:nor) = (4.d0*pi_const)**2
if(.not.allocated(ylm)) allocate(ylm(0:lax*(lax+2),nor))
do ior = 1,nor
CALL ylm4(lax,orvec,ylm(0,ior))
CALL ylm4(lax,orvec,ylm(:,ior))
do la1 = lan,lax ; do mu1 = -la1,la1
lamu = la1*(la1+1)+mu1
write(98,'(3i5,2f12.8)') la1,mu1,lamu,ylm(lamu,ior)
......
......@@ -49,7 +49,7 @@
rg=matmul(kr(:,j),cell%bmat)
CALL ylm4(&
& atoms%lmaxd,rg,&
& ylm(1,j) )!keep
& ylm(:,j) )!keep
ENDDO
ylm = conjg( ylm )
......
......@@ -31,7 +31,7 @@ math/SphBessel.f90
math/DoubleFactorial.f90
math/ExpSave.f90
math/intgr.F90
math/ylm4.f90
math/ylm4.F90
)
if (FLEUR_USE_FFTMKL)
set(fleur_F90 ${fleur_F90} math/mkl_dfti.f90)
......
MODULE m_ylm
IMPLICIT NONE
private
PRIVATE
!---> check whether or not normalizations are needed
REAL, ALLOCATABLE, SAVE :: ynorm(:)
INTEGER, SAVE :: lmaxd = -1 ! initial value
public ylm4,ylmnorm_init
PUBLIC ylm4,ylmnorm_init
#ifdef _CUDA
PUBLIC ylm4_dev
#endif
CONTAINS
#ifdef _CUDA
ATTRIBUTES(device) SUBROUTINE ylm4_dev(lmax,v,ylm)
!************************************************************
! generate the spherical harmonics for the vector v
! using a stable upward recursion in l. (see notes
! by m. weinert.)
! m.weinert january 1982
! modified by R. Podloucky (added in ynorm); July 1989
! cleaned up mw 1995
!
! modified to make use of f90 constructs. note that
! the normalization is an internal subroutine and hence
! can only be called from here. also, no need to dimension
! arrays for ynorm, done dynamically. mw 1999
!************************************************************
! USE m_juDFT
INTEGER, VALUE, INTENT (IN) :: lmax
REAL, DEVICE, INTENT (IN) :: v(3)
COMPLEX, DEVICE, INTENT (OUT):: ylm( (lmax+1)**2 )
REAL, PARAMETER :: small = 1.0e-12
INTEGER :: l,lm0,m
REAL :: fac,x,y,z,xy,r,rxy,cth,sth,cph,sph,cph2
REAL,allocatable,device :: c(:),s(:)
REAL,allocatable,device :: p(:,:)
COMPLEX :: ylms
REAL,allocatable,device :: ynorm_dev(:)
REAL a,cd,fpi
allocate(c(0:lmax),s(0:lmax),ynorm_dev((lmax+1)**2))
allocate(p(0:lmax,0:lmax))
!---> calculate norm
fpi = 4.0 * 3.1415926535897932
DO l=0,lmax
lm0 = l*(l+1) + 1
cd=1.0
a = sqrt( (2*l+1)/fpi )
ynorm_dev(lm0) = a
DO m=1,l
cd=cd/( (l+1-m)*(l+m) )
ynorm_dev(lm0+m) = a*sqrt(cd)
ynorm_dev(lm0-m) = ( (-1.0)**m )*ynorm_dev(lm0+m)
ENDDO
ENDDO
!---> calculate sin and cos of theta and phi
x = v(1)
y = v(2)
z = v(3)
xy = x*x + y*y
r = sqrt(xy + z*z)
rxy = sqrt(xy)
IF (r.GT.small) THEN
cth = z/r
sth = rxy/r
ELSE
sth = 0.0
cth = 1.0
ENDIF
IF(rxy.GT.small) THEN
cph = x/rxy
sph = y/rxy
ELSE
cph = 1.0
sph = 0.0
ENDIF
!---> generate associated legendre functions for m.ge.0
fac = 1.0
!---> loop over m values
DO m = 0, lmax - 1
fac = -(m+m-1)*fac
p(m,m) = fac
p(m+1,m) = (m+m+1)*cth*fac
!---> recurse upward in l
DO l = m+2, lmax
p(l,m)=((l+l-1)*cth*p(l-1,m)-(l+m-1)*p(l-2,m))/(l-m)
ENDDO
fac = fac*sth
ENDDO
p(lmax,lmax) = -(lmax+lmax-1)*fac
!---> determine sin and cos of phi
s(0) = 0.0
s(1) = sph
c(0) = 1.0
c(1) = cph
cph2 = cph+cph
DO m = 2, lmax
s(m) = cph2*s(m-1)-s(m-2)
c(m) = cph2*c(m-1)-c(m-2)
ENDDO
!---> multiply in the normalization factors
DO l=0,lmax
ylm(l*(l+1)+1) = ynorm_dev(l*(l+1)+1)*cmplx(p(l,0),0.0)
ENDDO
DO m = 1, lmax
DO l = m, lmax
lm0 = l*(l+1)+1
ylms = p(l,m)*cmplx(c(m),s(m))
ylm(lm0+m) = ynorm_dev(lm0+m)*ylms
ylm(lm0-m) = conjg(ylms)*ynorm_dev(lm0-m)
ENDDO
ENDDO
END SUBROUTINE ylm4_dev
#endif
SUBROUTINE ylm4(lmax,v,ylm)
!************************************************************
! generate the spherical harmonics for the vector v
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment