diff --git a/docs/applications/cylinder_configurations.ipynb b/docs/applications/cylinder_configurations.ipynb
index bea5692b..260a5c82 100644
--- a/docs/applications/cylinder_configurations.ipynb
+++ b/docs/applications/cylinder_configurations.ipynb
@@ -37,6 +37,58 @@
     "from visualisation import show_dislocation"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Accessing Key Dislocation Properties\n",
+    "As defined in a dimensionless (alat=1) frame, each kind of dislocation has several key properties which describe the geometric layout of the atomic cell, and also the position and orientation of the dislocation. Such properties are stored in matscipy as attributes of each dislocation class, such that they can be accessed without actually needing to construct a dislocation.\n",
+    "\n",
+    "As an example:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from matscipy.dislocation import BCCMixed111Dislocation\n",
+    "\n",
+    "print(\"For a BCC 1/2<111>{110} Mixed Dislocation\")\n",
+    "print()\n",
+    "\n",
+    "print(\"The dislocation is oriented in a cell described by the miller indices: \")\n",
+    "print(BCCMixed111Dislocation.axes)\n",
+    "\n",
+    "print(\"Dimensionless burgers vector is: \", BCCMixed111Dislocation.burgers_dimensionless)\n",
+    "\n",
+    "print(\"Dislocation core will be at fractional coords \", BCCMixed111Dislocation.unit_cell_core_position_dimensionless, \" within a unit cell\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Available Dislocation Systems\n",
+    "The available dislocation systems are:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matscipy.dislocation as disloc_mod\n",
+    "import inspect\n",
+    "# Find all classes in matscipy.dislocation which inherit from the Abstact Base Class CubicCrystalDislocation\n",
+    "for name, attr in disloc_mod.__dict__.items():\n",
+    "    if inspect.isclass(attr):\n",
+    "        if issubclass(attr, disloc_mod.CubicCrystalDislocation) and attr not in [disloc_mod.CubicCrystalDislocation, disloc_mod.CubicCrystalDissociatedDislocation]:\n",
+    "            print(name)"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
diff --git a/docs/applications/multispecies_dislocations.ipynb b/docs/applications/multispecies_dislocations.ipynb
new file mode 100644
index 00000000..e3f89837
--- /dev/null
+++ b/docs/applications/multispecies_dislocations.ipynb
@@ -0,0 +1,251 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Multi-species Dislocation Systems\n",
+    "The previous section discussed constructing cylindrical dislocation cells, using single-species systems as case studies. The dislocation classes in `matscipy.dislocation` can also be applied to systems which have much more chemical complexity. However, the following restrictions apply:\n",
+    "1. The desired system can be expressed in a cubic lattice\n",
+    "2. The desired system shares on-lattice geometry with an existing \"base\" crystal structure\n",
+    "\n",
+    "As an example, let's take Zincblende GaAs and compare with a diamond lattice:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from ase.build import bulk\n",
+    "from ase.lattice.cubic import Diamond\n",
+    "import numpy as np\n",
+    "import nglview # this import is necessary for rendering of the 3D view\n",
+    "from visualisation import show_dislocation, interactive_view\n",
+    "from ase.visualize import view\n",
+    "\n",
+    "# Data from https://doi.org/10.1080/08927022.2011.602975\n",
+    "alat = 11.306/2\n",
+    "C11 = 98.47\n",
+    "C12 = 15.25\n",
+    "C44 = 57.89\n",
+    "\n",
+    "GaAs = bulk(\"GaAs\", crystalstructure=\"zincblende\", cubic=True, a=alat)\n",
+    "\n",
+    "diamond = Diamond(\"C\", latticeconstant=alat)\n",
+    "\n",
+    "gaas_pos = GaAs.get_scaled_positions()\n",
+    "dia_pos = diamond.get_scaled_positions()\n",
+    "\n",
+    "# Sort the coords, as bulk uses a different atom order\n",
+    "gaas_idxs = np.lexsort((gaas_pos[:, 0], gaas_pos[:, 1], gaas_pos[:, 2]))\n",
+    "dia_idxs = np.lexsort((dia_pos[:, 0], dia_pos[:, 1], dia_pos[:, 2]))\n",
+    "print(\"GaAs Fractional Coordinates\")\n",
+    "print(gaas_pos[gaas_idxs, :])\n",
+    "print()\n",
+    "print(\"Diamond Fractional Coordinates\")\n",
+    "print(dia_pos[dia_idxs, :])\n",
+    "\n",
+    "interactive_view(GaAs, name=\"GaAs Bulk\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can see that the fractional coordinates of the Zincblende GaAs are the same as the Diamond structure we generated. We can therefore say that Zincblende shares geometry with Diamond. The displacement field used to generate dislocation structures is agnostic to atomic chemistry, and so we can model Zincblende dislocations as if they were dislocations in a Diamond crystal with the lattice constant and elastic properties of the Zincblende crystal.\n",
+    "\n",
+    "To build a dislocation with this GaAs bulk, we can simply pass it as an argument in place of the lattice constant:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from matscipy.dislocation import DiamondGlide90degreePartial\n",
+    "\n",
+    "disloc = DiamondGlide90degreePartial(GaAs, C11, C12, C44)\n",
+    "\n",
+    "GaAs_bulk, GaAs_dislocation = disloc.build_cylinder(radius=3 * alat)\n",
+    "\n",
+    "show_dislocation(GaAs_dislocation, scale=0.35, d_name=\"GaAs 1/6<112> 90 degree partial\", \n",
+    "                 diamond_structure=True, CNA_color=False, add_bonds=True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Additional Dislocations in Multispecies Systems\n",
+    "Dislocations in multispecies systems have additiional complexity over those in single-species crystals, due to the breaking of some symmetries caused by the added chemical complexity. For our Zincblende GaAs example, this means that we can have two different forms of some of our dislocations: $\\alpha$ (As-terminated), and $\\beta$ (Ga-terminated).\n",
+    "\n",
+    "$\\alpha$-$90^\\circ$ Partial Dislocation in GaAs:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# alpha-90 degree dislocation\n",
+    "GaAs = bulk(\"GaAs\", crystalstructure=\"zincblende\", cubic=True, a=alat)\n",
+    "symbols = np.array(GaAs.get_chemical_symbols())\n",
+    "# Swap Ga <-> As to get other alpha dislocation\n",
+    "new_symbols = symbols.copy()\n",
+    "new_symbols[symbols == \"Ga\"] = \"As\"\n",
+    "new_symbols[symbols == \"As\"] = \"Ga\"\n",
+    "GaAs.set_chemical_symbols(new_symbols)\n",
+    "disloc = DiamondGlide90degreePartial(GaAs, C11, C12, C44)\n",
+    "\n",
+    "GaAs_bulk, GaAs_dislocation = disloc.build_cylinder(radius=3 * alat)\n",
+    "view = show_dislocation(GaAs_dislocation, scale=0.35, d_name=\"GaAs 1/6<112> Alpha-90 degree partial\", \n",
+    "                 diamond_structure=True, CNA_color=False, add_bonds=True)\n",
+    "\n",
+    "view.control.zoom(0.8)\n",
+    "view"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "$\\beta$-$90^\\circ$ Partial Dislocation in GaAs:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "GaAs = bulk(\"GaAs\", crystalstructure=\"zincblende\", cubic=True, a=alat)\n",
+    "# beta dislocation is the one done by default, for the bulk built by ase.build.bulk\n",
+    "disloc = DiamondGlide90degreePartial(GaAs, C11, C12, C44)\n",
+    "\n",
+    "GaAs_bulk, GaAs_dislocation = disloc.build_cylinder(radius=3 * alat)\n",
+    "view = show_dislocation(GaAs_dislocation, scale=0.35, d_name=\"GaAs 1/6<112> Beta-90 degree partial\", \n",
+    "                 diamond_structure=True, CNA_color=False, add_bonds=True)\n",
+    "\n",
+    "view.control.zoom(0.8)\n",
+    "view"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Dislocations in Disordered Systems\n",
+    "Generating dislocation structures for systems which have some chemical disorder is also possible, with some caveats:\n",
+    "1. The generation routines will only work reliably when given on-lattice bulk structures with cubic cells\n",
+    "2. The generated structure is based on a periodic replication of the base unit_cell, thus won't have true disorder\n",
+    "3. The displacement field applied to the bulk structure does not depend on atomic species - known lattice distortions caused by the disorder will not be captured\n",
+    "4. The current code only allows $C_{11}$, $C_{12}$, and $C_{44}$ to be specified, so more complex elastic effects (e.g. $C_{11} \\neq C_{22}$) cannot currently be captured\n",
+    "\n",
+    "With these in mind, a recommended workflow would be to start by generating the dislocation system for an ordered system, but using the elastic constants of the disordered system. From there, the disorder can be applied by simply changing the chemcical symbols of the atoms to match the target composition.\n",
+    "\n",
+    "To show a worked example of this, consider the alloy $\\text{In}_{0.5} \\text{Ga}_{0.5} \\text{As}$. This should form in a Zincblende structure, where the Ga sites from the previous GaAs bulk are now 50% occupied by In.\n",
+    "\n",
+    "In order to model this, we can generate a dislocation for GaAs, using the lattice constant and elastic properties of $\\text{In}_{0.5} \\text{Ga}_{0.5} \\text{As}$.\n",
+    "\n",
+    "NOTE: Whilst disordered systems like this $\\text{In}_{0.5} \\text{Ga}_{0.5} \\text{As}$ example should have off-lattice distortions in the bulk state, it is heavily recommended that the dislocation structures are generated using an on-lattice crystal. This is because the off-lattice structure will interact differently with the continuum displacement field, which could lead to overlapping/extremely close atoms, or generally incorrect core structures. The off-lattice distortions should ideally be found by relaxing the dislocation structure."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from ase.build import bulk\n",
+    "from matscipy.dislocation import DiamondGlide90degreePartial\n",
+    "\n",
+    "# Data from https://doi.org/10.1080/08927022.2011.602975\n",
+    "alat = 11.2402/2\n",
+    "C11 = 120.31\n",
+    "C12 = 55.87\n",
+    "C44 = 58.26\n",
+    "\n",
+    "GaAs = bulk(\"GaAs\", crystalstructure=\"zincblende\", cubic=True, a=alat)\n",
+    "\n",
+    "\n",
+    "disloc = DiamondGlide90degreePartial(GaAs, C11, C12, C44)\n",
+    "\n",
+    "GaAs_bulk, GaAs_dislocation = disloc.build_cylinder(radius=3 * alat)\n",
+    "\n",
+    "show_dislocation(GaAs_dislocation, scale=0.35, d_name=\"GaAs 1/6<112> Beta-90 degree partial\", \n",
+    "                 diamond_structure=True, CNA_color=False, add_bonds=True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "From this GaAs structure, we then do Monte Carlo sampling to introduce 50% Indium to the structure:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "np.random.seed(1)\n",
+    "\n",
+    "species = np.array(GaAs_dislocation.get_chemical_symbols())\n",
+    "\n",
+    "Ga_idxs = np.argwhere(species == \"Ga\")[:, 0]\n",
+    "\n",
+    "# Choose random half of the idxs to be In\n",
+    "In_idxs = sorted(np.random.choice(Ga_idxs, size=int(Ga_idxs.shape[0]/2), replace=False))\n",
+    "\n",
+    "# Introduce the chemical disorder in In-Ga sites\n",
+    "species[In_idxs] = \"In\"\n",
+    "\n",
+    "InGaAs_bulk = GaAs_bulk.copy()\n",
+    "InGaAs_bulk.set_chemical_symbols(species)\n",
+    "\n",
+    "InGaAs_dislocation = GaAs_dislocation.copy()\n",
+    "InGaAs_dislocation.set_chemical_symbols(species)\n",
+    "\n",
+    "view = show_dislocation(InGaAs_dislocation, scale=0.35, d_name=\"InGaAs 1/6<112> Beta-90 degree partial\", \n",
+    "                        diamond_structure=True, CNA_color=False, add_bonds=True)\n",
+    "\n",
+    "# In and Ga have almost the same default colors in nglview\n",
+    "# so we add another component with the In atoms in red to\n",
+    "# see the chemical disorder better\n",
+    "In_ats = InGaAs_dislocation[In_idxs]\n",
+    "c = view.add_component(nglview.ASEStructure(In_ats),\n",
+    "                       default_representation=False)\n",
+    "c.add_spacefill(radius=1.0, color=\"red\")\n",
+    "view"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.12"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/docs/applications/plasticity.rst b/docs/applications/plasticity.rst
index 8d880457..c2c2dc41 100644
--- a/docs/applications/plasticity.rst
+++ b/docs/applications/plasticity.rst
@@ -17,4 +17,5 @@ Tutorials:
 
 .. toctree::
 
-    cylinder_configurations.ipynb
\ No newline at end of file
+    cylinder_configurations.ipynb
+    multispecies_dislocations.ipynb
\ No newline at end of file
diff --git a/docs/applications/visualisation.py b/docs/applications/visualisation.py
index e37590d5..1069346e 100644
--- a/docs/applications/visualisation.py
+++ b/docs/applications/visualisation.py
@@ -1,5 +1,6 @@
 import numpy as np
-
+# interactive visualisation  inside the notebook with nglview
+from nglview import show_ase, ASEStructure
 from ovito.io.ase import ase_to_ovito
 from ovito.modifiers import CommonNeighborAnalysisModifier, IdentifyDiamondModifier
 from ovito.pipeline import StaticSource, Pipeline
@@ -42,8 +43,30 @@ def get_structure_types(structure, diamond_structure=False):
 
     return atom_labels, structure_names, hex_colors
 
-# interactive visualisation  inside the notebook with nglview
-from nglview import show_ase, ASEStructure
+# custom tooltip for nglview to avoid showing molecule and residue names
+# that are not relevant for our type of structures
+tooltip_js = """
+                this.stage.mouseControls.add('hoverPick', (stage, pickingProxy) => {
+                    let tooltip = this.stage.tooltip;
+                    if(pickingProxy && pickingProxy.atom && !pickingProxy.bond){
+                        let atom = pickingProxy.atom;
+                        if (atom.structure.name.length > 0){
+                            tooltip.innerText = atom.atomname + " atom: " + atom.structure.name;
+                        } else {
+                            tooltip.innerText = atom.atomname + " atom";
+                        }
+                    } else if (pickingProxy && pickingProxy.bond){
+                        let bond = pickingProxy.bond;
+                        if (bond.structure.name.length > 0){
+                        tooltip.innerText = bond.atom1.atomname + "-" + bond.atom2.atomname + " bond: " + bond.structure.name;
+                        } else {
+                            tooltip.innerText = bond.atom1.atomname + "-" + bond.atom2.atomname + " bond";
+                        }
+                    } else if (pickingProxy && pickingProxy.unitcell){
+                        tooltip.innerText = "Unit cell";
+                    }
+                });
+                """
 
 def add_dislocation(view, system, name, color=[0, 1, 0], x_shift=0.0):
     '''Add dislocation line to the view as a cylinder and two cones.
@@ -71,7 +94,7 @@ def add_dislocation(view, system, name, color=[0, 1, 0], x_shift=0.0):
                         name)
     
 
-def show_dislocation(system, scale=0.5, diamond_structure=False, add_bonds=False,
+def show_dislocation(system, scale=0.5, CNA_color=True, diamond_structure=False, add_bonds=False,
                      d_name='', d_color=[0, 1, 0], partial_distance=None):
 
     atom_labels, structure_names, colors = get_structure_types(system, 
@@ -89,11 +112,17 @@ def show_dislocation(system, scale=0.5, diamond_structure=False, add_bonds=False
         mask = atom_labels == structure_type
         component = view.add_component(ASEStructure(system[mask]), 
                                        default_representation=False, name=str(structure_names[structure_type]))
-        if add_bonds:
-            component.add_ball_and_stick(color=colors[structure_type], radiusType='covalent', radiusScale=scale)
+        if CNA_color:
+            if add_bonds:
+                component.add_ball_and_stick(color=colors[structure_type], radiusType='covalent', radiusScale=scale)
+            else:
+                component.add_spacefill(color=colors[structure_type], radiusType='covalent', radiusScale=scale)
         else:
-            component.add_spacefill(color=colors[structure_type], radiusType='covalent', radiusScale=scale)
-        
+            if add_bonds:
+                component.add_ball_and_stick(radiusType='covalent', radiusScale=scale)
+            else:
+                component.add_spacefill(radiusType='covalent', radiusScale=scale)
+                    
     component.add_unitcell()
 
     if partial_distance is None:
@@ -111,20 +140,25 @@ def show_dislocation(system, scale=0.5, diamond_structure=False, add_bonds=False
 
     view._remote_call("setSize", target="Widget", args=["700px", "400px"])
     view.center()
-    
-    tooltip_js = """
-                 this.stage.mouseControls.add('hoverPick', (stage, pickingProxy) => {
-                        let tooltip = this.stage.tooltip;
-                        if(pickingProxy && pickingProxy.atom && !pickingProxy.bond){
-                            let atom = pickingProxy.atom;
-                            tooltip.innerText = atom.atomname + " atom: " + atom.structure.name;
-                        } else if (pickingProxy && pickingProxy.bond){
-                            let bond = pickingProxy.bond;
-                            tooltip.innerText = bond.atom1.atomname + "-" + bond.atom2.atomname + " bond: " + bond.structure.name;
-                        } else if (pickingProxy && pickingProxy.unitcell){
-                            tooltip.innerText = "Unit cell";
-                        }
-                    });
-                 """
+
+    view._js(tooltip_js)
+    return view
+
+
+def interactive_view(system, scale=0.5, name=""):
+    view = show_ase(system)
+    view._remove_representation()
+    component = view.add_component(ASEStructure(system), 
+                                   default_representation=False, name=name)
+    component.add_spacefill(radiusType='covalent', radiusScale=scale)
+    view.add_unitcell()
+    view.update_spacefill(radiusType='covalent',
+                          radiusScale=scale)
+
+    view.camera = 'orthographic'
+    view.parameters = {"clipDist": 0}
+
+    view.center()
+    view._remote_call("setSize", target="Widget", args=["300px", "300px"])
     view._js(tooltip_js)
     return view
\ No newline at end of file
diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index fc23479b..5ea63993 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -26,6 +26,8 @@
 
 import numpy as np
 
+from abc import ABCMeta
+
 from scipy.optimize import minimize
 
 from ase.lattice.cubic import (BodyCenteredCubic, FaceCenteredCubic,
@@ -40,6 +42,7 @@
 
 from matscipy.neighbours import neighbour_list, mic
 from matscipy.elasticity import fit_elastic_constants
+from matscipy.utils import validate_cubic_cell
 
 
 def make_screw_cyl(alat, C11, C12, C44,
@@ -2148,11 +2151,27 @@ def check_duplicates(atoms, distance=0.1):
     return duplicates.astype(np.bool)
 
 
-class CubicCrystalDislocation:
-    def __init__(self, unit_cell, alat, C11, C12, C44, axes, burgers,
-                 unit_cell_core_position=None,
-                 parity=None, glide_distance=None, n_planes=None,
-                 self_consistent=None):
+class CubicCrystalDislocation(metaclass=ABCMeta):
+    '''
+    Abstract base class for modelling a single dislocation
+    '''
+
+    # Mandatory Attributes of CubicCrystalDislocation with no defaults
+    # These should be set by the child dislocation class, or by CubicCrystalDislocation.__init__
+    # (see https://stackoverflow.com/questions/472000/usage-of-slots for more details on __slots__)
+    __slots__ = ("burgers_dimensionless", "unit_cell_core_position_dimensionless",
+                 "glide_distance_dimensionless", "crystalstructure", "axes",
+                 "C11", "C12", "C44", "alat", "unit_cell")
+
+    # Attributes with defaults
+    # These may be overridden by child classes
+    parity = np.zeros(2)
+    n_planes = 3
+    self_consistent = True
+    pbc = [True, True, True]
+    stroh = None
+
+    def __init__(self, a, C11, C12, C44, symbol="W"):
         """
         This class represents a dislocation in a cubic crystal
 
@@ -2162,64 +2181,95 @@ def __init__(self, unit_cell, alat, C11, C12, C44, axes, burgers,
 
         Parameters
         ----------
-        unit_cell : unit cell to build the dislocation configuration
-        alat : lattice constant
-        C11 : elastic constants
+        a : lattice constant OR cubic ase.atoms.Atoms object
+            Lattice constant passed to ase.lattice.cubic constructor or
+            Atoms object used by ase.build.cut to get unit cell in correct orientation
+        C11 : float
+            Elastic Constants
         C12
         C44
-        axes : cell axes (b is normally along z direction)
-        burgers : burgers vector of the dislocation
-        unit_cell_core_position : dislocation core position in the unit cell
-                                  used to shift atomic positions to
-                                  make the dislocation core the center
-                                  of the cell
-        parity
-        glide_distance : distance to the next equivalent
-                         core position in the glide direction
+        symbol : str
+            Chemical symbol used to construct unit cell (if "a" is a lattice constant)
+
+
+        Attributes
+        ------------------------------
+        alat : float
+            lattice parameter (in A)
+        unit_cell : ase.atoms.Atoms object
+            bulk cell used to generate dislocations
+        axes : np.array of float
+            unit cell axes (b is normally along z direction)
+        burgers_dimensionless : np.array of float
+            burgers vector of the dislocation class, in a dimensionless alat=1 system
+        burgers : np.array of float
+            burgers vector of the atomistic dislocation
+        unit_cell_core_position : np.array of float
+            dislocation core position in the unit cell used to shift atomic positions to
+            make the dislocation core the center of the cell
+        parity : np.array
+        glide_distance : float
+            distance (in A) to the next equivalent core position in the glide direction
         n_planes : int
             number of non equivalent planes in z direction
         self_consistent : float
             default value for the displacement calculation
+        crystalstructure : str
+            Name of basic structure defining atomic geometry ("fcc", "bcc", "diamond")
         """
 
-        self.unit_cell = unit_cell.copy()
-        self.alat = alat
+        # Create copies of immutable class attributes
+        # Prevents side effects when E.G. changing burgers vector of an instance
+        # Also makes changing cls.var and instance.var work as expected for these variables
+        self.axes = self.axes.copy()
+        self.burgers_dimensionless = self.burgers_dimensionless.copy()
+        self.unit_cell_core_position_dimensionless = self.unit_cell_core_position_dimensionless.copy()
+        self.parity = self.parity.copy()
+
         self.C11 = C11
         self.C12 = C12
         self.C44 = C44
+        
+        self.alat, self.unit_cell = validate_cubic_cell(a, symbol, self.axes, self.crystalstructure, self.pbc)
 
-        self.axes = axes
-        self.burgers = burgers
-        if unit_cell_core_position is None:
-            unit_cell_core_position = np.zeroes(3)
-        self.unit_cell_core_position = unit_cell_core_position
-        if parity is None:
-            parity = np.zeros(2, dtype=int)
-        self.parity = parity
-        if glide_distance is None:
-            glide_distance = 0.0
-        self.glide_distance = glide_distance
-        if n_planes is None:
-            n_planes = 3
-        self.n_planes = n_planes
-        if self_consistent is None:
-            self_consistent = True
-        self.self_consistent = self_consistent
-
-        self.stroh = None
+        self.glide_distance = self.alat * self.glide_distance_dimensionless
 
     def init_stroh(self):
-
         from atomman import ElasticConstants
         from atomman.defect import Stroh
         c = ElasticConstants(C11=self.C11, C12=self.C12, C44=self.C44)
         self.stroh = Stroh(c, burgers=self.burgers, axes=self.axes)
 
-    def set_burgers(self, burgers):
-        self.burgers = burgers
+    # @property & @var.setter decorators used to ensure var and var_dimensionless don't get out of sync
+    @property
+    def burgers(self):
+        return self.burgers_dimensionless * self.alat
+    
+    @burgers.setter
+    def burgers(self, burgers):
+        self.burgers_dimensionless = burgers / self.alat
         if self.stroh is None:
             self.init_stroh()
 
+    def set_burgers(self, burgers):
+        self.burgers = burgers
+
+    @property
+    def unit_cell_core_position(self):
+        return self.unit_cell_core_position_dimensionless * self.alat
+    
+    @unit_cell_core_position.setter
+    def unit_cell_core_position(self, position):
+        self.unit_cell_core_position_dimensionless = position / self.alat
+
+    @property
+    def glide_distance(self):
+        return self.glide_distance_dimensionless * self.alat
+    
+    @glide_distance.setter
+    def glide_distance(self, distance):
+        self.glide_distance_dimensionless = distance / self.alat
+
     def plot_unit_cell(self, ms=250, ax=None):
         import matplotlib.pyplot as plt
 
@@ -2509,193 +2559,23 @@ def build_impurity_cylinder(self, disloc, impurity, radius,
         return impurities_disloc
 
 
-class BCCScrew111Dislocation(CubicCrystalDislocation):
-    def __init__(self, alat, C11, C12, C44, symbol='W'):
-        axes = np.array([[1, 1, -2],
-                         [-1, 1, 0],
-                         [1, 1, 1]])
-        burgers = alat * np.array([1, 1, 1]) / 2.0
-        unit_cell_core_position = alat * np.array([np.sqrt(6.)/6.0,
-                                                   np.sqrt(2.)/6.0, 0])
-        parity = [0, 0]
-        unit_cell = BodyCenteredCubic(directions=axes.tolist(),
-                                      size=(1, 1, 1), symbol=symbol,
-                                      pbc=True,
-                                      latticeconstant=alat)
-        glide_distance = alat * np.linalg.norm(axes[0]) / 3.0
-        super().__init__(unit_cell, alat, C11, C12, C44,
-                         axes, burgers, unit_cell_core_position, parity,
-                         glide_distance)
-
-
-class BCCEdge111Dislocation(CubicCrystalDislocation):
-    def __init__(self, alat, C11, C12, C44, symbol='W'):
-        axes = np.array([[1, 1, 1],
-                         [1, -1, 0],
-                         [1, 1, -2]])
-        burgers = alat * np.array([1, 1, 1]) / 2.0
-        unit_cell_core_position = alat * np.array([(1.0/3.0) * np.sqrt(3.0)/2.0,
-                                                  0.25 * np.sqrt(2.0), 0])
-        parity = [0, 0]
-        unit_cell = BodyCenteredCubic(directions=axes.tolist(),
-                                      size=(1, 1, 1), symbol=symbol,
-                                      pbc=True,
-                                      latticeconstant=alat)
-        glide_distance = np.linalg.norm(burgers) / 3.0
-        n_planes = 6
-        super().__init__(unit_cell, alat, C11, C12, C44,
-                         axes, burgers, unit_cell_core_position, parity,
-                         glide_distance, n_planes=n_planes)
-
-
-class BCCMixed111Dislocation(CubicCrystalDislocation):
-    def __init__(self, alat, C11, C12, C44, symbol='W'):
-        axes = np.array([[1, -1, -2],
-                         [1, 1, 0],
-                         [1, -1, 1]])
-        burgers = alat * np.array([1, -1, -1]) / 2.0
-        parity = [0, 0]
-        unit_cell = BodyCenteredCubic(directions=axes.tolist(),
-                                      size=(1, 1, 1), symbol=symbol,
-                                      pbc=True,
-                                      latticeconstant=alat)
-
-        # middle of the right edge of the first upward triangle
-        core_position = (unit_cell.positions[1] +
-                         unit_cell.positions[2]) / 2.0
-
-        unit_cell_core_position = np.array([core_position[0],
-                                            core_position[1], 0])
-
-        glide_distance = alat * np.linalg.norm(axes[0]) / 3.0
-        super().__init__(unit_cell, alat, C11, C12, C44,
-                         axes, burgers, unit_cell_core_position, parity,
-                         glide_distance)
-
-
-class BCCEdge100Dislocation(CubicCrystalDislocation):
-    def __init__(self, alat, C11, C12, C44, symbol='W'):
-        axes = np.array([[1, 0, 0],
-                         [0, 0, -1],
-                         [0, 1, 0]])
-        burgers = alat * np.array([1, 0, 0])
-        unit_cell_core_position = alat * np.array([0.25,
-                                                   0.25, 0])
-        parity = [0, 0]
-        unit_cell = BodyCenteredCubic(directions=axes.tolist(),
-                                      size=(1, 1, 1), symbol=symbol,
-                                      pbc=True,
-                                      latticeconstant=alat)
-        glide_distance = alat
-        n_planes = 2
-        super().__init__(unit_cell, alat, C11, C12, C44,
-                         axes, burgers, unit_cell_core_position, parity,
-                         glide_distance, n_planes=n_planes)
-
-
-class BCCEdge100110Dislocation(CubicCrystalDislocation):
-    def __init__(self, alat, C11, C12, C44, symbol='W'):
-        axes = np.array([[1, 0, 0],
-                         [0, 1, 1],
-                         [0, -1, 1]])
-        burgers = alat * np.array([1, 0, 0])
-        unit_cell_core_position = alat * np.array([0.5,
-                                                   np.sqrt(2.) / 4.0, 0])
-        parity = [0, 0]
-        unit_cell = BodyCenteredCubic(directions=axes.tolist(),
-                                      size=(1, 1, 1), symbol=symbol,
-                                      pbc=True,
-                                      latticeconstant=alat)
-        glide_distance = 0.5 * alat
-        n_planes = 2
-        super().__init__(unit_cell, alat, C11, C12, C44,
-                         axes, burgers, unit_cell_core_position, parity,
-                         glide_distance, n_planes=n_planes)
-
-
-class DiamondGlide30degreePartial(CubicCrystalDislocation):
-    def __init__(self, alat, C11, C12, C44, symbol='C'):
-        axes = np.array([[1, 1, -2],
-                         [1, 1, 1],
-                         [1, -1, 0]])
-
-        burgers = alat * np.array([1, -2, 1.]) / 6.
-
-        disloCenterX = 0.5 * (alat * np.linalg.norm(axes[0])) / 6.0
-        # 1/4 + 1/2 * (1/3 - 1/4) - to be in the middle of the glide set
-        disloCenterY = 7.0 * (alat * np.linalg.norm(axes[1])) / 24.0
-
-        unit_cell_core_position = np.array([disloCenterX,
-                                            disloCenterY, 0])
-
-        parity = [0, 0]
-
-        unit_cell = Diamond(symbol, directions=axes.tolist(),
-                            pbc=(False, False, True),
-                            latticeconstant=alat)
-
-        glide_distance = alat * np.linalg.norm(axes[0]) / 4.0
-
-        n_planes = 2
-        # There is very small distance between
-        # atomic planes in glide configuration.
-        # Due to significant anisotropy application of the self consistent
-        # displacement field leads to deformation of the atomic planes.
-        # This leads to the cut plane crossing one of the atomic planes and
-        # thus breaking the stacking fault.
-        self_consistent = False
-        super().__init__(unit_cell, alat, C11, C12, C44,
-                         axes, burgers, unit_cell_core_position, parity,
-                         glide_distance, n_planes=n_planes,
-                         self_consistent=self_consistent)
-
-
-class DiamondGlide90degreePartial(CubicCrystalDislocation):
-    def __init__(self, alat, C11, C12, C44, symbol='C'):
-        axes = np.array([[1, 1, -2],
-                         [1, 1, 1],
-                         [1, -1, 0]])
-
-        burgers = alat * np.array([1., 1., -2.]) / 6.
-
-        disloCenterX = 0.5 * (alat * np.linalg.norm(axes[0])) / 6.0
-        # 1/4 + 1/2 * (1/3 - 1/4) - to be in the middle of the glide set
-        disloCenterY = 7.0 * (alat * np.linalg.norm(axes[1])) / 24.0
-
-        unit_cell_core_position = np.array([disloCenterX,
-                                            disloCenterY, 0])
-
-        parity = [0, 0]
-
-        unit_cell = Diamond(symbol, directions=axes.tolist(),
-                            pbc=(False, False, True),
-                            latticeconstant=alat)
-
-        glide_distance = alat * np.linalg.norm(axes[0]) / 4.0
-
-        n_planes = 2
-        # There is very small distance between
-        # atomic planes in glide configuration.
-        # Due to significant anisotropy application of the self consistent
-        # displacement field leads to deformation of the atomic planes.
-        # This leads to the cut plane crossing one of the atomic planes and
-        # thus breaking the stacking fault.
-        self_consistent = False
-        super().__init__(unit_cell, alat, C11, C12, C44,
-                         axes, burgers, unit_cell_core_position, parity,
-                         glide_distance, n_planes=n_planes,
-                         self_consistent=self_consistent)
-
+class CubicCrystalDissociatedDislocation(CubicCrystalDislocation, metaclass=ABCMeta):
+    '''
+    Abstract base class for modelling dissociated dislocation systems
+    '''
+    # Inherits all slots from CubicCrystalDislocation as well
+    __slots__ = ("left_dislocation", "right_dislocation")
 
-class CubicCrystalDissociatedDislocation(CubicCrystalDislocation):
-    def __init__(self, left_dislocation, right_dislocation, burgers):
-        """This class represents a dissociated dislocation in a cubic crystal
-        with burgers vercor b = b_left + b_right.
+    # Space for overriding the burgers vectors from cls.left_dislocation.burgers_dimensionless
+    new_left_burgers = None
+    new_right_burgers = None
+    def __init__(self, a, C11, C12, C44, symbol="W"):
+        """
+        This class represents a dissociated dislocation in a cubic crystal
+        with burgers vector b = b_left + b_right.
 
         Args:
-            left_dislocation (CubicCrystalDislocation): dislocation with b_left
-            right_dislocation (CubicCrystalDislocation): dislocation with b_right
-            burgers (ndarray of 3 floats): resulting burgers vector
+            identical to CubicCrystalDislocation
 
         Raises:
             ValueError: If resulting burgers vector
@@ -2704,10 +2584,34 @@ def __init__(self, left_dislocation, right_dislocation, burgers):
             ValueError: If one of the properties of
                         left and righ dislocations are not the same.
         """
+
+        self.left_dislocation = self.left_dislocation(a, C11, C12, C44, symbol)
+        self.right_dislocation = self.right_dislocation(a, C11, C12, C44, symbol)
+
+        # Change disloc burgers vectors, if requested
+        if self.new_left_burgers is not None:
+            self.left_dislocation.burgers_dimensionless = self.new_left_burgers.copy()
+        if self.new_right_burgers is not None:
+            self.right_dislocation.burgers_dimensionless = self.new_right_burgers.copy()
+        
+        # Set self.attrs based on left disloc attrs
+        left_dislocation = self.left_dislocation
+        right_dislocation = self.right_dislocation
+        self.crystalstructure = left_dislocation.crystalstructure
+        self.axes = left_dislocation.axes.copy()
+        self.unit_cell_core_position_dimensionless = left_dislocation.unit_cell_core_position_dimensionless.copy()
+        self.parity = left_dislocation.parity
+        self.glide_distance_dimensionless = left_dislocation.glide_distance_dimensionless
+        self.n_planes = left_dislocation.n_planes
+        self.self_consistent = left_dislocation.self_consistent
+
+        super().__init__(a, C11, C12, C44, symbol)
+
+        # Validation of disloc inputs
         try:
             np.testing.assert_almost_equal(left_dislocation.burgers +
                                            right_dislocation.burgers,
-                                           burgers)
+                                           self.burgers)
         except AssertionError as error:
             print(error)
             raise ValueError("Burgers vectors of left and right disloctions" +
@@ -2749,21 +2653,6 @@ def __init__(self, left_dislocation, right_dislocation, burgers):
             raise ValueError("Parameters of left and right" +
                              "partials must be the same")
 
-        self.left_dislocation = left_dislocation
-        self.right_dislocation = right_dislocation
-
-        super().__init__(left_dislocation.unit_cell,
-                         left_dislocation.alat,
-                         left_dislocation.C11,
-                         left_dislocation.C12,
-                         left_dislocation.C44,
-                         left_dislocation.axes,
-                         burgers,
-                         unit_cell_core_position=left_dislocation.unit_cell_core_position,
-                         parity=left_dislocation.parity,
-                         glide_distance=right_dislocation.glide_distance,
-                         n_planes=left_dislocation.n_planes,
-                         self_consistent=left_dislocation.self_consistent)
 
     def build_cylinder(self, radius, partial_distance=0,
                        core_position=np.array([0., 0., 0.]),
@@ -2827,148 +2716,148 @@ def displacements(self, bulk_positions, center,
         return left_u + right_u
 
 
-class DiamondGlideScrew(CubicCrystalDissociatedDislocation):
-    def __init__(self, alat, C11, C12, C44, symbol='C'):
-
-        # aiming for the resulting burgers vector
-        burgers = alat * np.array([1, -1, 0]) / 2.
-
-        # 30 degree
-        burgers_left = alat * np.array([2., -1., -1.]) / 6.
-
-        left30 = DiamondGlide30degreePartial(alat, C11, C12, C44,
-                                             symbol=symbol)
-        left30.set_burgers(burgers_left)
-        # another 30 degree
-        # burgers_right = alat * np.array([1, -2, 1.]) / 6. - default value
-        right30 = DiamondGlide30degreePartial(alat, C11, C12, C44,
-                                              symbol=symbol)
-
-        super().__init__(left30, right30, burgers)
+class BCCScrew111Dislocation(CubicCrystalDislocation):
+    crystalstructure = "bcc"
+    axes = np.array([[1, 1, -2],
+                     [-1, 1, 0],
+                     [1, 1, 1]])
+    burgers_dimensionless = np.array([1, 1, 1]) / 2.0
+    unit_cell_core_position_dimensionless = np.array([np.sqrt(6.)/6.0, np.sqrt(2.)/6.0, 0])
+    glide_distance_dimensionless = np.sqrt(6) / 3.0
 
 
-class DiamondGlide60Degree(CubicCrystalDissociatedDislocation):
-    def __init__(self, alat, C11, C12, C44, symbol='C'):
+class BCCEdge111Dislocation(CubicCrystalDislocation):
+    crystalstructure = "bcc"
+    axes = np.array([[1, 1, 1],
+                     [1, -1, 0],
+                     [1, 1, -2]])
+    burgers_dimensionless = np.array([1, 1, 1]) / 2.0
+    unit_cell_core_position_dimensionless =  np.array([(1.0/3.0) * np.sqrt(3.0)/2.0, 0.25 * np.sqrt(2.0), 0])
+    glide_distance_dimensionless = np.sqrt(3) / 3.0
+    n_planes = 6
 
-        # aiming for the resulting burgers vector
-        burgers = alat * np.array([1, 0, -1]) / 2.
 
-        # 30 degree
-        burgers_left = alat * np.array([2., -1., -1.]) / 6.
-        left30 = DiamondGlide30degreePartial(alat, C11, C12, C44,
-                                             symbol=symbol)
-        left30.set_burgers(burgers_left)
-        # 90 degree
-        # burgers_right = alat * np.array([1, 1, -2.]) / 6. - default value
-        right90 = DiamondGlide90degreePartial(alat, C11, C12, C44,
-                                              symbol=symbol)
+class BCCMixed111Dislocation(CubicCrystalDislocation):
+    crystalstructure = "bcc"
+    axes = np.array([[1, -1, -2],
+                     [1, 1, 0],
+                     [1, -1, 1]])
+    burgers_dimensionless = np.array([1, -1, -1]) / 2.0
+    glide_distance_dimensionless = np.sqrt(6) / 3.0
+    # middle of the right edge of the first upward triangle
+    # half way between (1/6, 1/2, 0) and (1/3, 0, 0) in fractional coords
+    unit_cell_core_position_dimensionless = np.array([1/6, 1/6, 0])
 
-        super().__init__(left30, right90, burgers)
 
+class BCCEdge100Dislocation(CubicCrystalDislocation):
+    crystalstructure = "bcc"
+    axes = np.array([[1, 0, 0],
+                     [0, 0, -1],
+                     [0, 1, 0]])
+    burgers_dimensionless = np.array([1, 0, 0])
+    unit_cell_core_position_dimensionless = np.array([1/4, 1/4, 0])
+    glide_distance_dimensionless = 1.0
+    n_planes = 2
 
-class FCCScrewShockleyPartial(CubicCrystalDislocation):
-    def __init__(self, alat, C11, C12, C44, symbol='Fe'):
-        axes = np.array([[1, 1, -2],
-                         [1, 1, 1],
-                         [1, -1, 0]])
 
-        burgers = alat * np.array([1, -2, 1.]) / 6.
+class BCCEdge100110Dislocation(CubicCrystalDislocation):
+    crystalstructure = "bcc"
+    axes = np.array([[1, 0, 0],
+                     [0, 1, 1],
+                     [0, -1, 1]])
+    burgers_dimensionless = np.array([1, 0, 0])
+    unit_cell_core_position_dimensionless = np.array([0.5, np.sqrt(2) / 4.0, 0])
+    glide_distance_dimensionless = 0.5
+    n_planes = 2
 
-        parity = [0, 0]
 
-        unit_cell = FaceCenteredCubic(symbol, directions=axes.tolist(),
-                                      pbc=(False, False, True),
-                                      latticeconstant=alat)
+class DiamondGlide30degreePartial(CubicCrystalDislocation):
+    crystalstructure="diamond"
+    axes = np.array([[1, 1, -2],
+                     [1, 1, 1],
+                     [1, -1, 0]])
+    burgers_dimensionless = np.array([1, -2, 1.]) / 6.
+    glide_distance_dimensionless = np.sqrt(6) / 4.0
+    # 1/4 + 1/2 * (1/3 - 1/4) - to be in the middle of the glide set
+    unit_cell_core_position_dimensionless = np.array([np.sqrt(6)/12, 7*np.sqrt(3)/24, 0])
+    n_planes = 2
+    # There is very small distance between
+    # atomic planes in glide configuration.
+    # Due to significant anisotropy application of the self consistent
+    # displacement field leads to deformation of the atomic planes.
+    # This leads to the cut plane crossing one of the atomic planes and
+    # thus breaking the stacking fault.
+    self_consistent = False
 
-        # put the dislocation in the centroid of the first triangle in xy plane
-        # get sorting of the atoms in the distance of atoms in xy plane
-        sorted_indices = np.argsort(np.linalg.norm(unit_cell.positions[:, :2],
-                                                   axis=1))
-        # centroid coordinates are simply
-        # mean of x y coordinates of the trianlge
-        (disloCenterX,
-         disloCenterY,
-         _) = unit_cell.positions[sorted_indices[:3]].mean(axis=0)
 
-        unit_cell_core_position = np.array([disloCenterX,
-                                            disloCenterY, 0])
+class DiamondGlide90degreePartial(CubicCrystalDislocation):
+    crystalstructure = "diamond"
+    axes = np.array([[1, 1, -2],
+                     [1, 1, 1],
+                     [1, -1, 0]])
+    # 1/4 + 1/2 * (1/3 - 1/4) - to be in the middle of the glide set
+    unit_cell_core_position_dimensionless = np.array([np.sqrt(6)/12, 7 * np.sqrt(3) / 24, 0])
+    burgers_dimensionless = np.array([1., 1., -2.]) / 6.
+    glide_distance_dimensionless = np.sqrt(6)/4
+    n_planes = 2
+    # There is very small distance between
+    # atomic planes in glide configuration.
+    # Due to significant anisotropy application of the self consistent
+    # displacement field leads to deformation of the atomic planes.
+    # This leads to the cut plane crossing one of the atomic planes and
+    # thus breaking the stacking fault.
+    self_consistent = False
 
-        glide_distance = alat * np.linalg.norm(axes[0]) / 4.0
 
-        n_planes = 2
-        self_consistent = True
-        super().__init__(unit_cell, alat, C11, C12, C44,
-                         axes, burgers, unit_cell_core_position, parity,
-                         glide_distance, n_planes=n_planes,
-                         self_consistent=self_consistent)
+class DiamondGlideScrew(CubicCrystalDissociatedDislocation):
+    burgers_dimensionless = np.array([1, -1, 0]) / 2
+    new_left_burgers = np.array([2., -1., -1.]) / 6
+    left_dislocation = DiamondGlide30degreePartial
+    right_dislocation = DiamondGlide30degreePartial
 
 
-class FCCScrew110Dislocation(CubicCrystalDissociatedDislocation):
-    def __init__(self, alat, C11, C12, C44, symbol='Fe'):
+class DiamondGlide60Degree(CubicCrystalDissociatedDislocation):
+    burgers_dimensionless = np.array([1, 0, -1]) / 2
+    new_left_burgers = np.array([2., -1., -1.]) / 6
+    left_dislocation = DiamondGlide30degreePartial
+    right_dislocation = DiamondGlide90degreePartial
 
-        # aiming for the resulting burgers vector
-        burgers = alat * np.array([1, -1, 0]) / 2.
 
-        # Shockley partial
-        burgers_left = alat * np.array([2., -1., -1.]) / 6.
-        left_shockley = FCCScrewShockleyPartial(alat, C11, C12, C44,
-                                                symbol=symbol)
-        left_shockley.set_burgers(burgers_left)
-        # another Shockley partial
-        # burgers_right = alat * np.array([1, -2, 1.]) / 6. - default value
-        right_shockley = FCCScrewShockleyPartial(alat, C11, C12, C44,
-                                                 symbol=symbol)
-        super().__init__(left_shockley, right_shockley, burgers)
+class FCCScrewShockleyPartial(CubicCrystalDislocation):
+    crystalstructure="fcc"
+    axes = np.array([[1, 1, -2],
+                     [1, 1, 1],
+                     [1, -1, 0]])
+    burgers_dimensionless = np.array([1, -2, 1.]) / 6.
+    glide_distance_dimensionless = np.sqrt(6) / 4.0
+    n_planes = 2
+    unit_cell_core_position_dimensionless = np.array([5/6, 1/9, 0])
 
 
 class FCCEdgeShockleyPartial(CubicCrystalDislocation):
-    def __init__(self, alat, C11, C12, C44, symbol='Fe'):
-        axes = np.array([[1, -1, 0],
-                         [1, 1, 1],
-                         [-1, -1, 2]])
-
-        burgers = alat * np.array([1, -2, 1.]) / 6.
+    crystalstructure = "fcc"
+    axes = np.array([[1, -1, 0],
+                     [1, 1, 1],
+                     [-1, -1, 2]])
+    burgers_dimensionless = np.array([1, -2, 1.]) / 6
+    unit_cell_core_position_dimensionless = np.array([0, 1/6 , 0])
+    glide_distance_dimensionless = np.sqrt(2)/4
+    n_planes = 6
 
-        parity = [0, 0]
 
-        unit_cell = FaceCenteredCubic(symbol, directions=axes.tolist(),
-                                      pbc=(False, False, True),
-                                      latticeconstant=alat)
-
-        disloCenterX = 0.0
-        # middle between two (111) planes
-        disloCenterY = unit_cell.cell[1][1] / 6.0
-
-        unit_cell_core_position = np.array([disloCenterX,
-                                            disloCenterY, 0])
-
-        glide_distance = alat * np.linalg.norm(axes[0]) / 4.0
-
-        n_planes = 6
-        self_consistent = True
-        super().__init__(unit_cell, alat, C11, C12, C44,
-                         axes, burgers, unit_cell_core_position, parity,
-                         glide_distance, n_planes=n_planes,
-                         self_consistent=self_consistent)
+class FCCScrew110Dislocation(CubicCrystalDissociatedDislocation):
+    burgers_dimensionless = np.array([1, -1, 0]) / 2
+    new_left_burgers = np.array([2., -1., -1.]) / 6
+    left_dislocation = FCCScrewShockleyPartial
+    right_dislocation = FCCScrewShockleyPartial
 
 
 class FCCEdge110Dislocation(CubicCrystalDissociatedDislocation):
-    def __init__(self, alat, C11, C12, C44, symbol='Fe'):
-
-        # aiming for the resulting burgers vector
-        burgers = alat * np.array([1, -1, 0]) / 2.
-
-        # Shockley partial
-        burgers_left = alat * np.array([2., -1., -1.]) / 6.
-        left_shockley = FCCEdgeShockleyPartial(alat, C11, C12, C44,
-                                               symbol=symbol)
-        left_shockley.set_burgers(burgers_left)
-        # another Shockley partial
-        # burgers_right = alat * np.array([1, -2, 1.]) / 6. - default value
-        right_shockley = FCCEdgeShockleyPartial(alat, C11, C12, C44,
-                                                symbol=symbol)
-
-        super().__init__(left_shockley, right_shockley, burgers)
+    crystalstructure = "fcc"
+    burgers_dimensionless = np.array([1, -1, 0]) / 2
+    new_left_burgers = np.array([2., -1., -1.]) / 6
+    left_dislocation = FCCEdgeShockleyPartial
+    right_dislocation = FCCEdgeShockleyPartial
 
 
 class FixedLineAtoms:
diff --git a/matscipy/meson.build b/matscipy/meson.build
index a8e6f488..83e93c5b 100644
--- a/matscipy/meson.build
+++ b/matscipy/meson.build
@@ -37,7 +37,8 @@ python_sources = [
     'surface.py',
     'visualise.py',
     'cauchy_born.py',
-    'surface_reconstruction.py'
+    'surface_reconstruction.py',
+    'utils.py'
 ]
 
 # Install pure Python
diff --git a/matscipy/utils.py b/matscipy/utils.py
new file mode 100644
index 00000000..9be8efc9
--- /dev/null
+++ b/matscipy/utils.py
@@ -0,0 +1,107 @@
+import numpy as np
+
+def validate_cubic_cell(a, symbol="w", axes=None, crystalstructure=None, pbc=True):
+    '''
+    Provide uniform interface for generating rotated atoms objects through two main methods:
+
+    For lattice constant + symbol + crystalstructure, build a valid cubic bulk rotated into the frame defined by axes
+    For cubic bulk atoms object + crystalstructure, validate structure matches expected cubic bulk geometry, and rotate to frame defined by axes
+
+    a: float or ase.atoms
+        EITHER lattice constant (in A), or the cubic bulk structure
+    symbol: str
+        Elemental symbol, passed to relevant crystalstructure generator when a is float
+    axes: np.array
+        Axes transform to apply to bulk cell
+    crystalstructure: str
+        Base Structure of bulk system, currently supported: fcc, bcc, diamond
+    pbc: list of bool
+        Periodic Boundary Conditions in x, y, z
+    '''
+
+    from ase.lattice.cubic import FaceCenteredCubic, BodyCenteredCubic, Diamond
+    from ase.atoms import Atoms
+    from ase.build import cut, rotate
+    from warnings import warn
+
+    constructors = {
+        "fcc": FaceCenteredCubic,
+        "bcc": BodyCenteredCubic,
+        "diamond": Diamond
+    }
+
+    # Choose correct ase.lattice.cubic constructor given crystalstructure
+    try:
+        cell_builder = constructors[crystalstructure.lower()]
+    except (KeyError, AttributeError) as e:
+        if type(e) == KeyError:
+            # KeyError when crystalstructure not a key in constructors
+            raise ValueError(f"Crystal Structure {crystalstructure} not in accepted list: {list(constructors.keys())}")
+        else:
+            # AttributeError when type(crystalstructure) doesn't support .lower() (not a valid input type)
+            raise TypeError(f"crystalstructure should be one of {constructors.keys()}")
+
+
+    if np.issubdtype(type(a), np.floating) or np.issubdtype(type(a), np.integer):
+        # Reproduce legacy behaviour with a==alat
+        alat = a
+            
+        unit_cell = cell_builder(symbol, directions=axes.tolist(),
+                                 pbc=pbc,
+                                 latticeconstant=alat)
+    elif type(a) == Atoms:
+        # New behaviour for arbitrary cubic unit cells (Zincblende, L12, ...)
+        alat = a.cell[0, 0]
+
+        ats = a.copy()
+
+        if crystalstructure is not None:
+            # Try to validate that "a" matches expected structure given by crystalstructure
+            tol = 1e-3
+            ref_ats = cell_builder("C", directions=np.eye(3).astype(int).tolist(),
+                                latticeconstant=alat)
+
+            # Check that number of atoms in ats is an integer multiple of ref_ats
+            # (ats has to be the cubic primitive, or a supercell of the cubic primitive)
+            # Also check that the structure is cubic
+            rel_size = len(ats) / len(ref_ats)
+            skip_fractional_check = False
+            try:
+                # Integer multiple of ats test
+                assert abs(rel_size - np.floor(rel_size)) < tol
+
+                # Cubic cell test: cell = a*I_3x3
+                assert np.allclose(ats.cell[:, :], np.eye(3) * ats.cell[0, 0], atol=tol)
+            except AssertionError:
+                # Test failed, skip next test + warn user
+                skip_fractional_check = True
+                
+            # Check fractional coords match expectation
+            frac_match = True
+            sup_size = int(rel_size**(1/3))
+            if not skip_fractional_check:
+                ref_supercell = ref_ats * (sup_size, sup_size, sup_size)
+
+                try:
+                    apos = np.sort(ats.get_scaled_positions(), axis=0)
+                    rpos = np.sort(ref_supercell.get_scaled_positions(), axis=0)
+                    assert np.allclose(apos, rpos, atol=tol)
+                except:
+                    frac_match = False
+
+            if not frac_match or skip_fractional_check:
+                warn(f"Input bulk does not appear to match bulk {crystalstructure}.", stacklevel=2)
+
+            alat /= sup_size # Account for larger supercells having multiple internal core sites
+
+        ats = cut(ats, a=axes[0, :], b=axes[1, :], c=axes[2, :])
+        rotate(ats, ats.cell[0, :].copy(), [1, 0, 0], ats.cell[1, :].copy(), [0, 1, 0])
+        unit_cell = ats.copy()
+
+        unit_cell.set_pbc(pbc)
+    
+    else:
+        # "a" variable not valid
+        raise TypeError("type(a) is not a float, or an Atoms object")
+
+    return alat, unit_cell
diff --git a/paper/paper.md b/paper/paper.md
index 3d7003f4..0a9a2142 100644
--- a/paper/paper.md
+++ b/paper/paper.md
@@ -111,7 +111,7 @@ Within materials science, the package has different application domains:
 
 - **Elasticity.** Solids respond to small external loads through a reversible elastic response. The strength of the response is characterized by the elastic moduli. `matscipy.elasticity` implements functions for computing elastic moduli from small deformation that consider potential symmetries of the underlying atomic system, in particular for crystals. The implementation also includes estimates of uncertainty on elastic moduli - either from a least-squares error, or from a Bayesian treatment if stress uncertainty is supplied. `matscipy` also implements analytic calculation of elastic moduli for some interatomic potentials, described in more detail below. The computation of elastic moduli is a prerequisite for multi-scale modelling of materials, as they are the most basic parameters of continuum material models. `matscipy` was used to study finite-pressure elastic constants and structural stability in crystals [@Griesser2023crystal] and glasses [@Griesser2023glass].
 
-- **Plasticity.** For large loads, solids can respond with irreversible deformation. One form of irreversibility is plasticity, that is carried by extended defects, the dislocations, in crystals. The module `matscipy.dislocation` implements tools for studying structure and movement of dislocations. Construction and analysis of model atomic systems is implemented for compact and dissociated screw, as well as edge dislocations in cubic crystals. The implementation supports ideal straight as well as kinked dislocations. Some of the dislocation functionality requires the `atomman` and/or `OVITO` packages as optional dependencies [@AtomMan;@Stukowski2009]. The module was employed in a study of interaction of impurities with screw dislocations in tungsten [@Grigorev2020;@Grigorev2023]. 
+- **Plasticity.** For large loads, solids can respond with irreversible deformation. One form of irreversibility is plasticity, that is carried by extended defects, the dislocations, in crystals. The module `matscipy.dislocation` implements tools for studying structure and movement of dislocations. Construction and analysis of model atomic systems is implemented for compact and dissociated screw, as well as edge dislocations in cubic crystals. The implementation supports ideal straight as well as kinked dislocations. Some of the dislocation functionality requires the `atomman` and/or `OVITO` packages as optional dependencies [@AtomMan;@Stukowski2009]. The toolkit can be applied to a wide range of single- and multi-component ordered systems, and could be used as an initial starting point for modelling dislocations in systems with chemical disorder. The module was employed in a study of interaction of impurities with screw dislocations in tungsten [@Grigorev2020;@Grigorev2023]. 
 
 - **Fracture mechanics.** Cracking is the process of generating new surface area by splitting the material apart. The module `matscipy.fracture_mechanics` provides functionality for calculating continuum linear elastic displacement fields near crack tips, including support for anisotropy in the elastic response [@Sih1965]. The module also implements generation of atomic structures that are deformed according to this near-tip field. This functionality has been used to quantify lattice trapping, which is the pinning of cracks due to the discreteness of the atomic lattice, and to compare simulations with experimental measurements of crack speeds in silicon [@Kermode2015]. There is also support for flexible boundary conditions in fracture simulations using the formalism proposed by Sinclair [@Sinclair1975], where the finite atomistic domain is coupled to an infinite elastic continuum. Finally, we provide an extension of this approach to give a flexible boundary scheme that uses numerical continuation to obtain full solution paths for cracks [@Buze2021].
 
diff --git a/tests/test_dislocation.py b/tests/test_dislocation.py
index 38e80823..231e953f 100644
--- a/tests/test_dislocation.py
+++ b/tests/test_dislocation.py
@@ -30,6 +30,8 @@
 from ase.calculators.lammpslib import LAMMPSlib
 from matscipy.calculators.eam import EAM
 
+from ase.build import bulk as ase_bulk
+
 test_dir = os.path.dirname(os.path.realpath(__file__))
 
 try:
@@ -410,13 +412,18 @@ def test_slice_long_dislo(self):
         self.assertEqual(len(sliced_left_kink), kink_length * 3 - 1)
 
     def check_disloc(self, cls, ref_angle, structure="BCC", test_u=True,
-                     burgers=0.5 * np.array([1.0, 1.0, 1.0]), tol=10.0):
+                     burgers=0.5 * np.array([1.0, 1.0, 1.0]), tol=10.0, gen_bulk=False):
         alat = 3.14339177996466
         C11 = 523.0266819809012
         C12 = 202.1786296941397
         C44 = 160.88179872237012
 
-        d = cls(alat, C11, C12, C44, symbol="Al")
+        if gen_bulk:
+            a = ase_bulk("Al", structure.lower(), alat, cubic=True)
+        else:
+            a = alat
+
+        d = cls(a, C11, C12, C44, symbol="Al")
         bulk, disloc = d.build_cylinder(20.0)
 
         # test that assigning non default symbol worked
@@ -493,6 +500,10 @@ def test_90degree_diamond_partial_dislocation(self,):
         self.check_disloc(sd.DiamondGlide90degreePartial, 90.0,
                           structure="Diamond",
                           burgers=(1.0 / 6.0) * np.array([2.0, 1.0, 1.0]))
+        self.check_disloc(sd.DiamondGlide90degreePartial, 90.0,
+                          structure="Diamond",
+                          burgers=(1.0 / 6.0) * np.array([2.0, 1.0, 1.0]),
+                          gen_bulk=True)
 
     @unittest.skipIf("atomman" not in sys.modules or
                      "ovito" not in sys.modules,
@@ -731,6 +742,5 @@ def test_gamma_line(self):
         np.testing.assert_almost_equal(E, target_E, decimal=3)
         np.testing.assert_almost_equal(shift, target_shift, decimal=3)
 
-
 if __name__ == '__main__':
     unittest.main()
diff --git a/tests/test_utils.py b/tests/test_utils.py
new file mode 100644
index 00000000..5ca80626
--- /dev/null
+++ b/tests/test_utils.py
@@ -0,0 +1,90 @@
+import os
+import unittest
+import matscipytest
+import sys
+
+import matscipy.utils as utils_mod
+import numpy as np
+
+
+test_dir = os.path.dirname(os.path.realpath(__file__))
+
+
+class TestDislocation(matscipytest.MatSciPyTestCase):
+    """Class to store test for utils.py module."""
+
+    def test_validate_cubic_cell(self):
+            from ase.build import bulk as build_bulk
+            import warnings
+
+            warnings.simplefilter("always")
+            
+            bulk_fcc = build_bulk("Cu", crystalstructure="fcc", cubic=True)
+            bulk_bcc = build_bulk("W", crystalstructure="bcc", cubic=True)
+            bulk_dia = build_bulk("C", crystalstructure="diamond", cubic=True)
+
+            bulks = {
+                "fcc": bulk_fcc,
+                "bcc": bulk_bcc,
+                "diamond": bulk_dia
+            }
+
+            axes = np.eye(3).astype(int)
+
+            # Check legacy behaviour works
+            for structure in bulks.keys():
+                bulk = bulks[structure]
+                alat = bulk.cell[0, 0]
+
+                alat_valid, struct = utils_mod.validate_cubic_cell(alat, symbol="Cu", axes=axes, crystalstructure=structure)
+
+                assert alat_valid == alat
+                bpos = np.sort(bulk.get_scaled_positions(), axis=0)
+                spos = np.sort(struct.get_scaled_positions(), axis=0)
+                assert np.allclose(struct.cell[:, :], bulk.cell[:, :])
+                assert np.allclose(bpos, spos)
+
+            # Check new functionality accepted
+            for structure in bulks.keys():
+                bulk = bulks[structure]
+                alat = bulk.cell[0, 0]
+                alat_valid, struct = utils_mod.validate_cubic_cell(bulk, symbol="Cu", axes=axes, crystalstructure=structure)
+
+                assert alat_valid == alat
+                bpos = np.sort(bulk.get_scaled_positions(), axis=0)
+                spos = np.sort(struct.get_scaled_positions(), axis=0)
+                assert np.allclose(struct.cell[:, :], bulk.cell[:, :])
+                assert np.allclose(bpos, spos)
+
+            # Check Fail/Warn states
+
+            # "a" validation
+            self.assertRaises(TypeError, utils_mod.validate_cubic_cell, "Not a bulk or lattice constant", symbol="Cu", axes=axes, crystalstructure="fcc")
+            
+            # Crystalstructure validation
+            self.assertRaises(ValueError, utils_mod.validate_cubic_cell, bulk_fcc, symbol="Cu", axes=axes, crystalstructure="made up crystalstructure")
+            self.assertRaises(TypeError, utils_mod.validate_cubic_cell, bulk_fcc, symbol="Cu", axes=axes, crystalstructure=["not the ", "right datatype"])
+
+            # Bulk structure validation
+            self.assertWarns(UserWarning, utils_mod.validate_cubic_cell, bulk_fcc, symbol="Cu", axes=axes, crystalstructure="bcc")
+            self.assertWarns(UserWarning, utils_mod.validate_cubic_cell, bulk_fcc, symbol="Cu", axes=axes, crystalstructure="diamond")
+            self.assertWarns(UserWarning, utils_mod.validate_cubic_cell, bulk_bcc, symbol="Cu", axes=axes, crystalstructure="fcc")
+
+            ats = bulk_fcc.copy()
+            del ats[0]
+            # FCC with vacancy is not bulk FCC!
+            self.assertWarns(UserWarning, utils_mod.validate_cubic_cell, ats, symbol="Cu", axes=axes, crystalstructure="fcc")
+
+            
+            # Test w/ supercells of bulk
+            utils_mod.validate_cubic_cell(bulk_fcc * (2, 2, 2), symbol="Cu", axes=axes, crystalstructure="fcc")
+
+            # Test w/ more complex crystal
+            nacl = build_bulk("NaCl", crystalstructure="zincblende", a=5.8, cubic=True)
+
+            utils_mod.validate_cubic_cell(nacl, symbol="Cu", axes=axes, crystalstructure="diamond")
+
+            warnings.simplefilter("default")
+
+if __name__ == '__main__':
+    unittest.main()