Added scalar action phi^4

Check Norm2 output (Complex type assumption)

lib/qcd/QCD.h

@@ -506,6 +506,9 @@ namespace QCD {
 #include <Grid/qcd/utils/SUnTwoIndex.h>
 #include <Grid/qcd/representations/hmc_types.h>
 
+// Scalar field
+#include <Grid/qcd/utils/ScalarObjs.h>
+
 #include <Grid/qcd/action/Actions.h>
 
 #include <Grid/qcd/smearing/Smearing.h>

lib/qcd/action/ActionBase.h

@@ -81,13 +81,13 @@ struct ActionLevel {
 };
 */
 
-template <class GaugeField, class Repr = NoHirep >
+template <class Field, class Repr = NoHirep >
 struct ActionLevel {
  public:
   unsigned int multiplier; 
 
   // Fundamental repr actions separated because of the smearing
-  typedef Action<GaugeField>* ActPtr;
+  typedef Action<Field>* ActPtr;
 
   // construct a tuple of vectors of the actions for the corresponding higher
   // representation fields
@@ -97,9 +97,6 @@ struct ActionLevel {
 
   std::vector<ActPtr>& actions;
 
-  // Temporary conversion between ActionLevel and ActionLevelHirep
-  //ActionLevelHirep(ActionLevel<GaugeField>& AL ):actions(AL.actions), multiplier(AL.multiplier){}
-
   ActionLevel(unsigned int mul = 1) : actions(std::get<0>(actions_hirep)), multiplier(mul) {
     // initialize the hirep vectors to zero.
     //apply(this->resize, actions_hirep, 0); //need a working resize
@@ -110,10 +107,10 @@ struct ActionLevel {
 
 
 
-  template < class Field >
-  void push_back(Action<Field>* ptr) {
+  template < class GenField >
+  void push_back(Action<GenField>* ptr) {
     // insert only in the correct vector
-    std::get< Index < Field, action_hirep_types>::value >(actions_hirep).push_back(ptr);
+    std::get< Index < GenField, action_hirep_types>::value >(actions_hirep).push_back(ptr);
   };
 
  
@@ -142,9 +139,6 @@ struct ActionLevel {
 };
 
 
-//template <class GaugeField>
-//using ActionSet = std::vector<ActionLevel<GaugeField> >;
-
 template <class GaugeField, class R>
 using ActionSet = std::vector<ActionLevel<GaugeField, R> >;
 

lib/qcd/action/Actions.h

@@ -60,6 +60,12 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #include <Grid/qcd/action/gauge/WilsonGaugeAction.h>
 #include <Grid/qcd/action/gauge/PlaqPlusRectangleAction.h>
 
+////////////////////////////////////////////
+// Scalar Actions
+////////////////////////////////////////////
+#include <Grid/qcd/action/scalar/scalarImpl.h>
+#include <Grid/qcd/action/scalar/ScalarAction.h>
+
 namespace Grid {
 namespace QCD {
 
@@ -90,6 +96,12 @@ typedef SymanzikGaugeAction<ConjugateGimplR>        ConjugateSymanzikGaugeAction
 typedef SymanzikGaugeAction<ConjugateGimplF>        ConjugateSymanzikGaugeActionF;
 typedef SymanzikGaugeAction<ConjugateGimplD>        ConjugateSymanzikGaugeActionD;
 
+
+  typedef ScalarAction<ScalarImplR>                 ScalarActionR;
+  typedef ScalarAction<ScalarImplF>                 ScalarActionF;
+  typedef ScalarAction<ScalarImplD>                 ScalarActionD;
+  
+  
 }}
 
 ////////////////////////////////////////////////////////////////////////////////////////////////////

lib/qcd/action/gauge/GaugeImpl.h

@@ -47,6 +47,8 @@ template <class Gimpl> class WilsonLoops;
   typedef typename GImpl::SiteLink SiteGaugeLink;
 
 #define INHERIT_FIELD_TYPES(Impl) \
+  typedef typename Impl::Simd Simd;           \
+  typedef typename Impl::SiteField SiteField; \
   typedef typename Impl::Field Field; 
 
 // hard codes the exponential approximation in the template
@@ -88,7 +90,9 @@ public:
     }
   }
 
-   
+   static inline Field projectForce(Field& P){
+    return Ta(P);
+   }
   
   static inline void update_field(Field& P, Field& U, double ep){
     for (int mu = 0; mu < Nd; mu++) {
@@ -117,7 +121,7 @@ public:
   static inline void TepidConfiguration(GridParallelRNG &pRNG, Field &U) {
     SU<Nc>::TepidConfiguration(pRNG, U);
   }
-  
+
   static inline void ColdConfiguration(GridParallelRNG &pRNG, Field &U) {
     SU<Nc>::ColdConfiguration(pRNG, U);
   }

lib/qcd/action/scalar/ScalarAction.h

@@ -0,0 +1,77 @@
+/*************************************************************************************
+
+  Grid physics library, www.github.com/paboyle/Grid
+
+  Source file: ./lib/qcd/action/gauge/WilsonGaugeAction.h
+
+  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
+  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 SCALAR_ACTION_H
+#define SCALAR_ACTION_H
+
+namespace Grid {
+namespace QCD {
+
+////////////////////////////////////////////////////////////////////////
+// Wilson Gauge Action .. should I template the Nc etc..
+////////////////////////////////////////////////////////////////////////
+
+template <class Impl>
+class ScalarAction : public Action<typename Impl::Field> {
+ public:
+  INHERIT_FIELD_TYPES(Impl);
+
+ private:
+  RealD mass_square;
+  RealD lambda;
+
+ public:
+  ScalarAction(RealD ms, RealD l) : mass_square(ms), lambda(l){};
+
+  virtual void refresh(const Field &U,
+                       GridParallelRNG &pRNG){};  // noop as no pseudoferms
+
+  virtual RealD S(const Field &p) {
+    return (mass_square * 0.5 + Nd) * ScalarObs<Impl>::sumphisquared(p) +
+           (lambda / 24.) * ScalarObs<Impl>::sumphifourth(p) +
+           ScalarObs<Impl>::sumphider(p);
+  };
+
+  virtual void deriv(const Field &p,
+                     Field &force) {
+    Field tmp(p._grid);
+    Field p2(p._grid);
+    ScalarObs<Impl>::phisquared(p2, p);
+    tmp = -(Cshift(p, 0, -1) + Cshift(p, 0, 1));
+    for (int mu = 1; mu < Nd; mu++) tmp -= Cshift(p, mu, -1) + Cshift(p, mu, 1);
+
+    force=+(mass_square + 2. * Nd) * p + (lambda / 6.) * p2 * p + tmp;
+  };
+};
+}
+}
+
+#endif

lib/qcd/action/scalar/scalarImpl.h

@@ -0,0 +1,56 @@
+#ifndef SCALAR_IMPL
+#define SCALAR_IMPL
+
+
+namespace Grid {
+namespace QCD {
+
+template <class S>
+class ScalarImplTypes {
+ public:
+  typedef S Simd;
+
+  template <typename vtype>
+  using iImplField = iScalar<iScalar<iScalar<vtype> > >;
+
+  typedef iImplField<Simd> SiteField;
+
+  
+  typedef Lattice<SiteField> Field;
+
+  static inline void generate_momenta(Field& P, GridParallelRNG& pRNG){
+  	gaussian(pRNG, P);
+  }
+
+  static inline Field projectForce(Field& P){return P;}
+
+  static inline void update_field(Field& P, Field& U, double ep){
+  	U += P*ep;
+  }
+
+  static inline RealD FieldSquareNorm(Field& U){
+  	return (- sum(trace(U*U))/2.0);
+  }
+
+  static inline void HotConfiguration(GridParallelRNG &pRNG, Field &U) {
+    gaussian(pRNG, U);
+  }
+
+  static inline void TepidConfiguration(GridParallelRNG &pRNG, Field &U) {
+    gaussian(pRNG, U);
+  }
+
+  static inline void ColdConfiguration(GridParallelRNG &pRNG, Field &U) {
+    U = 1.0;
+  }
+
+};
+
+
+typedef ScalarImplTypes<vReal> ScalarImplR;
+typedef ScalarImplTypes<vRealF> ScalarImplF;
+typedef ScalarImplTypes<vRealD> ScalarImplD;
+
+}}
+
+#endif

lib/qcd/hmc/GenericHMCrunner.h

@@ -34,16 +34,21 @@ namespace QCD {
 
 // Virtual Class for HMC specific for gauge theories
 // implement a specific theory by defining the BuildTheAction
-template <class Implementation, class RepresentationsPolicy = NoHirep>
+template <class Implementation,
+          class IOCheckpointer = BinaryHmcCheckpointer<Implementation>,
+          class RepresentationsPolicy = NoHirep>
 class BinaryHmcRunnerTemplate {
  public:
   INHERIT_FIELD_TYPES(Implementation);
+  typedef Implementation ImplPolicy;
 
   enum StartType_t { ColdStart, HotStart, TepidStart, CheckpointStart };
 
   ActionSet<Field, RepresentationsPolicy> TheAction;
   // Add here a vector of HmcObservable 
   // that can be injected from outside  
+  std::vector< HmcObservable<typename Implementation::Field> > ObservablesList;
+
 
   GridCartesian *UGrid;
   GridCartesian *FGrid;
@@ -119,8 +124,7 @@ class BinaryHmcRunnerTemplate {
     std::string format = std::string("IEEE64BIG");
     std::string conf_prefix = std::string("ckpoint_lat");
     std::string rng_prefix = std::string("ckpoint_rng");
-    BinaryHmcCheckpointer<Implementation> Checkpoint(conf_prefix, rng_prefix,
-                                                     SaveInterval, format);
+    IOCheckpointer Checkpoint(conf_prefix, rng_prefix, SaveInterval, format);
 
     HMCparameters HMCpar;
     HMCpar.StartTrajectory = StartTraj;
@@ -154,20 +158,31 @@ class BinaryHmcRunnerTemplate {
     SmearingPolicy.set_Field(U);
 
     HybridMonteCarlo<IntegratorType> HMC(HMCpar, MDynamics, sRNG, pRNG, U);
-    HMC.AddObservable(&Checkpoint);
+    //HMC.AddObservable(&Checkpoint);
+
+    for (int obs = 0; obs < ObservablesList.size(); obs++)
+      HMC.AddObservable(&ObservablesList[obs]); 
 
     // Run it
     HMC.evolve();
   }
 };
 
+// These are for gauge fields
 typedef BinaryHmcRunnerTemplate<PeriodicGimplR> BinaryHmcRunner;
 typedef BinaryHmcRunnerTemplate<PeriodicGimplF> BinaryHmcRunnerF;
 typedef BinaryHmcRunnerTemplate<PeriodicGimplD> BinaryHmcRunnerD;
 
+typedef BinaryHmcRunnerTemplate<PeriodicGimplR, NerscHmcCheckpointer<PeriodicGimplR> > NerscTestHmcRunner;
+
+
 template <class RepresentationsPolicy>
 using BinaryHmcRunnerTemplateHirep =
     BinaryHmcRunnerTemplate<PeriodicGimplR, RepresentationsPolicy>;
+
+
+
+typedef BinaryHmcRunnerTemplate<ScalarImplR, BinaryHmcCheckpointer<ScalarImplR>, ScalarFields> ScalarBinaryHmcRunner;
 }
 }
 #endif

lib/qcd/hmc/NerscCheckpointer.h

@@ -36,7 +36,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 namespace Grid{
   namespace QCD{
     
-    // Only for gauge fields
+    // Only for Gauge fields
     template<class Gimpl> 
     class NerscHmcCheckpointer : public HmcObservable<typename Gimpl::GaugeField> {
     private:
@@ -44,12 +44,13 @@ namespace Grid{
       std::string rngStem;
       int SaveInterval;
     public:
-      INHERIT_GIMPL_TYPES(Gimpl);
+      INHERIT_GIMPL_TYPES(Gimpl);// 
 
-      NerscHmcCheckpointer(std::string cf, std::string rn,int savemodulo) {
+      NerscHmcCheckpointer(std::string cf, std::string rn,int savemodulo, std::string format = "") {
         configStem  = cf;
         rngStem     = rn;
         SaveInterval= savemodulo;
+        // format is fixed to IEEE64BIG for NERSC
       };
 
       void TrajectoryComplete(int traj, GaugeField &U, GridSerialRNG &sRNG, GridParallelRNG & pRNG )

lib/qcd/hmc/integrators/Integrator.h

@@ -128,7 +128,7 @@ class Integrator {
                 << "Smearing (on/off): " << as[level].actions.at(a)->is_smeared
                 << std::endl;
       if (as[level].actions.at(a)->is_smeared) Smearer.smeared_force(force);
-      force = Ta(force);
+      force = FieldImplementation::projectForce(force); // Ta for gauge fields
       std::cout << GridLogIntegrator
                 << "Force average: " << norm2(force) / (U._grid->gSites())
                 << std::endl;
@@ -201,9 +201,10 @@ class Integrator {
 
   // Initialization of momenta and actions
   void refresh(Field& U, GridParallelRNG& pRNG) {
+    assert(P._grid == U._grid);
     std::cout << GridLogIntegrator << "Integrator refresh\n";
     FieldImplementation::generate_momenta(P, pRNG);
-
+ 
     // Update the smeared fields, can be implemented as observer
     // necessary to keep the fields updated even after a reject
     // of the Metropolis

lib/qcd/representations/fundamental.h

@@ -34,8 +34,21 @@ class FundamentalRep {
 
 };
 
+  template<class Field> 
+  class EmptyRep {
+  public:
+    typedef Field LatticeField;
+
+    explicit EmptyRep(GridBase* grid) {} //do nothing
+    void update_representation(const LatticeField& Uin) {} // do nothing
+    LatticeField RtoFundamentalProject(const LatticeField& in, Real scale = 1.0) const{}// do nothing
+  };
+    
+
+  
 typedef	 FundamentalRep<Nc> FundamentalRepresentation;
 
+  
 }
 }
 

lib/qcd/representations/hmc_types.h

@@ -4,6 +4,8 @@
 #include <Grid/qcd/representations/adjoint.h>
 #include <Grid/qcd/representations/two_index.h>
 #include <Grid/qcd/representations/fundamental.h>
+#include <Grid/qcd/action/scalar/scalarImpl.h>
+
 #include <tuple>
 #include <utility>
 
@@ -39,13 +41,17 @@ class Representations {
   int size() { return tuple_size; }
 
   // update the fields
+  // fields in the main representation always the first in the list
+  // get the field type
+  typedef typename std::tuple_element<0,Representation_Fields>::type LatticeSourceField;
+  
   template <std::size_t I = 0>
   inline typename std::enable_if<(I == tuple_size), void>::type update(
-      LatticeGaugeField& U) {}
+      LatticeSourceField& U) {}
 
   template <std::size_t I = 0>
   inline typename std::enable_if<(I < tuple_size), void>::type update(
-      LatticeGaugeField& U) {
+      LatticeSourceField& U) {
     std::get<I>(rep).update_representation(U);
     update<I + 1>(U);
   }
@@ -55,6 +61,7 @@ class Representations {
 };
 
 typedef Representations<FundamentalRepresentation> NoHirep;
+typedef Representations<EmptyRep<typename ScalarImplR::Field> > ScalarFields;
 
 // Helper classes to access the elements
 // Strips the first N parameters from the tuple

lib/qcd/smearing/GaugeConfiguration.h

@@ -11,24 +11,23 @@ namespace Grid {
 namespace QCD {
 
   //trivial class for no smearing
-  template< class Gimpl >
+  template< class Impl >
 class NoSmearing {
 public:
-  INHERIT_GIMPL_TYPES(Gimpl);
+  INHERIT_FIELD_TYPES(Impl);
 
-  GaugeField*
-      ThinLinks;
+  Field* ThinField;
 
-  NoSmearing(): ThinLinks(NULL) {}
+  NoSmearing(): ThinField(NULL) {}
 
-  void set_Field(GaugeField& U) { ThinLinks = &U; }
+  void set_Field(Field& U) { ThinField = &U; }
 
-  void smeared_force(GaugeField& SigmaTilde) const {}
+  void smeared_force(Field&) const {}
 
-  GaugeField& get_SmearedU() { return *ThinLinks; }
+  Field& get_SmearedU() { return *ThinField; }
 
-  GaugeField& get_U(bool smeared = false) {
-    return *ThinLinks;
+  Field& get_U(bool smeared = false) {
+    return *ThinField;
   }
 
 };

lib/qcd/utils/ScalarObjs.h

@@ -0,0 +1,96 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid
+
+    Source file: ./lib/qcd/utils/WilsonLoops.h
+
+    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
+    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 SCALAR_OBJS_H
+#define SCALAR_OBJS_H
+namespace Grid {
+namespace QCD {
+
+// Scalar field obs
+template <class Impl>
+class ScalarObs {
+ public:
+  //////////////////////////////////////////////////
+  // squared field
+  //////////////////////////////////////////////////
+
+  static void phisquared(typename Impl::Field &fsq,
+                         const typename Impl::Field &f) {
+    fsq = f * f;
+  }
+  //////////////////////////////////////////////////
+  // phi^4 interaction term
+  //////////////////////////////////////////////////
+
+  static void phifourth(typename Impl::Field &fsq,
+                        const typename Impl::Field &f) {
+    fsq = f * f * f * f;
+  }
+
+  //////////////////////////////////////////////////
+  // phi(x)phi(x+mu)
+  //////////////////////////////////////////////////
+
+  static void phider(typename Impl::Field &fsq,
+                     const typename Impl::Field &f) {
+    fsq = Cshift(f, 0, -1) * f;
+    for (int mu = 1; mu < Nd; mu++) fsq += Cshift(f, mu, -1) * f;
+  }
+
+  //////////////////////////////////////////////////
+  // Vol sum of the previous obs.
+  //////////////////////////////////////////////////
+
+  static RealD sumphider(const typename Impl::Field &f) {
+    typename Impl::Field tmp(f._grid);
+    tmp = Cshift(f, 0, -1) * f;
+    for (int mu = 1; mu < Nd; mu++) {
+      tmp += Cshift(f, mu, -1) * f;
+    }
+    return -sum(trace(tmp));
+  }
+
+  static RealD sumphisquared(const typename Impl::Field &f) {
+    typename Impl::Field tmp(f._grid);
+    tmp = f * f;
+    return sum(trace(tmp));
+  }
+
+  static RealD sumphifourth(const typename Impl::Field &f) {
+    typename Impl::Field tmp(f._grid);
+    phifourth(tmp, f);
+    return sum(trace(tmp));
+  }
+};
+}
+}
+
+#endif