Direct evaluation of SAP integrals

[kshitij-05]

Introduces generalized Gaussian potential operators (q_gau and op_q_gau_op)

Operator::q_gau and Operator::op_q_gau_op use a unified core evaluator parameterized by a list of Gaussian primitives {(alpha, coefficient)} per center. All standard nuclear models are special cases:

Model	Primitives {alpha, coefficient}	# primitives
Point nuclear	{INFINITY, 1.0}	1
Gaussian nuclear	{xi(Z), 1.0}	1
erf-attenuated	{omega*omega, 1.0}	1
erfc-attenuated	{INFINITY, 1.0}, {omega*omega, -1.0}	2
erfx	{INFINITY, sigma}, {omega*omega, -(sigma-lambda)}	2
SAP	{INFINITY, 1.0}, {a1, c1/q}, {a2, c2/q}, ...	1 + n
SAP + finite nuclear	{xi(Z), 1.0}, {a1, c1/q}, {a2, c2/q}, ...	1 + n

[Copilot]

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::sap with corresponding operator_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.

[evaleev]

n/a

[Copilot]

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_gau plus 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 wires q_gau through 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.

[evaleev]

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 / GaussianPotentialCentersData type hierarchy is well-designed. Using shared_ptr for same-Z sharing and nullptr for 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() and fortran_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

1. Is the opVop component reordering intentional as a breaking change, or should the σpVσp_Descr ordering be left alone for backward compatibility?
2. The q_gau_gm_eval calls fm_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?
3. 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?

[evaleev]

also: need to update the wiki