Debugging force term

lib/qcd/action/fermion/FermionOperatorImpl.h

@@ -254,8 +254,22 @@ namespace QCD {
       GaugeLinkField link(mat._grid);
       link = TraceIndex<SpinIndex>(outerProduct(Btilde,A)); 
       PokeIndex<LorentzIndex>(mat,link,mu);
-    }   
+    }  
+    
+    inline void outerProductImpl(PropagatorField &mat, const FermionField &B, const FermionField &A){
+      mat = outerProduct(B,A); 
+    }  
+
+    inline void TraceSpinImpl(GaugeLinkField &mat, PropagatorField&P) {
+      mat = TraceIndex<SpinIndex>(P); 
+    }
       
+    inline void extractLinkField(std::vector<GaugeLinkField> &mat, DoubledGaugeField &Uds){
+      for (int mu = 0; mu < Nd; mu++)
+      mat[mu] = PeekIndex<LorentzIndex>(Uds, mu);
+    }
+
+
     inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField &Atilde,int mu){
       
       int Ls=Btilde._grid->_fdimensions[0];
@@ -373,6 +387,19 @@ class DomainWallVec5dImpl :  public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepres
     assert(0);
   }
 
+  inline void outerProductImpl(PropagatorField &mat, const FermionField &Btilde, const FermionField &A){
+    assert(0);
+  } 
+
+  inline void TraceSpinImpl(GaugeLinkField &mat, PropagatorField&P) {
+    assert(0);
+  }
+
+  inline void extractLinkField(std::vector<GaugeLinkField> &mat, DoubledGaugeField &Uds){
+    assert(0);
+  }
+
+
   inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField &Atilde, int mu) {
 
     assert(0);
@@ -611,6 +638,25 @@ class GparityWilsonImpl : public ConjugateGaugeImpl<GaugeImplTypes<S, Nrepresent
    return;
  }
       
+ inline void outerProductImpl(PropagatorField &mat, const FermionField &Btilde, const FermionField &A){
+   //mat = outerProduct(Btilde, A);
+   assert(0);
+  }
+
+  inline void TraceSpinImpl(GaugeLinkField &mat, PropagatorField&P) {
+    assert(0);
+    /*
+    auto tmp = TraceIndex<SpinIndex>(P);
+    parallel_for(auto ss = tmp.begin(); ss < tmp.end(); ss++) {
+      mat[ss]() = tmp[ss](0, 0) + conjugate(tmp[ss](1, 1));
+    }
+    */
+  }
+
+  inline void extractLinkField(std::vector<GaugeLinkField> &mat, DoubledGaugeField &Uds){
+    assert(0);
+  }
+  
  inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField &Atilde, int mu) {
 
    int Ls = Btilde._grid->_fdimensions[0];
@@ -751,8 +797,8 @@ class StaggeredImpl : public PeriodicGaugeImpl<GaugeImplTypes<S, Representation:
       GaugeLinkField link(mat._grid);
       link = TraceIndex<SpinIndex>(outerProduct(Btilde,A)); 
       PokeIndex<LorentzIndex>(mat,link,mu);
-    }   
+    } 
-      
+          
     inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField &Atilde,int mu){
       assert (0); 
       // Must never hit

lib/qcd/action/fermion/WilsonCloverFermion.cc

@@ -43,11 +43,15 @@ RealD WilsonCloverFermion<Impl>::M(const FermionField &in, FermionField &out)
 
   // Wilson term
   out.checkerboard = in.checkerboard;
-  this->Dhop(in, out, DaggerNo);
+  //this->Dhop(in, out, DaggerNo);
 
   // Clover term
   Mooee(in, temp);
 
+  //hack
+  out = zero;
+
+
   out += temp;
   return norm2(out);
 }
@@ -59,11 +63,14 @@ RealD WilsonCloverFermion<Impl>::Mdag(const FermionField &in, FermionField &out)
 
   // Wilson term
   out.checkerboard = in.checkerboard;
-  this->Dhop(in, out, DaggerYes);
+  //this->Dhop(in, out, DaggerYes);
 
   // Clover term
   MooeeDag(in, temp);
 
+   //hack
+  out = zero;
+
   out += temp;
   return norm2(out);
 }
@@ -84,7 +91,7 @@ void WilsonCloverFermion<Impl>::ImportGauge(const GaugeField &_Umu)
   WilsonLoops<Impl>::FieldStrength(Ez, _Umu, Tdir, Zdir);
 
   // Compute the Clover Operator acting on Colour and Spin
-  CloverTerm = fillCloverYZ(Bx);
+  CloverTerm  = fillCloverYZ(Bx);
   CloverTerm += fillCloverXZ(By);
   CloverTerm += fillCloverXY(Bz);
   CloverTerm += fillCloverXT(Ex);
@@ -223,23 +230,6 @@ void WilsonCloverFermion<Impl>::MooeeInternal(const FermionField &in, FermionFie
 
 } // MooeeInternal
 
-// Derivative parts
-template <class Impl>
-void WilsonCloverFermion<Impl>::MDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag)
-{
-
-  GaugeField tmp(mat._grid);
-
-  conformable(U._grid, V._grid);
-  conformable(U._grid, mat._grid);
-
-  mat.checkerboard = U.checkerboard;
-  tmp.checkerboard = U.checkerboard;
-
-  this->DhopDeriv(mat, U, V, dag);
-  MooDeriv(tmp, U, V, dag);
-  mat += tmp;
-}
 
 // Derivative parts
 template <class Impl>

lib/qcd/action/fermion/WilsonCloverFermion.h

@@ -67,15 +67,18 @@ public:
                                                             CloverTermOdd(&Hgrid),
                                                             CloverTermInvEven(&Hgrid),
                                                             CloverTermInvOdd(&Hgrid),
-                                                            CloverTermDagEven(&Hgrid),    
+                                                            CloverTermDagEven(&Hgrid),
-                                                            CloverTermDagOdd(&Hgrid),     
+                                                            CloverTermDagOdd(&Hgrid),
-                                                            CloverTermInvDagEven(&Hgrid), 
+                                                            CloverTermInvDagEven(&Hgrid),
-                                                            CloverTermInvDagOdd(&Hgrid)   
+                                                            CloverTermInvDagOdd(&Hgrid)
   {
     csw = _csw;
     assert(Nd == 4); // require 4 dimensions
 
-    if (csw == 0) std::cout << GridLogWarning << "Initializing WilsonCloverFermion with csw = 0" << std::endl; 
+    if (csw == 0)
+      std::cout << GridLogWarning << "Initializing WilsonCloverFermion with csw = 0" << std::endl;
+
+    ImportGauge(_Umu);
   }
 
   virtual RealD M(const FermionField &in, FermionField &out);
@@ -87,16 +90,127 @@ public:
   virtual void MooeeInvDag(const FermionField &in, FermionField &out);
   virtual void MooeeInternal(const FermionField &in, FermionField &out, int dag, int inv);
 
-  virtual void MDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag);
+  //virtual void MDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag);
   virtual void MooDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag);
   virtual void MeeDeriv(GaugeField &mat, const FermionField &U, const FermionField &V, int dag);
 
   void ImportGauge(const GaugeField &_Umu);
 
+  // Derivative parts unpreconditioned pseudofermions
+  void MDeriv(GaugeField &force, const FermionField &X, const FermionField &Y, int dag)
+  {
+    conformable(X._grid, Y._grid);
+    conformable(X._grid, force._grid);
+    GaugeLinkField force_mu(force._grid), lambda(force._grid);
+    GaugeField clover_force(force._grid);
+    PropagatorField Lambda(force._grid);
+
+    // Here we are hitting some performance issues:
+    // need to extract the components of the DoubledGaugeField
+    // for each call
+    // Possible solution
+    // Create a vector object to store them? (cons: wasting space)
+    std::vector<GaugeLinkField> U(Nd, this->Umu._grid);
+    
+    Impl::extractLinkField(U, this->Umu);
+    
+    force = zero;
+    // Derivative of the Wilson hopping term
+    //this->DhopDeriv(force, X, Y, dag);
+
+    ///////////////////////////////////////////////////////////
+    // Clover term derivative
+    ///////////////////////////////////////////////////////////
+    Impl::outerProductImpl(Lambda, X, Y); 
+
+    Gamma::Algebra sigma[] = {
+        Gamma::Algebra::SigmaXY,
+        Gamma::Algebra::SigmaXZ,
+        Gamma::Algebra::SigmaXT,
+        Gamma::Algebra::MinusSigmaXY,
+        Gamma::Algebra::SigmaYZ,
+        Gamma::Algebra::SigmaYT,
+        Gamma::Algebra::MinusSigmaXZ,
+        Gamma::Algebra::MinusSigmaYZ,
+        Gamma::Algebra::SigmaZT,
+        Gamma::Algebra::MinusSigmaXT,
+        Gamma::Algebra::MinusSigmaYT,
+        Gamma::Algebra::MinusSigmaZT};
+
+    /*
+      sigma_{\mu \nu}=
+      | 0         sigma[0]  sigma[1]  sigma[2] |
+      | sigma[3]    0       sigma[4]  sigma[5] |
+      | sigma[6]  sigma[7]     0      sigma[8] |
+      | sigma[9]  sigma[10] sigma[11]   0      |
+    */
+
+    int count = 0;
+    clover_force  = zero;
+    for (int mu = 0; mu < 4; mu++)
+    {
+      force_mu = zero;
+      for (int nu = 0; nu < 4; nu++)
+      {
+        if (mu == nu) continue;
+        PropagatorField Slambda = Gamma(sigma[count]) * Lambda;
+        Impl::TraceSpinImpl(lambda, Slambda); //traceSpin
+        force_mu += Cmunu(U, lambda, mu, nu);
+        count++;
+      }
+
+      pokeLorentz(clover_force, U[mu] * force_mu, mu);
+    }
+    clover_force *= csw / 8.;
+    force += clover_force;
+  }
+
+  // Computing C_{\mu \nu}(x) as in Eq.(B.39) in Zbigniew Sroczynski's PhD thesis
+  GaugeLinkField Cmunu(std::vector<GaugeLinkField> &U, GaugeLinkField &lambda, int mu, int nu)
+  {
+    conformable(lambda._grid, U[0]._grid);
+    GaugeLinkField out(lambda._grid), tmp(lambda._grid);
+
+    // insertion in upper staple
+    // please check redundancy of shift operations
+
+    // C1+
+    tmp = lambda * U[nu];
+    out = Impl::ShiftStaple(Impl::CovShiftForward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu);
+
+    // C2+
+    tmp = U[mu] * Impl::CovShiftIdentityForward(adj(lambda), mu);
+    out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(tmp, mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu);
+
+    // C3+
+    tmp = U[nu] * Impl::CovShiftIdentityForward(adj(lambda), nu);
+    out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(tmp, nu))), mu);
+
+    // C4+
+    out += Impl::ShiftStaple(Impl::CovShiftForward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, Impl::CovShiftIdentityBackward(U[nu], nu))), mu) * lambda;
+
+    // insertion in lower staple
+    // C1-
+    out -= Impl::ShiftStaple(lambda, mu) * Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu);
+
+    // C2-
+    tmp = adj(lambda) * U[nu];
+    out -= Impl::ShiftStaple(Impl::CovShiftBackward(tmp, nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu);
+
+    // C3-
+    tmp = lambda * U[nu];
+    out -= Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, tmp)), mu);
+
+    // C4-
+    out -= Impl::ShiftStaple(Impl::CovShiftBackward(U[nu], nu, Impl::CovShiftBackward(U[mu], mu, U[nu])), mu) * lambda;
+
+    return out;
+  }
+
 private:
   // here fixing the 4 dimensions, make it more general?
 
-  RealD csw;                                                               // Clover coefficient
+  RealD csw;                                                 // Clover coefficient
   CloverFieldType CloverTerm, CloverTermInv;                 // Clover term
   CloverFieldType CloverTermEven, CloverTermOdd;             // Clover term EO
   CloverFieldType CloverTermInvEven, CloverTermInvOdd;       // Clover term Inv EO

tests/forces/Test_wilson_force.cc

@@ -50,7 +50,12 @@ int main (int argc, char ** argv)
   std::vector<int> seeds({1,2,3,4});
 
   GridParallelRNG          pRNG(&Grid);
-  pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
+  std::vector<int> vrand(4);
+  std::srand(std::time(0));
+  std::generate(vrand.begin(), vrand.end(), std::rand);
+  std::cout << GridLogMessage << vrand << std::endl;
+  pRNG.SeedFixedIntegers(vrand);
+  //pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
 
   LatticeFermion phi        (&Grid); gaussian(pRNG,phi);
   LatticeFermion Mphi       (&Grid); 

tests/forces/Test_wilsonclover_force.cc

@@ -0,0 +1,209 @@
+/*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./tests/Test_wilson_force.cc
+
+    Copyright (C) 2015
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+    *************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Grid.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+int main(int argc, char **argv)
+{
+  Grid_init(&argc, &argv);
+
+  std::vector<int> latt_size = GridDefaultLatt();
+  std::vector<int> simd_layout = GridDefaultSimd(Nd, vComplex::Nsimd());
+  std::vector<int> mpi_layout = GridDefaultMpi();
+
+  GridCartesian Grid(latt_size, simd_layout, mpi_layout);
+  GridRedBlackCartesian RBGrid(&Grid);
+
+  int threads = GridThread::GetThreads();
+  std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
+
+  std::vector<int> seeds({1, 2, 3, 4});
+
+  GridParallelRNG pRNG(&Grid);
+  std::vector<int> vrand(4);
+  std::srand(std::time(0));
+  std::generate(vrand.begin(), vrand.end(), std::rand);
+  std::cout << GridLogMessage << vrand << std::endl;
+  pRNG.SeedFixedIntegers(vrand);
+
+  LatticeFermion phi(&Grid);
+  gaussian(pRNG, phi);
+  LatticeFermion Mphi(&Grid);
+  LatticeFermion MphiPrime(&Grid);
+
+  LatticeGaugeField U(&Grid);
+
+  SU3::HotConfiguration(pRNG,U);
+  //SU3::ColdConfiguration(pRNG, U);
+
+  ////////////////////////////////////
+  // Unmodified matrix element
+  ////////////////////////////////////
+  RealD mass = -4.0; //kills the diagonal term
+  Real csw = 1.0;
+  WilsonCloverFermionR Dw(U, Grid, RBGrid, mass, csw);
+
+  Dw.M(phi, Mphi);
+  ComplexD S = innerProduct(Mphi, Mphi); // Action : pdag MdagM p
+
+  // get the deriv of phidag MdagM phi with respect to "U"
+  LatticeGaugeField UdSdU(&Grid);
+  LatticeGaugeField tmp(&Grid);
+
+  Dw.MDeriv(tmp, Mphi, phi, DaggerNo);  UdSdU  = tmp;
+  Dw.MDeriv(tmp, phi, Mphi, DaggerYes); UdSdU += tmp;
+  // Take the traceless antihermitian component
+  UdSdU = Ta(UdSdU);
+  
+  ////////////////////////////////////
+  // Modify the gauge field a little
+  ////////////////////////////////////
+  RealD dt = 0.0001;
+  RealD Hmom = 0.0;
+  RealD Hmomprime = 0.0;
+  RealD Hmompp = 0.0;
+  LatticeColourMatrix mommu(&Grid);
+  LatticeColourMatrix forcemu(&Grid);
+  LatticeGaugeField mom(&Grid);
+  LatticeGaugeField Uprime(&Grid);
+
+  for (int mu = 0; mu < Nd; mu++) {
+    // Traceless antihermitian momentum; gaussian in lie alg
+    SU3::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu);
+    Hmom -= real(sum(trace(mommu * mommu)));
+    PokeIndex<LorentzIndex>(mom, mommu, mu);
+  }
+
+  parallel_for(int ss=0;ss<mom._grid->oSites();ss++){
+    for (int mu = 0; mu < Nd; mu++) 
+      Uprime[ss]._internal[mu] = ProjectOnGroup(Exponentiate(mom[ss]._internal[mu], dt, 12) * U[ss]._internal[mu]);
+  } 
+
+  std::cout << GridLogMessage << "Initial mom hamiltonian is " << Hmom << std::endl;
+  
+  // New action
+  Dw.ImportGauge(Uprime);
+  Dw.M(phi, MphiPrime);
+
+  ComplexD Sprime = innerProduct(MphiPrime, MphiPrime);
+
+  //////////////////////////////////////////////
+  // Use derivative to estimate dS
+  //////////////////////////////////////////////
+
+  ///////////////////////////////////////////////////////
+  std::cout << GridLogMessage << "Antihermiticity tests - 1 " << std::endl;
+  for (int mu = 0; mu < Nd; mu++)
+  {
+    mommu = PeekIndex<LorentzIndex>(mom, mu);
+    std::cout << GridLogMessage << " Mommu  " << norm2(mommu) << std::endl;
+    mommu = mommu + adj(mommu);
+    std::cout << GridLogMessage << " Test: Mommu + Mommudag " << norm2(mommu) << std::endl;
+    mommu = PeekIndex<LorentzIndex>(UdSdU, mu);
+    std::cout << GridLogMessage << " dsdumu  " << norm2(mommu) << std::endl;
+    mommu = mommu + adj(mommu);
+    std::cout << GridLogMessage << " Test: dsdumu + dag  " << norm2(mommu) << std::endl;
+    std::cout << "" << std::endl;
+  }
+  ////////////////////////////////////////////////////////
+
+  LatticeComplex dS(&Grid);
+  dS = zero;
+  LatticeComplex dSmom(&Grid);
+  dSmom = zero;
+  LatticeComplex dSmom2(&Grid);
+  dSmom2 = zero;
+
+
+  // need for this??? 
+  // ultimately it is just a 2.0 factor in UdSdU
+  for (int mu = 0; mu < Nd; mu++)
+  {
+    mommu = PeekIndex<LorentzIndex>(UdSdU, mu); // P_mu = 
+    mommu = Ta(mommu) * 2.0;                    // Mom = (P_mu - P_mu^dag) - trace(P_mu - P_mu^dag)
+    PokeIndex<LorentzIndex>(UdSdU, mommu, mu);  // UdSdU_mu = Mom
+  }
+
+  std::cout << GridLogMessage<< "Antihermiticity tests - 2 " << std::endl;
+  for (int mu = 0; mu < Nd; mu++)
+  {
+    mommu = PeekIndex<LorentzIndex>(mom, mu);
+    std::cout << GridLogMessage << " Mommu  " << norm2(mommu) << std::endl;
+    mommu = mommu + adj(mommu);
+    std::cout << GridLogMessage << " Mommu + Mommudag " << norm2(mommu) << std::endl;
+    mommu = PeekIndex<LorentzIndex>(UdSdU, mu);
+    std::cout << GridLogMessage << " dsdumu  " << norm2(mommu) << std::endl;
+    mommu = mommu + adj(mommu);
+    std::cout << GridLogMessage << " dsdumu + dag  " << norm2(mommu) << std::endl;
+    std::cout << "" << std::endl;
+  }
+  /////////////////////////////////////////////////////
+
+
+  for (int mu = 0; mu < Nd; mu++)
+  {
+    forcemu = PeekIndex<LorentzIndex>(UdSdU, mu);
+    mommu = PeekIndex<LorentzIndex>(mom, mu);
+
+    // Update PF action density
+    dS = dS + trace(mommu * forcemu) * dt;
+
+    dSmom = dSmom - trace(mommu * forcemu) * dt;
+    dSmom2 = dSmom2 - trace(forcemu * forcemu) * (0.25 * dt * dt);
+
+    // Update mom action density
+    mommu = mommu + forcemu * (dt * 0.5);
+
+    Hmomprime -= real(sum(trace(mommu * mommu)));
+  }
+
+  ComplexD dSpred = sum(dS);
+  ComplexD dSm = sum(dSmom);
+  ComplexD dSm2 = sum(dSmom2);
+
+  std::cout << GridLogMessage << "Initial mom hamiltonian is " << Hmom << std::endl;
+  std::cout << GridLogMessage << "Final   mom hamiltonian is " << Hmomprime << std::endl;
+  std::cout << GridLogMessage << "Delta   mom hamiltonian is " << Hmomprime - Hmom << std::endl;
+
+  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 << "predict dS    " << dSpred << std::endl;
+  std::cout << GridLogMessage << "dSm " << dSm << std::endl;
+  std::cout << GridLogMessage << "dSm2" << dSm2 << std::endl;
+
+  std::cout << GridLogMessage << "Total dS    " << Hmomprime - Hmom + Sprime - S << std::endl;
+
+  assert(fabs(real(Sprime - S - dSpred)) < 1.0);
+
+  std::cout << GridLogMessage << "Done" << std::endl;
+  Grid_finalize();
+}