Computing energy of Guassian basis set
First, these Coulombic interaction integrals can be simplified to $ \langle \mu\nu|\frac{1}{r_{12}}|\lambda \sigma \rangle $, which is the common two electron integral in the AO basis (your atom-centered gaussian basis functions).
You are correct that energies are only well defined for eigenfunctions of the Hamiltonian, and thus basis set functions (the choice of which is arbitrary) do not have a well-defined energy. Therefore, the most straightforward workaround for this is to calculate the two electron integrals in the AO basis, forming the rank 4 tensor of two electron integrals, $g(\mu,\nu,\rho,\sigma) = \langle \mu\nu|\frac{1}{r_{12}}|\lambda \sigma \rangle $. Then, you must (approximately) solve the Schrödinger equation with some level of electronic structure theory (e.g. CCSD, MP2, etc.) to solve for the molecular orbital coefficient matrix, C. The a$^{th}$ column of the matrix C is the weightings of the AO basis (atom-centered gaussian basis functions) to give the a$^{th}$ molecular orbital with well-defined energy $\epsilon_a$.
Finally, you should rotate $g(\mu,\nu,\rho,\sigma)$ with the coefficient matrix C to yield $g(a,b,c,d)$:
$$g(a,b,c,d)=C_{\mu a}C_{\nu b}g(\mu,\nu,\rho,\sigma)C_{\rho c}C_{\sigma d}$$
Note that, in this form, this change of basis scales as $N^8$, where N is the number of basis functions used, but can be implemented efficiently (by doing one rotation at a time) to scale as $N^5$. See the Crawdad programming tutorial on MP2 (site seems to be down right now) for more discussion on tensor rotations. Although it was done for MP2, it works identically for you.