Reunitarise. Complete the HMC and integrator changes.

TODO

@@ -1,45 +1,44 @@
+RECENT
+---------------
+
+  - Clean up HMC                                                             -- DONE
+  - LorentzScalar<GaugeField> gets Gauge link type (cleaner).                -- DONE
+  - Simplified the integrators a bit.                                        -- DONE
+  - Multi-timescale looks broken and operating on single timescale for now.  -- DONE
+  - pass GaugeField as template param.                        -- DONE
+  - Reunitarise                                               -- DONE
+  - Force Gradient                                            -- DONE
+  - Prefer "RefreshInternal" or such like to "init" in naming -- DONE
+  - Parallel io improvements                                  -- DONE
+  - Plaquette and link trace checks into nersc reader from the Grid_nersc_io.cc test. -- DONE
+
 TODO:
-
+---------------
-=> Clean up HMC
+Policies:
-- Reunitarise
+* Link smearing/boundary conds; Policy class based implementation ; framework more in place
-- Reversibility test.
-- Link smearing/boundary conds; Policy class based implementation
-
-- Integrators
-  - Force Gradient              
-  - Sign of force term.
-  - Prefer "RefreshInternal" or such like to "init" in naming
-  - Rename "Ta" as too unclear
-  DONE: -- pass GaugeField as template param.
-        -- LorentzScalar<GaugeField> gets Gauge link type (cleaner).
-        -- Simplified the integrators a bit.
-        - Multi-timescale looks broken and operating on single timescale for now.
-          FIXED
-
 * Support different boundary conditions (finite temp, chem. potential ... )
 * Support different fermion representations? 
   - contained entirely within the integrator presently
+  - Sign of force term.
+  - Rename "Ta" as too unclear
 
-=> Lanczos
+- Reversibility test.
 
- - Rectangle gauge actions.
+- Lanczos
+
+- Rectangle gauge actions.
   Iwasaki,
   Symanzik,
   ... etc...
 
 - Prepare multigrid for HMC.
 - Alternate setup schemes.
-
+* Support for ILDG --- ugly, not done
-* Parallel io improvements DONE
-  - move Plaquette and link trace checks into nersc reader from the Grid_nersc_io.cc test. DONE
-  * Support for ILDG --- ugly, not done
-
 * Flavour matrices?
-
 * FFTnD ?
 
 ================================================================
-*** Hacks and bug fixes to clean up and Audits
+* Hacks and bug fixes to clean up and Audits
 ================================================================
 
 *  Extract/merge/set cleanup ; too many variants; rationalise and call simpler ones
@@ -65,7 +64,6 @@ Insert/Extract
 
 * Thread scaling tests Xeon, XeonPhi
 
-
 Not sure of status of this -- reverify. Things are working nicely now though.
 
 * Make the Tensor types and Complex etc... play more nicely.
@@ -84,6 +82,10 @@ Not sure of status of this -- reverify. Things are working nicely now though.
     template specialize the scalar scalar scalar sum and SliceSum,  on the basis of being
     pure scalar.
 
+======================================================================
+======================================================================
+======================================================================
+======================================================================
 Done: Cayley, Partial , ContFrac force terms.
 
 DONE

lib/lattice/Lattice_ET.h

@@ -178,6 +178,7 @@ GridUnopClass(UnaryConj,conjugate(a));
 GridUnopClass(UnaryTrace,trace(a));
 GridUnopClass(UnaryTranspose,transpose(a));
 GridUnopClass(UnaryTa,Ta(a));
+GridUnopClass(UnaryProjectOnGroup,ProjectOnGroup(a));
 GridUnopClass(UnaryReal,real(a));
 GridUnopClass(UnaryImag,imag(a));
 GridUnopClass(UnaryToReal,toReal(a));
@@ -290,6 +291,7 @@ GRID_DEF_UNOP(conjugate,UnaryConj);
 GRID_DEF_UNOP(trace,UnaryTrace);
 GRID_DEF_UNOP(transpose,UnaryTranspose);
 GRID_DEF_UNOP(Ta,UnaryTa);
+GRID_DEF_UNOP(ProjectOnGroup,UnaryProjectOnGroup);
 GRID_DEF_UNOP(real,UnaryReal);
 GRID_DEF_UNOP(imag,UnaryImag);
 GRID_DEF_UNOP(toReal,UnaryToReal);

lib/qcd/hmc/integrators/Integrator.h

@@ -73,7 +73,6 @@ namespace Grid{
       //      void register_observers();
       //      void notify_observers();
 
-
       void update_P(GaugeField&U, int level,double ep){
 	t_P[level]+=ep;
 	update_P(P,U,level,ep);
@@ -82,6 +81,7 @@ namespace Grid{
 	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);
@@ -108,18 +108,12 @@ namespace Grid{
 	  auto Umu=PeekIndex<LorentzIndex>(U, mu);
 	  auto Pmu=PeekIndex<LorentzIndex>(Mom, mu);
 	  Umu = expMat(Pmu, ep, Params.Nexp)*Umu;
+	  ProjectOnGroup(Umu);
 	  PokeIndex<LorentzIndex>(U, Umu, mu);
 	}
       }
       
-      /*
+      virtual void step (GaugeField& U,int level, int first,int last)=0;
-	friend void Algorithm::step (GaugeField& U, 
-				   int level, 
-				   std::vector<int>& clock,
-				   Integrator<GaugeField,Algorithm>* Integ);
-      */
-
-      virtual void step (GaugeField& U,int level, std::vector<int>& clock)=0;
 
     public:
 
@@ -178,23 +172,28 @@ namespace Grid{
 
       void integrate(GaugeField& U){
 
-	std::vector<int> clock;
+	// reset the clocks
+	t_U=0;
+	for(int level=0; level<as.size(); ++level){
+	  t_P[level]=0;
+	}	
 
-	clock.resize(as.size(),0);
-
-	// All the clock stuff is removed if we pass first, last to the step down the way
 	for(int step=0; step< Params.MDsteps; ++step){   // MD step
 	  int first_step = (step==0);
 	  int  last_step = (step==Params.MDsteps-1);
-	  this->step(U,0,clock);
+	  this->step(U,0,first_step,last_step);
 	}
 
-	// Check the clocks all match
+	// Check the clocks all match on all levels
 	for(int level=0; level<as.size(); ++level){
 	  assert(fabs(t_U - t_P[level])<1.0e-6); // must be the same
 	  std::cout<<GridLogMessage<<" times["<<level<<"]= "<<t_P[level]<< " " << t_U <<std::endl;
 	}	
 
+	// and that we indeed got to the end of the trajectory
+	assert(fabs(t_U-Params.trajL) < 1.0e-6);
+
+
       }
     };
     

lib/qcd/hmc/integrators/Integrator_algorithm.h

@@ -68,12 +68,12 @@ namespace Grid{
 
       typedef LeapFrog<GaugeField> Algorithm;
 
-    LeapFrog(GridBase* grid, 
+      LeapFrog(GridBase* grid, 
-	     IntegratorParameters Par,
+	       IntegratorParameters Par,
-	     ActionSet<GaugeField> & Aset): Integrator<GaugeField>(grid,Par,Aset) {};
+	       ActionSet<GaugeField> & Aset): Integrator<GaugeField>(grid,Par,Aset) {};
 
 
-      void step (GaugeField& U, int level, std::vector<int>& clock){
+      void step (GaugeField& U, int level,int _first, int _last){
 
 	int fl = this->as.size() -1;
 	// level  : current level
@@ -81,34 +81,27 @@ namespace Grid{
 	// eps    : current step size
 	
 	// Get current level step size
-	int fin = 2*this->Params.MDsteps;
+        RealD eps = this->Params.stepsize;
-	for(int l=0; l<=level; ++l) fin*= this->as[l].multiplier;
-	fin = fin-1;
-
-	RealD eps = this->Params.stepsize;
 	for(int l=0; l<=level; ++l) eps/= this->as[l].multiplier;
 	
 	int multiplier = this->as[level].multiplier;
 	for(int e=0; e<multiplier; ++e){
 
-	  int first_step,last_step;
+	  int first_step = _first && (e==0);
-
+	  int last_step  = _last  && (e==multiplier-1);
-	  first_step = (clock[level]==0);
 
 	  if(first_step){    // initial half step
-	    this->update_P(U, level,eps/2.0);	    ++clock[level];
+	    this->update_P(U, level,eps/2.0);
 	  }
 
 	  if(level == fl){          // lowest level
 	    this->update_U(U, eps);
 	  }else{                 // recursive function call
-	    this->step(U, level+1,clock);
+	    this->step(U, level+1,first_step,last_step);
 	  }
 
-	  last_step  = (clock[level]==fin);
 	  int mm = last_step ? 1 : 2;
 	  this->update_P(U, level,mm*eps/2.0);	    
-	  clock[level]+=mm;
 
 	}
       }
@@ -124,7 +117,7 @@ namespace Grid{
 		   IntegratorParameters Par,
 		   ActionSet<GaugeField> & Aset): Integrator<GaugeField>(grid,Par,Aset) {};
 
-      void step (GaugeField& U, int level, std::vector<int>& clock){
+      void step (GaugeField& U, int level, int _first,int _last){
 
 	// level  : current level
 	// fl     : final level
@@ -132,43 +125,38 @@ namespace Grid{
 
 	int fl = this->as.size() -1;
 
-	RealD eps = this->Params.stepsize;
+	RealD eps = this->Params.stepsize*2.0;                              
-	
+	for(int l=0; l<=level; ++l) eps/= 2.0*this->as[l].multiplier;   
-	for(int l=0; l<=level; ++l) eps/= 2.0*this->as[l].multiplier;
+
-	
+	// Nesting:  2xupdate_U of size eps/2
-	// which is final half step
+	// Next level is eps/2/multiplier
-	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;
+	  int first_step = _first && (e==0);
-
+	  int last_step  = _last  && (e==multiplier-1);
-	  first_step = (clock[level]==0);
 
 	  if(first_step){    // initial half step 
-	    this->update_P(U,level,lambda*eps);   ++clock[level];
+	    this->update_P(U,level,lambda*eps);
 	  }
 	  
 	  if(level == fl){          // lowest level 
 	    this->update_U(U,0.5*eps);
 	  }else{                 // recursive function call 
-	    this->step(U,level+1,clock);
+	    this->step(U,level+1,first_step,0);
 	  }
 	  
-	  this->update_P(U,level,(1.0-2.0*lambda)*eps); ++clock[level];
+	  this->update_P(U,level,(1.0-2.0*lambda)*eps);
 	  
 	  if(level == fl){          // lowest level 
 	    this->update_U(U,0.5*eps);
 	  }else{                 // recursive function call 
-	    this->step(U,level+1,clock);
+	    this->step(U,level+1,0,last_step);
 	  }    
 	  
-	  last_step  = (clock[level]==fin);
 	  int mm = (last_step) ? 1 : 2;
-	  this->update_P(U,level,lambda*eps*mm); clock[level]+=mm;
+	  this->update_P(U,level,lambda*eps*mm);
 
 	}
       }
@@ -207,50 +195,43 @@ namespace Grid{
 	this->update_P(Ufg,level,ep);
       }
 
-      void step (GaugeField& U, int level, std::vector<int>& clock){
+      void step (GaugeField& U, int level, int _first,int _last){
 
-	RealD eps = this->Params.stepsize;
+	RealD eps = this->Params.stepsize*2.0;                              
 	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);
+	  int first_step = _first && (e==0);
+	  int last_step  = _last  && (e==multiplier-1);
 
 	  if(first_step){    // initial half step 
-	    this->update_P(U,level,lambda*eps);   ++clock[level];
+	    this->update_P(U,level,lambda*eps);
 	  }
 	  
 	  if(level == fl){          // lowest level 
 	    this->update_U(U,0.5*eps);
 	  }else{                 // recursive function call 
-	    this->step(U,level+1,clock);
+	    this->step(U,level+1,first_step,0);
 	  }
 	  
-	  this->FG_update_P(U,level,2*Chi/((1.0-2.0*lambda)*eps),(1.0-2.0*lambda)*eps); ++clock[level];
+	  this->FG_update_P(U,level,2*Chi/((1.0-2.0*lambda)*eps),(1.0-2.0*lambda)*eps);
 	  
 	  if(level == fl){          // lowest level 
 	    this->update_U(U,0.5*eps);
 	  }else{                 // recursive function call 
-	    this->step(U,level+1,clock);
+	    this->step(U,level+1,0,last_step);
 	  }    
 	  
-	  last_step  = (clock[level]==fin);
 	  int mm = (last_step) ? 1 : 2;
-	  this->update_P(U,level,lambda*eps*mm); clock[level]+=mm;
+	  this->update_P(U,level,lambda*eps*mm); 
 
 	}
       }

scripts/hmc.sh

@@ -7,10 +7,14 @@ echo
 echo $SWEEPS thermalised sweeps
 echo
 plaq=`grep Plaq $LOG | tail -n $SWEEPS | awk '{ S=S+$10} END { print S/NR} ' `
-echo Plaquette: $plaq
+plaqe=`grep Plaq $LOG | tail -n $SWEEPS | awk '{ S=S+$10 ; SS=SS+$10*$10 } END { print sqrt( (SS/NR - S*S/NR/NR)/NR) } ' `
+echo "Plaquette: $plaq (${plaqe})"
 echo
+
+dHv=`grep dH $LOG | tail -n $SWEEPS | awk '{ S=S+$10 ; SS=SS+$10*$10 } END { print sqrt(SS/NR) } ' `
 edH=`grep dH $LOG | tail -n $SWEEPS | awk '{ S=S+exp(-$10)} END { print S/NR} '`
 echo "<e-dH>: $edH"
+echo "<rms dH>: $dHv"
 
 TRAJ=`grep Acc $LOG | wc -l`
 ACC=`grep Acc $LOG | grep ACCE | wc -l`
@@ -20,7 +24,7 @@ echo "Acceptance $PACC %  $ACC / $TRAJ "
 
 grep Plaq $LOG | awk '{ print $10 }' | uniq > plaq.dat
 grep dH $LOG | awk '{ print $10 }' > dH.dat
-echo set yrange [-1:1] > plot.gnu
+echo set yrange [-0.1:0.6] > plot.gnu
 echo set terminal 'pdf' >> plot.gnu
 echo "set output 'plaq.${LOG}.pdf'" >> plot.gnu
 echo "plot 'plaq.dat' w l, 'dH.dat' w l " >> plot.gnu

tests/Test_hmc_EOWilsonFermionGauge.cc

@@ -53,7 +53,7 @@ int main (int argc, char ** argv)
   //typedef MinimumNorm2<LatticeGaugeField>  IntegratorAlgorithm;// change here to change the algorithm
   typedef ForceGradient<LatticeGaugeField>  IntegratorAlgorithm;// change here to change the algorithm
 
-  IntegratorParameters MDpar(10);
+  IntegratorParameters MDpar(16);
   IntegratorAlgorithm MDynamics(&Fine,MDpar, FullSet);
 
   // Create HMC

tests/Test_quenched_update.cc

@@ -73,8 +73,10 @@ int main (int argc, char ** argv)
 	}
 
 	PokeIndex<LorentzIndex>(Umu,link,mu);
-	//reunitarise link;
 	
+	//reunitarise link;
+	ProjectOnGroup(Umu);
+
       }
 
     }