Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1import numpy as np
3# NB! This module was ported from a 4 year old CamposASE2 module.
5"""Bayesian Error Estimation
7For details, see: "Bayesian Error Estimation in Density Functional
8Theory", J. J. Mortensen, K. Kaasbjerg, S. L. Frederiksen,
9J. K. Norskov, J. P. Sethna, K. W. Jacobsen, Phys. Rev. Lett. 95,
10216401 (2005)."""
12# T
13# cost(c) = cost0 + 0.5 * (c - c0) H (c - c0)
14#
16# Cost function minimum value:
17cost0 = 3.4660625596
19# Best fit parameters:
20c0 = np.array([1.000787451, 0.1926284063, 1.896191546])
22# Hessian:
23# H = np.array([[ 1.770035168e+03, -3.732470432e+02, -2.105836167e+02],
24# [-3.732470432e+02, 1.188857209e+02, 6.054102443e+01],
25# [-2.105836167e+02, 6.054102443e+01, 3.211200293e+01]])
26#
27# 0.5 * np * T = cost0 (np=3: number of parameters)
28T = cost0 * 2 / 3
31def make_ensemble(N=1000, rng=np.random):
32 M = np.array([(0.066, -0.812, 1.996),
33 (0.055, 0.206, 0.082),
34 (-0.034, 0.007, 0.004)])
35 alpha = rng.normal(0.0, 1.0, (N, 3))
36 return c0 + np.dot(alpha, M)
39def get_ensemble_energies(atoms, ensemble):
40 if hasattr(atoms, 'get_calculator'):
41 coefs = atoms.calc.get_ensemble_coefficients()
42 else:
43 coefs = atoms
44 return coefs[0] + np.dot(ensemble, coefs[1:])