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

Solution 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}

Solution 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