6
$\begingroup$

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

$\endgroup$
  • 2
    $\begingroup$ 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. $\endgroup$ – pentavalentcarbon May 26 '18 at 14:54
  • 1
    $\begingroup$ 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. $\endgroup$ – Erik Kjellgren May 26 '18 at 18:06
  • 1
    $\begingroup$ 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. $\endgroup$ – pentavalentcarbon May 27 '18 at 21:33
  • 1
    $\begingroup$ @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. $\endgroup$ – pentavalentcarbon May 28 '18 at 1:22
  • 2
    $\begingroup$ @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. $\endgroup$ – pentavalentcarbon May 28 '18 at 3:06

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.