C C grdint, a script to integrate 3-dim grid function C (i.e. LDOS, RHO, DRHO, etc.) written in SIESTA by subr. iorho C over a region of choice (a sphere centered at a given atom). C C !!! -------------------- IMPORTANT ---------------------- !!! C compile this code with the same compiler switches as Siesta, C in what regards using single/double precision, C otherwise reading the data from unformatted files C written by iorho.F might be spoiled. C Don't say you haven't been warned. C !!! C C Written by Andrei Postnikov, Mar 2007 C postnikov@univ-metz.fr C program grdint implicit none integer ii1,ii2,io1 parameter (ii1=11,ii2=12,io1=14) integer mesh(3),nspin,is,ii,jj,iat,nat,nz,iix,iiy,iiz, . ind,mn,ialloc,ishift,jshift,kshift,limit,nrat parameter (limit=1) ! tried translations along each lattice vector character inpfil*60,outfil*60,syslab*30,suffix*6,owrite*1 logical filexist double precision b2ang,cc_bohr(3,3),diff,coort(3),cell(3,3), . dummy,modu,rmesh(3),small,volume,srad,srad2, . dist2,charge,spin parameter (b2ang=0.529177) ! Bohr to Angstroem parameter (small=1.0d-5) integer, allocatable :: ityp(:),iz(:) double precision, allocatable :: coor_bohr(:,:) real, allocatable :: func(:,:,:,:) ! NB! single precision, as in iorho.F real sum(8) ! NB! single precision, as in iorho.F C C string manipulation functions in Fortran used below: C len_trim(string): returns the length of string C without trailing blank characters, C trim(string) : returns the string with railing blanks removed C char(integer) : returns the character in the specified position C of computer's ASCII table, i.e. char(49)=1 write (6,701,advance="no") 701 format(" Specify SystemLabel (or 'siesta' if none): ") read (5,*) syslab C -- handle XV file: ---------------------------- inpfil = trim(syslab)//'.XV' open (ii1,file=inpfil,form='formatted',status='old',err=801) C -- check number of atoms in the XV file, for allocation: do ii=1,3 read (ii1,'(3x,3f18.9)',end=801,err=801) (dummy,jj=1,3) enddo read (ii1,*,end=802,err=802) nat allocate (ityp(nat)) allocate (iz(nat)) allocate (coor_bohr(1:3,1:nat),STAT=ialloc) if (ialloc.ne.0) then write (6,*) ' Fails to allocate space for ',nat,' atoms.' stop endif rewind (ii1) C -- read in translation vectors: do ii=1,3 read (ii1,'(3x,3f18.9)') (cc_bohr(jj,ii),jj=1,3) enddo read (ii1,*) nat do iat=1,nat read (ii1,'(i3,i6,3f18.9)',end=803,err=803) . ityp(iat),iz(iat),coor_bohr(:,iat) enddo close (ii1) C C --- Look for grid data files to include: 103 write (6,706,advance="no") read (5,*) suffix inpfil = trim(syslab)//'.'//trim(suffix) open (ii2,file=inpfil,form='unformatted',status='old',err=806) write (6,*) 'Found and opened: ',trim(inpfil) read (ii2,err=807) cell C check if cell vectors match those from .XV, just for the case... diff = 0.0 do ii=1,3 do jj=1,3 diff = diff + (cell(ii,jj)-cc_bohr(ii,jj))**2 enddo enddo if (diff.gt.small) then write (6,*) ' cell vectors from XV: ' write (6,'(3f12.6)') (cc_bohr(ii,1),ii=1,3) write (6,'(3f12.6)') (cc_bohr(ii,2),ii=1,3) write (6,'(3f12.6)') (cc_bohr(ii,3),ii=1,3) write (6,*) ' and from grid file: ' write (6,'(3f12.6)') (cell(ii,1),ii=1,3) write (6,'(3f12.6)') (cell(ii,2),ii=1,3) write (6,'(3f12.6)') (cell(ii,3),ii=1,3) write (6,*) ' differ!' stop endif read (ii2,err=808) mesh, nspin allocate ( . func(1:mesh(3),1:mesh(2),1:mesh(1),1:nspin), . STAT=ialloc) if (ialloc.ne.0) then write (6,*) ' Fails to allocate space for ', . mesh(1)*mesh(2)*mesh(3),' grid points,', . ' nspin =',nspin stop endif write (6,*) 'mesh = (',mesh,'), nspin=',nspin C microvolume per one mesh point: volume = ( cell(1,1)*cell(2,2)*cell(3,3) + + cell(1,2)*cell(2,3)*cell(3,1) + + cell(1,3)*cell(2,1)*cell(3,2) - - cell(1,1)*cell(2,3)*cell(3,2) - - cell(1,2)*cell(2,1)*cell(3,3) - - cell(1,3)*cell(2,2)*cell(3,1) ) / / float( mesh(1)*mesh(2)*mesh(3) ) do is=1,nspin do iiz=1,mesh(3) do iiy=1,mesh(2) read (ii2,err=809,end=810) . (func(iiz,iiy,iix,is),iix=1,mesh(1)) enddo enddo enddo close (ii2) C --- open output file: outfil = trim(syslab)//'.SPH' inquire (file=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=outfil,form='formatted',status='REPLACE') else write (6,*) ' Then rename is first. Bye...' stop endif else open (io1,file=outfil,form='formatted',status='NEW') endif write(io1,201) C -- loop over atoms; for each one integrate and write onto io1 101 write (6,*) ' Integrate around atom No.: (0 to quit)' read (5,*,err=101) nrat if (nrat.eq.0) then deallocate (func) close (io1) stop ! regular termination elseif (nrat.lt.0.or.nrat.gt.nat) then goto 101 endif 102 write (6,*) ' Over sphere of radius (in Bohr) :' read (5,*,err=102) srad if (srad.lt.0.0) goto 102 srad2 = srad*srad sum(:) = 0.d0 ! accumulates spin-resolved charge density over sphere do ishift=-limit,limit do jshift=-limit,limit do kshift=-limit,limit do iiz = 1,mesh(3) ! or, 0 to mesh(1)-1 - check ! do iiy = 1,mesh(2) do iix = 1,mesh(1) do jj=1,3 coort(jj) = (float(iix)/float(mesh(1))+ishift)*cell(jj,1) + + (float(iiy)/float(mesh(2))+jshift)*cell(jj,2) + + (float(iiz)/float(mesh(3))+kshift)*cell(jj,3) enddo dist2 = (coort(1) - coor_bohr(1,nrat))**2 + + (coort(2) - coor_bohr(2,nrat))**2 + + (coort(3) - coor_bohr(3,nrat))**2 if (dist2.lt.srad2) then C write(io1,203)iix,iiy,iiz,coort,dist2, C . (func(iiz,iiy,iix,is),is=1,2) do is = 1,nspin sum(is) = sum(is) + func(iiz,iiy,iix,is)*volume enddo endif enddo ! do iix enddo enddo enddo enddo enddo C --- print integrated values: charge = sum(1) if (nspin.eq.1) then charge = sum(1) spin = 0.d0 else charge = sum(1) + sum(2) spin = sum(1) - sum(2) endif write(io1,202) nrat,(coor_bohr(jj,nrat),jj=1,3), . srad,charge,spin goto 101 201 format(' Atom',4x,11('-'),' at ',11('-'),4x,'Raduis', . 5x,'Charge',4x,'Spin') 202 format(i4,3x,'(',3f9.4,' )',f8.4,2x,2f9.4) 203 format(3i4,3f8.4,' dist2=',f12.4,' func=',1p2e12.5) 706 format (' Add grid property (LDOS, RHO, ...;', . ' or BYE if none): ') 801 write (6,*) ' End/Error reading XV for cell vector ',ii stop 802 write (6,*) ' End/Error reading XV for number of atoms line' stop 803 write (6,*) ' End/Error reading XV for atom number ',iat stop 806 write (6,*) ' A wild guess! There is no file ', . trim(inpfil),'; close XSF and quit.' close (io1) stop 807 write (6,*) ' Error reading cell vectors' stop 808 write (6,*) ' Error reading n1,n2,n3,nspin ' stop 809 write (6,*) ' Error reading function values on the grid' stop 810 write (6,*) ' Unexpected end of data for iiy=',iiy, . ' iiz=',iiz,' is=',is stop 811 write (6,*) ' Error opening file ', . trim(outfil),' as new unformatted' stop end