Calculating atomization energies

Assuming I want to compute the atomization energy of methane with PyTorch which has an atomization energy of -0.631 Hartree. I compute the total energy with

eng ='b3lyp/6-311G(2DF,P)')

and I use the reference energies for the atoms in Hartree

H: -0.500273
C: -37.846772

By subtracting these reference energies I get an atomization energy of -0.687 Hartree which is a large difference to the expected result.

What am I missing?

I’m not familiar with PyTorch. You’re doing both the molecular and the atomic calculations at the same level of theory? Do you need to do open-shell and adjust the spin multiplicity for your atomic calcs? looks like a good paper to consult.

I meant Psi4 for Python instead of PyTorch, sry for the typo.

Maybe using the same level of theory solves the problem. How would my input file look like, when I want to calculate the energy of a single H atom?

For the the input

H 0.0 0.0 0.0

I get the error

Fatal Error: RHF: RHF reference is only for singlets.

That’s where multiplicity and open-shell come in. In that paper they use UHF for all the atomic calcs (set reference uhf). And you’ll probably have to adjust multiplicity, too. One place to start is +1 for mult (molecule {\n0 2\n<coordinates>\n} for neutral doublet).

I still get the error. My molecule string is

mol_str = '\n0 2\nH 0 0 0\n'

which I feed into

geom = psi4.geometry(mol_str)


eng ='b3lyp/'6-311G(2DF,P) mol=geom)


psi4.set_options({'reference': 'uhf'})

You’ll want to paste the whole input surrounded by triple backticks. For one thing, doesn’t take a kwarg mol but rather molecule, and set_options had better be before the energy call.

Maybe I did not state it clearly enough before, but I use the python module psi4 inside a python script (import psi4) and I pasted the exact lines of code that I use and that work well for other molecules.

Of course I set up set_options before the energy call.

import psi4

hatom = psi4.geometry("""
0 2

psi4.set_options({'reference': 'uhf'})

eng ='b3lyp/6-311G(2DF,P)', molecule=hatom)


yields -0.502157332116662
you’ve got a misplaced quote in the energy call

Thanks, that worked.

I ran the above code lines in ipython and I think it didn’t work for me in ipython because my paths were mixed with an older version of psi4.