# 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 = psi4.energy('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?

https://aip.scitation.org/doi/pdf/10.1063/1.481544 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 https://github.com/psi4/psi4/blob/master/psi4/src/psi4/libscf_solver/sad.cc#L175-L180 +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)
``````

then

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

using

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

You’ll want to paste the whole input surrounded by triple backticks. For one thing, `psi4.energy()` 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
H
""")

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

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

print(eng)
``````

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.