How to calculate Zak phase from numerical wavefunctions with arbitray phase? In numerical calculations, an arbitrary gauge or phase attached to a wavefunction at a particular $\mathbf{k}$ places an obstacle in calculating the Berry connection $$\mathcal{A}(\mathbf{k})=\langle u_{n}(\mathbf{k})|\nabla_{\mathbf{k}}u_{n}(\mathbf{k})\rangle$$ because of the derivative. In 2D, for instance, using Berry curvature formulae can overcome this when calculating the Chern number.
However, for a 1D system, is there a way to calculate the (quantized, when there's appropriate symmetry) Zak phase $$\gamma =\oint_\mathrm{BZ}d\mathbf{k}\mathcal{A}(\mathbf{k})$$ without being affected by the arbitray gauge of wavefunctions?
 A: You can do this in the same way that Wannier charge centers are computed (see papers by Vanderbilt).
Suppose $\mathcal{C}$ is some closed path in $\mathbf{k}$-space (e.g. a 1D BZ). We'll define the Berry phase of the $n$th band along $\mathcal{C}$ as (note the imaginary unit, absent in your definition):
$$\gamma_n(\mathcal{C}) = \mathrm{i}\oint_{\mathcal{C}} \langle u_{n\mathbf{k}}|\boldsymbol{\nabla}_{\mathbf{k}}u_{n\mathbf{k}}\rangle.$$
The discrete formulation can be obtained by using e.g. forward differences and eliminating gauge invariances by cleverly taking logarithms of 1 + small terms (or by parallel transport reasoning). If we suppose the path is discretized into (not necessarily equidistant) $\mathbf{k}_i$ steps with $i=1,\ldots, N$ and $\mathbf{k}_{N+1} \equiv \mathbf{k}_1$, the end result is:
$$\gamma_n(\mathcal{C}) = \mathrm{Im}\log \prod_{i=1}^N \langle u_{n\mathbf{k_i}}|u_{n\mathbf{k}_{i+1}}\rangle$$
You can view this as the product (i.e. phase summation) of $N$ small rotations of the eigenvector's phase as it's transported along $\mathcal{C}$; the $\mathrm{Im}\log$-part merely picks out the phase.
If $\mathcal{C}$ is a non-contractible path in the BZ along a reciprocal lattice vector $\mathbf{G}$, it is desirable to enforce a periodic gauge, in which case one would take $u_{n\mathbf{k}_{N+1}}(\mathbf{r}) = u_{n\mathbf{k}_1}(\mathbf{r})\mathrm{e}^{-\mathrm{i}\mathbf{G}\cdot\mathbf{r}}$.
Z2Pack is a tool which implements this (to a much greater degree of generality...). That's also a good starting point for further reading.
