by $$\begin{aligned}
    \label{eq:f_e}
    \mathbf f^e_{i} = g \sum_{j, \alpha} t_{ij}^{(0)} \hat{\mathbf e}_{ij} \left[ {\rm Tr}(\rho_G \, c^\dagger_{i,\alpha}c^{\;}_{j\alpha} ) + \mbox{c.c.}\right],\end{aligned}$$ where $\rho_G = \mathcal{P}_G \, \rho^{\rm qp} \, \mathcal{P}_G$ is the Gutzwiller electron density matrix [@sandri13], $\rho^{\rm qp}$ is the density matrix of quasi-particles, and $\mathcal{P} = \prod_i \mathcal{P}_i$ is the Gutzwiller projector. As in Gutzwiller’s original approach, the local projectors $\mathcal{P}_i$ are introduced to reduce the on-site double-occupancy variationally [@gutzwiller63; @gutzwiller64; @gutzwiller65]. In the modern formulation of the Gutzwiller approximation (GA), these projectors can be conveniently parameterized using the so-called $\Phi$-matrix [@lanata12], $\mathcal{P}_i \equiv \sum_{\alpha, \beta} \Phi_{i, \alpha\beta} /(\Pi_{i, \beta})^{-1/2} |\alpha\rangle \langle \beta|$, where $\Pi_{i, \beta}$ is the occupation probability of local atomic state $|\beta\rangle$ for the uncorrelated state. For single-band Hubbard model, the $\Phi$-matrix is diagonal $$\begin{aligned}
    \Phi_i = \left[ \begin{array}{cccc}
    e_i & 0 & 0 & 0 \\
    0 & p_{i, \uparrow} & 0 & 0 \\
    0 & 0 & p_{i, \downarrow} & 0 \\
    0 & 0 & 0 & d_i
    \end{array} \right],\end{aligned}$$ and its elements represent the empty, singly occupied (spin up and down), and doubly occupied electron state of the $i$-th atom [@lanata15; @lanata17]. As discussed in Sec. \[sec:intro\], they can also be 