Repository: libAtoms/matscipy (language: python)

Issue:
Support for celldisp in neighbor list
ASE's atoms object has a new celldisp parameters that specifies the origin of the cell. The neighbor list will break is it is present and not (0,0,0).


Discussion / hints:
Apparently celldisp is not the origin of the cell but only used for visualization purposes. We should probably defer implementing this until a proper cell origin is implemented in ASE.

See: https://gitlab.com/ase/ase/issues/197

Fix patch:
diff --git a/ChangeLog.md b/ChangeLog.md
index 4ad7a264..c37c7307 100644
--- a/ChangeLog.md
+++ b/ChangeLog.md
@@ -6,6 +6,8 @@ v0.9.0 (not yet released)
 
 - Scripts are now installable CLI
 - Updated documentation packages
+- Added tool for generating dislocation quadrupoles for any `CubicCrystalDislocation` or `CubicCrystalDissociatedDislocation`
+- Refactor of dislocation internals & testing to allow for easy integration of new dislocation types and displacement solver methods.
 
 v0.8.0 (12Apr23) 
 ----------------
diff --git a/docs/applications/cylinder_configurations.ipynb b/docs/applications/cylinder_configurations.ipynb
index 8de9fa09..933344ce 100644
--- a/docs/applications/cylinder_configurations.ipynb
+++ b/docs/applications/cylinder_configurations.ipynb
@@ -15,26 +15,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "application/vnd.jupyter.widget-view+json": {
-       "model_id": "928a58090bfa46b29cf864362bcbfdc9",
-       "version_major": 2,
-       "version_minor": 0
-      },
-      "text/plain": []
-     },
-     "metadata": {},
-     "output_type": "display_data"
-    }
-   ],
+   "outputs": [],
    "source": [
     "import numpy as np\n",
-    "import nglview # this import is necessary for rendering of the 3D view\n",
-    "from visualisation import show_dislocation"
+    "import nglview # this import is necessary for rendering of the 3D view"
    ]
   },
   {
@@ -77,29 +63,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "BCCScrew111Dislocation\n",
-      "BCCEdge111Dislocation\n",
-      "BCCMixed111Dislocation\n",
-      "BCCEdge100Dislocation\n",
-      "BCCEdge100110Dislocation\n",
-      "DiamondGlide30degreePartial\n",
-      "DiamondGlide90degreePartial\n",
-      "DiamondGlideScrew\n",
-      "DiamondGlide60Degree\n",
-      "FCCScrewShockleyPartial\n",
-      "FCCEdgeShockleyPartial\n",
-      "FCCScrew110Dislocation\n",
-      "FCCEdge110Dislocation\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "import matscipy.dislocation as disloc_mod\n",
     "import inspect\n",
@@ -124,123 +90,13 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": null,
    "metadata": {
     "tags": [
      "hide-output"
     ]
    },
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "      Step     Time          Energy         fmax\n",
-      "*Force-consistent energies used in optimization.\n",
-      "FIRE:    0 15:29:59      -17.792788*       1.5616\n",
-      "FIRE:    1 15:29:59      -17.796007*       1.2330\n",
-      "FIRE:    2 15:29:59      -17.798904*       0.6354\n",
-      "FIRE:    3 15:29:59      -17.799853*       0.2223\n",
-      "FIRE:    4 15:29:59      -17.799882*       0.1976\n",
-      "FIRE:    5 15:29:59      -17.799928*       0.1509\n",
-      "FIRE:    6 15:29:59      -17.799970*       0.0870\n",
-      "FIRE:    7 15:29:59      -17.799991*       0.0130\n",
-      "FIRE:    8 15:29:59      -17.799981*       0.0629\n",
-      "FIRE:    9 15:29:59      -17.799981*       0.0611\n",
-      "FIRE:   10 15:29:59      -17.799982*       0.0576\n",
-      "FIRE:   11 15:29:59      -17.799984*       0.0525\n",
-      "FIRE:   12 15:29:59      -17.799986*       0.0458\n",
-      "FIRE:   13 15:29:59      -17.799988*       0.0379\n",
-      "FIRE:   14 15:29:59      -17.799989*       0.0289\n",
-      "FIRE:   15 15:29:59      -17.799991*       0.0191\n",
-      "FIRE:   16 15:29:59      -17.799991*       0.0076\n",
-      "FIRE:   17 15:29:59      -17.799992*       0.0052\n",
-      "FIRE:   18 15:29:59      -17.799992*       0.0052\n",
-      "FIRE:   19 15:29:59      -17.799992*       0.0051\n",
-      "FIRE:   20 15:29:59      -17.799992*       0.0049\n",
-      "FIRE:   21 15:29:59      -17.799992*       0.0047\n",
-      "FIRE:   22 15:29:59      -17.799992*       0.0044\n",
-      "FIRE:   23 15:29:59      -17.799992*       0.0041\n",
-      "FIRE:   24 15:29:59      -17.799992*       0.0038\n",
-      "FIRE:   25 15:29:59      -17.799992*       0.0034\n",
-      "FIRE:   26 15:29:59      -17.799992*       0.0028\n",
-      "FIRE:   27 15:29:59      -17.799992*       0.0022\n",
-      "FIRE:   28 15:29:59      -17.799992*       0.0015\n",
-      "FIRE:   29 15:29:59      -17.799992*       0.0006\n",
-      "FIRE:   30 15:30:00      -17.799992*       0.0003\n",
-      "FIRE:   31 15:30:00      -17.799992*       0.0003\n",
-      "FIRE:   32 15:30:00      -17.799992*       0.0003\n",
-      "FIRE:   33 15:30:00      -17.799992*       0.0003\n",
-      "FIRE:   34 15:30:00      -17.799992*       0.0003\n",
-      "FIRE:   35 15:30:00      -17.799992*       0.0003\n",
-      "FIRE:   36 15:30:00      -17.799992*       0.0003\n",
-      "FIRE:   37 15:30:00      -17.799992*       0.0003\n",
-      "FIRE:   38 15:30:00      -17.799992*       0.0002\n",
-      "FIRE:   39 15:30:00      -17.799992*       0.0002\n",
-      "FIRE:   40 15:30:00      -17.799992*       0.0002\n",
-      "FIRE:   41 15:30:00      -17.799992*       0.0001\n",
-      "FIRE:   42 15:30:00      -17.799992*       0.0001\n",
-      "Fitting C_11\n",
-      "Strain array([-0.02, -0.01,  0.  ,  0.01,  0.02])\n",
-      "Stress array([-1.18143011e+01, -5.57729763e+00, -2.88261126e-04,  4.98286691e+00,\n",
-      "        9.85079614e+00]) GPa\n",
-      "Cij (gradient) / GPa    :     538.9035896451805\n",
-      "Error in Cij / GPa      :     16.584381718776033\n",
-      "Correlation coefficient :     0.9985824334709763\n",
-      "Setting C11 (1) to 3.363572 +/- 0.103512\n",
-      "\n",
-      "\n",
-      "Fitting C_21\n",
-      "Strain array([-0.02, -0.01,  0.  ,  0.01,  0.02])\n",
-      "Stress array([-5.18171556e+00, -2.29563243e+00, -2.88261126e-04,  1.85153957e+00,\n",
-      "        3.74132626e+00]) GPa\n",
-      "Cij (gradient) / GPa    :     219.93255649330877\n",
-      "Error in Cij / GPa      :     12.450223131451542\n",
-      "Correlation coefficient :     0.9952274688373499\n",
-      "Setting C21 (7) to 1.372711 +/- 0.077708\n",
-      "\n",
-      "\n",
-      "Fitting C_31\n",
-      "Strain array([-0.02, -0.01,  0.  ,  0.01,  0.02])\n",
-      "Stress array([-5.18171556e+00, -2.29563243e+00, -2.88261126e-04,  1.85153957e+00,\n",
-      "        3.74132626e+00]) GPa\n",
-      "Cij (gradient) / GPa    :     219.93255649330868\n",
-      "Error in Cij / GPa      :     12.450223131451247\n",
-      "Correlation coefficient :     0.9952274688373501\n",
-      "Updating C31 (7) with value 1.372711 +/- 0.077708\n",
-      "\n",
-      "\n",
-      "Fitting C_44\n",
-      "Strain array([-0.02, -0.01,  0.  ,  0.01,  0.02])\n",
-      "Stress array([-4.37426370e+00, -1.88654647e+00,  9.30635418e-15,  1.43563268e+00,\n",
-      "        2.90478138e+00]) GPa\n",
-      "Cij (gradient) / GPa    :     178.80269302570218\n",
-      "Error in Cij / GPa      :     12.702352202927047\n",
-      "Correlation coefficient :     0.9925146216927534\n",
-      "Setting C44 (4) to 1.115999 +/- 0.079282\n",
-      "\n",
-      "\n",
-      "[[b C11 b C12 b C12 b     b     b    ]\n",
-      " [b C12 b C11 b C12 b     b     b    ]\n",
-      " [b C12 b C12 b C11 b     b     b    ]\n",
-      " [b     b     b     b C44 b     b    ]\n",
-      " [b     b     b     b     b C44 b    ]\n",
-      " [b     b     b     b     b     b C44]]\n",
-      "\n",
-      " = \n",
-      "\n",
-      "[[538.9  219.93 219.93   0.     0.     0.  ]\n",
-      " [219.93 538.9  219.93   0.     0.     0.  ]\n",
-      " [219.93 219.93 538.9    0.     0.     0.  ]\n",
-      " [  0.     0.     0.   178.8    0.     0.  ]\n",
-      " [  0.     0.     0.     0.   178.8    0.  ]\n",
-      " [  0.     0.     0.     0.     0.   178.8 ]]\n",
-      "C_11 = 538.90 +/- 16.58 GPa\n",
-      "C_12 = 219.93 +/- 12.45 GPa\n",
-      "C_44 = 178.80 +/- 12.70 GPa\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "from matscipy.dislocation import get_elastic_constants\n",
     "# the calculator to provide forces and energies from the potential\n",
@@ -253,17 +109,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "3.143 (Angstrom), 538.90, 219.93, 178.80 (GPa)\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "print(f\"{alat:.3f} (Angstrom), {C11:.2f}, {C12:.2f}, {C44:.2f} (GPa)\")"
    ]
@@ -286,22 +134,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Cell orientation:\n",
-      "[[ 1  1 -2]\n",
-      " [-1  1  0]\n",
-      " [ 1  1  1]]\n",
-      "Burgers vector:\n",
-      "[1.57169453 1.57169453 1.57169453]\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "from matscipy.dislocation import BCCScrew111Dislocation\n",
     "\n",
@@ -323,26 +158,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 31,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "disloc SCF 0 |d1-d2|_inf = 0.001400723332078968\n",
-      "disloc SCF 1 |d1-d2|_inf = 6.704394990092455e-06\n",
-      "disloc SCF 2 |d1-d2|_inf = 3.2089362583809233e-08\n",
-      "\n",
-      "Cell vectors:\n",
-      "[[61.59759392  0.          0.        ]\n",
-      " [ 0.         62.235928    0.        ]\n",
-      " [ 0.          0.          2.72225477]]\n",
-      "\n",
-      "Burgers vector lenth: 2.72 Angstrom\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "W_screw_bulk, W_screw_dislo = W_screw.build_cylinder(radius=20)\n",
     "\n",
@@ -356,7 +174,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Now we can look at the created structure with the interactive view. Hovering your mouse over an atom will show chemical symbol of the atom and the structure identified by CNA algorithm. Dislocation is shown with an arrow corresponding to the dislocation line vector $\\vec l$. Hovering over the dislocation will display the dislocation name. Moreover, you can use your mouse to control the view:\n",
+    "Now we can use `CubicCrystalDislocation.view_cyl` to look at the created structure with an NGL interactive view. Hovering your mouse over an atom will show chemical symbol of the atom and the structure identified by CNA algorithm. Dislocation is shown with an arrow corresponding to the dislocation line vector $\\vec l$. Hovering over the dislocation will display the dislocation name. Moreover, you can use your mouse to control the view:\n",
     "\n",
     "- Translation: right click + drag\n",
     "- Rotation: left click + drag\n",
@@ -367,26 +185,11 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 32,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "application/vnd.jupyter.widget-view+json": {
-       "model_id": "df6d21351f824d99bd4ec008cc1417bc",
-       "version_major": 2,
-       "version_minor": 0
-      },
-      "text/plain": [
-       "NGLWidget()"
-      ]
-     },
-     "metadata": {},
-     "output_type": "display_data"
-    }
-   ],
+   "outputs": [],
    "source": [
-    "show_dislocation(W_screw_dislo, d_name=\"1/2<111> screw\")"
+    "W_screw.view_cyl(W_screw_dislo)"
    ]
   },
   {
@@ -402,27 +205,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 38,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "application/vnd.jupyter.widget-view+json": {
-       "model_id": "3645f686ef254210bd8007af69ca5ea2",
-       "version_major": 2,
-       "version_minor": 0
-      },
-      "text/plain": [
-       "NGLWidget()"
-      ]
-     },
-     "metadata": {},
-     "output_type": "display_data"
-    }
-   ],
+   "outputs": [],
    "source": [
     "longer_W_screw_dislo = W_screw_dislo * [1, 1, 10]\n",
-    "show_dislocation(longer_W_screw_dislo, d_name=\"1/2<111> screw\")"
+    "W_screw.view_cyl(longer_W_screw_dislo)"
    ]
   },
   {
@@ -438,22 +226,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Cell orientation:\n",
-      "[[ 1  1  1]\n",
-      " [ 1 -1  0]\n",
-      " [ 1  1 -2]]\n",
-      "Burgers vector:\n",
-      "[1.57169453 1.57169453 1.57169453]\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "from matscipy.dislocation import BCCEdge111Dislocation\n",
     "\n",
@@ -468,31 +243,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 34,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "disloc SCF 0 |d1-d2|_inf = 0.17503717189385382\n",
-      "disloc SCF 1 |d1-d2|_inf = 0.030669651063355258\n",
-      "disloc SCF 2 |d1-d2|_inf = 0.004376489640269432\n",
-      "disloc SCF 3 |d1-d2|_inf = 0.0008807112000284167\n",
-      "disloc SCF 4 |d1-d2|_inf = 0.00015452837692353505\n",
-      "disloc SCF 5 |d1-d2|_inf = 2.6112647029508107e-05\n",
-      "disloc SCF 6 |d1-d2|_inf = 5.7487804006051135e-06\n",
-      "disloc SCF 7 |d1-d2|_inf = 7.771108773702018e-07\n",
-      "\n",
-      "Cell vectors:\n",
-      "[[54.44509545  0.          0.        ]\n",
-      " [ 0.         53.34508114  0.        ]\n",
-      " [ 0.          0.          7.69969924]]\n",
-      "\n",
-      "Burgers vector lenth: 2.72 Angstrom\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "W_edge_bulk, W_edge_dislo = W_edge.build_cylinder(radius=15)\n",
     "\n",
@@ -511,26 +264,11 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 33,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "application/vnd.jupyter.widget-view+json": {
-       "model_id": "863a61925a3c4d8ca4de13786507f13c",
-       "version_major": 2,
-       "version_minor": 0
-      },
-      "text/plain": [
-       "NGLWidget()"
-      ]
-     },
-     "metadata": {},
-     "output_type": "display_data"
-    }
-   ],
+   "outputs": [],
    "source": [
-    "show_dislocation(W_edge_dislo, scale=0.25, d_name=\"1/2<111> edge\")"
+    "W_edge.view_cyl(W_edge_dislo, scale=0.25)"
    ]
   },
   {
@@ -544,22 +282,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Cell orientation:\n",
-      "[[ 1 -1 -2]\n",
-      " [ 1  1  0]\n",
-      " [ 1 -1  1]]\n",
-      "Burgers vector:\n",
-      "[ 1.57169453 -1.57169453 -1.57169453]\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "from matscipy.dislocation import BCCMixed111Dislocation\n",
     "\n",
@@ -574,45 +299,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 10,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "disloc SCF 0 |d1-d2|_inf = 0.10539187767041772\n",
-      "disloc SCF 1 |d1-d2|_inf = 0.018105716735450528\n",
-      "disloc SCF 2 |d1-d2|_inf = 0.002853744720166418\n",
-      "disloc SCF 3 |d1-d2|_inf = 0.00044617599789692486\n",
-      "disloc SCF 4 |d1-d2|_inf = 6.968641396387643e-05\n",
-      "disloc SCF 5 |d1-d2|_inf = 1.0882499856867955e-05\n",
-      "disloc SCF 6 |d1-d2|_inf = 1.6994203833053945e-06\n",
-      "disloc SCF 7 |d1-d2|_inf = 2.653822296916353e-07\n",
-      "\n",
-      "Cell vectors:\n",
-      "[[61.59759392  0.          0.        ]\n",
-      " [ 0.         62.235928    0.        ]\n",
-      " [ 0.          0.          2.72225477]]\n",
-      "\n",
-      "Burgers vector length: 2.72 Angstrom\n"
-     ]
-    },
-    {
-     "data": {
-      "application/vnd.jupyter.widget-view+json": {
-       "model_id": "c8e63ead26fc406590892d6e8897f0d0",
-       "version_major": 2,
-       "version_minor": 0
-      },
-      "text/plain": [
-       "NGLWidget()"
-      ]
-     },
-     "metadata": {},
-     "output_type": "display_data"
-    }
-   ],
+   "outputs": [],
    "source": [
     "W_mixed_bulk, W_mixed_dislo = W_mixed.build_cylinder(radius=20)\n",
     "\n",
@@ -621,7 +310,7 @@
     "\n",
     "print(f\"\\nBurgers vector length: {np.linalg.norm(W_mixed.burgers):.2f} Angstrom\")\n",
     "\n",
-    "show_dislocation(W_mixed_dislo, scale=0.5, d_name=\"1/2<111> mixed\")"
+    "W_mixed.view_cyl(W_mixed_dislo, scale=0.5)"
    ]
   },
   {
@@ -641,22 +330,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 11,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Cell orientation:\n",
-      "[[ 1  0  0]\n",
-      " [ 0  1  1]\n",
-      " [ 0 -1  1]]\n",
-      "Burgers vector:\n",
-      "[3.14338905 0.         0.        ]\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "from matscipy.dislocation import BCCEdge100110Dislocation\n",
     "\n",
@@ -671,43 +347,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 12,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "disloc SCF 0 |d1-d2|_inf = 0.256011923011219\n",
-      "disloc SCF 1 |d1-d2|_inf = 0.06016875789278883\n",
-      "disloc SCF 2 |d1-d2|_inf = 0.008688239875985548\n",
-      "disloc SCF 3 |d1-d2|_inf = 0.002036859690120879\n",
-      "disloc SCF 4 |d1-d2|_inf = 0.00025500641555943615\n",
-      "disloc SCF 5 |d1-d2|_inf = 8.302194247189476e-05\n",
-      "disloc SCF 6 |d1-d2|_inf = 1.6579146270331857e-05\n",
-      "disloc SCF 7 |d1-d2|_inf = 1.7272275882884713e-06\n",
-      "disloc SCF 8 |d1-d2|_inf = 6.753932630942927e-07\n"
-     ]
-    },
-    {
-     "data": {
-      "application/vnd.jupyter.widget-view+json": {
-       "model_id": "97a6a8aa8305419b802a254b534dcf62",
-       "version_major": 2,
-       "version_minor": 0
-      },
-      "text/plain": [
-       "NGLWidget()"
-      ]
-     },
-     "metadata": {},
-     "output_type": "display_data"
-    }
-   ],
+   "outputs": [],
    "source": [
     "W_100110_edge_bulk, W_100110_edge_dislo = W_100110_edge.build_cylinder(radius=17)\n",
-    "show_dislocation(W_100110_edge_dislo, scale=0.3, \n",
-    "                 d_name=\"<100>{110} edge\", d_color=[1, 0, 1])"
+    "W_100110_edge.view_cyl(W_100110_edge_dislo, scale=0.3, line_color=[1, 0, 1])"
    ]
   },
   {
@@ -721,22 +366,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Cell orientation:\n",
-      "[[ 1  0  0]\n",
-      " [ 0  0 -1]\n",
-      " [ 0  1  0]]\n",
-      "Burgers vector:\n",
-      "[3.14338905 0.         0.        ]\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "from matscipy.dislocation import BCCEdge100Dislocation\n",
     "\n",
@@ -751,42 +383,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 14,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "disloc SCF 0 |d1-d2|_inf = 0.29723896079013645\n",
-      "disloc SCF 1 |d1-d2|_inf = 0.03893895819760007\n",
-      "disloc SCF 2 |d1-d2|_inf = 0.011189961805871729\n",
-      "disloc SCF 3 |d1-d2|_inf = 0.001492901691747324\n",
-      "disloc SCF 4 |d1-d2|_inf = 0.00046542738716470744\n",
-      "disloc SCF 5 |d1-d2|_inf = 8.162136711258039e-05\n",
-      "disloc SCF 6 |d1-d2|_inf = 1.3593036747305831e-05\n",
-      "disloc SCF 7 |d1-d2|_inf = 4.370960947186475e-06\n",
-      "disloc SCF 8 |d1-d2|_inf = 5.007713447316431e-07\n"
-     ]
-    },
-    {
-     "data": {
-      "application/vnd.jupyter.widget-view+json": {
-       "model_id": "f3c5f9df4a574110902f9da11ff7372f",
-       "version_major": 2,
-       "version_minor": 0
-      },
-      "text/plain": [
-       "NGLWidget()"
-      ]
-     },
-     "metadata": {},
-     "output_type": "display_data"
-    }
-   ],
+   "outputs": [],
    "source": [
     "W_100_edge_bulk, W_100_edge_dislo = W_100_edge.build_cylinder(radius=20)\n",
-    "show_dislocation(W_100_edge_dislo, scale=0.4, d_name=\"<100>{001} edge\", d_color=[1, 0, 1])"
+    "W_100_edge.view_cyl(W_100_edge_dislo, scale=0.4, line_color=[1, 0, 1])"
    ]
   },
   {
@@ -800,81 +402,13 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 15,
+   "execution_count": null,
    "metadata": {
     "tags": [
      "hide-output"
     ]
    },
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "      Step     Time          Energy         fmax\n",
-      "*Force-consistent energies used in optimization.\n",
-      "FIRE:    0 15:30:04      -17.801046*       0.0001\n",
-      "Fitting C_11\n",
-      "Strain array([-0.02, -0.01,  0.  ,  0.01,  0.02])\n",
-      "Stress array([-5.12800014e+00, -2.50333320e+00, -1.41366167e-04,  2.37961959e+00,\n",
-      "        4.63379565e+00]) GPa\n",
-      "Cij (gradient) / GPa    :     244.0654437182867\n",
-      "Error in Cij / GPa      :     4.217968506243491\n",
-      "Correlation coefficient :     0.9995522932872427\n",
-      "Setting C11 (1) to 1.523337 +/- 0.026326\n",
-      "\n",
-      "\n",
-      "Fitting C_21\n",
-      "Strain array([-0.02, -0.01,  0.  ,  0.01,  0.02])\n",
-      "Stress array([-3.37539886e+00, -1.58459192e+00, -1.41366167e-04,  1.38619385e+00,\n",
-      "        2.58171226e+00]) GPa\n",
-      "Cij (gradient) / GPa    :     148.85008009740656\n",
-      "Error in Cij / GPa      :     6.776719382842348\n",
-      "Correlation coefficient :     0.9969053432709588\n",
-      "Setting C21 (7) to 0.929049 +/- 0.042297\n",
-      "\n",
-      "\n",
-      "Fitting C_31\n",
-      "Strain array([-0.02, -0.01,  0.  ,  0.01,  0.02])\n",
-      "Stress array([-3.37539886e+00, -1.58459192e+00, -1.41366167e-04,  1.38619385e+00,\n",
-      "        2.58171226e+00]) GPa\n",
-      "Cij (gradient) / GPa    :     148.850080097406\n",
-      "Error in Cij / GPa      :     6.776719382842323\n",
-      "Correlation coefficient :     0.9969053432709588\n",
-      "Updating C31 (7) with value 0.929049 +/- 0.042297\n",
-      "\n",
-      "\n",
-      "Fitting C_44\n",
-      "Strain array([-0.02, -0.01,  0.  ,  0.01,  0.02])\n",
-      "Stress array([-2.51031965e+00, -1.25545671e+00, -2.54901650e-17,  1.25652507e+00,\n",
-      "        2.51470655e+00]) GPa\n",
-      "Cij (gradient) / GPa    :     125.62034171207866\n",
-      "Error in Cij / GPa      :     0.0380968152678079\n",
-      "Correlation coefficient :     0.9999998620414655\n",
-      "Setting C44 (4) to 0.784061 +/- 0.000238\n",
-      "\n",
-      "\n",
-      "[[b C11 b C12 b C12 b     b     b    ]\n",
-      " [b C12 b C11 b C12 b     b     b    ]\n",
-      " [b C12 b C12 b C11 b     b     b    ]\n",
-      " [b     b     b     b C44 b     b    ]\n",
-      " [b     b     b     b     b C44 b    ]\n",
-      " [b     b     b     b     b     b C44]]\n",
-      "\n",
-      " = \n",
-      "\n",
-      "[[244.07 148.85 148.85   0.     0.     0.  ]\n",
-      " [148.85 244.07 148.85   0.     0.     0.  ]\n",
-      " [148.85 148.85 244.07   0.     0.     0.  ]\n",
-      " [  0.     0.     0.   125.62   0.     0.  ]\n",
-      " [  0.     0.     0.     0.   125.62   0.  ]\n",
-      " [  0.     0.     0.     0.     0.   125.62]]\n",
-      "C_11 = 244.07 +/- 4.22 GPa\n",
-      "C_12 = 148.85 +/- 6.78 GPa\n",
-      "C_44 = 125.62 +/- 0.04 GPa\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "eam_calc = EAM(\"../../tests/FeCuNi.eam.alloy\")\n",
     "\n",
@@ -884,17 +418,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "3.520 (Angstrom), 244.07, 148.85, 125.62 (GPa)\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "print(f\"{alat:.3f} (Angstrom), {C11:.2f}, {C12:.2f}, {C44:.2f} (GPa)\")"
    ]
@@ -908,22 +434,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 17,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Cell orientation:\n",
-      "[[ 1  1 -2]\n",
-      " [ 1  1  1]\n",
-      " [ 1 -1  0]]\n",
-      "Burgers vector:\n",
-      "[ 1.76 -1.76  0.  ]\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "from matscipy.dislocation import FCCScrew110Dislocation\n",
     "\n",
@@ -938,47 +451,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 18,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "disloc SCF 0 |d1-d2|_inf = 0.021791458261984697\n",
-      "disloc SCF 1 |d1-d2|_inf = 0.0002936041882874109\n",
-      "disloc SCF 2 |d1-d2|_inf = 4.8928735380371347e-05\n",
-      "disloc SCF 3 |d1-d2|_inf = 1.2696678054050815e-06\n",
-      "disloc SCF 4 |d1-d2|_inf = 1.0764587210498888e-07\n",
-      "disloc SCF 0 |d1-d2|_inf = 0.020321133865591423\n",
-      "disloc SCF 1 |d1-d2|_inf = 0.0006134137472214689\n",
-      "disloc SCF 2 |d1-d2|_inf = 5.9406489367308524e-05\n",
-      "disloc SCF 3 |d1-d2|_inf = 2.948341477265748e-06\n",
-      "disloc SCF 4 |d1-d2|_inf = 8.340714410803862e-08\n",
-      "\n",
-      "Cell vectors:\n",
-      "[[60.35542726  0.          0.        ]\n",
-      " [ 0.         60.96818843  0.        ]\n",
-      " [ 0.          0.          2.48901587]]\n",
-      "\n",
-      "Burgers vector length: 2.49 Angstrom\n"
-     ]
-    },
-    {
-     "data": {
-      "application/vnd.jupyter.widget-view+json": {
-       "model_id": "c4f9f14bc6c84ecebe6c05988046b9d6",
-       "version_major": 2,
-       "version_minor": 0
-      },
-      "text/plain": [
-       "NGLWidget()"
-      ]
-     },
-     "metadata": {},
-     "output_type": "display_data"
-    }
-   ],
+   "outputs": [],
    "source": [
     "Ni_screw_bulk, Ni_screw_dislo = Ni_screw.build_cylinder(radius=20)\n",
     "\n",
@@ -988,9 +463,7 @@
     "\n",
     "print(f\"\\nBurgers vector length: {np.linalg.norm(Ni_screw.burgers):.2f} Angstrom\")\n",
     "\n",
-    "show_dislocation(Ni_screw_dislo, \n",
-    "                 d_name=\"1/2<110> screw\", \n",
-    "                 d_color=[0, 0, 1])"
+    "Ni_screw.view_cyl(Ni_screw_dislo, line_color=[0, 0, 1])"
    ]
   },
   {
@@ -1008,48 +481,13 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 35,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "disloc SCF 0 |d1-d2|_inf = 0.02179145826198453\n",
-      "disloc SCF 1 |d1-d2|_inf = 0.0002936041882874109\n",
-      "disloc SCF 2 |d1-d2|_inf = 4.892873538059339e-05\n",
-      "disloc SCF 3 |d1-d2|_inf = 1.2696678051752307e-06\n",
-      "disloc SCF 4 |d1-d2|_inf = 1.0764587193845543e-07\n",
-      "disloc SCF 0 |d1-d2|_inf = 0.020321133865591506\n",
-      "disloc SCF 1 |d1-d2|_inf = 0.0006134137472214689\n",
-      "disloc SCF 2 |d1-d2|_inf = 5.9406489367308524e-05\n",
-      "disloc SCF 3 |d1-d2|_inf = 2.948341477515548e-06\n",
-      "disloc SCF 4 |d1-d2|_inf = 8.340714444804442e-08\n",
-      "\n",
-      "Expected partial distance: 10.8 Angstrom\n"
-     ]
-    },
-    {
-     "data": {
-      "application/vnd.jupyter.widget-view+json": {
-       "model_id": "edb77f22bec7476bb0cce5beb1723915",
-       "version_major": 2,
-       "version_minor": 0
-      },
-      "text/plain": [
-       "NGLWidget()"
-      ]
-     },
-     "metadata": {},
-     "output_type": "display_data"
-    }
-   ],
+   "outputs": [],
    "source": [
     "Ni_screw_bulk, Ni_screw_dislo = Ni_screw.build_cylinder(radius=20, partial_distance=5)\n",
     "print(f\"\\nExpected partial distance: {5 * Ni_screw.glide_distance:.1f} Angstrom\")\n",
-    "show_dislocation(Ni_screw_dislo, \n",
-    "                 partial_distance=5 * Ni_screw.glide_distance,\n",
-    "                 d_name=\"1/6<112> Shockley partial screw\")"
+    "Ni_screw.view_cyl(Ni_screw_dislo)"
    ]
   },
   {
@@ -1063,22 +501,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 20,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Cell orientation:\n",
-      "[[ 1 -1  0]\n",
-      " [ 1  1  1]\n",
-      " [-1 -1  2]]\n",
-      "Burgers vector:\n",
-      "[ 1.76 -1.76  0.  ]\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "from matscipy.dislocation import FCCEdge110Dislocation\n",
     "\n",
@@ -1093,47 +518,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 21,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "disloc SCF 0 |d1-d2|_inf = 0.04310715442453525\n",
-      "disloc SCF 1 |d1-d2|_inf = 0.0031359890399249857\n",
-      "disloc SCF 2 |d1-d2|_inf = 0.00018078113282300745\n",
-      "disloc SCF 3 |d1-d2|_inf = 1.4501349156037513e-05\n",
-      "disloc SCF 4 |d1-d2|_inf = 9.978781161934513e-07\n",
-      "disloc SCF 0 |d1-d2|_inf = 0.03462421627571943\n",
-      "disloc SCF 1 |d1-d2|_inf = 0.0022123322463467043\n",
-      "disloc SCF 2 |d1-d2|_inf = 0.00021564683900149317\n",
-      "disloc SCF 3 |d1-d2|_inf = 1.1838612354972411e-05\n",
-      "disloc SCF 4 |d1-d2|_inf = 6.324895580611667e-07\n",
-      "\n",
-      "Cell vectors:\n",
-      "[[64.71441261  0.          0.        ]\n",
-      " [ 0.         60.96818843  0.        ]\n",
-      " [ 0.          0.          4.31110195]]\n",
-      "\n",
-      "Burgers vector length: 2.49 Angstrom\n"
-     ]
-    },
-    {
-     "data": {
-      "application/vnd.jupyter.widget-view+json": {
-       "model_id": "1329680a67d3459ca901a4ccb7c66955",
-       "version_major": 2,
-       "version_minor": 0
-      },
-      "text/plain": [
-       "NGLWidget()"
-      ]
-     },
-     "metadata": {},
-     "output_type": "display_data"
-    }
-   ],
+   "outputs": [],
    "source": [
     "Ni_edge_bulk, Ni_edge_dislo = Ni_edge.build_cylinder(radius=20)\n",
     "\n",
@@ -1142,52 +529,18 @@
     "\n",
     "print(f\"\\nBurgers vector length: {np.linalg.norm(Ni_edge.burgers):.2f} Angstrom\")\n",
     "\n",
-    "show_dislocation(Ni_edge_dislo, d_name=\"1/2<110> edge\", d_color=[0, 0, 1])"
+    "Ni_edge.view_cyl(Ni_edge_dislo, line_color=[0, 0, 1])"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 22,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "disloc SCF 0 |d1-d2|_inf = 0.04310715442453525\n",
-      "disloc SCF 1 |d1-d2|_inf = 0.0031359890399249857\n",
-      "disloc SCF 2 |d1-d2|_inf = 0.00018078113282300745\n",
-      "disloc SCF 3 |d1-d2|_inf = 1.4501349155759957e-05\n",
-      "disloc SCF 4 |d1-d2|_inf = 9.978781162003902e-07\n",
-      "disloc SCF 0 |d1-d2|_inf = 0.03462421627571943\n",
-      "disloc SCF 1 |d1-d2|_inf = 0.002212332246347065\n",
-      "disloc SCF 2 |d1-d2|_inf = 0.00021564683900149317\n",
-      "disloc SCF 3 |d1-d2|_inf = 1.1838612354972411e-05\n",
-      "disloc SCF 4 |d1-d2|_inf = 6.324895580611667e-07\n",
-      "Expected partial distance: 12.4 Angstrom\n"
-     ]
-    },
-    {
-     "data": {
-      "application/vnd.jupyter.widget-view+json": {
-       "model_id": "cbffabd9c6eb418f91e58d021cd4912e",
-       "version_major": 2,
-       "version_minor": 0
-      },
-      "text/plain": [
-       "NGLWidget()"
-      ]
-     },
-     "metadata": {},
-     "output_type": "display_data"
-    }
-   ],
+   "outputs": [],
    "source": [
     "Ni_edge_bulk, Ni_edge_dislo = Ni_edge.build_cylinder(radius=20, partial_distance=10)\n",
     "print(f\"\\nExpected partial distance: {10 * Ni_edge.glide_distance:.1f} Angstrom\")\n",
-    "show_dislocation(Ni_edge_dislo, \n",
-    "                 partial_distance=10 * Ni_edge.glide_distance, \n",
-    "                 d_name=\"1/2<110> Shockley partial edge\")"
+    "Ni_edge.view_cyl(Ni_edge_dislo)"
    ]
   },
   {
@@ -1201,113 +554,13 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 23,
+   "execution_count": null,
    "metadata": {
     "tags": [
      "hide-output"
     ]
    },
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "      Step     Time          Energy         fmax\n",
-      "*Force-consistent energies used in optimization.\n",
-      "FIRE:    0 15:30:06      -34.692786*       0.0922\n",
-      "FIRE:    1 15:30:06      -34.692742*       0.1881\n",
-      "FIRE:    2 15:30:06      -34.692797*       0.0452\n"
-     ]
-    },
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "FIRE:    3 15:30:06      -34.692771*       0.1323\n",
-      "FIRE:    4 15:30:06      -34.692781*       0.1071\n",
-      "FIRE:    5 15:30:06      -34.692794*       0.0616\n",
-      "FIRE:    6 15:30:06      -34.692800*       0.0043\n",
-      "FIRE:    7 15:30:06      -34.692795*       0.0537\n",
-      "FIRE:    8 15:30:06      -34.692796*       0.0511\n",
-      "FIRE:    9 15:30:06      -34.692796*       0.0461\n",
-      "FIRE:   10 15:30:06      -34.692798*       0.0390\n",
-      "FIRE:   11 15:30:06      -34.692799*       0.0299\n",
-      "FIRE:   12 15:30:06      -34.692799*       0.0195\n",
-      "FIRE:   13 15:30:06      -34.692800*       0.0081\n",
-      "FIRE:   14 15:30:06      -34.692800*       0.0036\n",
-      "FIRE:   15 15:30:06      -34.692800*       0.0036\n",
-      "FIRE:   16 15:30:06      -34.692800*       0.0035\n",
-      "FIRE:   17 15:30:06      -34.692800*       0.0034\n",
-      "FIRE:   18 15:30:06      -34.692800*       0.0032\n",
-      "FIRE:   19 15:30:06      -34.692800*       0.0030\n",
-      "FIRE:   20 15:30:06      -34.692800*       0.0028\n",
-      "FIRE:   21 15:30:06      -34.692800*       0.0025\n",
-      "FIRE:   22 15:30:06      -34.692800*       0.0022\n",
-      "FIRE:   23 15:30:06      -34.692800*       0.0017\n",
-      "FIRE:   24 15:30:06      -34.692800*       0.0013\n",
-      "FIRE:   25 15:30:06      -34.692800*       0.0007\n",
-      "FIRE:   26 15:30:06      -34.692800*       0.0000\n",
-      "Fitting C_11\n",
-      "Strain array([-0.02, -0.01,  0.  ,  0.01,  0.02])\n",
-      "Stress array([-4.18437928e+00, -2.05313584e+00, -2.57946679e-05,  1.97559670e+00,\n",
-      "        3.87433910e+00]) GPa\n",
-      "Cij (gradient) / GPa    :     201.46169296718193\n",
-      "Error in Cij / GPa      :     2.6470882834179283\n",
-      "Correlation coefficient :     0.999741134311275\n",
-      "Setting C11 (1) to 1.257425 +/- 0.016522\n",
-      "\n",
-      "\n",
-      "Fitting C_21\n",
-      "Strain array([-0.02, -0.01,  0.  ,  0.01,  0.02])\n",
-      "Stress array([-1.13686951e+00, -5.40672589e-01, -2.57946679e-05,  4.89041361e-01,\n",
-      "        9.30285346e-01]) GPa\n",
-      "Cij (gradient) / GPa    :     51.640236590323106\n",
-      "Error in Cij / GPa      :     1.7644327851076396\n",
-      "Correlation coefficient :     0.9982534295801976\n",
-      "Setting C21 (7) to 0.322313 +/- 0.011013\n",
-      "\n",
-      "\n",
-      "Fitting C_31\n",
-      "Strain array([-0.02, -0.01,  0.  ,  0.01,  0.02])\n",
-      "Stress array([-1.13686951e+00, -5.40672589e-01, -2.57946679e-05,  4.89041361e-01,\n",
-      "        9.30285346e-01]) GPa\n",
-      "Cij (gradient) / GPa    :     51.64023659032299\n",
-      "Error in Cij / GPa      :     1.764432785107523\n",
-      "Correlation coefficient :     0.9982534295801978\n",
-      "Updating C31 (7) with value 0.322313 +/- 0.011013\n",
-      "\n",
-      "\n",
-      "Fitting C_44\n",
-      "Strain array([-0.02, -0.01,  0.  ,  0.01,  0.02])\n",
-      "Stress array([-2.43911684e+00, -1.19995048e+00,  2.16046225e-14,  1.16236368e+00,\n",
-      "        2.28858547e+00]) GPa\n",
-      "Cij (gradient) / GPa    :     118.17718786807852\n",
-      "Error in Cij / GPa      :     1.2857535435707688\n",
-      "Correlation coefficient :     0.9998224896147394\n",
-      "Setting C44 (4) to 0.737604 +/- 0.008025\n",
-      "\n",
-      "\n",
-      "[[b C11 b C12 b C12 b     b     b    ]\n",
-      " [b C12 b C11 b C12 b     b     b    ]\n",
-      " [b C12 b C12 b C11 b     b     b    ]\n",
-      " [b     b     b     b C44 b     b    ]\n",
-      " [b     b     b     b     b C44 b    ]\n",
-      " [b     b     b     b     b     b C44]]\n",
-      "\n",
-      " = \n",
-      "\n",
-      "[[201.46  51.64  51.64   0.     0.     0.  ]\n",
-      " [ 51.64 201.46  51.64   0.     0.     0.  ]\n",
-      " [ 51.64  51.64 201.46   0.     0.     0.  ]\n",
-      " [  0.     0.     0.   118.18   0.     0.  ]\n",
-      " [  0.     0.     0.     0.   118.18   0.  ]\n",
-      " [  0.     0.     0.     0.     0.   118.18]]\n",
-      "C_11 = 201.46 +/- 2.65 GPa\n",
-      "C_12 = 51.64 +/- 1.76 GPa\n",
-      "C_44 = 118.18 +/- 1.29 GPa\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "from matscipy.calculators.manybody.explicit_forms.stillinger_weber import StillingerWeber,\\\n",
     "                                                                Holland_Marder_PRL_80_746_Si\n",
@@ -1320,39 +573,18 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 24,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "5.431 (Angstrom), 201.46, 51.64, 118.18 (GPa)\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "print(f\"{alat:.3f} (Angstrom), {C11:.2f}, {C12:.2f}, {C44:.2f} (GPa)\")"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 25,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Cell orientation:\n",
-      "[[ 1  1 -2]\n",
-      " [ 1  1  1]\n",
-      " [ 1 -1  0]]\n",
-      "Burgers vector:\n",
-      "[ 2.71547466 -2.71547466  0.        ]\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "from matscipy.dislocation import DiamondGlideScrew\n",
     "\n",
@@ -1374,37 +606,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 36,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "\n",
-      "Cell vectors:\n",
-      "[[66.51527324  0.          0.        ]\n",
-      " [ 0.         75.25344122  0.        ]\n",
-      " [ 0.          0.          3.84026109]]\n",
-      "\n",
-      "Burgers vector length: 3.84 Angstrom\n"
-     ]
-    },
-    {
-     "data": {
-      "application/vnd.jupyter.widget-view+json": {
-       "model_id": "f7355f892a57430e9b62045a228e9ac3",
-       "version_major": 2,
-       "version_minor": 0
-      },
-      "text/plain": [
-       "NGLWidget()"
-      ]
-     },
-     "metadata": {},
-     "output_type": "display_data"
-    }
-   ],
+   "outputs": [],
    "source": [
     "Si_screw_bulk, Si_screw_dislo = Si_screw.build_cylinder(radius=20)\n",
     "\n",
@@ -1413,12 +617,9 @@
     "\n",
     "print(f\"\\nBurgers vector length: {np.linalg.norm(Si_screw.burgers):.2f} Angstrom\")\n",
     "\n",
-    "show_dislocation(Si_screw_dislo, \n",
-    "                 diamond_structure=True, \n",
-    "                 scale=0.3, \n",
+    "Si_screw.view_cyl(Si_screw_dislo, scale=0.3, line_color=[0, 0, 1],\n",
     "                 add_bonds=True, # bonds make it a bit easier to see the structure\n",
-    "                 d_name=\"1/2<110> screw\", \n",
-    "                 d_color=[0, 0, 1])"
+    "                 )"
    ]
   },
   {
@@ -1432,40 +633,14 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 27,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Expected partial distance: 16.6 Angstrom\n"
-     ]
-    },
-    {
-     "data": {
-      "application/vnd.jupyter.widget-view+json": {
-       "model_id": "ee85e17cbfb74e9fab64ef50de7c4517",
-       "version_major": 2,
-       "version_minor": 0
-      },
-      "text/plain": [
-       "NGLWidget()"
-      ]
-     },
-     "metadata": {},
-     "output_type": "display_data"
-    }
-   ],
+   "outputs": [],
    "source": [
     "Si_screw_bulk, Si_screw_dislo = Si_screw.build_cylinder(radius=20, partial_distance=5)\n",
     "\n",
     "print(f\"Expected partial distance: {5 * Si_screw.glide_distance:.1f} Angstrom\")\n",
-    "show_dislocation(Si_screw_dislo, \n",
-    "                 diamond_structure=True, # bonds make it a bit easier to see the structure\n",
-    "                 scale=0.3, add_bonds=True, \n",
-    "                 partial_distance=5 * Si_screw.glide_distance, \n",
-    "                 d_name=\"1/6<112> 30 degree partial screw\") "
+    "Si_screw.view_cyl(Si_screw_dislo, scale=0.3, add_bonds=True, line_color=[1, 0, 1]) "
    ]
   },
   {
@@ -1481,22 +656,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 28,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Cell orientation:\n",
-      "[[ 1  1 -2]\n",
-      " [ 1  1  1]\n",
-      " [ 1 -1  0]]\n",
-      "Burgers vector:\n",
-      "[ 2.71547466 -2.71547466  0.        ]\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "from matscipy.dislocation import DiamondGlide60Degree\n",
     "\n",
@@ -1511,37 +673,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 29,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "\n",
-      "Cell vectors:\n",
-      "[[66.51527324  0.          0.        ]\n",
-      " [ 0.         75.25344122  0.        ]\n",
-      " [ 0.          0.          3.84026109]]\n",
-      "\n",
-      "Burgers vector length: 3.84 Angstrom\n"
-     ]
-    },
-    {
-     "data": {
-      "application/vnd.jupyter.widget-view+json": {
-       "model_id": "8d7eb594109f426e99d12a4a0d60125b",
-       "version_major": 2,
-       "version_minor": 0
-      },
-      "text/plain": [
-       "NGLWidget()"
-      ]
-     },
-     "metadata": {},
-     "output_type": "display_data"
-    }
-   ],
+   "outputs": [],
    "source": [
     "Si_60_degree_screw_bulk, \\\n",
     "Si_60_degree_screw_dislo = Si_60_degree_screw.build_cylinder(radius=20)\n",
@@ -1551,62 +685,24 @@
     "\n",
     "print(f\"\\nBurgers vector length: {np.linalg.norm(Si_60_degree_screw.burgers):.2f} Angstrom\")\n",
     "\n",
-    "show_dislocation(Si_60_degree_screw_dislo, \n",
-    "                 diamond_structure=True, \n",
-    "                 scale=0.3, add_bonds=True,\n",
-    "                 d_name=\"1/2<110> 60 degree screw\", \n",
-    "                 d_color=[0, 0, 1])"
+    "Si_60_degree_screw.view_cyl(Si_60_degree_screw_dislo, scale=0.3, add_bonds=True,\n",
+    "                            line_color=[0, 0, 1])"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 30,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Expected partial distance: 16.6 Angstrom\n"
-     ]
-    },
-    {
-     "data": {
-      "application/vnd.jupyter.widget-view+json": {
-       "model_id": "8844142c06074a1ea1d1ceae83239da3",
-       "version_major": 2,
-       "version_minor": 0
-      },
-      "text/plain": [
-       "NGLWidget()"
-      ]
-     },
-     "metadata": {},
-     "output_type": "display_data"
-    }
-   ],
+   "outputs": [],
    "source": [
     "Si_60_degree_screw_bulk, \\\n",
     "Si_60_degree_screw_dislo = Si_60_degree_screw.build_cylinder(radius=20, \n",
     "                                                             partial_distance=5)\n",
     "\n",
     "print(f\"Expected partial distance: {5 * Si_60_degree_screw.glide_distance:.1f} Angstrom\")\n",
-    "show_dislocation(Si_60_degree_screw_dislo, \n",
-    "                 diamond_structure=True, \n",
-    "                 scale=0.3, \n",
-    "                 add_bonds=True, \n",
-    "                 d_color=[1, 0, 1],\n",
-    "                 partial_distance=5 * Si_60_degree_screw.glide_distance, \n",
-    "                 d_name=[\"1/6<112> 90 degree partial screw\", \n",
-    "                         \"1/6<112> 30 degree partial screw\"])"
+    "Si_60_degree_screw.view_cyl(Si_60_degree_screw_dislo, scale=0.3, add_bonds=True, \n",
+    "                                line_color=[[1, 0, 1], [0, 1, 1]])"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
@@ -1628,7 +724,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.15"
+   "version": "3.10.12"
   }
  },
  "nbformat": 4,
diff --git a/docs/applications/gamma_surface_advanced.ipynb b/docs/applications/gamma_surface_advanced.ipynb
index 80b080bf..1e897fef 100644
--- a/docs/applications/gamma_surface_advanced.ipynb
+++ b/docs/applications/gamma_surface_advanced.ipynb
@@ -26,7 +26,6 @@
     "from matscipy.calculators.manybody.explicit_forms.tersoff_brenner import \\\n",
     "                                   TersoffBrenner, Brenner_PRB_42_9458_C_I\n",
     "from matscipy.calculators.manybody import Manybody\n",
-    "from visualisation import show_dislocation\n",
     "\n",
     "calc = Manybody(**TersoffBrenner(Brenner_PRB_42_9458_C_I))\n",
     "\n",
@@ -37,11 +36,7 @@
     "\n",
     "C_screw_bulk, C_screw_dislo = C_screw.build_cylinder(radius=25, partial_distance=20)\n",
     "\n",
-    "view = show_dislocation(C_screw_dislo, \n",
-    "                 diamond_structure=True, \n",
-    "                 scale=0.6, \n",
-    "                 partial_distance=20 * C_screw.glide_distance, \n",
-    "                 d_name=\"1/6<112> 30 degree partial screw\") \n",
+    "view = C_screw.view_cyl(C_screw_dislo, scale=0.6) \n",
     "\n",
     "view.control.spin([0, 1, 0], np.math.pi) # flip along y-axis to align with SF visualisation\n",
     "view.control.zoom(0.7)\n",
@@ -138,7 +133,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.15"
+   "version": "3.10.12"
   }
  },
  "nbformat": 4,
diff --git a/docs/applications/multispecies_dislocations.ipynb b/docs/applications/multispecies_dislocations.ipynb
index 68fc4ab7..7236e480 100644
--- a/docs/applications/multispecies_dislocations.ipynb
+++ b/docs/applications/multispecies_dislocations.ipynb
@@ -22,7 +22,7 @@
     "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 visualisation import interactive_view\n",
     "from ase.visualize import view\n",
     "\n",
     "# Data from https://doi.org/10.1080/08927022.2011.602975\n",
@@ -71,8 +71,7 @@
     "\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)"
+    "disloc.view_cyl(GaAs_dislocation, scale=0.35, CNA_color=False, add_bonds=True)"
    ]
   },
   {
@@ -102,9 +101,8 @@
     "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, \n",
-    "                        d_name=\"GaAs 1/6<112> Alpha-90 degree partial\", \n",
-    "                        diamond_structure=True, CNA_color=False, add_bonds=True)\n",
+    "view = disloc.view_cyl(GaAs_dislocation, scale=0.35, CNA_color=False, add_bonds=True,\n",
+    "                        disloc_names=[\"GaAs 1/6<112> Alpha-90 degree partial\"])\n",
     "\n",
     "view.control.zoom(0.8)\n",
     "view"
@@ -128,9 +126,8 @@
     "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, \n",
-    "                        d_name=\"GaAs 1/6<112> Beta-90 degree partial\", \n",
-    "                        diamond_structure=True, CNA_color=False, add_bonds=True)\n",
+    "view = disloc.view_cyl(GaAs_dislocation, scale=0.35, CNA_color=False, add_bonds=True,\n",
+    "                        disloc_names=[\"GaAs 1/6<112> Beta-90 degree partial\"])\n",
     "\n",
     "view.control.zoom(0.8)\n",
     "view"
@@ -144,16 +141,15 @@
     "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",
+    "3. The displacement field applied to the bulk structure does not depend on atomic species - lattice distortions caused by the disorder will not 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",
+    "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 chemical 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",
+    "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}$. This can be done by passing the 6x6 Voigt elasticity matrix to the `C=` argument, rather than passing just `C11`, `C12`, and `C44` as was done previously.   \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."
+    "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 built by the dislocaion classes."
    ]
   },
   {
@@ -167,20 +163,30 @@
     "\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",
+    "# Reproduce data from Table 2 for In0.5 Ga0.5 As\n",
+    "C11, C22, C33, C12, C13, C23, C44, C55, C66 = \\\n",
+    "    120.31, 120.35, 121.12, 55.87, 57.24, 57.29, 58.26, 58.33, 58.25\n",
+    "\n",
+    "# Construct 6x6 elasticity tensor\n",
+    "C = np.array([\n",
+    "    [C11, C12, C13, 0,   0,   0],\n",
+    "    [C12, C22, C23, 0,   0,   0],\n",
+    "    [C13, C23, C33, 0,   0,   0],\n",
+    "    [0,   0,   0,   C44, 0,   0],\n",
+    "    [0,   0,   0,   0,   C55, 0],\n",
+    "    [0,   0,   0,   0,   0,   C66]\n",
+    "])\n",
     "\n",
     "GaAs = bulk(\"GaAs\", crystalstructure=\"zincblende\", cubic=True, a=alat)\n",
     "\n",
     "\n",
-    "disloc = DiamondGlide90degreePartial(GaAs, C11, C12, C44)\n",
+    "disloc = DiamondGlide90degreePartial(GaAs, C=C)\n",
     "\n",
     "GaAs_bulk, GaAs_dislocation = disloc.build_cylinder(radius=3 * alat)\n",
     "\n",
-    "show_dislocation(GaAs_dislocation, scale=0.35, \n",
-    "                 d_name=\"GaAs 1/6<112> Beta-90 degree partial\", \n",
-    "                 diamond_structure=True, CNA_color=False, add_bonds=True)"
+    "disloc.view_cyl(GaAs_dislocation, scale=0.35, CNA_color=False, add_bonds=True,\n",
+    "                        disloc_names=[\"GaAs 1/6<112> Beta-90 degree partial\"])"
    ]
   },
   {
@@ -216,9 +222,8 @@
     "InGaAs_dislocation = GaAs_dislocation.copy()\n",
     "InGaAs_dislocation.set_chemical_symbols(species)\n",
     "\n",
-    "view = show_dislocation(InGaAs_dislocation, scale=0.35, \n",
-    "                        d_name=\"InGaAs 1/6<112> Beta-90 degree partial\", \n",
-    "                        diamond_structure=True, CNA_color=False, add_bonds=True)\n",
+    "view = disloc.view_cyl(GaAs_dislocation, scale=0.35, CNA_color=False, add_bonds=True,\n",
+    "                        disloc_names=[\"InGaAs 1/6<112> Alpha-90 degree partial\"])\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",
@@ -247,7 +252,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.15"
+   "version": "3.10.12"
   }
  },
  "nbformat": 4,
diff --git a/docs/applications/plasticity.rst b/docs/applications/plasticity.rst
index 4aef0a89..7055924e 100644
--- a/docs/applications/plasticity.rst
+++ b/docs/applications/plasticity.rst
@@ -21,5 +21,6 @@ Tutorials:
 
     cylinder_configurations.ipynb
     multispecies_dislocations.ipynb
+    quadrupole_dislocations.ipynb
     gamma_surfaces.ipynb
     gamma_surface_advanced.ipynb
\ No newline at end of file
diff --git a/docs/applications/quadrupole_dislocations.ipynb b/docs/applications/quadrupole_dislocations.ipynb
new file mode 100644
index 00000000..ce55f393
--- /dev/null
+++ b/docs/applications/quadrupole_dislocations.ipynb
@@ -0,0 +1,276 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Dislocation Quadrupoles\n",
+    "Dislocation Quadrupoles are structures designed with the goal of enabling the study of dislocation systems in periodic cells. In order for the system to be periodic, the cell must have a net Burgers vector of zero. To achieve this, we place two dislocations with Burgers vectors of equal magnitude and opposite sign (i.e. $+\\textbf{b}$ and $-\\textbf{b}$) in the structure.\n",
+    "\n",
+    "The code is designed to allow the generation of these structures in as small a structure is required for a given separation of the dislocations along the glide direction.\n",
+    "\n",
+    "To start, lets look at the $90^\\circ$ partial dislocation in Diamond, using a dislocation cylinder approach documented in earlier docs and the same potential by [D. Holland and M. Marder](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.80.746):"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from matscipy.dislocation import DiamondGlide90degreePartial, get_elastic_constants\n",
+    "from matscipy.calculators.manybody.explicit_forms.stillinger_weber import StillingerWeber,\\\n",
+    "                                                                Holland_Marder_PRL_80_746_Si\n",
+    "from matscipy.calculators.manybody import Manybody\n",
+    "calc = Manybody(**StillingerWeber(Holland_Marder_PRL_80_746_Si))\n",
+    "\n",
+    "alat, C11, C12, C44 = get_elastic_constants(calculator=calc, symbol=\"Si\", verbose=False)\n",
+    "\n",
+    "\n",
+    "Si_disloc = DiamondGlide90degreePartial(alat, C11, C12, C44, symbol=\"Si\")\n",
+    "\n",
+    "bulk, disloc = Si_disloc.build_cylinder(radius=20)\n",
+    "\n",
+    "Si_disloc.view_cyl(disloc, scale=0.3, add_bonds=True)\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To build dislocation quadrupoles, we can use the `Quadrupole` class. The `__init__()` method of the class takes a dislocation class (e.g. `DiamondGlide90degreePartial`), followed by args and/or kwargs used to initialise that dislocation class (i.e, `C11`, `C12`, `C44`, and either `alat` or a bulk structure). From here, there are many class methods which generate different kinds of structures. Let's first look at the basic quadrupole structure, using `build_quadrupole`"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from matscipy.dislocation import Quadrupole\n",
+    "\n",
+    "quad = Quadrupole(DiamondGlide90degreePartial, alat, C11, C12, C44, symbol=\"Si\")\n",
+    "\n",
+    "quad_bulk, quad_disloc = quad.build_quadrupole(glide_separation=4)\n",
+    "\n",
+    "view = quad.view_quad(quad_disloc, scale=0.3, add_bonds=True)\n",
+    "view.control.zoom(0.2)\n",
+    "view"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "`glide_separation` is the most important argument to `build_quadrupole`, as it controls how far apart (in terms of glide distances) the quadrupole cores are, and thus controls the size of the quadrupole cell. \n",
+    "\n",
+    "By default, the cell is constructed such that each dislocation core is surrounded by four cores of opposite sign, and the periodicity of the system tries to make the distances between the central core and it's opposites as equal as possible (i.e. that the cores form two sublattices with square packing)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Dislocation Glide in quadrupoles\n",
+    "As quadrupoles are extremely useful for modelling dislocations using plane-wave DFT, it can be convenient to be able to set up initial guesses for complex processes such as dislocation glide. In this scenario, we assume that the full infinite dislocation line glides in unison, ignoring the true \"kink\" process.\n",
+    "\n",
+    "We can use the function `build_glide_quadrupoles` to construct a set of images for this system, which can optionally model the glide of either the \"left\" ($+\\mathbf{b}$) or \"right\" ($-\\mathbf{b}$) dislocation cores, or both at once."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "num_images = 11\n",
+    "\n",
+    "glide_quads = quad.build_glide_quadrupoles(nims=num_images, \n",
+    "                                            glide_left=True, # Allow left dislocation glide\n",
+    "                                            glide_right=True, # Allow right dislocation glide\n",
+    "                                            glide_separation=6,\n",
+    "                                            verbose=False)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To visualise the glide structures, we will combine ase's `plot_atoms` to convert a structure to a matplotlib plot, and then use FuncAnimation to animate the glide structures: "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def animate_glide(images, diamond=True):\n",
+    "    from ase.visualize.plot import plot_atoms\n",
+    "    import matplotlib.pyplot as plt\n",
+    "    from matplotlib.animation import FuncAnimation\n",
+    "    from matscipy.utils import get_structure_types\n",
+    "    from visualisation import show_HTML\n",
+    "\n",
+    "    fig, ax = plt.subplots(figsize=(10, 15))\n",
+    "    ax.axis(\"off\")\n",
+    "\n",
+    "    # Add extra reps of start and end points for clarity\n",
+    "    anim_images = [images[0]] * 3 + images + [images[-1]] * 3\n",
+    "\n",
+    "    def plot_frame(framedata):\n",
+    "        ax.clear()\n",
+    "        # Plot an individual frame of the animation \n",
+    "        framenum, atoms = framedata\n",
+    "\n",
+    "        # get CNA colours to enhance plot\n",
+    "        atom_labels, struct_names, colors = get_structure_types(atoms, \n",
+    "                                                                diamond_structure=diamond)\n",
+    "        atom_colors = [colors[atom_label] for atom_label in atom_labels]\n",
+    "\n",
+    "        plot_atoms(atoms, ax=ax, colors=atom_colors)\n",
+    "\n",
+    "\n",
+    "    animation = FuncAnimation(fig, plot_frame, frames=enumerate(anim_images),\n",
+    "                                save_count=len(anim_images),\n",
+    "                                init_func=lambda: None,\n",
+    "                                interval=200)\n",
+    "    \n",
+    "    # Need to convert animation to HTML in order for it to be visible on the docs\n",
+    "    return show_HTML(animation)\n",
+    "\n",
+    "animate_glide(glide_quads)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Quadrupole Kink\n",
+    "Dislocation kink is commonly a more energetically accessible way for a dislocation to move in the glide direction. Rather than the full dislocation line gliding as one, the kink mechanism involves a small section of the dislocation line nucleating out by a glide vector, forming a pair of kinks. This pair of kinks can then diffuse along the direction of the dislocation line, advancing the line locally by a single glide vector.\n",
+    "\n",
+    "In the quadrupole cells, we model an infinite chain of single kinks, where the dislocation line advances by one glide direction every periodic repetition of the structure in the Z direction. This is distinct from the kink-pair mechanism, which would need two kinks (which migrate the dislocation line in opposite directions), however the single kink structures can be useful in making kink behaviour accessible to the small scales required by expensive methods such as DFT.\n",
+    "\n",
+    "To actually construct kink structures from a quadrupole object, we use the `build_kink_quadrupole()` method of `Quadrupole`. The main two important parameters to define are `invert_direction`, which changes the directionality of the kinks (some systems have distinct \"left\" kinks and \"right\" kinks), and `z_reps`, which both changes the size of the final structure, as well as the initial guess for the length of the kink.\n",
+    "\n",
+    ":::{note}\n",
+    "Because the kink advances the system by a single glide vector, and cells are typically multiple glide vectors long, the structure will not initially be periodic. To solve this, some atomic layers are removed (as controlled by `layer_tol`). In multispecies systems, this may disrupt the atomic composition of the end structure.\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "kink_quad = quad.build_kink_quadrupole(z_reps=5,\n",
+    "                                       glide_separation=4,\n",
+    "                                       verbose=False, layer_tol=1e-2)\n",
+    "\n",
+    "quad.view_quad(kink_quad, scale=0.6, add_bonds=False)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Quadrupoles in other systems\n",
+    "The code has many other features to the general purpose methods described above, which are intended to support some common choices and desires which arise in other dislocation systems. Below are a few curated examples of dislocation systems, and how the output can be modified by additional parameters.\n",
+    "\n",
+    "### BCC Screw in W\n",
+    "Taking the perfect screw dislocation example in BCC W (and using the same Embedded Atom Potential from [Marinica _et. al._ 2013 paper](http://dx.doi.org/10.1088/0953-8984/25/39/395502) (version EAM4) as was used in previous dislocation documentation), we can see that the default behaviour of the code is to produce what can be called an \"easy-hard\" quadrupole, which is to say that one dislocation is in the \"easy\" core (with 3 atoms in white), and the other dislocation is in the \"hard\" core (with 6 atoms in white)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from matscipy.calculators.eam import EAM\n",
+    "from matscipy.dislocation import BCCScrew111Dislocation\n",
+    "eam_calc = EAM(\"../../tests/w_eam4.fs\")\n",
+    "alat, C11, C12, C44 = get_elastic_constants(calculator=eam_calc, symbol=\"W\", verbose=False)\n",
+    "\n",
+    "quad = Quadrupole(BCCScrew111Dislocation, alat, C11, C12, C44, symbol=\"W\")\n",
+    "\n",
+    "W_bulk, W_quad = quad.build_quadrupole(glide_separation=8, verbose=False)\n",
+    "\n",
+    "quad.view_quad(W_quad)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Whilst this kind of structure contains a lot of information about both kinds of cores, typically the \"easy\" core is much lower in energy, and thus it may be desirable to model an \"easy-easy\" quadrupole. To do this, we can just add an offset to the left core (half the glide distance in this case), using the `left_offset` parameter. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "W_bulk, W_quad = quad.build_quadrupole(glide_separation=8, verbose=False,\n",
+    "                                       left_offset=np.array([quad.glide_distance/2, 0, 0]))\n",
+    "\n",
+    "quad.view_quad(W_quad)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Dissociation of Screw dislocations in FCC Ni Quadrupoles\n",
+    "Dissociation of dislocation cores is another phenomenon that me be desirable in a quadrupole structure, especially as quadrupoles of just the partials can't give much insight into the physics controlling when the partial cores are very close. \n",
+    "\n",
+    "We can build quadrupole structures containing these dissociated cores by passing the `partial_distance` argument to `build_quadrupole`, similarly to how we would when building a cylinder for that dissociated system."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from matscipy.dislocation import FCCScrew110Dislocation\n",
+    "\n",
+    "eam_calc = EAM(\"../../tests/FeCuNi.eam.alloy\")\n",
+    "\n",
+    "# the function accepts any ASE type of calculator\n",
+    "alat, C11, C12, C44 = get_elastic_constants(calculator=eam_calc, symbol=\"Ni\", verbose=False)\n",
+    "\n",
+    "quad = Quadrupole(FCCScrew110Dislocation, alat, C11, C12, C44, symbol=\"Ni\")\n",
+    "\n",
+    "bulk, Ni_quad = quad.build_quadrupole(glide_separation=14, partial_distance=4, \n",
+    "                                      disp_tol=1e-2) # Lower the tolerance for the periodic displacement convergence\n",
+    "\n",
+    "quad.view_quad(Ni_quad)"
+   ]
+  }
+ ],
+ "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/visualisation.py b/docs/applications/visualisation.py
index 769adc41..873585f1 100644
--- a/docs/applications/visualisation.py
+++ b/docs/applications/visualisation.py
@@ -5,9 +5,10 @@
 
 import numpy as np
 # interactive visualisation  inside the notebook with nglview
-from nglview import show_ase, ASEStructure
 import matplotlib.pyplot as plt
 from IPython.display import HTML
+from nglview import show_ase, ASEStructure
+
 def show_HTML(anim):
     '''
     Convert matplotlib.animation.FuncAnimation 
@@ -25,155 +26,32 @@ def show_HTML(anim):
         </div>
         '''
     plt.close(fig=plt.gcf())
-    return HTML(output_html)
-
-
-# Get the results of Ovito Common Neighbor Analysis 
-# https://www.ovito.org/docs/current/reference/pipelines/modifiers/common_neighbor_analysis.html
-# and Identify Diamond modifier
-# https://www.ovito.org/docs/current/reference/pipelines/modifiers/identify_diamond.html
-# for better visualisation of the dislocation core
-# it will be identified as "other" structure type
-def get_structure_types(structure, diamond_structure=False):
-    """Get the results of Common Neighbor Analysis and 
-        Identify Diamond modifiers from Ovito
-    Args:
-        structure (ase.atoms): input structure
-    Returns:
-        atom_labels (array of ints): per atom labels of the structure types
-        structure_names (list of strings): names of the structure types
-        colors (list of strings): colors of the structure types in hex format
-    """
-    from ovito.io.ase import ase_to_ovito
-    from ovito.modifiers import CommonNeighborAnalysisModifier, IdentifyDiamondModifier
-    from ovito.pipeline import StaticSource, Pipeline
-    ovito_structure = structure.copy()
-    if "fix_mask" in ovito_structure.arrays:
-        del ovito_structure.arrays["fix_mask"]
-    
-    if diamond_structure:
-        modifier = IdentifyDiamondModifier()
-    else:
-        modifier = CommonNeighborAnalysisModifier() 
-    
-    data = ase_to_ovito(ovito_structure)
-    pipeline = Pipeline(source=StaticSource(data=data))
-    pipeline.modifiers.append(modifier)
-    data = pipeline.compute()
-
-    atom_labels = data.particles['Structure Type'].array
-
-    structure_names = [structure.name for structure in modifier.structures]
-    colors = [structure.color for structure in modifier.structures]
-    hex_colors = ['#{:02x}{:02x}{:02x}'.format(int(r*255), int(g*255), int(b*255)) for r, g, b in colors] 
+    return HTML(output_html)    
 
-    return atom_labels, structure_names, hex_colors
-
-# 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";
+def interactive_view(system, scale=0.5, name=""):
+    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";
                         }
-                    } 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.
-    The cylinder is hollow by default so the second cylinder is needed to close it.
-    In case partial distance is provided, two dislocation lines are added and are shifter accordingly.
-    '''
+                    });
+                    """
     
-    center = system.positions.mean(axis=0)
-    view.shape.add_cylinder((center[0] + x_shift, center[1], -2.0), 
-                            (center[0] + x_shift, center[1], system.cell[2][2] - 0.5),
-                            color,
-                            [0.3],
-                            name)
-    
-    view.shape.add_cone((center[0] + x_shift, center[1], -2.0), 
-                        (center[0] + x_shift, center[1], 0.5),
-                        color,
-                        [0.3],
-                        name)
-       
-    view.shape.add_cone((center[0] + x_shift, center[1], system.cell[2][2] - 0.5), 
-                        (center[0] + x_shift, center[1], system.cell[2][2] + 1.0),
-                        color,
-                        [0.55],
-                        name)
-    
-
-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, 
-                                                               diamond_structure=diamond_structure)
-
-    view = show_ase(system)
-    view.hide([0])
-    
-    if add_bonds: # add bonds between all atoms to have bonds between structures
-        component = view.add_component(ASEStructure(system), default_representation=False, name='between structures')
-        component.add_ball_and_stick(cylinderOnly=True, radiusType='covalent', radiusScale=scale, aspectRatio=0.1)
-    
-    for structure_type in np.unique(atom_labels):
-        # every structure type is a different component
-        mask = atom_labels == structure_type
-        component = view.add_component(ASEStructure(system[mask]), 
-                                       default_representation=False, name=str(structure_names[structure_type]))
-        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:
-            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:
-        add_dislocation(view, system, d_name + ' dislocation line', d_color)
-    else:
-        if type(d_name) == list:
-            add_dislocation(view, system, d_name[0] + ' dislocation line', d_color, x_shift= -partial_distance / 2.0)
-            add_dislocation(view, system, d_name[1] + ' dislocation line', d_color, x_shift= partial_distance / 2.0)
-        else:
-            add_dislocation(view, system, d_name + ' dislocation line', d_color, x_shift= -partial_distance / 2.0)
-            add_dislocation(view, system, d_name + ' dislocation line', d_color, x_shift= partial_distance / 2.0)
-
-    view.camera = 'orthographic'
-    view.parameters = {"clipDist": 0}
-
-    #view._remote_call("setSize", target="Widget", args=["400px", "300px"])
-    view.layout.width = '100%'
-    view.layout.height = '300px'
-    view.center()
-
-    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), 
diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index 08ff46e6..067a8bab 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -27,6 +27,7 @@
 import numpy as np
 
 from abc import ABCMeta
+import warnings
 
 from scipy.optimize import minimize
 
@@ -35,7 +36,7 @@
 from ase.constraints import FixAtoms, StrainFilter
 from ase.optimize import FIRE
 from ase.optimize.precon import PreconLBFGS
-from ase.build import bulk
+from ase.build import bulk, stack
 from ase.calculators.lammpslib import LAMMPSlib
 from ase.units import GPa  # unit conversion
 from ase.io import read
@@ -43,8 +44,8 @@
 from matscipy.neighbours import neighbour_list, mic
 from matscipy.elasticity import fit_elastic_constants
 from matscipy.elasticity import Voigt_6x6_to_full_3x3x3x3
-from matscipy.elasticity import cubic_to_Voigt_6x6
-from matscipy.utils import validate_cubic_cell
+from matscipy.elasticity import cubic_to_Voigt_6x6, coalesce_elastic_constants
+from matscipy.utils import validate_cubic_cell, points_in_polygon2D
 
 
 def make_screw_cyl(alat, C11, C12, C44,
@@ -2068,15 +2069,15 @@ def ovito_dxa_straight_dislo_info(disloc, structure="BCC", replicate_z=3):
     if 'fix_mask' in dxa_disloc.arrays:
         del dxa_disloc.arrays['fix_mask']
 
-    input_crystal_structures = {"BCC": DislocationAnalysisModifier.Lattice.BCC,
-                                "FCC": DislocationAnalysisModifier.Lattice.FCC,
-                                "Diamond": DislocationAnalysisModifier.Lattice.CubicDiamond}
+    input_crystal_structures = {"bcc": DislocationAnalysisModifier.Lattice.BCC,
+                                "fcc": DislocationAnalysisModifier.Lattice.FCC,
+                                "diamond": DislocationAnalysisModifier.Lattice.CubicDiamond}
 
     data = ase_to_ovito(dxa_disloc)
     pipeline = Pipeline(source=StaticSource(data=data))
     pipeline.modifiers.append(ReplicateModifier(num_z=replicate_z))
     dxa = DislocationAnalysisModifier(
-          input_crystal_structure=input_crystal_structures[structure])
+          input_crystal_structure=input_crystal_structures[structure.lower()])
     pipeline.modifiers.append(dxa)
 
     data = pipeline.compute()
@@ -2164,11 +2165,11 @@ class AnisotropicDislocation:
     Theory of Dislocations, 2nd ed. Similar to class `CubicCrystalDislocation`.
     """
 
-    def __init__(self, C11, C12, C44, axes,
-                 slip_plane, disloc_line, burgers):
+    def __init__(self, axes, slip_plane, disloc_line, burgers,
+                 C11=None, C12=None, C44=None, C=None):
         """
         Setup a dislocation in a cubic crystal. 
-        C11, C12 and C44 are the elastic constants in the Catersian geometry.
+        C11, C12 and C44 are the elastic constants in the Cartesian geometry.
         The arg 'axes' is a 3x3 array containing axes of frame of reference (eg.crack system).
         The dislocation parameters required are the slip plane normal ('slip_plane'), the 
         dislocation line direction ('disloc_line') and the Burgers vector ('burgers').
@@ -2188,14 +2189,15 @@ def __init__(self, C11, C12, C44, axes,
         n = np.einsum('ij,j', A, n)
         xi = np.einsum('ij,j', A, xi)
         b = np.einsum('ij,j', A, burgers)
+
+        
         
         # define elastic constant in Voigt notation
-        C_v = cubic_to_Voigt_6x6(C11, C12, C44)
-        # convert Voigt to tensor
-        C = Voigt_6x6_to_full_3x3x3x3(C_v)
+        Cijkl = coalesce_elastic_constants(C11, C12, C44, C, convention="Cijkl")
+
         # rotate elastic matrix
         cijkl = np.einsum('ig,jh,ghmn,km,ln', \
-            A, A, C, A, A)
+            A, A, Cijkl, A, A)
     
         # solve the Stroh sextic formalism: the same as Stroh.py from atomman
         # this part can also be replaced by calling Stroh.py from atomman
@@ -2339,7 +2341,7 @@ class CubicCrystalDislocation(metaclass=ABCMeta):
     # (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")
+                 "C", "alat", "unit_cell", "name")
 
     # Attributes with defaults
     # These may be overridden by child classes
@@ -2350,23 +2352,42 @@ class CubicCrystalDislocation(metaclass=ABCMeta):
     stroh = None
     ADstroh = None
 
-    def __init__(self, a, C11, C12, C44, symbol="W"):
+    # Available solvers for atomic displacements
+    avail_methods = [
+        "atomman", "adsl"
+        ]
+    
+    # Attempt to load the atomman module once, 
+    try:
+        import atomman as atm
+        atomman = atm
+    except ImportError:
+        atomman = None
+
+    def __init__(self, a, C11=None, C12=None, C44=None, C=None, symbol="W"):
         """
         This class represents a dislocation in a cubic crystal
 
         The dislocation is defined by the crystal unit cell,
-        elastic constants C11, C12 and C44, crystal axes,
+        elastic constants, crystal axes,
         burgers vector and optional shift and parity vectors.
 
+        The crystal used as the base for constructing dislocations can be defined 
+        by a lattice constant and a chemical symbol. For multi-species systems, an ase 
+        Atoms object defining the cubic unit cell can instead be used.
+
+        Elastic constants can be defined either through defining C11, C12 and C44, or by passing 
+        the full 6x6 elasticity matrix
+
         Parameters
         ----------
         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
+        C11, C12, C44 : float
+            Elastic Constants (in GPa) for cubic symmetry
+        C : np.array
+            Full 6x6 elasticity matrix (in GPa)
         symbol : str
             Chemical symbol used to construct unit cell (if "a" is a lattice constant)
 
@@ -2400,37 +2421,112 @@ def __init__(self, a, C11, C12, C44, symbol="W"):
         # 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
+        # Skip these for dissociated dislocations, as these operate via classmethod properties
+        if not issubclass(self.__class__, CubicCrystalDissociatedDislocation):
+            self.axes = self.axes.copy()
+            self.unit_cell_core_position_dimensionless = self.unit_cell_core_position_dimensionless.copy()
+            self.parity = self.parity.copy()
         
-        self.alat, self.unit_cell = validate_cubic_cell(a, symbol, self.axes, self.crystalstructure, self.pbc)
+            self.alat, self.unit_cell = validate_cubic_cell(a, symbol, self.axes, self.crystalstructure, self.pbc)
+
+            # Empty dict for solver methods (e.g. AnisotropicDislocation.displacements)
+            self.solvers = {
+                m : None for m in self.avail_methods
+            }
+
+        # Sort out elasticity matrix into 6x6 convention (as we know the system is cubic)
+        self.C = coalesce_elastic_constants(C11, C12, C44, C, convention="Cij")
+
+    def init_solver(self, method="atomman"):
+        '''
+        Run the correct solver initialiser
+
+        Solver init should take only the self arg, and set self.solver[method] to point to the function 
+        to calculate displacements
+
+        e.g. self.solvers["atomman"] = self.stroh.displacement
+        '''
+        solver_initters = {
+            "atomman" : self.init_atomman,
+            "adsl" : self.init_adsl
+        }
+
+        # Execute init routine
+        solver_initters[method]()
+    
+    def get_solvers(self, method="atomman"):
+        '''
+        Get list of dislocation displacement solvers
+
+        As CubicCrystalDislocation models a single core, there is only one solver,
+        therefore solvers is a len(1) list
+
+        Parameters
+        ----------
+        use_atomman: bool
+            Use the Stroh solver included in atomman (requires atomman package) to
+            solve for displacements, or use the AnisotropicDislocation class
 
-        self.glide_distance = self.alat * self.glide_distance_dimensionless
+        Returns
+        -------
+        solvers: list
+            List of functions f(positions) = displacements
+        '''
+        has_atomman = False if self.atomman is None else True
+
+        if type(method) != str:
+            if has_atomman:
+                warnings.warn("non-string method arg, falling back to method='atomman'")
+                method = "atomman"
+            else:
+                warnings.warn("non-string method arg, falling back to method='adsl'")
+                method = "adsl"
+
+        if method.lower() not in self.avail_methods:
+            # Passed an invalid method
+            if has_atomman:
+                warnings.warn(f"method='{method}' not one of {self.avail_methods}, falling back to method='atomman'")
+                method = "atomman"
+            else:
+                warnings.warn(f"method='{method}' not one of {self.avail_methods}, falling back to method='adsl'")
+                method = "adsl"
+
+        if method.lower() == "atomman" and not has_atomman:
+            # Atomman import failed
+            warnings.warn("Import of atomman failed, falling back to method='adsl'")
+            method="adsl"
+
+        if self.solvers[method] is None:
+            # Method needs initialising
+            self.init_solver(method)
+        
+        return [self.solvers[method]]
 
-    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 init_atomman(self):
+        '''
+        Init atomman stroh solution solver (requires atomman)
+        '''
+        c = self.atomman.ElasticConstants(Cij=self.C)
+        self.stroh = self.atomman.defect.Stroh(c, burgers=self.burgers, axes=self.axes)
 
-    def init_anisotropic_dislocation(self):
+        self.solvers["atomman"] = self.stroh.displacement
 
+    def init_adsl(self):
+        '''
+        Init adsl (Anisotropic DiSLocation) solver
+        '''
         # Extract consistent parameters
         axes = self.axes
         slip_plane = axes[1].copy() 
         disloc_line = axes[2].copy() 
 
         # Setup AnistoropicDislocation object
-        self.ADstroh = AnisotropicDislocation(self.C11, self.C12, self.C44, 
-                                           axes, slip_plane, disloc_line, 
-                                           self.burgers)
-
+        self.ADstroh = AnisotropicDislocation(axes, slip_plane, disloc_line, 
+                                           self.burgers, C=self.C)
+        
+        self.solvers["adsl"] = self.ADstroh.displacement
 
     # @property & @var.setter decorators used to ensure var and var_dimensionless don't get out of sync
     @property
@@ -2440,12 +2536,16 @@ def burgers(self):
     @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
 
+    def invert_burgers(self):
+        '''
+        Modify dislocation to produce same dislocation with opposite burgers vector
+        '''
+        self.burgers_dimensionless *= -1
+
     @property
     def unit_cell_core_position(self):
         return self.unit_cell_core_position_dimensionless * self.alat
@@ -2462,6 +2562,140 @@ def glide_distance(self):
     def glide_distance(self, distance):
         self.glide_distance_dimensionless = distance / self.alat
 
+    def _build_supercell(self, targetx, targety):
+        '''
+        Build supercell in 2D from self.unit_cell, with target dimensions (targetx, targety), specified in Angstrom.
+        Supercell is constructed to have cell lengths >= the target value (i.e. cell[0, 0] >= targetx)
+
+        Parameters
+        ----------
+        targetx: float
+            Target length (in Ang) of the cell in the "x" direction (i.e. cell[0, 0])
+        targety: float
+            Target length (in Ang) of the cell in the "y" direction (i.e. cell[1, 0])
+
+        
+        Returns a supercell of self.unit_cell with cell[0, 0] >= targetx and cell[1, 1] >= targety
+        '''
+
+        base_cell = self.unit_cell
+
+        xreps = np.ceil(targetx / (base_cell.cell[0, 0] - 1e-2)).astype(int)
+        yreps = np.ceil(targety / (base_cell.cell[1, 1] - 1e-2)).astype(int)
+
+        sup = base_cell.copy() * (xreps, yreps, 1)
+        return sup
+
+    def _build_bulk_cyl(self, radius, core_positions, fix_rad, extension,
+                        self_consistent, method, verbose):
+        '''
+        Build bulk cylinder config from args supplied by self.build_cylinder
+
+        Parameters
+        ----------
+        radius: float
+            Radius for a radial mask for cyl config
+        core_positions: np.array
+            Positions of all dislocation cores
+            Should be of shape (ncores, 3)
+            Can add fictitious core positions in order to extend cyl as needed
+            e.g. for glide configs, etc
+        fix_rad: float
+            Radius for fixed atoms constraints
+        extension: np.array
+            (N, 3) array of vector extensions from each core.
+            Will add a fictitious core at position 
+            core_positions[i, :] + extension[i, :], for the 
+            purposes of adding extra atoms
+
+        Returns
+        -------
+        cyl: ase.Atoms
+            Bulk cyl, cut down based on radius, complete
+            with FixAtoms constraints
+        '''
+        from matscipy.utils import radial_mask_from_polygon2D
+
+
+        if self_consistent is None:
+            self_consistent = self.self_consistent
+
+        if len(core_positions.shape) == 1:
+            core_positions = core_positions[np.newaxis, :]
+
+        if len(extension.shape) == 1:
+            extension = extension[np.newaxis, :]
+
+        if np.all(extension.shape != core_positions.shape):            
+            if extension.shape[0] < core_positions.shape[0] and \
+                extension.shape[1] == core_positions.shape[1]:
+                # Too few extensions passed, assume that the extensions are for the first N cores
+                # Pad extension with zeros
+                extension = np.vstack([extension, np.zeros((
+                    core_positions.shape[0] - extension.shape[0], core_positions.shape[1]
+                ))])
+            else:
+                raise ValueError(f"{extension.shape=} does not match {core_positions.shape=}")
+
+        # Only care about non-zero extensions
+        exts = extension + core_positions
+
+        xmax = np.max([np.max(core_positions[:, 0]), np.max(exts[:, 0])])
+        xmin = np.min([np.min(core_positions[:, 0]), np.min(exts[:, 0])])
+
+        ymax = np.max([np.max(core_positions[:, 1]), np.max(exts[:, 1])])
+        ymin = np.min([np.min(core_positions[:, 1]), np.min(exts[:, 1])])
+
+        xlen = xmax - xmin + 2*radius
+        ylen = ymax - ymin + 2*radius
+
+        sup = self._build_supercell(xlen, ylen)
+
+        center = np.diag(sup.cell) / 2
+
+        # Shift everything so that the dislocations are reasonably central
+
+        shift = center - 0.5 * (np.array([xmax, ymax, 0]) + np.array([xmin, ymin, 0]))
+
+        pos = sup.get_positions()
+        pos += shift
+        sup.set_positions(pos)
+        sup.wrap()
+
+        new_core_positions = core_positions.copy()
+
+        new_core_positions += shift
+
+        # Mask out extensions that are all zeros
+        nonzero_exts = np.sum(np.abs(extension), axis=-1) > 0.0
+
+        exts = exts[nonzero_exts, :]
+
+        mask_positions = np.vstack([new_core_positions, exts + shift])
+
+        idxs = np.argsort(mask_positions, axis=0)
+        mask_positions = np.take_along_axis(mask_positions, idxs, axis=0)
+        
+        mask = radial_mask_from_polygon2D(sup.get_positions(), mask_positions, radius, inner=True)
+
+        cyl = sup[mask]
+
+        # disloc is a copy of bulk with displacements applied
+        disloc = cyl.copy()
+
+        disloc.positions += self.displacements(cyl.positions, new_core_positions,
+                                               self_consistent=self_consistent, method=method,
+                                               verbose=verbose)
+
+        if fix_rad:
+            fix_mask = ~radial_mask_from_polygon2D(cyl.get_positions(), mask_positions, radius - fix_rad, inner=True)
+
+            fix_atoms = FixAtoms(mask=fix_mask)
+            disloc.set_constraint(fix_atoms)
+            disloc.arrays["fix_mask"] = fix_mask
+
+        return cyl, disloc, new_core_positions
+
     def plot_unit_cell(self, ms=250, ax=None):
         import matplotlib.pyplot as plt
 
@@ -2492,139 +2726,161 @@ def plot_unit_cell(self, ms=250, ax=None):
         ax.set_xlabel(r"$\AA$")
         ax.set_ylabel(r"$\AA$")
 
-    def self_consistent_displacements(self, solver, bulk_positions, center, 
-                                      tol=1e-6, max_iter=100, verbose=True):
+    def self_consistent_displacements(self, solvers, bulk_positions, core_positions, 
+                                      tol=1e-6, max_iter=100, verbose=True, mixing=0.8):
+        '''
+        Compute dislocation displacements self-consistently, with max_iter capping the number of iterations
+        Each dislocation core uses a separate solver, which computes the displacements associated with positions 
+        relative to that core (i.e. that the core is treated as lying at [0, 0, 0])
+
+        Parameters
+        ----------
+        solvers: list
+            List of functions, solvers[i] computes the displacement field for the dislocation
+            defined by core_positions[i, :]
+        bulk_positions: np.array
+            (N, 3) array of positions of atoms in the bulk configuration
+        core_positions: np.array
+            Positions of each dislocation core (shape = (ncores, 3))
+        tol: float
+            Displacement tolerance for the convergence of the self-consistent loop
+        max_iter: int
+            Maximum number of iterations of the self-consistent cycle to run before throwing 
+            an exception
+        verbose: bool
+            Enable/Disable printing progress of the self-consistent cycle each iteration
+
+        Returns
+        -------
+        displacements: np.array
+        shape (N, 3) array of displacements to apply to bulk_positions to form
+        the dislocation structure
+
+        Raises a RuntimeError if max_iter is reached without convergence
+        '''
+        if len(core_positions.shape) == 1:
+            core_positions = core_positions[np.newaxis, :]
+
+        ncores = core_positions.shape[0]
+
+        disp1 = np.zeros_like(bulk_positions)
+        for i in range(ncores):
+            disp1 += solvers[i](bulk_positions - core_positions[i, :]).real
         
-        disp1 = np.real(solver(bulk_positions - center))
         if max_iter == 0:
             return disp1
         
         # Self-consistent calculation
         res = np.inf
-        i = 0
-        while res > tol:
+        for i in range(max_iter):
             disloc_positions = bulk_positions + disp1
-            disp2 = np.real(solver(disloc_positions - center))
+            
+            disp2 = np.zeros_like(bulk_positions)
+            for j in range(ncores):
+                disp2 += solvers[j](disloc_positions - core_positions[j, :]).real
+
+            # Add mixing of previous iteration
+            disp2 = mixing * disp2 + (1-mixing) * disp1
+
             res = np.abs(disp1 - disp2).max()
             disp1 = disp2
             if verbose:
                 print('disloc SCF', i, '|d1-d2|_inf =', res)
-            i += 1
-            if i > max_iter:
-                raise RuntimeError('Self-consistency' +
-                                   f'did not converge in {max_iter} cycles')
-        return disp2
-
-    def displacements(self, bulk_positions, center, use_atomman=True,
-                      self_consistent=True, tol=1e-6, max_iter=100, verbose=True):
-        
-        if use_atomman:
-            # Use stroh solver from atomman package
-            if self.stroh is None:
-                self.init_stroh()
-            solver = self.stroh.displacement
+            if res < tol:
+                return disp2
+            
+        raise RuntimeError('Self-consistency ' +
+                            f'did not converge in {max_iter} cycles')
+
+    def displacements(self, bulk_positions, core_positions, method="atomman",
+                      self_consistent=True, tol=1e-6, max_iter=100, verbose=True, mixing=0.5):
+        '''
+        Compute dislocation displacements self-consistently, with max_iter capping the number of iterations
+        Each dislocation core uses a separate solver, which computes the displacements associated with positions 
+        relative to that core (i.e. that the core is treated as lying at [0, 0, 0])
 
-        else:
-            # Use AnistoropicDislocation class from matscipy.dislocation
-            if self.ADstroh is None:
-                self.init_anisotropic_dislocation()
-            solver = self.ADstroh.displacement
+        Parameters
+        ----------
+        bulk_positions: np.array
+            (N, 3) array of positions of atoms in the bulk configuration
+        core_positions: np.array
+            Positions of each dislocation core (shape = (ncores, 3))
+        use_atomman: bool
+            Use the Stroh solver included in atomman (requires atomman package) to
+            solve for displacements, or use the AnisotropicDislocation class
+        self_consistent: bool
+            Whether to detemine the dislocation displacements in a self-consistent manner
+            (self_consistent=False is equivalent to max_iter=0)
+        tol: float
+            Displacement tolerance for the convergence of the self-consistent loop
+        max_iter: int
+            Maximum number of iterations of the self-consistent cycle to run before throwing 
+            an exception
+        verbose: bool
+            Enable/Disable printing progress of the self-consistent cycle each iteration
+
+        Returns
+        -------
+        displacements: np.array
+        shape (N, 3) array of displacements to apply to bulk_positions to form
+        the dislocation structure
+
+        Raises a RuntimeError if max_iter is reached without convergence
+        '''
+        
+        solvers = self.get_solvers(method)
 
         if not self_consistent:
             max_iter = 0
 
-        disp = self.self_consistent_displacements(solver, bulk_positions, center, 
-                                                  tol, max_iter, verbose)
+        disp = self.self_consistent_displacements(solvers, bulk_positions, core_positions, 
+                                                  tol, max_iter, verbose, mixing)
 
         return disp
     
 
     def build_cylinder(self, radius,
                        core_position=np.array([0., 0., 0.]),
-                       extension=np.array([0., 0., 0.]),
-                       fix_width=10.0, self_consistent=None):
-
-        if self_consistent is None:
-            self_consistent = self.self_consistent
-
-        extent = np.array([2 * (radius + fix_width),
-                           2 * (radius + fix_width), 1.])
-        repeat = np.ceil(extent / np.diag(self.unit_cell.cell)).astype(int)
-
-        # if the extension and core position is
-        # within the unit cell, do not add extra unit cells
-        repeat_extension = np.floor(2.0 * extension /
-                                    np.diag(self.unit_cell.cell)).astype(int)
-        repeat_core_position = np.floor(2.0 * core_position /
-                                        np.diag(self.unit_cell.cell)).astype(int)
-
-        extra_repeat = np.stack((repeat_core_position,
-                                 repeat_extension)).max(axis=0)
-
-        repeat += extra_repeat
-
-        repeat[2] = 1  # exactly one cell in the periodic direction
-
-        # ensure correct parity in x and y directions
-        if repeat[0] % 2 != self.parity[0]:
-            repeat[0] += 1
-        if repeat[1] % 2 != self.parity[1]:
-            repeat[1] += 1
-
-        bulk = self.unit_cell * repeat
-        # in order to get center from an atom to the desired position
-        # we have to move the atoms in the opposite direction
-        bulk.positions -= self.unit_cell_core_position
+                       extension=np.array([0, 0, 0]),
+                       fix_width=10.0, self_consistent=None,
+                       method="atomman",
+                       verbose=True):
+        '''
+        Build dislocation cylinder for single dislocation system
 
-        center = np.diag(bulk.cell) / 2
-        shifted_center = center + core_position
-
-        cylinder_mask = get_centering_mask(bulk, radius,
-                                           core_position, extension)
-
-        # add square borders for the case of large extension or core position
-        x, y, _ = bulk.positions.T
-        x_mask = x - center[0] < extension[0] + core_position[0]
-        x_mask = x_mask * (x - center[0] > 0)
-        y_mask = np.abs(y - center[1]) < radius
-        square_mask = y_mask & x_mask
-
-        final_mask = square_mask | cylinder_mask
-
-        bulk = bulk[final_mask]
-
-        # disloc is a copy of bulk with displacements applied
-        disloc = bulk.copy()
-
-        disloc.positions += self.displacements(bulk.positions, shifted_center,
-                                               self_consistent=self_consistent)
-
-        r = np.sqrt(((bulk.positions[:, [0, 1]]
-                      - center[[0, 1]])**2).sum(axis=1))
-
-        fix_mask = r > radius - fix_width
-
-        shifted_r = np.sqrt(((bulk.positions[:, [0, 1]] -
-                              shifted_center[[0, 1]]) ** 2).sum(axis=1))
-
-        shifted_fix_max = shifted_r > radius - fix_width
-        extension = np.array(extension)
-        extended_center = center + extension
-        extended_r = np.sqrt(((bulk.positions[:, [0, 1]] -
-                               extended_center[[0, 1]]) ** 2).sum(axis=1))
-        extended_fix_max = extended_r > radius - fix_width
-        final_fix_mask = fix_mask & shifted_fix_max & extended_fix_max
-
-        x, y, _ = bulk.positions.T
-        x_mask = x - center[0] < extension[0] + core_position[0]
-        x_mask = x_mask * (x - center[0] > 0)
-        y_mask = np.abs(y - center[1]) > radius - fix_width
-
-        # change mask only between the centers of the cylinders
-        final_fix_mask[x_mask] = y_mask[x_mask]
-
-        disloc.set_array('fix_mask', final_fix_mask)
-        disloc.set_constraint(FixAtoms(mask=final_fix_mask))
+        Parameters
+        ----------
+        radius: float
+            Radius of bulk surrounding the dislocation core
+        core_position: np.array
+            Vector offset of the core position from default site
+        extension: np.array
+            Add extra bulk to the system, to also surround a fictitious core
+            that is at position core_position + extension
+        fix_width: float
+            Defines a region to apply the FixAtoms ase constraint to
+            Fixed region is given by (radius - fix_width) <= r <= radius,
+            where the position r is measured from the dislocation core
+            (or from the glide path, if extra_glide_widths is given)
+        self_consistent: bool
+            Controls whether the displacement field used to construct the dislocation is converged 
+            self-consistently. If None (default), the value of self.self_consistent is used
+        use_atomman: bool
+            Use the Stroh solver included in atomman (requires atomman package) to
+            solve for displacements, or use the AnisotropicDislocation class
+        '''
+
+        core_positions = np.array([
+            core_position + self.unit_cell_core_position
+        ])
+
+        bulk, disloc, core_positions = self._build_bulk_cyl(radius, core_positions, fix_width,
+                                                            extension, self_consistent, method, verbose)
+
+        disloc.info["core_positions"] = [list(core_positions[0, :])]
+        disloc.info["burgers_vectors"] = [list(self.burgers)]
+        disloc.info["dislocation_types"] = [self.name]
+        disloc.info["dislocation_classes"] = [str(self.__class__)]
 
         # adding vacuum and centering breaks consistency
         # of displacement =  dislo.positions - bulk.positions
@@ -2646,6 +2902,7 @@ def build_glide_configurations(self, radius,
 
         _, disloc_fin = self.build_cylinder(radius,
                                             core_position=final_core_position,
+                                            extension=-final_core_position,
                                             **kwargs)
         if average_positions:
             # get the fixed atoms constrain
@@ -2772,6 +3029,165 @@ def build_impurity_cylinder(self, disloc, impurity, radius,
         impurities_disloc.cell = disloc.cell
 
         return impurities_disloc
+    
+    @staticmethod
+    def view_cyl(system, scale=0.5, CNA_color=True, add_bonds=False,
+                     line_color=[0, 1, 0], disloc_names=None, hide_arrows=False):
+        '''
+        NGLview-based visualisation tool for structures generated from the dislocation class
+
+        Parameters
+        ----------
+        system: ase.Atoms
+            Dislocation system to view
+        scale: float
+            Adjust radiusScale of add_ball_and_stick (add_bonds=True) or add_spacefill (add_bonds=False)
+        CNA_color: bool
+            Turn on atomic/bond colours based on Common Neighbour Analysis structure identification
+        add_bonds: bool
+            Add atomic bonds
+        line_color: list or list of lists
+            [R, G, B] values for the colour of dislocation lines
+            if a list of lists, line_color[i] defines the colour of the ith dislocation (see structure.info["dislocation_types"])
+        disloc_names: list of str
+            Manual override for the automatic naming of the dislocation lines 
+            (see structure.info["dislocation_types"] for defaults)
+        hide_arrows: bool
+            Hide arrows for placed on the dislocation lines
+        '''
+
+        def _plot_disloc_line(view, disloc_pos, disloc_name, z_length, color=[0, 1, 0]):
+            '''Add dislocation line to the view as a cylinder and two cones.
+            The cylinder is hollow by default so the second cylinder is needed to close it.
+            In case partial distance is provided, two dislocation lines are added and are shifter accordingly.
+            
+            Parameters
+            ----------
+            view: NGLview
+                NGLview plot of the atomic structure
+            disloc_pos: list or array
+                Position of the dislocation core
+            color: list
+                [R, G, B] colour of the dislocation line
+            '''
+
+            view.shape.add_cylinder((disloc_pos[0], disloc_pos[1], -2.0), 
+                                    (disloc_pos[0], disloc_pos[1], z_length - 0.5),
+                                    color,
+                                    [0.3],
+                                    disloc_name)
+            
+            view.shape.add_cone((disloc_pos[0], disloc_pos[1], -2.0), 
+                                (disloc_pos[0], disloc_pos[1], 0.5),
+                                color,
+                                [0.3],
+                                disloc_name)
+            
+            view.shape.add_cone((disloc_pos[0], disloc_pos[1], z_length - 0.5), 
+                                (disloc_pos[0], disloc_pos[1], z_length + 1.0),
+                                color,
+                                [0.55],
+                                disloc_name)
+    
+        from nglview import show_ase, ASEStructure
+        from matscipy.utils import get_structure_types
+        # custom tooltip for nglview to avoid showing molecule and residue names
+        # that are not relevant for dislocation 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";
+                            }
+                        });
+                        """
+        
+        
+        # Check system contains all keys in structure.info we need to plot the dislocations
+        for expected_key in ["core_positions", "burgers_vectors", "dislocation_types", "dislocation_classes"]:
+            
+            if not expected_key in system.info.keys():
+                raise RuntimeError(f"{expected_key} not found in system.info, regenerate system from self.build_cylinder()")
+        
+        # Validate line_color arg
+        if type(line_color) in [list, np.array, tuple]:
+            if type(line_color[0]) in [list, np.array, tuple]:
+                # Given an RGB value per dislocation
+                colours = line_color
+            else:
+                # Assume we're given a single RBG value for all dislocs
+                colours = [line_color] * len(system.info["dislocation_types"])
+                
+        # Check if system contains a diamond dislocation (e.g. DiamondGlideScrew, Diamond90degreePartial)
+        diamond_structure = "Diamond" in system.info["dislocation_classes"][0]
+
+        atom_labels, structure_names, colors = get_structure_types(system, 
+                                                                diamond_structure=diamond_structure)
+        if disloc_names is None:
+            disloc_names = system.info["dislocation_types"]
+
+        disloc_positions = system.info["core_positions"]
+
+        if len(disloc_names) != len(disloc_positions):
+            raise RuntimeError("Number of dislocation positions is not equal to number of dislocation types/names")
+
+        view = show_ase(system)
+        view.hide([0])
+        
+        if add_bonds: # add bonds between all atoms to have bonds between structures
+            component = view.add_component(ASEStructure(system), default_representation=False, name='between structures')
+            component.add_ball_and_stick(cylinderOnly=True, radiusType='covalent', radiusScale=scale, aspectRatio=0.1)
+        
+        # Add structure and struct colours
+        for structure_type in np.unique(atom_labels):
+            # every structure type is a different component
+            mask = atom_labels == structure_type
+            component = view.add_component(ASEStructure(system[mask]), 
+                                        default_representation=False, name=str(structure_names[structure_type]))
+            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:
+                if add_bonds:
+                    component.add_ball_and_stick(radiusType='covalent', radiusScale=scale)
+                else:
+                    component.add_spacefill(radiusType='covalent', radiusScale=scale)
+                        
+        component.add_unitcell()
+
+        # Plot dislocation lines
+        z_length = system.cell[2, 2]
+
+        if not hide_arrows:
+            for i in range(len(disloc_positions)):
+                _plot_disloc_line(view, disloc_positions[i], disloc_names[i] + " dislocation line", z_length, colours[i])
+
+        view.camera = 'orthographic'
+        view.parameters = {"clipDist": 0}
+
+        #view._remote_call("setSize", target="Widget", args=["400px", "300px"])
+        view.layout.width = '100%'
+        view.layout.height = '300px'
+        view.center()
+
+        view._js(tooltip_js)
+        return view
 
 
 class CubicCrystalDissociatedDislocation(CubicCrystalDislocation, metaclass=ABCMeta):
@@ -2784,7 +3200,7 @@ class CubicCrystalDissociatedDislocation(CubicCrystalDislocation, metaclass=ABCM
     # 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"):
+    def __init__(self, a, C11=None, C12=None, C44=None, C=None, symbol="W"):
         """
         This class represents a dissociated dislocation in a cubic crystal
         with burgers vector b = b_left + b_right.
@@ -2800,8 +3216,12 @@ def __init__(self, a, C11, C12, C44, symbol="W"):
                         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)
+        if not (isinstance(self.left_dislocation, CubicCrystalDislocation) or \
+                isinstance(self.left_dislocation, CubicCrystalDissociatedDislocation)):
+            self.left_dislocation = self.left_dislocation(a, C11, C12, C44, C, symbol)
+        if not (isinstance(self.right_dislocation, CubicCrystalDislocation) or \
+                isinstance(self.right_dislocation, CubicCrystalDissociatedDislocation)):
+            self.right_dislocation = self.right_dislocation(a, C11, C12, C44, C, symbol)
 
         # Change disloc burgers vectors, if requested
         if self.new_left_burgers is not None:
@@ -2812,15 +3232,8 @@ def __init__(self, a, C11, C12, C44, symbol="W"):
         # 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)
+        super().__init__(a, C11, C12, C44, C, symbol)
 
         # Validation of disloc inputs
         try:
@@ -2829,16 +3242,14 @@ def __init__(self, a, C11, C12, C44, symbol="W"):
                                            self.burgers)
         except AssertionError as error:
             print(error)
-            raise ValueError("Burgers vectors of left and right disloctions" +
+            raise ValueError("Burgers vectors of left and right dislocations " +
                              "do not add up to the desired vector")
 
         # checking that parameters of
         # left and right dislocations are the same
         try:
             assert left_dislocation.alat == right_dislocation.alat
-            assert left_dislocation.C11 == right_dislocation.C11
-            assert left_dislocation.C12 == right_dislocation.C12
-            assert left_dislocation.C44 == right_dislocation.C44
+            assert np.allclose(left_dislocation.C, right_dislocation.C)
 
             np.testing.assert_equal(left_dislocation.unit_cell.get_chemical_symbols(),
                                     right_dislocation.unit_cell.get_chemical_symbols())
@@ -2868,11 +3279,97 @@ def __init__(self, a, C11, C12, C44, symbol="W"):
             raise ValueError("Parameters of left and right" +
                              "partials must be the same")
 
+    # classmethod properties that get props from left dislocation
+    # Used so e.g. cls.crystalstructure is setup prior to __init__
+    # as is the case with CubicCrystalDislocation
+
+    @classmethod
+    @property
+    def crystalstructure(cls):
+        return cls.left_dislocation.crystalstructure
+    
+    @classmethod
+    @property
+    def axes(cls):
+        return cls.left_dislocation.axes
+    
+    @classmethod
+    @property
+    def unit_cell_core_position_dimensionless(cls):
+        return cls.left_dislocation.unit_cell_core_position_dimensionless
+    
+    @classmethod
+    @property
+    def parity(cls):
+        return cls.left_dislocation.parity
+    
+    
+    @classmethod
+    @property
+    def n_planes(cls):
+        return cls.left_dislocation.n_planes
+    
+    @classmethod
+    @property
+    def self_consistent(cls):
+        return cls.left_dislocation.self_consistent
+    
+    @classmethod
+    @property
+    def glide_distance_dimensionless(cls):
+        return cls.left_dislocation.glide_distance_dimensionless
+    
+    @property
+    def alat(self):
+        return self.left_dislocation.alat
+    
+    @property
+    def unit_cell(self):
+        return self.left_dislocation.unit_cell
+
+    def invert_burgers(self):
+        '''
+        Modify dislocation to produce same dislocation with opposite burgers vector
+        '''
+        self.burgers_dimensionless *= -1
+
+        self.left_dislocation.invert_burgers()
+        self.right_dislocation.invert_burgers()
+
+        # Flip dislocations
+        self.left_dislocation, self.right_dislocation = self.right_dislocation, self.left_dislocation
+
+    def get_solvers(self, method="atomman"):
+        '''
+        Overload of CubicCrystalDislocation.get_solvers
+        Get list of dislocation displacement solvers
+
+        Parameters
+        ----------
+        use_atomman: bool
+            Use the Stroh solver included in atomman (requires atomman package) to
+            solve for displacements, or use the AnisotropicDislocation class
+
+        Returns
+        -------
+        solvers: list
+            List of functions f(positions) = displacements
+        '''
+        solvers = self.left_dislocation.get_solvers(method)
+
+        solvers.extend(
+            self.right_dislocation.get_solvers(method)
+        )
+
+        return solvers
 
     def build_cylinder(self, radius, partial_distance=0,
                        core_position=np.array([0., 0., 0.]),
-                       extension=np.array([0., 0., 0.]),
-                       fix_width=10.0, self_consistent=None):
+                       extension=np.array([[0., 0., 0.],
+                                          [0., 0., 0.]]),
+                       fix_width=10.0, self_consistent=None,
+                       method="atomman",
+                       verbose=True):
         """
         Overloaded function to make dissociated dislocations.
         Partial distance is provided as an integer to define number
@@ -2886,49 +3383,515 @@ def build_cylinder(self, radius, partial_distance=0,
             distance between partials (SF length) in number of glide distances.
             Default is 0 -> non dissociated dislocation
             with b = b_left + b_right is produced
+        use_atomman: bool
+            Use the Stroh solver included in atomman (requires atomman package) to
+            solve for displacements, or use the AnisotropicDislocation class
+        extension: np.array
+            Shape (2, 3) array giving additional extension vectors from each dislocation core.
+            Used to add extra bulk, e.g. to set up glide configurations.
         """
-
-        if self_consistent is None:
-            self_consistent = self.self_consistent
-
+        
         partial_distance_Angstrom = np.array(
             [self.glide_distance * partial_distance, 0.0, 0.0])
+            
+
+        core_positions = np.array([
+            core_position + self.unit_cell_core_position,
+            core_position + self.unit_cell_core_position + partial_distance_Angstrom
+        ])
+
+        bulk, disloc, core_positions = self._build_bulk_cyl(radius, core_positions, fix_width, extension,
+                                                            self_consistent, method, verbose)
+
+        if partial_distance > 0:
+            # Specify left & right dislocation separately
+            disloc.info["core_positions"] = [list(core_positions[0, :]), list(core_positions[1, :])]
+            disloc.info["burgers_vectors"] = [list(self.left_dislocation.burgers), list(self.right_dislocation.burgers)]
+            disloc.info["dislocation_types"] = [self.left_dislocation.name, self.right_dislocation.name]
+            disloc.info["dislocation_classes"] = [str(self.left_dislocation.__class__), str(self.right_dislocation.__class__)]
+        else:
+            # Perfect, non-dissociated dislocation, only show values for single dislocation
+            disloc.info["core_positions"] = [list(core_positions[0, :])]
+            disloc.info["burgers_vectors"] = [list(self.burgers)]
+            disloc.info["dislocation_types"] = [self.name]
+            disloc.info["dislocation_classes"] = [str(self.__class__)]
 
-        bulk, disloc = self.left_dislocation.build_cylinder(radius,
-                                                            extension=extension + partial_distance_Angstrom,
-                                                            core_position=core_position,
-                                                            fix_width=fix_width,
-                                                            self_consistent=self_consistent)
+        return bulk, disloc
 
-        _, disloc_right = self.right_dislocation.build_cylinder(radius,
-                                                                core_position=core_position + partial_distance_Angstrom,
-                                                                extension=extension,
-                                                                fix_width=fix_width,
-                                                                self_consistent=self_consistent)
+class CubicCrystalDislocationQuadrupole(CubicCrystalDissociatedDislocation):
+    burgers_dimensionless = np.zeros(3)
+    def __init__(self, disloc_class, *args, **kwargs):
+        '''
+        Initialise dislocation quadrupole class
 
-        u_right = disloc_right.positions - bulk.positions
-        disloc.positions += u_right
+        Arguments
+        ---------
+        disloc_class: Subclass of CubicCrystalDislocation
+            Dislocation class to create (e.g. DiamondGlide90DegreePartial)
 
-        return bulk, disloc
+        *args, **kwargs
+            Parameters fed to CubicCrystalDislocation.__init__()
+        '''
 
-    def displacements(self, bulk_positions, center,
-                      partial_distance=0, **kwargs):
-        """Overloaded function to provide correct displacements
-        for the dissociated dislocation.
-        Partial distance is provided as an integer to define number
-        of glide distances between two partials.
-        """
+        if isinstance(disloc_class, CubicCrystalDislocation):
+            disloc_cls = disloc_class.__class__
+        else:
+            disloc_cls = disloc_class
 
-        partial_distance_Angstrom = np.array(
-            [self.glide_distance * partial_distance, 0.0, 0.0])
+        self.left_dislocation = disloc_cls(*args, **kwargs)
+        self.right_dislocation = disloc_cls(*args, **kwargs)
+
+        self.left_dislocation.invert_burgers()
+        
+        super().__init__(*args, **kwargs)
+
+        self.glides_per_unit_cell = np.floor((self.unit_cell.cell[0, 0] + self.unit_cell.cell[1, 0]) 
+                                              / (self.glide_distance - 1e-2)).astype(int)
+
+    # Overload all of the classmethod properties from CubicCrystalDissociatedDislocation, as these cannot be
+    # known at class level (don't know which disloc to make a quadrupole of, so can't know what axes should be)
+        
+    @property
+    def crystalstructure(self):
+        return self.left_dislocation.crystalstructure
+    
+    @property
+    def axes(self):
+        return self.left_dislocation.axes
+    
+    @property
+    def unit_cell_core_position_dimensionless(self):
+        return self.left_dislocation.unit_cell_core_position_dimensionless
+    
+    @property
+    def parity(self):
+        return self.left_dislocation.parity
+
+    @property
+    def n_planes(self):
+        return self.left_dislocation.n_planes
+    
+    @property
+    def self_consistent(self):
+        return self.left_dislocation.self_consistent
+    
+    @property
+    def glide_distance_dimensionless(self):
+        return self.left_dislocation.glide_distance_dimensionless
+
+
+    def periodic_displacements(self, positions, v1, v2, core_positions, disp_tol=1e-3, max_neighs=60, 
+                               verbose="periodic", **kwargs):
+        '''
+        Calculate the stroh displacements for the periodic structure defined by 2D lattice vectors v1 & v2
+        given the core positions.
+
+        positions: np.array
+            Atomic positions to calculate displacements at
+        v1, v2: np.array
+            Lattice vectors defining periodicity in XY plane (dislocation line is in Z)
+        core_positions: np.array
+            positions of the all dislocation cores (including all partial cores)
+        disp_tol: float
+            Displacement tolerance controlling the termination of displacement convergence
+        max_neighs: int
+            Maximum number of nth-nearest neighbour cells to consider
+            max_neighs=0 is just the current cell, max_neighs=1 also includes nearest-neighbour cells, ...
+        verbose: bool, str or None
+            Controls verbosity of the displacement solving
+            False, None, or "off" : No printing
+            "periodic" (default) : Prints status of the periodic displacement convergence
+            True : Also prints status of self-consistent solutions for each dislocation
+                    (i.e. verbose=True for CubicCrystalDislocation.displacements())
+        **kwargs
+            Other keyword arguments fed to self.displacements() (e.g. method="adsl")
+        returns
+        -------
+        disps: np.array
+            displacements
+        '''
+
+        ncores = core_positions.shape[0]
+
+        if type(verbose) == str:
+            if verbose.lower() == "off":
+                periodic_verbose = False
+                disloc_verbose = False
+            elif verbose.lower() == "periodic":
+                periodic_verbose = True
+                disloc_verbose = False
+            else:
+                warnings.warn(f"Unrecognised string verbosity '{verbose}'")
+                periodic_verbose = True
+                disloc_verbose = False
+        elif type(verbose) == bool:
+            if verbose:
+                periodic_verbose = True
+                disloc_verbose = True
+            else:
+                periodic_verbose = False
+                disloc_verbose = False
+        elif verbose is None:
+            periodic_verbose = False
+            disloc_verbose = False
+        else:
+            warnings.warn(f"verbose argument '{verbose}' of type {type(verbose)} not understood.")
+            periodic_verbose = True
+            disloc_verbose = False
+
+        if "self_consistent" in kwargs:
+            # Disable disloc verbosity if no self-consistent solve of disloc displacements
+            # CubicCrystalDislocation.displacements doesn't print unless SCF is turned on 
+            disloc_verbose = disloc_verbose * bool(kwargs["self_consistent"])
+            self_consistent = bool(kwargs[self_consistent])
+        else:
+            kwargs["self_consistent"] = self.self_consistent
+
+        displacements = np.zeros_like(positions)
+        disp_update = np.zeros_like(positions)
+        
+        for neigh_idx in range(max_neighs):
+            disp_update = np.zeros_like(positions)
+            for ix in range(-neigh_idx, neigh_idx+1):
+                for iy in range(-neigh_idx, neigh_idx+1):
+                    if ix < neigh_idx and ix > -neigh_idx and iy < neigh_idx and iy > -neigh_idx:
+                        # Looping over cells covered by previous iteration, skip
+                        continue
+
+                    # In a new cell
+                    if disloc_verbose:
+                        print(f"Displacements for cell {ix} {iy}")
+
+                    core_pos = core_positions.copy()
+                    for i in range(ncores):
+                        core_pos[i, :] += ix * v1 + iy * v2
+                    disp_update += self.displacements(
+                        positions, core_pos,
+                        verbose=disloc_verbose, **kwargs)
+
+            displacements += disp_update
+            max_disp_change = np.max(np.linalg.norm(disp_update, axis=-1))
+
+            if periodic_verbose:
+                print(f"Periodic displacement iteration {neigh_idx} -> max(|dr|) = {np.round(max_disp_change, 4)}")
+
+            if max_disp_change < disp_tol and ix > 0:
+                
+                if periodic_verbose:
+                    print(f"Periodic displacements converged to r_tol={disp_tol} after {neigh_idx} iterations")
+
+                return displacements
+        warnings.warn("Periodic quadrupole displacments did not converge!")
+        return displacements
+
+
+    def build_quadrupole(self, glide_separation=4, partial_distance=0, extension=0, verbose='periodic', 
+                         left_offset=None, right_offset=None, **kwargs):
+        '''
+        Build periodic quadrupole cell
+
+        Parameters
+        ----------
+        glide_separation: int
+            Number of glide distances separating the two dislocation cores
+        partial_distance: int or list of int
+            Distance between dissociated partial dislocations
+            If a list, specify a partial distance for each dislocation
+            If an int, specify a partial distance for both dislocations
+        extension: int
+            Argument to modify the minimum length of the cell in the glide (x) direction.
+            Cell will always be at least 
+            (2 * glide_separation) + partial_distance + extension glide vectors long
+            Useful for setting up images for glide barrier calculations, as 
+            start and end structures must be the same size for the NEB method
+        verbose: bool, str, or None
+            Verbosity value to be fed to CubicCrystalDislocationQuadrupole.periodic_displacements
+        left_offset, right_offset: np.array
+            Translational offset (in Angstrom) for the left and right dislocation cores 
+            (i.e. for the +b and -b cores)
+        
+        '''
+
+        if left_offset is None:
+            left_offset = np.zeros(3)
+        if right_offset is None:
+            right_offset = np.zeros(3)
+
+        if glide_separation < 1:
+            raise RuntimeError("glide_separation should be >= 1")
+        elif glide_separation < 3:
+            msg = "glide_separation is very small. Resulting structure may be very unstable."
+            warnings.warn(msg)
+
+        assert extension >= 0
+        assert partial_distance >= 0
+
+        core_separation = glide_separation * self.glide_distance
+        core_vec = np.array([core_separation, 0, 0])
+
+        partial_vec = np.array([partial_distance * self.glide_distance, 0, 0])
+        # Additional space to add to the cell
+        extra_extension = np.array([extension + partial_distance, 0, 0]) * self.glide_distance
+        
+        # Replicate the unit cell enough times to fill the target cell
+        cell = self.unit_cell.cell[:, :].copy()
+
+        # Integer number of glide distances which fit in a single unit cell
+        glides_per_unit_cell_x = self.glides_per_unit_cell
+        glides_per_unit_cell_y = np.floor(np.linalg.norm(cell[1, :]) / (self.glide_distance - 1e-2)).astype(int)
+
+        # Number of unit cells in x (glide) direction
+        xreps = np.ceil((2*glide_separation + extension + 2*partial_distance) / glides_per_unit_cell_x).astype(int)
+
+
+        # Number of unit cells in y direction
+        yreps = np.ceil(glide_separation / glides_per_unit_cell_y).astype(int)
+
+        quad_bulk = self.unit_cell.copy() * (xreps, yreps, 1)
+        cell = quad_bulk.cell[:, :].copy()
+
+
+        # New cell based on old cell vectors
+        # Rhomboid shape enclosing both cores, + any stacking fault
+        new_cell = np.array([
+            cell[1, :] + 0.5 * cell[0, :],
+            -cell[1, :] + 0.5 *  cell[0, :],
+            cell[2, :]
+        ])
+
+        if self.crystalstructure == "bcc" and np.floor(glide_separation).astype(int) % 3 == 1:
+            new_cell[0, 0] -= self.glide_distance * 3/2
+            new_cell[1, 0] += self.glide_distance * 3/2
+            
+        quad_bulk.set_cell(new_cell)
+        quad_bulk.wrap()
+        pos = quad_bulk.get_positions()
+
+        quad_disloc = quad_bulk.copy()
+
+        # Get the core positions, translate everything so the cores are central in the cell
+        lens = np.sum(new_cell, axis=0)
+        core_pos_1 = lens/2 - 0.5 * core_vec - 0.5 * partial_vec
+        core_pos_2 = core_pos_1 + core_vec
+        pos += core_pos_1 - self.left_dislocation.unit_cell_core_position
+
+        core_pos_1 += left_offset
+        core_pos_2 += right_offset
+
+        core_positions = np.array([
+            core_pos_1,
+            core_pos_2
+        ])
+
+        if partial_distance > 0.0:
+            partial_core_pos = core_positions.copy()
+            partial_core_pos[:, 0] += partial_vec[0]
+            old_core_positions = core_positions
+            
+            core_positions = np.array([
+                old_core_positions[0, :],
+                partial_core_pos[0, :],
+                old_core_positions[1, :],
+                partial_core_pos[1, :]
+            ])
+
+        quad_bulk.set_positions(pos)
+        quad_bulk.wrap()
+        pos = quad_bulk.get_positions()
+
+        # Apply disloc displacements to quad_disloc
+        disps = self.periodic_displacements(pos, new_cell[0, :], new_cell[1, :], core_positions, 
+                                            verbose=verbose, **kwargs)
+
+        pos += disps
+        quad_disloc.set_positions(pos)
+        quad_disloc.wrap()
+
+        if partial_distance > 0 and isinstance(self.right_dislocation, CubicCrystalDissociatedDislocation):
+            # Dissociated dislocation quadrupole
+            # 4 dislocations in total
+            # "Left" dissociated dislocation
+            pos1 = core_pos_1
+            pos2 = core_pos_1 + np.array([partial_distance * self.glide_distance, 0, 0])
+            # "Right" dissociated dislocation
+            pos3 = core_pos_2
+            pos4 = core_pos_2 + np.array([partial_distance * self.glide_distance, 0, 0])
+                                         
+            quad_disloc.info["core_positions"] = [list(pos1), list(pos2), list(pos3), list(pos4)]
+
+            # Quadrupole is a "dissociated dissociated" dislocation
+            # 4 partial dislocs are .L.L, .L.R, .R.L, .R.R
+            # (L = left_dislocation, R = right_dislocation)
+
+            for key, prop in [
+                ("burgers_vectors", "burgers"),
+                ("dislocation_types", "name"),
+                ("dislocation_classes", "__class__")
+            ]:
+                quad_disloc.info[key] = [
+                    getattr(self.left_dislocation.left_dislocation, prop),
+                    getattr(self.left_dislocation.right_dislocation, prop),
+                    getattr(self.right_dislocation.left_dislocation, prop),
+                    getattr(self.right_dislocation.right_dislocation, prop)
+                ]
+
+            # Type conversions
+            quad_disloc.info["burgers_vectors"] = [list(burgers) for burgers in quad_disloc.info["burgers_vectors"]]
+            quad_disloc.info["dislocation_classes"] = [str(name) for name in quad_disloc.info["dislocation_classes"]]
+        else:
+            # Perfect, non-dissociated dislocation, only show values for single dislocation
+            
+            quad_disloc.info["core_positions"] = [list(core_pos_1), list(core_pos_2)]
+            quad_disloc.info["burgers_vectors"] = [list(self.left_dislocation.burgers), list(self.right_dislocation.burgers)]
+            quad_disloc.info["dislocation_types"] = [self.left_dislocation.name, self.right_dislocation.name]
+            quad_disloc.info["dislocation_classes"] = [str(self.left_dislocation.__class__), str(self.right_dislocation.__class__)]
+
+
+
+        return quad_bulk, quad_disloc
+    
+    def build_glide_quadrupoles(self, nims, invert_direction=False, glide_left=True, glide_right=True, 
+                                left_offset=None, right_offset=None, *args, **kwargs):
+        '''
+        Construct a sequence of quadrupole structures providing an initial guess of the dislocation glide
+        trajectory
+
+        Parameters
+        ----------
+        nims: int
+            Number of images to generate
+        invert_direction: bool
+            Invert the direction of the glide
+            invert_direction=False (default) glides in the +x direction
+            invert_direction=True glides in the -x direction
+        glide_left, glide_right: bool
+            Flags for toggling whether the left/right cores
+            are allowed to glide
+        *args, **kwargs
+            Fed to self.build_quadrupole()
+        
+        '''
+        
+        glide_offsets = np.linspace(0, self.glide_distance, nims, endpoint=True)
+
+        if invert_direction:
+            glide_offsets *= -1
+
+        if left_offset is not None:
+            _loff = left_offset
+        else:
+            _loff = np.zeros(3)
+
+        if right_offset is not None:
+            _roff = right_offset
+        else:
+            _roff = np.zeros(3)
+
+        images = []
+
+        for i in range(nims):
+            left_offset = np.array([glide_offsets[i], 0, 0]) if glide_left else np.zeros(3)
+            right_offset = np.array([glide_offsets[i], 0, 0]) if glide_right else np.zeros(3)
+
+            left_offset += _loff
+            right_offset += _roff
+
+            images.append(self.build_quadrupole(
+                left_offset=left_offset,
+                right_offset=right_offset,
+                *args,
+                **kwargs
+            )[1])
+        return images
+    
+    def build_kink_quadrupole(self, z_reps, layer_tol=1e-1, invert_direction=False, *args, **kwargs):
+        '''
+        Construct a quadrupole structure providing an initial guess of the dislocation kink
+        mechanism. Produces a periodic array of kinks, where each kink is 
+        z_reps * self.unit_cell[2, 2] Ang long
+
+        Parameters
+        ----------
+        z_reps: int
+            Number of replications of self.unit_cell to use per kink
+            Should be >1
+        layer_tol: float
+            Tolerance for trying to determine the top atomic layer of the cell
+            (required in order to complete periodicity)
+            Top layer is defined as any atom with z component >= max(atom_z) - layer_tol
+        invert_direction: bool
+            Invert the direction of the glide
+            invert_direction=False (default) kink in the +x direction
+            invert_direction=True kink in the -x direction
+        *args, **kwargs
+            Fed to self.build_quadrupole() & self.build_glide_quadrupoles
+        
+        '''
+        assert z_reps > 1
+
+        direction = -1 if invert_direction else 1
 
-        left_u = self.left_dislocation.displacements(bulk_positions, center,
-                                                     **kwargs)
-        right_u = self.right_dislocation.displacements(bulk_positions,
-                                                       center + partial_distance_Angstrom,
-                                                       **kwargs)
+        # Need both cores to move for the infinite kink structure
+        reps = self.build_glide_quadrupoles(z_reps, glide_left=True, glide_right=True, 
+                                            invert_direction=invert_direction, *args, **kwargs)
+
+        kink_struct = reps[0]
+
+        for image in reps[1:]:
+            kink_struct = stack(kink_struct, image)
+        
+        cell = kink_struct.cell[:, :]
+
+        cell[2, 0] += self.glide_distance * direction
+
+        kink_struct.set_cell(cell)
+        kink_struct.wrap()
+
+        glide_parity = (-direction) % self.glides_per_unit_cell
+
+        if glide_parity:
+            # Cell won't be periodic if multiple glides needed to complete periodicity
+            # glide_parity gives number of layers required to remove
+
+            for i in range(glide_parity):
+                atom_heights = kink_struct.get_positions()[:, 2]
+                top_atom = np.max(atom_heights)
+
+                layer_mask = atom_heights >= top_atom - layer_tol
+
+                avg_layer_pos = np.average(atom_heights[layer_mask])
+
+                kink_struct = kink_struct[~layer_mask]
+                cell[2, 2] = avg_layer_pos
+                kink_struct.set_cell(cell)
+                kink_struct.wrap()
+        return kink_struct
+    
+    def view_quad(self, system, *args, **kwargs):
+        '''
+        Specialised wrapper of view_cyl for showing quadrupoles, as the quadrupole cell 
+        causes the plot to be rotated erroneously
+        
+        Takes the same args and kwargs as view_cyl
+        '''
+
+        view = self.view_cyl(system, hide_arrows=True, *args, **kwargs)
+
+        cell = system.cell[:, :]
+        
+        # Default rotation has cell[0, :] pointing in x
+        # Rotate such that cell[0, :] has angle theta to x
+        # cell[1, :] has -theta to x
+        dot = (cell[0, :] @ cell[1, :]) / (np.linalg.norm(cell[0, :]) * np.linalg.norm(cell[1, :]))
+        rot_angle = -np.arccos(dot)/2 + np.pi
 
-        return left_u + right_u
+        # Rotate system by rot_angle about [0, 0, 1]
+        view.control.spin([0, 0, 1], rot_angle)
+        return view
+   
+# TODO: If non-cubic dislocation classes are implemented, need to build an
+# interface to make "Quadrupole" work for both
+Quadrupole = CubicCrystalDislocationQuadrupole
 
 
 class BCCScrew111Dislocation(CubicCrystalDislocation):
@@ -2939,6 +3902,7 @@ class BCCScrew111Dislocation(CubicCrystalDislocation):
     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
+    name = "1/2<111> screw"
 
 
 class BCCEdge111Dislocation(CubicCrystalDislocation):
@@ -2950,6 +3914,7 @@ class BCCEdge111Dislocation(CubicCrystalDislocation):
     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
+    name = "1/2<111> edge"
 
 class BCCEdge111barDislocation(CubicCrystalDislocation):
     crystalstructure = "bcc"
@@ -2960,6 +3925,7 @@ class BCCEdge111barDislocation(CubicCrystalDislocation):
     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 = 1
+    name = "1/2<11-1> edge"
 
 class BCCMixed111Dislocation(CubicCrystalDislocation):
     crystalstructure = "bcc"
@@ -2971,6 +3937,7 @@ class BCCMixed111Dislocation(CubicCrystalDislocation):
     # 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])
+    name = "1/2<111> mixed"
 
 
 class BCCEdge100Dislocation(CubicCrystalDislocation):
@@ -2982,6 +3949,7 @@ class BCCEdge100Dislocation(CubicCrystalDislocation):
     unit_cell_core_position_dimensionless = np.array([1/4, 1/4, 0])
     glide_distance_dimensionless = 1.0
     n_planes = 2
+    name = "<100>{110} edge"
 
 
 class BCCEdge100110Dislocation(CubicCrystalDislocation):
@@ -2993,6 +3961,7 @@ class BCCEdge100110Dislocation(CubicCrystalDislocation):
     unit_cell_core_position_dimensionless = np.array([0.5, np.sqrt(2) / 4.0, 0])
     glide_distance_dimensionless = 0.5
     n_planes = 2
+    name = "<100>{001} edge"
 
 
 class DiamondGlide30degreePartial(CubicCrystalDislocation):
@@ -3012,6 +3981,7 @@ class DiamondGlide30degreePartial(CubicCrystalDislocation):
     # This leads to the cut plane crossing one of the atomic planes and
     # thus breaking the stacking fault.
     self_consistent = False
+    name = "1/6<112> 30 degree partial"
 
 
 class DiamondGlide90degreePartial(CubicCrystalDislocation):
@@ -3031,6 +4001,7 @@ class DiamondGlide90degreePartial(CubicCrystalDislocation):
     # This leads to the cut plane crossing one of the atomic planes and
     # thus breaking the stacking fault.
     self_consistent = False
+    name = "1/6<112> 90 degree partial"
 
 
 class DiamondGlideScrew(CubicCrystalDissociatedDislocation):
@@ -3038,6 +4009,7 @@ class DiamondGlideScrew(CubicCrystalDissociatedDislocation):
     new_left_burgers = np.array([2., -1., -1.]) / 6
     left_dislocation = DiamondGlide30degreePartial
     right_dislocation = DiamondGlide30degreePartial
+    name = "1/2<110> screw"
 
 
 class DiamondGlide60Degree(CubicCrystalDissociatedDislocation):
@@ -3045,6 +4017,7 @@ class DiamondGlide60Degree(CubicCrystalDissociatedDislocation):
     new_left_burgers = np.array([2., -1., -1.]) / 6
     left_dislocation = DiamondGlide30degreePartial
     right_dislocation = DiamondGlide90degreePartial
+    name = "1/2<110> 60 degree screw"
 
 
 class FCCScrewShockleyPartial(CubicCrystalDislocation):
@@ -3056,6 +4029,7 @@ class FCCScrewShockleyPartial(CubicCrystalDislocation):
     glide_distance_dimensionless = np.sqrt(6) / 4.0
     n_planes = 2
     unit_cell_core_position_dimensionless = np.array([5/6, 1/9, 0])
+    name = "1/6<112> screw Shockley partial"
 
 
 class FCCEdgeShockleyPartial(CubicCrystalDislocation):
@@ -3067,6 +4041,7 @@ class FCCEdgeShockleyPartial(CubicCrystalDislocation):
     unit_cell_core_position_dimensionless = np.array([0, 1/6 , 0])
     glide_distance_dimensionless = np.sqrt(2)/4
     n_planes = 6
+    name = "1/2<110> edge Shockley partial"
 
 
 class FCCScrew110Dislocation(CubicCrystalDissociatedDislocation):
@@ -3074,6 +4049,7 @@ class FCCScrew110Dislocation(CubicCrystalDissociatedDislocation):
     new_left_burgers = np.array([2., -1., -1.]) / 6
     left_dislocation = FCCScrewShockleyPartial
     right_dislocation = FCCScrewShockleyPartial
+    name = "1/2<110>{111} screw"
 
 
 class FCCEdge110Dislocation(CubicCrystalDissociatedDislocation):
@@ -3082,6 +4058,7 @@ class FCCEdge110Dislocation(CubicCrystalDissociatedDislocation):
     new_left_burgers = np.array([2., -1., -1.]) / 6
     left_dislocation = FCCEdgeShockleyPartial
     right_dislocation = FCCEdgeShockleyPartial
+    name = "1/2<110> edge"
 
 
 class FixedLineAtoms:
diff --git a/matscipy/elasticity.py b/matscipy/elasticity.py
index b7ca2cdb..68b9cd58 100644
--- a/matscipy/elasticity.py
+++ b/matscipy/elasticity.py
@@ -208,6 +208,60 @@ def cubic_to_Voigt_6x6(C11, C12, C44):
                      [  0,  0,  0,  0,C44,  0],
                      [  0,  0,  0,  0,  0,C44]])
 
+def coalesce_elastic_constants(C11=None, C12=None, C44=None, C=None, convention="Cij"):
+    '''
+    Convert either all of C11, C12, C44, or 6x6 C, or 3x3x3x3 C to the form defined by convention.
+
+    convention="cubic": returns C11, C12, C44
+    convention="Cij": returns 6x6 C matrix
+    convention="Cijkl": returns 3x3x3x3 C tensor
+    
+    '''
+    # Aim to convert to 3x3x3x3 first, as this contains the maximal information
+    Cijkl = None
+
+    if C is not None:
+        # Prefer to use full matrix, as this contains more info
+        C = np.array(C)
+        if C.shape == (6, 6):
+            Cijkl = Voigt_6x6_to_full_3x3x3x3(C)
+
+
+        elif C.shape == (3, 3, 3, 3):
+            Cijkl = C.copy()
+        else:
+            # C not of the correct shape
+            warnings.warn("C was not the correct shape. Expected shape (6, 6) or (3, 3, 3, 3)")
+        
+        if Cijkl is not None:
+            # Check if C11, C12, or C44 have also been passed
+            # Ignore these inputs, but warn user
+            has_other_elast = [var is not None for var in [C11, C12, C44]]
+
+            if any(has_other_elast):
+                for i, var in enumerate(["C11", "C12", "C44"]):
+                    if has_other_elast[i]:
+                        warnings.warn(f"Ignoring argument {var} as C was passed")
+
+    if Cijkl is None:
+        # Didn't get correct info from C, look at C11, C12, C44
+        if C11 is not None and C12 is not None and C44 is not None:
+            Cijkl = Voigt_6x6_to_full_3x3x3x3(cubic_to_Voigt_6x6(C11, C12, C44))
+        else:
+            raise RuntimeError("Elastic Constants not correctly defined. "+
+                                "Pass either C, or all of C11, C12, C44.")
+
+    if convention.lower() == "cijkl":
+        # 3x3x3x3 convention
+        return Cijkl
+    elif convention.lower() == "cij":
+        # 6x6 convention
+        return full_3x3x3x3_to_Voigt_6x6(Cijkl)
+    elif convention.lower() == "cubic":
+        return Voigt_6x6_to_cubic(full_3x3x3x3_to_Voigt_6x6(Cijkl))
+    
+    raise ValueError(f"{convention=} not one of 'cubic', 'Cij', 'Cijkl")
+
 def _invariants(s, syy=None, szz=None, syz=None, sxz=None, sxy=None,
                full_3x3_to_Voigt_6=full_3x3_to_Voigt_6_stress):
     """
diff --git a/matscipy/gamma_surface.py b/matscipy/gamma_surface.py
index 11a26b4d..a5642854 100644
--- a/matscipy/gamma_surface.py
+++ b/matscipy/gamma_surface.py
@@ -128,7 +128,7 @@ def __init__(self, a, surface_direction, glide_direction=None,
                                                 symbol=symbol)
         self.offset *= alat
 
-    def _vec_to_miller(self, vec, latex=True):
+    def _vec_to_miller(self, vec, latex=True, brackets="["):
         '''
         Converts np.array vec to string miller notation
 
@@ -137,7 +137,32 @@ def _vec_to_miller(self, vec, latex=True):
             e.g. (100) or (11-2)
         latex: bool
             Use LaTeX expressions to construct representation
+        brackets: str
+            Bracketing style. 
+            Options are "(" or ")" for "()", "{" or "}" for "{}",
+             or "[" or "]" for "[]"
+            Other options can be supplied by passing a len(2) list of str
         '''
+        # Sort out bracketing
+        bracket_types = [
+            [["(", ")", "()"], ["(", ")"]],
+            [["{", "}", r"{}"], ["{", "}"]],
+            [["[", "]", "[]"], ["[", "]"]]
+        ]
+        found = False
+        for key, val in bracket_types:
+            if brackets in key:
+                lbrac = val[0]
+                rbrac = val[1]
+                found = True
+        if not found:
+            if len(brackets) ==2:
+                lbrac = brackets[0]
+                rbrac = brackets[1]
+            else:
+                lbrac = "("
+                rbrac = ")"
+
         int_vec = vec.astype(int)
         l = []
         for item in int_vec:
@@ -145,7 +170,7 @@ def _vec_to_miller(self, vec, latex=True):
                 l.extend(r"$\overline{" + str(-item) + "}$")
             else:
                 l.extend(str(item))
-        return "(" + "".join(l) + ")"
+        return lbrac + "".join(l) + rbrac
 
     def _gen_cellmove_images(self, base_struct, 
                              nx, ny, x_points, y_points, vacuum=0.0):
@@ -451,7 +476,7 @@ def plot_energy_densities(self, Es=None, ax=None,
             units = "J/m$^2$"
         else:
             mul = 1
-            units = "eV/${\AA}^2$"
+            units = "eV/$ \mathrm{A}^2$"
 
         if ax is None:
 
@@ -471,10 +496,10 @@ def plot_energy_densities(self, Es=None, ax=None,
         my_ax.set_xlim([0, np.linalg.norm(self.x_disp)])
         my_ax.set_ylim([0, np.linalg.norm(self.y_disp)])
 
-        my_ax.set_xlabel("Glide in " + self._vec_to_miller(self.surf_directions["x"]) + " ($\AA$)")
-        my_ax.set_ylabel("Glide in " + self._vec_to_miller(self.surf_directions["y"]) + " ($\AA$)")
+        my_ax.set_xlabel("Glide in " + self._vec_to_miller(self.surf_directions["x"]) + " ($ \mathrm{A}$)")
+        my_ax.set_ylabel("Glide in " + self._vec_to_miller(self.surf_directions["y"]) + " ($ \mathrm{A}$)")
 
-        my_ax.set_title(self._vec_to_miller(self.surf_directions["z"]) +
+        my_ax.set_title(self._vec_to_miller(self.surf_directions["z"], brackets="()") +
                         " Gamma Surface\nSurface Separation = {:0.2f} A".format(self.surface_separation))
 
         return fig, my_ax
@@ -633,7 +658,7 @@ def plot_energy_densities(self, Es=None, ax=None, si=False):
             units = "J/m$^2$"
         else:
             mul = 1
-            units = "eV/${\AA}^2$"
+            units = "eV/$ \mathrm{A}^2$"
 
         if ax is None:
 
@@ -644,9 +669,10 @@ def plot_energy_densities(self, Es=None, ax=None, si=False):
             fig = my_ax.get_figure()
 
         my_ax.plot(np.linalg.norm(self.y_disp) * np.arange(self.ny)/(self.ny-1), self.Es[0, :] * mul)
-        my_ax.set_xlabel("Glide distance in "+ self._vec_to_miller(self.surf_directions["y"]) + " ($\AA$)")
+        my_ax.set_xlabel("Glide distance in "+ self._vec_to_miller(self.surf_directions["y"]) + " ($ \mathrm{A}$)")
         my_ax.set_ylabel(f"Energy Density ({units})")
-        title = self._vec_to_miller(self.surf_directions["z"]) + " Stacking Fault\n" + "Surface Separation = {:0.2f}".format(self.surface_separation) + "${\AA}$"
+        title = self._vec_to_miller(self.surf_directions["z"], brackets="()") + " Stacking Fault\n" + \
+                "Surface Separation = {:0.2f}".format(self.surface_separation) + "$ \mathrm{A}$"
         my_ax.set_title(title)
         return fig, my_ax
 
diff --git a/matscipy/utils.py b/matscipy/utils.py
index 0fad8d8b..287a8d62 100644
--- a/matscipy/utils.py
+++ b/matscipy/utils.py
@@ -126,7 +126,7 @@ def complete_basis(v1, v2=None, normalise=False, nmax=5, tol=1E-6):
          (Searches for vectors perpendicular to v1 with small integer indices)
         If v2 is given and is not orthogonal to v1, raises a warning
     normalise: bool
-        return an orthonormal basis, rather than just orthogonal
+        return an float orthonormal basis, rather than integer orthogonal basis
     nmax: int
         Maximum integer index for v2 search
     tol: float
@@ -136,6 +136,7 @@ def complete_basis(v1, v2=None, normalise=False, nmax=5, tol=1E-6):
     -------
     V1, V2, V3: np.arrays
         Complete orthogonal basis, optionally normalised
+        dtype of arrays is int with normalise=False, float with normalise=True
     '''
 
     def _v2_search(v1, nmax, tol):
@@ -159,8 +160,8 @@ def _v2_search(v1, nmax, tol):
 
                     if np.abs(np.dot(v1, test_vec)) < tol:
                         return test_vec
-                    # No nice integer vector found!
-                    raise RuntimeError(f"Could not automatically find an integer basis from basis vector {v1}")
+        # No nice integer vector found!
+        raise RuntimeError(f"Could not automatically find an integer basis from basis vector {v1}")
 
     V1 = np.array(v1).copy().astype(int)
 
@@ -246,4 +247,171 @@ def get_structure_types(structure, diamond_structure=False):
     colors = [structure.color for structure in modifier.structures]
     hex_colors = ['#{:02x}{:02x}{:02x}'.format(int(r*255), int(g*255), int(b*255)) for r, g, b in colors] 
 
-    return atom_labels, structure_names, hex_colors
\ No newline at end of file
+    return atom_labels, structure_names, hex_colors
+
+def line_intersect_2D(p1, p2, x1, x2):
+    '''
+    Test if 2D finite line defined by points p1 & p2 intersects with the finite line defined by x1 & x2.
+    Essentially a Python conversion of the ray casting algorithm suggested in: 
+    https://stackoverflow.com/questions/217578/how-can-i-determine-whether-a-2d-point-is-within-a-polygon
+
+    Returns True if lines intersect (including if they are collinear)
+
+    p1, p2, x1, x2: np.array
+        2D points defining lines
+    '''
+
+    # Convert from pointwise p1, p2 form to ax + by + c = 0 form
+    a_p = p2[1] - p1[1]
+    b_p = p1[0] - p2[0]
+    c_p = (p2[0] * p1[1]) - (p1[0] * p2[1])
+
+    # Sub in x1, x2
+    d1 = (a_p * x1[0]) + (b_p * x1[1]) + c_p
+    d2 = (a_p * x2[0]) + (b_p * x2[1]) + c_p
+
+    # Both points above the line
+    if (d1 > 0 and d2 > 0):
+        return False
+    # Both points below the line
+    if (d1 < 0 and d2 < 0):
+        return False
+
+    # One point above, one point below the *infinite* line defined by p1, p2
+    # Instead of testing for finite line, switch p <-> x and repeat
+    # (Both finite lines intersect if both sets of points intersect the infinite
+    # lines defined by the other set of points)
+
+    # Convert from pointwise x1, x2 form to ax + by + c = 0 form
+    a_x = x2[1] - x1[1]
+    b_x = x1[0] - x2[0]
+    c_x = (x2[0] * x1[1]) - (x1[0] * x2[1])
+
+    # Sub in p1, p2
+    d1 = (a_x * p1[0]) + (b_x * p1[1]) + c_x
+    d2 = (a_x * p2[0]) + (b_x * p2[1]) + c_x
+
+    # Both points above the line
+    if (d1 > 0 and d2 > 0):
+        return False
+    # Both points below the line
+    if (d1 < 0 and d2 < 0): 
+        return False
+
+    # 2 possibilities left: the lines intersect at a single point, or the lines are collinear,
+    # both result in an intersection, so return True
+    return True
+
+def points_in_polygon2D(p, poly_points):
+    '''
+    Test if points lies within the closed 2D polygon defined by poly_points
+    Uses ray casting algorithm, as suggested by:
+    https://stackoverflow.com/questions/217578/how-can-i-determine-whether-a-2d-point-is-within-a-polygon
+
+    p: np.array
+        (n, 2) array of points to test
+    poly_points: np.array
+        Points defining the 2D polygon
+
+    Returns
+    -------
+    mask: np.array
+        Boolean mask. True for points inside the poylgon
+    '''
+    if len(p.shape) == 1:
+        # Single point provided
+        points = p.copy()[np.newaxis, :]
+    else:
+        points = p.copy()
+
+    points += 1E-3
+
+    # Ensure polygon is closed
+    if np.any(poly_points[0, :] != poly_points[-1, :]):
+        poly_points = np.append(poly_points, poly_points[0, :][np.newaxis, :], axis=0)
+
+    npoints = points.shape[0]
+
+    mask = np.zeros(npoints, dtype=bool)
+
+    # Get random point that is definitely outside the polygon
+    a, b = 10 + np.random.random(size=2) * 10
+
+    test_point = np.array([a * np.max(poly_points[:, 0]), b * np.max(poly_points[:, 1])])
+
+    for i in range(npoints):
+        intersections = 0
+        for j in range(poly_points.shape[0] - 1):
+            intersections += line_intersect_2D(poly_points[j, :2], poly_points[j+1, :2], test_point, points[i, :2])
+        if intersections % 2:
+            # Even number of intersections, point is inside polygon
+            mask[i] = True
+    return mask.astype(bool)
+
+def get_distance_from_polygon2D(test_points:np.array, polygon_points:np.array) -> np.array:
+    '''
+    Get shortest distance between a test point and a polygon defined by polygon_points
+        (i.e. the shortest distance between each point and the lines of the polygon)
+    
+    Uses formula from https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line#Line_defined_by_two_points
+
+    test_points: np.array
+        2D Points to get distances for
+    polygon_points: np.array
+        Coordinates defining the points of the polygon
+    '''
+
+    def get_dist(p, v, w):
+        '''
+        Gets distance between point p, and the line segment defined by v and w
+        From https://stackoverflow.com/questions/849211/shortest-distance-between-a-point-and-a-line-segment
+        '''
+        denominator = np.linalg.norm(v - w)
+
+        if denominator == 0.0:
+            # point 1 and point 2 are the same, return distance between test_point and poly_point_1
+            return np.linalg.norm(p - v)
+        
+        # Use t on [0, 1] to parameterize moving on the line segment defined by poly points
+        t = np.max([0, np.min([1, np.dot(p - v, w - v)/denominator**2])])
+
+        # Closest point on line segment
+        projection = v + t * (w - v)
+
+        return np.linalg.norm(p - projection)
+        
+    distances = np.zeros(test_points.shape[0])
+
+    N = polygon_points.shape[0]
+    for i in range(test_points.shape[0]):
+        test_point = test_points[i, :]
+
+        distances[i] = min([get_dist(test_point, polygon_points[j, :], polygon_points[(j+1)%N, :]) for j in range(N)])
+    
+    return distances
+
+def radial_mask_from_polygon2D(test_points:np.array, polygon_points:np.array, radius:float, inner:bool=True) -> np.array:
+    '''
+    Get a boolean mask of all test_points within a radius of any edge of the polygon defined by polygon_points
+
+    test_points: np.array
+        2D Points to get mask for
+    polygon_points: np.array
+        Coordinates defining the points of the polygon
+    radius: float
+        Radius to use as cutoff for mask
+    inner: bool
+        Whether test_points inside the polygon should always be masked as True, regardless of radius
+    '''
+
+    distances = get_distance_from_polygon2D(test_points, polygon_points)
+
+    outer_mask = distances <= radius 
+
+    if inner:
+        inner_mask = points_in_polygon2D(test_points, polygon_points)
+
+        full_mask = (inner_mask + outer_mask).astype(bool)
+    else:
+        full_mask = outer_mask.astype(bool)
+    return full_mask
\ No newline at end of file
diff --git a/pyproject.toml b/pyproject.toml
index 748859ff..6cf0ac57 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -27,6 +27,7 @@ dependencies = [
 [project.optional-dependencies]
 test = [
     "pytest",
+    "pytest-subtests",
     "sympy",
     "atomman",
     "ovito"


Test patch:
diff --git a/tests/conftest.py b/tests/conftest.py
index 4f42425a..326c83d8 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -1,8 +1,33 @@
 import os.path
 
 import pytest
+import inspect
+import numpy as np
+
 
 @pytest.fixture
 def datafile_directory():
     """Absolute path to the `tests` directory (which also contains the data files)"""
-    return os.path.dirname(__file__)
\ No newline at end of file
+    return os.path.dirname(__file__)
+
+
+def pytest_make_parametrize_id(config, val, argname):
+    '''
+    Define how pytest should label parametrized tests
+    '''
+    
+    if inspect.isclass(val):
+        # Strip full path (e.g. matscipy.dislocation.X)
+        # And also the <class 'x'> wrapper
+        v = str(val).split(".")[-1][:-2]
+
+        return f"{argname}={v}"
+    
+    elif type(val) in [bool, str] or \
+        np.issubdtype(type(val), np.integer) or \
+        np.issubdtype(type(val), np.floating):
+        # Common rule for nice class names
+        return f"{argname}={val}"
+
+    # Allow pytest default behaviour
+    return None
\ No newline at end of file
diff --git a/tests/matscipytest.py b/tests/matscipytest.py
index 5125fca4..57314d14 100644
--- a/tests/matscipytest.py
+++ b/tests/matscipytest.py
@@ -45,13 +45,15 @@
 import logging
 from io import StringIO
 import numpy as np
+import pytest
+import pytest_subtests
 
 def string_to_array(s):
     return np.loadtxt(StringIO(s)).T
 
-class MatSciPyTestCase(unittest.TestCase):
+class MatSciPyBaseTest():
     """
-    Subclass of unittest.TestCase with extra methods for comparing arrays and dictionaries
+    Base Class of methods for comparing arrays and dictionaries
     """
 
     def assertDictionariesEqual(self, d1, d2, skip_keys=[], ignore_case=True):
@@ -123,9 +125,44 @@ def assertArrayAlmostEqual(self, a, b, tol=1e-7):
     def assertAtomsAlmostEqual(self, a, b, tol=1e-7):
         self.assertArrayAlmostEqual(a.positions, b.positions, tol)
         self.assertArrayAlmostEqual(a.numbers, b.numbers)
-        self.assertArrayAlmostEqual(a._cell, b._cell)
+        self.assertArrayAlmostEqual(a.cell[:, :], b.cell[:, :])
         self.assertArrayAlmostEqual(a._pbc, b._pbc)
 
+
+class MatSciPyTestCase(MatSciPyBaseTest, unittest.TestCase):
+    '''
+    unittest.TestCase class for unittests requiring assertAtomsAlmostEqual etc.
+    Define methods like self.fail, self.subTest, etc through unittest.TestCase
+    '''
+    pass
+
+class MatSciPyTestFixture(MatSciPyBaseTest):
+    '''
+    Pytest-only variant of test framework, which doesn't depend on unittest.TestCase
+    (and so it can be @pytest.mark.parameterize'd)
+
+
+    self.subTest calls need to be replaced with pytest-subtests style calls
+    (i.e. including subtests as a function arg, and calling subtest.test())
+    See https://blog.ganssle.io/articles/2020/04/subtests-in-python.html for more info
+    '''
+
+    @staticmethod
+    def fail(msg=""):
+        return pytest.fail(msg)
+    
+    @staticmethod
+    def skipTest(msg=""):
+        return pytest.skip(msg)
+    
+    def assertNotAlmostEqual(self, a, b, tol=1e-6, msg=""):
+        if np.max(np.abs(a - b)) < tol:
+            return self.fail(msg)
+
+    def assertTrue(self, a, msg=""):
+        if not a == True:
+            self.fail(msg)
+
 def skip(f):
     """
     Decorator which can be used to skip unit tests
diff --git a/tests/test_dislocation.py b/tests/test_dislocation.py
index a6e3a1a3..0f4928c2 100644
--- a/tests/test_dislocation.py
+++ b/tests/test_dislocation.py
@@ -22,6 +22,7 @@
 import os
 import unittest
 import matscipytest
+import pytest
 import sys
 
 import matscipy.dislocation as sd
@@ -29,6 +30,10 @@
 
 from ase.calculators.lammpslib import LAMMPSlib
 from matscipy.calculators.eam import EAM
+from matscipy.dislocation import DiamondGlide90degreePartial, get_elastic_constants
+from matscipy.calculators.manybody.explicit_forms.stillinger_weber import StillingerWeber,\
+                                                                Holland_Marder_PRL_80_746_Si
+from matscipy.calculators.manybody import Manybody
 
 from ase.build import bulk as ase_bulk
 
@@ -55,6 +60,32 @@
 except ImportError:
     print("ovito not found: skipping some tests")
 
+# Generator functions for lists of cubic dislocations
+def cubic_dislocs():
+    import inspect
+
+    # Get all cubic dislocation classes
+    # So new disloc types are automatically added to testing framework
+    _cubic_dislocs = [item for name, item in sd.__dict__.items() if inspect.isclass(item) and issubclass(item, sd.CubicCrystalDislocation)]
+    _cubic_dislocs = [item for item in _cubic_dislocs if item not in [sd.CubicCrystalDislocation, 
+                                                                    sd.CubicCrystalDissociatedDislocation,
+                                                                    sd.CubicCrystalDislocationQuadrupole]]
+    return _cubic_dislocs
+
+def cubic_perfect_dislocs():
+    _cubic_dislocs = cubic_dislocs()
+
+    _cubic_perfect_dislocs = [item for item in _cubic_dislocs if not issubclass(item, sd.CubicCrystalDissociatedDislocation)]
+    
+    return _cubic_perfect_dislocs
+
+def cubic_dissociated_dislocs():
+    _cubic_dislocs = cubic_dislocs()
+
+    _cubic_dissociated_dislocs = [item for item in _cubic_dislocs if issubclass(item, sd.CubicCrystalDissociatedDislocation)]
+    
+    return _cubic_dissociated_dislocs
+
 
 class TestDislocation(matscipytest.MatSciPyTestCase):
     """Class to store test for dislocation.py module."""
@@ -107,7 +138,7 @@ def test_elastic_constants_EAM(self):
         pot_path = os.path.join(test_dir, pot_name)
         calc_EAM = EAM(pot_path)
         obtained_values = sd.get_elastic_constants(calculator=calc_EAM,
-                                                   delta=1.0e-3)
+                                                   delta=1.0e-3, verbose=False)
 
         self.assertArrayAlmostEqual(obtained_values, target_values, tol=1e-2)
 
@@ -136,7 +167,7 @@ def test_elastic_constants_lammpslib(self):
                            log_file="lammps.log")
 
         obtained_values = sd.get_elastic_constants(calculator=lammps,
-                                                   delta=1.0e-3)
+                                                   delta=1.0e-3, verbose=False)
 
         os.remove("lammps.log")
 
@@ -219,7 +250,7 @@ def test_differential_displacement(self):
         fig = sd.show_NEB_configurations([disloc_ini, disloc_fin], bulk_ini,
                                          xyscale=5.0, show=False)
         print("'dd_test.png' will be created: check the displacement map")
-        fig.savefig("dd_test.png")
+        #fig.savefig("dd_test.png")
 
     @unittest.skipIf("atomman" not in sys.modules,
                      'requires Stroh solution from atomman to run')
@@ -411,272 +442,6 @@ def test_slice_long_dislo(self):
         # to length of kink_length * 3 - 1 (for left kink)
         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, gen_bulk=False):
-        alat = 3.14339177996466
-        C11 = 523.0266819809012
-        C12 = 202.1786296941397
-        C44 = 160.88179872237012
-
-        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
-        assert np.unique(bulk.get_chemical_symbols()) == "Al"
-        assert np.unique(disloc.get_chemical_symbols()) == "Al"
-
-        assert len(bulk) == len(disloc)
-
-        if test_u:
-            # test the consistency
-            # displacement = disloc.positions - bulk.positions
-            stroh_displacement = d.displacements(bulk.positions,
-                                                 np.diag(bulk.cell) / 2.0,
-                                                 self_consistent=d.self_consistent)
-
-            displacement = disloc.positions - bulk.positions
-
-            np.testing.assert_array_almost_equal(displacement,
-                                                 stroh_displacement)
-
-        results = sd.ovito_dxa_straight_dislo_info(disloc, structure=structure)
-        assert len(results) == 1
-        position, b, line, angle = results[0]
-        self.assertArrayAlmostEqual(np.abs(b), burgers)  # signs can change
-
-        err = angle - ref_angle
-        print(f'angle = {angle} ref_angle = {ref_angle} err = {err}')
-        assert abs(err) < tol
-
-    @unittest.skipIf("atomman" not in sys.modules or
-                     "ovito" not in sys.modules,
-                     "requires atomman and ovito")
-    def test_screw_dislocation(self):
-        self.check_disloc(sd.BCCScrew111Dislocation, 0.0)
-
-    @unittest.skipIf("atomman" not in sys.modules or
-                     "ovito" not in sys.modules,
-                     "requires atomman and ovito")
-    def test_edge_dislocation(self):
-        self.check_disloc(sd.BCCEdge111Dislocation, 90.0)
-
-    @unittest.skipIf("atomman" not in sys.modules or
-                     "ovito" not in sys.modules,
-                     "requires atomman and ovito")
-    def test_edge111bar_dislocation(self):
-        self.check_disloc(sd.BCCEdge111barDislocation, 90.0)
-
-    @unittest.skipIf("atomman" not in sys.modules or
-                     "ovito" not in sys.modules,
-                     "requires atomman and ovito")
-    def test_edge100_dislocation(self,):
-        self.check_disloc(sd.BCCEdge100Dislocation, 90.0,
-                          burgers=np.array([1.0, 0.0, 0.0]))
-
-    @unittest.skipIf("atomman" not in sys.modules or
-                     "ovito" not in sys.modules,
-                     "requires atomman and ovito")
-    def test_edge100110_dislocation(self,):
-        self.check_disloc(sd.BCCEdge100110Dislocation, 90.0,
-                          burgers=np.array([1.0, 0.0, 0.0]))
-
-    @unittest.skipIf("atomman" not in sys.modules or
-                     "ovito" not in sys.modules,
-                     "requires atomman and ovito")
-    def test_mixed_dislocation(self):
-        self.check_disloc(sd.BCCMixed111Dislocation, 70.5)
-
-    @unittest.skipIf("atomman" not in sys.modules or
-                     "ovito" not in sys.modules,
-                     "requires atomman and ovito")
-    def test_30degree_diamond_partial_dislocation(self,):
-        self.check_disloc(sd.DiamondGlide30degreePartial, 30.0,
-                          structure="Diamond",
-                          burgers=(1.0 / 6.0) * np.array([1.0, 2.0, 1.0]))
-
-    @unittest.skipIf("atomman" not in sys.modules or
-                     "ovito" not in sys.modules,
-                     "requires atomman and ovito")
-    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,
-                     "requires atomman and ovito")
-    def test_screw_diamond_dislocation(self,):
-        self.check_disloc(sd.DiamondGlideScrew, 0.0,
-                          structure="Diamond",
-                          burgers=(1.0 / 2.0) * np.array([0.0, 1.0, 1.0]))
-
-    @unittest.skipIf("atomman" not in sys.modules or
-                     "ovito" not in sys.modules,
-                     "requires atomman and ovito")
-    def test_60degree_diamond_dislocation(self,):
-        self.check_disloc(sd.DiamondGlide60Degree, 60.0,
-                          structure="Diamond",
-                          burgers=(1.0 / 2.0) * np.array([1.0, 0.0, 1.0]))
-
-    @unittest.skipIf("atomman" not in sys.modules or
-                     "ovito" not in sys.modules,
-                     "requires atomman and ovito")
-    def test_fcc_screw_shockley_partial_dislocation(self,):
-        self.check_disloc(sd.FCCScrewShockleyPartial, 30.0,
-                          structure="FCC",
-                          burgers=(1.0 / 6.0) * np.array([1.0, 2.0, 1.0]))
-
-    @unittest.skipIf("atomman" not in sys.modules or
-                     "ovito" not in sys.modules,
-                     "requires atomman and ovito")
-    def test_fcc_screw110_dislocation(self,):
-        self.check_disloc(sd.FCCScrew110Dislocation, 0.0,
-                          structure="FCC",
-                          burgers=(1.0 / 2.0) * np.array([0.0, 1.0, 1.0]))
-
-    @unittest.skipIf("atomman" not in sys.modules or
-                     "ovito" not in sys.modules,
-                     "requires atomman and ovito")
-    def test_fcc_edge_shockley_partial_dislocation(self,):
-        self.check_disloc(sd.FCCEdgeShockleyPartial, 60.0,
-                          structure="FCC",
-                          burgers=(1.0 / 6.0) * np.array([1.0, 2.0, 1.0]))
-
-    @unittest.skipIf("atomman" not in sys.modules or
-                     "ovito" not in sys.modules,
-                     "requires atomman and ovito")
-    def test_fcc_edge110_partial_dislocation(self,):
-        self.check_disloc(sd.FCCEdge110Dislocation, 90.0,
-                          structure="FCC",
-                          burgers=(1.0 / 2.0) * np.array([1.0, 1.0, 0.0]))
-
-    def check_glide_configs(self, cls, structure="BCC"):
-        alat = 3.14339177996466
-        C11 = 523.0266819809012
-        C12 = 202.1786296941397
-        C44 = 160.88179872237012
-
-        d = cls(alat, C11, C12, C44)
-        bulk, disloc_ini, disloc_fin = d.build_glide_configurations(radius=40)
-
-        assert len(bulk) == len(disloc_ini)
-        assert len(disloc_ini) == len(disloc_fin)
-
-        assert all(disloc_ini.get_array("fix_mask") ==
-                   disloc_fin.get_array("fix_mask"))
-
-        results = sd.ovito_dxa_straight_dislo_info(disloc_ini,
-                                                   structure=structure)
-        assert len(results) == 1
-        ini_x_position = results[0][0][0]
-
-        results = sd.ovito_dxa_straight_dislo_info(disloc_fin,
-                                                   structure=structure)
-        assert len(results) == 1
-        fin_x_position = results[0][0][0]
-        # test that difference between initial and final positions are
-        # roughly equal to glide distance.
-        # Since the configurations are unrelaxed dxa gives
-        # a very rough estimation (especially for edge dislocations)
-        # thus tolerance is taken to be ~1 Angstrom
-        np.testing.assert_almost_equal(fin_x_position - ini_x_position,
-                                       d.glide_distance, decimal=0)
-
-    @unittest.skipIf("atomman" not in sys.modules or
-                     "ovito" not in sys.modules,
-                     "requires atomman and ovito")
-    def test_screw_glide(self):
-        self.check_glide_configs(sd.BCCScrew111Dislocation)
-
-    @unittest.skipIf("atomman" not in sys.modules or
-                     "ovito" not in sys.modules,
-                     "requires atomman and ovito")
-    def test_edge_glide(self):
-        self.check_glide_configs(sd.BCCEdge111Dislocation)
-
-    @unittest.skipIf("atomman" not in sys.modules or
-                     "ovito" not in sys.modules,
-                     "requires atomman and ovito")
-    def test_mixed_glide(self):
-        self.check_glide_configs(sd.BCCMixed111Dislocation)
-
-    @unittest.skipIf("atomman" not in sys.modules or
-                     "ovito" not in sys.modules,
-                     "requires atomman and ovito")
-    def test_edge100_glide(self):
-        self.check_glide_configs(sd.BCCEdge100Dislocation)
-
-    @unittest.skipIf("atomman" not in sys.modules or
-                     "ovito" not in sys.modules,
-                     "requires atomman and ovito")
-    def test_edge100110_glide(self):
-        self.check_glide_configs(sd.BCCEdge100110Dislocation)
-
-    @unittest.skipIf("atomman" not in sys.modules or
-                     "ovito" not in sys.modules,
-                     "requires atomman and ovito")
-    def test_30degree_diamond_partial_glide(self):
-        self.check_glide_configs(sd.DiamondGlide30degreePartial,
-                                 structure="Diamond")
-
-    @unittest.skipIf("atomman" not in sys.modules or
-                     "ovito" not in sys.modules,
-                     "requires atomman and ovito")
-    def test_90degree_diamond_partial_glide(self):
-        self.check_glide_configs(sd.DiamondGlide90degreePartial,
-                                 structure="Diamond")
-
-    @unittest.skipIf("atomman" not in sys.modules or
-                     "ovito" not in sys.modules,
-                     "requires atomman and ovito")
-    def test_screw_diamond_glide(self):
-        self.check_glide_configs(sd.DiamondGlideScrew,
-                                 structure="Diamond")
-
-    @unittest.skipIf("atomman" not in sys.modules or
-                     "ovito" not in sys.modules,
-                     "requires atomman and ovito")
-    def test_60degree_diamond_glide(self):
-        self.check_glide_configs(sd.DiamondGlide60Degree,
-                                 structure="Diamond")
-
-    @unittest.skipIf("atomman" not in sys.modules or
-                     "ovito" not in sys.modules,
-                     "requires atomman and ovito")
-    def test_fcc_screw_shockley_partial_glide(self):
-        self.check_glide_configs(sd.FCCScrewShockleyPartial,
-                                 structure="FCC")
-
-    @unittest.skipIf("atomman" not in sys.modules or
-                     "ovito" not in sys.modules,
-                     "requires atomman and ovito")
-    def test_fcc_screw_110_glide(self):
-        self.check_glide_configs(sd.FCCScrew110Dislocation,
-                                 structure="FCC")
-
-    @unittest.skipIf("atomman" not in sys.modules or
-                     "ovito" not in sys.modules,
-                     "requires atomman and ovito")
-    def test_fcc_edge_shockley_partial_glide(self):
-        self.check_glide_configs(sd.FCCEdgeShockleyPartial,
-                                 structure="FCC")
-
-    @unittest.skipIf("atomman" not in sys.modules or
-                     "ovito" not in sys.modules,
-                     "requires atomman and ovito")
-    def test_fcc_edge_dislocation_glide(self):
-        self.check_glide_configs(sd.FCCEdge110Dislocation,
-                                 structure="FCC")
-
     def test_fixed_line_atoms(self):
 
         from ase.build import bulk
@@ -760,12 +525,12 @@ def check_anisotropic_disloc(self, cls, ref_angle, structure="BCC",
 
         # Setup CubicCrystalDislocation object, and get displacments
         ccd = cls(a0, C11, C12, C44, symbol="Fe")
-        bulk, sc_disloc1 = ccd.build_cylinder(20.0)
+        bulk, sc_disloc1 = ccd.build_cylinder(20.0, verbose=False)
         center = np.diag(bulk.cell) / 2
-        stroh_disp = ccd.displacements(bulk.positions, center, use_atomman=True, self_consistent=False)
+        stroh_disp = ccd.displacements(bulk.positions, center, method="atomman", self_consistent=False, verbose=False)
 
         # Get displacements using AnistoropicDislocation class
-        adsl_disp = ccd.displacements(bulk.positions, center, use_atomman=False, self_consistent=False)
+        adsl_disp = ccd.displacements(bulk.positions, center, method="adsl", self_consistent=False, verbose=False)
 
         # Setup the dislcation from the AnistoropicDislocation object
         disloc = bulk.copy()
@@ -814,48 +579,41 @@ def check_anisotropic_disloc(self, cls, ref_angle, structure="BCC",
             grad2D_stroh = np.transpose(grad2D_stroh_T)
 
             # Find 2D gradient tensor from AnistoropicDislocation object
-            grad2D_adsl = ccd.ADstroh.deformation_gradient(bulk, center)
+            grad2D_adsl = ccd.ADstroh.deformation_gradient(bulk.positions, center)
 
             # Check gradients
             np.testing.assert_array_almost_equal(grad2D_adsl, grad2D_stroh)
         
-
-    @unittest.skip('Reenable in #206')
     @unittest.skipIf("atomman" not in sys.modules or
                      "ovito" not in sys.modules,
                      "requires atomman and ovito")
     def test_screw_dislocation_adsl(self):
         self.check_anisotropic_disloc(sd.BCCScrew111Dislocation, 0.0)
 
-    @unittest.skip('Reenable in #206')
     @unittest.skipIf("atomman" not in sys.modules or
                      "ovito" not in sys.modules,
                      "requires atomman and ovito")
     def test_edge_dislocation_adsl(self):
         self.check_anisotropic_disloc(sd.BCCEdge111Dislocation, 90.0)
 
-    @unittest.skip('Reenable in #206')
     @unittest.skipIf("atomman" not in sys.modules or
                      "ovito" not in sys.modules,
                      "requires atomman and ovito")
     def test_edge111bar_dislocation_adsl(self):
         self.check_anisotropic_disloc(sd.BCCEdge111barDislocation, 90.0)
 
-    @unittest.skip('Reenable in #206')
     @unittest.skipIf("atomman" not in sys.modules or
                      "ovito" not in sys.modules,
                      "requires atomman and ovito")
     def test_edge100_dislocation_adsl(self,):
         self.check_anisotropic_disloc(sd.BCCEdge100Dislocation, 90.0)
 
-    @unittest.skip('Reenable in #206')
     @unittest.skipIf("atomman" not in sys.modules or
                      "ovito" not in sys.modules,
                      "requires atomman and ovito")
     def test_edge100110_dislocation_adsl(self,):
         self.check_anisotropic_disloc(sd.BCCEdge100110Dislocation, 90.0)
 
-    @unittest.skip('Reenable in #206')
     @unittest.skipIf("atomman" not in sys.modules or
                      "ovito" not in sys.modules,
                      "requires atomman and ovito")
@@ -863,5 +621,331 @@ def test_mixed_dislocation_adsl(self):
         self.check_anisotropic_disloc(sd.BCCMixed111Dislocation, 70.5)
 
 
+test_props = {
+    # Calculators for several bulk properties
+    "diamond" : {
+        "symbol" : "Si",
+        "props" : sd.get_elastic_constants(calculator=Manybody(**StillingerWeber(Holland_Marder_PRL_80_746_Si)), 
+                                           symbol="Si", verbose=False)
+    },
+
+    "fcc" : {
+        "symbol" : "Ni",
+        "props" : sd.get_elastic_constants(calculator=EAM(test_dir + os.sep + "FeCuNi.eam.alloy"), 
+                                           symbol="Ni", verbose=False)
+    },
+
+    "bcc" : {
+        "symbol" : "W",
+        "props" : sd.get_elastic_constants(calculator=EAM(test_dir + os.sep + "w_eam4.fs"), 
+                                           symbol="W", verbose=False)
+    }
+}
+
+class BaseTestCubicCrystalDislocation(matscipytest.MatSciPyTestFixture):
+
+    has_atomman = "atomman" in sys.modules
+    has_ovito = "ovito" in sys.modules
+    
+    def set_up_cls(self, disloc_cls):
+        self.test_cls = disloc_cls
+
+        self.structure = self.test_cls.crystalstructure
+
+        self.symbol = test_props[self.structure.lower()]["symbol"]
+
+        self.alat, self.C11, self.C12, self.C44 = test_props[self.structure.lower()]["props"]
+
+        self.default_method = "atomman" if self.has_atomman else "adsl"
+
+
+        if issubclass(self.test_cls, sd.CubicCrystalDissociatedDislocation):
+            self.ncores = 2
+        else:
+            self.ncores = 1
+
+    @pytest.mark.parametrize("gen_bulk", [True, False])
+    def test_dislocation_cylinder(self, disloc, gen_bulk, 
+                                  test_u=True, tol=10.0):
+        '''
+        Test construction of CubicCrystalDislocation.build_cylinder
+        '''
+        self.set_up_cls(disloc)
+
+        if gen_bulk:
+            a = ase_bulk(self.symbol, self.structure.lower(), self.alat, cubic=True)
+        else:
+            a = self.alat
+
+        d = self.test_cls(a, self.C11, self.C12, self.C44, symbol=self.symbol)
+        bulk, disloc = d.build_cylinder(20.0, method=self.default_method, 
+                                        self_consistent=d.self_consistent, verbose=False)
+
+        # test that assigning non default symbol worked
+        assert np.unique(bulk.get_chemical_symbols()) == self.symbol
+        assert np.unique(disloc.get_chemical_symbols()) == self.symbol
+
+        assert len(bulk) == len(disloc)
+
+        if test_u:
+            # test the consistency
+            # displacement = disloc.positions - bulk.positions
+            stroh_displacement = d.displacements(bulk.positions,
+                                                 np.array(disloc.info["core_positions"]*self.ncores),
+                                                 self_consistent=d.self_consistent,
+                                                 verbose=False,
+                                                 method=self.default_method)
+
+            displacement = disloc.positions - bulk.positions
+
+            np.testing.assert_array_almost_equal(displacement,
+                                                 stroh_displacement)
+    
+    @pytest.mark.parametrize("gen_bulk", [True, False])
+    def test_cylinder_ovito_dxa(self, disloc, gen_bulk, 
+                                test_u=True, tol=10.0):
+        self.set_up_cls(disloc)
+            
+        if not self.has_ovito:
+            self.skipTest("ovito module not installed")
+
+        if gen_bulk:
+            a = ase_bulk(self.symbol, self.structure.lower(), self.alat, cubic=True)
+        else:
+            a = self.alat
+
+        d = self.test_cls(a, self.C11, self.C12, self.C44, symbol=self.symbol)
+        bulk, disloc = d.build_cylinder(20.0, method=self.default_method, verbose=False)
+
+        results = sd.ovito_dxa_straight_dislo_info(disloc, structure=self.structure)
+
+        assert len(results) == 1
+        position, b, line, angle = results[0]
+
+        # Sign flips and symmetry groups can cause differences in direct comparisons
+        self.assertArrayAlmostEqual(np.sort(np.abs(b)), np.sort(np.abs(d.burgers_dimensionless)))
+
+        # norm_burgers = d.burgers_dimensionless / np.linalg.norm(d.burgers_dimensionless)
+        # norm_axis = d.axes[:, 2] / np.linalg.norm(d.axes[:, 2])
+        # ref_angle = np.arccos(np.dot(norm_burgers, norm_axis))
+        # err = angle - ref_angle
+        # print(f'angle = {angle} ref_angle = {ref_angle} err = {err}')
+        # assert abs(err) < tol
+        
+    def test_glide_configs(self, disloc, subtests):  
+        self.set_up_cls(disloc)      
+
+        d = self.test_cls(self.alat, self.C11, self.C12, self.C44, symbol=self.symbol)
+        bulk, disloc_ini, disloc_fin = d.build_glide_configurations(radius=40, method=self.default_method, verbose=False)
+
+        assert len(bulk) == len(disloc_ini)
+        assert len(disloc_ini) == len(disloc_fin)
+
+        assert all(disloc_ini.get_array("fix_mask") ==
+                   disloc_fin.get_array("fix_mask"))
+
+        with subtests.test("Check glide configs against Ovito DXA"):
+            if not self.has_ovito:
+                self.skipTest("ovito module not installed")
+
+            results = sd.ovito_dxa_straight_dislo_info(disloc_ini,
+                                                    structure=self.structure)
+            assert len(results) == 1
+            ini_x_position = results[0][0][0]
+
+            results = sd.ovito_dxa_straight_dislo_info(disloc_fin,
+                                                    structure=self.structure)
+            assert len(results) == 1
+            fin_x_position = results[0][0][0]
+            # test that difference between initial and final positions are
+            # roughly equal to glide distance.
+            # Since the configurations are unrelaxed dxa gives
+            # a very rough estimation (especially for edge dislocations)
+            # thus tolerance is taken to be ~1 Angstrom
+            np.testing.assert_almost_equal(fin_x_position - ini_x_position,
+                                        d.glide_distance, decimal=0)
+        
+    def test_methods(self, disloc, subtests):
+        self.set_up_cls(disloc)
+
+        methods = self.test_cls.avail_methods
+
+        d = self.test_cls(self.alat, self.C11, self.C12, self.C44, symbol=self.symbol)
+        bulk, _ = d.build_cylinder(radius=20.0, verbose=False)
+
+        # Base method to compare against
+        base_method=self.default_method
+
+        methods = [method for method in methods if method != base_method]
+
+        core_positions = np.zeros((self.ncores, 3))
+
+        base_displacements = d.displacements(bulk.positions, core_positions, method=base_method, verbose=False)
+
+        for i, method in enumerate(methods):
+            
+            with subtests.test(f"Validating {method} method against {base_method}", i=i):
+                if method == "atomman" and not self.has_atomman:
+                    self.skipTest("atomman not installed")
+                method_disps = d.displacements(bulk.positions, core_positions, method=method, verbose=False)
+
+                try:
+                    np.testing.assert_array_almost_equal(method_disps, base_displacements)
+                except AssertionError as e:
+
+                    raise AssertionError(f"Displacements from {method} did not match {base_method}")
+
+@pytest.mark.parametrize("disloc", cubic_perfect_dislocs())
+class TestCubicCrystalDislocation(BaseTestCubicCrystalDislocation):
+    pass
+
+@pytest.mark.parametrize("disloc", cubic_dissociated_dislocs())
+class TestCubicCrystalDissociatedDislocation(BaseTestCubicCrystalDislocation):
+    @pytest.mark.parametrize("gen_bulk", [True, False])
+    def test_dissociated_disloc(self, disloc, gen_bulk):
+        self.set_up_cls(disloc)
+
+        if gen_bulk:
+            a = ase_bulk(self.symbol, self.structure.lower(), self.alat, cubic=True)
+        else:
+            a = self.alat
+
+        d = self.test_cls(a, self.C11, self.C12, self.C44, symbol=self.symbol)
+        
+        partial_dist = 5.0        
+        bulk, disloc = d.build_cylinder(40.0, partial_distance=partial_dist, verbose=False)
+
+        # Check good agreement of displacements method with sum of 
+        # left and right displacements
+        # May not be equal, due to self_consistent
+
+        core_positions = np.array([
+            [0, 0, 0],
+            [partial_dist, 0, 0]
+        ])
+        
+        # Left + Right displacements
+        full_disps = d.displacements(bulk.positions, core_positions, 
+                                      self_consistent=d.self_consistent, verbose=False)
+        # Left only
+        left_disps = d.left_dislocation.displacements(bulk.positions, core_positions[0, :], 
+                                      self_consistent=d.self_consistent, verbose=False)
+        # Right only
+        right_disps = d.right_dislocation.displacements(bulk.positions, core_positions[1, :], 
+                                      self_consistent=d.self_consistent, verbose=False)
+        
+        diff = left_disps + right_disps - full_disps
+
+        max_err = np.max(np.abs(diff))
+
+        # Percentage tolerance of the error, w/ respect to the lattice parameter
+        percent_tol = 1.0
+
+        assert (max_err / self.alat) <= (percent_tol / 100.0)
+
+class BaseTestCubicCrystalDislocationQuadrupole(matscipytest.MatSciPyTestFixture):
+    has_atomman = "atomman" in sys.modules
+    has_ovito = "ovito" in sys.modules
+    
+    def set_up_cls(self, disloc):
+        self.test_cls = disloc
+
+        self.structure = self.test_cls.crystalstructure
+
+        self.symbol = test_props[self.structure.lower()]["symbol"]
+
+        self.alat, self.C11, self.C12, self.C44 = test_props[self.structure.lower()]["props"]
+
+        self.default_method = "atomman" if self.has_atomman else "adsl"
+
+        if issubclass(self.test_cls, sd.CubicCrystalDissociatedDislocation):
+            self.ncores = 2
+        else:
+            self.ncores = 1
+
+    def test_quadrupole_struct(self, disloc, subtests):
+        '''
+        Validation that build_quadrupole runs without errors
+        '''
+        self.set_up_cls(disloc)
+
+        d = sd.Quadrupole(self.test_cls, self.alat, self.C11, self.C12, self.C44, symbol=self.symbol)
+
+        with subtests.test("Check build_quadrupole basic functionality"):
+
+            bulk, quad = d.build_quadrupole(glide_separation=4, verbose=False)
+
+            self.assertEqual(len(bulk), len(quad))
+
+        with subtests.test("Test build_quadrupole offset & extension args"):
+            bulk, quad = d.build_quadrupole(glide_separation=4, extension=2,
+                                            left_offset=np.array([d.glide_distance/2, 0, 0]),
+                                            verbose=False)
+
+
+    def test_glide_configs(self, disloc, subtests):
+        self.set_up_cls(disloc)
+
+        d = sd.Quadrupole(self.test_cls, self.alat, self.C11, self.C12, self.C44, symbol=self.symbol)
+
+
+        with subtests.test("Check build_glide_quadrupoles basic functionality"):
+            configs = d.build_glide_quadrupoles(nims=2, glide_separation=4, glide_left=True, glide_right=False,
+                                                verbose=False)
+
+        with subtests.test("Check build_glide_quadrupoles against equivalent build_quadrupoles calls"):
+            bulk, ini_quad = d.build_quadrupole(glide_separation=4,
+                                            verbose=False)
+            
+            bulk, fin_quad = d.build_quadrupole(glide_separation=4,
+                                            left_offset=np.array([d.glide_distance, 0, 0]),
+                                            verbose=False)
+            
+            self.assertAtomsAlmostEqual(configs[0], ini_quad)
+            self.assertAtomsAlmostEqual(configs[1], fin_quad)
+
+    def test_single_kink_quadrupole(self, disloc):
+        '''
+        Validate that generation of single kink structures runs without errors
+        '''
+        self.set_up_cls(disloc)
+        d = sd.Quadrupole(self.test_cls, self.alat, self.C11, self.C12, self.C44, symbol=self.symbol)
+
+
+        kink_struct = d.build_kink_quadrupole(z_reps=4, glide_separation=5)
+
+        # Check that the z vector has been tilted to get an infinite array of periodic single kinks
+        z_vec = kink_struct.cell[2, :]
+
+        vals = np.abs(z_vec[:2])
+
+        self.assertNotAlmostEqual(vals[0], 0.0)
+
+
+
+@pytest.mark.parametrize("disloc", cubic_perfect_dislocs())
+class TestCubicCrystalDislocationQuadrupole(BaseTestCubicCrystalDislocationQuadrupole):
+    pass
+
+
+@pytest.mark.parametrize("disloc", cubic_dissociated_dislocs())
+class TestCubicCrystalDissociatedDislocationQuadrupole(BaseTestCubicCrystalDislocationQuadrupole):
+    
+    def test_dissociated_quadrupole(self, disloc):
+        '''
+        Check execution with no errors
+        '''
+        self.set_up_cls(disloc)
+
+        d = sd.Quadrupole(self.test_cls, self.alat, self.C11, self.C12, self.C44, symbol=self.symbol)
+
+
+        # Cell has to be quite big to also have partial distances in there
+        bulk, quad = d.build_quadrupole(glide_separation=8, partial_distance=2, verbose=False)
+
+
+
+
+
 if __name__ == '__main__':
-    unittest.main()
+    unittest.main()
\ No newline at end of file
diff --git a/tests/test_utils.py b/tests/test_utils.py
index 5ca80626..ace5bad3 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -10,81 +10,271 @@
 test_dir = os.path.dirname(os.path.realpath(__file__))
 
 
-class TestDislocation(matscipytest.MatSciPyTestCase):
-    """Class to store test for utils.py module."""
+class TestUtils(matscipytest.MatSciPyTestCase):
+    """Class to store unittests of utils.py module."""
 
     def test_validate_cubic_cell(self):
-            from ase.build import bulk as build_bulk
-            import warnings
+        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"])
+        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")
 
-            # 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")
+        warnings.simplefilter("default")
 
-            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")
 
+    def test_complete_basis(self):
+
+        def test_case(vecs, check=True):
+            tol = 1e-6
+
+            v1, v2, v3 = utils_mod.complete_basis(*vecs, normalise=True)
             
-            # Test w/ supercells of bulk
-            utils_mod.validate_cubic_cell(bulk_fcc * (2, 2, 2), symbol="Cu", axes=axes, crystalstructure="fcc")
+            if check:
+                # Check orthogonal basis
+                assert v1.T @ v2 < tol
+                assert v2.T @ v3 < tol
+                assert v3.T @ v1 < tol
+
+                # Check norms of vecs
+                assert np.linalg.norm(v1) - 1 < tol
+                assert np.linalg.norm(v2) - 1 < tol
+                assert np.linalg.norm(v3) - 1 < tol
+
+                # Check handedness
+                assert np.cross(v1, v2).T @ v3 > 0
+
+                # Check dtype
+                assert v1.dtype == v2.dtype == v3.dtype
+                assert np.issubdtype(v1.dtype, np.floating)
+
+            v1, v2, v3 = utils_mod.complete_basis(*vecs, normalise=False)
+
+            if check:
+                # Check dtype for non-normalised mode
+                assert v1.dtype == v2.dtype == v3.dtype
+                assert np.issubdtype(v1.dtype, np.integer)
+
+
+        v1 = np.array([1, 1, 1])
+        v2 = np.array([1, 1, -2])
+
+        # Single vec given
+        test_case([v1])
+
+        # Two vecs
+        test_case([v1, v2])
+
+        v1 = np.array([1, 0, 0])
+        
+        # Non-orthoganal pair
+        # Will fail assertions, so bypass
+        self.assertWarns(RuntimeWarning, test_case, [v1, v2], check=False)
+
+
+        v1 = np.array([13, 19, 37])
+        # v1 is three prime numbers, unlikely to have a "neat" orthogonal vector
+        # with max value less than the default nmax
+        self.assertRaises(RuntimeError, test_case, [v1])
+
+    def test_line_intersect_2d(self):
+
+        p1 = np.array([0, 0])
+        p2 = np.array([2, 0])
+
+        # "Plus" pattern
+        x1 = np.array([1, 1])
+        x2 = np.array([1, -1])
+
+        assert utils_mod.line_intersect_2D(p1, p2, x1, x2) == True
+
+        # Lines don't cross ever
+        x2 = np.array([2, 1])
+        assert utils_mod.line_intersect_2D(p1, p2, x1, x2) == False
+
+        # Finite lines don't meet, but infinite ones would
+        x2 = np.array([2, 0.5])
+        assert utils_mod.line_intersect_2D(p1, p2, x1, x2) == False
+
+        # Lines meet at one of the points
+        x1 = np.array([0, 1])
+        x2 = np.array([0, -1])
+        assert utils_mod.line_intersect_2D(p1, p2, x1, x2) == True
+
+        # Lines end at the same point
+        x2 = p1
+        assert utils_mod.line_intersect_2D(p1, p2, x1, x2) == True
+
+
+    def test_points_in_polygon2D(self):
+
+        # Simple square shape
+        polygon = np.array([
+            [0, 0],
+            [1, 0],
+            [1, 1],
+            [0, 1]
+        ])
+
+        inside = np.array([0.5, 0.5])
+        outside = np.array([0.5, 4])
+
+        assert utils_mod.points_in_polygon2D(inside, polygon) == True
+        assert utils_mod.points_in_polygon2D(outside, polygon) == False
+
+        # Complex polygon
+
+        polygon = np.array([
+            [0, 0],
+            [5, 0],
+            [5, 5],
+            [3, 5],
+            [3, 4],
+            [4, 4],
+            [4, 1],
+            [1, 1],
+            [1, 4],
+            [2, 4],
+            [2, 5],
+            [0, 5]
+        ])
+
+        inside = np.array([
+            [0.5, 0.5],
+            [1.5, 4.5]
+        ])
+
+
+        outside = np.array([
+            [6, 0],
+            [0, 6],
+            [2.5, 4.5],
+            [2.5, 2.5],
+            [3.5, 3.5]
+        ])
+
+        assert np.all(utils_mod.points_in_polygon2D(inside, polygon) == True)
+        assert np.all(utils_mod.points_in_polygon2D(outside, polygon) == False)
+
+    def test_get_distance_from_polygon(self):
+        polygon = np.array([
+            [0, 0],
+            [1, 0],
+            [1, 1],
+            [0, 1]
+        ])
+
+        points = np.array([
+            [2, 0], # 1 away from point
+            [0.5, 0], # touching an edge of the polygon
+            [0.5, 0.5] # In the centre, 0.5 from a line
+        ])
+
+        expected_dists = np.array([1, 0, 0.5])
+
+        dists = utils_mod.get_distance_from_polygon2D(points, polygon)
+
+        assert np.allclose(dists, expected_dists)
+
+    def test_radial_mask_from_polygon2D(self):
+        polygon = np.array([
+            [0, 0],
+            [1, 0],
+            [1, 1],
+            [0, 1]
+        ])
+        
+        # Get 100 random points
+        points = np.random.rand(100, 2)
+
+        # Check that points_in_polygon2D is reproduced
+        assert np.all(utils_mod.radial_mask_from_polygon2D(points, polygon, radius=0.0, inner=True) ==
+                      utils_mod.points_in_polygon2D(points, polygon))
+        
+
+        radius = 0.2
+        inner = False
+
+        points = np.array([
+            [0.5, 0.5], # Inside polygon, outside radius from points/lines
+            [0.5, 0.1], # Close to a line
+            [2, 0.5] # Outside polygon, far from lines
+        ])
+
+        target_mask = np.array([False, True, False])
 
-            # Test w/ more complex crystal
-            nacl = build_bulk("NaCl", crystalstructure="zincblende", a=5.8, cubic=True)
+        assert np.all(utils_mod.radial_mask_from_polygon2D(points, polygon, radius=radius, 
+                                                           inner=inner) == target_mask)
+        
+        inner = True
+        target_mask = np.array([True, True, False])
+        assert np.all(utils_mod.radial_mask_from_polygon2D(points, polygon, radius=radius, 
+                                                           inner=inner) == target_mask)
 
-            utils_mod.validate_cubic_cell(nacl, symbol="Cu", axes=axes, crystalstructure="diamond")
 
-            warnings.simplefilter("default")
 
 if __name__ == '__main__':
     unittest.main()