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
  • 1
    $\begingroup$ This thesis discusses in Chapters 4 and 5 how to use only the nonredundant integrals in the density contraction. Not directly related to the 2e-integral part of your question, but the Fock matrix, $J$, and $K$ are also symmetric, so if you have a way to store them all lower/upper triangular, you can add them without having to explicitly symmetrize $H_\text{core}$. $\endgroup$
    – Tyberius
    Commented Jan 4, 2022 at 20:23
  • $\begingroup$ +1 and welcome to our new community! Thank you so much for contributing your question here and we hope to see much more of you in the future !!! Did you delete your similar question on Chem.SE? Also, you currently have an un-registered account and I'd recommend to register it so that you don't lose anything. $\endgroup$ Commented Jan 4, 2022 at 22:50
  • $\begingroup$ @NikeDattani I've now registered, trying to figure out how to delete the Chem.SE post since the registration link for Chem.SE timed out and now the site doesn't seem to think I'm the author. Thank you for the welcome! $\endgroup$
    – user4763
    Commented Jan 5, 2022 at 2:53
  • $\begingroup$ @Rob Your account still shows as unregistered. If you wound up making multiple accounts, you can get them merged by contacting the CMs here. Also, your attempt at the optimized code maybe better placed into an additional answer, rather than as an edit to the question. This will make it clearer to readers what the initial question was. $\endgroup$
    – Tyberius
    Commented Jan 5, 2022 at 20:44

1 Answer 1


AFAIK the standard way is to obtain the memory efficiency at the expense of some code repetition in the innermost loop. A preliminary example goes like:

  1. Construct a nested loop that only loops over the non-redundant (i,j,k,l) quartets. Say by
do i =1,norb
    do j =1,i
        do k =1,i
            do l =1,k
                if(k==i .and. l>j) cycle
  1. In the innermost loop, you can use the non-symmetrized integral e2ints(i,j,k,l), either by reading it from disk on-the-fly or by fetching it from an array in the memory. Enumerate all possible ways that e2ints(i,j,k,l) (and its symmetry-equivalent analogues, e.g. e2ints(i,j,l,k), e2ints(k,l,i,j) etc.) participates in the Fock building. You should come up with a maximum of 8 occurrences of this integral in the Coulomb term, and another 8 occurrences in the exchange term. When i==j or k==l or (i==k .and. j==l) etc., the number of terms are even fewer. You can either write a select case clause to deal with all the possible cases, or write a general case code that assumes (i,j,k,l) are mutually different, and scale the density_mat(*,*)*e2ints(i,j,k,l) contributions by a symmetrization factor (e.g. density_mat(k,l)*e2ints(i,j,k,l) should be scaled by 0.5 when k and l are equal) to avoid over-counting. The former approach has fewer FLOPs but requires 100~200 lines of additional code, while the latter approach has more FLOPs but gives a more elegant code.

The above approach is already essentially optimal in terms of memory. It can then be optimized for computational time, as follows:

  1. Some additions can be saved by taking advantage of the transpose symmetry of density_mat. For example mat2int(i,j) = mat2int(i,j) + density_mat(k,l)*e2ints(i,j,k,l); mat2int(i,j) = mat2int(i,j) + density_mat(l,k)*e2ints(i,j,k,l) can be reduced to mat2int(i,j) = mat2int(i,j) + 2.d0*density_mat(k,l)*e2ints(i,j,k,l).
  2. Similarly, some computation can be saved by taking advantage of the transpose symmetry of mat2int, but this is a more subtle issue. The trick is to: (1) whenever a term is being added to mat2int(i,j) and an equivalent term is being added to mat2int(j,i) (where i is not equal to j), only perform the former calculation; (2) whenever a term is being added to mat2int(i,i), scale the term by 0.5; (3) after all (i,j,k,l) quartets have been processed, add one line mat2int = mat2int + transpose(mat2int).
  3. Skip the quartet (i,j,k,l) if the absolute value of e2ints(i,j,k,l) is sufficiently small.

The above tricks (plus OpenMP parallelization) should already give you a code that rivals the best codes in terms of efficiency (but note that the contraction of integrals with density matrix elements is usually not the rate-determining step; rather, the bottleneck is usually the calculation of the integrals and/or the I/O of the FCIDUMP file), at least for molecules that are small enough such that the size of the FCIDUMP file is still reasonable. However, expect that your code expands to ~100 lines, and expect that most of the lines are copy-pasted from a previous line and then edited (thus error-prone). There are also some other optimizations that I did not mention, because they only make a difference when the molecule is so large that the FCIDUMP file cannot be stored on a typical disk.

Good luck!

  • $\begingroup$ Thank you! For just your first tip on optimizing for memory but not yet for time complexity, I've edited the post and written out at the very bottom what I think the nested loop should look like based on your description. Would you mind telling me if I've got it right? $\endgroup$
    – user4763
    Commented Jan 5, 2022 at 2:49
  • $\begingroup$ @Rob You have to do similar things for the permutation between i and j, and also the permutation between ij and kl. What I mentioned in my post are examples, not an exhaustive list of all the operations you need to do. That's why the final code can actually reach ~100 lines. Besides, in the loop you can only reference e2ints(i,j,k,l), because since you have not symmetrized e2ints(i,j,k,l), e2ints(i,l,k,j) does not exist (btw you have to switch to another data structure, as a 4-d array does not allow you to exploit the memory savings). $\endgroup$
    – wzkchem5
    Commented Jan 6, 2022 at 9:23
  • $\begingroup$ For example, instead of mat2int(i,j) = mat2int(i,j) + density_mat(k,l)*(2*e2ints(i,j,k,l) - e2ints(i,l,k,j)), you write two lines: mat2int(i,j) = mat2int(i,j) + 2.d0*density_mat(k,l)*e2ints(i,j,k,l), and mat2int(i,l) = mat2int(i,l) - density_mat(k,j)*e2ints(i,j,k,l). Note how the indices are set so that the same computation is performed but only e2ints(i,j,k,l) (instead of its symmetry equivalent counterparts) is referenced. Then perform the laborious but elementary step of including all other symmetry-equivalent operations and scaling them by appropriate factors, and you're done. $\endgroup$
    – wzkchem5
    Commented Jan 6, 2022 at 9:27

You must log in to answer this question.