diff --git a/src/libnptm/Commons.cpp b/src/libnptm/Commons.cpp
index 81cb853bfcb1f757cab084a1f587a620c9b14f0d..d19e62b58843443c3c1540d43f353b05e55bd025 100644
--- a/src/libnptm/Commons.cpp
+++ b/src/libnptm/Commons.cpp
@@ -19,234 +19,248 @@
 using namespace std;
 
 C1::C1(int ns, int l_max, int *_nshl, int *_iog) {
-	nlmmt = 2 * (l_max * (l_max + 2));
-	nsph = ns;
-	lm = l_max;
-
-	rmi = new complex<double>*[lm];
-	rei = new complex<double>*[lm];
-	for (int ri = 0; ri < lm; ri++) {
-		rmi[ri] = new complex<double>[nsph]();
-		rei[ri] = new complex<double>[nsph]();
-	}
-	w = new complex<double>*[nlmmt];
-	for (int wi = 0; wi < nlmmt; wi++) w[wi] = new complex<double>[4]();
-	vints = new complex<double>*[nsph];
-	rc = new double*[nsph];
-	nshl = new int[nsph]();
-	iog = new int[nsph]();
-	for (int vi = 0; vi < nsph; vi++) {
-		rc[vi] = new double[_nshl[vi]]();
-		vints[vi] = new complex<double>[16]();
-		nshl[vi] = _nshl[vi];
-		iog[vi] = _iog[vi];
-	}
-	fsas = new complex<double>[nsph]();
-	sscs = new double[nsph]();
-	sexs = new double[nsph]();
-	sabs = new double[nsph]();
-	sqscs = new double[nsph]();
-	sqexs = new double[nsph]();
-	sqabs = new double[nsph]();
-	gcsv = new double[nsph]();
-	rxx = new double[nsph]();
-	ryy = new double[nsph]();
-	rzz = new double[nsph]();
-	ros = new double[nsph]();
-
-	sas = new complex<double>**[nsph];
-	for (int si = 0; si < nsph; si++) {
-		sas[si] = new complex<double>*[2];
-		sas[si][0] = new complex<double>[2]();
-		sas[si][1] = new complex<double>[2]();
-	}
+  nlmmt = 2 * (l_max * (l_max + 2));
+  nsph = ns;
+  lm = l_max;
+
+  rmi = new complex<double>*[lm];
+  rei = new complex<double>*[lm];
+  for (int ri = 0; ri < lm; ri++) {
+    rmi[ri] = new complex<double>[nsph]();
+    rei[ri] = new complex<double>[nsph]();
+  }
+  w = new complex<double>*[nlmmt];
+  for (int wi = 0; wi < nlmmt; wi++) w[wi] = new complex<double>[4]();
+  int configurations = 0;
+  for (int ci = 1; ci <= nsph; ci++) {
+    if (_iog[ci - 1] >= ci) configurations++;
+  }
+  vints = new complex<double>*[nsph];
+  rc = new double*[configurations];
+  nshl = new int[configurations]();
+  iog = new int[nsph]();
+  int conf_index = 0;
+  for (int vi = 0; vi < nsph; vi++) {
+    vints[vi] = new complex<double>[16]();
+    iog[vi] = _iog[vi];
+    if (iog[vi] >= vi + 1) {
+      nshl[conf_index] = _nshl[conf_index];
+      rc[conf_index] = new double[_nshl[conf_index]]();
+      conf_index++;
+    }
+  }
+  fsas = new complex<double>[nsph]();
+  sscs = new double[nsph]();
+  sexs = new double[nsph]();
+  sabs = new double[nsph]();
+  sqscs = new double[nsph]();
+  sqexs = new double[nsph]();
+  sqabs = new double[nsph]();
+  gcsv = new double[nsph]();
+  rxx = new double[nsph]();
+  ryy = new double[nsph]();
+  rzz = new double[nsph]();
+  ros = new double[nsph]();
+
+  sas = new complex<double>**[nsph];
+  for (int si = 0; si < nsph; si++) {
+    sas[si] = new complex<double>*[2];
+    sas[si][0] = new complex<double>[2]();
+    sas[si][1] = new complex<double>[2]();
+  }
 }
 
 C1::~C1() {
-	delete[] rmi;
-	delete[] rei;
-	for (int wi = nlmmt - 1; wi > -1; wi--) delete[] w[wi];
-	for (int vi = nsph - 1; vi > - 1; vi--) {
-		delete[] rc[vi];
-		delete[] vints[vi];
-	}
-	delete[] rc;
-	delete[] vints;
-	for (int si = nsph - 1; si > -1; si--) {
-		delete[] sas[si][1];
-		delete[] sas[si][0];
-		delete[] sas[si];
-	}
-	delete[] sas;
-	delete[] fsas;
-	delete[] sscs;
-	delete[] sexs;
-	delete[] sabs;
-	delete[] sqscs;
-	delete[] sqexs;
-	delete[] sqabs;
-	delete[] gcsv;
-	delete[] rxx;
-	delete[] ryy;
-	delete[] rzz;
-	delete[] ros;
-	delete[] iog;
-	delete[] nshl;
+  delete[] rmi;
+  delete[] rei;
+  for (int wi = nlmmt - 1; wi > -1; wi--) delete[] w[wi];
+  int conf_index = 0;
+  for (int ci = 1; ci <= nsph; ci++) {
+    if (iog[ci] >= ci) {
+      delete[] rc[conf_index];
+      conf_index++;
+    }
+  }
+  delete[] rc;
+  for (int vi = nsph - 1; vi > - 1; vi--) {
+    delete[] vints[vi];
+  }
+  delete[] vints;
+  for (int si = nsph - 1; si > -1; si--) {
+    delete[] sas[si][1];
+    delete[] sas[si][0];
+    delete[] sas[si];
+  }
+  delete[] sas;
+  delete[] fsas;
+  delete[] sscs;
+  delete[] sexs;
+  delete[] sabs;
+  delete[] sqscs;
+  delete[] sqexs;
+  delete[] sqabs;
+  delete[] gcsv;
+  delete[] rxx;
+  delete[] ryy;
+  delete[] rzz;
+  delete[] ros;
+  delete[] iog;
+  delete[] nshl;
 }
 
 C1_AddOns::C1_AddOns(C4 *c4) {
-	nsph = c4->nsph;
-	lmpo = c4->lmpo;
-	nlemt = 2 * c4->nlem;
-	vh = new complex<double>[(nsph * nsph - 1) * c4->litpo]();
-	vj0 = new complex<double>[nsph * c4->lmtpo]();
-	vj = new complex<double>[1](); // QUESTION: is 1 really enough for a general case?
-	vyhj = new complex<double>[(nsph * nsph - 1) * c4->litpos]();
-	vyj0 = new complex<double>[nsph * c4->lmtpos]();
-	am0m = new complex<double>*[nlemt];
-	for (int ai = 0; ai < nlemt; ai++) {
-		am0m[ai] = new complex<double>[nlemt]();
-	}
-	vint = new complex<double>[16](); // QUESTION: is dimension 16 generally fixed?
-	vintm = new complex<double>[16]();
-	vintt = new complex<double>[16]();
-	vints = new complex<double>*[nsph];
-	for (int vi = 0; vi < nsph; vi++) vints[vi] = new complex<double>[16]();
-	fsac = new complex<double>*[2];
-	sac = new complex<double>*[2];
-	fsacm = new complex<double>*[2];
-	for (int fi = 0; fi < 2; fi++) {
-		fsac[fi] = new complex<double>[2]();
-		sac[fi] = new complex<double>[2]();
-		fsacm[fi] = new complex<double>[2]();
-	}
-	scscp = new complex<double>[2]();
-	ecscp = new complex<double>[2]();
-	scscpm = new complex<double>[2]();
-	ecscpm = new complex<double>[2]();
-	allocate_vectors(c4);
-	sscs = new double[nsph]();
-	ecscm = new double[2]();
-	scscm = new double[2]();
-	scsc = new double[2]();
-	ecsc = new double[2]();
+  nsph = c4->nsph;
+  lmpo = c4->lmpo;
+  nlemt = 2 * c4->nlem;
+  vh = new complex<double>[(nsph * nsph - 1) * c4->litpo]();
+  vj0 = new complex<double>[nsph * c4->lmtpo]();
+  vj = new complex<double>[1](); // QUESTION: is 1 really enough for a general case?
+  vyhj = new complex<double>[(nsph * nsph - 1) * c4->litpos]();
+  vyj0 = new complex<double>[nsph * c4->lmtpos]();
+  am0m = new complex<double>*[nlemt];
+  for (int ai = 0; ai < nlemt; ai++) {
+    am0m[ai] = new complex<double>[nlemt]();
+  }
+  vint = new complex<double>[16](); // QUESTION: is dimension 16 generally fixed?
+  vintm = new complex<double>[16]();
+  vintt = new complex<double>[16]();
+  vints = new complex<double>*[nsph];
+  for (int vi = 0; vi < nsph; vi++) vints[vi] = new complex<double>[16]();
+  fsac = new complex<double>*[2];
+  sac = new complex<double>*[2];
+  fsacm = new complex<double>*[2];
+  for (int fi = 0; fi < 2; fi++) {
+    fsac[fi] = new complex<double>[2]();
+    sac[fi] = new complex<double>[2]();
+    fsacm[fi] = new complex<double>[2]();
+  }
+  scscp = new complex<double>[2]();
+  ecscp = new complex<double>[2]();
+  scscpm = new complex<double>[2]();
+  ecscpm = new complex<double>[2]();
+  allocate_vectors(c4);
+  sscs = new double[nsph]();
+  ecscm = new double[2]();
+  scscm = new double[2]();
+  scsc = new double[2]();
+  ecsc = new double[2]();
 }
 
 C1_AddOns::~C1_AddOns() {
-	delete[] sscs;
-	delete[] vyj0;
-	delete[] vyhj;
-	for (int ii = lmpo - 1; ii > -1; ii--) delete[] ind3j[ii];
-	delete[] ind3j;
-	delete[] v3j0;
-	delete[] vh;
-	delete[] vj0;
-	for (int ai = nlemt - 1; ai > -1; ai--) {
-		delete[] am0m[ai];
-	}
-	delete am0m;
-	delete[] vint;
-	delete[] vintm;
-	delete[] vintt;
-	for (int vi = nsph - 1; vi > -1; vi--) delete[] vints[vi];
-	delete[] vints;
-	for (int fi = 1; fi > -1; fi--) {
-		delete[] fsac[fi];
-		delete[] sac[fi];
-		delete[] fsacm[fi];
-	}
-	delete[] fsac;
-	delete[] sac;
-	delete[] fsacm;
-	delete[] scscp;
-	delete[] ecscp;
-	delete[] scscpm;
-	delete[] ecscpm;
-	delete[] ecscm;
-	delete[] scscm;
-	delete[] scsc;
-	delete[] ecsc;
+  delete[] sscs;
+  delete[] vyj0;
+  delete[] vyhj;
+  for (int ii = lmpo - 1; ii > -1; ii--) delete[] ind3j[ii];
+  delete[] ind3j;
+  delete[] v3j0;
+  delete[] vh;
+  delete[] vj0;
+  for (int ai = nlemt - 1; ai > -1; ai--) {
+    delete[] am0m[ai];
+  }
+  delete am0m;
+  delete[] vint;
+  delete[] vintm;
+  delete[] vintt;
+  for (int vi = nsph - 1; vi > -1; vi--) delete[] vints[vi];
+  delete[] vints;
+  for (int fi = 1; fi > -1; fi--) {
+    delete[] fsac[fi];
+    delete[] sac[fi];
+    delete[] fsacm[fi];
+  }
+  delete[] fsac;
+  delete[] sac;
+  delete[] fsacm;
+  delete[] scscp;
+  delete[] ecscp;
+  delete[] scscpm;
+  delete[] ecscpm;
+  delete[] ecscm;
+  delete[] scscm;
+  delete[] scsc;
+  delete[] ecsc;
 }
 
 void C1_AddOns::allocate_vectors(C4 *c4) {
-	// This calculates the size of v3j0
-	int lm = (c4->li > c4->le) ? c4->li : c4->le;
-	const int nv3j = c4->nv3j;
-	v3j0 = new double[nv3j]();
-	ind3j = new int*[lmpo];
-	for (int ii = 0; ii < lmpo; ii++) ind3j[ii] = new int[c4->lm]();
-	// This calculates the size of vyhj
-	int ivy = (nsph * nsph - 1) * c4->litpos;
-	vyhj = new complex<double>[ivy]();
-	// This calculates the size of vyj0
-	ivy = c4->lmtpos * c4->nsph;
-	vyj0 = new complex<double>[ivy]();
+  // This calculates the size of v3j0
+  int lm = (c4->li > c4->le) ? c4->li : c4->le;
+  const int nv3j = c4->nv3j;
+  v3j0 = new double[nv3j]();
+  ind3j = new int*[lmpo];
+  for (int ii = 0; ii < lmpo; ii++) ind3j[ii] = new int[c4->lm]();
+  // This calculates the size of vyhj
+  int ivy = (nsph * nsph - 1) * c4->litpos;
+  vyhj = new complex<double>[ivy]();
+  // This calculates the size of vyj0
+  ivy = c4->lmtpos * c4->nsph;
+  vyj0 = new complex<double>[ivy]();
 }
 
 C2::C2(int ns, int nl, int npnt, int npntts) {
-	nsph = ns;
-	int max_n = (npnt > npntts) ? npnt : npntts;
-	nhspo = 2 * max_n - 1;
-	ris = new complex<double>[nhspo]();
-	dlri = new complex<double>[nhspo]();
-	vkt = new complex<double>[nsph]();
-	dc0 = new complex<double>[nl]();
-	vsz = new double[nsph]();
+  nsph = ns;
+  int max_n = (npnt > npntts) ? npnt : npntts;
+  nhspo = 2 * max_n - 1;
+  ris = new complex<double>[nhspo]();
+  dlri = new complex<double>[nhspo]();
+  vkt = new complex<double>[nsph]();
+  dc0 = new complex<double>[nl]();
+  vsz = new double[nsph]();
 }
 
 C2::~C2() {
-	delete[] ris;
-	delete[] dlri;
-	delete[] vkt;
-	delete[] dc0;
-	delete[] vsz;
+  delete[] ris;
+  delete[] dlri;
+  delete[] vkt;
+  delete[] dc0;
+  delete[] vsz;
 }
 
 C3::C3() {
-	tsas = new complex<double>*[2];
-	tsas[0] = new complex<double>[2];
-	tsas[1] = new complex<double>[2];
-	tfsas = complex<double>(0.0, 0.0);
-	gcs = 0.0;
-	scs = 0.0;
-	ecs = 0.0;
-	acs = 0.0;
+  tsas = new complex<double>*[2];
+  tsas[0] = new complex<double>[2];
+  tsas[1] = new complex<double>[2];
+  tfsas = complex<double>(0.0, 0.0);
+  gcs = 0.0;
+  scs = 0.0;
+  ecs = 0.0;
+  acs = 0.0;
 }
 
 C3::~C3() {
-	delete[] tsas[1];
-	delete[] tsas[0];
-	delete[] tsas;
+  delete[] tsas[1];
+  delete[] tsas[0];
+  delete[] tsas;
 }
 
 C6::C6(int lmtpo) {
-	rac3j = new double[lmtpo]();
+  rac3j = new double[lmtpo]();
 }
 
 C6::~C6() {
-	delete[] rac3j;
+  delete[] rac3j;
 }
 
 C9::C9(int ndi, int nlem, int ndit, int nlemt) {
-	gis = new complex<double>*[ndi];
-	gls = new complex<double>*[ndi];
-	for (int gi = 0; gi < ndi; gi++) {
-		gis[gi] = new complex<double>[nlem]();
-		gls[gi] = new complex<double>[nlem]();
-	}
-	sam = new complex<double>*[ndit];
-	for (int si = 0; si < ndit; si++) sam[si] = new complex<double>[nlemt]();
-	gis_size_0 = ndi;
-	sam_size_0 = ndit;
+  gis = new complex<double>*[ndi];
+  gls = new complex<double>*[ndi];
+  for (int gi = 0; gi < ndi; gi++) {
+    gis[gi] = new complex<double>[nlem]();
+    gls[gi] = new complex<double>[nlem]();
+  }
+  sam = new complex<double>*[ndit];
+  for (int si = 0; si < ndit; si++) sam[si] = new complex<double>[nlemt]();
+  gis_size_0 = ndi;
+  sam_size_0 = ndit;
 }
 
 C9::~C9() {
-	for (int gi = gis_size_0 - 1; gi > -1; gi--) {
-		delete[] gis[gi];
-		delete[] gls[gi];
-	}
-	delete[] gis;
-	delete[] gls;
-	for (int si = sam_size_0 - 1; si > -1; si--) delete[] sam[si];
-	delete[] sam;
+  for (int gi = gis_size_0 - 1; gi > -1; gi--) {
+    delete[] gis[gi];
+    delete[] gls[gi];
+  }
+  delete[] gis;
+  delete[] gls;
+  for (int si = sam_size_0 - 1; si > -1; si--) delete[] sam[si];
+  delete[] sam;
 }
diff --git a/src/libnptm/Configuration.cpp b/src/libnptm/Configuration.cpp
index e7ed703baf80b284252177a8732fe28930616c8c..04cdb5c03a0930affa172afb46d7db960166914b 100644
--- a/src/libnptm/Configuration.cpp
+++ b/src/libnptm/Configuration.cpp
@@ -1067,6 +1067,12 @@ void ScattererConfiguration::write_formatted(string file_name) {
 					       );
       break;
     }
+    // Clean memory
+    delete[] xi_vec;
+    delete[] pu_vec;
+    delete[] ev_vec;
+    delete[] wn_vec;
+    delete[] wl_vec;
   } else { // idfc < 0, Dielectric functions are at XIP and XI is scale for dimensions
     double pu, wn;
     xi_vec = scale_vec;