Added TM fermions for DSDR etc..

lib/qcd/action/fermion/WilsonFermion.cc

@@ -292,6 +292,7 @@ PARALLEL_FOR_LOOP
  
   FermOpTemplateInstantiate(WilsonFermion);
 
+
 }}
 
 

lib/qcd/action/fermion/WilsonFermion.h

@@ -42,10 +42,12 @@ namespace Grid {
       /////////////////////////////////////////////////////////
       void Meooe(const FermionField &in, FermionField &out) ;
       void MeooeDag(const FermionField &in, FermionField &out) ;
-      void Mooee(const FermionField &in, FermionField &out) ;
-      void MooeeDag(const FermionField &in, FermionField &out) ;
-      void MooeeInv(const FermionField &in, FermionField &out) ;
-      void MooeeInvDag(const FermionField &in, FermionField &out) ;
+
+      // allow override for twisted mass and clover
+      virtual void Mooee(const FermionField &in, FermionField &out) ;
+      virtual void MooeeDag(const FermionField &in, FermionField &out) ;
+      virtual void MooeeInv(const FermionField &in, FermionField &out) ;
+      virtual void MooeeInvDag(const FermionField &in, FermionField &out) ;
 
       ////////////////////////
       // Derivative interface

lib/qcd/action/fermion/WilsonTMFermion.cc

@@ -0,0 +1,71 @@
+#include <Grid.h>
+
+namespace Grid {
+namespace QCD {
+
+    /*
+     * BF sequence
+     *
+      void bfmbase<Float>::MooeeInv(Fermion_t psi, 
+			       Fermion_t chi, 
+			      int dag, int cb)
+
+    double m    = this->mass;
+    double tm   = this->twistedmass;
+    double mtil = 4.0+this->mass;
+
+    double sq = mtil*mtil + tm*tm;
+
+    double a = mtil/sq;
+    double b = -tm /sq;
+    if(dag) b=-b;
+    axpibg5x(chi,psi,a,b);
+
+      void bfmbase<Float>::Mooee(Fermion_t psi, 
+			   Fermion_t chi, 
+			   int dag,int cb)
+    double a = 4.0+this->mass;
+    double b = this->twistedmass;
+    if(dag) b=-b;
+    axpibg5x(chi,psi,a,b);
+    */
+
+  template<class Impl>
+  void WilsonTMFermion<Impl>::Mooee(const FermionField &in, FermionField &out) {
+    RealD a = 4.0+this->mass;
+    RealD b = this->mu;
+    out.checkerboard = in.checkerboard;
+    axpibg5x(out,in,a,b);
+  }
+  template<class Impl>
+  void WilsonTMFermion<Impl>::MooeeDag(const FermionField &in, FermionField &out) {
+    RealD a = 4.0+this->mass;
+    RealD b = -this->mu;
+    out.checkerboard = in.checkerboard;
+    axpibg5x(out,in,a,b);
+  }
+  template<class Impl>
+  void WilsonTMFermion<Impl>::MooeeInv(const FermionField &in, FermionField &out) {
+    RealD m    = this->mass;
+    RealD tm   = this->mu;
+    RealD mtil = 4.0+this->mass;
+    RealD sq   = mtil*mtil+tm*tm;
+    RealD a    = mtil/sq;
+    RealD b    = -tm /sq;
+    axpibg5x(out,in,a,b);
+  }
+  template<class Impl>
+  void WilsonTMFermion<Impl>::MooeeInvDag(const FermionField &in, FermionField &out) {
+    RealD m    = this->mass;
+    RealD tm   = this->mu;
+    RealD mtil = 4.0+this->mass;
+    RealD sq   = mtil*mtil+tm*tm;
+    RealD a    = mtil/sq;
+    RealD b    = tm /sq;
+    axpibg5x(out,in,a,b);
+  }
+
+  FermOpTemplateInstantiate(WilsonTMFermion);
+
+}
+}

lib/qcd/action/fermion/WilsonTMFermion.h

@@ -0,0 +1,47 @@
+#ifndef  GRID_QCD_WILSON_TM_FERMION_H
+#define  GRID_QCD_WILSON_TM_FERMION_H
+
+#include <Grid.h>
+
+namespace Grid {
+
+  namespace QCD {
+
+    template<class Impl>
+    class WilsonTMFermion : public WilsonFermion<Impl>
+    {
+    public:
+     INHERIT_IMPL_TYPES(Impl);
+    public:
+     RealD mu; // TwistedMass parameter
+
+      virtual void   Instantiatable(void) {};
+      // Constructors
+      WilsonTMFermion(GaugeField &_Umu,
+		    GridCartesian         &Fgrid,
+		    GridRedBlackCartesian &Hgrid, 
+		    RealD _mass,
+		    RealD _mu,
+		    const ImplParams &p= ImplParams()
+		      ) :
+	WilsonFermion<Impl>(_Umu,
+			    Fgrid,
+			    Hgrid,
+			    _mass,p)
+
+      {
+	mu = _mu;
+      }
+
+    };
+    // allow override for twisted mass and clover
+    virtual void Mooee(const FermionField &in, FermionField &out) ;
+    virtual void MooeeDag(const FermionField &in, FermionField &out) ;
+    virtual void MooeeInv(const FermionField &in, FermionField &out) ;
+    virtual void MooeeInvDag(const FermionField &in, FermionField &out) ;
+    
+
+  }
+}
+
+#endif

lib/qcd/utils/LinalgUtils.h

@@ -9,6 +9,25 @@ namespace QCD{
 //These routines support five-D chiral fermions and contain s-subslice indexing 
 //on the 5d (rb4d) checkerboarded lattices
 ////////////////////////////////////////////////////////////////////////
+
+template<class vobj> 
+void axibg5x(Lattice<vobj> &z,const Lattice<vobj> &x,RealD a,RealD b)
+{
+  z.checkerboard = x.checkerboard;
+  conformable(x,z);
+
+  GridBase *grid=x._grid;
+
+  Gamma G5(Gamma::Gamma5);
+PARALLEL_FOR_LOOP
+  for(int ss=0;ss<grid->oSites();ss++){
+    vobj tmp;
+    tmp = a*x._odata[ss];
+    tmp = tmp + G5*(b*timesI(x._odata[ss]));
+    vstream(z._odata[ss],tmp);
+  }
+}
+
 template<class vobj> 
 void axpby_ssp(Lattice<vobj> &z, RealD a,const Lattice<vobj> &x,RealD b,const Lattice<vobj> &y,int s,int sp)
 {

scripts/hmc.sh

@@ -24,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 [-0.1:0.6] > plot.gnu
+echo set yrange [-0.2:1.0] > 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_main.cc

@@ -35,7 +35,7 @@ int main (int argc, char ** argv)
 #ifdef AVX512
  for(int omp=128;omp<236;omp+=16){
 #else
- for(int omp=1;omp<8;omp*=20){
+ for(int omp=1;omp<2;omp*=20){
 #endif
 
 #ifdef OMP
@@ -43,6 +43,9 @@ int main (int argc, char ** argv)
 #endif 
 
   for(int lat=8;lat<=16;lat+=40){
+
+    std::cout << "Lat "<<lat<<std::endl;
+
     latt_size[0] = lat;
     latt_size[1] = lat;
     latt_size[2] = lat;
@@ -53,8 +56,19 @@ int main (int argc, char ** argv)
     GridRedBlackCartesian rbFine(latt_size,simd_layout,mpi_layout);
     GridParallelRNG       FineRNG(&Fine);
     GridSerialRNG       SerialRNG;
+
     FineRNG.SeedRandomDevice();
 
+    std::cout <<"SerialRNG" << SerialRNG._generators[0] <<std::endl;
+    std::stringstream output(std::stringstream::out|std::stringstream::binary);
+    output <<SerialRNG._generators[0]<<std::endl;
+    std::cout << output.str();
+
+    {
+      std::ofstream of("rngstate",std::ios::out|std::ios::binary);
+      of << SerialRNG._generators[0];
+    }
+
     LatticeColourMatrix Foo(&Fine);
     LatticeColourMatrix Bar(&Fine);
 
@@ -302,7 +316,6 @@ int main (int argc, char ** argv)
     }
 
     FooBar = Bar;
- 
     /*
     { 
       std::vector<int> coor(4);
@@ -315,6 +328,7 @@ int main (int argc, char ** argv)
     */
     lex_sites(Foo);
 
+
     Integer mm[4];
     mm[0]=1;
     mm[1]=Fine._rdimensions[0];
@@ -330,8 +344,12 @@ int main (int argc, char ** argv)
       
     }
 
-    Bar = zero;
-    Bar = where(lex<Integer(10),Foo,Bar);
+
+
+    //    Bar = zero;
+    //    Bar = where(lex<Integer(10),Foo,Bar);
+
+
     cout << "peeking sites..\n";
     {
       std::vector<int> coor(4);