Computing energy of Guassian basis set Actually, I'd like to compute decay rate using :

$$ \Gamma = 2\pi \ | \langle ab|V|cd \rangle|^2 \delta\ (\epsilon_f-\epsilon_i) $$

where a, b, c, d are the MOs and V is coulomb interaction, and $\epsilon_f$ and $\epsilon_i$ are final and initial energy of the system.

I am interested in computing this value regarding to a special atom (for example atom A) which means MOs should be expandes in terms of basis set:

$$ \Gamma = 2\pi \ \sum_{\mu, \nu,\lambda, \sigma \ \in A  }C_{\mu,a}\ C_{\nu,b}C_{\lambda,c}C_{\sigma,d}|\langle \mu\nu|V|\lambda \sigma \rangle|^2 \delta \ (\epsilon_{\mu}-\epsilon_{\nu}-\epsilon_{\lambda}+\epsilon_{\sigma}) $$

where $\mu, \nu, \lambda, \sigma$ are states in Gaussian basis sets, and $\epsilon$ are corresponding energies.

It seems a straightforward calculation, but I am completely confused because I don't know how I can calculate the energy of basis sets! I think the energy of the basis set is not meaningful.

Somebody, suggests using a numerical atomic orbital basis set to compute this kind of decay rate. (for numerical atomic orbital, the energy of states is meaningful) but I also don't know how to compute this kind of basis set!

Is there any suggestion to resolve this problem?
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.