Merge branch 'feature/gparity_HMC' of https://github.com/paboyle/Grid into feature/gparity_HMC

Grid/tensors/Tensor_outer.h

@@ -45,17 +45,6 @@ accelerator_inline RR outerProduct(const RR &l, const RR& r)
   return l*r;
 }
 
-template<class CC,IfComplex<CC> = 0>
-accelerator_inline CC outerProduct(const CC &l, const CC& r)
-{
-  return l*conj(r);
-}
-template<class RR,IfReal<RR> = 0>
-accelerator_inline RR outerProduct(const RR &l, const RR& r)
-{
-  return l*r;
-}
-
 template<class l,class r,int N> accelerator_inline
 auto outerProduct (const iVector<l,N>& lhs,const iVector<r,N>& rhs) -> iMatrix<decltype(outerProduct(lhs._internal[0],rhs._internal[0])),N>
 {

scripts/hmc.sh

@@ -1,19 +1,27 @@
 #!/bin/bash
 
 LOG=$1
-SWEEPS=`grep dH $LOG | wc -l`
+SWEEPS=`grep dH.= $LOG | wc -l`
-SWEEPS=`expr $SWEEPS - 80`
+SWEEPS=`expr $SWEEPS - 100`
 echo
 echo $SWEEPS thermalised sweeps
 echo
-plaq=`grep Plaq $LOG | tail -n $SWEEPS | awk '{ S=S+$10} END { print S/NR} ' `
+plaq=`grep Plaq $LOG | tail -n $SWEEPS | awk '{ S=S+$12} END { print S/NR} ' `
-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) } ' `
+plaqe=`grep Plaq $LOG | tail -n $SWEEPS | awk '{ S=S+$12 ; SS=SS+$12*$12 } 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) } ' `
+grep  Plaq $LOG | tail -n $SWEEPS | awk '{ S=S+$12/20; if(NR%20==0){ print NR/20, " ", S; S=0;} } '  > plaq.binned
-edH=`grep dH $LOG | tail -n $SWEEPS | awk '{ S=S+exp(-$10)} END { print S/NR} '`
+
-echo "<e-dH>: $edH"
+plaq=`cat plaq.binned  | awk '{ S=S+$2} END { print S/NR} ' `
+plaqe=`cat plaq.binned | awk '{ S=S+$2 ; SS=SS+$2*$2 } END { print sqrt( (SS/NR - S*S/NR/NR)/NR) } ' `
+echo "Binned Plaquette: $plaq (${plaqe})"
+echo
+
+dHv=`grep dH.= $LOG | tail -n $SWEEPS | awk '{ S=S+$16 ; SS=SS+$16*$16 } END { print sqrt(SS/NR) } ' `
+edH=`grep dH.= $LOG | tail -n $SWEEPS | awk '{ S=S+exp(-$16)} END { print S/NR} '`
+dedH=`grep dH.= $LOG | tail -n $SWEEPS | awk '{ S=S+exp(-$16); SS=SS+exp(-$16)*exp(-$16)} END { print sqrt( (SS/NR - S*S/NR/NR)/NR) } '`
+echo "<e-dH>: $edH (${dedH})"
 echo "<rms dH>: $dHv"
 
 TRAJ=`grep Acc $LOG | wc -l`
@@ -22,12 +30,13 @@ PACC=`expr  100 \* ${ACC} / ${TRAJ} `
 echo
 echo "Acceptance $PACC %  $ACC / $TRAJ "
 
-grep Plaq $LOG | awk '{ print $10 }' | uniq > plaq.dat
+grep Plaq $LOG | awk '{ print $12 }' | uniq > plaq.dat
-grep dH $LOG | awk '{ print $10 }' > dH.dat
+grep dH.= $LOG | awk '{ print $16 }' > dH.dat
-echo set yrange [-0.2:1.0] > plot.gnu
+echo set yrange [0.58:0.60] > plot.gnu
 echo set terminal 'pdf' >> plot.gnu
+echo "f(x) =0.588" >> plot.gnu
 echo "set output 'plaq.${LOG}.pdf'" >> plot.gnu
-echo "plot 'plaq.dat' w l, 'dH.dat' w l " >> plot.gnu
+echo "plot 'plaq.dat' w l, f(x) " >> plot.gnu
 echo
 gnuplot plot.gnu >& gnu.errs
 open plaq.${LOG}.pdf

tests/forces/Test_gpdwf_force.cc

@@ -71,8 +71,10 @@ int main (int argc, char ** argv)
   RealD mass=0.01; 
   RealD M5=1.8; 
 
-  const int nu = 3;
+  const int nu = 1;
-  std::vector<int> twists(Nd,0);  twists[nu] = 1;
+  std::vector<int> twists(Nd,0);
+  twists[nu] = 1;
+  twists[3] = 1;
   GparityDomainWallFermionR::ImplParams params;  params.twists = twists;
   GparityDomainWallFermionR Ddwf(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,params);
   Ddwf.M   (phi,Mphi);
@@ -91,16 +93,28 @@ int main (int argc, char ** argv)
   ////////////////////////////////////
   // Modify the gauge field a little 
   ////////////////////////////////////
-  RealD dt = 0.0001;
+  RealD dt = 0.01;
 
+  LatticeColourMatrix zz(UGrid); zz=Zero();
   LatticeColourMatrix mommu(UGrid); 
   LatticeColourMatrix forcemu(UGrid); 
   LatticeGaugeField mom(UGrid); 
   LatticeGaugeField Uprime(UGrid); 
 
+  const int Lnu=latt_size[nu];
+  Lattice<iScalar<vInteger> > coor(UGrid);
+  LatticeCoordinate(coor,nu);
   for(int mu=0;mu<Nd;mu++){
 
-    SU<Nc>::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
+    // Traceless antihermitian momentum; gaussian in lie alg
+    SU<Nc>::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu);
+    if(0){
+      if(mu==nu){
+	mommu=where(coor==Lnu-1,mommu,zz);
+      } else {
+	mommu=Zero();
+      }
+    }
 
     PokeIndex<LorentzIndex>(mom,mommu,mu);
 
@@ -125,6 +139,12 @@ int main (int argc, char ** argv)
 
   ComplexD Sprime    = innerProduct(MphiPrime   ,MphiPrime);
 
+
+  LatticeComplex lip(FGrid); lip=localInnerProduct(Mphi,Mphi);
+  LatticeComplex lipp(FGrid); lipp=localInnerProduct(MphiPrime,MphiPrime);
+  LatticeComplex dip(FGrid); dip = lipp - lip;
+  std::cout << " dip "<<dip<<std::endl;
+  
   //////////////////////////////////////////////
   // Use derivative to estimate dS
   //////////////////////////////////////////////

tests/forces/Test_gpwilson_force.cc

@@ -64,8 +64,12 @@ int main (int argc, char ** argv)
   ////////////////////////////////////
   RealD mass=0.01; 
 
-  const int nu = 3;
+  const int nu = 1;
-  std::vector<int> twists(Nd,0);  twists[nu] = 1;
+  const int Lnu=latt_size[nu];
+
+  std::vector<int> twists(Nd,0);
+  twists[nu] = 1;
+  twists[3]=1;
   GparityWilsonFermionR::ImplParams params;  params.twists = twists;
   GparityWilsonFermionR Wil(U,*UGrid,*UrbGrid,mass,params);
   Wil.M   (phi,Mphi);
@@ -87,17 +91,28 @@ int main (int argc, char ** argv)
   RealD dt = 0.01;
 
   LatticeColourMatrix mommu(UGrid); 
+  LatticeColourMatrix zz(UGrid);
   LatticeColourMatrix forcemu(UGrid); 
   LatticeGaugeField mom(UGrid); 
   LatticeGaugeField Uprime(UGrid); 
 
+  
+  Lattice<iScalar<vInteger> > coor(UGrid);
+  LatticeCoordinate(coor,nu);
+  zz=Zero();
   for(int mu=0;mu<Nd;mu++){
 
     // Traceless antihermitian momentum; gaussian in lie alg
-    SU<Nc>::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); 
+    SU<Nc>::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu);
-
+    if(0){
+      if(mu==nu){
+	mommu=where(coor==Lnu-1,mommu,zz);
+      } else {
+	mommu=Zero();
+      }
+    }
     PokeIndex<LorentzIndex>(mom,mommu,mu);
-
+    
     // fourth order exponential approx
     autoView( mom_v, mom, CpuRead);
     autoView( U_v , U, CpuRead);
@@ -130,6 +145,10 @@ int main (int argc, char ** argv)
     mommu=Ta(mommu)*2.0;
     PokeIndex<LorentzIndex>(UdSdU,mommu,mu);
   }
+  LatticeComplex lip(UGrid); lip=localInnerProduct(Mphi,Mphi);
+  LatticeComplex lipp(UGrid); lipp=localInnerProduct(MphiPrime,MphiPrime);
+  LatticeComplex dip(UGrid); dip = lipp - lip;
+  std::cout << " dip "<<dip<<std::endl;
 
   LatticeComplex dS(UGrid); dS = Zero();
   for(int mu=0;mu<Nd;mu++){
@@ -139,12 +158,14 @@ int main (int argc, char ** argv)
     // Update PF action density
     dS = dS+trace(mommu*forcemu)*dt;
   }
-
+  std::cout << "mommu"<<mommu<<std::endl;
+  std::cout << "dS" << dS<<std::endl;
+  
   ComplexD dSpred    = sum(dS);
 
   std::cout << GridLogMessage << " S      "<<S<<std::endl;
   std::cout << GridLogMessage << " Sprime "<<Sprime<<std::endl;
-  std::cout << GridLogMessage << "dS      "<<Sprime-S<<std::endl;
+  std::cout << GridLogMessage << "Delta S "<<Sprime-S<<std::endl;
   std::cout << GridLogMessage << "predict dS    "<< dSpred <<std::endl;
 
   assert( fabs(real(Sprime-S-dSpred)) < 2.0 ) ;

tests/hmc/Test_hmc_GparityIwasakiGauge.cc

@@ -58,7 +58,7 @@ int main(int argc, char **argv) {
   CheckpointerParameters CPparams;  
   CPparams.config_prefix = "ckpoint_EODWF_lat";
   CPparams.rng_prefix = "ckpoint_EODWF_rng";
-  CPparams.saveInterval = 5;
+  CPparams.saveInterval = 1;
   CPparams.format = "IEEE64BIG";
   
   TheHMC.Resources.LoadNerscCheckpointer(CPparams);
@@ -79,7 +79,7 @@ int main(int argc, char **argv) {
   // that have a complex construction
   // standard
   RealD beta = 2.6 ;
-  const int nu = 3;
+  const int nu = 1;
   std::vector<int> twists(Nd,0);
   twists[nu] = 1;
   ConjugateGimplD::setDirections(twists);