# Calculation of atomization energies¶

Warning: mainstream DFT is unable to describe correctly the electronic states of isolated atoms (especially of transition metals). See doi:10.1063/1.2723118. This usually manifests itself as SCF convergence problems. Please consult literature before reporting such problems on the mailing lists.

The following script will calculate the atomization energy of a hydrogen molecule:

```# web-page: atomization.txt

from ase import Atoms
from ase.parallel import paropen as open
from gpaw import GPAW, PW, FermiDirac

a = 10.  # Size of unit cell (Angstrom)
c = a / 2

# Hydrogen atom:
atom = Atoms('H',
positions=[(c, c, c)],
magmoms=[0],
cell=(a, a + 0.0001, a + 0.0002))  # Break cell symmetry

# gpaw calculator:
calc = GPAW(mode=PW(),
xc='PBE',
hund=True,
eigensolver='rmm-diis',  # This solver can parallelize over bands
occupations=FermiDirac(0.0, fixmagmom=True),
txt='H.out',
)
atom.calc = calc

e1 = atom.get_potential_energy()
calc.write('H.gpw')

# Hydrogen molecule:
d = 0.74  # Experimental bond length
molecule = Atoms('H2',
positions=([c - d / 2, c, c],
[c + d / 2, c, c]),
cell=(a, a, a))

calc.set(txt='H2.out')
calc.set(hund=False)  # No hund rule for molecules

molecule.calc = calc
e2 = molecule.get_potential_energy()
calc.write('H2.gpw')

fd = open('atomization.txt', 'w')
print(f'  hydrogen atom energy:     {e1:5.2f} eV', file=fd)
print(f'  hydrogen molecule energy: {e2:5.2f} eV', file=fd)
print(f'  atomization energy:       {2 * e1 - e2:5.2f} eV', file=fd)
fd.close()
```

First, an `Atoms` object containing one hydrogen atom with a magnetic moment of one, is created. Next, a GPAW calculator is created. The calculator will do a calculation using the PBE exchange-correlation functional, and output from the calculation will be written to a file `H.out`. The calculator is hooked on to the `atom` object, and the energy is calculated (the `e1` variable). Finally, the result of the calculation (wavefunctions, densities, …) is saved to a file.

The `molecule` object is defined, holding the hydrogen molecule at the experimental lattice constant. The calculator used for the atom calculation is used again for the molecule caluclation - only the filename for the output file needs to be changed to `H2.out`. We extract the energy into the `e2` variable.

If we run this script, we get the following result:

```  hydrogen atom energy:     -1.07 eV
hydrogen molecule energy: -6.64 eV
atomization energy:        4.50 eV
```

According to Blaha et al. 1, an all-electron calculation with PBE gives an atomization energy of 4.54 eV, which is in perfect agreement with our PAW result.

The energy of the spin polarized hydrogen atom is -1.09 eV. If we do the calculation for the atom with `magmom=0` and `hund=False`, then we get almost 0 eV. This number should converge to exactly zero for a very large cell and a very high grid-point density, because the energy of a non spin-polarized hydrogen atom is the reference energy.

1

S. Kurth, J. P. Perdew, and P. Blaha, Int. J. Quantum Chem. 75, 889 (1999)