C C rho2xy, a script to scan the 3-dim grid function C (i.e. LDOS, RHO, DRHO, etc.) written in SIESTA by subr. iorho C along a straight line, using linear 4-point interpolation 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, Dec 2009 C postnikov@univ-metz.fr C program rho2xy implicit none integer ii1,ii2,io1 parameter (ii1=11,ii2=12,io1=14) integer mesh0(3),mesh1(3),ip,nspin,is,ii,jj,n1,n2,n3, . iat,nat,nz,iix,iiy,iiz,ind,mn,ibox,nbox, . ixmax,iymax,izmax,ixmin,iymin,izmin,ialloc, . ishift,jshift,kshift,i1,i2,i3 integer limit,maxa,maxo,maxuo,maxnh,maxna,il,ia,nrela(3),idim parameter (limit=5) ! tried translations along each lattice vector character inpfil*60,outfil*60,syslab*30,suffix*6, . unitlab*1,labunit*9,labbox*1,owrite*1 logical unitb,charge,waves,filexist double precision b2ang,cc_bohr(3,3),cc_ang(3,3),cc_inv(3,3), . coort(3),obeg(3),oend(3),leng,xx, . cell(3,3),dum,rmaxo,rela,modu,rmesh(3),drela(3) parameter (b2ang=0.529177) ! Bohr to Angstroem integer, allocatable :: ityp(:),iz(:) double precision, allocatable :: mass(:),coor_ang(:,:) character(len=2), allocatable :: label(:) real, allocatable :: func(:) ! NB! single precision, as in iorho.F real fintp,fmax,fmin ! NB! single precision, as in iorho.F external test_xv,read_xv,inver3,intpl04 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 write (6,701,advance="no") 701 format(" Specify SystemLabel (or 'siesta' if none): ") read (5,*) syslab inpfil = syslab(1:len_trim(syslab))//'.XV' open (ii1,file=inpfil,form='formatted',status='old',err=801) call test_xv(ii1,nat) allocate (ityp(nat)) allocate (iz(nat)) allocate (mass(nat)) allocate (label(nat)) allocate (coor_ang(1:3,1:nat),STAT=ialloc) if (ialloc.ne.0) then write (6,*) ' Fails to allocate space for ',nat,' atoms.' stop endif call read_xv(ii1,nat,ityp,iz,cc_ang,mass,label,coor_ang) call inver3(cc_ang,cc_inv) close (ii1) C --- set up and fill output box: call makelin(obeg,oend,leng) C --- open output file: outfil = syslab(1:len_trim(syslab))//'.RHO_XY' inquire (file=outfil, exist=filexist) if (filexist) then write (6,*) ' File ',outfil(1:len_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 C write (6,704) 102 write (6,705,advance="no") read (5,*) n1 if (n1.le.1) then write (6,*) ' Number must be greater than 1, try again.' goto 102 endif idim = 1 ! dimensionalty of output grid C --- Look for grid data files to include: 103 write (6,706,advance="no") read (5,*) suffix inpfil = syslab(1:len_trim(syslab))// . '.'//suffix(1:len_trim(suffix)) open (ii2,file=inpfil,form='unformatted',status='old',err=806) write (6,*) 'Found and opened: ',inpfil(1:len_trim(inpfil)) read (ii2,err=807) cell read (ii2,err=808) mesh0, nspin allocate (func(1:mesh0(1)*mesh0(2)*mesh0(3)),STAT=ialloc) if (ialloc.ne.0) then write (6,*) ' Fails to allocate space for ', . mesh0(1)*mesh0(2)*mesh0(3),' grid points.' stop endif write (6,*) 'mesh0 = (',mesh0,'), nspin=',nspin do is=1,nspin write (io1,707) inpfil,is,n1,obeg,oend,suffix ind=0 do iiz=1,mesh0(3) do iiy=1,mesh0(2) read (ii2,err=809,end=810) (func(ind+iix),iix=1,mesh0(1)) ind = ind + mesh0(1) enddo enddo C --- loop over mesh points do i1=1,n1 do ii=1,3 rmesh(ii) = obeg(ii) + (oend(ii)-obeg(ii))*(i1-1)/(n1-1) enddo xx = leng*(i1-1)/(n1-1) C the current X value along the path, starting from zero. C - rmesh(1..3) are absolute Cartesian coordinates C of the mesh point (i1) in Angstroem C Find its relative coordinates on the unit cell grid: do ii=1,3 rela = 0.0 do jj=1,3 rela = rela + cc_inv(ii,jj)*rmesh(jj) enddo C take modulo to make sure that it falls within [0,1]: modu = modulo( rela*mesh0(ii), dble(mesh0(ii)) ) nrela(ii) = floor(modu) + 1 drela(ii) = modu - nrela(ii) + 1 enddo C - mesh point rmesh(1..3) falls within the grid microcell C originating at the grid point nrela(1..3), C its relative coordinates within this microcell are drela(1..3) C Select neighboring grid points and make the interpolation: call intpl04 (func(1),fintp, . mesh0(1),mesh0(2),mesh0(3),nrela,drela) write (io1,'(2f18.8)') xx, fintp enddo ! do i1 = enddo ! do is=1,nspin deallocate (func) close (ii2) goto 103 201 format (3f20.8) 202 format (i4,3f20.8) 203 format (3i6) 204 format (3f12.7) 205 format (1p,6e13.6) 704 format (" Now define the grid. ") 705 format (" Enter number of grid points along the path: ") 706 format (' Add grid property (LDOS, RHO, ...;', . ' or BYE if none): ') 707 format ('# trace grid file ',a60,/'# spin ',i2,/ . '#',i6,' points from (',3f10.6,' )',/ . '#',16x,'to (',3f10.6,' ), endpoints included'/ . '#',8x,'dist (Ang)',10x,a6) 801 write (6,*) ' Error opening file ', . inpfil(1:len_trim(inpfil)),' as old formatted' stop 806 write (6,*) ' A wild guess! There is no file ', . inpfil(1:len_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 ', . outfil(1:len_trim(outfil)),' as new unformatted' stop end C C............................................................... C subroutine test_xv(ii1,nat) C C reads from ii1 (XV file) and returns the number of atoms C implicit none integer ii1,ii,jj,nat double precision dummy rewind ii1 do ii=1,3 read (ii1,'(3x,3f18.9)',end=801,err=801) (dummy,jj=1,3) enddo read (ii1,'(i12)',end=802,err=802) nat return 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 end C C............................................................... C subroutine read_xv(ii1,nat,ityp,iz,cc_ang,mass,label,coor_ang) C C reads again from ii1 (XV file) to the end, C returns cell vectors and coordinates C transformed from Bohr (as in XV) into Angstroem (as for XCrysDen) C implicit none integer ii1,nat,nat1,iat,ii,jj,ityp(nat),iz(nat) double precision tau(3,3),mass(nat),coor_ang(3,nat),amass(110), . b2ang,cc_bohr(3,3),cc_ang(3,3),coor_bohr(3) parameter (b2ang=0.529177) ! Bohr to Angstroem character*2 label(nat),alab(110) C data amass(1:40) / 1.00, 4.00, 6.94, 9.01, 10.81, . 12.01, 14.01, 16.00, 19.00, 20.18, ! Ne . 23.99, 24.31, 26.98, 28.09, 30.97, . 32.07, 35.45, 39.95, 39.10, 40.08, ! Ca . 44.96, 47.88, 50.94, 52.00, 54.94, . 55.85, 58.93, 58.69, 63.35, 65.39, ! Zn . 69.72, 72.61, 74.92, 78.96, 79.80, . 83.80, 85.47, 87.62, 88.91, 91.22 / ! Zr data amass(41:110) / . 92.91, 95.94, 97.91, 101.07, 102.91, . 106.42, 107.87, 112.41, 114.82, 118.71, ! Sn . 121.76, 127.60, 126.90, 131.29, 132.91, . 137.33, 138.91, 140.12, 140.91, 144.24, ! Nd . 146.92, 150.36, 151.96, 157.35, 158.93, . 162.50, 164.93, 167.26, 168.93, 173.04, ! Yb . 174.97, 178.49, 180.95, 183.84, 186.21, . 190.23, 192.22, 195.08, 196.97, 200.59, ! Hg . 204.38, 207.20, 208.98, 208.98, 209.99, . 222.02, 223.02, 226.03, 227.03, 232.04, ! Th . 231.04, 238.03, 237.05, 244.06, 243.06, . 247.07, 247.07, 251.08, 252.08, 257.10, ! Fm . 258.10, 259.10, 262.11, 261.11, 262.11, . 263.12, 262.12, 265.13, 266.13, 271.00 / ! Ds data alab / ' H','He','Li','Be',' B',' C',' N',' O',' F','Ne', . 'Na','Mg','Al','Si',' P',' S','Cl','Ar',' K','Ca', . 'Sc','Ti',' V','Cr','Mn','Fe','Co','Ni','Cu','Zn', . 'Ga','Ge','As','Se','Br','Kr','Rb','Sr',' Y','Zr', . 'Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd','In','Sn', . 'Sb','Te',' I','Xe','Cs','Ba','La','Ce','Pr','Nd', . 'Pm','Sm','Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb', . 'Lu','Hf','Ta',' W','Re','Os','Ir','Pt','Au','Hg', . 'Tl','Pb','Bi','Po','At','Rn','Fr','Ra','Ac','Th', . 'Pa',' U','Np','Pu','Am','Cm','Bk','Cf','Es','Fm', . 'Md','No','Lr','Rf','Db','Sg','Bh','Hs','Mt','Ds' / C rewind (ii1) C -- read in translation vectors, convert into Ang: do ii=1,3 read (ii1,101) (cc_bohr(jj,ii),jj=1,3) enddo cc_ang = cc_bohr*b2ang 101 format(3x,3f18.9) read (ii1,'(i12)') nat1 if (nat1.ne.nat) then C check if the same as returned by test_xv : write (6,*) ' Number of atoms from first and second ', . ' reading of XV differ !!' stop endif do iat=1,nat read (ii1,103,end=801,err=801) ityp(iat),iz(iat),coor_bohr 103 format(i3,i6,3f18.9) mass(iat)=amass(iz(iat)) label(iat)=alab(iz(iat)) do ii=1,3 coor_ang(ii,iat) = coor_bohr(ii)*b2ang enddo enddo return 801 write (6,*) ' End/Error reading XV for atom number ',iat stop end C 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............................................................... C subroutine opnout(io1,outfil) C C Opens output file, preventing its accidental overrid C implicit none integer io1 character outfil*60,owrite*1 logical filexist inquire (file=outfil, exist=filexist) if (filexist) then write (6,*) ' File ',outfil(1:len_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',err=801) else write (6,*) ' Then rename is first. Bye...' stop endif else open (io1,file=outfil,form='formatted',status='NEW',err=802) endif return 801 write (6,*) ' Error replacing formatted file ', . outfil(1:len_trim(outfil)) stop 802 write (6,*) ' Error opening file ', . outfil(1:len_trim(outfil)),' as new formatted' stop end C C............................................................... C subroutine makelin(obeg,oend,leng) C C asks for start and end of the linear path C C Input: none (interactive in/out) C Output: obeg - first point of path C oend - end point of path C leng - length of path in Angstroem C implicit none integer ii double precision obeg(3),oend(3),leng,b2ang parameter (b2ang=0.529177) ! Bohr to Angstroem character unitlab*1,labunit*4 logical unitb write (6,702) 101 write (6,703,advance="no") read (5,*) unitlab if (unitlab.eq.'B'.or.unitlab.eq.'b') then unitb = .true. labunit = 'Bohr' elseif (unitlab.eq.'A'.or.unitlab.eq.'a') then unitb = .false. labunit = 'Ang ' else write (6,*) ' Sorry, no third choice.' goto 101 endif write (6,704,advance="no") labunit read (5,*) (obeg(ii),ii=1,3) write (6,705,advance="no") labunit read (5,*) (oend(ii),ii=1,3) if (unitb) then C transform anything into Ang, as is standard in XCrysden: obeg=obeg*b2ang oend=oend*b2ang endif leng = sqrt( (oend(1)-obeg(1))**2 + . (oend(2)-obeg(2))**2 + . (oend(3)-obeg(3))**2 ) return 702 format(" Now define the path for scanning the mesh property.",/ . " We'll define it by the start and end points,", . " each given by three coordinates. They can be in Bohr or Ang.") 703 format (" Would you use Bohr (B) or Ang (A) ? ") 704 format(' Enter start point in ',a4,' : ') 705 format(' Enter end point in ',a4,' : ') end C C ........................................................... C subroutine intpl04 (func,fintp,mesh01,mesh02,mesh03,nrela,drela) C C 4-points linear interpolation on a 3-dim. grid. C The argument point (drela) falls into the microcell (nrela). C 0 < drela(i) < 1. Of eight edge points of the microcell, C keep the four closest, and perform a linear interpolation C over this tetrahedron. C C mesh01, C mesh02, C mesh03 : divisions along X-, Y- and Z- box sides C nrela(1..3) : index of microcell C drela(1..3) : relative coordinates of point in the microcell C func : function to be interpolated C fintp : interpolation result C integer mesh01,mesh02,nrela(3),inear(3),id(3), . igrid0(3),igrid1(3),igrid2(3),igrid3(3) integer ii,ix double precision drela(3),xx(3) real F0,F1,F2,F3 real func(mesh01,mesh02,mesh03),fintp ! NB! single precision C C select 4 interpolation points: inear=int(2*drela) ! inear(1,2,3): closest corner of mcrocell xx = drela - float(inear) C xx(1..3) are coordinates of interpolation point in its microcell, C relatvely to the closest microcell corner. Therefore C xx(1..3) can be between -1 and 1. id = int(sign(1.d0,xx)) C id(1..3) are values of xx(1..3) extended to the next integer C value, i.e. either -1 or 1. Therefore the nearest grid point C inear() and its three next ones displaced by id() do hopefully C make the tetrahedron which contains the mesh point xx(). C Now we want to address function values in four grid points, C igrid0(), igrid1(), igrid2(), igrid3(). They may correspond to C a displacement beyond the initial grid booundaries. Take care C of this introducing mod(): igrid0(1) = modulo(nrela(1)+inear(1)-1,mesh01)+1 igrid0(2) = modulo(nrela(2)+inear(2)-1,mesh02)+1 igrid0(3) = modulo(nrela(3)+inear(3)-1,mesh03)+1 igrid1(1) = modulo(nrela(1)+inear(1)-1+id(1),mesh01)+1 igrid1(2) = modulo(nrela(2)+inear(2)-1,mesh02)+1 igrid1(3) = modulo(nrela(3)+inear(3)-1,mesh03)+1 igrid2(1) = modulo(nrela(1)+inear(1)-1,mesh01)+1 igrid2(2) = modulo(nrela(2)+inear(2)-1+id(2),mesh02)+1 igrid2(3) = modulo(nrela(3)+inear(3)-1,mesh03)+1 igrid3(1) = modulo(nrela(1)+inear(1)-1,mesh01)+1 igrid3(2) = modulo(nrela(2)+inear(2)-1,mesh02)+1 igrid3(3) = modulo(nrela(3)+inear(3)-1+id(3),mesh03)+1 F0 = func( igrid0(1),igrid0(2),igrid0(3) ) F1 = func( igrid1(1),igrid1(2),igrid1(3) ) F2 = func( igrid2(1),igrid2(2),igrid2(3) ) F3 = func( igrid3(1),igrid3(2),igrid3(3) ) C recover interpolated function value: fintp = F0 + (F1-F0)*abs(xx(1)) + + (F2-F0)*abs(xx(2)) + + (F3-F0)*abs(xx(3)) return end