diff --git a/src/cluster/cluster.cpp b/src/cluster/cluster.cpp
index 11d710569030a2a30c465daa30e1f7f87d195c65..e7da6b93ad2f1bd4fb51d352bfc6195beecf3d28 100644
--- a/src/cluster/cluster.cpp
+++ b/src/cluster/cluster.cpp
@@ -83,6 +83,10 @@
 #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...
@@ -771,12 +775,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
   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);
-    }
-  }
+  write_dcomplex_matrix(outam0, cid->am, ndit, ndit);
   outam0->write_to_disk(outam0_name);
   delete outam0;
 #endif
@@ -788,12 +787,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
   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);
-    }
-  }
+  write_dcomplex_matrix(outam1, cid->am, ndit, ndit, " %5d %5d (%17.8lE,%17.8lE)\n", 1);
   outam1->write_to_disk(outam1_name);
   delete outam1;
 #endif
@@ -816,12 +810,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
   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);
-    }
-  }
+  write_dcomplex_matrix(outam2, cid->am, ndit, ndit);
   outam2->write_to_disk(outam2_name);
   delete outam2;
 #endif
@@ -850,12 +839,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf
   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);
-    }
-  }
+  write_dcomplex_matrix(outam3, cid->am, ndit, ndit);
   outam3->write_to_disk(outam3_name);
   delete outam3;
 #endif
diff --git a/src/include/utils.h b/src/include/utils.h
index 8a3c372213dc6a063e48adfcacf28324d5d52bd8..82779beb33fa14026b8a20c3cbf9f47a01a0c65c 100644
--- a/src/include/utils.h
+++ b/src/include/utils.h
@@ -24,11 +24,17 @@
 
 /*! \brief Write a double complex matrix to a text file.
  *
- * \param file_name: `const string&` Name of the file to be written.
+ * \param af: `VirtualAsciiFile *` Pointer to an existing VirtualAsciiFile.
  * \param mat: `dcomplex **` Pointer to the matrix.
- * \param rows: `np_int` Number of rows in the matrix.
- * \param columns: `np_int` Number of columns in the matrix.
+ * \param rows: `int` Number of rows in the matrix.
+ * \param columns: `int` Number of columns in the matrix.
+ * \param format: `const string&` Format of the line (default is \" %5d %5d (%17.8lE,%17.8lE)\n\")
+ * \param first_index: `int` Index of the first element (default is 1, i.e. base 1 FORTRAN array notation)
  */
-int write_dcomplex_matrix(const std::string& file_name, dcomplex **mat, np_int rows, np_int columns);
+int write_dcomplex_matrix(
+			  VirtualAsciiFile *af, dcomplex **mat, int rows,
+			  int columns, const std::string& format=" %5d %5d (%17.8lE,%17.8lE)\n",
+			  int first_index=1
+);
 
 #endif
diff --git a/src/libnptm/utils.cpp b/src/libnptm/utils.cpp
index d7514d6dfbfc7e22093d11da5804d1b79d800dc3..f989f07565b33e18d93ed6390675468a63b41fb3 100644
--- a/src/libnptm/utils.cpp
+++ b/src/libnptm/utils.cpp
@@ -40,33 +40,43 @@
  * \brief Definition of auxiliary code utilities.
  */
 
-#include <cstdio>
+#include <hdf5.h>
 #ifndef INCLUDE_TYPES_H_
 #include "../include/types.h"
 #endif
+
+#ifndef INCLUDE_ERRORS_H_
+#include "../include/errors.h"
+#endif
+
+#ifndef INCLUDE_LIST_H_
+#include "../include/List.h"
+#endif
+
+#ifndef INCLUDE_FILE_IO_H_
+#include "../include/file_io.h"
+#endif
+
 #ifndef INCLUDE_UTILS_H_
 #include "../include/utils.h"
 #endif
 
 using namespace std;
 
-int write_dcomplex_matrix(const std::string& file_name, dcomplex **mat, np_int rows, np_int columns) {
+int write_dcomplex_matrix(
+			  VirtualAsciiFile *af, dcomplex **mat, int rows, int columns,
+			  const std::string& format, int first_index
+) {
   int result = 0;
-  FILE* output = fopen(file_name.c_str(), "w");
-  if (output) {
-    string str_format = "%5d%5d (%13.5lE,%13.5lE)\n";
-    if (sizeof(np_int) > sizeof(int)) {
-      str_format = "%5ld%5ld (%13.5lE,%13.5lE)\n";
-    }
-    for (np_int i = 0; i < rows; i++) {
-      for (np_int j = 0; j < columns; j++) {
-	fprintf(output, str_format.c_str(), i, j, real(mat[i][j]), imag(mat[i][j]));
-      }
+  char virtual_line[256];
+  for (int i=0; i < rows; i++) {
+    for (int j = 0; j < columns; j++) {
+      sprintf(
+	      virtual_line, format.c_str(), i + first_index, j + first_index,
+	      real(mat[i][j]), imag(mat[i][j])
+      );
+      af->append_line(virtual_line);
     }
-    fclose(output);
-  } else {
-    // Could not open the output file.
-    result = 1;
   }
   return result;
 }