diff --git a/src/IO.f90 b/src/IO.f90
index be8982aeb..1029e3360 100644
--- a/src/IO.f90
+++ b/src/IO.f90
@@ -463,10 +463,6 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
 
 !--------------------------------------------------------------------------------------------------
 ! math errors
-    case (400)
-      msg = 'matrix inversion error'
-    case (401)
-      msg = 'error in Eigenvalue calculation'
     case (402)
       msg = 'invalid orientation specified'
 
diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90
index 36bbaeda3..dcd2c98cb 100644
--- a/src/grid/spectral_utilities.f90
+++ b/src/grid/spectral_utilities.f90
@@ -713,7 +713,6 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
     allocate(s_reduced,mold = c_reduced)
     call math_invert(s_reduced, errmatinv, c_reduced)                                               ! invert reduced stiffness
     if (any(IEEE_is_NaN(s_reduced))) errmatinv = .true.
-    if (errmatinv) call IO_error(error_ID=400,ext_msg='utilities_maskedCompliance')
 
 !--------------------------------------------------------------------------------------------------
 ! check if inversion was successful
@@ -725,7 +724,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
       write(IO_STDOUT,trim(formatString),advance='no') ' C * S (load) ', &
                                                              transpose(matmul(c_reduced,s_reduced))
       write(IO_STDOUT,trim(formatString),advance='no') ' S (load) ', transpose(s_reduced)
-      if(errmatinv) call IO_error(error_ID=400,ext_msg='utilities_maskedCompliance')
+      if(errmatinv) error stop 'matrix inversion error'
     endif
     temp99_real = reshape(unpack(reshape(s_reduced,[size_reduced**2]),reshape(mask,[81]),0.0_pReal),[9,9])
   else
diff --git a/src/lattice.f90 b/src/lattice.f90
index 6cfa41fef..08385eac7 100644
--- a/src/lattice.f90
+++ b/src/lattice.f90
@@ -2209,11 +2209,11 @@ function equivalent_nu(C,assumption) result(nu)
       / 9.0_pReal
   elseif(IO_lc(assumption) == 'reuss') then
     call math_invert(S,error,C)
-    if(error) call IO_error(0)
+    if(error) error stop 'matrix inversion failed'
     K = 1.0_pReal &
       / (S(1,1)+S(2,2)+S(3,3) +2.0_pReal*(S(1,2)+S(2,3)+S(1,3)))
   else
-    call IO_error(0)
+    error stop 'invalid assumption'
     K = 0.0_pReal
   endif
 
@@ -2241,11 +2241,11 @@ function equivalent_mu(C,assumption) result(mu)
        / 15.0_pReal
   elseif(IO_lc(assumption) == 'reuss') then
     call math_invert(S,error,C)
-    if(error) call IO_error(0)
+    if(error) error stop 'matrix inversion failed'
     mu = 15.0_pReal &
        / (4.0_pReal*(S(1,1)+S(2,2)+S(3,3)) -4.0_pReal*(S(1,2)+S(2,3)+S(1,3)) +3.0_pReal*(S(4,4)+S(5,5)+S(6,6)))
   else
-    call IO_error(0)
+    error stop 'invalid assumption'
     mu = 0.0_pReal
   endif
 
diff --git a/src/math.f90 b/src/math.f90
index 163f4df6a..b01cf9e26 100644
--- a/src/math.f90
+++ b/src/math.f90
@@ -499,7 +499,7 @@ function math_invSym3333(A)
   call dgetrf(6,6,temp66,6,ipiv6,ierr_i)
   call dgetri(6,temp66,6,ipiv6,work,size(work,1),ierr_f)
   if (ierr_i /= 0 .or. ierr_f /= 0) then
-    call IO_error(400, ext_msg = 'math_invSym3333')
+    error stop 'matrix inversion error'
   else
     math_invSym3333 = math_66toSym3333(temp66)
   endif
diff --git a/src/rotations.f90 b/src/rotations.f90
index ea4a8a9d8..888e73762 100644
--- a/src/rotations.f90
+++ b/src/rotations.f90
@@ -640,13 +640,13 @@ function om2ax(om) result(ax)
     ax(1:3) = [ 0.0_pReal, 0.0_pReal, 1.0_pReal ]
   else
     call dgeev('N','V',3,om_,3,Wr,Wi,devNull,3,VR,3,work,size(work,1),ierr)
-    if (ierr /= 0) call IO_error(401,ext_msg='Error in om2ax: DGEEV return not zero')
+    if (ierr /= 0) error stop 'LAPACK error'
 #if defined(__GFORTRAN__) &&  __GNUC__<9 || defined(__INTEL_COMPILER) && INTEL_COMPILER<1800 || defined(__PGI)
     i = maxloc(merge(1,0,cEq(cmplx(Wr,Wi,pReal),cmplx(1.0_pReal,0.0_pReal,pReal),tol=1.0e-14_pReal)),dim=1)
 #else
     i = findloc(cEq(cmplx(Wr,Wi,pReal),cmplx(1.0_pReal,0.0_pReal,pReal),tol=1.0e-14_pReal),.true.,dim=1) !find eigenvalue (1,0)
 #endif
-    if (i == 0) call IO_error(401,ext_msg='Error in om2ax Real: eigenvalue not found')
+    if (i == 0) error stop 'om2ax conversion failed'
     ax(1:3) = VR(1:3,i)
     where (                dNeq0([om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)])) &
       ax(1:3) = sign(ax(1:3),-P *[om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)])