+ U \sum_i n_{i,\uparrow} n_{i,\downarrow} \nonumber \\
    + \frac{K_0}{2} \sum_i |\mathbf u_i|^2 + \frac{K_1}{2} \sum_{\langle ij \rangle} \left[\hat{\mathbf e}_{ij} \cdot (\mathbf u_j - \mathbf u_i) \right]^2 + \sum_i \frac{|\mathbf p_i|^2}{2m}. \quad\end{aligned}$$ Here $c^\dagger_{i, \alpha}$ is the creation operator of electron with spin $\alpha = \uparrow, \downarrow$ at site-$i$, $n_{i, \alpha} = c^\dagger_{i,\alpha} c^{\;}_{i, \alpha}$ is the electron number operator, $U$ is the Hubbard repulsion parameter, $K$ is a elastic constant, $\mathbf p_i$ is the momentum operator, $m$ is the mass of the atom, $\mathbf u_i$ denotes the displacement vector of the $i$-th site, i.e. $\mathbf r_i = \mathbf r^{(0)}_i + \mathbf u_i$, and $\hat{\mathbf e}_{ij}$ is a unit vector pointing from site-$i$ to $j$. We assume the following dependence of hopping integral on the displacements $$\begin{aligned}
    \label{eq:t_ij}
    t_{ij} = t(\mathbf u_i - \mathbf u_j) = t^{(0)}_{ij} \left[1 + g \, \hat{\mathbf e}_{ij} \cdot (\mathbf u_j - \mathbf u_i) \right],\end{aligned}$$ where $t^{(0)}_{ij}$ is the bare hopping constant, and $g$ is the electron-phonon coupling. The Hamiltonian (\[eq:H\]) with $t_{ij}$ given by Eq. (\[eq:t\_ij\]) is called the Peierls-Hubbard model [@mazumdar83; @hirsch83], which has served as the basic platform for studying the interplay between the Peierls instability and electron correlation [@tang88; @fehske92; @yuan02]. Here 