Abstract. Using a density matrix approach to Gutzwiller method, we present a formalism to treat ab-initio multiband Tight-Binding Hamiltonians including local Coulomb interaction in a solid, like, for e.g., the degenerate Hubbard model. We first derive the main results of our method: starting from the density matrix of the non-interacting state, we build a multi-configurational variational wave function. The probabilities of atomic configurations are the variational parameters of the method. The kinetic energy contributions are renormalized whereas the interaction contributions are exactly calculated. A renormalization of effective on-site levels, in contrast to the usual one-band Gutzwiller approach, is derived. After minimization with respect to the variational parameters, the approximate ground state is obtained, providing the equilibrium properties of a material. Academic models will illustrate the key points of our approach. Finally, as this method is not restricted to parametrized Tight-Binding Hamiltonians, it can be performed from first principles level by the use of the so-called "Linearized Muffin Tin Orbitals" technique. To avoid double counting of the repulsion, one subtracts the average interaction, already taken into account in this density functional theory within local density approximation (DFT-LDA) based band structures method and one adds an interaction part "a la Hubbard". Our method can be seen as an improvement of the more popular LDA+U method as the density-density correlations are treated beyond a standard mean field approach. First application to Plutonium will be presented with peculiar attention to the equilibrium volume, and investigations for other densities will be discussed.