The two-electron four-center integral of the homogeneous solid spherical harmonic Gaussian-type functions (GTF’s), r(2n+l)Y(lm)((r) over cap)exp(-alpha r(2)), has been evaluated analytically by decomposing it into a linear combination of two-center integrals through coincidence of centers. The two-electron two-center integrals are integrated analytically through the Fourier transformation convolution theorem. A compact integration formula is obtained for a general two-electron irregular solid spherical harmonic operator [4 pi/(2L + 1)]Y-1/2(LM)((r) over cap(12))/r(12)((L+1)). This formula is applied to evaluate two-center integrals of the Coulomb repulsion, the spin-other-orbit interaction and the spin-spin interaction by letting L=0, 1, and 2, respectively. The integration results are in terms of the L(n)’+(l’+1/2)(sigma R(2))R(l’)Y(l)’(m)’((R) over cap)exp(-sigma R(2)), of the relative nuclear coordinate plus one error-type F-function term. One-electron multicenter integrals have also been evaluated through Fourier transformation convolution theorem. These are the integrals of two and three center overlap, three-center multipole moments r(L)Y(LM)((r) over cap) (e.g., L=1 and 2 for dipole and quadrupole moment, respectively), and three-center irregular solid spherical harmonics Y-LM((r) over cap)/r((L+1)) (e.g., L=0, 1, and 2 for nuclear attraction, spin-orbit interaction, and electron spin-nuclear-spin interaction, respectively). All of the integral results are given explicitly in terms of the relative nuclear coordinates and the Gaussian exponential parameters. Explicit expressions of derivatives can be easily obtained to generate force constants, to optimize geometries and to optimize Gaussian exponential parameters.