  We describe and implement a first-principles algorithm based on
maximally-localized Wannier functions for calculating the shift-current
response of piezoelectric crystals in the independent-particle approximation.
The proposed algorithm presents several advantages over existing ones,
including full gauge invariance, low computational cost, and a correct
treatment of the optical matrix elements with nonlocal pseudopotentials.
Band-truncation errors are avoided by a careful formulation of $k\cdot p$
perturbation theory within the subspace of wannierized bands. The needed
ingredients are the matrix elements of the Hamiltonian and of the position
operator in the Wannier basis, which are readily available at the end of the
wannierization step. If the off-diagonal matrix elements of the position
operator are discarded, our expressions reduce to the ones that have been used
in recent tight-binding calculations of the shift current. We find that this
`diagonal' approximation can introduce sizeable errors, highlighting the
importance of carefully embedding the tight-binding model in real space for an
accurate description of the charge transfer that gives rise to the shift
current.
