diff --git a/extras/Hadrons/A2AVectors.hpp b/extras/Hadrons/A2AVectors.hpp
index 4cbc7f2a..6e865203 100644
--- a/extras/Hadrons/A2AVectors.hpp
+++ b/extras/Hadrons/A2AVectors.hpp
@@ -219,4 +219,4 @@ void A2AVectorsSchurDiagTwo<FImpl>::makeHighModeW5D(FermionField &wout_4d,
 
 END_HADRONS_NAMESPACE
 
-#endif // A2A_Vectors_hpp_
\ No newline at end of file
+#endif // A2A_Vectors_hpp_
diff --git a/extras/Hadrons/DilutedNoise.hpp b/extras/Hadrons/DilutedNoise.hpp
index 16561c99..250f8445 100644
--- a/extras/Hadrons/DilutedNoise.hpp
+++ b/extras/Hadrons/DilutedNoise.hpp
@@ -188,4 +188,4 @@ void TimeDilutedSpinColorDiagonalNoise<FImpl>::generateNoise(GridParallelRNG &rn
 
 END_HADRONS_NAMESPACE
 
-#endif // Hadrons_DilutedNoise_hpp_
\ No newline at end of file
+#endif // Hadrons_DilutedNoise_hpp_
diff --git a/extras/Hadrons/Global.cc b/extras/Hadrons/Global.cc
index a173bd8b..1e747cda 100644
--- a/extras/Hadrons/Global.cc
+++ b/extras/Hadrons/Global.cc
@@ -192,7 +192,7 @@ void Hadrons::printTimeProfile(const std::map<std::string, GridTime> &timing,
     }
     for (auto &rt: rtiming)
     {
-        LOG(Message) << std::right << std::setw(width) << rt.second << ": " 
+        LOG(Message) << std::setw(width) << rt.second << ": " 
                      << rt.first << " us (" << std::fixed 
                      << std::setprecision(1) 
                      << static_cast<double>(rt.first)/dtotal*100 << "%)"
diff --git a/extras/Hadrons/Module.cc b/extras/Hadrons/Module.cc
index f0f3def1..8e915f6b 100644
--- a/extras/Hadrons/Module.cc
+++ b/extras/Hadrons/Module.cc
@@ -78,12 +78,21 @@ void ModuleBase::startTimer(const std::string &name)
 GridTime ModuleBase::getTimer(const std::string &name)
 {
     GridTime t;
-
+    
     if (!name.empty())
     {
-        stopTimer(name);
-        t = timer_.at(name).Elapsed();
-        startTimer(name);
+        try
+        {
+            bool running = timer_.at(name).isRunning();
+
+            if (running) stopTimer(name);
+            t = timer_.at(name).Elapsed();
+            if (running) startTimer(name);
+        }
+        catch (std::out_of_range &)
+        {
+            t = GridTime::zero();
+        }
     }
     else
     {
@@ -93,6 +102,11 @@ GridTime ModuleBase::getTimer(const std::string &name)
     return t;
 }
 
+double ModuleBase::getDTimer(const std::string &name)
+{
+    return static_cast<double>(getTimer(name).count());
+}
+
 void ModuleBase::startCurrentTimer(const std::string &name)
 {
     if (!name.empty())
diff --git a/extras/Hadrons/Module.hpp b/extras/Hadrons/Module.hpp
index 8f260c19..02a20e5f 100644
--- a/extras/Hadrons/Module.hpp
+++ b/extras/Hadrons/Module.hpp
@@ -129,8 +129,6 @@ if (env().getGrid()->IsBoss() and !ioStem.empty())\
     }\
 }
 
-#define MODULE_TIMER(name) this->startCurrentTimer(name)
-
 /******************************************************************************
  *                            Module class                                    *
  ******************************************************************************/
@@ -166,6 +164,7 @@ public:
     // timers
     void                            startTimer(const std::string &name);
     GridTime                        getTimer(const std::string &name);
+    double                          getDTimer(const std::string &name);
     void                            startCurrentTimer(const std::string &name);
     void                            stopTimer(const std::string &name);
     void                            stopCurrentTimer(void);
diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp
index 5cce335a..86d50755 100644
--- a/extras/Hadrons/Modules.hpp
+++ b/extras/Hadrons/Modules.hpp
@@ -1,5 +1,40 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules.hpp
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+Author: Guido Cossu <guido.cossu@ed.ac.uk>
+Author: James Harrison <jch1g10@soton.ac.uk>
+Author: Lanny91 <andrew.lawson@gmail.com>
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: Vera Guelpers <vmg1n14@soton.ac.uk>
+Author: fionnoh <fionnoh@gmail.com>
+Author: pretidav <david.preti@csic.es>
+
+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/Hadrons/Modules/MContraction/Baryon.hpp>
 #include <Grid/Hadrons/Modules/MContraction/A2AMesonField.hpp>
+#include <Grid/Hadrons/Modules/MContraction/A2AMesonFieldKernels.hpp>
 #include <Grid/Hadrons/Modules/MContraction/Meson.hpp>
 #include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
 #include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp>
diff --git a/extras/Hadrons/Modules/MContraction/A2AMesonField.cc b/extras/Hadrons/Modules/MContraction/A2AMesonField.cc
index ffc45e22..3d85acb5 100644
--- a/extras/Hadrons/Modules/MContraction/A2AMesonField.cc
+++ b/extras/Hadrons/Modules/MContraction/A2AMesonField.cc
@@ -6,6 +6,7 @@ Source file: extras/Hadrons/Modules/MContraction/A2AMesonField.cc
 
 Copyright (C) 2015-2018
 
+Author: Antonin Portelli <antonin.portelli@me.com>
 Author: paboyle <paboyle@ph.ed.ac.uk>
 
 This program is free software; you can redistribute it and/or modify
diff --git a/extras/Hadrons/Modules/MContraction/A2AMesonField.hpp b/extras/Hadrons/Modules/MContraction/A2AMesonField.hpp
index e0818874..568985bc 100644
--- a/extras/Hadrons/Modules/MContraction/A2AMesonField.hpp
+++ b/extras/Hadrons/Modules/MContraction/A2AMesonField.hpp
@@ -34,8 +34,8 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #include <Grid/Hadrons/Module.hpp>
 #include <Grid/Hadrons/ModuleFactory.hpp>
 #include <Grid/Hadrons/A2AVectors.hpp>
-#include <Grid/Eigen/unsupported/CXX11/Tensor>
 #include <Grid/Hadrons/Modules/MSolver/A2AVectors.hpp>
+#include <Grid/Hadrons/Modules/MContraction/A2AMesonFieldKernels.hpp>
 
 BEGIN_HADRONS_NAMESPACE
 
@@ -44,9 +44,6 @@ BEGIN_HADRONS_NAMESPACE
  ******************************************************************************/
 BEGIN_MODULE_NAMESPACE(MContraction)
 
-typedef std::pair<Gamma::Algebra, Gamma::Algebra> GammaPair;
-
-
 class A2AMesonFieldPar: Serializable
 {
   public:
@@ -66,7 +63,6 @@ class TA2AMesonField : public Module<A2AMesonFieldPar>
 public:
     FERM_TYPE_ALIASES(FImpl,);
     SOLVER_TYPE_ALIASES(FImpl,);
-    typedef Eigen::TensorMap<Eigen::Tensor<Complex, 5, Eigen::RowMajor>> MesonField;
 public:
     // constructor
     TA2AMesonField(const std::string name);
@@ -144,24 +140,24 @@ void TA2AMesonField<FImpl>::setup(void)
     mom_.clear();
     if (par().gammas == "all")
     {
-    gamma_ = {
-        Gamma::Algebra::Gamma5,
-        Gamma::Algebra::Identity,    
-        Gamma::Algebra::GammaX,
-        Gamma::Algebra::GammaY,
-        Gamma::Algebra::GammaZ,
-        Gamma::Algebra::GammaT,
-        Gamma::Algebra::GammaXGamma5,
-        Gamma::Algebra::GammaYGamma5,
-        Gamma::Algebra::GammaZGamma5,
-        Gamma::Algebra::GammaTGamma5,
-        Gamma::Algebra::SigmaXY,
-        Gamma::Algebra::SigmaXZ,
-        Gamma::Algebra::SigmaXT,
-        Gamma::Algebra::SigmaYZ,
-        Gamma::Algebra::SigmaYT,
-        Gamma::Algebra::SigmaZT
-    };
+        gamma_ = {
+            Gamma::Algebra::Gamma5,
+            Gamma::Algebra::Identity,    
+            Gamma::Algebra::GammaX,
+            Gamma::Algebra::GammaY,
+            Gamma::Algebra::GammaZ,
+            Gamma::Algebra::GammaT,
+            Gamma::Algebra::GammaXGamma5,
+            Gamma::Algebra::GammaYGamma5,
+            Gamma::Algebra::GammaZGamma5,
+            Gamma::Algebra::GammaTGamma5,
+            Gamma::Algebra::SigmaXY,
+            Gamma::Algebra::SigmaXZ,
+            Gamma::Algebra::SigmaXT,
+            Gamma::Algebra::SigmaYZ,
+            Gamma::Algebra::SigmaYT,
+            Gamma::Algebra::SigmaZT
+        };
     }
     else
     {
@@ -227,7 +223,7 @@ void TA2AMesonField<FImpl>::execute(void)
 
     if (!hasPhase_)
     {
-        MODULE_TIMER("Momentum phases");
+        startTimer("Momentum phases");
         for (unsigned int j = 0; j < nmom; ++j)
         {
             Complex           i(0.0,1.0);
@@ -243,9 +239,9 @@ void TA2AMesonField<FImpl>::execute(void)
             ph[j] = exp((Real)(2*M_PI)*i*ph[j]);
         }
         hasPhase_ = true;
+        stopTimer("Momentum phases");
     }
     
-
     //////////////////////////////////////////////////////////////////////////
     // i,j   is first  loop over SchurBlock factors reusing 5D matrices
     // ii,jj is second loop over cacheBlock factors for high perf contractoin
@@ -256,12 +252,8 @@ void TA2AMesonField<FImpl>::execute(void)
     double flops = 0.0;
     double bytes = 0.0;
     double vol   = env().getVolume();
-    double t_schur=0;
     double t_contr=0;
-    double t_int_0=0;
-    double t_int_1=0;
-    double t_int_2=0;
-    double t_int_3=0;
+
 
     envGetTmp(Vector<Complex>, mfBuf);
     envGetTmp(Vector<Complex>, mfCache);
@@ -277,9 +269,6 @@ void TA2AMesonField<FImpl>::execute(void)
         int N_ii = MIN(N_i-i,block);
         int N_jj = MIN(N_j-j,block);
 
-        t_schur-=usecond();
-        t_schur+=usecond();
-
         LOG(Message) << "Meson field block " 
                     << j/block + NBlock_j*i/block + 1 
                     << "/" << NBlock_i*NBlock_j << " [" << i <<" .. " 
@@ -296,17 +285,17 @@ void TA2AMesonField<FImpl>::execute(void)
             int N_jjj = MIN(N_jj-jj,cacheBlock);
             MesonField mfCacheBlock(mfCache.data(),nmom,ngamma,nt,N_iii,N_jjj);    
 
-            t_contr-=usecond();
-            makeBlock(mfCacheBlock, &w[i+ii], &v[j+jj], gamma_, ph, 
-                      env().getNd() - 1, t_int_0, t_int_1, t_int_2, t_int_3);
-            t_contr+=usecond();
+            startTimer("contraction: total");
+            makeMesonFieldBlock(mfCacheBlock, &w[i+ii], &v[j+jj], gamma_, ph, 
+                                env().getNd() - 1, this);
+            stopTimer("contraction: total");
             
             // flops for general N_c & N_s
             flops += vol * ( 2 * 8.0 + 6.0 + 8.0*nmom) * N_iii*N_jjj*ngamma;
             bytes += vol * (12.0 * sizeof(Complex) ) * N_iii*N_jjj
                 +  vol * ( 2.0 * sizeof(Complex) *nmom ) * N_iii*N_jjj* ngamma;
 
-            MODULE_TIMER("Cache copy");
+            startTimer("cache copy");
             for(int iii=0;iii< N_iii;iii++)
             for(int jjj=0;jjj< N_jjj;jjj++)
             for(int m =0;m< nmom;m++)
@@ -315,223 +304,44 @@ void TA2AMesonField<FImpl>::execute(void)
             {
                 mfBlock(m,g,t,ii+iii,jj+jjj) = mfCacheBlock(m,g,t,iii,jjj);
             }
+            stopTimer("cache copy");
         }
 
         // IO
+        if (!par().output.empty())
+        {
         double blockSize, ioTime;
         
-        MODULE_TIMER("IO");
+            LOG(Message) << "Writing block to disk" << std::endl;
+            ioTime = -getDTimer("IO: write block");
+            startTimer("IO: total");
         for(int m = 0; m < nmom; m++)
         for(int g = 0; g < ngamma; g++)
         {
             if ((i == 0) and (j == 0))
             {
+                    startTimer("IO: file creation");
                 initFile(m, g);
+                    stopTimer("IO: file creation");
             }
+                startTimer("IO: write block");
             saveBlock(mfBlock, m, g, i, j);
+                stopTimer("IO: write block");
         }
-        blockSize = static_cast<double>(nmom*ngamma*nt*N_ii*N_jj*sizeof(Complex));
-        ioTime    = static_cast<double>(this->getTimer("IO").count());
-        LOG(Message) << "HDF5 IO " << blockSize/ioTime*1.0e6/1024/1024
+            stopTimer("IO: total");
+            blockSize  = static_cast<double>(nmom*ngamma*nt*N_ii*N_jj*sizeof(Complex));
+            ioTime    += getDTimer("IO: write block");
+            LOG(Message) << "HDF5 IO done " << blockSize/ioTime*1.0e6/1024/1024
                      << " MB/s" << std::endl;
     }
+    }
 
     double nodes    = env().getGrid()->NodeCount();
-    double t_kernel = t_int_0 + t_int_1;
+    double t_kernel = getDTimer("contraction: colour trace & mom.")
+                      + getDTimer("contraction: local space sum");
 
-    LOG(Message) << "Perf " << flops/(t_kernel)/1.0e3/nodes << " Gflop/s/node "  << std::endl;
-    LOG(Message) << "Perf " << bytes/(t_kernel)/1.0e3/nodes << " GB/s/node "  << std::endl;
-}
-
-//////////////////////////////////////////////////////////////////////////////////
-// Cache blocked arithmetic routine
-// Could move to Grid ???
-//////////////////////////////////////////////////////////////////////////////////
-template <typename FImpl>
-void TA2AMesonField<FImpl>::makeBlock(MesonField &mat, 
-					                  const FermionField *lhs_wi,
-					                  const FermionField *rhs_vj,
-					                  std::vector<Gamma::Algebra> gamma,
-					                  const std::vector<LatticeComplex> &mom,
-					                  int orthogdim,
-					                  double &t0,
-					                  double &t1,
-					                  double &t2,
-					                  double &t3) 
-{
-    typedef typename FImpl::SiteSpinor vobj;
-
-    typedef typename vobj::scalar_object sobj;
-    typedef typename vobj::scalar_type scalar_type;
-    typedef typename vobj::vector_type vector_type;
-
-    typedef iSpinMatrix<vector_type> SpinMatrix_v;
-    typedef iSpinMatrix<scalar_type> SpinMatrix_s;
-    
-    int Lblock = mat.dimension(3); 
-    int Rblock = mat.dimension(4);
-
-    GridBase *grid = lhs_wi[0]._grid;
-    
-    const int    Nd = grid->_ndimension;
-    const int Nsimd = grid->Nsimd();
-
-    int Nt     = grid->GlobalDimensions()[orthogdim];
-    int Ngamma = gamma.size();
-    int Nmom   = mom.size();
-
-    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
-    int MFrvol = rd*Lblock*Rblock*Nmom;
-    int MFlvol = ld*Lblock*Rblock*Nmom;
-
-    Vector<SpinMatrix_v > lvSum(MFrvol);
-    parallel_for (int r = 0; r < MFrvol; r++)
-    {
-        lvSum[r] = zero;
-    }
-
-    Vector<SpinMatrix_s > lsSum(MFlvol);             
-    parallel_for (int r = 0; r < MFlvol; r++)
-    {
-        lsSum[r]=scalar_type(0.0);
-    }
-
-    int e1=    grid->_slice_nblock[orthogdim];
-    int e2=    grid->_slice_block [orthogdim];
-    int stride=grid->_slice_stride[orthogdim];
-
-    t0-=usecond();
-    MODULE_TIMER("Colour trace * mom.");
-    // Nested parallelism would be ok
-    // Wasting cores here. Test case r
-    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_wi[i]._odata[ss]);
-
-                for(int j=0;j<Rblock;j++)
-                {
-                    SpinMatrix_v vv;
-                    auto right = rhs_vj[j]._odata[ss];
-
-                    for(int s1=0;s1<Ns;s1++)
-                    for(int s2=0;s2<Ns;s2++)
-                    {
-                        vv()(s1,s2)() = left()(s2)(0) * right()(s1)(0)
-                                        + left()(s2)(1) * right()(s1)(1)
-                                        + left()(s2)(2) * right()(s1)(2);
-                    }
-                    
-                    // After getting the sitewise product do the mom phase loop
-                    int base = Nmom*i+Nmom*Lblock*j+Nmom*Lblock*Rblock*r;
-
-                    for ( int m=0;m<Nmom;m++)
-                    {
-                        int idx = m+base;
-                        auto phase = mom[m]._odata[ss];
-                        mac(&lvSum[idx],&vv,&phase);
-                    }
-                }
-            }
-        }
-    }
-    t0+=usecond();
-
-    // Sum across simd lanes in the plane, breaking out orthog dir.
-    MODULE_TIMER("Local space sum");
-    t1-=usecond();
-    parallel_for(int rt=0;rt<rd;rt++)
-    {
-        std::vector<int> icoor(Nd);
-        std::vector<SpinMatrix_s> extracted(Nsimd);               
-
-        for(int i=0;i<Lblock;i++)
-        for(int j=0;j<Rblock;j++)
-        for(int m=0;m<Nmom;m++)
-        {
-
-            int ij_rdx = m+Nmom*i+Nmom*Lblock*j+Nmom*Lblock*Rblock*rt;
-
-            extract(lvSum[ij_rdx],extracted);
-            for(int idx=0;idx<Nsimd;idx++)
-            {
-                grid->iCoorFromIindex(icoor,idx);
-
-                int ldx    = rt+icoor[orthogdim]*rd;
-                int ij_ldx = m+Nmom*i+Nmom*Lblock*j+Nmom*Lblock*Rblock*ldx;
-
-                lsSum[ij_ldx]=lsSum[ij_ldx]+extracted[idx];
-            }
-        }
-    }
-    t1+=usecond();
-    assert(mat.dimension(0) == Nmom);
-    assert(mat.dimension(1) == Ngamma);
-    assert(mat.dimension(2) == Nt);
-    t2-=usecond();
-
-    // ld loop and local only??
-    MODULE_TIMER("Spin trace");
-    int pd = grid->_processors[orthogdim];
-    int pc = grid->_processor_coor[orthogdim];
-    parallel_for_nest2(int lt=0;lt<ld;lt++)
-    {
-        for(int pt=0;pt<pd;pt++)
-        {
-            int t = lt + pt*ld;
-            if (pt == pc)
-            {
-                for(int i=0;i<Lblock;i++)
-                for(int j=0;j<Rblock;j++)
-                for(int m=0;m<Nmom;m++)
-                {
-                    int ij_dx = m+Nmom*i + Nmom*Lblock * j + Nmom*Lblock * Rblock * lt;
-
-                    for(int mu=0;mu<Ngamma;mu++)
-                    {
-                        // this is a bit slow
-                        mat(m,mu,t,i,j) = trace(lsSum[ij_dx]*Gamma(gamma[mu]));
-                    }
-                }
-            } 
-            else 
-            { 
-                const scalar_type zz(0.0);
-
-                for(int i=0;i<Lblock;i++)
-                for(int j=0;j<Rblock;j++)
-                for(int mu=0;mu<Ngamma;mu++)
-                for(int m=0;m<Nmom;m++)
-                {
-                    mat(m,mu,t,i,j) =zz;
-                }
-            }
-        }
-    }
-    t2+=usecond();
-    ////////////////////////////////////////////////////////////////////
-    // This global sum is taking as much as 50% of time on 16 nodes
-    // Vector size is 7 x 16 x 32 x 16 x 16 x sizeof(complex) = 2MB - 60MB depending on volume
-    // Healthy size that should suffice
-    ////////////////////////////////////////////////////////////////////
-    t3-=usecond();
-    MODULE_TIMER("Global sum");
-    grid->GlobalSumVector(&mat(0,0,0,0,0),Nmom*Ngamma*Nt*Lblock*Rblock);
-    t3+=usecond();
+    LOG(Message) << "Perf " << flops/t_kernel/1.0e3/nodes << " Gflop/s/node "  << std::endl;
+    LOG(Message) << "Perf " << bytes/t_kernel/1.0e3/nodes << " GB/s/node "  << std::endl;
 }
 
 // IO
diff --git a/extras/Hadrons/Modules/MContraction/A2AMesonFieldKernels.hpp b/extras/Hadrons/Modules/MContraction/A2AMesonFieldKernels.hpp
new file mode 100644
index 00000000..4d06f796
--- /dev/null
+++ b/extras/Hadrons/Modules/MContraction/A2AMesonFieldKernels.hpp
@@ -0,0 +1,226 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MContraction/A2AMesonFieldKernels.hpp
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#ifndef Hadrons_MContraction_A2AMesonFieldKernels_hpp_
+#define Hadrons_MContraction_A2AMesonFieldKernels_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Eigen/unsupported/CXX11/Tensor>
+
+BEGIN_HADRONS_NAMESPACE
+
+BEGIN_MODULE_NAMESPACE(MContraction)
+
+typedef Eigen::TensorMap<Eigen::Tensor<Complex, 5, Eigen::RowMajor>> MesonField;
+
+////////////////////////////////////////////////////////////////////////////////
+// Cache blocked arithmetic routine
+// Could move to Grid ???
+////////////////////////////////////////////////////////////////////////////////
+template <typename Field>
+void makeMesonFieldBlock(MesonField &mat, 
+                         const Field *lhs_wi,
+                         const Field *rhs_vj,
+                         std::vector<Gamma::Algebra> gamma,
+                         const std::vector<LatticeComplex> &mom,
+                         int orthogdim,
+                         ModuleBase *caller = nullptr) 
+{
+    typedef typename Field::vector_object vobj;
+    typedef typename vobj::scalar_object  sobj;
+    typedef typename vobj::scalar_type    scalar_type;
+    typedef typename vobj::vector_type    vector_type;
+
+    typedef iSpinMatrix<vector_type> SpinMatrix_v;
+    typedef iSpinMatrix<scalar_type> SpinMatrix_s;
+    
+    int Lblock = mat.dimension(3); 
+    int Rblock = mat.dimension(4);
+
+    GridBase *grid = lhs_wi[0]._grid;
+    
+    const int    Nd = grid->_ndimension;
+    const int Nsimd = grid->Nsimd();
+
+    int Nt     = grid->GlobalDimensions()[orthogdim];
+    int Ngamma = gamma.size();
+    int Nmom   = mom.size();
+
+    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
+    int MFrvol = rd*Lblock*Rblock*Nmom;
+    int MFlvol = ld*Lblock*Rblock*Nmom;
+
+    Vector<SpinMatrix_v > lvSum(MFrvol);
+    parallel_for (int r = 0; r < MFrvol; r++)
+    {
+        lvSum[r] = zero;
+    }
+
+    Vector<SpinMatrix_s > lsSum(MFlvol);             
+    parallel_for (int r = 0; r < MFlvol; r++)
+    {
+        lsSum[r]=scalar_type(0.0);
+    }
+
+    int e1=    grid->_slice_nblock[orthogdim];
+    int e2=    grid->_slice_block [orthogdim];
+    int stride=grid->_slice_stride[orthogdim];
+
+    if (caller) caller->startTimer("contraction: colour trace & mom.");
+    // Nested parallelism would be ok
+    // Wasting cores here. Test case r
+    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_wi[i]._odata[ss]);
+
+                for(int j=0;j<Rblock;j++)
+                {
+                    SpinMatrix_v vv;
+                    auto right = rhs_vj[j]._odata[ss];
+
+                    for(int s1=0;s1<Ns;s1++)
+                    for(int s2=0;s2<Ns;s2++)
+                    {
+                        vv()(s1,s2)() = left()(s2)(0) * right()(s1)(0)
+                                        + left()(s2)(1) * right()(s1)(1)
+                                        + left()(s2)(2) * right()(s1)(2);
+                    }
+                    
+                    // After getting the sitewise product do the mom phase loop
+                    int base = Nmom*i+Nmom*Lblock*j+Nmom*Lblock*Rblock*r;
+
+                    for ( int m=0;m<Nmom;m++)
+                    {
+                        int idx = m+base;
+                        auto phase = mom[m]._odata[ss];
+                        mac(&lvSum[idx],&vv,&phase);
+                    }
+                }
+            }
+        }
+    }
+    if (caller) caller->stopTimer("contraction: colour trace & mom.");
+
+    // Sum across simd lanes in the plane, breaking out orthog dir.
+    if (caller) caller->startTimer("contraction: local space sum");
+    parallel_for(int rt=0;rt<rd;rt++)
+    {
+        std::vector<int> icoor(Nd);
+        std::vector<SpinMatrix_s> extracted(Nsimd);               
+
+        for(int i=0;i<Lblock;i++)
+        for(int j=0;j<Rblock;j++)
+        for(int m=0;m<Nmom;m++)
+        {
+
+            int ij_rdx = m+Nmom*i+Nmom*Lblock*j+Nmom*Lblock*Rblock*rt;
+
+            extract(lvSum[ij_rdx],extracted);
+            for(int idx=0;idx<Nsimd;idx++)
+            {
+                grid->iCoorFromIindex(icoor,idx);
+
+                int ldx    = rt+icoor[orthogdim]*rd;
+                int ij_ldx = m+Nmom*i+Nmom*Lblock*j+Nmom*Lblock*Rblock*ldx;
+
+                lsSum[ij_ldx]=lsSum[ij_ldx]+extracted[idx];
+            }
+        }
+    }
+    if (caller) caller->stopTimer("contraction: local space sum");
+
+    // ld loop and local only??
+    if (caller) caller->startTimer("contraction: spin trace");
+    int pd = grid->_processors[orthogdim];
+    int pc = grid->_processor_coor[orthogdim];
+    parallel_for_nest2(int lt=0;lt<ld;lt++)
+    {
+        for(int pt=0;pt<pd;pt++)
+        {
+            int t = lt + pt*ld;
+            if (pt == pc)
+            {
+                for(int i=0;i<Lblock;i++)
+                for(int j=0;j<Rblock;j++)
+                for(int m=0;m<Nmom;m++)
+                {
+                    int ij_dx = m+Nmom*i + Nmom*Lblock * j + Nmom*Lblock * Rblock * lt;
+
+                    for(int mu=0;mu<Ngamma;mu++)
+                    {
+                        // this is a bit slow
+                        mat(m,mu,t,i,j) = trace(lsSum[ij_dx]*Gamma(gamma[mu]));
+                    }
+                }
+            } 
+            else 
+            { 
+                const scalar_type zz(0.0);
+
+                for(int i=0;i<Lblock;i++)
+                for(int j=0;j<Rblock;j++)
+                for(int mu=0;mu<Ngamma;mu++)
+                for(int m=0;m<Nmom;m++)
+                {
+                    mat(m,mu,t,i,j) =zz;
+                }
+            }
+        }
+    }
+    if (caller) caller->stopTimer("contraction: spin trace");
+    ////////////////////////////////////////////////////////////////////
+    // This global sum is taking as much as 50% of time on 16 nodes
+    // Vector size is 7 x 16 x 32 x 16 x 16 x sizeof(complex) = 2MB - 60MB depending on volume
+    // Healthy size that should suffice
+    ////////////////////////////////////////////////////////////////////
+    if (caller) caller->startTimer("contraction: global sum");
+    grid->GlobalSumVector(&mat(0,0,0,0,0),Nmom*Ngamma*Nt*Lblock*Rblock);
+    if (caller) caller->stopTimer("contraction: global sum");
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif //Hadrons_MContraction_A2AMesonField_hpp_
diff --git a/extras/Hadrons/Modules/MFermion/FreeProp.cc b/extras/Hadrons/Modules/MFermion/FreeProp.cc
index 7b3a0ce5..0ff07887 100644
--- a/extras/Hadrons/Modules/MFermion/FreeProp.cc
+++ b/extras/Hadrons/Modules/MFermion/FreeProp.cc
@@ -6,6 +6,7 @@ Source file: extras/Hadrons/Modules/MFermion/FreeProp.cc
 
 Copyright (C) 2015-2018
 
+Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Vera Guelpers <vmg1n14@soton.ac.uk>
 
 This program is free software; you can redistribute it and/or modify
diff --git a/extras/Hadrons/Modules/MGauge/UnitEm.cc b/extras/Hadrons/Modules/MGauge/UnitEm.cc
index 6b3be47e..5c076592 100644
--- a/extras/Hadrons/Modules/MGauge/UnitEm.cc
+++ b/extras/Hadrons/Modules/MGauge/UnitEm.cc
@@ -6,6 +6,7 @@ Source file: extras/Hadrons/Modules/MGauge/UnitEm.cc
 
 Copyright (C) 2015-2018
 
+Author: Antonin Portelli <antonin.portelli@me.com>
 Author: James Harrison <jch1g10@soton.ac.uk>
 
 This program is free software; you can redistribute it and/or modify
diff --git a/extras/Hadrons/Modules/MGauge/UnitEm.hpp b/extras/Hadrons/Modules/MGauge/UnitEm.hpp
index a9cb9dfa..c6c89bb0 100644
--- a/extras/Hadrons/Modules/MGauge/UnitEm.hpp
+++ b/extras/Hadrons/Modules/MGauge/UnitEm.hpp
@@ -6,6 +6,7 @@ Source file: extras/Hadrons/Modules/MGauge/UnitEm.hpp
 
 Copyright (C) 2015-2018
 
+Author: Antonin Portelli <antonin.portelli@me.com>
 Author: James Harrison <jch1g10@soton.ac.uk>
 
 This program is free software; you can redistribute it and/or modify
diff --git a/extras/Hadrons/Modules/MScalar/ScalarVP.cc b/extras/Hadrons/Modules/MScalar/ScalarVP.cc
index d560ec5a..b43b0c9f 100644
--- a/extras/Hadrons/Modules/MScalar/ScalarVP.cc
+++ b/extras/Hadrons/Modules/MScalar/ScalarVP.cc
@@ -6,6 +6,7 @@ Source file: extras/Hadrons/Modules/MScalar/ScalarVP.cc
 
 Copyright (C) 2015-2018
 
+Author: Antonin Portelli <antonin.portelli@me.com>
 Author: James Harrison <jch1g10@soton.ac.uk>
 
 This program is free software; you can redistribute it and/or modify
diff --git a/extras/Hadrons/Modules/MScalar/ScalarVP.hpp b/extras/Hadrons/Modules/MScalar/ScalarVP.hpp
index e1b2dc3d..191348d3 100644
--- a/extras/Hadrons/Modules/MScalar/ScalarVP.hpp
+++ b/extras/Hadrons/Modules/MScalar/ScalarVP.hpp
@@ -6,6 +6,7 @@ Source file: extras/Hadrons/Modules/MScalar/ScalarVP.hpp
 
 Copyright (C) 2015-2018
 
+Author: Antonin Portelli <antonin.portelli@me.com>
 Author: James Harrison <jch1g10@soton.ac.uk>
 
 This program is free software; you can redistribute it and/or modify
diff --git a/extras/Hadrons/Modules/MScalar/VPCounterTerms.cc b/extras/Hadrons/Modules/MScalar/VPCounterTerms.cc
index 86867788..4c3fab85 100644
--- a/extras/Hadrons/Modules/MScalar/VPCounterTerms.cc
+++ b/extras/Hadrons/Modules/MScalar/VPCounterTerms.cc
@@ -6,6 +6,7 @@ Source file: extras/Hadrons/Modules/MScalar/VPCounterTerms.cc
 
 Copyright (C) 2015-2018
 
+Author: Antonin Portelli <antonin.portelli@me.com>
 Author: James Harrison <jch1g10@soton.ac.uk>
 
 This program is free software; you can redistribute it and/or modify
diff --git a/extras/Hadrons/Modules/MScalar/VPCounterTerms.hpp b/extras/Hadrons/Modules/MScalar/VPCounterTerms.hpp
index a2ccbb56..92d37686 100644
--- a/extras/Hadrons/Modules/MScalar/VPCounterTerms.hpp
+++ b/extras/Hadrons/Modules/MScalar/VPCounterTerms.hpp
@@ -6,6 +6,7 @@ Source file: extras/Hadrons/Modules/MScalar/VPCounterTerms.hpp
 
 Copyright (C) 2015-2018
 
+Author: Antonin Portelli <antonin.portelli@me.com>
 Author: James Harrison <jch1g10@soton.ac.uk>
 
 This program is free software; you can redistribute it and/or modify
diff --git a/extras/Hadrons/Modules/MSolver/A2AVectors.hpp b/extras/Hadrons/Modules/MSolver/A2AVectors.hpp
index 6d150499..7022b90b 100644
--- a/extras/Hadrons/Modules/MSolver/A2AVectors.hpp
+++ b/extras/Hadrons/Modules/MSolver/A2AVectors.hpp
@@ -178,7 +178,7 @@ void TA2AVectors<FImpl, Pack>::execute(void)
     {
         auto &epack  = envGet(Pack, par().eigenPack);
 
-        MODULE_TIMER("V low mode");
+        startTimer("V low mode");
         LOG(Message) << "V vector i = " << il << " (low mode)" << std::endl;
         if (Ls == 1)
         {
@@ -189,7 +189,8 @@ void TA2AVectors<FImpl, Pack>::execute(void)
             envGetTmp(FermionField, f5);
             a2a.makeLowModeV5D(v[il], f5, epack.evec[il], epack.eval[il]);
         }
-        MODULE_TIMER("W low mode");
+        stopTimer("V low mode");
+        startTimer("W low mode");
         LOG(Message) << "W vector i = " << il << " (low mode)" << std::endl;
         if (Ls == 1)
         {
@@ -199,13 +200,14 @@ void TA2AVectors<FImpl, Pack>::execute(void)
         {
             envGetTmp(FermionField, f5);
             a2a.makeLowModeW5D(w[il], f5, epack.evec[il], epack.eval[il]);
-        }   
+        }
+        stopTimer("W low mode");
     }
 
     // High modes
     for (unsigned int ih = 0; ih < noise.size(); ih++)
     {
-        MODULE_TIMER("V high mode");
+        startTimer("V high mode");
         LOG(Message) << "V vector i = " << Nl_ + ih
                      << " (" << ((Nl_ > 0) ? "high " : "") 
                      << "stochastic mode)" << std::endl;
@@ -218,7 +220,8 @@ void TA2AVectors<FImpl, Pack>::execute(void)
             envGetTmp(FermionField, f5);
             a2a.makeHighModeV5D(v[Nl_ + ih], f5, noise[ih]);
         }
-        MODULE_TIMER("W high mode");
+        stopTimer("V high mode");
+        startTimer("W high mode");
         LOG(Message) << "W vector i = " << Nl_ + ih
                      << " (" << ((Nl_ > 0) ? "high " : "") 
                      << "stochastic mode)" << std::endl;
@@ -231,6 +234,7 @@ void TA2AVectors<FImpl, Pack>::execute(void)
             envGetTmp(FermionField, f5);
             a2a.makeHighModeW5D(w[Nl_ + ih], f5, noise[ih]);
         }
+        stopTimer("W high mode");
     }
 }
 
diff --git a/extras/Hadrons/Modules/MUtilities/RandomVectors.cc b/extras/Hadrons/Modules/MUtilities/RandomVectors.cc
index cc122c90..7da32322 100644
--- a/extras/Hadrons/Modules/MUtilities/RandomVectors.cc
+++ b/extras/Hadrons/Modules/MUtilities/RandomVectors.cc
@@ -1,3 +1,30 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MUtilities/RandomVectors.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
 #include <Grid/Hadrons/Modules/MUtilities/RandomVectors.hpp>
 
 using namespace Grid;
diff --git a/extras/Hadrons/Modules/MUtilities/RandomVectors.hpp b/extras/Hadrons/Modules/MUtilities/RandomVectors.hpp
index 590ef224..4da84506 100644
--- a/extras/Hadrons/Modules/MUtilities/RandomVectors.hpp
+++ b/extras/Hadrons/Modules/MUtilities/RandomVectors.hpp
@@ -1,3 +1,30 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MUtilities/RandomVectors.hpp
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
 #ifndef Hadrons_MUtilities_RandomVectors_hpp_
 #define Hadrons_MUtilities_RandomVectors_hpp_
 
diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc
index 6bc07618..eef12b0d 100644
--- a/extras/Hadrons/modules.inc
+++ b/extras/Hadrons/modules.inc
@@ -61,6 +61,7 @@ modules_cc =\
 modules_hpp =\
   Modules/MContraction/Baryon.hpp \
   Modules/MContraction/A2AMesonField.hpp \
+  Modules/MContraction/A2AMesonFieldKernels.hpp \
   Modules/MContraction/Meson.hpp \
   Modules/MContraction/WeakHamiltonian.hpp \
   Modules/MContraction/WeakHamiltonianNonEye.hpp \