diff --git a/.ci_support/environment.yml b/.ci_support/environment.yml
index 09d356591..89feaec6b 100644
--- a/.ci_support/environment.yml
+++ b/.ci_support/environment.yml
@@ -10,10 +10,10 @@ dependencies:
 - numpy =1.26.3
 - phonopy =2.21.0
 - plotly =5.18.0
-- pymatgen =2023.12.18
+- pymatgen =2024.1.27
 - pyscal =2.10.18
-- pyxtal =0.6.1
+- pyxtal =0.6.2
 - scikit-learn =1.4.0
-- scipy =1.11.4
-- spglib =2.2.0
+- scipy =1.12.0
+- spglib =2.3.0
 - sqsgenerator =0.3
diff --git a/.ci_support/environment_mini.yml b/.ci_support/environment_mini.yml
index ad9d0d49f..ff2cb899f 100644
--- a/.ci_support/environment_mini.yml
+++ b/.ci_support/environment_mini.yml
@@ -4,4 +4,4 @@ dependencies:
 - ase =3.22.1
 - coverage
 - numpy =1.26.3
-- scipy =1.11.4
+- scipy =1.12.0
diff --git a/pyproject.toml b/pyproject.toml
index 23bd4a720..8d97aa129 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -26,7 +26,7 @@ classifiers = [
 dependencies = [
     "ase==3.22.1",
     "numpy==1.26.3",
-    "scipy==1.11.4",
+    "scipy==1.12.0",
 ]
 dynamic = ["version"]
 
@@ -38,23 +38,23 @@ Repository = "https://github.com/pyiron/structuretoolkit"
 [project.optional-dependencies]
 grainboundary = [
     "aimsgb==1.1.0",
-    "pymatgen==2023.12.18",
+    "pymatgen==2024.1.27",
 ]
 pyscal = ["pyscal2==2.10.18"]
 nglview = ["nglview==3.1.1"]
 matplotlib = ["matplotlib==3.8.2"]
 plotly = ["plotly==5.18.0"]
 clusters = ["scikit-learn==1.4.0"]
-symmetry = ["spglib==2.2.0"]
+symmetry = ["spglib==2.3.0"]
 surface = [
-    "spglib==2.2.0",
-    "pymatgen==2023.12.18",
+    "spglib==2.3.0",
+    "pymatgen==2024.1.27",
 ]
 phonopy = [
     "phonopy==2.21.0",
-    "spglib==2.2.0",
+    "spglib==2.3.0",
 ]
-pyxtal = ["pyxtal==0.6.1"]
+pyxtal = ["pyxtal==0.6.2"]
 
 [tool.setuptools.packages.find]
 include = ["structuretoolkit*"]
diff --git a/structuretoolkit/analyse/strain.py b/structuretoolkit/analyse/strain.py
index 3632c9dda..f9ebb6aaa 100644
--- a/structuretoolkit/analyse/strain.py
+++ b/structuretoolkit/analyse/strain.py
@@ -6,7 +6,6 @@
 
 
 class Strain:
-
     """
     Calculate local strain of each atom following the Lagrangian strain tensor:
 
diff --git a/structuretoolkit/analyse/symmetry.py b/structuretoolkit/analyse/symmetry.py
index bebeeae37..2da660fee 100644
--- a/structuretoolkit/analyse/symmetry.py
+++ b/structuretoolkit/analyse/symmetry.py
@@ -24,7 +24,6 @@
 
 
 class Symmetry(dict):
-
     """
 
     Return a class for operations related to box symmetries. Main attributes:
diff --git a/structuretoolkit/visualize.py b/structuretoolkit/visualize.py
index a65ab4849..4774c4f2b 100644
--- a/structuretoolkit/visualize.py
+++ b/structuretoolkit/visualize.py
@@ -120,6 +120,7 @@ def plot3d(
     elif mode == "plotly":
         return _plot3d_plotly(
             structure=structure,
+            show_cell=show_cell,
             camera=camera,
             particle_size=particle_size,
             select_atoms=select_atoms,
@@ -143,8 +144,20 @@ def plot3d(
         raise ValueError("plot method not recognized")
 
 
+def _get_box_skeleton(cell):
+    lines_dz = np.stack(np.meshgrid(*3 * [[0, 1]], indexing="ij"), axis=-1)
+    # eight corners of a unit cube, paired as four z-axis lines
+
+    all_lines = np.reshape(
+        [np.roll(lines_dz, i, axis=-1) for i in range(3)], (-1, 2, 3)
+    )
+    # All 12 two-point lines on the unit square
+    return all_lines @ cell
+
+
 def _plot3d_plotly(
     structure,
+    show_cell=True,
     scalar_field=None,
     select_atoms=None,
     particle_size=1.0,
@@ -177,6 +190,7 @@ def _plot3d_plotly(
     """
     try:
         import plotly.express as px
+        import plotly.graph_objects as go
     except ModuleNotFoundError:
         raise ModuleNotFoundError("plotly not installed - use plot3d instead")
     if select_atoms is None:
@@ -196,6 +210,13 @@ def _plot3d_plotly(
             scale=particle_size / (0.1 * structure.get_volume() ** (1 / 3)),
         ),
     )
+    if show_cell:
+        data = fig.data
+        for lines in _get_box_skeleton(structure.cell):
+            fig = px.line_3d(**{xx: vv for xx, vv in zip(["x", "y", "z"], lines.T)})
+            fig.update_traces(line_color="#000000")
+            data = fig.data + data
+        fig = go.Figure(data=data)
     fig.layout.scene.camera.projection.type = camera
     rot = _get_orientation(view_plane).T
     rot[0, :] *= distance_from_camera * 1.25
@@ -206,6 +227,7 @@ def _plot3d_plotly(
     fig.update_layout(scene_camera=angle)
     fig.update_traces(marker=dict(line=dict(width=0.1, color="DarkSlateGrey")))
     fig.update_scenes(aspectmode="data")
+    fig.update_layout(legend={"itemsizing": "constant"})
     return fig
 
 
diff --git a/tests/test_visualize.py b/tests/test_visualize.py
index 5abc2ae58..c649865a8 100644
--- a/tests/test_visualize.py
+++ b/tests/test_visualize.py
@@ -4,7 +4,7 @@
 
 import unittest
 import numpy as np
-from structuretoolkit.visualize import _get_flattened_orientation
+from structuretoolkit.visualize import _get_flattened_orientation, _get_box_skeleton
 
 
 class TestAtoms(unittest.TestCase):
@@ -21,6 +21,25 @@ def test_get_flattened_orientation(self):
         R = np.array(_get_flattened_orientation(R, 1)).reshape(4, 4)
         self.assertAlmostEqual(np.linalg.det(R), 1)
 
+    def test_get_frame(self):
+        frame = _get_box_skeleton(np.eye(3))
+        self.assertLessEqual(
+            np.unique(frame.reshape(-1, 6), axis=0, return_counts=True)[1].max(),
+            1
+        )
+        dx, counts = np.unique(
+            np.diff(frame, axis=-2).squeeze().astype(int), axis=0, return_counts=True
+        )
+        self.assertEqual(
+            dx.ptp(), 1, msg="Frames not drawn along the nearest edges"
+        )
+        msg = (
+            "There must be four lines along each direction"
+            + " (4 x [1, 0, 0], 4 x [0, 1, 0] and 4 x [0, 0, 1])"
+        )
+        self.assertEqual(counts.min(), 4, msg=msg)
+        self.assertEqual(counts.max(), 4, msg=msg)
+
 
 if __name__ == "__main__":
     unittest.main()