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