Added ForceGradient integrator. dH dropped so seems to work. Will only

believe it is right once I have pulled a dt^4 error scaling plot out.

TODO

@@ -2,6 +2,7 @@ TODO:
 
 => Clean up HMC
 - Reunitarise
+- Reversibility test.
 - Link smearing/boundary conds; Policy class based implementation
 
 - Integrators

lib/qcd/hmc/integrators/Integrator.h

@@ -76,28 +76,22 @@ namespace Grid{
 
       void update_P(GaugeField&U, int level,double ep){
 	t_P[level]+=ep;
-	for(int a=0; a<as[level].actions.size(); ++a){
+	update_P(P,U,level,ep);
-	  GaugeField force(U._grid);
-	  as[level].actions.at(a)->deriv(U,force);
-	  P = P - force*ep;
-	}
 
 	std::cout<<GridLogMessage;
 	for(int l=0; l<level;++l) std::cout<<"   ";	    
 	std::cout<<"["<<level<<"] P " << " dt "<< ep <<" : t_P "<< t_P[level] <<std::endl;
-
+      }
+      void update_P(GaugeField &Mom,GaugeField&U, int level,double ep){
+	for(int a=0; a<as[level].actions.size(); ++a){
+	  GaugeField force(U._grid);
+	  as[level].actions.at(a)->deriv(U,force);
+	  Mom = Mom - force*ep;
+	}
       }
 
       void update_U(GaugeField&U, double ep){
-	//rewrite exponential to deal automatically  with the lorentz index?
+	update_U(P,U,ep);
-	//	GaugeLinkField Umu(U._grid);
-	//	GaugeLinkField Pmu(U._grid);
-	for (int mu = 0; mu < Nd; mu++){
-	  auto Umu=PeekIndex<LorentzIndex>(U, mu);
-	  auto Pmu=PeekIndex<LorentzIndex>(P, mu);
-	  Umu = expMat(Pmu, ep, Params.Nexp)*Umu;
-	  PokeIndex<LorentzIndex>(U, Umu, mu);
-	}
 
 	t_U+=ep;
 	int fl = levels-1;
@@ -106,6 +100,17 @@ namespace Grid{
 	std::cout<<"["<<fl<<"] U " << " dt "<< ep <<" : t_U "<< t_U <<std::endl;
 
       }
+      void update_U(GaugeField &Mom, GaugeField&U, double ep){
+	//rewrite exponential to deal automatically  with the lorentz index?
+	//	GaugeLinkField Umu(U._grid);
+	//	GaugeLinkField Pmu(U._grid);
+	for (int mu = 0; mu < Nd; mu++){
+	  auto Umu=PeekIndex<LorentzIndex>(U, mu);
+	  auto Pmu=PeekIndex<LorentzIndex>(Mom, mu);
+	  Umu = expMat(Pmu, ep, Params.Nexp)*Umu;
+	  PokeIndex<LorentzIndex>(U, Umu, mu);
+	}
+      }
       
       /*
 	friend void Algorithm::step (GaugeField& U, 

lib/qcd/hmc/integrators/Integrator_algorithm.h

@@ -85,7 +85,7 @@ namespace Grid{
 	for(int l=0; l<=level; ++l) fin*= this->as[l].multiplier;
 	fin = fin-1;
 
-	double eps = this->Params.stepsize;
+	RealD eps = this->Params.stepsize;
 	for(int l=0; l<=level; ++l) eps/= this->as[l].multiplier;
 	
 	int multiplier = this->as[level].multiplier;
@@ -116,7 +116,7 @@ namespace Grid{
 
     template<class GaugeField> class MinimumNorm2 : public Integrator<GaugeField> {
     private:
-      const double lambda = 0.1931833275037836;
+      const RealD lambda = 0.1931833275037836;
 
     public:
 
@@ -132,7 +132,7 @@ namespace Grid{
 
 	int fl = this->as.size() -1;
 
-	double eps = this->Params.stepsize;
+	RealD eps = this->Params.stepsize;
 	
 	for(int l=0; l<=level; ++l) eps/= 2.0*this->as[l].multiplier;
 	
@@ -174,6 +174,88 @@ namespace Grid{
       }
     };
 
+
+    template<class GaugeField> class ForceGradient : public Integrator<GaugeField> {
+    private:
+      const RealD lambda = 1.0/6.0;;
+      const RealD chi    = 1.0/72.0;
+      const RealD xi     = 0.0;
+      const RealD theta  = 0.0;
+    public:
+
+    ForceGradient(GridBase* grid, 
+		  IntegratorParameters Par,
+		  ActionSet<GaugeField> & Aset): Integrator<GaugeField>(grid,Par,Aset) {};
+
+
+      void FG_update_P(GaugeField&U, int level,double fg_dt,double ep){
+	GaugeField Ufg(U._grid);
+	GaugeField Pfg(U._grid);
+	Ufg = U;
+	Pfg = zero;
+	std::cout << GridLogMessage << "FG update "<<fg_dt<<" "<<ep<<std::endl;
+	// prepare_fg; no prediction/result cache for now
+	// could relax CG stopping conditions for the 
+	// derivatives in the small step since the force gets multiplied by
+	// a tiny dt^2 term relative to main force.
+	//
+	// Presently 4 force evals, and should have 3, so 1.33x too expensive.
+	// could reduce this with sloppy CG to perhaps 1.15x too expensive
+	// even without prediction.
+	this->update_P(Pfg,Ufg,level,1.0); 
+	this->update_U(Pfg,Ufg,fg_dt);
+	this->update_P(Ufg,level,ep);
+      }
+
+      void step (GaugeField& U, int level, std::vector<int>& clock){
+
+	RealD eps = this->Params.stepsize;
+	for(int l=0; l<=level; ++l) eps/= 2.0*this->as[l].multiplier;
+
+	RealD Theta = theta*eps*eps*eps;
+	RealD Chi   = chi*eps*eps*eps;
+
+	int fl = this->as.size() -1;
+	
+	// which is final half step
+	int fin = this->as[0].multiplier;
+	for(int l=1; l<=level; ++l) fin*= 2.0*this->as[l].multiplier;
+	fin = 3*this->Params.MDsteps*fin -1;
+
+	int multiplier = this->as[level].multiplier;
+
+	for(int e=0; e<multiplier; ++e){       // steps per step
+
+	  int first_step,last_step;
+
+	  first_step = (clock[level]==0);
+
+	  if(first_step){    // initial half step 
+	    this->update_P(U,level,lambda*eps);   ++clock[level];
+	  }
+	  
+	  if(level == fl){          // lowest level 
+	    this->update_U(U,0.5*eps);
+	  }else{                 // recursive function call 
+	    this->step(U,level+1,clock);
+	  }
+	  
+	  this->FG_update_P(U,level,2*Chi/((1.0-2.0*lambda)*eps),(1.0-2.0*lambda)*eps); ++clock[level];
+	  
+	  if(level == fl){          // lowest level 
+	    this->update_U(U,0.5*eps);
+	  }else{                 // recursive function call 
+	    this->step(U,level+1,clock);
+	  }    
+	  
+	  last_step  = (clock[level]==fin);
+	  int mm = (last_step) ? 1 : 2;
+	  this->update_P(U,level,lambda*eps*mm); clock[level]+=mm;
+
+	}
+      }
+    };
+
   }
 }
 

tests/Test_hmc_EOWilsonFermionGauge.cc

@@ -49,9 +49,10 @@ int main (int argc, char ** argv)
   FullSet.push_back(Level2);
 
   // Create integrator
-  //  typedef LeapFrog  IntegratorAlgorithm;// change here to change the algorithm
+  //typedef LeapFrog<LatticeGaugeField>  IntegratorAlgorithm;// change here to change the algorithm
-  //  IntegratorParameters MDpar(12,40,1.0);
+  //typedef MinimumNorm2<LatticeGaugeField>  IntegratorAlgorithm;// change here to change the algorithm
-  typedef MinimumNorm2<LatticeGaugeField>  IntegratorAlgorithm;// change here to change the algorithm
+  typedef ForceGradient<LatticeGaugeField>  IntegratorAlgorithm;// change here to change the algorithm
+
   IntegratorParameters MDpar(10);
   IntegratorAlgorithm MDynamics(&Fine,MDpar, FullSet);