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. Specifically for the second-order dc photocurrent, it includes all of the following contributions-shift current, gyration current, (magnetic) injection current, Berry curvature dipole, and other Fermi surface contributions, instead of only a part of them as in most previous ab initio methods. 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 simulate the charge/spin dc and ac photocurrent of various semiconductors and metals, including GaAs, graphene-hBN heterostructure, monolayer WS2, monolayer GeS, bilayer antiferromagnetic MnBi2Te4, and topological Weyl semimetal RhSi. 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 h<overline>F(2 ), with F(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 nondegenerate perturbation theory, the numerical errors of optical susceptibilities of bilayer antiferromagnetic MnBi2Te4 with the PT symmetry can be reduced by one to two 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.