HMC for Wilson Gauge action works

Fixed bug in momenta generation

lib/qcd/action/gauge/WilsonGaugeAction.h

@@ -27,16 +27,13 @@ namespace Grid{
 	//not optimal implementation FIXME
 	//extend Ta to include Lorentz indexes
 	RealD factor = 0.5*beta/RealD(Nc);
-	std::cout << "Debug force Wilson  "<< factor << "\n";
 	MatrixField dSdU_mu(U._grid);
 	MatrixField Umu(U._grid);
 	for (int mu=0; mu < Nd; mu++){
 	  Umu = PeekIndex<LorentzIndex>(U,mu);
 	  // Staple in direction mu
 	  WilsonLoops<MatrixField,GaugeField>::Staple(dSdU_mu,U,mu);
-	  std::cout << "Staple norm ["<<mu<<"] = "<< norm2(dSdU_mu) <<"\n";
 	  dSdU_mu = Ta(Umu*adj(dSdU_mu))*factor;
-	  std::cout << "Force norm ["<<mu<<"] = "<< norm2(dSdU_mu) << "  : Umu "<< norm2(Umu)<<"\n";
 	  pokeLorentz(dSdU, dSdU_mu, mu);
 	}
       };

lib/qcd/hmc/HMC.cc

@@ -7,9 +7,9 @@ namespace Grid{
 	// FIXME fill this constructor  now just default values
 	  
 	////////////////////////////// Default values
-	Nsweeps             = 10;
+	Nsweeps             = 100;
-	TotalSweeps         = 10;
+	TotalSweeps         = 20;
-	ThermalizationSteps = 1;
+	ThermalizationSteps = 20;
 	StartingConfig      = 0;
 	SaveInterval        = 1;
 	Filename_prefix     = "Conf_";

lib/qcd/hmc/HMC.h

@@ -10,6 +10,7 @@
 
 #include <string>
 #include <memory>
+#include <unistd.h>
 
 namespace Grid{
   namespace QCD{
@@ -55,12 +56,12 @@ namespace Grid{
 
 	MD.init(U,pRNG);     // set U and initialize P and phi's 
 	RealD H0 = MD.S(U);     // current state            
-	std::cout<<"Total H_before = "<< H0 << "\n";
+	std::cout<<"Total H before = "<< H0 << "\n";
       
 	MD.integrate(U,0);
       
 	RealD H1 = MD.S(U);     // updated state            
-	std::cout<<"Total H_after = "<< H1 << "\n";
+	std::cout<<"Total H after = "<< H1 << "\n";
       
 	return (H1-H0);
       }
@@ -84,19 +85,14 @@ namespace Grid{
 
       void evolve(LatticeLorentzColourMatrix& Uin){
 	Real DeltaH;
-	Real timer;
       
 	// Thermalizations
 	for(int iter=1; iter <= Params.ThermalizationSteps; ++iter){
 	  std::cout << "-- # Thermalization step = "<< iter <<  "\n";
 	
 	  DeltaH = evolve_step(Uin);
-	
-	  std::cout<< "[Timing] Trajectory time (s) : "<< timer/1000.0 << "\n";
 	  std::cout<< " dH = "<< DeltaH << "\n";
 	
-	  // Update matrix
-	  //Uin = MD->get_U();  //accept every time
 	}
 
 	// Actual updates (evolve a copy Ucopy then copy back eventually)
@@ -109,13 +105,9 @@ namespace Grid{
 	  DeltaH = evolve_step(Ucopy);
 		
 	  if(metropolis_test(DeltaH)) Uin = Ucopy;
-	  // need sync?	
+	  //need sync?
 	}
-
-      
       }
-
-
     };
     
   }// QCD

lib/qcd/hmc/integrators/Integrator_base.cc

@@ -15,10 +15,12 @@ namespace Grid{
 
     void MDutils::generate_momenta_su3(LatticeLorentzColourMatrix& P,GridParallelRNG& pRNG){
       LatticeColourMatrix Pmu(P._grid);
+      Pmu = zero;
       for(int mu=0;mu<Nd;mu++){
 	SU3::GaussianLieAlgebraMatrix(pRNG, Pmu);
 	pokeLorentz(P, Pmu, mu);
       }
+      
     }
 
 

lib/qcd/hmc/integrators/Integrator_base.h

@@ -71,14 +71,13 @@ namespace Grid{
 
 
       void update_U(LatticeLorentzColourMatrix&U, double ep){
-	//rewrite exponential to deal with the lorentz index?
+	//rewrite exponential to deal automatically  with the lorentz index?
 	LatticeColourMatrix Umu(U._grid);
 	LatticeColourMatrix Pmu(U._grid);
 	for (int mu = 0; mu < Nd; mu++){
 	  Umu=peekLorentz(U, mu);
 	  Pmu=peekLorentz(*P, mu);
-	  std::cout << "U norm ["<<mu<<"] = "<< norm2(Umu) << "  : Pmu "<< norm2(Pmu)<<"\n";
+	  Umu = expMat(Pmu, ep)*Umu;
-	  Umu = expMat(Pmu, Complex(ep, 0.0))*Umu;
 	  pokeLorentz(U, Umu, mu);
 	}
 
@@ -105,14 +104,10 @@ namespace Grid{
 		GridParallelRNG& pRNG){
 	std::cout<< "Integrator init\n";
 	if (!P){
-	  std::unique_ptr<LatticeLorentzColourMatrix> Pnew(new LatticeLorentzColourMatrix(U._grid));
+	  std::unique_ptr<LatticeLorentzColourMatrix> Pnew(new LatticeLorentzColourMatrix(U._grid));	  
 	  P = std::move(Pnew);
-
 	}
 	MDutils::generate_momenta(*P,pRNG);
-	
-       
-	
 	for(int level=0; level< as.size(); ++level){
 	  for(int actionID=0; actionID<as.at(level).size(); ++actionID){
 	    as[level].at(actionID)->init(U, pRNG);
@@ -128,7 +123,7 @@ namespace Grid{
 	// Momenta
 	for (int mu=0; mu <Nd; mu++){
 	  LatticeColourMatrix Pmu = peekLorentz(*P, mu);
-	  Hloc -= trace(Pmu*adj(Pmu));
+	  Hloc -= trace(Pmu*Pmu);
 	}
 	Complex Hsum = sum(Hloc);
 	

lib/qcd/utils/SUn.h

@@ -545,21 +545,21 @@ Note that in step D setting B ~ X - A and using B in place of A in step E will g
 
   static void GaussianLieAlgebraMatrix(GridParallelRNG     &pRNG,LatticeMatrix &out,double scale=1.0){
     GridBase *grid = out._grid;
-    LatticeComplex ca (grid);
+    LatticeReal ca (grid);
-    LatticeMatrix  lie(grid);
     LatticeMatrix  la (grid);
     Complex ci(0.0,scale);
     Matrix ta;
 
     out=zero;
     for(int a=0;a<generators();a++){
-
       gaussian(pRNG,ca); 
       generator(a,ta);
-
+      
-      la=ci*ca*ta;
+      la=toComplex(ca)*ci*ta;
-      out = lie+la; // e^{i la ta}
+   
+      out += la; 
     }
+
   }
 
 

lib/tensors/Tensor_Ta.h

@@ -69,9 +69,7 @@ namespace Grid {
 	for(int c2=0;c2<N;c2++)
 	  inner += innerProduct(ret._internal[c1][c2],ret._internal[c1][c2]);
 
-	//nrm = 1.0/sqrt(Reduce(toReal(inner)));
 	nrm = rsqrt(inner);
-
 	for(int c2=0;c2<N;c2++)
 	  ret._internal[c1][c2]*= nrm;
       

tests/Make.inc

@@ -1,5 +1,5 @@
 
-bin_PROGRAMS = Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_contfrac_cg Test_contfrac_even_odd Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_fpgcr Test_dwf_hdcr Test_gamma Test_GaugeAction Test_hmc Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_quenched_update Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd
+bin_PROGRAMS = Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_contfrac_cg Test_contfrac_even_odd Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_dwf_fpgcr Test_dwf_hdcr Test_gamma Test_GaugeAction Test_hmc_WilsonGauge Test_lie_generators Test_main Test_multishift_sqrt Test_nersc_io Test_quenched_update Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd
 
 
 Test_cayley_cg_SOURCES=Test_cayley_cg.cc
@@ -78,8 +78,8 @@ Test_GaugeAction_SOURCES=Test_GaugeAction.cc
 Test_GaugeAction_LDADD=-lGrid
 
 
-Test_hmc_SOURCES=Test_hmc.cc
+Test_hmc_WilsonGauge_SOURCES=Test_hmc_WilsonGauge.cc
-Test_hmc_LDADD=-lGrid
+Test_hmc_WilsonGauge_LDADD=-lGrid
 
 
 Test_lie_generators_SOURCES=Test_lie_generators.cc

tests/Test_hmc_WilsonGauge.cc

@@ -26,7 +26,7 @@ int main (int argc, char ** argv)
   pRNG.SeedRandomDevice();
   LatticeLorentzColourMatrix     U(&Fine);
 
-  SU3::ColdConfiguration(pRNG, U);
+  SU3::HotConfiguration(pRNG, U);
  
 
   // simplify template?
@@ -39,7 +39,7 @@ int main (int argc, char ** argv)
   FullSet.push_back(Level1);
 
   // Create integrator
-  IntegratorParameters MDpar(12,50,1.0);
+  IntegratorParameters MDpar(12,30,1.0);
   std::vector<int> rel ={1};
   Integrator<LeapFrog> MDleapfrog(MDpar, FullSet,rel);