diff = divA - 2.0 * m_q * ps_ps - 2.0*mpps_ps ;
diff = divA - 2.0 * m_q * ps_ps - 2.0*mpps_ps + 2.0*q_bar_q;
void check_dwf_ward_identity(const multi1d<LatticeColorMatrix>& u,			     const multi1d<LatticePropagator>& p5d,			     const LatticePropagator& q_q,			     const LatticePropagator& q_mp_q,			     const Real& m_q){  QDPIO::cout<<"check_dwf_ward_identity: Checking the chiral Ward Identity...";  QDPIO::cout<<endl ;  LatticeComplex divA ;  divA = 0.0 ;  for(int mu(0);mu<Nd;mu++){    LatticeComplex tt ;    dwf_conserved_axial_ps_corr(tt,u,p5d,mu);    divA += tt - shift(tt,BACKWARD,mu) ;   }  LatticeComplex mpps_ps, ps_ps, contact ;  mpps_ps = localNorm2(q_mp_q) ;  ps_ps = localNorm2(q_q) ;  //contact =   LatticeComplex diff ;  diff = divA - 2.0 * m_q * ps_ps - 2.0*mpps_ps ;  multi1d<Double> corr ;       SftMom trick(0,false,3) ;  corr = sumMulti(localNorm2(diff), trick.getSet());  QDPIO::cout<<"check_dwf_ward_identity: Ward Identity violation per timeslice: "<<endl;  for(int t(0);t<corr.size(); t++){    QDPIO::cout<<"        "<<t<<"                "<< corr[t]<<endl ;  }      QDPIO::cout<<"check_dwf_ward_identity: Ward Identity violation: ";  QDPIO::cout<<norm2(diff)<<endl ;  QDPIO::cout<<"check_dwf_ward_identity: |divA|^2    : "<< norm2(divA)<<endl;  QDPIO::cout<<"check_dwf_ward_identity: |ps_ps|^2   : "<< norm2(ps_ps)<<endl;  QDPIO::cout<<"check_dwf_ward_identity: |mpps_ps|^2 : "<< norm2(mpps_ps)<<endl;}