portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-09 23:45:36 +00:00
Code Releases Activity

2indx antisymm representation of sp2n

Browse Source
This commit is contained in:
Alessandro Lupo 2021-11-04 18:27:35 +00:00
parent 4044536eea
commit 88bdd4344b
34 changed files with 687 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
Grid/Makefile.am
View File

@ -66,6 +66,7 @@ if BUILD_FERMION_REPS
extra_sources+=$(SP_FERMION_FILES)
extra_sources+=$(ADJ_FERMION_FILES)
extra_sources+=$(TWOIND_FERMION_FILES)
extra_sources+=$(SP_TWOIND_FILES)
endif
lib_LIBRARIES = libGrid.a

13
Grid/qcd/action/fermion/Fermion.h
View File

@ -115,11 +115,22 @@ typedef WilsonFermion<WilsonImplR> WilsonFermionR;
typedef WilsonFermion<WilsonImplF> WilsonFermionF;
typedef WilsonFermion<WilsonImplD> WilsonFermionD;
// sp
typedef WilsonFermion<SpWilsonImplR> SpWilsonFermionR;
typedef WilsonFermion<SpWilsonImplF> SpWilsonFermionF;
typedef WilsonFermion<SpWilsonImplD> SpWilsonFermionD;
typedef WilsonFermion<SpWilsonTwoIndexAntiSymmetricImplR> SpWilsonTwoIndexAntiSymmetricFermionR;
typedef WilsonFermion<SpWilsonTwoIndexAntiSymmetricImplF> SpWilsonTwoIndexAntiSymmetricFermionF;
typedef WilsonFermion<SpWilsonTwoIndexAntiSymmetricImplD> SpWilsonTwoIndexAntiSymmetricFermionD;
typedef WilsonFermion<SpWilsonTwoIndexSymmetricImplR> SpWilsonTwoIndexSymmetricFermionR;
typedef WilsonFermion<SpWilsonTwoIndexSymmetricImplF> SpWilsonTwoIndexSymmetricFermionF;
typedef WilsonFermion<SpWilsonTwoIndexSymmetricImplD> SpWilsonTwoIndexSymmetricFermionD;
// end sp
//typedef WilsonFermion<WilsonImplRL> WilsonFermionRL;
//typedef WilsonFermion<WilsonImplFH> WilsonFermionFH;
//typedef WilsonFermion<WilsonImplDF> WilsonFermionDF;

10
Grid/qcd/action/fermion/WilsonImpl.h
View File

@ -273,9 +273,13 @@ typedef WilsonImpl<vComplex, SpFundamentalRepresentation, CoeffReal > SpWilsonI
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
typedef WilsonImpl<vComplex, SpTwoIndexAntiSymmetricRepresentation, CoeffReal > SpWilsonTwoIndexAntiSymmetricImplR; // Real.. whichever prec
typedef WilsonImpl<vComplexF, SpTwoIndexAntiSymmetricRepresentation, CoeffReal > SpWilsonTwoIndexAntiSymmetricImplF; // Float
typedef WilsonImpl<vComplexD, SpTwoIndexAntiSymmetricRepresentation, CoeffReal > SpWilsonTwoIndexAntiSymmetricImplD; // Double
typedef WilsonImpl<vComplex, SpTwoIndexSymmetricRepresentation, CoeffReal > SpWilsonTwoIndexSymmetricImplR; // Real.. whichever prec
typedef WilsonImpl<vComplexF, SpTwoIndexSymmetricRepresentation, CoeffReal > SpWilsonTwoIndexSymmetricImplF; // Float
typedef WilsonImpl<vComplexD, SpTwoIndexSymmetricRepresentation, CoeffReal > SpWilsonTwoIndexSymmetricImplD; // Double
NAMESPACE_END(Grid);

1
Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexAntiSymmetricImplD/WilsonCloverFermionInstantiationSpWilsonTwoIndexAntiSymmetricImplD.cc Symbolic link
View File

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

1
Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexAntiSymmetricImplD/WilsonFermionInstantiationSpWilsonTwoIndexAntiSymmetricImplD.cc Symbolic link
View File

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

1
Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexAntiSymmetricImplD/WilsonKernelsInstantiationSpWilsonTwoIndexAntiSymmetricImplD.cc Symbolic link
View File

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

1
Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexAntiSymmetricImplD/WilsonTMFermionInstantiationSpWilsonTwoIndexAntiSymmetricImplD.cc Symbolic link
View File

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

1
Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexAntiSymmetricImplD/impl.h Normal file
View File

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

1
Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexAntiSymmetricImplF/WilsonCloverFermionInstantiationSpWilsonTwoIndexAntiSymmetricImplF.cc Symbolic link
View File

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

1
Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexAntiSymmetricImplF/WilsonFermionInstantiationSpWilsonTwoIndexAntiSymmetricImplF.cc Symbolic link
View File

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

1
Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexAntiSymmetricImplF/WilsonKernelsInstantiationSpWilsonTwoIndexAntiSymmetricImplF.cc Symbolic link
View File

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

1
Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexAntiSymmetricImplF/WilsonTMFermionInstantiationSpWilsonTwoIndexAntiSymmetricImplF.cc Symbolic link
View File

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

1
Grid/qcd/action/fermion/instantiation/SpWilsonTwoIndexAntiSymmetricImplF/impl.h Normal file
View File

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4
Grid/qcd/action/fermion/instantiation/generate_instantiations.sh
View File

@ -17,6 +17,10 @@ WILSON_IMPL_LIST=" \
WilsonTwoIndexSymmetricImplD \
WilsonTwoIndexAntiSymmetricImplF \
WilsonTwoIndexAntiSymmetricImplD \
SpWilsonTwoIndexAntiSymmetricImplF \
SpWilsonTwoIndexAntiSymmetricImplD \
SpWilsonTwoIndexSymmetricImplF \
SpWilsonTwoIndexSymmetricImplD \
GparityWilsonImplF \
GparityWilsonImplD "

1
Grid/qcd/representations/Representations.h
View File

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

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

@ -4,6 +4,8 @@
#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/sp_two_index.h>
#include <Grid/qcd/action/scalar/ScalarImpl.h>
#include <tuple>

100
Grid/qcd/representations/sp_two_index.h Normal file
View File

@ -0,0 +1,100 @@
/*
* Policy classes for the HMC
* Authors: Guido Cossu, David Preti
*/
#ifndef SP2N2INDEX_H_H
#define SP2N2INDEX_H_H
NAMESPACE_BEGIN(Grid);
/*
* This is an helper class for the HMC
* Should contain only the data for the two index representations
* 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 SpTwoIndexRep {
public:
// typdef to be used by the Representations class in HMC to get the
// 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 bool isFundamental = false;
static const int nsp = Nc / 2;
LatticeField U;
explicit SpTwoIndexRep(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 TwoIndexRep
// (U)_{(ij)(lk)} = tr [ adj(e^(ij)) U e^(lk) transpose(U) ]
conformable(U, Uin);
U = Zero();
LatticeColourMatrix tmp(Uin.Grid());
Vector<typename Sp<nsp>::Matrix> eij(Dimension);
for (int a = 0; a < Dimension; a++)
Sp_TwoIndex<ncolour, S>::base(a, eij[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);
}
}
LatticeGaugeField RtoFundamentalProject(const LatticeField &in,
Real scale = 1.0) const {
LatticeGaugeField out(in.Grid());
out = Zero();
for (int mu = 0; mu < Nd; mu++) {
LatticeColourMatrix out_mu(in.Grid()); // fundamental representation
LatticeMatrix in_mu = peekLorentz(in, mu);
out_mu = Zero();
typename Sp<nsp>::LatticeAlgebraVector h(in.Grid());
projectOnAlgebra(h, in_mu, double(Nc + 2 * S)); // factor T(r)/T(fund)
FundamentalLieAlgebraMatrix(h, out_mu); // apply scale only once
pokeLorentz(out, out_mu, mu); // should be 2 for sp4 as. ok
}
return out;
}
private:
void projectOnAlgebra(typename Sp<nsp>::LatticeAlgebraVector &h_out,
const LatticeMatrix &in, Real scale = 1.0) const {
Sp_TwoIndex<ncolour, S>::projectOnAlgebra(h_out, in, scale);
}
void FundamentalLieAlgebraMatrix(
typename Sp<nsp>::LatticeAlgebraVector &h,
typename Sp<nsp>::LatticeMatrix &out, Real scale = 1.0) const {
Sp<nsp>::FundamentalLieAlgebraMatrix(h, out, scale);
}
};
typedef SpTwoIndexRep<Nc, Symmetric> SpTwoIndexSymmetricRepresentation;
typedef SpTwoIndexRep<Nc, AntiSymmetric> SpTwoIndexAntiSymmetricRepresentation;
NAMESPACE_END(Grid);
#endif

27
Grid/qcd/utils/Sp2n.h
View File

@ -19,6 +19,7 @@ public:
static const int Dimension = ncolour*2;
static const int AlgebraDimension = ncolour*(2*ncolour +1);
static int su2subgroups(void) { return (ncolour * (ncolour - 1)) / 2; }
static const int nnsp = ncolour;
template <typename vtype>
@ -368,6 +369,21 @@ public:
out *= ci;
}
static void FundamentalLieAlgebraMatrix(const LatticeAlgebraVector &h,
LatticeMatrix &out,
Real scale = 1.0) {
conformable(h, out);
GridBase *grid = out.Grid();
LatticeMatrix la(grid);
Matrix ta;
out = Zero();
for (int a = 0; a < AlgebraDimension; a++) {
generator(a, ta);
la = peekColour(h, a) * timesI(ta) * scale;
out += la;
}
}
template <typename LatticeMatrixType>
@ -388,7 +404,16 @@ public:
}
}
static void projectOnAlgebra(LatticeAlgebraVector &h_out, const LatticeMatrix &in, Real scale = 1.0) {
conformable(h_out, in);
h_out = Zero();
Matrix Ta;
for (int a = 0; a < AlgebraDimension; a++) {
generator(a, Ta);
pokeColour(h_out, - 2.0 * (trace(timesI(Ta) * in)) * scale, a);
}
}
template <typename GaugeField>

320
Grid/qcd/utils/Sp2nTwoIndex.h Normal file
View File

@ -0,0 +1,320 @@
////////////////////////////////////////////////////////////////////////
//
// * 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_SP2N2INDEX_H
#define QCD_UTIL_SP2N2INDEX_H
NAMESPACE_BEGIN(Grid);
//enum SpTwoIndexSymmetry { S = 1, AS = -1 };
//inline Real delta(int a, int b) { return (a == b) ? 1.0 : 0.0; }
template <int ncolour, TwoIndexSymmetry S>
class Sp_TwoIndex : public Sp<Config_nsp> {
public:
static const int nsp = ncolour/2;
static const int Dimension = nsp * (ncolour + S) - 1 ;
static const int NumGenerators = Sp<nsp>::AlgebraDimension;
template <typename vtype>
using iSp2nTwoIndexMatrix = iScalar<iScalar<iMatrix<vtype, Dimension> > >;
typedef iSp2nTwoIndexMatrix<Complex> TIMatrix;
typedef iSp2nTwoIndexMatrix<ComplexF> TIMatrixF;
typedef iSp2nTwoIndexMatrix<ComplexD> TIMatrixD;
typedef iSp2nTwoIndexMatrix<vComplex> vTIMatrix;
typedef iSp2nTwoIndexMatrix<vComplexF> vTIMatrixF;
typedef iSp2nTwoIndexMatrix<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 iSp2nMatrix = iScalar<iScalar<iMatrix<vtype, ncolour> > >;
typedef iSp2nMatrix<Complex> Matrix;
typedef iSp2nMatrix<ComplexF> MatrixF;
typedef iSp2nMatrix<ComplexD> MatrixD;
template <class cplx>
static void base(int Index, iSp2nMatrix<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++) {
if (i==nsp)
{
j = j+1;
}
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, iSp2nMatrix<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, iSp2nMatrix<cplx> &eij) {
eij = Zero();
//std::cout << GridLogMessage << " doing i j = " << i << j << std::endl;
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 ;
}
}
}
else if (i != nsp+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;
}
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, iSp2nTwoIndexMatrix<cplx> &i2indTa) {
Vector<typename Sp<nsp>::template iSp2nMatrix<cplx> > ta(
NumGenerators);
Vector<typename Sp<nsp>::template iSp2nMatrix<cplx> > eij(Dimension);
typename Sp<nsp>::template iSp2nMatrix<cplx> tmp;
i2indTa = Zero();
for (int a = 0; a < NumGenerators; a++)
Sp<nsp>::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 Sp<nsp>::template iSp2nMatrix<cplx> tmp1 =
tmp * eij[b];
Complex iTr = TensorRemove(timesI(trace(tmp1)));
i2indTa()()(a, b) = iTr;
}
}
}
static void printGenerators(void) {
for (int gen = 0; gen < NumGenerators; 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 < NumGenerators; 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 < NumGenerators; 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 < NumGenerators; a++) {
for (int b = 0; b < NumGenerators; 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 Sp<nsp>::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 < NumGenerators; 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 Sp<nsp>::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 < NumGenerators; a++) {
generator(a, i2indTa);
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
static void projector(typename Sp<nsp>::LatticeAlgebraVector &h_out,
const LatticeTwoIndexMatrix &in, Real scale = 1.0) {
conformable(h_out, in);
// to store the generators
static std::vector<TIMatrix> i2indTa(NumGenerators);
h_out = Zero();
static bool precalculated = false;
if (!precalculated) {
precalculated = true;
for (int a = 0; a < NumGenerators; 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 < NumGenerators; a++) {
auto tmp = real(trace(i2indTa[a] * in)) * coefficient;
pokeColour(h_out, tmp, a);
}
}
};
// Some useful type names
typedef Sp_TwoIndex<Nc, Symmetric> SpTwoIndexSymmMatrices;
typedef Sp_TwoIndex<Nc, AntiSymmetric> SpTwoIndexAntiSymmMatrices;
typedef Sp_TwoIndex<2, Symmetric> Sp2TwoIndexSymm;
typedef Sp_TwoIndex<4, Symmetric> Sp4TwoIndexSymm;
typedef Sp_TwoIndex<2, AntiSymmetric> Sp2TwoIndexAntiSymm;
typedef Sp_TwoIndex<4, AntiSymmetric> Sp4TwoIndexAntiSymm;
NAMESPACE_END(Grid);
#endif

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

@ -12,6 +12,7 @@
#include <Grid/qcd/utils/Sp2n.h>
#include <Grid/qcd/utils/SUnAdjoint.h>
#include <Grid/qcd/utils/SUnTwoIndex.h>
#include <Grid/qcd/utils/Sp2nTwoIndex.h>
// All-to-all contraction kernels that touch the
// internal lattice structure

9
configure.ac
View File

@ -181,13 +181,18 @@ AC_ARG_ENABLE([Nc],
case ${ac_Nc} in
2)
AC_DEFINE([Config_Nc],[2],[Gauge group Nc]);;
AC_DEFINE([Config_Nc],[2],[Gauge group Nc])
AC_DEFINE([Config_nsp],[1],[n of Sp2n]);;
3)
AC_DEFINE([Config_Nc],[3],[Gauge group Nc]);;
4)
AC_DEFINE([Config_Nc],[4],[Gauge group Nc]);;
AC_DEFINE([Config_Nc],[4],[Gauge group Nc])
AC_DEFINE([Config_nsp],[2],[n of Sp2n]);;
5)
AC_DEFINE([Config_Nc],[5],[Gauge group Nc]);;
8)
AC_DEFINE([Config_Nc],[8],[Gauge group Nc])
AC_DEFINE([Config_nsp],[4],[n of Sp2n]);;
*)
AC_MSG_ERROR(["Unsupport gauge group choice Nc = ${ac_Nc}"]);;
esac

2
scripts/filelist
View File

@ -16,6 +16,7 @@ GP_FERMION_FILES=` find . -name '*.cc' -path '*/instantiation/*' -path '*/ins
ADJ_FERMION_FILES=` find . -name '*.cc' -path '*/instantiation/*' -path '*/instantiation/WilsonAdj*' `
TWOIND_FERMION_FILES=`find . -name '*.cc' -path '*/instantiation/*' -path '*/instantiation/WilsonTwoIndex*'`
SP_FERMION_FILES=`find . -name '*.cc' -path '*/instantiation/*' -path '*/instantiation/SpWilsonImpl*'`
SP_TWOIND_FILES=`find . -name '*.cc' -path '*/instantiation/*' -path '*/instantiation/SpWilsonTwo*'`
HPPFILES=`find . -type f -name '*.hpp'`
echo HFILES=$HFILES $HPPFILES > Make.inc
@ -29,6 +30,7 @@ echo GP_FERMION_FILES=$GP_FERMION_FILES >> Make.inc
echo ADJ_FERMION_FILES=$ADJ_FERMION_FILES >> Make.inc
echo TWOIND_FERMION_FILES=$TWOIND_FERMION_FILES >> Make.inc
echo SP_FERMION_FILES=$SP_FERMION_FILES >> Make.inc
echo SP_TWOIND_FILES=$SP_TWOIND_FILES >> Make.inc
# tests Make.inc
cd $home/tests

102
tests/sp2n/Test_hmc_Sp_WF_2_Fund_3_2AS.cc Normal file
View File

@ -0,0 +1,102 @@
#include <Grid/Grid.h>
int main(int argc, char **argv) {
using namespace Grid;
typedef Representations< SpFundamentalRepresentation, SpTwoIndexAntiSymmetricRepresentation > TheRepresentations;
Grid_init(&argc, &argv);
typedef GenericSp2nHMCRunnerHirep<TheRepresentations, MinimumNorm2> HMCWrapper;
typedef SpWilsonTwoIndexAntiSymmetricImplR TwoIndexFermionImplPolicy;
typedef SpWilsonTwoIndexAntiSymmetricFermionR TwoIndexFermionAction;
typedef typename TwoIndexFermionAction::FermionField TwoIndexFermionField;
typedef SpWilsonImplR FundFermionImplPolicy; // ok
typedef SpWilsonFermionR FundFermionAction; // ok
typedef typename FundFermionAction::FermionField FundFermionField;
//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
HMCWrapper TheHMC;
TheHMC.Resources.AddFourDimGrid("gauge");
// Checkpointer definition
CheckpointerParameters CPparams;
CPparams.config_prefix = "ckpoint_lat";
CPparams.rng_prefix = "ckpoint_rng";
CPparams.saveInterval = 5;
CPparams.format = "IEEE64BIG";
TheHMC.Resources.LoadNerscCheckpointer(CPparams);
RNGModuleParameters RNGpar;
RNGpar.serial_seeds = "1 2 3 4 5";
RNGpar.parallel_seeds = "6 7 8 9 10";
TheHMC.Resources.SetRNGSeeds(RNGpar);
// Construct observables
typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
TheHMC.Resources.AddObservable<PlaqObs>();
typedef PolyakovMod<HMCWrapper::ImplPolicy> PolyakovObs;
TheHMC.Resources.AddObservable<PolyakovObs>();
RealD beta = 6 ;
SymplWilsonGaugeActionR Waction(beta);
auto GridPtr = TheHMC.Resources.GetCartesian();
auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
SpFundamentalRepresentation::LatticeField fundU(GridPtr);
SpTwoIndexAntiSymmetricRepresentation::LatticeField asU(GridPtr);
//LatticeGaugeField U(GridPtr);
RealD Fundmass = -0.71;
RealD ASmass = -0.71;
std::vector<Complex> boundary = {-1,-1,-1,-1};
FundFermionAction::ImplParams bc(boundary);
TwoIndexFermionAction::ImplParams bbc(boundary);
FundFermionAction FundFermOp(fundU, *GridPtr, *GridRBPtr, Fundmass, bbc);
TwoIndexFermionAction TwoIndexFermOp(asU, *GridPtr, *GridRBPtr, ASmass, bbc);
ConjugateGradient<FundFermionField> fCG(1.0e-8, 2000, false);
ConjugateGradient<TwoIndexFermionField> asCG(1.0e-8, 2000, false);
OneFlavourRationalParams Params(1.0e-6, 64.0, 2000, 1.0e-6, 16);
TwoFlavourPseudoFermionAction<FundFermionImplPolicy> fundNf2(FundFermOp, fCG, fCG);
TwoFlavourPseudoFermionAction<TwoIndexFermionImplPolicy> asNf2(TwoIndexFermOp, asCG, asCG);
OneFlavourRationalPseudoFermionAction<TwoIndexFermionImplPolicy> asNf1(TwoIndexFermOp,Params);
fundNf2.is_smeared = false;
asNf2.is_smeared = false;
asNf1.is_smeared = false;
ActionLevel<HMCWrapper::Field, TheRepresentations > Level1(1);
Level1.push_back(&fundNf2);
Level1.push_back(&asNf2);
Level1.push_back(&asNf1);
ActionLevel<HMCWrapper::Field, TheRepresentations > Level2(4);
Level2.push_back(&Waction);
TheHMC.TheAction.push_back(Level1);
TheHMC.TheAction.push_back(Level2);
TheHMC.Parameters.MD.MDsteps = 28;
TheHMC.Parameters.MD.trajL = 1.0;
TheHMC.ReadCommandLine(argc, argv);
TheHMC.Run();
Grid_finalize();
}

82
tests/sp2n/Test_hmc_Sp_Wilson2ASFermionGauge.cc Normal file
View File

@ -0,0 +1,82 @@
#include <Grid/Grid.h>
int main(int argc, char **argv) {
using namespace Grid;
typedef Representations< SpFundamentalRepresentation, SpTwoIndexAntiSymmetricRepresentation > TheRepresentations;
Grid_init(&argc, &argv);
typedef GenericSp2nHMCRunnerHirep<TheRepresentations, MinimumNorm2> HMCWrapper;
typedef SpWilsonTwoIndexAntiSymmetricImplR FermionImplPolicy;
typedef SpWilsonTwoIndexAntiSymmetricFermionR FermionAction;
typedef typename FermionAction::FermionField FermionField;
//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
HMCWrapper TheHMC;
TheHMC.Resources.AddFourDimGrid("gauge");
// Checkpointer definition
CheckpointerParameters CPparams;
CPparams.config_prefix = "ckpoint_lat";
CPparams.rng_prefix = "ckpoint_rng";
CPparams.saveInterval = 5;
CPparams.format = "IEEE64BIG";
TheHMC.Resources.LoadNerscCheckpointer(CPparams);
RNGModuleParameters RNGpar;
RNGpar.serial_seeds = "1 2 3 4 5";
RNGpar.parallel_seeds = "6 7 8 9 10";
TheHMC.Resources.SetRNGSeeds(RNGpar);
// Construct observables
typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
TheHMC.Resources.AddObservable<PlaqObs>();
RealD beta = 6.7 ;
SymplWilsonGaugeActionR Waction(beta);
auto GridPtr = TheHMC.Resources.GetCartesian();
auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
SpTwoIndexAntiSymmetricRepresentation::LatticeField U(GridPtr);
//LatticeGaugeField U(GridPtr);
RealD mass = -0.85;
std::vector<Complex> boundary = {-1,-1,-1,-1};
FermionAction::ImplParams bc(boundary);
FermionAction FermOp(U, *GridPtr, *GridRBPtr, mass, bc);
ConjugateGradient<FermionField> CG(1.0e-8, 2000, false);
TwoFlavourPseudoFermionAction<FermionImplPolicy> Nf2(FermOp, CG, CG);
Nf2.is_smeared = false;
std::cout << GridLogMessage << "mass " << mass << std::endl;
ActionLevel<HMCWrapper::Field, TheRepresentations > Level1(1);
Level1.push_back(&Nf2);
ActionLevel<HMCWrapper::Field, TheRepresentations > Level2(4);
Level2.push_back(&Waction);
TheHMC.TheAction.push_back(Level1);
TheHMC.TheAction.push_back(Level2);
TheHMC.Parameters.MD.MDsteps = 16;
TheHMC.Parameters.MD.trajL = 1.0;
TheHMC.ReadCommandLine(argc, argv);
TheHMC.Run();
Grid_finalize();
}