diff --git a/Grid/qcd/utils/A2Autils.h b/Grid/qcd/utils/A2Autils.h
index f6db38be..97188ffe 100644
--- a/Grid/qcd/utils/A2Autils.h
+++ b/Grid/qcd/utils/A2Autils.h
@@ -68,16 +68,6 @@ public:
         const std::vector<ComplexField> &emB1,
         int orthogdim, double *t_kernel = nullptr, double *t_gsum = nullptr);
 
-  template <typename TensorType>
-  static void LeptonField(TensorType &mat,
-        const FermionField *lhs_wi,
-        const FermionField *rhs_vj,
-        const std::vector<ComplexField> &lepton02,
-        const std::vector<ComplexField> &lepton03,
-        const std::vector<ComplexField> &lepton12,
-        const std::vector<ComplexField> &lepton13,
-        int orthogdim, double *t_kernel = nullptr, double *t_gsum = nullptr);
-
   static void ContractWWVV(std::vector<PropagatorField> &WWVV,
 			   const Eigen::Tensor<ComplexD,3> &WW_sd,
 			   const FermionField *vs,
@@ -819,204 +809,6 @@ void A2Autils<FImpl>::AslashField(TensorType &mat,
     if (t_gsum) *t_gsum += usecond();
 }
 
-// "Lepton" field w_i(x)^dag * L_mu * g^L_mu * v_j(x)
-//
-// With V-A current:
-//
-// g^L_mu = g_mu * (1 - g5)
-//
-// then in spin space
-//
-//                 ( 0	0	L_02	L_03	)
-// L_mu * g^L_mu = ( 0	0	L_12	L_13 	)
-//                 ( 0	0       0       0	)
-//                 ( 0	0 	0       0	)
-//
-//
-// With
-//
-// L_02 =  2*L3 + 2*i*L2
-// L_03 = -2*L1 + 2*i*L0
-// L_12 =  2*L1 + 2*i*L0
-// L_13 =  2*L3 - 2*i*L2
-//
-// (where mu = (x,y,z,t) )
-//
-template <class FImpl>
-template <typename TensorType>
-void A2Autils<FImpl>::LeptonField(TensorType &mat,
-          const FermionField *lhs_wi,
-          const FermionField *rhs_vj,
-          const std::vector<ComplexField> &lepton02,
-          const std::vector<ComplexField> &lepton03,
-          const std::vector<ComplexField> &lepton12,
-          const std::vector<ComplexField> &lepton13,
-          int orthogdim, double *t_kernel, double *t_gsum)
-{
-    typedef typename FermionField::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;
-    typedef iSinglet<vector_type>    Singlet_v;
-    typedef iSinglet<scalar_type>    Singlet_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 Nlep = lepton02.size();
-    assert(lepton03.size() == Nlep);
-    assert(lepton12.size() == Nlep);
-    assert(lepton13.size() == Nlep);
-
-    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*Nlep;
-    int MFlvol = ld*Lblock*Rblock*Nlep;
-
-    Vector<vector_type> lvSum(MFrvol);
-    parallel_for (int r = 0; r < MFrvol; r++)
-    {
-        lvSum[r] = zero;
-    }
-
-    Vector<scalar_type> 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];
-
-    // Nested parallelism would be ok
-    // Wasting cores here. Test case r
-    if (t_kernel) *t_kernel = -usecond();
-    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 = Nlep*i+Nlep*Lblock*j+Nlep*Lblock*Rblock*r;
-
-                    for ( int m=0;m<Nlep;m++)
-                    {
-                        int idx  = m+base;
-
-                        auto lep02  = lepton02[m]._odata[ss];
-                        auto lep12  = lepton12[m]._odata[ss];
-                        auto lep03  = lepton03[m]._odata[ss];
-                        auto lep13  = lepton13[m]._odata[ss];
-
-			lvSum[idx] += vv()(2,0)()*lep02()()() + vv()(3,0)()*lep03()()()
-					+ vv()(2,1)()*lep12()()() + vv()(3,1)()*lep13()()();
-                    }
-                }
-            }
-        }
-    }
-    // Sum across simd lanes in the plane, breaking out orthog dir.
-    parallel_for(int rt=0;rt<rd;rt++)
-    {
-        std::vector<int> icoor(Nd);
-        std::vector<scalar_type> extracted(Nsimd);
-
-        for(int i=0;i<Lblock;i++)
-        for(int j=0;j<Rblock;j++)
-        for(int m=0;m<Nlep;m++)
-        {
-
-            int ij_rdx = m+Nlep*i+Nlep*Lblock*j+Nlep*Lblock*Rblock*rt;
-
-            extract<vector_type,scalar_type>(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+Nlep*i+Nlep*Lblock*j+Nlep*Lblock*Rblock*ldx;
-
-                lsSum[ij_ldx]=lsSum[ij_ldx]+extracted[idx];
-            }
-        }
-    }
-    if (t_kernel) *t_kernel += usecond();
-
-    // ld loop and local only??
-    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<Nlep;m++)
-                {
-                    int ij_dx = m+Nlep*i + Nlep*Lblock * j + Nlep*Lblock * Rblock * lt;
-
-                    mat(m,0,t,i,j) = lsSum[ij_dx];
-                }
-            }
-            else
-            {
-                const scalar_type zz(0.0);
-
-                for(int i=0;i<Lblock;i++)
-                for(int j=0;j<Rblock;j++)
-                for(int m=0;m<Nlep;m++)
-                {
-                    mat(m,0,t,i,j) = zz;
-                }
-            }
-        }
-    }
-    if (t_gsum) *t_gsum = -usecond();
-    grid->GlobalSumVector(&mat(0,0,0,0,0),Nlep*Nt*Lblock*Rblock);
-    if (t_gsum) *t_gsum += usecond();
-}
-
-
 ////////////////////////////////////////////
 // Schematic thoughts about more generalised four quark insertion
 //
diff --git a/Hadrons/Modules.hpp b/Hadrons/Modules.hpp
index b0aae2a8..898cb532 100644
--- a/Hadrons/Modules.hpp
+++ b/Hadrons/Modules.hpp
@@ -3,7 +3,6 @@
 #include <Hadrons/Modules/MContraction/A2AAslashField.hpp>
 #include <Hadrons/Modules/MContraction/A2ALoop.hpp>
 #include <Hadrons/Modules/MContraction/A2AMesonField.hpp>
-#include <Hadrons/Modules/MContraction/A2ALeptonField.hpp>
 #include <Hadrons/Modules/MContraction/Meson.hpp>
 #include <Hadrons/Modules/MContraction/DiscLoop.hpp>
 #include <Hadrons/Modules/MContraction/Gamma3pt.hpp>
diff --git a/Hadrons/Modules/MContraction/A2ALeptonField.cc b/Hadrons/Modules/MContraction/A2ALeptonField.cc
deleted file mode 100644
index 3cc3d86a..00000000
--- a/Hadrons/Modules/MContraction/A2ALeptonField.cc
+++ /dev/null
@@ -1,34 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: Hadrons/Modules/MContraction/A2ALeptonField.cc
-
-Copyright (C) 2015-2019
-
-Author: Vera Guelpers <Vera.Guelpers@ed.ac.uk>
-
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-the Free Software Foundation; either version 2 of the License, or
-(at your option) any later version.
-
-This program is distributed in the hope that it will be useful,
-but WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-GNU General Public License for more details.
-
-You should have received a copy of the GNU General Public License along
-with this program; if not, write to the Free Software Foundation, Inc.,
-51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-See the full license in the file "LICENSE" in the top level distribution directory
-*************************************************************************************/
-/*  END LEGAL */
-#include <Hadrons/Modules/MContraction/A2ALeptonField.hpp>
-
-using namespace Grid;
-using namespace Hadrons;
-using namespace MContraction;
-
-template class Grid::Hadrons::MContraction::TA2ALeptonField<FIMPL, PhotonR>;
diff --git a/Hadrons/Modules/MContraction/A2ALeptonField.hpp b/Hadrons/Modules/MContraction/A2ALeptonField.hpp
deleted file mode 100644
index 7e041cf5..00000000
--- a/Hadrons/Modules/MContraction/A2ALeptonField.hpp
+++ /dev/null
@@ -1,283 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: Hadrons/Modules/MContraction/A2ALeptonField.cc
-
-Copyright (C) 2015-2019
-
-Author: Vera Guelpers <Vera.Guelpers@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_A2ALeptonField_hpp_
-#define Hadrons_MContraction_A2ALeptonField_hpp_
-
-#include <Hadrons/Global.hpp>
-#include <Hadrons/Module.hpp>
-#include <Hadrons/ModuleFactory.hpp>
-#include <Hadrons/A2AMatrix.hpp>
-
-#ifndef MLF_IO_TYPE
-#define MLF_IO_TYPE ComplexF
-#endif
-
-BEGIN_HADRONS_NAMESPACE
-
-/******************************************************************************
- *                         A2ALeptonField                                 *
- ******************************************************************************/
-BEGIN_MODULE_NAMESPACE(MContraction)
-
-class A2ALeptonFieldPar: Serializable
-{
-public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(A2ALeptonFieldPar,
-                                    int, cacheBlock,
-                                    int, block,
-                                    std::string, left,
-                                    std::string, right,
-                                    std::string, output,
-                                    std::vector<std::string>, lepton);
-};
-
-class A2ALeptonFieldMetadata: Serializable
-{
-public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(A2ALeptonFieldMetadata,
-                                    std::string, leptonName,
-				    int, sidx1,
-				    int, sidx2);
-};
-
-template <typename T, typename FImpl>
-class LeptonFieldKernel: public A2AKernel<T, typename FImpl::FermionField>
-{
-public:
-    typedef typename FImpl::FermionField FermionField;
-public:
-    LeptonFieldKernel(const std::vector<LatticeComplex> &lepton02,
-                      const std::vector<LatticeComplex> &lepton03,
-                      const std::vector<LatticeComplex> &lepton12,
-                      const std::vector<LatticeComplex> &lepton13,
-                      GridBase *grid)
-    : lepton02_(lepton02), lepton03_(lepton03), lepton12_(lepton12), lepton13_(lepton13), grid_(grid)
-    {
-        vol_ = 1.;
-        for (auto &d: grid_->GlobalDimensions())
-        {
-            vol_ *= d;
-        }
-    }
-
-    virtual ~LeptonFieldKernel(void) = default;
-    virtual void operator()(A2AMatrixSet<T> &m, const FermionField *left, 
-                            const FermionField *right,
-                            const unsigned int orthogDim, double &t)
-    {
-        A2Autils<FImpl>::LeptonField(m, left, right, lepton02_, lepton03_, lepton12_, lepton13_, orthogDim, &t);
-    }
-
-    virtual double flops(const unsigned int blockSizei, const unsigned int blockSizej)
-    {
-        return 0.;
-    }
-
-    virtual double bytes(const unsigned int blockSizei, const unsigned int blockSizej)
-    {
-        return 0.;
-    }
-private:
-    const std::vector<LatticeComplex> &lepton02_, &lepton03_, &lepton12_, &lepton13_;
-    GridBase                          *grid_;
-    double                            vol_;
-};
-
-template <typename FImpl, typename PhotonImpl>
-class TA2ALeptonField: public Module<A2ALeptonFieldPar>
-{
-public:
-    FERM_TYPE_ALIASES(FImpl,);
-    typedef A2AMatrixBlockComputation<Complex, 
-                                      FermionField, 
-                                      A2ALeptonFieldMetadata, 
-                                      MLF_IO_TYPE> Computation;
-    typedef LeptonFieldKernel<Complex, FImpl> Kernel;
-public:
-    // constructor
-    TA2ALeptonField(const std::string name);
-    // destructor
-    virtual ~TA2ALeptonField(void) {};
-    // dependency relation
-    virtual std::vector<std::string> getInput(void);
-    virtual std::vector<std::string> getOutput(void);
-    // setup
-    virtual void setup(void);
-    // execution
-    virtual void execute(void);
-};
-
-MODULE_REGISTER_TMP(A2ALeptonField, ARG(TA2ALeptonField<FIMPL, PhotonR>), MContraction);
-
-/******************************************************************************
- *                 TA2ALeptonField implementation                             *
- ******************************************************************************/
-// constructor /////////////////////////////////////////////////////////////////
-template <typename FImpl, typename PhotonImpl>
-TA2ALeptonField<FImpl, PhotonImpl>::TA2ALeptonField(const std::string name)
-: Module<A2ALeptonFieldPar>(name)
-{}
-
-// dependencies/products ///////////////////////////////////////////////////////
-template <typename FImpl, typename PhotonImpl>
-std::vector<std::string> TA2ALeptonField<FImpl, PhotonImpl>::getInput(void)
-{
-    std::vector<std::string> in = par().lepton;
-    in.push_back(par().left);
-    in.push_back(par().right);
-
-    return in;
-}
-
-template <typename FImpl, typename PhotonImpl>
-std::vector<std::string> TA2ALeptonField<FImpl, PhotonImpl>::getOutput(void)
-{
-    std::vector<std::string> out = {};
-    
-    return out;
-}
-
-// setup ///////////////////////////////////////////////////////////////////////
-template <typename FImpl, typename PhotonImpl>
-void TA2ALeptonField<FImpl, PhotonImpl>::setup(void)
-{
-    envTmp(Computation, "computation", 1, envGetGrid(FermionField), 
-           env().getNd() - 1, par().lepton.size()*8, 1, par().block, 
-           par().cacheBlock, this);
-    envTmp(std::vector<LatticeComplex>, "Lmu", 1, 
-           4, envGetGrid(LatticeComplex));
-    envTmp(std::vector<LatticeComplex>, "L02", 1, 
-           par().lepton.size()*8, envGetGrid(LatticeComplex));
-    envTmp(std::vector<LatticeComplex>, "L03", 1, 
-           par().lepton.size()*8, envGetGrid(LatticeComplex)); 
-    envTmp(std::vector<LatticeComplex>, "L12", 1, 
-           par().lepton.size()*8, envGetGrid(LatticeComplex));
-    envTmp(std::vector<LatticeComplex>, "L13", 1, 
-           par().lepton.size()*8, envGetGrid(LatticeComplex));
-    envTmpLat(PropagatorField, "prop_buf");
-
-}
-
-// execution ///////////////////////////////////////////////////////////////////
-template <typename FImpl, typename PhotonImpl>
-void TA2ALeptonField<FImpl, PhotonImpl>::execute(void)
-{
-    auto &left  = envGet(std::vector<FermionField>, par().left);
-    auto &right = envGet(std::vector<FermionField>, par().right);
-
-    int nt         = env().getDim().back();
-    int N_i        = left.size();
-    int N_j        = right.size();
-    int block      = par().block;
-    int cacheBlock = par().cacheBlock;
-
-    LOG(Message) << "Computing all-to-all lepton insertion fields" << std::endl;
-    LOG(Message) << "Left: '" << par().left << "' Right: '" << par().right << "'" << std::endl;
-    LOG(Message) << "leptons: " << std::endl;
-    for (auto &name: par().lepton)
-    {
-        LOG(Message) << "  " << name << std::endl;
-    }
-    LOG(Message) << "lepton insertion field size: " << nt << "*" << N_i << "*" << N_j 
-                 << " (filesize " << sizeString(nt*N_i*N_j*sizeof(MLF_IO_TYPE)) 
-                 << "/spin index/lepton)" << std::endl;
-
-    envGetTmp(std::vector<ComplexField>, L02);
-    envGetTmp(std::vector<ComplexField>, L03);
-    envGetTmp(std::vector<ComplexField>, L12);
-    envGetTmp(std::vector<ComplexField>, L13);
-
-    unsigned int ii = 0;
-    envGetTmp(PropagatorField, prop_buf);
-    envGetTmp(std::vector<LatticeComplex>, Lmu);
-    for (unsigned int i = 0; i < par().lepton.size(); ++i)
-    {
-	auto &lepton = envGet(PropagatorField, par().lepton[i]);
-	// need only half of the spin indices (s1s2) of the lepton field
-	// the other ones are zero because of the left-handed current
-	for (unsigned int s1 = 0; s1 < 2; ++s1) 
-	for (unsigned int s2 = 0; s2 < 4; ++s2)
-	{  	
-	    for (unsigned int mu = 0; mu < 4; ++mu)
-	    {
-		prop_buf = GammaL(Gamma::gmu[mu]) * lepton;
-		//lepon is unit matix in color space, so just pick one diagonal enty
-		Lmu[mu] = peekColour(peekSpin(prop_buf,s1,s2),0,0);
-	    }
-	    // build the required combinations of lepton fields for the a2a kernel
-	    // [g^L_mu * L_mu]_{ab} = L_{ab} (with mu = (x,y,z,t))
-	    L02[ii] = 2.0*Lmu[3] + 2.0*timesI(Lmu[2]);
-	    L03[ii] = 2.0*timesI(Lmu[0]) - 2.0*Lmu[1];
-	    L12[ii] = 2.0*timesI(Lmu[0]) + 2.0*Lmu[1];
-	    L13[ii] = 2.0*Lmu[3] - 2.0*timesI(Lmu[2]);
-	    ii++;
-    	}
-    }
-
-    auto ionameFn = [this](const unsigned int index, const unsigned int dummy)
-    {
-	//index = 8*l + 4*sindex1 + sindex2
-	unsigned int sindex = index % 8;
-	unsigned int sindex2 = sindex % 4;
- 	unsigned int sindex1 = (sindex - sindex2)/4;
-	unsigned int l = (index - sindex)/8;
-	return par().lepton[l] + "_" + std::to_string(sindex1) + std::to_string(sindex2);
-    };
-
-    auto filenameFn = [this, &ionameFn](const unsigned int index, const unsigned int dummy)
-    {
-        return par().output + "." + std::to_string(vm().getTrajectory()) 
-               + "/" + ionameFn(index, dummy) + ".h5";
-    };
-
-    auto metadataFn = [this](const unsigned int index, const unsigned int dummy)
-    {
-	A2ALeptonFieldMetadata md;
-
-	unsigned int sindex = index % 8;
-	unsigned int sindex2 = sindex % 4;
- 	unsigned int sindex1 = (sindex - sindex2)/4;
-	unsigned int l = (index - sindex)/8;
-        md.leptonName = par().lepton[l];
-	md.sidx1 = sindex1;
-	md.sidx2 = sindex2;
-
-	return md;
-    };
-
-    // executing computation
-    Kernel kernel(L02, L03, L12, L13, envGetGrid(FermionField));
-    envGetTmp(Computation, computation);
-    computation.execute(left, right, kernel, ionameFn, filenameFn, metadataFn);
-}
-
-END_MODULE_NAMESPACE
-
-END_HADRONS_NAMESPACE
-
-#endif // Hadrons_MContraction_A2ALeptonField_hpp_
diff --git a/Hadrons/modules.inc b/Hadrons/modules.inc
index b428b5eb..3b9e6563 100644
--- a/Hadrons/modules.inc
+++ b/Hadrons/modules.inc
@@ -7,7 +7,6 @@ modules_cc =\
   Modules/MContraction/WeakNonEye3pt.cc \
   Modules/MContraction/A2AAslashField.cc \
   Modules/MContraction/A2AMesonField.cc \
-  Modules/MContraction/A2ALeptonField.cc \
   Modules/MContraction/DiscLoop.cc \
   Modules/MContraction/Gamma3pt.cc \
   Modules/MFermion/FreeProp.cc \
@@ -73,7 +72,6 @@ modules_hpp =\
   Modules/MContraction/A2AAslashField.hpp \
   Modules/MContraction/A2ALoop.hpp \
   Modules/MContraction/A2AMesonField.hpp \
-  Modules/MContraction/A2ALeptonField.hpp \
   Modules/MContraction/Meson.hpp \
   Modules/MContraction/DiscLoop.hpp \
   Modules/MContraction/Gamma3pt.hpp \