From f941c4ee180aa20b3f3a24a939018357457e5bbf Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Sun, 29 Oct 2017 11:43:33 +0000 Subject: [PATCH] Clover term force ok --- lib/qcd/action/fermion/WilsonCloverFermion.cc | 9 +- lib/qcd/action/fermion/WilsonCloverFermion.h | 24 ++-- lib/qcd/utils/WilsonLoops.h | 2 +- tests/forces/Test_wilsonclover_force.cc | 118 +++++++++++++++--- 4 files changed, 118 insertions(+), 35 deletions(-) diff --git a/lib/qcd/action/fermion/WilsonCloverFermion.cc b/lib/qcd/action/fermion/WilsonCloverFermion.cc index 2159fffc..e678835a 100644 --- a/lib/qcd/action/fermion/WilsonCloverFermion.cc +++ b/lib/qcd/action/fermion/WilsonCloverFermion.cc @@ -48,10 +48,7 @@ RealD WilsonCloverFermion::M(const FermionField &in, FermionField &out) // Clover term Mooee(in, temp); - //hack - out = zero; - - + out= zero; out += temp; return norm2(out); } @@ -68,9 +65,7 @@ RealD WilsonCloverFermion::Mdag(const FermionField &in, FermionField &out) // Clover term MooeeDag(in, temp); - //hack - out = zero; - + out=zero; out += temp; return norm2(out); } diff --git a/lib/qcd/action/fermion/WilsonCloverFermion.h b/lib/qcd/action/fermion/WilsonCloverFermion.h index d8a42129..402a9a7e 100644 --- a/lib/qcd/action/fermion/WilsonCloverFermion.h +++ b/lib/qcd/action/fermion/WilsonCloverFermion.h @@ -122,6 +122,8 @@ public: // Clover term derivative /////////////////////////////////////////////////////////// Impl::outerProductImpl(Lambda, X, Y); + //std::cout << "Lambda:" << Lambda << std::endl; + Gamma::Algebra sigma[] = { Gamma::Algebra::SigmaXY, @@ -153,16 +155,18 @@ public: for (int nu = 0; nu < 4; nu++) { if (mu == nu) continue; - PropagatorField Slambda = Gamma(sigma[count]) * Lambda; - Impl::TraceSpinImpl(lambda, Slambda); //traceSpin - force_mu += Cmunu(U, lambda, mu, nu); + PropagatorField Slambda = Gamma(sigma[count]) * Lambda; // sigma checked + Impl::TraceSpinImpl(lambda, Slambda); // traceSpin ok + force_mu -= Cmunu(U, lambda, mu, nu); // checked count++; } pokeLorentz(clover_force, U[mu] * force_mu, mu); } - clover_force *= csw / 8.; + clover_force *= csw; force += clover_force; + + } // Computing C_{\mu \nu}(x) as in Eq.(B.39) in Zbigniew Sroczynski's PhD thesis @@ -170,20 +174,19 @@ public: { conformable(lambda._grid, U[0]._grid); GaugeLinkField out(lambda._grid), tmp(lambda._grid); - // insertion in upper staple // please check redundancy of shift operations - + // C1+ tmp = lambda * U[nu]; out = Impl::ShiftStaple(Impl::CovShiftForward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu); - + // C2+ - tmp = U[mu] * Impl::CovShiftIdentityForward(adj(lambda), mu); + tmp = U[mu] * Impl::ShiftStaple(adj(lambda), mu); out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(tmp, mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu); - + // C3+ - tmp = U[nu] * Impl::CovShiftIdentityForward(adj(lambda), nu); + tmp = U[nu] * Impl::ShiftStaple(adj(lambda), nu); out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(tmp, nu))), mu); // C4+ @@ -259,6 +262,7 @@ private: PARALLEL_FOR_LOOP for (int i = 0; i < CloverTerm._grid->oSites(); i++) { + T._odata[i]()(0, 0) = timesMinusI(F._odata[i]()()); T._odata[i]()(1, 1) = timesI(F._odata[i]()()); T._odata[i]()(2, 2) = timesMinusI(F._odata[i]()()); diff --git a/lib/qcd/utils/WilsonLoops.h b/lib/qcd/utils/WilsonLoops.h index fe813298..86609ffc 100644 --- a/lib/qcd/utils/WilsonLoops.h +++ b/lib/qcd/utils/WilsonLoops.h @@ -327,7 +327,7 @@ static void StapleMult(GaugeMat &staple, const GaugeLorentz &Umu, int mu) { static void FieldStrength(GaugeMat &FS, const GaugeLorentz &Umu, int mu, int nu){ // Fmn +--<--+ Ut +--<--+ // | | | | - // (x)+-->--+ +-->--+(x) + // (x)+-->--+ +-->--+(x) - h.c. // | | | | // +--<--+ +--<--+ diff --git a/tests/forces/Test_wilsonclover_force.cc b/tests/forces/Test_wilsonclover_force.cc index c99cfa98..82adb8ab 100644 --- a/tests/forces/Test_wilsonclover_force.cc +++ b/tests/forces/Test_wilsonclover_force.cc @@ -1,6 +1,6 @@ /************************************************************************************* - Grid physics library, www.github.com/paboyle/Grid + Grid physics library, www.github.com/paboyle/Grid Source file: ./tests/Test_wilson_force.cc @@ -45,14 +45,17 @@ int main(int argc, char **argv) int threads = GridThread::GetThreads(); std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl; - std::vector seeds({1, 2, 3, 4}); + 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); gaussian(pRNG, phi); @@ -61,16 +64,53 @@ int main(int argc, char **argv) LatticeGaugeField U(&Grid); - SU3::HotConfiguration(pRNG,U); +/* + 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); + //////////////////////////////////// // Unmodified matrix element //////////////////////////////////// RealD mass = -4.0; //kills the diagonal term Real csw = 1.0; WilsonCloverFermionR Dw(U, Grid, RBGrid, mass, csw); - + Dw.ImportGauge(U); Dw.M(phi, Mphi); ComplexD S = innerProduct(Mphi, Mphi); // Action : pdag MdagM p @@ -78,11 +118,23 @@ int main(int argc, char **argv) LatticeGaugeField UdSdU(&Grid); LatticeGaugeField tmp(&Grid); - Dw.MDeriv(tmp, Mphi, phi, DaggerNo); UdSdU = tmp; - Dw.MDeriv(tmp, phi, Mphi, DaggerYes); UdSdU += tmp; + //////////////////////////////////////////// + Dw.MDeriv(tmp, Mphi, phi, DaggerNo); + UdSdU = tmp; + Dw.MDeriv(tmp, phi, Mphi, DaggerYes); + UdSdU += tmp; + ///////////////////////////////////////////// + // Take the traceless antihermitian component - UdSdU = Ta(UdSdU); - + //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 //////////////////////////////////// @@ -90,28 +142,63 @@ int main(int argc, char **argv) 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++) { // 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++) + 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++) + { + 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); + } + } 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); @@ -143,16 +230,14 @@ int main(int argc, char **argv) dSmom2 = zero; - // need for this??? - // ultimately it is just a 2.0 factor in UdSdU for (int mu = 0; mu < Nd; mu++) { - mommu = PeekIndex(UdSdU, mu); // P_mu = + mommu = PeekIndex(UdSdU, mu); // P_mu = mommu = Ta(mommu) * 2.0; // Mom = (P_mu - P_mu^dag) - trace(P_mu - P_mu^dag) PokeIndex(UdSdU, mommu, mu); // UdSdU_mu = Mom } - std::cout << GridLogMessage<< "Antihermiticity tests - 2 " << std::endl; + std::cout << GridLogMessage << "Antihermiticity tests - 2 " << std::endl; for (int mu = 0; mu < Nd; mu++) { mommu = PeekIndex(mom, mu); @@ -167,7 +252,6 @@ int main(int argc, char **argv) } ///////////////////////////////////////////////////// - for (int mu = 0; mu < Nd; mu++) { forcemu = PeekIndex(UdSdU, mu);