From 96f1d1b828f613a7d2278a441cc86ee59cfd2526 Mon Sep 17 00:00:00 2001
From: paboyle <paboyle@ph.ed.ac.uk>
Date: Mon, 10 Oct 2016 23:46:45 +0100
Subject: [PATCH] Debugged Domain wall and Overlap feynman rules (infinite Ls,
 finite mass).

---
 lib/qcd/action/fermion/WilsonFermion5D.cc | 184 ++++++++++------------
 1 file changed, 81 insertions(+), 103 deletions(-)

diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc
index 28f02b67..ea3e8535 100644
--- a/lib/qcd/action/fermion/WilsonFermion5D.cc
+++ b/lib/qcd/action/fermion/WilsonFermion5D.cc
@@ -393,92 +393,95 @@ void WilsonFermion5D<Impl>::DW(const FermionField &in, FermionField &out,int dag
 }
 
 template<class Impl>
-void WilsonFermion5D<Impl>::MomentumSpacePropagatorHt(FermionField &out,const FermionField &in) 
+void WilsonFermion5D<Impl>::MomentumSpacePropagatorHt(FermionField &out,const FermionField &in, RealD mass) 
 {
-    // what type LatticeComplex 
-    GridBase *_grid = _FourDimGrid;
-    conformable(_grid,out._grid);
+  // what type LatticeComplex 
+  GridBase *_grid = _FourDimGrid;
+  conformable(_grid,out._grid);
+  
+  typedef typename FermionField::vector_type vector_type;
+  typedef typename FermionField::scalar_type ScalComplex;
+  typedef iSinglet<ScalComplex> Tcomplex;
+  typedef Lattice<iSinglet<vector_type> > LatComplex;
+  
+  Gamma::GammaMatrix Gmu [] = {
+    Gamma::GammaX,
+    Gamma::GammaY,
+    Gamma::GammaZ,
+    Gamma::GammaT
+  };
 
-    typedef typename FermionField::vector_type vector_type;
-    typedef typename FermionField::scalar_type ScalComplex;
-    typedef iSinglet<ScalComplex> Tcomplex;
-    typedef Lattice<iSinglet<vector_type> > LatComplex;
+  std::vector<int> latt_size   = _grid->_fdimensions;
 
-    Gamma::GammaMatrix Gmu [] = {
-      Gamma::GammaX,
-      Gamma::GammaY,
-      Gamma::GammaZ,
-      Gamma::GammaT
-    };
+  
+  FermionField   num  (_grid); num  = zero;
 
-    std::vector<int> latt_size   = _grid->_fdimensions;
+  LatComplex    sk(_grid);  sk = zero;
+  LatComplex    sk2(_grid); sk2= zero;
+  LatComplex    W(_grid); W= zero;
+  LatComplex    a(_grid); a= zero;
+  LatComplex    one  (_grid); one = ScalComplex(1.0,0.0);
+  LatComplex denom(_grid); denom= zero;
+  LatComplex cosha(_grid); 
+  LatComplex kmu(_grid); 
+  LatComplex Wea(_grid); 
+  LatComplex Wema(_grid); 
 
-
-    LatComplex    sk(_grid);  sk = zero;
-    LatComplex    sk2(_grid); sk2= zero;
-
-    LatComplex    W(_grid); W= zero;
-    LatComplex    a(_grid); a= zero;
-    LatComplex    cosha(_grid); a= zero;
-
-    LatComplex     one  (_grid);     one = ScalComplex(1.0,0.0);
-
-    FermionField   num  (_grid); num  = zero;
-    LatComplex denom(_grid); denom= zero;
-    LatComplex kmu(_grid); 
-    ScalComplex ci(0.0,1.0);
-
-    for(int mu=0;mu<Nd;mu++) {
-
-      LatticeCoordinate(kmu,mu);
-
-      RealD TwoPiL =  M_PI * 2.0/ latt_size[mu];
-
-      kmu = TwoPiL * kmu;
-
-      sk2 = sk2 + 2.0*sin(kmu*0.5)*sin(kmu*0.5);
-      sk  = sk  + sin(kmu)*sin(kmu); 
-
-      num = num - sin(kmu)*ci*(Gamma(Gmu[mu])*in);
-
-    }
-
-    W = one - M5 - sk2;
-    cosha =  (one + W*W + sk) / (W*2.0);
+  ScalComplex ci(0.0,1.0);
+  
+  for(int mu=0;mu<Nd;mu++) {
     
-    for(int idx=0;idx<_grid->lSites();idx++){
-      Tcomplex cc;
-      RealD sgn;
-      std::vector<int> lcoor(Nd);
-      _grid->LocalIndexToLocalCoor(idx,lcoor);
-      peekLocalSite(cc,cosha,lcoor);
-      sgn= ::fabs(real(cc));
-      cc = ScalComplex(::acosh(sgn),0.0);
-      pokeLocalSite(cc,a,lcoor);
-    }
+    LatticeCoordinate(kmu,mu);
+    
+    RealD TwoPiL =  M_PI * 2.0/ latt_size[mu];
+    
+    kmu = TwoPiL * kmu;
+    
+    sk2 = sk2 + 2.0*sin(kmu*0.5)*sin(kmu*0.5);
+    sk  = sk  +     sin(kmu)    *sin(kmu); 
+    
+    num = num - sin(kmu)*ci*(Gamma(Gmu[mu])*in);
+    
+  }
+  
+  W = one - M5 + sk2;
 
-  std::cout << "Ht mom space num" << norm2(num)<<std::endl;// num ok
-  std::cout << "Ht mom space W  " << norm2(W)<<std::endl;// w ok
-  std::cout << "Ht mom space a  " << norm2(a)<<std::endl;// a ok
-  denom= ( exp(a) ); 
-  std::cout << "Ht mom space den" << norm2(denom)<<std::endl;
-  denom= ( exp(a) * W  ); 
-  std::cout << "Ht mom space den W" << norm2(denom)<<std::endl;
-  denom= ( exp(a) * W - one ); 
-  std::cout << "Ht mom space den W- 1" << norm2(denom)<<std::endl;
-  denom=where(sk==ScalComplex(0,0),one,denom);
-  std::cout << "Ht mom space den W- 1" << denom<<std::endl;
-  std::cout << "Ht mom space one " << norm2(one)<<std::endl;
-  denom= one/denom;
-  std::cout << "Ht mom space div " << norm2(denom)<<std::endl;
-  out = num*denom;
-  std::cout << "Ht mom space out" << norm2(out)<<std::endl;
+  ////////////////////////////////////////////
+  // Cosh alpha -> alpha
+  ////////////////////////////////////////////
+  cosha =  (one + W*W + sk) / (W*2.0);
+
+  // FIXME Need a Lattice acosh
+  for(int idx=0;idx<_grid->lSites();idx++){
+    std::vector<int> lcoor(Nd);
+    Tcomplex cc;
+    RealD sgn;
+    _grid->LocalIndexToLocalCoor(idx,lcoor);
+    peekLocalSite(cc,cosha,lcoor);
+    assert((double)real(cc)>=1.0);
+    assert(fabs((double)imag(cc))<=1.0e-15);
+    cc = ScalComplex(::acosh(real(cc)),0.0);
+    pokeLocalSite(cc,a,lcoor);
+  }
+  
+  Wea = ( exp( a) * W  ); 
+  Wema= ( exp(-a) * W  ); 
+  
+  num   = num + ( one - Wema ) * mass * in;
+  denom= ( Wea - one ) + mass*mass * (one - Wema); 
+  out = num/denom;
 }
 
 template<class Impl>
-void WilsonFermion5D<Impl>::MomentumSpacePropagatorHw(FermionField &out,const FermionField &in) 
+void WilsonFermion5D<Impl>::MomentumSpacePropagatorHw(FermionField &out,const FermionField &in,RealD mass) 
 {
-    // what type LatticeComplex 
+    Gamma::GammaMatrix Gmu [] = {
+      Gamma::GammaX,
+      Gamma::GammaY,
+      Gamma::GammaZ,
+      Gamma::GammaT
+    };
+
     GridBase *_grid = _FourDimGrid;
     conformable(_grid,out._grid);
 
@@ -487,16 +490,9 @@ void WilsonFermion5D<Impl>::MomentumSpacePropagatorHw(FermionField &out,const Fe
 
     typedef Lattice<iSinglet<vector_type> > LatComplex;
 
-    Gamma::GammaMatrix Gmu [] = {
-      Gamma::GammaX,
-      Gamma::GammaY,
-      Gamma::GammaZ,
-      Gamma::GammaT
-    };
 
     std::vector<int> latt_size   = _grid->_fdimensions;
 
-
     LatComplex    sk(_grid);  sk = zero;
     LatComplex    sk2(_grid); sk2= zero;
 
@@ -524,35 +520,17 @@ void WilsonFermion5D<Impl>::MomentumSpacePropagatorHw(FermionField &out,const Fe
       num = num - sin(kmu)*ci*(Gamma(Gmu[mu])*in);
 
     }
+    num = num + mass * in ;
 
     b_k = sk2 - M5;
      
     w_k = sqrt(sk + b_k*b_k);
 
-    std::cout << "Hw M5 " << M5<<std::endl;
-    std::cout << "Hw mom space NUM" << norm2(num)<<std::endl;
-    std::cout << "Hw mom space BK" << norm2(b_k)<<std::endl;
-    std::cout << "Hw mom space BK" << b_k<<std::endl;
-    std::cout << "Hw mom space WK" << norm2(w_k)<<std::endl;
-    std::cout << "Hw mom space WK" << w_k<<std::endl;
-
-    denom= ( w_k + b_k ) * (2.0 * M5);
-    denom=where(sk==ScalComplex(0,0),one,denom);
-
-    std::cout << "Hw mom space DEMON" << norm2(denom)<<std::endl;
-
-    std::cout << "Hw mom space DEMON" << denom<<std::endl;
-    
-    std::cout << "Hw mom space one" << norm2(one)<<std::endl;
+    denom= ( w_k + b_k + mass*mass) ;
 
     denom= one/denom;
-
-    std::cout << "Hw mom space div " << norm2(denom)<<std::endl;
-
-    std::cout << "denom " << denom<<std::endl;
-
     out = num*denom;
-    std::cout << "Hw mom space OUT" << norm2(out)<<std::endl;
+
 }