How to compute 2-electron integral for Hartree-Fock code?
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-&gt;agch', D4, Vprim, D4)

# 2-electron integral Option 2: transform indices one-by-one
V = Vprim.copy()
V = numpy.einsum('ba,bcde-&gt;acde', D2, V)    # i
V = numpy.einsum('abcd,be-&gt;aecd', V, D2)    # j
V = numpy.einsum('da,bcde-&gt;bcae', D2, V)    # k
V = numpy.einsum('abcd,de-&gt;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)  &amp; (\phi_1\phi_1|\phi_1\phi_2) &amp;  (\phi_1\phi_1|\phi_2\phi_1) &amp; (\phi_1\phi_1|\phi_2\phi_2)\\      (\phi_1\phi_2|\phi_1\phi_1)  &amp; (\phi_1\phi_2|\phi_1\phi_2) &amp;  (\phi_1\phi_2|\phi_2\phi_1) &amp; (\phi_1\phi_2|\phi_2\phi_2) \\      (\phi_2\phi_1|\phi_1\phi_1)  &amp; (\phi_2\phi_1|\phi_1\phi_2) &amp;  (\phi_2\phi_1|\phi_2\phi_1) &amp; (\phi_2\phi_1|\phi_2\phi_2) \\      (\phi_2\phi_2|\phi_1\phi_1)  &amp; (\phi_2\phi_2|\phi_1\phi_2) &amp;  (\phi_2\phi_2|\phi_2\phi_1) &amp; (\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} &amp; 0*\mathbf{D_2} \\
d_{12}*\mathbf{D_2} &amp; 0*\mathbf{D_2} \\
d_{13}*\mathbf{D_2} &amp; 0*\mathbf{D_2} \\
0*\mathbf{D_2} &amp; d_{21}*\mathbf{D_2} \\
0*\mathbf{D_2} &amp; d_{22}*\mathbf{D_2} \\
0*\mathbf{D_2} &amp; d_{23}*\mathbf{D_2} \\
\end{bmatrix} 
\end{equation}