Real space techniques in the theory of the chemical bond may in fact be understood as means to obtain chemical information from reduced density matrices (RDMs). Much work has been devoted to examine the one particle density, but the renewed interest in electron correlation effects on bonding has progressively shifted research towards further order densities. Up to date, however, not much work has been carried out to determine higher order RDMs from standard electronic structure packages beyond the single-determinant level. We show in this work a possible efficient route to this end. We discuss how to compute RDMs from multi-determinant wavefunctions, how to construct from them nth order cumulant density matrices (n-CDMs), and how to make profit of the recursivity and extensivity of the latter to define a full hierarchy of bonding indices through exhaustive partitions of the physical space. This hierarchy includes a full set of one- to n-center quantities (multi-center bonding indices) which may be conveniently decomposed into one-electron components, each of them partnered with a one-electron function (natural adaptive orbital, NAdO). Illustrative results are obtained and discussed for a set of simple test molecules.