diff --git a/Grid/serialisation/BaseIO.h b/Grid/serialisation/BaseIO.h
index dd15e7da..e96a2f9a 100644
--- a/Grid/serialisation/BaseIO.h
+++ b/Grid/serialisation/BaseIO.h
@@ -33,6 +33,7 @@ Author: Guido Cossu <guido.cossu@ed.ac.uk>
 #include <type_traits>
 #include <Grid/tensors/Tensors.h>
 #include <Grid/serialisation/VectorUtils.h>
+#include <Grid/Eigen/unsupported/CXX11/Tensor>
 
 namespace Grid {
   // Abstract writer/reader classes ////////////////////////////////////////////
@@ -60,6 +61,9 @@ namespace Grid {
     void write(const std::string &s, const iVector<U, N> &output);
     template <typename U, int N>
     void write(const std::string &s, const iMatrix<U, N> &output);
+    template <typename Scalar_, int NumIndices_>
+    void write(const std::string &s, const Eigen::Tensor<Scalar_, NumIndices_, Eigen::RowMajor> &output);
+    
     void         scientificFormat(const bool set);
     bool         isScientific(void);
     void         setPrecision(const unsigned int prec);
@@ -162,6 +166,14 @@ namespace Grid {
     upcast->writeDefault(s, tensorToVec(output));
   }
 
+  template <typename T>
+  template <typename Scalar_, int NumIndices_>
+  void Writer<T>::write(const std::string &s, const Eigen::Tensor<Scalar_, NumIndices_, Eigen::RowMajor> &output)
+  {
+    //upcast->writeDefault(s, tensorToVec(output));
+    std::cout << "I really should add code to write Eigen::Tensor ..." << std::endl;
+  }
+  
   template <typename T>
   void Writer<T>::scientificFormat(const bool set)
   {
diff --git a/Hadrons/Modules/MDistil/BContraction.hpp b/Hadrons/Modules/MDistil/BContraction.hpp
index c0843273..7fc5fc7c 100644
--- a/Hadrons/Modules/MDistil/BContraction.hpp
+++ b/Hadrons/Modules/MDistil/BContraction.hpp
@@ -151,8 +151,8 @@ void TBContraction<FImpl>::execute(void)
        Gamma::Algebra::SigmaYT, // C = i gamma_2 gamma_4
        Gamma::Algebra::GammaYGamma5, // i gamma_4 C gamma_5 = i gamma_2 gamma_5
     };
-    std::vector<Complex> factor23 = {(0.,-1.),(0.,1.),(0.,1.)};
-
+  std::vector<Complex> factor23{{0.,-1.},{0.,1.},{0.,1.}};
+    if((0)){
     for (int i1=0 ; i1 < N_1 ; i1++){
       for (int i2=0 ; i2 < N_2 ; i2++){
         for (int i3=0 ; i3 < N_3 ; i3++){
@@ -175,7 +175,7 @@ void TBContraction<FImpl>::execute(void)
 		  // this should be outerProduct??? Does not work.
                   for (int isl=0 ; isl < 4 ; isl++){
                     for (int isr=0 ; isr < 4 ; isr++){
-		      diquark()(isl,isr)() = factor23[0]*tmp22s()(isl)(),tmp333()(isr)();
+		      diquark()(isl,isr)() = factor23[0]*tmp22s()(isl)()*tmp333()(isr)();
 		    }
 		  }
 		  // Is there a way to compute gamma * SpinMatrix (left component)???
@@ -203,6 +203,43 @@ void TBContraction<FImpl>::execute(void)
       Bindex = 0 + N_1*(0 + N_2*(0 + N_3*(0+Nmom*t)));
       std::cout << "BaryonField(t=" << t << ") = " << BField[Bindex] << std::endl;
     }
+    } // end if 0
+    
+    // ONLY THIS IS CORRECT?
+    std::vector<SpinVector> BField2(Nmom*Nt*N_1*N_2*N_3);    
+    Complex diquark2;
+    for (int i1=0 ; i1 < N_1 ; i1++){
+      for (int i2=0 ; i2 < N_2 ; i2++){
+        for (int i3=0 ; i3 < N_3 ; i3++){
+          for (int imom=0 ; imom < Nmom ; imom++){
+            for (int t=0 ; t < Nt ; t++){
+              Bindex = i1 + N_1*(i2 + N_2*(i3 + N_3*(imom+Nmom*t)));
+	      ExtractSliceLocal(tmp1,one[i1],0,t,3);
+              parallel_for (unsigned int sU = 0; sU < grid3d->oSites(); ++sU)
+              {
+                for (int ie=0 ; ie < 6 ; ie++){
+		  // Why does peekColour not work????
+                  for (int is=0 ; is < 4 ; is++){
+	            tmp11s()(is)() = tmp11[sU]()(is)(epsilon[ie][0]);
+	            tmp22s()(is)() = tmp22[sU]()(is)(epsilon[ie][1]);
+	            tmp33s()(is)() = tmp33[sU]()(is)(epsilon[ie][2]);
+		  }
+		  tmp333 = Gamma(gamma23_[0])*tmp33s;
+                  diquark2 = factor23[0]*innerProduct(tmp22s,tmp333);
+                  BField2[Bindex]+=(double)epsilon_sgn[ie]*tmp11s*diquark2;
+		}
+  	      }
+            }
+	  }
+	}
+      }
+    }
+    for (int is=0 ; is < 4 ; is++){
+      for (int t=0 ; t < Nt ; t++){
+        Bindex = 0 + N_1*(0 + N_2*(0 + N_3*(0+Nmom*t)));
+        std::cout << "BaryonField(is=" << is << ",t=" << t << ") = " << BField2[Bindex]()(is)() << std::endl;
+      }
+    }
 }
 
 END_MODULE_NAMESPACE
diff --git a/tests/hadrons/Test_hadrons_distil.cc b/tests/hadrons/Test_hadrons_distil.cc
index 33f883ee..1529eb15 100644
--- a/tests/hadrons/Test_hadrons_distil.cc
+++ b/tests/hadrons/Test_hadrons_distil.cc
@@ -422,6 +422,51 @@ bool DebugEigenTest()
 
   return true;
 }
+
+typedef iMatrix<Complex,7> OddBall;
+
+// From Test_serialisation.cc
+class myclass: Serializable {
+public:
+  GRID_SERIALIZABLE_CLASS_MEMBERS(myclass
+                                  , OddBall, critter
+                                  , SpinColourVector, scv
+                                  , SpinColourMatrix, scm
+                                  );
+};
+
+template <typename W, typename R, typename O>
+bool ioTest(const std::string &filename, const O &object, const std::string &name)
+{
+  // writer needs to be destroyed so that writing physically happens
+  {
+    W writer(filename);
+    write(writer, "testobject", object);
+  }
+  
+  /*R    reader(filename);
+  O    buf;
+  bool good;
+  
+  read(reader, "testobject", buf);
+  good = (object == buf);
+  std::cout << name << " IO test: " << (good ? "success" : "failure");
+  std::cout << std::endl;
+  return good;*/
+  return true;
+}
+
+bool DebugIOTest(void) {
+  OddBall critter;
+  ioTest<Hdf5Writer, Hdf5Reader, OddBall>("iotest_oddball.h5", critter, "OddBall");
+  SpinColourMatrix scm;
+  ioTest<Hdf5Writer, Hdf5Reader, SpinColourMatrix>("iotest_matrix.h5", scm, "SpinColourMatrix");
+  SpinColourVector scv;
+  ioTest<Hdf5Writer, Hdf5Reader, SpinColourVector>("iotest_vector.h5", scv, "SpinColourVector");
+  myclass o;
+  ioTest<Hdf5Writer, Hdf5Reader, myclass>("iotest_object.h5", o, "myclass_object_instance_name");
+  return true;
+}
 #endif
 
 int main(int argc, char *argv[])
@@ -430,7 +475,8 @@ int main(int argc, char *argv[])
   // Debug only - test of Eigen::Tensor
   std::cout << "sizeof(std::streamsize) = " << sizeof(std::streamsize) << std::endl;
   std::cout << "sizeof(Eigen::Index) = " << sizeof(Eigen::Index) << std::endl;
-  if( DebugEigenTest() ) return 0;
+  //if( DebugEigenTest() ) return 0;
+  if(DebugIOTest()) return 0;
 #endif
 
   // Decode command-line parameters. 1st one is which test to run