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

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,
            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 = calc.new(hund=False,  # no hund rule for molecules
                txt='H2.out')

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

with open('atomization.txt', 'w') as fd:
    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)

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.08 eV
  hydrogen molecule energy: -6.65 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.