Repository: libAtoms/matscipy (language: python)

Issue:
Adding mixed mode NCFlex on top of 3D NCFlex
First pull request to reconcile all of the NCFlex changes made on different branches together, with a view to make a pull request between the ParallelNCFlex branch and master.


Discussion / hints:

Fix patch:
diff --git a/matscipy/cauchy_born.py b/matscipy/cauchy_born.py
index 2918ec73..ae7acab7 100644
--- a/matscipy/cauchy_born.py
+++ b/matscipy/cauchy_born.py
@@ -1599,3 +1599,9 @@ def load_regression_model(self):
         """Loads the regression model from file: CB_model.npy"""
         self.RM = RegressionModel()
         self.RM.load()
+
+    def get_model(self):
+        return self.RM.model
+    
+    def set_model(self, model):
+        self.RM.model = model
\ No newline at end of file
diff --git a/matscipy/fracture_mechanics/clusters.py b/matscipy/fracture_mechanics/clusters.py
index c771c84e..aa87a91b 100644
--- a/matscipy/fracture_mechanics/clusters.py
+++ b/matscipy/fracture_mechanics/clusters.py
@@ -24,8 +24,163 @@
 
 from ase.lattice.cubic import Diamond, FaceCenteredCubic, SimpleCubic, BodyCenteredCubic
 from matscipy.neighbours import neighbour_list
+import ase.io
+####
+
+def get_alpha_period(cryst):
+    pos = cryst.get_positions()
+    sx, sy, sz = cryst.cell.diagonal()
+    xpos = pos[:,0] - sx/2
+    ypos = pos[:,1] - sy/2
+    #find the closest y atoms by finding which atoms lie within 1e-2 of the min
+    #filter out all atoms with x less than 0
+    xmask = xpos>0
+    closest_y_mask = np.abs(ypos)<(np.min(np.abs(ypos[xmask]))+(1e-2))
+    #find the x positions of these atoms
+    closest_x = xpos[closest_y_mask&xmask]
+    #sort these atoms and find the largest x gap
+    sorted_x = np.sort(closest_x)
+    diffs = np.diff(sorted_x)
+    alpha_period = np.sum(np.unique(np.round(np.diff(sorted_x),decimals=4)))
+    return alpha_period
+
+def generate_3D_structure(cryst_2D,nzlayer,el,a0,lattice,crack_surface,crack_front,shift=np.array([0.0,0.0,0.0]),cb=None,switch_sublattices=False):
+    alpha_period = get_alpha_period(cryst_2D)
+    #make a single layer of cell 
+    single_cell = lattice(el,a0,[1,1,1],crack_surface,crack_front,cb=cb,shift=shift,switch_sublattices=switch_sublattices)
+    ase.io.write('single_cell.xyz',single_cell)
+    cell_x_period = single_cell.get_cell()[0,0]
+    print('NUM KINK PER CELL,', cell_x_period/alpha_period)
+    #original x,y dimensions
+    big_cryst = cryst_2D*[1,1,nzlayer]
+    og_cell = big_cryst.get_cell()
+    h = og_cell[2,2]
+    theta = np.arctan(cell_x_period/h)
+    pos = big_cryst.get_positions()
+    xpos = pos[:,0]
+    zpos = pos[:,2]
+    cell_xlength = og_cell[0,0]
+    ztantheta = zpos*np.tan(theta)
+    mask = xpos>(cell_xlength - ztantheta)
+    pos[mask,0] = pos[mask,0] - cell_xlength
+    big_cryst.set_positions(pos)
+    new_cell = og_cell.copy()
+    new_cell[2,0] = -cell_x_period
+    ase.io.write('big_cryst_pre_rotation.xyz',big_cryst)
+    big_cryst.set_cell(new_cell)
+    ase.io.write('big_cryst.xyz',big_cryst)
+
+    # Define the rotation matrix
+    R = np.array([[np.cos(theta), 0, np.sin(theta)],
+            [0, 1, 0],
+            [-np.sin(theta), 0, np.cos(theta)]])
+
+    # Get the cell of big_cryst
+    cell = big_cryst.get_cell()
+
+    # Rotate the cell using the rotation matrix
+    rotated_cell = np.dot(cell, R.T)
+
+    # Set the new cell of big_cryst
+    big_cryst.set_cell(rotated_cell,scale_atoms=True)
+    return big_cryst, theta
+
+def generate_3D_cubic_111(cryst_2D, nzlayer, el, a0, lattice, crack_surface, crack_front, shift=np.array([0.0,0.0,0.0]), cb=None, switch_sublattices=False):
+    """
+    Generate a kink-periodic cell, using the high symmetry of the 111 cubic surface
+    to reduce the number of kinks in the cell
+
+    Parameters:
+    -----------
+    cryst_2D : ASE atoms object
+        The 2D cryst structure for the crack to use as a template for the 3D structure.
+    nzlayer : int
+        The number of layers in the z direction.
+    el : str
+        The element symbol to use for the atoms.
+    a0 : float
+        The lattice constant.
+    lattice : function
+        The function to use for generating the lattice.
+    crack_surface : str
+        The surface on which the crack is located.
+    crack_front : str
+        The direction of the crack front.
+    shift : numpy.ndarray, optional
+        The shift vector to apply to the lattice. Default is [0.0, 0.0, 0.0].
+    cb : float or None, optional
+        The concentration of vacancies to introduce in the lattice. Default is None.
+    switch_sublattices : bool, optional
+        Whether to switch the sublattices. Default is False.
+
+    Returns:
+    --------
+    big_cryst : ase.Atoms
+        The 3D cubic crystal structure with a (111) surface and a crack.
+    """
+    # in this case, one can make use of the high symmetry of the 111 surface
+    # in order to reduce the number of kinks
+    alpha_period = get_alpha_period(cryst_2D)
+    #make a single layer of cell 
+    single_cell = lattice(el,a0,[1,1,1],crack_surface,crack_front,cb=cb,shift=shift,switch_sublattices=switch_sublattices)
+    #ase.io.write('single_cell.xyz',single_cell)
+    cell_x_period = single_cell.get_cell()[0,0]
+    nkink_per_cell = int(np.round((cell_x_period/alpha_period)))
+    print('NUM KINK PER CELL,', nkink_per_cell)
+    #original x,y dimensions
+    
+    if nkink_per_cell == 2:
+        big_cryst = cryst_2D*[1,1,nzlayer+1]
+        #ase.io.write('big_cryst_larger.xyz',big_cryst)
+        big_cryst.translate([0,0,-0.01])
+        big_cryst.wrap()
+        #ase.io.write('big_cryst_translated.xyz',big_cryst)
+        cell_x_period = cell_x_period/2
+        # remove the extra half layer in z
+        extended_cell = big_cryst.get_cell()
+        extended_cell[2,2] -= (single_cell.get_cell()[2,2]/2)
+        big_cryst.set_cell(extended_cell)
+        #ase.io.write('big_cryst_smaller_cell.xyz',big_cryst)
+        pos = big_cryst.get_positions()
+        mask = pos[:,2]>extended_cell[2,2]
+        del big_cryst[mask]
+        #ase.io.write('big_cryst_atoms_removed.xyz',big_cryst)
+        #big_cryst.wrap()
+    else:
+        big_cryst = cryst_2D*[1,1,nzlayer]
+
+    og_cell = big_cryst.get_cell()
+    h = og_cell[2,2]
+    theta = np.arctan(cell_x_period/h)
+    pos = big_cryst.get_positions()
+    xpos = pos[:,0]
+    zpos = pos[:,2]
+    cell_xlength = og_cell[0,0]
+    ztantheta = zpos*np.tan(theta)
+    mask = xpos>(cell_xlength - ztantheta)
+    pos[mask,0] = pos[mask,0] - cell_xlength
+    big_cryst.set_positions(pos)
+    new_cell = og_cell.copy()
+    new_cell[2,0] = -cell_x_period
+    #ase.io.write('big_cryst_pre_rotation.xyz',big_cryst)
+    big_cryst.set_cell(new_cell)
+    #ase.io.write('big_cryst.xyz',big_cryst)
+
+    # Define the rotation matrix
+    R = np.array([[np.cos(theta), 0, np.sin(theta)],
+            [0, 1, 0],
+            [-np.sin(theta), 0, np.cos(theta)]])
+
+    # Get the cell of big_cryst
+    cell = big_cryst.get_cell()
+
+    # Rotate the cell using the rotation matrix
+    rotated_cell = np.dot(cell, R.T)
+
+    # Set the new cell of big_cryst
+    big_cryst.set_cell(rotated_cell,scale_atoms=True)
+    return big_cryst, theta
 
-###
 
 def set_groups(a, n, skin_x, skin_y, central_x=-1./2, central_y=-1./2,
                invert_central=False):
@@ -59,7 +214,7 @@ def set_groups(a, n, skin_x, skin_y, central_x=-1./2, central_y=-1./2,
 
     a.set_array('groups', g)
 
-def set_regions(cryst, r_I, cutoff, r_III, extended_far_field=False,extended_region_I=False,exclude_surface=False):
+def set_regions(cryst, r_I, cutoff, r_III, extended_far_field=False,extended_region_I=False,exclude_surface=False,sort_type='r_theta_z'):
     sx, sy, sz = cryst.cell.diagonal()
     x, y = cryst.positions[:, 0], cryst.positions[:, 1]
     cx, cy = sx/2, sy/2
@@ -120,8 +275,41 @@ def set_regions(cryst, r_I, cutoff, r_III, extended_far_field=False,extended_reg
     # keep only cylinder defined by regions I - IV
     cryst = cryst[regionI | regionII | regionIII | regionIV]
 
-    # order by radial distance from tip
-    order = r[regionI | regionII | regionIII | regionIV ].argsort()
+    if sort_type =='radial':
+
+        print('Warning: Using old method of sorting by radial distance from tip in x-y plane.')
+
+        # order by radial distance from tip
+        order = r[regionI | regionII | regionIII | regionIV ].argsort()
+    
+    elif sort_type == 'region':
+
+        # sort by regions, retaining original bulk order within each region
+        region = cryst.arrays['region']
+        region_numbers = np.unique(region) # returns sorted region numbers
+        order = np.array([ index for n in region_numbers for index in np.where(region==n)[0] ])
+
+    elif sort_type == 'r_theta_z':
+            # sort by r, theta, z
+            sx, sy, sz = cryst.cell.diagonal()
+            x, y, z = cryst.positions[:, 0], cryst.positions[:, 1], cryst.positions[:, 2]
+            cx, cy = sx/2, sy/2
+            r = np.sqrt((x - cx)**2 + (y - cy)**2)
+            # get r from the centre of the cell
+            theta = np.arctan2(y-cy, x-cx)
+            #first, discretely bin r
+            r_sorted = np.sort(r)
+            r_sorted_index = np.argsort(r)
+            r_diff = np.diff(r_sorted)
+            for i,r_diff_val in enumerate(r_diff):
+                if r_diff_val<1e-3:
+                    r[r_sorted_index[i+1]] = r[r_sorted_index[i]]
+        
+            order = np.lexsort((z, theta, r))
+
+    else:
+        raise ValueError('sort_type must be one of "radial", "region", or "r_theta_z"')
+    
     cryst = cryst[order]
     return cryst
 
diff --git a/matscipy/fracture_mechanics/crack.py b/matscipy/fracture_mechanics/crack.py
index 1ead5a7a..a265b269 100644
--- a/matscipy/fracture_mechanics/crack.py
+++ b/matscipy/fracture_mechanics/crack.py
@@ -140,9 +140,9 @@ def _init_crack(self):
         self.mu2_q1 = self.mu2 * self.q1
 
 
-    def displacements(self, r, theta, k):
+    def displacements(self, r, theta, kI, kII=0.0):
         """
-        Displacement field in mode I fracture. Positions are passed in cylinder
+        Displacement field in mode I/II fracture. Positions are passed in cylinder
         coordinates.
 
         Parameters
@@ -162,20 +162,26 @@ def displacements(self, r, theta, k):
             Displacements normal to the plane of the crack.
         """
 
-        h1 = k * np.sqrt(2.0*r/math.pi)
+        h1 = kI * np.sqrt(2.0*r/math.pi)
 
         h2 = np.sqrt( np.cos(theta) + self.mu2*np.sin(theta) )
         h3 = np.sqrt( np.cos(theta) + self.mu1*np.sin(theta) )
 
-        u = h1*( self.inv_mu1_mu2*( self.mu1_p2*h2 - self.mu2_p1*h3 ) ).real
-        v = h1*( self.inv_mu1_mu2*( self.mu1_q2*h2 - self.mu2_q1*h3 ) ).real
+        u_I = h1*( self.inv_mu1_mu2*( self.mu1_p2*h2 - self.mu2_p1*h3 ) ).real
+        v_I = h1*( self.inv_mu1_mu2*( self.mu1_q2*h2 - self.mu2_q1*h3 ) ).real
 
+        h4 = kII*np.sqrt(2.0*r/math.pi)
+        u_II = h4*(self.inv_mu1_mu2*(self.p2*h2 - self.p1*h3)).real
+        v_II = h4*(self.inv_mu1_mu2*(self.q2*h2 - self.q1*h3)).real
+
+        u = u_I + u_II
+        v = v_I + v_II
         return u, v
 
 
-    def deformation_gradient(self, r, theta, k):
+    def deformation_gradient(self, r, theta, kI, kII=0.0):
         """
-        Deformation gradient tensor in mode I fracture. Positions are passed in
+        Deformation gradient tensor in mode I/II fracture. Positions are passed in
         cyclinder coordinates.
 
         Parameters
@@ -202,17 +208,31 @@ def deformation_gradient(self, r, theta, k):
             perpendicular to the plane.
         """
 
-        f = k / np.sqrt(2*math.pi*r)
+        f_I = kI / np.sqrt(2*math.pi*r)
 
         h1 = (self.mu1*self.mu2)*self.inv_mu1_mu2
         h2 = np.sqrt( np.cos(theta) + self.mu2*np.sin(theta) )
         h3 = np.sqrt( np.cos(theta) + self.mu1*np.sin(theta) )
 
-        du_dx = f*( self.inv_mu1_mu2*( self.mu1_p2/h2 - self.mu2_p1/h3 ) ).real
-        du_dy = f*( h1*( self.p2/h2 - self.p1/h3 ) ).real
+        du_dx_I = f_I*( self.inv_mu1_mu2*( self.mu1_p2/h2 - self.mu2_p1/h3 ) ).real
+        du_dy_I = f_I*( h1*( self.p2/h2 - self.p1/h3 ) ).real
+
+        dv_dx_I = f_I*( self.inv_mu1_mu2*( self.mu1_q2/h2 - self.mu2_q1/h3 ) ).real
+        dv_dy_I = f_I*( h1*( self.q2/h2 - self.q1/h3 ) ).real
+
+        f_II = kII / np.sqrt(2*math.pi*r)
+        h4 = self.inv_mu1_mu2
+        du_dx_II = f_II*( self.inv_mu1_mu2*( self.p2/h2 - self.p1/h3 ) ).real
+        du_dy_II = f_II*( h4*( (self.mu2*self.p2)/h2 - (self.mu1*self.p1)/h3 ) ).real
+
+        dv_dx_II = f_II*( self.inv_mu1_mu2*( self.q2/h2 - self.q1/h3 ) ).real
+        dv_dy_II = f_II*( h4*( (self.mu2*self.q2)/h2 - (self.mu1*self.q1)/h3 ) ).real
+
+        du_dx = du_dx_I + du_dx_II
+        du_dy = du_dy_I + du_dy_II
+        dv_dx = dv_dx_I + dv_dx_II
+        dv_dy = dv_dy_I + dv_dy_II
 
-        dv_dx = f*( self.inv_mu1_mu2*( self.mu1_q2/h2 - self.mu2_q1/h3 ) ).real
-        dv_dy = f*( h1*( self.q2/h2 - self.q1/h3 ) ).real
 
         # We need to add unity matrix to turn this into the deformation gradient
         # tensor.
@@ -221,9 +241,9 @@ def deformation_gradient(self, r, theta, k):
 
         return np.transpose([[du_dx, du_dy], [dv_dx, dv_dy]])
 
-    def stresses(self, r, theta, k):
+    def stresses(self, r, theta, kI, kII=0.0):
         """
-        Stress field in mode I fracture. Positions are passed in cylinder
+        Stress field in mode I/II fracture. Positions are passed in cylinder
         coordinates.
 
         Parameters
@@ -247,15 +267,26 @@ def stresses(self, r, theta, k):
             Off-diagonal component of the stress tensor.
         """
 
-        f = k / np.sqrt(2.0*math.pi*r)
+        f_I = kI / np.sqrt(2.0*math.pi*r)
 
         h1 = (self.mu1*self.mu2)*self.inv_mu1_mu2
         h2 = np.sqrt( np.cos(theta) + self.mu2*np.sin(theta) )
         h3 = np.sqrt( np.cos(theta) + self.mu1*np.sin(theta) )
 
-        sig_x  = f*(h1*(self.mu2/h2 - self.mu1/h3)).real
-        sig_y  = f*(self.inv_mu1_mu2*(self.mu1/h2 - self.mu2/h3)).real
-        sig_xy = f*(h1*(1/h3 - 1/h2)).real
+        sig_x_I  = f_I*(h1*(self.mu2/h2 - self.mu1/h3)).real
+        sig_y_I  = f_I*(self.inv_mu1_mu2*(self.mu1/h2 - self.mu2/h3)).real
+        sig_xy_I = f_I*(h1*(1/h3 - 1/h2)).real
+
+        f_II = kII / np.sqrt(2.0*math.pi*r)
+        h4 = self.inv_mu1_mu2
+
+        sig_x_II = f_II*(h4*((self.mu2**2)/h2 - (self.mu1**2)/h3)).real
+        sig_y_II = f_II*(h4*((1)/h2 - (1)/h3)).real
+        sig_xy_II = f_II*(h4*((self.mu1)/h3 - (self.mu2)/h2)).real
+
+        sig_x = sig_x_I + sig_x_II
+        sig_y = sig_y_I + sig_y_II
+        sig_xy = sig_xy_I + sig_xy_II
 
         return sig_x, sig_y, sig_xy
 
@@ -305,7 +336,7 @@ def k1gsqG(self):
 
 ###
 
-def displacement_residuals(r0, crack, x, y, ref_x, ref_y, k, power=1):
+def displacement_residuals(r0, crack, x, y, ref_x, ref_y, kI, kII=0, power=1):
     """
     Return actual displacement field minus ideal displacement field,
     divided by r**alpha.
@@ -317,7 +348,7 @@ def displacement_residuals(r0, crack, x, y, ref_x, ref_y, k, power=1):
     dy = ref_y - y0
     abs_dr = np.sqrt(dx*dx+dy*dy)
     theta = np.arctan2(dy, dx)
-    u2x, u2y = crack.displacements_from_cylinder_coordinates(abs_dr, theta, k)
+    u2x, u2y = crack.displacements_from_cylinder_coordinates(abs_dr, theta, kI, kII=kII)
     if abs(power) < 1e-12:
         power_of_abs_dr = 1.0
     else:
@@ -325,8 +356,8 @@ def displacement_residuals(r0, crack, x, y, ref_x, ref_y, k, power=1):
     return (u1x - u2x)/power_of_abs_dr, (u1y - u2y)/power_of_abs_dr
 
 
-def displacement_residual(r0, crack, x, y, ref_x, ref_y, k, mask=None, power=1):
-    dux, duy = displacement_residuals(r0, crack, x, y, ref_x, ref_y, k,
+def displacement_residual(r0, crack, x, y, ref_x, ref_y, kI, kII=0, mask=None, power=1):
+    dux, duy = displacement_residuals(r0, crack, x, y, ref_x, ref_y, kI, kII=kII,
                                       power=power)
     if mask is None:
         return np.transpose([dux, duy]).flatten()
@@ -334,8 +365,7 @@ def displacement_residual(r0, crack, x, y, ref_x, ref_y, k, mask=None, power=1):
         return np.transpose([dux[mask], duy[mask]]).flatten()
 
 
-def deformation_gradient_residuals(r0, crack, x, y, cur, ref_x, ref_y, ref, k,
-                                   cutoff):
+def deformation_gradient_residuals(r0, crack, x, y, cur, ref_x, ref_y, ref, kI, cutoff, kII=0):
     """
     Return actual displacement field minus ideal displacement field.
     """
@@ -352,15 +382,15 @@ def deformation_gradient_residuals(r0, crack, x, y, cur, ref_x, ref_y, ref, k,
     F1 = np.array([[F1xx, F1xy], [F1yx, F1yy]]).T
 
     F2 = crack.deformation_gradient(ref.positions[:, 0], ref.positions[:, 1],
-                                    x0, y0, k)
+                                    x0, y0, kI, kII=kII)
 
     return F1-F2
 
 
-def deformation_gradient_residual(r0, crack, x, y, cur, ref_x, ref_y, ref, k,
-                                  cutoff, mask=None):
+def deformation_gradient_residual(r0, crack, x, y, cur, ref_x, ref_y, ref, kI,
+                                  cutoff, kII=0, mask=None):
     dF = deformation_gradient_residuals(r0, crack, x, y, cur, ref_x, ref_y, ref,
-                                        k, cutoff)
+                                        kI, cutoff, kII=kII)
     if mask is None:
         return dF.flatten()
     else:
@@ -450,23 +480,23 @@ def k1gsqG(self):
         return self.crack.k1gsqG()
 
 
-    def displacements_from_cylinder_coordinates(self, r, theta, k):
+    def displacements_from_cylinder_coordinates(self, r, theta, kI, kII=0):
         """
-        Displacement field in mode I fracture from cylindrical coordinates.
+        Displacement field in mode I/II fracture from cylindrical coordinates.
         """
-        return self.crack.displacements(r, theta, k)
+        return self.crack.displacements(r, theta, kI, kII=kII)
 
 
-    def displacements_from_cartesian_coordinates(self, dx, dy, k):
+    def displacements_from_cartesian_coordinates(self, dx, dy, kI, kII=0):
         """
-        Displacement field in mode I fracture from cartesian coordinates.
+        Displacement field in mode I/II fracture from cartesian coordinates.
         """
         abs_dr = np.sqrt(dx*dx+dy*dy)
         theta = np.arctan2(dy, dx)
-        return self.displacements_from_cylinder_coordinates(abs_dr, theta, k)
+        return self.displacements_from_cylinder_coordinates(abs_dr, theta, kI, kII=kII)
 
 
-    def displacements(self, ref_x, ref_y, x0, y0, k):
+    def displacements(self, ref_x, ref_y, x0, y0, kI, kII=0):
         """
         Displacement field for a list of cartesian positions.
 
@@ -492,26 +522,26 @@ def displacements(self, ref_x, ref_y, x0, y0, k):
         """
         dx = ref_x - x0
         dy = ref_y - y0
-        return self.displacements_from_cartesian_coordinates(dx, dy, k)
+        return self.displacements_from_cartesian_coordinates(dx, dy, kI, kII=kII)
 
 
-    def deformation_gradient_from_cylinder_coordinates(self, r, theta, k):
+    def deformation_gradient_from_cylinder_coordinates(self, r, theta, kI, kII=0):
         """
         Displacement field in mode I fracture from cylindrical coordinates.
         """
-        return self.crack.deformation_gradient(r, theta, k)
+        return self.crack.deformation_gradient(r, theta, kI, kII=kII)
 
 
-    def deformation_gradient_from_cartesian_coordinates(self, dx, dy, k):
+    def deformation_gradient_from_cartesian_coordinates(self, dx, dy, kI, kII=0):
         """
         Displacement field in mode I fracture from cartesian coordinates.
         """
         abs_dr = np.sqrt(dx*dx+dy*dy)
         theta = np.arctan2(dy, dx)
-        return self.deformation_gradient_from_cylinder_coordinates(abs_dr, theta, k)
+        return self.deformation_gradient_from_cylinder_coordinates(abs_dr, theta, kI, kII=kII)
 
 
-    def deformation_gradient(self, ref_x, ref_y, x0, y0, k):
+    def deformation_gradient(self, ref_x, ref_y, x0, y0, kI, kII=0):
         """
         Deformation gradient for a list of cartesian positions.
 
@@ -537,9 +567,9 @@ def deformation_gradient(self, ref_x, ref_y, x0, y0, k):
         """
         dx = ref_x - x0
         dy = ref_y - y0
-        return self.deformation_gradient_from_cartesian_coordinates(dx, dy, k)
+        return self.deformation_gradient_from_cartesian_coordinates(dx, dy, kI, kII=kII)
 
-    def crack_tip_position(self, x, y, ref_x, ref_y, x0, y0, k, mask=None,
+    def crack_tip_position(self, x, y, ref_x, ref_y, x0, y0, kI, kII=0, mask=None,
                            residual_func=displacement_residual,
                            method='Powell', return_residuals=False):
         """
@@ -587,29 +617,29 @@ def crack_tip_position(self, x, y, ref_x, ref_y, x0, y0, k, mask=None,
             mask = np.ones(len(x), dtype=bool)
         if method == 'leastsq':
             (x1, y1), ier = leastsq(residual_func, (x0, y0),
-                                    args=(self, x, y, ref_x, ref_y, k, mask))
+                                    args=(self, x, y, ref_x, ref_y, kI, kII, mask))
             if ier not in [ 1, 2, 3, 4 ]:
                 raise RuntimeError('Could not find crack tip')
         else:
             opt = minimize(lambda *args: (residual_func(*args)**2).sum(),
-                           (x0, y0), args=(self, x, y, ref_x, ref_y, k, mask),
+                           (x0, y0), args=(self, x, y, ref_x, ref_y, kI, kII, mask),
                            method=method)
             if not opt.success:
                 raise RuntimeError('Could not find crack tip. Reason: {}'
                                    .format(opt.message))
             x1, y1 = opt.x
         if return_residuals:
-            return x1, y1, residual_func((x0, y0), self, x, y, ref_x, ref_y, k)
+            return x1, y1, residual_func((x0, y0), self, x, y, ref_x, ref_y, kI, kII=kII)
         else:
             return x1, y1
 
 
-    def _residual_y(self, y0, x0, x, y, ref_x, ref_y, k, mask):
-        dux, duy = self.displacement_residuals(x, y, ref_x, ref_y, x0, y0, k)
+    def _residual_y(self, y0, x0, x, y, ref_x, ref_y, kI, mask, kII=0):
+        dux, duy = self.displacement_residuals(x, y, ref_x, ref_y, x0, y0, kI, kII=kII)
         return dux[mask]*dux[mask]+duy[mask]*duy[mask]
 
 
-    def crack_tip_position_y(self, x, y, ref_x, ref_y, x0, y0, k, mask=None):
+    def crack_tip_position_y(self, x, y, ref_x, ref_y, x0, y0, kI, kII=0, mask=None):
         """
         Return an estimate of the y-coordinate of the real crack tip position
         assuming the stress intensity factor is k.
@@ -641,7 +671,7 @@ def crack_tip_position_y(self, x, y, ref_x, ref_y, x0, y0, k, mask=None):
         if mask is None:
             mask = np.ones(len(a), dtype=bool)
         ( y0, ), ier = leastsq(self._residual_y, y0,
-                          args=(x0, x, y, ref_x, ref_y, k, mask))
+                          args=(x0, x, y, ref_x, ref_y, kI, kII, mask))
         if ier not in [ 1, 2, 3, 4 ]:
             raise RuntimeError('Could not find crack tip')
         return y0
@@ -651,7 +681,7 @@ def scale_displacements(self, x, y, ref_x, ref_y, old_k, new_k):
         """
         Rescale atomic positions from stress intensity factor old_k to the new
         stress intensity factor new_k. This is useful for extrapolation of
-        relaxed positions.
+        relaxed positions. Note - this only works for pure mode I fracture
 
         Parameters
         ----------
@@ -678,23 +708,23 @@ def scale_displacements(self, x, y, ref_x, ref_y, old_k, new_k):
         return ref_x + new_k/old_k*(x-ref_x), ref_y + new_k/old_k*(y-ref_y)
 
 
-    def stresses_from_cylinder_coordinates(self, r, theta, k):
+    def stresses_from_cylinder_coordinates(self, r, theta, kI, kII=0):
         """
         Stress field in mode I fracture from cylindrical coordinates
         """
-        return self.crack.stresses(r, theta, k)
+        return self.crack.stresses(r, theta, kI, kII=kII)
 
 
-    def stresses_from_cartesian_coordinates(self, dx, dy, k):
+    def stresses_from_cartesian_coordinates(self, dx, dy, kI, kII=0):
         """
         Stress field in mode I fracture from cartesian coordinates
         """
         abs_dr = np.sqrt(dx*dx+dy*dy)
         theta = np.arctan2(dy, dx)
-        return self.stresses_from_cylinder_coordinates(abs_dr, theta, k)
+        return self.stresses_from_cylinder_coordinates(abs_dr, theta, kI, kII=kII)
 
 
-    def stresses(self, ref_x, ref_y, x0, y0, k):
+    def stresses(self, ref_x, ref_y, x0, y0, kI, kII=0):
         """
         Stress field for a list of cartesian positions
 
@@ -724,7 +754,7 @@ def stresses(self, ref_x, ref_y, x0, y0, k):
         """
         dx = ref_x - x0
         dy = ref_y - y0
-        sig_x, sig_y, sig_xy = self.stresses_from_cartesian_coordinates(dx, dy, k)
+        sig_x, sig_y, sig_xy = self.stresses_from_cartesian_coordinates(dx, dy, kI, kII=kII)
         return sig_x, sig_y, sig_xy
 
 
@@ -735,11 +765,11 @@ class SinclairCrack:
     Sinclair, J. E. The Influence of the Interatomic Force Law and of Kinks on the
         Propagation of Brittle Cracks. Philos. Mag. 31, 647–671 (1975)
     """
-    def __init__(self, crk, cryst, calc, k, alpha=0.0, vacuum=6.0,
+    def __init__(self, crk, cryst, calc, k, kI=None, kII=0, alpha=0.0, vacuum=6.0,
                  variable_alpha=True, variable_k=False,
                  alpha_scale=None, k_scale=None,
                  extended_far_field=False,rI=0.0,rIII=0.0,cutoff=0.0,
-                 incl_rI_f_alpha=False):
+                 incl_rI_f_alpha=False,is_3D=False,theta_3D=0.0,cont_k='k1'):
         """
 
         Parameters
@@ -747,7 +777,9 @@ def __init__(self, crk, cryst, calc, k, alpha=0.0, vacuum=6.0,
         crk - instance of CubicCrystalCrack
         cryst - Atoms object representing undeformed crystal
         calc - ASE-compatible calculator
-        k - stress intensity factor
+        k - stress intensity factor in mode I (deprecated name for kI)
+        kI - stress intensity factor in mode I
+        kII - stress intensity factor in mode II
         alpha - crack tip position
         vacuum - amount of vacuum to add to unit cell
         variable_alpha - if True, crack tip position can vary (flexible BCs)
@@ -780,9 +812,23 @@ def __init__(self, crk, cryst, calc, k, alpha=0.0, vacuum=6.0,
         self.k_scale = k_scale
         self.extended_far_field = extended_far_field
 
+
+
         self.u = np.zeros((self.N1, 3))
         self.alpha = alpha
-        self.k = k
+
+        self.is_3D = is_3D
+        if is_3D:
+            self.theta_3D = theta_3D
+            self.beta = self.alpha/np.sin(theta_3D)
+
+        if kI is not None:
+            self.kI = kI
+        else:
+            self.kI = k
+        
+        self.kII = kII
+        
         self.rI = rI
         self.rIII = rIII
         self.cutoff = cutoff
@@ -806,7 +852,7 @@ def __init__(self, crk, cryst, calc, k, alpha=0.0, vacuum=6.0,
 
         a0 = self.atoms.copy()
         self.x0 = a0.get_positions()
-        #self.E0 = self.calc.get_potential_energies(a0)[self.regionI_II].sum()
+        self.E0 = self.calc.get_potential_energies(a0)[self.regionI_II].sum()
         a0_II_III = a0[self.regionII | self.regionIII]
         f0bar = self.calc.get_forces(a0_II_III)
         self.f0bar = f0bar[a0_II_III.arrays['region'] == 2]
@@ -816,8 +862,19 @@ def __init__(self, crk, cryst, calc, k, alpha=0.0, vacuum=6.0,
         self.f_alpha_correction = 0
         self.alpha0 = alpha
         self.force_calls = 0
+        self.cont_k = cont_k
+        self.follow_G_contour = False
 
 
+    #alias self.k to self.kI such as not to break old code
+    @property
+    def k(self):
+        return self.kI
+    
+    @k.setter
+    def k(self, value):
+        self.kI = value
+
     def pack(self, u, alpha, k):
         dofs = list(u.reshape(-1))
         if self.variable_alpha:
@@ -842,17 +899,72 @@ def unpack(self, x, reshape=False, defaults=None):
         if self.variable_k:
             k = float(x[offset])
         else:
-            k = defaults.get('k', self.k)
+            if self.cont_k == 'k1':
+                k = defaults.get('kI', self.kI)
+            elif self.cont_k == 'k2':
+                k = defaults.get('kII', self.kII)
         return u, alpha, k
 
     def get_dofs(self):
-        return self.pack(self.u, self.alpha, self.k)
+        if self.is_3D:
+            if self.cont_k == 'k1':
+                return self.pack(self.u,self.beta,self.kI)
+            elif self.cont_k == 'k2':
+                return self.pack(self.u,self.alpha,self.kII)
+        else:
+            if self.cont_k == 'k1':
+                return self.pack(self.u,self.alpha,self.kI)
+            elif self.cont_k == 'k2':
+                return self.pack(self.u, self.alpha,self.kII)
 
     def set_dofs(self, x):
-        self.u[:], self.alpha, self.k = self.unpack(x, reshape=True)
-        print('alpha',self.alpha)
+        #if set_dofs gets passed an array that contains k1 and k2, it should
+        #set the one that isn't being used for continuation, and remove it from x
+        #before passing it forwards
+        #get the total number of dofs one would expect
+        prev_v_k = self.variable_k
+        prev_v_alph = self.variable_alpha
+        self.variable_k = True
+        self.variable_alpha = True
+        ndof = len(self)
+        self.variable_k = prev_v_k
+        self.variable_alpha = prev_v_alph
+        if len(x) == ndof+1: #array contains k1 and k2
+            if self.cont_k == 'k1':
+                #set kII and delete it from x
+                idx = -1
+                self.kII = x[idx]
+            elif self.cont_k == 'k2':
+                #set kI and delete it from x
+                idx = -2
+                self.kI = x[idx]
+            x = np.delete(x,idx)
+        
+        if self.is_3D:
+            if self.cont_k == 'k1':
+                self.u[:], self.beta, self.kI = self.unpack(x, reshape=True)
+            elif self.cont_k == 'k2':
+                self.u[:], self.beta, self.kII = self.unpack(x, reshape=True)
+            self.alpha = self.beta*np.sin(self.theta_3D)
+            print('beta',self.beta)
+            print('alpha',self.alpha)
+        else:
+            if self.cont_k == 'k1':
+                self.u[:], self.alpha, self.kI = self.unpack(x, reshape=True)
+            elif self.cont_k == 'k2':
+                self.u[:], self.alpha, self.kII = self.unpack(x, reshape=True)
+            print('alpha',self.alpha)
         if self.variable_k:
-            print('k',self.k)
+            print('kI',self.kI)
+            print('kII',self.kII)
+        
+        #adjust non-continuation K to follow constant G contour
+        
+        if self.follow_G_contour:
+            if self.cont_k == 'k2':
+                self.kI = np.sqrt(self.kI0**2 + self.kII0**2 - self.kII**2)
+            elif self.cont_k == 'k1':
+                self.kII = np.sqrt(self.kI0**2 + self.kII0**2 - self.kI**2)
         #self.k = self.k-(np.abs(self.alpha-self.alpha0)*self.k1g)
         self.update_atoms()
 
@@ -864,25 +976,31 @@ def __len__(self):
             N_dofs += 1
         return N_dofs
 
-    def u_cle(self, alpha=None):
+    def u_cle(self, alpha=None, kI=None, kII=None):
         """
         Returns CLE displacement solution at current (alpha, K)
 
-        Note that this is NOT pre-multipled by the stress intensity factor k
+        Note that this is now pre-multipled by the stress intensity factor k
         """
         if alpha is None:
             alpha = self.alpha
+        if kI is None:
+            kI = self.kI
+        if kII is None:
+            kII = self.kII
+        
         tip_x = self.cryst.cell.diagonal()[0] / 2.0 + alpha
         tip_y = self.cryst.cell.diagonal()[1] / 2.0
         ux, uy = self.crk.displacements(self.cryst.positions[:, 0],
                                         self.cryst.positions[:, 1],
-                                        tip_x, tip_y, 1.0)
+                                        tip_x, tip_y, self.kI, self.kII)
         u = np.c_[ux, uy, np.zeros_like(ux)] # convert to 3D field
         return u
 
 
     def fit_cle(self, r_fit=20.0, variable_alpha=True, variable_k=True, x0=None,
                 grid=None):
+        #Fit the CLE field - note that this only works for pure mode I fracture
         def residuals(x, mask):
             idx = 0
             if variable_alpha:
@@ -894,10 +1012,10 @@ def residuals(x, mask):
                 k = x[idx]
                 idx += 1
             else:
-                k = self.k
+                k = self.kI
             u = np.zeros((len(self.atoms), 3))
             u[self.regionI] = self.u
-            du = (self.k * self.u_cle(self.alpha) + u - k * self.u_cle(alpha))
+            du = (self.u_cle(self.alpha,kI=self.kI) + u - self.u_cle(alpha,kI=k))
             return du[mask, :].reshape(-1)
 
         mask = self.r < r_fit
@@ -906,7 +1024,7 @@ def residuals(x, mask):
             if variable_alpha:
                 x0.append(self.alpha)
             if variable_k:
-                x0.append(self.k)
+                x0.append(self.kI)
         if grid:
             alpha_grid, k_grid = grid
             vals = np.zeros((len(alpha_grid), len(k_grid)))
@@ -921,31 +1039,37 @@ def residuals(x, mask):
                 raise RuntimeError('CLE fit failed')
         return res
 
-    def get_deformation_gradient(self,x,y,k=None,de=0):
+    def get_deformation_gradient(self,x,y,kI=None,kII=None,de=0):
         alpha = self.alpha + de #take finite differences if necessary
-        if k is None:
-            k = self.k
+        if kI is None:
+            kI = self.kI
+        if kII is None:
+            kII = self.kII
         tip_x = self.cryst.cell.diagonal()[0] / 2.0 + alpha
         tip_y = self.cryst.cell.diagonal()[1] / 2.0
-        dg = self.crk.deformation_gradient(x, y ,tip_x, tip_y, k)
+        dg = self.crk.deformation_gradient(x, y ,tip_x, tip_y, kI, kII)
         return dg
 
     def set_shiftmask(self,radial_dist):
         self.shiftmask = self.r>radial_dist
     
-    def update_atoms(self):
+    def update_atoms(self,use_alpha_3D=False):
         """
-        Update self.atoms from degrees of freedom (self.u, self.alpha, self.k)
+        Update self.atoms from degrees of freedom (self.u, self.alpha, self.kI, self.kII)
         """
+        if self.is_3D and (not use_alpha_3D):
+            #set self.alpha according to self.beta
+            self.alpha = self.beta*np.sin(self.theta_3D)
         self.atoms.set_pbc([False, False, True])
         self.atoms.calc = self.calc
 
-        self.atoms.info['k'] = self.k
+        self.atoms.info['kI'] = self.kI
+        self.atoms.info['kII'] = self.kII
         self.atoms.info['alpha'] = self.alpha
 
         # x = x_cryst + K * u_cle + u
         self.atoms.positions[:, :] = self.cryst.positions
-        self.atoms.positions[:, :] += self.k * self.u_cle()
+        self.atoms.positions[:, :] += self.u_cle()
         self.atoms.positions[self.regionI, :] += self.u
         if self.crk.cauchy_born is not None: #if the crack has a multilattice cauchy-born object
             #get rotation matrix 
@@ -955,7 +1079,7 @@ def update_atoms(self):
             # very important to pass cryst rather than atoms here as the displacement gradient field
             # is found from the positions of the original atoms, not the deformed atoms
             shifts = self.crk.cauchy_born.predict_shifts(A,self.cryst,\
-                F_func=self.get_deformation_gradient, coordinates='cart2D',method='regression',k=self.k)
+                F_func=self.get_deformation_gradient, coordinates='cart2D',method='regression',kI=self.kI, kII=self.kII)
             print('done!')
             #apply shifts
             self.crk.cauchy_born.apply_shifts(self.atoms,shifts,mask=self.shiftmask)
@@ -978,12 +1102,13 @@ def get_f_alpha_correction(self):
         self.atoms.set_pbc([False, False, True])
         self.atoms.calc = self.calc
 
-        self.atoms.info['k'] = self.k
+        self.atoms.info['kI'] = self.kI
+        self.atoms.info['kII'] = self.kII
         self.atoms.info['alpha'] = self.alpha
 
         # x = x_cryst + K * u_cle + u
         self.atoms.positions[:, :] = self.cryst.positions
-        self.atoms.positions[:, :] += self.k * self.u_cle()
+        self.atoms.positions[:, :] += self.u_cle()
         if self.crk.cauchy_born is not None: #if the crack has a multilattice cauchy-born object
             #get rotation matrix 
             A = np.transpose(self.crk.RotationMatrix)
@@ -992,7 +1117,7 @@ def get_f_alpha_correction(self):
             # very important to pass cryst rather than atoms here as the displacement gradient field
             # is found from the positions of the original atoms, not the deformed atoms
             shifts = self.crk.cauchy_born.predict_shifts(A,self.cryst,\
-                F_func=self.get_deformation_gradient, coordinates='cart2D',method='regression',k=self.k)
+                F_func=self.get_deformation_gradient, coordinates='cart2D',method='regression',kI=self.kI,kII=self.kII)
             print('done!')
             #apply shifts
             self.crk.cauchy_born.apply_shifts(self.atoms,shifts,mask=self.shiftmask)
@@ -1010,7 +1135,8 @@ def set_atoms(self, atoms):
             self.alpha = atoms.info['alpha']
         else:
             self.alpha = 0.0
-        self.k = atoms.info['k']
+        self.kI = atoms.info['kI']
+        self.kII = atoms.info['kII']
         self.u[:] = np.zeros(3, self.N1)
         self.update_atoms()  # now we have same u_cle in atoms and self.atoms
         min_len = min(N1_in, self.N1)
@@ -1027,13 +1153,13 @@ def get_crack_tip_force(self, forces=None, mask=None, full_array_output=False):
         tip_y = self.cryst.cell.diagonal()[1] / 2.0
         dg = self.crk.deformation_gradient(self.cryst.positions[:, 0],
                                            self.cryst.positions[:, 1],
-                                           tip_x, tip_y, self.k)
+                                           tip_x, tip_y, self.kI, self.kII)
         V = np.zeros((len(self.cryst), 3))
         V[:, 0] = -(dg[:, 0, 0] - 1.0)
         V[:, 1] = -(dg[:, 0, 1])
         if self.crk.cauchy_born is not None:
             A = np.transpose(self.crk.RotationMatrix)
-            nu_grad = self.crk.cauchy_born.get_shift_gradients(A,self.cryst,F_func=self.get_deformation_gradient, coordinates='cart2D',k=self.k)
+            nu_grad = self.crk.cauchy_born.get_shift_gradients(A,self.cryst,F_func=self.get_deformation_gradient, coordinates='cart2D',kI=self.kI,kII=self.kII)
             if self.shiftmask is not None:
                 V[self.shiftmask]+=nu_grad[self.shiftmask]
             else:
@@ -1043,13 +1169,13 @@ def get_crack_tip_force(self, forces=None, mask=None, full_array_output=False):
         # V_fd = np.zeros((len(self.cryst), 3))
         # u, v = self.crk.displacements(self.cryst.positions[:, 0],
         #                               self.cryst.positions[:, 1],
-        #                               tip_x, tip_y, self.k)
+        #                               tip_x, tip_y, self.kI)
         # xp = self.cryst.positions[:, 0].copy()
         # for i in range(len(self.cryst)):
         #     xp[i] += eps
         #     up, vp = self.crk.displacements(xp,
         #                                     self.cryst.positions[:, 1],
-        #                                     tip_x, tip_y, self.k)
+        #                                     tip_x, tip_y, self.kI)
         #     xp[i] -= eps
         #     V_fd[i, 0] = - (up[i] - u[i]) / eps
         #     V_fd[i, 1] = - (vp[i] - v[i]) / eps
@@ -1114,9 +1240,14 @@ def get_forces(self, x1=None, xdot1=None, ds=None, forces=None, mask=None):
         if self.variable_alpha:
             f_alpha = self.get_crack_tip_force(forces, mask=mask)
             self.f_alpha_correction = 0#self.get_f_alpha_correction()
-            F.append((f_alpha-self.f_alpha_correction))
-            print('f_alpha',f_alpha)
-            #print('f_alpha_corrected',f_alpha-self.f_alpha_correction)
+            if self.is_3D:
+                f_beta = f_alpha*np.sin(self.theta_3D)
+                F.append(f_beta)
+                print('f_beta',f_beta)
+            else:
+                F.append((f_alpha-self.f_alpha_correction))
+                print('f_alpha',f_alpha)
+                #print('f_alpha_corrected',f_alpha-self.f_alpha_correction)
         if self.variable_k:
             f_k = self.get_k_force(x1, xdot1, ds)
             F.append(f_k)
@@ -1169,30 +1300,40 @@ def update_precon(self, x, F=None):
         Fu, Falpha, Fk = self.unpack(F)
         Pf_1 = spilu(P_1).solve(Fu)
         
+        #fix this for 3D
         if self.variable_alpha:
             dalpha = 0.001
             alph_before = self.alpha
             self.alpha += dalpha
-            self.update_atoms()
+            self.update_atoms(use_alpha_3D=True)
             F_after = self.get_forces(*self.precon_args)
             dF_dalpha = (F_after-F)/dalpha
+            if self.is_3D:
+                dF_dalpha *= np.sin(self.theta_3D)
             print('df_dalpha,',dF_dalpha)
             self.alpha = alph_before
-            self.update_atoms()
+            self.update_atoms(use_alpha_3D=True)
 
         print(self.variable_k)
         if self.variable_k:
             #currently this part isn't used and also doesn't really work
             dk = 0.0001
-            k_before = self.k
-            self.k += dk*self.k1g
+            if self.cont_k == 'k1':
+                k_before = self.kI
+                self.kI += dk*self.k1g
+            elif self.cont_k == 'k2':
+                k_before = self.kII
+                self.kII += dk*self.k2g
             self.update_atoms()
             F_after = self.get_forces(*self.precon_args)
             dF_dk = (F_after-F)/(dk*self.k1g)
             dF_dk[-1] = -self.dfk_dk_approx(self.precon_args[1])
             print(F_after[-1],F[-1],F_after[-1]-F[-1],dk*self.k1g)
             print('df_dk',dF_dk)
-            self.k = k_before
+            if self.cont_k == 'k1':
+                self.kI = k_before
+            elif self.cont_k == 'k2':
+                self.kII = k_before
             self.update_atoms()
 
         # extend diagonal of preconditioner for additional DoFs
@@ -1216,9 +1357,14 @@ def update_precon(self, x, F=None):
             I = np.append(I,offset)
             J = np.append(J,offset)
             Z = np.append(Z,-dF_dalpha[offset])
-            print('-dfdalpha',-dF_dalpha[offset])
+            if self.is_3D:
+                print('-dfdbeta',-dF_dalpha[offset])
+            else:
+                print('-dfdalpha',-dF_dalpha[offset])
             offset += 1
 
+
+        #TODO EITHER DELETE OR FIX THIS
         if self.variable_k:
             #add bottom edge of hessian, i = offset, j = 0:offset-1
             I = np.append(I,[offset for i in range(offset)])
@@ -1464,7 +1610,8 @@ def get_potential_energy(self):
         # print(f'E1={E1} E2={E2} total E={E1 + E2}')
         return E1 + E2
 
-    def rescale_k(self, new_k):
+    def rescale_k(self, new_kI):
+        # rescale_k, in the case of mode I fracture
         ref_x = self.cryst.positions[:, 0]
         ref_y = self.cryst.positions[:, 1]
         ref_z = self.cryst.positions[:, 2]
@@ -1473,11 +1620,11 @@ def rescale_k(self, new_k):
         x, y, z = self.atoms.get_positions().T
 
         # rescale full displacement field (CLE + atomistic corrector)
-        x = ref_x + new_k / self.k * (x - ref_x)
-        y = ref_y + new_k / self.k * (y - ref_y)
+        x = ref_x + new_kI / self.kI * (x - ref_x)
+        y = ref_y + new_kI / self.kI * (y - ref_y)
 
-        self.k = new_k
-        u_cle = new_k * self.u_cle()  # CLE solution at new K
+        self.kI = new_kI
+        u_cle = self.u_cle(kI=new_kI)  # CLE solution at new K
         self.u[:] = np.c_[x - u_cle[:, 0] - ref_x,
                           y - u_cle[:, 1] - ref_y,
                           z - ref_z][self.regionI, :]
@@ -1493,10 +1640,42 @@ def arc_length_continuation(self, x0, x1, N=10, ds=0.01, ftol=1e-2,
                                 cos_alpha_min=0.9,parallel=False,
                                 pipe_output=None,data_queue=None,
                                 kill_confirm_queue=None,
-                                allow_alpha_backtracking=False):
+                                allow_alpha_backtracking=False,
+                                follow_G_contour=False):
         import h5py
         assert self.variable_k  # only makes sense if K can vary
 
+        #check to see if x0 and x1 contain kI and kII, if they do
+        #then remove whichever one isn't being used for continuation
+        if len(self)+1 == len(x0):
+            if self.cont_k == 'k1':
+                idx = -1
+                self.kII = x0[idx]
+            elif self.cont_k == 'k2':
+                idx = -2
+                self.kI = x0[idx]
+            x0 = np.delete(x0,idx)
+
+        if len(self)+1 == len(x1):
+            if self.cont_k == 'k1':
+                idx = -1
+                self.kII = x1[idx]
+            elif self.cont_k == 'k2':
+                idx = -2
+                self.kI = x1[idx]
+            x1 = np.delete(x1,idx)
+                
+        #if following G contour, get the initial values of kI and kII
+        self.follow_G_contour = follow_G_contour
+        if self.follow_G_contour:
+            if self.cont_k == 'k1':
+                self.kI0 = x0[-1]
+                self.kII0 = self.kII
+            elif self.cont_k == 'k2':
+                self.kII0 = x0[-1]
+                self.kI0 = self.kI
+            print('og vals', self.kI0, self.kII0)
+
         if continuation:
             xdot1 = self.get_xdot(x0, x1, ds)
         else:
@@ -1517,14 +1696,19 @@ def arc_length_continuation(self, x0, x1, N=10, ds=0.01, ftol=1e-2,
                         if 'x' in hf.keys():
                             x_traj = hf['x']
                         else:
-                            x_traj = hf.create_dataset('x', (0, len(self)),
-                                                    maxshape=(None, len(self)),
+                            if self.variable_k:
+                                dof_length = len(self) + 1 #write both kI and kII to file
+                            else:
+                                dof_length = len(self)
+                            x_traj = hf.create_dataset('x', (0, dof_length),
+                                                    maxshape=(None, dof_length),
                                                     compression='gzip')
                             x_traj.attrs['ds'] = ds
                             x_traj.attrs['ftol'] = ftol
                             x_traj.attrs['direction'] = direction
                             x_traj.attrs['traj_interval'] = traj_interval
                             #row = x_traj.shape[0]
+                    break
                 except:
                     print('failed to access file 1')
                     time.sleep(1)
@@ -1579,7 +1763,9 @@ def arc_length_continuation(self, x0, x1, N=10, ds=0.01, ftol=1e-2,
                             with h5py.File(traj_file, 'a') as hf:
                                 x_traj = hf['x']
                                 x_traj.resize((x_traj.shape[0] + 1, x_traj.shape[1]))
-                                x_traj[-1, :] = x2
+                                print(f'writing k1 {self.kI}, k2 {self.kII}')
+                                #print(np.append(x2[:-1],[self.kI,self.kII]))
+                                x_traj[-1, :] = np.append(x2[:-1],[self.kI,self.kII])
                                 row += 1
                                 break
                         except:
@@ -1683,10 +1869,10 @@ def animate(self, x, k1g, regions='12', rzoom=np.inf, cutoff=2.8, frames=None, c
         region = a.arrays['region']
 
         i = 0
-        ax1.plot(x[:, -2], x[:, -1] / k1g, 'b-')
-        (blob,) = ax1.plot([x[i, -2]], [x[i, -1] / k1g], 'rx', mew=5, ms=20)
+        ax1.plot(x[:, -3], x[:, -2] / k1g, 'b-')
+        (blob,) = ax1.plot([x[i, -3]], [x[i, -2] / k1g], 'rx', mew=5, ms=20)
         ax1.set_xlabel(r'Crack position $\alpha$')
-        ax1.set_ylabel(r'Stress intensity factor $K/K_{G}$');
+        ax1.set_ylabel(r'Stress intensity factor $KI/KI_{G}$');
 
         self.set_dofs(x[i, :])
         if cutoff is None:
@@ -1786,7 +1972,7 @@ def write_atoms_to_file(self, fname):
 
             #get strains applied
             A = np.transpose(self.crk.RotationMatrix)
-            E,R = self.crk.cauchy_born.evaluate_F_or_E(A,self.cryst,F_func=self.get_deformation_gradient, coordinates='cart2D',k=self.k)
+            E,R = self.crk.cauchy_born.evaluate_F_or_E(A,self.cryst,F_func=self.get_deformation_gradient, coordinates='cart2D',k=self.kI)
             for i in range(len(crack_atoms)-1):
                 E[i,:,:] = np.transpose(A)@E[i,:,:]@A #transform back to lab frame
             #store values of E
@@ -2030,6 +2216,55 @@ def isotropic_modeI_crack_tip_displacement_field(K, G, nu, r, t,
 
     return u, v
 
+def isotropic_modeII_crack_tip_displacement_field(K, G, nu, r, t,
+                                                 stress_state=PLANE_STRAIN):
+    """
+    Compute Irwin singular crack tip displacement field for mode II fracture.
+
+    Parameters
+    ----------
+    K : float
+        Mode I stress intensity factor. Units should match units of `G` and `r`.
+    G : float
+        Shear modulus. Units should match units of `K` and `r`.
+    nu : float
+        Poisson ratio.
+    r : array_like
+        Radial distances from crack tip. Can be a multidimensional
+        array to evaluate stress field on a grid.
+    t : array_like
+        Angles from horizontal line y=0 ahead of crack tip,
+        measured anticlockwise. Should have same shape as `r`.
+    stress_state : str
+        One of"plane stress" or "plane strain". Used if xyz_only=False to
+        determine zz stresses.
+
+    Returns
+    -------
+    u : array
+    v : array
+        Displacements. Same shape as `r` and `t`.
+    """
+
+    if r.shape != t.shape:
+        raise ValueError('Shapes of radial and angular arrays "r" and "t" '
+                         'must match.')
+
+    if stress_state == PLANE_STRAIN:
+        kappa = 3-4*nu
+    elif stress_state == PLANE_STRESS:
+        kappa = (3.-nu)/(1.+nu)
+    else:
+        raise ValueError('"stress_state" should be either "{0}" or "{1}".'
+            .format(PLANE_STRAIN, PLANE_STRESS))
+
+    #use the lawn book form for now, as I do not have the form for mode II that matches above to hand
+    # Form in Lawn book is equivalent:
+    radial = (K/(4*G))*np.sqrt(r/(2.*math.pi))
+    u = radial*((2*kappa + 3)*np.sin(t/2) + np.sin(3*t/2))
+    v = -radial*((2*kappa - 3)*np.cos(t/2) + np.cos(3*t/2))
+
+    return u, v
 
 class IsotropicStressField(object):
     """
diff --git a/matscipy/surface_reconstruction.py b/matscipy/surface_reconstruction.py
index 56fe8dae..e3b7f5cd 100644
--- a/matscipy/surface_reconstruction.py
+++ b/matscipy/surface_reconstruction.py
@@ -62,7 +62,7 @@ def __init__(self, el, a0, calc, directions, surf_dir, lattice='diamond', multil
                 self.el, self.a0, self.calc, lattice=self.lattice)
             self.multilattice = True
 
-    def map_surface(self, fmax=0.0001, layers=6, cutoff=10, shift=0, switch=False):
+    def map_surface(self, fmax=0.0001, layers=6, cutoff=10, shift=0, switch=False,invert_dirs=[False,False,False]):
         """Map the relaxation of the crystal surface a certain number of layers deep.
         Parameters
         ----------
@@ -82,6 +82,11 @@ def map_surface(self, fmax=0.0001, layers=6, cutoff=10, shift=0, switch=False):
         switch : bool
             Whether or not to switch the sublattices around in the mapped structure.
             (Essentially, the correct surface can always be found from some combination of shift and switch)
+        invert_dirs : bool array
+            Whether to invert the components of the mapped shifts. This is useful in the case where a surface which originally has mirror
+            symmetry in a certain plane can break this symmetry during reconstruction. This is the case for the 110 silicon reconstructed surface,
+            where the surface breaks mirror symmetry in the {001} type plane during reconstruction. Two possible surfaces can form, related through
+            diad symmetry. To map the other surface, one has to switch sublattices and invert the mapped displacements along the <001> type direction.
         """
         # build a single cell to determine how many atomic surface
         # layers there are in a single unit cell made by ASE
@@ -102,10 +107,10 @@ def map_surface(self, fmax=0.0001, layers=6, cutoff=10, shift=0, switch=False):
             #check if multilattice atoms are coplanar
             diff = single_cell.get_positions()[self.cb.lattice1mask][0, self.surf_dir] - single_cell.get_positions()[self.cb.lattice2mask][0, self.surf_dir]
             print(diff)
-            if np.abs(diff)<0.001:
-                self.split_cb_pair = False
-            else:
-                self.split_cb_pair = True
+            #if np.abs(diff)<0.001:
+            #    self.split_cb_pair = False
+            #else:
+            #    self.split_cb_pair = True
         else:
             # otherwise, just look at all the layers directly
             n_layer_per_cell = len(np.unique(single_cell.get_positions()[
@@ -165,6 +170,11 @@ def map_surface(self, fmax=0.0001, layers=6, cutoff=10, shift=0, switch=False):
         print('POS SHIFT',pos_shift)
         # print('fafter',slab.get_forces())
 
+        #invert measured shifts if required
+        for direc in range(len(invert_dirs)):
+            if invert_dirs[direc]:
+                pos_shift[:,direc] *= -1
+        
         # for a multi-lattice:
         if self.multilattice:
             # if this is a multi-lattice
@@ -472,7 +482,7 @@ def map_pandey_111(self, fmax=0.0001, layers=6, cutoff=10,orientation=0,shift=0,
         # map out the pandey reconstruction on a diamond structure 111 surface
         # can then apply it later (this requires a seperate function as the unit cell of 2x1 reconstruction
         # breaks the symmetry of the surface)
-        self.split_cb_pair = True
+        #self.split_cb_pair = True
         #111 type surface always splits a cb pair
         ny = 1
         nx = 1
@@ -491,14 +501,14 @@ def map_pandey_111(self, fmax=0.0001, layers=6, cutoff=10,orientation=0,shift=0,
                     directions=self.pandey_dirs
                     )
 
-
         self.cb.set_sublattices(a,self.U)
         if switch:
             self.cb.switch_sublattices(a)
         sorted_z_vals = np.sort(a.get_positions()[self.cb.lattice1mask][:,2])
         self.inter_surf_dist = sorted_z_vals[2] - sorted_z_vals[0]
-        a.translate([0,0,-shift*self.inter_surf_dist])
-
+        #a.translate([0,+0.001,-shift*self.inter_surf_dist]) #REMOVE SHIFT FROM HERE
+        a.translate([0.01,0.01,0]) 
+        a.wrap()
         print('INTER SURF DIST',self.inter_surf_dist)
         #seed the 2x1 reconstruction on the top plane
         sx, sy, sz = a.get_cell().diagonal()
@@ -517,14 +527,16 @@ def map_pandey_111(self, fmax=0.0001, layers=6, cutoff=10,orientation=0,shift=0,
                                                     np.logical_not(mask))]
         a.set_distance(top1, top2, bondlen)
         a.set_distance(topA, topB, bondlen)
-
+        #ase.io.write('intermediate_atoms.xyz',a)
         x, y, z = a.positions.T
         y1 = (y[top1]+y[top2])/2
         yA = (y[topA]+y[topB])/2
         y[top1] += yA-y1+a.cell[1,1]/2
         y[top2] += yA-y1+a.cell[1,1]/2
+        #x[top1] += a.cell[0,0]/2
+        #x[topB] -= a.cell[0,0]/2
         a.set_positions(np.transpose([x, y, z]))
-        a.wrap()
+        #a.wrap()
 
         ase.io.write('0_bulk.xyz',bulk)
         ase.io.write('1_bulk.xyz',a)
@@ -547,6 +559,20 @@ def map_pandey_111(self, fmax=0.0001, layers=6, cutoff=10,orientation=0,shift=0,
         opt_a.run(fmax=fmax)
         ase.io.write('2_bulk.xyz',a)
         pos_after_unrotated = a.get_positions()
+        print(pos_after_unrotated[:,1]<0)
+        print(len(pos_after_unrotated[:,1]<0))
+        print(pos_after_unrotated[:,1]>a.cell[1,1])
+        print(len(pos_after_unrotated[:,1]>a.cell[1,1]))
+        atoms_outside_mask = (pos_after_unrotated[:,1]<0) | (pos_after_unrotated[:,1]>a.cell[1,1])
+        pushed_out_atoms = a[atoms_outside_mask]
+        wrapped=False
+        if len(pushed_out_atoms)>0:
+            print('detected atoms to be wrapped')
+            wrapped=True
+            a.wrap()
+            pos_after_unrotated = a.get_positions()
+            ase.io.write('3_bulk.xyz',a)
+        
         pos_diff_unrotated = pos_after_unrotated-pos_initial_unrotated
         #print('pos diff unrotated',pos_diff_unrotated)
         #map out the layers in the original structure
@@ -557,6 +583,69 @@ def map_pandey_111(self, fmax=0.0001, layers=6, cutoff=10,orientation=0,shift=0,
             self.identify_pandey_layers(bulk,surface_coords,read_from_atoms=True,frame_dirs=self.pandey_dirs)
         self.surf_dir = tmp
 
+        """
+        #The following is code written in order to minimise the amount that atoms
+        move during the relaxation- this is not necessary when there is no chance 
+        that the position of xlim always lies at the edge of a pandey unit cell
+        #also, i think actually that the atoms do end up in positions that minimise the travel
+        distance anyway
+
+        a.new_array('atomno',np.zeros(len(a),dtype=int))
+        a_atomno = a.arrays['atomno']
+        a_atomno[layer_mask_set_lattice1_atom1[:,0]] = 1
+        a_atomno[layer_mask_set_lattice1_atom2[:,0]] = 2
+        a_atomno[layer_mask_set_lattice2_atom1[:,0]] = 3
+        a_atomno[layer_mask_set_lattice2_atom2[:,0]] = 4
+        ase.io.write('4_bulk.xyz',a)
+
+        if wrapped:
+            print('DETECTED WRAPPED AS TRUE')
+            sublattice1_layer1_atom1_index = np.where(layer_mask_set_lattice1_atom1[:,0])[0]
+            sublattice1_layer1_atom2_index = np.where(layer_mask_set_lattice1_atom2[:,0])[0]
+            sublattice2_layer1_atom1_index = np.where(layer_mask_set_lattice2_atom1[:,0])[0]
+            sublattice2_layer1_atom2_index = np.where(layer_mask_set_lattice2_atom2[:,0])[0]
+            print('sublattice1_layer1_atom1_index',sublattice1_layer1_atom1_index)
+            print('sublattice1_layer1_atom2_index',sublattice1_layer1_atom2_index)
+            print('sublattice2_layer1_atom1_index',sublattice2_layer1_atom1_index)
+            print('sublattice2_layer1_atom2_index',sublattice2_layer1_atom2_index)
+
+
+            comparison_sublattice_1_list = [sublattice1_layer1_atom1_index,sublattice1_layer1_atom2_index]
+            comparison_sublattice_2_list = [sublattice2_layer1_atom1_index,sublattice2_layer1_atom2_index]
+            print('comparison sublattice 1 list',comparison_sublattice_1_list)
+            print('comparison sublattice 2 list',comparison_sublattice_2_list)
+
+            index_switch_list = []
+            for i in range(2):
+                print(i)
+                dist_sublattice1 = []
+                dist_sublattice2 = []
+                #this picks a single atom on each sublattice in the wrapped relaxed structure and compares it to each of the two
+                #atoms of the same sublattice in the initial structure.
+                for j in range(2):
+                    dist_sublattice1.append(np.abs(pos_after_unrotated[comparison_sublattice_1_list[i],1]-pos_initial_unrotated[comparison_sublattice_1_list[j],1]))
+                    dist_sublattice2.append(np.abs(pos_after_unrotated[comparison_sublattice_2_list[i],1]-pos_initial_unrotated[comparison_sublattice_2_list[j],1]))
+                print('dist sublattice 1',dist_sublattice1)
+                print('dist sublattice 2',dist_sublattice2)
+                #the next step is to check which of the two distances is smaller.
+                #if it is the one where i=j, then there is no problem and nothing needs to change
+                #if it is the one where i!=j, then the atom index needs to be switched - such that
+                #the atom in the RELAXED structure that was indexed by i is now indexed by j
+                min_ind_sublattice1 = np.argmin(np.abs(dist_sublattice1))
+                min_ind_sublattice2 = np.argmin(np.abs(dist_sublattice2))
+                if min_ind_sublattice1 != i:
+                    index_switch_list.append([comparison_sublattice_1_list[i],comparison_sublattice_1_list[min_ind_sublattice1]])
+
+            a_before = a.copy()
+            atomno_copy = a_atomno.copy()
+            print(a)
+            print('index swtich list',index_switch_list)
+            for i_list in index_switch_list:
+                a[i_list[0][0]].position = a_before[i_list[1][0]].position.copy()
+                a[i_list[1][0]].position = a_before[i_list[0][0]].position.copy()
+            ase.io.write('5_bulk.xyz',a)
+        """
+
         #relaxation array lattice 1 atom 1
         self.ral1a1 = np.zeros([layers,3])
         #relaxation array lattice 1 atom 2
@@ -582,15 +671,16 @@ def map_pandey_111(self, fmax=0.0001, layers=6, cutoff=10,orientation=0,shift=0,
         #print('with no rotation',np.transpose(self.A)@self.U)
         #print('with rotation',R)
         for layer in range(layers):
+            print(layer)
+            print(len(pos_diff_unrotated[layer_mask_set_lattice2_atom2[:,layer]]))
+
             self.ral1a1[layer,:] = R@pos_diff_unrotated[layer_mask_set_lattice1_atom1[:,layer]][0,:]
             self.ral1a2[layer,:] = R@pos_diff_unrotated[layer_mask_set_lattice1_atom2[:,layer]][0,:]
             self.ral2a1[layer,:] = R@pos_diff_unrotated[layer_mask_set_lattice2_atom1[:,layer]][0,:]
             self.ral2a2[layer,:] = R@pos_diff_unrotated[layer_mask_set_lattice2_atom2[:,layer]][0,:]
-        
         print(self.ral1a1,self.ral1a2,self.ral2a1,self.ral2a2)
 
 
-
     def identify_pandey_layers(self,
             atoms,
             surface_coords,
@@ -677,7 +767,8 @@ def identify_pandey_layers(self,
         OR[1,1] = np.cos(2*(np.pi/3)*orientation)
         OR[2,2] = 1
         
-        #T - lab frame to lattice fram
+        #OR - Orientation rotation tensor
+        #T - lab frame to lattice frame
         #U - pandey map frame to lattice frame
         R = OR@np.transpose(self.U)@T
         print(self.U,T,R)
@@ -687,29 +778,36 @@ def identify_pandey_layers(self,
         layer_mask_set_lattice2_atom1 = np.zeros_like(layer_mask_set_lattice1,dtype=bool)
         layer_mask_set_lattice2_atom2 = np.zeros_like(layer_mask_set_lattice1,dtype=bool)
         
-        #get the first layer 1 atom in the direction we need
-        layer1 = atoms[layer_mask_set_lattice1[:,0]]
-        layer1pos = layer1.get_positions()
+        #find an exclusion zone using the position of the leftmost lattice1
+        #atom in the third layer (this is where we cut the cell in map_pandey)
+        layer3 = atoms[layer_mask_set_lattice1[:,2]]
+        layer3pos = layer3.get_positions()
         #rotate
-        rotated_layer1 = np.zeros_like(layer1pos)
-        for j in range(np.shape(layer1pos)[0]):
-            rotated_layer1[j,:] = R@layer1pos[j,:]
+        rotated_layer3 = np.zeros_like(layer3pos)
+        for j in range(np.shape(layer3pos)[0]):
+            rotated_layer3[j,:] = R@layer3pos[j,:]
         
-        leftpos = np.min(rotated_layer1[:,1])
-        print('LEFTPOS',leftpos)
-        print('R',R)
-        shiftvec = np.transpose(R)@np.array([0,-leftpos,0])
-        #shift all positions left, add any negative ones back to the end
-        #(to make sure that labelling is consistent with mapping)
-        ase.io.write('1_atoms_shift.xyz',atoms)
-        #atoms.translate(shiftvec)
-        #atoms.wrap()
-        ase.io.write('2_atoms_shift.xyz',atoms)
-        #print('SHIFTVEC',shiftvec)
-
-        #we need to shift the whole unitcell so the left most 
-        print('SEARCH DIR,',search_dir)
+        #find the leftmost point of layer 3
+        leftpos = np.min(rotated_layer3[:,1])
+        pos = atoms.get_positions()
+        rotated_pos = np.zeros_like(pos)
+        for j in range(np.shape(pos)[0]):
+            rotated_pos[j,:] = R@pos[j,:]
+        
+        mask = rotated_pos[:,1]<leftpos-0.01
+        print(np.where(mask))
+        for i in range(np.shape(layer_mask_set_lattice1)[1]):
+            print(np.where(layer_mask_set_lattice1[:,i]))
+            print(np.where(layer_mask_set_lattice1[:,i]^(mask&layer_mask_set_lattice1[:,i])))
+            print('lattice2')
+            print(np.where(layer_mask_set_lattice2[:,i]))
+            print(np.where(layer_mask_set_lattice2[:,i]^(mask&layer_mask_set_lattice2[:,i])))
+
+            #TODO uncomment these lines? why did I even comment them? this really needs testing
+            #layer_mask_set_lattice1[:,i] = layer_mask_set_lattice1[:,i]^(mask&layer_mask_set_lattice1[:,i])
+            #layer_mask_set_lattice2[:,i] = layer_mask_set_lattice2[:,i]^(mask&layer_mask_set_lattice2[:,i])
 
+        print('SEARCH DIR,',search_dir)
         for i in range(np.shape(layer_mask_set_lattice1)[1]):
             #print(np.shape(layer_mask_set_lattice1))
             lattice1atoms = atoms[layer_mask_set_lattice1[:,i]]
@@ -790,7 +888,7 @@ def identify_pandey_layers(self,
         #print(layer_mask_set_lattice1_atom1[:,layer])
         try:
             atoms.new_array('atomno',np.zeros(len(atoms),dtype=int))
-        except:
+        except: #TODO add correct exception here (for the case where the array already exists)
             pass
         atomno = atoms.arrays['atomno']
         for layer in range(self.layers):
@@ -898,4 +996,494 @@ def apply_pandey_111(self,
 
         atoms.wrap()
 
+    def map_si_110_3x1(self, fmax=0.0001, layers=6, cutoff=10,switch=False,orientation=0,shift=0,permute=False):
+        """Map the silicon 110 3x1 relaxation of the crystal.
+        ------------------------------
+        fmax : float
+            Force tolerance for relaxation
+        layers : int
+            Number of layers deep to map the free surface
+        cutoff : float
+            Cutoff of the potential being used - used to decide the number of
+            atoms beneath the mapped layers that need to be modelled 
+        shift : float
+            Amount to shift the bulk structure before the surface is mapped
+            such that the top layer of surface is physically meaningful. One
+            use case of this would be to adjust the top layer of atoms being mapped
+            in diamond such that it matches the top layer of atoms seen on the cleavage
+            plane.
+        switch : bool
+            Whether or not to switch the sublattices around in the mapped structure.
+            (Essentially, the correct surface can always be found from some combination of shift and switch)
+        orientation : int, 1 or 0
+            Whether to flip the reconstruction by 180 degrees. This is necessary for the lower surface.
+        """
+        # map out the pandey reconstruction on a diamond structure 111 surface
+        # can then apply it later (this requires a seperate function as the unit cell of 3x1 reconstruction
+        # breaks the symmetry of the surface)
+        #self.split_cb_pair = True
+        #111 type surface always splits a cb pair
+        ny = 3
+        nx = 1
+        nz = 3*int(layers)
+        self.layers = layers
+        self.si_110_3x1_dirs = [ [0,0,1], [1,-1,0], [1,1,0] ]
+        #build the reconstruction frame rotation matrix
+        U = np.zeros([3, 3])
+        for i, direction in enumerate(self.si_110_3x1_dirs):
+            direction = np.array(direction)
+            U[:, i] = direction / np.linalg.norm(direction)
+        self.U = U  # reconstruction frame rotation matrix
+        a = Diamond(self.el,
+                    size             = [1, ny, nz],
+                    latticeconstant  = self.a0,
+                    directions=self.si_110_3x1_dirs
+                    )
+
+        self.cb.set_sublattices(a,self.U)
+        if switch:
+            self.cb.switch_sublattices(a)
+        sorted_z_vals = np.sort(a.get_positions()[self.cb.lattice1mask][:,2])
+        self.inter_surf_dist = sorted_z_vals[3] - sorted_z_vals[0]
+        a.translate([-0.01,0.0,-shift*self.inter_surf_dist]) #REMOVE SHIFT FROM HERE
+        #a.translate([0.01,0.01,-1]) 
+        a.wrap()
+        print('INTER SURF DIST',self.inter_surf_dist)
+        #seed the 3x1 reconstruction on the top plane
+        pos_initial_unrotated = a.get_positions()
+
+        bulk = a.copy()
+        top_atom_mask = a.positions[:,2]>(np.max(a.positions[:,2])-self.inter_surf_dist)
+        top_atom_indices = np.where(top_atom_mask)[0]
+        sorted_indices = top_atom_indices[np.argsort(a.positions[top_atom_indices,1])]
+        
+        if permute:
+            a.positions[sorted_indices[2],2] += -1
+            a.positions[sorted_indices[3],2] += -1
+            a.positions[sorted_indices[5],2] += -1
+            a.positions[sorted_indices[0],2] += -1
+        else:
+            a.positions[sorted_indices[1],2] += -1
+            a.positions[sorted_indices[2],2] += -1
+            a.positions[sorted_indices[4],2] += -1
+            a.positions[sorted_indices[5],2] += -1
+
+
+
+
+        ase.io.write('0_bulk.xyz',bulk)
+        ase.io.write('1_bulk.xyz',a)
+        #relax the structure to get the 3x1 111 reconstruction
+        # run an optimisation to relax the surface
+
+        #add some vacuum
+        cell = a.get_cell()
+        # vacuum along self.surf_dir axis (surface normal)
+        cell[2, :] *= 2
+        a.set_cell(cell, scale_atoms=False)
+        a.calc = self.calc
+
+        # fix bottom atoms to stop lower surface reconstruction
+        #print('POS INITIAL UNROTATED',pos_initial_unrotated)
+        mask = pos_initial_unrotated[:,2] < cutoff
+        #print(pos_initial_unrotated[:,2])
+        a.set_constraint(FixAtoms(mask=mask))
+        opt_a = PreconLBFGS(a)
+        opt_a.run(fmax=fmax)
+        ase.io.write('2_bulk.xyz',a)
+        pos_after_unrotated = a.get_positions()
+        print(pos_after_unrotated[:,1]<0)
+        print(len(pos_after_unrotated[:,1]<0))
+        print(pos_after_unrotated[:,1]>a.cell[1,1])
+        print(len(pos_after_unrotated[:,1]>a.cell[1,1]))
+        atoms_outside_mask = (pos_after_unrotated[:,1]<0) | (pos_after_unrotated[:,1]>a.cell[1,1])
+        pushed_out_atoms = a[atoms_outside_mask]
+        wrapped=False
+        if len(pushed_out_atoms)>0:
+            print('detected atoms to be wrapped')
+            wrapped=True
+            #a.wrap()
+            pos_after_unrotated = a.get_positions()
+            ase.io.write('3_bulk.xyz',a)
+        
+        pos_diff_unrotated = pos_after_unrotated-pos_initial_unrotated
+        #print('pos diff unrotated',pos_diff_unrotated)
+        #map out the layers in the original structure
+        surface_coords = pos_initial_unrotated[np.argmax(pos_initial_unrotated[:,2]),:]+[0,0,0.01]
+        tmp = self.surf_dir
+        self.surf_dir = 2
+        layer_mask_set_lattice1_atom1, layer_mask_set_lattice1_atom2, layer_mask_set_lattice1_atom3, layer_mask_set_lattice2_atom1, layer_mask_set_lattice2_atom2,\
+            layer_mask_set_lattice2_atom3 = self.identify_si_110_layers(bulk,surface_coords,read_from_atoms=True,frame_dirs=self.si_110_3x1_dirs,switch=switch)
+        self.surf_dir = tmp
+
+        #relaxation array lattice 1 atom 1
+        self.ral1a1 = np.zeros([layers,3])
+        #relaxation array lattice 1 atom 2
+        self.ral1a2 = np.zeros([layers,3])
+        #relaxation array lattice 1 atom 3
+        self.ral1a3 = np.zeros([layers,3])
+        #relaxation array lattice 2 atom 1
+        self.ral2a1 = np.zeros([layers,3])
+        #relaxation array lattice 2 atom 2
+        self.ral2a2 = np.zeros([layers,3])
+        #relaxation array lattice 2 atom 3
+        self.ral2a3 = np.zeros([layers,3])
+
+        #print(np.shape(pos_diff_unrotated))
+        #print(np.shape(layer_mask_set_lattice1_atom1))
+        #generate rotation matrix to get from pandey frame to lab frame 
+        #orientation rotation tensor
+        OR = np.zeros([3,3])
+        OR[0,0] = np.cos((np.pi)*orientation)
+        OR[0,1] = -np.sin((np.pi)*orientation)
+        OR[1,0] = np.sin((np.pi)*orientation)
+        OR[1,1] = np.cos((np.pi)*orientation)
+        OR[2,2] = 1
+        #print('OR',OR)
+
+        R = np.transpose(self.A)@self.U@np.transpose(OR)
+        #print('with no rotation',np.transpose(self.A)@self.U)
+        #print('with rotation',R)
+        for layer in range(layers):
+            print(layer)
+            print(len(pos_diff_unrotated[layer_mask_set_lattice2_atom2[:,layer]]))
+
+            self.ral1a1[layer,:] = R@pos_diff_unrotated[layer_mask_set_lattice1_atom1[:,layer]][0,:]
+            self.ral1a2[layer,:] = R@pos_diff_unrotated[layer_mask_set_lattice1_atom2[:,layer]][0,:]
+            self.ral1a3[layer,:] = R@pos_diff_unrotated[layer_mask_set_lattice1_atom3[:,layer]][0,:]
+            self.ral2a1[layer,:] = R@pos_diff_unrotated[layer_mask_set_lattice2_atom1[:,layer]][0,:]
+            self.ral2a2[layer,:] = R@pos_diff_unrotated[layer_mask_set_lattice2_atom2[:,layer]][0,:]
+            self.ral2a3[layer,:] = R@pos_diff_unrotated[layer_mask_set_lattice2_atom3[:,layer]][0,:]
+        print(self.ral1a1,self.ral1a2,self.ral2a1,self.ral2a2)
+
+
+    def identify_si_110_layers(self,
+            atoms,
+            surface_coords,
+            xlim=None,
+            ylim=None,
+            zlim=None,
+            search_dir=-1,
+            atoms_for_cb=None,
+            read_from_atoms=False,
+            frame_dirs=None,
+            orientation=0,
+            switch=False):
+        """Function for identifying the layers adjacent to a surface, in a supplied atoms structure.
+        This allows a mapped relaxation to be easily applied to a surface. This implementation is specific
+        to the 111 Pandey reconstruction, as more atoms need to be identified due to the breaking of surface
+        symmetry.
+
+        Parameters
+        ----------
+        atoms : ASE atoms object
+            Atoms object which contains surface to relax
+        surface_coords : array_like
+            Coordinates of a point that lies on the free surface (or within one layer distance above it)
+        xlim : array_like
+            [x lower, x upper], x range to apply relaxation
+        ylim : array_like
+            [y lower, y upper], y range to apply relaxation
+        zlim : array_like
+            [z lower, z upper], z range to apply relaxation
+        search_dir : int
+            -1: surface is below point provided, 1: surface is above point provided
+        atoms_for_cb : ASE atoms object
+            atoms object which allows the easy identification of sub-lattices in a multi-lattice
+            crystal using the CubicCauchyBorn set_sublattices function. The difference between this and atoms may
+            be periodic boundary conditions, for example.
+        read_from_atoms : bool
+            Whether or not to read in the cauchy-born sublattices from the atoms structure. This can be done if the
+            atoms object contains arrays called lattice1mask and lattice2mask
+        frame_dirs : 
+            The frame in which the Pandey reconstruction was mapped out - the atoms in the supplied frame will need to
+            be rotated to this frame before the reconstruction can be applied
+        orientation : 0,1 or 2
+            Which of the three different orientations of the Pandey reconstructions are being applied to the surface
+
+
+        Returns
+        -------
+        layer_mask_set : array_like, bool
+            A 2D array dimensions [nlayers,natoms] where a row i is a mask for
+            'atoms' corresponding to a layer i deep from the surface. In the case of
+            a multi-lattice, two masks are returned, one for each of the sublattices on each layer.
+        """
+        # this function needs to identify the 4 sets of atoms per layer
+        # that reconstruct (breaking symmetry)
+        # it does this by taking an atoms structure, reading the sublattices
+        # identifying the layers down from a point and then in each layer
+        # identifying the 4 unique types of atom
+        # to do this, it takes the positions of all the atoms of one sublattice
+        # rotates these position vectors from the supplied frame to the mapped reconstruction frame
+        # and then goes from left to right, shifting th
+        
+        #get the layer masks
+
+        layer_mask_set_lattice1, layer_mask_set_lattice2 = \
+            self.identify_layers(atoms,surface_coords,read_from_atoms=read_from_atoms,xlim=xlim,ylim=ylim,zlim=zlim,atoms_for_cb=atoms_for_cb,search_dir=search_dir)
+
+        #rotate the lattice structure to match the pandey directions
+        #build the rotation matrix
+        #this is equivalent to transforming from the coordinate system
+        #given by 'dirs' (the lab frame) to the lattice frame, and then from there
+        #to the frame where the pandey reconstruction is defined.
+
+        if frame_dirs is not None:
+            T = np.zeros([3, 3])
+            for i, direction in enumerate(frame_dirs):
+                direction = np.array(direction)
+                T[:, i] = direction / np.linalg.norm(direction)
+        else:
+            T = self.A
+
+        OR = np.zeros([3,3])
+        OR[0,0] = np.cos((np.pi)*orientation)
+        OR[0,1] = -np.sin((np.pi)*orientation)
+        OR[1,0] = np.sin((np.pi)*orientation)
+        OR[1,1] = np.cos((np.pi)*orientation)
+        OR[2,2] = 1
+        
+        #OR - Orientation rotation tensor
+        #T - lab frame to lattice frame
+        #U - lattice map frame to lattice frame
+        R = OR@np.transpose(self.U)@T
+        print(self.U,T,R)
+
+        layer_mask_set_lattice1_atom1 = np.zeros_like(layer_mask_set_lattice1,dtype=bool)
+        layer_mask_set_lattice1_atom2 = np.zeros_like(layer_mask_set_lattice1,dtype=bool)
+        layer_mask_set_lattice1_atom3 = np.zeros_like(layer_mask_set_lattice1,dtype=bool)
+        layer_mask_set_lattice2_atom1 = np.zeros_like(layer_mask_set_lattice1,dtype=bool)
+        layer_mask_set_lattice2_atom2 = np.zeros_like(layer_mask_set_lattice1,dtype=bool)
+        layer_mask_set_lattice2_atom3 = np.zeros_like(layer_mask_set_lattice1,dtype=bool)
+        
+        #find an exclusion zone using the position of the leftmost lattice1
+        #atom in the first layer (this is where we cut the cell in map_pandey)
+        if switch:
+            layer1 = atoms[layer_mask_set_lattice2[:,0]]
+        else:
+            layer1 = atoms[layer_mask_set_lattice1[:,0]]
+        
+        layer1pos = layer1.get_positions()
+        #rotate
+        rotated_layer1 = np.zeros_like(layer1pos)
+        for j in range(np.shape(layer1pos)[0]):
+            rotated_layer1[j,:] = R@layer1pos[j,:]
+        
+        #find the leftmost point of layer 1
+        leftpos = np.min(rotated_layer1[:,1])
+        pos = atoms.get_positions()
+        rotated_pos = np.zeros_like(pos)
+        for j in range(np.shape(pos)[0]):
+            rotated_pos[j,:] = R@pos[j,:]
+        
+        mask = rotated_pos[:,1]<leftpos-0.01
+        print(np.where(mask))
+        for i in range(np.shape(layer_mask_set_lattice1)[1]):
+            print(np.where(layer_mask_set_lattice1[:,i]))
+            print(np.where(layer_mask_set_lattice1[:,i]^(mask&layer_mask_set_lattice1[:,i])))
+            print('lattice2')
+            print(np.where(layer_mask_set_lattice2[:,i]))
+            print(np.where(layer_mask_set_lattice2[:,i]^(mask&layer_mask_set_lattice2[:,i])))
+            #layer_mask_set_lattice1[:,i] = layer_mask_set_lattice1[:,i]^(mask&layer_mask_set_lattice1[:,i])
+            #layer_mask_set_lattice2[:,i] = layer_mask_set_lattice2[:,i]^(mask&layer_mask_set_lattice2[:,i])
+        print('SEARCH DIR,',search_dir)
+        for i in range(np.shape(layer_mask_set_lattice1)[1]):
+            #print(np.shape(layer_mask_set_lattice1))
+            lattice1atoms = atoms[layer_mask_set_lattice1[:,i]]
+            pos_lattice_1_atoms = lattice1atoms.get_positions()
+            #rotate
+            rotated_pos_lattice1 = np.zeros_like(pos_lattice_1_atoms)
+            for j in range(np.shape(pos_lattice_1_atoms)[0]):
+                rotated_pos_lattice1[j,:] = R@pos_lattice_1_atoms[j,:]
 
+            #sort by the y position of the atoms in the rotated frame
+            #that is used to map the pandey reconstruction
+            #print('LAYER MASK LATTICE 1',layer_mask_set_lattice1[:,i])
+            atom_indices = np.where(layer_mask_set_lattice1[:,i] == True)[0]
+            pos_indices = np.argsort(rotated_pos_lattice1[:,1])
+            n_unique = 0
+            #if the search direction is 1 (looking at a lower surface), run this backwards
+            if search_dir == 1:
+                k = len(pos_indices)-1
+            else:
+                k = 0
+            compare_pos = rotated_pos_lattice1[:,1][pos_indices[k]]
+            for j in range(len(pos_indices)):
+                if search_dir == 1:
+                    k = len(pos_indices)-j-1
+                else:
+                    print('here')
+                    k = j
+                #build masks which hold the different types of atoms in
+                #the symmetry broken reconstruction (using the y axis in the rotated
+                # frame to differentiate them from eachother)
+                print(k,search_dir)
+                print(-search_dir*rotated_pos_lattice1[:,1][pos_indices[k]],(-search_dir*(compare_pos-(search_dir*0.01))))
+                if (-search_dir*rotated_pos_lattice1[:,1][pos_indices[k]])>(-search_dir*(compare_pos-(search_dir*0.01))):
+                #if rotated_pos_lattice1[:,1][pos_indices[k]]<(compare_pos-0.01):
+                    n_unique += 1
+                    compare_pos = rotated_pos_lattice1[:,1][pos_indices[k]]
+                if n_unique%3 == 0:
+                    layer_mask_set_lattice1_atom1[:,i][atom_indices[pos_indices[k]]] = True
+                elif n_unique%3 == 1:
+                    layer_mask_set_lattice1_atom2[:,i][atom_indices[pos_indices[k]]] = True
+                elif n_unique%3 == 2:
+                    layer_mask_set_lattice1_atom3[:,i][atom_indices[pos_indices[k]]] = True
+            
+            #repeat for lattice 2
+            lattice2atoms = atoms[layer_mask_set_lattice2[:,i]]
+            pos_lattice_2_atoms = lattice2atoms.get_positions()
+
+            #rotate
+            rotated_pos_lattice2 = np.zeros_like(pos_lattice_2_atoms)
+            for j in range(np.shape(pos_lattice_2_atoms)[0]):
+                rotated_pos_lattice2[j,:] = R@pos_lattice_2_atoms[j,:]
+            #sort by the y position of the atoms in the rotated frame
+            #that is used to map the pandey reconstruction
+            atom_indices = np.where(layer_mask_set_lattice2[:,i] == True)[0]
+            #print(atom_indices)
+            pos_indices = np.argsort(rotated_pos_lattice2[:,1])
+            n_unique = 0
+            if search_dir == 1:
+                k = len(pos_indices)-1
+            else:
+                k = 0
+            compare_pos = rotated_pos_lattice2[:,1][pos_indices[k]]
+            for j in range(len(pos_indices)):
+                if search_dir == 1:
+                    k = len(pos_indices)-j-1
+                else:
+                    k = j
+                #build masks which hold the different types of atoms in
+                #the symmetry broken reconstruction (using the y axis in the rotated
+                # frame to differentiate them from eachother)
+                if (-search_dir*rotated_pos_lattice2[:,1][pos_indices[k]])>(-search_dir*(compare_pos-(search_dir*0.01))):
+                    n_unique += 1
+                    compare_pos = rotated_pos_lattice2[:,1][pos_indices[k]]
+                if n_unique%3 == 0:
+                    layer_mask_set_lattice2_atom1[:,i][atom_indices[pos_indices[k]]] = True
+                elif n_unique%3 == 1:
+                    layer_mask_set_lattice2_atom2[:,i][atom_indices[pos_indices[k]]] = True
+                elif n_unique%3 == 2:
+                    layer_mask_set_lattice2_atom3[:,i][atom_indices[pos_indices[k]]] = True
+
+        
+        #print(layer_mask_set_lattice1_atom1[:,layer])
+        try:
+            atoms.new_array('atomno',np.zeros(len(atoms),dtype=int))
+        except: #TODO add correct exception here (for the case where the array already exists)
+            pass
+        atomno = atoms.arrays['atomno']
+        for layer in range(self.layers):
+            atomno[layer_mask_set_lattice1_atom1[:,layer]] = 1
+            atomno[layer_mask_set_lattice1_atom2[:,layer]] = 2
+            atomno[layer_mask_set_lattice2_atom1[:,layer]] = 3
+            atomno[layer_mask_set_lattice2_atom2[:,layer]] = 4
+        return layer_mask_set_lattice1_atom1, layer_mask_set_lattice1_atom2, layer_mask_set_lattice1_atom3, layer_mask_set_lattice2_atom1, layer_mask_set_lattice2_atom2, layer_mask_set_lattice2_atom3
+
+
+    def apply_si_110(self,
+            atoms,
+            surface_coords,
+            cb=None,
+            xlim=None,
+            ylim=None,
+            zlim=None,
+            search_dir=-1,
+            atoms_for_cb=None,
+            read_from_atoms=False,
+            orientation=0,
+            switch=False):
+        """Function which applies the mapped out Si deformation to the surface layers
+        in a supplied atoms structure.
+        Parameters
+        ----------
+        atoms : ASE atoms object
+            Atoms object which contains surface to relax
+        surface_coords : array_like
+            Coordinates of a point that lies on the free surface (or within one layer distance above it)
+        xlim : array_like
+            [x lower, x upper], x range to apply relaxation
+        ylim : array_like
+            [y lower, y upper], y range to apply relaxation
+        zlim : array_like
+            [z lower, z upper], z range to apply relaxation
+        search_dir : int
+            -1: surface is below point provided, 1: surface is above point provided
+        atoms_for_cb : ASE atoms object
+            atoms object which allows the easy identification of sub-lattices in a multi-lattice
+            crystal using the CubicCauchyBorn set_sublattices function. The difference between this and atoms may
+            be periodic boundary conditions, for example.
+        read_from_atoms : bool
+            Whether or not to read in the cauchy-born sublattices from the atoms structure. This can be done if the
+            atoms object contains arrays called lattice1mask and lattice2mask
+        orientation : 0,1 or 2
+            Which of the three different orientations of the Pandey reconstructions are being applied to the surface
+        """
+        self.eval_cb = cb
+
+        #apply the reconstruction on a diamond structure 110 surface within some set of layers
+        #get masks for atoms
+        layer_mask_set_lattice1_atom1, layer_mask_set_lattice1_atom2, layer_mask_set_lattice1_atom3, layer_mask_set_lattice2_atom1, layer_mask_set_lattice2_atom2,\
+            layer_mask_set_lattice2_atom3 = \
+        self.identify_si_110_layers(atoms,surface_coords,read_from_atoms=read_from_atoms,
+                                    xlim=xlim,ylim=ylim,zlim=zlim,search_dir=search_dir,
+                                    atoms_for_cb=atoms_for_cb,orientation=orientation,switch=switch)
+        ase.io.write(f'{self.directions[self.surf_dir]}1.xyz',atoms)
+
+        for layer in range(self.layers):
+            lattice_1_atom_1shift = self.ral1a1[layer, :]
+            lattice_1_atom_2shift = self.ral1a2[layer, :]
+            lattice_1_atom_3shift = self.ral1a3[layer, :]
+            lattice_2_atom_1shift = self.ral2a1[layer, :]
+            lattice_2_atom_2shift = self.ral2a2[layer, :]
+            lattice_2_atom_3shift = self.ral2a3[layer, :]
+            print(f'atom 1, layer{layer}',lattice_1_atom_1shift)
+            print(f'atom 2, layer{layer}',lattice_1_atom_2shift)
+            print(f'atom 3, layer{layer}',lattice_1_atom_3shift)
+            print(f'atom 4, layer{layer}',lattice_2_atom_1shift)
+            print(f'atom 5, layer{layer}',lattice_2_atom_2shift)
+            print(f'atom 6, layer{layer}',lattice_2_atom_3shift)
+            print(self.surf_dir)
+            #invert component based on surf_dir if necessary
+            
+            if search_dir == 1:
+                #TODO does this rotation matrix make this non-general? what was I thinking? this needs testing
+                R = np.array([[-1,0,0],
+                            [0,-1,0],
+                            [0,0,1]])
+                lattice_1_atom_1shift = R@ \
+                    lattice_1_atom_1shift
+                lattice_1_atom_2shift = R@ \
+                    lattice_1_atom_2shift
+                lattice_1_atom_3shift = R@ \
+                    lattice_1_atom_3shift
+                lattice_2_atom_1shift = R@ \
+                    lattice_2_atom_1shift
+                lattice_2_atom_2shift = R@ \
+                    lattice_2_atom_2shift
+                lattice_2_atom_3shift = R@ \
+                    lattice_2_atom_3shift
+            
+            #apply shifts to different atom types
+            atoms.positions[:,
+                            :][layer_mask_set_lattice1_atom1[:,
+                                                        layer]] += lattice_1_atom_1shift
+            atoms.positions[:,
+                            :][layer_mask_set_lattice1_atom2[:,
+                                                        layer]] += lattice_1_atom_2shift
+            atoms.positions[:,
+                            :][layer_mask_set_lattice1_atom3[:,
+                                                        layer]] += lattice_1_atom_3shift
+            atoms.positions[:,
+                            :][layer_mask_set_lattice2_atom1[:,
+                                                        layer]] += lattice_2_atom_1shift
+            atoms.positions[:,
+                            :][layer_mask_set_lattice2_atom2[:,
+                                                        layer]] += lattice_2_atom_2shift
+            atoms.positions[:,
+                            :][layer_mask_set_lattice2_atom3[:,
+                                                        layer]] += lattice_2_atom_3shift
+
+        atoms.wrap()
diff --git a/scripts/fracture_mechanics/parallel_NCFlex.py b/scripts/fracture_mechanics/parallel_NCFlex.py
index 21ac568d..b3f296fb 100644
--- a/scripts/fracture_mechanics/parallel_NCFlex.py
+++ b/scripts/fracture_mechanics/parallel_NCFlex.py
@@ -168,11 +168,10 @@ def walk(x0,x1,direction,pipe_output,sc_dict):
 
     #time.sleep(10)
     #go idle
-    data_queue.put([os.getpid(),x0,direction])
     status_queue.put([os.getpid(),'idle'],block=False)
     delete_calc(sc)
 
-def get_opt_K_alpha(walk_positions,trail_positions,already_killed_pids):    
+def get_opt_K_alpha(walk_positions,trail_positions,currently_walking_pids):    
     print('GETTING NEW OPT K ALPHA')
     #get a random value of alpha that's not yet been 
     #searched, as well as corresponding estimate for K from
@@ -183,7 +182,7 @@ def get_opt_K_alpha(walk_positions,trail_positions,already_killed_pids):
 
     alpha_covered = 0
     for pid in walk_positions:
-        if pid not in already_killed_pids:
+        if pid in currently_walking_pids:
             [[k_lead,alpha_lead],direction] = walk_positions[pid]
             [k_trail,alpha_trail] = trail_positions[pid]
             alpha_covered += (alpha_lead-alpha_trail)*direction
@@ -335,16 +334,23 @@ def main(K_range,alpha_range):
         sc_dict[pid] = sc
     print(sc_dict)
 
-    #launch half of the workers on a non-blocking initial search
+    #launch workers on a non-blocking initial search
+    #depending on search direction
     #generate the grid of initial points for exploration
+    print('explore direction',explore_direction)
+    if (explore_direction == 1) or (explore_direction ==-1):
+        launch_num = num_processors
+    elif (explore_direction == 0):
+        launch_num = int(num_processors/2)
+
     num_branches = int(np.floor((2*(alpha_range[1]-alpha_range[0]))/alpha_period))
     if num_branches == 0:
-        num_K_vals = int(num_processors/2)
+        num_K_vals = launch_num
     else:
-        num_K_vals = int(np.floor((num_processors/2)/num_branches))
+        num_K_vals = int(np.floor(launch_num/num_branches))
     
     num_K_val_array = [num_K_vals for i in range(num_branches)]
-    sum_diff = int(num_processors/2) - np.sum(num_K_val_array)
+    sum_diff = launch_num - np.sum(num_K_val_array)
     for i in range(sum_diff):
         num_K_val_array[i] += 1
     
@@ -369,12 +375,12 @@ def main(K_range,alpha_range):
 
     it_num = 0
     killed_num = 0
-    already_killed_pids = []
+    currently_walking_pids = []
     curve_explored = False
     percentage_covered = 0
     time.sleep(5)
     while not curve_explored:
-        print(f'already killed pids, {already_killed_pids}')
+        print(f'currently walking pids, {currently_walking_pids}')
         it_num += 1
         #first, check the status queue to update any worker statuses. 
         #get PID list
@@ -410,11 +416,14 @@ def main(K_range,alpha_range):
         
         #print('checking to launch new searches')
         #if there's unnaccounted for idle processes, launch new searches
-        num_new_searches = int(np.floor(idle_num-search_num))
+        if explore_direction == 0:
+            num_new_searches = int(np.floor(idle_num-search_num))
+        elif (explore_direction == 1) or (explore_direction ==-1):
+            num_new_searches = idle_num
         if num_new_searches>0:
             for i in range(num_new_searches):
                 print('LAUNCHING A NEW SEARCH')
-                new_K, new_alpha = get_opt_K_alpha(walk_positions,trail_positions,already_killed_pids)
+                new_K, new_alpha = get_opt_K_alpha(walk_positions,trail_positions,currently_walking_pids)
                 print('INITIAL K, ALPHA OF NEW SEARCH:',new_K,new_alpha)
                 worker_pool.apply_async(search, args=(new_K,new_alpha,sc_dict))
                 time.sleep(1)
@@ -442,7 +451,7 @@ def main(K_range,alpha_range):
                 if len(x0) == 2:
                     #this indicates that the search failed. Immediately restart it
                     print('search failed, restarting')
-                    new_K, new_alpha = get_opt_K_alpha(walk_positions,trail_positions,already_killed_pids)
+                    new_K, new_alpha = get_opt_K_alpha(walk_positions,trail_positions,currently_walking_pids)
                     worker_pool.apply_async(search, args=(new_K,new_alpha,sc_dict))
                     continue
 
@@ -454,7 +463,7 @@ def main(K_range,alpha_range):
 
                 #if found alpha is out of range, start a new search
                 if (found_alpha<alpha_range[0] or found_alpha>alpha_range[1]):
-                    new_K, new_alpha = get_opt_K_alpha(walk_positions,trail_positions,already_killed_pids)
+                    new_K, new_alpha = get_opt_K_alpha(walk_positions,trail_positions,currently_walking_pids)
                     worker_pool.apply_async(search, args=(new_K,new_alpha,sc_dict))
                     continue
 
@@ -470,37 +479,62 @@ def main(K_range,alpha_range):
                         #is at least 0.1 big. 
                         if (alpha_lead>comp_alpha_range[0]) and (alpha_lead<comp_alpha_range[1]):
                             print('alpha detected in already searched range, starting new search')
-                            new_K, new_alpha = get_opt_K_alpha(walk_positions,trail_positions,already_killed_pids)
+                            new_K, new_alpha = get_opt_K_alpha(walk_positions,trail_positions,currently_walking_pids)
                             worker_pool.apply_async(search, args=(new_K,new_alpha,sc_dict))
                             search_restarted=True
                             break
 
                 if search_restarted:
                     continue
+                
+                if explore_direction == 0:
+                    #create two new pipes
+                    pipes = []
+                    pipes.append(multiprocessing.Pipe())
+                    pipes.append(multiprocessing.Pipe())
+
+                    #start walk jobs
+                    worker_pool.apply_async(walk,args=(x0,x1,-1,pipes[0][0],sc_dict))
+                    worker_pool.apply_async(walk,args=(x0,x1,1,pipes[1][0],sc_dict))
+
+                    #send pipeIDs down pipes
+                    pipes[0][1].send(0)
+                    pipes[1][1].send(1)
+
+                    #assign resulting pipe communication to correct pid
+                    queue_empty = False
+                    while not queue_empty:
+                        try:
+                            [pipe_id,pid] = pipe_setup_queue.get(block=True,timeout=3)
+                            pipe_dict[pid] = pipes[pipe_id][1]
+                            #add initial positions as trails to dict
+                            trail_positions[pid] = [x0[-1],x0[-2]] #[K,alpha]
+                            currently_walking_pids.append(pid)
+                        except Empty:
+                            queue_empty = True
+                elif (explore_direction == 1) or (explore_direction ==-1):
+                    #create new pipe
+                    pipes = []
+                    pipes.append(multiprocessing.Pipe())
+
+                    #start walk job
+                    worker_pool.apply_async(walk,args=(x0,x1,explore_direction,pipes[0][0],sc_dict))
+
+                    #send pipeID down pipe
+                    pipes[0][1].send(0)
+
+                    #assign resulting pipe communication to correct pid
+                    queue_empty = False
+                    while not queue_empty:
+                        try:
+                            [pipe_id,pid] = pipe_setup_queue.get(block=True,timeout=3)
+                            pipe_dict[pid] = pipes[pipe_id][1]
+                            #add initial positions as trails to dict
+                            trail_positions[pid] = [x0[-1],x0[-2]] #[K,alpha]
+                            currently_walking_pids.append(pid)
+                        except Empty:
+                            queue_empty = True
 
-                #create two new pipes
-                pipes = []
-                pipes.append(multiprocessing.Pipe())
-                pipes.append(multiprocessing.Pipe())
-
-                #start walk jobs
-                worker_pool.apply_async(walk,args=(x0,x1,-1,pipes[0][0],sc_dict))
-                worker_pool.apply_async(walk,args=(x0,x1,1,pipes[1][0],sc_dict))
-
-                #send pipeIDs down pipes
-                pipes[0][1].send(0)
-                pipes[1][1].send(1)
-
-                #assign resulting pipe communication to correct pid
-                queue_empty = False
-                while not queue_empty:
-                    try:
-                        [pipe_id,pid] = pipe_setup_queue.get(block=True,timeout=3)
-                        pipe_dict[pid] = pipes[pipe_id][1]
-                        #add initial positions as trails to dict
-                        trail_positions[pid] = [x0[-1],x0[-2]] #[K,alpha]
-                    except Empty:
-                        queue_empty = True
 
                 
 
@@ -527,7 +561,7 @@ def main(K_range,alpha_range):
             kill_pids = []
             for item in items:
                 [pid, x, direction] = item
-                if pid in already_killed_pids:
+                if pid not in currently_walking_pids:
                     #if the pid that sent the message
                     #has already had a command sent to kill
                     #but has not yet been killed
@@ -589,7 +623,7 @@ def main(K_range,alpha_range):
             if len(kill_pids)>0:
                 for pid in kill_pids:
                     pipe_dict[pid].send(1) #send kill command
-                    already_killed_pids.append(pid)
+                    currently_walking_pids.remove(pid)
                     killed_num += 1
                     trail_positions[f'killed{killed_num}'] = trail_positions[pid]
                     walk_positions[f'killed{killed_num}'] = walk_positions[pid]
@@ -604,7 +638,6 @@ def main(K_range,alpha_range):
         while not queue_empty:
             try:
                 pid = kill_confirm_queue.get(block=False)
-                already_killed_pids.remove(pid)
                 del trail_positions[pid]
                 del walk_positions[pid]
             except Empty:
@@ -613,7 +646,7 @@ def main(K_range,alpha_range):
         #check if the search is finished and calculate %
         total_alpha_covered = 0
         for pid in walk_positions:
-            if pid in already_killed_pids:
+            if pid not in currently_walking_pids:
                 #need to avoid double counting when a process terminates
                 continue
             [[k_lead,alpha_lead],direction] = walk_positions[pid]
@@ -679,6 +712,7 @@ def main(K_range,alpha_range):
     cutoff = parameter('cutoff')
     dk = parameter('dk', 1e-4)
     dalpha = parameter('dalpha', 1e-1)
+    explore_direction = parameter('explore_direction', 0)
 
     if cb == 'None':
         cb = None
@@ -711,7 +745,7 @@ def main(K_range,alpha_range):
     
     #get k1g
     k1g = crk.k1g(parameter('surface_energy'))
-
+    print('griffthk1,',k1g)
     cluster = params.cluster.copy() 
     if crk.cauchy_born is not None:
         crk.cauchy_born.set_sublattices(cluster,np.transpose(crk.RotationMatrix),read_from_atoms=True)
diff --git a/scripts/fracture_mechanics/parallel_NCFlex_seperate_curves.py b/scripts/fracture_mechanics/parallel_NCFlex_seperate_curves.py
new file mode 100644
index 00000000..03c6b020
--- /dev/null
+++ b/scripts/fracture_mechanics/parallel_NCFlex_seperate_curves.py
@@ -0,0 +1,438 @@
+import matscipy; print(matscipy.__file__)
+
+import os
+import numpy as np
+
+import h5py
+import ase.io
+from ase.units import GPa
+from ase.constraints import FixAtoms
+#from ase.optimize import LBFGSLineSearch
+from ase.optimize.precon import PreconLBFGS
+from ase.optimize import LBFGS
+
+from matscipy import parameter
+from matscipy.elasticity import  fit_elastic_constants
+from matscipy.fracture_mechanics.crack import CubicCrystalCrack, SinclairCrack
+
+from scipy.optimize import fsolve, fminbound, minimize
+from scipy.optimize.nonlin import NoConvergence
+
+import sys
+sys.path.insert(0, '.')
+
+import multiprocessing
+import params
+import time
+
+from queue import Empty
+from copy import deepcopy
+
+
+
+def init_worker():
+    global data_queue
+    print('worker initialised')
+    data_queue.put(os.getpid())
+
+def assign_sc(sc):
+    return [os.getpid(),sc]
+
+def assign_calc(sc):
+    #originally was passing around sc objects, but some
+    #calculators aren't pickleable
+    calc = params.new_calculator_instance()
+    sc.calc = calc
+    sc.atoms.calc = calc
+    sc.cryst.calc = calc
+    if sc.crk.cauchy_born is not None:
+        sc.crk.cauchy_born.calc = calc
+
+def delete_calc(sc):
+    del sc.calc
+    del sc.atoms.calc
+    del sc.cryst.calc
+    if sc.crk.cauchy_born is not None:
+        del sc.crk.cauchy_born
+
+def search(K0,alpha0,sc_dict):
+    global status_queue, search_queue
+    #set status as 'searching'
+
+    status_queue.put([os.getpid(),'searching'],block=False)
+    sc = sc_dict[os.getpid()]
+    assign_calc(sc)
+
+    print(f'pid {os.getpid()} sc is {sc} sc calc is {sc.calc}')
+    sc.variable_k = False
+    k1g = sc.k1g
+    sc.k = K0*k1g
+    sc.alpha = alpha0
+    sc.update_atoms()
+
+    converged = False
+    try:
+        sc.optimize(ftol=0.0005, steps=1000,method='ode12r',precon=True)
+    except:
+        print('did not converge fully in ode12r section, proceeding with krylov')
+
+    try:
+        sc.optimize(ftol=fmax, steps=100,method='krylov')
+        converged = True
+    except:
+        print('No convergence')
+
+    if converged:
+        k_x0 = sc.k/k1g
+        x0 = np.r_[sc.get_dofs(), k_x0 * k1g]
+        alpha_x0 = sc.alpha
+
+        print(f'Rescaling K_I from {sc.k} to {sc.k + dk * k1g}')
+        k_x1 = k_x0 + dk
+        print(k0)
+        sc.k = k_x1*k1g
+        sc.update_atoms()
+        sc.optimize(ftol=0.0001, steps=100,method='krylov')
+        x1 = np.r_[sc.get_dofs(), k_x1 * k1g]
+        alpha_x1 = sc.alpha
+        print(f'k0={k_x0}, k1={k_x1} --> alpha0={alpha_x0}, alpha1={alpha_x1}')
+        
+    else:
+        #return sentinal values
+        #FIXME really this should do something more drastic
+        x0 = [-10000,-10000]
+        x1 = [-10000,-10000]
+
+    #place result in search queue (blocking)
+    print('Braodcasting search result')
+    search_queue.put([os.getpid(),[x0,x1]])
+
+    delete_calc(sc)
+    #at end,
+    #status_queue.put([os.getpid(),'idle'],block=False)
+
+def walk(x0,x1,direction,sc_dict):
+    global data_queue, status_queue
+
+    status_queue.put([os.getpid(),'walking'],block=False)
+    sc = sc_dict[os.getpid()]
+    assign_calc(sc)
+    #set status as 'walking'
+    print('here1')
+    print(f'pid {os.getpid()} sc is {sc} sc calc is {sc.calc}')
+
+    #start NCFlex, passing the data queue and the pipe output
+    # data queue sends back data to main
+    # pipe output allows main to kill the process
+    print('starting ncflex...')
+    sc.variable_k = True
+    precon = False
+    print('lammps lib process',os.getpid(),sc.calc)
+    traj_file_name = f'curve_alph_{np.round(x0[-2],decimals=3)}_K_{np.round(x0[-1]/sc.k1g,decimals=3)}_dir_{direction}.h5'
+    sc.arc_length_continuation(x0, x1, N=nsteps,
+                            ds=ds, ds_max=0.05, ftol=fmax, max_steps=10,
+                            direction=direction,
+                            continuation=continuation,
+                            traj_file=traj_file_name,
+                            traj_interval=traj_interval,
+                            precon=precon,opt_method='krylov',parallel=False,
+                            allow_alpha_backtracking=allow_alpha_backtracking)
+
+    #time.sleep(10)
+    #go idle
+    status_queue.put([os.getpid(),'idle'],block=False)
+    delete_calc(sc)
+
+
+def main(K_range,alpha_range):
+    global data_queue, sc, status_queue,search_queue
+
+    #kill old sc calculator
+    del sc.calc
+    del sc.atoms.calc
+    del sc.cryst.calc
+    if sc.crk.cauchy_born is not None:
+        del sc.crk.cauchy_born.calc
+
+    #create copies of sc for each process
+    sc_array = [deepcopy(sc) for i in range(num_processors)]
+    print(sc_array)
+
+    #now, start worker pool
+    worker_pool = multiprocessing.Pool(num_processors,initializer=init_worker)
+
+    #set up dictionary to hold each of the sc objects that workers need
+    sc_dict = {}
+
+    #set up dictionary to hold worker status
+    worker_status = {}
+
+    #get PID list
+    queue_empty = False
+    while not queue_empty:
+        try:
+            pid = data_queue.get(block=True,timeout=5)
+            worker_status[pid] = 'idle'
+        except Empty:
+            queue_empty = True
+
+    #assign scs
+    print('assigning SCs')
+    for i in range(num_processors):
+        #make a new calculator instance
+        calc = params.new_calculator_instance()
+        [pid,sc] = worker_pool.apply(assign_sc, args=(tuple([sc_array[i]])))
+        sc_dict[pid] = sc
+    print(sc_dict)
+
+    #launch workers on a non-blocking initial search
+    #depending on search direction
+    #generate the grid of initial points for exploration
+    print('explore direction',explore_direction)
+    if (explore_direction == 1) or (explore_direction ==-1):
+        launch_num = num_processors
+    elif (explore_direction == 0):
+        launch_num = int(num_processors/2)
+
+    num_branches = int(np.floor((2*(alpha_range[1]-alpha_range[0]))/alpha_period))
+    if num_branches == 0:
+        num_K_vals = launch_num
+    else:
+        num_K_vals = int(np.floor(launch_num/num_branches))
+    
+    num_K_val_array = [num_K_vals for i in range(num_branches)]
+    sum_diff = launch_num - np.sum(num_K_val_array)
+    for i in range(sum_diff):
+        num_K_val_array[i] += 1
+    
+    #FIXME this assumes that alpha_range > alpha_period/4 (should be usually)
+
+    init_alpha_vals = np.array([alpha_range[0] + (alpha_period/4) + ((alpha_period/2)*i) for i in range(num_branches)])
+
+    #alpha_range_buffer = (alpha_range[1]-alpha_range[0])*(1/(num_processors/2))
+    #alphas = np.linspace(alpha_range[0]+alpha_range_buffer,alpha_range[1]-alpha_range_buffer,int(num_processors/2))
+    #Ks = np.random.uniform(K_range[0],K_range[1],int(num_processors/2))
+
+    #launch searches
+    for i, alpha in enumerate(init_alpha_vals):
+        if num_K_val_array[i] == 0:
+            continue
+        else:
+            K_step = (K_range[1]-K_range[0])/num_K_val_array[i]
+            initial_K_step = np.random.uniform(0,K_step)
+            init_K_vals = np.array([K_range[0] + initial_K_step + K_step*p for p in range(num_K_val_array[i])])
+            for K in init_K_vals:
+                worker_pool.apply_async(search, args=(K,alpha,sc_dict))
+
+    it_num = 0
+    killed_num = 0
+    curve_explored = False
+    percentage_covered = 0
+    time.sleep(5)
+    while not curve_explored:
+        it_num += 1
+        #first, check the status queue to update any worker statuses. 
+        #get PID list
+        queue_empty = False
+        #print('checking status queue')
+        while not queue_empty:
+            try:
+                [pid,status] = status_queue.get(block=False)
+                worker_status[pid] = status
+            except Empty:
+                queue_empty = True
+        
+        #Write worker status to file and then launch any new workers if necessary
+        #print('writing to file')
+        if it_num%1 == 0:
+            with open('worker_status.txt', 'w') as f:
+                for pid in worker_status:
+                    f.write(f'{pid} : {worker_status[pid]}\n')
+        #print('checking worker status')
+        #check what workers are doing
+        search_num = 0
+        walk_num = 0
+        idle_num = 0
+        for pid in worker_status:
+            if worker_status[pid] == 'idle':
+                idle_num += 1
+            elif worker_status[pid] == 'searching':
+                search_num += 1
+            elif worker_status[pid] == 'walking':
+                walk_num += 1
+        if it_num % 100 == 0:
+            print(f'idle: {idle_num}, search, {search_num}, walk {walk_num}')
+        
+        #print('checking to launch new searches')
+        #if there's unnaccounted for idle processes, launch new searches
+        if explore_direction == 0:
+            num_new_searches = int(np.floor(idle_num-search_num))
+        elif (explore_direction == 1) or (explore_direction ==-1):
+            num_new_searches = idle_num
+        if num_new_searches>0:
+            for i in range(num_new_searches):
+                print('LAUNCHING A NEW SEARCH')
+                new_K = np.random.uniform(K_range[0],K_range[1])
+                new_alpha = np.random.uniform(alpha_range[0],alpha_range[1])
+                print('INITIAL K, ALPHA OF NEW SEARCH:',new_K,new_alpha)
+                worker_pool.apply_async(search, args=(new_K,new_alpha,sc_dict))
+                time.sleep(1)
+            
+
+        #print('checking for finished searches')
+        #check for any finished searches in the search_queue
+        #for any finished searches, start two walk jobs from the found position in each
+        #direction (if the search position lies in the alpha range). If it doesn't, start a new search.
+        queue_empty = False
+        search_results = []
+        while not queue_empty:
+            try:
+                item = search_queue.get(block=False)
+                search_results.append(item)
+            except Empty:
+                queue_empty = True
+
+        if len(search_results)>0:
+            print('found finished search!')
+            print(search_results)
+            for res in search_results:
+                [pid,[x0,x1]] = res
+
+                if len(x0) == 2:
+                    #this indicates that the search failed. Immediately restart it
+                    print('search failed, restarting')
+                    new_K = np.random.uniform(K_range[0],K_range[1])
+                    new_alpha = np.random.uniform(alpha_range[0],alpha_range[1])
+                    worker_pool.apply_async(search, args=(new_K,new_alpha,sc_dict))
+                    continue
+
+                found_alpha = x0[-2]
+
+                #set worker status to idle here rather than in function to prevent
+                #new searches being started erroneously
+                worker_status[pid] = 'idle'
+
+                #if found alpha is out of range, start a new search
+                if (found_alpha<alpha_range[0] or found_alpha>alpha_range[1]):
+                    new_K = np.random.uniform(K_range[0],K_range[1])
+                    new_alpha = np.random.uniform(alpha_range[0],alpha_range[1])
+                    worker_pool.apply_async(search, args=(new_K,new_alpha,sc_dict))
+                    continue
+                
+                if explore_direction == 0:
+                    #start walk jobs
+                    worker_pool.apply_async(walk,args=(x0,x1,-1,sc_dict))
+                    worker_pool.apply_async(walk,args=(x0,x1,1,sc_dict))
+
+
+                elif (explore_direction == 1) or (explore_direction ==-1):
+                    #start walk job
+                    worker_pool.apply_async(walk,args=(x0,x1,explore_direction,sc_dict))
+
+    
+    #kill all workers outside of while loop
+    worker_pool.close()
+
+
+
+if __name__ == '__main__':
+    
+    #multiprocessing configuration
+    num_processors = parameter('num_processors')
+    K_range = parameter('K_range')
+    alpha_range = parameter('alpha_range')
+    
+    # Set start method to fork so 
+    # each process inherits information from the parent process
+    # this includes the queue through which information is passed
+    multiprocessing.set_start_method('fork')
+    global data_queue, pipe_setup_queue, status_queue, search_queue, kill_confirm_queue
+
+    data_queue = multiprocessing.Queue()
+    status_queue = multiprocessing.Queue()
+    search_queue = multiprocessing.Queue()
+
+    #read parameters in from params file
+    calc = parameter('calc')
+    fmax = parameter('fmax', 1e-3)
+    max_opt_steps = parameter('max_opt_steps', 100)
+    max_steps = parameter('max_arc_steps', 10)
+    vacuum = parameter('vacuum', 10.0)
+    flexible = parameter('flexible', True)
+    continuation = parameter('continuation', False)
+    ds = parameter('ds', 1e-2)
+    nsteps = parameter('nsteps', 10)
+    a0 = parameter('a0') # lattice constant
+    k0 = parameter('k0', 1.0)
+    extended_far_field = parameter('extended_far_field', False)
+    alpha0 = parameter('alpha0', 0.0) # initial guess for crack position
+    dump = parameter('dump', False)
+    precon = parameter('precon', False)
+    traj_file = parameter('traj_file', 'x_traj.h5')
+    restart_file = parameter('restart_file', traj_file)
+    traj_interval = parameter('traj_interval', 1)
+    direction = parameter('direction', +1)
+    ds_max = parameter('ds_max', 0.1)
+    ds_min = parameter('ds_min', 1e-6)
+    ds_aggressiveness=parameter('ds_aggressiveness', 1.25)
+    opt_method=parameter('opt_method', 'krylov')
+    cb = parameter('cb', 'None')
+    rI = parameter('r_I')
+    rIII = parameter('r_III')
+    cutoff = parameter('cutoff')
+    dk = parameter('dk', 1e-4)
+    dalpha = parameter('dalpha', 1e-1)
+    explore_direction = parameter('explore_direction', 0)
+    allow_alpha_backtracking = parameter('allow_alpha_backtracking',False)
+
+    if cb == 'None':
+        cb = None
+    
+    cryst = params.cryst.copy()
+
+    pos = cryst.get_positions()
+    sx, sy, sz = cryst.cell.diagonal()
+    xpos = pos[:,0] - sx/2
+    ypos = pos[:,1] - sy/2
+
+    #find the closest y atoms by finding which atoms lie within 1e-2 of the min
+    #filter out all atoms with x less than 0
+    xmask = xpos>0
+    closest_y_mask = np.abs(ypos)<(np.min(np.abs(ypos[xmask]))+(1e-2))
+
+
+    #find the x positions of these atoms
+    closest_x = xpos[closest_y_mask&xmask]
+
+    #sort these atoms and find the largest x gap
+    sorted_x = np.sort(closest_x)
+    diffs = np.diff(sorted_x)
+    alpha_period = np.sum(np.unique(np.round(np.diff(sorted_x),decimals=4)))
+    print('alpha_period',alpha_period)
+    # setup the crack
+    crk = CubicCrystalCrack(parameter('crack_surface'),
+                    parameter('crack_front'),
+                    C=parameter('C')/GPa,cauchy_born=cb) #TODO change back to cb
+    
+    #get k1g
+    k1g = crk.k1g(parameter('surface_energy'))
+    print('griffthk1,',k1g)
+    cluster = params.cluster.copy() 
+    if crk.cauchy_born is not None:
+        crk.cauchy_born.set_sublattices(cluster,np.transpose(crk.RotationMatrix),read_from_atoms=True)
+
+    global sc
+
+    sc = SinclairCrack(crk, cluster, calc, k0 * k1g,
+                    alpha=alpha0,
+                    vacuum=vacuum,
+                    variable_alpha=flexible,
+                    extended_far_field=extended_far_field,rI=rI,
+                    rIII=rIII,cutoff=cutoff,incl_rI_f_alpha=True)
+    sc.k1g = k1g
+
+    sc.variable_alpha = True #True
+    sc.variable_k = True
+
+    main(K_range,alpha_range)
+
diff --git a/scripts/fracture_mechanics/parallel_mode_II_NCFlex.py b/scripts/fracture_mechanics/parallel_mode_II_NCFlex.py
new file mode 100644
index 00000000..703586cc
--- /dev/null
+++ b/scripts/fracture_mechanics/parallel_mode_II_NCFlex.py
@@ -0,0 +1,195 @@
+import matscipy; print(matscipy.__file__)
+from matscipy import parameter
+
+from mpi4py import MPI
+from matscipy.fracture_mechanics.crack import CubicCrystalCrack, SinclairCrack
+
+from ase.units import GPa
+import numpy as np
+
+from copy import deepcopy
+
+import h5py
+
+comm = MPI.COMM_WORLD
+rank = comm.Get_rank()
+import params
+import os
+
+import time
+
+num_processors = parameter('num_processors')
+
+#read parameters in from params file
+calc = parameter('calc')
+fmax = parameter('fmax', 1e-3)
+max_opt_steps = parameter('max_opt_steps', 100)
+max_steps = parameter('max_arc_steps', 10)
+vacuum = parameter('vacuum', 10.0)
+flexible = parameter('flexible', True)
+continuation = parameter('continuation', False)
+ds = parameter('ds', 1e-2)
+nsteps = parameter('nsteps', 10)
+a0 = parameter('a0') # lattice constant
+k0 = parameter('k0', 1.0)
+extended_far_field = parameter('extended_far_field', False)
+alpha0 = parameter('alpha0', 0.0) # initial guess for crack position
+dump = parameter('dump', False)
+precon = parameter('precon', False)
+traj_file = parameter('traj_file', 'x_traj.h5')
+restart_file = parameter('restart_file', traj_file)
+traj_interval = parameter('traj_interval', 1)
+direction = parameter('direction', +1)
+ds_max = parameter('ds_max', 0.1)
+ds_min = parameter('ds_min', 1e-6)
+ds_aggressiveness=parameter('ds_aggressiveness', 1.25)
+opt_method=parameter('opt_method', 'krylov')
+cb = parameter('cb', 'None')
+rI = parameter('r_I')
+rIII = parameter('r_III')
+cutoff = parameter('cutoff')
+dk = parameter('dk', 1e-4)
+dalpha = parameter('dalpha', 1e-1)
+explore_direction = parameter('explore_direction', 0)
+allow_alpha_backtracking = parameter('allow_alpha_backtracking',False)
+k1_curve_file = parameter('k1_curve_file')
+num_data_points = parameter('num_data_points')
+follow_G_contour = parameter('follow_G_contour',False)
+folder_name = parameter('folder_name','data')
+if cb == 'None':
+    cb = None
+    
+if cb is not None:
+    if rank == 0:
+        #cb.initial_regression_fit()
+        #cb.save_regression_model()
+        cb.load_regression_model()
+        cb_mod = cb.get_model()
+        #communicate cb model to all other processors using mpi broadcast
+    else:
+        cb_mod = 0
+    cb_mod = comm.bcast(cb_mod, root=0)
+    cb.set_model(cb_mod)
+
+cryst = params.cryst.copy()
+
+pos = cryst.get_positions()
+sx, sy, sz = cryst.cell.diagonal()
+xpos = pos[:,0] - sx/2
+ypos = pos[:,1] - sy/2
+
+crk = CubicCrystalCrack(parameter('crack_surface'),
+                parameter('crack_front'),
+                C=parameter('C')/GPa,cauchy_born=cb)
+
+
+k1g = crk.k1g(parameter('surface_energy'))
+print('griffthk1,',k1g)
+cluster = params.cluster.copy() 
+if crk.cauchy_born is not None:
+    crk.cauchy_born.set_sublattices(cluster,np.transpose(crk.RotationMatrix),read_from_atoms=True)
+
+
+sc = SinclairCrack(crk, cluster, calc, k0 * k1g,
+                alpha=alpha0,
+                vacuum=vacuum,
+                variable_alpha=flexible,
+                extended_far_field=extended_far_field,rI=rI,
+                rIII=rIII,cutoff=cutoff,incl_rI_f_alpha=True)
+sc.k1g = k1g
+sc.variable_alpha = True #True
+sc.variable_k = True
+sc.cont_k = 'k2'
+
+
+if rank == 0:
+    #now read in the initial data file and start to send over data
+    hf = h5py.File(k1_curve_file, 'r')
+    x_traj = hf['x']
+    x = x_traj[:,:]
+    hf.close()
+
+    total_data_points = np.shape(x)[0]
+    sample_frac = int(total_data_points/num_data_points)
+    x = x[::sample_frac,:]
+
+    data_points_per_proc = int(np.floor(np.shape(x)[0]/num_processors))
+    #k2 distance to walk size
+
+    num_data_point_array = [data_points_per_proc for i in range(num_processors)]
+    sum_diff = np.shape(x)[0] - np.sum(num_data_point_array)
+    for i in range(sum_diff):
+        num_data_point_array[i] += 1
+    print(num_data_point_array)
+    time.sleep(3)
+    #now communicate over all data
+    for proc_num in range(1,num_processors):
+        print(proc_num)
+        comm.send(num_data_point_array[proc_num], dest=proc_num, tag=12)
+        comm.send(x[sum(num_data_point_array[:proc_num]):sum(num_data_point_array[:proc_num+1]),:], dest=proc_num, tag=10)
+        comm.send(os.getpid(),dest=proc_num,tag=13)
+    #the first set of data (proc_num 0) is assigned to rank 0 (this processor)
+    num_data_points = num_data_point_array[0]
+    x = x[0:num_data_points,:]
+    pid = os.getpid()
+    
+else:
+    num_data_points = comm.recv(source=0, tag=12)
+    x = comm.recv(source=0, tag=10)
+
+    #recieve pid
+    pid = comm.recv(source=0, tag=13)
+    
+for i in range(num_data_points):
+    assert sc.cont_k == 'k2'
+    #pull kII0 from x
+    sc.variable_k = True
+    x0 = x[i,:]
+    sc.set_dofs(x0)
+    sc.variable_k = False
+    converged=False
+    dkcurr = dk
+    for attempt in range(10):
+        kII0 = x[i,-1]
+        k_x0 = kII0/sc.k1g
+        alpha_x0 = x[i,-3]
+        print(f'Rescaling K_II from {kII0} to {kII0 + dkcurr}')
+        k_x1 = k_x0 + dkcurr
+        sc.kII = k_x1*sc.k1g
+        if follow_G_contour:
+            print('before',sc.kI)
+            sc.kI = np.sqrt(((k_x0*sc.k1g)**2 + sc.kI**2) - (k_x1*sc.k1g)**2)
+            print('after',sc.kI)
+            time.sleep(5)
+        
+        sc.update_atoms()
+        print('starting search...')
+        try:
+            sc.optimize(ftol=0.0001, steps=100,method='krylov')
+            converged=True
+            break
+        except RuntimeError:
+            print('No convergence, retrying with smaller k step')
+            dkcurr = dkcurr/2
+    
+    if converged == False:
+        print('No convergence, skipping this point')
+        continue
+    x1 = np.r_[sc.get_dofs(), k_x1 * sc.k1g]
+    alpha_x1 = sc.alpha
+    print(f'k0={k_x0}, k1={k_x1} --> alpha0={alpha_x0}, alpha1={alpha_x1}')
+
+    #start that continuation
+    sc.variable_k = True
+    traj_file_name = f'{folder_name}/main_process_{pid}_sub_process_{rank}_point_{i}.h5'
+    print('starting ncflex...')
+    sc.arc_length_continuation(x0, x1, N=nsteps,
+                    ds=ds, ds_max=0.05, ftol=fmax, max_steps=10,
+                    direction=explore_direction,
+                    continuation=continuation,
+                    traj_file=traj_file_name,
+                    traj_interval=traj_interval,
+                    precon=precon,opt_method='krylov',parallel=False,
+                    allow_alpha_backtracking=allow_alpha_backtracking,
+                    follow_G_contour=follow_G_contour)
+


Test patch:
diff --git a/scripts/fracture_mechanics/parallel_rI_convergence_test.py b/scripts/fracture_mechanics/parallel_rI_convergence_test.py
new file mode 100644
index 00000000..4d00ab00
--- /dev/null
+++ b/scripts/fracture_mechanics/parallel_rI_convergence_test.py
@@ -0,0 +1,190 @@
+import matscipy; print(matscipy.__file__)
+from matscipy import parameter
+
+from mpi4py import MPI
+from matscipy.fracture_mechanics.crack import CubicCrystalCrack, SinclairCrack
+
+from ase.units import GPa
+import numpy as np
+
+import h5py
+
+comm = MPI.COMM_WORLD
+rank = comm.Get_rank()
+import params
+import os
+
+import time
+
+num_processors = parameter('num_processors')
+#read parameters in from params file
+calc = parameter('calc')
+fmax = parameter('fmax', 1e-3)
+cb = parameter('cb', 'None')
+r_I_vals = parameter('r_I_vals')
+r_III_vals = parameter('r_III_vals')
+cutoff = parameter('cutoff')
+k1_curve_file = parameter('k1_curve_file')
+num_data_points = parameter('num_data_points')
+follow_G_contour = parameter('follow_G_contour',False)
+folder_name = parameter('folder_name','data')
+a0 = parameter('a0') # lattice constant
+k0 = parameter('k0', 1.0)
+alpha0 = parameter('alpha0', 0.0) # initial guess for crack position
+traj_file = parameter('traj_file', 'x_traj.h5')
+extended_far_field = parameter('extended_far_field', False)
+vacuum = parameter('vacuum', 10.0)
+flexible = parameter('flexible', True)
+cont_k = parameter('cont_k','k1')
+clusters = parameter('clusters')#get full set of clusters being used
+
+if cb == 'None':
+    cb = None
+    
+if cb is not None:
+    if rank == 0:
+        #cb.initial_regression_fit()
+        #cb.save_regression_model()
+        cb.load_regression_model()
+        cb_mod = cb.get_model()
+        #communicate cb model to all other processors using mpi broadcast
+    else:
+        cb_mod = 0
+    cb_mod = comm.bcast(cb_mod, root=0)
+    cb.set_model(cb_mod)
+
+crk = CubicCrystalCrack(parameter('crack_surface'),
+                parameter('crack_front'),
+                C=parameter('C')/GPa,cauchy_born=cb)
+
+k1g = crk.k1g(parameter('surface_energy'))
+print('griffthk1,',k1g)
+
+if rank == 0:
+    #now read in the initial data file and start to send over data
+    hf = h5py.File(k1_curve_file, 'r')
+    x_traj = hf['x']
+    x = x_traj[:,:]
+    hf.close()
+
+    total_data_points = np.shape(x)[0]
+    sample_frac = int(total_data_points/num_data_points)
+    x = x[::sample_frac,:]
+
+    data_points_per_proc = int(np.floor(np.shape(x)[0]/num_processors))
+    #k2 distance to walk size
+
+    num_data_point_array = [data_points_per_proc for i in range(num_processors)]
+    sum_diff = np.shape(x)[0] - np.sum(num_data_point_array)
+    for i in range(sum_diff):
+        num_data_point_array[i] += 1
+    print(num_data_point_array)
+    time.sleep(3)
+    #now communicate over all data
+    for proc_num in range(1,num_processors):
+        print(proc_num)
+        comm.send(num_data_point_array[proc_num], dest=proc_num, tag=12)
+        comm.send(x[sum(num_data_point_array[:proc_num]):sum(num_data_point_array[:proc_num+1]),:], dest=proc_num, tag=10)
+        comm.send(os.getpid(),dest=proc_num,tag=13)
+    #the first set of data (proc_num 0) is assigned to rank 0 (this processor)
+    num_data_points = num_data_point_array[0]
+    x = x[0:num_data_points,:]
+    pid = os.getpid()
+else:
+    num_data_points = comm.recv(source=0, tag=12)
+    x = comm.recv(source=0, tag=10)
+
+    #recieve pid
+    pid = comm.recv(source=0, tag=13)
+
+for dp in range(num_data_points):
+    #each process creates a file titled f'{pid}_{x[dp][-1]}'
+    #this file contains the data for the current data point
+    #the file is created in the folder f'{folder_name}'
+
+    #create folder if it doesn't exist
+    if not os.path.exists(folder_name):
+        os.mkdir(folder_name)
+    #create file, with names rouneded to 3dp
+
+    file_name = f'{folder_name}/{pid}_rank_{rank}_point_{dp}.h5'
+    #create file
+    hf = h5py.File(file_name,'w')
+    for i, curr_cluster in enumerate(clusters):
+        print(r_I_vals[i],r_III_vals[i])
+        cluster = curr_cluster.copy()
+        if crk.cauchy_born is not None:
+            crk.cauchy_born.set_sublattices(cluster,np.transpose(crk.RotationMatrix),read_from_atoms=True)
+
+        sc = SinclairCrack(crk, cluster, calc, k0 * k1g,
+                    alpha=alpha0,
+                    vacuum=vacuum,
+                    variable_alpha=flexible,
+                    extended_far_field=extended_far_field,rI=r_I_vals[i],
+                    rIII=r_III_vals[i],cutoff=cutoff,incl_rI_f_alpha=True)
+        
+        sc.write_atoms_to_file()
+        sc.k1g = k1g
+        sc.variable_alpha = True #True
+        sc.variable_k=True
+        sc.cont_k = cont_k
+
+        if i == 0:
+            #on the first iteration, just set prev_relaxed to the initial config
+            #and pass
+            x_initial = x[dp]
+            sc.set_dofs(x_initial[:])
+            #sc.write_atoms_to_file()
+            #check if both kI and kII are in x[dp], if they are then
+            #remove whichever is involved in continuation
+            if len(x_initial) == len(sc)+1:
+                if sc.cont_k == 'k1':
+                    idx = -1
+                    sc.kII = x_initial[idx]
+                    kII0 = sc.kII
+                elif sc.cont_k == 'k2':
+                    idx = -2
+                    sc.kI = x_initial[idx]
+                    kI0 = sc.kI
+                x_initial = np.delete(x_initial,idx)
+            
+            prev_u, prev_alph, prev_k = sc.unpack(x_initial[:],reshape=True)
+            continue
+
+        #now take sc.atoms, mask out the inner r_I_vals[i-1] atoms, and set to
+        sx, sy, sz = cluster.cell.diagonal()
+        pos_rI = cluster.get_positions()[sc.regionI]
+        xcoord, ycoord = pos_rI[:, 0], pos_rI[:, 1]
+        cx, cy = sx/2, sy/2
+        r = np.sqrt((xcoord - cx)**2 + (ycoord - cy)**2)
+        mask = r<r_I_vals[i-1]
+        sc.alpha = prev_alph
+        if sc.cont_k == 'k1':
+            sc.kI = prev_k
+            sc.kII = kII0
+        elif sc.cont_k == 'k2':
+            sc.kII = prev_k
+            sc.kI = kI0
+        
+        sc.u[mask] = prev_u
+        sc.update_atoms()
+        
+        #sc.write_atoms_to_file()
+        sc.variable_k=False
+        #now relax
+        try:
+            sc.optimize(ftol=0.0001, steps=100,method='krylov')
+        except RuntimeError:
+            #proceed to next data point
+            break
+        sc.write_atoms_to_file()
+        #get out data
+        sc.variable_k = True
+
+        x_new = sc.get_dofs()
+        prev_u, prev_alph, prev_k = sc.unpack(x_new,reshape=True)
+
+        #write x_new to a new file dataset called f'r_I_vals[i]'
+        with h5py.File(file_name,'a') as hf:
+            hf.create_dataset(f'r_I_{r_I_vals[i]}',data=x_new, compression='gzip')
+        
\ No newline at end of file
diff --git a/scripts/fracture_mechanics/parallel_rI_convergence_test_post_process.py b/scripts/fracture_mechanics/parallel_rI_convergence_test_post_process.py
new file mode 100644
index 00000000..fdcd44a5
--- /dev/null
+++ b/scripts/fracture_mechanics/parallel_rI_convergence_test_post_process.py
@@ -0,0 +1,59 @@
+import h5py
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+import params
+from matscipy import parameter
+
+folder_name = parameter('folder_name','data')
+pid = parameter('data_load_pid')
+traj_files = [file for file in os.listdir(folder_name) if pid in file]
+
+data_sets = {}
+#now read each file and append the data to the data_sets for each key
+
+for file in traj_files:
+    hf = h5py.File(os.path.join(folder_name,file),'r')
+    keys = list(hf.keys())
+    for key in keys:
+        #if key is not in data sets, add it
+        if key not in data_sets:
+            data_sets[key] = []
+        data_sets[key].append(hf[key][:])
+    hf.close()
+print(data_sets)
+
+#now save each data set to new h5 files in the folder
+#set keys to data_sets.keys()
+keys = data_sets.keys()
+for key in keys:
+    hf = h5py.File(os.path.join(folder_name,f'{pid}_fin_{key}.h5'),'w')
+    hf.create_dataset('x',(0, len(data_sets[key][0])),maxshape=(None,len(data_sets[key][0])),compression='gzip')
+    #write each entry of data_sets[key] to the file
+    for data in data_sets[key]:
+        hf['x'].resize((hf['x'].shape[0]+1),axis=0)
+        hf['x'][-1,:] = data
+    hf.close()
+
+
+##open all traj files in folder that have just been made and plot the k-alpha curves
+#traj_files = [file for file in os.listdir(folder_name) if f'{pid}_fin' in file]
+#for i,traj_file in enumerate(traj_files):
+#    hf = h5py.File(os.path.join(folder_name,traj_file),'r')
+#    x_traj = hf['x'][:]
+#    #sort x_traj by second to last column
+#    x_traj = x_traj[x_traj[:,-2].argsort()]
+#    Ks = x_traj[:,-1]
+#    print(len(Ks))
+#    alphas = x_traj[:,-2]
+#    hf.close()
+#    print(Ks,alphas)
+#    plt.plot(Ks,alphas,label=traj_file)###
+
+#plt.legend()
+#plt.xlabel('K')
+#plt.ylabel('alpha')
+#plt.title('NCFlex diamond (1-10) solution curve')
+#plt.tight_layout()
+#plt.grid()
+#plt.savefig('rI_conv_curve.png')
\ No newline at end of file
diff --git a/tests/test_build_3D_crack.py b/tests/test_build_3D_crack.py
new file mode 100644
index 00000000..b134dc83
--- /dev/null
+++ b/tests/test_build_3D_crack.py
@@ -0,0 +1,81 @@
+import unittest
+import numpy as np
+from ase.build import bulk
+from ase.lattice.cubic import Diamond
+from ase.optimize.precon import PreconLBFGS
+from ase.constraints import ExpCellFilter
+from matscipy.cauchy_born import CubicCauchyBorn
+from matscipy.calculators.manybody.explicit_forms.stillinger_weber import StillingerWeber, Holland_Marder_PRL_80_746_Si
+from matscipy.calculators.manybody import Manybody
+from matscipy.fracture_mechanics.clusters import generate_3D_cubic_111,diamond,set_regions,generate_3D_structure
+import matscipytest
+import ase.io
+class TestBuild3DPeriodicFractureCell(matscipytest.MatSciPyTestCase):
+    """
+    Test to build 3D kink-periodic fracture cell. 
+    """
+
+    def setUp(self):
+        el = 'Si'  # chemical element
+        si = Diamond(el, latticeconstant=5.43, size=[1, 1, 1],
+                     directions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]])
+        calc = Manybody(**StillingerWeber(Holland_Marder_PRL_80_746_Si))
+        si.calc = calc
+        self.calc = calc
+        ecf = ExpCellFilter(si)
+        opt = PreconLBFGS(ecf)
+        opt.run(fmax=1e-6, smax=1e-6)
+        # print(si.get_positions())
+        a0 = si.cell[0, 0]
+        # print(a0)
+        self.cb = CubicCauchyBorn(el, a0, calc, lattice=Diamond)
+        self.unit_cell = si
+        self.a0 = a0
+        self.el = el
+    
+    def build_3D_cluster(self,build_func,nzlayer,crack_surface,crack_front):
+
+        A = np.zeros([3, 3])
+        directions = [np.cross(crack_surface,crack_front),crack_surface,crack_front]
+        for i, direction in enumerate(directions):
+            direction = np.array(direction)
+            A[:, i] = direction / np.linalg.norm(direction)
+    
+        r_III = 40
+        cutoff = 6
+        padding = 3
+        r_I = 20
+        [ax,ay,az] = diamond(self.el, self.a0, [1,1,1], crack_surface, crack_front).cell.lengths()
+        n = [2 * (int(np.ceil(((r_III + 2*cutoff)/ ax)))+padding)+1, 2 * (int(np.ceil((r_III + 2*cutoff)/ ay))+padding+1), 1]
+
+        cryst = diamond(self.el, self.a0, n, crack_surface, crack_front, cb=self.cb, shift=np.array([0.0,0.25,0.0]),switch_sublattices=True)
+        self.cb.set_sublattices(cryst,A,read_from_atoms=True)
+
+        cryst, theta = build_func(cryst,nzlayer,self.el,self.a0,diamond,crack_surface,crack_front,cb=self.cb,shift=np.array([0.0,0.25,0.0]),switch_sublattices=True)
+
+        #create cluster from this surface relaxed version
+        cluster = set_regions(cryst, r_I, cutoff, r_III)  # carve circular cluster
+
+        #assert that forces between 5 and 10 angstroms are all small
+        sx, sy, sz = cluster.cell.diagonal()
+        x, y = cluster.positions[:, 0], cluster.positions[:, 1]
+        cx, cy = sx/2, sy/2
+        r = np.sqrt((x - cx)**2 + (y - cy)**2)
+        
+        mask = (r>5) & (r<10)
+        cluster.calc = self.calc
+        forces = cluster.get_forces()
+        #ase.io.write('cluster.xyz',cluster)
+        assert np.max(np.linalg.norm(forces[mask],axis=1))<0.5
+
+    def test_build_cubic_111_periodic_cell(self):
+        nzlayer = 5
+        crack_surface = np.array([1, 1, 1])
+        crack_front = np.array([1, -1, 0])
+        self.build_3D_cluster(generate_3D_cubic_111,nzlayer,crack_surface,crack_front)
+
+    def test_build_3D_kink_periodic_cell(self):
+        nzlayer = 5
+        crack_surface = np.array([1, 1, 1])
+        crack_front = np.array([1, -1, 0])
+        self.build_3D_cluster(generate_3D_structure,nzlayer,crack_surface,crack_front)
diff --git a/tests/test_cauchy_born_corrector.py b/tests/test_cauchy_born_corrector.py
index de014f1b..976113f6 100644
--- a/tests/test_cauchy_born_corrector.py
+++ b/tests/test_cauchy_born_corrector.py
@@ -301,7 +301,7 @@ def test_regression_model_gradient_F(self):
         coordinates = 'cart3D'
         #print('MAKING RAW PREDICION WITH NO ADDED EPS')
         atoms, shifts, shift_err_before, A = self.model_prediction(
-            dirs, eps, method='regression', F_func=func, coordinates=coordinates, atol=1e-4, returnvals=True)
+            dirs, eps, method='regression', F_func=func, coordinates=coordinates, atol=2e-4, returnvals=True)
         
         #find the gradient using the function finite differences
         de = 1e-5
diff --git a/tests/test_cluster_stable_sort.py b/tests/test_cluster_stable_sort.py
new file mode 100644
index 00000000..63f0304f
--- /dev/null
+++ b/tests/test_cluster_stable_sort.py
@@ -0,0 +1,123 @@
+#write a unit test which tests the r_theta_z sort in cluster.py
+
+import unittest
+import numpy as np
+from ase.build import bulk
+from ase.lattice.cubic import Diamond
+from ase.optimize.precon import PreconLBFGS
+from ase.constraints import ExpCellFilter
+from scipy.linalg import sqrtm
+import ase.io
+import matscipytest
+from matscipy.fracture_mechanics.clusters import diamond, set_groups, set_regions, get_alpha_period,bcc
+
+
+class TestStableSortCluster(matscipytest.MatSciPyTestCase):
+    """
+    Tests of stable sort order when building clusters for fracture
+    testing.
+    """
+    def build_cluster(self,lattice,el,a0,rI,rIII,cutoff,crack_surface,crack_front,shift=np.array([0.0,0.0,0.0])):
+        padding = 3
+        [ax,ay,az] = lattice(el, a0, [1,1,1], crack_surface, crack_front).cell.lengths()
+        n = [2 * (int(np.ceil(((rIII + 2*cutoff)/ ax)))+padding)+1, 2 * (int(np.ceil((rIII + 2*cutoff)/ ay))+padding+1), 1]
+        print('N VALS', n)  
+        cryst = lattice(el, a0, n, crack_surface, crack_front, shift=shift)
+        cluster = set_regions(cryst, rI, cutoff, rIII,sort_type='r_theta_z')
+        print(cluster)
+        #ase.io.write(f'{el}_cluster.xyz',cluster)
+        return cluster
+
+        """    def check_theta_close(self,arr1,arr2,tol):
+        #test if two theta arrays are close, including when theta jumps by pi
+        diff_arr = arr2-arr1
+        wrong_indices = np.where(np.abs(diff_arr)>tol)[0]
+        pi_indices = np.where(np.abs((np.abs(diff_arr))-np.pi)<tol)[0]
+        if len(pi_indices)>0:
+            arr1 = np.delete(arr1, pi_indices)
+            arr2 = np.delete(arr2, pi_indices)
+        self.assertArrayAlmostEqual(arr1,arr2,tol=tol)"""
+
+    def stable_sort_test(self,lattice,el,a0,crack_surface,crack_front,shift=np.array([0.0,0.0,0.0])):
+        #set r_I vals between 24 and 64
+        r_I_vals = np.arange(24, 64, 10)
+        #set r_III vals between 44 and 84
+        r_III_vals = np.arange(44, 84, 10)
+        #set cutoff to 6
+        cutoff = 6
+        #set up largest cluster
+        big_cluster = self.build_cluster(lattice,el,a0,r_I_vals[-1],r_III_vals[-1],cutoff,crack_surface,crack_front,shift=shift)
+        big_cluster_pos = big_cluster.positions
+        #get radial big cluster positions, as well as theta and z
+        sx, sy, sz = big_cluster.cell.diagonal()
+        x, y, z_big = big_cluster_pos[:, 0], big_cluster_pos[:, 1], big_cluster_pos[:, 2]
+        cx, cy = sx/2, sy/2
+        big_cluster_r = np.sqrt((x - cx)**2 + (y - cy)**2)
+        big_cluster_theta = np.arctan2(y-cy,x-cx)
+        #big_cluster.new_array('r', np.zeros(len(big_cluster), dtype=float))
+        #big_cluster.new_array('theta',np.zeros(len(big_cluster)),dtype=float)
+        #cluster.new_array('index',np.zeros(len(cluster)),dtype=float)
+        #r_array = big_cluster.arrays['r']
+        #theta_array = big_cluster.arrays['theta']
+        #r_array[:]=big_cluster_r[:]
+        #theta_array[:]=big_cluster_theta[:]
+        #ase.io.write(f'{el}_big_cluster.xyz',big_cluster)
+
+        #now set up each smaller cluster, and check the r, theta, z of the 
+        # first atoms of big_cluster match those of each small cluster
+        for i in range(len(r_I_vals)-1):
+            cluster = self.build_cluster(lattice,el,a0,r_I_vals[i],r_III_vals[i],cutoff,crack_surface,crack_front,shift=shift)
+            small_cluster_pos = cluster.positions
+            #get radial small cluster positions (from x and y)
+            sx, sy, sz = cluster.cell.diagonal()
+            x, y, z_small = small_cluster_pos[:, 0], small_cluster_pos[:, 1], small_cluster_pos[:, 2]
+            cx, cy = sx/2, sy/2
+            small_cluster_r = np.sqrt((x - cx)**2 + (y - cy)**2)
+            small_cluster_theta = np.arctan2(y-cy,x-cx)
+            cluster.new_array('r', np.zeros(len(cluster), dtype=float))
+            cluster.new_array('theta',np.zeros(len(cluster)),dtype=float)
+            #cluster.new_array('index',np.zeros(len(cluster)),dtype=float)
+            #r_array = cluster.arrays['r']
+            #theta_array = cluster.arrays['theta']
+            #r_array[:]=small_cluster_r[:]
+            #theta_array[:]=small_cluster_theta[:]
+            #ase.io.write(f'{el}_small_cluster.xyz',cluster)
+            #assert that the first atoms of the big cluster match those of the small cluster
+            self.assertArrayAlmostEqual(big_cluster_r[:len(cluster)],small_cluster_r[:],tol=1e-5)
+            #print(big_cluster_theta[7],small_cluster_theta[7])
+            self.assertArrayAlmostEqual(big_cluster_theta[:len(cluster)],small_cluster_theta[:],tol=1e-5)
+            self.assertArrayAlmostEqual(z_small,z_big[:len(cluster)],tol=1e-5)
+        
+    def test_stable_sort_silicon(self):
+        #first, test silicon
+        el = 'Si'
+        a0= 5.43
+        crack_surface = np.array([1, 1, 1]) 
+        crack_front = np.array([1, -1, 0])
+        self.stable_sort_test(diamond,el,a0,crack_surface,crack_front)
+
+        crack_surface = np.array([1, 1, 0]) 
+        crack_front = np.array([0,0,1])
+        self.stable_sort_test(diamond,el,a0,crack_surface,crack_front,shift=np.array([0.0,0.25,0.0]))
+
+
+    
+    def test_stable_sort_iron(self):
+        #next, test iron
+        el = 'Fe'
+        a0 = 2.87
+        crack_surface = np.array([1, 0, 0])
+        crack_front = np.array([0, 0, 1])
+        self.stable_sort_test(bcc,el,a0,crack_surface,crack_front)
+
+        crack_surface = np.array([1,1,0])
+        crack_front = np.array([0,0,1])
+        self.stable_sort_test(bcc,el,a0,crack_surface,crack_front)
+
+        crack_surface = np.array([1,1,0])
+        crack_front = np.array([1,-1,0])
+        self.stable_sort_test(bcc,el,a0,crack_surface,crack_front)
+
+        crack_surface = np.array([1,0,0])
+        crack_front = np.array([0,1,1])
+        self.stable_sort_test(bcc,el,a0,crack_surface,crack_front)
\ No newline at end of file
diff --git a/tests/test_cubic_crystal_crack.py b/tests/test_cubic_crystal_crack.py
index 3a063a1d..da45d877 100644
--- a/tests/test_cubic_crystal_crack.py
+++ b/tests/test_cubic_crystal_crack.py
@@ -62,7 +62,7 @@
 from matscipy.elasticity import Voigt_6x6_to_cubic
 from matscipy.fracture_mechanics.crack import CubicCrystalCrack
 from matscipy.fracture_mechanics.crack import \
-    isotropic_modeI_crack_tip_displacement_field
+    isotropic_modeI_crack_tip_displacement_field, isotropic_modeII_crack_tip_displacement_field
 from matscipy.fracture_mechanics.idealbrittlesolid import IdealBrittleSolid
 
 ###
@@ -112,6 +112,7 @@ def test_isotropic_near_field_solution(self):
 
         k = crack.crack.k1g(1.0)
 
+        #test mode 1
         u, v = crack.crack.displacements(r, theta, k)
         ref_u, ref_v = isotropic_modeI_crack_tip_displacement_field(k, C44, nu,
                                                                     r, theta)
@@ -119,6 +120,26 @@ def test_isotropic_near_field_solution(self):
         self.assertArrayAlmostEqual(u, ref_u)
         self.assertArrayAlmostEqual(v, ref_v)
 
+        #test mode 2
+        u, v = crack.crack.displacements(r, theta, kI=0,kII=k)
+        ref_u, ref_v = isotropic_modeII_crack_tip_displacement_field(k, C44, nu,
+                                                                    r, theta)
+
+        self.assertArrayAlmostEqual(u, ref_u)
+        self.assertArrayAlmostEqual(v, ref_v)
+
+        #test mixed mode
+        u, v = crack.crack.displacements(r, theta, kI=k/2,kII=k/2)
+        ref_u_I, ref_v_I = isotropic_modeI_crack_tip_displacement_field(k/2, C44, nu,
+                                                                    r, theta)
+        ref_u_II, ref_v_II = isotropic_modeII_crack_tip_displacement_field(k/2, C44, nu,
+                                                                    r, theta)
+        ref_u = ref_u_I + ref_u_II
+        ref_v = ref_v_I + ref_v_II
+
+        self.assertArrayAlmostEqual(u, ref_u)
+        self.assertArrayAlmostEqual(v, ref_v)
+
     def test_anisotropic_near_field_solution(self):
         """
         Run an atomistic calculation of a harmonic solid and compare to
@@ -180,6 +201,7 @@ def test_anisotropic_near_field_solution(self):
                 self.assertTrue(np.max(residual[mask]) < 0.2)
 
     def test_consistency_of_deformation_gradient_and_stress(self):
+        #test pure mode 1 first
         for C11, C12, C44, surface_energy, k1 in self.materials:
             crack = CubicCrystalCrack([1,0,0], [0,1,0], C11, C12, C44)
             k = crack.k1g(surface_energy)*k1
@@ -203,8 +225,60 @@ def test_consistency_of_deformation_gradient_and_stress(self):
                 self.assertAlmostEqual(eps[0, 0], eps_xx)
                 self.assertAlmostEqual(eps[1, 1], eps_yy)
                 self.assertAlmostEqual(eps[0, 1], eps_xy)
+        
+        #now test mode 2
+        for C11, C12, C44, surface_energy, k1 in self.materials:
+            crack = CubicCrystalCrack([1,0,0], [0,1,0], C11, C12, C44)
+            k = crack.k1g(surface_energy)*k1
+            for i in range(10):
+                x = np.random.uniform(-10, 10)
+                y = np.random.uniform(-10, 10)
+                F = crack.deformation_gradient(x, y, 0, 0, 0, kII=k)
+                eps = (F+F.T)/2-np.eye(2)
+                sig_x, sig_y, sig_xy = crack.stresses(x, y, 0, 0, 0,kII=k)
+                eps_xx = crack.crack.a11*sig_x + \
+                         crack.crack.a12*sig_y + \
+                         crack.crack.a16*sig_xy
+                eps_yy = crack.crack.a12*sig_x + \
+                         crack.crack.a22*sig_y + \
+                         crack.crack.a26*sig_xy
+                # Factor 1/2 because elastic constants and matrix product are
+                # expressed in Voigt notation.
+                eps_xy = (crack.crack.a16*sig_x + \
+                          crack.crack.a26*sig_y + \
+                          crack.crack.a66*sig_xy)/2
+                self.assertAlmostEqual(eps[0, 0], eps_xx)
+                self.assertAlmostEqual(eps[1, 1], eps_yy)
+                self.assertAlmostEqual(eps[0, 1], eps_xy)
+        
+        #now test mixed mode
+        for C11, C12, C44, surface_energy, k1 in self.materials:
+            crack = CubicCrystalCrack([1,0,0], [0,1,0], C11, C12, C44)
+            k = crack.k1g(surface_energy)*k1
+            for i in range(10):
+                x = np.random.uniform(-10, 10)
+                y = np.random.uniform(-10, 10)
+                F = crack.deformation_gradient(x, y, 0, 0, k/2, k/2)
+                eps = (F+F.T)/2-np.eye(2)
+                sig_x, sig_y, sig_xy = crack.stresses(x, y, 0, 0, k/2, k/2)
+                eps_xx = crack.crack.a11*sig_x + \
+                         crack.crack.a12*sig_y + \
+                         crack.crack.a16*sig_xy
+                eps_yy = crack.crack.a12*sig_x + \
+                         crack.crack.a22*sig_y + \
+                         crack.crack.a26*sig_xy
+                # Factor 1/2 because elastic constants and matrix product are
+                # expressed in Voigt notation.
+                eps_xy = (crack.crack.a16*sig_x + \
+                          crack.crack.a26*sig_y + \
+                          crack.crack.a66*sig_xy)/2
+                self.assertAlmostEqual(eps[0, 0], eps_xx)
+                self.assertAlmostEqual(eps[1, 1], eps_yy)
+                self.assertAlmostEqual(eps[0, 1], eps_xy)
+
 
     def test_consistency_of_deformation_gradient_and_displacement(self):
+        #test mode 1 first
         eps = 1e-3
         for C11, C12, C44, surface_energy, k1 in self.materials:
             crack = CubicCrystalCrack([1,0,0], [0,1,0], C11, C12, C44)
@@ -228,7 +302,54 @@ def test_consistency_of_deformation_gradient_and_displacement(self):
                 F_num = np.transpose([[du_dx, du_dy], [dv_dx, dv_dy]])
                 self.assertArrayAlmostEqual(F, F_num, tol=1e-4)
 
+        #now test mode 2 
+        for C11, C12, C44, surface_energy, k1 in self.materials:
+            crack = CubicCrystalCrack([1,0,0], [0,1,0], C11, C12, C44)
+            k = crack.k1g(surface_energy)*k1
+            for i in range(10):
+                x = i+1
+                y = i+1
+                F = crack.deformation_gradient(x, y, 0, 0, 0, k)
+                # Finite difference approximation of deformation gradient
+                u, v = crack.displacements(x, y, 0, 0, 0, k)
+                ux0, vx0 = crack.displacements(x-eps, y, 0, 0, 0, k)
+                uy0, vy0 = crack.displacements(x, y-eps, 0, 0, 0, k)
+                ux1, vx1 = crack.displacements(x+eps, y, 0, 0, 0, k)
+                uy1, vy1 = crack.displacements(x, y+eps, 0, 0, 0, k)
+                du_dx = (ux1-ux0)/(2*eps)
+                du_dy = (uy1-uy0)/(2*eps)
+                dv_dx = (vx1-vx0)/(2*eps)
+                dv_dy = (vy1-vy0)/(2*eps)
+                du_dx += np.ones_like(du_dx)
+                dv_dy += np.ones_like(dv_dy)
+                F_num = np.transpose([[du_dx, du_dy], [dv_dx, dv_dy]])
+                self.assertArrayAlmostEqual(F, F_num, tol=1e-4)
+        
+        #now test mixed mode
+        for C11, C12, C44, surface_energy, k1 in self.materials:
+            crack = CubicCrystalCrack([1,0,0], [0,1,0], C11, C12, C44)
+            k = crack.k1g(surface_energy)*k1
+            for i in range(10):
+                x = i+1
+                y = i+1
+                F = crack.deformation_gradient(x, y, 0, 0, k/2, k/2)
+                # Finite difference approximation of deformation gradient
+                u, v = crack.displacements(x, y, 0, 0, k/2, k/2)
+                ux0, vx0 = crack.displacements(x-eps, y, 0, 0, k/2, k/2)
+                uy0, vy0 = crack.displacements(x, y-eps, 0, 0, k/2, k/2)
+                ux1, vx1 = crack.displacements(x+eps, y, 0, 0, k/2, k/2)
+                uy1, vy1 = crack.displacements(x, y+eps, 0, 0, k/2, k/2)
+                du_dx = (ux1-ux0)/(2*eps)
+                du_dy = (uy1-uy0)/(2*eps)
+                dv_dx = (vx1-vx0)/(2*eps)
+                dv_dy = (vy1-vy0)/(2*eps)
+                du_dx += np.ones_like(du_dx)
+                dv_dy += np.ones_like(dv_dy)
+                F_num = np.transpose([[du_dx, du_dy], [dv_dx, dv_dy]])
+                self.assertArrayAlmostEqual(F, F_num, tol=1e-4)
+
     def test_elastostatics(self):
+        #check for pure mode I
         eps = 1e-3
         for C11, C12, C44, surface_energy, k1 in self.materials:
             crack = CubicCrystalCrack([1,0,0], [0,1,0], C11, C12, C44)
@@ -247,6 +368,26 @@ def test_elastostatics(self):
                 # equilibrium)
                 self.assertAlmostEqual(divsx, 0.0, places=3)
                 self.assertAlmostEqual(divsy, 0.0, places=3)
+                eps = 1e-3
+
+        #also check for mixed mode
+        for C11, C12, C44, surface_energy, k1 in self.materials:
+            crack = CubicCrystalCrack([1,0,0], [0,1,0], C11, C12, C44)
+            k = crack.k1g(surface_energy)*k1
+            for i in range(10):
+                x = i+1
+                y = i+1
+                # Finite difference approximation of the stress divergence
+                sxx0x, syy0x, sxy0x = crack.stresses(x-eps, y, 0, 0, k/3,(2*k)/3)
+                sxx0y, syy0y, sxy0y = crack.stresses(x, y-eps, 0, 0, k/3,(2*k)/3)
+                sxx1x, syy1x, sxy1x = crack.stresses(x+eps, y, 0, 0, k/3,(2*k)/3)
+                sxx1y, syy1y, sxy1y = crack.stresses(x, y+eps, 0, 0, k/3,(2*k)/3)
+                divsx = (sxx1x-sxx0x)/(2*eps) + (sxy1y-sxy0y)/(2*eps)
+                divsy = (sxy1x-sxy0x)/(2*eps) + (syy1y-syy0y)/(2*eps)
+                # Check that divergence of stress is zero (elastostatic
+                # equilibrium)
+                self.assertAlmostEqual(divsx, 0.0, places=3)
+                self.assertAlmostEqual(divsy, 0.0, places=3)
 
     def test_J_integral(self):
         if quadrature is None: