Hi,
I have written a small fortran program to calculate yz averages as funciton of
x ect. based on the
grid2cube.f in Siesta 2.0.2 utils directory. Note that for reading unformatted
files exactly the same compiler has to be used as for siesta.
Regards,
Ulrich.
----- Original Message -----
From: root
[mailto:[email protected]]
To: [email protected]
Sent: Wed, 20 Apr
2011 05:35:33 +0200
Subject: [SIESTA-L] The average potential along certain
direction
> Hi all,
>
> Does anyone has experience of plotting the average potential along certain
> direction after a SIESTA run? I found that XCrysden can plot the 3D real
> space potential profile, however what I want is the average potential
> changes along a certain direction. Any suggestion for this?
>
> Thanks a lot.
>
---
Max-Planck-Institut für Eisenforschung GmbH
Max-Planck-Straße 1
D-40237 Düsseldorf
Handelsregister B 2533
Amtsgericht Düsseldorf
Geschäftsführung
Prof. Dr. Jörg Neugebauer
Prof. Dr. Dierk Raabe
Prof. Dr. Martin Stratmann
Dipl.-Kfm. Herbert Wilk
Ust.-Id.-Nr.: DE 11 93 58 514
Steuernummer: 105 5891 1000
-------------------------------------------------
program line_av
c****************************************************************************
c LINE_AV version 1.1.1
c
c written by P. Ulrich Biedermann March 2009
C last changed 30.3.2009
c This program LINE_AV reads files written by SIESTA with data stored on the
c mesh used by DHSCF, e.g., RHO (.RHO), Delta RHO (.DRHO), Electrostatic
c Potential (.VH) and total potential (.VT). It computes the averages over
c the yz planes of the data and writes them as function of x to file.
c systemlabel.TASK.line
c And so on for y and z directions.
c
c based on GRID2CUBE written by P. Ordejon, August 2002
c Last modification: June 2003
! GRID2CUBE is part of the SIESTA package.
!
! Copyright (c) Fundacion General Universidad Autonoma de Madrid:
! E.Artacho, J.Gale, A.Garcia, J.Junquera, P.Ordejon, D.Sanchez-Portal
! and J.M.Soler, 1996-2006.
!
! Use of this software constitutes agreement with the full conditions
! given in the SIESTA license, as signed by all legitimate users.
!
c****************************************************************************
c CAVEATS:
c
c This version only works with ORTHOGONAL cells, where the cell vectors
c are in the X, Y and Z directions. Extension to non-orthogonal cells
c might be done in future versions.
c****************************************************************************
c USAGE:
c
c This program reads files generated from SIESTA, with information of
c physical quantities on the real space grid, such as charge densities
c (filename.RHO, filename.DRHO, filename.LDOS) or potentials (filename.VH,
c filename.VT) and and computes a profile along the z coordinate
c
c The program needs two input files:
c
c 1) Main input file, read by standard input. A sample of input file is:
c
c --- begin input file ---
c h2o
c rho
c 4.0 6.0 5.0
c 2
c unformatted
c --- end input file ---
c
c where:
c - The first line is the label of the system, as in SIESTA SystemLabel.
c Files will be searched as SystemLabel.* (in the example, h2o.RHO).
c - The second line is the task, which should be one of the following:
c rho, drho, ldos, vh or vt (in lowercase!!).
c - The third line is a shift of the origin of coordinates (in bohr).
c This is useful, for instance, when plotting molecules, if the
c center of the molecule is close to the simulation cell walls.
c In that case, the molecule would appear split, due to the
c contributions from multiple images. Shifting the molecule will
c make it to appear centered in the image, and make it look
c as a single object
c (currently not implemented)
c - The forth line is an integer (nskip) that specifies the density of
c grid points in the output. For each direction, only one of every
c 'nskip' points from the input will be written to the output.
c nskip = 1 will yield an output with the same number of grid points
c than the input; nskip = N will provide an output N**3 less dense
c than the input. This option is provided because some
c visualization programs can not handle the large memory involved
c in large systems with very fine grids, so the grid must be
c made coarser.
c - The fifth line indicates if the grid data file SystemLabel.TASK
c is formatted or unformatted (the latter being the standard option
c in SIESTA)
c
c
c 2) SystemLabel.TASK file: this is a file generated by SIESTA, with
c the values of the appropriate quantity on the grid.
c In example above: h2o.RHO. You should copy it from the directory
c with your SIESTA output files.
c
c
c The program generates some informative output on standard output,
c and writes one or two files with the converted grid data.
c If your SIESTA calculation was not spin polarized, there will be only
c one output file: SystemLabel.TASK.cube with the information in Gaussian
c CUBE format.
c If your SIESTA calculation was spin polarized, there will be two
c output files: SystemLabel.TASK.UP.line and SystemLabel.TASK.DN.line,
c with the information of the physical quantity for the two spin components,
c in Gaussian CUBE format (except for the Hartree potential, VH, which
c has the same value for both spins, and therefore is only plotted once
c in SystemLabel.VH.cube).
c****************************************************************************
implicit none
integer maxg, nskip
parameter (maxg = 300)
integer ipt, isp, ix, iy, iz, i, ip, natoms, np,
. mesh(3), nspin, Ind, id(3), iix, iiy,
. iiz, ii, length, lb
character sysname*70, fnamein*75, fnameout(2)*75,
. fnamexv*75, paste*74, task*5, fform*12
real rho(maxg,maxg,maxg,2), line(maxg,3)
double precision cell(3,3), cm(3), rt(3),
. delta(3), dr(3), residual
external paste, lb
c ---------------------------------------------------------------------------
read(5,*) sysname
read(5,*) task
read(5,*) rt(1),rt(2),rt(3)
read(5,*) nskip
read(5,*) fform
fnamexv = paste(sysname,'.XV')
if (task .eq. 'rho') then
fnamein = paste(sysname,'.RHO')
else if (task .eq. 'drho') then
fnamein = paste(sysname,'.DRHO')
else if (task .eq. 'ldos') then
fnamein = paste(sysname,'.LDOS')
else if (task .eq. 'vt') then
fnamein = paste(sysname,'.VT')
else if (task .eq. 'vh') then
fnamein = paste(sysname,'.VH')
else
write(6,*) 'Wrong task'
write(6,*) 'Accepted values: rho, drho, ldos, vh, vt'
write(6,*) '(in lower case!!!!)'
stop
endif
length = lb(fnamein)
write(6,*)
write(6,*) 'Reading grid data from file ',fnamein(1:length)
c read function from the 3D grid --------------------------------------------
open( unit=1, file=fnamein, form=fform, status='old' )
if (fform .eq. 'unformatted') then
read(1) cell
else if (fform .eq. 'formatted') then
do ix=1,3
read(1,*) (cell(iy,ix),iy=1,3)
enddo
else
stop 'ERROR: last input line must be formatted or unformatted'
endif
write(6,*)
write(6,*) 'Cell vectors'
write(6,*)
write(6,*) cell(1,1),cell(2,1),cell(3,1)
write(6,*) cell(1,2),cell(2,2),cell(3,2)
write(6,*) cell(1,3),cell(2,3),cell(3,3)
residual = 0.0d0
do ix=1,3
do iy=ix+1,3
residual = residual + cell(ix,iy)**2
enddo
enddo
if (residual .gt. 1.0d-6) then
write(6,*)
write(6,*) 'ERROR: this progam can only handle orthogonal cells'
write(6,*) ' with vectors pointing in the X, Y and Z directions'
stop
endif
if (fform .eq. 'unformatted') then
read(1) mesh, nspin
else
read(1,*) mesh, nspin
endif
write(6,*)
write(6,*) 'Grid mesh: ',mesh(1),'x',mesh(2),'x',mesh(3)
write(6,*)
write(6,*) 'nspin = ',nspin
write(6,*)
do ix=1,3
dr(ix)=cell(ix,ix)/mesh(ix)
enddo
np = mesh(1) * mesh(2) * mesh(3)
if (np .gt. maxg*maxg*maxg) stop
. 'grid2d: Parameter MAXG too small'
C do isp=1,nspin
C Ind=0
C if (fform .eq. 'unformatted') then
C do iz=1,mesh(3)
C do iy=1,mesh(2)
C read(1) (rho(Ind+ix,isp),ix=1,mesh(1))
C Ind=Ind+mesh(1)
C enddo
C enddo
C else
C do iz=1,mesh(3)
C do iy=1,mesh(2)
C read(1,'(e15.6)') (rho(Ind+ix,isp),ix=1,mesh(1))
C Ind=Ind+mesh(1)
C enddo
C enddo
C endif
C enddo
if (fform .eq. 'unformatted') then
do isp=1,nspin
do iz=1,mesh(3)
do iy=1,mesh(2)
read(1) (rho(ix,iy,iz,isp),ix=1,mesh(1))
enddo
enddo
enddo
else
do isp=1,nspin
do iz=1,mesh(3)
do iy=1,mesh(2)
read(1,'(e15.6)') (rho(ix,iy,iz,isp),ix=1,mesh(1))
enddo
enddo
enddo
endif
C x-direction average
if (nspin .eq. 1) then
fnameout(1) = paste(fnamein,'.line.x')
else if (nspin .eq. 2) then
fnameout(1) = paste(fnamein,'.UP.line.x')
fnameout(2) = paste(fnamein,'.DN.line.x')
else
stop 'nspin must be either 1 or 2'
endif
do isp=1,nspin
length = lb(fnameout(isp))
write(6,*) 'Writing LINE file ',fnameout(isp)(1:length)
C open( unit=2, file=fnameout(isp), status='new', form='formatted')
open( unit=2, file=fnameout(isp), form='formatted')
length = lb(fnameout(isp))
write(2,*) fnameout(isp)(1:length)
do ix=1,mesh(1),nskip
line(ix,1)=0
do iy=1,mesh(2),nskip
do iz=1,mesh(3),nskip
line(ix,1)=line(ix,1)+rho(ix,iy,iz,isp)
enddo
enddo
line(ix,1)=line(ix,1)*nskip*nskip/mesh(2)/mesh(3)
write(2,'(f10.6,e13.5)')
. (ix-1)*cell(1,1)/mesh(1),line(ix,1)
enddo
close(2)
enddo
write(6,*)
C y-direction average
if (nspin .eq. 1) then
fnameout(1) = paste(fnamein,'.line.y')
else if (nspin .eq. 2) then
fnameout(1) = paste(fnamein,'.UP.line.y')
fnameout(2) = paste(fnamein,'.DN.line.y')
else
stop 'nspin must be either 1 or 2'
endif
do isp=1,nspin
length = lb(fnameout(isp))
write(6,*) 'Writing LINE file ',fnameout(isp)(1:length)
C open( unit=2, file=fnameout(isp), status='new', form='formatted')
open( unit=2, file=fnameout(isp), form='formatted')
length = lb(fnameout(isp))
write(2,*) fnameout(isp)(1:length)
do iy=1,mesh(2),nskip
line(iy,2)=0
do ix=1,mesh(1),nskip
do iz=1,mesh(3),nskip
line(iy,2)=line(iy,2)+rho(ix,iy,iz,isp)
enddo
enddo
line(iy,2)=line(iy,2)*nskip*nskip/mesh(1)/mesh(3)
write(2,'(f10.6,e13.5)')
. (iy-1)*cell(2,2)/mesh(2),line(iy,2)
enddo
close(2)
enddo
write(6,*)
C z-direction average
if (nspin .eq. 1) then
fnameout(1) = paste(fnamein,'.line.z')
else if (nspin .eq. 2) then
fnameout(1) = paste(fnamein,'.UP.line.z')
fnameout(2) = paste(fnamein,'.DN.line.z')
else
stop 'nspin must be either 1 or 2'
endif
do isp=1,nspin
length = lb(fnameout(isp))
write(6,*) 'Writing LINE file ',fnameout(isp)(1:length)
C open( unit=2, file=fnameout(isp), status='new', form='formatted')
open( unit=2, file=fnameout(isp), form='formatted')
length = lb(fnameout(isp))
write(2,*) fnameout(isp)(1:length)
do iz=1,mesh(3),nskip
line(iz,3)=0
do ix=1,mesh(1),nskip
do iy=1,mesh(2),nskip
line(iz,3)=line(iz,3)+rho(ix,iy,iz,isp)
enddo
enddo
line(iz,3)=line(iz,3)*nskip*nskip/mesh(1)/mesh(2)
write(2,'(f10.5,e13.5)')
. (iz-1)*cell(3,3)/mesh(3),line(iz,3)
enddo
close(2)
enddo
write(6,*)
end
CHARACTER*(*) FUNCTION PASTE( STR1, STR2 )
C CONCATENATES THE STRINGS STR1 AND STR2 REMOVING BLANKS IN BETWEEN
C Writen by Jose M. Soler
CHARACTER*(*) STR1, STR2
DO 10 L = LEN( STR1 ), 1, -1
IF (STR1(L:L) .NE. ' ') GOTO 20
10 CONTINUE
20 PASTE = STR1(1:L)//STR2
END
INTEGER FUNCTION LB ( STR1 )
C RETURNS THE SIZE IF STRING STR1 WITH BLANKS REMOVED
C Writen by P. Ordejon from Soler's paste.f
CHARACTER*(*) STR1
DO 10 L = LEN( STR1 ), 1, -1
IF (STR1(L:L) .NE. ' ') GOTO 20
10 CONTINUE
20 LB = L
END
