diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml
index e915b6fb958115732dcfc1e9bddbce9b881d4dda..37cf4a28f6744df507aa827042c053274af94b39 100644
--- a/.github/workflows/python-package.yml
+++ b/.github/workflows/python-package.yml
@@ -11,11 +11,13 @@ on:
 
 jobs:
   build:
-
-    runs-on: ubuntu-latest
     strategy:
+      fail-fast: false
       matrix:
         python-version: [3.6, 3.7, 3.8]
+        os: [ubuntu-latest, macos-11.0]
+
+    runs-on: ${{ matrix.os }}
 
     steps:
     - uses: actions/checkout@v2
@@ -25,23 +27,20 @@ jobs:
       with:
         python-version: ${{ matrix.python-version }}
 
-    - name: Get pip cache dir
-      id: pip-cache
+    - name: Install pre-install dependencies
       run: |
-        echo "::set-output name=dir::$(pip cache dir)"
+        python -m pip install --upgrade pip
+        pip install poetry wheel
+        poetry config virtualenvs.in-project true
 
-    - name: pip cache
+    - name: Cache venv
       uses: actions/cache@v2
       with:
-        path: ${{ steps.pip-cache.outputs.dir }}
-        key: ${{ runner.os }}-pip-${{ hashFiles('pyproject.toml') }}
-        restore-keys: |
-          ${{ runner.os }}-pip-
+        path: .venv
+        key: ${{ matrix.os }}-${{ matrix.python-version }}-venv-${{ hashFiles('poetry.lock') }}
 
-    - name: Install dependencies
+    - name: Install project dependencies
       run: |
-        python -m pip install --upgrade pip
-        pip install poetry wheel
         poetry install --no-root --extras docs --extras test
 
     - name: Test with pytest
@@ -59,3 +58,7 @@ jobs:
     - name: Check that test coverage is above 80%
       run: |
         poetry run coverage report --fail-under=80
+
+    - name: Test that distributable package builds
+      run: |
+        poetry build
diff --git a/.gitignore b/.gitignore
index 5ad6c26c8599dd3e7bc1c783061c35923ed5edba..6f36e7c0fc782b69a73352a4a555eca34e5ef4a9 100644
--- a/.gitignore
+++ b/.gitignore
@@ -15,6 +15,7 @@
 # Development and Testing
 .cache
 .coverage
+.python-version
 /cover/
 /dist/
 /docs/build/
diff --git a/README.md b/README.md
index 978c28ac3f13783e1ba039997faf6e1df58d3462..774ca1750db316566854426b533bc725ed9f37db 100644
--- a/README.md
+++ b/README.md
@@ -35,6 +35,13 @@ PyCGTOOL requires Python 3.6 or higher and may be installed using pip:
 pip install pycgtool
 ```
 
+### MDTraj on macOS
+
+On some versions macOS, with some versions of the Clang compiler, MDTraj may fail to load GROMACS XTC simulation trajectories.
+If you encounter this issue, make sure you have the latest version of MDTraj.
+
+For more information see [MDTraj/#1572](https://github.com/mdtraj/mdtraj/issues/1572).
+
 ## Usage
 
 Input to PyCGTOOL is GROMACS GRO and XTC files, along with two custom files: MAP and BND.  These files provide the atomistic-to-CG mapping and bonded topology respectively.  Example files are present in the [test/data](https://github.com/jag1g13/pycgtool/tree/master/test/data) directory.  The format of these files is described in the [full documentation](https://pycgtool.readthedocs.io/en/master/index.html).
diff --git a/pycgtool/__main__.py b/pycgtool/__main__.py
index dcd24128f415302815306d76ead07383f2398507..bdef110be40aa9127effdd70a9bbdf0e74375eb9 100755
--- a/pycgtool/__main__.py
+++ b/pycgtool/__main__.py
@@ -7,8 +7,6 @@ import pathlib
 import sys
 import typing
 
-import tqdm
-
 from .frame import Frame
 from .mapping import Mapping
 from .bondset import BondSet
@@ -76,15 +74,7 @@ def mapping_loop(frame: Frame, config) -> typing.Tuple[Frame, Mapping]:
     mapping = Mapping(config.map, config, itp_filename=config.itp)
 
     cg_frame = mapping.apply(frame)
-    cg_frame.save(get_output_filepath('gro', config))
-
-    logger.info("Beginning analysis of %d frames", config.end - config.begin)
-
-    if config.xtc:
-        # Main loop - perform mapping on every frame in trajectory
-        for _ in tqdm.trange(config.begin, config.end):
-            frame.next_frame()
-            mapping.apply(frame, cg_frame=cg_frame)
+    cg_frame.save(get_output_filepath('gro', config), frame_number=0)
 
     return cg_frame, mapping
 
@@ -99,14 +89,8 @@ def full_run(config):
     frame = Frame(
         topology_file=config.gro,
         trajectory_file=config.xtc,  # May be None
-        frame_start=config.begin)
-
-    try:
-        config.end = min(config.end, frame.numframes)
-
-    except TypeError:
-        # Value of args.end was None
-        config.end = frame.numframes
+        frame_start=config.begin,
+        frame_end=config.end)
 
     if config.map:
         cg_frame, mapping = mapping_loop(frame, config)
diff --git a/pycgtool/frame.py b/pycgtool/frame.py
index 511926d9351e54ca8e6512728d8659828c458153..d3bfe79f8d096915e813876b4ae1b2bc985dbd52 100644
--- a/pycgtool/frame.py
+++ b/pycgtool/frame.py
@@ -5,6 +5,7 @@ Both Frame and Residue are iterable. Residue is indexable with either atom numbe
 """
 
 import logging
+import pathlib
 import typing
 
 import mdtraj
@@ -34,25 +35,35 @@ class NonMatchingSystemError(ValueError):
 class Frame:
     """Load and store data from a simulation trajectory."""
     def __init__(self,
-                 topology_file: typing.Optional[str] = None,
-                 trajectory_file: typing.Optional[str] = None,
-                 frame_start: int = 0):
-        self.frame_number = frame_start
-
-        if topology_file:
+                 topology_file: typing.Optional[typing.Union[pathlib.Path,
+                                                             str]] = None,
+                 trajectory_file: typing.Optional[typing.Union[pathlib.Path,
+                                                               str]] = None,
+                 frame_start: int = 0,
+                 frame_end: typing.Optional[int] = None):
+        """Load a simulation trajectory.
+
+        :param topology_file: File containing system topology
+        :param trajectory_file: File containing simulation trajectory
+        :param frame_start: First frame of trajectory to use
+        :param frame_end: Last frame of trajectory to use
+        """
+        if topology_file is not None:
             try:
-                self._trajectory = mdtraj.load(topology_file)
+                self._trajectory = mdtraj.load(str(topology_file))
                 self._topology = self._trajectory.topology
 
-                if trajectory_file:
+                if trajectory_file is not None:
                     try:
-                        self._trajectory = mdtraj.load(trajectory_file,
+                        self._trajectory = mdtraj.load(str(trajectory_file),
                                                        top=self._topology)
 
+                        self._slice_trajectory(frame_start, frame_end)
+
                     except ValueError as exc:
                         raise NonMatchingSystemError from exc
 
-                self._load_trajectory_frame(self.frame_number)
+                self._load_trajectory()
 
             except OSError as exc:
                 if 'no loader' in str(exc):
@@ -64,39 +75,44 @@ class Frame:
             # No topology - we're probably building a CG frame
             self._topology = mdtraj.Topology()
 
-    def _load_trajectory_frame(self, frame_number) -> None:
-        """Load a trajectory frame into the frame attributes."""
+    def _slice_trajectory(
+            self,
+            frame_start: int = 0,
+            frame_end: typing.Optional[int] = None) -> mdtraj.Trajectory:
+        """Slice input simulation trajectory to a subset of frames.
+
+        :param frame_start: First frame of trajectory to use
+        :param frame_end: Last frame of trajectory to use
+        """
+        traj = self._trajectory
+
+        if frame_start != 0:
+            if frame_end is not None:
+                traj = traj[frame_start:frame_end]
+
+            else:
+                traj = traj[frame_start:]
+
+        return traj
+
+    def _load_trajectory(self) -> None:
+        """Load a trajectory into the frame attributes."""
         # Improve performance by not double indexing repeatedly
         traj = self._trajectory
 
-        xyz_frame = traj.xyz[frame_number]
         for atom in self._topology.atoms:
-            atom.coords = xyz_frame[atom.index]
+            atom.coords = traj.xyz[:, atom.index]
 
         # TODO handle non-cubic boxes
         try:
-            self.unitcell_lengths = traj.unitcell_lengths[frame_number]
-            self.unitcell_angles = traj.unitcell_angles[frame_number]
+            self.unitcell_lengths = traj.unitcell_lengths
+            self.unitcell_angles = traj.unitcell_angles
 
         except TypeError:
             self.unitcell_lengths = None
             self.unitcell_angles = None
 
-        self.time = traj.time[frame_number]
-
-    def next_frame(self) -> bool:
-        """Load the next trajectory frame into the frame attributes.
-
-        :return: Loading next frame was successful?
-        """
-        try:
-            self._load_trajectory_frame(self.frame_number + 1)
-
-        except IndexError:
-            return False
-
-        self.frame_number += 1
-        return True
+        self.time = traj.time
 
     @property
     def residues(self):
@@ -122,10 +138,9 @@ class Frame:
         return self._topology.n_atoms
 
     @property
-    def numframes(self) -> int:
-        """Number of frames in the frame trajectory."""
-        # The MDTraj trajectory has the topology file as frame 0
-        return self._trajectory.n_frames - 1
+    def n_frames(self) -> int:
+        """Number of frames in the trajectory."""
+        return self._trajectory.n_frames
 
     def add_residue(self,
                     name,
@@ -168,28 +183,39 @@ class Frame:
 
         return self._topology.add_atom(name, element, residue)
 
-    def save(self, filename, **kwargs) -> None:
+    def save(self,
+             filename: str,
+             frame_number: typing.Optional[int] = None,
+             **kwargs) -> None:
         """Write trajctory to file.
 
         :param filename: Name of output file
         """
-        self._trajectory.save(str(filename), **kwargs)
+        traj = self._trajectory
+        if frame_number is not None:
+            traj = traj.slice(frame_number)
 
-    def add_frame_to_trajectory(self) -> None:
-        """Add a new trajectory frame from the values stored as attributes on this frame."""
+        traj.save(str(filename), **kwargs)
+
+    def build_trajectory(self) -> None:
+        """Build an MDTraj trajectory from atom coordinates and the values stored as attributes on this frame."""
         xyz = np.array([atom.coords for atom in self._topology.atoms])
 
+        # We currently have axes: 0 - each atom, 1 - timestep, 2 - xyz coords
+        # Need to convert to axes: 0 - timestep, 1 - each atom, 2 - xyz coords
+        xyz = xyz.swapaxes(0, 1)
+
         optional_values = {
-            attr: None if getattr(self, attr) is None else [getattr(self, attr)]
+            attr: getattr(self, attr, None)
             for attr in {'time', 'unitcell_lengths', 'unitcell_angles'}
         }  # yapf: disable
 
-        new_frame = mdtraj.Trajectory(xyz,
-                                      topology=self._topology,
-                                      **optional_values)
+        new_trajectory = mdtraj.Trajectory(xyz,
+                                           topology=self._topology,
+                                           **optional_values)
 
         if hasattr(self, '_trajectory'):
-            self._trajectory += new_frame
+            self._trajectory += new_trajectory
 
         else:
-            self._trajectory = new_frame
+            self._trajectory = new_trajectory
diff --git a/pycgtool/mapping.py b/pycgtool/mapping.py
index 1197d501f8c0ba39d0f7f5c1076bf138403ba01d..7046654d48de7aec3d11e78ea76e02df3d50db9b 100644
--- a/pycgtool/mapping.py
+++ b/pycgtool/mapping.py
@@ -378,7 +378,6 @@ class Mapping:
             # Frame needs initialising
             cg_frame = self._cg_frame_setup(frame.residues)
 
-        cg_frame.frame_number = frame.frame_number
         cg_frame.time = frame.time
         cg_frame.unitcell_lengths = frame.unitcell_lengths
         cg_frame.unitcell_angles = frame.unitcell_angles
@@ -417,7 +416,7 @@ class Mapping:
                                                      bmap.weights,
                                                      unitcell_lengths)
 
-        cg_frame.add_frame_to_trajectory()
+        cg_frame.build_trajectory()
 
         return cg_frame
 
@@ -437,7 +436,7 @@ def calc_coords_weight(ref_coords, coords, weights, box=None):
         vectors -= box * np.rint(vectors / box)
 
     # Reshape weights array to match atom positions
-    weights_t = weights[np.newaxis].T
+    weights_t = weights[np.newaxis, np.newaxis].T
     result = np.sum(weights_t * vectors, axis=0)
     result += ref_coords
     return result
diff --git a/pycgtool/util.py b/pycgtool/util.py
index 5b6b2c375cd63bc3e658e58bc08342df68aa523c..82aff29f650f8bbff2b95f78b2e26a550119a073 100644
--- a/pycgtool/util.py
+++ b/pycgtool/util.py
@@ -1,6 +1,4 @@
-"""
-This module contains some general purpose utility functions used in PyCGTOOL.
-"""
+"""This module contains some general purpose utility functions used in PyCGTOOL."""
 
 import filecmp
 import functools
@@ -293,21 +291,37 @@ def any_starts_with(iterable: typing.Iterable, char: str) -> bool:
     return any(map(recurse, iterable))
 
 
-def compare_trajectories(*trajectory_files: typing.Iterable[typing.Union[
-    str, pathlib.Path]],
-                         topology_file: typing.Union[str, pathlib.Path],
-                         rtol: float = 0.001) -> bool:
+def load_optional_topology(
+    trajfile: typing.Union[str, pathlib.Path],
+    topfile: typing.Optional[typing.Union[str, pathlib.Path]] = None
+) -> mdtraj.Trajectory:
+    """Load an MDTraj trajectory with or without a separate topology file.
+
+    :param trajfile: Trajectory file
+    :param topfile: Corresponding topology file
+    :return: MDTraj trajectory
+    """
+    if topfile is None:
+        return mdtraj.load(str(trajfile))
+
+    return mdtraj.load(str(trajfile), top=str(topfile))
+
+
+def compare_trajectories(
+        *trajectory_files: typing.Iterable[typing.Union[str, pathlib.Path]],
+        topology_file: typing.Optional[typing.Union[str, pathlib.Path]] = None,
+        rtol: float = 0.001) -> bool:
     """Compare multiple trajectory files for equality.
-    
+
     :param trajectory_files: Paths of trajectory file to compare
     :param topology_file: Topology file path
     :param coords_only: Only compare coordinates from trajectory - e.g. not box size
     """
-    ref_traj = mdtraj.load(str(trajectory_files[0]), top=str(topology_file))
+    ref_traj = load_optional_topology(trajectory_files[0], topology_file)
 
     for traj_file in trajectory_files[1:]:
         try:
-            traj = mdtraj.load(str(traj_file), top=str(topology_file))
+            traj = load_optional_topology(traj_file, topology_file)
 
         except ValueError as exc:
             raise ValueError('Topology does not match') from exc
diff --git a/test/test_bondset.py b/test/test_bondset.py
index 08d56298f2f7102107b0f9d427cf33fc3a754cc2..aa34785b9403dc91d49d76a6f71c6762cf3f0ba0 100644
--- a/test/test_bondset.py
+++ b/test/test_bondset.py
@@ -2,6 +2,8 @@ import unittest
 
 import math
 import os
+import pathlib
+import tempfile
 
 from pycgtool.bondset import BondSet
 from pycgtool.frame import Frame
@@ -31,6 +33,9 @@ class DummyBond:
 
 # TODO add setup and teardown to put all files in a tmp directory
 class BondSetTest(unittest.TestCase):
+    base_dir = pathlib.Path(__file__).absolute().parent
+    data_dir = base_dir.joinpath('data')
+
     # Columns are: eqm value, standard fc, MARTINI defaults fc
     invert_test_ref_data = [
         ( 0.220474419132,  72658.18163, 1250),
@@ -104,8 +109,6 @@ class BondSetTest(unittest.TestCase):
         mapping = Mapping('test/data/sugar.map', DummyOptions)
 
         cg_frame = mapping.apply(frame)
-        while frame.next_frame():
-            cg_frame = mapping.apply(frame, cg_frame=cg_frame)
 
         measure.apply(cg_frame)
         measure.boltzmann_invert()
@@ -121,8 +124,6 @@ class BondSetTest(unittest.TestCase):
         mapping = Mapping('test/data/sugar.map', DefaultOptions)
 
         cg_frame = mapping.apply(frame)
-        while frame.next_frame():
-            cg_frame = mapping.apply(frame, cg_frame=cg_frame)
 
         measure.apply(cg_frame)
         measure.boltzmann_invert()
@@ -140,8 +141,6 @@ class BondSetTest(unittest.TestCase):
         mapping = Mapping('test/data/sugar.map', DummyOptions)
 
         cg_frame = mapping.apply(frame)
-        while frame.next_frame():
-            cg_frame = mapping.apply(frame, cg_frame=cg_frame)
 
         measure.apply(cg_frame)
         measure.boltzmann_invert()
@@ -180,17 +179,23 @@ class BondSetTest(unittest.TestCase):
         measure = BondSet('test/data/sugar.bnd', DummyOptions)
         frame = Frame('test/data/sugar.gro', 'test/data/sugar.xtc')
         mapping = Mapping('test/data/sugar.map', DummyOptions)
-        cg_frame = mapping.apply(frame)
 
-        while frame.next_frame():
-            cg_frame = mapping.apply(frame, cg_frame=cg_frame)
+        cg_frame = mapping.apply(frame)
 
         measure.apply(cg_frame)
         measure.boltzmann_invert()
-        measure.write_itp('sugar_out.itp', mapping)
 
-        self.assertTrue(cmp_file_whitespace_float('sugar_out.itp', 'test/data/sugar_out.itp',
-                                                  rtol=0.005, verbose=True))
+        with tempfile.TemporaryDirectory() as tmpdir:
+            tmp_path = pathlib.Path(tmpdir)
+
+            measure.write_itp(tmp_path.joinpath('sugar_out.itp'), mapping)
+
+            self.assertTrue(
+                cmp_file_whitespace_float(
+                    tmp_path.joinpath('sugar_out.itp'),
+                    self.data_dir.joinpath('sugar_out.itp'),
+                    rtol=0.005,
+                    verbose=True))
 
     def test_full_itp_vsites(self):
         """Test full operation to output of .itp file for molecule with vsites."""
@@ -200,17 +205,23 @@ class BondSetTest(unittest.TestCase):
         measure = BondSet('test/data/martini3/naphthalene.bnd', options)
         frame = Frame('test/data/martini3/naphthalene.gro')
         mapping = Mapping('test/data/martini3/naphthalene.map', options)
+
         cg_frame = mapping.apply(frame)
 
         measure.apply(cg_frame)
         measure.boltzmann_invert()
-        measure.write_itp('naphthalene_out.itp', mapping)
 
-        self.assertTrue(
-            cmp_file_whitespace_float('naphthalene_out.itp',
-                                      'test/data/martini3/naphthalene_out.itp',
-                                      rtol=0.005,
-                                      verbose=True))
+        with tempfile.TemporaryDirectory() as tmpdir:
+            tmp_path = pathlib.Path(tmpdir)
+            measure.write_itp(tmp_path.joinpath('naphthalene_out.itp'),
+                              mapping)
+
+            self.assertTrue(
+                cmp_file_whitespace_float(
+                    tmp_path.joinpath('naphthalene_out.itp'),
+                    self.data_dir.joinpath('martini3/naphthalene_out.itp'),
+                    rtol=0.005,
+                    verbose=True))
 
     def test_duplicate_atoms_in_bond(self):
         with self.assertRaises(ValueError):
@@ -220,19 +231,21 @@ class BondSetTest(unittest.TestCase):
         measure = BondSet('test/data/sugar.bnd', DummyOptions)
         frame = Frame('test/data/sugar.gro', 'test/data/sugar.xtc')
         mapping = Mapping('test/data/sugar.map', DummyOptions)
-        cg_frame = mapping.apply(frame)
 
-        while frame.next_frame():
-            cg_frame = mapping.apply(frame, cg_frame=cg_frame)
+        cg_frame = mapping.apply(frame)
 
         measure.apply(cg_frame)
         measure.boltzmann_invert()
         measure.dump_values()
 
+        data_dir = pathlib.Path('test/data')
         filenames = ('ALLA_length.dat', 'ALLA_angle.dat', 'ALLA_dihedral.dat')
         for filename in filenames:
-            self.assertTrue(cmp_file_whitespace_float(os.path.join('test/data', filename), filename,
-                                                      rtol=0.008, verbose=True))
+            self.assertTrue(
+                cmp_file_whitespace_float(data_dir.joinpath(filename),
+                                          filename,
+                                          rtol=0.008,
+                                          verbose=True))
             os.remove(filename)
 
     def test_get_lines_for_bond_dump(self):
diff --git a/test/test_frame.py b/test/test_frame.py
index fd340ddf252a084e54690dd90894bbe1f278fb71..0b073668c2c79affb47cd20c1edf8fa5f981de8e 100644
--- a/test/test_frame.py
+++ b/test/test_frame.py
@@ -1,11 +1,14 @@
 import filecmp
 import os
+import pathlib
 import unittest
+import tempfile
 
 import numpy as np
 
 from pycgtool.frame import NonMatchingSystemError, UnsupportedFormatException
 from pycgtool.frame import Frame
+from pycgtool import util
 
 
 def try_remove(filename) -> None:
@@ -18,6 +21,9 @@ def try_remove(filename) -> None:
 
 # TODO add ITP parsing tests
 class FrameTest(unittest.TestCase):
+    base_dir = pathlib.Path(__file__).absolute().parent
+    data_dir = base_dir.joinpath('data')
+
     def check_reference_topology(self, frame, skip_names=True):
         self.assertEqual(663, frame.natoms)
 
@@ -33,8 +39,8 @@ class FrameTest(unittest.TestCase):
     def check_reference_frame(self, frame):
         self.check_reference_topology(frame)
 
-        atom0_coords = np.array([0.696, 1.330, 1.211])
-        box_vectors = np.array([1.89868, 1.89868, 1.89868])
+        atom0_coords = np.array([[0.696, 1.330, 1.211]])
+        box_vectors = np.array([[1.89868, 1.89868, 1.89868]])
 
         np.testing.assert_allclose(atom0_coords, frame.atom(0).coords)
         np.testing.assert_allclose(box_vectors, frame.unitcell_lengths,
@@ -43,18 +49,38 @@ class FrameTest(unittest.TestCase):
     def check_reference_trajectory(self, frame):
         self.check_reference_topology(frame)
 
-        atom0_coords_array = np.array([[1.176, 1.152, 1.586],
-                                       [1.122, 1.130, 1.534]])
-
-        box_vectors_array = np.array([[1.9052, 1.9052, 1.9052],
-                                      [1.90325272, 1.90325272, 1.90325272]])
+        atom0_coords = np.array([
+            [1.176     , 1.1520001 , 1.5860001 ],
+            [1.1220001 , 1.13      , 1.534     ],
+            [1.0580001 , 1.1620001 , 1.462     ],
+            [0.91600007, 1.276     , 1.5580001 ],
+            [0.73200005, 1.1240001 , 1.286     ],
+            [0.64000005, 1.2160001 , 1.258     ],
+            [0.632     , 1.312     , 1.2520001 ],
+            [0.606     , 1.284     , 1.246     ],
+            [0.582     , 1.312     , 1.1600001 ],
+            [0.68200004, 1.22      , 1.25      ],
+            [0.69600004, 1.33      , 1.21      ],
+        ], dtype=np.float32)  # yapf: disable
+
+        unitcell_lengths = np.array([
+            [1.9052   , 1.9052   , 1.9052   ],
+            [1.9032527, 1.9032527, 1.9032527],
+            [1.9040661, 1.9040661, 1.9040661],
+            [1.896811 , 1.896811 , 1.896811 ],
+            [1.8985983, 1.8985983, 1.8985983],
+            [1.9033976, 1.9033976, 1.9033976],
+            [1.8904614, 1.8904614, 1.8904614],
+            [1.9013108, 1.9013108, 1.9013108],
+            [1.8946321, 1.8946321, 1.8946321],
+            [1.898091 , 1.898091 , 1.898091 ],
+            [1.898684 , 1.898684 , 1.898684 ],
+        ], dtype=np.float32)  # yapf: disable
 
-        for _, (atom0_coords, box_vectors) in enumerate(
-                zip(atom0_coords_array, box_vectors_array)):
-            np.testing.assert_allclose(atom0_coords, frame.atom(0).coords)
-            np.testing.assert_allclose(box_vectors, frame.unitcell_lengths,
-                                       rtol=1e-4)  # PDB files are f9.3
-            frame.next_frame()
+        np.testing.assert_allclose(atom0_coords, frame.atom(0).coords)
+        np.testing.assert_allclose(unitcell_lengths,
+                                   frame.unitcell_lengths,
+                                   rtol=1e-4)  # PDB files are f9.3
 
     def test_frame_create(self):
         Frame()
@@ -85,25 +111,43 @@ class FrameTest(unittest.TestCase):
 
     def test_frame_output_gro(self):
         frame = Frame('test/data/water.gro')
-        frame.save('water-out.gro')
-        self.assertTrue(filecmp.cmp('test/data/water.gro', 'water-out.gro'))
 
-        try_remove('water-out.gro')
+        with tempfile.TemporaryDirectory() as tmpdir:
+            tmp_path = pathlib.Path(tmpdir)
+
+            frame.save(tmp_path.joinpath('water-out.gro'))
+            self.assertTrue(
+                filecmp.cmp(self.data_dir.joinpath('water.gro'),
+                            tmp_path.joinpath('water-out.gro')))
 
     def test_frame_read_xtc_numframes(self):
         frame = Frame('test/data/water.gro', 'test/data/water.xtc')
-        self.assertEqual(10, frame.numframes)
+        self.assertEqual(11, frame.n_frames)
 
     def test_frame_read_xtc(self):
         frame = Frame('test/data/water.gro', 'test/data/water.xtc')
 
         self.check_reference_trajectory(frame)
 
-    def test_frame_write_xtc(self):
-        try_remove('water_test2.xtc')
+    def test_frame_accept_path(self):
+        """Test that Frame accepts :class:`pathlib.Path` arguments."""
+        _ = Frame(self.data_dir.joinpath('water.gro'),
+                  self.data_dir.joinpath('water.xtc'))
 
-        frame = Frame('test/data/water.gro', 'test/data/water.xtc')
-        frame.save('water_test2.xtc')
+    def test_frame_write_xtc(self):
+        """Test that :class:`Frame` can save a trajectory file."""
+        frame = Frame(self.data_dir.joinpath('water.gro'),
+                      self.data_dir.joinpath('water.xtc'))
+
+        with tempfile.TemporaryDirectory() as tmpdir:
+            tmp_path = pathlib.Path(tmpdir)
+
+            frame.save(tmp_path.joinpath('water.xtc'))
+            self.assertTrue(
+                util.compare_trajectories(
+                    self.data_dir.joinpath('water.xtc'),
+                    tmp_path.joinpath('water.xtc'),
+                    topology_file=self.data_dir.joinpath('water.gro')))
 
     def test_raise_nonmatching_system_mdtraj(self):
         with self.assertRaises(NonMatchingSystemError):
diff --git a/test/test_mapping.py b/test/test_mapping.py
index fa1d6caaf2182fb7c1192ede1d1ed39f4e0221d6..d6f987bab9759a17efc19e35e02823ee9a0dd404 100644
--- a/test/test_mapping.py
+++ b/test/test_mapping.py
@@ -110,7 +110,7 @@ class MappingTest(unittest.TestCase):
         mapping = Mapping("test/data/two.map", DummyOptions)
         cg_frame = mapping.apply(frame)
 
-        np.testing.assert_allclose(np.array([1.5, 1.5, 1.5]),
+        np.testing.assert_allclose(np.array([[1.5, 1.5, 1.5]]),
                                    cg_frame.residue(0).atom(0).coords)
 
     def test_virtual_mapping_weights_geom(self):
@@ -119,7 +119,7 @@ class MappingTest(unittest.TestCase):
         mapping = Mapping("test/data/martini3/four.map", DummyOptions)
         cg_frame = mapping.apply(frame)
 
-        np.testing.assert_allclose(np.array([2.5, 2.5, 2.5]), 
+        np.testing.assert_allclose(np.array([[2.5, 2.5, 2.5]]), 
                                    cg_frame.residue(0).atom(2).coords)
 
     def test_mapping_weights_mass(self):
@@ -130,7 +130,7 @@ class MappingTest(unittest.TestCase):
         mapping = Mapping("test/data/two.map", options, itp_filename="test/data/two.itp")
         cg_frame = mapping.apply(frame)
 
-        np.testing.assert_allclose(np.array([2., 2., 2.]),
+        np.testing.assert_allclose(np.array([[2., 2., 2.]]),
                                    cg_frame.residue(0).atom(0).coords)
 
     def test_virtual_mapping_weights_mass(self):
@@ -141,7 +141,7 @@ class MappingTest(unittest.TestCase):
         mapping = Mapping("test/data/martini3/four.map", options, itp_filename="test/data/martini3/four.itp")
         cg_frame = mapping.apply(frame)
 
-        np.testing.assert_allclose(np.array([3.0, 3.0, 3.0]),
+        np.testing.assert_allclose(np.array([[3.0, 3.0, 3.0]]),
                                    cg_frame.residue(0).atom(2).coords)
 
     def test_mapping_weights_guess_mass(self):
@@ -152,7 +152,7 @@ class MappingTest(unittest.TestCase):
         mapping = Mapping("test/data/two.map", options)
         cg_frame = mapping.apply(frame)
 
-        np.testing.assert_allclose(np.array([1.922575,  1.922575,  1.922575], dtype=np.float32),
+        np.testing.assert_allclose(np.array([[1.922575,  1.922575,  1.922575]], dtype=np.float32),
                                    cg_frame.residue(0).atom(0).coords, rtol=0.01)
 
     def test_virtual_mapping_weights_guess_mass(self):
@@ -163,7 +163,7 @@ class MappingTest(unittest.TestCase):
         mapping = Mapping("test/data/martini3/four.map", options)
         cg_frame = mapping.apply(frame)
 
-        np.testing.assert_allclose(np.array([2.83337, 2.83337, 2.83337], dtype=np.float32),
+        np.testing.assert_allclose(np.array([[2.83337, 2.83337, 2.83337]], dtype=np.float32),
                                    cg_frame.residue(0).atom(2).coords, rtol=0.01)
 
     def test_mapping_weights_first(self):
@@ -174,7 +174,7 @@ class MappingTest(unittest.TestCase):
         mapping = Mapping("test/data/two.map", options, itp_filename="test/data/two.itp")
         cg_frame = mapping.apply(frame)
 
-        np.testing.assert_allclose(np.array([1., 1., 1.]),
+        np.testing.assert_allclose(np.array([[1., 1., 1.]]),
                                    cg_frame.residue(0).atom(0).coords)
 
     def test_mapping_itp_multi(self):
diff --git a/test/test_pycgtool.py b/test/test_pycgtool.py
index 820449650cbc7669f61cf25af9c6eef862f9e27d..a2b9813a5c6ee9e4132a33aeaa47a2e6f758b75a 100644
--- a/test/test_pycgtool.py
+++ b/test/test_pycgtool.py
@@ -91,7 +91,6 @@ class PycgtoolTest(unittest.TestCase):
                 util.compare_trajectories(
                     self.data_dir.joinpath('sugar_out.gro'),
                     tmp_path.joinpath('out.gro'),
-                    topology_file=self.data_dir.joinpath('sugar_out.gro'),
                     rtol=0.001))
 
     def test_full_no_traj(self):
@@ -109,7 +108,6 @@ class PycgtoolTest(unittest.TestCase):
                 util.compare_trajectories(
                     self.data_dir.joinpath('sugar_out.gro'),
                     tmp_path.joinpath('out.gro'),
-                    topology_file=self.data_dir.joinpath('sugar_out.gro'),
                     rtol=0.001))
 
     def test_measure_only(self):
diff --git a/test/test_util.py b/test/test_util.py
index 8c3b1ef3d51bc1a4cf0cc5c3bbff4c75cd0d844c..7e5cbc7350185f67b7ba2ebc7ba95a62605df48c 100644
--- a/test/test_util.py
+++ b/test/test_util.py
@@ -172,16 +172,14 @@ class CompareTrajectoryTest(unittest.TestCase):
     def test_compare_trajectory_single(self):
         self.assertTrue(util.compare_trajectories(
             'test/data/sugar.gro',
-            'test/data/sugar.gro',
-            topology_file='test/data/sugar.gro'
+            'test/data/sugar.gro'
         ))
 
     def test_compare_trajectory_single_false(self):
         with self.assertRaises(ValueError):
             util.compare_trajectories(
                 'test/data/sugar.gro',
-                'test/data/water.gro',
-                topology_file='test/data/sugar.gro'
+                'test/data/water.gro'
             )
 
     def test_compare_trajectory(self):