diff --git a/extras/Hadrons/AllToAllVectors.hpp b/extras/Hadrons/AllToAllVectors.hpp
index 7ddc0e81..de97051f 100644
--- a/extras/Hadrons/AllToAllVectors.hpp
+++ b/extras/Hadrons/AllToAllVectors.hpp
@@ -17,77 +17,53 @@ class A2AModesSchurDiagTwo
     A2AModesSchurDiagTwo(void) = default;
     virtual ~A2AModesSchurDiagTwo(void) = default;
 
-    void low_mode_v(Matrix &_Matrix, const Field &evec, const RealD &eval, Field &vout)
+    void Doo(Matrix &action, const Field &in, Field &out)
     {
+        Field tmp(in._grid);
 
-        GridBase *grid = _Matrix.RedBlackGrid();
+        action.MooeeInv(in, out);
+        action.Meooe(out, tmp);
+        action.MooeeInv(tmp, out);
+        action.Meooe(out, tmp);
 
-        Field src_o(grid);
-        Field sol_e(grid);
-        Field sol_o(grid);
-        Field tmp(grid);
-
-        pickCheckerboard(Odd, src_o, evec);
-        pickCheckerboard(Even, sol_e, vout);
-        pickCheckerboard(Odd, sol_o, vout);
-
-        /////////////////////////////////////////////////////
-        // v_ie = -(1/eval_i)* Mee Meo MooInv evec_i
-        /////////////////////////////////////////////////////
-        _Matrix.MooeeInv(src_o, tmp);
-        assert(tmp.checkerboard == Odd);
-        _Matrix.Meooe(tmp, sol_e);
-        assert(sol_e.checkerboard == Even);
-        _Matrix.Mooee(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
-        /////////////////////////////////////////////////////
-        _Matrix.MooeeInv(src_o, tmp);
-        assert(tmp.checkerboard == Odd);
-        sol_o = -(1.0 / eval) * tmp;
-        assert(sol_o.checkerboard == Odd);
-
-        setCheckerboard(vout, sol_e);
-        assert(sol_e.checkerboard == Even);
-        setCheckerboard(vout, sol_o);
-        assert(sol_o.checkerboard == Odd);
+        axpy(out, -1.0, tmp, in);
     }
 
-    void low_mode_w(Matrix &_Matrix, const Field &evec, const RealD &eval, Field &wout)
+    void low_mode_v(Matrix &action, const Field &evec, const RealD &eval, Field &vout, bool return_5d = true)
     {
-        GridBase *grid = _Matrix.RedBlackGrid();
-        SchurDiagTwoOperator<Matrix, Field> _HermOpEO(_Matrix);
 
+        GridBase *grid = action.RedBlackGrid();
         Field src_o(grid);
         Field sol_e(grid);
         Field sol_o(grid);
         Field tmp(grid);
 
-        GridBase *fgrid = _Matrix.Grid();
+        GridBase *fgrid = action.Grid();
         Field tmp_out(fgrid);
 
-        pickCheckerboard(Odd, src_o, evec);
-        pickCheckerboard(Even, sol_e, wout);
-        pickCheckerboard(Odd, sol_o, wout);
+        src_o = evec;
+        src_o.checkerboard = Odd;
+        pickCheckerboard(Even, sol_e, tmp_out);
+        pickCheckerboard(Odd, sol_o, tmp_out);
 
         /////////////////////////////////////////////////////
-        // w_ie = MeeInvDag MoeDag Doo evec_i
+        // v_ie = -(1/eval_i) * MeeInv Meo MooInv evec_i
         /////////////////////////////////////////////////////
-        _HermOpEO.Mpc(src_o, sol_e);
-        assert(sol_e.checkerboard == Odd);
-        _Matrix.MeooeDag(sol_e, tmp);
+        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);
-        _Matrix.MooeeInvDag(tmp, sol_e);
+        sol_e = (-1.0 / eval) * tmp;
         assert(sol_e.checkerboard == Even);
 
         /////////////////////////////////////////////////////
-        // w_io = Doo evec_i
+        // v_io = (1/eval_i) * MooInv evec_i
         /////////////////////////////////////////////////////
-        _HermOpEO.Mpc(src_o, sol_o);
+        action.MooeeInv(src_o, tmp);
+        assert(tmp.checkerboard == Odd);
+        sol_o = (1.0 / eval) * tmp;
         assert(sol_o.checkerboard == Odd);
 
         setCheckerboard(tmp_out, sol_e);
@@ -95,21 +71,85 @@ class A2AModesSchurDiagTwo
         setCheckerboard(tmp_out, sol_o);
         assert(sol_o.checkerboard == Odd);
 
-        _Matrix.Dminus(tmp_out, wout);
+        this->return_dim(action, tmp_out, vout, return_5d);
     }
 
-    void high_mode_v(Matrix &_Matrix, std::function<void(Field &, const Field &)> &Solver, const Field &source, Field &vout)
+    void low_mode_w(Matrix &action, const Field &evec, const RealD &eval, Field &wout, bool return_5d = true)
     {
-        GridBase *fgrid = _Matrix.Grid();
+        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_out(fgrid);
+        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
+        /////////////////////////////////////////////////////
+        Doo(action, 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
+        /////////////////////////////////////////////////////
+        Doo(action, 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, tmp_out);
+        this->return_dim(action, tmp_out, wout, return_5d);
+    }
+
+    void high_mode_v(Matrix &action, std::function<void(Field &, const Field &)> &Solver, const Field &source, Field &vout, bool return_5d = true)
+    {
+        GridBase *fgrid = action.Grid();
         Field tmp(fgrid);
+        Field tmp_out(fgrid);
 
-        _Matrix.Dminus(source, tmp);
-        Solver(vout, tmp);
+        action.Dminus(source, tmp);
+        Solver(tmp_out, source); // Note: Solver is Solver(out, in)
+        this->return_dim(action, tmp_out, vout, return_5d);
     }
 
-    void high_mode_w(Matrix &_Matrix, const Field &source, Field &wout)
+    void high_mode_w(Matrix &action, const Field &source4d, Field &wout, bool return_5d = true)
     {
-        wout = source;
+        // GridBase *fgrid = action.Grid();
+        // Field tmp_out(fgrid);
+
+        // tmp_out = source;
+        // this->return_dim(action, tmp_out, wout, return_5d);
+        wout = source4d;
+    }
+
+    void return_dim(Matrix &action, const Field &in, Field &out, bool return_5d)
+    {
+        if (return_5d)
+        {
+            out = in;
+        }
+        else
+        {
+            action.ExportPhysicalFermionSolution(in, out);
+        }
     }
 };
 
@@ -123,7 +163,7 @@ class A2ALMSchurDiagTwo : public A2AModesSchurDiagTwo<Field, Matrix>
   private:
     const std::vector<Field> &evec;
     const std::vector<RealD> &eval;
-    Matrix                   &action;
+    Matrix &action;
 
   public:
     A2ALMSchurDiagTwo(const std::vector<Field> &_evec, const std::vector<RealD> &_eval, Matrix &_action) : evec(_evec), eval(_eval), action(_action){};
@@ -138,10 +178,10 @@ template <class FineField, class CoarseField, class Matrix>
 class A2ALMSchurDiagTwoCoarse : public A2AModesSchurDiagTwo<FineField, Matrix>
 {
   private:
-    const std::vector<FineField>   &subspace;
+    const std::vector<FineField> &subspace;
     const std::vector<CoarseField> &evec_coarse;
-    const std::vector<RealD>       &eval_coarse;
-    Matrix                         &action;
+    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)
@@ -164,58 +204,10 @@ template <class Field, class Matrix>
 class A2AHMSchurDiagTwo : virtual public A2AModesSchurDiagTwo<Field, Matrix>
 {
   public:
-    void operator()(Matrix &_Matrix, std::function<void(Field &, const Field &)> &Solver, const Field &source, Field &vout, Field &wout)
+    void operator()(Matrix &action, std::function<void(Field &, const Field &)> &Solver, const Field &source, Field &vout, Field &wout)
     {
-        this->high_mode_v(_Matrix, Solver, source, vout);
-        this->high_mode_w(_Matrix, source, wout);
-    }
-};
-
-template <class Field, class Matrix>
-class A2AVectorsReturnHigh : public A2AModesSchurDiagTwo<Field, Matrix>
-{
-  private:
-    Matrix &action;
-    std::function<void(Field &, const Field &)> &Solver;
-    const int Nh, Ls;
-
-  public:
-    std::vector<Field> w_high, v_high;
-    A2AVectorsReturnHigh(Matrix &_action,
-                         std::function<void(Field &, const Field &)> &_Solver,
-                         const int _Nh, const int _Ls)
-        : action(_action),
-          Solver(_Solver),
-          Nh(_Nh), Ls(_Ls)
-    {
-        GridBase *fgrid = action.Grid();
-        resize(Nh, fgrid);
-    };
-
-    void resize(const size_t size, GridBase *grid)
-    {
-        w_high.resize(size, grid);
-        v_high.resize(size, grid);
-    }
-
-    void high_modes(Field &source, int i)
-    {
-        this->high_mode_v(action, Solver, source, v_high[i]);
-        this->high_mode_w(action, source, w_high[i]);
-    }
-
-    void operator()(int i, Field &vout, Field &wout)
-    {
-        if (Ls > 1)
-        {
-            vout = v_high[i];
-            wout = w_high[i];
-        }
-        else
-        {
-            action.ExportPhysicalFermionSolution(v_high[i], vout);
-            action.ExportPhysicalFermionSolution(w_high[i], wout);
-        }
+        this->high_mode_v(action, Solver, source, vout);
+        this->high_mode_w(action, source, wout);
     }
 };
 
@@ -227,26 +219,36 @@ template <class Field, class Matrix>
 class A2AVectorsReturn : public A2AModesSchurDiagTwo<Field, Matrix>
 {
   private:
-    const std::vector<Field> &evec;
-    const std::vector<RealD> &eval;
+    const std::vector<Field> *evec;
+    const std::vector<RealD> *eval;
     Matrix &action;
     std::function<void(Field &, const Field &)> &Solver;
-    const int Nl, Nh, Ls;
-
-  public:
+    const int Nl, Nh;
+    const bool return_5d;
     std::vector<Field> w_high, v_high;
 
-    A2AVectorsReturn(const std::vector<Field> &_evec, const std::vector<RealD> &_eval,
+  public:
+    A2AVectorsReturn(const std::vector<Field> *_evec, const std::vector<RealD> *_eval,
                      Matrix &_action,
                      std::function<void(Field &, const Field &)> &_Solver,
-                     const int _Nl, const int _Nh, const int _Ls)
+                     const int _Nl, const int _Nh,
+                     const bool _return_5d)
         : evec(_evec), eval(_eval),
           action(_action),
           Solver(_Solver),
-          Nl(_Nl), Nh(_Nh), Ls(_Ls)
+          Nl(_Nl), Nh(_Nh),
+          return_5d(_return_5d)
     {
-        GridBase *fgrid = action.Grid();
-        resize(Nh, fgrid);
+        GridBase *grid;
+        if (return_5d)
+        {
+            grid = action.Grid();
+        }
+        else
+        {
+            grid = action.GaugeGrid();
+        }
+        resize(Nh, grid);
     };
 
     void resize(const size_t size, GridBase *grid)
@@ -255,38 +257,25 @@ class A2AVectorsReturn : public A2AModesSchurDiagTwo<Field, Matrix>
         v_high.resize(size, grid);
     }
 
-    void high_modes(Field &source, int i)
+    void high_modes(Field &source5d, Field &source4d, int i)
     {
-        this->high_mode_v(action, Solver, source, v_high[i]);
-        this->high_mode_w(action, source, w_high[i]);
+        LOG(Message) << "A2A high modes for i = " << i << std::endl;
+        this->high_mode_v(action, Solver, source5d, v_high[i], return_5d);
+        this->high_mode_w(action, source4d, w_high[i], return_5d);
     }
 
     void operator()(int i, Field &vout, Field &wout)
     {
-
-        GridBase *fgrid = action.Grid();
-        Field vtmp(fgrid);
-        Field wtmp(fgrid);
         if (i < Nl)
-        {   
-            this->low_mode_v(action, evec[i], eval[i], vtmp);
-            this->low_mode_w(action, evec[i], eval[i], wtmp);
+        {
+            LOG(Message) << "A2A low modes for i = " << i << std::endl;
+            this->low_mode_v(action, evec->at(i), eval->at(i), vout, return_5d);
+            this->low_mode_w(action, evec->at(i), eval->at(i), wout, return_5d);
         }
         else
         {
-            vtmp = v_high[i-Nl];
-            wtmp = w_high[i-Nl];
-        }
-
-        if (Ls > 1)
-        {
-            vout = vtmp;
-            wout = wtmp;
-        }
-        else
-        {
-            action.ExportPhysicalFermionSolution(vtmp, vout);
-            action.ExportPhysicalFermionSolution(wtmp, wout);
+            vout = v_high[i - Nl];
+            wout = w_high[i - Nl];
         }
     }
 };
diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp
index e46bf01b..f81c01da 100644
--- a/extras/Hadrons/Modules.hpp
+++ b/extras/Hadrons/Modules.hpp
@@ -1,52 +1,53 @@
-#include <Grid/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp>
-#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp>
-#include <Grid/Hadrons/Modules/MContraction/DiscLoop.hpp>
-#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp>
-#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
-#include <Grid/Hadrons/Modules/MContraction/WardIdentity.hpp>
-#include <Grid/Hadrons/Modules/MContraction/Gamma3pt.hpp>
-#include <Grid/Hadrons/Modules/MContraction/Baryon.hpp>
-#include <Grid/Hadrons/Modules/MContraction/Meson.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/MGauge/Unit.hpp>
-#include <Grid/Hadrons/Modules/MSource/Wall.hpp>
-#include <Grid/Hadrons/Modules/MSource/Point.hpp>
-#include <Grid/Hadrons/Modules/MSource/Z2.hpp>
 #include <Grid/Hadrons/Modules/MSource/SeqConserved.hpp>
 #include <Grid/Hadrons/Modules/MSource/SeqGamma.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/MScalarSUN/TwoPointNPR.hpp>
-#include <Grid/Hadrons/Modules/MScalarSUN/TrMag.hpp>
-#include <Grid/Hadrons/Modules/MScalarSUN/Grad.hpp>
-#include <Grid/Hadrons/Modules/MScalarSUN/EMT.hpp>
-#include <Grid/Hadrons/Modules/MScalarSUN/TwoPoint.hpp>
-#include <Grid/Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp>
-#include <Grid/Hadrons/Modules/MScalarSUN/TrKinetic.hpp>
-#include <Grid/Hadrons/Modules/MScalarSUN/Div.hpp>
-#include <Grid/Hadrons/Modules/MScalarSUN/ShiftProbe.hpp>
-#include <Grid/Hadrons/Modules/MScalarSUN/StochFreeField.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/MScalar/ChargedProp.hpp>
-#include <Grid/Hadrons/Modules/MScalar/Scalar.hpp>
-#include <Grid/Hadrons/Modules/MScalar/FreeProp.hpp>
-#include <Grid/Hadrons/Modules/MIO/LoadNersc.hpp>
-#include <Grid/Hadrons/Modules/MIO/LoadEigenPack.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/MSink/Smear.hpp>
+#include <Grid/Hadrons/Modules/MSink/Point.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/MFermion/GaugeProp.hpp>
-#include <Grid/Hadrons/Modules/MFermion/FreeProp.hpp>
-#include <Grid/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp>
-#include <Grid/Hadrons/Modules/MSolver/RBPrecCG.hpp>
-#include <Grid/Hadrons/Modules/MSolver/A2AVectors.hpp>
-#include <Grid/Hadrons/Modules/MAction/Wilson.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/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/MAction/ZMobiusDWF.hpp>
+#include <Grid/Hadrons/Modules/MAction/Wilson.hpp>
 #include <Grid/Hadrons/Modules/MAction/WilsonClover.hpp>
 #include <Grid/Hadrons/Modules/MAction/DWF.hpp>
-#include <Grid/Hadrons/Modules/MSink/Point.hpp>
-#include <Grid/Hadrons/Modules/MSink/Smear.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/Baryon.hpp>
+#include <Grid/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp>
+#include <Grid/Hadrons/Modules/MContraction/A2AMeson.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/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/StochEm.hpp>
+#include <Grid/Hadrons/Modules/MGauge/FundtoHirep.hpp>
+#include <Grid/Hadrons/Modules/MGauge/Unit.hpp>
+#include <Grid/Hadrons/Modules/MGauge/Random.hpp>
diff --git a/extras/Hadrons/Modules/MSolver/A2AVectors.hpp b/extras/Hadrons/Modules/MSolver/A2AVectors.hpp
index e4690b13..107d1edb 100644
--- a/extras/Hadrons/Modules/MSolver/A2AVectors.hpp
+++ b/extras/Hadrons/Modules/MSolver/A2AVectors.hpp
@@ -18,7 +18,7 @@ class A2AVectorsPar: Serializable
 {
 public:
   GRID_SERIALIZABLE_CLASS_MEMBERS(A2AVectorsPar,
-                                  int, Ls,
+                                  bool, return_5d,
                                   int, Nl,
                                   int, N,
                                   std::vector<std::string>, sources,
@@ -39,7 +39,6 @@ class TA2AVectors : public Module<A2AVectorsPar>
 
     typedef A2AModesSchurDiagTwo<typename FImpl::FermionField, FMat> A2ABase;
     typedef A2AVectorsReturn<typename FImpl::FermionField, FMat> A2AReturn;
-    typedef A2AVectorsReturnHigh<typename FImpl::FermionField, FMat> A2AReturnHigh;
 
   public:
     // constructor
@@ -56,7 +55,7 @@ class TA2AVectors : public Module<A2AVectorsPar>
 
   private:
     unsigned int Ls_;
-    std::string retName_;
+    std::string retName_, whighName_, vhighName_;
 };
 
 MODULE_REGISTER_TMP(A2AVectors, ARG(TA2AVectors<FIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>), MSolver);
@@ -70,13 +69,15 @@ template <typename FImpl, int nBasis>
 TA2AVectors<FImpl, nBasis>::TA2AVectors(const std::string name)
 : Module<A2AVectorsPar>(name)
 , retName_ (name + "_ret")
+, whighName_ (name + "_whigh")
+, vhighName_ (name + "_vhigh")
 {}
 
 // dependencies/products ///////////////////////////////////////////////////////
 template <typename FImpl, int nBasis>
 std::vector<std::string> TA2AVectors<FImpl, nBasis>::getInput(void)
 {
-    std::vector<std::string> in = {par().action, par().solver + "_subtract"};
+    std::vector<std::string> in = {par().action, par().solver, par().solver + "_subtract"};
 
     int n = par().sources.size();
 
@@ -103,39 +104,57 @@ void TA2AVectors<FImpl, nBasis>::setup(void)
     int N = par().N;
     int Nl = par().Nl;
     int Nh = N - Nl;
-    int Ls = par().Ls;
+    bool return_5d = par().return_5d;
+    int Ls, Ls_;
+
     Ls_ = env().getObjectLs(par().solver + "_subtract");
+    auto &solver = envGet(SolverFn, par().solver + "_subtract");
+    if (!(Nl > 0))
+    {
+        Ls_ = env().getObjectLs(par().solver);
+        auto &solver = envGet(SolverFn, par().solver);
+    }
+
+    if (return_5d)
+    {
+        Ls = Ls_;
+    }
+    else
+    {
+        Ls = 1;
+    }
 
     auto &action = envGet(FMat, par().action);
-    auto &solver = envGet(SolverFn, par().solver + "_subtract");
-
-    GridBase *fgrid = action.Grid();
 
     envTmpLat(FermionField, "ferm_src", Ls_);
     envTmpLat(FermionField, "tmp");
 
+    std::vector<FermionField> *evec;
+    const std::vector<RealD> *eval;
+
     if (Nl > 0)
     {
         // Low modes
         auto &epack = envGet(EPack, par().eigenpack);
 
-        LOG(Message) << "using a2a with eigenpack '"
-                     << par().eigenpack << "' ("
-                     << epack.evec.size() << " modes)" << std::endl;
-
-        envCreateDerived(A2ABase, A2AReturn, retName_, Ls,
-                         epack.evec, epack.eval,
-                         action,
-                         solver,
-                         Nl, Nh, Ls);
+        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) << "using a2a with high modes only" << std::endl;
-        envCreateDerived(A2ABase, A2AReturnHigh, retName_, Ls, action,
-                         solver,
-                         Nh, Ls);
+        LOG(Message) << "Creating a2a vectors " << getName() <<
+                     " using " << Nh << " high modes only." << std::endl;
     }
+    envCreateDerived(A2ABase, A2AReturn, retName_, Ls,
+                     evec, eval,
+                     action,
+                     solver,
+                     Nl, Nh,
+                     return_5d);
 }
 
 // execution ///////////////////////////////////////////////////////////////////
@@ -145,30 +164,41 @@ void TA2AVectors<FImpl, nBasis>::execute(void)
     auto &action = envGet(FMat, par().action);
 
     int Nt = env().getDim(Tp);
-    int Nl = par().Nl; // Number of low modes
     int Nc = FImpl::Dimension;
+    int Ls_;
+    int Nl = par().Nl;
+    Ls_ = env().getObjectLs(par().solver + "_subtract");
+    if (!(Nl > 0))
+    {
+        Ls_ = env().getObjectLs(par().solver);
+    }
+
+    auto &a2areturn = envGetDerived(A2ABase, A2AReturn, retName_);
 
     // High modes
     auto sources = par().sources;
+    int Nsrc = par().sources.size();
 
     envGetTmp(FermionField, ferm_src);
     envGetTmp(FermionField, tmp);
+
     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 < Nt; T++)
+        for (unsigned int T = 0; T < Nsrc; T++)
         {
-            auto &fullSrc = envGet(PropagatorField, sources[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, fullSrc, s, c);
+                    PropToFerm<FImpl>(ferm_src, prop_src, s, c);
                 }
                 else
                 {
-                    PropToFerm<FImpl>(tmp, fullSrc, s, c);
+                    PropToFerm<FImpl>(tmp, prop_src, s, c);
                     action.ImportPhysicalFermionSource(tmp, ferm_src);
                 }
             }
@@ -181,23 +211,13 @@ void TA2AVectors<FImpl, nBasis>::execute(void)
                 }
                 else
                 {
-                    PropToFerm<FImpl>(ferm_src, fullSrc, s, c);
+                    PropToFerm<FImpl>(ferm_src, prop_src, s, c);
                 }
             }
-
-            if (Nl > 0)
-            {
-                auto a2areturn = envGetDerived(A2ABase, A2AReturn, retName_);
-                a2areturn.high_modes(ferm_src, N_count);
-            }
-            else
-            {
-                auto a2areturnhigh = envGetDerived(A2ABase, A2AReturnHigh, retName_);
-                a2areturnhigh.high_modes(ferm_src, N_count);
-            }
+            LOG(Message) << "a2areturn.high_modes Ncount = " << N_count << std::endl;
+            a2areturn.high_modes(ferm_src, tmp, N_count);
             N_count++;
         }
-    
 }
 END_MODULE_NAMESPACE
 
diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc
index aecfdd5a..9bdb48cc 100644
--- a/extras/Hadrons/modules.inc
+++ b/extras/Hadrons/modules.inc
@@ -1,105 +1,107 @@
 modules_cc =\
-  Modules/MContraction/WeakHamiltonianEye.cc \
-  Modules/MContraction/Baryon.cc \
-  Modules/MContraction/WeakNeutral4ptDisc.cc \
-  Modules/MContraction/DiscLoop.cc \
-  Modules/MContraction/WeakHamiltonianNonEye.cc \
-  Modules/MContraction/Meson.cc \
-  Modules/MContraction/WardIdentity.cc \
-  Modules/MContraction/Gamma3pt.cc \
-  Modules/MGauge/Unit.cc \
-  Modules/MGauge/FundtoHirep.cc \
-  Modules/MGauge/StochEm.cc \
-  Modules/MGauge/Random.cc \
-  Modules/MSource/Wall.cc \
   Modules/MSource/SeqConserved.cc \
-  Modules/MSource/Z2.cc \
   Modules/MSource/SeqGamma.cc \
   Modules/MSource/Point.cc \
-  Modules/MUtilities/TestSeqConserved.cc \
-  Modules/MUtilities/TestSeqGamma.cc \
-  Modules/MLoop/NoiseLoop.cc \
-  Modules/MScalarSUN/TrPhi.cc \
-  Modules/MScalarSUN/TrMag.cc \
-  Modules/MScalarSUN/Grad.cc \
-  Modules/MScalarSUN/TransProj.cc \
-  Modules/MScalarSUN/StochFreeField.cc \
-  Modules/MScalarSUN/TwoPoint.cc \
-  Modules/MScalarSUN/TwoPointNPR.cc \
-  Modules/MScalarSUN/EMT.cc \
-  Modules/MScalarSUN/ShiftProbe.cc \
-  Modules/MScalarSUN/TrKinetic.cc \
-  Modules/MScalarSUN/Div.cc \
-  Modules/MScalarSUN/TimeMomProbe.cc \
-  Modules/MScalar/FreeProp.cc \
-  Modules/MScalar/ChargedProp.cc \
-  Modules/MIO/LoadCoarseEigenPack.cc \
+  Modules/MSource/Z2.cc \
+  Modules/MSource/Wall.cc \
+  Modules/MSink/Point.cc \
+  Modules/MSink/Smear.cc \
   Modules/MIO/LoadNersc.cc \
   Modules/MIO/LoadEigenPack.cc \
+  Modules/MIO/LoadCoarseEigenPack.cc \
   Modules/MIO/LoadBinary.cc \
-  Modules/MFermion/FreeProp.cc \
-  Modules/MFermion/GaugeProp.cc \
-  Modules/MSolver/LocalCoherenceLanczos.cc \
-  Modules/MSolver/A2AVectors.cc \
-  Modules/MSolver/RBPrecCG.cc \
-  Modules/MAction/WilsonClover.cc \
-  Modules/MAction/Wilson.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/ZMobiusDWF.cc \
-  Modules/MSink/Point.cc \
-  Modules/MSink/Smear.cc
+  Modules/MAction/Wilson.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/Gamma3pt.cc \
+  Modules/MContraction/WeakNeutral4ptDisc.cc \
+  Modules/MContraction/Meson.cc \
+  Modules/MScalar/ChargedProp.cc \
+  Modules/MScalar/FreeProp.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/StochEm.cc \
+  Modules/MGauge/FundtoHirep.cc
 
 modules_hpp =\
-  Modules/MContraction/WeakNeutral4ptDisc.hpp \
-  Modules/MContraction/WeakHamiltonianEye.hpp \
-  Modules/MContraction/DiscLoop.hpp \
-  Modules/MContraction/WeakHamiltonianNonEye.hpp \
-  Modules/MContraction/WeakHamiltonian.hpp \
-  Modules/MContraction/WardIdentity.hpp \
-  Modules/MContraction/Gamma3pt.hpp \
-  Modules/MContraction/Baryon.hpp \
-  Modules/MContraction/Meson.hpp \
-  Modules/MGauge/Random.hpp \
-  Modules/MGauge/FundtoHirep.hpp \
-  Modules/MGauge/StochEm.hpp \
-  Modules/MGauge/Unit.hpp \
-  Modules/MSource/Wall.hpp \
-  Modules/MSource/Point.hpp \
-  Modules/MSource/Z2.hpp \
   Modules/MSource/SeqConserved.hpp \
   Modules/MSource/SeqGamma.hpp \
-  Modules/MUtilities/TestSeqGamma.hpp \
-  Modules/MUtilities/TestSeqConserved.hpp \
-  Modules/MLoop/NoiseLoop.hpp \
-  Modules/MScalarSUN/TwoPointNPR.hpp \
-  Modules/MScalarSUN/TrMag.hpp \
-  Modules/MScalarSUN/Grad.hpp \
-  Modules/MScalarSUN/EMT.hpp \
-  Modules/MScalarSUN/TwoPoint.hpp \
-  Modules/MScalarSUN/TimeMomProbe.hpp \
-  Modules/MScalarSUN/TrKinetic.hpp \
-  Modules/MScalarSUN/Div.hpp \
-  Modules/MScalarSUN/ShiftProbe.hpp \
-  Modules/MScalarSUN/StochFreeField.hpp \
-  Modules/MScalarSUN/TrPhi.hpp \
-  Modules/MScalarSUN/Utils.hpp \
-  Modules/MScalarSUN/TransProj.hpp \
-  Modules/MScalar/ChargedProp.hpp \
-  Modules/MScalar/Scalar.hpp \
-  Modules/MScalar/FreeProp.hpp \
-  Modules/MIO/LoadNersc.hpp \
-  Modules/MIO/LoadEigenPack.hpp \
+  Modules/MSource/Z2.hpp \
+  Modules/MSource/Point.hpp \
+  Modules/MSource/Wall.hpp \
+  Modules/MSink/Smear.hpp \
+  Modules/MSink/Point.hpp \
   Modules/MIO/LoadBinary.hpp \
+  Modules/MIO/LoadEigenPack.hpp \
   Modules/MIO/LoadCoarseEigenPack.hpp \
-  Modules/MFermion/GaugeProp.hpp \
-  Modules/MFermion/FreeProp.hpp \
-  Modules/MSolver/LocalCoherenceLanczos.hpp \
-  Modules/MSolver/RBPrecCG.hpp \
-  Modules/MSolver/A2AVectors.hpp \
-  Modules/MAction/Wilson.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/TimeMomProbe.hpp \
+  Modules/MScalarSUN/Div.hpp \
+  Modules/MScalarSUN/TrMag.hpp \
+  Modules/MScalarSUN/EMT.hpp \
   Modules/MAction/ZMobiusDWF.hpp \
+  Modules/MAction/Wilson.hpp \
   Modules/MAction/WilsonClover.hpp \
   Modules/MAction/DWF.hpp \
-  Modules/MSink/Point.hpp \
-  Modules/MSink/Smear.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/Baryon.hpp \
+  Modules/MContraction/WeakNeutral4ptDisc.hpp \
+  Modules/MContraction/A2AMeson.hpp \
+  Modules/MScalar/Scalar.hpp \
+  Modules/MScalar/FreeProp.hpp \
+  Modules/MScalar/ChargedProp.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/StochEm.hpp \
+  Modules/MGauge/FundtoHirep.hpp \
+  Modules/MGauge/Unit.hpp \
+  Modules/MGauge/Random.hpp