Merge branch 'sp2n/unify_twoindex' of github.com:LupoA/Grid into sp2n/unify_twoindex

Grid/qcd/representations/sp_two_index.h

@@ -27,7 +27,7 @@ public:
   // types for the higher representation fields
   typedef typename Sp_TwoIndex<ncolour, S>::LatticeTwoIndexMatrix LatticeMatrix;
   typedef typename Sp_TwoIndex<ncolour, S>::LatticeTwoIndexField LatticeField;
-  static const int Dimension = (ncolour * (ncolour + S) / 2) - 1;
+  static const int Dimension = (ncolour * (ncolour + S) / 2) + S;
   static const bool isFundamental = false;
   //static const int nsp = Nc / 2;
   LatticeField U;

Grid/qcd/utils/GaugeGroupTwoIndex.h

@@ -85,7 +85,7 @@ class GaugeGroupTwoIndex : public GaugeGroup<ncolour, group_name> {
   template <class cplx>
   static void base(int Index, iGroupMatrix<cplx> &eij) {
     // returns (e)^(ij)_{kl} necessary for change of base U_F -> U_R
-    assert(Index < NumGenerators);
+    assert(Index < Dimension);
     eij = Zero();
 
     // for the linearisation of the 2 indexes
@@ -96,9 +96,10 @@ class GaugeGroupTwoIndex : public GaugeGroup<ncolour, group_name> {
       for (int i = 1; i < ncolour; i++) {
         for (int j = 0; j < i; j++) {
           a[counter][0] = i;
-          j = i == ngroup ? j + 1
-                          : j;  // this will only ever trigger for Sp because
-                                // ngroup == ncolour is out of range for SU
+            if (j==0 && ngroup == ncolour/2 && i==ngroup+j) {
+                //std::cout << "skipping" << std::endl; // for Sp2n this vanishes identically.
+                j = j+1;
+            }
           a[counter][1] = j;
           counter++;
         }
@@ -328,6 +329,7 @@ typedef SU_TwoIndex<3, Symmetric> SU3TwoIndexSymm;
 typedef SU_TwoIndex<4, Symmetric> SU4TwoIndexSymm;
 typedef SU_TwoIndex<5, Symmetric> SU5TwoIndexSymm;
 
+typedef SU_TwoIndex<2, AntiSymmetric> SU2TwoIndexAntiSymm;
 typedef SU_TwoIndex<3, AntiSymmetric> SU3TwoIndexAntiSymm;
 typedef SU_TwoIndex<4, AntiSymmetric> SU4TwoIndexAntiSymm;
 typedef SU_TwoIndex<5, AntiSymmetric> SU5TwoIndexAntiSymm;

tests/sp2n/Test_2as_base.cc

@@ -0,0 +1,110 @@
+#include <Grid/Grid.h>
+
+#define verbose 0
+
+using namespace Grid;
+
+int main(int argc, char** argv) {
+  Grid_init(&argc, &argv);
+
+  const int this_nc = 4;
+  const int this_n = this_nc/2;
+  const int this_irrep_dim = Sp_TwoIndex<this_nc, AntiSymmetric>::Dimension;
+  const int this_algebra_dim = Sp<this_nc>::AlgebraDimension;
+  typedef Sp_TwoIndex<this_nc, AntiSymmetric>::iGroupMatrix<Complex> Matrix;
+  typedef Sp_TwoIndex<this_nc, AntiSymmetric>::iGroupTwoIndexMatrix<Complex> ASMatrix;
+  
+  Matrix Omega;
+  Matrix eij_a;
+  Matrix eij_b;
+  Matrix eij_c;
+  Matrix e_sum;
+  Omega = Zero();
+  for (int i = 0; i < this_n; i++)
+  {
+    Omega()()(i, this_n + i) = 1.;
+    Omega()()(this_n + i, i) = -1;
+  }
+    std::cout << "Omega " << Omega << std::endl;
+
+  RealD realA;
+  RealD realB;
+  std::cout << GridLogMessage << "2as dimension is " << this_irrep_dim << std::endl;
+  std::cout << GridLogMessage << "algebra dimension is " << this_algebra_dim << std::endl;
+  realA = Sp_TwoIndex<this_nc, AntiSymmetric>::Dimension + Sp_TwoIndex<this_nc, Symmetric>::Dimension;
+  realB = Sp<this_nc>::Dimension*Sp<this_nc>::Dimension;
+  assert ( realA == realB);
+
+  std::cout << GridLogMessage << "checking base is antisymmetric " << std::endl;
+  for (int a=0; a < this_irrep_dim; a++)
+  {
+      Sp_TwoIndex<this_nc, AntiSymmetric>::base(a, eij_c);
+      e_sum = eij_c + transpose(eij_c);
+      std::cout << GridLogMessage << "e_ab + e_ab^T " << norm2(e_sum) << std::endl;
+      assert(norm2(e_sum) < 1e-8);
+        
+  }
+  std::cout << GridLogMessage << "Checking Tr (e^(ab) Omega ) = 0 and Tr (e^(ab) e^(cd) = delta^((ab)(cd)) ) " << std::endl;
+  for (int a=0; a < Sp_TwoIndex<this_nc, AntiSymmetric>::Dimension; a++) {
+      Sp_TwoIndex<this_nc, AntiSymmetric>::base(a, eij_a);
+      realA = norm2(trace(Omega*eij_a));
+      std::cout << GridLogMessage << "Omega trace for (ab) = " << a << std::endl;
+      assert(realA == 0);
+      for (int b=0; b < Sp_TwoIndex<this_nc, AntiSymmetric>::Dimension; b++) {
+          Sp_TwoIndex<this_nc, AntiSymmetric>::base(b, eij_b);
+          auto d_ab = TensorRemove(trace(eij_a * eij_b));
+#if verbose
+          std::cout << GridLogMessage << "Tr( e_{ab=" << a << "} e_{cd=" << b << "} ) = " << d_ab << std::endl;
+#endif
+          std::cout << GridLogMessage << "Orthonormality for (ab) = " << a << std::endl;
+          if (a==b) {
+              assert(real(d_ab)+1 < 1e-8);
+              assert(imag(d_ab) < 1e-8);
+          } else {
+              assert(real(d_ab) < 1e-8);
+              assert(imag(d_ab) < 1e-8);
+          }
+      }
+  }
+
+  int sum = 0;
+  int sum_im = 0;
+  Vector<Matrix> ta_fund(this_algebra_dim);
+  Vector<Matrix> eij(this_irrep_dim);
+  Matrix tmp_l;
+  Matrix tmp_r;
+  for (int n = 0; n < this_algebra_dim; n++)
+  {
+      Sp<this_nc>::generator(n, ta_fund[n]);
+  }
+  for (int a = 0; a < this_irrep_dim; a++)
+  {
+      Sp_TwoIndex<this_nc, AntiSymmetric>::base(a, eij[a]);
+  }
+  for (int gen_id = 0; gen_id < this_algebra_dim; gen_id++)
+  {
+      Complex iTr;
+      sum = 0;
+      sum_im = 0;
+      std::cout << GridLogMessage <<  "generator number " << gen_id << std::endl;
+      for (int a = 0; a < this_irrep_dim; a++)
+      {
+      
+          tmp_l = adj(eij[a])*ta_fund[gen_id]*eij[a];
+          tmp_r = adj(eij[a])*eij[a]*transpose(ta_fund[gen_id]);
+#if verbose
+          std::cout << GridLogMessage << " as_indx = " << a << " eDag T_F e = " << std::endl << tmp_l << std::endl;
+          std::cout << GridLogMessage << " as_indx = " << a << " eDag e T_F^T = " << std::endl << tmp_r << std::endl;
+#endif
+          std::cout << GridLogMessage << " as_indx = " << a << " Tr(sum) = " << TensorRemove(trace(tmp_l+tmp_r)) << std::endl;
+          sum += real(TensorRemove(trace(tmp_l+tmp_r)));
+          sum_im += imag(TensorRemove(trace(tmp_l+tmp_r)));
+      }
+      std::cout << GridLogMessage << "re-evaluated trace of the generator " << gen_id << " is " << sum << " " << sum_im << std::endl;
+      assert ( sum < 1e-8) ;
+      assert ( sum_im < 1e-8) ;
+
+  }
+  
+  Grid_finalize();
+}

tests/sp2n/Test_sp2n_lie_gen.cc

@@ -4,10 +4,10 @@ using namespace Grid;
 
 int main(int argc, char** argv) {
   Grid_init(&argc, &argv);
-    
+
   std::cout << GridLogMessage << "*********************************************"
               << std::endl;
-  std::cout << GridLogMessage << "* Generators for Sp(2)" << std::endl;
+  std::cout << GridLogMessage << "* Generators for Sp(2) (print and test)" << std::endl;
   std::cout << GridLogMessage << "*********************************************"
               << std::endl;
       
@@ -16,7 +16,7 @@ int main(int argc, char** argv) {
   
   std::cout << GridLogMessage << "*********************************************"
             << std::endl;
-  std::cout << GridLogMessage << "* Generators for Sp(4)" << std::endl;
+  std::cout << GridLogMessage << "* Generators for Sp(4) (print and test)" << std::endl;
   std::cout << GridLogMessage << "*********************************************"
             << std::endl;
     
@@ -25,47 +25,42 @@ int main(int argc, char** argv) {
     
   std::cout << GridLogMessage << "*********************************************"
             << std::endl;
-  std::cout << GridLogMessage << "* Generators for Sp(6)" << std::endl;
+  std::cout << GridLogMessage << "* Generators for Sp(6) (test)" << std::endl;
   std::cout << GridLogMessage << "*********************************************"
             << std::endl;
-  
-  Sp6::printGenerators();
+
   Sp6::testGenerators();
-    
+
   std::cout << GridLogMessage << "*********************************************"
               << std::endl;
-  std::cout << GridLogMessage << "* Generators for Sp(8)" << std::endl;
+  std::cout << GridLogMessage << "* Generators for Sp(8) (test)" << std::endl;
   std::cout << GridLogMessage << "*********************************************"
               << std::endl;
     
-  Sp8::printGenerators();
-  Sp8::testGenerators(); 
+  Sp8::testGenerators();
     
   std::cout << GridLogMessage << "*********************************************"
             << std::endl;
-  std::cout << GridLogMessage << "* Generators for Sp(4) TwoIndex AntiSymmetric" << std::endl;
+  std::cout << GridLogMessage << "* Generators for Sp(4) TwoIndexAS (test)" << std::endl;
   std::cout << GridLogMessage << "*********************************************"
             << std::endl;
 
-  Sp4TwoIndexAntiSymm::printGenerators();
-  Sp4TwoIndexAntiSymm::testGenerators();
+  Sp_TwoIndex<4, AntiSymmetric>:::testGenerators();
     
   std::cout << GridLogMessage << "*********************************************"
             << std::endl;
-  std::cout << GridLogMessage << "* Generators for Sp(6) TwoIndex AntiSymmetric" << std::endl;
+  std::cout << GridLogMessage << "* Generators for Sp(6) TwoIndexAS (test)" << std::endl;
   std::cout << GridLogMessage << "*********************************************"
             << std::endl;
-    
-  SU_TwoIndex<6, AntiSymmetric>::printGenerators();
-  SU_TwoIndex<6, AntiSymmetric>::testGenerators();
+
+  Sp_TwoIndex<6, AntiSymmetric>::testGenerators();
 
   std::cout << GridLogMessage << "*********************************************"
             << std::endl;
-  std::cout << GridLogMessage << "* Generators for Sp(8) TwoIndex AntiSymmetric" << std::endl;
+  std::cout << GridLogMessage << "* Generators for Sp(8) TwoIndexAS (test)" << std::endl;
   std::cout << GridLogMessage << "*********************************************"
             << std::endl;
-  
-  Sp_TwoIndex<8, AntiSymmetric>::printGenerators();
+
   Sp_TwoIndex<8, AntiSymmetric>::testGenerators();
 
   Grid_finalize();