Dear Siesta users,in view of some modifications of the PDOS format in recent
versions of Siesta,my script fmpdos.f (for extraction of atom-resolved partial
densities of states)required modification. The updated version (hopefully able
to manage bothold and new formats) is in attachment, and also available
athttp://www.home.uni-osnabrueck.de/apostnik/Software/fmpdos.fThank you very
much in advance for testing and sending me bug reports. Best regardsAndrei
Postnikov----- James Sifuna <sifunaja...@gmail.com> a écrit :
>Hello,
I have successfully done the dos and pdos calculations,
I need to extract the pdos of atom 1 with n=4, l=0, m=0 etc. but i get an
error..
siesta-4.0.2/Util/Contrib/APostnikov/fmpdos
Input file name (PDOS):
SrTiO3.PDOS
Output file name :
Sr.dat
Extract data for atom index (enter atom NUMBER, or 0 to select all),
or for all atoms of given species (enter its chemical LABEL):
1
Extract data for n= ... (0 for all n ):
0
Unknown identifier in PDOS file, line 6
-0.699999252E+02
Best wishes,
J.
program fmPDOS
!
! extracts partial densities of states from PDOS file of SIESTA,
! taking care of flexible format in index, atom_index, etc. fields
! Andrei Postnikov, University Paul Verlaine - Metz (Jul 2010)
! University of Lorraine (Jan 2019)
! andrei.postni...@univ-lorraine.fr
!
implicit none
integer ii1,io1
parameter (ii1=11,io1=12)
integer nt,nmax,i0,ispin,nspin,norbs,it,npts,nene,is,idos,l,lref,
. m,mode,mref,zeta,zvalue,atind,n,nref,indref,index,nline,
. iquoted,iparsed
parameter (nmax=10000) ! max. number of energy values in DOS
double precision ene(nmax),dos(nmax,4),dos1(4),dparsed,efermi
character inpfil*60,outfil*60,string*80,llabel*80,rlabel*80,
. species*6,squoted*6,chlab*6,owrite*1,s1*1,polar*6,
. dtail*7
logical filexist,addos,nptdef,ene_fo,isferm
external iquoted,squoted,iparsed,dparsed
do it=1,nmax
do is=1,2
dos(it,is)=0.d0
enddo
enddo
nptdef = .false. ! We don't know whether the file contains <npoints>
isferm = .false. ! Show the Fermi energy if included in the PDOS
! --- open input and outpute files: ----------------------
write (6,*) ' Input file name (PDOS):'
read (5,*) inpfil
inquire (file=trim(inpfil), exist=filexist)
if (filexist) then
open (ii1,file=trim(inpfil),status='old',form='formatted')
else
write (6,*) ' File does not exist!'
stop
endif
write (6,*) ' Output file name :'
read (5,*) outfil
inquire (file=trim(outfil), exist=filexist)
if (filexist) then
write (6,*) ' File ',trim(outfil),' exists. Overwrite? (Y/N)'
read (5,*) owrite
if (owrite.eq.'Y'.or.owrite.eq.'y') then
open (io1,file=trim(outfil),form='formatted',status='REPLACE')
else
write (6,*) ' Then rename is first. Bye...'
stop
endif
else
open (io1,file=trim(outfil),form='formatted',status='NEW')
endif
write (io1,309) trim(inpfil)
! --- select projected DOS which are needed: -------------
100 continue
write (6,*) ' Extract data for atom index',
. ' (enter atom NUMBER, or 0 to select all),'
write (6,*) ' or for all atoms of given species',
. ' (enter its chemical LABEL):'
read (5,*) string
! check first character of string, whether it is decimal
i0 = ichar(string(1:1))
if (i0.ge.48.and.i0.le.57) then
! 1st character passed is a number;
! assume the whole string is atom number
mode = 1
read (string,*,err=401) indref
write (io1,301)indref
elseif ((i0.ge.65.and.i0.le.90).or.(i0.ge.97.and.i0.le.122)) then
! 1st character passed is a character [A-Z,a-z];
! assume the whole string is atom label
mode = 2
read (string,*,err=402) species
write (io1,302) trim(species)
else
! 1st character is a special symbol; presumably an error
write (6,306) i0
goto 100
endif
write (6,*) ' Extract data for n= ... (0 for all n ):'
read (5,*) nref
if (nref.ne.0) then
write (6,*) ' Extract data for l= ... (-1 for all l ):'
read (5,*) lref
if (lref.ne.-1) then
write (6,*) ' Extract data for m= ... (9 for all m ):'
read (5,*) mref
else
mref=9
endif
else
lref=-1
mref=9
endif
write (io1,308) nref,lref,mref
ene_fo = .false. ! are we in the energy reading mode ?
addos = .false. ! add this DOS to the result ?
nene=0 ! counter for energy values
nline = 0 ! global conter for lines in PDOS file
! --- read and analyze the PDOS file line by line: -----------
11 continue
nline = nline+1
read (ii1,'(a80)',err=201,end=202) string
if (string(1:6).eq.'<pdos>') goto 11
if (string(1:7).eq.'<nspin>') then
llabel = '<nspin>'
rlabel = '</nspin>'
nspin = iparsed(string,llabel,rlabel)
if (nspin.eq.8) then
write (6,*) ' nspin=8 changed to nspin=4'
nspin=4
endif
goto 11
elseif (string(1:11).eq.'<norbitals>') then
llabel = '<norbitals>'
rlabel = '</norbitals>'
norbs = iparsed(string,llabel,rlabel)
goto 11
elseif (string(1:13).eq.'<fermi_energy') then
llabel = '<fermi_energy units="eV">' ! no other options
! than units="eV" is so far provided in pdos.F
rlabel = '</fermi_energy>'
efermi = dparsed(string,llabel,rlabel)
isferm = .true.
goto 11
elseif (string(1:9).eq.'<npoints>') then
llabel = '<npoints>'
rlabel = '</npoints>'
nptdef = .true.
npts = iparsed(string,llabel,rlabel)
goto 11
elseif (string(1:14).eq.'<energy_values') then
ene_fo = .true. ! Ready for reading energy values that must follow
nene=0
goto 11
elseif (string(1:16).eq.'</energy_values>') then
ene_fo = .false. ! List of energies closes
nt=nene
if (nt.gt.nmax) then
write(6,*)' nt=',nt,' > nmax=',nmax,', must be increased'
stop
endif
if (nptdef.and.nt.ne.npts) then
write(6,*)' nt=',nt,' differs from npoints=',npts
stop
endif
goto 11
elseif (string(1:8).eq.'<orbital') then ! new orbital follows:
addos = .false. ! will be set to .true. below if needed
goto 11
elseif (string(2:7).eq.'index=') then ! orbital index line,
index = iquoted(string,nline) ! not used for selection
goto 11
elseif (string(2:12).eq.'atom_index=') then ! atom index line
atind = iquoted(string,nline) ! (atom Nr); check if we need it
if (mode.eq.1.and.(atind.eq.indref.or.indref.eq.0)) addos=.true.
goto 11
elseif (string(2:9).eq.'species=') then ! atom species line
chlab = squoted(string,nline) ! (atom label); check if we need it
! write (6,*) ' Compare chlab=',trim(chlab),
! . '-- and species=',trim(species),'---'
if (mode.eq.2.and.(trim(chlab).eq.trim(species))) addos=.true.
goto 11
elseif (string(2:3).eq.'Z=') then ! Z-value,
zvalue = iquoted(string,nline) ! not used for selection
goto 11
elseif (string(2:10).eq.'position=') then ! atom position,
! not used for selection
goto 11
elseif (string(2:3).eq.'n=') then ! n-value:
n = iquoted(string,nline)
if (nref.ne.0.and.n.ne.nref) addos=.false.
goto 11
elseif (string(2:3).eq.'l=') then ! l-value:
l = iquoted(string,nline)
if (lref.ne.-1.and.l.ne.lref) addos=.false.
goto 11
elseif (string(2:3).eq.'m=') then ! m-value:
m = iquoted(string,nline)
if (mref.ne.9.and.m.ne.mref) addos=.false.
goto 11
elseif (string(2:3).eq.'z=') then ! zeta-value: not used for selection
zeta = iquoted(string,nline) ! (all zetas are included into summation)
goto 11
elseif (string(2:3).eq.'P=') then ! whether polarisation orbital
polar = squoted(string,nline)
! write(io1,*)' polar=',polar
if (trim(polar).eq.'true') addos=.false. ! skip P="true"
! write(io1,*)' addos=',addos
goto 11
elseif (string(1:1).eq.'>') then ! closing the orbital block
! write (io1,*) n,l,m,zeta,polar,' - Orbital block closed'
! write (io1,*) ' addos=',addos
goto 11
elseif (string(1:6).eq.'<data>') then ! PDOS block follows:
if (addos) write (io1,303) atind,n,l,m,zeta
do it=1,nt
read (ii1,*) (dos1(ispin),ispin=1,nspin) ! read always...
if (addos) then
do ispin=1,nspin
dos(it,ispin)=dos(it,ispin)+dos1(ispin) ! add only if 'addos'
enddo
endif
enddo
read (ii1,'(a7)',err=201,end=202) dtail
if (dtail.ne.'</data>') then ! check that <data> block
write (6,307) nline,dtail,nt ! correctly ends with </data>
stop
endif
goto 11
elseif (string(1:10).eq.'</orbital>') then ! orbital ends:
goto 11
elseif (string(1:7).eq.'</pdos>') then ! regular end of everything
goto 19
elseif (ene_fo) then ! expect energy value...
nene=nene+1
read (string,*) ene(nene)
goto 11
else
write(6,*) ' Unknown identifier in PDOS file, line ',nline
write(6,*) string
stop
endif
19 continue
close (ii1)
C --- write down accumulated DOS values: -----------------
if (isferm) write (io1,310) efermi
write (io1,304)
do it=1,nt
write (io1,305) ene(it),(dos(it,ispin),ispin=1,nspin)
enddo
close (io1)
stop
201 write(6,*) 'Error reading PDOS file, line ',nline
stop
202 write(6,*) 'Unexpected end of PDOS file, line ',nline
stop
401 write (6,*) ' Error reading ',trim(string),' as numeric'
goto 100
402 write (6,*) ' Error reading ',trim(string),' as string'
goto 100
301 format('# summed up for atom index ',i5)
302 format('# summed up over atom species: ',(a))
303 format('# Add data for atom_index =',i4,', n,l,m,z=',4i3)
304 format('#',/,'# Energy',10x,'spin 1',8x,'spin2',/,'#')
305 format(f13.7,4f15.8)
306 format(' Illegal first character (ASCII =',i4,').',/
. ' Atom number must start from [0-9],',
. ' atom label - from [A-Z,a-z]. Try again...')
307 format(' Data mismatch in line ',i6,' :',/
. ' read in -|',a7,'|- instead of expected </data>',
. ' after ',i4,' data lines')
309 format('#',/,'# partial DOS extracted from file ',60a)
310 format('#',/,'# Fermi energy :',f13.8,' eV')
308 format('# selecting orbital keys n=',i2,', l=',i2,', m=',i2,/
. '# The summation is over all zetas (see z= below),',/
. '# neglecting polarisation orbitals (retain orbitals',
. ' with P="false" only)',/'#')
end
C
C ......................................................................
C
integer function iquoted(string,nline)
C
C returns integer value contained in string between two delimiters
C
implicit none
integer lquote,rquote,nline
character*80 string
character*1 delim
data delim /'"'/
lquote = scan(string,delim)
rquote = scan(string,delim,back=.true.)
if (lquote.lt.rquote) then
read (string(lquote+1:rquote-1),*) iquoted
else
write (6,*) ' Error locating quotes in line ',nline
write (6,*) ' lquote, rquote =',lquote, rquote
write (6,*) string
stop
endif
return
end
C
C ......................................................................
C
character*6 function squoted(string,nline)
C
C returns substring contained in string between two delimiters
C
implicit none
integer lquote,rquote,nline
character*80 string
character*1 delim
data delim /'"'/
lquote = scan(string,delim)
rquote = scan(string,delim,back=.true.)
if (lquote.lt.rquote) then
read (string(lquote+1:rquote-1),*) squoted
else
write (6,*) ' Error locating quotes in line ',nline
write (6,*) ' lquote, rquote =',lquote, rquote
write (6,*) string
stop
endif
return
end
C
C ......................................................................
C
integer function iparsed(string,llabel,rlabel)
C
C reads in an integer value from the field of string 'string'
C which is situated between two substrings 'llabel' and 'rlabel'.
implicit none
integer first,last
character*80 string
character*80 llabel,rlabel
C --- position in 'string' just after appearance of substring 'llabel'
first = scan(string,llabel(1:len_trim(llabel))) +
+ len_trim(llabel)
C --- position in 'string' just before the appearance of substring 'rlabel'
last = scan(string,rlabel(1:len_trim(rlabel)),back=.true.) -
- len_trim(rlabel)
if (first.le.last) then
read (string(first:last),*,err=801) iparsed
return
else
write (6,*) ' Fail to locate integer field in the string :'
write (6,*) string
write (6,*) ' first, last =',first,last
stop
endif
801 write (6,*) ' Fail to read integer number from positions ',
. first,' through ',last,' of the string :'
write (6,*) string
stop
end
C
C ......................................................................
C
double precision function dparsed(string,llabel,rlabel)
C
C reads in an integer value from the field of string 'string'
C which is situated between two substrings 'llabel' and 'rlabel'.
implicit none
integer first,last
character*80 string
character*80 llabel,rlabel
C --- position in 'string' just after appearance of substring 'llabel'
first = scan(string,llabel(1:len_trim(llabel))) +
+ len_trim(llabel)
C --- position in 'string' just before the appearance of substring 'rlabel'
last = scan(string,rlabel(1:len_trim(rlabel)),back=.true.) -
- len_trim(rlabel)
if (first.le.last) then
read (string(first:last),*,err=801) dparsed
return
else
write (6,*) ' Fail to locate dp field in the string :'
write (6,*) string
write (6,*) ' first, last =',first,last
stop
endif
801 write (6,*) ' Fail to read dp number from positions ',
. first,' through ',last,' of the string :'
write (6,*) string
stop
end
C
C ......................................................................
