diff --git a/.gitignore b/.gitignore
index deb19d0..9980d03 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,6 +1,10 @@
 *.o
 *.mod
 *.dbg
+*.cmdx
+*.cmod
+*.ilm
+*.stb
 .vscode
 Slater*
 Updates*
diff --git a/Makefile b/Makefile
index 3af3823..3673340 100644
--- a/Makefile
+++ b/Makefile
@@ -50,7 +50,7 @@ EXEC := $(BIN_DIR)/cMaponiA3_test_3x3_3 \
 all: $(EXEC)
 
 clean:
-	@rm -vrf $(OBJ_DIR) *.dbg
+	@rm -vrf $(OBJ_DIR) *.dbg *.cmdx *.cmod *.ilm *.stb
 
 distclean: clean
 	@rm -vrf $(BIN_DIR) \
diff --git a/smvars.sh b/smvars.sh
index 447f7bb..06cc242 100644
--- a/smvars.sh
+++ b/smvars.sh
@@ -2,12 +2,13 @@
 
 ENV=$1
 
+## Set Sherman-Morrison root dir
 PWD=$(pwd)
 SRCDIR=$(dirname $BASH_SOURCE)
 case $SRCDIR in
-	/*) SMROOT=$SRCDIR ;; ## sourced from absolute path
-	*) ## sourced from absolute path
-		if [[ $SRCDIR = . ]] ## check if already in root
+	/*) SMROOT=$SRCDIR ;;
+	*)
+		if [[ $SRCDIR = . ]]
 		then
 			SMROOT=$PWD
 		else
@@ -17,6 +18,7 @@ case $SRCDIR in
 esac
 export SMROOT
 
+## Set environment for hdf5-tools and Makefile
 case $ENV in
   intel)
 	export HDF5_CXX=icpc
@@ -40,6 +42,7 @@ case $ENV in
 	;;
 esac
 
+## Export path, but only once
 if [[ -z $SMVARS ]]
 then
 	export PATH=$SMROOT/bin:$PATH
diff --git a/src/SM_MaponiA3.cpp b/src/SM_MaponiA3.cpp
index 0d55188..df4540f 100644
--- a/src/SM_MaponiA3.cpp
+++ b/src/SM_MaponiA3.cpp
@@ -102,7 +102,7 @@ void MaponiA3(double *Slater_inv, unsigned int Dim,
     cout << "beta = 1 + ylk[" << l << "][" << p[l+1] << "][" << component << "] = " << beta << endl;
     cout << endl;
     #endif
-    if (fabs(beta) < 1e-6) {
+    if (fabs(beta) < 1e-3) {
       cerr << "Break-down occured." << endl;
     }
 
diff --git a/src/SM_Standard.cpp b/src/SM_Standard.cpp
index 4082bd2..8ac8b72 100644
--- a/src/SM_Standard.cpp
+++ b/src/SM_Standard.cpp
@@ -4,10 +4,16 @@
 #include "Helpers.hpp"
 #include <iostream>
 
+// Set common break-down threshold
+double threshold() {
+  double const threshold = 1e-3;
+  return threshold;
+}
 
 // Naïve Sherman Morrison
 void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
-        double *Updates, unsigned int *Updates_index) {
+         double *Updates, unsigned int *Updates_index)
+{
 
   std::cerr << "Called SM1 with " << N_updates << " updates" << std::endl;
   double C[Dim];
@@ -15,30 +21,37 @@ void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
 
   unsigned int l = 0;
   // For each update
-  while (l < N_updates) {
+  while (l < N_updates)
+  {
     // C = A^{-1} x U_l
-    for (unsigned int i = 0; i < Dim; i++) {
+    for (unsigned int i = 0; i < Dim; i++)
+    {
       C[i] = 0;
-      for (unsigned int j = 0; j < Dim; j++) {
+      for (unsigned int j = 0; j < Dim; j++)
+      {
         C[i] += Slater_inv[i * Dim + j] * Updates[l * Dim + j];
       }
     }
 
     // Denominator
     double den = 1 + C[Updates_index[l] - 1];
-    if (fabs(den) < 1e-6) {
-      std::cerr << "Breakdown condition triggered at " << Updates_index[l]  << std::endl;
+    if (fabs(den) < threshold())
+    {
+      std::cerr << "Breakdown condition triggered at " << Updates_index[l] << std::endl;
     }
     double iden = 1 / den;
 
     // D = v^T x A^{-1}
-    for (unsigned int j = 0; j < Dim; j++) {
+    for (unsigned int j = 0; j < Dim; j++)
+    {
       D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j];
     }
 
     // A^{-1} = A^{-1} - C x D / den
-    for (unsigned int i = 0; i < Dim; i++) {
-      for (unsigned int j = 0; j < Dim; j++) {
+    for (unsigned int i = 0; i < Dim; i++)
+    {
+      for (unsigned int j = 0; j < Dim; j++)
+      {
         double update = C[i] * D[j] * iden;
         Slater_inv[i * Dim + j] -= update;
       }
@@ -51,52 +64,61 @@ void SM1(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
 // Sherman Morrison, with J. Slagel splitting
 // http://hdl.handle.net/10919/52966
 void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
-        double *Updates, unsigned int *Updates_index) {
+         double *Updates, unsigned int *Updates_index)
+{
 
   std::cerr << "Called SM2 with updates " << N_updates << std::endl;
   double C[Dim];
   double D[Dim];
 
-  double later_updates[Dim*N_updates];
+  double later_updates[Dim * N_updates];
   unsigned int later_index[N_updates];
   unsigned int later = 0;
 
   unsigned int l = 0;
   // For each update
-  while (l < N_updates) {
+  while (l < N_updates)
+  {
     // C = A^{-1} x U_l
-    for (unsigned int i = 0; i < Dim; i++) {
+    for (unsigned int i = 0; i < Dim; i++)
+    {
       C[i] = 0;
-      for (unsigned int j = 0; j < Dim; j++) {
+      for (unsigned int j = 0; j < Dim; j++)
+      {
         C[i] += Slater_inv[i * Dim + j] * Updates[l * Dim + j];
       }
     }
 
     // Denominator
     double den = 1 + C[Updates_index[l] - 1];
-    if (fabs(den) < 1e-6) {
-      std::cerr << "Breakdown condition triggered at " << Updates_index[l]  << std::endl;
+    if (fabs(den) < threshold())
+    {
+      std::cerr << "Breakdown condition triggered at " << Updates_index[l] << std::endl;
 
       // U_l = U_l / 2 (do the split)
-      for (unsigned int j = 0; j < Dim; j++) {
-        later_updates[later*Dim+j] = Updates[l*Dim+j]/2.0;
+      for (unsigned int j = 0; j < Dim; j++)
+      {
+        later_updates[later * Dim + j] = Updates[l * Dim + j] / 2.0;
         C[j] /= 2.0;
       }
       later_index[later] = Updates_index[l];
       later++;
 
-      den = 1+C[Updates_index[l]-1];
+      den = 1 + C[Updates_index[l] - 1];
     }
     double iden = 1 / den;
 
     // D = v^T x A^{-1}
-    for (unsigned int j = 0; j < Dim; j++) {
+    for (unsigned int j = 0; j < Dim; j++)
+    {
       D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j];
     }
 
     // A^{-1} = A^{-1} - C x D / den
-    for (unsigned int i = 0; i < Dim; i++) {
-      for (unsigned int j = 0; j < Dim; j++) {
+    for (unsigned int i = 0; i < Dim; i++)
+    {
+      for (unsigned int j = 0; j < Dim; j++)
+      {
         double update = C[i] * D[j] * iden;
         Slater_inv[i * Dim + j] -= update;
       }
@@ -104,41 +126,48 @@ void SM2(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
     l += 1;
   }
 
-  if (later > 0) {
+  if (later > 0)
+  {
     SM2(Slater_inv, Dim, later, later_updates, later_index);
   }
 }
 
 // Sherman Morrison, leaving zero denominators for later
 void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
-        double *Updates, unsigned int *Updates_index) {
+         double *Updates, unsigned int *Updates_index)
+{
 
   std::cerr << "Called SM3 with updates " << N_updates << std::endl;
   double C[Dim];
   double D[Dim];
 
-  double later_updates[Dim*N_updates];
+  double later_updates[Dim * N_updates];
   unsigned int later_index[N_updates];
   unsigned int later = 0;
 
   unsigned int l = 0;
   // For each update
-  while (l < N_updates) {
+  while (l < N_updates)
+  {
     // C = A^{-1} x U_l
-    for (unsigned int i = 0; i < Dim; i++) {
+    for (unsigned int i = 0; i < Dim; i++)
+    {
       C[i] = 0;
-      for (unsigned int j = 0; j < Dim; j++) {
+      for (unsigned int j = 0; j < Dim; j++)
+      {
         C[i] += Slater_inv[i * Dim + j] * Updates[l * Dim + j];
       }
     }
 
     // Denominator
     double den = 1 + C[Updates_index[l] - 1];
-    if (fabs(den) < 1e-6) {
-      std::cerr << "Breakdown condition triggered at " << Updates_index[l]  << std::endl;
+    if (fabs(den) < threshold())
+    {
+      std::cerr << "Breakdown condition triggered at " << Updates_index[l] << std::endl;
 
-      for (unsigned int j = 0; j < Dim; j++) {
-        later_updates[later*Dim+j] = Updates[l*Dim+j];
+      for (unsigned int j = 0; j < Dim; j++)
+      {
+        later_updates[later * Dim + j] = Updates[l * Dim + j];
       }
       later_index[later] = Updates_index[l];
       later++;
@@ -148,13 +177,16 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
     double iden = 1 / den;
 
     // D = v^T x A^{-1}
-    for (unsigned int j = 0; j < Dim; j++) {
+    for (unsigned int j = 0; j < Dim; j++)
+    {
       D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j];
     }
 
     // A^{-1} = A^{-1} - C x D / den
-    for (unsigned int i = 0; i < Dim; i++) {
-      for (unsigned int j = 0; j < Dim; j++) {
+    for (unsigned int i = 0; i < Dim; i++)
+    {
+      for (unsigned int j = 0; j < Dim; j++)
+      {
         double update = C[i] * D[j] * iden;
         Slater_inv[i * Dim + j] -= update;
       }
@@ -162,7 +194,8 @@ void SM3(double *Slater_inv, unsigned int Dim, unsigned int N_updates,
     l += 1;
   }
 
-  if (later > 0) {
+  if (later > 0)
+  {
     SM3(Slater_inv, Dim, later, later_updates, later_index);
   }
 }
diff --git a/src/det.irp.f b/src/det.irp.f
index e83d64a..3189986 100644
--- a/src/det.irp.f
+++ b/src/det.irp.f
@@ -452,7 +452,13 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
   !DIR$ ASSUME (LDS >= $n)
   integer                        :: i,j
   double precision :: zj, zj1, zj2, zj3
+
+  101 format ('Updating column ', I2)
+  102 format ('#BREAKDOWN_OCCURED: TH=1.d-3')
+
+  write(502, '("In subroutine DET_UPDATE21()")')
   
+  write(502,101) l
   do i=1,$n
     u(i) = m(i) - S(i,l)
   enddo
@@ -469,9 +475,9 @@ subroutine det_update$n(n,LDS,m,l,S,S_inv,d)
   d_inv = 1.d0/d
   d = d+zj
   lambda = d*d_inv
-  if (dabs(lambda) < 1.d-6) then
+  if (dabs(lambda) < 1.d-3) then
     d = 0.d0
-    write(502,"('#BREAKDOWN_OCCURED')")
+    write(502,102)
     return
   end if
   
@@ -970,7 +976,7 @@ subroutine det_update_general(n,LDS,m,l,S,S_inv,d)
   !DIR$ ASSUME (MOD(LDS,$IRP_ALIGN/8) == 0)
   !DIR$ ASSUME (n>150)
 
-  integer          :: i,j,n4
+  integer                        :: i,j,n4
   double precision :: zl
 
   !DIR$ NOPREFETCH
@@ -1168,7 +1174,7 @@ END_PROVIDER
 
   !! Some usefull formats for output
   10001 format ('#START_PACKET')
-  10008 format ('#CYCLE_ID: ', I4)
+  10008 format ('#CYCLE_ID: ', I5)
   10000 format ('#SLATER_MATRIX_DIM: ', I3)
   10002 format ('#NUPDATES: ', I2)
   10003 format ('#SLATER_MATRIX: (i (outer), j (inner)), slater_matrix_alpha(i,j), ddet * slater_matrix_alpha_inv_det(i,j) / ddet')
@@ -1177,8 +1183,9 @@ END_PROVIDER
   10006 format ('#COL_UPDATE_COMP_(',I0.2,'): ', E23.15)
   10007 format ('#END_PACKET',/)
 
-  open (unit = 501, file = "dataset.dat") !! slightly cleaner output
+  open (unit = 501, file = "dataset.dat") !! Only Slater matrices plus updates
   open (unit = 502, file = "dataset.fulltrace.dat")
+  write(502,'("In subroutine BLD_DET_ALPHA_VALUE_CURR(...)")')
 
   if (ifirst == 0) then !! If this is the first time we enter this subroutine
     ifirst = 1
@@ -1191,7 +1198,6 @@ END_PROVIDER
   PROVIDE mo_value
   if (det_i /= det_alpha_order(1) ) then
     ! write(*,*) "det_i: ", det_i
-!  if (det_i == -1 ) then
 
     !! determin number of updates, the updates and 
     n_to_do = 0
@@ -1225,38 +1231,6 @@ END_PROVIDER
       end do
     end do
  
-!  print n_to_do (number of updates to do)
-!  to_do (array of the columns that need to be swapped)
-!  mo_list_alpha_curr (list of orbitals to build the current determinant)
-!  mo_list_alpha_prev (list of the previous determinant)
-!
-!  slater_matrix_alpha (slater matrix that needs to be inverted. This is the initial matrix)
-!  slater_matrix_alpha_inv_det = inverse of the slater matrix divided by the determinant (ddet: l. 1216)
-!
-!
-!
-!       1 2 3  4
-!       --------
-!       2 4 6  8
-!       2 4 10 8
-!    
-!     mo_list(:)       (2,4,10,8)
-!     mo_list_prev(:)  (2,4,6,8)
-!     n_todo           1
-!     to_do            (3)
-!
-!
-!       1 2 3  4
-!       --------
-!       2 4 6  8
-!       2 5 10 8
-!    
-!     mo_list(:)       (2,5,10,8)
-!     mo_list_prev(:)  (2,4,6,8)
-!     n_todo           2
-!     to_do            (2,3)
-
-
     ddet = 0.d0
     
     if (n_to_do < shiftl(elec_alpha_num,1)) then
@@ -1291,6 +1265,7 @@ END_PROVIDER
               slater_matrix_alpha,                                     &
               slater_matrix_alpha_inv_det,                             &
               ddet)
+          write(502, '("BACK in subroutine BLD_DET_ALPHA_VALUE_CURR(...)")')
           if (ddet /= 0.d0) then
             det_alpha_value_curr = ddet
           else
@@ -1305,15 +1280,9 @@ END_PROVIDER
         endif
       enddo
     endif
-
-  write(501,10007)
-  write(502,10007)
-  cycle_id = cycle_id + 1
   
   else
-
     ddet = 0.d0
-
   endif
 
     
@@ -1337,12 +1306,24 @@ END_PROVIDER
         slater_matrix_alpha_inv_det(k,i) = mo_value(i,mo_list_alpha_curr(k))
       enddo
     enddo
+    if (det_i /= det_alpha_order(1)) then
+      write(502,"('#NO_INVERSE_FOUND: ')") cycle_id
+    else
+      write(502,"('#FIRST_DETERMINANT_IN_SUPERCYCLE')")
+    endif
     call invert(slater_matrix_alpha_inv_det,elec_alpha_num_8,elec_alpha_num,ddet)
-    
   endif
+
+  write(501,10007)
+  write(502,10007)
+  cycle_id = cycle_id + 1
+
   ASSERT (ddet /= 0.d0)
   
   det_alpha_value_curr = ddet 
+
+
+
 END_PROVIDER
 
  BEGIN_PROVIDER [ double precision, det_beta_value_curr ]
diff --git a/tests/run_test.sh b/tests/run_test.sh
new file mode 100755
index 0000000..01c425b
--- /dev/null
+++ b/tests/run_test.sh
@@ -0,0 +1,41 @@
+#!/usr/bin/env bash
+
+## IMPORTANT: start cycle will always be 1 because of the intention of this
+## script and the way it is written now. In fact, var. $START is not used at
+## the moment.
+
+## Call: $ ./run_test.sh <{sm1|sm2|sm3|maponia3}> <start cycle> <stop cycle>
+## Example: ./run_test.sh sm2 1 500
+
+TEST=test_h5
+ALGO=$1
+START=$2
+STOP=$3
+BLOCKSIZE=329
+
+BLOCKS=$((STOP / BLOCKSIZE))
+LEFT=$((STOP % BLOCKSIZE))
+
+if [[ $LEFT -ne 0 ]]
+then
+    BSTART=0
+    BSTOP=0
+    LAST=0
+    for ((i=1; i<=$BLOCKS; i++))
+    do
+        BSTART=$(((i-1)*BLOCKSIZE+1))
+        BSTOP=$((i*BLOCKSIZE-1))
+        $TEST $ALGO $BSTART $BSTOP
+        LAST=$i
+    done
+    LSTART=$((LAST*BLOCKSIZE+1))
+    LSTOP=$((LSTART+LEFT-1))
+    $TEST $ALGO $LSTART $LSTOP
+else
+    for ((i=1; i<=$BLOCKS; i++))
+    do
+        BSTART=$(((i-1)*BLOCKSIZE+1))
+        BSTOP=$((i*BLOCKSIZE-1))
+        $TEST $ALGO $BSTART $BSTOP
+    done
+fi
diff --git a/tests/test_h5.cpp b/tests/test_h5.cpp
index 0e38a43..b7ad2e1 100644
--- a/tests/test_h5.cpp
+++ b/tests/test_h5.cpp
@@ -13,7 +13,7 @@ using namespace H5;
 const H5std_string FILE_NAME( "dataset.hdf5" );
 
 double residual_max(double * A, unsigned int Dim) {
-  double max= 0.0;
+  double max = 0.0;
   for (unsigned int i = 0; i < Dim; i++) {
     for (unsigned int j = 0; j < Dim; j++) {
       double delta = (A[i * Dim + j] - (i == j));
diff --git a/tools/stats.py b/tools/stats.py
new file mode 100755
index 0000000..f685547
--- /dev/null
+++ b/tools/stats.py
@@ -0,0 +1,26 @@
+#!/usr/bin/env python
+
+import sys
+import numpy as np
+
+filename = sys.argv[1]
+
+Cycle, Rmax, R2 = np.loadtxt(filename, usecols = (3, 4, 5), unpack = True)
+
+print ("========================================================")
+print ("Statistics for file: ", filename)
+print ("========================================================")
+print ("Max element residual")
+print ("--------------------")
+print ("MINIMUM, Cycle#: ", np.nanmin(Rmax), ", ", Cycle[np.nanargmin(Rmax)])
+print ("MAXIMUM, Cycle#: ", np.nanmax(Rmax), ", " , Cycle[np.nanargmax(Rmax)])
+print ("MEAN: ", np.nanmean(Rmax))
+print ("STD: ", np.nanstd(Rmax))
+print ("")
+print ("Frobenius norm squared residual")
+print ("--------------------------------")
+print ("MINIMUM, Cycle#: ", np.nanmin(R2), ", ", Cycle[np.nanargmin(R2)])
+print ("MAXIMUM, Cycle#: ", np.nanmax(R2), ", " , Cycle[np.nanargmax(R2)])
+print ("MEAN: ", np.nanmean(R2))
+print ("STD: ", np.nanstd(R2))
+print ("")