diff --git a/src/loadGadget.cpp b/src/loadGadget.cpp
index 4cb3d7d..763d2c6 100644
--- a/src/loadGadget.cpp
+++ b/src/loadGadget.cpp
@@ -65,95 +65,152 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags)
 {
   SimuData *data;
   int p, n;
-  UnformattedRead f(fname);
+  UnformattedRead *f;
   GadgetHeader h;
 
+  if (id >= 0) {
+    int numDigits = 1;
+    int mul = 1;
+    while (mul < id)
+      {
+	mul *= 10;
+	numDigits++;
+      }
+    size_t len = numDigits+2+strlen(fname);
+    char *out_fname = new char[numDigits+2+strlen(fname)];
+
+    if (snprintf(out_fname, len, "%s.%d", fname, id) != len)
+      abort();
+     
+    f = new UnformattedRead(out_fname);
+    if (f == 0)
+      return 0;
+
+    delete out_fname;
+
+  } else {
+
+    f = new UnformattedRead(fname);
+    if (f == 0)
+      return 0;
+
+  }
+
   data = new SimuData;
+  if (data == 0) {
+    delete f;
+    return 0;
+  }
 
-  f.beginCheckpoint();
+  f->beginCheckpoint();
   for (int i = 0; i < 6; i++)
-    h.npart[i] = f.readInt32();
+    h.npart[i] = f->readInt32();
   for (int i = 0; i < 6; i++)    
-    h.mass[i] = f.readReal64();
-  h.time = f.readReal64();
-  h.redshift = f.readReal64();
-  h.flag_sfr = f.readInt32();
-  h.flag_feedback = f.readInt32();
+    h.mass[i] = f->readReal64();
+  h.time = f->readReal64();
+  h.redshift = f->readReal64();
+  h.flag_sfr = f->readInt32();
+  h.flag_feedback = f->readInt32();
   for (int i = 0; i < 6; i++)
-    h.npartTotal[i] = f.readInt32();
-  h.flag_cooling = f.readInt32();
-  h.num_files = f.readInt32();
-  h.BoxSize = f.readReal64();
-  h.Omega0 = f.readReal64();
-  h.OmegaLambda = f.readReal64();
-  h.HubbleParam = f.readReal64();
-  f.endCheckpoint(true);
+    h.npartTotal[i] = f->readInt32();
+  h.flag_cooling = f->readInt32();
+  h.num_files = f->readInt32();
+  data->BoxSize = h.BoxSize = f->readReal64();
+  h.Omega0 = f->readReal64();
+  h.OmegaLambda = f->readReal64();
+  h.HubbleParam = f->readReal64();
+  f->endCheckpoint(true);
 
-  long NumPart = 0;
+  long NumPart = 0, NumPartTotal = 0;
   for(int k=0; k<5; k++)
-    NumPart += h.npart[k];
+    {
+      NumPart += h.npart[k];
+      NumPartTotal += h.npartTotal[k];
+    }
   data->NumPart = NumPart;
+  data->TotalNumPart = NumPartTotal;
 
   if (loadflags & NEED_POSITION) {
       
-    for (int i = 0; i < 3; i++)
+    for (int i = 0; i < 3; i++) {
       data->Pos[i] = new float[data->NumPart];
-    
-    f.beginCheckpoint();
+      if (data->Pos[i] == 0) {
+	delete data;
+	return 0;
+      }
+    }
+
+    f->beginCheckpoint();
     for(int k = 0, p = 0; k < 5; k++) {
       for(int n = 0; n < h.npart[k]; n++) {
-	data->Pos[p][0] = f.readReal32();
-	data->Pos[p][1] = f.readReal32();
-	data->Pos[p][2] = f.readReal32();
+	data->Pos[p][0] = f->readReal32();
+	data->Pos[p][1] = f->readReal32();
+	data->Pos[p][2] = f->readReal32();
 	p++;
       }
     }
-    f.endCheckpoint();
+    f->endCheckpoint();
     
   } else {
     // Skip positions
-    f.skip(NumPart * 3 + 2*4);
+    f->skip(NumPart * 3 * sizeof(float) + 2*4);
   }
   
   if (loadflags & NEED_VELOCITY) {
     for (int i = 0; i < 3; i++)
-      data->Vel[i] = new float[data->NumPart];
-    
-    f.beginCheckpoint();
+      {
+	data->Vel[i] = new float[data->NumPart];
+	if (data->Vel[i] == 0) 
+	  {
+	    delete data;
+	    return 0;
+	  }
+      }
+
+    f->beginCheckpoint();
     for(int k = 0, p = 0; k < 5; k++) {
       for(int n = 0; n < h.npart[k]; n++) {
-	data->Vel[p][0] = f.readReal32();
-	data->Vel[p][1] = f.readReal32();
-	data->Vel[p][2] = f.readReal32();
+	data->Vel[p][0] = f->readReal32();
+	data->Vel[p][1] = f->readReal32();
+	data->Vel[p][2] = f->readReal32();
 	p++;
       }
     }
-    f.endCheckpoint();
+    f->endCheckpoint();
     
     // TODO: FIX THE UNITS OF THESE FUNKY VELOCITIES !!!
   } else {
     // Skip velocities
-    f.skip(NumPart*3+2*4);
+    f->skip(NumPart*3*sizeof(float)+2*4);
   }
 
   // Skip ids
   if (loadflags & NEED_GADGET_ID) {
-    f.beginCheckpoint();
+    f->beginCheckpoint();
+    data->Id = new int[data->NumPart];
+    if (data->Id == 0)
+      {
+	delete data;
+	return 0;
+      }
+
     for(int k = 0, p = 0; k < 6; k++)
       {
 	for(int n = 0; n < h.npart[k]; n++)
 	  {
-	    data->Id[p] = f.readInt32();
+	    data->Id[p] = f->readInt32();
 	    p++;
 	  }
       }
-    f.endCheckpoint();
+    f->endCheckpoint();
   } else {
-    f.skip(2*4);
+    f->skip(2*4);
     for (int k = 0; k < 6; k++)
-      f.skip(h.npart[k]);
+      f->skip(h.npart[k]*4);
   }
 
+  delete f;
+
   return data;
 }
 
diff --git a/src/loadSimu.hpp b/src/loadSimu.hpp
index 1215a88..965f266 100644
--- a/src/loadSimu.hpp
+++ b/src/loadSimu.hpp
@@ -16,6 +16,7 @@ namespace CosmoTool
     float time;
 
     long NumPart;
+    long TotalNumPart;
     int *Id;
     float *Pos[3];
     float *Vel[3];