Bayesian Error Estimation for RPBE

Here is a trick for making the Bayesian Error Estimation (BEE) with the RPBE functional. Just edit the lines in ASE and GPAW codes by adding RPBE as an exception.

To find the needed files, run

find ./ -name ""

In ase/dft/ change one line:

class BEEFEnsemble:

            if self.xc in ['BEEF-vdW', 'BEEF', 'PBE', 'RPBE']: # add RPBE
                self.beef_type = 'beefvdw'

In gpaw/xc/ add two lines:

class BEEFEnsemble:
    """BEEF ensemble error estimation."""
    def __init__(self, calc):

        # determine functional and read parameters
        self.xc = self.calc.get_xc_functional()
        if self.xc == 'BEEF-vdW':
            self.bee_type = 1
        elif self.xc == 'RPBE': # catch the RPBE exchange functional
            self.bee_type = 1   # assign BEEF coefficients the RBPE

Below we use BEEF-vdW, RPBE, and PBE dimensionless density (n) with gradient (s) and apply BEEF coefficients (E₀, ΔEᵢ) to evaluate the BEE as the standard deviation for the ensemble total energies with the variable enhancement factor (F(s,θᵢ)).

from ase import Atoms
from ase.dft.bee import BEEFEnsemble
from ase.parallel import parprint
from gpaw import GPAW
import time

for xc in ['BEEF-vdW','RPBE','PBE']:
    start_time = time.time()

    h2 = Atoms('H2',[[0.,0.,0.],[0.,0.,0.741]]) #exp. bond length
    cell = h2.get_cell()

    calc = GPAW(xc=xc,txt='H2_{0}.txt'.format(xc))
    h2.calc = calc
    e_h2 = h2.get_potential_energy()
    ens = BEEFEnsemble(calc)
    de_h2 = ens.get_ensemble_energies()
    del h2, calc, ens

    h = Atoms('H')
    calc = GPAW(xc=xc,txt='H_{0}.txt'.format(xc), hund=True)
    h.calc = calc
    e_h = h.get_potential_energy()
    ens = BEEFEnsemble(calc)
    de_h = ens.get_ensemble_energies()
    del h, calc, ens

    E_bind = 2*e_h - e_h2
    dE_bind = 2*de_h[:] - de_h2[:]
    dE_bind = dE_bind.std()
    parpting('{0} functional'.format(xc))
    parprint('Time: {0} s'.format(round(time.time()-start_time,0)))
    parprint('E_bind: {0} eV'.format(round(E_bind,4)))
    parprint('Error bar {0} eV'.format(round(dE_bind,4)))