Many-body eigenstates beyond the gaussian approximation can be constructed in
terms of local integrals of motion (IOM), although their actual computation has
been until now a daunting task. We present a new practical computation of IOMS
based on displacement transformations. It represents a general and systematic
way to extend Hartree-Fock and configuration interaction theories to higher
order. Our method combines minimization of energy and energy variance of a
reference state with exact diagonalization. We show that our implementation is
able to perform ground state calculations with high precision for relatively
large systems. Since it keeps track of the IMO's forming a reference state, our
method is particularly efficient dealing with excited states, both in accuracy
and the number of different states that can be constructed