How to calculate "heat of formation" using Vasp?

Message

ferdows · #1 Post by **ferdows** » Mon Dec 14, 2009 2:37 am

Hi,

I am trying to calculate heat of formation using Vasp. As a test case I have chosen the molecule: C6H6.

This is how I am going:

# All calculations are done at zero temperature, not considering thermal contribution at the moment.

# PBE with ENCUT = 500 for all calculations.

# Obtained the total electronic energy for C6H6. It is:

e0(c6h6) = -75.985 eV

# Obtained the electronic energy for H atom from H2 (standard state of H):

e0(H) = e0(H2)/2 = -3.385 eV

# Standard state for C is graphite. Optmized the graphite structure with 4 carbon atom unit cell. Then got the electronic energy of carbon from graphite as:

e0(C) = e0(graphite cell)/4 = 9.210 eV.

# Then use the following equation for heat of formation of C6H6:

delH(C6H6) = e0(c6h6) - 6*e0(C) - 6*e0(H) = -9.47 kcal/mol

Problem is: from Gaussian calculation I'm getting

delH(C6H6) = 30.1 kcal/mol

There is a big difference.

Am I calculating the heat of formation from Vasp correctly? Am I missing something? Please let me know. Any help will be greatly appreciated.

Thank you.

alex · #2 Post by **alex** » Mon Dec 14, 2009 9:17 am

Hi,
some things to check:
- is your basis set converged/large enough to give good results?
- do you use the same functional?
- did you use the "good" pseudopotiantials with small r_c?
- have you got some symmetry issue in the c6h6 calculation?
- experiment is delta_H = +82.9 kJ/mol, so both are off by about 10 kcal/mol assuming you have a typo in the vasp calc.
- did you take electronic energy in Gaussian or did you do a frequency calculation with thermochemistry afterwards ...

cheers

alex

panda · #3 Post by **panda** » Mon Dec 14, 2009 4:54 pm

yes your choice of pseudopotential in the case of H matters, see manual. Also as alex said did you do opt+freq in Gaussian? Either way you are off a bit...

ferdows · #4 Post by **ferdows** » Mon Dec 14, 2009 10:33 pm

[quote author= +82.9 kJ/mol, so both are off by about 10 kcal/mol assuming you have a typo in the vasp calc.
- did you take electronic energy in Gaussian or did you do a frequency calculation with thermochemistry afterwards ...

cheers

alex[/quote]

alex · #5 Post by **alex** » Tue Dec 15, 2009 7:57 am

Hi,

about 1) basic thermodynamics: product - reactant ... check sone textbook on that (very strong recommondation!)
about 2) normally I'd guess that total energy differencies cover more than 90% of the reaction enthalpies/free energies. However, at high T vibrational & entropy contributions may have a significant impact on that ...
- symmetry may force your wavefunction in a not favoured state. In the very beginning of OUTCAR (grep for ymmetr) the point group is plotted. Go then to your final DOS and check if this is as expected.
typo: you have -9 kcal for vasp and +30 for gaussian ...
experiment is gas phase (source webbook.nist.gov)
sign mismatch: indicates you are making a _huge_ mistake somewhere!
what is your gaussian basis set?
the difference between well converged calculations (e.g. wrt. basis set size) of both programs should be neglegible!

Btw. Your graphite is wrong. I haven't seen it in the first place. E_tot has to be negtive ...

Good luck!

alex

ferdows · #6 Post by **ferdows** » Wed Dec 16, 2009 2:00 am

Thank you, Alex. It is very helpful.

No, graphite is not wrong, it is just a typo. The value for
e0(C) is -9.21 eV.

Thermal energy is not the issue at the moment. Everything is at zero temperature. I will worry about thermal energy later.

Symmetry is most likely not the issue. I have done calculations for several other molecules. And the discrepancy in the Gaussian and Vasp results are there for all the molecules.

Gaussian calculations are done with: B3LYP/6-311+G(3df,2p)

Only thing I can think of now is to add the zero-point energy term (due to vibration) in addition to electronic energy. And that makes quite a bit of difference. The value of the heat of formation for C6H6 becomes 6.993 kcal/mol from -9.467 kcal/mol when I add the zero-point energy.

Overall, form your reply, it looks to me my basic procedure (and equation) to calculate heat of formation using Vasp is OK. There can be something off in the details.

Is there any Vasp references/documentation available for heat of formation calculations for compunds? Or any good paper where people used Vasp to calculate heat of formation of compounds?

Gaussian has excellent documentations for their thermochemistry calculations but I could not find any such documentation for Vasp.

[ Edited Wed Dec 16 2009, 03:02AM ]

alex · #7 Post by **alex** » Wed Dec 16, 2009 9:12 am

Hi,

you are actually comparing Gaussian B3LYP calculations with standard DFT of VASP? And you really wonder?

alex

ferdows · #8 Post by **ferdows** » Wed Dec 16, 2009 10:59 pm

Sorry, I did not understand what did you mean by "standard DFT of Vasp". I used PBE functional in my Vasp calculations.

Vasp with PBE and Gaussian with B3LYP are supposed to give very differernt results? Could you please elaborate on this? I do not have much experience on the comparison between different XC functionals.

Thanks for your help.

alex · #9 Post by **alex** » Thu Dec 17, 2009 6:18 am

Yep. The term standard DFT often refers to DFT without exact exchange, like PW91 or PBE. B3LYP however contains exact exchange and functionals different from PBE. So your results are _not_ supposed to be the same. Try PBE in Gaussian and compare afterwards. If then everthing is the same, you might go on.

Cheers,

Alex

ferdows · #10 Post by **ferdows** » Thu Dec 17, 2009 6:51 pm

Thank you very much, Alex. It was very helpful. Much appreciated.

panda · #11 Post by **panda** » Wed Dec 30, 2009 10:38 pm

There is LD(U+A) for exchange and correlation in VASP, it would be interesting to try this, although I have not done so myself as I am a bit skeptical...Electron correlation is a big issue in some systems...

vasp · #12 Post by **vasp** » Wed Jan 13, 2010 5:21 pm

Hi,
The way how you calculated the energy for a single carbon atom is intriguing for me. I have recently calculated the energy of carbon atom and got the energy around -1.0 plus or minus .15. depending on the parameters.
Your calculated energy of 9.21 by dividing the energy of graphene unit cell by four is very different. May be we are calculating different energies? Please explain how dividing the formation energy of four carbon atoms by four gives you that of the individual atom?

forsdan · #13 Post by **forsdan** » Wed Jan 13, 2010 5:43 pm

@vasp:

If one uses the graphite as reference phase to build the molecules then -9.21 eV/atom is the energy for one carbon atom in that environment with respect to the non-polarized carbon atom. The energy -9.21 eV/atom is therefore the correct energy to use in the evaluation of the formation energy when the carbon atoms is taken from the graphite phase. Your calculations are probably for a free carbon atom. There is a significant difference.

Cheers,
/Dan

[ Edited Wed Jan 13 2010, 06:58PM ]