diff --git a/nautilus/bounds/__init__.py b/nautilus/bounds/__init__.py
index 0c93117..5ed4b09 100644
--- a/nautilus/bounds/__init__.py
+++ b/nautilus/bounds/__init__.py
@@ -1,8 +1,10 @@
 """Modules implementing various multi-dimensional bounds."""
 
 from .basic import UnitCube, Ellipsoid, UnitCubeEllipsoidMixture
-from .neural import NeuralBound, NautilusBound
+from .nautilus import NautilusBound
+from .neural import NeuralBound
+from .periodic import PhaseShift
 from .union import Union
 
-__all__ = ['UnitCube', 'Ellipsoid', 'UnitCubeEllipsoidMixture', 'Union',
-           'NeuralBound', 'NautilusBound']
+__all__ = ['Ellipsoid', 'NautilusBound', 'NeuralBound',
+           'PhaseShift', 'Union', 'UnitCube', 'UnitCubeEllipsoidMixture']
diff --git a/nautilus/bounds/nautilus.py b/nautilus/bounds/nautilus.py
new file mode 100644
index 0000000..d4282a8
--- /dev/null
+++ b/nautilus/bounds/nautilus.py
@@ -0,0 +1,394 @@
+"""Module implementing the nautilus bound."""
+
+import numpy as np
+from functools import partial
+from threadpoolctl import threadpool_limits
+
+from .basic import Ellipsoid, UnitCubeEllipsoidMixture
+from .neural import NeuralBound
+from .periodic import PhaseShift
+from .union import Union
+
+
+class NautilusBound():
+    """Union of multiple non-overlapping neural network-based bounds.
+
+    The bound is the overlap of the union of multiple neural network-based
+    bounds and the unit hypercube.
+
+    Attributes
+    ----------
+    n_dim : int
+        Number of dimensions.
+    neural_bounds : list
+        List of the individual neural network-based bounds.
+    outer_bound : Union
+        Outer bound used for sampling.
+    rng : None or numpy.random.Generator
+        Determines random number generation.
+    points : numpy.ndarray
+        Points that a call to `sample` will return next.
+    n_sample : int
+        Number of points sampled from the outer bound.
+    n_reject : int
+        Number of points rejected due to not falling into the neural
+        network-based bounds.
+    """
+
+    @classmethod
+    def compute(cls, points, log_l, log_l_min, log_v_target,
+                enlarge_per_dim=1.1, n_points_min=None, split_threshold=100,
+                periodic=None, n_networks=4, neural_network_kwargs={},
+                pool=None, rng=None):
+        """Compute a union of multiple neural network-based bounds.
+
+        Parameters
+        ----------
+        points : numpy.ndarray with shape (m, n)
+            A 2-D array where each row represents a point.
+        log_l : numpy.ndarray of length m
+            Likelihood of each point.
+        log_l_min : float
+            Target likelihood threshold of the bound.
+        log_v_target : float
+            Expected volume of the bound. Used for multi-ellipsoidal
+            decomposition.
+        enlarge_per_dim : float, optional
+            Along each dimension, the ellipsoid of the outer bound is enlarged
+            by this factor. Default is 1.1.
+        n_points_min : int or None, optional
+            The minimum number of points each ellipsoid should have.
+            Effectively, ellipsoids with less than twice that number will not
+            be split further. If None, uses `n_points_min = n_dim + 1`. Default
+            is None.
+        split_threshold: float, optional
+            Threshold used for splitting the multi-ellipsoidal bound used for
+            sampling. If the volume of the bound is larger than
+            `split_threshold` times the target volume, the multi-ellipsiodal
+            bound is split further, if possible. Default is 100.
+        periodic : numpy.ndarray or None
+            Indices of the parameters that are periodic.
+        n_networks : int, optional
+            Number of networks used in the emulator. Default is 4.
+        neural_network_kwargs : dict, optional
+            Non-default keyword arguments passed to the constructor of
+            MLPRegressor.
+        pool : nautilus.pool.NautilusPool or None, optional
+            Pool used for parallel processing. Default is None.
+        rng : None or numpy.random.Generator, optional
+            Determines random number generation. Default is None.
+
+        Returns
+        -------
+        bound : NautilusBound
+            The bound.
+
+        """
+        bound = cls()
+        bound.n_dim = points.shape[1]
+
+        if periodic is not None:
+            bound.shift = PhaseShift.compute(points[log_l >= log_l_min],
+                                             periodic)
+            points = bound.shift.transform(points)
+        else:
+            bound.shift = None
+
+        bound.neural_bounds = []
+
+        multi_ellipsoid = Union.compute(
+            points[log_l >= log_l_min], enlarge_per_dim=enlarge_per_dim,
+            n_points_min=n_points_min, bound_class=Ellipsoid,
+            rng=rng)
+
+        while multi_ellipsoid.split(allow_overlap=False):
+            pass
+
+        for ellipsoid in multi_ellipsoid.bounds:
+            select = ellipsoid.contains(points)
+            bound.neural_bounds.append(NeuralBound.compute(
+                points[select], log_l[select], log_l_min,
+                enlarge_per_dim=enlarge_per_dim, n_networks=n_networks,
+                neural_network_kwargs=neural_network_kwargs, pool=pool,
+                rng=rng))
+
+        bound.outer_bound = Union.compute(
+            points[log_l >= log_l_min], enlarge_per_dim=enlarge_per_dim,
+            n_points_min=n_points_min, bound_class=UnitCubeEllipsoidMixture,
+            rng=rng)
+
+        # If the single bounding ellipsoid is too large, split ellipsoids
+        # further.
+        while bound.outer_bound.log_v - log_v_target > np.log(
+                split_threshold * enlarge_per_dim**points.shape[1]):
+            if not bound.outer_bound.split():
+                break
+
+        # If the ellipsoid union is still too large, check whether some
+        # ellipsoids have too low densities and should be dropped.
+        while bound.outer_bound.log_v - log_v_target > np.log(
+                split_threshold * enlarge_per_dim**points.shape[1]):
+            if not bound.outer_bound.trim():
+                break
+
+        if rng is None:
+            bound.rng = np.random.default_rng()
+        else:
+            bound.rng = rng
+
+        bound.points = np.zeros((0, points.shape[1]))
+        bound.n_sample = 0
+        bound.n_reject = 0
+
+        return bound
+
+    def contains(self, points):
+        """Check whether points are contained in the bound.
+
+        Parameters
+        ----------
+        points : numpy.ndarray
+            A 1-D or 2-D array containing single point or a collection of
+            points. If more than one-dimensional, each row represents a point.
+
+        Returns
+        -------
+        in_bound : bool or numpy.ndarray
+            Bool or array of bools describing for each point it is contained in
+            the bound.
+
+        """
+        if self.shift is not None:
+            points = self.shift.transform(points)
+        in_bound = self.outer_bound.contains(points)
+        if len(self.neural_bounds) > 0:
+            in_bound = in_bound & np.any(
+                [bound.contains(points) for bound in self.neural_bounds],
+                axis=0)
+        return in_bound
+
+    @threadpool_limits.wrap(limits=1)
+    def _reset_and_sample(self, n_points=100, rng=None):
+        """Reset the bound, sample points internally and return the result.
+
+        Parameters
+        ----------
+        n_points : int, optional
+            How many points to draw. Default is 100.
+        rng : None or numpy.random.Generator, optional
+            Determines random number generation. If None, random number
+            generation is not reset. Default is None.
+
+        Returns
+        -------
+        bound : NautilusBound
+            The bound.
+
+        """
+        self.reset(rng=rng)
+        self.sample(n_points=n_points, return_points=False)
+        return self
+
+    def sample(self, n_points=100, return_points=True, pool=None):
+        """Sample points from the the bound.
+
+        Parameters
+        ----------
+        n_points : int, optional
+            How many points to draw. Default is 100.
+        return_points : bool, optional
+            If True, return sampled points. Otherwise, sample internally until
+            at least `n_points` are saved.
+        pool : nautilus.pool.NautilusPool or None, optional
+            Pool used for parallel processing. Default is None.
+
+        Returns
+        -------
+        points : numpy.ndarray
+            Points as two-dimensional array of shape (n_points, n_dim).
+
+        """
+        if len(self.points) < n_points:
+            if pool is None:
+                while len(self.points) < n_points:
+                    n_sample = 1000
+                    points = self.outer_bound.sample(n_sample)
+                    in_bound = np.any([bound.contains(points) for bound in
+                                       self.neural_bounds], axis=0)
+                    points = points[in_bound]
+                    self.points = np.vstack([self.points, points])
+                    self.n_sample += n_sample
+                    self.n_reject += n_sample - len(points)
+            else:
+                n_jobs = pool.size
+                n_points_per_job = (
+                    (max(n_points - len(self.points), 10000)) // n_jobs) + 1
+                func = partial(self._reset_and_sample, n_points_per_job)
+                rngs = [np.random.default_rng(seed) for seed in
+                        np.random.SeedSequence(self.rng.integers(
+                            2**32 - 1)).spawn(n_jobs)]
+                bounds = pool.map(func, rngs)
+                for bound in bounds:
+                    self.points = np.vstack([self.points, bound.points])
+                    self.n_sample += bound.n_sample
+                    self.n_reject += bound.n_reject
+                    self.outer_bound.n_sample += bound.outer_bound.n_sample
+                    self.outer_bound.n_reject += bound.outer_bound.n_reject
+
+        if return_points:
+            points = self.points[:n_points]
+            self.points = self.points[n_points:]
+            if self.shift is not None:
+                points = self.shift.transform(points, inverse=True)
+            return points
+
+    @property
+    def log_v(self):
+        """Return the natural log of the volume.
+
+        Returns
+        -------
+        log_v : float
+            An estimate of the natural log of the volume. Will become more
+            accurate as more points are sampled.
+
+        """
+        if self.n_sample == 0:
+            self.sample(return_points=False)
+
+        return self.outer_bound.log_v + np.log(
+            1.0 - self.n_reject / self.n_sample)
+
+    @property
+    def n_ell(self):
+        """Return the number of ellipsoids in the bound.
+
+        Returns
+        -------
+        n_ell : int
+            The number of ellipsoids.
+        """
+        return np.sum([np.any(~bound.dim_cube) for bound in
+                       self.outer_bound.bounds])
+
+    @property
+    def n_net(self):
+        """Return the number of neural networks in the bound.
+
+        Returns
+        -------
+        n_net : int
+            The number of neural networks.
+        """
+        if self.neural_bounds[0].emulator is not None:
+            return len(self.neural_bounds) * len(
+                self.neural_bounds[0].emulator.neural_networks)
+        else:
+            return 0
+
+    def write(self, group):
+        """Write the bound to an HDF5 group.
+
+        Parameters
+        ----------
+        group : h5py.Group
+            HDF5 group to write to.
+
+        """
+        group.attrs['type'] = 'NautilusBound'
+        group.attrs['n_dim'] = self.n_dim
+
+        if self.shift is not None:
+            self.shift.write(group.create_group('shift'))
+
+        group.attrs['n_neural_bounds'] = len(self.neural_bounds)
+
+        for i, neural_bound in enumerate(self.neural_bounds):
+            neural_bound.write(group.create_group('neural_bound_{}'.format(i)))
+
+        self.outer_bound.write(group.create_group('outer_bound'))
+
+        group.create_dataset('points', data=self.points,
+                             maxshape=(None, self.n_dim))
+        group.attrs['n_sample'] = self.n_sample
+        group.attrs['n_reject'] = self.n_reject
+
+    def update(self, group):
+        """Update bound information previously written to an HDF5 group.
+
+        Parameters
+        ----------
+        group : h5py.Group
+            HDF5 group to write to.
+
+        """
+        group.attrs['n_sample'] = self.n_sample
+        group.attrs['n_reject'] = self.n_reject
+        self.outer_bound.update(group['outer_bound'])
+        group['points'].resize(self.points.shape)
+        group['points'][...] = self.points
+
+    @classmethod
+    def read(cls, group, rng=None):
+        """Read the bound from an HDF5 group.
+
+        Parameters
+        ----------
+        group : h5py.Group
+            HDF5 group to write to.
+        rng : None or numpy.random.Generator, optional
+            Determines random number generation. Default is None.
+
+        Returns
+        -------
+        bound : NautilusBound
+            The bound.
+
+        """
+        bound = cls()
+
+        if rng is None:
+            bound.rng = np.random.default_rng()
+        else:
+            bound.rng = rng
+
+        bound.n_dim = group.attrs['n_dim']
+
+        if 'shift' in group:
+            bound.shift = PhaseShift.read(group['shift'])
+        else:
+            bound.shift = None
+
+        bound.neural_bounds = []
+        i = 0
+        while 'neural_bound_{}'.format(i) in group:
+            bound.neural_bounds.append(NeuralBound.read(
+                group['neural_bound_{}'.format(i)],
+                rng=bound.rng))
+            i += 1
+
+        bound.outer_bound = Union.read(
+            group['outer_bound'], rng=rng)
+
+        bound.points = np.array(group['points'])
+        bound.n_sample = group.attrs['n_sample']
+        bound.n_reject = group.attrs['n_reject']
+
+        return bound
+
+    def reset(self, rng=None):
+        """Reset random number generation and any progress, if applicable.
+
+        Parameters
+        ----------
+        rng : None or numpy.random.Generator, optional
+            Determines random number generation. If None, random number
+            generation is not reset. Default is None.
+
+        """
+        self.points = np.zeros((0, self.n_dim))
+        self.n_sample = 0
+        self.n_reject = 0
+        self.outer_bound.reset(rng)
+        if rng is not None:
+            self.rng = rng
diff --git a/nautilus/bounds/neural.py b/nautilus/bounds/neural.py
index 57933c7..4c71293 100644
--- a/nautilus/bounds/neural.py
+++ b/nautilus/bounds/neural.py
@@ -1,13 +1,9 @@
 """Module implementing multi-dimensional neural network-based bounds."""
 
 import numpy as np
-from functools import partial
 from scipy.stats import rankdata
-from threadpoolctl import threadpool_limits
 
-from .basic import Ellipsoid, UnitCubeEllipsoidMixture
-from .union import Union
-from .periodic import PhaseShift
+from .basic import Ellipsoid
 from ..neural import NeuralNetworkEmulator
 
 
@@ -174,383 +170,3 @@ def read(cls, group, rng=None):
             bound.emulator = None
 
         return bound
-
-
-class NautilusBound():
-    """Union of multiple non-overlapping neural network-based bounds.
-
-    The bound is the overlap of the union of multiple neural network-based
-    bounds and the unit hypercube.
-
-    Attributes
-    ----------
-    n_dim : int
-        Number of dimensions.
-    neural_bounds : list
-        List of the individual neural network-based bounds.
-    outer_bound : Union
-        Outer bound used for sampling.
-    rng : None or numpy.random.Generator
-        Determines random number generation.
-    points : numpy.ndarray
-        Points that a call to `sample` will return next.
-    n_sample : int
-        Number of points sampled from the outer bound.
-    n_reject : int
-        Number of points rejected due to not falling into the neural
-        network-based bounds.
-    """
-
-    @classmethod
-    def compute(cls, points, log_l, log_l_min, log_v_target,
-                enlarge_per_dim=1.1, n_points_min=None, split_threshold=100,
-                periodic=None, n_networks=4, neural_network_kwargs={},
-                pool=None, rng=None):
-        """Compute a union of multiple neural network-based bounds.
-
-        Parameters
-        ----------
-        points : numpy.ndarray with shape (m, n)
-            A 2-D array where each row represents a point.
-        log_l : numpy.ndarray of length m
-            Likelihood of each point.
-        log_l_min : float
-            Target likelihood threshold of the bound.
-        log_v_target : float
-            Expected volume of the bound. Used for multi-ellipsoidal
-            decomposition.
-        enlarge_per_dim : float, optional
-            Along each dimension, the ellipsoid of the outer bound is enlarged
-            by this factor. Default is 1.1.
-        n_points_min : int or None, optional
-            The minimum number of points each ellipsoid should have.
-            Effectively, ellipsoids with less than twice that number will not
-            be split further. If None, uses `n_points_min = n_dim + 1`. Default
-            is None.
-        split_threshold: float, optional
-            Threshold used for splitting the multi-ellipsoidal bound used for
-            sampling. If the volume of the bound is larger than
-            `split_threshold` times the target volume, the multi-ellipsiodal
-            bound is split further, if possible. Default is 100.
-        periodic : numpy.ndarray or None
-            Indices of the parameters that are periodic.
-        n_networks : int, optional
-            Number of networks used in the emulator. Default is 4.
-        neural_network_kwargs : dict, optional
-            Non-default keyword arguments passed to the constructor of
-            MLPRegressor.
-        pool : nautilus.pool.NautilusPool or None, optional
-            Pool used for parallel processing. Default is None.
-        rng : None or numpy.random.Generator, optional
-            Determines random number generation. Default is None.
-
-        Returns
-        -------
-        bound : NautilusBound
-            The bound.
-
-        """
-        bound = cls()
-        bound.n_dim = points.shape[1]
-
-        if periodic is not None:
-            bound.shift = PhaseShift.compute(points[log_l >= log_l_min],
-                                             periodic)
-            points = bound.shift.transform(points)
-        else:
-            bound.shift = None
-
-        bound.neural_bounds = []
-
-        multi_ellipsoid = Union.compute(
-            points[log_l >= log_l_min], enlarge_per_dim=enlarge_per_dim,
-            n_points_min=n_points_min, bound_class=Ellipsoid,
-            rng=rng)
-
-        while multi_ellipsoid.split(allow_overlap=False):
-            pass
-
-        for ellipsoid in multi_ellipsoid.bounds:
-            select = ellipsoid.contains(points)
-            bound.neural_bounds.append(NeuralBound.compute(
-                points[select], log_l[select], log_l_min,
-                enlarge_per_dim=enlarge_per_dim, n_networks=n_networks,
-                neural_network_kwargs=neural_network_kwargs, pool=pool,
-                rng=rng))
-
-        bound.outer_bound = Union.compute(
-            points[log_l >= log_l_min], enlarge_per_dim=enlarge_per_dim,
-            n_points_min=n_points_min, bound_class=UnitCubeEllipsoidMixture,
-            rng=rng)
-
-        # If the single bounding ellipsoid is too large, split ellipsoids
-        # further.
-        while bound.outer_bound.log_v - log_v_target > np.log(
-                split_threshold * enlarge_per_dim**points.shape[1]):
-            if not bound.outer_bound.split():
-                break
-
-        # If the ellipsoid union is still too large, check whether some
-        # ellipsoids have too low densities and should be dropped.
-        while bound.outer_bound.log_v - log_v_target > np.log(
-                split_threshold * enlarge_per_dim**points.shape[1]):
-            if not bound.outer_bound.trim():
-                break
-
-        if rng is None:
-            bound.rng = np.random.default_rng()
-        else:
-            bound.rng = rng
-
-        bound.points = np.zeros((0, points.shape[1]))
-        bound.n_sample = 0
-        bound.n_reject = 0
-
-        return bound
-
-    def contains(self, points):
-        """Check whether points are contained in the bound.
-
-        Parameters
-        ----------
-        points : numpy.ndarray
-            A 1-D or 2-D array containing single point or a collection of
-            points. If more than one-dimensional, each row represents a point.
-
-        Returns
-        -------
-        in_bound : bool or numpy.ndarray
-            Bool or array of bools describing for each point it is contained in
-            the bound.
-
-        """
-        if self.shift is not None:
-            points = self.shift.transform(points)
-        in_bound = self.outer_bound.contains(points)
-        if len(self.neural_bounds) > 0:
-            in_bound = in_bound & np.any(
-                [bound.contains(points) for bound in self.neural_bounds],
-                axis=0)
-        return in_bound
-
-    @threadpool_limits.wrap(limits=1)
-    def _reset_and_sample(self, n_points=100, rng=None):
-        """Reset the bound, sample points internally and return the result.
-
-        Parameters
-        ----------
-        n_points : int, optional
-            How many points to draw. Default is 100.
-        rng : None or numpy.random.Generator, optional
-            Determines random number generation. If None, random number
-            generation is not reset. Default is None.
-
-        Returns
-        -------
-        bound : NautilusBound
-            The bound.
-
-        """
-        self.reset(rng=rng)
-        self.sample(n_points=n_points, return_points=False)
-        return self
-
-    def sample(self, n_points=100, return_points=True, pool=None):
-        """Sample points from the the bound.
-
-        Parameters
-        ----------
-        n_points : int, optional
-            How many points to draw. Default is 100.
-        return_points : bool, optional
-            If True, return sampled points. Otherwise, sample internally until
-            at least `n_points` are saved.
-        pool : nautilus.pool.NautilusPool or None, optional
-            Pool used for parallel processing. Default is None.
-
-        Returns
-        -------
-        points : numpy.ndarray
-            Points as two-dimensional array of shape (n_points, n_dim).
-
-        """
-        if len(self.points) < n_points:
-            if pool is None:
-                while len(self.points) < n_points:
-                    n_sample = 1000
-                    points = self.outer_bound.sample(n_sample)
-                    in_bound = np.any([bound.contains(points) for bound in
-                                       self.neural_bounds], axis=0)
-                    points = points[in_bound]
-                    self.points = np.vstack([self.points, points])
-                    self.n_sample += n_sample
-                    self.n_reject += n_sample - len(points)
-            else:
-                n_jobs = pool.size
-                n_points_per_job = (
-                    (max(n_points - len(self.points), 10000)) // n_jobs) + 1
-                func = partial(self._reset_and_sample, n_points_per_job)
-                rngs = [np.random.default_rng(seed) for seed in
-                        np.random.SeedSequence(self.rng.integers(
-                            2**32 - 1)).spawn(n_jobs)]
-                bounds = pool.map(func, rngs)
-                for bound in bounds:
-                    self.points = np.vstack([self.points, bound.points])
-                    self.n_sample += bound.n_sample
-                    self.n_reject += bound.n_reject
-                    self.outer_bound.n_sample += bound.outer_bound.n_sample
-                    self.outer_bound.n_reject += bound.outer_bound.n_reject
-
-        if return_points:
-            points = self.points[:n_points]
-            self.points = self.points[n_points:]
-            if self.shift is not None:
-                points = self.shift.transform(points, inverse=True)
-            return points
-
-    @property
-    def log_v(self):
-        """Return the natural log of the volume.
-
-        Returns
-        -------
-        log_v : float
-            An estimate of the natural log of the volume. Will become more
-            accurate as more points are sampled.
-
-        """
-        if self.n_sample == 0:
-            self.sample(return_points=False)
-
-        return self.outer_bound.log_v + np.log(
-            1.0 - self.n_reject / self.n_sample)
-
-    @property
-    def n_ell(self):
-        """Return the number of ellipsoids in the bound.
-
-        Returns
-        -------
-        n_ell : int
-            The number of ellipsoids.
-        """
-        return np.sum([np.any(~bound.dim_cube) for bound in
-                       self.outer_bound.bounds])
-
-    @property
-    def n_net(self):
-        """Return the number of neural networks in the bound.
-
-        Returns
-        -------
-        n_net : int
-            The number of neural networks.
-        """
-        if self.neural_bounds[0].emulator is not None:
-            return len(self.neural_bounds) * len(
-                self.neural_bounds[0].emulator.neural_networks)
-        else:
-            return 0
-
-    def write(self, group):
-        """Write the bound to an HDF5 group.
-
-        Parameters
-        ----------
-        group : h5py.Group
-            HDF5 group to write to.
-
-        """
-        group.attrs['type'] = 'NautilusBound'
-        group.attrs['n_dim'] = self.n_dim
-        group.attrs['n_neural_bounds'] = len(self.neural_bounds)
-
-        for i, neural_bound in enumerate(self.neural_bounds):
-            neural_bound.write(group.create_group('neural_bound_{}'.format(i)))
-
-        self.outer_bound.write(group.create_group('outer_bound'))
-
-        group.create_dataset('points', data=self.points,
-                             maxshape=(None, self.n_dim))
-        group.attrs['n_sample'] = self.n_sample
-        group.attrs['n_reject'] = self.n_reject
-
-    def update(self, group):
-        """Update bound information previously written to an HDF5 group.
-
-        Parameters
-        ----------
-        group : h5py.Group
-            HDF5 group to write to.
-
-        """
-        group.attrs['n_sample'] = self.n_sample
-        group.attrs['n_reject'] = self.n_reject
-        self.outer_bound.update(group['outer_bound'])
-        group['points'].resize(self.points.shape)
-        group['points'][...] = self.points
-
-    @classmethod
-    def read(cls, group, rng=None):
-        """Read the bound from an HDF5 group.
-
-        Parameters
-        ----------
-        group : h5py.Group
-            HDF5 group to write to.
-        rng : None or numpy.random.Generator, optional
-            Determines random number generation. Default is None.
-
-        Returns
-        -------
-        bound : NautilusBound
-            The bound.
-
-        """
-        bound = cls()
-
-        if rng is None:
-            bound.rng = np.random.default_rng()
-        else:
-            bound.rng = rng
-
-        bound.n_dim = group.attrs['n_dim']
-
-        if 'shift' in group:
-            bound.shift = PhaseShift.read(group['shift'])
-        else:
-            bound.shift = None
-
-        bound.neural_bounds = []
-        i = 0
-        while 'neural_bound_{}'.format(i) in group:
-            bound.neural_bounds.append(NeuralBound.read(
-                group['neural_bound_{}'.format(i)],
-                rng=bound.rng))
-            i += 1
-
-        bound.outer_bound = Union.read(
-            group['outer_bound'], rng=rng)
-
-        bound.points = np.array(group['points'])
-        bound.n_sample = group.attrs['n_sample']
-        bound.n_reject = group.attrs['n_reject']
-
-        return bound
-
-    def reset(self, rng=None):
-        """Reset random number generation and any progress, if applicable.
-
-        Parameters
-        ----------
-        rng : None or numpy.random.Generator, optional
-            Determines random number generation. If None, random number
-            generation is not reset. Default is None.
-
-        """
-        self.points = np.zeros((0, self.n_dim))
-        self.n_sample = 0
-        self.n_reject = 0
-        self.outer_bound.reset(rng)
-        if rng is not None:
-            self.rng = rng
diff --git a/nautilus/bounds/periodic.py b/nautilus/bounds/periodic.py
index b7dd497..f77c84f 100644
--- a/nautilus/bounds/periodic.py
+++ b/nautilus/bounds/periodic.py
@@ -36,8 +36,6 @@ def compute(cls, points, periodic):
 
         """
         bound = cls()
-        if periodic is None:
-            periodic = np.zeros(0, dtype=bool)
         bound.periodic = periodic
         bound.centers = np.zeros(len(periodic))
 
@@ -67,13 +65,10 @@ def transform(self, points, inverse=False):
             Transformed points.
 
         """
-        if not np.any(self.periodic):
-            return points
-
         points_t = np.copy(points)
         for dim in self.periodic:
             points_t[:, dim] = (points_t[:, dim] + (-1 if inverse else +1) *
-                                (self.centers[dim] + 0.5)) % 1
+                                (-self.centers[dim] + 0.5)) % 1
         return points_t
 
     def write(self, group):
diff --git a/tests/test_bounds.py b/tests/test_bounds.py
index 844e15c..617f73f 100644
--- a/tests/test_bounds.py
+++ b/tests/test_bounds.py
@@ -142,7 +142,7 @@ def test_ellipsoid_volume(points_on_hypersphere_boundary):
         n_dim = points_on_hypersphere_boundary.shape[1]
         assert np.isclose(
             ell.log_v, np.log(enlarge_per_dim**n_dim * np.pi**(n_dim / 2) /
-                               gamma(n_dim / 2 + 1)))
+                              gamma(n_dim / 2 + 1)))
 
 
 def test_ellipsoid_transform(random_points_from_hypersphere):
@@ -311,6 +311,22 @@ def test_neural_bound_contains(random_points_from_hypercube):
     assert np.mean(log_l[in_bound] > log_l_min) >= 0.9
 
 
+def test_phase_shift():
+    # Test whether the phase-shift correctly centers periodic dimensions.
+
+    np.random.seed(0)
+    for i in range(100):
+        n_points = 10  # 100
+        n_dim = 2  # 10
+        points = (np.random.random(size=(n_points, n_dim)) * 0.1 +
+                  np.random.random(size=n_dim)) % 1
+        shift = bounds.PhaseShift.compute(points, np.arange(n_dim // 2))
+        assert np.amin(shift.transform(points)[:, :n_dim // 2]) >= 0.45
+        assert np.amax(shift.transform(points)[:, :n_dim // 2]) <= 0.55
+        assert np.allclose(points, shift.transform(shift.transform(
+            points), inverse=True), rtol=0, atol=1e-12)
+
+
 def test_nautilus_bound_sample_and_contains(random_points_from_hypercube):
     # Test whether the nautilus sampling and boundary work as expected.
 
diff --git a/tests/test_io.py b/tests/test_io.py
index 67e01d1..b9ffccb 100644
--- a/tests/test_io.py
+++ b/tests/test_io.py
@@ -6,8 +6,8 @@
 
 from nautilus import Sampler
 from nautilus.bounds import (
-    UnitCube, Ellipsoid, Union, UnitCubeEllipsoidMixture, NeuralBound,
-    NautilusBound)
+    Ellipsoid, NautilusBound, NeuralBound, PhaseShift, Union, UnitCube,
+    UnitCubeEllipsoidMixture)
 from nautilus.neural import NeuralNetworkEmulator
 
 
@@ -31,9 +31,9 @@ def test_neural_io(h5py_group):
                   emulator_read.predict(points))
 
 
-@pytest.mark.parametrize("bound_class", [UnitCube, Ellipsoid, Union,
-                                         UnitCubeEllipsoidMixture, NeuralBound,
-                                         NautilusBound])
+@pytest.mark.parametrize("bound_class", [Ellipsoid, NautilusBound, NeuralBound,
+                                         Union, UnitCube,
+                                         UnitCubeEllipsoidMixture])
 @pytest.mark.parametrize("rng_sync", [True, False])
 def test_bounds_io(h5py_group, bound_class, rng_sync):
     # Test that we can write and read a bound correctly. In particular, also
@@ -61,7 +61,10 @@ def test_bounds_io(h5py_group, bound_class, rng_sync):
             args = (points, log_l, log_l_min, np.log(0.5))
         kwargs = dict(n_networks=1, pool=None)
 
-    bound_write = bound_class.compute(*args, **kwargs, rng=rng)
+    if bound_class != PhaseShift:
+        kwargs['rng'] = rng
+
+    bound_write = bound_class.compute(*args, **kwargs)
     if bound_class == Union:
         bound_write.split()
 
@@ -111,7 +114,9 @@ def test_bounds_io(h5py_group, bound_class, rng_sync):
 @pytest.mark.parametrize("n_like_max", [np.inf, 500])
 @pytest.mark.parametrize("discard_exploration", [True, False])
 @pytest.mark.parametrize("n_networks", [0, 1])
-def test_sampler_io(blobs, n_like_max, discard_exploration, n_networks):
+@pytest.mark.parametrize("periodic", [None, np.arange(1)])
+def test_sampler_io(blobs, n_like_max, discard_exploration, n_networks,
+                    periodic):
     # Test that we can write and read a sampler correctly. In particular, also
     # test that the random number generator is correctly set after writing and
     # reading. Also make sure that the sampler can print out the progress.
@@ -126,13 +131,13 @@ def likelihood(x):
             return -np.linalg.norm(x - 0.5) * 0.001
 
     sampler_write = Sampler(prior, likelihood, n_dim=2, n_live=100,
-                            n_networks=n_networks, filepath='test.hdf5',
-                            resume=False, seed=0)
+                            n_networks=n_networks, periodic=periodic,
+                            filepath='test.hdf5', resume=False, seed=0)
     sampler_write.run(f_live=0.45, n_eff=1000, n_like_max=n_like_max,
                       discard_exploration=discard_exploration)
     sampler_read = Sampler(prior, likelihood, n_dim=2, n_live=100,
-                           n_networks=n_networks, filepath='test.hdf5',
-                           resume=True)
+                           n_networks=n_networks, periodic=periodic,
+                           filepath='test.hdf5', resume=True)
     assert sampler_write.log_z == sampler_read.log_z
 
     sampler_write.run(f_live=0.45, n_eff=5000,
diff --git a/tests/test_sampler.py b/tests/test_sampler.py
index 7e99dad..6faa8f7 100644
--- a/tests/test_sampler.py
+++ b/tests/test_sampler.py
@@ -371,3 +371,27 @@ def likelihood(x):
         # plateaus.
         assert np.all(np.isclose(sampler.shell_log_l_min[1:],
                       np.arange(len(sampler.bounds) - 1) + 2, rtol=0))
+
+
+@pytest.mark.parametrize("periodic", [False, True])
+def test_sampler_periodic(periodic):
+    # Test that the periodic boundary conditions work correctly. In particular,
+    # the sampler shouldn't split modes extending over a boundary.
+
+    n_dim = 2
+
+    def prior(x):
+        return x
+
+    def likelihood(x):
+        return multivariate_normal.logpdf(
+            np.abs(x - 0.5), mean=[0.5] * n_dim, cov=0.1)
+
+    sampler = Sampler(prior, likelihood, n_dim,
+                      periodic=np.arange(n_dim) if periodic else None,
+                      n_networks=0, seed=0)
+    sampler.run(verbose=True)
+    points, log_w, log_l = sampler.posterior()
+
+    for bound in sampler.bounds[1:]:
+        assert len(bound.neural_bounds) == 1 if periodic else 4