diff --git a/src/scripts/model_maker.py b/src/scripts/model_maker.py
index 9e6bb1cc6abd60b47beee5519c5096bb6db8f741..bd31e186678b3c14a8191df724e04e4be7ee2ab9 100755
--- a/src/scripts/model_maker.py
+++ b/src/scripts/model_maker.py
@@ -168,16 +168,17 @@ def load_model(model_file):
             print("ERROR: declared array of optical constants does not match configurations!")
             return (None, None)
         else:
-            sconf['dielec_id'] = [
-                [ 0 for j in range(max_layers)] for i in range(sconf['configurations'])
-            ]
-            for i in range(sconf['configurations']):
-                if (len(model['particle_settings']['dielec_id'][i]) != sconf['nshl'][i]):
-                    print("ERROR: Declared materials in type %d do not match number of layers!"%i)
-                    return (None, None)
-                else:
-                    for j in range(sconf['nshl'][i]):
-                        sconf['dielec_id'][i][j] = float(model['particle_settings']['dielec_id'][i][j])
+            if (sconf['idfc'] == 0):
+                sconf['dielec_id'] = [
+                    [ 0 for j in range(max_layers)] for i in range(sconf['configurations'])
+                ]
+                for i in range(sconf['configurations']):
+                    if (len(model['particle_settings']['dielec_id'][i]) != sconf['nshl'][i]):
+                        print("ERROR: Declared materials in type %d do not match number of layers!"%i)
+                        return (None, None)
+                    else:
+                        for j in range(sconf['nshl'][i]):
+                            sconf['dielec_id'][i][j] = float(model['particle_settings']['dielec_id'][i][j])
         if (model['radiation_settings']['scale_name'] == "XI"):
             sconf['insn'] = 1
             sconf['nxi'] = 1 + int((sconf['xi_end'] - sconf['xi_start']) / sconf['xi_step'])
@@ -203,7 +204,7 @@ def load_model(model_file):
                     for di in range(max_layers):
                         for dj in range(sconf['configurations']):
                             if (len(model['material_settings']['diel_const'][dj]) / 2 != sconf['nshl'][dj]):
-                                print("ERROR: dielectric constants for type %dj do not match number of layers!"%(dj + 1))
+                                print("ERROR: dielectric constants for type %d do not match number of layers!"%(dj + 1))
                                 return (None, None)
                             else:
                                 sconf['rdc0'][di][dj][0] = float(model['material_settings']['diel_const'][dj][2 * di])
@@ -289,7 +290,8 @@ def load_model(model_file):
         gconf['inpol'] = 0 if str_polar == "LINEAR" else 1
         gconf['npnt'] = int(model['geometry_settings']['npnt'])
         gconf['npntts'] = int(model['geometry_settings']['npntts'])
-        gconf['iavm'] = int(model['geometry_settings']['iavm'])
+        if (gconf['application'] != "SPHERE"):
+            gconf['iavm'] = int(model['geometry_settings']['iavm'])
         gconf['isam'] = int(model['geometry_settings']['isam'])
         gconf['jwtm'] = int(model['output_settings']['jwtm'])
         gconf['th'] = float(model['geometry_settings']['in_th_start'])
@@ -304,6 +306,9 @@ def load_model(model_file):
         gconf['phs'] = float(model['geometry_settings']['sc_ph_start'])
         gconf['phsstp'] = float(model['geometry_settings']['sc_ph_step'])
         gconf['phslst'] = float(model['geometry_settings']['sc_ph_end'])
+        gconf['vec_sph_x'] = [0.0 for i in range(gconf['nsph'])]
+        gconf['vec_sph_y'] = [0.0 for i in range(gconf['nsph'])]
+        gconf['vec_sph_z'] = [0.0 for i in range(gconf['nsph'])]
         if (gconf['application'] != "SPHERE" or gconf['nsph'] != 1):
             # TODO: put here a test to allow for empty vectors in case of random generation
             if (len(model['geometry_settings']['x_coords']) != gconf['nsph']):
@@ -315,13 +320,10 @@ def load_model(model_file):
             if (len(model['geometry_settings']['z_coords']) != gconf['nsph']):
                 print("ERROR: Z coordinates do not match the number of spheres!")
                 return (None, None)
-        gconf['vec_sph_x'] = [0.0 for i in range(gconf['nsph'])]
-        gconf['vec_sph_y'] = [0.0 for i in range(gconf['nsph'])]
-        gconf['vec_sph_z'] = [0.0 for i in range(gconf['nsph'])]
-        for si in range(gconf['nsph']):
-            gconf['vec_sph_x'][si] = float(model['geometry_settings']['x_coords'][si])
-            gconf['vec_sph_y'][si] = float(model['geometry_settings']['y_coords'][si])
-            gconf['vec_sph_z'][si] = float(model['geometry_settings']['z_coords'][si])
+            for si in range(gconf['nsph']):
+                gconf['vec_sph_x'][si] = float(model['geometry_settings']['x_coords'][si])
+                gconf['vec_sph_y'][si] = float(model['geometry_settings']['y_coords'][si])
+                gconf['vec_sph_z'][si] = float(model['geometry_settings']['z_coords'][si])
     else: # model is None
         print("ERROR: could not parse " + model_file + "!")
     return (sconf, gconf)
@@ -507,7 +509,7 @@ def write_legacy_sconf(conf):
     output = open(out_file, 'w')
     str_line = " {0:3d}{1:3d}\n".format(nsph, ies)
     output.write(str_line)
-    str_line = " {0:12.7E} {1:12.7E} {2:12.7E} {3:2d} {4:7d} {5:4d} {6:3d}\n".format(
+    str_line = " {0:12.7E} {1:12.7E} {2:12.7E} {3:2d} {4:4d} {5:4d} {6:3d}\n".format(
         exdc, wp, xip, idfc, nxi, instpc, xi_flag
     )
     output.write(str_line)
diff --git a/test_data/cluster/config_dev.yml b/test_data/cluster/config_dev.yml
new file mode 100644
index 0000000000000000000000000000000000000000..5a9fd27de76b94b81cc2d43522f82c8deaa39a7f
--- /dev/null
+++ b/test_data/cluster/config_dev.yml
@@ -0,0 +1,146 @@
+system_settings:
+  # Limit on host RAM use in Gb (0 for no configuration limit)
+  max_host_ram : 0
+  # Limit on GPU RAM use in Gb ( 0 for no configuration limit)
+  max_gpu_ram  : 0
+
+input_settings:
+  # Folder to write the code input configuration files
+  input_folder : "."
+  # Name of the scatterer description file
+  spheres_file : "DEDFB"
+  # Name of the geometry description file
+  geometry_file: "DCLU"
+  
+output_settings:
+  # Folder for the code output storage
+  output_folder: "test_subdir"
+  # Name of the main output file
+  output_name  : "c_OCLU"
+  # Requested output formats
+  formats      : [ "LEGACY", "HDF5" ]
+  # Index of the scale for transition matrix output
+  jwtm         : 1
+
+particle_settings:
+  # What application to use (SPHERE | CLUSTER | INCLUSION)
+  application : "CLUSTER"
+  # Number of spheres
+  n_spheres   : 4
+  # Number of sphere types
+  n_types     : 4
+  # Vector of sphere type identifiers (what type is each sphere)
+  sph_types   : [ 1, 2, 3, 4 ]
+  # Vector of layers in types (how many layers in each type)
+  n_layers    : [ 1, 1, 1, 1 ]
+  # Spherical monomer radii in m (one size for each type)
+  radii       : [ 5.0e-8, 4.0e-8, 7.5e-8, 3.0e-8 ]
+  # Layer fractional radii (one per layer in each type)
+  rad_frac    : [ [ 1.0 ], [ 1.0 ], [ 1.0 ], [ 1.0 ] ]
+  # Index of the dielectric constants (one per layer in each type)
+  #
+  # 1 is first file in `dielec_file`, 2 is second ...
+  dielec_id   : [ [ 1 ], [ 1 ], [ 1 ], [ 1 ] ]
+
+material_settings:
+  diel_flag   : -1
+  # External medium dielectric constant
+  extern_diel : 2.3104e0
+  # Dielectric constant files folder
+  dielec_path : "."
+  # List of dielectric constant files (used if diel_flag = 0)
+  dielec_file : [ ]
+  # Dielectric constant files format (same for all files)
+  dielec_fmt  : [ "CSV" ]
+  # Matching method between optical constants and radiation wavelengths
+  #
+  # INTERPOLATE: the constants are interpolated on wavelengths
+  # GRID: only the wavelengths with defined constants are computed
+  #
+  match_mode  : "GRID"
+  # Reference dielectric constants (used if diel_flag = -1)
+  #
+  # One real and one imaginary part for each layer in each type.
+  diel_const  : [
+    [ 3.25e00, 0.00e00 ],
+    [ 3.25e00, 0.00e00 ],
+    [ 3.25e00, 0.00e00 ],
+    [ 3.25e00, 0.00e00 ]
+    ]
+
+radiation_settings:
+  # Radiation field polarization (LINEAR | CIRCULAR)
+  polarization: "LINEAR"
+  # First scale to be used
+  scale_start : 1.0e00
+  # Last scale to be used
+  scale_end   : 2.0e00
+  # Calculation step (overridden if `match_mode` is GRID)
+  scale_step  : 1.0e00
+  # Peak Omega
+  wp          : 1.372e15
+  # Peak scale
+  xip         : 1.0e00
+  # Define scale explicitly (0) or in equal steps (1)
+  step_flag   : 0
+  # Type of scaling variable
+  scale_name  : "XI"
+
+geometry_settings:
+  # Maximum internal field expansion
+  li          : 8
+  # Maximum external field expansion (not used by SPHERE)
+  le          : 8
+  # Number of transition layer integration points
+  npnt        : 149
+  # Number of non transition layer integration points
+  npntts      : 300
+  # Averaging mode
+  iavm        : 0
+  # Meridional plane flag
+  isam        : 0
+  # Starting incidence azimuth angle
+  in_th_start : 79.0
+  # Incidence azimuth angle incremental step
+  in_th_step  : 10.0
+  # Ending incidence azimuth angle
+  in_th_end   : 89.0
+  # Starting incidence elevation angle
+  in_ph_start : 0.0
+  # Incidence elevation angle incremental step
+  in_ph_step  : 10.0
+  # Ending incidence elevation angle
+  in_ph_end   : 10.0
+  # Starting scattered azimuth angle
+  sc_th_start : 34.0
+  # Scattered azimuth angle incremental step
+  sc_th_step  : 15.0
+  # Ending scattered azimuth angle
+  sc_th_end   : 49.0
+  # Starting scattered elevation angle
+  sc_ph_start : 5.0
+  # Scattered elevation angle incremental step
+  sc_ph_step  : 5.0
+  # Ending scattered elevation angle
+  sc_ph_end   : 10.0
+  # Vector of sphere X coordinates (one per sphere or empty for random)
+  x_coords    : [
+    0.00e00,
+    0.00e00,
+    1.18882e-07,
+    -5.656855e-08,
+    ]
+  # Vector of sphere Y coordinates (one per sphere or empty for random)
+  y_coords    : [
+    3.00e-08,
+    3.00e-08,
+    3.00e-08,
+    -2.656855e-08
+    ]
+  # Vector of sphere Z coordinates (one per sphere or empty for random)
+  z_coords    : [
+    0.00e00,
+    9.00e-08,
+    3.8627e-08,
+    0.00e00
+    ]
diff --git a/test_data/config_example.yml b/test_data/config_example.yml
index 4a01e2d0b9295448f24e68de49bca9aaa19ed96a..7a9553c4ed41efd52f31984c5afcfdd81ceb09fc 100644
--- a/test_data/config_example.yml
+++ b/test_data/config_example.yml
@@ -10,13 +10,13 @@ input_settings:
   # Name of the scatterer description file
   spheres_file : "DEDFB"
   # Name of the geometry description file
-  geometry_file: "DSPH"
+  geometry_file: "DCLU"
   
 output_settings:
   # Folder for the code output storage
   output_folder: "test_subdir"
   # Name of the main output file
-  output_name  : "c_OSPH"
+  output_name  : "c_OCLU"
   # Requested output formats
   formats      : [ "LEGACY", "HDF5" ]
   # Index of the scale for transition matrix output
@@ -24,64 +24,73 @@ output_settings:
 
 particle_settings:
   # What application to use (SPHERE | CLUSTER | INCLUSION)
-  application : "SPHERE"
+  application : "CLUSTER"
   # Number of spheres
-  n_spheres   : 1
+  n_spheres   : 4
   # Number of sphere types
-  n_types     : 1
+  n_types     : 4
   # Vector of sphere type identifiers (what type is each sphere)
-  sph_types   : [ 1 ]
+  sph_types   : [ 1, 2, 3, 4 ]
   # Vector of layers in types (how many layers in each type)
-  n_layers    : [ 2 ]
+  n_layers    : [ 1, 1, 1, 1 ]
   # Spherical monomer radii in m (one size for each type)
-  radii       : [ 2.5e-7 ]
+  radii       : [ 5.0e-8, 4.0e-8, 7.5e-8, 3.0e-8 ]
   # Layer fractional radii (one per layer in each type)
-  rad_frac    : [ [ 0.5, 1.0 ] ]
+  rad_frac    : [ [ 1.0 ], [ 1.0 ], [ 1.0 ], [ 1.0 ] ]
   # Index of the dielectric constants (one per layer in each type)
   #
   # 1 is first file in `dielec_file`, 2 is second ...
-  dielec_id   : [ [ 1, 2 ] ]
+  dielec_id   : [ [ 1 ], [ 1 ], [ 1 ], [ 1 ] ]
 
 material_settings:
-  diel_flag   : 0
+  diel_flag   : -1
   # External medium dielectric constant
-  extern_diel : 1.0e0
+  extern_diel : 2.3104e0
   # Dielectric constant files folder
-  dielec_path : "../"
-  # List of dielectric constant files
-  dielec_file : [ "eps_draine_long_Si", "eps_ashok_long_C" ]
+  dielec_path : "."
+  # List of dielectric constant files (used if diel_flag = 0)
+  dielec_file : [ ]
   # Dielectric constant files format (same for all files)
   dielec_fmt  : [ "CSV" ]
   # Matching method between optical constants and radiation wavelengths
   #
-  # interpolate: the constants are interpolated on wavelengths
-  # grid: only the wavelengths with defined constants are computed
+  # INTERPOLATE: the constants are interpolated on wavelengths
+  # GRID: only the wavelengths with defined constants are computed
   #
   match_mode  : "GRID"
+  # Reference dielectric constants (used if diel_flag = -1)
+  #
+  # One real and one imaginary part for each layer in each type.
+  diel_const  : [
+    [ 3.25e00, 0.00e00 ],
+    [ 3.25e00, 0.00e00 ],
+    [ 3.25e00, 0.00e00 ],
+    [ 3.25e00, 0.00e00 ]
+    ]
 
 radiation_settings:
   # Radiation field polarization (LINEAR | CIRCULAR)
   polarization: "LINEAR"
   # First scale to be used
-  scale_start : 1.0e-7
+  scale_start : 1.0e00
   # Last scale to be used
-  scale_end   : 1.0e-6
+  scale_end   : 2.0e00
   # Calculation step (overridden if `match_mode` is GRID)
-  scale_step  : 5.0e-9
+  scale_step  : 1.0e00
   # Peak Omega
-  wp          : 3.0e8
+  wp          : 1.372e15
   # Peak scale
-  xip         : 1.0e0
+  xip         : 1.0e00
   # Define scale explicitly (0) or in equal steps (1)
   step_flag   : 0
-  # Type of scaling variable (only wavelength supported, for now)
-  scale_name  : "WAVELENGTH"
+  # Type of scaling variable
+  scale_name  : "XI"
 
 geometry_settings:
   # Maximum internal field expansion
-  li          : 20
+  li          : 8
   # Maximum external field expansion (not used by SPHERE)
-  le          : 20
+  le          : 8
   # Number of transition layer integration points
   npnt        : 149
   # Number of non transition layer integration points
@@ -91,33 +100,47 @@ geometry_settings:
   # Meridional plane flag
   isam        : 0
   # Starting incidence azimuth angle
-  in_th_start : 0.0
+  in_th_start : 79.0
   # Incidence azimuth angle incremental step
-  in_th_step  : 0.0
+  in_th_step  : 10.0
   # Ending incidence azimuth angle
-  in_th_end   : 0.0
+  in_th_end   : 89.0
   # Starting incidence elevation angle
   in_ph_start : 0.0
   # Incidence elevation angle incremental step
-  in_ph_step  : 0.0
+  in_ph_step  : 10.0
   # Ending incidence elevation angle
-  in_ph_end   : 0.0
+  in_ph_end   : 10.0
   # Starting scattered azimuth angle
-  sc_th_start : 0.0
+  sc_th_start : 34.0
   # Scattered azimuth angle incremental step
-  sc_th_step  : 0.0
+  sc_th_step  : 15.0
   # Ending scattered azimuth angle
-  sc_th_end   : 0.0
+  sc_th_end   : 49.0
   # Starting scattered elevation angle
-  sc_ph_start : 0.0
+  sc_ph_start : 5.0
   # Scattered elevation angle incremental step
-  sc_ph_step  : 0.0
+  sc_ph_step  : 5.0
   # Ending scattered elevation angle
-  sc_ph_end   : 0.0
+  sc_ph_end   : 10.0
   # Vector of sphere X coordinates (one per sphere or empty for random)
-  x_coords    : []
+  x_coords    : [
+    0.00e00,
+    0.00e00,
+    1.18882e-07,
+    -5.656855e-08,
+    ]
   # Vector of sphere Y coordinates (one per sphere or empty for random)
-  y_coords    : []
+  y_coords    : [
+    3.00e-08,
+    3.00e-08,
+    3.00e-08,
+    -2.656855e-08
+    ]
   # Vector of sphere Z coordinates (one per sphere or empty for random)
-  z_coords    : []
-  
\ No newline at end of file
+  z_coords    : [
+    0.00e00,
+    9.00e-08,
+    3.8627e-08,
+    0.00e00
+    ]
diff --git a/test_data/sphere/config_dev.yml b/test_data/sphere/config_dev.yml
new file mode 100644
index 0000000000000000000000000000000000000000..ac2d0e569d393da9238bee1562ef9f60964551ec
--- /dev/null
+++ b/test_data/sphere/config_dev.yml
@@ -0,0 +1,126 @@
+system_settings:
+  # Limit on host RAM use in Gb (0 for no configuration limit)
+  max_host_ram : 0
+  # Limit on GPU RAM use in Gb ( 0 for no configuration limit)
+  max_gpu_ram  : 0
+
+input_settings:
+  # Folder to write the code input configuration files
+  input_folder : "."
+  # Name of the scatterer description file
+  spheres_file : "DEDFB"
+  # Name of the geometry description file
+  geometry_file: "DSPH"
+  
+output_settings:
+  # Folder for the code output storage
+  output_folder: "."
+  # Name of the main output file
+  output_name  : "c_OSPH"
+  # Requested output formats
+  formats      : [ "LEGACY", "HDF5" ]
+  # Index of the scale for transition matrix output
+  jwtm         : 1
+
+particle_settings:
+  # What application to use (SPHERE | CLUSTER | INCLUSION)
+  application : "SPHERE"
+  # Number of spheres
+  n_spheres   : 1
+  # Number of sphere types
+  n_types     : 1
+  # Vector of sphere type identifiers (what type is each sphere)
+  sph_types   : [ 1 ]
+  # Vector of layers in types (how many layers in each type)
+  n_layers    : [ 2 ]
+  # Spherical monomer radii in m (one size for each type)
+  radii       : [ 4.0e-8 ]
+  # Layer fractional radii (one per layer in each type)
+  rad_frac    : [ [ 0.75, 1.0 ] ]
+  # Index of the dielectric constants (one per layer in each type)
+  #
+  # 1 is first file in `dielec_file`, 2 is second ...
+  dielec_id   : [ ]
+
+material_settings:
+  diel_flag   : -1
+  # External medium dielectric constant
+  extern_diel : 1.7689e+00
+  # Dielectric constant files folder
+  dielec_path : "."
+  # List of dielectric constant files (used if diel_flag = 0)
+  dielec_file : [ ]
+  # Dielectric constant files format (same for all files)
+  dielec_fmt  : [ "CSV" ]
+  # Matching method between optical constants and radiation wavelengths
+  #
+  # INTERPOLATE: the constants are interpolated on wavelengths
+  # GRID: only the wavelengths with defined constants are computed
+  #
+  match_mode  : "GRID"
+  # Reference dielectric constants (used if diel_flag = -1)
+  #
+  # One real and one imaginary part for each layer in each type.
+  diel_const  : [
+    [ 2.25e00, 2.56e-03, 2.25e00, 2.56e-03 ]
+    ]
+
+radiation_settings:
+  # Radiation field polarization (LINEAR | CIRCULAR)
+  polarization: "LINEAR"
+  # First scale to be used
+  scale_start : 1.0e+00
+  # Last scale to be used
+  scale_end   : 2.0e+00
+  # Calculation step (overridden if `match_mode` is GRID)
+  scale_step  : 1.0e+00
+  # Peak Omega
+  wp          : 1.7715748e+15
+  # Peak scale
+  xip         : 1.0e+00
+  # Define scale explicitly (0) or in equal steps (1)
+  step_flag   : 0
+  # Type of scaling variable
+  scale_name  : "XI"
+
+geometry_settings:
+  # Maximum internal field expansion
+  li          : 8
+  # Maximum external field expansion (not used by SPHERE)
+  le          : 8
+  # Number of transition layer integration points
+  npnt        : 149
+  # Number of non transition layer integration points
+  npntts      : 300
+  # Meridional plane flag
+  isam        : 0
+  # Starting incidence azimuth angle
+  in_th_start : 79.0
+  # Incidence azimuth angle incremental step
+  in_th_step  : 10.0
+  # Ending incidence azimuth angle
+  in_th_end   : 89.0
+  # Starting incidence elevation angle
+  in_ph_start : 0.0
+  # Incidence elevation angle incremental step
+  in_ph_step  : 10.0
+  # Ending incidence elevation angle
+  in_ph_end   : 10.0
+  # Starting scattered azimuth angle
+  sc_th_start : 34.0
+  # Scattered azimuth angle incremental step
+  sc_th_step  : 15.0
+  # Ending scattered azimuth angle
+  sc_th_end   : 49.0
+  # Starting scattered elevation angle
+  sc_ph_start : 5.0
+  # Scattered elevation angle incremental step
+  sc_ph_step  : 5.0
+  # Ending scattered elevation angle
+  sc_ph_end   : 10.0
+  # Vector of sphere X coordinates (one per sphere or empty for random)
+  x_coords    : [ ]
+  # Vector of sphere Y coordinates (one per sphere or empty for random)
+  y_coords    : [ ]
+  # Vector of sphere Z coordinates (one per sphere or empty for random)
+  z_coords    : [ ]