by using electron density as the basic variable. DFT transformed many sciences and technologies, but finding accurate, parameter-free approximations to its exchange-correlation energy functional that avoid self-interaction and capture strong electron correlation has remained difficult. One-particle density-matrix (1-DM) functional theories [@Gilbert1975] have one more degree of complexity. In them, the 1-DM is often represented by its eigenvalues, the occupation numbers, and the corresponding eigenvectors, the natural spin-orbitals (NSOs), e.g. [@Gritsenko2005; @Lathiotakis2008]. While avoiding the mean-field form of the 1-DM of DFT [@Kohn1965], the approximations to the exchange-correlation functional of the 1-DM have difficulties like those of the DFT approximations. Two-particle density-matrix (2-DM) functional theories, e.g. [@Garrod1964], are less reduced, an advantage. The ground-state energy is a known, explicit functional of the 2-DM in Coulombic systems. However, while necessary and sufficient conditions are known for the $N$-representability of the 1-DM, no such conditions exist for the 2-DM [@Coleman1963a]; reconstructing the $N$-particle wave function from the 2-DM is a QMA-hard problem [@Liu2007]. Nevertheless, major progress has been made towards necessary conditions for $N$-representability that can be systematically refined [@Zhao2004; @Mazziotti2012]. While not variational, the resulting calculations are almost as accurate as full configuration interaction (FCI) calculations [@Zhao2004; @Nakata2001; @Gidofalvi2008]. Their computational cost 