From 2a87d69a39e8a3ee939ba318761523b27379cb0d Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Mon, 19 Oct 2026 14:02:09 -0500 Subject: [PATCH] claude's first version --- README.md | 199 +++++++--- examples/eri_error_explorer/README.md | 83 +++++ .../eri_error_explorer/eri_error_explorer.py | 340 ++++++++++++++++++ 3 files changed, 573 insertions(+), 49 deletions(-) create mode 100644 examples/eri_error_explorer/README.md create mode 100644 examples/eri_error_explorer/eri_error_explorer.py diff --git a/README.md b/README.md index 74c5fa25..9645a634 100644 --- a/README.md +++ b/README.md @@ -19,58 +19,159 @@ Integrals [![Nightly Workflow](https://github.com/NWChemEx/Integrals/actions/workflows/nightly.yaml/badge.svg)](https://github.com/NWChemEx/Integrals/actions/workflows/nightly.yaml) -Repository for the generic integral API and implementations for specific integral libraries. +Repository for the generic integral API and implementations for specific +integral libraries. Integrals is a PluginPlay plugin; its modules are written +in C++ and are also usable from Python. -Building Integrals ------------------- +Installing a Pre-Built Wheel +---------------------------- -Integrals is built using -[CPP](https://github.com/CMakePackagingProject/CMakePackagingProject.git). -Assuming you have already installed CPP, that you are on a sane Unix-like -computer, and you are willing to let Integrals build all dependencies, then -the following will suffice to build Integrals: +If you just want to use Integrals from Python, install the published wheel. +The wheels are built for Python 3.12: +```.sh +python3.12 -m venv venv +source venv/bin/activate +pip install nwchemex-integrals ``` + +The rest of this file describes building Integrals from source. + +Prerequisites +------------- + +- A C++20 compiler. CI uses GCC 14 or Clang 18 on Linux and AppleClang on + macOS. On macOS do **not** use GCC: the published wheels of Integrals' + dependencies are built with AppleClang (libc++), and GCC (libstdc++) cannot + link against them. The build stops with an error if you try. +- CMake 3.14 or newer. Ninja is optional but recommended. +- Python 3 with its development headers, and pip. Use Python 3.12 if you want + the build to reuse the published wheels of the dependencies (see + [How dependencies are found](#how-dependencies-are-found)). +- MPI, Boost (headers only) and Eigen3. The build can not obtain these for you. + +On Ubuntu: + +```.sh +sudo apt-get install g++-14 cmake ninja-build libopenmpi-dev openmpi-bin \ + libboost-dev libeigen3-dev +``` + +On macOS, with Homebrew: + +```.sh +brew install cmake ninja open-mpi boost eigen +``` + +How dependencies are found +-------------------------- + +Integrals depends on [SimDE](https://github.com/NWChemEx/SimDE) (and through it +the rest of the NWChemEx stack: PluginPlay, Chemist, TensorWrapper, +ParallelZone, Utilities, ...) and on [Libint](https://github.com/evaleev/libint). +The build system comes from [NWXCMake](https://github.com/NWChemEx/NWXCMake). +You do not need to install any of these yourself: + +- **NWXCMake** is taken from the `nwchemex-nwxcmake` pip package if it is + installed in the Python CMake finds; otherwise it is downloaded from GitHub. +- **NWChemEx dependencies** (SimDE and everything below it) are taken from + their `nwchemex-*` wheels if they are installed in the active Python + environment. Otherwise they are downloaded from the `master` branch on GitHub + and built as part of this build, which takes considerably longer. +- **Libint** is taken from an existing installation if CMake can find one + (`find_package(Libint2)`, e.g., via `CMAKE_PREFIX_PATH`). Otherwise the + v2.11.0 release is downloaded and built, which takes a while. + +The quickest builds come from a Python 3.12 virtual environment that has the +dependency wheels installed, which is what the instructions below set up. + +Building with pip +----------------- + +This builds Integrals and installs it, including the Python module, into the +active virtual environment: + +```.sh +git clone https://github.com/NWChemEx/Integrals.git +cd Integrals +python3.12 -m venv venv +source venv/bin/activate +pip install scikit-build-core setuptools-scm nwchemex-nwxcmake nwchemex-simde +pip install . --no-build-isolation +``` + +The third `pip install` puts Integrals' build requirements (listed under +`[build-system]` in `pyproject.toml`) into the environment. Together with +`--no-build-isolation`, this makes Integrals link against the dependency wheels +that stay installed, rather than copies in a temporary build environment that +pip deletes afterwards. + +For development, do an editable install instead: + +```.sh +pip install -e ".[dev]" --no-build-isolation +``` + +This turns on `DEVELOPER_SETUP` (see the table below), so it also builds the +tests, uses a Debug build and installs the pre-commit hooks. The build tree is +`build/`. + +Building with CMake +------------------- + +Integrals can also be built with CMake directly. Activating the virtual +environment from the previous section first is optional, but if the dependency +wheels are installed there, CMake reuses them: + +```.sh git clone https://github.com/NWChemEx/Integrals.git -cd SDE -cmake -H. -Bbuild -DCPP_GITHUB_TOKEN= \ - -DMPI_ROOT= \ - -DCMAKE_PREFIX_PATH= \ - -DCMAKE_INSTALL_PREFIX= -cd build -cmake --build . -#May need to run as an admin depending on where you are installing -cmake --build . --target install +cd Integrals +cmake -B build -G Ninja \ + -DPython_EXECUTABLE="$(command -v python3)" \ + -DBUILD_TESTING=ON \ + -DCMAKE_INSTALL_PREFIX= +cmake --build build --parallel +cmake --install build ``` -The build process is not capable of building MPI so you will have to provide a -path to a known installation via the `MPI_ROOT` variable. The GitHub token is -necessary because, at the moment, Chemist, TAMM, SDE, and Utilities are -private repositories (instructions for generating a token are -[here](https://help.github.com/articles/creating-a-personal-access-token-for-the-command-line/)). - -For finer-grained control over the build we direct the reader to the more -thorough CPP build instructions located -[here](https://cmakepackagingproject.readthedocs.io/en/latest/end_user/quick_start.html) - -and note that Integrals depends on several other projects: - -- [utilities](https://github.com/NWChemEx/Utilities) -- [SDE](https://github.com/NWChemEx/SDE) - - [bphash](https://github.com/bennybp/BPHash) - - [pybind11](https://github.com/pybind/pybind11) - - Requires a development version of Python - - [cereal](https://github.com/USCiLab/cereal) -- [Chemist](https://github.com/NWChemEx/Chemist) -- [TAMM](https://github.com/NWChemEx/TAMM) - - [HPTT](https://github.com/ajaypanyala/hptt) - - [MSGSL](https://github.com/Microsoft/GSL) - - [LAPACKE](http://www.netlib.org/lapack/) - - [Eigen3](https://github.com/eigenteam/eigen-git-mirror) - - [GlobalArrays](https://github.com/GlobalArrays/ga) - - MPI -- [LibInt](https://github.com/evaleev/libint) -- (For testing only)[Catch2](https://github.com/catchorg/Catch2) - -We recommend building Libint separately for the moment because a proper -autotools build (the only autotools variant supported by CPP) builds the -generator and many more integrals than are necessary. + +Setting `Python_EXECUTABLE` makes sure CMake uses the Python of your active +environment (and therefore its wheels) rather than some other Python it finds +first. Commonly used options are: + +| Option | Default | Description | +|--------------------------------------|---------|-----------------------------------------------------------------------------------------------------| +| `BUILD_TESTING` | `OFF` | Build the C++ unit tests and register the Python tests with CTest. | +| `BUILD_PYBIND11_BINDINGS` | `ON` | Build the `integrals` Python module. | +| `ENABLE_SIGMA` | `OFF` | Use [Sigma](https://github.com/QCUncertainty/sigma) for uncertainty quantification (UQ modules). | +| `INTEGRALS_ENABLE_OPENMP` | `OFF` | Thread the integral evaluation with OpenMP. On macOS, install `libomp` and add its prefix to `CMAKE_PREFIX_PATH`. | +| `DEVELOPER_SETUP` | `OFF` | Debug build, turns on `BUILD_TESTING`, and installs the pre-commit hooks. | +| `ENABLE_XHOST` | `ON` | Tune dependencies for this machine's CPU. Turn off if the build must run on other machines. | +| `NWX_ECOSYSTEM_FROM_SOURCE` | `OFF` | Ignore installed `nwchemex-*` wheels and build every NWChemEx dependency from GitHub. | +| `FETCHCONTENT_SOURCE_DIR_` | unset | Use a local checkout of dependency `` (e.g. `FETCHCONTENT_SOURCE_DIR_SIMDE=/path/to/SimDE`). | + +If you change which Python or which dependency wheels the build should use, +delete `build/` and configure again; earlier choices are cached. + +Running the tests +----------------- + +After a CMake build with `-DBUILD_TESTING=ON` (or an editable pip install), +run the C++ and Python tests with: + +```.sh +ctest --test-dir build --output-on-failure +``` + +After a pip install, the Python tests can also be run directly: + +```.sh +pip install pytest +pytest tests/python +``` + +Examples +-------- + +The [examples](examples) directory contains small programs showing how to use +Integrals. `examples/primitive_error_models` is built along with the tests; +the Python examples have their own `README.md`. diff --git a/examples/eri_error_explorer/README.md b/examples/eri_error_explorer/README.md new file mode 100644 index 00000000..3c355057 --- /dev/null +++ b/examples/eri_error_explorer/README.md @@ -0,0 +1,83 @@ + + +ERI Error Explorer +================== + +`eri_error_explorer.py` explores how the basis set parameters, the separation +between centers, and the screening tolerance affect the error in the +four-center electron repulsion integrals (ERIs). See the docstring at the top +of the script for details. + +Setting up +---------- + +1. Build Integrals, with its Python bindings, as described in the + [top-level README](../../README.md). Building with pip in a Python 3.12 + virtual environment is the simplest option. +2. In the same environment, install ChemCache (the script takes its STO-3G + basis set from ChemCache) and NumPy: + + ```.sh + pip install nwchemex-chemcache numpy + ``` + + ChemCache and Integrals must share one copy of SimDE and the libraries + below it. That is the case if Integrals was built against the dependency + wheels installed in this environment, which the top-level README's + instructions do. If you built Integrals with `NWX_ECOSYSTEM_FROM_SOURCE=ON`, + or in an environment without the wheels, build ChemCache from source the + same way instead of installing its wheel. +3. If you built Integrals with CMake and did not install it into the + environment, add the build directory to `PYTHONPATH` so Python can find the + `integrals` module: + + ```.sh + export PYTHONPATH=/path/to/Integrals/build:$PYTHONPATH + ``` + + This step is not needed after a pip install. + +Running +------- + +From the root of the Integrals repository: + +```.sh +python examples/eri_error_explorer/eri_error_explorer.py +``` + +With the default settings this takes a few seconds. The script prints the +sampled shells, one line per shell quartet and tolerance, and a summary per +separation and tolerance. + +Changing what is explored +------------------------- + +The parameters are hard-coded near the top of the script. Edit them and run +the script again: + +| Parameter | Default | Meaning | +|-------------------------|--------------------|-----------------------------------------------------------------------------------| +| `SHELL_SAMPLES` | `[(1, 0), (8, 2)]` | `(Z, i)` pairs: shell `i` of the STO-3G basis set for atomic number `Z` (H 1s, O 2p). | +| `R_MIN`, `R_MAX` | `1.0`, `5.0` | Range of center separations, in bohr. | +| `N_SEPARATIONS` | `3` | Number of evenly spaced separations in that range. | +| `TOLERANCES` | `[1e-6, 1e-8]` | Screening tolerances to study. | +| `REFERENCE_TOLERANCE` | `1e-25` | Tolerance used to compute the reference ("exact") ERIs. | + +Each sampled shell is placed on all four corners of a regular tetrahedron +whose edge length is the separation, so the number of shells, and therefore of +shell quartets, grows quickly with the number of samples. diff --git a/examples/eri_error_explorer/eri_error_explorer.py b/examples/eri_error_explorer/eri_error_explorer.py new file mode 100644 index 00000000..9d5f2f85 --- /dev/null +++ b/examples/eri_error_explorer/eri_error_explorer.py @@ -0,0 +1,340 @@ +# +# Copyright 2026 NWChemEx-Project +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +"""Explores how screening introduces error into four-center ERIs. + +This example looks at how three knobs impact the error in the four-center +electron repulsion integrals (ERIs): + +1. The basis set parameters, i.e., which STO-3G shells are used. +2. The separation between the centers the shells sit on. +3. The screening tolerance (the "Threshold" input of the "ERI4" module). + +The sampled shells are placed on every vertex of a regular tetrahedron whose +edge length is the center separation, i.e., every pair of centers is the same +distance apart. The ERIs are then computed one shell quartet at a time, only +for symmetry-unique quartets. Every quartet is computed once at +``REFERENCE_TOLERANCE``, which defines the "exact" ERIs, and once per +screening tolerance. The error is the difference between the screened and the +reference ERIs. + +Exploring a different region of parameter space is done by changing the +hard-coded parameters below and re-running the module. +""" + +import math + +import chemcache +import integrals +import numpy as np +import parallelzone as pz +import pluginplay as pp +import simde +from chemist import Electron, PointD +from chemist.basis_set import AOBasisSetD, AtomicBasisSetD, ShellD +from chemist.braket import BraKet, _BraKet_ERI4 +from chemist.qm_operator import Coulomb +from chemist.wavefunction import AOs + +# ----------------------------------------------------------------------------- +# Parameters controlling the run +# ----------------------------------------------------------------------------- + +#: Shells to sample. Each entry is (Z, i) where i is the index of the shell in +#: ChemCache's STO-3G atomic basis set for the atom with atomic number Z. The +#: default is the H 1s and the O 2p shell. +SHELL_SAMPLES = [(1, 0), (8, 2)] + +#: The center separations (in bohr) are ``N_SEPARATIONS`` evenly spaced values +#: from ``R_MIN`` to ``R_MAX`` (inclusive). +R_MIN = 1.0 +R_MAX = 5.0 +N_SEPARATIONS = 3 + +#: Screening tolerances to study. +TOLERANCES = [1.0e-6, 1.0e-8] + +#: Screening tolerance used to compute the reference ("exact") ERIs. +REFERENCE_TOLERANCE = 1.0e-25 + +# ----------------------------------------------------------------------------- +# Implementation +# ----------------------------------------------------------------------------- + +_L_LABELS = "spdfghik" + +#: Number of centers, i.e., vertices of a tetrahedron. +_N_CENTERS = 4 + + +class SampledShell: + """The parameters of a shell taken from the STO-3G basis set. + + The center of a shell is not part of the shell's parameters, so that the + same parameters can be placed on any center. + + :param label: A human-readable name for the shell, e.g., "O[2] p". + :type label: str + :param Z: The atomic number of the atom the shell came from. + :type Z: int + :param am: The angular momentum of the shell. + :type am: int + :param pure: Whether the shell is pure or Cartesian. + :type pure: chemist.ShellType + :param coefs: The contraction coefficients of the shell. + :type coefs: list(float) + :param exps: The exponents of the primitives in the shell. + :type exps: list(float) + """ + + def __init__(self, label, Z, am, pure, coefs, exps): + self.label = label + self.Z = Z + self.am = am + self.pure = pure + self.coefs = coefs + self.exps = exps + + def at(self, center, center_label): + """Places the shell on ``center``. + + :param center: The (x, y, z) coordinates of the center. + :type center: tuple(float) + :param center_label: A name for the center, used for printing. + :type center_label: str + + :return: A single-shell AO space and a label for it. + :rtype: tuple(chemist.wavefunction.AOs, str) + """ + x, y, z = center + shell = ShellD(self.pure, self.am, self.coefs, self.exps, x, y, z) + atom = AtomicBasisSetD("sto-3g", self.Z, PointD(x, y, z), [shell]) + bs = AOBasisSetD() + bs.add_center(atom) + return AOs(bs), f"{self.label}@{center_label}" + + +def get_sto3g_shell(mm, Z, i): + """Pulls the parameters of shell ``i`` of atom ``Z`` out of ChemCache. + + :param mm: A ModuleManager with ChemCache loaded. + :type mm: pluginplay.ModuleManager + :param Z: The atomic number of the atom. + :type Z: int + :param i: The index of the shell in the atom's STO-3G basis set. + :type i: int + + :return: The parameters of the requested shell. + :rtype: SampledShell + """ + atom = mm.run_as(simde.AtomicBasisSetFromZ(), "sto-3g atomic basis", Z) + shell = atom.at(i) + prims = [shell.primitive(p) for p in range(shell.n_primitives())] + coefs = [prim.coefficient for prim in prims] + exps = [prim.exponent for prim in prims] + symbol = mm.run_as(simde.SymbolFromZ(), "Symbol from Z", Z) + label = f"{symbol}[{i}] {_L_LABELS[shell.l]}" + return SampledShell(label, Z, shell.l, shell.pure, coefs, exps) + + +def tetrahedron_centers(r): + """The vertices of a regular tetrahedron with edge length ``r``. + + :param r: The edge length of the tetrahedron. + :type r: float + + :return: The (x, y, z) coordinates of the four vertices. + :rtype: list(tuple(float)) + """ + scale = r / (2.0 * math.sqrt(2.0)) + signs = [(1, 1, 1), (1, -1, -1), (-1, 1, -1), (-1, -1, 1)] + return [tuple(scale * s for s in sign) for sign in signs] + + +def unique_quartets(n): + """Yields the symmetry-unique shell quartets for ``n`` shells. + + A quartet (i, j, k, l) is symmetry unique if i >= j, k >= l, and the + pair index of (i, j) is not less than that of (k, l). + + :param n: The number of shells. + :type n: int + + :return: A generator over the unique quartets. + :rtype: generator(tuple(int)) + """ + pairs = [(i, j) for i in range(n) for j in range(i + 1)] + for ij, (i, j) in enumerate(pairs): + for k, m in pairs[: ij + 1]: + yield i, j, k, m + + +def make_eri4_modules(mm): + """Makes one copy of the "ERI4" module per screening tolerance. + + Modules are locked once they run, so rather than changing the threshold + of "ERI4" between runs, each tolerance gets its own copy. + + :param mm: A ModuleManager with Integrals loaded. + :type mm: pluginplay.ModuleManager + + :return: The key of the reference module and a dictionary from each + tolerance to the key of its module. + :rtype: tuple(str, dict(float, str)) + """ + ref_key = "ERI4 reference" + mm.copy_module("ERI4", ref_key) + mm.change_input(ref_key, "Threshold", REFERENCE_TOLERANCE) + + tol_keys = {} + for tol in TOLERANCES: + tol_keys[tol] = f"ERI4 tol={tol:.1e}" + mm.copy_module("ERI4", tol_keys[tol]) + mm.change_input(tol_keys[tol], "Threshold", tol) + return ref_key, tol_keys + + +def compute_eri4(mm, key, aos): + """Computes the ERIs for one shell quartet. + + :param mm: A ModuleManager with Integrals loaded. + :type mm: pluginplay.ModuleManager + :param key: The key of the module to use. + :type key: str + :param aos: The single-shell AO spaces of the quartet. + :type aos: tuple(chemist.wavefunction.AOs) + + :return: The (n_i, n_j, n_k, n_l) block of ERIs. + :rtype: numpy.ndarray + """ + v_ee = Coulomb(Electron(), Electron()) + braket = BraKet((aos[0], aos[1]), v_ee, (aos[2], aos[3])) + # N.B. C++ modules currently fail to convert the erased BraKet to the + # type ERI4 wants, so we hand them the concrete type. + t = mm.run_as(simde.ERI4(), key, _BraKet_ERI4(braket)) + return np.array(t) + + +def print_shells(shells): + """Prints the parameters of the sampled shells. + + :param shells: The sampled shells. + :type shells: list(SampledShell) + """ + print("Sampled STO-3G shells:") + for shell in shells: + print(f" {shell.label} (Z = {shell.Z}, l = {shell.am}, {shell.pure})") + print(" exponents :", " ".join(f"{e:.10e}" for e in shell.exps)) + print(" coefficients :", " ".join(f"{c:.10e}" for c in shell.coefs)) + print() + + +def print_summary(r, tol, records): + """Summarizes the errors for one separation and tolerance. + + :param r: The center separation. + :type r: float + :param tol: The screening tolerance. + :type tol: float + :param records: For each quartet, a tuple of the quartet's label, the + number of distinct centers, and the absolute error block. + :type records: list(tuple(str, int, numpy.ndarray)) + """ + print(f"Summary for R = {r:.4f}, tol = {tol:.1e}") + header = ( + f" {'centers':>7} {'quartets':>8} {'max|err|':>12} {'sum|err|':>12}" + f" worst quartet" + ) + print(header) + for n_centers in range(1, _N_CENTERS + 1): + subset = [rec for rec in records if rec[1] == n_centers] + if not subset: + continue + worst = max(subset, key=lambda rec: rec[2].max()) + total = sum(rec[2].sum() for rec in subset) + print( + f" {n_centers:>7} {len(subset):>8} {worst[2].max():12.4e}" + f" {total:12.4e} {worst[0]}" + ) + print() + + +def main(): + rv = pz.runtime.RuntimeView() + + # N.b. passing None for the cache turns memoization off + mm = pp.ModuleManager(rv, None) + mm.load_plugin(chemcache.ChemCache()) + mm.load_plugin(integrals.Integrals()) + + ref_key, tol_keys = make_eri4_modules(mm) + shells = [get_sto3g_shell(mm, Z, i) for Z, i in SHELL_SAMPLES] + print_shells(shells) + + row = "{:>8} {:>8} {:>3} {:>12} {:>12} {:>12} {}" + header = ( + "R", + "tol", + "nc", + "max|ERI|", + "max|err|", + "max|err|/tol", + "quartet", + ) + for r in np.linspace(R_MIN, R_MAX, N_SEPARATIONS): + centers = tetrahedron_centers(r) + + # Each element is (AOs, label, center index) + basis = [] + for c, center in enumerate(centers): + for shell in shells: + aos, label = shell.at(center, f"C{c}") + basis.append((aos, label, c)) + + print(f"R = {r:.4f} bohr ({len(basis)} shells)") + print(row.format(*header)) + + records = {tol: [] for tol in TOLERANCES} + for quartet in unique_quartets(len(basis)): + aos = tuple(basis[s][0] for s in quartet) + labels = [basis[s][1] for s in quartet] + label = f"({labels[0]} {labels[1]}|{labels[2]} {labels[3]})" + n_centers = len({basis[s][2] for s in quartet}) + + ref = compute_eri4(mm, ref_key, aos) + max_eri = np.abs(ref).max() + for tol, key in tol_keys.items(): + error = np.abs(compute_eri4(mm, key, aos) - ref) + records[tol].append((label, n_centers, error)) + print( + row.format( + f"{r:.4f}", + f"{tol:.1e}", + n_centers, + f"{max_eri:.4e}", + f"{error.max():.4e}", + f"{error.max() / tol:.4e}", + label, + ) + ) + print() + + for tol in TOLERANCES: + print_summary(r, tol, records[tol]) + + +if __name__ == "__main__": + main()