Almost all functionals that are currently used in density matrix functional theory have been created by some a priori ansatz that generates approximations to the second-order reduced density matrix (2RDM). In this paper, a more consistent approach is used: we analyze the 2RDMs (in the natural orbital basis) of rather accurate multi-reference configuration interaction expansions for several small molecules (CH4, NH3, H2O, FH, and N2) and use the knowledge gained to generate new functionals. The analysis shows that a geminal-like structure is present in the 2RDMs, even though no geminal theory has been applied from the onset. It is also shown that the leading non-geminal dynamical correlation contributions are generated by a specific set of double excitations. The corresponding determinants give rise to non-JKL (non Coulomb/Exchange like) multipole-multipole dispersive attractive terms between geminals. Due to the proximity of the geminals, these dispersion terms are large and cannot be omitted, proving pure JKL functionals to be essentially deficient. A second correction emerges from the observation that the "normal" geminal-like exchange between geminals breaks down when one breaks multiple bonds. This problem can be fixed by doubling the exchange between bond broken geminals, effectively restoring the often physically correct high-spin configurations on the bond broken fragments. Both of these corrections have been added to the commonly used antisymmetrized product of strongly orthogonal geminals functional. The resulting non-JKL functional Extended Löwdin-Shull Dynamical-Multibond is capable of reproducing complete active space self-consistent field curves, in which one active orbital is used for each valence electron.