portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-10-29 19:14:33 +00:00
Code Releases Activity

Merge pull request #22 from LupoA/sp2n/unify_twoindex

Browse Source 
Unify TwoIndex
This commit is contained in:
chillenzer
2023-05-18 14:55:02 +00:00
committed by GitHub
parent aa9df63a05 a5125e23f4
commit d169c275b6
18 changed files with 319 additions and 443 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexSymmetricImplD/WilsonCloverFermionInstantiationSpWilsonTwoIndexSymmetricImplD.cc
View File

@@ -1 +0,0 @@
../WilsonCloverFermionInstantiation.cc.master

1
Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexSymmetricImplD/WilsonFermionInstantiationSpWilsonTwoIndexSymmetricImplD.cc
View File

@@ -1 +0,0 @@
../WilsonFermionInstantiation.cc.master

1
Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexSymmetricImplD/WilsonKernelsInstantiationSpWilsonTwoIndexSymmetricImplD.cc
View File

@@ -1 +0,0 @@
../WilsonKernelsInstantiation.cc.master

1
Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexSymmetricImplD/WilsonTMFermionInstantiationSpWilsonTwoIndexSymmetricImplD.cc
View File

@@ -1 +0,0 @@
../WilsonTMFermionInstantiation.cc.master

1
Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexSymmetricImplD/impl.h
View File

@@ -1 +0,0 @@
#define IMPLEMENTATION SpWilsonTwoIndexSymmetricImplD

1
Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexSymmetricImplF/WilsonCloverFermionInstantiationSpWilsonTwoIndexSymmetricImplF.cc
View File

@@ -1 +0,0 @@
../WilsonCloverFermionInstantiation.cc.master

1
Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexSymmetricImplF/WilsonFermionInstantiationSpWilsonTwoIndexSymmetricImplF.cc
View File

@@ -1 +0,0 @@
../WilsonFermionInstantiation.cc.master

1
Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexSymmetricImplF/WilsonKernelsInstantiationSpWilsonTwoIndexSymmetricImplF.cc
View File

@@ -1 +0,0 @@
../WilsonKernelsInstantiation.cc.master

1
Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexSymmetricImplF/WilsonTMFermionInstantiationSpWilsonTwoIndexSymmetricImplF.cc
View File

@@ -1 +0,0 @@
../WilsonTMFermionInstantiation.cc.master

1
Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexSymmetricImplF/impl.h
View File

@@ -1 +0,0 @@
#define IMPLEMENTATION SpWilsonTwoIndexSymmetricImplF

2
Grid/qcd/representations/sp_two_index.h
View File

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

217
Grid/qcd/utils/Sp2nTwoIndex.h → Grid/qcd/utils/GaugeGroupTwoIndex.h
View File

@@ -22,71 +22,84 @@
// Authors: David Preti, Guido Cossu // Authors: David Preti, Guido Cossu
#ifndef QCD_UTIL_SP2N2INDEX_H #ifndef QCD_UTIL_SUN2INDEX_H
#define QCD_UTIL_SP2N2INDEX_H #define QCD_UTIL_SUN2INDEX_H
NAMESPACE_BEGIN(Grid); NAMESPACE_BEGIN(Grid);
//enum SpTwoIndexSymmetry { S = 1, AS = -1 }; enum TwoIndexSymmetry { Symmetric = 1, AntiSymmetric = -1 };
//inline Real delta(int a, int b) { return (a == b) ? 1.0 : 0.0; } inline Real delta(int a, int b) { return (a == b) ? 1.0 : 0.0; }
template <int ncolour, TwoIndexSymmetry S> template <int ncolour, TwoIndexSymmetry S, class group_name>
class Sp_TwoIndex : public Sp<ncolour> { class GaugeGroupTwoIndex : public GaugeGroup<ncolour, group_name> {
public: public:
static const int nsp = ncolour/2; static_assert(
static const int Dimension = nsp * (ncolour + S) - 1 ; std::is_same<group_name, GroupName::Sp>::value ? S != Symmetric : true,
static const int NumGenerators = Sp<ncolour>::AlgebraDimension; "The symmetric two-index representation of Sp(2N) does not work "
"currently. If you want to use it, you need to implement the equivalent "
"of Eq. (27) and (28) from https://doi.org/10.48550/arXiv.2202.05516.");
// The chosen convention is that we are taking ncolour to be N in SU<N> but 2N
// in Sp(2N). ngroup is equal to N for SU but 2N/2 = N for Sp(2N).
static_assert(std::is_same<group_name, GroupName::SU>::value or
std::is_same<group_name, GroupName::Sp>::value,
"ngroup is only implemented for SU and Sp currently.");
static const int ngroup =
std::is_same<group_name, GroupName::SU>::value ? ncolour : ncolour / 2;
static const int Dimension = std::is_same<group_name, GroupName::SU>::value
? ncolour * (ncolour + S) / 2
: ngroup * (ncolour + S) + S;
static const int NumGenerators =
GaugeGroup<ncolour, group_name>::AlgebraDimension;
template <typename vtype> template <typename vtype>
using iSp2nTwoIndexMatrix = iScalar<iScalar<iMatrix<vtype, Dimension> > >; using iGroupTwoIndexMatrix = iScalar<iScalar<iMatrix<vtype, Dimension> > >;
typedef iSp2nTwoIndexMatrix<Complex> TIMatrix; typedef iGroupTwoIndexMatrix<Complex> TIMatrix;
typedef iSp2nTwoIndexMatrix<ComplexF> TIMatrixF; typedef iGroupTwoIndexMatrix<ComplexF> TIMatrixF;
typedef iSp2nTwoIndexMatrix<ComplexD> TIMatrixD; typedef iGroupTwoIndexMatrix<ComplexD> TIMatrixD;
typedef iSp2nTwoIndexMatrix<vComplex> vTIMatrix; typedef iGroupTwoIndexMatrix<vComplex> vTIMatrix;
typedef iSp2nTwoIndexMatrix<vComplexF> vTIMatrixF; typedef iGroupTwoIndexMatrix<vComplexF> vTIMatrixF;
typedef iSp2nTwoIndexMatrix<vComplexD> vTIMatrixD; typedef iGroupTwoIndexMatrix<vComplexD> vTIMatrixD;
typedef Lattice<vTIMatrix> LatticeTwoIndexMatrix; typedef Lattice<vTIMatrix> LatticeTwoIndexMatrix;
typedef Lattice<vTIMatrixF> LatticeTwoIndexMatrixF; typedef Lattice<vTIMatrixF> LatticeTwoIndexMatrixF;
typedef Lattice<vTIMatrixD> LatticeTwoIndexMatrixD; typedef Lattice<vTIMatrixD> LatticeTwoIndexMatrixD;
typedef Lattice<iVector<iScalar<iMatrix<vComplex, Dimension> >, Nd> > typedef Lattice<iVector<iScalar<iMatrix<vComplex, Dimension> >, Nd> >
LatticeTwoIndexField; LatticeTwoIndexField;
typedef Lattice<iVector<iScalar<iMatrix<vComplexF, Dimension> >, Nd> > typedef Lattice<iVector<iScalar<iMatrix<vComplexF, Dimension> >, Nd> >
LatticeTwoIndexFieldF; LatticeTwoIndexFieldF;
typedef Lattice<iVector<iScalar<iMatrix<vComplexD, Dimension> >, Nd> > typedef Lattice<iVector<iScalar<iMatrix<vComplexD, Dimension> >, Nd> >
LatticeTwoIndexFieldD; LatticeTwoIndexFieldD;
template <typename vtype> template <typename vtype>
using iSp2nMatrix = iScalar<iScalar<iMatrix<vtype, ncolour> > >; using iGroupMatrix = iScalar<iScalar<iMatrix<vtype, ncolour> > >;
typedef iSp2nMatrix<Complex> Matrix; typedef iGroupMatrix<Complex> Matrix;
typedef iSp2nMatrix<ComplexF> MatrixF; typedef iGroupMatrix<ComplexF> MatrixF;
typedef iSp2nMatrix<ComplexD> MatrixD; typedef iGroupMatrix<ComplexD> MatrixD;
template <class cplx> template <class cplx>
static void base(int Index, iSp2nMatrix<cplx> &eij) { static void base(int Index, iGroupMatrix<cplx> &eij) {
// returns (e)^(ij)_{kl} necessary for change of base U_F -> U_R // returns (e)^(ij)_{kl} necessary for change of base U_F -> U_R
assert(Index < NumGenerators); assert(Index < Dimension);
eij = Zero(); eij = Zero();
// for the linearisation of the 2 indexes // for the linearisation of the 2 indexes
static int a[ncolour * (ncolour - 1) / 2][2]; // store the a <-> i,j static int a[ncolour * (ncolour - 1) / 2][2]; // store the a <-> i,j
static bool filled = false; static bool filled = false;
if (!filled) { if (!filled) {
int counter = 0; int counter = 0;
for (int i = 1; i < ncolour; i++) { for (int i = 1; i < ncolour; i++) {
for (int j = 0; j < i; j++) { for (int j = 0; j < i; j++) {
a[counter][0] = i;
if (i==nsp) if (j==0 && ngroup == ncolour/2 && i==ngroup+j) {
{ //std::cout << "skipping" << std::endl; // for Sp2n this vanishes identically.
j = j+1; j = j+1;
} }
a[counter][0] = i;
a[counter][1] = j; a[counter][1] = j;
counter++; counter++;
} }
@@ -95,69 +108,65 @@ public:
} }
if (Index < ncolour * (ncolour - 1) / 2) { if (Index < ncolour * (ncolour - 1) / 2) {
baseOffDiagonal(a[Index][0], a[Index][1], eij); baseOffDiagonal(a[Index][0], a[Index][1], eij, group_name());
} else { } else {
baseDiagonal(Index, eij); baseDiagonal(Index, eij);
} }
} }
template <class cplx> template <class cplx>
static void baseDiagonal(int Index, iSp2nMatrix<cplx> &eij) { static void baseDiagonal(int Index, iGroupMatrix<cplx> &eij) {
eij = Zero(); eij = Zero();
eij()()(Index - ncolour * (ncolour - 1) / 2, eij()()(Index - ncolour * (ncolour - 1) / 2,
Index - ncolour * (ncolour - 1) / 2) = 1.0; Index - ncolour * (ncolour - 1) / 2) = 1.0;
} }
template <class cplx> template <class cplx>
static void baseOffDiagonal(int i, int j, iSp2nMatrix<cplx> &eij) { static void baseOffDiagonal(int i, int j, iGroupMatrix<cplx> &eij,
GroupName::Sp) {
eij = Zero(); eij = Zero();
//std::cout << GridLogMessage << " doing i j = " << i << j << std::endl; RealD tmp;
RealD tmp;
if ( (i == nsp + j) && (1 <= j) && (j < nsp) )
{
for (int k = 0; k < nsp; k++)
{
if (k < j)
{
//std::cout << GridLogMessage << "b=N+a 1"<< std::endl;
//std::cout << GridLogMessage << "j i "<< j << i << std::endl;
tmp = sqrt( 2*j*(j+1) );
tmp = 1/tmp;
tmp *= std::sqrt(2.0);
eij()()(k,k+nsp) = tmp;
eij()()(k+nsp,k) = - tmp;
}
if (k == j)
{
//std::cout << GridLogMessage << "b=N+a 2"<< std::endl;
//std::cout << GridLogMessage << "j i "<< j << i << std::endl;
tmp = sqrt(2*j*(j+1) );
tmp = -j/tmp;
tmp *= std::sqrt(2.0);
eij()()(k,k+nsp) = tmp ;
eij()()(k+nsp,k) = - tmp ;
}
}
if ((i == ngroup + j) && (1 <= j) && (j < ngroup)) {
for (int k = 0; k < ngroup; k++) {
if (k < j) {
tmp = sqrt(2 * j * (j + 1));
tmp = 1 / tmp;
tmp *= std::sqrt(2.0);
eij()()(k, k + ngroup) = tmp;
eij()()(k + ngroup, k) = -tmp;
} }
if (k == j) {
tmp = sqrt(2 * j * (j + 1));
else if (i != nsp+j) tmp = -j / tmp;
{ tmp *= std::sqrt(2.0);
for (int k = 0; k < ncolour; k++) eij()()(k, k + ngroup) = tmp;
for (int l = 0; l < ncolour; l++) eij()()(k + ngroup, k) = -tmp;
{
eij()()(l, k) = delta(i, k) * delta(j, l) + S * delta(j, k) * delta(i, l);
}
} }
}
}
else if (i != ngroup + j) {
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;
}
template <class cplx>
static void baseOffDiagonal(int i, int j, iGroupMatrix<cplx> &eij,
GroupName::SU) {
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); RealD nrm = 1. / std::sqrt(2.0);
eij = eij * nrm; eij = eij * nrm;
@@ -174,23 +183,21 @@ public:
} }
template <class cplx> template <class cplx>
static void generator(int Index, iSp2nTwoIndexMatrix<cplx> &i2indTa) { static void generator(int Index, iGroupTwoIndexMatrix<cplx> &i2indTa) {
Vector<iSp2nMatrix<cplx> > ta( Vector<iGroupMatrix<cplx> > ta(NumGenerators);
NumGenerators); Vector<iGroupMatrix<cplx> > eij(Dimension);
Vector<iSp2nMatrix<cplx> > eij(Dimension); iGroupMatrix<cplx> tmp;
iSp2nMatrix<cplx> tmp;
i2indTa = Zero(); i2indTa = Zero();
for (int a = 0; a < NumGenerators; a++) for (int a = 0; a < NumGenerators; a++)
Sp<ncolour>::generator(a, ta[a]); GaugeGroup<ncolour, group_name>::generator(a, ta[a]);
for (int a = 0; a < Dimension; a++) base(a, eij[a]); for (int a = 0; a < Dimension; a++) base(a, eij[a]);
for (int a = 0; a < Dimension; a++) { for (int a = 0; a < Dimension; a++) {
tmp = transpose(ta[Index]) * adj(eij[a]) + adj(eij[a]) * ta[Index]; tmp = transpose(ta[Index]) * adj(eij[a]) + adj(eij[a]) * ta[Index];
for (int b = 0; b < Dimension; b++) { for (int b = 0; b < Dimension; b++) {
iSp2nMatrix<cplx> tmp1 = iGroupMatrix<cplx> tmp1 = tmp * eij[b];
tmp * eij[b];
Complex iTr = TensorRemove(timesI(trace(tmp1))); Complex iTr = TensorRemove(timesI(trace(tmp1)));
i2indTa()()(a, b) = iTr; i2indTa()()(a, b) = iTr;
} }
@@ -239,14 +246,21 @@ public:
Complex Tr = -TensorRemove(trace(i2indTa * i2indTb)); Complex Tr = -TensorRemove(trace(i2indTa * i2indTb));
std::cout << GridLogMessage << "a=" << a << "b=" << b << "Tr=" << Tr std::cout << GridLogMessage << "a=" << a << "b=" << b << "Tr=" << Tr
<< std::endl; << std::endl;
if (a == b) {
assert(imag(Tr) < 1e-8);
assert(real(Tr) - ((ncolour + S * 2) * 0.5) < 1e-8);
} else {
assert(imag(Tr) < 1e-8);
assert(real(Tr) < 1e-8);
}
} }
} }
std::cout << GridLogMessage << std::endl; std::cout << GridLogMessage << std::endl;
} }
static void TwoIndexLieAlgebraMatrix( static void TwoIndexLieAlgebraMatrix(
const typename Sp<ncolour>::LatticeAlgebraVector &h, const typename GaugeGroup<ncolour, group_name>::LatticeAlgebraVector &h,
LatticeTwoIndexMatrix &out, Real scale = 1.0) { LatticeTwoIndexMatrix &out, Real scale = 1.0) {
conformable(h, out); conformable(h, out);
GridBase *grid = out.Grid(); GridBase *grid = out.Grid();
LatticeTwoIndexMatrix la(grid); LatticeTwoIndexMatrix la(grid);
@@ -264,8 +278,8 @@ public:
// Projects the algebra components // Projects the algebra components
// of a lattice matrix ( of dimension ncol*ncol -1 ) // of a lattice matrix ( of dimension ncol*ncol -1 )
static void projectOnAlgebra( static void projectOnAlgebra(
typename Sp<ncolour>::LatticeAlgebraVector &h_out, typename GaugeGroup<ncolour, group_name>::LatticeAlgebraVector &h_out,
const LatticeTwoIndexMatrix &in, Real scale = 1.0) { const LatticeTwoIndexMatrix &in, Real scale = 1.0) {
conformable(h_out, in); conformable(h_out, in);
h_out = Zero(); h_out = Zero();
TIMatrix i2indTa; TIMatrix i2indTa;
@@ -273,14 +287,15 @@ public:
// 2/(Nc +/- 2) for the normalization of the trace in the two index rep // 2/(Nc +/- 2) for the normalization of the trace in the two index rep
for (int a = 0; a < NumGenerators; a++) { for (int a = 0; a < NumGenerators; a++) {
generator(a, i2indTa); generator(a, i2indTa);
pokeColour(h_out, real(trace(i2indTa * in)) * coefficient, a); pokeColour(h_out, real(trace(i2indTa * in)) * coefficient, a);
} }
} }
// a projector that keeps the generators stored to avoid the overhead of // a projector that keeps the generators stored to avoid the overhead of
// recomputing them // recomputing them
static void projector(typename Sp<ncolour>::LatticeAlgebraVector &h_out, static void projector(
const LatticeTwoIndexMatrix &in, Real scale = 1.0) { typename GaugeGroup<ncolour, group_name>::LatticeAlgebraVector &h_out,
const LatticeTwoIndexMatrix &in, Real scale = 1.0) {
conformable(h_out, in); conformable(h_out, in);
// to store the generators // to store the generators
static std::vector<TIMatrix> i2indTa(NumGenerators); static std::vector<TIMatrix> i2indTa(NumGenerators);
@@ -292,7 +307,7 @@ public:
} }
Real coefficient = Real coefficient =
-2.0 / (ncolour + 2 * S) * scale; // 2/(Nc +/- 2) for the normalization -2.0 / (ncolour + 2 * S) * scale; // 2/(Nc +/- 2) for the normalization
// of the trace in the two index rep // of the trace in the two index rep
for (int a = 0; a < NumGenerators; a++) { for (int a = 0; a < NumGenerators; a++) {
@@ -302,18 +317,34 @@ public:
} }
}; };
template <int ncolour, TwoIndexSymmetry S>
using SU_TwoIndex = GaugeGroupTwoIndex<ncolour, S, GroupName::SU>;
// Some useful type names // 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<2, AntiSymmetric> SU2TwoIndexAntiSymm;
typedef SU_TwoIndex<3, AntiSymmetric> SU3TwoIndexAntiSymm;
typedef SU_TwoIndex<4, AntiSymmetric> SU4TwoIndexAntiSymm;
typedef SU_TwoIndex<5, AntiSymmetric> SU5TwoIndexAntiSymm;
template <int ncolour, TwoIndexSymmetry S>
using Sp_TwoIndex = GaugeGroupTwoIndex<ncolour, S, GroupName::Sp>;
typedef Sp_TwoIndex<Nc, Symmetric> SpTwoIndexSymmMatrices; typedef Sp_TwoIndex<Nc, Symmetric> SpTwoIndexSymmMatrices;
typedef Sp_TwoIndex<Nc, AntiSymmetric> SpTwoIndexAntiSymmMatrices; typedef Sp_TwoIndex<Nc, AntiSymmetric> SpTwoIndexAntiSymmMatrices;
typedef Sp_TwoIndex<2, Symmetric> Sp2TwoIndexSymm; typedef Sp_TwoIndex<2, Symmetric> Sp2TwoIndexSymm;
typedef Sp_TwoIndex<4, Symmetric> Sp4TwoIndexSymm; typedef Sp_TwoIndex<4, Symmetric> Sp4TwoIndexSymm;
typedef Sp_TwoIndex<2, AntiSymmetric> Sp2TwoIndexAntiSymm;
typedef Sp_TwoIndex<4, AntiSymmetric> Sp4TwoIndexAntiSymm; typedef Sp_TwoIndex<4, AntiSymmetric> Sp4TwoIndexAntiSymm;
NAMESPACE_END(Grid); NAMESPACE_END(Grid);
#endif #endif

271
Grid/qcd/utils/SUnTwoIndex.h
View File

@@ -1,271 +0,0 @@
////////////////////////////////////////////////////////////////////////
//
// * Two index representation generators
//
// * Normalisation for the fundamental generators:
// trace ta tb = 1/2 delta_ab = T_F delta_ab
// T_F = 1/2 for SU(N) groups
//
//
// base for NxN two index (anti-symmetric) matrices
// normalized to 1 (d_ij is the kroenecker delta)
//
// (e^(ij)_{kl} = 1 / sqrt(2) (d_ik d_jl +/- d_jk d_il)
//
// Then the generators are written as
//
// (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_BEGIN(Grid);
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 = iScalar<iScalar<iMatrix<vtype, Dimension> > >;
typedef iSUnTwoIndexMatrix<Complex> TIMatrix;
typedef iSUnTwoIndexMatrix<ComplexF> TIMatrixF;
typedef iSUnTwoIndexMatrix<ComplexD> TIMatrixD;
typedef iSUnTwoIndexMatrix<vComplex> vTIMatrix;
typedef iSUnTwoIndexMatrix<vComplexF> vTIMatrixF;
typedef iSUnTwoIndexMatrix<vComplexD> vTIMatrixD;
typedef Lattice<vTIMatrix> LatticeTwoIndexMatrix;
typedef Lattice<vTIMatrixF> LatticeTwoIndexMatrixF;
typedef Lattice<vTIMatrixD> LatticeTwoIndexMatrixD;
typedef Lattice<iVector<iScalar<iMatrix<vComplex, Dimension> >, Nd> >
LatticeTwoIndexField;
typedef Lattice<iVector<iScalar<iMatrix<vComplexF, Dimension> >, Nd> >
LatticeTwoIndexFieldF;
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 base(int Index, iSUnMatrix<cplx> &eij) {
// returns (e)^(ij)_{kl} necessary for change of base U_F -> U_R
assert(Index < NumGenerators);
eij = Zero();
// 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;
counter++;
}
}
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<iSUnMatrix<cplx> > ta(ncolour * ncolour - 1);
Vector<iSUnMatrix<cplx> > eij(Dimension);
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++) {
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 < ncolour * ncolour - 1; gen++) {
TIMatrix i2indTa;
generator(gen, i2indTa);
std::cout << GridLogMessage << "Nc = " << ncolour << " t_" << gen
<< std::endl;
std::cout << GridLogMessage << i2indTa << std::endl;
}
}
static void testGenerators(void) {
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);
}
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();
LatticeTwoIndexMatrix la(grid);
TIMatrix i2indTa;
out = Zero();
for (int a = 0; a < ncolour * ncolour - 1; a++) {
generator(a, i2indTa);
la = peekColour(h, a) * i2indTa;
out += la;
}
out *= scale;
}
// Projects the algebra components
// 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 i2indTa;
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 < ncolour * ncolour - 1; a++) {
generator(a, i2indTa);
pokeColour(h_out, real(trace(i2indTa * in)) * coefficient, a);
}
}
// a projector that keeps the generators stored to avoid the overhead of
// recomputing them
static void projector(typename SU<ncolour>::LatticeAlgebraVector &h_out,
const LatticeTwoIndexMatrix &in, Real scale = 1.0) {
conformable(h_out, in);
// 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 < 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
for (int a = 0; a < ncolour * ncolour - 1; a++) {
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<2, AntiSymmetric> SU2TwoIndexAntiSymm;
typedef SU_TwoIndex<3, AntiSymmetric> SU3TwoIndexAntiSymm;
typedef SU_TwoIndex<4, AntiSymmetric> SU4TwoIndexAntiSymm;
typedef SU_TwoIndex<5, AntiSymmetric> SU5TwoIndexAntiSymm;
NAMESPACE_END(Grid);
#endif

3
Grid/qcd/utils/Utils.h
View File

@@ -10,8 +10,7 @@
// Include representations // Include representations
#include <Grid/qcd/utils/GaugeGroup.h> #include <Grid/qcd/utils/GaugeGroup.h>
#include <Grid/qcd/utils/SUnAdjoint.h> #include <Grid/qcd/utils/SUnAdjoint.h>
#include <Grid/qcd/utils/SUnTwoIndex.h> #include <Grid/qcd/utils/GaugeGroupTwoIndex.h>
#include <Grid/qcd/utils/Sp2nTwoIndex.h>
// All-to-all contraction kernels that touch the // All-to-all contraction kernels that touch the
// internal lattice structure // internal lattice structure

2
tests/core/Test_lie_generators.cc
View File

@@ -35,7 +35,7 @@ directory
#include <Grid/qcd/utils/GaugeGroup.h> #include <Grid/qcd/utils/GaugeGroup.h>
#include <Grid/qcd/utils/SUnAdjoint.h> #include <Grid/qcd/utils/SUnAdjoint.h>
#include <Grid/qcd/utils/SUnTwoIndex.h> #include <Grid/qcd/utils/GaugeGroupTwoIndex.h>
#include <Grid/qcd/representations/adjoint.h> #include <Grid/qcd/representations/adjoint.h>
#include <Grid/qcd/representations/two_index.h> #include <Grid/qcd/representations/two_index.h>

110
tests/sp2n/Test_2as_base.cc Normal file
View File

@@ -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();
}

38
tests/sp2n/Test_hmc_Sp_Wilson2ASFermionGauge.cc
View File

@@ -3,11 +3,14 @@
int main(int argc, char **argv) { int main(int argc, char **argv) {
using namespace Grid; using namespace Grid;
typedef Representations< SpFundamentalRepresentation, SpTwoIndexAntiSymmetricRepresentation > TheRepresentations; typedef Representations<SpFundamentalRepresentation,
SpTwoIndexAntiSymmetricRepresentation>
TheRepresentations;
Grid_init(&argc, &argv); Grid_init(&argc, &argv);
typedef GenericSp2nHMCRunnerHirep<TheRepresentations, MinimumNorm2> HMCWrapper; typedef GenericSp2nHMCRunnerHirep<TheRepresentations, MinimumNorm2>
HMCWrapper;
typedef SpWilsonTwoIndexAntiSymmetricImplR FermionImplPolicy; typedef SpWilsonTwoIndexAntiSymmetricImplR FermionImplPolicy;
typedef SpWilsonTwoIndexAntiSymmetricFermionR FermionAction; typedef SpWilsonTwoIndexAntiSymmetricFermionR FermionAction;
typedef typename FermionAction::FermionField FermionField; typedef typename FermionAction::FermionField FermionField;
@@ -36,24 +39,22 @@ int main(int argc, char **argv) {
typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs; typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
TheHMC.Resources.AddObservable<PlaqObs>(); TheHMC.Resources.AddObservable<PlaqObs>();
RealD beta = 6.7 ; RealD beta = 6.7;
SpWilsonGaugeActionR Waction(beta); SpWilsonGaugeActionR Waction(beta);
auto GridPtr = TheHMC.Resources.GetCartesian(); auto GridPtr = TheHMC.Resources.GetCartesian();
auto GridRBPtr = TheHMC.Resources.GetRBCartesian(); auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
SpTwoIndexAntiSymmetricRepresentation::LatticeField U(GridPtr); SpTwoIndexAntiSymmetricRepresentation::LatticeField U(GridPtr);
//LatticeGaugeField U(GridPtr); // LatticeGaugeField U(GridPtr);
RealD mass = -0.115; RealD mass = -0.115;
std::vector<Complex> boundary = {-1, -1, -1, -1};
std::vector<Complex> boundary = {-1,-1,-1,-1};
FermionAction::ImplParams bc(boundary); FermionAction::ImplParams bc(boundary);
FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, bc); FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, bc);
ConjugateGradient<FermionField> CG(1.0e-8, 2000, false); ConjugateGradient<FermionField> CG(1.0e-8, 2000, false);
TwoFlavourPseudoFermionAction<FermionImplPolicy> Nf2(FermOp, CG, CG); TwoFlavourPseudoFermionAction<FermionImplPolicy> Nf2(FermOp, CG, CG);
@@ -61,21 +62,20 @@ int main(int argc, char **argv) {
Nf2.is_smeared = false; Nf2.is_smeared = false;
std::cout << GridLogMessage << "mass " << mass << std::endl; std::cout << GridLogMessage << "mass " << mass << std::endl;
ActionLevel<HMCWrapper::Field, TheRepresentations > Level1(1); ActionLevel<HMCWrapper::Field, TheRepresentations> Level1(1);
Level1.push_back(&Nf2); Level1.push_back(&Nf2);
ActionLevel<HMCWrapper::Field, TheRepresentations > Level2(4); ActionLevel<HMCWrapper::Field, TheRepresentations> Level2(4);
Level2.push_back(&Waction); Level2.push_back(&Waction);
TheHMC.TheAction.push_back(Level1); TheHMC.TheAction.push_back(Level1);
TheHMC.TheAction.push_back(Level2); TheHMC.TheAction.push_back(Level2);
TheHMC.Parameters.MD.MDsteps = 16; TheHMC.Parameters.MD.MDsteps = 16;
TheHMC.Parameters.MD.trajL = 1.0; TheHMC.Parameters.MD.trajL = 1.0;
TheHMC.ReadCommandLine(argc, argv);
TheHMC.Run();
TheHMC.ReadCommandLine(argc, argv);
TheHMC.Run();
Grid_finalize(); Grid_finalize();
} }

39
tests/sp2n/Test_sp2n_lie_gen.cc
View File

@@ -5,14 +5,9 @@ using namespace Grid;
int main(int argc, char** argv) { int main(int argc, char** argv) {
Grid_init(&argc, &argv); Grid_init(&argc, &argv);
//std::vector<int> latt({4, 4, 4, 8});
//GridCartesian* grid = SpaceTimeGrid::makeFourDimGrid(
//latt, GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi());
//GridRedBlackCartesian* rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(grid);
std::cout << GridLogMessage << "*********************************************" std::cout << GridLogMessage << "*********************************************"
<< std::endl; << 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::cout << GridLogMessage << "*********************************************"
<< std::endl; << std::endl;
@@ -21,7 +16,7 @@ int main(int argc, char** argv) {
std::cout << GridLogMessage << "*********************************************" std::cout << GridLogMessage << "*********************************************"
<< std::endl; << 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::cout << GridLogMessage << "*********************************************"
<< std::endl; << std::endl;
@@ -30,21 +25,43 @@ int main(int argc, char** argv) {
std::cout << GridLogMessage << "*********************************************" std::cout << GridLogMessage << "*********************************************"
<< std::endl; << 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::cout << GridLogMessage << "*********************************************"
<< std::endl; << std::endl;
Sp6::printGenerators();
Sp6::testGenerators(); Sp6::testGenerators();
std::cout << GridLogMessage << "*********************************************" std::cout << GridLogMessage << "*********************************************"
<< std::endl; << 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::cout << GridLogMessage << "*********************************************"
<< std::endl; << std::endl;
Sp8::printGenerators();
Sp8::testGenerators(); Sp8::testGenerators();
std::cout << GridLogMessage << "*********************************************"
<< std::endl;
std::cout << GridLogMessage << "* Generators for Sp(4) TwoIndexAS (test)" << std::endl;
std::cout << GridLogMessage << "*********************************************"
<< std::endl;
Sp_TwoIndex<4, AntiSymmetric>::testGenerators();
std::cout << GridLogMessage << "*********************************************"
<< std::endl;
std::cout << GridLogMessage << "* Generators for Sp(6) TwoIndexAS (test)" << std::endl;
std::cout << GridLogMessage << "*********************************************"
<< std::endl;
Sp_TwoIndex<6, AntiSymmetric>::testGenerators();
std::cout << GridLogMessage << "*********************************************"
<< std::endl;
std::cout << GridLogMessage << "* Generators for Sp(8) TwoIndexAS (test)" << std::endl;
std::cout << GridLogMessage << "*********************************************"
<< std::endl;
Sp_TwoIndex<8, AntiSymmetric>::testGenerators();
Grid_finalize(); Grid_finalize();
} }