diff --git a/Grid/algorithms/iterative/ConjugateGradient.h b/Grid/algorithms/iterative/ConjugateGradient.h
index ff4ba8ac..8a6a4bae 100644
--- a/Grid/algorithms/iterative/ConjugateGradient.h
+++ b/Grid/algorithms/iterative/ConjugateGradient.h
@@ -150,13 +150,13 @@ class ConjugateGradient : public OperatorFunction<Field> {
 	std::cout << GridLogMessage << "\tTrue residual " << true_residual<<std::endl;
 	std::cout << GridLogMessage << "\tTarget " << Tolerance << std::endl;
 
-        std::cout << GridLogMessage << "Time breakdown "<<std::endl;
-	std::cout << GridLogMessage << "\tElapsed    " << SolverTimer.Elapsed() <<std::endl;
-	std::cout << GridLogMessage << "\tMatrix     " << MatrixTimer.Elapsed() <<std::endl;
-	std::cout << GridLogMessage << "\tLinalg     " << LinalgTimer.Elapsed() <<std::endl;
-	std::cout << GridLogMessage << "\tInner      " << InnerTimer.Elapsed() <<std::endl;
-	std::cout << GridLogMessage << "\tAxpyNorm   " << AxpyNormTimer.Elapsed() <<std::endl;
-	std::cout << GridLogMessage << "\tLinearComb " << LinearCombTimer.Elapsed() <<std::endl;
+        std::cout << GridLogPerformance << "Time breakdown "<<std::endl;
+	std::cout << GridLogPerformance << "\tElapsed    " << SolverTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tMatrix     " << MatrixTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tLinalg     " << LinalgTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tInner      " << InnerTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tAxpyNorm   " << AxpyNormTimer.Elapsed() <<std::endl;
+	std::cout << GridLogPerformance << "\tLinearComb " << LinearCombTimer.Elapsed() <<std::endl;
 
         if (ErrorOnNoConverge) assert(true_residual / Tolerance < 10000.0);
 
diff --git a/Grid/qcd/QCD.h b/Grid/qcd/QCD.h
index e1e38e82..7e9b2da0 100644
--- a/Grid/qcd/QCD.h
+++ b/Grid/qcd/QCD.h
@@ -90,17 +90,20 @@ namespace QCD {
     // That probably makes for GridRedBlack4dCartesian grid.
 
     // s,sp,c,spc,lc
-    template<typename vtype> using iSinglet                   = iScalar<iScalar<iScalar<vtype> > >;
-    template<typename vtype> using iSpinMatrix                = iScalar<iMatrix<iScalar<vtype>, Ns> >;
-    template<typename vtype> using iColourMatrix              = iScalar<iScalar<iMatrix<vtype, Nc> > > ;
-    template<typename vtype> using iSpinColourMatrix          = iScalar<iMatrix<iMatrix<vtype, Nc>, Ns> >;
-    template<typename vtype> using iLorentzColourMatrix       = iVector<iScalar<iMatrix<vtype, Nc> >, Nd > ;
-    template<typename vtype> using iDoubleStoredColourMatrix  = iVector<iScalar<iMatrix<vtype, Nc> >, Nds > ;
-    template<typename vtype> using iSpinVector                = iScalar<iVector<iScalar<vtype>, Ns> >;
-    template<typename vtype> using iColourVector              = iScalar<iScalar<iVector<vtype, Nc> > >;
-    template<typename vtype> using iSpinColourVector          = iScalar<iVector<iVector<vtype, Nc>, Ns> >;
-    template<typename vtype> using iHalfSpinVector            = iScalar<iVector<iScalar<vtype>, Nhs> >;
-    template<typename vtype> using iHalfSpinColourVector      = iScalar<iVector<iVector<vtype, Nc>, Nhs> >;
+
+    template<typename vtype> using iSinglet                     = iScalar<iScalar<iScalar<vtype> > >;
+    template<typename vtype> using iSpinMatrix                  = iScalar<iMatrix<iScalar<vtype>, Ns> >;
+    template<typename vtype> using iColourMatrix                = iScalar<iScalar<iMatrix<vtype, Nc> > > ;
+    template<typename vtype> using iSpinColourMatrix            = iScalar<iMatrix<iMatrix<vtype, Nc>, Ns> >;
+    template<typename vtype> using iLorentzColourMatrix         = iVector<iScalar<iMatrix<vtype, Nc> >, Nd > ;
+    template<typename vtype> using iDoubleStoredColourMatrix    = iVector<iScalar<iMatrix<vtype, Nc> >, Nds > ;
+    template<typename vtype> using iSpinVector                  = iScalar<iVector<iScalar<vtype>, Ns> >;
+    template<typename vtype> using iColourVector                = iScalar<iScalar<iVector<vtype, Nc> > >;
+    template<typename vtype> using iSpinColourVector            = iScalar<iVector<iVector<vtype, Nc>, Ns> >;
+    template<typename vtype> using iHalfSpinVector              = iScalar<iVector<iScalar<vtype>, Nhs> >;
+    template<typename vtype> using iHalfSpinColourVector        = iScalar<iVector<iVector<vtype, Nc>, Nhs> >;
+    template<typename vtype> using iSpinColourSpinColourMatrix  = iScalar<iMatrix<iMatrix<iMatrix<iMatrix<vtype, Nc>, Ns>, Nc>, Ns> >;
+
 
     template<typename vtype> using iGparitySpinColourVector       = iVector<iVector<iVector<vtype, Nc>, Ns>, Ngp >;
     template<typename vtype> using iGparityHalfSpinColourVector   = iVector<iVector<iVector<vtype, Nc>, Nhs>, Ngp >;
@@ -127,10 +130,28 @@ namespace QCD {
     typedef iSpinColourMatrix<Complex  >    SpinColourMatrix;
     typedef iSpinColourMatrix<ComplexF >    SpinColourMatrixF;
     typedef iSpinColourMatrix<ComplexD >    SpinColourMatrixD;
-
+    
     typedef iSpinColourMatrix<vComplex >    vSpinColourMatrix;
     typedef iSpinColourMatrix<vComplexF>    vSpinColourMatrixF;
     typedef iSpinColourMatrix<vComplexD>    vSpinColourMatrixD;
+    
+    // SpinColourSpinColour matrix
+    typedef iSpinColourSpinColourMatrix<Complex  >    SpinColourSpinColourMatrix;
+    typedef iSpinColourSpinColourMatrix<ComplexF >    SpinColourSpinColourMatrixF;
+    typedef iSpinColourSpinColourMatrix<ComplexD >    SpinColourSpinColourMatrixD;
+
+    typedef iSpinColourSpinColourMatrix<vComplex >    vSpinColourSpinColourMatrix;
+    typedef iSpinColourSpinColourMatrix<vComplexF>    vSpinColourSpinColourMatrixF;
+    typedef iSpinColourSpinColourMatrix<vComplexD>    vSpinColourSpinColourMatrixD;
+
+    // SpinColourSpinColour matrix
+    typedef iSpinColourSpinColourMatrix<Complex  >    SpinColourSpinColourMatrix;
+    typedef iSpinColourSpinColourMatrix<ComplexF >    SpinColourSpinColourMatrixF;
+    typedef iSpinColourSpinColourMatrix<ComplexD >    SpinColourSpinColourMatrixD;
+
+    typedef iSpinColourSpinColourMatrix<vComplex >    vSpinColourSpinColourMatrix;
+    typedef iSpinColourSpinColourMatrix<vComplexF>    vSpinColourSpinColourMatrixF;
+    typedef iSpinColourSpinColourMatrix<vComplexD>    vSpinColourSpinColourMatrixD;
 
     // LorentzColour
     typedef iLorentzColourMatrix<Complex  > LorentzColourMatrix;
@@ -229,6 +250,9 @@ namespace QCD {
     typedef Lattice<vSpinColourMatrixF>     LatticeSpinColourMatrixF;
     typedef Lattice<vSpinColourMatrixD>     LatticeSpinColourMatrixD;
 
+    typedef Lattice<vSpinColourSpinColourMatrix>      LatticeSpinColourSpinColourMatrix;
+    typedef Lattice<vSpinColourSpinColourMatrixF>     LatticeSpinColourSpinColourMatrixF;
+    typedef Lattice<vSpinColourSpinColourMatrixD>     LatticeSpinColourSpinColourMatrixD;
 
     typedef Lattice<vLorentzColourMatrix>  LatticeLorentzColourMatrix;
     typedef Lattice<vLorentzColourMatrixF> LatticeLorentzColourMatrixF;
diff --git a/Grid/qcd/action/fermion/CayleyFermion5D.cc b/Grid/qcd/action/fermion/CayleyFermion5D.cc
index d4ed5a7c..a95ea4a0 100644
--- a/Grid/qcd/action/fermion/CayleyFermion5D.cc
+++ b/Grid/qcd/action/fermion/CayleyFermion5D.cc
@@ -68,6 +68,26 @@ void CayleyFermion5D<Impl>::ExportPhysicalFermionSolution(const FermionField &so
   ExtractSlice(exported4d, tmp, 0, 0);
 }
 template<class Impl>  
+void CayleyFermion5D<Impl>::P(const FermionField &psi, FermionField &chi)
+{
+  int Ls= this->Ls;
+  chi=zero;
+  for(int s=0;s<Ls;s++){
+    axpby_ssp_pminus(chi,1.0,chi,1.0,psi,s,s);
+    axpby_ssp_pplus (chi,1.0,chi,1.0,psi,s,(s+1)%Ls);
+  }
+}
+template<class Impl>  
+void CayleyFermion5D<Impl>::Pdag(const FermionField &psi, FermionField &chi)
+{
+  int Ls= this->Ls;
+  chi=zero;
+  for(int s=0;s<Ls;s++){
+    axpby_ssp_pminus(chi,1.0,chi,1.0,psi,s,s);
+    axpby_ssp_pplus (chi,1.0,chi,1.0,psi,s,(s-1+Ls)%Ls);
+  }
+}
+template<class Impl>  
 void CayleyFermion5D<Impl>::ExportPhysicalFermionSource(const FermionField &solution5d,FermionField &exported4d)
 {
   int Ls = this->Ls;
@@ -465,9 +485,13 @@ void CayleyFermion5D<Impl>::SetCoefficientsInternal(RealD zolo_hi,std::vector<Co
     
   double bpc = b+c;
   double bmc = b-c;
+  _b = b;
+  _c = c;
+  _gamma  = gamma; // Save the parameters so we can change mass later.
+  _zolo_hi= zolo_hi;
   for(int i=0; i < Ls; i++){
     as[i] = 1.0;
-    omega[i] = gamma[i]*zolo_hi; //NB reciprocal relative to Chroma NEF code
+    omega[i] = _gamma[i]*_zolo_hi; //NB reciprocal relative to Chroma NEF code
     assert(omega[i]!=Coeff_t(0.0));
     bs[i] = 0.5*(bpc/omega[i] + bmc);
     cs[i] = 0.5*(bpc/omega[i] - bmc);
diff --git a/Grid/qcd/action/fermion/CayleyFermion5D.h b/Grid/qcd/action/fermion/CayleyFermion5D.h
index 3bf2a8b6..203a015e 100644
--- a/Grid/qcd/action/fermion/CayleyFermion5D.h
+++ b/Grid/qcd/action/fermion/CayleyFermion5D.h
@@ -93,6 +93,17 @@ namespace Grid {
       virtual void ImportPhysicalFermionSource(const FermionField &input4d,FermionField &imported5d);
       virtual void ImportUnphysicalFermion(const FermionField &solution5d, FermionField &exported4d);
 
+      ///////////////////////////////////////////////////////////////
+      // Support for MADWF tricks
+      ///////////////////////////////////////////////////////////////
+      RealD Mass(void) { return mass; };
+      void  SetMass(RealD _mass) { 
+	mass=_mass; 
+	SetCoefficientsInternal(_zolo_hi,_gamma,_b,_c);  // Reset coeffs
+      } ;
+      void  P(const FermionField &psi, FermionField &chi);
+      void  Pdag(const FermionField &psi, FermionField &chi);
+
       /////////////////////////////////////////////////////
       // Instantiate different versions depending on Impl
       /////////////////////////////////////////////////////
@@ -139,6 +150,12 @@ namespace Grid {
       //    protected:
       RealD mass;
 
+      // Save arguments to SetCoefficientsInternal
+      std::vector<Coeff_t> _gamma;
+      RealD                _zolo_hi;
+      RealD                _b;
+      RealD                _c;
+
       // Cayley form Moebius (tanh and zolotarev)
       std::vector<Coeff_t> omega;
       std::vector<Coeff_t> bs;    // S dependent coeffs
diff --git a/Grid/qcd/action/fermion/Fermion.h b/Grid/qcd/action/fermion/Fermion.h
index 2a008cb7..77a4681f 100644
--- a/Grid/qcd/action/fermion/Fermion.h
+++ b/Grid/qcd/action/fermion/Fermion.h
@@ -80,12 +80,24 @@ Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
 ///////////////////////////////////////////////////////////////////////////////
 #include <Grid/qcd/action/fermion/g5HermitianLinop.h>
 
+///////////////////////////////////////////////////////////////////////////////
+// Fourier accelerated Pauli Villars inverse support
+///////////////////////////////////////////////////////////////////////////////
+#include <Grid/qcd/action/fermion/WilsonTMFermion5D.h>   
+
+////////////////////////////////////////////////////////////////////////////////
+// Move this group to a DWF specific tools/algorithms subdir? 
+////////////////////////////////////////////////////////////////////////////////
+#include <Grid/qcd/action/fermion/FourierAcceleratedPV.h>
+#include <Grid/qcd/action/fermion/PauliVillarsInverters.h>
+#include <Grid/qcd/action/fermion/Reconstruct5Dprop.h>
+#include <Grid/qcd/action/fermion/MADWF.h>
+
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 // More maintainable to maintain the following typedef list centrally, as more "impl" targets
 // are added, (e.g. extension for gparity, half precision project in comms etc..)
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 
-
 // Cayley 5d
 namespace Grid {
   namespace QCD {
diff --git a/Grid/qcd/action/fermion/FermionOperatorImpl.h b/Grid/qcd/action/fermion/FermionOperatorImpl.h
index e6f6e136..b0ffa90e 100644
--- a/Grid/qcd/action/fermion/FermionOperatorImpl.h
+++ b/Grid/qcd/action/fermion/FermionOperatorImpl.h
@@ -141,6 +141,7 @@ namespace QCD {
   ////////////////////////////////////////////////////////////////////////
   
 #define INHERIT_FIMPL_TYPES(Impl)\
+  typedef Impl Impl_t;							\
   typedef typename Impl::FermionField           FermionField;		\
   typedef typename Impl::PropagatorField     PropagatorField;		\
   typedef typename Impl::DoubledGaugeField DoubledGaugeField;		\
diff --git a/Grid/qcd/action/fermion/FourierAcceleratedPV.h b/Grid/qcd/action/fermion/FourierAcceleratedPV.h
new file mode 100644
index 00000000..d6196eee
--- /dev/null
+++ b/Grid/qcd/action/fermion/FourierAcceleratedPV.h
@@ -0,0 +1,237 @@
+
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/qcd/action/fermion/FourierAcceleratedPV.h
+
+    Copyright (C) 2015
+
+Author: Christoph Lehner (lifted with permission by Peter Boyle, brought back to Grid)
+Author: Peter Boyle <pabobyle@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 */
+#pragma once
+namespace Grid {
+namespace QCD {
+
+  template<typename M>
+    void get_real_const_bc(M& m, RealD& _b, RealD& _c) {
+    ComplexD b,c;
+    b=m.bs[0];
+    c=m.cs[0];
+    std::cout << GridLogMessage << "b=" << b << ", c=" << c << std::endl;
+    for (size_t i=1;i<m.bs.size();i++) {
+      assert(m.bs[i] == b);
+      assert(m.cs[i] == c);
+    }
+    assert(b.imag() == 0.0);
+    assert(c.imag() == 0.0);
+    _b = b.real();
+    _c = c.real();
+  }
+
+
+template<typename Vi, typename M, typename G>
+class FourierAcceleratedPV {
+ public:
+
+  ConjugateGradient<Vi> &cg;
+  M& dwfPV;
+  G& Umu;
+  GridCartesian* grid5D;
+  GridRedBlackCartesian* gridRB5D;
+  int group_in_s;
+
+  FourierAcceleratedPV(M& _dwfPV, G& _Umu, ConjugateGradient<Vi> &_cg, int _group_in_s = 2) 
+   : dwfPV(_dwfPV), Umu(_Umu), cg(_cg), group_in_s(_group_in_s) 
+  {
+    assert( dwfPV.FermionGrid()->_fdimensions[0] % (2*group_in_s) == 0);
+    grid5D = QCD::SpaceTimeGrid::makeFiveDimGrid(2*group_in_s, (GridCartesian*)Umu._grid);
+    gridRB5D = QCD::SpaceTimeGrid::makeFiveDimRedBlackGrid(2*group_in_s, (GridCartesian*)Umu._grid);
+  }
+
+  void rotatePV(const Vi& _src, Vi& dst, bool forward) const {
+
+    GridStopWatch gsw1, gsw2;
+
+    typedef typename Vi::scalar_type Coeff_t;
+    int Ls = dst._grid->_fdimensions[0];
+
+    Vi _tmp(dst._grid);
+    double phase = M_PI / (double)Ls;
+    Coeff_t bzero(0.0,0.0);
+
+    FFT theFFT((GridCartesian*)dst._grid);
+
+    if (!forward) {
+      gsw1.Start();
+      for (int s=0;s<Ls;s++) {
+	Coeff_t a(::cos(phase*s),-::sin(phase*s));
+	axpby_ssp(_tmp,a,_src,bzero,_src,s,s);
+      }
+      gsw1.Stop();
+
+      gsw2.Start();
+      theFFT.FFT_dim(dst,_tmp,0,FFT::forward);
+      gsw2.Stop();
+
+    } else {
+
+      gsw2.Start();
+      theFFT.FFT_dim(_tmp,_src,0,FFT::backward);
+      gsw2.Stop();
+
+      gsw1.Start();
+      for (int s=0;s<Ls;s++) {
+	Coeff_t a(::cos(phase*s),::sin(phase*s));
+	axpby_ssp(dst,a,_tmp,bzero,_tmp,s,s);
+      }
+      gsw1.Stop();
+    }
+
+    std::cout << GridLogMessage << "Timing rotatePV: " << gsw1.Elapsed() << ", " << gsw2.Elapsed() << std::endl;
+
+  }
+
+  void pvInv(const Vi& _src, Vi& _dst) const {
+
+    std::cout << GridLogMessage << "Fourier-Accelerated Outer Pauli Villars"<<std::endl;
+
+    typedef typename Vi::scalar_type Coeff_t;
+    int Ls = _dst._grid->_fdimensions[0];
+
+    GridStopWatch gswT;
+    gswT.Start();
+
+    RealD b,c;
+    get_real_const_bc(dwfPV,b,c);
+    RealD M5 = dwfPV.M5;
+    
+    // U(true) Rightinv TMinv U(false) = Minv
+
+    Vi _src_diag(_dst._grid);
+    Vi _src_diag_slice(dwfPV.GaugeGrid());
+    Vi _dst_diag_slice(dwfPV.GaugeGrid());
+    Vi _src_diag_slices(grid5D);
+    Vi _dst_diag_slices(grid5D);
+    Vi _dst_diag(_dst._grid);
+
+    rotatePV(_src,_src_diag,false);
+
+    // now do TM solves
+    Gamma G5(Gamma::Algebra::Gamma5);
+
+    GridStopWatch gswA, gswB;
+
+    gswA.Start();
+
+    typedef typename M::Impl_t Impl;
+    //WilsonTMFermion<Impl> tm(x.Umu,*x.UGridF,*x.UrbGridF,0.0,0.0,solver_outer.parent.par.wparams_f);
+    std::vector<RealD> vmass(grid5D->_fdimensions[0],0.0);
+    std::vector<RealD> vmu(grid5D->_fdimensions[0],0.0);
+
+    WilsonTMFermion5D<Impl> tm(Umu,*grid5D,*gridRB5D,
+			   *(GridCartesian*)dwfPV.GaugeGrid(),
+			   *(GridRedBlackCartesian*)dwfPV.GaugeRedBlackGrid(),
+			   vmass,vmu);
+    
+    //SchurRedBlackDiagTwoSolve<Vi> sol(cg);
+    SchurRedBlackDiagMooeeSolve<Vi> sol(cg); // same performance as DiagTwo
+    gswA.Stop();
+
+    gswB.Start();
+
+    for (int sgroup=0;sgroup<Ls/2/group_in_s;sgroup++) {
+
+      for (int sidx=0;sidx<group_in_s;sidx++) {
+
+	int s = sgroup*group_in_s + sidx;
+	int sprime = Ls-s-1;
+
+	RealD phase = M_PI / (RealD)Ls * (2.0 * s + 1.0);
+	RealD cosp = ::cos(phase);
+	RealD sinp = ::sin(phase);
+	RealD denom = b*b + c*c + 2.0*b*c*cosp;
+	RealD mass = -(b*b*M5 + c*(1.0 - cosp + c*M5) + b*(-1.0 + cosp + 2.0*c*cosp*M5))/denom;
+	RealD mu = (b+c)*sinp/denom;
+
+	vmass[2*sidx + 0] = mass;
+	vmass[2*sidx + 1] = mass;
+	vmu[2*sidx + 0] = mu;
+	vmu[2*sidx + 1] = -mu;
+
+      }
+
+      tm.update(vmass,vmu);
+
+      for (int sidx=0;sidx<group_in_s;sidx++) {
+
+	int s = sgroup*group_in_s + sidx;
+	int sprime = Ls-s-1;
+
+	ExtractSlice(_src_diag_slice,_src_diag,s,0);
+	InsertSlice(_src_diag_slice,_src_diag_slices,2*sidx + 0,0);
+
+	ExtractSlice(_src_diag_slice,_src_diag,sprime,0);
+	InsertSlice(_src_diag_slice,_src_diag_slices,2*sidx + 1,0);
+
+      }
+
+      GridStopWatch gsw;
+      gsw.Start();
+      _dst_diag_slices = zero; // zero guess
+      sol(tm,_src_diag_slices,_dst_diag_slices);
+      gsw.Stop();
+      std::cout << GridLogMessage << "Solve[sgroup=" << sgroup << "] completed in " << gsw.Elapsed() << ", " << gswA.Elapsed() << std::endl;
+
+      for (int sidx=0;sidx<group_in_s;sidx++) {
+
+	int s = sgroup*group_in_s + sidx;
+	int sprime = Ls-s-1;
+
+	RealD phase = M_PI / (RealD)Ls * (2.0 * s + 1.0);
+	RealD cosp = ::cos(phase);
+	RealD sinp = ::sin(phase);
+
+	// now rotate with inverse of
+	Coeff_t pA = b + c*cosp;
+	Coeff_t pB = - Coeff_t(0.0,1.0)*c*sinp;
+	Coeff_t pABden = pA*pA - pB*pB;
+	// (pA + pB * G5) * (pA - pB*G5) = (pA^2 - pB^2)
+      
+	ExtractSlice(_dst_diag_slice,_dst_diag_slices,2*sidx + 0,0);
+	_dst_diag_slice = (pA/pABden) * _dst_diag_slice - (pB/pABden) * (G5 * _dst_diag_slice);
+	InsertSlice(_dst_diag_slice,_dst_diag,s,0);
+	
+	ExtractSlice(_dst_diag_slice,_dst_diag_slices,2*sidx + 1,0);
+	_dst_diag_slice = (pA/pABden) * _dst_diag_slice + (pB/pABden) * (G5 * _dst_diag_slice);
+	InsertSlice(_dst_diag_slice,_dst_diag,sprime,0);
+      }
+    }
+    gswB.Stop();
+
+    rotatePV(_dst_diag,_dst,true);
+
+    gswT.Stop();
+    std::cout << GridLogMessage << "PV completed in " << gswT.Elapsed() << " (Setup: " << gswA.Elapsed() << ", s-loop: " << gswB.Elapsed() << ")" << std::endl;
+  }
+
+};
+}}
diff --git a/Grid/qcd/action/fermion/MADWF.h b/Grid/qcd/action/fermion/MADWF.h
new file mode 100644
index 00000000..064b13a8
--- /dev/null
+++ b/Grid/qcd/action/fermion/MADWF.h
@@ -0,0 +1,193 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/algorithms/iterative/MADWF.h
+
+    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 */
+#pragma once
+
+namespace Grid {
+namespace QCD {
+
+template <class Fieldi, class Fieldo,IfNotSame<Fieldi,Fieldo> X=0>
+inline void convert(const Fieldi &from,Fieldo &to) 
+{
+  precisionChange(to,from);
+}
+template <class Fieldi, class Fieldo,IfSame<Fieldi,Fieldo> X=0>
+inline void convert(const Fieldi &from,Fieldo &to) 
+{
+  to=from;
+}
+
+template<class Matrixo,class Matrixi,class PVinverter,class SchurSolver, class Guesser> 
+class MADWF 
+{
+ private:
+  typedef typename Matrixo::FermionField FermionFieldo;
+  typedef typename Matrixi::FermionField FermionFieldi;
+
+  PVinverter  & PauliVillarsSolvero;// For the outer field
+  SchurSolver & SchurSolveri;       // For the inner approx field
+  Guesser     & Guesseri;           // To deflate the inner approx solves
+
+  Matrixo & Mato;                   // Action object for outer
+  Matrixi & Mati;                   // Action object for inner
+
+  RealD target_resid;
+  int   maxiter;
+ public:
+
+  MADWF(Matrixo &_Mato,
+	Matrixi &_Mati, 
+	PVinverter &_PauliVillarsSolvero, 
+	SchurSolver &_SchurSolveri,
+	Guesser & _Guesseri,
+	RealD resid,
+	int _maxiter) :
+
+  Mato(_Mato),Mati(_Mati),
+    SchurSolveri(_SchurSolveri),
+    PauliVillarsSolvero(_PauliVillarsSolvero),Guesseri(_Guesseri)
+  {   
+    target_resid=resid;
+    maxiter     =_maxiter; 
+  };
+
+  void operator() (const FermionFieldo &src4,FermionFieldo &sol5)
+  {
+    std::cout << GridLogMessage<< " ************************************************" << std::endl;
+    std::cout << GridLogMessage<< "  MADWF-like algorithm                           " << std::endl;
+    std::cout << GridLogMessage<< " ************************************************" << std::endl;
+
+    FermionFieldi    c0i(Mati.GaugeGrid()); // 4d 
+    FermionFieldi    y0i(Mati.GaugeGrid()); // 4d
+    FermionFieldo    c0 (Mato.GaugeGrid()); // 4d 
+    FermionFieldo    y0 (Mato.GaugeGrid()); // 4d
+
+    FermionFieldo    A(Mato.FermionGrid()); // Temporary outer
+    FermionFieldo    B(Mato.FermionGrid()); // Temporary outer
+    FermionFieldo    b(Mato.FermionGrid()); // 5d source
+
+    FermionFieldo    c(Mato.FermionGrid()); // PVinv source; reused so store
+    FermionFieldo    defect(Mato.FermionGrid()); // 5d source
+
+    FermionFieldi   ci(Mati.FermionGrid()); 
+    FermionFieldi   yi(Mati.FermionGrid()); 
+    FermionFieldi   xi(Mati.FermionGrid()); 
+    FermionFieldi srci(Mati.FermionGrid()); 
+    FermionFieldi   Ai(Mati.FermionGrid()); 
+
+    RealD m=Mati.Mass();
+
+    ///////////////////////////////////////
+    //Import source, include Dminus factors
+    ///////////////////////////////////////
+    Mato.ImportPhysicalFermionSource(src4,b); 
+    std::cout << GridLogMessage << " src4 " <<norm2(src4)<<std::endl;
+    std::cout << GridLogMessage << " b    " <<norm2(b)<<std::endl;
+
+    defect = b;
+    sol5=zero;
+    for (int i=0;i<maxiter;i++) {
+
+      ///////////////////////////////////////
+      // Set up c0 from current defect
+      ///////////////////////////////////////
+      PauliVillarsSolvero(Mato,defect,A);
+      Mato.Pdag(A,c);
+      ExtractSlice(c0, c, 0 , 0);
+
+      ////////////////////////////////////////////////
+      // Solve the inner system with surface term c0
+      ////////////////////////////////////////////////
+      ci = zero;  
+      convert(c0,c0i); // Possible precison change
+      InsertSlice(c0i,ci,0, 0);
+
+      // Dwm P y = Dwm x = D(1) P (c0,0,0,0)^T
+      Mati.P(ci,Ai);
+      Mati.SetMass(1.0);      Mati.M(Ai,srci);      Mati.SetMass(m);
+      SchurSolveri(Mati,srci,xi,Guesseri); 
+      Mati.Pdag(xi,yi);
+      ExtractSlice(y0i, yi, 0 , 0);
+      convert(y0i,y0); // Possible precision change
+
+      //////////////////////////////////////
+      // Propagate solution back to outer system
+      // Build Pdag PV^-1 Dm P [-sol4,c2,c3... cL]
+      //////////////////////////////////////
+      c0 = - y0;
+      InsertSlice(c0, c, 0   , 0);
+
+      /////////////////////////////
+      // Reconstruct the bulk solution Pdag PV^-1 Dm P 
+      /////////////////////////////
+      Mato.P(c,B);
+      Mato.M(B,A);
+      PauliVillarsSolvero(Mato,A,B);
+      Mato.Pdag(B,A);
+
+      //////////////////////////////
+      // Reinsert surface prop
+      //////////////////////////////
+      InsertSlice(y0,A,0,0);
+
+      //////////////////////////////
+      // Convert from y back to x 
+      //////////////////////////////
+      Mato.P(A,B);
+
+      //         sol5' = sol5 + M^-1 defect
+      //               = sol5 + M^-1 src - M^-1 M sol5  ...
+      sol5 = sol5 + B;
+      std::cout << GridLogMessage << "***************************************" <<std::endl;
+      std::cout << GridLogMessage << " Sol5 update "<<std::endl;
+      std::cout << GridLogMessage << "***************************************" <<std::endl;
+      std::cout << GridLogMessage << " Sol5 now "<<norm2(sol5)<<std::endl;
+      std::cout << GridLogMessage << " delta    "<<norm2(B)<<std::endl;
+
+       // New defect  = b - M sol5
+       Mato.M(sol5,A);
+       defect = b - A;
+
+       std::cout << GridLogMessage << " defect   "<<norm2(defect)<<std::endl;
+
+       double resid = ::sqrt(norm2(defect) / norm2(b));
+       std::cout << GridLogMessage << "Residual " << i << ": " << resid  << std::endl;
+       std::cout << GridLogMessage << "***************************************" <<std::endl;
+
+       if (resid < target_resid) {
+	 return;
+       }
+    }
+
+    std::cout << GridLogMessage << "MADWF : Exceeded maxiter "<<std::endl;
+    assert(0);
+
+  }
+
+};
+
+}}
diff --git a/Grid/qcd/action/fermion/PauliVillarsInverters.h b/Grid/qcd/action/fermion/PauliVillarsInverters.h
new file mode 100644
index 00000000..b10640e3
--- /dev/null
+++ b/Grid/qcd/action/fermion/PauliVillarsInverters.h
@@ -0,0 +1,95 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/algorithms/iterative/SchurRedBlack.h
+
+    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 */
+#pragma once
+
+namespace Grid {
+namespace QCD {
+
+template<class Field>
+class PauliVillarsSolverUnprec
+{
+ public:
+  ConjugateGradient<Field> & CG;
+  PauliVillarsSolverUnprec(  ConjugateGradient<Field> &_CG) : CG(_CG){};
+
+  template<class Matrix>
+  void operator() (Matrix &_Matrix,const Field &src,Field &sol)
+  {
+    RealD m = _Matrix.Mass();
+    Field A  (_Matrix.FermionGrid());
+
+    MdagMLinearOperator<Matrix,Field> HermOp(_Matrix);
+
+    _Matrix.SetMass(1.0);
+    _Matrix.Mdag(src,A);
+    CG(HermOp,A,sol);
+    _Matrix.SetMass(m);
+  };
+};
+
+template<class Field,class SchurSolverType>
+class PauliVillarsSolverRBprec
+{
+ public:
+  SchurSolverType & SchurSolver;
+  PauliVillarsSolverRBprec( SchurSolverType &_SchurSolver) : SchurSolver(_SchurSolver){};
+
+  template<class Matrix>
+  void operator() (Matrix &_Matrix,const Field &src,Field &sol)
+  {
+    RealD m = _Matrix.Mass();
+    Field A  (_Matrix.FermionGrid());
+
+    _Matrix.SetMass(1.0);
+    SchurSolver(_Matrix,src,sol);
+    _Matrix.SetMass(m);
+  };
+};
+
+template<class Field,class GaugeField>
+class PauliVillarsSolverFourierAccel
+{
+ public:
+  GaugeField      & Umu;
+  ConjugateGradient<Field> & CG;
+
+  PauliVillarsSolverFourierAccel(GaugeField &_Umu,ConjugateGradient<Field> &_CG) :  Umu(_Umu), CG(_CG)
+  {
+  };
+
+  template<class Matrix>
+  void operator() (Matrix &_Matrix,const Field &src,Field &sol)
+  {
+    FourierAcceleratedPV<Field, Matrix, typename Matrix::GaugeField > faPV(_Matrix,Umu,CG) ;
+    faPV.pvInv(src,sol);
+  };
+};
+
+
+}
+}
diff --git a/Grid/qcd/action/fermion/Reconstruct5Dprop.h b/Grid/qcd/action/fermion/Reconstruct5Dprop.h
new file mode 100644
index 00000000..6862c5ee
--- /dev/null
+++ b/Grid/qcd/action/fermion/Reconstruct5Dprop.h
@@ -0,0 +1,135 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/algorithms/iterative/SchurRedBlack.h
+
+    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 */
+#pragma once
+
+namespace Grid {
+namespace QCD {
+
+template<class Field,class PVinverter> class Reconstruct5DfromPhysical {
+ private:
+  PVinverter & PauliVillarsSolver;
+ public:
+
+ /////////////////////////////////////////////////////
+ // First cut works, 10 Oct 2018.
+ //
+ // Must form a plan to get this into production for Zmobius acceleration
+ // of the Mobius exact AMA corrections.
+ //
+ // TODO : understand absence of contact term in eqns in Hantao's thesis
+ //        sol4 is contact term subtracted, but thesis & Brower's paper suggests not.
+ //
+ // Step 1: Localise PV inverse in a routine. [DONE]
+ // Step 2: Schur based PV inverse            [DONE]
+ // Step 3: Fourier accelerated PV inverse    [DONE]
+ //
+ /////////////////////////////////////////////////////
+ 
+  Reconstruct5DfromPhysical(PVinverter &_PauliVillarsSolver) 
+    : PauliVillarsSolver(_PauliVillarsSolver) 
+  { 
+  };
+
+
+   template<class Matrix>
+   void PV(Matrix &_Matrix,const Field &src,Field &sol)
+   {
+     RealD m = _Matrix.Mass();
+     _Matrix.SetMass(1.0);
+     _Matrix.M(src,sol);
+     _Matrix.SetMass(m);
+   }
+   template<class Matrix>
+   void PVdag(Matrix &_Matrix,const Field &src,Field &sol)
+   {
+     RealD m = _Matrix.Mass();
+     _Matrix.SetMass(1.0);
+     _Matrix.Mdag(src,sol);
+     _Matrix.SetMass(m);
+   }
+  template<class Matrix>
+  void operator() (Matrix & _Matrix,const Field &sol4,const Field &src4, Field &sol5){
+
+    int Ls =  _Matrix.Ls;
+
+    Field psi4(_Matrix.GaugeGrid());
+    Field psi(_Matrix.FermionGrid());
+    Field A  (_Matrix.FermionGrid());
+    Field B  (_Matrix.FermionGrid());
+    Field c  (_Matrix.FermionGrid());
+
+    typedef typename Matrix::Coeff_t Coeff_t;
+
+    std::cout << GridLogMessage<< " ************************************************" << std::endl;
+    std::cout << GridLogMessage<< " Reconstruct5Dprop: c.f. MADWF algorithm         " << std::endl;
+    std::cout << GridLogMessage<< " ************************************************" << std::endl;
+
+    ///////////////////////////////////////
+    //Import source, include Dminus factors
+    ///////////////////////////////////////
+    _Matrix.ImportPhysicalFermionSource(src4,B); 
+
+    ///////////////////////////////////////
+    // Set up c from src4
+    ///////////////////////////////////////
+    PauliVillarsSolver(_Matrix,B,A);
+    _Matrix.Pdag(A,c);
+
+    //////////////////////////////////////
+    // Build Pdag PV^-1 Dm P [-sol4,c2,c3... cL]
+    //////////////////////////////////////
+    psi4 = - sol4;
+    InsertSlice(psi4, psi, 0   , 0);
+    for (int s=1;s<Ls;s++) {
+      ExtractSlice(psi4,c,s,0);
+       InsertSlice(psi4,psi,s,0);
+    }
+
+    /////////////////////////////
+    // Pdag PV^-1 Dm P 
+    /////////////////////////////
+    _Matrix.P(psi,B);
+    _Matrix.M(B,A);
+    PauliVillarsSolver(_Matrix,A,B);
+    _Matrix.Pdag(B,A);
+
+    //////////////////////////////
+    // Reinsert surface prop
+    //////////////////////////////
+    InsertSlice(sol4,A,0,0);
+
+    //////////////////////////////
+    // Convert from y back to x 
+    //////////////////////////////
+    _Matrix.P(A,sol5);
+    
+  }
+};
+
+}
+}
diff --git a/Grid/qcd/action/fermion/WilsonTMFermion5D.h b/Grid/qcd/action/fermion/WilsonTMFermion5D.h
new file mode 100644
index 00000000..6c75d908
--- /dev/null
+++ b/Grid/qcd/action/fermion/WilsonTMFermion5D.h
@@ -0,0 +1,155 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./lib/qcd/action/fermion/WilsonTMFermion5D.h
+
+    Copyright (C) 2015
+
+Author: paboyle <paboyle@ph.ed.ac.uk> ; NB Christoph did similar in GPT
+
+    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 */
+#pragma once 
+
+#include <Grid/qcd/action/fermion/FermionCore.h>
+#include <Grid/qcd/action/fermion/WilsonFermion.h>
+
+
+namespace Grid {
+
+  namespace QCD {
+    
+    template<class Impl>
+      class WilsonTMFermion5D : public WilsonFermion5D<Impl>
+      {
+      public:
+	INHERIT_IMPL_TYPES(Impl);
+      public:
+
+	virtual void   Instantiatable(void) {};
+
+	// Constructors
+        WilsonTMFermion5D(GaugeField &_Umu,
+			  GridCartesian         &Fgrid,
+			  GridRedBlackCartesian &Frbgrid, 
+			  GridCartesian         &Ugrid,
+			  GridRedBlackCartesian &Urbgrid, 
+			  const std::vector<RealD> _mass,
+			  const std::vector<RealD> _mu,
+			  const ImplParams &p= ImplParams()
+			  ) :
+	WilsonFermion5D<Impl>(_Umu,
+			      Fgrid,
+			      Frbgrid,
+			      Ugrid,
+			      Urbgrid,
+			      4.0,p)
+	
+	  {
+	    update(_mass,_mu);
+	  }
+
+	virtual void Meooe(const FermionField &in, FermionField &out) {
+	  if (in.checkerboard == Odd) {
+	    this->DhopEO(in, out, DaggerNo);
+	  } else {
+	    this->DhopOE(in, out, DaggerNo);
+	  }
+	}
+
+	virtual void MeooeDag(const FermionField &in, FermionField &out) {
+	  if (in.checkerboard == Odd) {
+	    this->DhopEO(in, out, DaggerYes);
+	  } else {
+	    this->DhopOE(in, out, DaggerYes);
+	  }
+	}	
+	
+	// allow override for twisted mass and clover
+	virtual void Mooee(const FermionField &in, FermionField &out) {
+	  out.checkerboard = in.checkerboard;
+	  //axpibg5x(out,in,a,b); // out = a*in + b*i*G5*in
+	  for (int s=0;s<(int)this->mass.size();s++) {
+	    ComplexD a = 4.0+this->mass[s];
+	    ComplexD b(0.0,this->mu[s]);
+	    axpbg5y_ssp(out,a,in,b,in,s,s);
+	  }
+	}
+
+	virtual void MooeeDag(const FermionField &in, FermionField &out) {
+	  out.checkerboard = in.checkerboard;
+	  for (int s=0;s<(int)this->mass.size();s++) {
+	    ComplexD a = 4.0+this->mass[s];
+	    ComplexD b(0.0,-this->mu[s]);
+	    axpbg5y_ssp(out,a,in,b,in,s,s);
+	  }
+	}
+	virtual void MooeeInv(const FermionField &in, FermionField &out) {
+	  for (int s=0;s<(int)this->mass.size();s++) {
+	    RealD m    = this->mass[s];
+	    RealD tm   = this->mu[s];
+	    RealD mtil = 4.0+this->mass[s];
+	    RealD sq   = mtil*mtil+tm*tm;
+	    ComplexD a    = mtil/sq;
+	    ComplexD b(0.0, -tm /sq);
+	    axpbg5y_ssp(out,a,in,b,in,s,s);
+	  }
+	}
+	virtual void MooeeInvDag(const FermionField &in, FermionField &out) {
+	  for (int s=0;s<(int)this->mass.size();s++) {
+	    RealD m    = this->mass[s];
+	    RealD tm   = this->mu[s];
+	    RealD mtil = 4.0+this->mass[s];
+	    RealD sq   = mtil*mtil+tm*tm;
+	    ComplexD a    = mtil/sq;
+	    ComplexD b(0.0,tm /sq);
+	    axpbg5y_ssp(out,a,in,b,in,s,s);
+	  }
+	}
+
+	virtual RealD M(const FermionField &in, FermionField &out) {
+	  out.checkerboard = in.checkerboard;
+	  this->Dhop(in, out, DaggerNo);
+	  FermionField tmp(out._grid);
+	  for (int s=0;s<(int)this->mass.size();s++) {
+	    ComplexD a = 4.0+this->mass[s];
+	    ComplexD b(0.0,this->mu[s]);
+	    axpbg5y_ssp(tmp,a,in,b,in,s,s);
+	  }
+	  return axpy_norm(out, 1.0, tmp, out);
+	}
+	
+	// needed for fast PV
+	void update(const std::vector<RealD>& _mass, const std::vector<RealD>& _mu) {
+	  assert(_mass.size() == _mu.size());
+	  assert(_mass.size() == this->FermionGrid()->_fdimensions[0]);
+	  this->mass = _mass;
+	  this->mu = _mu;
+	}
+	
+      private:
+	std::vector<RealD> mu;
+	std::vector<RealD> mass;
+	
+      };
+   
+    typedef WilsonTMFermion5D<WilsonImplF> WilsonTMFermion5DF; 
+    typedef WilsonTMFermion5D<WilsonImplD> WilsonTMFermion5DD; 
+
+}}
diff --git a/Hadrons/A2AVectors.hpp b/Hadrons/A2AVectors.hpp
index 423a5f08..ff32ff5e 100644
--- a/Hadrons/A2AVectors.hpp
+++ b/Hadrons/A2AVectors.hpp
@@ -36,7 +36,7 @@ See the full license in the file "LICENSE" in the top level distribution directo
 BEGIN_HADRONS_NAMESPACE
 
 /******************************************************************************
- *               Classes to generate V & W all-to-all vectors                 *
+ *                 Class to generate V & W all-to-all vectors                 *
  ******************************************************************************/
 template <typename FImpl>
 class A2AVectorsSchurDiagTwo
@@ -70,6 +70,42 @@ private:
     SchurDiagTwoOperator<FMat, FermionField> op_;
 };
 
+/******************************************************************************
+ *                  Methods for V & W all-to-all vectors I/O                  *
+ ******************************************************************************/
+class A2AVectorsIo
+{
+public:
+    struct Record: Serializable
+    {
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Record,
+                                        unsigned int, index);
+        Record(void): index(0) {}
+    };
+public:
+    template <typename Field>
+    static void write(const std::string fileStem, std::vector<Field> &vec, 
+                      const bool multiFile, const int trajectory = -1);
+    template <typename Field>
+    static void read(std::vector<Field> &vec, const std::string fileStem,
+                     const bool multiFile, const int trajectory = -1);
+private:
+    static inline std::string vecFilename(const std::string stem, const int traj, 
+                                          const bool multiFile)
+    {
+        std::string t = (traj < 0) ? "" : ("." + std::to_string(traj));
+
+        if (multiFile)
+        {
+            return stem + t;
+        }
+        else
+        {
+            return stem + t + ".bin";
+        }
+    }
+};
+
 /******************************************************************************
  *               A2AVectorsSchurDiagTwo template implementation               *
  ******************************************************************************/
@@ -217,6 +253,90 @@ void A2AVectorsSchurDiagTwo<FImpl>::makeHighModeW5D(FermionField &wout_4d,
     }
 }
 
+/******************************************************************************
+ *               all-to-all vectors I/O template implementation               *
+ ******************************************************************************/
+template <typename Field>
+void A2AVectorsIo::write(const std::string fileStem, std::vector<Field> &vec, 
+                         const bool multiFile, const int trajectory)
+{
+    Record       record;
+    GridBase     *grid = vec[0]._grid;
+    ScidacWriter binWriter(grid->IsBoss());
+    std::string  filename = vecFilename(fileStem, multiFile, trajectory);
+
+    if (multiFile)
+    {
+        std::string fullFilename;
+
+        for (unsigned int i = 0; i < vec.size(); ++i)
+        {
+            fullFilename = filename + "/elem" + std::to_string(i) + ".bin";
+
+            LOG(Message) << "Writing vector " << i << std::endl;
+            makeFileDir(fullFilename, grid);
+            binWriter.open(fullFilename);
+            record.index = i;
+            binWriter.writeScidacFieldRecord(vec[i], record);
+            binWriter.close();
+        }
+    }
+    else
+    {
+        makeFileDir(filename, grid);
+        binWriter.open(filename);
+        for (unsigned int i = 0; i < vec.size(); ++i)
+        {
+            LOG(Message) << "Writing vector " << i << std::endl;
+            record.index = i;
+            binWriter.writeScidacFieldRecord(vec[i], record);
+        }
+        binWriter.close();
+    }
+}
+
+template <typename Field>
+void A2AVectorsIo::read(std::vector<Field> &vec, const std::string fileStem, 
+                        const bool multiFile, const int trajectory)
+{
+    Record       record;
+    ScidacReader binReader;
+    std::string  filename = vecFilename(fileStem, multiFile, trajectory);
+
+    if (multiFile)
+    {
+        std::string fullFilename;
+
+        for (unsigned int i = 0; i < vec.size(); ++i)
+        {
+            fullFilename = filename + "/elem" + std::to_string(i) + ".bin";
+
+            LOG(Message) << "Reading vector " << i << std::endl;
+            binReader.open(fullFilename);
+            binReader.readScidacFieldRecord(vec[i], record);
+            binReader.close();
+            if (record.index != i)
+            {
+                HADRONS_ERROR(Io, "vector index mismatch");
+            }
+        }
+    }
+    else
+    {
+        binReader.open(filename);
+        for (unsigned int i = 0; i < vec.size(); ++i)
+        {
+            LOG(Message) << "Reading vector " << i << std::endl;
+            binReader.readScidacFieldRecord(vec[i], record);
+            if (record.index != i)
+            {
+                HADRONS_ERROR(Io, "vector index mismatch");
+            }
+        }
+        binReader.close();
+    }
+}
+
 END_HADRONS_NAMESPACE
 
 #endif // A2A_Vectors_hpp_
diff --git a/Hadrons/DiskVector.hpp b/Hadrons/DiskVector.hpp
index f3023a6e..1fddbe36 100644
--- a/Hadrons/DiskVector.hpp
+++ b/Hadrons/DiskVector.hpp
@@ -151,6 +151,12 @@ class EigenDiskVector: public DiskVectorBase<EigenDiskVectorMat<T>>
 public:
     using DiskVectorBase<EigenDiskVectorMat<T>>::DiskVectorBase;
     typedef EigenDiskVectorMat<T> Matrix;
+public:
+    T operator()(const unsigned int i, const Eigen::Index j,
+                 const Eigen::Index k) const
+    {
+        return (*this)[i](j, k);
+    }
 private:
     virtual void load(EigenDiskVectorMat<T> &obj, const std::string filename) const
     {
diff --git a/Hadrons/Modules.hpp b/Hadrons/Modules.hpp
index 6d56ea2d..a866c5cd 100644
--- a/Hadrons/Modules.hpp
+++ b/Hadrons/Modules.hpp
@@ -16,6 +16,7 @@
 #include <Hadrons/Modules/MSource/Wall.hpp>
 #include <Hadrons/Modules/MSource/Z2.hpp>
 #include <Hadrons/Modules/MSource/SeqConserved.hpp>
+#include <Hadrons/Modules/MSource/Momentum.hpp>
 #include <Hadrons/Modules/MSink/Smear.hpp>
 #include <Hadrons/Modules/MSink/Point.hpp>
 #include <Hadrons/Modules/MSolver/MixedPrecisionRBPrecCG.hpp>
@@ -27,6 +28,7 @@
 #include <Hadrons/Modules/MGauge/StoutSmearing.hpp>
 #include <Hadrons/Modules/MGauge/Unit.hpp>
 #include <Hadrons/Modules/MGauge/Random.hpp>
+#include <Hadrons/Modules/MGauge/GaugeFix.hpp>
 #include <Hadrons/Modules/MGauge/FundtoHirep.hpp>
 #include <Hadrons/Modules/MGauge/StochEm.hpp>
 #include <Hadrons/Modules/MNoise/TimeDilutedSpinColorDiagonal.hpp>
@@ -40,6 +42,9 @@
 #include <Hadrons/Modules/MScalar/ScalarVP.hpp>
 #include <Hadrons/Modules/MScalar/Scalar.hpp>
 #include <Hadrons/Modules/MScalar/ChargedProp.hpp>
+#include <Hadrons/Modules/MNPR/Bilinear.hpp>
+#include <Hadrons/Modules/MNPR/Amputate.hpp>
+#include <Hadrons/Modules/MNPR/FourQuark.hpp>
 #include <Hadrons/Modules/MAction/DWF.hpp>
 #include <Hadrons/Modules/MAction/MobiusDWF.hpp>
 #include <Hadrons/Modules/MAction/Wilson.hpp>
@@ -60,6 +65,7 @@
 #include <Hadrons/Modules/MScalarSUN/TrKinetic.hpp>
 #include <Hadrons/Modules/MIO/LoadEigenPack.hpp>
 #include <Hadrons/Modules/MIO/LoadNersc.hpp>
+#include <Hadrons/Modules/MIO/LoadA2AVectors.hpp>
 #include <Hadrons/Modules/MIO/LoadCosmHol.hpp>
 #include <Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp>
 #include <Hadrons/Modules/MIO/LoadBinary.hpp>
diff --git a/Hadrons/Modules/MAction/DWF.cc b/Hadrons/Modules/MAction/DWF.cc
index 6cda84ac..38d25cb9 100644
--- a/Hadrons/Modules/MAction/DWF.cc
+++ b/Hadrons/Modules/MAction/DWF.cc
@@ -32,4 +32,6 @@ using namespace Hadrons;
 using namespace MAction;
 
 template class Grid::Hadrons::MAction::TDWF<FIMPL>;
+#ifdef GRID_DEFAULT_PRECISION_DOUBLE
 template class Grid::Hadrons::MAction::TDWF<FIMPLF>;
+#endif
diff --git a/Hadrons/Modules/MAction/DWF.hpp b/Hadrons/Modules/MAction/DWF.hpp
index 67cfeb1b..257782a1 100644
--- a/Hadrons/Modules/MAction/DWF.hpp
+++ b/Hadrons/Modules/MAction/DWF.hpp
@@ -73,7 +73,9 @@ protected:
 };
 
 MODULE_REGISTER_TMP(DWF, TDWF<FIMPL>, MAction);
+#ifdef GRID_DEFAULT_PRECISION_DOUBLE
 MODULE_REGISTER_TMP(DWFF, TDWF<FIMPLF>, MAction);
+#endif
 
 /******************************************************************************
  *                        DWF template implementation                         *
diff --git a/Hadrons/Modules/MAction/MobiusDWF.cc b/Hadrons/Modules/MAction/MobiusDWF.cc
index 41683920..879452d8 100644
--- a/Hadrons/Modules/MAction/MobiusDWF.cc
+++ b/Hadrons/Modules/MAction/MobiusDWF.cc
@@ -32,4 +32,6 @@ using namespace Hadrons;
 using namespace MAction;
 
 template class Grid::Hadrons::MAction::TMobiusDWF<FIMPL>;
+#ifdef GRID_DEFAULT_PRECISION_DOUBLE
 template class Grid::Hadrons::MAction::TMobiusDWF<FIMPLF>;
+#endif
diff --git a/Hadrons/Modules/MAction/MobiusDWF.hpp b/Hadrons/Modules/MAction/MobiusDWF.hpp
index 22964c9a..01223267 100644
--- a/Hadrons/Modules/MAction/MobiusDWF.hpp
+++ b/Hadrons/Modules/MAction/MobiusDWF.hpp
@@ -72,7 +72,9 @@ public:
 };
 
 MODULE_REGISTER_TMP(MobiusDWF, TMobiusDWF<FIMPL>, MAction);
+#ifdef GRID_DEFAULT_PRECISION_DOUBLE
 MODULE_REGISTER_TMP(MobiusDWFF, TMobiusDWF<FIMPLF>, MAction);
+#endif
 
 /******************************************************************************
  *                      TMobiusDWF implementation                             *
diff --git a/Hadrons/Modules/MAction/ScaledDWF.cc b/Hadrons/Modules/MAction/ScaledDWF.cc
index 70eafac5..7008bf5d 100644
--- a/Hadrons/Modules/MAction/ScaledDWF.cc
+++ b/Hadrons/Modules/MAction/ScaledDWF.cc
@@ -32,4 +32,6 @@ using namespace Hadrons;
 using namespace MAction;
 
 template class Grid::Hadrons::MAction::TScaledDWF<FIMPL>;
+#ifdef GRID_DEFAULT_PRECISION_DOUBLE
 template class Grid::Hadrons::MAction::TScaledDWF<FIMPLF>;
+#endif
diff --git a/Hadrons/Modules/MAction/ScaledDWF.hpp b/Hadrons/Modules/MAction/ScaledDWF.hpp
index c890f0e9..f7a09eef 100644
--- a/Hadrons/Modules/MAction/ScaledDWF.hpp
+++ b/Hadrons/Modules/MAction/ScaledDWF.hpp
@@ -71,7 +71,9 @@ public:
 };
 
 MODULE_REGISTER_TMP(ScaledDWF, TScaledDWF<FIMPL>, MAction);
+#ifdef GRID_DEFAULT_PRECISION_DOUBLE
 MODULE_REGISTER_TMP(ScaledDWFF, TScaledDWF<FIMPLF>, MAction);
+#endif
 
 /******************************************************************************
  *                      TScaledDWF implementation                             *
diff --git a/Hadrons/Modules/MAction/Wilson.cc b/Hadrons/Modules/MAction/Wilson.cc
index 4ce3e7ef..1e801ed6 100644
--- a/Hadrons/Modules/MAction/Wilson.cc
+++ b/Hadrons/Modules/MAction/Wilson.cc
@@ -32,4 +32,6 @@ using namespace Hadrons;
 using namespace MAction;
 
 template class Grid::Hadrons::MAction::TWilson<FIMPL>;
+#ifdef GRID_DEFAULT_PRECISION_DOUBLE
 template class Grid::Hadrons::MAction::TWilson<FIMPLF>;
+#endif
diff --git a/Hadrons/Modules/MAction/Wilson.hpp b/Hadrons/Modules/MAction/Wilson.hpp
index a9327c2f..093449e6 100644
--- a/Hadrons/Modules/MAction/Wilson.hpp
+++ b/Hadrons/Modules/MAction/Wilson.hpp
@@ -71,7 +71,9 @@ protected:
 };
 
 MODULE_REGISTER_TMP(Wilson, TWilson<FIMPL>, MAction);
+#ifdef GRID_DEFAULT_PRECISION_DOUBLE
 MODULE_REGISTER_TMP(WilsonF, TWilson<FIMPLF>, MAction);
+#endif
 
 /******************************************************************************
  *                     TWilson template implementation                        *
diff --git a/Hadrons/Modules/MAction/WilsonClover.cc b/Hadrons/Modules/MAction/WilsonClover.cc
index 2c5c0e66..eed1582c 100644
--- a/Hadrons/Modules/MAction/WilsonClover.cc
+++ b/Hadrons/Modules/MAction/WilsonClover.cc
@@ -32,4 +32,6 @@ using namespace Hadrons;
 using namespace MAction;
 
 template class Grid::Hadrons::MAction::TWilsonClover<FIMPL>;
+#ifdef GRID_DEFAULT_PRECISION_DOUBLE
 template class Grid::Hadrons::MAction::TWilsonClover<FIMPLF>;
+#endif
diff --git a/Hadrons/Modules/MAction/WilsonClover.hpp b/Hadrons/Modules/MAction/WilsonClover.hpp
index 349abe84..0b78bb55 100644
--- a/Hadrons/Modules/MAction/WilsonClover.hpp
+++ b/Hadrons/Modules/MAction/WilsonClover.hpp
@@ -75,7 +75,9 @@ public:
 };
 
 MODULE_REGISTER_TMP(WilsonClover, TWilsonClover<FIMPL>, MAction);
+#ifdef GRID_DEFAULT_PRECISION_DOUBLE
 MODULE_REGISTER_TMP(WilsonCloverF, TWilsonClover<FIMPLF>, MAction);
+#endif
 
 /******************************************************************************
  *                    TWilsonClover template implementation                   *
diff --git a/Hadrons/Modules/MAction/ZMobiusDWF.cc b/Hadrons/Modules/MAction/ZMobiusDWF.cc
index ef8e4799..609b76cc 100644
--- a/Hadrons/Modules/MAction/ZMobiusDWF.cc
+++ b/Hadrons/Modules/MAction/ZMobiusDWF.cc
@@ -32,4 +32,6 @@ using namespace Hadrons;
 using namespace MAction;
 
 template class Grid::Hadrons::MAction::TZMobiusDWF<ZFIMPL>;
+#ifdef GRID_DEFAULT_PRECISION_DOUBLE
 template class Grid::Hadrons::MAction::TZMobiusDWF<ZFIMPLF>;
+#endif
diff --git a/Hadrons/Modules/MAction/ZMobiusDWF.hpp b/Hadrons/Modules/MAction/ZMobiusDWF.hpp
index f7959127..40b04d7f 100644
--- a/Hadrons/Modules/MAction/ZMobiusDWF.hpp
+++ b/Hadrons/Modules/MAction/ZMobiusDWF.hpp
@@ -73,7 +73,9 @@ public:
 };
 
 MODULE_REGISTER_TMP(ZMobiusDWF, TZMobiusDWF<ZFIMPL>, MAction);
+#ifdef GRID_DEFAULT_PRECISION_DOUBLE
 MODULE_REGISTER_TMP(ZMobiusDWFF, TZMobiusDWF<ZFIMPLF>, MAction);
+#endif
 
 /******************************************************************************
  *                     TZMobiusDWF implementation                             *
diff --git a/Hadrons/Modules/MContraction/A2AAslashField.cc b/Hadrons/Modules/MContraction/A2AAslashField.cc
index eb76932e..65c49198 100644
--- a/Hadrons/Modules/MContraction/A2AAslashField.cc
+++ b/Hadrons/Modules/MContraction/A2AAslashField.cc
@@ -1,3 +1,30 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MContraction/A2AAslashField.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+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 <Hadrons/Modules/MContraction/A2AAslashField.hpp>
 
 using namespace Grid;
diff --git a/Hadrons/Modules/MContraction/A2AAslashField.hpp b/Hadrons/Modules/MContraction/A2AAslashField.hpp
index a1d64d3d..8b99692b 100644
--- a/Hadrons/Modules/MContraction/A2AAslashField.hpp
+++ b/Hadrons/Modules/MContraction/A2AAslashField.hpp
@@ -1,3 +1,30 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MContraction/A2AAslashField.hpp
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+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 */
 #ifndef Hadrons_MContraction_A2AAslashField_hpp_
 #define Hadrons_MContraction_A2AAslashField_hpp_
 
diff --git a/Hadrons/Modules/MGauge/GaugeFix.cc b/Hadrons/Modules/MGauge/GaugeFix.cc
new file mode 100644
index 00000000..53aa16da
--- /dev/null
+++ b/Hadrons/Modules/MGauge/GaugeFix.cc
@@ -0,0 +1,36 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MGauge/GaugeFix.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+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 <Hadrons/Modules/MGauge/GaugeFix.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MGauge;
+
+template class Grid::Hadrons::MGauge::TGaugeFix<GIMPL>;
diff --git a/Hadrons/Modules/MGauge/GaugeFix.hpp b/Hadrons/Modules/MGauge/GaugeFix.hpp
new file mode 100644
index 00000000..ece8c19d
--- /dev/null
+++ b/Hadrons/Modules/MGauge/GaugeFix.hpp
@@ -0,0 +1,135 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MGauge/GaugeFix.hpp
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+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 */
+
+#ifndef Hadrons_MGaugeFix_hpp_
+#define Hadrons_MGaugeFix_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+#include <Grid/qcd/utils/GaugeFix.h>
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                              Fix gauge                                    *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MGauge)
+
+class GaugeFixPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(GaugeFixPar,
+                                    std::string, gauge,
+                                    Real,  alpha,
+                                    int, maxiter, 
+                                    Real, Omega_tol, 
+                                    Real, Phi_tol,
+                                    bool, Fourier);
+};
+
+template <typename GImpl>
+class TGaugeFix: public Module<GaugeFixPar>
+{
+public:
+    GAUGE_TYPE_ALIASES(GImpl,);
+public:
+    // constructor
+    TGaugeFix(const std::string name);
+    // destructor
+    virtual ~TGaugeFix(void) {};
+    // dependencies/products
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_TMP(GaugeFix, TGaugeFix<GIMPL>, MGauge);
+
+/******************************************************************************
+*                            TGaugeFix implementation                             *
+******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename GImpl>
+TGaugeFix<GImpl>::TGaugeFix(const std::string name)
+: Module<GaugeFixPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename GImpl>
+std::vector<std::string> TGaugeFix<GImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().gauge};
+    return in;
+}
+
+template <typename GImpl>
+std::vector<std::string> TGaugeFix<GImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename GImpl>
+void TGaugeFix<GImpl>::setup(void)
+{
+    envCreateLat(GaugeField, getName());
+}
+
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename GImpl>
+void TGaugeFix<GImpl>::execute(void)
+//Loads the gauge and fixes it
+{
+    std::cout << "executing" << std::endl;
+    LOG(Message) << "Fixing the Gauge" << std::endl;
+    LOG(Message) << par().gauge << std::endl;
+    auto &U     = envGet(GaugeField, par().gauge);
+    auto &Umu   = envGet(GaugeField, getName());
+    LOG(Message) << "Gauge Field fetched" << std::endl;
+    //do we allow maxiter etc to be user set?
+    Real alpha     = par().alpha;
+    int  maxiter   = par().maxiter;
+    Real Omega_tol = par().Omega_tol;
+    Real Phi_tol   = par().Phi_tol;
+    bool Fourier   = par().Fourier;
+    FourierAcceleratedGaugeFixer<PeriodicGimplR>::SteepestDescentGaugeFix(U,alpha,maxiter,Omega_tol,Phi_tol,Fourier);
+    Umu = U;
+    LOG(Message) << "Gauge Fixed" << std::endl;
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MGaugeFix_hpp_
diff --git a/Hadrons/Modules/MIO/LoadA2AVectors.cc b/Hadrons/Modules/MIO/LoadA2AVectors.cc
new file mode 100644
index 00000000..7a40a6f5
--- /dev/null
+++ b/Hadrons/Modules/MIO/LoadA2AVectors.cc
@@ -0,0 +1,34 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MIO/LoadA2AVectors.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+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 <Hadrons/Modules/MIO/LoadA2AVectors.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MIO;
+
+template class Grid::Hadrons::MIO::TLoadA2AVectors<FIMPL>;
diff --git a/Hadrons/Modules/MIO/LoadA2AVectors.hpp b/Hadrons/Modules/MIO/LoadA2AVectors.hpp
new file mode 100644
index 00000000..5b194c16
--- /dev/null
+++ b/Hadrons/Modules/MIO/LoadA2AVectors.hpp
@@ -0,0 +1,120 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MIO/LoadA2AVectors.hpp
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+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 */
+#ifndef Hadrons_MIO_LoadA2AVectors_hpp_
+#define Hadrons_MIO_LoadA2AVectors_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+#include <Hadrons/A2AVectors.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                    Module to load all-to-all vectors                       *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MIO)
+
+class LoadA2AVectorsPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(LoadA2AVectorsPar,
+                                    std::string,  filestem,
+                                    bool,         multiFile,
+                                    unsigned int, size);
+};
+
+template <typename FImpl>
+class TLoadA2AVectors: public Module<LoadA2AVectorsPar>
+{
+public:
+    FERM_TYPE_ALIASES(FImpl,);
+public:
+    // constructor
+    TLoadA2AVectors(const std::string name);
+    // destructor
+    virtual ~TLoadA2AVectors(void) {};
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_TMP(LoadA2AVectors, TLoadA2AVectors<FIMPL>, MIO);
+
+/******************************************************************************
+ *                      TLoadA2AVectors implementation                        *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TLoadA2AVectors<FImpl>::TLoadA2AVectors(const std::string name)
+: Module<LoadA2AVectorsPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TLoadA2AVectors<FImpl>::getInput(void)
+{
+    std::vector<std::string> in;
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TLoadA2AVectors<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TLoadA2AVectors<FImpl>::setup(void)
+{
+    envCreate(std::vector<FermionField>, getName(), 1, par().size, 
+              envGetGrid(FermionField));
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TLoadA2AVectors<FImpl>::execute(void)
+{
+    auto &vec = envGet(std::vector<FermionField>, getName());
+
+    A2AVectorsIo::read(vec, par().filestem, par().multiFile, vm().getTrajectory());
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MIO_LoadA2AVectors_hpp_
diff --git a/Hadrons/Modules/MIO/LoadEigenPack.cc b/Hadrons/Modules/MIO/LoadEigenPack.cc
index f0ae63c3..28fdeb01 100644
--- a/Hadrons/Modules/MIO/LoadEigenPack.cc
+++ b/Hadrons/Modules/MIO/LoadEigenPack.cc
@@ -32,4 +32,6 @@ using namespace Hadrons;
 using namespace MIO;
 
 template class Grid::Hadrons::MIO::TLoadEigenPack<FermionEigenPack<FIMPL>>;
+#ifdef GRID_DEFAULT_PRECISION_DOUBLE
 template class Grid::Hadrons::MIO::TLoadEigenPack<FermionEigenPack<FIMPL, FIMPLF>>;
+#endif
diff --git a/Hadrons/Modules/MIO/LoadEigenPack.hpp b/Hadrons/Modules/MIO/LoadEigenPack.hpp
index 2302b15d..016675c9 100644
--- a/Hadrons/Modules/MIO/LoadEigenPack.hpp
+++ b/Hadrons/Modules/MIO/LoadEigenPack.hpp
@@ -72,7 +72,9 @@ public:
 };
 
 MODULE_REGISTER_TMP(LoadFermionEigenPack, TLoadEigenPack<FermionEigenPack<FIMPL>>, MIO);
+#ifdef GRID_DEFAULT_PRECISION_DOUBLE
 MODULE_REGISTER_TMP(LoadFermionEigenPackIo32, ARG(TLoadEigenPack<FermionEigenPack<FIMPL, FIMPLF>>), MIO);
+#endif
 
 /******************************************************************************
  *                    TLoadEigenPack implementation                           *
diff --git a/Hadrons/Modules/MNPR/Amputate.cc b/Hadrons/Modules/MNPR/Amputate.cc
new file mode 100644
index 00000000..ec7c5940
--- /dev/null
+++ b/Hadrons/Modules/MNPR/Amputate.cc
@@ -0,0 +1,36 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MNPR/Amputate.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+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 <Hadrons/Modules/MNPR/Amputate.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MNPR;
+
+template class Grid::Hadrons::MNPR::TAmputate<FIMPL,FIMPL>;
+
diff --git a/Hadrons/Modules/MNPR/Amputate.hpp b/Hadrons/Modules/MNPR/Amputate.hpp
new file mode 100644
index 00000000..953bb354
--- /dev/null
+++ b/Hadrons/Modules/MNPR/Amputate.hpp
@@ -0,0 +1,200 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MNPR/Amputate.hpp
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+Author: Julia Kettle J.R.Kettle-2@sms.ed.ac.uk
+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 */
+
+#ifndef Hadrons_Amputate_hpp_
+#define Hadrons_Amputate_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+#include <Grid/Eigen/LU>
+//#include <Grid/qcd/utils/PropagatorUtils.h>
+//#include <Grid/serialisation/Serialisation.h>
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                                TAmputate                                       *
+        Performs bilinear contractions of the type tr[g5*adj(Sout)*g5*G*Sin]
+        Suitable for non exceptional momenta
+******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MNPR)
+
+class AmputatePar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(AmputatePar,
+                                    std::string,    Sin, //need to make this a propogator type?
+                                    std::string,    Sout, //same
+                                    std::string,    vertex,
+                                    std::string,    pin,
+                                    std::string,    pout,
+                                    std::string,    output,
+                                    std::string,    input);
+};
+
+template <typename FImpl1, typename FImpl2>
+class TAmputate: public Module<AmputatePar>
+{
+public:
+    FERM_TYPE_ALIASES(FImpl1, 1);
+    FERM_TYPE_ALIASES(FImpl2, 2);
+    class Result: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
+                                        std::vector<Complex>, Vamp,
+                                        ); 
+    };
+public:
+    // constructor
+    TAmputate(const std::string name);
+    // destructor
+    virtual ~TAmputate(void) {};
+    // dependencies/products
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    virtual SpinColourMatrix invertspincolmat(SpinColourMatrix &scmat);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_TMP(Amputate, ARG(TAmputate<FIMPL, FIMPL>), MNPR);
+
+/******************************************************************************
+ *                           TAmputate implementation                            *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2>
+TAmputate<FImpl1, FImpl2>::TAmputate(const std::string name)
+: Module<AmputatePar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2>
+std::vector<std::string> TAmputate<FImpl1, FImpl2>::getInput(void)
+{
+    std::vector<std::string> input = {par().Sin, par().Sout, par().vertex};
+    
+    return input;
+}
+
+template <typename FImpl1, typename FImpl2>
+std::vector<std::string> TAmputate<FImpl1, FImpl2>::getOutput(void)
+{
+    std::vector<std::string> output = {getName()};
+    
+    
+    return output;
+}
+
+// Invert spin colour matrix using Eigen
+template <typename Fimpl1, typename Fimpl2>
+SpinColourMatrix TAmputate<Fimpl1, Fimpl2>::invertspincolmat(SpinColourMatrix &scmat)
+{
+    Eigen::MatrixXcf scmat_2d(Ns*Nc,Ns*Nc);
+    for(int ic=0; ic<Nc; ic++){
+    for(int jc=0; jc<Nc; jc++){
+        for(int is=0; is<Ns; is++){
+        for(int js=0; js<Ns; js++){
+            scmat_2d(Ns*ic+is,Ns*jc+js) = scmat()(is,js)(ic,jc);
+        }}
+    }}      
+    Eigen::MatrixXcf scmat_2d_inv = scmat_2d.inverse();
+    SpinColourMatrix scmat_inv;
+    for(int ic=0; ic<Nc; ic++){
+    for(int jc=0; jc<Nc; jc++){
+        for(int is=0; is<Ns; is++){
+        for(int js=0; js<Ns; js++){
+            scmat_inv()(is,js)(ic,jc) = scmat_2d_inv(Ns*ic+is,Ns*jc+js);
+        }}
+    }}      
+    return scmat_inv;
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2>
+void TAmputate<FImpl1, FImpl2>::execute(void)
+{
+    LOG(Message) << "Computing bilinear amputations '" << getName() << "' using"
+                 << " momentum '" << par().Sin << "' and '" << par().Sout << "'"
+                 << std::endl;
+    BinaryWriter                    writer(par().output);
+    PropagatorField1                &Sin = *env().template getObject<PropagatorField1>(par().Sin); //Do these have the phases taken into account?? Don't think so. FIX
+    PropagatorField2                &Sout = *env().template getObject<PropagatorField2>(par().Sout);
+    std::vector<int>                pin  = strToVec<int>(par().pin), pout = strToVec<int>(par().pout);
+    std::vector<Real>               latt_size(pin.begin(), pin.end()); 
+    LatticeComplex                  pdotxin(env().getGrid()), pdotxout(env().getGrid()), coor(env().getGrid());
+    LOG(Message) << "Propagators set up " << std::endl;
+    std::vector<SpinColourMatrix>   vertex; // Let's read from file here
+    Gamma                           g5(Gamma::Algebra::Gamma5);
+    Result                          result;
+    LOG(Message) << "reading file - "  << par().input << std::endl;
+    BinaryReader                    reader(par().input); 
+    Complex                         Ci(0.0,1.0);
+
+    std::string svertex;
+    read(reader,"vertex", vertex);
+    LOG(Message) << "vertex read" << std::endl;
+
+    pdotxin=zero;
+    pdotxout=zero;
+    for (unsigned int mu = 0; mu < 4; ++mu)
+    {
+        Real TwoPiL =  M_PI * 2.0/ latt_size[mu];
+        LatticeCoordinate(coor,mu);
+        pdotxin = pdotxin +(TwoPiL * pin[mu]) * coor;
+        pdotxout= pdotxout +(TwoPiL * pout[mu]) * coor;
+    }
+    Sin = Sin*exp(-Ci*pdotxin); //phase corrections
+    Sout = Sout*exp(-Ci*pdotxout);
+
+    SpinColourMatrix Sin_mom = sum(Sin);
+    SpinColourMatrix Sout_mom = sum(Sout);
+    LOG(Message) << "summed over lattice" << std::endl;
+   
+    LOG(Message) << "Lattice -> spincolourmatrix conversion" << std::endl;
+
+    SpinColourMatrix Sin_inv = invertspincolmat(Sin_mom);
+    SpinColourMatrix Sout_inv = invertspincolmat(Sout_mom);
+    LOG(Message) << "Inversions done" << std::endl;
+
+    result.Vamp.resize(Gamma::nGamma/2);
+    for( int mu=0; mu < Gamma::nGamma/2; mu++){
+        Gamma::Algebra gam = mu;
+        result.Vamp[mu] = 1/12.0*trace(adj(Gamma(mu*2+1))*g5*Sout_inv*g5*vertex[mu]*Sin_inv);
+        LOG(Message) << "Vamp[" << mu << "] - " << result.Vamp[mu] << std::endl;
+        }
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Amputate_hpp_
diff --git a/Hadrons/Modules/MNPR/Bilinear.cc b/Hadrons/Modules/MNPR/Bilinear.cc
new file mode 100644
index 00000000..c5b38d37
--- /dev/null
+++ b/Hadrons/Modules/MNPR/Bilinear.cc
@@ -0,0 +1,36 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MNPR/Bilinear.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+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 <Hadrons/Modules/MNPR/Bilinear.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MNPR;
+
+template class Grid::Hadrons::MNPR::TBilinear<FIMPL,FIMPL>;
+
diff --git a/Hadrons/Modules/MNPR/Bilinear.hpp b/Hadrons/Modules/MNPR/Bilinear.hpp
new file mode 100644
index 00000000..66cd22a6
--- /dev/null
+++ b/Hadrons/Modules/MNPR/Bilinear.hpp
@@ -0,0 +1,225 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MNPR/Bilinear.hpp
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+Author: Julia Kettle J.R.Kettle-2@sms.ed.ac.uk
+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 */
+
+#ifndef Hadrons_Bilinear_hpp_
+#define Hadrons_Bilinear_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+//#include <Grid/qcd/utils/PropagatorUtils.h>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                                TBilinear                                       *
+        Performs bilinear contractions of the type tr[g5*adj(Sout)*g5*G*Sin]
+        Suitable for non exceptional momenta in Rome-Southampton NPR
+******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MNPR)
+
+class BilinearPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(BilinearPar,
+                                    std::string,    Sin,
+                                    std::string,    Sout,
+                                    std::string,    pin,
+                                    std::string,    pout,
+                                    std::string,    output);
+};
+
+template <typename FImpl1, typename FImpl2>
+class TBilinear: public Module<BilinearPar>
+{
+public:
+    FERM_TYPE_ALIASES(FImpl1, 1);
+    FERM_TYPE_ALIASES(FImpl2, 2);
+    class Result: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Result, 
+                                        std::vector<SpinColourMatrix>, bilinear);
+    };
+public:
+    // constructor
+    TBilinear(const std::string name);
+    // destructor
+    virtual ~TBilinear(void) {};
+    // dependencies/products
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    //LatticeSpinColourMatrix PhaseProps(LatticeSpinColourMatrix S, std::vector<Real> p);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_TMP(Bilinear, ARG(TBilinear<FIMPL, FIMPL>), MNPR);
+
+/******************************************************************************
+ *                           TBilinear implementation                            *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2>
+TBilinear<FImpl1, FImpl2>::TBilinear(const std::string name)
+: Module<BilinearPar>(name)
+{}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2>
+void TBilinear<FImpl1, FImpl2>::setup(void)
+{
+    //env().template registerLattice<LatticeSpinColourMatrix>(getName());
+    //env().template registerObject<SpinColourMatrix>(getName());
+}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2>
+std::vector<std::string> TBilinear<FImpl1, FImpl2>::getInput(void)
+{
+    std::vector<std::string> input = {par().Sin, par().Sout};
+    
+    return input;
+}
+
+template <typename FImpl1, typename FImpl2>
+std::vector<std::string> TBilinear<FImpl1, FImpl2>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+/*
+/////Phase propagators//////////////////////////
+template <typename FImpl1, typename FImpl2>
+LatticeSpinColourMatrix TBilinear<FImpl1, FImpl2>::PhaseProps(LatticeSpinColourMatrix S, std::vector<Real> p)
+{
+    GridBase *grid = S._grid;
+    LatticeComplex      pdotx(grid),  coor(grid);
+    std::vector<int>   latt_size = grid->_fdimensions; 
+    Complex             Ci(0.0,1.0);
+    pdotx=zero;
+    for (unsigned int mu = 0; mu < 4; ++mu)
+    {
+        Real TwoPiL =  M_PI * 2.0/ latt_size[mu];
+        LatticeCoordinate(coor,mu);
+        pdotx = pdotx +(TwoPiL * p[mu]) * coor;
+    }
+    S = S*exp(-Ci*pdotx);
+    return S;
+}
+*/
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2>
+void TBilinear<FImpl1, FImpl2>::execute(void)
+{
+/**************************************************************************
+
+Compute the bilinear vertex needed for the NPR.
+V(G) = sum_x  [ g5 * adj(S'(x,p2)) * g5 * G * S'(x,p1) ]_{si,sj,ci,cj}
+G is one of the 16 gamma vertices [I,gmu,g5,g5gmu,sig(mu,nu)]
+
+        * G
+       / \
+    p1/   \p2
+     /     \
+    /       \
+
+Returns a spin-colour matrix, with indices si,sj, ci,cj
+
+Conventions:
+p1 - incoming momenta
+p2 - outgoing momenta
+q = (p1-p2)
+**************************************************************************/
+
+    LOG(Message) << "Computing bilinear contractions '" << getName() << "' using"
+                 << " momentum '" << par().Sin << "' and '" << par().Sout << "'"
+                 << std::endl;
+     
+    BinaryWriter             writer(par().output);
+    
+
+    // Propogators
+    LatticeSpinColourMatrix     &Sin = *env().template getObject<LatticeSpinColourMatrix>(par().Sin);
+    LatticeSpinColourMatrix     &Sout = *env().template getObject<LatticeSpinColourMatrix>(par().Sout);
+    LatticeComplex              pdotxin(env().getGrid()), pdotxout(env().getGrid()), coor(env().getGrid());
+    // momentum on legs
+    std::vector<Real>           pin  = strToVec<Real>(par().pin), pout = strToVec<Real>(par().pout);
+    std::vector<Real>           latt_size(pin.begin(), pin.end()); 
+    //bilinears
+    LatticeSpinColourMatrix     bilinear_x(env().getGrid());
+    SpinColourMatrix            bilinear;
+    Gamma                       g5(Gamma::Algebra::Gamma5);
+    Result                      result;
+    Complex                     Ci(0.0,1.0);
+
+    //
+
+    pdotxin=zero;
+    pdotxout=zero;
+    for (unsigned int mu = 0; mu < 4; ++mu)
+    {
+        Real TwoPiL =  M_PI * 2.0/ latt_size[mu];
+        LatticeCoordinate(coor,mu);
+        pdotxin = pdotxin +(TwoPiL * pin[mu]) * coor;
+        pdotxout= pdotxout +(TwoPiL * pout[mu]) * coor;
+    }
+    Sin = Sin*exp(-Ci*pdotxin); //phase corrections
+    Sout = Sout*exp(-Ci*pdotxout);
+    
+    ////Set up gamma vector//////////////////////////
+    std::vector<Gamma> gammavector;
+    for( int i=0; i<Gamma::nGamma; i++){
+        Gamma::Algebra gam = i;
+        gammavector.push_back(Gamma(gam));
+    }
+    result.bilinear.resize(Gamma::nGamma);
+    /////////////////////////////////////////////////
+    //LatticeSpinMatrix temp = g5*Sout;
+    ////////Form Vertex//////////////////////////////
+    for (int i=0; i < Gamma::nGamma; i++){
+        bilinear_x = g5*adj(Sout)*g5*gammavector[i]*Sin; 
+        result.bilinear[i] = sum(bilinear_x); //sum over lattice sites
+    }
+    //////////////////////////////////////////////////
+    write(writer, par().output, result.bilinear);
+    LOG(Message) << "Complete. Writing results to " << par().output << std:: endl;
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Bilinear_hpp_
diff --git a/Hadrons/Modules/MNPR/FourQuark.cc b/Hadrons/Modules/MNPR/FourQuark.cc
new file mode 100644
index 00000000..1943c25b
--- /dev/null
+++ b/Hadrons/Modules/MNPR/FourQuark.cc
@@ -0,0 +1,36 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MNPR/FourQuark.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+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 <Hadrons/Modules/MNPR/FourQuark.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MNPR;
+
+template class Grid::Hadrons::MNPR::TFourQuark<FIMPL,FIMPL>;
+
diff --git a/Hadrons/Modules/MNPR/FourQuark.hpp b/Hadrons/Modules/MNPR/FourQuark.hpp
new file mode 100644
index 00000000..f961a366
--- /dev/null
+++ b/Hadrons/Modules/MNPR/FourQuark.hpp
@@ -0,0 +1,274 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MNPR/FourQuark.hpp
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+Author: Julia Kettle J.R.Kettle-2@sms.ed.ac.uk
+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 */
+
+#ifndef Hadrons_FourQuark_hpp_
+#define Hadrons_FourQuark_hpp_
+
+#include <typeinfo>
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+#include <Grid/serialisation/Serialisation.h>
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                                TFourQuark                                       *
+        Performs fourquark contractions of the type tr[g5*adj(Sout)*g5*G*Sin]
+        Suitable for non exceptional momenta
+******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MNPR)
+
+class FourQuarkPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(FourQuarkPar,
+                                    std::string,    Sin, //need to make this a propogator type?
+                                    std::string,    Sout, //same
+                                    std::string,    pin,
+                                    std::string,    pout,
+                                    bool,           fullbasis,
+                                    std::string,    output);
+};
+
+template <typename FImpl1, typename FImpl2>
+class TFourQuark: public Module<FourQuarkPar>
+{
+public:
+    FERM_TYPE_ALIASES(FImpl1, 1);
+    FERM_TYPE_ALIASES(FImpl2, 2);
+    class Result: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
+                                        std::vector<SpinColourSpinColourMatrix>, fourquark);
+    };
+public:
+    // constructor
+    TFourQuark(const std::string name);
+    // destructor
+    virtual ~TFourQuark(void) {};
+    // dependencies/products
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void tensorprod(LatticeSpinColourSpinColourMatrix &lret, LatticeSpinColourMatrix a, LatticeSpinColourMatrix b);
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_TMP(FourQuark, ARG(TFourQuark<FIMPL, FIMPL>), MNPR);
+
+/******************************************************************************
+ *                           TFourQuark implementation                            *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2>
+TFourQuark<FImpl1, FImpl2>::TFourQuark(const std::string name)
+: Module<FourQuarkPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2>
+std::vector<std::string> TFourQuark<FImpl1, FImpl2>::getInput(void)
+{
+    std::vector<std::string> input = {par().Sin, par().Sout};
+    
+    return input;
+}
+
+template <typename FImpl1, typename FImpl2>
+std::vector<std::string> TFourQuark<FImpl1, FImpl2>::getOutput(void)
+{
+    std::vector<std::string> output = {getName()};
+    
+    return output;
+}
+
+
+template <typename FImpl1, typename FImpl2>
+void TFourQuark<FImpl1, FImpl2>::tensorprod(LatticeSpinColourSpinColourMatrix &lret, LatticeSpinColourMatrix a, LatticeSpinColourMatrix b)
+{
+#if 0
+            parallel_for(auto site=lret.begin();site<lret.end();site++) {
+                for (int si; si < 4; ++si){
+                for(int sj; sj <4; ++sj){
+                    for (int ci; ci < 3; ++ci){
+                    for (int cj; cj < 3; ++cj){
+                        for (int sk; sk < 4; ++sk){
+                        for(int sl; sl <4; ++sl){
+                            for (int ck; ck < 3; ++ck){
+                            for (int cl; cl < 3; ++cl){
+                        lret[site]()(si,sj)(ci,cj)(sk,sl)(ck,cl)=a[site]()(si,sj)(ci,cj)*b[site]()(sk,sl)(ck,cl);
+                            }}
+                        }}
+                    }}
+                }}
+        }
+#else 
+            // FIXME ; is there a general need for this construct ? In which case we should encapsulate the
+            //         below loops in a helper function.
+            //LOG(Message) << "sp co mat a is - " << a << std::endl;
+            //LOG(Message) << "sp co mat b is - " << b << std::endl;
+            parallel_for(auto site=lret.begin();site<lret.end();site++) {
+            vTComplex left;
+                for(int si=0; si < Ns; ++si){
+                for(int sj=0; sj < Ns; ++sj){
+                    for (int ci=0; ci < Nc; ++ci){
+                    for (int cj=0; cj < Nc; ++cj){
+                      //LOG(Message) << "si, sj, ci, cj -  " << si << ", " << sj  << ", "<< ci  << ", "<< cj << std::endl;
+                      left()()() = a[site]()(si,sj)(ci,cj);
+                      //LOG(Message) << left << std::endl;
+                      lret[site]()(si,sj)(ci,cj)=left()*b[site]();
+                    }}
+                }}
+            }
+#endif      
+}
+
+
+
+
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2>
+void TFourQuark<FImpl1, FImpl2>::setup(void)
+{
+    envCreateLat(LatticeSpinColourMatrix, getName());
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2>
+void TFourQuark<FImpl1, FImpl2>::execute(void)
+{
+
+/*********************************************************************************
+
+TFourQuark : Creates the four quark vertex required for the NPR of four-quark ops
+
+V_{Gamma_1,Gamma_2} = sum_x [ ( g5 * adj(S'(x,p2)) * g5 * G1 * S'(x,p1) )_ci,cj;si,sj x ( g5 * adj(S'(x,p2)) * g5 * G2 S'(x,p1) )_ck,cl;sk,cl ]
+
+Create a bilinear vertex for G1 and G2  the spin and colour indices are kept free. Where there are 16 potential Gs.
+We then find the outer product of V1 and V2, keeping the spin and colour indices uncontracted
+Then this is summed over the lattice coordinate
+Result is a SpinColourSpinColourMatrix - with 4 colour and 4 spin indices. 
+We have up to 256 of these including the offdiag (G1 != G2).
+
+        \         /
+         \p1   p1/
+          \     /
+           \   /
+         G1 * * G2
+           /   \
+          /     \
+         /p2   p2\
+        /         \
+
+*********************************************************************************/
+
+
+
+
+    LOG(Message) << "Computing fourquark contractions '" << getName() << "' using"
+                 << " momentum '" << par().Sin << "' and '" << par().Sout << "'"
+                 << std::endl;
+    
+    BinaryWriter             writer(par().output);
+    
+    PropagatorField1                            &Sin = *env().template getObject<PropagatorField1>(par().Sin);
+    PropagatorField2                            &Sout = *env().template getObject<PropagatorField2>(par().Sout);
+    std::vector<Real>                           pin  = strToVec<Real>(par().pin), pout = strToVec<Real>(par().pout);
+    bool                                        fullbasis = par().fullbasis;
+    Gamma                                       g5(Gamma::Algebra::Gamma5);
+    Result                                      result;
+    std::vector<Real>                           latt_size(pin.begin(), pin.end());
+    LatticeComplex                              pdotxin(env().getGrid()), pdotxout(env().getGrid()), coor(env().getGrid());
+    LatticeSpinColourMatrix                     bilinear_mu(env().getGrid()), bilinear_nu(env().getGrid());
+    LatticeSpinColourSpinColourMatrix           lret(env().getGrid()); 
+    Complex                         Ci(0.0,1.0);
+
+    //Phase propagators
+    //Sin = Grid::QCD::PropUtils::PhaseProps(Sin,pin);
+    //Sout = Grid::QCD::PropUtils::PhaseProps(Sout,pout);
+    
+    //find p.x for in and out so phase can be accounted for in propagators
+    pdotxin=zero;
+    pdotxout=zero;
+    for (unsigned int mu = 0; mu < 4; ++mu)
+    {
+        Real TwoPiL =  M_PI * 2.0/ latt_size[mu];
+        LatticeCoordinate(coor,mu);
+        pdotxin = pdotxin +(TwoPiL * pin[mu]) * coor;
+        pdotxout= pdotxout +(TwoPiL * pout[mu]) * coor;
+    }
+    Sin = Sin*exp(-Ci*pdotxin); //phase corrections
+    Sout = Sout*exp(-Ci*pdotxout);
+
+
+    //Set up Gammas 
+    std::vector<Gamma> gammavector;
+     for( int i=1; i<Gamma::nGamma; i+=2){
+         Gamma::Algebra gam = i;
+         gammavector.push_back(Gamma(gam));
+       }
+    
+    lret = zero;
+    if (fullbasis == true){ // all combinations of mu and nu
+        result.fourquark.resize(Gamma::nGamma/2*Gamma::nGamma/2);
+        for( int mu=0; mu<Gamma::nGamma/2; mu++){ 
+            bilinear_mu = g5*adj(Sout)*g5*gammavector[mu]*Sin;
+            for ( int nu=0; nu<Gamma::nGamma; nu++){
+                LatticeSpinColourMatrix     bilinear_nu(env().getGrid());
+                bilinear_nu = g5*adj(Sout)*g5*gammavector[nu]*Sin;
+                LOG(Message) << "bilinear_nu for nu = " << nu << " is - " << bilinear_mu << std::endl;
+                result.fourquark[mu*Gamma::nGamma/2 + nu] = zero;
+                tensorprod(lret,bilinear_mu,bilinear_nu);
+                result.fourquark[mu*Gamma::nGamma/2 + nu] = sum(lret);
+            }
+        }
+    } else {
+        result.fourquark.resize(Gamma::nGamma/2);
+        for ( int mu=0; mu<1; mu++){
+        //for( int mu=0; mu<Gamma::nGamma/2; mu++ ){
+            bilinear_mu = g5*adj(Sout)*g5*gammavector[mu]*Sin;
+            //LOG(Message) << "bilinear_mu for mu = " << mu << " is - " << bilinear_mu << std::endl;
+            result.fourquark[mu] = zero;
+            tensorprod(lret,bilinear_mu,bilinear_mu); //tensor outer product
+            result.fourquark[mu] = sum(lret);
+        }
+    }
+    write(writer, "fourquark", result.fourquark);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_FourQuark_hpp_
diff --git a/Hadrons/Modules/MSolver/A2AVectors.hpp b/Hadrons/Modules/MSolver/A2AVectors.hpp
index ee7cfa30..47ae6aee 100644
--- a/Hadrons/Modules/MSolver/A2AVectors.hpp
+++ b/Hadrons/Modules/MSolver/A2AVectors.hpp
@@ -51,7 +51,9 @@ public:
                                   std::string, noise,
                                   std::string, action,
                                   std::string, eigenPack,
-                                  std::string, solver);
+                                  std::string, solver,
+                                  std::string, output,
+                                  bool,        multiFile);
 };
 
 template <typename FImpl, typename Pack>
@@ -236,6 +238,13 @@ void TA2AVectors<FImpl, Pack>::execute(void)
         }
         stopTimer("W high mode");
     }
+
+    // I/O if necessary
+    if (!par().output.empty())
+    {
+        A2AVectorsIo::write(par().output + "_w", w, par().multiFile, vm().getTrajectory());
+        A2AVectorsIo::write(par().output + "_v", v, par().multiFile, vm().getTrajectory());
+    }
 }
 
 END_MODULE_NAMESPACE
diff --git a/Hadrons/Modules/MSource/Momentum.cc b/Hadrons/Modules/MSource/Momentum.cc
new file mode 100644
index 00000000..9bcf65ae
--- /dev/null
+++ b/Hadrons/Modules/MSource/Momentum.cc
@@ -0,0 +1,36 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MSource/Momentum.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+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 <Hadrons/Modules/MSource/Momentum.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MSource;
+
+template class Grid::Hadrons::MSource::TMomentum<FIMPL>;
+
diff --git a/Hadrons/Modules/MSource/Momentum.hpp b/Hadrons/Modules/MSource/Momentum.hpp
new file mode 100644
index 00000000..8336ea5c
--- /dev/null
+++ b/Hadrons/Modules/MSource/Momentum.hpp
@@ -0,0 +1,149 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MSource/Momentum.hpp
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+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 */
+#ifndef Hadrons_Momentum_hpp_
+#define Hadrons_Momentum_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/* 
+Plane Wave source
+-----------------
+src_x = e^i2pi/L * p *position
+*/
+
+/******************************************************************************
+ *                          Plane Wave source                                 *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MSource)
+
+class MomentumPar: Serializable
+{
+public:
+//What is meant by serializable in this context
+GRID_SERIALIZABLE_CLASS_MEMBERS(MomentumPar,
+std::string, mom);
+};
+
+
+template <typename FImpl>
+class TMomentum: public Module<MomentumPar>
+{
+public:
+FERM_TYPE_ALIASES(FImpl,);
+public:
+// constructor
+TMomentum(const std::string name);
+// destructor
+virtual ~TMomentum(void) {};
+// dependency relation
+virtual std::vector<std::string> getInput(void);
+virtual std::vector<std::string> getOutput(void);
+// setup
+virtual void setup(void);
+// execution
+virtual void execute(void);
+};
+
+MODULE_REGISTER_TMP(Momentum, TMomentum<FIMPL>, MSource);
+//MODULE_REGISTER_NS(Momentum, TMomentum, MSource);
+
+/******************************************************************************
+*                       TMomentum template implementation                     *
+******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TMomentum<FImpl>::TMomentum(const std::string name)
+: Module<MomentumPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TMomentum<FImpl>::getInput(void)
+{
+    std::vector<std::string> in;
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TMomentum<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    return out;
+}
+
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TMomentum<FImpl>::setup(void)
+{
+    envCreateLat(PropagatorField, getName());
+}
+
+
+//execution//////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TMomentum<FImpl>::execute(void)
+{
+    LOG(Message) << "Generating planewave momentum source with momentum " << par().mom << std::endl;
+    //what does this env do?
+    PropagatorField &src = envGet(PropagatorField, getName());
+    Lattice<iScalar<vInteger>> t(env().getGrid());
+    LatticeComplex             C(env().getGrid()), coor(env().getGrid());
+    std::vector<Real>          p;
+    std::vector<Real> latt_size(GridDefaultLatt().begin(), GridDefaultLatt().end()); 
+    Complex                    i(0.0,1.0);
+
+    LOG(Message) << " " << std::endl;
+    //get the momentum from parameters
+    p  = strToVec<Real>(par().mom);
+    C = zero;
+    LOG(Message) << "momentum converted from string - " << std::to_string(p[0]) <<std::to_string(p[1]) <<std::to_string(p[2]) <<   std::to_string(p[3]) << std::endl;
+    for(int mu=0;mu<4;mu++){
+    Real TwoPiL =  M_PI * 2.0/ latt_size[mu];
+    LatticeCoordinate(coor,mu);
+    C = C +(TwoPiL * p[mu]) * coor;
+    }
+    C = exp(C*i);
+    LOG(Message) << "exponential of pdotx taken " << std::endl;
+    src = src + C;
+    LOG(Message) << "source created" << std::endl;
+
+}
+
+
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Momentum_hpp_
diff --git a/Hadrons/TimerArray.cc b/Hadrons/TimerArray.cc
index c99ed41f..2b85bc7e 100644
--- a/Hadrons/TimerArray.cc
+++ b/Hadrons/TimerArray.cc
@@ -1,3 +1,30 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/TimerArray.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+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 <Hadrons/TimerArray.hpp>
 
 using namespace Grid;
diff --git a/Hadrons/TimerArray.hpp b/Hadrons/TimerArray.hpp
index 477f11cd..77cc2b8c 100644
--- a/Hadrons/TimerArray.hpp
+++ b/Hadrons/TimerArray.hpp
@@ -1,3 +1,30 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/TimerArray.hpp
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+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 */
 #ifndef Hadrons_TimerArray_hpp_
 #define Hadrons_TimerArray_hpp_
 
diff --git a/Hadrons/Utilities/EigenPackCast.cc b/Hadrons/Utilities/EigenPackCast.cc
index 7247c92d..8f2755e3 100644
--- a/Hadrons/Utilities/EigenPackCast.cc
+++ b/Hadrons/Utilities/EigenPackCast.cc
@@ -1,3 +1,30 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Utilities/EigenPackCast.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+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 <Hadrons/EigenPack.hpp>
 #include <Hadrons/Environment.hpp>
 
diff --git a/Hadrons/Utilities/HadronsXmlRun.cc b/Hadrons/Utilities/HadronsXmlRun.cc
index c6ea7dc5..9b8fc901 100644
--- a/Hadrons/Utilities/HadronsXmlRun.cc
+++ b/Hadrons/Utilities/HadronsXmlRun.cc
@@ -2,7 +2,7 @@
 
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: Hadrons/HadronsXmlRun.cc
+Source file: Hadrons/Utilities/HadronsXmlRun.cc
 
 Copyright (C) 2015-2018
 
diff --git a/Hadrons/modules.inc b/Hadrons/modules.inc
index 74138238..c716b470 100644
--- a/Hadrons/modules.inc
+++ b/Hadrons/modules.inc
@@ -11,6 +11,7 @@ modules_cc =\
   Modules/MContraction/Gamma3pt.cc \
   Modules/MFermion/FreeProp.cc \
   Modules/MFermion/GaugeProp.cc \
+  Modules/MSource/Momentum.cc \
   Modules/MSource/Point.cc \
   Modules/MSource/Wall.cc \
   Modules/MSource/SeqConserved.cc \
@@ -28,6 +29,7 @@ modules_cc =\
   Modules/MGauge/StochEm.cc \
   Modules/MGauge/Random.cc \
   Modules/MGauge/FundtoHirep.cc \
+  Modules/MGauge/GaugeFix.cc \
   Modules/MNoise/TimeDilutedSpinColorDiagonal.cc \
   Modules/MUtilities/RandomVectors.cc \
   Modules/MUtilities/TestSeqGamma.cc \
@@ -38,6 +40,9 @@ modules_cc =\
   Modules/MScalar/VPCounterTerms.cc \
   Modules/MScalar/ChargedProp.cc \
   Modules/MScalar/ScalarVP.cc \
+  Modules/MNPR/Amputate.cc \
+  Modules/MNPR/Bilinear.cc \
+  Modules/MNPR/FourQuark.cc \
   Modules/MAction/Wilson.cc \
   Modules/MAction/MobiusDWF.cc \
   Modules/MAction/ZMobiusDWF.cc \
@@ -59,7 +64,8 @@ modules_cc =\
   Modules/MIO/LoadBinary.cc \
   Modules/MIO/LoadNersc.cc \
   Modules/MIO/LoadCoarseEigenPack.cc \
-  Modules/MIO/LoadCosmHol.cc
+  Modules/MIO/LoadCosmHol.cc \
+  Modules/MIO/LoadA2AVectors.cc
 
 modules_hpp =\
   Modules/MContraction/Baryon.hpp \
@@ -80,6 +86,7 @@ modules_hpp =\
   Modules/MSource/Wall.hpp \
   Modules/MSource/Z2.hpp \
   Modules/MSource/SeqConserved.hpp \
+  Modules/MSource/Momentum.hpp \
   Modules/MSink/Smear.hpp \
   Modules/MSink/Point.hpp \
   Modules/MSolver/MixedPrecisionRBPrecCG.hpp \
@@ -91,6 +98,7 @@ modules_hpp =\
   Modules/MGauge/StoutSmearing.hpp \
   Modules/MGauge/Unit.hpp \
   Modules/MGauge/Random.hpp \
+  Modules/MGauge/GaugeFix.hpp \
   Modules/MGauge/FundtoHirep.hpp \
   Modules/MGauge/StochEm.hpp \
   Modules/MNoise/TimeDilutedSpinColorDiagonal.hpp \
@@ -104,6 +112,9 @@ modules_hpp =\
   Modules/MScalar/ScalarVP.hpp \
   Modules/MScalar/Scalar.hpp \
   Modules/MScalar/ChargedProp.hpp \
+  Modules/MNPR/Bilinear.hpp \
+  Modules/MNPR/Amputate.hpp \
+  Modules/MNPR/FourQuark.hpp \
   Modules/MAction/DWF.hpp \
   Modules/MAction/MobiusDWF.hpp \
   Modules/MAction/Wilson.hpp \
@@ -124,6 +135,7 @@ modules_hpp =\
   Modules/MScalarSUN/TrKinetic.hpp \
   Modules/MIO/LoadEigenPack.hpp \
   Modules/MIO/LoadNersc.hpp \
+  Modules/MIO/LoadA2AVectors.hpp \
   Modules/MIO/LoadCosmHol.hpp \
   Modules/MIO/LoadCoarseEigenPack.hpp \
   Modules/MIO/LoadBinary.hpp
diff --git a/benchmarks/Benchmark_IO.cc b/benchmarks/Benchmark_IO.cc
index 479ae037..58e0340d 100644
--- a/benchmarks/Benchmark_IO.cc
+++ b/benchmarks/Benchmark_IO.cc
@@ -1,108 +1,48 @@
-#include <Grid/Grid.h>
-#ifdef HAVE_LIME
 
-using namespace std;
-using namespace Grid;
-using namespace Grid::QCD;
+#include "Benchmark_IO.hpp"
 
-#define MSG cout << GridLogMessage
-#define SEP \
-"============================================================================="
 #ifndef BENCH_IO_LMAX
 #define BENCH_IO_LMAX 40
 #endif
 
-typedef function<void(const string, LatticeFermion &)> WriterFn;
-typedef function<void(LatticeFermion &, const string)> ReaderFn;
+using namespace Grid;
+using namespace QCD;
 
-string filestem(const int l)
+std::string filestem(const int l)
 {
-  return "iobench_l" + to_string(l);
-}
-
-void limeWrite(const string filestem, LatticeFermion &vec)
-{
-  emptyUserRecord record;
-  ScidacWriter    binWriter(vec._grid->IsBoss());
-
-  binWriter.open(filestem + ".bin");
-  binWriter.writeScidacFieldRecord(vec, record);
-  binWriter.close();
-}
-
-void limeRead(LatticeFermion &vec, const string filestem)
-{
-  emptyUserRecord record;
-  ScidacReader    binReader;
-
-  binReader.open(filestem + ".bin");
-  binReader.readScidacFieldRecord(vec, record);
-  binReader.close();
-}
-
-void writeBenchmark(const int l, const WriterFn &write)
-{
-  auto                      mpi  = GridDefaultMpi();
-  auto                      simd = GridDefaultSimd(Nd, vComplex::Nsimd());
-  vector<int>               latt = {l*mpi[0], l*mpi[1], l*mpi[2], l*mpi[3]};
-  unique_ptr<GridCartesian> gPt(SpaceTimeGrid::makeFourDimGrid(latt, simd, mpi));
-  GridCartesian             *g = gPt.get();
-  GridParallelRNG           rng(g);
-  LatticeFermion            vec(g);
-  emptyUserRecord           record;
-  ScidacWriter              binWriter(g->IsBoss());
-
-  cout << "-- Local volume " << l << "^4" << endl;
-  random(rng, vec);
-  write(filestem(l), vec);
-}
-
-void readBenchmark(const int l, const ReaderFn &read)
-{
-  auto                      mpi  = GridDefaultMpi();
-  auto                      simd = GridDefaultSimd(Nd, vComplex::Nsimd());
-  vector<int>               latt = {l*mpi[0], l*mpi[1], l*mpi[2], l*mpi[3]};
-  unique_ptr<GridCartesian> gPt(SpaceTimeGrid::makeFourDimGrid(latt, simd, mpi));
-  GridCartesian             *g = gPt.get();
-  LatticeFermion            vec(g);
-  emptyUserRecord           record;
-  ScidacReader              binReader;
-
-  cout << "-- Local volume " << l << "^4" << endl;
-  read(vec, filestem(l));
+  return "iobench_l" + std::to_string(l);
 }
 
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
 
-  auto simd = GridDefaultSimd(Nd,vComplex::Nsimd());
-  auto mpi  = GridDefaultMpi();
-
   int64_t threads = GridThread::GetThreads();
-  MSG << "Grid is setup to use " << threads << " threads" << endl;
-  MSG << SEP << endl;
-  MSG << "Benchmark Lime write" << endl;
-  MSG << SEP << endl;
+  MSG << "Grid is setup to use " << threads << " threads" << std::endl;
+  MSG << SEP << std::endl;
+  MSG << "Benchmark Lime write" << std::endl;
+  MSG << SEP << std::endl;
   for (int l = 4; l <= BENCH_IO_LMAX; l += 2)
   {
-    writeBenchmark(l, limeWrite);
+    auto             mpi  = GridDefaultMpi();
+    std::vector<int> latt = {l*mpi[0], l*mpi[1], l*mpi[2], l*mpi[3]};
+
+    std::cout << "-- Local volume " << l << "^4" << std::endl;
+    writeBenchmark<LatticeFermion>(latt, filestem(l), limeWrite<LatticeFermion>);
   }
 
-  MSG << "Benchmark Lime read" << endl;
-  MSG << SEP << endl;
+  MSG << "Benchmark Lime read" << std::endl;
+  MSG << SEP << std::endl;
   for (int l = 4; l <= BENCH_IO_LMAX; l += 2)
   {
-    readBenchmark(l, limeRead);
+    auto             mpi  = GridDefaultMpi();
+    std::vector<int> latt = {l*mpi[0], l*mpi[1], l*mpi[2], l*mpi[3]};
+
+    std::cout << "-- Local volume " << l << "^4" << std::endl;
+    readBenchmark<LatticeFermion>(latt, filestem(l), limeRead<LatticeFermion>);
   }
 
   Grid_finalize();
 
   return EXIT_SUCCESS;
 }
-#else
-int main (int argc, char ** argv)
-{
-  return EXIT_SUCCESS;
-}
-#endif
diff --git a/benchmarks/Benchmark_IO.hpp b/benchmarks/Benchmark_IO.hpp
new file mode 100644
index 00000000..9565c928
--- /dev/null
+++ b/benchmarks/Benchmark_IO.hpp
@@ -0,0 +1,107 @@
+#ifndef Benchmark_IO_hpp_
+#define Benchmark_IO_hpp_
+
+#include <Grid/Grid.h>
+
+#define MSG std::cout << GridLogMessage
+#define SEP \
+"============================================================================="
+
+namespace Grid {
+
+template <typename Field>
+using WriterFn = std::function<void(const std::string, Field &)> ;
+template <typename Field>
+using ReaderFn = std::function<void(Field &, const std::string)>;
+
+template <typename Field>
+void limeWrite(const std::string filestem, Field &vec)
+{
+  emptyUserRecord   record;
+  QCD::ScidacWriter binWriter(vec._grid->IsBoss());
+
+  binWriter.open(filestem + ".bin");
+  binWriter.writeScidacFieldRecord(vec, record);
+  binWriter.close();
+}
+
+template <typename Field>
+void limeRead(Field &vec, const std::string filestem)
+{
+  emptyUserRecord   record;
+  QCD::ScidacReader binReader;
+
+  binReader.open(filestem + ".bin");
+  binReader.readScidacFieldRecord(vec, record);
+  binReader.close();
+}
+
+inline void makeGrid(std::shared_ptr<GridBase> &gPt, 
+                     const std::shared_ptr<GridCartesian> &gBasePt,
+                     const unsigned int Ls = 1, const bool rb = false)
+{
+  if (rb)
+  {
+    if (Ls > 1)
+    {
+      gPt.reset(QCD::SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, gBasePt.get()));
+    }
+    else
+    {
+      gPt.reset(QCD::SpaceTimeGrid::makeFourDimRedBlackGrid(gBasePt.get()));
+    }
+  }
+  else
+  {
+    if (Ls > 1)
+    {
+        gPt.reset(QCD::SpaceTimeGrid::makeFiveDimGrid(Ls, gBasePt.get()));
+    }
+    else
+    {
+        gPt = gBasePt;
+    }
+  }
+}
+
+template <typename Field>
+void writeBenchmark(const std::vector<int> &latt, const std::string filename,
+                    const WriterFn<Field> &write, 
+                    const unsigned int Ls = 1, const bool rb = false)
+{
+  auto                           mpi  = GridDefaultMpi();
+  auto                           simd = GridDefaultSimd(latt.size(), Field::vector_type::Nsimd());
+  std::shared_ptr<GridCartesian> gBasePt(QCD::SpaceTimeGrid::makeFourDimGrid(latt, simd, mpi));
+  std::shared_ptr<GridBase>      gPt;
+
+  makeGrid(gPt, gBasePt, Ls, rb);
+
+  GridBase                       *g = gPt.get();
+  GridParallelRNG                rng(g);
+  Field                          vec(g);
+
+  random(rng, vec);
+  write(filename, vec);
+}
+
+template <typename Field>
+void readBenchmark(const std::vector<int> &latt, const std::string filename,
+                   const ReaderFn<Field> &read, 
+                   const unsigned int Ls = 1, const bool rb = false)
+{
+  auto                           mpi  = GridDefaultMpi();
+  auto                           simd = GridDefaultSimd(latt.size(), Field::vector_type::Nsimd());
+  std::shared_ptr<GridCartesian> gBasePt(QCD::SpaceTimeGrid::makeFourDimGrid(latt, simd, mpi));
+  std::shared_ptr<GridBase>      gPt;
+
+  makeGrid(gPt, gBasePt, Ls, rb);
+
+  GridBase                       *g = gPt.get();
+  Field                          vec(g);
+
+  read(vec, filename);
+}
+
+}
+
+#endif // Benchmark_IO_hpp_
diff --git a/benchmarks/Benchmark_IO_vs_dir.cc b/benchmarks/Benchmark_IO_vs_dir.cc
new file mode 100644
index 00000000..15f6194f
--- /dev/null
+++ b/benchmarks/Benchmark_IO_vs_dir.cc
@@ -0,0 +1,79 @@
+#include "Benchmark_IO.hpp"
+
+#define MSG std::cout << GridLogMessage
+#define SEP \
+"============================================================================="
+
+using namespace Grid;
+using namespace QCD;
+
+int main (int argc, char ** argv)
+{
+  std::vector<std::string> dir;
+  unsigned int             Ls;
+  bool                     rb;
+  if (argc < 4)
+  {
+    std::cerr << "usage: " << argv[0] << " <Ls> <RB {0|1}> <dir1> [<dir2> ... <dirn>] [Grid options]";
+    std::cerr << std::endl;
+  }
+  Ls = std::stoi(argv[1]);
+  rb = (std::string(argv[2]) == "1");
+  for (unsigned int i = 3; i < argc; ++i)
+  {
+    std::string a = argv[i];
+
+    if (a[0] != '-')
+    {
+      dir.push_back(std::string(argv[i]));
+    }
+    else
+    {
+      break;
+    }
+  }
+  Grid_init(&argc,&argv);
+
+
+  int64_t threads = GridThread::GetThreads();
+  MSG << "Grid is setup to use " << threads << " threads" << std::endl;
+  MSG << SEP << std::endl;
+  MSG << "Benchmark double precision Lime write" << std::endl;
+  MSG << SEP << std::endl;
+  for (auto &d: dir)
+  {
+    MSG << "-- Directory " << d << std::endl;
+    writeBenchmark<LatticeFermion>(GridDefaultLatt(), d + "/ioBench", limeWrite<LatticeFermion>, Ls, rb);
+  }
+
+  MSG << SEP << std::endl;
+  MSG << "Benchmark double precision Lime read" << std::endl;
+  MSG << SEP << std::endl;
+  for (auto &d: dir)
+  {
+    MSG << "-- Directory " << d << std::endl;
+    readBenchmark<LatticeFermion>(GridDefaultLatt(), d + "/ioBench", limeRead<LatticeFermion>, Ls, rb);
+  }
+
+  MSG << SEP << std::endl;
+  MSG << "Benchmark single precision Lime write" << std::endl;
+  MSG << SEP << std::endl;
+  for (auto &d: dir)
+  {
+    MSG << "-- Directory " << d << std::endl;
+    writeBenchmark<LatticeFermionF>(GridDefaultLatt(), d + "/ioBench", limeWrite<LatticeFermionF>, Ls, rb);
+  }
+
+  MSG << SEP << std::endl;
+  MSG << "Benchmark single precision Lime read" << std::endl;
+  MSG << SEP << std::endl;
+  for (auto &d: dir)
+  {
+    MSG << "-- Directory " << d << std::endl;
+    readBenchmark<LatticeFermionF>(GridDefaultLatt(), d + "/ioBench", limeRead<LatticeFermionF>, Ls, rb);
+  }
+
+  Grid_finalize();
+
+  return EXIT_SUCCESS;
+}
diff --git a/tests/debug/Test_cayley_cg.cc b/tests/debug/Test_cayley_cg.cc
index 57fe7717..5150268f 100644
--- a/tests/debug/Test_cayley_cg.cc
+++ b/tests/debug/Test_cayley_cg.cc
@@ -1,5 +1,4 @@
     /*************************************************************************************
-
     Grid physics library, www.github.com/paboyle/Grid 
 
     Source file: ./tests/Test_cayley_cg.cc
@@ -27,6 +26,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
     *************************************************************************************/
     /*  END LEGAL */
 #include <Grid/Grid.h>
+#include <Grid/qcd/action/fermion/Reconstruct5Dprop.h>
 
 using namespace std;
 using namespace Grid;
@@ -46,6 +46,7 @@ struct scal {
 
 template<class What> 
 void  TestCGinversions(What & Ddwf, 
+		       LatticeGaugeField &Umu,
 		       GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
 		       GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
 		       RealD mass, RealD M5,
@@ -75,6 +76,25 @@ void  TestCGprec(What & Ddwf,
 		 GridParallelRNG *RNG4,
 		 GridParallelRNG *RNG5);
 
+template<class What> 
+void  TestReconstruct5D(What & Ddwf, 
+			LatticeGaugeField &Umu,
+			GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
+			GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
+			RealD mass, RealD M5,
+			GridParallelRNG *RNG4,
+			GridParallelRNG *RNG5);
+
+template<class What,class WhatF> 
+void  TestReconstruct5DFA(What & Ddwf, 
+			  WhatF & DdwfF, 
+			  LatticeGaugeField &Umu,
+			  GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
+			  GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
+			  RealD mass, RealD M5,
+			  GridParallelRNG *RNG4,
+			  GridParallelRNG *RNG5);
+
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
@@ -83,63 +103,104 @@ int main (int argc, char ** argv)
   std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
 
   const int Ls=8;
-  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), 
+								   GridDefaultSimd(Nd,vComplex::Nsimd()),
+								   GridDefaultMpi());
   GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
   GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
   GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
 
 
+  GridCartesian         * UGridF   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), 
+								    GridDefaultSimd(Nd,vComplexF::Nsimd()),
+								    GridDefaultMpi());
+  GridRedBlackCartesian * UrbGridF = SpaceTimeGrid::makeFourDimRedBlackGrid(UGridF);
+  GridCartesian         * FGridF   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGridF);
+  GridRedBlackCartesian * FrbGridF = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGridF);
+
+
   std::vector<int> seeds4({1,2,3,4});
   std::vector<int> seeds5({5,6,7,8});
   GridParallelRNG          RNG5(FGrid);  RNG5.SeedFixedIntegers(seeds5);
   GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds4);
 
   LatticeGaugeField Umu(UGrid);
+  LatticeGaugeFieldF UmuF(UGridF);
   SU3::HotConfiguration(RNG4,Umu);
+  precisionChange(UmuF,Umu);
   std::vector<LatticeColourMatrix> U(4,UGrid);
 
   RealD mass=0.1;
   RealD M5  =1.8;
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   std::cout<<GridLogMessage <<"DomainWallFermion test"<<std::endl;
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
-  TestCGinversions<DomainWallFermionR>(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  DomainWallFermionF DdwfF(UmuF,*FGridF,*FrbGridF,*UGridF,*UrbGridF,mass,M5);
+  TestCGinversions<DomainWallFermionR>(Ddwf,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestReconstruct5DFA<DomainWallFermionR,DomainWallFermionF>(Ddwf,DdwfF,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
 
   RealD b=1.5;// Scale factor b+c=2, b-c=1
   RealD c=0.5;
   std::vector<ComplexD> gamma(Ls,ComplexD(1.0,0.0));
 
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   std::cout<<GridLogMessage <<"MobiusFermion test"<<std::endl;
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   MobiusFermionR Dmob(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c);
-  TestCGinversions<MobiusFermionR>(Dmob,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  MobiusFermionF DmobF(UmuF,*FGridF,*FrbGridF,*UGridF,*UrbGridF,mass,M5,b,c);
+  TestCGinversions<MobiusFermionR>(Dmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestReconstruct5DFA<MobiusFermionR,MobiusFermionF>(Dmob,DmobF,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
 
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   std::cout<<GridLogMessage <<"ZMobiusFermion test"<<std::endl;
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   ZMobiusFermionR ZDmob(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,gamma,b,c);
-  TestCGinversions<ZMobiusFermionR>(ZDmob,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestCGinversions<ZMobiusFermionR>(ZDmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestReconstruct5D<ZMobiusFermionR>(ZDmob,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
 
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   std::cout<<GridLogMessage <<"MobiusZolotarevFermion test"<<std::endl;
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   MobiusZolotarevFermionR Dzolo(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c,0.1,2.0);
-  TestCGinversions<MobiusZolotarevFermionR>(Dzolo,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestCGinversions<MobiusZolotarevFermionR>(Dzolo,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestReconstruct5D<MobiusZolotarevFermionR>(Dzolo,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
 
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   std::cout<<GridLogMessage <<"ScaledShamirFermion test"<<std::endl;
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   ScaledShamirFermionR Dsham(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,2.0);
-  TestCGinversions<ScaledShamirFermionR>(Dsham,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  ScaledShamirFermionF DshamF(UmuF,*FGridF,*FrbGridF,*UGridF,*UrbGridF,mass,M5,2.0);
+  TestCGinversions<ScaledShamirFermionR>(Dsham,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestReconstruct5DFA<ScaledShamirFermionR,ScaledShamirFermionF>(Dsham,DshamF,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
 
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   std::cout<<GridLogMessage <<"ShamirZolotarevFermion test"<<std::endl;
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   ShamirZolotarevFermionR Dshamz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,2.0);
-  TestCGinversions<ShamirZolotarevFermionR>(Dshamz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestCGinversions<ShamirZolotarevFermionR>(Dshamz,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestReconstruct5D<ShamirZolotarevFermionR>(Dshamz,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
 
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   std::cout<<GridLogMessage <<"OverlapWilsonCayleyTanhFermion test"<<std::endl;
-  OverlapWilsonCayleyTanhFermionR Dov(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0);
-  TestCGinversions<OverlapWilsonCayleyTanhFermionR>(Dov,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
+  OverlapWilsonCayleyTanhFermionR Dov (Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0);
+  OverlapWilsonCayleyTanhFermionF DovF(UmuF,*FGridF,*FrbGridF,*UGridF,*UrbGridF,mass,M5,1.0);
+  TestCGinversions<OverlapWilsonCayleyTanhFermionR>(Dov,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestReconstruct5DFA<OverlapWilsonCayleyTanhFermionR,OverlapWilsonCayleyTanhFermionF>(Dov,DovF,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
 
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   std::cout<<GridLogMessage <<"OverlapWilsonCayleyZolotarevFermion test"<<std::endl;
+  std::cout<<GridLogMessage <<"======================"<<std::endl;
   OverlapWilsonCayleyZolotarevFermionR Dovz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,2.0);
-  TestCGinversions<OverlapWilsonCayleyZolotarevFermionR>(Dovz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestCGinversions<OverlapWilsonCayleyZolotarevFermionR>(Dovz,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
+  TestReconstruct5D<OverlapWilsonCayleyZolotarevFermionR>(Dovz,Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
 
   Grid_finalize();
 }
 template<class What> 
 void  TestCGinversions(What & Ddwf, 
+		       LatticeGaugeField &Umu,
 		       GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
 		       GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
 		       RealD mass, RealD M5,
@@ -154,6 +215,7 @@ void  TestCGinversions(What & Ddwf,
   TestCGschur<What>(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5);
 }
 
+
 template<class What> 
 void  TestCGunprec(What & Ddwf, 
 		   GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
@@ -189,6 +251,147 @@ void  TestCGprec(What & Ddwf,
   CG(HermOpEO,src_o,result_o);
 }
 
+template<class What> 
+void  TestReconstruct5D(What & Ddwf, 
+			LatticeGaugeField & Umu,
+			GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
+			GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
+			RealD mass, RealD M5,
+			GridParallelRNG *RNG4,
+			GridParallelRNG *RNG5)
+{
+  LatticeFermion src4   (UGrid); random(*RNG4,src4);
+  LatticeFermion res4   (UGrid); res4 = zero;
+
+  LatticeFermion src   (FGrid);
+  LatticeFermion src_NE(FGrid);
+  LatticeFermion result(FGrid);
+  LatticeFermion result_rec(FGrid);
+  LatticeFermion result_madwf(FGrid);
+
+  MdagMLinearOperator<What,LatticeFermion> HermOp(Ddwf);
+  double Resid = 1.0e-12;
+  double Residi = 1.0e-6;
+  ConjugateGradient<LatticeFermion> CG(Resid,10000);
+  ConjugateGradient<LatticeFermion> CGi(Residi,10000);
+
+  Ddwf.ImportPhysicalFermionSource(src4,src);
+  Ddwf.Mdag(src,src_NE);
+  CG(HermOp,src_NE,result);
+
+  Ddwf.ExportPhysicalFermionSolution(result, res4);
+
+  Ddwf.M(result,src_NE);
+  src_NE = src_NE - src;
+  std::cout <<GridLogMessage<< " True residual is " << norm2(src_NE)<<std::endl;
+
+  std::cout <<GridLogMessage<< " Reconstructing " <<std::endl;
+
+  ////////////////////////////
+  // RBprec PV inverse
+  ////////////////////////////
+  typedef LatticeFermion Field;
+  typedef SchurRedBlackDiagTwoSolve<Field> SchurSolverType; 
+  typedef SchurRedBlackDiagTwoSolve<Field> SchurSolverTypei; 
+  typedef PauliVillarsSolverRBprec<Field,SchurSolverType> PVinverter;
+  SchurSolverType SchurSolver(CG);
+  PVinverter      PVinverse(SchurSolver);
+
+  Reconstruct5DfromPhysical<LatticeFermion,PVinverter> reconstructor(PVinverse);
+
+  reconstructor(Ddwf,res4,src4,result_rec);
+
+  std::cout <<GridLogMessage << "Result     "<<norm2(result)<<std::endl;
+  std::cout <<GridLogMessage << "Result_rec "<<norm2(result_rec)<<std::endl;
+
+  result_rec = result_rec - result;
+  std::cout <<GridLogMessage << "Difference "<<norm2(result_rec)<<std::endl;
+
+  //////////////////////////////
+  // Now try MADWF
+  //////////////////////////////
+  SchurSolverTypei SchurSolveri(CGi);
+  ZeroGuesser<LatticeFermion> Guess;
+  MADWF<What,What,PVinverter,SchurSolverTypei,ZeroGuesser<LatticeFermion> > 
+    madwf(Ddwf,Ddwf,PVinverse,SchurSolveri,Guess,Resid,10);
+  
+  madwf(src4,result_madwf);
+  result_madwf = result_madwf - result;
+  std::cout <<GridLogMessage << "Difference "<<norm2(result_madwf)<<std::endl;
+
+
+}
+template<class What,class WhatF> 
+void  TestReconstruct5DFA(What & Ddwf, 
+			  WhatF & DdwfF, 
+			  LatticeGaugeField & Umu,
+			  GridCartesian         * FGrid,	       GridRedBlackCartesian * FrbGrid,
+			  GridCartesian         * UGrid,	       GridRedBlackCartesian * UrbGrid,
+			  RealD mass, RealD M5,
+			  GridParallelRNG *RNG4,
+			  GridParallelRNG *RNG5)
+{
+  LatticeFermion src4   (UGrid); random(*RNG4,src4);
+  LatticeFermion res4   (UGrid); res4 = zero;
+
+  LatticeFermion src   (FGrid);
+  LatticeFermion src_NE(FGrid);
+  LatticeFermion result(FGrid);
+  LatticeFermion result_rec(FGrid);
+  LatticeFermion result_madwf(FGrid);
+
+  MdagMLinearOperator<What,LatticeFermion> HermOp(Ddwf);
+  double Resid = 1.0e-12;
+  double Residi = 1.0e-5;
+  ConjugateGradient<LatticeFermion> CG(Resid,10000);
+  ConjugateGradient<LatticeFermionF> CGi(Residi,10000);
+
+  Ddwf.ImportPhysicalFermionSource(src4,src);
+  Ddwf.Mdag(src,src_NE);
+  CG(HermOp,src_NE,result);
+
+  Ddwf.ExportPhysicalFermionSolution(result, res4);
+
+  Ddwf.M(result,src_NE);
+  src_NE = src_NE - src;
+  std::cout <<GridLogMessage<< " True residual is " << norm2(src_NE)<<std::endl;
+
+  std::cout <<GridLogMessage<< " Reconstructing " <<std::endl;
+
+  ////////////////////////////
+  // Fourier accel PV inverse
+  ////////////////////////////
+  typedef LatticeFermion Field;
+  typedef LatticeFermionF FieldF;
+  typedef SchurRedBlackDiagTwoSolve<FieldF> SchurSolverTypei; 
+  typedef PauliVillarsSolverFourierAccel<LatticeFermion,LatticeGaugeField> PVinverter;
+  PVinverter PVinverse(Umu,CG);
+
+  Reconstruct5DfromPhysical<LatticeFermion,PVinverter> reconstructor(PVinverse);
+
+  reconstructor(Ddwf,res4,src4,result_rec);
+
+  std::cout <<GridLogMessage << "Result     "<<norm2(result)<<std::endl;
+  std::cout <<GridLogMessage << "Result_rec "<<norm2(result_rec)<<std::endl;
+
+  result_rec = result_rec - result;
+  std::cout <<GridLogMessage << "Difference "<<norm2(result_rec)<<std::endl;
+
+  //////////////////////////////
+  // Now try MADWF
+  //////////////////////////////
+  SchurSolverTypei SchurSolver(CGi);
+  ZeroGuesser<LatticeFermionF> Guess;
+  MADWF<What,WhatF,PVinverter,SchurSolverTypei,ZeroGuesser<LatticeFermionF> > 
+    madwf(Ddwf,DdwfF,PVinverse,SchurSolver,Guess,Resid,10);
+  
+  madwf(src4,result_madwf);
+  result_madwf = result_madwf - result;
+  std::cout <<GridLogMessage << "Difference "<<norm2(result_madwf)<<std::endl;
+
+}
+
+
 
 template<class What> 
 void  TestCGschur(What & Ddwf,