From a37109bcd109dd614d763c9728b1e09fef847ab2 Mon Sep 17 00:00:00 2001
From: thomas-rocke <tom@rocke.co.uk>
Date: Fri, 10 May 2024 15:18:07 +0100
Subject: [PATCH 01/46] MAINT: Tweaks to disloc kink + glide

---
 matscipy/dislocation.py | 152 +++++++++++++++++++++++++++++++---------
 1 file changed, 120 insertions(+), 32 deletions(-)

diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index 1565dbb6..f642cd72 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -2587,7 +2587,7 @@ def _build_supercell(self, targetx, targety):
         return sup
 
     def _build_bulk_cyl(self, radius, core_positions, fix_rad, extension,
-                        self_consistent, method, verbose):
+                        self_consistent, method, verbose, cyl_mask=None):
         '''
         Build bulk cylinder config from args supplied by self.build_cylinder
 
@@ -2607,6 +2607,8 @@ def _build_bulk_cyl(self, radius, core_positions, fix_rad, extension,
             Will add a fictitious core at position 
             core_positions[i, :] + extension[i, :], for the 
             purposes of adding extra atoms
+        cyl_mask: array of bool
+            Optional ovverride for atomic mask to convert supercell into cyl
 
         Returns
         -------
@@ -2671,12 +2673,15 @@ def _build_bulk_cyl(self, radius, core_positions, fix_rad, extension,
 
         exts = exts[nonzero_exts, :]
 
-        mask_positions = np.vstack([new_core_positions, exts + shift])
+        if cyl_mask is None:
+            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)
+            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)
+        else:
+            mask = cyl_mask
 
         cyl = sup[mask]
 
@@ -2687,14 +2692,14 @@ def _build_bulk_cyl(self, radius, core_positions, fix_rad, extension,
                                                self_consistent=self_consistent, method=method,
                                                verbose=verbose)
 
+        fix_mask = ~radial_mask_from_polygon2D(cyl.get_positions(), mask_positions, radius - fix_rad, inner=True)
+        
         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
+        return cyl, disloc, new_core_positions, cyl_mask, fix_mask
 
     def plot_unit_cell(self, ms=250, ax=None):
         import matplotlib.pyplot as plt
@@ -2839,12 +2844,17 @@ def displacements(self, bulk_positions, core_positions, method="atomman",
         return disp
     
 
-    def build_cylinder(self, radius,
+    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,
+                       fix_width=10.0, 
+                       self_consistent=None,
                        method="atomman",
-                       verbose=True):
+                       verbose=True, 
+                       return_cyl_mask=False,
+                       return_fix_mask=False,
+                       cyl_mask=None):
         '''
         Build dislocation cylinder for single dislocation system
 
@@ -2865,17 +2875,35 @@ def build_cylinder(self, radius,
         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
+        method: str
+            Displacement solver method name. See self.avail_methods for allowed values
+            method = atomman requires atomman package
+        return_cyl_mask: bool
+            Whether to additionally return the atoms mask
+        return_fix_mask: bool
+            Whether to additionally return the FixAtoms mask
+        cyl_mask: np.array
+            Enforce a specific cyl mask top convert supercell to cyl
+
+        Returns:
+        --------
+        bulk: ase Atoms object
+            Cylinder of bulk, for a formation energy reference
+        disloc: ase Atoms object:
+            Cylinder containing dislocation
+        cyl_mask: np.array
+            If return_cyl_mask=True, return the mask used to convert suqare supercells 
+            to cylinders of correct shape
+        fix_mask: np.array
+            If return_fix_mask=True, return the mask used for the FixAtoms constraint
         '''
 
         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)
+        bulk, disloc, core_positions, cyl_mask, fix_mask = self._build_bulk_cyl(radius, core_positions, fix_width,
+                                                            extension, self_consistent, method, verbose, cyl_mask=cyl_mask)
 
         disloc.info["core_positions"] = [list(core_positions[0, :])]
         disloc.info["burgers_vectors"] = [list(self.burgers)]
@@ -2889,20 +2917,40 @@ def build_cylinder(self, radius,
         # without spoiling the displacement
         # disloc.center(vacuum=2 * fix_width, axis=(0, 1))
 
-        return bulk, disloc
+        out = [bulk, disloc]
+
+        if return_cyl_mask:
+            out.append(cyl_mask)
+        if return_fix_mask:
+            out.append(fix_mask)
+
+        return out
 
     def build_glide_configurations(self, radius,
                                    average_positions=False, **kwargs):
 
+        # Figure out number of dislocation cores
+        if type(self) == CubicCrystalDislocationQuadrupole or issubclass(self.__class__, CubicCrystalDislocationQuadrupole):
+            n_cores = 4
+        elif issubclass(self.__class__, CubicCrystalDissociatedDislocation):
+            n_cores = 2
+        else:
+            n_cores = 1
+
         final_core_position = np.array([self.glide_distance, 0.0, 0.0])
 
-        bulk_ini, disloc_ini = self.build_cylinder(radius,
-                                                   extension=final_core_position,
+        extensions = np.array([final_core_position for i in range(n_cores)])
+
+        bulk_ini, disloc_ini, cyl_mask = self.build_cylinder(radius,
+                                                   core_position=np.zeros(3),
+                                                   extension=extensions,
+                                                   return_cyl_mask=True,
                                                    **kwargs)
 
         _, disloc_fin = self.build_cylinder(radius,
                                             core_position=final_core_position,
-                                            extension=-final_core_position,
+                                            extension=-extensions,
+                                            cyl_mask=cyl_mask,
                                             **kwargs)
         if average_positions:
             # get the fixed atoms constrain
@@ -3363,13 +3411,19 @@ def get_solvers(self, method="atomman"):
 
         return solvers
 
-    def build_cylinder(self, radius, partial_distance=0,
+    def build_cylinder(self, 
+                       radius, 
+                       partial_distance=0,
                        core_position=np.array([0., 0., 0.]),
                        extension=np.array([[0., 0., 0.],
                                           [0., 0., 0.]]),
-                       fix_width=10.0, self_consistent=None,
+                       fix_width=10.0, 
+                       self_consistent=None,
                        method="atomman",
-                       verbose=True):
+                       verbose=True, 
+                       return_cyl_mask=False,
+                       return_fix_mask=False,
+                       cyl_mask=None):
         """
         Overloaded function to make dissociated dislocations.
         Partial distance is provided as an integer to define number
@@ -3383,12 +3437,30 @@ 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.
+        method: str
+            Displacement solver method name. See self.avail_methods for allowed values
+            method = atomman requires atomman package
+        return_cyl_mask: bool
+            Whether to additionally return the atoms mask
+        return_fix_mask: bool
+            Whether to additionally return the FixAtoms mask
+        cyl_mask: np.array
+            Enforce a specific cyl mask top convert supercell to cyl
+
+        Returns:
+        --------
+        bulk: ase Atoms object
+            Cylinder of bulk, for a formation energy reference
+        disloc: ase Atoms object:
+            Cylinder containing dislocation
+        cyl_mask: np.array
+            If return_cyl_mask=True, return the mask used to convert suqare supercells 
+            to cylinders of correct shape
+        fix_mask: np.array
+            If return_fix_mask=True, return the mask used for the FixAtoms constraint
         """
         
         partial_distance_Angstrom = np.array(
@@ -3400,8 +3472,8 @@ def build_cylinder(self, radius, partial_distance=0,
             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)
+        bulk, disloc, core_positions, cyl_mask, fix_mask = self._build_bulk_cyl(radius, core_positions, fix_width, extension,
+                                                            self_consistent, method, verbose, cyl_mask=cyl_mask)
 
         if partial_distance > 0:
             # Specify left & right dislocation separately
@@ -3416,7 +3488,15 @@ def build_cylinder(self, radius, partial_distance=0,
             disloc.info["dislocation_types"] = [self.name]
             disloc.info["dislocation_classes"] = [str(self.__class__)]
 
-        return bulk, disloc
+        
+        out = [bulk, disloc]
+
+        if return_cyl_mask:
+            out.append(cyl_mask)
+        if return_fix_mask:
+            out.append(fix_mask)
+
+        return out
 
 class CubicCrystalDislocationQuadrupole(CubicCrystalDissociatedDislocation):
     burgers_dimensionless = np.zeros(3)
@@ -3831,11 +3911,15 @@ def build_kink_quadrupole(self, z_reps, layer_tol=1e-1, invert_direction=False,
         # 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]
 
+        # TODO: Replace by passing bulk through build_glide_quadrupoles
+        base_bulk = self.build_quadrupole(*args, **kwargs)[0]
+        bulk = base_bulk.copy()
+        
         for image in reps[1:]:
             kink_struct = stack(kink_struct, image)
+            bulk = stack(bulk, base_bulk)
         
         cell = kink_struct.cell[:, :]
 
@@ -3844,6 +3928,9 @@ def build_kink_quadrupole(self, z_reps, layer_tol=1e-1, invert_direction=False,
         kink_struct.set_cell(cell)
         kink_struct.wrap()
 
+        bulk.set_cell(cell)
+        bulk.wrap()
+
         glide_parity = (-direction) % self.glides_per_unit_cell
 
         if glide_parity:
@@ -3851,7 +3938,7 @@ def build_kink_quadrupole(self, z_reps, layer_tol=1e-1, invert_direction=False,
             # glide_parity gives number of layers required to remove
 
             for i in range(glide_parity):
-                atom_heights = kink_struct.get_positions()[:, 2]
+                atom_heights = bulk.get_positions()[:, 2]
                 top_atom = np.max(atom_heights)
 
                 layer_mask = atom_heights >= top_atom - layer_tol
@@ -3859,6 +3946,7 @@ def build_kink_quadrupole(self, z_reps, layer_tol=1e-1, invert_direction=False,
                 avg_layer_pos = np.average(atom_heights[layer_mask])
 
                 kink_struct = kink_struct[~layer_mask]
+                bulk = bulk[~layer_mask]
                 cell[2, 2] = avg_layer_pos
                 kink_struct.set_cell(cell)
                 kink_struct.wrap()

From e73932320ca81f044cb1fbc3ced7e5612dcb8d6e Mon Sep 17 00:00:00 2001
From: thomas-rocke <tom@rocke.co.uk>
Date: Mon, 13 May 2024 14:45:52 +0100
Subject: [PATCH 02/46] ENH: Find smaller representation of cut structure

---
 matscipy/utils.py | 88 +++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 88 insertions(+)

diff --git a/matscipy/utils.py b/matscipy/utils.py
index 287a8d62..d2c6b09c 100644
--- a/matscipy/utils.py
+++ b/matscipy/utils.py
@@ -1,5 +1,6 @@
 import warnings
 import numpy as np
+from ase.utils.structure_comparator import SymmetryEquivalenceCheck
 
 def validate_cubic_cell(a, symbol="w", axes=None, crystalstructure=None, pbc=True):
     '''
@@ -7,6 +8,7 @@ def validate_cubic_cell(a, symbol="w", axes=None, crystalstructure=None, pbc=Tru
 
     For lattice constant + symbol + crystalstructure, build a valid cubic bulk rotated into the frame defined by axes
     For cubic bulk atoms object + crystalstructure, validate structure matches expected cubic bulk geometry, and rotate to frame defined by axes
+    Also, if cubic atoms object is supplied, search for a more condensed representation than is provided by ase.build.cut
 
     a: float or ase.atoms
         EITHER lattice constant (in A), or the cubic bulk structure
@@ -112,6 +114,92 @@ def validate_cubic_cell(a, symbol="w", axes=None, crystalstructure=None, pbc=Tru
 
     return alat, unit_cell
 
+
+def find_condensed_repr(atoms, precision=2):
+    '''
+    Search for a condensed representation of atoms
+    Equivalent to attempting to undo a supercell
+
+    atoms: ase Atoms object
+        structure to condense
+    precision: int
+        Number of decimal places to use in determining whether scaled positions are equal
+
+    returns a condensed copy of atoms, if such a condensed representation is found. Else returns a copy of atoms
+    '''
+    ats = atoms.copy()
+
+    for axis in range(2, -1, -1):
+        ats = find_condensed_repr_along_axis(ats, axis, precision)
+
+    return ats
+
+def find_condensed_repr_along_axis(atoms, axis=-1, precision=2):
+    '''
+    Search for a condensed representation of atoms about axis. 
+    Essentially an inverse to taking a supercell.
+
+    atoms: ase Atoms object
+        structure to condense
+    axis: int
+        axis about which to find a condensed repr
+        axis=-1 essentially will invert a (1, 1, x) supercell
+    precision: int
+        Number of decimal places to use in determining whether scaled positions are equal
+
+    returns a condensed copy of atoms, if such a condensed representation is found. Else returns a copy of atoms
+    '''
+    comp = SymmetryEquivalenceCheck()
+
+    cart_directions = [
+        [1, 2], # axis=0
+        [2, 0], # axis=1
+        [0, 1] # axis=2
+    ]
+
+    dirs = cart_directions[axis]
+
+    ats = atoms.copy()
+    
+    # Find all atoms which are in line with the 0th atom
+    p = np.round(ats.get_scaled_positions(), precision)
+
+    p_diff = (p - p[0, :]) % 1
+
+    matches = np.argwhere((p_diff[:, dirs[0]] < 10**(-precision)) * (p_diff[:, dirs[1]] < 10**(-precision)))[:, 0]
+    min_off = np.inf
+    ret_struct = ats
+
+    for match in matches:
+        if match == 0:
+            # skip i=0
+            continue
+
+        # Fractional change in positions
+        dz = (p[0, axis] - p[match, axis]) % 1.0
+        # Create a test atoms object and cut cell along in axis
+        # Test whether test atoms is equivalent to original structure
+        test_ats = ats.copy()
+        test_cell = test_ats.cell[:, :].copy()
+        test_cell[axis, :] *= dz
+        test_ats.set_cell(test_cell)
+        del_mask = test_ats.get_scaled_positions(wrap=False)[:, axis] > 1.0 - 10**(-precision)
+        test_ats = test_ats[~del_mask]
+
+        # Create a supercell of test ats, and see if it matches the original atoms
+        test_sup = [1] * 3
+        test_sup[axis] = int(1/dz)
+
+        is_equiv = comp.compare(ats, test_ats * tuple(test_sup))
+
+        if is_equiv:
+            if dz < min_off:
+                ret_struct = test_ats.copy()
+                min_off = dz
+
+    return ret_struct
+
+
 def complete_basis(v1, v2=None, normalise=False, nmax=5, tol=1E-6):
     '''
     Generate a complete (v1, v2, v3) orthogonal basis in 3D from v1 and an optional v2

From e279234be2cb428e93c4451a4ec1cde0f157c26b Mon Sep 17 00:00:00 2001
From: thomas-rocke <tom@rocke.co.uk>
Date: Mon, 13 May 2024 14:47:57 +0100
Subject: [PATCH 03/46] ENH: Apply structure representation reduction to
 dislocation unit_cell

---
 matscipy/utils.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/matscipy/utils.py b/matscipy/utils.py
index d2c6b09c..3680f591 100644
--- a/matscipy/utils.py
+++ b/matscipy/utils.py
@@ -62,6 +62,7 @@ def validate_cubic_cell(a, symbol="w", axes=None, crystalstructure=None, pbc=Tru
         alat = a.cell[0, 0]
 
         ats = a.copy()
+        
 
         if crystalstructure is not None:
             # Try to validate that "a" matches expected structure given by crystalstructure
@@ -104,6 +105,10 @@ def validate_cubic_cell(a, symbol="w", axes=None, crystalstructure=None, pbc=Tru
 
         ats = cut(ats, a=axes[0, :], b=axes[1, :], c=axes[2, :])
         rotate(ats, ats.cell[0, :].copy(), [1, 0, 0], ats.cell[1, :].copy(), [0, 1, 0])
+        
+        # cut and rotate can result in a supercell structure. Attempt to find smaller repr 
+        ats = find_condensed_repr(ats)
+        
         unit_cell = ats.copy()
 
         unit_cell.set_pbc(pbc)

From 37d904f779c7f480ee172562ebc41f15e92156c4 Mon Sep 17 00:00:00 2001
From: thomas-rocke <tom@rocke.co.uk>
Date: Tue, 14 May 2024 13:22:41 +0100
Subject: [PATCH 04/46] WIP: Fixes to kink quadrupoles

---
 matscipy/dislocation.py | 63 ++++++++++++++++++++++++++---------------
 1 file changed, 40 insertions(+), 23 deletions(-)

diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index f642cd72..0fc421c3 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -3881,7 +3881,7 @@ def build_glide_quadrupoles(self, nims, invert_direction=False, glide_left=True,
             )[1])
         return images
     
-    def build_kink_quadrupole(self, z_reps, layer_tol=1e-1, invert_direction=False, *args, **kwargs):
+    def build_kink_quadrupole(self, z_reps, layer_decimal_precision=2, 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 
@@ -3911,45 +3911,62 @@ def build_kink_quadrupole(self, z_reps, layer_tol=1e-1, invert_direction=False,
         # 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]
 
         # TODO: Replace by passing bulk through build_glide_quadrupoles
         base_bulk = self.build_quadrupole(*args, **kwargs)[0]
+
+
+        for rep in reps:
+            cell = rep.cell[:, :].copy()
+            cell[2, 0] += self.glide_distance * direction / z_reps
+            rep.set_cell(cell)
+            #rep.wrap()
+
+        base_bulk.set_cell(cell)
+        #base_bulk.wrap()
         bulk = base_bulk.copy()
+        kink_struct = reps[0]
+
+        # Figure out layer spacing
+        p = bulk.get_scaled_positions()[:, 2]
+
+        layer_pos = np.unique(np.round(p, layer_decimal_precision))
+
+        layer_seps = np.diff(layer_pos)[::-1] * bulk.cell[2, 2]
         
+        n_layers = len(layer_seps)
+
         for image in reps[1:]:
             kink_struct = stack(kink_struct, image)
             bulk = stack(bulk, base_bulk)
         
-        cell = kink_struct.cell[:, :]
-
-        cell[2, 0] += self.glide_distance * direction
+        cell = kink_struct.cell[:, :].copy()
+        
+        glide_parity = self.glides_per_unit_cell % 2
+        if glide_parity:
+            # Cell won't be periodic if multiple glides needed to complete periodicity
+            # glide_parity gives number of layers required to remove
 
-        kink_struct.set_cell(cell)
-        kink_struct.wrap()
+            atom_heights = kink_struct.get_positions()[:, 2]
 
-        bulk.set_cell(cell)
-        bulk.wrap()
+            # Get lower bound for atomic layer, to ensure we remove all atoms
+            layer_dist = layer_seps[i%n_layers] + 0.5 * layer_seps[(i+1)%n_layers]
 
-        glide_parity = (-direction) % self.glides_per_unit_cell
+            layer_mask = atom_heights >= cell[2, 2] - layer_dist
 
-        if glide_parity:
-            # Cell won't be periodic if multiple glides needed to complete periodicity
-            # glide_parity gives number of layers required to remove
+            print(i, len(kink_struct), np.sum(layer_mask))
 
-            for i in range(glide_parity):
-                atom_heights = bulk.get_positions()[:, 2]
-                top_atom = np.max(atom_heights)
+            del kink_struct[layer_mask]
+            del bulk[layer_mask]
 
-                layer_mask = atom_heights >= top_atom - layer_tol
+            # Shorten z cell vector to shave off layer_sep[i%n_layer] from cell[2, 2] 
+            cell_frac = layer_seps[i%n_layers] / cell[2, 2]
 
-                avg_layer_pos = np.average(atom_heights[layer_mask])
+            cell[2, :] -= cell[2, :] * cell_frac
+            kink_struct.set_cell(cell)
+            bulk.set_cell(cell)
 
-                kink_struct = kink_struct[~layer_mask]
-                bulk = bulk[~layer_mask]
-                cell[2, 2] = avg_layer_pos
-                kink_struct.set_cell(cell)
-                kink_struct.wrap()
+        kink_struct.wrap()        
         return kink_struct
     
     def view_quad(self, system, *args, **kwargs):

From 062c1d46e42667ff80cce5c35a5158dce119ea78 Mon Sep 17 00:00:00 2001
From: thomas-rocke <tom@rocke.co.uk>
Date: Fri, 17 May 2024 13:08:57 +0100
Subject: [PATCH 05/46] ENH: Coordination number calculation

---
 matscipy/neighbours.py | 15 +++++++++++++++
 1 file changed, 15 insertions(+)

diff --git a/matscipy/neighbours.py b/matscipy/neighbours.py
index b9eaa638..d20f483c 100644
--- a/matscipy/neighbours.py
+++ b/matscipy/neighbours.py
@@ -940,3 +940,18 @@ def find_common_neighbours(i_n, j_n, nat):
         cnl_i1_i2[slice1, 0] = i1
         cnl_i1_i2[slice1, 1] = i_n_2[slice2]
     return cnl_i1_i2, cnl_j1, nl_index_i1_j1, nl_index_i2_j1
+
+
+def coordination(structure, cutoff):
+    '''
+    Get coordination number of each atom in structure
+    based on the number of neighbours within cutoff
+
+    structure: ase Atoms object
+        Structure to take coordination of
+    cutoff: float
+        Radial cutoff used to determine neighbours
+    '''
+
+    i = neighbour_list("i", structure, cutoff)
+    return np.bincount(i)
\ No newline at end of file

From d985ec4e7866b53f6454c633e43d829c3f1345bf Mon Sep 17 00:00:00 2001
From: thomas-rocke <tom@rocke.co.uk>
Date: Sun, 19 May 2024 15:51:07 +0100
Subject: [PATCH 06/46] WIP: Working Quad Kinks

---
 matscipy/dislocation.py | 185 ++++++++++++++++++++++++++++------------
 1 file changed, 131 insertions(+), 54 deletions(-)

diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index 0fc421c3..d0d72bdd 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -41,7 +41,7 @@
 from ase.units import GPa  # unit conversion
 from ase.io import read
 
-from matscipy.neighbours import neighbour_list, mic
+from matscipy.neighbours import neighbour_list, mic, coordination
 from matscipy.elasticity import fit_elastic_constants
 from matscipy.elasticity import Voigt_6x6_to_full_3x3x3x3
 from matscipy.elasticity import cubic_to_Voigt_6x6, coalesce_elastic_constants
@@ -2158,6 +2158,7 @@ def check_duplicates(atoms, distance=0.1):
     # print(f"found {duplicates.sum()} duplicates")
     return duplicates.astype(np.bool)
 
+
 class AnisotropicDislocation:
     """
     Displacement and displacement gradient field of straight dislocation 
@@ -3881,7 +3882,113 @@ def build_glide_quadrupoles(self, nims, invert_direction=False, glide_left=True,
             )[1])
         return images
     
-    def build_kink_quadrupole(self, z_reps, layer_decimal_precision=2, invert_direction=False, *args, **kwargs):
+    def _get_kink_quad_mask(self, ref_bulk, direction, precision=3):
+        '''
+        Get an atom mask and tilted cell which generates a perfect bulk structure
+        Mask is needed in building kinks to ensure continuity of atomic basis
+
+        ref bulk: ase Atoms
+            Bulk structure to use to generate mask
+        direction: int
+            Direction of kink (+1 or -1)
+        precision: int
+            Scaled precision to use when determining a step size
+            Each step is self.alat * 10**_precision Ang
+            
+
+        returns a boolean mask, and a new cell
+
+        '''
+
+        if self.left_dislocation.__class__ == BCCEdge111barDislocation:
+            raise RuntimeError("Kink Quadrupoles do not currently work for BCCEdge111barDislocation")
+
+        bulk = ref_bulk.copy()
+        tilt_bulk = ref_bulk.copy()
+        cell = tilt_bulk.cell[:, :]
+        cell[2, 0] += direction * self.glide_distance
+
+        tilt_bulk.set_cell(cell, scale_atoms=False)
+        tilt_bulk.wrap()
+
+        p = ref_bulk.get_positions()
+
+        step_size = self.alat * 10**(-precision)
+
+        # Shift bulk if atoms too close to the border
+        if any(p[:, 2] < step_size):
+            p[:, 2] += step_size
+            tilt_bulk.set_positions(p)
+            tilt_bulk.wrap()
+            p = tilt_bulk.get_positions()
+
+    
+        # 1st, 2nd, & 3rd nearest neighbour distances
+        # for each crytal structure
+        neigh_dists = {
+        "bcc" : [0.88, 1.01, 1.42],
+        "fcc" : [0.72, 1.01, 1.23],
+        "diamond" : [0.44, 0.72, 0.84]
+        }
+        cutoffs = neigh_dists[self.crystalstructure]
+
+        bulk_coord = [np.min(coordination(bulk, cutoff * self.alat)) for cutoff in cutoffs]
+
+        full_mask = np.ones((len(tilt_bulk)), dtype=bool)
+        atom_heights = np.round(p, precision)[:, 2]
+
+        for i in range(int(cell[2, 2] / step_size)):
+            # Check if tilt_bulk has same 1st, 2nd, & 3rd neighbour coordination
+            # as normal bulk
+            # TODO: This currently does not work for BCCEdge111barDislocation
+            if np.product([
+                all(coordination(tilt_bulk[full_mask], cutoff * self.alat) == bulk_coord[j])
+                for j, cutoff in enumerate(cutoffs)
+                ]):
+                return full_mask, cell
+                
+            # Else, trim the top part of the cell and retry
+            # (atoms outside the cell are also trimmed)
+            cell[2, 2] -= step_size
+            atom_mask = atom_heights > cell[2,2]
+
+            full_mask[atom_mask] = False
+
+            tilt_bulk.set_cell(cell, scale_atoms=False)
+
+        # If we end up here, all layers removed before we found a match. Raise error
+        raise RuntimeError("Could not find a valid periodic kink cell.")
+
+    def _smooth_kink_displacements(self, kink1, kink2, ref_bulk, smoothing_width, base_cell):
+        '''
+        Use a sigmoid
+        '''
+
+        cell = base_cell
+
+        dr1 = mic(kink1.positions - ref_bulk.positions, cell)
+        dr2 = mic(kink2.positions - ref_bulk.positions, cell)
+
+        dr_forward = mic(dr2 - dr1, cell)
+
+
+        kink = ref_bulk.copy()
+
+        x = ref_bulk.positions[:, 2]
+        h = kink.cell[2, 2]
+        sw = smoothing_width
+
+        # ReLU interpolation between kink1 & kink2 displacements 
+        smooth_facs = (x > h - sw) * (x + sw - h) / sw
+
+        smoothed_disp = dr1 + smooth_facs[:, np.newaxis] * dr_forward
+
+        kink.positions += smoothed_disp
+
+        return kink
+
+    def build_kink_quadrupole(self, z_reps=2, layer_decimal_precision=3, invert_direction=False, smooth_width=None,
+                                *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 
@@ -3904,70 +4011,40 @@ def build_kink_quadrupole(self, z_reps, layer_decimal_precision=2, invert_direct
             Fed to self.build_quadrupole() & self.build_glide_quadrupoles
         
         '''
-        assert z_reps > 1
 
         direction = -1 if invert_direction else 1
 
-        # 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)
+        base_bulk, _ = self.build_quadrupole(*args, **kwargs)
 
-        # TODO: Replace by passing bulk through build_glide_quadrupoles
-        base_bulk = self.build_quadrupole(*args, **kwargs)[0]
+        base_quad1, base_quad2 = self.build_glide_quadrupoles(2, glide_left=True, glide_right=True, 
+                                                 invert_direction=invert_direction, *args, **kwargs)
 
-
-        for rep in reps:
-            cell = rep.cell[:, :].copy()
-            cell[2, 0] += self.glide_distance * direction / z_reps
-            rep.set_cell(cell)
-            #rep.wrap()
-
-        base_bulk.set_cell(cell)
-        #base_bulk.wrap()
-        bulk = base_bulk.copy()
-        kink_struct = reps[0]
-
-        # Figure out layer spacing
-        p = bulk.get_scaled_positions()[:, 2]
-
-        layer_pos = np.unique(np.round(p, layer_decimal_precision))
-
-        layer_seps = np.diff(layer_pos)[::-1] * bulk.cell[2, 2]
+        kink_struct1 = base_quad1.copy() * (1, 1, z_reps)
+        kink_struct2 = base_quad2.copy() * (1, 1, z_reps)
         
-        n_layers = len(layer_seps)
-
-        for image in reps[1:]:
-            kink_struct = stack(kink_struct, image)
-            bulk = stack(bulk, base_bulk)
-        
-        cell = kink_struct.cell[:, :].copy()
-        
-        glide_parity = self.glides_per_unit_cell % 2
-        if glide_parity:
-            # Cell won't be periodic if multiple glides needed to complete periodicity
-            # glide_parity gives number of layers required to remove
-
-            atom_heights = kink_struct.get_positions()[:, 2]
-
-            # Get lower bound for atomic layer, to ensure we remove all atoms
-            layer_dist = layer_seps[i%n_layers] + 0.5 * layer_seps[(i+1)%n_layers]
+        bulk = base_bulk.copy() * (1, 1, z_reps)
 
-            layer_mask = atom_heights >= cell[2, 2] - layer_dist
+        bulk.wrap()
+        kink_struct1.wrap()
+        kink_struct2.wrap()
 
-            print(i, len(kink_struct), np.sum(layer_mask))
+        # Cell won't always be periodic, make sure we end up with something that is
+        mask, cell = self._get_kink_quad_mask(bulk, direction, layer_decimal_precision)
+        bulk = bulk[mask]
+        bulk.set_cell(cell, scale_atoms=False)
 
-            del kink_struct[layer_mask]
-            del bulk[layer_mask]
+        kink_struct1 = kink_struct1[mask]
+        kink_struct1.set_cell(cell, scale_atoms=False)
 
-            # Shorten z cell vector to shave off layer_sep[i%n_layer] from cell[2, 2] 
-            cell_frac = layer_seps[i%n_layers] / cell[2, 2]
+        kink_struct2 = kink_struct2[mask]
+        kink_struct2.set_cell(cell, scale_atoms=False)
 
-            cell[2, :] -= cell[2, :] * cell_frac
-            kink_struct.set_cell(cell)
-            bulk.set_cell(cell)
+        if smooth_width is None:
+            smooth_width = 0.5 * self.unit_cell.cell[2, 2]
 
-        kink_struct.wrap()        
-        return kink_struct
+        # Smooth displacements across the periodic boundary
+        kink = self._smooth_kink_displacements(kink_struct1, kink_struct2, bulk, smooth_width, base_bulk.cell[:, :])
+        return kink
     
     def view_quad(self, system, *args, **kwargs):
         '''

From 58582b6f701a4e63ced41632bc410ab9e9e215b2 Mon Sep 17 00:00:00 2001
From: thomas-rocke <tom@rocke.co.uk>
Date: Mon, 20 May 2024 15:08:33 +0100
Subject: [PATCH 07/46] ENH: Add extra fixed_points args when building disloc
 cyls

---
 matscipy/dislocation.py | 73 ++++++++++++++++++++++++-----------------
 1 file changed, 43 insertions(+), 30 deletions(-)

diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index d0d72bdd..f3c3827f 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -2588,7 +2588,7 @@ def _build_supercell(self, targetx, targety):
         return sup
 
     def _build_bulk_cyl(self, radius, core_positions, fix_rad, extension,
-                        self_consistent, method, verbose, cyl_mask=None):
+                        fixed_points, self_consistent, method, verbose, cyl_mask=None, **kwargs):
         '''
         Build bulk cylinder config from args supplied by self.build_cylinder
 
@@ -2643,11 +2643,21 @@ def _build_bulk_cyl(self, radius, core_positions, fix_rad, extension,
         # 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])])
+        # Mask out extensions that are all zeros
+        nonzero_exts = np.sum(np.abs(extension), axis=-1) > 0.0
+
+        exts = exts[nonzero_exts, :]
+
+        fictitious_core_positions = exts
+
+        if fixed_points is not None:
+            fictitious_core_positions = np.vstack([fictitious_core_positions, fixed_points])
+
+        xmax = np.max([np.max(core_positions[:, 0]), np.max(fictitious_core_positions[:, 0])])
+        xmin = np.min([np.min(core_positions[:, 0]), np.min(fictitious_core_positions[:, 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])])
+        ymax = np.max([np.max(core_positions[:, 1]), np.max(fictitious_core_positions[:, 1])])
+        ymin = np.min([np.min(core_positions[:, 1]), np.min(fictitious_core_positions[:, 1])])
 
         xlen = xmax - xmin + 2*radius
         ylen = ymax - ymin + 2*radius
@@ -2669,13 +2679,8 @@ def _build_bulk_cyl(self, radius, core_positions, fix_rad, extension,
 
         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, :]
-
         if cyl_mask is None:
-            mask_positions = np.vstack([new_core_positions, exts + shift])
+            mask_positions = np.vstack([new_core_positions, fictitious_core_positions + shift])
 
             idxs = np.argsort(mask_positions, axis=0)
             mask_positions = np.take_along_axis(mask_positions, idxs, axis=0)
@@ -2691,7 +2696,7 @@ def _build_bulk_cyl(self, radius, core_positions, fix_rad, extension,
 
         disloc.positions += self.displacements(cyl.positions, new_core_positions,
                                                self_consistent=self_consistent, method=method,
-                                               verbose=verbose)
+                                               verbose=verbose, **kwargs)
 
         fix_mask = ~radial_mask_from_polygon2D(cyl.get_positions(), mask_positions, radius - fix_rad, inner=True)
         
@@ -2849,13 +2854,15 @@ def build_cylinder(self,
                        radius,
                        core_position=np.array([0., 0., 0.]),
                        extension=np.array([0, 0, 0]),
+                       fixed_points=None,
                        fix_width=10.0, 
                        self_consistent=None,
                        method="atomman",
                        verbose=True, 
                        return_cyl_mask=False,
                        return_fix_mask=False,
-                       cyl_mask=None):
+                       cyl_mask=None,
+                       **kwargs):
         '''
         Build dislocation cylinder for single dislocation system
 
@@ -2868,6 +2875,9 @@ def build_cylinder(self,
         extension: np.array
             Add extra bulk to the system, to also surround a fictitious core
             that is at position core_position + extension
+        fixed_points: np.array
+            Shape (N, 3) array of points defining extra bulk.
+            Similar to extension arg, but is not relative to core positions 
         fix_width: float
             Defines a region to apply the FixAtoms ase constraint to
             Fixed region is given by (radius - fix_width) <= r <= radius,
@@ -2904,7 +2914,7 @@ def build_cylinder(self,
         ])
 
         bulk, disloc, core_positions, cyl_mask, fix_mask = self._build_bulk_cyl(radius, core_positions, fix_width,
-                                                            extension, self_consistent, method, verbose, cyl_mask=cyl_mask)
+                                                            extension, fixed_points, self_consistent, method, verbose, cyl_mask=cyl_mask, **kwargs)
 
         disloc.info["core_positions"] = [list(core_positions[0, :])]
         disloc.info["burgers_vectors"] = [list(self.burgers)]
@@ -2927,30 +2937,23 @@ def build_cylinder(self,
 
         return out
 
-    def build_glide_configurations(self, radius,
-                                   average_positions=False, **kwargs):
+    def build_glide_configurations(self, radius, average_positions=False, **kwargs):
 
-        # Figure out number of dislocation cores
-        if type(self) == CubicCrystalDislocationQuadrupole or issubclass(self.__class__, CubicCrystalDislocationQuadrupole):
-            n_cores = 4
-        elif issubclass(self.__class__, CubicCrystalDissociatedDislocation):
-            n_cores = 2
-        else:
-            n_cores = 1
+        initial_core_position = np.array([0, 0, 0])
 
-        final_core_position = np.array([self.glide_distance, 0.0, 0.0])
+        final_core_position = np.array([self.glide_distance, 0, 0])
 
-        extensions = np.array([final_core_position for i in range(n_cores)])
+        fixed_points = np.vstack([initial_core_position, final_core_position])
 
         bulk_ini, disloc_ini, cyl_mask = self.build_cylinder(radius,
-                                                   core_position=np.zeros(3),
-                                                   extension=extensions,
+                                                   core_position=initial_core_position,
+                                                   fixed_points=fixed_points,
                                                    return_cyl_mask=True,
                                                    **kwargs)
 
         _, disloc_fin = self.build_cylinder(radius,
                                             core_position=final_core_position,
-                                            extension=-extensions,
+                                            fixed_points=fixed_points,
                                             cyl_mask=cyl_mask,
                                             **kwargs)
         if average_positions:
@@ -3079,6 +3082,11 @@ def build_impurity_cylinder(self, disloc, impurity, radius,
 
         return impurities_disloc
     
+    def build_kink_cyl(self, kink_map=[0, 1], *args, **kwargs) -> None:
+        map = np.array(kink_map, dtype=int)
+        map -= np.min(map)
+        range = np.max(map)
+    
     @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):
@@ -3418,13 +3426,15 @@ def build_cylinder(self,
                        core_position=np.array([0., 0., 0.]),
                        extension=np.array([[0., 0., 0.],
                                           [0., 0., 0.]]),
+                       fixed_points=None,
                        fix_width=10.0, 
                        self_consistent=None,
                        method="atomman",
                        verbose=True, 
                        return_cyl_mask=False,
                        return_fix_mask=False,
-                       cyl_mask=None):
+                       cyl_mask=None,
+                       **kwargs):
         """
         Overloaded function to make dissociated dislocations.
         Partial distance is provided as an integer to define number
@@ -3441,6 +3451,9 @@ def build_cylinder(self,
         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.
+        fixed_points: np.array
+            Shape (N, 3) array of points defining extra bulk.
+            Similar to extension arg, but is not relative to core positions 
         method: str
             Displacement solver method name. See self.avail_methods for allowed values
             method = atomman requires atomman package
@@ -3474,7 +3487,7 @@ def build_cylinder(self,
         ])
 
         bulk, disloc, core_positions, cyl_mask, fix_mask = self._build_bulk_cyl(radius, core_positions, fix_width, extension,
-                                                            self_consistent, method, verbose, cyl_mask=cyl_mask)
+                                                            fixed_points, self_consistent, method, verbose, cyl_mask=cyl_mask, **kwargs)
 
         if partial_distance > 0:
             # Specify left & right dislocation separately

From 4b1591fa63e47cf5185b0369ac8ad2cd3a507d55 Mon Sep 17 00:00:00 2001
From: thomas-rocke <tom@rocke.co.uk>
Date: Mon, 20 May 2024 15:27:40 +0100
Subject: [PATCH 08/46] ENH: Disloc Kink Cyl generation

---
 matscipy/dislocation.py | 24 +++++++++++++++++++++++-
 1 file changed, 23 insertions(+), 1 deletion(-)

diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index f3c3827f..6b08e03a 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -3082,10 +3082,32 @@ def build_impurity_cylinder(self, disloc, impurity, radius,
 
         return impurities_disloc
     
-    def build_kink_cyl(self, kink_map=[0, 1], *args, **kwargs) -> None:
+    def build_kink_cyl(self, kink_map=[0, 1], *args, **kwargs):
         map = np.array(kink_map, dtype=int)
         map -= np.min(map)
         range = np.max(map)
+
+        unique_kinks = np.sort(np.unique(map))
+
+        glide_structs = {}
+
+        fixed_points = np.array([
+            [0, 0, 0],
+            [self.glide_distance * range, 0, 0]
+        ])
+
+        for kink_pos in unique_kinks:
+            glide_structs[kink_pos] = self.build_cylinder(*args, 
+                                                          fixed_points=fixed_points, 
+                                                          core_position=np.array([kink_pos * self.glide_distance, 0 , 0]),
+                                                          **kwargs)[1]
+            
+        kink_cyl = glide_structs[map[0]].copy()
+
+        for i in map[1:]:
+            kink_cyl = stack(kink_cyl, glide_structs[i])
+
+        return kink_cyl
     
     @staticmethod
     def view_cyl(system, scale=0.5, CNA_color=True, add_bonds=False,

From 634a19c583d977da925b3819b6c3262d5ad6522a Mon Sep 17 00:00:00 2001
From: thomas-rocke <tom@rocke.co.uk>
Date: Mon, 20 May 2024 16:46:48 +0100
Subject: [PATCH 09/46] MAINT: Skip single kink quad test

---
 tests/test_dislocation.py | 9 +++------
 1 file changed, 3 insertions(+), 6 deletions(-)

diff --git a/tests/test_dislocation.py b/tests/test_dislocation.py
index 0f4928c2..251a7268 100644
--- a/tests/test_dislocation.py
+++ b/tests/test_dislocation.py
@@ -904,6 +904,7 @@ def test_glide_configs(self, disloc, subtests):
             self.assertAtomsAlmostEqual(configs[0], ini_quad)
             self.assertAtomsAlmostEqual(configs[1], fin_quad)
 
+    @pytest.skip()
     def test_single_kink_quadrupole(self, disloc):
         '''
         Validate that generation of single kink structures runs without errors
@@ -912,14 +913,10 @@ def test_single_kink_quadrupole(self, 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)
+        kink_struct = d.build_kink_quadrupole(z_reps=4, glide_separation=4)
 
         # 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)
+        self.assertNotAlmostEqual(kink_struct.cell[2, 0], 0.0)
 
 
 

From a866e14a731abce3502f9bddb4b9e59946f7b226 Mon Sep 17 00:00:00 2001
From: thomas-rocke <tom@rocke.co.uk>
Date: Mon, 20 May 2024 16:51:26 +0100
Subject: [PATCH 10/46] MAINT: Add pytest timeout of 300s through
 pytest-timeout

---
 .github/workflows/tests.yml | 2 +-
 pyproject.toml              | 1 +
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index d78f64b5..9bfa494e 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -41,7 +41,7 @@ jobs:
       - name: pytest
         run: |
           cd tests
-          pytest -v --junitxml=report.xml --durations=20
+          pytest -v --junitxml=report.xml --durations=20 --timeout=300
 
       - name: report
         uses: mikepenz/action-junit-report@v2
diff --git a/pyproject.toml b/pyproject.toml
index 6cf0ac57..df7c32b8 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -28,6 +28,7 @@ dependencies = [
 test = [
     "pytest",
     "pytest-subtests",
+    "pytest-timeout",
     "sympy",
     "atomman",
     "ovito"

From 87e21ed0d680084682e0e0634dea45d4b79c602f Mon Sep 17 00:00:00 2001
From: thomas-rocke <tom@rocke.co.uk>
Date: Mon, 20 May 2024 16:52:21 +0100
Subject: [PATCH 11/46] BUG: pytest.skip -> pytest.mark.skip

---
 tests/test_dislocation.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_dislocation.py b/tests/test_dislocation.py
index 251a7268..dad5ba9e 100644
--- a/tests/test_dislocation.py
+++ b/tests/test_dislocation.py
@@ -904,7 +904,7 @@ def test_glide_configs(self, disloc, subtests):
             self.assertAtomsAlmostEqual(configs[0], ini_quad)
             self.assertAtomsAlmostEqual(configs[1], fin_quad)
 
-    @pytest.skip()
+    @pytest.mark.skip()
     def test_single_kink_quadrupole(self, disloc):
         '''
         Validate that generation of single kink structures runs without errors

From a62c43ec847d65c039f2d5509a3a1c3be00c4761 Mon Sep 17 00:00:00 2001
From: thomas-rocke <tom@rocke.co.uk>
Date: Tue, 21 May 2024 00:46:26 +0100
Subject: [PATCH 12/46] BUG: Fix fixed_points not scaling off
 unit_cell_core_position

---
 matscipy/dislocation.py | 91 ++++++++++++++++++++++-------------------
 1 file changed, 49 insertions(+), 42 deletions(-)

diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index 6b08e03a..3cd74280 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -2585,6 +2585,7 @@ def _build_supercell(self, targetx, targety):
         yreps = np.ceil(targety / (base_cell.cell[1, 1] - 1e-2)).astype(int)
 
         sup = base_cell.copy() * (xreps, yreps, 1)
+        sup.wrap()
         return sup
 
     def _build_bulk_cyl(self, radius, core_positions, fix_rad, extension,
@@ -2619,45 +2620,53 @@ def _build_bulk_cyl(self, radius, core_positions, fix_rad, extension,
         '''
         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
-
-        # Mask out extensions that are all zeros
-        nonzero_exts = np.sum(np.abs(extension), axis=-1) > 0.0
-
-        exts = exts[nonzero_exts, :]
+        # Validate core position, extension, and fixed_points args
+        errs = []
+        for argname, arg in [["core_position", core_positions], ["extension", extension], ["fixed_points", fixed_points]]:
+            if type(arg) == np.ndarray:
+                if arg.shape[-1] != 3:
+                    # Needs to be 3d vectors
+                    errs.append(f"Argument {argname} misspecified. Should be an array of shape (N, 3)")
+                if len(arg.shape) == 1:
+                    # Convert all arrays to 2D
+                    arg.reshape((np.newaxis, arg.shape[-1]))
+            elif arg is not None:
+                # non-array, and not None arg provided, which is not allowed
+                errs.append(f"Argument {argname} misspecified. Should be an array of shape (N, 3)")
+        if len(errs):
+            raise RuntimeError("\n".join(errs))
+
+        fictitious_core_positions = core_positions.copy()
+
+        if extension is not None:
+            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=}")
 
-        fictitious_core_positions = exts
+            # Only care about non-zero extensions
+            fictitious_core_positions = np.vstack([fictitious_core_positions, extension + core_positions])
 
         if fixed_points is not None:
-            fictitious_core_positions = np.vstack([fictitious_core_positions, fixed_points])
+            fictitious_core_positions = np.vstack([fictitious_core_positions, fixed_points + self.unit_cell_core_position])
+        
+        fictitious_core_positions = np.unique(fictitious_core_positions, axis=0)
 
-        xmax = np.max([np.max(core_positions[:, 0]), np.max(fictitious_core_positions[:, 0])])
-        xmin = np.min([np.min(core_positions[:, 0]), np.min(fictitious_core_positions[:, 0])])
+        # Get bounds for supercell box
+        xmax = np.max([np.max(fictitious_core_positions[:, 0])])
+        xmin = np.min([np.min(fictitious_core_positions[:, 0])])
 
-        ymax = np.max([np.max(core_positions[:, 1]), np.max(fictitious_core_positions[:, 1])])
-        ymin = np.min([np.min(core_positions[:, 1]), np.min(fictitious_core_positions[:, 1])])
+        ymax = np.max([np.max(fictitious_core_positions[:, 1])])
+        ymin = np.min([np.min(fictitious_core_positions[:, 1])])
 
         xlen = xmax - xmin + 2*radius
         ylen = ymax - ymin + 2*radius
@@ -2667,7 +2676,6 @@ def _build_bulk_cyl(self, radius, core_positions, fix_rad, extension,
         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()
@@ -2679,13 +2687,12 @@ def _build_bulk_cyl(self, radius, core_positions, fix_rad, extension,
 
         new_core_positions += shift
 
-        if cyl_mask is None:
-            mask_positions = np.vstack([new_core_positions, fictitious_core_positions + shift])
+        # Sort positions from smallest to largest
+        idxs = np.argsort(fictitious_core_positions, axis=0)
+        fictitious_core_positions = np.take_along_axis(fictitious_core_positions, idxs, axis=0) + 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)
+        if cyl_mask is None:    
+            mask = radial_mask_from_polygon2D(sup.get_positions(), fictitious_core_positions, radius, inner=True)
         else:
             mask = cyl_mask
 
@@ -2698,7 +2705,8 @@ def _build_bulk_cyl(self, radius, core_positions, fix_rad, extension,
                                                self_consistent=self_consistent, method=method,
                                                verbose=verbose, **kwargs)
 
-        fix_mask = ~radial_mask_from_polygon2D(cyl.get_positions(), mask_positions, radius - fix_rad, inner=True)
+        # Generate FixAtoms constraint
+        fix_mask = ~radial_mask_from_polygon2D(cyl.get_positions(), fictitious_core_positions, radius - fix_rad, inner=True)
         
         if fix_rad:
             fix_atoms = FixAtoms(mask=fix_mask)
@@ -2853,7 +2861,7 @@ def displacements(self, bulk_positions, core_positions, method="atomman",
     def build_cylinder(self, 
                        radius,
                        core_position=np.array([0., 0., 0.]),
-                       extension=np.array([0, 0, 0]),
+                       extension=None,
                        fixed_points=None,
                        fix_width=10.0, 
                        self_consistent=None,
@@ -2986,7 +2994,7 @@ def build_glide_configurations(self, radius, average_positions=False, **kwargs):
     def build_impurity_cylinder(self, disloc, impurity, radius,
                                 imp_symbol="H",
                                 core_position=np.array([0., 0., 0.]),
-                                extension=np.array([0., 0., 0.]),
+                                extension=None,
                                 self_consistent=False,
                                 extra_bulk_at_core=False,
                                 core_radius=0.5,
@@ -3101,7 +3109,6 @@ def build_kink_cyl(self, kink_map=[0, 1], *args, **kwargs):
                                                           fixed_points=fixed_points, 
                                                           core_position=np.array([kink_pos * self.glide_distance, 0 , 0]),
                                                           **kwargs)[1]
-            
         kink_cyl = glide_structs[map[0]].copy()
 
         for i in map[1:]:

From f17a9d90c7a28c336eaf469cc4ca3987fa8ce637 Mon Sep 17 00:00:00 2001
From: thomas-rocke <tom@rocke.co.uk>
Date: Thu, 23 May 2024 16:12:18 +0100
Subject: [PATCH 13/46] MAINT: Promote _smooth_kink_displacements to CCD from
 CCDQ

---
 matscipy/dislocation.py | 58 +++++++++++++++++++++--------------------
 1 file changed, 30 insertions(+), 28 deletions(-)

diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index 3cd74280..40744f08 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -3089,6 +3089,36 @@ def build_impurity_cylinder(self, disloc, impurity, radius,
         impurities_disloc.cell = disloc.cell
 
         return impurities_disloc
+
+    def _smooth_kink_displacements(self, kink1, kink2, ref_bulk, smoothing_width, base_cell):
+        '''
+        Use ReLU interpolation to smooth the displacment fields across a kink boundary between kink1 and kink2.
+        Returns a copy of kink1, but with the displacement field at ref_bulk.positions[:, 2] > base_cell[2, 2] - smoothing width
+        linearly interpolated. 
+        '''
+
+        cell = base_cell
+
+        dr1 = mic(kink1.positions - ref_bulk.positions, cell)
+        dr2 = mic(kink2.positions - ref_bulk.positions, cell)
+
+        dr_forward = mic(dr2 - dr1, cell)
+
+
+        kink = ref_bulk.copy()
+
+        x = ref_bulk.positions[:, 2]
+        h = kink.cell[2, 2]
+        sw = smoothing_width
+
+        # ReLU interpolation between kink1 & kink2 displacements 
+        smooth_facs = (x > h - sw) * (x + sw - h) / sw
+
+        smoothed_disp = dr1 + smooth_facs[:, np.newaxis] * dr_forward
+
+        kink.positions += smoothed_disp
+
+        return kink
     
     def build_kink_cyl(self, kink_map=[0, 1], *args, **kwargs):
         map = np.array(kink_map, dtype=int)
@@ -4001,34 +4031,6 @@ def _get_kink_quad_mask(self, ref_bulk, direction, precision=3):
         # If we end up here, all layers removed before we found a match. Raise error
         raise RuntimeError("Could not find a valid periodic kink cell.")
 
-    def _smooth_kink_displacements(self, kink1, kink2, ref_bulk, smoothing_width, base_cell):
-        '''
-        Use a sigmoid
-        '''
-
-        cell = base_cell
-
-        dr1 = mic(kink1.positions - ref_bulk.positions, cell)
-        dr2 = mic(kink2.positions - ref_bulk.positions, cell)
-
-        dr_forward = mic(dr2 - dr1, cell)
-
-
-        kink = ref_bulk.copy()
-
-        x = ref_bulk.positions[:, 2]
-        h = kink.cell[2, 2]
-        sw = smoothing_width
-
-        # ReLU interpolation between kink1 & kink2 displacements 
-        smooth_facs = (x > h - sw) * (x + sw - h) / sw
-
-        smoothed_disp = dr1 + smooth_facs[:, np.newaxis] * dr_forward
-
-        kink.positions += smoothed_disp
-
-        return kink
-
     def build_kink_quadrupole(self, z_reps=2, layer_decimal_precision=3, invert_direction=False, smooth_width=None,
                                 *args, **kwargs):
         '''

From 863fea85ef82971889571b5629a2987a7e0ab5fb Mon Sep 17 00:00:00 2001
From: thomas-rocke <tom@rocke.co.uk>
Date: Thu, 23 May 2024 16:54:55 +0100
Subject: [PATCH 14/46] ENH: Displacmeent smoothing for cyl kinks

---
 matscipy/dislocation.py | 37 ++++++++++++++++++++++++++-----------
 1 file changed, 26 insertions(+), 11 deletions(-)

diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index 40744f08..0873771d 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -3120,29 +3120,44 @@ def _smooth_kink_displacements(self, kink1, kink2, ref_bulk, smoothing_width, ba
 
         return kink
     
-    def build_kink_cyl(self, kink_map=[0, 1], *args, **kwargs):
-        map = np.array(kink_map, dtype=int)
-        map -= np.min(map)
-        range = np.max(map)
+    def build_kink_cyl(self, kink_map=[0, 1], smooth_width=None, *args, **kwargs):
+        kmap = np.array(kink_map, dtype=int)
+        kmap -= np.min(kmap)
+        krange = np.max(kmap)
 
-        unique_kinks = np.sort(np.unique(map))
+        unique_kinks = np.sort(np.unique(kmap))
 
         glide_structs = {}
 
         fixed_points = np.array([
             [0, 0, 0],
-            [self.glide_distance * range, 0, 0]
+            [self.glide_distance * krange, 0, 0]
         ])
 
         for kink_pos in unique_kinks:
-            glide_structs[kink_pos] = self.build_cylinder(*args, 
+            ref_bulk, glide_structs[kink_pos] = self.build_cylinder(*args, 
                                                           fixed_points=fixed_points, 
                                                           core_position=np.array([kink_pos * self.glide_distance, 0 , 0]),
-                                                          **kwargs)[1]
-        kink_cyl = glide_structs[map[0]].copy()
+                                                          **kwargs)
+        
+        kink_structs = [glide_structs[midx].copy() for midx in kmap]
+
+        # Smooth displacements across kink boundaries
+        if smooth_width is None:
+            smooth_width = 0.5 * self.unit_cell.cell[2, 2]
+
+        N = len(kink_structs)
+
+        smoothed_kink_structs = []
+
+        for i in range(N):
+            smoothed_kink_structs.append(
+                self._smooth_kink_displacements(kink_structs[i], kink_structs[(i+1) % N], ref_bulk, smooth_width, ref_bulk.cell[:, :])
+            )
 
-        for i in map[1:]:
-            kink_cyl = stack(kink_cyl, glide_structs[i])
+        kink_cyl = smoothed_kink_structs[0].copy()
+        for i in range(1, N):
+            kink_cyl = stack(kink_cyl, smoothed_kink_structs[i])
 
         return kink_cyl
     

From d06a5e7f6b94a8c42022908ab586e4df4406589e Mon Sep 17 00:00:00 2001
From: thomas-rocke <tom@rocke.co.uk>
Date: Thu, 23 May 2024 17:14:48 +0100
Subject: [PATCH 15/46] ENH: Multi-kink quads

---
 matscipy/dislocation.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index 0873771d..85fc8b79 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -4046,7 +4046,7 @@ def _get_kink_quad_mask(self, ref_bulk, direction, precision=3):
         # If we end up here, all layers removed before we found a match. Raise error
         raise RuntimeError("Could not find a valid periodic kink cell.")
 
-    def build_kink_quadrupole(self, z_reps=2, layer_decimal_precision=3, invert_direction=False, smooth_width=None,
+    def build_kink_quadrupole(self, z_reps=2, layer_decimal_precision=3, n_kink=1, invert_direction=False, smooth_width=None,
                                 *args, **kwargs):
         '''
         Construct a quadrupole structure providing an initial guess of the dislocation kink
@@ -4088,7 +4088,7 @@ def build_kink_quadrupole(self, z_reps=2, layer_decimal_precision=3, invert_dire
         kink_struct2.wrap()
 
         # Cell won't always be periodic, make sure we end up with something that is
-        mask, cell = self._get_kink_quad_mask(bulk, direction, layer_decimal_precision)
+        mask, cell = self._get_kink_quad_mask(bulk, direction * n_kink, layer_decimal_precision)
         bulk = bulk[mask]
         bulk.set_cell(cell, scale_atoms=False)
 

From a1d76ec8b9c681ffd507b20b7acbd3e84d59a9a0 Mon Sep 17 00:00:00 2001
From: thomas-rocke <tom@rocke.co.uk>
Date: Tue, 28 May 2024 11:47:12 +0100
Subject: [PATCH 16/46] ENH: Kink networks in quadrupoles

---
 matscipy/dislocation.py | 72 ++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 71 insertions(+), 1 deletion(-)

diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index 85fc8b79..53a245c6 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -3988,7 +3988,8 @@ def _get_kink_quad_mask(self, ref_bulk, direction, precision=3):
         '''
 
         if self.left_dislocation.__class__ == BCCEdge111barDislocation:
-            raise RuntimeError("Kink Quadrupoles do not currently work for BCCEdge111barDislocation")
+            #raise RuntimeError("Kink Quadrupoles do not currently work for BCCEdge111barDislocation")
+            pass
 
         bulk = ref_bulk.copy()
         tilt_bulk = ref_bulk.copy()
@@ -4104,6 +4105,75 @@ def build_kink_quadrupole(self, z_reps=2, layer_decimal_precision=3, n_kink=1, i
         # Smooth displacements across the periodic boundary
         kink = self._smooth_kink_displacements(kink_struct1, kink_struct2, bulk, smooth_width, base_bulk.cell[:, :])
         return kink
+
+    def build_kink_quadrupole_network(self, kink_map=[0, 1], layer_decimal_precision=3, smooth_width=None,
+                                *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
+        
+        '''
+
+        kmap = np.array(kink_map, dtype=int)
+        kmap -= np.min(kmap)
+        krange = np.max(kmap)
+
+        unique_kinks = np.sort(np.unique(kmap))
+
+        glide_structs = {}
+
+        for kink_pos in unique_kinks:
+            off = np.array([kink_pos * self.glide_distance, 0 , 0])
+            ref_bulk, glide_structs[kink_pos] = self.build_quadrupole(*args, 
+                                                          left_offset=off, right_offset=off,
+                                                          **kwargs)
+        
+        kink_structs = [glide_structs[midx].copy() for midx in kmap]
+
+        # Smooth displacements across kink boundaries
+        if smooth_width is None:
+            smooth_width = 0.5 * self.unit_cell.cell[2, 2]
+
+        N = len(kink_structs)
+
+        smoothed_kink_structs = []
+        # Extra copy of last struct, so last one gets correct unsmoothed displacements
+        kink_structs.append(kink_structs[-1].copy())
+        for i in range(N):
+            smoothed_kink_structs.append(
+                self._smooth_kink_displacements(kink_structs[i], kink_structs[i+1], ref_bulk, smooth_width, ref_bulk.cell[:, :])
+            )
+
+        direction = kmap[-1]
+
+        bulk = ref_bulk.copy() * (1, 1, N)
+        kink = smoothed_kink_structs[0].copy()
+        for i in range(1, N):
+            kink = stack(kink, smoothed_kink_structs[i])
+
+        # Cell won't always be periodic, make sure we end up with something that is
+        mask, cell = self._get_kink_quad_mask(bulk, direction, layer_decimal_precision)
+
+        kink = kink[mask]
+        kink.set_cell(cell, scale_atoms=False)
+        return kink
     
     def view_quad(self, system, *args, **kwargs):
         '''

From 15220f2f07258e6d64acd8744c1a1047f205943c Mon Sep 17 00:00:00 2001
From: thomas-rocke <tom@rocke.co.uk>
Date: Tue, 28 May 2024 12:15:03 +0100
Subject: [PATCH 17/46] WIP: Start of disloc mobility docs

---
 docs/applications/disloc_mobility.ipynb | 55 +++++++++++++++++++++++++
 1 file changed, 55 insertions(+)
 create mode 100644 docs/applications/disloc_mobility.ipynb

diff --git a/docs/applications/disloc_mobility.ipynb b/docs/applications/disloc_mobility.ipynb
new file mode 100644
index 00000000..84fba5ad
--- /dev/null
+++ b/docs/applications/disloc_mobility.ipynb
@@ -0,0 +1,55 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Modelling the Mobility of Dislocations with `matscipy`\n",
+    "\n",
+    "The dislocation modelling toolkit provided in `matscipy.dislocation` includes functions to assist in characterising the mobility of dislocations in various systems. The two supported mechanisms are glide and kink, which both allow the dislocation to move in the glide direction.\n",
+    "\n",
+    "## Dislocation Glide\n",
+    "Dislocation Glide is where the entire dislocation line migrates across a glide barrier to advance the dislocation. These kinds of structures are easy to generate both in cylindrical and in quadrupolar structures.\n",
+    "\n",
+    "\n",
+    "### Glide in Dislocation Cylinders\n",
+    "As an example of modelling dislocation glide in cylindrical cells, let's look at the BCC screw dislocation:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from matscipy.dislocation import get_elastic_constants, BCCScrew111Dislocation\n",
+    "# the calculator to provide forces and energies from the potential\n",
+    "from matscipy.calculators.eam import EAM\n",
+    "eam_calc = EAM(\"../../tests/w_eam4.fs\")\n",
+    "\n",
+    "# the function accepts any ASE type of calculator\n",
+    "alat, C11, C12, C44 = get_elastic_constants(calculator=eam_calc, symbol=\"W\", verbose=False)\n",
+    "\n",
+    "W_screw = BCCScrew111Dislocation(alat, C11, C12, C44, symbol=\"W\")\n",
+    "\n",
+    "bulk, screw = W_screw.build_dislocation_cylinder(radius=30)\n",
+    "W_screw_bulk, W_screw_dislo = W_screw.build_cylinder(radius=20)\n",
+    "W_screw.view_cyl(W_screw_dislo)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "language_info": {
+   "name": "python"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

From a4bde2dbd98647c38e2f0ea31cb655de6be58253 Mon Sep 17 00:00:00 2001
From: thomas-rocke <tom@rocke.co.uk>
Date: Thu, 30 May 2024 16:23:57 +0100
Subject: [PATCH 18/46] BUG: Fix issue with SymmetryEquivalenceCheck modifying
 ats

---
 matscipy/utils.py | 22 ++++++++++++++--------
 1 file changed, 14 insertions(+), 8 deletions(-)

diff --git a/matscipy/utils.py b/matscipy/utils.py
index 3680f591..ff2175a8 100644
--- a/matscipy/utils.py
+++ b/matscipy/utils.py
@@ -133,7 +133,6 @@ def find_condensed_repr(atoms, precision=2):
     returns a condensed copy of atoms, if such a condensed representation is found. Else returns a copy of atoms
     '''
     ats = atoms.copy()
-
     for axis in range(2, -1, -1):
         ats = find_condensed_repr_along_axis(ats, axis, precision)
 
@@ -165,37 +164,44 @@ def find_condensed_repr_along_axis(atoms, axis=-1, precision=2):
     dirs = cart_directions[axis]
 
     ats = atoms.copy()
+    ats.wrap()
+
+    # Choose an origin atom, closest to cell origin
+    origin_idx = np.argmin(np.linalg.norm(ats.positions, axis=-1))
     
-    # Find all atoms which are in line with the 0th atom
+    # Find all atoms which are in line with the origin_idx th atom
     p = np.round(ats.get_scaled_positions(), precision)
 
-    p_diff = (p - p[0, :]) % 1
+    p_diff = (p - p[origin_idx, :]) % 1
 
     matches = np.argwhere((p_diff[:, dirs[0]] < 10**(-precision)) * (p_diff[:, dirs[1]] < 10**(-precision)))[:, 0]
     min_off = np.inf
     ret_struct = ats
 
     for match in matches:
-        if match == 0:
-            # skip i=0
+        if match == origin_idx:
+            # skip i=origin_idx
             continue
 
         # Fractional change in positions
-        dz = (p[0, axis] - p[match, axis]) % 1.0
+        dz = (p[origin_idx, axis] - p[match, axis]) % 1.0
         # Create a test atoms object and cut cell along in axis
         # Test whether test atoms is equivalent to original structure
         test_ats = ats.copy()
         test_cell = test_ats.cell[:, :].copy()
         test_cell[axis, :] *= dz
         test_ats.set_cell(test_cell)
-        del_mask = test_ats.get_scaled_positions(wrap=False)[:, axis] > 1.0 - 10**(-precision)
+
+        sp = np.round(test_ats.get_scaled_positions(wrap=False), precision)
+
+        del_mask = sp[:, axis] > 1.0 - 10**(-precision)
         test_ats = test_ats[~del_mask]
 
         # Create a supercell of test ats, and see if it matches the original atoms
         test_sup = [1] * 3
         test_sup[axis] = int(1/dz)
 
-        is_equiv = comp.compare(ats, test_ats * tuple(test_sup))
+        is_equiv = comp.compare(ats.copy(), test_ats * tuple(test_sup))
 
         if is_equiv:
             if dz < min_off:

From 09846da3e949705205e9b4bdeb133cd1c049db5d Mon Sep 17 00:00:00 2001
From: Thomas Rocke <45517493+thomas-rocke@users.noreply.github.com>
Date: Thu, 13 Jun 2024 11:46:21 +0100
Subject: [PATCH 19/46] MAINT: Convert to using np.atleast_2d

---
 matscipy/dislocation.py | 33 +++++++++++++++++++++++++--------
 1 file changed, 25 insertions(+), 8 deletions(-)

diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index 53a245c6..1c4249bb 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -2618,24 +2618,37 @@ def _build_bulk_cyl(self, radius, core_positions, fix_rad, extension,
             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
 
-        # Validate core position, extension, and fixed_points args
-        errs = []
-        for argname, arg in [["core_position", core_positions], ["extension", extension], ["fixed_points", fixed_points]]:
+        def enforce_arr_shape(argname, arg, errs):
+            '''
+            Enforce correct array dimensionality for 
+            
+            '''
             if type(arg) == np.ndarray:
                 if arg.shape[-1] != 3:
                     # Needs to be 3d vectors
                     errs.append(f"Argument {argname} misspecified. Should be an array of shape (N, 3)")
                 if len(arg.shape) == 1:
                     # Convert all arrays to 2D
-                    arg.reshape((np.newaxis, arg.shape[-1]))
+                    arg = np.atleast_2d(arg)
             elif arg is not None:
                 # non-array, and not None arg provided, which is not allowed
                 errs.append(f"Argument {argname} misspecified. Should be an array of shape (N, 3)")
+            return arg, errs
+
+
+        from matscipy.utils import radial_mask_from_polygon2D
+
+        if self_consistent is None:
+            self_consistent = self.self_consistent
+
+        # Validate core position, extension, and fixed_points args
+        errs = []
+        
+        core_positions, errs = enforce_arr_shape("core_positions", core_positions, errs)
+        extension, errs = enforce_arr_shape("extension", extension, errs)
+        fixed_points, errs = enforce_arr_shape("fixed_points", fixed_points, errs)
+
         if len(errs):
             raise RuntimeError("\n".join(errs))
 
@@ -2953,6 +2966,10 @@ def build_glide_configurations(self, radius, average_positions=False, **kwargs):
 
         fixed_points = np.vstack([initial_core_position, final_core_position])
 
+        if "partial_distance" in kwargs.keys():
+            # Account for dissociated glide
+            fixed_points[1, 0] += kwargs["partial_distance"] * self.glide_distance
+
         bulk_ini, disloc_ini, cyl_mask = self.build_cylinder(radius,
                                                    core_position=initial_core_position,
                                                    fixed_points=fixed_points,

From 1146543c77ef5f853ee9c6c134cb29a46a6de30a Mon Sep 17 00:00:00 2001
From: Thomas Rocke <45517493+thomas-rocke@users.noreply.github.com>
Date: Wed, 31 Jul 2024 15:38:01 +0100
Subject: [PATCH 20/46] WIP: Disloc mobility docs

---
 docs/applications/disloc_mobility.ipynb | 238 +++++++++++++++++++++++-
 1 file changed, 228 insertions(+), 10 deletions(-)

diff --git a/docs/applications/disloc_mobility.ipynb b/docs/applications/disloc_mobility.ipynb
index 84fba5ad..214c59ed 100644
--- a/docs/applications/disloc_mobility.ipynb
+++ b/docs/applications/disloc_mobility.ipynb
@@ -13,7 +13,7 @@
     "\n",
     "\n",
     "### Glide in Dislocation Cylinders\n",
-    "As an example of modelling dislocation glide in cylindrical cells, let's look at the BCC screw dislocation:"
+    "As an example of modelling dislocation glide in cylindrical cells, let's look at the Diamond Screw dislocation in Si, using the [Holland and Marder](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.80.746) potential:"
    ]
   },
   {
@@ -22,19 +22,223 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from matscipy.dislocation import get_elastic_constants, BCCScrew111Dislocation\n",
+    "from matscipy.dislocation import get_elastic_constants, DiamondGlideScrew\n",
     "# the calculator to provide forces and energies from the potential\n",
-    "from matscipy.calculators.eam import EAM\n",
-    "eam_calc = EAM(\"../../tests/w_eam4.fs\")\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",
     "# the function accepts any ASE type of calculator\n",
-    "alat, C11, C12, C44 = get_elastic_constants(calculator=eam_calc, symbol=\"W\", verbose=False)\n",
+    "alat, C11, C12, C44 = get_elastic_constants(calculator=calc, symbol=\"Si\", verbose=False)\n",
     "\n",
-    "W_screw = BCCScrew111Dislocation(alat, C11, C12, C44, symbol=\"W\")\n",
+    "Si_screw = DiamondGlideScrew(alat, C11, C12, C44, symbol=\"Si\")\n",
     "\n",
-    "bulk, screw = W_screw.build_dislocation_cylinder(radius=30)\n",
-    "W_screw_bulk, W_screw_dislo = W_screw.build_cylinder(radius=20)\n",
-    "W_screw.view_cyl(W_screw_dislo)\n"
+    "Si_screw_bulk, Si_screw_dislo = Si_screw.build_cylinder(radius=20)\n",
+    "Si_screw.view_cyl(Si_screw_dislo)\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To generate the iniial and final structures for a glide event, we can use the `build_glide_configurations` function. We can also then use the `ase.neb` tools to smoothly interpolate an approximate glide path:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "bulk, glide_ini, glide_fin = Si_screw.build_glide_configurations(radius=20)\n",
+    "\n",
+    "from ase.neb import NEB\n",
+    "\n",
+    "nims = 5\n",
+    "\n",
+    "glide_neb = NEB([glide_ini.copy() for i in range(nims-1)] + [glide_fin.copy()])\n",
+    "\n",
+    "glide_neb.interpolate(method=\"idpp\", apply_constraint=True)"
+   ]
+  },
+  {
+   "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, radii=None):\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, 10))\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, radii=radii)\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_neb.images, radii=0.3)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This is also possible in the dissociated case:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "bulk, glide_ini, glide_fin = Si_screw.build_glide_configurations(radius=20, partial_distance=5)\n",
+    "\n",
+    "nims = 5\n",
+    "\n",
+    "glide_neb = NEB([glide_ini.copy() for i in range(nims-1)] + [glide_fin.copy()])\n",
+    "\n",
+    "glide_neb.interpolate(method=\"idpp\", apply_constraint=True)\n",
+    "\n",
+    "animate_glide(glide_neb.images, radii=0.3)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Glide in Dislocation 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": [
+    "from matscipy.dislocation import Quadrupole\n",
+    "\n",
+    "Si_quad = Quadrupole(DiamondGlideScrew, alat, C11, C12, C44, symbol=\"Si\")\n",
+    "\n",
+    "num_images = 5\n",
+    "\n",
+    "glide_quads = Si_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)\n",
+    "\n",
+    "animate_glide(glide_quads)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "num_images = 5\n",
+    "\n",
+    "glide_quads = Si_quad.build_glide_quadrupoles(nims=num_images, \n",
+    "                                            glide_left=False, # Prevent left dislocation glide\n",
+    "                                            glide_right=True, # Allow right dislocation glide\n",
+    "                                            partial_distance=2, # Dissociated Glide\n",
+    "                                            glide_separation=6,\n",
+    "                                            verbose=False)\n",
+    "\n",
+    "animate_glide(glide_quads)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Dislocation Kink\n",
+    "Dislocation kink is often the preferred mechanism for migration of the dislocation line. It involves a short segment of the dislocation line hopping by one glide vector, which then provides a low barrier for the rest of the dislocation line to migrate.\n",
+    "\n",
+    "Here the space of structures to explore is greater, due to the possibility of entire networks of kinks forming, and also due to the possibility for the dislocation to kink out by more than one glide.\n",
+    "\n",
+    "### Kink maps\n",
+    "Dislocation kink in `matscipy` is controlled by a kink map, which controls the position of the dislocation line as a function of the vertical coordinate. The kink map is a list of integers, where the values give a line position in units of glide distances, and each element corresponds to an extra cell in z.\n",
+    "\n",
+    "A kink map of `[0, 0, 1, 1]` means that the dislocation line initially starts at position zero for two repetitions, and then kinks out by one glide distance for two repetitions.\n",
+    "\n",
+    "Both dislocation cylinders and quadrupoles have support for this kind of kink map, but the boundaries are treated differently. In dislocation cylinders (using `build_kink_cyl`), the dislocation line is wrapped back around, and so there is an extra kink back to position zero in the generated structures. In quadrupoles (using `build_kink_quadrupole_network`), the cell is modified such that the dislocation position at the top of the cell becomes the new dislocation line position at the bottom of the cell. Therefore, in dislocation qudrupoles there would be no boundary.\n",
+    "\n",
+    "Since the kink map is just a list of integers, it can be more convenient to exploit the concatenation overloading of list addition, and specify kink maps in a form similar to `[0] * 5 + [1] * 5`.\n",
+    "\n",
+    "### Kinks in cylinders"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "Si_kink = Si_screw.build_kink_cyl(\n",
+    "                    radius=20,\n",
+    "                    kink_map= [0] * 5 + [1] * 5\n",
+    "                )\n",
+    "\n",
+    "Si_screw.view_cyl(Si_kink)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Kink in quadrupoles"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "Si_quad_kink = Si_quad.build_kink_quadrupole_network(\n",
+    "    glide_separation=8,\n",
+    "    kink_map=[0]*5 + [1]*5\n",
+    ")\n",
+    "\n",
+    "Si_quad.view_quad(Si_quad_kink)"
    ]
   },
   {
@@ -46,8 +250,22 @@
   }
  ],
  "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
   "language_info": {
-   "name": "python"
+   "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,

From 6afd0cf2b9c4b604fab6e62e5018d533e60db801 Mon Sep 17 00:00:00 2001
From: Tom <thomas.rocke@warwick.ac.uk>
Date: Wed, 7 Aug 2024 06:50:39 +0100
Subject: [PATCH 21/46] MAINT: classmethod property -> classproperty

---
 matscipy/dislocation.py | 63 +++++++++++++++++++++++++++++------------
 matscipy/utils.py       | 54 ++++++++++++++++++++++-------------
 2 files changed, 80 insertions(+), 37 deletions(-)

diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index 1c4249bb..4626fa46 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -45,7 +45,7 @@
 from matscipy.elasticity import fit_elastic_constants
 from matscipy.elasticity import Voigt_6x6_to_full_3x3x3x3
 from matscipy.elasticity import cubic_to_Voigt_6x6, coalesce_elastic_constants
-from matscipy.utils import validate_cubic_cell, points_in_polygon2D
+from matscipy.utils import validate_cubic_cell, points_in_polygon2D, classproperty
 
 
 def make_screw_cyl(alat, C11, C12, C44,
@@ -2610,7 +2610,7 @@ def _build_bulk_cyl(self, radius, core_positions, fix_rad, extension,
             core_positions[i, :] + extension[i, :], for the 
             purposes of adding extra atoms
         cyl_mask: array of bool
-            Optional ovverride for atomic mask to convert supercell into cyl
+            Optional override for atomic mask to convert supercell into cyl
 
         Returns
         -------
@@ -3137,7 +3137,27 @@ def _smooth_kink_displacements(self, kink1, kink2, ref_bulk, smoothing_width, ba
 
         return kink
     
-    def build_kink_cyl(self, kink_map=[0, 1], smooth_width=None, *args, **kwargs):
+    def build_kink_cyl(self, kink_map=None, smooth_width=None, *args, **kwargs):
+        """
+        Build a cylindrical cell with a dislocation kink network, defined by kink_map
+
+        kink_map: iterable of ints
+            Map of the location of the dislocation core in units of the glide vector
+            Default is a kink map of [0, 1]
+            See examples for more details.
+        smooth_width: float
+            Size (in Ang) of the region for displacement smoothing at each kink site.
+            Larger smoothing width assumes a broader kink structure.
+            Default is 0.5 * self.unit_cell.cell[2, 2]
+
+        *args, **kwargs
+            Extra arguments sent to self.build_cylinder
+        
+        """
+        # Deal with default kink_map value
+        if kink_map is None:
+            kink_map = [0, 1]
+
         kmap = np.array(kink_map, dtype=int)
         kmap -= np.min(kmap)
         krange = np.max(kmap)
@@ -3431,39 +3451,46 @@ def __init__(self, a, C11=None, C12=None, C44=None, C=None, symbol="W"):
     # Used so e.g. cls.crystalstructure is setup prior to __init__
     # as is the case with CubicCrystalDislocation
 
-    @classmethod
-    @property
+    # @classmethod
+    # @property
+    @classproperty
     def crystalstructure(cls):
         return cls.left_dislocation.crystalstructure
     
-    @classmethod
-    @property
+    # @classmethod
+    # @property
+    @classproperty
     def axes(cls):
         return cls.left_dislocation.axes
     
-    @classmethod
-    @property
+    # @classmethod
+    # @property
+    @classproperty
     def unit_cell_core_position_dimensionless(cls):
         return cls.left_dislocation.unit_cell_core_position_dimensionless
     
-    @classmethod
-    @property
+    # @classmethod
+    # @property
+    @classproperty
     def parity(cls):
         return cls.left_dislocation.parity
     
     
-    @classmethod
-    @property
+    # @classmethod
+    # @property
+    @classproperty
     def n_planes(cls):
         return cls.left_dislocation.n_planes
     
-    @classmethod
-    @property
+    # @classmethod
+    # @property
+    @classproperty
     def self_consistent(cls):
         return cls.left_dislocation.self_consistent
     
-    @classmethod
-    @property
+    # @classmethod
+    # @property
+    @classproperty
     def glide_distance_dimensionless(cls):
         return cls.left_dislocation.glide_distance_dimensionless
     
@@ -4046,7 +4073,7 @@ def _get_kink_quad_mask(self, ref_bulk, direction, precision=3):
             # Check if tilt_bulk has same 1st, 2nd, & 3rd neighbour coordination
             # as normal bulk
             # TODO: This currently does not work for BCCEdge111barDislocation
-            if np.product([
+            if np.prod([
                 all(coordination(tilt_bulk[full_mask], cutoff * self.alat) == bulk_coord[j])
                 for j, cutoff in enumerate(cutoffs)
                 ]):
diff --git a/matscipy/utils.py b/matscipy/utils.py
index ff2175a8..a86046e8 100644
--- a/matscipy/utils.py
+++ b/matscipy/utils.py
@@ -1,9 +1,24 @@
+import functools
 import warnings
 import numpy as np
 from ase.utils.structure_comparator import SymmetryEquivalenceCheck
 
-def validate_cubic_cell(a, symbol="w", axes=None, crystalstructure=None, pbc=True):
+
+class classproperty:
+    '''
+    Decorator class to replace classmethod property decorators
     '''
+    def __init__(self, method):
+        self.method = method
+        functools.update_wrapper(self, method)
+    def __get__(self, obj, cls=None):
+        if cls is None:
+            cls = type(obj)
+        return self.method(cls)
+
+
+def validate_cubic_cell(a, symbol="w", axes=None, crystalstructure=None, pbc=True):
+    """
     Provide uniform interface for generating rotated atoms objects through two main methods:
 
     For lattice constant + symbol + crystalstructure, build a valid cubic bulk rotated into the frame defined by axes
@@ -20,7 +35,7 @@ def validate_cubic_cell(a, symbol="w", axes=None, crystalstructure=None, pbc=Tru
         Base Structure of bulk system, currently supported: fcc, bcc, diamond
     pbc: list of bool
         Periodic Boundary Conditions in x, y, z
-    '''
+    """
 
     from ase.lattice.cubic import FaceCenteredCubic, BodyCenteredCubic, Diamond
     from ase.atoms import Atoms
@@ -121,7 +136,7 @@ def validate_cubic_cell(a, symbol="w", axes=None, crystalstructure=None, pbc=Tru
 
 
 def find_condensed_repr(atoms, precision=2):
-    '''
+    """
     Search for a condensed representation of atoms
     Equivalent to attempting to undo a supercell
 
@@ -131,7 +146,7 @@ def find_condensed_repr(atoms, precision=2):
         Number of decimal places to use in determining whether scaled positions are equal
 
     returns a condensed copy of atoms, if such a condensed representation is found. Else returns a copy of atoms
-    '''
+    """
     ats = atoms.copy()
     for axis in range(2, -1, -1):
         ats = find_condensed_repr_along_axis(ats, axis, precision)
@@ -139,7 +154,7 @@ def find_condensed_repr(atoms, precision=2):
     return ats
 
 def find_condensed_repr_along_axis(atoms, axis=-1, precision=2):
-    '''
+    """
     Search for a condensed representation of atoms about axis. 
     Essentially an inverse to taking a supercell.
 
@@ -152,7 +167,7 @@ def find_condensed_repr_along_axis(atoms, axis=-1, precision=2):
         Number of decimal places to use in determining whether scaled positions are equal
 
     returns a condensed copy of atoms, if such a condensed representation is found. Else returns a copy of atoms
-    '''
+    """
     comp = SymmetryEquivalenceCheck()
 
     cart_directions = [
@@ -172,7 +187,8 @@ def find_condensed_repr_along_axis(atoms, axis=-1, precision=2):
     # Find all atoms which are in line with the origin_idx th atom
     p = np.round(ats.get_scaled_positions(), precision)
 
-    p_diff = (p - p[origin_idx, :]) % 1
+    # MIC distance vector from origin atom, in scaled coordinates
+    p_diff = (p - p[origin_idx, :]) % 1.0
 
     matches = np.argwhere((p_diff[:, dirs[0]] < 10**(-precision)) * (p_diff[:, dirs[1]] < 10**(-precision)))[:, 0]
     min_off = np.inf
@@ -212,7 +228,7 @@ def find_condensed_repr_along_axis(atoms, axis=-1, precision=2):
 
 
 def complete_basis(v1, v2=None, normalise=False, nmax=5, tol=1E-6):
-    '''
+    """
     Generate a complete (v1, v2, v3) orthogonal basis in 3D from v1 and an optional v2
 
     (V1, V2, V3) is always right-handed.
@@ -236,7 +252,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):
         for i in range(nmax):
@@ -349,7 +365,7 @@ def get_structure_types(structure, diamond_structure=False):
     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
@@ -358,7 +374,7 @@ def line_intersect_2D(p1, p2, x1, x2):
 
     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]
@@ -402,7 +418,7 @@ def line_intersect_2D(p1, p2, x1, x2):
     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
@@ -416,7 +432,7 @@ def points_in_polygon2D(p, poly_points):
     -------
     mask: np.array
         Boolean mask. True for points inside the poylgon
-    '''
+    """
     if len(p.shape) == 1:
         # Single point provided
         points = p.copy()[np.newaxis, :]
@@ -448,7 +464,7 @@ def points_in_polygon2D(p, poly_points):
     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)
     
@@ -458,13 +474,13 @@ def get_distance_from_polygon2D(test_points:np.array, polygon_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:
@@ -490,7 +506,7 @@ def get_dist(p, v, w):
     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
@@ -501,7 +517,7 @@ def radial_mask_from_polygon2D(test_points:np.array, polygon_points:np.array, ra
         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)
 

From 5da21cb78af23159c9cbdab7a77a53849f162657 Mon Sep 17 00:00:00 2001
From: Tom <thomas.rocke@warwick.ac.uk>
Date: Wed, 7 Aug 2024 06:51:03 +0100
Subject: [PATCH 22/46] DOC: Tweaks to disloc mobility docs

---
 docs/applications/disloc_mobility.ipynb | 82 ++++++++++++++++++++-----
 1 file changed, 68 insertions(+), 14 deletions(-)

diff --git a/docs/applications/disloc_mobility.ipynb b/docs/applications/disloc_mobility.ipynb
index 214c59ed..98890c76 100644
--- a/docs/applications/disloc_mobility.ipynb
+++ b/docs/applications/disloc_mobility.ipynb
@@ -5,11 +5,14 @@
    "metadata": {},
    "source": [
     "# Modelling the Mobility of Dislocations with `matscipy`\n",
-    "\n",
-    "The dislocation modelling toolkit provided in `matscipy.dislocation` includes functions to assist in characterising the mobility of dislocations in various systems. The two supported mechanisms are glide and kink, which both allow the dislocation to move in the glide direction.\n",
+    "The dislocation modelling toolkit provided in `matscipy.dislocation` includes functions to assist in characterising the mobility of dislocations in various systems. Currently, only motion in the glide direction is supported.\n",
     "\n",
     "## Dislocation Glide\n",
-    "Dislocation Glide is where the entire dislocation line migrates across a glide barrier to advance the dislocation. These kinds of structures are easy to generate both in cylindrical and in quadrupolar structures.\n",
+    "Dislocation Glide is where the dislocation line migrates in the glide direction. There are two main mechanisms for attempting to describe dislocation glide dynamics.\n",
+    "\n",
+    "The first mechanism is the (double) kink mechanism. Starting from a perfectly straight dislocation line, a small segment  nucleates out in the glide directionforming a pair of dislocation kinks. These kinks can then migrate along the dislocation line.\n",
+    "\n",
+    "A second, simpler mechanism is to consider the entire dislocation line moving at once - this is much less physical than the double kink mechanism, but can be a good first approximation. The 2D nature of this mechanism can also make it much simpler to study, rather than having to deal with the fully 3D kinks.\n",
     "\n",
     "\n",
     "### Glide in Dislocation Cylinders\n",
@@ -18,9 +21,50 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 1,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "273fa26b287d4c61a24ff07d20735b7f",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": []
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/home/europa/.local/lib/python3.10/site-packages/matscipy/dislocation.py:484: FutureWarning: Import StrainFilter from ase.filters\n",
+      "  sf = StrainFilter(unit_cell)\n"
+     ]
+    },
+    {
+     "ename": "ImportError",
+     "evalue": "libQt6Core.so.6: cannot open shared object file: No such file or directory",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mImportError\u001b[0m                               Traceback (most recent call last)",
+      "Cell \u001b[0;32mIn[1], line 14\u001b[0m\n\u001b[1;32m     11\u001b[0m Si_screw \u001b[38;5;241m=\u001b[39m DiamondGlideScrew(alat, C11, C12, C44, symbol\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mSi\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m     13\u001b[0m Si_screw_bulk, Si_screw_dislo \u001b[38;5;241m=\u001b[39m Si_screw\u001b[38;5;241m.\u001b[39mbuild_cylinder(radius\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m20\u001b[39m)\n\u001b[0;32m---> 14\u001b[0m \u001b[43mSi_screw\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mview_cyl\u001b[49m\u001b[43m(\u001b[49m\u001b[43mSi_screw_dislo\u001b[49m\u001b[43m)\u001b[49m\n",
+      "File \u001b[0;32m~/.local/lib/python3.10/site-packages/matscipy/dislocation.py:3306\u001b[0m, in \u001b[0;36mCubicCrystalDislocation.view_cyl\u001b[0;34m(system, scale, CNA_color, add_bonds, line_color, disloc_names, hide_arrows)\u001b[0m\n\u001b[1;32m   3303\u001b[0m \u001b[38;5;66;03m# Check if system contains a diamond dislocation (e.g. DiamondGlideScrew, Diamond90degreePartial)\u001b[39;00m\n\u001b[1;32m   3304\u001b[0m diamond_structure \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mDiamond\u001b[39m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;129;01min\u001b[39;00m system\u001b[38;5;241m.\u001b[39minfo[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mdislocation_classes\u001b[39m\u001b[38;5;124m\"\u001b[39m][\u001b[38;5;241m0\u001b[39m]\n\u001b[0;32m-> 3306\u001b[0m atom_labels, structure_names, colors \u001b[38;5;241m=\u001b[39m \u001b[43mget_structure_types\u001b[49m\u001b[43m(\u001b[49m\u001b[43msystem\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\n\u001b[1;32m   3307\u001b[0m \u001b[43m                                                        \u001b[49m\u001b[43mdiamond_structure\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdiamond_structure\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m   3308\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m disloc_names \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m   3309\u001b[0m     disloc_names \u001b[38;5;241m=\u001b[39m system\u001b[38;5;241m.\u001b[39minfo[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mdislocation_types\u001b[39m\u001b[38;5;124m\"\u001b[39m]\n",
+      "File \u001b[0;32m~/.local/lib/python3.10/site-packages/matscipy/utils.py:342\u001b[0m, in \u001b[0;36mget_structure_types\u001b[0;34m(structure, diamond_structure)\u001b[0m\n\u001b[1;32m    331\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mget_structure_types\u001b[39m(structure, diamond_structure\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mFalse\u001b[39;00m):\n\u001b[1;32m    332\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"Get the results of Common Neighbor Analysis and \u001b[39;00m\n\u001b[1;32m    333\u001b[0m \u001b[38;5;124;03m        Identify Diamond modifiers from Ovito\u001b[39;00m\n\u001b[1;32m    334\u001b[0m \u001b[38;5;124;03m        (Requires Ovito python module)\u001b[39;00m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    340\u001b[0m \u001b[38;5;124;03m        colors (list of strings): colors of the structure types in hex format\u001b[39;00m\n\u001b[1;32m    341\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[0;32m--> 342\u001b[0m     \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01movito\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mio\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mase\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m ase_to_ovito\n\u001b[1;32m    343\u001b[0m     \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01movito\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mmodifiers\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m CommonNeighborAnalysisModifier, IdentifyDiamondModifier\n\u001b[1;32m    344\u001b[0m     \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01movito\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mpipeline\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m StaticSource, Pipeline\n",
+      "File \u001b[0;32m~/.local/lib/python3.10/site-packages/ovito/__init__.py:21\u001b[0m\n\u001b[1;32m     18\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01movito\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mnonpublic\u001b[39;00m\n\u001b[1;32m     20\u001b[0m \u001b[38;5;66;03m# Load all extension modules.\u001b[39;00m\n\u001b[0;32m---> 21\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01movito\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01m_extensions\u001b[39;00m\n\u001b[1;32m     23\u001b[0m \u001b[38;5;66;03m# The public symbols of this root module:\u001b[39;00m\n\u001b[1;32m     24\u001b[0m __all__ \u001b[38;5;241m=\u001b[39m [\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mversion\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mversion_string\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mscene\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mScene\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124menable_logging\u001b[39m\u001b[38;5;124m'\u001b[39m]\n",
+      "File \u001b[0;32m~/.local/lib/python3.10/site-packages/ovito/_extensions/__init__.py:6\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[38;5;66;03m# Load all plugin extension modules under the _extensions package.\u001b[39;00m\n\u001b[1;32m      5\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m modinfo \u001b[38;5;129;01min\u001b[39;00m pkgutil\u001b[38;5;241m.\u001b[39mwalk_packages(__path__, \u001b[38;5;18m__name__\u001b[39m \u001b[38;5;241m+\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124m.\u001b[39m\u001b[38;5;124m'\u001b[39m):\n\u001b[0;32m----> 6\u001b[0m     \u001b[43mimportlib\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mimport_module\u001b[49m\u001b[43m(\u001b[49m\u001b[43mmodinfo\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mname\u001b[49m\u001b[43m)\u001b[49m\n",
+      "File \u001b[0;32m/usr/lib/python3.10/importlib/__init__.py:126\u001b[0m, in \u001b[0;36mimport_module\u001b[0;34m(name, package)\u001b[0m\n\u001b[1;32m    124\u001b[0m             \u001b[38;5;28;01mbreak\u001b[39;00m\n\u001b[1;32m    125\u001b[0m         level \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m \u001b[38;5;241m1\u001b[39m\n\u001b[0;32m--> 126\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43m_bootstrap\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_gcd_import\u001b[49m\u001b[43m(\u001b[49m\u001b[43mname\u001b[49m\u001b[43m[\u001b[49m\u001b[43mlevel\u001b[49m\u001b[43m:\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mpackage\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mlevel\u001b[49m\u001b[43m)\u001b[49m\n",
+      "File \u001b[0;32m~/.local/lib/python3.10/site-packages/ovito/_extensions/anari.py:2\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[38;5;66;03m# Load dependencies.\u001b[39;00m\n\u001b[0;32m----> 2\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01movito\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01m_extensions\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mpyscript\u001b[39;00m\n\u001b[1;32m      4\u001b[0m \u001b[38;5;66;03m# Load the C extension module.\u001b[39;00m\n\u001b[1;32m      5\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01movito\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mplugins\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mAnariRendererPython\u001b[39;00m\n",
+      "File \u001b[0;32m~/.local/lib/python3.10/site-packages/ovito/_extensions/pyscript.py:2\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[38;5;66;03m# Load the C extension module.\u001b[39;00m\n\u001b[0;32m----> 2\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01movito\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mplugins\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mPyScript\u001b[39;00m\n\u001b[1;32m      4\u001b[0m \u001b[38;5;66;03m# Load class add-ons.\u001b[39;00m\n\u001b[1;32m      5\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01movito\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01m_scene_class\u001b[39;00m\n",
+      "File \u001b[0;32m~/.local/lib/python3.10/site-packages/ovito/plugins/__init__.py:36\u001b[0m\n\u001b[1;32m     26\u001b[0m     warnings\u001b[38;5;241m.\u001b[39mwarn(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mIncompatible version of the Qt cross-platform framework detected!\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124mThis version of the OVITO Python module is based on Qt6 (loaded via the PySide6 bindings module), \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m     27\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mbut bindings for old Qt5 are already loaded at this point (through PyQt5 or PySide2 imports preceding the import of OVITO). To avoid library version conflicts, please make sure the rest of \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m     28\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124myour application uses Qt6 too instead of Qt5. \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m     31\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mIn addition, it may help to set the environment variable QT_API=pyside6 to force third-party packages (e.g. matplotlib) to load Qt6 instead of Qt5. \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m     32\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mIf you have any questions, please contact support@ovito.org.\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m     34\u001b[0m \u001b[38;5;66;03m# Load all the Qt bindings first before OVITO's own C++ modules get loaded.\u001b[39;00m\n\u001b[1;32m     35\u001b[0m \u001b[38;5;66;03m# This ensures that the right Qt shared libraries needed by OVITO are already loaded into the process when running in a system Python interpreter.\u001b[39;00m\n\u001b[0;32m---> 36\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01movito\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mqt_compat\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m QtCore\n\u001b[1;32m     37\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01movito\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mqt_compat\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m QtGui\n\u001b[1;32m     38\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01movito\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mqt_compat\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m QtWidgets\n",
+      "File \u001b[0;32m~/.local/lib/python3.10/site-packages/ovito/qt_compat.py:11\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[38;5;66;03m# Import Shiboken module first to preload shiboken6.dll, which is needed by the Qt DLLs.\u001b[39;00m\n\u001b[1;32m      9\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mshiboken6\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mshiboken\u001b[39;00m \n\u001b[0;32m---> 11\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mPySide6\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m QtCore, QtGui, QtWidgets, QtNetwork, QtXml\n\u001b[1;32m     13\u001b[0m \u001b[38;5;66;03m# The QtOpenGLWidgets module is not included in the embedded PySide6 installation of OVITO Pro 3.7,\u001b[39;00m\n\u001b[1;32m     14\u001b[0m \u001b[38;5;66;03m# but it should be loaded when using the PyPI version of PySide6.\u001b[39;00m\n\u001b[1;32m     15\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n",
+      "\u001b[0;31mImportError\u001b[0m: libQt6Core.so.6: cannot open shared object file: No such file or directory"
+     ]
+    }
+   ],
    "source": [
     "from matscipy.dislocation import get_elastic_constants, DiamondGlideScrew\n",
     "# the calculator to provide forces and energies from the potential\n",
@@ -42,7 +86,9 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "To generate the iniial and final structures for a glide event, we can use the `build_glide_configurations` function. We can also then use the `ase.neb` tools to smoothly interpolate an approximate glide path:"
+    "To generate the initial and final structures for a full glide event, we can use the `build_glide_configurations` function. The structures contain straight dislocation lines, separated by the dislocation glide distance. These structures will be similar to those produced by a call to `build_cylinder`, except extra bulk is added (creating a pill shape, not a circle) in order to make the initial and final configurations symmetric.\n",
+    "\n",
+    "We can also then use the `ase.neb` tools to smoothly interpolate an approximate glide path, which allows us to generate structures for the simpler dislocation glide mechanism discussed above."
    ]
   },
   {
@@ -141,9 +187,9 @@
    "metadata": {},
    "source": [
     "### Glide in Dislocation 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",
+    "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.\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.`"
+    "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."
    ]
   },
   {
@@ -192,16 +238,24 @@
     "## Dislocation Kink\n",
     "Dislocation kink is often the preferred mechanism for migration of the dislocation line. It involves a short segment of the dislocation line hopping by one glide vector, which then provides a low barrier for the rest of the dislocation line to migrate.\n",
     "\n",
-    "Here the space of structures to explore is greater, due to the possibility of entire networks of kinks forming, and also due to the possibility for the dislocation to kink out by more than one glide.\n",
+    "Here the space of structures to explore is greater, due to the 3D nature of dislocation kink. Kink nucleation is also possible on segments of an already kinked-out dislocation line, which can lead to a full network of dislocation kinks.  \n",
     "\n",
     "### Kink maps\n",
-    "Dislocation kink in `matscipy` is controlled by a kink map, which controls the position of the dislocation line as a function of the vertical coordinate. The kink map is a list of integers, where the values give a line position in units of glide distances, and each element corresponds to an extra cell in z.\n",
+    "Dislocation kink in `matscipy` is controlled by a kink map, which controls the position of the dislocation line as a function of the vertical coordinate. The kink map is a list of integers, where the values give a line position in units of glide distances, and each element corresponds to an extra cell in z (replication of `disloc.unit_cell`).\n",
+    "\n",
+    "A kink map of `[0, 0, 1, 1]` means that the dislocation line initially starts at position zero for two repetitions, and then moves out by one glide distance for two repetitions.\n",
+    "\n",
+    "Both dislocation cylinders and quadrupoles have support for this kind of kink map, but the periodic boundaries are treated differently. \n",
+    "\n",
+    "In dislocation cylinders (using `build_kink_cyl`), the periodicity in z is enforced by requiring that the dislocation line returns back to the initial position across the periodic boundary. The example kink map of `[0, 0, 1, 1]` will therefore have a single kink out in the center of the cell, and a single kink back at the periodic boundary (the dislocation line will go like `0, 0, 1, 1, 0, 0, 1, 1, ...`).\n",
     "\n",
-    "A kink map of `[0, 0, 1, 1]` means that the dislocation line initially starts at position zero for two repetitions, and then kinks out by one glide distance for two repetitions.\n",
+    "In quadrupoles (using `build_kink_quadrupole_network`), the cell is modified such that the dislocation position at the top of the cell becomes the new dislocation line position at the bottom of the cell. This means that for quadrupoles an extra kink will not be created, and that the example map of `[0, 0, 1, 1]` will create only one kink in the center of the cell.\n",
     "\n",
-    "Both dislocation cylinders and quadrupoles have support for this kind of kink map, but the boundaries are treated differently. In dislocation cylinders (using `build_kink_cyl`), the dislocation line is wrapped back around, and so there is an extra kink back to position zero in the generated structures. In quadrupoles (using `build_kink_quadrupole_network`), the cell is modified such that the dislocation position at the top of the cell becomes the new dislocation line position at the bottom of the cell. Therefore, in dislocation qudrupoles there would be no boundary.\n",
+    ":::{note}\n",
+    "For the cell shift to always produce the correct crystalstructure, some atoms need to be deleted for some values of `kink_map[-1]`. This means that some kink maps may not conserve stoichiometry for some multispecies systems.\n",
+    ":::\n",
     "\n",
-    "Since the kink map is just a list of integers, it can be more convenient to exploit the concatenation overloading of list addition, and specify kink maps in a form similar to `[0] * 5 + [1] * 5`.\n",
+    "Since the kink map is just a list of integers, it can be more convenient to exploit list addition, and specify kink maps in a form similar to `[0] * 5 + [1] * 5`, which would be identical to the input of `[0, 0, 0, 0, 0, 1, 1, 1, 1, 1]`.\n",
     "\n",
     "### Kinks in cylinders"
    ]
@@ -265,7 +319,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.12"
+   "version": "3.10.14"
   }
  },
  "nbformat": 4,

From 9cabb6e8f3247aa0e8721e255aa8c3633088048d Mon Sep 17 00:00:00 2001
From: Tom <thomas.rocke@warwick.ac.uk>
Date: Wed, 7 Aug 2024 06:51:39 +0100
Subject: [PATCH 23/46] DOC: Proper note syntax in multispecies dislocs docs

---
 docs/applications/multispecies_dislocations.ipynb | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/docs/applications/multispecies_dislocations.ipynb b/docs/applications/multispecies_dislocations.ipynb
index 7236e480..b9e20ce7 100644
--- a/docs/applications/multispecies_dislocations.ipynb
+++ b/docs/applications/multispecies_dislocations.ipynb
@@ -149,7 +149,9 @@
     "\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 built by the dislocaion classes."
+    ":::{note}\n",
+    "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 dislocation classes.\n",
+    ":::"
    ]
   },
   {

From 1c9239963cbe4d86474cbeea3654ecf9c057372a Mon Sep 17 00:00:00 2001
From: Tom <thomas.rocke@warwick.ac.uk>
Date: Wed, 7 Aug 2024 16:27:21 +0100
Subject: [PATCH 24/46] ENH: Separate kink cyl gen into distinct steps

---
 matscipy/dislocation.py | 63 +++++++++++++++++++++++++++++++++++++----
 1 file changed, 58 insertions(+), 5 deletions(-)

diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index 4626fa46..56c471bd 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -3137,13 +3137,14 @@ def _smooth_kink_displacements(self, kink1, kink2, ref_bulk, smoothing_width, ba
 
         return kink
     
-    def build_kink_cyl(self, kink_map=None, smooth_width=None, *args, **kwargs):
+    def build_kink_glide_structs(self, kink_map=None, *args, **kwargs):
         """
-        Build a cylindrical cell with a dislocation kink network, defined by kink_map
+        Build a reference bulk and a minimal set of glide structures required to build a kink structure of the given kink map (default [0, 1]).
+        Does not actually build the kink structure, this is handled by self.kink_from_glide_cyls.
 
-        kink_map: iterable of ints
+        kink_map: iterable of ints or None
             Map of the location of the dislocation core in units of the glide vector
-            Default is a kink map of [0, 1]
+            Default (kink_map=None) is a kink map of [0, 1]
             See examples for more details.
         smooth_width: float
             Size (in Ang) of the region for displacement smoothing at each kink site.
@@ -3162,6 +3163,8 @@ def build_kink_cyl(self, kink_map=None, smooth_width=None, *args, **kwargs):
         kmap -= np.min(kmap)
         krange = np.max(kmap)
 
+        assert len(kmap) > 0
+
         unique_kinks = np.sort(np.unique(kmap))
 
         glide_structs = {}
@@ -3176,9 +3179,38 @@ def build_kink_cyl(self, kink_map=None, smooth_width=None, *args, **kwargs):
                                                           fixed_points=fixed_points, 
                                                           core_position=np.array([kink_pos * self.glide_distance, 0 , 0]),
                                                           **kwargs)
+        return ref_bulk, glide_structs
+    
+    def kink_from_glide_cyls(self, kink_map, ref_bulk, glide_structs, smooth_width=None):
+        """
+        Build a kink cylinder cell from the given kink map, using the structures contained in glide_structs
+
+        kink_map: iterable of ints or None
+            Map of the location of the dislocation core in units of the glide vector
+            Default (kink_map=None) is a kink map of [0, 1]
+            See examples for more details.
+        ref_bulk: ase Atoms object
+            Reference bulk structure, as returned by self.build_kink_glide_structs
+        glide_structs: list of ase Atoms
+            Glide structures e.g. those produced by self.build_kink_glide_structs.
+            The kink structure is constructed using glide_structs[kink_map[i]] for the ith cell.
+        smooth_width: float
+            Size (in Ang) of the region for displacement smoothing at each kink site.
+            Larger smoothing width assumes a broader kink structure.
+            Default is 0.5 * self.unit_cell.cell[2, 2]
         
-        kink_structs = [glide_structs[midx].copy() for midx in kmap]
+        """
+        # Deal with default kink_map value
+        if kink_map is None:
+            kink_map = [0, 1]
+
+        kmap = np.array(kink_map, dtype=int)
+        kmap -= np.min(kmap)
 
+        assert np.max(kmap) == len(glide_structs)
+
+        kink_structs = [glide_structs[midx].copy() for midx in kmap]
+        
         # Smooth displacements across kink boundaries
         if smooth_width is None:
             smooth_width = 0.5 * self.unit_cell.cell[2, 2]
@@ -3198,6 +3230,27 @@ def build_kink_cyl(self, kink_map=None, smooth_width=None, *args, **kwargs):
 
         return kink_cyl
     
+    def build_kink_cyl(self, kink_map=None, smooth_width=None, *args, **kwargs):
+        """
+        Build a cylindrical cell with a dislocation kink network, defined by kink_map
+
+        kink_map: iterable of ints or None
+            Map of the location of the dislocation core in units of the glide vector
+            Default (kink_map=None) is a kink map of [0, 1]
+            See examples for more details.
+        smooth_width: float
+            Size (in Ang) of the region for displacement smoothing at each kink site.
+            Larger smoothing width assumes a broader kink structure.
+            Default is 0.5 * self.unit_cell.cell[2, 2]
+
+        *args, **kwargs
+            Extra arguments sent to self.build_cylinder
+        
+        """
+        ref_bulk, glide_structs = self.build_kink_glide_structs(kink_map, *args, **kwargs)
+
+        return self.kink_from_glide_cyls(kink_map, ref_bulk, glide_structs, smooth_width)
+    
     @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):

From e4f6cbac35cac4490fe5e57c44854fb61c181e9a Mon Sep 17 00:00:00 2001
From: thomas-rocke <tom@rocke.co.uk>
Date: Sat, 17 Aug 2024 14:26:56 +0100
Subject: [PATCH 25/46] ENH: Quad kink in separate steps

---
 matscipy/dislocation.py | 224 ++++++++++++++++++++++++++++------------
 1 file changed, 158 insertions(+), 66 deletions(-)

diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index 56c471bd..c1b604dc 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -3146,15 +3146,12 @@ def build_kink_glide_structs(self, kink_map=None, *args, **kwargs):
             Map of the location of the dislocation core in units of the glide vector
             Default (kink_map=None) is a kink map of [0, 1]
             See examples for more details.
-        smooth_width: float
-            Size (in Ang) of the region for displacement smoothing at each kink site.
-            Larger smoothing width assumes a broader kink structure.
-            Default is 0.5 * self.unit_cell.cell[2, 2]
 
         *args, **kwargs
             Extra arguments sent to self.build_cylinder
         
         """
+
         # Deal with default kink_map value
         if kink_map is None:
             kink_map = [0, 1]
@@ -3167,7 +3164,8 @@ def build_kink_glide_structs(self, kink_map=None, *args, **kwargs):
 
         unique_kinks = np.sort(np.unique(kmap))
 
-        glide_structs = {}
+        # Make an empty list of the correct length
+        glide_structs = [0] * (krange+1)
 
         fixed_points = np.array([
             [0, 0, 0],
@@ -3175,13 +3173,37 @@ def build_kink_glide_structs(self, kink_map=None, *args, **kwargs):
         ])
 
         for kink_pos in unique_kinks:
+            print(kink_pos, krange)
             ref_bulk, glide_structs[kink_pos] = self.build_cylinder(*args, 
                                                           fixed_points=fixed_points, 
                                                           core_position=np.array([kink_pos * self.glide_distance, 0 , 0]),
                                                           **kwargs)
+            
+        # Deal with small displacements in boundary conditions caused by differences in disloc positions
+        # by taking the average of the extremes
+        fix_mask = glide_structs[0].arrays["fix_mask"]
+
+        ini_pos = glide_structs[0].positions[fix_mask, :]
+        fin_pos = glide_structs[-1].positions[fix_mask]
+        av_pos = (ini_pos + fin_pos) / 2
+
+        for struct in glide_structs:
+            # Guarantee that the fix masks are all identical
+            struct.arrays["fix_mask"] = fix_mask
+
+            p = struct.positions
+            p[fix_mask, :] = av_pos
+            struct.set_positions(p)
+
+            # Reapply the FixAtoms constraints
+            struct.set_constraint()
+            struct.set_constraint(
+                FixAtoms(mask=fix_mask)
+            )
+
         return ref_bulk, glide_structs
     
-    def kink_from_glide_cyls(self, kink_map, ref_bulk, glide_structs, smooth_width=None):
+    def build_kink_from_glide_cyls(self, kink_map, ref_bulk, glide_structs, smooth_width=None):
         """
         Build a kink cylinder cell from the given kink map, using the structures contained in glide_structs
 
@@ -3207,7 +3229,7 @@ def kink_from_glide_cyls(self, kink_map, ref_bulk, glide_structs, smooth_width=N
         kmap = np.array(kink_map, dtype=int)
         kmap -= np.min(kmap)
 
-        assert np.max(kmap) == len(glide_structs)
+        assert np.max(kmap) == len(glide_structs) - 1
 
         kink_structs = [glide_structs[midx].copy() for midx in kmap]
         
@@ -3225,12 +3247,20 @@ def kink_from_glide_cyls(self, kink_map, ref_bulk, glide_structs, smooth_width=N
             )
 
         kink_cyl = smoothed_kink_structs[0].copy()
+
         for i in range(1, N):
             kink_cyl = stack(kink_cyl, smoothed_kink_structs[i])
 
+        # Add the correct FixAtoms constraint by concatenating glide struct fix masks in the correct order
+        fix_mask = np.array(sum([list(glide_structs[midx].arrays["fix_mask"]) for midx in kmap], []), dtype=bool)
+        kink_cyl.arrays["fix_mask"] = fix_mask
+        kink_cyl.set_constraint(
+            FixAtoms(mask=fix_mask)
+        )
+
         return kink_cyl
     
-    def build_kink_cyl(self, kink_map=None, smooth_width=None, *args, **kwargs):
+    def build_kink_cylinder(self, kink_map=None, smooth_width=None, *args, **kwargs):
         """
         Build a cylindrical cell with a dislocation kink network, defined by kink_map
 
@@ -3249,7 +3279,7 @@ def build_kink_cyl(self, kink_map=None, smooth_width=None, *args, **kwargs):
         """
         ref_bulk, glide_structs = self.build_kink_glide_structs(kink_map, *args, **kwargs)
 
-        return self.kink_from_glide_cyls(kink_map, ref_bulk, glide_structs, smooth_width)
+        return self.build_kink_from_glide_cyls(kink_map, ref_bulk, glide_structs, smooth_width)
     
     @staticmethod
     def view_cyl(system, scale=0.5, CNA_color=True, add_bonds=False,
@@ -4077,9 +4107,8 @@ def _get_kink_quad_mask(self, ref_bulk, direction, precision=3):
             Direction of kink (+1 or -1)
         precision: int
             Scaled precision to use when determining a step size
-            Each step is self.alat * 10**_precision Ang
+            Each step is self.alat * 10**-precision Ang
             
-
         returns a boolean mask, and a new cell
 
         '''
@@ -4144,7 +4173,107 @@ def _get_kink_quad_mask(self, ref_bulk, direction, precision=3):
         # If we end up here, all layers removed before we found a match. Raise error
         raise RuntimeError("Could not find a valid periodic kink cell.")
 
-    def build_kink_quadrupole(self, z_reps=2, layer_decimal_precision=3, n_kink=1, invert_direction=False, smooth_width=None,
+    def build_kink_quadrupole_glide_structs(self, kink_map=None, *args, **kwargs):
+        """
+        Build a reference quadrupole bulk and a minimal set of glide structures required to build a kink structure of the given kink map (default [0, 1]).
+        Does not actually build the kink structure, this is handled by self.kink_from_glide_cyls.
+
+        kink_map: iterable of ints or None
+            Map of the location of the dislocation core in units of the glide vector
+            Default (kink_map=None) is a kink map of [0, 1]
+            See examples for more details.
+
+        *args, **kwargs
+            Extra arguments sent to self.build_quadrupole
+        """
+        if kink_map is None:
+            kink_map = [0, 1]
+
+        kmap = np.array(kink_map, dtype=int)
+        kmap -= np.min(kmap)
+        krange = np.max(kmap)
+
+        assert len(kmap) > 0
+
+        unique_kinks = np.sort(np.unique(kmap))
+
+        glide_structs = [0] * (krange+1)
+
+        for kink_pos in unique_kinks:
+            off = np.array([kink_pos * self.glide_distance, 0 , 0])
+            ref_bulk, glide_structs[kink_pos] = self.build_quadrupole(*args, 
+                                                          left_offset=off, right_offset=off,
+                                                          **kwargs)
+        return ref_bulk, glide_structs
+
+    def build_kink_quadrupole_from_glide_structs(self, kink_map, ref_bulk, glide_structs, smooth_width=None, layer_decimal_precision=3):
+        """
+        Build a kink cylinder cell from the given kink map, using the structures contained in glide_structs
+
+        kink_map: iterable of ints or None
+            Map of the location of the dislocation core in units of the glide vector
+            Default (kink_map=None) is a kink map of [0, 1]
+            See examples for more details.
+        ref_bulk: ase Atoms object
+            Reference bulk structure, as returned by self.build_kink_quadrupole_glide_structs
+        glide_structs: list of ase Atoms
+            Glide structures e.g. those produced by self.build_kink_quadrupole_glide_structs.
+            The kink structure is constructed using glide_structs[kink_map[i]] for the ith cell.
+        smooth_width: float
+            Size (in Ang) of the region for displacement smoothing at each kink site.
+            Larger smoothing width assumes a broader kink structure.
+            Default is 0.5 * self.unit_cell.cell[2, 2]
+        layer_decimal_precision: float
+            Tolerance for trying to determine the top atomic layer of the cell
+            (required in order to ensure complete periodicity)
+            Top layer is defined as any atom with z component >= max(atom_z) - 10**(-layer_tol)
+            in fractional coordinates
+        
+        """
+        # Deal with default kink_map value
+        if kink_map is None:
+            kink_map = [0, 1]
+
+        kmap = np.array(kink_map, dtype=int)
+        kmap -= np.min(kmap)
+
+        assert np.max(kmap) == len(glide_structs) - 1
+
+        kink_structs = [glide_structs[midx].copy() for midx in kmap]
+        
+        # Smooth displacements across kink boundaries
+        if smooth_width is None:
+            smooth_width = 0.5 * self.unit_cell.cell[2, 2]
+
+        N = len(kink_structs)
+
+        smoothed_kink_structs = []
+
+        for i in range(N):
+            smoothed_kink_structs.append(
+                self._smooth_kink_displacements(kink_structs[i], kink_structs[(i+1) % N], ref_bulk, smooth_width, ref_bulk.cell[:, :])
+            )
+
+        kink_quad = smoothed_kink_structs[0].copy()
+
+        for i in range(1, N):
+            kink_quad = stack(kink_quad, smoothed_kink_structs[i])
+
+        # Add correct cell tilt
+        direction = kmap[-1]
+
+        bulk = ref_bulk.copy() * (1, 1, N)
+
+        # Cell won't always be periodic, make sure we end up with something that is
+        mask, cell = self._get_kink_quad_mask(bulk, direction, layer_decimal_precision)
+
+        kink_quad = kink_quad[mask]
+        kink_quad.set_cell(cell, scale_atoms=False)
+
+        return kink_quad
+    
+    #TODO: Fix this
+    def build_minimal_kink_quadrupole(self, n_kink=1, layer_decimal_precision=3, invert_direction=False, smooth_width=None,
                                 *args, **kwargs):
         '''
         Construct a quadrupole structure providing an initial guess of the dislocation kink
@@ -4203,7 +4332,7 @@ def build_kink_quadrupole(self, z_reps=2, layer_decimal_precision=3, n_kink=1, i
         kink = self._smooth_kink_displacements(kink_struct1, kink_struct2, bulk, smooth_width, base_bulk.cell[:, :])
         return kink
 
-    def build_kink_quadrupole_network(self, kink_map=[0, 1], layer_decimal_precision=3, smooth_width=None,
+    def build_kink_quadrupole(self, kink_map=[0, 1], layer_decimal_precision=3, smooth_width=None,
                                 *args, **kwargs):
         '''
         Construct a quadrupole structure providing an initial guess of the dislocation kink
@@ -4212,65 +4341,28 @@ def build_kink_quadrupole_network(self, kink_map=[0, 1], layer_decimal_precision
 
         Parameters
         ----------
-        z_reps: int
-            Number of replications of self.unit_cell to use per kink
-            Should be >1
-        layer_tol: float
+        kink_map: iterable of ints or None
+            Map of the location of the dislocation core in units of the glide vector
+            Default (kink_map=None) is a kink map of [0, 1]
+            See examples for more details.
+        layer_decimal_precision: 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
+            (required in order to ensure complete periodicity)
+            Top layer is defined as any atom with z component >= max(atom_z) - 10**(-layer_tol)
+            in fractional coordinates
+        smooth_width: float
+            Size (in Ang) of the region for displacement smoothing at each kink site.
+            Larger smoothing width assumes a broader kink structure.
+            Default is 0.5 * self.unit_cell.cell[2, 2]
         *args, **kwargs
-            Fed to self.build_quadrupole() & self.build_glide_quadrupoles
+            Fed to self.build_quadrupole()
         
         '''
 
-        kmap = np.array(kink_map, dtype=int)
-        kmap -= np.min(kmap)
-        krange = np.max(kmap)
-
-        unique_kinks = np.sort(np.unique(kmap))
-
-        glide_structs = {}
+        ref_bulk, glide_structs = self.build_kink_quadrupole_glide_structs(kink_map, *args, **kwargs)
 
-        for kink_pos in unique_kinks:
-            off = np.array([kink_pos * self.glide_distance, 0 , 0])
-            ref_bulk, glide_structs[kink_pos] = self.build_quadrupole(*args, 
-                                                          left_offset=off, right_offset=off,
-                                                          **kwargs)
-        
-        kink_structs = [glide_structs[midx].copy() for midx in kmap]
-
-        # Smooth displacements across kink boundaries
-        if smooth_width is None:
-            smooth_width = 0.5 * self.unit_cell.cell[2, 2]
-
-        N = len(kink_structs)
-
-        smoothed_kink_structs = []
-        # Extra copy of last struct, so last one gets correct unsmoothed displacements
-        kink_structs.append(kink_structs[-1].copy())
-        for i in range(N):
-            smoothed_kink_structs.append(
-                self._smooth_kink_displacements(kink_structs[i], kink_structs[i+1], ref_bulk, smooth_width, ref_bulk.cell[:, :])
-            )
-
-        direction = kmap[-1]
-
-        bulk = ref_bulk.copy() * (1, 1, N)
-        kink = smoothed_kink_structs[0].copy()
-        for i in range(1, N):
-            kink = stack(kink, smoothed_kink_structs[i])
-
-        # Cell won't always be periodic, make sure we end up with something that is
-        mask, cell = self._get_kink_quad_mask(bulk, direction, layer_decimal_precision)
-
-        kink = kink[mask]
-        kink.set_cell(cell, scale_atoms=False)
-        return kink
+        return self.build_kink_quadrupole_from_glide_structs(kink_map, ref_bulk, glide_structs,
+                                                             smooth_width, layer_decimal_precision)
     
     def view_quad(self, system, *args, **kwargs):
         '''

From 6d97fe70159ef65d1a071b5495b1e223e5ff2d34 Mon Sep 17 00:00:00 2001
From: thomas-rocke <tom@rocke.co.uk>
Date: Sat, 17 Aug 2024 14:33:47 +0100
Subject: [PATCH 26/46] MAINT: Non-mutable default for kink_map

---
 matscipy/dislocation.py | 24 ++++++++++++------------
 1 file changed, 12 insertions(+), 12 deletions(-)

diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index c1b604dc..fa4e9d2e 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -3203,19 +3203,19 @@ def build_kink_glide_structs(self, kink_map=None, *args, **kwargs):
 
         return ref_bulk, glide_structs
     
-    def build_kink_from_glide_cyls(self, kink_map, ref_bulk, glide_structs, smooth_width=None):
+    def build_kink_from_glide_cyls(self, ref_bulk, glide_structs, kink_map=None, smooth_width=None):
         """
         Build a kink cylinder cell from the given kink map, using the structures contained in glide_structs
 
-        kink_map: iterable of ints or None
-            Map of the location of the dislocation core in units of the glide vector
-            Default (kink_map=None) is a kink map of [0, 1]
-            See examples for more details.
         ref_bulk: ase Atoms object
             Reference bulk structure, as returned by self.build_kink_glide_structs
         glide_structs: list of ase Atoms
             Glide structures e.g. those produced by self.build_kink_glide_structs.
             The kink structure is constructed using glide_structs[kink_map[i]] for the ith cell.
+        kink_map: iterable of ints or None
+            Map of the location of the dislocation core in units of the glide vector
+            Default (kink_map=None) is a kink map of [0, 1]
+            See examples for more details.
         smooth_width: float
             Size (in Ang) of the region for displacement smoothing at each kink site.
             Larger smoothing width assumes a broader kink structure.
@@ -3711,7 +3711,7 @@ def build_cylinder(self,
         if return_fix_mask:
             out.append(fix_mask)
 
-        return out
+        return tuple(out)
 
 class CubicCrystalDislocationQuadrupole(CubicCrystalDissociatedDislocation):
     burgers_dimensionless = np.zeros(3)
@@ -4206,19 +4206,19 @@ def build_kink_quadrupole_glide_structs(self, kink_map=None, *args, **kwargs):
                                                           **kwargs)
         return ref_bulk, glide_structs
 
-    def build_kink_quadrupole_from_glide_structs(self, kink_map, ref_bulk, glide_structs, smooth_width=None, layer_decimal_precision=3):
+    def build_kink_quadrupole_from_glide_structs(self, ref_bulk, glide_structs, kink_map=None, smooth_width=None, layer_decimal_precision=3):
         """
         Build a kink cylinder cell from the given kink map, using the structures contained in glide_structs
 
-        kink_map: iterable of ints or None
-            Map of the location of the dislocation core in units of the glide vector
-            Default (kink_map=None) is a kink map of [0, 1]
-            See examples for more details.
         ref_bulk: ase Atoms object
             Reference bulk structure, as returned by self.build_kink_quadrupole_glide_structs
         glide_structs: list of ase Atoms
             Glide structures e.g. those produced by self.build_kink_quadrupole_glide_structs.
             The kink structure is constructed using glide_structs[kink_map[i]] for the ith cell.
+        kink_map: iterable of ints or None
+            Map of the location of the dislocation core in units of the glide vector
+            Default (kink_map=None) is a kink map of [0, 1]
+            See examples for more details.
         smooth_width: float
             Size (in Ang) of the region for displacement smoothing at each kink site.
             Larger smoothing width assumes a broader kink structure.
@@ -4332,7 +4332,7 @@ def build_minimal_kink_quadrupole(self, n_kink=1, layer_decimal_precision=3, inv
         kink = self._smooth_kink_displacements(kink_struct1, kink_struct2, bulk, smooth_width, base_bulk.cell[:, :])
         return kink
 
-    def build_kink_quadrupole(self, kink_map=[0, 1], layer_decimal_precision=3, smooth_width=None,
+    def build_kink_quadrupole(self, kink_map=None, layer_decimal_precision=3, smooth_width=None,
                                 *args, **kwargs):
         '''
         Construct a quadrupole structure providing an initial guess of the dislocation kink

From 2f9c95fb6abf139ba50fc44eb72631e1afe16fa3 Mon Sep 17 00:00:00 2001
From: thomas-rocke <tom@rocke.co.uk>
Date: Tue, 20 Aug 2024 15:29:27 +0100
Subject: [PATCH 27/46] BUG: Allow CCDD build_cylinder to specify both core
 positions

---
 matscipy/dislocation.py | 21 ++++++++++++++++-----
 1 file changed, 16 insertions(+), 5 deletions(-)

diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index fa4e9d2e..ce3b46aa 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -24,6 +24,7 @@
 
 """Tools for studying structure and movement of dislocations."""
 
+from ctypes import ArgumentError
 import numpy as np
 
 from abc import ABCMeta
@@ -3681,11 +3682,21 @@ def build_cylinder(self,
         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
-        ])
+        if np.atleast_2d(core_position).shape[0] == 1:
+            # Core positions for only one core specified
+            # Use partial distance argument to find second core
+            core_positions = np.array([
+                core_position + self.unit_cell_core_position,
+                core_position + self.unit_cell_core_position + partial_distance_Angstrom
+            ])
+        else:
+            # Two core positions specified
+            if partial_distance != 0:
+                # Strange input specifying both core positions, plus a partial distance
+                raise ArgumentError("Core positions for both cores and a partial distance were both specified.\n" + 
+                                    "Either specify both core positions directly with the core_position argument\n" + 
+                                    "or specify the position of the left core, and a partial distance")
+            core_positions = core_position + self.unit_cell_core_position
 
         bulk, disloc, core_positions, cyl_mask, fix_mask = self._build_bulk_cyl(radius, core_positions, fix_width, extension,
                                                             fixed_points, self_consistent, method, verbose, cyl_mask=cyl_mask, **kwargs)

From d72edc1668ef2d7baa5bcebebe35235685ddf060 Mon Sep 17 00:00:00 2001
From: thomas-rocke <tom@rocke.co.uk>
Date: Tue, 20 Aug 2024 15:43:53 +0100
Subject: [PATCH 28/46] BUG: Fix arg order in build_kink_cylinder

---
 matscipy/dislocation.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index ce3b46aa..21ff9821 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -3174,7 +3174,6 @@ def build_kink_glide_structs(self, kink_map=None, *args, **kwargs):
         ])
 
         for kink_pos in unique_kinks:
-            print(kink_pos, krange)
             ref_bulk, glide_structs[kink_pos] = self.build_cylinder(*args, 
                                                           fixed_points=fixed_points, 
                                                           core_position=np.array([kink_pos * self.glide_distance, 0 , 0]),
@@ -3280,7 +3279,7 @@ def build_kink_cylinder(self, kink_map=None, smooth_width=None, *args, **kwargs)
         """
         ref_bulk, glide_structs = self.build_kink_glide_structs(kink_map, *args, **kwargs)
 
-        return self.build_kink_from_glide_cyls(kink_map, ref_bulk, glide_structs, smooth_width)
+        return self.build_kink_from_glide_cyls(ref_bulk, glide_structs, kink_map, smooth_width)
     
     @staticmethod
     def view_cyl(system, scale=0.5, CNA_color=True, add_bonds=False,

From 9203e6a7d8fd479af83bffadddc6dbddcc8bff4c Mon Sep 17 00:00:00 2001
From: thomas-rocke <tom@rocke.co.uk>
Date: Tue, 20 Aug 2024 16:18:06 +0100
Subject: [PATCH 29/46] MAINT: Re-add exception for BCC Edge 111bar quad kink
 construction

---
 matscipy/dislocation.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index 21ff9821..14e86102 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -4124,7 +4124,7 @@ def _get_kink_quad_mask(self, ref_bulk, direction, precision=3):
         '''
 
         if self.left_dislocation.__class__ == BCCEdge111barDislocation:
-            #raise RuntimeError("Kink Quadrupoles do not currently work for BCCEdge111barDislocation")
+            raise RuntimeError("Kink Quadrupoles do not currently work for BCCEdge111barDislocation")
             pass
 
         bulk = ref_bulk.copy()
@@ -4371,7 +4371,7 @@ def build_kink_quadrupole(self, kink_map=None, layer_decimal_precision=3, smooth
 
         ref_bulk, glide_structs = self.build_kink_quadrupole_glide_structs(kink_map, *args, **kwargs)
 
-        return self.build_kink_quadrupole_from_glide_structs(kink_map, ref_bulk, glide_structs,
+        return self.build_kink_quadrupole_from_glide_structs(ref_bulk, glide_structs, kink_map,
                                                              smooth_width, layer_decimal_precision)
     
     def view_quad(self, system, *args, **kwargs):

From 368ec850a9e94e3671bab5baa6f85f2675981918 Mon Sep 17 00:00:00 2001
From: thomas-rocke <tom@rocke.co.uk>
Date: Tue, 20 Aug 2024 16:23:04 +0100
Subject: [PATCH 30/46] MAINT: Fix broken test (Caused by finding min repr of
 rotated cell)

---
 tests/test_gamma_surface.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/tests/test_gamma_surface.py b/tests/test_gamma_surface.py
index 18b363a1..21953c72 100644
--- a/tests/test_gamma_surface.py
+++ b/tests/test_gamma_surface.py
@@ -27,10 +27,11 @@ def test_gamma_surface(self):
 
     def test_stacking_fault(self):
         surface = StackingFault(self.at0, np.array([1, 1, 0]), np.array([-1, 1, 0]))
-        surface.generate_images(9, z_reps=1, path_ylims=[0, 0.5])
+        surface.generate_images(9, z_reps=1, path_ylims=[0, 1])
         surface.relax_images(self.model, self.fmax)
         Es = surface.get_energy_densities(self.model)[0, :]
-        assert np.allclose(np.max(Es), 0.27185793519964985)
+        print(np.max(Es))
+        assert np.allclose(np.max(Es), 0.2628608260182023)
 
     def test_disloc_stacking_fault(self):
         from matscipy.dislocation import DiamondGlideScrew

From 00db06fb3eb57eb4731b64d693601ae429b8fc9c Mon Sep 17 00:00:00 2001
From: Thomas Rocke <45517493+thomas-rocke@users.noreply.github.com>
Date: Thu, 22 Aug 2024 15:31:57 +0100
Subject: [PATCH 31/46] DOC: Glide Screw -> 60deg for mobility docs

---
 docs/applications/disloc_mobility.ipynb | 65 +++++--------------------
 1 file changed, 12 insertions(+), 53 deletions(-)

diff --git a/docs/applications/disloc_mobility.ipynb b/docs/applications/disloc_mobility.ipynb
index 98890c76..c28fe8ca 100644
--- a/docs/applications/disloc_mobility.ipynb
+++ b/docs/applications/disloc_mobility.ipynb
@@ -21,52 +21,11 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "application/vnd.jupyter.widget-view+json": {
-       "model_id": "273fa26b287d4c61a24ff07d20735b7f",
-       "version_major": 2,
-       "version_minor": 0
-      },
-      "text/plain": []
-     },
-     "metadata": {},
-     "output_type": "display_data"
-    },
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "/home/europa/.local/lib/python3.10/site-packages/matscipy/dislocation.py:484: FutureWarning: Import StrainFilter from ase.filters\n",
-      "  sf = StrainFilter(unit_cell)\n"
-     ]
-    },
-    {
-     "ename": "ImportError",
-     "evalue": "libQt6Core.so.6: cannot open shared object file: No such file or directory",
-     "output_type": "error",
-     "traceback": [
-      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
-      "\u001b[0;31mImportError\u001b[0m                               Traceback (most recent call last)",
-      "Cell \u001b[0;32mIn[1], line 14\u001b[0m\n\u001b[1;32m     11\u001b[0m Si_screw \u001b[38;5;241m=\u001b[39m DiamondGlideScrew(alat, C11, C12, C44, symbol\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mSi\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m     13\u001b[0m Si_screw_bulk, Si_screw_dislo \u001b[38;5;241m=\u001b[39m Si_screw\u001b[38;5;241m.\u001b[39mbuild_cylinder(radius\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m20\u001b[39m)\n\u001b[0;32m---> 14\u001b[0m \u001b[43mSi_screw\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mview_cyl\u001b[49m\u001b[43m(\u001b[49m\u001b[43mSi_screw_dislo\u001b[49m\u001b[43m)\u001b[49m\n",
-      "File \u001b[0;32m~/.local/lib/python3.10/site-packages/matscipy/dislocation.py:3306\u001b[0m, in \u001b[0;36mCubicCrystalDislocation.view_cyl\u001b[0;34m(system, scale, CNA_color, add_bonds, line_color, disloc_names, hide_arrows)\u001b[0m\n\u001b[1;32m   3303\u001b[0m \u001b[38;5;66;03m# Check if system contains a diamond dislocation (e.g. DiamondGlideScrew, Diamond90degreePartial)\u001b[39;00m\n\u001b[1;32m   3304\u001b[0m diamond_structure \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mDiamond\u001b[39m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;129;01min\u001b[39;00m system\u001b[38;5;241m.\u001b[39minfo[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mdislocation_classes\u001b[39m\u001b[38;5;124m\"\u001b[39m][\u001b[38;5;241m0\u001b[39m]\n\u001b[0;32m-> 3306\u001b[0m atom_labels, structure_names, colors \u001b[38;5;241m=\u001b[39m \u001b[43mget_structure_types\u001b[49m\u001b[43m(\u001b[49m\u001b[43msystem\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\n\u001b[1;32m   3307\u001b[0m \u001b[43m                                                        \u001b[49m\u001b[43mdiamond_structure\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdiamond_structure\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m   3308\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m disloc_names \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m   3309\u001b[0m     disloc_names \u001b[38;5;241m=\u001b[39m system\u001b[38;5;241m.\u001b[39minfo[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mdislocation_types\u001b[39m\u001b[38;5;124m\"\u001b[39m]\n",
-      "File \u001b[0;32m~/.local/lib/python3.10/site-packages/matscipy/utils.py:342\u001b[0m, in \u001b[0;36mget_structure_types\u001b[0;34m(structure, diamond_structure)\u001b[0m\n\u001b[1;32m    331\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mget_structure_types\u001b[39m(structure, diamond_structure\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mFalse\u001b[39;00m):\n\u001b[1;32m    332\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"Get the results of Common Neighbor Analysis and \u001b[39;00m\n\u001b[1;32m    333\u001b[0m \u001b[38;5;124;03m        Identify Diamond modifiers from Ovito\u001b[39;00m\n\u001b[1;32m    334\u001b[0m \u001b[38;5;124;03m        (Requires Ovito python module)\u001b[39;00m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    340\u001b[0m \u001b[38;5;124;03m        colors (list of strings): colors of the structure types in hex format\u001b[39;00m\n\u001b[1;32m    341\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[0;32m--> 342\u001b[0m     \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01movito\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mio\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mase\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m ase_to_ovito\n\u001b[1;32m    343\u001b[0m     \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01movito\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mmodifiers\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m CommonNeighborAnalysisModifier, IdentifyDiamondModifier\n\u001b[1;32m    344\u001b[0m     \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01movito\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mpipeline\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m StaticSource, Pipeline\n",
-      "File \u001b[0;32m~/.local/lib/python3.10/site-packages/ovito/__init__.py:21\u001b[0m\n\u001b[1;32m     18\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01movito\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mnonpublic\u001b[39;00m\n\u001b[1;32m     20\u001b[0m \u001b[38;5;66;03m# Load all extension modules.\u001b[39;00m\n\u001b[0;32m---> 21\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01movito\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01m_extensions\u001b[39;00m\n\u001b[1;32m     23\u001b[0m \u001b[38;5;66;03m# The public symbols of this root module:\u001b[39;00m\n\u001b[1;32m     24\u001b[0m __all__ \u001b[38;5;241m=\u001b[39m [\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mversion\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mversion_string\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mscene\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mScene\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124menable_logging\u001b[39m\u001b[38;5;124m'\u001b[39m]\n",
-      "File \u001b[0;32m~/.local/lib/python3.10/site-packages/ovito/_extensions/__init__.py:6\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[38;5;66;03m# Load all plugin extension modules under the _extensions package.\u001b[39;00m\n\u001b[1;32m      5\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m modinfo \u001b[38;5;129;01min\u001b[39;00m pkgutil\u001b[38;5;241m.\u001b[39mwalk_packages(__path__, \u001b[38;5;18m__name__\u001b[39m \u001b[38;5;241m+\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124m.\u001b[39m\u001b[38;5;124m'\u001b[39m):\n\u001b[0;32m----> 6\u001b[0m     \u001b[43mimportlib\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mimport_module\u001b[49m\u001b[43m(\u001b[49m\u001b[43mmodinfo\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mname\u001b[49m\u001b[43m)\u001b[49m\n",
-      "File \u001b[0;32m/usr/lib/python3.10/importlib/__init__.py:126\u001b[0m, in \u001b[0;36mimport_module\u001b[0;34m(name, package)\u001b[0m\n\u001b[1;32m    124\u001b[0m             \u001b[38;5;28;01mbreak\u001b[39;00m\n\u001b[1;32m    125\u001b[0m         level \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m \u001b[38;5;241m1\u001b[39m\n\u001b[0;32m--> 126\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43m_bootstrap\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_gcd_import\u001b[49m\u001b[43m(\u001b[49m\u001b[43mname\u001b[49m\u001b[43m[\u001b[49m\u001b[43mlevel\u001b[49m\u001b[43m:\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mpackage\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mlevel\u001b[49m\u001b[43m)\u001b[49m\n",
-      "File \u001b[0;32m~/.local/lib/python3.10/site-packages/ovito/_extensions/anari.py:2\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[38;5;66;03m# Load dependencies.\u001b[39;00m\n\u001b[0;32m----> 2\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01movito\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01m_extensions\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mpyscript\u001b[39;00m\n\u001b[1;32m      4\u001b[0m \u001b[38;5;66;03m# Load the C extension module.\u001b[39;00m\n\u001b[1;32m      5\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01movito\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mplugins\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mAnariRendererPython\u001b[39;00m\n",
-      "File \u001b[0;32m~/.local/lib/python3.10/site-packages/ovito/_extensions/pyscript.py:2\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[38;5;66;03m# Load the C extension module.\u001b[39;00m\n\u001b[0;32m----> 2\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01movito\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mplugins\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mPyScript\u001b[39;00m\n\u001b[1;32m      4\u001b[0m \u001b[38;5;66;03m# Load class add-ons.\u001b[39;00m\n\u001b[1;32m      5\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01movito\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01m_scene_class\u001b[39;00m\n",
-      "File \u001b[0;32m~/.local/lib/python3.10/site-packages/ovito/plugins/__init__.py:36\u001b[0m\n\u001b[1;32m     26\u001b[0m     warnings\u001b[38;5;241m.\u001b[39mwarn(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mIncompatible version of the Qt cross-platform framework detected!\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124mThis version of the OVITO Python module is based on Qt6 (loaded via the PySide6 bindings module), \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m     27\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mbut bindings for old Qt5 are already loaded at this point (through PyQt5 or PySide2 imports preceding the import of OVITO). To avoid library version conflicts, please make sure the rest of \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m     28\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124myour application uses Qt6 too instead of Qt5. \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m     31\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mIn addition, it may help to set the environment variable QT_API=pyside6 to force third-party packages (e.g. matplotlib) to load Qt6 instead of Qt5. \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m     32\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mIf you have any questions, please contact support@ovito.org.\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m     34\u001b[0m \u001b[38;5;66;03m# Load all the Qt bindings first before OVITO's own C++ modules get loaded.\u001b[39;00m\n\u001b[1;32m     35\u001b[0m \u001b[38;5;66;03m# This ensures that the right Qt shared libraries needed by OVITO are already loaded into the process when running in a system Python interpreter.\u001b[39;00m\n\u001b[0;32m---> 36\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01movito\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mqt_compat\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m QtCore\n\u001b[1;32m     37\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01movito\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mqt_compat\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m QtGui\n\u001b[1;32m     38\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01movito\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mqt_compat\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m QtWidgets\n",
-      "File \u001b[0;32m~/.local/lib/python3.10/site-packages/ovito/qt_compat.py:11\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[38;5;66;03m# Import Shiboken module first to preload shiboken6.dll, which is needed by the Qt DLLs.\u001b[39;00m\n\u001b[1;32m      9\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mshiboken6\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mshiboken\u001b[39;00m \n\u001b[0;32m---> 11\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mPySide6\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m QtCore, QtGui, QtWidgets, QtNetwork, QtXml\n\u001b[1;32m     13\u001b[0m \u001b[38;5;66;03m# The QtOpenGLWidgets module is not included in the embedded PySide6 installation of OVITO Pro 3.7,\u001b[39;00m\n\u001b[1;32m     14\u001b[0m \u001b[38;5;66;03m# but it should be loaded when using the PyPI version of PySide6.\u001b[39;00m\n\u001b[1;32m     15\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n",
-      "\u001b[0;31mImportError\u001b[0m: libQt6Core.so.6: cannot open shared object file: No such file or directory"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
-    "from matscipy.dislocation import get_elastic_constants, DiamondGlideScrew\n",
+    "from matscipy.dislocation import get_elastic_constants, DiamondGlide60Degree\n",
     "# the calculator to provide forces and energies from the potential\n",
     "from matscipy.calculators.manybody.explicit_forms.stillinger_weber import StillingerWeber,\\\n",
     "                                                                Holland_Marder_PRL_80_746_Si\n",
@@ -76,10 +35,10 @@
     "# the function accepts any ASE type of calculator\n",
     "alat, C11, C12, C44 = get_elastic_constants(calculator=calc, symbol=\"Si\", verbose=False)\n",
     "\n",
-    "Si_screw = DiamondGlideScrew(alat, C11, C12, C44, symbol=\"Si\")\n",
+    "Si_disloc = DiamondGlide60Degree(alat, C11, C12, C44, symbol=\"Si\")\n",
     "\n",
-    "Si_screw_bulk, Si_screw_dislo = Si_screw.build_cylinder(radius=20)\n",
-    "Si_screw.view_cyl(Si_screw_dislo)\n"
+    "Si_disloc_bulk, Si_disloc_dislo = Si_disloc.build_cylinder(radius=20)\n",
+    "Si_disloc.view_cyl(Si_disloc_dislo)\n"
    ]
   },
   {
@@ -97,7 +56,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "bulk, glide_ini, glide_fin = Si_screw.build_glide_configurations(radius=20)\n",
+    "bulk, glide_ini, glide_fin = Si_disloc.build_glide_configurations(radius=20)\n",
     "\n",
     "from ase.neb import NEB\n",
     "\n",
@@ -171,7 +130,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "bulk, glide_ini, glide_fin = Si_screw.build_glide_configurations(radius=20, partial_distance=5)\n",
+    "bulk, glide_ini, glide_fin = Si_disloc.build_glide_configurations(radius=20, partial_distance=5)\n",
     "\n",
     "nims = 5\n",
     "\n",
@@ -200,7 +159,7 @@
    "source": [
     "from matscipy.dislocation import Quadrupole\n",
     "\n",
-    "Si_quad = Quadrupole(DiamondGlideScrew, alat, C11, C12, C44, symbol=\"Si\")\n",
+    "Si_quad = Quadrupole(DiamondGlide60Degree, alat, C11, C12, C44, symbol=\"Si\")\n",
     "\n",
     "num_images = 5\n",
     "\n",
@@ -266,12 +225,12 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "Si_kink = Si_screw.build_kink_cyl(\n",
+    "Si_kink = Si_disloc.build_kink_cyl(\n",
     "                    radius=20,\n",
     "                    kink_map= [0] * 5 + [1] * 5\n",
     "                )\n",
     "\n",
-    "Si_screw.view_cyl(Si_kink)"
+    "Si_disloc.view_cyl(Si_kink)"
    ]
   },
   {
@@ -319,7 +278,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.14"
+   "version": "3.10.12"
   }
  },
  "nbformat": 4,

From 6ba85a1ac0522e03b3a57d50f47dc8b075148823 Mon Sep 17 00:00:00 2001
From: Thomas Rocke <45517493+thomas-rocke@users.noreply.github.com>
Date: Tue, 3 Sep 2024 18:51:17 +0100
Subject: [PATCH 32/46] ENH: Disloc DXA for view_cyl

---
 docs/applications/disloc_mobility.ipynb |   6 +-
 matscipy/dislocation.py                 | 380 ++++++++++++++++++++----
 2 files changed, 317 insertions(+), 69 deletions(-)

diff --git a/docs/applications/disloc_mobility.ipynb b/docs/applications/disloc_mobility.ipynb
index c28fe8ca..b7753051 100644
--- a/docs/applications/disloc_mobility.ipynb
+++ b/docs/applications/disloc_mobility.ipynb
@@ -38,7 +38,7 @@
     "Si_disloc = DiamondGlide60Degree(alat, C11, C12, C44, symbol=\"Si\")\n",
     "\n",
     "Si_disloc_bulk, Si_disloc_dislo = Si_disloc.build_cylinder(radius=20)\n",
-    "Si_disloc.view_cyl(Si_disloc_dislo)\n"
+    "Si_disloc.view_cyl(Si_disloc_dislo, mode=\"dxa\")\n"
    ]
   },
   {
@@ -184,7 +184,7 @@
     "                                            glide_left=False, # Prevent left dislocation glide\n",
     "                                            glide_right=True, # Allow right dislocation glide\n",
     "                                            partial_distance=2, # Dissociated Glide\n",
-    "                                            glide_separation=6,\n",
+    "                                            glide_separation=8,\n",
     "                                            verbose=False)\n",
     "\n",
     "animate_glide(glide_quads)"
@@ -225,7 +225,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "Si_kink = Si_disloc.build_kink_cyl(\n",
+    "Si_kink = Si_disloc.build_kink_cylinder(\n",
     "                    radius=20,\n",
     "                    kink_map= [0] * 5 + [1] * 5\n",
     "                )\n",
diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index 14e86102..2a61e4f4 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -2118,7 +2118,6 @@ def ovito_dxa_straight_dislo_info(disloc, structure="BCC", replicate_z=3):
 
     return results
 
-
 def get_centering_mask(atoms, radius,
                        core_position=[0., 0., 0.],
                        extension=[0., 0., 0.],):
@@ -3281,9 +3280,309 @@ def build_kink_cylinder(self, kink_map=None, smooth_width=None, *args, **kwargs)
 
         return self.build_kink_from_glide_cyls(ref_bulk, glide_structs, kink_map, smooth_width)
     
-    @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):
+    def ovito_dxa(self, disloc):
+        """
+        General interface for ovito dxa analysis.
+        dxa requires 3b thick cells at minimum, so uses ovito replicate if given structure is insufficient
+
+        Parameters
+        ----------
+        disloc: ase.Atoms
+            Atoms object containing the atomic configuration to analyse
+
+        Returns
+        -------
+        data: ovito DataCollection object. 
+            See https://www.ovito.org/docs/current/python/modules/ovito_data.html#ovito.data.DataCollection
+            data.dislocations.segments is a list of all discovered dislocation segments
+
+        """
+        from ovito.io.ase import ase_to_ovito
+        from ovito.modifiers import ReplicateModifier, DislocationAnalysisModifier
+        from ovito.pipeline import StaticSource, Pipeline
+
+        dxa_disloc = disloc.copy()
+        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}
+        
+        data = ase_to_ovito(dxa_disloc)
+        pipeline = Pipeline(source=StaticSource(data=data))
+        
+        # Ensure disloc structure is at least 3b thick
+        cell_scale = int(disloc.cell[2, 2] // self.unit_cell.cell[2, 2])
+        
+        if cell_scale < 3:
+            replicate_z = int(np.ceil(3 / cell_scale))
+        else:
+            replicate_z = 1
+        
+        pipeline.modifiers.append(ReplicateModifier(num_z=replicate_z))
+        
+        # Setup DXA pipeline
+        dxa = DislocationAnalysisModifier(
+                input_crystal_structure=input_crystal_structures[self.crystalstructure.lower()])
+        pipeline.modifiers.append(dxa)
+        
+        data = pipeline.compute()
+
+        return data
+
+
+    def _plot_CLE_disloc_line(self, view, disloc, z_length, disloc_names=None, 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: ase Atoms
+            Dislocation structure
+        disloc_names: str
+            Custom names for the dislocations 
+        color: list
+            [R, G, B] colour of the dislocation line
+        '''
+
+        # 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 disloc.info.keys():
+                warnings.warn(f"{expected_key} not found in system.info, skipping plot of dislocation line")
+        
+        
+        # Validate line_color arg
+        if type(color) in [list, np.array, tuple]:
+            if type(color[0]) in [list, np.array, tuple]:
+                # Given an RGB value per dislocation
+                colours = color
+            else:
+                # Assume we're given a single RBG value for all dislocs
+                colours = [color] * len(disloc.info["dislocation_types"])
+        
+        if disloc_names is None:
+            disloc_names = [name + " dislocation line" for name in disloc.info["dislocation_types"]]
+
+        disloc_pos = disloc.info["core_positions"]
+
+        if type(disloc_names) in (list, tuple, np.array):
+            if len(disloc_names) != len(disloc_pos):
+                raise RuntimeError("Number of dislocation positions is not equal to number of dislocation types/names")
+
+        for i in range(len(disloc_pos)):
+            view.shape.add_cylinder((disloc_pos[i][0], disloc_pos[i][1], -2.0), 
+                                    (disloc_pos[i][0], disloc_pos[i][1], z_length - 0.5),
+                                    colours,
+                                    [0.3],
+                                    disloc_names[i])
+            
+            view.shape.add_cone((disloc_pos[i][0], disloc_pos[i][1], -2.0), 
+                                (disloc_pos[i][0], disloc_pos[i][1], 0.5),
+                                colours,
+                                [0.3],
+                                disloc_names[i])
+            
+            view.shape.add_cone((disloc_pos[i][0], disloc_pos[i][1], z_length - 0.5), 
+                                (disloc_pos[i][0], disloc_pos[i][1], z_length + 1.0),
+                                colours,
+                                [0.55],
+                                disloc_names[i])
+        return view
+
+    def _plot_DXA_disloc_line(self, view, disloc, disloc_names=None, color=None):
+        def inf_line_transform(p, z_max):
+            """
+            Find the dislocation line spanning 0<=p[:, 2]<=z_max from points which may lie outside the cell
+
+            Parameters
+            ----------
+            p: np.array
+                Points to transform
+            z_max: float
+                cell[2, 2] component.
+
+            Returns
+            -------
+            points: np.array
+                Transformed points
+            
+            """
+            points = []
+            
+            # Start by translating by z_max (periodicity) such that there are points below the cell
+            n = np.ceil(np.min(p[:, 2]) / z_max)
+            p_off = n * z_max
+            p[:, 2] -= p_off
+
+            # Now start from the point with negative z closest to z=0
+            z_tmp = p[:, 2].copy()
+            z_tmp[z_tmp > 0] = 100
+            argstart = np.argmin(np.abs(z_tmp))
+
+            # Next point will be above z=0
+            # Generate the line segment from z=0 to this point
+            p1 = p[argstart, :]
+            p2 = p[argstart+1, :]
+            dp = p2 - p1
+            dz = dp[2]
+
+            p_intercept = p1 + dp * (0 - p1[2]) / dz # z=0 intercept point
+            points.append(p_intercept.copy())
+
+            for i in range(argstart+1, p.shape[0]):
+                p2 = p[i, :]
+                p1 = p[i-1, :]
+
+                if p2[2] > z_max:
+                    # Current point is off the end of the cell
+                    # Generate the line segment between p1 and this intercept
+                    dp = p2 - p1
+                    dz = dp[2]
+                
+                    p_intercept = p1 + dp * (z_max - p1[2]) / dz # z=z_max intercept point
+                    points.append(p_intercept)
+
+                    # Can stop now, as we have points spanning the full cell
+                    return np.array(points)
+                    
+                points.append(p2)
+            
+            if p[-1, 2] < z_max:
+                # Missing the line segment(s) wrapping over the periodic boundary
+                for i in range(argstart+1):
+                    p1 = p[i, :].copy()
+                    p2 = p[i+1, :].copy()
+                    p1[2] += z_max
+                    p2[2] += z_max
+
+                    points.append(p1)
+                    
+                    if p2[2] > z_max:
+                        # Again crossing the boundary
+                        # Generate line segment up to z_max
+                        dp = p2 - p1
+                        dz = dp[2]
+                    
+                        p_intercept = p1 + dp * (z_max - p1[2]) / dz # z=z_max intercept point
+                        points.append(p_intercept)
+            
+                        # Can stop now, as we have points spanning the full cell
+                        return np.array(points)
+                return np.array(points)
+        from fractions import Fraction
+        
+        # Ovito colours for dislocation lines (as of v 3.10.1) 
+        disloc_colours = {
+            "diamond" : {
+                "default" : [230, 51, 51], # Red
+                "1/2 1/2 0" : [51, 51, 255], # Dark Blue; Screw, 60 degree
+                "1/3 1/6 1/6" : [0, 255, 0], # Green; 30 & 90 degree Partials
+                "1/6 1/6 0" : [255, 0, 255], # Pink
+                "1/3 1/3 1/3" : [0, 255, 255] # Light Blue
+            }
+        }
+
+        z_max = disloc.cell[2, 2]
+        
+        data = self.ovito_dxa(disloc)
+
+        nseg = len(data.dislocations.segments)
+
+        if disloc_names is not None:
+            if len(disloc_names) != nseg:
+                warnings.warn(f"{len(disloc_names)} names specified, but found {nseg} dislocation segments. Skipping plot of dislocation line")
+                return view
+        
+        if color is not None:
+            if len(color) != nseg and (len(color) != 3 or (type(color[0]) not in (int, float))):
+                # Not enough colours specified, and colour is not a universal RGB
+                warnings.warn(f"{len(disloc_names)} names specified, but found {nseg} dislocation segments. Skipping plot of dislocation line")
+                return view
+            elif (len(color) == 3 and (type(color[0]) in (int, float))):
+                # Universal color for all segments
+                color = [color for i in range(nseg)]
+
+        for idx, segment in enumerate(data.dislocations.segments):
+            b = segment.true_burgers_vector
+            points = segment.points
+
+            if segment.is_infinite_line:
+                # Infinite dislocation line
+                # Points are likely to leave the cell
+                # Especially if we had to replicate
+                points = inf_line_transform(points, z_max)
+    
+                # Make sure points are unique to 2dp
+                _, idxs = np.unique(np.round(points.copy(), 2), return_index=True, axis=0)
+                points = points[idxs, :]
+                idxs = np.argsort(points[:, 2])
+                points = points[idxs, :]
+    
+            b_sorted = sorted(np.abs(b))[::-1]
+    
+            # Use burgers vector to get a key for the dislocation line colour
+            # Sorting and abs remove the rotational dependance, reducing the number of possible keys
+            colour_key = " ".join([str(Fraction(elem).limit_denominator(10)) for elem in b_sorted])
+
+            if color is None:
+                if colour_key in disloc_colours[self.crystalstructure.lower()]:
+                    colour = disloc_colours[self.crystalstructure.lower()][colour_key]
+                else:
+                    colour = disloc_colours[self.crystalstructure.lower()]["default"]
+                colour = np.array(colour) / 255
+            else:
+                colour = color[idx]
+
+            if disloc_names is None:
+                seg_name = f"b=[" + " ".join([str(Fraction(elem).limit_denominator(10)) for elem in b]) + "]"
+            else:
+                seg_name = disloc_names[idx]
+
+            for i in range(points.shape[0] - 1):
+                P1 = points[i, :]
+                P2 = points[i+1, :]
+            
+                view.shape.add_cylinder(P1, P2, colour, 0.3, seg_name)
+
+            # Add a cone on the end of the dislocation line, based on the burgers vector
+            # Helps to visualise dislocation quadrupoles
+
+            if segment.spatial_burgers_vector @ np.array([1, 1, 1]) >= 0:
+                # "Positive" burgers vector, arrow on top
+                p1 = points[-2, :]
+                p2 = points[-1, :]
+            else:
+                # "Negative" burgers vector, arrow on bottom
+                p1 = points[1, :]
+                p2 = points[0, :]
+
+            if not ((p2[2] - 0) **2 < 1e-3 or (p2[2] - z_max)**2 < 1e-3):
+                # p2 is not on either border, therefore the current segment shouldn't have an arrow drawn
+                continue
+
+            startpoint = p2
+
+            extension = 3.0
+            cone_length = 1.0
+            
+            dp = p2 - p1
+            dz = dp[2]
+        
+            endpoint = p2 + dp * (extension) / np.abs(dz)
+            coneendpoint = p2 + dp * (extension + cone_length) / np.abs(dz)
+
+            view.shape.add_cylinder(startpoint, endpoint, colour, 0.3, seg_name)
+            view.shape.add_cone(endpoint, coneendpoint, colour, 0.6, seg_name)
+        return view
+
+    def view_cyl(self, system, scale=0.5, CNA_color=True, add_bonds=False,
+                     line_color=[0, 1, 0], disloc_names=None, hide_arrows=False,
+                     mode="dxa"):
         '''
         NGLview-based visualisation tool for structures generated from the dislocation class
 
@@ -3305,41 +3604,12 @@ def view_cyl(system, scale=0.5, CNA_color=True, add_bonds=False,
             (see structure.info["dislocation_types"] for defaults)
         hide_arrows: bool
             Hide arrows for placed on the dislocation lines
+        mode: str
+            Mode for plotting dislocation line
+            mode="dxa" uses ovito dxa to plot an exact line
+            mode="cle" uses the cle solution used to generate the structure (if the dislocation line is straight)
         '''
 
-        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
@@ -3366,36 +3636,12 @@ def _plot_disloc_line(view, disloc_pos, disloc_name, z_length, color=[0, 1, 0]):
                             }
                         });
                         """
-        
-        
-        # 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]
+        diamond_structure = "diamond" in self.crystalstructure.lower()
 
         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])
         
@@ -3422,12 +3668,14 @@ def _plot_disloc_line(view, disloc_pos, disloc_name, z_length, color=[0, 1, 0]):
                         
         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])
+            if mode.lower() == "cle":
+                # Plot CLE dislocation lines
+                z_length = system.cell[2, 2]
+                self._plot_CLE_disloc_line(view, system, z_length, disloc_names, line_color)
+            elif mode.lower() == "dxa":
+                self._plot_DXA_disloc_line(view, system, disloc_names, line_color)
 
         view.camera = 'orthographic'
         view.parameters = {"clipDist": 0}

From e9874e0b343e6033d6da1540507b6193a0698d18 Mon Sep 17 00:00:00 2001
From: thomas-rocke <tom@rocke.co.uk>
Date: Tue, 3 Sep 2024 22:40:13 +0100
Subject: [PATCH 33/46] ENH: Fix build_minimal_kink_quadrupole

---
 matscipy/dislocation.py | 149 ++++++++++++++++++++--------------------
 1 file changed, 73 insertions(+), 76 deletions(-)

diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index 2a61e4f4..939f562c 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -212,7 +212,7 @@ def make_screw_cyl(alat, C11, C12, C44,
 def make_edge_cyl(alat, C11, C12, C44,
                   cylinder_r=10, cutoff=5.5,
                   symbol='W'):
-    '''
+    """
     makes edge dislocation using atomman library
 
     cylinder_r - radius of cylinder of unconstrained atoms around the
@@ -223,7 +223,7 @@ def make_edge_cyl(alat, C11, C12, C44,
     symbol : string
         Symbol of the element to pass to ase.lattice.cubic.SimpleCubicFactory
         default is "W" for tungsten
-    '''
+    """
     from atomman import ElasticConstants
     from atomman.defect import Stroh
     # Create a Stroh object with junk data
@@ -2333,9 +2333,9 @@ def deformation_gradient(self, bulk_positions, center=np.array([0.0,0.0,0.0]), r
 
 
 class CubicCrystalDislocation(metaclass=ABCMeta):
-    '''
+    """
     Abstract base class for modelling a single dislocation
-    '''
+    """
 
     # Mandatory Attributes of CubicCrystalDislocation with no defaults
     # These should be set by the child dislocation class, or by CubicCrystalDislocation.__init__
@@ -2441,14 +2441,14 @@ def __init__(self, a, C11=None, C12=None, C44=None, C=None, symbol="W"):
         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
@@ -2458,7 +2458,7 @@ def init_solver(self, method="atomman"):
         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,
@@ -2474,7 +2474,7 @@ def get_solvers(self, method="atomman"):
         -------
         solvers: list
             List of functions f(positions) = displacements
-        '''
+        """
         has_atomman = False if self.atomman is None else True
 
         if type(method) != str:
@@ -2506,18 +2506,18 @@ def get_solvers(self, method="atomman"):
         return [self.solvers[method]]
 
     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)
 
         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() 
@@ -2542,9 +2542,9 @@ 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
@@ -2564,7 +2564,7 @@ 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)
 
@@ -2577,7 +2577,7 @@ def _build_supercell(self, targetx, targety):
 
         
         Returns a supercell of self.unit_cell with cell[0, 0] >= targetx and cell[1, 1] >= targety
-        '''
+        """
 
         base_cell = self.unit_cell
 
@@ -2590,7 +2590,7 @@ def _build_supercell(self, targetx, targety):
 
     def _build_bulk_cyl(self, radius, core_positions, fix_rad, extension,
                         fixed_points, self_consistent, method, verbose, cyl_mask=None, **kwargs):
-        '''
+        """
         Build bulk cylinder config from args supplied by self.build_cylinder
 
         Parameters
@@ -2617,13 +2617,13 @@ def _build_bulk_cyl(self, radius, core_positions, fix_rad, extension,
         cyl: ase.Atoms
             Bulk cyl, cut down based on radius, complete
             with FixAtoms constraints
-        '''
+        """
 
         def enforce_arr_shape(argname, arg, errs):
-            '''
+            """
             Enforce correct array dimensionality for 
             
-            '''
+            """
             if type(arg) == np.ndarray:
                 if arg.shape[-1] != 3:
                     # Needs to be 3d vectors
@@ -2760,7 +2760,7 @@ def plot_unit_cell(self, ms=250, ax=None):
 
     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])
@@ -2789,7 +2789,7 @@ def self_consistent_displacements(self, solvers, bulk_positions, core_positions,
         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, :]
 
@@ -2826,7 +2826,7 @@ def self_consistent_displacements(self, solvers, bulk_positions, core_positions,
 
     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])
@@ -2858,7 +2858,7 @@ def displacements(self, bulk_positions, core_positions, method="atomman",
         the dislocation structure
 
         Raises a RuntimeError if max_iter is reached without convergence
-        '''
+        """
         
         solvers = self.get_solvers(method)
 
@@ -2884,7 +2884,7 @@ def build_cylinder(self,
                        return_fix_mask=False,
                        cyl_mask=None,
                        **kwargs):
-        '''
+        """
         Build dislocation cylinder for single dislocation system
 
         Parameters
@@ -2928,7 +2928,7 @@ def build_cylinder(self,
             to cylinders of correct shape
         fix_mask: np.array
             If return_fix_mask=True, return the mask used for the FixAtoms constraint
-        '''
+        """
 
         core_positions = np.array([
             core_position + self.unit_cell_core_position
@@ -3108,11 +3108,11 @@ def build_impurity_cylinder(self, disloc, impurity, radius,
         return impurities_disloc
 
     def _smooth_kink_displacements(self, kink1, kink2, ref_bulk, smoothing_width, base_cell):
-        '''
+        """
         Use ReLU interpolation to smooth the displacment fields across a kink boundary between kink1 and kink2.
         Returns a copy of kink1, but with the displacement field at ref_bulk.positions[:, 2] > base_cell[2, 2] - smoothing width
         linearly interpolated. 
-        '''
+        """
 
         cell = base_cell
 
@@ -3333,7 +3333,7 @@ def ovito_dxa(self, disloc):
 
 
     def _plot_CLE_disloc_line(self, view, disloc, z_length, disloc_names=None, color=[0, 1, 0]):
-        '''Add dislocation line to the view as a cylinder and two cones.
+        """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.
         
@@ -3347,7 +3347,7 @@ def _plot_CLE_disloc_line(self, view, disloc, z_length, disloc_names=None, color
             Custom names for the dislocations 
         color: list
             [R, G, B] colour of the dislocation line
-        '''
+        """
 
         # 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"]:
@@ -3583,7 +3583,7 @@ def inf_line_transform(p, z_max):
     def view_cyl(self, system, scale=0.5, CNA_color=True, add_bonds=False,
                      line_color=[0, 1, 0], disloc_names=None, hide_arrows=False,
                      mode="dxa"):
-        '''
+        """
         NGLview-based visualisation tool for structures generated from the dislocation class
 
         Parameters
@@ -3608,7 +3608,7 @@ def view_cyl(self, system, scale=0.5, CNA_color=True, add_bonds=False,
             Mode for plotting dislocation line
             mode="dxa" uses ovito dxa to plot an exact line
             mode="cle" uses the cle solution used to generate the structure (if the dislocation line is straight)
-        '''
+        """
 
         from nglview import show_ase, ASEStructure
         from matscipy.utils import get_structure_types
@@ -3690,9 +3690,9 @@ def view_cyl(self, system, scale=0.5, CNA_color=True, add_bonds=False,
 
 
 class CubicCrystalDissociatedDislocation(CubicCrystalDislocation, metaclass=ABCMeta):
-    '''
+    """
     Abstract base class for modelling dissociated dislocation systems
-    '''
+    """
     # Inherits all slots from CubicCrystalDislocation as well
     __slots__ = ("left_dislocation", "right_dislocation")
 
@@ -3834,9 +3834,9 @@ 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()
@@ -3846,7 +3846,7 @@ def invert_burgers(self):
         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
 
@@ -3860,7 +3860,7 @@ def get_solvers(self, method="atomman"):
         -------
         solvers: list
             List of functions f(positions) = displacements
-        '''
+        """
         solvers = self.left_dislocation.get_solvers(method)
 
         solvers.extend(
@@ -3974,7 +3974,7 @@ def build_cylinder(self,
 class CubicCrystalDislocationQuadrupole(CubicCrystalDissociatedDislocation):
     burgers_dimensionless = np.zeros(3)
     def __init__(self, disloc_class, *args, **kwargs):
-        '''
+        """
         Initialise dislocation quadrupole class
 
         Arguments
@@ -3984,7 +3984,7 @@ def __init__(self, disloc_class, *args, **kwargs):
 
         *args, **kwargs
             Parameters fed to CubicCrystalDislocation.__init__()
-        '''
+        """
 
         if isinstance(disloc_class, CubicCrystalDislocation):
             disloc_cls = disloc_class.__class__
@@ -4035,7 +4035,7 @@ def glide_distance_dimensionless(self):
 
     def periodic_displacements(self, positions, v1, v2, core_positions, disp_tol=1e-3, max_neighs=15, 
                                verbose="periodic", **kwargs):
-        '''
+        """
         Calculate the stroh displacements for the periodic structure defined by 2D lattice vectors v1 & v2
         given the core positions.
 
@@ -4062,7 +4062,7 @@ def periodic_displacements(self, positions, v1, v2, core_positions, disp_tol=1e-
         -------
         disps: np.array
             displacements
-        '''
+        """
 
         ncores = core_positions.shape[0]
 
@@ -4137,7 +4137,7 @@ def periodic_displacements(self, positions, v1, v2, core_positions, disp_tol=1e-
 
     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
@@ -4160,7 +4160,7 @@ def build_quadrupole(self, glide_separation=4, partial_distance=0, extension=0,
             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)
@@ -4302,7 +4302,7 @@ def build_quadrupole(self, glide_separation=4, partial_distance=0, extension=0,
     
     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
 
@@ -4320,7 +4320,7 @@ def build_glide_quadrupoles(self, nims, invert_direction=False, glide_left=True,
         *args, **kwargs
             Fed to self.build_quadrupole()
         
-        '''
+        """
         
         glide_offsets = np.linspace(0, self.glide_distance, nims, endpoint=True)
 
@@ -4355,7 +4355,7 @@ def build_glide_quadrupoles(self, nims, invert_direction=False, glide_left=True,
         return images
     
     def _get_kink_quad_mask(self, ref_bulk, direction, precision=3):
-        '''
+        """
         Get an atom mask and tilted cell which generates a perfect bulk structure
         Mask is needed in building kinks to ensure continuity of atomic basis
 
@@ -4369,7 +4369,7 @@ def _get_kink_quad_mask(self, ref_bulk, direction, precision=3):
             
         returns a boolean mask, and a new cell
 
-        '''
+        """
 
         if self.left_dislocation.__class__ == BCCEdge111barDislocation:
             raise RuntimeError("Kink Quadrupoles do not currently work for BCCEdge111barDislocation")
@@ -4530,50 +4530,47 @@ def build_kink_quadrupole_from_glide_structs(self, ref_bulk, glide_structs, kink
 
         return kink_quad
     
-    #TODO: Fix this
-    def build_minimal_kink_quadrupole(self, n_kink=1, layer_decimal_precision=3, invert_direction=False, smooth_width=None,
+    def build_minimal_kink_quadrupole(self, n_kink=1, layer_decimal_precision=3, smooth_width=None,
                                 *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
+        <= 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
+        n_kink: int
+            Number of glide distances to kink out.
+            n_kink=1 for single kink in +x, n_kink=-2 for kink by two glide distances in -x.
+        layer_decimal_precision: 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
+            (required in order to ensure complete periodicity)
+            Top layer is defined as any atom with z component >= max(atom_z) - 10**(-layer_tol)
+            in fractional coordinates
+        smooth_width: float
+            Size (in Ang) of the region for displacement smoothing at each kink site.
+            Larger smoothing width assumes a broader kink structure.
+            Default is 0.5 * self.unit_cell.cell[2, 2]
         *args, **kwargs
-            Fed to self.build_quadrupole() & self.build_glide_quadrupoles
-        
-        '''
-
-        direction = -1 if invert_direction else 1
+            Fed to self.build_quadrupole()
+        """
 
         base_bulk, _ = self.build_quadrupole(*args, **kwargs)
 
         base_quad1, base_quad2 = self.build_glide_quadrupoles(2, glide_left=True, glide_right=True, 
-                                                 invert_direction=invert_direction, *args, **kwargs)
+                                                 invert_direction=False, *args, **kwargs)
 
-        kink_struct1 = base_quad1.copy() * (1, 1, z_reps)
-        kink_struct2 = base_quad2.copy() * (1, 1, z_reps)
+        kink_struct1 = base_quad1.copy()
+        kink_struct2 = base_quad2.copy()
         
-        bulk = base_bulk.copy() * (1, 1, z_reps)
+        bulk = base_bulk.copy()
 
         bulk.wrap()
         kink_struct1.wrap()
         kink_struct2.wrap()
 
         # Cell won't always be periodic, make sure we end up with something that is
-        mask, cell = self._get_kink_quad_mask(bulk, direction * n_kink, layer_decimal_precision)
+        mask, cell = self._get_kink_quad_mask(bulk, n_kink, layer_decimal_precision)
         bulk = bulk[mask]
         bulk.set_cell(cell, scale_atoms=False)
 
@@ -4592,7 +4589,7 @@ def build_minimal_kink_quadrupole(self, n_kink=1, layer_decimal_precision=3, inv
 
     def build_kink_quadrupole(self, kink_map=None, layer_decimal_precision=3, smooth_width=None,
                                 *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
@@ -4615,7 +4612,7 @@ def build_kink_quadrupole(self, kink_map=None, layer_decimal_precision=3, smooth
         *args, **kwargs
             Fed to self.build_quadrupole()
         
-        '''
+        """
 
         ref_bulk, glide_structs = self.build_kink_quadrupole_glide_structs(kink_map, *args, **kwargs)
 
@@ -4623,12 +4620,12 @@ def build_kink_quadrupole(self, kink_map=None, layer_decimal_precision=3, smooth
                                                              smooth_width, layer_decimal_precision)
     
     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)
 

From 45a56803f9dd0933fe79e88d91cf29a891b12e01 Mon Sep 17 00:00:00 2001
From: Thomas Rocke <45517493+thomas-rocke@users.noreply.github.com>
Date: Wed, 4 Sep 2024 11:59:41 +0100
Subject: [PATCH 34/46] DOC: Disloc mobility docs

---
 docs/applications/disloc_mobility.ipynb | 63 +++++++++++++++++++++++--
 1 file changed, 59 insertions(+), 4 deletions(-)

diff --git a/docs/applications/disloc_mobility.ipynb b/docs/applications/disloc_mobility.ipynb
index b7753051..51c9b086 100644
--- a/docs/applications/disloc_mobility.ipynb
+++ b/docs/applications/disloc_mobility.ipynb
@@ -237,7 +237,8 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "### Kink in quadrupoles"
+    "### Kink in quadrupoles\n",
+    "Like with dislocation cylinders, we can build networks of kinks in dislocation quadrupoles:"
    ]
   },
   {
@@ -246,12 +247,19 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "Si_quad_kink = Si_quad.build_kink_quadrupole_network(\n",
+    "Si_quad_kink = Si_quad.build_kink_quadrupole(\n",
     "    glide_separation=8,\n",
     "    kink_map=[0]*5 + [1]*5\n",
     ")\n",
     "\n",
-    "Si_quad.view_quad(Si_quad_kink)"
+    "Si_quad.view_quad(Si_quad_kink, mode=\"dxa\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "There is also another routine `build_minimal_kink_quadrupole` for building only the smallest possible kink structures. This is where the kink happens in a single burgers vector cell. Here, the `n_kink` parameter controls the number of glide distances covered by the kink, and the direction of the kink: `n_kink=2` builds a compressed version of the kink map `[0, 2]`, and `n_kink=-1` constructs a compressed version of `[0, -1]`."
    ]
   },
   {
@@ -259,7 +267,54 @@
    "execution_count": null,
    "metadata": {},
    "outputs": [],
-   "source": []
+   "source": [
+    "Si_quad_kink = Si_quad.build_minimal_kink_quadrupole(\n",
+    "    glide_separation=8,\n",
+    "    n_kink=1\n",
+    ")\n",
+    "\n",
+    "Si_quad.view_quad(Si_quad_kink, mode=\"dxa\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Improving the kink structures\n",
+    "So far, we have only used the Continuum Linear Elastic (CLE) solutions when building kink structures. As the kink structures are built by interpolating between glide structures, we can get a better approximation of the kink if we relax these glide structures before building the kink. \n",
+    "\n",
+    "To do this, we can replace a call to `build_kink_cylinder` with a call to `build_kink_glide_structs` to build the required glide structures, and then `build_kink_from_glide_cyls` to actually construct the kink. For quadrupoles, the equivalent is replacing `build_kink_quadrupole` with `build_kink_quadrupole_glide_structs` and `build_kink_quadrupole_from_glide_structs`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "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",
+    "from ase.optimize.precon import PreconLBFGS\n",
+    "\n",
+    "\n",
+    "calc = Manybody(**StillingerWeber(Holland_Marder_PRL_80_746_Si))\n",
+    "\n",
+    "kink_map = [0, 0, 0, 1, 1, 1]\n",
+    "\n",
+    "ref_bulk, glide_structs = Si_disloc.build_kink_glide_structs(kink_map, radius=40)\n",
+    "\n",
+    "# glide_structs has a length of 2, as only two unique values in the kink map\n",
+    "\n",
+    "for struct in glide_structs:\n",
+    "    struct.calc = calc\n",
+    "    opt = PreconLBFGS(struct)\n",
+    "    opt.run(1e-1)\n",
+    "\n",
+    "Si_kink = Si_disloc.build_kink_from_glide_cyls(ref_bulk, glide_structs, kink_map)\n",
+    "\n",
+    "Si_disloc.view_cyl(Si_kink)"
+   ]
   }
  ],
  "metadata": {

From d25f4fc237f745bf4431695270f0b48118e609c6 Mon Sep 17 00:00:00 2001
From: Thomas Rocke <45517493+thomas-rocke@users.noreply.github.com>
Date: Wed, 4 Sep 2024 14:20:35 +0100
Subject: [PATCH 35/46] DOC: Move quadrupole glide + kink to
 disloc_mobility.ipynb

---
 .../quadrupole_dislocations.ipynb             | 104 ------------------
 1 file changed, 104 deletions(-)

diff --git a/docs/applications/quadrupole_dislocations.ipynb b/docs/applications/quadrupole_dislocations.ipynb
index ce55f393..20a4aee6 100644
--- a/docs/applications/quadrupole_dislocations.ipynb
+++ b/docs/applications/quadrupole_dislocations.ipynb
@@ -67,110 +67,6 @@
     "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": {},

From c3960782138180e3d74909f481935dc3478b0836 Mon Sep 17 00:00:00 2001
From: Thomas Rocke <45517493+thomas-rocke@users.noreply.github.com>
Date: Wed, 4 Sep 2024 14:44:58 +0100
Subject: [PATCH 36/46] MAINT: Formatting of plasticity docs toc

---
 docs/applications/disloc_mobility.ipynb         | 2 +-
 docs/applications/plasticity.rst                | 5 +++--
 docs/applications/quadrupole_dislocations.ipynb | 7 ++++---
 3 files changed, 8 insertions(+), 6 deletions(-)

diff --git a/docs/applications/disloc_mobility.ipynb b/docs/applications/disloc_mobility.ipynb
index 51c9b086..fb43d1f4 100644
--- a/docs/applications/disloc_mobility.ipynb
+++ b/docs/applications/disloc_mobility.ipynb
@@ -280,7 +280,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Improving the kink structures\n",
+    "### Improving the kink structures\n",
     "So far, we have only used the Continuum Linear Elastic (CLE) solutions when building kink structures. As the kink structures are built by interpolating between glide structures, we can get a better approximation of the kink if we relax these glide structures before building the kink. \n",
     "\n",
     "To do this, we can replace a call to `build_kink_cylinder` with a call to `build_kink_glide_structs` to build the required glide structures, and then `build_kink_from_glide_cyls` to actually construct the kink. For quadrupoles, the equivalent is replacing `build_kink_quadrupole` with `build_kink_quadrupole_glide_structs` and `build_kink_quadrupole_from_glide_structs`."
diff --git a/docs/applications/plasticity.rst b/docs/applications/plasticity.rst
index 7055924e..ec83148b 100644
--- a/docs/applications/plasticity.rst
+++ b/docs/applications/plasticity.rst
@@ -18,9 +18,10 @@ Tutorials:
 ----------
 
 .. toctree::
-
+    
     cylinder_configurations.ipynb
-    multispecies_dislocations.ipynb
     quadrupole_dislocations.ipynb
+    disloc_mobility.ipynb
+    multispecies_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
index 20a4aee6..bd762b2f 100644
--- a/docs/applications/quadrupole_dislocations.ipynb
+++ b/docs/applications/quadrupole_dislocations.ipynb
@@ -5,6 +5,7 @@
    "metadata": {},
    "source": [
     "# Dislocation Quadrupoles\n",
+    "## Quadrupoles in Si Diamond\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",
@@ -74,7 +75,7 @@
     "## 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",
+    "### Screw in BCC 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)."
    ]
   },
@@ -120,8 +121,8 @@
    "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",
+    "### Dissociated Screw dislocations in FCC Ni\n",
+    "Dissociation of dislocation cores is another phenomenon that may 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."
    ]

From 3b2d446a2647edf7b3a37d3a5f6071f91fa7a8dd Mon Sep 17 00:00:00 2001
From: Thomas Rocke <45517493+thomas-rocke@users.noreply.github.com>
Date: Mon, 9 Sep 2024 14:31:00 +0100
Subject: [PATCH 37/46] ENH: More general kink_maps now work

---
 matscipy/dislocation.py | 74 ++++++++++++++++++++++-------------------
 1 file changed, 40 insertions(+), 34 deletions(-)

diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index 939f562c..add796d8 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -3165,19 +3165,25 @@ def build_kink_glide_structs(self, kink_map=None, *args, **kwargs):
         unique_kinks = np.sort(np.unique(kmap))
 
         # Make an empty list of the correct length
-        glide_structs = [0] * (krange+1)
+        glide_structs = [0] * len(unique_kinks)
+        
+        struct_map = []
 
+            
         fixed_points = np.array([
             [0, 0, 0],
             [self.glide_distance * krange, 0, 0]
         ])
 
-        for kink_pos in unique_kinks:
-            ref_bulk, glide_structs[kink_pos] = self.build_cylinder(*args, 
+        for i in range(len(kmap)):
+            struct_map.append(np.argmax(unique_kinks == kmap[i]))
+
+        for i, kink_pos in enumerate(unique_kinks):
+            ref_bulk, glide_structs[i] = self.build_cylinder(*args, 
                                                           fixed_points=fixed_points, 
                                                           core_position=np.array([kink_pos * self.glide_distance, 0 , 0]),
                                                           **kwargs)
-            
+
         # Deal with small displacements in boundary conditions caused by differences in disloc positions
         # by taking the average of the extremes
         fix_mask = glide_structs[0].arrays["fix_mask"]
@@ -3200,9 +3206,9 @@ def build_kink_glide_structs(self, kink_map=None, *args, **kwargs):
                 FixAtoms(mask=fix_mask)
             )
 
-        return ref_bulk, glide_structs
+        return ref_bulk, glide_structs, struct_map
     
-    def build_kink_from_glide_cyls(self, ref_bulk, glide_structs, kink_map=None, smooth_width=None):
+    def build_kink_from_glide_cyls(self, ref_bulk, glide_structs, struct_map, smooth_width=None):
         """
         Build a kink cylinder cell from the given kink map, using the structures contained in glide_structs
 
@@ -3211,26 +3217,18 @@ def build_kink_from_glide_cyls(self, ref_bulk, glide_structs, kink_map=None, smo
         glide_structs: list of ase Atoms
             Glide structures e.g. those produced by self.build_kink_glide_structs.
             The kink structure is constructed using glide_structs[kink_map[i]] for the ith cell.
-        kink_map: iterable of ints or None
-            Map of the location of the dislocation core in units of the glide vector
-            Default (kink_map=None) is a kink map of [0, 1]
-            See examples for more details.
+        struct_map: iterable of ints
+            Map of glide_structs back onto kink profile
+            Obtained by a call to build_kink_glide_structs
         smooth_width: float
             Size (in Ang) of the region for displacement smoothing at each kink site.
             Larger smoothing width assumes a broader kink structure.
             Default is 0.5 * self.unit_cell.cell[2, 2]
         
         """
-        # Deal with default kink_map value
-        if kink_map is None:
-            kink_map = [0, 1]
+        assert np.max(struct_map) == len(glide_structs) - 1
 
-        kmap = np.array(kink_map, dtype=int)
-        kmap -= np.min(kmap)
-
-        assert np.max(kmap) == len(glide_structs) - 1
-
-        kink_structs = [glide_structs[midx].copy() for midx in kmap]
+        kink_structs = [glide_structs[midx].copy() for midx in struct_map]
         
         # Smooth displacements across kink boundaries
         if smooth_width is None:
@@ -3251,7 +3249,7 @@ def build_kink_from_glide_cyls(self, ref_bulk, glide_structs, kink_map=None, smo
             kink_cyl = stack(kink_cyl, smoothed_kink_structs[i])
 
         # Add the correct FixAtoms constraint by concatenating glide struct fix masks in the correct order
-        fix_mask = np.array(sum([list(glide_structs[midx].arrays["fix_mask"]) for midx in kmap], []), dtype=bool)
+        fix_mask = np.array(sum([list(glide_structs[midx].arrays["fix_mask"]) for midx in struct_map], []), dtype=bool)
         kink_cyl.arrays["fix_mask"] = fix_mask
         kink_cyl.set_constraint(
             FixAtoms(mask=fix_mask)
@@ -3276,9 +3274,9 @@ def build_kink_cylinder(self, kink_map=None, smooth_width=None, *args, **kwargs)
             Extra arguments sent to self.build_cylinder
         
         """
-        ref_bulk, glide_structs = self.build_kink_glide_structs(kink_map, *args, **kwargs)
-
-        return self.build_kink_from_glide_cyls(ref_bulk, glide_structs, kink_map, smooth_width)
+        ref_bulk, glide_structs, struct_map = self.build_kink_glide_structs(kink_map, *args, **kwargs)
+        
+        return self.build_kink_from_glide_cyls(ref_bulk, glide_structs, struct_map, smooth_width=smooth_width)
     
     def ovito_dxa(self, disloc):
         """
@@ -4455,16 +4453,21 @@ def build_kink_quadrupole_glide_structs(self, kink_map=None, *args, **kwargs):
 
         unique_kinks = np.sort(np.unique(kmap))
 
-        glide_structs = [0] * (krange+1)
+        struct_map = []
 
-        for kink_pos in unique_kinks:
+        glide_structs = [0] * len(unique_kinks)
+
+        for i in range(len(kmap)):
+            struct_map.append(np.argmax(unique_kinks == kmap[i]))
+
+        for i, kink_pos in enumerate(unique_kinks):
             off = np.array([kink_pos * self.glide_distance, 0 , 0])
-            ref_bulk, glide_structs[kink_pos] = self.build_quadrupole(*args, 
+            ref_bulk, glide_structs[i] = self.build_quadrupole(*args, 
                                                           left_offset=off, right_offset=off,
                                                           **kwargs)
-        return ref_bulk, glide_structs
+        return ref_bulk, glide_structs, struct_map
 
-    def build_kink_quadrupole_from_glide_structs(self, ref_bulk, glide_structs, kink_map=None, smooth_width=None, layer_decimal_precision=3):
+    def build_kink_quadrupole_from_glide_structs(self, ref_bulk, glide_structs, kink_map, struct_map, smooth_width=None, layer_decimal_precision=3):
         """
         Build a kink cylinder cell from the given kink map, using the structures contained in glide_structs
 
@@ -4475,8 +4478,10 @@ def build_kink_quadrupole_from_glide_structs(self, ref_bulk, glide_structs, kink
             The kink structure is constructed using glide_structs[kink_map[i]] for the ith cell.
         kink_map: iterable of ints or None
             Map of the location of the dislocation core in units of the glide vector
-            Default (kink_map=None) is a kink map of [0, 1]
-            See examples for more details.
+            Use the same kink_map as the call to build_kink_quadrupole_glide_structs
+        struct_map: iterable of ints
+            Map of glide_structs back onto kink profile
+            Obtained by a call to build_kink_quadrupole_glide_structs
         smooth_width: float
             Size (in Ang) of the region for displacement smoothing at each kink site.
             Larger smoothing width assumes a broader kink structure.
@@ -4488,6 +4493,7 @@ def build_kink_quadrupole_from_glide_structs(self, ref_bulk, glide_structs, kink
             in fractional coordinates
         
         """
+
         # Deal with default kink_map value
         if kink_map is None:
             kink_map = [0, 1]
@@ -4495,9 +4501,9 @@ def build_kink_quadrupole_from_glide_structs(self, ref_bulk, glide_structs, kink
         kmap = np.array(kink_map, dtype=int)
         kmap -= np.min(kmap)
 
-        assert np.max(kmap) == len(glide_structs) - 1
+        assert np.max(struct_map) == len(glide_structs) - 1
 
-        kink_structs = [glide_structs[midx].copy() for midx in kmap]
+        kink_structs = [glide_structs[midx].copy() for midx in struct_map]
         
         # Smooth displacements across kink boundaries
         if smooth_width is None:
@@ -4614,9 +4620,9 @@ def build_kink_quadrupole(self, kink_map=None, layer_decimal_precision=3, smooth
         
         """
 
-        ref_bulk, glide_structs = self.build_kink_quadrupole_glide_structs(kink_map, *args, **kwargs)
+        ref_bulk, glide_structs, struct_map = self.build_kink_quadrupole_glide_structs(kink_map, *args, **kwargs)
 
-        return self.build_kink_quadrupole_from_glide_structs(ref_bulk, glide_structs, kink_map,
+        return self.build_kink_quadrupole_from_glide_structs(ref_bulk, glide_structs, kink_map, struct_map,
                                                              smooth_width, layer_decimal_precision)
     
     def view_quad(self, system, *args, **kwargs):

From 2f11c971a56d4bc49227e91402d3e5906830976c Mon Sep 17 00:00:00 2001
From: Thomas Rocke <45517493+thomas-rocke@users.noreply.github.com>
Date: Mon, 9 Sep 2024 15:23:28 +0100
Subject: [PATCH 38/46] MAINT: Kink generation functions return reference bulk

---
 matscipy/dislocation.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index add796d8..9a402680 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -3255,7 +3255,7 @@ def build_kink_from_glide_cyls(self, ref_bulk, glide_structs, struct_map, smooth
             FixAtoms(mask=fix_mask)
         )
 
-        return kink_cyl
+        return bulk, kink_cyl
     
     def build_kink_cylinder(self, kink_map=None, smooth_width=None, *args, **kwargs):
         """
@@ -4591,7 +4591,7 @@ def build_minimal_kink_quadrupole(self, n_kink=1, layer_decimal_precision=3, smo
 
         # Smooth displacements across the periodic boundary
         kink = self._smooth_kink_displacements(kink_struct1, kink_struct2, bulk, smooth_width, base_bulk.cell[:, :])
-        return kink
+        return bulk, kink
 
     def build_kink_quadrupole(self, kink_map=None, layer_decimal_precision=3, smooth_width=None,
                                 *args, **kwargs):

From 2b2dd38a6463811f05c8756a96551ededad69d86 Mon Sep 17 00:00:00 2001
From: Thomas Rocke <45517493+thomas-rocke@users.noreply.github.com>
Date: Mon, 9 Sep 2024 16:09:00 +0100
Subject: [PATCH 39/46] MAINT: Quad kink also returns bulk

---
 matscipy/dislocation.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index 9a402680..cadfa03b 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -4371,7 +4371,6 @@ def _get_kink_quad_mask(self, ref_bulk, direction, precision=3):
 
         if self.left_dislocation.__class__ == BCCEdge111barDislocation:
             raise RuntimeError("Kink Quadrupoles do not currently work for BCCEdge111barDislocation")
-            pass
 
         bulk = ref_bulk.copy()
         tilt_bulk = ref_bulk.copy()
@@ -4534,7 +4533,7 @@ def build_kink_quadrupole_from_glide_structs(self, ref_bulk, glide_structs, kink
         kink_quad = kink_quad[mask]
         kink_quad.set_cell(cell, scale_atoms=False)
 
-        return kink_quad
+        return bulk, kink_quad
     
     def build_minimal_kink_quadrupole(self, n_kink=1, layer_decimal_precision=3, smooth_width=None,
                                 *args, **kwargs):

From 74587ae389bf56b05533bd4c047b3a368bb6386c Mon Sep 17 00:00:00 2001
From: Thomas Rocke <45517493+thomas-rocke@users.noreply.github.com>
Date: Mon, 9 Sep 2024 17:14:31 +0100
Subject: [PATCH 40/46] BUG: Fix kink building not properly reading struct_map

---
 matscipy/dislocation.py | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index cadfa03b..94dd5e54 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -3234,13 +3234,13 @@ def build_kink_from_glide_cyls(self, ref_bulk, glide_structs, struct_map, smooth
         if smooth_width is None:
             smooth_width = 0.5 * self.unit_cell.cell[2, 2]
 
-        N = len(kink_structs)
+        N = len(struct_map)
 
         smoothed_kink_structs = []
 
         for i in range(N):
             smoothed_kink_structs.append(
-                self._smooth_kink_displacements(kink_structs[i], kink_structs[(i+1) % N], ref_bulk, smooth_width, ref_bulk.cell[:, :])
+                self._smooth_kink_displacements(kink_structs[struct_map[i]], kink_structs[struct_map[(i+1) % N]], ref_bulk, smooth_width, ref_bulk.cell[:, :])
             )
 
         kink_cyl = smoothed_kink_structs[0].copy()
@@ -3255,7 +3255,7 @@ def build_kink_from_glide_cyls(self, ref_bulk, glide_structs, struct_map, smooth
             FixAtoms(mask=fix_mask)
         )
 
-        return bulk, kink_cyl
+        return ref_bulk * (1, 1, N), kink_cyl
     
     def build_kink_cylinder(self, kink_map=None, smooth_width=None, *args, **kwargs):
         """
@@ -4508,13 +4508,13 @@ def build_kink_quadrupole_from_glide_structs(self, ref_bulk, glide_structs, kink
         if smooth_width is None:
             smooth_width = 0.5 * self.unit_cell.cell[2, 2]
 
-        N = len(kink_structs)
+        N = len(struct_map)
 
         smoothed_kink_structs = []
 
         for i in range(N):
             smoothed_kink_structs.append(
-                self._smooth_kink_displacements(kink_structs[i], kink_structs[(i+1) % N], ref_bulk, smooth_width, ref_bulk.cell[:, :])
+                self._smooth_kink_displacements(kink_structs[struct_map[i]], kink_structs[struct_map[(i+1) % N]], ref_bulk, smooth_width, ref_bulk.cell[:, :])
             )
 
         kink_quad = smoothed_kink_structs[0].copy()
@@ -4530,6 +4530,9 @@ def build_kink_quadrupole_from_glide_structs(self, ref_bulk, glide_structs, kink
         # Cell won't always be periodic, make sure we end up with something that is
         mask, cell = self._get_kink_quad_mask(bulk, direction, layer_decimal_precision)
 
+        bulk = bulk[mask]
+        bulk.set_cell(cell)
+        
         kink_quad = kink_quad[mask]
         kink_quad.set_cell(cell, scale_atoms=False)
 

From 6a7f5f16eee8634f30fa9756be92c757e0f39261 Mon Sep 17 00:00:00 2001
From: Thomas Rocke <45517493+thomas-rocke@users.noreply.github.com>
Date: Mon, 9 Sep 2024 17:14:47 +0100
Subject: [PATCH 41/46] TEST: Unittests for cyl & quad kinks

---
 tests/test_dislocation.py | 105 +++++++++++++++++++++++++++++++++++---
 1 file changed, 97 insertions(+), 8 deletions(-)

diff --git a/tests/test_dislocation.py b/tests/test_dislocation.py
index dad5ba9e..a53a6fb9 100644
--- a/tests/test_dislocation.py
+++ b/tests/test_dislocation.py
@@ -794,6 +794,41 @@ def test_methods(self, disloc, subtests):
                 except AssertionError as e:
 
                     raise AssertionError(f"Displacements from {method} did not match {base_method}")
+                
+    def test_kink_round_trip(self, disloc, subtests):
+        self.set_up_cls(disloc)
+        d = self.test_cls(self.alat, self.C11, self.C12, self.C44, symbol=self.symbol)
+
+        kink_map = [0, 1]
+
+        ref_bulk, glide_structs, struct_map = d.build_kink_glide_structs(kink_map=kink_map, radius=30)
+        bulk, kink1 = d.build_kink_from_glide_cyls(ref_bulk, glide_structs, struct_map)
+
+        bulk, kink2 = d.build_kink_cylinder(kink_map=kink_map, radius=30)
+
+        assert len(kink1) == len(kink2)
+        assert len(bulk) == len(kink2)
+
+        np.testing.assert_array_almost_equal(kink1.positions, kink2.positions)
+    
+    def test_kink_equiv_maps(self, disloc, subtests):
+        self.set_up_cls(disloc)
+        d = self.test_cls(self.alat, self.C11, self.C12, self.C44, symbol=self.symbol)
+
+        equiv_kinks = [
+            [[0, 1], [-1, 0]],
+            [[0, 2], [-1, 1]]
+        ]
+
+        for kmap1, kmap2 in equiv_kinks:
+            with subtests.test(f"Comparing equivalent kink maps {kmap1} and {kmap2}"):
+                bulk, kc1 = d.build_kink_cylinder(kink_map=kmap1, radius=30)
+                bulk, kc2 = d.build_kink_cylinder(kink_map=kmap2, radius=30)
+
+                assert len(kc1) == len(kc2)
+                assert len(bulk) == len(kc2)
+
+                np.testing.assert_array_almost_equal(kc1.positions, kc2.positions)
 
 @pytest.mark.parametrize("disloc", cubic_perfect_dislocs())
 class TestCubicCrystalDislocation(BaseTestCubicCrystalDislocation):
@@ -904,19 +939,73 @@ def test_glide_configs(self, disloc, subtests):
             self.assertAtomsAlmostEqual(configs[0], ini_quad)
             self.assertAtomsAlmostEqual(configs[1], fin_quad)
 
-    @pytest.mark.skip()
-    def test_single_kink_quadrupole(self, disloc):
-        '''
-        Validate that generation of single kink structures runs without errors
-        '''
+    def test_quad_kink_round_trip(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)
+
+        kink_map = [0, 1]
+
+        is_111bar = disloc.name == "1/2<11-1> edge" and disloc.crystalstructure.lower() == "bcc"
+
+        if is_111bar:
+            # Check RuntimeError is raised for BCC 111bar disloc
+            with pytest.raises(RuntimeError):
+                bulk, kc1 = d.build_kink_quadrupole(kink_map=[0, 1], glide_separation=8)
+        else:
+            ref_bulk, glide_structs, struct_map = d.build_kink_quadrupole_glide_structs(kink_map=kink_map, glide_separation=8)
+            bulk, kink1 = d.build_kink_quadrupole_from_glide_structs(ref_bulk, glide_structs, kink_map, struct_map)
+            
+            bulk, kink2 = d.build_kink_quadrupole(kink_map=kink_map, glide_separation=8)
+            
+
+            assert len(kink1) == len(kink2)
+            assert len(bulk) == len(kink2)
+
+            np.testing.assert_array_almost_equal(kink1.positions, kink2.positions)
+    
+    def test_quad_kink_equiv_maps(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)
 
+        equiv_kinks = [
+            [[0, 1], [-1, 0]],
+            [[0, 2], [-1, 1]]
+        ]
 
-        kink_struct = d.build_kink_quadrupole(z_reps=4, glide_separation=4)
+        is_111bar = disloc.name == "1/2<11-1> edge" and disloc.crystalstructure.lower() == "bcc"
 
-        # Check that the z vector has been tilted to get an infinite array of periodic single kinks
-        self.assertNotAlmostEqual(kink_struct.cell[2, 0], 0.0)
+        if is_111bar:
+            # Check RuntimeError is raised for BCC 111bar disloc
+            with pytest.raises(RuntimeError):
+                bulk, kc1 = d.build_kink_quadrupole(kink_map=[0, 1], glide_separation=8)
+        else:
+            for kmap1, kmap2 in equiv_kinks:
+                with subtests.test(f"Comparing equivalent kink maps {kmap1} and {kmap2}"):
+                    bulk, kc1 = d.build_kink_quadrupole(kink_map=kmap1, glide_separation=8)
+                    bulk, kc2 = d.build_kink_quadrupole(kink_map=kmap2, glide_separation=8)
+                    
+                    assert len(kc1) == len(kc2)
+                    assert len(bulk) == len(kc2)
+
+                    np.testing.assert_array_almost_equal(kc1.positions, kc2.positions)
+
+    def test_quad_minimal_kink(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)
+
+        n_kinks = np.arange(-2, 2)
+
+        is_111bar = disloc.name == "1/2<11-1> edge" and disloc.crystalstructure.lower() == "bcc"
+
+        if is_111bar:
+            # Check RuntimeError is raised for BCC 111bar disloc
+            with pytest.raises(RuntimeError):
+                bulk, kc1 = d.build_kink_quadrupole(kink_map=[0, 1], glide_separation=8)
+        else:
+            for n_kink in n_kinks:
+                with subtests.test(f"Building minimal {n_kink=} cell"):
+                    bulk, kink = d.build_minimal_kink_quadrupole(n_kink, glide_separation=8)
+                    assert len(bulk) == len(kink)
 
 
 

From f834c4ff3abba46aebd4c3cfa1e1872ae3ae4529 Mon Sep 17 00:00:00 2001
From: Thomas Rocke <45517493+thomas-rocke@users.noreply.github.com>
Date: Mon, 9 Sep 2024 17:15:35 +0100
Subject: [PATCH 42/46] DOC: Updating disloc mobility docs to work with
 returned bulk

---
 docs/applications/disloc_mobility.ipynb | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/docs/applications/disloc_mobility.ipynb b/docs/applications/disloc_mobility.ipynb
index fb43d1f4..b1ec51b2 100644
--- a/docs/applications/disloc_mobility.ipynb
+++ b/docs/applications/disloc_mobility.ipynb
@@ -225,7 +225,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "Si_kink = Si_disloc.build_kink_cylinder(\n",
+    "Si_bulk, Si_kink = Si_disloc.build_kink_cylinder(\n",
     "                    radius=20,\n",
     "                    kink_map= [0] * 5 + [1] * 5\n",
     "                )\n",
@@ -247,7 +247,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "Si_quad_kink = Si_quad.build_kink_quadrupole(\n",
+    "Si_quad_bulk, Si_quad_kink = Si_quad.build_kink_quadrupole(\n",
     "    glide_separation=8,\n",
     "    kink_map=[0]*5 + [1]*5\n",
     ")\n",
@@ -268,7 +268,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "Si_quad_kink = Si_quad.build_minimal_kink_quadrupole(\n",
+    "Si_quad_bulk, Si_quad_kink = Si_quad.build_minimal_kink_quadrupole(\n",
     "    glide_separation=8,\n",
     "    n_kink=1\n",
     ")\n",
@@ -302,7 +302,7 @@
     "\n",
     "kink_map = [0, 0, 0, 1, 1, 1]\n",
     "\n",
-    "ref_bulk, glide_structs = Si_disloc.build_kink_glide_structs(kink_map, radius=40)\n",
+    "ref_bulk, glide_structs, struct_map = Si_disloc.build_kink_glide_structs(kink_map, radius=40)\n",
     "\n",
     "# glide_structs has a length of 2, as only two unique values in the kink map\n",
     "\n",
@@ -311,7 +311,7 @@
     "    opt = PreconLBFGS(struct)\n",
     "    opt.run(1e-1)\n",
     "\n",
-    "Si_kink = Si_disloc.build_kink_from_glide_cyls(ref_bulk, glide_structs, kink_map)\n",
+    "Si_bulk, Si_kink = Si_disloc.build_kink_from_glide_cyls(ref_bulk, glide_structs, kink_map, struct_map)\n",
     "\n",
     "Si_disloc.view_cyl(Si_kink)"
    ]

From 4f87a8e0c3836fd74915bbf05858a203f0b29b2b Mon Sep 17 00:00:00 2001
From: Thomas Rocke <45517493+thomas-rocke@users.noreply.github.com>
Date: Tue, 10 Sep 2024 11:32:08 +0100
Subject: [PATCH 43/46] ENH: DXA colours for BCC & FCC

---
 matscipy/dislocation.py | 20 ++++++++++++++++++--
 1 file changed, 18 insertions(+), 2 deletions(-)

diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index 94dd5e54..4cd9eb6b 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -3474,7 +3474,9 @@ def inf_line_transform(p, z_max):
                 return np.array(points)
         from fractions import Fraction
         
-        # Ovito colours for dislocation lines (as of v 3.10.1) 
+        # Ovito colours for dislocation lines (as of v 3.10.1)
+        # Keys are abs of burgers as fractions, sorted in descending order
+        # (to collapse rotational symmetries) 
         disloc_colours = {
             "diamond" : {
                 "default" : [230, 51, 51], # Red
@@ -3482,6 +3484,20 @@ def inf_line_transform(p, z_max):
                 "1/3 1/6 1/6" : [0, 255, 0], # Green; 30 & 90 degree Partials
                 "1/6 1/6 0" : [255, 0, 255], # Pink
                 "1/3 1/3 1/3" : [0, 255, 255] # Light Blue
+            },
+            "bcc" : {
+                "default" : [230, 51, 51], # Red
+                "1/2 1/2 1/2" : [0, 255, 0], # Green; Screw/Edge 111
+                "1 0 0" : [255, 77, 204], # Pink; Edge 100
+                "1 1 0" : [51, 128, 255] # Blue
+            },
+            "fcc" : {
+                "default" : [230, 51, 51], # Red
+                "1/2 1/2 0" : [51, 51, 255], # Blue; Screw/Edge 110
+                "1/3 1/6 1/6" : [0, 255, 0], # Green; Shockley Partials
+                "1/6 1/6 0" : [255, 0, 255], # Pink; Stair Rod
+                "1/3 0 0" : [255, 255, 0], # Yellow; Hirth
+                "1/3 1/3 1/3" : [0, 255, 255] # Cyan; Frank
             }
         }
 
@@ -4532,7 +4548,7 @@ def build_kink_quadrupole_from_glide_structs(self, ref_bulk, glide_structs, kink
 
         bulk = bulk[mask]
         bulk.set_cell(cell)
-        
+
         kink_quad = kink_quad[mask]
         kink_quad.set_cell(cell, scale_atoms=False)
 

From e5d9addd7a69cbb1bb9a69acaccb54008de83ac0 Mon Sep 17 00:00:00 2001
From: Thomas Rocke <45517493+thomas-rocke@users.noreply.github.com>
Date: Tue, 10 Sep 2024 11:51:03 +0100
Subject: [PATCH 44/46] BUG: Fix issue with kink structure improvement docs

---
 docs/applications/disloc_mobility.ipynb | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/docs/applications/disloc_mobility.ipynb b/docs/applications/disloc_mobility.ipynb
index b1ec51b2..cab1cbac 100644
--- a/docs/applications/disloc_mobility.ipynb
+++ b/docs/applications/disloc_mobility.ipynb
@@ -300,7 +300,7 @@
     "\n",
     "calc = Manybody(**StillingerWeber(Holland_Marder_PRL_80_746_Si))\n",
     "\n",
-    "kink_map = [0, 0, 0, 1, 1, 1]\n",
+    "kink_map = [0] * 5 + [1] * 5\n",
     "\n",
     "ref_bulk, glide_structs, struct_map = Si_disloc.build_kink_glide_structs(kink_map, radius=40)\n",
     "\n",
@@ -308,10 +308,10 @@
     "\n",
     "for struct in glide_structs:\n",
     "    struct.calc = calc\n",
-    "    opt = PreconLBFGS(struct)\n",
+    "    opt = PreconLBFGS(struct, logfile=None)\n",
     "    opt.run(1e-1)\n",
     "\n",
-    "Si_bulk, Si_kink = Si_disloc.build_kink_from_glide_cyls(ref_bulk, glide_structs, kink_map, struct_map)\n",
+    "Si_bulk, Si_kink = Si_disloc.build_kink_from_glide_cyls(ref_bulk, glide_structs, struct_map)\n",
     "\n",
     "Si_disloc.view_cyl(Si_kink)"
    ]

From e3f653edc5a6676aef9c8ae42645c53e54283e10 Mon Sep 17 00:00:00 2001
From: Thomas Rocke <45517493+thomas-rocke@users.noreply.github.com>
Date: Tue, 10 Sep 2024 12:34:11 +0100
Subject: [PATCH 45/46] BUG: Fix kink structures not being built correctly

---
 matscipy/dislocation.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index 4cd9eb6b..c00a31a5 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -3234,13 +3234,13 @@ def build_kink_from_glide_cyls(self, ref_bulk, glide_structs, struct_map, smooth
         if smooth_width is None:
             smooth_width = 0.5 * self.unit_cell.cell[2, 2]
 
-        N = len(struct_map)
+        N = len(kink_structs)
 
         smoothed_kink_structs = []
 
         for i in range(N):
             smoothed_kink_structs.append(
-                self._smooth_kink_displacements(kink_structs[struct_map[i]], kink_structs[struct_map[(i+1) % N]], ref_bulk, smooth_width, ref_bulk.cell[:, :])
+                self._smooth_kink_displacements(kink_structs[i], kink_structs[(i+1) % N], ref_bulk, smooth_width, ref_bulk.cell[:, :])
             )
 
         kink_cyl = smoothed_kink_structs[0].copy()
@@ -4524,13 +4524,13 @@ def build_kink_quadrupole_from_glide_structs(self, ref_bulk, glide_structs, kink
         if smooth_width is None:
             smooth_width = 0.5 * self.unit_cell.cell[2, 2]
 
-        N = len(struct_map)
+        N = len(kink_structs)
 
         smoothed_kink_structs = []
 
         for i in range(N):
             smoothed_kink_structs.append(
-                self._smooth_kink_displacements(kink_structs[struct_map[i]], kink_structs[struct_map[(i+1) % N]], ref_bulk, smooth_width, ref_bulk.cell[:, :])
+                self._smooth_kink_displacements(kink_structs[i], kink_structs[(i+1) % N], ref_bulk, smooth_width, ref_bulk.cell[:, :])
             )
 
         kink_quad = smoothed_kink_structs[0].copy()

From 8da8bdb693b137995229a43f89322030df571869 Mon Sep 17 00:00:00 2001
From: Thomas Rocke <45517493+thomas-rocke@users.noreply.github.com>
Date: Tue, 10 Sep 2024 16:11:11 +0100
Subject: [PATCH 46/46] DOC: Reduce verbosity of kink quad output in docs

---
 docs/applications/disloc_mobility.ipynb | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/docs/applications/disloc_mobility.ipynb b/docs/applications/disloc_mobility.ipynb
index cab1cbac..c4274543 100644
--- a/docs/applications/disloc_mobility.ipynb
+++ b/docs/applications/disloc_mobility.ipynb
@@ -216,7 +216,7 @@
     "\n",
     "Since the kink map is just a list of integers, it can be more convenient to exploit list addition, and specify kink maps in a form similar to `[0] * 5 + [1] * 5`, which would be identical to the input of `[0, 0, 0, 0, 0, 1, 1, 1, 1, 1]`.\n",
     "\n",
-    "### Kinks in cylinders"
+    "### Kink in cylinders"
    ]
   },
   {
@@ -249,7 +249,8 @@
    "source": [
     "Si_quad_bulk, Si_quad_kink = Si_quad.build_kink_quadrupole(\n",
     "    glide_separation=8,\n",
-    "    kink_map=[0]*5 + [1]*5\n",
+    "    kink_map=[0]*5 + [1]*5,\n",
+    "    verbose=False\n",
     ")\n",
     "\n",
     "Si_quad.view_quad(Si_quad_kink, mode=\"dxa\")"
@@ -270,7 +271,8 @@
    "source": [
     "Si_quad_bulk, Si_quad_kink = Si_quad.build_minimal_kink_quadrupole(\n",
     "    glide_separation=8,\n",
-    "    n_kink=1\n",
+    "    n_kink=1,\n",
+    "    verbose=False\n",
     ")\n",
     "\n",
     "Si_quad.view_quad(Si_quad_kink, mode=\"dxa\")"