[Pw_forum] Re: doing QHA calculations (contd)

[Eduardo Ariel Menendez P]

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
!