I am trying to numerically compute the Berry Curvature for a generic quadratic Bosonic Hamiltonian of the form $$H = \sum_{ij} A_{ij} b_{i}^\dagger b_j + \frac{1}{2} \sum_{ij}\left( B_{ij} b_i b_j + \text{H.c.}\right).$$ After an appriate Fourier transform and Bogoliubov transformation, the Hamiltonian for the $n^{th}$ band can be written as
$$H_{n} = \sum_{\mathbf k} E(\mathbf k) \alpha_\mathbf k^\dagger \alpha_{\mathbf k}$$
for some bosonic operators
$$\alpha_{\mathbf k} := \sum_{j} \left[C_j(\mathbf k) b_j(\mathbf k) + D_j(\mathbf k) b_j^\dagger(\mathbf k) \right]$$
which satisfy $[\alpha_{\mathbf k},\alpha_\mathbf{k'}^\dagger] = \delta_{\mathbf{k}\mathbf{k'}}$, where $b_j(\mathbf k)$ is the $j^{th}$ bosonic annihilation operator in momentum space. A standard method for computing the Berry Curvature was introduced by Fukui et. al.
Fukui, Hatsugai, and Suzuki: Chern Numbers in Discretized Brillouin Zone: Efficient Method of Computing (Spin) Hall Conductances J. Phys. Soc. Jpn. 74, pp 1674-1677 (2005). https://arxiv.org/abs/cond-mat/0503172 .
Which describes computing the Berry Curvature in terms of so called $U(1)$ link variables
$$U_\mu({\mathbf{k}}) := \frac{\langle n(\mathbf{k})|n(\mathbf{k}+ \delta\mathbf{k}_\mu) \rangle}{|\langle n(\mathbf{k})|n(\mathbf{k}+ \delta\mathbf{k}_\mu) \rangle|}$$ where $\delta\mathbf{k}_\mu$ is a small vector that points in the $\mu^\text{th}$ direction in reciprocal space. The Berry Curvature is then approximated as $$F_{12}(\mathbf{k}) = \ln U_1(\mathbf{k}) U_2(\mathbf{k}+ \delta\mathbf{k}_1)U_1(\mathbf{k}+ \delta\mathbf{k}_2)^{-1}U_2(\mathbf{k})^{-1}.$$
In my context, I have computed the energies and corresponding Bogoliubov operators numerically, and can specify the energy eigenstates as $$|n(\mathbf{k})\rangle = \alpha^\dagger_{\mathbf k} |0\rangle$$ where $|0\rangle$ is the vacuum state. In this case however, it seems that $$\langle n(\mathbf{k})|n(\mathbf{k}+ \delta\mathbf{k}_\mu) \rangle = \langle 0|\alpha_{\mathbf k} \alpha^\dagger_{\mathbf{k}+ \delta\mathbf{k}_\mu} |0\rangle = \langle 0| \alpha^\dagger_{\mathbf{k}+ \delta\mathbf{k}_\mu} \alpha_{\mathbf k}|0\rangle = 0 $$ for any finite translation of the momentum vector by $\delta \mathbf k_\mu$. How can I proceed to compute the curvature numerically?