# [Pw_forum] help me: how to calculate free energy

hailin yu yuhailin_79 at msn.com
Tue Aug 22 16:39:32 CEST 2006

Dear Stefano B:
I have writen a programme to compute the free energy according to
Equ(1)(PRB,57,10421(1998)). And I have checked the normalization of the DOS
and found the integral of the DOS  equal to 5.932386, agreement with 3*2 (2
atomic per unit cell, is unit cell or primitive cell?? ). I also have using
Equ(3)( Chem. Phys. Lett. 417,272-276(2005)) to compute free energy and
attainted the same results with using Equ(1)(PRB,57,10421(1998)). but our
results is still different with the previous work (PRB,57,10421(1998)). At
the range of 0~400 K, our results ( free energy) decrease nearly 6 eV,
while the previous work is about 0.26 eV. Where are wrong?

Thank you very much!
yuhailin
integer i
double precision hi,k,T,Z,F,dos
double precision aSn(197),aSndos(197)

open(1,file="aSndos.dat")
open(2,file="aSn-FreeEnergy.dat")

hi= 1.055  !
k=1.38	   ! the Boltzman's constant

do i=1,197
read(1,*)aSn(i),aSndos(i) ! the unit of omega is CM-1
dos=dos+aSndos(i)
enddo
print *,dos

do T=1.0,400.0,1.0
F=0.0
do i=2,197
Z=(hi*aSn(i)*0.029979*10.0)/(2*k*T) ! the value of (h*omega)/(2*K*T);
0.029979 is the coefficient of CM-1 ~THz
F=F+k*6.24151*(10.0E-5)*T*log(exp(z)-exp(-z))*aSndos(i) ! the sum of
k*T*ln[2*sinh(h*omega/2*k*T)]
!the unit of F is eV;
1J=6.24151*10E18 eV
enddo
write(2,*)T,F
enddo

close(1)
close(2)

end

>From: Stefano Baroni <baroni at sissa.it>
>To: pw_forum at pwscf.org
>Subject: Re: [Pw_forum] help me: how to calculate free energy
>Date: Tue, 22 Aug 2006 13:09:32 +0200
>
>
>On Aug 22, 2006, at 5:14 AM, hailin yu wrote:
>
>>Dear all:
>>       In recently, I have calculated the phonon dispersion and  phonon
>>density of states of alpha-Sn using PWSCF. And our results  are very
>>agreement with previously work(PRB 57,10421(1998), phonon  dispersion and
>>phonon density of states).
>
>very good
>
>>But when I calucalte the free energy, our results are very  different from
>>the previously work, and our results are too large.  So I want to konw
>>there have any programme to calculate the free  energy F(T) when attainted
>>the phonon density of states.
>
>The vibrational free energy can ben calculated from the vibrational  DOS
>using Eq. (1) of the above paper and replacing  the sum over  normal modes
>with an \int n(\omega) d\omega. This results in a simple  one-dimensional
>integral over a finite range. Are you sure you want  us to send you a
>10-lines piece of code doing just this?
>
>If you have already written down your own 10-lines code and still  find
>wrong results, you may want to:
>
>1) Check that your 10-lines code give correct results when used to
>integrate functions that you can integrate analitically. Beware that,  due
>to Van Hove singularities, the DOS may have to be sampled at a  rather
>large number of points in order to obtain sensible numerical  results.
>
>2) Check the normalization of the DOS. In order for the \int n (\omega)
>d\omega to give the same result as in Eq. (1), The integral  of the DOS
>must be equal to the number of normal modes per unit cell,  i.e. 3N, N
>being the number of atoms. With this normalization, Eq.  (1) will give the
>free energy per unit cell.
>
>Hoper this helps.
>
>Stefano B.
>
>
>>
>>Thanks a lot!
>>
>>yuhailin
>>
>>_________________________________________________________________
>>http://messenger.msn.click-url.com/go/onm00200471ave/direct/01/
>>
>>_______________________________________________
>>Pw_forum mailing list
>>Pw_forum at pwscf.org
>>http://www.democritos.it/mailman/listinfo/pw_forum
>
>---
>Stefano Baroni - SISSA  &  DEMOCRITOS National Simulation Center -  Trieste
>[+39] 040 3787 406 (tel) -528 (fax) / stefanobaroni (skype)
>
>Please, if possible, don't  send me MS Word or PowerPoint attachments
>Why? See:  http://www.gnu.org/philosophy/no-word-attachments.html
>
>
>

_________________________________________________________________
Check the weather nationwide with MSN Search: Try it now!
http://search.msn.com/results.aspx?q=weather&FORM=WLMTAG