Repository navigation
Add truncated Coulomb support to libint and expose it to Engine - #434
Conversation
|
Added another implementation. It seems to perform well, both in speed and accuracy. |
|
Thanks for V2 — we tested it independently and it looks like the right way forward. What we measured (standalone harness, V2 sources from fc257fa, Apple M-series, single core):
So we'd like to drop V1 entirely — the ~107 MB of generated tables, Request: to keep the tables out of the repository history, the branch needs to be rewritten, not just have V1 deleted in a later commit. The target is roughly:
You're welcome to do this yourself; alternatively, may we rewrite and force-push Follow-ups we'd like before merging (happy to help with these):
This comment (analysis, tests and text) was agent-generated and human-reviewed. |
|
Yes, please feel free to rewrite and force-push Thank you very much for the independent tests — I am really glad to see that V2 reaches the double-precision VRR floor. Also, V1 can be very safely deleted. I previously intend to keep that for benchmarking, now there's no need for V1. |
|
Thanks — we'll go ahead with the rewrite. One more item for the V2 review: performance on hardware without FMA. On an Ivy Bridge Xeon (E5-2690 v2, AVX but no FMA), V2 is ~25× slower than on Apple M-series, and 7× slower than a table-based evaluator we prototyped for comparison:
Without a hardware FMA instruction, glibc computes The replacement above is not a fix, since the rounding analysis relies on fused operations in places. Options: an FMA-free formulation (e.g. Dekker/Veltkamp two-product, or an explicitly widened error budget as This comment (analysis, tests and text) was agent-generated and human-reviewed. |
|
Correction to the FMA comment above: a default That case is now addressed in 62c8e16 (260a866 before the rebase onto master): Also pushed: the This comment (analysis, tests and text) was agent-generated and human-reviewed. |
…= T_crit The Chebyshev table has T_crit/delta = 819 intervals covering [0, T_crit). At exactly T = T_crit the interval index int(T/delta) is 819, one past the end of the table, so the interpolation read out of bounds. Add a test of values at and around the endpoint against the long-double reference.
Operator::truncated_coulomb_v2 takes the cutoff radius Rc (in the same length units as Shell::O) as its parameter; Rc = 0 gives zero integrals and Rc = infinity ordinary Coulomb integrals. The radial auxiliary integrals G_m(R, T) are evaluated without interpolation tables, in binary64 only, with explicit error budgets (os_core_ints::truncated_coulomb_v2_gm_eval). Tests compare G_m against independent 100-digit references generated by tools/generate_truncated_coulomb_reference.py (mpmath), and Engine integrals against the Coulomb and zero limits, finite differences, and braket/center permutations.
With the table-based evaluator removed there is a single implementation: Operator::truncated_coulomb, os_core_ints::truncated_coulomb_gm_eval, libint2/truncated_coulomb.h and src/lib/libint/truncated_coulomb.cc.
Without a compiler-emitted FMA instruction (FP_FAST_FMA undefined, e.g. x86-64 builds without -mfma), std::fma is a library call, emulated in software on CPUs without FMA: ~150 ns per call on Ivy Bridge, ~90% of the evaluation time. Replace the 26 calls by - difference_of_square(T, R): T - R^2 with Dekker's exact product, which keeps the accuracy of the displacement near R = sqrt(T); - outward_mul_add(): rigorous bounds with one outward step per rounding; - mul_add(): plain a*b+c elsewhere, as taylor_add already did. With FP_FAST_FMA the results are unchanged. LIBINT2_TRUNCATED_COULOMB_FAST_FMA=0 forces the fallback for testing. Ivy Bridge (E5-2690 v2), -march=x86-64, mmax=16: 11890 -> 990 ns per call. Integral-level error ((L0|00) by VRR over random geometries across the cutoff transition, L <= 16) is unchanged with either path.
The G_m tests check absolute errors. Add a test of what matters for integrals: (L0|00), L <= 16, by the Obara-Saika VRR from evaluated vs. reference G_m, relative to the larger of the truncated and the Coulomb integral, within 8x the double-precision VRR floor 2^(L/2) eps. The 96 geometries (tools/generate_truncated_coulomb_reference.py --vrr-output) are physical primitive pairs across the cutoff transition: exponents in [1e-2, 1e6], Rc in [1, 100], sqrt(T) - R in [-16, 16]; the first three are the worst cases of a 20000-geometry survey. The reference G_m are checked against 40 extra digits relative to max(|G_m|, F_m), and agree bit for bit with an independent Boost.Multiprecision evaluation. The test fails for G_m with a relative error of 1e-12, for Coulomb F_m, and for G_m replaced by F_m wherever |G_m - F_m| < eps F_0.
…cratch The evaluator keeps mutable work storage (series coefficients, recurrence and Miller buffers), so it does not fit GenericGmEval, whose wrapped evaluators must be stateless; that needed an eval_unchecked hook in GenericGmEval and two CoreEvalScratch specializations. Instead os_core_ints::truncated_coulomb_gm_eval is now the core_eval_type itself, with instance(), max_m(), precision() and eval(Gm, rho, T, mmax, cutoff, scratch); Engine passes its CoreEvalScratch, which is the evaluator's Scratch. The kernel has one entry point, detail::truncated_coulomb::evaluate(); its resource limits are internal. Requests are validated on every call (a few comparisons); results are no longer scanned for non-finite values. A requested precision below double epsilon (e.g. from a long double Engine) is raised to it instead of throwing. Public headers stay C++11.
The kernel's error bounds rely on infinities as sentinels, compensated sums, Dekker's exact product and outward rounding. Value-unsafe modes break them: the Intel icx export build, whose default is -fp-model=fast, failed the integral test at large R, and a kernel built with -ffast-math fails both truncated-Coulomb test cases with or without FMA. Compile truncated_coulomb.cc with value-safe options appended after the user's flags (-fp-model=precise for icx, /fp:precise for MSVC-compatible compilers, -fno-fast-math for GCC and Clang), and fail the build if the file still sees __FAST_MATH__ or __FINITE_MATH_ONLY__.
…r R = sqrt(T) - Accept the F_m and zero limits only when their error bound is below precision * F_n(T) (rigorous lower bound on F_n via lgamma). The absolute test replaced high-order G_m at large T by values without a correct digit (e.g. G_16 = -2.8e-62 vs F_16 = 8.2e-71 at T = 1e5); no measurable cost. - For R >= sqrt(T), fall back to the erfc-difference form of G0 when the F0-based form cancels (ratio up to 3 at R = sqrt(T)), instead of running radius and Taylor series that cannot converge there: 3.7x faster in that band with unchanged accuracy. - Tests: G_m relative to max(|G_m|, F_m) just beyond the former switch distance (fails before this change); Engine (ss|ss) against the closed form, checking cutoff units and the sqrt(rho) scaling independently of the kernel.
GCC contracts a*b+c into a fused multiply-add across statements by default (-ffp-contract=fast, unaffected by -fno-fast-math). Fusing a product into a following TwoSum makes the compensated sums inexact; g++ -O2 -mfma does so for a product followed by TwoSum. The kernel calls std::fma itself where it wants one, so pass -ffp-contract=off to GCC, Clang and icx.
… T < n+1/2
For 1 <~ T < n + 1/2 the downward branch needs Q_m that only Miller's
boundary solve supplied, at 3-6 us on Apple M-series and 9-13 us on Ivy
Bridge (no FMA); it dominated the mean cost at mmax >= 8.
Evaluate all orders directly from the Taylor expansion about T = 0,
G_m = F_m(T) - exp(-x) sum_j (-T)^j / j! P_{m+j}(x) / (2m+2j+1),
with x = R^2 and P_k = 1F1(-k; 1/2; x), in double-word arithmetic. The
absolute terms grow like exp(T), so binary64 would lose log10 exp(T) digits;
double-word keeps the error near u^2 exp(T). (1/2)_k P_k obeys a
division-free three-term recurrence; double-word tables supply
1/((1/2)_k (2k+1)) and 1/k!. Used for 0 < T < min(n + 1/2, 12), outside the
small-R radius-series regime.
Each call bounds the truncation and the double-word rounding of the weights,
the coefficients and the accumulation, and estimates the propagated
recurrence error to first order from its Green's function: a second solution
V and the closed-form Casoratian i! (1/2)_i. A term-by-term bound suits the
oscillating regime; where both solutions grow (x >~ 4k) it overestimates by
up to 1e9, and a nested form with positive prefix sums is evaluated instead.
If the total exceeds precision/4, or the order exceeds the tables (n + degree
>= 160) or the magnitudes could overflow, the existing path is used.
Physical-distribution kernel mean (Apple M-series), before -> after:
mmax 8: 776 -> 223 ns (p99 7.1 -> 1.1 us); 16: 1384 -> 529 ns; 32: 2540 ->
1586 ns. The estimate costs ~15% of the new path. Accuracy against 300-digit
references is unchanged at mmax 8, 16 and 32; the fuzz maximum absolute error
stays 1.8e-16. A new test covers the top of the T range at order 32, a
large-x point that needs the nested estimate, and agreement with the other
paths on the same inputs.
The export sets C++11 as its minimum standard. Options{precision} is not
aggregate initialization in C++11, whose aggregates cannot have default
member initializers; give Options a constructor. The kernel now compiles
with -std=c++11 -pedantic.
Adds
Operator::truncated_coulomb, the spherically truncated Coulomb kernel θ(Rc − r12)/r12 (the truncated-Coulomb formulation used in CP2K), with the cutoff radius Rc as its parameter, in the length units ofShell::O. Rc = 0 gives zero integrals and Rc = ∞ the Coulomb integrals.The auxiliary integrals Gm(R, T), with R = Rc√ρ, are evaluated by
os_core_ints::truncated_coulomb_gm_eval(libint2/truncated_coulomb.h,src/lib/libint/truncated_coulomb.cc). There are no interpolation tables and all arithmetic is binary64, with an explicit error budget on every path. The target is absolute precision, like the Boys and Ten-no evaluators. The Coulomb and zero limits are used only where they are exact relative to Fm(T).Commits (on master
0ede003f)ab3908dcFmEval_Chebyshev7: use the asymptotic branch at the table endpoint T = Tcrit. An out-of-bounds read, found along the way.a314bd1eAdd the truncated Coulomb operator. This is the V2 evaluator, with Engine support and tests against 100-digit mpmath references.0a78605eDrop the_v2suffix.62c8e165Usestd::fmaonly where it is a hardware instruction. On CPUs without FMA it was emulated in software.668acbd3Integral-level test: (L0|00), L ≤ 16, by Obara–Saika VRR over 96 physical geometries across the cutoff transition, against stored references.58834295Make the evaluator an Engine core evaluator with per-Engine scratch.15a87a87Compile the kernel with IEEE floating-point semantics. The Intel icx default is-fp-model=fast.27072883Relative Coulomb/zero limits, and a well-conditioned G0 near R = √T.7a9433e6Compile the kernel without FMA contraction. GCC fuses products into TwoSum by default.98656217A double-word Taylor expansion about T = 0 replaces Miller's algorithm for 0 < T < min(n + ½, 12). Every call bounds its own error and falls back to the previous path if the bound misses the target.725fe382Compile as C++11.The PR's earlier table-based evaluator (V1, with ~107 MB of generated tables) is not in the history.
Accuracy
Performance
These timings are for the kernel alone (ns per call, one core) on physical primitive pairs. Exponents are log-uniform in [10−2, 106], Rc is log-uniform in [1, 30] bohr and |PQ| is uniform in [0, 3Rc]. The three kernels compared are:
27072883, before the origin expansion.725fe382.Apple M-series:
Intel Xeon E5-2690 v2 (Ivy Bridge, no FMA), g++ 9.4,
-O3 -march=native:At mmax ≥ 16 the p99 and much of the mean come from calls with 12 < T < mmax + ½. The origin expansion covers only T ≤ 12, so these calls still use Miller's algorithm. On M-series they are 2.5% of calls and 35% of the time at mmax = 16, and 7.9% of calls and 57% of the time at mmax = 32.
History
The original description, by RTChou67:
The formulation documents are truncated-coulomb.pdf (V1, since removed) and truncated-coulomb-v2.pdf (V2, the basis of this PR).