I believe the answer to your question is summarized well by the documentation of Psi4, an open source quantum chemistry package which can do similar calculations. I am fairly certain this is the problem you're running into but I could be wrong.
It is important to know that PSI4, like any other quantum chemistry program, does not compute the usual enthalpies, entropies, or Gibbs free energies of formation provided by most reference books. Instead, quantum chemistry programs compute “absolute” thermodynamic properties relative to infinitely separated nuclei and electrons, not “formation” values relative to elements in their standard states. If you are computing thermodynamic differences, like a reaction enthalpy computed as the enthalpy of the products minus the enthalpy of the reactants, then these “absolute” enthalpies are perfectly valid and usable. However, they cannot be mixed and matched with enthalpies of formation from reference books, since the zero of energy is not the same. Additionally, the “thermal energies” reported in kcal/mol are the finite-temperature corrections to the electronic total energy, and not the overall thermal energies themselves. If in doubt, use the reported Total Energies in Hartree/particle.
Link to documentation
In other words, you're value may be perfectly reasonable once you use it to calculate a change in energy, but as an absolute energy both values are fairly meaningless because they both pick a reasonable, but arbitrary, zero of energy.
This is also why you're internal energy can be negative. This simply means it has a smaller internal energy than whatever the reference internal energy is. I didn't check your work, but it is possible to have a negative energy depending on the reference value.
Also, here is the very thorough documentation (with examples!) from Gaussian. They also note,
This is generally not the same as the output from Gaussian for a calculation on an isolated gas-phase atom. These values are referenced to the standard states of the elements. For example, the value for hydrogen atom is 1.01 kcal/mol. This is ($H_{H_2}^{\circ}(298K)- H_{H_2}^{\circ}(0K))/2$, and not ($H_{H}^{\circ}(298K)- H_{H}^{\circ}(0K))/2$ , which is what Gaussian calculates.
So you're going to have the same problem. Read through there and I think it tells you how to make the correction. I think you just need to inclue some zero-point energies.
