  We present a self-consistent numerical approach to solve the Gutzwiller
variational problem for general multi-band models with arbitrary on-site
interaction. The proposed method generalizes and improves the procedure derived
by Deng et al., Phys. Rev. B. 79 075114 (2009), overcoming the restriction to
density-density interaction without increasing the complexity of the
computational algorithm. Our approach drastically reduces the problem of the
high-dimensional Gutzwiller minimization by mapping it to a minimization only
in the variational density matrix, in the spirit of the Levy and Lieb
formulation of DFT. For fixed density the Gutzwiller renormalization matrix is
determined as a fixpoint of a proper functional, whose evaluation only requires
ground-state calculations of matrices defined in the Gutzwiller variational
space. Furthermore, the proposed method is able to account for the symmetries
of the variational function in a controlled way, reducing the number of
variational parameters. After a detailed description of the method we present
calculations for multi-band Hubbard models with full (rotationally invariant)
Hund's rule on-site interaction. Our analysis shows that the numerical
algorithm is very efficient, stable and easy to implement. For these reasons
this method is particularly suitable for first principle studies -- e.g., in
combination with DFT -- of many complex real materials, where the full
intra-atomic interaction is important to obtain correct results.
