diff --git a/Grid/tensors/Tensor_outer.h b/Grid/tensors/Tensor_outer.h index f45babeb..a32a2a91 100644 --- a/Grid/tensors/Tensor_outer.h +++ b/Grid/tensors/Tensor_outer.h @@ -45,17 +45,6 @@ accelerator_inline RR outerProduct(const RR &l, const RR& r) return l*r; } -template = 0> -accelerator_inline CC outerProduct(const CC &l, const CC& r) -{ - return l*conj(r); -} -template = 0> -accelerator_inline RR outerProduct(const RR &l, const RR& r) -{ - return l*r; -} - template accelerator_inline auto outerProduct (const iVector& lhs,const iVector& rhs) -> iMatrix { diff --git a/scripts/hmc.sh b/scripts/hmc.sh index c2cd6593..4d7191ff 100755 --- a/scripts/hmc.sh +++ b/scripts/hmc.sh @@ -1,19 +1,27 @@ #!/bin/bash LOG=$1 -SWEEPS=`grep dH $LOG | wc -l` -SWEEPS=`expr $SWEEPS - 80` +SWEEPS=`grep dH.= $LOG | wc -l` +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} ' ` -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) } ' ` +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+$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) } ' ` -edH=`grep dH $LOG | tail -n $SWEEPS | awk '{ S=S+exp(-$10)} END { print S/NR} '` -echo ": $edH" +grep Plaq $LOG | tail -n $SWEEPS | awk '{ S=S+$12/20; if(NR%20==0){ print NR/20, " ", S; S=0;} } ' > plaq.binned + +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 ": $edH (${dedH})" echo ": $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 dH $LOG | awk '{ print $10 }' > dH.dat -echo set yrange [-0.2:1.0] > plot.gnu +grep Plaq $LOG | awk '{ print $12 }' | uniq > plaq.dat +grep dH.= $LOG | awk '{ print $16 }' > dH.dat +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 diff --git a/tests/forces/Test_gpdwf_force.cc b/tests/forces/Test_gpdwf_force.cc index d6744080..b1c0dd80 100644 --- a/tests/forces/Test_gpdwf_force.cc +++ b/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; - std::vector twists(Nd,0); twists[nu] = 1; + const int nu = 1; + std::vector 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 > coor(UGrid); + LatticeCoordinate(coor,nu); for(int mu=0;mu::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg + // Traceless antihermitian momentum; gaussian in lie alg + SU::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); + if(0){ + if(mu==nu){ + mommu=where(coor==Lnu-1,mommu,zz); + } else { + mommu=Zero(); + } + } PokeIndex(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 "< twists(Nd,0); twists[nu] = 1; + const int nu = 1; + const int Lnu=latt_size[nu]; + + std::vector 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 > coor(UGrid); + LatticeCoordinate(coor,nu); + zz=Zero(); for(int mu=0;mu::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); - + SU::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); + if(0){ + if(mu==nu){ + mommu=where(coor==Lnu-1,mommu,zz); + } else { + mommu=Zero(); + } + } PokeIndex(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(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 "< twists(Nd,0); twists[nu] = 1; ConjugateGimplD::setDirections(twists);