Dear Karolina:thanks for your interest; please try the attached code(a single
fortran source file, including all the subroutines called).Just compile it.It
is not separately documented, but you just run it and answer the
questions.Please ask me if something is not clear or goes wrong,best
regardsAndrei-- prof. Andrei Postnikov -- tel. +33-372749149 -- University of
Lorraine - Laboratoire de Chimie/Physique - A2MCICPM, 1 Bd Arago - BP 95823,
F-57078 Metz Cedex 03,
France------------------------------------------------------------------------
----- Karolina Milowska <karolina.milow...@gmail.com> a écrit :
>Dear All,
I would like to plot a phonon DOS resolved over particular atoms in my system,
in a similar fashion to the standard projected DOS. How may I extract this
information in practice?
I have found a presentation from 2010:
http://www.home.uni-osnabrueck.de/apostnik/Lectures/APostnikov-Phonons.pdf
which shows what I would like to get. On slide 21, there is an information
about utility - phdos. Where can I find it?
Kind regards,
Karolina
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 apost...@uos.de
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
C delta = 5.0
delta = 2.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)
102 format(14x,i6)
! 103 format('Frequency =',f13.6)
103 format(14x,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),zzsum
character*2 label(nat),lab_s
data lab_s/'++'/
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),lab_s
write (io2,204)
do ii=1,npoints
xx = xmin+(ii-1)*xstep
zzsum = 0.d0
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
zzsum = zzsum + zz(idif)
enddo
write (io2,205) xx,(zz(idif),idif=1,nodif),zzsum
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
