From 64f58cec244b6b4ed6cdea4714641ff83307d094 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Sat, 27 Sep 2025 00:36:47 -0400 Subject: [PATCH 01/14] Add screening and test to evaluate_stress_tensor --- gbasis/evals/stress_tensor.py | 26 ++++++++++- tests/test_stress_tensor.py | 83 ++++++++++++++++++++++++++++++----- 2 files changed, 96 insertions(+), 13 deletions(-) diff --git a/gbasis/evals/stress_tensor.py b/gbasis/evals/stress_tensor.py index 6a33ac51..8fd4b46c 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, @@ -8,7 +9,16 @@ # TODO: need to be tested against reference -def evaluate_stress_tensor(one_density_matrix, basis, points, alpha=1, beta=0, transform=None): +def evaluate_stress_tensor( + one_density_matrix, + basis, + points, + alpha=1, + beta=0, + transform=None, + screen_basis=True, + tol_screen=1e-8, +): r"""Return the stress tensor evaluated at the given coordinates. Stress tensor is defined here as: @@ -70,6 +80,14 @@ def evaluate_stress_tensor(one_density_matrix, basis, points, alpha=1, beta=0, t beta : {int, float} Second parameter of the stress tensor. Default value is 0. + screen_basis : bool, optional + Whether to screen out points with negligible contributions. Default value is True + (enable screening). + tol_screen : float + Screening tolerance for excluding evaluations. Points with values below this tolerance + will not be evaluated (they will be set to zero). Internal computed quantities that + affect the results below this tolerance will also be ignored to speed up the + evaluation. Default value is 1e-8. Returns ------- @@ -99,6 +117,8 @@ def evaluate_stress_tensor(one_density_matrix, basis, points, alpha=1, beta=0, t basis, points, transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, ) if alpha != 1: output[i, j] += (1 - alpha) * evaluate_deriv_reduced_density_matrix( @@ -108,6 +128,8 @@ def evaluate_stress_tensor(one_density_matrix, basis, points, alpha=1, beta=0, t basis, points, transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, ) if i == j and beta != 0: output[i, j] -= ( @@ -118,6 +140,8 @@ def evaluate_stress_tensor(one_density_matrix, basis, points, alpha=1, beta=0, t basis, points, transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, ) ) output[j, i] = output[i, j] diff --git a/tests/test_stress_tensor.py b/tests/test_stress_tensor.py index a9ef3f56..3361a584 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, @@ -15,7 +16,9 @@ from utils import find_datafile, HortonContractions -def test_evaluate_stress_tensor(): +@pytest.mark.parametrize("screen_basis", [True, False]) +@pytest.mark.parametrize("tol_screen", [1e-8]) +def test_evaluate_stress_tensor(screen_basis, tol_screen): """Test gbasis.evals.stress_tensor.evaluate_stress_tensor.""" basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem")) coords = np.array([[0, 0, 0]]) @@ -29,10 +32,47 @@ def test_evaluate_stress_tensor(): with pytest.raises(TypeError): evaluate_stress_tensor(np.identity(40), basis, points, 1, 0j, np.identity(40)) - test_a = evaluate_stress_tensor(np.identity(40), basis, points, 0, 0, np.identity(40)) - test_b = evaluate_stress_tensor(np.identity(40), basis, points, 1, 0, np.identity(40)) - test_c = evaluate_stress_tensor(np.identity(40), basis, points, 1, 2, np.identity(40)) - test_d = evaluate_stress_tensor(np.identity(40), basis, points, 0.5, 2, np.identity(40)) + test_a = evaluate_stress_tensor( + np.identity(40), + basis, + points, + 0, + 0, + np.identity(40), + screen_basis=screen_basis, + tol_screen=tol_screen, + ) + test_b = evaluate_stress_tensor( + np.identity(40), + basis, + points, + 1, + 0, + np.identity(40), + screen_basis=screen_basis, + tol_screen=tol_screen, + ) + test_c = evaluate_stress_tensor( + np.identity(40), + basis, + points, + 1, + 2, + np.identity(40), + screen_basis=screen_basis, + tol_screen=tol_screen, + ) + test_d = evaluate_stress_tensor( + np.identity(40), + basis, + points, + 0.5, + 2, + np.identity(40), + screen_basis=screen_basis, + tol_screen=tol_screen, + ) + # compute reference without screening for i in range(3): for j in range(3): orders_i = np.array([0, 0, 0]) @@ -41,10 +81,22 @@ def test_evaluate_stress_tensor(): orders_j[j] += 1 temp1 = evaluate_deriv_reduced_density_matrix( - orders_i, orders_j, np.identity(40), basis, points, np.identity(40) + orders_i, + orders_j, + np.identity(40), + basis, + points, + np.identity(40), + screen_basis=False, ) temp2 = evaluate_deriv_reduced_density_matrix( - orders_j, orders_i, np.identity(40), basis, points, np.identity(40) + orders_j, + orders_i, + np.identity(40), + basis, + points, + np.identity(40), + screen_basis=False, ) temp3 = evaluate_deriv_reduced_density_matrix( orders_i + orders_j, @@ -53,6 +105,7 @@ def test_evaluate_stress_tensor(): basis, points, np.identity(40), + screen_basis=False, ) temp4 = evaluate_deriv_reduced_density_matrix( orders_i + orders_j, @@ -61,16 +114,22 @@ def test_evaluate_stress_tensor(): basis, points, np.identity(40), + screen_basis=False, ) if i == j: - temp5 = evaluate_density_laplacian(np.identity(40), basis, points, np.identity(40)) + temp5 = evaluate_density_laplacian( + np.identity(40), basis, points, np.identity(40), screen_basis=False + ) else: temp5 = 0 - assert np.allclose(test_a[:, i, j], 0.5 * temp3 + 0.5 * temp4) - assert np.allclose(test_b[:, i, j], -0.5 * temp1 - 0.5 * temp2) - assert np.allclose(test_c[:, i, j], -0.5 * temp1 - 0.5 * temp2 - temp5) + # check that non screened reference matches screened result within tol_screen + assert np.allclose(test_a[:, i, j], 0.5 * temp3 + 0.5 * temp4, atol=tol_screen) + assert np.allclose(test_b[:, i, j], -0.5 * temp1 - 0.5 * temp2, atol=tol_screen) + assert np.allclose(test_c[:, i, j], -0.5 * temp1 - 0.5 * temp2 - temp5, atol=tol_screen) assert np.allclose( - test_d[:, i, j], -0.25 * temp1 - 0.25 * temp2 + 0.25 * temp3 + 0.25 * temp4 - temp5 + test_d[:, i, j], + -0.25 * temp1 - 0.25 * temp2 + 0.25 * temp3 + 0.25 * temp4 - temp5, + atol=tol_screen, ) From 5d0d6221707053f163c1951f3f0ac98d8f953590 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Sat, 25 Oct 2025 02:14:24 -0400 Subject: [PATCH 02/14] Set screening of integrals to False by default --- gbasis/integrals/angular_momentum.py | 4 ++-- gbasis/integrals/kinetic_energy.py | 4 ++-- gbasis/integrals/moment.py | 2 +- gbasis/integrals/momentum.py | 4 ++-- gbasis/integrals/overlap.py | 4 ++-- gbasis/integrals/overlap_asymm.py | 2 +- 6 files changed, 10 insertions(+), 10 deletions(-) diff --git a/gbasis/integrals/angular_momentum.py b/gbasis/integrals/angular_momentum.py index 9f6503bc..7a8413ca 100644 --- a/gbasis/integrals/angular_momentum.py +++ b/gbasis/integrals/angular_momentum.py @@ -60,7 +60,7 @@ class AngularMomentumIntegral(BaseTwoIndexSymmetric): @staticmethod def construct_array_contraction( - contractions_one, contractions_two, screen_basis=True, tol_screen=1e-8 + contractions_one, contractions_two, screen_basis=False, tol_screen=1e-8 ): """Return the integrals over the angular momentum operator of the given contractions. @@ -183,7 +183,7 @@ def construct_array_contraction( return -1j * np.transpose(output, (3, 2, 4, 1, 0)) -def angular_momentum_integral(basis, transform=None, screen_basis=True, tol_screen=1e-8): +def angular_momentum_integral(basis, transform=None, screen_basis=False, tol_screen=1e-8): r"""Return the integral over :math:`hat{L}` of the given basis set. .. math:: diff --git a/gbasis/integrals/kinetic_energy.py b/gbasis/integrals/kinetic_energy.py index 3276abee..bed75d8b 100644 --- a/gbasis/integrals/kinetic_energy.py +++ b/gbasis/integrals/kinetic_energy.py @@ -55,7 +55,7 @@ class KineticEnergyIntegral(BaseTwoIndexSymmetric): @staticmethod def construct_array_contraction( - contractions_one, contractions_two, screen_basis=True, tol_screen=1e-8 + contractions_one, contractions_two, screen_basis=False, tol_screen=1e-8 ): """Return the evaluations of the given contractions at the given coordinates. @@ -146,7 +146,7 @@ def construct_array_contraction( return -0.5 * np.sum(output, axis=0) -def kinetic_energy_integral(basis, transform=None, screen_basis=True, tol_screen=1e-8): +def kinetic_energy_integral(basis, transform=None, screen_basis=False, tol_screen=1e-8): r"""Return kinetic energy integral of the given basis set. .. math:: diff --git a/gbasis/integrals/moment.py b/gbasis/integrals/moment.py index 682c95d3..49bdcf80 100644 --- a/gbasis/integrals/moment.py +++ b/gbasis/integrals/moment.py @@ -188,7 +188,7 @@ class if the first four indices correspond to the segmented contraction and the def moment_integral( - basis, moment_coord, moment_orders, transform=None, screen_basis=True, tol_screen=1e-8 + basis, moment_coord, moment_orders, transform=None, screen_basis=False, tol_screen=1e-8 ): r"""Return moment integral of the given basis set. diff --git a/gbasis/integrals/momentum.py b/gbasis/integrals/momentum.py index 86a2fc37..e8496992 100644 --- a/gbasis/integrals/momentum.py +++ b/gbasis/integrals/momentum.py @@ -57,7 +57,7 @@ class MomentumIntegral(BaseTwoIndexSymmetric): @staticmethod def construct_array_contraction( - contractions_one, contractions_two, screen_basis=True, tol_screen=1e-8 + contractions_one, contractions_two, screen_basis=False, tol_screen=1e-8 ): """Return the integrals over the momentum operator of the given contractions. @@ -140,7 +140,7 @@ def construct_array_contraction( return -1j * np.transpose(output, (1, 2, 3, 4, 0)) -def momentum_integral(basis, transform=None, screen_basis=True, tol_screen=1e-8): +def momentum_integral(basis, transform=None, screen_basis=False, tol_screen=1e-8): r"""Return integral over momentum operator of the given basis set. .. math:: diff --git a/gbasis/integrals/overlap.py b/gbasis/integrals/overlap.py index 7ee3c9ce..0b767170 100644 --- a/gbasis/integrals/overlap.py +++ b/gbasis/integrals/overlap.py @@ -53,7 +53,7 @@ class Overlap(BaseTwoIndexSymmetric): @staticmethod def construct_array_contraction( - contractions_one, contractions_two, screen_basis=True, tol_screen=1e-8 + contractions_one, contractions_two, screen_basis=False, tol_screen=1e-8 ): """Return the evaluations of the given contractions at the given coordinates. @@ -134,7 +134,7 @@ def construct_array_contraction( )[0] -def overlap_integral(basis, transform=None, screen_basis=True, tol_screen=1e-8): +def overlap_integral(basis, transform=None, screen_basis=False, tol_screen=1e-8): r"""Return overlap integral of the given basis set. .. math:: diff --git a/gbasis/integrals/overlap_asymm.py b/gbasis/integrals/overlap_asymm.py index 0a4c80a6..7e96d4a5 100644 --- a/gbasis/integrals/overlap_asymm.py +++ b/gbasis/integrals/overlap_asymm.py @@ -65,7 +65,7 @@ class OverlapAsymmetric(BaseTwoIndexAsymmetric): def overlap_integral_asymmetric( - basis_one, basis_two, transform_one=None, transform_two=None, screen_basis=True, tol_screen=1e-8 + basis_one, basis_two, transform_one=None, transform_two=None, screen_basis=False, tol_screen=1e-8 ): r"""Return overlap integrals between two basis sets. From 8e6bb37237973a522717876b62d93602ce8f5260 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Sat, 27 Sep 2025 00:46:34 -0400 Subject: [PATCH 03/14] Add screening and test to evaluate_ehrenfest_force --- gbasis/evals/stress_tensor.py | 29 ++++++++++- tests/test_stress_tensor.py | 90 ++++++++++++++++++++++++++++++----- 2 files changed, 104 insertions(+), 15 deletions(-) diff --git a/gbasis/evals/stress_tensor.py b/gbasis/evals/stress_tensor.py index 8fd4b46c..12125601 100644 --- a/gbasis/evals/stress_tensor.py +++ b/gbasis/evals/stress_tensor.py @@ -16,7 +16,7 @@ def evaluate_stress_tensor( alpha=1, beta=0, transform=None, - screen_basis=True, + screen_basis=False, tol_screen=1e-8, ): r"""Return the stress tensor evaluated at the given coordinates. @@ -149,7 +149,16 @@ def evaluate_stress_tensor( # TODO: need to be tested against reference -def evaluate_ehrenfest_force(one_density_matrix, basis, points, alpha=1, beta=0, transform=None): +def evaluate_ehrenfest_force( + one_density_matrix, + basis, + points, + alpha=1, + beta=0, + transform=None, + screen_basis=False, + tol_screen=1e-8, +): r"""Return the Ehrenfest force. Ehrenfest force is the negative of the divergence of the stress tensor: @@ -207,6 +216,14 @@ def evaluate_ehrenfest_force(one_density_matrix, basis, points, alpha=1, beta=0, beta : {int, float} Second parameter of the stress tensor. Default value is 0. + screen_basis : bool, optional + Whether to screen out points with negligible contributions. Default value is True + (enable screening). + tol_screen : float + Screening tolerance for excluding evaluations. Points with values below this tolerance + will not be evaluated (they will be set to zero). Internal computed quantities that + affect the results below this tolerance will also be ignored to speed up the + evaluation. Default value is 1e-8. Returns ------- @@ -235,6 +252,8 @@ def evaluate_ehrenfest_force(one_density_matrix, basis, points, alpha=1, beta=0, basis, points, transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, ) if alpha != 1: output[i] -= (1 - alpha) * evaluate_deriv_reduced_density_matrix( @@ -244,6 +263,8 @@ def evaluate_ehrenfest_force(one_density_matrix, basis, points, alpha=1, beta=0, basis, points, transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, ) if alpha != 0.5: output[i] -= (1 - 2 * alpha) * evaluate_deriv_reduced_density_matrix( @@ -253,6 +274,8 @@ def evaluate_ehrenfest_force(one_density_matrix, basis, points, alpha=1, beta=0, basis, points, transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, ) if beta != 0: output[i] += ( @@ -264,6 +287,8 @@ def evaluate_ehrenfest_force(one_density_matrix, basis, points, alpha=1, beta=0, basis, points, transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, ) ) return output.T diff --git a/tests/test_stress_tensor.py b/tests/test_stress_tensor.py index 3361a584..6eda85a3 100644 --- a/tests/test_stress_tensor.py +++ b/tests/test_stress_tensor.py @@ -133,7 +133,9 @@ def test_evaluate_stress_tensor(screen_basis, tol_screen): ) -def test_evaluate_ehrenfest_force(): +@pytest.mark.parametrize("screen_basis", [True, False]) +@pytest.mark.parametrize("tol_screen", [1e-8]) +def test_evaluate_ehrenfest_force(screen_basis, tol_screen): """Test gbasis.evals.stress_tensor.evaluate_ehrenfest_force.""" basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem")) coords = np.array([[0, 0, 0]]) @@ -147,16 +149,53 @@ def test_evaluate_ehrenfest_force(): with pytest.raises(TypeError): evaluate_ehrenfest_force(np.identity(40), basis, points, 1, 0j, np.identity(40)) - test_a = evaluate_ehrenfest_force(np.identity(40), basis, points, 0, 0, np.identity(40)) - test_b = evaluate_ehrenfest_force(np.identity(40), basis, points, 1, 0, np.identity(40)) - test_c = evaluate_ehrenfest_force(np.identity(40), basis, points, 0.5, 0, np.identity(40)) - test_d = evaluate_ehrenfest_force(np.identity(40), basis, points, 0, 2, np.identity(40)) + test_a = evaluate_ehrenfest_force( + np.identity(40), + basis, + points, + 0, + 0, + np.identity(40), + screen_basis=screen_basis, + tol_screen=tol_screen, + ) + test_b = evaluate_ehrenfest_force( + np.identity(40), + basis, + points, + 1, + 0, + np.identity(40), + screen_basis=screen_basis, + tol_screen=tol_screen, + ) + test_c = evaluate_ehrenfest_force( + np.identity(40), + basis, + points, + 0.5, + 0, + np.identity(40), + screen_basis=screen_basis, + tol_screen=tol_screen, + ) + test_d = evaluate_ehrenfest_force( + np.identity(40), + basis, + points, + 0, + 2, + np.identity(40), + screen_basis=screen_basis, + tol_screen=tol_screen, + ) for j in range(3): ref_a = np.zeros(points.shape[0]) ref_b = np.zeros(points.shape[0]) ref_c = np.zeros(points.shape[0]) ref_d = np.zeros(points.shape[0]) + # compute reference without screening orders_j = np.array([0, 0, 0]) orders_j[j] += 1 for i in range(3): @@ -164,10 +203,22 @@ def test_evaluate_ehrenfest_force(): orders_i[i] += 1 temp1 = evaluate_deriv_reduced_density_matrix( - 2 * orders_i, orders_j, np.identity(40), basis, points, np.identity(40) + 2 * orders_i, + orders_j, + np.identity(40), + basis, + points, + np.identity(40), + screen_basis=False, ) temp2 = evaluate_deriv_reduced_density_matrix( - orders_i, orders_i + orders_j, np.identity(40), basis, points, np.identity(40) + orders_i, + orders_i + orders_j, + np.identity(40), + basis, + points, + np.identity(40), + screen_basis=False, ) temp3 = evaluate_deriv_reduced_density_matrix( 2 * orders_i + orders_j, @@ -176,22 +227,35 @@ def test_evaluate_ehrenfest_force(): basis, points, np.identity(40), + screen_basis=False, ) temp4 = evaluate_deriv_reduced_density_matrix( - orders_i + orders_j, orders_i, np.identity(40), basis, points, np.identity(40) + orders_i + orders_j, + orders_i, + np.identity(40), + basis, + points, + np.identity(40), + screen_basis=False, ) temp5 = evaluate_deriv_density( - 2 * orders_i + orders_j, np.identity(40), basis, points, np.identity(40) + 2 * orders_i + orders_j, + np.identity(40), + basis, + points, + np.identity(40), + screen_basis=False, ) ref_a += temp3 + temp4 ref_b += -temp1 - temp2 ref_c += -0.5 * temp1 - 0.5 * temp2 + 0.5 * temp3 + 0.5 * temp4 ref_d += temp3 + temp4 - temp5 - assert np.allclose(test_a[:, j], -ref_a) - assert np.allclose(test_b[:, j], -ref_b) - assert np.allclose(test_c[:, j], -ref_c) - assert np.allclose(test_d[:, j], -ref_d) + # check that non screened reference matches screened result within tol_screen + assert np.allclose(test_a[:, j], -ref_a, atol=tol_screen) + assert np.allclose(test_b[:, j], -ref_b, atol=tol_screen) + assert np.allclose(test_c[:, j], -ref_c, atol=tol_screen) + assert np.allclose(test_d[:, j], -ref_d, atol=tol_screen) def test_evaluate_ehrenfest_hessian(): From 72ab5d34dcff1359d7e2d37b535791d6e163213e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Sat, 27 Sep 2025 00:55:21 -0400 Subject: [PATCH 04/14] Add screening and test to evaluate_ehrenfest_hessian --- gbasis/evals/stress_tensor.py | 24 ++++++++++ tests/test_stress_tensor.py | 85 ++++++++++++++++++++++++++++++----- 2 files changed, 97 insertions(+), 12 deletions(-) diff --git a/gbasis/evals/stress_tensor.py b/gbasis/evals/stress_tensor.py index 12125601..6d8dfd79 100644 --- a/gbasis/evals/stress_tensor.py +++ b/gbasis/evals/stress_tensor.py @@ -303,6 +303,8 @@ def evaluate_ehrenfest_hessian( beta=0, transform=None, symmetric=False, + screen_basis=False, + tol_screen=1e-8, ): r"""Return the Ehrenfest Hessian. @@ -375,6 +377,14 @@ def evaluate_ehrenfest_hessian( Flag for symmetrizing the Hessian. If True, then the Hessian is symmetrized by averaging it with its transpose. Default value is False. + screen_basis : bool, optional + Whether to screen out points with negligible contributions. Default value is True + (enable screening). + tol_screen : float + Screening tolerance for excluding evaluations. Points with values below this tolerance + will not be evaluated (they will be set to zero). Internal computed quantities that + affect the results below this tolerance will also be ignored to speed up the + evaluation. Default value is 1e-8. Returns ------- @@ -405,6 +415,8 @@ def evaluate_ehrenfest_hessian( basis, points, transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, ) output[i, j] += alpha * evaluate_deriv_reduced_density_matrix( 2 * orders_one, @@ -413,6 +425,8 @@ def evaluate_ehrenfest_hessian( basis, points, transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, ) if alpha != 1: output[i, j] -= (1 - alpha) * evaluate_deriv_reduced_density_matrix( @@ -422,6 +436,8 @@ def evaluate_ehrenfest_hessian( basis, points, transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, ) output[i, j] -= (1 - alpha) * evaluate_deriv_reduced_density_matrix( 2 * orders_one + orders_two, @@ -430,6 +446,8 @@ def evaluate_ehrenfest_hessian( basis, points, transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, ) if alpha != 0.5: output[i, j] -= (1 - 2 * alpha) * evaluate_deriv_reduced_density_matrix( @@ -439,6 +457,8 @@ def evaluate_ehrenfest_hessian( basis, points, transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, ) output[i, j] -= (1 - 2 * alpha) * evaluate_deriv_reduced_density_matrix( orders_one + orders_two, @@ -447,6 +467,8 @@ def evaluate_ehrenfest_hessian( basis, points, transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, ) if beta != 0: output[i, j] += ( @@ -458,6 +480,8 @@ def evaluate_ehrenfest_hessian( basis, points, transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, ) ) if symmetric: diff --git a/tests/test_stress_tensor.py b/tests/test_stress_tensor.py index 6eda85a3..da45c979 100644 --- a/tests/test_stress_tensor.py +++ b/tests/test_stress_tensor.py @@ -258,7 +258,9 @@ def test_evaluate_ehrenfest_force(screen_basis, tol_screen): assert np.allclose(test_d[:, j], -ref_d, atol=tol_screen) -def test_evaluate_ehrenfest_hessian(): +@pytest.mark.parametrize("screen_basis", [True, False]) +@pytest.mark.parametrize("tol_screen", [1e-8]) +def test_evaluate_ehrenfest_hessian(screen_basis, tol_screen): """Test gbasis.evals.stress_tensor.evaluate_ehrenfest_hessian.""" basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem")) coords = np.array([[0, 0, 0]]) @@ -274,10 +276,48 @@ def test_evaluate_ehrenfest_hessian(): with pytest.raises(TypeError): evaluate_ehrenfest_hessian(np.identity(40), basis, points, 1, 0j, np.identity(40)) - test_a = evaluate_ehrenfest_hessian(np.identity(40), basis, points, 0, 0, np.identity(40)) - test_b = evaluate_ehrenfest_hessian(np.identity(40), basis, points, 1, 0, np.identity(40)) - test_c = evaluate_ehrenfest_hessian(np.identity(40), basis, points, 0.5, 0, np.identity(40)) - test_d = evaluate_ehrenfest_hessian(np.identity(40), basis, points, 0, 2, np.identity(40)) + test_a = evaluate_ehrenfest_hessian( + np.identity(40), + basis, + points, + 0, + 0, + np.identity(40), + screen_basis=screen_basis, + tol_screen=tol_screen, + ) + test_b = evaluate_ehrenfest_hessian( + np.identity(40), + basis, + points, + 1, + 0, + np.identity(40), + screen_basis=screen_basis, + tol_screen=tol_screen, + ) + test_c = evaluate_ehrenfest_hessian( + np.identity(40), + basis, + points, + 0.5, + 0, + np.identity(40), + screen_basis=screen_basis, + tol_screen=tol_screen, + ) + test_d = evaluate_ehrenfest_hessian( + np.identity(40), + basis, + points, + 0, + 2, + np.identity(40), + screen_basis=screen_basis, + tol_screen=tol_screen, + ) + + # compute reference without screening for j in range(3): for k in range(3): ref_a = np.zeros(points.shape[0]) @@ -300,6 +340,7 @@ def test_evaluate_ehrenfest_hessian(): basis, points, np.identity(40), + screen_basis=False, ) temp2 = evaluate_deriv_reduced_density_matrix( 2 * orders_i, @@ -308,6 +349,7 @@ def test_evaluate_ehrenfest_hessian(): basis, points, np.identity(40), + screen_basis=False, ) temp3 = evaluate_deriv_reduced_density_matrix( orders_i + orders_k, @@ -316,6 +358,7 @@ def test_evaluate_ehrenfest_hessian(): basis, points, np.identity(40), + screen_basis=False, ) temp4 = evaluate_deriv_reduced_density_matrix( orders_i, @@ -324,6 +367,7 @@ def test_evaluate_ehrenfest_hessian(): basis, points, np.identity(40), + screen_basis=False, ) temp5 = evaluate_deriv_reduced_density_matrix( 2 * orders_i + orders_j + orders_k, @@ -332,6 +376,7 @@ def test_evaluate_ehrenfest_hessian(): basis, points, np.identity(40), + screen_basis=False, ) temp6 = evaluate_deriv_reduced_density_matrix( 2 * orders_i + orders_j, @@ -340,6 +385,7 @@ def test_evaluate_ehrenfest_hessian(): basis, points, np.identity(40), + screen_basis=False, ) temp7 = evaluate_deriv_reduced_density_matrix( orders_i + orders_j + orders_k, @@ -348,6 +394,7 @@ def test_evaluate_ehrenfest_hessian(): basis, points, np.identity(40), + screen_basis=False, ) temp8 = evaluate_deriv_reduced_density_matrix( orders_i + orders_j, @@ -356,6 +403,7 @@ def test_evaluate_ehrenfest_hessian(): basis, points, np.identity(40), + screen_basis=False, ) temp9 = evaluate_deriv_density( 2 * orders_i + orders_j + orders_k, @@ -363,6 +411,7 @@ def test_evaluate_ehrenfest_hessian(): basis, points, np.identity(40), + screen_basis=False, ) ref_a += temp5 + temp6 + temp7 + temp8 @@ -378,18 +427,30 @@ def test_evaluate_ehrenfest_hessian(): + 0.5 * temp8 ) ref_d += temp3 + temp4 + temp5 + temp6 - temp9 - assert np.allclose(test_a[:, j, k], -ref_a) - assert np.allclose(test_b[:, j, k], -ref_b) - assert np.allclose(test_c[:, j, k], -ref_c) - assert np.allclose(test_d[:, j, k], -ref_d) + assert np.allclose(test_a[:, j, k], -ref_a, atol=tol_screen) + assert np.allclose(test_b[:, j, k], -ref_b, atol=tol_screen) + assert np.allclose(test_c[:, j, k], -ref_c, atol=tol_screen) + assert np.allclose(test_d[:, j, k], -ref_d, atol=tol_screen) + # check symmetry if requested assert np.allclose( evaluate_ehrenfest_hessian( - np.identity(40), basis, points, 0, 0, np.identity(40), symmetric=True + np.identity(40), + basis, + points, + 0, + 0, + np.identity(40), + symmetric=True, + screen_basis=False, ), ( - evaluate_ehrenfest_hessian(np.identity(40), basis, points, 0, 0, np.identity(40)) + evaluate_ehrenfest_hessian( + np.identity(40), basis, points, 0, 0, np.identity(40), screen_basis=False + ) + np.swapaxes( - evaluate_ehrenfest_hessian(np.identity(40), basis, points, 0, 0, np.identity(40)), + evaluate_ehrenfest_hessian( + np.identity(40), basis, points, 0, 0, np.identity(40), screen_basis=False + ), 1, 2, ) From 6644798ae89077efd81246fbb51ed2beb8ce6f44 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Sat, 27 Sep 2025 00:55:21 -0400 Subject: [PATCH 05/14] Add screening and test to evaluate_ehrenfest_hessian --- gbasis/evals/stress_tensor.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/gbasis/evals/stress_tensor.py b/gbasis/evals/stress_tensor.py index 6d8dfd79..76994c07 100644 --- a/gbasis/evals/stress_tensor.py +++ b/gbasis/evals/stress_tensor.py @@ -81,8 +81,7 @@ def evaluate_stress_tensor( Second parameter of the stress tensor. Default value is 0. screen_basis : bool, optional - Whether to screen out points with negligible contributions. Default value is True - (enable screening). + Whether to screen out points with negligible contributions. Default value is False. tol_screen : float Screening tolerance for excluding evaluations. Points with values below this tolerance will not be evaluated (they will be set to zero). Internal computed quantities that @@ -217,8 +216,7 @@ def evaluate_ehrenfest_force( Second parameter of the stress tensor. Default value is 0. screen_basis : bool, optional - Whether to screen out points with negligible contributions. Default value is True - (enable screening). + Whether to screen out points with negligible contributions. Default value is False. tol_screen : float Screening tolerance for excluding evaluations. Points with values below this tolerance will not be evaluated (they will be set to zero). Internal computed quantities that @@ -378,8 +376,7 @@ def evaluate_ehrenfest_hessian( If True, then the Hessian is symmetrized by averaging it with its transpose. Default value is False. screen_basis : bool, optional - Whether to screen out points with negligible contributions. Default value is True - (enable screening). + Whether to screen out points with negligible contributions. Default value is False. tol_screen : float Screening tolerance for excluding evaluations. Points with values below this tolerance will not be evaluated (they will be set to zero). Internal computed quantities that From c4662c4cf3fb2c987dd80c3a7e73a73c31f81e98 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Sat, 25 Oct 2025 03:11:07 -0400 Subject: [PATCH 06/14] Set screening default to False and update docstrings --- gbasis/evals/density.py | 40 +++++++++++----------------- gbasis/evals/eval.py | 10 +++---- gbasis/evals/eval_deriv.py | 8 +++--- gbasis/integrals/angular_momentum.py | 4 +-- gbasis/integrals/kinetic_energy.py | 4 +-- gbasis/integrals/moment.py | 6 ++--- gbasis/integrals/momentum.py | 4 +-- gbasis/integrals/overlap.py | 4 +-- gbasis/integrals/overlap_asymm.py | 2 +- 9 files changed, 36 insertions(+), 46 deletions(-) diff --git a/gbasis/evals/density.py b/gbasis/evals/density.py index a72772ec..0cdd96e1 100644 --- a/gbasis/evals/density.py +++ b/gbasis/evals/density.py @@ -69,7 +69,7 @@ def evaluate_density( points, transform=None, threshold=1.0e-8, - screen_basis=True, + screen_basis=False, tol_screen=1e-8, ): r"""Return the density of the given basis set at the given points. @@ -104,8 +104,7 @@ def evaluate_density( The absolute value below which negative density values are acceptable. Any negative density value with an absolute value smaller than this threshold will be set to zero. screen_basis : bool, optional - Whether to screen out points with negligible contributions. Default value is True - (enable screening). + Whether to screen out points with negligible contributions. Default value is False. tol_screen : float Screening tolerance for excluding evaluations. Points with values below this tolerance will not be evaluated (they will be set to zero). Internal computed quantities that @@ -137,7 +136,7 @@ def evaluate_deriv_reduced_density_matrix( points, transform=None, deriv_type="general", - screen_basis=True, + screen_basis=False, tol_screen=1e-8, ): r"""Return the derivative of the first-order reduced density matrix at the given points. @@ -192,8 +191,7 @@ def evaluate_deriv_reduced_density_matrix( and "direct" makes reference to specific implementation of first and second order derivatives for generalized contraction (_eval_first_second_order_deriv_contractions()). screen_basis : bool, optional - Whether to screen out points with negligible contributions. Default value is True - (enable screening). + Whether to screen out points with negligible contributions. Default value is False. tol_screen : float Screening tolerance for excluding evaluations. Points with values below this tolerance will not be evaluated (they will be set to zero). Internal computed quantities that @@ -241,7 +239,7 @@ def evaluate_deriv_density( points, transform=None, deriv_type="general", - screen_basis=True, + screen_basis=False, tol_screen=1e-8, ): r"""Return the derivative of density of the given transformed basis set at the given points. @@ -291,8 +289,7 @@ def evaluate_deriv_density( and "direct" makes reference to specific implementation of first and second order derivatives for generalized contraction (_eval_first_second_order_deriv_contractions()). screen_basis : bool, optional - Whether to screen out points with negligible contributions. Default value is True - (enable screening). + Whether to screen out points with negligible contributions. Default value is False. tol_screen : float Screening tolerance for excluding evaluations. Points with values below this tolerance will not be evaluated (they will be set to zero). Internal computed quantities that @@ -357,7 +354,7 @@ def evaluate_density_gradient( points, transform=None, deriv_type="general", - screen_basis=True, + screen_basis=False, tol_screen=1e-8, ): r"""Return the gradient of the density evaluated at the given points. @@ -397,8 +394,7 @@ def evaluate_density_gradient( and "direct" makes reference to specific implementation of first and second order derivatives for generalized contraction (_eval_first_second_order_deriv_contractions()). screen_basis : bool, optional - Whether to screen out points with negligible contributions. Default value is True - (enable screening). + Whether to screen out points with negligible contributions. Default value is False. tol_screen : float Screening tolerance for excluding evaluations. Points with values below this tolerance will not be evaluated (they will be set to zero). Internal computed quantities that @@ -448,7 +444,7 @@ def evaluate_density_laplacian( points, transform=None, deriv_type="general", - screen_basis=True, + screen_basis=False, tol_screen=1e-8, ): r"""Return the Laplacian of the density evaluated at the given points. @@ -486,8 +482,7 @@ def evaluate_density_laplacian( and "direct" makes reference to specific implementation of first and second order derivatives for generalized contraction (_eval_first_second_order_deriv_contractions()). screen_basis : bool, optional - Whether to screen out points with negligible contributions. Default value is True - (enable screening). + Whether to screen out points with negligible contributions. Default value is False. tol_screen : float Screening tolerance for excluding evaluations. Points with values below this tolerance will not be evaluated (they will be set to zero). Internal computed quantities that @@ -565,7 +560,7 @@ def evaluate_density_hessian( points, transform=None, deriv_type="general", - screen_basis=True, + screen_basis=False, tol_screen=1e-8, ): r"""Return the Hessian of the density evaluated at the given points. @@ -611,8 +606,7 @@ def evaluate_density_hessian( and "direct" makes reference to specific implementation of first and second order derivatives for generalized contraction (_eval_first_second_order_deriv_contractions()). screen_basis : bool, optional - Whether to screen out points with negligible contributions. Default value is True - (enable screening). + Whether to screen out points with negligible contributions. Default value is False. tol_screen : float Screening tolerance for excluding evaluations. Points with values below this tolerance will not be evaluated (they will be set to zero). Internal computed quantities that @@ -719,7 +713,7 @@ def evaluate_posdef_kinetic_energy_density( transform=None, deriv_type="general", threshold=1.0e-8, - screen_basis=True, + screen_basis=False, tol_screen=1e-8, ): r"""Return evaluations of positive definite kinetic energy density at the given points. @@ -765,8 +759,7 @@ def evaluate_posdef_kinetic_energy_density( The absolute value below which negative density values are acceptable. Any negative density value with an absolute value smaller than this threshold will be set to zero. screen_basis : bool, optional - Whether to screen out points with negligible contributions. Default value is True - (enable screening). + Whether to screen out points with negligible contributions. Default value is False. tol_screen : float Screening tolerance for excluding evaluations. Points with values below this tolerance will not be evaluated (they will be set to zero). Internal computed quantities that @@ -808,7 +801,7 @@ def evaluate_general_kinetic_energy_density( alpha, transform=None, deriv_type="general", - screen_basis=True, + screen_basis=False, tol_screen=1e-8, ): r"""Return evaluations of general form of the kinetic energy density at the given points. @@ -846,8 +839,7 @@ def evaluate_general_kinetic_energy_density( and "direct" makes reference to specific implementation of first and second order derivatives for generalized contraction (_eval_first_second_order_deriv_contractions()). screen_basis : bool, optional - Whether to screen out points with negligible contributions. Default value is True - (enable screening). + Whether to screen out points with negligible contributions. Default value is False. tol_screen : float Screening tolerance for excluding evaluations. Points with values below this tolerance will not be evaluated (they will be set to zero). Internal computed quantities that diff --git a/gbasis/evals/eval.py b/gbasis/evals/eval.py index de06e235..6c2fa110 100644 --- a/gbasis/evals/eval.py +++ b/gbasis/evals/eval.py @@ -55,7 +55,7 @@ class Eval(BaseOneIndex): """ @staticmethod - def construct_array_contraction(contractions, points, screen_basis=True, tol_screen=1e-8): + def construct_array_contraction(contractions, points, screen_basis=False, tol_screen=1e-8): r"""Return the evaluations of the given contractions at the given coordinates. Parameters @@ -69,8 +69,7 @@ def construct_array_contraction(contractions, points, screen_basis=True, tol_scr Rows correspond to the points and columns correspond to the :math:`x, y, \text{and} z` components. screen_basis : bool, optional - Whether to screen out points with negligible contributions. Default value is True - (enable screening). + Whether to screen out points with negligible contributions. Default value is False. tol_screen : float Screening tolerance for excluding evaluations. Points with values below this tolerance will not be evaluated (they will be set to zero). Internal computed quantities that @@ -145,7 +144,7 @@ def construct_array_contraction(contractions, points, screen_basis=True, tol_scr return output -def evaluate_basis(basis, points, transform=None, screen_basis=True, tol_screen=1e-8): +def evaluate_basis(basis, points, transform=None, screen_basis=False, tol_screen=1e-8): r"""Evaluate the basis set in the given coordinate system at the given points. Parameters @@ -164,8 +163,7 @@ def evaluate_basis(basis, points, transform=None, screen_basis=True, tol_screen= and index 0 of the array for contractions. Default is no transformation. screen_basis : bool, optional - Whether to screen out points with negligible contributions. Default value is True - (enable screening). + Whether to screen out points with negligible contributions. Default value is False. tol_screen : float Screening tolerance for excluding evaluations. Points with values below this tolerance will not be evaluated (they will be set to zero). Internal computed quantities that diff --git a/gbasis/evals/eval_deriv.py b/gbasis/evals/eval_deriv.py index 5bd8edcd..9d51e252 100644 --- a/gbasis/evals/eval_deriv.py +++ b/gbasis/evals/eval_deriv.py @@ -58,7 +58,7 @@ class EvalDeriv(BaseOneIndex): @staticmethod def construct_array_contraction( - contractions, points, orders, deriv_type="general", screen_basis=True, tol_screen=1e-8 + contractions, points, orders, deriv_type="general", screen_basis=False, tol_screen=1e-8 ): r"""Return the array associated with a set of contracted Cartesian Gaussians. @@ -82,7 +82,7 @@ def construct_array_contraction( contraction (_eval_first_second_order_deriv_contractions()). screen_basis : bool, optional Whether to screen out points that are too far from the contraction center. Default value - is True (enable screening). + False. tol_screen : float Screening tolerance for excluding points. This value, together with the minimum contraction parameters, determines a cutoff distance. Points @@ -177,7 +177,7 @@ def evaluate_deriv_basis( orders, transform=None, deriv_type="general", - screen_basis=True, + screen_basis=False, tol_screen=1e-8, ): r"""Evaluate the derivative of the basis set in the given coordinate system at the given points. @@ -220,7 +220,7 @@ def evaluate_deriv_basis( derivatives for generalized contraction (_eval_first_second_order_deriv_contractions()). screen_basis : bool, optional Whether to screen out points that are too far from the contractions center to contribute to - the result above the given tolerance. Default value is True (enable screening). + the result above the given tolerance. Default value is False. tol_screen : float Screening tolerance for excluding points. This value, together with the minimum contraction parameters, determines a cutoff distance. Points diff --git a/gbasis/integrals/angular_momentum.py b/gbasis/integrals/angular_momentum.py index 7a8413ca..ff01dc72 100644 --- a/gbasis/integrals/angular_momentum.py +++ b/gbasis/integrals/angular_momentum.py @@ -73,7 +73,7 @@ def construct_array_contraction( Contracted Cartesian Gaussians (of the same shell) associated with the second index of the angular momentum integral. screen_basis : bool, optional - A toggle to enable or disable screening. Default value is True to enable screening. + A toggle to enable or disable screening. Default value is False. tol_screen : float, optional The tolerance used for screening overlap integrals. `tol_screen` is combined with the minimum contraction exponents to compute a cutoff which is compared to the distance @@ -218,7 +218,7 @@ def angular_momentum_integral(basis, transform=None, screen_basis=False, tol_scr and index 0 of the array for contractions. Default is no transformation. screen_basis : bool, optional - A toggle to enable or disable screening. Default value is True to enable screening. + A toggle to enable or disable screening. Default value is False. tol_screen : float, optional The tolerance used for screening overlap integrals. `tol_screen` is combined with the minimum contraction exponents to compute a cutoff which is compared to the distance diff --git a/gbasis/integrals/kinetic_energy.py b/gbasis/integrals/kinetic_energy.py index bed75d8b..eb016633 100644 --- a/gbasis/integrals/kinetic_energy.py +++ b/gbasis/integrals/kinetic_energy.py @@ -68,7 +68,7 @@ def construct_array_contraction( Contracted Cartesian Gaussians (of the same shell) associated with the second index of the kinetic energy integral. screen_basis : bool, optional - A toggle to enable or disable screening. Default value is True to enable screening. + A toggle to enable or disable screening. Default value is False. tol_screen : float, optional The tolerance used for screening overlap integrals. `tol_screen` is combined with the minimum contraction exponents to compute a cutoff which is compared to the distance @@ -174,7 +174,7 @@ def kinetic_energy_integral(basis, transform=None, screen_basis=False, tol_scree and index 0 of the array for contractions. Default is no transformation. screen_basis : bool, optional - A toggle to enable or disable screening. Default value is True to enable screening. + A toggle to enable or disable screening. Default value is False. tol_screen : float, optional The tolerance used for screening overlap integrals. `tol_screen` is combined with the minimum contraction exponents to compute a cutoff which is compared to the distance diff --git a/gbasis/integrals/moment.py b/gbasis/integrals/moment.py index 49bdcf80..e2f9a6d3 100644 --- a/gbasis/integrals/moment.py +++ b/gbasis/integrals/moment.py @@ -56,7 +56,7 @@ def construct_array_contraction( contractions_two, moment_coord, moment_orders, - screen_basis=True, + screen_basis=False, tol_screen=1e-8, ): """Return the evaluations of the given contractions at the given coordinates. @@ -76,7 +76,7 @@ def construct_array_contraction( Note that a two dimensional array must be given, even if there is only one set of orders of the moment. screen_basis : bool, optional - A toggle to enable or disable screening. Default value is True to enable screening. + A toggle to enable or disable screening. Default value is False. tol_screen : float, optional The tolerance used for screening overlap integrals. `tol_screen` is combined with the minimum contraction exponents to compute a cutoff which is compared to the distance @@ -217,7 +217,7 @@ def moment_integral( and index 0 of the array for contractions. Default is no transformation. screen_basis : bool, optional - A toggle to enable or disable screening. Default value is True to enable screening. + A toggle to enable or disable screening. Default value is False. tol_screen : float, optional The tolerance used for screening overlap integrals. `tol_screen` is combined with the minimum contraction exponents to compute a cutoff which is compared to the distance diff --git a/gbasis/integrals/momentum.py b/gbasis/integrals/momentum.py index e8496992..b505a0b8 100644 --- a/gbasis/integrals/momentum.py +++ b/gbasis/integrals/momentum.py @@ -70,7 +70,7 @@ def construct_array_contraction( Contracted Cartesian Gaussians (of the same shell) associated with the second index of the momentum energy integral. screen_basis : bool, optional - A toggle to enable or disable screening. Default value is True to enable screening. + A toggle to enable or disable screening. Default value is False. tol_screen : float, optional The tolerance used for screening overlap integrals. `tol_screen` is combined with the minimum contraction exponents to compute a cutoff which is compared to the distance @@ -166,7 +166,7 @@ def momentum_integral(basis, transform=None, screen_basis=False, tol_screen=1e-8 Default is no transformation. Default is no transformation. screen_basis : bool, optional - A toggle to enable or disable screening. Default value is True to enable screening. + A toggle to enable or disable screening. Default value is False. tol_screen : float, optional The tolerance used for screening overlap integrals. `tol_screen` is combined with the minimum contraction exponents to compute a cutoff which is compared to the distance diff --git a/gbasis/integrals/overlap.py b/gbasis/integrals/overlap.py index 0b767170..d5b7b69f 100644 --- a/gbasis/integrals/overlap.py +++ b/gbasis/integrals/overlap.py @@ -66,7 +66,7 @@ def construct_array_contraction( Contracted Cartesian Gaussians (of the same shell) associated with the second index of the overlap. screen_basis : bool, optional - A toggle to enable or disable screening. Default value is True to enable screening. + A toggle to enable or disable screening. Default value is False. tol_screen : float, optional The tolerance used for screening overlap integrals. `tol_screen` is combined with the minimum contraction exponents to compute a cutoff which is compared to the distance @@ -154,7 +154,7 @@ def overlap_integral(basis, transform=None, screen_basis=False, tol_screen=1e-8) and index 0 of the array for contractions. Default is no transformation. screen_basis : bool, optional - A toggle to enable or disable screening. Default value is True to enable screening. + A toggle to enable or disable screening. Default value is False. tol_screen : float, optional The tolerance used for screening overlap integrals. `tol_screen` is combined with the minimum contraction exponents to compute a cutoff which is compared to the distance diff --git a/gbasis/integrals/overlap_asymm.py b/gbasis/integrals/overlap_asymm.py index 7e96d4a5..28f61f9d 100644 --- a/gbasis/integrals/overlap_asymm.py +++ b/gbasis/integrals/overlap_asymm.py @@ -94,7 +94,7 @@ def overlap_integral_asymmetric( and index 0 of the array for contractions. Default is no transformation. screen_basis : bool, optional - A toggle to enable or disable screening. Default value is True to enable screening. + A toggle to enable or disable screening. Default value is False. tol_screen : float, optional The tolerance used for screening overlap integrals. `tol_screen` is combined with the minimum contraction exponents to compute a cutoff which is compared to the distance From 9e75b415c0860632be1069eec13b55703160faa2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Sat, 27 Sep 2025 03:06:04 -0400 Subject: [PATCH 07/14] Add screening to point_charge --- gbasis/integrals/point_charge.py | 134 ++++++++++++++++++++++++++----- 1 file changed, 115 insertions(+), 19 deletions(-) diff --git a/gbasis/integrals/point_charge.py b/gbasis/integrals/point_charge.py index 294f1bd2..e8ce88a7 100644 --- a/gbasis/integrals/point_charge.py +++ b/gbasis/integrals/point_charge.py @@ -1,7 +1,9 @@ """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 +from gbasis.screening import is_two_index_overlap_screened, get_points_mask_for_contraction import numpy as np from scipy.special import hyp1f1 # pylint: disable=E0611 @@ -117,7 +119,13 @@ def boys_func(orders, weighted_dist): @classmethod def construct_array_contraction( - cls, contractions_one, contractions_two, points_coords, points_charge + cls, + contractions_one, + contractions_two, + points_coords, + points_charge, + screen_basis=False, + tol_screen=1e-8, ): r"""Return point charge interaction integral for the given contractions and point charges. @@ -135,6 +143,13 @@ def construct_array_contraction( :math:`x, y, \text{and} z` components. points_charge : np.ndarray(N) Charge of each point charge. + screen_basis : bool, optional + Whether to enable screening of the basis functions. Default value is False. + tol_screen : float, optional + Screening tolerance for excluding evaluations. Integrals with values below this tolerance + will not be evaluated (they will be set to zero). Internal computed quantities that + affect the results below this tolerance will also be ignored to speed up the + evaluation. Default value is 1e-8. Returns ------- @@ -198,7 +213,26 @@ def construct_array_contraction( "`points_charge`." ) - # TODO: Overlap screening + # Overlap screening + if screen_basis and is_two_index_overlap_screened( + contractions_one, contractions_two, tol_screen + ): + # Shape of output array: + # axis 0 : index for segmented contractions of contraction one + # axis 1 : angular momentum vector of contraction one (in the same order as angmoms_a) + # axis 2 : index for segmented contractions of contraction two + # axis 3 : angular momentum vector of contraction two (in the same order as angmoms_b) + # axis 4 : point charge + return np.zeros( + ( + contractions_one.num_seg_cont, + len(contractions_one.angmom_components_cart), + contractions_two.num_seg_cont, + len(contractions_two.angmom_components_cart), + points_coords.shape[0], + ), + dtype=np.float64, + ) coord_a = contractions_one.coord angmom_a = contractions_one.angmom @@ -221,18 +255,66 @@ def construct_array_contraction( coeffs_a, coeffs_b = coeffs_b, coeffs_a ab_swapped = True - integrals = _compute_one_elec_integrals( - points_coords, - cls.boys_func, - coord_a, - angmom_a, - exps_a, - coeffs_a, - coord_b, - angmom_b, - exps_b, - coeffs_b, - ) + if screen_basis: + pts_mask_one = get_points_mask_for_contraction( + contractions_one, points_coords, deriv_order=0, tol_screen=tol_screen + ) + pts_mask_two = get_points_mask_for_contraction( + contractions_two, points_coords, deriv_order=0, tol_screen=tol_screen + ) + pts_mask = np.logical_and(pts_mask_one, pts_mask_two) + + # NOTE: Ordering convention for horizontal recursion of integrals + # axis 0 : a_x (size: angmom_a + 1) + # axis 1 : a_y (size: angmom_a + 1) + # axis 2 : a_z (size: angmom_a + 1) + # axis 3 : b_x (size: angmom_b + 1) + # axis 4 : b_y (size: angmom_b + 1) + # axis 5 : b_z (size: angmom_b + 1) + # axis 6 : point charge (size: N) + # axis 7 : index for segmented contractions of contraction a (size: M_a) + # axis 8 : index for segmented contractions of contraction b (size: M_b) + integrals = np.zeros( + ( + angmom_a + 1, + angmom_a + 1, + angmom_a + 1, + angmom_b + 1, + angmom_b + 1, + angmom_b + 1, + points_coords.shape[0], + coeffs_a.shape[1], + coeffs_b.shape[1], + ), + dtype=np.float64, + ) + if np.any(pts_mask): + integrals[:, :, :, :, :, :, pts_mask, :, :] = _compute_one_elec_integrals( + points_coords[pts_mask], + cls.boys_func, + coord_a, + angmom_a, + exps_a, + coeffs_a, + coord_b, + angmom_b, + exps_b, + coeffs_b, + ) + + else: + integrals = _compute_one_elec_integrals( + points_coords, + cls.boys_func, + coord_a, + angmom_a, + exps_a, + coeffs_a, + coord_b, + angmom_b, + exps_b, + coeffs_b, + ) integrals = np.transpose(integrals, (7, 0, 1, 2, 8, 3, 4, 5, 6)) angmoms_a_x, angmoms_a_y, angmoms_a_z = angmoms_a.T @@ -266,7 +348,9 @@ def construct_array_contraction( return output -def point_charge_integral(basis, points_coords, points_charge, transform=None): +def point_charge_integral( + basis, points_coords, points_charge, transform=None, screen_basis=False, tol_screen=1e-8 +): r"""Return the point-charge interaction integrals of basis set in the given coordinate systems. .. math:: @@ -293,6 +377,13 @@ def point_charge_integral(basis, points_coords, points_charge, transform=None): Transformation is applied to the left, i.e. the sum is over the index 1 of `transform` and index 0 of the array for contractions. Default is no transformation. + screen_basis : bool, optional + Whether to enable screening of the basis functions. Default value is False. + tol_screen : float, optional + Screening tolerance for excluding evaluations. Integrals with values below this tolerance + will not be evaluated (they will be set to zero). Internal computed quantities that + affect the results below this tolerance will also be ignored to speed up the + evaluation. Default value is 1e-8. Returns ------- @@ -305,19 +396,24 @@ def point_charge_integral(basis, points_coords, points_charge, transform=None): """ coord_type = [ct for ct in [shell.coord_type for shell in basis]] + kwargs = {"tol_screen": tol_screen, "screen_basis": screen_basis} if transform is not None: return PointChargeIntegral(basis).construct_array_lincomb( - transform, coord_type, points_coords=points_coords, points_charge=points_charge + transform, + coord_type, + points_coords=points_coords, + points_charge=points_charge, + **kwargs, ) if all(ct == "cartesian" for ct in coord_type): return PointChargeIntegral(basis).construct_array_cartesian( - points_coords=points_coords, points_charge=points_charge + points_coords=points_coords, points_charge=points_charge, **kwargs ) if all(ct == "spherical" for ct in coord_type): return PointChargeIntegral(basis).construct_array_spherical( - points_coords=points_coords, points_charge=points_charge + points_coords=points_coords, points_charge=points_charge, **kwargs ) return PointChargeIntegral(basis).construct_array_mix( - coord_type, points_coords=points_coords, points_charge=points_charge + coord_type, points_coords=points_coords, points_charge=points_charge, **kwargs ) From a012b4609bedc35c50c3e40d4c98a6aee2b3fa5e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Sun, 26 Oct 2025 02:20:11 -0400 Subject: [PATCH 08/14] Add minor grammar fix to documentation --- gbasis/evals/eval_deriv.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gbasis/evals/eval_deriv.py b/gbasis/evals/eval_deriv.py index 9d51e252..d524acfa 100644 --- a/gbasis/evals/eval_deriv.py +++ b/gbasis/evals/eval_deriv.py @@ -82,7 +82,7 @@ def construct_array_contraction( contraction (_eval_first_second_order_deriv_contractions()). screen_basis : bool, optional Whether to screen out points that are too far from the contraction center. Default value - False. + is False. tol_screen : float Screening tolerance for excluding points. This value, together with the minimum contraction parameters, determines a cutoff distance. Points From 047124d16007352f7b72bdb68b73ece3526563be Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Sat, 27 Sep 2025 03:06:50 -0400 Subject: [PATCH 09/14] Add screening tests to point_charge --- tests/test_point_charge.py | 22 +++++++++++++++++----- 1 file changed, 17 insertions(+), 5 deletions(-) diff --git a/tests/test_point_charge.py b/tests/test_point_charge.py index 0b4fd6a7..7d0819ca 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 @@ -30,7 +31,8 @@ def test_boys_func(): assert np.allclose(test, ref.reshape(10, 20, 30)) -def test_construct_array_contraction(): +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) +def test_construct_array_contraction(screen_basis, tol_screen): """Test gbasis.integrals.point_charge.PointChargeIntegral.construct_array_contraction.""" coord_one = np.array([0.5, 1, 1.5]) test_one = GeneralizedContractionShell( @@ -58,7 +60,9 @@ def test_construct_array_contraction(): PointChargeIntegral.construct_array_contraction(test_one, test_two, coord, np.array([0, 1])) assert np.allclose( - PointChargeIntegral.construct_array_contraction(test_one, test_two, coord, charge), + PointChargeIntegral.construct_array_contraction( + test_one, test_two, coord, charge, screen_basis=screen_basis, tol_screen=tol_screen + ), ( 2 * np.pi @@ -71,6 +75,7 @@ def test_construct_array_contraction(): * 3 * (-1) ), + atol=tol_screen, ) test_one = GeneralizedContractionShell( @@ -88,7 +93,9 @@ def test_construct_array_contraction(): for i in range(2) ] assert np.allclose( - PointChargeIntegral.construct_array_contraction(test_one, test_two, coord, charge).ravel(), + PointChargeIntegral.construct_array_contraction( + test_one, test_two, coord, charge, screen_basis=screen_basis, tol_screen=tol_screen + ).ravel(), ( ((coord_wac - coord_one) * v_000_000[0] - (coord_wac - coord[0]) * v_000_000[1]) * (2 * 0.1 / np.pi) ** (3 / 4) @@ -98,6 +105,7 @@ def test_construct_array_contraction(): * 3 * (-1) ), + atol=tol_screen, ) test_one = GeneralizedContractionShell( @@ -107,7 +115,9 @@ def test_construct_array_contraction(): 1, np.array([1.5, 2, 3]), np.array([3.0]), np.array([0.2]), "spherical" ) assert np.allclose( - PointChargeIntegral.construct_array_contraction(test_one, test_two, coord, charge).ravel(), + PointChargeIntegral.construct_array_contraction( + test_one, test_two, coord, charge, screen_basis=screen_basis, tol_screen=tol_screen + ).ravel(), ( ((coord_wac - coord_two) * v_000_000[0] - (coord_wac - coord[0]) * v_000_000[1]) * (2 * 0.1 / np.pi) ** (3 / 4) @@ -117,6 +127,7 @@ def test_construct_array_contraction(): * 3 * (-1) ), + atol=tol_screen, ) @@ -130,7 +141,8 @@ def test_point_charge_cartesian(): points_charge = np.random.rand(5) assert np.allclose( point_charge_obj.construct_array_cartesian( - points_coords=points_coords, points_charge=points_charge + points_coords=points_coords, + points_charge=points_charge, ), point_charge_integral(basis, points_coords=points_coords, points_charge=points_charge), ) From 4f5108c0d647e5c1ae5c7d78097c68330c25b7b0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Sat, 27 Sep 2025 03:20:23 -0400 Subject: [PATCH 10/14] Add screening tests to angular momentum integrals --- gbasis/evals/electrostatic_potential.py | 1 + tests/test_angular_momentum.py | 56 ++++++++++++++++++------- 2 files changed, 42 insertions(+), 15 deletions(-) diff --git a/gbasis/evals/electrostatic_potential.py b/gbasis/evals/electrostatic_potential.py index 4d850975..6cdd4209 100644 --- a/gbasis/evals/electrostatic_potential.py +++ b/gbasis/evals/electrostatic_potential.py @@ -136,6 +136,7 @@ def electrostatic_potential( raise TypeError( "`coord_type` must be a list/tuple of the strings 'spherical' or 'cartesian'." ) + # compute Hartree potential (no screening since it is cumulative for all points) hartree_potential = point_charge_integral( basis, points, -np.ones(points.shape[0]), transform=transform ) diff --git a/tests/test_angular_momentum.py b/tests/test_angular_momentum.py index f915e26c..b75de4c5 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, @@ -11,7 +12,9 @@ from utils import find_datafile -def test_angular_momentum_construct_array_contraction(): +@pytest.mark.parametrize("screen_basis", [True, False]) +@pytest.mark.parametrize("tol_screen", [1e-8]) +def test_angular_momentum_construct_array_contraction(screen_basis, tol_screen): """Test integrals.angular_momentum.angular_momentumIntegral.construct_array_contraction.""" test_one = GeneralizedContractionShell( 1, np.array([0.5, 1, 1.5]), np.array([1.0, 2.0]), np.array([0.1, 0.01]), "spherical" @@ -346,7 +349,9 @@ def test_angular_momentum_construct_array_contraction(): int_pz_dyz = int_pz_dyz.squeeze() int_pz_dzz = int_pz_dzz.squeeze() - test = AngularMomentumIntegral.construct_array_contraction(test_one, test_two) + test = AngularMomentumIntegral.construct_array_contraction( + test_one, test_two, screen_basis=screen_basis, tol_screen=tol_screen + ) assert test.shape == (1, 3, 1, 6, 3) assert np.allclose( test.squeeze(), @@ -358,6 +363,7 @@ def test_angular_momentum_construct_array_contraction(): [int_pz_dxx, int_pz_dxy, int_pz_dxz, int_pz_dyy, int_pz_dyz, int_pz_dzz], ] ), + atol=tol_screen, ) with pytest.raises(TypeError): @@ -366,50 +372,70 @@ def test_angular_momentum_construct_array_contraction(): AngularMomentumIntegral.construct_array_contraction(None, test_two) -def test_angular_momentum_integral_cartesian(): +@pytest.mark.parametrize("screen_basis", [True, False]) +@pytest.mark.parametrize("tol_screen", [1e-8]) +def test_angular_momentum_integral_cartesian(screen_basis, tol_screen): """Test gbasis.integrals.angular_momentum.angular_momentum_integral_cartesian.""" basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "cartesian") angular_momentum_integral_obj = AngularMomentumIntegral(basis) assert np.allclose( - angular_momentum_integral_obj.construct_array_cartesian(), - angular_momentum_integral(basis), + angular_momentum_integral_obj.construct_array_cartesian( + screen_basis=screen_basis, tol_screen=tol_screen + ), + angular_momentum_integral(basis, screen_basis=False), + atol=tol_screen, ) -def test_angular_momentum_integral_spherical(): +@pytest.mark.parametrize("screen_basis", [True, False]) +@pytest.mark.parametrize("tol_screen", [1e-8]) +def test_angular_momentum_integral_spherical(screen_basis, tol_screen): """Test gbasis.integrals.angular_momentum.angular_momentum_integral_spherical.""" basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "spherical") angular_momentum_integral_obj = AngularMomentumIntegral(basis) assert np.allclose( - angular_momentum_integral_obj.construct_array_spherical(), - angular_momentum_integral(basis), + angular_momentum_integral_obj.construct_array_spherical( + screen_basis=screen_basis, tol_screen=tol_screen + ), + angular_momentum_integral(basis, screen_basis=False), + atol=tol_screen, ) -def test_angular_momentum_integral_mix(): +@pytest.mark.parametrize("screen_basis", [True, False]) +@pytest.mark.parametrize("tol_screen", [1e-8]) +def test_angular_momentum_integral_mix(screen_basis, tol_screen): """Test gbasis.integrals.angular_momentum.angular_momentum_integral_mix.""" basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), ["spherical"] * 8) angular_momentum_integral_obj = AngularMomentumIntegral(basis) assert np.allclose( - angular_momentum_integral_obj.construct_array_mix(["spherical"] * 8), - angular_momentum_integral(basis), + angular_momentum_integral_obj.construct_array_mix( + ["spherical"] * 8, screen_basis=screen_basis, tol_screen=tol_screen + ), + angular_momentum_integral(basis, screen_basis=False), + atol=tol_screen, ) -def test_angular_momentum_integral_lincomb(): +@pytest.mark.parametrize("screen_basis", [True, False]) +@pytest.mark.parametrize("tol_screen", [1e-8]) +def test_angular_momentum_integral_lincomb(screen_basis, tol_screen): """Test gbasis.integrals.angular_momentum.angular_momentum_integral_lincomb.""" basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "spherical") angular_momentum_integral_obj = AngularMomentumIntegral(basis) transform = np.random.rand(14, 18) assert np.allclose( - angular_momentum_integral_obj.construct_array_lincomb(transform, ["spherical"]), - angular_momentum_integral(basis, transform), + angular_momentum_integral_obj.construct_array_lincomb( + transform, ["spherical"], screen_basis=screen_basis, tol_screen=tol_screen + ), + angular_momentum_integral(basis, transform, screen_basis=False), + atol=tol_screen, ) @@ -423,7 +449,7 @@ def test_angular_momentum_screening_accuracy(precision): contraction = make_contractions(basis_dict, atsymbols, atcoords, "cartesian") # the screening tolerance needs to be 1e-4 times the desired precision - tol_screen = precision * 1e-4 + tol_screen = precision * 1e-2 angular_momentum = angular_momentum_integral(contraction, tol_screen=tol_screen) angular_momentum_no_screen = angular_momentum_integral(contraction, screen_basis=False) assert np.allclose(angular_momentum, angular_momentum_no_screen, atol=precision) From a2dfc21518d798a53bb664c08046e4215be217ae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Sun, 26 Oct 2025 03:54:14 -0400 Subject: [PATCH 11/14] Simplify density tests parametrization --- tests/test_density.py | 30 ++++++++++-------------------- 1 file changed, 10 insertions(+), 20 deletions(-) diff --git a/tests/test_density.py b/tests/test_density.py index 844d2e0c..8a69f648 100644 --- a/tests/test_density.py +++ b/tests/test_density.py @@ -60,8 +60,7 @@ def test_evaluate_density_using_evaluated_orbs(): evaluate_density_using_evaluated_orbs(density_mat, orb_eval) -@pytest.mark.parametrize("screen_basis", [True, False]) -@pytest.mark.parametrize("tol_screen", [1e-8]) +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) def test_evaluate_density(screen_basis, tol_screen): """Test gbasis.evals.density.evaluate_density.""" basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) @@ -83,8 +82,7 @@ def test_evaluate_density(screen_basis, tol_screen): ) -@pytest.mark.parametrize("screen_basis", [True, False]) -@pytest.mark.parametrize("tol_screen", [1e-8]) +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) @pytest.mark.parametrize("deriv_type", ["general", "direct"]) def test_evaluate_deriv_density(screen_basis, tol_screen, deriv_type): """Test gbasis.evals.density.evaluate_deriv_density.""" @@ -272,8 +270,7 @@ def test_evaluate_deriv_density(screen_basis, tol_screen, deriv_type): ) -@pytest.mark.parametrize("screen_basis", [True, False]) -@pytest.mark.parametrize("tol_screen", [1e-8]) +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) @pytest.mark.parametrize("deriv_type", ["general", "direct"]) def test_evaluate_density_gradient(screen_basis, tol_screen, deriv_type): """Test gbasis.evals.density.evaluate_density_gradient.""" @@ -338,8 +335,7 @@ def test_evaluate_density_gradient(screen_basis, tol_screen, deriv_type): ) -@pytest.mark.parametrize("screen_basis", [True, False]) -@pytest.mark.parametrize("tol_screen", [1e-8]) +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) def test_evaluate_density_horton(screen_basis, tol_screen): """Test gbasis.evals.density.evaluate_density against result from HORTON. @@ -370,8 +366,7 @@ def test_evaluate_density_horton(screen_basis, tol_screen): assert np.allclose(dens, horton_density, atol=tol_screen) -@pytest.mark.parametrize("screen_basis", [True, False]) -@pytest.mark.parametrize("tol_screen", [1e-8]) +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) @pytest.mark.parametrize("deriv_type", ["general", "direct"]) def test_evaluate_density_gradient_horton(screen_basis, tol_screen, deriv_type): """Test gbasis.evals.density.evaluate_density_gradient against result from HORTON. @@ -405,8 +400,7 @@ def test_evaluate_density_gradient_horton(screen_basis, tol_screen, deriv_type): ) -@pytest.mark.parametrize("screen_basis", [True, False]) -@pytest.mark.parametrize("tol_screen", [1e-8]) +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) @pytest.mark.parametrize("deriv_type", ["general", "direct"]) def test_evaluate_hessian_deriv_horton(screen_basis, tol_screen, deriv_type): """Test gbasis.evals.density.evaluate_density_hessian against result from HORTON. @@ -447,8 +441,7 @@ def test_evaluate_hessian_deriv_horton(screen_basis, tol_screen, deriv_type): ) -@pytest.mark.parametrize("screen_basis", [True, False]) -@pytest.mark.parametrize("tol_screen", [1e-8]) +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) @pytest.mark.parametrize("deriv_type", ["general", "direct"]) def test_evaluate_laplacian_deriv_horton(screen_basis, tol_screen, deriv_type): """Test gbasis.evals.density.evaluate_density_laplacian against result from HORTON. @@ -483,8 +476,7 @@ def test_evaluate_laplacian_deriv_horton(screen_basis, tol_screen, deriv_type): ) -@pytest.mark.parametrize("screen_basis", [True, False]) -@pytest.mark.parametrize("tol_screen", [1e-8]) +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) @pytest.mark.parametrize("deriv_type", ["general", "direct"]) def test_evaluate_posdef_kinetic_energy_density(screen_basis, tol_screen, deriv_type): """Test evaluate_posdef_kinetic_energy_density against results from HORTON. @@ -519,8 +511,7 @@ def test_evaluate_posdef_kinetic_energy_density(screen_basis, tol_screen, deriv_ assert np.allclose(dens, horton_density_kinetic_density, atol=tol_screen) -@pytest.mark.parametrize("screen_basis", [True, False]) -@pytest.mark.parametrize("tol_screen", [1e-8]) +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) @pytest.mark.parametrize("deriv_type", ["general", "direct"]) def test_evaluate_general_kinetic_energy_density_horton(screen_basis, tol_screen, deriv_type): """Test evaluate_general_kinetic_energy_density against results from HORTON. @@ -560,8 +551,7 @@ def test_evaluate_general_kinetic_energy_density_horton(screen_basis, tol_screen ) -@pytest.mark.parametrize("screen_basis", [True, False]) -@pytest.mark.parametrize("tol_screen", [1e-8]) +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) @pytest.mark.parametrize("deriv_type", ["general", "direct"]) def test_evaluate_general_kinetic_energy_density(screen_basis, tol_screen, deriv_type): """Test density.evaluate_general_kinetic_energy_density.""" From 2567ff1efb9e93e5855f3be8cf68d40df2939ede Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Sun, 28 Sep 2025 04:45:09 -0400 Subject: [PATCH 12/14] Add screening to nuclear_electron_attraction --- .../integrals/nuclear_electron_attraction.py | 21 +++++++++++++++++-- 1 file changed, 19 insertions(+), 2 deletions(-) diff --git a/gbasis/integrals/nuclear_electron_attraction.py b/gbasis/integrals/nuclear_electron_attraction.py index 00e2241e..4efd5180 100644 --- a/gbasis/integrals/nuclear_electron_attraction.py +++ b/gbasis/integrals/nuclear_electron_attraction.py @@ -1,9 +1,12 @@ """Module for computing the nuclear electron attraction.""" + from gbasis.integrals.point_charge import point_charge_integral import numpy as np -def nuclear_electron_attraction_integral(basis, nuclear_coords, nuclear_charges, transform=None): +def nuclear_electron_attraction_integral( + basis, nuclear_coords, nuclear_charges, transform=None, screen_basis=False, tol_screen=1e-8 +): """Return the nuclear electron attraction integrals of the basis set in the Cartesian form. .. math:: @@ -27,6 +30,13 @@ def nuclear_electron_attraction_integral(basis, nuclear_coords, nuclear_charges, Transformation is applied to the left, i.e. the sum is over the index 1 of `transform` and index 0 of the array for contractions. Default is no transformation. + screen_basis : bool, optional + Whether to enable screening of the basis functions. Default value is False. + tol_screen : float, optional + Screening tolerance for excluding evaluations. Integrals with values below this tolerance + will not be evaluated (they will be set to zero). Internal computed quantities that + affect the results below this tolerance will also be ignored to speed up the + evaluation. Default value is 1e-8. Returns ------- @@ -37,6 +47,13 @@ def nuclear_electron_attraction_integral(basis, nuclear_coords, nuclear_charges, """ return np.sum( - point_charge_integral(basis, nuclear_coords, nuclear_charges, transform=transform), + point_charge_integral( + basis, + nuclear_coords, + nuclear_charges, + transform=transform, + screen_basis=screen_basis, + tol_screen=tol_screen, + ), axis=2, ) From 2350cfce1302342f4711179ce55ea8a04f2783d4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Sun, 26 Oct 2025 04:31:16 -0400 Subject: [PATCH 13/14] Test screening for nuclear_electron_attraction --- tests/test_nuclear_electron_attraction.py | 54 ++++++++++++++++------- 1 file changed, 38 insertions(+), 16 deletions(-) diff --git a/tests/test_nuclear_electron_attraction.py b/tests/test_nuclear_electron_attraction.py index 0e67871b..1840ed3e 100644 --- a/tests/test_nuclear_electron_attraction.py +++ b/tests/test_nuclear_electron_attraction.py @@ -1,17 +1,21 @@ """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 import numpy as np from utils import find_datafile, HortonContractions +import pytest -def test_nuclear_electron_attraction_horton_anorcc_hhe(): +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) +def test_nuclear_electron_attraction_horton_anorcc_hhe(screen_basis, tol_screen): """Test nuclear_electron_attraciton.nuclear_electron_attraction_cartesian with HORTON. The test case is diatomic with H and He separated by 0.8 angstroms with basis set ANO-RCC. """ + kwargs = {"screen_basis": screen_basis, "tol_screen": tol_screen} basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem")) # NOTE: used HORTON's conversion factor for angstroms to bohr coords = np.array([[0, 0, 0], [0.8 * 1.0 / 0.5291772083, 0, 0]]) @@ -20,17 +24,20 @@ def test_nuclear_electron_attraction_horton_anorcc_hhe(): horton_nucattract = np.load(find_datafile("data_horton_hhe_cart_nucattract.npy")) assert np.allclose( - nuclear_electron_attraction_integral(basis, coords, np.array([1, 2])), + nuclear_electron_attraction_integral(basis, coords, np.array([1, 2]), **kwargs), horton_nucattract, + atol=tol_screen, ) -def test_nuclear_electron_attraction_horton_anorcc_bec(): +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) +def test_nuclear_electron_attraction_horton_anorcc_bec(screen_basis, tol_screen): """Test nuclear_electron_attraciton.nuclear_electron_attraction_cartesian with HORTON. The test case is diatomic with B and C separated by 1.0 angstroms with basis set ANO-RCC. """ + kwargs = {"screen_basis": screen_basis, "tol_screen": tol_screen} basis_dict = parse_nwchem(find_datafile("data_anorcc.nwchem")) # NOTE: used HORTON's conversion factor for angstroms to bohr coords = np.array([[0, 0, 0], [1.0 * 1.0 / 0.5291772083, 0, 0]]) @@ -39,65 +46,80 @@ def test_nuclear_electron_attraction_horton_anorcc_bec(): horton_nucattract = np.load(find_datafile("data_horton_bec_cart_nucattract.npy")) assert np.allclose( - nuclear_electron_attraction_integral(basis, coords, np.array([4, 6])), + nuclear_electron_attraction_integral(basis, coords, np.array([4, 6]), **kwargs), horton_nucattract, + atol=tol_screen, ) -def test_nuclear_electron_attraction_cartesian(): +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) +def test_nuclear_electron_attraction_cartesian(screen_basis, tol_screen): """Test gbasis.integrals.nuclear_electron_attraction.nuclear_electron_attraction_cartesian.""" + kwargs = {"screen_basis": screen_basis, "tol_screen": tol_screen} basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "cartesian") nuclear_coords = np.random.rand(5, 3) nuclear_charges = np.random.rand(5) - ref = point_charge_integral(basis, nuclear_coords, nuclear_charges) + ref = point_charge_integral(basis, nuclear_coords, nuclear_charges, **kwargs) assert np.allclose( ref[:, :, 0] + ref[:, :, 1] + ref[:, :, 2] + ref[:, :, 3] + ref[:, :, 4], - nuclear_electron_attraction_integral(basis, nuclear_coords, nuclear_charges), + nuclear_electron_attraction_integral(basis, nuclear_coords, nuclear_charges, **kwargs), + atol=1e-10, ) -def test_nuclear_electron_attraction_spherical(): +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) +def test_nuclear_electron_attraction_spherical(screen_basis, tol_screen): """Test gbasis.integrals.nuclear_electron_attraction.nuclear_electron_attraction_spherical.""" + kwargs = {"screen_basis": screen_basis, "tol_screen": tol_screen} basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "spherical") nuclear_coords = np.random.rand(5, 3) nuclear_charges = np.random.rand(5) - ref = point_charge_integral(basis, nuclear_coords, nuclear_charges) + ref = point_charge_integral(basis, nuclear_coords, nuclear_charges, **kwargs) assert np.allclose( ref[:, :, 0] + ref[:, :, 1] + ref[:, :, 2] + ref[:, :, 3] + ref[:, :, 4], - nuclear_electron_attraction_integral(basis, nuclear_coords, nuclear_charges), + nuclear_electron_attraction_integral(basis, nuclear_coords, nuclear_charges, **kwargs), + atol=1e-10, ) -def test_nuclear_electron_attraction_mix(): +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) +def test_nuclear_electron_attraction_mix(screen_basis, tol_screen): """Test gbasis.integrals.nuclear_electron_attraction.nuclear_electron_attraction_mix.""" + kwargs = {"screen_basis": screen_basis, "tol_screen": tol_screen} basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), ["spherical"] * 8) nuclear_coords = np.random.rand(5, 3) nuclear_charges = np.random.rand(5) - ref = point_charge_integral(basis, nuclear_coords, nuclear_charges) + ref = point_charge_integral(basis, nuclear_coords, nuclear_charges, **kwargs) assert np.allclose( ref[:, :, 0] + ref[:, :, 1] + ref[:, :, 2] + ref[:, :, 3] + ref[:, :, 4], - nuclear_electron_attraction_integral(basis, nuclear_coords, nuclear_charges), + nuclear_electron_attraction_integral(basis, nuclear_coords, nuclear_charges, **kwargs), + atol=1e-10, ) -def test_nuclear_electron_attraction_lincomb(): +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) +def test_nuclear_electron_attraction_lincomb(screen_basis, tol_screen): """Test gbasis.integrals.nuclear_electron_attraction.nuclear_electron_attraction_lincomb.""" + kwargs = {"screen_basis": screen_basis, "tol_screen": tol_screen} basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "spherical") nuclear_coords = np.random.rand(5, 3) nuclear_charges = np.random.rand(5) transform = np.random.rand(14, 18) - ref = point_charge_integral(basis, nuclear_coords, nuclear_charges, transform=transform) + ref = point_charge_integral( + basis, nuclear_coords, nuclear_charges, transform=transform, **kwargs + ) assert np.allclose( ref[:, :, 0] + ref[:, :, 1] + ref[:, :, 2] + ref[:, :, 3] + ref[:, :, 4], nuclear_electron_attraction_integral( - basis, nuclear_coords, nuclear_charges, transform=transform + basis, nuclear_coords, nuclear_charges, transform=transform, **kwargs ), + atol=1e-10, ) From 508f16b638acc463b735d17efc162352d26b64d4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marco=20Mart=C3=ADnez=20Gonz=C3=A1lez?= Date: Sun, 26 Oct 2025 04:59:51 -0400 Subject: [PATCH 14/14] Fix inconsitency in a comment, simplify tests parametrizations for angular_momentum --- tests/test_angular_momentum.py | 17 ++++++----------- 1 file changed, 6 insertions(+), 11 deletions(-) diff --git a/tests/test_angular_momentum.py b/tests/test_angular_momentum.py index b75de4c5..b23fd688 100644 --- a/tests/test_angular_momentum.py +++ b/tests/test_angular_momentum.py @@ -12,8 +12,7 @@ from utils import find_datafile -@pytest.mark.parametrize("screen_basis", [True, False]) -@pytest.mark.parametrize("tol_screen", [1e-8]) +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) def test_angular_momentum_construct_array_contraction(screen_basis, tol_screen): """Test integrals.angular_momentum.angular_momentumIntegral.construct_array_contraction.""" test_one = GeneralizedContractionShell( @@ -372,8 +371,7 @@ def test_angular_momentum_construct_array_contraction(screen_basis, tol_screen): AngularMomentumIntegral.construct_array_contraction(None, test_two) -@pytest.mark.parametrize("screen_basis", [True, False]) -@pytest.mark.parametrize("tol_screen", [1e-8]) +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) def test_angular_momentum_integral_cartesian(screen_basis, tol_screen): """Test gbasis.integrals.angular_momentum.angular_momentum_integral_cartesian.""" basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) @@ -388,8 +386,7 @@ def test_angular_momentum_integral_cartesian(screen_basis, tol_screen): ) -@pytest.mark.parametrize("screen_basis", [True, False]) -@pytest.mark.parametrize("tol_screen", [1e-8]) +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) def test_angular_momentum_integral_spherical(screen_basis, tol_screen): """Test gbasis.integrals.angular_momentum.angular_momentum_integral_spherical.""" basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) @@ -405,8 +402,7 @@ def test_angular_momentum_integral_spherical(screen_basis, tol_screen): ) -@pytest.mark.parametrize("screen_basis", [True, False]) -@pytest.mark.parametrize("tol_screen", [1e-8]) +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) def test_angular_momentum_integral_mix(screen_basis, tol_screen): """Test gbasis.integrals.angular_momentum.angular_momentum_integral_mix.""" basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) @@ -422,8 +418,7 @@ def test_angular_momentum_integral_mix(screen_basis, tol_screen): ) -@pytest.mark.parametrize("screen_basis", [True, False]) -@pytest.mark.parametrize("tol_screen", [1e-8]) +@pytest.mark.parametrize("screen_basis, tol_screen", [(True, 1e-8), (False, 1e-12)]) def test_angular_momentum_integral_lincomb(screen_basis, tol_screen): """Test gbasis.integrals.angular_momentum.angular_momentum_integral_lincomb.""" basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem")) @@ -448,7 +443,7 @@ def test_angular_momentum_screening_accuracy(precision): atcoords = np.array([[0, 0, 0], [1, 1, 1], [2, 2, 2]]) contraction = make_contractions(basis_dict, atsymbols, atcoords, "cartesian") - # the screening tolerance needs to be 1e-4 times the desired precision + # the screening tolerance needs to be 1e-2 times the desired precision tol_screen = precision * 1e-2 angular_momentum = angular_momentum_integral(contraction, tol_screen=tol_screen) angular_momentum_no_screen = angular_momentum_integral(contraction, screen_basis=False)