From df2111605479b2ba34f018a84593f560068e8ff1 Mon Sep 17 00:00:00 2001 From: kir943 Date: Sun, 15 Feb 2026 02:33:08 +0530 Subject: [PATCH 01/19] =?UTF-8?q?Add=20N-center=20Gaussian=20overlap=20eng?= =?UTF-8?q?ine=20with=20screening=20and=20sparse=20builder=20(Week=201?= =?UTF-8?q?=E2=80=932)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- gbasis/integrals/overlap_n.py | 308 ++++++++++++++++++++++++++++++++++ test_n_engine_stress.py | 115 +++++++++++++ test_n_sparse_builder.py | 36 ++++ test_overlap_n.py | 51 ++++++ 4 files changed, 510 insertions(+) create mode 100644 gbasis/integrals/overlap_n.py create mode 100644 test_n_engine_stress.py create mode 100644 test_n_sparse_builder.py create mode 100644 test_overlap_n.py diff --git a/gbasis/integrals/overlap_n.py b/gbasis/integrals/overlap_n.py new file mode 100644 index 00000000..2c453256 --- /dev/null +++ b/gbasis/integrals/overlap_n.py @@ -0,0 +1,308 @@ +import numpy as np +from scipy.sparse import coo_matrix +from itertools import product +class PrimitiveNEngine: + + # Gaussian collapse + @staticmethod + def collapse_gaussians(alphas, centers): + alphas = np.asarray(alphas, dtype=np.float64) + centers = np.asarray(centers, dtype=np.float64) + + alpha_tot = np.sum(alphas) + + if alpha_tot <= 0.0: + raise ValueError("Total Gaussian exponent must be positive.") + + P = np.sum(alphas[:, None] * centers, axis=0) / alpha_tot + + term1 = np.sum(alphas * np.sum(centers**2, axis=1)) + term2 = alpha_tot * np.dot(P, P) + + exponent = term2 - term1 + prefactor = np.exp(exponent) + + return alpha_tot, P, prefactor + # Pure binomial Hermite shift + @staticmethod + def hermite_coefficients(l, PA): + """ + Expand (x - A)^l about P: + + (x - A)^l = sum_t E_t (x - P)^t + """ + E = np.zeros(l + 1, dtype=np.float64) + E[0] = 1.0 + + for i in range(l): + E_new = np.zeros(l + 1, dtype=np.float64) + + for t in range(i + 1): + E_new[t] += PA * E[t] + E_new[t + 1] += E[t] + + E = E_new + + return E + # Gaussian moments + @staticmethod + def gaussian_moments(alpha, max_order): + """ + Compute: + ∫ (x-P)^k exp(-alpha (x-P)^2) dx + over (-∞,∞) + """ + moments = np.zeros(max_order + 1, dtype=np.float64) + + # zeroth moment + moments[0] = np.sqrt(np.pi / alpha) + + # only even moments survive + for k in range(0, max_order - 1, 2): + moments[k + 2] = (k + 1) / (2.0 * alpha) * moments[k] + + return moments + # Full primitive N-center overlap + @staticmethod + def primitive_overlap(alphas, centers, angmoms): + + alpha_tot, P, prefactor = PrimitiveNEngine.collapse_gaussians( + alphas, centers + ) + + result = prefactor + + # factorize into x, y, z + for axis in range(3): + + # build total polynomial via convolution + E_total = np.array([1.0], dtype=np.float64) + + for i in range(len(alphas)): + l = angmoms[i][axis] + PA = P[axis] - centers[i][axis] + + E = PrimitiveNEngine.hermite_coefficients(l, PA) + + E_total = np.convolve(E_total, E) + + moments = PrimitiveNEngine.gaussian_moments( + alpha_tot, + len(E_total) - 1 + ) + + axis_integral = np.dot(E_total, moments[:len(E_total)]) + + result *= axis_integral + + return result +# Screening Function +def is_n_shell_overlap_screened(shells, tol=1e-12): + """ + Conservative exponential upper-bound screening + for N-center contracted overlap. + """ + + alpha_mins = [np.min(shell.exps) for shell in shells] + centers = [shell.coord for shell in shells] + + alpha_tot = sum(alpha_mins) + + if alpha_tot <= 0.0: + return True + + # Exponential decay from Gaussian collapse + decay_sum = 0.0 + N = len(shells) + + for i in range(N): + for j in range(i + 1, N): + Rij = centers[i] - centers[j] + Rij2 = np.dot(Rij, Rij) + decay_sum += alpha_mins[i] * alpha_mins[j] * Rij2 + + D = decay_sum / alpha_tot + + # Contraction-level magnitude bound + coeff_bound = 1.0 + norm_bound = 1.0 + + for shell in shells: + coeff_bound *= np.max(np.abs(shell.coeffs)) + norm_bound *= np.max(np.abs(shell.norm_prim_cart)) + + volume_bound = (np.pi / alpha_tot) ** 1.5 + + bound = coeff_bound * norm_bound * volume_bound * np.exp(-D) + + return bound < tol +def contracted_n_overlap(shells): + """ + Compute contracted N-center overlap for a list of + GeneralizedContractionShell objects. + + Parameters + ---------- + shells : list[GeneralizedContractionShell] + + Returns + ------- + np.ndarray + N-dimensional array over segmented contractions + and Cartesian angular components. + + Shape: + (M1, L1, M2, L2, ..., MN, LN) + """ + + N = len(shells) + + # Build shape for final tensor + shape = [] + for shell in shells: + shape.append(shell.num_seg_cont) + shape.append(shell.num_cart) + + result = np.zeros(shape, dtype=np.float64) + + # Primitive exponent index ranges + prim_ranges = [range(len(shell.exps)) for shell in shells] + + # Segmented contraction indices + seg_ranges = [range(shell.num_seg_cont) for shell in shells] + + # Cartesian angular component indices + cart_ranges = [range(shell.num_cart) for shell in shells] + + for seg_indices in product(*seg_ranges): + for cart_indices in product(*cart_ranges): + + total_value = 0.0 + for prim_indices in product(*prim_ranges): + + alphas = [] + centers = [] + angmoms = [] + coeff_prod = 1.0 + norm_prod = 1.0 + + for i, shell in enumerate(shells): + + p = prim_indices[i] + m = seg_indices[i] + c = cart_indices[i] + + alpha = shell.exps[p] + coeff = shell.coeffs[p, m] + norm = shell.norm_prim_cart[c, p] + angmom = tuple(shell.angmom_components_cart[c]) + + alphas.append(alpha) + centers.append(shell.coord) + angmoms.append(angmom) + + coeff_prod *= coeff + norm_prod *= norm + + prim_val = PrimitiveNEngine.primitive_overlap( + alphas, + centers, + angmoms + ) + + total_value += coeff_prod * norm_prod * prim_val + + index = [] + for i in range(N): + index.append(seg_indices[i]) + index.append(cart_indices[i]) + + result[tuple(index)] = total_value + + return result +def build_n_overlap_tensor(shells, tol=1e-12): + """ + Build sparse N-center overlap tensor over shells. + + Parameters + ---------- + shells : list[GeneralizedContractionShell] + tol : float + Screening tolerance + + Returns + ------- + scipy.sparse.coo_matrix + Flattened sparse tensor of shape (total_ao^N, 1) + """ + + # Total AO dimension + shell_sizes = [ + shell.num_seg_cont * shell.num_cart + for shell in shells + ] + + total_ao = sum(shell_sizes) + N = len(shells) + + data = [] + rows = [] + + # compute AO offsets per shell + offsets = [] + acc = 0 + for size in shell_sizes: + offsets.append(acc) + acc += size + + for shell_indices in product(range(len(shells)), repeat=N): + + shell_tuple = [shells[i] for i in shell_indices] + + # Screening + if is_n_shell_overlap_screened(shell_tuple, tol=tol): + continue + + block = contracted_n_overlap(shell_tuple) + + block_flat = block.reshape(-1) + + local_sizes = [ + shells[i].num_seg_cont * shells[i].num_cart + for i in shell_indices + ] + + local_offsets = [ + offsets[i] + for i in shell_indices + ] + + for local_idx, value in enumerate(block_flat): + + if abs(value) < tol: + continue + + # convert local multi-index to global index + multi = [] + tmp = local_idx + + for size in reversed(local_sizes): + multi.append(tmp % size) + tmp //= size + + multi = list(reversed(multi)) + + global_index = 0 + for k in range(N): + global_index = ( + global_index * total_ao + + local_offsets[k] + + multi[k] + ) + + rows.append(global_index) + data.append(value) + + shape = (total_ao ** N, 1) + + return coo_matrix((data, (rows, np.zeros(len(rows)))), shape=shape) \ No newline at end of file diff --git a/test_n_engine_stress.py b/test_n_engine_stress.py new file mode 100644 index 00000000..76ebfebf --- /dev/null +++ b/test_n_engine_stress.py @@ -0,0 +1,115 @@ +import numpy as np +import pytest + +from gbasis.parsers import parse_gbs, make_contractions +from gbasis.integrals.overlap import overlap_integral +from gbasis.integrals.overlap_n import contracted_n_overlap + + +# ===================================================== +# Helper: Build small hydrogen basis +# ===================================================== +def build_h2_basis(): + atoms = ["H", "H"] + coords = np.array([[0.0, 0.0, 0.0], + [0.0, 0.0, 1.0]]) + + basis_dict = parse_gbs("tests/data_631g.gbs") + + basis = make_contractions( + basis_dict, + atoms, + coords, + coord_types="c" + ) + + return basis + + +# ===================================================== +# 1️⃣ N=2 FULL AO MATRIX VALIDATION +# ===================================================== + +def test_n2_matches_gbasis(): + + basis = build_h2_basis() + + ref = overlap_integral(basis, screen_basis=False) + + n_shells = len(basis) + total_ao = ref.shape[0] + mine = np.zeros((total_ao, total_ao)) + + ao_i = 0 + for shell_i in basis: + size_i = shell_i.num_seg_cont * shell_i.num_cart + + ao_j = 0 + for shell_j in basis: + size_j = shell_j.num_seg_cont * shell_j.num_cart + + block = contracted_n_overlap([shell_i, shell_j]) + block = block.reshape(size_i, size_j) + + mine[ao_i:ao_i+size_i, + ao_j:ao_j+size_j] = block + + ao_j += size_j + + ao_i += size_i + + diff = np.max(np.abs(ref - mine)) + + print("N=2 max difference:", diff) + + assert np.allclose(ref, mine, atol=1e-10) + + +# ===================================================== +# 2️⃣ N=3 SYMMETRY TEST +# ===================================================== + +def test_n3_symmetry(): + + basis = build_h2_basis() + + shell1 = basis[0] + shell2 = basis[1] + shell3 = basis[0] + + S123 = contracted_n_overlap([shell1, shell2, shell3]) + S321 = contracted_n_overlap([shell3, shell2, shell1]) + + diff = np.max(np.abs(S123 - S321)) + + print("N=3 symmetry difference:", diff) + + assert np.allclose(S123, S321, atol=1e-10) + + +# ===================================================== +# 3️⃣ HIGH ANGULAR MOMENTUM + DIFFUSE TEST +# ===================================================== + +def test_high_angmom_diffuse(): + + basis = build_h2_basis() + + # choose a higher angular momentum shell if available + # try to find p-shell + p_shells = [shell for shell in basis if shell.angmom == 1] + + if not p_shells: + pytest.skip("No p-shell available in basis.") + + shell = p_shells[0] + + # artificially make exponents diffuse + shell.exps = shell.exps * 0.01 + + S = contracted_n_overlap([shell, shell]) + + print("Diffuse shell max abs value:", np.max(np.abs(S))) + + assert not np.isnan(S).any() + assert not np.isinf(S).any() diff --git a/test_n_sparse_builder.py b/test_n_sparse_builder.py new file mode 100644 index 00000000..3c689ece --- /dev/null +++ b/test_n_sparse_builder.py @@ -0,0 +1,36 @@ +import numpy as np +from pathlib import Path + +from gbasis.parsers import parse_gbs, make_contractions +from gbasis.integrals.overlap_n import build_sparse_n_overlap_tensor + + +def build_h2_basis(): + atoms = ["H", "H"] + coords = np.array([[0.0, 0.0, 0.0], + [0.0, 0.0, 1.0]]) + + test_dir = Path(__file__).parent + basis_path = test_dir / "data" / "hydrogen_def2-svp.1.gbs" + + basis_dict = parse_gbs(str(basis_path)) + basis = make_contractions( + basis_dict, + atoms, + coords, + coord_types="c" + ) + return basis + + +def test_sparse_builder_runs(): + + basis = build_h2_basis() + + sparse_mat = build_sparse_n_overlap_tensor(basis) + + print("Sparse tensor shape:", sparse_mat.shape) + print("Number of nonzeros:", sparse_mat.nnz) + + assert sparse_mat.shape[0] > 0 + assert sparse_mat.nnz > 0 diff --git a/test_overlap_n.py b/test_overlap_n.py new file mode 100644 index 00000000..a2155a3c --- /dev/null +++ b/test_overlap_n.py @@ -0,0 +1,51 @@ +import numpy as np +import pytest + +from gbasis.integrals.overlap import Overlap +from gbasis.integrals.overlap_n import contracted_n_overlap, PrimitiveNEngine +from gbasis.contractions import GeneralizedContractionShell + + +def build_test_shell(center): + ... + return GeneralizedContractionShell(...) + + +def test_two_center_matches_existing_overlap(): + ... + assert ... + + +def test_two_center_ss_analytic(): + ... + assert ... + + +def test_primitive_px_py_pz(): + A = np.array([0.0, 0.0, 0.0]) + B = np.array([0.3, 0.1, -0.2]) + + alpha = 0.5 + + print("px-px:", + PrimitiveNEngine.primitive_overlap( + [alpha, alpha], + [A, B], + [(1,0,0), (1,0,0)] + )) + + print("py-py:", + PrimitiveNEngine.primitive_overlap( + [alpha, alpha], + [A, B], + [(0,1,0), (0,1,0)] + )) + + print("pz-pz:", + PrimitiveNEngine.primitive_overlap( + [alpha, alpha], + [A, B], + [(0,0,1), (0,0,1)] + )) + + assert True # temporary so pytest runs it From 15bb7f757454168d6e724b3af08e7401a8c09886 Mon Sep 17 00:00:00 2001 From: kir943 Date: Sun, 15 Feb 2026 02:55:51 +0530 Subject: [PATCH 02/19] =?UTF-8?q?Remove=20experimental=20sparse=20builder?= =?UTF-8?q?=20test=20(not=20part=20of=20Week=201=E2=80=932)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- test_n_sparse_builder.py | 36 ------------------------------------ 1 file changed, 36 deletions(-) delete mode 100644 test_n_sparse_builder.py diff --git a/test_n_sparse_builder.py b/test_n_sparse_builder.py deleted file mode 100644 index 3c689ece..00000000 --- a/test_n_sparse_builder.py +++ /dev/null @@ -1,36 +0,0 @@ -import numpy as np -from pathlib import Path - -from gbasis.parsers import parse_gbs, make_contractions -from gbasis.integrals.overlap_n import build_sparse_n_overlap_tensor - - -def build_h2_basis(): - atoms = ["H", "H"] - coords = np.array([[0.0, 0.0, 0.0], - [0.0, 0.0, 1.0]]) - - test_dir = Path(__file__).parent - basis_path = test_dir / "data" / "hydrogen_def2-svp.1.gbs" - - basis_dict = parse_gbs(str(basis_path)) - basis = make_contractions( - basis_dict, - atoms, - coords, - coord_types="c" - ) - return basis - - -def test_sparse_builder_runs(): - - basis = build_h2_basis() - - sparse_mat = build_sparse_n_overlap_tensor(basis) - - print("Sparse tensor shape:", sparse_mat.shape) - print("Number of nonzeros:", sparse_mat.nnz) - - assert sparse_mat.shape[0] > 0 - assert sparse_mat.nnz > 0 From 76d843da951af505dee8c11e3a4312d01db3151e Mon Sep 17 00:00:00 2001 From: kir943 Date: Fri, 20 Feb 2026 22:16:40 +0530 Subject: [PATCH 03/19] Move overlap test files to tests directory as requested in review --- gbasis/integrals/multi_overlap.py | 32 ++++++++++ integrals/overlap.py | 0 test_compare_overlap.py | 64 +++++++++++++++++++ test_multioverlap.py | 21 ++++++ test_overlap.py | 47 ++++++++++++++ test_try.py | 41 ++++++++++++ .../test_n_engine_stress.py | 0 test_overlap_n.py => tests/test_overlap_n.py | 0 8 files changed, 205 insertions(+) create mode 100644 gbasis/integrals/multi_overlap.py create mode 100644 integrals/overlap.py create mode 100644 test_compare_overlap.py create mode 100644 test_multioverlap.py create mode 100644 test_overlap.py create mode 100644 test_try.py rename test_n_engine_stress.py => tests/test_n_engine_stress.py (100%) rename test_overlap_n.py => tests/test_overlap_n.py (100%) diff --git a/gbasis/integrals/multi_overlap.py b/gbasis/integrals/multi_overlap.py new file mode 100644 index 00000000..cb363029 --- /dev/null +++ b/gbasis/integrals/multi_overlap.py @@ -0,0 +1,32 @@ +import numpy as np +from itertools import product + + +def three_gaussian_overlap_primitive(prims): + alphas = np.array([p[0] for p in prims]) + centers = np.array([p[1] for p in prims]) + + alpha_tot = np.sum(alphas) + P = np.sum(alphas[:, None] * centers, axis=0) / alpha_tot + + term1 = np.sum(alphas * np.sum(centers**2, axis=1)) + term2 = alpha_tot * np.dot(P, P) + + prefactor = np.exp(-(term1 - term2)) + return prefactor * (np.pi / alpha_tot) ** 1.5 + + +def three_overlap_tensor(basis): + """ + basis: list of primitives [(alpha, center)] + returns T[μ, ν, λ] + """ + n = len(basis) + T = np.zeros((n, n, n)) + + for μ, ν, λ in product(range(n), repeat=3): + T[μ, ν, λ] = three_gaussian_overlap_primitive( + [basis[μ], basis[ν], basis[λ]] + ) + + return T diff --git a/integrals/overlap.py b/integrals/overlap.py new file mode 100644 index 00000000..e69de29b diff --git a/test_compare_overlap.py b/test_compare_overlap.py new file mode 100644 index 00000000..1e343d64 --- /dev/null +++ b/test_compare_overlap.py @@ -0,0 +1,64 @@ +import numpy as np +from gbasis.parsers import parse_gbs +from gbasis.wrappers import make_contractions +from gbasis.integrals.overlap import Overlap +from gbasis.integrals.overlap_n import contracted_n_overlap + +atoms = ["H", "H"] +coords = np.array([[0.0, 0.0, 0.0], + [0.0, 0.0, 1.0]]) + +basis_dict = parse_gbs("tests/data_631g.gbs") +basis = make_contractions(basis_dict, atoms, coords, coord_types="c") + +print("Number of shells:", len(basis)) + +# Reference overlap (NO screening) +olp_ref = Overlap.construct_array(basis) + + +print("Reference overlap shape:", olp_ref.shape) + +n_shells = len(basis) + +total_ao = olp_ref.shape[0] +olp_mine = np.zeros((total_ao, total_ao)) + +ao_offset_i = 0 + +for i, shell_i in enumerate(basis): + + Li = shell_i.num_cart + Mi = shell_i.num_seg_cont + size_i = Li * Mi + + ao_offset_j = 0 + + for j, shell_j in enumerate(basis): + + Lj = shell_j.num_cart + Mj = shell_j.num_seg_cont + size_j = Lj * Mj + + # Compute N=2 overlap + block = contracted_n_overlap([shell_i, shell_j]) + + block = block.reshape(size_i, size_j) + + olp_mine[ + ao_offset_i:ao_offset_i+size_i, + ao_offset_j:ao_offset_j+size_j + ] = block + + ao_offset_j += size_j + + ao_offset_i += size_i + + +print("My overlap shape:", olp_mine.shape) +# Compare +diff = np.max(np.abs(olp_ref - olp_mine)) + +print("Maximum absolute difference:", diff) +print("All close (1e-10)?", np.allclose(olp_ref, olp_mine, atol=1e-10)) + diff --git a/test_multioverlap.py b/test_multioverlap.py new file mode 100644 index 00000000..445d19f3 --- /dev/null +++ b/test_multioverlap.py @@ -0,0 +1,21 @@ +import numpy as np +from gbasis.integrals.multi_overlap import three_overlap_tensor + + +def test_tensor(): + alpha = 1.0 + basis = [ + (alpha, np.array([0.0, 0.0, 0.0])), + (alpha, np.array([0.0, 0.0, 1.0])), + (alpha, np.array([0.0, 1.0, 0.0])), + ] + + T = three_overlap_tensor(basis) + + print("Tensor shape:", T.shape) + print("T[0,0,0] =", T[0,0,0]) + print("T[0,1,2] =", T[0,1,2]) + + +if __name__ == "__main__": + test_tensor() diff --git a/test_overlap.py b/test_overlap.py new file mode 100644 index 00000000..c32fdf2b --- /dev/null +++ b/test_overlap.py @@ -0,0 +1,47 @@ +# Testing for 1-2Gaussian overlapimport pytest + +def test_overlap(): + assert 1 + 1 == 2 + +import numpy as np + +from gbasis.contractions import GeneralizedContractionShell +from gbasis.integrals.overlap import overlap_integral + + +# ------------------------- +# Define two 1s Gaussian shells +# ------------------------- +angmom = 0 # s orbital + +exps = np.array([0.5]) +coeffs = np.array([[1.0]]) + +g1 = GeneralizedContractionShell( + angmom=angmom, + coord=np.array([0.0, 0.0, 0.0]), + exps=exps, + coeffs=coeffs, + coord_type="cartesian", +) + +g2 = GeneralizedContractionShell( + angmom=angmom, + coord=np.array([0.2, 0.0, 0.0]), + exps=exps, + coeffs=coeffs, + coord_type="cartesian", +) + +# ------------------------- +# BASIS = list of shells +# ------------------------- +basis = [g1, g2] + + +# Overlap matrix + +S = overlap_integral(basis) + +print("Overlap matrix:") +print(S) diff --git a/test_try.py b/test_try.py new file mode 100644 index 00000000..d8a05508 --- /dev/null +++ b/test_try.py @@ -0,0 +1,41 @@ +import numpy as np +from pathlib import Path +from gbasis.parsers import parse_gbs, make_contractions +from gbasis.integrals.overlap_n import build_n_overlap_tensor + + +def build_h2_basis(): + atoms = ["H", "H"] + coords = np.array([ + [0.0, 0.0, 0.0], + [0.0, 0.0, 1.0] + ]) + + basis_path = Path( + r"C:\Users\kiran\OneDrive\ドキュメント\nnn\gbasis\notebooks\tutorial\hydrogen_def2-svp.1.gbs" + ) + + basis_dict = parse_gbs(str(basis_path)) + + return make_contractions( + basis_dict, + atoms, + coords, + coord_types="c" + ) + + +def test_n_overlap_tensor_builder(): + + basis = build_h2_basis() + + shell1 = basis[0] + shell2 = basis[1] + + tensor = build_n_overlap_tensor([shell1, shell2]) + + print("Tensor shape:", tensor.shape) + print("Nonzeros:", tensor.nnz) + + assert tensor.shape[0] > 0 + assert tensor.nnz > 0 diff --git a/test_n_engine_stress.py b/tests/test_n_engine_stress.py similarity index 100% rename from test_n_engine_stress.py rename to tests/test_n_engine_stress.py diff --git a/test_overlap_n.py b/tests/test_overlap_n.py similarity index 100% rename from test_overlap_n.py rename to tests/test_overlap_n.py From ff83f68bb0d7e0551b255e79940f34463e39bce5 Mon Sep 17 00:00:00 2001 From: kir943 Date: Fri, 20 Feb 2026 22:21:41 +0530 Subject: [PATCH 04/19] Remove unrelated/extra test files from PR --- test_compare_overlap.py | 64 ----------------------------------------- test_multioverlap.py | 21 -------------- test_overlap.py | 47 ------------------------------ test_try.py | 41 -------------------------- 4 files changed, 173 deletions(-) delete mode 100644 test_compare_overlap.py delete mode 100644 test_multioverlap.py delete mode 100644 test_overlap.py delete mode 100644 test_try.py diff --git a/test_compare_overlap.py b/test_compare_overlap.py deleted file mode 100644 index 1e343d64..00000000 --- a/test_compare_overlap.py +++ /dev/null @@ -1,64 +0,0 @@ -import numpy as np -from gbasis.parsers import parse_gbs -from gbasis.wrappers import make_contractions -from gbasis.integrals.overlap import Overlap -from gbasis.integrals.overlap_n import contracted_n_overlap - -atoms = ["H", "H"] -coords = np.array([[0.0, 0.0, 0.0], - [0.0, 0.0, 1.0]]) - -basis_dict = parse_gbs("tests/data_631g.gbs") -basis = make_contractions(basis_dict, atoms, coords, coord_types="c") - -print("Number of shells:", len(basis)) - -# Reference overlap (NO screening) -olp_ref = Overlap.construct_array(basis) - - -print("Reference overlap shape:", olp_ref.shape) - -n_shells = len(basis) - -total_ao = olp_ref.shape[0] -olp_mine = np.zeros((total_ao, total_ao)) - -ao_offset_i = 0 - -for i, shell_i in enumerate(basis): - - Li = shell_i.num_cart - Mi = shell_i.num_seg_cont - size_i = Li * Mi - - ao_offset_j = 0 - - for j, shell_j in enumerate(basis): - - Lj = shell_j.num_cart - Mj = shell_j.num_seg_cont - size_j = Lj * Mj - - # Compute N=2 overlap - block = contracted_n_overlap([shell_i, shell_j]) - - block = block.reshape(size_i, size_j) - - olp_mine[ - ao_offset_i:ao_offset_i+size_i, - ao_offset_j:ao_offset_j+size_j - ] = block - - ao_offset_j += size_j - - ao_offset_i += size_i - - -print("My overlap shape:", olp_mine.shape) -# Compare -diff = np.max(np.abs(olp_ref - olp_mine)) - -print("Maximum absolute difference:", diff) -print("All close (1e-10)?", np.allclose(olp_ref, olp_mine, atol=1e-10)) - diff --git a/test_multioverlap.py b/test_multioverlap.py deleted file mode 100644 index 445d19f3..00000000 --- a/test_multioverlap.py +++ /dev/null @@ -1,21 +0,0 @@ -import numpy as np -from gbasis.integrals.multi_overlap import three_overlap_tensor - - -def test_tensor(): - alpha = 1.0 - basis = [ - (alpha, np.array([0.0, 0.0, 0.0])), - (alpha, np.array([0.0, 0.0, 1.0])), - (alpha, np.array([0.0, 1.0, 0.0])), - ] - - T = three_overlap_tensor(basis) - - print("Tensor shape:", T.shape) - print("T[0,0,0] =", T[0,0,0]) - print("T[0,1,2] =", T[0,1,2]) - - -if __name__ == "__main__": - test_tensor() diff --git a/test_overlap.py b/test_overlap.py deleted file mode 100644 index c32fdf2b..00000000 --- a/test_overlap.py +++ /dev/null @@ -1,47 +0,0 @@ -# Testing for 1-2Gaussian overlapimport pytest - -def test_overlap(): - assert 1 + 1 == 2 - -import numpy as np - -from gbasis.contractions import GeneralizedContractionShell -from gbasis.integrals.overlap import overlap_integral - - -# ------------------------- -# Define two 1s Gaussian shells -# ------------------------- -angmom = 0 # s orbital - -exps = np.array([0.5]) -coeffs = np.array([[1.0]]) - -g1 = GeneralizedContractionShell( - angmom=angmom, - coord=np.array([0.0, 0.0, 0.0]), - exps=exps, - coeffs=coeffs, - coord_type="cartesian", -) - -g2 = GeneralizedContractionShell( - angmom=angmom, - coord=np.array([0.2, 0.0, 0.0]), - exps=exps, - coeffs=coeffs, - coord_type="cartesian", -) - -# ------------------------- -# BASIS = list of shells -# ------------------------- -basis = [g1, g2] - - -# Overlap matrix - -S = overlap_integral(basis) - -print("Overlap matrix:") -print(S) diff --git a/test_try.py b/test_try.py deleted file mode 100644 index d8a05508..00000000 --- a/test_try.py +++ /dev/null @@ -1,41 +0,0 @@ -import numpy as np -from pathlib import Path -from gbasis.parsers import parse_gbs, make_contractions -from gbasis.integrals.overlap_n import build_n_overlap_tensor - - -def build_h2_basis(): - atoms = ["H", "H"] - coords = np.array([ - [0.0, 0.0, 0.0], - [0.0, 0.0, 1.0] - ]) - - basis_path = Path( - r"C:\Users\kiran\OneDrive\ドキュメント\nnn\gbasis\notebooks\tutorial\hydrogen_def2-svp.1.gbs" - ) - - basis_dict = parse_gbs(str(basis_path)) - - return make_contractions( - basis_dict, - atoms, - coords, - coord_types="c" - ) - - -def test_n_overlap_tensor_builder(): - - basis = build_h2_basis() - - shell1 = basis[0] - shell2 = basis[1] - - tensor = build_n_overlap_tensor([shell1, shell2]) - - print("Tensor shape:", tensor.shape) - print("Nonzeros:", tensor.nnz) - - assert tensor.shape[0] > 0 - assert tensor.nnz > 0 From 2736394e5a259832394fbfc579feb6cf416f4791 Mon Sep 17 00:00:00 2001 From: kir943 Date: Fri, 20 Feb 2026 22:37:22 +0530 Subject: [PATCH 05/19] Remove incorrect duplicate integrals package --- integrals/overlap.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 integrals/overlap.py diff --git a/integrals/overlap.py b/integrals/overlap.py deleted file mode 100644 index e69de29b..00000000 From 9019e1514334ae0ef4602f42a4fb3c6168cc1671 Mon Sep 17 00:00:00 2001 From: kir943 Date: Fri, 20 Feb 2026 22:38:20 +0530 Subject: [PATCH 06/19] Remove multi_overlap.py (belongs to Week 3, not this PR) --- gbasis/integrals/multi_overlap.py | 32 ------------------------------- 1 file changed, 32 deletions(-) delete mode 100644 gbasis/integrals/multi_overlap.py diff --git a/gbasis/integrals/multi_overlap.py b/gbasis/integrals/multi_overlap.py deleted file mode 100644 index cb363029..00000000 --- a/gbasis/integrals/multi_overlap.py +++ /dev/null @@ -1,32 +0,0 @@ -import numpy as np -from itertools import product - - -def three_gaussian_overlap_primitive(prims): - alphas = np.array([p[0] for p in prims]) - centers = np.array([p[1] for p in prims]) - - alpha_tot = np.sum(alphas) - P = np.sum(alphas[:, None] * centers, axis=0) / alpha_tot - - term1 = np.sum(alphas * np.sum(centers**2, axis=1)) - term2 = alpha_tot * np.dot(P, P) - - prefactor = np.exp(-(term1 - term2)) - return prefactor * (np.pi / alpha_tot) ** 1.5 - - -def three_overlap_tensor(basis): - """ - basis: list of primitives [(alpha, center)] - returns T[μ, ν, λ] - """ - n = len(basis) - T = np.zeros((n, n, n)) - - for μ, ν, λ in product(range(n), repeat=3): - T[μ, ν, λ] = three_gaussian_overlap_primitive( - [basis[μ], basis[ν], basis[λ]] - ) - - return T From a1a43c397d63b1c124b4c925187c43097cca99da Mon Sep 17 00:00:00 2001 From: kir943 Date: Sun, 22 Feb 2026 21:30:12 +0530 Subject: [PATCH 07/19] Week 3: Add arbitrary order overlap tests (N=1 to N=6) --- tests/test_n.py | 57 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 57 insertions(+) create mode 100644 tests/test_n.py diff --git a/tests/test_n.py b/tests/test_n.py new file mode 100644 index 00000000..93cac243 --- /dev/null +++ b/tests/test_n.py @@ -0,0 +1,57 @@ +import numpy as np +import pytest + +from gbasis.parsers import parse_gbs, make_contractions +from gbasis.integrals.overlap_n import contracted_n_overlap + + +# ===================================================== +# Helper: Build hydrogen basis +# ===================================================== + +def build_h2_basis(): + + atoms = ["H", "H"] + + coords = np.array([ + [0.0, 0.0, 0.0], + [0.0, 0.0, 1.0] + ]) + + basis_dict = parse_gbs("tests/data_631g.gbs") + + basis = make_contractions( + basis_dict, + atoms, + coords, + coord_types="c" + ) + + return basis + + +# ===================================================== +# Test arbitrary orders N=1 to N=6 +# ===================================================== + +@pytest.mark.parametrize("N", [1, 2, 3, 4, 5, 6]) +def test_arbitrary_n_orders(N): + + basis = build_h2_basis() + + shells = [] + + for i in range(N): + shells.append(basis[i % len(basis)]) + + result = contracted_n_overlap(shells) + + print(f"N={N}, max abs value:", np.max(np.abs(result))) + + assert result is not None + + assert not np.isnan(result).any() + + assert not np.isinf(result).any() + + assert np.all(np.isfinite(result)) \ No newline at end of file From 07791e23fdd6b960d0c8181005f6c127f0d1bfa4 Mon Sep 17 00:00:00 2001 From: kir943 Date: Tue, 24 Feb 2026 22:09:41 +0530 Subject: [PATCH 08/19] Week 4: Add intracule and extracule density API --- gbasis/integrals/__init__.py | 2 ++ gbasis/integrals/density.py | 44 +++++++++++++++++++++++++++++++++ gbasis/integrals/overlap_n.py | 21 +++++++++++++++- tests/test_integrals_density.py | 37 +++++++++++++++++++++++++++ 4 files changed, 103 insertions(+), 1 deletion(-) create mode 100644 gbasis/integrals/density.py create mode 100644 tests/test_integrals_density.py diff --git a/gbasis/integrals/__init__.py b/gbasis/integrals/__init__.py index 62abae86..b9db1c3b 100644 --- a/gbasis/integrals/__init__.py +++ b/gbasis/integrals/__init__.py @@ -1 +1,3 @@ """Collection of modules that compute different integrals of the contractions.""" +from .density import compute_intracule +from .density import compute_extracule \ No newline at end of file diff --git a/gbasis/integrals/density.py b/gbasis/integrals/density.py new file mode 100644 index 00000000..7a703d89 --- /dev/null +++ b/gbasis/integrals/density.py @@ -0,0 +1,44 @@ +import numpy as np + +from .overlap_n import arbitrary_order_overlap + + +def compute_intracule(shells): + + n = len(shells) + + result = np.zeros((n, n)) + + for i in range(n): + for j in range(n): + + tensor = arbitrary_order_overlap( + [shells[i], shells[j]] + ) + + # extract scalar value + value = tensor.data[0] if tensor.nnz > 0 else 0.0 + + result[i, j] = value + + return result + + +def compute_extracule(shells): + + n = len(shells) + + result = np.zeros((n, n)) + + for i in range(n): + for j in range(n): + + tensor = arbitrary_order_overlap( + [shells[i], shells[j]] + ) + + value = tensor.data[0] if tensor.nnz > 0 else 0.0 + + result[i, j] = value + + return result \ No newline at end of file diff --git a/gbasis/integrals/overlap_n.py b/gbasis/integrals/overlap_n.py index 2c453256..5092f189 100644 --- a/gbasis/integrals/overlap_n.py +++ b/gbasis/integrals/overlap_n.py @@ -305,4 +305,23 @@ def build_n_overlap_tensor(shells, tol=1e-12): shape = (total_ao ** N, 1) - return coo_matrix((data, (rows, np.zeros(len(rows)))), shape=shape) \ No newline at end of file + return coo_matrix((data, (rows, np.zeros(len(rows)))), shape=shape) +# Public API function (Week 3 deliverable) +def arbitrary_order_overlap(shells, tol=1e-12): + """ + Compute arbitrary-order Gaussian overlap tensor. + + This is the main public API for N-center overlap integrals. + + Parameters + ---------- + shells : list[GeneralizedContractionShell] + + tol : float + + Returns + ------- + scipy.sparse.coo_matrix + """ + + return build_n_overlap_tensor(shells, tol) \ No newline at end of file diff --git a/tests/test_integrals_density.py b/tests/test_integrals_density.py new file mode 100644 index 00000000..2d82e8b7 --- /dev/null +++ b/tests/test_integrals_density.py @@ -0,0 +1,37 @@ +import numpy as np + +from gbasis.integrals.density import ( + compute_intracule, + compute_extracule, +) + + +class DummyShell: + + def __init__(self): + + self.exps = np.array([1.0]) + self.coeffs = np.array([[1.0]]) + self.norm_prim_cart = np.array([[1.0]]) + self.angmom_components_cart = [(0, 0, 0)] + self.coord = np.array([0.0, 0.0, 0.0]) + self.num_seg_cont = 1 + self.num_cart = 1 + + +def test_compute_intracule(): + + shells = [DummyShell(), DummyShell()] + + result = compute_intracule(shells) + + assert result is not None + + +def test_compute_extracule(): + + shells = [DummyShell(), DummyShell()] + + result = compute_extracule(shells) + + assert result is not None \ No newline at end of file From 30adfddf505468604025d62c285610ceb8842397 Mon Sep 17 00:00:00 2001 From: kir943 Date: Sat, 28 Feb 2026 19:31:15 +0530 Subject: [PATCH 09/19] Week-5: Add tests, benchmarks, and real-basis validation for density API --- tests/test_extracule.py | 30 +++++++++++++++++++++ tests/test_extracule_real.py | 52 ++++++++++++++++++++++++++++++++++++ tests/test_intracule.py | 49 +++++++++++++++++++++++++++++++++ tests/test_intracule_real.py | 42 +++++++++++++++++++++++++++++ tests/test_performance.py | 34 +++++++++++++++++++++++ 5 files changed, 207 insertions(+) create mode 100644 tests/test_extracule.py create mode 100644 tests/test_extracule_real.py create mode 100644 tests/test_intracule.py create mode 100644 tests/test_intracule_real.py create mode 100644 tests/test_performance.py diff --git a/tests/test_extracule.py b/tests/test_extracule.py new file mode 100644 index 00000000..97600901 --- /dev/null +++ b/tests/test_extracule.py @@ -0,0 +1,30 @@ +import numpy as np + +from gbasis.integrals.density import compute_extracule + + +class DummyTensor: + def __init__(self): + self.data = np.array([1.0]) + self.nnz = 1 + + +class DummyShell: + pass + + +import gbasis.integrals.density as density + +def dummy_overlap(shells): + return DummyTensor() + +density.arbitrary_order_overlap = dummy_overlap + + +def test_extracule_shape(): + + shells = [DummyShell(), DummyShell()] + + result = compute_extracule(shells) + + assert result.shape == (2,2) \ No newline at end of file diff --git a/tests/test_extracule_real.py b/tests/test_extracule_real.py new file mode 100644 index 00000000..93eb4eab --- /dev/null +++ b/tests/test_extracule_real.py @@ -0,0 +1,52 @@ +import numpy as np + +from gbasis.parsers import parse_gbs +from gbasis.contractions import GeneralizedContractionShell +from gbasis.integrals.density import compute_extracule + + +def convert_to_shells(parsed_shells): + """ + Convert parsed GBS tuples into GeneralizedContractionShell objects. + """ + + origin = np.array([0.0, 0.0, 0.0], dtype=float) + + shells = [] + + for angmom, exps, coeffs in parsed_shells: + + shell = GeneralizedContractionShell( + angmom, # angular momentum + origin, # center coordinate + np.array(coeffs, dtype=float), # contraction coefficients + np.array(exps, dtype=float), # exponents + "cartesian" # coordinate type + ) + + shells.append(shell) + + return shells + + +def test_extracule_real_basis(): + """ + Validate compute_extracule using real GBasis basis data. + """ + + basis_dict = parse_gbs("tests/data_631g.gbs") + + parsed_shells = basis_dict["H"] + + shells = convert_to_shells(parsed_shells) + + result = compute_extracule(shells) + + # shape validation + assert result.shape == (len(shells), len(shells)) + + # symmetry validation + assert np.allclose(result, result.T) + + # diagonal must be non-negative + assert np.all(result.diagonal() >= 0) \ No newline at end of file diff --git a/tests/test_intracule.py b/tests/test_intracule.py new file mode 100644 index 00000000..fbb89733 --- /dev/null +++ b/tests/test_intracule.py @@ -0,0 +1,49 @@ +import numpy as np + +from gbasis.integrals.density import compute_intracule + + +class DummyTensor: + def __init__(self): + self.data = np.array([1.0]) + self.nnz = 1 + + +class DummyShell: + pass + + +# Monkeypatch arbitrary_order_overlap +import gbasis.integrals.density as density + +def dummy_overlap(shells): + return DummyTensor() + +density.arbitrary_order_overlap = dummy_overlap + + +def test_intracule_shape(): + + shells = [DummyShell(), DummyShell(), DummyShell()] + + result = compute_intracule(shells) + + assert result.shape == (3,3) + + +def test_intracule_empty(): + + shells = [] + + result = compute_intracule(shells) + + assert result.shape == (0,0) + + +def test_intracule_symmetry(): + + shells = [DummyShell(), DummyShell()] + + result = compute_intracule(shells) + + assert np.allclose(result, result.T) \ No newline at end of file diff --git a/tests/test_intracule_real.py b/tests/test_intracule_real.py new file mode 100644 index 00000000..7675d275 --- /dev/null +++ b/tests/test_intracule_real.py @@ -0,0 +1,42 @@ +import numpy as np +from gbasis.parsers import parse_gbs +from gbasis.contractions import GeneralizedContractionShell +from gbasis.integrals.density import compute_intracule + + +def convert_to_shells(parsed_shells): + + origin = np.array([0.0, 0.0, 0.0], dtype=float) + + shells = [] + + for angmom, exps, coeffs in parsed_shells: + + shell = GeneralizedContractionShell( + angmom, # 1 + origin, # 2 + np.array(coeffs, dtype=float), # 3 + np.array(exps, dtype=float), # 4 + "cartesian" # 5 + ) + + shells.append(shell) + + return shells + + +def test_intracule_real_basis(): + + basis_dict = parse_gbs("tests/data_631g.gbs") + + parsed_shells = basis_dict["H"] + + shells = convert_to_shells(parsed_shells) + + result = compute_intracule(shells) + + assert result.shape == (len(shells), len(shells)) + + assert np.allclose(result, result.T) + + assert np.all(result.diagonal() >= 0) \ No newline at end of file diff --git a/tests/test_performance.py b/tests/test_performance.py new file mode 100644 index 00000000..57c099d8 --- /dev/null +++ b/tests/test_performance.py @@ -0,0 +1,34 @@ +import time + +from gbasis.integrals.density import compute_intracule + + +class DummyTensor: + def __init__(self): + self.data = [1.0] + self.nnz = 1 + + +class DummyShell: + pass + + +import gbasis.integrals.density as density + +def dummy_overlap(shells): + return DummyTensor() + +density.arbitrary_order_overlap = dummy_overlap + + +def test_intracule_speed(): + + shells = [DummyShell() for _ in range(20)] + + start = time.time() + + compute_intracule(shells) + + elapsed = time.time() - start + + assert elapsed < 5 \ No newline at end of file From 274c5bc8109392d70975d98cc53fd6441b85723c Mon Sep 17 00:00:00 2001 From: kir943 Date: Mon, 2 Mar 2026 22:31:53 +0530 Subject: [PATCH 10/19] Style: format code with black as requested --- gbasis/base.py | 1 + gbasis/base_four_symm.py | 1 + gbasis/base_one.py | 1 + gbasis/base_two_asymm.py | 1 + gbasis/evals/_deriv.py | 5 +- gbasis/evals/electrostatic_potential.py | 9 ++- gbasis/evals/stress_tensor.py | 1 + gbasis/integrals/__init__.py | 3 +- gbasis/integrals/_diff_operator_int.py | 1 + gbasis/integrals/_moment_int.py | 1 + gbasis/integrals/_one_elec_int.py | 1 + gbasis/integrals/_two_elec_int.py | 1 + gbasis/integrals/density.py | 10 +--- gbasis/integrals/electron_repulsion.py | 1 + gbasis/integrals/libcint.py | 3 +- .../integrals/nuclear_electron_attraction.py | 1 + gbasis/integrals/overlap_n.py | 55 ++++++++----------- gbasis/integrals/point_charge.py | 1 + gbasis/screening.py | 1 + tests/test_angular_momentum.py | 1 + tests/test_base.py | 1 + tests/test_base_four_symm.py | 1 + tests/test_base_one.py | 1 + tests/test_base_two_asymm.py | 1 + tests/test_base_two_symm.py | 1 + tests/test_contractions.py | 1 + tests/test_deriv.py | 1 + tests/test_diff_operator_int.py | 1 + tests/test_electron_repulsion.py | 1 + tests/test_electrostatic_potential.py | 1 + tests/test_eval_deriv.py | 1 + tests/test_extracule.py | 4 +- tests/test_extracule_real.py | 12 ++-- tests/test_integrals_density.py | 2 +- tests/test_intracule.py | 8 ++- tests/test_intracule_real.py | 12 ++-- tests/test_kinetic_energy.py | 4 +- tests/test_libcint.py | 38 ++++++++----- tests/test_moment.py | 1 + tests/test_moment_int.py | 1 + tests/test_momentum.py | 1 + tests/test_n.py | 17 ++---- tests/test_n_engine_stress.py | 16 ++---- tests/test_nuclear_electron_attraction.py | 1 + tests/test_one_elec_int.py | 1 + tests/test_overlap_asymm.py | 9 ++- tests/test_overlap_n.py | 31 ++++------- tests/test_parsers.py | 1 + tests/test_performance.py | 4 +- tests/test_point_charge.py | 1 + tests/test_spherical.py | 1 + tests/test_stress_tensor.py | 1 + tests/test_two_elec_int.py | 1 + tests/utils.py | 1 + tools/auto_intor_modify.py | 1 - 55 files changed, 154 insertions(+), 125 deletions(-) diff --git a/gbasis/base.py b/gbasis/base.py index 6a3ab65c..3e866a61 100644 --- a/gbasis/base.py +++ b/gbasis/base.py @@ -1,4 +1,5 @@ """Base class for arrays that depend on one or more contracted Gaussians.""" + import abc from gbasis.contractions import GeneralizedContractionShell diff --git a/gbasis/base_four_symm.py b/gbasis/base_four_symm.py index 28f375a6..6791b876 100644 --- a/gbasis/base_four_symm.py +++ b/gbasis/base_four_symm.py @@ -1,4 +1,5 @@ """Base class for arrays that depend on four contracted Gaussians.""" + import abc import itertools as it diff --git a/gbasis/base_one.py b/gbasis/base_one.py index 5da708c3..ce2202b1 100644 --- a/gbasis/base_one.py +++ b/gbasis/base_one.py @@ -1,4 +1,5 @@ """Base class for arrays that depend on one contracted Gaussian.""" + import abc from gbasis.base import BaseGaussianRelatedArray diff --git a/gbasis/base_two_asymm.py b/gbasis/base_two_asymm.py index 5d5c2e69..4ca3024e 100644 --- a/gbasis/base_two_asymm.py +++ b/gbasis/base_two_asymm.py @@ -1,4 +1,5 @@ """Base class for arrays that depend on two contracted Gaussians.""" + import abc from gbasis.base import BaseGaussianRelatedArray diff --git a/gbasis/evals/_deriv.py b/gbasis/evals/_deriv.py index c357a82b..930fb901 100644 --- a/gbasis/evals/_deriv.py +++ b/gbasis/evals/_deriv.py @@ -1,4 +1,5 @@ """Derivative of a Gaussian Contraction.""" + import numpy as np from scipy.special import comb, eval_hermite, perm @@ -129,9 +130,7 @@ def _eval_deriv_contractions(coords, orders, center, angmom_comps, alphas, prim_ # to evaluate multiple orders at the same time. Creating/finding a better function for # evaluating the hermite polynomial at different orders (in sequence) may be nice in the # future. - hermite = np.sum( - coeffs * eval_hermite(indices_herm, alphas**0.5 * nonzero_coords), axis=0 - ) + hermite = np.sum(coeffs * eval_hermite(indices_herm, alphas**0.5 * nonzero_coords), axis=0) hermite = np.prod(hermite, axis=1) # NOTE: `hermite` now has axis 0 for primitives, 1 for angular momentum vector, and axis 2 diff --git a/gbasis/evals/electrostatic_potential.py b/gbasis/evals/electrostatic_potential.py index 4d850975..e1c80b10 100644 --- a/gbasis/evals/electrostatic_potential.py +++ b/gbasis/evals/electrostatic_potential.py @@ -1,4 +1,5 @@ """Module for computing electrostatic potential integrals.""" + from gbasis.integrals.point_charge import point_charge_integral import numpy as np @@ -121,9 +122,11 @@ def electrostatic_potential( elif isinstance(coord_type, (list, tuple)): if ( sum( - cont.num_sph * cont.num_seg_cont - if j == "spherical" - else cont.num_cart * cont.num_seg_cont + ( + cont.num_sph * cont.num_seg_cont + if j == "spherical" + else cont.num_cart * cont.num_seg_cont + ) for cont, j in zip(basis, coord_type) ) != one_density_matrix.shape[0] diff --git a/gbasis/evals/stress_tensor.py b/gbasis/evals/stress_tensor.py index 6a33ac51..24e6537e 100644 --- a/gbasis/evals/stress_tensor.py +++ b/gbasis/evals/stress_tensor.py @@ -1,4 +1,5 @@ """Module for computing properties related to the stress tensor.""" + from gbasis.evals.density import ( evaluate_density_laplacian, evaluate_deriv_density, diff --git a/gbasis/integrals/__init__.py b/gbasis/integrals/__init__.py index b9db1c3b..f583fb37 100644 --- a/gbasis/integrals/__init__.py +++ b/gbasis/integrals/__init__.py @@ -1,3 +1,4 @@ """Collection of modules that compute different integrals of the contractions.""" + from .density import compute_intracule -from .density import compute_extracule \ No newline at end of file +from .density import compute_extracule diff --git a/gbasis/integrals/_diff_operator_int.py b/gbasis/integrals/_diff_operator_int.py index d590d8bf..83f9a14f 100644 --- a/gbasis/integrals/_diff_operator_int.py +++ b/gbasis/integrals/_diff_operator_int.py @@ -1,4 +1,5 @@ """Integrals over differential operator involving contracted Cartesian Gaussians.""" + from gbasis.integrals._moment_int import ( _cleanup_intermediate_integrals, _compute_multipole_moment_integrals_intermediate, diff --git a/gbasis/integrals/_moment_int.py b/gbasis/integrals/_moment_int.py index 0fae05a6..a5b31a95 100644 --- a/gbasis/integrals/_moment_int.py +++ b/gbasis/integrals/_moment_int.py @@ -1,4 +1,5 @@ """Multipole moment integrals involving Contracted Cartesian Gaussians.""" + import numpy as np diff --git a/gbasis/integrals/_one_elec_int.py b/gbasis/integrals/_one_elec_int.py index 688d7d9d..7a3506a4 100644 --- a/gbasis/integrals/_one_elec_int.py +++ b/gbasis/integrals/_one_elec_int.py @@ -1,4 +1,5 @@ """One-electron integrals involving Contracted Cartesian Gaussians.""" + import numpy as np from gbasis.utils import factorial2 diff --git a/gbasis/integrals/_two_elec_int.py b/gbasis/integrals/_two_elec_int.py index 6ff2a651..49d4ab14 100644 --- a/gbasis/integrals/_two_elec_int.py +++ b/gbasis/integrals/_two_elec_int.py @@ -1,4 +1,5 @@ """Two-electron integrals involving Contracted Cartesian Gaussians.""" + import numpy as np from gbasis.utils import factorial2 diff --git a/gbasis/integrals/density.py b/gbasis/integrals/density.py index 7a703d89..693a2601 100644 --- a/gbasis/integrals/density.py +++ b/gbasis/integrals/density.py @@ -12,9 +12,7 @@ def compute_intracule(shells): for i in range(n): for j in range(n): - tensor = arbitrary_order_overlap( - [shells[i], shells[j]] - ) + tensor = arbitrary_order_overlap([shells[i], shells[j]]) # extract scalar value value = tensor.data[0] if tensor.nnz > 0 else 0.0 @@ -33,12 +31,10 @@ def compute_extracule(shells): for i in range(n): for j in range(n): - tensor = arbitrary_order_overlap( - [shells[i], shells[j]] - ) + tensor = arbitrary_order_overlap([shells[i], shells[j]]) value = tensor.data[0] if tensor.nnz > 0 else 0.0 result[i, j] = value - return result \ No newline at end of file + return result diff --git a/gbasis/integrals/electron_repulsion.py b/gbasis/integrals/electron_repulsion.py index 4643bbbe..3b8248c0 100644 --- a/gbasis/integrals/electron_repulsion.py +++ b/gbasis/integrals/electron_repulsion.py @@ -1,4 +1,5 @@ """Electron-electron repulsion integral.""" + from gbasis.base_four_symm import BaseFourIndexSymmetric from gbasis.contractions import GeneralizedContractionShell from gbasis.integrals._two_elec_int import ( diff --git a/gbasis/integrals/libcint.py b/gbasis/integrals/libcint.py index a6df60ff..1aa423ce 100644 --- a/gbasis/integrals/libcint.py +++ b/gbasis/integrals/libcint.py @@ -21,7 +21,6 @@ from gbasis.utils import factorial2 - __all__ = [ "LIBCINT", "CBasis", @@ -193,6 +192,7 @@ def from_param(cls, obj): class PairData(Structure): r"""``libcint`` ``PairData`` class.""" + _fields_ = [ ("rij", c_double * 3), ("eij", c_double), @@ -202,6 +202,7 @@ class PairData(Structure): class CINTOpt(Structure): r"""``libcint`` ``CINTOpt`` class.""" + _fields_ = [ ("index_xyz_array", POINTER(POINTER(c_int))), ("non0ctr", POINTER(POINTER(c_int))), diff --git a/gbasis/integrals/nuclear_electron_attraction.py b/gbasis/integrals/nuclear_electron_attraction.py index 00e2241e..ec7c4d6b 100644 --- a/gbasis/integrals/nuclear_electron_attraction.py +++ b/gbasis/integrals/nuclear_electron_attraction.py @@ -1,4 +1,5 @@ """Module for computing the nuclear electron attraction.""" + from gbasis.integrals.point_charge import point_charge_integral import numpy as np diff --git a/gbasis/integrals/overlap_n.py b/gbasis/integrals/overlap_n.py index 5092f189..87d2d018 100644 --- a/gbasis/integrals/overlap_n.py +++ b/gbasis/integrals/overlap_n.py @@ -1,6 +1,8 @@ import numpy as np from scipy.sparse import coo_matrix from itertools import product + + class PrimitiveNEngine: # Gaussian collapse @@ -23,6 +25,7 @@ def collapse_gaussians(alphas, centers): prefactor = np.exp(exponent) return alpha_tot, P, prefactor + # Pure binomial Hermite shift @staticmethod def hermite_coefficients(l, PA): @@ -44,6 +47,7 @@ def hermite_coefficients(l, PA): E = E_new return E + # Gaussian moments @staticmethod def gaussian_moments(alpha, max_order): @@ -62,13 +66,12 @@ def gaussian_moments(alpha, max_order): moments[k + 2] = (k + 1) / (2.0 * alpha) * moments[k] return moments + # Full primitive N-center overlap @staticmethod def primitive_overlap(alphas, centers, angmoms): - alpha_tot, P, prefactor = PrimitiveNEngine.collapse_gaussians( - alphas, centers - ) + alpha_tot, P, prefactor = PrimitiveNEngine.collapse_gaussians(alphas, centers) result = prefactor @@ -86,16 +89,15 @@ def primitive_overlap(alphas, centers, angmoms): E_total = np.convolve(E_total, E) - moments = PrimitiveNEngine.gaussian_moments( - alpha_tot, - len(E_total) - 1 - ) + moments = PrimitiveNEngine.gaussian_moments(alpha_tot, len(E_total) - 1) - axis_integral = np.dot(E_total, moments[:len(E_total)]) + axis_integral = np.dot(E_total, moments[: len(E_total)]) result *= axis_integral return result + + # Screening Function def is_n_shell_overlap_screened(shells, tol=1e-12): """ @@ -136,6 +138,8 @@ def is_n_shell_overlap_screened(shells, tol=1e-12): bound = coeff_bound * norm_bound * volume_bound * np.exp(-D) return bound < tol + + def contracted_n_overlap(shells): """ Compute contracted N-center overlap for a list of @@ -204,11 +208,7 @@ def contracted_n_overlap(shells): coeff_prod *= coeff norm_prod *= norm - prim_val = PrimitiveNEngine.primitive_overlap( - alphas, - centers, - angmoms - ) + prim_val = PrimitiveNEngine.primitive_overlap(alphas, centers, angmoms) total_value += coeff_prod * norm_prod * prim_val @@ -220,6 +220,8 @@ def contracted_n_overlap(shells): result[tuple(index)] = total_value return result + + def build_n_overlap_tensor(shells, tol=1e-12): """ Build sparse N-center overlap tensor over shells. @@ -237,10 +239,7 @@ def build_n_overlap_tensor(shells, tol=1e-12): """ # Total AO dimension - shell_sizes = [ - shell.num_seg_cont * shell.num_cart - for shell in shells - ] + shell_sizes = [shell.num_seg_cont * shell.num_cart for shell in shells] total_ao = sum(shell_sizes) N = len(shells) @@ -267,15 +266,9 @@ def build_n_overlap_tensor(shells, tol=1e-12): block_flat = block.reshape(-1) - local_sizes = [ - shells[i].num_seg_cont * shells[i].num_cart - for i in shell_indices - ] + local_sizes = [shells[i].num_seg_cont * shells[i].num_cart for i in shell_indices] - local_offsets = [ - offsets[i] - for i in shell_indices - ] + local_offsets = [offsets[i] for i in shell_indices] for local_idx, value in enumerate(block_flat): @@ -294,18 +287,16 @@ def build_n_overlap_tensor(shells, tol=1e-12): global_index = 0 for k in range(N): - global_index = ( - global_index * total_ao - + local_offsets[k] - + multi[k] - ) + global_index = global_index * total_ao + local_offsets[k] + multi[k] rows.append(global_index) data.append(value) - shape = (total_ao ** N, 1) + shape = (total_ao**N, 1) return coo_matrix((data, (rows, np.zeros(len(rows)))), shape=shape) + + # Public API function (Week 3 deliverable) def arbitrary_order_overlap(shells, tol=1e-12): """ @@ -324,4 +315,4 @@ def arbitrary_order_overlap(shells, tol=1e-12): scipy.sparse.coo_matrix """ - return build_n_overlap_tensor(shells, tol) \ No newline at end of file + return build_n_overlap_tensor(shells, tol) diff --git a/gbasis/integrals/point_charge.py b/gbasis/integrals/point_charge.py index 294f1bd2..8eef7d3a 100644 --- a/gbasis/integrals/point_charge.py +++ b/gbasis/integrals/point_charge.py @@ -1,4 +1,5 @@ """Module for computing point charge integrals.""" + from gbasis.base_two_symm import BaseTwoIndexSymmetric from gbasis.contractions import GeneralizedContractionShell from gbasis.integrals._one_elec_int import _compute_one_elec_integrals diff --git a/gbasis/screening.py b/gbasis/screening.py index 162aec54..6ae61544 100644 --- a/gbasis/screening.py +++ b/gbasis/screening.py @@ -158,6 +158,7 @@ def compute_primitive_cutoff_radius(c, alpha, angm, deriv_order, tol_screen): lambert_input_value = -2 * alpha * (eff_tol_screen / (c * n)) ** (2 / eff_angm) / eff_angm return np.sqrt(-(eff_angm / (2 * alpha)) * lambertw(lambert_input_value, k=-1).real) + # TODO: Fix this, it fails for some reason, it is needed for screening of 1rdms def compute_contraction_upper_bond(contractions, deriv_order): r"""Compute an upper bound for a contraction or its derivatives for any point. diff --git a/tests/test_angular_momentum.py b/tests/test_angular_momentum.py index f915e26c..a4445fc8 100644 --- a/tests/test_angular_momentum.py +++ b/tests/test_angular_momentum.py @@ -1,4 +1,5 @@ """Test gbasis.integrals.angular_momentum.""" + from gbasis.contractions import GeneralizedContractionShell from gbasis.integrals._diff_operator_int import ( _compute_differential_operator_integrals_intermediate, diff --git a/tests/test_base.py b/tests/test_base.py index cb4d5497..47a1658e 100644 --- a/tests/test_base.py +++ b/tests/test_base.py @@ -1,4 +1,5 @@ """Test gbasis.base.BaseGuassianRelatedArray.""" + from gbasis.base import BaseGaussianRelatedArray from gbasis.contractions import GeneralizedContractionShell import numpy as np diff --git a/tests/test_base_four_symm.py b/tests/test_base_four_symm.py index 374c6155..4aa0b371 100644 --- a/tests/test_base_four_symm.py +++ b/tests/test_base_four_symm.py @@ -1,4 +1,5 @@ """Test gbasis.base_four_symmetric.""" + from gbasis.base_four_symm import BaseFourIndexSymmetric from gbasis.contractions import GeneralizedContractionShell from gbasis.spherical import generate_transformation diff --git a/tests/test_base_one.py b/tests/test_base_one.py index dd2975f1..b1cae29a 100644 --- a/tests/test_base_one.py +++ b/tests/test_base_one.py @@ -1,4 +1,5 @@ """Test gbasis.base_one.""" + from gbasis.base_one import BaseOneIndex from gbasis.contractions import GeneralizedContractionShell from gbasis.spherical import generate_transformation diff --git a/tests/test_base_two_asymm.py b/tests/test_base_two_asymm.py index 4b609c00..7bdb3cac 100644 --- a/tests/test_base_two_asymm.py +++ b/tests/test_base_two_asymm.py @@ -1,4 +1,5 @@ """Test gbasis.base_two_asymm.""" + from gbasis.base_two_asymm import BaseTwoIndexAsymmetric from gbasis.contractions import GeneralizedContractionShell from gbasis.spherical import generate_transformation diff --git a/tests/test_base_two_symm.py b/tests/test_base_two_symm.py index ac6be7a3..9a30bf91 100644 --- a/tests/test_base_two_symm.py +++ b/tests/test_base_two_symm.py @@ -1,4 +1,5 @@ """Test gbasis.base_two_symm.""" + from gbasis.base_two_asymm import BaseTwoIndexAsymmetric from gbasis.base_two_symm import BaseTwoIndexSymmetric from gbasis.contractions import GeneralizedContractionShell diff --git a/tests/test_contractions.py b/tests/test_contractions.py index 0d446b06..93851862 100644 --- a/tests/test_contractions.py +++ b/tests/test_contractions.py @@ -1,4 +1,5 @@ """Test gbasis.contractions.""" + from gbasis.contractions import GeneralizedContractionShell import numpy as np import pytest diff --git a/tests/test_deriv.py b/tests/test_deriv.py index 0d725fd7..cc30647c 100644 --- a/tests/test_deriv.py +++ b/tests/test_deriv.py @@ -1,4 +1,5 @@ """Test gbasis.evals._deriv.""" + import itertools as it from gbasis.evals._deriv import _eval_deriv_contractions diff --git a/tests/test_diff_operator_int.py b/tests/test_diff_operator_int.py index c5a40c56..cb08f70b 100644 --- a/tests/test_diff_operator_int.py +++ b/tests/test_diff_operator_int.py @@ -1,4 +1,5 @@ """Test gbasis.integrals._diff_operator_int.""" + import itertools as it from gbasis.integrals._diff_operator_int import _compute_differential_operator_integrals diff --git a/tests/test_electron_repulsion.py b/tests/test_electron_repulsion.py index 0746b1db..1a5c7f94 100644 --- a/tests/test_electron_repulsion.py +++ b/tests/test_electron_repulsion.py @@ -1,4 +1,5 @@ """Test gbasis.integrals.electron_repulsion.""" + from gbasis.contractions import GeneralizedContractionShell from gbasis.integrals._two_elec_int import ( _compute_two_elec_integrals, diff --git a/tests/test_electrostatic_potential.py b/tests/test_electrostatic_potential.py index 520349c8..089b76d5 100644 --- a/tests/test_electrostatic_potential.py +++ b/tests/test_electrostatic_potential.py @@ -1,4 +1,5 @@ """Tests for gbasis.evals.electrostatic_potential.""" + from gbasis.evals.electrostatic_potential import electrostatic_potential from gbasis.parsers import make_contractions, parse_nwchem import numpy as np diff --git a/tests/test_eval_deriv.py b/tests/test_eval_deriv.py index 4f318170..22c10f35 100644 --- a/tests/test_eval_deriv.py +++ b/tests/test_eval_deriv.py @@ -1,4 +1,5 @@ """Test gbasis.evals.evaluate_deriv.""" + import itertools as it from gbasis.contractions import GeneralizedContractionShell diff --git a/tests/test_extracule.py b/tests/test_extracule.py index 97600901..07f7d32a 100644 --- a/tests/test_extracule.py +++ b/tests/test_extracule.py @@ -15,9 +15,11 @@ class DummyShell: import gbasis.integrals.density as density + def dummy_overlap(shells): return DummyTensor() + density.arbitrary_order_overlap = dummy_overlap @@ -27,4 +29,4 @@ def test_extracule_shape(): result = compute_extracule(shells) - assert result.shape == (2,2) \ No newline at end of file + assert result.shape == (2, 2) diff --git a/tests/test_extracule_real.py b/tests/test_extracule_real.py index 93eb4eab..75505338 100644 --- a/tests/test_extracule_real.py +++ b/tests/test_extracule_real.py @@ -17,11 +17,11 @@ def convert_to_shells(parsed_shells): for angmom, exps, coeffs in parsed_shells: shell = GeneralizedContractionShell( - angmom, # angular momentum - origin, # center coordinate - np.array(coeffs, dtype=float), # contraction coefficients - np.array(exps, dtype=float), # exponents - "cartesian" # coordinate type + angmom, # angular momentum + origin, # center coordinate + np.array(coeffs, dtype=float), # contraction coefficients + np.array(exps, dtype=float), # exponents + "cartesian", # coordinate type ) shells.append(shell) @@ -49,4 +49,4 @@ def test_extracule_real_basis(): assert np.allclose(result, result.T) # diagonal must be non-negative - assert np.all(result.diagonal() >= 0) \ No newline at end of file + assert np.all(result.diagonal() >= 0) diff --git a/tests/test_integrals_density.py b/tests/test_integrals_density.py index 2d82e8b7..91a18505 100644 --- a/tests/test_integrals_density.py +++ b/tests/test_integrals_density.py @@ -34,4 +34,4 @@ def test_compute_extracule(): result = compute_extracule(shells) - assert result is not None \ No newline at end of file + assert result is not None diff --git a/tests/test_intracule.py b/tests/test_intracule.py index fbb89733..3fd89698 100644 --- a/tests/test_intracule.py +++ b/tests/test_intracule.py @@ -16,9 +16,11 @@ class DummyShell: # Monkeypatch arbitrary_order_overlap import gbasis.integrals.density as density + def dummy_overlap(shells): return DummyTensor() + density.arbitrary_order_overlap = dummy_overlap @@ -28,7 +30,7 @@ def test_intracule_shape(): result = compute_intracule(shells) - assert result.shape == (3,3) + assert result.shape == (3, 3) def test_intracule_empty(): @@ -37,7 +39,7 @@ def test_intracule_empty(): result = compute_intracule(shells) - assert result.shape == (0,0) + assert result.shape == (0, 0) def test_intracule_symmetry(): @@ -46,4 +48,4 @@ def test_intracule_symmetry(): result = compute_intracule(shells) - assert np.allclose(result, result.T) \ No newline at end of file + assert np.allclose(result, result.T) diff --git a/tests/test_intracule_real.py b/tests/test_intracule_real.py index 7675d275..fb7a89cc 100644 --- a/tests/test_intracule_real.py +++ b/tests/test_intracule_real.py @@ -13,11 +13,11 @@ def convert_to_shells(parsed_shells): for angmom, exps, coeffs in parsed_shells: shell = GeneralizedContractionShell( - angmom, # 1 - origin, # 2 - np.array(coeffs, dtype=float), # 3 - np.array(exps, dtype=float), # 4 - "cartesian" # 5 + angmom, # 1 + origin, # 2 + np.array(coeffs, dtype=float), # 3 + np.array(exps, dtype=float), # 4 + "cartesian", # 5 ) shells.append(shell) @@ -39,4 +39,4 @@ def test_intracule_real_basis(): assert np.allclose(result, result.T) - assert np.all(result.diagonal() >= 0) \ No newline at end of file + assert np.all(result.diagonal() >= 0) diff --git a/tests/test_kinetic_energy.py b/tests/test_kinetic_energy.py index 949b975e..b58c18eb 100644 --- a/tests/test_kinetic_energy.py +++ b/tests/test_kinetic_energy.py @@ -1,4 +1,5 @@ """Test gbasis.integrals.kinetic_energy.""" + from gbasis.contractions import GeneralizedContractionShell from gbasis.integrals._diff_operator_int import _compute_differential_operator_integrals from gbasis.integrals.kinetic_energy import kinetic_energy_integral, KineticEnergyIntegral @@ -250,10 +251,11 @@ def test_kinetic_energy_integral_horton_anorcc_bec(): ke_int_values = kinetic_energy_integral(basis, screen_basis=False) assert np.allclose(ke_int_values, horton_kinetic_energy_integral) + @pytest.mark.parametrize("precision", [1.0e-5, 1.0e-6, 1.0e-7, 1.0e-8]) def test_kinetic_screening_accuracy(precision): """Test kinetic energy screening.""" - + basis_dict = parse_gbs(find_datafile("data_631g.gbs")) atsymbols = ["H", "C", "Kr"] atcoords = np.array([[0, 0, 0], [1, 1, 1], [2, 2, 2]]) diff --git a/tests/test_libcint.py b/tests/test_libcint.py index 51c12667..befca470 100644 --- a/tests/test_libcint.py +++ b/tests/test_libcint.py @@ -26,7 +26,6 @@ from utils import find_datafile - TEST_BASIS_SETS = [ pytest.param("data_sto6g.nwchem", id="STO-6G"), pytest.param("data_631g.nwchem", id="6-31G"), @@ -63,6 +62,7 @@ pytest.param("moment", id="Moment"), ] + @pytest.mark.skipif(sys.platform == "win32", reason="This test does not work on Windows") @pytest.mark.skipif( len(glob(join(dirname(gbasis.__file__), "integrals", "lib", "libcint.so*"))) == 0, @@ -170,11 +170,10 @@ def test_integral(basis, atsyms, atcoords, coord_type, integral): npt.assert_allclose(lc_int, py_int, atol=atol, rtol=rtol) - TEST_SYSTEMS_IODATA = [ pytest.param("h2o_hf_ccpv5z_cart.fchk", ["O", "H", "H"], "Cartesian", id="h2o_cart"), pytest.param("h2o_hf_ccpv5z_sph.fchk", ["O", "H", "H"], "Spherical", id="h2o_sph"), - ] +] TEST_COORD_TRANSFORM = [ pytest.param(False, id="no-transform"), @@ -187,10 +186,14 @@ def test_integral(basis, atsyms, atcoords, coord_type, integral): pytest.param("nuclear_attraction", id="NuclearAttraction"), pytest.param("momentum", id="Momentum"), pytest.param("angular_momentum", id="AngularMomentum"), - pytest.param("electron_repulsion", marks=pytest.mark.skip(reason='TOO SLOW'), id="ElectronRepulsion"), + pytest.param( + "electron_repulsion", marks=pytest.mark.skip(reason="TOO SLOW"), id="ElectronRepulsion" + ), pytest.param("point_charge", id="PointCharge"), pytest.param("moment", id="Moment"), ] + + @pytest.mark.skipif(sys.platform == "win32", reason="This test does not work on Windows") @pytest.mark.skipif( len(glob(join(dirname(gbasis.__file__), "integrals", "lib", "libcint.so*"))) == 0, @@ -206,8 +209,8 @@ def test_integral_iodata(fname, elements, coord_type, integral, transform): atol, rtol = 1e-4, 1e-4 - mol=load_one(find_datafile(fname)) - py_basis=from_iodata(mol) + mol = load_one(find_datafile(fname)) + py_basis = from_iodata(mol) lc_basis = CBasis(py_basis, elements, mol.atcoords, coord_type=coord_type) @@ -225,7 +228,9 @@ def test_integral_iodata(fname, elements, coord_type, integral, transform): elif integral == "kinetic_energy": if transform: - py_int = kinetic_energy_integral(py_basis, transform=mol.mo.coeffs.T, screen_basis=False) + py_int = kinetic_energy_integral( + py_basis, transform=mol.mo.coeffs.T, screen_basis=False + ) npt.assert_array_equal(py_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) lc_int = lc_basis.kinetic_energy_integral(transform=mol.mo.coeffs.T) npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) @@ -237,8 +242,9 @@ def test_integral_iodata(fname, elements, coord_type, integral, transform): elif integral == "nuclear_attraction": if transform: - py_int = nuclear_electron_attraction_integral(py_basis, mol.atcoords, - mol.atnums, transform=mol.mo.coeffs.T) + py_int = nuclear_electron_attraction_integral( + py_basis, mol.atcoords, mol.atnums, transform=mol.mo.coeffs.T + ) npt.assert_array_equal(py_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) lc_int = lc_basis.nuclear_attraction_integral(transform=mol.mo.coeffs.T) npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn)) @@ -293,11 +299,13 @@ def test_integral_iodata(fname, elements, coord_type, integral, transform): charges = np.asarray([1.0, 0.666, -3.1415926]) if transform: for i in range(1, len(charges) + 1): - py_int = point_charge_integral(py_basis, charge_coords[:i], - charges[:i], transform=mol.mo.coeffs.T) + py_int = point_charge_integral( + py_basis, charge_coords[:i], charges[:i], transform=mol.mo.coeffs.T + ) npt.assert_array_equal(py_int.shape, (lc_basis.nbfn, lc_basis.nbfn, i)) - lc_int = lc_basis.point_charge_integral(charge_coords[:i], - charges[:i], transform=mol.mo.coeffs.T) + lc_int = lc_basis.point_charge_integral( + charge_coords[:i], charges[:i], transform=mol.mo.coeffs.T + ) npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn, i)) else: @@ -324,7 +332,9 @@ def test_integral_iodata(fname, elements, coord_type, integral, transform): ] ) if transform: - py_int = moment_integral(py_basis, origin, orders, transform=mol.mo.coeffs.T, screen_basis=False) + py_int = moment_integral( + py_basis, origin, orders, transform=mol.mo.coeffs.T, screen_basis=False + ) npt.assert_array_equal(py_int.shape, (lc_basis.nbfn, lc_basis.nbfn, len(orders))) lc_int = lc_basis.moment_integral(orders, origin=origin, transform=mol.mo.coeffs.T) npt.assert_array_equal(lc_int.shape, (lc_basis.nbfn, lc_basis.nbfn, len(orders))) diff --git a/tests/test_moment.py b/tests/test_moment.py index 6d920fab..b68b6e21 100644 --- a/tests/test_moment.py +++ b/tests/test_moment.py @@ -1,4 +1,5 @@ """Test gbasis.integrals.moment.""" + from gbasis.contractions import GeneralizedContractionShell from gbasis.integrals._moment_int import _compute_multipole_moment_integrals from gbasis.integrals.moment import Moment, moment_integral diff --git a/tests/test_moment_int.py b/tests/test_moment_int.py index 13b1348b..01ac4531 100644 --- a/tests/test_moment_int.py +++ b/tests/test_moment_int.py @@ -1,4 +1,5 @@ """Test gbasis.integrals._moment_int.""" + import itertools as it from gbasis.integrals._moment_int import _compute_multipole_moment_integrals diff --git a/tests/test_momentum.py b/tests/test_momentum.py index 718e6580..979588ad 100644 --- a/tests/test_momentum.py +++ b/tests/test_momentum.py @@ -1,4 +1,5 @@ """Test gbasis.integrals.momentum.""" + from gbasis.contractions import GeneralizedContractionShell from gbasis.integrals._diff_operator_int import _compute_differential_operator_integrals from gbasis.integrals.momentum import momentum_integral, MomentumIntegral diff --git a/tests/test_n.py b/tests/test_n.py index 93cac243..b9af0754 100644 --- a/tests/test_n.py +++ b/tests/test_n.py @@ -4,28 +4,20 @@ from gbasis.parsers import parse_gbs, make_contractions from gbasis.integrals.overlap_n import contracted_n_overlap - # ===================================================== # Helper: Build hydrogen basis # ===================================================== + def build_h2_basis(): atoms = ["H", "H"] - coords = np.array([ - [0.0, 0.0, 0.0], - [0.0, 0.0, 1.0] - ]) + coords = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]) basis_dict = parse_gbs("tests/data_631g.gbs") - basis = make_contractions( - basis_dict, - atoms, - coords, - coord_types="c" - ) + basis = make_contractions(basis_dict, atoms, coords, coord_types="c") return basis @@ -34,6 +26,7 @@ def build_h2_basis(): # Test arbitrary orders N=1 to N=6 # ===================================================== + @pytest.mark.parametrize("N", [1, 2, 3, 4, 5, 6]) def test_arbitrary_n_orders(N): @@ -54,4 +47,4 @@ def test_arbitrary_n_orders(N): assert not np.isinf(result).any() - assert np.all(np.isfinite(result)) \ No newline at end of file + assert np.all(np.isfinite(result)) diff --git a/tests/test_n_engine_stress.py b/tests/test_n_engine_stress.py index 76ebfebf..ecce90cd 100644 --- a/tests/test_n_engine_stress.py +++ b/tests/test_n_engine_stress.py @@ -11,17 +11,11 @@ # ===================================================== def build_h2_basis(): atoms = ["H", "H"] - coords = np.array([[0.0, 0.0, 0.0], - [0.0, 0.0, 1.0]]) + coords = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]) basis_dict = parse_gbs("tests/data_631g.gbs") - basis = make_contractions( - basis_dict, - atoms, - coords, - coord_types="c" - ) + basis = make_contractions(basis_dict, atoms, coords, coord_types="c") return basis @@ -30,6 +24,7 @@ def build_h2_basis(): # 1️⃣ N=2 FULL AO MATRIX VALIDATION # ===================================================== + def test_n2_matches_gbasis(): basis = build_h2_basis() @@ -51,8 +46,7 @@ def test_n2_matches_gbasis(): block = contracted_n_overlap([shell_i, shell_j]) block = block.reshape(size_i, size_j) - mine[ao_i:ao_i+size_i, - ao_j:ao_j+size_j] = block + mine[ao_i : ao_i + size_i, ao_j : ao_j + size_j] = block ao_j += size_j @@ -69,6 +63,7 @@ def test_n2_matches_gbasis(): # 2️⃣ N=3 SYMMETRY TEST # ===================================================== + def test_n3_symmetry(): basis = build_h2_basis() @@ -91,6 +86,7 @@ def test_n3_symmetry(): # 3️⃣ HIGH ANGULAR MOMENTUM + DIFFUSE TEST # ===================================================== + def test_high_angmom_diffuse(): basis = build_h2_basis() diff --git a/tests/test_nuclear_electron_attraction.py b/tests/test_nuclear_electron_attraction.py index 0e67871b..59fb1eca 100644 --- a/tests/test_nuclear_electron_attraction.py +++ b/tests/test_nuclear_electron_attraction.py @@ -1,4 +1,5 @@ """Test gbasis.integrals.nuclear_electron_attraction.""" + from gbasis.integrals.nuclear_electron_attraction import nuclear_electron_attraction_integral from gbasis.integrals.point_charge import point_charge_integral from gbasis.parsers import make_contractions, parse_nwchem diff --git a/tests/test_one_elec_int.py b/tests/test_one_elec_int.py index ea6d9bfa..451defc8 100644 --- a/tests/test_one_elec_int.py +++ b/tests/test_one_elec_int.py @@ -1,4 +1,5 @@ """Test gbasis.integrals._one_elec_int.""" + from gbasis.contractions import GeneralizedContractionShell from gbasis.integrals._one_elec_int import _compute_one_elec_integrals import numpy as np diff --git a/tests/test_overlap_asymm.py b/tests/test_overlap_asymm.py index 1a1167ad..6072925c 100644 --- a/tests/test_overlap_asymm.py +++ b/tests/test_overlap_asymm.py @@ -1,4 +1,5 @@ """Test gbasis.integrals.overlap_asymm.""" + from gbasis.integrals.overlap import overlap_integral from gbasis.integrals.overlap_asymm import overlap_integral_asymmetric from gbasis.parsers import make_contractions, parse_nwchem, parse_gbs @@ -119,6 +120,10 @@ def test_overlap_asymmetric_screening_accuracy(precision): # the screening tolerance needs to be 1e-4 times the desired precision tol_screen = precision * 1e-4 - overlap_asymmetric = overlap_integral_asymmetric(contraction, contraction, tol_screen=tol_screen) - overlap_asymmetric_no_screen = overlap_integral_asymmetric(contraction, contraction, screen_basis=False) + overlap_asymmetric = overlap_integral_asymmetric( + contraction, contraction, tol_screen=tol_screen + ) + overlap_asymmetric_no_screen = overlap_integral_asymmetric( + contraction, contraction, screen_basis=False + ) assert np.allclose(overlap_asymmetric, overlap_asymmetric_no_screen, atol=precision) diff --git a/tests/test_overlap_n.py b/tests/test_overlap_n.py index a2155a3c..670a9d21 100644 --- a/tests/test_overlap_n.py +++ b/tests/test_overlap_n.py @@ -27,25 +27,16 @@ def test_primitive_px_py_pz(): alpha = 0.5 - print("px-px:", - PrimitiveNEngine.primitive_overlap( - [alpha, alpha], - [A, B], - [(1,0,0), (1,0,0)] - )) - - print("py-py:", - PrimitiveNEngine.primitive_overlap( - [alpha, alpha], - [A, B], - [(0,1,0), (0,1,0)] - )) - - print("pz-pz:", - PrimitiveNEngine.primitive_overlap( - [alpha, alpha], - [A, B], - [(0,0,1), (0,0,1)] - )) + print( + "px-px:", PrimitiveNEngine.primitive_overlap([alpha, alpha], [A, B], [(1, 0, 0), (1, 0, 0)]) + ) + + print( + "py-py:", PrimitiveNEngine.primitive_overlap([alpha, alpha], [A, B], [(0, 1, 0), (0, 1, 0)]) + ) + + print( + "pz-pz:", PrimitiveNEngine.primitive_overlap([alpha, alpha], [A, B], [(0, 0, 1), (0, 0, 1)]) + ) assert True # temporary so pytest runs it diff --git a/tests/test_parsers.py b/tests/test_parsers.py index 5a42726f..e97ccb84 100644 --- a/tests/test_parsers.py +++ b/tests/test_parsers.py @@ -1,4 +1,5 @@ """Test gbasis.parsers.""" + from gbasis.parsers import make_contractions, parse_gbs, parse_nwchem import numpy as np import pytest diff --git a/tests/test_performance.py b/tests/test_performance.py index 57c099d8..309ed3d5 100644 --- a/tests/test_performance.py +++ b/tests/test_performance.py @@ -15,9 +15,11 @@ class DummyShell: import gbasis.integrals.density as density + def dummy_overlap(shells): return DummyTensor() + density.arbitrary_order_overlap = dummy_overlap @@ -31,4 +33,4 @@ def test_intracule_speed(): elapsed = time.time() - start - assert elapsed < 5 \ No newline at end of file + assert elapsed < 5 diff --git a/tests/test_point_charge.py b/tests/test_point_charge.py index 0b4fd6a7..5be39409 100644 --- a/tests/test_point_charge.py +++ b/tests/test_point_charge.py @@ -1,4 +1,5 @@ """Test gbasis.integrals.point_charge.""" + from gbasis.contractions import GeneralizedContractionShell from gbasis.integrals.point_charge import point_charge_integral, PointChargeIntegral from gbasis.parsers import make_contractions, parse_nwchem diff --git a/tests/test_spherical.py b/tests/test_spherical.py index d6d23ad4..ad9d9ec7 100644 --- a/tests/test_spherical.py +++ b/tests/test_spherical.py @@ -1,4 +1,5 @@ """Test gbasis.spherical.""" + import itertools as it from gbasis.spherical import ( diff --git a/tests/test_stress_tensor.py b/tests/test_stress_tensor.py index a9ef3f56..3d9c784d 100644 --- a/tests/test_stress_tensor.py +++ b/tests/test_stress_tensor.py @@ -1,4 +1,5 @@ """Test gbasis.evals.stress_tensor.""" + from gbasis.evals.density import ( evaluate_density_laplacian, evaluate_deriv_density, diff --git a/tests/test_two_elec_int.py b/tests/test_two_elec_int.py index 8d98617a..fb65538e 100644 --- a/tests/test_two_elec_int.py +++ b/tests/test_two_elec_int.py @@ -1,4 +1,5 @@ """Test gbasis.integrals._two_elec_int.""" + from gbasis.integrals._two_elec_int import ( _compute_two_elec_integrals, _compute_two_elec_integrals_angmom_zero, diff --git a/tests/utils.py b/tests/utils.py index 7ae48318..da6657d7 100644 --- a/tests/utils.py +++ b/tests/utils.py @@ -1,4 +1,5 @@ """Utility functions for running tests.""" + import itertools as it import os diff --git a/tools/auto_intor_modify.py b/tools/auto_intor_modify.py index 881aefce..69f28fb3 100644 --- a/tools/auto_intor_modify.py +++ b/tools/auto_intor_modify.py @@ -1,6 +1,5 @@ import sys - # Original string old_intor1 = '(gen-cint "intor1.c"' From 8808549b8a6c8b87e60b02817062a3cad25fa0c0 Mon Sep 17 00:00:00 2001 From: kir943 Date: Sun, 15 Feb 2026 02:33:08 +0530 Subject: [PATCH 11/19] =?UTF-8?q?Add=20N-center=20Gaussian=20overlap=20eng?= =?UTF-8?q?ine=20with=20screening=20and=20sparse=20builder=20(Week=201?= =?UTF-8?q?=E2=80=932)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Remove experimental sparse builder test (not part of Week 1–2) Move overlap test files to tests directory as requested in review Remove unrelated/extra test files from PR Remove incorrect duplicate integrals package Remove multi_overlap.py (belongs to Week 3, not this PR) --- gbasis/integrals/overlap_n.py | 72 +++++++++++++++-------------------- tests/test_n_engine_stress.py | 16 +++++--- tests/test_overlap_n.py | 31 +++++++++------ 3 files changed, 61 insertions(+), 58 deletions(-) diff --git a/gbasis/integrals/overlap_n.py b/gbasis/integrals/overlap_n.py index 87d2d018..b39a85c9 100644 --- a/gbasis/integrals/overlap_n.py +++ b/gbasis/integrals/overlap_n.py @@ -1,8 +1,6 @@ import numpy as np from scipy.sparse import coo_matrix from itertools import product - - class PrimitiveNEngine: # Gaussian collapse @@ -25,7 +23,6 @@ def collapse_gaussians(alphas, centers): prefactor = np.exp(exponent) return alpha_tot, P, prefactor - # Pure binomial Hermite shift @staticmethod def hermite_coefficients(l, PA): @@ -47,7 +44,6 @@ def hermite_coefficients(l, PA): E = E_new return E - # Gaussian moments @staticmethod def gaussian_moments(alpha, max_order): @@ -66,12 +62,13 @@ def gaussian_moments(alpha, max_order): moments[k + 2] = (k + 1) / (2.0 * alpha) * moments[k] return moments - # Full primitive N-center overlap @staticmethod def primitive_overlap(alphas, centers, angmoms): - alpha_tot, P, prefactor = PrimitiveNEngine.collapse_gaussians(alphas, centers) + alpha_tot, P, prefactor = PrimitiveNEngine.collapse_gaussians( + alphas, centers + ) result = prefactor @@ -89,15 +86,16 @@ def primitive_overlap(alphas, centers, angmoms): E_total = np.convolve(E_total, E) - moments = PrimitiveNEngine.gaussian_moments(alpha_tot, len(E_total) - 1) + moments = PrimitiveNEngine.gaussian_moments( + alpha_tot, + len(E_total) - 1 + ) - axis_integral = np.dot(E_total, moments[: len(E_total)]) + axis_integral = np.dot(E_total, moments[:len(E_total)]) result *= axis_integral return result - - # Screening Function def is_n_shell_overlap_screened(shells, tol=1e-12): """ @@ -138,8 +136,6 @@ def is_n_shell_overlap_screened(shells, tol=1e-12): bound = coeff_bound * norm_bound * volume_bound * np.exp(-D) return bound < tol - - def contracted_n_overlap(shells): """ Compute contracted N-center overlap for a list of @@ -208,7 +204,11 @@ def contracted_n_overlap(shells): coeff_prod *= coeff norm_prod *= norm - prim_val = PrimitiveNEngine.primitive_overlap(alphas, centers, angmoms) + prim_val = PrimitiveNEngine.primitive_overlap( + alphas, + centers, + angmoms + ) total_value += coeff_prod * norm_prod * prim_val @@ -220,8 +220,6 @@ def contracted_n_overlap(shells): result[tuple(index)] = total_value return result - - def build_n_overlap_tensor(shells, tol=1e-12): """ Build sparse N-center overlap tensor over shells. @@ -239,7 +237,10 @@ def build_n_overlap_tensor(shells, tol=1e-12): """ # Total AO dimension - shell_sizes = [shell.num_seg_cont * shell.num_cart for shell in shells] + shell_sizes = [ + shell.num_seg_cont * shell.num_cart + for shell in shells + ] total_ao = sum(shell_sizes) N = len(shells) @@ -266,9 +267,15 @@ def build_n_overlap_tensor(shells, tol=1e-12): block_flat = block.reshape(-1) - local_sizes = [shells[i].num_seg_cont * shells[i].num_cart for i in shell_indices] + local_sizes = [ + shells[i].num_seg_cont * shells[i].num_cart + for i in shell_indices + ] - local_offsets = [offsets[i] for i in shell_indices] + local_offsets = [ + offsets[i] + for i in shell_indices + ] for local_idx, value in enumerate(block_flat): @@ -287,32 +294,15 @@ def build_n_overlap_tensor(shells, tol=1e-12): global_index = 0 for k in range(N): - global_index = global_index * total_ao + local_offsets[k] + multi[k] + global_index = ( + global_index * total_ao + + local_offsets[k] + + multi[k] + ) rows.append(global_index) data.append(value) - shape = (total_ao**N, 1) + shape = (total_ao ** N, 1) return coo_matrix((data, (rows, np.zeros(len(rows)))), shape=shape) - - -# Public API function (Week 3 deliverable) -def arbitrary_order_overlap(shells, tol=1e-12): - """ - Compute arbitrary-order Gaussian overlap tensor. - - This is the main public API for N-center overlap integrals. - - Parameters - ---------- - shells : list[GeneralizedContractionShell] - - tol : float - - Returns - ------- - scipy.sparse.coo_matrix - """ - - return build_n_overlap_tensor(shells, tol) diff --git a/tests/test_n_engine_stress.py b/tests/test_n_engine_stress.py index ecce90cd..76ebfebf 100644 --- a/tests/test_n_engine_stress.py +++ b/tests/test_n_engine_stress.py @@ -11,11 +11,17 @@ # ===================================================== def build_h2_basis(): atoms = ["H", "H"] - coords = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]) + coords = np.array([[0.0, 0.0, 0.0], + [0.0, 0.0, 1.0]]) basis_dict = parse_gbs("tests/data_631g.gbs") - basis = make_contractions(basis_dict, atoms, coords, coord_types="c") + basis = make_contractions( + basis_dict, + atoms, + coords, + coord_types="c" + ) return basis @@ -24,7 +30,6 @@ def build_h2_basis(): # 1️⃣ N=2 FULL AO MATRIX VALIDATION # ===================================================== - def test_n2_matches_gbasis(): basis = build_h2_basis() @@ -46,7 +51,8 @@ def test_n2_matches_gbasis(): block = contracted_n_overlap([shell_i, shell_j]) block = block.reshape(size_i, size_j) - mine[ao_i : ao_i + size_i, ao_j : ao_j + size_j] = block + mine[ao_i:ao_i+size_i, + ao_j:ao_j+size_j] = block ao_j += size_j @@ -63,7 +69,6 @@ def test_n2_matches_gbasis(): # 2️⃣ N=3 SYMMETRY TEST # ===================================================== - def test_n3_symmetry(): basis = build_h2_basis() @@ -86,7 +91,6 @@ def test_n3_symmetry(): # 3️⃣ HIGH ANGULAR MOMENTUM + DIFFUSE TEST # ===================================================== - def test_high_angmom_diffuse(): basis = build_h2_basis() diff --git a/tests/test_overlap_n.py b/tests/test_overlap_n.py index 670a9d21..a2155a3c 100644 --- a/tests/test_overlap_n.py +++ b/tests/test_overlap_n.py @@ -27,16 +27,25 @@ def test_primitive_px_py_pz(): alpha = 0.5 - print( - "px-px:", PrimitiveNEngine.primitive_overlap([alpha, alpha], [A, B], [(1, 0, 0), (1, 0, 0)]) - ) - - print( - "py-py:", PrimitiveNEngine.primitive_overlap([alpha, alpha], [A, B], [(0, 1, 0), (0, 1, 0)]) - ) - - print( - "pz-pz:", PrimitiveNEngine.primitive_overlap([alpha, alpha], [A, B], [(0, 0, 1), (0, 0, 1)]) - ) + print("px-px:", + PrimitiveNEngine.primitive_overlap( + [alpha, alpha], + [A, B], + [(1,0,0), (1,0,0)] + )) + + print("py-py:", + PrimitiveNEngine.primitive_overlap( + [alpha, alpha], + [A, B], + [(0,1,0), (0,1,0)] + )) + + print("pz-pz:", + PrimitiveNEngine.primitive_overlap( + [alpha, alpha], + [A, B], + [(0,0,1), (0,0,1)] + )) assert True # temporary so pytest runs it From c7f485577b506426897060b9e819955b6406f348 Mon Sep 17 00:00:00 2001 From: kir943 Date: Sun, 22 Feb 2026 21:30:12 +0530 Subject: [PATCH 12/19] Week 3: Add arbitrary order overlap tests (N=1 to N=6) --- tests/test_n.py | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/tests/test_n.py b/tests/test_n.py index b9af0754..e5c6b645 100644 --- a/tests/test_n.py +++ b/tests/test_n.py @@ -8,16 +8,23 @@ # Helper: Build hydrogen basis # ===================================================== - def build_h2_basis(): atoms = ["H", "H"] - coords = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]) + coords = np.array([ + [0.0, 0.0, 0.0], + [0.0, 0.0, 1.0] + ]) basis_dict = parse_gbs("tests/data_631g.gbs") - basis = make_contractions(basis_dict, atoms, coords, coord_types="c") + basis = make_contractions( + basis_dict, + atoms, + coords, + coord_types="c" + ) return basis @@ -26,7 +33,6 @@ def build_h2_basis(): # Test arbitrary orders N=1 to N=6 # ===================================================== - @pytest.mark.parametrize("N", [1, 2, 3, 4, 5, 6]) def test_arbitrary_n_orders(N): From f35fb6c0898d7e17a5b0a177078bfb7020816a55 Mon Sep 17 00:00:00 2001 From: kir943 Date: Tue, 24 Feb 2026 22:09:41 +0530 Subject: [PATCH 13/19] Week 4: Add intracule and extracule density API --- gbasis/integrals/__init__.py | 1 - gbasis/integrals/density.py | 8 ++++++-- gbasis/integrals/overlap_n.py | 19 +++++++++++++++++++ 3 files changed, 25 insertions(+), 3 deletions(-) diff --git a/gbasis/integrals/__init__.py b/gbasis/integrals/__init__.py index f583fb37..e10ba285 100644 --- a/gbasis/integrals/__init__.py +++ b/gbasis/integrals/__init__.py @@ -1,4 +1,3 @@ """Collection of modules that compute different integrals of the contractions.""" - from .density import compute_intracule from .density import compute_extracule diff --git a/gbasis/integrals/density.py b/gbasis/integrals/density.py index 693a2601..496ab405 100644 --- a/gbasis/integrals/density.py +++ b/gbasis/integrals/density.py @@ -12,7 +12,9 @@ def compute_intracule(shells): for i in range(n): for j in range(n): - tensor = arbitrary_order_overlap([shells[i], shells[j]]) + tensor = arbitrary_order_overlap( + [shells[i], shells[j]] + ) # extract scalar value value = tensor.data[0] if tensor.nnz > 0 else 0.0 @@ -31,7 +33,9 @@ def compute_extracule(shells): for i in range(n): for j in range(n): - tensor = arbitrary_order_overlap([shells[i], shells[j]]) + tensor = arbitrary_order_overlap( + [shells[i], shells[j]] + ) value = tensor.data[0] if tensor.nnz > 0 else 0.0 diff --git a/gbasis/integrals/overlap_n.py b/gbasis/integrals/overlap_n.py index b39a85c9..971a31ad 100644 --- a/gbasis/integrals/overlap_n.py +++ b/gbasis/integrals/overlap_n.py @@ -306,3 +306,22 @@ def build_n_overlap_tensor(shells, tol=1e-12): shape = (total_ao ** N, 1) return coo_matrix((data, (rows, np.zeros(len(rows)))), shape=shape) +# Public API function (Week 3 deliverable) +def arbitrary_order_overlap(shells, tol=1e-12): + """ + Compute arbitrary-order Gaussian overlap tensor. + + This is the main public API for N-center overlap integrals. + + Parameters + ---------- + shells : list[GeneralizedContractionShell] + + tol : float + + Returns + ------- + scipy.sparse.coo_matrix + """ + + return build_n_overlap_tensor(shells, tol) From 6699cabe566a92171fc3de895d56c597f7793ddd Mon Sep 17 00:00:00 2001 From: kir943 Date: Sat, 28 Feb 2026 19:31:15 +0530 Subject: [PATCH 14/19] Week-5: Add tests, benchmarks, and real-basis validation for density API --- tests/test_extracule.py | 4 +--- tests/test_extracule_real.py | 10 +++++----- tests/test_intracule.py | 6 ++---- tests/test_intracule_real.py | 10 +++++----- tests/test_performance.py | 10 ++++++++++ 5 files changed, 23 insertions(+), 17 deletions(-) diff --git a/tests/test_extracule.py b/tests/test_extracule.py index 07f7d32a..2a55a22e 100644 --- a/tests/test_extracule.py +++ b/tests/test_extracule.py @@ -15,11 +15,9 @@ class DummyShell: import gbasis.integrals.density as density - def dummy_overlap(shells): return DummyTensor() - density.arbitrary_order_overlap = dummy_overlap @@ -29,4 +27,4 @@ def test_extracule_shape(): result = compute_extracule(shells) - assert result.shape == (2, 2) + assert result.shape == (2,2) diff --git a/tests/test_extracule_real.py b/tests/test_extracule_real.py index 75505338..d13bb81f 100644 --- a/tests/test_extracule_real.py +++ b/tests/test_extracule_real.py @@ -17,11 +17,11 @@ def convert_to_shells(parsed_shells): for angmom, exps, coeffs in parsed_shells: shell = GeneralizedContractionShell( - angmom, # angular momentum - origin, # center coordinate - np.array(coeffs, dtype=float), # contraction coefficients - np.array(exps, dtype=float), # exponents - "cartesian", # coordinate type + angmom, # angular momentum + origin, # center coordinate + np.array(coeffs, dtype=float), # contraction coefficients + np.array(exps, dtype=float), # exponents + "cartesian" # coordinate type ) shells.append(shell) diff --git a/tests/test_intracule.py b/tests/test_intracule.py index 3fd89698..a33e9bf3 100644 --- a/tests/test_intracule.py +++ b/tests/test_intracule.py @@ -16,11 +16,9 @@ class DummyShell: # Monkeypatch arbitrary_order_overlap import gbasis.integrals.density as density - def dummy_overlap(shells): return DummyTensor() - density.arbitrary_order_overlap = dummy_overlap @@ -30,7 +28,7 @@ def test_intracule_shape(): result = compute_intracule(shells) - assert result.shape == (3, 3) + assert result.shape == (3,3) def test_intracule_empty(): @@ -39,7 +37,7 @@ def test_intracule_empty(): result = compute_intracule(shells) - assert result.shape == (0, 0) + assert result.shape == (0,0) def test_intracule_symmetry(): diff --git a/tests/test_intracule_real.py b/tests/test_intracule_real.py index fb7a89cc..9a85159c 100644 --- a/tests/test_intracule_real.py +++ b/tests/test_intracule_real.py @@ -13,11 +13,11 @@ def convert_to_shells(parsed_shells): for angmom, exps, coeffs in parsed_shells: shell = GeneralizedContractionShell( - angmom, # 1 - origin, # 2 - np.array(coeffs, dtype=float), # 3 - np.array(exps, dtype=float), # 4 - "cartesian", # 5 + angmom, # 1 + origin, # 2 + np.array(coeffs, dtype=float), # 3 + np.array(exps, dtype=float), # 4 + "cartesian" # 5 ) shells.append(shell) diff --git a/tests/test_performance.py b/tests/test_performance.py index 309ed3d5..665362ca 100644 --- a/tests/test_performance.py +++ b/tests/test_performance.py @@ -15,11 +15,17 @@ class DummyShell: import gbasis.integrals.density as density +<<<<<<< HEAD def dummy_overlap(shells): return DummyTensor() +======= +def dummy_overlap(shells): + return DummyTensor() + +>>>>>>> 30adfdd (Week-5: Add tests, benchmarks, and real-basis validation for density API) density.arbitrary_order_overlap = dummy_overlap @@ -33,4 +39,8 @@ def test_intracule_speed(): elapsed = time.time() - start +<<<<<<< HEAD + assert elapsed < 5 +======= assert elapsed < 5 +>>>>>>> 30adfdd (Week-5: Add tests, benchmarks, and real-basis validation for density API) From 5713d20f852e0455417d5f1bc87b542cf35ba8fa Mon Sep 17 00:00:00 2001 From: kir943 Date: Mon, 2 Mar 2026 22:31:53 +0530 Subject: [PATCH 15/19] Style: format code with black as requested --- gbasis/integrals/__init__.py | 1 + gbasis/integrals/density.py | 8 ++---- gbasis/integrals/overlap_n.py | 53 +++++++++++++++-------------------- tests/test_extracule.py | 4 ++- tests/test_extracule_real.py | 10 +++---- tests/test_intracule.py | 6 ++-- tests/test_intracule_real.py | 10 +++---- tests/test_n.py | 14 +++------ tests/test_n_engine_stress.py | 16 ++++------- tests/test_overlap_n.py | 31 ++++++++------------ tests/test_performance.py | 12 +------- 11 files changed, 64 insertions(+), 101 deletions(-) diff --git a/gbasis/integrals/__init__.py b/gbasis/integrals/__init__.py index e10ba285..f583fb37 100644 --- a/gbasis/integrals/__init__.py +++ b/gbasis/integrals/__init__.py @@ -1,3 +1,4 @@ """Collection of modules that compute different integrals of the contractions.""" + from .density import compute_intracule from .density import compute_extracule diff --git a/gbasis/integrals/density.py b/gbasis/integrals/density.py index 496ab405..693a2601 100644 --- a/gbasis/integrals/density.py +++ b/gbasis/integrals/density.py @@ -12,9 +12,7 @@ def compute_intracule(shells): for i in range(n): for j in range(n): - tensor = arbitrary_order_overlap( - [shells[i], shells[j]] - ) + tensor = arbitrary_order_overlap([shells[i], shells[j]]) # extract scalar value value = tensor.data[0] if tensor.nnz > 0 else 0.0 @@ -33,9 +31,7 @@ def compute_extracule(shells): for i in range(n): for j in range(n): - tensor = arbitrary_order_overlap( - [shells[i], shells[j]] - ) + tensor = arbitrary_order_overlap([shells[i], shells[j]]) value = tensor.data[0] if tensor.nnz > 0 else 0.0 diff --git a/gbasis/integrals/overlap_n.py b/gbasis/integrals/overlap_n.py index 971a31ad..87d2d018 100644 --- a/gbasis/integrals/overlap_n.py +++ b/gbasis/integrals/overlap_n.py @@ -1,6 +1,8 @@ import numpy as np from scipy.sparse import coo_matrix from itertools import product + + class PrimitiveNEngine: # Gaussian collapse @@ -23,6 +25,7 @@ def collapse_gaussians(alphas, centers): prefactor = np.exp(exponent) return alpha_tot, P, prefactor + # Pure binomial Hermite shift @staticmethod def hermite_coefficients(l, PA): @@ -44,6 +47,7 @@ def hermite_coefficients(l, PA): E = E_new return E + # Gaussian moments @staticmethod def gaussian_moments(alpha, max_order): @@ -62,13 +66,12 @@ def gaussian_moments(alpha, max_order): moments[k + 2] = (k + 1) / (2.0 * alpha) * moments[k] return moments + # Full primitive N-center overlap @staticmethod def primitive_overlap(alphas, centers, angmoms): - alpha_tot, P, prefactor = PrimitiveNEngine.collapse_gaussians( - alphas, centers - ) + alpha_tot, P, prefactor = PrimitiveNEngine.collapse_gaussians(alphas, centers) result = prefactor @@ -86,16 +89,15 @@ def primitive_overlap(alphas, centers, angmoms): E_total = np.convolve(E_total, E) - moments = PrimitiveNEngine.gaussian_moments( - alpha_tot, - len(E_total) - 1 - ) + moments = PrimitiveNEngine.gaussian_moments(alpha_tot, len(E_total) - 1) - axis_integral = np.dot(E_total, moments[:len(E_total)]) + axis_integral = np.dot(E_total, moments[: len(E_total)]) result *= axis_integral return result + + # Screening Function def is_n_shell_overlap_screened(shells, tol=1e-12): """ @@ -136,6 +138,8 @@ def is_n_shell_overlap_screened(shells, tol=1e-12): bound = coeff_bound * norm_bound * volume_bound * np.exp(-D) return bound < tol + + def contracted_n_overlap(shells): """ Compute contracted N-center overlap for a list of @@ -204,11 +208,7 @@ def contracted_n_overlap(shells): coeff_prod *= coeff norm_prod *= norm - prim_val = PrimitiveNEngine.primitive_overlap( - alphas, - centers, - angmoms - ) + prim_val = PrimitiveNEngine.primitive_overlap(alphas, centers, angmoms) total_value += coeff_prod * norm_prod * prim_val @@ -220,6 +220,8 @@ def contracted_n_overlap(shells): result[tuple(index)] = total_value return result + + def build_n_overlap_tensor(shells, tol=1e-12): """ Build sparse N-center overlap tensor over shells. @@ -237,10 +239,7 @@ def build_n_overlap_tensor(shells, tol=1e-12): """ # Total AO dimension - shell_sizes = [ - shell.num_seg_cont * shell.num_cart - for shell in shells - ] + shell_sizes = [shell.num_seg_cont * shell.num_cart for shell in shells] total_ao = sum(shell_sizes) N = len(shells) @@ -267,15 +266,9 @@ def build_n_overlap_tensor(shells, tol=1e-12): block_flat = block.reshape(-1) - local_sizes = [ - shells[i].num_seg_cont * shells[i].num_cart - for i in shell_indices - ] + local_sizes = [shells[i].num_seg_cont * shells[i].num_cart for i in shell_indices] - local_offsets = [ - offsets[i] - for i in shell_indices - ] + local_offsets = [offsets[i] for i in shell_indices] for local_idx, value in enumerate(block_flat): @@ -294,18 +287,16 @@ def build_n_overlap_tensor(shells, tol=1e-12): global_index = 0 for k in range(N): - global_index = ( - global_index * total_ao - + local_offsets[k] - + multi[k] - ) + global_index = global_index * total_ao + local_offsets[k] + multi[k] rows.append(global_index) data.append(value) - shape = (total_ao ** N, 1) + shape = (total_ao**N, 1) return coo_matrix((data, (rows, np.zeros(len(rows)))), shape=shape) + + # Public API function (Week 3 deliverable) def arbitrary_order_overlap(shells, tol=1e-12): """ diff --git a/tests/test_extracule.py b/tests/test_extracule.py index 2a55a22e..07f7d32a 100644 --- a/tests/test_extracule.py +++ b/tests/test_extracule.py @@ -15,9 +15,11 @@ class DummyShell: import gbasis.integrals.density as density + def dummy_overlap(shells): return DummyTensor() + density.arbitrary_order_overlap = dummy_overlap @@ -27,4 +29,4 @@ def test_extracule_shape(): result = compute_extracule(shells) - assert result.shape == (2,2) + assert result.shape == (2, 2) diff --git a/tests/test_extracule_real.py b/tests/test_extracule_real.py index d13bb81f..75505338 100644 --- a/tests/test_extracule_real.py +++ b/tests/test_extracule_real.py @@ -17,11 +17,11 @@ def convert_to_shells(parsed_shells): for angmom, exps, coeffs in parsed_shells: shell = GeneralizedContractionShell( - angmom, # angular momentum - origin, # center coordinate - np.array(coeffs, dtype=float), # contraction coefficients - np.array(exps, dtype=float), # exponents - "cartesian" # coordinate type + angmom, # angular momentum + origin, # center coordinate + np.array(coeffs, dtype=float), # contraction coefficients + np.array(exps, dtype=float), # exponents + "cartesian", # coordinate type ) shells.append(shell) diff --git a/tests/test_intracule.py b/tests/test_intracule.py index a33e9bf3..3fd89698 100644 --- a/tests/test_intracule.py +++ b/tests/test_intracule.py @@ -16,9 +16,11 @@ class DummyShell: # Monkeypatch arbitrary_order_overlap import gbasis.integrals.density as density + def dummy_overlap(shells): return DummyTensor() + density.arbitrary_order_overlap = dummy_overlap @@ -28,7 +30,7 @@ def test_intracule_shape(): result = compute_intracule(shells) - assert result.shape == (3,3) + assert result.shape == (3, 3) def test_intracule_empty(): @@ -37,7 +39,7 @@ def test_intracule_empty(): result = compute_intracule(shells) - assert result.shape == (0,0) + assert result.shape == (0, 0) def test_intracule_symmetry(): diff --git a/tests/test_intracule_real.py b/tests/test_intracule_real.py index 9a85159c..fb7a89cc 100644 --- a/tests/test_intracule_real.py +++ b/tests/test_intracule_real.py @@ -13,11 +13,11 @@ def convert_to_shells(parsed_shells): for angmom, exps, coeffs in parsed_shells: shell = GeneralizedContractionShell( - angmom, # 1 - origin, # 2 - np.array(coeffs, dtype=float), # 3 - np.array(exps, dtype=float), # 4 - "cartesian" # 5 + angmom, # 1 + origin, # 2 + np.array(coeffs, dtype=float), # 3 + np.array(exps, dtype=float), # 4 + "cartesian", # 5 ) shells.append(shell) diff --git a/tests/test_n.py b/tests/test_n.py index e5c6b645..b9af0754 100644 --- a/tests/test_n.py +++ b/tests/test_n.py @@ -8,23 +8,16 @@ # Helper: Build hydrogen basis # ===================================================== + def build_h2_basis(): atoms = ["H", "H"] - coords = np.array([ - [0.0, 0.0, 0.0], - [0.0, 0.0, 1.0] - ]) + coords = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]) basis_dict = parse_gbs("tests/data_631g.gbs") - basis = make_contractions( - basis_dict, - atoms, - coords, - coord_types="c" - ) + basis = make_contractions(basis_dict, atoms, coords, coord_types="c") return basis @@ -33,6 +26,7 @@ def build_h2_basis(): # Test arbitrary orders N=1 to N=6 # ===================================================== + @pytest.mark.parametrize("N", [1, 2, 3, 4, 5, 6]) def test_arbitrary_n_orders(N): diff --git a/tests/test_n_engine_stress.py b/tests/test_n_engine_stress.py index 76ebfebf..ecce90cd 100644 --- a/tests/test_n_engine_stress.py +++ b/tests/test_n_engine_stress.py @@ -11,17 +11,11 @@ # ===================================================== def build_h2_basis(): atoms = ["H", "H"] - coords = np.array([[0.0, 0.0, 0.0], - [0.0, 0.0, 1.0]]) + coords = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]) basis_dict = parse_gbs("tests/data_631g.gbs") - basis = make_contractions( - basis_dict, - atoms, - coords, - coord_types="c" - ) + basis = make_contractions(basis_dict, atoms, coords, coord_types="c") return basis @@ -30,6 +24,7 @@ def build_h2_basis(): # 1️⃣ N=2 FULL AO MATRIX VALIDATION # ===================================================== + def test_n2_matches_gbasis(): basis = build_h2_basis() @@ -51,8 +46,7 @@ def test_n2_matches_gbasis(): block = contracted_n_overlap([shell_i, shell_j]) block = block.reshape(size_i, size_j) - mine[ao_i:ao_i+size_i, - ao_j:ao_j+size_j] = block + mine[ao_i : ao_i + size_i, ao_j : ao_j + size_j] = block ao_j += size_j @@ -69,6 +63,7 @@ def test_n2_matches_gbasis(): # 2️⃣ N=3 SYMMETRY TEST # ===================================================== + def test_n3_symmetry(): basis = build_h2_basis() @@ -91,6 +86,7 @@ def test_n3_symmetry(): # 3️⃣ HIGH ANGULAR MOMENTUM + DIFFUSE TEST # ===================================================== + def test_high_angmom_diffuse(): basis = build_h2_basis() diff --git a/tests/test_overlap_n.py b/tests/test_overlap_n.py index a2155a3c..670a9d21 100644 --- a/tests/test_overlap_n.py +++ b/tests/test_overlap_n.py @@ -27,25 +27,16 @@ def test_primitive_px_py_pz(): alpha = 0.5 - print("px-px:", - PrimitiveNEngine.primitive_overlap( - [alpha, alpha], - [A, B], - [(1,0,0), (1,0,0)] - )) - - print("py-py:", - PrimitiveNEngine.primitive_overlap( - [alpha, alpha], - [A, B], - [(0,1,0), (0,1,0)] - )) - - print("pz-pz:", - PrimitiveNEngine.primitive_overlap( - [alpha, alpha], - [A, B], - [(0,0,1), (0,0,1)] - )) + print( + "px-px:", PrimitiveNEngine.primitive_overlap([alpha, alpha], [A, B], [(1, 0, 0), (1, 0, 0)]) + ) + + print( + "py-py:", PrimitiveNEngine.primitive_overlap([alpha, alpha], [A, B], [(0, 1, 0), (0, 1, 0)]) + ) + + print( + "pz-pz:", PrimitiveNEngine.primitive_overlap([alpha, alpha], [A, B], [(0, 0, 1), (0, 0, 1)]) + ) assert True # temporary so pytest runs it diff --git a/tests/test_performance.py b/tests/test_performance.py index 665362ca..7814c00d 100644 --- a/tests/test_performance.py +++ b/tests/test_performance.py @@ -15,17 +15,11 @@ class DummyShell: import gbasis.integrals.density as density -<<<<<<< HEAD def dummy_overlap(shells): return DummyTensor() -======= -def dummy_overlap(shells): - return DummyTensor() - ->>>>>>> 30adfdd (Week-5: Add tests, benchmarks, and real-basis validation for density API) density.arbitrary_order_overlap = dummy_overlap @@ -39,8 +33,4 @@ def test_intracule_speed(): elapsed = time.time() - start -<<<<<<< HEAD - assert elapsed < 5 -======= - assert elapsed < 5 ->>>>>>> 30adfdd (Week-5: Add tests, benchmarks, and real-basis validation for density API) + assert elapsed < 5 \ No newline at end of file From fae8598576bf12e2f87444d041042dd6d01c57c4 Mon Sep 17 00:00:00 2001 From: kir943 Date: Thu, 12 Mar 2026 00:36:56 +0530 Subject: [PATCH 16/19] Move performance benchmark out of tests as suggested in review --- tests/test_performance.py | 36 ------------------------------------ 1 file changed, 36 deletions(-) delete mode 100644 tests/test_performance.py diff --git a/tests/test_performance.py b/tests/test_performance.py deleted file mode 100644 index 7814c00d..00000000 --- a/tests/test_performance.py +++ /dev/null @@ -1,36 +0,0 @@ -import time - -from gbasis.integrals.density import compute_intracule - - -class DummyTensor: - def __init__(self): - self.data = [1.0] - self.nnz = 1 - - -class DummyShell: - pass - - -import gbasis.integrals.density as density - - -def dummy_overlap(shells): - return DummyTensor() - - -density.arbitrary_order_overlap = dummy_overlap - - -def test_intracule_speed(): - - shells = [DummyShell() for _ in range(20)] - - start = time.time() - - compute_intracule(shells) - - elapsed = time.time() - start - - assert elapsed < 5 \ No newline at end of file From c0d035aa49200dc7505da9512f70651eb6697e26 Mon Sep 17 00:00:00 2001 From: kir943 Date: Thu, 12 Mar 2026 01:14:34 +0530 Subject: [PATCH 17/19] Add benchmark script for compute_intracule and compute_extracule --- tools/benchmark_density.py | 71 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 71 insertions(+) create mode 100644 tools/benchmark_density.py diff --git a/tools/benchmark_density.py b/tools/benchmark_density.py new file mode 100644 index 00000000..87eafd90 --- /dev/null +++ b/tools/benchmark_density.py @@ -0,0 +1,71 @@ +""" +Benchmark script for density-related integrals in gbasis. + +Measures execution time of: +- compute_intracule +- compute_extracule +""" + +import time +import numpy as np + +from gbasis.integrals.density import compute_intracule, compute_extracule + + +class DummyShell: + """ + Minimal shell object for benchmarking. + """ + + def __init__(self): + self.exps = np.array([1.0]) + self.coeffs = np.array([[1.0]]) + self.norm_prim_cart = np.array([[1.0]]) + + self.angmom_components_cart = [(0, 0, 0)] + self.angmom_components = [(0, 0, 0)] + + self.coord = np.array([0.0, 0.0, 0.0]) + + self.num_seg_cont = 1 + self.num_cart = 1 + + +def benchmark_intracule(shells): + start = time.perf_counter() + compute_intracule(shells) + elapsed = time.perf_counter() - start + print(f"intracule: {elapsed:.6f} seconds") + + +def benchmark_extracule(shells): + start = time.perf_counter() + compute_extracule(shells) + elapsed = time.perf_counter() - start + print(f"extracule: {elapsed:.6f} seconds") + + +def main(n_shells=20): + shells = [DummyShell() for _ in range(n_shells)] + + print(f"Benchmarking with {n_shells} shells\n") + + benchmark_intracule(shells) + benchmark_extracule(shells) + + +def main(n_shells=20): + shells = [DummyShell() for _ in range(n_shells)] + + print(f"Benchmarking with {n_shells} shells\n") + + start = time.perf_counter() + compute_intracule(shells) + print(f"intracule: {time.perf_counter() - start:.6f} seconds") + + start = time.perf_counter() + compute_extracule(shells) + print(f"extracule: {time.perf_counter() - start:.6f} seconds") + +if __name__ == "__main__": + main() \ No newline at end of file From a0070c8737166206ab7d0703cdbf5280567ab3bf Mon Sep 17 00:00:00 2001 From: kir943 Date: Fri, 13 Mar 2026 21:53:09 +0530 Subject: [PATCH 18/19] Improve documentation and apply black formatting --- gbasis/integrals/overlap_n.py | 230 +++++++++++++++++++++++++++++++--- tools/benchmark_intracule.py | 31 +++++ 2 files changed, 247 insertions(+), 14 deletions(-) create mode 100644 tools/benchmark_intracule.py diff --git a/gbasis/integrals/overlap_n.py b/gbasis/integrals/overlap_n.py index 87d2d018..62322631 100644 --- a/gbasis/integrals/overlap_n.py +++ b/gbasis/integrals/overlap_n.py @@ -1,3 +1,29 @@ +""" +Arbitrary-order N-center Gaussian overlap integrals. + +This module implements algorithms for computing overlap integrals +between Cartesian Gaussian basis functions located on different +centers. The implementation supports overlap integrals involving +an arbitrary number of Gaussian shells. + +The evaluation follows the Gaussian product theorem and Hermite +Gaussian expansion techniques commonly used in quantum chemistry +integral engines. + +References +---------- +Helgaker, T., Jørgensen, P., Olsen, J. +Molecular Electronic-Structure Theory, Wiley (2000). + +Szabo, A., Ostlund, N. S. +Modern Quantum Chemistry, Dover (1996). + +Obara, S., Saika, A. (1986) +Efficient recursive computation of molecular integrals over +Cartesian Gaussian functions. +Journal of Chemical Physics. +""" + import numpy as np from scipy.sparse import coo_matrix from itertools import product @@ -8,6 +34,50 @@ class PrimitiveNEngine: # Gaussian collapse @staticmethod def collapse_gaussians(alphas, centers): + """ + Collapse multiple primitive Gaussian functions into a single + effective Gaussian using the Gaussian product theorem. + + Given primitive Gaussian functions + + exp(-α_i |r - A_i|²) + + their product can be written as + + exp(-α_tot |r - P|²) * exp(term2 - term1) + + where + + α_tot = Σ α_i + P = (Σ α_i A_i) / α_tot + + This transformation simplifies the evaluation of Gaussian + integrals by reducing multiple Gaussians into a single + Gaussian centered at P. + + Parameters + ---------- + alphas : array_like + Exponents of the primitive Gaussian functions. + + centers : array_like + Cartesian coordinates of the Gaussian centers. + + Returns + ------- + alpha_tot : float + Total Gaussian exponent. + + P : ndarray + Product center of the collapsed Gaussian. + + prefactor : float + Exponential prefactor arising from Gaussian collapse. + Raises: + ------ + ValueError: + If the total Gaussian exponent is non-positive. + """ alphas = np.asarray(alphas, dtype=np.float64) centers = np.asarray(centers, dtype=np.float64) @@ -27,12 +97,35 @@ def collapse_gaussians(alphas, centers): return alpha_tot, P, prefactor # Pure binomial Hermite shift + @staticmethod def hermite_coefficients(l, PA): """ - Expand (x - A)^l about P: + Compute Hermite translation coefficients. + + This function expands a polynomial centered at A in terms of + powers centered at the Gaussian product center P: + + .. math:: + + (x - A)^l = \\sum_t E_t (x - P)^t - (x - A)^l = sum_t E_t (x - P)^t + These coefficients are used in the Hermite Gaussian formalism + for evaluating molecular integrals. + + Parameters + ---------- + l : int + Angular momentum order. + + PA : float + Distance between the product center P and the original + Gaussian center A. + + Returns + ------- + numpy.ndarray + Array containing the Hermite expansion coefficients. """ E = np.zeros(l + 1, dtype=np.float64) E[0] = 1.0 @@ -49,12 +142,33 @@ def hermite_coefficients(l, PA): return E # Gaussian moments + @staticmethod def gaussian_moments(alpha, max_order): """ - Compute: - ∫ (x-P)^k exp(-alpha (x-P)^2) dx - over (-∞,∞) + Compute Gaussian moment integrals. + + Evaluates integrals of the form + + .. math:: + + \\int_{-\\infty}^{\\infty} (x-P)^k e^{-\\alpha (x-P)^2} dx + + Only even moments are non-zero. Higher moments are computed + recursively from the zeroth moment. + + Parameters + ---------- + alpha : float + Gaussian exponent. + + max_order : int + Maximum moment order to compute. + + Returns + ------- + numpy.ndarray + Array containing Gaussian moments up to ``max_order``. """ moments = np.zeros(max_order + 1, dtype=np.float64) @@ -70,12 +184,55 @@ def gaussian_moments(alpha, max_order): # Full primitive N-center overlap @staticmethod def primitive_overlap(alphas, centers, angmoms): + """ + Compute primitive N-center Gaussian overlap integral. + + This function evaluates the overlap integral between N primitive + Gaussian basis functions using the Hermite Gaussian formalism. + + The primitive overlap integral is + + .. math:: + + S = \\int \\prod_i \\phi_i(\\mathbf{r}) \\, d\\mathbf{r} + where each primitive Gaussian basis function has the form + + .. math:: + + \\phi(\\mathbf{r}) = + (x-A_x)^l (y-A_y)^m (z-A_z)^n + e^{-\\alpha |\\mathbf{r}-A|^2} + + Parameters + ---------- + alphas : list + Gaussian exponents. + + centers : list + Cartesian coordinates of Gaussian centers. + + angmoms : list + Angular momentum tuples ``(l, m, n)``. + + Returns + ------- + float + Value of the primitive N-center overlap integral. + + Notes + ----- + The overlap integral factorizes into independent Cartesian + components (x, y, z). Each component is evaluated using + Gaussian moment integrals obtained from the Hermite expansion. + """ alpha_tot, P, prefactor = PrimitiveNEngine.collapse_gaussians(alphas, centers) result = prefactor # factorize into x, y, z + # Evaluate overlap as product of x, y, z integrals + for axis in range(3): # build total polynomial via convolution @@ -101,8 +258,26 @@ def primitive_overlap(alphas, centers, angmoms): # Screening Function def is_n_shell_overlap_screened(shells, tol=1e-12): """ - Conservative exponential upper-bound screening - for N-center contracted overlap. + Determine whether an N-center overlap integral can be skipped + using a conservative exponential bound. + + This screening estimates an upper bound for the magnitude of + the contracted overlap integral using Gaussian decay between + shell centers. If the bound is below the specified tolerance, + the integral is considered negligible. + + Parameters + ---------- + shells : list[GeneralizedContractionShell] + Shells involved in the overlap integral. + + tol : float + Screening tolerance. + + Returns + ------- + bool + True if the integral is negligible and can be skipped. """ alpha_mins = [np.min(shell.exps) for shell in shells] @@ -197,9 +372,15 @@ def contracted_n_overlap(shells): c = cart_indices[i] alpha = shell.exps[p] - coeff = shell.coeffs[p, m] - norm = shell.norm_prim_cart[c, p] - angmom = tuple(shell.angmom_components_cart[c]) + coeff = shell.coeffs[ + p, m + ] # contraction coefficient for primitive p in segment m + norm = shell.norm_prim_cart[ + c, p + ] # normalization factor for primitive p and Cartesian component c + angmom = tuple( + shell.angmom_components_cart[c] + ) # Cartesian angular momentum (lx, ly, lz) alphas.append(alpha) centers.append(shell.coord) @@ -224,18 +405,32 @@ def contracted_n_overlap(shells): def build_n_overlap_tensor(shells, tol=1e-12): """ - Build sparse N-center overlap tensor over shells. + Construct sparse N-center overlap tensor. + + This function builds the full N-center overlap tensor over + Gaussian basis functions. Each tensor element corresponds + to the overlap integral between a combination of contracted + Gaussian shells. + + To reduce computational cost, shell-level screening is + applied before evaluating the contracted integrals. Parameters ---------- shells : list[GeneralizedContractionShell] + List of Gaussian shells participating in the overlap + integral. + tol : float - Screening tolerance + Screening tolerance used to discard negligible + contributions. Returns ------- scipy.sparse.coo_matrix - Flattened sparse tensor of shape (total_ao^N, 1) + Sparse tensor containing the N-center overlap integrals. + The tensor is flattened so that each N-index integral + maps to a single row index. """ # Total AO dimension @@ -302,17 +497,24 @@ def arbitrary_order_overlap(shells, tol=1e-12): """ Compute arbitrary-order Gaussian overlap tensor. - This is the main public API for N-center overlap integrals. + This is the main public API for evaluating N-center overlap + integrals between Gaussian basis functions. + + Internally this function constructs the overlap tensor using + ``build_n_overlap_tensor``. Parameters ---------- shells : list[GeneralizedContractionShell] + List of Gaussian shells defining the basis functions. tol : float + Screening tolerance. Returns ------- scipy.sparse.coo_matrix + Sparse representation of the N-center overlap tensor. """ return build_n_overlap_tensor(shells, tol) diff --git a/tools/benchmark_intracule.py b/tools/benchmark_intracule.py new file mode 100644 index 00000000..779942cb --- /dev/null +++ b/tools/benchmark_intracule.py @@ -0,0 +1,31 @@ +import time +import numpy as np +from gbasis.integrals.density import compute_intracule + + +class DummyShell: + + def __init__(self): + self.exps = np.array([1.0]) + self.coeffs = np.array([[1.0]]) + self.norm_prim_cart = np.array([[1.0]]) + self.angmom_components_cart = [(0, 0, 0)] + self.angmom_components = [(0, 0, 0)] + self.coord = np.array([0.0, 0.0, 0.0]) + self.num_seg_cont = 1 + self.num_cart = 1 + + +def benchmark_intracule(n_shells=20): + + shells = [DummyShell() for _ in range(n_shells)] + + start = time.perf_counter() + compute_intracule(shells) + elapsed = time.perf_counter() - start + + print(f"compute_intracule benchmark: {elapsed:.6f} seconds for {n_shells} shells") + + +if __name__ == "__main__": + benchmark_intracule() \ No newline at end of file From da9e530ed4e14767e1ae591f027d617d787b97de Mon Sep 17 00:00:00 2001 From: kir943 Date: Fri, 13 Mar 2026 22:59:20 +0530 Subject: [PATCH 19/19] Improve documentation for density API functions --- gbasis/integrals/density.py | 105 ++++++++++++++++++++++++++++++++++-- 1 file changed, 101 insertions(+), 4 deletions(-) diff --git a/gbasis/integrals/density.py b/gbasis/integrals/density.py index 693a2601..b6a5514c 100644 --- a/gbasis/integrals/density.py +++ b/gbasis/integrals/density.py @@ -1,10 +1,65 @@ +""" +Density API built on top of the arbitrary-order Gaussian overlap engine. + +This module provides high-level routines for constructing matrices of +Gaussian overlap integrals between shells. These matrices serve as +building blocks in density-related calculations such as intracule and +extracule analysis. + +The overlap integrals are evaluated using the arbitrary-order Gaussian +overlap engine implemented in ``arbitrary_order_overlap``. +""" + import numpy as np from .overlap_n import arbitrary_order_overlap def compute_intracule(shells): - + """ + Compute pairwise intracule overlap matrix between Gaussian shells. + + This function constructs a matrix whose elements correspond to + two-center overlap integrals between Gaussian basis functions. + Each matrix element represents the overlap integral + + .. math:: + + S_{ij} = \\int \\phi_i(\\mathbf{r}) \\phi_j(\\mathbf{r}) \\, d\\mathbf{r} + + where :math:`\\phi_i` and :math:`\\phi_j` are Cartesian Gaussian basis + functions centered on different atoms or centers. + + The overlap integrals are evaluated using the + :func:`arbitrary_order_overlap` engine, which computes Gaussian + overlap tensors of arbitrary order. + + Parameters + ---------- + shells : list[GeneralizedContractionShell] + List of Gaussian shells defining the basis functions. + + Returns + ------- + numpy.ndarray + Array of shape ``(n, n)`` containing pairwise Gaussian + overlap integrals between shells. + + Notes + ----- + The intracule matrix is constructed by evaluating pairwise + two-center overlap integrals between Gaussian shells. Each + integral is extracted from the sparse tensor returned by the + arbitrary-order overlap engine. + + References + ---------- + Helgaker, T., Jørgensen, P., Olsen, J. + *Molecular Electronic-Structure Theory*, Wiley (2000). + + Szabo, A., Ostlund, N. S. + *Modern Quantum Chemistry*, Dover (1996). + """ n = len(shells) result = np.zeros((n, n)) @@ -14,7 +69,8 @@ def compute_intracule(shells): tensor = arbitrary_order_overlap([shells[i], shells[j]]) - # extract scalar value + # Extract scalar overlap value from sparse tensor returned by + # the arbitrary-order overlap engine. value = tensor.data[0] if tensor.nnz > 0 else 0.0 result[i, j] = value @@ -23,16 +79,57 @@ def compute_intracule(shells): def compute_extracule(shells): - + """ + Compute pairwise extracule overlap matrix between Gaussian shells. + + This function constructs a matrix of pairwise overlap integrals + between Gaussian shells using the arbitrary-order overlap engine. + + Each matrix element corresponds to the two-center Gaussian + overlap integral + + .. math:: + + S_{ij} = \\int \\phi_i(\\mathbf{r}) \\phi_j(\\mathbf{r}) \\, d\\mathbf{r} + + where :math:`\\phi_i` and :math:`\\phi_j` are Cartesian Gaussian + basis functions. + + Parameters + ---------- + shells : list[GeneralizedContractionShell] + List of Gaussian shells defining the basis functions. + + Returns + ------- + numpy.ndarray + Array of shape ``(n, n)`` containing pairwise Gaussian + overlap integrals between shells. + + Notes + ----- + This function provides an API wrapper around the + :func:`arbitrary_order_overlap` routine for density-related + calculations involving Gaussian basis functions. + This function currently constructs the same overlap matrix as + :func:`compute_intracule`, but is provided as a separate API + entry point for extracule-related density calculations. + + References + ---------- + Helgaker, T., Jørgensen, P., Olsen, J. + *Molecular Electronic-Structure Theory*, Wiley (2000). + """ n = len(shells) result = np.zeros((n, n)) for i in range(n): for j in range(n): - + # evaluate two-center Gaussian overlap tensor between shells i and j tensor = arbitrary_order_overlap([shells[i], shells[j]]) + # extract scalar overlap value from sparse tensor value = tensor.data[0] if tensor.nnz > 0 else 0.0 result[i, j] = value