diff --git a/src/include/Commons.h b/src/include/Commons.h
index 1a1de770a39d1e4e35cbccfc7222e073860d26f9..c74427f1cf6ac45f4ccb19163f4a087c6a2d27da 100644
--- a/src/include/Commons.h
+++ b/src/include/Commons.h
@@ -58,6 +58,10 @@ class C1 {
 protected:
   //! \brief Maximum order of field expansion.
   int lm;
+  //! \brief Maximum order of internal field expansion.
+  int li;
+  //! \brief Maximum order of external field expansion.
+  int le;
   //! \brief Contiguous RMI vector.
   dcomplex *vec_rmi;
   //! \brief Contiguous REI vector.
@@ -70,6 +74,8 @@ protected:
 public:
   //! \brief Number of spheres.
   int nsph;
+  //! \brief NLEMT = 2 * LE * (LE + 2).
+  int nlemt;
   //! \brief NLMMT = 2 * LM * (LM + 2).
   int nlmmt;
   //! \brief Number of configurations
diff --git a/src/libnptm/Commons.cpp b/src/libnptm/Commons.cpp
index 0cecd8e615cbdc552561ac0c14cdf2cf631ad12f..8539e0b2210ac8fd067f0aef5fa3d2fc7b463a85 100644
--- a/src/libnptm/Commons.cpp
+++ b/src/libnptm/Commons.cpp
@@ -36,25 +36,26 @@
 
 C1::C1(GeometryConfiguration *gconf, ScattererConfiguration *sconf) {
   lm = gconf->l_max;
-  int li = gconf->li;
-  int le = gconf->le;
+  li = gconf->li;
+  le = gconf->le;
   if (lm == 0) {
     lm = (li > le) ? li : le;
   }
   nsph = gconf->number_of_spheres;
   nlmmt = 2 * (lm * (lm + 2));
+  nlemt = 2 * (le * (le + 2));
 
-  vec_rmi = new dcomplex[nsph * lm]();
-  vec_rei = new dcomplex[nsph * lm]();
-  rmi = new dcomplex*[lm];
-  rei = new dcomplex*[lm];
-  for (int ri = 0; ri < lm; ri++) {
+  vec_rmi = new dcomplex[nsph * li]();
+  vec_rei = new dcomplex[nsph * li]();
+  rmi = new dcomplex*[li];
+  rei = new dcomplex*[li];
+  for (int ri = 0; ri < li; ri++) {
     rmi[ri] = &(vec_rmi[nsph * ri]);
     rei[ri] = &(vec_rei[nsph * ri]);
   }
-  vec_w = new dcomplex[4 * nlmmt]();
-  w = new dcomplex*[nlmmt];
-  for (int wi = 0; wi < nlmmt; wi++) w[wi] = &(vec_w[4 * wi]);
+  vec_w = new dcomplex[4 * nlemt]();
+  w = new dcomplex*[nlemt];
+  for (int wi = 0; wi < nlemt; wi++) w[wi] = &(vec_w[4 * wi]);
   configurations = sconf->configurations;
   vint = new dcomplex[16]();
   vec_vints = new dcomplex[nsph * 16]();
@@ -108,28 +109,31 @@ C1::C1(GeometryConfiguration *gconf, ScattererConfiguration *sconf) {
 
 C1::C1(const C1& rhs) {
   nlmmt = rhs.nlmmt;
+  nlemt = rhs.nlemt;
   nsph = rhs.nsph;
   lm = rhs.lm;
+  li = rhs.li;
+  le = rhs.le;
 
-  vec_rmi = new dcomplex[nsph * lm]();
-  vec_rei = new dcomplex[nsph * lm]();
-  for (long li=0; li<(nsph*lm); li++) {
-    vec_rmi[li] = rhs.vec_rmi[li];
-    vec_rei[li] = rhs.vec_rei[li];
-  }
-  rmi = new dcomplex*[lm];
-  rei = new dcomplex*[lm];
-  for (int ri = 0; ri < lm; ri++) {
+  vec_rmi = new dcomplex[nsph * li]();
+  vec_rei = new dcomplex[nsph * li]();
+  for (long ili=0; ili < (nsph * li); ili++) {
+    vec_rmi[ili] = rhs.vec_rmi[ili];
+    vec_rei[ili] = rhs.vec_rei[ili];
+  }
+  rmi = new dcomplex*[li];
+  rei = new dcomplex*[li];
+  for (int ri = 0; ri < li; ri++) {
     rmi[ri] = &(vec_rmi[nsph * ri]);
     rei[ri] = &(vec_rei[nsph * ri]);
     /*! The contents were already copied via vec_rmi and vec_rei */
   }
-  vec_w = new dcomplex[4 * nlmmt]();
-  for (long li=0; li<(4*nlmmt); li++) {
-    vec_w[li] = rhs.vec_w[li];
+  vec_w = new dcomplex[4 * nlemt]();
+  for (long ili=0; ili < (4*nlemt); ili++) {
+    vec_w[ili] = rhs.vec_w[ili];
   }    
-  w = new dcomplex*[nlmmt];
-  for (int wi = 0; wi < nlmmt; wi++) {
+  w = new dcomplex*[nlemt];
+  for (int wi = 0; wi < nlemt; wi++) {
     w[wi] = &(vec_w[4 * wi]);
     /*! The contents were already copied via vec_w */
   }
@@ -198,23 +202,26 @@ C1::C1(const C1& rhs) {
 #ifdef MPI_VERSION
 C1::C1(const mixMPI *mpidata) {
   MPI_Bcast(&nlmmt, 1, MPI_INT, 0, MPI_COMM_WORLD);
+  MPI_Bcast(&nlemt, 1, MPI_INT, 0, MPI_COMM_WORLD);
   MPI_Bcast(&nsph, 1, MPI_INT, 0, MPI_COMM_WORLD);
   MPI_Bcast(&lm, 1, MPI_INT, 0, MPI_COMM_WORLD);
-  vec_rmi = new dcomplex[nsph * lm]();
-  vec_rei = new dcomplex[nsph * lm]();
-  MPI_Bcast(vec_rmi, nsph*lm, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
-  MPI_Bcast(vec_rei, nsph*lm, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
-  rmi = new dcomplex*[lm];
-  rei = new dcomplex*[lm];
-  for (int ri = 0; ri < lm; ri++) {
+  MPI_Bcast(&li, 1, MPI_INT, 0, MPI_COMM_WORLD);
+  MPI_Bcast(&le, 1, MPI_INT, 0, MPI_COMM_WORLD);
+  vec_rmi = new dcomplex[nsph * li]();
+  vec_rei = new dcomplex[nsph * li]();
+  MPI_Bcast(vec_rmi, nsph*li, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+  MPI_Bcast(vec_rei, nsph*li, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+  rmi = new dcomplex*[li];
+  rei = new dcomplex*[li];
+  for (int ri = 0; ri < li; ri++) {
     rmi[ri] = &(vec_rmi[nsph * ri]);
     rei[ri] = &(vec_rei[nsph * ri]);
     /*! The contents were copied already via vec_rmi and vec_rei */
   }
-  vec_w = new dcomplex[4 * nlmmt]();
-  MPI_Bcast(vec_w, 4*nlmmt, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
-  w = new dcomplex*[nlmmt];
-  for (int wi = 0; wi < nlmmt; wi++) {
+  vec_w = new dcomplex[4 * nlemt]();
+  MPI_Bcast(vec_w, 4*nlemt, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+  w = new dcomplex*[nlemt];
+  for (int wi = 0; wi < nlemt; wi++) {
     w[wi] = &(vec_w[4 * wi]);
     /*! The contents were copied already via vec_w */
   }
@@ -279,11 +286,14 @@ C1::C1(const mixMPI *mpidata) {
 
 void C1::mpibcast(const mixMPI *mpidata) {
   MPI_Bcast(&nlmmt, 1, MPI_INT, 0, MPI_COMM_WORLD);
+  MPI_Bcast(&nlemt, 1, MPI_INT, 0, MPI_COMM_WORLD);
   MPI_Bcast(&nsph, 1, MPI_INT, 0, MPI_COMM_WORLD);
   MPI_Bcast(&lm, 1, MPI_INT, 0, MPI_COMM_WORLD);
-  MPI_Bcast(vec_rmi, nsph*lm, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
-  MPI_Bcast(vec_rei, nsph*lm, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
-  MPI_Bcast(vec_w, 4*nlmmt, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+  MPI_Bcast(&li, 1, MPI_INT, 0, MPI_COMM_WORLD);
+  MPI_Bcast(&le, 1, MPI_INT, 0, MPI_COMM_WORLD);
+  MPI_Bcast(vec_rmi, nsph*li, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+  MPI_Bcast(vec_rei, nsph*li, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+  MPI_Bcast(vec_w, 4*nlemt, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
   MPI_Bcast(&configurations, 1, MPI_INT, 0, MPI_COMM_WORLD);
   MPI_Bcast(vint, 16, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
   MPI_Bcast(vec_vints, nsph*16, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);