ij \rangle$ goes over nearest neighbouring sites only, and using second-quantized notation where $a^\dagger_\alpha/a_\alpha$ creates/annihilates a fermion in a single-particle state $\alpha = \{i\sigma\}$ with $i$ the lattice label and $\sigma$ the spin projection ($\sigma = \pm\frac{1}{2}$), see [*e.g.*]{} [@bijbel; @fetter_walecka]. The 2D version is of particular interest because it is thought to exhibit high-temperature superconductivity [@anderson; @scalapino]. This simple model has a rich phase diagram and is therefore challenging to solve. This is due to the competition between the delocalizing first term (see Eq. \[hub\_ham\]), which allows for hopping from one site to its neighbours, and the second term, which locally repels electrons that are on the same site. A lot of numerical studies have been carried out for this system, using various methods such as quantum Monte Carlo [@qmc_2dhub], and the density matrix renormalization group (DMRG) [@xiang]. For an excellent overview of numerical studies performed on the 2D Hubbard model, we refer to the review by Scalapino [@scalapino]. In spite of all these numerical studies, surprisingly few reference results exist for the ground-state energy on finite size lattices. An earlier v2DM study focussed on the $4\times4$ lattice and was limited to the half-filled case [@nakata_2dhub]. In this paper 