Adding some documentation for HMC

lib/qcd/action/Actions.h

@@ -72,6 +72,10 @@ namespace QCD {
 typedef WilsonGaugeAction<PeriodicGimplR>          WilsonGaugeActionR;
 typedef WilsonGaugeAction<PeriodicGimplF>          WilsonGaugeActionF;
 typedef WilsonGaugeAction<PeriodicGimplD>          WilsonGaugeActionD;
+typedef VariableWilsonGaugeAction<PeriodicGimplR>          VariableWilsonGaugeActionR;
+typedef VariableWilsonGaugeAction<PeriodicGimplF>          VariableWilsonGaugeActionF;
+typedef VariableWilsonGaugeAction<PeriodicGimplD>          VariableWilsonGaugeActionD;
+
 typedef PlaqPlusRectangleAction<PeriodicGimplR>    PlaqPlusRectangleActionR;
 typedef PlaqPlusRectangleAction<PeriodicGimplF>    PlaqPlusRectangleActionF;
 typedef PlaqPlusRectangleAction<PeriodicGimplD>    PlaqPlusRectangleActionD;

lib/qcd/action/gauge/WilsonGaugeAction.h

@@ -1,86 +1,260 @@
-    /*************************************************************************************
+/*************************************************************************************
 
-    Grid physics library, www.github.com/paboyle/Grid 
+Grid physics library, www.github.com/paboyle/Grid
 
-    Source file: ./lib/qcd/action/gauge/WilsonGaugeAction.h
+Source file: ./lib/qcd/action/gauge/WilsonGaugeAction.h
 
-    Copyright (C) 2015
+Copyright (C) 2015
 
 Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 Author: neo <cossu@post.kek.jp>
 Author: paboyle <paboyle@ph.ed.ac.uk>
 
-    This program is free software; you can redistribute it and/or modify
+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
+it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
+the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
+(at your option) any later version.
 
-    This program is distributed in the hope that it will be useful,
+This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
+but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
+GNU General Public License for more details.
 
-    You should have received a copy of the GNU General Public License along
+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.,
+with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 
-    See the full license in the file "LICENSE" in the top level distribution directory
+See the full license in the file "LICENSE" in the top level distribution
-    *************************************************************************************/
+directory
-    /*  END LEGAL */
+*************************************************************************************/
+/*  END LEGAL */
 #ifndef QCD_WILSON_GAUGE_ACTION_H
 #define QCD_WILSON_GAUGE_ACTION_H
 
-namespace Grid{
+namespace Grid {
-  namespace QCD{
+namespace QCD {
 
-    ////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////
-    // Wilson Gauge Action .. should I template the Nc etc..
+// Wilson Gauge Action .. should I template the Nc etc..
-    ////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////
-    template<class Gimpl>
+template <class Gimpl>
-    class WilsonGaugeAction : public Action<typename Gimpl::GaugeField> {
+class WilsonGaugeAction : public Action<typename Gimpl::GaugeField> {
-    public:
+ public:
+  INHERIT_GIMPL_TYPES(Gimpl);
 
-      INHERIT_GIMPL_TYPES(Gimpl);
+  //      typedef LorentzScalar<GaugeField> GaugeLinkField;
 
-      //      typedef LorentzScalar<GaugeField> GaugeLinkField;
+ private:
+  RealD beta;
 
-    private:
+ public:
-      RealD beta;
+  WilsonGaugeAction(RealD b) : beta(b){};
-    public:
-    WilsonGaugeAction(RealD b):beta(b){};
 
-      virtual void refresh(const GaugeField &U, GridParallelRNG& pRNG) {}; // noop as no pseudoferms
+  virtual void refresh(const GaugeField &U,
+                       GridParallelRNG &pRNG){};  // noop as no pseudoferms
 
-      virtual RealD S(const GaugeField &U) {
+  virtual RealD S(const GaugeField &U) {
-	RealD plaq = WilsonLoops<Gimpl>::avgPlaquette(U);
+    RealD plaq = WilsonLoops<Gimpl>::avgPlaquette(U);
-	RealD vol = U._grid->gSites();
+    RealD vol = U._grid->gSites();
-	RealD action=beta*(1.0 -plaq)*(Nd*(Nd-1.0))*vol*0.5;
+    RealD action = beta * (1.0 - plaq) * (Nd * (Nd - 1.0)) * vol * 0.5;
-	return action;
+    return action;
-      };
+  };
 
-      virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
+  virtual void deriv(const GaugeField &U, GaugeField &dSdU) {
-	//not optimal implementation FIXME
+    // not optimal implementation FIXME
-	//extend Ta to include Lorentz indexes
+    // extend Ta to include Lorentz indexes
 
-	//RealD factor = 0.5*beta/RealD(Nc);
+    // RealD factor = 0.5*beta/RealD(Nc);
-	RealD factor = 0.5*beta/RealD(Nc);
+    RealD factor = 0.5 * beta / RealD(Nc);
 
-	GaugeLinkField Umu(U._grid);
+    GaugeLinkField Umu(U._grid);
-	GaugeLinkField dSdU_mu(U._grid);
+    GaugeLinkField dSdU_mu(U._grid);
-	for (int mu=0; mu < Nd; mu++){
+    for (int mu = 0; mu < Nd; mu++) {
+      Umu = PeekIndex<LorentzIndex>(U, mu);
 
-	  Umu = PeekIndex<LorentzIndex>(U,mu);
+      // Staple in direction mu
+      WilsonLoops<Gimpl>::Staple(dSdU_mu, U, mu);
+      dSdU_mu = Ta(Umu * dSdU_mu) * factor;
       
-	  // Staple in direction mu
+      PokeIndex<LorentzIndex>(dSdU, dSdU_mu, mu);
-	  WilsonLoops<Gimpl>::Staple(dSdU_mu,U,mu);
+    }
-	  dSdU_mu = Ta(Umu*dSdU_mu)*factor;
+  };
-	  PokeIndex<LorentzIndex>(dSdU, dSdU_mu, mu);
+};
-	}
-      };
-    };
 
+template <class Gimpl>
+class VariableWilsonGaugeAction : public Action<typename Gimpl::GaugeField> {
+ public:
+  INHERIT_GIMPL_TYPES(Gimpl);
+
+ private:
+  std::vector<RealD> b_bulk;// bulk couplings
+  std::vector<RealD> b_xdim;//extra dimension couplings
+  GridBase *grid;
+  LatticeComplex beta_xdim;
+  LatticeComplex beta_xdim_shifted;
+  LatticeComplex beta_bulk;
+
+ public:
+  VariableWilsonGaugeAction(std::vector<RealD> bulk, std::vector<RealD> xdim,
+                            GridBase *_grid)
+      : b_bulk(bulk),
+        b_xdim(xdim),
+        grid(_grid),
+        beta_xdim(grid),
+        beta_xdim_shifted(grid),
+        beta_bulk(grid)
+        {
+          //check that the grid is ok
+          //todo
+          int Ndim = Nd;//change later
+
+          LatticeComplex temp(grid);
+
+          Lattice<iScalar<vInteger> > coor(grid);
+
+          LatticeCoordinate(coor, Ndim - 1);
+
+          int Nex = grid->_fdimensions[Ndim - 1];
+          assert(b_bulk.size() == Nex);
+          assert(b_xdim.size() == Nex);
+
+          beta_xdim = zero;
+          for (int tau = 0; tau < Nex; tau++) {
+            temp = b_xdim[tau];
+            beta_xdim = where(coor == tau, temp, beta_xdim);
+          }
+
+          beta_xdim_shifted = Cshift(beta_xdim, Ndim - 1, -1);
+
+          beta_bulk = zero;
+          for (int tau = 0; tau < Nex; tau++) {
+            temp = b_bulk[tau];
+            beta_bulk = where(coor == tau, temp, beta_bulk);
+          }
+        };
+
+  virtual void refresh(const GaugeField &U,
+                       GridParallelRNG &pRNG){};  // noop as no pseudoferms
+
+  virtual RealD S(const GaugeField &Umu) {
+    int Ndim = Nd; // change later for generality
+    conformable(grid, Umu._grid);
+
+    std::vector<GaugeLinkField> U(Ndim, grid);
+
+    for (int mu = 0; mu < Ndim; mu++) {
+      U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
+    }
+
+    LatticeComplex dirPlaq(grid);
+    LatticeComplex Plaq(grid);
+
+    RealD OneOnNc = 1.0 / Real(Nc);
+
+    /////////////
+    // Lower dim plaquettes
+    /////////////
+    Plaq = zero;
+    for (int mu = 1; mu < Ndim - 1; mu++) {
+      for (int nu = 0; nu < mu; nu++) {
+        WilsonLoops<Gimpl>::traceDirPlaquette(dirPlaq, U, mu, nu);
+        Plaq = Plaq + (1.0 - dirPlaq * OneOnNc) * beta_bulk;
+      }
+    }
+
+    /////////////
+    // Extra dimension
+    /////////////
+    {
+      int mu = Ndim - 1;
+      for (int nu = 0; nu < mu; nu++) {
+        WilsonLoops<Gimpl>::traceDirPlaquette(dirPlaq, U, mu, nu);
+        Plaq = Plaq + (1.0 - dirPlaq * OneOnNc) * beta_xdim;
+      }
+    }
+
+    TComplex Tp = sum(Plaq);
+    Complex p = TensorRemove(Tp);
+    RealD action = p.real();
+    return action;
+  };
+
+  virtual void deriv(const GaugeField &U, GaugeField &dSdU) {
+    // not optimal implementation FIXME
+    // extend Ta to include Lorentz indexes
+
+    // for the higher dimension plaquettes take the upper plaq of the
+    // 4d slice and multiply by beta[s] and the lower and multiply by beta[s-1]
+    // todo
+    // derivative of links mu = 0, ... Nd-1 inside plaq (mu,5)
+    // for these I need upper and lower staples separated
+    // each multiplied with their own beta
+    // derivative of links mu = 5
+    // living on the same slice, share the same beta
+
+
+conformable(grid,U._grid);
+int Ndim = Nd;  // change later
+RealD factor = 0.5 / RealD(Nc);
+
+GaugeLinkField Umu(grid);
+GaugeLinkField dSdU_mu(grid);
+GaugeLinkField staple(grid);
+
+for (int mu = 0; mu < Ndim; mu++) {
+  Umu = PeekIndex<LorentzIndex>(U, mu);
+  dSdU_mu = zero;
+
+  for (int nu = 0; nu < Ndim; nu++) {
+    if (nu != mu) {
+      if ((mu < (Ndim - 1)) && (nu < (Ndim - 1))) {
+        // Spacelike case apply beta space
+        WilsonLoops<Gimpl>::Staple(staple, U, mu, nu);
+        staple = staple * beta_bulk;
+        dSdU_mu += staple;
+
+      } else if (mu == (Ndim - 1)) {
+        // nu space; mu time link
+        assert(nu < (Ndim - 1));
+        assert(mu == (Ndim - 1));
+
+        // mu==tau dir link deriv, nu spatial
+        WilsonLoops<Gimpl>::Staple(staple, U, mu, nu);
+        staple = staple * beta_xdim;
+        dSdU_mu += staple;
+
+      } else {
+        assert(mu < (Ndim - 1));
+        assert(nu == (Ndim - 1));
+
+        // nu time; mu space link
+
+        // staple forwards in tau
+        WilsonLoops<Gimpl>::StapleUpper(staple, U, mu, nu);
+        staple = staple * beta_xdim;
+        dSdU_mu += staple;
+
+        // staple backwards in tau
+        WilsonLoops<Gimpl>::StapleLower(staple, U, mu, nu);
+        staple = staple * beta_xdim_shifted;
+        dSdU_mu += staple;
+      }
+    }
   }
+
+  dSdU_mu = Ta(Umu * dSdU_mu) * factor;
+  PokeIndex<LorentzIndex>(dSdU, dSdU_mu, mu);
+  }
+
+
+
+  };
+};
+
+
+
+}
 }
 
 #endif

lib/qcd/hmc/UsingHMC.md

@@ -0,0 +1,107 @@
+Using HMC in Grid version 0.5.1
+
+These are the instructions to use the Generalised HMC on Grid version 0.5.1.
+Disclaimer: GRID is still under active development so any information here can be changed in future releases.
+
+
+Command line options
+===================
+(relevant file GenericHMCrunner.h)
+The initial configuration can be changed at the command line using 
+--StartType <your choice>
+valid choices, one among these
+HotStart, ColdStart, TepidStart, CheckpointStart
+default: HotStart
+
+example
+./My_hmc_exec  --StartType HotStart
+
+The CheckpointStart option uses the prefix for the configurations and rng seed files defined in your executable and the initial configuration is specified by
+--StartTrajectory <integer>
+default: 0
+
+The number of trajectories for a specific run are specified at command line by
+--Trajectories <integer>
+default: 1
+
+The number of thermalization steps (i.e. steps when the Metropolis acceptance check is turned off) is specified by
+--Thermalizations <integer>
+default: 10
+
+
+Any other parameter is defined in the source for the executable.
+
+HMC controls
+===========
+
+The lines 
+
+  std::vector<int> SerSeed({1, 2, 3, 4, 5});
+  std::vector<int> ParSeed({6, 7, 8, 9, 10});
+
+define the seeds for the serial and the parallel RNG.
+
+The line 
+
+  TheHMC.MDparameters.set(20, 1.0);// MDsteps, traj length
+
+declares the number of molecular dynamics steps and the total trajectory length.
+
+
+Actions
+======
+
+Action names are defined in the file
+lib/qcd/Actions.h
+
+Gauge actions list:
+
+WilsonGaugeActionR;
+WilsonGaugeActionF;
+WilsonGaugeActionD;
+PlaqPlusRectangleActionR;
+PlaqPlusRectangleActionF;
+PlaqPlusRectangleActionD;
+IwasakiGaugeActionR;
+IwasakiGaugeActionF;
+IwasakiGaugeActionD;
+SymanzikGaugeActionR;
+SymanzikGaugeActionF;
+SymanzikGaugeActionD;
+
+
+ConjugateWilsonGaugeActionR;
+ConjugateWilsonGaugeActionF;
+ConjugateWilsonGaugeActionD;
+ConjugatePlaqPlusRectangleActionR;
+ConjugatePlaqPlusRectangleActionF;
+ConjugatePlaqPlusRectangleActionD;
+ConjugateIwasakiGaugeActionR;
+ConjugateIwasakiGaugeActionF;
+ConjugateIwasakiGaugeActionD;
+ConjugateSymanzikGaugeActionR;
+ConjugateSymanzikGaugeActionF;
+ConjugateSymanzikGaugeActionD;
+
+
+ScalarActionR;
+ScalarActionF;
+ScalarActionD;
+
+
+each of these action accept one single parameter at creation time (beta).
+Example for creating a Symanzik action with beta=4.0
+
+	SymanzikGaugeActionR(4.0)
+
+The suffixes R,F,D in the action names refer to the Real
+(the precision is defined at compile time by the --enable-precision flag in the configure),
+Float and Double, that force the precision of the action to be 32, 64 bit respectively.
+
+
+
+
+
+
+
+

lib/qcd/utils/WilsonLoops.h

@@ -86,8 +86,7 @@ public:
   // sum over all x,y,z,t and over all planes of plaquette
   //////////////////////////////////////////////////
   static RealD sumPlaquette(const GaugeLorentz &Umu) {
-    std::vector<GaugeMat> U(4, Umu._grid);
+    std::vector<GaugeMat> U(Nd, Umu._grid);
-
     for (int mu = 0; mu < Nd; mu++) {
       U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
     }
@@ -95,11 +94,12 @@ public:
     LatticeComplex Plaq(Umu._grid);
 
     sitePlaquette(Plaq, U);
-
     TComplex Tp = sum(Plaq);
     Complex p = TensorRemove(Tp);
     return p.real();
   }
+
+
   //////////////////////////////////////////////////
   // average over all x,y,z,t and over all planes of plaquette
   //////////////////////////////////////////////////
@@ -114,7 +114,7 @@ public:
   // average over traced single links
   //////////////////////////////////////////////////
   static RealD linkTrace(const GaugeLorentz &Umu) {
-    std::vector<GaugeMat> U(4, Umu._grid);
+    std::vector<GaugeMat> U(Nd, Umu._grid);
 
     LatticeComplex Tr(Umu._grid);
     Tr = zero;
@@ -139,7 +139,7 @@ public:
 
     GridBase *grid = Umu._grid;
 
-    std::vector<GaugeMat> U(4, grid);
+    std::vector<GaugeMat> U(Nd, grid);
     for (int d = 0; d < Nd; d++) {
       U[d] = PeekIndex<LorentzIndex>(Umu, d);
     }
@@ -233,7 +233,7 @@ public:
     if (nu != mu) {
       GridBase *grid = Umu._grid;
 
-      std::vector<GaugeMat> U(4, grid);
+      std::vector<GaugeMat> U(Nd, grid);
       for (int d = 0; d < Nd; d++) {
         U[d] = PeekIndex<LorentzIndex>(Umu, d);
       }
@@ -256,6 +256,37 @@ public:
     }
   }
 
+  //////////////////////////////////////////////////
+  // 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) {
+    staple = zero;
+
+    if (nu != mu) {
+      GridBase *grid = Umu._grid;
+
+      std::vector<GaugeMat> U(Nd, grid);
+      for (int d = 0; d < Nd; d++) {
+        U[d] = PeekIndex<LorentzIndex>(Umu, d);
+      }
+
+      // mu
+      // ^
+      // |__>  nu
+
+      //  __
+      // |
+      // |__
+      //
+      //
+      staple += Gimpl::ShiftStaple(
+          Gimpl::CovShiftBackward(U[nu], nu,
+                                  Gimpl::CovShiftBackward(U[mu], mu, U[nu])),
+          mu);
+    }
+  }
+
   //////////////////////////////////////////////////////
   // Similar to above for rectangle is required
   //////////////////////////////////////////////////////
@@ -375,8 +406,8 @@ public:
         //             |___ ___|
         //
 
-        //	tmp= Staple2x1* Cshift(U[mu],mu,-2);
+        //  tmp= Staple2x1* Cshift(U[mu],mu,-2);
-        //	Stap+= Cshift(tmp,mu,1) ;
+        //  Stap+= Cshift(tmp,mu,1) ;
         Stap += Cshift(Staple2x1, mu, 1) * Cshift(U[mu], mu, -1);
         ;
 

tests/hmc/Make.inc

@@ -1,5 +1,5 @@
-tests: Test_hmc_EODWFRatio_Binary Test_hmc_EODWFRatio Test_hmc_EODWFRatio_Gparity Test_hmc_EOWilsonFermionGauge Test_hmc_EOWilsonRatio Test_hmc_GparityIwasakiGauge Test_hmc_GparityWilsonGauge Test_hmc_IwasakiGauge Test_hmc_RectGauge Test_hmc_ScalarAction Test_hmc_WilsonAdjointFermionGauge Test_hmc_WilsonFermionGauge_Binary Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_hmc_WilsonMixedRepresentationsFermionGauge Test_hmc_WilsonRatio Test_hmc_WilsonTwoIndexSymmetricFermionGauge Test_multishift_sqrt Test_remez Test_rhmc_EOWilson1p1 Test_rhmc_EOWilsonRatio Test_rhmc_Wilson1p1 Test_rhmc_WilsonRatio
+tests: Test_hmc_EODWFRatio_Binary Test_hmc_EODWFRatio Test_hmc_EODWFRatio_Gparity Test_hmc_EOWilsonFermionGauge Test_hmc_EOWilsonRatio Test_hmc_GparityIwasakiGauge Test_hmc_GparityWilsonGauge Test_hmc_IwasakiGauge Test_hmc_RectGauge Test_hmc_ScalarAction Test_hmc_WilsonAdjointFermionGauge Test_hmc_WilsonFermionGauge_Binary Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge_Binary Test_hmc_WilsonGauge Test_hmc_WilsonMixedRepresentationsFermionGauge Test_hmc_WilsonRatio Test_hmc_WilsonTwoIndexSymmetricFermionGauge Test_multishift_sqrt Test_remez Test_rhmc_EOWilson1p1 Test_rhmc_EOWilsonRatio Test_rhmc_Wilson1p1 Test_rhmc_WilsonRatio
-EXTRA_PROGRAMS = Test_hmc_EODWFRatio_Binary Test_hmc_EODWFRatio Test_hmc_EODWFRatio_Gparity Test_hmc_EOWilsonFermionGauge Test_hmc_EOWilsonRatio Test_hmc_GparityIwasakiGauge Test_hmc_GparityWilsonGauge Test_hmc_IwasakiGauge Test_hmc_RectGauge Test_hmc_ScalarAction Test_hmc_WilsonAdjointFermionGauge Test_hmc_WilsonFermionGauge_Binary Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge Test_hmc_WilsonMixedRepresentationsFermionGauge Test_hmc_WilsonRatio Test_hmc_WilsonTwoIndexSymmetricFermionGauge Test_multishift_sqrt Test_remez Test_rhmc_EOWilson1p1 Test_rhmc_EOWilsonRatio Test_rhmc_Wilson1p1 Test_rhmc_WilsonRatio
+EXTRA_PROGRAMS = Test_hmc_EODWFRatio_Binary Test_hmc_EODWFRatio Test_hmc_EODWFRatio_Gparity Test_hmc_EOWilsonFermionGauge Test_hmc_EOWilsonRatio Test_hmc_GparityIwasakiGauge Test_hmc_GparityWilsonGauge Test_hmc_IwasakiGauge Test_hmc_RectGauge Test_hmc_ScalarAction Test_hmc_WilsonAdjointFermionGauge Test_hmc_WilsonFermionGauge_Binary Test_hmc_WilsonFermionGauge Test_hmc_WilsonGauge_Binary Test_hmc_WilsonGauge Test_hmc_WilsonMixedRepresentationsFermionGauge Test_hmc_WilsonRatio Test_hmc_WilsonTwoIndexSymmetricFermionGauge Test_multishift_sqrt Test_remez Test_rhmc_EOWilson1p1 Test_rhmc_EOWilsonRatio Test_rhmc_Wilson1p1 Test_rhmc_WilsonRatio
 
 Test_hmc_EODWFRatio_Binary_SOURCES=Test_hmc_EODWFRatio_Binary.cc
 Test_hmc_EODWFRatio_Binary_LDADD=-lGrid
@@ -40,6 +40,9 @@ Test_hmc_WilsonFermionGauge_Binary_LDADD=-lGrid
 Test_hmc_WilsonFermionGauge_SOURCES=Test_hmc_WilsonFermionGauge.cc
 Test_hmc_WilsonFermionGauge_LDADD=-lGrid
 
+Test_hmc_WilsonGauge_Binary_SOURCES=Test_hmc_WilsonGauge_Binary.cc
+Test_hmc_WilsonGauge_Binary_LDADD=-lGrid
+
 Test_hmc_WilsonGauge_SOURCES=Test_hmc_WilsonGauge.cc
 Test_hmc_WilsonGauge_LDADD=-lGrid
 

tests/hmc/Test_hmc_EODWFRatio_Binary.cc

@@ -78,7 +78,8 @@ class HmcRunner : public BinaryHmcRunner {
     LatticeGaugeField  U(UGrid);
 
     // Gauge action
-    IwasakiGaugeActionR Iaction(4.0);
+    double beta = 4.0;
+    IwasakiGaugeActionR Iaction(beta);
 
     Real mass = 0.04;
     Real pv   = 1.0;
@@ -87,7 +88,9 @@ class HmcRunner : public BinaryHmcRunner {
     FermionAction DenOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,scale);
     FermionAction NumOp(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,pv,M5,scale);
   
-    ConjugateGradient<FermionField>  CG(1.0e-8,10000);
+    double StoppingCondition = 1.0e-8;
+    double MaxCGIterations = 10000;
+    ConjugateGradient<FermionField>  CG(StoppingCondition,MaxCGIterations);
     TwoFlavourEvenOddRatioPseudoFermionAction<ImplPolicy> Nf2(NumOp, DenOp,CG,CG);
   
     // Set smearing (true/false), default: false
@@ -98,6 +101,9 @@ class HmcRunner : public BinaryHmcRunner {
     ActionLevel<Field> Level1(1);
     Level1.push_back(&Nf2);
 
+    // this level will integrate with a
+    // step that is 4 times finer
+    // than the previous level
     ActionLevel<Field> Level2(4);
     Level2.push_back(&Iaction);
 

tests/hmc/Test_hmc_WilsonGauge_Binary.cc

@@ -0,0 +1,136 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./tests/Test_hmc_WilsonFermionGauge.cc
+
+Copyright (C) 2015
+
+Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
+Author: neo <cossu@post.kek.jp>
+Author: Guido Cossu <guido.cossu@ed.ac.uk>
+
+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 */
+#include <Grid/Grid.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+namespace Grid {
+namespace QCD {
+
+class HMCRunnerParameters : Serializable {
+ public:
+  GRID_SERIALIZABLE_CLASS_MEMBERS(HMCRunnerParameters,
+                                  double, beta,
+                                  double, mass,
+                                  int, MaxCGIterations,
+                                  double, StoppingCondition,
+                                  bool, smearedAction,
+                                  int, SaveInterval,
+                                  std::string, format,
+                                  std::string, conf_prefix,
+                                  std::string, rng_prefix,
+                                  double, rho,
+                                  int, SmearingLevels,
+                                  );
+
+  HMCRunnerParameters() {}
+};
+
+// Derive from the BinaryHmcRunner (templated for gauge fields)
+class HmcRunner : public BinaryHmcRunner {
+ public:
+  void BuildTheAction(int argc, char **argv)
+
+  {
+    typedef WilsonImplR ImplPolicy;
+    typedef WilsonFermionR FermionAction;
+    typedef typename FermionAction::FermionField FermionField;
+
+    UGrid = SpaceTimeGrid::makeFourDimGrid(
+        GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()),
+        GridDefaultMpi());
+    UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+
+    FGrid = UGrid;
+    FrbGrid = UrbGrid;
+
+    // temporarily need a gauge field
+    LatticeGaugeField U(UGrid);
+
+    // Gauge action
+    int Ls = UGrid->_gdimensions[Nd - 1];
+    std::vector<RealD> betat(Ls,5);
+    std::vector<RealD> betas(Ls);
+    betas={5,6,6,5};
+    std:cout << "Betas:" << betas << std::endl;
+    VariableWilsonGaugeActionR Waction(betas, betat, UGrid);
+    //WilsonGaugeActionR Waction(5.6);
+
+    // Collect actions
+    ActionLevel<Field> Level1(1);
+    Level1.push_back(&Waction);
+    TheAction.push_back(Level1);
+
+    // Add observables
+    int SaveInterval = 1;
+    std::string format = std::string("IEEE64BIG");
+    std::string conf_prefix = std::string("ckpoint_lat");
+    std::string rng_prefix = std::string("ckpoint_rng");
+    BinaryHmcCheckpointer<BinaryHmcRunner::ImplPolicy> Checkpoint(
+        conf_prefix, rng_prefix, SaveInterval, format);
+    // Can implement also a specific function in the hmcrunner
+    // AddCheckpoint (...) that takes the same parameters + a string/tag
+    // defining the type of the checkpointer
+    // with tags can be implemented by overloading and no ifs
+    // Then force all checkpoint to have few common functions
+    // return an object that is then passed to the Run function
+
+    PlaquetteLogger<BinaryHmcRunner::ImplPolicy> PlaqLog(
+        std::string("Plaquette"));
+    ObservablesList.push_back(&PlaqLog);
+    ObservablesList.push_back(&Checkpoint);
+
+    Run(argc, argv, Checkpoint);  // no smearing
+  };
+};
+}
+}
+
+int main(int argc, char **argv) {
+  Grid_init(&argc, &argv);
+
+  int threads = GridThread::GetThreads();
+  std::cout << GridLogMessage << "Grid is setup to use " << threads
+            << " threads" << std::endl;
+
+  HmcRunner TheHMC;
+
+  // Seeds for the random number generators
+  std::vector<int> SerSeed({1, 2, 3, 4, 5});
+  std::vector<int> ParSeed({6, 7, 8, 9, 10});
+  TheHMC.RNGSeeds(SerSeed, ParSeed);
+
+  TheHMC.MDparameters.set(20, 1.0);// MDsteps, traj length
+
+  TheHMC.BuildTheAction(argc, argv);
+}