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.