I've made some small modifications of SRC_tetra, such that it can
calculate the DOS only for tetrahedra which contain the desired k-point.
a) Check case.klist and select the k-point of your choice (eg. Gamma as
number 1, or XX for any other selection.
b) modify case.int as: (3rd line)
PbTiO3 #Title
-0.91780 0.002 .9915901704 0.003 #Emin, DE, Emax, Gauss-Broad
11 N 0.000 KSEL=1 #Number of DOS-cases,G/L/B
broadening (Ry), KSELect=XX
Using the keyword KSEL= directs tetra to integrate only those tetrahedra
which contain the selected k-point.
A negative number (-XX) (or no KSEL=) switches this off and the regular
integration happens.
Depending on the density of your -mesh, it will cover a certain volume
of the BZ (a few percent ...).
NOTE: At the moment you cannot add this up for 2 different k-points,
since X1 and X2 may share the same tetrahedron and thus double counting
occurs. It would be straight forward, however, to specify several
k-points at the same time and use all tetrahedra which contain at least
one of the selected k-points ....
On 8/28/19 8:59 AM, Peter Blaha wrote:
Not with standard wien2k tools.
However, you can write your own little programs, eg.
reading the case.qtl file for an arbitrary k-mesh of your choice and do
a "root sampling", i.e. divide your energy mesh into some intervals and
count the number of eigenvalues (times the multiplicity of the k-point)
in each E-interval. Plot the corresponding histogram.
An alternative is to modify tetra such, that instead of summing over all
tetrahedra, it sums only over those involving k-points of your choice.
For this you need to understand the case.kgen file, which lists the
vertices of all tetrahedra and you select only thos which contain the
desired k-points.
On 8/27/19 6:12 PM, Jens Biegert wrote:
Dear All,
I was wondering wether I can get the DOS, or better PDOS, for a
specified k-path. Say, instead of an integration over the whole
Brillouin zone, just (like for the bandstructure) along gamma to K to
M to gamma.
If not possible, can I restrict 3D k-space integration close to e.g.
the K point such as to get the PDOS close to K?
I saw some previous answers that seem to hint at this not being
possible, but it was not clear to me.
Best regards,
Jens Biegert
_______________________________________________
Wien mailing list
[email protected]
http://zeus.theochem.tuwien.ac.at/mailman/listinfo/wien
SEARCH the MAILING-LIST at:
http://www.mail-archive.com/[email protected]/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: [email protected] WIEN2k: http://www.wien2k.at
WWW: http://www.imc.tuwien.ac.at/TC_Blaha
--------------------------------------------------------------------------
PROGRAM TETRA
!
! D E N S I T Y O F S T A T E S
! ACTUAL VERSION with BLOECHL SCHEME
!
use reallocate
INCLUDE 'param.inc'
CHARACTER *7 A7
CHARACTER *9 A9,B9
CHARACTER *20 TITLE
CHARACTER*120,allocatable :: TEXT(:)
CHARACTER *70 SYSTEM
CHARACTER*11 FORM,STATUS
CHARACTER*180 FNAME,fname1,fname2
CHARACTER*80 errfn
character*2 fnameupdn,fnum2
character*1 fnum1 ,abrsw
!
real*4,allocatable :: RNUMB(:,:), DENSTY(:,:),dsum(:)
real*4 fermden(MG)
COMMON /NCOUNT/ NNOC
real*4,pointer :: EBS(:,:), FC(:,:,:)
CHARACTER *14 DOSTYP(MG)
DIMENSION IDOS(MG,2)
real*4,allocatable :: QTL(:,:),xqtl(:,:)
real*4,allocatable :: ehelp(:),denbr(:,:)
integer,allocatable :: isumdos(:,:)
logical rxes,rxesw,enefile
common /rxescom/ rxes,rxesw
!
DATA PI /3.141592654/, EEF /0.0/
!--------------------------------------------------------------------
! call getarg(2,fname)
rxes=.false.
rxesw=.false.
enefile=.false.
CALL GTFNAM(fname,ERRFN)
! if(fname.eq.' ') call getarg(1,fname)
OPEN(1,FILE=fname,STATUS='OLD',ERR=8000)
8003 READ(1,*,END=8001) IUNIT,FNAME,STATUS,FORM,IRECL
OPEN(IUNIT,FILE=FNAME,STATUS=STATUS,FORM=FORM,ERR=8002)
if(iunit.eq.7) then
do i=180,5,-1
if(fname(i:i).eq.'1') then
nnchar=i-1
fname1(1:i-1)=fname(1:i-1)
fnameupdn=''
GOTO 8003
else if(fname(i:i).eq.'p') then
nnchar=i-3
fname1(1:i-3)=fname(1:i-3)
fnameupdn='up'
GOTO 8003
else if(fname(i:i).eq.'n') then
nnchar=i-3
fname1(1:i-3)=fname(1:i-3)
fnameupdn='dn'
GOTO 8003
endif
enddo
endif
if(iunit.eq.108) then
enefile=.true.
kspin=2
if(fnameupdn.eq.'') kspin=1
endif
if(iunit.eq.8) then
rxes=.true.
read(8,'(a)')system
write(6,'("RXES spectroscopy for",a )') system
endif
if(iunit.eq.9) then
rxesw=.true.
write(6,'("RXES spectroscopy weight file case.rxes calculated")')
endif
GOTO 8003
8000 WRITE(*,*) ' ERROR IN OPENING TETRA.DEF !!!!'
STOP 'TETRA.DEF'
8002 WRITE(*,*) ' ERROR IN OPENING UNIT:',IUNIT
WRITE(*,*) ' FILENAME: ',FNAME,' STATUS: ',STATUS,' FORM:',FORM
STOP 'OPEN FAILED'
8001 CONTINUE
CLOSE (1)
!----------------------------------------------------------------------
READ(5,1) TITLE
IAV=0
NPRINT=1
! MI=1000
eef=999.
! EEF=0.
! NST number of k points
! IAV 0 NO PLOT OTHERWISE DOS AVERAGED OVER IAV VALUES
! NPRINT 0 NO DOS PRINTED
! 1 DOS AND INTEGRATED DOS PRINTED
! NYMIN LOWER BAND INDEX
! NYMAX UPPER BAND INDEX
! ZEL NUMBER OF ELECTRONS
! EMIN,DE,EMAX ENERGY GRID FOR DOS
!
! DOS1 DOS: GES, S, P, D, EG, T2G, F, FOR ATOM NR.1
! DOS2 DOS: S, ... F, OUT FOR ATOM NR.2
!
!***** FIRST PRINTED ENERGY ON DOS PLOT IS
!***** E0 + (IAV-1) * DE / 2 + 4 * DE * IAV
!***** E.G.IAV=5 FIRST ENERGY E0 + 22*DE
!***** NEXT AT 25*DE HIGHER
!
! read(5,*) nst
read(3,*) nst
rewind 3
! if(nst.gt.nktot) then
! write(*,*) ' nktot',nktot,' .lt. nst=',nst
! stop
! end if
! READ(5,*) NYMIN,NYMAX
! read(5,*) ZEL
nymin=1
! nymax=nbtot
prev=0.d0
! prev: Gaussian broadening parameter in Ry!!!
read(5,*,ERR=4848) EMIN,DE,EMAX,prev
4848 continue
! if(nymax.gt.nbtot) then
! write(*,*) ' nbtot',nbtot,' .lt. nymax=',nymax
! stop
! end if
ibrsw=-1 !use tetra by default
br=0.01
brlor=0.01
read(5,'(a)') system
READ(system,*,end=123,err=123) NNDOS,abrsw,br
if(abrsw.eq.'G') then
ibrsw=1
else if(abrsw.eq.'L') then
ibrsw=0
brlor=br
else if(abrsw.eq.'B') then
ibrsw=2
READ(system,*,err=123,end=123) NNDOS,abrsw,br,brlor
endif
123 continue
do i=1,65
if(system(i:i).eq.'#') exit
if(system(i:i+4).eq.'KSEL=') then
read(system(i+5:70),*) kselect
print*,'DOS for tetrahedra around K=',kselect
exit
endif
enddo
if(ibrsw.ge.0) then
print*,'using broadening method instead of TETRA',ibrsw,br,brlor
write(6,*)'using broadening method instead of TETRA',ibrsw,br,brlor
endif
if(nndos.gt.mg) stop 'mg too small, TETRA does not support more the 51 cases'
!
if(nndos.gt.1 .and. enefile) then
STOP 'NNDOS must be 1 when using case.energy file'
endif
! reading case.energy instead of case.qtl
if(enefile) then
779 read(109,'(A)',end=997) system
if(system(1:4).ne.':FER') goto 779
do i=1,69
if(system(i:i).eq.'=') exit
enddo
read(system(i+1:70),*) EFERM
eef=EFERM
spin=1.d0
if(KSPIN.eq.2) spin=0.5d0
KSO=0
SO=1.d0
NSORT=1
KLMAX=1
allocate( TEXT(0:Nsort+1),QTL(Nsort+1,klmax+1),xqtl(nsort,lxdos2*(lxdos2+1)) )
IDOS(1,1)=0
IDOS(1,2)=0
DOSTYP(1)='TOTAL'
NRBAND=99999999
emax1=-999.d0
nbtot=1000
allocate (EBS(Nst,NBTOT), FC(MG,Nst,NBTOT))
do k=1,nst
777 read(108,'(3e19.10,i10,6x,i6)',err=777,end=999) xk,xk,zk,numk,ne
NRBAND=min(NRBAND,ne)
if(ne.gt.nbtot) then
nbtot=ne*1.5
call doreallocate(ebs, nst, nbtot)
call doreallocate(fc, mg, nst, nbtot)
endif
do i=1,ne
read(108,*) i1,ebs(numk,i)
enddo
enddo
do i=1,NRBAND
BMIN=10000.
BMAX=-10000.
do k=1,nst
BMIN=AMIN1(ebs(k,i),BMIN)
BMAX=AMAX1(ebs(k,i),BMAX)
enddo
WRITE(6,14) i,BMIN,BMAX
emax1=amax1(bmin,emax1)
enddo
nymax=NRBAND
if(emax.gt.emax1) then
write(6,*) ' EMAX reduced due to lower HIGHEST BAND-minimum'
emax=emax1
endif
!
fc=1.d0*spin
endif
! end reading case.energy
WRITE(6,3) TITLE
WRITE(6,4) IAV,NPRINT
! WRITE(6,5) NYMIN,NYMAX
WRITE(6,7) NNDOS
!
nnsum_dos=0
if(enefile) goto 200 ! read case.energy instead of case.qtl
REWIND 4
READ(4,1000) SYSTEM
READ(4,1010) ALAT1,ALAT2,ALAT3,EFERM
READ(4,1020) MINWAV,MAXWAV,KSPIN,NSORT,kso,kLmax
if(klmax.eq.0) klmax=14
allocate( TEXT(0:Nsort+1),QTL(Nsort+1,klmax+1),xqtl(nsort,lxdos2*(lxdos2+1)) )
! INCREASE NSORT BY ONE FOR THE INTERSTITIAL DOS (OUT)
NSORT=NSORT+1
WRITE(6,2000) SYSTEM
WRITE(6,2010) ALAT1,ALAT2,ALAT3,EFERM
WRITE(6,2020) MINWAV,MAXWAV,KSPIN,NSORT
eef=eferm
SPIN=1.d0
IF(KSPIN.EQ.2) SPIN=0.5d0
so=1.d0
if(kso.eq.2) then
so=2.d0
spin=spin/2.d0
endif
!
! 106 KAT=NST*NSORT
!
DO 107 ISORT=1,NSORT
IF(ISORT.EQ.NSORT) GOTO 107
READ(4,1030) JATOM,NEQU,isplit,text(isort)
WRITE(6,2030) JATOM,NEQU,isplit,text(isort)
107 CONTINUE
nndos_save=nndos
DO 90 JCASE=1,NNDOS
READ(5,*,err=998) IATOM,ICOLUM !,DOSTYP(JCASE)
IDOS(JCASE,1)=IATOM
IDOS(JCASE,2)=ICOLUM
call parsetext(iatom,icolum,dostyp(jcase),text(iatom))
IF(IATOM.EQ.0) THEN
IDOS(JCASE,2)=0
ENDIF
WRITE(6,25) JCASE,IDOS(JCASE,1),IDOS(JCASE,2),DOSTYP(JCASE)
90 CONTINUE
nnsum_dos=0
read(5,'(a)',end=91) system
if(system(1:3).ne.'SUM') goto 91
read(system(5:70),*,ERR=91,END=91) nnsum_dos,nnsum_dos_max
if(nnsum_dos.gt.7) then
nnsum_dos=7
print*,':WARNING: maximal NNSUM_DOS=7, input truncated'
endif
WRITE(6,*) 'We will add ',nnsum_dos,' DOS-cases together:'
allocate(isumdos(nnsum_dos,nnsum_dos_max),dsum(nnsum_dos))
isumdos=0
do i=1,nnsum_dos
read(5,'(a)') system
READ(system,*,end=92,err=92) (isumdos(i,j),j=1,nnsum_dos_max)
92 continue
write(6,'("adding:",99i3)') (isumdos(i,j),j=1,nnsum_dos_max)
enddo
91 continue
! in spin.pol SO case set total DOS to DOS of interstital (=normso)
if(kso.eq.2.and.kspin.eq.2) then
DO JCASE=1,NNDOS
if(IDOS(JCASE,1).eq.0) then
IDOS(JCASE,1)=NSORT
IDOS(JCASE,2)=1
ENDIF
enddo
endif
text(nsort)=' '
!
emax1=-999.d0
nbtot=1000
allocate (EBS(Nst,NBTOT), FC(MG,Nst,NBTOT))
ny=0
115 NY=NY+1
if(ny.gt.nbtot) then
nbtot=nbtot*2
call doreallocate(ebs, nst, nbtot)
call doreallocate(fc, mg, nst, nbtot)
endif
read(4,'(a)',end=141) title
READ(title,1041,err=1042) NRBAND
goto 1043
1042 READ(title,1040) NRBAND
1043 continue
BMIN=10000.
BMAX=-10000.
DO 120 I=1,NST
DO 113 ISORT=1,NSORT
READ(4,1050,err=114) EVAL,ISRT,( QTL(ISORT,J), J=1,klmax+1 )
! WRITE(6,1050) EVAL,ISRT,( QTL(ISORT,J), J=1,13 )
if(text(isort)(13:20).eq.'xdos(i,j') then
if(lxdos.ne.3) stop 'LXDOS must be 3 for ISPLIT 99 or 88'
read(4,208)(xqtl(isort,i1),i1=1,(lxdos2)*(lxdos2+1))
else if(text(isort)(13:20).eq.'xdos(i,i') then
read(4,208)(xqtl(isort,i1),i1=1,lxdos2)
endif
goto 113
114 print*,'qtl-reading error',eval,isrt,( QTL(ISORT,J), J=1,klmax+1 ),'BAND',nrband,'K=',i,'ISORT=',isort
stop 'error reading qtl-file'
113 CONTINUE
EBS(I,NY)=EVAL
BMIN=AMIN1(EVAL,BMIN)
BMAX=AMAX1(EVAL,BMAX)
DO 201 NDOS=1,NNDOS
IAT=IDOS(NDOS,1)
JAT=IDOS(NDOS,2)
if(iat.gt.NSORT) stop 'iat in input too large'
IF(IAT.EQ.0) FC(NDOS,I,NY)=1.d0*SPIN
! IF(IAT.GT.0.and.jat.le.14)then
IF(IAT.GT.0.and.jat.le.klmax+1)then
FC(NDOS,I,NY)=QTL(IAT,JAT)*SPIN*so
else if(IAT.GT.0) then
FC(NDOS,I,NY)=xQTL(IAT,JAT-100)*SPIN
endif
201 continue
120 CONTINUE
!
WRITE(6,14) NRBAND,BMIN,BMAX
emax1=amax1(bmin,emax1)
goto 115
141 nymax=ny-1
if(emax.gt.emax1) then
write(6,*) ' EMAX reduced due to lower HIGHEST BAND-minimum'
emax=emax1
endif
200 CONTINUE
!
JEND=1 + (EMAX-EMIN)/DE
if(jend.lt.1) stop 'Error: Energy mesh not ok'
allocate (RNUMB(jend,MG), DENSTY(jend,MG))
allocate (ehelp(jend),denbr(jend,mg))
DENSTY=0.0
RNUMB=0.0
write(6,55) EMIN,DE,EMAX
if(rxesw) then
! if(jend.ne.1) then
! print*, "For RXESW switch exactly ONE energy must be specified in case.int"
! stop 'Error: RXESW-select only one energy'
! endif
! if(nndos.ne.1) then
! print*, "For RXESW switch only ONE DOS must be specified in case.int"
! stop 'Error: RXESW-select only one DOS'
! endif
write(9,'("Energy",2f6.3," atom,column",7(i4,i3))') emin,emax,(IDOS(j,1),IDOS(j,2),j=1,nndos)
write(6,'(/,"RXES weight file case.rxes calculated", &
& " for Atom",i4," col=",i3," Energy=",2f8.4,/)') &
IDOS(1,1),IDOS(1,2),emin,emax
endif
CALL ARBDOS (0,NNDOS,1,JEND,EMIN,EMAX,NYMIN,NYMAX,DE,ebs,fc,nst,rnumb,densty,kselect)
if(nst.eq.1.and.ibrsw.eq.-1) then ! for one k-point calculate the DOS with broadening
ibrsw=1
if(br.lt.1.d-6) then
br=0.002
! print*,' No broadening specified, increasing broadening to',br
endif
print*,'Only 1 k-point, switching to Gauss-broadening with',br
Write(6,*) 'Only 1 k-point, switching to Gauss-broadening with',br
endif
if(ibrsw.ne.-1) then
CALL broadDOS (ibrsw,br,brlor,NNDOS,1,JEND,EMIN,EMAX,NYMIN,NYMAX,DE,ebs,fc,nst,densty)
endif
!
!.....WRITE OUT RESULTS TO DOS FILES
IFILE=6
ilmin=1
do 295 j=1,jend
if(emin+(j-1)*de.ge.eef) then
if(idos(1,1).eq.0) zel=rnumb(j,1)
jef=j
goto 300
endif
295 continue
300 CONTINUE
IFILE=IFILE+1
if(ifile.gt.7) then
close(ifile-1)
close(ifile-1+50)
fname=''
fname2=''
if(ifile-6.lt.10) then
write(fnum1,'(i1)') ifile-6
fname=fname1(1:nnchar)//fnum1//fnameupdn
fname2=fname1(1:nnchar)//fnum1//'ev'//fnameupdn
else
write(fnum2,'(i2)') ifile-6
fname=fname1(1:nnchar)//fnum2//fnameupdn
fname2=fname1(1:nnchar)//fnum2//'ev'//fnameupdn
endif
print*,'openfilename:',fname
OPEN(ifile,FILE=fname,STATUS='unknown')
OPEN(ifile+50,FILE=fname2,STATUS='unknown')
endif
IL=MIN0(7,NNDOS)
ilmax=il+ilmin-1
WRITE(IFILE,15) TITLE,EEF,IL,JEND,prev
WRITE(IFILE,217) (DOSTYP(IX),IX=ilmin,ILmax)
WRITE(IFILE+50,15) TITLE,0.0,IL,JEND,prev
WRITE(IFILE+50,217) (DOSTYP(IX),IX=ilmin,ILmax)
WRITE(6,15) TITLE,EEF,IL,JEND,prev
if(idos(1,1).eq.0) then
! write(6,44) zel
if(spin.eq.1.d0) then
write(6,*) 'DOS in states/Ry'
else
write(6,*) 'DOS in states/Ry/spin'
endif
endif
!
DO J=1,JEND
Ehelp(j)=EMIN+(J-1)*DE
enddo
do IND=ILMIN,ILMAX
call Broad(prev,jend,ehelp,DENSTY(1,IND),denbr(1,ind))
enddo
write(6,44) (rnumb(jef,ix),ix=ilmin,ilmax)
WRITE(6,117) (IDOS(IX,1),DOSTYP(IX),IX=ilmin,ILmax)
!
DO J=1,JEND
ENRG=EMIN+(J-1)*DE
WRITE(6,16) ENRG,(DENSTY(J,IND)*2.d0,RNUMB(J,IND),IND=ILMIN,ILMAX)
WRITE(IFILE,116) ENRG,(DENbr(J,IND)*2.d0,IND=ILMIN,ILMAX)
WRITE(IFILE+50,116) (ENRG-eef)*13.605698,(DENbr(J,IND)/13.605698*2.d0,IND=ILMIN,ILMAX)
IF(ENRG.GE.eferm.AND.enrg-DE.LE.eferm.and.j.gt.1) THEN
DO ind=ilmin,ilmax
fermden(ind)=2.d0*((eferm-enrg+de)*DENSTY(J,IND)+(enrg-eferm)*DENSTY(J-1,IND))/de
ENDDO
ENDIF
ENDDO
240 CONTINUE
NNDOS=NNDOS-7
WRITE(6,*) '******** EF and DOS at fermi level *******'
WRITE(6,'(f9.5,7(F9.2,9x))') eferm,(fermDEN(ind),ind=ILMIN,ilmax)
! AM (2.80)
DO ind=ilmin,ilmax
! fermDEN(ind)=fermDEN(ind)/2.d0*(8.617D-5)**2*(pi**2)/3.d0*2.6255D6/27.212d0*1000.d0
fermDEN(ind)=fermDEN(ind)*(pi**2)/3.d0*(6.3336303d-6)**2*6.0221d23*13.605698d0*1.60219d-19*1000.d0
ENDDO
WRITE(6,*) 'Gamma in mJ/(mol cell K**2). (Divide by number of formula units in cell to get it per mole only)'
WRITE(6,'(A9,7(F9.2,9x))') 'Cv/T ',(fermDEN(ind),ind=ILMIN,ilmax)
ilmin=ilmin+7
IF(NNDOS.GT.0) GOTO 300
! summing up some DOS columns
if(nnsum_dos.gt.0) then
do i=1,nnsum_dos
WRITE(6,1176) i,(isumdos(i,i1),i1=1,nnsum_dos_max)
write(dostyp(i),'(7i2)') (isumdos(i,i1),i1=1,min(7,nnsum_dos_max))
enddo
WRITE(6,1177) (i1,i1=1,nnsum_dos)
fname=fname1(1:nnchar)//'sum'//fnameupdn
fname2=fname1(1:nnchar)//'sumev'//fnameupdn
print*,'openfilename:',fname
OPEN(ifile,FILE=fname,STATUS='unknown')
OPEN(ifile+50,FILE=fname2,STATUS='unknown')
WRITE(IFILE,15) TITLE,EEF,IL,JEND,prev
WRITE(IFILE,217) (DOSTYP(IX),IX=1,nnsum_dos)
WRITE(IFILE+50,15) TITLE,0.0,IL,JEND,prev
WRITE(IFILE+50,217) (DOSTYP(IX),IX=1,nnsum_dos)
DO J=1,JEND
ENRG=EMIN+(J-1)*DE
dsum=0.d0
do i=1,nnsum_dos
do i1=1,nnsum_dos_max
if(isumdos(i,i1).eq.0) exit
! dsum(i)=dsum(i)+DENSTY(J,isumdos(i,i1))*2.d0
dsum(i)=dsum(i)+DENbr(J,isumdos(i,i1))*2.d0
enddo
enddo
WRITE(6,116) ENRG,(Dsum(i),i=1,nnsum_dos)
WRITE(IFILE,116) ENRG,(Dsum(i),i=1,nnsum_dos)
WRITE(IFILE+50,116) (ENRG-eef)*13.605698,(Dsum(i)/13.605698,i=1,nnsum_dos)
enddo
endif
!DOS at Fermi
STOP ' LEGAL END TETRA'
997 print*, "ERROR when reading EF from case.scf2"
STOP 'Error when reading case.scf2'
998 print*, "ERROR when reading the DOS-cases in case.int"
STOP 'Error when reading case.int'
999 print*, "ERROR when reading case.energy"
STOP 'Error when reading case.energy'
!
1 FORMAT(A20)
! 1 FORMAT(A20/A4,/I2/I2)
2 FORMAT(2I3/F10.5/3F10.5)
3 FORMAT(1X,A20/)
4 FORMAT(1X,'IAV',25X,':',I3,/,1X,'NPRINT',22X,':',I3)
44 Format(/,1X,'NUMBER OF ELECTRONS UP TO EF:',/,':CTOTAL: ',7(F10.4,9x))
5 FORMAT(1X,'LOWER AND UPPER BAND-INDEX',2X,':',2I5)
55 format(1X,'EMIN, DE, EMAX',14X,':',3F10.5,/)
6 FORMAT(I5)
7 FORMAT(1X,I2,' CASES FOR DOS',12X,': ATOM L')
13 FORMAT(8X,I2)
14 FORMAT(' BAND LIMITS OF BAND',I4,' ARE',2F10.5)
15 FORMAT('# ',A20,/,'#EF=',F10.5,5X,'NDOS=',I2,5X,'NENRG=',I5,' Gaussian bradening:',f8.5)
16 FORMAT(f9.5,7(F9.2,f9.4))
116 FORMAT(f10.5,7f14.8)
117 FORMAT('# ENERGY',1X,7(1x,I3,1X,A14))
1175 FORMAT(//,'DOS summed up for following columns:')
1176 FORMAT('Sum',i3,':',1X,7I3)
1177 FORMAT('# ENERGY',1X,7(5x,I3,6X))
217 FORMAT('# ENERGY',2X,7(A14,1x))
18 FORMAT(' ATOM NO.',I3,4X,A6,'-CHARGE',F10.5)
19 FORMAT(A72/17X,3F9.5,18X,F7.4/5X,I3,9X,I2,11X,I3)
20 FORMAT(A7,F4.1,1X,A9,F8.5,A9,I3)
221 FORMAT(F10.5,I3,1X,F9.5,3X,6F9.5)
21 FORMAT(F10.5,I4,1X,F10.5,5X,12F10.5)
22 FORMAT(/,1X,A72/17X,3F9.5,18X,F7.4/5X,I3,9X,I2,11X,I3)
23 FORMAT(1X,A7,F4.1,1X,A9,F8.5,A9,I3)
24 FORMAT(3X,I2,2X,I3,3X,A6)
25 FORMAT(1X,'CASE',I2,' : ATOM NUMBER',I3,' COLUMN READ',I3,' DOSTYPE=',A14)
208 FORMAT((10F10.5))
1000 FORMAT(1X,A70,/)
1010 FORMAT(1X,15X,3F8.5,16X,F10.5)
1020 FORMAT(1X,I4,10X,I4,8X,I1,7X,I3,8x,i2,6x,i3)
1030 FORMAT(6X,I3,7X,I2,9x,i2,2x,a120)
1040 FORMAT(6X,I4)
1041 FORMAT(5X,I5)
1050 FORMAT(F10.5,I3,F8.5,3X,(40F8.5))
2000 FORMAT(/,1X,A70)
2010 FORMAT(1X,'LATTICE CONST.=',3F8.5,3X,'FERMI ENERGY=',F10.5)
2020 FORMAT(1X,I4,' < NMAT < ',I4,3X,'SPIN=',I1,3X,'NATO=',I4)
2030 FORMAT(1X,'JATOM',I3,2X,'MULT=',I2,2x,'ISPLIT=',i2,2x,a120)
END
!.......................................................................
SUBROUTINE ADOS (NS,NFU,NEMIN,NEMAX,ebs,fc,nst,RNUMB,DENSTY,SNUMB,nlast,kselect)
! ** ARBDOS READS EIGENVALUES AND EIGENFUNCTIONS FROM FILE 8 **
! ** **
! ** SUBROUTINES USED: **
! ** NOCULC **
! ** **
! ** INPUT AND OUTPUT: **
! ** NFU: NUMBER OF EIGENFUNCTIONS PER K-POINT *
! ** NNOC: NUMBER OF ALLREADY FILLED TETRAHEDRA **
! ** NEMAX: NUMBER OF LAST BAND **
! ** NEMIN: NUMBER OF FIRST BAND-1 **
! ** NUMK: NUMBER OF K-POINTS **
! ** NTT NUMBER OF DIFFERENT TETRAHEDRA **
! ** **
INCLUDE 'param.inc'
COMMON /NCOUNT/ NNOC
real*4 EBS(Nst,*), FC(MG,Nst,*)
real*4 RNUMB(nlast,MF), DENSTY(nlast,MF)
real*4 SNUMB(nlast,MF)
COMMON /EMICRO/ D(4), F(MF,4), V
logical rxes,rxesw
common /rxescom/ rxes,rxesw
real*4, allocatable :: rxesweight(:,:)
real*4, allocatable :: dosold(:,:),dosold1(:,:),dosnew(:)
logical doit
PARAMETER (NITT = 505)
DIMENSION ITTFL(nitt)
DATA IWTET /0/
NNOC=0
NWRIT=100
READ(3,1234)Ndim,NTT,V1,nWRIT,NREC
1234 format(2i10,e20.12,2i10)
! READ(12,1101) NDIM,NTT,V1
WRITE (6,*) 'NUMBER OF K-POINTS:',NDIM
write(6,*) 'NUMBER OF TETRAHEDRONS:',NTT
if(rxesw ) then
allocate (dosold1(nlast,ns+nfu),dosnew(ns+nfu))
dosold1=0.0
endif
if(rxes ) then
allocate (rxesweight(ntt,ns+nfu))
allocate (dosold(nlast,ns+nfu),dosold1(nlast,ns+nfu),dosnew(ns+nfu))
dosold=0.0
dosold1=0.0
read(8,*) (j1,j2,(rxesweight(j,j3),j3=1,ns+nfu),j=1,ntt)
endif
! NTTW=(NTT/NWRIT)+1
!.....correction by Nelhiebel
if((ntt/nwrit)*nwrit.eq.ntt) then
Nttw=NTT/NWRIT
else
Nttw=NTT/NWRIT+1
end if
!.....
ind=0
dostot=0.d0
volumetot=0.d0
volumepart=0.d0
DO 30 N=1,NTTW
IF (N.EQ.NTTW) THEN
KMAX=NTT-(N-1)*NWRIT
ELSE
KMAX=NWRIT
ENDIF
!cad
if(5*kmax.gt.NITT) then
write(*,*) 'increase parameter NITT in ados!'
write(*,*) 'the value should be at least ',5*kmax
STOP 'NITT in ados too small'
endif
!cad
! READ(12,1102) (ITTFL(K),K=1,KMAX*5)
READ(3,1235) (ITTFL(K),K=1,KMAX*5)
1235 format(6i10)
DO 30 K=1,KMAX
V=FLOAT(ITTFL((K-1)*5+1))*V1
ind=ind+1
doit=.true.
if(kselect.gt.0) doit=.false.
DO 31 JJ=NEMIN,NEMAX
DO 20 KP=1,4
KPP=ITTFL((K-1)*5+KP+1)
if(kpp.eq.kselect.and.kselect.gt.0) doit=.true.
DO 10 I=1,NFU
10 F(I,KP)=FC(I,KPP,JJ)
20 D(KP)=EBS(KPP,JJ)
if(jj.eq.nemin) volumetot=volumetot+v
if(.not.doit) exit
if(jj.eq.nemin.and.kselect.gt.0) print*,'doing tetrahedra',(n-1)*nwrit+k
if(jj.eq.nemin) volumepart=volumepart+v
CALL NOCULC (NS, RNUMB,DENSTY,SNUMB,nlast)
if(rxesw .and. (jj.eq.nemax)) then
do l=1,nfu
dosnew(l)=0.0
do i=1,nlast ! loop over energy-mesh
dosnew(l)=dosnew(l)+densty(i,l+ns)-dosold1(i,l)
dosold1(i,l)=densty(i,l+ns)
enddo
dostot=dostot+dosnew(l)
enddo
write(9,'(2i6,7e20.12)') n,k,(dosnew(l),l=1,nfu)
endif
if(rxes) then
DO L=1,NFU ! loop over cases
do i=1,nlast ! loop over energy-mesh
! if(i.eq.600) print*,n,k,jj,rxesweight(ind,1),dosold(i,1)
dosnew(l)=densty(i,L+NS)-dosold(i,L+NS)
densty(i,L+NS)=dosold(i,L+NS)+dosnew(l)*rxesweight(ind,L+ns)
dosold(i,L+NS)=densty(i,L+NS)
enddo
enddo
endif
31 CONTINUE
30 CONTINUE
print*,'covered volume (%)',volumepart/ volumetot*100.d0
write(6,*) 'dostot in rxesw', dostot
RETURN
END
SUBROUTINE ARBDOS (NFLOW,NFUN1,NFIRST,NLAST,EMIN,EMAX,NEMIN,NEMAX &
,DET,ebs,fc,nst,RNUMB,DENSTY,kselect)
! NEW TETRA.-SUBR. (P.BLOECHL, OKA, O. JEPSEN, M. ALOUANI
!**********************************************************************
!* *
!* DENSITY OF STATES PROGRAM *
!* CALCULATES THE PROJECTED DENSITIES OF STATES, CORRESPONDING *
!* TO THE BAND ENERGIES AND ANGULAR MOMENTUM WEIGHTS, STORED ON *
!* THE FILE 8 *
!*********************************************************************
INCLUDE 'param.inc'
real*4 EBS(Nst,*), FC(MG,Nst,*)
real*4 RNUMB(nlast,mg), DENSTY(nlast,mg)
real*4,allocatable :: SNUMB(:,:)
COMMON /EME/ EMIN1, EMAX1, EFACTR, ESTEP, NO, NFUN, NU
!
allocate (SNUMB(nlast,MG))
snumb=0.d0
ESTEP=DET
EFACTR=1.E0/DET
NFUN=NFUN1
EMIN1=EMIN
EMAX1=EMAX
NU=NFIRST
NO=NLAST
WRITE (6,10) EMIN,EMAX,EFACTR,ESTEP,NEMIN,NEMAX,NU,NO
10 FORMAT (1H ,' EMIN=',F10.5,' EMAX=',F10.5,' EFACTR=',F16.8, &
' ESTEP =',F10.5/' ENERGY BAND',I5,' THROUGH',I5,' ENERGY CHANNEL:', &
I5,' TO',I5)
CALL ADOS (NFLOW,NFUN1,NEMIN,NEMAX,ebs,fc,nst, RNUMB,DENSTY,SNUMB,nlast,kselect)
!
NU1=NU+1
DO 30 L=1,NFUN1
DO 20 I=NU1,NO
SNUMB(I,L+NFLOW)=SNUMB(I,L+NFLOW)+SNUMB(I-1,L+NFLOW)
20 RNUMB(I,L+NFLOW)=(RNUMB(I,L+NFLOW)+SNUMB(I,L+NFLOW))*2.
30 CONTINUE
! DO 40 L=1,NFUN1
! DO 40 I=NU,NO
! TDOS(I)=TDOS(I)+DENSTY(I,L+NFLOW)
! 40 TNOS(I)=TNOS(I)+RNUMB(I,L+NFLOW)
RETURN
! 50 WRITE (6,60) Nlast
! 60 FORMAT (' ARRAYS TOO SMALL. Nlast=',I7)
! STOP
END
_______________________________________________
Wien mailing list
[email protected]
http://zeus.theochem.tuwien.ac.at/mailman/listinfo/wien
SEARCH the MAILING-LIST at:
http://www.mail-archive.com/[email protected]/index.html
