----< 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)--