diff --git a/gbasis/evals/density.py b/gbasis/evals/density.py
index 82341804..820ed05a 100644
--- a/gbasis/evals/density.py
+++ b/gbasis/evals/density.py
@@ -62,7 +62,15 @@ def evaluate_density_using_evaluated_orbs(one_density_matrix, orb_eval):
     return np.sum(density, axis=0)
 
 
-def evaluate_density(one_density_matrix, basis, points, transform=None, threshold=1.0e-8):
+def evaluate_density(
+    one_density_matrix,
+    basis,
+    points,
+    transform=None,
+    threshold=1.0e-8,
+    screen_basis=True,
+    tol_screen=1e-8,
+):
     r"""Return the density of the given basis set at the given points.
 
     .. math::
@@ -94,6 +102,13 @@ def evaluate_density(one_density_matrix, basis, points, transform=None, threshol
     threshold : float, optional
         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
+        A toggle to enable or disable screening. Default value is True (enable screening).
+    tol_screen : float, optional
+        The tolerance used for screening overlap integrals. `tol_screen` is combined with the
+        minimum contraction exponents to compute a cutoff radius which is compared to the distance
+        between the points and the contraction centers to decide whether the basis function
+        should be evaluated or set to zero at that point.
 
     Returns
     -------
@@ -101,7 +116,9 @@ def evaluate_density(one_density_matrix, basis, points, transform=None, threshol
         Density evaluated at `N` grid points.
 
     """
-    orb_eval = evaluate_basis(basis, points, transform=transform)
+    orb_eval = evaluate_basis(
+        basis, points, transform=transform, screen_basis=screen_basis, tol_screen=tol_screen
+    )
     output = evaluate_density_using_evaluated_orbs(one_density_matrix, orb_eval)
     # Fix #117: check magnitude of small negative density values, then use clip to remove them
     min_output = np.min(output)
@@ -489,7 +506,7 @@ def evaluate_density_hessian(
     r"""Return the Hessian of the density evaluated at the given points.
 
     .. math::
-        
+
         H[\rho(\mathbf{r})]
         =
         \begin{bmatrix}
diff --git a/gbasis/evals/eval.py b/gbasis/evals/eval.py
index bd14b8d6..2e3817a0 100644
--- a/gbasis/evals/eval.py
+++ b/gbasis/evals/eval.py
@@ -2,6 +2,7 @@
 from gbasis.base_one import BaseOneIndex
 from gbasis.contractions import GeneralizedContractionShell
 from gbasis.evals._deriv import _eval_deriv_contractions
+from gbasis.screening import get_points_mask_for_contraction
 import numpy as np
 
 
@@ -53,7 +54,7 @@ class Eval(BaseOneIndex):
     """
 
     @staticmethod
-    def construct_array_contraction(contractions, points):
+    def construct_array_contraction(contractions, points, screen_basis=True, tol_screen=1e-8):
         r"""Return the evaluations of the given contractions at the given coordinates.
 
         Parameters
@@ -105,13 +106,29 @@ def construct_array_contraction(contractions, points):
         angmom_comps = contractions.angmom_components_cart
         center = contractions.coord
         norm_prim_cart = contractions.norm_prim_cart
-        output = _eval_deriv_contractions(
-            points, np.zeros(3), center, angmom_comps, alphas, prim_coeffs, norm_prim_cart
+
+        # if screening is not requested, evaluate all points
+        if not screen_basis:
+            return _eval_deriv_contractions(
+                points, np.zeros(3), center, angmom_comps, alphas, prim_coeffs, norm_prim_cart
+            )
+
+        # default case, screen points that are too far from the contraction center
+        points_mask = get_points_mask_for_contraction(contractions, points, tol_screen=tol_screen)
+        # reconstruct the array with correct shape
+        L = angmom_comps.shape[0]
+        M = prim_coeffs.shape[1]
+        N = points.shape[0]
+        output = np.zeros((M, L, N), dtype=np.float64)
+
+        # fill non-screened points in the output array
+        output[:,:, points_mask] = _eval_deriv_contractions(
+            points[points_mask], np.zeros(3), center, angmom_comps, alphas, prim_coeffs, norm_prim_cart
         )
         return output
 
 
-def evaluate_basis(basis, points, transform=None):
+def evaluate_basis(basis, points, transform=None, screen_basis=True, tol_screen=1e-8):
     r"""Evaluate the basis set in the given coordinate system at the given points.
 
     Parameters
@@ -129,6 +146,13 @@ def evaluate_basis(basis, points, 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
+        A toggle to enable or disable screening. Default value is True (enable screening).
+    tol_screen : float, optional
+        The tolerance used for screening overlap integrals. `tol_screen` is combined with the
+        minimum contraction exponents to compute a cutoff radius which is compared to the distance
+        between the points and the contraction centers to decide whether the basis function
+        should be evaluated or set to zero at that point.
 
     Returns
     -------
@@ -141,11 +165,12 @@ def evaluate_basis(basis, points, 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 Eval(basis).construct_array_lincomb(transform, coord_type, points=points)
+        return Eval(basis).construct_array_lincomb(transform, coord_type, points=points, **kwargs)
     if all(ct == "cartesian" for ct in coord_type):
-        return Eval(basis).construct_array_cartesian(points=points)
+        return Eval(basis).construct_array_cartesian(points=points, **kwargs)
     if all(ct == "spherical" for ct in coord_type):
-        return Eval(basis).construct_array_spherical(points=points)
-    return Eval(basis).construct_array_mix(coord_type, points=points)
+        return Eval(basis).construct_array_spherical(points=points, **kwargs)
+    return Eval(basis).construct_array_mix(coord_type, points=points, **kwargs)
diff --git a/gbasis/screening.py b/gbasis/screening.py
index 8a214705..1ace911f 100644
--- a/gbasis/screening.py
+++ b/gbasis/screening.py
@@ -1,6 +1,8 @@
 """1 and 2-index screening functions"""
 
 import numpy as np
+from scipy.special import lambertw
+from gbasis.utils import factorial2
 
 
 def is_two_index_overlap_screened(contractions_one, contractions_two, tol_screen):
@@ -42,3 +44,98 @@ def is_two_index_overlap_screened(contractions_one, contractions_two, tol_screen
     cutoff = np.sqrt(-(alpha_a + alpha_b) / (alpha_a * alpha_b) * np.log(tol_screen))
     # integrals are screened if centers are further apart than the cutoff
     return np.linalg.norm(r_12) > cutoff
+
+
+def get_points_mask_for_contraction(contractions, points, tol_screen):
+    r"""Return a boolean mask indicating which points should be screened for a contraction shell
+
+    A point is considered screened if it lies farther from the contraction center than a cutoff
+    radius computed from the contraction parameters and the screening tolerance
+    (see :func:`compute_primitive_cutoff_radius`).
+
+    Parameters
+    ----------
+    contractions : GeneralizedContractionShell
+        Contracted Cartesian Gaussians (of the same shell) for which the mask is computed.
+    points : np.ndarray of shape (N, 3)
+        Cartesian coordinates of the points in space (in atomic units) where the
+        basis functions are evaluated. Rows correspond to points; columns
+        correspond to the :math:`x`, :math:`y`, and :math:`z` components.
+    tol_screen : float
+        Screening tolerance for excluding points. This value, together with the
+        minimum contraction parameters, determines a cutoff distance. Points
+        farther than the cutoff are excluded from contraction evaluation.
+
+    Returns
+    -------
+    np.ndarray of bool, shape (N,)
+        Boolean mask where `False` marks points to be screened out.
+    """
+    angm = contractions.angmom
+    # reshape exponents for broadcasting
+    exps = contractions.exps[:, np.newaxis]  # shape (K, 1)
+    # use absolute value (indicating magnitude) of primitive contributions
+    coeffs = np.abs(contractions.coeffs)  # shape (K, M)
+
+    # compute cutoff radius for all primitives in all contractions
+    r_cuts = compute_primitive_cutoff_radius(coeffs, exps, angm, tol_screen)
+
+    # pick the maximum radius over all primitives and contractions
+    cutoff_radius = np.max(r_cuts)
+
+    # screen out points beyond the cutoff radius
+    points_r = np.linalg.norm(points - contractions.coord, axis=1)
+    return points_r <= cutoff_radius
+
+
+def compute_primitive_cutoff_radius(c, alpha, angm, tol_screen):
+    r"""Compute the cutoff radius for a primitive Gaussian.
+
+    The cutoff radius is the maximum distance from the center of the primitive Gaussian at which the
+    radial density remains above a given tolerance :math:`\epsilon`. This radius is computed by
+    solving the equation:
+
+    .. math::
+
+        r^{2\ell} \chi(r)^2 = \epsilon^2
+
+    where :math:`\chi(r)` is the radial part of the primitive Gaussian, defined as:
+
+    .. math::
+
+        \chi(r) = c \, n \, e^{-\alpha r^2}
+
+    Here, :math:`c` is the coefficient of the primitive Gaussian, :math:`\alpha` is its exponent,
+    :math:`\ell` is the angular momentum quantum number, and :math:`n` is the normalization factor
+    given by:
+
+    .. math::
+
+        n = \left( \frac{2 \alpha}{\pi} \right)^{\frac{1}{4}}
+            \frac{(4 \alpha)^{\frac{\ell}{2}}}{\sqrt{(2\ell + 1)!!}}
+
+    Parameters
+    ----------
+    c : float
+        Coefficient of the primitive Gaussian.
+    alpha : float
+        Exponent :math:`\alpha` of the primitive Gaussian.
+    angm : int
+        Angular momentum quantum number :math:`\ell` of the primitive Gaussian.
+    dens_tol : float
+        Radial density tolerance :math:`\epsilon` for computing the cutoff radius.
+
+    Returns
+    -------
+    float
+        The cutoff radius for the primitive Gaussian.
+    """
+    # Compute normalization factor n for the primitive Gaussian
+    n = (2 * alpha / np.pi) ** 0.25 * (4 * alpha) ** (angm / 2) / np.sqrt(factorial2(2 * angm + 1))
+    # special case for angular momentum 0. Solution found using logarithm
+    if angm == 0:
+        return np.sqrt(-np.log(tol_screen / (c * n)) / alpha)
+    # general case for angular momentum > 0. Solution found in terms of the Lambert W function
+    # W_{-1} branch corresponds to the outermost solution
+    lambert_input_value = -2 * alpha * (tol_screen / (c * n)) ** (2 / angm) / angm
+    return np.sqrt(-(angm / (2 * alpha)) * lambertw(lambert_input_value, k=-1).real)
diff --git a/tests/test_density.py b/tests/test_density.py
index 0b12eeda..15d07a81 100644
--- a/tests/test_density.py
+++ b/tests/test_density.py
@@ -68,7 +68,7 @@ def test_evaluate_density():
     density += density.T
     points = np.random.rand(10, 3)
 
-    evaluate_orbs = evaluate_basis(basis, points, transform)
+    evaluate_orbs = evaluate_basis(basis, points, transform, screen_basis=False)
     dens = evaluate_density(density, basis, points, transform)
     assert np.all(dens >= 0.0)
     assert np.allclose(dens, np.einsum("ij,ik,jk->k", density, evaluate_orbs, evaluate_orbs))
@@ -89,12 +89,12 @@ def test_evaluate_deriv_density():
             "ij,ik,jk->k",
             density,
             evaluate_deriv_basis(basis, points, np.array([1, 0, 0]), transform),
-            evaluate_basis(basis, points, transform),
+            evaluate_basis(basis, points, transform, screen_basis=False),
         )
         + np.einsum(
             "ij,ik,jk->k",
             density,
-            evaluate_basis(basis, points, transform),
+            evaluate_basis(basis, points, transform, screen_basis=False),
             evaluate_deriv_basis(basis, points, np.array([1, 0, 0]), transform),
         ),
     )
@@ -105,12 +105,12 @@ def test_evaluate_deriv_density():
             "ij,ik,jk->k",
             density,
             evaluate_deriv_basis(basis, points, np.array([0, 1, 0]), transform),
-            evaluate_basis(basis, points, transform),
+            evaluate_basis(basis, points, transform, screen_basis=False),
         )
         + np.einsum(
             "ij,ik,jk->k",
             density,
-            evaluate_basis(basis, points, transform),
+            evaluate_basis(basis, points, transform, screen_basis=False),
             evaluate_deriv_basis(basis, points, np.array([0, 1, 0]), transform),
         ),
     )
@@ -121,12 +121,12 @@ def test_evaluate_deriv_density():
             "ij,ik,jk->k",
             density,
             evaluate_deriv_basis(basis, points, np.array([0, 0, 1]), transform),
-            evaluate_basis(basis, points, transform),
+            evaluate_basis(basis, points, transform, screen_basis=False),
         )
         + np.einsum(
             "ij,ik,jk->k",
             density,
-            evaluate_basis(basis, points, transform),
+            evaluate_basis(basis, points, transform, screen_basis=False),
             evaluate_deriv_basis(basis, points, np.array([0, 0, 1]), transform),
         ),
     )
@@ -235,36 +235,36 @@ def test_evaluate_density_gradient():
                     "ij,ik,jk->k",
                     density,
                     evaluate_deriv_basis(basis, points, np.array([1, 0, 0]), transform),
-                    evaluate_basis(basis, points, transform),
+                    evaluate_basis(basis, points, transform, screen_basis=False),
                 )
                 + np.einsum(
                     "ij,ik,jk->k",
                     density,
-                    evaluate_basis(basis, points, transform),
+                    evaluate_basis(basis, points, transform, screen_basis=False),
                     evaluate_deriv_basis(basis, points, np.array([1, 0, 0]), transform),
                 ),
                 np.einsum(
                     "ij,ik,jk->k",
                     density,
                     evaluate_deriv_basis(basis, points, np.array([0, 1, 0]), transform),
-                    evaluate_basis(basis, points, transform),
+                    evaluate_basis(basis, points, transform, screen_basis=False),
                 )
                 + np.einsum(
                     "ij,ik,jk->k",
                     density,
-                    evaluate_basis(basis, points, transform),
+                    evaluate_basis(basis, points, transform, screen_basis=False),
                     evaluate_deriv_basis(basis, points, np.array([0, 1, 0]), transform),
                 ),
                 np.einsum(
                     "ij,ik,jk->k",
                     density,
                     evaluate_deriv_basis(basis, points, np.array([0, 0, 1]), transform),
-                    evaluate_basis(basis, points, transform),
+                    evaluate_basis(basis, points, transform, screen_basis=False),
                 )
                 + np.einsum(
                     "ij,ik,jk->k",
                     density,
-                    evaluate_basis(basis, points, transform),
+                    evaluate_basis(basis, points, transform, screen_basis=False),
                     evaluate_deriv_basis(basis, points, np.array([0, 0, 1]), transform),
                 ),
             ]
diff --git a/tests/test_density_direct.py b/tests/test_density_direct.py
index a95397b0..618763a5 100644
--- a/tests/test_density_direct.py
+++ b/tests/test_density_direct.py
@@ -69,7 +69,7 @@ def test_evaluate_density():
     density += density.T
     points = np.random.rand(10, 3)
 
-    evaluate_orbs = evaluate_basis(basis, points, transform)
+    evaluate_orbs = evaluate_basis(basis, points, transform, screen_basis=False)
     assert np.allclose(
         evaluate_density(density, basis, points, transform),
         np.einsum("ij,ik,jk->k", density, evaluate_orbs, evaluate_orbs),
@@ -93,12 +93,12 @@ def test_evaluate_deriv_density():
             "ij,ik,jk->k",
             density,
             evaluate_deriv_basis(basis, points, np.array([1, 0, 0]), transform),
-            evaluate_basis(basis, points, transform),
+            evaluate_basis(basis, points, transform, screen_basis=False),
         )
         + np.einsum(
             "ij,ik,jk->k",
             density,
-            evaluate_basis(basis, points, transform),
+            evaluate_basis(basis, points, transform, screen_basis=False),
             evaluate_deriv_basis(basis, points, np.array([1, 0, 0]), transform),
         ),
     )
@@ -111,12 +111,12 @@ def test_evaluate_deriv_density():
             "ij,ik,jk->k",
             density,
             evaluate_deriv_basis(basis, points, np.array([0, 1, 0]), transform),
-            evaluate_basis(basis, points, transform),
+            evaluate_basis(basis, points, transform, screen_basis=False),
         )
         + np.einsum(
             "ij,ik,jk->k",
             density,
-            evaluate_basis(basis, points, transform),
+            evaluate_basis(basis, points, transform, screen_basis=False),
             evaluate_deriv_basis(basis, points, np.array([0, 1, 0]), transform),
         ),
     )
@@ -129,12 +129,12 @@ def test_evaluate_deriv_density():
             "ij,ik,jk->k",
             density,
             evaluate_deriv_basis(basis, points, np.array([0, 0, 1]), transform),
-            evaluate_basis(basis, points, transform),
+            evaluate_basis(basis, points, transform, screen_basis=False),
         )
         + np.einsum(
             "ij,ik,jk->k",
             density,
-            evaluate_basis(basis, points, transform),
+            evaluate_basis(basis, points, transform, screen_basis=False),
             evaluate_deriv_basis(basis, points, np.array([0, 0, 1]), transform),
         ),
     )
@@ -245,36 +245,36 @@ def test_evaluate_density_gradient():
                     "ij,ik,jk->k",
                     density,
                     evaluate_deriv_basis(basis, points, np.array([1, 0, 0]), transform),
-                    evaluate_basis(basis, points, transform),
+                    evaluate_basis(basis, points, transform, screen_basis=False),
                 )
                 + np.einsum(
                     "ij,ik,jk->k",
                     density,
-                    evaluate_basis(basis, points, transform),
+                    evaluate_basis(basis, points, transform, screen_basis=False),
                     evaluate_deriv_basis(basis, points, np.array([1, 0, 0]), transform),
                 ),
                 np.einsum(
                     "ij,ik,jk->k",
                     density,
                     evaluate_deriv_basis(basis, points, np.array([0, 1, 0]), transform),
-                    evaluate_basis(basis, points, transform),
+                    evaluate_basis(basis, points, transform, screen_basis=False),
                 )
                 + np.einsum(
                     "ij,ik,jk->k",
                     density,
-                    evaluate_basis(basis, points, transform),
+                    evaluate_basis(basis, points, transform, screen_basis=False),
                     evaluate_deriv_basis(basis, points, np.array([0, 1, 0]), transform),
                 ),
                 np.einsum(
                     "ij,ik,jk->k",
                     density,
                     evaluate_deriv_basis(basis, points, np.array([0, 0, 1]), transform),
-                    evaluate_basis(basis, points, transform),
+                    evaluate_basis(basis, points, transform, screen_basis=False),
                 )
                 + np.einsum(
                     "ij,ik,jk->k",
                     density,
-                    evaluate_basis(basis, points, transform),
+                    evaluate_basis(basis, points, transform, screen_basis=False),
                     evaluate_deriv_basis(basis, points, np.array([0, 0, 1]), transform),
                 ),
             ]
diff --git a/tests/test_eval.py b/tests/test_eval.py
index 2bcb3a88..7c681adf 100644
--- a/tests/test_eval.py
+++ b/tests/test_eval.py
@@ -2,7 +2,7 @@
 from gbasis.contractions import GeneralizedContractionShell
 from gbasis.evals._deriv import _eval_deriv_contractions
 from gbasis.evals.eval import Eval, evaluate_basis
-from gbasis.parsers import make_contractions, parse_nwchem
+from gbasis.parsers import make_contractions, parse_nwchem, parse_gbs
 from gbasis.utils import factorial2
 import numpy as np
 import pytest
@@ -60,7 +60,7 @@ def test_evaluate_basis_cartesian():
     evaluate_obj = Eval(basis)
     assert np.allclose(
         evaluate_obj.construct_array_cartesian(points=np.array([[0, 0, 0]])),
-        evaluate_basis(basis, np.array([[0, 0, 0]])),
+        evaluate_basis(basis, np.array([[0, 0, 0]]), screen_basis=False),
     )
 
 
@@ -73,21 +73,21 @@ def test_evaluate_basis_spherical():
     evaluate_obj = Eval(basis)
     assert np.allclose(
         evaluate_obj.construct_array_cartesian(points=np.array([[0, 0, 0]])),
-        evaluate_basis(basis, np.array([[0, 0, 0]])),
+        evaluate_basis(basis, np.array([[0, 0, 0]]), screen_basis=False),
     )
     # p orbitals are zero at center
     basis = make_contractions(basis_dict, ["Li"], np.array([[0, 0, 0]]), "spherical")
     evaluate_obj = Eval(basis)
     assert np.allclose(
         evaluate_obj.construct_array_cartesian(points=np.array([[0, 0, 0]])),
-        evaluate_basis(basis, np.array([[0, 0, 0]])),
+        evaluate_basis(basis, np.array([[0, 0, 0]]), screen_basis=False),
     )
 
     basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "spherical")
     evaluate_obj = Eval(basis)
     assert np.allclose(
         evaluate_obj.construct_array_spherical(points=np.array([[1, 1, 1]])),
-        evaluate_basis(basis, np.array([[1, 1, 1]])),
+        evaluate_basis(basis, np.array([[1, 1, 1]]), screen_basis=False),
     )
 
 
@@ -101,8 +101,8 @@ def test_evaluate_basis_mix():
         basis_dict, ["H"], np.array([[0, 0, 0]]), ["spherical"]
     )
     assert np.allclose(
-        evaluate_basis(spherical_basis, np.array([[0, 0, 0]])),
-        evaluate_basis(spherical_basis_list, np.array([[0, 0, 0]])),
+        evaluate_basis(spherical_basis, np.array([[0, 0, 0]]), screen_basis=False),
+        evaluate_basis(spherical_basis_list, np.array([[0, 0, 0]]), screen_basis=False),
     )
 
     cartesian_basis = make_contractions(basis_dict, ["H"], np.array([[0, 0, 0]]), "cartesian")
@@ -110,8 +110,8 @@ def test_evaluate_basis_mix():
         basis_dict, ["H"], np.array([[0, 0, 0]]), ["cartesian"]
     )
     assert np.allclose(
-        evaluate_basis(cartesian_basis, np.array([[0, 0, 0]])),
-        evaluate_basis(cartesian_basis_list, np.array([[0, 0, 0]])),
+        evaluate_basis(cartesian_basis, np.array([[0, 0, 0]]), screen_basis=False),
+        evaluate_basis(cartesian_basis_list, np.array([[0, 0, 0]]), screen_basis=False),
     )
 
     spherical_basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "spherical")
@@ -119,8 +119,8 @@ def test_evaluate_basis_mix():
         basis_dict, ["Kr"], np.array([[0, 0, 0]]), ["spherical"] * 8
     )
     assert np.allclose(
-        evaluate_basis(spherical_basis, np.array([[1, 1, 1]])),
-        evaluate_basis(spherical_basis_list, np.array([[1, 1, 1]])),
+        evaluate_basis(spherical_basis, np.array([[1, 1, 1]]), screen_basis=False),
+        evaluate_basis(spherical_basis_list, np.array([[1, 1, 1]]), screen_basis=False),
     )
 
     cartesian_basis = make_contractions(basis_dict, ["Kr"], np.array([[0, 0, 0]]), "cartesian")
@@ -128,8 +128,8 @@ def test_evaluate_basis_mix():
         basis_dict, ["Kr"], np.array([[0, 0, 0]]), ["cartesian"] * 8
     )
     assert np.allclose(
-        evaluate_basis(cartesian_basis, np.array([[1, 1, 1]])),
-        evaluate_basis(cartesian_basis_list, np.array([[1, 1, 1]])),
+        evaluate_basis(cartesian_basis, np.array([[1, 1, 1]]), screen_basis=False),
+        evaluate_basis(cartesian_basis_list, np.array([[1, 1, 1]]), screen_basis=False),
     )
 
 
@@ -143,7 +143,7 @@ def test_evaluate_basis_lincomb():
         evaluate_obj.construct_array_lincomb(
             transform, ["spherical"], points=np.array([[1, 1, 1]])
         ),
-        evaluate_basis(basis, np.array([[1, 1, 1]]), transform=transform),
+        evaluate_basis(basis, np.array([[1, 1, 1]]), transform=transform, screen_basis=False),
     )
 
 
@@ -166,8 +166,12 @@ def test_evaluate_basis_horton():
     grid_x, grid_y, grid_z = np.meshgrid(grid_1d, grid_1d, grid_1d)
     grid_3d = np.vstack([grid_x.ravel(), grid_y.ravel(), grid_z.ravel()]).T
 
-    assert np.allclose(evaluate_basis(cartesian_basis, grid_3d), horton_eval_cart.T)
-    assert np.allclose(evaluate_basis(spherical_basis, grid_3d), horton_eval_sph.T)
+    assert np.allclose(
+        evaluate_basis(cartesian_basis, grid_3d, screen_basis=False), horton_eval_cart.T
+    )
+    assert np.allclose(
+        evaluate_basis(spherical_basis, grid_3d, screen_basis=False), horton_eval_sph.T
+    )
 
 
 def test_evaluate_basis_pyscf():
@@ -195,21 +199,47 @@ def test_evaluate_basis_pyscf():
     pyscf_eval_cart = gto.eval_gto(mol, "GTOval_cart", grid_3d)
 
     # s orbitals
-    assert np.allclose(evaluate_basis(cartesian_basis, grid_3d)[:6], pyscf_eval_cart.T[:6])
-    assert np.allclose(evaluate_basis(cartesian_basis, grid_3d)[46:53], pyscf_eval_cart.T[46:53])
-    assert np.allclose(evaluate_basis(spherical_basis, grid_3d)[:6], pyscf_eval_sph.T[:6])
-    assert np.allclose(evaluate_basis(spherical_basis, grid_3d)[40:47], pyscf_eval_sph.T[40:47])
+    assert np.allclose(
+        evaluate_basis(cartesian_basis, grid_3d, screen_basis=False)[:6], pyscf_eval_cart.T[:6]
+    )
+    assert np.allclose(
+        evaluate_basis(cartesian_basis, grid_3d, screen_basis=False)[46:53],
+        pyscf_eval_cart.T[46:53],
+    )
+    assert np.allclose(
+        evaluate_basis(spherical_basis, grid_3d, screen_basis=False)[:6], pyscf_eval_sph.T[:6]
+    )
+    assert np.allclose(
+        evaluate_basis(spherical_basis, grid_3d, screen_basis=False)[40:47], pyscf_eval_sph.T[40:47]
+    )
     # p orbitals
-    assert np.allclose(evaluate_basis(cartesian_basis, grid_3d)[6:18], pyscf_eval_cart.T[6:18])
-    assert np.allclose(evaluate_basis(cartesian_basis, grid_3d)[53:65], pyscf_eval_cart.T[53:65])
-    assert np.allclose(evaluate_basis(spherical_basis, grid_3d)[6:18], pyscf_eval_sph.T[6:18])
-    assert np.allclose(evaluate_basis(spherical_basis, grid_3d)[47:59], pyscf_eval_sph.T[47:59])
+    assert np.allclose(
+        evaluate_basis(cartesian_basis, grid_3d, screen_basis=False)[6:18], pyscf_eval_cart.T[6:18]
+    )
+    assert np.allclose(
+        evaluate_basis(cartesian_basis, grid_3d, screen_basis=False)[53:65],
+        pyscf_eval_cart.T[53:65],
+    )
+    assert np.allclose(
+        evaluate_basis(spherical_basis, grid_3d, screen_basis=False)[6:18], pyscf_eval_sph.T[6:18]
+    )
+    assert np.allclose(
+        evaluate_basis(spherical_basis, grid_3d, screen_basis=False)[47:59], pyscf_eval_sph.T[47:59]
+    )
     # d orbitals are off by some constant for the cartesian case
-    assert np.allclose(evaluate_basis(spherical_basis, grid_3d)[18:33], pyscf_eval_sph.T[18:33])
-    assert np.allclose(evaluate_basis(spherical_basis, grid_3d)[59:74], pyscf_eval_sph.T[59:74])
+    assert np.allclose(
+        evaluate_basis(spherical_basis, grid_3d, screen_basis=False)[18:33], pyscf_eval_sph.T[18:33]
+    )
+    assert np.allclose(
+        evaluate_basis(spherical_basis, grid_3d, screen_basis=False)[59:74], pyscf_eval_sph.T[59:74]
+    )
     # f orbitals are off by some constant for the cartesian case
-    assert np.allclose(evaluate_basis(spherical_basis, grid_3d)[33:40], pyscf_eval_sph.T[33:40])
-    assert np.allclose(evaluate_basis(spherical_basis, grid_3d)[74:88], pyscf_eval_sph.T[74:88])
+    assert np.allclose(
+        evaluate_basis(spherical_basis, grid_3d, screen_basis=False)[33:40], pyscf_eval_sph.T[33:40]
+    )
+    assert np.allclose(
+        evaluate_basis(spherical_basis, grid_3d, screen_basis=False)[74:88], pyscf_eval_sph.T[74:88]
+    )
 
 
 @pytest.mark.xfail
@@ -237,15 +267,39 @@ def test_evaluate_basis_pyscf_cart_norm():
 
     # d orbitals are all off by some scalar factor
     assert np.allclose(
-        evaluate_basis(basis, grid_3d, coord_type="cartesian")[18:36], pyscf_eval_cart.T[18:36]
+        evaluate_basis(basis, grid_3d, coord_type="cartesian", screen_basis=False)[18:36],
+        pyscf_eval_cart.T[18:36],
     )
     assert np.allclose(
-        evaluate_basis(basis, grid_3d, coord_type="cartesian")[65:83], pyscf_eval_cart.T[65:83]
+        evaluate_basis(basis, grid_3d, coord_type="cartesian", screen_basis=False)[65:83],
+        pyscf_eval_cart.T[65:83],
     )
     # f orbitals are all off by some scalar factor
     assert np.allclose(
-        evaluate_basis(basis, grid_3d, coord_type="cartesian")[36:46], pyscf_eval_cart.T[36:46]
+        evaluate_basis(basis, grid_3d, coord_type="cartesian", screen_basis=False)[36:46],
+        pyscf_eval_cart.T[36:46],
     )
     assert np.allclose(
-        evaluate_basis(basis, grid_3d, coord_type="cartesian")[83:103], pyscf_eval_cart.T[83:103]
+        evaluate_basis(basis, grid_3d, coord_type="cartesian", screen_basis=False)[83:103],
+        pyscf_eval_cart.T[83:103],
     )
+
+
+@pytest.mark.parametrize("precision", [1.0e-5, 1.0e-6, 1.0e-7, 1.0e-8])
+def test_evaluate_basis_screening_accuracy(precision):
+    """Test basis set evaluation screening."""
+
+    basis_dict = parse_gbs(find_datafile("data_631g.gbs"))
+    atsymbols = ["H", "C", "Kr"]
+    atcoords = np.array([[0, 0, 0], [1, 1, 1], [2, 2, 2]])
+    contraction = make_contractions(basis_dict, atsymbols, atcoords, "cartesian")
+
+    grid_1d = np.linspace(-2, 2, num=5)
+    grid_x, grid_y, grid_z = np.meshgrid(grid_1d, grid_1d, grid_1d)
+    grid_3d = np.vstack([grid_x.ravel(), grid_y.ravel(), grid_z.ravel()]).T
+
+    #  the screening tolerance needs to be 1e-4 times the desired precision
+    tol_screen = precision * 1e-4
+    basis_evaluation = evaluate_basis(contraction, grid_3d, tol_screen=tol_screen)
+    basis_evaluation_no_screen = evaluate_basis(contraction, grid_3d, screen_basis=False)
+    assert np.allclose(basis_evaluation, basis_evaluation_no_screen, atol=precision)
diff --git a/tests/test_screening.py b/tests/test_screening.py
new file mode 100644
index 00000000..85bbed88
--- /dev/null
+++ b/tests/test_screening.py
@@ -0,0 +1,68 @@
+"""Test gbasis.screening"""
+import numpy as np
+import pytest
+from gbasis.parsers import make_contractions, parse_nwchem
+from gbasis.screening import is_two_index_overlap_screened
+from utils import find_datafile
+from gbasis.utils import factorial2
+from gbasis.screening import compute_primitive_cutoff_radius
+
+
+def get_atom_contractions_data(atsym, atcoords):
+    """Get the STO-6G contractions for a given atom symbol and coordinates."""
+    basis_dict = parse_nwchem(find_datafile("data_sto6g.nwchem"))
+    basis = make_contractions(basis_dict, [atsym], atcoords, "cartesian")
+    return basis
+
+
+@pytest.mark.parametrize("bond_length", [0, 0.999, 2.0, 4.0, 8.0, 50.0, 100.0])
+@pytest.mark.parametrize("tol_screen", [1e-4, 1e-8, 1e-12])
+def test_is_two_index_overlap_screened(bond_length, tol_screen):
+    contractions_one = get_atom_contractions_data("H", np.array([[0, 0, 0]]))
+    contractions_two = get_atom_contractions_data("O", np.array([[0, 0, bond_length]]))
+    screen_pairs_list = []
+    screening_cutoffs = []
+    for h_shell in contractions_one:
+        for o_shell in contractions_two:
+            screen_pairs_list.append(is_two_index_overlap_screened(h_shell, o_shell, tol_screen))
+            alpha_a = min(h_shell.exps)
+            alpha_b = min(o_shell.exps)
+            cutoff = np.sqrt(-(alpha_a + alpha_b) / (alpha_a * alpha_b) * np.log(tol_screen))
+            screening_cutoffs.append(cutoff)
+
+    # bonds too close, not any pairs should be screened
+    if bond_length < 1.0:
+        assert not any(screen_pairs_list)
+    # bonds too far, all pairs should be screened
+    elif bond_length > 30.0:
+        assert all(screen_pairs_list)
+    # intermediate bond lengths, some pairs should be screened
+    else:
+        ref_screened = np.array(screening_cutoffs) < bond_length
+        assert np.all(
+            np.array(screen_pairs_list) == ref_screened
+        ), "Screening results do not match the expected values based on cutoff distances."
+
+
+@pytest.mark.parametrize("angm", [0, 1, 2, 3])
+@pytest.mark.parametrize("alpha", [0.5, 1.0, 2.0])
+@pytest.mark.parametrize("coeff", [0.1, 0.4, 0.8])
+@pytest.mark.parametrize("tol_screen", [1e-4, 1e-8])
+def test_compute_primitive_cutoff_radius(angm, alpha, coeff, tol_screen):
+    """Test the computation of the primitive cutoff radius."""
+
+    def compute_primitive_value(r, c, alpha, angm):
+        """Compute the primitive value at the given radius."""
+        n = (
+            (2 * alpha / np.pi) ** 0.25
+            * (4 * alpha) ** (angm / 2)
+            / np.sqrt(factorial2(2 * angm + 1))
+        )
+        return c * n * np.exp(-alpha * r**2)
+
+    cutoff_r = compute_primitive_cutoff_radius(coeff, alpha, angm, tol_screen)
+    val_over_cutoff = compute_primitive_value(cutoff_r + 1e-3, coeff, alpha, angm)
+
+    assert (
+        val_over_cutoff < tol_screen
+    ), f"Value {val_over_cutoff} at r={cutoff_r + 1e-3} is not below the tolerance {tol_screen}"