diff --git a/tests/Test_meson_field.cc b/tests/Test_meson_field.cc
new file mode 100644
index 00000000..b108ee57
--- /dev/null
+++ b/tests/Test_meson_field.cc
@@ -0,0 +1,91 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: tests/core/Test_meson_field.cc
+
+Copyright (C) 2015-2018
+
+Author: Felix Erben <felix.erben@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
+*************************************************************************************/
+
+#include <Grid/Grid.h>
+#include <Grid/qcd/utils/A2Autils.h>
+
+using namespace Grid;
+
+typedef typename DomainWallFermionR::ComplexField ComplexField;
+typedef typename DomainWallFermionR::FermionField FermionField;
+
+int main(int argc, char *argv[])
+{
+  // initialization
+  Grid_init(&argc, &argv);
+  std::cout << GridLogMessage << "Grid initialized" << std::endl;
+ 
+  Coordinate latt_size   = GridDefaultLatt();
+  Coordinate simd_layout = GridDefaultSimd(4, vComplex::Nsimd());
+  Coordinate mpi_layout  = GridDefaultMpi();
+  GridCartesian    grid(latt_size,simd_layout,mpi_layout);
+
+  // MesonField lhs and rhs vectors
+  int mfDim = 10;
+  std::vector<FermionField> phi(mfDim,&grid);
+  std::vector<FermionField> rho(mfDim,&grid);
+  // Gamma matrices used in the contraction
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
+  };
+
+  // momentum phases e^{ipx}
+  std::vector<std::vector<double>> momenta = {
+	  {0.,0.,0.},
+	  {1.,0.,0.},
+	  {1.,1.,0.},
+	  {1.,1.,1.},
+	  {2.,0.,0.}
+  };
+  std::vector<ComplexField> phases(momenta.size(),&grid);
+  ComplexField coor(&grid);
+  Complex Ci(0.0,1.0);
+  for (unsigned int j = 0; j < momenta.size(); ++j)
+  {
+    phases[j] = Zero();
+    for(unsigned int mu = 0; mu < momenta[j].size(); mu++)
+    {
+      LatticeCoordinate(coor, mu);
+      phases[j] = phases[j] + momenta[j][mu]/GridDefaultLatt()[mu]*coor;
+    }
+    phases[j] = exp((Real)(2*M_PI)*Ci*phases[j]);
+  }
+
+  Eigen::Tensor<ComplexD,5, Eigen::RowMajor> mf;//(momenta.size(),12,GridDefaultLatt()[3],10,10);
+
+  //execute meson field routine
+  //A2Autils<WilsonImplR>::MesonField(mf,phi,phi,Gmu,phases,3);
+
+  // epilogue
+  std::cout << GridLogMessage << "Grid is finalizing now" << std::endl;
+  Grid_finalize();
+  
+  return EXIT_SUCCESS;
+}