  We present a general ab initio method based on Wannier functions using the
covariant derivative for simulating the photocurrent in solids. The method is
widely applicable to charge/spin DC and AC photocurrent at any perturbation
levels in both semiconductors and metals for both linearly and circularly
polarized light. This is because the method is theoretically complete (within
the relaxation time approximation), that is to say, it includes all intraband,
interband and their cross terms. It is also free from the degeneracy issue,
i.e., applicable to arbitrary band structures with arbitrary numbers of
degenerate bands. We apply the method to various semiconductors and metals,
including GaAs, graphene-hBN heterostructure, monolayer WS2, a 2D ferroelectric
material - monolayer GeS, bilayer anti-ferromagnetic MnBi2Te4 and topological
Weyl semimetal RhSi, to simulate their charge and/or spin, DC and/or AC
photocurrent. Our theoretical results are in agreement with previous
theoretical works. Our numerical tests of GaAs, WS2 and GeS suggest setting the
degeneracy threshold in the conventional method as \hbar\Gamma^{2}, with
\Gamma^{2} the relaxation rate of the off-diagonal elements of the density
matrix between two states with close energies. We find that compared with the
conventional Wannier-function-based method using non-dgenerate perturbation
theory, the numerical errors of optical susceptibilities of bilayer
anti-ferromagnetic MnBi2Te4 with the PT symmetry can be reduced by 1-2 orders
of magnitude by our method for circularly polarized light. Our method provides
a universal computational tool for reliable and accurate predictions of
abundant weak-field photocurrent phenomena in disparate materials.
