[Pw_forum] Re: doing QHA calculations (contd)
Eduardo Ariel Menendez P
emenendez at macul.ciencias.uchile.cl
Mon Aug 21 23:34:48 CEST 2006
Hi Tuhina,
Maybe you can find useful this personal modification of fqha.f90. The
program is small, and it is easy to modify (if you know fortran).
Comments:
I think the format of the phonon DOS file is not exactly what
is read by fqha.x (see the commented line "! read (1,*) ndiv, emax, de".
There is no first line like "ndiv, emax, de" in the DOS file).
Moreover, the first frequency (read as nu(1)) is negative in the few DOS
files that I have generated, then the logarithm log(1.0d0-exp(-a3*nu(i)/T)
is NaN (as Eyvaz pointed). I just jumped and began the loop by
i=2, and I hope it will not harm, as DOS(0)=0.
Regards
Eduardo
!
! Copyright (C) 2001-2003 PWSCF group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
! Modified by E Menendez, 21 Ago 2006
! Calculate Free Energy F
! Given phonon DOS, calculate F at various temperatures
!
program fqha
!
implicit none
integer, parameter:: ndivx=10000
real(8) :: dos(ndivx),nu(ndivx),T,a1,a2,a3,Ftot,norm,F0
real(8) :: de, emax
integer :: i,ndiv
character(len=256) :: filename
!
a1=0.5d0/13.6058d0/8065.5d0
a2=8.617d-5/13.6058d0 ! Boltzmann constant kB in Ry/K
a3=1.0d0/8065.5d0/8.617d-5 ! Boltzmann constant kB in cm^-1/K
!
write (*,'(a,$)') ' file containing the dos >>>'
read(*,'(a)') filename
open(unit=1,file=filename,status='old')
! read (1,*) ndiv, emax, de
ndiv=0
i=1
120 continue
read(1,*,err=110,end=110) nu(i),dos(i)
! ndiv=ndiv+1
! write(*,*) nu(i),dos(i),i
i=i+1
! if (abs(nu(i)-de*(i-0.5d0)).gt.1.0d-6) stop ' wrong grid'
goto 120
110 continue
close(1)
ndiv=i-1
if (ndiv.gt.ndivx) stop ' ndivx too small'
write (*,*) 'ndiv=',ndiv
emax=nu(ndiv)
de=(emax-nu(1))/(ndiv-1)
write (*,'(a,$)') ' file for the Free energy >>>'
read(*,'(a)') filename
open(unit=1,file=filename,status='new')
!
1 continue
write (*,'(a,$)') ' temperature >>>'
read (*,*,err=10,end=10) T
Ftot=0.0d0
F0=0.0d0
norm=0.d0
do i=2,ndiv
! I found nu(1) in the dos file was negative,then I changed
! the initial 'i' from 1 to 2.
F0=F0+dos(i)*a1*nu(i)
if (T.gt.0.d0) Ftot=Ftot+dos(i)*a2*T*log(1.0d0-exp(-a3*nu(i)/T))
norm=norm+dos(i)
enddo
Ftot=(F0+Ftot)*de
norm=norm*de
write(1,*) T,Ftot
write(*,*) T,Ftot
write(*,*) norm,F0
!
go to 1
10 close(1)
!
stop
end program fqha
!
More information about the Pw_forum
mailing list