From 35e10a1159ce6f8c05217d682ce5f85db1c9e2c2 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Fri, 3 Oct 2025 12:17:13 -0400 Subject: [PATCH] Changes for Nd=3 --- Grid/GridQCDcore.h | 1 + Grid/cartesian/Cartesian_base.h | 1 + Grid/qcd/QCD.h | 3 +- Grid/qcd/action/fermion/DomainWallVec5dImpl.h | 4 +- Grid/qcd/action/fermion/Fermion.h | 9 ++++ Grid/qcd/action/fermion/FermionCore.h | 3 +- Grid/qcd/action/fermion/FermionOperatorImpl.h | 6 +++ Grid/qcd/action/fermion/GparityWilsonImpl.h | 6 +-- .../action/fermion/ImprovedStaggeredFermion.h | 2 + .../fermion/ImprovedStaggeredFermion5D.h | 2 + .../action/fermion/NaiveStaggeredFermion.h | 2 + Grid/qcd/action/fermion/StaggeredImpl.h | 6 +-- Grid/qcd/action/fermion/WilsonFermion.h | 2 + Grid/qcd/action/fermion/WilsonFermion5D.h | 2 + Grid/qcd/action/fermion/WilsonImpl.h | 2 +- Grid/qcd/action/fermion/WilsonTMFermion5D.h | 12 ++--- .../DomainWallEOFAFermionImplementation.h | 2 +- .../WilsonCloverFermionImplementation.h | 14 +++--- .../WilsonFermion5DImplementation.h | 36 ++++++++------- .../WilsonFermionImplementation.h | 10 ++--- .../WilsonKernelsHandGparityImplementation.h | 2 +- ...ImprovedStaggeredFermion5DInstantiation.cc | 26 ++++++++++- .../ImprovedStaggeredFermionInstantiation.cc | 23 +++++++++- .../NaiveStaggeredFermionInstantiation.cc | 24 +++++++++- .../WilsonFermion5DInstantiation.cc | 24 +++++++++- .../WilsonFermionInstantiation.cc | 22 ++++++++- .../instantiation/generate_instantiations.sh | 18 ++++++-- Grid/qcd/smearing/WilsonFlow.h | 4 +- Grid/qcd/utils/WilsonLoops.h | 45 +++++++++++++------ Grid/util/Init.cc | 4 +- configure.ac | 17 +++++++ 31 files changed, 256 insertions(+), 78 deletions(-) diff --git a/Grid/GridQCDcore.h b/Grid/GridQCDcore.h index 065b62cd..db63f565 100644 --- a/Grid/GridQCDcore.h +++ b/Grid/GridQCDcore.h @@ -37,6 +37,7 @@ Author: paboyle #include #include #include +#include // depends on Gparity #include #include NAMESPACE_CHECK(GridQCDCore); diff --git a/Grid/cartesian/Cartesian_base.h b/Grid/cartesian/Cartesian_base.h index 66400787..d029b9a7 100644 --- a/Grid/cartesian/Cartesian_base.h +++ b/Grid/cartesian/Cartesian_base.h @@ -90,6 +90,7 @@ public: // Checkerboarding interface is virtual and overridden by // GridCartesian / GridRedBlackCartesian //////////////////////////////////////////////////////////////// + virtual int Icosahedral(void) { return 0;} virtual int CheckerBoarded(int dim) =0; virtual int CheckerBoard(const Coordinate &site)=0; virtual int CheckerBoardDestination(int source_cb,int shift,int dim)=0; diff --git a/Grid/qcd/QCD.h b/Grid/qcd/QCD.h index dbedfa7c..3325c4b1 100644 --- a/Grid/qcd/QCD.h +++ b/Grid/qcd/QCD.h @@ -49,7 +49,7 @@ static constexpr int Tm = 7; static constexpr int Nc=Config_Nc; static constexpr int Ns=4; -static constexpr int Nd=4; +static constexpr int Nd=Config_Nd; static constexpr int Nhs=2; // half spinor static constexpr int Nds=8; // double stored gauge field static constexpr int Ngp=2; // gparity index range @@ -75,6 +75,7 @@ static constexpr int InverseYes=1; //typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; const int SpinorIndex = 2; +const int PauliIndex = 2; //TensorLevel counts from the bottom! template struct isSpinor { static constexpr bool value = (SpinorIndex==T::TensorLevel); }; diff --git a/Grid/qcd/action/fermion/DomainWallVec5dImpl.h b/Grid/qcd/action/fermion/DomainWallVec5dImpl.h index 0c8a0930..3135b7b4 100644 --- a/Grid/qcd/action/fermion/DomainWallVec5dImpl.h +++ b/Grid/qcd/action/fermion/DomainWallVec5dImpl.h @@ -123,10 +123,10 @@ public: GaugeGrid->LocalIndexToLocalCoor(lidx, lcoor); peekLocalSite(ScalarUmu, Umu_v, lcoor); - for (int mu = 0; mu < 4; mu++) ScalarUds(mu) = ScalarUmu(mu); + for (int mu = 0; mu < Nd; mu++) ScalarUds(mu) = ScalarUmu(mu); peekLocalSite(ScalarUmu, Uadj_v, lcoor); - for (int mu = 0; mu < 4; mu++) ScalarUds(mu + 4) = ScalarUmu(mu); + for (int mu = 0; mu < Nd; mu++) ScalarUds(mu + Nd) = ScalarUmu(mu); pokeLocalSite(ScalarUds, Uds_v, lcoor); }); diff --git a/Grid/qcd/action/fermion/Fermion.h b/Grid/qcd/action/fermion/Fermion.h index a3d96d9b..28165915 100644 --- a/Grid/qcd/action/fermion/Fermion.h +++ b/Grid/qcd/action/fermion/Fermion.h @@ -85,6 +85,15 @@ NAMESPACE_CHECK(DomainWall); #include #include NAMESPACE_CHECK(Overlap); + + +/////////////////////////////////////////////////////////////////////////////// +// Two spin wilson fermion based +/////////////////////////////////////////////////////////////////////////////// + +#include +NAMESPACE_CHECK(TwoSpinWilson); + /////////////////////////////////////////////////////////////////////////////// // G5 herm -- this has to live in QCD since dirac matrix is not in the broader sector of code /////////////////////////////////////////////////////////////////////////////// diff --git a/Grid/qcd/action/fermion/FermionCore.h b/Grid/qcd/action/fermion/FermionCore.h index 6745032e..fa6ec072 100644 --- a/Grid/qcd/action/fermion/FermionCore.h +++ b/Grid/qcd/action/fermion/FermionCore.h @@ -41,8 +41,9 @@ NAMESPACE_CHECK(Compressor); NAMESPACE_CHECK(FermionOperatorImpl); #include NAMESPACE_CHECK(FermionOperator); -#include //used by all wilson type fermions +#include //used by all wilson type fermions #include //used by all wilson type fermions +#include //used for 3D fermions, pauli in place of Dirac NAMESPACE_CHECK(Kernels); #endif diff --git a/Grid/qcd/action/fermion/FermionOperatorImpl.h b/Grid/qcd/action/fermion/FermionOperatorImpl.h index 56aaca12..6593fda9 100644 --- a/Grid/qcd/action/fermion/FermionOperatorImpl.h +++ b/Grid/qcd/action/fermion/FermionOperatorImpl.h @@ -180,6 +180,12 @@ NAMESPACE_CHECK(ImplGparityWilson); #include NAMESPACE_CHECK(ImplStaggered); +///////////////////////////////////////////////////////////////////////////// +// Two component spinor Wilson action for 3d / Boston +///////////////////////////////////////////////////////////////////////////// +#include +NAMESPACE_CHECK(ImplTwoSpinWilson); + ///////////////////////////////////////////////////////////////////////////// // Single flavour one component spinors with colour index. 5d vec ///////////////////////////////////////////////////////////////////////////// diff --git a/Grid/qcd/action/fermion/GparityWilsonImpl.h b/Grid/qcd/action/fermion/GparityWilsonImpl.h index 8017bc76..1ab4e4f2 100644 --- a/Grid/qcd/action/fermion/GparityWilsonImpl.h +++ b/Grid/qcd/action/fermion/GparityWilsonImpl.h @@ -274,7 +274,7 @@ public: autoView( Uds_v , Uds, CpuWrite); autoView( Utmp_v, Utmp, CpuWrite); thread_foreach(ss,Utmp_v,{ - Uds_v[ss](0)(mu+4) = Utmp_v[ss](); + Uds_v[ss](0)(mu+Nd) = Utmp_v[ss](); }); } Utmp = Uconj; @@ -286,7 +286,7 @@ public: autoView( Uds_v , Uds, CpuWrite); autoView( Utmp_v, Utmp, CpuWrite); thread_foreach(ss,Utmp_v,{ - Uds_v[ss](1)(mu+4) = Utmp_v[ss](); + Uds_v[ss](1)(mu+Nd) = Utmp_v[ss](); }); } } @@ -320,7 +320,7 @@ public: } Uconj = conjugate(*Upoke); - pokeGparityDoubledGaugeField(Uds, *Upoke, Uconj, mu + 4); + pokeGparityDoubledGaugeField(Uds, *Upoke, Uconj, mu + Nd); } } diff --git a/Grid/qcd/action/fermion/ImprovedStaggeredFermion.h b/Grid/qcd/action/fermion/ImprovedStaggeredFermion.h index f7655f24..1c05d398 100644 --- a/Grid/qcd/action/fermion/ImprovedStaggeredFermion.h +++ b/Grid/qcd/action/fermion/ImprovedStaggeredFermion.h @@ -36,6 +36,8 @@ public: static const std::vector directions; static const std::vector displacements; static const int npoint = 16; + static std::vector MakeDirections(void); + static std::vector MakeDisplacements(void); }; template diff --git a/Grid/qcd/action/fermion/ImprovedStaggeredFermion5D.h b/Grid/qcd/action/fermion/ImprovedStaggeredFermion5D.h index 2641a6b8..192645c6 100644 --- a/Grid/qcd/action/fermion/ImprovedStaggeredFermion5D.h +++ b/Grid/qcd/action/fermion/ImprovedStaggeredFermion5D.h @@ -40,6 +40,8 @@ public: static const std::vector directions; static const std::vector displacements; const int npoint = 16; + static std::vector MakeDirections(void); + static std::vector MakeDisplacements(void); }; template diff --git a/Grid/qcd/action/fermion/NaiveStaggeredFermion.h b/Grid/qcd/action/fermion/NaiveStaggeredFermion.h index 9ec6be90..551beac5 100644 --- a/Grid/qcd/action/fermion/NaiveStaggeredFermion.h +++ b/Grid/qcd/action/fermion/NaiveStaggeredFermion.h @@ -36,6 +36,8 @@ public: static const std::vector directions; static const std::vector displacements; static const int npoint = 8; + static std::vector MakeDirections(void); + static std::vector MakeDisplacements(void); }; template diff --git a/Grid/qcd/action/fermion/StaggeredImpl.h b/Grid/qcd/action/fermion/StaggeredImpl.h index f44d12f4..b5aa03f0 100644 --- a/Grid/qcd/action/fermion/StaggeredImpl.h +++ b/Grid/qcd/action/fermion/StaggeredImpl.h @@ -141,9 +141,9 @@ public: Udag = Udag *phases; InsertGaugeField(Uds,U,mu); - InsertGaugeField(Uds,Udag,mu+4); + InsertGaugeField(Uds,Udag,mu+Nd); // PokeIndex(Uds, U, mu); - // PokeIndex(Uds, Udag, mu + 4); + // PokeIndex(Uds, Udag, mu + Nd); // 3 hop based on thin links. Crazy huh ? U = PeekIndex(Uthin, mu); @@ -156,7 +156,7 @@ public: UUUdag = UUUdag *phases; InsertGaugeField(UUUds,UUU,mu); - InsertGaugeField(UUUds,UUUdag,mu+4); + InsertGaugeField(UUUds,UUUdag,mu+Nd); } } diff --git a/Grid/qcd/action/fermion/WilsonFermion.h b/Grid/qcd/action/fermion/WilsonFermion.h index 16320a93..d64d2f15 100644 --- a/Grid/qcd/action/fermion/WilsonFermion.h +++ b/Grid/qcd/action/fermion/WilsonFermion.h @@ -38,6 +38,8 @@ public: static int MortonOrder; static const std::vector directions; static const std::vector displacements; + static std::vector MakeDirections(void); + static std::vector MakeDisplacements(void); static const int npoint = 8; }; diff --git a/Grid/qcd/action/fermion/WilsonFermion5D.h b/Grid/qcd/action/fermion/WilsonFermion5D.h index 9dadc866..4c2bcf9c 100644 --- a/Grid/qcd/action/fermion/WilsonFermion5D.h +++ b/Grid/qcd/action/fermion/WilsonFermion5D.h @@ -62,6 +62,8 @@ public: static const std::vector directions; static const std::vector displacements; static constexpr int npoint = 8; + static std::vector MakeDirections(void); + static std::vector MakeDisplacements(void); }; template diff --git a/Grid/qcd/action/fermion/WilsonImpl.h b/Grid/qcd/action/fermion/WilsonImpl.h index 07248160..bee2b8f8 100644 --- a/Grid/qcd/action/fermion/WilsonImpl.h +++ b/Grid/qcd/action/fermion/WilsonImpl.h @@ -166,7 +166,7 @@ public: U = adj(Cshift(U, mu, -1)); U = where(coor == 0, conjugate(phase) * U, U); - PokeIndex(Uds, U, mu + 4); + PokeIndex(Uds, U, mu + Nd); } } diff --git a/Grid/qcd/action/fermion/WilsonTMFermion5D.h b/Grid/qcd/action/fermion/WilsonTMFermion5D.h index 982e722a..a5fb4a20 100644 --- a/Grid/qcd/action/fermion/WilsonTMFermion5D.h +++ b/Grid/qcd/action/fermion/WilsonTMFermion5D.h @@ -56,7 +56,7 @@ class WilsonTMFermion5D : public WilsonFermion5D Frbgrid, Ugrid, Urbgrid, - 4.0,p) + Nd*1.0,p) { update(_mass,_mu); @@ -83,7 +83,7 @@ class WilsonTMFermion5D : public WilsonFermion5D out.Checkerboard() = in.Checkerboard(); //axpibg5x(out,in,a,b); // out = a*in + b*i*G5*in for (int s=0;s<(int)this->mass.size();s++) { - ComplexD a = 4.0+this->mass[s]; + ComplexD a = Nd*1.0+this->mass[s]; ComplexD b(0.0,this->mu[s]); axpbg5y_ssp(out,a,in,b,in,s,s); } @@ -92,7 +92,7 @@ class WilsonTMFermion5D : public WilsonFermion5D virtual void MooeeDag(const FermionField &in, FermionField &out) { out.Checkerboard() = in.Checkerboard(); for (int s=0;s<(int)this->mass.size();s++) { - ComplexD a = 4.0+this->mass[s]; + ComplexD a = Nd*1.0+this->mass[s]; ComplexD b(0.0,-this->mu[s]); axpbg5y_ssp(out,a,in,b,in,s,s); } @@ -101,7 +101,7 @@ class WilsonTMFermion5D : public WilsonFermion5D for (int s=0;s<(int)this->mass.size();s++) { RealD m = this->mass[s]; RealD tm = this->mu[s]; - RealD mtil = 4.0+this->mass[s]; + RealD mtil = Nd*1.0+this->mass[s]; RealD sq = mtil*mtil+tm*tm; ComplexD a = mtil/sq; ComplexD b(0.0, -tm /sq); @@ -112,7 +112,7 @@ class WilsonTMFermion5D : public WilsonFermion5D for (int s=0;s<(int)this->mass.size();s++) { RealD m = this->mass[s]; RealD tm = this->mu[s]; - RealD mtil = 4.0+this->mass[s]; + RealD mtil = Nd*1.0+this->mass[s]; RealD sq = mtil*mtil+tm*tm; ComplexD a = mtil/sq; ComplexD b(0.0,tm /sq); @@ -126,7 +126,7 @@ class WilsonTMFermion5D : public WilsonFermion5D this->Dhop(in, out, DaggerNo); FermionField tmp(out.Grid()); for (int s=0;s<(int)this->mass.size();s++) { - ComplexD a = 4.0+this->mass[s]; + ComplexD a = Nd*1.0+this->mass[s]; ComplexD b(0.0,this->mu[s]); axpbg5y_ssp(tmp,a,in,b,in,s,s); } diff --git a/Grid/qcd/action/fermion/implementation/DomainWallEOFAFermionImplementation.h b/Grid/qcd/action/fermion/implementation/DomainWallEOFAFermionImplementation.h index 53b44ca2..47074fd4 100644 --- a/Grid/qcd/action/fermion/implementation/DomainWallEOFAFermionImplementation.h +++ b/Grid/qcd/action/fermion/implementation/DomainWallEOFAFermionImplementation.h @@ -240,7 +240,7 @@ void DomainWallEOFAFermion::SetCoefficientsInternal(RealD zolo_hi, std::ve this->ceo.resize(Ls); for(int i=0; ibee[i] = 4.0 - this->M5 + 1.0; + this->bee[i] = Nd*1.0 - this->M5 + 1.0; this->cee[i] = 1.0; } diff --git a/Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h b/Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h index 2a7e7535..18144356 100644 --- a/Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h +++ b/Grid/qcd/action/fermion/implementation/WilsonCloverFermionImplementation.h @@ -61,7 +61,7 @@ WilsonCloverFermion::WilsonCloverFermion(GaugeField& diag_mass = _mass + 1.0 + (Nd - 1) * (clover_anisotropy.nu / clover_anisotropy.xi_0); } else { csw_r = _csw_r * 0.5; - diag_mass = 4.0 + _mass; + diag_mass = Nd*1.0 + _mass; } csw_t = _csw_t * 0.5; @@ -297,9 +297,9 @@ void WilsonCloverFermion::MDeriv(GaugeField &force, const F { if (mu == nu) continue; - + RealD factor; - if (nu == 4 || mu == 4) + if (nu == (Nd-1) || mu == (Nd-1)) // This was a bug - surely mu/nu is NEVER 4 but rather (Nd-1)=3 ?? { factor = 2.0 * csw_t; } @@ -307,9 +307,11 @@ void WilsonCloverFermion::MDeriv(GaugeField &force, const F { factor = 2.0 * csw_r; } - PropagatorField Slambda = Gamma(sigma[count]) * Lambda; // sigma checked - Impl::TraceSpinImpl(lambda, Slambda); // traceSpin ok - force_mu -= factor*CloverHelpers::Cmunu(U, lambda, mu, nu); // checked + if ( mu < Nd && nu < Nd ) { // Allow to restrict range to Nd=3, but preserve orders of SigmaMuNu in table by counting ALL + PropagatorField Slambda = Gamma(sigma[count]) * Lambda; // sigma checked + Impl::TraceSpinImpl(lambda, Slambda); // traceSpin ok + force_mu -= factor*CloverHelpers::Cmunu(U, lambda, mu, nu); // checked + } count++; } diff --git a/Grid/qcd/action/fermion/implementation/WilsonFermion5DImplementation.h b/Grid/qcd/action/fermion/implementation/WilsonFermion5DImplementation.h index 9598552f..75fbcaa6 100644 --- a/Grid/qcd/action/fermion/implementation/WilsonFermion5DImplementation.h +++ b/Grid/qcd/action/fermion/implementation/WilsonFermion5DImplementation.h @@ -63,10 +63,10 @@ WilsonFermion5D::WilsonFermion5D(GaugeField &_Umu, Dirichlet(0) { // some assertions - assert(FiveDimGrid._ndimension==5); - assert(FourDimGrid._ndimension==4); - assert(FourDimRedBlackGrid._ndimension==4); - assert(FiveDimRedBlackGrid._ndimension==5); + assert(FiveDimGrid._ndimension==Nd+1); + assert(FourDimGrid._ndimension==Nd); + assert(FourDimRedBlackGrid._ndimension==Nd); + assert(FiveDimRedBlackGrid._ndimension==Nd+1); assert(FiveDimRedBlackGrid._checker_dim==1); // Don't checker the s direction // extent of fifth dim and not spread out @@ -76,7 +76,7 @@ WilsonFermion5D::WilsonFermion5D(GaugeField &_Umu, assert(FiveDimRedBlackGrid._processors[0] ==1); // Other dimensions must match the decomposition of the four-D fields - for(int d=0;d<4;d++){ + for(int d=0;d::WilsonFermion5D(GaugeField &_Umu, if ( p.dirichlet.size() == Nd+1) { Coordinate block = p.dirichlet; - if ( block[0] || block[1] || block[2] || block[3] || block[4] ){ - Dirichlet = 1; - std::cout << GridLogMessage << " WilsonFermion: non-trivial Dirichlet condition "<< block << std::endl; - std::cout << GridLogMessage << " WilsonFermion: partial Dirichlet "<< p.partialDirichlet << std::endl; - Block = block; + for(int d=0;d::WilsonFermion5D(GaugeField &_Umu, assert(FiveDimGrid._simd_layout[0] ==nsimd); assert(FiveDimRedBlackGrid._simd_layout[0]==nsimd); - for(int d=0;d<4;d++){ + for(int d=0;d::DhopDir(const FermionField &in, FermionField &out,in // assert( (dir>=0)&&(dir<4) ); //must do x,y,z or t; int skip = (disp==1) ? 0 : 1; - int dirdisp = dir+skip*4; - int gamma = dir+(1-skip)*4; + int dirdisp = dir+skip*Nd; + int gamma = dir+(1-skip)*Nd; Compressor compressor(DaggerNo); Stencil.HaloExchange(in,compressor); @@ -483,7 +485,7 @@ void WilsonFermion5D::DW(const FermionField &in, FermionField &out,int dag { out.Checkerboard()=in.Checkerboard(); Dhop(in,out,dag); // -0.5 is included - axpy(out,4.0-M5,in,out); + axpy(out,Nd*1.0-M5,in,out); } template void WilsonFermion5D::Meooe(const FermionField &in, FermionField &out) @@ -509,7 +511,7 @@ template void WilsonFermion5D::Mooee(const FermionField &in, FermionField &out) { out.Checkerboard() = in.Checkerboard(); - typename FermionField::scalar_type scal(4.0 + M5); + typename FermionField::scalar_type scal(Nd*1.0 + M5); out = scal * in; } @@ -524,7 +526,7 @@ template void WilsonFermion5D::MooeeInv(const FermionField &in, FermionField &out) { out.Checkerboard() = in.Checkerboard(); - out = (1.0/(4.0 + M5))*in; + out = (1.0/(Nd*1.0 + M5))*in; } template @@ -635,7 +637,7 @@ void WilsonFermion5D::MomentumSpacePropagatorHt_5d(FermionField &out,const A = one / (abs(W) * sinha * 2.0) * one / (sinhaLs * 2.0); F = eaLs * (one - Wea + (Wema - one) * mass*mass); F = F + emaLs * (Wema - one + (one - Wea) * mass*mass); - F = F - abs(W) * sinha * 4.0 * mass; + F = F - abs(W) * sinha * (Nd* 1.0) * mass; Bpp = (A/F) * (ema2Ls - one) * (one - Wema) * (one - mass*mass * one); Bmm = (A/F) * (one - ea2Ls) * (one - Wea) * (one - mass*mass * one); diff --git a/Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h b/Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h index 8c58f692..313d1a95 100644 --- a/Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h +++ b/Grid/qcd/action/fermion/implementation/WilsonFermionImplementation.h @@ -63,7 +63,7 @@ WilsonFermion::WilsonFermion(GaugeField &_Umu, GridCartesian &Fgrid, if (anisotropyCoeff.isAnisotropic){ diag_mass = mass + 1.0 + (Nd-1)*(anisotropyCoeff.nu / anisotropyCoeff.xi_0); } else { - diag_mass = 4.0 + mass; + diag_mass = Nd*1.0 + mass; } int vol4; @@ -354,8 +354,8 @@ void WilsonFermion::DhopDir(const FermionField &in, FermionField &out, int Stencil.HaloExchange(in, compressor); int skip = (disp == 1) ? 0 : 1; - int dirdisp = dir + skip * 4; - int gamma = dir + (1 - skip) * 4; + int dirdisp = dir + skip * Nd; + int gamma = dir + (1 - skip) * Nd; DhopDirCalc(in, out, dirdisp, gamma, DaggerNo); }; @@ -370,8 +370,8 @@ void WilsonFermion::DhopDirAll(const FermionField &in, std::vector distance = st._distances[DIR]; \ sl = st._simd_layout[direction]; \ inplace_twist = 0; \ - if(SE->_around_the_world && st.parameters.twists[DIR % 4]){ \ + if(SE->_around_the_world && st.parameters.twists[DIR % Nd]){ \ if(sl == 1){ \ g = (F+1) % 2; \ }else{ \ diff --git a/Grid/qcd/action/fermion/instantiation/ImprovedStaggeredFermion5DInstantiation.cc b/Grid/qcd/action/fermion/instantiation/ImprovedStaggeredFermion5DInstantiation.cc index 986460b0..596cf74c 100644 --- a/Grid/qcd/action/fermion/instantiation/ImprovedStaggeredFermion5DInstantiation.cc +++ b/Grid/qcd/action/fermion/instantiation/ImprovedStaggeredFermion5DInstantiation.cc @@ -32,8 +32,30 @@ Author: Peter Boyle NAMESPACE_BEGIN(Grid); // S-direction is INNERMOST and takes no part in the parity. -const std::vector ImprovedStaggeredFermion5DStatic::directions({1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4}); -const std::vector ImprovedStaggeredFermion5DStatic::displacements({1, 1, 1, 1, -1, -1, -1, -1, 3, 3, 3, 3, -3, -3, -3, -3}); +const std::vector ImprovedStaggeredFermion5DStatic::directions(ImprovedStaggeredFermion5DStatic::MakeDirections()); +const std::vector ImprovedStaggeredFermion5DStatic::displacements(ImprovedStaggeredFermion5DStatic::MakeDisplacements()); +std::vector ImprovedStaggeredFermion5DStatic::MakeDirections(void) +{ + std::vector directions(4*Nd); + for(int d=0;d ImprovedStaggeredFermion5DStatic::MakeDisplacements(void) +{ + std::vector displacements(4*Nd); + for(int d=0;d ImprovedStaggeredFermionStatic::directions({0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3}); const std::vector ImprovedStaggeredFermionStatic::displacements({1, 1, 1, 1, -1, -1, -1, -1, 3, 3, 3, 3, -3, -3, -3, -3}); - +std::vector ImprovedStaggeredFermionStatic::MakeDirections(void) +{ + std::vector directions(4*Nd); + for(int d=0;d ImprovedStaggeredFermionStatic::MakeDisplacements(void) +{ + std::vector displacements(4*Nd); + for(int d=0;d NaiveStaggeredFermionStatic::directions({0, 1, 2, 3, 0, 1, 2, 3}); -const std::vector NaiveStaggeredFermionStatic::displacements({1, 1, 1, 1, -1, -1, -1, -1}); +//const std::vector NaiveStaggeredFermionStatic::directions({0, 1, 2, 3, 0, 1, 2, 3}); +//const std::vector NaiveStaggeredFermionStatic::displacements({1, 1, 1, 1, -1, -1, -1, -1}); +const std::vector NaiveStaggeredFermionStatic::directions(NaiveStaggeredFermionStatic::MakeDirections()); +const std::vector NaiveStaggeredFermionStatic::displacements(NaiveStaggeredFermionStatic::MakeDisplacements()); +std::vector NaiveStaggeredFermionStatic::MakeDirections(void) +{ + std::vector directions(4*Nd); + for(int d=0;d NaiveStaggeredFermionStatic::MakeDisplacements(void) +{ + std::vector displacements(4*Nd); + for(int d=0;d WilsonFermion5DStatic::directions ({1,2,3,4, 1, 2, 3, 4}); -const std::vector WilsonFermion5DStatic::displacements({1,1,1,1,-1,-1,-1,-1}); + +const std::vector WilsonFermion5DStatic::directions (WilsonFermion5DStatic::MakeDirections()); +const std::vector WilsonFermion5DStatic::displacements(WilsonFermion5DStatic::MakeDisplacements()); + +std::vector WilsonFermion5DStatic::MakeDirections (void) +{ + std::vector directions(2*Nd); + for(int d=0;d WilsonFermion5DStatic::MakeDisplacements(void) +{ + std::vector displacements(2*Nd); + for(int d=0;d WilsonFermionStatic::directions({0, 1, 2, 3, 0, 1, 2, 3}); -const std::vector WilsonFermionStatic::displacements({1, 1, 1, 1, -1, -1, -1, -1}); +const std::vector WilsonFermionStatic::directions(WilsonFermionStatic::MakeDirections()); +const std::vector WilsonFermionStatic::displacements(WilsonFermionStatic::MakeDisplacements()); int WilsonFermionStatic::HandOptDslash; +std::vector WilsonFermionStatic::MakeDirections (void) +{ + std::vector directions(2*Nd); + for(int d=0;d WilsonFermionStatic::MakeDisplacements(void) +{ + std::vector displacements(2*Nd); + for(int d=0;d::energyDensityCloverleaf(const RealD t, const GaugeF LatticeComplexD R(U.Grid()); R = Zero(); - for(int mu=0;mu<3;mu++){ - for(int nu=mu+1;nu<4;nu++){ + for(int mu=0;mu::FieldStrength(F, U, mu, nu); R = R + trace(F*F); } diff --git a/Grid/qcd/utils/WilsonLoops.h b/Grid/qcd/utils/WilsonLoops.h index 7466f4bf..f9dbdb37 100644 --- a/Grid/qcd/utils/WilsonLoops.h +++ b/Grid/qcd/utils/WilsonLoops.h @@ -179,20 +179,17 @@ public: ////////////////////////////////////////////////// // average over all x,y,z the temporal loop ////////////////////////////////////////////////// - static ComplexD avgPolyakovLoop(const GaugeField &Umu) { //assume Nd=4 + static ComplexD avgPolyakovLoop(const GaugeField &Umu) { GaugeMat Ut(Umu.Grid()), P(Umu.Grid()); ComplexD out; - int T = Umu.Grid()->GlobalDimensions()[3]; - int X = Umu.Grid()->GlobalDimensions()[0]; - int Y = Umu.Grid()->GlobalDimensions()[1]; - int Z = Umu.Grid()->GlobalDimensions()[2]; - - Ut = peekLorentz(Umu,3); //Select temporal direction + uint64_t vol = Umu.Grid()->gSites(); + int T = Umu.Grid()->GlobalDimensions()[Nd-1]; + Ut = peekLorentz(Umu,Nd-1); //Select temporal direction P = Ut; for (int t=1;tgSites(); - return p.real() / vol / (4.0 * Nc ) ; + return p.real() / vol / (Nd * Nc ) ; }; ////////////////////////////////////////////////// @@ -740,6 +737,7 @@ public: //cf https://arxiv.org/pdf/hep-lat/9701012.pdf Eq 6 //output is the charge by timeslice: sum over timeslices to obtain the total static std::vector TimesliceTopologicalChargeMxN(const GaugeLorentz &U, int M, int N){ + // Audit: 4D epsilon is hard coded assert(Nd == 4); std::vector > F(Nd,std::vector(Nd,nullptr)); //Note F_numu = - F_munu @@ -829,6 +827,25 @@ public: return out; } + //Compute the 5Li topological charge density + static std::vector TopologicalChargeDensity5Li(const GaugeLorentz &U){ + + static const int exts[5][2] = { {1,1}, {2,2}, {1,2}, {1,3}, {3,3} }; + std::vector > loops = TimesliceTopologicalCharge5LiContributions(U); + + double c5=1./20.; + double c4=1./5.-2.*c5; + double c3=(-64.+640.*c5)/45.; + double c2=(1-64.*c5)/9.; + double c1=(19.-55.*c5)/9.; + + int Lt = loops[0].size(); + std::vector out(Lt,0.); + for(int t=0;t Qt = TimesliceTopologicalCharge5Li(U); Real Q = 0.; @@ -1455,7 +1472,7 @@ public: ////////////////////////////////////////////////// static Real sumWilsonLoop(const GaugeLorentz &Umu, const int R1, const int R2) { - std::vector U(4, Umu.Grid()); + std::vector U(Nd, Umu.Grid()); for (int mu = 0; mu < Umu.Grid()->_ndimension; mu++) { U[mu] = PeekIndex(Umu, mu); @@ -1474,7 +1491,7 @@ public: ////////////////////////////////////////////////// static Real sumTimelikeWilsonLoop(const GaugeLorentz &Umu, const int R1, const int R2) { - std::vector U(4, Umu.Grid()); + std::vector U(Nd, Umu.Grid()); for (int mu = 0; mu < Umu.Grid()->_ndimension; mu++) { U[mu] = PeekIndex(Umu, mu); @@ -1492,8 +1509,8 @@ public: // sum over all x,y,z,t and over all planes of spatial Wilson loop ////////////////////////////////////////////////// static Real sumSpatialWilsonLoop(const GaugeLorentz &Umu, - const int R1, const int R2) { - std::vector U(4, Umu.Grid()); + const int R1, const int R2) { + std::vector U(Nd, Umu.Grid()); for (int mu = 0; mu < Umu.Grid()->_ndimension; mu++) { U[mu] = PeekIndex(Umu, mu); diff --git a/Grid/util/Init.cc b/Grid/util/Init.cc index 43d959c0..d8d34373 100644 --- a/Grid/util/Init.cc +++ b/Grid/util/Init.cc @@ -187,8 +187,8 @@ void GridParseLayout(char **argv,int argc, Coordinate &latt_c, Coordinate &mpi_c) { - auto mpi =std::vector({1,1,1,1}); - auto latt=std::vector({8,8,8,8}); + auto mpi =std::vector(1,Nd); + auto latt=std::vector(8,Nd); GridThread::SetMaxThreads(); diff --git a/configure.ac b/configure.ac index 9664a675..1471a0d1 100644 --- a/configure.ac +++ b/configure.ac @@ -198,6 +198,8 @@ AC_ARG_ENABLE([Nc], [ac_Nc=${enable_Nc}], [ac_Nc=3]) case ${ac_Nc} in + 1) + AC_DEFINE([Config_Nc],[1],[Gauge group Nc]);; 2) AC_DEFINE([Config_Nc],[2],[Gauge group Nc]);; 3) @@ -211,6 +213,21 @@ case ${ac_Nc} in *) AC_MSG_ERROR(["Unsupport gauge group choice Nc = ${ac_Nc}"]);; esac +############### Nd +AC_ARG_ENABLE([Nd], + [AS_HELP_STRING([--enable-Nd=2|3|4],[enable default LGT dimension])], + [ac_Nd=${enable_Nd}], [ac_Nd=4]) + +case ${ac_Nd} in + 2) + AC_DEFINE([Config_Nd],[2],[Gauge field dimension Nd]);; + 3) + AC_DEFINE([Config_Nd],[3],[Gauge field dimension Nd]);; + 4) + AC_DEFINE([Config_Nd],[4],[Gauge field dimension Nd]);; + *) + AC_MSG_ERROR(["Unsupport gauge group choice Nc = ${ac_Nc}"]);; +esac ############### Symplectic group AC_ARG_ENABLE([Sp],