Added support for the Two index Symmetric and Antisymmetric representations

Tested for HMC convergence: OK
Added also a test file showing an example for mixed representations

lib/qcd/QCD.h

@@ -45,7 +45,7 @@ namespace QCD {
     static const int Zm = 6;
     static const int Tm = 7;
 
-    static const int Nc=3;
+    static const int Nc=2;
     static const int Ns=4;
     static const int Nd=4;
     static const int Nhs=2; // half spinor
@@ -500,9 +500,10 @@ namespace QCD {
 #include <Grid/qcd/utils/LinalgUtils.h>
 #include <Grid/qcd/utils/CovariantCshift.h>
 
+// Include representations 	
 #include <Grid/qcd/utils/SUn.h>
 #include <Grid/qcd/utils/SUnAdjoint.h>
-	
+#include <Grid/qcd/utils/SUnTwoIndex.h>
 #include <Grid/qcd/representations/hmc_types.h>
 
 #include <Grid/qcd/action/Actions.h>

lib/qcd/action/Actions.h

@@ -120,6 +120,10 @@ typedef SymanzikGaugeAction<ConjugateGimplD>        ConjugateSymanzikGaugeAction
   template class A<WilsonAdjImplF>; \
   template class A<WilsonAdjImplD>; 
 
+#define TwoIndexFermOpTemplateInstantiate(A) \
+  template class A<WilsonTwoIndexSymmetricImplF>; \
+  template class A<WilsonTwoIndexSymmetricImplD>; 
+
 #define FermOp5dVecTemplateInstantiate(A) \
   template class A<DomainWallVec5dImplF>;	\
   template class A<DomainWallVec5dImplD>;	\
@@ -180,6 +184,10 @@ typedef WilsonFermion<WilsonAdjImplR> WilsonAdjFermionR;
 typedef WilsonFermion<WilsonAdjImplF> WilsonAdjFermionF;
 typedef WilsonFermion<WilsonAdjImplD> WilsonAdjFermionD;
 
+typedef WilsonFermion<WilsonTwoIndexSymmetricImplR> WilsonTwoIndexSymmetricFermionR;
+typedef WilsonFermion<WilsonTwoIndexSymmetricImplF> WilsonTwoIndexSymmetricFermionF;
+typedef WilsonFermion<WilsonTwoIndexSymmetricImplD> WilsonTwoIndexSymmetricFermionD;
+
 typedef WilsonTMFermion<WilsonImplR> WilsonTMFermionR;
 typedef WilsonTMFermion<WilsonImplF> WilsonTMFermionF;
 typedef WilsonTMFermion<WilsonImplD> WilsonTMFermionD;

lib/qcd/action/fermion/FermionOperatorImpl.h

@@ -523,6 +523,10 @@ namespace Grid {
     typedef WilsonImpl<vComplexF, AdjointRepresentation > WilsonAdjImplF;  // Float
     typedef WilsonImpl<vComplexD, AdjointRepresentation > WilsonAdjImplD;  // Double
 
+    typedef WilsonImpl<vComplex,  TwoIndexSymmetricRepresentation > WilsonTwoIndexSymmetricImplR;   // Real.. whichever prec
+    typedef WilsonImpl<vComplexF, TwoIndexSymmetricRepresentation > WilsonTwoIndexSymmetricImplF;  // Float
+    typedef WilsonImpl<vComplexD, TwoIndexSymmetricRepresentation > WilsonTwoIndexSymmetricImplD;  // Double
+
     typedef DomainWallVec5dImpl<vComplex ,Nc> DomainWallVec5dImplR; // Real.. whichever prec
     typedef DomainWallVec5dImpl<vComplexF,Nc> DomainWallVec5dImplF; // Float
     typedef DomainWallVec5dImpl<vComplexD,Nc> DomainWallVec5dImplD; // Double

lib/qcd/action/fermion/WilsonFermion.cc

@@ -309,6 +309,7 @@ void WilsonFermion<Impl>::DhopInternal(StencilImpl &st, LebesgueOrder &lo,
 
 FermOpTemplateInstantiate(WilsonFermion);
 AdjointFermOpTemplateInstantiate(WilsonFermion);
+TwoIndexFermOpTemplateInstantiate(WilsonFermion);
 GparityFermOpTemplateInstantiate(WilsonFermion);
 }
 }

lib/qcd/action/fermion/WilsonKernels.cc

@@ -526,6 +526,7 @@ void WilsonKernels<Impl>::DiracOptDhopDir(
 
 FermOpTemplateInstantiate(WilsonKernels);
 AdjointFermOpTemplateInstantiate(WilsonKernels);
+TwoIndexFermOpTemplateInstantiate(WilsonKernels);
 
 }}
 

lib/qcd/representations/hmc_types.h

@@ -2,6 +2,7 @@
 #define HMC_TYPES_H
 
 #include <Grid/qcd/representations/adjoint.h>
+#include <Grid/qcd/representations/two_index.h>
 #include <Grid/qcd/representations/fundamental.h>
 #include <tuple>
 #include <utility>

lib/qcd/representations/two_index.h

@@ -1,22 +1,28 @@
 /*
  *  Policy classes for the HMC
- *  Author: Guido Cossu
+ *  Authors: Guido Cossu, David Preti
 */
 
-#ifndef ADJOINT_H
+#ifndef SUN2INDEX_H_H
-#define ADJOINT_H
+#define SUN2INDEX_H_H
 
 namespace Grid {
 namespace QCD {
 
 /*
  * This is an helper class for the HMC
-* Should contain only the data for the adjoint representation
+ * Should contain only the data for the two index representations
-* and the facility to convert from the fundamental -> adjoint
+ * and the facility to convert from the fundamental -> two index
+ * The templated parameter TwoIndexSymmetry choses between the 
+ * symmetric and antisymmetric representations
+ * 
+ * There is an 
+ * enum TwoIndexSymmetry { Symmetric = 1, AntiSymmetric = -1 };
+ * in the SUnTwoIndex.h file
  */
 
 template <int ncolour, TwoIndexSymmetry S>
-class TwoIndexSymmetricRep {
+class TwoIndexRep {
  public:
   // typdef to be used by the Representations class in HMC to get the
   // types for the higher representation fields
@@ -26,28 +32,31 @@ class TwoIndexSymmetricRep {
 
   LatticeField U;
 
-  explicit TwoIndexSymmetricRep(GridBase *grid) : U(grid) {}
+  explicit TwoIndexRep(GridBase *grid) : U(grid) {}
 
   void update_representation(const LatticeGaugeField &Uin) {
     std::cout << GridLogDebug << "Updating TwoIndex representation\n";
     // Uin is in the fundamental representation
-    // get the U in AdjointRep
+    // get the U in TwoIndexRep
-    // (U)(ij)_(lk) = 
+    // (U)_{(ij)(lk)} = tr [ adj(e^(ij)) U e^(lk) transpose(U) ]
-    // e^a = 
     conformable(U, Uin);
     U = zero;
     LatticeColourMatrix tmp(Uin._grid);
 
-    Vector<typename SU<ncolour>::Matrix> ta(Dimension);
+    Vector<typename SU<ncolour>::Matrix> eij(Dimension);
 
-    // FIXME probably not very efficient to get all the generators
+    for (int a = 0; a < Dimension; a++)
-    // everytime
+      SU_TwoIndex<ncolour, S>::base(a, eij[a]);
-    for (int a = 0; a < Dimension; a++) SU<ncolour>::generator(a, ta[a]);
 
     for (int mu = 0; mu < Nd; mu++) {
       auto Uin_mu = peekLorentz(Uin, mu);
       auto U_mu = peekLorentz(U, mu);
-
+      for (int a = 0; a < Dimension; a++) {
+        tmp = transpose(Uin_mu) * adj(eij[a]) * Uin_mu;
+        for (int b = 0; b < Dimension; b++)
+          pokeColour(U_mu, trace(tmp * eij[b]), a, b);
+      }
+      pokeLorentz(U, U_mu, mu);
     }
   }
 
@@ -87,5 +96,4 @@ class TwoIndexSymmetricRep {
 typedef TwoIndexRep<Nc, AntiSymmetric> TwoIndexAntiSymmetricRepresentation;
 }
 }
-
 #endif

lib/qcd/utils/SUnTwoIndex.h

@@ -1,6 +1,3 @@
-#ifndef QCD_UTIL_SUNADJOINT_H
-#define QCD_UTIL_SUNADJOINT_H
-
 ////////////////////////////////////////////////////////////////////////
 //
 // * Two index representation generators
@@ -17,23 +14,33 @@
 //
 //   Then the generators are written as
 //
-//   (iT^(ij))_lk = i 
+//   (iT_a)^(ij)(lk) = i * ( tr[e^(ij)^dag e^(lk) T^trasp_a] +
+//   tr[e^(lk)e^(ij)^dag T_a] )  //
+//   
 //
 ////////////////////////////////////////////////////////////////////////
 
+// Authors: David Preti, Guido Cossu
+
+#ifndef QCD_UTIL_SUN2INDEX_H
+#define QCD_UTIL_SUN2INDEX_H
+
+
 namespace Grid {
 namespace QCD {
 
 enum TwoIndexSymmetry { Symmetric = 1, AntiSymmetric = -1 };
 
+inline Real delta(int a, int b) { return (a == b) ? 1.0 : 0.0; }
+
 template <int ncolour, TwoIndexSymmetry S>
 class SU_TwoIndex : public SU<ncolour> {
  public:
   static const int Dimension = ncolour * (ncolour + S) / 2;
+  static const int NumGenerators = SU<ncolour>::AdjointDimension;
 
   template <typename vtype>
-      using iSUnTwoIndexMatrix =
+  using iSUnTwoIndexMatrix = iScalar<iScalar<iMatrix<vtype, Dimension> > >;
-	iScalar<iScalar<iMatrix<vtype, Dimension > > >;
 
   typedef iSUnTwoIndexMatrix<Complex> TIMatrix;
   typedef iSUnTwoIndexMatrix<ComplexF> TIMatrixF;
@@ -43,9 +50,9 @@ namespace Grid {
   typedef iSUnTwoIndexMatrix<vComplexF> vTIMatrixF;
   typedef iSUnTwoIndexMatrix<vComplexD> vTIMatrixD;
 
-      typedef Lattice<vAMatrix>  LatticeTwoIndexMatrix;
+  typedef Lattice<vTIMatrix> LatticeTwoIndexMatrix;
-      typedef Lattice<vAMatrixF> LatticeTwoIndexMatrixF;
+  typedef Lattice<vTIMatrixF> LatticeTwoIndexMatrixF;
-      typedef Lattice<vAMatrixD> LatticeTwoIndexMatrixD;
+  typedef Lattice<vTIMatrixD> LatticeTwoIndexMatrixD;
 
   typedef Lattice<iVector<iScalar<iMatrix<vComplex, Dimension> >, Nd> >
       LatticeTwoIndexField;
@@ -54,108 +61,213 @@ namespace Grid {
   typedef Lattice<iVector<iScalar<iMatrix<vComplexD, Dimension> >, Nd> >
       LatticeTwoIndexFieldD;
 
+  template <typename vtype>
+  using iSUnMatrix = iScalar<iScalar<iMatrix<vtype, ncolour> > >;
 
-
+  typedef iSUnMatrix<Complex> Matrix;
+  typedef iSUnMatrix<ComplexF> MatrixF;
+  typedef iSUnMatrix<ComplexD> MatrixD;
 
   template <class cplx>
-      static void generator(int Index, iSUnTwoIndexMatrix<cplx> &iTwoIdxTa) {
+  static void base(int Index, iSUnMatrix<cplx> &eij) {
-	// returns i(T)^(ij) necessary for the projectors
+    // returns (e)^(ij)_{kl} necessary for change of base U_F -> U_R
-	// see definitions above
+    assert(Index < NumGenerators);
-	iTwoIdxTa = zero;
+    eij = zero;
-	Vector<typename SU<ncolour>::template iSUnMatrix<cplx> > tij(Dimension);
-	typename SU<ncolour>::template iSUnMatrix<cplx> tmp;
 
+    // for the linearisation of the 2 indexes 
+    static int a[ncolour * (ncolour - 1) / 2][2]; // store the a <-> i,j
+    static bool filled = false;
+    if (!filled) {
+      int counter = 0;
+      for (int i = 1; i < ncolour; i++) {
+        for (int j = 0; j < i; j++) {
+          a[counter++][0] = i;
+          a[counter++][1] = j;
+        }
+      }
+      filled = true;
+    }
+
+    if (Index < ncolour * (ncolour - 1) / 2) {
+      baseOffDiagonal(a[Index][0], a[Index][1], eij);
+    } else {
+      baseDiagonal(Index, eij);
+    }
+  }
+
+  template <class cplx>
+  static void baseDiagonal(int Index, iSUnMatrix<cplx> &eij) {
+    eij = zero;
+    eij()()(Index - ncolour * (ncolour - 1) / 2,
+            Index - ncolour * (ncolour - 1) / 2) = 1.0;
+  }
+
+  template <class cplx>
+  static void baseOffDiagonal(int i, int j, iSUnMatrix<cplx> &eij) {
+    eij = zero;
+    for (int k = 0; k < ncolour; k++)
+      for (int l = 0; l < ncolour; l++)
+        eij()()(l, k) = delta(i, k) * delta(j, l) +
+                        S * delta(j, k) * delta(i, l);
+
+    RealD nrm = 1. / std::sqrt(2.0);
+    eij = eij * nrm;
+  }
+
+  static void printBase(void) {
+    for (int gen = 0; gen < Dimension; gen++) {
+      Matrix tmp;
+      base(gen, tmp);
+      std::cout << GridLogMessage << "Nc = " << ncolour << " t_" << gen
+                << std::endl;
+      std::cout << GridLogMessage << tmp << std::endl;
+    }
+  }
+
+  template <class cplx>
+  static void generator(int Index, iSUnTwoIndexMatrix<cplx> &i2indTa) {
+    Vector<typename SU<ncolour>::template iSUnMatrix<cplx> > ta(
+        ncolour * ncolour - 1);
+    Vector<typename SU<ncolour>::template iSUnMatrix<cplx> > eij(Dimension);
+    typename SU<ncolour>::template iSUnMatrix<cplx> tmp;
+    i2indTa = zero;
+    
+    for (int a = 0; a < ncolour * ncolour - 1; a++)
+      SU<ncolour>::generator(a, ta[a]);
+    
+    for (int a = 0; a < Dimension; a++) base(a, eij[a]);
 
     for (int a = 0; a < Dimension; a++) {
-
+      tmp = transpose(ta[Index]) * adj(eij[a]) + adj(eij[a]) * ta[Index];
+      for (int b = 0; b < Dimension; b++) {
+        typename SU<ncolour>::template iSUnMatrix<cplx> tmp1 =
+            tmp * eij[b]; 
+        Complex iTr = TensorRemove(timesI(trace(tmp1)));
+        i2indTa()()(a, b) = iTr;
+      }
     }
   }
 
   static void printGenerators(void) {
-	for (int gen = 0; gen < Dimension; gen++) {
+    for (int gen = 0; gen < ncolour * ncolour - 1; gen++) {
-	  AMatrix ta;
+      TIMatrix i2indTa;
-	  generator(gen, ta);
+      generator(gen, i2indTa);
       std::cout << GridLogMessage << "Nc = " << ncolour << " t_" << gen
                 << std::endl;
-	  std::cout << GridLogMessage << ta << std::endl;
+      std::cout << GridLogMessage << i2indTa << std::endl;
     }
   }
 
   static void testGenerators(void) {
-	TIMatrix TwoIndexTa;
+    TIMatrix i2indTa, i2indTb;
+    std::cout << GridLogMessage << "2IndexRep - Checking if traceless"
+              << std::endl;
+    for (int a = 0; a < ncolour * ncolour - 1; a++) {
+      generator(a, i2indTa);
+      std::cout << GridLogMessage << a << std::endl;
+      assert(norm2(trace(i2indTa)) < 1.0e-6);
+    }
+    std::cout << GridLogMessage << std::endl;
 
+    std::cout << GridLogMessage << "2IndexRep - Checking if antihermitean"
+              << std::endl;
+    for (int a = 0; a < ncolour * ncolour - 1; a++) {
+      generator(a, i2indTa);
+      std::cout << GridLogMessage << a << std::endl;
+      assert(norm2(adj(i2indTa) + i2indTa) < 1.0e-6);
     }
 
-      static void TwoIndexLieAlgebraMatrix(const typename SU<ncolour>::LatticeAlgebraVector &h,
+    std::cout << GridLogMessage << std::endl;
+    std::cout << GridLogMessage
+              << "2IndexRep - Checking Tr[Ta*Tb]=delta(a,b)*(N +- 2)/2"
+              << std::endl;
+    for (int a = 0; a < ncolour * ncolour - 1; a++) {
+      for (int b = 0; b < ncolour * ncolour - 1; b++) {
+        generator(a, i2indTa);
+        generator(b, i2indTb);
+
+        // generator returns iTa, so we need a minus sign here
+        Complex Tr = -TensorRemove(trace(i2indTa * i2indTb));
+        std::cout << GridLogMessage << "a=" << a << "b=" << b << "Tr=" << Tr
+                  << std::endl;
+      }
+    }
+    std::cout << GridLogMessage << std::endl;
+  }
+
+  static void TwoIndexLieAlgebraMatrix(
+      const typename SU<ncolour>::LatticeAlgebraVector &h,
       LatticeTwoIndexMatrix &out, Real scale = 1.0) {
     conformable(h, out);
     GridBase *grid = out._grid;
-	LatticeAdjMatrix la(grid);
+    LatticeTwoIndexMatrix la(grid);
-	TIMatrix iTa;
+    TIMatrix i2indTa;
 
     out = zero;
-	for (int a = 0; a < Dimension; a++) {
+    for (int a = 0; a < ncolour * ncolour - 1; a++) {
-	  generator(a, iTa);
+      generator(a, i2indTa);
-	  la = peekColour(h, a) * iTa;
+      la = peekColour(h, a) * i2indTa;
       out += la;
     }
     out *= scale;
   }
 
-      // Projects the algebra components a lattice matrix (of dimension ncol*ncol -1 )
+  // Projects the algebra components 
-      static void projectOnAlgebra(typename SU<ncolour>::LatticeAlgebraVector &h_out, const LatticeTwoIndexMatrix &in, Real scale = 1.0) {
+  // of a lattice matrix ( of dimension ncol*ncol -1 )
+  static void projectOnAlgebra(
+      typename SU<ncolour>::LatticeAlgebraVector &h_out,
+      const LatticeTwoIndexMatrix &in, Real scale = 1.0) {
     conformable(h_out, in);
     h_out = zero;
-	TIMatrix iTa;
+    TIMatrix i2indTa;
-	Real coefficient = - 2.0/(ncolour + 2*S) * scale;// 2/(Nc +/- 2) for the normalization of the trace in the two index rep
+    Real coefficient = -2.0 / (ncolour + 2 * S) * scale;
-
+    // 2/(Nc +/- 2) for the normalization of the trace in the two index rep
-	for (int a = 0; a < Dimension; a++) {
+    for (int a = 0; a < ncolour * ncolour - 1; a++) {
-	  generator(a, iTa);
+      generator(a, i2indTa);
-	  auto tmp = real(trace(iTa * in)) * coefficient;
+      auto tmp = real(trace(i2indTa * in)) * coefficient;
       pokeColour(h_out, tmp, a);
     }
   }
 
-      // a projector that keeps the generators stored to avoid the overhead of recomputing them 
+  // a projector that keeps the generators stored to avoid the overhead of
-      static void projector(typename SU<ncolour>::LatticeAlgebraVector &h_out, const LatticeTwoIndexMatrix &in, Real scale = 1.0) {
+  // recomputing them
+  static void projector(typename SU<ncolour>::LatticeAlgebraVector &h_out,
+                        const LatticeTwoIndexMatrix &in, Real scale = 1.0) {
     conformable(h_out, in);
-	static std::vector<TIMatrix> iTa(Dimension);  // to store the generators
+    // to store the generators
+    static std::vector<TIMatrix> i2indTa(ncolour * ncolour -1); 
     h_out = zero;
     static bool precalculated = false;
     if (!precalculated) {
       precalculated = true;
-	  for (int a = 0; a < Dimension; a++) generator(a, iTa[a]);
+      for (int a = 0; a < ncolour * ncolour - 1; a++) generator(a, i2indTa[a]);
     }
 
-	Real coefficient = - 2.0/(ncolour + 2*S) * scale; // 2/(Nc +/- 2) for the normalization of the trace in the two index rep
+    Real coefficient =
+        -2.0 / (ncolour + 2 * S) * scale;  // 2/(Nc +/- 2) for the normalization
+                                           // of the trace in the two index rep
 
-	for (int a = 0; a < Dimension; a++) {
+    for (int a = 0; a < ncolour * ncolour - 1; a++) {
-	  auto tmp = real(trace(iTa[a] * in)) * coefficient; 
+      auto tmp = real(trace(i2indTa[a] * in)) * coefficient;
       pokeColour(h_out, tmp, a);
     }
   }
-
-
 };
 
-
-
-
 // Some useful type names
+typedef SU_TwoIndex<Nc, Symmetric> TwoIndexSymmMatrices;
+typedef SU_TwoIndex<Nc, AntiSymmetric> TwoIndexAntiSymmMatrices;
+
 typedef SU_TwoIndex<2, Symmetric> SU2TwoIndexSymm;
 typedef SU_TwoIndex<3, Symmetric> SU3TwoIndexSymm;
 typedef SU_TwoIndex<4, Symmetric> SU4TwoIndexSymm;
 typedef SU_TwoIndex<5, Symmetric> SU5TwoIndexSymm;
 
-    typedef SU_TwoIndex<Nc, Symmetric> TwoIndexSymmMatrices;
-
 typedef SU_TwoIndex<2, AntiSymmetric> SU2TwoIndexAntiSymm;
 typedef SU_TwoIndex<3, AntiSymmetric> SU3TwoIndexAntiSymm;
 typedef SU_TwoIndex<4, AntiSymmetric> SU4TwoIndexAntiSymm;
 typedef SU_TwoIndex<5, AntiSymmetric> SU5TwoIndexAntiSymm;
 
-    typedef SU_TwoIndex<Nc, AntiSymmetric> TwoIndexAntiSymmMatrices;
-
 
 }
 }

tests/core/Test_lie_generators.cc

@@ -8,6 +8,7 @@ Copyright (C) 2015
 
 Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: Guido Cossu <guido.cossu@ed.ac.uk>
 
 This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
@@ -30,9 +31,13 @@ directory
 #include <Grid/Grid.h>
 
 #include <Grid/qcd/utils/CovariantCshift.h>
+
 #include <Grid/qcd/utils/SUn.h>
 #include <Grid/qcd/utils/SUnAdjoint.h>
+#include <Grid/qcd/utils/SUnTwoIndex.h>
+
 #include <Grid/qcd/representations/adjoint.h>
+#include <Grid/qcd/representations/two_index.h>
 #include <Grid/qcd/utils/WilsonLoops.h>
 
 using namespace std;
@@ -79,13 +84,6 @@ int main(int argc, char** argv) {
   SU4::testGenerators();
   SU4Adjoint::testGenerators();
 
-  //  std::cout<<GridLogMessage<<"*********************************************"<<std::endl;
-  //  std::cout<<GridLogMessage<<"* Generators for SU(5)"<<std::endl;
-  //  std::cout<<GridLogMessage<<"*********************************************"<<std::endl;
-  //  SU5::printGenerators();
-  //  SU5::testGenerators();
-
-
   // Projectors 
   GridParallelRNG gridRNG(grid);
   gridRNG.SeedRandomDevice();
@@ -112,12 +110,14 @@ int main(int argc, char** argv) {
   AdjointRep<Nc> AdjRep(grid);
 
   // AdjointRepresentation has the predefined number of colours Nc
-  Representations<FundamentalRepresentation, AdjointRepresentation> RepresentationTypes(grid);  
+  Representations<FundamentalRepresentation, AdjointRepresentation, TwoIndexSymmetricRepresentation> RepresentationTypes(grid);  
+
   
   LatticeGaugeField U(grid), V(grid);
   SU<Nc>::HotConfiguration<LatticeGaugeField>(gridRNG, U);
   SU<Nc>::HotConfiguration<LatticeGaugeField>(gridRNG, V);
 
+  // Adjoint representation
   // Test group structure
   // (U_f * V_f)_r = U_r * V_r
   LatticeGaugeField UV(grid);
@@ -147,7 +147,7 @@ int main(int argc, char** argv) {
   }
 
   typename AdjointRep<Nc>::LatticeField Diff_check = UVr - UrVr;
-  std::cout << GridLogMessage << "Group structure SU("<<Nc<<") check difference : " << norm2(Diff_check) << std::endl;
+  std::cout << GridLogMessage << "Group structure SU("<<Nc<<") check difference (Adjoint representation) : " << norm2(Diff_check) << std::endl;
 
   // Check correspondence of algebra and group transformations
   // Create a random vector
@@ -161,7 +161,7 @@ int main(int argc, char** argv) {
   SU<Nc>::LatticeAlgebraVector h_adj2(grid);
   SU_Adjoint<Nc>::projectOnAlgebra(h_adj2, Ar);
   SU<Nc>::LatticeAlgebraVector h_diff = h_adj - h_adj2;
-  std::cout << GridLogMessage << "Projections structure check vector difference : " << norm2(h_diff) << std::endl;
+  std::cout << GridLogMessage << "Projections structure check vector difference (Adjoint representation) : " << norm2(h_diff) << std::endl;
 
   // Exponentiate
   typename AdjointRep<Nc>::LatticeMatrix Uadj(grid);
@@ -210,5 +210,323 @@ int main(int argc, char** argv) {
   typename AdjointRep<Nc>::LatticeMatrix Diff_check_mat = Ur0 - Uadj;
   std::cout << GridLogMessage << "Projections structure check group difference : " << norm2(Diff_check_mat) << std::endl;
 
+
+
+
+  // TwoIndexRep tests
+
+  std::cout << GridLogMessage << "*********************************************"
+            << std::endl;
+  std::cout << GridLogMessage << "*********************************************"
+            << std::endl;
+
+
+  
+  std::cout << GridLogMessage << "* eS^{ij} base for SU(2)" << std::endl;
+  std::cout << GridLogMessage << "*********************************************"
+            << std::endl;
+  std::cout << GridLogMessage << "Dimension of Two Index Symmetric representation: "<< SU2TwoIndexSymm::Dimension << std::endl;
+  SU2TwoIndexSymm::printBase();
+      std::cout << GridLogMessage << "*********************************************"
+            << std::endl;
+        std::cout << GridLogMessage << "Generators of Two Index Symmetric representation: "<< SU2TwoIndexSymm::Dimension << std::endl;
+  SU2TwoIndexSymm::printGenerators();
+        std::cout << GridLogMessage << "Test of Two Index Symmetric Generators: "<< SU2TwoIndexSymm::Dimension << std::endl;
+  SU2TwoIndexSymm::testGenerators();
+  std::cout << GridLogMessage << "*********************************************"
+            << std::endl;
+
+
+  
+  std::cout << GridLogMessage << "*********************************************"
+            << std::endl;
+  std::cout << GridLogMessage << "* eAS^{ij} base for SU(2)" << std::endl;
+  
+  std::cout << GridLogMessage << "*********************************************"
+            << std::endl;
+  std::cout << GridLogMessage << "Dimension of Two Index anti-Symmetric representation: "<< SU2TwoIndexAntiSymm::Dimension << std::endl;
+  SU2TwoIndexAntiSymm::printBase();
+      std::cout << GridLogMessage << "*********************************************"
+            << std::endl;
+        std::cout << GridLogMessage << "Dimension of Two Index anti-Symmetric representation: "<< SU2TwoIndexAntiSymm::Dimension << std::endl;
+  SU2TwoIndexAntiSymm::printGenerators();
+  std::cout << GridLogMessage << "Test of Two Index anti-Symmetric Generators: "<< SU2TwoIndexAntiSymm::Dimension << std::endl;
+  SU2TwoIndexAntiSymm::testGenerators();
+  
+  
+  
+  std::cout << GridLogMessage << "*********************************************"
+      << std::endl;
+  std::cout << GridLogMessage << "Test for the Two Index Symmetric projectors"
+      << std::endl;
+  // Projectors 
+  SU3TwoIndexSymm::LatticeTwoIndexMatrix Gauss2(grid);
+  random(gridRNG,Gauss2);
+  
+  std::cout << GridLogMessage << "Start projectOnAlgebra" << std::endl;
+  SU3TwoIndexSymm::projectOnAlgebra(ha, Gauss2);
+  std::cout << GridLogMessage << "end projectOnAlgebra" << std::endl;
+  std::cout << GridLogMessage << "Start projector" << std::endl;
+  SU3TwoIndexSymm::projector(hb, Gauss2);
+  std::cout << GridLogMessage << "end projector" << std::endl;
+  
+  std::cout << GridLogMessage << "ReStart projector" << std::endl;
+  SU3TwoIndexSymm::projector(hb, Gauss2);
+  std::cout << GridLogMessage << "end projector" << std::endl;
+  SU3::LatticeAlgebraVector diff2 = ha - hb;
+  std::cout << GridLogMessage << "Difference: " << norm2(diff) << std::endl;
+  std::cout << GridLogMessage << "*********************************************"
+      << std::endl;
+
+  
+    std::cout << GridLogMessage << "*********************************************"
+      << std::endl;
+  std::cout << GridLogMessage << "Test for the Two index anti-Symmetric projectors"
+      << std::endl;
+  // Projectors
+  SU3TwoIndexAntiSymm::LatticeTwoIndexMatrix Gauss2a(grid);
+  random(gridRNG,Gauss2a);
+  
+  std::cout << GridLogMessage << "Start projectOnAlgebra" << std::endl;
+  SU3TwoIndexAntiSymm::projectOnAlgebra(ha, Gauss2a);
+  std::cout << GridLogMessage << "end projectOnAlgebra" << std::endl;
+  std::cout << GridLogMessage << "Start projector" << std::endl;
+  SU3TwoIndexAntiSymm::projector(hb, Gauss2a);
+  std::cout << GridLogMessage << "end projector" << std::endl;
+  
+  std::cout << GridLogMessage << "ReStart projector" << std::endl;
+  SU3TwoIndexAntiSymm::projector(hb, Gauss2a);
+  std::cout << GridLogMessage << "end projector" << std::endl;
+  SU3::LatticeAlgebraVector diff2a = ha - hb;
+  std::cout << GridLogMessage << "Difference: " << norm2(diff2a) << std::endl;
+  std::cout << GridLogMessage << "*********************************************"
+      << std::endl;
+    
+  
+  std::cout << GridLogMessage << "Two index Symmetric: Checking Group Structure"
+      << std::endl;
+  // Testing HMC representation classes
+  TwoIndexRep< Nc, Symmetric > TIndexRep(grid);
+
+  // Test group structure
+  // (U_f * V_f)_r = U_r * V_r
+  LatticeGaugeField U2(grid), V2(grid);
+  SU<Nc>::HotConfiguration<LatticeGaugeField>(gridRNG, U2);
+  SU<Nc>::HotConfiguration<LatticeGaugeField>(gridRNG, V2);
+  
+  LatticeGaugeField UV2(grid);
+  UV2 = zero;
+  for (int mu = 0; mu < Nd; mu++) {
+    SU<Nc>::LatticeMatrix Umu2 = peekLorentz(U2,mu);
+    SU<Nc>::LatticeMatrix Vmu2 = peekLorentz(V2,mu);
+    pokeLorentz(UV2,Umu2*Vmu2, mu);
+  }
+  
+  TIndexRep.update_representation(UV2);
+  typename TwoIndexRep< Nc, Symmetric >::LatticeField UVr2 = TIndexRep.U;  // (U_f * V_f)_r
+  
+  TIndexRep.update_representation(U2);
+  typename TwoIndexRep< Nc, Symmetric >::LatticeField Ur2 = TIndexRep.U;  // U_r
+  
+  TIndexRep.update_representation(V2);
+  typename TwoIndexRep< Nc, Symmetric >::LatticeField Vr2 = TIndexRep.U;  // V_r
+  
+  typename TwoIndexRep< Nc, Symmetric >::LatticeField Ur2Vr2(grid);
+  Ur2Vr2 = zero;
+  for (int mu = 0; mu < Nd; mu++) {
+    typename TwoIndexRep< Nc, Symmetric >::LatticeMatrix Urmu2 = peekLorentz(Ur2,mu);
+    typename TwoIndexRep< Nc, Symmetric >::LatticeMatrix Vrmu2 = peekLorentz(Vr2,mu);
+    pokeLorentz(Ur2Vr2,Urmu2*Vrmu2, mu);
+  }
+  
+  typename TwoIndexRep< Nc, Symmetric >::LatticeField Diff_check2 = UVr2 - Ur2Vr2;
+  std::cout << GridLogMessage << "Group structure SU("<<Nc<<") check difference (Two Index Symmetric): " << norm2(Diff_check2) << std::endl;
+
+  
+  // Check correspondence of algebra and group transformations
+  // Create a random vector
+  SU<Nc>::LatticeAlgebraVector h_sym(grid);
+  typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix Ar_sym(grid);
+  random(gridRNG,h_sym);
+  h_sym = real(h_sym);
+  SU_TwoIndex<Nc,Symmetric>::TwoIndexLieAlgebraMatrix(h_sym,Ar_sym);
+  
+  // Re-extract h_sym
+  SU<Nc>::LatticeAlgebraVector h_sym2(grid);
+  SU_TwoIndex< Nc, Symmetric>::projectOnAlgebra(h_sym2, Ar_sym);
+  SU<Nc>::LatticeAlgebraVector h_diff_sym = h_sym - h_sym2;
+  std::cout << GridLogMessage << "Projections structure check vector difference (Two Index Symmetric): " << norm2(h_diff_sym) << std::endl;
+
+  
+  // Exponentiate
+  typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix U2iS(grid);
+  U2iS  = expMat(Ar_sym, 1.0, 16);
+  
+  typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix uno2iS(grid);
+  uno2iS = 1.0;
+  // Check matrix U2iS, must be real orthogonal
+  typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix Ucheck2iS = U2iS - conjugate(U2iS);
+  std::cout << GridLogMessage << "Reality check: " << norm2(Ucheck2iS)
+      << std::endl;
+  
+  Ucheck2iS = U2iS * adj(U2iS) - uno2iS;
+  std::cout << GridLogMessage << "orthogonality check 1: " << norm2(Ucheck2iS)
+      << std::endl;
+  Ucheck2iS = adj(U2iS) * U2iS - uno2iS;
+  std::cout << GridLogMessage << "orthogonality check 2: " << norm2(Ucheck2iS)
+      << std::endl;
+  
+  
+  
+  // Construct the fundamental matrix in the group
+  SU<Nc>::LatticeMatrix Af_sym(grid);
+  SU<Nc>::FundamentalLieAlgebraMatrix(h_sym,Af_sym);
+  SU<Nc>::LatticeMatrix Ufund2(grid);
+  Ufund2  = expMat(Af_sym, 1.0, 16);
+  SU<Nc>::LatticeMatrix UnitCheck2(grid);
+  UnitCheck2 = Ufund2 * adj(Ufund2) - uno_f;
+  std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck2)
+      << std::endl;
+  UnitCheck2 = adj(Ufund2) * Ufund2 - uno_f;
+  std::cout << GridLogMessage << "unitarity check 2: " << norm2(UnitCheck2)
+      << std::endl;
+  
+
+  // Tranform to the 2Index Sym representation
+  U = zero; // fill this with only one direction
+  pokeLorentz(U,Ufund2,0); // the representation transf acts on full gauge fields
+  
+  TIndexRep.update_representation(U);
+  Ur2 = TIndexRep.U;  // U_r  
+  typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix Ur02 = peekLorentz(Ur2,0); // this should be the same as U2iS
+  
+  typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix Diff_check_mat2 = Ur02 - U2iS;
+  std::cout << GridLogMessage << "Projections structure check group difference (Two Index Symmetric): " << norm2(Diff_check_mat2) << std::endl;
+  
+
+
+
+  if (SU2TwoIndexAntiSymm::Dimension != 1){
+
+  std::cout << GridLogMessage << "*********************************************"
+      << std::endl;
+    
+  
+  std::cout << GridLogMessage << "Two Index anti-Symmetric: Check Group Structure"
+      << std::endl;
+  // Testing HMC representation classes
+  TwoIndexRep< Nc, AntiSymmetric > TIndexRepA(grid);
+
+
+  // Test group structure
+  // (U_f * V_f)_r = U_r * V_r
+  LatticeGaugeField U2A(grid), V2A(grid);
+  SU<Nc>::HotConfiguration<LatticeGaugeField>(gridRNG, U2A);
+  SU<Nc>::HotConfiguration<LatticeGaugeField>(gridRNG, V2A);
+  
+  LatticeGaugeField UV2A(grid);
+  UV2A = zero;
+  for (int mu = 0; mu < Nd; mu++) {
+    SU<Nc>::LatticeMatrix Umu2A = peekLorentz(U2,mu);
+    SU<Nc>::LatticeMatrix Vmu2A = peekLorentz(V2,mu);
+    pokeLorentz(UV2A,Umu2A*Vmu2A, mu);
+  }
+  
+  TIndexRep.update_representation(UV2A);
+  typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField UVr2A = TIndexRepA.U;  // (U_f * V_f)_r
+  
+  TIndexRep.update_representation(U2A);
+  typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Ur2A = TIndexRepA.U;  // U_r
+  
+  TIndexRep.update_representation(V2A);
+  typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Vr2A = TIndexRepA.U;  // V_r
+  
+  typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Ur2Vr2A(grid);
+  Ur2Vr2A = zero;
+  for (int mu = 0; mu < Nd; mu++) {
+    typename TwoIndexRep< Nc, AntiSymmetric >::LatticeMatrix Urmu2A = peekLorentz(Ur2A,mu);
+    typename TwoIndexRep< Nc, AntiSymmetric >::LatticeMatrix Vrmu2A = peekLorentz(Vr2A,mu);
+    pokeLorentz(Ur2Vr2A,Urmu2A*Vrmu2A, mu);
+  }
+  
+  typename TwoIndexRep< Nc, AntiSymmetric >::LatticeField Diff_check2A = UVr2A - Ur2Vr2A;
+  std::cout << GridLogMessage << "Group structure SU("<<Nc<<") check difference (Two Index anti-Symmetric): " << norm2(Diff_check2A) << std::endl;
+
+  
+  // Check correspondence of algebra and group transformations
+  // Create a random vector
+  SU<Nc>::LatticeAlgebraVector h_Asym(grid);
+  typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Ar_Asym(grid);
+  random(gridRNG,h_Asym);
+  h_Asym = real(h_Asym);
+  SU_TwoIndex< Nc, AntiSymmetric>::TwoIndexLieAlgebraMatrix(h_Asym,Ar_Asym);
+  
+  // Re-extract h_sym
+  SU<Nc>::LatticeAlgebraVector h_Asym2(grid);
+  SU_TwoIndex< Nc, AntiSymmetric>::projectOnAlgebra(h_Asym2, Ar_Asym);
+  SU<Nc>::LatticeAlgebraVector h_diff_Asym = h_Asym - h_Asym2;
+  std::cout << GridLogMessage << "Projections structure check vector difference (Two Index anti-Symmetric): " << norm2(h_diff_Asym) << std::endl;
+
+  
+  // Exponentiate
+  typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix U2iAS(grid);
+  U2iAS  = expMat(Ar_Asym, 1.0, 16);
+  
+  typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix uno2iAS(grid);
+  uno2iAS = 1.0;
+  // Check matrix U2iS, must be real orthogonal
+  typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Ucheck2iAS = U2iAS - conjugate(U2iAS);
+  std::cout << GridLogMessage << "Reality check: " << norm2(Ucheck2iAS)
+      << std::endl;
+  
+  Ucheck2iAS = U2iAS * adj(U2iAS) - uno2iAS;
+  std::cout << GridLogMessage << "orthogonality check 1: " << norm2(Ucheck2iAS)
+      << std::endl;
+  Ucheck2iAS = adj(U2iAS) * U2iAS - uno2iAS;
+  std::cout << GridLogMessage << "orthogonality check 2: " << norm2(Ucheck2iAS)
+      << std::endl;
+  
+  
+  
+  // Construct the fundamental matrix in the group
+  SU<Nc>::LatticeMatrix Af_Asym(grid);
+  SU<Nc>::FundamentalLieAlgebraMatrix(h_Asym,Af_Asym);
+  SU<Nc>::LatticeMatrix Ufund2A(grid);
+  Ufund2A  = expMat(Af_Asym, 1.0, 16);
+  SU<Nc>::LatticeMatrix UnitCheck2A(grid);
+  UnitCheck2A = Ufund2A * adj(Ufund2A) - uno_f;
+  std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck2A)
+      << std::endl;
+  UnitCheck2A = adj(Ufund2A) * Ufund2A - uno_f;
+  std::cout << GridLogMessage << "unitarity check 2: " << norm2(UnitCheck2A)
+      << std::endl;
+  
+
+  // Tranform to the 2Index Sym representation
+  U = zero; // fill this with only one direction
+  pokeLorentz(U,Ufund2A,0); // the representation transf acts on full gauge fields
+  
+  TIndexRepA.update_representation(U);
+  Ur2A = TIndexRepA.U;  // U_r  
+  typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Ur02A = peekLorentz(Ur2A,0); // this should be the same as U2iS
+  
+  typename TwoIndexRep< Nc, AntiSymmetric>::LatticeMatrix Diff_check_mat2A = Ur02A - U2iAS;
+  std::cout << GridLogMessage << "Projections structure check group difference (Two Index anti-Symmetric): " << norm2(Diff_check_mat2A) << std::endl;
+  
+} else  {
+  std::cout << GridLogMessage << "Skipping Two Index anti-Symmetric tests "
+                                 "because representation is trivial (dim = 1)"
+            << std::endl;
+}
+
+
+
+
+
+
+
+
+
+  
   Grid_finalize();
 }

tests/hmc/Make.inc

@@ -1,5 +1,5 @@
-tests: Test_hmc_EODWFRatio Test_hmc_EODWFRatio_Gparity Test_hmc_EOWilsonFermionGauge Test_hmc_EOWilsonRatio Test_hmc_GparityIwasakiGauge Test_hmc_GparityWilsonGauge Test_hmc_IwasakiGauge Test_hmc_RectGauge Test_hmc_WilsonAdjointFermionGauge Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_hmc_WilsonRatio Test_multishift_sqrt Test_remez Test_rhmc_EOWilson1p1 Test_rhmc_EOWilsonRatio Test_rhmc_Wilson1p1 Test_rhmc_WilsonRatio
+tests: Test_hmc_EODWFRatio Test_hmc_EODWFRatio_Gparity Test_hmc_EOWilsonFermionGauge Test_hmc_EOWilsonRatio Test_hmc_GparityIwasakiGauge Test_hmc_GparityWilsonGauge Test_hmc_IwasakiGauge Test_hmc_RectGauge Test_hmc_WilsonAdjointFermionGauge Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_hmc_WilsonMixedRepresentationsFermionGauge Test_hmc_WilsonRatio Test_hmc_WilsonTwoIndexSymmetricFermionGauge Test_multishift_sqrt Test_remez Test_rhmc_EOWilson1p1 Test_rhmc_EOWilsonRatio Test_rhmc_Wilson1p1 Test_rhmc_WilsonRatio
-EXTRA_PROGRAMS = Test_hmc_EODWFRatio Test_hmc_EODWFRatio_Gparity Test_hmc_EOWilsonFermionGauge Test_hmc_EOWilsonRatio Test_hmc_GparityIwasakiGauge Test_hmc_GparityWilsonGauge Test_hmc_IwasakiGauge Test_hmc_RectGauge Test_hmc_WilsonAdjointFermionGauge Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_hmc_WilsonRatio Test_multishift_sqrt Test_remez Test_rhmc_EOWilson1p1 Test_rhmc_EOWilsonRatio Test_rhmc_Wilson1p1 Test_rhmc_WilsonRatio
+EXTRA_PROGRAMS = Test_hmc_EODWFRatio Test_hmc_EODWFRatio_Gparity Test_hmc_EOWilsonFermionGauge Test_hmc_EOWilsonRatio Test_hmc_GparityIwasakiGauge Test_hmc_GparityWilsonGauge Test_hmc_IwasakiGauge Test_hmc_RectGauge Test_hmc_WilsonAdjointFermionGauge Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_hmc_WilsonMixedRepresentationsFermionGauge Test_hmc_WilsonRatio Test_hmc_WilsonTwoIndexSymmetricFermionGauge Test_multishift_sqrt Test_remez Test_rhmc_EOWilson1p1 Test_rhmc_EOWilsonRatio Test_rhmc_Wilson1p1 Test_rhmc_WilsonRatio
 
 Test_hmc_EODWFRatio_SOURCES=Test_hmc_EODWFRatio.cc
 Test_hmc_EODWFRatio_LDADD=-lGrid
@@ -34,9 +34,15 @@ Test_hmc_WilsonFermionGauge_LDADD=-lGrid
 Test_hmc_WilsonGauge_SOURCES=Test_hmc_WilsonGauge.cc
 Test_hmc_WilsonGauge_LDADD=-lGrid
 
+Test_hmc_WilsonMixedRepresentationsFermionGauge_SOURCES=Test_hmc_WilsonMixedRepresentationsFermionGauge.cc
+Test_hmc_WilsonMixedRepresentationsFermionGauge_LDADD=-lGrid
+
 Test_hmc_WilsonRatio_SOURCES=Test_hmc_WilsonRatio.cc
 Test_hmc_WilsonRatio_LDADD=-lGrid
 
+Test_hmc_WilsonTwoIndexSymmetricFermionGauge_SOURCES=Test_hmc_WilsonTwoIndexSymmetricFermionGauge.cc
+Test_hmc_WilsonTwoIndexSymmetricFermionGauge_LDADD=-lGrid
+
 Test_multishift_sqrt_SOURCES=Test_multishift_sqrt.cc
 Test_multishift_sqrt_LDADD=-lGrid
 

tests/hmc/Test_hmc_WilsonMixedRepresentationsFermionGauge.cc

@@ -0,0 +1,113 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./tests/Test_hmc_WilsonAdjointFermionGauge.cc
+
+Copyright (C) 2015
+
+Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: neo <cossu@post.kek.jp>
+Author: paboyle <paboyle@ph.ed.ac.uk>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
+#include "Grid/Grid.h"
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+namespace Grid {
+namespace QCD {
+
+// Here change the allowed (higher) representations
+typedef Representations< FundamentalRepresentation, AdjointRepresentation , TwoIndexSymmetricRepresentation> TheRepresentations;
+
+class HmcRunner : public NerscHmcRunnerHirep< TheRepresentations > {
+ public:
+  void BuildTheAction(int argc, char **argv)
+
+  {
+    typedef WilsonAdjImplR AdjImplPolicy; // gauge field implemetation for the pseudofermions
+    typedef WilsonAdjFermionR AdjFermionAction; // type of lattice fermions (Wilson, DW, ...)
+    typedef WilsonTwoIndexSymmetricImplR SymmImplPolicy; 
+    typedef WilsonTwoIndexSymmetricFermionR SymmFermionAction; 
+
+ 
+    typedef typename AdjFermionAction::FermionField AdjFermionField;
+    typedef typename SymmFermionAction::FermionField SymmFermionField;
+
+    UGrid = SpaceTimeGrid::makeFourDimGrid(
+        GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()),
+        GridDefaultMpi());
+    UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+
+    FGrid = UGrid;
+    FrbGrid = UrbGrid;
+
+    // temporarily need a gauge field
+    //LatticeGaugeField U(UGrid);
+    AdjointRepresentation::LatticeField UA(UGrid);
+    TwoIndexSymmetricRepresentation::LatticeField US(UGrid);
+
+    // Gauge action
+    WilsonGaugeActionR Waction(2.25);
+
+    Real adjoint_mass = -0.1;
+    Real symm_mass = -0.5;
+    AdjFermionAction AdjFermOp(UA, *FGrid, *FrbGrid, adjoint_mass);
+    SymmFermionAction SymmFermOp(US, *FGrid, *FrbGrid, symm_mass);
+
+    ConjugateGradient<AdjFermionField> CG_adj(1.0e-8, 10000, false);
+    ConjugateGradient<SymmFermionField> CG_symm(1.0e-8, 10000, false);
+
+    // Pass two solvers: one for the force computation and one for the action
+    TwoFlavourPseudoFermionAction<AdjImplPolicy> Nf2_Adj(AdjFermOp, CG_adj, CG_adj);
+    TwoFlavourPseudoFermionAction<SymmImplPolicy> Nf2_Symm(SymmFermOp, CG_symm, CG_symm);
+
+    // Collect actions
+    ActionLevelHirep<LatticeGaugeField, TheRepresentations > Level1(1);
+    Level1.push_back(&Nf2_Adj);
+    Level1.push_back(&Nf2_Symm);
+
+    ActionLevelHirep<LatticeGaugeField, TheRepresentations > Level2(4);
+    Level2.push_back(&Waction);
+
+    TheAction.push_back(Level1);
+    TheAction.push_back(Level2);
+
+    Run(argc, argv);
+  };
+};
+}
+}
+
+int main(int argc, char **argv) {
+  Grid_init(&argc, &argv);
+
+  int threads = GridThread::GetThreads();
+  std::cout << GridLogMessage << "Grid is setup to use " << threads
+            << " threads" << std::endl;
+
+  HmcRunner TheHMC;
+
+  TheHMC.BuildTheAction(argc, argv);
+}

tests/hmc/Test_hmc_WilsonTwoIndexSymmetricFermionGauge.cc

@@ -0,0 +1,103 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./tests/Test_hmc_WilsonAdjointFermionGauge.cc
+
+Copyright (C) 2015
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: neo <cossu@post.kek.jp>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution
+directory
+*************************************************************************************/
+/*  END LEGAL */
+#include "Grid/Grid.h"
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+namespace Grid {
+namespace QCD {
+
+// Here change the allowed (higher) representations
+typedef Representations< FundamentalRepresentation, TwoIndexSymmetricRepresentation > TheRepresentations;
+
+
+class HmcRunner : public NerscHmcRunnerHirep< TheRepresentations > {
+ public:
+  void BuildTheAction(int argc, char **argv)
+
+  {
+    typedef WilsonTwoIndexSymmetricImplR ImplPolicy; // gauge field implemetation for the pseudofermions
+    typedef WilsonTwoIndexSymmetricFermionR FermionAction; // type of lattice fermions (Wilson, DW, ...)
+    typedef typename FermionAction::FermionField FermionField;
+
+    UGrid = SpaceTimeGrid::makeFourDimGrid(
+        GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()),
+        GridDefaultMpi());
+    UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+
+    FGrid = UGrid;
+    FrbGrid = UrbGrid;
+
+    // temporarily need a gauge field
+    TwoIndexSymmetricRepresentation::LatticeField U(UGrid);
+
+    // Gauge action
+    WilsonGaugeActionR Waction(2.0);
+
+    Real mass = -0.0;
+    FermionAction FermOp(U, *FGrid, *FrbGrid, mass);
+
+    ConjugateGradient<FermionField> CG(1.0e-8, 10000, false);
+
+    // Pass two solvers: one for the force computation and one for the action
+    TwoFlavourPseudoFermionAction<ImplPolicy> Nf2(FermOp, CG, CG);
+
+    // Set smearing (true/false), default: false
+    Nf2.is_smeared = false;
+
+    // Collect actions
+    ActionLevelHirep<LatticeGaugeField, TheRepresentations > Level1(1);
+    Level1.push_back(&Nf2);
+
+    ActionLevelHirep<LatticeGaugeField, TheRepresentations > Level2(4);
+    Level2.push_back(&Waction);
+
+    TheAction.push_back(Level1);
+    TheAction.push_back(Level2);
+
+    Run(argc, argv);
+  };
+};
+}
+}
+
+int main(int argc, char **argv) {
+  Grid_init(&argc, &argv);
+
+  int threads = GridThread::GetThreads();
+  std::cout << GridLogMessage << "Grid is setup to use " << threads
+            << " threads" << std::endl;
+
+  HmcRunner TheHMC;
+
+  TheHMC.BuildTheAction(argc, argv);
+}