# Two-electron integral algorithm

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).

• Semi-related, have you thought about using a lookup table for the Boys function? Your solution is probably faster than the analytic one, but some experts once told me about 30% of the time in calculating (ss|ss) is spent there. May 26, 2018 at 14:54
• Thanks for your observation. I have thought about making my Boys function faster through look up tabels as you suggest. This is also the recommended method in Molecular Electronic-Structure Theory. I think that right now I am mostly concerned about figuring out how I want the broad structure of the code. May 26, 2018 at 18:06
• I agree. As far as I know, the way the Boys function is incorporated never changes; it’s just called as a function, as you do, so any performance improvements are localized. May 27, 2018 at 21:33
• @user1271772 Being slow is not the intention, but probably poking fun at the fact that being fast enough for practical computation is very difficult. Avoiding recomputation of intermediates is the primary form of reducing the number of FLOPs, since you can't change the scaling complexity. May 28, 2018 at 1:22
• @user1271772 In that case I apologize. That document is not actually helpful because it is too high-level; using an analogy, it's like comparing Python with C that includes assembly in some places. The starting place is usually rsc.anu.edu.au/~pgill/papers/045Review.pdf, and in particular the section on "drivers". This is difficult because your intermediate is stored in a fixed-dimension Fortran array meant to behave as a hash table, where the key is an integer formed from basis set and shell pair info. In Python, a dictionary avoids some of that. May 28, 2018 at 3:06

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.