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?

• I'm not sure how the partitioning to atom $A$ will work. Would all AO basis functions $\mu, \nu, \lambda, \sigma$ have to be centered on $A$? Just the ones on one side of the integral? (The same is true for the answer by @jezzo) Apr 15, 2020 at 19:18
• In the summation, one just considers the terms which localized on atom A. Apr 15, 2020 at 19:33

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.