Commit b44aa77e authored by Alexander Neukirchen's avatar Alexander Neukirchen
Browse files

Merge branch 'juPhon/DFPT' into 'develop'

Juphon/dfpt to dev

See merge request fleur/fleur!203
parents 8895136e 68ad957b
......@@ -252,12 +252,12 @@ CONTAINS
+ denCoeffsOffdiag%du21(l,itype)*( g(j,1,l,2)*f(j,1,l,1) +g(j,2,l,2)*f(j,2,l,1) )&
+ denCoeffsOffdiag%dd21(l,itype)*( g(j,1,l,2)*g(j,1,l,1) +g(j,2,l,2)*g(j,2,l,1) )
!rho21(j,0,itype) = rho21(j,0,itype)+ conjg(cs)/(atoms%neq(itype)*sfp_const)
rho21=CONJG(cs)/(atoms%neq(itype)*sfp_const)
rho(j,0,itype,3)=rho(j,0,itype,3)+REAL(rho21)
rho(j,0,itype,4)=rho(j,0,itype,4)-aimag(rho21)
rho21 = cs/(atoms%neq(itype)*sfp_const)
rho(j,0,itype,3) = rho(j,0,itype,3) + REAL(rho21)
rho(j,0,itype,4) = rho(j,0,itype,4) + AIMAG(rho21)
IF (l.LE.input%lResMax) THEN
moments%rhoLRes(j,0,llp,itype,3) = moments%rhoLRes(j,0,llp,itype,3)+REAL(conjg(cs)/(atoms%neq(itype)*sfp_const))
moments%rhoLRes(j,0,llp,itype,4) = moments%rhoLRes(j,0,llp,itype,4)-AIMAG(conjg(cs)/(atoms%neq(itype)*sfp_const))
moments%rhoLRes(j,0,llp,itype,3) = moments%rhoLRes(j,0,llp,itype,3)+ REAL(cs/(atoms%neq(itype)*sfp_const))
moments%rhoLRes(j,0,llp,itype,4) = moments%rhoLRes(j,0,llp,itype,4)+ AIMAG(cs/(atoms%neq(itype)*sfp_const))
END IF
ENDDO
ENDDO
......@@ -276,12 +276,12 @@ CONTAINS
+ g(j,2,lp,2)*f(j,2,l,1) )+ denCoeffsOffdiag%ddnmt21(llp,lh,itype)*(g(j,1,lp,2)*g(j,1,l,1)&
+ g(j,2,lp,2)*g(j,2,l,1) )
!rho21(j,lh,itype)= rho21(j,lh,itype)+ CONJG(cs)/atoms%neq(itype)
rho21=CONJG(cs)/atoms%neq(itype)
rho(j,lh,itype,3)=rho(j,lh,itype,3)+REAL(rho21)
rho(j,lh,itype,4)=rho(j,lh,itype,4)-aimag(rho21)
rho21 = cs/atoms%neq(itype)
rho(j,lh,itype,3) = rho(j,lh,itype,3) + REAL(rho21)
rho(j,lh,itype,4) = rho(j,lh,itype,4) + AIMAG(rho21)
IF ((l.LE.input%lResMax).AND.(lp.LE.input%lResMax)) THEN
moments%rhoLRes(j,lh,llpb,itype,3)= moments%rhoLRes(j,lh,llpb,itype,3) + REAL(conjg(cs)/atoms%neq(itype))
moments%rhoLRes(j,lh,llpb,itype,4)= moments%rhoLRes(j,lh,llpb,itype,4) - AIMAG(conjg(cs)/atoms%neq(itype))
moments%rhoLRes(j,lh,llpb,itype,3)= moments%rhoLRes(j,lh,llpb,itype,3) + REAL(cs/atoms%neq(itype))
moments%rhoLRes(j,lh,llpb,itype,4)= moments%rhoLRes(j,lh,llpb,itype,4) + AIMAG(cs/atoms%neq(itype))
END IF
ENDDO
ENDDO
......
......@@ -41,7 +41,6 @@ MODULE m_rhonmt21
! .. Local Scalars ..
COMPLEX coef, cil, coef1
COMPLEX, PARAMETER :: mi = (0.0,-1.0)
COMPLEX :: temp(ne)
#include"cpp_double.h"
......@@ -81,7 +80,7 @@ MODULE m_rhonmt21
IF (lplow.GT.lphi) CYCLE m_loop
DO lp = lplow, lphi,2
cil = mi**(l-lp)
cil = ImagUnit**(lp-l)
coef1 = cil * sphhar%clnu(jmem,lh,ns)
lmp = lp*(lp+1) + mp
......
This diff is collapsed.
......@@ -53,7 +53,7 @@ CONTAINS
!$OMP SHARED(fmpi, stars, bbmat, gvecPr, gvec, kvecPr, kvec) &
!$OMP SHARED(nvPr, nv, iTkin, fact, l_smat, l_fullj, vpw, hmat, smat) &
!$OMP PRIVATE(ikGPr, ikG, ikG0, gPrG, gInd, th, ts, phase, bvecPr, bvec, r2)
DO ikG = fmpi%n_rank + 1, nv, fmpi%n_size
DO ikG = fmpi%n_rank + 1, nv, fmpi%n_size
ikG0 = (ikG-1) / fmpi%n_size + 1
DO ikGPr = 1, MERGE(nvPr, MIN(ikG, nvPr), l_fullj)
gPrG = fact * (gvecPr(:, ikGPr) - gvec(:, ikG))
......@@ -94,8 +94,8 @@ CONTAINS
hmat%data_r(ikGPr, ikG0) = REAL(th)
smat%data_r(ikGPr, ikG0) = REAL(ts)
ELSE
hmat%data_c(ikGPr, ikG0) = th ! This corresponds to the
smat%data_c(ikGPr, ikG0) = ts ! conjugation in eigen.
hmat%data_c(ikGPr, ikG0) = th
smat%data_c(ikGPr, ikG0) = ts
END IF
END DO
END DO
......
......@@ -86,9 +86,9 @@ CONTAINS
!This is for collinear calculations: the (1,1) element of the matrices is all
!that is needed and allocated
CALL hsmt_sph(n,atoms,fmpi,ilSpinPr,input,nococonv,1,1,chi_one,lapw,enpara%el0,td%e_shift(n,ilSpinPr),usdus,fjgj,smat(1,1),hmat(1,1),.FALSE.,.FALSE.,cell%bmat,1)
CALL hsmt_sph(n,atoms,fmpi,ilSpinPr,input,nococonv,1,1,chi_one,lapw,enpara%el0,td%e_shift(n,ilSpinPr),usdus,fjgj,smat(1,1),hmat(1,1),.FALSE.,.FALSE.)
CALL hsmt_nonsph(n,fmpi,sym,atoms,ilSpinPr,ilSpin,1,1,chi_one,noco,nococonv,cell,lapw,td,fjgj,hmat(1,1),.FALSE.)
CALL hsmt_lo(input,atoms,sym,cell,fmpi,noco,nococonv,lapw,usdus,td,fjgj,n,chi_one,ilSpinPr,ilSpin,igSpinPr,igSpin,hmat(1,1),.FALSE.,smat(1,1))
CALL hsmt_lo(input,atoms,sym,cell,fmpi,noco,nococonv,lapw,usdus,td,fjgj,n,chi_one,ilSpinPr,ilSpin,igSpinPr,igSpin,hmat(1,1),.FALSE.,.FALSE.,smat=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
......@@ -98,9 +98,9 @@ CONTAINS
!initialize the non-LO part of hmat_tmp matrix with zeros
CALL hsmt_nonsph(n,fmpi,sym,atoms,ilSpinPr,ilSpinPr,1,1,chi_one,noco,nococonv,cell,lapw,td,fjgj,hmat_tmp,.TRUE.)
!initialize the smat_tmp matrix with zeros
CALL hsmt_sph(n,atoms,fmpi,ilSpinPr,input,nococonv,1,1,chi_one,lapw,enpara%el0,td%e_shift(n,ilSpinPr),usdus,fjgj,smat_tmp,hmat_tmp,.TRUE.,.FALSE.,cell%bmat,1)
CALL hsmt_sph(n,atoms,fmpi,ilSpinPr,input,nococonv,1,1,chi_one,lapw,enpara%el0,td%e_shift(n,ilSpinPr),usdus,fjgj,smat_tmp,hmat_tmp,.TRUE.,.FALSE.)
!initialize the LO part of the matrices with zeros
CALL hsmt_lo(input,atoms,sym,cell,fmpi,noco,nococonv,lapw,usdus,td,fjgj,n,chi_one,ilSpinPr,ilSpin,igSpinPr,igSpin,hmat_tmp,.TRUE.,smat_tmp)
CALL hsmt_lo(input,atoms,sym,cell,fmpi,noco,nococonv,lapw,usdus,td,fjgj,n,chi_one,ilSpinPr,ilSpin,igSpinPr,igSpin,hmat_tmp,.TRUE.,.FALSE.,smat=smat_tmp)
CALL hsmt_spinor(ilSpinPr,n,nococonv,chi)
CALL timestart("hsmt_distspins")
CALL hsmt_distspins(chi,smat_tmp,smat)
......@@ -121,11 +121,11 @@ CONTAINS
DO igSpin=1,2
IF (ilSpinPr==ilSpin) THEN !local diagonal spin
CALL hsmt_sph(n,atoms,fmpi,ilSpinPr,input,nococonv,igSpinPr,igSpin,chi(igSpinPr,igSpin),&
lapw,enpara%el0,td%e_shift(n,ilSpinPr),usdus,fjgj,smat(igSpinPr,igSpin),hmat(igSpinPr,igSpin),.FALSE.,.FALSE.,cell%bmat,1)
lapw,enpara%el0,td%e_shift(n,ilSpinPr),usdus,fjgj,smat(igSpinPr,igSpin),hmat(igSpinPr,igSpin),.FALSE.,.FALSE.)
CALL hsmt_nonsph(n,fmpi,sym,atoms,ilSpinPr,ilSpin,igSpinPr,igSpin,chi(igSpinPr,igSpin),noco,nococonv,cell,&
lapw,td,fjgj,hmat(igSpinPr,igSpin),.FALSE.)
CALL hsmt_lo(input,atoms,sym,cell,fmpi,noco,nococonv,lapw,usdus,td,fjgj,&
n,chi(igSpinPr,igSpin),ilSpinPr,ilSpin,igSpinPr,igSpin,hmat(igSpinPr,igSpin),.FALSE.,smat(igSpinPr,igSpin))
n,chi(igSpinPr,igSpin),ilSpinPr,ilSpin,igSpinPr,igSpin,hmat(igSpinPr,igSpin),.FALSE.,.FALSE.,smat=smat(igSpinPr,igSpin))
ELSE
IF (any(noco%l_unrestrictMT).OR.noco%l_spinoffd_ldau(n)) call hsmt_mtNocoPot_offdiag(n,input,fmpi,sym,atoms,noco,nococonv,cell,lapw,usdus,td,fjgj,igSpinPr,igSpin,hmat_tmp,hmat)
IF (any(noco%l_constrained)) CALL hsmt_offdiag(n,atoms,fmpi,nococonv,lapw,td,usdus,fjgj,ilSpinPr,ilSpin,igSpinPr,igSpin,hmat)
......
......@@ -20,8 +20,8 @@ CONTAINS
! This translates to (e.g. order = 1 == acof):
! fjgj = \frac{4\pi}{\sqrt{\Omega}W} * (udot_{l}(R_{\alpha}) * K * j_{l})^{'}(K R_{\alpha})
! - udot_{l}^{'}(R_{\alpha}) * j_{l})(K R_{\alpha})
! abCoeffs(lm, k) = CONJG(c_ph(k,igSpin)) *
! ylm(lm, k) *
! abCoeffs(lm, k) = c_ph(k,igSpin) *
! CONJG(ylm(lm, k)) *
! fjgj%fj(k,l,ilSpin,igSpin)
! A necessary factor i**l is omitted here and instead calculated where
! the coefficients are used.
......
......@@ -16,7 +16,7 @@ MODULE m_hsmt_lo
PRIVATE
PUBLIC hsmt_lo
CONTAINS
SUBROUTINE hsmt_lo(Input,Atoms,Sym,Cell,fmpi,Noco,nococonv,Lapw,Ud,Tlmplm,FjGj,N,Chi,ilSpinPr,ilSpin,igSpinPr,igSpin,Hmat,set0,Smat)
SUBROUTINE hsmt_lo(Input,Atoms,Sym,Cell,fmpi,Noco,nococonv,Lapw,Ud,Tlmplm,FjGj,N,Chi,ilSpinPr,ilSpin,igSpinPr,igSpin,Hmat,set0,l_fullj,Smat,lapwq)
USE m_hlomat
USE m_slomat
USE m_setabc1lo
......@@ -35,7 +35,8 @@ CONTAINS
TYPE(t_usdus),INTENT(IN) :: ud
TYPE(t_tlmplm),INTENT(IN) :: tlmplm
TYPE(t_fjgj),INTENT(IN) :: fjgj
LOGICAL,INTENT(IN) :: set0 !if true, initialize the LO-part of the matrices with zeros
LOGICAL,INTENT(IN) :: l_fullj, set0 !if true, initialize the LO-part of the matrices with zeros
TYPE(t_lapw),OPTIONAL,INTENT(IN) :: lapwq
CLASS(t_mat),INTENT(INOUT)::hmat
CLASS(t_mat),INTENT(INOUT),OPTIONAL::smat
......@@ -119,14 +120,25 @@ CONTAINS
!---> hamiltonian matrix, if they are used for this atom.
IF (ilSpinPr==ilSpin) THEN
IF (.NOT.PRESENT(smat)) CALL judft_error("Bug in hsmt_lo, called without smat")
CALL slomat(input,atoms,sym,fmpi,lapw,cell,nococonv,n,na,&
ilSpinPr,ud, alo1(:,ilSpinPr),blo1(:,ilSpinPr),clo1(:,ilSpinPr),fjgj,&
igSpinPr,igSpin,chi,smat)
IF (.NOT.PRESENT(smat).AND..NOT.PRESENT(lapwq)) CALL judft_error("Bug in hsmt_lo, called without smat")
IF (PRESENT(lapwq)) THEN
CALL slomat(input,atoms,sym,fmpi,lapw,cell,nococonv,n,na,&
ilSpinPr,ud, alo1(:,ilSpinPr),blo1(:,ilSpinPr),clo1(:,ilSpinPr),fjgj,&
igSpinPr,igSpin,chi,smat,l_fullj,lapwq)
ELSE
CALL slomat(input,atoms,sym,fmpi,lapw,cell,nococonv,n,na,&
ilSpinPr,ud, alo1(:,ilSpinPr),blo1(:,ilSpinPr),clo1(:,ilSpinPr),fjgj,&
igSpinPr,igSpin,chi,smat,l_fullj)
END IF
END IF
CALL timestart("hlomat")
CALL hlomat(input,atoms,fmpi,lapw,ud,tlmplm,sym,cell,noco,nococonv,ilSpinPr,ilSpin,&
n,na,fjgj,alo1,blo1,clo1,igSpinPr,igSpin,chi,hmat)
IF (PRESENT(lapwq)) THEN
CALL hlomat(input,atoms,fmpi,lapw,ud,tlmplm,sym,cell,noco,nococonv,ilSpinPr,ilSpin,&
n,na,fjgj,alo1,blo1,clo1,igSpinPr,igSpin,chi,hmat,l_fullj,lapwq)
ELSE
CALL hlomat(input,atoms,fmpi,lapw,ud,tlmplm,sym,cell,noco,nococonv,ilSpinPr,ilSpin,&
n,na,fjgj,alo1,blo1,clo1,igSpinPr,igSpin,chi,hmat,l_fullj)
END IF
CALL timestop("hlomat")
END IF
END IF
......
......@@ -43,7 +43,7 @@ CONTAINS
chi_one=1.0
!The spin2,1 matrix is calculated(real part of potential)
CALL hsmt_nonsph(n,fmpi,sym,atoms,2,1,igSpinPr,igSpin,chi_one,noco,nococonv,cell,lapw,td,fjgj,hmat_tmp,.TRUE.)
CALL hsmt_lo(input,atoms,sym,cell,fmpi,noco,nococonv,lapw,ud,td,fjgj,n,chi_one,2,1,igSpinPr,igSpin,hmat_tmp,.TRUE.)
CALL hsmt_lo(input,atoms,sym,cell,fmpi,noco,nococonv,lapw,ud,td,fjgj,n,chi_one,2,1,igSpinPr,igSpin,hmat_tmp,.TRUE.,.FALSE.)
!call hmat_tmp%u2l()
CALL hsmt_spinor(3,n,nococonv,chi) !spinor for off-diagonal part
CALL hsmt_distspins(chi,hmat_tmp,hmat)
......@@ -57,7 +57,7 @@ CONTAINS
!The spin1,2 matrix is calculated(imag part of potential)
!chi_one=CMPLX(0.,1.)
CALL hsmt_nonsph(n,fmpi,sym,atoms,1,2,igSpinPr,igSpin,chi_one,noco,nococonv,cell,lapw,td,fjgj,hmat_tmp,.TRUE.)
CALL hsmt_lo(input,atoms,sym,cell,fmpi,noco,nococonv,lapw,ud,td,fjgj,n,chi_one,1,2,igSpinPr,igSpin,hmat_tmp,.TRUE.)
CALL hsmt_lo(input,atoms,sym,cell,fmpi,noco,nococonv,lapw,ud,td,fjgj,n,chi_one,1,2,igSpinPr,igSpin,hmat_tmp,.TRUE.,.FALSE.)
!call hmat_tmp%u2l()
CALL hsmt_spinor(4,n,nococonv,chi)
......
......@@ -46,6 +46,7 @@ CONTAINS
INTEGER :: nn, na, ab_size, l, ll, size_ab_select
INTEGER :: size_data_c, size_ab, size_ab2 !these data-dimensions are not available so easily in openacc, hence we store them
INTEGER :: ikGPr, ikG
REAL :: rchi
COMPLEX :: cchi
LOGICAL :: l_samelapw
......
......@@ -18,7 +18,7 @@ MODULE m_hsmt_sph
CONTAINS
SUBROUTINE hsmt_sph_acc(n,atoms,fmpi,isp,input,nococonv,igSpinPr,igSpin,chi,lapw,el,e_shift,usdus,fjgj,smat,hmat,set0,l_fullj,bmat,idir,lapwq)
SUBROUTINE hsmt_sph_acc(n,atoms,fmpi,isp,input,nococonv,igSpinPr,igSpin,chi,lapw,el,e_shift,usdus,fjgj,smat,hmat,set0,l_fullj,lapwq)
USE m_constants, ONLY : fpi_const, tpi_const
USE m_types
USE m_hsmt_fjgj
......@@ -35,12 +35,11 @@ CONTAINS
TYPE(t_fjgj), INTENT(IN) :: fjgj
CLASS(t_mat), INTENT(INOUT) :: smat, hmat
LOGICAL, INTENT(IN) :: l_fullj, set0 !if true, initialize the smat matrix with zeros
REAL, INTENT(IN) :: bmat(3, 3)
TYPE(t_lapw), OPTIONAL, INTENT(IN) :: lapwq
! Scalar Arguments
INTEGER, INTENT(IN) :: n, isp, igSpinPr, igSpin, idir
INTEGER, INTENT(IN) :: n, isp, igSpinPr, igSpin
COMPLEX, INTENT(IN) :: chi
! Array Arguments
......@@ -63,7 +62,6 @@ CONTAINS
REAL :: dot, fct, fct2
COMPLEX :: cfac
COMPLEX :: pref(3)
LOGICAL :: l_samelapw
......@@ -87,7 +85,7 @@ CONTAINS
qssAdd = MERGE(-nococonv%qss/2, +nococonv%qss/2, igSpin.EQ.1)
qssAddPr = MERGE(-nococonv%qss/2, +nococonv%qss/2, igSpinPr.EQ.1)
!$acc data &
!$acc& copyin(igSpin,igSpinPr,n,fleg1,fleg2,isp,fl2p1,el,e_shift,chi,qssAdd,qssAddPr,l_fullj,bmat,idir)&
!$acc& copyin(igSpin,igSpinPr,n,fleg1,fleg2,isp,fl2p1,el,e_shift,chi,qssAdd,qssAddPr,l_fullj)&
!$acc& copyin(lapw,lapwPr,atoms,fmpi,input,usdus)&
!$acc& copyin(lapw%nv,lapw%gvec,lapw%gk,lapwPr%nv,lapwPr%gvec,lapwPr%gk,lapw%bkpt,lapwPr%bkpt)&
!$acc& copyin(atoms%lmax,atoms%rmt,atoms%lnonsph,atoms%neq,atoms%taual)&
......@@ -101,7 +99,7 @@ CONTAINS
!$acc loop gang
DO ikG = fmpi%n_rank+1, lapw%nv(igSpin), fmpi%n_size
!$acc loop vector independent&
!$acc & PRIVATE(ikGPr,ikG0,ski,plegend,tnn,vechelps,vechelph,xlegend,fjkiln,gjkiln,ddnln,elall,l3,l,fct,fct2,cph_re,cph_im,cfac,pref,dot)
!$acc & PRIVATE(ikGPr,ikG0,ski,plegend,tnn,vechelps,vechelph,xlegend,fjkiln,gjkiln,ddnln,elall,l3,l,fct,fct2,cph_re,cph_im,cfac,dot)
DO ikGPr = 1, MERGE(lapwPr%nv(igSpinPr),MIN(ikG,lapwPr%nv(igSpinPr)),l_fullj)
ikG0 = (ikG-1)/fmpi%n_size + 1
ski = lapw%gvec(:,ikG,igSpin) + qssAdd(:) + lapw%bkpt
......@@ -171,11 +169,11 @@ CONTAINS
END DO ! nn
cfac = CMPLX(cph_re,cph_im)
! Prefactor: i * (k + G + qssAdd - k' - G' - qssAdd')
IF (l_fullj) THEN
pref = ImagUnit * MATMUL(ski(1:3) - lapwPr%gvec(1:3,ikGPr,igSpinPr) - qssAddPr(1:3) - lapwPr%bkpt, bmat)
cfac = pref(idir) * cfac
END IF
! ! Prefactor: i * (k + G + qssAdd - k' - G' - qssAdd')
! IF (l_fullj) THEN
! pref = ImagUnit * MATMUL(ski(1:3) - lapwPr%gvec(1:3,ikGPr,igSpinPr) - qssAddPr(1:3) - lapwPr%bkpt, bmat)
! cfac = pref(idir) * cfac
! END IF
IF (smat%l_real) THEN
IF (set0) THEN
......@@ -209,7 +207,7 @@ CONTAINS
RETURN
END SUBROUTINE hsmt_sph_acc
SUBROUTINE hsmt_sph_cpu(n,atoms,fmpi,isp,input,nococonv,igSpinPr,igSpin,chi,lapw,el,e_shift,usdus,fjgj,smat,hmat,set0,l_fullj,bmat,idir,lapwq)
SUBROUTINE hsmt_sph_cpu(n,atoms,fmpi,isp,input,nococonv,igSpinPr,igSpin,chi,lapw,el,e_shift,usdus,fjgj,smat,hmat,set0,l_fullj,lapwq)
USE m_constants, ONLY : fpi_const, tpi_const
USE m_types
USE m_hsmt_fjgj
......@@ -226,12 +224,11 @@ CONTAINS
TYPE(t_fjgj), INTENT(IN) :: fjgj
CLASS(t_mat), INTENT(INOUT) :: smat, hmat
LOGICAL, INTENT(IN) :: l_fullj, set0 !if true, initialize the smat matrix with zeros
REAL, INTENT(IN) :: bmat(3, 3)
TYPE(t_lapw), OPTIONAL, INTENT(IN) :: lapwq
! Scalar Arguments
INTEGER, INTENT(IN) :: n, isp, igSpinPr, igSpin, idir
INTEGER, INTENT(IN) :: n, isp, igSpinPr, igSpin
COMPLEX, INTENT(IN) :: chi
! Array Arguments
......@@ -258,7 +255,7 @@ CONTAINS
REAL, ALLOCATABLE :: cph_re(:), cph_im(:)
REAL, ALLOCATABLE :: dot(:), fct(:), fct2(:)
COMPLEX, ALLOCATABLE :: cfac(:), pref(:,:)
COMPLEX, ALLOCATABLE :: cfac(:)
INTEGER, PARAMETER :: NVEC = 128
......@@ -281,13 +278,13 @@ CONTAINS
!$OMP PARALLEL DEFAULT(NONE)&
!$OMP SHARED(lapw,lapwPr,atoms,nococonv,fmpi,input,usdus,smat,hmat)&
!$OMP SHARED(igSpin,igSpinPr,n,fleg1,fleg2,fjgj,isp,fl2p1,el,e_shift,chi,set0,l_fullj,bmat,idir)&
!$OMP SHARED(igSpin,igSpinPr,n,fleg1,fleg2,fjgj,isp,fl2p1,el,e_shift,chi,set0,l_fullj)&
!$OMP PRIVATE(ikG0,ikG,ski,ikGPr,kj_off,kj_vec,plegend,xlegend,l,l3,kj_end,qssAdd,qssAddPr,fct2)&
!$OMP PRIVATE(cph_re,cph_im,cfac,dot,pref,nn,tnn,fjkiln,gjkiln)&
!$OMP PRIVATE(cph_re,cph_im,cfac,dot,nn,tnn,fjkiln,gjkiln)&
!$OMP PRIVATE(w1,apw_lo1,apw_lo2,ddnln,elall,fct)&
!$OMP PRIVATE(VecHelpS,VecHelpH,NVEC_rem)
ALLOCATE(cph_re(NVEC),cph_im(NVEC),cfac(NVEC))
ALLOCATE(dot(NVEC),pref(NVEC,3),fct(NVEC),fct2(NVEC))
ALLOCATE(dot(NVEC),fct(NVEC),fct2(NVEC))
ALLOCATE(plegend(NVEC,0:2))
ALLOCATE(xlegend(NVEC))
ALLOCATE(VecHelpS(NVEC),VecHelpH(NVEC))
......@@ -367,15 +364,15 @@ CONTAINS
DO jv = 0, NVEC_rem-1
ikGPr = jv + kj_off
dot(jv+1) = DOT_PRODUCT(ski(1:3) - lapwPr%gvec(1:3,ikGPr,igSpinPr) - qssAddPr(1:3) - lapwPr%bkpt, tnn(1:3))
pref(jv+1,1:3) = ImagUnit * MATMUL(ski(1:3) - lapwPr%gvec(1:3,ikGPr,igSpinPr) - qssAddPr(1:3) - lapwPr%bkpt, bmat)
!pref(jv+1,1:3) = ImagUnit * MATMUL(ski(1:3) - lapwPr%gvec(1:3,ikGPr,igSpinPr) - qssAddPr(1:3) - lapwPr%bkpt, bmat)
END DO ! ikGPr
cph_re(:NVEC_REM) = cph_re(:NVEC_REM) + COS(dot(:NVEC_REM))
cph_im(:NVEC_REM) = cph_im(:NVEC_REM) + SIN(dot(:NVEC_REM))
cfac(:NVEC_REM) = CMPLX(cph_re(:NVEC_REM),cph_im(:NVEC_REM))
! Prefactor: i * (k + G + qssAdd - k' - G' - qssAdd')
IF (l_fullj) THEN
cfac(:NVEC_REM) = pref(:NVEC_REM, idir) * cfac(:NVEC_REM)
END IF
! ! Prefactor: i * (k + G + qssAdd - k' - G' - qssAdd')
! IF (l_fullj) THEN
! cfac(:NVEC_REM) = pref(:NVEC_REM, idir) * cfac(:NVEC_REM)
! END IF
! IF (igSpinPr.NE.igSpin) cph_im=-cph_im
END DO ! nn
......
......@@ -13,7 +13,7 @@ MODULE m_slomat
CONTAINS
SUBROUTINE slomat(input,atoms,sym,fmpi,lapw,cell,nococonv,ntyp,na,&
isp,ud, alo1,blo1,clo1,fjgj,&
igSpinPr,igSpin,chi,smat)
igSpinPr,igSpin,chi,smat,l_fullj,lapwq)
!***********************************************************************
! locol stores the number of columns already processed; on parallel
! computers this decides, whether the LO-contribution is
......@@ -48,19 +48,40 @@ CONTAINS
REAL, INTENT (IN) :: alo1(atoms%nlod),blo1(atoms%nlod),clo1(atoms%nlod)
TYPE(t_usdus),INTENT(IN) :: ud
CLASS(t_mat),INTENT(INOUT) :: smat
LOGICAL, INTENT(IN) :: l_fullj
TYPE(t_lapw), OPTIONAL, INTENT(IN) :: lapwq
! Local Scalars
REAL :: con,dotp,fact1,fact2,fact3,fl2p1
INTEGER :: invsfct,k ,l,lo,lop,lp,nkvec,nkvecp,kp,i
INTEGER :: locol,lorow
LOGICAL :: l_samelapw
COMPLEX, ALLOCATABLE :: cphPr(:), cph(:)
TYPE(t_lapw) :: lapwPr
l_samelapw = .FALSE.
IF (.NOT.PRESENT(lapwq)) l_samelapw = .TRUE.
IF (.NOT.l_samelapw) THEN
lapwPr = lapwq
ELSE
lapwPr = lapw
END IF
COMPLEX, ALLOCATABLE :: cph(:,:)
ALLOCATE(cph(MAXVAL(lapw%nv)))
ALLOCATE(cphPr(MAXVAL(lapwPr%nv)))
ALLOCATE(cph(MAXVAL(lapw%nv),2))
! TODO: Introduce the logic for different lapw and the full rectangular
! instead of triangular construction...
DO i=MIN(igSpin,igSpinPr),MAX(igSpin,igSpinPr)
CALL lapw%phase_factors(i,atoms%taual(:,na),nococonv%qss,cph(:,i))
END DO
CALL lapw%phase_factors(igSpinPr,atoms%taual(:,na),nococonv%qss,cphPr)
IF (l_samelapw) THEN
cph = cphPr
ELSE
CALL lapw%phase_factors(igSpin,atoms%taual(:,na),nococonv%qss,cph)
END IF
IF ((sym%invsat(na) == 0) .OR. (sym%invsat(na) == 1)) THEN
! If this atom is the first of two atoms related by inversion, the
......@@ -75,9 +96,9 @@ CONTAINS
con = fpi_const/SQRT(cell%omtil)* ((atoms%rmt(ntyp))**2)/2.0
!$acc kernels present(fjgj,fjgj%fj,fjgj%gj,smat,smat%data_c,smat%data_r)&
!$acc & copyin(l,lapw,lapw%kvec(:,:,na),ud,clo1(:),dotp,cph(:,:),atoms,lapw%index_lo(:,na),lapw%gk(:,:,:)) &
!$acc copyin(ud%dulon(:,ntyp,isp),ud%ddn(:,ntyp,isp),ud%uulon(:,ntyp,isp),ud%uloulopn(:,:,ntyp,isp),blo1(:)) &
!$acc copyin(atoms%nlo(ntyp),lapw%nv(:),atoms%llo(:,ntyp),alo1(:), fmpi, fmpi%n_size, fmpi%n_rank)&
!$acc & copyin(l,lapw,lapw%kvec(:,:,na),,lapwPr,lapwPr%kvec(:,:,na)ud,clo1(:),dotp,cph(:),cphPr(:),atoms,lapw%index_lo(:,na),lapw%gk(:,:,:)) &
!$acc copyin(lapwPr%index_lo(:,na),lapwPr%gk(:,:,:),ud%dulon(:,ntyp,isp),ud%ddn(:,ntyp,isp),ud%uulon(:,ntyp,isp),ud%uloulopn(:,:,ntyp,isp),blo1(:)) &
!$acc copyin(atoms%nlo(ntyp),lapw%nv(:),lapwPr%nv(:),atoms%llo(:,ntyp),alo1(:), fmpi, fmpi%n_size, fmpi%n_rank)&
!$acc default(none)
DO lo = 1,atoms%nlo(ntyp) !loop over all LOs for this atom
......@@ -90,6 +111,31 @@ CONTAINS
& + 2*clo1(lo)*ud%dulon(lo,ntyp,isp)) + &
& clo1(lo)*clo1(lo) )
DO nkvec = 1,invsfct* (2*l+1) !Each LO can have several functions
! lop-kG part (l_fullj). |====|
IF (l_fullj) THEN
lorow = lapwPr%nv(igSpinPr)+lapwPr%index_lo(lo,na)+nkvec
kp = lapwPr%kvec(nkvec,lo,na)
!$acc loop gang private(fact2,dotp,kp)
DO k = fmpi%n_rank + 1, lapw%nv(igSpin), fmpi%n_size
fact2 = con * fl2p1 * ( &
& fjgj%fj(k,l,isp,igSpin)*(alo1(lo) &
& + clo1(lo)*ud%uulon(lo,ntyp,isp)) + &
& fjgj%gj(k,l,isp,igSpin)*(blo1(lo)*ud%ddn(l,ntyp,isp) &
& + clo1(lo)*ud%dulon(lo,ntyp,isp)) )
dotp = dot_PRODUCT(lapw%gk(:,k,igSpin),lapwPr%gk(:,kp,igSpinPr))
IF (smat%l_real) THEN
smat%data_r(lorow,k) = smat%data_r(kp,locol) &
& + chi * invsfct * fact2 * legpol(atoms%llo(lo,ntyp),dotp) &
& * CONJG(cphPr(kp)) * cph(k)
ELSE
smat%data_c(lorow,k) = smat%data_c(kp,locol) &
& + chi * invsfct * fact2 * legpol(atoms%llo(lo,ntyp),dotp) &
& * CONJG(cphPr(kp)) * cph(k)
END IF
END DO
!$acc end loop
END IF
locol = lapw%nv(igSpin)+lapw%index_lo(lo,na)+nkvec !this is the column of the matrix
IF (MOD(locol-1,fmpi%n_size) == fmpi%n_rank) THEN
locol=(locol-1)/fmpi%n_size+1 !this is the column in local storage!
......@@ -97,27 +143,28 @@ CONTAINS
! Calculate the overlap matrix elements with the regular Flapw
! basis functions
!$acc loop gang private(fact2,dotp,kp)
DO kp = 1,lapw%nv(igSpinPr)
DO kp = 1,lapwPr%nv(igSpinPr)
fact2 = con * fl2p1 * ( &
& fjgj%fj(kp,l,isp,igSpinPr)*(alo1(lo) &
& + clo1(lo)*ud%uulon(lo,ntyp,isp)) + &
& fjgj%gj(kp,l,isp,igSpinPr)*(blo1(lo)*ud%ddn(l,ntyp,isp) &
& + clo1(lo)*ud%dulon(lo,ntyp,isp)) )
dotp = dot_PRODUCT(lapw%gk(:,k,igSpin),lapw%gk(:,kp,igSpinPr))
dotp = dot_PRODUCT(lapw%gk(:,k,igSpin),lapwPr%gk(:,kp,igSpinPr))
IF (smat%l_real) THEN
smat%data_r(kp,locol) = smat%data_r(kp,locol) &
& + chi * invsfct * fact2 * legpol(atoms%llo(lo,ntyp),dotp) &
& * CONJG(cph(kp,igSpinPr)) * cph(k,igSpin)
& * CONJG(cphPr(kp)) * cph(k)
ELSE
smat%data_c(kp,locol) = smat%data_c(kp,locol) &
& + chi * invsfct * fact2 * legpol(atoms%llo(lo,ntyp),dotp) &
& * CONJG(cph(kp,igSpinPr)) * cph(k,igSpin)
& * CONJG(cphPr(kp)) * cph(k)
END IF
END DO
!$acc end loop
! Calculate the overlap matrix elements with other local orbitals
! of the same atom, if they have the same l
DO lop = 1, MERGE(lo-1,atoms%nlo(ntyp),igSpinPr==igSpin)
DO lop = 1, MERGE(lo-1,atoms%nlo(ntyp),igSpinPr==igSpin.AND..NOT.l_fullj)
IF (lop==lo) CYCLE !Do later
lp = atoms%llo(lop,ntyp)
IF (l == lp) THEN
......@@ -130,17 +177,18 @@ CONTAINS
& + blo1(lo)*ud%dulon(lop,ntyp,isp) &
& + clo1(lo)*ud%uloulopn(lop,lo,ntyp,isp)) )
DO nkvecp = 1,invsfct* (2*lp+1)
kp = lapw%kvec(nkvecp,lop,na)
lorow=lapw%nv(igSpinPr)+lapw%index_lo(lop,na)+nkvecp
dotp = dot_PRODUCT(lapw%gk(:,k,igSpin),lapw%gk(:,kp,igSpinPr))
kp = lapwPr%kvec(nkvecp,lop,na)
lorow=lapwPr%nv(igSpinPr)+lapwPr%index_lo(lop,na)+nkvecp
dotp = dot_PRODUCT(lapw%gk(:,k,igSpin),lapwPr%gk(:,kp,igSpinPr))
IF (smat%l_real) THEN
smat%data_r(lorow,locol) = smat%data_r(lorow,locol) &
& + chi * invsfct * fact3 * legpol(atoms%llo(lo,ntyp),dotp) &
& * CONJG(cph(kp,igSpinPr)) * cph(k,igSpin)
& * CONJG(cphPr(kp)) * cph(k)
ELSE
smat%data_c(lorow,locol) = smat%data_c(lorow,locol) &
& + chi * invsfct * fact3 * legpol(atoms%llo(lo,ntyp),dotp) &
& * CONJG(cph(kp,igSpinPr)) * cph(k,igSpin)
& * CONJG(cphPr(kp)) * cph(k)
END IF
END DO
END IF
......@@ -149,17 +197,18 @@ CONTAINS
! orbital with itself
lop=lo
DO nkvecp = 1,MERGE(nkvec,invsfct* (2*l+1),igSpinPr==igSpin)
kp = lapw%kvec(nkvecp,lo,na)
lorow = lapw%nv(igSpinPr)+lapw%index_lo(lo,na)+nkvecp
dotp = dot_PRODUCT(lapw%gk(:,k,igSpin),lapw%gk(:,kp,igSpinPr))
kp = lapwPr%kvec(nkvecp,lo,na)
lorow = lapwPr%nv(igSpinPr)+lapwPr%index_lo(lo,na)+nkvecp
dotp = dot_PRODUCT(lapw%gk(:,k,igSpin),lapwPr%gk(:,kp,igSpinPr))
IF (smat%l_real) THEN
smat%data_r(lorow,locol) = smat%data_r(lorow,locol) &
& + chi * invsfct * fact1 * legpol(l,dotp) &
& * CONJG(cph(kp, igSpinPr)) * cph(k, igSpin)
& * CONJG(cphPr(kp)) * cph(k)
ELSE
smat%data_c(lorow,locol) = smat%data_c(lorow,locol) &
& + chi * invsfct * fact1 * legpol(l,dotp) &
& * CONJG(cph(kp, igSpinPr)) * cph(k, igSpin)
& * CONJG(cphPr(kp)) * cph(k)
END IF
END DO
END IF ! mod(locol-1,n_size) = nrank
......
MODULE m_tlmplm
USE m_judft
IMPLICIT NONE
!*********************************************************************
! sets up the local Hamiltonian, i.e. the Hamiltonian in the
! l',m',l,m,u- basis which is independent from k!
!*********************************************************************
USE m_judft
IMPLICIT NONE
CONTAINS
SUBROUTINE tlmplm(n,sphhar,atoms,sym,enpara,nococonv,&
ilSpinPr,ilSpin,iSpinV,fmpi,v,vx,input,hub1inp,hub1data,td,ud,alpha_hybrid,lh0,one,v1)
USE m_constants
USE m_intgr, ONLY : intgr3
USE m_genMTBasis
USE m_tlo
USE m_gaunt, ONLY: gaunt1
USE m_types
IMPLICIT NONE
TYPE(t_input),INTENT(IN) :: input
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_enpara),INTENT(IN) :: enpara
TYPE(t_nococonv),INTENT(IN) :: nococonv