The steps to compute the winding number are the following:
We discretize the BZ of the system. We need to take into account that k can take values that depend on the size of the system. Therefore, if L is the number of supercells of our system, K takes the following values = $\frac{2\pi m}{L} + \frac{\pi}{L}$ where $m = \{1, ..., L\}$. We have to make sure that the system is large enough, hence the discretization of BZ is small enough. We obtain a list $\{k_0, k_2, ..., k_{L-1}\}$. To be sure that we are not missing a phase, we can add a small value $\delta$ to the sequence. For the infinite system, we can choose any discretization.
For every $k_i$, we take the Hamitlonian at $k$ and $k_{i+1}$ and diagonalize it. We want to take all the eigenvectors that correspond to energies below zero (or above zero). Since the matrix is $2q\times 2q$ and the spectrum is symmetric, we will take $q$ eigenvectors. If we take ALL the eigenvectors, we should obtain a winding number of $0$. When we reach the last value of $k_i$, we need to impose periodic boundary conditions and therefore connect it to the first value of the chain.
For every $k_i$, we create a $q\times q$ tensor, $U^{(i)}$ such that: $$ U^{(i)}_{m,n} = \langle \psi_m(k_{i+1})\vert\psi_n (k_i) \rangle $$ Here $\psi_{m(n)}$ is the $m (n)$-th eigenvector of the Hamiltonian.
Finally, we will compute: $$ \nu = \frac{1}{\pi} \text{Im} \left[ \log \left(\prod_{i=1}^s \frac{\vert U^{(i)}\vert}{\sqrt{\vert U^{(i)} \vert \vert U^{(i)} \vert^*}}\right)\right] $$
For the long-range systems, it is not possible to obtain half-integer winding numbers if we close the loop. Nevertheless, without closing the loop we obtain a Gauge dependent quantity.
Here we compute the winding number as follows:
$$ w = Tr\left[Q_{BA}[X, Q_{AB}]\right] $$where $Tr$ is the trace per unit cell. Here, $Q = P_{+}-P_{-}$:
$$ P_{\pm} = \sum_{\pm} \mid \psi_{\pm} >< \psi_{\pm}\mid $$Then, $Q_{AB} = \Gamma_A Q \Gamma_B$, $Q_{BA} = \Gamma_B Q \Gamma_A$, where:
$$ \Gamma_A = \begin{pmatrix} 1 & 0 & 0 & 0 & .. \\ 0 & 0 & 0 & 0 & .. \\ 0 & 0 & 1 & 0 & .. \\ 0 & 0 & 0 & 0 & .. \\ \end{pmatrix}. $$$$ \Gamma_B = \begin{pmatrix} 0 & 0 & 0 & 0 & .. \\ 0 & 1 & 0 & 0 & .. \\ 0 & 0 & 0 & 0 & .. \\ 0 & 0 & 0 & 1 & .. \\ \end{pmatrix}. $$Computing $[X, Q_{AB}]$ is easy for an open system, but not for a system with periodic boundary conditions. Therefore, we use the steps described in this paper.
$$ [X, Q_{AB}] = \sum_{\lambda \neq 1} c_{\lambda} \lambda^X Q \lambda^{-X}, $$where we have a system in the domain $-N \leq x \leq N$, $\lambda$ are the solutions to $z^{2N+1} = 1$, $c_{\lambda} = \frac{\lambda^{N+1}}{1-\lambda}$
First, one has to find the basis where the chiral Hamiltonian has an off-diagonal form $H = \begin{pmatrix} 0 & h \\ h^\dagger &0\end{pmatrix}$ and the chiral operator has the form $\Gamma = \begin{pmatrix} \mathbb{1} & 0 \\ 0 & -\mathbb{1} \end{pmatrix}$. For this model, the latter can be done by transforming the Hamiltonian $H$ as follows
$$ H' = R_y H R_y^\dagger $$where $R_y = e^{i\sigma_y \frac{\pi}{2}}$. Moreover, we perform a rotation in order to rearrange the basis to the following form
$$ \chi_k' = \left(c_{k,0}, ..., c_{k, q-1}, c^\dagger_{-k,0}, ..., c^\dagger_{-k, q-1}\right)^T. $$Then, we take the upper-right block, $h$. The winding number will be given by:
$$ w = \int_0^{2\pi} \frac{1}{2\pi i} Tr\left[h^{-1}\delta_k h\right], $$or in the discrete version:
$$ w = \sum_k \frac{1}{L} Tr\left[h^{-1}\delta_k h\right]. $$We can consider the finite Hamiltonian, for which we have to choose the values of $k$ accordingly, and then take the numerical derivative of the Hamiltonian. The other opcion that is used for the paper is to take the infinite Hamiltonian and its the analytical derivative.
2.2 Thouless-Kohmoto-Nightingale-Nijs algorthm (TKNN)
In order to compute the winding number for the AAH edge states, we use the the Thouless-Kohmoto-Nightingale-Nijs (TKNN) algorithm, where they define the Berry curvature as follows
$$ \Omega_{k_x k_y} = i\sum_{\substack{E_{\alpha}<E_F \\ E_{\beta}>E_F}} \frac{\partial_{k_x} H_{\alpha\beta}\partial_{k_y} H_{\beta\alpha}-\partial_{\partial k_y} H_{\alpha\beta}\partial_{k_x} H_{\beta\alpha}}{(E_{\alpha}-E_{\beta})^2}, $$where $k_x$ and $k_y$ are quasi-momenta of the Brillouin zone. Recall that here, $k_x = k$ and $k_y = \phi$. We choose the value of the Fermi energy $E_F$ to be located inside the gap that we want to characterize. Note that
$$ \partial_{ k_x} H_{\alpha\beta} = \langle\psi_\alpha\vert \partial_{k_x} H \vert \psi_\beta\rangle $$The Chern number is defined as
$$ C= \frac{1}{2\pi i}\int_\text{BZ} \Omega_{k_x k_y}. $$We here discretize the integral as a Riemann sum $$ C = \frac{1}{2\pi i}\sum_{k_x,k_y} \Omega_{k_x k_y} \Delta k_x \Delta k_y. $$