diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index e23f6c2..0e72460 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -35,10 +35,49 @@ jobs:
           pip install .
       - name: Test
         run: |
-          python -m pytest
+          python -m pytest -m "not lammps"
       - name: Install optional dependencies
         run: |
           pip install -r tests/requirements-optin.txt
       - name: Test (optin)
         run: |
-          python -m pytest
+          python -m pytest -m "not lammps"
+  test-lammps:
+    name: unit test lammps [python-${{ matrix.python-version }}, lammps-${{ matrix.lammps-version }}]
+    needs: test
+    runs-on: ubuntu-latest
+    defaults:
+      run:
+        shell: bash -el {0}
+    strategy:
+      fail-fast: false
+      matrix:
+        python-version: ["3.9", "3.10", "3.11", "3.12", "3.13"]
+        lammps-version: ["2023.08.02", "2024.08.29"]
+        exclude:
+          - python-version: "3.12"
+            lammps-version: "2023.08.02"
+          - python-version: "3.13"
+            lammps-version: "2023.08.02"
+          - python-version: "3.13"
+            lammps-version: "2024.08.29"
+    steps:
+      - name: Checkout code
+        uses: actions/checkout@v4
+      - name: Setup Python ${{ matrix.python-version }}
+        uses: conda-incubator/setup-miniconda@v3
+        with:
+          channels: conda-forge,defaults
+          channel-priority: strict
+          python-version: ${{ matrix.python-version }}
+      - name: Install LAMMPS
+        run: |
+          conda install "lammps=${{ matrix.lammps-version }}"
+      - name: Install
+        run: |
+          python -m pip install pip --upgrade pip
+          pip install -r tests/requirements.txt -r tests/requirements-optin.txt
+          pip install .
+      - name: Test
+        run: |
+          python -m pytest -m lammps
diff --git a/pyproject.toml b/pyproject.toml
index 129d004..c03f183 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -11,3 +11,6 @@ target-version = ["py39"]
 
 [tool.isort]
 profile = "black"
+
+[tool.pytest.ini_options]
+markers = ["lammps"]
diff --git a/tests/conftest.py b/tests/conftest.py
index 084fc56..b7118d4 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -6,12 +6,12 @@
 
 @pytest.fixture
 def orthorhombic():
-    return lammpsio.Box([-5.0, -10.0, 0.0], [1.0, 10.0, 8.0])
+    return lammpsio.Box([-5.0, -10.0, -1.0], [2.0, 10.0, 8.0])
 
 
 @pytest.fixture
 def triclinic():
-    return lammpsio.Box([-5.0, -10.0, 0.0], [1.0, 10.0, 8.0], [1.0, -2.0, 0.5])
+    return lammpsio.Box([-5.0, -10.0, -1.0], [2.0, 10.0, 8.0], [2, -2.0, 0.5])
 
 
 @pytest.fixture(params=[lf("orthorhombic"), lf("triclinic")])
diff --git a/tests/test_box.py b/tests/test_box.py
index 111065f..aaa612c 100644
--- a/tests/test_box.py
+++ b/tests/test_box.py
@@ -4,8 +4,8 @@
 
 def test_orthorhombic(orthorhombic):
     box = orthorhombic
-    assert numpy.allclose(box.low, [-5, -10, 0])
-    assert numpy.allclose(box.high, [1, 10, 8])
+    assert numpy.allclose(box.low, [-5, -10, -1])
+    assert numpy.allclose(box.high, [2, 10, 8])
     assert box.tilt is None
 
     box.low = [0, 0, 0]
@@ -21,9 +21,9 @@ def test_orthorhombic(orthorhombic):
 
 def test_triclinic(triclinic):
     box = triclinic
-    assert numpy.allclose(box.low, [-5, -10, 0])
-    assert numpy.allclose(box.high, [1, 10, 8])
-    assert numpy.allclose(box.tilt, [1.0, -2.0, 0.5])
+    assert numpy.allclose(box.low, [-5, -10, -1])
+    assert numpy.allclose(box.high, [2, 10, 8])
+    assert numpy.allclose(box.tilt, [2.0, -2.0, 0.5])
 
     box.low = [0, 0, 0]
     assert numpy.allclose(box.low, [0, 0, 0])
diff --git a/tests/test_data_file.py b/tests/test_data_file.py
index 539f336..be4b94a 100644
--- a/tests/test_data_file.py
+++ b/tests/test_data_file.py
@@ -3,6 +3,13 @@
 
 import lammpsio
 
+try:
+    import lammps
+
+    has_lammps = True
+except ImportError:
+    has_lammps = False
+
 
 @pytest.mark.parametrize("shuffle_ids", [False, True])
 @pytest.mark.parametrize("atom_style", ["atomic", "molecular", "charge", "full"])
@@ -183,7 +190,7 @@ def test_data_file_topology(snap_8, tmp_path, shuffle_ids):
     assert filename.exists
 
     snap_2 = data.read()
-    # test none topology features of a snapshot with topology
+    # test non-topology features of a snapshot with topology
     assert snap_2.N == snap_8.N
     if shuffle_ids:
         assert snap_2.has_id()
@@ -249,3 +256,304 @@ def test_data_file_topology(snap_8, tmp_path, shuffle_ids):
     assert numpy.allclose(snap_2.impropers.typeid, snap_8.impropers.typeid)
     assert snap_2.impropers.has_members()
     assert numpy.allclose(snap_2.impropers.members, snap_8.impropers.members)
+
+
+@pytest.mark.lammps
+@pytest.mark.skipif(not has_lammps, reason="lammps not installed")
+@pytest.mark.parametrize("shuffle_ids", [False, True])
+@pytest.mark.parametrize("atom_style", ["atomic", "molecular", "charge", "full"])
+def test_data_file_min_lammps(tmp_path, snap, atom_style, shuffle_ids):
+    if shuffle_ids:
+        snap.id = [3, 1, 2]
+
+    # write the data file with default values using lammpsio
+    filename = tmp_path / "atoms.data"
+    assert not snap.has_image()
+    assert not snap.has_velocity()
+    lammpsio.DataFile.create(filename, snap, atom_style)
+    assert filename.exists
+
+    # read it in LAMMPS and write it out
+    lmps = lammps.lammps(cmdargs=["-log", "none", "-nocite"])
+    cmds = [f"atom_style {atom_style}"]
+    cmds += [f"read_data {filename}"]
+    cmds += ["mass 1 1.0"]
+    cmds += [f"write_data {filename} nocoeff nofix"]
+    lmps.commands_list(cmds)
+    lmps.close()
+
+    # read it back in and check
+    snap_2 = lammpsio.DataFile(filename).read()
+    assert snap_2.N == snap.N
+    assert snap_2.step is None
+    assert numpy.allclose(snap_2.box.low, snap.box.low)
+    assert numpy.allclose(snap_2.box.high, snap.box.high)
+    if snap.box.tilt is not None:
+        assert numpy.allclose(snap_2.box.tilt, snap.box.tilt)
+    else:
+        assert snap_2.box.tilt is None
+    if shuffle_ids:
+        assert snap_2.has_id()
+        assert numpy.allclose(snap_2.id, snap.id)
+    else:
+        assert not snap_2.has_id()
+    assert snap_2.has_position()
+    assert numpy.allclose(snap_2.position, 0)
+    assert snap_2.has_typeid()
+    assert numpy.allclose(snap_2.typeid, 1)
+
+    if atom_style in ("molecular", "full"):
+        assert snap_2.has_molecule()
+        assert numpy.allclose(snap_2.molecule, 0)
+    else:
+        assert not snap_2.has_molecule()
+    if atom_style in ("charge", "full"):
+        assert snap_2.has_charge()
+        assert numpy.allclose(snap_2.charge, 0)
+    else:
+        assert not snap_2.has_charge()
+
+
+@pytest.mark.lammps
+@pytest.mark.skipif(not has_lammps, reason="lammps not installed")
+@pytest.mark.parametrize("shuffle_ids", [False, True])
+@pytest.mark.parametrize("set_style", [True, False])
+@pytest.mark.parametrize("atom_style", ["atomic", "molecular", "charge", "full"])
+def test_data_file_all_lammps(snap, atom_style, set_style, shuffle_ids, tmp_path):
+    if shuffle_ids:
+        snap.id = [3, 1, 2]
+
+    # write the data file with nondefault values using lammpsio
+    snap.position = [[0.1, 0.2, 0.3], [-0.4, -0.5, -0.6], [0.7, 0.8, 0.9]]
+    snap.image = [[1, 2, 3], [-4, -5, -6], [7, 8, 9]]
+    snap.velocity = [[-3, -2, -1], [6, 5, 4], [9, 8, 7]]
+    snap.typeid = [2, 1, 2]
+    snap.mass = [3, 2, 3]
+    if atom_style in ("molecular", "full"):
+        snap.molecule = [2, 0, 1]
+    if atom_style in ("charge", "full"):
+        snap.charge = [-1, 0, 1]
+
+    filename = tmp_path / "atoms.data"
+
+    # create the data file using lammpsio
+    lammpsio.DataFile.create(filename, snap, atom_style if set_style else None)
+    assert filename.exists
+
+    lmps = lammps.lammps(cmdargs=["-log", "none", "-nocite"])
+    cmds = [f"atom_style {atom_style}"]
+    cmds += [f"read_data {filename}"]
+    cmds += [f"write_data {filename}"]
+    lmps.commands_list(cmds)
+    lmps.close()
+
+    # read it back in and check
+    snap_2 = lammpsio.DataFile(filename).read()
+    assert snap_2.N == snap.N
+    assert snap_2.step is None
+    assert numpy.allclose(snap_2.box.low, snap.box.low)
+    assert numpy.allclose(snap_2.box.high, snap.box.high)
+    if snap.box.tilt is not None:
+        assert numpy.allclose(snap_2.box.tilt, snap.box.tilt)
+    else:
+        assert snap_2.box.tilt is None
+    if shuffle_ids:
+        assert snap_2.has_id()
+        assert numpy.allclose(snap_2.id, snap.id)
+    else:
+        assert not snap_2.has_id()
+    assert snap_2.has_position()
+    assert numpy.allclose(snap_2.position, snap.position)
+    assert snap_2.has_velocity()
+    assert numpy.allclose(snap_2.velocity, snap.velocity)
+    assert snap_2.has_image()
+    assert numpy.allclose(snap_2.image, snap.image)
+    assert snap_2.has_typeid()
+    assert numpy.allclose(snap_2.typeid, snap.typeid)
+    assert snap_2.has_mass()
+    assert numpy.allclose(snap_2.mass, snap.mass)
+    if atom_style in ("molecular", "full"):
+        assert snap_2.has_molecule()
+        assert numpy.allclose(snap_2.molecule, snap.molecule)
+    else:
+        assert not snap_2.has_molecule()
+    if atom_style in ("charge", "full"):
+        assert snap_2.has_charge()
+        assert numpy.allclose(snap_2.charge, snap.charge)
+    else:
+        assert not snap_2.has_charge()
+
+
+@pytest.mark.lammps
+@pytest.mark.skipif(not has_lammps, reason="lammps not installed")
+@pytest.mark.parametrize("shuffle_ids", [False, True])
+def test_data_file_topology_lammps(snap_8, tmp_path, shuffle_ids):
+    # set ids to be assigned
+    # Note: LAMMPS reorders the ids of topology objects when ids are shuffled
+    # We only check that topology objects typeid and members survive the round
+    # trip, but in future, we could check the IDs if this behavior were changed.
+    if shuffle_ids:
+        particle_id = [1, 5, 2, 6, 3, 7, 4, 8]
+        bond_id = [1, 4, 2, 5, 3, 6]
+        angle_id = [1, 3, 2, 4]
+        dihedral_id = [2, 1]
+        improper_id = [2, 1]
+    else:
+        particle_id = [1, 2, 3, 4, 5, 6, 7, 8]
+        bond_id = [1, 2, 3, 4, 5, 6]
+        angle_id = [1, 2, 3, 4]
+        dihedral_id = [1, 2]
+        improper_id = [1, 2]
+
+    # particle information
+    snap_8.id = particle_id
+    snap_8.typeid = [1, 1, 1, 1, 2, 2, 2, 2]
+    snap_8.position = [
+        [0, 0, 0],
+        [0.1, 0.1, 0.1],
+        [0.2, 0.2, 0.2],
+        [0.3, 0.3, 0.3],
+        [1, 1, 1],
+        [1.1, 1.1, 1.1],
+        [1.2, 1.2, 1.2],
+        [1.3, 1.3, 1.3],
+    ]
+    snap_8.mass = [1, 1, 1, 1, 2, 2, 2, 2]
+
+    # bond information
+    snap_8.bonds = lammpsio.topology.Bonds(N=6, num_types=2)
+    snap_8.bonds.id = bond_id
+    snap_8.bonds.typeid = [1, 2, 1, 2, 1, 2]
+    snap_8.bonds.members = [
+        [1, 2],
+        [2, 3],
+        [3, 4],
+        [5, 6],
+        [6, 7],
+        [7, 8],
+    ]
+
+    # angle information
+    snap_8.angles = lammpsio.topology.Angles(N=4, num_types=2)
+    snap_8.angles.id = angle_id
+    snap_8.angles.typeid = [1, 2, 2, 1]
+    snap_8.angles.members = [
+        [1, 2, 3],
+        [2, 3, 4],
+        [5, 6, 7],
+        [6, 7, 8],
+    ]
+
+    # dihedral information
+    snap_8.dihedrals = lammpsio.topology.Dihedrals(N=2, num_types=2)
+    snap_8.dihedrals.id = dihedral_id
+    snap_8.dihedrals.typeid = [1, 2]
+    snap_8.dihedrals.members = [
+        [1, 2, 3, 4],
+        [5, 6, 7, 8],
+    ]
+
+    # improper information
+    snap_8.impropers = lammpsio.topology.Impropers(N=2, num_types=2)
+    snap_8.impropers.id = improper_id
+    snap_8.impropers.typeid = [1, 2]
+    snap_8.impropers.members = [
+        [1, 2, 3, 4],
+        [5, 6, 7, 8],
+    ]
+
+    filename = tmp_path / "atoms.data"
+    # create the data file using lammpsio
+    lammpsio.DataFile.create(filename, snap_8)
+
+    # read it in LAMMPS and write it out
+    lmps = lammps.lammps(cmdargs=["-log", "none", "-nocite"])
+    cmds = ["atom_style molecular"]
+    cmds += [f"read_data {filename}"]
+    cmds += [f"write_data {filename}"]
+    lmps.commands_list(cmds)
+    lmps.close()
+
+    # read it back in and check
+    snap_2 = lammpsio.DataFile(filename).read()
+    # test non-topology features of a snapshot with topology
+    assert snap_2.N == snap_8.N
+    if shuffle_ids:
+        assert snap_2.has_id()
+        assert numpy.allclose(snap_2.id, snap_8.id)
+    else:
+        assert not snap_2.has_id()
+
+    assert numpy.allclose(snap_2.id, snap_8.id)
+    assert snap_2.has_typeid()
+    assert numpy.allclose(snap_2.typeid, snap_8.typeid)
+    assert snap_2.has_position()
+    assert numpy.allclose(snap_2.position, snap_8.position)
+    assert snap_2.has_mass()
+    assert numpy.allclose(snap_2.mass, snap_8.mass)
+
+    # test bonds
+    assert snap_2.bonds.N == snap_8.bonds.N
+    assert snap_2.bonds.has_typeid()
+    assert snap_2.bonds.has_members()
+
+    # test that types and bond members survive round trip
+    # Note: This sorting only works because each member only appears in the first
+    # column one time.
+    snap_8_bond_data = numpy.column_stack((snap_8.bonds.members, snap_8.bonds.typeid))
+    snap_2_bond_data = numpy.column_stack((snap_2.bonds.members, snap_2.bonds.typeid))
+    assert numpy.allclose(
+        snap_2_bond_data[snap_2_bond_data[:, 0].argsort()],
+        snap_8_bond_data[snap_8_bond_data[:, 0].argsort()],
+    )
+
+    # test angles
+    assert snap_2.angles.N == snap_8.angles.N
+    assert snap_2.angles.has_typeid()
+    assert snap_2.angles.has_members()
+
+    # test that types and angle members survive round trip
+    snap_8_angle_data = numpy.column_stack(
+        (snap_8.angles.members, snap_8.angles.typeid)
+    )
+    snap_2_angle_data = numpy.column_stack(
+        (snap_2.angles.members, snap_2.angles.typeid)
+    )
+    assert numpy.allclose(
+        snap_2_angle_data[snap_2_angle_data[:, 0].argsort()],
+        snap_8_angle_data[snap_8_angle_data[:, 0].argsort()],
+    )
+
+    # test dihedrals
+    assert snap_2.dihedrals.N == snap_8.dihedrals.N
+    assert snap_2.dihedrals.has_typeid()
+    assert snap_2.dihedrals.has_members()
+
+    # test that types and dihedral members survive round trip
+    snap_8_dihedral_data = numpy.column_stack(
+        (snap_8.dihedrals.members, snap_8.dihedrals.typeid)
+    )
+    snap_2_dihedral_data = numpy.column_stack(
+        (snap_2.dihedrals.members, snap_2.dihedrals.typeid)
+    )
+    assert numpy.allclose(
+        snap_2_dihedral_data[snap_2_dihedral_data[:, 0].argsort()],
+        snap_8_dihedral_data[snap_8_dihedral_data[:, 0].argsort()],
+    )
+
+    # test impropers
+    assert snap_2.impropers.N == snap_8.impropers.N
+    assert snap_2.impropers.has_typeid()
+    assert snap_2.impropers.has_members()
+
+    # test that types and improper members survive round trip
+    snap_8_improper_data = numpy.column_stack(
+        (snap_8.impropers.members, snap_8.impropers.typeid)
+    )
+    snap_2_improper_data = numpy.column_stack(
+        (snap_2.impropers.members, snap_2.impropers.typeid)
+    )
+    assert numpy.allclose(
+        snap_2_improper_data[snap_2_improper_data[:, 0].argsort()],
+        snap_8_improper_data[snap_8_improper_data[:, 0].argsort()],
+    )
diff --git a/tests/test_dump_file.py b/tests/test_dump_file.py
index d236d81..5d0a161 100644
--- a/tests/test_dump_file.py
+++ b/tests/test_dump_file.py
@@ -12,6 +12,13 @@
 except ModuleNotFoundError:
     has_pyzstd = False
 
+try:
+    import lammps
+
+    has_lammps = True
+except ImportError:
+    has_lammps = False
+
 
 @pytest.mark.parametrize("sort_ids", [False, True])
 @pytest.mark.parametrize("shuffle_ids", [False, True])
@@ -318,3 +325,225 @@ def test_faulty_dump_schema(tmp_path, schema):
     with pytest.raises(IOError):
         for snap in traj:
             pass
+
+
+@pytest.mark.lammps
+@pytest.mark.parametrize("sort_ids", [False, True])
+@pytest.mark.parametrize("shuffle_ids", [False, True])
+@pytest.mark.parametrize("compression_extension", ["", ".gz", ".zst"])
+@pytest.mark.skipif(not has_lammps, reason="lammps not installed")
+def test_dump_file_min_lammps(
+    snap, compression_extension, shuffle_ids, sort_ids, tmp_path
+):
+    if not has_pyzstd and compression_extension == ".zst":
+        pytest.skip("pyzstd not installed")
+
+    # create file with 2 snapshots with defaults, changing N & step
+    snap.mass = [1, 1, 1]
+    snap.step = 0
+    snap_2 = lammpsio.Snapshot(snap.N + 2, snap.box, snap.step + 1)
+    snap_2.mass = [1, 1, 1, 1, 1]
+    snaps = [snap, snap_2]
+    if shuffle_ids:
+        snap.id = snap.id[::-1]
+        # LAMMPS create_atoms gives atom smallest available id
+        snap_2.id = numpy.append(snap.id, [4, 5])
+
+    # create file with first snapshot
+    filename_data = tmp_path / "atoms.data"
+    lammpsio.DataFile.create(filename_data, snap)
+    assert filename_data.exists
+
+    # create dump with 2 snapshots in LAMMPS
+    filename = tmp_path / f"atoms.lammpstrj{compression_extension}"
+    lmps = lammps.lammps(cmdargs=["-log", "none", "-nocite"])
+
+    cmds = []
+    cmds += ["units lj"]
+    cmds += ["atom_style atomic"]
+    cmds += ["boundary p p p"]
+    cmds += [f"read_data {filename_data}"]
+    cmds += ["timestep 1.0"]
+
+    # Dump setup to capture both frames
+    cmds += [f"dump dmp all custom 1 {filename} id x y z"]
+    cmds += ["dump_modify dmp first yes"]
+
+    # Dump initial state frame 0
+    cmds += ["run 0"]
+
+    # Add two atoms of type 1 at 0, 0, 0
+    cmds += ["create_atoms 1 single 0.0 0.0 0.0"]
+    cmds += ["create_atoms 1 single 0.0 0.0 0.0"]
+
+    # dump frame 1
+    cmds += ["run 1"]
+
+    lmps.commands_list(cmds)
+    lmps.close()
+
+    assert filename.exists
+
+    # read it back in and check snapshots
+    f = lammpsio.DumpFile(filename, sort_ids=sort_ids)
+    assert filename.exists
+    assert len(f) == 2
+    for read_snap, snap in zip(f, snaps):
+        assert read_snap.N == snap.N
+        assert read_snap.step == snap.step
+        assert numpy.allclose(read_snap.box.low, snap.box.low)
+        assert numpy.allclose(read_snap.box.high, snap.box.high)
+        if snap.box.tilt is not None:
+            assert numpy.allclose(read_snap.box.tilt, snap.box.tilt)
+        else:
+            assert read_snap.box.tilt is None
+        if shuffle_ids:
+            assert read_snap.has_id()
+            if sort_ids:
+                assert numpy.allclose(read_snap.id, numpy.arange(1, snap.N + 1))
+            else:
+                if snap.N == 3:
+                    assert numpy.allclose(
+                        read_snap.id, numpy.arange(1, snap.N + 1)[::-1]
+                    )
+                elif snap.N == 5:
+                    assert numpy.allclose(read_snap.id, [3, 2, 1, 4, 5])
+        else:
+            assert not read_snap.has_id()
+        assert read_snap.has_position()
+        assert numpy.allclose(read_snap.position, 0)
+        assert not read_snap.has_image()
+        assert not read_snap.has_velocity()
+        assert not read_snap.has_typeid()
+        assert numpy.allclose(read_snap.mass, 1)
+        assert not read_snap.has_molecule()
+        assert not read_snap.has_charge()
+
+
+@pytest.mark.lammps
+@pytest.mark.parametrize("sort_ids", [False, True])
+@pytest.mark.parametrize("shuffle_ids", [False, True])
+@pytest.mark.parametrize("compression_extension", ["", ".gz", ".zst"])
+@pytest.mark.skipif(not has_lammps, reason="lammps not installed")
+def test_dump_file_all_lammps(  #
+    snap, compression_extension, shuffle_ids, sort_ids, tmp_path
+):
+    if not has_pyzstd and compression_extension == ".zst":
+        pytest.skip("pyzstd not installed")
+
+    snap.step = 0
+    snap.position = [[0.1, 0.2, 0.3], [-0.4, -0.5, -0.6], [0.7, 0.8, 0.9]]
+    snap.image = [[1, 2, 3], [-4, -5, -6], [7, 8, 9]]
+    snap.velocity = [[-3, -2, -1], [6, 5, 4], [9, 8, 7]]
+    snap.typeid = [2, 1, 2]
+    snap.mass = [3, 2, 3]
+    snap.molecule = [1, 2, 3]
+    snap.charge = [-1, 0, 1]
+
+    snap_2 = lammpsio.Snapshot(snap.N, snap.box, snap.step + 1)
+    snap_2.step = 1
+    snap_2.position = snap.position[::-1]
+    snap_2.image = snap.image[::-1]
+    snap_2.velocity = snap.velocity[::-1]
+    snap_2.typeid = snap.typeid[::-1]
+    snap_2.mass = snap.mass[::-1]
+    snap_2.molecule = snap.molecule[::-1]
+    snap_2.charge = snap.charge[::-1]
+    snaps = [snap, snap_2]
+
+    if shuffle_ids:
+        for s in snaps:
+            s.id = s.id[::-1]
+        if sort_ids:
+            order = numpy.arange(snap.N)[::-1]
+        else:
+            order = numpy.arange(snap.N)
+    else:
+        order = numpy.arange(snap.N)
+
+    # create file with first snapshot
+    filename_data = tmp_path / "atoms.data"
+    lammpsio.DataFile.create(filename_data, snap)
+    assert filename_data.exists
+
+    # create dump with 2 snapshots in LAMMPS
+    filename = tmp_path / f"atoms.lammpstrj{compression_extension}"
+    lmps = lammps.lammps(cmdargs=["-log", "none", "-nocite"])
+
+    cmds = []
+    cmds += ["units lj"]
+    cmds += ["atom_style full"]
+    cmds += ["boundary p p p"]
+    cmds += [f"read_data {filename_data}"]
+    cmds += ["timestep 1.0"]
+
+    # Dump setup to capture both frames
+    cmds += [
+        f"dump dmp all custom 1 {filename} id type x y z ix iy iz vx vy vz mass mol q"
+    ]
+    cmds += ["dump_modify dmp first yes"]
+
+    # Dump initial state (frame 0)
+    cmds += ["run 0"]
+
+    # Manually reverse all atom data to match snap_2
+    for i in range(snap_2.N):
+        cmds += [
+            f"set atom {snap_2.id[i]} x {snap_2.position[i, 0]} "
+            f"y {snap_2.position[i, 1]} z {snap_2.position[i, 2]}"
+        ]
+        cmds += [
+            f"set atom {snap_2.id[i]} vx {snap_2.velocity[i, 0]} "
+            f"vy {snap_2.velocity[i, 1]} vz {snap_2.velocity[i, 2]}"
+        ]
+        cmds += [f"set atom {snap_2.id[i]} mol {snap_2.molecule[i]}"]
+        cmds += [f"set atom {snap_2.id[i]} type {snap_2.typeid[i]}"]
+        cmds += [f"set atom {snap_2.id[i]} charge {snap_2.charge[i]}"]
+        cmds += [
+            f"set atom {snap_2.id[i]} image {snap_2.image[i, 0]} "
+            f"{snap_2.image[i, 1]} {snap_2.image[i, 2]}"
+        ]
+
+    # Dump reversed state (frame 1)
+    cmds += ["run 1"]
+
+    lmps.commands_list(cmds)
+    lmps.close()
+
+    assert filename.exists
+
+    # read it back in and check snapshots
+    f = lammpsio.DumpFile(filename, sort_ids=sort_ids)
+    assert filename.exists
+    assert len(f) == 2
+    for read_snap, snap in zip(f, snaps):
+        assert read_snap.N == snap.N
+        assert read_snap.step == snap.step
+        assert numpy.allclose(read_snap.box.low, snap.box.low)
+        assert numpy.allclose(read_snap.box.high, snap.box.high)
+        if snap.box.tilt is not None:
+            assert numpy.allclose(read_snap.box.tilt, snap.box.tilt)
+        else:
+            assert read_snap.box.tilt is None
+        if shuffle_ids:
+            assert read_snap.has_id()
+            if sort_ids:
+                assert numpy.allclose(read_snap.id, numpy.arange(1, snap.N + 1))
+            else:
+                assert numpy.allclose(read_snap.id, numpy.arange(1, snap.N + 1)[::-1])
+        else:
+            assert not read_snap.has_id()
+        assert read_snap.has_position()
+        assert numpy.allclose(read_snap.position, snap.position[order])
+        assert read_snap.has_image()
+        assert numpy.allclose(read_snap.image, snap.image[order])
+        assert read_snap.has_velocity()
+        assert numpy.allclose(read_snap.velocity, snap.velocity[order])
+        assert read_snap.has_typeid()
+        assert numpy.allclose(read_snap.typeid, snap.typeid[order])
+        assert read_snap.has_mass()
+        assert numpy.allclose(read_snap.mass, snap.mass[order])
+        assert read_snap.has_molecule()
+        assert numpy.allclose(read_snap.molecule, snap.molecule[order])
+        assert read_snap.has_charge()
+        assert numpy.allclose(read_snap.charge, snap.charge[order])