No subject
Fri Nov 19 11:58:34 CET 2010
----< local rotation matrix >-----
New local rotation matrix in global orthogonal system
new x new y new z
LOCAL ROT MATRIX: -0.7643046 0.6448554 0.0000000
0.6042329 0.7161574 0.3493153
0.2252578 0.2669833-0.9370053
---< another choice of local rotation matrix >---
New local rotation matrix in global orthogonal system
new x new y new z
LOCAL ROT MATRIX: 0.7643046-0.6448554 0.0000000
-0.6042329-0.7161574 0.3493153
-0.2252578-0.2669833-0.9370053
(Compared with previous one, newx --> -newx, newy --> -newy)
These two chocies should give same decomposition of d-states into five
different orbitals. But, qtl gives completely different decomposition
for these two choices.
The problem only arises when local coordinates are specified in case.inq.
If I modify local rot matrix in case.struct and execute qtl program,
the problem does not appear.
In my understanding, the problem arises from the definition
of local rotation matrix in QTL code.
When genrating local rot matrix in QTL, column and row are interchanged compared with
lapw2 and other program. When I modified QTL code
and used same defnition as in lapw2, I got expected results.
Attached file is modified qtlmain.f. (v10.1, tested only for non spin-orbit case)
Best regards,
--
National Institute for Materials Science (NIMS)
Computational Materials Science Center (CMSC)
Masao Arai
--Boundary_(ID_mCpo2DQOnEfdFVv6S7Rfrw)
Content-type: text/x-fortran; charset=utf-8; name=qtlmain.f
Content-transfer-encoding: 7BIT
Content-disposition: attachment; filename=qtlmain.f
PROGRAM QTLMAIN
! Program QTL calculates atomic-site projected densities of states.
! For description of program see 'QTL-Technical Report'
USE param
USE struct
USE case
USE sym2
USE com
USE abc
IMPLICIT REAL*8 (A-H,O-Z)
complex*16 trmat(-7:7,-7:7),dimag,qtl1(24)
dimension rcf(14,14),xx(3),xz(3),acf(14),bcf(14)
CHARACTER*4 adum
CHARACTER*5 MODSYM,CHAR
CHARACTER*10 KNAME
CHARACTER*11 STATUS,FORM
CHARACTER*67 ERRMSG
CHARACTER*80 DEFFN, ERRFN
CHARACTER*80 FNAME,FNAME1,enefn,enefnd,qtlfn
CHARACTER*80 VECFN
CHARACTER*161 TEXT
REAL*8,ALLOCATABLE :: qtle(:,:),eb(:),suml(:)
real*8 qtle1(6144) ! JR: for merging files with popmat=88 or 99
INTEGER,ALLOCATABLE :: ieband(:),kL1(:)
INTEGER IPROC,popmat,qsplit,symmetrize
LOGICAL SO,cmplx
COMMON /GENER/ BR1(3,3),BR2(3,3)
COMMON /RADFU/ RRAD1(NRAD,0:LMAX2),RADE1(NRAD,0:LMAX2), &
RRAD2(NRAD,0:LMAX2),RADE2(NRAD,0:LMAX2)
COMMON /PROC/ VECFN(2)
COMMON /IPROC/ IPROC
DIMENSION S(3),saveiz(3,3)
DIMENSION xnew(3), ynew(3), znew(3)
DATA SO /.false./,cmplx /.false./
!-----------------------------------------------------------------------
dimag=(0.,1.)
CALL GTFNAM(DEFFN,ERRFN,IPROC)
CALL ERRFLG(ERRFN,'Error in QTL')
nspin1=1
! nspin1=2 changed 5.4.08 PN
OPEN (1,FILE=DEFFN,STATUS='OLD',ERR=910)
iso=1
iso2=1
10 CONTINUE
READ (1,*,END=20,ERR=960) IUNIT,FNAME,STATUS,FORM,IRECL
OPEN (IUNIT,FILE=FNAME,STATUS=STATUS,FORM=FORM,ERR=920)
if(iunit.eq.9) VECFN(1)=FNAME
if(iunit.eq.10) VECFN(2)=FNAME
if(iunit.eq.18) then
do i=80,5,-1
if(fname(i:i).ne.' ') then
if(fname(i-2:i).eq.'pup') nspin1=2
if(fname(i-2:i).eq.'pdn') nspin1=2 ! added 10.07.08 PN
! if(fname(i-2:i).eq.'vsp') nspin1=1 changed 5.4.08 PN
endif
enddo
endif
if(iunit==31) qtlfn=fname
if(iunit==59) enefnd=fname
if(iunit==60) enefn=fname
if(iunit.eq.7) then
read(iunit,123,end=12)adum
123 format(A5)
cmplx=.true.
12 continue
end if
if(iunit.eq.9) then
do i=80,1,-1
if(fname(i:i).ne.' ') then
if(fname(i-1:i).eq.'so') then
SO=.true.
iso=2
iso2=1
cmplx=.true.
endif
! JR: if(fname(i-3:i).eq.'soup') then
if((fname(i-3:i).eq.'soup').or.(fname(i-3:i).eq.'sodn')) then
SO=.true.
iso=2
iso2=2
cmplx=.true.
endif
goto 10
endif
enddo
endif
GOTO 10
20 CONTINUE
CLOSE (1)
!.....READ STRUCT
CALL init_struct
write(6,*)'******Program QTL: projected DOS or population matrix******'
815 format(' Compound:',a60)
if (cmplx) then
write(6,*)'CALCULATION WITH COMPLEX VECTORS'
endif
if (so)then
if(iso2.eq.2)write(6,*)'CALCULATION WITH SOC-UP AND DN VECTORS REQUIRED'
if(iso2.eq.1)write(6,*)'CALCULATION WITH SOC VECTORS REQUIRED'
endif
IF(IPROC.GT.0)THEN
write(6,*)'Running QTL on ',IPROC,' processors'
write(6,*)' '
do is=1,2
call mknam(FNAME,VECFN(is),1)
iunit=8+is
status='UNKNOWN'
form='UNFORMATTED'
OPEN (IUNIT,FILE=FNAME,STATUS=STATUS,FORM=FORM,ERR=920)
end do
ELSE
write(6,*)'Running QTL in single processor mode'
write(6,*)' '
END IF
! find Fermi energy using scf2 file
666 format(a4)
667 format(38x,f10.5)
13 read(8,666,end=14)adum
if(adum.eq.':NOE')then
read(8,*)
read(8,667)EF
goto 14
endif
goto 13
14 continue ! Fermi energy end
READ(5,*) EMIN,EMAX
read(5,*)natom
ALLOCATE(CF(NDIM2,NDIM2,natom,4),ISUM(0:ndim2,natom,4),nl(natom))
allocate(eb(nat),qtle(nat,34),suml(nat))
allocate(ll(natom,0:4),iatom(natom),nm(natom))
887 format(' Projected density of states calculation for',i3,' atoms')
write(6,887)natom
write(21,887)natom
write(6,877)Emin,Emax
write(21,877)Emin,Emax
877 format(' Energy window:',f10.4,' < E < ',f10.3)
CALL init_sym2(iord)
WRITE(6,800)
WRITE(6,805) TITLE
WRITE(6,810) LATTIC
WRITE(6,820) AA,BB,CC
WRITE(6,840) NAT
WRITE(6,850) IREL
write(30,4030)title
4030 format(' Ordering of DOS in QTL file for: ',a60,/)
CALL LATGEN
write(6,*)'IORD=',IORD
DO I=1,NATOM
READ(5,*)IATOM(I),qsplit,symmetrize,loro
read(5,*)nL(i),(LL(i,j),j=1,nL(i))
!Pavel fix for SO-vector, but nonrel. basis
if(SO)then ! s-o calculation, but nonrelativistic basis
if(qsplit.lt.80)then ! exclude density matrix calculation
if(qsplit.ne.-1)then ! exclude relativistic basis
if(qsplit.ne.0)then ! exclude relativistic basis
if(qsplit.ne.6)then ! exclude user defined basis
iso=1
endif
endif
endif
endif
! iso=1
! if((qsplit.eq.-1).or.(qsplit.eq.0))iso=2
endif
!pb qsplit=-2 uses the isplit parameter of the struct file
if(qsplit.eq.-2) then
if(isplit(iatom(i)).eq.8) then
qsplit=2
elseif(isplit(iatom(i)).eq.2) then
qsplit=5
elseif(isplit(iatom(i)).eq.4) then
qsplit=4
elseif(isplit(iatom(i)).eq.-2) then
qsplit=3
else
qsplit=2
endif
endif
nm(i)=qsplit
write(6,888)i,iatom(i),qsplit,(LL(i,j),j=1,nL(i))
write(21,888)i,iatom(i),qsplit,(LL(i,j),j=1,nL(i))
888 format(' atom',i3,'; type ',I3,'; qsplit=',i3,'; for L=',4i3)
if(symmetrize.ne.0)then
nsymop=iord
write(6,*)'Symmetrization over eq. k-points is performed'
write(21,*)'Symmetrization over eq. k-points is performed'
else
nsymop=1
write(21,*)'Symmetrization over eq. k-points is not performed'
write(6,*)'Symmetrization over eq. k-points is not performed'
write(6,*)'allowed for invariant DOS'
end if
ia=iatom(i)
if(loro.gt.0)then ! change of local rotation matrix
read(5,*)(xz(j),j=1,3) ! new z-axis expressed in unit cell vectors
write(6,120)(xz(j),j=1,3)
write(21,120)(xz(j),j=1,3)
120 format(' New z axis || ',3f9.4)
call angle(xz,thetaz,phiz)
znew(1) = sin(thetaz)*cos(phiz)
znew(2) = sin(thetaz)*sin(phiz)
znew(3) = cos(thetaz)
if(loro.eq.1)then ! only z-axis fixed
xnew(1) = 0. ! new x perpendicular to new z
ddd=abs(znew(3)**2-1.)
if(ddd.gt.0.000001)then
xnew(1) = -znew(2)
xnew(2) = znew(1)
else
xnew(1) = 1.
xnew(2) = 0.
endif
ddd=0. ! normalize new x
do j=1,3
ddd=ddd+xnew(j)**2
enddo
do j=1,3
xnew(j) = xnew(j)/sqrt(ddd)
enddo
elseif(loro.eq.2)then ! also new x-axis fixed
read(5,*)(xx(j),j=1,3)
write(6,121)(xx(j),j=1,3)
write(21,121)(xx(j),j=1,3)
121 format(' New x axis || ',3f9.4)
call angle(xx,thetax,phix)
xnew(1) = sin(thetax)*cos(phix)
xnew(2) = sin(thetax)*sin(phix)
xnew(3) = cos(thetax)
! check orthogonality of new x and z axes
124 check=0.
do j=1,3
check=check+xnew(j)*znew(j)
enddo
if(abs(check).gt.0.00001)then
write(6,'(" new x and z axes are not orthogonal",f8.5,f8.2," degrees")') &
abs(check),acos(check)*180./acos(-1.d0)
print*,'new x and z axes are not orthogonal for atom',ia
print*,'angle:',acos(check)*180./acos(-1.d0),"degrees"
write(6,'(" vector x",3f10.5)') xnew(1:3)
write(6,'(" vector z",3f10.5)') znew(1:3)
if(abs(check).gt.0.2d0) stop 'x,z vectors not orthogonal'
! x' = x + a (zx) with x'.z=0
check1=xnew(1) * znew(1)
check2=xnew(2) * znew(2)
check3=xnew(3) * znew(3)
xz1=xnew(1) - znew(1)
xz2=xnew(2) - znew(2)
xz3=xnew(3) - znew(3)
ang=-(check1+check2+check3)/(xz1*znew(1) + xz2*znew(2) + xz3*znew(3))
xnew(1) = xnew(1) + ang * xz1
xnew(2) = xnew(2) + ang * xz2
xnew(3) = xnew(3) + ang * xz3
ddd=0. ! normalize new x
do j=1,3
ddd=ddd+xnew(j)**2
enddo
do j=1,3
xnew(j) = xnew(j)/sqrt(ddd)
enddo
print*, 'small non-orthogonality corrected '
write(6,*) 'small non-orthogonality corrected to:'
write(6,'(" vector x:",3f10.5)') xnew(1:3)
goto 124
endif
endif
ynew(1) = znew(2)*xnew(3) - znew(3)*xnew(2)
ynew(2) = znew(3)*xnew(1) - znew(1)*xnew(3)
ynew(3) = znew(1)*xnew(2) - znew(2)*xnew(1)
rotloc(1,1:3,ia) = xnew(1:3)
rotloc(2,1:3,ia) = ynew(1:3)
rotloc(3,1:3,ia) = znew(1:3)
write(6,*)' New local rotation matrix in global orthogonal system'
write(6,*)' new x new y new z'
write(21,*)' New local rotation matrix in global orthogonal system'
write(21,*)' new x new y new z'
write(6,1013)((rotloc(j1,j,ia),j1=1,3),j=1,3) !written as in .struct
write(21,1013)((rotloc(j1,j,ia),j1=1,3),j=1,3)
! write(6,1013)((rotloc(j,j1,ia),j1=1,3),j=1,3) !written as in .struct
! write(21,1013)((rotloc(j,j1,ia),j1=1,3),j=1,3)
1013 format('LOCAL ROT MATRIX: ',3f10.7,/,20x,3f10.7,/,20x,3f10.7)
endif !change of loc.rot. end
if(qsplit.lt.6.and.qsplit.gt.0)then ! unitary trafo for DOS calculation
do j=1,nL(i) ! determined by program
L=LL(i,j)
call transf(L,qsplit,trmat) ! determine unitary transformation
if(iso.eq.1)then ! no s-o coupling
do m1=-L,L
ma=m1+L+1
do m2=-L,L
mb=m2+L+1
cf(ma,mb,i,j)=trmat(m1,m2)
enddo
enddo
else ! s-o coupling starts
nlm=2*L+1
! if(qsplit.gt.0)then ! trafo matrix for ms=1/2 = matrix (ms=-1/2)
do m1=-L,L
do m2=-L,L
cf(m1+nlm,m2+nlm,i,j)=trmat(m1,m2)
cf(m1+2*nlm,m2+2*nlm,i,j)=trmat(m1,m2)
enddo
enddo
endif
enddo
else if(qsplit.lt.1)then ! relativistic basis
do j=1,nL(i) ! determined by program
L=LL(i,j)
nlm=2*L+1
call trafoso(L,rcf)
do ma=1,2*nlm
do mb=1,2*nlm
cf(ma,mb,i,j)=rcf(ma,mb)
enddo
enddo
enddo
! endif
! endif
elseif(qsplit.eq.6)then ! unitary transformation read from files 31+ia
if(31+i.gt.58) stop 'only up to 27 atoms allowed for qsplit=6'
do j=1,nL(i)
L=LL(i,j)
write(6,101)L
101 format(' L=',i2,'. Unitary transformation read from input')
nlm=2*L+1
do ma=1,iso*nlm
read(31+i,*)(acf(mb),bcf(mb),mb=1,iso*nlm)
do mb=1,iso*nlm
cf(ma,mb,i,j)=acf(mb)+dimag*bcf(mb)
enddo
enddo
write(6,*)' Real part of unitary matrix'
do m1=1,iso*nlm
write(6,555)(dble(cf(m1,m2,i,j)),m2=1,iso*nlm)
enddo
write(6,*)' Imaginary part of unitary matrix'
do m1=1,iso*nlm
write(6,555)(imag(cf(m1,m2,i,j)),m2=1,iso*nlm)
enddo
enddo
555 format(7f9.4)
endif ! Unitary trafo for DOS end
popmat=0
if(qsplit.gt.80)then ! Population matrix start
popmat=qsplit
if((qsplit.eq.88).or.(qsplit.eq.99))then
write(6,*)' Population matrix for TELNES'
write(21,*)' Population matrix for TELNES'
endif
nbl=0
do j=1,nL(i)
nbl=nbl+iso*(2*LL(i,j)+1)
enddo
nbl=nbl*(nbl+1)+2
if((qsplit.eq.87).or.(qsplit.eq.88))then
write(6,111)(LL(i,j),j=1,nL(i))
write(21,111)(LL(i,j),j=1,nL(i))
111 format(' Population matrix diagonal in L for L=',4i3)
elseif((qsplit.eq.98).or.(qsplit.eq.99))then
if(70+i.gt.98) stop 'only up to 27 atoms supported for qsplit=98/99'
write(6,110)(LL(i,j),j=1,nL(i))
write(21,110)(LL(i,j),j=1,nL(i))
110 format(' Complete population matrix with <L1|L2> terms for L=',4i3)
112 format(' Number of cases for TETRA=',i4)
write(6,112)nbl
! begin to write input files for tetra (to be completed in outp)
write(70+i,*)' Population matrix data for tetra from QTL'
dE=0.0001
gauss=0.
write(70+i,776)Emin,dE,Emax,gauss
write(70+i,777)nbl
776 format(4f12.6,' Emin,dE,Emax,dgauss')
777 format(i5,' number of cases for tetra')
endif
endif
END DO ! Loop over atoms end
!....reading *.inso
if (so) then
do i=1,3
read(4,*)
end do
read(4,*)s(1),s(2),s(3)
126 format(' Magnetic system with s-o coupling; M theta, phi:',2f8.4)
call angle(s,theta,fi)
write(6,126)theta,fi
write(21,126)theta,fi
!.... test of collinearity of the z-axis with spin quantization axis
cost=cos(theta)
cosf=cos(fi)
sint=sin(theta)
sinf=sin(fi)
xz(1)=sint*cosf
xz(2)=sint*sinf
xz(3)=cost
do i=1,natom
ia=iatom(i)
sum=rotloc(1,3,ia)*xz(1)+rotloc(2,3,ia)*xz(2)+rotloc(3,3,ia)*xz(3)
! if(abs(sum-1.).gt.0.00001.and.(nm(i).eq.0.or.nm(i).eq.-1))then ! old z not along M, new rotloc
! if(abs(sum-1.).gt.0.00001)then ! old z not along M, new rotloc
if(nm(i).eq.0.or.nm(i).eq.-1)then ! old z not along M, new rotloc
rotloc(1,1,i)=cosf*cost
rotloc(1,2,i)=-sinf
rotloc(1,3,i)=cosf*sint
rotloc(2,1,i)=sinf*cost
rotloc(2,2,i)=cosf
rotloc(2,3,i)=sinf*sint
rotloc(3,1,i)=-sint
rotloc(3,2,i)=0
rotloc(3,3,i)=cost
write(6,*)' New local rotation matrix with z || M'
write(21,*)' New local rotation matrix with z || M'
write(6,1013)((rotloc(j1,j,i),j1=1,3),j=1,3) ! written as in .struct
write(21,1013)((rotloc(j1,j,i),j1=1,3),j=1,3)
! write(6,1013)((rotloc(j,j1,i),j1=1,3),j=1,3) ! written as in .struct
! write(21,1013)((rotloc(j,j1,i),j1=1,3),j=1,3)
endif
end do
end if
IF (nsymop.gt.1) THEN
if (so) then
CALL SYM(THETA,FI)
nsymop=iord
endif
ELSE ! put identity as first operation and exchange
saveiz(1:3,1:3)=iz(1:3,1:3,1)
do ind=1,iord
if(iz(1,1,ind).eq.1.d0.and.iz(2,2,ind).eq.1d0.and.iz(3,3,ind).eq.1.d0.and. &
iz(2,1,ind).eq.0.d0.and.iz(3,1,ind).eq.0.d0.and.iz(1,2,ind).eq.0.d0.and. &
iz(3,2,ind).eq.0.d0.and.iz(1,3,ind).eq.0.d0.and.iz(2,3,ind).eq.0.d0) then
do i=1,3
do j=1,3
iz(i,j,1)=0
if (i.eq.j) iz(i,j,1)=1
iz(i,j,ind)=saveiz(i,j)
end do
end do
! print*,'identity found as',ind,' operation'
endif
enddo
END IF
! find nkpt, nmat and nume in energy file
k=0
iloop=0
778 continue
if(IPROC.GT.0) then
iloop=iloop+1
close(59)
close(60)
call mknam(FNAME,enefn,ILOOP)
open(60,FILE=FNAME,STATUS='old',ERR=951)
call mknam(FNAME,enefnd,ILOOP)
open(59,FILE=FNAME,STATUS='unknown',ERR=951)
endif
jspin=1 !pb
DO I=1,NAT
READ(60,'(f9.5)') EMIST
READ(60,'(f9.5)') EMIST
ENDDO
DO I=1,NAT
READ(59,'(f9.5)',END=1005) EMIST
READ(59,'(f9.5)',END=1005) EMIST
ENDDO
JSPIN=2
1005 CONTINUE
DO
READ(60,'(3e19.12,a10,2i6)',IOSTAT=ios) SS,T,Z,KNAME,N,NEn
IF (ios /= 0) EXIT
k=k+1
nmat=MAX(n,nmat)
nume=MAX(nen,nume)
DO ii=1,nen
READ(60,*) NUM,E1
ENDDO
IF(jspin.EQ.2) THEN
READ(59,'(3e19.12,a10,2i6)',IOSTAT=ios) SS,T,Z,KNAME,N,NEn
nmat=MAX(n,nmat)
nume=MAX(nen,nume)
DO ii=1,nen
READ(59,*) NUM,E1
ENDDO
ENDIF
ENDDO
IF(ILOOP.LT.IPROC) goto 778
nkpt=k
REWIND(59)
REWIND(60)
CALL init_com(nume,nkpt)
CALL init_abc(nume,nmat,lmax2,ndim2,nrf)
if(nm(1).gt.50)then ! for pop. mat. calc.read weights calculated by lapw2
call readw ! which are used to integrate dmat over energy
write(6,*)' Weights for energy integration read'
endif
itape=10
jtape=18
call l2main(cmplx,nsymop,popmat,qtlfn,so)
WRITE(16,3000) TITLE(1:79)
WRITE(16,3010) AA,BB,CC,ef
if (.not. so) then
WRITE(16,3020) MINWAV,MAXWAV,nspin1,nat,iso-1,kLmax !pb
else
WRITE(16,3020) MINWAV,MAXWAV,nspin1,nat,2,kLmax !pb
endif
icase=1
do jatom=1,nat
if(icase.le.natom)then
ia=iatom(icase)
else
ia=-1
endif
if(ia.eq.jatom)then ! atom was selected
if(nm(icase).eq.-1)then
text=' projected DOS in relativistic |J,MJ> basis'
else if(nm(icase).eq.0)then
text=' |J,MJ> basis, L-1/2; L+1/2 sums'
else if(nm(icase).eq.1)then
text=' projected DOS in |L,M> basis'
else if(nm(icase).eq.2)then
text=' projected DOS in real basis'
else if(nm(icase).eq.3)then
text=' projected DOS in axial symmetry, [001] symmetry axis'
else if(nm(icase).eq.5)then
text=' projected DOS in cubic symmetry'
else if(nm(icase).eq.6)then
text=' projected DOS, user own unitary transformation'
else if(nm(icase).eq.88)then
text='tot,0,1,2,3,xdos(i,i),i=1,lxdos2)'
else if(nm(icase).eq.99)then
text='tot,0,1,2,3,xdos(i,j),j=1,i),i=1,lxdos2)'
! JR: headers for high-precision pop-mat calculation
else if(nm(icase).eq.87)then
write(text,'("NL=",i3,":",4i3)') nl(icase), (ll(icase,j),j=1,nl(icase))
else if(nm(icase).eq.98)then
write(text,'("NL=",i3,":",4i3)') nl(icase), (ll(icase,j),j=1,nl(icase))
! --------------------------------------------------
endif
if(nm(icase).lt.6)then ! for DOS write ordering of DOS on file 30
write(30,4031)ia
!pb put good text-header into case.qtl
text='tot,'
ltxt=5
do j=1,nL(icase)
call qtltext(LL(icase,j),nm(icase),text,ltxt)
enddo
write(30,*)
4031 format(' atom',i4,' ordering of projected DOS ')
endif
! JR: if((popmat.eq.87).or.(popmat.eq.98))then
! JR: WRITE(16,3031)jatom,MULT(jatom)
! JR: else
WRITE(16,3030)jatom,MULT(jatom),nm(icase),text
! JR: endif
icase=icase+1
else ! atom was not selected
text=' atom not selected for QTL calculation'
WRITE(16,3030)jatom,MULT(jatom),0,text(1:40)
endif
enddo
! JR: if(popmat.eq.0)then ! reorganize atoms qtl files to single file
if(popmat.ne.87 .and. popmat.ne.98)then ! reorganize atoms qtl files to single file
987 format(a5)
allocate (ieband(nband+1),kL1(nat))
! JR: we know already the number of bands and k-points!
! rewind(31)
! ib=0
! nene=0
! read(31,987)char
!143 continue
! nene=nene+1
! read(31,987,end=142)char
! if(char.eq.' BAND')then
! ib=ib+1
! ieband(ib)=nene/2
! nene=0
! endif
! goto 143
!142 continue
! ib=ib+1
! ieband(ib)=nene/2
! close(31)
! ------------------------------------------------------
do icase=1,natom
! jcase=30+icase
! rewind(jcase)
call mknam(FNAME,qtlfn,Icase)
open(1000+icase,FILE=FNAME,STATUS='old',FORM='formatted',ERR=951)
enddo
do ib=1,nband
write(16,3040)ib
! JR: do ie=1,ieband(ib)
! JR: the number of k-points is the same for all bands.... or?
do ie=1,nkpt
icase=1
sumlout=1.d0
do jatom=1,nat
kL=0
if(icase.le.natom)then
ia=iatom(icase)
else
ia=-1
endif
if(ia.eq.jatom)then ! atom was selected
if(ie.eq.1)read(1000+icase,987)char
! JR: merge also for qsplit=88 or 99
! JR: read(1000+icase,1050) Eb(icase),ja,kL,sumL(icase),(qtle(icase,jL),jL=1,kL)
! JR: kL1(icase)=kL
! JR: read(1000+icase,1050) Eb1,ja,kL,sumL1
208 FORMAT(10F10.5)
if(popmat.eq.0) then
read(1000+icase,1050) Eb(icase),ja,kL,sumL(icase),(qtle(icase,jL),jL=1,kL)
kL1(icase)=kL
read(1000+icase,1050) Eb1,ja,kL,sumi
else
read(1000+icase,1050) Eb(icase),ja,kL,sumL(icase),(qtle(icase,jL),jL=1,nl(icase))
if(popmat.eq.88) then
kL1(icase)=kL
elseif(popmat.eq.99) then
kL1(icase)=2*kL
else
stop 'ERROR: unknown popmat mode'
endif
read(1000+icase,208) (qtle1(jL),jL=1,kL1(icase))
read(1000+icase,1050) Eb1,ja,kL,sumi
endif
! --------------------------------
icase=icase+1
endif
enddo
icase=1
sumlout=1.d0
do jatom=1,nat
kL=0
if(icase.le.natom)then
ia=iatom(icase)
else
ia=-1
endif
if(ia.eq.jatom)then ! atom was selected
sumlout=sumlout-suml(icase)
kL=kL1(icase)
! JR: merge also for qsplit=88 or 99
! write(16,1051) Eb(icase),jatom,sumL(icase),(qtle(icase,jL),jL=1,kL)
if(popmat.eq.0) then
write(16,1051) Eb(icase),jatom,sumL(icase),(qtle(icase,jL),jL=1,kL)
else
write(16,1051) Eb(icase),jatom,sumL(icase),(qtle(icase,jL),jL=1,nl(icase))
write(16,208) (qtle1(jL),jL=1,kL)
endif
! ----------------------------------
icase=icase+1
else
sumi=0.
write(16,1051) Eb(1),jatom,sumi
endif
enddo
write(16,1051) Eb1,nat+1,sumLout
enddo
enddo
write(30,1052)nat+1
1052 format(' Data for interstital DOS correspond to atom index ',i4)
endif
1050 FORMAT(F10.5,I3,I5,F8.5,3X,39F8.5)
1051 FORMAT(F10.5,I3,F8.5,3X,39F8.5)
3000 FORMAT(A80,/)
3010 FORMAT(1X,'LATTICE CONST.=',3F8.4,3X,'FERMI ENERGY=',F10.5)
3020 FORMAT(I5,' < NMAT <',I5,3X,'SPIN=',I1,3X,'NAT=',I3, &
6x,'SO',i2,' KLmax',i3) !pb
3031 FORMAT(1X,'JATOM',I3,2X,'MULT=',I2,60x)
3030 FORMAT(1X,'JATOM',I3,2X,'MULT=',I2,2x,'ISPLIT=',i2,2x,a)
3040 FORMAT(1X,'BAND:',I4)
CALL CPUTIM(TTIME)
!.....CALCULATE CPUTIME REQUIRED
WRITE(6,2000)
WRITE(6,2010) TTIME
CALL ERRCLR(ERRFN)
STOP ' QTL END'
!
! error handling
!
910 INFO = 1
!
! 'QTL.def' couldn't be opened
!
WRITE (ERRMSG,9000) FNAME
CALL OUTERR('QTL',ERRMSG)
GOTO 999
920 INFO = 2
!
! file FNAME couldn't be opened
!
WRITE (ERRMSG,9010) IUNIT
CALL OUTERR('QTL',ERRMSG)
WRITE (ERRMSG,9020) FNAME
CALL OUTERR('QTL',ERRMSG)
WRITE (ERRMSG,9030) STATUS, FORM
CALL OUTERR('QTL',ERRMSG)
GOTO 999
951 INFO = 5
write(6,*) 'error:',FNAME
CALL OUTERR('QTL','file open error')
GOTO 999
960 INFO = 7
!
! Error reading file 'QTL.def'
!
WRITE (ERRMSG,9040) FNAME
CALL OUTERR('QTL',ERRMSG)
GOTO 999
999 STOP 'QTL - Error'
!
!
800 FORMAT(////,30X,50(1H-),/,33X,'S T R U C T U R A L ', &
'I N F O R M A T I O N',/,30X,50(1H-),//)
805 FORMAT(3X,'SUBSTANCE',20X,'= ',A80,/)
810 FORMAT(3X,'LATTICE',22X,'= ',A4)
820 FORMAT(3X,'LATTICE CONSTANTS ARE',8X,'= ',3F12.7)
840 FORMAT(3X,'NUMBER OF ATOMS IN UNITCELL = ',I3)
850 FORMAT(3X,'MODE OF CALCULATION IS',7X,'= ',A4)
1003 FORMAT(A4)
1004 FORMAT(2A5)
2000 FORMAT(//,3X,'=====>>> CPU TIME SUMMARY',/)
2010 FORMAT(12X,'TOTAL : ',F8.1)
9000 FORMAT('can''t open definition file ',A40)
9010 FORMAT('can''t open unit: ',I2)
9020 FORMAT(' filename: ',A50)
9030 FORMAT(' status: ',A,' form: ',A)
9040 FORMAT('Error reading file: ',A47)
END
--Boundary_(ID_mCpo2DQOnEfdFVv6S7Rfrw)--
More information about the Wien
mailing list