Dear Miguel and Rudy (and all),
now that we are discussing this topic, I also have a script
to convert RHO and other grid properties into XCrysden format.
I attach a Fortran source and a README file.
Your suggestions and criticisms are welcome.
(Miguel: you did not attach your script in your mail,
at least it didn't go to the list).
For whoever is interested, I also have a tool to plot Fermi surfaces
from SIESTA data with XCrysden.
As you may figure out, a tool to transfer wavefunctions with XCrysden is
not yet ready.
Best regards,
Andrei Postnikov
+-- Dr. habil. Andrei Postnikov ----- Tel. +49-541-969.2377 -- Fax .2351 ---+
| Institut fuer Festkoerperforschung - FZ Juelich, D-52425 Juelich, Germany |
+-- [EMAIL PROTECTED] --------- http://www.home.uni-osnabrueck.de/apostnik/ --+
On Mon, 25 Apr 2005, Miguel Pruneda wrote:
| Dear Rudy (and all),
|
| As far as I know, the format "cube" is only valid for orthogonal cells. There
| are other options to plot .RHO (and .DRHO, .VH, .VT, etc). Some of them are
| included in the siesta/Utils directory. Check the mailing-list history,
because
| this is a frequent topic of discussion.
|
| You can also use the format of XCrysden (called xfs) to plot 3D surfaces. I
| attach a script similar to grid2cube.f, written in f77, that you can use in
the
| same way as grid2cube. You can download XCrysden at: http://www.xcrysden.org/
|
| Cheers,
| Miguel
|
| Missatge citat per Rudy Coquet <[EMAIL PROTECTED]>:
|
| > Dear SIESTA users,
| >
| > I would like to know if there is an easy way to convert a .RHO file to a
| > cube file? Usually grid2cube works fine but I have a non-orthogonal
| > cell so I cannot use it.
| >
| > Thank you,
| >
| > Rudy Coquet
| > Cardiff Universtity, UK.
| >
|
|
|
|
| -------------------------------------------------
| This mail sent through IMP: http://horde.org/imp/
|
|
C
C rho2xsf, a script to transform 3-dim grid function
C (i.e. LDOS, RHO, DRHO, etc.) written in SIESTA
C by subr. iorho
C into a grid for XCrysden
C
C !!! -------------------- IMPORTANT ---------------------- !!!
C compile this code with the same compiler switches as Siesta,
C otherwise reading the data from unformatted files will be spoiled.
C
C Written by Andrei Postnikov, Mar 2005
C [EMAIL PROTECTED]
C
program rho2xsf
implicit none
integer MPTS,mesh0(3),mesh1(3),ip,nspin,is,ii,jj,
. iat,nat,ityp,nz,ix,iy,iz,ind,mn,
. ixmax,iymax,izmax,ixmin,iymin,izmin
parameter (MPTS=10000)
character inpfil*60,outfil*60,syslab*30,suffix*6
real*8 cell(3,3),coord(3),b2ang
parameter (b2ang=0.529177) ! Bohr to Angstroem
real func(MPTS),fdum ! NB! single precision, as in iorho.F
real fmax,fmin
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
write (6,701)
701 format(" Specify SystemLabel (or 'siesta' if none): ",$)
read (5,*) syslab
inpfil = syslab(1:len_trim(syslab))//'.XV'
open (11,file=inpfil,form='formatted',status='old',err=801)
write (6,*) 'Found and opened: ',inpfil
C reads from .XV and writes into .XSF file in XCrysden format:
outfil = syslab(1:len_trim(syslab))//'.XSF'
open (12,file=outfil,form='formatted',status='new',err=802)
write (6,*) 'Opened as new: ',outfil
C --- write crystal structure data:
write (12,'(a7)') 'CRYSTAL'
write (12,'(a7)') 'PRIMVEC'
do ii=1,3
read (11,*,end=803,err=803) (cell(ii,jj),jj=1,3)
write (12,204) (cell(ii,jj)*b2ang,jj=1,3)
enddo
write (12,'(a9)') 'PRIMCOORD'
read (11,*,end=804,err=804) nat
write (12,'(i4,a3)') nat,' 1'
do iat=1,nat
read (11,*,end=805,err=805) ityp, nz, (coord(ii),ii=1,3)
write (12,201) nz, (coord(ii)*b2ang,ii=1,3)
enddo
close (11)
C --- finished with .XV; now look for grid data files to include:
C
101 write (6,702)
702 format (' Add grid property (LDOS, RHO, ...;',
. ' or BYE if none): ',$)
read (5,*) suffix
inpfil = syslab(1:len_trim(syslab))//
. '.'//suffix(1:len_trim(suffix))
open (11,file=inpfil,form='unformatted',status='old',err=806)
write (6,*) 'Found and opened: ',inpfil(1:len_trim(inpfil))
read (11,err=807) cell
read (11,err=808) mesh0, nspin
C
C --- Here, introdruce mesh step (keep every mn's point along
C each grid dimension).
11 write (6,'(a18,3i5)') 'Grid divisions: ',mesh0
write (6,'(a39,$)') " Keep for XSF file every N'th point, N="
read (5,*) mn
if (mn.le.0) then
write (6,*) ' N must be positive, try again'
goto 11
endif
if (mod(mesh0(1),mn).ne.0.or.
. mod(mesh0(2),mn).ne.0.or.
. mod(mesh0(3),mn).ne.0) then
write (6,*) ' Must be a common divisor of all three, try again'
goto 11
endif
mesh1(1)=mesh0(1)/mn
mesh1(2)=mesh0(2)/mn
mesh1(3)=mesh0(3)/mn
if (mesh1(1).gt.MPTS) then
write (6,*) ' Number of grid points along 1st cell vector =',
. mesh1(1),' exceeds MPTS=',MPTS
stop
endif
write (12,'(a23)') 'BEGIN_BLOCK_DATAGRID_3D'
write (12,*) 'DATA_from:'//inpfil(1:len_trim(inpfil))
do is=1,nspin
fmax= -9.999E+10
fmin= 9.999E+10
write (12,*) 'BEGIN_DATAGRID_3D_'//
. suffix(1:len_trim(suffix))//':spin_'//char(48+is)
write (12,203) mesh1(1)+1, mesh1(2)+1, mesh1(3)+1
write (12,204) 0., 0., 0.
write (12,204) ((cell(ii,jj)*b2ang,jj=1,3),ii=1,3)
C
C --- for the following it is important that XCrysDen uses
C general grids (i.e. those with redundant (translated) points
C and Siesta - periodic grids. (See explanation in the XCrysDen
C documentation, www.xcrysden.org
C -> Documentation
C -> Specification of the XSF Format
C -> Specification of 2D and 3D scalar-fields
C -> General vs. Periodic grids
C Therefore it is necessary to copy the data
C from ix=1 into ix=mesh1(1)+1
C and similarly for iy and iz. This is organized below
C using two temp files: unit (13) saves a vector of data
C [ix = 1,...,mesh1(1)] for iy=1 and releases it after a loop in iy;
C unit (14) saves a plane of data
C [iy=1,...,mesh1(2)+1,[ix=1,...,mesh1(1)]] for iz=1
C and releases it after the loop in iz.
C
C --- read data as written by subr. iorho:
open (13,status='scratch',form='unformatted')
open (14,status='scratch',form='unformatted')
do iz=1,mesh0(3)
if (iz.eq.1) rewind (14)
if (mod(iz-1,mn).ne.0) then
C --- read and forget this iz-layer :
do iy=1,mesh0(2)
read (11,err=809,end=810) (fdum,ix=1,mesh0(1))
enddo
else
C --- this iz-layer contains useful data; run through iy :
do iy=1,mesh0(2)
if (mod(iy-1,mn).ne.0) then
C read and forget this iy-vector :
read (11,err=809,end=810) (fdum,ix=1,mesh0(1))
else
C --- read and keep this iy-vector, or, more precisely,
C keep its 1, 1+mn, 1+2*mn, etc. element :
read (11,err=809,end=810) (func(ix), ! main
. (fdum,ii=1,mn-1),ix=1,mesh1(1)) ! data read
write (12,205) (func(ix),ix=1,mesh1(1)),func(1)
C --- mark max. and min. values on the grid ... (just for info):
do ix=1,mesh1(1)
if (func(ix).gt.fmax) then
ixmax=ix
iymax=(iy-1)/mn+1
izmax=(iz-1)/mn+1
fmax=func(ix)
endif
if (func(ix).lt.fmin) then
ixmin=ix
iymin=(iy-1)/mn+1
izmin=(iz-1)/mn+1
fmin=func(ix)
endif
enddo
C --- ... up to here
if (iy.eq.1) then
rewind (13)
write(13) (func(ix),ix=1,mesh1(1))
endif
if (iz.eq.1) write(14) (func(ix),ix=1,mesh1(1))
endif ! (mod(iy-1,mn).ne.0)
enddo ! do iy=1,mesh0(2)
C copy all (ix) values from iy=1 vector to iy=mesh1(2)+1 :
rewind (13)
read(13) (func(ix),ix=1,mesh1(1))
write (12,205) (func(ix),ix=1,mesh1(1)),func(1)
if (iz.eq.1) write(14) (func(ix),ix=1,mesh1(1))
endif ! (mod(iz-1,mn).ne.0)
enddo ! do iz=1,mesh0(3)
C --- copy all (ix,iy) values from iz=1 plane to iz=mesh1(3)+1 :
rewind (14)
do iy=1,mesh1(2)+1
read (14) (func(ix),ix=1,mesh1(1))
write (12,205) (func(ix),ix=1,mesh1(1)),func(1)
enddo
close (13)
close (14)
write (12,*) 'END_DATAGRID_3D'
write (6,206) is,'max',fmax,ixmax,iymax,izmax
write (6,206) is,'min',fmin,ixmin,iymin,izmin
enddo ! do is=1,nspin
write (12,'(a21)') 'END_BLOCK_DATAGRID_3D'
close (11)
goto 101
201 format (i4,3f20.8)
203 format (3i6)
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
804 write (6,*) ' End/Error reading XV for number of atoms line'
stop
805 write (6,*) ' End/Error reading XV for atom number ',iat
stop
806 write (6,*) ' A wild guess! There is no file ',
. inpfil(1:len_trim(inpfil)),'; close XSF and quit.'
close (12)
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 iy=',iy,
. ' iz=',iz,' is=',is
stop
end
This directory contains routines written by Andrei Postnikov
(Forschungszentrum Juelich; [EMAIL PROTECTED])
for transforming SIESTA properties files into input files
of the XCrySDen package. XCrySDen is being developed
by Anton Kokalj at the Jozef Stefan Institute in Ljubljana, Slovenia.
XCrySDen is released under the GNU General Public License,
its web site is http://www.xcrysden.org/
It offers a straightforward installation under Unix/Linux,
easy user's interface and gascinating graphics, enabling to manipulate
crystal structures and their 2-dim. and 3-dim. functions as 2-dim. plots
or 3-dim. isosurfaces. Its input data format is simple and well documented
on the XCrySDen site
[ Home -> Documentation -> Description of XCrySDen's XSF format ] .
The need of 3-dimensional manipulation of Siesta data falls into
four categories:
1. Viewing of crystal structures.
2. Viewing data stored on the grid spanning the simulation cell:
RHO (switch SaveRho );
IORHO (switch SaveDeltaRho );
VH (switch SaveElectrostaticPotential );
VT (switch SaveTotalPotential );
IOCH (switch SaveIonicCharge );
TOCH (switch SaveTotalCharge );
LDOS (block LocalDensityOfStates )
These data are written in subr. iorho (iorho.F) unformatted,
and post-processing recommended by the program plrho. The present script
rho2xcf offers a convenient alternative.
3. Viewing data stored as decomposed over basis functions,
like wavefunctions, ...
4. Viewing Fermi surfaces.
--------------------- TASKS 1) AND 2) --------------------------------
are done using the script rho2xcf, which works
interactively and must be started without passing any parameters to it.
If you need just to vizualize crystal/ molecular structure from
a Siesta run, you need to have .XV file which will be written into the
.XCF file in the correct format. For vizualizing 3-dim. properties
one needs to have a corresponding property file(s) in the work directory.
Note that many properties may be put into a single .XCF file, if needed.
Note that rho2xcf makes unformatted read-in, according to the sequence of
(unformatted) records done in iorho.F. Therefore in the future - in case of
problems - one might need to check the consistency of record sequences.
Moreover it is ESSENTIAL to compile rho2xcf.f using the same compiler
and parameters as for Siesta (apart from parallelization issues), to ensure
the consistency of unformatted record lengths: there are variables defined
as real and as real*8 in both rho2xcf.f and iorho.F, and compiler settings
may influence the actually used word lengths, e.g. imposing double precision
throughout.
The program rho2xcf has to be simply started in the directory which contains
the .XV and .LDOS etc. files; the rest is interactive. Since the grid
might be too dense for plotting purposes (overkill in mesh points,
resulting in very large formatted files and slow graphic manipulation
by XCrySDen), an option is provided to keep every N'th point of grid
values, where N is a common divisor of three FTT mesh divisions.
Moreover it is important to understand that XCrySDen uses general grids,
i.e. those with redundant (translated) points at the cell edges,
and Siesta - periodic grids, with no redundant points.
See explanation in the XCrySDen documentation, www.xcrysden.org
-> Documentation
-> Specification of the XSF Format
-> Specification of 2D and 3D scalar-fields
-> General vs. Periodic grids
This is taken care of in rho2xcf by copying the function values of the grid:
f(1,iy,iz) into f(Nx+1,iy,iz),
f(ix,1,iz) into f(ix,Ny+1,iz),
f(ix,iy,1) into f(ix,iy,Nz+1). The grid points #1 and #N?+1 lie then
in the (equivalent) limiting planes of the simulation cell.
Note that XCrySDen permits plotting isosurfaces on non-orthogonal grids.
The way to use the generated .XCF file:
Choose in the XCrySDen Menu: File -> Open Structure -> Open XCF
This reads in crystal structure information and shows the atoms of a unit cell.
If there is a scalar field (3-dim. property) an the XCF file,
the submenu "Tools" of the XCrySDen main menu shows an active entry
"Data Grid", which opens a map of grid data available in he XCF file.
Choose a Block (if there are many, like LDOS, RHO etc.) to plot,
and within the Block one can activate-deactivate Sub-blocks and sum their
contents up with arbitrary "Multiply factor"s (in the windows).
rho2xsf makes two blocks whenever there are two spins in the Siesta
property files. Note that XCrySDen activates by default only the first
Sub-block. One may wish to activate both (for summing over spins)
and set on of "Multiply factor"s to -1.0 (for analyzing e.g. spin density).
------------- TASK 4 (Fermi surface plotting) -------------------------
is done using the script ene2bxsf, which manipulates the eigenvalues
previously calculated by Siesta on a (fine enough) k-mesh.
The script must be started in the directory which contains
the files from a Siesta calculation; after inquiring for the Siesta
SystemLabel, the script opens the following files:
...XV (to read the lattice vectors),
...KP (to read the list of irreducible k-points),
...EIG (to read the eigenvalues over k-points and bands).
The resulting file(s)
...BXSF (two for a spin-polarized calculation)
contain(s) the ordered list of eigenvalues, band by band, over a full
k-mesh, consistently with the rules for constructing bandgrids
for XCrySDen, see www.xcrysden.org
-> Documentation
-> Specification of the XSF Format
-> Bandgrids (visualization of Fermi surfaces).
The visualization proceeds by invoking XCrySDen,
File -> Open Structure -> Open BXSF
then confirm (or change) the Fermi energy,
select bands for Fermi surface drawing.
The access to the graphic menu for visual control
is via the right mouse button. See details in the XCrySDen Documentation
-> Description of XCrySDen's Fermi surface viewer.
I find it a bit annoying that the option
Display -> Depth Cuing
is acctivated by default,
and apparently there is no way to change the width (or color) of
the Brillouin zone wire, but one can live with this.
Two important issues:
1. The origin of a bandgrid for XCrySDen must be ai Gamma, therefore
the k-mesh used for the calculation should not be shifted. This is
actually tested in ene2bxsf, as the list of irreducible k-points is mapped
onto the general bandgrid.
2. The resulting .BXSF file contains only the bands which are actually crossing
the Fermi level, whose value is read from the first line of .EIG.
The subsequent tuning of the Fermi energy is possible in the "Isolevel"
window of the "XCrySDen-Fermi Surface" program. However, if you need
to choose other band(s) for potting you must change the Fermi energy
in the first line of .EIG and re-run the ene2bxsf script.
