On Wed, 31 May 2006, Sharat Chandra wrote:
| Hi
| Can any one please tell me how to calculate the phonon DOS using Vibra
| module? Has any one developed a code for carrying out the calculation?
Hallo Sharat:
yes I have a tool to calculate (local) phonon DOS,
however it is of limited usefulness:
it assumes a single q-point in a presumable large supercell,
and smears out the peaks at frequencies (and with weights)
as provided by the "vibrator".
Please find enclosed the source code and README.
Best regards,
Andrei
+-- 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 phdos, a script to calculte phonon density of states
C in a supercell,
C using phonon eigenvectors provided by "vibrator"
C Written by Andrei Postnikov, Oct 2005 Vers_0.2
C XV read format bug fixed May 2006
C [EMAIL PROTECTED]
C
program phdos
implicit none
integer ii1,ii2,io1,io2,ivmin,ivmax
parameter (ii1=11,ii2=12,io1=14,io2=15)
integer ialloc,iat,nat,iatmin,iatmax,mode,nodif,idif,npoints
integer, allocatable :: ityp(:),iz(:),jat(:,:),jtyp(:),nna(:)
double precision tau(3,3),qvec(3),qq(3),delta
double precision, allocatable :: mass(:),coor(:,:),
. freq(:),evr(:,:,:),evi(:,:,:),wwsum(:,:),zz(:)
character inpfil*60,outfil*60,syslab*30,qlab*1
character*2, allocatable :: label(:)
external test_xv,read_xv,read_ev,full_dos,qres_dos
C
C string manipulation functions in Fortran used below:
C len_trim(string): returns the length of string
C without trailing blank characters,
C --- read lattice vectors and atom positions from the .XV file:
write (6,701)
701 format(' Specify SystemLabel of .XV file: ',$)
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 (jtyp(nat))
allocate (nna(nat))
allocate (jat(nat,nat))
allocate (coor(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,tau,mass,label,coor)
close (ii1)
C --- count different types in the XV file. Note that they can be changed
C arbitrarily (by hand) in the XV file as compared to the Siesta run.
C E.g. to select some atoms of type=4 as different, assign them type=14.
C They will then get a separate column in output files.
nodif=0
do 98 iat=1,nat
if (nodif.gt.0) then
do idif=1,nodif
if (ityp(iat).eq.jtyp(idif)) then ! add atom to existing type
nna(idif)=nna(idif)+1
jat(idif,nna(idif))=iat
goto 98
endif
enddo
endif
nodif=nodif+1 ! add a new type
jtyp(nodif)=ityp(iat)
nna(nodif)=1
jat(nodif,1)=iat
98 continue
write (6,"(i5,' different types found in the XV file:')") nodif
do idif=1,nodif
write (6,"(i5,' atoms of type ',i5)") nna(idif),jtyp(idif)
enddo
C --- read and store frequencies /eigenvectors from the .vector file,
C q =(0 0 0) only :
C allocate for all modes, 1 through nat*3 :
ivmin=1
ivmax=nat*3
C (change ivmin, ivmax above if want to restrict to less modes)
allocate (freq(ivmin:ivmax))
allocate (evr(1:3,1:nat,ivmin:ivmax))
allocate (evi(1:3,1:nat,ivmin:ivmax))
allocate (zz(1:nodif))
allocate (wwsum(1:nodif,ivmin:ivmax),STAT=ialloc)
if (ialloc.ne.0) then
write (6,*) ' Fails to allocate space for vibration modes'
stop
endif
write (6,702)
702 format(' Specify SystemLabel of .vectors file: ',$)
read (5,*) syslab
inpfil = syslab(1:len_trim(syslab))//'.vectors'
open (ii2,file=inpfil,form='formatted',status='old',err=801)
call read_ev(ii2,nat,qq,ivmin,ivmax,evr,evi,freq)
write (6,*) 'opened and read ',inpfil
close (ii2)
C ---------------------------------------------------
write (6,*) 'You have two options: Total density of states ',
. '(sum of squares of eigenvectors),'
write (6,*) 'or Q-projected density of states ',
. '(convolution with a given Q-vector).'
101 write (6,703)
703 format(' Do you want total density of states (T) ',
. ' or convolution (Q) ? : ',$)
read (5,*,err=101,end=101) qlab
if (qlab.eq.'T'.or.qlab.eq.'t') then
mode=1
else if (qlab.eq.'Q'.or.qlab.eq.'q') then
mode=2
write (6,*) 'Type three numbers qx qy qz for convolution. ',
. 'The units are 1/Bohr.'
write (6,*) 'So for the X point of simple cubic lattice with ',
. 'A=5 Bohr type: 0.1 0 0'
102 write (6,704)
704 format (' Enter qx qy qz : ',$)
read (5,*,err=102,end=102) qvec
else
goto 101
endif
C ---------------------------------------------------
C --- output file(s) .WS? contain the same information as .VS?,
C but smeared, as an (continuous) functions of frequency.
C Fix smearing parameter and No. of steps, to avoid many questions:
C delta = 10.0
delta = 5.0
npoints = 2000
C
if (mode.eq.1) then
outfil = syslab(1:len_trim(syslab))//'.VST'
open (io1,file=outfil,form='formatted',status='new',err=802)
call full_peaks(io1,nat,jat,label,qq,nodif,jtyp,nna,
. ivmin,ivmax,freq,evr,evi,wwsum)
close (io1)
outfil = syslab(1:len_trim(syslab))//'.WST'
open (io2,file=outfil,form='formatted',status='new',err=802)
call full_smear(delta,npoints,io2,ivmin,ivmax,nat,jat,qq,
. nodif,jtyp,label,freq,wwsum,zz)
close (io2)
else if (mode.eq.2) then
outfil = syslab(1:len_trim(syslab))//'.VSQ'
open (io1,file=outfil,form='formatted',status='new',err=802)
call qres_peaks(io1,nat,jat,label,qq,nodif,jtyp,nna,coor,
. ivmin,ivmax,freq,evr,evi,qvec,wwsum)
close (io1)
outfil = syslab(1:len_trim(syslab))//'.WSQ'
open (io2,file=outfil,form='formatted',status='new',err=802)
call qres_smear(delta,npoints,io2,ivmin,ivmax,nat,jat,qq,
. nodif,jtyp,label,freq,wwsum,qvec,zz)
close (io2)
endif
C
deallocate (ityp,iz,mass,label,jtyp,nna,jat,coor)
deallocate (freq,evr,wwsum,zz)
stop
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(inpfil)),' as new'
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,nat
double precision dummy
rewind ii1
read (ii1,101) (dummy,ii=1,9)
101 format(3x,3f18.9)
C read (ii1,'(i13)') nat
read (ii1,*) nat
return
end
C
C...............................................................
C
subroutine read_xv(ii1,nat,ityp,iz,tau,mass,label,coor)
C
C reads again from ii1 (XV file) to the end
C
implicit none
integer ii1,nat,nat1,iat,ii,ityp(nat),iz(nat)
double precision tau(3,3),mass(nat),coor(3,nat),amass(110)
character*2 label(nat),alab(110)
C
data amass / 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
. 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)
read (ii1,101) (tau(ii,1),ii=1,3)
read (ii1,101) (tau(ii,2),ii=1,3)
read (ii1,101) (tau(ii,3),ii=1,3)
101 format(3x,3f18.9)
C read (ii1,'(i13)') nat1
read (ii1,*) 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) ityp(iat),iz(iat),(coor(ii,iat),ii=1,3)
103 format(i3,i6,3f18.9)
mass(iat)=amass(iz(iat))
label(iat)=alab(iz(iat))
enddo
return
end
C
C...............................................................
C
subroutine read_ev(ii2,nat,qq,ivmin,ivmax,evr,evi,freq)
C
C reads selected modes only from ii2 ('vector' file),
C check for consistency, and keeps selected modes only
C (ivmin to ivmax)
C
implicit none
integer ii2,ii,nat,iat,ivmin,ivmax,iev,idum
double precision qq(3),dummy,small,freq(ivmin:ivmax),
. evr(3,nat,ivmin:ivmax),evi(3,nat,ivmin:ivmax)
data small/10.e-4/
C
rewind (ii2)
read (ii2,100,err=301) ! leading line
read (ii2,101,err=301) (qq(ii),ii=1,3)
if (sqrt(qq(1)**2+qq(2)**2+qq(3)**2).gt.small) then
write (6,*) ' Attention: q.ne.0 in .vectors file!'
C stop
endif
do iev=1,nat*3
read (ii2,102,err=302,end=401) idum
if (idum.ne.iev) then
write (6,*) ' Stop in read_ev: expected eigenvector ',iev,
. ' but found ',idum
stop
endif
if (iev.lt.ivmin.or.iev.gt.ivmax) then ! -- skip these modes
read (ii2,103,err=303) dummy
read (ii2,100) ! Eigenvector, real part follows
do iat=1,nat
read (ii2,104,err=304) (dummy,ii=1,3)
enddo
read (ii2,100) ! Eigenvector, imag part follows
do iat=1,nat
read (ii2,104,err=305) (dummy,ii=1,3)
enddo
else ! --- selected modes...
read (ii2,103,err=303) freq(iev)
read (ii2,100) ! Eigenvector, real part follows
do iat=1,nat
read (ii2,104,err=304) (evr(ii,iat,iev),ii=1,3)
enddo
read (ii2,100) ! Eigenvector, imag part follows
do iat=1,nat
read (ii2,104,err=305) (evi(ii,iat,iev),ii=1,3)
enddo
endif
enddo
return
301 print *,' Error in 1st or 2d line of .vector file'
stop
302 print *,' Error reading Eigenvector line, last iev=',iev
stop
303 print *,' Error reading Frequency line, iev=',iev
stop
304 print *,' Error reading real part of eigenvector ',iev,
. ' iat=',iat
stop
305 print *,' Error reading imag part of eigenvector ',iev,
. ' iat=',iat
stop
401 print *,' End of file while expecting eigenvector ',iev
stop
100 format()
101 format(15x,3f12.6)
102 format('Eigenvector =',i6)
103 format('Frequency =',f13.6)
104 format(3e12.4)
end
C
C...............................................................
C
subroutine trans_1(MAT,MEV,nat,nev,nolab,nna,jat,idisp,
. lab,evr,evi,freq,qvec,wwr,wwi)
C Fourier transform eigenvectors,
C sum up over types and cartesian directions
C and take to square
C
implicit none
real*8 twopi
parameter (twopi=6.283185307)
integer MAT,MEV,nat,i,j,iev,nev,nolab,ilab,nrat,
. idisp(3,MAT),nna(MAT),jat(MAT,MAT)
real*8 evr(3,MAT,MEV),evi(3,MAT,MEV),freq(MEV),
. qvec(3),wwr(3,MAT),wwi(3,MAT),ww(3),arg,sinarg,cosarg
character*6 lab(MAT)
do iev=1,nev
write (6,204) iev, freq(iev)
204 format(' iev=',i4,' frequency = ',f10.4)
do ilab=1,nolab
do i=1,3
wwr(i,ilab)=0.d0
wwi(i,ilab)=0.d0
enddo
do j=1,nna(ilab) ! loop over atoms within the same type
nrat = jat(ilab,j)
arg=idisp(1,nrat)*qvec(1) +
+ idisp(2,nrat)*qvec(2) +
+ idisp(3,nrat)*qvec(3)
cosarg = cos(twopi*arg)
sinarg = sin(twopi*arg)
do i=1,3
wwr(i,ilab) = wwr(i,ilab) +
+ evr(i,nrat,iev)*cosarg - evi(i,nrat,iev)*sinarg
wwi(i,ilab) = wwi(i,ilab) +
+ evr(i,nrat,iev)*sinarg + evi(i,nrat,iev)*cosarg
enddo
enddo
do i=1,3
ww(i) = wwr(i,ilab)*wwr(i,ilab) + wwi(i,ilab)*wwi(i,ilab)
enddo
write (6,205) ilab, lab(ilab), ww
205 format (i6,a6,3f12.6)
enddo
enddo
return
end
C
C...............................................................
C
subroutine full_peaks(io1,nat,jat,label,qq,nodif,jtyp,nna,
. ivmin,ivmax,freq,evr,evi,wwsum)
C
C Takes square of eigenvectors,
C sums up over all atoms in type
C
implicit none
integer io1,iev,nodif,idif,jtyp(nat),nna(nat),
. nat,jat(nat,nat),nrat,ii,jj,ivmin,ivmax
double precision freq(ivmin:ivmax),
. evr(3,nat,ivmin:ivmax),evi(3,nat,ivmin:ivmax),
. wwsum(nodif,ivmin:ivmax),qq(3)
character*2 label(nat)
write (io1,201) qq
write (io1,202) (jtyp(idif),idif=1,nodif)
write (io1,203) (label(jat(idif,1)),idif=1,nodif)
write (io1,204)
C
do iev=ivmin,ivmax
do idif=1,nodif
wwsum(idif,iev)=0.d0
do jj=1,nna(idif) ! loop over atoms within the same type
nrat = jat(idif,jj)
do ii=1,3
wwsum(idif,iev) = wwsum(idif,iev) +
+ evr(ii,nrat,iev)**2 + evi(ii,nrat,iev)**2
enddo
enddo
enddo
write (io1,205) freq(iev),(wwsum(idif,iev),idif=1,nodif)
enddo
return
201 format('# Sum of squares of eigenvectors for Q =(',3f10.6,' )')
202 format('# Freq ',2x,100(2x,i4,4x))
203 format('# (cm-1) ',2x,100(4x,a2,4x))
204 format('#')
205 format (f10.3,100(f10.4))
end
C
C...............................................................
C
subroutine qres_peaks(io1,nat,jat,label,qq,nodif,jtyp,nna,coor,
. ivmin,ivmax,freq,evr,evi,qvec,wwsum)
C
C Fourier transforms eigenvectors,
C sums up over types and cartesian directions
C and takes to square
C
implicit none
integer io1,iev,nodif,idif,jtyp(nat),nna(nat),
. nat,jat(nat,nat),nrat,ii,jj,ivmin,ivmax
double precision freq(ivmin:ivmax),
. evr(3,nat,ivmin:ivmax),evi(3,nat,ivmin:ivmax),
. wwsum(nodif,ivmin:ivmax),
. wwr(3),wwi(3),coor(3,nat),
. qq(3),qvec(3),twopi,arg,sinarg,cosarg
character*2 label(nat)
parameter (twopi=6.283185307)
write (io1,201) qq,qvec
write (io1,202) (jtyp(idif),idif=1,nodif)
write (io1,203) (label(jat(idif,1)),idif=1,nodif)
write (io1,204)
C
do iev=ivmin,ivmax
do idif=1,nodif
wwsum(idif,iev)=0.d0
do ii=1,3
wwr(ii)=0.d0
wwi(ii)=0.d0
enddo
do jj=1,nna(idif) ! loop over atoms within the same type
nrat = jat(idif,jj)
arg=coor(1,nrat)*qvec(1) +
+ coor(2,nrat)*qvec(2) +
+ coor(3,nrat)*qvec(3)
cosarg = cos(twopi*arg)
sinarg = sin(twopi*arg)
do ii=1,3
wwr(ii) = wwr(ii) +
+ evr(ii,nrat,iev)*cosarg - evi(ii,nrat,iev)*sinarg
wwi(ii) = wwi(ii) +
+ evr(ii,nrat,iev)*sinarg + evi(ii,nrat,iev)*cosarg
enddo
enddo ! jj=1,nna(idif)
do ii=1,3
wwsum(idif,iev) = wwsum(idif,iev)+wwr(ii)**2+wwi(ii)**2
enddo
enddo ! do idif=1,nodif
write (io1,205) freq(iev),(wwsum(idif,iev),idif=1,nodif)
enddo
return
201 format('# Squared eigenvectors for Q =(',3f10.6,' )',/
. '# convoluted with qvec =(',3f10.6,' )')
202 format('# Freq ',2x,100(2x,i4,4x))
203 format('# (cm-1) ',2x,100(4x,a2,4x))
204 format('#')
205 format (f10.3,100(f10.4))
end
C
C...............................................................
C
subroutine full_smear(delta,npoints,io2,ivmin,ivmax,nat,jat,
. qq,nodif,jtyp,label,freq,wwsum,zz)
C Smear the peaks of phonon density of states
C (frequencies in "freq", squared eigenvectors in "wwsum")
C with a width parameter delta
C
implicit none
integer npoints,io2,iev,ivmin,ivmax,nodif,idif,ii,
. jtyp(nodif),nat,jat(nat,nat)
double precision freq(ivmin:ivmax),wwsum(nodif,ivmin:ivmax),
. zz(nodif),xx,delta,xmin,xmax,xstep,broad,qq(3)
character*2 label(nat)
external broad
C determine limits and mesh, to avoid many questions:
xmin = freq(ivmin) - (freq(ivmax)-freq(ivmin))*0.15
xmax = freq(ivmax) + (freq(ivmax)-freq(ivmin))*0.15
xstep=(xmax-xmin)/(npoints-1)
write (io2,201) qq,delta
write (io2,202) (jtyp(idif),idif=1,nodif)
write (io2,203) (label(jat(idif,1)),idif=1,nodif)
write (io2,204)
do ii=1,npoints
xx = xmin+(ii-1)*xstep
do idif=1,nodif
zz(idif)=0.d0
do iev=ivmin,ivmax
C --- discard acoustic mode of zero frequency in total DOS
if ( freq(iev).gt.1.d0 ) then
zz(idif)=zz(idif) +
+ broad(xx-freq(iev),delta)*wwsum(idif,iev)
endif
enddo
enddo
write (io2,205) xx,(zz(idif),idif=1,nodif)
enddo
return
201 format('# ',/
. '# Total phonon DOS for supercell Q = (',3f10.6,' )',/
. '# Sum of squares of eigenvectors, smeared with ',f8.4)
202 format('# Freq ',2x,100(2x,i4,4x))
203 format('# (cm-1) ',2x,100(4x,a2,4x))
204 format('#')
205 format (f10.4,100(f10.4))
end
C
C...............................................................
C
subroutine qres_smear(delta,npoints,io2,ivmin,ivmax,nat,jat,
. qq,nodif,jtyp,label,freq,wwsum,qvec,zz)
C Smear the peaks of phonon density of states
C (frequencies in "freq", squared eigenvectors in "wwsum")
C with a width parameter delta
C
implicit none
integer npoints,io2,iev,ivmin,ivmax,nodif,idif,ii,
. jtyp(nodif),nat,jat(nat,nat)
double precision freq(ivmin:ivmax),wwsum(nodif,ivmin:ivmax),
. zz(nodif),xx,delta,xmin,xmax,xstep,broad,qq(3),qvec(3)
character*2 label(nat)
external broad
C determine limits and mesh, to avoid many questions:
xmin = freq(ivmin) - (freq(ivmax)-freq(ivmin))*0.15
xmax = freq(ivmax) + (freq(ivmax)-freq(ivmin))*0.15
xstep=(xmax-xmin)/(npoints-1)
write (io2,201) qq,qvec,delta
write (io2,202) (jtyp(idif),idif=1,nodif)
write (io2,203) (label(jat(idif,1)),idif=1,nodif)
write (io2,204)
do ii=1,npoints
xx = xmin+(ii-1)*xstep
do idif=1,nodif
zz(idif)=0.d0
do iev=ivmin,ivmax
zz(idif)=zz(idif) +
+ broad(xx-freq(iev),delta)*wwsum(idif,iev)
enddo
enddo
write (io2,205) xx,(zz(idif),idif=1,nodif)
enddo
return
201 format('# ',/
. '# Squared eigenvectors for Q = (',3f10.6,' )',/
. '# convoluted with qvec = (',3f10.6,' )',/
. '# smeared with ',f8.4 )
202 format('# Freq ',2x,100(2x,i4,4x))
203 format('# (cm-1) ',2x,100(4x,a2,4x))
204 format('#')
205 format (f10.4,100(f10.4))
end
C
C...............................................................
C
function broad(xx,dd)
implicit none
double precision broad,xx,dd
C Explicitly defined broadening function:
C e.g. Lorentz:
broad=dd/3.1415926536/(xx*xx+dd*dd)
return
end
This directory contains routines written by Andrei Postnikov
([EMAIL PROTECTED])
for post-processing the results of siesta and/or
(SIESTA HOME)/Util/Vibra/Vibra/vibrator,
with the aim to vizualize, or better understand, the results
of lattice dynamics simulations.
Standalone programs followed by detailed explanations:
1. bonds.f - makes a table and a plot of bond lengths
between atoms of given types in supercell.
Reads only XV file, not the "vibrator" results.
2. phdos.f - phonon density of states over supercell
3. zerofc.f - zero out force constants beyond a given length
--- bonds.f ----------------------------------------------------------
- analyzis of interatomic distances between atoms of given two types
(which of course can be both the same), reading the atom coordinates
from the XV file and applying translations. The input file may
have an arbitrary name (which is asked for), but it must have
the structure of the XV file. One can change the atom type labels
in the input file, in order to select atoms within a type, the bonding
of which is of interest. I.e. it suffice to change the type
from, say, "3" to, say, "13" for those atoms whithin the type 3
whose bonds we want to analyze separately.
The script asks for (full) input file name, two atom types
to measure the distances between, and a bond length range
(min, max), in Ang. Finally, the input file name.
The program applies lattice translations along a fixed number
of displacements (fixed to 2 in the code; can be changed)
along each lattice vector, and produces a table of bond lengths,
which ndicates also running numbers of two involved atoms and the
lattice translation. Finally, the distribution of lengths is added
after this table as a continuous function, assuming an artificial
smearing of bond lengths. The discritization step of this
continuous function and the width parameter of the smearing are
fixed at the beginning of the code and can be changed.
--- phdos.f ----------------------------------------------------------
- calculation of phonon density of states after a supercell calculation
by siesta and then vibrator. Only single q-value is supported so far,
i.e. no q-space integration is performed. However, the supercell presumably
contain many atoms, which fall into different types, so that the local
density of states (according to different atom types) makes sense.
The density of states is in fact a collection of delta-peaks
at frequencies as calculated by vibrator, with the weights equal to
corresponding parts of the eigenvector squared.
The script offers a choice between the straightforward density of states
and the one involving convolution with a wave of given Q.
The idea behind is to project possibly complicated vibration modes
of a supercell onto lattice waves into a reference system of
the underlying single cell, and so recover something like
the phonon spectral density - see details in PRB 71, 115206 (2005).
---- Input files:
.XV - as from Siesta, but it can be manually modified to change
the attribution of atoms to different types. According to different
types found in this file, the local densities of states will be numbered.
For example, in order to split the atoms of type 3 into two different
species for the density of states, modify a to-be-split-off atom type
to, e.g., 13.
.vectors - as generated from vibrator.
The script asks for the calculation mode, -
full phonon DOS or q-projected, - and, in the latter case, for the
q vector to project upon.
---- Output files from a total DOS calculation:
.VST - discrete frequency values and corresponding weights,
.WST - the same, smeared and as continuous frequency function.
The frequence range and smearing parameter (10 inv.cm)
are set in the code:
program phdos
...
C Fix smearing parameter and No. of steps, to avoid many questions:
delta = 10.0
npoints = 2000
subroutine full_smear(...)
...
C determine limits and mesh, to avoid many questions:
xmin = freq(ivmin) - (freq(ivmax)-freq(ivmin))*0.15
xmax = freq(ivmax) + (freq(ivmax)-freq(ivmin))*0.15
xstep=(xmax-xmin)/(npoints-1)
Moreover, note that the smeared function supresses the low-frequency
acoustic mode:
subroutine full_smear(...)
C --- discard acoustic mode of zero frequency in total DOS
if ( freq(iev).gt.1.d0 ) then
zz(idif)=zz(idif) +
+ broad(xx-freq(iev),delta)*wwsum(idif,iev)
endif
---- Output files from a projected-q calculation:
.VSQ - discrete frequency values and corresponding weights,
.WSQ - the same, smeared and as continuous frequency function.
Differently from the total DOS option, a small frequencies-contribution
is not suppressed, because we might want to keep them in the spectral
function. The projection q value is indicated in the file heading,
otherwise the files structure as the same as for the total DOS.
--- zerofc.f ---------------------------------------------------------
discards the interactions beyond a given separation between atoms
in a .FC file. The idea is, to control the effect of different interactions
on the vibration spectrum.
---- Input files:
.XV
.FC - as from Siesta
---- Output file:
a replica of FC file, with certain elements set to zero,
and the meaning (atom pair; Cartesian displacement) appended
at the end of each line.
