/* * Warning: This code illustrative only: not well tested, and not meant for production use * without regression / tests being applied */ #include using namespace std; using namespace Grid; RealD LLscale =1.0; RealD LCscale =1.0; template class CovariantLaplacianCshift : public SparseMatrixBase { public: INHERIT_GIMPL_TYPES(Gimpl); GridBase *grid; GaugeField U; CovariantLaplacianCshift(GaugeField &_U) : grid(_U.Grid()), U(_U) { }; virtual GridBase *Grid(void) { return grid; }; virtual void M (const Field &in, Field &out) { out=Zero(); for(int mu=0;mu(U, mu); // NB: Inefficent out = out - Gimpl::CovShiftForward(Umu,mu,in); out = out - Gimpl::CovShiftBackward(Umu,mu,in); out = out + 2.0*in; } }; virtual void Mdag (const Field &in, Field &out) { M(in,out);}; // Laplacian is hermitian virtual void Mdiag (const Field &in, Field &out) {assert(0);}; // Unimplemented need only for multigrid virtual void Mdir (const Field &in, Field &out,int dir, int disp){assert(0);}; // Unimplemented need only for multigrid virtual void MdirAll (const Field &in, std::vector &out) {assert(0);}; // Unimplemented need only for multigrid }; void MakePhase(Coordinate mom,LatticeComplex &phase) { GridBase *grid = phase.Grid(); auto latt_size = grid->GlobalDimensions(); ComplexD ci(0.0,1.0); phase=Zero(); LatticeComplex coor(phase.Grid()); for(int mu=0;mu void GaussianSmear(LatticeGaugeField &U,Field &unsmeared,Field &smeared) { typedef CovariantLaplacianCshift Laplacian_t; Laplacian_t Laplacian(U); Integer Iterations = 40; Real width = 2.0; Real coeff = (width*width) / Real(4*Iterations); Field tmp(U.Grid()); smeared=unsmeared; // chi = (1-p^2/2N)^N kronecker for(int n = 0; n < Iterations; ++n) { Laplacian.M(smeared,tmp); smeared = smeared - coeff*tmp; std::cout << " smear iter " << n<<" " < void MasslessFreePropagator(Action &D,LatticePropagator &source,LatticePropagator &propagator) { GridBase *UGrid = source.Grid(); GridBase *FGrid = D.FermionGrid(); bool fiveD = true; //calculate 5d free propagator RealD mass = D.Mass(); LatticeFermion src4 (UGrid); LatticeFermion result4 (UGrid); LatticeFermion result5(FGrid); LatticeFermion src5(FGrid); LatticePropagator prop5(FGrid); for(int s=0;s(src4,source,s,c); D.ImportPhysicalFermionSource(src4,src5); D.FreePropagator(src5,result5,mass,true); std::cout<(prop5,result5,s,c); FermToProp(propagator,result4,s,c); } } LatticePropagator Vector_mu(UGrid); LatticeComplex VV (UGrid); std::vector sumVV; Gamma::Algebra GammaV[3] = { Gamma::Algebra::GammaX, Gamma::Algebra::GammaY, Gamma::Algebra::GammaZ }; for( int mu=0;mu<3;mu++ ) { Gamma gV(GammaV[mu]); D.ContractConservedCurrent(prop5,prop5,Vector_mu,source,Current::Vector,mu); VV = trace(gV*Vector_mu); // (local) Vector-Vector conserved current sliceSum(VV,sumVV,Tdir); int Nt = sumVV.size(); for(int t=0;t void MasslessFreePropagator1(Action &D,LatticePropagator &source,LatticePropagator &propagator) { bool fiveD = false; //calculate 4d free propagator RealD mass = D.Mass(); GridBase *UGrid = source.Grid(); LatticeFermion src4 (UGrid); LatticeFermion result4 (UGrid); for(int s=0;s(src4,source,s,c); D.FreePropagator(src4,result4,mass,false); FermToProp(propagator,result4,s,c); } } } template void Solve(Action &D,LatticePropagator &source,LatticePropagator &propagator) { GridBase *UGrid = D.GaugeGrid(); GridBase *FGrid = D.FermionGrid(); LatticeFermion src4 (UGrid); LatticeFermion src5 (FGrid); LatticeFermion result5(FGrid); LatticeFermion result4(UGrid); LatticePropagator prop5(FGrid); ConjugateGradient CG(1.0e-10,100000); SchurRedBlackDiagMooeeSolve schur(CG); ZeroGuesser ZG; // Could be a DeflatedGuesser if have eigenvectors for(int s=0;s(src4,source,s,c); D.ImportPhysicalFermionSource(src4,src5); result5=Zero(); schur(D,src5,result5,ZG); std::cout<(prop5,result5,s,c); FermToProp(propagator,result4,s,c); } } LatticePropagator Axial_mu(UGrid); LatticePropagator Vector_mu(UGrid); LatticeComplex PA (UGrid); LatticeComplex VV (UGrid); LatticeComplex PJ5q(UGrid); LatticeComplex PP (UGrid); std::vector sumPA; std::vector sumVV; std::vector sumPP; std::vector sumPJ5q; Gamma g5(Gamma::Algebra::Gamma5); D.ContractConservedCurrent(prop5,prop5,Axial_mu,source,Current::Axial,Tdir); PA = trace(g5*Axial_mu); // Pseudoscalar-Axial conserved current sliceSum(PA,sumPA,Tdir); int Nt{static_cast(sumPA.size())}; for(int t=0;t >, data); }; void MesonTrace(std::string file,LatticePropagator &q1,LatticePropagator &q2,LatticeComplex &phase) { const int nchannel=4; Gamma::Algebra Gammas[nchannel][2] = { {Gamma::Algebra::GammaXGamma5,Gamma::Algebra::GammaXGamma5}, {Gamma::Algebra::GammaYGamma5,Gamma::Algebra::GammaYGamma5}, {Gamma::Algebra::GammaZGamma5,Gamma::Algebra::GammaZGamma5}, {Gamma::Algebra::Identity,Gamma::Algebra::Identity} }; LatticeComplex meson_CF(q1.Grid()); MesonFile MF; for(int ch=0;ch meson_T; sliceSum(meson_CF,meson_T, Tdir); int nt=meson_T.size(); std::vector corr(nt); for(int t=0;t seeds4({1,2,3,4}); // GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); LatticeGaugeField Umu(UGrid); std::string config; RealD M5=atof(getenv("M5")); RealD mq = atof(getenv("mass")); int tadpole = atof(getenv("tadpole")); std::vector masses({ mq} ); // u/d, s, c ?? if( argc > 1 && argv[1][0] != '-' ) { std::cout<::ColdConfiguration(Umu); config="ColdConfig"; // RealD P=1.0; // Don't scale // RealD P=0.6388238 // 32Ifine // RealD P=0.6153342; // 64I RealD P=0.5871119; // 48I RealD u0 = sqrt(sqrt(P)); RealD w0 = 1 - M5; std::cout< boundary = {1,1,1,-1}; FermionActionR::ImplParams Params(boundary); RealD b=1.5; RealD c=0.5; std::cout< PointProps(nmass,UGrid); // std::vector FreeProps(nmass,UGrid); // LatticePropagator delta(UGrid); for(int m=0;m