I'm trying to write a restricted Hartree Fock code in Fortran that reads in a file of zeroth-iteration 1 and 2 electron integrals (FCIDUMP format) and uses them to do the SCF procedure from an initial guess of a density matrix. The thing I'm struggling with is the 8-fold symmetry of the 2-electron integrals and how I can leverage it to construct my Coulomb and Exchange matrices efficiently.
FCIDUMP files use the symmetries of 1 and 2 electron integrals to save space, i.e. since $H_{\text{core},p,q}$ = $H_{\text{core},q,p}$, it will only print matrix elements of $H_{core}$ where $p \geq q$. Thus, reading in the matrix elements from FCIDUMP makes a lower triangular matrix which must be symmetrized along the diagonal to get the $H_{core}$ matrix. Easy enough:
h_core = e1ints
do q =1,norb
do p =1,q
h_core(p,q) = h_core(q,p) !lower triangular to symmetric
end do
end do
Where norb is the number of spatial orbitals (i.e. the dimension of the Fock matrix is norb*norb) and e1ints is the lower triangular matrix I've constructed from FCIDUMP. Of course, it's more memory-efficient to keep the matrix lower-diagonal, but since $H_{core}$ must be added to the Fock matrix eventually, I don't think there's any way around this.
The same printing methodology holds true for the 2 electron integrals, where the symmetry is $${\small (pq|rs) = (pq|sr) = (qp|sr) = (qp|rs) = (rs|pq) = (rs|qp) = (sr|qp) = (sr|pq)\tag{1}}$$
Where I'm getting tripped up is
- how much more time and memory symmetrizing a "lower triangular" 4-dimensional electron repulsion integral (ERI) tensor will take
- generally thinking about how to code the following across eight-fold symmetry instead of two-fold.
The general formula for the elements of the Coulomb and Exchange matrices, respectively, are:
$$J_{pq} = \sum_{r,s=1}^{\text{norb}} (pr|qs)P_{rs}\tag{2}$$
$$K_{pq} = \sum_{r,s=1}^{\text{norb}} (pr|sq)P_{rs}\tag{3}$$
where $P$ is the density matrix. Since of the $\text{norb}^4$ elements of the full ERI tensor, I only have $$\frac{\frac{\text{norb}(\text{norb}+1)}{2}\cdot(\frac{\text{norb}(\text{norb}+1)}{2}+1)}{2}\tag{4}$$ of them (for a $n\times n$ triangular matrix, there are $\frac{n(n+1)}{2}$ nonzero matrix elements).
How can I efficiently use this symmetry to construct the $J$ and $K$ matrices without making the ERI tensor huge by symmetrizing it? Since the tensor is not used explicitly in the construction of the Fock matrix, I don't think it makes sense to symmetrize it - I'm just at a loss for how to get around this. There are codes on GitHub that deal with FCIDUMP files, but after many hours of searching, I unfortunately haven't been able to find one that uses these files in the way I need to. With a fully symmetrized ERI tensor, I could do something like this, but there must be a more efficient way that avoids symmetrizing, right?
do i =1,norb
do j =1,norb
mat2int(i,j) = 0.0 !Coulomb + Exchange together
do k =1,norb
do l =1,norb
mat2int(i,j) = mat2int(i,j) + density_mat(k,l)*(2*e2ints_sym(i,j,k,l) - e2ints_sym(i,l,k,j))
end do
end do
end do
end do
