On Sun, 26 Jun 2005, xianghjun wrote:
| Dear Andrei Postnikov,
| I have tried rho2xsf, and it works well.
| Now I want to plot the Fermi surface using Xcrysden.
| I learn that you have a program named ene2bxsf to do this job.
| Would you please send a copy of ene2bxsf to me?
| Thanks a lot.
|
| Best regards,
| xianghjun
Dear Xianghjun,
please find in the attachment a code to calculate Fermi surface file
in Xcrysden. I did not check it extensively, but it seems to work
also with the version 1.4 (that with an additional option
to plot several sheets of the Fermi surface in the same plot).
Please let me know if you encounter any problems.
Be sure that you use many enough k-points without shifting.
Good luck,
Andrei
C
C ene2bxsf, a script to make the band-XSF file (for plotting
C the Fermi surface) from siesta KP and EIG files.
C The k-mesh must be generated without shift
C (requirement of the BXSF format).
C
C Written by Andrei Postnikov, Mar 2005
C [EMAIL PROTECTED]
C
program ene2bxsf
implicit none
integer MKP,MNB
parameter (MKP=8000,MNB=500)
integer ii1,ii2,ii3,io1
parameter (ii1=11,ii2=12,ii3=13,io1=14)
integer ii,jj,nbands,nbmin,nbmax,nkp,ikp,ndum,iis,
. mdiv1,mdiv2,mdiv3,ndiv(3,MKP),ind,irrek(MKP),nspin,is,
. iband,iik,idiv1,idiv2,idiv3,itry1,itry2,itry3,ib,
. homo(2),lumo(2),ik,in2,in3
character inpfil*60,outfil*60,syslab*30,suffix*6
real*8 cell(3,3),efermi,twopi,rcell(3,3),small,relmin(3),
. relk(3,MKP),kkp(3),eneb(MNB,2),enek(MKP),dum
equivalence (eneb,enek)
parameter (small=1.0d-04)
C
C string manipulation functions in Fortran used below:
C len_trim(string): returns the length of string
C without trailing blank characters,
C char(integer) : returns the character in the specified position
C of computer's ASCII table, i.e. char(49)=1
C
twopi = 8.d0*atan(1.d0)
write (6,701)
701 format(" Specify SystemLabel (or 'siesta' if none): ",$)
read (5,*) syslab
C --- open .XV and .EIG :
inpfil = syslab(1:len_trim(syslab))//'.XV'
open (ii1,file=inpfil,form='formatted',status='old',err=801)
write (6,*) 'Found and opened: ',inpfil
inpfil = syslab(1:len_trim(syslab))//'.EIG'
open (ii2,file=inpfil,form='formatted',status='old',err=801)
write (6,*) 'Found and opened: ',inpfil
C --- read Fermi energy and total number of bands from .EIG :
read (ii2,*) efermi
read (ii2,*) nbands, nspin, nkp
if (nbands.gt.MNB) then
write (6,*) ' nbands=',nbands,' .gt. MNB=',MNB
stop
endif
if (nspin.ne.1.and.nspin.ne.2) then
write (6,*) 'A problem encountered: nspin=',nspin
stop
endif
if (nkp.gt.MKP) then
write (6,*) ' nkp=',nkp,' .gt. MKP=',MKP
stop
endif
C --- finds bands crossing the Fermi energy:
do is=1,nspin
homo(is)=0 ! higest (partially?) occupied band
lumo(is)=nbands+1 ! lowest (partially?) unoccupied band
enddo
write (6,*) ' nkp=',nkp,' nbands=',nbands
do ik=1,nkp
read (ii2,"(i5,10f12.5,/,(5x,10f12.5))") iik, ! written in
. ((eneb(ib,is),ib=1,nbands),is=1,min(nspin,2)) ! ioeif.f
if (iik.ne.ik) then
write (6,*) ' iik=',iik,'.ne. ik=',ik,' for spin ',is
stop
endif
do is=1,nspin
do ib=1,nbands
if (eneb(ib,is).lt.efermi.and.ib.gt.homo(is)) homo(is) = ib
if (eneb(ib,is).gt.efermi.and.ib.lt.lumo(is)) lumo(is) = ib
enddo
enddo
enddo
do is=1,nspin
write (6,*) ' is=',is,' homo, lumo=',homo(is),lumo(is)
if (homo(is).lt.lumo(is)) write (6,201) is, homo(is),lumo(is)
enddo
C --- read cell vectors from .XV, convert to Ang, find reciprocal:
do ii=1,3
read (ii1,*,end=803,err=803) (cell(ii,jj),jj=1,3)
enddo
close (ii1)
call inver3(cell,rcell)
C write (6,*) ' cell:'
C do ii=1,3
C write(6,'(3f15.6)') (cell(ii,jj),jj=1,3)
C enddo
C write (6,*) ' rcell:'
C do ii=1,3
C write(6,'(3f15.6)') (rcell(ii,jj),jj=1,3)
C enddo
C --- open .KP as old:
inpfil = syslab(1:len_trim(syslab))//'.KP'
open (ii3,file=inpfil,form='formatted',status='old',err=801)
write (6,*) 'Found and opened: ',inpfil
C --- read k-points from the .KP file and recover their fractional
C coordinates in terms of reciprocal vectors:
read (ii3,*) nkp
relmin(:)=1.d0/small
do ikp=1,nkp
read (ii3,*) iik,(kkp(jj),jj=1,3)
if (iik.ne.ikp) then
write (6,*) ' a mess in KP list: read in iik=',iik,
. ', but expected ikp=',ikp
stop
endif
C --- find relative coordinates of k-point:
do ii=1,3
relk(ii,ikp)=( cell(ii,1)*kkp(1) +
+ cell(ii,2)*kkp(2) +
+ cell(ii,3)*kkp(3) )/twopi
if (abs(relk(ii,ikp)).lt.relmin(ii).and.
. abs(relk(ii,ikp)).gt.small)
. relmin(ii)=abs(relk(ii,ikp))
enddo
C write(6,'(i5,3f15.6)') ikp,(relk(jj,ikp),jj=1,3)
enddo
close (ii3)
C write (6,*) ' relmin=',relmin
C --- relmin(ii) is now the smallest relative coordinate of k points
C along the reciprocal vector (ii). Good chance that its inverse
C is number of divisions.
mdiv1=1.d0/relmin(1)+small
mdiv2=1.d0/relmin(2)+small
mdiv3=1.d0/relmin(3)+small
write (6,*) ' mdiv1,2,3=',mdiv1,mdiv2,mdiv3
if ((mdiv1+1)*(mdiv2+1)*(mdiv3+1).gt.MKP) then
write (6,*) ' MKP needs to be at least ',
. (mdiv1+1)*(mdiv2+1)*(mdiv3+1)
stop
endif
C --- decifer all k-point coordinates as ndiv(ii,ikp)/mdiv_ii
do ikp=1,nkp
do ii=1,3
if (relk(ii,ikp).ge.0.d0) then
ndiv(ii,ikp)=relk(ii,ikp)/relmin(ii)+small
else
ndiv(ii,ikp)=(relk(ii,ikp)+1.d0)/relmin(ii)+small
endif
enddo
C write(6,203) ikp, ndiv(1,ikp),mdiv1,
C . ndiv(2,ikp),mdiv2,
C . ndiv(3,ikp),mdiv3
203 format(' ikp=',i4,': (',3(i3,'/',i3,3x),' )')
enddo
C --- attribute irreducible k-points to k-points on the grid:
ind=0
do idiv3 = 0,mdiv3
do idiv2 = 0,mdiv2
do idiv1 = 0,mdiv1
ind = ind + 1 ! global address on the mesh
itry1=mod(idiv1,mdiv1)
itry2=mod(idiv2,mdiv2)
itry3=mod(idiv3,mdiv3)
C write (6,302) idiv1,idiv2,idiv3,itry1,itry2,itry3
302 format (' Search idiv1,2,3=',3i5,' itry1,2,3=',3i5)
do 12 ikp=1,nkp
if (ndiv(1,ikp).ne.itry1) goto 12
if (ndiv(2,ikp).ne.itry2) goto 12
if (ndiv(3,ikp).ne.itry3) goto 12
irrek(ind) = ikp
C write (6,*) ' Found: ikp=',ikp
goto 14
12 enddo
C --- haven't find anything; try inversion:
itry1=mod(mdiv1-idiv1,mdiv1)
itry2=mod(mdiv2-idiv2,mdiv2)
itry3=mod(mdiv3-idiv3,mdiv3)
C write (6,303) itry1,itry2,itry3
303 format (' --- try inversed: itry1,2,3=',3i5)
do 13 ikp=1,nkp
C write (6,*) ' ikp:',ikp,' ndiv=',(ndiv(ii,ikp),ii=1,3)
if (ndiv(1,ikp).ne.itry1) goto 13
if (ndiv(2,ikp).ne.itry2) goto 13
if (ndiv(3,ikp).ne.itry3) goto 13
irrek(ind) = ikp
C write (6,*) ' Found: ikp=',ikp
goto 14
13 enddo
C --- haven't try anything with inversion as well;
write (6,*) ' No irreducible k-point found ',
. ' for k-grid point (',idiv1,'/',mdiv1,', ',
. idiv2,'/',mdiv2,', ',idiv3,'/',mdiv3,' ) .'
write (6,*) ' Are you sure your k-grid was without shift?'
stop
14 continue
enddo
enddo
enddo
C --- The writing sequence (into two BXSF files if two spins)
do is=1,nspin
if (nspin.eq.1) then
outfil = syslab(1:len_trim(syslab))//'.BXSF'
else
outfil =
. syslab(1:len_trim(syslab))//'_spin_'//char(48+is)//'.BXSF'
endif
open (io1,file=outfil,form='formatted',status='new',err=802)
write (6,*) 'Opened as new: ',outfil
write (io1, "(a10)") 'BEGIN_INFO'
write (io1, "(a5)") ' # '
write (io1, "(a33)") ' # Band-XCRYSDEN-Structure-File'
write (io1, "(a5)") ' # '
write (io1, "(a16,f10.4)") ' Fermi energy: ',efermi
write (io1, "(a8,/)") 'END_INFO'
write (io1, "(a23)") 'BEGIN_BLOCK_BANDGRID_3D'
write (io1, *) ' ',syslab(1:len_trim(syslab))
write (io1, "(a19,a1)") ' BANDGRID_3D_spin_',char(48+is)
write (io1, "(i5)") homo(is)-lumo(is)+1 ! No. of bands
write (io1, "(3i5)") mdiv1+1, mdiv2+1, mdiv3+1
write (io1, "(3f16.8)") 0.0, 0.0, 0.0
write (io1, "(3f16.8)") (rcell(jj,1)*twopi,jj=1,3)
write (io1, "(3f16.8)") (rcell(jj,2)*twopi,jj=1,3)
write (io1, "(3f16.8)") (rcell(jj,3)*twopi,jj=1,3)
do iband=lumo(is),homo(is)
write (io1, '(a7,i5)') ' BAND:',iband
C --- again read band energies of the needed band and spin from .ENE
C and write them in the correct order into .BXSF
rewind (ii2)
read (ii2,*) dum
read (ii2,*) ndum, ndum, ndum
C --- read in all energy values over all k points for the given spin band:
do ik=1,nkp
read (ii2,"(i5,10f12.5,/,(5x,10f12.5))") iik,
. ((dum,ib=1,nbands),iis=1,is-1), ! dummy read prev. spin, if any
. (dum,ib=1,iband-1),enek(ik),(dum,ib=iband+1,nbands),
. ((dum,ib=1,nbands),iis=is+1,nspin) ! dummy read next spin, if any
enddo
C --- write into .BXSF file:
write (io1, "(7f11.5)")
. (enek(irrek(ii)),ii=1,(mdiv1+1)*(mdiv2+1)*(mdiv3+1))
enddo ! do iband=lumo(is),homo(is)
write (io1, '(a17)') ' END_BANDGRID_3D'
write (io1, '(a21)') 'END_BLOCK_BANDGRID_3D'
close (io1)
enddo ! do is=1,ispin
close (ii2)
stop
201 format (' spin',i2,': band gap between bands ',i5,' and ',i5)
204 format (3f12.7)
205 format (1p,6e13.6)
C 205 format (1p,8e10.3)
206 format (' For is=',i1,': ',a3,'. grid value =',1p,e12.5,
. ' at ix,iy,iz=',3i4)
801 write (6,*) ' Error opening file ',
. inpfil(1:len_trim(inpfil)),' as old formatted'
stop
802 write (6,*) ' Error opening file ',
. outfil(1:len_trim(outfil)),' as new formatted'
stop
803 write (6,*) ' End/Error reading XV for cell vector ',ii
stop
end
C -----------------------------------
C
subroutine inver3(a,b)
C
C Inverts a 3x3 matrix
implicit none
integer i
double precision a(3,3), b(3,3), c
b(1,1) = a(2,2)*a(3,3) - a(3,2)*a(2,3)
b(1,2) = a(3,2)*a(1,3) - a(1,2)*a(3,3)
b(1,3) = a(1,2)*a(2,3) - a(2,2)*a(1,3)
b(2,1) = a(2,3)*a(3,1) - a(3,3)*a(2,1)
b(2,2) = a(3,3)*a(1,1) - a(1,3)*a(3,1)
b(2,3) = a(1,3)*a(2,1) - a(2,3)*a(1,1)
b(3,1) = a(2,1)*a(3,2) - a(3,1)*a(2,2)
b(3,2) = a(3,1)*a(1,2) - a(1,1)*a(3,2)
b(3,3) = a(1,1)*a(2,2) - a(2,1)*a(1,2)
do i = 1, 3
c=1.d0/(a(1,i)*b(i,1) + a(2,i)*b(i,2) + a(3,i)*b(i,3) )
b(i,1)=b(i,1)*c
b(i,2)=b(i,2)*c
b(i,3)=b(i,3)*c
enddo
end
C
C ...........................................................
