diff --git a/Grid/tensors/Tensor_Ta.h b/Grid/tensors/Tensor_Ta.h
index 15120ffd..cc21378e 100644
--- a/Grid/tensors/Tensor_Ta.h
+++ b/Grid/tensors/Tensor_Ta.h
@@ -212,58 +212,52 @@ accelerator_inline iMatrix<vtype,N> ProjectOnSpGroup(const iMatrix<vtype,N> &arg
   iMatrix<vtype,N> ret(arg);
   vtype nrm;
   vtype inner;
-  vtype tmp;
   
   for(int c1=0;c1<N/2;c1++)
   {
       
-      for (int b=0; b<c1; b++)                  // remove the b-rows from U_c1
-      {
-          decltype(ret._internal[b][b]*ret._internal[b][b]) pr;
-          decltype(ret._internal[b][b]*ret._internal[b][b]) prn;
-          zeroit(pr);
-          zeroit(prn);
+    for (int b=0; b<c1; b++)                  // remove the b-rows from U_c1
+    {
+      decltype(ret._internal[b][b]*ret._internal[b][b]) pr;
+      decltype(ret._internal[b][b]*ret._internal[b][b]) prn;
+      zeroit(pr);
+      zeroit(prn);
           
-          for(int c=0; c<N; c++)
-          {
-              pr += conjugate(ret._internal[c1][c])*ret._internal[b][c];        // <U_c1 | U_b >
-              prn += conjugate(ret._internal[c1][c])*ret._internal[b+N/2][c];   // <U_c1 | U_{b+N} >
-          }
+      for(int c=0; c<N; c++)
+      {
+        pr += conjugate(ret._internal[c1][c])*ret._internal[b][c];        // <U_c1 | U_b >
+        prn += conjugate(ret._internal[c1][c])*ret._internal[b+N/2][c];   // <U_c1 | U_{b+N} >
+      }
        
 
-          for(int c=0; c<N; c++)
-          {
-              ret._internal[c1][c] -= (conjugate(pr) * ret._internal[b][c] + conjugate(prn) * ret._internal[b+N/2][c] );    //  U_c1 -= (  <U_c1 | U_b > U_b + <U_c1 | U_{b+N} > U_{b+N}  )
-          }
+      for(int c=0; c<N; c++)
+      {
+        ret._internal[c1][c] -= (conjugate(pr) * ret._internal[b][c] + conjugate(prn) * ret._internal[b+N/2][c] );    //  U_c1 -= (  <U_c1 | U_b > U_b + <U_c1 | U_{b+N} > U_{b+N}  )
       }
+    }
     
-      zeroit(inner);
-      for(int c2=0;c2<N;c2++)
-      {
-          inner += innerProduct(ret._internal[c1][c2],ret._internal[c1][c2]);
-      }
+    zeroit(inner);
+    for(int c2=0;c2<N;c2++)
+    {
+      inner += innerProduct(ret._internal[c1][c2],ret._internal[c1][c2]);
+    }
       
-      nrm = sqrt(inner);
-      nrm = 1.0/nrm;
-      for(int c2=0;c2<N;c2++)
-      {
-          ret._internal[c1][c2]*= nrm;
-      }
+    nrm = sqrt(inner);
+    nrm = 1.0/nrm;
+    for(int c2=0;c2<N;c2++)
+    {
+      ret._internal[c1][c2]*= nrm;
+    }
       
-
-      for(int c2=0;c2<N/2;c2++)
-      {
-          tmp = conjugate(ret._internal[c1][c2]);       // (up-left)* of the old matrix
-          ret._internal[c1+N/2][c2+N/2] = tmp;          // down right in the new matrix = (up-left)* of the old matrix
-      }
-         
-      for(int c2=N/2;c2<N;c2++)
-      {
-          tmp = conjugate(ret._internal[c1][c2]);   // (up-right)* of the old
-          ret._internal[c1+N/2][c2-N/2] = -tmp;     // down left in the new matrix = -(up-right)* of the old
-      }
-    
+    for(int c2=0;c2<N/2;c2++)
+    {
+      ret._internal[c1+N/2][c2+N/2] = conjugate(ret._internal[c1][c2]);          // down right in the new matrix = (up-left)* of the old matrix
+    }
       
+    for(int c2=N/2;c2<N;c2++)
+    {
+      ret._internal[c1+N/2][c2-N/2] = -conjugate(ret._internal[c1][c2]);;     // down left in the new matrix = -(up-right)* of the old
+    }
   }
   return ret;
 }
diff --git a/tests/sp2n/Test_project_on_Sp.cc b/tests/sp2n/Test_project_on_Sp.cc
index e514106f..2f8dd32b 100644
--- a/tests/sp2n/Test_project_on_Sp.cc
+++ b/tests/sp2n/Test_project_on_Sp.cc
@@ -2,11 +2,121 @@
 
 using namespace Grid;
 
+int SpGroupQuiz (const LatticeColourMatrixD Uin)
+{
+    double vol = Uin.Grid()->gSites();
+    const int nsp = Nc / 2;
+    LatticeColourMatrixD Omega(Uin.Grid());
+    Sp<Nc>::Omega(Omega);
+    LatticeColourMatrixD aux(Uin.Grid());
+    LatticeColourMatrixD identity(Uin.Grid());
+    Complex i(0., 1.);
+    identity = 1;
+    
+    std::cout << GridLogMessage << "Check matrix is non-zero " << std::endl;
+    assert(norm2(Uin) > 1e-8);
+    aux = Uin*adj(Uin) - identity;
+    std::cout << GridLogMessage << "U adjU - 1 = " << norm2(aux)/vol << std::endl;
+    assert(norm2(aux) < 1e-8);
+    aux = Omega - (Uin * Omega * transpose(Uin));
+    std::cout << GridLogMessage << "Omega - U Omega transpose(U) = " << norm2(aux)/vol << std::endl;
+    assert(norm2(aux) < 1e-8);
+    //Sp<Nc>::OmegaInvariance(Uin)
+    std::cout << GridLogMessage << "Checking the structure is " << std::endl;
+    std::cout << GridLogMessage << "U  =  (   W    X   )  " << std::endl;
+    std::cout << GridLogMessage << "      (  -X^*  W^* )  " << std::endl;
+    for (int c1 = 0; c1 < nsp; c1++) //check on W
+    {
+        for (int c2 = 0; c2 < nsp; c2++)
+        {
+            auto W = PeekIndex<ColourIndex>(Uin,c1,c2);
+            auto Wstar =  PeekIndex<ColourIndex>(Uin,c1+nsp,c2+nsp);
+            auto Ww = conjugate( Wstar );
+            auto amizero = sum(W - Ww);
+            auto amizeroo = TensorRemove(amizero);
+            assert(  amizeroo.real() < 10e-6 );
+            amizeroo *= i;
+            assert(  amizeroo.real() < 10e-6 );
+        }
+    }
+    for (int c1 = 0; c1 < nsp ; c1++)
+    {
+        for (int c2 = 0; c2 < nsp; c2++)
+        {
+            auto X = PeekIndex<ColourIndex>(Uin,c1,c2+nsp);
+            auto minusXstar = PeekIndex<ColourIndex>(Uin,c1+nsp,c2);
+            auto minusXx = conjugate(minusXstar);
+            auto amizero = sum (X + minusXx);
+            auto amizeroo = TensorRemove(amizero);
+            assert(  amizeroo.real() < 10e-6 );
+            amizeroo *= i;
+            assert(  amizeroo.real() < 10e-6 );
+        }
+    }
+    std::cout << GridLogMessage << "|Det| = " << norm2( Determinant(Uin) ) / vol << std::endl;
+    assert( norm2( Determinant(Uin) ) / vol - 1 < 1e-8);
+    return 0;
+}
+
+int SpAntiHermitianAlgebraQuiz (const LatticeColourMatrixD Uin)
+{
+    double vol = Uin.Grid()->gSites();
+    const int nsp = Nc / 2;
+    LatticeColourMatrixD Omega(Uin.Grid());
+    Sp<Nc>::Omega(Omega);
+    LatticeColourMatrixD aux(Uin.Grid());
+    LatticeColourMatrixD identity(Uin.Grid());
+    Complex i(0., 1.);
+    identity = 1;
+    
+    std::cout << GridLogMessage << "Check matrix is non-zero " << std::endl;
+    assert(norm2(Uin) > 1e-8);
+    aux = Uin - adj(Uin);
+    std::cout << GridLogMessage << "SpTa ::: T - Tda = " << norm2(aux)/vol << std::endl;
+    aux = Uin + adj(Uin);
+    std::cout << GridLogMessage << "SpTa ::: T + Tda = " << norm2(aux)/vol << std::endl;
+    assert( norm2(aux) - 1 < 1e-8);
+    std::cout << GridLogMessage << "Check that Omega T Omega + conj(T) = 0 " << std::endl;
+    aux = Omega*Uin*Omega + conjugate(Uin);
+    assert( norm2(aux) < 1e-8);
+    std::cout << GridLogMessage << "Checking the structure is " << std::endl;
+    std::cout << GridLogMessage << "U  =  (   W    X   )  " << std::endl;
+    std::cout << GridLogMessage << "      (  -X^*  W^* )  " << std::endl;
+    for (int c1 = 0; c1 < nsp; c1++) //check on W
+    {
+        for (int c2 = 0; c2 < nsp; c2++)
+        {
+            auto W = PeekIndex<ColourIndex>(Uin,c1,c2);
+            auto Wstar =  PeekIndex<ColourIndex>(Uin,c1+nsp,c2+nsp);
+            auto Ww = conjugate( Wstar );
+            auto amizero = sum(W - Ww);
+            auto amizeroo = TensorRemove(amizero);
+            assert(  amizeroo.real() < 10e-6 );
+            amizeroo *= i;
+            assert(  amizeroo.real() < 10e-6 );
+        }
+    }
+    for (int c1 = 0; c1 < nsp ; c1++)
+    {
+        for (int c2 = 0; c2 < nsp; c2++)
+        {
+            auto X = PeekIndex<ColourIndex>(Uin,c1,c2+nsp);
+            auto minusXstar = PeekIndex<ColourIndex>(Uin,c1+nsp,c2);
+            auto minusXx = conjugate(minusXstar);
+            auto amizero = sum (X + minusXx);
+            auto amizeroo = TensorRemove(amizero);
+            assert(  amizeroo.real() < 10e-6 );
+            amizeroo *= i;
+            assert(  amizeroo.real() < 10e-6 );
+        }
+    }
+    return 0;
+}
+
 int main (int argc, char **argv)
 {
     Grid_init(&argc,&argv);
     
-    
     Coordinate latt_size   = GridDefaultLatt();
     Coordinate simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
     Coordinate mpi_layout  = GridDefaultMpi();
@@ -14,27 +124,18 @@ int main (int argc, char **argv)
     GridCartesian             Grid(latt_size,simd_layout,mpi_layout);
     GridRedBlackCartesian     RBGrid(&Grid);
     
-    
     LatticeGaugeField Umu(&Grid);
     LatticeColourMatrixD U(&Grid);
-    LatticeColourMatrixD Up(&Grid);
     LatticeColourMatrixD aux(&Grid);
     LatticeColourMatrixD identity(&Grid);
-    
-    // Will test resimplectification-related functionalities (from ProjectOnGaugeGroup, ProjectOnSpGroup, ProjectGn) and projection on the algebra (from ProjectSp2nAlgebra)
-    // we work with matrices with positive determinant so detU = 1 even if in principle ProjectOnGaugeGroup and ProjectOnSpGroup allow for detU=-1
-    // so the checks will be the same for the three functions
-    // NB only ProjectGn is the proper simplectification function
-    
-    
-    const int nsp = Nc / 2;
-    
+    LatticeColourMatrixD Omega(&Grid);
+    Sp<Nc>::Omega(Omega);
+
     identity = 1.0;
-    RealD epsilon = 0.01;
-    RealD Delta = 666.;
+    RealD epsilon = 0.1;
     Complex i(0., 1.);
-    RealD u = 0.;
     double vol = Umu.Grid()->gSites();
+    const int nsp = Nc / 2;
     
     std::vector<int> pseeds({1,2,3,4,5});
     GridParallelRNG  pRNG(&Grid); pRNG.SeedFixedIntegers(pseeds);
@@ -42,236 +143,76 @@ int main (int argc, char **argv)
     SU<Nc>::HotConfiguration(pRNG,Umu);
     U = PeekIndex<LorentzIndex>(Umu,0);
     
-    aux = U*adj(U) - identity;
     std::cout << GridLogMessage << "Starting with random SUn matrix " << std::endl;
-    std::cout << GridLogMessage << "Unitary check " << std::endl;
+    aux = U*adj(U) - identity;
     std::cout <<GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl;
     assert ( norm2(aux) < 1e-8 );
-    std::cout <<GridLogMessage << std::endl;
-    if (Nc != 2)
+    if (Nc != 2)    // Sp2 = SU2
     {
-        std::cout << GridLogMessage << "This matrix should not leave Omega invariant, expect a warning" << std::endl;
+        std::cout << GridLogMessage <<  "Checking matrix does NOT leave Omega invariant, to avoid a trivial test " << std::endl;
+        aux = Omega - (U * Omega * transpose(U));
+        assert ( norm2(aux) > 1e-8 );
     }
-    Sp<Nc>::OmegaInvariance(U);
-    std::cout <<GridLogMessage << std::endl;
-    
-    U = U + epsilon*identity;
-    aux = U*adj(U) - identity;
     
+    U = U + epsilon*identity - i*identity;
     std::cout << GridLogMessage << "Unitary matrix deformed " << std::endl;
+    auto det =  sum( Determinant(U) );
+    std::cout << GridLogMessage << "Re(Det) = " << real(det) / vol << std::endl;
+    det = det*i;
+    std::cout << GridLogMessage << "Im(Det) = " << real(det) / vol << std::endl;
+    aux = U*adj(U) - identity;
     std::cout << GridLogMessage << "now U adjU - 1 = " << norm2(aux) << std::endl;
     
-    std::cout <<GridLogMessage << std::endl;
-    std::cout << GridLogMessage << "#   #   #   #" << std::endl;
-    std::cout << GridLogMessage << "Group" << std::endl;
-    std::cout << GridLogMessage << "#   #   #   #" << std::endl;
-    std::cout <<GridLogMessage << std::endl;
-    //  Testing ProjectOnSpGroup
+    //  ProjectOnSpGroup
     std::cout << GridLogMessage << "Testing ProjectOnSpGroup" << std::endl;
     std::cout << GridLogMessage << "Apply ProjectOnSpGroup to deformed matrix" << std::endl;
     U = ProjectOnSpGroup(U);
-
-    aux = U*adj(U) - identity;
-    std::cout << GridLogMessage << "Unitary check after ProjectOnSpGroup " << std::endl;
-    std::cout << GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl;
-    assert( norm2(aux) < 1e-8);
-    
-    // actual sp2n check
-    std::cout << GridLogMessage << "Checking Omega invariance after ProjectOnSpGroup" << std::endl;
-    Sp<Nc>::OmegaInvariance(U);     // no assertion here, but the next check will kill us if we are not simplectic
-    
-    // checks on elements
-    
-    std::cout << GridLogMessage << "Checking the structure is " << std::endl;
-    std::cout << GridLogMessage << "U  =  (   W    X   )  " << std::endl;
-    std::cout << GridLogMessage << "      (  -X^*  W^* )  " << std::endl;
-    std::cout <<GridLogMessage << std::endl;
-    for (int c1 = 0; c1 < nsp; c1++) //check on W
-    {
-        for (int c2 = 0; c2 < nsp; c2++)
-        {
-            auto W = PeekIndex<ColourIndex>(U,c1,c2);
-            auto Wstar =  PeekIndex<ColourIndex>(U,c1+nsp,c2+nsp);
-            auto Ww = conjugate( Wstar );
-            auto amizero = sum(W - Ww);
-            auto amizeroo = TensorRemove(amizero);
-            assert(  amizeroo.real() < 10e-6 );
-            amizeroo *= i;
-            assert(  amizeroo.real() < 10e-6 );
-        }
-        
-    }
-
-    for (int c1 = 0; c1 < nsp ; c1++)
-    {
-        for (int c2 = 0; c2 < nsp; c2++)
-        {
-            auto X = PeekIndex<ColourIndex>(U,c1,c2+nsp);
-            auto minusXstar = PeekIndex<ColourIndex>(U,c1+nsp,c2);
-            auto minusXx = conjugate(minusXstar);
-            auto amizero = sum (X + minusXx);
-            auto amizeroo = TensorRemove(amizero);
-            assert(  amizeroo.real() < 10e-6 );
-            amizeroo *= i;
-            assert(  amizeroo.real() < 10e-6 );
-        }
-    }
+    std::cout << GridLogMessage << "Run checks:" << std::endl;
+    SpGroupQuiz(U);
+    det =  sum( Determinant(U) );
+    std::cout << GridLogMessage << "Re(Det) after ProjectOnSpGroup (nothing to assert) = " << real(det) / vol << std::endl;
+    det = det*i;
+    std::cout << GridLogMessage << "Im(Det) after ProjectOnSpGroup (nothing to assert) = " << real(det) / vol << std::endl;
     
+    //  ProjectOnGaugeGroup
     SU<Nc>::HotConfiguration(pRNG,Umu); //refresh
-    U = PeekIndex<LorentzIndex>(Umu,1);
-    
+    U = PeekIndex<LorentzIndex>(Umu,0);
     std::cout << GridLogMessage << "Testing ProjectOnGaugeGroup" << std::endl;
-    U = U + Delta*identity;
+    U = U + epsilon*identity - i*identity;
     std::cout << GridLogMessage << "Apply ProjectOnGaugeGroup to deformed matrix" << std::endl;
     Sp<Nc>::ProjectOnGaugeGroup(U);
-    aux = U*adj(U) - identity;
-    std::cout << GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl;
-    assert( norm2(aux) < 1e-8);
-    Sp<Nc>::OmegaInvariance(U);
-    
-    std::cout << GridLogMessage << "Checking the structure is " << std::endl;
-    std::cout << GridLogMessage << "U  =  (   W    X   )  " << std::endl;
-    std::cout << GridLogMessage << "      (  -X^*  W^* )  " << std::endl;
-    
-    for (int c1 = 0; c1 < nsp; c1++) //check on W
-    {
-        for (int c2 = 0; c2 < nsp; c2++)
-        {
-            auto W = PeekIndex<ColourIndex>(U,c1,c2);
-            auto Wstar =  PeekIndex<ColourIndex>(U,c1+nsp,c2+nsp);
-            auto Ww = conjugate( Wstar );
-            auto amizero = sum(W - Ww);
-            auto amizeroo = TensorRemove(amizero);
-            assert(  amizeroo.real() < 10e-6 );
-            amizeroo *= i;
-            assert(  amizeroo.real() < 10e-6 );
-        }
-        
-    }
-
-    for (int c1 = 0; c1 < nsp ; c1++)
-    {
-        for (int c2 = 0; c2 < nsp; c2++)
-        {
-            auto X = PeekIndex<ColourIndex>(U,c1,c2+nsp);
-            auto minusXstar = PeekIndex<ColourIndex>(U,c1+nsp,c2);
-            auto minusXx = conjugate(minusXstar);
-            auto amizero = sum (X + minusXx);
-            auto amizeroo = TensorRemove(amizero);
-            assert(  amizeroo.real() < 10e-6 );
-            amizeroo *= i;
-            assert(  amizeroo.real() < 10e-6 );
-        }
-    }
+    std::cout << GridLogMessage <<"Run checks:" << std::endl;
+    SpGroupQuiz(U);
+    det =  sum( Determinant(U) );
+    std::cout << GridLogMessage << "Re(Det) after ProjectOnGaugeGroup (nothing to assert) = " << real(det) / vol << std::endl;
+    det = det*i;
+    std::cout << GridLogMessage << "Im(Det) after ProjectOnGaugeGroup (nothing to assert) = " << real(det) / vol << std::endl;
     
+    //  ProjectGn
+    SU<Nc>::HotConfiguration(pRNG,Umu); //refresh
+    U = PeekIndex<LorentzIndex>(Umu,0);
     std::cout << GridLogMessage << "Testing ProjectGn" << std::endl;
-    U = U + Delta*identity;
+    U = U + epsilon*identity - i*identity;
     std::cout << GridLogMessage << "Apply ProjectGn to deformed matrix" << std::endl;
     Sp<Nc>::ProjectGn(U);
-    aux = U*adj(U) - identity;
-    std::cout << GridLogMessage << "U adjU - 1 = " << norm2(aux) << std::endl;
-    assert( norm2(aux) < 1e-8);
-    std::cout << GridLogMessage << "Det after ProjectGn = " << norm2( Determinant(U) ) / vol << std::endl;
-    assert( norm2(aux) - 1 < 1e-8);
-    Sp<Nc>::OmegaInvariance(U);
-    
-    std::cout << GridLogMessage << "Checking the structure is " << std::endl;
-    std::cout << GridLogMessage << "U  =  (   W    X   )  " << std::endl;
-    std::cout << GridLogMessage << "      (  -X^*  W^* )  " << std::endl;
-    std::cout <<GridLogMessage << std::endl;
-    for (int c1 = 0; c1 < nsp; c1++) //check on W
-    {
-        for (int c2 = 0; c2 < nsp; c2++)
-        {
-            auto W = PeekIndex<ColourIndex>(U,c1,c2);
-            auto Wstar =  PeekIndex<ColourIndex>(U,c1+nsp,c2+nsp);
-            auto Ww = conjugate( Wstar );
-            auto amizero = sum(W - Ww);
-            auto amizeroo = TensorRemove(amizero);
-            assert(  amizeroo.real() < 10e-6 );
-            amizeroo *= i;
-            assert(  amizeroo.real() < 10e-6 );
-        }
-        
-    }
-
-    for (int c1 = 0; c1 < nsp ; c1++)
-    {
-        for (int c2 = 0; c2 < nsp; c2++)
-        {
-            auto X = PeekIndex<ColourIndex>(U,c1,c2+nsp);
-            auto minusXstar = PeekIndex<ColourIndex>(U,c1+nsp,c2);
-            auto minusXx = conjugate(minusXstar);
-            auto amizero = sum (X + minusXx);
-            auto amizeroo = TensorRemove(amizero);
-            assert(  amizeroo.real() < 10e-6 );
-            amizeroo *= i;
-            assert(  amizeroo.real() < 10e-6 );
-        }
-    }
-    std::cout <<GridLogMessage << std::endl;
-    std::cout << GridLogMessage << "#   #   #   #" << std::endl;
-    std::cout << GridLogMessage << "Algebra" << std::endl;
-    std::cout << GridLogMessage << "#   #   #   #" << std::endl;
-    std::cout <<GridLogMessage << std::endl;
-    std::cout << GridLogMessage << "Testing SpTa" << std::endl;
+    std::cout << GridLogMessage << "Run checks:" << std::endl;
+    SpGroupQuiz(U);
+    det =  sum( Determinant(U) );
+    std::cout << GridLogMessage << "Re(Det) after ProjectGn (constrained) = " << real(det) / vol << std::endl;
+    assert ( (real(det) / vol) - 1 < 1e-8 );
+    det = det*i;
+    std::cout << GridLogMessage << "Im(Det) after ProjectGn (constrained) = " << real(det) / vol << std::endl;
+    assert ( real(det) / vol < 1e-8 );
     
+    //  SpTa
     SU<Nc>::HotConfiguration(pRNG,Umu);//refresh
-    U = PeekIndex<LorentzIndex>(Umu,2);
-
-    U = U + Delta*identity;
-    std::cout << GridLogMessage << "Matrix deformed " << std::endl;
+    U = PeekIndex<LorentzIndex>(Umu,0);
+    std::cout << GridLogMessage << "Testing SpTa" << std::endl;
+    U = U + epsilon*identity - i*identity;
     std::cout << GridLogMessage << "Apply SpTa to deformed matrix" << std::endl;
     U = SpTa(U);
-    
-    aux = U - adj(U);
-    std::cout << GridLogMessage << "SpTa ::: T - Tda = " << norm2(aux) << std::endl;
-    aux = U + adj(U);
-    std::cout << GridLogMessage << "SpTa ::: T + Tda = " << norm2(aux) << std::endl;
-    assert( norm2(aux) - 1 < 1e-8);
-    
-    std::cout << GridLogMessage << "Check that Omega T Omega + conj(T) = 0 " << std::endl;
-
-    LatticeColourMatrixD Omega(&Grid);
-    Sp<Nc>::Omega(Omega);
-    aux = Omega*U*Omega + conjugate(U);
-    assert( norm2(aux) < 1e-8);
-    
-    
-    std::cout << GridLogMessage << "Checking the structure is " << std::endl;
-    std::cout << GridLogMessage << "U  =  (   W    X   )  " << std::endl;
-    std::cout << GridLogMessage << "      (  -X^*  W^* )  " << std::endl;
-    std::cout <<GridLogMessage << std::endl;
-    for (int c1 = 0; c1 < nsp; c1++) //check on W
-    {
-        for (int c2 = 0; c2 < nsp; c2++)
-        {
-            auto W = PeekIndex<ColourIndex>(U,c1,c2);
-            auto Wstar =  PeekIndex<ColourIndex>(U,c1+nsp,c2+nsp);
-            auto Ww = conjugate( Wstar );
-            auto amizero = sum(W - Ww);
-            auto amizeroo = TensorRemove(amizero);
-            assert(  amizeroo.real() < 10e-6 );
-            amizeroo *= i;
-            assert(  amizeroo.real() < 10e-6 );
-        }
-        
-    }
-
-    for (int c1 = 0; c1 < nsp ; c1++)
-    {
-        for (int c2 = 0; c2 < nsp; c2++)
-        {
-            auto X = PeekIndex<ColourIndex>(U,c1,c2+nsp);
-            auto minusXstar = PeekIndex<ColourIndex>(U,c1+nsp,c2);
-            auto minusXx = conjugate(minusXstar);
-            auto amizero = sum (X + minusXx);
-            auto amizeroo = TensorRemove(amizero);
-            assert(  amizeroo.real() < 10e-6 );
-            amizeroo *= i;
-            assert(  amizeroo.real() < 10e-6 );
-        }
-    }
+    std::cout << GridLogMessage << "Run checks:" << std::endl;
+    SpAntiHermitianAlgebraQuiz(U);
         
     Grid_finalize();