diff --git a/structuretoolkit/analyse/distance.py b/structuretoolkit/analyse/distance.py
index 9cee723da..c8ccfc03c 100644
--- a/structuretoolkit/analyse/distance.py
+++ b/structuretoolkit/analyse/distance.py
@@ -3,7 +3,13 @@
 import numpy as np
 
 
-def get_distances_array(structure: Atoms, p1: Optional[np.ndarray] = None, p2: Optional[np.ndarray] = None, mic: bool = True, vectors: bool = False):
+def get_distances_array(
+    structure: Atoms,
+    p1: Optional[np.ndarray] = None,
+    p2: Optional[np.ndarray] = None,
+    mic: bool = True,
+    vectors: bool = False,
+):
     """
     Return distance matrix of every position in p1 with every position in
     p2. If p2 is not set, it is assumed that distances between all
diff --git a/structuretoolkit/analyse/dscribe.py b/structuretoolkit/analyse/dscribe.py
index 5c23183fd..68b14875d 100644
--- a/structuretoolkit/analyse/dscribe.py
+++ b/structuretoolkit/analyse/dscribe.py
@@ -11,7 +11,7 @@ def soap_descriptor_per_atom(
     sigma: Optional[float] = 1.0,
     rbf: str = "gto",
     weighting: Optional[np.ndarray] = None,
-    average: str ="off",
+    average: str = "off",
     compression: dict = {"mode": "off", "species_weighting": None},
     species: Optional[list] = None,
     periodic: bool = True,
diff --git a/structuretoolkit/analyse/neighbors.py b/structuretoolkit/analyse/neighbors.py
index 1f4896ac8..ee7200b68 100644
--- a/structuretoolkit/analyse/neighbors.py
+++ b/structuretoolkit/analyse/neighbors.py
@@ -130,7 +130,12 @@ def copy(self):
         new_neigh._positions = self._positions.copy()
         return new_neigh
 
-    def _reshape(self, value: np.ndarray, key: Optional[str] = None, ref_vector: Optional[np.ndarray] = None):
+    def _reshape(
+        self,
+        value: np.ndarray,
+        key: Optional[str] = None,
+        ref_vector: Optional[np.ndarray] = None,
+    ):
         if value is None:
             raise ValueError("Neighbors not initialized yet")
         if key is None:
@@ -227,7 +232,9 @@ def _get_wrapped_indices(self) -> np.ndarray:
             return np.arange(len(self._ref_structure.positions))
         return self._wrapped_indices
 
-    def _get_wrapped_positions(self, positions: np.ndarray, distance_buffer: float = 1.0e-12):
+    def _get_wrapped_positions(
+        self, positions: np.ndarray, distance_buffer: float = 1.0e-12
+    ):
         if not self.wrap_positions:
             return np.asarray(positions)
         x = np.array(positions).copy()
@@ -319,7 +326,10 @@ def _get_vectors(
         return vectors
 
     def _estimate_num_neighbors(
-        self, num_neighbors: Optional[int] = None, cutoff_radius: float = np.inf, width_buffer: float = 1.2
+        self,
+        num_neighbors: Optional[int] = None,
+        cutoff_radius: float = np.inf,
+        width_buffer: float = 1.2,
     ):
         """
 
@@ -356,7 +366,10 @@ def _estimate_num_neighbors(
         return num_neighbors
 
     def _estimate_width(
-        self, num_neighbors: Optional[int] = None, cutoff_radius: float = np.inf, width_buffer: float = 1.2
+        self,
+        num_neighbors: Optional[int] = None,
+        cutoff_radius: float = np.inf,
+        width_buffer: float = 1.2,
     ):
         """
 
@@ -456,7 +469,13 @@ def _check_width(self, width: float, pbc: list[bool, bool, bool]):
                 return True
         return False
 
-    def get_spherical_harmonics(self, l: np.ndarray, m: np.ndarray, cutoff_radius: float = np.inf, rotation: Optional[np.ndarray] = None) -> np.ndarray:
+    def get_spherical_harmonics(
+        self,
+        l: np.ndarray,
+        m: np.ndarray,
+        cutoff_radius: float = np.inf,
+        rotation: Optional[np.ndarray] = None,
+    ) -> np.ndarray:
         """
         Args:
             l (int/numpy.array): Degree of the harmonic (int); must have ``l >= 0``.
@@ -497,7 +516,9 @@ def get_spherical_harmonics(self, l: np.ndarray, m: np.ndarray, cutoff_radius: f
             within_cutoff, axis=-1
         )
 
-    def get_steinhardt_parameter(self, l: np.ndarray, cutoff_radius: float = np.inf) -> np.ndarray:
+    def get_steinhardt_parameter(
+        self, l: np.ndarray, cutoff_radius: float = np.inf
+    ) -> np.ndarray:
         """
         Args:
             l (int/numpy.array): Order of Steinhardt parameter
@@ -790,7 +811,10 @@ def get_global_shells(
         return self._reshape(shells, key=mode)
 
     def get_shell_matrix(
-        self, chemical_pair: Optional[list] = None, cluster_by_distances: bool = False, cluster_by_vecs: bool = False
+        self,
+        chemical_pair: Optional[list] = None,
+        cluster_by_distances: bool = False,
+        cluster_by_vecs: bool = False,
     ):
         """
         Shell matrices for pairwise interaction. Note: The matrices are always symmetric, meaning if you
@@ -851,7 +875,9 @@ def get_shell_matrix(
                 )
         return shell_matrix
 
-    def find_neighbors_by_vector(self, vector: np.ndarray, return_deviation: bool = False) -> np.ndarray:
+    def find_neighbors_by_vector(
+        self, vector: np.ndarray, return_deviation: bool = False
+    ) -> np.ndarray:
         """
         Args:
             vector (list/np.ndarray): vector by which positions are translated (and neighbors are searched)
@@ -1052,7 +1078,12 @@ def __probe_cluster(self, c_count: int, neighbors: list, id_list: list):
                     self.__probe_cluster(c_count, nbrs, id_list)
 
     # TODO: combine with corresponding routine in plot3d
-    def get_bonds(self, radius: float = np.inf, max_shells: Optional[int] = None, prec: float = 0.1):
+    def get_bonds(
+        self,
+        radius: float = np.inf,
+        max_shells: Optional[int] = None,
+        prec: float = 0.1,
+    ):
         """
 
         Args:
@@ -1110,7 +1141,7 @@ def get_volume_of_n_sphere_in_p_norm(n: int = 3, p: int = 2) -> float:
 
 def get_neighbors(
     structure: Atoms,
-    num_neighbors: int =12,
+    num_neighbors: int = 12,
     tolerance: int = 2,
     id_list: Optional[list] = None,
     cutoff_radius: float = np.inf,
diff --git a/structuretoolkit/analyse/phonopy.py b/structuretoolkit/analyse/phonopy.py
index c5683bddc..bc587fce5 100644
--- a/structuretoolkit/analyse/phonopy.py
+++ b/structuretoolkit/analyse/phonopy.py
@@ -17,7 +17,9 @@
 __date__ = "Sep 1, 2018"
 
 
-def get_equivalent_atoms(structure: Atoms, symprec: float = 1e-5, angle_tolerance: float = -1.0):
+def get_equivalent_atoms(
+    structure: Atoms, symprec: float = 1e-5, angle_tolerance: float = -1.0
+):
     """
     Args: (read phonopy.structure.spglib for more details)
         symprec:
diff --git a/structuretoolkit/analyse/pyscal.py b/structuretoolkit/analyse/pyscal.py
index 0d30e6ae2..e6119807c 100644
--- a/structuretoolkit/analyse/pyscal.py
+++ b/structuretoolkit/analyse/pyscal.py
@@ -64,7 +64,9 @@ def get_steinhardt_parameters(
         return sysq
 
 
-def get_centro_symmetry_descriptors(structure: Atoms, num_neighbors: int = 12) -> np.ndarray:
+def get_centro_symmetry_descriptors(
+    structure: Atoms, num_neighbors: int = 12
+) -> np.ndarray:
     """
     Analyse centrosymmetry parameter
 
@@ -183,7 +185,9 @@ def get_diamond_structure_descriptors(
         )
 
 
-def get_adaptive_cna_descriptors(structure: Atoms, mode: str = "total", ovito_compatibility: bool = False) -> np.ndarray:
+def get_adaptive_cna_descriptors(
+    structure: Atoms, mode: str = "total", ovito_compatibility: bool = False
+) -> np.ndarray:
     """
     Use common neighbor analysis
 
@@ -253,7 +257,7 @@ def get_voronoi_volumes(structure: Atoms) -> np.ndarray:
 
 def find_solids(
     structure: Atoms,
-    neighbor_method: str ="cutoff",
+    neighbor_method: str = "cutoff",
     cutoff: float = 0.0,
     bonds: float = 0.5,
     threshold: float = 0.5,
diff --git a/structuretoolkit/analyse/spatial.py b/structuretoolkit/analyse/spatial.py
index 2db4f54e0..c7c1ce8f0 100644
--- a/structuretoolkit/analyse/spatial.py
+++ b/structuretoolkit/analyse/spatial.py
@@ -28,7 +28,9 @@
 __date__ = "Sep 1, 2017"
 
 
-def get_mean_positions(positions: np.ndarray, cell: np.ndarray, pbc: np.ndarray, labels: np.ndarray) -> np.ndarray:
+def get_mean_positions(
+    positions: np.ndarray, cell: np.ndarray, pbc: np.ndarray, labels: np.ndarray
+) -> np.ndarray:
     """
     This function calculates the average position(-s) across periodic boundary conditions according
     to the labels
@@ -57,7 +59,9 @@ def get_mean_positions(positions: np.ndarray, cell: np.ndarray, pbc: np.ndarray,
     return mean_positions
 
 
-def create_gridpoints(structure: Atoms, n_gridpoints_per_angstrom: int = 5) -> np.ndarray:
+def create_gridpoints(
+    structure: Atoms, n_gridpoints_per_angstrom: int = 5
+) -> np.ndarray:
     cell = get_vertical_length(structure=structure)
     n_points = (n_gridpoints_per_angstrom * cell).astype(int)
     positions = np.meshgrid(
@@ -67,12 +71,16 @@ def create_gridpoints(structure: Atoms, n_gridpoints_per_angstrom: int = 5) -> n
     return np.einsum("ji,nj->ni", structure.cell, positions)
 
 
-def remove_too_close(positions: np.ndarray, structure: Atoms, min_distance: float = 1) -> np.ndarray:
+def remove_too_close(
+    positions: np.ndarray, structure: Atoms, min_distance: float = 1
+) -> np.ndarray:
     neigh = get_neighborhood(structure=structure, positions=positions, num_neighbors=1)
     return positions[neigh.distances.flatten() > min_distance]
 
 
-def set_to_high_symmetry_points(positions: np.ndarray, structure: Atoms, neigh, decimals: int = 4) -> np.ndarray:
+def set_to_high_symmetry_points(
+    positions: np.ndarray, structure: Atoms, neigh, decimals: int = 4
+) -> np.ndarray:
     for _ in range(10):
         neigh = neigh.get_neighborhood(positions)
         dx = np.mean(neigh.vecs, axis=-2)
@@ -87,7 +95,14 @@ def set_to_high_symmetry_points(positions: np.ndarray, structure: Atoms, neigh,
     raise ValueError("High symmetry points could not be detected")
 
 
-def cluster_by_steinhardt(positions: np.ndarray, neigh, l_values: List[int], q_eps: float, var_ratio: float, min_samples: int) -> np.ndarray:
+def cluster_by_steinhardt(
+    positions: np.ndarray,
+    neigh,
+    l_values: List[int],
+    q_eps: float,
+    var_ratio: float,
+    min_samples: int,
+) -> np.ndarray:
     """
     Clusters candidate positions via Steinhardt parameters and the variance in distances to host atoms.
 
@@ -240,7 +255,9 @@ def __init__(
         self._positions = None
         self.structure = structure
 
-    def run_workflow(self, positions: Optional[np.ndarray] = None, steps: int = -1) -> np.ndarray:
+    def run_workflow(
+        self, positions: Optional[np.ndarray] = None, steps: int = -1
+    ) -> np.ndarray:
         if positions is None:
             positions = self.initial_positions.copy()
         for ii, ww in enumerate(self.workflow):
@@ -470,7 +487,10 @@ def get_layers(
 
 
 def get_voronoi_vertices(
-    structure: Atoms, epsilon: float = 2.5e-4, distance_threshold: float = 0, width_buffer: float = 10.0
+    structure: Atoms,
+    epsilon: float = 2.5e-4,
+    distance_threshold: float = 0,
+    width_buffer: float = 10.0,
 ) -> np.ndarray:
     """
     Get voronoi vertices of the box.
@@ -578,7 +598,11 @@ def get_delaunay_neighbors(structure: Atoms, width_buffer: float = 10.0) -> np.n
 
 
 def get_cluster_positions(
-    structure: Atoms, positions: Optional[np.ndarray] = None, eps: float = 1.0, buffer_width: Optional[float] =None, return_labels: bool = False
+    structure: Atoms,
+    positions: Optional[np.ndarray] = None,
+    eps: float = 1.0,
+    buffer_width: Optional[float] = None,
+    return_labels: bool = False,
 ) -> np.ndarray:
     """
     Cluster positions according to the distances. Clustering algorithm uses DBSCAN:
diff --git a/structuretoolkit/analyse/strain.py b/structuretoolkit/analyse/strain.py
index 55ac4ffd1..a65dcd32a 100644
--- a/structuretoolkit/analyse/strain.py
+++ b/structuretoolkit/analyse/strain.py
@@ -26,7 +26,11 @@ class Strain:
     """
 
     def __init__(
-        self, structure: Atoms, ref_structure: Atoms, num_neighbors: Optional[int] = None, only_bulk_type: bool = False
+        self,
+        structure: Atoms,
+        ref_structure: Atoms,
+        num_neighbors: Optional[int] = None,
+        only_bulk_type: bool = False,
     ):
         """
 
@@ -71,7 +75,9 @@ def _nullify_non_bulk(self) -> np.ndarray:
             self.structure.analyse.pyscal_cna_adaptive(mode="str") != self.crystal_phase
         )
 
-    def _get_perpendicular_unit_vectors(self, vec: np.ndarray, vec_axis: Optional[np.ndarray] = None) -> np.ndarray:
+    def _get_perpendicular_unit_vectors(
+        self, vec: np.ndarray, vec_axis: Optional[np.ndarray] = None
+    ) -> np.ndarray:
         if vec_axis is not None:
             vec_axis = self._get_safe_unit_vectors(vec_axis)
             vec = np.array(
@@ -80,7 +86,9 @@ def _get_perpendicular_unit_vectors(self, vec: np.ndarray, vec_axis: Optional[np
         return self._get_safe_unit_vectors(vec)
 
     @staticmethod
-    def _get_safe_unit_vectors(vectors: np.ndarray, minimum_value: float = 1.0e-8) -> np.ndarray:
+    def _get_safe_unit_vectors(
+        vectors: np.ndarray, minimum_value: float = 1.0e-8
+    ) -> np.ndarray:
         v = np.linalg.norm(vectors, axis=-1)
         v += (v < minimum_value) * minimum_value
         return np.einsum("...i,...->...i", vectors, 1 / v)
@@ -94,7 +102,12 @@ def _get_angle(self, v: np.ndarray, w: np.ndarray) -> np.ndarray:
             prod[np.absolute(prod) > 1] = np.sign(prod)[np.absolute(prod) > 1]
         return np.arccos(prod)
 
-    def _get_rotation_from_vectors(self, vec_before: np.ndarray, vec_after: np.ndarray, vec_axis: Optional[np.ndarray] = None) -> np.ndarray:
+    def _get_rotation_from_vectors(
+        self,
+        vec_before: np.ndarray,
+        vec_after: np.ndarray,
+        vec_axis: Optional[np.ndarray] = None,
+    ) -> np.ndarray:
         v = self._get_perpendicular_unit_vectors(vec_before, vec_axis)
         w = self._get_perpendicular_unit_vectors(vec_after, vec_axis)
         if vec_axis is None:
@@ -132,7 +145,9 @@ def rotations(self) -> np.ndarray:
         return self._rotations
 
     @staticmethod
-    def _get_best_match_indices(coords: np.ndarray, ref_coord: np.ndarray) -> np.ndarray:
+    def _get_best_match_indices(
+        coords: np.ndarray, ref_coord: np.ndarray
+    ) -> np.ndarray:
         distances = np.linalg.norm(
             coords[:, :, None, :] - ref_coord[None, None, :, :], axis=-1
         )
@@ -195,7 +210,7 @@ def strain(self) -> np.ndarray:
 def get_strain(
     structure: Atoms,
     ref_structure: Atoms,
-    num_neighbors: Optional[int] =None,
+    num_neighbors: Optional[int] = None,
     only_bulk_type: bool = False,
     return_object: bool = False,
 ):
diff --git a/structuretoolkit/analyse/symmetry.py b/structuretoolkit/analyse/symmetry.py
index e3f100362..19e699f4b 100644
--- a/structuretoolkit/analyse/symmetry.py
+++ b/structuretoolkit/analyse/symmetry.py
@@ -230,7 +230,9 @@ def symmetrize_vectors(
             np.einsum("ijk->jki", v_reshaped)[self.permutations],
         ).reshape(np.shape(vectors)) / len(self["rotations"])
 
-    def _get_spglib_cell(self, use_elements: Optional[bool] = None, use_magmoms: Optional[bool] = None) -> tuple:
+    def _get_spglib_cell(
+        self, use_elements: Optional[bool] = None, use_magmoms: Optional[bool] = None
+    ) -> tuple:
         lattice = np.array(self._structure.get_cell(), dtype="double", order="C")
         positions = np.array(
             self._structure.get_scaled_positions(wrap=False), dtype="double", order="C"
@@ -326,7 +328,10 @@ def spacegroup(self) -> dict:
         }
 
     def get_primitive_cell(
-        self, standardize: bool = False, use_elements: Optional[bool] = None, use_magmoms: Optional[bool] = None
+        self,
+        standardize: bool = False,
+        use_elements: Optional[bool] = None,
+        use_magmoms: Optional[bool] = None,
     ) -> Atoms:
         """
         Get primitive cell of a given structure.
diff --git a/structuretoolkit/common/helper.py b/structuretoolkit/common/helper.py
index fab2dd730..99a3246a1 100644
--- a/structuretoolkit/common/helper.py
+++ b/structuretoolkit/common/helper.py
@@ -6,7 +6,11 @@
 
 
 def get_extended_positions(
-    structure: Atoms, width: float, return_indices: bool =False, norm_order: int = 2, positions: Optional[np.ndarray] = None
+    structure: Atoms,
+    width: float,
+    return_indices: bool = False,
+    norm_order: int = 2,
+    positions: Optional[np.ndarray] = None,
 ):
     """
     Get all atoms in the boundary around the supercell which have a distance
@@ -79,7 +83,9 @@ def get_vertical_length(structure: Atoms, norm_order: int = 2):
     )
 
 
-def get_wrapped_coordinates(structure: Atoms, positions: np.ndarray, epsilon: float=1.0e-8) -> np.ndarray:
+def get_wrapped_coordinates(
+    structure: Atoms, positions: np.ndarray, epsilon: float = 1.0e-8
+) -> np.ndarray:
     """
     Return coordinates in wrapped in the periodic cell
 
@@ -153,7 +159,9 @@ def get_average_of_unique_labels(labels: np.ndarray, values: np.ndarray) -> floa
     return mean_values
 
 
-def center_coordinates_in_unit_cell(structure: Atoms, origin: float = 0.0, eps: float = 1e-4) -> Atoms:
+def center_coordinates_in_unit_cell(
+    structure: Atoms, origin: float = 0.0, eps: float = 1e-4
+) -> Atoms:
     """
     Wrap atomic coordinates within the supercell.
 
@@ -174,7 +182,9 @@ def center_coordinates_in_unit_cell(structure: Atoms, origin: float = 0.0, eps:
     return structure
 
 
-def apply_strain(structure: Atoms, epsilon: float, return_box: bool = False, mode: str = "linear"):
+def apply_strain(
+    structure: Atoms, epsilon: float, return_box: bool = False, mode: str = "linear"
+):
     """
     Apply a given strain on the structure. It applies the matrix `F` in the manner:
 
diff --git a/structuretoolkit/common/phonopy.py b/structuretoolkit/common/phonopy.py
index c37cae9e7..c6915950d 100644
--- a/structuretoolkit/common/phonopy.py
+++ b/structuretoolkit/common/phonopy.py
@@ -1,5 +1,6 @@
 from ase.atoms import Atoms
 
+
 def phonopy_to_atoms(ph_atoms) -> Atoms:
     """
     Convert Phonopy Atoms to ASE-like Atoms
diff --git a/structuretoolkit/visualize.py b/structuretoolkit/visualize.py
index 420e49bfb..4f0f67078 100644
--- a/structuretoolkit/visualize.py
+++ b/structuretoolkit/visualize.py
@@ -32,7 +32,7 @@ def plot3d(
     select_atoms: Optional[np.ndarray] = None,
     background: str = "white",
     color_scheme: Optional[str] = None,
-    colors: Optional[np.ndarray] =None,
+    colors: Optional[np.ndarray] = None,
     scalar_field: Optional[np.ndarray] = None,
     scalar_start: Optional[float] = None,
     scalar_end: Optional[float] = None,
@@ -512,7 +512,14 @@ def _plot3d_ase(
     return view
 
 
-def _ngl_write_cell(a1: float, a2: float, a3: float, f1: float = 90.0, f2: float = 90.0, f3: float = 90.0):
+def _ngl_write_cell(
+    a1: float,
+    a2: float,
+    a3: float,
+    f1: float = 90.0,
+    f2: float = 90.0,
+    f3: float = 90.0,
+):
     """
     Writes a PDB-formatted line to represent the simulation cell.
 
@@ -568,7 +575,9 @@ def _ngl_write_atom(
     )
 
 
-def _ngl_write_structure(elements: np.ndarray, positions: np.ndarray, cell: np.ndarray) -> str:
+def _ngl_write_structure(
+    elements: np.ndarray, positions: np.ndarray, cell: np.ndarray
+) -> str:
     """
     Turns structure information into a NGLView-readable protein-database-formatted string.
 
@@ -604,7 +613,9 @@ def _ngl_write_structure(elements: np.ndarray, positions: np.ndarray, cell: np.n
     return pdb_str
 
 
-def _atomic_number_to_radius(atomic_number: int, shift: float =0.2, slope: float = 0.1, scale: float = 1.0) -> float:
+def _atomic_number_to_radius(
+    atomic_number: int, shift: float = 0.2, slope: float = 0.1, scale: float = 1.0
+) -> float:
     """
     Give the atomic radius for plotting, which scales like the root of the atomic number.
 
@@ -621,7 +632,11 @@ def _atomic_number_to_radius(atomic_number: int, shift: float =0.2, slope: float
 
 
 def _add_colorscheme_spacefill(
-    view, elements: np.ndarray, atomic_numbers: np.ndarray, particle_size: float, scheme: str = "element"
+    view,
+    elements: np.ndarray,
+    atomic_numbers: np.ndarray,
+    particle_size: float,
+    scheme: str = "element",
 ):
     """
     Set NGLView spacefill parameters according to a color-scheme.
@@ -651,7 +666,9 @@ def _add_colorscheme_spacefill(
     return view
 
 
-def _add_custom_color_spacefill(view, atomic_numbers: np.ndarray, particle_size: float, colors: np.ndarray):
+def _add_custom_color_spacefill(
+    view, atomic_numbers: np.ndarray, particle_size: float, colors: np.ndarray
+):
     """
     Set NGLView spacefill parameters according to per-atom colors.
 
@@ -674,7 +691,12 @@ def _add_custom_color_spacefill(view, atomic_numbers: np.ndarray, particle_size:
     return view
 
 
-def _scalars_to_hex_colors(scalar_field: np.ndarray, start: Optional[float] = None, end: Optional[float] = None, cmap=None):
+def _scalars_to_hex_colors(
+    scalar_field: np.ndarray,
+    start: Optional[float] = None,
+    end: Optional[float] = None,
+    cmap=None,
+):
     """
     Convert scalar values to hex codes using a colormap.
 
@@ -753,7 +775,9 @@ def _get_orientation(view_plane: np.ndarray) -> np.ndarray:
     ).T
 
 
-def _get_flattened_orientation(view_plane: np.ndarray, distance_from_camera: float) -> list:
+def _get_flattened_orientation(
+    view_plane: np.ndarray, distance_from_camera: float
+) -> list:
     """
     A helper method to plot3d, which generates a rotation matrix from the input `view_plane`, and returns a
     flattened list of len = 16. This flattened list becomes the input argument to `view.contol.orient`.