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
CXV 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
