Repository: libAtoms/matscipy (language: python)

Issue:
WIP: Bring adaptive-cont up-to-date with master and add Cauchy-Born shift and Surface Reconstruction code



Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index 27f88c06..dbbf5b3d 100644
--- a/.gitignore
+++ b/.gitignore
@@ -64,3 +64,5 @@ docs/_build/
 
 #MAC DS_Store
 .DS_Store
+.eggs/README.txt
+.mesonpy-native-file.ini
diff --git a/examples/fracture_mechanics/sinclair_crack/params_Fe.py b/examples/fracture_mechanics/sinclair_crack/params_Fe.py
new file mode 100644
index 00000000..184670b1
--- /dev/null
+++ b/examples/fracture_mechanics/sinclair_crack/params_Fe.py
@@ -0,0 +1,57 @@
+import ase.io
+import numpy as np
+from matscipy.fracture_mechanics.clusters import bcc, set_regions
+
+# Interaction potential
+#from ase.calculators.lammpslib import LAMMPSlib
+#cmds = ["pair_style eam/fs", "pair_coeff * * Fe_2_eam.fs Fe"]
+#calc = LAMMPSlib(lmpcmds=cmds, log_file='test.log')
+#from matscipy.calculators.eam import EAM
+#calc = EAM("eam/pot/Fe_2_eam.fs")
+from pyjulip import ACE1
+calc = ACE1("ace/pot/Fe_ace_N4_D18_R5.json")
+
+# material properties
+el = 'Fe'
+elastic_symmetry = 'cubic'
+from matscipy.dislocation import get_elastic_constants
+a0, C11, C12, C44 = get_elastic_constants(calculator=calc, symbol=el)
+
+# crack system
+r_I = 30.0
+r_III = 60.0
+cutoff = 6.0 # for region 4, beyond r_III
+crack_surface = [1, 0, 0]
+crack_front = [0, 0, 1]
+
+# crack surface energy
+#from matscipy.fracture_mechanics.clusters import find_surface_energy
+from matscipy.surface import find_surface_energy
+surface_energy = find_surface_energy(el,calc,a0,'bcc100',unit='0.1J/m^2')[0]
+
+# simulation control
+fmax = 1e-4 #5e-4
+ds = 0.001
+opt_method = 'ode12r' # currently only for flex1 & flex2, arc_cont uses krylov
+#precon = True # turn on precon for arc_cont using ode12r
+max_opt_steps = 500 # number of ode12r/krylov/sg steps for flex1/flex2
+max_arc_steps = 10 # number of corrector steps in the arc-continuation
+nsteps = 1000 # number of arc continuation steps
+circular_regions = True
+extended_far_field = True # set to True to inslude region III
+# restart tags
+#continuation = True
+#restart_index = 0
+
+# set circular regions I-II-III-IV
+if circular_regions:
+    n = [2 * int((r_III + cutoff)/ a0), 2 * int((r_III + cutoff)/ a0) - 1, 1]
+    print('n=', n)
+
+    # Setup crack system and regions I, II, III and IV
+    cryst = bcc(el, a0, n, crack_surface, crack_front)
+    cluster = set_regions(cryst, r_I, cutoff, r_III)  # carve circular cluster
+
+# Write initial configs
+ase.io.write('cryst.cfg', cryst)
+ase.io.write('cluster.cfg', cluster)
diff --git a/matscipy/cauchy_born.py b/matscipy/cauchy_born.py
new file mode 100644
index 00000000..33b11942
--- /dev/null
+++ b/matscipy/cauchy_born.py
@@ -0,0 +1,1370 @@
+import numpy as np
+from ase.lattice.cubic import Diamond
+
+from ase.optimize import LBFGS
+
+from scipy.linalg import polar, sqrtm
+
+from itertools import permutations
+
+from scipy.stats.qmc import LatinHypercube, scale
+
+import ase
+class RegressionModel:
+    """Simple regression model object wrapper for np.linalg.lstsq for
+    predicting shifts."""
+
+    def fit(self, phi, values):
+        """Fit a simple least-squares regression model to data. Saves output
+        vector to self.model.
+
+        Parameters
+        ----------
+        phi : array_like, 2D
+            Design matrix of shape [number of data points, number of basis
+            functions]. Contains each basis function evaluated at
+            each data point.
+        values : array_like
+            value of function to be fitted evaluated at each data point
+        """
+        self.model = np.linalg.lstsq(phi, values,rcond=None)[0]
+
+    def predict(self, phi):
+        """Evaluate the simple least-squares regression model.
+
+        Parameters
+        ----------
+        phi : array_like, 2D
+            Design matrix of shape [number of data points,
+            number of basis functions].
+            Contains each basis function evaluated at each data point.
+
+        Returns
+        -------
+        predictions : array_like
+            1D vector of least squares model predictions for simple model.
+        """
+
+        return phi @ self.model
+
+
+class CubicCauchyBorn:
+    """Corrector for the cauchy-born prediction of atomistic positions in
+    multi-lattices subject to continuum strain This model exploits the symmetry
+    triad down the 111 axis that all cubic crystals possess by making models
+    that only fit a single component of the shift-correction vector.
+
+    The other 2 components can be obtained by calling the same model
+    with strain states that are rotated about the 111 axis.
+    """
+
+    def __init__(self, el, a0, calc, lattice=Diamond):
+        """Parameters
+        ----------
+        el : string
+            ASE chemical element that constitutes lattice
+        a0 : float
+            Lattice constant for cubic lattice
+        calc : ASE calculator object
+            ASE calculator.
+        lattice : ASE lattice builder function
+            ASE function from ase.lattice.cubic
+        """
+        self.el = el
+        self.a0 = a0
+        self.calc = calc
+        self.lattice = lattice
+        self.regression_Fitted = False  # status of model fit
+        self.lattice1mask = None  # mask for the lattice 1 atoms
+        self.lattice2mask = None  # mask for the lattice 2 atoms
+
+    def set_sublattices(self, atoms, A):
+        """Apply a small strain to all atoms in the supplied atoms structure
+        and determine which atoms belong to which sublattice using forces. NOTE
+        as this method is based on forces, it does not work in cases where atom
+        forces are high already, This means it should be applied to an
+        unstrained version of the atoms structure and any sublattices assigned
+        adjacent to free surfaces cannot be trusted. This function updates
+        self.lattice1mask and self.lattice2mask with the measured sublattices.
+
+        Parameters
+        ----------
+        atoms : ASE atoms object
+            atoms structure with atoms on seperate
+            sublattices that need assigning
+        A : 3x3 numpy array, floats
+            rotation matrix of the form [x^T,y^T,z^T],
+            where x^T,y^T,z^T are normalised column vectors
+            of the lab frame directions expressed
+            in terms of the crystal lattice directions.
+        """
+
+        # currently, the idea for this is to apply a small strain to atoms
+        # measure the force difference, and set all
+        # the forces which are positive to sublattice
+        # 1, and all those which are negative to sublattice 2.
+
+        # generate U (right hand stretch tensor) to apply
+        U_voigt = np.array([1.001, 1.003, 1.002, 0.006, 0.002, 0.004])
+        U = np.zeros([3, 3])
+        U[0, 0] = U_voigt[0]
+        U[1, 1] = U_voigt[1]
+        U[2, 2] = U_voigt[2]
+        U[1, 2], U[2, 1] = U_voigt[3], U_voigt[3]
+        U[0, 2], U[2, 0] = U_voigt[4], U_voigt[4]
+        U[0, 1], U[1, 0] = U_voigt[5], U_voigt[5]
+        atoms_copy = atoms.copy()
+
+        # generate a copy of the atoms and stretch the cell, scaling atoms.
+        # get forces before and after and compare these
+        atoms_copy.calc = self.calc
+        f_before = atoms_copy.get_forces()
+        cell = atoms_copy.get_cell()
+        scaled_cell = U @ cell
+        atoms_copy.set_cell(scaled_cell, scale_atoms=True)
+        f_after = atoms_copy.get_forces()
+
+        force_diff = f_after - f_before
+
+        # transform force_diff to lattice frame
+        force_diff_lattice = np.zeros_like(force_diff)
+        for i in range(len(atoms)):
+            force_diff_lattice[i, :] = A @ force_diff[i, :]
+
+        # generate the lattice masks, checking to
+        # ensure that if by chance no shift occurs in one
+        # direction due to the applied strains
+        # that the other directions are used instead and things don't break
+        lattice1mask = force_diff_lattice[:, 0] > 0
+        lattice2mask = force_diff_lattice[:, 0] < 0
+        # if they're all false or true, try the y component
+        if all(element == lattice1mask[0] for element in lattice1mask):
+            lattice1mask = force_diff_lattice[:, 1] > 0
+            lattice2mask = force_diff_lattice[:, 1] < 0
+            # if they're all false, try the z component
+            if all(element == lattice1mask[0] for element in lattice1mask):
+                lattice1mask = force_diff_lattice[:, 2] > 0
+                lattice2mask = force_diff_lattice[:, 2] < 0
+
+        # set the masks for which atoms are in each of the sublattices
+        self.lattice1mask = lattice1mask
+        self.lattice2mask = lattice2mask
+
+    def fit_taylor(self, de=1e-4):
+        """Fit a simple taylor-expansion type model to predict the cauchy-born
+        shift corrector at a given applied strain using finite differences.
+        Sets self.grad_f to the vector of first derivatives of the first shift
+        component with each of the 6 strain components and self.hess_f to the
+        hessian of the first shift component with the strain components. NOTE
+        this model is generally worse than the model provided later in
+        initial_regression_fit.
+
+        Parameters
+        ----------
+        de : float
+            tolerance for finite difference gradient calculations
+        """
+        calc = self.calc
+
+        def f_gl(E_vec, calc, unitcell):
+            """Function that calculates the first component of the cauchy-born
+            shift corrector for a given green-lagrange strain tensor.
+
+            Parameters
+            ----------
+            E_vec : array_like
+                Green-Lagrange strain state to apply in voigt notation
+            calc : ASE calculator object
+                ASE calculator to use to find shift
+            unitcell : ASE atoms object
+                Single unit cell of crystal to apply
+                strain vector to and to find shift
+
+            Returns
+            -------
+            shift_diff[0] : float
+                first component of cauchy-born shift correction
+                for applied strain state
+            """
+            # green lagrange version of function
+            # turn E into matrix
+            E = np.zeros([3, 3])
+            E[0, 0] = E_vec[0]
+            E[1, 1] = E_vec[1]
+            E[2, 2] = E_vec[2]
+            E[2, 1], E[1, 2] = E_vec[3], E_vec[3]
+            E[2, 0], E[0, 2] = E_vec[4], E_vec[4]
+            E[0, 1], E[1, 0] = E_vec[5], E_vec[5]
+            # get U^2
+            Usqr = 2 * E + np.eye(3)
+            # square root matrix to get U
+            U = sqrtm(Usqr, disp=True)
+
+            # this is just the symmetric stretch tensor, exactly what we need.
+            x = U
+
+            initial_shift = np.zeros([3])
+            relaxed_shift = np.zeros([3])
+
+            # build cell
+            unitcell_copy = unitcell.copy()
+            unitcell_copy.calc = self.calc
+            cell = unitcell_copy.get_cell()
+            cell_rescale = np.transpose(x[:, :] @ cell)
+            unitcell_copy.set_cell(cell_rescale, scale_atoms=True)
+            initial_shift[:] = unitcell_copy.get_positions()[1] - unitcell_copy.get_positions()[0]
+
+            # relax cell
+            opt = LBFGS(unitcell_copy, logfile=None)
+            opt.run(fmax=1e-10)
+            relaxed_shift[:] = unitcell_copy.get_positions()[1] - unitcell_copy.get_positions()[0]
+
+            # get shift
+            shift_diff = relaxed_shift - initial_shift
+            # print(shift_diff)
+            # THIS CAN BE SPED UP
+            # TECHNICALLY, EACH CALL OF THIS FUNCTION CAN BE USED
+            # TO COMPUTE 3 DERIVATIVES RATHER THAN JUST 1
+            # IF YOU EXPLOIT THE SYMMETRY OF CUBIC SYSTEMS
+            # WOULD MAKE THINGS 3x FASTER, BUT WOULD REQUIRE
+            # *SOME* CODING
+            return shift_diff[0]  # return 1 component
+
+        def find_grad(eps, i, *args):
+            """"Compute gradient of first shift component by varying a single
+            component of eps using finite differences.
+
+            Parameters
+            ----------
+            eps : array_like
+                Green-Lagrange strain state to vary components of in
+                Voigt notation
+            i : integer
+                Component of eps to find gradient of shift with respect to
+            *args: additional arguments accepted by shift evaluation function.
+
+            Returns
+            -------
+            $$\frac{\\partial nu_{1}}{\\partial {\\epsilon_{i}}}$$ : float
+            """
+            eps_up = eps.copy()
+            eps_down = eps.copy()
+            eps_up[i] += self.de
+            eps_down[i] -= self.de
+            return (f_gl(eps_up, *args) - f_gl(eps_down, *args)) / (2 * self.de)
+
+        def find_second_deriv(eps, i, j, *args):
+            """"Compute mixed 2nd derivative of first shift component by
+            varying two components of eps using finite differences. Essentially
+            computes a single term in the Hessian.
+
+            Parameters
+            ----------
+            eps : array_like
+                Green-Lagrange strain state to vary components
+                of in Voigt notation
+            i : integer
+                First component of eps to find
+                gradient of shift with respect to
+            j : integer
+                Second component of eps to find
+                gradient of shift with respect to
+            *args: additional arguments accepted by shift evaluation function.
+
+            Returns
+            -------
+            $$\frac{\\partial^{2} nu_{1}}{\\partial \\epsilon_{i}\\partial \\epsilon_{j}}$$ : float
+            """
+            eps_up = eps.copy()
+            eps_down = eps.copy()
+            eps_up[j] += self.de
+            eps_down[j] -= self.de
+            return (find_grad(eps_up, i, calc, unitcell) -
+                    find_grad(eps_down, i, calc, unitcell)) / (2 * self.de)
+
+        self.de = de
+        unitcell = self.lattice(
+            self.el, latticeconstant=self.a0, size=[1, 1, 1])
+        unitcell.set_pbc([True, True, True])
+
+        # compute all the terms we need:
+        grad_f = np.zeros([6])
+        hess_f = np.zeros([6, 6])
+        eps_vec = np.zeros([6])
+        for i in range(6):
+            grad_f[i] = find_grad(eps_vec, i, calc, unitcell)
+            # finds the gradient of the gradient of f with e[i] with respect to
+            # each strain component e[j]
+            for j in range(6 - i):
+                hess_f[i, j] = find_second_deriv(eps_vec, i, j, calc, unitcell)
+                # include relevant symmetries
+                hess_f[j, i] = hess_f[i, j]
+
+        self.grad_f = grad_f
+        self.hess_f = hess_f
+
+    def save_taylor(self):
+        """Save the results of the taylor expansion in to numpy readable files
+        rather than re-calculating each time."""
+        try:
+            np.save('grad_f.npy', self.grad_f)
+            np.save('hess_f.npy', self.hess_f)
+        except AttributeError:
+            print('Need to fit a model with fit taylor before it can be saved!')
+
+    def load_taylor(self):
+        """Read the results of the taylor expansion rather than re-calculating
+        each time."""
+        try:
+            self.grad_f = np.load('grad_f.npy')
+            self.hess_f = np.load('hess_f.npy')
+        except FileNotFoundError:
+            print('No files found containing model parameters - make sure to call save_taylor() to save the model.')
+
+    def tensor_field_3D_from_atoms(
+            self,
+            atoms,
+            func,
+            coordinates='cart3D',
+            cell=None,
+            *args,
+            **kwargs):
+        """Generate a 3D tensor field (F or E) for an atoms
+        object in a given coordinate system, using a provided
+        function
+
+        Parameters
+        ----------
+        atoms : ASE atoms object
+            Atoms object where atoms sit at the coordinates needed
+            to generate the tensor field from func.
+        func : Function
+            Function to calculate the tensor field
+            from a set of coordinates. The system of coordinates that
+            the function accepts is given by 'coordinates' and the
+            function also accepts anything extra specified in *args
+            and **kwargs. Function must return the field for the atoms
+            in form of a numpy array shape [natoms,ndims,ndims],
+            where output is specified in cartesian coordinates.
+        coordinates: string
+            The coordinates system that func accepts. Must
+            be 'cylind2D', 'cylind3D', 'spherical', 'cart2D', 'cart3D'
+        cell : array_like, optional
+            An optional argument allowing the user to specify
+            the cell parameters used to define the coordinate system.
+
+        Returns
+        -------
+        field3D : array_like
+            The 3D tensor field for all
+            the atoms in the system, expressed in the lab frame.
+        """
+        natoms = len(atoms)
+        if cell is None:
+            sx, sy, sz = atoms.cell.diagonal()
+        else:
+            [sx, sy, sz] = cell
+
+        if coordinates == 'cylind2D':
+            x, y = atoms.positions[:, 0], atoms.positions[:, 1]
+            cx, cy = sx / 2, sy / 2
+            r = np.sqrt((x - cx)**2 + (y - cy)**2)
+            theta = np.arctan2((y - cy), (x - cx))
+            TF_2D = func(r, theta, *args, **kwargs)
+            # pad F into a 3x3 matrix
+            TF_3D = (np.zeros([natoms, 3, 3]))
+            TF_3D[:, 0:2, 0:2] = TF_2D[:, :, :]
+            TF_3D[:, 2, 2] = 1
+
+        elif coordinates == 'cylind3D':
+            x, y, z = atoms.positions[:,
+                                      0], atoms.positions[:, 1], atoms.positions[:, 2]
+            cx, cy, cz = sx / 2, sy / 2, sz / 2
+            r = np.sqrt((x - cx)**2 + (y - cy)**2)
+            theta = np.arctan2((y - cy), (x - cx))
+            TF_3D = func(r, theta, z, *args, **kwargs)
+
+        elif coordinates == 'spherical':
+            raise NotImplementedError
+
+        elif coordinates == 'cart2D':
+            x, y = atoms.positions[:, 0], atoms.positions[:, 1]
+            TF_2D = func(x, y, *args, **kwargs)
+            TF_3D = (np.zeros([natoms, 3, 3]))
+            TF_3D[:, 0:2, 0:2] = TF_2D[:, :, :]
+            TF_3D[:, 2, 2] = 1
+
+        elif coordinates == 'cart3D':
+            x, y, z = atoms.positions[:,
+                                      0], atoms.positions[:, 1], atoms.positions[:, 2]
+            TF_3D = func(x, y, z, *args, **kwargs)
+
+        else:
+            print(
+                'Coordinate system should be cart2D, cart3D, cylind2D, cylind3D, spherical')
+            raise NotImplementedError
+
+        return TF_3D
+
+    def evaluate_F_or_E(
+            self,
+            A,
+            atoms,
+            F_func=None,
+            E_func=None,
+            cell=None,
+            coordinates='cart3D',
+            *args,
+            **kwargs):
+        """Get the deformation gradient tensor field or Green-Lagrange strain
+        tensor field for a system of atoms in the lab frame, depending on which
+        function is provided. From these, find the Green-Lagrange strain tensor
+        field in the lattice frame.
+
+        Parameters
+        ----------
+        A : 3x3 numpy array, floats
+            rotation matrix of the form [x^T,y^T,z^T],
+            where x^T,y^T,z^T are normalised column vectors
+            of the lab frame directions expressed
+            in terms of the crystal lattice directions.
+        atoms : ASE atoms object
+            Atoms object where atoms sit at the coordinates needed to
+            generate the deformation gradient tensor field from F_func.
+        F_func : Function
+            Function to calculate the deformation gradient tensor field
+            from a set of atomic coordinates. The system of coordinates
+            that the function accepts is given by 'coordinates' and the
+            function also accepts anything extra specified in *args and
+            **kwargs. Function must return the deformation gradient field
+            for the atoms in form of a numpy array shape [natoms,ndims,ndims],
+            where F is specified in cartesian coordinates. The returned
+            tensor field must be in the lab frame Optional, but one of
+            F_func or E_func must be provided.
+        E_func : Function
+            Function to calculate the Green-Lagrange tensor field
+            from a set of atomic coordinates. The system of coordinates
+            that the function accepts is given by 'coordinates'
+            and the function also accepts anything extra specified
+            in *args and **kwargs. Function must return the
+            Green-Lagrange tensor field for the atoms in form of a
+            numpy array shape [natoms,ndims,ndims], where E is
+            specified in cartesian coordinates, and the returned
+            tensor field is in the lab frame. Optional, but one of
+            F_func or E_func must be provided.
+        cell : array_like, optional
+            An optional argument allowing the user to specify
+            the cell dimensions which will be used to find the
+            centre of the cell. If not provided, the centre
+            of the cell will be found from the atoms object provided.
+        coordinates: string
+            Specifies the coordinate system that the function
+            (F_func or E_func) provided accepts. The coordinates
+            of the atoms in atoms will be converted to this
+            coordinate system before being passed to the function
+            as a vector. Default is 3D cartesian coordinates.
+            Other options are cart2D, cart3D, cylind2D, cylind3D,
+            spherical
+
+
+
+        Returns
+        -------
+        E : array_like
+            The 3D Green-Lagrange tensor field for all the
+            atoms in the system, expressed in the lattice frame.
+        R : array_like
+            The rotation tensor component of the deformation
+            gradient tensor field applied to the system,
+            according to the polar decomposition F = RU.
+            See
+            https://en.wikipedia.org/wiki/Finite_strain_theory#Polar_decomposition_of_the_deformation_gradient_tensor.
+        """
+        if F_func is not None:
+            if E_func is not None:
+                print('Need to only provide one of E or F, not both')
+                raise ValueError
+            return self.evaluate_F(
+                A,
+                atoms,
+                F_func,
+                cell=cell,
+                coordinates=coordinates,
+                *args,
+                **kwargs)
+        elif E_func is not None:
+            return self.evaluate_E(
+                A,
+                atoms,
+                E_func,
+                cell=cell,
+                coordinates=coordinates,
+                *args,
+                **kwargs)
+        else:
+            print('Error, please provide one of E or F')
+            raise ValueError
+
+    def evaluate_F(
+            self,
+            A,
+            atoms,
+            F_func,
+            cell=None,
+            coordinates='cart3D',
+            *args,
+            **kwargs):
+        """Get the deformation gradient tensor field for a system of atoms in the lab frame.
+        From this, find the Green-Lagrange strain tensor field in the lattice frame.
+
+        Parameters
+        ----------
+        A : 3x3 numpy array, floats
+            rotation matrix of the form [x^T,y^T,z^T],
+            where x^T,y^T,z^T are normalised column vectors
+            of the lab frame directions expressed
+            in terms of the crystal lattice directions.
+        atoms : ASE atoms object
+            Atoms object where atoms sit at the coordinates needed to
+            generate the deformation gradient tensor field from F_func.
+        F_func : Function
+            Function to calculate the deformation gradient tensor field
+            from a set of atomic coordinates. The system of coordinates
+            that the function accepts is given by 'coordinates' and the
+            function also accepts anything extra specified in *args and
+            **kwargs. Function must return the deformation gradient field
+            for the atoms in form of a numpy array shape [natoms,ndims,ndims],
+            where F is specified in cartesian coordinates. The returned
+            tensor field must be in the lab frame.
+        cell : array_like, optional
+            An optional argument allowing the user to specify
+            the cell dimensions which will be used to find the
+            centre of the cell. If not provided, the centre
+            of the cell will be found from the atoms object provided.
+        coordinates: string
+            Specifies the coordinate system that the function
+            (F_func or E_func) provided accepts. The coordinates
+            of the atoms in atoms will be converted to this
+            coordinate system before being passed to the function
+            as a vector. Default is 3D cartesian coordinates.
+            Other options are cart2D, cart3D, cylind2D, cylind3D,
+            spherical
+        Returns
+        -------
+        E : array_like
+            The 3D Green-Lagrange tensor field for all the
+            atoms in the system, expressed in the lattice frame.
+        R : array_like
+            The rotation tensor component of the deformation
+            gradient tensor field applied to the system,
+            according to the polar decomposition F = RU.
+            See
+            https://en.wikipedia.org/wiki/Finite_strain_theory#Polar_decomposition_of_the_deformation_gradient_tensor.
+        """
+
+        natoms = len(atoms)
+        F_3D = self.tensor_field_3D_from_atoms(
+            atoms,
+            func=F_func,
+            cell=cell,
+            coordinates=coordinates,
+            *args,
+            **kwargs)
+
+        # transform F from lab frame to lattice frame,
+        # find right handed stretch tensor U and rotation matrix
+        # R via polar decomposition, then use these to find
+        # green-lagrange strain tensors E.
+
+        Fprime = np.zeros_like(F_3D)
+        R = np.zeros_like(F_3D)
+        U = np.zeros_like(F_3D)
+        E = np.zeros_like(F_3D)
+        for i in range(natoms):
+            Fprime[i, :, :] = A @ F_3D[i, :, :] @ np.transpose(A)
+            R[i, :, :], U[i, :, :] = polar(Fprime[i, :, :])
+            E[i, :, :] = 0.5 * \
+                (U[i, :, :] @ (np.transpose(U[i, :, :])) - np.eye(3))
+        return E, R
+
+    def evaluate_E(
+            self,
+            A,
+            atoms,
+            E_func,
+            cell=None,
+            coordinates='cart3D',
+            *args,
+            **kwargs):
+        """Get the Green-Lagrange strain tensor field for a system of atoms in the lab frame
+        directly from a provided function.
+
+        Parameters
+        ----------
+        A : 3x3 numpy array, floats
+            rotation matrix of the form [x^T,y^T,z^T],
+            where x^T,y^T,z^T are normalised column vectors
+            of the lab frame directions expressed
+            in terms of the crystal lattice directions.
+        atoms : ASE atoms object
+            Atoms object where atoms sit at the coordinates needed to
+            generate the deformation gradient tensor field from F_func.
+        E_func : Function
+            Function to calculate the Green-Lagrange tensor field
+            from a set of atomic coordinates. The system of coordinates
+            that the function accepts is given by 'coordinates'
+            and the function also accepts anything extra specified
+            in *args and **kwargs. Function must return the
+            Green-Lagrange tensor field for the atoms in form of a
+            numpy array shape [natoms,ndims,ndims], where E is
+            specified in cartesian coordinates, and the returned
+            tensor field is in the lab frame.
+        cell : array_like, optional
+            An optional argument allowing the user to specify
+            the cell dimensions which will be used to find the
+            centre of the cell. If not provided, the centre
+            of the cell will be found from the atoms object provided.
+        coordinates: string
+            Specifies the coordinate system that the function
+            (F_func or E_func) provided accepts. The coordinates
+            of the atoms in atoms will be converted to this
+            coordinate system before being passed to the function
+            as a vector. Default is 3D cartesian coordinates.
+            Other options are cart2D, cart3D, cylind2D, cylind3D,
+            spherical
+        Returns
+        -------
+        E : array_like
+            The 3D Green-Lagrange tensor field for all the
+            atoms in the system, expressed in the lattice frame.
+        R : array_like
+            The rotation tensor component of the deformation
+            gradient tensor field applied to the system,
+            according to the polar decomposition F = RU.
+            See
+            https://en.wikipedia.org/wiki/Finite_strain_theory#Polar_decomposition_of_the_deformation_gradient_tensor.
+        """
+        natoms = len(atoms)
+        # in the case that only the E function is provided, find and rotate the
+        # Green-Lagrange strain tensor field directly
+        E_3D_lab = self.tensor_field_3D_from_atoms(
+            atoms,
+            func=E_func,
+            cell=cell,
+            coordinates=coordinates,
+            *args,
+            **kwargs)
+        E = np.zeros_like(E_3D_lab)
+        R = np.zeros_like(E_3D_lab)
+        # for a specified strain, R is just I for all atoms
+        for i in range(natoms):
+            E[i, :, :] = A @ E_3D_lab[i, :, :] @ np.transpose(A)
+            R[i, :, :] = np.eye(3)
+        return E, R
+
+    def predict_shifts(
+            self,
+            A,
+            atoms,
+            F_func=None,
+            E_func=None,
+            coordinates='cart3D',
+            method='regression',
+            *args,
+            **kwargs):
+        """Get the deformation gradient tensor field or Green-Lagrange strain
+        tensor field for a system of atoms in the lab frame, depending on which
+        function is provided. From these, find the Green-Lagrange strain tensor
+        field in the lattice frame.
+
+        Parameters
+        ----------
+        A : 3x3 numpy array, floats
+            rotation matrix of the form [x^T,y^T,z^T], where x^T,y^T,z^T
+            are normalised column vectors of the lab frame directions
+            expressed in terms of the crystal lattice directions.
+        atoms : ASE atoms object
+            Atoms object where atoms sit at the coordinates needed to
+            generate the deformation gradient tensor field from F_func.
+        F_func : Function
+            Function to calculate the deformation gradient tensor field
+            from a set of atomic coordinates. The system of coordinates
+            that the function accepts is given by 'coordinates' and the
+            function also accepts anything extra specified in *args and
+            **kwargs. Function must return the deformation gradient field
+            for the atoms in form of a numpy array shape [natoms,ndims,ndims],
+            where F is specified in cartesian coordinates. The returned
+            tensor field must be in the lab frame Optional, but one of
+            F_func or E_func must be provided.
+        E_func : Function
+            Function to calculate the Green-Lagrange tensor field
+            from a set of atomic coordinates. The system of coordinates
+            that the function accepts is given by 'coordinates' and
+            the function also accepts anything extra specified in *args
+            and **kwargs. Function must return the Green-Lagrange tensor
+            field for the atoms in form of a numpy array shape
+            [natoms,ndims,ndims], where E is specified in cartesian
+            coordinates, and the returned tensor field is in the lab frame.
+             Optional, but one of F_func or E_func must be provided.
+        coordinates : string
+            Specifies the coordinate system that the function
+            (F_func or E_func) provided accepts. The coordinates
+            of the atoms in atoms will be converted to this coordinate
+            system before being passed to the function as a vector.
+            Default is 3D cartesian coordinates. Other options are
+            cart2D, cart3D, cylind2D, cylind3D, spherical
+        method : string
+            Specifies which method to use when evaluating shifts.
+            Options are 'taylor' (for the taylor expansion model)
+            and 'regression' for the regression model.
+            'regression' by default.
+        Returns
+        -------
+        shifts : array-like
+            The cauchy-born corrector shift predictions made by the model
+        """
+        E, R = self.evaluate_F_or_E(
+            A, atoms, F_func=F_func, E_func=E_func,
+            coordinates=coordinates, *args, **kwargs)
+        natoms = len(atoms)
+        # get the cauchy born shifts unrotated
+        shifts_no_rr = self.evaluate_shift_model(E, method=method)
+        shifts = np.zeros_like(shifts_no_rr)
+
+        # rotate the cauchy shifts both by the rotation induced by F
+        # and to get them back into the lab frame
+        for i in range(natoms):
+            shifts[i, :] = np.transpose(A) @ R[i, :, :] @ shifts_no_rr[i, :]
+        return shifts
+
+    def evaluate_shift_model(self, E, method='regression'):
+        """Evaluate a fitted shift model and return the shift predictions for
+        each atom.
+
+        Parameters
+        ----------
+        E : array_like
+            list of Green-Lagrange strain tensors defining the strain state
+            on each atom, of the form [natoms,3,3]
+        method : string
+            Specifies which method to use when evaluating shifts.
+            Options are 'taylor' (for the taylor expansion model)
+            and 'regression' for the regression model. 'regression' by default.
+        Returns
+        -------
+        shifts : array-like
+            The cauchy-born corrector shift predictions made by the model
+            specified in 'method'
+        """
+        # takes E in the form of a list of matrices [natoms, :, :]
+        # define function for evaluating model
+
+        def permutation(strain_vec, perm_shift):
+            """Rotate a Green-Lagrange strain tensor in Voigt notation about
+            the 111 axis a number of times given by perm_shift. This
+            corresponds to permuting the top three components and bottom three
+            components seperately. For example, the vector [1,2,3,4,5,6]
+            becomes [2,3,1,5,6,4] under this operation. Due to the triad 111
+            symmetry in diamond, feeding these rotated strain states to the
+            shift prediction model allows the prediction of the full shift
+            vector.
+
+            Parameters
+            ----------
+            strain_vec : array_like
+                2D array of size [natoms,6], where each row is a Green-Lagrange
+                strain tensor written in Voigt notation.
+            perm_shift : integer
+                The permeutation number wanted, with 1 corresponding
+                to a 120 degree rotation about 111 (these strain
+                vectors give the y component of shift prediction when
+                fed to the model) and 2 corresponding to a 240 degree
+                rotation about 111 (these strain vectors give the z
+                component of shift prediction when fed to the model)
+
+            Returns
+            -------
+            strain_vec_perm : array_like
+                Rotated version of strain_vec
+            """
+            strain_vec_perm = np.zeros_like(strain_vec)
+            for i in range(3):
+                perm_num = (i - perm_shift) % 3
+                strain_vec_perm[:, perm_num] = strain_vec[:, i]
+                strain_vec_perm[:, perm_num + 3] = strain_vec[:, i + 3]
+            return strain_vec_perm
+
+        def evaluate_shift_regression(eps):
+            """Evaluate the regression model a given set of Green-Lagrange
+            strain vectors.
+
+            Parameters
+            ----------
+            eps : array_like
+                2D array of size [natoms,6], where each row is a Green-Lagrange
+                strain tensor written in Voigt notation.
+
+            Returns
+            -------
+            predictions : array_like
+                prediction of Cauchy-Born shift correctors from
+                regression model.
+            """
+
+            # Get the rotated versions of the strain tensors
+            epsx = eps
+            epsy = permutation(eps, 1)
+            epsz = permutation(eps, 2)
+            predictions = np.zeros([np.shape(eps)[0], 3])
+
+            # Predict the shifts using the fitted regression model
+            predictions[:, 0] = self.RM.predict(
+                self.basis_function_evaluation(epsx))
+            predictions[:, 1] = self.RM.predict(
+                self.basis_function_evaluation(epsy))
+            predictions[:, 2] = self.RM.predict(
+                self.basis_function_evaluation(epsz))
+            return predictions
+
+        def eval_tay_model(self, eps):
+            """Evaluate the taylor model a given set of Green-Lagrange strain
+            vectors.
+
+            Parameters
+            ----------
+            eps : array_like
+                2D array of size [natoms,6], where each row is a Green-Lagrange
+                strain tensor written in Voigt notation.
+
+            Returns
+            -------
+            predictions : array_like
+                prediction of Cauchy-Born shift correctors from taylor model.
+            """
+            grad_f = self.grad_f
+            hess_f = self.hess_f
+            # takes E as a 2D array of voigt vectors, with number of
+            # rows=number atoms
+            # Note that the permutation is computed from observing
+            # how the matrix components permute
+            # when transforming the strain matrix by the 120 degree
+            # rotation matrices about [111],
+            # such that the axis are permuted in the order
+            # 1 ---> 2 ---> 3
+
+            def eval_eps(eps_vec, grad_f, hess_f):
+                """Function to explicitly evaluate the taylor expansion model.
+
+                Parameters
+                ----------
+                eps_vec : array_like
+                    2D array of Green-Lagrange strain states written in
+                    Voigt notation.
+
+                Returns
+                -------
+                prediction : float
+                    Prediction of shift components from strain vectors
+                    using taylor expansion model.
+                """
+                term_1 = grad_f @ eps_vec
+                term_2 = (1 / 2) * (np.transpose(eps_vec) @ hess_f @ eps_vec)
+
+                return term_1 + term_2
+
+            epsx = eps
+            epsy = permutation(eps, 1)
+            epsz = permutation(eps, 2)
+            predictions = np.zeros([np.shape(eps)[0], 3])
+            for i in range(np.shape(eps)[0]):
+                predictions[i, 0] = eval_eps(epsx[i, :], grad_f, hess_f)
+                predictions[i, 1] = eval_eps(epsy[i, :], grad_f, hess_f)
+                predictions[i, 2] = eval_eps(epsz[i, :], grad_f, hess_f)
+
+            return predictions
+
+        # turn E into voigt vector
+        E_voigt = np.zeros([np.shape(E)[0], 6])
+        E_voigt[:, 0] = E[:, 0, 0]
+        E_voigt[:, 1] = E[:, 1, 1]
+        E_voigt[:, 2] = E[:, 2, 2]
+        E_voigt[:, 3] = E[:, 1, 2]
+        E_voigt[:, 4] = E[:, 0, 2]
+        E_voigt[:, 5] = E[:, 0, 1]
+
+        # return the shift vectors
+        if method == 'taylor':
+            return eval_tay_model(self, E_voigt)
+        elif method == 'regression':
+            return evaluate_shift_regression(E_voigt)
+        else:
+            print(
+                'Error! Can only predict cauchy born shift with implemented methods - "taylor" or "regression"')
+            raise NotImplementedError
+
+    def apply_shifts(self, atoms, shifts, mask=None):
+        """Apply predicted Cauchy-Born corrector shifts to atoms object with an
+        optional mask.
+
+        Parameters
+        --------
+        atoms : ASE atoms object
+            Atoms object to apply Cauchy-Born corrector shifts too.
+        shifts : array_like
+            Cauchy-Born corrector shifts to apply
+        mask : array_like, bool, optional
+            Optional mask which specifies which atoms to apply shifts to.
+            All atoms are shifted by default.
+        """
+        # mask covers atoms that we want to be shifted
+        if mask is not None:
+            lattice1mask = np.logical_and(self.lattice1mask, mask)
+            lattice2mask = np.logical_and(self.lattice2mask, mask)
+        else:
+            lattice1mask = self.lattice1mask
+            lattice2mask = self.lattice2mask
+        atoms.positions[lattice1mask] += -0.5 * shifts[lattice1mask]
+        atoms.positions[lattice2mask] += 0.5 * shifts[lattice2mask]
+
+    def get_cb_error(self, atoms, mask=None, forces=None):
+        """Get an approximate quantifier of the error in an atoms object due to
+        poor or not implemented Cauchy-Born corrector shifts.
+
+        Parameters
+        ----------
+        atoms : ASE atoms object
+            ASE atoms for which to quantify the Cauchy-Born
+            corrector related error.
+        mask : array_like, bool
+            Mask for which atoms to include in the error
+            quantification process.
+        forces : array_like
+            Forces on atoms. Can be provided in order
+            to stop this function from carrying out a
+            potentially expensive atoms.get_forces() operation.
+        Returns
+        -------
+        predictions : cb_err
+            A quantification of the Cauchy-Born shift
+            related error in an atoms object, which
+            can be used as an objective to minimise.
+        """
+
+        atoms.set_calculator = self.calc
+        if forces is None:
+            forces = atoms.get_forces()
+
+        if mask is not None:
+            lattice1mask = np.logical_and(self.lattice1mask, mask)
+            lattice2mask = np.logical_and(self.lattice2mask, mask)
+        else:
+            lattice1mask = self.lattice1mask
+            lattice2mask = self.lattice2mask
+
+        lattice1_force_vec = forces[lattice1mask]
+        lattice2_force_vec = forces[lattice2mask]
+        lattice1_force = np.sum(lattice1_force_vec)
+        lattice2_force = np.sum(lattice2_force_vec)
+
+        # find the total forces on each sublattice and subtract them
+        cb_err = np.abs(lattice1_force - lattice2_force)
+
+        return cb_err  # we should be trying to minimise this
+
+    def initial_regression_fit(self, initial_samples=10):
+        """Perform an initial fit of the regression model with some number of
+        initial samples.
+
+        Parameters
+        ----------
+        initial_samples : int
+            The number of Latin Hypercube samples to fit the initial
+            Cauchy-Born regression model to.
+        """
+        # initialise the regression model
+        self.RM = RegressionModel()
+        # build the initial simple dataset to begin fitting with
+        # take samples in the plus/minus 1% strain hypercube in strain space
+        sampler = LatinHypercube(d=6)
+        samples = sampler.random(n=initial_samples)
+        lbounds = [-0.01, -0.01, -0.01, -0.01, -0.01, -0.01]
+        ubounds = [0.01, 0.01, 0.01, 0.01, 0.01, 0.01]
+        samples = scale(samples, lbounds, ubounds)
+        self.regression_data = samples
+        self.regression_phi = None
+        E_vecs, shifts = self.get_data_points(samples)
+        # get design matrix
+        phi = self.basis_function_evaluation(E_vecs)
+        self.regression_phi = phi
+        self.regression_shifts = shifts
+        print('Performing initial regression fit for CB corrector.....')
+        self.RM.fit(self.regression_phi, self.regression_shifts.flatten())
+        print('Fit completed')
+
+    def check_for_refit(
+            self,
+            A,
+            atoms,
+            forces,
+            tol=1e-3,
+            mask=None,
+            F_func=None,
+            E_func=None,
+            coordinates='cart3D',
+            refit_points=10,
+            err_vec=None,
+            *args,
+            **kwargs):
+        """
+        Check and perform a refit of the Cauchy-Born shift
+        regression model based on the error metric defined above,
+        adding additional strain-states to the dataset from the
+        highest force atoms.
+
+        Parameters
+        ----------
+        A : 3x3 numpy array, floats
+            Rotation matrix of the form [x^T,y^T,z^T], where x^T,
+            y^T, z^T are normalised column vectors of the lab frame
+            directions expressed in terms of the crystal lattice
+            directions.
+        atoms : ASE atoms object
+            Atoms object where atoms sit at the coordinates needed
+            to generate the deformation gradient tensor field from
+            F_func.
+        forces : array_like
+            The forces on the atoms.
+        tol : float
+            The acceptable level of Cauchy-Born error per atom.
+        mask : array_like, bool
+            Mask for which atoms to include in the error quantification
+            process.
+        F_func : Function
+            Function to calculate the deformation gradient tensor field
+            from a set of atomic coordinates. The system of coordinates
+            that the function accepts is given by 'coordinates' and the
+            function also accepts anything extra specified in *args and
+            **kwargs. Function must return the deformation gradient field
+            for the atoms in form of a numpy array shape [natoms,ndims,ndims],
+            where F is specified in cartesian coordinates. The returned tensor
+            field must be in the lab frame. Optional, but one of F_func or
+            E_func must be provided.
+        E_func : Function
+            Function to calculate the Green-Lagrange tensor field
+            from a set of atomic coordinates. The system of
+            coordinates that the function accepts is given by
+            'coordinates' and the function also accepts anything
+            extra specified in *args and **kwargs. Function
+            must return the Green-Lagrange tensor field for
+            the atoms in form of a numpy array shape [natoms,ndims,ndims],
+            where E is specified in cartesian coordinates,
+            and the returned tensor field is in the lab frame.
+            Optional, but one of F_func or E_func must be provided.
+        coordinates : string
+            Specifies the coordinate system that the function
+            (F_func or E_func) provided accepts. The coordinates
+            of the atoms in atoms will be converted to this
+            coordinate system before being passed to the function
+            as a vector. Default is 3D cartesian coordinates.
+            Other options are cart2D, cart3D, cylind2D,
+            cylind3D, spherical
+        refit_points : int
+            The number of strain states from high force atoms to add to the
+            dataset for refitting.
+        err_vec : array_like
+            A vector of the errors on each atom which can be used in
+            place of the absolute force on each atom as a criteria for
+            selecting the 10 highest error strain states to re-fit
+            the Cauchy-Born model to.
+
+        Returns
+        -------
+        refit : integer, 1 or 0
+            0 if no refit took place, 1 if it did.
+        cb_err_per_atom : float
+            The computed Cauchy-Born shift error metric per atom
+        """
+        # tolerance is error per atom
+        # multiply by number of atoms to get full tolerance
+        if mask is not None:
+            tol = tol * len(atoms[mask])
+        else:
+            tol = tol * len(atoms)
+        print('tol', tol)
+
+        cb_err = self.get_cb_error(atoms, forces=forces, mask=mask)
+        print('cb_err', cb_err)
+        if abs(cb_err) < tol:  # do not refit in this case
+            return 0, cb_err / len(atoms)
+
+        # otherwise, find the lattice 1 and lattice 2 atoms with the highest
+        # contribution to the cauchy-born errors,
+        # add these to the training dataset and refit the model.
+        if mask is not None:
+            lattice1mask = np.logical_and(self.lattice1mask, mask)
+            lattice2mask = np.logical_and(self.lattice2mask, mask)
+        else:
+            lattice1mask = self.lattice1mask
+            lattice2mask = self.lattice2mask
+
+        if err_vec is None:
+            # get force magnitude vectors, only do this when no vector of
+            # errors is provided
+            force_mag = np.linalg.norm(forces, axis=1)
+        else:
+            force_mag = err_vec
+
+        atom_numbers = np.array([i for i in range(len(atoms))])
+
+        # work out how many points need refitting from each sublattice
+        l1points = int(refit_points / 2)
+        l2points = refit_points - l1points
+
+        # get the indices of the highest force atoms from each sublattice and
+        # select these for refitting.
+        force_mag_lattice1 = force_mag[lattice1mask]
+        atom_nos_lattice1 = atom_numbers[lattice1mask]
+        top_idx = np.argpartition(force_mag_lattice1, -l1points)[-l1points:]
+        high_force_atoms1 = atom_nos_lattice1[top_idx]
+
+        # get 5 largest lattice 2 force indices
+        force_mag_lattice2 = force_mag[lattice2mask]
+        atom_nos_lattice2 = atom_numbers[lattice2mask]
+        top_idx = np.argpartition(force_mag_lattice2, -l2points)[-l2points:]
+        high_force_atoms2 = atom_nos_lattice2[top_idx]
+
+        high_force_atoms = np.concatenate(
+            (high_force_atoms1, high_force_atoms2))
+        # print('atoms selected:', high_force_atoms)
+        # print('forces on atoms:', err_vec[high_force_atoms])
+        reduced_atoms = atoms[high_force_atoms]
+
+        sx, sy, sz = atoms.cell.diagonal()
+        cell = [sx, sy, sz]
+
+        # evaluate E and R for reduced list of atoms
+        E, R = self.evaluate_F_or_E(A, reduced_atoms, F_func=F_func,
+                                    E_func=E_func, cell=cell,
+                                    coordinates=coordinates, *args, **kwargs)
+        # print('Es',E)
+        # we get back the strain vectors in the lattice frame
+        # turn E into voigt vectors
+        E_voigt = np.zeros([np.shape(E)[0], 6])
+        E_voigt[:, 0] = E[:, 0, 0]
+        E_voigt[:, 1] = E[:, 1, 1]
+        E_voigt[:, 2] = E[:, 2, 2]
+        E_voigt[:, 3] = E[:, 1, 2]
+        E_voigt[:, 4] = E[:, 0, 2]
+        E_voigt[:, 5] = E[:, 0, 1]
+
+        # pass set of Es to evaluate to refit_regression
+        self.refit_regression(atoms, E_voigt)
+        return 1, cb_err / len(atoms)
+
+    def refit_regression(self, atoms, E_voigt):
+        """Refit the regression model to a new set of data.
+
+        Parameters
+        ----------
+        atoms : ASE atoms object
+            Atoms object where atoms sit at the coordinates
+            needed to generate the deformation gradient tensor field
+            from F_func.
+        E_voigt : array_like
+            Numpy array of the Green-Lagrange strain tensors
+            (in Voigt notation) to refit to.
+        """
+        calc = self.calc
+        # take our given E voigt and evaluate the model
+        E_vecs, shifts = self.get_data_points(E_voigt)
+        # print('evecs,shifts',E_vecs,shifts)
+        # get design matrix
+        phi = self.basis_function_evaluation(E_vecs)
+        temp_predict = self.RM.predict(phi)
+        # print('predicted shifts', temp_predict)
+        if self.regression_phi is None:
+            self.regression_phi = phi
+            self.regression_shifts = shifts
+        else:
+            self.regression_phi = np.concatenate(
+                (self.regression_phi, phi), axis=0)
+            self.regression_shifts = np.concatenate(
+                (self.regression_shifts, shifts), axis=0)
+
+        print('High error detected, refitting CB corrector.....')
+        # refit the regression model to the given data
+        self.RM.fit(self.regression_phi, self.regression_shifts.flatten())
+        print('Re-fit completed')
+
+    def get_data_points(self, E_vec):
+        """Get the 'true' shifts (by building and relaxing a unit cell under
+        that strain state) using a series of supplied strain vectors. The
+        cauchy-born shift correction along all 3 axes can be measured from one
+        applied strain state and relaxation. The regression model only predicts
+        the component of the shift vector parallel to axis 1, but we can still
+        use all 3 predicted values for fitting if we find the rotated strain
+        tensors about 111 which would predict the shifts along axis 2 and 3.
+
+        Parameters
+        ----------
+        E_vec : array_like
+            Numpy array of the Green-Lagrange strain
+            tensors (in Voigt notation) to apply
+
+        Returns
+        -------
+        strains : array_like
+            A larger array including the supplied E_vec array,
+            as well as the versions of E_vec which have been
+            rotated about 120 degrees about the [111] axis.
+        shift_vals : array_like
+            An array of shift components along the 1 axis
+            corresponding to the matching strain state in 'strains'.
+        """
+        unitcell = self.lattice(
+            self.el, latticeconstant=self.a0, size=[1, 1, 1])
+        unitcell.set_pbc([True, True, True])
+
+        def permutation(strain_vec, perm_shift):
+            """Rotate a Green-Lagrange strain tensor in Voigt notation about
+            the 111 axis a number of times given by perm_shift. This
+            corresponds to permuting the top three components and bottom three
+            components seperately. For example, the vector [1,2,3,4,5,6]
+            becomes [2,3,1,5,6,4] under this operation. Due to the triad 111
+            symmetry in diamond, feeding these rotated strain states to the
+            shift prediction model allows the prediction of the full shift
+            vector.
+
+            Parameters
+            ----------
+            strain_vec : array_like
+                2D array of size [natoms,6], where each row is a Green-Lagrange
+                strain tensor written in Voigt notation.
+            perm_shift : integer
+                The permeutation number wanted, with 1 corresponding
+                to a 120 degree rotation about 111 (these strain vectors give
+                the y component of shift prediction when fed to the model)
+                and 2 corresponding to a 240 degree rotation
+                about 111(these strain vectors give the z component
+                of shift prediction when fed to the model)
+
+            Returns
+            -------
+            strain_vec_perm : array_like
+                Rotated version of strain_vec
+            """
+            strain_vec_perm = np.zeros_like(strain_vec)
+            for i in range(3):
+                perm_num = (i - perm_shift) % 3
+                strain_vec_perm[:, perm_num] = strain_vec[:, i]
+                strain_vec_perm[:, perm_num + 3] = strain_vec[:, i + 3]
+            return strain_vec_perm
+
+        ndp = np.shape(E_vec)[0]  # number of data points
+        shift_vals = np.zeros([3 * ndp, 1])
+        for i, E in enumerate(E_vec):
+            # get the Cauchy-Born shift correction
+            shift_diff = self.eval_shift(E, unitcell)
+            shift_vals[i] = shift_diff[0]
+            shift_vals[i + ndp] = shift_diff[1]
+            shift_vals[i + 2 * ndp] = shift_diff[2]
+
+        # Permute the supplied strain vectors to rotate about 111
+        eps = E_vec
+        epsx = eps
+        epsy = permutation(eps, 1)
+        epsz = permutation(eps, 2)
+
+        strains = np.zeros([3 * ndp, 6])
+        strains[0:ndp, :] = epsx
+        strains[ndp:2 * ndp, :] = epsy
+        strains[2 * ndp:3 * ndp, :] = epsz
+
+        return strains, shift_vals
+
+    def eval_shift(self, E_vec, unitcell):
+        """Function that calculates the cauchy-born shift corrector for a given
+        Green-Lagrange strain tensor.
+
+        Parameters
+        ----------
+        E_vec : array_like
+            Green-Lagrange strain state to apply in voigt notation
+        unitcell : ASE atoms object
+            Single unit cell of crystal to apply strain vector to and to find
+            shift
+
+        Returns
+        -------
+        shift_diff : array_like
+            Vector of cauchy-born shift correction for applied strain state
+        """
+        # green lagrange version of function
+        # turn E into matrix
+        E = np.zeros([3, 3])
+        E[0, 0] = E_vec[0]
+        E[1, 1] = E_vec[1]
+        E[2, 2] = E_vec[2]
+        E[2, 1], E[1, 2] = E_vec[3], E_vec[3]
+        E[2, 0], E[0, 2] = E_vec[4], E_vec[4]
+        E[0, 1], E[1, 0] = E_vec[5], E_vec[5]
+        # get U^2
+        Usqr = 2 * E + np.eye(3)
+        # square root matrix
+        U = sqrtm(Usqr, disp=True)
+
+        # this is just the symmetric stretch tensor, exactly what we need.
+        x = U
+        initial_shift = np.zeros([3])
+        relaxed_shift = np.zeros([3])
+
+        # build cell
+        unitcell_copy = unitcell.copy()
+        unitcell_copy.calc = self.calc
+        cell = unitcell_copy.get_cell()
+        cell_rescale = np.transpose(x[:, :] @ cell)
+        unitcell_copy.set_cell(cell_rescale, scale_atoms=True)
+        initial_shift[:] = unitcell_copy.get_positions()[1] - \
+            unitcell_copy.get_positions()[0]
+        # relax cell
+        opt = LBFGS(unitcell_copy, logfile=None)
+        opt.run(fmax=1e-10)
+        relaxed_shift[:] = unitcell_copy.get_positions()[1] - \
+            unitcell_copy.get_positions()[0]
+
+        shift_diff = relaxed_shift - initial_shift
+        # TRANSFER SHIFT_DIFF BACK INTO THE LATTICE FRAME, OUT OF THE DEFORMED
+        # FRAME
+        # back_transform = np.transpose(cell_rescale)
+        # shift_diff_transform = back_transform@shift_diff
+        # print('shift_diff',shift_diff,'shift_diff_transform',shift_diff_transform)
+        return shift_diff  # return all 3 components
+
+    def basis_function_evaluation(self, E_vecs):
+        """Function that evaluates the polynomial basis functions used in the
+        regression model on a set of strain vectors to build the design matrix
+        phi.
+
+        Parameters
+        ----------
+        E_vecs : array_like
+            Set of Green-Lagrange strain states in voigt notation to fit to
+
+        Returns
+        -------
+        phi : array_like
+            Design matrix composed of all of the basis functions evaluated at
+            every data point given in E_vecs.
+        """
+        # first basis function is just the e vector components evaluated, 2nd
+        # is all of the cross terms
+        triu_indices = np.triu_indices(6)  # get upper triangle indices
+        phi_2nd_term = np.zeros(
+            [np.shape(E_vecs)[0], np.shape(triu_indices)[1]])
+        for i, E in enumerate(E_vecs):
+            e_outer_p = np.outer(E, E)
+            phi_2nd_term[i, :] = e_outer_p[triu_indices]
+
+        # return the design matrix
+
+        return np.concatenate((E_vecs, phi_2nd_term), axis=1)
diff --git a/matscipy/dislocation.py b/matscipy/dislocation.py
index 48982754..fc23479b 100644
--- a/matscipy/dislocation.py
+++ b/matscipy/dislocation.py
@@ -32,6 +32,7 @@
                                Diamond, SimpleCubicFactory)
 from ase.constraints import FixAtoms, StrainFilter
 from ase.optimize import FIRE
+from ase.optimize.precon import PreconLBFGS
 from ase.build import bulk
 from ase.calculators.lammpslib import LAMMPSlib
 from ase.units import GPa  # unit conversion
diff --git a/matscipy/fracture_mechanics/clusters.py b/matscipy/fracture_mechanics/clusters.py
index 36d8bc85..1b8d5c78 100644
--- a/matscipy/fracture_mechanics/clusters.py
+++ b/matscipy/fracture_mechanics/clusters.py
@@ -23,6 +23,7 @@
 import numpy as np
 
 from ase.lattice.cubic import Diamond, FaceCenteredCubic, SimpleCubic, BodyCenteredCubic
+from matscipy.neighbours import neighbour_list
 
 ###
 
@@ -58,18 +59,57 @@ 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):
+def set_regions(cryst, r_I, cutoff, r_III, extended_far_field=False,extended_region_I=False,exclude_surface=False):
     sx, sy, sz = cryst.cell.diagonal()
     x, y = cryst.positions[:, 0], cryst.positions[:, 1]
     cx, cy = sx/2, sy/2
     r = np.sqrt((x - cx)**2 + (y - cy)**2)
 
-    # Regions I-III defined by radial distance from center
+    # Regions I and III defined by radial distance from center
     regionI = r < r_I
     regionII = (r >= r_I) & (r < (r_I + cutoff))
-    regionIII = (r >= (r_I + cutoff)) & (r < r_III)
+    regionIII = (r >= r_I+cutoff) & (r < r_III)
+    #regionIII = (r >= r_I) & (r < r_III)
     regionIV = (r >= r_III) & (r < (r_III + cutoff))
 
+    """    regionII = np.zeros(len(cryst), bool)
+    i, j = neighbour_list('ij', cryst, cutoff)
+    for idx in regionI.nonzero()[0]:
+        neighbs = j[i == idx] 
+        mask = np.zeros(len(cryst), bool)
+        mask[neighbs] = True # include all neighbours of atom `idx`
+        # print(f'adding {mask.sum()} neigbours of atom {idx} to regionII')
+        mask[regionI] = False # exclude those in region I already
+        regionII[mask] = True # add to region I
+        regionIII[mask] = False # remove from region III"""
+
+    if exclude_surface or extended_region_I:
+        #build a mask of the material surface based on the following criteria:
+        # - the atoms x coordinate should lie between rI and 
+        #  rI + 2*cutoff (if non extended) and rIII + cutoff if extended.
+        # - the atoms y coordinate should lie between +- cutoff
+        if extended_far_field:
+            x_criteria = ((x-cx)<-(r_I-cutoff))&\
+                    (((x-cx)>-((r_III)+(cutoff))))
+        else:
+            x_criteria = ((x-cx)<-((r_I-cutoff)))&\
+                (((x-cx)>-((r_I)+((2*cutoff)))))
+        y_criteria = ((y-cy)>-(cutoff))&\
+                ((y-cy)<(cutoff))
+        surface_mask = np.logical_and(x_criteria,y_criteria)
+    if extended_region_I:
+        #re-do all the region logic so they are correct
+        #(and followed by an xor to exclude new region I atoms.)
+        regionI = regionI|surface_mask
+        regionII = regionII^np.logical_and(regionII,surface_mask)
+        regionIII = regionIII^np.logical_and(regionIII,surface_mask)
+        regionIV = regionIV^np.logical_and(regionIV,surface_mask)
+    
+    elif exclude_surface:
+        regionII = regionII^np.logical_and(regionII,surface_mask)
+        regionIII = regionIII^np.logical_and(regionIII,surface_mask)
+        regionIV = regionIV|(surface_mask^np.logical_and(surface_mask,regionI)) #add the surface mask to region IV
+    
     cryst.new_array('region', np.zeros(len(cryst), dtype=int))
     region = cryst.arrays['region']
     region[regionI]  = 1
@@ -130,3 +170,4 @@ def bcc(*args, **kwargs):
 def sc(*args, **kwargs):
     kwargs['lattice'] = SimpleCubic
     return cluster(*args, **kwargs)
+
diff --git a/matscipy/fracture_mechanics/crack.py b/matscipy/fracture_mechanics/crack.py
index fb3564f9..6a80aef4 100644
--- a/matscipy/fracture_mechanics/crack.py
+++ b/matscipy/fracture_mechanics/crack.py
@@ -36,14 +36,22 @@
 
 import ase.units as units
 from ase.optimize.precon import Exp
+#from ase.optimize.ode import ode12r
 
 from matscipy.atomic_strain import atomic_strain
 from matscipy.elasticity import (rotate_elastic_constants,
                                  rotate_cubic_elastic_constants,
                                  Voigt_6_to_full_3x3_stress)
 from matscipy.surface import MillerDirection, MillerPlane
-
 from matscipy.neighbours import neighbour_list
+from matscipy.optimize import ode12r
+
+import matplotlib.pyplot as plt
+from ase.optimize.precon import Exp, PreconLBFGS
+from ase.optimize import LBFGS
+from ase.constraints import FixAtoms
+from ase import Atom
+import ase.io
 
 ###
 
@@ -53,7 +61,13 @@
 
 MPa_sqrt_m = 1e6*units.Pascal*np.sqrt(units.m)
 
-###
+#wrapper to count function calls for writing
+def counted(f):
+    def wrapped(*args):
+        wrapped.calls += 1
+        return f(*args,str(wrapped.calls))
+    wrapped.calls = 0
+    return wrapped
 
 class RectilinearAnisotropicCrack:
     """
@@ -359,7 +373,8 @@ class CubicCrystalCrack:
     """
 
     def __init__(self, crack_surface, crack_front, C11=None, C12=None,
-                 C44=None, stress_state=PLANE_STRAIN, C=None, Crot=None):
+                 C44=None, stress_state=PLANE_STRAIN, C=None, Crot=None,
+                 cauchy_born=None):
         """
         Initialize a crack in a cubic crystal with elastic constants C11, C12
         and C44 (or optionally a full 6x6 elastic constant matrix C).
@@ -405,6 +420,12 @@ def __init__(self, crack_surface, crack_front, C11=None, C12=None,
                                         S6[0, 5], S6[1, 5], S6[2, 5],
                                         S6[5, 5])
 
+        self.third_dir = third_dir
+        self.crack_surface = crack_surface
+        self.crack_front = crack_front
+        self.RotationMatrix = A
+        self.cauchy_born = cauchy_born
+
 
     def k1g(self, surface_energy):
         """
@@ -715,7 +736,7 @@ class SinclairCrack:
     def __init__(self, crk, cryst, calc, k, alpha=0.0, vacuum=6.0,
                  variable_alpha=True, variable_k=False,
                  alpha_scale=None, k_scale=None,
-                 extended_far_field=False):
+                 extended_far_field=False, rI=0.0, rIII=0.0, cutoff=0.0):
         """
 
         Parameters
@@ -755,7 +776,10 @@ def __init__(self, crk, cryst, calc, k, alpha=0.0, vacuum=6.0,
         self.u = np.zeros((self.N1, 3))
         self.alpha = alpha
         self.k = k
-
+        self.rI = rI
+        self.rIII = rIII
+        self.cutoff = cutoff
+        self.shiftmask = None
         # check atoms are sorted by distance from centre so we can use N1,N2,N3
         tip_x = cryst.cell.diagonal()[0] / 2.0
         tip_y = cryst.cell.diagonal()[1] / 2.0
@@ -770,13 +794,15 @@ 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]
 
         self.precon = None
         self.precon_count = 0
+        self.f_alpha_correction = 0
+        self.alpha0 = alpha
 
     def pack(self, u, alpha, k):
         dofs = list(u.reshape(-1))
@@ -810,6 +836,10 @@ def get_dofs(self):
 
     def set_dofs(self, x):
         self.u[:], self.alpha, self.k = self.unpack(x, reshape=True)
+        print('alpha',self.alpha)
+        if self.variable_k:
+            print('k',self.k)
+        #self.k = self.k-(np.abs(self.alpha-self.alpha0)*self.k1g)
         self.update_atoms()
 
     def __len__(self):
@@ -836,6 +866,7 @@ def u_cle(self, alpha=None):
         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):
         def residuals(x, mask):
@@ -876,6 +907,12 @@ def residuals(x, mask):
                 raise RuntimeError('CLE fit failed')
         return res
 
+    def get_deformation_gradient(self,r,theta,k):
+        return self.crk.crack.deformation_gradient(r, theta, k)
+
+    def set_shiftmask(self,radial_dist):
+        self.shiftmask = self.r>radial_dist
+    
     def update_atoms(self):
         """
         Update self.atoms from degrees of freedom (self.u, self.alpha, self.k)
@@ -890,12 +927,57 @@ def update_atoms(self):
         self.atoms.positions[:, :] = self.cryst.positions
         self.atoms.positions[:, :] += self.k * 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 
+            A = np.transpose(self.crk.RotationMatrix)
+            #find shifts
+            print('finding shifts.....')
+            # 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='cylind2D',method='regression',k=self.k)
+            print('done!')
+            #apply shifts
+            self.crk.cauchy_born.apply_shifts(self.atoms,shifts,mask=self.shiftmask)
+
 
         # add vacuum
         self.atoms.cell = self.cryst.cell
         self.atoms.cell[0, 0] += self.vacuum
         self.atoms.cell[1, 1] += self.vacuum
 
+    def get_f_alpha_correction(self):
+        #function which gets the f_alpha contribution from regions II and III without u being applied
+        #'update atoms' but don't add u
+        self.atoms.set_pbc([False, False, True])
+        self.atoms.calc = self.calc
+
+        self.atoms.info['k'] = self.k
+        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()
+        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)
+            #find shifts
+            print('finding shifts.....')
+            # 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='cylind2D',method='regression',k=self.k)
+            print('done!')
+            #apply shifts
+            self.crk.cauchy_born.apply_shifts(self.atoms,shifts,mask=self.shiftmask)
+        if self.extended_far_field:
+            mask = self.regionII|self.regionIII
+        else:
+            mask = self.regionII
+        f_alpha = self.get_crack_tip_force(mask=mask)
+        self.update_atoms()
+        return f_alpha
+
     def set_atoms(self, atoms):
         N1_in = (atoms.arrays['region'] == 1).sum()
         if 'alpha' in atoms.info:
@@ -907,13 +989,13 @@ def set_atoms(self, atoms):
         self.update_atoms()  # now we have same u_cle in atoms and self.atoms
         min_len = min(N1_in, self.N1)
         # FIXME this assumes stable sort order for atoms and self.atoms
-        u = atoms.positions[:min_len] - self.atoms.positions[:min_len]
+        u = atoms.positions - self.atoms.positions
         shift = np.diag(self.atoms.cell)/2 - np.diag(atoms.cell)/2
         u += shift
-        self.u[:min_len] = u
+        self.u = u[self.regionI]
         self.update_atoms()
 
-    def get_crack_tip_force(self, forces=None, mask=None):
+    def get_crack_tip_force(self, forces=None, mask=None, full_array_output=False):
         # V_alpha = -\nabla_1 U_CLE(alpha)
         tip_x = self.cryst.cell.diagonal()[0] / 2.0 + self.alpha
         tip_y = self.cryst.cell.diagonal()[1] / 2.0
@@ -947,6 +1029,10 @@ def get_crack_tip_force(self, forces=None, mask=None):
             mask = self.regionII
             if self.extended_far_field:
                 mask = self.regionII | self.regionIII
+        if full_array_output is True:
+            reduced_forces = forces[mask,:]
+            reduced_V = V[mask,:]
+            return np.tensordot(forces[mask, :], V[mask, :]), np.array([np.dot(reduced_forces[i, :], reduced_V[i, :]) for i in range(np.shape(reduced_forces)[0])])
         return np.tensordot(forces[mask, :], V[mask, :])
 
     def get_xdot(self, x1, x2, ds=None):
@@ -989,10 +1075,13 @@ def get_forces(self, x1=None, xdot1=None, ds=None, forces=None, mask=None):
         F = list(forces[self.regionI, :].reshape(-1))
         if self.variable_alpha:
             f_alpha = self.get_crack_tip_force(forces, mask=mask)
-            F.append(f_alpha)
+            self.f_alpha_correction = 0.0 #self.get_f_alpha_correction()
+            F.append((f_alpha-self.f_alpha_correction))
+            print('f_alpha',f_alpha)
         if self.variable_k:
             f_k = self.get_k_force(x1, xdot1, ds)
             F.append(f_k)
+
         return np.array(F)
 
     def update_precon(self, x, F=None):
@@ -1057,7 +1146,7 @@ def update_precon(self, x, F=None):
             Pf = self.P_ilu.solve(F)
             print(f'norm(F) = {np.linalg.norm(F)}, norm(P^-1 F) = {np.linalg.norm(Pf)}')
             Pfu, Pfalpha, Pfk = self.unpack(Pf)
-            print(f'|P^-1 f_I| = {np.linalg.norm(Pfu, np.inf)}, P^-1 f_alpha = {Pfalpha}')
+            print(f'|P^-1 f_I| = {np.linalg.norm(Pfu, np.inf)}, P^-1 f_alpha = {Pfalpha}, P^-1 f_k = {Pfk}')
 
     def get_precon(self, x, F):
         self.update_precon(x, F)
@@ -1192,19 +1281,54 @@ def cg2_jacobian(x):
                                 'maxiter': steps,
                                 'jac_options': {'inner_M': M}},
                        callback=log)
+            
+        elif method == 'ode12r':
+
+            class ODEResult:
+                def __init__(self, success, x, nit):
+                    self.success = success
+                    self.x = x
+                    self.nit = nit
+
+            if args is None:
+                # args==None -> variable_k=False (flex1 and flex2)
+                x, nit = ode12r(residuals, x0, verbose=2, fmax=ftol, steps=steps)
+
+            else:
+                # args!=None -> variable_k=True (arc_continuation)
+
+                def oderes(f,x):
+                    return np.linalg.norm(f, np.inf)
+
+                if precon:
+                    def odeprecon(f,x):
+                        M = self.get_precon(x,f)
+                        Fp = M@f
+                        Rp = oderes(f,x)
+                        return Fp, Rp # returns preconditioned forces and residual
+                else:
+                    def odeprecon(f,x):
+                        Rp = oderes(f,x)
+                        return f, Rp #this precon function doesn't do anything
+                    
+                x, nit = ode12r(residuals, x0, args, verbose=1, fmax=ftol, steps=steps, apply_precon=odeprecon, residual=oderes, maxtol=np.inf) 
+            
+            res = ODEResult(True, x, nit)
+
         else:
             raise RuntimeError(f'unknown method {method}')
 
         if res.success:
             self.set_dofs(res.x)
+            return res.nit
         else:
             self.atoms.write('no_convergence.xyz')
             raise RuntimeError(f"no convergence of scipy optimizer {method}")
 
     def get_potential_energy(self):
         # E1: energy of region I and II atoms
-        E = self.atoms.get_potential_energies()[self.regionI_II].sum()
-        E1 = E - self.E0
+        #E = self.atoms.get_potential_energies()[self.regionI_II].sum()
+        #E1 = E - self.E0
 
         # E2: energy of far-field (region III)
         regionII_III = self.regionII | self.regionIII
@@ -1237,10 +1361,14 @@ def rescale_k(self, new_k):
                           z - ref_z][self.regionI, :]
 
     def arc_length_continuation(self, x0, x1, N=10, ds=0.01, ftol=1e-2,
-                                direction=1, steps=100,
+                                direction=1, max_steps=10,
                                 continuation=False, traj_file='x_traj.h5',
                                 traj_interval=1,
-                                precon=False):
+                                precon=False,
+                                ds_max=np.inf, ds_min=0,
+                                ds_aggressiveness=2,
+                                opt_method='krylov',
+                                cos_alpha_min=0.9):
         import h5py
         assert self.variable_k  # only makes sense if K can vary
 
@@ -1254,7 +1382,7 @@ def arc_length_continuation(self, x0, x1, N=10, ds=0.01, ftol=1e-2,
             _, alphadot1, _ = self.unpack(xdot1)
             if direction * np.sign(alphadot1) < 0:
                 xdot1 = -xdot1
-
+        
         row = 0
         with h5py.File(traj_file, 'a') as hf:
             if 'x' in hf.keys():
@@ -1269,21 +1397,47 @@ def arc_length_continuation(self, x0, x1, N=10, ds=0.01, ftol=1e-2,
                 x_traj.attrs['traj_interval'] = traj_interval
             row = x_traj.shape[0]
 
-        for i in range(N):
+        i = 0
+        while i < N:
             x2 = x1 + ds * xdot1
             print(f'ARC LENGTH step={i} ds={ds}, k1 = {x1[-1]:.3f}, k2 = {x2[-1]:.3f}, '
                   f' |F| = {np.linalg.norm(self.get_forces(x1=x1, xdot1=xdot1, ds=ds)):.4f}')
             self.set_dofs(x2)
-            self.optimize(ftol, steps, args=(x1, xdot1, ds), precon=precon)
+            try:
+                if opt_method == 'krylov':
+                    num_steps = self.optimize(ftol, max_steps, args=(x1, xdot1, ds), precon=precon, method=opt_method,verbose=2)
+                elif opt_method == 'ode12r':
+                    num_steps = self.optimize(ftol, max_steps, args=[x1, xdot1, ds], precon=precon, method=opt_method,verbose=2)
+                print(f'Corrector converged in {num_steps}/{max_steps} steps')
+            except NoConvergence:
+                if ds < ds_min:
+                    print(f'Corrector failed to converge even with ds={ds}<{ds_min}. Aborting job.')
+                    break
+                else:
+                    ds /= ds_aggressiveness
+                    print(f'Corrector failed to converge, reducing step size to ds={ds}')
+                    continue
+
             x2 = self.get_dofs()
             xdot2 = self.get_xdot(x1, x2, ds)
 
-            # monitor sign of \dot{alpha} and flip if necessary
             if self.variable_alpha:
-                _, alphadot2, _ = self.unpack(xdot2)
+                # monitor sign of \dot{alpha} and flip if necessary
+                #udot1, alphadot1, kdot1 = self.unpack(xdot1)
+                udot2, alphadot2, kdot2 = self.unpack(xdot2)
                 if direction * np.sign(alphadot2) < 0:
-                    xdot2 = -xdot2
-
+                    xdot2 = -xdot2  
+
+            # cos_alpha = np.dot(xdot1, xdot2) / np.linalg.norm(xdot1) / np.linalg.norm(xdot2)
+            # print(f'cos_alpha = {cos_alpha}')
+            # if cos_alpha < cos_alpha_min:
+            #     print("Angle between subsequent predictors too large "
+            #           f"(cos(alpha) = {cos_alpha} < {cos_alpha_min}). "
+            #            "Restart with smaller step size.")
+            #     ds *= 0.5
+            #     continue
+
+            # Write new relaxed point to the traj file
             if i % traj_interval == 0:
                 for nattempt in range(1000):
                     try:
@@ -1299,12 +1453,20 @@ def arc_length_continuation(self, x0, x1, N=10, ds=0.01, ftol=1e-2,
                         continue
                 else:
                     raise IOError("ran out of attempts to access trajectory file")
-
+            
+            # Update variables
             x1[:] = x2
             xdot1[:] = xdot2
 
-    def plot(self, ax=None, regions='1234', styles=None, bonds=None, cutoff=2.8,
-             tip=False, atoms_args=None, bonds_args=None, tip_args=None):
+            # Update the stepsize
+            ds *= ( 1 + ds_aggressiveness * ((max_steps - num_steps) / (max_steps - 1)) ** 2 )
+            ds = min(ds_max, ds) 
+
+            # Increase iteration number
+            i += 1
+
+    def plot(self, ax=None, regions='1234', rzoom=np.inf, bonds=None, cutoff=2.8, tip=False, 
+             regions_styles=None, atoms_args=None, bonds_args=None, tip_args=None):
         import matplotlib.pyplot as plt
         from matplotlib.collections import LineCollection
 
@@ -1317,17 +1479,20 @@ def plot(self, ax=None, regions='1234', styles=None, bonds=None, cutoff=2.8,
 
         if ax is None:
             fig, ax = plt.subplots()
+
         if isinstance(regions, str):
             regions = [int(r) for r in regions]
         a = self.atoms
         region = a.arrays['region']
-        if styles is None:
-            styles = ['bo', 'ko', 'r.', 'rx']
         plot_elements = []
-        for i, fmt in zip(regions, styles):
-            (p,) = ax.plot(a.positions[region == i, 0],
-                           a.positions[region == i, 1], **atoms_args)
+        if regions_styles is None:
+            regions_styles = ['bo', 'ko', 'r.', 'rx'][:len(regions)] #upto 4 regions
+        for i, fmt in zip(regions, regions_styles):
+            iatoms = np.array([ i for i,rad in enumerate(self.r[region==i]) if rad<rzoom ]) #zoom
+            (p,) = ax.plot(a.positions[region == i, 0][iatoms],
+                           a.positions[region == i, 1][iatoms], **atoms_args)
             plot_elements.append(p)
+
         if bonds:
             if isinstance(bonds, bool):
                 bonds = regions
@@ -1340,6 +1505,7 @@ def plot(self, ax=None, regions='1234', styles=None, bonds=None, cutoff=2.8,
             lc = LineCollection(lines, **bonds_args)
             ax.add_collection(lc)
             plot_elements.append(lc)
+
         if tip:
             (tip,) = ax.plot(self.cryst.cell[0, 0] / 2.0 + self.alpha,
                              self.cryst.cell[1, 1] / 2.0, **tip_args)
@@ -1348,10 +1514,11 @@ def plot(self, ax=None, regions='1234', styles=None, bonds=None, cutoff=2.8,
         ax.set_aspect('equal')
         ax.set_xticks([])
         ax.set_yticks([])
+
         return plot_elements
 
-    def animate(self, x, k1g, regions='12', cutoff=2.8, frames=None,
-                callback=None):
+    def animate(self, x, k1g, regions='12', rzoom=np.inf, cutoff=2.8, frames=None, callback=None,
+                plot_tip=True, regions_styles=None, atoms_args=None, bonds_args=None, tip_args=None):
         import matplotlib.pyplot as plt
         from matplotlib.animation import FuncAnimation
         fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 12))
@@ -1370,9 +1537,17 @@ def animate(self, x, k1g, regions='12', cutoff=2.8, frames=None,
         ax1.set_ylabel(r'Stress intensity factor $K/K_{G}$');
 
         self.set_dofs(x[i, :])
-        plot_elements = self.plot(ax2, regions=regions, bonds=regions, tip=True)
-        tip = plot_elements.pop(-1)
-        lc = plot_elements.pop(-1)
+        if cutoff is None:
+            # Do not plot bonds (eg. for metals)
+            plot_elements = self.plot(ax2, regions=regions, rplot=rzoom, bonds=False, cutoff=cutoff, tip=plot_tip, 
+                                      regions_styles=regions_styles, atoms_args=atoms_args, bonds_args=bonds_args, tip_args=tip_args)
+            tip = plot_elements.pop(-1)
+            lc = None # placeholder, no bond lines
+        else:
+            plot_elements = self.plot(ax2, regions=regions, rzoom=rzoom, bonds=regions, cutoff=cutoff, tip=plot_tip, 
+                                      regions_styles=regions_styles, atoms_arg=atoms_args, bonds_args=bonds_args, tip_args=tip_args)
+            tip = plot_elements.pop(-1)
+            lc = plot_elements.pop(-1)
 
         if callback:
             callback(ax1, ax2)
@@ -1387,11 +1562,12 @@ def frame(idx):
                 p.set_data(
                     [a.positions[region == r, 0], a.positions[region == r, 1]])
             # update bonds - requires a neighbour list recalculation
-            i, j = neighbour_list('ij', a, cutoff)
-            i, j = np.array([(I, J) for I, J in zip(i, j) if
-                             region[I] in regions and region[J] in regions]).T
-            lines = list(zip(a.positions[i, 0:2], a.positions[j, 0:2]))
-            lc.set_segments(lines)
+            if cutoff is not None:
+                i, j = neighbour_list('ij', a, cutoff)
+                i, j = np.array([(I, J) for I, J in zip(i, j) if
+                                region[I] in regions and region[J] in regions]).T
+                lines = list(zip(a.positions[i, 0:2], a.positions[j, 0:2]))
+                lc.set_segments(lines)
             # update crack tip indicator
             tip.set_data([self.cryst.cell[0, 0] / 2.0 + x[idx, -2]],
                          [self.cryst.cell[1, 1] / 2.0])
@@ -1399,8 +1575,204 @@ def frame(idx):
 
         if frames is None:
             frames = range(0, len(x), 100)
+
         return FuncAnimation(fig, frame, frames)
 
+    @counted
+    def write_atoms_to_file(self, fname):
+        alpha_at = Atom('Au')
+        crack_atoms = self.atoms.copy()
+
+        r = crack_atoms.arrays['region']
+        r[self.regionI] = 1
+        r[self.regionII] = 2
+        r[self.regionIII] = 3
+        r[self.regionIV] = 4
+
+        crack_atoms.calc = self.calc
+        forces = crack_atoms.get_forces()
+        crack_atoms.append(alpha_at)
+        crack_atoms[-1].position = [self.cryst.cell[0, 0] / 2.0 + self.alpha,
+                             self.cryst.cell[1, 1] / 2.0,0]
+        crack_atoms.new_array('fx', np.zeros(len(crack_atoms), dtype=float))
+        crack_atoms.new_array('fy', np.zeros(len(crack_atoms), dtype=float))
+        crack_atoms.new_array('fz', np.zeros(len(crack_atoms), dtype=float))
+        crack_atoms.new_array('ftot', np.zeros(len(crack_atoms), dtype=float))
+        crack_atoms.new_array('logabsftot', np.zeros(len(crack_atoms), dtype=float))
+        crack_atoms.new_array('falphacomponents',np.zeros(len(crack_atoms)),dtype =float)
+        crack_atoms.new_array('logabsfalphacomp',np.zeros(len(crack_atoms)),dtype =float)
+        crack_atoms.new_array('exx', np.zeros(len(crack_atoms),dtype=float))
+        crack_atoms.new_array('eyy',np.zeros(len(crack_atoms)),dtype =float)
+        crack_atoms.new_array('ezz',np.zeros(len(crack_atoms)),dtype =float)
+        crack_atoms.new_array('eyz',np.zeros(len(crack_atoms)),dtype =float)
+        crack_atoms.new_array('exz',np.zeros(len(crack_atoms)),dtype =float)
+        crack_atoms.new_array('exy',np.zeros(len(crack_atoms)),dtype =float)
+
+        crack_atoms.new_array('ux', np.zeros(len(crack_atoms), dtype=float))
+        crack_atoms.new_array('uy', np.zeros(len(crack_atoms), dtype=float))
+        crack_atoms.new_array('uz', np.zeros(len(crack_atoms), dtype=float))
+        
+
+        fx = crack_atoms.arrays['fx']
+        fy = crack_atoms.arrays['fy']
+        fz = crack_atoms.arrays['fz']
+        ux = crack_atoms.arrays['ux']
+        uy = crack_atoms.arrays['uy']
+        uz = crack_atoms.arrays['uz']
+        ftot = crack_atoms.arrays['ftot']
+        logabsftot = crack_atoms.arrays['logabsftot']
+        falphacomponents = crack_atoms.arrays['falphacomponents']
+        logabsfalphacomp = crack_atoms.arrays['logabsfalphacomp']
+        store_strain=False
+        if store_strain:
+            exx = crack_atoms.arrays['exx']
+            eyy = crack_atoms.arrays['eyy']
+            ezz = crack_atoms.arrays['ezz']
+            eyz = crack_atoms.arrays['eyz']
+            exz = crack_atoms.arrays['exz']
+            exy = crack_atoms.arrays['exy']
+
+            #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='cylind2D',k=self.k)
+            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
+            exx[0:len(crack_atoms)-1] = E[:,0,0]
+            eyy[0:len(crack_atoms)-1] = E[:,1,1]
+            ezz[0:len(crack_atoms)-1] = E[:,2,2]
+            eyz[0:len(crack_atoms)-1] = E[:,1,2]
+            exz[0:len(crack_atoms)-1] = E[:,0,2]
+            exy[0:len(crack_atoms)-1] = E[:,0,1]
+        #store other things
+        fx[0:len(crack_atoms)-1] = (forces[:,0])
+        fy[0:len(crack_atoms)-1] = (forces[:,1])
+        fz[0:len(crack_atoms)-1] = (forces[:,2])
+        ftot[0:len(crack_atoms)-1] = (np.linalg.norm(forces,ord=2,axis=1))
+        ux[0:len(crack_atoms)-1][self.regionI] = (self.u[:,0])
+        uy[0:len(crack_atoms)-1][self.regionI] = (self.u[:,1])
+        uz[0:len(crack_atoms)-1][self.regionI] = (self.u[:,2])
+
+        mask = np.append(self.regionII,np.array([False],dtype=bool)) #extra false added as we have added alpha to array
+        crack_tip_force_mask = self.regionII
+        if self.extended_far_field:
+            mask = np.append((self.regionII | self.regionIII),np.array([False],dtype=bool))
+            crack_tip_force_mask = self.regionII | self.regionIII
+        
+        falpha, falphas = self.get_crack_tip_force(mask = crack_tip_force_mask,full_array_output=True)
+        logfalphas = np.log10(np.abs(falphas))
+        falphacomponents[mask] = falphas
+        logabsfalphacomp[mask] = logfalphas
+        fx[len(crack_atoms)-1] = falpha
+        ftot[len(crack_atoms)-1] = falpha
+        logabsftot[:] = np.log10(np.abs(ftot))
+        ase.io.write(fname+'.xyz',crack_atoms)
+    
+    def strain_err(self,cutoff,seperate_surface=False):
+        '''Function that returns the atomistic corrector strain error Dv for each atom using the norm of the difference
+        between the corrector on each atom and all those around it within some cutoff. Also has an option to return states
+        adjacent to the free surface of the crack seperately for comparison.
+
+        Parameters
+        ----------
+        cutoff : float
+            cutoff around each atom to use to find the strain error norm.
+        seperate_surface : bool
+            whether or not to return the surface atoms in a seperate array.
+
+        Returns
+        -------
+        r : array of atom radial distance from the centre of the crack system
+
+        dv : norm strain error for each atom in r.
+        """
+        
+        '''
+        #want to get all neighbours in region I
+        I, J = neighbour_list('ij', self.atoms[self.regionI], cutoff)
+        #print(I)
+        #print(J)
+        v = self.u
+        dv = np.linalg.norm(v[I, :] - v[J, :], axis=1)
+        r = self.r[self.regionI][I]
+        mask = r<(self.rI-cutoff)
+        if seperate_surface:
+            mid_cell_x = self.cryst.cell.diagonal()[0] / 2.0
+            mid_cell_y = self.cryst.cell.diagonal()[1] / 2.0
+            I_positions = self.cryst.get_positions()[I]
+            x_criteria = ((I_positions[:,0] - mid_cell_x) < 0.5)
+            y_criteria = ((I_positions[:,1] - mid_cell_y) < ((1.5 * self.cutoff)))&\
+            ((I_positions[:,1] - mid_cell_y) > (-1.5 * self.cutoff))
+            surface_mask_full = np.logical_and(x_criteria,y_criteria)
+            bulk_mask_full = np.logical_not(surface_mask_full)
+            surface_mask = np.logical_and(surface_mask_full,mask)
+            bulk_mask = np.logical_and(bulk_mask_full,mask)
+            return r[surface_mask],dv[surface_mask],r[bulk_mask],dv[bulk_mask]
+        else:
+            return r[mask],dv[mask]
+
+        
+    def convergence_line_plot(self,num=0):
+        def gen_mask(r,rval,dx,dr,cx,cy):
+            #for a line going vertically up
+            #so r is essentiall y here
+            rcriteria = (r > rval)&(r<(rval+dr))
+            xcriteria = (x-cx > 0)&(x-cx<dx)
+            return np.logical_and(rcriteria,xcriteria)
+        crack_atoms = self.atoms.copy()
+        crack_atoms.calc = self.calc
+        forces = crack_atoms.get_forces()
+
+        sx, sy, sz = crack_atoms.cell.diagonal()
+        x, y = crack_atoms.positions[:, 0], crack_atoms.positions[:, 1]
+        cx, cy = sx/2, sy/2
+        r = np.sqrt((x - cx)**2 + (y - cy)**2)
+
+        dr = 0.2 #Angstroms
+        dx = 1 #angstroms
+        regionIvals = np.arange(0,self.rI,dr)
+        other_region_vals = np.arange(self.rI,self.rIII,dr)
+        full_list = np.concatenate((regionIvals,other_region_vals))
+        U = np.zeros([len(regionIvals),3])
+        F = np.zeros([len(full_list),3])
+        for i,rval in enumerate(regionIvals):
+            mask = gen_mask(r,rval,dx,dr,cx,cy)
+            if len(forces[mask])==0:
+                U[i,:] = np.zeros([1,3])
+                F[i,:] = np.zeros([1,3])
+            else:
+                print(f'non 0 at r={rval}')
+                print('vals',self.u[mask[self.regionI]])
+                print('mean',np.mean(self.u[mask[self.regionI]],axis=0))
+                U[i,:] = np.mean(self.u[mask[self.regionI]],axis=0)
+                F[i,:] = np.mean(forces[mask],axis=0)
+
+        for i,rval in enumerate(other_region_vals):
+            mask = gen_mask(r,rval,dx,dr,cx,cy)
+            k = i+len(regionIvals)
+            print(forces[mask])
+            if len(forces[mask])==0:
+                F[k,:] = np.zeros([1,3])
+            else:
+                print(f'non 0 at r={rval}')
+                F[k,:] = np.mean(forces[mask],axis=0)
+
+        plt.figure()
+        for i in range(3):
+            plt.plot(regionIvals,U[:,i])
+        plt.xlabel('r')
+        plt.ylabel('atomistic corrector U at r')
+        plt.legend(['Ux','Uy','Uz'])
+        plt.savefig(f'AtomisticCorrector{num}.png')
+        
+        plt.figure()
+        for i in range(3):
+            plt.plot(full_list,F[:,i])
+        plt.xlabel('r')
+        plt.ylabel('Force on atoms at r')
+        plt.legend(['Fx','Fy','Fz'])
+        plt.savefig(f'ForcesOnAtoms{num}.png')
+
 
 
 def isotropic_modeI_crack_tip_stress_field(K, r, t, xy_only=True,
diff --git a/matscipy/meson.build b/matscipy/meson.build
index bd5a2ba3..a8e6f488 100644
--- a/matscipy/meson.build
+++ b/matscipy/meson.build
@@ -28,13 +28,16 @@ python_sources = [
     'numerical.py',
     'numpy_tricks.py',
     'opls.py',
+    'optimize.py',
     'precon.py',
     'pressurecoupling.py',
     'rings.py',
     'socketcalc.py',
     'spatial_correlation_function.py',
     'surface.py',
-    'visualise.py'
+    'visualise.py',
+    'cauchy_born.py',
+    'surface_reconstruction.py'
 ]
 
 # Install pure Python
diff --git a/matscipy/optimize.py b/matscipy/optimize.py
new file mode 100644
index 00000000..60a4b382
--- /dev/null
+++ b/matscipy/optimize.py
@@ -0,0 +1,197 @@
+
+# This function is temporarily here, for use in NCFlex.
+# Eventually meant to be pushed to ase.optimize.ode
+
+import numpy as np
+from ase.optimize.sciopt import SciPyOptimizer, OptimizerConvergenceError
+
+def ode12r(f, X0, args=None, h=None, verbose=1, fmax=1e-6, maxtol=1e3, steps=100,
+           rtol=1e-1, C1=1e-2, C2=2.0, hmin=1e-10, extrapolate=3,
+           callback=None, apply_precon=None, converged=None, residual=None):
+    """
+    Adaptive ODE solver, which uses 1st and 2nd order approximations to
+    estimate local error and choose a new step length. The original
+    implementation in ase.optimize.ode, and has been modified to this 
+    matscipy version to ensure compatibility with the NCFlex routines.
+
+    This optimizer is described in detail in:
+
+                S. Makri, C. Ortner and J. R. Kermode, J. Chem. Phys.
+                150, 094109 (2019)
+                https://dx.doi.org/10.1063/1.5064465
+
+    Parameters
+    ----------
+
+    f : function
+        function returning driving force on system
+    X0 : 1-dimensional array
+        initial value of degrees of freedom
+    h : float
+        step size, if None an estimate is used based on ODE12
+    verbose: int
+        verbosity level. 0=no output, 1=log output (default), 2=debug output.
+    fmax : float
+        convergence tolerance for residual force
+    maxtol: float
+        terminate if reisdual exceeds this value
+    rtol : float
+        relative tolerance
+    C1 : float
+        sufficient contraction parameter
+    C2 : float
+        residual growth control (Inf means there is no control)
+    hmin : float
+        minimal allowed step size
+    extrapolate : int
+        extrapolation style (3 seems the most robust)
+    callback : function
+        optional callback function to call after each update step
+    apply_precon: function
+        apply a apply_preconditioner to the optimisation
+    converged: function
+        alternative function to check convergence, rather than
+        using a simple norm of the forces.
+    residual: function
+        compute the residual from the current forces
+
+    Returns
+    -------
+
+    X: array
+        final value of degrees of freedom
+    """
+
+    X = X0
+    if args is None:
+        Fn = f(X)
+    else:
+        #x1, xdot1, ds = args 
+        #Fn = f(X, x1, xdot1, ds)
+        Fn = f(X, *args)
+        # For debugging
+        print('Last 4 forces:', Fn[-4:])
+
+    if callback is None:
+        def callback(X):
+            pass
+    callback(X)
+
+    if residual is None:
+        def residual(F, X):
+            return np.linalg.norm(F, np.inf)
+    Rn = residual(Fn, X)
+
+    if apply_precon is None:
+        def apply_precon(F, X):
+            return F, residual(F, X)
+    Fp, Rp = apply_precon(Fn, X)
+    # For debugging
+    if args is not None:
+        print('Last 4 precon forces:',Fp[-4:])
+
+    def log(*args):
+        if verbose >= 1:
+            print(*args)
+
+    def debug(*args):
+        if verbose >= 2:
+            print(*args)
+
+    if converged is None:
+        def converged(F, X):
+            return residual(F, X) <= fmax
+
+    if converged(Fn, X):
+        log("ODE12r terminates successfully after 0 iterations")
+        return X, 0
+    if Rn >= maxtol:
+        raise OptimizerConvergenceError(f"ODE12r: Residual {Rn} is too large "
+                                        "at iteration 0")
+
+    # computation of the initial step
+    r = residual(Fp, X)  # pick the biggest force
+    if h is None:
+        h = 0.5 * rtol ** 0.5 / r  # Chose a stepsize based on that force
+        h = max(h, hmin)  # Make sure the step size is not too big
+
+    for nit in range(1, steps):
+        Xnew = X + h * Fp  # Pick a new position
+        #Fn_new = f(Xnew) if args is None else f(Xnew, x1, xdot1, ds) # Calculate the new forces at this position 
+        Fn_new = f(Xnew) if args is None else f(Xnew, *args) # Calculate the new forces at this position
+        # For debugging
+        if args is not None:
+            print('Last 4 forces:',Fn_new[-4:])
+        Rn_new = residual(Fn_new, Xnew)
+        Fp_new, Rp_new = apply_precon(Fn_new, Xnew)
+        # For debugging
+        if args is not None:
+            print('Last 4 precon forces:',Fp_new[-4:])
+
+        e = 0.5 * h * (Fp_new - Fp)  # Estimate the area under the forces curve
+        err = np.linalg.norm(e, np.inf)  # Error estimate
+
+        # Accept step if residual decreases sufficiently and/or error acceptable
+        accept = ((Rp_new <= Rp * (1 - C1 * h)) or
+                  ((Rp_new <= Rp * C2) and err <= rtol))
+
+        # Pick an extrapolation scheme for the system & find new increment
+        y = Fp - Fp_new
+        if extrapolate == 1:  # F(xn + h Fp)
+            h_ls = h * (Fp @ y) / (y @ y)
+        elif extrapolate == 2:  # F(Xn + h Fp)
+            h_ls = h * (Fp @ Fp_new) / (Fp @ y + 1e-10)
+        elif extrapolate == 3:  # min | F(Xn + h Fp) |
+            h_ls = h * (Fp @ y) / (y @ y + 1e-10)
+        else:
+            raise ValueError(f'invalid extrapolate value: {extrapolate}. '
+                             'Must be 1, 2 or 3')
+        if np.isnan(h_ls) or h_ls < hmin:  # Rejects if increment is too small
+            h_ls = np.inf
+
+        h_err = h * 0.5 * np.sqrt(rtol / err)
+
+        # Accept the step and do the update
+        if accept:
+            X = Xnew
+            Rn = Rn_new
+            Fn = Fn_new
+            Fp = Fp_new
+            Rp = Rp_new
+            callback(X)
+
+            # We check the residuals again
+            if Rn >= maxtol:
+                raise OptimizerConvergenceError(
+                    f"ODE12r: Residual {Rn} is too "
+                    f"large at iteration number {nit}")
+
+            if converged(Fn, X):
+                log(f"ODE12r: terminates successfully "
+                    f"after {nit} iterations.")
+                return X, nit
+
+            # Compute a new step size.
+            # Based on the extrapolation and some other heuristics
+            h = max(0.25 * h,
+                    min(4 * h, h_err, h_ls))  # Log steep-size analytic results
+
+            debug(f"ODE12r:      accept: new h = {h}, |F| = {Rp}")
+            debug(f"ODE12r:                hls = {h_ls}")
+            debug(f"ODE12r:               herr = {h_err}")
+        else:
+            # Compute a new step size.
+            h = max(0.1 * h, min(0.25 * h, h_err,
+                                 h_ls))
+            debug(f"ODE12r:      reject: new h = {h}")
+            debug(f"ODE12r:               |Fnew| = {Rp_new}")
+            debug(f"ODE12r:               |Fold| = {Rp}")
+            debug(f"ODE12r:        |Fnew|/|Fold| = {Rp_new/Rp}")
+
+        # abort if step size is too small
+        if abs(h) <= hmin:
+            raise OptimizerConvergenceError('ODE12r terminates unsuccessfully'
+                                            f' Step size {h} too small')
+
+    raise OptimizerConvergenceError(f'ODE12r terminates unsuccessfully after '
+                                    f'{steps} iterations.')
diff --git a/matscipy/surface.py b/matscipy/surface.py
index 2dd85be6..dc1abee5 100644
--- a/matscipy/surface.py
+++ b/matscipy/surface.py
@@ -26,6 +26,7 @@
 
 import numpy as np
 from numpy.linalg import norm, inv
+from ase import units
 
 def gcd(a, b):
     """Calculate the greatest common divisor of a and b"""
@@ -281,3 +282,56 @@ def make_unit_slab(unit_cell, axes):
 
     sup.set_scaled_positions(sup.get_scaled_positions())
     return sup    
+
+
+def find_surface_energy(symbol,calc,a0,surface,size=(8,1,1),vacuum=10,fmax=0.0001,unit='0.1J/m^2'):
+
+    # Import required lattice builder
+    if surface.startswith('bcc'):
+        from ase.lattice.cubic import BodyCenteredCubic as lattice_builder
+    elif surface.startswith('fcc'):
+        from ase.lattice.cubic import FaceCenteredCubic as lattice_builder #untested
+    elif surface.startswith('diamond'):
+        from ase.lattice.cubic import Diamond as lattice_builder #untested
+    ## Append other lattice builders here
+    else:
+        print('Error: Unsupported lattice ordering.')
+
+    # Set orthogonal directions for cell axes
+    if surface.endswith('100'):
+        directions=[[1,0,0], [0,1,0], [0,0,1]] #tested for bcc
+    elif surface.endswith('110'):
+        directions=[[1,1,0], [-1,1,0], [0,0,1]] #tested for bcc
+    elif surface.endswith('111'):
+        directions=[[1,1,1], [-2,1,1],[0,-1,1]] #tested for bcc
+    ## Append other cell axis options here
+    else:
+        print('Error: Unsupported surface orientation.')
+    
+    # Make bulk and slab with same number of atoms (size)
+    bulk = lattice_builder(directions=directions, size=size, symbol=symbol, latticeconstant=a0, pbc=(1,1,1))
+    cell = bulk.get_cell() ; cell[0,:] *=2 # vacuum along x axis (surface normal)
+    slab = bulk.copy() ; slab.set_cell(cell)
+    
+    # Optimize the geometries
+    from ase.optimize import LBFGSLineSearch
+    bulk.calc = calc ; opt_bulk = LBFGSLineSearch(bulk) ; opt_bulk.run(fmax=fmax)
+    slab.calc = calc ; opt_slab = LBFGSLineSearch(slab) ; opt_slab.run(fmax=fmax)
+
+    # Find surface energy
+    import numpy as np
+    Ebulk = bulk.get_potential_energy() ; Eslab = slab.get_potential_energy()
+    area = np.linalg.norm(np.cross(slab.get_cell()[1,:],slab.get_cell()[2,:]))
+    gamma_ase = (Eslab - Ebulk)/(2*area)
+
+    # Convert to required units
+    if unit == 'ASE':
+        return [gamma_ase,'ase_units']
+    else:
+        gamma_SI = (gamma_ase / units.J ) * (units.m)**2
+        if unit =='J/m^2':
+            return gamma_SI
+        elif unit == '0.1J/m^2':
+            return 10*gamma_SI # units required for the fracture code
+        else:
+            print('Error: Unsupported unit of surface energy.')
diff --git a/matscipy/surface_reconstruction.py b/matscipy/surface_reconstruction.py
new file mode 100644
index 00000000..bcb12a6a
--- /dev/null
+++ b/matscipy/surface_reconstruction.py
@@ -0,0 +1,410 @@
+import numpy as np
+from ase.lattice.cubic import Diamond, FaceCenteredCubic, SimpleCubic, BodyCenteredCubic
+
+from ase.optimize import LBFGS
+
+from ase.constraints import FixAtoms
+
+from matscipy.cauchy_born import CubicCauchyBorn
+from ase.optimize.precon import PreconLBFGS
+import ase
+
+
+class SurfaceReconstruction:
+    """Object for mapping and applying a surface reconstruction in simple and multi-lattices.
+
+    The relaxation of a given free surface of a crystal at a given orientation can be mapped and
+    saved. This relaxation can then be applied to atoms in a different ASE atoms object.
+    This is extremely useful in the case of fracture mechanics, where surface relaxation
+    can be mapped to the internal crack surfaces of an atoms object.
+    """
+
+    def __init__(self, el, a0, calc, directions, surf_dir, lattice='diamond', multilattice=False):
+        """Parameters
+        ----------
+        el : string
+            Element crystal is composed of in ASE format
+        a0 : float
+            Lattice constant of crystal
+        calc : ASE calculator object
+            Calculator for surface relaxation
+        directions : list
+            The x, y and z lab frame directions expressed in the lattice frame.
+            e.g directions=[[1,1,1], [0,-1,1], [-2,1,1]]
+        surf_dir : int
+            Direction which is normal to the free surface. Can be 0, 1 or 2 for
+            lab x, lab y or lab z (as defined by directions).
+        lattice : ASE ase.lattice.cubic function
+            Callable function which can be used to generate the crystal.
+        """
+        # directions is a list of 3 vectors specifiying the surface of interest
+        # in the orientation specified by the problem. surf_dir says which axis is free surface
+        # e.g directions=[[1,1,1], [0,-1,1], [-2,1,1]]
+        self.el = el
+        self.a0 = a0
+        self.calc = calc
+        self.directions = directions
+        self.surf_dir = surf_dir
+        self.inter_planar_dist = a0/np.linalg.norm(directions[surf_dir])
+        # build the rotation matrix from the directions
+        A = np.zeros([3, 3])
+        for i, direction in enumerate(directions):
+            direction = np.array(direction)
+            A[:, i] = direction / np.linalg.norm(direction)
+        self.A = A  # rotation matrix
+        self.eval_cb = None
+        self.lattice = lattice
+        self.cb = None
+        if multilattice:
+            self.cb = CubicCauchyBorn(
+                self.el, self.a0, self.calc, lattice=self.lattice)
+
+    def map_surface(self, fmax=0.0001, layers=6, cutoff=10, shift=0):
+        """Map the relaxation of the crystal surface a certain number of layers deep.
+        Parameters
+        ----------
+        fmax : float
+            Force tolerance for relaxation
+        layers : int
+            Number of layers deep to map the free surface
+        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.
+        """
+        # build a single cell to determine how many atomic surface
+        # layers there are in a single unit cell made by ASE
+        single_cell = self.lattice(
+            directions=self.directions,
+            size=[1, 1, 1],
+            symbol=self.el,
+            latticeconstant=self.a0,
+            pbc=(1, 1, 1))
+        # ase.io.write(f'{self.directions[self.surf_dir]}_unitcell.xyz',single_cell)
+
+        if self.cb is not None:
+            # for a multilattice, get number of layers by just looking
+            # at one of the sublattices
+            self.cb.set_sublattices(single_cell, self.A)
+            n_layer_per_cell = len(np.unique(single_cell.get_positions()[
+                                   self.cb.lattice1mask][:, self.surf_dir].round(decimals=6)))
+        else:
+            # otherwise, just look at all the layers directly
+            n_layer_per_cell = len(np.unique(single_cell.get_positions()[
+                                   :, self.surf_dir].round(decimals=6)))
+
+        # print('nlayer per cell', n_layer_per_cell)
+        # ase.io.write('single_surface_cell.xyz',single_cell)
+
+        # set the number of atomic unit cells to map
+        # as twice the number of layers to map + a number of layers
+        # equal to the potential cutoff for fixing
+        height = 2*int(layers/n_layer_per_cell) + \
+            int(np.round(cutoff/self.inter_planar_dist))
+        size = [1, 1, 1]
+        size[self.surf_dir] *= height
+        self.layers = layers
+
+        # build the bulk structure, and the slab structure which has vacuum in
+        # surface direction
+        bulk = self.lattice(
+            directions=self.directions,
+            size=size,
+            symbol=self.el,
+            latticeconstant=self.a0,
+            pbc=(1, 1, 1))
+
+        # shift the cell slightly and wrap to stop dodgy surfaces when vacuum added
+        # as well as to add on any user-defined shift
+        shift_vec = np.zeros([3])
+        shift_vec[self.surf_dir] += 0.01 + shift
+        bulk.positions += shift_vec
+        bulk.wrap()
+
+        cell = bulk.get_cell()
+        # vacuum along self.surf_dir axis (surface normal)
+        cell[self.surf_dir, :] *= 2
+        slab = bulk.copy()
+        slab.set_cell(cell, scale_atoms=False)
+
+        # get the mask for the bottom surface and fix it
+        pos_before = slab.get_positions()
+        # fix bottom atoms to stop lower surface reconstruction
+        mask = pos_before[:, self.surf_dir] < cutoff
+        slab.set_constraint(FixAtoms(mask=mask))
+        slab.calc = self.calc
+        # print('fbefore',slab.get_forces())
+        # ase.io.write('0_slab.xyz',slab)
+
+        # run an optimisation to relax the surface
+        opt_slab = PreconLBFGS(slab)
+        opt_slab.run(fmax=fmax)
+        # ase.io.write('1_slab.xyz',slab)
+        pos_after = slab.get_positions()
+
+        # measure the shift between the final position and the intial position
+        pos_shift = pos_after - pos_before
+        # print('fafter',slab.get_forces())
+
+        # for a multi-lattice:
+        if self.cb is not None:
+            # if this is a multi-lattice
+
+            # set which atoms are in each multilattice based on the bulk structure. This is important because
+            # different relaxations happen for each multilattice atom in each
+            # atomic layer near the surface
+            self.cb.set_sublattices(bulk, self.A)
+
+            # split the full lattice into the two multilattices
+            total_layers = n_layer_per_cell*height
+            n_atoms_per_layer = int(len(slab) / total_layers)
+
+            # get the indices of the atoms based on how close they are to the
+            # surface in order, for each sub lattice
+            sorted_surf_indices_lattice1 = np.argsort(
+                pos_before[self.cb.lattice1mask][:, self.surf_dir])
+            sorted_surf_indices_lattice2 = np.argsort(
+                pos_before[self.cb.lattice2mask][:, self.surf_dir])
+
+            # split these into layers and save them to an array.
+            # each layer holds the indices of the atoms it contains [note - these are the indices of the array
+            # pos_shift[self.cb.latticemask] - not pos_shift]
+            surf_layer_list_lattice1 = np.zeros(
+                [total_layers, int(n_atoms_per_layer / 2)], dtype=int)
+            surf_layer_list_lattice2 = np.zeros(
+                [total_layers, int(n_atoms_per_layer / 2)], dtype=int)
+            for i in range(total_layers):
+                surf_layer_list_lattice1[i, :] = sorted_surf_indices_lattice1[i * int(
+                    n_atoms_per_layer / 2):(i + 1) * int(n_atoms_per_layer / 2)]
+                surf_layer_list_lattice2[i, :] = sorted_surf_indices_lattice2[i * int(
+                    n_atoms_per_layer / 2):(i + 1) * int(n_atoms_per_layer / 2)]
+
+            # get the inter-layer distance in the direction we care about.
+            self.inter_surf_dist = np.abs(pos_before[self.cb.lattice1mask][surf_layer_list_lattice1[1, :]][0, self.surf_dir]
+                                          - pos_before[self.cb.lattice1mask][surf_layer_list_lattice1[0, :]][0, self.surf_dir])
+
+            # now go through the top layers that we are applying a surface relaxation to and save the shifts of
+            # the lattice 1 and lattice 2 atoms in these layers in the form
+            # [x,y,z] for each layer
+
+            self.relaxation_array_lattice1 = np.zeros([layers, 3])
+            self.relaxation_array_lattice2 = np.zeros([layers, 3])
+            for layer in range(layers):
+                pos_shifts_lattice1 = pos_shift[self.cb.lattice1mask][
+                    surf_layer_list_lattice1[total_layers - layer - 1, :]]
+                pos_shifts_lattice2 = pos_shift[self.cb.lattice2mask][
+                    surf_layer_list_lattice2[total_layers - layer - 1, :]]
+
+                self.relaxation_array_lattice1[layer,
+                                               :] = pos_shifts_lattice1[0, :]
+                self.relaxation_array_lattice2[layer,
+                                               :] = pos_shifts_lattice2[0, :]
+
+        else:
+            # for a non-multilattice
+            # get the total number of in the surface from the number of layers per cell
+            total_layers = n_layer_per_cell*height
+            n_atoms_per_layer = int(len(slab) / total_layers)
+            # print(n_atoms_per_layer)
+            # print(len(slab))
+            # print(pos_shift)
+            # get the indices of the atoms sorted in height order
+            sorted_surf_indices = np.argsort(pos_before[:, self.surf_dir])
+
+            # split these into layers and save them to an array.
+            # each layer holds the indices of the atoms it contains [note - these are the indices of the array
+            # pos_shift[] - not pos_shift]
+            surf_layer_list = np.zeros(
+                [total_layers, n_atoms_per_layer], dtype=int)
+            for i in range(total_layers):
+                surf_layer_list[i, :] = sorted_surf_indices[i *
+                                                            n_atoms_per_layer:((i + 1) * n_atoms_per_layer)]
+            # print(surf_layer_list)
+            self.inter_surf_dist = np.abs(
+                pos_before[surf_layer_list[1, :]][0, self.surf_dir] - pos_before[surf_layer_list[0, :]][0, self.surf_dir])
+            self.relaxation_array = np.zeros([layers, 3])
+            for layer in range(layers):
+                pos_shifts_by_layer = pos_shift[surf_layer_list[total_layers - layer - 1, :]]
+                self.relaxation_array[layer, :] = pos_shifts_by_layer[0, :]
+            # print(self.relaxation_array)
+
+    def identify_layers(
+            self,
+            atoms,
+            surface_coords,
+            xlim=None,
+            ylim=None,
+            zlim=None,
+            search_dir=-1,
+            atoms_for_cb=None):
+        """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.
+        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.
+
+        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.
+        """
+        # function for identifying the layers adjacent to a surface, where the top of the surface
+        # is found near coordinates surface_coords
+        pos = atoms.get_positions()
+
+        # get the masks for the limits imposed by xlim,ylim and zlim
+        if xlim is not None:
+            xmask = (pos[:, 0] < xlim[1]) & (pos[:, 0] > xlim[0])
+        else:
+            xmask = np.ones([len(atoms)], dtype=bool)
+
+        if ylim is not None:
+            ymask = (pos[:, 1] < ylim[1]) & (pos[:, 1] > ylim[0])
+        else:
+            ymask = np.ones([len(atoms)], dtype=bool)
+
+        if zlim is not None:
+            zmask = (pos[:, 2] < zlim[1]) & (pos[:, 2] > zlim[0])
+        else:
+            zmask = np.ones([len(atoms)], dtype=bool)
+
+        lim_mask = xmask & ymask & zmask
+
+        # now, get the mask for the provided atoms structure
+        # for each individual surface layer
+        layer_mask_set = np.zeros([len(atoms), self.layers], dtype=bool)
+        for layer in range(self.layers):
+            # this is the crucial line, essentially finds a mask for the
+            # atoms that lie within a +- 1 layer distance from the starting point.
+            # we then shift down or up a layer each iteration, depending on
+            # search_dir
+            if search_dir == -1:
+                # this corresponds to searching downwards for surface layers
+                layer_mask = np.logical_and((pos[:, self.surf_dir] <= surface_coords[self.surf_dir] - (
+                    layer * self.inter_surf_dist)), (pos[:, self.surf_dir] > surface_coords[self.surf_dir] - ((layer + 1) * self.inter_surf_dist)))
+            elif search_dir == 1:
+                # this corresponds to searching upward for surface layers
+                layer_mask = np.logical_and((pos[:, self.surf_dir] >= surface_coords[self.surf_dir] + (
+                    layer * self.inter_surf_dist)), (pos[:, self.surf_dir] < surface_coords[self.surf_dir] + ((layer + 1) * self.inter_surf_dist)))
+            else:
+                print(
+                    'search_dir should be -1 for a downward surface, 1 for an upward surface!')
+                raise ValueError
+            layer_mask = layer_mask & lim_mask
+
+            # add the layer mask found here to the full set
+            layer_mask_set[:, layer] = layer_mask
+
+        if self.eval_cb is None:
+            return layer_mask_set
+        else:
+            # for a multi-lattice, also have to account for shifts between multilattice atoms
+            # in layers, so need to refine the masks further to be the lattice1
+            # and lattice2 atoms in each layer
+            if atoms_for_cb is None:
+                atoms_for_cb = atoms
+            # note that to find the sublattices, need to use a version of the
+            # atoms structure with no free surface
+            if self.eval_cb.lattice1mask is None:
+                self.eval_cb.set_sublattices(atoms_for_cb, self.A)
+            layer_mask_set_lattice1 = np.zeros_like(layer_mask_set, dtype=bool)
+            layer_mask_set_lattice2 = np.zeros_like(layer_mask_set, dtype=bool)
+            for layer in range(self.layers):
+                layer_mask_set_lattice1[:,
+                                        layer] = layer_mask_set[:,
+                                                                layer] & self.eval_cb.lattice1mask
+                layer_mask_set_lattice2[:,
+                                        layer] = layer_mask_set[:,
+                                                                layer] & self.eval_cb.lattice2mask
+            return layer_mask_set_lattice1, layer_mask_set_lattice2
+
+    def apply_surface_shift(
+            self,
+            atoms,
+            surface_coords,
+            cb=None,
+            xlim=None,
+            ylim=None,
+            zlim=None,
+            search_dir=-1,
+            atoms_for_cb=None):
+        """Function which applies the mapped out 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.
+        """
+        # This finds the layers of atoms that need shifting and applies the
+        # surface shift accordingly.
+        self.eval_cb = cb
+        # print(self.eval_cb)
+        if self.eval_cb is not None:
+            layer_mask_set_lattice1, layer_mask_set_lattice2 = self.identify_layers(
+                atoms, surface_coords, xlim=xlim, ylim=ylim, zlim=zlim, search_dir=search_dir, atoms_for_cb=atoms_for_cb)
+            for layer in range(self.layers):
+                # print('diff',self.relaxation_array_lattice1[layer,:])
+                # print('initial',atoms.positions[:, :][layer_mask_set_lattice1[:,layer]])
+                lattice_1_shift = self.relaxation_array_lattice1[layer, :]
+                lattice_2_shift = self.relaxation_array_lattice2[layer, :]
+
+                # if the surface is going upward, invert surf_dir component of
+                # applied shift
+                lattice_1_shift[self.surf_dir] = (-search_dir) * \
+                    lattice_1_shift[self.surf_dir]
+                lattice_2_shift[self.surf_dir] = (-search_dir) * \
+                    lattice_2_shift[self.surf_dir]
+                atoms.positions[:,
+                                :][layer_mask_set_lattice1[:,
+                                                           layer]] += lattice_1_shift
+                atoms.positions[:,
+                                :][layer_mask_set_lattice2[:,
+                                                           layer]] += lattice_2_shift
+                # print('after',atoms.positions[:, :][layer_mask_set_lattice1[:,layer]])
+
+        else:
+            layer_mask_set = self.identify_layers(
+                atoms,
+                surface_coords,
+                xlim=xlim,
+                ylim=ylim,
+                zlim=zlim,
+                search_dir=search_dir)
+            for layer in range(self.layers):
+                shift = self.relaxation_array[layer, :]
+                shift[self.surf_dir] = (-search_dir) * shift[self.surf_dir]
+                atoms.positions[:, :][layer_mask_set[:, layer]
+                                      ] += shift
diff --git a/scripts/fracture_mechanics/sinclair_continuation.py b/scripts/fracture_mechanics/sinclair_continuation.py
index 7364a20a..081d0da4 100644
--- a/scripts/fracture_mechanics/sinclair_continuation.py
+++ b/scripts/fracture_mechanics/sinclair_continuation.py
@@ -27,6 +27,7 @@
 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 matscipy import parameter
@@ -34,15 +35,16 @@
 from matscipy.fracture_mechanics.crack import CubicCrystalCrack, SinclairCrack
 
 from scipy.optimize import fsolve, fminbound
+from scipy.optimize.nonlin import NoConvergence
 
 import sys
 sys.path.insert(0, '.')
 
 import params
-
 calc = parameter('calc')
 fmax = parameter('fmax', 1e-3)
-max_steps = parameter('max_steps', 100)
+max_opt_steps = parameter('max_opt_steps', 100)
+max_arc_steps = parameter('max_arc_steps', 10)
 vacuum = parameter('vacuum', 10.0)
 flexible = parameter('flexible', True)
 continuation = parameter('continuation', False)
@@ -56,45 +58,58 @@
 precon = parameter('precon', False)
 prerelax = parameter('prerelax', 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')
 
-# compute elastic constants
+
+# make copies of initial configs
 cryst = params.cryst.copy()
+cluster = params.cluster.copy()
+
+# compute elastic constants
 cryst.pbc = True
 cryst.calc = calc
-C, C_err = fit_elastic_constants(cryst,
-                                 symmetry=parameter('elastic_symmetry',
-                                                    'triclinic'))
+C, C_err = fit_elastic_constants(cryst, symmetry=parameter('elastic_symmetry','triclinic'))
 
+# setup the crack
 crk = CubicCrystalCrack(parameter('crack_surface'),
                         parameter('crack_front'),
                         Crot=C/GPa)
 
-# Get Griffith's k1.
+# get Griffith's stess intensity factor
 k1g = crk.k1g(parameter('surface_energy'))
 print('Griffith k1 = %f' % k1g)
 
-cluster = params.cluster.copy()
-
-if continuation and not os.path.exists(traj_file):
+# check for restart files if continuation=True
+if continuation and not os.path.exists(restart_file):
     continuation = False
-
 if continuation:
-    with h5py.File(traj_file, 'r') as hf:
+
+    # works only if restart file contains atleast two consecutive flex solutions
+    with h5py.File(restart_file, 'r') as hf:
         x = hf['x']
         restart_index = parameter('restart_index', x.shape[0] - 1)
         x0 = x[restart_index-1, :]
         x1 = x[restart_index, :]
-        k0 = x0[-1] / k1g
+        k0 = x0[-1] / k1g ; k1 = x1[-1] / k1g
+        alpha0 = x0[-2] / k1g ; alpha1 = x1[-2] / k1g
+        # print restart info
+        print(f'Restarting from k0={k0}, k1={k1} --> alpha0={alpha0}, alpha1={alpha1}')
+
+else:
 
-sc = SinclairCrack(crk, cluster, calc, k0 * k1g,
-                   alpha=alpha0,
+    # setup Sinclair boundary conditions with variable_alpha and no variable_k, for approx solution
+    sc = SinclairCrack(crk, cluster, calc, 
+                   k0*k1g, alpha=alpha0,
                    vacuum=vacuum,
                    variable_alpha=flexible,
                    extended_far_field=extended_far_field)
 
-if not continuation:
     traj = None
     dk = parameter('dk', 1e-4)
     dalpha = parameter('dalpha', 1e-3)
@@ -105,33 +120,43 @@
         a = ase.io.read(f'k_{int(k0 * 1000):04d}.xyz')
         sc.set_atoms(a)
 
+    # setup mask fpr relevant regions
     mask = sc.regionI | sc.regionII
     if extended_far_field:
         mask = sc.regionI | sc.regionII | sc.regionIII
 
-    # first use the CLE-only approximation`: define a function f_alpha0(k, alpha)
-    def f(k, alpha):
-        sc.k = k * k1g
-        sc.alpha = alpha
-        sc.update_atoms()
-        return sc.get_crack_tip_force(mask=mask)
-
-    # identify approximate range of stable k
-    alpha_range = parameter('alpha_range', np.linspace(-a0, a0, 100))
-    k = k0  # initial guess for k
-    alpha_k = []
-    for alpha in alpha_range:
-        # look for solution to f(k, alpha) = 0 close to alpha = alpha
-        (k,) = fsolve(f, k, args=(alpha,))
-        print(f'alpha={alpha:.3f} k={k:.3f} ')
-        alpha_k.append((alpha, k))
-    alpha_k = np.array(alpha_k)
-    kmin, kmax = alpha_k[:, 1].min(), alpha_k[:, 1].max()
-    print(f'Estimated stable K range is {kmin} < k / k_G < {kmax}')
+    
+    # estimate the stable K range
+    k_range= parameter('k_range', 'Unknown') 
+    if isinstance(k_range,list) and len(k_range)==2:
+        kmin, kmax = k_range
+        print(f'Using K range from params file: {kmin} < k / k_G < {kmax}')
+    else:
+        print('No k_range=[kmin,kmax] given in params.')
+        print('Running CLE-only approximation to estimate stable K range.')
+        # first use the CLE-only approximation`: define a function f_alpha0(k, alpha)
+        def f(k, alpha):
+            sc.k = k * k1g
+            sc.alpha = alpha
+            sc.update_atoms()
+            return sc.get_crack_tip_force(mask=mask)
+
+        # identify approximate range of stable k
+        alpha_range = parameter('alpha_range', np.linspace(-a0, a0, 20))
+        k = k0  # initial guess for k
+        alpha_k = []
+        for alpha in alpha_range:
+            # look for solution to f(k, alpha) = 0 close to alpha = alpha
+            (k,) = fsolve(f, k, args=(alpha,))
+            print(f'alpha={alpha:.3f} k={k:.3f} ')
+            alpha_k.append((alpha, k))
+        alpha_k = np.array(alpha_k)
+        kmin, kmax = alpha_k[:, 1].min(), alpha_k[:, 1].max()
+        print(f'Estimated stable K range is {kmin} < k / k_G < {kmax}')
 
     # define a function to relax with static scheme at a given value of k
-    # note that we use a looser fmax of 1e-3 and reduced max steps of 50
-    def g(k):
+    traj = open("traj.xyz", "w")
+    def g(k, do_abs=True):
         print(f'Static minimisation with k={k}, alpha={alpha0}.')
         sc.k = k * k1g
         sc.alpha = alpha0
@@ -142,56 +167,70 @@ def g(k):
         atoms.calc = sc.calc
         atoms.set_constraint(FixAtoms(mask=~sc.regionI))
         opt = PreconLBFGS(atoms, logfile=None)
-        opt.run(fmax=1e-5)
+        #opt = LBFGSLineSearch(atoms)
+        opt.run(fmax=1e-3)
+        atoms.write(traj, format="extxyz")
         sc.set_atoms(atoms)
         f_alpha = sc.get_crack_tip_force(mask=mask)
         print(f'Static minimisation with k={k}, alpha={alpha0} --> f_alpha={f_alpha}')
-        return abs(f_alpha)
+        if do_abs:
+            f_alpha = abs(f_alpha)
+        return f_alpha
 
     # minimise g(k) in [kmin, kmax]
     kopt, falpha_min, ierr, funccalls = fminbound(g, kmin, kmax, xtol=1e-8, full_output=True)
     print(f'Brent minimisation yields f_alpha={falpha_min} for k = {kopt} after {funccalls} calls')
+    traj.close()
 
-    # re-optimize, first with static scheme, then flexible scheme
-    sc.k = kopt * k1g
+    # Find flex1
+    # first optimize with static scheme
+    #sc.k = kopt * k1g
+    k0 = kopt
+    sc.rescale_k(k0 * k1g)
     sc.alpha = alpha0
     sc.variable_alpha = False
-    sc.optimize(ftol=1e-3, steps=max_steps)
-
-    # finally, we revert to target fmax precision
-    sc.variable_alpha = flexible
-    sc.optimize(ftol=fmax, steps=max_steps)
-
+    sc.optimize(ftol=1e-3, steps=max_opt_steps,method=opt_method)
+    # then revert to target fmax precision and optimize with flexible scheme
+    sc.variable_alpha = flexible #True
+    sc.optimize(ftol=fmax, steps=max_opt_steps,method=opt_method)
+    # save flex1
     sc.atoms.write('x0.xyz')
     x0 = np.r_[sc.get_dofs(), k0 * k1g]
     alpha0 = sc.alpha
 
-    # obtain a second solution x1 = (u_1, alpha_1, k_1) where
+    # Find flex2: a second solution x1 = (u_1, alpha_1, k_1) where
     # k_1 ~= k_0 and alpha_1 ~= alpha_0
-
+    # Increase k, and rescale Us accordingly
     print(f'Rescaling K_I from {sc.k} to {sc.k + dk * k1g}')
     k1 = k0 + dk
     sc.rescale_k(k1 * k1g)
-
-
+    # optimize at target fmax precision with flexible scheme
+    sc.variable_alpha = flexible #True
+    sc.optimize(ftol=fmax, steps=max_opt_steps,method=opt_method)
+    # save flex2
     sc.atoms.write('x1.xyz')
     x1 = np.r_[sc.get_dofs(), k1 * k1g]
-    # check crack tip didn't jump too far
     alpha1 = sc.alpha
+
+    # check crack tip didn't jump too far
     print(f'k0={k0}, k1={k1} --> alpha0={alpha0}, alpha1={alpha1}')
     assert abs(alpha1 - alpha0) < dalpha
 
-scv = SinclairCrack(crk, cluster, calc, k0 * k1g,
-                    alpha=alpha0,
+# setup new crack with variable_k for full solution
+scv = SinclairCrack(crk, cluster, calc, 
+                    k0*k1g, alpha=alpha0,
                     vacuum=vacuum,
-                    variable_alpha=flexible,
-                    variable_k=True,
+                    variable_alpha=flexible, variable_k=True,
                     extended_far_field=extended_far_field)
 
+# run full arc length continuation
 scv.arc_length_continuation(x0, x1, N=nsteps,
-                            ds=ds, ftol=fmax, steps=max_steps,
-                            direction=direction,
-                            continuation=continuation,
-                            traj_file=traj_file,
-                            traj_interval=traj_interval,
-                            precon=precon)
+                        ds=ds, ftol=fmax, max_steps=max_arc_steps,
+                        direction=direction,
+                        continuation=continuation,
+                        traj_file=traj_file,
+                        traj_interval=traj_interval,
+                        precon=precon,
+                        ds_max=ds_max, ds_min=ds_min,
+                        ds_aggressiveness=ds_aggressiveness,
+                        opt_method='krylov') # 'ode12r' preconditioning needs debugging
\ No newline at end of file


Test patch:
diff --git a/tests/test_cauchy_born_corrector.py b/tests/test_cauchy_born_corrector.py
new file mode 100644
index 00000000..606604cb
--- /dev/null
+++ b/tests/test_cauchy_born_corrector.py
@@ -0,0 +1,238 @@
+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 scipy.linalg import sqrtm
+import ase.io
+import matscipytest
+
+
+class TestPredictCauchyBornShifts(matscipytest.MatSciPyTestCase):
+    """
+    Tests of Cauchy-Born shift prediction in the case of multilattices.
+
+    We test with a cubic Silicon diamond multi-lattice
+    """
+
+    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
+        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
+
+    def test_set_sublattice(self):
+        A = np.eye(3)
+        self.cb.set_sublattices(self.unit_cell, A)
+        assert np.array_equal(self.cb.lattice1mask, np.array(
+            [True, False, True, False, True, False, True, False], dtype=bool))
+        assert np.array_equal(self.cb.lattice2mask, np.array(
+            [False, True, False, True, False, True, False, True], dtype=bool))
+
+    def test_fit_taylor_model(self):
+        self.cb.fit_taylor()
+        # print(self.cb.grad_f)
+        # print(self.cb.hess_f)
+        grad_f = np.array([0.00000000e+00, 0.00000000e+00, 0.00000000e+00, -
+                          1.09995082e+00, -2.22044605e-12, -3.33066907e-12])
+        hess_f = np.array([[0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 2.03958794e-01, 0.00000000e+00, 1.11022302e-08],
+                           [0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
+                               9.12311482e-01, -5.55111512e-09, 0.00000000e+00],
+                           [0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
+                               9.12311465e-01, 0.00000000e+00, 0.00000000e+00],
+                           [2.03958794e-01, 9.12311493e-01, 9.12311465e-01,
+                            0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
+                           [0.00000000e+00, -5.55111512e-09, 0.00000000e+00,
+                            0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
+                           [1.11022302e-08, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00]])
+        assert np.allclose(self.cb.grad_f, grad_f, atol=1e-6)
+        assert np.allclose(self.cb.hess_f, hess_f, atol=1e-6)
+
+    def E_cart3D(self, x, y, z, eps=None):
+        E = np.zeros([len(x), 3, 3])
+        E[:, 0, 0] = eps[0]
+        E[:, 1, 1] = eps[1]
+        E[:, 2, 2] = eps[2]
+        E[:, 1, 2], E[:, 2, 1] = eps[3], eps[3]
+        E[:, 0, 2], E[:, 2, 0] = eps[4], eps[4]
+        E[:, 0, 1], E[:, 1, 0] = eps[5], eps[5]
+        return E
+
+    def E_cylind3D(self, r, theta, z, eps=None):
+        E = np.zeros([len(r), 3, 3])
+        E[:, 0, 0] = eps[0]
+        E[:, 1, 1] = eps[1]
+        E[:, 2, 2] = eps[2]
+        E[:, 1, 2], E[:, 2, 1] = eps[3], eps[3]
+        E[:, 0, 2], E[:, 2, 0] = eps[4], eps[4]
+        E[:, 0, 1], E[:, 1, 0] = eps[5], eps[5]
+        return E
+
+    def F_cart3D(self, x, y, z, eps=None):
+        E = np.zeros([len(x), 3, 3])
+        E[:, 0, 0] = eps[0]
+        E[:, 1, 1] = eps[1]
+        E[:, 2, 2] = eps[2]
+        E[:, 1, 2], E[:, 2, 1] = eps[3], eps[3]
+        E[:, 0, 2], E[:, 2, 0] = eps[4], eps[4]
+        E[:, 0, 1], E[:, 1, 0] = eps[5], eps[5]
+        U = np.zeros_like(E)
+        for i in range(np.shape(E)[0]):
+            U[i, :, :] = sqrtm((2*E[i, :, :])+np.eye(3))
+        return U  # with no rigid rotation, U is F
+
+    def F_cylind3D(self, r, theta, z, eps=None):
+        E = np.zeros([len(r), 3, 3])
+        E[:, 0, 0] = eps[0]
+        E[:, 1, 1] = eps[1]
+        E[:, 2, 2] = eps[2]
+        E[:, 1, 2], E[:, 2, 1] = eps[3], eps[3]
+        E[:, 0, 2], E[:, 2, 0] = eps[4], eps[4]
+        E[:, 0, 1], E[:, 1, 0] = eps[5], eps[5]
+        U = np.zeros_like(E)
+        for i in range(np.shape(E)[0]):
+            U[i, :, :] = sqrtm((2*E[i, :, :])+np.eye(3))
+        return U  # with no rigid rotation, U is F
+
+    def model_prediction(self, dirs, eps, method, E_func=None, F_func=None, coordinates='cart3D', atol=1.5e-5, returnvals=False):
+        # print(dirs,coordinates)
+        # Get rotation matrix from dirs (A)
+        A = np.zeros([3, 3])
+        for i, direction in enumerate(dirs):
+            direction = np.array(direction)
+            direction = (1/np.linalg.norm(direction))*direction
+            A[:, i] = direction
+
+        # generate rotated unitcell (according to dirs)
+        atoms = Diamond(self.cb.el, latticeconstant=self.cb.a0,
+                        size=[1, 1, 1], directions=dirs)
+        atoms.set_pbc([True, True, True])
+        atoms.calc = self.cb.calc
+        atoms_copy = atoms.copy()
+        self.cb.set_sublattices(atoms, A)
+        shift_diff_true = self.cb.eval_shift(
+            eps, atoms)  # lab frame shift diff
+        if E_func is not None:
+            shifts = self.cb.predict_shifts(
+                A, atoms, E_func=E_func, eps=eps, method=method, coordinates=coordinates)
+            shift_diff_predicted = shifts[0, :]
+        elif F_func is not None:
+            shifts = self.cb.predict_shifts(
+                A, atoms, F_func=F_func, eps=eps, method=method, coordinates=coordinates)
+            shift_diff_predicted = shifts[0, :]
+        # print(shift_diff_predicted,shift_diff_true)
+        # print(shift_diff_predicted-shift_diff_true)
+        assert np.allclose(shift_diff_predicted, shift_diff_true, atol=atol)
+        if returnvals:
+            return atoms_copy, shifts, np.linalg.norm(shift_diff_predicted-shift_diff_true), A
+
+    def test_taylor_model_E(self):
+        self.cb.fit_taylor()
+        eps_vals = [np.array([0, 0, 0, 0.0001, 0, 0]),
+                    np.array([0.0001, 0, 0.0001, 0.0001, 0.0001, 0]),
+                    np.array([0, 0, 0, 0.0001, 0, 0]),
+                    np.array([0, 0, 0, 0.0001, 0, 0])]
+
+        dir_vals = [[[1, 0, 0], [0, 1, 0], [0, 0, 1]],
+                    [[1, 1, 0], [-1, 1, 0], np.cross([1, 1, 0], [-1, 1, 0])],
+                    [[1, 1, 1], [-2, 1, 1], np.cross([1, 1, 1], [-2, 1, 1])],
+                    [[1, 1, 0], [-1, 1, 0], np.cross([1, 1, 0], [-1, 1, 0])]]
+
+        funcs = [self.E_cart3D, self.E_cart3D, self.E_cart3D, self.E_cylind3D]
+        coords = ['cart3D', 'cart3D', 'cart3D', 'cylind3D']
+        for i in range(len(funcs)):
+            self.model_prediction(
+                dir_vals[i], eps_vals[i], method='taylor', E_func=funcs[i], coordinates=coords[i], atol=1e-7)
+
+    def test_taylor_model_F(self):
+        self.cb.fit_taylor()
+        eps_vals = [np.array([0, 0, 0, 0.0001, 0, 0]),
+                    np.array([0.0001, 0, 0.0001, 0.0001, 0.0001, 0]),
+                    np.array([0, 0, 0, 0.0001, 0, 0]),
+                    np.array([0, 0, 0, 0.0001, 0, 0])]
+
+        dir_vals = [[[1, 0, 0], [0, 1, 0], [0, 0, 1]],
+                    [[1, 1, 0], [-1, 1, 0], np.cross([1, 1, 0], [-1, 1, 0])],
+                    [[1, 1, 1], [-2, 1, 1], np.cross([1, 1, 1], [-2, 1, 1])],
+                    [[1, 1, 0], [-1, 1, 0], np.cross([1, 1, 0], [-1, 1, 0])]]
+
+        funcs = [self.F_cart3D, self.F_cart3D, self.F_cart3D, self.F_cylind3D]
+        coords = ['cart3D', 'cart3D', 'cart3D', 'cylind3D']
+        for i in range(len(funcs)):
+            self.model_prediction(
+                dir_vals[i], eps_vals[i], method='taylor', F_func=funcs[i], coordinates=coords[i], atol=1e-7)
+
+    def test_fit_regression_model(self):
+        self.cb.initial_regression_fit()
+
+    def test_regression_model_E(self):
+        self.cb.initial_regression_fit()
+        eps_vals = [np.array([0, 0, 0, 0.0001, 0, 0]),
+                    np.array([0.0001, 0, 0.0001, 0.0001, 0.0001, 0]),
+                    np.array([0, 0, 0, 0.0001, 0, 0]),
+                    np.array([0, 0, 0, 0.0001, 0, 0])]
+
+        dir_vals = [[[1, 0, 0], [0, 1, 0], [0, 0, 1]],
+                    [[1, 1, 0], [-1, 1, 0], np.cross([1, 1, 0], [-1, 1, 0])],
+                    [[1, 1, 1], [-2, 1, 1], np.cross([1, 1, 1], [-2, 1, 1])],
+                    [[1, 1, 0], [-1, 1, 0], np.cross([1, 1, 0], [-1, 1, 0])]]
+
+        funcs = [self.E_cart3D, self.E_cart3D, self.E_cart3D, self.E_cylind3D]
+        coords = ['cart3D', 'cart3D', 'cart3D', 'cylind3D']
+        for i in range(len(funcs)):
+            self.model_prediction(
+                dir_vals[i], eps_vals[i], method='regression', E_func=funcs[i], coordinates=coords[i], atol=1e-6)
+
+    def test_regression_model_F(self):
+        self.cb.initial_regression_fit()
+        eps_vals = [np.array([0, 0, 0, 0.0001, 0, 0]),
+                    np.array([0.0001, 0, 0.0001, 0.0001, 0.0001, 0]),
+                    np.array([0, 0, 0, 0.0001, 0, 0]),
+                    np.array([0, 0, 0, 0.0001, 0, 0])]
+
+        dir_vals = [[[1, 0, 0], [0, 1, 0], [0, 0, 1]],
+                    [[1, 1, 0], [-1, 1, 0], np.cross([1, 1, 0], [-1, 1, 0])],
+                    [[1, 1, 1], [-2, 1, 1], np.cross([1, 1, 1], [-2, 1, 1])],
+                    [[1, 1, 0], [-1, 1, 0], np.cross([1, 1, 0], [-1, 1, 0])]]
+
+        funcs = [self.F_cart3D, self.F_cart3D, self.F_cart3D, self.F_cylind3D]
+        coords = ['cart3D', 'cart3D', 'cart3D', 'cylind3D']
+        for i in range(len(funcs)):
+            self.model_prediction(
+                dir_vals[i], eps_vals[i], method='regression', F_func=funcs[i], coordinates=coords[i], atol=1e-6)
+
+    def test_regression_model_refit(self):
+        self.cb.initial_regression_fit()
+        eps = np.array([0.01, 0, 0.01, 0.01, 0, -0.01])
+        dirs = [[1, 1, 1], [-2, 1, 1], np.cross([1, 1, 1], [-2, 1, 1])]
+        func = self.E_cart3D
+        coordinates = 'cart3D'
+        atoms, shifts, shift_err_before, A = self.model_prediction(
+            dirs, eps, method='regression', E_func=func, coordinates=coordinates, atol=1e-4, returnvals=True)
+        atoms_copy = atoms.copy()
+        atoms_copy.calc = self.cb.calc
+        self.cb.apply_shifts(atoms_copy, shifts)
+        forces_after = atoms_copy.get_forces()
+        converged, err = self.cb.check_for_refit(
+            A, atoms_copy, forces_after, E_func=func, eps=eps, tol=1e-4, refit_points=2)
+        atoms, shifts, shift_err_after, A = self.model_prediction(
+            dirs, eps, method='regression', E_func=func, coordinates=coordinates, atol=1e-4, returnvals=True)
+        assert shift_err_after < shift_err_before
+        # print(shift_err_before,shift_err_after)
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/tests/test_surface_reconstruction.py b/tests/test_surface_reconstruction.py
new file mode 100644
index 00000000..2141284d
--- /dev/null
+++ b/tests/test_surface_reconstruction.py
@@ -0,0 +1,180 @@
+import unittest
+import numpy as np
+from matscipy.cauchy_born import CubicCauchyBorn
+from ase.build import bulk
+from ase.lattice.cubic import Diamond, FaceCenteredCubic
+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.explicit_forms.tersoff_brenner import TersoffBrenner, Erhart_PRB_71_035211_Si
+from matscipy.calculators.manybody import Manybody
+from matscipy.surface_reconstruction import SurfaceReconstruction
+import ase
+from matscipy.calculators.eam import EAM
+import matscipytest
+
+
+class TestSurfaceReconstructionMap(matscipytest.MatSciPyTestCase):
+    """
+    Tests of Cauchy-Born shift prediction in the case of multilattices.
+
+    Test uses EAM for gold (Au-Grochola-JCP05.eam.alloy)
+    and 
+    """
+
+    def setupmultilattice(self):
+        el = 'Si'  # chemical element
+        unitcell = Diamond(el, latticeconstant=5.43, size=[1, 1, 1],
+                           directions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]])
+        # model = "c_v3.yace"
+        # d2_table = "d2.table"
+        # cmds = ['pair_style hybrid/overlay pace table linear 10000',
+        # 'pair_coeff * * pace c_v3.yace C', 'pair_coeff * * table d2.table D2 9.0']
+        # calc = LAMMPSlib(amendments=["mass 1 12.101"],lmpcmds = cmds,log_file='test.log',keep_alive=True)
+        calc = Manybody(**TersoffBrenner(Erhart_PRB_71_035211_Si))
+        unitcell.calc = calc
+        ecf = ExpCellFilter(unitcell)
+        opt = PreconLBFGS(ecf)
+        opt.run(fmax=0.0001, smax=0.0001)
+        a0 = unitcell.cell[0, 0]
+        self.cb = CubicCauchyBorn(el, a0, calc, lattice=Diamond)
+        self.el = el
+        self.a0 = a0
+        self.calc = calc
+
+    def test_surface_mapping_multilattice(self):
+        self.setupmultilattice()
+        dir_vals = [[[1, 1, 0], [-1, 1, 0], [0, 0, 1]],
+                    [[1, 1, 1], [-2, 1, 1], [0, -1, 1]]]
+        surf_dirs = [0, 1, 2]
+        for dirs in dir_vals:
+            for surf_dir in surf_dirs:
+                # print(f'testing relaxation of {dirs[surf_dir]} surface')
+                sr = SurfaceReconstruction(
+                    self.el, self.a0, self.calc, dirs, surf_dir, lattice=Diamond, multilattice=True)
+                sr.map_surface(fmax=0.0005, shift=0)
+                self.check_surface_mapping_multilattice(
+                    sr, dirs, Diamond, surf_dir)
+
+    def check_surface_mapping_multilattice(self, sr, directions, lattice, surf_dir, shift=0):
+        cb = CubicCauchyBorn(self.el, self.a0, self.calc, lattice=lattice)
+
+        size = [1, 1, 1]
+        size[surf_dir] *= 20
+
+        shift_vec = np.zeros([3])
+        shift_vec[surf_dir] += 0.01 + shift
+
+        bulk = lattice(directions=directions, size=size,
+                       symbol=self.el, latticeconstant=self.a0, pbc=(1, 1, 1))
+        calc = self.calc
+        bulk.positions += 0.01+shift
+        bulk.wrap()
+
+        cell = bulk.get_cell()
+        cell[surf_dir, :] *= 2  # vacuum along surface axis (surface normal)
+        slab = bulk.copy()
+        slab.set_cell(cell, scale_atoms=False)
+        pos = slab.get_positions()
+
+        sorted_surf_indices = np.argsort(pos[:, surf_dir])
+
+        # apply the surface reconstruction to the top surface only
+        surface_atom = np.argmax(pos[:, surf_dir])
+        surface_coords = pos[surface_atom, :]+np.array([0.01, 0.01, 0.01])
+
+        slab.set_calculator(calc)
+        initial_force = slab.get_forces()
+        # ase.io.write(f'{directions[surf_dir]}_before.xyz',slab)
+        sr.apply_surface_shift(slab, surface_coords, cb=cb, xlim=None,
+                               ylim=None, zlim=None, search_dir=-1, atoms_for_cb=bulk)
+        final_force = slab.get_forces()
+        # ase.io.write(f'{directions[surf_dir]}_after.xyz',slab)
+
+        # get the atoms nearest to the top surface
+        # apply the surface reconstruction to the bottom surface
+        natoms = len(slab)
+        maxfbefore = np.max(
+            np.abs(initial_force[sorted_surf_indices[int(natoms/2):natoms]]))
+        # print('max force before:',maxfbefore)
+        maxfafter = np.max(
+            np.abs(final_force[sorted_surf_indices[int(natoms/2):natoms]]))
+        # print('max force after:',maxfafter)
+        assert maxfafter < (maxfbefore/10)
+
+    def setupsinglelattice(self):
+        el = 'Au'  # chemical element
+        unitcell = FaceCenteredCubic(el, latticeconstant=4.07, size=[1, 1, 1],
+                                     directions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]])
+        calc = EAM('Au-Grochola-JCP05.eam.alloy')
+        unitcell.calc = calc
+        ecf = ExpCellFilter(unitcell)
+        opt = PreconLBFGS(ecf)
+        opt.run(fmax=0.0001, smax=0.0001)
+        a0 = unitcell.cell[0, 0]
+        self.el = el
+        self.a0 = a0
+        self.calc = calc
+
+    def test_surface_mapping_single_lattice(self):
+        self.setupsinglelattice()
+        dir_vals = [[[1, 1, 0], [-1, 1, 0], [0, 0, 1]],
+                    [[1, 1, 1], [-2, 1, 1], [0, -1, 1]]]
+        surf_dirs = [0, 1, 2]
+        for dirs in dir_vals:
+            for surf_dir in surf_dirs:
+                print(f'testing relaxation of {dirs[surf_dir]} surface')
+                sr = SurfaceReconstruction(
+                    self.el, self.a0, self.calc, dirs, surf_dir, lattice=FaceCenteredCubic)
+                sr.map_surface(fmax=0.0005, shift=0, layers=20)
+                self.check_surface_mapping_single_lattice(
+                    sr, dirs, FaceCenteredCubic, surf_dir)
+
+    def check_surface_mapping_single_lattice(self, sr, directions, lattice, surf_dir, shift=0):
+        size = [1, 1, 1]
+        size[surf_dir] *= 20
+
+        shift_vec = np.zeros([3])
+        shift_vec[surf_dir] += 0.01 + shift
+
+        bulk = lattice(directions=directions, size=size,
+                       symbol=self.el, latticeconstant=self.a0, pbc=(1, 1, 1))
+        calc = self.calc
+        bulk.positions += 0.01+shift
+        bulk.wrap()
+
+        cell = bulk.get_cell()
+        cell[surf_dir, :] *= 2  # vacuum along surface axis (surface normal)
+        slab = bulk.copy()
+        slab.set_cell(cell, scale_atoms=False)
+        pos = slab.get_positions()
+
+        sorted_surf_indices = np.argsort(pos[:, surf_dir])
+
+        # apply the surface reconstruction to the top surface only
+        surface_atom = np.argmax(pos[:, surf_dir])
+        surface_coords = pos[surface_atom, :]+np.array([0.01, 0.01, 0.01])
+
+        slab.set_calculator(calc)
+        initial_force = slab.get_forces()
+        # ase.io.write(f'{directions[surf_dir]}0.xyz',slab)
+        sr.apply_surface_shift(slab, surface_coords,
+                               xlim=None, ylim=None, zlim=None, search_dir=-1)
+        final_force = slab.get_forces()
+        # ase.io.write(f'{directions[surf_dir]}1.xyz',slab)
+
+        # get the atoms nearest to the top surface
+        # apply the surface reconstruction to the bottom surface
+        natoms = len(slab)
+        maxfbefore = np.max(
+            np.abs(initial_force[sorted_surf_indices[int(natoms/2):natoms]]))
+        print('max force before:', maxfbefore)
+        maxfafter = np.max(
+            np.abs(final_force[sorted_surf_indices[int(natoms/2):natoms]]))
+        print('max force after:', maxfafter)
+        assert maxfafter < (maxfbefore/10)
+
+
+if __name__ == '__main__':
+    unittest.main()