From df2111605479b2ba34f018a84593f560068e8ff1 Mon Sep 17 00:00:00 2001 From: kir943 Date: Sun, 15 Feb 2026 02:33:08 +0530 Subject: [PATCH 1/7] =?UTF-8?q?Add=20N-center=20Gaussian=20overlap=20engin?= =?UTF-8?q?e=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 2/7] =?UTF-8?q?Remove=20experimental=20sparse=20builder=20?= =?UTF-8?q?test=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 3/7] 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 4/7] 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 5/7] 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 6/7] 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 7/7] 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