diff --git a/docs/cli/interatomic.md b/docs/cli/calculators.md
similarity index 86%
rename from docs/cli/interatomic.md
rename to docs/cli/calculators.md
index 8a96647a..a27d44fb 100644
--- a/docs/cli/interatomic.md
+++ b/docs/cli/calculators.md
@@ -4,7 +4,7 @@
 
 This command generates an averaged EAM potential for an alloy. The basic theory behind averaged alloy EAM potentials is
 described in
-[Varvenne, Luque,  Nöhring, Curtin, Phys. Rev. B 93, 104201 (2016)](https://doi.org/10.1103/PhysRevB.93.104201). The
+[Varvenne, Luque, Nöhring, Curtin, Phys. Rev. B 93, 104201 (2016)](https://doi.org/10.1103/PhysRevB.93.104201). The
 basic usage of the command is
 
 ```bash
diff --git a/docs/cli/fracture_mechanics.md b/docs/cli/fracture_mechanics.md
new file mode 100644
index 00000000..d1075eec
--- /dev/null
+++ b/docs/cli/fracture_mechanics.md
@@ -0,0 +1,67 @@
+# Fracture mechanics
+
+## `matscipy-quasistatic-crack`
+
+This command carries out a quasistatic K-field-controlled crack calculation on a small cluster. The outer boundaries of
+the cluster are fixed according to the near field solution of linear elastic fracture mechanics, considering elastic
+anisotropy, as described in [Sih, Paris, Irwin, Int. J. Fract. Mech. 1, 189 (1965)](https://doi.org/10.1007/BF00186854).
+The crack advances by stepwise updates of $K_\textrm{I}$ and and crack tip position.
+
+The command does not take any arguments but can be configured using a `params.py` file in the current directory.
+An example `params.py` file follows:
+```Python
+import numpy as np
+
+import ase.io
+
+from matscipy.fracture_mechanics.clusters import diamond, set_groups
+
+from atomistica import TersoffScr, Tersoff_PRB_39_5566_Si_C__Scr
+
+###
+# Interaction potential
+calc = TersoffScr(**Tersoff_PRB_39_5566_Si_C__Scr)
+
+# Fundamental material properties
+el = 'C'
+a0 = 3.57
+
+# Elastic constants: either specify explicitly or compute automatically
+compute_elastic_constants = True
+# C11 = 1220.   # GPa
+# C12 = -3.    # GPa
+# C44 = 535.    # GPa
+
+# Surface energy
+surface_energy = 2.7326 * 10  # GPa*A = 0.1 J/m^2
+
+# Crack system
+crack_surface = [1, 1, 1]  # Normal of crack face
+crack_front = [1, -1, 0]  # Direction of crack front
+# bond = (10, 11)
+bondlength = 1.7  # Bond length for crack tip detection
+bulk_nn = 4  # Number of nearest neighbors in the bulk
+
+vacuum = 6.0  # Vacuum surrounding the cracked cluster
+
+# Simulation control
+nsteps = 31
+# Increase stress intensity factor
+k1 = np.linspace(0.8, 1.2, nsteps)
+# Don't move crack tip
+tip_dx = np.zeros_like(k1)
+tip_dz = np.zeros_like(k1)
+
+fmax = 0.05  # Tolerance for quasistatic optimizer
+
+# The actual crack system
+n = [1, 1, 1]
+skin_x, skin_y = 1, 1
+cryst = diamond(el, a0, n, crack_surface, crack_front)
+set_groups(cryst, n, skin_x, skin_y)  # Outer fixed atoms
+ase.io.write('cryst.xyz', cryst, format='extxyz')  # Dump initial crack system (without notch)
+```
+
+## `matscipy-sinclair-continuation`
+
+## `matscipy-sinclair-crack`
\ No newline at end of file
diff --git a/docs/cli/index.rst b/docs/cli/index.rst
index 08065d71..7b37e5c4 100644
--- a/docs/cli/index.rst
+++ b/docs/cli/index.rst
@@ -14,5 +14,6 @@ whether a default parameter was used, to the log.
 
    diffusion
    electrochemistry
-   interatomic
+   fracture_mechanics
+   calculators
    structure
\ No newline at end of file
diff --git a/matscipy/cli/fracture_mechanics/__init__.py b/matscipy/cli/fracture_mechanics/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/matscipy/tool/meson.build b/matscipy/cli/fracture_mechanics/meson.build
similarity index 53%
rename from matscipy/tool/meson.build
rename to matscipy/cli/fracture_mechanics/meson.build
index 0efaf425..5a8491e0 100644
--- a/matscipy/tool/meson.build
+++ b/matscipy/cli/fracture_mechanics/meson.build
@@ -3,11 +3,14 @@
 # Pure Python sources
 python_sources = [
     '__init__.py',
-    'plot.py'
+    'quasistatic_crack.py',
+    'setup_crack.py',
+    'sinclair_continuation.py',
+    'sinclair_crack.py'
 ]
 
 # Install pure Python
 python.install_sources(
     python_sources,
-    subdir: 'matscipy/tool'
+    subdir: 'matscipy/cli/fracture_mechanics'
 )
diff --git a/matscipy/cli/fracture_mechanics/quasistatic_crack.py b/matscipy/cli/fracture_mechanics/quasistatic_crack.py
new file mode 100644
index 00000000..965918b5
--- /dev/null
+++ b/matscipy/cli/fracture_mechanics/quasistatic_crack.py
@@ -0,0 +1,177 @@
+#
+# Copyright 2014-2015, 2021 Lars Pastewka (U. Freiburg)
+#           2015 James Kermode (Warwick U.)
+#
+# matscipy - Materials science with Python at the atomic-scale
+# https://github.com/libAtoms/matscipy
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+
+# USAGE:
+#
+# Code imports the file 'params.py' from current working directory. params.py
+# contains simulation parameters. Some parameters can be omitted, see below.
+#
+# Parameters
+# ----------
+# calc : ase.Calculator
+#     Calculator object for energy and force computation.
+# tip_x0 : float
+#     Initial x-position of crack tip.
+# tip_y0 : float
+#     Initial y-position of crack tip.
+# tip_dx : array-like
+#     Displacement of tip in x-direction during run. x- and y-positions will be
+#     optimized self-consistently if omitted.
+# tip_dy : array-like
+#     Displacement of tip in y-direction during run. Position will be optimized
+#     self-consistently if omitted.
+
+import sys
+
+import numpy as np
+
+import ase.io
+import ase.constraints
+import ase.optimize
+
+from matscipy import parameter
+from matscipy.logger import screen
+
+from setup_crack import setup_crack
+from ase.optimize.precon import PreconLBFGS
+
+###
+
+def main():
+    # Atom types used for outputting the crack tip position.
+    ACTUAL_CRACK_TIP = 'Au'
+    FITTED_CRACK_TIP = 'Ag'
+
+    ###
+
+    logger = screen
+
+    ###
+
+    a, cryst, crk, k1g, tip_x0, tip_y0, bond1, bond2, boundary_mask, \
+        boundary_mask_bulk, tip_mask = setup_crack(logger=logger)
+    ase.io.write('notch.xyz', a, format='extxyz')
+
+    # Global parameters
+    basename = parameter('basename', 'quasistatic_crack')
+    calc = parameter('calc')
+    fmax = parameter('fmax', 0.01)
+
+    # Determine simulation control
+    k1_list = parameter('k1')
+    old_k1 = k1_list[0]
+    nsteps = len(k1_list)
+    tip_dx_list = parameter('tip_dx', np.zeros(nsteps))
+    tip_dy_list = parameter('tip_dy', np.zeros(nsteps))
+
+    # Run crack calculation.
+    tip_x = tip_x0
+    tip_y = tip_y0
+    a.set_calculator(calc)
+    for i, (k1, tip_dx, tip_dy) in enumerate(zip(k1_list, tip_dx_list,
+                                                 tip_dy_list)):
+        logger.pr('=== k1 = {0}*k1g, tip_dx = {1}, tip_dy = {2} ===' \
+                  .format(k1, tip_dx, tip_dy))
+        if tip_dx is None or tip_dy is None:
+            #
+            # Optimize crack tip position
+            #
+            old_y = tip_y + 1.0
+            old_x = tip_x + 1.0
+            while abs(tip_x - old_x) > 1e-6 and abs(tip_y - old_y) > 1e-6:
+                b = cryst.copy()
+                ux, uy = crk.displacements(cryst.positions[:, 0], cryst.positions[:, 1],
+                                           tip_x, tip_y, k * k1g)
+                b.positions[:, 0] += ux
+                b.positions[:, 1] += uy
+                a.set_constraint(None)
+                a.positions[boundary_mask] = b.positions[boundary_mask]
+                a.set_constraint(ase.constraints.FixAtoms(mask=boundary_mask))
+                logger.pr('Optimizing positions...')
+                opt = ase.optimize.FIRE(a, logfile=None)
+                opt.run(fmax=fmax)
+                logger.pr('...done. Converged within {0} steps.' \
+                          .format(opt.get_number_of_steps()))
+
+                old_x = tip_x
+                old_y = tip_y
+                tip_x, tip_y = crk.crack_tip_position(a.positions[:, 0],
+                                                      a.positions[:, 1],
+                                                      cryst.positions[:, 0],
+                                                      cryst.positions[:, 1],
+                                                      tip_x, tip_y, k * k1g,
+                                                      mask=mask)
+        else:
+            #
+            # Do not optimize tip position.
+            #
+            tip_x = tip_x0 + tip_dx
+            tip_y = tip_y0 + tip_dy
+            logger.pr('Setting crack tip position to {0} {1}' \
+                      .format(tip_x, tip_y))
+            # Scale strain field and optimize crack
+            a.set_constraint(None)
+            x, y = crk.scale_displacements(a.positions[:len(cryst), 0],
+                                           a.positions[:len(cryst), 1],
+                                           cryst.positions[:, 0],
+                                           cryst.positions[:, 1],
+                                           old_k1, k1)
+            a.positions[:len(cryst), 0] = x
+            a.positions[:len(cryst), 1] = y
+            # Optimize atoms in center
+            a.set_constraint(ase.constraints.FixAtoms(mask=boundary_mask))
+            logger.pr('Optimizing positions...')
+            opt = PreconLBFGS(a)
+            opt.run(fmax=fmax)
+            logger.pr('...done. Converged within {0} steps.' \
+                      .format(opt.get_number_of_steps()))
+            old_k1 = k1
+
+        # Output optimized configuration ot file. Include the mask array in
+        # output, so we know which atoms were used in fitting.
+        a.set_array('atoms_used_for_fitting_crack_tip', tip_mask)
+
+        # Output a file that contains the target crack tip (used for
+        # the displacmenets of the boundary atoms) and the fitted crack tip
+        # positions. The target crack tip is marked by a Hydrogen atom.
+        b = a.copy()
+        b += ase.Atom(ACTUAL_CRACK_TIP, (tip_x, tip_y, b.cell[2, 2] / 2))
+
+        # Measure the true (fitted) crack tip position.
+        try:
+            measured_tip_x, measured_tip_y = \
+                crk.crack_tip_position(a.positions[:, 0], a.positions[:, 1],
+                                       cryst.positions[:, 0], cryst.positions[:, 1],
+                                       tip_x, tip_y, k1 * k1g, mask=tip_mask)
+            measured_tip_x %= a.cell[0][0]
+            measured_tip_y %= a.cell[0][0]
+        except:
+            measured_tip_x = 0.0
+            measured_tip_y = 0.0
+
+        # The fitted crack tip is marked by a Helium atom.
+        b += ase.Atom(FITTED_CRACK_TIP, (measured_tip_x, measured_tip_y,
+                                         b.cell[2, 2] / 2))
+
+        b.info['bond_length'] = a.get_distance(bond1, bond2)
+        b.info['energy'] = a.get_potential_energy()
+        b.info['cell_origin'] = [0, 0, 0]
+        ase.io.write('%s_%4.4i.xyz' % (basename, i), b, format='extxyz')
diff --git a/scripts/fracture_mechanics/setup_crack.py b/matscipy/cli/fracture_mechanics/setup_crack.py
similarity index 83%
rename from scripts/fracture_mechanics/setup_crack.py
rename to matscipy/cli/fracture_mechanics/setup_crack.py
index 9525e62b..07661869 100644
--- a/scripts/fracture_mechanics/setup_crack.py
+++ b/matscipy/cli/fracture_mechanics/setup_crack.py
@@ -18,6 +18,7 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
+
 import numpy as np
 
 import ase.optimize
@@ -30,6 +31,7 @@
 from matscipy.logger import screen
 from matscipy.neighbours import neighbour_list
 
+
 ###
 
 def setup_crack(logger=screen):
@@ -56,11 +58,11 @@ def setup_crack(logger=screen):
                                          fmax=elastic_fmax)
         log_file.close()
         logger.pr('Measured elastic constants (in GPa):')
-        logger.pr(np.round(C*10/GPa)/10)
+        logger.pr(np.round(C * 10 / GPa) / 10)
 
         crk = crack.CubicCrystalCrack(parameter('crack_surface'),
                                       parameter('crack_front'),
-                                      Crot=C/GPa)
+                                      Crot=C / GPa)
     else:
         if has_parameter('C'):
             crk = crack.CubicCrystalCrack(parameter('crack_surface'),
@@ -72,17 +74,16 @@ def setup_crack(logger=screen):
                                           parameter('C11'), parameter('C12'),
                                           parameter('C44'))
 
-
     logger.pr('Elastic constants used for boundary condition (in GPa):')
-    logger.pr(np.round(crk.C*10)/10)
+    logger.pr(np.round(crk.C * 10) / 10)
 
     # Get Griffith's k1.
     k1g = crk.k1g(parameter('surface_energy'))
     logger.pr('Griffith k1 = %f' % k1g)
 
     # Apply initial strain field.
-    tip_x = parameter('tip_x', cryst.cell.diagonal()[0]/2)
-    tip_y = parameter('tip_y', cryst.cell.diagonal()[1]/2)
+    tip_x = parameter('tip_x', cryst.cell.diagonal()[0] / 2)
+    tip_y = parameter('tip_y', cryst.cell.diagonal()[1] / 2)
 
     bondlength = parameter('bondlength', 2.7)
     bulk_nn = parameter('bulk_nn', 4)
@@ -97,19 +98,19 @@ def setup_crack(logger=screen):
         # Get surface atoms of cluster with crack
         a = hydrogenate(cryst, bondlength, parameter('XH_bondlength'), b=a)
         g = a.get_array('groups')
-        g[a.numbers==1] = -1
+        g[a.numbers == 1] = -1
         a.set_array('groups', g)
         cryst = a.copy()
 
     k1 = parameter('k1')
     try:
-      k1 = k1[0]
+        k1 = k1[0]
     except:
-      pass
-    ux, uy = crk.displacements(cryst.positions[:,0], cryst.positions[:,1],
-                               tip_x, tip_y, k1*k1g)
-    a.positions[:len(cryst),0] += ux
-    a.positions[:len(cryst),1] += uy
+        pass
+    ux, uy = crk.displacements(cryst.positions[:, 0], cryst.positions[:, 1],
+                               tip_x, tip_y, k1 * k1g)
+    a.positions[:len(cryst), 0] += ux
+    a.positions[:len(cryst), 1] += uy
 
     # Center notched configuration in simulation cell and ensure enough vacuum.
     oldr = a[0].position.copy()
@@ -124,7 +125,7 @@ def setup_crack(logger=screen):
         parameter('bond', crack.find_tip_coordination(a, bondlength=bondlength, bulk_nn=bulk_nn))
 
     if parameter('center_crack_tip_on_bond', False):
-        tip_x, tip_y, dummy = (a.positions[bond1]+a.positions[bond2])/2
+        tip_x, tip_y, dummy = (a.positions[bond1] + a.positions[bond2]) / 2
 
     # Hydrogenate?
     coord = np.bincount(neighbour_list('i', a, bondlength), minlength=len(a))
@@ -134,19 +135,19 @@ def setup_crack(logger=screen):
         g = a.get_array('groups')
         gcryst = cryst.get_array('groups')
         coord = a.get_array('coord')
-        g[coord!=4] = -1
-        gcryst[coord!=4] = -1
+        g[coord != 4] = -1
+        gcryst[coord != 4] = -1
         a.set_array('groups', g)
         cryst.set_array('groups', gcryst)
 
     if hydrogenate_flag and hydrogenate_crack_face_flag:
         # Get surface atoms of cluster with crack
-        exclude = np.logical_and(a.get_array('groups')==1, coord!=4)
+        exclude = np.logical_and(a.get_array('groups') == 1, coord != 4)
         a.set_array('exclude', exclude)
         a = hydrogenate(cryst, bondlength, parameter('XH_bondlength'), b=a,
                         exclude=exclude)
         g = a.get_array('groups')
-        g[a.numbers==1] = -1
+        g[a.numbers == 1] = -1
         a.set_array('groups', g)
         basename = parameter('basename', 'energy_barrier')
         ase.io.write('{0}_hydrogenated.xyz'.format(basename), a,
@@ -162,15 +163,15 @@ def setup_crack(logger=screen):
     gcryst = cryst.get_array('groups')
 
     logger.pr('Opening bond {0}--{1}, initial bond length {2}'.
-          format(bond1, bond2, a.get_distance(bond1, bond2, mic=True)))
+              format(bond1, bond2, a.get_distance(bond1, bond2, mic=True)))
 
     # centre vertically on the opening bond
     if parameter('center_cell_on_bond', True):
-      a.translate([0., a.cell[1,1]/2.0 -
-                      (a.positions[bond1, 1] +
-                       a.positions[bond2, 1])/2.0, 0.])
+        a.translate([0., a.cell[1, 1] / 2.0 -
+                     (a.positions[bond1, 1] +
+                      a.positions[bond2, 1]) / 2.0, 0.])
 
     a.info['bond1'] = bond1
     a.info['bond2'] = bond2
 
-    return a, cryst, crk, k1g, tip_x, tip_y, bond1, bond2, g==0, gcryst==0, g==1
+    return a, cryst, crk, k1g, tip_x, tip_y, bond1, bond2, g == 0, gcryst == 0, g == 1
diff --git a/matscipy/cli/fracture_mechanics/sinclair_continuation.py b/matscipy/cli/fracture_mechanics/sinclair_continuation.py
new file mode 100644
index 00000000..fbb61e40
--- /dev/null
+++ b/matscipy/cli/fracture_mechanics/sinclair_continuation.py
@@ -0,0 +1,246 @@
+#
+# Copyright 2020-2021 James Kermode (Warwick U.)
+#
+# matscipy - Materials science with Python at the atomic-scale
+# https://github.com/libAtoms/matscipy
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+
+import matscipy
+
+import os
+import numpy as np
+
+import h5py
+import ase.io
+from ase.units import GPa
+from ase.constraints import FixAtoms
+# from ase.optimize import LBFGSLineSearch
+from ase.optimize.precon import PreconLBFGS
+
+from matscipy import parameter
+from matscipy.elasticity import fit_elastic_constants
+from matscipy.fracture_mechanics.crack import CubicCrystalCrack, SinclairCrack
+
+from scipy.optimize import fsolve, fminbound
+
+import sys
+
+def main():
+    print(matscipy.__file__)
+
+    sys.path.insert(0, '../../../scripts/fracture_mechanics')
+
+    import params
+
+    calc = parameter('calc')
+    fmax = parameter('fmax', 1e-3)
+    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)
+    ds = parameter('ds', 1e-2)
+    nsteps = parameter('nsteps', 10)
+    a0 = parameter('a0')  # lattice constant
+    k0 = parameter('k0', 1.0)
+    extended_far_field = parameter('extended_far_field', False)
+    alpha0 = parameter('alpha0', 0.0)  # initial guess for crack position
+    dump = parameter('dump', False)
+    precon = parameter('precon', False)
+    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')
+
+    # 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'))
+
+    # setup the crack
+    crk = CubicCrystalCrack(parameter('crack_surface'),
+                            parameter('crack_front'),
+                            Crot=C / GPa)
+
+    # get Griffith's stess intensity factor
+    k1g = crk.k1g(parameter('surface_energy'))
+    print('Griffith k1 = %f' % k1g)
+
+    # check for restart files if continuation=True
+    if continuation and not os.path.exists(restart_file):
+        continuation = False
+    if continuation:
+
+        # 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;
+            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:
+
+        # 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)
+
+        traj = None
+        dk = parameter('dk', 1e-4)
+        dalpha = parameter('dalpha', 1e-3)
+
+        # reuse output from sinclair_crack.py if possible
+        if os.path.exists(f'k_{int(k0 * 1000):04d}.xyz'):
+            print(f'Reading atoms from k_{int(k0 * 1000):04d}.xyz')
+            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
+
+        # 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
+        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
+            sc.variable_alpha = False
+            sc.variable_k = False
+            sc.update_atoms()
+            atoms = sc.atoms.copy()
+            atoms.calc = sc.calc
+            atoms.set_constraint(FixAtoms(mask=~sc.regionI))
+            opt = PreconLBFGS(atoms, logfile=None)
+            # 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}')
+            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()
+
+        # 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_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
+
+        # 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]
+        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
+
+    # 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,
+                        extended_far_field=extended_far_field)
+
+    # run full arc length continuation
+    scv.arc_length_continuation(x0, x1, N=nsteps,
+                                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
diff --git a/matscipy/cli/fracture_mechanics/sinclair_crack.py b/matscipy/cli/fracture_mechanics/sinclair_crack.py
new file mode 100644
index 00000000..c89155b2
--- /dev/null
+++ b/matscipy/cli/fracture_mechanics/sinclair_crack.py
@@ -0,0 +1,149 @@
+#
+# Copyright 2014, 2020 James Kermode (Warwick U.)
+#           2020 Arnaud Allera (U. Lyon 1)
+#           2014 Lars Pastewka (U. Freiburg)
+#
+# matscipy - Materials science with Python at the atomic-scale
+# https://github.com/libAtoms/matscipy
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+
+import os
+import numpy as np
+
+import ase.io
+from ase.units import GPa
+from ase.constraints import FixAtoms
+from ase.optimize.precon import PreconLBFGS
+
+from matscipy import parameter
+from matscipy.elasticity import fit_elastic_constants
+from matscipy.fracture_mechanics.crack import CubicCrystalCrack, SinclairCrack
+
+import sys
+
+def main():
+    sys.path.insert(0, '../../../scripts/fracture_mechanics')
+    import params
+
+    calc = parameter('calc')
+    fmax = parameter('fmax', 1e-3)
+    vacuum = parameter('vacuum', 10.0)
+    flexible = parameter('flexible', True)
+    extended_far_field = parameter('extended_far_field', False)
+    k0 = parameter('k0', 1.0)
+    alpha0 = parameter('alpha0', 0.0)  # initial guess for crack position
+    dump = parameter('dump', False)
+    precon = parameter('precon', False)
+    prerelax = parameter('prerelax', False)
+    lbfgs = parameter('lbfgs', not flexible)  # use LBGS by default if not flexible
+
+    # compute elastic constants
+    cryst = params.cryst.copy()
+    cryst.pbc = True
+    cryst.calc = calc
+    C, C_err = fit_elastic_constants(cryst,
+                                     symmetry=parameter('elastic_symmetry',
+                                                        'triclinic'))
+
+    crk = CubicCrystalCrack(parameter('crack_surface'),
+                            parameter('crack_front'),
+                            Crot=C / GPa)
+
+    # Get Griffith's k1.
+    k1g = crk.k1g(parameter('surface_energy'))
+    print('Griffith k1 = %f' % k1g)
+
+    cluster = params.cluster.copy()
+    sc = SinclairCrack(crk, cluster, calc, k0 * k1g,
+                       alpha=alpha0,
+                       variable_alpha=flexible,
+                       vacuum=vacuum,
+                       extended_far_field=extended_far_field)
+
+    nsteps = parameter('nsteps')
+    k1_range = parameter('k1_range', np.linspace(0.8, 1.2, nsteps))
+
+    ks = list(k1_range)
+    ks_out = []
+    alphas = []
+
+    max_steps = parameter('max_steps', 10)
+    fit_alpha = parameter('fit_alpha', False)
+
+    for i, k in enumerate(ks):
+        restart_file = parameter('restart_file', f'k_{int(k * 1000):04d}.xyz')
+        if os.path.exists(restart_file):
+            print(f'Restarting from {restart_file}')
+            if restart_file.endswith('.xyz'):
+                a = ase.io.read(restart_file)
+                sc.set_atoms(a)
+            elif restart_file.endswith('.txt'):
+                x = np.loadtxt(restart_file)
+                if len(x) == len(sc):
+                    sc.set_dofs(x)
+                elif len(x) == len(sc) + 1:
+                    sc.variable_k = True
+                    sc.set_dofs(x)
+                    sc.variable_k = False
+                elif len(x) == len(sc) + 2:
+                    sc.variable_alpha = True
+                    sc.variable_k = True
+                    sc.set_dofs(x)
+                    sc.variable_alpha = False
+                    sc.variable_k = False
+                else:
+                    raise RuntimeError('cannot guess how to restart with'
+                                       f' {len(x)} variables for {len(self)} DoFs')
+            else:
+                raise ValueError(f"don't know how to restart from {restart_file}")
+        else:
+            sc.rescale_k(k * k1g)
+
+        if fit_alpha:
+            sc.alpha, = sc.fit_cle(variable_alpha=True, variable_k=False)
+            print(f'Fitted value of alpha: {sc.alpha}')
+        print(f'k = {sc.k / k1g} * k1g')
+        print(f'alpha = {sc.alpha}')
+
+        if lbfgs:
+            print('Optimizing with LBFGS')
+            atoms = sc.atoms.copy()
+            atoms.calc = sc.calc
+            atoms.set_constraint(FixAtoms(mask=np.logical_not(sc.regionI)))
+            opt = PreconLBFGS(atoms)
+            opt.run(fmax=fmax)
+            sc.set_atoms(atoms)
+        else:
+            if prerelax:
+                print('Pre-relaxing with Conjugate-Gradients')
+                sc.optimize(ftol=1e-5, steps=max_steps, dump=dump,
+                            method='cg')
+
+            sc.optimize(fmax, steps=max_steps, dump=dump,
+                        precon=precon, method='krylov')
+
+        if flexible:
+            print(f'Optimized alpha = {sc.alpha:.3f}')
+        ks_out.append(k)
+        alphas.append(sc.alpha)
+
+        a = sc.atoms
+        a.get_forces()
+        ase.io.write(f'k_{int(k * 1000):04d}.xyz', a)
+        with open('x.txt', 'a') as fh:
+            np.savetxt(fh, [sc.get_dofs()])
+
+    np.savetxt('k_vs_alpha.txt', np.c_[ks_out, alphas])
diff --git a/matscipy/tool/__init__.py b/matscipy/tool/__init__.py
deleted file mode 100644
index 62485b6d..00000000
--- a/matscipy/tool/__init__.py
+++ /dev/null
@@ -1,24 +0,0 @@
-#
-# Copyright 2014-2015, 2021 Lars Pastewka (U. Freiburg)
-#           2015-2016 Adrien Gola (KIT)
-#           2014 James Kermode (Warwick U.)
-#
-# matscipy - Materials science with Python at the atomic-scale
-# https://github.com/libAtoms/matscipy
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 2 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
-
-# Generic stuff may go here.
-
diff --git a/matscipy/tool/plot.py b/matscipy/tool/plot.py
deleted file mode 100644
index b5181e6c..00000000
--- a/matscipy/tool/plot.py
+++ /dev/null
@@ -1,201 +0,0 @@
-#
-# Copyright 2016, 2021 Lars Pastewka (U. Freiburg)
-#           2016 Adrien Gola (KIT)
-#
-# matscipy - Materials science with Python at the atomic-scale
-# https://github.com/libAtoms/matscipy
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 2 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
-
-
-"""
- Usage : python plot.py [filename]
- 
- This tool script allow to quicly plot data from text files arranged in column with space as separator. It uses a GUI through Tkinter
- "log.lammp" can be read and splited according to different "run" present in the log file
-                                                                                            
- Adrien Gola // 4.11.2015                                                                        
-"""
-
-import numpy as np
-import matplotlib
-matplotlib.rcParams['backend'] = "Qt4Agg"
-import matplotlib.pyplot as p 
-
-import os,sys
-import os.path
-
-import Tkinter as tk
-import ttk
-
-# ------------------------------------
-# ------- Function definition --------
-# ------------------------------------
-def close_windows():
-    root.destroy()
-
-def plot_button():
-    global xlabel,ylabel
-    x = d[:,Xvar.get()]
-    y = d[:,Yvar.get()]
-  
-    if xlabel.get() != "X-axis label":
-        hx = xlabel.get()
-    else:
-        hx = h[Xvar.get()]
-    if ylabel.get() != "Y-axis label":
-        hy = ylabel.get()
-    else:
-        hy = h[Yvar.get()]
-  
-    p.plot(x, y, 'o')
-    p.xlabel(hx)
-    p.ylabel(hy)
-    p.show() 
-  
-def process_file():
-    global variables_frame,h,d,Xvar,Yvar,skipline,flag_loglammps
-    selected_file = sorted(outlogfiles)[Fvar.get()]
-    if flag_loglammps:
-        h = open(selected_file,'r').readlines()[0].strip('#').strip().split(" ") # headers list 
-        d = np.loadtxt(selected_file, skiprows=1) # data  
-    else:
-        h = filter(None,open(selected_file,'r').readlines()[0+int(skipline.get())].strip('#').strip().split(" ")) # headers list 
-        d = np.loadtxt(selected_file, skiprows=1+int(skipline.get())) # data   
-    try:
-        for child in variables_frame.winfo_children():
-            child.destroy()
-    except:
-        pass
-    #Xvar
-    tk.Label(variables_frame,text="Select X-axis variable to plot",anchor='n').grid(row=next_row+1)
-    Xvar = tk.IntVar()
-    for i,x in enumerate(h):
-        x_row=i+2
-        tk.Radiobutton(variables_frame,text=x,variable=Xvar, value=i).grid(row=next_row+x_row)
-    #Yvar
-    tk.Label(variables_frame,text="Select Y-axis variable to plot",anchor='n').grid(row=next_row+1,column=1)
-    Yvar = tk.IntVar()
-    for j,x in enumerate(h):
-        y_row=j+2
-        tk.Radiobutton(variables_frame,text=x,variable=Yvar, value=j).grid(row=next_row+y_row,column=1)
-
-
-mypath = os.getcwd()
-# -------------------------------------------------------
-# ------- Reading and spliting the main log file --------
-# -------------------------------------------------------
-flag_loglammps=0
-# --- file selection ---
-try :
-    in_file = sys.argv[1]
-except:
-    print("Usage : plot_general.py file_name")
-    quit()
-
-if in_file == "log.lammps":
-    flag_loglammps=1
-    data = open("log.lammps",'r').readlines()  # select log.lammps as main log file if it exists
-      
-    # --- spliting into "out.log.*" files ---
-    j = 0
-    flag = 0
-    flag_out = 0
-    step = ""
-    for lines in data:
-        lines = lines.strip()
-        if lines.startswith("run:"):
-            step = "."+lines.strip().strip("run: ")
-        elif lines.startswith("Memory usage per processor"):
-          if len(str(j)) == 1:
-              J = "0"+str(j)
-          else:
-              J = str(j)
-          out = open("out.log.%s%s-plt-tmp"%(J,step),'w')
-          flag_out = 1
-          flag = 1
-        elif flag and lines.startswith("Loop") or lines.startswith("kill"):
-            flag = 0
-            j+=1
-            out.close()
-            step = ""
-        elif flag:
-            out.write(lines+'\n')
-    if flag_out:
-        out.close()
-    outlogfiles = [ f for f in os.listdir(mypath) if os.path.isfile(os.path.join(mypath,f)) and f.startswith("out.log") and f.endswith("plt-tmp")] # list of "out.log.*" files
-else:
-    outlogfiles = [in_file]
-    pass
-  
-
-
-
-
-###################
-### MAIN script ###
-###################
-
-root=tk.Tk()
-# Positioning frames
-button_frame = tk.Frame(root,heigh=5)
-button_frame.pack(side='right')
-
-source_frame = tk.Frame(root)
-source_frame.pack()
-
-labels_frame = tk.Frame(root)
-labels_frame.pack(side="bottom")
-
-variables_frame = tk.Frame(root)
-variables_frame.pack(side="bottom")
-
-# Filling frames
-tk.Label(source_frame,text="Select data source",anchor='n').grid(row=0)
-Fvar = tk.IntVar()
-for i,x in enumerate(outlogfiles):
-    next_row=i+1
-    tk.Radiobutton(source_frame,text=x,variable=Fvar, value=i, command = process_file).grid(row=next_row)
-    if not flag_loglammps:
-        skipline = tk.Entry(source_frame,width=5,justify="center")
-        skipline.grid(row=next_row,column=1)
-        skipline.insert(0,0)
-    else:
-        skipline=0
-next_row+=1  
-ttk.Separator(source_frame,orient='horizontal').grid(row=next_row,columnspan=2,sticky="ew") 
-
-# Insert X,Y label text field
-ttk.Separator(labels_frame,orient='horizontal').grid(columnspan=2,sticky="ew") 
-xlabel = tk.Entry(labels_frame,width=50,justify="center")
-xlabel.grid(row=next_row+1,column=0,columnspan=2)
-xlabel.insert(0,"X-axis label")
-ylabel = tk.Entry(labels_frame,width=50,justify="center")
-ylabel.grid(row=next_row+2,column=0,columnspan=2)
-ylabel.insert(0,"Y-axis label")
-
-close = tk.Button(button_frame, text = "Close",width=10, command = close_windows)
-close.grid(row=1)
-plot = tk.Button(button_frame, text = "Plot",width=10, command = plot_button)
-plot.grid(row=0)
-
-# show GUI
-root.mainloop()
-
-
-# Cleaning tmp files
-if flag_loglammps:
-    for f in outlogfiles:
-        os.remove(f)
diff --git a/pyproject.toml b/pyproject.toml
index ab976ec9..88b69b46 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -70,7 +70,7 @@ cli = [
 # Generic
 matscipy-rms = "matscipy.cli.diffusion.rms:main"
 
-# Glasses
+# Structure generation
 matscipy-quench = "matscipy.cli.glasses.quench:main"
 
 # Electrochemistry
@@ -78,6 +78,11 @@ matscipy-continuous2discrete = "matscipy.cli.electrochemistry.continuous2discret
 matscipy-poisson-nernst-planck = "matscipy.cli.electrochemistry.poisson_nernst_planck_solver:main"
 matscipy-stericify = "matscipy.cli.electrochemistry.stericify:main"
 
+# Fracture mechanics
+matscipy-quasistatic-crack = "matscipy.cli.fracture_mechanics.quasistatic_crack:main"
+matscipy-sinclair-continuation = "matscipy.cli.fracture_mechanics.sinclair_continuation:main"
+matscipy-sinclair-crack = "matscipy.cli.fracture_mechanics.sinclair_crack:main"
+
 # Interatomic potentials
 matscipy-average-eam-potential = "matscipy.cli.calculators.average_eam_potential:main"
 
diff --git a/scripts/fracture_mechanics/energy_barrier.py b/scripts/fracture_mechanics/energy_barrier.py
index 992ece83..f18abfbb 100644
--- a/scripts/fracture_mechanics/energy_barrier.py
+++ b/scripts/fracture_mechanics/energy_barrier.py
@@ -46,26 +46,21 @@
 # Code imports the file 'params.py' from current working directory. params.py
 # contains simulation parameters. Some parameters can be omitted, see below.
 
-from math import sqrt
-
 import os
-import sys
 
 import numpy as np
 
-import ase
 from ase.constraints import FixConstraint
 import ase.io
 import ase.optimize
 from ase.data import atomic_numbers
-from ase.units import GPa
 from ase.geometry import find_mic
 
 import matscipy.fracture_mechanics.crack as crack
 from matscipy import parameter
 from matscipy.logger import screen
 
-from setup_crack import setup_crack
+from matscipy.cli.fracture_mechanics.setup_crack import setup_crack
 
 ###
 
diff --git a/scripts/fracture_mechanics/neb_crack_tip.py b/scripts/fracture_mechanics/neb_crack_tip.py
index 69b8048f..5ce3a63f 100644
--- a/scripts/fracture_mechanics/neb_crack_tip.py
+++ b/scripts/fracture_mechanics/neb_crack_tip.py
@@ -46,11 +46,9 @@
 # contains simulation parameters. Some parameters can be omitted, see below.
 
 import os
-import sys
 
 import numpy as np
 
-import ase
 import ase.constraints
 import ase.io
 #import ase.optimize.precon
@@ -58,14 +56,13 @@
 from ase.neb import NEB
 
 from ase.data import atomic_numbers
-from ase.units import GPa
 
 import matscipy.fracture_mechanics.crack as crack
 from matscipy import parameter
 from matscipy.logger import screen
 from matscipy.io import savetbl
 
-from setup_crack import setup_crack
+from matscipy.cli.fracture_mechanics.setup_crack import setup_crack
 
 from scipy.integrate import cumtrapz
 
diff --git a/scripts/fracture_mechanics/optimize_crack_tip.py b/scripts/fracture_mechanics/optimize_crack_tip.py
index d4ca8069..93198d62 100644
--- a/scripts/fracture_mechanics/optimize_crack_tip.py
+++ b/scripts/fracture_mechanics/optimize_crack_tip.py
@@ -45,23 +45,18 @@
 # Code imports the file 'params.py' from current working directory. params.py
 # contains simulation parameters. Some parameters can be omitted, see below.
 
-import os
-import sys
-
 import numpy as np
 
-import ase
 import ase.constraints
 import ase.io
 import ase.optimize
 from ase.data import atomic_numbers
-from ase.units import GPa
 
 import matscipy.fracture_mechanics.crack as crack
 from matscipy import parameter
 from matscipy.logger import screen
 
-from setup_crack import setup_crack
+from matscipy.cli.fracture_mechanics.setup_crack import setup_crack
 
 ###
 
diff --git a/scripts/fracture_mechanics/quasistatic_crack.py b/scripts/fracture_mechanics/quasistatic_crack.py
deleted file mode 100644
index 8730f1dd..00000000
--- a/scripts/fracture_mechanics/quasistatic_crack.py
+++ /dev/null
@@ -1,206 +0,0 @@
-#
-# Copyright 2014-2015, 2021 Lars Pastewka (U. Freiburg)
-#           2015 James Kermode (Warwick U.)
-#
-# matscipy - Materials science with Python at the atomic-scale
-# https://github.com/libAtoms/matscipy
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 2 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
-
-# ======================================================================
-# matscipy - Python materials science tools
-# https://github.com/libAtoms/matscipy
-#
-# Copyright (2014) James Kermode, King's College London
-#                  Lars Pastewka, Karlsruhe Institute of Technology
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 2 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-# ======================================================================
-
-# USAGE:
-#
-# Code imports the file 'params.py' from current working directory. params.py
-# contains simulation parameters. Some parameters can be omitted, see below.
-#
-# Parameters
-# ----------
-# calc : ase.Calculator
-#     Calculator object for energy and force computation.
-# tip_x0 : float
-#     Initial x-position of crack tip.
-# tip_y0 : float
-#     Initial y-position of crack tip.
-# tip_dx : array-like
-#     Displacement of tip in x-direction during run. x- and y-positions will be
-#     optimized self-consistently if omitted.
-# tip_dy : array-like
-#     Displacement of tip in y-direction during run. Position will be optimized
-#     self-consistently if omitted.
-
-import os
-import sys
-
-import numpy as np
-
-import ase
-import ase.io
-import ase.constraints
-import ase.optimize
-from ase.units import GPa
-
-import matscipy.fracture_mechanics.crack as crack
-from matscipy import parameter
-from matscipy.logger import screen
-
-from setup_crack import setup_crack
-from ase.optimize.precon import PreconLBFGS
-
-###
-
-sys.path += ['.', '..']
-import params
-
-###
-
-# Atom types used for outputting the crack tip position.
-ACTUAL_CRACK_TIP = 'Au'
-FITTED_CRACK_TIP = 'Ag'
-
-###
-
-logger = screen
-
-###
-
-a, cryst, crk, k1g, tip_x0, tip_y0, bond1, bond2, boundary_mask, \
-    boundary_mask_bulk, tip_mask = setup_crack(logger=logger)
-ase.io.write('notch.xyz', a, format='extxyz')   
-
-# Global parameters
-basename = parameter('basename', 'quasistatic_crack')
-calc = parameter('calc')
-fmax = parameter('fmax', 0.01)
-
-# Determine simulation control
-k1_list = parameter('k1')
-old_k1 = k1_list[0]
-nsteps = len(k1_list)
-tip_dx_list = parameter('tip_dx', np.zeros(nsteps))
-tip_dy_list = parameter('tip_dy', np.zeros(nsteps))
-
-# Run crack calculation.
-tip_x = tip_x0
-tip_y = tip_y0
-a.set_calculator(calc)
-for i, ( k1, tip_dx, tip_dy ) in enumerate(zip(k1_list, tip_dx_list,
-                                               tip_dy_list)):
-    logger.pr('=== k1 = {0}*k1g, tip_dx = {1}, tip_dy = {2} ===' \
-              .format(k1, tip_dx, tip_dy))
-    if tip_dx is None or tip_dy is None:
-        #
-        # Optimize crack tip position
-        #
-        old_y = tip_y+1.0
-        old_x = tip_x+1.0
-        while abs(tip_x-old_x) > 1e-6 and abs(tip_y-old_y) > 1e-6:
-            b = cryst.copy()
-            ux, uy = crk.displacements(cryst.positions[:,0], cryst.positions[:,1],
-                                       tip_x, tip_y, k*k1g)
-            b.positions[:,0] += ux
-            b.positions[:,1] += uy
-            a.set_constraint(None)
-            a.positions[boundary_mask] = b.positions[boundary_mask]
-            a.set_constraint(ase.constraints.FixAtoms(mask=boundary_mask))
-            logger.pr('Optimizing positions...')
-            opt = ase.optimize.FIRE(a, logfile=None)
-            opt.run(fmax=fmax)
-            logger.pr('...done. Converged within {0} steps.' \
-                      .format(opt.get_number_of_steps()))
-        
-            old_x = tip_x
-            old_y = tip_y
-            tip_x, tip_y = crk.crack_tip_position(a.positions[:,0],
-                                                  a.positions[:,1],
-                                                  cryst.positions[:,0],
-                                                  cryst.positions[:,1],
-                                                  tip_x, tip_y, k*k1g,
-                                                  mask=mask)
-    else:
-        #
-        # Do not optimize tip position.
-        #
-        tip_x = tip_x0+tip_dx
-        tip_y = tip_y0+tip_dy
-        logger.pr('Setting crack tip position to {0} {1}' \
-                  .format(tip_x, tip_y))
-        # Scale strain field and optimize crack
-        a.set_constraint(None)
-        x, y = crk.scale_displacements(a.positions[:len(cryst),0],
-                                       a.positions[:len(cryst),1],
-                                       cryst.positions[:,0],
-                                       cryst.positions[:,1],
-                                       old_k1, k1)
-        a.positions[:len(cryst),0] = x
-        a.positions[:len(cryst),1] = y
-        # Optimize atoms in center
-        a.set_constraint(ase.constraints.FixAtoms(mask=boundary_mask))
-        logger.pr('Optimizing positions...')
-        opt = PreconLBFGS(a)
-        opt.run(fmax=fmax)
-        logger.pr('...done. Converged within {0} steps.' \
-                  .format(opt.get_number_of_steps()))
-        old_k1 = k1
-
-    # Output optimized configuration ot file. Include the mask array in
-    # output, so we know which atoms were used in fitting.
-    a.set_array('atoms_used_for_fitting_crack_tip', tip_mask)
-
-    # Output a file that contains the target crack tip (used for
-    # the displacmenets of the boundary atoms) and the fitted crack tip
-    # positions. The target crack tip is marked by a Hydrogen atom.
-    b = a.copy()
-    b += ase.Atom(ACTUAL_CRACK_TIP, (tip_x, tip_y, b.cell[2, 2]/2))
-
-    # Measure the true (fitted) crack tip position.
-    try:
-        measured_tip_x, measured_tip_y = \
-            crk.crack_tip_position(a.positions[:,0], a.positions[:,1],
-                                   cryst.positions[:,0], cryst.positions[:,1],
-                                   tip_x, tip_y, k1*k1g, mask=tip_mask)
-        measured_tip_x %= a.cell[0][0]
-        measured_tip_y %= a.cell[0][0]
-    except:
-        measured_tip_x = 0.0
-        measured_tip_y = 0.0
-
-    # The fitted crack tip is marked by a Helium atom.
-    b += ase.Atom(FITTED_CRACK_TIP, (measured_tip_x, measured_tip_y,
-                                     b.cell[2, 2]/2))
-
-    b.info['bond_length'] = a.get_distance(bond1, bond2)
-    b.info['energy'] = a.get_potential_energy()
-    b.info['cell_origin'] = [0, 0, 0]
-    ase.io.write('%s_%4.4i.xyz' % (basename, i), b, format='extxyz')
\ No newline at end of file
diff --git a/scripts/fracture_mechanics/sinclair_continuation.py b/scripts/fracture_mechanics/sinclair_continuation.py
deleted file mode 100644
index 081d0da4..00000000
--- a/scripts/fracture_mechanics/sinclair_continuation.py
+++ /dev/null
@@ -1,236 +0,0 @@
-#
-# Copyright 2020-2021 James Kermode (Warwick U.)
-#
-# matscipy - Materials science with Python at the atomic-scale
-# https://github.com/libAtoms/matscipy
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 2 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
-
-import matscipy; print(matscipy.__file__)
-
-import os
-import numpy as np
-
-import h5py
-import ase.io
-from ase.units import GPa
-from ase.constraints import FixAtoms
-#from ase.optimize import LBFGSLineSearch
-from ase.optimize.precon import PreconLBFGS
-
-from matscipy import parameter
-from matscipy.elasticity import  fit_elastic_constants
-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_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)
-ds = parameter('ds', 1e-2)
-nsteps = parameter('nsteps', 10)
-a0 = parameter('a0') # lattice constant
-k0 = parameter('k0', 1.0)
-extended_far_field = parameter('extended_far_field', False)
-alpha0 = parameter('alpha0', 0.0) # initial guess for crack position
-dump = parameter('dump', False)
-precon = parameter('precon', False)
-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')
-
-
-# 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'))
-
-# setup the crack
-crk = CubicCrystalCrack(parameter('crack_surface'),
-                        parameter('crack_front'),
-                        Crot=C/GPa)
-
-# get Griffith's stess intensity factor
-k1g = crk.k1g(parameter('surface_energy'))
-print('Griffith k1 = %f' % k1g)
-
-# check for restart files if continuation=True
-if continuation and not os.path.exists(restart_file):
-    continuation = False
-if continuation:
-
-    # 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 ; 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:
-
-    # 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)
-
-    traj = None
-    dk = parameter('dk', 1e-4)
-    dalpha = parameter('dalpha', 1e-3)
-
-    # reuse output from sinclair_crack.py if possible
-    if os.path.exists(f'k_{int(k0 * 1000):04d}.xyz'):
-        print(f'Reading atoms from k_{int(k0 * 1000):04d}.xyz')
-        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
-
-    
-    # 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
-    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
-        sc.variable_alpha = False
-        sc.variable_k = False
-        sc.update_atoms()
-        atoms = sc.atoms.copy()
-        atoms.calc = sc.calc
-        atoms.set_constraint(FixAtoms(mask=~sc.regionI))
-        opt = PreconLBFGS(atoms, logfile=None)
-        #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}')
-        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()
-
-    # 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_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
-
-    # 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]
-    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
-
-# 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,
-                    extended_far_field=extended_far_field)
-
-# run full arc length continuation
-scv.arc_length_continuation(x0, x1, N=nsteps,
-                        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
diff --git a/scripts/fracture_mechanics/sinclair_crack.py b/scripts/fracture_mechanics/sinclair_crack.py
deleted file mode 100644
index bc2518a0..00000000
--- a/scripts/fracture_mechanics/sinclair_crack.py
+++ /dev/null
@@ -1,149 +0,0 @@
-#
-# Copyright 2014, 2020 James Kermode (Warwick U.)
-#           2020 Arnaud Allera (U. Lyon 1)
-#           2014 Lars Pastewka (U. Freiburg)
-#
-# matscipy - Materials science with Python at the atomic-scale
-# https://github.com/libAtoms/matscipy
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 2 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
-
-import os
-import numpy as np
-
-import ase.io
-from ase.units import GPa
-from ase.constraints import FixAtoms
-from ase.optimize.precon import PreconLBFGS
-
-from matscipy import parameter
-from matscipy.elasticity import fit_elastic_constants
-from matscipy.fracture_mechanics.crack import CubicCrystalCrack, SinclairCrack
-
-import sys
-
-sys.path.insert(0, '.')
-import params
-
-calc = parameter('calc')
-fmax = parameter('fmax', 1e-3)
-vacuum = parameter('vacuum', 10.0)
-flexible = parameter('flexible', True)
-extended_far_field = parameter('extended_far_field', False)
-k0 = parameter('k0', 1.0)
-alpha0 = parameter('alpha0', 0.0) # initial guess for crack position
-dump = parameter('dump', False)
-precon = parameter('precon', False)
-prerelax = parameter('prerelax', False)
-lbfgs = parameter('lbfgs', not flexible) # use LBGS by default if not flexible
-
-# compute elastic constants
-cryst = params.cryst.copy()
-cryst.pbc = True
-cryst.calc = calc
-C, C_err = fit_elastic_constants(cryst,
-                                 symmetry=parameter('elastic_symmetry',
-                                                    'triclinic'))
-
-crk = CubicCrystalCrack(parameter('crack_surface'),
-                        parameter('crack_front'),
-                        Crot=C/GPa)
-
-# Get Griffith's k1.
-k1g = crk.k1g(parameter('surface_energy'))
-print('Griffith k1 = %f' % k1g)
-
-
-cluster = params.cluster.copy()
-sc = SinclairCrack(crk, cluster, calc, k0 * k1g,
-                   alpha=alpha0,
-                   variable_alpha=flexible,
-                   vacuum=vacuum,
-                   extended_far_field=extended_far_field)
-
-nsteps = parameter('nsteps')
-k1_range = parameter('k1_range', np.linspace(0.8, 1.2, nsteps))
-
-ks = list(k1_range)
-ks_out = []
-alphas = []
-
-max_steps = parameter('max_steps', 10)
-fit_alpha = parameter('fit_alpha', False)
-
-for i, k in enumerate(ks):
-    restart_file = parameter('restart_file', f'k_{int(k * 1000):04d}.xyz')
-    if os.path.exists(restart_file):
-        print(f'Restarting from {restart_file}')
-        if restart_file.endswith('.xyz'):
-            a = ase.io.read(restart_file)
-            sc.set_atoms(a)
-        elif restart_file.endswith('.txt'):
-            x = np.loadtxt(restart_file)
-            if len(x) == len(sc):
-                sc.set_dofs(x)
-            elif len(x) == len(sc) + 1:
-                sc.variable_k = True
-                sc.set_dofs(x)
-                sc.variable_k = False
-            elif len(x) == len(sc) + 2:
-                sc.variable_alpha = True
-                sc.variable_k = True
-                sc.set_dofs(x)
-                sc.variable_alpha = False
-                sc.variable_k = False
-            else:
-                raise RuntimeError('cannot guess how to restart with'
-                                   f' {len(x)} variables for {len(self)} DoFs')
-        else:
-            raise ValueError(f"don't know how to restart from {restart_file}")
-    else:
-        sc.rescale_k(k * k1g)
-
-    if fit_alpha:
-        sc.alpha, = sc.fit_cle(variable_alpha=True, variable_k=False)
-        print(f'Fitted value of alpha: {sc.alpha}')
-    print(f'k = {sc.k / k1g} * k1g')
-    print(f'alpha = {sc.alpha}')
-
-    if lbfgs:
-        print('Optimizing with LBFGS')
-        atoms = sc.atoms.copy()
-        atoms.calc = sc.calc
-        atoms.set_constraint(FixAtoms(mask=np.logical_not(sc.regionI)))
-        opt = PreconLBFGS(atoms)
-        opt.run(fmax=fmax)
-        sc.set_atoms(atoms)
-    else:
-        if prerelax:
-            print('Pre-relaxing with Conjugate-Gradients')
-            sc.optimize(ftol=1e-5, steps=max_steps, dump=dump,
-                        method='cg')
-
-        sc.optimize(fmax, steps=max_steps, dump=dump,
-                    precon=precon, method='krylov')
-
-    if flexible:
-        print(f'Optimized alpha = {sc.alpha:.3f}')
-    ks_out.append(k)
-    alphas.append(sc.alpha)
-
-    a = sc.atoms
-    a.get_forces()
-    ase.io.write(f'k_{int(k * 1000):04d}.xyz', a)
-    with open('x.txt', 'a') as fh:
-        np.savetxt(fh, [sc.get_dofs()])
-
-np.savetxt('k_vs_alpha.txt', np.c_[ks_out, alphas])