Add truncated Coulomb support to libint and expose it to Engine

export/CMakeLists.txt.export

@@ -237,6 +237,28 @@ configure_file(
 libint2_install_path(_libint2_installed_datadir "${LIBINT2_INSTALL_DATADIR}" "\$\{_IMPORT_PREFIX\}")
 
 include(srclist.cmake)
+# Keep the truncated-Coulomb numerical kernel in its own translation unit,
+# compiled with IEEE floating-point semantics: its error bounds rely on
+# infinities, compensated sums and exact products, which value-unsafe modes
+# (-ffast-math, and icx's default -fp-model=fast) break. Contraction of a*b+c
+# into a fused multiply-add is disabled as well: GCC's default contracts
+# across statements, which makes compensated sums inexact; the kernel calls
+# std::fma itself where it wants one. The options come after
+# CMAKE_CXX_FLAGS on the command line, so they override them.
+if (MSVC)  # includes clang-cl and icx-cl
+    set(_libint2_ieee_fp_options "/fp:precise")
+elseif (CMAKE_CXX_COMPILER_ID STREQUAL "IntelLLVM")
+    set(_libint2_ieee_fp_options "-fp-model=precise;-ffp-contract=off")
+elseif (CMAKE_CXX_COMPILER_ID MATCHES "GNU|Clang")
+    set(_libint2_ieee_fp_options "-fno-fast-math;-ffp-contract=off")
+endif()
+foreach(_source IN LISTS LIBINT2_LIBRARY_CXX_SRC)
+    if (_source MATCHES "truncated_coulomb")
+        set_source_files_properties(${_source} PROPERTIES
+            SKIP_UNITY_BUILD_INCLUSION ON
+            COMPILE_OPTIONS "${_libint2_ieee_fp_options}")
+    endif()
+endforeach()
 list(LENGTH LIBINT2_LIBRARY_CXX_SRC _source_count)
 message(STATUS "Loading ${_source_count} library source files from LIBINT2_LIBRARY_CXX_SRC")
 message(DEBUG "LIBINT2_LIBRARY_CXX_SRC=${LIBINT2_LIBRARY_CXX_SRC}")

export/tests/CMakeLists.txt

@@ -85,6 +85,7 @@ if (LIBINT2_REQUIRE_CXX_API)
             unit/test-core.cc
             unit/test-core-ints.cc
             unit/test-molden.cc
+            unit/test-truncated-coulomb.cc
             unit/test-permute.cc
             unit/test-precision.cc
             unit/test-shell-order.cc

export/tests/unit/test-core-ints.cc

@@ -92,6 +92,25 @@ TEST_CASE("Boys reference values", "[core-ints][precision]") {
 
 #endif  // LIBINT_HAS_MPFR
 
+TEST_CASE("Boys interpolation endpoint", "[core-ints][precision]") {
+  const auto evaluator = libint2::FmEval_Chebyshev7<double>::instance(40);
+  // The table covers [0,117); its right endpoint uses the asymptotic branch.
+  for (double T :
+       {std::nextafter(117., 0.), 117., std::nextafter(117., 118.)}) {
+    double values[41];
+    long double reference[41];
+    libint2::FmEval_Reference2<long double>::eval(reference, T, 40);
+    for (int last : {0, 40}) {
+      evaluator->eval(values, T, last);
+      for (int m = 0; m <= last; ++m) {
+        INFO("T=" << T << " m=" << m);
+        CHECK(std::abs(values[m] - reference[m]) <=
+              std::numeric_limits<double>::epsilon());
+      }
+    }
+  }
+}
+
 TEST_CASE("Slater/Yukawa core integral values", "[core-ints]") {
   using scalar_type = libint2::scalar_type;