Dear Siesta community,
attached is a small tool designed to integrate a grid property
(I had primarily RHO in mind) over a sphere of given side,
centered at a given atom.
The motivation was to enable a ~ direct comparison with
other codes where e.g. local magnetic moments are printed out
relative to within a given sphere. The Mulliken analysis in Siesta
may provide gravely different numbers, and sometimes one would like to know
where the difference comes from, even as the numbers as such
do not have much sense.
The tool is standalone, interactive and "self-explaining".
The algorithm is exremely simplistic:
try all points on the grid within +/- one cell displaced along
each basis vector, retain all grid points which fall within
the given sphere, and sum up the corresponding values of grid property,
multiplied by the volume of grid microvolume. So this is in fact
a simple summation, but given a usually dense enough grid and
no high accuracy needed, "it works".
Output is just the charge and spin, with all non-collinearity neglected.
The tool is also available on my web page,
http://www.home.uni-osnabrueck.de/apostnik/download.html
Bug reports and suggestions are welcome.
Best regards,
Andrei Postnikov
+-- Dr. Andrei Postnikov ---- Tel. +33-387315873 ----- mobile +33-666784053 ---+
| Paul Verlaine University - Institute de Physique Electronique et Chimie, |
| Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz, France |
+-- [EMAIL PROTECTED] ------ http://www.home.uni-osnabrueck.de/apostnik/ --+
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 [EMAIL PROTECTED]
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
