From 3d0fe1537436cbd0e5545d5185c175e05e543693 Mon Sep 17 00:00:00 2001 From: Guido Cossu Date: Fri, 17 Mar 2017 16:14:57 +0900 Subject: [PATCH] Added topological charge measurement --- lib/qcd/QCD.h | 4 ++ lib/qcd/modules/ObservableModules.h | 16 ++++++++ lib/qcd/observables/hmc_observable.h | 1 + lib/qcd/observables/plaquette.h | 3 -- lib/qcd/observables/topological_charge.h | 47 +++++++++++++++++++++--- lib/qcd/utils/WilsonLoops.h | 32 +++++++++++----- tests/hmc/Test_hmc_WilsonGauge.cc | 2 + 7 files changed, 87 insertions(+), 18 deletions(-) diff --git a/lib/qcd/QCD.h b/lib/qcd/QCD.h index 8602fcb7..9d1bba0b 100644 --- a/lib/qcd/QCD.h +++ b/lib/qcd/QCD.h @@ -34,6 +34,10 @@ Author: paboyle namespace Grid{ namespace QCD { + static const int Xdir = 0; + static const int Ydir = 1; + static const int Zdir = 2; + static const int Tdir = 3; static const int Xp = 0; static const int Yp = 1; diff --git a/lib/qcd/modules/ObservableModules.h b/lib/qcd/modules/ObservableModules.h index 27efecc7..579fc1ec 100644 --- a/lib/qcd/modules/ObservableModules.h +++ b/lib/qcd/modules/ObservableModules.h @@ -94,6 +94,22 @@ class PlaquetteMod: public ObservableModule, NoParameters> PlaquetteMod(): ObsBase(NoParameters()){} }; +template < class Impl > +class TopologicalChargeMod: public ObservableModule, NoParameters>{ + typedef ObservableModule, NoParameters> ObsBase; + using ObsBase::ObsBase; // for constructors + + + + // acquire resource + virtual void initialize(){ + this->ObservablePtr.reset(new TopologicalCharge()); + } + public: + TopologicalChargeMod(): ObsBase(NoParameters()){} +}; + + }// QCD temporarily here diff --git a/lib/qcd/observables/hmc_observable.h b/lib/qcd/observables/hmc_observable.h index 446b3dd1..db629ce7 100644 --- a/lib/qcd/observables/hmc_observable.h +++ b/lib/qcd/observables/hmc_observable.h @@ -44,5 +44,6 @@ class HmcObservable { } // namespace Grid #include "plaquette.h" +#include "topological_charge.h" #endif // HMC_OBSERVABLE_H diff --git a/lib/qcd/observables/plaquette.h b/lib/qcd/observables/plaquette.h index 84e7be04..c848dcd9 100644 --- a/lib/qcd/observables/plaquette.h +++ b/lib/qcd/observables/plaquette.h @@ -36,9 +36,6 @@ namespace QCD { // this is only defined for a gauge theory template class PlaquetteLogger : public HmcObservable { - private: - std::string Stem; - public: // here forces the Impl to be of gauge fields // if not the compiler will complain diff --git a/lib/qcd/observables/topological_charge.h b/lib/qcd/observables/topological_charge.h index aa76be1b..6188b9e8 100644 --- a/lib/qcd/observables/topological_charge.h +++ b/lib/qcd/observables/topological_charge.h @@ -31,19 +31,54 @@ directory #define HMC_TOP_CHARGE_H namespace Grid { -namespace QCD ( +namespace QCD { // this is only defined for a gauge theory template class TopologicalCharge : public HmcObservable { - private: - std::string Stem; - public: + // here forces the Impl to be of gauge fields + // if not the compiler will complain + INHERIT_GIMPL_TYPES(Impl); + + // necessary for HmcObservable compatibility + typedef typename Impl::Field Field; + + void TrajectoryComplete(int traj, + Field &U, + GridSerialRNG &sRNG, + GridParallelRNG &pRNG) { + + // 4d topological charge + // Bx = -iF(y,z), By = -iF(z,y), Bz = -iF(x,y) + GaugeLinkField Bx(U._grid), By(U._grid), Bz(U._grid); + WilsonLoops::FieldStrength(Bx, U, Ydir, Zdir); + WilsonLoops::FieldStrength(By, U, Zdir, Xdir); + WilsonLoops::FieldStrength(Bz, U, Xdir, Ydir); + + // Ex = -iF(t,x), Ey = -iF(t,y), Ez = -iF(t,z) + GaugeLinkField Ex(U._grid), Ey(U._grid), Ez(U._grid); + WilsonLoops::FieldStrength(Ex, U, Tdir, Xdir); + WilsonLoops::FieldStrength(Ey, U, Tdir, Ydir); + WilsonLoops::FieldStrength(Ez, U, Tdir, Zdir); + + double coeff = 8.0/(32.0*M_PI*M_PI); + + LatticeComplex qfield = coeff*trace(Bx*Ex + By*Ey + Bz*Ez); + TComplex Tq = sum(qfield); + Real q = TensorRemove(Tq).real(); + + int def_prec = std::cout.precision(); + + std::cout << GridLogMessage + << std::setprecision(std::numeric_limits::digits10 + 1) + << "Topological Charge: [ " << traj << " ] "<< q << std::endl; + + std::cout.precision(def_prec); + } }; - -) +} } #endif // HMC_TOP_CHARGE_H diff --git a/lib/qcd/utils/WilsonLoops.h b/lib/qcd/utils/WilsonLoops.h index 5e0e2d5c..a610cfa4 100644 --- a/lib/qcd/utils/WilsonLoops.h +++ b/lib/qcd/utils/WilsonLoops.h @@ -228,8 +228,7 @@ public: // staple += Gimpl::ShiftStaple( Gimpl::CovShiftBackward(U[nu], nu, - Gimpl::CovShiftBackward(U[mu], mu, U[nu])), - mu); + Gimpl::CovShiftBackward(U[mu], mu, U[nu])), mu); } } } @@ -239,9 +238,6 @@ public: ////////////////////////////////////////////////// static void StapleUpper(GaugeMat &staple, const GaugeLorentz &Umu, int mu, int nu) { - - staple = zero; - if (nu != mu) { GridBase *grid = Umu._grid; @@ -259,7 +255,7 @@ public: // __| // - staple += Gimpl::ShiftStaple( + staple = Gimpl::ShiftStaple( Gimpl::CovShiftForward( U[nu], nu, Gimpl::CovShiftBackward( @@ -273,8 +269,6 @@ public: ////////////////////////////////////////////////// static void StapleLower(GaugeMat &staple, const GaugeLorentz &Umu, int mu, int nu) { - staple = zero; - if (nu != mu) { GridBase *grid = Umu._grid; @@ -292,13 +286,33 @@ public: // |__ // // - staple += Gimpl::ShiftStaple( + 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 = adj(Vup) - adj(Vdn); + GaugeMat u = PeekIndex(Umu, mu); // some redundant copies + GaugeMat vu = v*u; + FS = 0.25*Ta(u*v - Cshift(vu, mu, +1)); + } + + ////////////////////////////////////////////////////// // Similar to above for rectangle is required ////////////////////////////////////////////////////// diff --git a/tests/hmc/Test_hmc_WilsonGauge.cc b/tests/hmc/Test_hmc_WilsonGauge.cc index 337bbcaf..09b73cb8 100644 --- a/tests/hmc/Test_hmc_WilsonGauge.cc +++ b/tests/hmc/Test_hmc_WilsonGauge.cc @@ -66,7 +66,9 @@ int main(int argc, char **argv) { // Construct observables // here there is too much indirection typedef PlaquetteMod PlaqObs; + typedef TopologicalChargeMod QObs; TheHMC.Resources.AddObservable(); + TheHMC.Resources.AddObservable(); ////////////////////////////////////////////// /////////////////////////////////////////////////////////////