diff --git a/README.md b/README.md
index e5b7c3db0..b49035037 100644
--- a/README.md
+++ b/README.md
@@ -10,46 +10,54 @@ represent atomistic structures in python. The `structuretoolkit` is integrated i
 ## Disclaimer 
 The `structuretoolkit` is currently under development. 
 
-## Example 
+## Example
+
 ```python
 import structuretoolkit as stk
 from ase.build import bulk
 
 structure = bulk("Al", cubic=True)
-stk.analyse_cna_adaptive(structure)
+stk.analyse.get_adaptive_cna_descriptors(structure)
 stk.plot3d(structure)
 ```
 
 ## Features 
 ### Analysis
-* `get_neighbors`
-* `get_neighborhood`
-* `analyse_phonopy_equivalent_atoms`
-* `get_steinhardt_parameter_structure`
-* `analyse_centro_symmetry` 
-* `analyse_diamond_structure` 
-* `analyse_cna_adaptive` 
-* `analyse_voronoi_volume` 
-* `analyse_find_solids`
-* `get_mean_positions`
-* `get_average_of_unique_labels`
-* `get_interstitials`
-* `get_layers`
-* `get_voronoi_vertices`
-* `get_voronoi_neighbors`
-* `get_delaunay_neighbors`
-* `cluster_positions`
-* `get_strain`
+* `stk.analyse.get_neighbors()`
+* `stk.analyse.get_neighborhood()`
+* `stk.analyse.get_equivalent_atoms()`
+* `stk.analyse.get_steinhardt_parameters()`
+* `stk.analyse.get_centro_symmetry_descriptors()` 
+* `stk.analyse.get_diamond_structure_descriptors()` 
+* `stk.analyse.get_adaptive_cna_descriptors()` 
+* `stk.analyse.get_voronoi_volumes()` 
+* `stk.analyse.find_solids()`
+* `stk.analyse.get_mean_positions()`
+* `stk.analyse.get_average_of_unique_labels()`
+* `stk.analyse.get_interstitials()`
+* `stk.analyse.get_layers()`
+* `stk.analyse.get_voronoi_vertices()`
+* `stk.analyse.get_voronoi_neighbors()`
+* `stk.analyse.get_delaunay_neighbors()`
+* `stk.analyse.get_cluster_positions()`
+* `stk.analyse.get_strain()`
 
 ### Build
-* `grainboundary_build`
-* `grainboundary_info`
-* `get_sqs_structures`
-* `B2`
-* `C14`
-* `C15`
-* `C36`
-* `D03`
+* `stk.build.get_grainboundary_info()`
+* `stk.build.grainboundary()`
+* `stk.build.high_index_surface()`
+* `stk.build.get_high_index_surface_info()`
+* `stk.build.sqs_structures()`
+* `stk.build.B2()`
+* `stk.build.C14()`
+* `stk.build.C15()`
+* `stk.build.C36()`
+* `stk.build.D03()`
 
 ### Visualize 
-* `plot3d`
+* `stk.visualize.plot3d()`
+
+### Common 
+* `stk.common.ase_to_pymatgen()`
+* `stk.common.pymatgen_to_ase()`
+* `stk.common.ase_to_pyscal()`
diff --git a/setup.py b/setup.py
index aa2a38fe2..962389d01 100644
--- a/setup.py
+++ b/setup.py
@@ -30,15 +30,20 @@
     keywords='pyiron',
     packages=find_packages(exclude=["*tests*", "*docs*", "*binder*", "*conda*", "*notebooks*", "*.ci_support*"]),
     install_requires=[
-        'aimsgb==0.1.1',
         'ase==3.22.1',
-        'matplotlib==3.7.0',
-        'numpy==1.24.2',
-        'phonopy==2.17.1',
-        'pymatgen==2022.11.7',
-        'scipy==1.10.0',
-        'scikit-learn==1.2.1',
-        'spglib==2.0.2',
+        'matplotlib==3.7.0',  # ase already requires matplotlib
+        'numpy==1.24.2',  # ase already requires numpy
+        'scipy==1.10.0',  # ase already requires scipy
     ],
+    extras_require={
+        "grainboundary": ['aimsgb==0.1.1', 'pymatgen==2022.11.7'],
+        "pyscal": ['pyscal2==2.10.20'],
+        "nglview": ['nglview==3.0.4'],
+        "plotly": ['plotly==5.14.1'],
+        "clusters": ['scikit-learn==1.2.1'],
+        "symmetry": ['spglib==2.0.2'],
+        "surface": ['spglib==2.0.2', 'pymatgen==2022.11.7'],
+        "phonopy": ['phonopy==2.17.1', 'spglib==2.0.2'],
+    },
     cmdclass=versioneer.get_cmdclass(),
 )
diff --git a/structuretoolkit/__init__.py b/structuretoolkit/__init__.py
index 0a114bac8..666658d3d 100644
--- a/structuretoolkit/__init__.py
+++ b/structuretoolkit/__init__.py
@@ -1,16 +1,16 @@
 # Analyse
-from structuretoolkit.analyse.distance import get_distances_array, find_mic
-from structuretoolkit.analyse.neighbors import get_neighbors, get_neighborhood
-from structuretoolkit.analyse.phonopy import analyse_phonopy_equivalent_atoms
-from structuretoolkit.analyse.pyscal import (
-    get_steinhardt_parameter_structure,
-    analyse_centro_symmetry,
-    analyse_diamond_structure,
-    analyse_cna_adaptive,
-    analyse_voronoi_volume,
-    analyse_find_solids,
-)
-from structuretoolkit.analyse.spatial import (
+from structuretoolkit.analyse import (
+    get_distances_array,
+    find_mic,
+    get_neighbors,
+    get_neighborhood,
+    get_equivalent_atoms,
+    get_steinhardt_parameters,
+    get_centro_symmetry_descriptors,
+    get_diamond_structure_descriptors,
+    get_adaptive_cna_descriptors,
+    get_voronoi_volumes,
+    find_solids,
     get_mean_positions,
     get_average_of_unique_labels,
     get_interstitials,
@@ -18,19 +18,47 @@
     get_voronoi_vertices,
     get_voronoi_neighbors,
     get_delaunay_neighbors,
-    cluster_positions,
+    get_cluster_positions,
+    get_strain,
+    get_symmetry,
+    # for backwards compatibility
+    get_cluster_positions as cluster_positions,
+    get_equivalent_atoms as analyse_phonopy_equivalent_atoms,
+    get_steinhardt_parameters as get_steinhardt_parameter_structure,
+    get_centro_symmetry_descriptors as analyse_centro_symmetry,
+    get_diamond_structure_descriptors as analyse_diamond_structure,
+    get_adaptive_cna_descriptors as analyse_cna_adaptive,
+    get_voronoi_volumes as analyse_voronoi_volume,
+    find_solids as analyse_find_solids,
 )
-from structuretoolkit.analyse.strain import get_strain
-from structuretoolkit.analyse.symmetry import get_symmetry
 
 # Build
-from structuretoolkit.build.aimsgb import grainboundary_build, grainboundary_info
-from structuretoolkit.build.compound import B2, C14, C15, C36, D03
-from structuretoolkit.build.sqs import get_sqs_structures
-from structuretoolkit.build.surface import high_index_surface, high_index_surface_info
+from structuretoolkit.build import (
+    grainboundary,
+    get_grainboundary_info,
+    B2,
+    C14,
+    C15,
+    C36,
+    D03,
+    sqs_structures,
+    high_index_surface,
+    get_high_index_surface_info,
+    # for backwards compatibility
+    grainboundary as grainboundary_build,
+    get_grainboundary_info as grainboundary_info,
+    sqs_structures as get_sqs_structures,
+    get_high_index_surface_info as high_index_surface_info,
+)
 
-# Other
-from structuretoolkit.helper import (
+# Visualize
+from structuretoolkit.visualize import plot3d
+
+# Common
+from structuretoolkit.common import (
+    ase_to_pymatgen,
+    pymatgen_to_ase,
+    ase_to_pyscal,
     get_atomic_numbers,
     get_extended_positions,
     get_vertical_length,
@@ -38,5 +66,5 @@
     select_index,
     center_coordinates_in_unit_cell,
     apply_strain,
+    SymmetryError,
 )
-from structuretoolkit.visualize import plot3d
diff --git a/structuretoolkit/analyse/__init__.py b/structuretoolkit/analyse/__init__.py
index e69de29bb..9fa35065c 100644
--- a/structuretoolkit/analyse/__init__.py
+++ b/structuretoolkit/analyse/__init__.py
@@ -0,0 +1,48 @@
+from structuretoolkit.analyse.distance import get_distances_array, find_mic
+from structuretoolkit.analyse.neighbors import get_neighbors, get_neighborhood
+from structuretoolkit.analyse.phonopy import get_equivalent_atoms
+from structuretoolkit.analyse.pyscal import (
+    get_steinhardt_parameters,
+    get_centro_symmetry_descriptors,
+    get_diamond_structure_descriptors,
+    get_adaptive_cna_descriptors,
+    get_voronoi_volumes,
+    find_solids,
+    ase_to_pyscal,
+)
+from structuretoolkit.analyse.spatial import (
+    get_mean_positions,
+    get_average_of_unique_labels,
+    get_interstitials,
+    get_layers,
+    get_voronoi_vertices,
+    get_voronoi_neighbors,
+    get_delaunay_neighbors,
+    get_cluster_positions,
+)
+from structuretoolkit.analyse.strain import get_strain
+
+
+def get_symmetry(
+    structure, use_magmoms=False, use_elements=True, symprec=1e-5, angle_tolerance=-1.0
+):
+    """
+    Args:
+        structure (ase.atoms.Atoms): Atomistic Structure object
+        use_magmoms (bool): Whether to consider magnetic moments (cf. get_initial_magnetic_moments())
+        use_elements (bool): If False, chemical elements will be ignored
+        symprec (float): Symmetry search precision
+        angle_tolerance (float): Angle search tolerance
+
+    Returns:
+        symmetry (:class:`structuretoolkit.analyse.symmetry.Symmetry`): Symmetry class
+    """
+    from structuretoolkit.analyse.symmetry import get_symmetry
+
+    return get_symmetry(
+        structure=structure,
+        use_magmoms=use_magmoms,
+        use_elements=use_elements,
+        symprec=symprec,
+        angle_tolerance=angle_tolerance,
+    )
diff --git a/structuretoolkit/analyse/neighbors.py b/structuretoolkit/analyse/neighbors.py
index 15eb906c1..1fa9daaeb 100644
--- a/structuretoolkit/analyse/neighbors.py
+++ b/structuretoolkit/analyse/neighbors.py
@@ -3,7 +3,6 @@
 # Distributed under the terms of "New BSD License", see the LICENSE file.
 
 import numpy as np
-from sklearn.cluster import AgglomerativeClustering
 from scipy.sparse import coo_matrix
 from scipy.special import gamma
 from scipy.spatial.transform import Rotation
@@ -11,7 +10,10 @@
 from scipy.spatial import cKDTree
 import warnings
 import itertools
-from structuretoolkit.helper import get_extended_positions, get_average_of_unique_labels
+from structuretoolkit.common.helper import (
+    get_extended_positions,
+    get_average_of_unique_labels,
+)
 
 __author__ = "Joerg Neugebauer, Sam Waseda"
 __copyright__ = (
@@ -905,6 +907,8 @@ def cluster_by_vecs(
                 `euclidean` is accepted.
 
         """
+        from sklearn.cluster import AgglomerativeClustering
+
         if distance_threshold is None and n_clusters is None:
             distance_threshold = np.min(self.filled.distances)
         dr = self.flattened.vecs
@@ -952,6 +956,8 @@ def cluster_by_distances(
                 obtained from the clustered vectors is used for the distance clustering. Otherwise
                 neigh.distances is used.
         """
+        from sklearn.cluster import AgglomerativeClustering
+
         if distance_threshold is None:
             distance_threshold = 0.1 * np.min(self.flattened.distances)
         dr = self.flattened.distances
diff --git a/structuretoolkit/analyse/phonopy.py b/structuretoolkit/analyse/phonopy.py
index 2d0247dc5..3b9c8d895 100644
--- a/structuretoolkit/analyse/phonopy.py
+++ b/structuretoolkit/analyse/phonopy.py
@@ -3,8 +3,6 @@
 # Distributed under the terms of "New BSD License", see the LICENSE file.
 
 import numpy as np
-from phonopy.structure.atoms import PhonopyAtoms
-import spglib as spg
 
 __author__ = "Osamu Waseda"
 __copyright__ = (
@@ -18,7 +16,7 @@
 __date__ = "Sep 1, 2018"
 
 
-def analyse_phonopy_equivalent_atoms(structure, symprec=1e-5, angle_tolerance=-1.0):
+def get_equivalent_atoms(structure, symprec=1e-5, angle_tolerance=-1.0):
     """
     Args: (read phonopy.structure.spglib for more details)
         symprec:
@@ -29,6 +27,9 @@ def analyse_phonopy_equivalent_atoms(structure, symprec=1e-5, angle_tolerance=-1
                 is used to judge symmetry.
 
     """
+    import spglib as spg
+    from phonopy.structure.atoms import PhonopyAtoms
+
     positions = structure.get_scaled_positions()
     cell = structure.cell
     types = structure.get_chemical_symbols()
diff --git a/structuretoolkit/analyse/pyscal.py b/structuretoolkit/analyse/pyscal.py
index 323ef8e7a..cde8ea2d1 100644
--- a/structuretoolkit/analyse/pyscal.py
+++ b/structuretoolkit/analyse/pyscal.py
@@ -4,8 +4,7 @@
 
 import numpy as np
 from ase.atoms import Atoms
-import pyscal.core as pc
-from sklearn import cluster
+from structuretoolkit.common.pyscal import ase_to_pyscal
 
 __author__ = "Sarath Menon, Jan Janssen"
 __copyright__ = (
@@ -19,7 +18,7 @@
 __date__ = "Nov 6, 2019"
 
 
-def get_steinhardt_parameter_structure(
+def get_steinhardt_parameters(
     structure,
     neighbor_method="cutoff",
     cutoff=0,
@@ -43,7 +42,7 @@ def get_steinhardt_parameter_structure(
         numpy.ndarray: (number of q's, number of atoms) shaped array of q parameters
         numpy.ndarray: If `clustering=True`, an additional per-atom array of cluster ids is also returned
     """
-    sys = ase_to_pyscal_system(structure)
+    sys = ase_to_pyscal(structure)
     q = (4, 6) if q is None else q
 
     sys.find_neighbors(method=neighbor_method, cutoff=cutoff)
@@ -53,6 +52,8 @@ def get_steinhardt_parameter_structure(
     sysq = np.array(sys.get_qvals(q, averaged=averaged))
 
     if n_clusters is not None:
+        from sklearn import cluster
+
         cl = cluster.KMeans(n_clusters=n_clusters)
 
         ind = cl.fit(list(zip(*sysq))).labels_
@@ -61,7 +62,7 @@ def get_steinhardt_parameter_structure(
         return sysq
 
 
-def analyse_centro_symmetry(structure, num_neighbors=12):
+def get_centro_symmetry_descriptors(structure, num_neighbors=12):
     """
     Analyse centrosymmetry parameter
 
@@ -72,11 +73,13 @@ def analyse_centro_symmetry(structure, num_neighbors=12):
     Returns:
         csm (list) : list of centrosymmetry parameter
     """
-    sys = ase_to_pyscal_system(structure)
+    sys = ase_to_pyscal(structure)
     return np.array(sys.calculate_centrosymmetry(nmax=num_neighbors))
 
 
-def analyse_diamond_structure(structure, mode="total", ovito_compatibility=False):
+def get_diamond_structure_descriptors(
+    structure, mode="total", ovito_compatibility=False
+):
     """
     Analyse diamond structure
 
@@ -93,7 +96,7 @@ def analyse_diamond_structure(structure, mode="total", ovito_compatibility=False
     Returns:
         (depends on `mode`)
     """
-    sys = ase_to_pyscal_system(structure)
+    sys = ase_to_pyscal(structure)
     diamond_dict = sys.identify_diamond()
 
     ovito_identifiers = [
@@ -178,7 +181,7 @@ def analyse_diamond_structure(structure, mode="total", ovito_compatibility=False
         )
 
 
-def analyse_cna_adaptive(structure, mode="total", ovito_compatibility=False):
+def get_adaptive_cna_descriptors(structure, mode="total", ovito_compatibility=False):
     """
     Use common neighbor analysis
 
@@ -195,7 +198,7 @@ def analyse_cna_adaptive(structure, mode="total", ovito_compatibility=False):
     Returns:
         (depends on `mode`)
     """
-    sys = ase_to_pyscal_system(structure)
+    sys = ase_to_pyscal(structure)
     if mode not in ["total", "numeric", "str"]:
         raise ValueError("Unsupported mode")
 
@@ -233,39 +236,20 @@ def analyse_cna_adaptive(structure, mode="total", ovito_compatibility=False):
             )
 
 
-def analyse_voronoi_volume(structure):
+def get_voronoi_volumes(structure):
     """
     Calculate the Voronoi volume of atoms
 
     Args:
         structure : (ase.atoms.Atoms): The structure to analyze.
     """
-    sys = ase_to_pyscal_system(structure)
+    sys = ase_to_pyscal(structure)
     sys.find_neighbors(method="voronoi")
     structure = sys.atoms
     return np.array([atom.volume for atom in structure])
 
 
-def ase_to_pyscal_system(structure):
-    """
-    Converts atoms to ase atoms and than to a pyscal system.
-    Also adds the pyscal publication.
-
-    Args:
-        structure (ase.atoms.Atoms): Structure to convert.
-
-    Returns:
-        Pyscal system: See the pyscal documentation.
-    """
-    sys = pc.System()
-    sys.read_inputfile(
-        structure,
-        format="ase",
-    )
-    return sys
-
-
-def analyse_find_solids(
+def find_solids(
     structure,
     neighbor_method="cutoff",
     cutoff=0,
@@ -296,7 +280,7 @@ def analyse_find_solids(
         int: number of solids,
         pyscal system: pyscal system when return_sys=True
     """
-    sys = ase_to_pyscal_system(structure)
+    sys = ase_to_pyscal(structure)
     sys.find_neighbors(method=neighbor_method, cutoff=cutoff)
     sys.find_solids(
         bonds=bonds,
diff --git a/structuretoolkit/analyse/spatial.py b/structuretoolkit/analyse/spatial.py
index f38df0422..05c71eecb 100644
--- a/structuretoolkit/analyse/spatial.py
+++ b/structuretoolkit/analyse/spatial.py
@@ -3,11 +3,10 @@
 # Distributed under the terms of "New BSD License", see the LICENSE file.
 
 import numpy as np
-from sklearn.cluster import AgglomerativeClustering, DBSCAN
 from scipy.sparse import coo_matrix
 from scipy.spatial import Voronoi, Delaunay
 from scipy.spatial import ConvexHull
-from structuretoolkit.helper import (
+from structuretoolkit.common.helper import (
     get_extended_positions,
     get_wrapped_coordinates,
     get_vertical_length,
@@ -261,6 +260,8 @@ def _kick_out_points(self, variance_buffer=0.01):
         self.positions = self.positions[variance < min_var + variance_buffer]
 
     def _cluster_points(self, eps=0.1):
+        from sklearn.cluster import DBSCAN
+
         if eps == 0:
             return
         extended_positions, indices = get_extended_positions(
@@ -428,6 +429,8 @@ def get_layers(
             "ij,i,nj->ni", mat, 1 / np.linalg.norm(mat, axis=-1), positions
         )
     if cluster_method is None:
+        from sklearn.cluster import AgglomerativeClustering
+
         cluster_method = AgglomerativeClustering(
             linkage="complete",
             n_clusters=None,
@@ -493,6 +496,8 @@ def get_voronoi_vertices(
     )
     xx = voro.vertices
     if distance_threshold > 0:
+        from sklearn.cluster import AgglomerativeClustering
+
         cluster = AgglomerativeClustering(
             linkage="single", distance_threshold=distance_threshold, n_clusters=None
         )
@@ -567,7 +572,7 @@ def get_delaunay_neighbors(structure, width_buffer: float = 10.0) -> np.ndarray:
     )
 
 
-def cluster_positions(
+def get_cluster_positions(
     structure, positions=None, eps=1, buffer_width=None, return_labels=False
 ):
     """
@@ -609,6 +614,8 @@ def cluster_positions(
         positions (numpy.ndarray): Mean positions
         label (numpy.ndarray): Labels of the positions (returned when `return_labels = True`)
     """
+    from sklearn.cluster import DBSCAN
+
     positions = structure.positions if positions is None else np.array(positions)
     if buffer_width is None:
         buffer_width = eps
diff --git a/structuretoolkit/analyse/strain.py b/structuretoolkit/analyse/strain.py
index eba6066f4..41c407689 100644
--- a/structuretoolkit/analyse/strain.py
+++ b/structuretoolkit/analyse/strain.py
@@ -1,7 +1,7 @@
 import numpy as np
 from scipy.spatial.transform import Rotation
 from structuretoolkit.analyse.neighbors import get_neighbors
-from structuretoolkit.analyse.pyscal import analyse_cna_adaptive
+from structuretoolkit.analyse.pyscal import get_adaptive_cna_descriptors
 
 
 class Strain:
@@ -138,7 +138,7 @@ def _get_best_match_indices(coords, ref_coord):
 
     @staticmethod
     def _get_majority_phase(structure):
-        cna = analyse_cna_adaptive(structure=structure)
+        cna = get_adaptive_cna_descriptors(structure=structure)
         return np.asarray([k for k in cna.keys()])[np.argmax([v for v in cna.values()])]
 
     @staticmethod
diff --git a/structuretoolkit/analyse/symmetry.py b/structuretoolkit/analyse/symmetry.py
index a1e2b6f31..3e4d7b653 100644
--- a/structuretoolkit/analyse/symmetry.py
+++ b/structuretoolkit/analyse/symmetry.py
@@ -7,8 +7,8 @@
 import spglib
 import ast
 
-import structuretoolkit.helper
-from structuretoolkit.helper import get_structure_indices
+import structuretoolkit.common.helper
+from structuretoolkit.common.error import SymmetryError
 
 __author__ = "Joerg Neugebauer, Sam Waseda"
 __copyright__ = (
@@ -22,10 +22,6 @@
 __date__ = "Sep 1, 2017"
 
 
-class SymmetryError(Exception):
-    pass
-
-
 class Symmetry(dict):
 
     """
@@ -243,11 +239,17 @@ def _get_spglib_cell(self, use_elements=None, use_magmoms=None):
             use_magmoms = self._use_magmoms
         if use_elements:
             numbers = np.array(
-                get_structure_indices(structure=self._structure), dtype="intc"
+                structuretoolkit.common.helper.get_structure_indices(
+                    structure=self._structure
+                ),
+                dtype="intc",
             )
         else:
             numbers = np.ones_like(
-                get_structure_indices(structure=self._structure), dtype="intc"
+                structuretoolkit.common.helper.get_structure_indices(
+                    structure=self._structure
+                ),
+                dtype="intc",
             )
         if use_magmoms:
             return (
@@ -354,7 +356,7 @@ def get_primitive_cell(
         new_structure = self._structure.copy()
         new_structure.cell = cell
         new_structure = new_structure[: len(indices)]
-        new_structure = structuretoolkit.helper.set_indices(
+        new_structure = structuretoolkit.common.helper.set_indices(
             structure=new_structure, indices=indices
         )
         new_structure.positions = positions
diff --git a/structuretoolkit/build/__init__.py b/structuretoolkit/build/__init__.py
index e69de29bb..fdc7eebf3 100644
--- a/structuretoolkit/build/__init__.py
+++ b/structuretoolkit/build/__init__.py
@@ -0,0 +1,10 @@
+from structuretoolkit.build.aimsgb import (
+    grainboundary,
+    get_grainboundary_info
+)
+from structuretoolkit.build.compound import B2, C14, C15, C36, D03
+from structuretoolkit.build.sqs import sqs_structures
+from structuretoolkit.build.surface import (
+    high_index_surface,
+    get_high_index_surface_info
+)
\ No newline at end of file
diff --git a/structuretoolkit/build/aimsgb.py b/structuretoolkit/build/aimsgb.py
index 176162c52..2842ca6a0 100644
--- a/structuretoolkit/build/aimsgb.py
+++ b/structuretoolkit/build/aimsgb.py
@@ -2,8 +2,7 @@
 # Copyright (c) Max-Planck-Institut für Eisenforschung GmbH - Computational Materials Design (CM) Department
 # Distributed under the terms of "New BSD License", see the LICENSE file.
 
-from aimsgb import GrainBoundary, Grain, GBInformation
-from pymatgen.io.ase import AseAtomsAdaptor
+from structuretoolkit.common.pymatgen import ase_to_pymatgen, pymatgen_to_ase
 
 __author__ = "Ujjal Saikia"
 __copyright__ = (
@@ -17,7 +16,7 @@
 __date__ = "Feb 26, 2021"
 
 
-def grainboundary_info(axis, max_sigma):
+def get_grainboundary_info(axis, max_sigma):
     """
     Provides a list of possible GB structures for a given rotational axis and upto the given maximum sigma value.
 
@@ -37,10 +36,11 @@ def grainboundary_info(axis, max_sigma):
     To construct the grain boundary select a GB plane and sigma value from the list and pass it to the
     GBBuilder.gb_build() function along with the rotational axis and initial bulk structure.
     """
+    from aimsgb import GBInformation
     return GBInformation(axis=axis, max_sigma=max_sigma)
 
 
-def grainboundary_build(
+def grainboundary(
     axis,
     sigma,
     plane,
@@ -76,8 +76,8 @@ def grainboundary_build(
     Returns:
         :class:`.Atoms`: final grain boundary structure
     """
-    adapter = AseAtomsAdaptor()
-    basis_pymatgen = adapter.get_structure(atoms=initial_struct, cls=None)
+    from aimsgb import GrainBoundary, Grain
+    basis_pymatgen = ase_to_pymatgen(structure=initial_struct)
     grain_init = Grain(
         basis_pymatgen.lattice, basis_pymatgen.species, basis_pymatgen.frac_coords
     )
@@ -90,7 +90,7 @@ def grainboundary_build(
         uc_b=uc_b,
     )
 
-    return adapter.get_atoms(
+    return pymatgen_to_ase(
         structure=gb_obj.build_gb(
             to_primitive=to_primitive,
             delete_layer=delete_layer,
diff --git a/structuretoolkit/build/sqs.py b/structuretoolkit/build/sqs.py
index 786dea08b..efd6136bc 100644
--- a/structuretoolkit/build/sqs.py
+++ b/structuretoolkit/build/sqs.py
@@ -5,7 +5,6 @@
 from ase.data import atomic_numbers
 from ase.atoms import Atoms
 import numpy as np
-from sqsgenerator import sqs_optimize
 from typing import Dict, Optional, Union, Iterable
 
 
@@ -101,7 +100,7 @@ def transpose(it):
     return zip(*it)
 
 
-def get_sqs_structures(
+def sqs_structures(
     structure: Atoms,
     mole_fractions: Dict[str, Union[float, int]],
     weights: Optional[Dict[int, float]] = None,
@@ -119,6 +118,7 @@ def get_sqs_structures(
     minimal: Optional[bool] = True,
     similar: Optional[bool] = True,
 ):
+    from sqsgenerator import sqs_optimize
     composition = mole_fractions_to_composition(mole_fractions, len(structure))
 
     settings = dict(
diff --git a/structuretoolkit/build/surface.py b/structuretoolkit/build/surface.py
index 3dc863bce..f460f5f38 100644
--- a/structuretoolkit/build/surface.py
+++ b/structuretoolkit/build/surface.py
@@ -1,11 +1,10 @@
 import numpy as np
 from ase.build import bulk, surface
-from structuretoolkit.analyse.symmetry import get_symmetry
-from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
-from pymatgen.io.ase import AseAtomsAdaptor
+from structuretoolkit.analyse import get_symmetry
+from structuretoolkit.common.pymatgen import ase_to_pymatgen, pymatgen_to_ase
 
 
-def high_index_surface_info(
+def get_high_index_surface_info(
         element,
         crystal_structure,
         lattice_constant,
@@ -136,13 +135,14 @@ def high_index_surface(
     Returns:
         slab: ase.atoms.Atoms instance Required surface
     """
+    from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
     basis = bulk(
         name=element,
         crystalstructure=crystal_structure,
         a=lattice_constant,
         cubic=True
     )
-    high_index_surface, _, _ = high_index_surface_info(
+    high_index_surface, _, _ = get_high_index_surface_info(
         element=element,
         crystal_structure=crystal_structure,
         lattice_constant=lattice_constant,
@@ -155,10 +155,9 @@ def high_index_surface(
         length_kink=length_kink,
     )
     surf = surface(basis, high_index_surface, layers, vacuum)
-    adapter = AseAtomsAdaptor()
-    sga = SpacegroupAnalyzer(adapter.get_structure(atoms=surf))
-    pmg_refined = sga.get_refined_structure()
-    slab = adapter.get_atoms(structure=pmg_refined)
+    slab = pymatgen_to_ase(
+        SpacegroupAnalyzer(ase_to_pymatgen(structure=surf)).get_refined_structure()
+    )
     slab.positions[:, 2] = slab.positions[:, 2] - np.min(slab.positions[:, 2])
     slab.set_pbc(True)
     return slab
diff --git a/structuretoolkit/common/__init__.py b/structuretoolkit/common/__init__.py
new file mode 100644
index 000000000..51e3a74f2
--- /dev/null
+++ b/structuretoolkit/common/__init__.py
@@ -0,0 +1,12 @@
+from structuretoolkit.common.pymatgen import ase_to_pymatgen, pymatgen_to_ase
+from structuretoolkit.common.pyscal import ase_to_pyscal
+from structuretoolkit.common.helper import (
+    get_atomic_numbers,
+    get_extended_positions,
+    get_vertical_length,
+    get_wrapped_coordinates,
+    select_index,
+    center_coordinates_in_unit_cell,
+    apply_strain,
+)
+from structuretoolkit.common.error import SymmetryError
diff --git a/structuretoolkit/common/error.py b/structuretoolkit/common/error.py
new file mode 100644
index 000000000..1fb557ccc
--- /dev/null
+++ b/structuretoolkit/common/error.py
@@ -0,0 +1,2 @@
+class SymmetryError(Exception):
+    pass
diff --git a/structuretoolkit/helper.py b/structuretoolkit/common/helper.py
similarity index 100%
rename from structuretoolkit/helper.py
rename to structuretoolkit/common/helper.py
diff --git a/structuretoolkit/common/pymatgen.py b/structuretoolkit/common/pymatgen.py
new file mode 100644
index 000000000..dffea57e4
--- /dev/null
+++ b/structuretoolkit/common/pymatgen.py
@@ -0,0 +1,12 @@
+def ase_to_pymatgen(structure):
+    from pymatgen.io.ase import AseAtomsAdaptor
+
+    adapter = AseAtomsAdaptor()
+    return adapter.get_structure(atoms=structure)
+
+
+def pymatgen_to_ase(structure):
+    from pymatgen.io.ase import AseAtomsAdaptor
+
+    adapter = AseAtomsAdaptor()
+    return adapter.get_atoms(structure=structure)
diff --git a/structuretoolkit/common/pyscal.py b/structuretoolkit/common/pyscal.py
new file mode 100644
index 000000000..9a010964d
--- /dev/null
+++ b/structuretoolkit/common/pyscal.py
@@ -0,0 +1,19 @@
+def ase_to_pyscal(structure):
+    """
+    Converts atoms to ase atoms and than to a pyscal system.
+    Also adds the pyscal publication.
+
+    Args:
+        structure (ase.atoms.Atoms): Structure to convert.
+
+    Returns:
+        Pyscal system: See the pyscal documentation.
+    """
+    import pyscal.core as pc
+
+    sys = pc.System()
+    sys.read_inputfile(
+        filename=structure,
+        format="ase",
+    )
+    return sys
diff --git a/structuretoolkit/visualize.py b/structuretoolkit/visualize.py
index 921f9afe3..1707981d0 100644
--- a/structuretoolkit/visualize.py
+++ b/structuretoolkit/visualize.py
@@ -4,9 +4,8 @@
 
 import numpy as np
 import warnings
-from matplotlib.colors import rgb2hex
 from scipy.interpolate import interp1d
-from structuretoolkit.helper import get_atomic_numbers
+from structuretoolkit.common.helper import get_atomic_numbers
 
 __author__ = "Joerg Neugebauer, Sudarsan Surendralal"
 __copyright__ = (
@@ -652,6 +651,8 @@ def _scalars_to_hex_colors(scalar_field, start=None, end=None, cmap=None):
     Returns:
         (list): The corresponding hex codes for each scalar value passed in.
     """
+    from matplotlib.colors import rgb2hex
+
     if start is None:
         start = np.amin(scalar_field)
     if end is None:
diff --git a/tests/test_aimsgb.py b/tests/test_aimsgb.py
index 4d9d76f9e..f4ab57cc1 100644
--- a/tests/test_aimsgb.py
+++ b/tests/test_aimsgb.py
@@ -16,9 +16,9 @@ def test_grain_thickness(self):
         axis = [0, 0, 1]
         sigma = 5
         plane = [1, 2, 0]
-        gb1 = stk.grainboundary_build(axis, sigma, plane, self.fcc_basis)  # Default thicknesses expected to be 1
+        gb1 = stk.build.grainboundary(axis, sigma, plane, self.fcc_basis)  # Default thicknesses expected to be 1
         uc_a, uc_b = 2, 3  # Make grains thicker
-        gb2 = stk.grainboundary_build(axis, sigma, plane, self.fcc_basis, uc_a=uc_a, uc_b=uc_b)
+        gb2 = stk.build.grainboundary(axis, sigma, plane, self.fcc_basis, uc_a=uc_a, uc_b=uc_b)
         self.assertEqual(
             ((uc_a + uc_b)/2)*len(gb1), len(gb2),
             msg="Expected structure to be bigger in proportion to grain thickness"
diff --git a/tests/test_analyse.py b/tests/test_analyse.py
index 8c12eb392..a9c3dc337 100644
--- a/tests/test_analyse.py
+++ b/tests/test_analyse.py
@@ -18,25 +18,25 @@ def test_get_layers(self):
         a_0 = 4
         struct = bulk(name='Al', a=a_0, crystalstructure='fcc', cubic=True).repeat(10)
         struct_pure = struct.copy()
-        layers = stk.get_layers(structure=struct)
+        layers = stk.analyse.get_layers(structure=struct)
         self.assertAlmostEqual(np.linalg.norm(layers-np.rint(2*struct.positions/a_0).astype(int)), 0)
         struct.append(Atom(symbol='C', position=np.random.random(3)))
         self.assertEqual(
-            np.linalg.norm(layers-stk.get_layers(structure=struct, id_list=stk.select_index(structure=struct, element='Al'))), 0
+            np.linalg.norm(layers-stk.analyse.get_layers(structure=struct, id_list=stk.common.select_index(structure=struct, element='Al'))), 0
         )
         self.assertEqual(
-            np.linalg.norm(layers-stk.get_layers(
+            np.linalg.norm(layers-stk.analyse.get_layers(
                 structure=struct,
-                id_list=stk.select_index(structure=struct, element='Al'),
+                id_list=stk.common.select_index(structure=struct, element='Al'),
                 wrap_atoms=False
             )), 0
         )
         with self.assertRaises(ValueError):
-            _ = stk.get_layers(structure=struct, distance_threshold=0)
+            _ = stk.analyse.get_layers(structure=struct, distance_threshold=0)
         with self.assertRaises(ValueError):
-            _ = stk.get_layers(structure=stk, id_list=[])
+            _ = stk.analyse.get_layers(structure=struct, id_list=[])
 
-        self.assertTrue(np.all(stk.get_layers(structure=struct) == stk.get_layers(
+        self.assertTrue(np.all(stk.analyse.get_layers(structure=struct) == stk.analyse.get_layers(
             structure=struct,
             cluster_method=AgglomerativeClustering(
                 linkage='complete',
@@ -46,7 +46,7 @@ def test_get_layers(self):
         ), "Overriding cluster method with default parameters does not return the same results.")
         self.assertTrue(
             np.all(
-                stk.get_layers(structure=struct_pure) == stk.get_layers(
+                stk.analyse.get_layers(structure=struct_pure) == stk.analyse.get_layers(
                     structure=struct_pure,
                     cluster_method=DBSCAN(eps=0.01)
                 )
@@ -56,57 +56,57 @@ def test_get_layers(self):
 
     def test_get_layers_other_planes(self):
         structure = bulk(name='Fe', a=3.5, crystalstructure='fcc', cubic=True).repeat(2)
-        layers = stk.get_layers(structure=structure, planes=[1, 1, 1])
+        layers = stk.analyse.get_layers(structure=structure, planes=[1, 1, 1])
         self.assertEqual(np.unique(layers).tolist(), [0, 1, 2, 3, 4])
 
     def test_get_layers_with_strain(self):
         structure = bulk(name='Fe', a=2.8, crystalstructure='bcc', cubic=True).repeat(2)
-        layers = stk.get_layers(structure=structure).tolist()
-        stk.apply_strain(structure=structure, epsilon=0.1*(np.random.random((3, 3))-0.5))
+        layers = stk.analyse.get_layers(structure=structure).tolist()
+        stk.common.apply_strain(structure=structure, epsilon=0.1*(np.random.random((3, 3))-0.5))
         self.assertEqual(
-            layers, stk.get_layers(structure=structure, planes=np.linalg.inv(structure.cell).T).tolist()
+            layers, stk.analyse.get_layers(structure=structure, planes=np.linalg.inv(structure.cell).T).tolist()
         )
 
     def test_get_layers_across_pbc(self):
         structure = bulk(name='Fe', a=2.8, crystalstructure='bcc', cubic=True).repeat(2)
-        layers = stk.get_layers(structure=structure)
+        layers = stk.analyse.get_layers(structure=structure)
         structure.cell[1, 0] += 0.01
-        structure = stk.center_coordinates_in_unit_cell(structure=structure)
-        self.assertEqual(len(np.unique(layers[stk.get_layers(structure=structure)[:, 0] == 0, 0])), 1)
+        structure = stk.common.center_coordinates_in_unit_cell(structure=structure)
+        self.assertEqual(len(np.unique(layers[stk.analyse.get_layers(structure=structure)[:, 0] == 0, 0])), 1)
 
     def test_pyscal_cna_adaptive(self):
         basis = Atoms(
             "FeFe", scaled_positions=[(0, 0, 0), (0.5, 0.5, 0.5)], cell=np.identity(3)
         )
         self.assertTrue(
-            stk.analyse_cna_adaptive(structure=basis)["bcc"] == 2
+            stk.analyse.get_adaptive_cna_descriptors(structure=basis)["bcc"] == 2
         )
 
     def test_pyscal_centro_symmetry(self):
         basis = bulk(name='Fe', a=2.8, crystalstructure='bcc', cubic=True)
         self.assertTrue(
-            all([np.isclose(v, 0.0) for v in stk.analyse_centro_symmetry(structure=basis, num_neighbors=8)])
+            all([np.isclose(v, 0.0) for v in stk.analyse.get_centro_symmetry_descriptors(structure=basis, num_neighbors=8)])
         )
 
     def test_get_voronoi_vertices(self):
         basis = bulk(name='Al', a=4, crystalstructure='fcc', cubic=True)
-        self.assertEqual(len(stk.get_voronoi_vertices(structure=basis)), 12)
-        self.assertEqual(len(stk.get_voronoi_vertices(structure=basis, distance_threshold=2)), 1)
+        self.assertEqual(len(stk.analyse.get_voronoi_vertices(structure=basis)), 12)
+        self.assertEqual(len(stk.analyse.get_voronoi_vertices(structure=basis, distance_threshold=2)), 1)
 
     def test_get_interstitials_bcc(self):
         bcc = bulk('Fe', cubic=True)
         x_octa_ref = bcc.positions[:, None, :]+0.5*bcc.cell[None, :, :]
         x_octa_ref = x_octa_ref.reshape(-1, 3)
-        x_octa_ref = stk.get_wrapped_coordinates(structure=bcc, positions=x_octa_ref)
-        int_octa = stk.get_interstitials(structure=bcc, num_neighbors=6)
+        x_octa_ref = stk.common.get_wrapped_coordinates(structure=bcc, positions=x_octa_ref)
+        int_octa = stk.analyse.get_interstitials(structure=bcc, num_neighbors=6)
         self.assertEqual(len(int_octa.positions), len(x_octa_ref))
         self.assertAlmostEqual(
             np.linalg.norm(
                 x_octa_ref[:, None, :]-int_octa.positions[None, :, :], axis=-1
             ).min(axis=0).sum(), 0
         )
-        int_tetra = stk.get_interstitials(structure=bcc, num_neighbors=4)
-        x_tetra_ref = stk.get_wrapped_coordinates(structure=bcc, positions=stk.get_voronoi_vertices(structure=bcc))
+        int_tetra = stk.analyse.get_interstitials(structure=bcc, num_neighbors=4)
+        x_tetra_ref = stk.common.get_wrapped_coordinates(structure=bcc, positions=stk.get_voronoi_vertices(structure=bcc))
         self.assertEqual(len(int_tetra.positions), len(x_tetra_ref))
         self.assertAlmostEqual(
             np.linalg.norm(
@@ -118,8 +118,8 @@ def test_get_interstitials_fcc(self):
         fcc = bulk('Al', cubic=True)
         a_0 = fcc.cell[0, 0]
         x_tetra_ref = 0.25*a_0*np.ones(3)*np.array([[1], [-1]])+fcc.positions[:, None, :]
-        x_tetra_ref = stk.get_wrapped_coordinates(structure=fcc, positions=x_tetra_ref).reshape(-1, 3)
-        int_tetra = stk.get_interstitials(structure=fcc, num_neighbors=4)
+        x_tetra_ref = stk.common.get_wrapped_coordinates(structure=fcc, positions=x_tetra_ref).reshape(-1, 3)
+        int_tetra = stk.analyse.get_interstitials(structure=fcc, num_neighbors=4)
         self.assertEqual(len(int_tetra.positions), len(x_tetra_ref))
         self.assertAlmostEqual(
             np.linalg.norm(
@@ -127,8 +127,8 @@ def test_get_interstitials_fcc(self):
             ).min(axis=0).sum(), 0
         )
         x_octa_ref = 0.5*a_0*np.array([1, 0, 0])+fcc.positions
-        x_octa_ref = stk.get_wrapped_coordinates(structure=fcc, positions=x_octa_ref)
-        int_octa = stk.get_interstitials(structure=fcc, num_neighbors=6)
+        x_octa_ref = stk.common.get_wrapped_coordinates(structure=fcc, positions=x_octa_ref)
+        int_octa = stk.analyse.get_interstitials(structure=fcc, num_neighbors=6)
         self.assertEqual(len(int_octa.positions), len(x_octa_ref))
         self.assertAlmostEqual(
             np.linalg.norm(x_octa_ref[:, None, :]-int_octa.positions[None, :, :], axis=-1).min(axis=0).sum(), 0
@@ -166,11 +166,11 @@ def test_strain(self):
         voro = Voronoi(structure.positions[:, :2])
         center = voro.vertices[np.linalg.norm(voro.vertices-structure.cell.diagonal()[:2]*0.5, axis=-1).argmin()]
         structure.positions[:, 2] += b/(2*np.pi)*np.arctan2(*(structure.positions[:, :2]-center).T[::-1])
-        structure = stk.center_coordinates_in_unit_cell(structure=structure)
+        structure = stk.common.center_coordinates_in_unit_cell(structure=structure)
         r_0 = 0.9*L/2
         r = np.linalg.norm(structure.positions[:, :2]-center, axis=-1)
         core_region = (r < r_0)*(r > 10)
-        strain = stk.get_strain(structure=structure, ref_structure=structure_bulk, num_neighbors=8)
+        strain = stk.analyse.get_strain(structure=structure, ref_structure=structure_bulk, num_neighbors=8)
         strain = strain[core_region]
         positions = structure.positions[core_region, :2]
         x = positions-center
@@ -183,20 +183,20 @@ def test_tessellations(self):
         structure_bulk = bulk('Fe', cubic=True)
         a_0 = structure_bulk.cell[0, 0]
         structure = structure_bulk.repeat(3)
-        self.assertAlmostEqual(np.linalg.norm(stk.find_mic(structure=structure, v=np.diff(
-            structure.positions[stk.get_delaunay_neighbors(structure=structure)], axis=-2
+        self.assertAlmostEqual(np.linalg.norm(stk.analyse.find_mic(structure=structure, v=np.diff(
+            structure.positions[stk.analyse.get_delaunay_neighbors(structure=structure)], axis=-2
         )), axis=-1).flatten().max(), a_0)
-        self.assertAlmostEqual(np.linalg.norm(stk.find_mic(structure=structure, v=np.diff(
-            structure.positions[stk.get_voronoi_neighbors(structure=structure)], axis=-2
+        self.assertAlmostEqual(np.linalg.norm(stk.analyse.find_mic(structure=structure, v=np.diff(
+            structure.positions[stk.analyse.get_voronoi_neighbors(structure=structure)], axis=-2
         )), axis=-1).flatten().max(), a_0)
 
     def test_cluster_positions(self):
         structure_bulk = bulk('Fe', cubic=True)
-        self.assertEqual(len(stk.cluster_positions(structure=structure_bulk)), len(structure_bulk))
+        self.assertEqual(len(stk.analyse.get_cluster_positions(structure=structure_bulk)), len(structure_bulk))
         positions = np.append(structure_bulk.positions, structure_bulk.positions, axis=0)
-        self.assertEqual(len(stk.cluster_positions(structure=structure_bulk, positions=positions)), len(structure_bulk))
+        self.assertEqual(len(stk.analyse.get_cluster_positions(structure=structure_bulk, positions=positions)), len(structure_bulk))
         self.assertEqual(
-            stk.cluster_positions(structure=structure_bulk, positions=np.zeros((2, 3)), return_labels=True)[1].tolist(),
+            stk.analyse.get_cluster_positions(structure=structure_bulk, positions=np.zeros((2, 3)), return_labels=True)[1].tolist(),
             [0, 0]
         )
 
diff --git a/tests/test_compound.py b/tests/test_compound.py
index bf22c3f91..4c8f44b75 100644
--- a/tests/test_compound.py
+++ b/tests/test_compound.py
@@ -11,23 +11,23 @@
 
 class TestCompound(unittest.TestCase):
     def test_B2(self):
-        structure = stk.B2('Fe', 'Al')
+        structure = stk.build.B2('Fe', 'Al')
         self.assertAlmostEqual(bulk('Fe', cubic=True).cell[0, 0], structure.cell[0, 0],
                                msg="Docstring claims lattice constant defaults to primary species")
         self.assertEqual(2, len(structure))
-        neigh = stk.get_neighbors(structure=structure, num_neighbors=8)
+        neigh = stk.analyse.get_neighbors(structure=structure, num_neighbors=8)
         symbols = np.array(structure.get_chemical_symbols())
         self.assertEqual(8, np.sum(symbols[neigh.indices[0]] == 'Al'),
                          msg="Expected the primary atom to have all secondary neighbors")
         self.assertEqual(8, np.sum(symbols[neigh.indices[1]] == 'Fe'),
                          msg="Expected the secondary atom to have all primary neighbors")
-        structure = stk.B2('Fe', 'Al', a=1)
+        structure = stk.build.B2('Fe', 'Al', a=1)
         self.assertTrue(np.allclose(np.diag(structure.cell.array), 1), "Expected cubic cell with specified size.")
 
     def test_C14(self):
         a_type = 'Mg'
         b_type = 'Cu'
-        c14 = stk.C14(a_type, b_type)
+        c14 = stk.build.C14(a_type, b_type)
         self.assertEqual(len(c14), 12, "Wrong number of atoms in C14 structure.")
         self.assertEqual(c14.get_chemical_formula(), "Cu8Mg4", "Wrong chemical formula.")
 
@@ -38,29 +38,29 @@ def test_C15(self):
 
         a_type = 'Mg'
         b_type = 'Cu'
-        structure = stk.C15(a_type, b_type)
+        structure = stk.build.C15(a_type, b_type)
 
         self.assertEqual(len(structure), 24, "Wrong number of atoms in C15 structure.")
         self.assertEqual(structure.get_chemical_formula(), "Cu16Mg8", "Wrong chemical formula.")
 
-        a_type_nn_distance = stk.get_neighbors(structure=bulk(a_type), num_neighbors=1).distances[0, 0]
+        a_type_nn_distance = stk.analyse.get_neighbors(structure=bulk(a_type), num_neighbors=1).distances[0, 0]
         self.assertAlmostEqual((4 / np.sqrt(3)) * a_type_nn_distance, structure.cell.array[0, 0],
                                msg="Default lattice constant should relate to NN distance of A-type element.")
 
-        unique_ids = np.unique(stk.get_symmetry(structure)['equivalent_atoms'])
+        unique_ids = np.unique(stk.analyse.get_symmetry(structure)['equivalent_atoms'])
         self.assertEqual(2, len(unique_ids), msg="Expected only A- and B1-type sites.")
         symbols = np.array(structure.get_chemical_symbols())
         a_id = unique_ids[np.argwhere(symbols[unique_ids] == a_type)[0, 0]]
         b_id = unique_ids[np.argwhere(symbols[unique_ids] == b_type)[0, 0]]
         unique_ids = [a_id, b_id]  # Now with guaranteed ordering
 
-        csa = stk.analyse_centro_symmetry(structure)[unique_ids]
+        csa = stk.analyse.get_centro_symmetry_descriptors(structure)[unique_ids]
         self.assertLess(1, csa[0], msg="A site for AB_2 C15 should be significantly non-centro-symmetric.")
         self.assertAlmostEqual(0, csa[1], msg="B site for AB_2 C15 should be nearly centro-symmetric.")
 
         num_a_neighs = 16
         num_b_neighs = 12
-        neigh = stk.get_neighbors(structure=structure, num_neighbors=num_a_neighs)
+        neigh = stk.analyse.get_neighbors(structure=structure, num_neighbors=num_a_neighs)
         a_neighs = neigh.indices[unique_ids[0]]
         b_neighs = neigh.indices[unique_ids[1], :num_b_neighs]
         symbols = np.array(structure.get_chemical_symbols())
@@ -72,15 +72,15 @@ def test_C15(self):
     def test_C36(self):
         a_type = 'Mg'
         b_type = 'Cu'
-        c36 = stk.C36(a_type, b_type)
+        c36 = stk.build.C36(a_type, b_type)
         self.assertEqual(len(c36), 24, "Wrong number of atoms in C36 structure.")
         self.assertEqual(c36.get_chemical_formula(), "Cu16Mg8", "Wrong chemical formula.")
 
     def test_D03(self):
         element_a, element_b = 'Al', 'Fe'
-        structure = stk.D03(element_a, element_b)
+        structure = stk.build.D03(element_a, element_b)
         symbols = np.array(structure.get_chemical_symbols())
-        neigh = stk.get_neighbors(structure=structure, num_neighbors=8)
+        neigh = stk.analyse.get_neighbors(structure=structure, num_neighbors=8)
 
         a_neighbors = neigh.indices[symbols == element_a]
         self.assertTrue(np.all(symbols[a_neighbors] == element_b), msg="A-type should only have B-type neighbors.")
diff --git a/tests/test_high_index_surface.py b/tests/test_high_index_surface.py
index 812a5a91b..d9b569909 100644
--- a/tests/test_high_index_surface.py
+++ b/tests/test_high_index_surface.py
@@ -4,7 +4,7 @@
 
 class TestHighIndexSurface(unittest.TestCase):
     def test_high_index_surface(self):
-        slab = stk.high_index_surface(
+        slab = stk.build.high_index_surface(
             element='Ni',
             crystal_structure='fcc',
             lattice_constant=3.526,
@@ -21,7 +21,7 @@ def test_high_index_surface(self):
         self.assertEqual(len(slab), 60)
 
     def test_high_index_surface_info(self):
-        h, s, k = stk.high_index_surface_info(
+        h, s, k = stk.build.get_high_index_surface_info(
             element='Ni',
             crystal_structure='fcc',
             lattice_constant=3.526,
@@ -38,7 +38,7 @@ def test_high_index_surface_info(self):
         self.assertEqual(len(k), 3)
         self.assertEqual(len(s), 3)
         with self.assertRaises(ValueError):
-            stk.high_index_surface_info(
+            stk.build.get_high_index_surface_info(
                 element='Ni',
                 crystal_structure='fcc',
                 lattice_constant=3.526,
diff --git a/tests/test_neighbors.py b/tests/test_neighbors.py
index 2926e616b..8cb1a34ec 100644
--- a/tests/test_neighbors.py
+++ b/tests/test_neighbors.py
@@ -22,7 +22,7 @@ def setUpClass(cls):
     def test_allow_ragged(self):
         struct = bulk('Al', a=4, cubic=True).repeat(10)
         del struct[0]
-        neigh = stk.get_neighbors(structure=struct, cutoff_radius=3, mode="ragged")
+        neigh = stk.analyse.get_neighbors(structure=struct, cutoff_radius=3, mode="ragged")
         self.assertTrue(neigh._mode["ragged"])
         self.assertTrue(isinstance(neigh.indices, list))
         indices = neigh.indices.copy()
@@ -33,14 +33,14 @@ def test_allow_ragged(self):
             struct.positions[neigh.indices] = struct.positions[neigh.indices]
         neigh._set_mode(neigh._allow_ragged_to_mode(True))
         self.assertTrue(np.array_equal(neigh.indices[0], indices[0]))
-        neigh = stk.get_neighbors(structure=struct, cutoff_radius=3, num_neighbors=None)
+        neigh = stk.analyse.get_neighbors(structure=struct, cutoff_radius=3, num_neighbors=None)
         self.assertFalse(neigh._mode["ragged"])
         self.assertTrue(isinstance(neigh.indices, np.ndarray))
 
     def test_getter_and_ragged(self):
         struct = bulk('Al', a=4, cubic=True).repeat(2)
         del struct[0]
-        neigh = stk.get_neighbors(structure=struct, cutoff_radius=3, mode='filled')
+        neigh = stk.analyse.get_neighbors(structure=struct, cutoff_radius=3, mode='filled')
         vecs = neigh.filled.vecs
         distances = neigh.filled.distances
         indices = neigh.filled.indices
@@ -50,7 +50,7 @@ def test_getter_and_ragged(self):
 
     def test_get_neighborhood_single(self):
         struct = bulk('Al', a=4, cubic=True)
-        neigh = stk.get_neighborhood(structure=struct, positions=np.random.random(3), cutoff_radius=3)
+        neigh = stk.analyse.get_neighborhood(structure=struct, positions=np.random.random(3), cutoff_radius=3)
         distances = neigh.distances.copy()
         neigh.allow_ragged = False
         self.assertTrue(np.array_equal(neigh.distances, distances))
@@ -65,7 +65,7 @@ def test_get_neighbors(self):
         cell += np.random.random((3, 3))-0.5
         struct.positions += np.random.random((len(struct), 3))-0.5
         struct.set_cell(cell, scale_atoms=True)
-        neigh = stk.get_neighbors(structure=struct)
+        neigh = stk.analyse.get_neighbors(structure=struct)
         self.assertAlmostEqual(
             np.absolute(neigh.distances-np.linalg.norm(neigh.vecs, axis=-1)).max(), 0
         )
@@ -91,7 +91,7 @@ def test_pbc_false(self):
         cell = struct.cell.copy()
         cell += np.random.random((3, 3))-0.5
         struct.set_cell(cell, scale_atoms=True)
-        neigh = stk.get_neighbors(structure=struct)
+        neigh = stk.analyse.get_neighbors(structure=struct)
         self.assertAlmostEqual(
             np.absolute(neigh.distances-np.linalg.norm(neigh.vecs, axis=-1)).max(), 0
         )
@@ -104,7 +104,7 @@ def test_pbc_false(self):
 
     def test_fe_large(self):
         struct = bulk('Fe', a=2.85, cubic=True).repeat(10)
-        neigh = stk.get_neighbors(structure=struct)
+        neigh = stk.analyse.get_neighbors(structure=struct)
         self.assertAlmostEqual(
             np.absolute(neigh.distances-np.linalg.norm(neigh.vecs, axis=-1)).max(), 0
         )
@@ -113,14 +113,14 @@ def test_fe_large(self):
 
     def test_fe_small(self):
         struct = bulk('Fe', a=2.85, cubic=True)
-        neigh = stk.get_neighbors(structure=struct)
+        neigh = stk.analyse.get_neighbors(structure=struct)
         self.assertAlmostEqual(neigh.vecs[neigh.shells == 1].sum(), 0)
         with self.assertRaises(ValueError):
-            _ = stk.get_neighbors(structure=struct, num_neighbors=None)
+            _ = stk.analyse.get_neighbors(structure=struct, num_neighbors=None)
 
     def test_al_large(self):
         struct = bulk('Al', a=4.04, cubic=True).repeat(10)
-        neigh = stk.get_neighbors(structure=struct)
+        neigh = stk.analyse.get_neighbors(structure=struct)
         self.assertAlmostEqual(
             np.absolute(neigh.distances-np.linalg.norm(neigh.vecs, axis=-1)).max(), 0
         )
@@ -129,17 +129,17 @@ def test_al_large(self):
 
     def test_al_small(self):
         struct = bulk('Al', a=4.04, cubic=True)
-        neigh = stk.get_neighbors(structure=struct)
+        neigh = stk.analyse.get_neighbors(structure=struct)
         self.assertAlmostEqual(
             np.absolute(neigh.distances-np.linalg.norm(neigh.vecs, axis=-1)).max(), 0
         )
         self.assertAlmostEqual(neigh.vecs[neigh.shells == 1].sum(), 0)
         with self.assertRaises(ValueError):
-            stk.get_neighbors(structure=struct, num_neighbors=0)
+            stk.analyse.get_neighbors(structure=struct, num_neighbors=0)
 
     def test_wrapped_positions(self):
         structure = bulk('Al', a=4, cubic=True).repeat(2)
-        neigh = stk.get_neighbors(structure=structure)
+        neigh = stk.analyse.get_neighbors(structure=structure)
         distances = neigh.distances
         new_positions = structure.positions+structure.cell.diagonal()*2
         self.assertFalse(
@@ -152,17 +152,17 @@ def test_wrapped_positions(self):
 
     def test_get_global_shells(self):
         structure = bulk('Al', a=4, cubic=True).repeat(2)
-        neigh = stk.get_neighbors(structure=structure)
+        neigh = stk.analyse.get_neighbors(structure=structure)
         self.assertTrue(np.array_equal(neigh.shells, neigh.get_global_shells()))
         structure += Atoms(symbols='C', positions=[[0, 0, 0.5*4]])
-        neigh = stk.get_neighbors(structure=structure)
+        neigh = stk.analyse.get_neighbors(structure=structure)
         self.assertFalse(np.array_equal(neigh.shells, neigh.get_global_shells()))
         structure = bulk('Al', a=4, cubic=True).repeat(2)
-        neigh = stk.get_neighbors(structure=structure)
+        neigh = stk.analyse.get_neighbors(structure=structure)
         shells = neigh.get_global_shells()
         structure.positions += 0.01*(np.random.random((len(structure), 3))-0.5)
-        stk.center_coordinates_in_unit_cell(structure=structure)
-        neigh = stk.get_neighbors(structure=structure)
+        stk.common.center_coordinates_in_unit_cell(structure=structure)
+        neigh = stk.analyse.get_neighbors(structure=structure)
         self.assertTrue(
             np.array_equal(
                 shells, neigh.get_global_shells(cluster_by_vecs=True, cluster_by_distances=True)
@@ -174,10 +174,10 @@ def test_get_global_shells(self):
 
     def test_get_local_shells(self):
         structure = bulk('Al', a=4, cubic=True).repeat(2)
-        neigh = stk.get_neighbors(structure=structure)
+        neigh = stk.analyse.get_neighbors(structure=structure)
         shells = neigh.get_local_shells()
         structure.positions += 0.01*(np.random.random((len(structure), 3))-0.5)
-        neigh = stk.get_neighbors(structure=structure)
+        neigh = stk.analyse.get_neighbors(structure=structure)
         self.assertTrue(
             np.array_equal(
                 shells, neigh.get_local_shells(cluster_by_vecs=True, cluster_by_distances=True)
@@ -190,7 +190,7 @@ def test_get_local_shells(self):
     def test_get_global_shells_ragged(self):
         structure = bulk('Al', a=4, cubic=True).repeat(2)
         del structure[0]
-        neigh = stk.get_neighbors(structure=structure, cutoff_radius=3.5, num_neighbors=None)
+        neigh = stk.analyse.get_neighbors(structure=structure, cutoff_radius=3.5, num_neighbors=None)
         self.assertEqual(np.sum(neigh.get_global_shells() == -1), 12)
         self.assertEqual(np.sum(neigh.get_global_shells(cluster_by_distances=True) == -1), 12)
         self.assertEqual(np.sum(neigh.get_global_shells(cluster_by_vecs=True) == -1), 12)
@@ -206,14 +206,14 @@ def test_get_global_shells_ragged(self):
     def test_get_local_shells_ragged(self):
         structure = bulk('Al', a=4, cubic=True).repeat(2)
         del structure[0]
-        neigh = stk.get_neighbors(structure=structure, cutoff_radius=3.5, num_neighbors=None)
+        neigh = stk.analyse.get_neighbors(structure=structure, cutoff_radius=3.5, num_neighbors=None)
         self.assertEqual(np.sum(neigh.shells == -1), 12)
         self.assertEqual(np.sum(neigh.get_local_shells(cluster_by_distances=True) == -1), 12)
         self.assertEqual(np.sum(neigh.get_local_shells(cluster_by_vecs=True) == -1), 12)
         self.assertEqual(
             np.sum(neigh.get_local_shells(cluster_by_distances=True, cluster_by_vecs=True) == -1), 12
         )
-        neigh = stk.get_neighbors(structure=structure, cutoff_radius=3.5, num_neighbors=None, mode='ragged')
+        neigh = stk.analyse.get_neighbors(structure=structure, cutoff_radius=3.5, num_neighbors=None, mode='ragged')
         self.assertEqual(np.sum([len(s) == 11 for s in neigh.shells]), 12)
         self.assertEqual(np.sum([len(s) == 11 for s in neigh.get_local_shells(cluster_by_distances=True)]), 12)
         self.assertEqual(np.sum([len(s) == 11 for s in neigh.get_local_shells(cluster_by_vecs=True)]), 12)
@@ -223,14 +223,14 @@ def test_get_shells_flattened(self):
         structure = bulk('Al', cubic=True).repeat(2)
         del structure[0]
         r = structure.cell[0, 0]*0.49
-        neigh = stk.get_neighbors(structure=structure, cutoff_radius=r, num_neighbors=None, mode='flattened')
+        neigh = stk.analyse.get_neighbors(structure=structure, cutoff_radius=r, num_neighbors=None, mode='flattened')
         self.assertEqual(len(np.unique(neigh.shells)), 1)
         self.assertEqual(len(neigh.shells), 360)
         self.assertEqual(len(np.unique(neigh.get_local_shells())), 1)
         self.assertEqual(len(neigh.get_local_shells()), 360)
         self.assertEqual(len(np.unique(neigh.get_global_shells())), 1)
         self.assertEqual(len(neigh.get_global_shells()), 360)
-        neigh = stk.get_neighbors(structure=structure, cutoff_radius=r, num_neighbors=None)
+        neigh = stk.analyse.get_neighbors(structure=structure, cutoff_radius=r, num_neighbors=None)
         self.assertEqual(len(np.unique(neigh.flattened.shells)), 1)
         self.assertEqual(len(neigh.flattened.shells), 360)
 
@@ -238,7 +238,7 @@ def test_get_distances_flattened(self):
         structure = bulk('Al', cubic=True).repeat(2)
         del structure[0]
         r = structure.cell[0, 0]*0.49
-        neigh = stk.get_neighbors(structure=structure, cutoff_radius=r, num_neighbors=None, mode='flattened')
+        neigh = stk.analyse.get_neighbors(structure=structure, cutoff_radius=r, num_neighbors=None, mode='flattened')
         self.assertAlmostEqual(np.std(neigh.distances), 0)
         self.assertEqual(len(neigh.distances), 360)
         self.assertEqual(neigh.vecs.shape, (360, 3, ))
@@ -247,13 +247,13 @@ def test_atom_numbers(self):
         structure = bulk('Al', cubic=True).repeat(2)
         del structure[0]
         r = structure.cell[0, 0]*0.49
-        neigh = stk.get_neighbors(structure=structure, cutoff_radius=r, num_neighbors=None, mode='filled')
+        neigh = stk.analyse.get_neighbors(structure=structure, cutoff_radius=r, num_neighbors=None, mode='filled')
         n = len(structure)
         self.assertEqual(neigh.atom_numbers.sum(), int(n*(n-1)/2*12))
-        neigh = stk.get_neighbors(structure=structure, cutoff_radius=r, num_neighbors=None, mode='ragged')
+        neigh = stk.analyse.get_neighbors(structure=structure, cutoff_radius=r, num_neighbors=None, mode='ragged')
         for i, (a, d) in enumerate(zip(neigh.atom_numbers, neigh.distances)):
             self.assertEqual(np.sum(a-len(d)*[i]), 0)
-        neigh = stk.get_neighbors(structure=structure, cutoff_radius=r, num_neighbors=None, mode='flattened')
+        neigh = stk.analyse.get_neighbors(structure=structure, cutoff_radius=r, num_neighbors=None, mode='flattened')
         labels, counts = np.unique(
             np.unique(neigh.atom_numbers, return_counts=True)[1], return_counts=True
         )
@@ -263,7 +263,7 @@ def test_atom_numbers(self):
     def test_get_shell_matrix(self):
         structure = bulk('Fe', a=2.83, cubic=True).repeat(2)
         structure.symbols[0] = 'Ni'
-        neigh = stk.get_neighbors(structure=structure, num_neighbors=8)
+        neigh = stk.analyse.get_neighbors(structure=structure, num_neighbors=8)
         mat = neigh.get_shell_matrix()
         self.assertEqual(mat[0].sum(), 8*len(structure))
         mat = neigh.get_shell_matrix(chemical_pair=['Fe', 'Ni'])
@@ -273,7 +273,7 @@ def test_get_shell_matrix(self):
 
     def test_cluster_analysis(self):
         basis = bulk('Al', a=4.2, cubic=True).repeat(10)
-        neigh = stk.get_neighbors(structure=basis, num_neighbors=100)
+        neigh = stk.analyse.get_neighbors(structure=basis, num_neighbors=100)
         key, counts = neigh.cluster_analysis(id_list=[0, 1], return_cluster_sizes=True)
         self.assertTrue(np.array_equal(key[1], [0, 1]))
         self.assertEqual(counts[0], 2)
@@ -285,7 +285,7 @@ def test_cluster_analysis(self):
 
     def test_get_bonds(self):
         basis = bulk('Al', a=4.2, cubic=True).repeat(5)
-        neigh = stk.get_neighbors(structure=basis, num_neighbors=20)
+        neigh = stk.analyse.get_neighbors(structure=basis, num_neighbors=20)
         bonds = neigh.get_bonds()
         self.assertTrue(np.array_equal(np.sort(bonds[0]['Al'][0]),
                         np.sort(neigh.indices[0, neigh.shells[0] == 1])))
@@ -295,7 +295,7 @@ def test_find_neighbors_by_vector(self):
                       scaled_positions=[(0, 0, 0), (0.5, 0.5, 0.5)],
                       cell=np.identity(3),
                       pbc=True)
-        neigh = stk.get_neighbors(structure=basis, num_neighbors=14)
+        neigh = stk.analyse.get_neighbors(structure=basis, num_neighbors=14)
         id_lst, dist = neigh.find_neighbors_by_vector([0, 0, 1],
                                                       return_deviation=True)
         self.assertEqual(len(np.unique(np.unique(id_lst, return_counts=True)[1])), 1)
@@ -305,7 +305,7 @@ def test_find_neighbors_by_vector(self):
 
     def test_get_distances_arbitrary_array(self):
         basis = bulk('Al', a=4.2, cubic=True).repeat(3)
-        neigh = stk.get_neighbors(structure=basis, cutoff_radius=3.5, num_neighbors=None)
+        neigh = stk.analyse.get_neighbors(structure=basis, cutoff_radius=3.5, num_neighbors=None)
         self.assertEqual(len(neigh.get_neighborhood(np.random.random(3), num_neighbors=12).indices), 12)
         self.assertEqual(
             len(neigh.get_neighborhood(np.random.random(3), num_neighbors=12).ragged.distances), 12
@@ -314,7 +314,7 @@ def test_get_distances_arbitrary_array(self):
             len(neigh.get_neighborhood(np.random.random(3), num_neighbors=12, cutoff_radius=3.5).ragged.vecs), 12
         )
         self.assertTrue(neigh.get_neighborhood(np.random.random((2, 3)), num_neighbors=12).vecs.shape == (2, 12, 3))
-        neigh = stk.get_neighbors(structure=basis, num_neighbors=50)
+        neigh = stk.analyse.get_neighbors(structure=basis, num_neighbors=50)
         self.assertTrue(neigh.get_neighborhood(np.random.random(3)).distances.shape == (50,))
         self.assertTrue(neigh.get_neighborhood(np.random.random((2, 3))).indices.shape == (2, 50))
         self.assertTrue(neigh.get_neighborhood(np.random.random((2, 2, 3))).vecs.shape == (2, 2, 50, 3))
@@ -327,13 +327,13 @@ def test_get_distances_arbitrary_array(self):
 
     def test_repr(self):
         basis = bulk('Al', a=4.2, cubic=True).repeat(3)
-        neigh = stk.get_neighbors(structure=basis, cutoff_radius=3.5, num_neighbors=None)
+        neigh = stk.analyse.get_neighbors(structure=basis, cutoff_radius=3.5, num_neighbors=None)
         self.assertTrue('each atom' in neigh.__repr__())
 
     def test_norm_order(self):
         a_0 = 2.8
         basis = bulk('Fe', a=a_0, cubic=True).repeat(10)
-        neigh = stk.get_neighbors(structure=basis, num_neighbors=None, norm_order=np.inf, cutoff_radius=a_0+0.01)
+        neigh = stk.analyse.get_neighbors(structure=basis, num_neighbors=None, norm_order=np.inf, cutoff_radius=a_0+0.01)
         self.assertEqual(len(neigh.indices[0]), 34)
         with self.assertRaises(ValueError):
             neigh.norm_order = 3
@@ -341,25 +341,25 @@ def test_norm_order(self):
     def test_chemical_symbols(self):
         basis = bulk('Fe', cubic=True)
         basis.symbols[0] = 'Ni'
-        neigh = stk.get_neighbors(structure=basis, num_neighbors=1)
+        neigh = stk.analyse.get_neighbors(structure=basis, num_neighbors=1)
         self.assertEqual(neigh.chemical_symbols[0, 0], 'Fe')
         self.assertEqual(neigh.chemical_symbols[1, 0], 'Ni')
         vacancy = bulk('Fe', cubic=True).repeat(4)
         del vacancy[0]
-        neigh = stk.get_neighbors(structure=vacancy, num_neighbors=None, cutoff_radius=3)
+        neigh = stk.analyse.get_neighbors(structure=vacancy, num_neighbors=None, cutoff_radius=3)
         self.assertEqual(neigh.chemical_symbols[0, -1], 'v')
 
     def test_steinhardt_parameters(self):
-        neigh = stk.get_neighbors(structure=bulk('Al'), num_neighbors=12)
+        neigh = stk.analyse.get_neighbors(structure=bulk('Al'), num_neighbors=12)
         # values obtained with pyscal
         self.assertAlmostEqual(0, neigh.get_steinhardt_parameter(2)[0])
         self.assertAlmostEqual(0.19094065395649323, neigh.get_steinhardt_parameter(4)[0])
         self.assertAlmostEqual(0.5745242597140696, neigh.get_steinhardt_parameter(6)[0])
-        neigh = stk.get_neighbors(structure=bulk('Mg', a=1, c=np.sqrt(8/3)), num_neighbors=12)
+        neigh = stk.analyse.get_neighbors(structure=bulk('Mg', a=1, c=np.sqrt(8/3)), num_neighbors=12)
         self.assertAlmostEqual(0, neigh.get_steinhardt_parameter(2)[0])
         self.assertAlmostEqual(0.097222222, neigh.get_steinhardt_parameter(4)[0])
         self.assertAlmostEqual(0.484761685, neigh.get_steinhardt_parameter(6)[0])
-        neigh = stk.get_neighbors(structure=bulk('Fe'), num_neighbors=14)
+        neigh = stk.analyse.get_neighbors(structure=bulk('Fe'), num_neighbors=14)
         self.assertAlmostEqual(0.03636964837266537, neigh.get_steinhardt_parameter(4)[0])
         self.assertAlmostEqual(0.5106882308569508, neigh.get_steinhardt_parameter(6)[0])
         self.assertRaises(ValueError, neigh.get_steinhardt_parameter, 2, 2)
@@ -367,35 +367,35 @@ def test_steinhardt_parameters(self):
     def test_numbers_of_neighbors(self):
         basis = bulk('Al', cubic=True).repeat(2)
         del basis[0]
-        neigh = stk.get_neighbors(structure=basis, num_neighbors=None, cutoff_radius=0.45*basis.cell[0,0])
+        neigh = stk.analyse.get_neighbors(structure=basis, num_neighbors=None, cutoff_radius=0.45*basis.cell[0,0])
         n, c = np.unique(neigh.numbers_of_neighbors, return_counts=True)
         self.assertEqual(n.tolist(), [11, 12])
         self.assertEqual(c.tolist(), [12, 19])
 
     def test_modes(self):
         basis = bulk('Al', cubic=True)
-        neigh = stk.get_neighbors(structure=basis)
+        neigh = stk.analyse.get_neighbors(structure=basis)
         self.assertTrue(neigh.mode == 'filled')
         with self.assertRaises(KeyError):
-            neigh = stk.get_neighbors(structure=basis, mode='random_key')
+            neigh = stk.analyse.get_neighbors(structure=basis, mode='random_key')
 
     def test_centrosymmetry(self):
         structure = bulk('Fe').repeat(4)
-        cs = stk.get_neighbors(structure=structure, num_neighbors=8).centrosymmetry
+        cs = stk.analyse.get_neighbors(structure=structure, num_neighbors=8).centrosymmetry
         self.assertAlmostEqual(cs.max(), 0)
         self.assertAlmostEqual(cs.min(), 0)
         structure.positions += 0.01*(2*np.random.random(structure.positions.shape)-1)
-        neigh = stk.get_neighbors(structure=structure, num_neighbors=8)
+        neigh = stk.analyse.get_neighbors(structure=structure, num_neighbors=8)
         self.assertGreater(neigh.centrosymmetry.min(), 0)
         self.assertTrue(
             np.allclose(
-                neigh.centrosymmetry, stk.analyse_centro_symmetry(structure=structure, num_neighbors=8)
+                neigh.centrosymmetry, stk.analyse.get_centro_symmetry_descriptors(structure=structure, num_neighbors=8)
             )
         )
 
     def test_get_all_pairs(self):
         structure = bulk('Fe').repeat(4)
-        neigh = stk.get_neighbors(structure=structure, num_neighbors=8)
+        neigh = stk.analyse.get_neighbors(structure=structure, num_neighbors=8)
         for n in [2, 4, 6]:
             pairs = neigh._get_all_possible_pairs(n)
             self.assertEqual(
diff --git a/tests/test_pyscal.py b/tests/test_pyscal.py
index eadb7e388..3c53f24d0 100644
--- a/tests/test_pyscal.py
+++ b/tests/test_pyscal.py
@@ -4,7 +4,6 @@
 
 import unittest
 import numpy as np
-import pyscal.core as pc
 from ase.build import bulk
 from ase.atoms import Atoms
 import structuretoolkit as stk
@@ -22,14 +21,13 @@ def test_simple_system(self):
         """
         Test a simple ase to pyscal conversion
         """
-        sysp = pc.System()
-        sysp.read_inputfile(self.structure, format="ase")
-        self.assertEqual(len(sysp.atoms), 256)
+        self.assertEqual(len(self.structure), 256)
+        self.assertEqual(len(stk.common.ase_to_pyscal(self.structure).atoms), 256)
 
     def test_steinhardt_parameters_returns(self):
-        self.assertEqual(2, len(stk.get_steinhardt_parameter_structure(self.structure)),
+        self.assertEqual(2, len(stk.analyse.get_steinhardt_parameters(self.structure)),
                          msg='Expected default return value to be a tuple of qs and cluster indices.')
-        self.assertIsInstance(stk.get_steinhardt_parameter_structure(self.structure, n_clusters=None), np.ndarray,
+        self.assertIsInstance(stk.analyse.get_steinhardt_parameters(self.structure, n_clusters=None), np.ndarray,
                               msg='Expected just the qs when no clustering is used.')
 
     def test_steinhardt_parameters_qs(self):
@@ -41,7 +39,7 @@ def test_steinhardt_parameters_qs(self):
 
         qtest = np.random.randint(2, 13, size=2)
 
-        qs, _ = stk.get_steinhardt_parameter_structure(self.structure, cutoff=0, n_clusters=2, q=qtest)
+        qs, _ = stk.analyse.get_steinhardt_parameters(self.structure, cutoff=0, n_clusters=2, q=qtest)
         for c, q in enumerate(qs):
             self.assertLess(np.abs(np.mean(q) - perfect_vals[qtest[c]-2]), 1E-3)
 
@@ -49,27 +47,27 @@ def test_steinhardt_parameters_clustering(self):
         noisy_structure = self.structure.copy()
         noisy_structure.positions += 0.5 * np.random.rand(*noisy_structure.positions.shape)
         n_clusters = 3
-        _, inds = stk.get_steinhardt_parameter_structure(noisy_structure, n_clusters=n_clusters)
+        _, inds = stk.analyse.get_steinhardt_parameters(noisy_structure, n_clusters=n_clusters)
         self.assertEqual(n_clusters, len(np.unique(inds)), msg='Expected to find one label for each cluster.')
 
     def test_centrosymmetry(self):
-        csm = stk.analyse_centro_symmetry(self.structure, num_neighbors=12)
+        csm = stk.analyse.get_centro_symmetry_descriptors(self.structure, num_neighbors=12)
         self.assertLess(np.mean(csm), 1E-5)
 
     def test_cna(self):
-        cna = stk.analyse_cna_adaptive(self.structure)
+        cna = stk.analyse.get_adaptive_cna_descriptors(self.structure)
         self.assertEqual(cna['fcc'], len(self.structure))
 
         rand = np.random.randint(0, len(self.structure))
 
-        cna = stk.analyse_cna_adaptive(self.structure, mode="numeric")
+        cna = stk.analyse.get_adaptive_cna_descriptors(self.structure, mode="numeric")
         self.assertEqual(cna[rand], 1)
 
-        cna = stk.analyse_cna_adaptive(self.structure, mode="str")
+        cna = stk.analyse.get_adaptive_cna_descriptors(self.structure, mode="str")
         self.assertEqual(cna[rand], "fcc")
 
     def test_volume(self):
-        vols = stk.analyse_voronoi_volume(self.structure)
+        vols = stk.analyse.get_voronoi_volumes(self.structure)
         self.assertLess(np.abs(np.mean(vols) - 16.0), 1E-3)
 
 
@@ -89,29 +87,29 @@ def test_steinhardt_parameters(self):
 
         qtest = np.random.randint(2, 13, size=2)
 
-        qs, _ = stk.get_steinhardt_parameter_structure(structure=self.al_fcc_4, cutoff=0, n_clusters=2, q=qtest)
+        qs, _ = stk.analyse.get_steinhardt_parameters(structure=self.al_fcc_4, cutoff=0, n_clusters=2, q=qtest)
         for c, q in enumerate(qs):
             self.assertLess(np.abs(np.mean(q) - perfect_vals[qtest[c]-2]), 1E-3)
 
     def test_analyse_pyscal_centro_symmetry(self):
-        self.assertTrue(all([np.isclose(v, 0.0) for v in stk.analyse_centro_symmetry(structure=self.al_fcc, num_neighbors=12)]))
-        self.assertTrue(all([np.isclose(v, 0.0) for v in stk.analyse_centro_symmetry(structure=self.fe_bcc, num_neighbors=8)]))
-        self.assertTrue(all([np.isclose(v, 8.7025) for v in stk.analyse_centro_symmetry(structure=self.ti_hcp, num_neighbors=12)]))
-        self.assertTrue(all([np.isclose(v, 14.742449) for v in stk.analyse_centro_symmetry(structure=self.si_dia, num_neighbors=4)]))
-        self.assertEqual(len(stk.analyse_centro_symmetry(structure=self.al_fcc)), len(self.al_fcc))
-        self.assertEqual(len(stk.analyse_centro_symmetry(structure=self.fe_bcc)), len(self.fe_bcc))
-        self.assertEqual(len(stk.analyse_centro_symmetry(structure=self.ti_hcp)), len(self.ti_hcp))
-        self.assertEqual(len(stk.analyse_centro_symmetry(structure=self.si_dia)), len(self.si_dia))
+        self.assertTrue(all([np.isclose(v, 0.0) for v in stk.analyse.get_centro_symmetry_descriptors(structure=self.al_fcc, num_neighbors=12)]))
+        self.assertTrue(all([np.isclose(v, 0.0) for v in stk.analyse.get_centro_symmetry_descriptors(structure=self.fe_bcc, num_neighbors=8)]))
+        self.assertTrue(all([np.isclose(v, 8.7025) for v in stk.analyse.get_centro_symmetry_descriptors(structure=self.ti_hcp, num_neighbors=12)]))
+        self.assertTrue(all([np.isclose(v, 14.742449) for v in stk.analyse.get_centro_symmetry_descriptors(structure=self.si_dia, num_neighbors=4)]))
+        self.assertEqual(len(stk.analyse.get_centro_symmetry_descriptors(structure=self.al_fcc)), len(self.al_fcc))
+        self.assertEqual(len(stk.analyse.get_centro_symmetry_descriptors(structure=self.fe_bcc)), len(self.fe_bcc))
+        self.assertEqual(len(stk.analyse.get_centro_symmetry_descriptors(structure=self.ti_hcp)), len(self.ti_hcp))
+        self.assertEqual(len(stk.analyse.get_centro_symmetry_descriptors(structure=self.si_dia)), len(self.si_dia))
 
     def test_analyse_pyscal_voronoi_volume(self):
-        self.assertAlmostEqual(np.mean(stk.analyse_voronoi_volume(structure=self.al_fcc)), 16.60753125)
-        self.assertAlmostEqual(np.mean(stk.analyse_voronoi_volume(structure=self.fe_bcc)), 11.8199515)
-        self.assertAlmostEqual(np.mean(stk.analyse_voronoi_volume(structure=self.ti_hcp)), 17.65294557)
-        self.assertAlmostEqual(np.mean(stk.analyse_voronoi_volume(structure=self.si_dia)), 20.01287587)
-        self.assertEqual(len(stk.analyse_voronoi_volume(structure=self.al_fcc)), len(self.al_fcc))
-        self.assertEqual(len(stk.analyse_voronoi_volume(structure=self.fe_bcc)), len(self.fe_bcc))
-        self.assertEqual(len(stk.analyse_voronoi_volume(structure=self.ti_hcp)), len(self.ti_hcp))
-        self.assertEqual(len(stk.analyse_voronoi_volume(structure=self.si_dia)), len(self.si_dia))
+        self.assertAlmostEqual(np.mean(stk.analyse.get_voronoi_volumes(structure=self.al_fcc)), 16.60753125)
+        self.assertAlmostEqual(np.mean(stk.analyse.get_voronoi_volumes(structure=self.fe_bcc)), 11.8199515)
+        self.assertAlmostEqual(np.mean(stk.analyse.get_voronoi_volumes(structure=self.ti_hcp)), 17.65294557)
+        self.assertAlmostEqual(np.mean(stk.analyse.get_voronoi_volumes(structure=self.si_dia)), 20.01287587)
+        self.assertEqual(len(stk.analyse.get_voronoi_volumes(structure=self.al_fcc)), len(self.al_fcc))
+        self.assertEqual(len(stk.analyse.get_voronoi_volumes(structure=self.fe_bcc)), len(self.fe_bcc))
+        self.assertEqual(len(stk.analyse.get_voronoi_volumes(structure=self.ti_hcp)), len(self.ti_hcp))
+        self.assertEqual(len(stk.analyse.get_voronoi_volumes(structure=self.si_dia)), len(self.si_dia))
 
     def test_analyse_pyscal_cna_adaptive(self):
         pyscal_keys = [
@@ -124,81 +122,81 @@ def test_analyse_pyscal_cna_adaptive(self):
             'CommonNeighborAnalysis.counts.BCC',
              'CommonNeighborAnalysis.counts.ICO'
         ]
-        res_dict_total = stk.analyse_cna_adaptive(structure=self.al_fcc, mode="total", ovito_compatibility=False)
+        res_dict_total = stk.analyse.get_adaptive_cna_descriptors(structure=self.al_fcc, mode="total", ovito_compatibility=False)
         self.assertEqual(sum([k in res_dict_total.keys() for k in pyscal_keys]), len(pyscal_keys))
         self.assertEqual(res_dict_total[pyscal_keys[1]], len(self.al_fcc))
-        res_dict_total = stk.analyse_cna_adaptive(structure=self.fe_bcc, mode="total", ovito_compatibility=False)
+        res_dict_total = stk.analyse.get_adaptive_cna_descriptors(structure=self.fe_bcc, mode="total", ovito_compatibility=False)
         self.assertEqual(sum([k in res_dict_total.keys() for k in pyscal_keys]), len(pyscal_keys))
         self.assertEqual(res_dict_total[pyscal_keys[3]], len(self.fe_bcc))
-        res_dict_total = stk.analyse_cna_adaptive(structure=self.ti_hcp, mode="total", ovito_compatibility=False)
+        res_dict_total = stk.analyse.get_adaptive_cna_descriptors(structure=self.ti_hcp, mode="total", ovito_compatibility=False)
         self.assertEqual(sum([k in res_dict_total.keys() for k in pyscal_keys]), len(pyscal_keys))
         self.assertEqual(res_dict_total[pyscal_keys[2]], len(self.ti_hcp))
-        res_dict_total = stk.analyse_cna_adaptive(structure=self.si_dia, mode="total", ovito_compatibility=False)
+        res_dict_total = stk.analyse.get_adaptive_cna_descriptors(structure=self.si_dia, mode="total", ovito_compatibility=False)
         self.assertEqual(sum([k in res_dict_total.keys() for k in pyscal_keys]), len(pyscal_keys))
         self.assertEqual(res_dict_total[pyscal_keys[0]], len(self.si_dia))
 
-        res_numeric = stk.analyse_cna_adaptive(structure=self.al_fcc, mode="numeric", ovito_compatibility=False)
+        res_numeric = stk.analyse.get_adaptive_cna_descriptors(structure=self.al_fcc, mode="numeric", ovito_compatibility=False)
         self.assertEqual(len(res_numeric), len(self.al_fcc))
         self.assertTrue(all([v == 1 for v in res_numeric]))
-        res_numeric = stk.analyse_cna_adaptive(structure=self.fe_bcc, mode="numeric", ovito_compatibility=False)
+        res_numeric = stk.analyse.get_adaptive_cna_descriptors(structure=self.fe_bcc, mode="numeric", ovito_compatibility=False)
         self.assertEqual(len(res_numeric), len(self.fe_bcc))
         self.assertTrue(all([v == 3 for v in res_numeric]))
-        res_numeric = stk.analyse_cna_adaptive(structure=self.ti_hcp, mode="numeric", ovito_compatibility=False)
+        res_numeric = stk.analyse.get_adaptive_cna_descriptors(structure=self.ti_hcp, mode="numeric", ovito_compatibility=False)
         self.assertEqual(len(res_numeric), len(self.ti_hcp))
         self.assertTrue(all([v == 2 for v in res_numeric]))
-        res_numeric = stk.analyse_cna_adaptive(structure=self.si_dia, mode="numeric", ovito_compatibility=False)
+        res_numeric = stk.analyse.get_adaptive_cna_descriptors(structure=self.si_dia, mode="numeric", ovito_compatibility=False)
         self.assertEqual(len(res_numeric), len(self.si_dia))
         self.assertTrue(all([v == 0 for v in res_numeric]))
 
-        res_str = stk.analyse_cna_adaptive(structure=self.al_fcc, mode="str", ovito_compatibility=False)
+        res_str = stk.analyse.get_adaptive_cna_descriptors(structure=self.al_fcc, mode="str", ovito_compatibility=False)
         self.assertEqual(len(res_str), len(self.al_fcc))
         self.assertTrue(all([v == 'fcc' for v in res_str]))
-        res_str = stk.analyse_cna_adaptive(structure=self.fe_bcc, mode="str", ovito_compatibility=False)
+        res_str = stk.analyse.get_adaptive_cna_descriptors(structure=self.fe_bcc, mode="str", ovito_compatibility=False)
         self.assertEqual(len(res_str), len(self.fe_bcc))
         self.assertTrue(all([v == 'bcc' for v in res_str]))
-        res_str = stk.analyse_cna_adaptive(structure=self.ti_hcp, mode="str", ovito_compatibility=False)
+        res_str = stk.analyse.get_adaptive_cna_descriptors(structure=self.ti_hcp, mode="str", ovito_compatibility=False)
         self.assertEqual(len(res_str), len(self.ti_hcp))
         self.assertTrue(all([v == 'hcp' for v in res_str]))
-        res_str = stk.analyse_cna_adaptive(structure=self.si_dia, mode="str", ovito_compatibility=False)
+        res_str = stk.analyse.get_adaptive_cna_descriptors(structure=self.si_dia, mode="str", ovito_compatibility=False)
         self.assertEqual(len(res_str), len(self.si_dia))
         self.assertTrue(all([v == 'others' for v in res_str]))
 
-        res_dict_total = stk.analyse_cna_adaptive(structure=self.al_fcc, mode="total", ovito_compatibility=True)
+        res_dict_total = stk.analyse.get_adaptive_cna_descriptors(structure=self.al_fcc, mode="total", ovito_compatibility=True)
         self.assertEqual(sum([k in res_dict_total.keys() for k in ovito_keys]), len(ovito_keys))
         self.assertEqual(res_dict_total[ovito_keys[1]], len(self.al_fcc))
-        res_dict_total = stk.analyse_cna_adaptive(structure=self.fe_bcc, mode="total", ovito_compatibility=True)
+        res_dict_total = stk.analyse.get_adaptive_cna_descriptors(structure=self.fe_bcc, mode="total", ovito_compatibility=True)
         self.assertEqual(sum([k in res_dict_total.keys() for k in ovito_keys]), len(ovito_keys))
         self.assertEqual(res_dict_total[ovito_keys[3]], len(self.fe_bcc))
-        res_dict_total = stk.analyse_cna_adaptive(structure=self.ti_hcp, mode="total", ovito_compatibility=True)
+        res_dict_total = stk.analyse.get_adaptive_cna_descriptors(structure=self.ti_hcp, mode="total", ovito_compatibility=True)
         self.assertEqual(sum([k in res_dict_total.keys() for k in ovito_keys]), len(ovito_keys))
         self.assertEqual(res_dict_total[ovito_keys[2]], len(self.ti_hcp))
-        res_dict_total = stk.analyse_cna_adaptive(structure=self.si_dia, mode="total", ovito_compatibility=True)
+        res_dict_total = stk.analyse.get_adaptive_cna_descriptors(structure=self.si_dia, mode="total", ovito_compatibility=True)
         self.assertEqual(sum([k in res_dict_total.keys() for k in ovito_keys]), len(ovito_keys))
         self.assertEqual(res_dict_total[ovito_keys[0]], len(self.si_dia))
 
-        res_numeric = stk.analyse_cna_adaptive(structure=self.al_fcc, mode="numeric", ovito_compatibility=True)
+        res_numeric = stk.analyse.get_adaptive_cna_descriptors(structure=self.al_fcc, mode="numeric", ovito_compatibility=True)
         self.assertEqual(len(res_numeric), len(self.al_fcc))
         self.assertTrue(all([v == 1 for v in res_numeric]))
-        res_numeric = stk.analyse_cna_adaptive(structure=self.fe_bcc, mode="numeric", ovito_compatibility=True)
+        res_numeric = stk.analyse.get_adaptive_cna_descriptors(structure=self.fe_bcc, mode="numeric", ovito_compatibility=True)
         self.assertEqual(len(res_numeric), len(self.fe_bcc))
         self.assertTrue(all([v == 3 for v in res_numeric]))
-        res_numeric = stk.analyse_cna_adaptive(structure=self.ti_hcp, mode="numeric", ovito_compatibility=True)
+        res_numeric = stk.analyse.get_adaptive_cna_descriptors(structure=self.ti_hcp, mode="numeric", ovito_compatibility=True)
         self.assertEqual(len(res_numeric), len(self.ti_hcp))
         self.assertTrue(all([v == 2 for v in res_numeric]))
-        res_numeric = stk.analyse_cna_adaptive(structure=self.si_dia, mode="numeric", ovito_compatibility=True)
+        res_numeric = stk.analyse.get_adaptive_cna_descriptors(structure=self.si_dia, mode="numeric", ovito_compatibility=True)
         self.assertEqual(len(res_numeric), len(self.si_dia))
         self.assertTrue(all([v == 0 for v in res_numeric]))
 
-        res_str = stk.analyse_cna_adaptive(structure=self.al_fcc, mode="str", ovito_compatibility=True)
+        res_str = stk.analyse.get_adaptive_cna_descriptors(structure=self.al_fcc, mode="str", ovito_compatibility=True)
         self.assertEqual(len(res_str), len(self.al_fcc))
         self.assertTrue(all([v == 'FCC' for v in res_str]))
-        res_str = stk.analyse_cna_adaptive(structure=self.fe_bcc, mode="str", ovito_compatibility=True)
+        res_str = stk.analyse.get_adaptive_cna_descriptors(structure=self.fe_bcc, mode="str", ovito_compatibility=True)
         self.assertEqual(len(res_str), len(self.fe_bcc))
         self.assertTrue(all([v == 'BCC' for v in res_str]))
-        res_str = stk.analyse_cna_adaptive(structure=self.ti_hcp, mode="str", ovito_compatibility=True)
+        res_str = stk.analyse.get_adaptive_cna_descriptors(structure=self.ti_hcp, mode="str", ovito_compatibility=True)
         self.assertEqual(len(res_str), len(self.ti_hcp))
         self.assertTrue(all([v == 'HCP' for v in res_str]))
-        res_str = stk.analyse_cna_adaptive(structure=self.si_dia, mode="str", ovito_compatibility=True)
+        res_str = stk.analyse.get_adaptive_cna_descriptors(structure=self.si_dia, mode="str", ovito_compatibility=True)
         self.assertEqual(len(res_str), len(self.si_dia))
         self.assertTrue(all([v == 'Other' for v in res_str]))
 
@@ -217,81 +215,81 @@ def test_analyse_pyscal_diamond_structure(self):
             'IdentifyDiamond.counts.HEX_DIAMOND_SECOND_NEIGHBOR',
             'IdentifyDiamond.counts.OTHER'
         ]
-        res_dict_total = stk.analyse_diamond_structure(structure=self.al_fcc, mode="total", ovito_compatibility=False)
+        res_dict_total = stk.analyse.get_diamond_structure_descriptors(structure=self.al_fcc, mode="total", ovito_compatibility=False)
         self.assertEqual(sum([k in res_dict_total.keys() for k in pyscal_keys]), len(pyscal_keys))
         self.assertEqual(res_dict_total[pyscal_keys[0]], len(self.al_fcc))
-        res_dict_total = stk.analyse_diamond_structure(structure=self.fe_bcc, mode="total", ovito_compatibility=False)
+        res_dict_total = stk.analyse.get_diamond_structure_descriptors(structure=self.fe_bcc, mode="total", ovito_compatibility=False)
         self.assertEqual(sum([k in res_dict_total.keys() for k in pyscal_keys]), len(pyscal_keys))
         self.assertEqual(res_dict_total[pyscal_keys[0]], len(self.fe_bcc))
-        res_dict_total = stk.analyse_diamond_structure(structure=self.ti_hcp, mode="total", ovito_compatibility=False)
+        res_dict_total = stk.analyse.get_diamond_structure_descriptors(structure=self.ti_hcp, mode="total", ovito_compatibility=False)
         self.assertEqual(sum([k in res_dict_total.keys() for k in pyscal_keys]), len(pyscal_keys))
         self.assertEqual(res_dict_total[pyscal_keys[0]], len(self.ti_hcp))
-        res_dict_total = stk.analyse_diamond_structure(structure=self.si_dia, mode="total", ovito_compatibility=False)
+        res_dict_total = stk.analyse.get_diamond_structure_descriptors(structure=self.si_dia, mode="total", ovito_compatibility=False)
         self.assertEqual(sum([k in res_dict_total.keys() for k in pyscal_keys]), len(pyscal_keys))
         self.assertEqual(res_dict_total[pyscal_keys[5]], len(self.si_dia))
 
-        res_numeric = stk.analyse_diamond_structure(structure=self.al_fcc, mode="numeric", ovito_compatibility=False)
+        res_numeric = stk.analyse.get_diamond_structure_descriptors(structure=self.al_fcc, mode="numeric", ovito_compatibility=False)
         self.assertEqual(len(res_numeric), len(self.al_fcc))
         self.assertTrue(all([v == 0 for v in res_numeric]))
-        res_numeric = stk.analyse_diamond_structure(structure=self.fe_bcc, mode="numeric", ovito_compatibility=False)
+        res_numeric = stk.analyse.get_diamond_structure_descriptors(structure=self.fe_bcc, mode="numeric", ovito_compatibility=False)
         self.assertEqual(len(res_numeric), len(self.fe_bcc))
         self.assertTrue(all([v == 0 for v in res_numeric]))
-        res_numeric = stk.analyse_diamond_structure(structure=self.ti_hcp, mode="numeric", ovito_compatibility=False)
+        res_numeric = stk.analyse.get_diamond_structure_descriptors(structure=self.ti_hcp, mode="numeric", ovito_compatibility=False)
         self.assertEqual(len(res_numeric), len(self.ti_hcp))
         self.assertTrue(all([v == 0 for v in res_numeric]))
-        res_numeric = stk.analyse_diamond_structure(structure=self.si_dia, mode="numeric", ovito_compatibility=False)
+        res_numeric = stk.analyse.get_diamond_structure_descriptors(structure=self.si_dia, mode="numeric", ovito_compatibility=False)
         self.assertEqual(len(res_numeric), len(self.si_dia))
         self.assertTrue(all([v == 5 for v in res_numeric]))
 
-        res_str = stk.analyse_diamond_structure(structure=self.al_fcc, mode="str", ovito_compatibility=False)
+        res_str = stk.analyse.get_diamond_structure_descriptors(structure=self.al_fcc, mode="str", ovito_compatibility=False)
         self.assertEqual(len(res_str), len(self.al_fcc))
         self.assertTrue(all([v == 'others' for v in res_str]))
-        res_str = stk.analyse_diamond_structure(structure=self.fe_bcc, mode="str", ovito_compatibility=False)
+        res_str = stk.analyse.get_diamond_structure_descriptors(structure=self.fe_bcc, mode="str", ovito_compatibility=False)
         self.assertEqual(len(res_str), len(self.fe_bcc))
         self.assertTrue(all([v == 'others' for v in res_str]))
-        res_str = stk.analyse_diamond_structure(structure=self.ti_hcp, mode="str", ovito_compatibility=False)
+        res_str = stk.analyse.get_diamond_structure_descriptors(structure=self.ti_hcp, mode="str", ovito_compatibility=False)
         self.assertEqual(len(res_str), len(self.ti_hcp))
         self.assertTrue(all([v == 'others' for v in res_str]))
-        res_str = stk.analyse_diamond_structure(structure=self.si_dia, mode="str", ovito_compatibility=False)
+        res_str = stk.analyse.get_diamond_structure_descriptors(structure=self.si_dia, mode="str", ovito_compatibility=False)
         self.assertEqual(len(res_str), len(self.si_dia))
         self.assertTrue(all([v == 'cubic diamond' for v in res_str]))
 
-        res_dict_total = stk.analyse_diamond_structure(structure=self.al_fcc, mode="total", ovito_compatibility=True)
+        res_dict_total = stk.analyse.get_diamond_structure_descriptors(structure=self.al_fcc, mode="total", ovito_compatibility=True)
         self.assertEqual(sum([k in res_dict_total.keys() for k in ovito_keys]), len(ovito_keys))
         self.assertEqual(res_dict_total[ovito_keys[6]], len(self.al_fcc))
-        res_dict_total = stk.analyse_diamond_structure(structure=self.fe_bcc, mode="total", ovito_compatibility=True)
+        res_dict_total = stk.analyse.get_diamond_structure_descriptors(structure=self.fe_bcc, mode="total", ovito_compatibility=True)
         self.assertEqual(sum([k in res_dict_total.keys() for k in ovito_keys]), len(ovito_keys))
         self.assertEqual(res_dict_total[ovito_keys[6]], len(self.fe_bcc))
-        res_dict_total = stk.analyse_diamond_structure(structure=self.ti_hcp, mode="total", ovito_compatibility=True)
+        res_dict_total = stk.analyse.get_diamond_structure_descriptors(structure=self.ti_hcp, mode="total", ovito_compatibility=True)
         self.assertEqual(sum([k in res_dict_total.keys() for k in ovito_keys]), len(ovito_keys))
         self.assertEqual(res_dict_total[ovito_keys[6]], len(self.ti_hcp))
-        res_dict_total = stk.analyse_diamond_structure(structure=self.si_dia, mode="total", ovito_compatibility=True)
+        res_dict_total = stk.analyse.get_diamond_structure_descriptors(structure=self.si_dia, mode="total", ovito_compatibility=True)
         self.assertEqual(sum([k in res_dict_total.keys() for k in ovito_keys]), len(ovito_keys))
         self.assertEqual(res_dict_total[ovito_keys[0]], len(self.si_dia))
 
-        res_numeric = stk.analyse_diamond_structure(structure=self.al_fcc, mode="numeric", ovito_compatibility=True)
+        res_numeric = stk.analyse.get_diamond_structure_descriptors(structure=self.al_fcc, mode="numeric", ovito_compatibility=True)
         self.assertEqual(len(res_numeric), len(self.al_fcc))
         self.assertTrue(all([v == 6 for v in res_numeric]))
-        res_numeric = stk.analyse_diamond_structure(structure=self.fe_bcc, mode="numeric", ovito_compatibility=True)
+        res_numeric = stk.analyse.get_diamond_structure_descriptors(structure=self.fe_bcc, mode="numeric", ovito_compatibility=True)
         self.assertEqual(len(res_numeric), len(self.fe_bcc))
         self.assertTrue(all([v == 6 for v in res_numeric]))
-        res_numeric = stk.analyse_diamond_structure(structure=self.ti_hcp, mode="numeric", ovito_compatibility=True)
+        res_numeric = stk.analyse.get_diamond_structure_descriptors(structure=self.ti_hcp, mode="numeric", ovito_compatibility=True)
         self.assertEqual(len(res_numeric), len(self.ti_hcp))
         self.assertTrue(all([v == 6 for v in res_numeric]))
-        res_numeric = stk.analyse_diamond_structure(structure=self.si_dia, mode="numeric", ovito_compatibility=True)
+        res_numeric = stk.analyse.get_diamond_structure_descriptors(structure=self.si_dia, mode="numeric", ovito_compatibility=True)
         self.assertEqual(len(res_numeric), len(self.si_dia))
         self.assertTrue(all([v == 0 for v in res_numeric]))
 
-        res_str = stk.analyse_diamond_structure(structure=self.al_fcc, mode="str", ovito_compatibility=True)
+        res_str = stk.analyse.get_diamond_structure_descriptors(structure=self.al_fcc, mode="str", ovito_compatibility=True)
         self.assertEqual(len(res_str), len(self.al_fcc))
         self.assertTrue(all([v == 'Other' for v in res_str]))
-        res_str = stk.analyse_diamond_structure(structure=self.fe_bcc, mode="str", ovito_compatibility=True)
+        res_str = stk.analyse.get_diamond_structure_descriptors(structure=self.fe_bcc, mode="str", ovito_compatibility=True)
         self.assertEqual(len(res_str), len(self.fe_bcc))
         self.assertTrue(all([v == 'Other' for v in res_str]))
-        res_str = stk.analyse_diamond_structure(structure=self.ti_hcp, mode="str", ovito_compatibility=True)
+        res_str = stk.analyse.get_diamond_structure_descriptors(structure=self.ti_hcp, mode="str", ovito_compatibility=True)
         self.assertEqual(len(res_str), len(self.ti_hcp))
         self.assertTrue(all([v == 'Other' for v in res_str]))
-        res_str = stk.analyse_diamond_structure(structure=self.si_dia, mode="str", ovito_compatibility=True)
+        res_str = stk.analyse.get_diamond_structure_descriptors(structure=self.si_dia, mode="str", ovito_compatibility=True)
         self.assertEqual(len(res_str), len(self.si_dia))
         self.assertTrue(all([v == 'Cubic diamond' for v in res_str]))
 
diff --git a/tests/test_strain.py b/tests/test_strain.py
index e21d6b917..4a6cc4a18 100644
--- a/tests/test_strain.py
+++ b/tests/test_strain.py
@@ -9,12 +9,12 @@ class TestAtoms(unittest.TestCase):
     @classmethod
     def setUpClass(cls):
         bulk_structure = bulk('Fe', cubic=True)
-        cls.strain = stk.get_strain(structure=bulk_structure, ref_structure=bulk_structure, return_object=True)
+        cls.strain = stk.analyse.get_strain(structure=bulk_structure, ref_structure=bulk_structure, return_object=True)
 
     def test_number_of_neighbors(self):
         self.assertEqual(self.strain.num_neighbors, 8)
         bulk_structure = bulk('Al', cubic=True)
-        strain = stk.get_strain(structure=bulk_structure, ref_structure=bulk_structure, return_object=True)
+        strain = stk.analyse.get_strain(structure=bulk_structure, ref_structure=bulk_structure, return_object=True)
         self.assertEqual(strain.num_neighbors, 12)
 
     def test_get_angle(self):
diff --git a/tests/test_symmetry.py b/tests/test_symmetry.py
index 6c77b0f1a..defd8d326 100644
--- a/tests/test_symmetry.py
+++ b/tests/test_symmetry.py
@@ -6,7 +6,6 @@
 import numpy as np
 from ase.build import bulk
 from ase.atoms import Atoms
-from structuretoolkit.analyse.symmetry import SymmetryError
 import structuretoolkit as stk
 
 
@@ -14,22 +13,22 @@ class TestAtoms(unittest.TestCase):
     def test_get_arg_equivalent_sites(self):
         a_0 = 4.0
         structure = bulk('Al', cubic=True, a=a_0).repeat(2)
-        sites = stk.get_wrapped_coordinates(structure=structure, positions=structure.positions + np.array([0, 0, 0.5 * a_0]))
+        sites = stk.common.get_wrapped_coordinates(structure=structure, positions=structure.positions + np.array([0, 0, 0.5 * a_0]))
         v_position = structure.positions[0]
         del structure[0]
         pairs = np.stack((
-            stk.get_symmetry(structure=structure).get_arg_equivalent_sites(sites),
-            np.unique(np.round(stk.get_distances_array(structure=structure, p1=v_position, p2=sites), decimals=2), return_inverse=True)[1]
+            stk.analyse.get_symmetry(structure=structure).get_arg_equivalent_sites(sites),
+            np.unique(np.round(stk.analyse.get_distances_array(structure=structure, p1=v_position, p2=sites), decimals=2), return_inverse=True)[1]
         ), axis=-1)
         unique_pairs = np.unique(pairs, axis=0)
         self.assertEqual(len(unique_pairs), len(np.unique(unique_pairs[:, 0])))
         with self.assertRaises(ValueError):
-            stk.get_symmetry(structure=structure).get_arg_equivalent_sites([0, 0, 0])
+            stk.analyse.get_symmetry(structure=structure).get_arg_equivalent_sites([0, 0, 0])
 
     def test_generate_equivalent_points(self):
         a_0 = 4
         structure = bulk('Al', cubic=True, a=a_0)
-        sym = stk.get_symmetry(structure)
+        sym = stk.analyse.get_symmetry(structure)
         self.assertEqual(
             len(structure),
             len(sym.generate_equivalent_points([0, 0, 0.5 * a_0]))
@@ -40,7 +39,7 @@ def test_generate_equivalent_points(self):
         y = np.tile(y, len(sym_x))
         sym_x = sym_x.reshape(-1, 3)
         xy = np.round(
-            [stk.get_neighborhood(structure, sym_x, num_neighbors=1).distances.flatten(), y],
+            [stk.analyse.get_neighborhood(structure, sym_x, num_neighbors=1).distances.flatten(), y],
             decimals=8
         )
         self.assertEqual(
@@ -52,55 +51,55 @@ def test_generate_equivalent_points(self):
     def test_get_symmetry(self):
         cell = 2.2 * np.identity(3)
         Al = Atoms("AlAl", positions=[(0, 0, 0), (0.5, 0.5, 0.5)], cell=cell, pbc=True).repeat(2)
-        self.assertEqual(len(set(stk.get_symmetry(structure=Al)["equivalent_atoms"])), 1)
-        self.assertEqual(len(stk.get_symmetry(structure=Al)["translations"]), 96)
+        self.assertEqual(len(set(stk.analyse.get_symmetry(structure=Al)["equivalent_atoms"])), 1)
+        self.assertEqual(len(stk.analyse.get_symmetry(structure=Al)["translations"]), 96)
         self.assertEqual(
-            len(stk.get_symmetry(structure=Al)["translations"]), len(stk.get_symmetry(structure=Al)["rotations"])
+            len(stk.analyse.get_symmetry(structure=Al)["translations"]), len(stk.analyse.get_symmetry(structure=Al)["rotations"])
         )
         cell = 2.2 * np.identity(3)
         Al = Atoms("AlAl", scaled_positions=[(0, 0, 0), (0.5, 0.5, 0.5)], cell=cell, pbc=True)
         v = np.random.rand(6).reshape(-1, 3)
-        self.assertAlmostEqual(np.linalg.norm(stk.get_symmetry(structure=Al).symmetrize_vectors(v)), 0)
+        self.assertAlmostEqual(np.linalg.norm(stk.analyse.get_symmetry(structure=Al).symmetrize_vectors(v)), 0)
         vv = np.random.rand(12).reshape(2, 2, 3)
-        for vvv in stk.get_symmetry(structure=Al).symmetrize_vectors(vv):
+        for vvv in stk.analyse.get_symmetry(structure=Al).symmetrize_vectors(vv):
             self.assertAlmostEqual(np.linalg.norm(vvv), 0)
         Al.positions[0, 0] += 0.01
-        w = stk.get_symmetry(structure=Al).symmetrize_vectors(v)
+        w = stk.analyse.get_symmetry(structure=Al).symmetrize_vectors(v)
         self.assertAlmostEqual(np.absolute(w[:, 0]).sum(), np.linalg.norm(w, axis=-1).sum())
 
     def test_get_symmetry_dataset(self):
         cell = 2.2 * np.identity(3)
         Al_sc = Atoms("AlAl", scaled_positions=[(0, 0, 0), (0.5, 0.5, 0.5)], cell=cell)
         Al_sc = Al_sc.repeat([2, 2, 2])
-        self.assertEqual(stk.get_symmetry(structure=Al_sc).info["number"], 229)
+        self.assertEqual(stk.analyse.get_symmetry(structure=Al_sc).info["number"], 229)
 
     def test_get_ir_reciprocal_mesh(self):
         cell = 2.2 * np.identity(3)
         Al_sc = Atoms("AlAl", scaled_positions=[(0, 0, 0), (0.5, 0.5, 0.5)], cell=cell)
-        self.assertEqual(len(stk.get_symmetry(structure=Al_sc).get_ir_reciprocal_mesh([3, 3, 3])[0]), 27)
+        self.assertEqual(len(stk.analyse.get_symmetry(structure=Al_sc).get_ir_reciprocal_mesh([3, 3, 3])[0]), 27)
 
     def test_get_primitive_cell(self):
         cell = 2.2 * np.identity(3)
         basis = Atoms("AlFe", scaled_positions=[(0, 0, 0), (0.5, 0.5, 0.5)], cell=cell)
         structure = basis.repeat([2, 2, 2])
-        sym = stk.get_symmetry(structure=structure)
+        sym = stk.analyse.get_symmetry(structure=structure)
         self.assertEqual(len(basis), len(sym.get_primitive_cell(standardize=True)))
-        self.assertEqual(stk.get_symmetry(structure=sym.get_primitive_cell()).spacegroup["Number"], 221)
+        self.assertEqual(stk.analyse.get_symmetry(structure=sym.get_primitive_cell()).spacegroup["Number"], 221)
 
     def test_get_equivalent_points(self):
         basis = Atoms("FeFe", positions=[[0.01, 0, 0], [0.5, 0.5, 0.5]], cell=np.identity(3))
-        arr = stk.get_symmetry(structure=basis).generate_equivalent_points([0, 0, 0.5])
+        arr = stk.analyse.get_symmetry(structure=basis).generate_equivalent_points([0, 0, 0.5])
         self.assertAlmostEqual(np.linalg.norm(arr - np.array([0.51, 0.5, 0]), axis=-1).min(), 0)
 
     def test_get_space_group(self):
         cell = 2.2 * np.identity(3)
         Al_sc = Atoms("AlAl", scaled_positions=[(0, 0, 0), (0.5, 0.5, 0.5)], cell=cell)
-        self.assertEqual(stk.get_symmetry(structure=Al_sc).spacegroup["InternationalTableSymbol"], "Im-3m")
-        self.assertEqual(stk.get_symmetry(structure=Al_sc).spacegroup["Number"], 229)
+        self.assertEqual(stk.analyse.get_symmetry(structure=Al_sc).spacegroup["InternationalTableSymbol"], "Im-3m")
+        self.assertEqual(stk.analyse.get_symmetry(structure=Al_sc).spacegroup["Number"], 229)
         cell = 4.2 * (0.5 * np.ones((3, 3)) - 0.5 * np.eye(3))
         Al_fcc = Atoms("Al", scaled_positions=[(0, 0, 0)], cell=cell)
-        self.assertEqual(stk.get_symmetry(structure=Al_fcc).spacegroup["InternationalTableSymbol"], "Fm-3m")
-        self.assertEqual(stk.get_symmetry(structure=Al_fcc).spacegroup["Number"], 225)
+        self.assertEqual(stk.analyse.get_symmetry(structure=Al_fcc).spacegroup["InternationalTableSymbol"], "Fm-3m")
+        self.assertEqual(stk.analyse.get_symmetry(structure=Al_fcc).spacegroup["Number"], 225)
         a = 3.18
         c = 1.623 * a
         cell = np.eye(3)
@@ -110,7 +109,7 @@ def test_get_space_group(self):
         cell[1, 1] = np.sqrt(3) * a / 2.0
         pos = np.array([[0.0, 0.0, 0.0], [1.0 / 3.0, 2.0 / 3.0, 1.0 / 2.0]])
         Mg_hcp = Atoms("Mg2", scaled_positions=pos, cell=cell)
-        self.assertEqual(stk.get_symmetry(structure=Mg_hcp).spacegroup["Number"], 194)
+        self.assertEqual(stk.analyse.get_symmetry(structure=Mg_hcp).spacegroup["Number"], 194)
         cell = np.eye(3)
         cell[0, 0] = a
         cell[2, 2] = c
@@ -124,16 +123,16 @@ def test_get_space_group(self):
             ]
         )
         Mg_hcp = Atoms("Mg4", scaled_positions=pos, cell=cell)
-        self.assertEqual(stk.get_symmetry(structure=Mg_hcp).spacegroup["Number"], 194)
+        self.assertEqual(stk.analyse.get_symmetry(structure=Mg_hcp).spacegroup["Number"], 194)
 
     def test_permutations(self):
         structure = bulk('Al', cubic=True).repeat(2)
         x_vacancy = structure.positions[0]
         del structure[0]
-        neigh = stk.get_neighborhood(structure=structure, positions=x_vacancy)
+        neigh = stk.analyse.get_neighborhood(structure=structure, positions=x_vacancy)
         vec = np.zeros_like(structure.positions)
         vec[neigh.indices[0]] = neigh.vecs[0]
-        sym = stk.get_symmetry(structure=structure)
+        sym = stk.analyse.get_symmetry(structure=structure)
         all_vectors = np.einsum('ijk,ink->inj', sym.rotations, vec[sym.permutations])
         for i, v in zip(neigh.indices, neigh.vecs):
             vec = np.zeros_like(structure.positions)
@@ -142,11 +141,11 @@ def test_permutations(self):
 
     def test_arg_equivalent_vectors(self):
         structure = bulk('Al', cubic=True).repeat(2)
-        self.assertEqual(np.unique(stk.get_symmetry(structure=structure).arg_equivalent_vectors).squeeze(), 0)
+        self.assertEqual(np.unique(stk.analyse.get_symmetry(structure=structure).arg_equivalent_vectors).squeeze(), 0)
         x_v = structure.positions[0]
         del structure[0]
-        arg_v = stk.get_symmetry(structure=structure).arg_equivalent_vectors
-        dx = stk.get_distances_array(structure=structure, p1=structure.positions, p2=x_v, vectors=True)
+        arg_v = stk.analyse.get_symmetry(structure=structure).arg_equivalent_vectors
+        dx = stk.analyse.get_distances_array(structure=structure, p1=structure.positions, p2=x_v, vectors=True)
         dx_round = np.round(np.absolute(dx), decimals=3)
         self.assertEqual(len(np.unique(dx_round + arg_v)), len(np.unique(arg_v)))
 
@@ -155,8 +154,8 @@ def test_error(self):
 
         structure = bulk('Al')
         structure += structure[-1]
-        with self.assertRaises(SymmetryError):
-            stk.get_symmetry(structure=structure)
+        with self.assertRaises(stk.common.SymmetryError):
+            stk.analyse.get_symmetry(structure=structure)
 
 
 if __name__ == "__main__":