Commit da702ccb authored by Gregor Michalicek's avatar Gregor Michalicek

Even more adaptions to the Wannier code provided by JPH

parent faabc5f1
......@@ -19,7 +19,7 @@ c ********************************************************
> lmax,rmsh,jri,pos,ngopr,ntypsy,nvd,
> omtil,amat,bmat,odi,ods,nlod,llod,nlo,llo,
> ff,gg,flo,acof,bcof,ccof,z,
> nv,k1,k2,k3,lmd,nbasfcn,
> nv,k1,k2,k3,lmd,nbasfcn,l_ss,qss,jspin,addnoco,
< xdnout)
c
USE m_types, ONLY : od_inp, od_sym
......@@ -32,8 +32,9 @@ C .. Scalar Arguments ..
INTEGER, INTENT (IN) :: lmaxd,jmtd,ntypd,natd,nmzd
INTEGER, INTENT (IN) :: iflag,n,na,iv,lmd,nv,nvd,nbasfcn
INTEGER, INTENT (IN) :: nq3,nvac,nmz,nmzxy,nq2,nop,nop2
LOGICAL, INTENT (IN) :: invs
REAL, INTENT (IN) :: z1,delz,omtil,bkpt(3)
LOGICAL, INTENT (IN) :: invs,l_ss
REAL, INTENT (IN) :: z1,delz,omtil,bkpt(3),qss(3)
INTEGER, INTENT (IN) :: jspin,addnoco
COMPLEX, INTENT (OUT):: xdnout
c-odim
TYPE (od_inp), INTENT (IN) :: odi
......@@ -69,6 +70,7 @@ C ..
C .. Local Arrays ..
COMPLEX sf2(n2d),sf3(n3d),ylm((lmaxd+1)**2)
REAL rcc(3),x(3),rcc2(3)
REAL bqpt(3)
C ..
C .. Intrinsic Functions ..
......@@ -80,7 +82,18 @@ C ..
c..define the factor e^{-ikr}
rcc2=matmul(bmat,p)/tpi_const
arg = -tpi*( bkpt(1)*rcc2(1) + bkpt(2)*rcc2(2) + bkpt(3)*rcc2(3) )
bqpt = 0.0
! if(l_ss.and.jspin.eq.1) then
! bqpt = -qss/2.0
! elseif(l_ss.and.jspin.eq.2) then
! bqpt = +qss/2.0
! endif
arg = -tpi*( (bkpt(1)+bqpt(1))*rcc2(1)
> + (bkpt(2)+bqpt(2))*rcc2(2)
> + (bkpt(3)+bqpt(3))*rcc2(3) )
arg1 = tpi*( bkpt(1)*rcc2(1) + bkpt(2)*rcc2(2) + bkpt(3)*rcc2(3) )
const2 = cmplx(cos(arg),sin(arg))
const3 = cmplx(cos(arg1),sin(arg1))
......@@ -100,7 +113,7 @@ c write (6,*) 'z(k)=', z(k)
arg = tpi*((k1(k))*rcc(1)+
+ (k2(k))*rcc(2)+
+ (k3(k))*rcc(3) )
xdnout = xdnout + z(k)*cmplx(cos(arg),sin(arg))*const
xdnout = xdnout + z(k+addnoco)*cmplx(cos(arg),sin(arg))*const
IF (((abs(p(1)-2.2).le.0.0001).and.(abs(p(2)).le.0.0001))
& .or.((abs(p(2)-2.2).le.0.0001).and.(abs(p(1)).le.0.0001)))then
c write (6,*) 'p(i)=',p(1:2)
......
......@@ -22,7 +22,7 @@ c****************************************************************
> l_bzsym,l_byindex,l_bynumber,l_byenergy,
> irreduc,odi,band_min,band_max,numbands,
> e1s,e2s,ef,l_paulimag,nkpt,
< nbnd,kpoints)
< nbnd,kpoints,l_gwf,iqpt)
use m_cdnread, only:cdn_read
use m_types,only:od_inp
......@@ -67,6 +67,12 @@ c****************************************************************
complex, allocatable :: paulimat(:,:,:,:)
real :: sum1,sum2,sum3
integer :: num_bands,err
! BEGIN QPOINTS
LOGICAL, INTENT(IN) :: l_gwf
INTEGER, INTENT(IN) :: iqpt
CHARACTER(len=12) :: fname
! END QPOINTS
real,parameter :: hartree=27.21138505
#if ( !defined(CPP_INVERSION) || defined(CPP_SOC) )
......@@ -74,22 +80,29 @@ c****************************************************************
#else
REAL :: z(nbasfcn,neigd)
#endif
! WRITE(*,*)'min',band_min,'max',band_max,'num',numbands
! WRITE(*,*)'wann_rw_eig: neigd',neigd
n_start=1
n_end=neigd
IF (l_gwf) THEN
WRITE(fname,'("_",i4.4,".eig")')iqpt
fname=spin12(jspin)//trim(fname)
ELSE
fname=spin12(jspin)//'.eig'
ENDIF
! WRITE(*,*)'wann_rw_eig.F: writing eig file ',fname
inquire (file=spin12(jspin)//'.eig',exist=l_eig)
inquire (file=fname,exist=l_eig)
if (l_eig) then
open (306,file=spin12(jspin)//'.eig',
open (306,file=fname,
& form='formatted',status='old')
rewind (306)
else
open (306,file=spin12(jspin)//'.eig',
open (306,file=fname,
& form='formatted',status='new')
endif!l_eig
write(*,*)"nkpt=",nkpt
if(l_paulimag.and..false.)then
num_bands=band_max-band_min+1
allocate( paulimat(num_bands,num_bands,3,nkpt),stat=err )
......@@ -125,10 +138,12 @@ c****************************************************************
> 0,isize,kptibz,jspin,nbasfcn,nlotot,
> l_ss,l_noco,nrec,irecl,
< nmat,nv,ello,evdu,epar,kveclo,
< k1,k2,k3,bkpt,wk,nbands,eig,z,cp_time)
< k1,k2,k3,bkpt,wk,nbands,eig,z,cp_time,66,
> l_gwf,iqpt)
nkbnd = 0
! WRITE(*,*)'wann_rw_eig: nbands',nbands
do i = 1,nbands
if((eig(i).ge.e1s .and. nkbnd.lt.numbands.and.l_bynumber).or.
& (eig(i).ge.e1s.and.eig(i).le.e2s.and.l_byenergy ).or.
......@@ -146,7 +161,8 @@ c****************************************************************
endif
endif
endif
enddo
enddo
! WRITE(*,*)'wann_rw_eig: nkbnd',nkbnd,'nbnd',nbnd
if (nkbnd.ge.nbnd) nbnd = nkbnd
if (odi%d1) kpoints(ikpt) = bkpt(3)
......@@ -170,7 +186,8 @@ c****************************************************************
> 0,isize,kptibz,jspin,nbasfcn,nlotot,
> l_ss,l_noco,nrec,irecl,
< nmat,nv,ello,evdu,epar,kveclo,
< k1,k2,k3,bkpt,wk,nbands,eig,z,cp_time)
< k1,k2,k3,bkpt,wk,nbands,eig,z,cp_time,66,
> l_gwf,iqpt)
nkbnd = 0
......@@ -195,20 +212,21 @@ c*********************************************************************
> irank,isize,kptibz,jspin,nbasfcn,nlotot,
> l_ss,l_noco,nrec,irecl,
< nmat,nv,ello,evdu,epar,kveclo,
< k1,k2,k3,bkpt,wk,nbands,eig,z,cp_time)
< k1,k2,k3,bkpt,wk,nbands,eig,z,cp_time,funit,
> l_gwf,iqpt)
use m_cdnread, only:cdn_read
implicit none
INTEGER, INTENT (IN) :: irank,isize,kptibz,nbasfcn,neigd,nlotot
INTEGER, INTENT (IN) :: nrec,nvd,jspd,jspin
INTEGER, INTENT (IN) :: lmaxd,ntypd,nlod
INTEGER, INTENT (IN) :: lmaxd,ntypd,nlod,iqpt
LOGICAL, INTENT (IN) :: l_ss,l_noco,l_gwf
#ifndef CPP_MPI
integer(8), intent (in) :: irecl
#else
include 'mpif.h'
INTEGER(KIND=MPI_OFFSET_KIND),intent(in) :: irecl
#endif
LOGICAL, INTENT (IN) :: l_ss,l_noco
INTEGER, INTENT (OUT) :: nbands,nmat
REAL, INTENT (OUT) :: wk
......@@ -225,7 +243,9 @@ c*********************************************************************
REAL, INTENT (OUT) :: z(nbasfcn,neigd)
#endif
character(len=12) :: fname
INTEGER, INTENT(IN) :: funit
character(len=12) :: fname,fend
integer :: aoff,d1,d10
integer :: cpu_index
logical :: l_eig
......@@ -234,6 +254,8 @@ c*********************************************************************
n_start=1
n_end=neigd
fend=''
if(l_gwf) write(fend,'("_",i4.4)')iqpt
do cpu_index=0,isize-1
......@@ -242,25 +264,28 @@ c*********************************************************************
if(cpu_index.eq.irank)then
#if 1==1
call judft_error("NOT IMPLEMENTED")
#else
call cdn_read(
> eig_id,nvd,jspd,
> irank,isize,kptibz,jspin,nbasfcn,
> l_noco,neigd,
> lmaxd,ntypd,nlod,neigd,nvd,jspd,
> irank,isize,kptibz,jspin,nbasfcn,nlotot,
> l_ss,l_noco,nrec,kptibz,funit,neigd,
> n_start,n_end,
< nmat,nv,ello,evdu,epar,kveclo,
< k1,k2,k3,bkpt,wk,nbands,eig,zmat)
#endif
< k1,k2,k3,bkpt,wk,nbands,eig,z,cp_time)
else !eigenvalues are stored in the file of another process
#if 1==2
d10=int((cpu_index+0.5)/10)
d1=mod(cpu_index,10)
aoff=iachar('1')-1
if(d10.ge.10)aoff=iachar('7')
fname='eig'//achar(d10+aoff)//achar(d1+iachar('1')-1)
fname=trim(fname)//trim(fend)
l_eig=.false.
!write(*,*)fname
c if(irank.eq.0) write(*,*) 'opening file=',fname
inquire(file=fname,exist=l_eig)
IF(.NOT.l_eig) CALL juDFT_error("one of eig files not found"
+ ,calledby ="wann_rw_eig")
......@@ -275,37 +300,35 @@ c*********************************************************************
< k1,k2,k3,bkpt,wk,nbands,eig,z,cp_time)
close(666)
#endif
endif !cpu_index==irank
endif !find process
enddo !cpu_index
enddo !cpu_index
end subroutine wann_mpi_read_eig
c*****************************************************************
c read in eig file
c*****************************************************************
subroutine wann_read_eig(
subroutine wann_read_eig(
> ntapwf,
> lmaxd,ntypd,nlod,neigd,nvd,jspd,
> irank,isize,kptibz,jspin,nbasfcn,nlotot,
> l_ss,l_noco,nrec,irecl,
< nmat,nv,ello,evdu,epar,kveclo,
< k1,k2,k3,bkpt,wk,nbands,eig,z,cp_time)
< k1,k2,k3,bkpt,wk,nbands,eig,z,cp_time,funit,
> l_gwf,iqpt)
use m_cdnread, only:cdn_read
implicit none
integer, intent (in) :: ntapwf
INTEGER, INTENT (IN) :: irank,isize,kptibz,nbasfcn,neigd,nlotot
INTEGER, INTENT (IN) :: nrec,nvd,jspd,jspin
INTEGER, INTENT (IN) :: lmaxd,ntypd,nlod
INTEGER, INTENT (IN) :: lmaxd,ntypd,nlod,iqpt
LOGICAL, INTENT (IN) :: l_ss,l_noco,l_gwf
#ifndef CPP_MPI
integer(8), intent (in) :: irecl
#else
include 'mpif.h'
INTEGER(KIND=MPI_OFFSET_KIND),intent(in) :: irecl
#endif
LOGICAL, INTENT (IN) :: l_ss,l_noco
INTEGER, INTENT (OUT) :: nbands,nmat
REAL, INTENT (OUT) :: wk
......@@ -322,6 +345,8 @@ c*****************************************************************
REAL, INTENT (OUT) :: z(nbasfcn,neigd)
#endif
INTEGER, INTENT(IN) :: funit !QPOINTS
integer :: n_start,n_end
n_start=1
......@@ -329,37 +354,24 @@ c*****************************************************************
#if( defined(CPP_MPI) && !defined(CPP_HDF) )
! MPI without HDF
#ifdef CPP_HYBEIG
c-----if hybrid functionals are used, and mpi is defined, mpi-io is
c-----used and all processors use the same eig-file
call cdn_read(
> lmaxd,ntypd,nlod,neigd,nvd,jspd,
> irank,isize,kptibz,jspin,nbasfcn,nlotot,
> l_ss,l_noco,nrec,kptibz,ntapwf,neigd,
> n_start,n_end,
< nmat,nv,ello,evdu,epar,kveclo,
< k1,k2,k3,bkpt,wk,nbands,eig,z,cp_time)
#else
call wann_mpi_read_eig(
call wann_mpi_read_eig(
> lmaxd,ntypd,nlod,neigd,nvd,jspd,
> irank,isize,kptibz,jspin,nbasfcn,nlotot,
> l_ss,l_noco,nrec,irecl,
< nmat,nv,ello,evdu,epar,kveclo,
< k1,k2,k3,bkpt,wk,nbands,eig,z,cp_time)
#endif
< k1,k2,k3,bkpt,wk,nbands,eig,z,cp_time,funit,
> l_gwf,iqpt)
#else
!MPI with HDF or no MPI
#if 1==1
call judft_error("NOT IMPLEMENTED")
#else
call cdn_read(
> lmaxd,ntypd,nlod,neigd,nvd,jspd,
> 0,isize,kptibz,jspin,nbasfcn,nlotot,
> l_ss,l_noco,nrec,kptibz,66,
> l_ss,l_noco,nrec,kptibz,funit,
> neigd,n_start,n_end,
< nmat,nv,ello,evdu,epar,kveclo,
< k1,k2,k3,bkpt,wk,nbands,eig,z,cp_time)
#endif
#endif
END SUBROUTINE wann_read_eig
END MODULE m_wann_rw_eig
......@@ -15,10 +15,11 @@ c Frank Freimuth, October 2006
c***********************************************************************
> llod,nntot,kdiff,lmax,
> ntype,ntypd,bbmat,bmat,
> nlod,nlo,llo,flo,f,g,jri,rmsh,dx,jmtd,
> nlod,nlo,llo,flo,flo_b,f,f_b,g,g_b,
> jri,rmsh,dx,jmtd,
> lmaxd,lmd,
< ujug,ujdg,djug,djdg,ujulog,djulog,
< ulojug,ulojdg,ulojulog)
< ulojug,ulojdg,ulojulog,l_q,sign_q)
use m_constants, only : pimach
use m_matmul , only : matmul3,matmul3r
use m_sphbes
......@@ -26,7 +27,7 @@ c***********************************************************************
use m_types, only : od_inp, od_sym
use m_intgr, only : intgr3
use m_gaunt, only: gaunt1
IMPLICIT NONE
integer, intent (in) :: llod
INTEGER, INTENT (IN) :: nntot,ntype,ntypd
......@@ -40,11 +41,17 @@ c***********************************************************************
integer, intent (in) :: nlo(:) !(ntypd)
integer, intent (in) :: llo(:,:) !(nlod,ntypd)
real, intent (in) :: f(:,:,:,0:) !(ntypd,jmtd,2,0:lmaxd)
real, intent (in) :: f_b(:,:,:,0:) !(ntypd,jmtd,2,0:lmaxd)
real, intent (in) :: g(:,:,:,0:) !(ntypd,jmtd,2,0:lmaxd)
real, intent (in) :: g_b(:,:,:,0:) !(ntypd,jmtd,2,0:lmaxd)
real, intent (in) :: flo(:,:,:,:) !(ntypd,jmtd,2,nlod)
real, intent (in) :: flo_b(:,:,:,:) !(ntypd,jmtd,2,nlod)
real, intent (in) :: rmsh(:,:) !(jmtd,ntypd)
real, intent (in) :: dx(:) !(ntypd)
logical, intent (in) :: l_q ! if true, we deal with q points
integer, intent (in) :: sign_q ! if we deal with q points, we might pick up an additional sign
complex, intent (out) :: ujug(0:,0:,:,:) !(0:lmd,0:lmd,1:ntype,1:nntot)
complex, intent (out) :: ujdg(0:,0:,:,:) !(0:lmd,0:lmd,1:ntype,1:nntot)
complex, intent (out) :: djug(0:,0:,:,:) !(0:lmd,0:lmd,1:ntype,1:nntot)
......@@ -56,7 +63,9 @@ c***********************************************************************
complex, intent (out) :: ulojulog(:,-llod:,:,-llod:,:,:) !(1:nlod,-llod:llod,1:nlod,-llod:llod,1:ntype,1:nntot)
real, allocatable :: djd(:,:,:),ujd(:,:,:),uju(:,:,:)
real, allocatable :: dju(:,:,:)
real, allocatable :: ujulo(:,:,:),djulo(:,:,:),ulojulo(:,:,:)
real, allocatable :: uloju(:,:,:),ulojd(:,:,:)
integer :: ikpt_b,i,lwn,n,lpp,lop,lo,l,lp
integer :: lmini,lmaxi,m,mp,llpp,mpp
integer :: lmpp,lminp,lmaxp,lm,lpmp
......@@ -64,15 +73,19 @@ c***********************************************************************
real :: jj(0:lmaxd,jmtd),x(jmtd)
real :: bkrot(3)
complex :: ylmpp((lmaxd+1)**2),factor,ic
complex :: factor2
ic = cmplx(0.,1.)
allocate( djd(0:lmaxd,0:lmaxd,0:lmaxd) )
allocate( ujd(0:lmaxd,0:lmaxd,0:lmaxd) )
allocate( dju(0:lmaxd,0:lmaxd,0:lmaxd) )
allocate( uju(0:lmaxd,0:lmaxd,0:lmaxd) )
allocate( ujulo(nlod,0:lmaxd,0:lmaxd) )
allocate( djulo(nlod,0:lmaxd,0:lmaxd) )
allocate( uloju(nlod,0:lmaxd,0:lmaxd) )
allocate( ulojd(nlod,0:lmaxd,0:lmaxd) )
allocate( ulojulo(nlod,nlod,0:lmaxd) )
......@@ -87,10 +100,12 @@ c***********************************************************************
ulojdg = cmplx(0.0,0.0)
ulojulog = cmplx(0.0,0.0)
do ikpt_b=1,nntot
bpt(:)=kdiff(:,ikpt_b)
rk = sqrt(dot_product(bpt,matmul(bbmat,bpt)))
!write(*,*)'ujugaunt rk',rk
do n=1,ntype
lwn = lmax(n)
c...generate the j_lpp(br) on the radial grid
......@@ -116,8 +131,8 @@ c..the gaunt conditions
ulojulo(lo,lop,lpp) = 0.
else
do i = 1,jri(n)
x(i) = ( flo(n,i,1,lo)*flo(n,i,1,lop)+
+ flo(n,i,2,lo)*flo(n,i,2,lop) )*jj(lpp,i)
x(i) = ( flo(n,i,1,lo)*flo_b(n,i,1,lop)+
+ flo(n,i,2,lo)*flo_b(n,i,2,lop) )*jj(lpp,i)
enddo
call intgr3(x,rmsh(1:,n),dx(n),jri(n),ulojulo(lo,lop,lpp))
endif
......@@ -136,23 +151,30 @@ c..gaunt conditions
+ (lpp.gt.lmaxi)) then
uju(l,lp,lpp) = 0.
ujd(l,lp,lpp) = 0.
dju(l,lp,lpp) = 0.
djd(l,lp,lpp) = 0.
else
do i = 1,jri(n)
x(i) = ( f(n,i,1,l)*f(n,i,1,lp)+
+ f(n,i,2,l)*f(n,i,2,lp) )*jj(lpp,i)
x(i) = ( f(n,i,1,l)*f_b(n,i,1,lp)+
+ f(n,i,2,l)*f_b(n,i,2,lp) )*jj(lpp,i)
enddo
call intgr3(x,rmsh(1:,n),dx(n),jri(n),uju(l,lp,lpp))
do i = 1,jri(n)
x(i) = ( f(n,i,1,l)*g(n,i,1,lp)+
+ f(n,i,2,l)*g(n,i,2,lp) )*jj(lpp,i)
x(i) = ( f(n,i,1,l)*g_b(n,i,1,lp)+
+ f(n,i,2,l)*g_b(n,i,2,lp) )*jj(lpp,i)
enddo
call intgr3(x,rmsh(1:,n),dx(n),jri(n),ujd(l,lp,lpp))
do i = 1,jri(n)
x(i) = ( g(n,i,1,l)*g(n,i,1,lp)+
+ g(n,i,2,l)*g(n,i,2,lp) )*jj(lpp,i)
x(i) = ( g(n,i,1,l)*f_b(n,i,1,lp)+
+ g(n,i,2,l)*f_b(n,i,2,lp) )*jj(lpp,i)
enddo
call intgr3(x,rmsh(1:,n),dx(n),jri(n),dju(l,lp,lpp))
do i = 1,jri(n)
x(i) = ( g(n,i,1,l)*g_b(n,i,1,lp)+
+ g(n,i,2,l)*g_b(n,i,2,lp) )*jj(lpp,i)
enddo
call intgr3(x,rmsh(1:,n),dx(n),jri(n),djd(l,lp,lpp))
endif
......@@ -171,17 +193,29 @@ c..gaunt conditions
+ (lpp.gt.lmaxi)) then
ujulo(lo,lp,lpp) = 0.
djulo(lo,lp,lpp) = 0.
uloju(lo,lp,lpp) = 0.
ulojd(lo,lp,lpp) = 0.
else
do i = 1,jri(n)
x(i) = ( flo(n,i,1,lo)*f(n,i,1,lp)+
+ flo(n,i,2,lo)*f(n,i,2,lp) )*jj(lpp,i)
x(i) = ( flo(n,i,1,lo)*f_b(n,i,1,lp)+
+ flo(n,i,2,lo)*f_b(n,i,2,lp) )*jj(lpp,i)
enddo
call intgr3(x,rmsh(1:,n),dx(n),jri(n),ujulo(lo,lp,lpp))
do i = 1,jri(n)
x(i) = ( flo(n,i,1,lo)*g(n,i,1,lp)+
+ flo(n,i,2,lo)*g(n,i,2,lp) )*jj(lpp,i)
x(i) = ( flo(n,i,1,lo)*g_b(n,i,1,lp)+
+ flo(n,i,2,lo)*g_b(n,i,2,lp) )*jj(lpp,i)
enddo
call intgr3(x,rmsh(1:,n),dx(n),jri(n),djulo(lo,lp,lpp))
do i = 1,jri(n)
x(i) = ( flo_b(n,i,1,lo)*f(n,i,1,lp)+
+ flo_b(n,i,2,lo)*f(n,i,2,lp) )*jj(lpp,i)
enddo
call intgr3(x,rmsh(1:,n),dx(n),jri(n),uloju(lo,lp,lpp))
do i = 1,jri(n)
x(i) = ( flo_b(n,i,1,lo)*g(n,i,1,lp)+
+ flo_b(n,i,2,lo)*g(n,i,2,lp) )*jj(lpp,i)
enddo
call intgr3(x,rmsh(1:,n),dx(n),jri(n),ulojd(lo,lp,lpp))
endif
enddo !lo
endif ! local orbitals
......@@ -191,7 +225,6 @@ c********************************************************************
c multiply with gaunt-coefficient (apw-apw)
c********************************************************************
bkrot=matmul(bpt,bmat)
!call cotra3(bpt,bkrot,bmat)
call ylm4(lwn,bkrot,ylmpp)
do l = 0,lwn
do m = -l,l
......@@ -201,21 +234,31 @@ c********************************************************************
lpmp=lp*(lp+1)+mp
do lpp = 0,lwn
llpp = lpp*(lpp+1)
mpp = mp - m
lmpp = llpp + mpp
lmini = abs(l-lpp)
lmaxi = l+lpp
if ((lmini.le.lp).and.(lp.le.lmaxi).and.
& (mod(l+lp+lpp,2).eq.0).and.(abs(mpp).LE.lpp))then
factor=conjg(ylmpp(lmpp+1))*(ic**(l+lpp-lp))*
* gaunt1(lp,lpp,l,mp,mpp,m,lmaxd)
c if(factor.ne.0 .and. uju(l,lp,lpp).ne.0)
c > write(*,*)lpp,lp,l,mp,m
if(l_q) then
factor=(sign_q**lpp)*factor ! additional sign for q points
endif
ujug(lpmp,lm,n,ikpt_b)=ujug(lpmp,lm,n,ikpt_b)+
+ factor*uju(l,lp,lpp)
ujdg(lpmp,lm,n,ikpt_b)=ujdg(lpmp,lm,n,ikpt_b)+
+ factor*ujd(l,lp,lpp)
djug(lpmp,lm,n,ikpt_b)=djug(lpmp,lm,n,ikpt_b)+
+ factor*ujd(lp,l,lpp)
+ factor*dju(l,lp,lpp)
djdg(lpmp,lm,n,ikpt_b)=djdg(lpmp,lm,n,ikpt_b)+
+ factor*djd(l,lp,lpp)
......@@ -248,6 +291,9 @@ c******************************************************************
lmpp = llpp + mpp
factor=conjg(ylmpp(lmpp+1))*(ic**(l+lpp-lp))*
* gaunt1(lp,lpp,l,mp,mpp,m,lmaxd)
if(l_q) then
factor=(sign_q**lpp)*factor ! additional sign for q points
endif
ujulog(lpmp,lo,m,n,ikpt_b)=ujulog(lpmp,lo,m,n,ikpt_b)+
+ factor*ujulo(lo,lp,lpp)
......@@ -261,11 +307,14 @@ c******************************************************************
lmpp = llpp + mpp
factor=conjg(ylmpp(lmpp+1))*(ic**(lp+lpp-l))*
* gaunt1(l,lpp,lp,m,mpp,mp,lmaxd)
if(l_q) then
factor=(sign_q**lpp)*factor
endif
ulojug(lpmp,lo,m,n,ikpt_b)=ulojug(lpmp,lo,m,n,ikpt_b)+
+ factor*ujulo(lo,lp,lpp)
+ factor*uloju(lo,lp,lpp)
ulojdg(lpmp,lo,m,n,ikpt_b)=ulojdg(lpmp,lo,m,n,ikpt_b)+
+ factor*djulo(lo,lp,lpp)
+ factor*ulojd(lo,lp,lpp)
endif
enddo ! lpp
enddo ! mp
......@@ -291,11 +340,19 @@ c*************************************************************
lmaxi = l+lpp
if ((lmini.le.lp).and.(lp.le.lmaxi).and.
& (mod(l+lp+lpp,2).eq.0).and.(abs(mpp).LE.lpp))then
ulojulog(lop,mp,lo,m,n,ikpt_b)=
factor= conjg(ylmpp(lmpp+1))*(ic**(l+lpp-lp))*
* gaunt1(lp,lpp,l,mp,mpp,m,lmaxd)
if(l_q) then
factor=(sign_q**lpp)*factor ! additional sign for q points
endif
ulojulog(lop,mp,lo,m,n,ikpt_b)=
= ulojulog(lop,mp,lo,m,n,ikpt_b)+
+ conjg(ylmpp(lmpp+1))*(ic**(l+lpp-lp))*
* gaunt1(lp,lpp,l,mp,mpp,m,lmaxd)*
* ulojulo(lo,lop,lpp)
+ ulojulo(lo,lop,lpp)*factor
! + conjg(ylmpp(lmpp+1))*(ic**(l+lpp-lp))*
! * gaunt1(lp,lpp,l,mp,mpp,m,lmaxd)*
! * ulojulo(lo,lop,lpp)
endif
enddo ! lpp
enddo ! mp lop
......@@ -306,6 +363,7 @@ c*************************************************************
enddo !ntype
enddo !ikpt_b
deallocate(djd,ujd,uju,ujulo,djulo,ulojulo)
deallocate(dju,uloju,ulojd)
end subroutine wann_ujugaunt
end module m_wann_ujugaunt
......@@ -114,7 +114,7 @@ c****************************************************************************
cccccccccccccccccc local variables cccccccccccccccccccc
type(t_wann) :: wann
integer :: lmd,nlotot,n,nmat,iter,ikpt,ikpt_b,nmat_b
integer :: addnoco,funbas,loplod
integer :: addnoco,addnoco2,funbas,loplod
integer :: noccbd,noccbd_b,nn,nkpts,i,jspin,j,l,i_rec,m,nwf,nwfp
integer :: jsp_start,jsp_end,nrec,nrec1,nbands,nbands_b
integer :: nodeu,noded,n_size,na,n_rank,nbnd,numbands
......@@ -335,8 +335,8 @@ c**************************************************************
open(713,file='kptsmap')
do i=1,fullnkpts
read(713,*)kpt,irreduc(i),mapkoper(i),shiftkpt(:,i)
IF(kpt/=i) CALL juDFT_error("kpt.ne.i",calledby ="wann_updown"
+ )
IF(kpt/=i) CALL juDFT_error("kpt.ne.i",
+ calledby ="wann_updown")
if(l_p0)print*,i,irreduc(i),mapkoper(i)
enddo
close(713)
......@@ -431,8 +431,8 @@ c..reading the proj.1 / proj.2 / proj file
rewind (203)
close (203)
else
CALL juDFT_error("no proj/proj.1/proj.2",calledby ="wann_updown"
+ )
CALL juDFT_error("no proj/proj.1/proj.2",
+ calledby ="wann_updown")
endif