  In the present paper, we propose an efficient numerical scheme for Gutzwiller
method for multi-band Hubbard models with general onsite Coulomb interaction.
Following the basic idea of Deng et al. [Phys. Rev. B 79, 075114 (2009)] and
extensions by Lanata et al. [Phys. Rev. B 85, 035133 (2012)], the ground state
is variationally determined through optimizing the total energy with respect to
the variational single particle density matrix (n0), which is called "outer
loop". In the corresponding "inner loop" where n0 is fixed, the non-interacting
wave function and the parameters contained in the Gutzwiller projector are
determined by a two-step iterative approach. All derivatives of the
implementation process have been analytically derived, which allows us to apply
some advanced minimization or root-searching algorithms for both the inner and
outer loops leading to the highly efficient convergence. In addition, an atomic
diagonalization method taking the point group symmetry into account has been
developed for the customized design of the Gutzwiller projector, making it
convenient to explore many interesting orders at a lower cost of computation.
As benchmarks, several different types of correlated models have been studied
utilizing the proposed method, which are in perfect agreement with the previous
results by DMFT and multi-orbital slave-boson mean field method. Compared with
the linear mixing method, the newton method with analytical derivatives shows
much faster convergence for the inner loop. As for the outer loop, the
minimization using analytical derivatives also shows much better stability and
efficiency compared with that using numerical derivatives.
