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

Question

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 1

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}

Answer 2

As was mentioned in the comments, this amounts to the same fundamental operation as the AO-to-MO transformation; they are both 4-index transformations.

(...) each contracted integral $(\mathbf{ab}|\mathbf{cd})$ is expressed as a sum of its component primitive integrals $[\mathbf{ab}|\mathbf{cd}]$ which, in turn, are computed individually, i.e.

$$ (\mathbf{ab}|\mathbf{cd}) = \sum_{i=1}^{K} \sum_{j=1}^{K} \sum_{k=1}^{K} \sum_{l=1}^{K} \mathrm{D}_{\mathrm{a}i} \mathrm{D}_{\mathrm{b}j} \mathrm{D}_{\mathrm{c}k} \mathrm{D}_{\mathrm{d}l} \left[ \mathbf{a}_{i} \mathbf{b}_{j} | \mathbf{c}_{k} \mathbf{d}_{l} \right] $$

For the naive implementation, assuming you have all the primitives already formed, this is the simplest approach. If you have an AO-to-MO transformation that works, replace the MO coefficients as an argument with the contraction coefficient matrix. To get the naive implementation working, it would probably be easiest to reshape your primitives into a 4-index tensor.

• Gill, Peter M. W. Molecular integrals over Gaussian basis functions. Adv. Quantum Chem. 1994, 25, 141--205