diff --git a/Grid/qcd/utils/GaugeGroup.h b/Grid/qcd/utils/GaugeGroup.h index ecd92780..d0cdec17 100644 --- a/Grid/qcd/utils/GaugeGroup.h +++ b/Grid/qcd/utils/GaugeGroup.h @@ -170,8 +170,7 @@ class GaugeGroup { } } - // reunitarise?? - template + template static void LieRandomize(GridParallelRNG &pRNG, LatticeMatrixType &out, double scale = 1.0) { GridBase *grid = out.Grid(); @@ -282,6 +281,7 @@ class GaugeGroup { PokeIndex(out, Umu, mu); } } + template static void ColdConfiguration(GaugeField &out) { typedef typename GaugeField::vector_type vector_type; @@ -294,6 +294,7 @@ class GaugeGroup { PokeIndex(out, Umu, mu); } } + template static void ColdConfiguration(GridParallelRNG &pRNG, GaugeField &out) { ColdConfiguration(out); @@ -303,6 +304,7 @@ class GaugeGroup { static void taProj(const LatticeMatrixType &in, LatticeMatrixType &out) { out = Ta(in); } + template static void taExp(const LatticeMatrixType &x, LatticeMatrixType &ex) { typedef typename LatticeMatrixType::scalar_type ComplexType; @@ -320,6 +322,42 @@ class GaugeGroup { ex = ex + xn * nfac; // x2/2!, x3/3!.... } } + + template // reunitarise, resimplectify... + static void ProjectOnGaugeGroup(Lattice >, Nd> > &U) { + ProjectOnGaugeGroup(U, group_name()); + } + + template // reunitarise, resimplectify... + static void ProjectOnGaugeGroup(Lattice > > > &Umu) { + ProjectOnGaugeGroup(Umu, group_name()); + } + + + template // reunitarise, resimplectify... previously ProjectSUn + static void ProjectGn(Lattice > > > &Umu) { + ProjectOnGaugeGroup(Umu); + auto det = Determinant(Umu); + + det = conjugate(det); + + for (int i = 0; i < N; i++) { + auto element = PeekIndex(Umu, N - 1, i); + element = element * det; + PokeIndex(Umu, element, Nc - 1, i); + } + } + + template // reunitarise, resimplectify... previously ProjectSUn + static void ProjectGn(Lattice >, Nd> > &U, GroupName::SU) { + // Reunitarise + for (int mu = 0; mu < Nd; mu++) { + auto Umu = PeekIndex(U, mu); + ProjectGn(Umu); + PokeIndex(U, Umu, mu); + } + } + }; template @@ -347,6 +385,8 @@ LatticeComplexD Determinant( }); return ret; } + + template static void ProjectSUn( Lattice > > > &Umu) { @@ -361,6 +401,7 @@ static void ProjectSUn( PokeIndex(Umu, element, Nc - 1, i); } } + template static void ProjectSUn( Lattice >, Nd> > &U) { @@ -373,6 +414,7 @@ static void ProjectSUn( PokeIndex(U, Umu, mu); } } + // Explicit specialisation for SU(3). // Explicit specialisation for SU(3). static void ProjectSU3( diff --git a/Grid/qcd/utils/SUn.impl b/Grid/qcd/utils/SUn.impl index eaf82877..0e4b6124 100644 --- a/Grid/qcd/utils/SUn.impl +++ b/Grid/qcd/utils/SUn.impl @@ -508,6 +508,20 @@ static void testGenerators(GroupName::SU) { std::cout << GridLogMessage << std::endl; } +template +static void ProjectOnGaugeGroup(Lattice > > > &Umu, GroupName::SU) { + Umu = ProjectOnGroup(Umu); + } + +template +static void ProjectOnGaugeGroup(Lattice >, Nd> > &U, GroupName::SU) { + // Reunitarise + for (int mu = 0; mu < Nd; mu++) { + auto Umu = PeekIndex(U, mu); + Umu = ProjectOnGroup(Umu); + } +} + /* * Fundamental rep gauge xform */ diff --git a/Grid/qcd/utils/SUnAdjoint.h b/Grid/qcd/utils/SUnAdjoint.h index c03036ef..3a5954ca 100644 --- a/Grid/qcd/utils/SUnAdjoint.h +++ b/Grid/qcd/utils/SUnAdjoint.h @@ -51,14 +51,16 @@ public: typedef Lattice >, Nd> > LatticeAdjFieldF; typedef Lattice >, Nd> > LatticeAdjFieldD; + template + using iSUnMatrix = iScalar > >; template static void generator(int Index, iSUnAdjointMatrix &iAdjTa) { // returns i(T_Adj)^index necessary for the projectors // see definitions above iAdjTa = Zero(); - Vector::template iSUnMatrix > ta(ncolour * ncolour - 1); - typename SU::template iSUnMatrix tmp; + Vector > ta(ncolour * ncolour - 1); + iSUnMatrix tmp; // FIXME not very efficient to get all the generators everytime for (int a = 0; a < Dimension; a++) SU::generator(a, ta[a]); @@ -66,7 +68,7 @@ public: for (int a = 0; a < Dimension; a++) { tmp = ta[a] * ta[Index] - ta[Index] * ta[a]; for (int b = 0; b < (ncolour * ncolour - 1); b++) { - typename SU::template iSUnMatrix tmp1 = + iSUnMatrix tmp1 = 2.0 * tmp * ta[b]; // 2.0 from the normalization Complex iTr = TensorRemove(timesI(trace(tmp1))); //iAdjTa()()(b, a) = iTr; diff --git a/Grid/qcd/utils/SUnTwoIndex.h b/Grid/qcd/utils/SUnTwoIndex.h index 50ddf5c1..5e283de7 100644 --- a/Grid/qcd/utils/SUnTwoIndex.h +++ b/Grid/qcd/utils/SUnTwoIndex.h @@ -126,10 +126,10 @@ public: template static void generator(int Index, iSUnTwoIndexMatrix &i2indTa) { - Vector::template iSUnMatrix > ta( + Vector > ta( ncolour * ncolour - 1); - Vector::template iSUnMatrix > eij(Dimension); - typename SU::template iSUnMatrix tmp; + Vector > eij(Dimension); + iSUnMatrix tmp; i2indTa = Zero(); for (int a = 0; a < ncolour * ncolour - 1; a++) @@ -140,7 +140,7 @@ public: 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::template iSUnMatrix tmp1 = + iSUnMatrix tmp1 = tmp * eij[b]; Complex iTr = TensorRemove(timesI(trace(tmp1))); i2indTa()()(a, b) = iTr; diff --git a/Grid/qcd/utils/Sp2n.impl b/Grid/qcd/utils/Sp2n.impl index c8ed2ddc..77da9ad0 100644 --- a/Grid/qcd/utils/Sp2n.impl +++ b/Grid/qcd/utils/Sp2n.impl @@ -253,6 +253,20 @@ static void testGenerators(GroupName::Sp) { } } +template +static void ProjectOnGaugeGroup(Lattice > > > &Umu, GroupName::Sp) { + Umu = ProjectOnSpGroup(Umu); + } + +template +static void ProjectOnGaugeGroup(Lattice >, Nd> > &U, GroupName::Sp) { + // Reunitarise + for (int mu = 0; mu < Nd; mu++) { + auto Umu = PeekIndex(U, mu); + Umu = ProjectOnSpGroup(Umu); + } +} + public: template static void OmegaInvariance(ColourMatrix &in) { diff --git a/tests/core/Test_lie_generators.cc b/tests/core/Test_lie_generators.cc index 372e1408..5b044061 100644 --- a/tests/core/Test_lie_generators.cc +++ b/tests/core/Test_lie_generators.cc @@ -28,6 +28,7 @@ See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ + #include #include @@ -42,11 +43,11 @@ directory using namespace std; using namespace Grid; - ; + int main(int argc, char** argv) { Grid_init(&argc, &argv); - + std::vector latt({4, 4, 4, 8}); GridCartesian* grid = SpaceTimeGrid::makeFourDimGrid( latt, GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi()); @@ -66,7 +67,7 @@ int main(int argc, char** argv) { std::cout << GridLogMessage << "*********************************************" << std::endl; - std::cout << GridLogMessage << "* Generators for SU(Nc" << std::endl; + std::cout << GridLogMessage << "* Generators for SU(3)" << std::endl; std::cout << GridLogMessage << "*********************************************" << std::endl; SU3::printGenerators(); @@ -87,22 +88,22 @@ int main(int argc, char** argv) { // Projectors GridParallelRNG gridRNG(grid); gridRNG.SeedFixedIntegers(std::vector({45,12,81,9})); - SU3Adjoint::LatticeAdjMatrix Gauss(grid); - SU3::LatticeAlgebraVector ha(grid); - SU3::LatticeAlgebraVector hb(grid); + SU_Adjoint::LatticeAdjMatrix Gauss(grid); + SU::LatticeAlgebraVector ha(grid); + SU::LatticeAlgebraVector hb(grid); random(gridRNG,Gauss); std::cout << GridLogMessage << "Start projectOnAlgebra" << std::endl; - SU3Adjoint::projectOnAlgebra(ha, Gauss); + SU_Adjoint::projectOnAlgebra(ha, Gauss); std::cout << GridLogMessage << "end projectOnAlgebra" << std::endl; std::cout << GridLogMessage << "Start projector" << std::endl; - SU3Adjoint::projector(hb, Gauss); + SU_Adjoint::projector(hb, Gauss); std::cout << GridLogMessage << "end projector" << std::endl; std::cout << GridLogMessage << "ReStart projector" << std::endl; - SU3Adjoint::projector(hb, Gauss); + SU_Adjoint::projector(hb, Gauss); std::cout << GridLogMessage << "end projector" << std::endl; - SU3::LatticeAlgebraVector diff = ha -hb; + SU::LatticeAlgebraVector diff = ha -hb; std::cout << GridLogMessage << "Difference: " << norm2(diff) << std::endl; @@ -114,8 +115,8 @@ int main(int argc, char** argv) { LatticeGaugeField U(grid), V(grid); - SU3::HotConfiguration(gridRNG, U); - SU3::HotConfiguration(gridRNG, V); + SU::HotConfiguration(gridRNG, U); + SU::HotConfiguration(gridRNG, V); // Adjoint representation // Test group structure @@ -123,8 +124,8 @@ int main(int argc, char** argv) { LatticeGaugeField UV(grid); UV = Zero(); for (int mu = 0; mu < Nd; mu++) { - SU3::LatticeMatrix Umu = peekLorentz(U,mu); - SU3::LatticeMatrix Vmu = peekLorentz(V,mu); + SU::LatticeMatrix Umu = peekLorentz(U,mu); + SU::LatticeMatrix Vmu = peekLorentz(V,mu); pokeLorentz(UV,Umu*Vmu, mu); } @@ -151,16 +152,16 @@ int main(int argc, char** argv) { // Check correspondence of algebra and group transformations // Create a random vector - SU3::LatticeAlgebraVector h_adj(grid); + SU::LatticeAlgebraVector h_adj(grid); typename AdjointRep::LatticeMatrix Ar(grid); random(gridRNG,h_adj); h_adj = real(h_adj); SU_Adjoint::AdjointLieAlgebraMatrix(h_adj,Ar); // Re-extract h_adj - SU3::LatticeAlgebraVector h_adj2(grid); + SU::LatticeAlgebraVector h_adj2(grid); SU_Adjoint::projectOnAlgebra(h_adj2, Ar); - SU3::LatticeAlgebraVector h_diff = h_adj - h_adj2; + SU::LatticeAlgebraVector h_diff = h_adj - h_adj2; std::cout << GridLogMessage << "Projections structure check vector difference (Adjoint representation) : " << norm2(h_diff) << std::endl; // Exponentiate @@ -183,14 +184,14 @@ int main(int argc, char** argv) { // Construct the fundamental matrix in the group - SU3::LatticeMatrix Af(grid); - SU3::FundamentalLieAlgebraMatrix(h_adj,Af); - SU3::LatticeMatrix Ufund(grid); + SU::LatticeMatrix Af(grid); + SU::FundamentalLieAlgebraMatrix(h_adj,Af); + SU::LatticeMatrix Ufund(grid); Ufund = expMat(Af, 1.0, 16); // Check unitarity - SU3::LatticeMatrix uno_f(grid); + SU::LatticeMatrix uno_f(grid); uno_f = 1.0; - SU3::LatticeMatrix UnitCheck(grid); + SU::LatticeMatrix UnitCheck(grid); UnitCheck = Ufund * adj(Ufund) - uno_f; std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck) << std::endl; @@ -260,20 +261,20 @@ int main(int argc, char** argv) { std::cout << GridLogMessage << "Test for the Two Index Symmetric projectors" << std::endl; // Projectors - SU3TwoIndexSymm::LatticeTwoIndexMatrix Gauss2(grid); + SU_TwoIndex::LatticeTwoIndexMatrix Gauss2(grid); random(gridRNG,Gauss2); std::cout << GridLogMessage << "Start projectOnAlgebra" << std::endl; - SU3TwoIndexSymm::projectOnAlgebra(ha, Gauss2); + SU_TwoIndex::projectOnAlgebra(ha, Gauss2); std::cout << GridLogMessage << "end projectOnAlgebra" << std::endl; std::cout << GridLogMessage << "Start projector" << std::endl; - SU3TwoIndexSymm::projector(hb, Gauss2); + SU_TwoIndex::projector(hb, Gauss2); std::cout << GridLogMessage << "end projector" << std::endl; std::cout << GridLogMessage << "ReStart projector" << std::endl; - SU3TwoIndexSymm::projector(hb, Gauss2); + SU_TwoIndex::projector(hb, Gauss2); std::cout << GridLogMessage << "end projector" << std::endl; - SU3::LatticeAlgebraVector diff2 = ha - hb; + SU::LatticeAlgebraVector diff2 = ha - hb; std::cout << GridLogMessage << "Difference: " << norm2(diff) << std::endl; std::cout << GridLogMessage << "*********************************************" << std::endl; @@ -284,20 +285,20 @@ int main(int argc, char** argv) { std::cout << GridLogMessage << "Test for the Two index anti-Symmetric projectors" << std::endl; // Projectors - SU3TwoIndexAntiSymm::LatticeTwoIndexMatrix Gauss2a(grid); + SU_TwoIndex::LatticeTwoIndexMatrix Gauss2a(grid); random(gridRNG,Gauss2a); std::cout << GridLogMessage << "Start projectOnAlgebra" << std::endl; - SU3TwoIndexAntiSymm::projectOnAlgebra(ha, Gauss2a); + SU_TwoIndex::projectOnAlgebra(ha, Gauss2a); std::cout << GridLogMessage << "end projectOnAlgebra" << std::endl; std::cout << GridLogMessage << "Start projector" << std::endl; - SU3TwoIndexAntiSymm::projector(hb, Gauss2a); + SU_TwoIndex::projector(hb, Gauss2a); std::cout << GridLogMessage << "end projector" << std::endl; std::cout << GridLogMessage << "ReStart projector" << std::endl; - SU3TwoIndexAntiSymm::projector(hb, Gauss2a); + SU_TwoIndex::projector(hb, Gauss2a); std::cout << GridLogMessage << "end projector" << std::endl; - SU3::LatticeAlgebraVector diff2a = ha - hb; + SU::LatticeAlgebraVector diff2a = ha - hb; std::cout << GridLogMessage << "Difference: " << norm2(diff2a) << std::endl; std::cout << GridLogMessage << "*********************************************" << std::endl; @@ -311,14 +312,14 @@ int main(int argc, char** argv) { // Test group structure // (U_f * V_f)_r = U_r * V_r LatticeGaugeField U2(grid), V2(grid); - SU3::HotConfiguration(gridRNG, U2); - SU3::HotConfiguration(gridRNG, V2); + SU::HotConfiguration(gridRNG, U2); + SU::HotConfiguration(gridRNG, V2); LatticeGaugeField UV2(grid); UV2 = Zero(); for (int mu = 0; mu < Nd; mu++) { - SU3::LatticeMatrix Umu2 = peekLorentz(U2,mu); - SU3::LatticeMatrix Vmu2 = peekLorentz(V2,mu); + SU::LatticeMatrix Umu2 = peekLorentz(U2,mu); + SU::LatticeMatrix Vmu2 = peekLorentz(V2,mu); pokeLorentz(UV2,Umu2*Vmu2, mu); } @@ -345,16 +346,16 @@ int main(int argc, char** argv) { // Check correspondence of algebra and group transformations // Create a random vector - SU3::LatticeAlgebraVector h_sym(grid); + SU::LatticeAlgebraVector h_sym(grid); typename TwoIndexRep< Nc, Symmetric>::LatticeMatrix Ar_sym(grid); random(gridRNG,h_sym); h_sym = real(h_sym); SU_TwoIndex::TwoIndexLieAlgebraMatrix(h_sym,Ar_sym); // Re-extract h_sym - SU3::LatticeAlgebraVector h_sym2(grid); + SU::LatticeAlgebraVector h_sym2(grid); SU_TwoIndex< Nc, Symmetric>::projectOnAlgebra(h_sym2, Ar_sym); - SU3::LatticeAlgebraVector h_diff_sym = h_sym - h_sym2; + SU::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; @@ -379,11 +380,11 @@ int main(int argc, char** argv) { // Construct the fundamental matrix in the group - SU3::LatticeMatrix Af_sym(grid); - SU3::FundamentalLieAlgebraMatrix(h_sym,Af_sym); - SU3::LatticeMatrix Ufund2(grid); + SU::LatticeMatrix Af_sym(grid); + SU::FundamentalLieAlgebraMatrix(h_sym,Af_sym); + SU::LatticeMatrix Ufund2(grid); Ufund2 = expMat(Af_sym, 1.0, 16); - SU3::LatticeMatrix UnitCheck2(grid); + SU::LatticeMatrix UnitCheck2(grid); UnitCheck2 = Ufund2 * adj(Ufund2) - uno_f; std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck2) << std::endl; @@ -421,14 +422,14 @@ int main(int argc, char** argv) { // Test group structure // (U_f * V_f)_r = U_r * V_r LatticeGaugeField U2A(grid), V2A(grid); - SU3::HotConfiguration(gridRNG, U2A); - SU3::HotConfiguration(gridRNG, V2A); + SU::HotConfiguration(gridRNG, U2A); + SU::HotConfiguration(gridRNG, V2A); LatticeGaugeField UV2A(grid); UV2A = Zero(); for (int mu = 0; mu < Nd; mu++) { - SU3::LatticeMatrix Umu2A = peekLorentz(U2,mu); - SU3::LatticeMatrix Vmu2A = peekLorentz(V2,mu); + SU::LatticeMatrix Umu2A = peekLorentz(U2,mu); + SU::LatticeMatrix Vmu2A = peekLorentz(V2,mu); pokeLorentz(UV2A,Umu2A*Vmu2A, mu); } @@ -455,16 +456,16 @@ int main(int argc, char** argv) { // Check correspondence of algebra and group transformations // Create a random vector - SU3::LatticeAlgebraVector h_Asym(grid); + SU::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 - SU3::LatticeAlgebraVector h_Asym2(grid); + SU::LatticeAlgebraVector h_Asym2(grid); SU_TwoIndex< Nc, AntiSymmetric>::projectOnAlgebra(h_Asym2, Ar_Asym); - SU3::LatticeAlgebraVector h_diff_Asym = h_Asym - h_Asym2; + SU::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; @@ -489,11 +490,11 @@ int main(int argc, char** argv) { // Construct the fundamental matrix in the group - SU3::LatticeMatrix Af_Asym(grid); - SU3::FundamentalLieAlgebraMatrix(h_Asym,Af_Asym); - SU3::LatticeMatrix Ufund2A(grid); + SU::LatticeMatrix Af_Asym(grid); + SU::FundamentalLieAlgebraMatrix(h_Asym,Af_Asym); + SU::LatticeMatrix Ufund2A(grid); Ufund2A = expMat(Af_Asym, 1.0, 16); - SU3::LatticeMatrix UnitCheck2A(grid); + SU::LatticeMatrix UnitCheck2A(grid); UnitCheck2A = Ufund2A * adj(Ufund2A) - uno_f; std::cout << GridLogMessage << "unitarity check 1: " << norm2(UnitCheck2A) << std::endl; @@ -526,7 +527,7 @@ int main(int argc, char** argv) { - Grid_finalize(); } + diff --git a/tests/core/Test_reunitarise.cc b/tests/core/Test_reunitarise.cc index 6644be1a..e421a2a8 100644 --- a/tests/core/Test_reunitarise.cc +++ b/tests/core/Test_reunitarise.cc @@ -26,6 +26,7 @@ Author: Peter Boyle See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ + #include using namespace std; @@ -35,7 +36,7 @@ using namespace Grid; int main (int argc, char ** argv) { Grid_init(&argc,&argv); - + std::vector latt({8,8,8,8}); GridCartesian * grid = SpaceTimeGrid::makeFourDimGrid(latt, GridDefaultSimd(Nd,vComplexD::Nsimd()), @@ -122,14 +123,14 @@ int main (int argc, char ** argv) std::cout << "Determinant defect before projection " <::HotConfiguration(pRNG_f, Umu); RealD checkTolerance = (getPrecision::value == 1) ? 1e-7 : 1e-15;