diff --git a/morfeus/buried_volume.py b/morfeus/buried_volume.py
index 3fbf969..b003918 100644
--- a/morfeus/buried_volume.py
+++ b/morfeus/buried_volume.py
@@ -115,6 +115,7 @@ class BuriedVolume:
     _density: float
     _excluded_atoms: set[int]
     _free_points: Array2DFloat
+    _method: str
     _octant_limits: dict[
         int, tuple[tuple[float, float], tuple[float, float], tuple[float, float]]
     ]
@@ -134,6 +135,7 @@ def __init__(
         density: float = 0.001,
         z_axis_atoms: Sequence[int] | None = None,
         xz_plane_atoms: Sequence[int] | None = None,
+        method: str = "projection",
     ) -> None:
         # Get center and and reortient coordinate system
         coordinates: Array2DFloat = np.array(coordinates)
@@ -178,7 +180,7 @@ def __init__(
         self._density = density
         self._all_coordinates = coordinates
 
-        # Converting element ids to atomic numbers if the are symbols
+        # Converting element ids to atomic numbers if they are symbols
         elements = convert_elements(elements, output="numbers")
 
         # Getting radii if they are not supplied
@@ -204,7 +206,12 @@ def __init__(
         self._excluded_atoms = set(excluded_atoms)
 
         # Compute buried volume
-        self._compute_buried_volume(center=center, radius=radius, density=density)
+        self._compute_buried_volume(
+            center=center, radius=radius, density=density, method=method
+        )
+
+        # Save method used
+        self._method = method
 
     def octant_analysis(self) -> "BuriedVolume":
         """Perform octant analysis of the buried volume."""
@@ -296,14 +303,16 @@ def octant_analysis(self) -> "BuriedVolume":
         return self
 
     def _compute_buried_volume(
-        self, center: ArrayLike1D, radius: float, density: float
+        self,
+        center: ArrayLike1D,
+        radius: float,
+        density: float,
+        method: str = "projection",
     ) -> None:
         """Compute buried volume."""
         center: Array1DFloat = np.array(center)
         # Construct sphere at metal center
-        sphere = Sphere(
-            center, radius, method="projection", density=density, filled=True
-        )
+        sphere = Sphere(center, radius, method=method, density=density, filled=True)
 
         # Prune sphere points which are within vdW radius of other atoms.
         tree = scipy.spatial.cKDTree(
@@ -393,6 +402,7 @@ def compute_distal_volume(
                 center=self._sphere.center,
                 radius=new_radius,
                 density=self._sphere.density,
+                method=self._method,
             )
             self.molecular_volume = temp_bv.buried_volume
             self.distal_volume = self.molecular_volume - self.buried_volume
diff --git a/morfeus/geometry.py b/morfeus/geometry.py
index 154872f..75e5c75 100644
--- a/morfeus/geometry.py
+++ b/morfeus/geometry.py
@@ -7,6 +7,8 @@
 
 import numpy as np
 from scipy.spatial.transform import Rotation
+import scipy.special
+import scipy.stats
 
 from morfeus.data import ANGSTROM_TO_BOHR
 from morfeus.typing import (
@@ -209,8 +211,8 @@ class Sphere:
         density: Area per point (Ų) for empty sphere
             and volume per point (ų) for filled sphere.
         filled: Whether a sphere with internal points should be constructed (works only
-            with method='projection')
-        method: Method for generating points: 'fibonacci', 'polar' or 'projection'
+            with method='projection' and method='sobol', which only works with filled)
+        method: Method for generating points: 'fibonacci', 'polar', 'projection' or 'sobol'
         radius: Radius (Å)
 
     Attributes:
@@ -252,6 +254,8 @@ def __init__(
             self.points = self._get_points_projected(density=density, filled=filled)
         elif method == "fibonacci":
             self.points = self._get_points_fibonacci(density=density)
+        elif method == "sobol" and filled:
+            self.points = self._get_points_sobol(density=density)
 
     @staticmethod
     def _get_cartesian_coordinates(
@@ -379,6 +383,39 @@ def _get_points_projected(
 
         return points
 
+    def _get_points_sobol(
+        self,
+        density: float,
+    ) -> Array2DFloat:
+        """Construct points on sphere surface by pseudorandom sampling.
+
+        Args:
+            density: Area per point (Ų) for empty sphere and volume
+                per point (ų) for filled sphere
+
+        Returns:
+            points: Array of surface points (Å)
+        """
+        # Calculate number of points from density of filled sphere and round to a power of 2.
+        n = self.volume / density * 6 / math.pi
+        m = int(round(math.log2(n)))
+
+        # Sample using a Sobol pseudorandom sequence
+        x = scipy.stats.norm.ppf(
+            scipy.stats.qmc.Sobol(d=3, scramble=True).random_base2(m)
+        )
+        ssq = np.sum(x**2, axis=1)
+
+        # Generate points based on gammainc distribution
+        r = self.radius
+        fr = r * scipy.special.gammainc(3 / 2, ssq / 2) ** (1 / 3) / np.sqrt(ssq)
+        points: Array2DFloat = np.multiply(x, np.tile(fr.reshape(2**m, 1), (1, 3)))
+
+        # Adjust with sphere center
+        points = points + self.center
+
+        return points
+
     def __repr__(self) -> str:
         return (
             f"{self.__class__.__name__}(center: {self.center}, "
diff --git a/morfeus/sasa.py b/morfeus/sasa.py
index 08921ca..44008cc 100644
--- a/morfeus/sasa.py
+++ b/morfeus/sasa.py
@@ -14,7 +14,13 @@
 from morfeus.geometry import Atom, Sphere
 from morfeus.io import read_geometry
 from morfeus.typing import Array1DFloat, Array2DFloat, ArrayLike1D, ArrayLike2D
-from morfeus.utils import convert_elements, get_radii, Import, requires_dependency
+from morfeus.utils import (
+    convert_elements,
+    get_connectivity_matrix,
+    get_radii,
+    Import,
+    requires_dependency,
+)
 
 if typing.TYPE_CHECKING:
     from matplotlib.colors import hex2color
@@ -132,31 +138,30 @@ def _calculate(self) -> None:
 
     def _determine_accessible_points(self) -> None:
         """Determine occluded and accessible points of each atom."""
-        # Based on distances to all other atoms (brute force).
-        for atom in self._atoms:
+        # Based on distances to other atoms.
+        radii: Array2DFloat = np.array([atom.radius for atom in self._atoms]).reshape(
+            -1, 1
+        )
+        radii_: Array1DFloat = np.array([atom.radius for atom in self._atoms])
+        coordinates: Array2DFloat = np.array(
+            [atom.coordinates for atom in self._atoms]
+        ).reshape(-1, 3)
+
+        # neighbor_masks[i] is a boolean mask that keeps neighbors of the ith atom
+        neighbor_masks = get_connectivity_matrix(
+            coordinates=coordinates, radii=radii_, scale_factor=1, matrix_type=bool
+        )
+
+        for i, atom in enumerate(self._atoms):
             # Construct sphere for atom
             sphere = Sphere(atom.coordinates, atom.radius, density=self._density)
             atom.points = sphere.points
 
-            # Select atoms that are at a distance less than the sum of radii
-            # !TODO can be vectorized
-            test_atoms = []
-            for test_atom in self._atoms:
-                if test_atom is not atom:
-                    distance = scipy.spatial.distance.euclidean(
-                        atom.coordinates, test_atom.coordinates
-                    )
-                    radii_sum = atom.radius + test_atom.radius
-                    if distance < radii_sum:
-                        test_atoms.append(test_atom)
-
-            # Select coordinates and radii for other atoms
-            test_coordinates = [test_atom.coordinates for test_atom in test_atoms]
-            test_radii = [test_atom.radius for test_atom in test_atoms]
-            test_radii: Array1DFloat = np.array(test_radii).reshape(-1, 1)
-
-            # Get distances to other atoms and subtract radii
-            if test_coordinates:
+            if np.any(neighbor_masks[i]):
+                # Select atoms that are at a distance less than the sum of radii
+                test_coordinates = coordinates[neighbor_masks[i]]
+                test_radii = radii[neighbor_masks[i]]
+
                 distances = scipy.spatial.distance.cdist(
                     test_coordinates, sphere.points
                 )
diff --git a/morfeus/utils.py b/morfeus/utils.py
index c10e409..c9431a9 100644
--- a/morfeus/utils.py
+++ b/morfeus/utils.py
@@ -7,7 +7,7 @@
 from importlib import import_module
 from numbers import Integral
 import shutil
-from typing import Any, cast, Literal, overload
+from typing import Any, cast, Literal, overload, Type
 
 import numpy as np
 from scipy.sparse import csr_matrix
@@ -340,8 +340,9 @@ def get_connectivity_matrix(
     radii: ArrayLike1D | None = None,
     radii_type: str = "pyykko",
     scale_factor: float = 1.2,
+    matrix_type: Type = int,
 ) -> Array2DInt:
-    """Get connectivity matrix from covalent radii.
+    """Get connectivity matrix from radii.
 
     Args:
         elements: Elements as atomic symbols or numbers
@@ -349,6 +350,7 @@ def get_connectivity_matrix(
         radii: Radii (Å)
         radii_type: Radii type: 'pyykko'
         scale_factor: Factor for scaling covalent radii
+        matrix_type: Type of the returned matrix: int
 
     Returns:
         connectivity_matrix: Connectivity matrix
@@ -357,17 +359,24 @@ def get_connectivity_matrix(
         RuntimeError: When neither elements nor radii given
     """
     coordinates: Array2DFloat = np.array(coordinates)
-    n_atoms = len(coordinates)
     if radii is None:
         if elements is None:
             raise RuntimeError("Either elements or radii needed.")
         elements = convert_elements(elements, output="numbers")
         radii = get_radii(elements, radii_type=radii_type)
     radii: Array1DFloat = np.array(radii)
-    distance_matrix = scipy.spatial.distance_matrix(coordinates, coordinates)
-    radii_matrix = np.add.outer(radii, radii) * scale_factor
-    connectivity_matrix = (distance_matrix < radii_matrix) - np.identity(
-        n_atoms
-    ).astype(int)
 
-    return connectivity_matrix
+    # The connectivity matrix is an upper triangular.
+    n = len(radii)
+    connectivity_matrix = np.zeros((n, n))
+    row, col = np.triu_indices(n, 1)
+
+    # Connected if distance is smaller than sum of radii * factor.
+    connectivity_matrix[row, col] = connectivity_matrix[
+        col, row
+    ] = scipy.spatial.distance.pdist(
+        coordinates
+    ) - scale_factor * scipy.spatial.distance.pdist(
+        radii.reshape(-1, 1), metric=lambda x, y: x + y
+    )
+    return (connectivity_matrix < 0).astype(matrix_type)
diff --git a/tests/test_buried_volume.py b/tests/test_buried_volume.py
index 37f1152..a3dd286 100644
--- a/tests/test_buried_volume.py
+++ b/tests/test_buried_volume.py
@@ -35,8 +35,21 @@ def test_reference(bv_data):
     percent_buried_volume_ref = float(data["buried_volume"])
     excluded_atoms = [int(i) for i in data["excluded_atoms"].split()]
     elements, coordinates = read_xyz(DATA_DIR / "xyz" / f"{idx}.xyz")
-    bv = BuriedVolume(elements, coordinates, 1, excluded_atoms=excluded_atoms)
+    bv_projection = BuriedVolume(
+        elements, coordinates, 1, excluded_atoms=excluded_atoms, method="projection"
+    )
+    bv_sobol = BuriedVolume(
+        elements, coordinates, 1, excluded_atoms=excluded_atoms, method="sobol"
+    )
 
     assert_almost_equal(
-        bv.fraction_buried_volume * 100, percent_buried_volume_ref, decimal=0
+        bv_projection.fraction_buried_volume * 100,
+        bv_sobol.fraction_buried_volume * 100,
+        decimal=0,
+    )
+    assert_almost_equal(
+        bv_projection.fraction_buried_volume * 100, percent_buried_volume_ref, decimal=0
+    )
+    assert_almost_equal(
+        bv_sobol.fraction_buried_volume * 100, percent_buried_volume_ref, decimal=0
     )