Hi,
Thank you very much for this bug report. I can confirm the problem.
It was a nasty bug, occurring only if you have
. non-spinpolarized spin-orbit calculations AND
. more than one atom/cell (the first atom is always ok).
The calculated partial charges (case.qtl) are all wrong except for the
first atom (not only the relativistic decomposition into 5/2 and 7/2,
..., but also the regular partial charges using "x qtl".
In addition there was a factor of two in the normalization, so that the
total DOS could be less than the partial one.
I include a few subroutines for SRC_qtl which should fix the problem.
Please verify.
Regards
Peter Blaha
On 06/19/2017 11:06 AM, Jindrich Kolorenc wrote:
Dear Wien2k developers,
I am trying to decompose the U 5f density of states in UGa2 to j=5/2 and
j=7/2 components and it appears that I am hitting a bug in the QTL
program. In the non-magnetic + SO calculation, the 5f DOS has two peaks,
the lower should be mostly j=5/2, the upper j=7/2. An earlier
calculation can be found in Fig.3 of Divis et al.,
https://doi.org/10.1103/PhysRevB.53.9658
The figure
http://www.fzu.cz/~kolorenc/wien2k_UGa2_bug/UGa2_f52_f72.png
shows what I get with Wien2k, version 16.1 (and 14.2 gives the same).
One can see that
1) it differs from the results of Divis at al.
2) the result depends on the direction of magnetization specified in
case.inso (001 or 100), which I believe it really should not when there
is no magnetization at all
3) the result for the direction 100 is clearly incorrect (the lower peak
certainly is not dominantly 7/2)
4) QTL gives long tail of the 5f DOS up to high energies whereas the 5f
DOS calculated after "lapw2 -so -qtl" does not have this tail.
Given the above, I think the problem is somewhere in interpreting the
lattice symmetry in QTL (here P6/mmm, group 191). My case.struct is
attached, the other input files were essentially straight from init_lapw
and initso_lapw. I use symmetrization in QTL just to be on the safe
side. I have also tried to run the calculations with the "ferromagnetic"
case.struct generated by initso_lapw (it reduces the number of symmetry
operations in the 100 case) but there seems to be no difference.
Any help leading toward a fix (or a workaround) would be really
appreciated.
Best regards,
Jindrich Kolorenc
_______________________________________________
Wien mailing list
Wien@zeus.theochem.tuwien.ac.at
http://zeus.theochem.tuwien.ac.at/mailman/listinfo/wien
SEARCH the MAILING-LIST at:
http://www.mail-archive.com/wien@zeus.theochem.tuwien.ac.at/index.html
--
P.Blaha
--------------------------------------------------------------------------
Peter BLAHA, Inst.f. Materials Chemistry, TU Vienna, A-1060 Vienna
Phone: +43-1-58801-165300 FAX: +43-1-58801-165982
Email: bl...@theochem.tuwien.ac.at WIEN2k: http://www.wien2k.at
WWW: http://www.imc.tuwien.ac.at/TC_Blaha
--------------------------------------------------------------------------
SUBROUTINE ATPAR (JATOM,LATOM,itape,jtape,is,ISPIN,popmat)
USE param
USE struct
use com
IMPLICIT REAL*8 (A-H,O-Z)
integer popmat
LOGICAL loor(0:lomax),lapw(0:lmax2),rlo(1:nloat, &
0:lomax)
DIMENSION emist(0:lomax,nloat),E(0:LMAX2),elo(0:LOMAX,nloat), &
pei(0:lmax2),ilo(0:lomax),vr(nrad)
COMMON /UHELP/ A(NRAD),B(NRAD),AE(NRAD),BE(NRAD)
COMMON /ATSPDT/ P(0:LMAX2,2,nrf),DP(0:LMAX2,2,nrf)
COMMON /RADFU/ RF1(NRAD,0:LMAX2,2,nrf),RF2(NRAD,0:LMAX2,2,nrf)
COMMON /RINTEG/ RI_MAT(0:lmax2,0:lmax2,nrf,nrf,2,2)
common /loabc/ alo(0:lomax,2,nloat,nrf)
common /lolog/ nlo,nlov,nlon,loor,ilo,lapw
save vr
!---------------------------------------------------------------------
2022 FORMAT(3X,4E19.12)
!.....READ TOTAL SPHERICAL POTENTIAL V(0,0) OF TAPEjtape=VSP
! NORM OF V(0,0)=V00(R)*R/SQRT(4.D0*PI)
if(iso2.eq.1.and.is.eq.1) then
READ(jtape,1980)
READ(jtape,2000)IDUMMY
READ(jtape,2031)
READ(jtape,2022)(VR(J),J=1,JRJ(JATOM))
READ(jtape,2031)
READ(jtape,2030)
DO 115 J=1,JRJ(JATOM)
115 VR(J)=VR(J)/2.0D0
endif
if (ispin.eq.2) then
if(is.eq.1)then
write(6,*)'ATPAR for spin up'
else
write(6,*)'ATPAR for spin down'
endif
end if
nlo=0
nlov=0
nlon=0
do i=0,lomax
ilo(i)=0
end do
DO 15 I=1,JATOM
READ(itape)E
READ(itape)elo
IF(i.EQ.jatom) THEN
DO l=0,lmax2
lapw(l)=.TRUE.
IF(e(l).GT.150.) THEN
e(l)=e(l)-200.d+0
lapw(l)=.FALSE.
ENDIF
ENDDO
ENDIF
rlo(:,:)=.FALSE.
do 15 l = 0,lomax
loor(l)=.FALSE.
DO k=1,nloat
IF (i.EQ.jatom) THEN
IF (elo(l,k).LT.(995.d+0)) THEN
ilo(l)=ilo(l)+1
nlo=nlo+((2*l+1))*mult(i)
IF(.NOT.lapw(l).AND.k.EQ.1) GOTO 666
IF(k.EQ.nloat) THEN
rlo(ilo(l),l)=.TRUE.
GOTO 666
ENDIF
loor(l)=.TRUE.
666 CONTINUE
ENDIF
ELSE
IF (elo(l,k).LT.(995.d+0)) nlov=nlov+((2*l+1))*mult(i)
ENDIF
ENDDO
15 continue
IF(JATOM.EQ.NAT) GOTO 30
DO 20 I=JATOM+1,NAT
READ(itape) EMIST
READ(itape) EMIST
do 20 l=0,lomax
do 20 k=1,nloat
if (emist(l,k).lt.(995.0D+0)) &
nlon=nlon+((2*l+1))*mult(i)
20 continue
30 CONTINUE
WRITE(6,7) ANAME(JATOM)
WRITE(6,5) E
WRITE(6,14)
!
DO 70 l=0,LMAX2
DELE=2.0D-3
DELEI=0.25D0/DELE
FL=L
EI=E(l)/2.0d0
! CALCULATE ENERGY-DERIVATIVE BY FINITE DIFFERENCE
! DELE IS THE UPWARD AND DOWNWARD ENERGY SHIFT IN HARTREES
!
E1=EI-DELE
CALL OUTWIN(REL,VR,r0(JATOM),DX(JATOM),jrj(JATOM),E1, &
FL,UVB,DUVB,NODEL,zz(jatom))
CALL rint13(A,B,A,B,OVLP,JATOM)
TRX=1.0D0/SQRT(OVLP)
IMAX=jrj(JATOM)
DO 45 M=1,IMAX
AE(M)=TRX*A(M)
BE(M)=TRX*B(M)
45 CONTINUE
UVB=TRX*UVB
DUVB=TRX*DUVB
E1=EI+DELE
CALL OUTWIN(REL,VR,r0(JATOM),DX(JATOM),jrj(JATOM),E1, &
FL,UVE,DUVE,NODE,zz(jatom))
CALL rint13(A,B,A,B,OVLP,JATOM)
TRX=1.0d0/SQRT(OVLP)
UVE=DELEI*(TRX*UVE-UVB)
DUVE=DELEI*(TRX*DUVE-DUVB)
IMAX=jrj(JATOM)
DO 50 M=1,IMAX
AE(M)=DELEI*(TRX*A(M)-AE(M))
BE(M)=DELEI*(TRX*B(M)-BE(M))
50 CONTINUE
!
! CALCULATE FUNCTION AT EI
!
CALL OUTWIN(REL,VR,r0(JATOM),DX(JATOM),jrj(JATOM),EI, &
FL,UV,DUV,NODES,zz(jatom))
CALL rint13(A,B,A,B,OVLP,JATOM)
TRX=1.0d0/SQRT(OVLP)
P(l,is,1)=TRX*UV
DP(l,is,1)=TRX*DUV
IMAX=jrj(JATOM)
DO 60 M=1,IMAX
A(M)=TRX*A(M)
60 B(M)=TRX*B(M)
!
! INSURE ORTHOGONALIZATION
!
CALL rint13(A,B,AE,BE,CROSS,JATOM)
TRY=-CROSS
IMAX=jrj(JATOM)
DO 55 M=1,IMAX
AE(M)=(AE(M)+TRY*A(M))
55 BE(M)=(BE(M)+TRY*B(M))
IMAX=jrj(JATOM)
DO 80 I=1,IMAX
RF1(I,l,is,1)=A(I)
RF2(I,l,is,1)=B(I)
RF1(I,l,is,2)=AE(I)
RF2(I,l,is,2)=BE(I)
80 CONTINUE
P(l,is,2)=UVE+TRY*P(l,is,1)
DP(l,is,2)=DUVE+TRY*DP(l,is,1)
CALL rint13(AE,BE,AE,BE,PEI(l),JATOM)
70 WRITE(6,8) L,P(l,is,1),DP(l,is,1),P(l,is,2),DP(l,is,2)
!
! nun fur lo
!
DO 170 l=0,lomax
irf=2
do 170 jlo=1,ilo(l)
if ((.not.lapw(l)).and.(jlo.eq.1)) goto 180
irf=irf+1
DELE=2.0D-3
DELEI=0.25D0/DELE
FL=L
EI=elo(l,jlo)/2.d0
!
! CALCULATE FUNCTION AT EI
IF(rlo(jlo,l)) THEN
ei=elo(l,nloat)/2.d0
kappa=l
CALL diracout(rel,vr,r0(jatom),dx(jatom),jrj(jatom), &
ei,fl,kappa,uv,duv,nodes,zz(jatom))
CALL dergl(a,b,r0(jatom),dx(jatom),jrj(jatom))
DO m = 1, jrj(jatom)
r_m = r0(jatom)*exp(dx(jatom)*(m-1))
b(m) = b(m)*r_m/(2.d0*clight+(elo(l,jlo)- &
2.d0*vr(m)/r_m)/(2.d0*clight))
b(m)=b(m)*clight
ENDDO
ELSE
CALL outwin(rel,vr,r0(jatom),dx(jatom),jrj(jatom), &
ei,fl,uv,duv,nodes,zz(jatom))
ENDIF
!
CALL rint13(A,B,A,B,OVLP,JATOM)
TRX=1.0d0/SQRT(OVLP)
P(l,is,irf)=TRX*UV
DP(l,is,irf)=TRX*DUV
IMAX=jrj(JATOM)
DO 160 M=1,IMAX
rf1(M,l,is,irf)=TRX*A(M)
160 rf2(M,l,is,irf)=TRX*B(M)
CALL rint13(rf1(1,l,is,1),rf2(1,l,is,1), &
rf1(1,l,is,irf),rf2(1,l,is,irf),pi12lo,JATOM)
CALL rint13(rf1(1,l,is,2),rf2(1,l,is,2), &
rf1(1,l,is,irf),rf2(1,l,is,irf),pe12lo,JATOM)
180 continue
call abcd(l,jatom,pei(l),pi12lo,pe12lo,is,jlo,lapw(l))
WRITE(6,*)
170 continue
!... CALCULATION OF RADIAL INTEGRALS
IF (ISPIN.EQ.2.AND.IS.EQ.1) GOTO 469
CALL RADINT(JATOM,ISPIN,popmat)
469 RETURN
2 FORMAT(5E14.7)
3 FORMAT(16I5)
4 FORMAT(I4,E16.7)
5 FORMAT(10X,' ENERGY PARAMETERS ARE',7F7.2)
6 FORMAT(10X,'E(',I2,2H)=,F10.4)
7 FORMAT(/10X,'ATOMIC PARAMETERS FOR ',A10/)
14 FORMAT(/11X,1HL,5X,4HU(R),10X, &
5HU'(R),9X,5HDU/DE,8X,6HDU'/DE,6X,7HNORM-U')
8 FORMAT(10X,I2,5E14.6,5X,3I2)
9 FORMAT(10X,'FOR L=',I2,' CORRECTION=',E14.5,' OVERLAP=',E14.5)
11 FORMAT(7F10.5)
13 FORMAT(////,':POS',i2.2,':',1x,'AT.NR.',I3,2X,'POSITION =', &
3F8.5,2X,'MULTIPLICITY =',I3)
1040 FORMAT(8I2)
1060 FORMAT(//,3X,'NOT EQUIV ATOM ',A10,' LOCAL ROTATION MATRIX')
1070 FORMAT(30X,3F10.5)
1080 FORMAT(13X,'EQUIV ATOM ',I3,3X,'POSITION: ',3F8.3)
1090 FORMAT(30X,3F10.5,5X,F10.5)
1980 FORMAT(3X)
2000 FORMAT(16X,I2//)
2030 FORMAT(///)
2031 FORMAT(/)
END
SUBROUTINE OUTP(ICASE,popmat,so)
USE param
USE struct
USE case
USE com
! use abc
IMPLICIT REAL*8 (A-H,O-Z)
CHARACTER*7 LDUMMY
CHARACTER*10 KNAME
CHARACTER*67 ERRMSG
logical so
integer popmat,naxial(0:3),nhexa(0:3),ncubic(0:3),iL(0:4)
! index to the lms,l'm's' component: indices L1,M1,S1,L2,M2,S2,re/im ... maxl fixed to 5
integer kkla(0:5,-5:5,1:2,0:5,-5:5,1:2,1:2)
DIMENSION qtl1(4,14),qtl2(34),suml1(4)
complex*16 qtle,qtle1(3072)
dimension e(1000)
data naxial/1,2,4,5/,nhexa/1,2,3,5/,ncubic/1,1,2,3/
data suml1/0.,0.,0.,0./
!--------------------------------------------------------------------
ic1=1
write(6,*)'NUMBER OF K-POINTS:',nkpt
nband=100000
do i=1,nkpt
write(6,*)i,nb(i),'bands'
if(nb(i).lt.nband)nband=nb(i)
end do
if(popmat.eq.0)then
!pb iL(lcase)=0
iL=0
do lcase=1,nl(icase)
l=ll(icase,lcase)
if((nm(icase).le.2).or.(nm(icase).eq.6))then
iL(lcase)=(2*l+1)*iso
if(nm(icase).eq.0)iL(lcase)=2 ! |J,MJ> basis, L-1/2,L+1/2 sums
else if(nm(icase).eq.3)then ! axial symmetry
iL(lcase)=naxial(L)*iso
else if(nm(icase).eq.4)then ! hexagonal symmetry
iL(lcase)=nhexa(L)*iso
else if(nm(icase).eq.5)then ! cubic symmetry
iL(lcase)=ncubic(L)*iso
endif
enddo
else
kkL=100 ! write rest of popmat input for TETRA on files 70+icase
do Lcase1=1,nl(icase)
L1=LL(icase,Lcase1)
do ms1=1,iso
do m1=-L1,L1
do Lcase2=1,nL(icase)
L2=LL(icase,Lcase2)
do ms2=1,iso
do m2=-L2,L2
if(ms1.gt.ms2)then
kkL=kkL+1
kkla(l1,m1,ms1,l2,m2,ms2,1) = kkl
kkL=kkL+1
kkla(l1,m1,ms1,l2,m2,ms2,2) = kkl
else if(ms1.eq.ms2)then
if(L1.gt.L2)then
kkL=kkL+1
kkla(l1,m1,ms1,l2,m2,ms2,1) = kkl
kkL=kkL+1
kkla(l1,m1,ms1,l2,m2,ms2,2) = kkl
else if(L1.eq.L2)then
if(m1.ge.m2)then
kkL=kkL+1
kkla(l1,m1,ms1,l2,m2,ms2,1) = kkl
kkL=kkL+1
kkla(l1,m1,ms1,l2,m2,ms2,2) = kkl
endif
endif
endif
enddo
enddo
enddo
enddo
enddo
enddo
556 format(2i5,' real',6i3,' L1,M1,Ms1,L2,M2,Ms2')
557 format(2i5,' imag',6i3,' L1,M1,Ms1,L2,M2,Ms2')
do ms1=1,iso
ks1=1
if(ms1.eq.2)ks1=-1
do ms2=1,iso
ks2=1
if(ms2.eq.2)ks2=-1
do Lcase1=1,nl(icase)
L1=LL(icase,Lcase1)
do m1=-L1,L1
do Lcase2=1,nL(icase)
L2=LL(icase,Lcase2)
do m2=-L2,L2
if(ms1.gt.ms2)then
write(70+icase,556)ic1, kkla(L1,m1,ms1,L2,m2,ms2,1), L1,m1,ks1,L2,m2,ks2
write(70+icase,557)ic1, kkla(L1,m1,ms1,L2,m2,ms2,2), L1,m1,ks1,L2,m2,ks2
else if(ms1.eq.ms2)then
if(L1.gt.L2)then
write(70+icase,556)ic1, kkla(L1,m1,ms1,L2,m2,ms2,1), L1,m1,ks1,L2,m2,ks2
write(70+icase,557)ic1, kkla(L1,m1,ms1,L2,m2,ms2,2), L1,m1,ks1,L2,m2,ks2
else if(L1.eq.L2)then
if(m1.ge.m2)then
write(70+icase,556)ic1, kkla(L1,m1,ms1,L2,m2,ms2,1), L1,m1,ks1,L2,m2,ks2
write(70+icase,557)ic1, kkla(L1,m1,ms1,L2,m2,ms2,2), L1,m1,ks1,L2,m2,ks2
endif
endif
endif
enddo
enddo
enddo
enddo
enddo
enddo ! end writing popmat input for tetra
endif
!.....WRITE QTL TO FILE 31
do 10 jb0=1,nband
rewind(29)
WRITE(31,2040) jb0
DO 11 MK=1,nkpt
READ(29) R,S,T,NV,NE,KNAME
if(popmat.eq.0)then
sumL=0.
do lcase=1,nl(icase)
sumL1(lcase)=0.
DO 20 JB=1,nb(mk)
READ(29) NUM,Ee
if(JB.eq.JB0)Eb=Ee
DO JL=1,iL(lcase)
READ(29)QTL
if(so .and. iso2.eq.1) qtl=qtl/2.d0
if(JB.eq.JB0)then
qtl1(lcase,jl)=qtl
sumL1(lcase)=sumL1(lcase)+qtl
sumL=sumL+qtl ! sum of all DOS's considered
endif
END DO
20 CONTINUE
enddo !end Lcase loop
kL=0 ! reorganize sumL1,qtL1 -> qtl2
do lcase=1,nL(icase)
kL=kL+1
qtl2(kL)=suml1(lcase)
L1=LL(icase,Lcase)
if(L1.gt.0)then
do 111 jL=1,iL(lcase)
if((L1.eq.1).and.(nm(icase).eq.5))goto 111 !do not write x+y+z in cubic symmetry
kL=kL+1
qtl2(kL)=qtl1(Lcase,jL)
111 continue
endif
enddo
!pb WRITE(31,1050) Eb,jatom,kL,(qtl2(jL),jL=1,kL),sumL
WRITE(31,1050) Eb,icase,kL,suml,(qtl2(jL),jL=1,kL)
kLmax=max(kLmax,kL) !pb
else ! popmat starts
! JR: we initialize sums here!
suml=0.
do lcase1=1,nl(icase)
suml1(lcase1)=0.
enddo
! ----------------------------
DO 21 JB=1,nb(mk)
READ(29) NUM,Ee
if(JB.eq.JB0)Eb=Ee
kl=0
! JR: suml=0.
if((popmat.eq.98).or.(popmat.eq.99))then ! complete matrix with <L1||L2> terms
do Lcase1=1,nl(icase)
! JR: suml1(lcase1)=0.
L1=LL(icase,Lcase1)
do ms1=1,iso
do m1=-L1,L1
do Lcase2=1,nL(icase)
L2=LL(icase,Lcase2)
do ms2=1,iso
do m2=-L2,L2
if(ms1.gt.ms2)then
read(29)qtle
! JR: this condition solves the fact, that telnes2 does not expect spin-crossterms in SO calculation
if((popmat.eq.98).or.((ms1.eq.1).and.(ms2.eq.1))) then
kL=kL+1
if(jb.eq.jb0) qtle1(kL)=qtle
endif
eLse if(ms1.eq.ms2)then
if(L1.gt.L2)then
read(29)qtle
! JR: this condition solves the fact, that telnes2 does not expect spin-crossterms in SO calculation
if((popmat.eq.98).or.((ms1.eq.1).and.(ms2.eq.1))) then
kL=kL+1
if(jb.eq.jb0) qtle1(kL)=qtle
endif
eLse if(L1.eq.L2)then
if(m1.ge.m2)then
read(29)qtle
! JR: this condition solves the fact, that telnes2 does not expect spin-crossterms in SO calculation
if((popmat.eq.98).or.((ms1.eq.1).and.(ms2.eq.1))) then
kL=kL+1
if(jb.eq.jb0) qtle1(kL)=qtle
if(m1.eq.m2)then
if(jb.eq.jb0) suml1(lcase1)=suml1(lcase1)+qtle
if(jb.eq.jb0) suml=suml+qtle
endif
endif
endif
endif
endif
enddo
enddo
enddo
enddo
enddo
enddo
else if((popmat.eq.87).or.(popmat.eq.88))then ! matrix with <L1||L1> diagonal terms only
do Lcase1=1,nl(icase)
! JR: suml1(lcase1)=0.
L1=LL(icase,Lcase1)
do ms1=1,iso
do m1=-L1,L1
do Lcase2=1,nL(icase)
L2=LL(icase,Lcase2)
do ms2=1,iso
do m2=-L2,L2
if(ms1.gt.ms2)then
read(29)qtle
if((L1.eq.L2).and.(m1.eq.m2))then
! JR: this condition solves the fact, that telnes2 does not expect spin-crossterms in SO calculation
if(popmat.eq.87) then
kL=kL+1
if(jb.eq.jb0) qtle1(kL)=qtle
endif
endif
else if(ms1.eq.ms2)then
if(L1.gt.L2)then
read(29)qtle
else if(L1.eq.L2)then
if(m1.ge.m2)then
read(29)qtle
if(m1.eq.m2)then
! JR: this condition solves the fact, that telnes2 does not expect spin-crossterms in SO calculation
if((popmat.eq.87).or.(ms1.eq.1)) then
kL=kL+1
if(jb.eq.jb0)then
qtle1(kL)=qtle
suml1(lcase1)=suml1(lcase1)+qtle
suml=suml+qtle
! write(*,'(i3,2f10.6,3x,5f10.6)') kl, qtle, suml, suml1(1), suml1(2), suml1(3), suml1(4)
endif
endif
endif
endif
endif
endif
enddo
enddo
enddo
enddo
enddo
enddo
endif
21 continue
! WRITE(31,1050) Eb,jatom,kL,1.,sumL,(suml1(jL),jL=1,4)
! JR: WRITE(31,1050) Eb,icase,kL,1.,sumL,(suml1(jL),jL=1,4)
WRITE(31,1050) Eb,icase,kL,sumL,(suml1(jL),jL=1,nl(icase))
if(popmat.eq.88)then
write(31,208)(dble(qtle1(jL)),jL=1,kL)
else if(popmat.eq.87)then
write(31,209)(dble(qtle1(jL)),jL=1,kL)
else if(popmat.eq.99)then
write(31,208)(qtle1(jL),jL=1,kL)
else if(popmat.eq.98)then
write(31,209)(qtle1(jL),jL=1,kL)
endif
209 format(10E16.8)
208 FORMAT(10F10.5)
endif ! end popmat vs qtl
sumL=0.
! WRITE(31,1050) Eb,jatom+1,1,sumL ! because of interstitial
WRITE(31,1050) Eb,icase+1,1,sumL ! because of interstitial
11 CONTINUE
10 CONTINUE
RETURN
!
! Error messages
!
CALL OUTERR('OUTP',ERRMSG)
STOP 'OUTP - Error'
!
1000 FORMAT(/,10X,3F8.4,1X,I5,I4,2X,A10)
1001 FORMAT(/,1X,' K-POINT:',3F8.4,1X,I5,I4,2X,A10)
1009 FORMAT(8X,I3,4X,F9.5,9X,F10.7,4x,i8)
1010 FORMAT(8X,I3,4X,F9.5,9X,F10.7,12x,5f8.5)
1011 FORMAT(1X,' BAND #',I3,' E=',F9.5,' WEIGHT=',F10.7,' NL=',i8,5f8.5)
1020 FORMAT(2e17.8)
1050 FORMAT(F10.5,I3,I5,F8.5,3X,39F8.5)
2040 FORMAT(1X,'BAND',I5)
END
SUBROUTINE L2MAIN(cmplx,nsymop,popmat,qtlfn,so)
USE param
USE struct
USE case
USE sym2
USE com
USE abc
IMPLICIT REAL*8 (A-H,O-Z)
CHARACTER *10 BNAME
CHARACTER *180 VECFN,FNAME,qtlfn
LOGICAL notcalc,cmplx , so
REAL*8,ALLOCATABLE :: a_real(:),fj(:,:),dfj(:,:)
REAL*8,ALLOCATABLE :: weight(:)
DIMENSION ilo(0:lomax)
integer popmat
COMPLEX*16 YL((LMAX2+1)*(LMAX2+1)),densmat(7,7)
COMPLEX*16 PHSHEL,CFAC,IMAG,CZERO,PH_SPIN(2)
COMMON /GENER/ BR1(3,3),BR2(3,3)
COMMON /ATSPDT/ P(0:LMAX2,2,nrf),DP(0:LMAX2,2,nrf)
COMMON /RADFU/ RF1(NRAD,0:LMAX2,2,nrf),RF2(NRAD,0:LMAX2,2,nrf)
COMMON /RINTEG/ RI_MAT(0:lmax2,0:lmax2,nrf,nrf,2,2)
COMMON /XA/ R(NRAD),BK(3),BKROT(3),BKRLOC(3)
logical loor(0:lomax),lapw(0:lmax2)
common /loabc/ alo(0:lomax,2,nloat,nrf)
common /lolog/ nlo,nlov,nlon,loor,ilo,lapw
COMMON /PROC/ VECFN(2)
COMMON /IPROC/ IPROC
!...............................
complex*16 h_yl(2*LMAX2+1,iblock)
complex*16 h_alyl(2*LMAX2+1,iblock,2)
complex*16 h_blyl(2*LMAX2+1,iblock,2)
!...............................
DATA CZERO/(0.0D0,0.0D0)/,IMAG/(0.0D0,1.0D0)/
!------------------------------------------------------------------
assign 2022 to iform1
2022 FORMAT(3X,4E19.12)
PI=ACOS(-1.0D0)
TWOPI=2.D0*PI
EDELT=0.0D0
TEST=0.0D0
ETOT=0.0D0
MAXWAV=0
MINWAV=100000
DO JK=1,NKPT
NB(JK)=0
END DO
do i=1,iord
do ii=1,3
do jj=1,3
tmat(ii,jj,i)=iz(ii,jj,i)
end do
end do
end do
ALLOCATE(a_real(nmat),fj(0:lmax2,nmat),dfj(0:lmax2,nmat))
allocate (weight(nume))
if(so .and. iso2.eq.1) print*, 'non-spinpolarized so calculation'
! ---------------------------------
! START LOOP FOR ALL ATOMS BY DO 50
! ---------------------------------
do is=1,iso2
READ(17+is,2032) ISCF
end do
LFIRST=1
icount=1
DO 50 JATOM=1,NAT
if(icount.gt.natom)return
if(SO)then
iso=2 ! s-o calculation, but nonrelativistic basis
if(nm(jatom).lt.80)then ! exclude density matrix calculation
if(nm(jatom).ne.-1)then ! exclude relativistic basis
if(nm(jatom).ne.0)then ! exclude relativistic basis
if(nm(jatom).ne.6)then ! exclude user defined basis
iso=1
endif
endif
endif
endif
endif
notcalc=.true.
if (jatom.eq.iatom(icount)) then
notcalc=.false.
icase=icount
icount=icount+1
end if
if(notcalc)then
write(6,*)' ATOM=',JATOM,' left out '
else
write(6,*)' Calculation for atom',jatom
if(nm(icase).gt.50)then
do lcase1=1,nl(icase)
L1=LL(icase,lcase1)
n1=(2*L1+1)*ISO
do lcase2=1,nl(icase)
L2=LL(icase,lcase2)
n2=(2*L2+1)*iso
do i1=1,n1
do i2=1,n2
dmat(L1,L2,i1,i2)=(0.d0,0.d0)
enddo
enddo
enddo
enddo
endif
endif
IF(JATOM.GT.1) LFIRST=LFIRST + MULT(JATOM-1)
ITAP=30+JATOM
! CALCULATE RADIAL FUNCTIONS U(R), UE(R), ...
do is=1,iso
itape=8+is
jtape=17+is
rewind(itape)
if(is.eq.2.and.iso2.eq.1) then ! fix for SO-coupling and non-spinpolarized
jtape=18
! rewind(jtape)
! READ(jtape,2032) ISCF
endif
CALL ATPAR(JATOM,LFIRST,itape,jtape,is,iso,popmat)
rewind(itape)
end do
if (notcalc) goto 50
FAC=4.0D0*PI*RMT(JATOM)**2/SQRT(VOL)
!................................
KKK=0
!.....READ IN K POINT AND BASIS VECTORS
iloop=0
788 continue
if (IPROC.GT.0) then
iloop=iloop+1
write(6,*)'ILOOP ',iloop
do is=1,iso
itape=8+is
close(itape)
call mknam(FNAME,VECFN(is),ILOOP)
write(6,*)'opening ',FNAME
open(itape,FILE=FNAME,STATUS='old',FORM='unformatted',ERR=950)
end do
endif
DO I=1,NAT
do is=1,iso
itape=8+is
READ(itape) EMIST
READ(itape) EMIST
enddo
ENDDO
4 continue
! call inispl this took too much CPU time
do lcase1=1,nl(icase)
L1=LL(icase,lcase1)
n1=(2*L1+1)*ISO
do lcase2=1,nl(icase)
L2=LL(icase,lcase2)
n2=(2*L2+1)*iso
do i1=1,n1
do i2=1,n2
do num=1,nume
xqtl(L1,L2,i1,i2,num)=(0.,0.)
enddo
enddo
enddo
enddo
enddo
!......reading from vector for both spins
DO 17 is=1,iso
itape=8+is
READ(itape,END=998) S,T,Z,BNAME,N,NE
IF(N.GT.MAXWAV) MAXWAV=N
IF(N.LT.MINWAV) MINWAV=N
if (is.eq.1) then
KKK=KKK+1
WRITE(29) S,T,Z,N,NE,BNAME
end if
205 FORMAT(/,1X,' K-POINT:',3F8.4,1X,I5,I4,2X,A10)
! IF(icase.EQ.1.and.is.eq.1) NK=NK+1
READ(itape) (KX(I),KY(I),KZ(I),I=1,N)
DO 5 I=1,N
BKX(I)=(S+KX(I))
BKY(I)=(T+KY(I))
BKZ(I)=(Z+KZ(I))
5 CONTINUE
CALL HARMON(N,BKX,BKY,BKZ,LMAX2,FJ,DFJ,RMT(JATOM))
NEMIN=1
NEMAX=0
14 READ(itape) NUM,exxx
E(NUM)=exxx
if (.not.cmplx) then
READ(itape) (A_real(I),I=1,N)
do i=1,n
a(i,num,is)=(1.d0,0.d0)*a_real(i)
end do
else
READ(itape) (A(I,NUM,is),I=1,N)
end if
weight(num)=weigh(kkk,num)
IF(E(NUM).LT.EMIN) NEMIN=NEMIN+1
IF(E(NUM).LT.EMAX) NEMAX=NEMAX+1
IF(NUM.EQ.NE) GOTO 16
GOTO 14
16 CONTINUE
NB(KKK)=NEMAX-NEMIN+1
17 CONTINUE
!.......reading from vect end
!.......sum over sym. operations
DO 777 MU=1,nsymop
! loop over L(jatom)
do lcase=1,nl(icase)
l=ll(icase,lcase)
CFAC=IMAG**L
DO IS=1,iso
PH_SPIN(is)=EXP((2*is-3)*IMAG*PHASE(MU)/2)
DO 9 I=1,2*LMAX2+1
DO 9 NUM=NEMIN,NEMAX
DO 9 irf=1,nrf
ALM(I,num,irf,is)=CZERO
9 CONTINUE
END DO
DO 120 iI=1,N-(nlo+nlon+nlov),iblock
i3=0
do 121 i=ii,min(ii+iblock-1,N-(nlo+nlon+nlov))
i3=i3+1
BK(1)=BKX(I)
BK(2)=BKY(I)
BK(3)=BKZ(I)
CALL ROTATE (BK,TMAT(1,1,mu),BKROT)
BK(1)=BKROT(1)*BR1(1,1)+BKROT(2)*BR1(1,2)+BKROT(3)*BR1(1,3)
BK(2)=BKROT(1)*BR1(2,1)+BKROT(2)*BR1(2,2)+BKROT(3)*BR1(2,3)
BK(3)=BKROT(1)*BR1(3,1)+BKROT(2)*BR1(3,2)+BKROT(3)*BR1(3,3)
CALL ROTATE (BK,ROTLOC(1,1,JATOM),BKRLOC)
CALL YLM (BKRLOC,LMAX2,YL)
ARG1=BKROT(1)*POS(1,LFIRST)*TWOPI
ARG2=BKROT(2)*POS(2,LFIRST)*TWOPI
ARG3=BKROT(3)*POS(3,LFIRST)*TWOPI
ARGT=(BKX(I)*TAU(1,MU)+BKY(I)*TAU(2,MU)+ &
BKZ(I)*TAU(3,MU))*TWOPI
PHSHEL=EXP(IMAG*(ARG1+ARG2+ARG3+ARGT))
DO M=1,2*L+1
IND_YL=M+L*L
h_yl(M,i3)=conjg(yl(ind_yl))*phshel
END DO
121 continue
DO 122 IS=1,iso
i3=0
do i=ii,min(ii+iblock-1,N-(nlo+nlon+nlov))
i3=i3+1
DO M=1,2*L+1
if (lapw(l)) then
h_ALYL(m,i3,is)=(DFJ(L,I)*P(l,is,2)-FJ(L,I)*DP(l,is,2))* &
h_yl(M,i3)*ph_spin(is)
h_BLYL(m,i3,is)=(FJ(L,I)*DP(l,is,1)-DFJ(L,I)*P(l,is,1))* &
h_yl(M,i3)*ph_spin(is)
else
h_ALYL(m,i3,is)=FJ(L,I)/P(l,is,1)/RMT(JATOM)**2* &
h_yl(M,i3)*ph_spin(is)
h_BLYL(m,i3,is)=(0.d0,0.d0)
end if
END DO
enddo
ibb=min(iblock,N-(nlo+nlon+nlov)-ii+1)
lda=2*LMAX2+1
ldc=lda
ldb=nmat
call zgemm('N','N',2*l+1,nemax-nemin+1,ibb,(1.d0,0.d0), &
h_alyl(1,1,is),lda,a(ii,nemin,is),ldb,(1.d0,0.d0), &
alm(1,nemin,1,is),ldc)
call zgemm('N','N',2*l+1,nemax-nemin+1,ibb,(1.d0,0.d0), &
h_blyl(1,1,is),lda,a(ii,nemin,is),ldb,(1.d0,0.d0), &
alm(1,nemin,2,is),ldc)
122 CONTINUE
120 CONTINUE
if (nlo.ne.0) then
call &
lomain (nemin,nemax,lfirst,lfirst,n,jatom, &
mu,L,iso)
end if
DO 230 M=1,2*L+1
DO 892 NUM=NEMIN,NEMAX
DO 892 is=1, iso
DO 892 irf=1,nrf
ALML(L,M,NUM,irf,is)=ALM(M,NUM,irf,is)*FAC*CFAC
892 CONTINUE
230 CONTINUE
enddo !end of atom L loop (alm)
do lcase1=1,nl(icase)
l1=ll(icase,lcase1)
if(popmat.eq.0)then
CALL XSPLT0(nemin,nemax,l1,l1,MU,iso,nsymop)
else
do lcase2=1,lcase1
l2=ll(icase,lcase2)
CALL XSPLT(nemin,nemax,l1,l2,MU,iso,weight,nsymop)
enddo
endif
enddo
777 CONTINUE !end of symmetry loop
if(popmat.eq.0)then
do lcase=1,nl(icase)
CALL QTL(nemin,nemax,icase,lcase,iso)
enddo
endif
CALL PSPLIT(nemin,nemax,icase,iso,popmat)
GOTO 4 !k-points end
998 CONTINUE !k-points end (jump)
if (iloop.lt.iproc) goto 788
close(31)
call mknam(FNAME,qtlfn,Icase)
open(31,FILE=FNAME,STATUS='unknown',FORM='formatted',ERR=951)
CALL OUTP(ICASE,popmat,so)
if(nm(icase).gt.50)then ! population matrix output
write(12,200)jatom
write(21,*)
write(21,200)jatom
write(6,*)
write(6,200)jatom
200 format(i5,' atom population matrix')
do lcase1=1,nl(icase)
l1=ll(icase,lcase1)
if(popmat.gt.90)then
low=1
else
low=lcase1
endif
do lcase2=low,lcase1
l2=ll(icase,lcase2)
N1=2*L1+1
N2=2*L2+1
if(popmat.gt.90)then
write(6,201)l1,l2
write(21,201)l1,l2
else
write(6,202)l1
write(21,202)l1
endif
201 format(' L1=',i2,', L2=',i2)
202 format(' L=',i2)
N1=2*L1+1
N2=2*L2+1
do is1=1,iso
do is2=is1,iso
do m1=-l1,l1
ly=L1+1+m1
ind1=L1+1+m1+N1*(is1-1)
do m2=-l2,l2
lpy=l2+1+m2
ind2=L2+1+m2+N2*(is2-1)
write(12,655)dmat(l1,l2,ind1,ind2),is1,is2,l1,l2,m1,m2
densmat(ly,lpy)=dmat(l1,l2,ind1,ind2)
654 format(6i3,2f14.8)
655 format(2f14.8,6i3,' Ms1,Ms2,L1,L2,ML1,ML2')
enddo
enddo
call OUTDMAT(L1,L2,N1,N2,IS1,IS2,ISO,densmat)
enddo
enddo
enddo
enddo
endif
! CLOSE(15)
rewind(29)
50 CONTINUE
! ....END LOOP OVER ALL ATOMS
REWIND(itape)
!
!********************************************************************PH*
!.... WRITE Output for the BR1 matrix and the ROTIJ Matrix in file case.rotlm:
! they are read by LectureROTIJ in lecture.f in SRC_elnes / or SRC_mdff/.
! These lines are in SRC_lapw2lm/l2main.frc: surch *PH*
! The number of the output file is 22. his extention must be .rotlm
794 CALL ROTDEF
DO JR=1, 3
WRITE(22, 101) BR1(1,JR), BR1(2,JR), BR1(3,JR)
ENDDO
INDEX = 0
DO LATOM=1, NAT
DO LATEQ=1, MULT(LATOM)
INDEX = INDEX + 1
WRITE(22, 100) LATOM, LATEQ, INDEX
DO JR=1, 3
WRITE(22, 101) (ROTIJ(JC,JR,INDEX),JC=1,3)
ENDDO
ENDDO
ENDDO
100 FORMAT('inequivalent atomnumber ',I3,' number ',I2,' total ',I4)
101 FORMAT(3F10.5)
!**********************************************************************PH*
RETURN
950 CALL OUTERR('l2main','error reading parallel vectors')
STOP 'L2main - Error'
951 CALL OUTERR('l2main','error opening qtl_XX files')
STOP 'L2main - Error'
!
2032 FORMAT(50X,I2,//)
END
_______________________________________________
Wien mailing list
Wien@zeus.theochem.tuwien.ac.at
http://zeus.theochem.tuwien.ac.at/mailman/listinfo/wien
SEARCH the MAILING-LIST at:
http://www.mail-archive.com/wien@zeus.theochem.tuwien.ac.at/index.html
