quantum chemistry - How to compute 2-electron integral for Hartree-Fock code?

I'm working on writing my own code for $\ce{H2}$ in an STO-3G basis set using Hartree-Fock (HF), and I am currently stuck on how to construct the two-electron integral matrix. I know how to evaluate the primitive two-electron integrals, but I can't figure out how to transform that into the $4\times4$ matrix of possible values for this system.

For my previous integrals, I determined the sum of the primitives by constructing a $6\times6$ matrix of the primitive integrals, and did a transformation using a $6\times2$ matrix of the coefficient values from the fit to the Slater function, for example: \begin{align} \mathbf{D} &= \begin{bmatrix} d_{11} & 0 \\ d_{12} & 0 \\ d_{13} & 0 \\ 0 & d_{21} \\ 0 & d_{22} \\ 0 & d_{23} \end{bmatrix}\\ \mathbf{S} &= \mathbf{D}' * \mathbf{S}_{prim} * \mathbf{D} \end{align}

I tried doing a similar technique for the two-electron integrals but with a $12\times12$ matrix, but the values didn't come out right. It seems that I need to make a $4$-dimensional array of all the primitive integrals, but I don't know how to evaluate that in order to determine the sum and hence the final

$4\times4$ matrix.

I've been using Szabo and Ostlund's Modern Quantum Chemistry book as my main reference.

Answer

You can also write this down in a similar way as for the one-electron integral. You already had:

\begin{equation} \mathbf{S} = \mathbf{D_2}' * \mathbf{S}_{\rm prim} * \mathbf{D_2} \end{equation}

where $\mathbf{S}$ is $(2\times2)$, $\mathbf{D_2}$ is $(6\times2)$ and $\mathbf{S}_{\rm prim}$ is $(6\times6)$.

For the 2-electron integral we now need

\begin{equation} \mathbf{V} = \mathbf{D_4}' * \mathbf{V}_{\rm prim} * \mathbf{D_4} \end{equation}

where $\mathbf{V}$ is $(2\times2\times2\times2)$, $\mathbf{D_4}$ is $(6\times2\times6\times2)$ and $\mathbf{V}_{\rm prim}$ is $(6\times6\times6\times6)$.

You can get $\mathbf{D_4}$ from $\mathbf{D_2}$: \begin{equation} \mathbf{D_4} = \mathbf{D_2} * \mathbf{D_2} \end{equation} where you do not sum over any index in order to get the 4-index tensor.

The equation for $\mathbf{V}$ is a bit more tricky, since you need to pay attention over which indices you sum (which is actually not specified by the notation). I add an example on how you can do this with Python/Numpy with Einstein summations.

# 1-electron integral
S = numpy.einsum('ba,bc,cd', D2, Sprim, D2)

# 2-electron integral Option 1: form (6x2x6x2) tensor
D4 = numpy.tensordot(D2, D2, 0)
V = numpy.einsum('badc,bedf,egfh->agch', D4, Vprim, D4)

# 2-electron integral Option 2: transform indices one-by-one
V = Vprim.copy()

V = numpy.einsum('ba,bcde->acde', D2, V)    # i
V = numpy.einsum('abcd,be->aecd', V, D2)    # j
V = numpy.einsum('da,bcde->bcae', D2, V)    # k
V = numpy.einsum('abcd,de->abce', V, D2)    # l

You will need to give D2, Sprim and Vprim as input.

(If anyone knows how to write this down in a formal notation, please edit.)

You were asking how to get $\mathbf{V}$ as a 2-index $(4\times4)$ matrix. For this we can simply reshape the above 4-index tensors. So $\mathbf{V}_{\rm prim}$ becomes $(36\times36)$. Now this means $\mathbf{D_4} $ needs to be $(36\times4)$, and not $(12\times12)$. To these reshaped 2-index matrices we can apply again the same transformation equation from above. Here is the code for my example:

# 2-electron integral Option 3: reshape to 2-index matrices
Vprim = Vprim.reshape((Nprim**2, Nprim**2))

D4 = numpy.swapaxes(D4, 1, 2)    # making sure to combine the correct axes
D4 = D4.reshape((Nprim**2, Ncntr**2))
V = numpy.einsum('ba,bc,cd', D4, Vprim, D4)

Now $\mathbf V$ will be: \begin{bmatrix} (\phi_1\phi_1|\phi_1\phi_1) & (\phi_1\phi_1|\phi_1\phi_2) & (\phi_1\phi_1|\phi_2\phi_1) & (\phi_1\phi_1|\phi_2\phi_2)\\ (\phi_1\phi_2|\phi_1\phi_1) & (\phi_1\phi_2|\phi_1\phi_2) & (\phi_1\phi_2|\phi_2\phi_1) & (\phi_1\phi_2|\phi_2\phi_2) \\ (\phi_2\phi_1|\phi_1\phi_1) & (\phi_2\phi_1|\phi_1\phi_2) & (\phi_2\phi_1|\phi_2\phi_1) & (\phi_2\phi_1|\phi_2\phi_2) \\ (\phi_2\phi_2|\phi_1\phi_1) & (\phi_2\phi_2|\phi_1\phi_2) & (\phi_2\phi_2|\phi_2\phi_1) & (\phi_2\phi_2|\phi_2\phi_2) \end{bmatrix}

And here is the structure of $\mathbf{D_4}$: \begin{equation} \mathbf{D_4} = \mathbf{D_2} * \mathbf{D_2} = \begin{bmatrix} d_{11}*\mathbf{D_2} & 0*\mathbf{D_2} \\ d_{12}*\mathbf{D_2} & 0*\mathbf{D_2} \\ d_{13}*\mathbf{D_2} & 0*\mathbf{D_2} \\ 0*\mathbf{D_2} & d_{21}*\mathbf{D_2} \\ 0*\mathbf{D_2} & d_{22}*\mathbf{D_2} \\ 0*\mathbf{D_2} & d_{23}*\mathbf{D_2} \\ \end{bmatrix} \end{equation}