gauge and fermion implementation for sp2n

Grid/qcd/action/fermion/Fermion.h

@@ -115,6 +115,11 @@ typedef WilsonFermion<WilsonImplR> WilsonFermionR;
 typedef WilsonFermion<WilsonImplF> WilsonFermionF;
 typedef WilsonFermion<WilsonImplD> WilsonFermionD;
 
+
+typedef WilsonFermion<SpWilsonImplR> SpWilsonFermionR;
+typedef WilsonFermion<SpWilsonImplF> SpWilsonFermionF;
+typedef WilsonFermion<SpWilsonImplD> SpWilsonFermionD;
+
 //typedef WilsonFermion<WilsonImplRL> WilsonFermionRL;
 //typedef WilsonFermion<WilsonImplFH> WilsonFermionFH;
 //typedef WilsonFermion<WilsonImplDF> WilsonFermionDF;

Grid/qcd/action/fermion/WilsonFermion.h

@@ -208,5 +208,7 @@ public:
 
 typedef WilsonFermion<WilsonImplF> WilsonFermionF;
 typedef WilsonFermion<WilsonImplD> WilsonFermionD;
+//typedef WilsonFermion<SpWilsonImplF> SpWilsonFermionF;
+//typedef WilsonFermion<SpWilsonImplD> SpWilsonFermionD;
 
 NAMESPACE_END(Grid);

Grid/qcd/action/fermion/WilsonImpl.h

@@ -267,6 +267,16 @@ typedef WilsonImpl<vComplex,  TwoIndexAntiSymmetricRepresentation, CoeffReal > W
 typedef WilsonImpl<vComplexF, TwoIndexAntiSymmetricRepresentation, CoeffReal > WilsonTwoIndexAntiSymmetricImplF;  // Float
 typedef WilsonImpl<vComplexD, TwoIndexAntiSymmetricRepresentation, CoeffReal > WilsonTwoIndexAntiSymmetricImplD;  // Double
 
+//sp 2n
+
+typedef WilsonImpl<vComplex,  SpFundamentalRepresentation, CoeffReal > SpWilsonImplR;  // Real.. whichever prec
+typedef WilsonImpl<vComplexF, SpFundamentalRepresentation, CoeffReal > SpWilsonImplF;  // Float
+typedef WilsonImpl<vComplexD, SpFundamentalRepresentation, CoeffReal > SpWilsonImplD;  // Double
+
+//typedef WilsonImpl<vComplex,  SpFundamentalRepresentation, CoeffComplex > SpZWilsonImplR; // Real.. whichever prec
+//typedef WilsonImpl<vComplexF, SpFundamentalRepresentation, CoeffComplex > SpZWilsonImplF; // Float
+//typedef WilsonImpl<vComplexD, SpFundamentalRepresentation, CoeffComplex > SpZWilsonImplD; // Double
+
 
 NAMESPACE_END(Grid);
 

Grid/qcd/action/gauge/Gauge.h

@@ -39,6 +39,9 @@ NAMESPACE_BEGIN(Grid);
 typedef WilsonGaugeAction<PeriodicGimplR>          WilsonGaugeActionR;
 typedef WilsonGaugeAction<PeriodicGimplF>          WilsonGaugeActionF;
 typedef WilsonGaugeAction<PeriodicGimplD>          WilsonGaugeActionD;
+typedef WilsonGaugeAction<SymplPeriodicGimplR>     SymplWilsonGaugeActionR;
+typedef WilsonGaugeAction<SymplPeriodicGimplF>     SymplWilsonGaugeActionF;
+typedef WilsonGaugeAction<SymplPeriodicGimplD>     SymplWilsonGaugeActionD;
 typedef PlaqPlusRectangleAction<PeriodicGimplR>    PlaqPlusRectangleActionR;
 typedef PlaqPlusRectangleAction<PeriodicGimplF>    PlaqPlusRectangleActionF;
 typedef PlaqPlusRectangleAction<PeriodicGimplD>    PlaqPlusRectangleActionD;

Grid/qcd/action/gauge/GaugeImplTypes.h

@@ -61,7 +61,8 @@ NAMESPACE_BEGIN(Grid);
   typedef typename Impl::Field Field;
 
 // hardcodes the exponential approximation in the template
-template <class S, int Nrepresentation = Nc, int Nexp = 12 > class GaugeImplTypes {
+//template <class S, int Nrepresentation = Nc, int Nexp = 12 > class GaugeImplTypes {
+template <class S, int Nrepresentation = Nc, int Nexp = 12, bool isSp2n = false > class GaugeImplTypes {
 public:
   typedef S Simd;
   typedef typename Simd::scalar_type scalar_type;
@@ -70,6 +71,7 @@ public:
   template <typename vtype> using iImplGaugeLink  = iScalar<iScalar<iMatrix<vtype, Nrepresentation> > >;
   template <typename vtype> using iImplGaugeField = iVector<iScalar<iMatrix<vtype, Nrepresentation> >, Nd>;
 
+
   typedef iImplScalar<Simd>     SiteComplex;
   typedef iImplGaugeLink<Simd>  SiteLink;
   typedef iImplGaugeField<Simd> SiteField;
@@ -117,8 +119,17 @@ public:
     //
     LinkField Pmu(P.Grid());
     Pmu = Zero();
-    for (int mu = 0; mu < Nd; mu++) {
+    for (int mu = 0; mu < Nd; mu++)
-      SU<Nrepresentation>::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu);
+    {
+        if (isSp2n == true)
+        {
+            const int nSp = Nrepresentation/2;
+            Sp<nSp>::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu);
+        } else
+        {
+        
+            SU<Nrepresentation>::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu);
+        }
       RealD scale = ::sqrt(HMC_MOMENTUM_DENOMINATOR) ;
       Pmu = Pmu*scale;
       PokeIndex<LorentzIndex>(P, Pmu, mu);
@@ -135,14 +146,21 @@ public:
     autoView(P_v,P,AcceleratorRead);
     accelerator_for(ss, P.Grid()->oSites(),1,{
       for (int mu = 0; mu < Nd; mu++) {
-        U_v[ss](mu) = ProjectOnGroup(Exponentiate(P_v[ss](mu), ep, Nexp) * U_v[ss](mu));
+          if (isSp2n == true)
+          {
+              U_v[ss](mu) = ProjectOnSpGroup(Exponentiate(P_v[ss](mu), ep, Nexp) * U_v[ss](mu));
+          } else
+          {
+              U_v[ss](mu) = ProjectOnGroup(Exponentiate(P_v[ss](mu), ep, Nexp) * U_v[ss](mu));
+          }
+        
       }
     });
    //auto end = std::chrono::high_resolution_clock::now();
    // diff += end - start;
    // std::cout << "Time to exponentiate matrix " << diff.count() << " s\n";
   }
-
+    
   static inline RealD FieldSquareNorm(Field& U){
     LatticeComplex Hloc(U.Grid());
     Hloc = Zero();
@@ -154,21 +172,77 @@ public:
     return Hsum.real();
   }
 
-  static inline void Project(Field &U) {
+  static inline void Project(Field &U)
-    ProjectSUn(U);
+  {
+      if (isSp2n == true)
+      {
+          ProjectSp2n(U);
+      } else
+      {
+          ProjectSUn(U);
+      }
   }
 
-  static inline void HotConfiguration(GridParallelRNG &pRNG, Field &U) {
+  static inline void HotConfiguration(GridParallelRNG &pRNG, Field &U)
-    SU<Nc>::HotConfiguration(pRNG, U);
+  {
+      SU<Nc>::HotConfiguration(pRNG, U);
   }
 
   static inline void TepidConfiguration(GridParallelRNG &pRNG, Field &U) {
     SU<Nc>::TepidConfiguration(pRNG, U);
   }
 
-  static inline void ColdConfiguration(GridParallelRNG &pRNG, Field &U) {
+  static inline void ColdConfiguration(GridParallelRNG &pRNG, Field &U)
-    SU<Nc>::ColdConfiguration(pRNG, U);
+  {
+      if (isSp2n == true)
+      {
+          const int nSp = Nrepresentation/2;
+          Sp<nSp>::ColdConfiguration(pRNG, U);
+      } else
+      {
+          SU<Nc>::ColdConfiguration(pRNG, U);
+      }
   }
+    
+  //sp2n... see sp2n.h
+/*
+    static inline void generate_sp2n_momenta(Field &P, GridSerialRNG & sRNG, GridParallelRNG &pRNG)
+    {
+      
+      const int nSp = Nrepresentation/2;
+      LinkField Pmu(P.Grid());
+      Pmu = Zero();
+      for (int mu = 0; mu < Nd; mu++) {
+        
+        Sp<nSp>::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu);
+        RealD scale = ::sqrt(HMC_MOMENTUM_DENOMINATOR) ;
+        Pmu = Pmu*scale;
+        PokeIndex<LorentzIndex>(P, Pmu, mu);
+      }
+    }
+    
+    static inline void update_sp2n_field(Field& P, Field& U, double ep){
+ 
+      autoView(U_v,U,AcceleratorWrite);
+      autoView(P_v,P,AcceleratorRead);
+      accelerator_for(ss, P.Grid()->oSites(),1,{
+        for (int mu = 0; mu < Nd; mu++) {
+          U_v[ss](mu) = ProjectOnSpGroup(Exponentiate(P_v[ss](mu), ep, Nexp) * U_v[ss](mu));
+        }
+      });
+    }
+    
+    static inline void SpProject(Field &U) {
+      ProjectSp2n(U);
+    }
+
+
+    static inline void ColdSpConfiguration(GridParallelRNG &pRNG, Field &U) {
+      const int nSp = Nrepresentation/2;
+      Sp<nSp>::ColdConfiguration(pRNG, U);
+    }
+*/
+
 };
 
 
@@ -176,10 +250,16 @@ typedef GaugeImplTypes<vComplex, Nc> GimplTypesR;
 typedef GaugeImplTypes<vComplexF, Nc> GimplTypesF;
 typedef GaugeImplTypes<vComplexD, Nc> GimplTypesD;
 
+typedef GaugeImplTypes<vComplex, Nc, 12, true> SymplGimplTypesR;
+typedef GaugeImplTypes<vComplexF, Nc, 12, true> SymplGimplTypesF;
+typedef GaugeImplTypes<vComplexD, Nc, 12, true> SymplGimplTypesD;
+
 typedef GaugeImplTypes<vComplex, SU<Nc>::AdjointDimension> GimplAdjointTypesR;
 typedef GaugeImplTypes<vComplexF, SU<Nc>::AdjointDimension> GimplAdjointTypesF;
 typedef GaugeImplTypes<vComplexD, SU<Nc>::AdjointDimension> GimplAdjointTypesD;
 
+
+
 NAMESPACE_END(Grid);
 
 #endif // GRID_GAUGE_IMPL_TYPES_H

Grid/qcd/action/gauge/GaugeImplementations.h

@@ -155,6 +155,11 @@ typedef ConjugateGaugeImpl<GimplTypesR> ConjugateGimplR; // Real.. whichever pre
 typedef ConjugateGaugeImpl<GimplTypesF> ConjugateGimplF; // Float
 typedef ConjugateGaugeImpl<GimplTypesD> ConjugateGimplD; // Double
 
+typedef PeriodicGaugeImpl<SymplGimplTypesR> SymplPeriodicGimplR; // Real.. whichever prec
+typedef PeriodicGaugeImpl<SymplGimplTypesF> SymplPeriodicGimplF; // Float
+typedef PeriodicGaugeImpl<SymplGimplTypesD> SymplPeriodicGimplD; // Double
+
+
 NAMESPACE_END(Grid);
 
 #endif

Grid/qcd/representations/Representations.h

@@ -4,6 +4,7 @@
 #include <Grid/qcd/representations/adjoint.h>
 #include <Grid/qcd/representations/two_index.h>
 #include <Grid/qcd/representations/fundamental.h>
+#include <Grid/qcd/representations/spfundamental.h>
 #include <Grid/qcd/representations/hmc_types.h>
 
 #endif

Grid/qcd/representations/spfundamental.h

@@ -0,0 +1,42 @@
+
+
+#ifndef SPFUNDAMENTAL_H
+#define SPFUNDAMENTAL_H
+
+NAMESPACE_BEGIN(Grid);
+
+/*
+ * This is an helper class for the HMC
+ * Empty since HMC updates already the fundamental representation 
+ */
+
+template <int ncolour>
+class SpFundamentalRep {
+public:
+  static const int Dimension = ncolour;
+  static const int nSp = ncolour/2;
+  static const bool isFundamental = true;
+
+  // typdef to be used by the Representations class in HMC to get the
+  // types for the higher representation fields
+  typedef typename Sp<nSp>::LatticeMatrix LatticeMatrix;
+  typedef LatticeGaugeField LatticeField;
+  
+  explicit SpFundamentalRep(GridBase* grid) {} //do nothing
+  void update_representation(const LatticeGaugeField& Uin) {} // do nothing
+
+  LatticeField RtoFundamentalProject(const LatticeField& in, Real scale = 1.0) const{
+    return (scale * in);
+  }
+
+};
+
+
+    
+
+  
+typedef	 SpFundamentalRep<Nc> SpFundamentalRepresentation;
+
+NAMESPACE_END(Grid);  
+
+#endif

Grid/qcd/utils/Sp2n.h

@@ -6,6 +6,13 @@ NAMESPACE_BEGIN(Grid);
 
 // Sp(2N)
 // ncolour = N
+
+// need to be careful with n and 2n
+// I am defining the Sp class for Sp(2n) to be such that the template variable ncolour
+// is n inside 2n and the typedef at the end of the file should eliminate possible confusion.
+
+// the other routines, like projectOnSp2n, N will be the dimension of the actual number of colors for consistency with the sun routines
+
 template <int ncolour>
 class Sp {
 public:
@@ -22,11 +29,31 @@ public:
     using iSp2nAlgebraVector = iScalar<iScalar<iVector<vtype, AlgebraDimension> > >;
     
     typedef iSp2nMatrix<Complex> Matrix;
+    typedef iSp2nMatrix<ComplexF> MatrixF;
+    typedef iSp2nMatrix<ComplexD> MatrixD;
     
     typedef iSp2nMatrix<vComplex> vMatrix;
+    typedef iSp2nMatrix<vComplexF> vMatrixF;
+    typedef iSp2nMatrix<vComplexD> vMatrixD;
     
+    typedef iSp2nAlgebraVector<Complex> AlgebraVector;
+    typedef iSp2nAlgebraVector<ComplexF> AlgebraVectorF;
+    typedef iSp2nAlgebraVector<ComplexD> AlgebraVectorD;
+    
+    typedef iSp2nAlgebraVector<vComplex> vAlgebraVector;
+    typedef iSp2nAlgebraVector<vComplexF> vAlgebraVectorF;
+    typedef iSp2nAlgebraVector<vComplexD> vAlgebraVectorD;
     
     typedef Lattice<vMatrix> LatticeMatrix;
+    typedef Lattice<vMatrixF> LatticeMatrixF;
+    typedef Lattice<vMatrixD> LatticeMatrixD;
+    
+    typedef Lattice<vAlgebraVector> LatticeAlgebraVector;
+    typedef Lattice<vAlgebraVectorF> LatticeAlgebraVectorF;
+    typedef Lattice<vAlgebraVectorD> LatticeAlgebraVectorD;
+    
+    
+    
     // Sp(2N) has N(2N+1) = 2N^2+N generators
     //
     // normalise the generators such that
@@ -322,75 +349,113 @@ public:
       }
     }
     
+ 
     
     
+/*    template <typename GaugeField>
+    static void HotSpConfiguration(GridParallelRNG &pRNG, GaugeField &out) {
+      typedef typename GaugeField::vector_type vector_type;
+      typedef iSp2nMatrix<vector_type> vMatrixType;
+      typedef Lattice<vMatrixType> LatticeMatrixType;
+
+      LatticeMatrixType Umu(out.Grid());
+      for (int mu = 0; mu < Nd; mu++) {
+        LieRandomize(pRNG, Umu, 1.0);   //def
+        PokeIndex<LorentzIndex>(out, Umu, mu);
+      }
+    }
+    template<typename GaugeField>
+    static void TepidSpConfiguration(GridParallelRNG &pRNG,GaugeField &out){
+      typedef typename GaugeField::vector_type vector_type;
+      typedef iSp2nMatrix<vector_type> vMatrixType;
+      typedef Lattice<vMatrixType> LatticeMatrixType;
+
+      LatticeMatrixType Umu(out.Grid());
+      for(int mu=0;mu<Nd;mu++){
+        LieRandomize(pRNG,Umu,0.01);    //def
+        PokeIndex<LorentzIndex>(out,Umu,mu);
+      }
+    } */
+    
+    template<typename GaugeField>
+    static void ColdConfiguration(GaugeField &out){
+      typedef typename GaugeField::vector_type vector_type;
+      typedef iSp2nMatrix<vector_type> vMatrixType;
+      typedef Lattice<vMatrixType> LatticeMatrixType;
+
+      LatticeMatrixType Umu(out.Grid());
+      Umu=1.0;
+      for(int mu=0;mu<Nd;mu++){
+        PokeIndex<LorentzIndex>(out,Umu,mu);
+      }
+    }
+    template<typename GaugeField>
+    static void ColdConfiguration(GridParallelRNG &pRNG,GaugeField &out){
+      ColdConfiguration(out);
+    }
+    
     
     
 }; // end of class Sp
     
     
 
+
+template<int N> //same of su
+ LatticeComplexD SpDeterminant(const Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu)
+ {
+   GridBase *grid=Umu.Grid();
+   auto lvol = grid->lSites();
+   LatticeComplexD ret(grid);
+
+   autoView(Umu_v,Umu,CpuRead);
+   autoView(ret_v,ret,CpuWrite);
+   thread_for(site,lvol,{
+     Eigen::MatrixXcd EigenU = Eigen::MatrixXcd::Zero(N,N);
+     Coordinate lcoor;
+     grid->LocalIndexToLocalCoor(site, lcoor);
+     iScalar<iScalar<iMatrix<ComplexD, N> > > Us;
+     peekLocalSite(Us, Umu_v, lcoor);
+     for(int i=0;i<N;i++){
+       for(int j=0;j<N;j++){
+     EigenU(i,j) = Us()()(i,j);
+       }}
+     ComplexD det = EigenU.determinant();
+     pokeLocalSite(det,ret_v,lcoor);
+   });
+   return ret;
+ }
+ 
+ template<int N>
+ static void ProjectSp2n(Lattice<iScalar<iScalar<iMatrix<vComplexD, N> > > > &Umu)
+ {
+   Umu      = ProjectOnSpGroup(Umu);
+   auto det = SpDeterminant(Umu); // ok ?
+
+   det = conjugate(det);
+
+   for(int i=0;i<N;i++){
+     auto element = PeekIndex<ColourIndex>(Umu,N-1,i);
+     element = element * det;
+     PokeIndex<ColourIndex>(Umu,element,Nc-1,i);
+   }
+ }
+ template<int N>
+ static void ProjectSp2n(Lattice<iVector<iScalar<iMatrix<vComplexD, N> >,Nd> > &U)
+ {
+   GridBase *grid=U.Grid();
+   for(int mu=0;mu<Nd;mu++){
+     auto Umu = PeekIndex<LorentzIndex>(U,mu);
+     Umu      = ProjectOnSpGroup(Umu);
+     ProjectSp2n(Umu);
+     PokeIndex<LorentzIndex>(U,Umu,mu);
+   }
+ }
+
+typedef Sp<1> Sp2;
 typedef Sp<2> Sp4;
 typedef Sp<3> Sp6;
 typedef Sp<4> Sp8;
 
 NAMESPACE_END(Grid);
 #endif
-
-
-    
-    
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-    
-    
-    
-/*     }
-     sigxy = lieIndex & 0x1;     // 1 if odd, 0 if even
-     su2Index = lieIndex >> 1 ;  // where to put the sigma_x(y)
-                                 //for the even(odd) lieindex, sigmax(y)
-     if (sigxy)
-         //put sigmay at su2index
-     else
-         //put sigmax    */

tests/core/test_sp2n_lie_gen.cc

@@ -2,7 +2,7 @@
 
 
 #include <Grid/Grid.h>
-#include <Grid/qcd/utils/SUn.h>
+#include <Grid/qcd/utils/Sp2n.h>
 
 using namespace Grid;
 
@@ -13,6 +13,15 @@ int main(int argc, char** argv) {
   //GridCartesian* grid = SpaceTimeGrid::makeFourDimGrid(
   //latt, GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi());
   //GridRedBlackCartesian* rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(grid);
+    
+    std::cout << GridLogMessage << "*********************************************"
+              << std::endl;
+    std::cout << GridLogMessage << "* Generators for Sp(2)" << std::endl;
+    std::cout << GridLogMessage << "*********************************************"
+              << std::endl;
+      
+    Sp2::printGenerators();
+    Sp2::testGenerators();
   
   std::cout << GridLogMessage << "*********************************************"
             << std::endl;