Two-electron integral algorithm

Question

I am writing a molecular integral code from scratch. Right now my code have the following structure:

loop over shellpair ab:
    loop over shellpair cd:
        loop over primitive ij in ab:
            calculate E_ij
            loop over primitive kl in cd:
                calculate E_kl
                calculate R_ij_kl
                calculate V_ab_cd_ij_kl
        contract V_ab_cd_ij_kl to V_ab_cd

The above pseudo-code have alot of implicit notation, so let me explain what it means before asking my question.

E_ij and E_kl is the expansions coefficients from the McMurchie Davidson scheme. R_ij_kl is the Hermite integrals also used in the Mcmurchie Davidson scheme. V_ab_cd_ij_kl is the cartesian integrals found by (note the equations below are not strictly correct):

$$V_{ab,cd,ij,kl}=\sum E_{kl}\sum E_{ij} R_{ij,kl}$$

And V_ab_cd is the contraction of the primitive to the basisfunction:

$$V_{ab,cd}=\sum_i\sum_j\sum_k\sum_l c_ic_jc_kc_lV_{ab,cd,ij,kl}$$

When calculating integrals of higher angular momentum (above $0$), I make sure to reuse common elements in E and R. I.e. when doing p-type integrals px, py and pz have the first element in common, and it is only calculated once in the recurrence relations of the McMurchie Davidson scheme.

For my implementation of the above pseudo-code even the (ss|ss) integrals are very slow. When doing benenze with 6-31G the (ss|ss) integrals uses 2.5 seconds, which is more than can be justified by the use of Python (accelerated with Numba) instead of a compiled language. I have thought about rearranging and saving more intermidiates in memory, i.e. doing something like:

 loop over shellpair ab:
    calculate E_ab (Here E_ab means all E_ij saved in memory)
    loop over shellpair cd:
        calculate E_cd (Here E_ab means all E_kl saved in memory)
        loop over primitive ij in ab:
            loop over primitive kl in cd:
                calculate R_ij_kl
                calculate V_ab_cd_ij_kl
        contract V_ab_cd_ij_kl to V_ab_cd

In the above by calculating E_ab and E_cd, recomputation of E_ij and E_kl is avoided. At this point I know that I lag some expertise in knowing what is worth it compared to memory usage.

My question is know, what is a proper way of implementing molecular integrals with regard of what intermediates to keep at which points. (I know thing like vertical recurrence relations exist, but these will for obvious reasons not speed on my (ss|ss) integrals).

Answer 1

There are many things to consider and improve your code. Many components in an ab initio electronic structure calculations need supervision to achieve a good performance. I am no expert on those calculations, but I will point you to some ideas and references. On the theoretical side, you could calculate E_ab for all shell pairs first. For more clarification on the order of operations and their scaling look at Figure 3 in Multi-electron integrals. The article is written by a well-recognized researcher in the field that wrote a monumental book where you can find all needed expressions.

On the computational side, I would say that maybe Numba is not able to optimize your code. Writing it with Numba is not a guarantee of speed. I just took a glimpse at the code and, without going deep in the code, I saw many if statements. Only a straightforward observation. As general advice, it is all about how data is moved and the flow of operations. Avoid allocations, try to rewrite as much as possible. Profile your code to find the bottlenecks. Try to write it as vector or matrix operations to be able to use the large registers of the processor. There are many tricks to do these operations fast, but sadly they are buried elsewhere.

I would like to end with a reference to the recent literature on the topic. I just came up with an article that I managed to find again now: Fast Evaluation of Two-Center Integrals over Gaussian Charge Distributions and Gaussian Orbitals with General Interaction Kernels.