Skip to content

fleur
fleur
Source
Select Git revision
...
Target
Select Git revision
 
Commits (20)
Hide whitespace changes
Inline Side-by-side
Showing with 254 additions and 149 deletions
eigen/eigen_hssetup.F90
View file @ 641590dc
......@@ -71,17 +71,19 @@ CONTAINS
ENDDO
CALL timestart("Interstitial part")
! CALL timestart("Interstitial part")
!Generate interstitial part of Hamiltonian
CALL hs_int(input,noco,stars,lapw,mpi,cell,isp,v%pw_w,smat,hmat)
CALL timestop("Interstitial part")
! CALL hs_int(input,noco,stars,lapw,mpi,cell,isp,v%pw_w,smat,hmat)
! CALL timestop("Interstitial part")
CALL timestart("MT part")
!MT-part of Hamiltonian. In case of noco, we need an loop over the local spin of the atoms
DO ispin=MERGE(1,isp,noco%l_noco),MERGE(2,isp,noco%l_noco)
CALL hsmt(atoms,sym,enpara,ispin,input,mpi,noco,cell,lapw,ud,td,smat,hmat)
ENDDO
CALL timestop("MT part")
STOP "hsmt_spherical only"
!Vacuum contributions
IF (input%film) THEN
CALL timestart("Vacuum part")
......
eigen/hsmt.F90
View file @ 641590dc
......@@ -73,10 +73,10 @@ CONTAINS
!that is needed and allocated
CALL hsmt_sph(n,atoms,mpi,ispin,input,noco,1,1,chi_one,lapw,enpara%el0,&
td%e_shift(n,ispin),usdus,fj(:,0:,ispin,:),gj(:,0:,ispin,:),smat(1,1),hmat(1,1))
CALL hsmt_nonsph(n,mpi,sym,atoms,ispin,1,1,chi_one,noco,cell,lapw,td,&
fj(:,0:,ispin,:),gj(:,0:,ispin,:),hmat(1,1))
CALL hsmt_lo(input,atoms,sym,cell,mpi,noco,lapw,usdus,td,fj(:,0:,ispin,:),gj(:,0:,ispin,:),&
n,chi_one,ispin,iintsp,jintsp,hmat(1,1),smat(1,1))
! CALL hsmt_nonsph(n,mpi,sym,atoms,ispin,1,1,chi_one,noco,cell,lapw,td,&
! fj(:,0:,ispin,:),gj(:,0:,ispin,:),hmat(1,1))
! CALL hsmt_lo(input,atoms,sym,cell,mpi,noco,lapw,usdus,td,fj(:,0:,ispin,:),gj(:,0:,ispin,:),&
! n,chi_one,ispin,iintsp,jintsp,hmat(1,1),smat(1,1))
ELSEIF(noco%l_noco.AND..NOT.noco%l_ss) THEN
!The NOCO but non-spinspiral setup follows:
!The Matrix-elements are first calculated in the local frame of the atom and
......
eigen/hsmt_sph.F90
View file @ 641590dc
......@@ -313,42 +313,43 @@ SUBROUTINE hsmt_sph_gpu(n,atoms,mpi,isp,input,noco,iintsp,jintsp,chi,lapw,el,e_s
END SUBROUTINE hsmt_sph_gpu
#endif
SUBROUTINE hsmt_sph_cpu(n,atoms,mpi,isp,input,noco,iintsp,jintsp,chi,lapw,el,e_shift,usdus,fj,gj,smat,hmat)
USE m_constants, ONLY : fpi_const,tpi_const
USE m_types
IMPLICIT NONE
TYPE(t_input),INTENT(IN) :: input
TYPE(t_mpi),INTENT(IN) :: mpi
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_lapw),INTENT(IN) :: lapw
TYPE(t_usdus),INTENT(IN) :: usdus
CLASS(t_mat),INTENT(INOUT) :: smat,hmat
! ..
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: n,isp,iintsp,jintsp
COMPLEX, INTENT(IN) :: chi
! ..
! .. Array Arguments ..
REAL, INTENT (IN) :: el(0:atoms%lmaxd,atoms%ntype,input%jspins)
REAL, INTENT (IN) :: e_shift!(atoms%ntype,input%jspins)
REAL, INTENT (IN) :: fj(:,0:,:),gj(:,0:,:)
! ..
! .. Local Scalars ..
REAL tnn(3), elall,fct,fct2,fjkiln,gjkiln,ddnln,ski(3)
REAL apw_lo1,apw_lo2,apw1,w1
COMPLEX capw1
INTEGER kii,ki,kj,l,nn,kj_end
! ..
! .. Local Arrays ..
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 :: VecHelpS(:),VecHelpH(:)
COMPLEX, ALLOCATABLE :: cph(:)
SUBROUTINE hsmt_sph_cpu(n,atoms,mpi,isp,input,noco,iintsp,jintsp,chi,lapw,el,e_shift,usdus,fj,gj,smat,hmat)
USE m_constants, ONLY : fpi_const,tpi_const
USE m_types
USE m_LegendrePoly
IMPLICIT NONE
TYPE(t_input),INTENT(IN) :: input
TYPE(t_mpi),INTENT(IN) :: mpi
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_lapw),INTENT(IN) :: lapw
TYPE(t_usdus),INTENT(IN) :: usdus
CLASS(t_mat),INTENT(INOUT) :: smat,hmat
! ..
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: n,isp,iintsp,jintsp
COMPLEX, INTENT(IN) :: chi
! ..
! .. Array Arguments ..
REAL, INTENT (IN) :: el(0:atoms%lmaxd,atoms%ntype,input%jspins)
REAL, INTENT (IN) :: e_shift!(atoms%ntype,input%jspins)
REAL, INTENT (IN) :: fj(:,0:,:),gj(:,0:,:)
! ..
! .. Local Scalars ..
REAL tnn(3), elall,fct,fct2,fjkiln,gjkiln,ddnln,ski(3)
REAL apw_lo1,apw_lo2,apw1,w1,plegend
COMPLEX capw1
INTEGER kii,ki,kj,l,nn
INTEGER kj_BlockSize, kj_BlockNum, kjb, kj_start, kj_end
! ..
! .. Local Arrays ..
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 :: gdot(:)
REAL, ALLOCATABLE :: VecHelpS(:),VecHelpH(:)
COMPLEX, ALLOCATABLE :: cph(:),cp_h(:,:)
LOGICAL apw(0:atoms%lmaxd)
CALL timestart("spherical setup")
......@@ -359,111 +360,123 @@ SUBROUTINE hsmt_sph_cpu(n,atoms,mpi,isp,input,noco,iintsp,jintsp,chi,lapw,el,e_s
fl2p1(l) = REAL(l+l+1)/fpi_const
fl2p1bt(l) = fl2p1(l)*0.5
END DO
!$OMP PARALLEL DEFAULT(NONE)&
!$OMP SHARED(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,kj_end,qssbti,qssbtj,fct2)&
!$OMP PRIVATE(cph,nn,tnn,fjkiln,gjkiln)&
!$OMP PRIVATE(w1,apw_lo1,apw_lo2,ddnln,elall,fct,apw1)&
!$OMP PRIVATE(capw1,VecHelpS,VecHelpH)
ALLOCATE(cph(MAXVAL(lapw%nv)))
ALLOCATE(plegend(MAXVAL(lapw%nv),0:atoms%lmaxd))
ALLOCATE(VecHelpS(MAXVAL(lapw%nv)),VecHelpH(MAXVAL(lapw%nv)))
plegend=0.0
plegend(:,0)=1.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)
DO ki = mpi%n_rank+1, lapw%nv(jintsp), mpi%n_size
kii=(ki-1)/mpi%n_size+1
ski = lapw%gvec(:,ki,jintsp) + qssbti
!---> legendre polynomials
DO kj = 1,ki
plegend(kj,1) = DOT_PRODUCT(lapw%gk(:,kj,iintsp),lapw%gk(:,ki,jintsp))
END DO
DO l = 1,atoms%lmax(n) - 1
plegend(:ki,l+1) = fleg1(l)*plegend(:ki,1)*plegend(:ki,l) - fleg2(l)*plegend(:ki,l-1)
END DO
!---> set up phase factors
cph = 0.0
DO nn = SUM(atoms%neq(:n-1))+1,SUM(atoms%neq(:n))
tnn = tpi_const*atoms%taual(:,nn)
DO kj = 1,ki
cph(kj) = cph(kj) +&
CMPLX(COS(DOT_PRODUCT(ski-lapw%gvec(:,kj,iintsp)-qssbtj,tnn)),&
SIN(DOT_PRODUCT(lapw%gvec(:,kj,iintsp)+qssbtj-ski,tnn)))
! IF (iintsp.NE.jintsp) cph(kj)=CONJG(cph(kj))
END DO
END DO
!---> update overlap and l-diagonal hamiltonian matrix
kj_end = MIN(ki,lapw%nv(iintsp))
VecHelpS = 0.d0
VecHelpH = 0.d0
DO l = 0,atoms%lmax(n)
fjkiln = fj(ki,l,jintsp)
gjkiln = gj(ki,l,jintsp)
!
IF (input%l_useapw) THEN
w1 = 0.5 * ( usdus%uds(l,n,isp)*usdus%dus(l,n,isp) + &
usdus%us(l,n,isp)*usdus%duds(l,n,isp) )
apw_lo1 = fl2p1(l) * 0.5 * atoms%rmt(n)**2 * ( gjkiln * w1 +&
fjkiln * usdus%us(l,n,isp) * usdus%dus(l,n,isp) )
apw_lo2 = fl2p1(l) * 0.5 * atoms%rmt(n)**2 * ( fjkiln * w1 +&
gjkiln * usdus%uds(l,n,isp) * usdus%duds(l,n,isp) )
!
ENDIF
ddnln = usdus%ddn(l,n,isp)
elall = el(l,n,isp)
IF (l<=atoms%lnonsph(n)) elall=elall-e_shift!(isp)
IF (smat%l_real) THEN
DO kj = 1,ki
fct = plegend(kj,l)*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) )
ALLOCATE(cp_h(MAXVAL(lapw%nv),SUM(atoms%neq(:n-1))+1:SUM(atoms%neq(:n))))
DO nn = SUM(atoms%neq(:n-1))+1,SUM(atoms%neq(:n))
CALL lapw%phase_factors(1,atoms%taual(:,nn),noco%qss,cp_h(:,nn))
ENDDO
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) * &
kj_BlockSize = 1000
!$OMP PARALLEL DEFAULT(NONE)&
!$OMP SHARED(cp_h,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,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(gdot(MAXVAL(lapw%nv)))
ALLOCATE(VecHelpS(MAXVAL(lapw%nv)),VecHelpH(MAXVAL(lapw%nv)))
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)
DO ki = mpi%n_rank+1, lapw%nv(jintsp), mpi%n_size
kii=(ki-1)/mpi%n_size+1
ski = lapw%gvec(:,ki,jintsp) + qssbti
kj_BlockNum = ki/kj_BlockSize + 1
DO kjb = 1, kj_BlockNum
kj_start = (kjb-1) * kj_BlockSize + 1
kj_end = MIN(ki,kjb * kj_BlockSize,lapw%nv(iintsp))
!---> legendre polynomials
DO kj = kj_start,kj_end
gdot(kj) = DOT_PRODUCT(lapw%gk(:,kj,iintsp),lapw%gk(:,ki,jintsp))
END DO
!---> set up phase factors
cph(kj_start:kj_end) = 0.0
DO nn = SUM(atoms%neq(:n-1))+1,SUM(atoms%neq(:n))
tnn = tpi_const*atoms%taual(:,nn)
DO kj = kj_start,kj_end
cph(kj) = cph(kj) +cp_h(ki,nn)/cp_h(kj,nn)
! CMPLX(COS(DOT_PRODUCT(ski-lapw%gvec(:,kj,iintsp)-qssbtj,tnn)),&
! SIN(DOT_PRODUCT(lapw%gvec(:,kj,iintsp)+qssbtj-ski,tnn)))
! IF (iintsp.NE.jintsp) cph(kj)=CONJG(cph(kj))
END DO
END DO
!---> update overlap and l-diagonal hamiltonian matrix
VecHelpS = 0.d0
VecHelpH = 0.d0
DO l = 0,atoms%lmax(n)
fjkiln = fj(ki,l,jintsp)
gjkiln = gj(ki,l,jintsp)
IF (input%l_useapw) THEN
w1 = 0.5 * ( usdus%uds(l,n,isp)*usdus%dus(l,n,isp) + &
usdus%us(l,n,isp)*usdus%duds(l,n,isp) )
apw_lo1 = fl2p1(l) * 0.5 * atoms%rmt(n)**2 * ( gjkiln * w1 +&
fjkiln * usdus%us(l,n,isp) * usdus%dus(l,n,isp) )
apw_lo2 = fl2p1(l) * 0.5 * atoms%rmt(n)**2 * ( fjkiln * w1 +&
gjkiln * usdus%uds(l,n,isp) * usdus%duds(l,n,isp) )
ENDIF
ddnln = usdus%ddn(l,n,isp)
elall = el(l,n,isp)
IF (l<=atoms%lnonsph(n)) elall=elall-e_shift!(isp)
IF (smat%l_real) THEN
DO kj = kj_start,kj_end
plegend = LegendrePoly_scalar(l,gdot(kj))
fct = plegend*fl2p1(l)*&
( fjkiln*fj(kj,l,iintsp) + gjkiln*gj(kj,l,iintsp)*ddnln )
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 * &
( 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
DO kj = 1,kj_end
fct = plegend(kj,l)*fl2p1(l)*&
hmat%data_r(kj,kii)=hmat%data_r(kj,kii) + apw1
ENDIF
!-APW
ENDDO
ELSE
plegend = LegendrePoly_scalar(l,gdot(kj))
DO kj = kj_start,kj_end
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) )
VecHelpS(kj) = VecHelpS(kj) + fct
VecHelpH(kj) = VecHelpH(kj) + fct*elall + fct2
IF (input%l_useapw) THEN
capw1 = cph(kj)*plegend(kj,l)&
* ( apw_lo1 * fj(kj,l,iintsp) + apw_lo2 * gj(kj,l,iintsp) )
hmat%data_c(kj,kii)=hmat%data_c(kj,kii) + capw1
ENDIF
END DO
ENDIF
!---> end loop over l
ENDDO
IF (.not.smat%l_real) THEN
smat%data_c(:kj_end,kii)=smat%data_c(:kj_end,kii) + chi*cph(:kj_end) * VecHelpS(:kj_end)
hmat%data_c(:kj_end,kii)=hmat%data_c(:kj_end,kii) + chi*cph(:kj_end) * VecHelpH(:kj_end)
ENDIF
!---> end loop over ki
ENDDO
!$OMP END DO
DEALLOCATE(plegend)
DEALLOCATE(cph)
DEALLOCATE(VecHelpS,VecHelpH)
!$OMP END PARALLEL
CALL timestop("spherical setup")
RETURN
END SUBROUTINE hsmt_sph_cpu
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 *&
( apw_lo1 * fj(kj,l,iintsp) + apw_lo2 * gj(kj,l,iintsp) )
hmat%data_c(kj,kii)=hmat%data_c(kj,kii) + capw1
ENDIF
END DO
ENDIF
!---> end loop over l
ENDDO
IF (.not.smat%l_real) THEN
smat%data_c(kj_start:kj_end,kii)=smat%data_c(kj_start:kj_end,kii) + chi*cph(kj_start:kj_end) * VecHelpS(kj_start:kj_end)
hmat%data_c(kj_start:kj_end,kii)=hmat%data_c(kj_start:kj_end,kii) + chi*cph(kj_start:kj_end) * VecHelpH(kj_start:kj_end)
ENDIF
!---> end loop over kj-blocks
ENDDO
!---> end loop over ki
ENDDO
!$OMP END DO
DEALLOCATE(gdot)
DEALLOCATE(cph)
DEALLOCATE(VecHelpS,VecHelpH)
!$OMP END PARALLEL
CALL timestop("spherical setup")
RETURN
END SUBROUTINE hsmt_sph_cpu
END MODULE m_hsmt_sph
main/fleur.F90
View file @ 641590dc
......@@ -233,11 +233,20 @@ CONTAINS
!!$ END IF
!---< gwf
INQUIRE(file="vtotal",exist=l_cont)
IF (l_cont) THEN
CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
0,results%ef,l_qfix,vtot,"vtotal")
ELSE
CALL timestart("generation of potential")
CALL vgen(hybrid,field,input,xcpot,DIMENSION,atoms,sphhar,stars,vacuum,sym,&
obsolete,cell,oneD,sliceplot,mpi,results,noco,inDen,vTot,vx,vCoul)
CALL timestop("generation of potential")
CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
0,-1.0,results%ef,.FALSE.,vtot,"vtotal")
stop "Potential written"
END IF
#ifdef CPP_MPI
CALL MPI_BARRIER(mpi%mpi_comm,ierr)
#endif
......
math/CMakeLists.txt
View file @ 641590dc
......@@ -32,6 +32,7 @@ math/DoubleFactorial.f90
math/ExpSave.f90
math/intgr.F90
math/ylm4.F90
math/LegendrePoly.f90
)
if (FLEUR_USE_FFTMKL)
set(fleur_F90 ${fleur_F90} math/mkl_dfti.f90)
......
math/LegendrePoly.f90 0 → 100644
View file @ 641590dc
module m_LegendrePoly
implicit none
contains
pure recursive function LegendrePoly(l,x) result(p)
implicit none
integer, intent(in) :: l
real, intent(in) :: x(:)
real :: p(size(x))
real, parameter :: one_128 = 1.0/128.0, one_256 = 1.0/256.0
select case(l)
case(0)
p = 1.0
case(1)
p = x
case(2)
p = 0.5 * ( 3 * x**2 - 1)
case(3)
p = 0.5 * ( 5 * x**3 - 3 * x)
case(4)
p = 0.125 * (35 * x**4 - 30 * x**2 + 3)
case(5)
p = 0.125 * (63 * x**5 - 70 * x**3 + 15 * x)
case(6)
p = 0.0625 * (231 * x**6 -315 * x**4 + 105 * x**2 - 5)
case(7)
p = 0.0625 * (429 * x**7 -693 * x**5 + 315 * x**3 - 35 * x)
case(8)
p = one_128 * (6435 * x**8 -12012 *x**6 + 6930* x**4 - 1260 * x**2 + 35)
case(9)
p = one_128 * (12155* x**9 -25740 *x**7 +18018* x**5 - 4620 * x**3 + 315 * x)
case(10)
p = one_256 * (46189* x**10 -109395*x**8 +90090* x**6 -30030 * x**4 +3465 * x**2 - 63)
case default
p = ( (2*l-1)*x*LegendrePoly(l-1,x) - (l-1)*LegendrePoly(l-2,x) ) / l
end select
end function LegendrePoly
pure recursive function LegendrePoly_scalar(l,x) result(p)
implicit none
integer, intent(in) :: l
real, intent(in) :: x
real :: p
real, parameter :: one_128 = 1.0/128.0, one_256 = 1.0/256.0
select case(l)
case(0)
p = 1.0
case(1)
p = x
case(2)
p = 0.5 * ( 3 * x**2 - 1)
case(3)
p = 0.5 * ( 5 * x**3 - 3 * x)
case(4)
p = 0.125 * (35 * x**4 - 30 * x**2 + 3)
case(5)
p = 0.125 * (63 * x**5 - 70 * x**3 + 15 * x)
case(6)
p = 0.0625 * (231 * x**6 -315 * x**4 + 105 * x**2 - 5)
case(7)
p = 0.0625 * (429 * x**7 -693 * x**5 + 315 * x**3 - 35 * x)
case(8)
p = one_128 * (6435 * x**8 -12012 *x**6 + 6930* x**4 - 1260 * x**2 + 35)
case(9)
p = one_128 * (12155* x**9 -25740 *x**7 +18018* x**5 - 4620 * x**3 + 315 * x)
case(10)
p = one_256 * (46189* x**10 -109395*x**8 +90090* x**6 -30030 * x**4 +3465 * x**2 - 63)
case default
p = ( (2*l-1)*x*LegendrePoly_scalar(l-1,x) - (l-1)*LegendrePoly_scalar(l-2,x) ) / l
end select
end function LegendrePoly_scalar
end module m_LegendrePoly
types/types_lapw.F90
View file @ 641590dc
......@@ -371,10 +371,12 @@ CONTAINS
INTEGER:: k
REAL:: th
!$OMP PARALLEL DO PRIVATE(k,th) SHARED(tau,cph,lapw,iintsp,qss) DEFAULT(none)
DO k = 1,lapw%nv(iintsp)
th= DOT_PRODUCT(lapw%gvec(:,k,iintsp)+(iintsp-1.5)*qss,tau)
cph(k) = CMPLX(COS(tpi_const*th),-SIN(tpi_const*th))
END DO
!$OMP END PARALLEL DO
END SUBROUTINE lapw_phase_factors
......@@ -426,6 +428,7 @@ CONTAINS
ENDIF
!---> set up phase factors
!$OMP PARALLEL DO PRIVATE(k,th) SHARED(rph,cph)
DO k = 1,lapw%nv(iintsp)
th= tpi_const*DOT_PRODUCT((/lapw%k1(k,iintsp),lapw%k2(k,iintsp),lapw%k3(k,iintsp)/)+qssbti,atoms%taual(:,na))
rph(k,iintsp) = COS(th)
......