Deriving Hartree from fundamental atomic units leads to error? I am trying to find how many newton meters are in a Hartree by using the following definition in terms of other physical constants:
$$E_h = \frac{\hbar^2}{m_ea_0^2} $$
The values of the other physical constants in SI units (according to Wikipedia and the NIST reference on constants) are:
\begin{align}
 m_e &= 9.109\ 383\ 7015(28) \times 10^{-31}\ \mathrm{kg}\\
 \hbar = \frac{h}{2\pi}; h &= 6.626\ 070\ 15  \times 10^{-34}\ \mathrm{kg} \frac{\rm m^2}{\rm s} \\
a_0 &= 5.291\ 772\ 109\ 03(80)\times 10^{−11}\ \mathrm m 
\end{align}
I cannot reproduce the answer from the NIST database, which says that
$$E_h = 4.359\ 744\ 722\ 2071(85)\times10^{−18}\ \mathrm J $$
However, the energy I get is $4.359\ 744\ 722\ 223\ 2755 \times 10^{-18} ... $
Can anyone reproduce their answer from the above numbers? Or alternatively, is there an explanation as to why the energy seems slightly off from the calculated value? (I am using the mpmath library in python with high precision. The code used to calculate and print the answer is below)
mpmath.mp.dps = 75
m_e = mpmath.mpf("9.1093837015e-31")
a_0 = mpmath.mpf("5.29177210903e-11")
h = mpmath.mpf("6.62607015e-34")
hartrees = mpmath.power((h/(mpmath.mpf(2)*mpmath.mp.pi)),2)/(m_e*a_0**2)
mpmath.nprint(hartrees,50)

 A: You are comparing
\begin{align}
E_h &= 4.359\ 744\ 722\ 2071(85)\times10^{−18}\ \mathrm J \\
\\
\frac{\hbar^2}{m_e a_0^2} &= 4.359\ 744\ 722\ 2232\ 755 \times 10^{-18}\ \mathrm J
\end{align}
Your final digits $\cdots 755$ are superfluous, smaller than the uncertainty.  Your disagreement is
$$
\frac{(\cdots 2232) - (\cdots 2017)}{85}
=1.9\sigma
$$
which is a little high.  However, you’re using a different expression for $E_h$ than the NIST website.  If I use their expression, I get
$$
\alpha^2 m_e c^2 = 4.359\ 744\ 722\ 1987 \times10^{-18}\ \mathrm J
$$
which is off by $-1.0\sigma$.
The NIST expression is probably preferable to yours, because $\alpha$ and $m_e$ are independent of each other, while the definition of $a_0$ includes the electron mass.  So to the extent that the electron mass is uncertain, that uncertainty multiplies your expression three times (once from $m_e$ and twice from $a_0^2$), while it enters the NIST expression only once.
If you really want to understand what’s happening here, you probably want to read the ugly details of CODATA’s enormous least-squares fit of all the data on physical constants; see the links in this other constants-related answer.
Note that you can get these results using regular double-precision arithmetic, rather than using a multiple-precision library.  (Multiple-precision libraries have their uses; they also have their quirks.) The problem from an arithmetic standpoint is that, in a world where the Hartree energy is a constant of nature, the measured values of $\alpha$ and $m_e$ are inconsistent with each other at the one-ish-sigma level. Dealing with those kinds of inconsistencies is the purpose of CODATA's least-squares fit to the entirety of the constants literature.
