In the tutorial on Calculation of atomization energies, we calculated the atomization energy for H2 using the PBE functional. In this tutorial, we wish to estimate the uncertainty of the calculated atomization energy.

In the following script, an ensemble of 1000 atomization energies are calculated (non-selfconsistently) with an ensemble of 1000 exchange-correlation functionals distributed according to their probability (see article [Mor05b] for details).

from __future__ import print_function

import numpy as np
from ase.utils.bee import get_ensemble_energies
from ase.parallel import paropen as open
from gpaw import GPAW 

atom = GPAW('H.gpw', txt=None).get_atoms()
molecule = GPAW('H2.gpw', txt=None).get_atoms()
e1 = atom.get_potential_energy()
e2 = molecule.get_potential_energy()
ea = 2 * e1 - e2

fd = open('ensemble_energies.txt', 'w')
print('PBE:', ea, 'eV', file=fd)

e1i = get_ensemble_energies(atom)
e2i = get_ensemble_energies(molecule)
eai = 2 * e1i - e2i

n = len(eai)
ea0 = np.sum(eai) / n
sigma = (np.sum((eai - ea0)**2) / n)**0.5
print('Best fit:', ea0, '+-', sigma, 'eV', file=fd)

fd = open('ensemble.dat', 'w')
for e in eai:
    print(e, file=fd)

The script produces a text file named ensemble_energies.txt.

Ensemble of atomization energies for H2 :

[Mor05b]J. J. Mortensen, K. Kaasbjerg, S. L. Frederiksen, J. K. Nørskov, J. P. Sethna, and K. W. Jacobsen, Phys. Rev. Lett. 95, 216401 (2005)