Commit 3d41caec authored by Matthias Redies's avatar Matthias Redies

calculate legendre polys on the fly

parent 53a57db4
......@@ -345,7 +345,7 @@ END SUBROUTINE hsmt_sph_gpu
REAL fleg1(0:atoms%lmaxd),fleg2(0:atoms%lmaxd),fl2p1(0:atoms%lmaxd)
REAL fl2p1bt(0:atoms%lmaxd)
REAL qssbti(3),qssbtj(3)
REAL, ALLOCATABLE :: plegend(:,:)
REAL, ALLOCATABLE :: gdot(:)
REAL, ALLOCATABLE :: VecHelpS(:),VecHelpH(:)
COMPLEX, ALLOCATABLE :: cph(:)
LOGICAL apw(0:atoms%lmaxd)
......@@ -363,16 +363,15 @@ END SUBROUTINE hsmt_sph_gpu
!$OMP PARALLEL DEFAULT(NONE)&
!$OMP SHARED(kj_BlockSize,lapw,atoms,noco,mpi,input,usdus,smat,hmat)&
!$OMP SHARED(jintsp,iintsp,n,fleg1,fleg2,fj,gj,isp,fl2p1,el,e_shift,fl2p1bt,chi)&
!$OMP PRIVATE(kii,ki,ski,kj,plegend,l,qssbti,qssbtj,fct2)&
!$OMP PRIVATE(kii,ki,ski,kj,gdot,l,qssbti,qssbtj,fct2,plegend)&
!$OMP PRIVATE(cph,nn,tnn,fjkiln,gjkiln)&
!$OMP PRIVATE(w1,apw_lo1,apw_lo2,ddnln,elall,fct,apw1)&
!$OMP PRIVATE(capw1,VecHelpS,VecHelpH)&
!$OMP PRIVATE(kj_start,kj_end,kj_BlockNum, kjb)
ALLOCATE(cph(MAXVAL(lapw%nv)))
ALLOCATE(plegend(MAXVAL(lapw%nv),0:atoms%lmaxd))
ALLOCATE(gdot(0:atoms%lmaxd))
ALLOCATE(VecHelpS(MAXVAL(lapw%nv)),VecHelpH(MAXVAL(lapw%nv)))
plegend=0.0
plegend(:,0)=1.0
gdot=0.0
qssbti=MERGE(- noco%qss/2,+ noco%qss/2,jintsp.EQ.1)
qssbtj=MERGE(- noco%qss/2,+ noco%qss/2,iintsp.EQ.1)
!$OMP DO SCHEDULE(DYNAMIC,1)
......@@ -386,11 +385,7 @@ END SUBROUTINE hsmt_sph_gpu
!---> legendre polynomials
DO kj = kj_start,kj_end
plegend(kj,1) = DOT_PRODUCT(lapw%gk(:,kj,iintsp),lapw%gk(:,ki,jintsp))
END DO
DO l = 1,atoms%lmax(n) - 1
plegend(kj_start:kj_end,l+1) = fleg1(l)*plegend(kj_start:kj_end,1)*plegend(kj_start:kj_end,l) -&
fleg2(l)*plegend(kj_start:kj_end,l-1)
gdot(kj,1) = DOT_PRODUCT(lapw%gk(:,kj,iintsp),lapw%gk(:,ki,jintsp))
END DO
!---> set up phase factors
cph(kj_start:kj_end) = 0.0
......@@ -424,31 +419,33 @@ END SUBROUTINE hsmt_sph_gpu
IF (l<=atoms%lnonsph(n)) elall=elall-e_shift!(isp)
IF (smat%l_real) THEN
DO kj = kj_start,kj_end
fct = plegend(kj,l)*fl2p1(l)*&
plegend = LegendrePoly(l,gdot(kj))
fct = plegend*fl2p1(l)*&
( fjkiln*fj(kj,l,iintsp) + gjkiln*gj(kj,l,iintsp)*ddnln )
fct2 = plegend(kj,l)*fl2p1bt(l) * ( fjkiln*gj(kj,l,iintsp) + gjkiln*fj(kj,l,iintsp) )
fct2 = plegend*fl2p1bt(l) * ( fjkiln*gj(kj,l,iintsp) + gjkiln*fj(kj,l,iintsp) )
smat%data_r(kj,kii)=smat%data_r(kj,kii)+REAL(cph(kj))*fct
hmat%data_r(kj,kii)=hmat%data_r(kj,kii) + REAL(cph(kj)) * ( fct * elall + fct2)
!+APW
IF (input%l_useapw) THEN
apw1 = REAL(cph(kj)) * plegend(kj,l) * &
apw1 = REAL(cph(kj)) * plegend * &
( apw_lo1 * fj(kj,l,iintsp) + apw_lo2 * gj(kj,l,iintsp) )
hmat%data_r(kj,kii)=hmat%data_r(kj,kii) + apw1
ENDIF
!-APW
ENDDO
ELSE
plegend = LegendrePoly(l,gdot(kj))
DO kj = kj_start,kj_end
fct = plegend(kj,l)*fl2p1(l)*&
fct = plegend*fl2p1(l)*&
( fjkiln*fj(kj,l,iintsp) + gjkiln*gj(kj,l,iintsp)*ddnln )
fct2 = plegend(kj,l)*fl2p1bt(l) * ( fjkiln*gj(kj,l,iintsp) + gjkiln*fj(kj,l,iintsp) )
fct2 = plegend*fl2p1bt(l) * ( fjkiln*gj(kj,l,iintsp) + gjkiln*fj(kj,l,iintsp) )
VecHelpS(kj) = VecHelpS(kj) + fct
VecHelpH(kj) = VecHelpH(kj) + fct*elall + fct2
IF (input%l_useapw) THEN
capw1 = cph(kj)*plegend(kj,l) *&
capw1 = cph(kj)*plegend *&
( apw_lo1 * fj(kj,l,iintsp) + apw_lo2 * gj(kj,l,iintsp) )
hmat%data_c(kj,kii)=hmat%data_c(kj,kii) + capw1
ENDIF
......@@ -467,7 +464,7 @@ END SUBROUTINE hsmt_sph_gpu
!---> end loop over ki
ENDDO
!$OMP END DO
DEALLOCATE(plegend)
DEALLOCATE(gdot)
DEALLOCATE(cph)
DEALLOCATE(VecHelpS,VecHelpH)
!$OMP END PARALLEL
......
......@@ -32,9 +32,6 @@ contains
p = x**9 - 36./17. * x**7 + 126./85. * x**5 - 84./221. * x**3 + 17./656.
case(10)
p = x**10 - 45./19. * x**8 + 630./323. * x**6 - 210./323. * x**4 + 106./1413. * x**2 - 1./733.
case default
write (*,*) "we don't support l this high (yet)"
stop 7
end select
end function LegendrePoly
......@@ -68,9 +65,6 @@ contains
p = x**9 - 36./17. * x**7 + 126./85. * x**5 - 84./221. * x**3 + 17./656.
case(10)
p = x**10 - 45./19. * x**8 + 630./323. * x**6 - 210./323. * x**4 + 106./1413. * x**2 - 1./733.
case default
write (*,*) "we don't support l this high (yet)"
stop 7
end select
end function LegendrePoly_scalar
......
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