diff --git a/build/Makefile.in b/build/Makefile.in
index aa9ce7570396f2958d25380e1d5d129521e2cacb..dfc4f1b8f9bb056082ad149aa35ba966718afbcd 100644
--- a/build/Makefile.in
+++ b/build/Makefile.in
@@ -456,6 +456,7 @@ CXXCPP = @CXXCPP@
 CXXDEPMODE = @CXXDEPMODE@
 CXXFLAGS = @CXXFLAGS@
 CYGPATH_W = @CYGPATH_W@
+DEBUGFLAGS = @DEBUGFLAGS@
 DEFS = @DEFS@
 DEPDIR = @DEPDIR@
 DLLTOOL = @DLLTOOL@
diff --git a/build/configure b/build/configure
index 5fa0692ffb1fd46bed6a5ab4ad5bd3d331f365a9..72b47b22d9fe3727ffb8a279f2c62efa15d53fae 100755
--- a/build/configure
+++ b/build/configure
@@ -669,6 +669,7 @@ OFFLOADFLAGS
 BUILDFORTRAN_FALSE
 BUILDFORTRAN_TRUE
 ENABLE_ILP64
+DEBUGFLAGS
 HDF5_LDFLAGS
 CXXCPP
 LT_SYS_LIBRARY_PATH
@@ -813,6 +814,7 @@ with_aix_soname
 with_gnu_ld
 with_sysroot
 enable_libtool_lock
+enable_debug
 enable_ilp64
 enable_fortran
 enable_offload
@@ -1479,6 +1481,7 @@ Optional Features:
   --enable-fast-install[=PKGS]
                           optimize for fast installation [default=yes]
   --disable-libtool-lock  avoid locking (might break parallel builds)
+  --enable-debug=FEATURE  enable debug of FEATURE [default=no]
   --enable-ilp64          enable 64-bit indexing [default=yes]
   --enable-fortran        enable legacy FORTRAN compilation [default=auto]
   --enable-offload        enable target offloading (requires g++ version >=
@@ -24936,6 +24939,25 @@ esac
 fi
 
 # Configure the optional features
+# Check whether --enable-debug was given.
+if test ${enable_debug+y}
+then :
+  enableval=$enable_debug;
+    if test "x$enableval" = "xAM"; then
+      DEBUGFLAGS="-DDEBUG_AM"
+
+    fi
+
+else case e in #(
+  e)
+    DEBUGFLAGS=""
+
+
+ ;;
+esac
+fi
+
+
 # Check whether --enable-ilp64 was given.
 if test ${enable_ilp64+y}
 then :
@@ -25647,7 +25669,7 @@ else case e in #(
  ;;
 esac
 fi
-CXXFLAGS="$CLANGFLAGS $OPTFLAGS -ggdb $OFFLOADFLAGS $USER_INCLUDE -I$HDF5_INCLUDE $OMPFLAGS $MPIFLAGS $LAPACKFLAGS $MAGMAFLAGS $NVTXFLAGS"
+CXXFLAGS="$CLANGFLAGS $OPTFLAGS -ggdb $DEBUGFLAGS $OFFLOADFLAGS $USER_INCLUDE -I$HDF5_INCLUDE $OMPFLAGS $MPIFLAGS $LAPACKFLAGS $MAGMAFLAGS $NVTXFLAGS"
 SUBDIRS="cluster libnptm sphere testing trapping"
 
 # Generate the output
diff --git a/build/configure.ac b/build/configure.ac
index 7ff1a35c695fef8c5c07baefcd69ee239d4e48eb..0cc9b84493d563d6f9466f27a2bf1d1b01612cbb 100644
--- a/build/configure.ac
+++ b/build/configure.ac
@@ -405,6 +405,19 @@ AS_IF(
 )
 
 # Configure the optional features
+AC_ARG_ENABLE(
+  [debug],
+  [AS_HELP_STRING([--enable-debug=FEATURE], [enable debug of FEATURE [default=no]])],
+  [
+    if test "x$enableval" = "xAM"; then
+      AC_SUBST([DEBUGFLAGS], ["-DDEBUG_AM"])
+    fi
+  ],
+  [
+    AC_SUBST([DEBUGFLAGS], [""])
+  ]
+)
+
 AC_ARG_ENABLE(
   [ilp64],
   [AS_HELP_STRING([--enable-ilp64], [enable 64-bit indexing [default=yes]])],
@@ -641,7 +654,7 @@ AS_IF(
   [AC_SUBST([OMPFLAGS], [""])],
   [AC_SUBST([OMPFLAGS], [$OMPFLAGS])]
 )
-CXXFLAGS="$CLANGFLAGS $OPTFLAGS -ggdb $OFFLOADFLAGS $USER_INCLUDE -I$HDF5_INCLUDE $OMPFLAGS $MPIFLAGS $LAPACKFLAGS $MAGMAFLAGS $NVTXFLAGS"
+CXXFLAGS="$CLANGFLAGS $OPTFLAGS -ggdb $DEBUGFLAGS $OFFLOADFLAGS $USER_INCLUDE -I$HDF5_INCLUDE $OMPFLAGS $MPIFLAGS $LAPACKFLAGS $MAGMAFLAGS $NVTXFLAGS"
 SUBDIRS="cluster libnptm sphere testing trapping"
 
 # Generate the output
diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index b1f0c0b79043b5e17b66ff98233679ac7788eee2..11d710569030a2a30c465daa30e1f7f87d195c65 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -83,10 +83,6 @@
 #include "../include/file_io.h"
 #endif
 
-#ifndef INCLUDE_UTILS_H_
-#include "../include/utils.h"
-#endif
-
 using namespace std;
 
 // I would like to put it all in a struct, but then I'd have to write a constructor for it, due to members defined as references, creating a worse nightmare than the one I'd like to simplify...
@@ -374,7 +370,7 @@ void cluster(const string& config_file, const string& data_file, const string& o
       int myjxi488startoffset = 0;
       int myMPIstride = ompnumthreads;
       int myMPIblock = ompnumthreads;
-      // Define here shared arrays of virtual ascii and binary files, so that thread 0 will be able to access the all later
+      // Define here shared arrays of virtual ascii and binary files, so that thread 0 will be able to access them all later
       VirtualAsciiFile **p_outarray = NULL;
       VirtualBinaryFile **vtppoanarray = NULL;
 
@@ -766,8 +762,41 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
   interval_start = chrono::high_resolution_clock::now();
 #ifdef USE_NVTX
   nvtxRangePush("Calculate inverted matrix");
+#endif
+#ifdef DEBUG_AM
+  /* now, before cms, output am to p_outam0 */
+  VirtualAsciiFile *outam0 = new VirtualAsciiFile();
+  string outam0_name = output_path + "/c_AM0_JXI" + to_string(jxi488) + ".txt";
+  sprintf(virtual_line, " AM matrix before CMS\n");
+  outam0->append_line(virtual_line);
+  sprintf(virtual_line, " I1+1   I2+1    Real    Imag\n");
+  outam0->append_line(virtual_line);
+  for (int iam=0; iam<ndit; iam++) {
+    for (int jam=0; jam<ndit; jam++) {
+      sprintf(virtual_line, " %d   %d   %17.8lE   %17.8lE\n", iam, jam, real(cid->am[iam][jam]), imag(cid->am[iam][jam]));
+      outam0->append_line(virtual_line);
+    }
+  }
+  outam0->write_to_disk(outam0_name);
+  delete outam0;
 #endif
   cms(cid->am, cid->c1, cid->c1ao, cid->c4, cid->c6);
+#ifdef DEBUG_AM
+  VirtualAsciiFile *outam1 = new VirtualAsciiFile();
+  string outam1_name = output_path + "/c_AM1_JXI" + to_string(jxi488) + ".txt";
+  sprintf(virtual_line, " AM matrix after CMS before LUCIN\n");
+  outam1->append_line(virtual_line);
+  sprintf(virtual_line, " I1+1   I2+1    Real    Imag\n");
+  outam1->append_line(virtual_line);
+  for (int iam=0; iam<ndit; iam++) {
+    for (int jam=0; jam<ndit; jam++) {
+      sprintf(virtual_line, " %d   %d   %17.8lE   %17.8lE\n", iam, jam, real(cid->am[iam][jam]), imag(cid->am[iam][jam]));
+      outam1->append_line(virtual_line);
+    }
+  }
+  outam1->write_to_disk(outam1_name);
+  delete outam1;
+#endif
 #ifdef USE_NVTX
   nvtxRangePop();
 #endif
@@ -775,15 +804,27 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
   elapsed = interval_end - interval_start;
   message = "INFO: matrix calculation for scale " + to_string(jxi488) + " took " + to_string(elapsed.count()) + "s.\n";
   logger->log(message);
-  if (jxi488 == 1) {
-    string am_file_name = "c_am_jxi_1.txt";
-    write_dcomplex_matrix(am_file_name, cid->am, 624, 624);
-  }
   interval_start = chrono::high_resolution_clock::now();
 #ifdef USE_NVTX
   nvtxRangePush("Invert the matrix");
 #endif
   invert_matrix(cid->am, ndit, jer, mxndm, cid->proc_device);
+#ifdef DEBUG_AM
+  VirtualAsciiFile *outam2 = new VirtualAsciiFile();
+  string outam2_name = output_path + "/c_AM2_JXI" + to_string(jxi488) + ".txt";
+  sprintf(virtual_line, " AM matrix after LUCIN before ZTM\n");
+  outam2->append_line(virtual_line);
+  sprintf(virtual_line, " I1+1   I2+1    Real    Imag\n");
+  outam2->append_line(virtual_line);
+  for (int iam=0; iam<ndit; iam++) {
+    for (int jam=0; jam<ndit; jam++) {
+      sprintf(virtual_line, " %d   %d   %17.8lE   %17.8lE\n", iam, jam, real(cid->am[iam][jam]), imag(cid->am[iam][jam]));
+      outam2->append_line(virtual_line);
+    }
+  }
+  outam2->write_to_disk(outam2_name);
+  delete outam2;
+#endif
 #ifdef USE_NVTX
   nvtxRangePop();
 #endif
@@ -802,6 +843,22 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
   nvtxRangePush("Average calculation");
 #endif
   ztm(cid->am, cid->c1, cid->c1ao, cid->c4, cid->c6, cid->c9);
+#ifdef DEBUG_AM
+  VirtualAsciiFile *outam3 = new VirtualAsciiFile();
+  string outam3_name = output_path + "/c_AM3_JXI" + to_string(jxi488) + ".txt";
+  sprintf(virtual_line, " AM matrix after ZTM\n");
+  outam3->append_line(virtual_line);
+  sprintf(virtual_line, " I1+1   I2+1    Real    Imag\n");
+  outam3->append_line(virtual_line);
+  for (int iam=0; iam<ndit; iam++) {
+    for (int jam=0; jam<ndit; jam++) {
+      sprintf(virtual_line, " %d   %d   %17.8lE   %17.8lE\n", iam, jam, real(cid->am[iam][jam]), imag(cid->am[iam][jam]));
+      outam3->append_line(virtual_line);
+    }
+  }
+  outam3->write_to_disk(outam3_name);
+  delete outam3;
+#endif
   if (idfc >= 0) {
     if (jxi488 == jwtm) {
       int nlemt = 2 * cid->c4->nlem;