  Based on ab initio software packages using nonorthogonal localized orbitals,
we develop a general scheme of calculating response functions. We test the
performance of this method by calculating nonlinear optical responses of
materials, like the shift current conductivity of monolayer WS2, and achieve
good agreement with previous calculations. This method bears many similarities
to Wannier interpolation, which requires a challenging optimization of Wannier
functions due to the conflicting requirements of orthogonality and
localization. Although computationally heavier compared to Wannier
interpolation, our procedure avoids the construction of Wannier functions and
thus enables automated high throughput calculations of linear and nonlinear
responses related to electrical, magnetic and optical material properties.
