[Pw_forum] help me: how to calculate free energy

[hailin yu]

Dear Stefano B:
      Thank you for your advice.
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
the follow is my pragramme
                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>
>Reply-To: pw_forum at pwscf.org
>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
>>
>>_________________________________________________________________
>>Express yourself instantly with MSN Messenger! Download today it's  FREE! 
>>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