  Although self-consistent multi-configuration methods have been used for
decades to address the description of atomic and molecular many-body systems,
only a few trials have been made in the context of nuclear structure. This work
aims at the development of such an approach to describe in a unified way
various types of correlations in nuclei, in a self-consistent manner where the
mean-field is improved as correlations are introduced. The goal is to reconcile
the usually set apart Shell-Model and Self-Consistent Mean-Field methods. This
approach is referred as "variational multiparticle-multihole configuration
mixing method". It is based on a double variational principle which yields a
set of two coupled equations that determine at the same time the expansion
coefficients of the many-body wave function and the single particle states. The
formalism is derived and discussed in a general context, starting from a
three-body Hamiltonian. Links to existing many-body techniques such as the
formalism of Green's functions are established. First applications are done
using the two-body D1S Gogny effective force. The numerical procedure is tested
on the $^{12}$C nucleus in order to study the convergence features of the
algorithm in different contexts. Ground state properties as well as
single-particle quantities are analyzed, and the description of the first $2^+$
state is examined. This study allows to validate our numerical algorithm and
leads to encouraging results. In order to test the method further, we will
realize in the second article of this series, a systematic description of more
nuclei and observables obtained by applying the newly-developed numerical
procedure with the same Gogny force. As raised in the present work,
applications of the variational multiparticle-multihole configuration mixing
method will however ultimately require the use of an extended and more
constrained Gogny force.
