diff --git a/matscipy/fracture_mechanics/crack.py b/matscipy/fracture_mechanics/crack.py
index 98512837..07829ae6 100644
--- a/matscipy/fracture_mechanics/crack.py
+++ b/matscipy/fracture_mechanics/crack.py
@@ -45,7 +45,7 @@
                                  Voigt_6_to_full_3x3_stress)
 from matscipy.surface import MillerDirection, MillerPlane
 from matscipy.neighbours import neighbour_list
-from matscipy.fracture_mechanics.solvers import ode12r
+from matscipy.optimize import ode12r
 
 import matplotlib.pyplot as plt
 from ase.optimize.precon import Exp, PreconLBFGS
@@ -1849,9 +1849,9 @@ def strain_err(self,cutoff,seperate_surface=False):
             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))
+            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)
diff --git a/matscipy/fracture_mechanics/meson.build b/matscipy/fracture_mechanics/meson.build
index 7a01e6df..9cef288d 100644
--- a/matscipy/fracture_mechanics/meson.build
+++ b/matscipy/fracture_mechanics/meson.build
@@ -7,8 +7,7 @@ python_sources = [
     'crack.py',
     'crackpathsel.py',
     'energy_release.py',
-    'idealbrittlesolid.py',
-    'solvers.py'
+    'idealbrittlesolid.py'
 ]
 
 # Install pure Python
diff --git a/matscipy/meson.build b/matscipy/meson.build
index dc597d26..a8e6f488 100644
--- a/matscipy/meson.build
+++ b/matscipy/meson.build
@@ -28,6 +28,7 @@ python_sources = [
     'numerical.py',
     'numpy_tricks.py',
     'opls.py',
+    'optimize.py',
     'precon.py',
     'pressurecoupling.py',
     'rings.py',
diff --git a/matscipy/fracture_mechanics/solvers.py b/matscipy/optimize.py
similarity index 100%
rename from matscipy/fracture_mechanics/solvers.py
rename to matscipy/optimize.py
diff --git a/matscipy/surface.py b/matscipy/surface.py
index 4fbbd271..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"""
@@ -327,11 +328,10 @@ def find_surface_energy(symbol,calc,a0,surface,size=(8,1,1),vacuum=10,fmax=0.000
     if unit == 'ASE':
         return [gamma_ase,'ase_units']
     else:
-        from ase import units
         gamma_SI = (gamma_ase / units.J ) * (units.m)**2
         if unit =='J/m^2':
-            return [gamma_SI,'J/m^2']
+            return gamma_SI
         elif unit == '0.1J/m^2':
-            return [10*gamma_SI,'0.1J/m^2'] # units required for the fracture code
+            return 10*gamma_SI # units required for the fracture code
         else:
             print('Error: Unsupported unit of surface energy.')