Conversation
…odifying the nuclear core intgrals.
… able to switch between SAP and SOAD guess
There was a problem hiding this comment.
Pull request overview
This PR adds a new 1-body operator (Operator::sap) to directly evaluate SAP (Superposition of Atomic Potentials) integrals by reusing the existing nuclear-attraction (“elecpot”) machinery with a modified core integral evaluator, plus helpers to load SAP primitive data from Gaussian94-format files.
Changes:
- Introduces
Operator::sapwith correspondingoperator_traits, core-evaluator (sap_gm_eval), and Engine plumbing. - Adds SAP basis parsing/loading helpers (
read_sap_basis_library,make_sap_prim_data) and refactors basis-path/Fortran-float helpers into free functions. - Extends tests and the Hartree–Fock sample to compute and use SAP integrals (including an optional SAP initial guess).
Reviewed changes
Copilot reviewed 9 out of 11 changed files in this pull request and generated 6 comments.
Show a summary per file
| File | Description |
|---|---|
| include/libint2/shell.h | Adds SAP primitive/container type aliases (SAPPrimitive, SAPElementData, SAPElementsData). |
| include/libint2/lcao/1body.h | Ensures 1-body integral driver sets engine params for Operator::sap. |
| include/libint2/engine.impl.h | Registers sap in operator list, routes sap through nuclear-like code paths, and wires sap core-eval into compute_primdata. |
| include/libint2/engine.h | Adds Operator::sap and operator_traits<Operator::sap>. |
| include/libint2/boys_fwd.h | Forward-declares sap_gm_eval. |
| include/libint2/boys.h | Implements sap_gm_eval and its scratch storage. |
| include/libint2/basis.h.in | Adds SAP basis reading/loading helpers and refactors basis_data_path / fortran_dfloats_to_efloats. |
| export/tests/unit/test-1body.cc | Adds a unit test validating a reference SAP matrix. |
| export/tests/hartree-fock/hartree-fock++.cc | Adds SAP-based initial guess option and SAP matrix construction helper. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
1. Locale — Added #ifdef _MSC_VER for "en-US" / "POSIX" in read_sap_basis_library 2. & 4. Assert → throw — Replaced assert with throw std::logic_error for non-s-type shells 3. lmax — Changed to orbital_basis.max_l() instead of LIBINT2_MAX_AM_elecpot 4. Atomic number cast — Added std::round + assert for integer charge validation 5. Division by zero — Added early return with zeroed Gm when q == 0
There was a problem hiding this comment.
Pull request overview
This PR introduces a generalized Gaussian-charge nuclear potential operator (Operator::q_gau) to enable direct evaluation of SAP integrals (and related finite/attenuated nuclear models) by reusing the existing nuclear attraction integral machinery with modified core integral evaluation.
Changes:
- Adds
Operator::q_gauplus supporting parameter types for per-center Gaussian potential expansions (including point charge, finite nuclear, SAP, erf/erfc/erfx). - Implements generalized core-integral evaluation (
q_gau_gm_eval) and wiresq_gauthrough the Engine/operator plumbing. - Adds helpers to build per-center potential data from atoms / SAP basis files, plus unit tests and a Hartree–Fock example hook for SAP guesses.
Reviewed changes
Copilot reviewed 10 out of 12 changed files in this pull request and generated 4 comments.
Show a summary per file
| File | Description |
|---|---|
| include/libint2/shell.h | Adds Gaussian potential primitive/center data types used by q_gau. |
| include/libint2/lcao/1body.h | Ensures compute_1body_ints sets params for q_gau like other nuclear-type operators. |
| include/libint2/engine.impl.h | Integrates q_gau into operator lists, parameter handling, and core integral evaluation dispatch. |
| include/libint2/engine.h | Adds Operator::q_gau and its operator_traits definition. |
| include/libint2/chemistry/elements.h | Adds Gaussian nuclear exponent utilities and isotope mass-number table. |
| include/libint2/boys_fwd.h | Forward-declares q_gau_gm_eval. |
| include/libint2/boys.h | Implements q_gau_gm_eval and scratch specialization. |
| include/libint2/basis.h.in | Adds SAP-basis reader + make_q_gau_data* helpers for building per-center potential data. |
| export/tests/unit/test-1body.cc | Adds unit tests comparing q_gau to existing nuclear/erf/erfc/erfx operators and SAP reference data. |
| export/tests/hartree-fock/hartree-fock++.cc | Adds optional SAP-based initial guess path via Operator::q_gau. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
…of `q_gau`. Supports `opVop`,`op(Vsap)op`, `op(erfV)op`, `op(erfcV)op`, `op(erfxV)op`, etc. with both finite and point nuclear models
evaleev
left a comment
evaleev
left a comment
There was a problem hiding this comment.
Review of PR #408: Direct evaluation of SAP integrals
Summary
This PR introduces Operator::q_gau (and its relativistic small-component variant Operator::op_q_gau_op) — a generalized Gaussian-charge nuclear potential operator that unifies point nuclear, finite (Gaussian) nuclear, SAP, erf/erfc/erfx models under a single operator. The design is clean: each center carries a shared_ptr<GaussianPotentialData> (a vector of {exponent, coefficient} primitives), and the core integral evaluator (q_gau_gm_eval) simply sums contributions. This is a significant generalization that subsumes all existing nuclear-type operators.
Architecture & Design
Strengths:
- The
GaussianPotentialPrimitive/GaussianPotentialData/GaussianPotentialCentersDatatype hierarchy is well-designed. Usingshared_ptrfor same-Z sharing andnullptrfor ghost atoms is practical. - The
make_q_gau_data*convenience functions cover all common use cases (point, finite, SAP, erf, erfc, erfx) with a consistent interface. - Excellent test coverage: unit tests verify q_gau matches every existing nuclear operator variant (nuclear, erf, erfc, erfx), plus Gaussian nuclear model and SAP with hardcoded reference data. The op_q_gau_op tests verify all 4 quaternionic components.
- Refactoring
basis_data_path()andfortran_dfloats_to_efloats()into free functions is a nice cleanup.
The quaternionic reordering in oper.h (0=scalar, 1=X, 2=Y, 3=Z instead of 0=scalar, 1=Z, 2=X, 3=Y) — is this intentional? This changes the component ordering for the existing opVop operator too, since σpVσp_Descr is shared. This could break downstream code that relies on the old ordering. The engine.h comment Component order: 0 = scalar, 1 = x, 2 = y, 3 = z is added for both opVop and op_q_gau_op, but if this is a change to opVop's ordering, it's an API break that needs calling out explicitly.
Issues
1. (Important) opVop component order change is an API break
The rename from pauli_index to quaternionic_index and reordering X/Y/Z in oper.h affects the generated recurrence code for opVop. Any existing code using opVop with index assumptions (1=Z, 2=X, 3=Y) will silently get wrong results. This should be documented as a breaking change, or the old operator's ordering should be preserved and only op_q_gau_op should use the new ordering.
2. (Medium) #include <map> in shell.h is unused
The shell.h header adds #include <map> but no std::map is used there. <map> is used in basis.h.in (in make_q_gau_data), which is fine. Remove from shell.h.
3. (Medium) SAP basis file parse error handling
In read_sap_basis_library(), after iss >> amlabel >> nprim >> rest and iss2 >> e >> c, stream extraction failures are not checked. On a malformed file, nprim or e/c could be indeterminate. Add stream-state checks after extraction.
4. (Minor) gaussian_nuclear_exponent_from_A lacks input validation
Unlike gaussian_nuclear_exponent(int Z) which validates Z ∈ [1,118], gaussian_nuclear_exponent_from_A(double A) accepts any value including A ≤ 0. Consider adding a guard.
5. (Minor) Ho isotope mass number
/* 67 Ho */ 162 — Holmium's most abundant (and only stable) isotope is Ho-165, not 162. 162 is Dy's value (the line above). Looks like a copy-paste error.
6. (Minor) Boilerplate in tests could be reduced
The gau_max_nprim calculation pattern is repeated ~8 times across tests. A small helper like max_nprim(const GaussianPotentialCentersData&) would clean this up.
7. (Cosmetic) T.resize(0,0) moved in hartree-fock++.cc
The kinetic energy matrix T deallocation was moved from after H = T + V to after the guess block. This means T stays alive through the guess computation. For the SAP guess this is intentional (since F = T + V_sap), but for SOAD it's now unnecessarily held. Very minor memory concern.
Questions for the author
- Is the
opVopcomponent reordering intentional as a breaking change, or should theσpVσp_Descrordering be left alone for backward compatibility? - The
q_gau_gm_evalcallsfm_eval_->eval()once per primitive per shell pair. For SAP bases with many primitives per element (the files show ~30+ primitives for heavy atoms), has the performance impact been benchmarked vs. the alternative of pre-contracting? - Your comment on the PR mentions SAP centers need explicit positions for QM/MM. Is the current design (positions come from the point_charges vector) sufficient, or does this need a follow-up?
|
also: need to update the wiki |
Introduces generalized Gaussian potential operators (
q_gauandop_q_gau_op)Operator::q_gauandOperator::op_q_gau_opuse a unified core evaluator parameterized by a list of Gaussian primitives{(alpha, coefficient)}per center. All standard nuclear models are special cases:{alpha, coefficient}{INFINITY, 1.0}{xi(Z), 1.0}{omega*omega, 1.0}{INFINITY, 1.0}, {omega*omega, -1.0}{INFINITY, sigma}, {omega*omega, -(sigma-lambda)}{INFINITY, 1.0}, {a1, c1/q}, {a2, c2/q}, ...{xi(Z), 1.0}, {a1, c1/q}, {a2, c2/q}, ...