/************************************************************************************* Grid physics library, www.github.com/paboyle/Grid Source file: ./lib/qcd/utils/WilsonLoops.h Copyright (C) 2015 Author: Azusa Yamaguchi Author: Peter Boyle Author: neo Author: paboyle This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ #ifndef QCD_UTILS_WILSON_LOOPS_H #define QCD_UTILS_WILSON_LOOPS_H NAMESPACE_BEGIN(Grid); // Common wilson loop observables template class WilsonLoops : public Gimpl { public: INHERIT_GIMPL_TYPES(Gimpl); typedef typename Gimpl::GaugeLinkField GaugeMat; typedef typename Gimpl::GaugeField GaugeLorentz; ////////////////////////////////////////////////// // directed plaquette oriented in mu,nu plane ////////////////////////////////////////////////// static void dirPlaquette(GaugeMat &plaq, const std::vector &U, const int mu, const int nu) { // Annoyingly, must use either scope resolution to find dependent base // class, // or this-> ; there is no "this" in a static method. This forces explicit // Gimpl scope // resolution throughout the usage in this file, and rather defeats the // purpose of deriving // from Gimpl. /* plaq = Gimpl::CovShiftBackward( U[mu], mu, Gimpl::CovShiftBackward( U[nu], nu, Gimpl::CovShiftForward(U[mu], mu, U[nu]))); */ // _ //|< _| plaq = Gimpl::CovShiftForward(U[mu],mu, Gimpl::CovShiftForward(U[nu],nu, Gimpl::CovShiftBackward(U[mu],mu, Gimpl::CovShiftIdentityBackward(U[nu], nu)))); } ////////////////////////////////////////////////// // trace of directed plaquette oriented in mu,nu plane ////////////////////////////////////////////////// static void traceDirPlaquette(ComplexField &plaq, const std::vector &U, const int mu, const int nu) { GaugeMat sp(U[0]._grid); dirPlaquette(sp, U, mu, nu); plaq = trace(sp); } ////////////////////////////////////////////////// // sum over all planes of plaquette ////////////////////////////////////////////////// static void sitePlaquette(ComplexField &Plaq, const std::vector &U) { ComplexField sitePlaq(U[0]._grid); Plaq = zero; for (int mu = 1; mu < Nd; mu++) { for (int nu = 0; nu < mu; nu++) { traceDirPlaquette(sitePlaq, U, mu, nu); Plaq = Plaq + sitePlaq; } } } ////////////////////////////////////////////////// // sum over all x,y,z,t and over all planes of plaquette ////////////////////////////////////////////////// static RealD sumPlaquette(const GaugeLorentz &Umu) { std::vector U(Nd, Umu._grid); // inefficient here for (int mu = 0; mu < Nd; mu++) { U[mu] = PeekIndex(Umu, mu); } ComplexField Plaq(Umu._grid); sitePlaquette(Plaq, U); auto Tp = sum(Plaq); auto p = TensorRemove(Tp); return p.real(); } ////////////////////////////////////////////////// // average over all x,y,z,t and over all planes of plaquette ////////////////////////////////////////////////// static RealD avgPlaquette(const GaugeLorentz &Umu) { RealD sumplaq = sumPlaquette(Umu); double vol = Umu._grid->gSites(); double faces = (1.0 * Nd * (Nd - 1)) / 2.0; return sumplaq / vol / faces / Nc; // Nd , Nc dependent... FIXME } ////////////////////////////////////////////////// // average over traced single links ////////////////////////////////////////////////// static RealD linkTrace(const GaugeLorentz &Umu) { std::vector U(Nd, Umu._grid); ComplexField Tr(Umu._grid); Tr = zero; for (int mu = 0; mu < Nd; mu++) { U[mu] = PeekIndex(Umu, mu); Tr = Tr + trace(U[mu]); } auto Tp = sum(Tr); auto p = TensorRemove(Tp); double vol = Umu._grid->gSites(); return p.real() / vol / 4.0 / 3.0; }; ////////////////////////////////////////////////// // the sum over all staples on each site in direction mu,nu ////////////////////////////////////////////////// static void Staple(GaugeMat &staple, const GaugeLorentz &Umu, int mu, int nu) { GridBase *grid = Umu._grid; std::vector U(Nd, grid); for (int d = 0; d < Nd; d++) { U[d] = PeekIndex(Umu, d); } staple = zero; if (nu != mu) { // mu // ^ // |__> nu // __ // | // __| // staple += Gimpl::ShiftStaple( Gimpl::CovShiftForward( U[nu], nu, Gimpl::CovShiftBackward( U[mu], mu, Gimpl::CovShiftIdentityBackward(U[nu], nu))), mu); // __ // | // |__ // // staple += Gimpl::ShiftStaple( Gimpl::CovShiftBackward(U[nu], nu, Gimpl::CovShiftBackward(U[mu], mu, U[nu])), mu); } } // For the force term static void StapleMult(GaugeMat &staple, const GaugeLorentz &Umu, int mu) { GridBase *grid = Umu._grid; std::vector U(Nd, grid); for (int d = 0; d < Nd; d++) { // this operation is taking too much time U[d] = PeekIndex(Umu, d); } staple = zero; GaugeMat tmp1(grid); GaugeMat tmp2(grid); for (int nu = 0; nu < Nd; nu++) { if (nu != mu) { // this is ~10% faster than the Staple tmp1 = Cshift(U[nu], mu, 1); tmp2 = Cshift(U[mu], nu, 1); staple += tmp1* adj(U[nu]*tmp2); tmp2 = adj(U[mu]*tmp1)*U[nu]; staple += Cshift(tmp2, nu, -1); } } staple = U[mu]*staple; } ////////////////////////////////////////////////// // the sum over all staples on each site ////////////////////////////////////////////////// static void Staple(GaugeMat &staple, const GaugeLorentz &Umu, int mu) { GridBase *grid = Umu._grid; std::vector U(Nd, grid); for (int d = 0; d < Nd; d++) { U[d] = PeekIndex(Umu, d); } staple = zero; for (int nu = 0; nu < Nd; nu++) { if (nu != mu) { // mu // ^ // |__> nu // __ // | // __| // staple += Gimpl::ShiftStaple( Gimpl::CovShiftForward( U[nu], nu, Gimpl::CovShiftBackward( U[mu], mu, Gimpl::CovShiftIdentityBackward(U[nu], nu))), mu); // __ // | // |__ // // staple += Gimpl::ShiftStaple( Gimpl::CovShiftBackward(U[nu], nu, Gimpl::CovShiftBackward(U[mu], mu, U[nu])), mu); } } } ////////////////////////////////////////////////// // the sum over all staples on each site in direction mu,nu, upper part ////////////////////////////////////////////////// static void StapleUpper(GaugeMat &staple, const GaugeLorentz &Umu, int mu, int nu) { if (nu != mu) { GridBase *grid = Umu._grid; std::vector U(Nd, grid); for (int d = 0; d < Nd; d++) { U[d] = PeekIndex(Umu, d);// some redundant copies } // mu // ^ // |__> nu // __ // | // __| // staple = Gimpl::ShiftStaple( Gimpl::CovShiftForward( U[nu], nu, Gimpl::CovShiftBackward( U[mu], mu, Gimpl::CovShiftIdentityBackward(U[nu], nu))), mu); } } ////////////////////////////////////////////////// // the sum over all staples on each site in direction mu,nu, lower part ////////////////////////////////////////////////// static void StapleLower(GaugeMat &staple, const GaugeLorentz &Umu, int mu, int nu) { if (nu != mu) { GridBase *grid = Umu._grid; std::vector U(Nd, grid); for (int d = 0; d < Nd; d++) { U[d] = PeekIndex(Umu, d);// some redundant copies } // mu // ^ // |__> nu // __ // | // |__ // // staple = Gimpl::ShiftStaple( Gimpl::CovShiftBackward(U[nu], nu, Gimpl::CovShiftBackward(U[mu], mu, U[nu])), mu); } } ////////////////////////////////////////////////////// // Field Strength ////////////////////////////////////////////////////// static void FieldStrength(GaugeMat &FS, const GaugeLorentz &Umu, int mu, int nu){ // Fmn +--<--+ Ut +--<--+ // | | | | // (x)+-->--+ +-->--+(x) // | | | | // +--<--+ +--<--+ GaugeMat Vup(Umu._grid), Vdn(Umu._grid); StapleUpper(Vup, Umu, mu, nu); StapleLower(Vdn, Umu, mu, nu); GaugeMat v = Vup - Vdn; GaugeMat u = PeekIndex(Umu, mu); // some redundant copies GaugeMat vu = v*u; FS = 0.25*Ta(u*v + Cshift(vu, mu, -1)); } static Real TopologicalCharge(GaugeLorentz &U){ // 4d topological charge assert(Nd==4); // Bx = -iF(y,z), By = -iF(z,y), Bz = -iF(x,y) GaugeMat Bx(U._grid), By(U._grid), Bz(U._grid); FieldStrength(Bx, U, Ydir, Zdir); FieldStrength(By, U, Zdir, Xdir); FieldStrength(Bz, U, Xdir, Ydir); // Ex = -iF(t,x), Ey = -iF(t,y), Ez = -iF(t,z) GaugeMat Ex(U._grid), Ey(U._grid), Ez(U._grid); FieldStrength(Ex, U, Tdir, Xdir); FieldStrength(Ey, U, Tdir, Ydir); FieldStrength(Ez, U, Tdir, Zdir); double coeff = 8.0/(32.0*M_PI*M_PI); ComplexField qfield = coeff*trace(Bx*Ex + By*Ey + Bz*Ez); auto Tq = sum(qfield); return TensorRemove(Tq).real(); } ////////////////////////////////////////////////////// // Similar to above for rectangle is required ////////////////////////////////////////////////////// static void dirRectangle(GaugeMat &rect, const std::vector &U, const int mu, const int nu) { rect = Gimpl::CovShiftForward( U[mu], mu, Gimpl::CovShiftForward(U[mu], mu, U[nu])) * // ->->| adj(Gimpl::CovShiftForward( U[nu], nu, Gimpl::CovShiftForward(U[mu], mu, U[mu]))); rect = rect + Gimpl::CovShiftForward( U[mu], mu, Gimpl::CovShiftForward(U[nu], nu, U[nu])) * // ->|| adj(Gimpl::CovShiftForward( U[nu], nu, Gimpl::CovShiftForward(U[nu], nu, U[mu]))); } static void traceDirRectangle(ComplexField &rect, const std::vector &U, const int mu, const int nu) { GaugeMat sp(U[0]._grid); dirRectangle(sp, U, mu, nu); rect = trace(sp); } static void siteRectangle(ComplexField &Rect, const std::vector &U) { ComplexField siteRect(U[0]._grid); Rect = zero; for (int mu = 1; mu < Nd; mu++) { for (int nu = 0; nu < mu; nu++) { traceDirRectangle(siteRect, U, mu, nu); Rect = Rect + siteRect; } } } ////////////////////////////////////////////////// // sum over all x,y,z,t and over all planes of plaquette ////////////////////////////////////////////////// static RealD sumRectangle(const GaugeLorentz &Umu) { std::vector U(Nd, Umu._grid); for (int mu = 0; mu < Nd; mu++) { U[mu] = PeekIndex(Umu, mu); } ComplexField Rect(Umu._grid); siteRectangle(Rect, U); auto Tp = sum(Rect); auto p = TensorRemove(Tp); return p.real(); } ////////////////////////////////////////////////// // average over all x,y,z,t and over all planes of plaquette ////////////////////////////////////////////////// static RealD avgRectangle(const GaugeLorentz &Umu) { RealD sumrect = sumRectangle(Umu); double vol = Umu._grid->gSites(); double faces = (1.0 * Nd * (Nd - 1)); // 2 distinct orientations summed return sumrect / vol / faces / Nc; // Nd , Nc dependent... FIXME } ////////////////////////////////////////////////// // the sum over all staples on each site ////////////////////////////////////////////////// static void RectStapleDouble(GaugeMat &U2, const GaugeMat &U, int mu) { U2 = U * Cshift(U, mu, 1); } //////////////////////////////////////////////////////////////////////////// // Hop by two optimisation strategy does not work nicely with Gparity. (could // do, // but need to track two deep where cross boundary and apply a conjugation). // Must differentiate this in Gimpl, and use Gimpl::isPeriodicGaugeField to do // so . //////////////////////////////////////////////////////////////////////////// static void RectStapleOptimised(GaugeMat &Stap, std::vector &U2, std::vector &U, int mu) { Stap = zero; GridBase *grid = U[0]._grid; GaugeMat Staple2x1(grid); GaugeMat tmp(grid); for (int nu = 0; nu < Nd; nu++) { if (nu != mu) { // Up staple ___ ___ // | | tmp = Cshift(adj(U[nu]), nu, -1); tmp = adj(U2[mu]) * tmp; tmp = Cshift(tmp, mu, -2); Staple2x1 = Gimpl::CovShiftForward(U[nu], nu, tmp); // Down staple // |___ ___| // tmp = adj(U2[mu]) * U[nu]; Staple2x1 += Gimpl::CovShiftBackward(U[nu], nu, Cshift(tmp, mu, -2)); // ___ ___ // | ___| // |___ ___| // Stap += Cshift(Gimpl::CovShiftForward(U[mu], mu, Staple2x1), mu, 1); // ___ ___ // |___ | // |___ ___| // // tmp= Staple2x1* Cshift(U[mu],mu,-2); // Stap+= Cshift(tmp,mu,1) ; Stap += Cshift(Staple2x1, mu, 1) * Cshift(U[mu], mu, -1); ; // -- // | | // // | | tmp = Cshift(adj(U2[nu]), nu, -2); tmp = Gimpl::CovShiftBackward(U[mu], mu, tmp); tmp = U2[nu] * Cshift(tmp, nu, 2); Stap += Cshift(tmp, mu, 1); // | | // // | | // -- tmp = Gimpl::CovShiftBackward(U[mu], mu, U2[nu]); tmp = adj(U2[nu]) * tmp; tmp = Cshift(tmp, nu, -2); Stap += Cshift(tmp, mu, 1); } } } static void RectStaple(GaugeMat &Stap, const GaugeLorentz &Umu, int mu) { RectStapleUnoptimised(Stap, Umu, mu); } static void RectStaple(const GaugeLorentz &Umu, GaugeMat &Stap, std::vector &U2, std::vector &U, int mu) { if (Gimpl::isPeriodicGaugeField()) { RectStapleOptimised(Stap, U2, U, mu); } else { RectStapleUnoptimised(Stap, Umu, mu); } } static void RectStapleUnoptimised(GaugeMat &Stap, const GaugeLorentz &Umu, int mu) { GridBase *grid = Umu._grid; std::vector U(Nd, grid); for (int d = 0; d < Nd; d++) { U[d] = PeekIndex(Umu, d); } Stap = zero; for (int nu = 0; nu < Nd; nu++) { if (nu != mu) { // __ ___ // | __ | // Stap += Gimpl::ShiftStaple( Gimpl::CovShiftForward( U[mu], mu, Gimpl::CovShiftForward( U[nu], nu, Gimpl::CovShiftBackward( U[mu], mu, Gimpl::CovShiftBackward( U[mu], mu, Gimpl::CovShiftIdentityBackward(U[nu], nu))))), mu); // __ // |__ __ | Stap += Gimpl::ShiftStaple( Gimpl::CovShiftForward( U[mu], mu, Gimpl::CovShiftBackward( U[nu], nu, Gimpl::CovShiftBackward( U[mu], mu, Gimpl::CovShiftBackward(U[mu], mu, U[nu])))), mu); // __ // |__ __ | Stap += Gimpl::ShiftStaple( Gimpl::CovShiftBackward( U[nu], nu, Gimpl::CovShiftBackward( U[mu], mu, Gimpl::CovShiftBackward( U[mu], mu, Gimpl::CovShiftForward(U[nu], nu, U[mu])))), mu); // __ ___ // |__ | Stap += Gimpl::ShiftStaple( Gimpl::CovShiftForward( U[nu], nu, Gimpl::CovShiftBackward( U[mu], mu, Gimpl::CovShiftBackward( U[mu], mu, Gimpl::CovShiftBackward(U[nu], nu, U[mu])))), mu); // -- // | | // // | | Stap += Gimpl::ShiftStaple( Gimpl::CovShiftForward( U[nu], nu, Gimpl::CovShiftForward( U[nu], nu, Gimpl::CovShiftBackward( U[mu], mu, Gimpl::CovShiftBackward( U[nu], nu, Gimpl::CovShiftIdentityBackward(U[nu], nu))))), mu); // | | // // | | // -- Stap += Gimpl::ShiftStaple( Gimpl::CovShiftBackward( U[nu], nu, Gimpl::CovShiftBackward( U[nu], nu, Gimpl::CovShiftBackward( U[mu], mu, Gimpl::CovShiftForward(U[nu], nu, U[nu])))), mu); } } } }; typedef WilsonLoops ColourWilsonLoops; typedef WilsonLoops U1WilsonLoops; typedef WilsonLoops SU2WilsonLoops; typedef WilsonLoops SU3WilsonLoops; NAMESPACE_END(Grid); #endif