DOS plot (1) Outline

Posted on : April 14, 2012 (Sat) 12:58:08

by Koji Kobashi

I recall that several years ago, Dr Kotani, who was then in Prof Akai's lab, presented a package to plot a DOS. One may still find it somewhere in the BBC. All you need is to input a following command in the terminal (example):

./gpd out/Cu_00_dos

To run this properly, one needs to prepare "gpdos" that has been compiled as follows:

g77 gpdos gpdos.f

Alternatively, one can also do the following:

g77 gpdos.f
mv a.out gpdos

Then an eps file of the DOS of "Cu_00_dos" is created that could be seen by a linux software "evince" (Document Viewer) or other softwares.

DOS plot (2) gpdos.f

Posted on : April 14, 2012 (Sat) 13:12:34

by Koji Kobashi

The original file of gpdos.f is as follows. Please note that Prof Akai owns the copyright. The program is "program dosplot" but the file name is gpdos.f. It consists of a Fortran 77 code with 159 lines.
-----------

program dosplt
c-----------------------------------------------------------------------
c construct DOS data from cpa98 output.
c coded by H.Akai, 22 Feb. 1998, Osaka
c-----------------------------------------------------------------------
implicit real*8 (a-h,o-z)
parameter(msex=201,icmpx=10,maxlin=10000,mmxl=4)
real*8 dat(mmxl*icmpx+3,msex)
character file*80,a*120,fmt*80,fmtx*80
logical exist
c--- get input file
50 write(*,'(a)')' input file name ?'
read(*,'(a)',end=10)file
call lftfil(file)
call chleng(file,len)
if(len .le. 0) go to 10
c
c If a directory of the same name exists, inquire can not
c ditect existence of a file correctly. Therefore one should
c use open statement with error return.
inquire(file=file,exist=exist)
if(.not. exist) then
write(*,'(a)')' file not found'
go to 50
endif
open(12,file=file)
c--- create output file
ip=index(file,'.')
if(ip .ge. 2) file=file(1:ip-1)
call lftfil(file)
call chleng(file,len)
if(len .le. 0) go to 10
c
c--- read-in files (unit=12) and find DOS data.
c
mse=0
mxl=0
icmp=0
do 70 line=1,maxlin
c --- get the number of meshes, l_max, and number of components
read(12,'(a)',end=10)a
if(a(5:28) .eq. 'meshr mse ng mxl') then
read (12,*)meshr,mse,ng,mxl
endif
if(a(4:8) .eq. 'ntyp='.and. a(13:17) .eq. 'natm='
& .and. a(22:27) .eq. 'ncmpx=') then
read(a(28:29),'(i2)')icmp
endif
c --- the last data before the DOS data begin are adopted
if(a(2:19) .eq. 'DOS of component 1') go to 100
70 continue
write(*,'(a)')' seems not to contain DOS data'
go to 10
100 write(*,'(a,3i3)')' mse,mxl,icmp=',mse,mxl,icmp
c --- check the sizes ---
if(mse .gt. msex) then
write(*,'(a,i3)')' ***error...mse>',msex
go to 10
endif
if(mxl .gt. mmxl) then
write(*,'(a,i2)')' ***error...mxl>',mmxl
go to 10
endif
if(icmp .gt. icmpx) then
write(*,'(a,i2)')' ***error...icmp>',icmpx
go to 10
endif
c --- Now the format statement can be constructed
write(fmtx,'(a,i1,a)')'(1x,f7.4,3x,',mxl,'f8.4)'
c write(*,*)'format=',fmtx
open(13,file=file(1:len)//'.plt',status='unknown')
c write(13,'(a)')'#!/usr/bin/gnuplot -persist'
write(13,'(a)')'#!/usr/bin/gnuplot'
write(13,'(2a)')'# set terminal postscript landscape noenhanced ',
& 'monochrome dashed defaultplex "Helvetica" 14'
write(13,'(a)')'# set output "dos.eps"'
write(13,'(a)')'set xzeroaxis lt -1 lw 0.1'
write(13,'(a)')'set yzeroaxis lt -1 lw 0.1'
write(13,'(a)')'set mxtics 5'
write(13,'(a)')'set mytics 5'
write(13,'(a)')'set xtics border mirror norotate 0.5'
write(13,'(a)')'set ytics border mirror norotate 10'
write(13,'(3a)')'set title "',file(1:len),
& '" 0.000000,0.000000 ""'
write(13,'(2a)')
& 'set xlabel "Energy relative to the Fermi energy (Ry)" ',
& '0.000000,0.000000 ""'
write(13,'(a)')
& 'set xrange [ * : * ] noreverse nowriteback'
write(13,'(a)')
& 'set ylabel "DOS (1/Ry/unit cell)" 0.000000,0.000000 ""'
write(13,'(a)')
& 'set yrange [ * : * ] noreverse nowriteback'
write(13,'(2a)')'plot "-" u 1:2 t "" w l, ',
& '"-" u 1:(-1*$2) t "" w l'
do 20 is=1,2
c write(*,'(1x,i3)')is
c write(*,'(1x,a)')file(1:len)//ext(is)
do 30 line=1,maxlin
if(a(2:19) .eq. 'DOS of component 1') then
c write(*,*)'DOS 1 met'
read(12,fmtx)((dat(l,k),l=3,mxl+3),k=1,mse)
c write(*,fmtx)((dat(l,k),l=3,mxl+3),k=1,mse)
60 continue
do 40 i=2,icmp
c write(*,*)'DOS',i,' met'
i0=mxl*(i-1)+4
read(12,'(//)')
40 read(12,fmtx)
& (dmy,(dat(l,k),l=i0,i0+mxl-1),k=1,mse)
endif
if(a(2:10) .eq. 'total DOS') then
c write(*,*)'total DOS met'
read(12,'(1x,f12.7,3x,f13.5)')(dat(1,k),dat(2,k),k=1,mse-1)
 dat(1,mse)=dat(1,mse-1)
 dat(2,mse)=dat(2,mse-1)
write(fmt,'(a,i2,a)')'((f8.4,f7.2,f8.4,',mxl*icmp,'f7.2))'
write(13,fmt)((dat(l,k),l=1,mxl*icmp+3),k=1,mse)
write(13,'(a)')'end'
go to 20
endif
30 read(12,'(a)',end=10)a
20 read(12,'(a)',end=10)a
write(13,'(a)')'pause -1'
write(13,'(a)')'# EOF'
close(13)
10 continue
end
subroutine lftfil(fil)
c---------------------------------------------------------------------
implicit real*8 (a-h,o-z)
character fil*(*)
n=len(fil)
j=0
is=1
do 30 i=1,n
is=i
if(llt(fil(i:i),'0') .or. llt('9',fil(i:i))) go to 40
30 continue
40 do 10 i=is,n
if(fil(i:i).eq.' ') go to 10
j=j+1
fil(j:j)=fil(i:i)
10 continue
if(j.ge.n) return
fil(j+1:n)=' '
return
end
subroutine chleng(a,ln)
c---------------------------------------------------------------------
character a*(*)
n=len(a)
do 10 i=n,1,-1
ln=i
if(a(i:i) .ne. ' '.and. ichar(a(i:i)) .gt. 27) return
10 continue
ln=0
return
end

DOS plot (3) gpdos.f for covalent crystals

Posted on : April 14, 2012 (Sat) 13:17:12

by Koji Kobashi

In case of crystals with covalent bonds only, one can use the following code for gpdos.f:
-------------------------

program dosplt
c-----------------------------------------------------------------------
c construct DOS data from cpa98 output.
c coded by H.Akai, 22 Feb. 1998, Osaka
c-----------------------------------------------------------------------
implicit real*8 (a-h,o-z)
parameter(msex=201,icmpx=10,maxlin=10000,mmxl=4)
real*8 dat(mmxl*icmpx+3,msex)
character file*80,a*120,fmt*80,fmtx*80
logical exist
c--- get input file
50 write(*,'(a)')' input file name ?'
read(*,'(a)',end=10)file
call lftfil(file)
call chleng(file,len)
if(len .le. 0) go to 10
c
c If a directory of the same name exists, inquire can not
c ditect existence of a file correctly. Therefore one should
c use open statement with error return.
inquire(file=file,exist=exist)
if(.not. exist) then
write(*,'(a)')' file not found'
go to 50
endif
open(12,file=file)
c--- create output file
ip=index(file,'.')
if(ip .ge. 2) file=file(1:ip-1)
call lftfil(file)
call chleng(file,len)
if(len .le. 0) go to 10
c
c--- read-in files (unit=12) and find DOS data.
c
mse=0
mxl=0
icmp=0
do 70 line=1,maxlin
c --- get the number of meshes, l_max, and number of components
read(12,'(a)',end=10)a
if(a(5:28) .eq. 'meshr mse ng mxl') then
read (12,*)meshr,mse,ng,mxl
endif
if(a(4:8) .eq. 'ntyp='.and. a(13:17) .eq. 'natm='
& .and. a(22:27) .eq. 'ncmpx=') then
read(a(28:29),'(i2)')icmp
endif
c --- the last data before the DOS data begin are adopted
if(a(2:19) .eq. 'DOS of component 1') go to 100
70 continue
write(*,'(a)')' seems not to contain DOS data'
go to 10
100 write(*,'(a,3i3)')' mse,mxl,icmp=',mse,mxl,icmp
c --- check the sizes ---
if(mse .gt. msex) then
write(*,'(a,i3)')' ***error...mse>',msex
go to 10
endif
if(mxl .gt. mmxl) then
write(*,'(a,i2)')' ***error...mxl>',mmxl
go to 10
endif
if(icmp .gt. icmpx) then
write(*,'(a,i2)')' ***error...icmp>',icmpx
go to 10
endif
c --- Now the format statement can be constructed
write(fmtx,'(a,i1,a)')'(1x,f7.4,3x,',mxl,'f10.4)'
c write(*,*)'format=',fmtx
open(13,file=file(1:len)//'.plt',status='unknown')
c write(13,'(a)')'#!/usr/bin/gnuplot -persist'
write(13,'(a)')'#!/usr/bin/gnuplot'
write(13,'(2a)')'set terminal postscript landscape noenhanced ',
& 'monochrome dashed defaultplex "Helvetica" 14'
c write(13,'(a)')'set output "dos.eps"'
write(13,'(a)')'set output "dos.eps"'
write(13,'(a)')'set border 15 lw 0.1'
write(13,'(a)')'set size ratio 0.7071'
write(13,'(a)')'set xzeroaxis lt -1 lw 0.1'
write(13,'(a)')'set yzeroaxis lt -1 lw 0.1'
write(13,'(a)')'set mxtics 2'
write(13,'(a)')'set mytics 2'
write(13,'(a)')'set xtics border mirror norotate 0.2'
write(13,'(a)')'set ytics border mirror norotate 10'
write(13,'(3a)')'set title "',file(1:len)
write(13,'(2a)')
& 'set xlabel "Energy relative to the Fermi energy (Ry)" '
write(13,'(a)')
& 'set xrange [ * : * ] noreverse nowriteback'
write(13,'(a)')
& 'set ylabel "DOS (1/Ry/unit cell)" '
write(13,'(a)')
& 'set yrange [ * : * ] noreverse nowriteback'
write(13,'(2a)')'plot "-" u 1:2 t "" w l 1'
do 20 is=1,2
c write(*,'(1x,i3)')is
c write(*,'(1x,a)')file(1:len)//ext(is)
do 30 line=1,maxlin
if(a(2:19) .eq. 'DOS of component 1') then
c write(*,*)'DOS 1 met'
read(12,fmtx)((dat(l,k),l=3,mxl+3),k=1,mse)
c write(*,fmtx)((dat(l,k),l=3,mxl+3),k=1,mse)
60 continue
do 40 i=2,icmp
c write(*,*)'DOS',i,' met'
i0=mxl*(i-1)+4
read(12,'(//)')
40 read(12,fmtx)
& (dmy,(dat(l,k),l=i0,i0+mxl-1),k=1,mse)
endif
if(a(2:10) .eq. 'total DOS') then
c write(*,*)'total DOS met'
read(12,'(1x,f12.7,3x,f13.5)')(dat(1,k),dat(2,k),k=1,mse-1)
dat(1,mse)=dat(1,mse-1)
dat(2,mse)=dat(2,mse-1)
write(fmt,'(a,i2,a)')'((f8.4,f7.2,f8.4,',mxl*icmp,'f7.2))'
write(13,fmt)((dat(l,k),l=1,mxl*icmp+3),k=1,mse)
write(13,'(a)')'end'
go to 20
endif
30 read(12,'(a)',end=10)a
20 read(12,'(a)',end=10)a
write(13,'(a)')'pause -1'
write(13,'(a)')'# EOF'
close(13)
10 continue
end
subroutine lftfil(fil)
c---------------------------------------------------------------------
implicit real*8 (a-h,o-z)
character fil*(*)
n=len(fil)
j=0
is=1
do 30 i=1,n
is=i
if(llt(fil(i:i),'0') .or. llt('9',fil(i:i))) go to 40
30 continue
40 do 10 i=is,n
if(fil(i:i).eq.' ') go to 10
j=j+1
fil(j:j)=fil(i:i)
10 continue
if(j.ge.n) return
fil(j+1:n)=' '
return
end
subroutine chleng(a,ln)
c---------------------------------------------------------------------
character a*(*)
n=len(a)
do 10 i=n,1,-1
ln=i
if(a(i:i) .ne. ' '.and. ichar(a(i:i)) .gt. 27) return
10 continue
ln=0
return
end

DOS plot (4) gpd (bash)

Posted on : April 14, 2012 (Sat) 13:20:38

by Koji Kobashi

The gpd file of bash reads as follows:

#!/bin/bash
echo $1|./gpdos
echo $1|awk '{z=split($1,name,".");file=name[1] ".plt";print "#!/bin/bash";print "chmod a+x " file;print "./" file}'>gpd.tmp
chmod a+x gpd.tmp
./gpd.tmp
rm gpd.tmp

DOS plot (5) gpd (csh)

Posted on : April 14, 2012 (Sat) 13:24:36

by Koji Kobashi

For the case of csh, use the following for gpd:

#!/bin/csh
echo $1|./gpdos
echo $1|awk '{z=split($1,name,".");file=name[1] ".plt";print "#\ !/bin/csh";print "chmod a+x " file;print "./" file}'>gpd.tmp
chmod a+x gpd.tmp
./gpd.tmp
rm gpd.tmp

DOS plot (6) DOS plot on Microsoft Excel

Posted on : April 14, 2012 (Sat) 13:35:57

by Koji Kobashi

Using Microsoft Excel to make a DOS plot is an alternative way. It is tedious but once one makes a macro, it is as convenient as the way using ./gpd XXX because one can change the energy units from Ry to eV, and also modify the plot as you like.

In the "out" folder, one can find a dos file in which there are lines as follows:

total DOS
-1.0447500 7.97723
-1.0342500 7.08628
-1.0237500 8.62280
-1.0132500 11.46989
-1.0027500 10.33481
-0.9922500 10.38610
-0.9817500 20.78855
-0.9712500 14.75674
... .....

You can copy the relevant data to an Excel sheet, and create a DOS graph. If you make and use a macro, it is quite convenient.

[Re:06] DOS plot (7) Current status

Posted on : July 22, 2014 (Tue) 13:06:38

by Koji Kobashi

In the current version of cpa2002v009c.tar.gz, Professor Akai added gpd.f in the source file so that all you need is to run (Fortran compiler) god in/god.f. Else (Fortran compiler) gpd.f creates a.out and you change the name to gpd. ./gpd out/{file name} generates the DOS plot. The detail of the procedure also is described in the release note of cpa2002v009c.tar.gz.
In the default code, the DOS plot includes spin-up and spin-down states so that for nonmagnetic materials, it may be modified.
Otherwise, copy and paste the data on Excel sheet, and make a DOS graph.