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 "bee.py"`

In ase/dft/bee.py change one line:

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

In gpaw/xc/bee.py 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
h2.center(vacuum=3)
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')
h.set_cell(cell)
h.center()
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)))
```