Relativistic 2e 3- and 4-center integrals#399
Open
kshitij-05 wants to merge 46 commits into
Open
Conversation
loriab
reviewed
Feb 9, 2026
loriab
left a comment
loriab
left a comment
Collaborator
There was a problem hiding this comment.
I just noticed this and saw some tweaks to propose. You might want to add the new class to INSTALL.md, too.
… unique am shell sets and phase change for this operator
…l differentiator when on MacOS
…t and more cleanup
…+ progress bar + sign fix - ShellQuartetSetPredicate: add braket-swap tiebreaker for bra_ket_coswappable operators (σpσpCoulombσpσp). When la+lb == lc+ld, use max(la,lb) <= lc to pick one canonical representative, reducing duplicate quartet generation. - Engine (engine.impl.h): update swap_braket logic for opop_coulomb_opop to match the new predicate tiebreaker. Add coupled-swap sign correction in the swap_braket branch (was missing — exposed by d-shell testing). - build_libint.cc: disable CSE (do_cse/condense_expr) for multi-component operators since their 16 components share no intermediates at the expression level. This eliminates the superlinear optimize_rr_out bottleneck (e.g., 8.8s → 71ms for (ss|ds) prerequisite DAG). - build_libint.cc: fix compilation when only LIBINT_INCLUDE_RKB_ERI is defined (without LIBINT_INCLUDE_ERI): extend #ifdef guards for build_TwoPRep_2b_2k, add forward declaration, move make_descr to detail namespace, use if constexpr for component descriptor construction. - buildtest.h: add CodeGenProgress spinner showing elapsed time, function count, and current task name on stderr during code generation. - int_am.cmake: fix typo in OPT_AM variable reference.
kshitij-05
force-pushed
the
kshitij/feature/2e_rkb_ints
branch
from
March 22, 2026 18:19
25265a6 to
41f1dad
Compare
kshitij-05
force-pushed
the
kshitij/feature/2e_rkb_ints
branch
from
April 1, 2026 16:26
025cdfb to
3194cad
Compare
Apply the same optimization as HRR: since differentiation of a Gaussian depends only on that shell's quanta (not spectators), generate code once per unique differentiated shell and pass spectator dims at the call site. Verified: 295,944 assertions pass (Coulomb + ERI derivs + RKB integrals).
handle_trivial_nodes() used default_dims() (hardcoded "1") before
adapt_dims_() provided correct runtime dims ("lowdim"/"highdim").
Pass localdims through optimize_rr_out → handle_trivial_nodes.
kshitij-05
force-pushed
the
kshitij/feature/2e_rkb_ints
branch
from
April 10, 2026 18:02
3194cad to
9b50b0b
Compare
…p files Add code-sharing overrides (generate_label, spfunction_call, adapt_dims_, generate_code) to CR_11_Coulombσpσp_11 and CR_11_σpσpCoulombσpσp_11. Shell quartets with the same quaternion component share a single function parameterized by highdim. Hand-emits element-wise loops to avoid TwoPRep particle-swap child deduplication issues in DAG-based codegen.
(μ ∂_a ν | 1/r12 | κ ∂_b λ) for a, b ∈ {x,y,z}. Needed for Gaunt LS
Fock where coulomb_opop's 4 folded outputs are rank-deficient. Only
2-fold bra↔ket symmetry (p1_p2_swappable, no within-side swap); uses
dedicated predicate la+lb <= lc+ld. Also: LIBINT2_SIMPLE_CORE_EVAL_CASE
macro to compactify Coulomb-family dispatch in engine.impl.h.
Scalar trace + 3 antisym + 5 sym-TL replaces raw 3*a+b dyadic. Bra<->ket swap: per-component sign flip on antisym (was: index transpose). p1_p2_swappable=true. Tests pass: 2,112,533 assertions.
12 components = 3 dipole directions × 4 Pauli quaternion (trace + 3 antisym), mirroring σpVσp's fold. Engine origin via set_params. Master integral list split via mpl::joint_view (boost list50 limit).
…rkb_ints # Conflicts: # export/tests/unit/test-1body.cc # include/libint2/engine.h # include/libint2/engine.impl.h # src/bin/libint/comp_1_σpVσp_1.h # src/bin/libint/oper.h
kshitij-05
changed the title
Contributor
There was a problem hiding this comment.
Pull request overview
This PR extends Libint’s code generator and C++ engine to support relativistic restricted kinetic balance (RKB) 2-electron 4-center and 3-center integrals (including Gaunt LS bilinear forms), and adds the one-body (σ·p) r (σ·p) (“oprop”) operator. It also introduces an optimized “DerivGaussV2” recurrence strategy and improves generator progress reporting.
Changes:
- Add new operator descriptors/types for RKB Coulomb-family two-body operators and the one-body
opropoperator, wiring them through the generator’s master lists and strategy selection. - Update code generation graph/DAG plumbing to propagate correct implicit dimensions into RR optimization, and add
DerivGaussV2to enable better code sharing for derivative-Gaussian expansions. - Extend the runtime
Engineto recognize the new operators, adjust braket canonicalization/sign conventions, and add unit tests + CMake configuration switches for RKB ERIs (4c/3c).
Reviewed changes
Copilot reviewed 9 out of 25 changed files in this pull request and generated 2 comments.
Show a summary per file
| File | Description |
|---|---|
| src/bin/libint/strategy.cc | Switch TwoPRep derivative RR strategy to DerivGaussV2; register new RKB operator master strategies. |
| src/bin/libint/rr.cc | Compute localdims earlier and pass dims into DAG optimization for correct trivial-node handling. |
| src/bin/libint/oper.h | Add descriptors for σpRσp, Coulombσpσp, σpσpCoulombσpσp, σpCoulombσp; rename quaternionic_index → quaternion_index for σpVσp. |
| src/bin/libint/master_rrs_list.h | Register new composite recurrence headers and typedefs, including DerivGaussV2. |
| src/bin/libint/master_ints_list.h | Add new integral-set types; split master integral list into two parts joined via joint_view. |
| src/bin/libint/gauss.cc | Make primitive labels filesystem-safe (case-insensitive FS) by adding an underscore marker for uncontracted. |
| src/bin/libint/dg.h | Thread ImplicitDimensions into optimize_rr_out() / handle_trivial_nodes(). |
| src/bin/libint/dg.cc | Implement dims-aware optimize_rr_out() and handle_trivial_nodes() behavior. |
| src/bin/libint/comp_deriv_gauss_v2.h | New optimized derivative-Gaussian recurrence relation with HRR-like code sharing. |
| src/bin/libint/comp_11_σpσpCoulombσpσp_11.h | New CR implementing 16-component (σ·p)(σ·p) (1/r12) (σ·p)(σ·p) decomposition and code emission. |
| src/bin/libint/comp_11_σpCoulombσp_11.h | New CR implementing 9-component SO(3) irrep decomposition for Gaunt LS “bilinear” operator. |
| src/bin/libint/comp_11_Coulombσpσp_11.h | New CR implementing 4-component quaternion decomposition for `(ab |
| src/bin/libint/comp_1_σpVσp_1.h | Align σpVσp descriptor API (quaternion_index) with naming changes. |
| src/bin/libint/comp_1_σpRσp_1.h | New CR implementing 12-component (σ·p) r_k (σ·p) fold. |
| src/bin/libint/buildtest.h | Add global codegen progress tracker and hook it into GenerateCode(). |
| src/bin/libint/build_libint.cc | Add RKB task wiring, component loops, canonicalization logic, and progress printing; template two-/three-center build helpers by operator. |
| include/libint2/engine.impl.h | Extend operator list, nargs=3 routing, braket canonicalization/sign rules, intrinsic deriv order handling, and core-eval dispatch for new operators. |
| include/libint2/engine.h | Add new Operator enum entries + traits (nopers, intrinsic_deriv_order) for RKB ops and oprop. |
| include/libint2/cxxapi.h | Allow building C++ API when only RKB ERIs are enabled as the 2-body class. |
| include/libint2/config.h.cmake.in | Add config macros for enabling RKB ERI and RKB ERI3 and their AM/OPT settings. |
| include/libint2.h | Update generated macro names (underscore normalization). |
| export/tests/unit/test-2body.cc | Add unit tests validating RKB Coulomb-family operators + 3-center paths and xx_xs alias behavior. |
| export/tests/unit/test-1body.cc | Add unit test validating oprop Hermiticity/sign conventions across bra↔ket. |
| CMakeLists.txt | Add CMake options for enabling RKB ERI / RKB ERI3 and their AM/optimization controls. |
| cmake/modules/int_am.cmake | Teach component generation to include RKB ERI and RKB ERI3 component sets; fix OPT_AM variable assignment. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
The (ss|ss) fast-path predicate in engine.impl.h used a single bitwise & between the deriv_order_ and intrinsic_deriv_order() boolean tests, which drops short-circuit evaluation. Use &&. (Copilot review.)
…lomb_op doc - INSTALL.md: document LIBINT2_ENABLE_RKB_ERI / RKB_ERI3 and their MAX_AM/OPT_AM options (Lori review). - CMakeLists.txt + config.h.cmake.in: state that RKB_ERI3/ERI3 MAX_AM set the density-fitting (auxiliary) basis AM on the single fitting center. - int_am.cmake: collapse the RKB ERI component-code token rkb_eri_ -> rkberi_ so the operator token is single, underscore-free (like eri/onebody/g12), avoiding ambiguity if codes are split on '_'. Public LIBINT2_*_RKB_ERI* options and the RKB_ERI class id are unchanged (Lori review). - engine.h: correct the Operator::op_coulomb_op docstring to describe the actual SO(3) irreducible 9-component layout, not a raw 3*a+b Cartesian dyadic (Copilot review).
operator_traits<op_coulomb_op> carried the same stale "raw 3*a+b Cartesian dyadic" description that the Operator enum doc had; correct it to the actual SO(3) irreducible 9-component layout (Scalar / Antisym / SymTL). Completes the Copilot docstring fix, which had only updated the enum copy.
- engine.h: default_braket(coulomb_op) now returns BraKet::xs_xx. coulomb_op has no 4-center kernel, so the rank-based xx_xx default silently dispatched to the plain-Coulomb `default` task and returned wrong integrals with no error in release builds. (Review finding #1.) - comp_11_σpσpCoulombσpσp_11.h / comp_11_Coulombσpσp_11.h: add a `default` branch to generate_code()'s child-count switch that throws instead of emitting a wrong/empty `target[hsi] = ...` expression if num_children() is ever an unexpected value (e.g. future child dedup). (Review finding #4.)
…, coulomb_op default-braket test - deriv_map.h (#2): DerivMapGenerator::instance() now throws a descriptive std::logic_error for unsupported (braket, deriv>0) combinations (e.g. xx_xs) instead of std::abort()-ing the host process. xx_xs derivatives remain unimplemented, but now fail catchably rather than crashing consumers. - comp_11_σpσpCoulombσpσp_11.h / comp_11_Coulombσpσp_11.h (#6): replace the guessed `max_nchildren = 100 // TODO` with the true bounds (9 and 3, derived from the densest fold component). This is a children_.reserve() hint only (generic_rr.h), so no behavior change. - test-2body.cc (#1): regression test that default_braket(coulomb_op)==xs_xx and that a coulomb_op engine built WITHOUT an explicit set() produces results identical to the explicit-xs_xx engine (guards the silent-wrong-integrals footgun fixed in the previous commit). Verified: low-AM (MAX_AM=2) compiler->export->library->test all green; RKB Coulomb case 7600->7653 assertions, full suite 2,114,807 / 33 cases, no regressions.
The four σ·p "element-wise fold" recurrence relations (Coulombσp, Coulombσpσp,
σpCoulombσp, σpσpCoulombσpσp) had byte-identical spfunction_call, adapt_dims_,
and generate_code scaffold logic copy-pasted four times; only the operator label
and the per-target RHS expression differ.
Extract the shared logic into three free-function helpers in the new
src/bin/libint/rkb_fold_codegen.h:
- rkb_fold_adapt_dims (the "highdim" dimension swap)
- rkb_fold_spfunction_call (call-site + total_dim emission)
- rkb_fold_generate_code (symbol assignment + function scaffold + the
target[hsi] = <rhs> loop)
Each comp header now delegates spfunction_call/adapt_dims_ in one line and keeps
only its operator-specific RHS computation (copy / nc-keyed / component-keyed)
before calling rkb_fold_generate_code. No CRTP/Instance()/inheritance changes —
the classes keep their exact GenericRecurrenceRelation base; the helpers are
chosen specifically to avoid touching the registry/hierarchy. This relies only
on the public RecurrenceRelation accessors (rr_target/rr_child/num_children/
label), passing the typed target_ where bra()/ket() iteration is needed.
Verified byte-identical: low-AM compiler->export->library->test reproduces the
exact pre-refactor counts (RKB Coulomb 7653 assertions, full suite 2,114,807 /
33 cases), so the generated code is unchanged.
kshitij-05
marked this pull request as ready for review
June 9, 2026 13:25
GenerateCode() (defined in buildtest.h) references the global CodeGenProgress g_progress, which was defined only in build_libint.cc. The standalone compiler test (src/bin/libint/test.cc -> libint-libcompiler-test, the check-libint2compiler target) includes buildtest.h and instantiates GenerateCode but does not link build_libint.cc, so it failed with an undefined `_g_progress` symbol — breaking CI (the macOS Repo job's check-libint2compiler step, which in turn starved the Export job of its compiler artifact). g_progress is an executable-scope global (build_libint defines its own); define it in test.cc as well so the compiler test links. buildtest.h is included only by these two executable mains and never by a library TU (GenerateCode is non-inline), so this cannot be moved into libint2-libcompiler.a without a duplicate-symbol clash — per-main definition is the correct shape. Verified: `cmake --build ... --target check-libint2compiler` now passes (coretest/build + coretest/run, 2/2).
The 3-center (xs_xx) dispatch in compute2 selected the ket stride `ket_lmax` using a hard-coded `LIBINT2_CENTER_DEPENDENT_MAX_AM_3eri`. That macro is 1 for plain 3-center ERI (whose fitting center reaches ERI3_MAX_AM while the paired centers only reach MAX_AM, so the build table is non-cubic and the ket stride must be hard_default_lmax_). But _initialize() sets hard_default_lmax_ per the *actual* task's center-dependence: for the RKB 3-center operators (3coulomb_opop, 3coulomb_op), which are NOT center-dependent, it is INT_MAX. Because the stride selection ignored the real operator and always trusted the 3eri macro, those RKB engines computed buildfnidx = (bra1.l*ket_lmax + ket1.l)*ket_lmax + ket2.l with ket_lmax = hard_default_lmax_ = INT_MAX, i.e. a wildly out-of-bounds index into the symmetric [hard_lmax_]^3 build table — a garbage function pointer and a segfault on the first shell triplet with l(ket1) > 0 (the (s|s,p) μ<ν case; the (s|s,s) case has buildfnidx 0 and happened to survive). Replace the hard-coded, deriv-switched 3eri logic with ket_lmax = min(hard_lmax_, hard_default_lmax_) which is correct for every 3-center task: for center-dependent 3eri it yields hard_default_lmax_ (the paired max, unchanged behaviour); for the non-center-dependent RKB operators it yields hard_lmax_, matching their cubic build-table stride. The μ<->ν canonicalization (swap_tket) was already correct; only the stride was wrong. Verified at MAX_AM=2 with RKB_ERI3 enabled: "RKB Coulomb integrals" 7653 assertions pass (3-center sections now run, previously segfaulted); full unit suite 20,034,963 / 33 cases pass, no regression to plain ERI3.
- cmake.yml: enable LIBINT2_ENABLE_RKB_ERI3=0 (RKB_ERI3_MAX_AM=2) in the shared BUILD_CONFIG so CI actually generates the 3-center RKB kernels and runs the "Coulombσpσp/Coulombσp 3-center" unit sections (previously RKB_ERI3 defaulted to OFF, so that code path was never exercised on CI). - test-2body.cc: compile-guard the three 3-center RKB sections on LIBINT2_TASK_EXISTS_3coulomb_opop / _3coulomb_op. When RKB_ERI3 is disabled, the 4-center coulomb_opop engine still constructs (so the lmax_exceeded catch does not fire), but compute() on the 3-center braket would dispatch to an absent kernel; the guard makes these sections skip cleanly instead of aborting. With RKB_ERI3 enabled (the new CI default) the guards are transparent and the sections run. Verified both ways at MAX_AM=2: RKB_ERI3 on -> full suite 20,034,963 / 33 pass (3-center sections run); RKB_ERI3 off -> 20,034,734 / 33 pass (3-center sections skipped, no abort).
The cxxapi.h availability guard accepted RKB_ERI (4-center) as a sufficient 2-body class but omitted RKB_ERI3 (3-center). A build with 1-body + only RKB_ERI3 enabled (no ERI/ERI3/ERI2/RKB_ERI) would therefore hard-error even though the library provides usable 2-body integrals. Add LIBINT_INCLUDE_RKB_ERI3 to the disjunction and the message. (Copilot review.)
…ss) fundamental names The macOS case-insensitive-filesystem differentiator (commit d4247b0) appends an extra underscore to primitive (non-contracted) function labels in CGF::label(), so the generated struct members, prereq reads, and LIBINT2_DEFINED__ macros for the (ss|ss) fundamentals are now _aB_s____0__s____1___TwoPRep_s____0__s____1___Ab__up_N (4 underscores). include/libint2.h's LIBINT_T_SS_EREP_SS macro was updated to match, but two static consumers still hardcoded the old 3-underscore name: - export/fortran/test-eri.cc - export/fortran/libint_f.F90 As a result every #ifdef LIBINT2_DEFINED__aB_s___0... guard evaluated false, the fundamental-fill code compiled out, and the prereq read uninitialized memory -> the "Fortran ERI" test produced garbage and failed in CI (it passed locally only because a fresh zeroed heap yields rel_error = 0/0 = NaN, a vacuous pass). Update both files to the 4-underscore primitive naming so the convention is applied consistently regardless of platform or Fortran/C++ API.
…tency The engine Operator enum uses ASCII names for the σ·p (RKB) operators (opVop, coulomb_opop, op_coulomb_op, opop_coulomb_opop, coulomb_op); the dipole operator was the lone Unicode outlier. Everything downstream of it was already ASCII opemuop (the descriptor label in oper.h, the BOOST_PP_NBODY_OPERATOR_LIST entry, the task label, and the generated LIBINT2_MAX_AM_opemuop macro), so the enum was the only mismatch. Rename Operator::σpeμσp -> Operator::opemuop (enum value, last_1body_oper, operator_traits specialization, the two engine.impl.h dispatch checks, and the test-1body.cc Engine construction). The compiler-side descriptor class σpeμσp_Descr / σpeμσpOper stays Unicode (matches σpVσp_Descr et al.); the math notation σpeμσp is retained in comments. External consumers now reference the ASCII libint2::Operator::opemuop rather than a Unicode identifier.
CR_1_σpVσp_1 and CR_1_σpeμσp_1 had near-identical constructors: both build the 6 first-derivative basis functions and the same 4-case Pauli-quaternion fold (σ_a σ_b = δ_ab + iε_abc σ_c), differing only in the inner child integral (electrostatic potential V vs Cartesian dipole r_k). Extract the shared structure into sigma_p_1body_fold() in the new rkb_fold_1body.h; each operator now supplies just its inner child via a make_child(Da, Db) lambda. No behavior change (same children, signs, and flop counts). This is the 1-body sibling of the 2-electron rkb_fold_codegen.h sharing already used by the comp_11_*σp*_11.h headers.
3-center RKB ERIs had no pure-solid-harmonic option for the density fitting (auxiliary) center, unlike plain ERI3. Add LIBINT2_RKB_ERI3_PURE_SH following the ERI3_PURE_SH counterpart (default OFF): - CMakeLists.txt: option_with_print + booleanize01 + LIBINT2_->LIBINT_ propagation. - config.h.cmake.in: #cmakedefine LIBINT_RKB_ERI3_PURE_SH. - int_am.cmake: build the RKB pure-SH (_pureamlist) fitting-center components and select them via the same NOT-PURE_SH fallback used by ERI3. - build_libint.cc: apply pure_sh() to the fitting center for the RKB (non-TwoPRep) operators in build_TwoPRep_1b_2k (target + label construction), and propagate the RKB_ERI3_PURE_SH macro to the generated library — structurally identical to the existing ERI3 path. (There is no 2-center RKB ERI, so no RKB_ERI2 counterpart is needed.) All codegen additions are gated on LIBINT_RKB_ERI3_PURE_SH (default OFF), so the default build is unchanged.
Document two known refactors in-place so they are tracked rather than re-discovered: - gauss.cc CGF/CGShell::label(): the primitive-disambiguation underscore leaks into C identifiers, not just filenames; note the preferred funcname->filename mapping. - build_libint.cc / engine.impl.h / strategy.cc: the per-operator std::is_same<OperType,...> / Operator::... branches (component counts, swappability, swap phase/remap, strategy specializations) should move onto operator-descriptor traits and be queried generically. Comments only; no behavior change.
The sigma_p_1body_fold trace component folds 3 first-derivative children and the antisym components fold 2, so max_nchildren is 3, not 100. Set the real bound on both 1-body RKB folds and drop the stale "TODO figure out", matching the child-count hardening already done on the 4-center σ·p folds.
The σpσpCoulombσpσp reference in the "Coulombσpσp and σpσpCoulombσpσp" section builds its 16 components from the didx(a,b,c,d)/D(a,b,c,d) lambdas, which superseded the ~40 hand-written four-derivative der_idx constants (xxxx..xyzz). Those constants were defined but never referenced. Remove them (49 lines). The nine two-derivative constants (d_xx..d_yx) are kept — the (LL|SS) coulomb_opop reference still uses them via eri_drrrr. No behavior change; reference math is unchanged.
The "σpeμσp correctness" test previously only checked Hermiticity (bra↔ket
swap symmetry) and non-triviality — it never compared against known values.
Redesign it to mirror the σpVσp ("W correctness") test: compute the 12
components over a (d|d) shell pair and compare each against a frozen
shellset_ref_standard array via LIBINT2_TEST_ONEBODY.
The reference values are the operator output, independently validated against
the Pauli-quaternion fold of the dipole's geometric derivatives
∂A_a ∂B_b (μ|r_k|ν) — obtained from emultipole1 deriv-2 and, as a cross-check,
finite differences of the deriv-0 dipole — agreeing to ~1e-6. So they are a
verified reference, not a self-referential snapshot. Components 1, 6, 11
(Pauli direction aligned with the dipole direction) vanish by symmetry.
Add .margin(1e-10) to LIBINT2_TEST_ONEBODY so those structurally-zero
components pass; 1e-10 is far below the magnitudes the other one-body value
tests (W/erf/nuclear) check, so their pass/fail is unchanged. Verified: full
unit suite green at both LIBINT2_ENABLE_ONEBODY=0 and =2.
|
Big thank you for implementing these. I paused working on amfX2C in Psi4 until I could write an atomic DHF solver. |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Implement 2-electron 4-center relativistic integrals with restricted kinetic balance condition (RKB).
Implement 2e 3-center relativistic integrals with RKB.
Also adds
σpeμσpoperator for small component dipole moment integrals