# 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.

• Welcome to chemistry.SE! If you have any questions about the policies of our community, please ‎visit the help center. I edited your post a little to improve the mathematical formatting (let me know if I messed anything up). In the future, if you have any questions on using mathjax, visit this page on the math SE. – Tyberius Jun 16 '17 at 3:42
• Welcome to Chemistry.SE! We have help pages on MathJax, too, no need to meander off to Mathematics. If you want to know more, please have a look here and here. We prefer to not use MathJax in the title field, see here for details. – Martin - マーチン Jun 16 '17 at 4:07
• I am not entirely certain, but transformation from primitive integrals to AO integrals boils down to a 4-index transformation just like the one from AO to MO?! In this case, you will need some sort of 4-index tensor representation, which is precisely what I used back when writing my first exercise SCF (we were provided with the AO integrals through a subroutine). – TAR86 Jun 16 '17 at 8:51
• Are the $d_{ij}$ the contraction coefficients 0.15432897, 0.53532814, 0.44463454? – pentavalentcarbon Jun 16 '17 at 18:11
• @pentavalentcarbon: yes – Feodoran Jun 16 '17 at 21:57

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}

• Here are some timings for the 3 options. The 2nd one is by far the fastest on my Core i7-3667U, probably because it requires the fewest memory options or the strides are optimal, but I didn't bother figuring out FLOPs/MOPs. – pentavalentcarbon Jun 20 '17 at 16:12
• numpy.einsum is not very efficient, since it is pure Python code. But one can easily replace these calls using dot and tensordot, which interface to BLAS. See here. – Feodoran Jun 20 '17 at 16:31

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.

• From yours and other peoples' comments, I definitely see now that I cannot avoid making a 4-dimensional tensor of the primitive gaussian integrals. The problem I am having is doing the above summation in a way that condenses all the primitives into a $4\times4$ matrix. This is what I'm trying to accomplish: – SKB0514 Jun 20 '17 at 0:45
• $\begin{bmatrix} (\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_1)\\ (\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_1\phi_2) \\ (\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_1) \\ (\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|\phi_2\phi_2) \end{bmatrix}$ – SKB0514 Jun 20 '17 at 0:53
• I extended my answer to get the $(4\times4)$ matrix. Your mistake was to assume you need a $(12\times12)$ transformation matrix, but it should be $(36\times4)$. – Feodoran Jun 20 '17 at 8:46