Adapt dc0 vector in ParticlDescriptor for multi-layer aggregate case

src/include/Commons.h

@@ -80,12 +80,12 @@ protected:
   int _nsph;
   //! \brief Maximum internal field expansion order.
   int _li;
+  //! \brief Maximum number of layers in known sphere types.
+  int _max_layers;
   //! \brief Number of different sphere types.
   int _num_configurations;
   //! \brief Total number of layers from all sphere types.
   int _num_layers;
-  //! \brief Space for different sphere types.
-  int _nl;
   //! \brief NHSPO = 2 * MAX(NPNT,NPNTTS) - 1
   int _nhspo;
   //! \brief Number of points for numerical integration in layered spheres.
@@ -184,12 +184,12 @@ public:
   const int &nsph = _nsph;
   //! \brief Read-only view of maximum internal field expansion order.
   const int &li = _li;
+  //! \brief Read-only view of the maximum number of layers in types.
+  const int &max_layers = _max_layers;
   //! \brief Read-only view of number of different sphere types.
   const int &num_configurations = _num_configurations;
   //! \brief Read-only view of total number of layers from all sphere types.
   const int &num_layers = _num_layers;
-  //! \brief Read-only view on the space for different sphere configurations.
-  const int &nl = _nl;
   //! \brief Read-only view of NHSPO.
   const int &nhspo = _nhspo;
   //! \brief Read-only view on number of points for numerical integration in layered spheres.

src/libnptm/Commons.cpp

@@ -84,7 +84,12 @@ ParticleDescriptor::ParticleDescriptor(GeometryConfiguration *gconf, ScattererCo
   gcs = 0.0;
   _num_configurations = sconf->configurations;
   _num_layers = (sconf->use_external_sphere) ? 1 : 0;
-  for (int nli = 0; nli < num_configurations; nli++) _num_layers += sconf->get_nshl(nli);
+  _max_layers = 1;
+  for (int nli = 0; nli < num_configurations; nli++) {
+    int nl = sconf->get_nshl(nli);
+    _num_layers += nl;
+    if (nl >  _max_layers) _max_layers = nl;
+  }
 
   vec_rmi = new dcomplex[_li * _nsph]();
   rmi = new dcomplex*[_li];
@@ -131,12 +136,12 @@ ParticleDescriptor::ParticleDescriptor(GeometryConfiguration *gconf, ScattererCo
   _npntts = gconf->npntts;
   int max_n = (npnt > npntts) ? npnt : npntts;
   _nhspo = 2 * max_n - 1;
-  _nl = sconf->configurations;
-  if (_nsph == 1 && _nl == 1) _nl = 5;
+  // _nl = sconf->configurations;
+  // if (_nsph == 1 && _nl == 1) _nl = 5;
   ris = new dcomplex[_nhspo]();
   dlri = new dcomplex[_nhspo]();
   vkt = new dcomplex[_nsph]();
-  dc0 = new dcomplex[_nl]();
+  dc0 = new dcomplex[_max_layers + 1]();
   vsz = new double[_nsph]();
   
   // >>> NEEDED BY SPHERE AND CLUSTER <<<
@@ -219,6 +224,7 @@ ParticleDescriptor::ParticleDescriptor(const ParticleDescriptor &rhs) {
   _ndit = rhs._ndit;
   _ndm = rhs._ndm;
   gcs = rhs.gcs;
+  _max_layers = rhs._max_layers;
   _num_configurations = rhs._num_configurations;
   _num_layers = rhs._num_layers;
 
@@ -271,7 +277,7 @@ ParticleDescriptor::ParticleDescriptor(const ParticleDescriptor &rhs) {
   _npnt = rhs._npnt;
   _npntts = rhs._npntts;
   _nhspo = rhs._nhspo;
-  _nl = rhs._nl;
+  _max_layers = rhs._max_layers;
   ris = new dcomplex[_nhspo]();
   dlri = new dcomplex[_nhspo]();
   for (int ri = 0; ri < _nhspo; ri++) {
@@ -284,8 +290,8 @@ ParticleDescriptor::ParticleDescriptor(const ParticleDescriptor &rhs) {
     vkt[vi] = rhs.vkt[vi];
     vsz[vi] = rhs.vsz[vi];
   }
-  dc0 = new dcomplex[_nl]();
-  for (int di = 0; di < _nl; di++) {
+  dc0 = new dcomplex[_max_layers]();
+  for (int di = 0; di < _max_layers; di++) {
     dc0[di] = rhs.dc0[di];
   }
   // >>> NEEDED BY SPHERE AND CLUSTER <<<
@@ -416,7 +422,7 @@ ParticleDescriptor::ParticleDescriptor(const mixMPI *mpidata) {
   MPI_Bcast(&_npnt, 1, MPI_INT, 0, MPI_COMM_WORLD);
   MPI_Bcast(&_npntts, 1, MPI_INT, 0, MPI_COMM_WORLD);
   MPI_Bcast(&_nhspo, 1, MPI_INT, 0, MPI_COMM_WORLD);
-  MPI_Bcast(&_nl, 1, MPI_INT, 0, MPI_COMM_WORLD);
+  MPI_Bcast(&_max_layers, 1, MPI_INT, 0, MPI_COMM_WORLD);
   ris = new dcomplex[_nhspo];
   MPI_Bcast(ris, _nhspo, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
   dlri = new dcomplex[_nhspo];
@@ -425,8 +431,8 @@ ParticleDescriptor::ParticleDescriptor(const mixMPI *mpidata) {
   MPI_Bcast(vkt, _nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
   vsz = new double[_nsph];
   MPI_Bcast(vsz, _nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
-  dc0 = new dcomplex[_nl];
-  MPI_Bcast(dc0, _nl, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+  dc0 = new dcomplex[_max_layers];
+  MPI_Bcast(dc0, _max_layers, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
 }
 
 void ParticleDescriptor::mpibcast(const mixMPI *mpidata) {
@@ -466,12 +472,12 @@ void ParticleDescriptor::mpibcast(const mixMPI *mpidata) {
   MPI_Bcast(&_npnt, 1, MPI_INT, 0, MPI_COMM_WORLD);
   MPI_Bcast(&_npntts, 1, MPI_INT, 0, MPI_COMM_WORLD);
   MPI_Bcast(&_nhspo, 1, MPI_INT, 0, MPI_COMM_WORLD);
-  MPI_Bcast(&_nl, 1, MPI_INT, 0, MPI_COMM_WORLD);
+  MPI_Bcast(&_max_layers, 1, MPI_INT, 0, MPI_COMM_WORLD);
   MPI_Bcast(ris, _nhspo, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
   MPI_Bcast(dlri, _nhspo, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
   MPI_Bcast(vkt, _nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
   MPI_Bcast(vsz, _nsph, MPI_DOUBLE, 0, MPI_COMM_WORLD);
-  MPI_Bcast(dc0, _nl, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
+  MPI_Bcast(dc0, _max_layers, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);
   // >>> NEEDED BY SPHERE AND CLUSTER <<< //
   if (_class_type == SPHERE_TYPE || _class_type == CLUSTER_TYPE) {
     MPI_Bcast(vec_sas, 4 * _nsph, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);