diff --git a/benchmarks/Benchmark_meson_field.cc b/benchmarks/Benchmark_meson_field.cc
new file mode 100644
index 00000000..442768d5
--- /dev/null
+++ b/benchmarks/Benchmark_meson_field.cc
@@ -0,0 +1,252 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./benchmarks/Benchmark_wilson.cc
+
+    Copyright (C) 2018
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: paboyle <paboyle@ph.ed.ac.uk>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+    *************************************************************************************/
+    /*  END LEGAL */
+#include <Grid/Grid.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+
+#include "Grid/util/Profiling.h"
+
+template<class vobj>
+void sliceInnerProductMesonField(std::vector< std::vector<ComplexD> > &mat, 
+				 const std::vector<Lattice<vobj> > &lhs,
+				 const std::vector<Lattice<vobj> > &rhs,
+				 int orthogdim) 
+{
+  typedef typename vobj::scalar_object sobj;
+  typedef typename vobj::scalar_type scalar_type;
+  typedef typename vobj::vector_type vector_type;
+  
+  int Lblock = lhs.size();
+  int Rblock = rhs.size();
+
+  GridBase *grid = lhs[0]._grid;
+  
+  const int    Nd = grid->_ndimension;
+  const int Nsimd = grid->Nsimd();
+  int Nt     = grid->GlobalDimensions()[orthogdim];
+
+  assert(mat.size()==Lblock*Rblock);
+  for(int t=0;t<mat.size();t++){
+    assert(mat[t].size()==Nt);
+  }
+
+  int fd=grid->_fdimensions[orthogdim];
+  int ld=grid->_ldimensions[orthogdim];
+  int rd=grid->_rdimensions[orthogdim];
+
+  // will locally sum vectors first
+  // sum across these down to scalars
+  // splitting the SIMD
+  std::vector<vector_type,alignedAllocator<vector_type> > lvSum(rd*Lblock*Rblock);
+  for (int r = 0; r < rd * Lblock * Rblock; r++){
+    lvSum[r] = zero;
+  }
+  std::vector<scalar_type > lsSum(ld*Lblock*Rblock,scalar_type(0.0));             
+
+  int e1=    grid->_slice_nblock[orthogdim];
+  int e2=    grid->_slice_block [orthogdim];
+  int stride=grid->_slice_stride[orthogdim];
+  
+  std::cout << GridLogMessage << " Entering first parallel loop "<<std::endl;
+  // Parallelise over t-direction doesn't expose as much parallelism as needed for KNL
+  parallel_for(int r=0;r<rd;r++){
+
+    int so=r*grid->_ostride[orthogdim]; // base offset for start of plane 
+
+    for(int n=0;n<e1;n++){
+      for(int b=0;b<e2;b++){
+	int ss= so+n*stride+b;
+	for(int i=0;i<Lblock;i++){
+	  auto left = conjugate(lhs[i]._odata[ss]);
+	  for(int j=0;j<Rblock;j++){
+	    int idx = i+Lblock*j+Lblock*Rblock*r;
+	    auto right = rhs[j]._odata[ss];
+#if 1	    
+	    vector_type vv = left()(0)(0) * right()(0)(0)
+	      +              left()(0)(1) * right()(0)(1)
+	      +              left()(0)(2) * right()(0)(2)
+              +              left()(1)(0) * right()(1)(0)
+	      +              left()(1)(1) * right()(1)(1)
+	      +              left()(1)(2) * right()(1)(2)
+              +              left()(2)(0) * right()(2)(0)
+	      +              left()(2)(1) * right()(2)(1)
+	      +              left()(2)(2) * right()(2)(2)
+              +              left()(3)(0) * right()(3)(0)
+	      +              left()(3)(1) * right()(3)(1)
+	      +              left()(3)(2) * right()(3)(2);
+#else
+	    vector_type vv = TensorRemove(innerProduct(left,right));
+#endif
+	    lvSum[idx]=lvSum[idx]+vv;
+	  }
+	}
+      }
+    }
+  }
+
+  std::cout << GridLogMessage << " Entering second parallel loop "<<std::endl;
+  // Sum across simd lanes in the plane, breaking out orthog dir.
+  parallel_for(int rt=0;rt<rd;rt++){
+
+    std::vector<int> icoor(Nd);
+
+    for(int i=0;i<Lblock;i++){
+    for(int j=0;j<Rblock;j++){
+
+      iScalar<vector_type> temp; 
+      std::vector<iScalar<scalar_type> > extracted(Nsimd);               
+
+      temp._internal = lvSum[i+Lblock*j+Lblock*Rblock*rt];
+
+      extract(temp,extracted);
+
+      for(int idx=0;idx<Nsimd;idx++){
+
+	grid->iCoorFromIindex(icoor,idx);
+
+	int ldx =rt+icoor[orthogdim]*rd;
+      
+	int ij_dx = i+Lblock*j+Lblock*Rblock*ldx;
+	lsSum[ij_dx]=lsSum[ij_dx]+extracted[idx]._internal;
+
+      }
+    }}
+  }
+
+  std::cout << GridLogMessage << " Entering non parallel loop "<<std::endl;
+  for(int t=0;t<fd;t++)
+  {
+    int pt = t / ld; // processor plane
+    int lt = t % ld;
+    for(int i=0;i<Lblock;i++){
+    for(int j=0;j<Rblock;j++){
+      if (pt == grid->_processor_coor[orthogdim]){
+        int ij_dx = i + Lblock * j + Lblock * Rblock * lt;
+        mat[i+j*Lblock][t] = lsSum[ij_dx];
+      }
+      else{
+        mat[i+j*Lblock][t] = scalar_type(0.0);
+      }
+    }}
+  }
+  std::cout << GridLogMessage << " Done "<<std::endl;
+  // defer sum over nodes.
+  return;
+}
+
+/*
+template void sliceInnerProductMesonField<SpinColourVector>(std::vector< std::vector<ComplexD> > &mat, 
+						   const std::vector<Lattice<SpinColourVector> > &lhs,
+						   const std::vector<Lattice<SpinColourVector> > &rhs,
+						   int orthogdim) ;
+*/
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  std::vector<int> latt_size   = GridDefaultLatt();
+  std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
+  std::vector<int> mpi_layout  = GridDefaultMpi();
+  GridCartesian               Grid(latt_size,simd_layout,mpi_layout);
+
+  int nt = latt_size[Tp];
+  uint64_t vol = 1;
+  for(int d=0;d<Nd;d++){
+    vol = vol*latt_size[d];
+  }
+  
+  std::vector<int> seeds({1,2,3,4});
+  GridParallelRNG          pRNG(&Grid);
+  pRNG.SeedFixedIntegers(seeds);
+
+
+  const int Nm = 32; // number of all modes (high + low)
+
+  std::vector<LatticeFermion> v(Nm,&Grid);
+  std::vector<LatticeFermion> w(Nm,&Grid);
+
+  for(int i=0;i<Nm;i++) { 
+    random(pRNG,v[i]);
+    random(pRNG,w[i]);
+  }
+
+  double flops = vol * (11.0 * 8.0 + 6.0) * Nm*Nm;
+  double byte  = vol * (12.0 * sizeof(Complex) ) * Nm*Nm;
+
+  std::vector<ComplexD> ip(nt);
+  std::vector<std::vector<ComplexD> > MesonFields   (Nm*Nm);
+  std::vector<std::vector<ComplexD> > MesonFieldsRef(Nm*Nm);
+
+  for(int i=0;i<Nm;i++) { 
+  for(int j=0;j<Nm;j++) { 
+    MesonFields   [i+j*Nm].resize(nt);
+    MesonFieldsRef[i+j*Nm].resize(nt);
+  }}
+
+  GridLogMessage.TimingMode(1);
+
+  std::cout<<GridLogMessage << "Running loop with sliceInnerProductVector"<<std::endl;
+  double t0 = usecond();
+  for(int i=0;i<Nm;i++) { 
+  for(int j=0;j<Nm;j++) { 
+    sliceInnerProductVector(ip, w[i],v[j],Tp);
+    for(int t=0;t<nt;t++){
+      MesonFieldsRef[i+j*Nm][t] = ip[t];
+    }
+  }}
+  double t1 = usecond();
+  std::cout<<GridLogMessage << "Done "<< flops/(t1-t0) <<" mflops " <<std::endl;
+  std::cout<<GridLogMessage << "Done "<< byte /(t1-t0) <<" MB/s " <<std::endl;
+
+  std::cout<<GridLogMessage << "Running loop with new code for Nt="<<nt<<std::endl;
+  double t2 = usecond();
+  sliceInnerProductMesonField(MesonFields,w,v,Tp);
+  double t3 = usecond();
+  std::cout<<GridLogMessage << "Done "<< flops/(t3-t2) <<" mflops " <<std::endl;
+  std::cout<<GridLogMessage << "Done "<< byte /(t3-t2) <<" MB/s " <<std::endl;
+
+  RealD err = 0;
+  ComplexD diff;
+
+  for(int i=0;i<Nm;i++) { 
+  for(int j=0;j<Nm;j++) { 
+    for(int t=0;t<nt;t++){
+      diff = MesonFields[i+Nm*j][t] - MesonFieldsRef[i+Nm*j][t];
+      err += real(diff*conj(diff));
+    }
+  }}
+  std::cout<<GridLogMessage << "Norm error "<< err <<std::endl;
+  
+  Grid_finalize();
+}
+
diff --git a/extras/Hadrons/AllToAllReduction.hpp b/extras/Hadrons/AllToAllReduction.hpp
new file mode 100644
index 00000000..7b2d2c1a
--- /dev/null
+++ b/extras/Hadrons/AllToAllReduction.hpp
@@ -0,0 +1,146 @@
+#ifndef A2A_Reduction_hpp_
+#define A2A_Reduction_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Environment.hpp>
+#include <Grid/Hadrons/Solver.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+////////////////////////////////////////////
+// A2A Meson Field Inner Product
+////////////////////////////////////////////
+
+template <class FermionField>
+void sliceInnerProductMesonField(std::vector<std::vector<ComplexD>> &mat,
+                                 const std::vector<Lattice<FermionField>> &lhs,
+                                 const std::vector<Lattice<FermionField>> &rhs,
+                                 int orthogdim)
+{
+    typedef typename FermionField::scalar_type scalar_type;
+    typedef typename FermionField::vector_type vector_type;
+
+    int Lblock = lhs.size();
+    int Rblock = rhs.size();
+
+    GridBase *grid = lhs[0]._grid;
+
+    const int Nd = grid->_ndimension;
+    const int Nsimd = grid->Nsimd();
+    int Nt = grid->GlobalDimensions()[orthogdim];
+
+    assert(mat.size() == Lblock * Rblock);
+    for (int t = 0; t < mat.size(); t++)
+    {
+        assert(mat[t].size() == Nt);
+    }
+
+    int fd = grid->_fdimensions[orthogdim];
+    int ld = grid->_ldimensions[orthogdim];
+    int rd = grid->_rdimensions[orthogdim];
+
+    // will locally sum vectors first
+    // sum across these down to scalars
+    // splitting the SIMD
+    std::vector<vector_type, alignedAllocator<vector_type>> lvSum(rd * Lblock * Rblock);
+    for(int r=0;r<rd * Lblock * Rblock;r++)
+    {
+        lvSum[r]=zero;
+    }
+    std::vector<scalar_type> lsSum(ld * Lblock * Rblock, scalar_type(0.0));
+
+    int e1 = grid->_slice_nblock[orthogdim];
+    int e2 = grid->_slice_block[orthogdim];
+    int stride = grid->_slice_stride[orthogdim];
+
+    // std::cout << GridLogMessage << " Entering first parallel loop " << std::endl;
+    // Parallelise over t-direction doesn't expose as much parallelism as needed for KNL
+    parallel_for(int r = 0; r < rd; r++)
+    {
+        int so = r * grid->_ostride[orthogdim]; // base offset for start of plane
+        for (int n = 0; n < e1; n++)
+        {
+            for (int b = 0; b < e2; b++)
+            {
+                int ss = so + n * stride + b;
+                for (int i = 0; i < Lblock; i++)
+                {
+                    auto left = conjugate(lhs[i]._odata[ss]);
+                    for (int j = 0; j < Rblock; j++)
+                    {
+                        int idx = i + Lblock * j + Lblock * Rblock * r;
+                        auto right = rhs[j]._odata[ss];
+                        vector_type vv = left()(0)(0) * right()(0)(0) 
+                                       + left()(0)(1) * right()(0)(1) 
+                                       + left()(0)(2) * right()(0)(2) 
+                                       + left()(1)(0) * right()(1)(0) 
+                                       + left()(1)(1) * right()(1)(1) 
+                                       + left()(1)(2) * right()(1)(2) 
+                                       + left()(2)(0) * right()(2)(0) 
+                                       + left()(2)(1) * right()(2)(1) 
+                                       + left()(2)(2) * right()(2)(2) 
+                                       + left()(3)(0) * right()(3)(0) 
+                                       + left()(3)(1) * right()(3)(1) 
+                                       + left()(3)(2) * right()(3)(2);
+
+                        lvSum[idx] = lvSum[idx] + vv;
+                    }
+                }
+            }
+        }
+    }
+
+    // std::cout << GridLogMessage << " Entering second parallel loop " << std::endl;
+    // Sum across simd lanes in the plane, breaking out orthog dir.
+    parallel_for(int rt = 0; rt < rd; rt++)
+    {
+        std::vector<int> icoor(Nd);
+        for (int i = 0; i < Lblock; i++)
+        {
+            for (int j = 0; j < Rblock; j++)
+            {
+                iScalar<vector_type> temp;
+                std::vector<iScalar<scalar_type>> extracted(Nsimd);
+                temp._internal = lvSum[i + Lblock * j + Lblock * Rblock * rt];
+                extract(temp, extracted);
+                for (int idx = 0; idx < Nsimd; idx++)
+                {
+                    grid->iCoorFromIindex(icoor, idx);
+                    int ldx = rt + icoor[orthogdim] * rd;
+                    int ij_dx = i + Lblock * j + Lblock * Rblock * ldx;
+                    lsSum[ij_dx] = lsSum[ij_dx] + extracted[idx]._internal;
+                }
+            }
+        }
+    }
+
+    // std::cout << GridLogMessage << " Entering non parallel loop " << std::endl;
+    for (int t = 0; t < fd; t++)
+    {
+        int pt = t/ld; // processor plane
+        int lt = t%ld;
+        for (int i = 0; i < Lblock; i++)
+        {
+            for (int j = 0; j < Rblock; j++)
+            {
+                if (pt == grid->_processor_coor[orthogdim])
+                {
+                    int ij_dx = i + Lblock * j + Lblock * Rblock * lt;
+                    mat[i + j * Lblock][t] = lsSum[ij_dx];
+                }
+                else
+                {
+                    mat[i + j * Lblock][t] = scalar_type(0.0);
+                }
+                
+            }
+        }
+    }
+    // std::cout << GridLogMessage << " Done " << std::endl;
+    // defer sum over nodes.
+    return;
+}
+
+END_HADRONS_NAMESPACE
+
+#endif // A2A_Reduction_hpp_
\ No newline at end of file
diff --git a/extras/Hadrons/AllToAllVectors.hpp b/extras/Hadrons/AllToAllVectors.hpp
new file mode 100644
index 00000000..df60b964
--- /dev/null
+++ b/extras/Hadrons/AllToAllVectors.hpp
@@ -0,0 +1,205 @@
+#ifndef A2A_Vectors_hpp_
+#define A2A_Vectors_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Environment.hpp>
+#include <Grid/Hadrons/Solver.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+////////////////////////////////
+// A2A Modes
+////////////////////////////////
+
+template <class Field, class Matrix, class Solver>
+class A2AModesSchurDiagTwo
+{
+  private:
+    const std::vector<Field> *evec;
+    const std::vector<RealD> *eval;
+    Matrix &action;
+    Solver &solver;
+    std::vector<Field> w_high_5d, v_high_5d, w_high_4d, v_high_4d;
+    const int Nl, Nh;
+    const bool return_5d;
+
+  public:
+    A2AModesSchurDiagTwo(const std::vector<Field> *_evec, const std::vector<RealD> *_eval,
+                         Matrix &_action,
+                         Solver &_solver,
+                         std::vector<Field> _w_high_5d, std::vector<Field> _v_high_5d,
+                         std::vector<Field> _w_high_4d, std::vector<Field> _v_high_4d,
+                         const int _Nl, const int _Nh,
+                         const bool _return_5d)
+                        : evec(_evec), eval(_eval),
+                        action(_action),
+                        solver(_solver),
+                        w_high_5d(_w_high_5d), v_high_5d(_v_high_5d),
+                        w_high_4d(_w_high_4d), v_high_4d(_v_high_4d),
+                        Nl(_Nl), Nh(_Nh),
+                        return_5d(_return_5d){};
+
+    void high_modes(Field &source_5d, Field &w_source_5d, Field &source_4d, int i)
+    {
+        int i5d;
+        LOG(Message) << "A2A high modes for i = " << i << std::endl;
+        i5d = 0;
+        if (return_5d) i5d = i;
+        this->high_mode_v(action, solver, source_5d, v_high_5d[i5d], v_high_4d[i]);
+        this->high_mode_w(w_source_5d, source_4d, w_high_5d[i5d], w_high_4d[i]);
+    }
+
+    void return_v(int i, Field &vout_5d, Field &vout_4d)
+    {
+        if (i < Nl)
+        {
+            this->low_mode_v(action, evec->at(i), eval->at(i), vout_5d, vout_4d);
+        }
+        else
+        {
+            vout_4d = v_high_4d[i - Nl];
+            if (!(return_5d)) i = Nl;
+            vout_5d = v_high_5d[i - Nl];
+        }
+    }
+    void return_w(int i, Field &wout_5d, Field &wout_4d)
+    {
+        if (i < Nl)
+        {
+            this->low_mode_w(action, evec->at(i), eval->at(i), wout_5d, wout_4d);
+        }
+        else
+        {
+            wout_4d = w_high_4d[i - Nl];
+            if (!(return_5d)) i = Nl;
+            wout_5d = w_high_5d[i - Nl];
+        }
+    }
+
+    void low_mode_v(Matrix &action, const Field &evec, const RealD &eval, Field &vout_5d, Field &vout_4d)
+    {
+        GridBase *grid = action.RedBlackGrid();
+        Field src_o(grid);
+        Field sol_e(grid);
+        Field sol_o(grid);
+        Field tmp(grid);
+
+        src_o = evec;
+        src_o.checkerboard = Odd;
+        pickCheckerboard(Even, sol_e, vout_5d);
+        pickCheckerboard(Odd, sol_o, vout_5d);
+
+        /////////////////////////////////////////////////////
+        // v_ie = -(1/eval_i) * MeeInv Meo MooInv evec_i
+        /////////////////////////////////////////////////////
+        action.MooeeInv(src_o, tmp);
+        assert(tmp.checkerboard == Odd);
+        action.Meooe(tmp, sol_e);
+        assert(sol_e.checkerboard == Even);
+        action.MooeeInv(sol_e, tmp);
+        assert(tmp.checkerboard == Even);
+        sol_e = (-1.0 / eval) * tmp;
+        assert(sol_e.checkerboard == Even);
+
+        /////////////////////////////////////////////////////
+        // v_io = (1/eval_i) * MooInv evec_i
+        /////////////////////////////////////////////////////
+        action.MooeeInv(src_o, tmp);
+        assert(tmp.checkerboard == Odd);
+        sol_o = (1.0 / eval) * tmp;
+        assert(sol_o.checkerboard == Odd);
+
+        setCheckerboard(vout_5d, sol_e);
+        assert(sol_e.checkerboard == Even);
+        setCheckerboard(vout_5d, sol_o);
+        assert(sol_o.checkerboard == Odd);
+
+        action.ExportPhysicalFermionSolution(vout_5d, vout_4d);
+    }
+
+    void low_mode_w(Matrix &action, const Field &evec, const RealD &eval, Field &wout_5d, Field &wout_4d)
+    {
+        GridBase *grid = action.RedBlackGrid();
+        SchurDiagTwoOperator<Matrix, Field> _HermOpEO(action);
+
+        Field src_o(grid);
+        Field sol_e(grid);
+        Field sol_o(grid);
+        Field tmp(grid);
+
+        GridBase *fgrid = action.Grid();
+        Field tmp_wout(fgrid);
+
+        src_o = evec;
+        src_o.checkerboard = Odd;
+        pickCheckerboard(Even, sol_e, tmp_wout);
+        pickCheckerboard(Odd, sol_o, tmp_wout);
+
+        /////////////////////////////////////////////////////
+        // w_ie = - MeeInvDag MoeDag Doo evec_i
+        /////////////////////////////////////////////////////
+        _HermOpEO.Mpc(src_o, tmp);
+        assert(tmp.checkerboard == Odd);
+        action.MeooeDag(tmp, sol_e);
+        assert(sol_e.checkerboard == Even);
+        action.MooeeInvDag(sol_e, tmp);
+        assert(tmp.checkerboard == Even);
+        sol_e = (-1.0) * tmp;
+
+        /////////////////////////////////////////////////////
+        // w_io = Doo evec_i
+        /////////////////////////////////////////////////////
+        _HermOpEO.Mpc(src_o, sol_o);
+        assert(sol_o.checkerboard == Odd);
+
+        setCheckerboard(tmp_wout, sol_e);
+        assert(sol_e.checkerboard == Even);
+        setCheckerboard(tmp_wout, sol_o);
+        assert(sol_o.checkerboard == Odd);
+
+        action.DminusDag(tmp_wout, wout_5d);
+
+        action.ExportPhysicalFermionSource(wout_5d, wout_4d);
+    }
+
+    void high_mode_v(Matrix &action, Solver &solver, const Field &source, Field &vout_5d, Field &vout_4d)
+    {
+        GridBase *fgrid = action.Grid();
+        solver(vout_5d, source); // Note: solver is solver(out, in)
+        action.ExportPhysicalFermionSolution(vout_5d, vout_4d);
+    }
+
+    void high_mode_w(const Field &w_source_5d, const Field &source_4d, Field &wout_5d, Field &wout_4d)
+    {
+        wout_5d = w_source_5d;
+        wout_4d = source_4d;
+    }
+};
+
+// TODO: A2A for coarse eigenvectors
+
+// template <class FineField, class CoarseField, class Matrix, class Solver>
+// class A2ALMSchurDiagTwoCoarse : public A2AModesSchurDiagTwo<FineField, Matrix, Solver>
+// {
+//   private:
+//     const std::vector<FineField> &subspace;
+//     const std::vector<CoarseField> &evec_coarse;
+//     const std::vector<RealD> &eval_coarse;
+//     Matrix &action;
+
+//   public:
+//     A2ALMSchurDiagTwoCoarse(const std::vector<FineField> &_subspace, const std::vector<CoarseField> &_evec_coarse, const std::vector<RealD> &_eval_coarse, Matrix &_action)
+//         : subspace(_subspace), evec_coarse(_evec_coarse), eval_coarse(_eval_coarse), action(_action){};
+
+//     void operator()(int i, FineField &vout, FineField &wout)
+//     {
+//         FineField prom_evec(subspace[0]._grid);
+//         blockPromote(evec_coarse[i], prom_evec, subspace);
+//         this->low_mode_v(action, prom_evec, eval_coarse[i], vout);
+//         this->low_mode_w(action, prom_evec, eval_coarse[i], wout);
+//     }
+// };
+
+END_HADRONS_NAMESPACE
+
+#endif // A2A_Vectors_hpp_
\ No newline at end of file
diff --git a/extras/Hadrons/Makefile.am b/extras/Hadrons/Makefile.am
index eb0bc3ad..abfcf9a8 100644
--- a/extras/Hadrons/Makefile.am
+++ b/extras/Hadrons/Makefile.am
@@ -14,6 +14,8 @@ libHadrons_a_SOURCES = \
 libHadrons_adir = $(pkgincludedir)/Hadrons
 nobase_libHadrons_a_HEADERS = \
 	$(modules_hpp)            \
+	AllToAllVectors.hpp       \
+	AllToAllReduction.hpp       \
 	Application.hpp           \
 	EigenPack.hpp             \
 	Environment.hpp           \
diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp
index 6bf5d3ed..58e9fec3 100644
--- a/extras/Hadrons/Modules.hpp
+++ b/extras/Hadrons/Modules.hpp
@@ -1,57 +1,60 @@
-#include <Grid/Hadrons/Modules/MContraction/Baryon.hpp>
-#include <Grid/Hadrons/Modules/MContraction/Meson.hpp>
-#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
-#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp>
-#include <Grid/Hadrons/Modules/MContraction/DiscLoop.hpp>
-#include <Grid/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp>
-#include <Grid/Hadrons/Modules/MContraction/Gamma3pt.hpp>
-#include <Grid/Hadrons/Modules/MContraction/WardIdentity.hpp>
-#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp>
-#include <Grid/Hadrons/Modules/MFermion/FreeProp.hpp>
-#include <Grid/Hadrons/Modules/MFermion/GaugeProp.hpp>
+#include <Grid/Hadrons/Modules/MSource/SeqConserved.hpp>
 #include <Grid/Hadrons/Modules/MSource/SeqGamma.hpp>
+#include <Grid/Hadrons/Modules/MSource/Z2.hpp>
 #include <Grid/Hadrons/Modules/MSource/Point.hpp>
 #include <Grid/Hadrons/Modules/MSource/Wall.hpp>
-#include <Grid/Hadrons/Modules/MSource/Z2.hpp>
-#include <Grid/Hadrons/Modules/MSource/SeqConserved.hpp>
 #include <Grid/Hadrons/Modules/MSink/Smear.hpp>
 #include <Grid/Hadrons/Modules/MSink/Point.hpp>
-#include <Grid/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp>
-#include <Grid/Hadrons/Modules/MSolver/RBPrecCG.hpp>
-#include <Grid/Hadrons/Modules/MGauge/UnitEm.hpp>
-#include <Grid/Hadrons/Modules/MGauge/StoutSmearing.hpp>
-#include <Grid/Hadrons/Modules/MGauge/Unit.hpp>
-#include <Grid/Hadrons/Modules/MGauge/Random.hpp>
-#include <Grid/Hadrons/Modules/MGauge/FundtoHirep.hpp>
-#include <Grid/Hadrons/Modules/MGauge/StochEm.hpp>
-#include <Grid/Hadrons/Modules/MUtilities/TestSeqGamma.hpp>
-#include <Grid/Hadrons/Modules/MUtilities/TestSeqConserved.hpp>
-#include <Grid/Hadrons/Modules/MLoop/NoiseLoop.hpp>
-#include <Grid/Hadrons/Modules/MScalar/FreeProp.hpp>
-#include <Grid/Hadrons/Modules/MScalar/VPCounterTerms.hpp>
-#include <Grid/Hadrons/Modules/MScalar/ScalarVP.hpp>
-#include <Grid/Hadrons/Modules/MScalar/Scalar.hpp>
-#include <Grid/Hadrons/Modules/MScalar/ChargedProp.hpp>
-#include <Grid/Hadrons/Modules/MAction/DWF.hpp>
-#include <Grid/Hadrons/Modules/MAction/MobiusDWF.hpp>
-#include <Grid/Hadrons/Modules/MAction/Wilson.hpp>
-#include <Grid/Hadrons/Modules/MAction/WilsonClover.hpp>
-#include <Grid/Hadrons/Modules/MAction/ZMobiusDWF.hpp>
-#include <Grid/Hadrons/Modules/MAction/ScaledDWF.hpp>
-#include <Grid/Hadrons/Modules/MScalarSUN/StochFreeField.hpp>
+#include <Grid/Hadrons/Modules/MIO/LoadBinary.hpp>
+#include <Grid/Hadrons/Modules/MIO/LoadEigenPack.hpp>
+#include <Grid/Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp>
+#include <Grid/Hadrons/Modules/MIO/LoadNersc.hpp>
+#include <Grid/Hadrons/Modules/MScalarSUN/Utils.hpp>
+#include <Grid/Hadrons/Modules/MScalarSUN/Grad.hpp>
+#include <Grid/Hadrons/Modules/MScalarSUN/TrPhi.hpp>
 #include <Grid/Hadrons/Modules/MScalarSUN/TwoPointNPR.hpp>
+#include <Grid/Hadrons/Modules/MScalarSUN/TwoPoint.hpp>
+#include <Grid/Hadrons/Modules/MScalarSUN/TransProj.hpp>
+#include <Grid/Hadrons/Modules/MScalarSUN/TrKinetic.hpp>
+#include <Grid/Hadrons/Modules/MScalarSUN/StochFreeField.hpp>
 #include <Grid/Hadrons/Modules/MScalarSUN/ShiftProbe.hpp>
-#include <Grid/Hadrons/Modules/MScalarSUN/Div.hpp>
 #include <Grid/Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp>
+#include <Grid/Hadrons/Modules/MScalarSUN/Div.hpp>
 #include <Grid/Hadrons/Modules/MScalarSUN/TrMag.hpp>
 #include <Grid/Hadrons/Modules/MScalarSUN/EMT.hpp>
-#include <Grid/Hadrons/Modules/MScalarSUN/TwoPoint.hpp>
-#include <Grid/Hadrons/Modules/MScalarSUN/TrPhi.hpp>
-#include <Grid/Hadrons/Modules/MScalarSUN/Utils.hpp>
-#include <Grid/Hadrons/Modules/MScalarSUN/TransProj.hpp>
-#include <Grid/Hadrons/Modules/MScalarSUN/Grad.hpp>
-#include <Grid/Hadrons/Modules/MScalarSUN/TrKinetic.hpp>
-#include <Grid/Hadrons/Modules/MIO/LoadEigenPack.hpp>
-#include <Grid/Hadrons/Modules/MIO/LoadNersc.hpp>
-#include <Grid/Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp>
-#include <Grid/Hadrons/Modules/MIO/LoadBinary.hpp>
+#include <Grid/Hadrons/Modules/MAction/ZMobiusDWF.hpp>
+#include <Grid/Hadrons/Modules/MAction/ScaledDWF.hpp>
+#include <Grid/Hadrons/Modules/MAction/Wilson.hpp>
+#include <Grid/Hadrons/Modules/MAction/WilsonClover.hpp>
+#include <Grid/Hadrons/Modules/MAction/MobiusDWF.hpp>
+#include <Grid/Hadrons/Modules/MAction/DWF.hpp>
+#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
+#include <Grid/Hadrons/Modules/MContraction/DiscLoop.hpp>
+#include <Grid/Hadrons/Modules/MContraction/Meson.hpp>
+#include <Grid/Hadrons/Modules/MContraction/WardIdentity.hpp>
+#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp>
+#include <Grid/Hadrons/Modules/MContraction/Gamma3pt.hpp>
+#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp>
+#include <Grid/Hadrons/Modules/MContraction/MesonFieldGamma.hpp>
+#include <Grid/Hadrons/Modules/MContraction/Baryon.hpp>
+#include <Grid/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp>
+#include <Grid/Hadrons/Modules/MContraction/A2AMeson.hpp>
+#include <Grid/Hadrons/Modules/MScalar/ScalarVP.hpp>
+#include <Grid/Hadrons/Modules/MScalar/Scalar.hpp>
+#include <Grid/Hadrons/Modules/MScalar/FreeProp.hpp>
+#include <Grid/Hadrons/Modules/MScalar/ChargedProp.hpp>
+#include <Grid/Hadrons/Modules/MScalar/VPCounterTerms.hpp>
+#include <Grid/Hadrons/Modules/MUtilities/TestSeqConserved.hpp>
+#include <Grid/Hadrons/Modules/MUtilities/TestSeqGamma.hpp>
+#include <Grid/Hadrons/Modules/MFermion/FreeProp.hpp>
+#include <Grid/Hadrons/Modules/MFermion/GaugeProp.hpp>
+#include <Grid/Hadrons/Modules/MSolver/A2AVectors.hpp>
+#include <Grid/Hadrons/Modules/MSolver/RBPrecCG.hpp>
+#include <Grid/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp>
+#include <Grid/Hadrons/Modules/MLoop/NoiseLoop.hpp>
+#include <Grid/Hadrons/Modules/MGauge/StoutSmearing.hpp>
+#include <Grid/Hadrons/Modules/MGauge/StochEm.hpp>
+#include <Grid/Hadrons/Modules/MGauge/FundtoHirep.hpp>
+#include <Grid/Hadrons/Modules/MGauge/Unit.hpp>
+#include <Grid/Hadrons/Modules/MGauge/Random.hpp>
+#include <Grid/Hadrons/Modules/MGauge/UnitEm.hpp>
diff --git a/extras/Hadrons/Modules/MContraction/A2AMeson.cc b/extras/Hadrons/Modules/MContraction/A2AMeson.cc
new file mode 100644
index 00000000..5e5c122b
--- /dev/null
+++ b/extras/Hadrons/Modules/MContraction/A2AMeson.cc
@@ -0,0 +1,8 @@
+#include <Grid/Hadrons/Modules/MContraction/A2AMeson.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MContraction;
+
+template class Grid::Hadrons::MContraction::TA2AMeson<FIMPL>;
+template class Grid::Hadrons::MContraction::TA2AMeson<ZFIMPL>;
\ No newline at end of file
diff --git a/extras/Hadrons/Modules/MContraction/A2AMeson.hpp b/extras/Hadrons/Modules/MContraction/A2AMeson.hpp
new file mode 100644
index 00000000..a13336ef
--- /dev/null
+++ b/extras/Hadrons/Modules/MContraction/A2AMeson.hpp
@@ -0,0 +1,207 @@
+#ifndef Hadrons_MContraction_A2AMeson_hpp_
+#define Hadrons_MContraction_A2AMeson_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+#include <Grid/Hadrons/AllToAllVectors.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         A2AMeson                                 *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MContraction)
+
+typedef std::pair<Gamma::Algebra, Gamma::Algebra> GammaPair;
+
+class A2AMesonPar : Serializable
+{
+  public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(A2AMesonPar,
+                                    int, Nl,
+                                    int, N,
+                                    std::string, A2A1,
+                                    std::string, A2A2,
+                                    std::string, gammas,
+                                    std::string, output);
+};
+
+template <typename FImpl>
+class TA2AMeson : public Module<A2AMesonPar>
+{
+  public:
+    FERM_TYPE_ALIASES(FImpl, );
+    SOLVER_TYPE_ALIASES(FImpl, );
+
+    typedef A2AModesSchurDiagTwo<typename FImpl::FermionField, FMat, Solver> A2ABase;
+
+    class Result : Serializable
+    {
+      public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
+                                        Gamma::Algebra, gamma_snk,
+                                        Gamma::Algebra, gamma_src,
+                                        std::vector<Complex>, corr);
+    };
+
+  public:
+    // constructor
+    TA2AMeson(const std::string name);
+    // destructor
+    virtual ~TA2AMeson(void){};
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    virtual void parseGammaString(std::vector<GammaPair> &gammaList);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER(A2AMeson, ARG(TA2AMeson<FIMPL>), MContraction);
+MODULE_REGISTER(ZA2AMeson, ARG(TA2AMeson<ZFIMPL>), MContraction);
+
+/******************************************************************************
+*                  TA2AMeson implementation                             *
+******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TA2AMeson<FImpl>::TA2AMeson(const std::string name)
+    : Module<A2AMesonPar>(name)
+{
+}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TA2AMeson<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().A2A1 + "_class", par().A2A2 + "_class"};
+    in.push_back(par().A2A1 + "_w_high_4d");
+    in.push_back(par().A2A2 + "_v_high_4d");
+
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TA2AMeson<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {};
+
+    return out;
+}
+
+template <typename FImpl>
+void TA2AMeson<FImpl>::parseGammaString(std::vector<GammaPair> &gammaList)
+{
+    gammaList.clear();
+    // Parse individual contractions from input string.
+    gammaList = strToVec<GammaPair>(par().gammas);
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TA2AMeson<FImpl>::setup(void)
+{
+    int nt = env().getDim(Tp);
+    int N = par().N;
+
+    int Ls_ = env().getObjectLs(par().A2A1 + "_class");
+
+    envTmp(std::vector<FermionField>, "w1", 1, N, FermionField(env().getGrid(1)));
+    envTmp(std::vector<FermionField>, "v1", 1, N, FermionField(env().getGrid(1)));
+    envTmpLat(FermionField, "tmpv_5d", Ls_);
+    envTmpLat(FermionField, "tmpw_5d", Ls_);
+
+    envTmp(std::vector<ComplexD>, "MF_x", 1, nt);
+    envTmp(std::vector<ComplexD>, "MF_y", 1, nt);
+    envTmp(std::vector<ComplexD>, "tmp", 1, nt);
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TA2AMeson<FImpl>::execute(void)
+{
+    LOG(Message) << "Computing A2A meson contractions" << std::endl;
+
+    Result result;
+    Gamma g5(Gamma::Algebra::Gamma5);
+    std::vector<GammaPair> gammaList;
+    int nt = env().getDim(Tp);
+
+    parseGammaString(gammaList);
+
+    result.gamma_snk = gammaList[0].first;
+    result.gamma_src = gammaList[0].second;
+    result.corr.resize(nt);
+
+    int Nl = par().Nl;
+    int N = par().N;
+    LOG(Message) << "N for A2A cont: " << N << std::endl;
+
+    envGetTmp(std::vector<ComplexD>, MF_x);
+    envGetTmp(std::vector<ComplexD>, MF_y);
+    envGetTmp(std::vector<ComplexD>, tmp);
+
+    for (unsigned int t = 0; t < nt; ++t)
+    {
+        tmp[t] = TensorRemove(MF_x[t] * MF_y[t] * 0.0);
+    }
+
+    Gamma gSnk(gammaList[0].first);
+    Gamma gSrc(gammaList[0].second);
+
+    auto &a2a1_fn = envGet(A2ABase, par().A2A1 + "_class");
+
+    envGetTmp(std::vector<FermionField>, w1);
+    envGetTmp(std::vector<FermionField>, v1);
+    envGetTmp(FermionField, tmpv_5d);
+    envGetTmp(FermionField, tmpw_5d);
+
+    LOG(Message) << "Finding v and w vectors for N =  " << N << std::endl;
+    for (int i = 0; i < N; i++)
+    {
+        a2a1_fn.return_v(i, tmpv_5d, v1[i]);
+        a2a1_fn.return_w(i, tmpw_5d, w1[i]);
+    }
+    LOG(Message) << "Found v and w vectors for N =  " << N << std::endl;
+    for (unsigned int i = 0; i < N; i++)
+    {
+        v1[i] = gSnk * v1[i];
+    }
+    int ty;
+    for (unsigned int i = 0; i < N; i++)
+    {
+        for (unsigned int j = 0; j < N; j++)
+        {
+            mySliceInnerProductVector(MF_x, w1[i], v1[j], Tp);
+            mySliceInnerProductVector(MF_y, w1[j], v1[i], Tp);
+            for (unsigned int t = 0; t < nt; ++t)
+            {
+                for (unsigned int tx = 0; tx < nt; tx++)
+                {
+                    ty = (tx + t) % nt;
+                    tmp[t] += TensorRemove((MF_x[tx]) * (MF_y[ty]));
+                }
+            }
+        }
+        if (i % 10 == 0)
+        {
+            LOG(Message) << "MF for i = " << i << " of " << N << std::endl;
+        }
+    }
+    double NTinv = 1.0 / static_cast<double>(nt);
+    for (unsigned int t = 0; t < nt; ++t)
+    {
+        result.corr[t] = NTinv * tmp[t];
+    }
+
+    saveResult(par().output, "meson", result);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MContraction_A2AMeson_hpp_
diff --git a/extras/Hadrons/Modules/MContraction/MesonFieldGamma.cc b/extras/Hadrons/Modules/MContraction/MesonFieldGamma.cc
new file mode 100644
index 00000000..9218d587
--- /dev/null
+++ b/extras/Hadrons/Modules/MContraction/MesonFieldGamma.cc
@@ -0,0 +1,8 @@
+#include <Grid/Hadrons/Modules/MContraction/MesonFieldGamma.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MContraction;
+
+template class Grid::Hadrons::MContraction::TMesonFieldGamma<FIMPL>;
+template class Grid::Hadrons::MContraction::TMesonFieldGamma<ZFIMPL>;
\ No newline at end of file
diff --git a/extras/Hadrons/Modules/MContraction/MesonFieldGamma.hpp b/extras/Hadrons/Modules/MContraction/MesonFieldGamma.hpp
new file mode 100644
index 00000000..6128534c
--- /dev/null
+++ b/extras/Hadrons/Modules/MContraction/MesonFieldGamma.hpp
@@ -0,0 +1,269 @@
+#ifndef Hadrons_MContraction_MesonFieldGamma_hpp_
+#define Hadrons_MContraction_MesonFieldGamma_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+#include <Grid/Hadrons/AllToAllVectors.hpp>
+#include <Grid/Hadrons/AllToAllReduction.hpp>
+#include <Grid/Grid_Eigen_Dense.h>
+#include <fstream>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         MesonFieldGamma                                 *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MContraction)
+
+class MesonFieldPar : Serializable
+{
+  public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(MesonFieldPar,
+                                    int, Nl,
+                                    int, N,
+                                    int, Nblock,
+                                    std::string, A2A1,
+                                    std::string, A2A2,
+                                    std::string, gammas,
+                                    std::string, output);
+};
+
+template <typename FImpl>
+class TMesonFieldGamma : public Module<MesonFieldPar>
+{
+  public:
+    FERM_TYPE_ALIASES(FImpl, );
+    SOLVER_TYPE_ALIASES(FImpl, );
+
+    typedef A2AModesSchurDiagTwo<typename FImpl::FermionField, FMat, Solver> A2ABase;
+
+    class Result : Serializable
+    {
+      public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
+                                        Gamma::Algebra, gamma,
+                                        std::vector<std::vector<std::vector<ComplexD>>>, MesonField);
+    };
+
+  public:
+    // constructor
+    TMesonFieldGamma(const std::string name);
+    // destructor
+    virtual ~TMesonFieldGamma(void){};
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    virtual void parseGammaString(std::vector<Gamma::Algebra> &gammaList);
+    virtual void vectorOfWs(std::vector<FermionField> &w, int i, int Nblock, FermionField &tmpw_5d, std::vector<FermionField> &vec_w);
+    virtual void vectorOfVs(std::vector<FermionField> &v, int j, int Nblock, FermionField &tmpv_5d, std::vector<FermionField> &vec_v);
+    virtual void gammaMult(std::vector<FermionField> &v, Gamma gamma);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER(MesonFieldGamma, ARG(TMesonFieldGamma<FIMPL>), MContraction);
+MODULE_REGISTER(ZMesonFieldGamma, ARG(TMesonFieldGamma<ZFIMPL>), MContraction);
+
+/******************************************************************************
+*                  TMesonFieldGamma implementation                             *
+******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TMesonFieldGamma<FImpl>::TMesonFieldGamma(const std::string name)
+    : Module<MesonFieldPar>(name)
+{
+}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TMesonFieldGamma<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().A2A1 + "_class", par().A2A2 + "_class"};
+    in.push_back(par().A2A1 + "_w_high_4d");
+    in.push_back(par().A2A2 + "_v_high_4d");
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TMesonFieldGamma<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {};
+
+    return out;
+}
+
+template <typename FImpl>
+void TMesonFieldGamma<FImpl>::parseGammaString(std::vector<Gamma::Algebra> &gammaList)
+{
+    gammaList.clear();
+    // Determine gamma matrices to insert at source/sink.
+    if (par().gammas.compare("all") == 0)
+    {
+        // Do all contractions.
+        for (unsigned int i = 1; i < Gamma::nGamma; i += 2)
+        {
+            gammaList.push_back(((Gamma::Algebra)i));
+        }
+    }
+    else
+    {
+        // Parse individual contractions from input string.
+        gammaList = strToVec<Gamma::Algebra>(par().gammas);
+    }
+}
+
+template <typename FImpl>
+void TMesonFieldGamma<FImpl>::vectorOfWs(std::vector<FermionField> &w, int i, int Nblock, FermionField &tmpw_5d, std::vector<FermionField> &vec_w)
+{
+    for (unsigned int ni = 0; ni < Nblock; ni++)
+    {
+        vec_w[ni] = w[i + ni];
+    }
+}
+
+template <typename FImpl>
+void TMesonFieldGamma<FImpl>::vectorOfVs(std::vector<FermionField> &v, int j, int Nblock, FermionField &tmpv_5d, std::vector<FermionField> &vec_v)
+{
+    for (unsigned int nj = 0; nj < Nblock; nj++)
+    {
+        vec_v[nj] = v[j+nj];
+    }
+}
+
+template <typename FImpl>
+void TMesonFieldGamma<FImpl>::gammaMult(std::vector<FermionField> &v, Gamma gamma)
+{
+    int Nblock = v.size();
+    for (unsigned int nj = 0; nj < Nblock; nj++)
+    {
+        v[nj] = gamma * v[nj];
+    }
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TMesonFieldGamma<FImpl>::setup(void)
+{
+    int nt = env().getDim(Tp);
+    int N = par().N;
+    int Nblock = par().Nblock;
+
+    int Ls_ = env().getObjectLs(par().A2A1 + "_class");
+
+    envTmpLat(FermionField, "tmpv_5d", Ls_);
+    envTmpLat(FermionField, "tmpw_5d", Ls_);
+
+    envTmp(std::vector<FermionField>, "w", 1, N, FermionField(env().getGrid(1)));
+    envTmp(std::vector<FermionField>, "v", 1, N, FermionField(env().getGrid(1)));
+
+    envTmp(Eigen::MatrixXcd, "MF", 1, Eigen::MatrixXcd::Zero(nt, N * N));
+
+    envTmp(std::vector<FermionField>, "w_block", 1, Nblock, FermionField(env().getGrid(1)));
+    envTmp(std::vector<FermionField>, "v_block", 1, Nblock, FermionField(env().getGrid(1)));
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TMesonFieldGamma<FImpl>::execute(void)
+{
+    LOG(Message) << "Computing A2A meson field for gamma = " << par().gammas << ", taking w from " << par().A2A1 << " and v from " << par().A2A2 << std::endl;
+
+    int N = par().N;
+    int nt = env().getDim(Tp);
+    int Nblock = par().Nblock;
+
+    std::vector<Result> result;
+    std::vector<Gamma::Algebra> gammaResultList;
+    std::vector<Gamma> gammaList;
+
+    parseGammaString(gammaResultList);
+    result.resize(gammaResultList.size());
+
+    Gamma g5(Gamma::Algebra::Gamma5);
+    gammaList.resize(gammaResultList.size(), g5);
+
+    for (unsigned int i = 0; i < result.size(); ++i)
+    {
+        result[i].gamma = gammaResultList[i];
+        result[i].MesonField.resize(N, std::vector<std::vector<ComplexD>>(N, std::vector<ComplexD>(nt)));
+
+        Gamma gamma(gammaResultList[i]);
+        gammaList[i] = gamma;
+    }
+
+    auto &a2a1 = envGet(A2ABase, par().A2A1 + "_class");
+    auto &a2a2 = envGet(A2ABase, par().A2A2 + "_class");
+
+    envGetTmp(FermionField, tmpv_5d);
+    envGetTmp(FermionField, tmpw_5d);
+
+    envGetTmp(std::vector<FermionField>, v);
+    envGetTmp(std::vector<FermionField>, w);
+    LOG(Message) << "Finding v and w vectors for N =  " << N << std::endl;
+    for (int i = 0; i < N; i++)
+    {
+        a2a2.return_v(i, tmpv_5d, v[i]);
+        a2a1.return_w(i, tmpw_5d, w[i]);
+    }
+    LOG(Message) << "Found v and w vectors for N =  " << N << std::endl;
+
+    std::vector<std::vector<ComplexD>> MesonField_ij;
+    LOG(Message) << "Before blocked MFs, Nblock = " << Nblock << std::endl;
+    envGetTmp(std::vector<FermionField>, v_block);
+    envGetTmp(std::vector<FermionField>, w_block);
+    MesonField_ij.resize(Nblock * Nblock, std::vector<ComplexD>(nt));
+
+    envGetTmp(Eigen::MatrixXcd, MF);
+
+    LOG(Message) << "Before blocked MFs, Nblock = " << Nblock << std::endl;
+    for (unsigned int i = 0; i < N; i += Nblock)
+    {
+        vectorOfWs(w, i, Nblock, tmpw_5d, w_block);
+        for (unsigned int j = 0; j < N; j += Nblock)
+        {
+            vectorOfVs(v, j, Nblock, tmpv_5d, v_block);
+            for (unsigned int k = 0; k < result.size(); k++)
+            {
+                gammaMult(v_block, gammaList[k]);
+                sliceInnerProductMesonField(MesonField_ij, w_block, v_block, Tp);
+                for (unsigned int nj = 0; nj < Nblock; nj++)
+                {
+                    for (unsigned int ni = 0; ni < Nblock; ni++)
+                    {
+                        MF.col((i + ni) + (j + nj) * N) = Eigen::VectorXcd::Map(&MesonField_ij[nj * Nblock + ni][0], MesonField_ij[nj * Nblock + ni].size());
+                    }
+                }
+            }
+        }
+        if (i % 10 == 0)
+        {
+            LOG(Message) << "MF for i = " << i << " of " << N << std::endl;
+        }
+    }
+    LOG(Message) << "Before Global sum, Nblock = " << Nblock << std::endl;
+    v_block[0]._grid->GlobalSumVector(MF.data(), MF.size());
+    LOG(Message) << "After Global sum, Nblock = " << Nblock << std::endl;
+    for (unsigned int i = 0; i < N; i++)
+    {
+        for (unsigned int j = 0; j < N; j++)
+        {
+            for (unsigned int k = 0; k < result.size(); k++)
+            {
+                for (unsigned int t = 0; t < nt; t++)
+                {
+                    result[k].MesonField[i][j][t] = MF.col(i + N * j)[t];
+                }
+            }
+        }
+    }
+    saveResult(par().output, "meson", result);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MContraction_MesonFieldGm_hpp_
diff --git a/extras/Hadrons/Modules/MSolver/A2AVectors.cc b/extras/Hadrons/Modules/MSolver/A2AVectors.cc
new file mode 100644
index 00000000..f72f405d
--- /dev/null
+++ b/extras/Hadrons/Modules/MSolver/A2AVectors.cc
@@ -0,0 +1,8 @@
+#include <Grid/Hadrons/Modules/MSolver/A2AVectors.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MSolver;
+
+template class Grid::Hadrons::MSolver::TA2AVectors<FIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>;
+template class Grid::Hadrons::MSolver::TA2AVectors<ZFIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>;
diff --git a/extras/Hadrons/Modules/MSolver/A2AVectors.hpp b/extras/Hadrons/Modules/MSolver/A2AVectors.hpp
new file mode 100644
index 00000000..9481f268
--- /dev/null
+++ b/extras/Hadrons/Modules/MSolver/A2AVectors.hpp
@@ -0,0 +1,253 @@
+#ifndef Hadrons_MSolver_A2AVectors_hpp_
+#define Hadrons_MSolver_A2AVectors_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+#include <Grid/Hadrons/Solver.hpp>
+#include <Grid/Hadrons/EigenPack.hpp>
+#include <Grid/Hadrons/AllToAllVectors.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         A2AVectors                                 *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MSolver)
+
+class A2AVectorsPar: Serializable
+{
+public:
+  GRID_SERIALIZABLE_CLASS_MEMBERS(A2AVectorsPar,
+                                  bool, return_5d,
+                                  int, Nl,
+                                  int, N,
+                                  std::vector<std::string>, sources,
+                                  std::string, action,
+                                  std::string, eigenPack,
+                                  std::string, solver);
+};
+
+template <typename FImpl, int nBasis>
+class TA2AVectors : public Module<A2AVectorsPar>
+{
+  public:
+    FERM_TYPE_ALIASES(FImpl,);
+    SOLVER_TYPE_ALIASES(FImpl,);
+
+    typedef FermionEigenPack<FImpl> EPack;
+    typedef CoarseFermionEigenPack<FImpl, nBasis> CoarseEPack;
+
+    typedef A2AModesSchurDiagTwo<typename FImpl::FermionField, FMat, Solver> A2ABase;
+
+  public:
+    // constructor
+    TA2AVectors(const std::string name);
+    // destructor
+    virtual ~TA2AVectors(void) {};
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getReference(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+
+  private:
+    unsigned int Ls_;
+    std::string className_;
+};
+
+MODULE_REGISTER_TMP(A2AVectors, ARG(TA2AVectors<FIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>), MSolver);
+MODULE_REGISTER_TMP(ZA2AVectors, ARG(TA2AVectors<ZFIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>), MSolver);
+
+/******************************************************************************
+ *                 TA2AVectors implementation                             *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl, int nBasis>
+TA2AVectors<FImpl, nBasis>::TA2AVectors(const std::string name)
+: Module<A2AVectorsPar>(name)
+, className_ (name + "_class")
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl, int nBasis>
+std::vector<std::string> TA2AVectors<FImpl, nBasis>::getInput(void)
+{
+    int Nl = par().Nl;
+    std::string sub_string = "";
+    if (Nl > 0) sub_string = "_subtract";
+    std::vector<std::string> in = {par().solver + sub_string};
+
+    int n = par().sources.size();
+
+    for (unsigned int t = 0; t < n; t += 1)
+    {
+        in.push_back(par().sources[t]);
+    }
+
+    return in;
+}
+
+template <typename FImpl, int nBasis>
+std::vector<std::string> TA2AVectors<FImpl, nBasis>::getReference(void)
+{
+    std::vector<std::string> ref = {par().action};
+
+    if (!par().eigenPack.empty())
+    {
+        ref.push_back(par().eigenPack);
+    }
+
+    return ref;
+}
+
+template <typename FImpl, int nBasis>
+std::vector<std::string> TA2AVectors<FImpl, nBasis>::getOutput(void)
+{
+    std::vector<std::string> out = {getName(), className_,
+                                    getName() + "_w_high_5d", getName() + "_v_high_5d",
+                                    getName() + "_w_high_4d", getName() + "_v_high_4d"};
+
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl, int nBasis>
+void TA2AVectors<FImpl, nBasis>::setup(void)
+{
+    int N = par().N;
+    int Nl = par().Nl;
+    int Nh = N - Nl;
+    bool return_5d = par().return_5d;
+    int Ls;
+
+    std::string sub_string = "";
+    if (Nl > 0) sub_string = "_subtract";
+    auto &solver = envGet(Solver, par().solver + sub_string);
+    Ls = env().getObjectLs(par().solver + sub_string);
+
+    auto &action = envGet(FMat, par().action);
+
+    envTmpLat(FermionField, "ferm_src", Ls);
+    envTmpLat(FermionField, "unphys_ferm", Ls);
+    envTmpLat(FermionField, "tmp");
+    envTmpLat(FermionField, "tmp2");
+
+    std::vector<FermionField> *evec;
+    const std::vector<RealD> *eval;
+
+    if (Nl > 0)
+    {
+        // Low modes
+        auto &epack = envGet(EPack, par().eigenPack);
+
+        LOG(Message) << "Creating a2a vectors " << getName() <<
+                     " using eigenpack '" << par().eigenPack << "' ("
+                     << epack.evec.size() << " modes)" <<
+                     " and " << Nh << " high modes." << std::endl;
+        evec = &epack.evec;
+        eval = &epack.eval;
+    }
+    else
+    {
+        LOG(Message) << "Creating a2a vectors " << getName() <<
+                     " using " << Nh << " high modes only." << std::endl;
+    }
+
+    int size_5d = 1;
+    if (return_5d) size_5d = Nh;
+    envCreate(std::vector<FermionField>, getName() + "_w_high_5d", Ls, size_5d, FermionField(env().getGrid(Ls)));
+    envCreate(std::vector<FermionField>, getName() + "_v_high_5d", Ls, size_5d, FermionField(env().getGrid(Ls)));
+    envCreate(std::vector<FermionField>, getName() + "_w_high_4d", 1, Nh, FermionField(env().getGrid(1)));
+    envCreate(std::vector<FermionField>, getName() + "_v_high_4d", 1, Nh, FermionField(env().getGrid(1)));
+
+    auto &w_high_5d = envGet(std::vector<FermionField>, getName() + "_w_high_5d");
+    auto &v_high_5d = envGet(std::vector<FermionField>, getName() + "_v_high_5d");
+    auto &w_high_4d = envGet(std::vector<FermionField>, getName() + "_w_high_4d");
+    auto &v_high_4d = envGet(std::vector<FermionField>, getName() + "_v_high_4d");
+
+    envCreate(A2ABase, className_, Ls,
+                     evec, eval,
+                     action,
+                     solver,
+                     w_high_5d, v_high_5d,
+                     w_high_4d, v_high_4d,
+                     Nl, Nh,
+                     return_5d);
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl, int nBasis>
+void TA2AVectors<FImpl, nBasis>::execute(void)
+{
+    auto &action = envGet(FMat, par().action);
+
+    int Nt = env().getDim(Tp);
+    int Nc = FImpl::Dimension;
+    int Ls_;
+    int Nl = par().Nl;
+
+    std::string sub_string = "";
+    if (Nl > 0) sub_string = "_subtract";
+    Ls_ = env().getObjectLs(par().solver + sub_string);
+
+    auto &a2areturn = envGet(A2ABase, className_);
+
+    // High modes
+    auto sources = par().sources;
+    int Nsrc = par().sources.size();
+
+    envGetTmp(FermionField, ferm_src);
+    envGetTmp(FermionField, unphys_ferm);
+    envGetTmp(FermionField, tmp);
+    envGetTmp(FermionField, tmp2);
+
+    int N_count = 0;
+    for (unsigned int s = 0; s < Ns; ++s)
+        for (unsigned int c = 0; c < Nc; ++c)
+        for (unsigned int T = 0; T < Nsrc; T++)
+        {
+            auto &prop_src = envGet(PropagatorField, sources[T]);
+            LOG(Message) << "A2A src for s = " << s << " , c = " << c << ", T = " << T << std::endl;
+            // source conversion for 4D sources
+            if (!env().isObject5d(sources[T]))
+            {
+                if (Ls_ == 1)
+                {
+                    PropToFerm<FImpl>(ferm_src, prop_src, s, c);
+                    tmp = ferm_src;
+                }
+                else
+                {
+                    PropToFerm<FImpl>(tmp, prop_src, s, c);
+                    action.ImportPhysicalFermionSource(tmp, ferm_src);
+                    action.ImportUnphysicalFermion(tmp, unphys_ferm);
+                }
+            }
+            // source conversion for 5D sources
+            else
+            {
+                if (Ls_ != env().getObjectLs(sources[T]))
+                {
+                    HADRONS_ERROR(Size, "Ls mismatch between quark action and source");
+                }
+                else
+                {
+                    PropToFerm<FImpl>(ferm_src, prop_src, s, c);
+                    action.ExportPhysicalFermionSolution(ferm_src, tmp);
+                    unphys_ferm = ferm_src;
+                }
+            }
+            LOG(Message) << "a2areturn.high_modes Ncount = " << N_count << std::endl;
+            a2areturn.high_modes(ferm_src, unphys_ferm, tmp, N_count);
+            N_count++;
+        }
+}
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MSolver_A2AVectors_hpp_
diff --git a/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp b/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp
index 31be621f..347371a3 100644
--- a/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp
+++ b/extras/Hadrons/Modules/MSolver/RBPrecCG.hpp
@@ -118,7 +118,7 @@ std::vector<std::string> TRBPrecCG<FImpl, nBasis>::getReference(void)
 template <typename FImpl, int nBasis>
 std::vector<std::string> TRBPrecCG<FImpl, nBasis>::getOutput(void)
 {
-    std::vector<std::string> out = {getName()};
+    std::vector<std::string> out = {getName(), getName() + "_subtract"};
     
     return out;
 }
@@ -168,19 +168,22 @@ void TRBPrecCG<FImpl, nBasis>::setup(void)
             guesser.reset(new FineGuesser(epack.evec, epack.eval));
         }
     }
-    auto solver = [&mat, guesser, this](FermionField &sol, 
-                                        const FermionField &source)
-    {
-        ConjugateGradient<FermionField>           cg(par().residual, 
-                                                     par().maxIteration);
-        HADRONS_DEFAULT_SCHUR_SOLVE<FermionField> schurSolver(cg);
-        
-        schurSolver(mat, source, sol, *guesser);
+    auto makeSolver = [&mat, guesser, this](bool subGuess) {
+        return [&mat, guesser, subGuess, this](FermionField &sol,
+                                     const FermionField &source) {
+            ConjugateGradient<FermionField> cg(par().residual,
+                                               par().maxIteration);
+            HADRONS_DEFAULT_SCHUR_SOLVE<FermionField> schurSolver(cg);
+            schurSolver.subtractGuess(subGuess);
+            schurSolver(mat, source, sol, *guesser);
+        };
     };
+    auto solver = makeSolver(false);
     envCreate(Solver, getName(), Ls, solver, mat);
+    auto solver_subtract = makeSolver(true);
+    envCreate(Solver, getName() + "_subtract", Ls, solver_subtract, mat);
 }
 
-
 // execution ///////////////////////////////////////////////////////////////////
 template <typename FImpl, int nBasis>
 void TRBPrecCG<FImpl, nBasis>::execute(void)
diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc
index 58b082d0..b321f2b7 100644
--- a/extras/Hadrons/modules.inc
+++ b/extras/Hadrons/modules.inc
@@ -1,115 +1,121 @@
 modules_cc =\
-  Modules/MContraction/WeakHamiltonianEye.cc \
-  Modules/MContraction/Baryon.cc \
-  Modules/MContraction/Meson.cc \
-  Modules/MContraction/WeakNeutral4ptDisc.cc \
-  Modules/MContraction/WeakHamiltonianNonEye.cc \
-  Modules/MContraction/WardIdentity.cc \
-  Modules/MContraction/DiscLoop.cc \
-  Modules/MContraction/Gamma3pt.cc \
-  Modules/MFermion/FreeProp.cc \
-  Modules/MFermion/GaugeProp.cc \
-  Modules/MSource/Point.cc \
-  Modules/MSource/Wall.cc \
   Modules/MSource/SeqConserved.cc \
   Modules/MSource/SeqGamma.cc \
+  Modules/MSource/Point.cc \
   Modules/MSource/Z2.cc \
+  Modules/MSource/Wall.cc \
   Modules/MSink/Point.cc \
   Modules/MSink/Smear.cc \
-  Modules/MSolver/RBPrecCG.cc \
-  Modules/MSolver/LocalCoherenceLanczos.cc \
-  Modules/MGauge/StoutSmearing.cc \
-  Modules/MGauge/Unit.cc \
-  Modules/MGauge/UnitEm.cc \
-  Modules/MGauge/StochEm.cc \
-  Modules/MGauge/Random.cc \
-  Modules/MGauge/FundtoHirep.cc \
-  Modules/MUtilities/TestSeqGamma.cc \
-  Modules/MUtilities/TestSeqConserved.cc \
-  Modules/MLoop/NoiseLoop.cc \
-  Modules/MScalar/FreeProp.cc \
-  Modules/MScalar/VPCounterTerms.cc \
-  Modules/MScalar/ChargedProp.cc \
-  Modules/MScalar/ScalarVP.cc \
-  Modules/MAction/Wilson.cc \
+  Modules/MIO/LoadNersc.cc \
+  Modules/MIO/LoadEigenPack.cc \
+  Modules/MIO/LoadCoarseEigenPack.cc \
+  Modules/MIO/LoadBinary.cc \
+  Modules/MScalarSUN/TwoPointNPR.cc \
+  Modules/MScalarSUN/TrPhi.cc \
+  Modules/MScalarSUN/StochFreeField.cc \
+  Modules/MScalarSUN/TrMag.cc \
+  Modules/MScalarSUN/Div.cc \
+  Modules/MScalarSUN/ShiftProbe.cc \
+  Modules/MScalarSUN/Grad.cc \
+  Modules/MScalarSUN/TwoPoint.cc \
+  Modules/MScalarSUN/TimeMomProbe.cc \
+  Modules/MScalarSUN/EMT.cc \
+  Modules/MScalarSUN/TransProj.cc \
+  Modules/MScalarSUN/TrKinetic.cc \
+  Modules/MAction/DWF.cc \
   Modules/MAction/MobiusDWF.cc \
   Modules/MAction/ZMobiusDWF.cc \
-  Modules/MAction/WilsonClover.cc \
-  Modules/MAction/DWF.cc \
+  Modules/MAction/Wilson.cc \
   Modules/MAction/ScaledDWF.cc \
-  Modules/MScalarSUN/TrPhi.cc \
-  Modules/MScalarSUN/Grad.cc \
-  Modules/MScalarSUN/TimeMomProbe.cc \
-  Modules/MScalarSUN/TrMag.cc \
-  Modules/MScalarSUN/TrKinetic.cc \
-  Modules/MScalarSUN/EMT.cc \
-  Modules/MScalarSUN/ShiftProbe.cc \
-  Modules/MScalarSUN/TransProj.cc \
-  Modules/MScalarSUN/StochFreeField.cc \
-  Modules/MScalarSUN/TwoPoint.cc \
-  Modules/MScalarSUN/TwoPointNPR.cc \
-  Modules/MScalarSUN/Div.cc \
-  Modules/MIO/LoadEigenPack.cc \
-  Modules/MIO/LoadBinary.cc \
-  Modules/MIO/LoadNersc.cc \
-  Modules/MIO/LoadCoarseEigenPack.cc
+  Modules/MAction/WilsonClover.cc \
+  Modules/MContraction/WeakHamiltonianNonEye.cc \
+  Modules/MContraction/WardIdentity.cc \
+  Modules/MContraction/WeakHamiltonianEye.cc \
+  Modules/MContraction/DiscLoop.cc \
+  Modules/MContraction/A2AMeson.cc \
+  Modules/MContraction/Baryon.cc \
+  Modules/MContraction/MesonFieldGamma.cc \
+  Modules/MContraction/Gamma3pt.cc \
+  Modules/MContraction/WeakNeutral4ptDisc.cc \
+  Modules/MContraction/Meson.cc \
+  Modules/MScalar/VPCounterTerms.cc \
+  Modules/MScalar/ChargedProp.cc \
+  Modules/MScalar/FreeProp.cc \
+  Modules/MScalar/ScalarVP.cc \
+  Modules/MUtilities/TestSeqGamma.cc \
+  Modules/MUtilities/TestSeqConserved.cc \
+  Modules/MFermion/FreeProp.cc \
+  Modules/MFermion/GaugeProp.cc \
+  Modules/MSolver/RBPrecCG.cc \
+  Modules/MSolver/LocalCoherenceLanczos.cc \
+  Modules/MSolver/A2AVectors.cc \
+  Modules/MLoop/NoiseLoop.cc \
+  Modules/MGauge/Unit.cc \
+  Modules/MGauge/Random.cc \
+  Modules/MGauge/UnitEm.cc \
+  Modules/MGauge/StochEm.cc \
+  Modules/MGauge/FundtoHirep.cc \
+  Modules/MGauge/StoutSmearing.cc
 
 modules_hpp =\
-  Modules/MContraction/Baryon.hpp \
-  Modules/MContraction/Meson.hpp \
-  Modules/MContraction/WeakHamiltonian.hpp \
-  Modules/MContraction/WeakHamiltonianNonEye.hpp \
-  Modules/MContraction/DiscLoop.hpp \
-  Modules/MContraction/WeakNeutral4ptDisc.hpp \
-  Modules/MContraction/Gamma3pt.hpp \
-  Modules/MContraction/WardIdentity.hpp \
-  Modules/MContraction/WeakHamiltonianEye.hpp \
-  Modules/MFermion/FreeProp.hpp \
-  Modules/MFermion/GaugeProp.hpp \
+  Modules/MSource/SeqConserved.hpp \
   Modules/MSource/SeqGamma.hpp \
+  Modules/MSource/Z2.hpp \
   Modules/MSource/Point.hpp \
   Modules/MSource/Wall.hpp \
-  Modules/MSource/Z2.hpp \
-  Modules/MSource/SeqConserved.hpp \
   Modules/MSink/Smear.hpp \
   Modules/MSink/Point.hpp \
-  Modules/MSolver/LocalCoherenceLanczos.hpp \
-  Modules/MSolver/RBPrecCG.hpp \
-  Modules/MGauge/UnitEm.hpp \
-  Modules/MGauge/StoutSmearing.hpp \
-  Modules/MGauge/Unit.hpp \
-  Modules/MGauge/Random.hpp \
-  Modules/MGauge/FundtoHirep.hpp \
-  Modules/MGauge/StochEm.hpp \
-  Modules/MUtilities/TestSeqGamma.hpp \
-  Modules/MUtilities/TestSeqConserved.hpp \
-  Modules/MLoop/NoiseLoop.hpp \
-  Modules/MScalar/FreeProp.hpp \
-  Modules/MScalar/VPCounterTerms.hpp \
-  Modules/MScalar/ScalarVP.hpp \
-  Modules/MScalar/Scalar.hpp \
-  Modules/MScalar/ChargedProp.hpp \
-  Modules/MAction/DWF.hpp \
-  Modules/MAction/MobiusDWF.hpp \
-  Modules/MAction/Wilson.hpp \
-  Modules/MAction/WilsonClover.hpp \
-  Modules/MAction/ZMobiusDWF.hpp \
-  Modules/MAction/ScaledDWF.hpp \
-  Modules/MScalarSUN/StochFreeField.hpp \
+  Modules/MIO/LoadBinary.hpp \
+  Modules/MIO/LoadEigenPack.hpp \
+  Modules/MIO/LoadCoarseEigenPack.hpp \
+  Modules/MIO/LoadNersc.hpp \
+  Modules/MScalarSUN/Utils.hpp \
+  Modules/MScalarSUN/Grad.hpp \
+  Modules/MScalarSUN/TrPhi.hpp \
   Modules/MScalarSUN/TwoPointNPR.hpp \
+  Modules/MScalarSUN/TwoPoint.hpp \
+  Modules/MScalarSUN/TransProj.hpp \
+  Modules/MScalarSUN/TrKinetic.hpp \
+  Modules/MScalarSUN/StochFreeField.hpp \
   Modules/MScalarSUN/ShiftProbe.hpp \
-  Modules/MScalarSUN/Div.hpp \
   Modules/MScalarSUN/TimeMomProbe.hpp \
+  Modules/MScalarSUN/Div.hpp \
   Modules/MScalarSUN/TrMag.hpp \
   Modules/MScalarSUN/EMT.hpp \
-  Modules/MScalarSUN/TwoPoint.hpp \
-  Modules/MScalarSUN/TrPhi.hpp \
-  Modules/MScalarSUN/Utils.hpp \
-  Modules/MScalarSUN/TransProj.hpp \
-  Modules/MScalarSUN/Grad.hpp \
-  Modules/MScalarSUN/TrKinetic.hpp \
-  Modules/MIO/LoadEigenPack.hpp \
-  Modules/MIO/LoadNersc.hpp \
-  Modules/MIO/LoadCoarseEigenPack.hpp \
-  Modules/MIO/LoadBinary.hpp
+  Modules/MAction/ZMobiusDWF.hpp \
+  Modules/MAction/ScaledDWF.hpp \
+  Modules/MAction/Wilson.hpp \
+  Modules/MAction/WilsonClover.hpp \
+  Modules/MAction/MobiusDWF.hpp \
+  Modules/MAction/DWF.hpp \
+  Modules/MContraction/WeakHamiltonian.hpp \
+  Modules/MContraction/DiscLoop.hpp \
+  Modules/MContraction/Meson.hpp \
+  Modules/MContraction/WardIdentity.hpp \
+  Modules/MContraction/WeakHamiltonianEye.hpp \
+  Modules/MContraction/Gamma3pt.hpp \
+  Modules/MContraction/WeakHamiltonianNonEye.hpp \
+  Modules/MContraction/MesonFieldGamma.hpp \
+  Modules/MContraction/Baryon.hpp \
+  Modules/MContraction/WeakNeutral4ptDisc.hpp \
+  Modules/MContraction/A2AMeson.hpp \
+  Modules/MScalar/ScalarVP.hpp \
+  Modules/MScalar/Scalar.hpp \
+  Modules/MScalar/FreeProp.hpp \
+  Modules/MScalar/ChargedProp.hpp \
+  Modules/MScalar/VPCounterTerms.hpp \
+  Modules/MUtilities/TestSeqConserved.hpp \
+  Modules/MUtilities/TestSeqGamma.hpp \
+  Modules/MFermion/FreeProp.hpp \
+  Modules/MFermion/GaugeProp.hpp \
+  Modules/MSolver/A2AVectors.hpp \
+  Modules/MSolver/RBPrecCG.hpp \
+  Modules/MSolver/LocalCoherenceLanczos.hpp \
+  Modules/MLoop/NoiseLoop.hpp \
+  Modules/MGauge/StoutSmearing.hpp \
+  Modules/MGauge/StochEm.hpp \
+  Modules/MGauge/FundtoHirep.hpp \
+  Modules/MGauge/Unit.hpp \
+  Modules/MGauge/Random.hpp \
+  Modules/MGauge/UnitEm.hpp
 
diff --git a/lib/algorithms/iterative/Deflation.h b/lib/algorithms/iterative/Deflation.h
index 316afe90..7c3e8e57 100644
--- a/lib/algorithms/iterative/Deflation.h
+++ b/lib/algorithms/iterative/Deflation.h
@@ -63,7 +63,7 @@ public:
 
   DeflatedGuesser(const std::vector<Field> & _evec,const std::vector<RealD> & _eval) : evec(_evec), eval(_eval) {};
 
-  virtual void operator()(const Field &src,Field &guess) { 
+  virtual void operator()(const Field &src,Field &guess) {
     guess = zero;
     assert(evec.size()==eval.size());
     auto N = evec.size();
@@ -71,6 +71,7 @@ public:
       const Field& tmp = evec[i];
       axpy(guess,TensorRemove(innerProduct(tmp,src)) / eval[i],tmp,guess);
     }
+    guess.checkerboard = src.checkerboard;
   }
 };
 
@@ -101,6 +102,7 @@ public:
       axpy(guess_coarse,TensorRemove(innerProduct(tmp,src_coarse)) / eval_coarse[i],tmp,guess_coarse);
     }
     blockPromote(guess_coarse,guess,subspace);
+    guess.checkerboard = src.checkerboard;
   };
 };
 
diff --git a/lib/algorithms/iterative/SchurRedBlack.h b/lib/algorithms/iterative/SchurRedBlack.h
index 091330b2..5abc4d9f 100644
--- a/lib/algorithms/iterative/SchurRedBlack.h
+++ b/lib/algorithms/iterative/SchurRedBlack.h
@@ -95,16 +95,26 @@ namespace Grid {
   private:
     OperatorFunction<Field> & _HermitianRBSolver;
     int CBfactorise;
+    bool subGuess;
   public:
 
     /////////////////////////////////////////////////////
     // Wrap the usual normal equations Schur trick
     /////////////////////////////////////////////////////
-  SchurRedBlackStaggeredSolve(OperatorFunction<Field> &HermitianRBSolver)  :
+  SchurRedBlackStaggeredSolve(OperatorFunction<Field> &HermitianRBSolver, const bool initSubGuess = false)  :
      _HermitianRBSolver(HermitianRBSolver) 
     { 
       CBfactorise=0;
+      subtractGuess(initSubGuess);
     };
+    void subtractGuess(const bool initSubGuess)
+    {
+      subGuess = initSubGuess;
+    }
+    bool isSubtractGuess(void)
+    {
+      return subGuess;
+    }
 
     template<class Matrix>
     void operator() (Matrix & _Matrix,const Field &in, Field &out){
@@ -150,9 +160,12 @@ namespace Grid {
       // Call the red-black solver
       //////////////////////////////////////////////////////////////
       std::cout<<GridLogMessage << "SchurRedBlackStaggeredSolver calling the Mpc solver" <<std::endl;
-      guess(src_o,sol_o);
+      guess(src_o, sol_o);
+      Mtmp = sol_o;
       _HermitianRBSolver(_HermOpEO,src_o,sol_o);  assert(sol_o.checkerboard==Odd);
       std::cout<<GridLogMessage << "SchurRedBlackStaggeredSolver called  the Mpc solver" <<std::endl;
+      // Fionn A2A boolean behavioural control
+      if (subGuess)        sol_o = sol_o-Mtmp;
 
       ///////////////////////////////////////////////////
       // sol_e = M_ee^-1 * ( src_e - Meo sol_o )...
@@ -167,11 +180,15 @@ namespace Grid {
       std::cout<<GridLogMessage << "SchurRedBlackStaggeredSolver inserted solution" <<std::endl;
 
       // Verify the unprec residual
-      _Matrix.M(out,resid); 
-      resid = resid-in;
-      RealD ns = norm2(in);
-      RealD nr = norm2(resid);
-      std::cout<<GridLogMessage << "SchurRedBlackStaggered solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
+      if ( ! subGuess ) {
+        _Matrix.M(out,resid); 
+        resid = resid-in;
+        RealD ns = norm2(in);
+        RealD nr = norm2(resid);
+        std::cout<<GridLogMessage << "SchurRedBlackStaggered solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
+      } else {
+        std::cout << GridLogMessage << "Guess subtracted after solve." << std::endl;
+      }
     }     
   };
   template<class Field> using SchurRedBlackStagSolve = SchurRedBlackStaggeredSolve<Field>;
@@ -184,15 +201,25 @@ namespace Grid {
   private:
     OperatorFunction<Field> & _HermitianRBSolver;
     int CBfactorise;
+    bool subGuess;
   public:
 
     /////////////////////////////////////////////////////
     // Wrap the usual normal equations Schur trick
     /////////////////////////////////////////////////////
-  SchurRedBlackDiagMooeeSolve(OperatorFunction<Field> &HermitianRBSolver,int cb=0)  :  _HermitianRBSolver(HermitianRBSolver) 
+  SchurRedBlackDiagMooeeSolve(OperatorFunction<Field> &HermitianRBSolver,int cb=0, const bool initSubGuess = false)  :  _HermitianRBSolver(HermitianRBSolver) 
   { 
     CBfactorise=cb;
+    subtractGuess(initSubGuess);
   };
+    void subtractGuess(const bool initSubGuess)
+    {
+      subGuess = initSubGuess;
+    }
+    bool isSubtractGuess(void)
+    {
+      return subGuess;
+    }
     template<class Matrix>
     void operator() (Matrix & _Matrix,const Field &in, Field &out){
       ZeroGuesser<Field> guess;
@@ -236,7 +263,10 @@ namespace Grid {
       //////////////////////////////////////////////////////////////
       std::cout<<GridLogMessage << "SchurRedBlack solver calling the MpcDagMp solver" <<std::endl;
       guess(src_o,sol_o);
+      Mtmp = sol_o;
       _HermitianRBSolver(_HermOpEO,src_o,sol_o);  assert(sol_o.checkerboard==Odd);
+      // Fionn A2A boolean behavioural control
+      if (subGuess)        sol_o = sol_o-Mtmp;
 
       ///////////////////////////////////////////////////
       // sol_e = M_ee^-1 * ( src_e - Meo sol_o )...
@@ -249,12 +279,16 @@ namespace Grid {
       setCheckerboard(out,sol_o); assert(  sol_o.checkerboard ==Odd );
 
       // Verify the unprec residual
-      _Matrix.M(out,resid); 
-      resid = resid-in;
-      RealD ns = norm2(in);
-      RealD nr = norm2(resid);
+      if ( ! subGuess ) {
+        _Matrix.M(out,resid); 
+        resid = resid-in;
+        RealD ns = norm2(in);
+        RealD nr = norm2(resid);
 
-      std::cout<<GridLogMessage << "SchurRedBlackDiagMooee solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
+        std::cout<<GridLogMessage << "SchurRedBlackDiagMooee solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
+      } else {
+        std::cout << GridLogMessage << "Guess subtracted after solve." << std::endl;
+      }
     }     
   };
 
@@ -267,16 +301,26 @@ namespace Grid {
   private:
     OperatorFunction<Field> & _HermitianRBSolver;
     int CBfactorise;
+    bool subGuess;
   public:
 
     /////////////////////////////////////////////////////
     // Wrap the usual normal equations Schur trick
     /////////////////////////////////////////////////////
-  SchurRedBlackDiagTwoSolve(OperatorFunction<Field> &HermitianRBSolver)  :
+  SchurRedBlackDiagTwoSolve(OperatorFunction<Field> &HermitianRBSolver, const bool initSubGuess = false)  :
      _HermitianRBSolver(HermitianRBSolver) 
     { 
-      CBfactorise=0;
+      CBfactorise = 0;
+      subtractGuess(initSubGuess);
     };
+    void subtractGuess(const bool initSubGuess)
+    {
+      subGuess = initSubGuess;
+    }
+    bool isSubtractGuess(void)
+    {
+      return subGuess;
+    }
 
     template<class Matrix>
     void operator() (Matrix & _Matrix,const Field &in, Field &out){
@@ -322,8 +366,11 @@ namespace Grid {
       std::cout<<GridLogMessage << "SchurRedBlack solver calling the MpcDagMp solver" <<std::endl;
 //      _HermitianRBSolver(_HermOpEO,src_o,sol_o);  assert(sol_o.checkerboard==Odd);
       guess(src_o,tmp);
+      Mtmp = tmp;
       _HermitianRBSolver(_HermOpEO,src_o,tmp);  assert(tmp.checkerboard==Odd);
-      _Matrix.MooeeInv(tmp,sol_o);        assert(  sol_o.checkerboard   ==Odd);
+      // Fionn A2A boolean behavioural control
+      if (subGuess)      tmp = tmp-Mtmp;
+      _Matrix.MooeeInv(tmp,sol_o);       assert(  sol_o.checkerboard   ==Odd);
 
       ///////////////////////////////////////////////////
       // sol_e = M_ee^-1 * ( src_e - Meo sol_o )...
@@ -336,12 +383,16 @@ namespace Grid {
       setCheckerboard(out,sol_o); assert(  sol_o.checkerboard ==Odd );
 
       // Verify the unprec residual
-      _Matrix.M(out,resid); 
-      resid = resid-in;
-      RealD ns = norm2(in);
-      RealD nr = norm2(resid);
+      if ( ! subGuess ) {
+        _Matrix.M(out,resid); 
+        resid = resid-in;
+        RealD ns = norm2(in);
+        RealD nr = norm2(resid);
 
-      std::cout<<GridLogMessage << "SchurRedBlackDiagTwo solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
+        std::cout<<GridLogMessage << "SchurRedBlackDiagTwo solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
+      } else {
+        std::cout << GridLogMessage << "Guess subtracted after solve." << std::endl;
+      }
     }     
   };
   ///////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -352,16 +403,26 @@ namespace Grid {
   private:
     LinearFunction<Field> & _HermitianRBSolver;
     int CBfactorise;
+    bool subGuess;
   public:
 
     /////////////////////////////////////////////////////
     // Wrap the usual normal equations Schur trick
     /////////////////////////////////////////////////////
-  SchurRedBlackDiagTwoMixed(LinearFunction<Field> &HermitianRBSolver)  :
+  SchurRedBlackDiagTwoMixed(LinearFunction<Field> &HermitianRBSolver, const bool initSubGuess = false)  :
      _HermitianRBSolver(HermitianRBSolver) 
     { 
       CBfactorise=0;
+      subtractGuess(initSubGuess);
     };
+    void subtractGuess(const bool initSubGuess)
+    {
+      subGuess = initSubGuess;
+    }
+    bool isSubtractGuess(void)
+    {
+      return subGuess;
+    }
 
     template<class Matrix>
     void operator() (Matrix & _Matrix,const Field &in, Field &out){
@@ -408,7 +469,10 @@ namespace Grid {
 //      _HermitianRBSolver(_HermOpEO,src_o,sol_o);  assert(sol_o.checkerboard==Odd);
 //      _HermitianRBSolver(_HermOpEO,src_o,tmp);  assert(tmp.checkerboard==Odd);
       guess(src_o,tmp);
-      _HermitianRBSolver(src_o,tmp);  assert(tmp.checkerboard==Odd);
+      Mtmp = tmp;
+      _HermitianRBSolver(_HermOpEO,src_o,tmp);  assert(tmp.checkerboard==Odd);
+      // Fionn A2A boolean behavioural control
+      if (subGuess)      tmp = tmp-Mtmp;
       _Matrix.MooeeInv(tmp,sol_o);        assert(  sol_o.checkerboard   ==Odd);
 
       ///////////////////////////////////////////////////
@@ -422,12 +486,16 @@ namespace Grid {
       setCheckerboard(out,sol_o); assert(  sol_o.checkerboard ==Odd );
 
       // Verify the unprec residual
-      _Matrix.M(out,resid); 
-      resid = resid-in;
-      RealD ns = norm2(in);
-      RealD nr = norm2(resid);
+      if ( ! subGuess ) {
+        _Matrix.M(out,resid); 
+        resid = resid-in;
+        RealD ns = norm2(in);
+        RealD nr = norm2(resid);
 
-      std::cout<<GridLogMessage << "SchurRedBlackDiagTwo solver true unprec resid "<< std::sqrt(nr/ns) <<" nr "<< nr <<" ns "<<ns << std::endl;
+        std::cout << GridLogMessage << "SchurRedBlackDiagTwo solver true unprec resid " << std::sqrt(nr / ns) << " nr " << nr << " ns " << ns << std::endl;
+      } else {
+        std::cout << GridLogMessage << "Guess subtracted after solve." << std::endl;
+      }
     }     
   };
 
diff --git a/lib/lattice/Lattice_reduction.h b/lib/lattice/Lattice_reduction.h
index 3be4b6cb..41351c98 100644
--- a/lib/lattice/Lattice_reduction.h
+++ b/lib/lattice/Lattice_reduction.h
@@ -39,7 +39,7 @@ template<class vobj> inline RealD norm2(const Lattice<vobj> &arg){
 
 // Double inner product
 template<class vobj>
-inline ComplexD innerProduct(const Lattice<vobj> &left,const Lattice<vobj> &right) 
+inline ComplexD innerProduct(const Lattice<vobj> &left,const Lattice<vobj> &right)
 {
   typedef typename vobj::scalar_type scalar_type;
   typedef typename vobj::vector_typeD vector_type;
@@ -274,6 +274,115 @@ template<class vobj> inline void sliceSum(const Lattice<vobj> &Data,std::vector<
   }
 }
 
+template<class vobj>
+static void mySliceInnerProductVector( std::vector<ComplexD> & result, const Lattice<vobj> &lhs,const Lattice<vobj> &rhs,int orthogdim)
+{
+  // std::cout << GridLogMessage << "Start mySliceInnerProductVector" << std::endl;
+
+  typedef typename vobj::scalar_type scalar_type;
+  std::vector<scalar_type> lsSum;
+  localSliceInnerProductVector(result, lhs, rhs, lsSum, orthogdim);
+  globalSliceInnerProductVector(result, lhs, lsSum, orthogdim);
+  // std::cout << GridLogMessage << "End mySliceInnerProductVector" << std::endl;
+}
+
+template <class vobj>
+static void localSliceInnerProductVector(std::vector<ComplexD> &result, const Lattice<vobj> &lhs, const Lattice<vobj> &rhs, std::vector<typename vobj::scalar_type> &lsSum, int orthogdim)
+{
+  // std::cout << GridLogMessage << "Start prep" << std::endl;
+  typedef typename vobj::vector_type   vector_type;
+  typedef typename vobj::scalar_type   scalar_type;
+  GridBase  *grid = lhs._grid;
+  assert(grid!=NULL);
+  conformable(grid,rhs._grid);
+
+  const int    Nd = grid->_ndimension;
+  const int Nsimd = grid->Nsimd();
+
+  assert(orthogdim >= 0);
+  assert(orthogdim < Nd);
+
+  int fd=grid->_fdimensions[orthogdim];
+  int ld=grid->_ldimensions[orthogdim];
+  int rd=grid->_rdimensions[orthogdim];
+  // std::cout << GridLogMessage << "Start alloc" << std::endl;
+
+  std::vector<vector_type,alignedAllocator<vector_type> > lvSum(rd); // will locally sum vectors first
+  lsSum.resize(ld,scalar_type(0.0));                    // sum across these down to scalars
+  std::vector<iScalar<scalar_type>> extracted(Nsimd);   // splitting the SIMD  
+  // std::cout << GridLogMessage << "End alloc" << std::endl;
+
+  result.resize(fd); // And then global sum to return the same vector to every node for IO to file
+  for(int r=0;r<rd;r++){
+    lvSum[r]=zero;
+  }
+
+  int e1=    grid->_slice_nblock[orthogdim];
+  int e2=    grid->_slice_block [orthogdim];
+  int stride=grid->_slice_stride[orthogdim];
+  // std::cout << GridLogMessage << "End prep" << std::endl;
+  // std::cout << GridLogMessage << "Start parallel inner product, _rd = " << rd << std::endl;
+  vector_type vv;
+  parallel_for(int r=0;r<rd;r++)
+  {
+
+    int so=r*grid->_ostride[orthogdim]; // base offset for start of plane 
+
+    for(int n=0;n<e1;n++){
+      for(int b=0;b<e2;b++){
+        int ss = so + n * stride + b;
+        vv = TensorRemove(innerProduct(lhs._odata[ss], rhs._odata[ss]));
+        lvSum[r] = lvSum[r] + vv;
+      }
+    }
+  }
+  // std::cout << GridLogMessage << "End parallel inner product" << std::endl;
+
+  // Sum across simd lanes in the plane, breaking out orthog dir.
+  std::vector<int> icoor(Nd);
+  for(int rt=0;rt<rd;rt++){
+
+    iScalar<vector_type> temp; 
+    temp._internal = lvSum[rt];
+    extract(temp,extracted);
+
+    for(int idx=0;idx<Nsimd;idx++){
+
+      grid->iCoorFromIindex(icoor,idx);
+
+      int ldx =rt+icoor[orthogdim]*rd;
+
+      lsSum[ldx]=lsSum[ldx]+extracted[idx]._internal;
+
+    }
+  }
+  // std::cout << GridLogMessage << "End sum over simd lanes" << std::endl;
+}
+template <class vobj>
+static void globalSliceInnerProductVector(std::vector<ComplexD> &result, const Lattice<vobj> &lhs, std::vector<typename vobj::scalar_type> &lsSum, int orthogdim)
+{
+  typedef typename vobj::scalar_type scalar_type;
+  GridBase *grid = lhs._grid;
+  int fd = result.size();
+  int ld = lsSum.size();
+  // sum over nodes.
+  std::vector<scalar_type> gsum;
+  gsum.resize(fd, scalar_type(0.0));
+  // std::cout << GridLogMessage << "Start of gsum[t] creation:" << std::endl;
+  for(int t=0;t<fd;t++){
+    int pt = t/ld; // processor plane
+    int lt = t%ld;
+    if ( pt == grid->_processor_coor[orthogdim] ) {
+      gsum[t]=lsSum[lt];
+    }
+  }
+  // std::cout << GridLogMessage << "End of gsum[t] creation:" << std::endl;
+  // std::cout << GridLogMessage << "Start of GlobalSumVector:" << std::endl;
+  grid->GlobalSumVector(&gsum[0], fd);
+  // std::cout << GridLogMessage << "End of GlobalSumVector:" << std::endl;
+
+  result = gsum;
+}
 template<class vobj>
 static void sliceInnerProductVector( std::vector<ComplexD> & result, const Lattice<vobj> &lhs,const Lattice<vobj> &rhs,int orthogdim) 
 {
diff --git a/lib/qcd/action/fermion/CayleyFermion5D.cc b/lib/qcd/action/fermion/CayleyFermion5D.cc
index e3b77390..d4ed5a7c 100644
--- a/lib/qcd/action/fermion/CayleyFermion5D.cc
+++ b/lib/qcd/action/fermion/CayleyFermion5D.cc
@@ -67,6 +67,33 @@ void CayleyFermion5D<Impl>::ExportPhysicalFermionSolution(const FermionField &so
   axpby_ssp_pplus (tmp, 1., tmp       , 1., solution5d, 0, Ls-1);
   ExtractSlice(exported4d, tmp, 0, 0);
 }
+template<class Impl>  
+void CayleyFermion5D<Impl>::ExportPhysicalFermionSource(const FermionField &solution5d,FermionField &exported4d)
+{
+  int Ls = this->Ls;
+  FermionField tmp(this->FermionGrid());
+  tmp = solution5d;
+  conformable(solution5d._grid,this->FermionGrid());
+  conformable(exported4d._grid,this->GaugeGrid());
+  axpby_ssp_pplus (tmp, 0., solution5d, 1., solution5d, 0, 0);
+  axpby_ssp_pminus(tmp, 1., tmp       , 1., solution5d, 0, Ls-1);
+  ExtractSlice(exported4d, tmp, 0, 0);
+}
+template<class Impl>
+void CayleyFermion5D<Impl>::ImportUnphysicalFermion(const FermionField &input4d,FermionField &imported5d)
+{
+  int Ls = this->Ls;
+  FermionField tmp(this->FermionGrid());
+  conformable(imported5d._grid,this->FermionGrid());
+  conformable(input4d._grid   ,this->GaugeGrid());
+  tmp = zero;
+  InsertSlice(input4d, tmp, 0   , 0);
+  InsertSlice(input4d, tmp, Ls-1, 0);
+  axpby_ssp_pplus (tmp, 0., tmp, 1., tmp, 0, 0);
+  axpby_ssp_pminus(tmp, 0., tmp, 1., tmp, Ls-1, Ls-1);
+  imported5d=tmp;
+}
+
 template<class Impl>  
 void CayleyFermion5D<Impl>::ImportPhysicalFermionSource(const FermionField &input4d,FermionField &imported5d)
 {
diff --git a/lib/qcd/action/fermion/CayleyFermion5D.h b/lib/qcd/action/fermion/CayleyFermion5D.h
index b370b09d..3bf2a8b6 100644
--- a/lib/qcd/action/fermion/CayleyFermion5D.h
+++ b/lib/qcd/action/fermion/CayleyFermion5D.h
@@ -89,7 +89,9 @@ namespace Grid {
       virtual void   Dminus(const FermionField &psi, FermionField &chi);
       virtual void   DminusDag(const FermionField &psi, FermionField &chi);
       virtual void ExportPhysicalFermionSolution(const FermionField &solution5d,FermionField &exported4d);
+      virtual void ExportPhysicalFermionSource(const FermionField &solution5d, FermionField &exported4d);
       virtual void ImportPhysicalFermionSource(const FermionField &input4d,FermionField &imported5d);
+      virtual void ImportUnphysicalFermion(const FermionField &solution5d, FermionField &exported4d);
 
       /////////////////////////////////////////////////////
       // Instantiate different versions depending on Impl
diff --git a/lib/qcd/action/fermion/FermionOperator.h b/lib/qcd/action/fermion/FermionOperator.h
index 77fdbec5..dc1ab924 100644
--- a/lib/qcd/action/fermion/FermionOperator.h
+++ b/lib/qcd/action/fermion/FermionOperator.h
@@ -162,10 +162,18 @@ namespace Grid {
       {
 	imported = input;
       };
+      virtual void ImportUnphysicalFermion(const FermionField &input,FermionField &imported)
+      {
+	imported=input;
+      };
       virtual void ExportPhysicalFermionSolution(const FermionField &solution,FermionField &exported)
       {
 	exported=solution;
       };
+      virtual void ExportPhysicalFermionSource(const FermionField &solution,FermionField &exported)
+      {
+	exported=solution;
+      };
     };
 
   }