Repository: GEMDAT-repos/GEMDAT (language: python)

Issue:
Consider dropping yapf for ruff formatter
Formatting is now built into ruff (https://astral.sh/blog/the-ruff-formatter) as part of [0.1.2](https://github.com/astral-sh/ruff/releases/tag/v0.1.2).

It's about 100x faster than yapf. Although it is in preview, it's considered production ready.


Discussion / hints:
Its very nice for new projects, not so much for existing ones: https://github.com/duqtools/duqtools/issues/669

The diff would be:  ` 31 files changed, 555 insertions(+), 732 deletions(-)`
I still want this 😅

Now that it is has come along a bit more, and I have used it in other repos, I also want to update GEMDAT to use ruff for formatting. I prefer the formatting style of ruff and it now has some configuration options if necessary (most importantly a toggle for single quotes)

It can replace yapf and nbqa (for linting notebooks) in the pre-commit config.

Any thoughts on this @SCiarella ?


Fix patch:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 50b9ecca..1f33abb3 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -19,13 +19,6 @@ repos:
     rev: v0.2.0
     hooks:
       - id: nbcheckorder
-  - repo: https://github.com/google/yapf
-    rev: v0.40.2
-    hooks:
-      - id: yapf
-        additional_dependencies:
-          - toml
-          - tomli
   - repo: https://github.com/myint/docformatter
     rev: v1.7.5
     hooks:
@@ -35,17 +28,11 @@ repos:
     hooks:
       - id: ruff
         args: [--fix]
+        types_or: [python, pyi, jupyter]
+      - id: ruff-format
+        types_or: [python, pyi, jupyter]
   - repo: https://github.com/pre-commit/mirrors-mypy
     rev: v1.10.0
     hooks:
       - id: mypy
         additional_dependencies: [matplotlib, MDAnalysis, numpy, pymatgen, rich, scikit-image, scipy]
-  - repo: https://github.com/nbQA-dev/nbQA
-    rev: 1.8.5
-    hooks:
-      - id: nbqa-yapf
-        additional_dependencies: [toml, yapf]
-        types_or: [jupyter]
-      - id: nbqa-ruff
-        args: [--fix, --ignore=E402]
-        types_or: [jupyter]
diff --git a/pyproject.toml b/pyproject.toml
index 225813d5..17a71019 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -99,6 +99,10 @@ known-first-party=["gemdat"]
 known-third-party = ["pymatgen"]
 required-imports = ["from __future__ import annotations"]
 
+[tool.ruff.format]
+quote-style = "single"
+indent-style = "space"
+
 [tool.mypy]
 files = ["src", "tests"]
 allow_redefinition = true
diff --git a/scripts/analyse_md.py b/scripts/analyse_md.py
index 2f046885..bfd40ea2 100644
--- a/scripts/analyse_md.py
+++ b/scripts/analyse_md.py
@@ -1,5 +1,6 @@
 """This module provides a code to deal with the original matlab code that
 Gemdat is based on."""
+
 from __future__ import annotations
 
 import matplotlib.pyplot as plt
@@ -16,7 +17,7 @@ def analyse_md(
     diff_elem: str,
     material: str,
     supercell: tuple[int, int, int] | None = None,
-    equil_time: float = 2.5E-12,
+    equil_time: float = 2.5e-12,
     diffusion_dimensions: int = 3,
     z_ion: int = 1,
     nr_parts: int = 10,
diff --git a/src/gemdat/caching.py b/src/gemdat/caching.py
index 36531dd0..56919c61 100644
--- a/src/gemdat/caching.py
+++ b/src/gemdat/caching.py
@@ -13,7 +13,6 @@ def weak_lru_cache(maxsize=128, typed=False):
     """
 
     def wrapper(func):
-
         @functools.lru_cache(maxsize, typed)
         def _func(_self, *args, **kwargs):
             return func(_self(), *args, **kwargs)
diff --git a/src/gemdat/collective.py b/src/gemdat/collective.py
index 01eb5140..0299e99f 100644
--- a/src/gemdat/collective.py
+++ b/src/gemdat/collective.py
@@ -27,12 +27,14 @@ class Collective:
         List of indices to collective jump events
     """
 
-    def __init__(self,
-                 jumps: Jumps,
-                 sites: Structure,
-                 lattice: Lattice,
-                 max_steps: int,
-                 max_dist: float = 1):
+    def __init__(
+        self,
+        jumps: Jumps,
+        sites: Structure,
+        lattice: Lattice,
+        max_steps: int,
+        max_dist: float = 1,
+    ):
         """Determine number of jumps which could show collective behaviour.
 
         Parameters
@@ -64,8 +66,7 @@ def _compute(self):
         max_dist = self.max_dist
 
         events = self.jumps.data
-        events = events.sort_values(['stop time', 'start time'],
-                                    ignore_index=True)
+        events = events.sort_values(['stop time', 'start time'], ignore_index=True)
 
         coll_jumps = []
         collective = []
@@ -73,7 +74,7 @@ def _compute(self):
 
         # Compare all pairs
         for i, event_i in events[:-1].iterrows():
-            for j, event_j in events[i + 1:].iterrows():
+            for j, event_j in events[i + 1 :].iterrows():
                 if event_j['start time'] - event_i['stop time'] > max_steps:
                     break
                 if event_i['start time'] - event_j['stop time'] > max_steps:
@@ -81,26 +82,28 @@ def _compute(self):
                 if event_i['atom index'] == event_j['atom index']:
                     continue
 
-                a = sites.frac_coords[[
-                    event_i['start site'], event_i['destination site']
-                ]]
-                b = sites.frac_coords[[
-                    event_j['start site'], event_j['destination site']
-                ]]
+                a = sites.frac_coords[
+                    [event_i['start site'], event_i['destination site']]
+                ]
+                b = sites.frac_coords[
+                    [event_j['start site'], event_j['destination site']]
+                ]
 
                 dists = lattice.get_all_distances(a, b)
 
                 if np.any(dists < max_dist):
                     collective.append((event_i, event_j))
                     coll_jumps.append(
-                        ((event_i['start site'], event_i['destination site']),
-                         (event_j['start site'], event_j['destination site'])))
+                        (
+                            (event_i['start site'], event_i['destination site']),
+                            (event_j['start site'], event_j['destination site']),
+                        )
+                    )
                     collective_matrix[i, j] = True
                     collective_matrix[j, i] = True
 
         # Get solo jumps from the collective matrix
-        self.n_solo_jumps = len(events) - np.any(collective_matrix,
-                                                 axis=0).sum()
+        self.n_solo_jumps = len(events) - np.any(collective_matrix, axis=0).sum()
 
         self.collective = collective
         self.coll_jumps = coll_jumps
@@ -119,10 +122,9 @@ def site_pair_count_matrix(self) -> np.ndarray:
         coll_jumps = self.coll_jumps
         site_pairs = self.site_pair_count_matrix_labels()
 
-        site_pair_count_matrix = np.zeros((len(site_pairs), len(site_pairs)),
-                                          dtype=int)
+        site_pair_count_matrix = np.zeros((len(site_pairs), len(site_pairs)), dtype=int)
 
-        for ((start_i, stop_i), (start_j, stop_j)) in coll_jumps:
+        for (start_i, stop_i), (start_j, stop_j) in coll_jumps:
             name_start_i = labels[start_i]
             name_stop_i = labels[stop_i]
             name_start_j = labels[start_j]
@@ -138,9 +140,7 @@ def site_pair_count_matrix(self) -> np.ndarray:
     @weak_lru_cache()
     def site_pair_count_matrix_labels(self) -> list:
         labels = self.sites.labels
-        return list({(label1, label2)
-                     for label1 in labels
-                     for label2 in labels})
+        return list({(label1, label2) for label1 in labels for label2 in labels})
 
     @weak_lru_cache()
     def multiple_collective(self) -> tuple[np.ndarray, np.ndarray]:
@@ -153,12 +153,16 @@ def multiple_collective(self) -> tuple[np.ndarray, np.ndarray]:
         multiple_collective : np.ndarray, np.ndarray
             Result of a np.unique on the collective jump pairs
         """
-        collective = np.array(dtype=[('start', int), ('stop', int)],
-                              object=[[(event_i['start site'],
-                                        event_i['destination site']),
-                                       (event_j['start site'],
-                                        event_j['destination site'])]
-                                      for event_i, event_j in self.collective])
+        collective = np.array(
+            dtype=[('start', int), ('stop', int)],
+            object=[
+                [
+                    (event_i['start site'], event_i['destination site']),
+                    (event_j['start site'], event_j['destination site']),
+                ]
+                for event_i, event_j in self.collective
+            ],
+        )
         # Sort jumps so equal jumps are orderd similarily
         collective = np.sort(collective, axis=1)
 
diff --git a/src/gemdat/io.py b/src/gemdat/io.py
index dd1fb7b8..c6b7e40e 100644
--- a/src/gemdat/io.py
+++ b/src/gemdat/io.py
@@ -1,4 +1,5 @@
 """This module contains functions to read and write data."""
+
 from __future__ import annotations
 
 import warnings
diff --git a/src/gemdat/jumps.py b/src/gemdat/jumps.py
index 4963bdaf..a485c17a 100644
--- a/src/gemdat/jumps.py
+++ b/src/gemdat/jumps.py
@@ -19,9 +19,9 @@
     from mypy_extensions import DefaultNamedArg
 
 
-def _generic_transitions_to_jumps(transitions: Transitions,
-                                  *,
-                                  minimal_residence: int = 0) -> pd.DataFrame:
+def _generic_transitions_to_jumps(
+    transitions: Transitions, *, minimal_residence: int = 0
+) -> pd.DataFrame:
     """Generic function to convert transition events to jumps.
 
     Parameters
@@ -52,13 +52,11 @@ def _generic_transitions_to_jumps(transitions: Transitions,
             if candidate_jump['atom index'] != event['atom index']:
                 jumps.append(candidate_jump)
                 candidate_jump = None
-            elif event['start time'] - candidate_jump[
-                    'stop time'] >= minimal_residence:
+            elif event['start time'] - candidate_jump['stop time'] >= minimal_residence:
                 jumps.append(candidate_jump)
                 candidate_jump = None
                 fromevent = None
-            elif candidate_jump['destination site'] != event[
-                    'destination site']:
+            elif candidate_jump['destination site'] != event['destination site']:
                 candidate_jump = None
 
         # Specify the start of a jump if we encounter one
@@ -96,8 +94,7 @@ def _generic_transitions_to_jumps(transitions: Transitions,
     jumps = pd.DataFrame(data=jumps)
 
     # remove jumps where start == destination
-    jumps = jumps[jumps['start site'] !=
-                  jumps['destination site']].reset_index()
+    jumps = jumps[jumps['start site'] != jumps['destination site']].reset_index()
 
     # remove old index
     del jumps['index']
@@ -107,15 +104,15 @@ def _generic_transitions_to_jumps(transitions: Transitions,
 
 
 class Jumps:
-
-    def __init__(self,
-                 transitions: Transitions,
-                 *,
-                 conversion_method: Callable[
-                     [Transitions,
-                      DefaultNamedArg(int, 'minimal_residence')],
-                     pd.DataFrame] = _generic_transitions_to_jumps,
-                 minimal_residence: int = 0):
+    def __init__(
+        self,
+        transitions: Transitions,
+        *,
+        conversion_method: Callable[
+            [Transitions, DefaultNamedArg(int, 'minimal_residence')], pd.DataFrame
+        ] = _generic_transitions_to_jumps,
+        minimal_residence: int = 0,
+    ):
         """Analyze transitions and classify them as jumps.
 
         Parameters
@@ -134,8 +131,7 @@ def __init__(self,
         self.trajectory = transitions.diff_trajectory
         self.sites = transitions.sites
         self.conversion_method = conversion_method
-        self.data = conversion_method(transitions,
-                                      minimal_residence=minimal_residence)
+        self.data = conversion_method(transitions, minimal_residence=minimal_residence)
 
     @property
     def n_jumps(self) -> int:
@@ -190,8 +186,7 @@ def jump_diffusivity(self, dimensions: int) -> float:
         pdist = lattice.get_all_distances(sites.frac_coords, sites.frac_coords)
 
         jump_diff = np.sum(pdist**2 * self.matrix())
-        jump_diff *= angstrom**2 / (2 * dimensions * self.n_floating *
-                                    total_time)
+        jump_diff *= angstrom**2 / (2 * dimensions * self.n_floating * total_time)
 
         jump_diff = FloatWithUnit(jump_diff, 'm^2 s^-1')
 
@@ -206,8 +201,9 @@ def matrix(self) -> np.ndarray:
         transitions_matrix : np.ndarray
             Square matrix with number of each transitions
         """
-        return _calculate_transitions_matrix(self.data,
-                                             n_sites=self.transitions.n_sites)
+        return _calculate_transitions_matrix(
+            self.data, n_sites=self.transitions.n_sites
+        )
 
     @weak_lru_cache()
     def collective(self, max_dist: float = 1) -> Collective:
@@ -264,21 +260,19 @@ def activation_energies(self, n_parts: int = 10) -> pd.DataFrame:
         atom_locations_parts = [
             part.atom_locations() for part in self.transitions.split(n_parts)
         ]
-        jumps_counter_parts = [
-            part.jumps_counter() for part in self.split(n_parts)
-        ]
+        jumps_counter_parts = [part.jumps_counter() for part in self.split(n_parts)]
         n_floating = self.n_floating
 
         for site_pair in self.site_pairs:
             site_start, site_stop = site_pair
 
-            n_jumps = np.array(
-                [part[site_pair] for part in jumps_counter_parts])
+            n_jumps = np.array([part[site_pair] for part in jumps_counter_parts])
 
             part_time = trajectory.total_time / n_parts
 
             atom_percentage = np.array(
-                [part[site_start] for part in atom_locations_parts])
+                [part[site_start] for part in atom_locations_parts]
+            )
 
             denom = atom_percentage * n_floating * part_time
 
@@ -288,11 +282,13 @@ def activation_energies(self, n_parts: int = 10) -> pd.DataFrame:
             if site_start == site_stop:
                 eff_rate /= 2
 
-            e_act_arr = -np.log(eff_rate / attempt_freq) * (
-                Boltzmann * temperature) / elementary_charge
+            e_act_arr = (
+                -np.log(eff_rate / attempt_freq)
+                * (Boltzmann * temperature)
+                / elementary_charge
+            )
 
-            dct[site_start, site_stop] = np.mean(e_act_arr), np.std(e_act_arr,
-                                                                    ddof=1)
+            dct[site_start, site_stop] = np.mean(e_act_arr), np.std(e_act_arr, ddof=1)
 
         df = pd.DataFrame(dct).T
         df.columns = ('energy', 'std')
@@ -308,8 +304,14 @@ def jumps_counter(self) -> Counter:
             Dictionary with number of jumps per sites combination
         """
         labels = self.sites.labels
-        jumps = Counter([(labels[i], labels[j]) for _, (
-            i, j) in self.data[['start site', 'destination site']].iterrows()])
+        jumps = Counter(
+            [
+                (labels[i], labels[j])
+                for _, (i, j) in self.data[
+                    ['start site', 'destination site']
+                ].iterrows()
+            ]
+        )
         return jumps
 
     def split(self, n_parts: int) -> list[Jumps]:
@@ -326,10 +328,7 @@ def split(self, n_parts: int) -> list[Jumps]:
         """
         parts = self.transitions.split(n_parts)
 
-        return [
-            Jumps(part, conversion_method=self.conversion_method)
-            for part in parts
-        ]
+        return [Jumps(part, conversion_method=self.conversion_method) for part in parts]
 
     @weak_lru_cache()
     def rates(self, n_parts: int = 10) -> pd.DataFrame:
@@ -363,24 +362,29 @@ def rates(self, n_parts: int = 10) -> pd.DataFrame:
     def plot_jumps_vs_distance(self, **kwargs):
         """See [gemdat.plots.jumps_vs_distance][] for more information."""
         from gemdat import plots
+
         return plots.jumps_vs_distance(jumps=self, **kwargs)
 
     def plot_jumps_vs_time(self, **kwargs):
         """See [gemdat.plots.jumps_vs_time][] for more information."""
         from gemdat import plots
+
         return plots.jumps_vs_time(jumps=self, **kwargs)
 
     def plot_collective_jumps(self, **kwargs):
         """See [gemdat.plots.collective_jumps][] for more information."""
         from gemdat import plots
+
         return plots.collective_jumps(jumps=self, **kwargs)
 
     def plot_jumps_3d(self, **kwargs):
         """See [gemdat.plots.jumps_3d][] for more information."""
         from gemdat import plots
+
         return plots.jumps_3d(jumps=self, **kwargs)
 
     def plot_jumps_3d_animation(self, **kwargs):
         """See [gemdat.plots.jumps_3d_animation][] for more information."""
         from gemdat import plots
+
         return plots.jumps_3d_animation(jumps=self, **kwargs)
diff --git a/src/gemdat/orientations.py b/src/gemdat/orientations.py
index 5d1dbfbd..aa1e45c4 100644
--- a/src/gemdat/orientations.py
+++ b/src/gemdat/orientations.py
@@ -26,6 +26,7 @@ class Orientations:
     vectors: np.ndarray
         Vectors representing orientation direction
     """
+
     trajectory: Trajectory
     center_type: str
     satellite_type: str
@@ -71,14 +72,20 @@ def _distances(self) -> np.ndarray:
         central_start_coord = self._trajectory_cent.base_positions
         satellite_start_coord = self._trajectory_sat.base_positions
         lattice = self.trajectory.get_lattice()
-        distance = np.array([[
-            lattice.get_all_distances(central, satellite)
-            for satellite in satellite_start_coord
-        ] for central in central_start_coord])
+        distance = np.array(
+            [
+                [
+                    lattice.get_all_distances(central, satellite)
+                    for satellite in satellite_start_coord
+                ]
+                for central in central_start_coord
+            ]
+        )
         return distance
 
-    def _matching_matrix(self, distance: np.ndarray,
-                         frac_coord_cent: np.ndarray) -> np.ndarray:
+    def _matching_matrix(
+        self, distance: np.ndarray, frac_coord_cent: np.ndarray
+    ) -> np.ndarray:
         """Determine which satellite atoms are close enough to central atom to
         be connected.
 
@@ -97,14 +104,14 @@ def _matching_matrix(self, distance: np.ndarray,
         match_criteria = 1.5 * np.min(distance)
 
         distance_match = np.where(distance < match_criteria, distance, 0)
-        matching_matrix = np.zeros((len(frac_coord_cent[0, :, 0]), 4),
-                                   dtype=int)
+        matching_matrix = np.zeros((len(frac_coord_cent[0, :, 0]), 4), dtype=int)
         for k in range(len(frac_coord_cent[0, :, 0])):
             matching_matrix[k, :] = np.where(distance_match[k, :] != 0)[0][:4]
         return matching_matrix
 
-    def _central_satellite_matrix(self, distance: np.ndarray,
-                                  frac_coord_cent: np.ndarray) -> np.ndarray:
+    def _central_satellite_matrix(
+        self, distance: np.ndarray, frac_coord_cent: np.ndarray
+    ) -> np.ndarray:
         """Get the combinations of central atoms and satellite atoms in a
         matrix.
 
@@ -126,8 +133,9 @@ def _central_satellite_matrix(self, distance: np.ndarray,
 
         # index_central_atoms = np.arange(self.nr_central_atoms)
         matching_matrix = self._matching_matrix(distance, frac_coord_cent)
-        combinations = np.array([(i, j) for i in index_central_atoms
-                                 for j in matching_matrix[i, :]])
+        combinations = np.array(
+            [(i, j) for i in index_central_atoms for j in matching_matrix[i, :]]
+        )
         return combinations
 
     def _fractional_directions(self, distance: np.ndarray) -> np.ndarray:
@@ -147,8 +155,7 @@ def _fractional_directions(self, distance: np.ndarray) -> np.ndarray:
         frac_coord_cent = self._trajectory_cent.positions
         frac_coord_sat = self._trajectory_sat.positions
 
-        combinations = self._central_satellite_matrix(distance,
-                                                      frac_coord_cent)
+        combinations = self._central_satellite_matrix(distance, frac_coord_cent)
 
         sat = frac_coord_sat[:, combinations[:, 1], :]
         cent = frac_coord_cent[:, combinations[:, 0], :]
@@ -168,14 +175,13 @@ def normalize(self) -> Orientations:
         -------
         Orientations
         """
-        vectors = self.vectors / np.linalg.norm(
-            self.vectors, axis=-1, keepdims=True)
+        vectors = self.vectors / np.linalg.norm(self.vectors, axis=-1, keepdims=True)
 
         return replace(self, in_vectors=vectors)
 
-    def symmetrize(self,
-                   sym_group: str | None = None,
-                   sym_ops: np.ndarray | None = None) -> Orientations:
+    def symmetrize(
+        self, sym_group: str | None = None, sym_ops: np.ndarray | None = None
+    ) -> Orientations:
         """Apply symmetry elements to the trajectory to improve statistics.
 
         One of `sym_group` and `sym_ops` must be supplied.
@@ -198,12 +204,13 @@ def symmetrize(self,
             sym_ops = sym_ops
         elif sym_group:
             g = PointGroup(sym_group)
-            sym_ops = np.array([
-                element.rotation_matrix for element in g.symmetry_ops
-            ]).transpose(1, 2, 0)
+            sym_ops = np.array(
+                [element.rotation_matrix for element in g.symmetry_ops]
+            ).transpose(1, 2, 0)
         else:
             raise ValueError(
-                'At least one of `sym_group` or `sym_ops` must be provided')
+                'At least one of `sym_group` or `sym_ops` must be provided'
+            )
 
         n_ts = self.vectors.shape[0]
         n_bonds = self.vectors.shape[1]
@@ -258,21 +265,23 @@ def autocorrelation(self):
     def plot_rectilinear(self, **kwargs):
         """See [gemdat.plots.rectilinear][] for more info."""
         from gemdat import plots
+
         return plots.rectilinear(orientations=self, **kwargs)
 
     def plot_bond_length_distribution(self, **kwargs):
         """See [gemdat.plots.bond_length_distribution][] for more info."""
         from gemdat import plots
+
         return plots.bond_length_distribution(orientations=self, **kwargs)
 
     def plot_autocorrelation(self, **kwargs):
         """See [gemdat.plots.unit_vector_autocorrelation][] for more info."""
         from gemdat import plots
+
         return plots.autocorrelation(orientations=self, **kwargs)
 
 
-def calculate_spherical_areas(shape: tuple[int, int],
-                              radius: float = 1) -> np.ndarray:
+def calculate_spherical_areas(shape: tuple[int, int], radius: float = 1) -> np.ndarray:
     """Calculate the areas of a section of a sphere, defined in spherical
     coordinates. Useful for normalization purposes.
 
@@ -296,9 +305,12 @@ def calculate_spherical_areas(shape: tuple[int, int],
 
     for i in range(shape[1]):
         for j in range(shape[0]):
-
-            areas[j, i] = (radius**2) * azimuthal_increment * np.sin(
-                np.deg2rad(elevation_angles[j])) * elevation_increment
-            #hacky way to get rid of singularity on poles
+            areas[j, i] = (
+                (radius**2)
+                * azimuthal_increment
+                * np.sin(np.deg2rad(elevation_angles[j]))
+                * elevation_increment
+            )
+            # hacky way to get rid of singularity on poles
             areas[0, :] = areas[-1, 0]
     return areas
diff --git a/src/gemdat/path.py b/src/gemdat/path.py
index 6e2665e5..69bd99ea 100644
--- a/src/gemdat/path.py
+++ b/src/gemdat/path.py
@@ -83,8 +83,7 @@ def wrapped_sites(self) -> list[tuple[int, int, int]]:
             Voxel coordinates wrapped to bounding box.
         """
         if not self.dims:
-            raise AttributeError(
-                f'Dimensions are needed for this method {self.dims=}')
+            raise AttributeError(f'Dimensions are needed for this method {self.dims=}')
         xdim, ydim, zdim = self.dims
         return [(x % xdim, y % xdim, z % xdim) for x, y, z in self.sites]
 
@@ -99,8 +98,7 @@ def frac_sites(self) -> np.ndarray:
             Fractional coordinates for sites
         """
         if not self.dims:
-            raise AttributeError(
-                f'Dimensions are needed for this method {self.dims=}')
+            raise AttributeError(f'Dimensions are needed for this method {self.dims=}')
         sites = self.wrapped_sites()
         return (np.array(sites) + 0.5) / np.array(self.dims)
 
@@ -123,7 +121,8 @@ def path_over_structure(
         frac_sites = np.array(self.frac_sites())
 
         nearest_structure_tree, nearest_structure_map = nearest_structure_reference(
-            structure)
+            structure
+        )
 
         # Get the indices of the nearest structure sites to the path sites
         nearest_structure_indices = [
@@ -149,17 +148,21 @@ def stop_site(self) -> tuple[int, int, int]:
     def plot_energy_along_path(self, **kwargs):
         """See [gemdat.plots.energy_along_path][] for more info."""
         from gemdat import plots
+
         return plots.energy_along_path(path=self, **kwargs)
 
     def plot_path_on_grid(self, **kwargs):
         """See [gemdat.plots.path_on_grid][] for more info."""
         from gemdat import plots
+
         return plots.path_on_grid(path=self, **kwargs)
 
 
-def free_energy_graph(F: np.ndarray | FreeEnergyVolume,
-                      max_energy_threshold: float = 1e20,
-                      diagonal: bool = True) -> nx.Graph:
+def free_energy_graph(
+    F: np.ndarray | FreeEnergyVolume,
+    max_energy_threshold: float = 1e20,
+    diagonal: bool = True,
+) -> nx.Graph:
     """Compute the graph of the free energy for networkx functions.
 
     Parameters
@@ -178,15 +181,30 @@ def free_energy_graph(F: np.ndarray | FreeEnergyVolume,
     """
 
     # Define possible movements in 3D space
-    movements = np.array([(1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0),
-                          (0, 0, 1), (0, 0, -1)])
+    movements = np.array(
+        [(1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1)]
+    )
     if diagonal:
-        diagonal_movements = np.array([(1, 1, 0), (-1, -1, 0), (1, -1, 0),
-                                       (-1, 1, 0), (1, 0, 1), (-1, 0, -1),
-                                       (1, 0, -1), (-1, 0, 1), (0, 1, 1),
-                                       (0, -1, -1), (0, 1, -1), (0, -1, 1),
-                                       (1, 1, 1), (-1, -1, -1), (1, -1, -1),
-                                       (-1, 1, 1)])
+        diagonal_movements = np.array(
+            [
+                (1, 1, 0),
+                (-1, -1, 0),
+                (1, -1, 0),
+                (-1, 1, 0),
+                (1, 0, 1),
+                (-1, 0, -1),
+                (1, 0, -1),
+                (-1, 0, 1),
+                (0, 1, 1),
+                (0, -1, -1),
+                (0, 1, -1),
+                (0, -1, 1),
+                (1, 1, 1),
+                (-1, -1, -1),
+                (1, -1, -1),
+                (-1, 1, 1),
+            ]
+        )
         movements = np.vstack((movements, diagonal_movements))
 
     G = nx.Graph()
@@ -208,16 +226,14 @@ def free_energy_graph(F: np.ndarray | FreeEnergyVolume,
                 else:
                     weight_exp = max_energy_threshold
 
-                G.add_edge(node,
-                           neighbor,
-                           weight=weight,
-                           weight_exp=weight_exp)
+                G.add_edge(node, neighbor, weight=weight, weight_exp=weight_exp)
 
     return G
 
 
-_PATHFINDING_METHODS = Literal['dijkstra', 'bellman-ford', 'minmax-energy',
-                               'dijkstra-exp', 'simple']
+_PATHFINDING_METHODS = Literal[
+    'dijkstra', 'bellman-ford', 'minmax-energy', 'dijkstra-exp', 'simple'
+]
 
 
 def calculate_path_difference(path1: list, path2: list) -> float:
@@ -249,8 +265,7 @@ def calculate_path_difference(path1: list, path2: list) -> float:
     return 1 - (shared_nodes / len(shortest))
 
 
-def _paths_too_similar(path: list, list_of_paths: list,
-                       min_diff: float) -> bool:
+def _paths_too_similar(path: list, list_of_paths: list, min_diff: float) -> bool:
     """Check if the path is too similar to the other paths.
 
     Parameters
@@ -321,10 +336,9 @@ def optimal_n_paths(
     list_of_paths = [best_path]
 
     # Compute the iterator over all the short paths
-    all_paths = nx.shortest_simple_paths(F_graph,
-                                         source=start,
-                                         target=stop,
-                                         weight='weight')
+    all_paths = nx.shortest_simple_paths(
+        F_graph, source=start, target=stop, weight='weight'
+    )
 
     # Attempt to find the Np shortest paths
     for idx, path in enumerate(all_paths):
@@ -383,17 +397,14 @@ def optimal_path(
     start = tuple(start)
     stop = tuple(stop)
 
-    optimal_path = nx.shortest_path(F_graph,
-                                    source=start,
-                                    target=stop,
-                                    weight=weight,
-                                    method=method)
+    optimal_path = nx.shortest_path(
+        F_graph, source=start, target=stop, weight=weight, method=method
+    )
 
     if method == 'minmax-energy':
-        optimal_path = _optimal_path_minmax_energy(F_graph,
-                                                   start=start,
-                                                   stop=stop,
-                                                   optimal_path=optimal_path)
+        optimal_path = _optimal_path_minmax_energy(
+            F_graph, start=start, stop=stop, optimal_path=optimal_path
+        )
     elif method not in ('dijkstra', 'bellman-ford', 'dijkstra-exp'):
         raise ValueError(f'Unknown method {method}')
 
@@ -444,8 +455,7 @@ def _optimal_path_minmax_energy(
             target=stop,
             weight='weight',
         )
-        minmax_energy = max(
-            [F_graph.nodes[node]['energy'] for node in pruned_path])
+        minmax_energy = max([F_graph.nodes[node]['energy'] for node in pruned_path])
 
         if minmax_energy < max_energy:
             optimal_path = pruned_path
diff --git a/src/gemdat/plots/__init__.py b/src/gemdat/plots/__init__.py
index 3dbb71a8..db5e5a9a 100644
--- a/src/gemdat/plots/__init__.py
+++ b/src/gemdat/plots/__init__.py
@@ -1,4 +1,5 @@
 """This module contains all the plots that Gemdat can generate."""
+
 from __future__ import annotations
 
 from .matplotlib import (
diff --git a/src/gemdat/plots/_shared.py b/src/gemdat/plots/_shared.py
index 7386ae02..b02d8944 100644
--- a/src/gemdat/plots/_shared.py
+++ b/src/gemdat/plots/_shared.py
@@ -13,7 +13,8 @@
 
 
 def _mean_displacements_per_element(
-        trajectory: Trajectory) -> dict[str, tuple[np.ndarray, np.ndarray]]:
+    trajectory: Trajectory,
+) -> dict[str, tuple[np.ndarray, np.ndarray]]:
     """Calculate mean displacements per element type.
 
     Helper for `displacement_per_atom`.
@@ -21,8 +22,7 @@ def _mean_displacements_per_element(
     species = trajectory.species
 
     grouped = defaultdict(list)
-    for sp, distances in zip(species,
-                             trajectory.distances_from_base_position()):
+    for sp, distances in zip(species, trajectory.distances_from_base_position()):
         grouped[sp.symbol].append(distances)
 
     means = {}
@@ -35,10 +35,9 @@ def _mean_displacements_per_element(
     return means
 
 
-def _orientations_to_histogram(orientations: Orientations,
-                               /,
-                               *,
-                               bins: int = 50) -> pd.DataFrame:
+def _orientations_to_histogram(
+    orientations: Orientations, /, *, bins: int = 50
+) -> pd.DataFrame:
     """Calculate bond lenth histogram from `Orientations`
 
     Helper for `bond_length_distribution`.
@@ -49,17 +48,19 @@ def _orientations_to_histogram(orientations: Orientations,
     hist, edges = np.histogram(bond_lengths, bins=bins, density=True)
     bin_centers = bin_centers = (edges[:-1] + edges[1:]) / 2
 
-    return pd.DataFrame({
-        'prob': hist,
-        'center': bin_centers,
-        'left_edge': edges[:-1],
-        'right_edge': edges[1:]
-    })
+    return pd.DataFrame(
+        {
+            'prob': hist,
+            'center': bin_centers,
+            'left_edge': edges[:-1],
+            'right_edge': edges[1:],
+        }
+    )
 
 
-def _fit_skewnorm_to_hist(df: pd.DataFrame,
-                          /,
-                          steps: int = 100) -> tuple[np.ndarray, np.ndarray]:
+def _fit_skewnorm_to_hist(
+    df: pd.DataFrame, /, steps: int = 100
+) -> tuple[np.ndarray, np.ndarray]:
     """Fit skewnorm to bond lenth histogram from `Orientations`
 
     Helper for `bond_length_distribution`.
@@ -67,10 +68,7 @@ def _fit_skewnorm_to_hist(df: pd.DataFrame,
     bin_centers = df['center']
     hist = df['prob']
 
-    params, covariance = curve_fit(skewnorm.pdf,
-                                   bin_centers,
-                                   hist,
-                                   p0=[1.5, 1, 1.5])
+    params, covariance = curve_fit(skewnorm.pdf, bin_centers, hist, p0=[1.5, 1, 1.5])
 
     x = np.linspace(min(bin_centers), max(bin_centers), steps)
     y = skewnorm.pdf(x, *params)
diff --git a/src/gemdat/plots/matplotlib/__init__.py b/src/gemdat/plots/matplotlib/__init__.py
index 98e9f516..336c3cce 100644
--- a/src/gemdat/plots/matplotlib/__init__.py
+++ b/src/gemdat/plots/matplotlib/__init__.py
@@ -1,5 +1,6 @@
 """This module contains all the plots that Gemdat can generate with
 matplotlib."""
+
 from __future__ import annotations
 
 from ._displacement_histogram import displacement_histogram
diff --git a/src/gemdat/plots/matplotlib/_autocorrelation.py b/src/gemdat/plots/matplotlib/_autocorrelation.py
index e2aa9f89..e2687bc8 100644
--- a/src/gemdat/plots/matplotlib/_autocorrelation.py
+++ b/src/gemdat/plots/matplotlib/_autocorrelation.py
@@ -3,8 +3,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-from gemdat.orientations import (
-    Orientations, )
+from gemdat.orientations import Orientations
 
 
 def autocorrelation(
@@ -49,11 +48,13 @@ def autocorrelation(
             ax.plot(tgrid, ac_i, lw=0.1, c=last_color, label=label)
 
     if show_shaded:
-        ax.fill_between(tgrid,
-                        ac_mean - ac_std,
-                        ac_mean + ac_std,
-                        alpha=0.2,
-                        label='Standard deviation')
+        ax.fill_between(
+            tgrid,
+            ac_mean - ac_std,
+            ac_mean + ac_std,
+            alpha=0.2,
+            label='Standard deviation',
+        )
 
     ax.set_xlabel('Time lag (ps)')
     ax.set_ylabel('Autocorrelation')
diff --git a/src/gemdat/plots/matplotlib/_bond_length_distribution.py b/src/gemdat/plots/matplotlib/_bond_length_distribution.py
index 59610c43..7a4a648c 100644
--- a/src/gemdat/plots/matplotlib/_bond_length_distribution.py
+++ b/src/gemdat/plots/matplotlib/_bond_length_distribution.py
@@ -6,9 +6,9 @@
 from .._shared import _orientations_to_histogram, _fit_skewnorm_to_hist
 
 
-def bond_length_distribution(*,
-                             orientations: Orientations,
-                             bins: int = 50) -> plt.Figure:
+def bond_length_distribution(
+    *, orientations: Orientations, bins: int = 50
+) -> plt.Figure:
     """Plot the bond length probability distribution.
 
     Parameters
@@ -35,12 +35,14 @@ def bond_length_distribution(*,
 
     ax.plot(x, y, 'r-', lw=2, label='Skewed Gaussian Fit')
 
-    ax.bar(x=centers,
-           height=prob,
-           width=bin_width * 0.8,
-           edgecolor='none',
-           alpha=0.7,
-           label='Probability')
+    ax.bar(
+        x=centers,
+        height=prob,
+        width=bin_width * 0.8,
+        edgecolor='none',
+        alpha=0.7,
+        label='Probability',
+    )
 
     ax.set_xlabel('Bond length (Å)')
     ax.set_ylabel(r'Probability density (Å$^{-1}$)')
diff --git a/src/gemdat/plots/matplotlib/_displacement_histogram.py b/src/gemdat/plots/matplotlib/_displacement_histogram.py
index af214536..44fffc39 100644
--- a/src/gemdat/plots/matplotlib/_displacement_histogram.py
+++ b/src/gemdat/plots/matplotlib/_displacement_histogram.py
@@ -20,8 +20,10 @@ def displacement_histogram(trajectory: Trajectory) -> plt.Figure:
     """
     fig, ax = plt.subplots()
     ax.hist(trajectory.distances_from_base_position()[:, -1])
-    ax.set(title='Displacement per element',
-           xlabel='Displacement (Å)',
-           ylabel='Nr. of atoms')
+    ax.set(
+        title='Displacement per element',
+        xlabel='Displacement (Å)',
+        ylabel='Nr. of atoms',
+    )
 
     return fig
diff --git a/src/gemdat/plots/matplotlib/_displacement_per_atom.py b/src/gemdat/plots/matplotlib/_displacement_per_atom.py
index 2b33ae76..89e43a4a 100644
--- a/src/gemdat/plots/matplotlib/_displacement_per_atom.py
+++ b/src/gemdat/plots/matplotlib/_displacement_per_atom.py
@@ -23,8 +23,6 @@ def displacement_per_atom(*, trajectory: Trajectory) -> plt.Figure:
     for distances in trajectory.distances_from_base_position():
         ax.plot(distances, lw=0.3)
 
-    ax.set(title='Displacement per site',
-           xlabel='Time step',
-           ylabel='Displacement (Å)')
+    ax.set(title='Displacement per site', xlabel='Time step', ylabel='Displacement (Å)')
 
     return fig
diff --git a/src/gemdat/plots/matplotlib/_displacement_per_element.py b/src/gemdat/plots/matplotlib/_displacement_per_element.py
index 626c7248..8d35dfd3 100644
--- a/src/gemdat/plots/matplotlib/_displacement_per_element.py
+++ b/src/gemdat/plots/matplotlib/_displacement_per_element.py
@@ -28,8 +28,8 @@ def displacement_per_element(*, trajectory: Trajectory) -> plt.Figure:
         ax.plot(mean, lw=0.3, label=symbol)
 
     ax.legend()
-    ax.set(title='Displacement per element',
-           xlabel='Time step',
-           ylabel='Displacement (Å)')
+    ax.set(
+        title='Displacement per element', xlabel='Time step', ylabel='Displacement (Å)'
+    )
 
     return fig
diff --git a/src/gemdat/plots/matplotlib/_energy_along_path.py b/src/gemdat/plots/matplotlib/_energy_along_path.py
index fb258303..f6844354 100644
--- a/src/gemdat/plots/matplotlib/_energy_along_path.py
+++ b/src/gemdat/plots/matplotlib/_energy_along_path.py
@@ -43,8 +43,7 @@ def energy_along_path(
         prev = nearest_sites[0]
         for i, site in enumerate(nearest_sites):
             if (site.coords != prev.coords).any() or i == 0:
-                sitecoord_xlabel.append(', '.join(f'{val:.1f}'
-                                                  for val in site.coords))
+                sitecoord_xlabel.append(', '.join(f'{val:.1f}' for val in site.coords))
                 site_xlabel.append(site.label)
             else:
                 sitecoord_xlabel.append('')
@@ -52,18 +51,17 @@ def energy_along_path(
 
             prev = site
 
-        non_empty_ticks = [
-            i for i, label in enumerate(sitecoord_xlabel) if label != ''
-        ]
+        non_empty_ticks = [i for i, label in enumerate(sitecoord_xlabel) if label != '']
 
         extra_ticks = non_empty_ticks.copy()
         extra_ticks.append(ax.get_xlim()[1])
-        centered_ticks = [(extra_ticks[i] + extra_ticks[i + 1]) / 2
-                          for i in range(len(extra_ticks) - 1)]
+        centered_ticks = [
+            (extra_ticks[i] + extra_ticks[i + 1]) / 2
+            for i in range(len(extra_ticks) - 1)
+        ]
 
         ax.set_xticks(centered_ticks)
-        ax.set_xticklabels([sitecoord_xlabel[i] for i in non_empty_ticks],
-                           rotation=45)
+        ax.set_xticklabels([sitecoord_xlabel[i] for i in non_empty_ticks], rotation=45)
 
         for start, stop in zip(non_empty_ticks[::2], non_empty_ticks[1::2]):
             ax.axvspan(start, stop, facecolor='lightgray', edgecolor='none')
@@ -71,17 +69,14 @@ def energy_along_path(
         ax_up = ax.twiny()
         ax_up.set_xlim(ax.get_xlim())
         ax_up.set_xticks(centered_ticks)
-        ax_up.set_xticklabels([site_xlabel[i] for i in non_empty_ticks],
-                              rotation=45)
+        ax_up.set_xticklabels([site_xlabel[i] for i in non_empty_ticks], rotation=45)
         ax_up.get_yaxis().set_visible(False)
 
     if other_paths:
         for idx, path in enumerate(other_paths):
             if path.energy is None:
                 raise ValueError('Pathway does not contain energy data')
-            ax.plot(range(len(path.energy)),
-                    path.energy,
-                    label=f'Alternative {idx+1}')
+            ax.plot(range(len(path.energy)), path.energy, label=f'Alternative {idx+1}')
 
         ax.legend(fontsize=8)
 
diff --git a/src/gemdat/plots/matplotlib/_frequency_vs_occurence.py b/src/gemdat/plots/matplotlib/_frequency_vs_occurence.py
index c89531f2..6b7d53a3 100644
--- a/src/gemdat/plots/matplotlib/_frequency_vs_occurence.py
+++ b/src/gemdat/plots/matplotlib/_frequency_vs_occurence.py
@@ -51,11 +51,14 @@ def frequency_vs_occurence(*, trajectory: Trajectory) -> plt.Figure:
             0,
             y_max,
             colors='red',
-            linestyles='dashed')
-
-    ax.set(title='Frequency vs Occurence',
-           xlabel='Frequency (Hz)',
-           ylabel='Occurrence (a.u.)')
+            linestyles='dashed',
+        )
+
+    ax.set(
+        title='Frequency vs Occurence',
+        xlabel='Frequency (Hz)',
+        ylabel='Occurrence (a.u.)',
+    )
 
     ax.set_ylim([0, y_max])
     ax.set_xlim([-0.1e13, 2.5e13])
diff --git a/src/gemdat/plots/matplotlib/_jumps_3d.py b/src/gemdat/plots/matplotlib/_jumps_3d.py
index 4a540519..43f9f911 100644
--- a/src/gemdat/plots/matplotlib/_jumps_3d.py
+++ b/src/gemdat/plots/matplotlib/_jumps_3d.py
@@ -27,7 +27,6 @@ def jumps_3d(*, jumps: Jumps) -> plt.Figure:
     sites = jumps.sites
 
     class LabelItems:
-
         def __init__(self, labels, coords):
             self.labels = labels
             self.coords = coords
@@ -43,15 +42,18 @@ def items(self):
 
     site_labels = LabelItems(jumps.sites.labels, coords)
 
-    xyz_labels = LabelItems('OABC', [[-0.1, -0.1, -0.1], [1.1, -0.1, -0.1],
-                                     [-0.1, 1.1, -0.1], [-0.1, -0.1, 1.1]])
+    xyz_labels = LabelItems(
+        'OABC',
+        [
+            [-0.1, -0.1, -0.1],
+            [1.1, -0.1, -0.1],
+            [-0.1, 1.1, -0.1],
+            [-0.1, -0.1, 1.1],
+        ],
+    )
 
     plotter.plot_lattice_vectors(lattice, ax=ax, linewidth=1)
-    plotter.plot_labels(xyz_labels,
-                        lattice=lattice,
-                        ax=ax,
-                        color='green',
-                        size=12)
+    plotter.plot_labels(xyz_labels, lattice=lattice, ax=ax, color='green', size=12)
     plotter.plot_points(coords, lattice=lattice, ax=ax)
 
     for i, j in zip(*np.triu_indices(len(coords), k=1)):
@@ -73,17 +75,9 @@ def items(self):
             lines = [(coord_i, coord_j)]
 
         for line in lines:
-            plotter.plot_path(line,
-                              lattice=lattice,
-                              ax=ax,
-                              color='red',
-                              linewidth=lw)
-
-    plotter.plot_labels(site_labels,
-                        lattice=lattice,
-                        ax=ax,
-                        color='black',
-                        size=8)
+            plotter.plot_path(line, lattice=lattice, ax=ax, color='red', linewidth=lw)
+
+    plotter.plot_labels(site_labels, lattice=lattice, ax=ax, color='black', size=8)
 
     ax.set(
         title='Jumps between sites',
diff --git a/src/gemdat/plots/matplotlib/_jumps_3d_animation.py b/src/gemdat/plots/matplotlib/_jumps_3d_animation.py
index 74031e95..0c474229 100644
--- a/src/gemdat/plots/matplotlib/_jumps_3d_animation.py
+++ b/src/gemdat/plots/matplotlib/_jumps_3d_animation.py
@@ -48,7 +48,6 @@ def jumps_3d_animation(
     trajectory = jumps.trajectory
 
     class LabelItems:
-
         def __init__(self, labels, coords):
             self.labels = labels
             self.coords = coords
@@ -65,24 +64,24 @@ def items(self):
     fig = plt.figure()
     ax = fig.add_subplot(projection='3d')
 
-    xyz_labels = LabelItems('OABC', [[-0.1, -0.1, -0.1], [1.1, -0.1, -0.1],
-                                     [-0.1, 1.1, -0.1], [-0.1, -0.1, 1.1]])
+    xyz_labels = LabelItems(
+        'OABC',
+        [
+            [-0.1, -0.1, -0.1],
+            [1.1, -0.1, -0.1],
+            [-0.1, 1.1, -0.1],
+            [-0.1, -0.1, 1.1],
+        ],
+    )
 
     plotter.plot_lattice_vectors(lattice, ax=ax, linewidth=1)
 
-    plotter.plot_labels(xyz_labels,
-                        lattice=lattice,
-                        ax=ax,
-                        color='green',
-                        size=12)
+    plotter.plot_labels(xyz_labels, lattice=lattice, ax=ax, color='green', size=12)
 
     assert len(ax.collections) == 0
-    plotter.plot_points(coords,
-                        lattice=lattice,
-                        ax=ax,
-                        s=50,
-                        color='white',
-                        edgecolor='black')
+    plotter.plot_points(
+        coords, lattice=lattice, ax=ax, s=50, color='white', edgecolor='black'
+    )
     points = ax.collections
 
     events = jumps.data.sort_values('start time', ignore_index=True)
@@ -100,11 +99,7 @@ def items(self):
 
         line = [coord_i, coord_j + image]
 
-        plotter.plot_path(line,
-                          lattice=lattice,
-                          ax=ax,
-                          color='red',
-                          linewidth=lw)
+        plotter.plot_path(line, lattice=lattice, ax=ax, color='red', linewidth=lw)
 
     lines = ax.lines[3:]
 
@@ -121,29 +116,27 @@ def update(frame_no):
         t_frame = t_start + (frame_no * skip)
 
         for i, event in events.iterrows():
-
             if event['start time'] > t_frame:
                 break
 
-            lw = max(maxwidth - decay * (t_frame - event['start time']),
-                     minwidth)
+            lw = max(maxwidth - decay * (t_frame - event['start time']), minwidth)
 
             line = lines[i]
             line.set_color('red')
             line.set_linewidth(lw)
 
             points[event['start site']].set_facecolor(
-                color_from[event['atom index'] % len(color_from)])
+                color_from[event['atom index'] % len(color_from)]
+            )
             points[event['destination site']].set_facecolor(
-                color_to[event['atom index'] % len(color_to)])
+                color_to[event['atom index'] % len(color_to)]
+            )
 
         start_time = event['start time']
         ax.set_title(f'T: {t_frame} | Next jump: {start_time}')
 
     n_frames = int((t_stop - t_start) / skip)
 
-    return animation.FuncAnimation(fig=fig,
-                                   func=update,
-                                   frames=n_frames,
-                                   interval=interval,
-                                   repeat=False)
+    return animation.FuncAnimation(
+        fig=fig, func=update, frames=n_frames, interval=interval, repeat=False
+    )
diff --git a/src/gemdat/plots/matplotlib/_jumps_vs_distance.py b/src/gemdat/plots/matplotlib/_jumps_vs_distance.py
index ba6b1057..ff05c1ea 100644
--- a/src/gemdat/plots/matplotlib/_jumps_vs_distance.py
+++ b/src/gemdat/plots/matplotlib/_jumps_vs_distance.py
@@ -48,8 +48,6 @@ def jumps_vs_distance(
 
     ax.bar(x, counts, width=(jump_res * 0.8))
 
-    ax.set(title='Jumps vs. Distance',
-           xlabel='Distance (Å)',
-           ylabel='Number of jumps')
+    ax.set(title='Jumps vs. Distance', xlabel='Distance (Å)', ylabel='Number of jumps')
 
     return fig
diff --git a/src/gemdat/plots/matplotlib/_jumps_vs_time.py b/src/gemdat/plots/matplotlib/_jumps_vs_time.py
index d0ac2a21..069dfe11 100644
--- a/src/gemdat/plots/matplotlib/_jumps_vs_time.py
+++ b/src/gemdat/plots/matplotlib/_jumps_vs_time.py
@@ -34,8 +34,6 @@ def jumps_vs_time(*, jumps: Jumps, binsize: int = 500) -> plt.Figure:
 
     ax.hist(jumps.data['stop time'], bins=bins, width=0.8 * binsize)
 
-    ax.set(title='Jumps vs. time',
-           xlabel='Time (steps)',
-           ylabel='Number of jumps')
+    ax.set(title='Jumps vs. time', xlabel='Time (steps)', ylabel='Number of jumps')
 
     return fig
diff --git a/src/gemdat/plots/matplotlib/_msd_per_element.py b/src/gemdat/plots/matplotlib/_msd_per_element.py
index 2655cfbe..870ee135 100644
--- a/src/gemdat/plots/matplotlib/_msd_per_element.py
+++ b/src/gemdat/plots/matplotlib/_msd_per_element.py
@@ -53,16 +53,20 @@ def msd_per_element(
                 ax.plot(t_values, traj, lw=0.1, c=last_color, label=label)
 
         if show_shaded:
-            ax.fill_between(t_values,
-                            msd_mean - msd_std,
-                            msd_mean + msd_std,
-                            color=last_color,
-                            alpha=0.2,
-                            label=f'{sp.symbol} std')
+            ax.fill_between(
+                t_values,
+                msd_mean - msd_std,
+                msd_mean + msd_std,
+                color=last_color,
+                alpha=0.2,
+                label=f'{sp.symbol} std',
+            )
 
     ax.legend()
-    ax.set(title='Mean squared displacement per element',
-           xlabel='Time lag (ps)',
-           ylabel='MSD (Å$^2$)')
+    ax.set(
+        title='Mean squared displacement per element',
+        xlabel='Time lag (ps)',
+        ylabel='MSD (Å$^2$)',
+    )
 
     return fig
diff --git a/src/gemdat/plots/matplotlib/_radial_distribution.py b/src/gemdat/plots/matplotlib/_radial_distribution.py
index fc75d37b..7de79559 100644
--- a/src/gemdat/plots/matplotlib/_radial_distribution.py
+++ b/src/gemdat/plots/matplotlib/_radial_distribution.py
@@ -29,8 +29,10 @@ def radial_distribution(rdfs: Iterable[RDFData]) -> plt.Figure:
     states = ', '.join({rdf.state for rdf in rdfs})
 
     ax.legend()
-    ax.set(title=f'Radial distribution function ({states})',
-           xlabel='Distance (Å)',
-           ylabel='Counts')
+    ax.set(
+        title=f'Radial distribution function ({states})',
+        xlabel='Distance (Å)',
+        ylabel='Counts',
+    )
 
     return fig
diff --git a/src/gemdat/plots/matplotlib/_rectilinear.py b/src/gemdat/plots/matplotlib/_rectilinear.py
index 9d08800a..7b3a5040 100644
--- a/src/gemdat/plots/matplotlib/_rectilinear.py
+++ b/src/gemdat/plots/matplotlib/_rectilinear.py
@@ -9,10 +9,12 @@
 )
 
 
-def rectilinear(*,
-                orientations: Orientations,
-                shape: tuple[int, int] = (90, 360),
-                normalize_histo: bool = True) -> plt.Figure:
+def rectilinear(
+    *,
+    orientations: Orientations,
+    shape: tuple[int, int] = (90, 360),
+    normalize_histo: bool = True,
+) -> plt.Figure:
     """Plot a rectilinear projection of a spherical function. This function
     uses the transformed trajectory.
 
diff --git a/src/gemdat/plots/matplotlib/_shape.py b/src/gemdat/plots/matplotlib/_shape.py
index 46c02d2c..adc15c3d 100644
--- a/src/gemdat/plots/matplotlib/_shape.py
+++ b/src/gemdat/plots/matplotlib/_shape.py
@@ -37,11 +37,13 @@ def shape(
     x_labels = ('X / Å', 'Y / Å', 'Z / Å')
     y_labels = ('Y / Å', 'Z / Å', 'X / Å')
 
-    fig, axes = plt.subplots(nrows=2,
-                             ncols=3,
-                             sharex=True,
-                             figsize=(12, 5),
-                             gridspec_kw={'height_ratios': (4, 1)})
+    fig, axes = plt.subplots(
+        nrows=2,
+        ncols=3,
+        sharex=True,
+        figsize=(12, 5),
+        gridspec_kw={'height_ratios': (4, 1)},
+    )
 
     distances = shape.distances()
     distances_sq = distances**2
diff --git a/src/gemdat/plots/matplotlib/_vibrational_amplitudes.py b/src/gemdat/plots/matplotlib/_vibrational_amplitudes.py
index 64e12cd2..c91b665c 100644
--- a/src/gemdat/plots/matplotlib/_vibrational_amplitudes.py
+++ b/src/gemdat/plots/matplotlib/_vibrational_amplitudes.py
@@ -30,8 +30,10 @@ def vibrational_amplitudes(*, trajectory: Trajectory) -> plt.Figure:
     y_gauss = stats.norm.pdf(x, 0, metrics.vibration_amplitude())
     ax.plot(x, y_gauss, 'r')
 
-    ax.set(title='Histogram of vibrational amplitudes with fitted Gaussian',
-           xlabel='Amplitude (Å)',
-           ylabel='Occurrence (a.u.)')
+    ax.set(
+        title='Histogram of vibrational amplitudes with fitted Gaussian',
+        xlabel='Amplitude (Å)',
+        ylabel='Occurrence (a.u.)',
+    )
 
     return fig
diff --git a/src/gemdat/plots/plotly/__init__.py b/src/gemdat/plots/plotly/__init__.py
index 66bf517c..6cafdf25 100644
--- a/src/gemdat/plots/plotly/__init__.py
+++ b/src/gemdat/plots/plotly/__init__.py
@@ -1,4 +1,5 @@
 """This module contains all the plots that Gemdat can generate with plotly."""
+
 from __future__ import annotations
 
 from ._autocorrelation import autocorrelation
diff --git a/src/gemdat/plots/plotly/_autocorrelation.py b/src/gemdat/plots/plotly/_autocorrelation.py
index 60cb2c67..2c764e1f 100644
--- a/src/gemdat/plots/plotly/_autocorrelation.py
+++ b/src/gemdat/plots/plotly/_autocorrelation.py
@@ -36,36 +36,39 @@ def autocorrelation(
     fig = go.Figure()
 
     if show_shaded:
-        error_y = {
-            'type': 'data',
-            'array': ac_std,
-            'width': 0.1,
-            'thickness': 0.1
-        }
+        error_y = {'type': 'data', 'array': ac_std, 'width': 0.1, 'thickness': 0.1}
     else:
         error_y = None
 
     fig.add_trace(
-        go.Scatter(x=t_values,
-                   y=ac_mean,
-                   error_y=error_y,
-                   name='FFT Autocorrelation',
-                   mode='lines',
-                   line={'width': 3},
-                   legendgroup='autocorr'))
+        go.Scatter(
+            x=t_values,
+            y=ac_mean,
+            error_y=error_y,
+            name='FFT Autocorrelation',
+            mode='lines',
+            line={'width': 3},
+            legendgroup='autocorr',
+        )
+    )
 
     if show_traces:
         for i, trace in enumerate(ac):
             fig.add_trace(
-                go.Scatter(x=t_values,
-                           y=trace,
-                           name=i,
-                           mode='lines',
-                           line={'width': 0.25},
-                           showlegend=False))
+                go.Scatter(
+                    x=t_values,
+                    y=trace,
+                    name=i,
+                    mode='lines',
+                    line={'width': 0.25},
+                    showlegend=False,
+                )
+            )
 
-    fig.update_layout(title='FFT Autocorrelation',
-                      xaxis_title='Time lag (ps)',
-                      yaxis_title='mean + std')
+    fig.update_layout(
+        title='FFT Autocorrelation',
+        xaxis_title='Time lag (ps)',
+        yaxis_title='mean + std',
+    )
 
     return fig
diff --git a/src/gemdat/plots/plotly/_bond_length_distribution.py b/src/gemdat/plots/plotly/_bond_length_distribution.py
index d5750858..809bd042 100644
--- a/src/gemdat/plots/plotly/_bond_length_distribution.py
+++ b/src/gemdat/plots/plotly/_bond_length_distribution.py
@@ -8,9 +8,9 @@
 from .._shared import _orientations_to_histogram, _fit_skewnorm_to_hist
 
 
-def bond_length_distribution(*,
-                             orientations: Orientations,
-                             bins: int = 50) -> go.Figure:
+def bond_length_distribution(
+    *, orientations: Orientations, bins: int = 50
+) -> go.Figure:
     """Plot the bond length probability distribution.
 
     Parameters
@@ -35,14 +35,15 @@ def bond_length_distribution(*,
     )
 
     fig.add_trace(
-        go.Scatter(x=x,
-                   y=y,
-                   name='Skewed Gaussian Fit',
-                   mode='lines',
-                   line={'width': 3}))
-
-    fig.update_layout(title='Bond length probability distribution',
-                      xaxis_title='Bond length (Å)',
-                      yaxis_title='Probability density (Å<sup>-1</sup>)')
+        go.Scatter(
+            x=x, y=y, name='Skewed Gaussian Fit', mode='lines', line={'width': 3}
+        )
+    )
+
+    fig.update_layout(
+        title='Bond length probability distribution',
+        xaxis_title='Bond length (Å)',
+        yaxis_title='Probability density (Å<sup>-1</sup>)',
+    )
 
     return fig
diff --git a/src/gemdat/plots/plotly/_collective_jumps.py b/src/gemdat/plots/plotly/_collective_jumps.py
index 6bbc310a..f15f4431 100644
--- a/src/gemdat/plots/plotly/_collective_jumps.py
+++ b/src/gemdat/plots/plotly/_collective_jumps.py
@@ -31,16 +31,10 @@ def collective_jumps(*, jumps: Jumps) -> go.Figure:
 
     ticks = list(range(len(labels)))
 
-    fig.update_layout(xaxis={
-        'tickmode': 'array',
-        'tickvals': ticks,
-        'ticktext': labels
-    },
-                      yaxis={
-                          'tickmode': 'array',
-                          'tickvals': ticks,
-                          'ticktext': labels
-                      },
-                      title='Cooperative jumps per jump-type combination')
+    fig.update_layout(
+        xaxis={'tickmode': 'array', 'tickvals': ticks, 'ticktext': labels},
+        yaxis={'tickmode': 'array', 'tickvals': ticks, 'ticktext': labels},
+        title='Cooperative jumps per jump-type combination',
+    )
 
     return fig
diff --git a/src/gemdat/plots/plotly/_density.py b/src/gemdat/plots/plotly/_density.py
index 8ad1258e..77d5e619 100644
--- a/src/gemdat/plots/plotly/_density.py
+++ b/src/gemdat/plots/plotly/_density.py
@@ -35,4 +35,5 @@ def density(
         Output as plotly figure
     """
     from ._plot3d import plot_3d
+
     return plot_3d(volume=volume, structure=structure, lattice=force_lattice)
diff --git a/src/gemdat/plots/plotly/_displacement_histogram.py b/src/gemdat/plots/plotly/_displacement_histogram.py
index 7b9f7828..7755dca1 100644
--- a/src/gemdat/plots/plotly/_displacement_histogram.py
+++ b/src/gemdat/plots/plotly/_displacement_histogram.py
@@ -22,18 +22,21 @@ def _trajectory_to_dataframe(trajectory: Trajectory) -> pd.DataFrame:
     """
     data = []
     for specie, distance in zip(
-            trajectory.species,
-            trajectory.distances_from_base_position()[:, -1]):
+        trajectory.species, trajectory.distances_from_base_position()[:, -1]
+    ):
         data.append((specie, round(distance)))
 
     df = pd.DataFrame(columns=['Element', 'Displacement'], data=data)
-    df = df.groupby(['Displacement', 'Element'
-                     ]).size().reset_index().rename(columns={0: 'count'})
+    df = (
+        df.groupby(['Displacement', 'Element'])
+        .size()
+        .reset_index()
+        .rename(columns={0: 'count'})
+    )
     return df
 
 
-def displacement_histogram(trajectory: Trajectory,
-                           n_parts: int = 1) -> go.Figure:
+def displacement_histogram(trajectory: Trajectory, n_parts: int = 1) -> go.Figure:
     """Plot histogram of total displacement at final timestep.
 
     Parameters
@@ -51,21 +54,16 @@ def displacement_histogram(trajectory: Trajectory,
     if n_parts == 1:
         df = _trajectory_to_dataframe(trajectory)
 
-        fig = px.bar(df,
-                     x='Displacement',
-                     y='count',
-                     color='Element',
-                     barmode='stack')
+        fig = px.bar(df, x='Displacement', y='count', color='Element', barmode='stack')
 
-        fig.update_layout(title='Displacement per element',
-                          xaxis_title='Displacement (Å)',
-                          yaxis_title='Nr. of atoms')
+        fig.update_layout(
+            title='Displacement per element',
+            xaxis_title='Displacement (Å)',
+            yaxis_title='Nr. of atoms',
+        )
     else:
         interval = np.linspace(0, len(trajectory) - 1, n_parts + 1)
-        dfs = [
-            _trajectory_to_dataframe(part)
-            for part in trajectory.split(n_parts)
-        ]
+        dfs = [_trajectory_to_dataframe(part) for part in trajectory.split(n_parts)]
 
         all_df = pd.concat(dfs)
 
@@ -75,17 +73,19 @@ def displacement_histogram(trajectory: Trajectory,
         std = grouped.std().reset_index().rename(columns={'count': 'std'})
         df = mean.merge(std, how='inner')
 
-        fig = px.bar(df,
-                     x='Displacement',
-                     y='mean',
-                     color='Element',
-                     error_y='std',
-                     barmode='group')
+        fig = px.bar(
+            df,
+            x='Displacement',
+            y='mean',
+            color='Element',
+            error_y='std',
+            barmode='group',
+        )
 
         fig.update_layout(
-            title=
-            f'Displacement per element after {int(interval[1]-interval[0])} timesteps',
+            title=f'Displacement per element after {int(interval[1]-interval[0])} timesteps',
             xaxis_title='Displacement (Å)',
-            yaxis_title='Nr. of atoms')
+            yaxis_title='Nr. of atoms',
+        )
 
     return fig
diff --git a/src/gemdat/plots/plotly/_displacement_per_atom.py b/src/gemdat/plots/plotly/_displacement_per_atom.py
index f80888ba..4ffb492a 100644
--- a/src/gemdat/plots/plotly/_displacement_per_atom.py
+++ b/src/gemdat/plots/plotly/_displacement_per_atom.py
@@ -25,14 +25,15 @@ def displacement_per_atom(*, trajectory: Trajectory) -> go.Figure:
 
     for i, distance in enumerate(distances):
         fig.add_trace(
-            go.Scatter(y=distance,
-                       name=i,
-                       mode='lines',
-                       line={'width': 1},
-                       showlegend=False))
-
-    fig.update_layout(title='Displacement per atom',
-                      xaxis_title='Time step',
-                      yaxis_title='Displacement (Å)')
+            go.Scatter(
+                y=distance, name=i, mode='lines', line={'width': 1}, showlegend=False
+            )
+        )
+
+    fig.update_layout(
+        title='Displacement per atom',
+        xaxis_title='Time step',
+        yaxis_title='Displacement (Å)',
+    )
 
     return fig
diff --git a/src/gemdat/plots/plotly/_displacement_per_element.py b/src/gemdat/plots/plotly/_displacement_per_element.py
index 33aaeeb0..9d62a650 100644
--- a/src/gemdat/plots/plotly/_displacement_per_element.py
+++ b/src/gemdat/plots/plotly/_displacement_per_element.py
@@ -25,29 +25,40 @@ def displacement_per_element(*, trajectory: Trajectory) -> go.Figure:
 
     for symbol, (mean, std) in displacements.items():
         fig.add_trace(
-            go.Scatter(y=mean,
-                       name=symbol + ' + std',
-                       mode='lines',
-                       line={'width': 3},
-                       legendgroup=symbol))
+            go.Scatter(
+                y=mean,
+                name=symbol + ' + std',
+                mode='lines',
+                line={'width': 3},
+                legendgroup=symbol,
+            )
+        )
         fig.add_trace(
-            go.Scatter(y=mean + std,
-                       name=symbol + ' + std',
-                       mode='lines',
-                       line={'width': 0},
-                       legendgroup=symbol,
-                       showlegend=False))
+            go.Scatter(
+                y=mean + std,
+                name=symbol + ' + std',
+                mode='lines',
+                line={'width': 0},
+                legendgroup=symbol,
+                showlegend=False,
+            )
+        )
         fig.add_trace(
-            go.Scatter(y=mean - std,
-                       name=symbol + ' + std',
-                       mode='lines',
-                       line={'width': 0},
-                       legendgroup=symbol,
-                       showlegend=False,
-                       fill='tonexty'))
-
-    fig.update_layout(title='Displacement per element',
-                      xaxis_title='Time step',
-                      yaxis_title='Displacement (Å)')
+            go.Scatter(
+                y=mean - std,
+                name=symbol + ' + std',
+                mode='lines',
+                line={'width': 0},
+                legendgroup=symbol,
+                showlegend=False,
+                fill='tonexty',
+            )
+        )
+
+    fig.update_layout(
+        title='Displacement per element',
+        xaxis_title='Time step',
+        yaxis_title='Displacement (Å)',
+    )
 
     return fig
diff --git a/src/gemdat/plots/plotly/_energy_along_path.py b/src/gemdat/plots/plotly/_energy_along_path.py
index 5ff745ec..cd6af709 100644
--- a/src/gemdat/plots/plotly/_energy_along_path.py
+++ b/src/gemdat/plots/plotly/_energy_along_path.py
@@ -33,11 +33,14 @@ def energy_along_path(
     fig = go.Figure()
 
     fig.add_trace(
-        go.Scatter(x=np.arange(len(path.energy)),
-                   y=path.energy,
-                   name='Optimal path',
-                   mode='lines',
-                   line={'width': 3}))
+        go.Scatter(
+            x=np.arange(len(path.energy)),
+            y=path.energy,
+            name='Optimal path',
+            mode='lines',
+            line={'width': 3},
+        )
+    )
 
     if structure:
         nearest_sites = path.path_over_structure(structure)
@@ -79,14 +82,17 @@ def energy_along_path(
     if other_paths:
         for idx, path in enumerate(other_paths):
             fig.add_trace(
-                go.Scatter(x=np.arange(len(path.energy)),
-                           y=path.energy,
-                           name=f'Alternative {idx+1}',
-                           mode='lines',
-                           line={'width': 1}))
-
-    fig.update_layout(title='Pathway',
-                      xaxis_title='Steps',
-                      yaxis_title='Free energy (eV)')
+                go.Scatter(
+                    x=np.arange(len(path.energy)),
+                    y=path.energy,
+                    name=f'Alternative {idx+1}',
+                    mode='lines',
+                    line={'width': 1},
+                )
+            )
+
+    fig.update_layout(
+        title='Pathway', xaxis_title='Steps', yaxis_title='Free energy (eV)'
+    )
 
     return fig
diff --git a/src/gemdat/plots/plotly/_frequency_vs_occurence.py b/src/gemdat/plots/plotly/_frequency_vs_occurence.py
index 0af7db28..1d0e31bd 100644
--- a/src/gemdat/plots/plotly/_frequency_vs_occurence.py
+++ b/src/gemdat/plots/plotly/_frequency_vs_occurence.py
@@ -38,14 +38,14 @@ def frequency_vs_occurence(*, trajectory: Trajectory) -> go.Figure:
     sum_freqs = np.sum(one_sided, axis=0)
     smoothed = np.convolve(sum_freqs, np.ones(51), 'same') / 51
     fig.add_trace(
-        go.Scatter(y=smoothed,
-                   x=f,
-                   mode='lines',
-                   line={
-                       'width': 3,
-                       'color': 'blue'
-                   },
-                   showlegend=False))
+        go.Scatter(
+            y=smoothed,
+            x=f,
+            mode='lines',
+            line={'width': 3, 'color': 'blue'},
+            showlegend=False,
+        )
+    )
 
     y_max = np.max(sum_freqs)
 
@@ -54,25 +54,23 @@ def frequency_vs_occurence(*, trajectory: Trajectory) -> go.Figure:
     if attempt_freq:
         fig.add_vline(x=attempt_freq, line={'width': 2, 'color': 'red'})
     if attempt_freq and attempt_freq_std:
-        fig.add_vline(x=attempt_freq + attempt_freq_std,
-                      line={
-                          'width': 2,
-                          'color': 'red',
-                          'dash': 'dash'
-                      })
-        fig.add_vline(x=attempt_freq - attempt_freq_std,
-                      line={
-                          'width': 2,
-                          'color': 'red',
-                          'dash': 'dash'
-                      })
-
-    fig.update_layout(title='Frequency vs Occurence',
-                      xaxis_title='Frequency (Hz)',
-                      yaxis_title='Occurrence (a.u.)',
-                      xaxis_range=[-0.1e13, 2.5e13],
-                      yaxis_range=[0, y_max],
-                      width=600,
-                      height=500)
+        fig.add_vline(
+            x=attempt_freq + attempt_freq_std,
+            line={'width': 2, 'color': 'red', 'dash': 'dash'},
+        )
+        fig.add_vline(
+            x=attempt_freq - attempt_freq_std,
+            line={'width': 2, 'color': 'red', 'dash': 'dash'},
+        )
+
+    fig.update_layout(
+        title='Frequency vs Occurence',
+        xaxis_title='Frequency (Hz)',
+        yaxis_title='Occurrence (a.u.)',
+        xaxis_range=[-0.1e13, 2.5e13],
+        yaxis_range=[0, y_max],
+        width=600,
+        height=500,
+    )
 
     return fig
diff --git a/src/gemdat/plots/plotly/_jumps_3d.py b/src/gemdat/plots/plotly/_jumps_3d.py
index fcad5ae5..ddc0f30a 100644
--- a/src/gemdat/plots/plotly/_jumps_3d.py
+++ b/src/gemdat/plots/plotly/_jumps_3d.py
@@ -22,4 +22,5 @@ def jumps_3d(*, jumps: Jumps) -> go.Figure:
         Output figure
     """
     from ._plot3d import plot_3d
+
     return plot_3d(jumps=jumps, structure=jumps.sites)
diff --git a/src/gemdat/plots/plotly/_jumps_vs_distance.py b/src/gemdat/plots/plotly/_jumps_vs_distance.py
index 3570ba77..9fae9033 100644
--- a/src/gemdat/plots/plotly/_jumps_vs_distance.py
+++ b/src/gemdat/plots/plotly/_jumps_vs_distance.py
@@ -11,10 +11,9 @@
     from gemdat import Jumps
 
 
-def jumps_vs_distance(*,
-                      jumps: Jumps,
-                      jump_res: float = 0.1,
-                      n_parts: int = 1) -> go.Figure:
+def jumps_vs_distance(
+    *, jumps: Jumps, jump_res: float = 0.1, n_parts: int = 1
+) -> go.Figure:
     """Plot jumps vs. distance histogram.
 
     Parameters
@@ -45,8 +44,7 @@ def jumps_vs_distance(*,
     data = []
     for transitions_part in jumps.split(n_parts=n_parts):
         counts = np.zeros_like(x)
-        for idx, n in zip(bin_idx.flatten(),
-                          transitions_part.matrix().flatten()):
+        for idx, n in zip(bin_idx.flatten(), transitions_part.matrix().flatten()):
             counts[idx] += n
         for idx in range(n_bins):
             if counts[idx] > 0:
@@ -62,14 +60,12 @@ def jumps_vs_distance(*,
     if n_parts == 1:
         fig = px.bar(df, x='Displacement', y='mean', barmode='stack')
     else:
-        fig = px.bar(df,
-                     x='Displacement',
-                     y='mean',
-                     error_y='std',
-                     barmode='stack')
+        fig = px.bar(df, x='Displacement', y='mean', error_y='std', barmode='stack')
 
-    fig.update_layout(title='Jumps vs. Distance',
-                      xaxis_title='Distance (Å)',
-                      yaxis_title='Number of jumps')
+    fig.update_layout(
+        title='Jumps vs. Distance',
+        xaxis_title='Distance (Å)',
+        yaxis_title='Number of jumps',
+    )
 
     return fig
diff --git a/src/gemdat/plots/plotly/_jumps_vs_time.py b/src/gemdat/plots/plotly/_jumps_vs_time.py
index 733718e1..7dc9e6dd 100644
--- a/src/gemdat/plots/plotly/_jumps_vs_time.py
+++ b/src/gemdat/plots/plotly/_jumps_vs_time.py
@@ -11,10 +11,7 @@
     from gemdat import Jumps
 
 
-def jumps_vs_time(*,
-                  jumps: Jumps,
-                  bins: int = 8,
-                  n_parts: int = 1) -> go.Figure:
+def jumps_vs_time(*, jumps: Jumps, bins: int = 8, n_parts: int = 1) -> go.Figure:
     """Plot jumps vs. distance histogram.
 
     Parameters
@@ -37,9 +34,10 @@ def jumps_vs_time(*,
 
     for jumps_part in jumps.split(n_parts=n_parts):
         data.append(
-            np.histogram(jumps_part.data['start time'],
-                         bins=bins,
-                         range=(0., maxlen))[0])
+            np.histogram(jumps_part.data['start time'], bins=bins, range=(0.0, maxlen))[
+                0
+            ]
+        )
 
     df = pd.DataFrame(data=data)
     columns = [binsize / 2 + binsize * col for col in range(bins)]
@@ -47,17 +45,18 @@ def jumps_vs_time(*,
     mean = [df[col].mean() for col in df.columns]
     std = [df[col].std() for col in df.columns]
 
-    df = pd.DataFrame(data=zip(columns, mean, std),
-                      columns=['time', 'count', 'std'])
+    df = pd.DataFrame(data=zip(columns, mean, std), columns=['time', 'count', 'std'])
 
     if n_parts > 1:
         fig = px.bar(df, x='time', y='count', error_y='std')
     else:
         fig = px.bar(df, x='time', y='count')
 
-    fig.update_layout(bargap=0.2,
-                      title='Jumps vs. time',
-                      xaxis_title='Time (steps)',
-                      yaxis_title='Number of jumps')
+    fig.update_layout(
+        bargap=0.2,
+        title='Jumps vs. time',
+        xaxis_title='Time (steps)',
+        yaxis_title='Number of jumps',
+    )
 
     return fig
diff --git a/src/gemdat/plots/plotly/_msd_per_element.py b/src/gemdat/plots/plotly/_msd_per_element.py
index 6ec26096..f8ef6bea 100644
--- a/src/gemdat/plots/plotly/_msd_per_element.py
+++ b/src/gemdat/plots/plotly/_msd_per_element.py
@@ -35,19 +35,21 @@ def msd_per_element(*, trajectory: Trajectory) -> go.Figure:
         t_values = np.arange(len(msd_mean)) * time_ps
 
         fig.add_trace(
-            go.Scatter(x=t_values,
-                       y=msd_mean,
-                       error_y=dict(type='data',
-                                    array=msd_std,
-                                    width=0.1,
-                                    thickness=0.1),
-                       name=f'{sp.symbol} mean+std',
-                       mode='lines',
-                       line={'width': 3},
-                       legendgroup=sp.symbol))
-
-    fig.update_layout(title='Mean squared displacement per element',
-                      xaxis_title='Time lag (ps)',
-                      yaxis_title=r'MSD (Å<sup>2</sup>)')
+            go.Scatter(
+                x=t_values,
+                y=msd_mean,
+                error_y=dict(type='data', array=msd_std, width=0.1, thickness=0.1),
+                name=f'{sp.symbol} mean+std',
+                mode='lines',
+                line={'width': 3},
+                legendgroup=sp.symbol,
+            )
+        )
+
+    fig.update_layout(
+        title='Mean squared displacement per element',
+        xaxis_title='Time lag (ps)',
+        yaxis_title=r'MSD (Å<sup>2</sup>)',
+    )
 
     return fig
diff --git a/src/gemdat/plots/plotly/_plot3d.py b/src/gemdat/plots/plotly/_plot3d.py
index 0206d218..4a2e60a9 100644
--- a/src/gemdat/plots/plotly/_plot3d.py
+++ b/src/gemdat/plots/plotly/_plot3d.py
@@ -34,10 +34,21 @@ def plot_lattice_vectors(lattice: Lattice, *, fig: go.Figure):
     abc = a + b + c
 
     for i, (start, stop) in enumerate(
-        ((org, a), (org, b), (org, c), (a, a + b), (a, a + c), (b, b + a),
-         (b, b + c), (c, c + a), (c, c + b), (abc, abc - a), (abc, abc - b),
-         (abc, abc - c))):
-
+        (
+            (org, a),
+            (org, b),
+            (org, c),
+            (a, a + b),
+            (a, a + c),
+            (b, b + a),
+            (b, b + c),
+            (c, c + a),
+            (c, c + b),
+            (abc, abc - a),
+            (abc, abc - b),
+            (abc, abc - c),
+        )
+    ):
         color = ('red', 'blue', 'green', 'black')[min(i, 3)]
         width = (5, 5, 5, 1)[min(i, 3)]
 
@@ -48,19 +59,15 @@ def plot_lattice_vectors(lattice: Lattice, *, fig: go.Figure):
                 z=(start[2], stop[2]),
                 mode='lines',
                 name=None,
-                line={
-                    'color': color,
-                    'width': width
-                },
+                line={'color': color, 'width': width},
                 showlegend=False,
-            ))
+            )
+        )
 
 
-def plot_points(points: np.ndarray,
-                labels: Sequence,
-                *,
-                fig: go.Figure,
-                point_size: int = 5):
+def plot_points(
+    points: np.ndarray, labels: Sequence, *, fig: go.Figure, point_size: int = 5
+):
     """Plot points using plotly.
 
     Parameters
@@ -86,25 +93,21 @@ def plot_points(points: np.ndarray,
         color = colors[label]
 
         fig.add_trace(
-            go.Scatter3d(x=[x],
-                         y=[y],
-                         z=[z],
-                         mode='markers',
-                         name=label,
-                         marker={
-                             'size': point_size,
-                             'color': color,
-                             'line': {
-                                 'width': 2.5
-                             }
-                         },
-                         showlegend=False))
-
-
-def plot_structure(structure: Structure,
-                   *,
-                   lattice: Lattice | None = None,
-                   fig: go.Figure):
+            go.Scatter3d(
+                x=[x],
+                y=[y],
+                z=[z],
+                mode='markers',
+                name=label,
+                marker={'size': point_size, 'color': color, 'line': {'width': 2.5}},
+                showlegend=False,
+            )
+        )
+
+
+def plot_structure(
+    structure: Structure, *, lattice: Lattice | None = None, fig: go.Figure
+):
     """Plot structure using plotly.
 
     Parameters
@@ -168,16 +171,19 @@ def plot_volume(
         cart_verts = lattice.get_cartesian_coords(verts)
 
         fig.add_trace(
-            go.Mesh3d(x=cart_verts[:, 0],
-                      y=cart_verts[:, 1],
-                      z=cart_verts[:, 2],
-                      i=faces[:, 0],
-                      j=faces[:, 1],
-                      k=faces[:, 2],
-                      name=f'{isoval=}',
-                      opacity=alphavals[i],
-                      color=colors[i],
-                      showlegend=False))
+            go.Mesh3d(
+                x=cart_verts[:, 0],
+                y=cart_verts[:, 1],
+                z=cart_verts[:, 2],
+                i=faces[:, 0],
+                j=faces[:, 1],
+                k=faces[:, 2],
+                name=f'{isoval=}',
+                opacity=alphavals[i],
+                color=colors[i],
+                showlegend=False,
+            )
+        )
 
 
 def plot_paths(
@@ -198,7 +204,7 @@ def plot_paths(
         Plotly figure to add paths too
     """
     if not isinstance(paths, Collection):
-        paths = (paths, )
+        paths = (paths,)
 
     path, *other_paths = paths
 
@@ -214,8 +220,7 @@ def plot_paths(
             size = 5
             color = None
 
-        x_path, y_path, z_path = lattice.get_cartesian_coords(
-            path.frac_sites()).T
+        x_path, y_path, z_path = lattice.get_cartesian_coords(path.frac_sites()).T
 
         fig.add_trace(
             go.Scatter3d(
@@ -228,10 +233,11 @@ def plot_paths(
                     'size': size,
                     'color': color,
                     'symbol': 'circle',
-                    'opacity': 0.9
+                    'opacity': 0.9,
                 },
                 name=name,
-            ))
+            )
+        )
 
 
 def plot_jumps(jumps: Jumps, *, fig: go.Figure):
@@ -278,14 +284,17 @@ def plot_jumps(jumps: Jumps, *, fig: go.Figure):
                     line_dash='dashdot' if any(image) != 0 else 'solid',
                     line_width=lw * 3,
                     line_color='black',
-                ))
+                )
+            )
 
 
-def update_layout(*,
-                  lattice: Lattice,
-                  fig: go.Figure,
-                  title: str = 'Gemdat 3D plot',
-                  zoom: float = 0.1):
+def update_layout(
+    *,
+    lattice: Lattice,
+    fig: go.Figure,
+    title: str = 'Gemdat 3D plot',
+    zoom: float = 0.1,
+):
     """Update layout, title, scene, etc for figure.
 
     Arguments
@@ -299,51 +308,48 @@ def update_layout(*,
     title : str
         Title of the plot
     """
-    fig.update_layout(title=title,
-                      scene={
-                          'aspectmode': 'manual',
-                          'aspectratio': {
-                              'x': lattice.a * zoom,
-                              'y': lattice.b * zoom,
-                              'z': lattice.c * zoom,
-                          },
-                          'xaxis_title': 'X (Å)',
-                          'yaxis_title': 'Y (Å)',
-                          'zaxis_title': 'Z (Å)'
-                      },
-                      legend={
-                          'orientation': 'h',
-                          'yanchor': 'bottom',
-                          'xanchor': 'left',
-                          'x': 0,
-                          'y': -0.1
-                      },
-                      showlegend=True,
-                      margin={
-                          'l': 0,
-                          'r': 0,
-                          'b': 0,
-                          't': 0
-                      },
-                      scene_camera={
-                          'projection': {
-                              'type': 'orthographic'
-                          },
-                          'eye': {
-                              'x': -lattice.a * 0.05,
-                              'y': -lattice.b * 0.2,
-                              'z': lattice.c * 0.15,
-                          }
-                      })
-
-
-def plot_3d(*,
-            volume: Volume | None = None,
-            structure: Structure | None = None,
-            paths: Pathway | list[Pathway] | None = None,
-            jumps: Jumps | None = None,
-            lattice: Lattice | None = None,
-            title: str = '3D plot') -> go.Figure:
+    fig.update_layout(
+        title=title,
+        scene={
+            'aspectmode': 'manual',
+            'aspectratio': {
+                'x': lattice.a * zoom,
+                'y': lattice.b * zoom,
+                'z': lattice.c * zoom,
+            },
+            'xaxis_title': 'X (Å)',
+            'yaxis_title': 'Y (Å)',
+            'zaxis_title': 'Z (Å)',
+        },
+        legend={
+            'orientation': 'h',
+            'yanchor': 'bottom',
+            'xanchor': 'left',
+            'x': 0,
+            'y': -0.1,
+        },
+        showlegend=True,
+        margin={'l': 0, 'r': 0, 'b': 0, 't': 0},
+        scene_camera={
+            'projection': {'type': 'orthographic'},
+            'eye': {
+                'x': -lattice.a * 0.05,
+                'y': -lattice.b * 0.2,
+                'z': lattice.c * 0.15,
+            },
+        },
+    )
+
+
+def plot_3d(
+    *,
+    volume: Volume | None = None,
+    structure: Structure | None = None,
+    paths: Pathway | list[Pathway] | None = None,
+    jumps: Jumps | None = None,
+    lattice: Lattice | None = None,
+    title: str = '3D plot',
+) -> go.Figure:
     """Plot 3d."""
     fig = go.Figure()
 
diff --git a/src/gemdat/plots/plotly/_radial_distribution.py b/src/gemdat/plots/plotly/_radial_distribution.py
index 58e3d661..382ed269 100644
--- a/src/gemdat/plots/plotly/_radial_distribution.py
+++ b/src/gemdat/plots/plotly/_radial_distribution.py
@@ -31,11 +31,14 @@ def radial_distribution(rdfs: Iterable[RDFData]) -> go.Figure:
                 name=rdf.symbol,
                 mode='lines',
                 # line={'width': 0.25}
-            ))
+            )
+        )
 
     states = ', '.join({rdf.state for rdf in rdfs})
-    fig.update_layout(title=f'Radial distribution function ({states})',
-                      xaxis_title='Distance (Å)',
-                      yaxis_title='Counts')
+    fig.update_layout(
+        title=f'Radial distribution function ({states})',
+        xaxis_title='Distance (Å)',
+        yaxis_title='Counts',
+    )
 
     return fig
diff --git a/src/gemdat/plots/plotly/_rectilinear.py b/src/gemdat/plots/plotly/_rectilinear.py
index 6b1e0980..c16326dc 100644
--- a/src/gemdat/plots/plotly/_rectilinear.py
+++ b/src/gemdat/plots/plotly/_rectilinear.py
@@ -10,10 +10,12 @@
 )
 
 
-def rectilinear(*,
-                orientations: Orientations,
-                shape: tuple[int, int] = (90, 360),
-                normalize_histo: bool = True) -> go.Figure:
+def rectilinear(
+    *,
+    orientations: Orientations,
+    shape: tuple[int, int] = (90, 360),
+    normalize_histo: bool = True,
+) -> go.Figure:
     """Plot a rectilinear projection of a spherical function. This function
     uses the transformed trajectory.
 
@@ -48,13 +50,17 @@ def rectilinear(*,
     phi = np.linspace(0, 360, axis_phi)
     theta = np.linspace(0, 180, axis_theta)
 
-    fig = go.Figure(data=go.Contour(x=phi,
-                                    y=theta,
-                                    z=hist,
-                                    colorbar={
-                                        'title': 'Areal probability',
-                                        'titleside': 'right',
-                                    }))
+    fig = go.Figure(
+        data=go.Contour(
+            x=phi,
+            y=theta,
+            z=hist,
+            colorbar={
+                'title': 'Areal probability',
+                'titleside': 'right',
+            },
+        )
+    )
 
     fig.update_layout(
         title='Rectilinear plot',
diff --git a/src/gemdat/plots/plotly/_vibrational_amplitudes.py b/src/gemdat/plots/plotly/_vibrational_amplitudes.py
index c226e718..c8356143 100644
--- a/src/gemdat/plots/plotly/_vibrational_amplitudes.py
+++ b/src/gemdat/plots/plotly/_vibrational_amplitudes.py
@@ -10,10 +10,9 @@
 from gemdat.trajectory import Trajectory
 
 
-def vibrational_amplitudes(*,
-                           trajectory: Trajectory,
-                           bins: int = 50,
-                           n_parts: int = 1) -> go.Figure:
+def vibrational_amplitudes(
+    *, trajectory: Trajectory, bins: int = 50, n_parts: int = 1
+) -> go.Figure:
     """Plot histogram of vibrational amplitudes with fitted Gaussian.
 
     Parameters
@@ -32,8 +31,7 @@ def vibrational_amplitudes(*,
     trajectories = trajectory.split(n_parts)
     single_metrics = SimulationMetrics(trajectory)
     metrics = [
-        SimulationMetrics(trajectory).amplitudes()
-        for trajectory in trajectories
+        SimulationMetrics(trajectory).amplitudes() for trajectory in trajectories
     ]
 
     max_amp = max(max(metric) for metric in metrics)
@@ -46,26 +44,22 @@ def vibrational_amplitudes(*,
 
     for metric in metrics:
         data.append(
-            np.histogram(metric,
-                         bins=bins,
-                         range=(min_amp, max_amp),
-                         density=True)[0])
+            np.histogram(metric, bins=bins, range=(min_amp, max_amp), density=True)[0]
+        )
 
     df = pd.DataFrame(data=data)
 
     # offset to middle of bar
     offset = (max_amp - min_amp) / (bins * 2)
 
-    columns = np.linspace(min_amp + offset,
-                          max_amp + offset,
-                          bins,
-                          endpoint=False)
+    columns = np.linspace(min_amp + offset, max_amp + offset, bins, endpoint=False)
 
     mean = [df[col].mean() for col in df.columns]
     std = [df[col].std() for col in df.columns]
 
-    df = pd.DataFrame(data=zip(columns, mean, std),
-                      columns=['amplitude', 'count', 'std'])
+    df = pd.DataFrame(
+        data=zip(columns, mean, std), columns=['amplitude', 'count', 'std']
+    )
 
     if n_parts == 1:
         fig = px.bar(df, x='amplitude', y='count')
@@ -79,6 +73,7 @@ def vibrational_amplitudes(*,
     fig.update_layout(
         title='Histogram of vibrational amplitudes with fitted Gaussian',
         xaxis_title='Amplitude (Å)',
-        yaxis_title='Occurrence (a.u.)')
+        yaxis_title='Occurrence (a.u.)',
+    )
 
     return fig
diff --git a/src/gemdat/rdf.py b/src/gemdat/rdf.py
index 167e2fe9..62da52d7 100644
--- a/src/gemdat/rdf.py
+++ b/src/gemdat/rdf.py
@@ -47,8 +47,7 @@ def _get_states(labels: list[str]) -> dict[int, str]:
     return states
 
 
-def _get_states_array(transitions: Transitions,
-                      labels: list[str]) -> np.ndarray:
+def _get_states_array(transitions: Transitions, labels: list[str]) -> np.ndarray:
     """Helper function to generate integer array of transition states."""
     states = _uniqify_labels(transitions.states, labels)
     states_prev = _uniqify_labels(transitions.states_prev(), labels)
@@ -63,9 +62,7 @@ def _get_symbol_indices(structure: Structure) -> dict[str, np.ndarray]:
     """Helper function to generate symbol indices."""
     symbols = structure.symbol_set
     return {
-        symbol:
-        np.argwhere([sp.symbol == symbol
-                     for sp in structure.species]).flatten()
+        symbol: np.argwhere([sp.symbol == symbol for sp in structure.species]).flatten()
         for symbol in symbols
     }
 
@@ -86,6 +83,7 @@ class RDFData:
         State that the floating species is in, e.g.
         the jump that it is making.
     """
+
     x: np.ndarray
     y: np.ndarray
     symbol: str
@@ -94,6 +92,7 @@ class RDFData:
     def plot(self, **kwargs):
         """See [gemdat.plots.radial_distribution][] for more info."""
         from gemdat import plots
+
         return plots.radial_distribution(rdfs=self, **kwargs)
 
 
@@ -138,13 +137,13 @@ def radial_distribution(
     bins = np.arange(0, max_dist + resolution, resolution)
     length = len(bins) + 1
 
-    rdfs: dict[tuple[str, str],
-               np.ndarray] = defaultdict(lambda: np.zeros(length, dtype=int))
+    rdfs: dict[tuple[str, str], np.ndarray] = defaultdict(
+        lambda: np.zeros(length, dtype=int)
+    )
 
     n_steps = len(trajectory)
 
     for i in track(range(n_steps), transient=True):
-
         t_coords = coords[i]
         t_sp_coords = sp_coords[i]
 
@@ -162,8 +161,7 @@ def radial_distribution(
 
             for symbol, symbol_idx in symbol_indices.items():
                 rdf_state = rdf[k_idx, symbol_idx].flatten()
-                rdfs[state_str, symbol] += np.bincount(rdf_state,
-                                                       minlength=length)
+                rdfs[state_str, symbol] += np.bincount(rdf_state, minlength=length)
 
     ret: dict[str, list[RDFData]] = {}
 
diff --git a/src/gemdat/shape.py b/src/gemdat/shape.py
index 8645bd14..7866f05f 100644
--- a/src/gemdat/shape.py
+++ b/src/gemdat/shape.py
@@ -28,6 +28,7 @@ class ShapeData:
     radius : float
         Maximum distance from center (Å)
     """
+
     site: PeriodicSite
     coords: np.ndarray
     radius: float
@@ -75,6 +76,7 @@ def centroid(self) -> np.ndarray:
     def plot(self, **kwargs):
         """See [gemdat.plots.shape][] for more info."""
         from gemdat import plots
+
         return plots.shape(shape=self, **kwargs)
 
 
@@ -88,8 +90,13 @@ class ShapeAnalyzer:
     for performing shape analysis.
     """
 
-    def __init__(self, *, sites: Collection[PeriodicSite], lattice: Lattice,
-                 spacegroup: SpaceGroup):
+    def __init__(
+        self,
+        *,
+        sites: Collection[PeriodicSite],
+        lattice: Lattice,
+        spacegroup: SpaceGroup,
+    ):
         """Set up shape analyzer from a collection of unique periodic sites,
         the lattice, and spacegroup.
 
@@ -107,7 +114,6 @@ def __init__(self, *, sites: Collection[PeriodicSite], lattice: Lattice,
         self.spacegroup = spacegroup
 
     def __repr__(self):
-
         def to_str(x):
             return f'{x:>10.6f}'
 
@@ -118,11 +124,13 @@ def to_str(x):
             symbol = self.spacegroup.symbol
 
         out = [
-            self.__class__.__name__, 'Spacegroup',
-            f'    {symbol} ({self.spacegroup.int_number})', 'Lattice',
+            self.__class__.__name__,
+            'Spacegroup',
+            f'    {symbol} ({self.spacegroup.int_number})',
+            'Lattice',
             f"    abc   : {' '.join(to_str(val) for val in self.lattice.abc)}",
             f"    angles: {' '.join(to_str(val) for val in self.lattice.angles)}",
-            f'Unique sites ({len(self.sites)})'
+            f'Unique sites ({len(self.sites)})',
         ]
         for site in self.sites:
             out.append(f'    {site!r}')
@@ -163,11 +171,9 @@ def from_structure(cls, structure: Structure):
         symmetrized_structure = sga.get_symmetrized_structure()
         return cls.from_symmetrized_structure(structure=symmetrized_structure)
 
-    def find_equivalent_positions(self,
-                                  *,
-                                  site: PeriodicSite,
-                                  positions: np.ndarray,
-                                  radius: float = 1.0) -> np.ndarray:
+    def find_equivalent_positions(
+        self, *, site: PeriodicSite, positions: np.ndarray, radius: float = 1.0
+    ) -> np.ndarray:
         """Cluster all symmetrically equivalent positions within sphere around
         `site`.
 
@@ -210,8 +216,7 @@ def find_equivalent_positions(self,
 
             # digitize differences to move all close positions to
             # same sphere around coordr
-            offsets = np.digitize(close - sym_coords, bins=[0.5, -0.4999999
-                                                            ]) - 1
+            offsets = np.digitize(close - sym_coords, bins=[0.5, -0.4999999]) - 1
             close += offsets
 
             inversed = op.inverse.operate_multi(close)
@@ -224,12 +229,13 @@ def find_equivalent_positions(self,
         cart_coords = self.lattice.get_cartesian_coords(centered)
         return cart_coords
 
-    def analyze_trajectory(self,
-                           trajectory: Trajectory,
-                           *,
-                           supercell: None
-                           | tuple[float, float, float] = None,
-                           radius: float = 1.0) -> list[ShapeData]:
+    def analyze_trajectory(
+        self,
+        trajectory: Trajectory,
+        *,
+        supercell: None | tuple[float, float, float] = None,
+        radius: float = 1.0,
+    ) -> list[ShapeData]:
         """Perform shape analysis on trajectory.
 
         Similar to [analyze_positions()][ShapeAnalyzer.analyze_positions]. Handles
@@ -269,10 +275,9 @@ def analyze_trajectory(self,
 
         return self.analyze_positions(positions=positions, radius=radius)
 
-    def analyze_positions(self,
-                          positions: np.ndarray,
-                          *,
-                          radius: float = 1.0) -> list[ShapeData]:
+    def analyze_positions(
+        self, positions: np.ndarray, *, radius: float = 1.0
+    ) -> list[ShapeData]:
         """Perform shape analysis on array of fractional coordinates.
 
         Parameters
@@ -296,9 +301,9 @@ def analyze_positions(self,
         shapes = []
 
         for site in self.sites:
-            eqv_coords = self.find_equivalent_positions(site=site,
-                                                        positions=positions,
-                                                        radius=radius)
+            eqv_coords = self.find_equivalent_positions(
+                site=site, positions=positions, radius=radius
+            )
 
             shape = ShapeData(
                 site=site,
@@ -310,9 +315,11 @@ def analyze_positions(self,
 
         return shapes
 
-    def shift_sites(self,
-                    vectors: Sequence[None | Sequence[float]],
-                    coords_are_cartesian: bool = True) -> ShapeAnalyzer:
+    def shift_sites(
+        self,
+        vectors: Sequence[None | Sequence[float]],
+        coords_are_cartesian: bool = True,
+    ) -> ShapeAnalyzer:
         """Shift `.unique_sites` by given vectors.
 
         Parameters
@@ -337,20 +344,22 @@ def shift_sites(self,
 
             coords = site.coords if coords_are_cartesian else site.frac_coords
 
-            new_site = PeriodicSite(site.species,
-                                    coords + offset,
-                                    self.lattice,
-                                    coords_are_cartesian=coords_are_cartesian,
-                                    label=site.label)
+            new_site = PeriodicSite(
+                site.species,
+                coords + offset,
+                self.lattice,
+                coords_are_cartesian=coords_are_cartesian,
+                label=site.label,
+            )
             new_sites.append(new_site)
 
-        return ShapeAnalyzer(sites=new_sites,
-                             spacegroup=self.spacegroup,
-                             lattice=self.lattice)
+        return ShapeAnalyzer(
+            sites=new_sites, spacegroup=self.spacegroup, lattice=self.lattice
+        )
 
-    def optimize_sites(self,
-                       shapes: Sequence[ShapeData],
-                       func: None | Callable = None) -> ShapeAnalyzer:
+    def optimize_sites(
+        self, shapes: Sequence[ShapeData], func: None | Callable = None
+    ) -> ShapeAnalyzer:
         """Optimize unique sites from shape objects.
 
         Note: This function does not take into account special positions.
diff --git a/src/gemdat/simulation_metrics.py b/src/gemdat/simulation_metrics.py
index 3e81ff4e..1e32723b 100644
--- a/src/gemdat/simulation_metrics.py
+++ b/src/gemdat/simulation_metrics.py
@@ -89,10 +89,11 @@ def tracer_diffusivity(self, *, dimensions: int = 3) -> FloatWithUnit:
             Tracer diffusivity in $m^2/s$
         """
         distances = self.trajectory.distances_from_base_position()
-        msd = np.mean(distances[:, -1]**2)  # Angstrom^2
+        msd = np.mean(distances[:, -1] ** 2)  # Angstrom^2
 
-        tracer_diff = (msd * angstrom**2) / (2 * dimensions *
-                                             self.trajectory.total_time)
+        tracer_diff = (msd * angstrom**2) / (
+            2 * dimensions * self.trajectory.total_time
+        )
 
         return FloatWithUnit(tracer_diff, 'm^2 s^-1')
 
@@ -134,14 +135,11 @@ def haven_ratio(self, *, dimensions: int = 3) -> float:
         haven_ratio : float
         """
         return self.tracer_diffusivity(
-            dimensions=dimensions) / self.tracer_diffusivity_center_of_mass(
-                dimensions=dimensions)
+            dimensions=dimensions
+        ) / self.tracer_diffusivity_center_of_mass(dimensions=dimensions)
 
     @weak_lru_cache()
-    def tracer_conductivity(self,
-                            *,
-                            z_ion: int,
-                            dimensions: int = 3) -> FloatWithUnit:
+    def tracer_conductivity(self, *, z_ion: int, dimensions: int = 3) -> FloatWithUnit:
         """Return tracer conductivity as S/m.
 
         Defined as: elementary_charge^2 * charge_ion^2 * diffusivity *
@@ -161,8 +159,9 @@ def tracer_conductivity(self,
         """
         temperature = self.trajectory.metadata['temperature']
         tracer_diff = self.tracer_diffusivity(dimensions=dimensions)
-        tracer_conduc = ((elementary_charge**2) * (z_ion**2) * tracer_diff *
-                         self.particle_density()) / (Boltzmann * temperature)
+        tracer_conduc = (
+            (elementary_charge**2) * (z_ion**2) * tracer_diff * self.particle_density()
+        ) / (Boltzmann * temperature)
 
         return FloatWithUnit(tracer_conduc, 'S m^-1')
 
@@ -249,9 +248,7 @@ def __init__(self, trajectories: list[Trajectory]):
         trajectories: list[Trajectory]
             Input trajectories
         """
-        self.metrics = [
-            SimulationMetrics(trajectory) for trajectory in trajectories
-        ]
+        self.metrics = [SimulationMetrics(trajectory) for trajectory in trajectories]
 
     def speed(self) -> tuple[np.ndarray, np.ndarray]:
         """Calculate mean speed and standard deviations.
@@ -282,8 +279,7 @@ def tracer_diffusivity(self, *, dimensions: int) -> u.ufloat:
             Tracer diffusivity in $m^2/s$, mean and standard deviation
         """
         diffusivities = [
-            metric.tracer_diffusivity(dimensions=dimensions)
-            for metric in self.metrics
+            metric.tracer_diffusivity(dimensions=dimensions) for metric in self.metrics
         ]
         mean_diffusivities = FloatWithUnit(np.mean(diffusivities), 'm^2 s^-1')
         std_diffusivities = FloatWithUnit(np.std(diffusivities), 'm^2 s^-1')
@@ -325,8 +321,7 @@ def vibration_amplitude(self) -> u.ufloat:
         """
         vibes = [metric.vibration_amplitude() for metric in self.metrics]
         mean_vibes = FloatWithUnit(np.mean(vibes), 'ang')
-        standard_vibes = FloatWithUnit(np.std(vibes),
-                                       'ang')  # Standard deviation
+        standard_vibes = FloatWithUnit(np.std(vibes), 'ang')  # Standard deviation
         return u.ufloat(mean_vibes, standard_vibes)
 
     def amplitudes(self) -> tuple[np.ndarray, np.ndarray]:
diff --git a/src/gemdat/trajectory.py b/src/gemdat/trajectory.py
index c4430c78..432f3cbe 100644
--- a/src/gemdat/trajectory.py
+++ b/src/gemdat/trajectory.py
@@ -81,10 +81,8 @@ def to_str(x):
             return f'{x:>10.6f}'
 
         outs.append('abc   : ' + ' '.join(to_str(i) for i in base.lattice.abc))
-        outs.append('angles: ' +
-                    ' '.join(to_str(i) for i in base.lattice.angles))
-        outs.append('pbc   : ' +
-                    ' '.join(str(p).rjust(10) for p in base.lattice.pbc))
+        outs.append('angles: ' + ' '.join(to_str(i) for i in base.lattice.angles))
+        outs.append('pbc   : ' + ' '.join(str(p).rjust(10) for p in base.lattice.pbc))
         if base._charge:
             outs.append(f'Overall Charge: {base._charge:+}')
         outs.append(f'Constant lattice ({self.constant_lattice})')
@@ -121,6 +119,7 @@ def to_volume(self, resolution: float = 0.2) -> Volume:
             Output volume
         """
         from gemdat.volume import trajectory_to_volume
+
         return trajectory_to_volume(self, resolution=resolution)
 
     @property
@@ -230,7 +229,8 @@ def from_vasprun(
         except ET.ParseError as e:
             raise Exception(
                 f'Error parsing {xml_file!r}, to parse incomplete data '
-                'adding `exception_on_bad_xml=False` might help') from e
+                'adding `exception_on_bad_xml=False` might help'
+            ) from e
 
         metadata = {'temperature': run.parameters['TEBEG']}
 
@@ -297,17 +297,19 @@ def center_of_mass(self) -> Trajectory:
         """Return trajectory with center of mass for positions."""
         weights = [s.atomic_mass for s in self.species]
 
-        positions_no_pbc = (self.base_positions +
-                            self.cumulative_displacements)
+        positions_no_pbc = self.base_positions + self.cumulative_displacements
 
-        center_of_mass = np.average(positions_no_pbc, axis=1,
-                                    weights=weights).reshape(-1, 1, 3)
+        center_of_mass = np.average(positions_no_pbc, axis=1, weights=weights).reshape(
+            -1, 1, 3
+        )
 
-        return self.__class__(species=['X'],
-                              coords=center_of_mass,
-                              lattice=self.get_lattice(),
-                              metadata=self.metadata,
-                              time_step=self.time_step)
+        return self.__class__(
+            species=['X'],
+            coords=center_of_mass,
+            lattice=self.get_lattice(),
+            metadata=self.metadata,
+            time_step=self.time_step,
+        )
 
     def drift(
         self,
@@ -338,8 +340,7 @@ def drift(
             displacements = self.filter(species=fixed_species).displacements
         elif floating_species:
             species = {
-                sp.symbol
-                for sp in self.species if sp.symbol not in floating_species
+                sp.symbol for sp in self.species if sp.symbol not in floating_species
             }
             displacements = self.filter(species=species).displacements
         else:
@@ -372,16 +373,19 @@ def apply_drift_correction(
         trajectory : Trajectory
             Ouput trajectory with positions corrected for drift
         """
-        drift = self.drift(fixed_species=fixed_species,
-                           floating_species=floating_species)
+        drift = self.drift(
+            fixed_species=fixed_species, floating_species=floating_species
+        )
 
-        return self.__class__(species=self.species,
-                              coords=self.displacements - drift,
-                              lattice=self.get_lattice(),
-                              metadata=self.metadata,
-                              coords_are_displacement=True,
-                              base_positions=self.base_positions,
-                              time_step=self.time_step)
+        return self.__class__(
+            species=self.species,
+            coords=self.displacements - drift,
+            lattice=self.get_lattice(),
+            metadata=self.metadata,
+            coords_are_displacement=True,
+            base_positions=self.base_positions,
+            time_step=self.time_step,
+        )
 
     def filter(self, species: str | Collection[str]) -> Trajectory:
         """Return trajectory with coordinates for given species only.
@@ -403,15 +407,15 @@ def filter(self, species: str | Collection[str]) -> Trajectory:
         new_coords = self.positions[:, idx]
         new_species = list(compress(self.species, idx))
 
-        return self.__class__(species=new_species,
-                              coords=new_coords,
-                              lattice=self.get_lattice(),
-                              metadata=self.metadata,
-                              time_step=self.time_step)
+        return self.__class__(
+            species=new_species,
+            coords=new_coords,
+            lattice=self.get_lattice(),
+            metadata=self.metadata,
+            time_step=self.time_step,
+        )
 
-    def split(self,
-              n_parts: int = 10,
-              equal_parts: bool = False) -> list[Trajectory]:
+    def split(self, n_parts: int = 10, equal_parts: bool = False) -> list[Trajectory]:
         """Split the trajectory in n similar parts.
 
         Parameters
@@ -426,9 +430,7 @@ def split(self,
         List[Trajectory]
         """
         interval = np.linspace(0, len(self) - 1, n_parts + 1, dtype=int)
-        subtrajectories = [
-            self[start:stop] for start, stop in pairwise(interval)
-        ]
+        subtrajectories = [self[start:stop] for start, stop in pairwise(interval)]
 
         if equal_parts:
             minsize = len(self)
@@ -439,9 +441,7 @@ def split(self,
                 minsize = min(minsize, size)
 
             # Trim all subtrajectories
-            subtrajectories = [
-                trajectory[0:minsize] for trajectory in subtrajectories
-            ]
+            subtrajectories = [trajectory[0:minsize] for trajectory in subtrajectories]
 
         return subtrajectories
 
@@ -460,14 +460,13 @@ def mean_squared_displacement(self) -> np.ndarray:
 
         # Autocorrelation term using FFT [https://doi.org/10.1051/sfn/201112010]:
         # - perform FFT, square it, and then perform inverse FFT
-        fft_result = np.fft.ifft(np.abs(np.fft.fft(pos, n=2 * n_times,
-                                                   axis=-2))**2,
-                                 axis=-2)
+        fft_result = np.fft.ifft(
+            np.abs(np.fft.fft(pos, n=2 * n_times, axis=-2)) ** 2, axis=-2
+        )
         # - keep only the first n_times elements
         fft_result = fft_result[:, :n_times, :].real
         # - sum over the coordinates and divide by the corresponding time window
-        S2 = np.sum(fft_result,
-                    axis=-1) / (n_times - np.arange(n_times)[None, :])
+        S2 = np.sum(fft_result, axis=-1) / (n_times - np.arange(n_times)[None, :])
 
         # Compute the sum of squared displacements
         D = np.square(pos).sum(axis=-1)
@@ -477,12 +476,11 @@ def mean_squared_displacement(self) -> np.ndarray:
         # - compute 2* the sum of squared positions
         double_sum_D = 2 * np.sum(D, axis=-1)[:, None]
         # - compute the cumulative sum of the sum of squares of the positions
-        cumsum_D = np.cumsum(np.insert(D[:, 0:-1], 0, 0, axis=-1) +
-                             np.flip(D, axis=-1),
-                             axis=-1)
+        cumsum_D = np.cumsum(
+            np.insert(D[:, 0:-1], 0, 0, axis=-1) + np.flip(D, axis=-1), axis=-1
+        )
         # - compute the first term in the MSD calculation
-        S1 = (double_sum_D - cumsum_D)[:, :-1] / (n_times -
-                                                  np.arange(n_times)[None, :])
+        S1 = (double_sum_D - cumsum_D)[:, :-1] / (n_times - np.arange(n_times)[None, :])
 
         msd = S1 - 2 * S2
         return msd
@@ -516,6 +514,7 @@ def transitions_between_sites(
         transitions: Transitions
         """
         from gemdat.transitions import Transitions
+
         return Transitions.from_trajectory(
             trajectory=self,
             sites=sites,
@@ -527,29 +526,35 @@ def transitions_between_sites(
     def plot_displacement_per_atom(self, **kwargs):
         """See [gemdat.plots.displacement_per_atom][] for more info."""
         from gemdat import plots
+
         return plots.displacement_per_atom(trajectory=self, **kwargs)
 
     def plot_displacement_per_element(self, **kwargs):
         """See [gemdat.plots.displacement_per_element][] for more info."""
         from gemdat import plots
+
         return plots.displacement_per_element(trajectory=self, **kwargs)
 
     def plot_msd_per_element(self, **kwargs):
         """See [gemdat.plots.msd_per_element][] for more info."""
         from gemdat import plots
+
         return plots.msd_per_element(trajectory=self, **kwargs)
 
     def plot_displacement_histogram(self, **kwargs):
         """See [gemdat.plots.displacement_histogram][] for more info."""
         from gemdat import plots
+
         return plots.displacement_histogram(trajectory=self, **kwargs)
 
     def plot_frequency_vs_occurence(self, **kwargs):
         """See [gemdat.plots.frequency_vs_occurence][] for more info."""
         from gemdat import plots
+
         return plots.frequency_vs_occurence(trajectory=self, **kwargs)
 
     def plot_vibrational_amplitudes(self, **kwargs):
         """See [gemdat.plots.vibrational_amplitudes][] for more info."""
         from gemdat import plots
+
         return plots.vibrational_amplitudes(trajectory=self, **kwargs)
diff --git a/src/gemdat/transitions.py b/src/gemdat/transitions.py
index dd9202c0..76006169 100644
--- a/src/gemdat/transitions.py
+++ b/src/gemdat/transitions.py
@@ -80,7 +80,7 @@ def from_trajectory(
         sites: Structure,
         floating_specie: str,
         site_radius: float | dict[str, float] | None = None,
-        site_inner_fraction: float = 1.,
+        site_inner_fraction: float = 1.0,
     ) -> Transitions:
         """Compute transitions between given sites for floating specie.
 
@@ -107,12 +107,14 @@ def from_trajectory(
 
         if site_radius is None:
             vibration_amplitude = SimulationMetrics(
-                diff_trajectory).vibration_amplitude()
+                diff_trajectory
+            ).vibration_amplitude()
 
             site_radius = _compute_site_radius(
                 trajectory=trajectory,
                 sites=sites,
-                vibration_amplitude=vibration_amplitude)
+                vibration_amplitude=vibration_amplitude,
+            )
 
         if isinstance(site_radius, float):
             site_radius = {'': site_radius}
@@ -130,8 +132,9 @@ def from_trajectory(
             site_inner_fraction=site_inner_fraction,
         )
 
-        events = _calculate_transition_events(atom_sites=states,
-                                              atom_inner_sites=inner_states)
+        events = _calculate_transition_events(
+            atom_sites=states, atom_inner_sites=inner_states
+        )
 
         obj = cls(
             sites=sites,
@@ -162,6 +165,7 @@ def jumps(self, minimal_residence: int = 0, **kwargs) -> Jumps:
         """
 
         from gemdat.jumps import Jumps
+
         return Jumps(self, **kwargs)
 
     @property
@@ -236,9 +240,9 @@ def occupancy(self) -> Structure:
         counts = counts / len(states)
         occupancies = dict(zip(unq, counts))
 
-        species = [{
-            site.specie.name: occupancies.get(i, 0)
-        } for i, site in enumerate(sites)]
+        species = [
+            {site.specie.name: occupancies.get(i, 0)} for i, site in enumerate(sites)
+        ]
 
         return Structure(
             lattice=sites.lattice,
@@ -285,8 +289,7 @@ def split(self, n_parts: int = 10) -> list[Transitions]:
         """
         split_states = np.array_split(self.states, n_parts)
         split_inner_states = np.array_split(self.inner_states, n_parts)
-        split_events = _split_transitions_events(self.events, self.n_states,
-                                                 n_parts)
+        split_events = _split_transitions_events(self.events, self.n_states, n_parts)
 
         split_trajectory = self.trajectory.split(n_parts)
         split_diff_trajectory = self.diff_trajectory.split(n_parts)
@@ -302,7 +305,8 @@ def split(self, n_parts: int = 10) -> list[Transitions]:
                     states=split_states[i],
                     inner_states=split_inner_states[i],
                     events=split_events[i],
-                ))
+                )
+            )
 
         return parts
 
@@ -333,15 +337,14 @@ def _calculate_transition_events(
         atom_site, atom_inner_site = site
 
         # Indices when atom jumps in or out of site
-        i, = np.nonzero((atom_site != np.roll(atom_site, shift=-1)))
+        (i,) = np.nonzero((atom_site != np.roll(atom_site, shift=-1)))
 
         # continue if atom does not jump
         if len(i) < 1:
             continue
 
         # Indices when atom jumps in or out of inner site
-        i2, = np.nonzero((atom_inner_site != np.roll(atom_inner_site,
-                                                     shift=-1)))
+        (i2,) = np.nonzero((atom_inner_site != np.roll(atom_inner_site, shift=-1)))
 
         # Drop last event if it is on the last timestep (side effect of np.roll)
         if i[-1] == len(atom_site) - 1:
@@ -352,29 +355,37 @@ def _calculate_transition_events(
         time = np.unique(np.concatenate((i, i2)))
 
         # Select the timestep just before the transition out of the site
-        transitions = np.vstack([
-            np.ones_like(time) * atom_index,
-            atom_site[time],
-            atom_site[time + 1],
-            atom_inner_site[time],
-            atom_inner_site[time + 1],
-            time,
-        ]).T
+        transitions = np.vstack(
+            [
+                np.ones_like(time) * atom_index,
+                atom_site[time],
+                atom_site[time + 1],
+                atom_inner_site[time],
+                atom_inner_site[time + 1],
+                time,
+            ]
+        ).T
 
         events.append(transitions)
 
     events = np.vstack(events)
-    events = pd.DataFrame(data=events,
-                          columns=[
-                              'atom index', 'start site', 'destination site',
-                              'start inner site', 'destination inner site',
-                              'time'
-                          ])
+    events = pd.DataFrame(
+        data=events,
+        columns=[
+            'atom index',
+            'start site',
+            'destination site',
+            'start inner site',
+            'destination inner site',
+            'time',
+        ],
+    )
     return events
 
 
-def _compute_site_radius(trajectory: Trajectory, sites: Structure,
-                         vibration_amplitude: float) -> float:
+def _compute_site_radius(
+    trajectory: Trajectory, sites: Structure, vibration_amplitude: float
+) -> float:
     """Calculate tolerance wihin which atoms are considered to be close to a
     site.
 
@@ -421,7 +432,8 @@ def _compute_site_radius(trajectory: Trajectory, sites: Structure,
 
             raise ValueError(
                 f'Crystallographic sites are too close together (expected: >{site_radius*2:.4f}, '
-                f'got: {min_dist:.4f} for {msg}')
+                f'got: {min_dist:.4f} for {msg}'
+            )
 
     return site_radius
 
@@ -430,7 +442,7 @@ def _calculate_atom_states(
     sites: Structure,
     trajectory: Trajectory,
     site_radius: dict[str, float],
-    site_inner_fraction: float = 1.,
+    site_inner_fraction: float = 1.0,
 ) -> np.ndarray:
     """Calculate nearest site for each atom coordinate in the trajectory.
 
@@ -462,11 +474,11 @@ def _calculate_atom_states(
     def _site_radius_iterator():
         for label, radius in site_radius.items():
             if label:
-                grouped = ((k, site) for k, site in enumerate(sites)
-                           if site.label == label)
+                grouped = (
+                    (k, site) for k, site in enumerate(sites) if site.label == label
+                )
                 key, site_group = zip(*grouped)
-                frac_coords = np.array(
-                    [site.frac_coords for site in site_group])
+                frac_coords = np.array([site.frac_coords for site in site_group])
                 yield frac_coords, np.array(key), radius
             else:
                 yield sites.frac_coords, None, radius
@@ -479,7 +491,8 @@ def _site_radius_iterator():
     traj_cart_coords = lattice.get_cartesian_coords(traj_frac_coords)
 
     periodic_tree: PeriodicKDTree = PeriodicKDTree(
-        box=np.array(lattice.parameters, dtype=np.float32))
+        box=np.array(lattice.parameters, dtype=np.float32)
+    )
     periodic_tree.set_coords(traj_cart_coords, cutoff=cutoff)
 
     shape = trajectory.positions.shape[0:2]
@@ -488,23 +501,21 @@ def _site_radius_iterator():
 
     for coords, key, radius in _site_radius_iterator():
         cart_coords = lattice.get_cartesian_coords(coords)
-        site_index = periodic_tree.search_tree(cart_coords,
-                                               radius * site_inner_fraction)
+        site_index = periodic_tree.search_tree(
+            cart_coords, radius * site_inner_fraction
+        )
 
         siteno, index = site_index.T
 
         if key is not None:
-            siteno = integer_remap(a=siteno,
-                                   key=key,
-                                   palette=np.unique(siteno))
+            siteno = integer_remap(a=siteno, key=key, palette=np.unique(siteno))
 
         atom_sites[index] = siteno
 
     return atom_sites.reshape(shape)
 
 
-def _calculate_transitions_matrix(events: pd.DataFrame,
-                                  n_sites: int) -> np.ndarray:
+def _calculate_transitions_matrix(events: pd.DataFrame, n_sites: int) -> np.ndarray:
     """Convert list of transition events to dense transitions matrix.
 
     Parameters
@@ -520,19 +531,21 @@ def _calculate_transitions_matrix(events: pd.DataFrame,
         Square matrix with number of each transitions
     """
     transitions = np.zeros((n_sites, n_sites), dtype=int)
-    idx, counts = np.unique(events[['start site', 'destination site']],
-                            return_counts=True,
-                            axis=0)
+    idx, counts = np.unique(
+        events[['start site', 'destination site']], return_counts=True, axis=0
+    )
     start_idx, stop_idx = idx.T
     transitions[start_idx, stop_idx] = counts
     return transitions
 
 
-def _split_transitions_events(events: pd.DataFrame,
-                              n_states: int,
-                              n_parts=10,
-                              split_key='time',
-                              dependent_keys='time') -> list[np.ndarray]:
+def _split_transitions_events(
+    events: pd.DataFrame,
+    n_states: int,
+    n_parts=10,
+    split_key='time',
+    dependent_keys='time',
+) -> list[np.ndarray]:
     """Split list of transition events into equal parts in time.
 
     Parameters
@@ -555,14 +568,12 @@ def _split_transitions_events(events: pd.DataFrame,
         The first dimension corresponds to `n_parts`.
     """
     if len(events) < n_parts:
-        raise ValueError(
-            f'Not enough transitions per part to split into {n_parts}')
+        raise ValueError(f'Not enough transitions per part to split into {n_parts}')
 
     bins = np.linspace(0, n_states + 1, n_parts + 1, dtype=int)
 
     parts = [
-        events[(events[split_key] >= start)
-               & (events[split_key] < stop)].copy()
+        events[(events[split_key] >= start) & (events[split_key] < stop)].copy()
         for start, stop in pairwise(bins)
     ]
 
diff --git a/src/gemdat/utils.py b/src/gemdat/utils.py
index cbd32b94..a672e0aa 100644
--- a/src/gemdat/utils.py
+++ b/src/gemdat/utils.py
@@ -1,4 +1,5 @@
 """This module connects generally useful utilties."""
+
 from __future__ import annotations
 
 import warnings
@@ -11,16 +12,21 @@
 from scipy.spatial import cKDTree
 
 # shortcut to test data
-VASPRUN = Path(__file__).parents[
-    2] / 'tests' / 'data' / 'short_simulation' / 'vasprun.xml'
-VASPCACHE_ORIENTATIONS = Path(__file__).parents[
-    2] / 'tests' / 'data' / 'short_simulation' / 'vasprun_rotations.cache'
+VASPRUN = (
+    Path(__file__).parents[2] / 'tests' / 'data' / 'short_simulation' / 'vasprun.xml'
+)
+VASPCACHE_ORIENTATIONS = (
+    Path(__file__).parents[2]
+    / 'tests'
+    / 'data'
+    / 'short_simulation'
+    / 'vasprun_rotations.cache'
+)
 
 DATA = files('gemdat') / 'data'
 
 
-def nearest_structure_reference(
-        structure: Structure) -> tuple[cKDTree, list[int]]:
+def nearest_structure_reference(structure: Structure) -> tuple[cKDTree, list[int]]:
     """Find distance and index of the nearest site of the structure for each
     voxel using a KD-tree.
 
@@ -41,8 +47,7 @@ def nearest_structure_reference(
     periodic_ids: list[int] = []
     images = np.mgrid[-1:2, -1:2, -1:2].reshape(3, -1).T
     for dx, dy, dz in images:
-        periodic_structure.extend(structure.frac_coords +
-                                  np.array([dx, dy, dz]))
+        periodic_structure.extend(structure.frac_coords + np.array([dx, dy, dz]))
 
         # store the id of the site in the original structure
         periodic_ids.extend(range(len(structure.cart_coords)))
@@ -95,9 +100,9 @@ def bfill(arr: np.ndarray, fill_val: int = -1, axis=-1) -> np.ndarray:
     return np.fliplr(ffill(np.fliplr(arr), fill_val=fill_val))
 
 
-def integer_remap(a: np.ndarray,
-                  key: np.ndarray,
-                  palette: np.ndarray | None = None) -> np.ndarray:
+def integer_remap(
+    a: np.ndarray, key: np.ndarray, palette: np.ndarray | None = None
+) -> np.ndarray:
     """Map integers in array `a` from `palette` -> `key`
 
     Parameters
@@ -157,10 +162,12 @@ def meanfreq(x: np.ndarray, fs: float = 1.0) -> np.ndarray:
     return mnfreq
 
 
-def is_lattice_similar(a: Lattice | Structure,
-                       b: Lattice | Structure,
-                       length_tol: float = 0.5,
-                       angle_tol: float = 1.0) -> bool:
+def is_lattice_similar(
+    a: Lattice | Structure,
+    b: Lattice | Structure,
+    length_tol: float = 0.5,
+    angle_tol: float = 1.0,
+) -> bool:
     """Return True if lattices are similar within given tolerance.
 
     Parameters
@@ -197,12 +204,14 @@ def warn_lattice_not_close(a: Lattice, b: Lattice):
     """Raises a userwarning if lattices are not close."""
     if not is_lattice_similar(a, b):
         warnings.warn(
-            'Lattices are not similar.'
-            f'a: {a.parameters}, b: {b.parameters}', UserWarning)
+            'Lattices are not similar.' f'a: {a.parameters}, b: {b.parameters}',
+            UserWarning,
+        )
 
 
-def _cart2sph(x: np.ndarray, y: np.ndarray,
-              z: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+def _cart2sph(
+    x: np.ndarray, y: np.ndarray, z: np.ndarray
+) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
     """Transform cartesian coordinates to spherical coordinates.
 
     Parameters
@@ -229,9 +238,9 @@ def _cart2sph(x: np.ndarray, y: np.ndarray,
     return az, el, r
 
 
-def cartesian_to_spherical(cart_coords: np.ndarray,
-                           *,
-                           degrees: bool = True) -> np.ndarray:
+def cartesian_to_spherical(
+    cart_coords: np.ndarray, *, degrees: bool = True
+) -> np.ndarray:
     """Trajectory from cartesian coordinates to spherical coordinates.
 
     Parameters
diff --git a/src/gemdat/volume.py b/src/gemdat/volume.py
index 518c79dc..d6a20119 100644
--- a/src/gemdat/volume.py
+++ b/src/gemdat/volume.py
@@ -41,6 +41,7 @@ class Volume:
     units : Unit | None
         Optional unit for the data
     """
+
     data: np.ndarray
     lattice: Lattice
     label: str = 'volume'
@@ -95,8 +96,7 @@ def from_volumetric_data(cls, volume: VolumetricData):
             lattice=volume.structure.lattice,
         )
 
-    def voxel_to_cart_coords(self,
-                             voxel: np.ndarray | list[Any]) -> np.ndarray:
+    def voxel_to_cart_coords(self, voxel: np.ndarray | list[Any]) -> np.ndarray:
         """Convert voxel coordinates to cartesian coordinates.
 
         Parameters
@@ -112,8 +112,7 @@ def voxel_to_cart_coords(self,
         frac_coords = self.voxel_to_frac_coords(voxel)
         return self.lattice.get_cartesian_coords(frac_coords)
 
-    def voxel_to_frac_coords(self,
-                             voxel: np.ndarray | list[Any]) -> np.ndarray:
+    def voxel_to_frac_coords(self, voxel: np.ndarray | list[Any]) -> np.ndarray:
         """Convert voxel coordinates to fractional coordinates.
 
         Parameters
@@ -207,11 +206,13 @@ def find_peaks(
 
         return coords[:, 0:3].astype(int)
 
-    def to_vasp_volume(self,
-                       structure: Structure,
-                       *,
-                       filename: str | None = None,
-                       other: list[Volume] | None = None) -> VolumetricData:
+    def to_vasp_volume(
+        self,
+        structure: Structure,
+        *,
+        filename: str | None = None,
+        other: list[Volume] | None = None,
+    ) -> VolumetricData:
         """Convert to vasp volume.
 
         Parameters
@@ -250,8 +251,9 @@ def to_vasp_volume(self,
 
         return vol_vasp
 
-    def _peaks_to_props(self, peaks: np.ndarray,
-                        background_level: float) -> list[RegionProperties]:
+    def _peaks_to_props(
+        self, peaks: np.ndarray, background_level: float
+    ) -> list[RegionProperties]:
         """Segment volume using watershed algorithm.
 
         Return regionprops.
@@ -282,7 +284,7 @@ def _props_to_frac_coords_centroid(
             for axis in (0, 1, 2):
                 dim = self.dims[axis]
                 if prop.image.shape[axis] == dim:
-                    sel = (coords[:, axis] < dim / 2)
+                    sel = coords[:, axis] < dim / 2
                     coords[sel, axis] += dim
 
             centroids.append(coords.mean(axis=0))
@@ -329,8 +331,7 @@ def to_structure(
         if peaks is None:
             peaks = self.find_peaks(**kwargs)
 
-        props = self._peaks_to_props(peaks=peaks,
-                                     background_level=background_level)
+        props = self._peaks_to_props(peaks=peaks, background_level=background_level)
 
         props_to_frac_coords = self._props_to_frac_coords_centroid
 
@@ -338,9 +339,11 @@ def to_structure(
 
         frac_coords = np.mod(frac_coords, 1)
 
-        structure = Structure(lattice=self.lattice,
-                              coords=frac_coords,
-                              species=[specie for _ in frac_coords])
+        structure = Structure(
+            lattice=self.lattice,
+            coords=frac_coords,
+            species=[specie for _ in frac_coords],
+        )
 
         structure.merge_sites(tol=0.1, mode='average')
 
@@ -363,8 +366,11 @@ def get_free_energy(
             Free energy in eV on the voxel grid
         """
         prob = self.probability()
-        free_energy = -temperature * physical_constants[
-            'Boltzmann constant in eV/K'][0] * np.log(prob)
+        free_energy = (
+            -temperature
+            * physical_constants['Boltzmann constant in eV/K'][0]
+            * np.log(prob)
+        )
 
         return FreeEnergyVolume(
             data=np.nan_to_num(free_energy),
@@ -374,22 +380,21 @@ def get_free_energy(
     def plot_3d(self, **kwargs):
         """See [gemdat.plots.plot_3d][] for more info."""
         from gemdat import plots
+
         return plots.plot_3d(volume=self, **kwargs)
 
 
 class FreeEnergyVolume(Volume):
-
     def free_energy_graph(self, **kwargs) -> nx.Graph:
         """Compute the graph of the free energy for networkx functions.
 
         See [gemdat.path.free_energy_graph][] for more info.
         """
         from .path import free_energy_graph
+
         return free_energy_graph(self.data, **kwargs)
 
-    def optimal_path(self,
-                     F_graph: nx.Graph | None = None,
-                     **kwargs) -> Pathway:
+    def optimal_path(self, F_graph: nx.Graph | None = None, **kwargs) -> Pathway:
         """Calculate the shortest cost-effective path using the desired method.
 
         Parameters
@@ -415,9 +420,9 @@ def optimal_path(self,
         path.dims = self.dims
         return path
 
-    def optimal_n_paths(self,
-                        F_graph: nx.Graph | None = None,
-                        **kwargs) -> list[Pathway]:
+    def optimal_n_paths(
+        self, F_graph: nx.Graph | None = None, **kwargs
+    ) -> list[Pathway]:
         """Calculate the n_paths shortest paths between two sites on the graph.
 
         See [gemdat.path.optimal_n_paths][] for more info.
@@ -439,6 +444,7 @@ def optimal_percolating_path(self, **kwargs) -> Pathway | None:
         See [gemdat.path.optimal_percolating_path][] for more info.
         """
         from .path import optimal_percolating_path
+
         return optimal_percolating_path(self, **kwargs)
 
 
@@ -484,11 +490,13 @@ def trajectory_to_volume(
     ybins = np.linspace(y0, y1, ny)[1:]
     zbins = np.linspace(z0, z1, nz)[1:]
 
-    digitized_coords = np.vstack([
-        np.digitize(coords[:, 0], bins=xbins),
-        np.digitize(coords[:, 1], bins=ybins),
-        np.digitize(coords[:, 2], bins=zbins),
-    ]).T
+    digitized_coords = np.vstack(
+        [
+            np.digitize(coords[:, 0], bins=xbins),
+            np.digitize(coords[:, 1], bins=ybins),
+            np.digitize(coords[:, 2], bins=zbins),
+        ]
+    ).T
 
     indices, counts = np.unique(digitized_coords, return_counts=True, axis=0)
     i, j, k = indices.T


Test patch:
diff --git a/tests/conftest.py b/tests/conftest.py
index 70c1853a..fadb5a60 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -10,23 +10,23 @@
 
 @pytest.fixture()
 def trajectory():
-    coords = np.array([
-        [[.2, .0, .0], [.0, .0, .5], [.0, .0, .5], [.0, .0, .5]],
-        [[.4, .0, .0], [.0, .0, .5], [.0, .0, .5], [.0, .0, .5]],
-        [[.6, .0, .0], [.0, .0, .5], [.0, .0, .5], [.0, .0, .5]],
-        [[.8, .0, .0], [.0, .0, .5], [.0, .0, .5], [.0, .0, .5]],
-        [[.1, .0, .0], [.0, .0, .5], [.0, .0, .5], [.0, .0, .5]],
-    ])
+    coords = np.array(
+        [
+            [[0.2, 0.0, 0.0], [0.0, 0.0, 0.5], [0.0, 0.0, 0.5], [0.0, 0.0, 0.5]],
+            [[0.4, 0.0, 0.0], [0.0, 0.0, 0.5], [0.0, 0.0, 0.5], [0.0, 0.0, 0.5]],
+            [[0.6, 0.0, 0.0], [0.0, 0.0, 0.5], [0.0, 0.0, 0.5], [0.0, 0.0, 0.5]],
+            [[0.8, 0.0, 0.0], [0.0, 0.0, 0.5], [0.0, 0.0, 0.5], [0.0, 0.0, 0.5]],
+            [[0.1, 0.0, 0.0], [0.0, 0.0, 0.5], [0.0, 0.0, 0.5], [0.0, 0.0, 0.5]],
+        ]
+    )
 
     return Trajectory(
-        species=[Species('B'),
-                 Species('Si'),
-                 Species('S'),
-                 Species('C')],
+        species=[Species('B'), Species('Si'), Species('S'), Species('C')],
         coords=coords,
         lattice=np.eye(3),
         metadata={'temperature': 123},
-        time_step=1)
+        time_step=1,
+    )
 
 
 @pytest.fixture()
diff --git a/tests/integration/conftest.py b/tests/integration/conftest.py
index c20852c3..c5098cd1 100644
--- a/tests/integration/conftest.py
+++ b/tests/integration/conftest.py
@@ -20,16 +20,20 @@
 def pytest_configure():
     pytest.vaspxml_available = pytest.mark.skipif(
         not VASP_XML.exists(),
-        reason=
-        ('Simulation data from vasprun.xml example is required for this test. '
-         'Run `git submodule init`/`update`, and extract using `tar -C tests/data/short_simulation '
-         '-xjf tests/data/short_simulation/vasprun.xml.bz2`'))
+        reason=(
+            'Simulation data from vasprun.xml example is required for this test. '
+            'Run `git submodule init`/`update`, and extract using `tar -C tests/data/short_simulation '
+            '-xjf tests/data/short_simulation/vasprun.xml.bz2`'
+        ),
+    )
     pytest.vasporicache_available = pytest.mark.skipif(
         not VASP_ORI_CACHE.exists(),
-        reason=
-        ('Simulation data from vasprun_rotations.cache example is required for this test. '
-         'Run `git submodule init`/`update`, and extract using `tar -C tests/data/short_simulation '
-         '-xjf tests/data/short_simulation/vasprun.xml.bz2`'))
+        reason=(
+            'Simulation data from vasprun_rotations.cache example is required for this test. '
+            'Run `git submodule init`/`update`, and extract using `tar -C tests/data/short_simulation '
+            '-xjf tests/data/short_simulation/vasprun.xml.bz2`'
+        ),
+    )
 
 
 @pytest.fixture(scope='module')
@@ -58,8 +62,9 @@ def structure():
 
 @pytest.fixture(scope='module')
 def vasp_transitions(vasp_traj, structure):
-    transitions = vasp_traj.transitions_between_sites(sites=structure,
-                                                      floating_specie='Li')
+    transitions = vasp_traj.transitions_between_sites(
+        sites=structure, floating_specie='Li'
+    )
     return transitions
 
 
@@ -72,8 +77,7 @@ def vasp_jumps(vasp_transitions):
 def vasp_rdf_data(vasp_traj, structure):
     # Shorten trajectory for faster test
     trajectory = vasp_traj[-1000:]
-    transitions = trajectory.transitions_between_sites(structure,
-                                                       floating_specie='Li')
+    transitions = trajectory.transitions_between_sites(structure, floating_specie='Li')
 
     rdfs = radial_distribution(
         transitions=transitions,
@@ -124,8 +128,7 @@ def vasp_F_graph(vasp_path_vol):
 def vasp_orientations(vasp_traj_orientations):
     central_atoms = 'S'
     satellite_atoms = 'O'
-    orientations = Orientations(vasp_traj_orientations, central_atoms,
-                                satellite_atoms)
+    orientations = Orientations(vasp_traj_orientations, central_atoms, satellite_atoms)
 
     return orientations
 
diff --git a/tests/integration/jumps_test.py b/tests/integration/jumps_test.py
index f45251e6..9c73f44a 100644
--- a/tests/integration/jumps_test.py
+++ b/tests/integration/jumps_test.py
@@ -16,23 +16,35 @@ class TestJumps:  # type: ignore
 
     def test_site_inner_fraction(self, vasp_traj, structure):
         transitions = vasp_traj.transitions_between_sites(
-            sites=structure, floating_specie='Li', site_inner_fraction=0.5)
+            sites=structure, floating_specie='Li', site_inner_fraction=0.5
+        )
         jumps = Jumps(transitions=transitions, minimal_residence=100)
 
         assert len(jumps.data) == 267
-        assert np.all(jumps.data[::100].to_numpy() == np.array([[
-            0, 94, 0, 282, 303
-        ], [15, 74, 8, 1271, 1286], [34, 49, 33, 3141, 3296]]))
+        assert np.all(
+            jumps.data[::100].to_numpy()
+            == np.array(
+                [
+                    [0, 94, 0, 282, 303],
+                    [15, 74, 8, 1271, 1286],
+                    [34, 49, 33, 3141, 3296],
+                ]
+            )
+        )
 
     def test_minimal_residency(self, vasp_traj, structure):
-        transitions = vasp_traj.transitions_between_sites(sites=structure,
-                                                          floating_specie='Li')
+        transitions = vasp_traj.transitions_between_sites(
+            sites=structure, floating_specie='Li'
+        )
         jumps = Jumps(transitions=transitions, minimal_residence=3)
 
         assert len(jumps.data) == 462
-        assert np.all(jumps.data[::200].to_numpy() == np.array([[
-            0, 94, 0, 282, 284
-        ], [18, 54, 24, 2937, 3015], [41, 41, 67, 849, 851]]))
+        assert np.all(
+            jumps.data[::200].to_numpy()
+            == np.array(
+                [[0, 94, 0, 282, 284], [18, 54, 24, 2937, 3015], [41, 41, 67, 849, 851]]
+            )
+        )
 
     def test_n_jumps(self, vasp_jumps):
         assert vasp_jumps.n_solo_jumps == 450
@@ -61,16 +73,13 @@ def test_activation_energies(self, vasp_jumps):
         assert isclose(row['std'], 0.00405952, abs_tol=1e-6)
 
     def test_jump_diffusivity(self, vasp_jumps):
-        assert isclose(vasp_jumps.jump_diffusivity(3),
-                       9.484382e-09,
-                       rel_tol=1e-6)
+        assert isclose(vasp_jumps.jump_diffusivity(3), 9.484382e-09, rel_tol=1e-6)
 
     def test_correlation_factor(self, vasp_traj, vasp_jumps):
         vasp_diff_traj = vasp_traj.filter('Li')
         metrics = SimulationMetrics(vasp_diff_traj)
         tracer_diff = metrics.tracer_diffusivity(dimensions=3)
-        correlation_factor = tracer_diff / vasp_jumps.jump_diffusivity(
-            dimensions=3)
+        correlation_factor = tracer_diff / vasp_jumps.jump_diffusivity(dimensions=3)
         assert isclose(correlation_factor, 0.165600, rel_tol=1e-6)
 
     def test_collective(self, vasp_jumps):
@@ -86,9 +95,14 @@ def test_coll_jumps(self, vasp_jumps):
         coll_jumps = collective.coll_jumps
         assert len(coll_jumps) == 6
         assert vasp_jumps.n_jumps == collective.n_solo_jumps + collective.n_coll_jumps
-        assert coll_jumps == [((49, 31), (33, 49)), ((54, 38), (24, 54)),
-                              ((42, 75), (64, 42)), ((60, 36), (68, 60)),
-                              ((36, 52), (80, 36)), ((2, 58), (92, 2))]
+        assert coll_jumps == [
+            ((49, 31), (33, 49)),
+            ((54, 38), (24, 54)),
+            ((42, 75), (64, 42)),
+            ((60, 36), (68, 60)),
+            ((36, 52), (80, 36)),
+            ((2, 58), (92, 2)),
+        ]
 
     def test_collective_matrix(self, vasp_jumps):
         collective = vasp_jumps.collective()
@@ -102,7 +116,16 @@ def test_multiple_collective(self, vasp_jumps):
         assert jumps.shape == (6, 2)
         assert counts.sum() == 6
         assert np.all(
-            jumps == np.array([[(2, 58), (92, 2)], [(24, 54), (
-                54, 38)], [(33, 49), (49, 31)], [(36, 52), (
-                    80, 36)], [(42, 75), (64, 42)], [(60, 36), (68, 60)]],
-                              dtype=[('start', '<i8'), ('stop', '<i8')]))
+            jumps
+            == np.array(
+                [
+                    [(2, 58), (92, 2)],
+                    [(24, 54), (54, 38)],
+                    [(33, 49), (49, 31)],
+                    [(36, 52), (80, 36)],
+                    [(42, 75), (64, 42)],
+                    [(60, 36), (68, 60)],
+                ],
+                dtype=[('start', '<i8'), ('stop', '<i8')],
+            )
+        )
diff --git a/tests/integration/orientations_test.py b/tests/integration/orientations_test.py
index 8d1e3be3..dbf70760 100644
--- a/tests/integration/orientations_test.py
+++ b/tests/integration/orientations_test.py
@@ -8,8 +8,12 @@
 from gemdat.orientations import calculate_spherical_areas
 
 matrix = np.array(
-    [[1 / 2**0.5, -1 / 6**0.5, 1 / 3**0.5],
-     [1 / 2**0.5, 1 / 6**0.5, -1 / 3**0.5], [0, 2 / 6**0.5, 1 / 3**0.5]], )
+    [
+        [1 / 2**0.5, -1 / 6**0.5, 1 / 3**0.5],
+        [1 / 2**0.5, 1 / 6**0.5, -1 / 3**0.5],
+        [0, 2 / 6**0.5, 1 / 3**0.5],
+    ],
+)
 
 
 @pytest.vasporicache_available  # type: ignore
@@ -17,8 +21,7 @@ def test_normalize(vasp_orientations):
     ret = vasp_orientations.normalize()
     assert isclose(ret.vectors.std(), 0.5773501929401034)
     first = ret.vectors[0, 0]
-    np.testing.assert_allclose(first, [0.300804, 0.433833, -0.849297],
-                               atol=1e-06)
+    np.testing.assert_allclose(first, [0.300804, 0.433833, -0.849297], atol=1e-06)
 
 
 @pytest.vasporicache_available  # type: ignore
@@ -26,8 +29,7 @@ def test_conventional(vasp_orientations):
     ret = vasp_orientations.transform(matrix=matrix)
     assert isclose(ret.vectors.std(), 0.8668510720176622)
     first = ret.vectors[0, 0]
-    np.testing.assert_allclose(first, [-0.712261, 1.378551, -0.213198],
-                               atol=1e-06)
+    np.testing.assert_allclose(first, [-0.712261, 1.378551, -0.213198], atol=1e-06)
 
 
 @pytest.vasporicache_available  # type: ignore
@@ -35,8 +37,7 @@ def test_symmetrize(vasp_orientations):
     ret = vasp_orientations.symmetrize(sym_group='m-3m')
     assert isclose(ret.vectors.std(), 0.8668511628167984)
     first = ret.vectors[0, 0]
-    np.testing.assert_allclose(first, [1.330221, -0.679495, -0.471138],
-                               atol=1e-06)
+    np.testing.assert_allclose(first, [1.330221, -0.679495, -0.471138], atol=1e-06)
 
 
 @pytest.vasporicache_available  # type: ignore
@@ -44,34 +45,36 @@ def test_normalize_symmetrize(vasp_orientations):
     ret = vasp_orientations.normalize().symmetrize(sym_group='m-3m')
     assert isclose(ret.vectors.std(), 0.577350269189626)
     first = ret.vectors[0, 0]
-    np.testing.assert_allclose(first, [0.849297, -0.433833, -0.300804],
-                               atol=1e-06)
+    np.testing.assert_allclose(first, [0.849297, -0.433833, -0.300804], atol=1e-06)
 
 
 @pytest.vasporicache_available  # type: ignore
 def test_normalize_conventional_symmetrize(vasp_orientations):
-    ret = vasp_orientations.normalize().transform(matrix=matrix).symmetrize(
-        sym_group='m-3m')
+    ret = (
+        vasp_orientations.normalize()
+        .transform(matrix=matrix)
+        .symmetrize(sym_group='m-3m')
+    )
     assert isclose(ret.vectors.std(), 0.577350269189626)
     first = ret.vectors[0, 0]
-    np.testing.assert_allclose(first, [0.136119, -0.880154, 0.454753],
-                               atol=1e-06)
+    np.testing.assert_allclose(first, [0.136119, -0.880154, 0.454753], atol=1e-06)
 
 
 @pytest.vasporicache_available  # type: ignore
 def test_conventional_normalize_symmetrize(vasp_orientations):
-    ret = vasp_orientations.transform(matrix=matrix).normalize().symmetrize(
-        sym_group='m-3m')
+    ret = (
+        vasp_orientations.transform(matrix=matrix)
+        .normalize()
+        .symmetrize(sym_group='m-3m')
+    )
     assert isclose(ret.vectors.std(), 0.577350269189626)
     first = ret.vectors[0, 0]
-    np.testing.assert_allclose(first, [0.136119, -0.880154, 0.454753],
-                               atol=1e-06)
+    np.testing.assert_allclose(first, [0.136119, -0.880154, 0.454753], atol=1e-06)
 
 
 @pytest.vasporicache_available  # type: ignore
 def test_fractional_coordinates(vasp_orientations):
-    direction = vasp_orientations._fractional_directions(
-        vasp_orientations._distances)
+    direction = vasp_orientations._fractional_directions(vasp_orientations._distances)
 
     assert isclose(direction.mean(), -2.8363494021277384e-05)
     assert isclose(direction[3, 12, 2], -0.10902666)
@@ -81,8 +84,7 @@ def test_fractional_coordinates(vasp_orientations):
 def test_matching_matrix(vasp_orientations):
     distance = np.array([[1.0, 2.0, 1.5, 0.5, 3.0], [2.0, 1.0, 0.5, 1.5, 3.0]])
     frac_coord_cent = np.array([[[0.3, 0.11, 0.78], [1.6, 1.0, 2.3]]])
-    matching_matrix = vasp_orientations._matching_matrix(
-        distance, frac_coord_cent)
+    matching_matrix = vasp_orientations._matching_matrix(distance, frac_coord_cent)
 
     assert isclose(matching_matrix[-1, -1], 2)
     assert isclose(matching_matrix.sum(), 20)
diff --git a/tests/integration/path_test.py b/tests/integration/path_test.py
index 380c90a7..7fd110a4 100644
--- a/tests/integration/path_test.py
+++ b/tests/integration/path_test.py
@@ -70,11 +70,13 @@ def test_nearest_structure_reference(vasp_path_vol, vasp_path):
 
 @pytest.vaspxml_available  # type: ignore
 def test_optimal_n_paths(vasp_F_graph):
-    paths = optimal_n_paths(F_graph=vasp_F_graph,
-                            start=(10, 4, 13),
-                            stop=(21, 3, 10),
-                            n_paths=3,
-                            min_diff=0.1)
+    paths = optimal_n_paths(
+        F_graph=vasp_F_graph,
+        start=(10, 4, 13),
+        stop=(21, 3, 10),
+        n_paths=3,
+        min_diff=0.1,
+    )
 
     assert len(paths) == 3
     assert isclose(sum(paths[-1].energy), 5.351758190646607)
diff --git a/tests/integration/plot_mpl_test.py b/tests/integration/plot_mpl_test.py
index 99d4c384..023d3b78 100644
--- a/tests/integration/plot_mpl_test.py
+++ b/tests/integration/plot_mpl_test.py
@@ -61,9 +61,13 @@ def test_jumps_3d_animation(vasp_jumps):
     plots.jumps_3d_animation(jumps=vasp_jumps, t_start=1000, t_stop=1001)
 
 
-@image_comparison2(baseline_images=[
-    'radial_distribution_1', 'radial_distribution_2', 'radial_distribution_3'
-])
+@image_comparison2(
+    baseline_images=[
+        'radial_distribution_1',
+        'radial_distribution_2',
+        'radial_distribution_3',
+    ]
+)
 def test_radial_distribution(vasp_rdf_data):
     assert len(vasp_rdf_data) == 3
     for rdfs in vasp_rdf_data.values():
@@ -91,8 +95,12 @@ def test_energy_along_path(vasp_path):
 @image_comparison2(baseline_images=['rectilinear'])
 def test_rectilinear(vasp_orientations):
     matrix = np.array(
-        [[1 / 2**0.5, -1 / 6**0.5, 1 / 3**0.5],
-         [1 / 2**0.5, 1 / 6**0.5, -1 / 3**0.5], [0, 2 / 6**0.5, 1 / 3**0.5]], )
+        [
+            [1 / 2**0.5, -1 / 6**0.5, 1 / 3**0.5],
+            [1 / 2**0.5, 1 / 6**0.5, -1 / 3**0.5],
+            [0, 2 / 6**0.5, 1 / 3**0.5],
+        ],
+    )
 
     orientations = vasp_orientations.normalize().transform(matrix=matrix)
     plots.rectilinear(orientations=orientations, normalize_histo=False)
diff --git a/tests/integration/plot_plotly_test.py b/tests/integration/plot_plotly_test.py
index 29f3f091..f3b95d7c 100644
--- a/tests/integration/plot_plotly_test.py
+++ b/tests/integration/plot_plotly_test.py
@@ -12,11 +12,7 @@
 from matplotlib.testing.compare import compare_images
 
 
-def assert_image_similar(fig,
-                         *,
-                         name: str,
-                         ext: str = 'png',
-                         rms: float = 0.0):
+def assert_image_similar(fig, *, name: str, ext: str = 'png', rms: float = 0.0):
     # Ensure same font is used on different machines (local/CI)
     fig.update_layout(
         font_family='Arial',
@@ -26,8 +22,7 @@ def assert_image_similar(fig,
     RESULTS_DIR = Path() / 'result_images' / Path(__file__).stem
     RESULTS_DIR.mkdir(exist_ok=True, parents=True)
 
-    EXPECTED_DIR = Path(__file__).parent / 'baseline_images' / Path(
-        __file__).stem
+    EXPECTED_DIR = Path(__file__).parent / 'baseline_images' / Path(__file__).stem
 
     filename = f'{name}.{ext}'
 
@@ -42,17 +37,19 @@ def assert_image_similar(fig,
 
     expected_link.symlink_to(expected)
 
-    err: dict[str, Any] = compare_images(expected=str(expected_link),
-                                         actual=str(actual),
-                                         tol=rms,
-                                         in_decorator=True)  # type: ignore
+    err: dict[str, Any] = compare_images(
+        expected=str(expected_link), actual=str(actual), tol=rms, in_decorator=True
+    )  # type: ignore
 
     if err:
         for key in ('actual', 'expected', 'diff'):
             err[key] = Path(err[key]).relative_to('.')
         raise AssertionError(
-            ('images not close (RMS {rms:.3f}):'
-             '\n\t{actual}\n\t{expected}\n\t{diff}'.format(**err)))
+            (
+                'images not close (RMS {rms:.3f}):'
+                '\n\t{actual}\n\t{expected}\n\t{diff}'.format(**err)
+            )
+        )
 
 
 def test_displacement_per_element(vasp_traj):
@@ -145,8 +142,12 @@ def test_energy_along_path(vasp_path):
 
 def test_rectilinear(vasp_orientations):
     matrix = np.array(
-        [[1 / 2**0.5, -1 / 6**0.5, 1 / 3**0.5],
-         [1 / 2**0.5, 1 / 6**0.5, -1 / 3**0.5], [0, 2 / 6**0.5, 1 / 3**0.5]], )
+        [
+            [1 / 2**0.5, -1 / 6**0.5, 1 / 3**0.5],
+            [1 / 2**0.5, 1 / 6**0.5, -1 / 3**0.5],
+            [0, 2 / 6**0.5, 1 / 3**0.5],
+        ],
+    )
 
     orientations = vasp_orientations.normalize().transform(matrix=matrix)
     fig = plots.rectilinear(orientations=orientations, normalize_histo=False)
@@ -155,8 +156,7 @@ def test_rectilinear(vasp_orientations):
 
 
 def test_bond_length_distribution(vasp_orientations):
-    fig = plots.bond_length_distribution(orientations=vasp_orientations,
-                                         bins=50)
+    fig = plots.bond_length_distribution(orientations=vasp_orientations, bins=50)
 
     assert_image_similar(fig, name='bond_length_distribution', rms=0.5)
 
diff --git a/tests/integration/trajectory_test.py b/tests/integration/trajectory_test.py
index 90d6e0fd..5b2e87db 100644
--- a/tests/integration/trajectory_test.py
+++ b/tests/integration/trajectory_test.py
@@ -93,8 +93,10 @@ def test_vibration_metrics(vasp_traj):
     assert speed.shape == (48, 3750)
     assert np.allclose(
         speed[::24, ::1000],
-        [[0., 0.01348809, -0.02682779, 0.00293911],
-         [0., -0.0059066, 0.01617687, 0.01080066]],
+        [
+            [0.0, 0.01348809, -0.02682779, 0.00293911],
+            [0.0, -0.0059066, 0.01617687, 0.01080066],
+        ],
     )
 
     attempt_freq, attempt_freq_std = metrics.attempt_frequency()
@@ -103,12 +105,17 @@ def test_vibration_metrics(vasp_traj):
     assert np.isclose(attempt_freq_std, 857338723421.6494)
 
     amplitudes = metrics.amplitudes()
-    assert amplitudes.shape == (6644, )
+    assert amplitudes.shape == (6644,)
     assert np.allclose(
         amplitudes[::1000],
         [
-            0.29206348, -0.11529727, 0.33333244, -1.05528599, 0.53086991,
-            -0.48080325, 0.61573288
+            0.29206348,
+            -0.11529727,
+            0.33333244,
+            -1.05528599,
+            0.53086991,
+            -0.48080325,
+            0.61573288,
         ],
     )
 
diff --git a/tests/integration/transitions_test.py b/tests/integration/transitions_test.py
index 12fb2356..2a0a0988 100644
--- a/tests/integration/transitions_test.py
+++ b/tests/integration/transitions_test.py
@@ -3,6 +3,7 @@
 
 VASP_XML=/home/stef/md-analysis-matlab-example/vasprun.xml pytest
 """
+
 from __future__ import annotations
 
 from math import isclose
@@ -40,7 +41,8 @@ def test_site_radius_dict(self, vasp_traj, structure):
 
         site_radius = {'A': 0.5, 'B': 0.6, 'C': 0.7}
         transitions = vasp_traj.transitions_between_sites(
-            sites=sites, floating_specie='Li', site_radius=site_radius)
+            sites=sites, floating_specie='Li', site_radius=site_radius
+        )
 
         assert transitions.states.sum() == 3445344
 
@@ -53,20 +55,23 @@ def test_atom_sites(self, vasp_transitions):
         states = vasp_transitions.states
         assert states.shape == (n_steps, n_diffusing)
         assert states.sum() == 6154859
-        assert_allclose(states[slice_],
-                        np.array([[94, -1], [94, -1], [0, 65], [0, 65]]))
+        assert_allclose(
+            states[slice_], np.array([[94, -1], [94, -1], [0, 65], [0, 65]])
+        )
 
         states_next = vasp_transitions.states_next()
         assert states_next.shape == (n_steps, n_diffusing)
         assert states_next.sum() == 8172006
-        assert_allclose(states_next[slice_],
-                        np.array([[94, 1], [94, 65], [0, 65], [0, 65]]))
+        assert_allclose(
+            states_next[slice_], np.array([[94, 1], [94, 65], [0, 65], [0, 65]])
+        )
 
         states_prev = vasp_transitions.states_prev()
         assert states_prev.shape == (n_steps, n_diffusing)
         assert states_prev.sum() == 8148552
-        assert_allclose(states_prev[slice_],
-                        np.array([[94, -1], [94, 65], [0, 65], [0, 65]]))
+        assert_allclose(
+            states_prev[slice_], np.array([[94, -1], [94, 65], [0, 65], [0, 65]])
+        )
 
     def test_n_floating(self, vasp_transitions):
         # https://github.com/GEMDAT-repos/GEMDAT/issues/252
@@ -83,8 +88,14 @@ def test_all_transitions_events(self, vasp_transitions):
         assert vasp_transitions.events.shape == (2105, 6)
         assert_allclose(
             events[::1000],
-            np.array([[0, 94, -1, 94, -1, 202], [22, 70, -1, 70, -1, 3472],
-                      [45, 10, -1, 10, -1, 2386]]))
+            np.array(
+                [
+                    [0, 94, -1, 94, -1, 202],
+                    [22, 70, -1, 70, -1, 3472],
+                    [45, 10, -1, 10, -1, 2386],
+                ]
+            ),
+        )
 
     def test_all_transitions_matrix(self, vasp_transitions):
         matrix = vasp_transitions.matrix()
@@ -92,19 +103,18 @@ def test_all_transitions_matrix(self, vasp_transitions):
         assert matrix.shape == (n_sites, n_sites)
         assert matrix.sum() == 2099
         assert_allclose(
-            np.argwhere(matrix)[::50],
-            np.array([[0, 95], [49, 95], [95, 5], [95, 55]]))
+            np.argwhere(matrix)[::50], np.array([[0, 95], [49, 95], [95, 5], [95, 55]])
+        )
 
     def test_transitions_parts(self, vasp_transitions):
         n_sites = vasp_transitions.n_sites
-        tp = np.stack(
-            [part.matrix() for part in vasp_transitions.split(n_parts=10)])
+        tp = np.stack([part.matrix() for part in vasp_transitions.split(n_parts=10)])
 
         assert tp.shape == (self.n_parts, n_sites, n_sites)
         assert tp.sum() == 2100
         assert_allclose(
-            np.argwhere(tp)[::500],
-            np.array([[0, 0, 95], [4, 23, 95], [8, 28, 95]]))
+            np.argwhere(tp)[::500], np.array([[0, 0, 95], [4, 23, 95], [8, 28, 95]])
+        )
 
     def test_occupancy(self, vasp_transitions):
         structure = vasp_transitions.occupancy()
@@ -118,14 +128,15 @@ def test_occupancy_parts(self, vasp_transitions):
         structures = [part.occupancy() for part in parts]
 
         values1 = [structure[0].species.num_atoms for structure in structures]
-        assert values1 == [
-            0.11066666666666666, 0.708, 0.408, 0.62, 0.42533333333333334
-        ]
+        assert values1 == [0.11066666666666666, 0.708, 0.408, 0.62, 0.42533333333333334]
 
         values2 = [structure.composition.num_atoms for structure in structures]
         assert values2 == [
-            36.653333333333336, 37.129333333333335, 37.98933333333333,
-            35.49199999999999, 35.43733333333334
+            36.653333333333336,
+            37.129333333333335,
+            37.98933333333333,
+            35.49199999999999,
+            35.43733333333334,
         ]
 
     def test_atom_locations(self, vasp_transitions):
@@ -137,19 +148,9 @@ def test_atom_locations_parts(self, vasp_transitions):
         dcts = [part.atom_locations() for part in parts]
 
         assert dcts == [
-            {
-                '48h': 0.7636111111111111
-            },
-            {
-                '48h': 0.7735277777777778
-            },
-            {
-                '48h': 0.7914444444444443
-            },
-            {
-                '48h': 0.7394166666666665
-            },
-            {
-                '48h': 0.7382777777777779
-            },
+            {'48h': 0.7636111111111111},
+            {'48h': 0.7735277777777778},
+            {'48h': 0.7914444444444443},
+            {'48h': 0.7394166666666665},
+            {'48h': 0.7382777777777779},
         ]
diff --git a/tests/orientations_test.py b/tests/orientations_test.py
index c353b10f..6b7ca94c 100644
--- a/tests/orientations_test.py
+++ b/tests/orientations_test.py
@@ -14,9 +14,9 @@
 
 
 def test_orientations_init(trajectory):
-    orientations = Orientations(trajectory=trajectory,
-                                center_type='B',
-                                satellite_type='Si')
+    orientations = Orientations(
+        trajectory=trajectory, center_type='B', satellite_type='Si'
+    )
 
     assert isinstance(orientations, Orientations)
     assert orientations.center_type == 'B'
@@ -25,40 +25,48 @@ def test_orientations_init(trajectory):
 
 
 def test_normalize(trajectory):
-    orientations = Orientations(trajectory=trajectory,
-                                center_type='B',
-                                satellite_type='Si',
-                                in_vectors=np.array([[1, 2, 2], [2, 2, 1]],
-                                                    dtype=float))
+    orientations = Orientations(
+        trajectory=trajectory,
+        center_type='B',
+        satellite_type='Si',
+        in_vectors=np.array([[1, 2, 2], [2, 2, 1]], dtype=float),
+    )
     ret = orientations.normalize()
-    assert_allclose(ret.vectors,
-                    np.array([[1 / 3, 2 / 3, 2 / 3], [2 / 3, 2 / 3, 1 / 3]]))
+    assert_allclose(
+        ret.vectors, np.array([[1 / 3, 2 / 3, 2 / 3], [2 / 3, 2 / 3, 1 / 3]])
+    )
 
 
 def test_conventional(trajectory):
-    orientations = Orientations(trajectory=trajectory,
-                                center_type='B',
-                                satellite_type='Si',
-                                in_vectors=np.array(
-                                    [[1, 0, 0], [0, 1, 0], [0, 0, 1]],
-                                    dtype=float))
+    orientations = Orientations(
+        trajectory=trajectory,
+        center_type='B',
+        satellite_type='Si',
+        in_vectors=np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]], dtype=float),
+    )
     matrix = np.eye(3) * [1, 2, 3]
     ret = orientations.transform(matrix=matrix)
     assert_allclose(ret.vectors, np.array([[1, 0, 0], [0, 2, 0], [0, 0, 3]]))
 
 
 def test_symmetrize(trajectory):
-    orientations = Orientations(trajectory=trajectory,
-                                center_type='B',
-                                satellite_type='Si',
-                                in_vectors=np.array([[[1, 0, 0]], [[0, 1, 0]]],
-                                                    dtype=float))
+    orientations = Orientations(
+        trajectory=trajectory,
+        center_type='B',
+        satellite_type='Si',
+        in_vectors=np.array([[[1, 0, 0]], [[0, 1, 0]]], dtype=float),
+    )
     sym_ops = np.array([[0, -1, 0], [1, 0, 0], [0, 0, -1]])
     ret = orientations.symmetrize(sym_ops=sym_ops)
     assert_allclose(
         ret.vectors,
-        np.array([[[0., 1., 0.], [-1., 0., 0.], [0., 0., -1.]],
-                  [[0., 1., 0.], [-1., 0., 0.], [0., 0., -1.]]]))
+        np.array(
+            [
+                [[0.0, 1.0, 0.0], [-1.0, 0.0, 0.0], [0.0, 0.0, -1.0]],
+                [[0.0, 1.0, 0.0], [-1.0, 0.0, 0.0], [0.0, 0.0, -1.0]],
+            ]
+        ),
+    )
 
 
 def test_orientations(orientations):
@@ -84,5 +92,7 @@ def test_fft_autocorrelation(trajectory):
     autocorr = fft_autocorrelation(trajectory.positions)
     assert isinstance(autocorr, np.ndarray)
     assert isclose(autocorr.mean(), 0.8142314269325723)
-    assert autocorr.shape == (trajectory.positions.shape[1],
-                              trajectory.positions.shape[0])
+    assert autocorr.shape == (
+        trajectory.positions.shape[1],
+        trajectory.positions.shape[0],
+    )
diff --git a/tests/plots_tests.py b/tests/plots_tests.py
index 45306ae5..b6ad3442 100644
--- a/tests/plots_tests.py
+++ b/tests/plots_tests.py
@@ -1,8 +1,10 @@
 def test_matplotlib_imports():
     from gemdat.plots import matplotlib
+
     assert matplotlib.__all__
 
 
 def test_plotly_imports():
     from gemdat.plots import plotly
+
     assert plotly.__all__
diff --git a/tests/segmentation_test.py b/tests/segmentation_test.py
index dd4104c4..cfceeae9 100644
--- a/tests/segmentation_test.py
+++ b/tests/segmentation_test.py
@@ -33,17 +33,19 @@ def test_watershed_horizontal():
         mask=mask,
     )
 
-    expected = np.array([
-        [0, 0, 0, 0, 0, 0, 0],
-        [0, 0, 0, 0, 0, 1, 1],
-        [1, 1, 1, 0, 0, 1, 1],
-        [1, 1, 1, 0, 0, 1, 1],
-        [1, 1, 1, 0, 0, 1, 1],
-        [1, 1, 1, 0, 0, 1, 1],
-        [1, 1, 1, 0, 0, 0, 0],
-        [0, 0, 0, 0, 0, 0, 0],
-        [0, 0, 0, 0, 0, 0, 0],
-    ])
+    expected = np.array(
+        [
+            [0, 0, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 1, 1],
+            [1, 1, 1, 0, 0, 1, 1],
+            [1, 1, 1, 0, 0, 1, 1],
+            [1, 1, 1, 0, 0, 1, 1],
+            [1, 1, 1, 0, 0, 1, 1],
+            [1, 1, 1, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0],
+        ]
+    )
 
     np.testing.assert_allclose(out, expected)
 
@@ -70,14 +72,16 @@ def test_watershed_vertical():
 
     out = watershed_pbc(data, markers, mask=mask)
 
-    expected = np.array([
-        [0, 0, 1, 1, 1, 1, 1, 0, 0],
-        [0, 0, 1, 1, 1, 1, 1, 0, 0],
-        [0, 0, 1, 1, 1, 1, 1, 0, 0],
-        [0, 0, 0, 0, 0, 0, 0, 0, 0],
-        [0, 0, 0, 0, 0, 0, 0, 0, 0],
-        [0, 1, 1, 1, 1, 1, 0, 0, 0],
-        [0, 1, 1, 1, 1, 1, 0, 0, 0],
-    ])
+    expected = np.array(
+        [
+            [0, 0, 1, 1, 1, 1, 1, 0, 0],
+            [0, 0, 1, 1, 1, 1, 1, 0, 0],
+            [0, 0, 1, 1, 1, 1, 1, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 0, 0],
+            [0, 1, 1, 1, 1, 1, 0, 0, 0],
+            [0, 1, 1, 1, 1, 1, 0, 0, 0],
+        ]
+    )
 
     np.testing.assert_allclose(out, expected)
diff --git a/tests/shape_test.py b/tests/shape_test.py
index bf6d629a..7c7f6519 100644
--- a/tests/shape_test.py
+++ b/tests/shape_test.py
@@ -18,11 +18,13 @@ def site():
 
 @pytest.fixture
 def shape(site):
-    coords = np.array([
-        [0, 3, 4],
-        [5, 0, 12],
-        [8, 15, 0],
-    ])
+    coords = np.array(
+        [
+            [0, 3, 4],
+            [5, 0, 12],
+            [8, 15, 0],
+        ]
+    )
     return ShapeData(site=site, coords=coords, radius=20)
 
 
@@ -55,10 +57,12 @@ def test_shape_analyzer(shape_analyzer):
 
 def test_shape_analyzer_analyze_positions(shape_analyzer):
     shapes = shape_analyzer.analyze_positions(
-        np.array([
-            (0.2, 0.2, 0.2),
-            (0.8, 0.8, 0.8),
-        ]),
+        np.array(
+            [
+                (0.2, 0.2, 0.2),
+                (0.8, 0.8, 0.8),
+            ]
+        ),
         radius=2,
     )
 
@@ -110,7 +114,7 @@ def test_shape_analyzer_optimize_sites(shape_analyzer):
     assert_allclose(site.frac_coords, (0.1, 0.2, 0.3))
 
     shape = ShapeData(site=site, coords=np.array([[2, 2, 2]]), radius=1)
-    shifted = shape_analyzer.optimize_sites((shape, ))
+    shifted = shape_analyzer.optimize_sites((shape,))
 
     site = shifted.sites[0]
     assert_allclose(site.frac_coords, (0.3, 0.4, 0.5))
diff --git a/tests/simulation_metrics_test.py b/tests/simulation_metrics_test.py
index fd15d56b..a6ffd1a5 100644
--- a/tests/simulation_metrics_test.py
+++ b/tests/simulation_metrics_test.py
@@ -9,11 +9,10 @@ def test_tracer_metrics(trajectory):
     diff_trajectory = trajectory.filter('B')
     metrics = SimulationMetrics(diff_trajectory)
 
-    assert (np.isclose(metrics.particle_density(), 1e30))
-    assert (np.isclose(metrics.mol_per_liter(), 1660.53906))
-    assert (np.isclose(metrics.tracer_diffusivity(dimensions=3), 2.666667e-10))
-    assert (np.isclose(metrics.tracer_conductivity(z_ion=1, dimensions=3),
-                       4.081278e-09))
+    assert np.isclose(metrics.particle_density(), 1e30)
+    assert np.isclose(metrics.mol_per_liter(), 1660.53906)
+    assert np.isclose(metrics.tracer_diffusivity(dimensions=3), 2.666667e-10)
+    assert np.isclose(metrics.tracer_conductivity(z_ion=1, dimensions=3), 4.081278e-09)
 
 
 def test_tracer_metrics_std(trajectory):
@@ -33,7 +32,7 @@ def test_vibration_metrics(trajectory):
     diff_trajectory = trajectory.filter('B')
     metrics = SimulationMetrics(diff_trajectory)
 
-    assert np.allclose(metrics.speed(), [[0., 0.2, 0.2, 0.2, 0.3]])
+    assert np.allclose(metrics.speed(), [[0.0, 0.2, 0.2, 0.2, 0.3]])
 
     attempt_freq, attempt_freq_std = metrics.attempt_frequency()
 
@@ -48,7 +47,7 @@ def test_vibration_metrics_std(trajectory):
     metrics = SimulationMetricsStd(diff_trajectory.split(2, equal_parts=True))
 
     speed = metrics.speed()
-    assert np.allclose(speed[0], [0., 0.2])
+    assert np.allclose(speed[0], [0.0, 0.2])
 
     amplitudes_mean, amplitudes_std = metrics.amplitudes()
     assert np.isclose(amplitudes_mean, np.array([0.2]))
diff --git a/tests/sites_test.py b/tests/sites_test.py
index 290247e6..9bea5e9f 100644
--- a/tests/sites_test.py
+++ b/tests/sites_test.py
@@ -7,21 +7,22 @@
 
 
 def test_transitions_matrix():
-    a = np.array([
-        [0, 0, 1, 10],
-        [0, 1, 0, 20],
-        [1, 2, 1, 5],
-        [1, 1, 0, 6],
-        [2, 2, 0, 30],
-    ])
+    a = np.array(
+        [
+            [0, 0, 1, 10],
+            [0, 1, 0, 20],
+            [1, 2, 1, 5],
+            [1, 1, 0, 6],
+            [2, 2, 0, 30],
+        ]
+    )
     a = pd.DataFrame(
-        data=a,
-        columns=['atom index', 'start site', 'destination site', 'time'])
+        data=a, columns=['atom index', 'start site', 'destination site', 'time']
+    )
 
     n = 3
     transitions = _calculate_transitions_matrix(a, n_sites=n)
 
     assert transitions.shape == (n, n)
     assert transitions.dtype == int
-    np.testing.assert_equal(transitions,
-                            np.array([[0, 1, 0], [2, 0, 0], [1, 1, 0]]))
+    np.testing.assert_equal(transitions, np.array([[0, 1, 0], [2, 0, 0], [1, 1, 0]]))
diff --git a/tests/trajectory_test.py b/tests/trajectory_test.py
index af3f8830..dc4b8e42 100644
--- a/tests/trajectory_test.py
+++ b/tests/trajectory_test.py
@@ -10,8 +10,10 @@
 def test_trajectory(trajectory):
     assert isinstance(trajectory, Trajectory)
     assert trajectory.species == [
-        Species('B'), Species('Si'),
-        Species('S'), Species('C')
+        Species('B'),
+        Species('Si'),
+        Species('S'),
+        Species('C'),
     ]
     assert trajectory.positions.shape == (5, 4, 3)
     assert trajectory.metadata == {'temperature': 123}
@@ -29,7 +31,7 @@ def test_slice(trajectory):
 def test_filter(trajectory):
     t = trajectory.filter('C')
     assert t.species == [Species('C')]
-    assert np.all(t.positions == [.0, .0, .5])
+    assert np.all(t.positions == [0.0, 0.0, 0.5])
 
 
 def test_get_lattice(trajectory):
@@ -61,13 +63,16 @@ def test_displacements_property(trajectory):
     trajectory = trajectory.filter(['B', 'C'])
     trajectory.to_positions()
 
-    assert_allclose(trajectory.displacements, [
-        [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],
-        [[0.2, 0.0, 0.0], [0.0, 0.0, 0.0]],
-        [[0.2, 0.0, 0.0], [0.0, 0.0, 0.0]],
-        [[0.2, 0.0, 0.0], [0.0, 0.0, 0.0]],
-        [[0.3, 0.0, 0.0], [0.0, 0.0, 0.0]],
-    ])
+    assert_allclose(
+        trajectory.displacements,
+        [
+            [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],
+            [[0.2, 0.0, 0.0], [0.0, 0.0, 0.0]],
+            [[0.2, 0.0, 0.0], [0.0, 0.0, 0.0]],
+            [[0.2, 0.0, 0.0], [0.0, 0.0, 0.0]],
+            [[0.3, 0.0, 0.0], [0.0, 0.0, 0.0]],
+        ],
+    )
 
     assert trajectory.coords_are_displacement
 
@@ -76,13 +81,15 @@ def test_positions_property(trajectory):
     trajectory.to_displacements()
 
     assert_allclose(
-        trajectory.filter(['B', 'S']).positions, [
+        trajectory.filter(['B', 'S']).positions,
+        [
             [[0.2, 0.0, 0.0], [0.0, 0.0, 0.5]],
             [[0.4, 0.0, 0.0], [0.0, 0.0, 0.5]],
             [[0.6, 0.0, 0.0], [0.0, 0.0, 0.5]],
             [[0.8, 0.0, 0.0], [0.0, 0.0, 0.5]],
             [[0.1, 0.0, 0.0], [0.0, 0.0, 0.5]],
-        ])
+        ],
+    )
 
     assert not trajectory.coords_are_displacement
 
@@ -102,21 +109,27 @@ def test_drift_correction(trajectory):
 
 def test_distances_from_base_position(trajectory):
     distances = trajectory.filter(['B', 'Si']).distances_from_base_position()
-    assert_allclose(distances, [
-        [0.0, 0.2, 0.4, 0.6, 0.9],
-        [0.0, 0.0, 0.0, 0.0, 0.0],
-    ])
+    assert_allclose(
+        distances,
+        [
+            [0.0, 0.2, 0.4, 0.6, 0.9],
+            [0.0, 0.0, 0.0, 0.0, 0.0],
+        ],
+    )
 
 
 def test_cumulative_displacements(trajectory):
     displacements = trajectory.filter(['B', 'C']).cumulative_displacements
-    assert_allclose(displacements, [
-        [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],
-        [[0.2, 0.0, 0.0], [0.0, 0.0, 0.0]],
-        [[0.4, 0.0, 0.0], [0.0, 0.0, 0.0]],
-        [[0.6, 0.0, 0.0], [0.0, 0.0, 0.0]],
-        [[0.9, 0.0, 0.0], [0.0, 0.0, 0.0]],
-    ])
+    assert_allclose(
+        displacements,
+        [
+            [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],
+            [[0.2, 0.0, 0.0], [0.0, 0.0, 0.0]],
+            [[0.4, 0.0, 0.0], [0.0, 0.0, 0.0]],
+            [[0.6, 0.0, 0.0], [0.0, 0.0, 0.0]],
+            [[0.9, 0.0, 0.0], [0.0, 0.0, 0.0]],
+        ],
+    )
 
 
 def test_filter_similar_names(trajectory):
@@ -133,8 +146,10 @@ def test_trajectory_extend(trajectory):
 
     assert isinstance(trajectory, Trajectory)
     assert len(trajectory) == 10
-    assert_allclose(trajectory.positions[:, 0, 0],
-                    [0.2, 0.4, 0.6, 0.8, 0.1, 0.2, 0.4, 0.6, 0.8, 0.1])
+    assert_allclose(
+        trajectory.positions[:, 0, 0],
+        [0.2, 0.4, 0.6, 0.8, 0.1, 0.2, 0.4, 0.6, 0.8, 0.1],
+    )
 
 
 def test_mean_squared_displacement(trajectory):
diff --git a/tests/utils_test.py b/tests/utils_test.py
index 32a77fea..9494d866 100644
--- a/tests/utils_test.py
+++ b/tests/utils_test.py
@@ -9,53 +9,63 @@
 
 @pytest.fixture
 def arr():
-    return np.array([
-        [5, -1, -1, 7, 2],
-        [3, -1, 1, 8, -1],
-        [4, 9, 6, -1, -1],
-    ])
+    return np.array(
+        [
+            [5, -1, -1, 7, 2],
+            [3, -1, 1, 8, -1],
+            [4, 9, 6, -1, -1],
+        ]
+    )
 
 
 def test_ffill(arr):
     ret = ffill(arr)
-    expected = np.array([
-        [5, 5, 5, 7, 2],
-        [3, 3, 1, 8, 8],
-        [4, 9, 6, 6, 6],
-    ])
+    expected = np.array(
+        [
+            [5, 5, 5, 7, 2],
+            [3, 3, 1, 8, 8],
+            [4, 9, 6, 6, 6],
+        ]
+    )
 
     assert_equal(ret, expected)
 
 
 def test_bfill(arr):
     ret = bfill(arr)
-    expected = np.array([
-        [5, 7, 7, 7, 2],
-        [3, 1, 1, 8, -1],
-        [4, 9, 6, -1, -1],
-    ])
+    expected = np.array(
+        [
+            [5, 7, 7, 7, 2],
+            [3, 1, 1, 8, -1],
+            [4, 9, 6, -1, -1],
+        ]
+    )
 
     assert_equal(ret, expected)
 
 
 def test_ffill_axis0(arr):
     ret = ffill(arr, axis=0)
-    expected = np.array([
-        [5, -1, -1, 7, 2],
-        [3, -1, 1, 8, 2],
-        [4, 9, 6, 8, 2],
-    ])
+    expected = np.array(
+        [
+            [5, -1, -1, 7, 2],
+            [3, -1, 1, 8, 2],
+            [4, 9, 6, 8, 2],
+        ]
+    )
 
     assert_equal(ret, expected)
 
 
 def test_bfill_axis0(arr):
     ret = bfill(arr, axis=0)
-    expected = np.array([
-        [5, 9, 1, 7, 2],
-        [3, 9, 1, 8, -1],
-        [4, 9, 6, -1, -1],
-    ])
+    expected = np.array(
+        [
+            [5, 9, 1, 7, 2],
+            [3, 9, 1, 8, -1],
+            [4, 9, 6, -1, -1],
+        ]
+    )
 
     assert_equal(ret, expected)
 
@@ -77,11 +87,13 @@ def test_meanfreq_single_timestep():
 
 
 def test_meanfreq():
-    x = np.array([
-        np.sin(np.linspace(0, 1, 6)),
-        np.sin(np.linspace(0, 2, 6)),
-        np.sin(np.linspace(0, 3, 6)),
-    ])
+    x = np.array(
+        [
+            np.sin(np.linspace(0, 1, 6)),
+            np.sin(np.linspace(0, 2, 6)),
+            np.sin(np.linspace(0, 3, 6)),
+        ]
+    )
     ret = meanfreq(x)
 
     expected = np.array([[0.2303359], [0.21308077], [0.17074241]])
@@ -96,20 +108,27 @@ def test_cartesian_to_spherical():
     xyz = np.expand_dims(np.stack((x, y, z), axis=-1), axis=0)
     # test in radians
     ret = cartesian_to_spherical(xyz, degrees=False)
-    expected = np.array([[[0, np.arccos(1 / np.sqrt(2)),
-                           np.sqrt(2)],
-                          [np.pi / 2,
-                           np.arccos(1 / np.sqrt(2)),
-                           np.sqrt(2)]]])
+    expected = np.array(
+        [
+            [
+                [0, np.arccos(1 / np.sqrt(2)), np.sqrt(2)],
+                [np.pi / 2, np.arccos(1 / np.sqrt(2)), np.sqrt(2)],
+            ]
+        ]
+    )
     assert_allclose(ret, expected, rtol=1e-5)
 
     ret = cartesian_to_spherical(xyz, degrees=True)
     expected = np.array(
-        [[[0, np.degrees(np.arccos(1 / np.sqrt(2))),
-           np.sqrt(2)],
-          [
-              np.degrees(np.pi / 2),
-              np.degrees(np.arccos(1 / np.sqrt(2))),
-              np.sqrt(2)
-          ]]])
+        [
+            [
+                [0, np.degrees(np.arccos(1 / np.sqrt(2))), np.sqrt(2)],
+                [
+                    np.degrees(np.pi / 2),
+                    np.degrees(np.arccos(1 / np.sqrt(2))),
+                    np.sqrt(2),
+                ],
+            ]
+        ]
+    )
     assert_allclose(ret, expected, rtol=1e-5)
diff --git a/tests/volume_test.py b/tests/volume_test.py
index ae46068e..a19c7e48 100644
--- a/tests/volume_test.py
+++ b/tests/volume_test.py
@@ -35,9 +35,9 @@ def test_volume_voxel_to_frac(volume, voxel, expected):
 
 TEST_DATA_FRAC_TO_VOX = (
     # coord, expected
-    ((0., 0., 0.), (0, 0, 0)),
+    ((0.0, 0.0, 0.0), (0, 0, 0)),
     ((0.5, 0.5, 0.5), (1, 2, 2)),
-    ((1., 1., 1.), (3, 4, 5)),
+    ((1.0, 1.0, 1.0), (3, 4, 5)),
 )
 
 
@@ -58,8 +58,8 @@ def test_volume_site_to_vox(volume, coord, expected):
 TEST_DATA_VOX_TO_CART = (
     # voxel, expected
     ((-0.5, -0.5, -0.5), (0, 0, 0)),
-    ((0, 0, 0), (1., 1., 1.)),
-    ((3, 4, 5), (7., 9., 11.)),
+    ((0, 0, 0), (1.0, 1.0, 1.0)),
+    ((3, 4, 5), (7.0, 9.0, 11.0)),
 )
 
 
@@ -70,16 +70,16 @@ def test_volume_vox_to_cart(volume, voxel, expected):
 
 
 def test_voxel_size(volume):
-    assert_allclose(volume.voxel_size, (2., 2., 2.))
+    assert_allclose(volume.voxel_size, (2.0, 2.0, 2.0))
 
 
 def test_to_vasp_volume(volume):
-    structure = Structure(lattice=volume.lattice,
-                          coords=[(0, 0, 0), (0.5, 0.5, 0.5)],
-                          species=['Si', 'Si'])
-    volume2 = Volume(data=volume.data,
-                     lattice=volume.lattice,
-                     label='free_energy')
+    structure = Structure(
+        lattice=volume.lattice,
+        coords=[(0, 0, 0), (0.5, 0.5, 0.5)],
+        species=['Si', 'Si'],
+    )
+    volume2 = Volume(data=volume.data, lattice=volume.lattice, label='free_energy')
     vol_data = volume.to_vasp_volume(structure=structure, other=[volume2])
     assert isinstance(vol_data, VolumetricData)
     assert set(vol_data.data.keys()) == {'total', 'free_energy'}