  Several algorithms have been proposed to calculate the spatial entanglement
spectrum from high order Renyi entropies. In this work we present an
alternative approach for computing the entanglement spectrum with quantum Monte
Carlo for both continuum and lattice Hamiltonians. This method provides direct
access to the matrix elements of the spatially reduced density matrix and we
determine an estimator that can be used in variational Monte Carlo as well as
other Monte Carlo methods. The algorithm is based on using a generalization of
the Swap operator, which can be extended to calculate a general class of
density matrices that can include combinations of spin, space, particle and
even momentum coordinates. We demonstrate the method by applying it to the
Hydrogen and Nitrogen molecules and describe for the first time how the spatial
entanglement spectrum encodes a covalent bond that includes all the many body
correlations.
