diff --git a/lib/qcd/action/fermion/WilsonCloverFermion.cc b/lib/qcd/action/fermion/WilsonCloverFermion.cc index 7f58f277..5e7e0034 100644 --- a/lib/qcd/action/fermion/WilsonCloverFermion.cc +++ b/lib/qcd/action/fermion/WilsonCloverFermion.cc @@ -71,10 +71,10 @@ void WilsonCloverFermion::ImportGauge(const GaugeField &_Umu) GridBase *grid = _Umu._grid; typename Impl::GaugeLinkField Bx(grid), By(grid), Bz(grid), Ex(grid), Ey(grid), Ez(grid); - // Compute the field strength terms - WilsonLoops::FieldStrength(Bx, _Umu, Ydir, Zdir); + // Compute the field strength terms mu>nu + WilsonLoops::FieldStrength(Bx, _Umu, Zdir, Ydir); WilsonLoops::FieldStrength(By, _Umu, Zdir, Xdir); - WilsonLoops::FieldStrength(Bz, _Umu, Xdir, Ydir); + WilsonLoops::FieldStrength(Bz, _Umu, Ydir, Xdir); WilsonLoops::FieldStrength(Ex, _Umu, Tdir, Xdir); WilsonLoops::FieldStrength(Ey, _Umu, Tdir, Ydir); WilsonLoops::FieldStrength(Ez, _Umu, Tdir, Zdir); @@ -86,7 +86,7 @@ void WilsonCloverFermion::ImportGauge(const GaugeField &_Umu) CloverTerm += fillCloverXT(Ex); CloverTerm += fillCloverYT(Ey); CloverTerm += fillCloverZT(Ez); - CloverTerm *= 0.5 * csw; // FieldStrength normalization? should be ( -i/8 ). Is it the anti-symmetric combination? + CloverTerm *= (0.5) * csw; int lvol = _Umu._grid->lSites(); @@ -232,7 +232,8 @@ void WilsonCloverFermion::MooeeInternal(const FermionField &in, FermionFie out = *Clover * in; } else { Clover = (inv) ? &CloverTermInv : &CloverTerm; - out = adj(*Clover) * in; + //out = adj(*Clover) * in; + out = adj(CloverTerm) * in; } diff --git a/lib/qcd/action/fermion/WilsonCloverFermion.h b/lib/qcd/action/fermion/WilsonCloverFermion.h index 7840af90..18386485 100644 --- a/lib/qcd/action/fermion/WilsonCloverFermion.h +++ b/lib/qcd/action/fermion/WilsonCloverFermion.h @@ -156,10 +156,10 @@ private: PARALLEL_FOR_LOOP for (int i = 0; i < CloverTerm._grid->oSites(); i++) { - T._odata[i]()(0, 1) = timesI(F._odata[i]()()); //fixed - T._odata[i]()(1, 0) = timesI(F._odata[i]()()); //fixed - T._odata[i]()(2, 3) = timesMinusI(F._odata[i]()()); //fixed - T._odata[i]()(3, 2) = timesMinusI(F._odata[i]()()); //fixed + T._odata[i]()(0, 1) = timesI(F._odata[i]()()); + T._odata[i]()(1, 0) = timesI(F._odata[i]()()); + T._odata[i]()(2, 3) = timesMinusI(F._odata[i]()()); + T._odata[i]()(3, 2) = timesMinusI(F._odata[i]()()); } return T; @@ -172,10 +172,10 @@ private: PARALLEL_FOR_LOOP for (int i = 0; i < CloverTerm._grid->oSites(); i++) { - T._odata[i]()(0, 1) = -(F._odata[i]()()); //fixed - T._odata[i]()(1, 0) = (F._odata[i]()()); //fixed - T._odata[i]()(2, 3) = (F._odata[i]()()); //fixed - T._odata[i]()(3, 2) = -(F._odata[i]()()); //fixed + T._odata[i]()(0, 1) = -(F._odata[i]()()); + T._odata[i]()(1, 0) = (F._odata[i]()()); + T._odata[i]()(2, 3) = (F._odata[i]()()); + T._odata[i]()(3, 2) = -(F._odata[i]()()); } return T; @@ -188,10 +188,10 @@ private: PARALLEL_FOR_LOOP for (int i = 0; i < CloverTerm._grid->oSites(); i++) { - T._odata[i]()(0, 0) = timesI(F._odata[i]()()); //fixed - T._odata[i]()(1, 1) = timesMinusI(F._odata[i]()()); //fixed - T._odata[i]()(2, 2) = timesMinusI(F._odata[i]()()); //fixed - T._odata[i]()(3, 3) = timesI(F._odata[i]()()); //fixed + T._odata[i]()(0, 0) = timesI(F._odata[i]()()); + T._odata[i]()(1, 1) = timesMinusI(F._odata[i]()()); + T._odata[i]()(2, 2) = timesMinusI(F._odata[i]()()); + T._odata[i]()(3, 3) = timesI(F._odata[i]()()); } return T; diff --git a/tests/qdpxx/Test_qdpxx_wilson.cc b/tests/qdpxx/Test_qdpxx_wilson.cc index 70a8b1bf..1e1f6a23 100644 --- a/tests/qdpxx/Test_qdpxx_wilson.cc +++ b/tests/qdpxx/Test_qdpxx_wilson.cc @@ -29,7 +29,7 @@ #include // Mass -double mq = 0.1; +double mq = 0.0; // Define Wilson Types typedef Grid::QCD::WilsonImplR::FermionField FermionField; @@ -274,7 +274,7 @@ public: p.Mass = _mq; p.clovCoeffR = QDP::Real(1.0); p.clovCoeffT = QDP::Real(1.0); - Real u0 = QDP::Real(0.0); + Real u0 = QDP::Real(1.0); Chroma::Handle> fbc(new Chroma::SimpleFermBC(bcs)); @@ -316,6 +316,8 @@ int main(int argc, char **argv) FermionField src(UGrid); FermionField res_chroma(UGrid); FermionField res_grid(UGrid); + FermionField only_wilson(UGrid); + FermionField difference(UGrid); std::vector ActionList({Wilson, WilsonClover}); std::vector ActionName({"Wilson", "WilsonClover"}); @@ -346,8 +348,19 @@ int main(int argc, char **argv) std::cout << "Norm of Grid " << ActionName[i] << " multiply " << Grid::norm2(res_grid) << std::endl; - res_chroma = res_chroma - res_grid; - std::cout << "Norm of difference " << Grid::norm2(res_chroma) << std::endl; + difference = res_chroma - res_grid; + std::cout << "Norm of difference " << Grid::norm2(difference) << std::endl; + + // Isolate Clover term + calc_grid(Wilson, Ug, src, only_wilson, dag);// Wilson term + res_grid -= only_wilson; + res_chroma -= only_wilson; + + std::cout << "Chroma:" << res_chroma << std::endl; + std::cout << "Grid :" << res_grid << std::endl; + + + } } @@ -416,7 +429,36 @@ void make_gauge(GaugeField &Umu, FermionField &src) Grid::GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); Grid::QCD::SU3::HotConfiguration(RNG4, Umu); - Grid::gaussian(RNG4, src); + + // Fermion field + //Grid::gaussian(RNG4, src); + Grid::QCD::SpinColourVector F; + Grid::Complex c; + + std::vector x(4); // 4d fermions + std::vector gd = src._grid->GlobalDimensions(); + + for (x[0] = 0; x[0] < gd[0]; x[0]++) + { + for (x[1] = 0; x[1] < gd[1]; x[1]++) + { + for (x[2] = 0; x[2] < gd[2]; x[2]++) + { + for (x[3] = 0; x[3] < gd[3]; x[3]++) + { + for (int sp = 0; sp < 1; sp++) + { + for (int j = 1; j < 2; j++)// colours + { + c = Grid::Complex(1.0, 0.0); + F()(sp)(j) = c; + } + } + Grid::pokeSite(F, src, x); + } + } + } + } } void calc_grid(ChromaAction action, Grid::QCD::LatticeGaugeField &Umu, Grid::QCD::LatticeFermion &src, Grid::QCD::LatticeFermion &res, int dag)