How to get a two-electron expectation value for atoms in PySCF?

Question

I am interested in calculating the following expectation value for atoms using a converged Slater deteminant from in DFT

$$ \left< \Psi \left| \sum_{i,j=1}^N r_i^a r_j^b\right| \Psi \right> \; ,\tag{1}$$

with $a,b$ being some small integers (let's say $a=b=1$) and $i,j$ run over each individual electron up to the total number of electrons $N$, $r_i$ being the operator corresponding to the distance of electron $i$ from the nucleus. I want to use PySCF because the code should contain all tools I need.

Originally I assumed I can just use the converged electron density to calculate the expectation values, but I believe I would lose the exchange terms as opposed to explicitly writing out the Slater determinant on both sides.

I've figured that I will be needing one and two electron integrals. For example if $a=1$ and $b=0$ I think the matrix element reduces to

$$ \left< \Psi \left| \sum_{i=1}^N r_i\right| \Psi \right> = \frac{1}{N} \left( \sum_n \left< \phi_n|r|\phi_n \right> + \sum_{m,n} \sum_n \left< \phi_m|r|\phi_n \right> \right)\; ,\tag{2}$$

which should correspond to the standard one electron integral with the Slater determinant ansatz. I also believe if $b \neq 0$ I should be getting a two electron integral. Based on my previous question, I also have a code that can calculate an arbitrary one-electron integral between any molecular orbitals. I don't know yet how to do it for two electron integrals.

However, I believe calculating matrix elements of Slater determinants should not be a rare problem, so I think there must be a standard way in PySCF of doing this.

What is the best approach of calculating the matrix element I am looking for?

Answer 1

First of all, the "diagonal" contribution of Eq. (1), $$ \langle \Psi | \sum_i^N r_i^{a+b} | \Psi \rangle \tag{3} $$ is the expectation value of a one-body operator, so you already know how to deal with it. The rest of my answer deals with the remaining term $2\langle \Psi | \sum_{i<j}^N r_i^a r_j^b | \Psi \rangle$.

The expectation value of an arbitrary two-body operator $\hat{G} = \sum_{i<j} \hat{g}(i,j)$ over Slater determinants is given by the Slater-Condon rules (regardless of whether there is only atom in the system, and regardless of the form of the operator except that it's two-body)

$$ \langle \Psi | \hat{G} | \Psi \rangle = \frac{1}{2} \sum_m \sum_n \left( \langle \phi_m(\vec{r}_1) \phi_n(\vec{r}_2) | \hat{g} | \phi_m(\vec{r}_1) \phi_n(\vec{r}_2) \rangle - \langle \phi_m(\vec{r}_1) \phi_n(\vec{r}_2) | \hat{g} | \phi_n(\vec{r}_1) \phi_m(\vec{r}_2) \rangle \right) \tag{4} $$

In your case, you'll have two-electron integrals of the form $\langle \phi_m(\vec{r}_1) \phi_n(\vec{r}_2) | r_1^a r_2^b | \phi_m(\vec{r}_1) \phi_n(\vec{r}_2) \rangle$ and $\langle \phi_m(\vec{r}_1) \phi_n(\vec{r}_2) | r_1^a r_2^b | \phi_n(\vec{r}_1) \phi_m(\vec{r}_2) \rangle$. Since the integrands are a product of a function of $\vec{r}_1$ and a function of $\vec{r}_2$, the integrals are separable, e.g.

$$ \langle \phi_m(\vec{r}_1) \phi_n(\vec{r}_2) | r_1^a r_2^b | \phi_m(\vec{r}_1) \phi_n(\vec{r}_2) \rangle = \langle \phi_m(\vec{r}_1) | r_1^a | \phi_m(\vec{r}_1) \rangle \langle \phi_n(\vec{r}_2) | r_2^b | \phi_n(\vec{r}_2) \rangle \tag{5} $$

By now, the problem is completely reduced to the calculation of one-electron integrals.

For arbitrary two-electron operators, the separation Eq. (5) is generally not possible. However, it is always possible to write the two-electron integral as a sum of products of one-electron integrals, say via the multipole expansion, although the sum may be an infinite series such that it can only be evaluated approximately. Even so, this separation is nevertheless sometimes used to speed up the evaluation of the two-electron term in e.g. DFT and TDDFT calculations, as in program packages like Dmol3 and BDF.