Calculating atomization energies

1

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?

2

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.

3

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.
4

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).

5

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'})
6

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.

7

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.

8 
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

9

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.