Numerical Berry curvature for bosons 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?
 A: There are multiple things to address here. First you would have reached the same conclusion whether the particles are bosons or fermions. In your reasoning as you defined $|n(\textbf{k})\rangle = \alpha^\dagger_\textbf{k} |0\rangle$, $\langle n(\textbf{k})|n(\textbf{k} + \delta \textbf{k})\rangle = 0$, whether $[\alpha_\textbf{k}, \alpha^\dagger_{\textbf{k}^\prime}] = \delta_{\textbf{k}\textbf{k}^\prime}$ or $\{\alpha_\textbf{k}, \alpha^\dagger_{\textbf{k}^\prime}\}=\delta_{\textbf{k}\textbf{k}^\prime}$. However I do think it's an interesting question to ask what's the meaning of the Berry curvature defined this way for a system of bosons. That's a different question though.
Before getting to the main confusion, I also wanna stress that you ran into a numerical problem. Mathematically you haven't shown that the $U(1)$ link is zero. You've shown that the numerator is zero, but you also have the denominator to be zero. So strictly speaking, one would need to work harder. However of course I understand, if you are running some sort of code this is bad.
Chiral Anomaly's answer address the confusion about the different Hilbert spaces. Here I'll try to be more explicit. The paper you cited is written in the context of condensed matter systems, that is periodic systems with some discreet transnational symmetries. For such periodic systems one always make use of Bloch's theorem and write the wavefunction as,
$$\psi^n_\textbf{k}(\textbf{r}) = e^{i\textbf{k} \cdot\textbf{r}} u^n_\textbf{k}(\textbf{r}),$$
where $u^n(\textbf{r})$ is a function that doesn't transform under the transnational symmetry of the system. It can be easily verified that
$$\langle \psi^n_\textbf{k}(\textbf{r})|\psi^{n^\prime}_{\textbf{k}^\prime}(\textbf{r})\rangle = \delta_{\textbf{k}\textbf{k}^\prime}\delta_{nn^\prime}.$$
Now in your question we have
$$\alpha^{n \dagger}_\textbf{k}|0\rangle = |\psi^n_\textbf{k} \rangle,$$
and with our choice of normalization above, this is consistent with your commutation relationships for $\alpha^n_{\textbf{k}}$, and $\alpha^{n\dagger}_{\textbf{k}}$, because remember that in general.
$$[\alpha^n_\textbf{k}, \alpha^{n\dagger}_{\textbf{k}^\prime}] = \langle \psi^n_\textbf{k}(\textbf{r})|\psi^{n^\prime}_{\textbf{k}^\prime}(\textbf{r})\rangle.$$
The wavefunction used in the paper you referred to, $|n(\textbf{k})\rangle$, is not the same. Rather with these periodic systems usually it's convenient to restrict oneself to only one unit cell. Here $|n(\textbf{k})\rangle$ would be defined within one unit cell of the system such that,
$$\langle r|n(\textbf{k})\rangle \equiv u_\textbf{k}^n(\textbf{r}).$$
What I mean by within a unit cell of the system is that when we see expressions like $\langle n(\textbf{k}) | n(\textbf{k}^\prime)\rangle $, the integration is only taken within one unit cell. Notice again that $| n(\textbf{k})\rangle$ is not what you have as $\alpha^n_{\textbf{k}}|0\rangle $.
Now the all important point is that even though $\langle \psi^n_\textbf{k}(\textbf{r})|\psi^{n}_{\textbf{k} + \delta \textbf{k}}(\textbf{r})\rangle = 0$, and hence the confusion, in general $\langle n_\textbf{k}(\textbf{r})|n_{\textbf{k} + \delta \textbf{k}}(\textbf{r})\rangle \neq 0$. I'm not fully aware how your code works for the calculation, but basically you need to extract the periodic part from your $\alpha^\dagger_{\textbf{k}}|0\rangle$, and then put that into the formula for the $U(1)$ link.
You might ask why use $|n(\textbf{k})\rangle$ in the definition of the $U(1)$ link and not $|\psi_{\textbf{k}}\rangle$. We already mentioned the awkward zero over zero situation. But also note that as one move through any path in the Brillouin zone, the Berry phase gained by $|\psi_{\textbf{k}}\rangle$ is the same as the berry phase gained by $|n(\textbf{k})\rangle$, and so it's usefull to define the $U(1)$ link in terms of $|n(\text{k})\rangle$.
A: If $|n(\mathbf{k})\rangle := \alpha^\dagger_\mathbf{k}|0\rangle$, then the Berry curvature is undefined: that's the dilemma that was highlighted in the question. The states $|n(\mathbf{k})\rangle$ typically used in Berry-curvature calculations cannot be the states $\alpha^\dagger_\mathbf{k}|0\rangle$ that were defined in the question.
The states $|n(\mathbf{k})\rangle$ typically used in Berry-curvature calculations (illustrated below) are eigenvectors of a matrix $H(\mathbf{k})$ that is distinct from the operator that was denoted $E(\mathbf{k})\alpha_\mathbf{k}^\dagger\alpha_\mathbf{k}$ in the question, but they are related: the coefficient $E(\mathbf{k})$ is an eigenvalue of the matrix $H(\mathbf{k})$ as well as being an eigenvalue of the operator $E(\mathbf{k})\alpha_\mathbf{k}^\dagger\alpha_\mathbf{k}$. However, the eigenvectors $|n(\mathbf{k})\rangle$ of $H(\mathbf{k})$ live in a different Hilbert space from the eigenstates of $E(\mathbf{k})\alpha_\mathbf{k}^\dagger\alpha_\mathbf{k}$.
To illustrate this, consider equation (6.45) in [1]:
$$
\newcommand{\da}{a^\dagger}
\newcommand{\db}{b^\dagger}
\newcommand{\la}{\langle}
\newcommand{\ra}{\rangle}
\newcommand{\pl}{\partial}
 H = \sum_\mathbf{k} \Big(\da(\mathbf{k})\ \db(\mathbf{k})\Big)
  H(\mathbf{k})\left(\begin{matrix} a(\mathbf{k})\\ b(\mathbf{k})\end{matrix}\right)
\tag{1}
$$
where $H(\mathbf{k})$ is a $2\times 2$ matrix. The Hamiltonian (1) can be derived from a special case of the Hamiltonian shown in the question, using appropriately-defined Fourier transforms.$^{[2]}$ This model has two bands.$^{[3]}$ In this case, the states $\alpha^\dagger_\mathbf{k}|0\ra$ described in the question are linear combinations of $\da(\mathbf{k})|0\ra$ and $\db(\mathbf{k})|0\ra$, where $|0\ra$ is the ground state of (1), and the coefficients in this linear combination are the components of one the eigenvectors of the matrix $H(\mathbf{k})$. Different eigenvectors of $H(\mathbf{k})$ correspond to different bands. The states $\alpha_\mathbf{k}^\dagger|0\ra$ are mutually orthogonal, as observed in the question, so we cannot define a Berry curvature for them. In contrast, the states $|n(\mathbf{k})\ra$ typically used in Berry-curvature calculations do not belong to the original Hilbert space at all. They belong to a different Hilbert space that is two-dimensional (in this two-band example), namely the Hilbert space in which the matrix $H(\mathbf{k})$ is defined. If we take the state $|n(\mathbf{k})\ra$ to be one of the eigenstates of this $2\times 2$ matrix, parameterized by $k$, then the states in this one-parameter family are typically not orthogonal to each other, so they can have a well-defined Berry curvature.
If this is the intent of the question, then the Berry curvature can be computed numerically using the approach in cond-mat/0503172, as reviewed in the question.

[1] http://www-personal.umich.edu/~sunkai/teaching/Fall_2014/Chapter6.pdf
[2] For the model considered in [1], a Bogoliubov transform is not needed, but it can be generalized so that a Bogoliubov transform is needed.
[3] The model can be generalized to an arbitrary number of bands, in which case $H(\mathbf{k})$ is generalized to a matrix of arbitrary size.
