Added topological charge measurement

lib/qcd/QCD.h

@@ -34,6 +34,10 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 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;

lib/qcd/modules/ObservableModules.h

@@ -94,6 +94,22 @@ class PlaquetteMod: public ObservableModule<PlaquetteLogger<Impl>, NoParameters>
   PlaquetteMod(): ObsBase(NoParameters()){}
 };
 
+template < class Impl >
+class TopologicalChargeMod: public ObservableModule<TopologicalCharge<Impl>, NoParameters>{
+  typedef ObservableModule<TopologicalCharge<Impl>, NoParameters> ObsBase;
+  using ObsBase::ObsBase; // for constructors
+
+
+
+  // acquire resource
+  virtual void initialize(){
+    this->ObservablePtr.reset(new TopologicalCharge<Impl>());
+  }
+  public:
+  TopologicalChargeMod(): ObsBase(NoParameters()){}
+};
+
+
 
 }// QCD temporarily here
 

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

lib/qcd/observables/plaquette.h

@@ -36,9 +36,6 @@ namespace QCD {
 // this is only defined for a gauge theory
 template <class Impl>
 class PlaquetteLogger : public HmcObservable<typename Impl::Field> {
- private:
-  std::string Stem;
-
  public:
   // here forces the Impl to be of gauge fields
   // if not the compiler will complain

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 Impl>
 class TopologicalCharge : public HmcObservable<typename Impl::Field> {
- 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<Impl>::FieldStrength(Bx, U, Ydir, Zdir);
+    WilsonLoops<Impl>::FieldStrength(By, U, Zdir, Xdir);
+    WilsonLoops<Impl>::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<Impl>::FieldStrength(Ex, U, Tdir, Xdir);
+    WilsonLoops<Impl>::FieldStrength(Ey, U, Tdir, Ydir);
+    WilsonLoops<Impl>::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<Real>::digits10 + 1)
+        << "Topological Charge: [ " << traj << " ] "<< q << std::endl;
+
+    std::cout.precision(def_prec);
+    }
 
 };
-
+}
-)
 }
 
 #endif  //  HMC_TOP_CHARGE_H

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])),
+                                    Gimpl::CovShiftBackward(U[mu], mu, U[nu])), mu);
-            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<LorentzIndex>(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
   //////////////////////////////////////////////////////

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<HMCWrapper::ImplPolicy> PlaqObs;
+  typedef TopologicalChargeMod<HMCWrapper::ImplPolicy> QObs;
   TheHMC.Resources.AddObservable<PlaqObs>();
+  TheHMC.Resources.AddObservable<QObs>();
   //////////////////////////////////////////////
 
   /////////////////////////////////////////////////////////////