\label{eq:H_SB}
    \mathcal{H}^{\rm SB}_i &=& \sum_{\alpha}  \frac{  \left(\Delta_{i, \alpha} \mathcal{M}_{\alpha} + \Delta^*_{i, \alpha} \mathcal{M}^{\dagger}_{\alpha} \right)}{\sqrt{n_{i\alpha} (1-n_{i\alpha})}} + U \,\mathcal{U} \nonumber \\
      & & + \sum_\alpha  \left( \Lambda_{i\alpha} + \mu_{i\alpha} \right) \mathcal{N}_\alpha + T \Theta_i,\end{aligned}$$ where $\Delta_{i, \alpha}$ encodes the local electron bonding $$\begin{aligned}
    \label{eq:Delta}
    \Delta_{i, \alpha} = \sum_j t_{ij} \mathcal{R}^*_{j,\alpha} {\rm Tr}(\rho^{\rm qp}  c^\dagger_{i\alpha} c^{\;}_{j \alpha} ).\end{aligned}$$ Explicit expressions for the various matrices $\mathcal{M}_\alpha$, $\mathcal{N}_\alpha$, $\mathcal{U}$ and $\Theta_\alpha$, as well as expression for the coefficient $\Lambda_{i, \alpha}$ are given in Appendix \[sec:GA\]. The Hamiltonian $\mathcal{H}^{\rm SB}$ essentially describes the competition between electron bonding favored by the first $\mathcal{M}$ term, and the Coulomb repulsion represented by the second $\mathcal{U}$ term. The Lagrangian multiplier $\mu_{i\alpha}$ appears in the third term acts as a chemical potential in order to satisfy the self-consistent Gutzwiller constraint.

The above two steps, namely solving the density matrix $\rho^{\rm qp}$ and variational matrices $\Phi_i$, have to be iterated until convergence is reached at every MD time-step. In particular, the Gutzwiller constraint $$\begin{aligned}
    {\rm Tr}(\Phi^\dagger_i \Phi^{\;}_i \hat{n}_{i, \sigma})  = {\rm Tr}(\rho^{\rm qp} \hat{n}_{i, \alpha}) = n_{i,\alpha},\end{aligned}$$ has to be satisfied. The bottleneck of the MD simulation here is solving the renormalized tight-binding Hamiltonian, which is required at every Gutzwiller iteration 