diff --git a/lib/qcd/action/fermion/WilsonCloverFermion.cc b/lib/qcd/action/fermion/WilsonCloverFermion.cc index e678835a..fff970a2 100644 --- a/lib/qcd/action/fermion/WilsonCloverFermion.cc +++ b/lib/qcd/action/fermion/WilsonCloverFermion.cc @@ -43,12 +43,11 @@ RealD WilsonCloverFermion::M(const FermionField &in, FermionField &out) // Wilson term out.checkerboard = in.checkerboard; - //this->Dhop(in, out, DaggerNo); + this->Dhop(in, out, DaggerNo); // Clover term Mooee(in, temp); - out= zero; out += temp; return norm2(out); } @@ -60,12 +59,11 @@ RealD WilsonCloverFermion::Mdag(const FermionField &in, FermionField &out) // Wilson term out.checkerboard = in.checkerboard; - //this->Dhop(in, out, DaggerYes); + this->Dhop(in, out, DaggerYes); // Clover term MooeeDag(in, temp); - out=zero; out += temp; return norm2(out); } diff --git a/lib/qcd/action/fermion/WilsonCloverFermion.h b/lib/qcd/action/fermion/WilsonCloverFermion.h index 402a9a7e..cd13b225 100644 --- a/lib/qcd/action/fermion/WilsonCloverFermion.h +++ b/lib/qcd/action/fermion/WilsonCloverFermion.h @@ -116,7 +116,7 @@ public: force = zero; // Derivative of the Wilson hopping term - //this->DhopDeriv(force, X, Y, dag); + this->DhopDeriv(force, X, Y, dag); /////////////////////////////////////////////////////////// // Clover term derivative diff --git a/tests/forces/Test_wilsonclover_force.cc b/tests/forces/Test_wilsonclover_force.cc index 82adb8ab..bcf67be4 100644 --- a/tests/forces/Test_wilsonclover_force.cc +++ b/tests/forces/Test_wilsonclover_force.cc @@ -48,13 +48,12 @@ int main(int argc, char **argv) std::vector seeds({1, 2, 30, 50}); GridParallelRNG pRNG(&Grid); - + std::vector vrand(4); std::srand(std::time(0)); std::generate(vrand.begin(), vrand.end(), std::rand); std::cout << GridLogMessage << vrand << std::endl; pRNG.SeedFixedIntegers(vrand); - //pRNG.SeedFixedIntegers(seeds); LatticeFermion phi(&Grid); @@ -64,50 +63,14 @@ int main(int argc, char **argv) LatticeGaugeField U(&Grid); -/* - std::vector x(4); // 4d fermions - std::vector gd = Grid.GlobalDimensions(); - Grid::QCD::SpinColourVector F; - Grid::Complex c; - - phi = zero; - for (x[0] = 0; x[0] < 1; x[0]++) - { - for (x[1] = 0; x[1] < 1; x[1]++) - { - for (x[2] = 0; x[2] < 1; x[2]++) - { - for (x[3] = 0; x[3] < 1; x[3]++) - { - for (int sp = 0; sp < 4; sp++) - { - for (int j = 0; j < 3; j++) // colours - { - F()(sp)(j) = Grid::Complex(0.0,0.0); - if (((sp == 0) && (j==0))) - { - c = Grid::Complex(1.0, 0.0); - F()(sp)(j) = c; - } - } - } - Grid::pokeSite(F, phi, x); - - } - } - } - } -*/ - std::vector site = {0, 0, 0, 0}; SU3::HotConfiguration(pRNG, U); - //SU3::ColdConfiguration(pRNG, U); - + //SU3::ColdConfiguration(pRNG, U);// Clover term zero //////////////////////////////////// // Unmodified matrix element //////////////////////////////////// - RealD mass = -4.0; //kills the diagonal term + RealD mass = 0.1; Real csw = 1.0; WilsonCloverFermionR Dw(U, Grid, RBGrid, mass, csw); Dw.ImportGauge(U); @@ -125,103 +88,42 @@ int main(int argc, char **argv) UdSdU += tmp; ///////////////////////////////////////////// - // Take the traceless antihermitian component - //UdSdU = Ta(UdSdU); - - //std::cout << UdSdU << std::endl; - //SU3::LatticeAlgebraVector hforce(&Grid); - LatticeColourMatrix mommu(&Grid); - //mommu = PeekIndex(UdSdU, 0); - //SU3::projectOnAlgebra(hforce, mommu); - //std::cout << hforce << std::endl; - //////////////////////////////////// // Modify the gauge field a little //////////////////////////////////// - RealD dt = 0.0001; + RealD dt = 0.00005; RealD Hmom = 0.0; RealD Hmomprime = 0.0; RealD Hmompp = 0.0; + LatticeColourMatrix mommu(&Grid); LatticeColourMatrix forcemu(&Grid); LatticeGaugeField mom(&Grid); LatticeGaugeField Uprime(&Grid); - - for (int mu = 0; mu < Nd; mu++) { + for (int mu = 0; mu < Nd; mu++) + { // Traceless antihermitian momentum; gaussian in lie alg SU3::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu); Hmom -= real(sum(trace(mommu * mommu))); PokeIndex(mom, mommu, mu); - } - /* - SU3::AlgebraVector h; - SU3::LatticeAlgebraVector hl(&Grid); - h()()(0) = 1.0; - hl = zero; - pokeSite(h, hl, site); - SU3::FundamentalLieAlgebraMatrix(hl, mommu); - mom = zero; - PokeIndex(mom, mommu, 0); - Hmom -= real(sum(trace(mommu * mommu))); - */ - /* - parallel_for(int ss=0;ssoSites();ss++){ - for (int mu = 0; mu < Nd; mu++) - Uprime[ss]._internal[mu] = ProjectOnGroup(Exponentiate(mom[ss]._internal[mu], dt, 12) * U[ss]._internal[mu]); - } -*/ - - for (int mu = 0; mu < Nd; mu++) - { - parallel_for(auto i = mom.begin(); i < mom.end(); i++) + parallel_for(int ss = 0; ss < mom._grid->oSites(); ss++) { - Uprime[i](mu) = U[i](mu); - Uprime[i](mu) += mom[i](mu) * U[i](mu) * dt; - Uprime[i](mu) += mom[i](mu) * mom[i](mu) * U[i](mu) * (dt * dt / 2.0); - Uprime[i](mu) += mom[i](mu) * mom[i](mu) * mom[i](mu) * U[i](mu) * (dt * dt * dt / 6.0); - Uprime[i](mu) += mom[i](mu) * mom[i](mu) * mom[i](mu) * mom[i](mu) * U[i](mu) * (dt * dt * dt * dt / 24.0); - Uprime[i](mu) += mom[i](mu) * mom[i](mu) * mom[i](mu) * mom[i](mu) * mom[i](mu) * U[i](mu) * (dt * dt * dt * dt * dt / 120.0); - Uprime[i](mu) += mom[i](mu) * mom[i](mu) * mom[i](mu) * mom[i](mu) * mom[i](mu) * mom[i](mu) * U[i](mu) * (dt * dt * dt * dt * dt * dt / 720.0); + Uprime[ss]._internal[mu] = ProjectOnGroup(Exponentiate(mom[ss]._internal[mu], dt, 12) * U[ss]._internal[mu]); } } std::cout << GridLogMessage << "Initial mom hamiltonian is " << Hmom << std::endl; // New action - LatticeGaugeField diff(&Grid); - diff = Uprime - U; - //std::cout << "Diff:" << diff << std::endl; Dw.ImportGauge(Uprime); Dw.M(phi, MphiPrime); - LatticeFermion DiffFermion(&Grid); - DiffFermion = MphiPrime - Mphi; - //std::cout << "DiffFermion:" << DiffFermion << std::endl; - //std::cout << "Mphi:" << Mphi << std::endl; - //std::cout << "MphiPrime:" << MphiPrime << std::endl; - ComplexD Sprime = innerProduct(MphiPrime, MphiPrime); ////////////////////////////////////////////// // Use derivative to estimate dS ////////////////////////////////////////////// - /////////////////////////////////////////////////////// - std::cout << GridLogMessage << "Antihermiticity tests - 1 " << std::endl; - for (int mu = 0; mu < Nd; mu++) - { - mommu = PeekIndex(mom, mu); - std::cout << GridLogMessage << " Mommu " << norm2(mommu) << std::endl; - mommu = mommu + adj(mommu); - std::cout << GridLogMessage << " Test: Mommu + Mommudag " << norm2(mommu) << std::endl; - mommu = PeekIndex(UdSdU, mu); - std::cout << GridLogMessage << " dsdumu " << norm2(mommu) << std::endl; - mommu = mommu + adj(mommu); - std::cout << GridLogMessage << " Test: dsdumu + dag " << norm2(mommu) << std::endl; - std::cout << "" << std::endl; - } - //////////////////////////////////////////////////////// - LatticeComplex dS(&Grid); dS = zero; LatticeComplex dSmom(&Grid); @@ -229,7 +131,6 @@ int main(int argc, char **argv) LatticeComplex dSmom2(&Grid); dSmom2 = zero; - for (int mu = 0; mu < Nd; mu++) { mommu = PeekIndex(UdSdU, mu); // P_mu = @@ -237,7 +138,7 @@ int main(int argc, char **argv) PokeIndex(UdSdU, mommu, mu); // UdSdU_mu = Mom } - std::cout << GridLogMessage << "Antihermiticity tests - 2 " << std::endl; + std::cout << GridLogMessage << "Antihermiticity tests" << std::endl; for (int mu = 0; mu < Nd; mu++) { mommu = PeekIndex(mom, mu); @@ -279,8 +180,8 @@ int main(int argc, char **argv) std::cout << GridLogMessage << " S " << S << std::endl; std::cout << GridLogMessage << " Sprime " << Sprime << std::endl; - std::cout << GridLogMessage << "dS " << Sprime - S << std::endl; - std::cout << GridLogMessage << "predict dS " << dSpred << std::endl; + std::cout << GridLogMessage << "dS (S' - S) :" << Sprime - S << std::endl; + std::cout << GridLogMessage << "predict dS (force) :" << dSpred << std::endl; std::cout << GridLogMessage << "dSm " << dSm << std::endl; std::cout << GridLogMessage << "dSm2" << dSm2 << std::endl;