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_{core,p,q}$ = $H_{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 (ie 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 (pq|rs) = (pq|sr) = (qp|sr) = (qp|rs) = (rs|pq) = (rs|qp) = (sr|qp) = (sr|pq). Where I'm getting tripped up is a) how much more time and memory symmetrizing a "lower triangular" 4-dimensional electron repulsion integral (ERI) tensor will take and b) 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}^{norb} (pr|qs)P_{rs}$

and $K_{pq} = \sum_{r,s=1}^{norb} (pr|sq)P_{rs}$

where $P$ is the density matrix. Since of the $norb^4$ elements of the full ERI tensor, I only have $\frac{\frac{norb*(norb+1)}{2} *(\frac{norb*(norb+1)}{2}+1)}{2}$ of them (for a nxn 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.

  • 1
    $\begingroup$ Suggest: by content and purpose of your question, consider to move your question to the sibling site of mattermodeling.stackexchange.com. $\endgroup$
    – Buttonwood
    Jan 3 at 6:32
  • $\begingroup$ For a given two electron integral work out whih elements of the Hamiltonian it contributes to and then add in each of those contributions. You are then finished with that integral - no need to store the full rank 4 tensor, just (at most) the rank 2 Hamiltonian. $\endgroup$
    – Ian Bush
    Jan 3 at 7:39
  • $\begingroup$ BTW do you mean (pq|rs) in your expression for your Coulomb integral? $\endgroup$
    – Ian Bush
    Jan 3 at 15:03
  • 1
    $\begingroup$ There's an entire tag for FCIDUMP on Matter Modeling SE (see the link given by @Buttonwood) then click "tags" then "FCIDUMP". $\endgroup$ Jan 3 at 22:30


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.