Notebooks · lib/hf/mdg_eri

Electron repulsion integrals on the kernels lib/cas generated

The production ERI member. It is lib/hf/md_eri's method, McMurchie-Davidson one shell quartet at a time with the sums over primitives inside two matrix products, with three changes:

Generated kernels. The R table of each primitive quartet comes from mdg_rfill and the E coefficients and the Hermite expansion of each primitive pair from mdg_pair, straight-line defs per angular-momentum class that lib/cas wrote (lib/hf/mdg_kernels, tools/cas/gen_integrals.viv). herm_rfill's loops recompute the cube offsets of every entry; the generated def has them as constants and keeps the levels in locals. The R values go into the matrix of the ket product through mdg_gather, generated per ket order with the ket offsets and signs as constants; mdg_gather_general reads them from a table and checks it.
Spherical rows in the store. For a spherical basis, each pair's E matrix is moved to spherical functions once, when the store is built: E_sph = (M_a (x) M_b)^T E_cart, M the c2s blocks of lib/hf/cart2sph. A quartet's block then comes out spherical, with fewer rows (5 for d, 7 for f, against 6 and 10), and no quartet is transformed afterwards; md_eri's mderi_sph multiplies every block by the two transforms.
The E matrices are kept, one record per pair (MdgE), where md_eri copies each one out of its flat store for every quartet it is in.

The store's directory and its per-pair headers are md_eri's (MDERI_*), and the shell groups of general contractions are md_eri's groups. Each pair also keeps the packed index of each of its rows, so writing a block is two loops over rows and columns (mdg_put), and the Hermite offset tables md_eri makes for every quartet are made once per build (mdg_htab).

Direct J and K (mdg_jk_screened, the production direct build of lib/hf/jk_direct) is lib/hf/md_direct's density-screened build on these quartets: the screening, the chunks and the eight contributions of each unique integral are md_direct's, and its header explains them. md_direct works over Cartesian functions, so a spherical run there moves each density to Cartesian functions and J and K back; here the quartets come out in the run's functions, so the density, J and K stay in them, with fewer functions per quartet and no transforms per build. The Schwarz factors must be those of the same functions (mdg_schwarz_basis).

Verification. md_eri stays as the general code and checks this module: mdg_screened against mderi_full, mdg_sph against mderi_sph and mdg_jk_screened against md_direct's mdjk_direct on the same molecules, with the tolerance of the registry, |v - median| / (1 + |median|) <= 1e-12, in lib/hf/tests/mdg_eri. The two do not run the same sums in the same order (the spherical transform moved into the store, the E entries multiplied in another order), so they agree to rounding, not bit for bit; the measured gap is in that test. The R kernels are bit for bit herm_rfill, and the tests here hold every class to it and every pair class to herm_efill.

Used with import hf.mdg_eri.

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

This program declares a record type (MdgE), and records run only in compiled C for now, so the Run buttons here are off. Open it in the workbench and use its C tab to compile it.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

[1]lines 57–821

from hf.opts import OPT_DIRECT_TOL, OPT_SCREEN_TOL
from hf.basis import (
    bas_atom_pos, bas_contract, bas_nao_cart, bas_nao_sph, bas_nc, bas_ncart, bas_nsph,
    bas_powers, bas_prim_exp, shell_atom, shell_cart0, shell_f0, shell_l, shell_nprim, shell_p0,
    shell_sph0,
)
from hf.eri_store import es_len8
from hf.hermite import herm_boys_fill, herm_esize, herm_etab, herm_nh, herm_rfill, herm_tuv
from hf.md_eri import (
    MDERI_BATCH, mderi_canon_pairs, MDERI_D_GA, MDERI_D_GB, MDERI_D_NA, MDERI_D_NB, MDERI_D_SA,
    MDERI_D_SB, MDERI_DIR_LEN, mderi_grid, mderi_groups, MDERI_HDR, mderi_hoff, MDERI_PI,
    MDERI_PP_LEN, mderi_qlist, MDERI_SCREEN_TOL, MDERI_TOP,
)
from hf.md_direct import MDJK_NCHUNK
from hf.mdg_kernels import mdg_gather, MDG_LMAX, mdg_pair, mdg_rfill, MDG_RMAX

#| One pair's E matrix, `[rows, primitive pairs * Hermite functions]`.
record MdgE { E: Arr }

# the directory entry md_eri leaves 0, which this store uses for where the
# pair's row index starts: per row `(i (i + 1)/2 + j, i, j)`, the first `-1`
# for a row whose pair is not canonical (`j > i` within one shell group)
MDG_D_IX = 7
MDG_IX_LEN = 3

#| The pair transform `(M_a (x) M_b)^T`, `[nsph_a nsph_b, ncart_a ncart_b]`,
#| from the c2s blocks of shells `sa` and `sb` in the block-diagonal `M`.
#@ 2543cdfa
def mdg_t2(M: Arr, shells: Arr, sa, sb) {
    na = bas_ncart(shell_l(shells, sa))
    nb = bas_ncart(shell_l(shells, sb))
    ma = bas_nsph(shell_l(shells, sa))
    mb = bas_nsph(shell_l(shells, sb))
    ca = shell_cart0(shells, sa)
    cb = shell_cart0(shells, sb)
    pa = shell_sph0(shells, sa)
    pb = shell_sph0(shells, sb)
    T = zeros(ma * mb, na * nb)
    for qa to ma {
        for qb to mb {
            for ia to na {
                for ib to nb {
                    T[qa * mb + qb, ia * nb + ib] = M[ca + ia, pa + qa] * M[cb + ib, pb + qb]
                }
            }
        }
    }
    return T
}

#| `mdg_pair` for a class with no generated kernel: the same entries from
#| `herm_etab` (lib/hf/hermite) in loops, into the flat `V`.
#@ 088d6645
def mdg_pair_general(V: Arr, la, lb, a, b, X, Y, Z) {
    PA = bas_powers(la)
    PB = bas_powers(lb)
    nb = len(PB)
    H = herm_tuv(la + lb)
    EX = herm_etab(la, lb, a, b, X)
    EY = herm_etab(la, lb, a, b, Y)

from hf.opts import OPT_DIRECT_TOL, OPT_SCREEN_TOL
from hf.basis import (
    bas_atom_pos, bas_contract, bas_nao_cart, bas_nao_sph, bas_nc, bas_ncart, bas_nsph,
    bas_powers, bas_prim_exp, shell_atom, shell_cart0, shell_f0, shell_l, shell_nprim, shell_p0,
    shell_sph0,
)
from hf.eri_store import es_len8
from hf.hermite import herm_boys_fill, herm_esize, herm_etab, herm_nh, herm_rfill, herm_tuv
from hf.md_eri import (
    MDERI_BATCH, mderi_canon_pairs, MDERI_D_GA, MDERI_D_GB, MDERI_D_NA, MDERI_D_NB, MDERI_D_SA,
    MDERI_D_SB, MDERI_DIR_LEN, mderi_grid, mderi_groups, MDERI_HDR, mderi_hoff, MDERI_PI,
    MDERI_PP_LEN, mderi_qlist, MDERI_SCREEN_TOL, MDERI_TOP,
)
from hf.md_direct import MDJK_NCHUNK
from hf.mdg_kernels import mdg_gather, MDG_LMAX, mdg_pair, mdg_rfill, MDG_RMAX

#| One pair's E matrix, `[rows, primitive pairs * Hermite functions]`.
record MdgE { E: Arr }

# the directory entry md_eri leaves 0, which this store uses for where the
# pair's row index starts: per row `(i (i + 1)/2 + j, i, j)`, the first `-1`
# for a row whose pair is not canonical (`j > i` within one shell group)
MDG_D_IX = 7
MDG_IX_LEN = 3

#| The pair transform `(M_a (x) M_b)^T`, `[nsph_a nsph_b, ncart_a ncart_b]`,
#| from the c2s blocks of shells `sa` and `sb` in the block-diagonal `M`.
#@ 2543cdfa
def mdg_t2(M: Arr, shells: Arr, sa, sb) {
    na = bas_ncart(shell_l(shells, sa))
    nb = bas_ncart(shell_l(shells, sb))
    ma = bas_nsph(shell_l(shells, sa))
    mb = bas_nsph(shell_l(shells, sb))
    ca = shell_cart0(shells, sa)
    cb = shell_cart0(shells, sb)
    pa = shell_sph0(shells, sa)
    pb = shell_sph0(shells, sb)
    T = zeros(ma * mb, na * nb)
    for qa to ma {
        for qb to mb {
            for ia to na {
                for ib to nb {
                    T[qa * mb + qb, ia * nb + ib] = M[ca + ia, pa + qa] * M[cb + ib, pb + qb]
                }
            }
        }
    }
    return T
}

#| `mdg_pair` for a class with no generated kernel: the same entries from
#| `herm_etab` (lib/hf/hermite) in loops, into the flat `V`.
#@ 088d6645
def mdg_pair_general(V: Arr, la, lb, a, b, X, Y, Z) {
    PA = bas_powers(la)
    PB = bas_powers(lb)
    nb = len(PB)
    H = herm_tuv(la + lb)
    EX = herm_etab(la, lb, a, b, X)
    EY = herm_etab(la, lb, a, b, Y)