Merge branch 'develop' into feature/distil

* develop: (29 commits)
  precision fix
  Updates after review with Peter.
  Wilson clover multi grid for lime lattice
  Recommendations for Traits classes
  Hadrons: uninitialised pointer fix (might have been harmless)
  Hadrons: beware of the nasty uninitialised twists
  Smearing test. Test on free field.
  Smearing for quark observables
  Smearing
  Hadrons: XML validator utility
  display relative norm during field IO norm check
  possibility to set a build number
  IO norm check on relative norm
  Output field norm check during IO
  Hadrons: random vector utility module I/O
  quieter initialisation
  fix patch command for eigen in bootstrap.sh
  Mres changes and gauge xform mat changes
  Hadrons: 32 bit I/O directly in Lanczos module
  Hadrons: copyright update
  ...

# Conflicts:
#	Grid/tensors/Tensor_traits.h
#	Hadrons/Modules.hpp
#	Hadrons/modules.inc

Hadrons/Modules/MAction/DWF.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MAction/DWF.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MAction/DWF.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MAction/DWF.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Lanny91 <andrew.lawson@gmail.com>
@@ -112,8 +112,6 @@ void TDWF<FImpl>::setup(void)
                  << par().mass << ", M5= " << par().M5 << " and Ls= "
                  << par().Ls << " using gauge field '" << par().gauge << "'"
                  << std::endl;
-    LOG(Message) << "Fermion boundary conditions: " << par().boundary
-                 << std::endl;
                  
     auto &U    = envGet(GaugeField, par().gauge);
     auto &g4   = *envGetGrid(FermionField);
@@ -121,8 +119,26 @@ void TDWF<FImpl>::setup(void)
     auto &g5   = *envGetGrid(FermionField, par().Ls);
     auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
     typename DomainWallFermion<FImpl>::ImplParams implParams;
-    implParams.boundary_phases = strToVec<Complex>(par().boundary);
-    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    if (!par().boundary.empty())
+    {
+        implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    }
+    if (!par().twist.empty())
+    {
+        implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    }
+    LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases
+                 << std::endl;
+    LOG(Message) << "Twists: " << implParams.twist_n_2pi_L
+                 << std::endl;
+    if (implParams.boundary_phases.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of boundary phase");
+    }
+    if (implParams.twist_n_2pi_L.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of twist");
+    }
     envCreateDerived(FMat, DomainWallFermion<FImpl>, getName(), par().Ls, U, g5,
                      grb5, g4, grb4, par().mass, par().M5, implParams);
 }

Hadrons/Modules/MAction/MobiusDWF.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MAction/MobiusDWF.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MAction/MobiusDWF.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MAction/MobiusDWF.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 
@@ -112,17 +112,33 @@ void TMobiusDWF<FImpl>::setup(void)
                  << ", b= " << par().b << ", c= " << par().c
                  << " using gauge field '" << par().gauge << "'"
                  << std::endl;
-    LOG(Message) << "Fermion boundary conditions: " << par().boundary
-                 << std::endl;
-
+                 
     auto &U    = envGet(GaugeField, par().gauge);
     auto &g4   = *envGetGrid(FermionField);
     auto &grb4 = *envGetRbGrid(FermionField);
     auto &g5   = *envGetGrid(FermionField, par().Ls);
     auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
     typename MobiusFermion<FImpl>::ImplParams implParams;
-    implParams.boundary_phases = strToVec<Complex>(par().boundary);
-    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    if (!par().boundary.empty())
+    {
+        implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    }
+    if (!par().twist.empty())
+    {
+        implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    }
+    LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases
+                 << std::endl;
+    LOG(Message) << "Twists: " << implParams.twist_n_2pi_L
+                 << std::endl;
+    if (implParams.boundary_phases.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of boundary phase");
+    }
+    if (implParams.twist_n_2pi_L.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of twist");
+    }
     envCreateDerived(FMat, MobiusFermion<FImpl>, getName(), par().Ls, U, g5,
                      grb5, g4, grb4, par().mass, par().M5, par().b, par().c,
                      implParams);

Hadrons/Modules/MAction/ScaledDWF.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MAction/ScaledDWF.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MAction/ScaledDWF.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MAction/ScaledDWF.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 
@@ -111,8 +111,6 @@ void TScaledDWF<FImpl>::setup(void)
                  << ", scale= " << par().scale
                  << " using gauge field '" << par().gauge << "'"
                  << std::endl;
-    LOG(Message) << "Fermion boundary conditions: " << par().boundary
-                 << std::endl;
 
     auto &U    = envGet(GaugeField, par().gauge);
     auto &g4   = *envGetGrid(FermionField);
@@ -120,8 +118,26 @@ void TScaledDWF<FImpl>::setup(void)
     auto &g5   = *envGetGrid(FermionField, par().Ls);
     auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
     typename ScaledShamirFermion<FImpl>::ImplParams implParams;
-    implParams.boundary_phases = strToVec<Complex>(par().boundary);
-    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    if (!par().boundary.empty())
+    {
+        implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    }
+    if (!par().twist.empty())
+    {
+        implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    }
+    LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases
+                 << std::endl;
+    LOG(Message) << "Twists: " << implParams.twist_n_2pi_L
+                 << std::endl;
+    if (implParams.boundary_phases.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of boundary phase");
+    }
+    if (implParams.twist_n_2pi_L.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of twist");
+    }
     envCreateDerived(FMat, ScaledShamirFermion<FImpl>, getName(), par().Ls, U, g5,
                      grb5, g4, grb4, par().mass, par().M5, par().scale,
                      implParams);

Hadrons/Modules/MAction/Wilson.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MAction/Wilson.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MAction/Wilson.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MAction/Wilson.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Lanny91 <andrew.lawson@gmail.com>
@@ -109,15 +109,31 @@ void TWilson<FImpl>::setup(void)
 {
     LOG(Message) << "Setting up Wilson fermion matrix with m= " << par().mass
                  << " using gauge field '" << par().gauge << "'" << std::endl;
-    LOG(Message) << "Fermion boundary conditions: " << par().boundary
-                 << std::endl;
                  
     auto &U      = envGet(GaugeField, par().gauge);
     auto &grid   = *envGetGrid(FermionField);
     auto &gridRb = *envGetRbGrid(FermionField);
     typename WilsonFermion<FImpl>::ImplParams implParams;
-    implParams.boundary_phases = strToVec<Complex>(par().boundary);
-    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    if (!par().boundary.empty())
+    {
+        implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    }
+    if (!par().twist.empty())
+    {
+        implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    }
+    LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases
+                 << std::endl;
+    LOG(Message) << "Twists: " << implParams.twist_n_2pi_L
+                 << std::endl;
+    if (implParams.boundary_phases.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of boundary phase");
+    }
+    if (implParams.twist_n_2pi_L.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of twist");
+    }
     envCreateDerived(FMat, WilsonFermion<FImpl>, getName(), 1, U, grid, gridRb,
                      par().mass, implParams);
 }

Hadrons/Modules/MAction/WilsonClover.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MAction/WilsonClover.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MAction/WilsonClover.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MAction/WilsonClover.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Guido Cossu <guido.cossu@ed.ac.uk>
@@ -112,17 +112,34 @@ void TWilsonClover<FImpl>::setup(void)
 {
     LOG(Message) << "Setting up Wilson clover fermion matrix with m= " << par().mass
                  << " using gauge field '" << par().gauge << "'" << std::endl;
-    LOG(Message) << "Fermion boundary conditions: " << par().boundary 
-                 << std::endl;
     LOG(Message) << "Clover term csw_r: " << par().csw_r
                  << " csw_t: " << par().csw_t
                  << std::endl;
+                 
     auto &U      = envGet(GaugeField, par().gauge);
     auto &grid   = *envGetGrid(FermionField);
     auto &gridRb = *envGetRbGrid(FermionField);
     typename WilsonCloverFermion<FImpl>::ImplParams implParams;
-    implParams.boundary_phases = strToVec<Complex>(par().boundary);
-    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    if (!par().boundary.empty())
+    {
+        implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    }
+    if (!par().twist.empty())
+    {
+        implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    }
+    LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases
+                 << std::endl;
+    LOG(Message) << "Twists: " << implParams.twist_n_2pi_L
+                 << std::endl;
+    if (implParams.boundary_phases.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of boundary phase");
+    }
+    if (implParams.twist_n_2pi_L.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of twist");
+    }
     envCreateDerived(FMat, WilsonCloverFermion<FImpl>, getName(), 1, U, grid,
                      gridRb, par().mass, par().csw_r, par().csw_t, 
                      par().clover_anisotropy, implParams); 

Hadrons/Modules/MAction/ZMobiusDWF.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MAction/ZMobiusDWF.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MAction/ZMobiusDWF.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MAction/ZMobiusDWF.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 
@@ -118,10 +118,7 @@ void TZMobiusDWF<FImpl>::setup(void)
     {
         LOG(Message) << "  omega[" << i << "]= " << par().omega[i] << std::endl;
     }
-    LOG(Message) << "Fermion boundary conditions: " << par().boundary
-                 << std::endl;
 
-    env().createGrid(par().Ls);
     auto &U    = envGet(GaugeField, par().gauge);
     auto &g4   = *envGetGrid(FermionField);
     auto &grb4 = *envGetRbGrid(FermionField);
@@ -129,8 +126,26 @@ void TZMobiusDWF<FImpl>::setup(void)
     auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
     auto omega = par().omega;
     typename ZMobiusFermion<FImpl>::ImplParams implParams;
-    implParams.boundary_phases = strToVec<Complex>(par().boundary);
-    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    if (!par().boundary.empty())
+    {
+        implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    }
+    if (!par().twist.empty())
+    {
+        implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    }
+    LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases
+                 << std::endl;
+    LOG(Message) << "Twists: " << implParams.twist_n_2pi_L
+                 << std::endl;
+    if (implParams.boundary_phases.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of boundary phase");
+    }
+    if (implParams.twist_n_2pi_L.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of twist");
+    }
     envCreateDerived(FMat, ZMobiusFermion<FImpl>, getName(), par().Ls, U, g5,
                      grb5, g4, grb4, par().mass, par().M5, omega,
                      par().b, par().c, implParams);

Hadrons/Modules/MContraction/A2AAslashField.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MContraction/A2AAslashField.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MContraction/A2AAslashField.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MContraction/A2AAslashField.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MContraction/A2ALoop.cc

@@ -2,9 +2,9 @@
 
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: Hadrons/Modules/MContraction/WardIdentity.cc
+Source file: Hadrons/Modules/MContraction/A2ALoop.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 
@@ -25,11 +25,10 @@ with this program; if not, write to the Free Software Foundation, Inc.,
 See the full license in the file "LICENSE" in the top level distribution directory
 *************************************************************************************/
 /*  END LEGAL */
-#include <Hadrons/Modules/MContraction/WardIdentity.hpp>
+#include <Hadrons/Modules/MContraction/A2ALoop.hpp>
 
 using namespace Grid;
 using namespace Hadrons;
 using namespace MContraction;
 
-template class Grid::Hadrons::MContraction::TWardIdentity<FIMPL>;
-
+template class Grid::Hadrons::MContraction::TA2ALoop<FIMPL>;

Hadrons/Modules/MContraction/A2ALoop.hpp

@@ -2,12 +2,11 @@
 
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: Hadrons/Modules/MLoop/NoiseLoop.hpp
+Source file: Hadrons/Modules/MContraction/A2ALoop.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
-Author: Lanny91 <andrew.lawson@gmail.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
@@ -26,9 +25,8 @@ with this program; if not, write to the Free Software Foundation, Inc.,
 See the full license in the file "LICENSE" in the top level distribution directory
 *************************************************************************************/
 /*  END LEGAL */
-
-#ifndef Hadrons_MLoop_NoiseLoop_hpp_
-#define Hadrons_MLoop_NoiseLoop_hpp_
+#ifndef Hadrons_MContraction_A2ALoop_hpp_
+#define Hadrons_MContraction_A2ALoop_hpp_
 
 #include <Hadrons/Global.hpp>
 #include <Hadrons/Module.hpp>
@@ -36,99 +34,90 @@ See the full license in the file "LICENSE" in the top level distribution directo
 
 BEGIN_HADRONS_NAMESPACE
 
-/*
- 
- Noise loop propagator
- -----------------------------
- * loop_x = q_x * adj(eta_x)
- 
- * options:
- - q = Result of inversion on noise source.
- - eta = noise source.
-
- */
-
-
 /******************************************************************************
- *                         NoiseLoop                                          *
+ *                  From closed loop from all-to-all vectors                  *
  ******************************************************************************/
-BEGIN_MODULE_NAMESPACE(MLoop)
+BEGIN_MODULE_NAMESPACE(MContraction)
 
-class NoiseLoopPar: Serializable
+class A2ALoopPar: Serializable
 {
 public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(NoiseLoopPar,
-                                    std::string, q,
-                                    std::string, eta);
+    GRID_SERIALIZABLE_CLASS_MEMBERS(A2ALoopPar,
+                                    std::string, left,
+                                    std::string, right);
 };
 
 template <typename FImpl>
-class TNoiseLoop: public Module<NoiseLoopPar>
+class TA2ALoop: public Module<A2ALoopPar>
 {
 public:
     FERM_TYPE_ALIASES(FImpl,);
 public:
     // constructor
-    TNoiseLoop(const std::string name);
+    TA2ALoop(const std::string name);
     // destructor
-    virtual ~TNoiseLoop(void) {};
+    virtual ~TA2ALoop(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
-protected:
     // setup
     virtual void setup(void);
     // execution
     virtual void execute(void);
 };
 
-MODULE_REGISTER_TMP(NoiseLoop, TNoiseLoop<FIMPL>, MLoop);
+MODULE_REGISTER_TMP(A2ALoop, TA2ALoop<FIMPL>, MContraction);
 
 /******************************************************************************
- *                 TNoiseLoop implementation                                  *
+ *                        TA2ALoop implementation                             *
  ******************************************************************************/
 // constructor /////////////////////////////////////////////////////////////////
 template <typename FImpl>
-TNoiseLoop<FImpl>::TNoiseLoop(const std::string name)
-: Module<NoiseLoopPar>(name)
+TA2ALoop<FImpl>::TA2ALoop(const std::string name)
+: Module<A2ALoopPar>(name)
 {}
 
 // dependencies/products ///////////////////////////////////////////////////////
 template <typename FImpl>
-std::vector<std::string> TNoiseLoop<FImpl>::getInput(void)
+std::vector<std::string> TA2ALoop<FImpl>::getInput(void)
 {
-    std::vector<std::string> in = {par().q, par().eta};
+    std::vector<std::string> in = {par().left, par().right};
     
     return in;
 }
 
 template <typename FImpl>
-std::vector<std::string> TNoiseLoop<FImpl>::getOutput(void)
+std::vector<std::string> TA2ALoop<FImpl>::getOutput(void)
 {
-    std::vector<std::string> out = {getName()};
+    std::vector<std::string> out = {};
     
     return out;
 }
 
 // setup ///////////////////////////////////////////////////////////////////////
 template <typename FImpl>
-void TNoiseLoop<FImpl>::setup(void)
+void TA2ALoop<FImpl>::setup(void)
 {
     envCreateLat(PropagatorField, getName());
 }
 
 // execution ///////////////////////////////////////////////////////////////////
 template <typename FImpl>
-void TNoiseLoop<FImpl>::execute(void)
+void TA2ALoop<FImpl>::execute(void)
 {
-    auto &loop = envGet(PropagatorField, getName());
-    auto &q    = envGet(PropagatorField, par().q);
-    auto &eta  = envGet(PropagatorField, par().eta);
-    loop = q*adj(eta);
+    auto &loop  = envGet(PropagatorField, getName());
+    auto &left  = envGet(std::vector<FermionField>, par().left);
+    auto &right = envGet(std::vector<FermionField>, par().right);
+
+    loop = zero;
+    for (unsigned int i = 0; i < left.size(); ++i)
+    {
+        loop += outerProduct(adj(left[i]), right[i]);
+    }
 }
 
 END_MODULE_NAMESPACE
 
 END_HADRONS_NAMESPACE
 
-#endif // Hadrons_MLoop_NoiseLoop_hpp_
+#endif // Hadrons_MContraction_A2ALoop_hpp_

Hadrons/Modules/MContraction/A2AMesonField.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MContraction/A2AMesonField.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: paboyle <paboyle@ph.ed.ac.uk>

Hadrons/Modules/MContraction/A2AMesonField.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MContraction/A2AMesonField.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>

Hadrons/Modules/MContraction/Baryon.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MContraction/Baryon.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MContraction/Baryon.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MContraction/Baryon.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Lanny91 <andrew.lawson@gmail.com>

Hadrons/Modules/MContraction/DiscLoop.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MContraction/DiscLoop.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MContraction/DiscLoop.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MContraction/DiscLoop.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Lanny91 <andrew.lawson@gmail.com>

Hadrons/Modules/MContraction/Gamma3pt.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MContraction/Gamma3pt.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MContraction/Gamma3pt.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MContraction/Gamma3pt.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Lanny91 <andrew.lawson@gmail.com>

Hadrons/Modules/MContraction/Meson.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MContraction/Meson.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MContraction/Meson.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MContraction/Meson.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Lanny91 <andrew.lawson@gmail.com>

Hadrons/Modules/MContraction/WardIdentity.hpp

@@ -1,224 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: Hadrons/Modules/MContraction/WardIdentity.hpp
-
-Copyright (C) 2015-2018
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-Author: Lanny91 <andrew.lawson@gmail.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_MContraction_WardIdentity_hpp_
-#define Hadrons_MContraction_WardIdentity_hpp_
-
-#include <Hadrons/Global.hpp>
-#include <Hadrons/Module.hpp>
-#include <Hadrons/ModuleFactory.hpp>
-
-BEGIN_HADRONS_NAMESPACE
-
-/*
-  Ward Identity contractions
- -----------------------------
- 
- * options:
- - q:          propagator, 5D if available (string)
- - action:     action module used for propagator solution (string)
- - mass:       mass of quark (double)
- - test_axial: whether or not to test PCAC relation.
-*/
-
-/******************************************************************************
- *                              WardIdentity                                  *
- ******************************************************************************/
-BEGIN_MODULE_NAMESPACE(MContraction)
-
-class WardIdentityPar: Serializable
-{
-public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(WardIdentityPar,
-                                    std::string, q,
-                                    std::string, action,
-                                    double,      mass,
-                                    bool,        test_axial);
-};
-
-template <typename FImpl>
-class TWardIdentity: public Module<WardIdentityPar>
-{
-public:
-    FERM_TYPE_ALIASES(FImpl,);
-public:
-    // constructor
-    TWardIdentity(const std::string name);
-    // destructor
-    virtual ~TWardIdentity(void) {};
-    // dependency relation
-    virtual std::vector<std::string> getInput(void);
-    virtual std::vector<std::string> getOutput(void);
-protected:
-    // setup
-    virtual void setup(void);
-    // execution
-    virtual void execute(void);
-private:
-    unsigned int Ls_;
-};
-
-MODULE_REGISTER_TMP(WardIdentity, TWardIdentity<FIMPL>, MContraction);
-
-/******************************************************************************
- *                     TWardIdentity implementation                           *
- ******************************************************************************/
-// constructor /////////////////////////////////////////////////////////////////
-template <typename FImpl>
-TWardIdentity<FImpl>::TWardIdentity(const std::string name)
-: Module<WardIdentityPar>(name)
-{}
-
-// dependencies/products ///////////////////////////////////////////////////////
-template <typename FImpl>
-std::vector<std::string> TWardIdentity<FImpl>::getInput(void)
-{
-    std::vector<std::string> in = {par().q, par().action};
-    
-    return in;
-}
-
-template <typename FImpl>
-std::vector<std::string> TWardIdentity<FImpl>::getOutput(void)
-{
-    std::vector<std::string> out = {};
-    
-    return out;
-}
-
-// setup ///////////////////////////////////////////////////////////////////////
-template <typename FImpl>
-void TWardIdentity<FImpl>::setup(void)
-{
-    Ls_ = env().getObjectLs(par().q);
-    if (Ls_ != env().getObjectLs(par().action))
-    {
-        HADRONS_ERROR(Size, "Ls mismatch between quark action and propagator");
-    }
-    envTmpLat(PropagatorField, "tmp");
-    envTmpLat(PropagatorField, "vector_WI");
-    if (par().test_axial)
-    {
-        envTmpLat(PropagatorField, "psi");
-        envTmpLat(LatticeComplex,  "PP");
-        envTmpLat(LatticeComplex,  "axial_defect");
-        envTmpLat(LatticeComplex,  "PJ5q");
-    }
-}
-
-// execution ///////////////////////////////////////////////////////////////////
-template <typename FImpl>
-void TWardIdentity<FImpl>::execute(void)
-{
-    LOG(Message) << "Performing Ward Identity checks for quark '" << par().q
-                 << "'." << std::endl;
-
-    auto  &q   = envGet(PropagatorField, par().q);
-    auto  &act = envGet(FMat, par().action);
-    Gamma g5(Gamma::Algebra::Gamma5);
-
-    // Compute D_mu V_mu, D here is backward derivative.
-    envGetTmp(PropagatorField, tmp);
-    envGetTmp(PropagatorField, vector_WI);
-    vector_WI    = zero;
-    for (unsigned int mu = 0; mu < Nd; ++mu)
-    {
-        act.ContractConservedCurrent(q, q, tmp, Current::Vector, mu);
-        tmp -= Cshift(tmp, mu, -1);
-        vector_WI += tmp;
-    }
-
-    // Test ward identity D_mu V_mu = 0;
-    LOG(Message) << "Vector Ward Identity check Delta_mu V_mu = " 
-                 << norm2(vector_WI) << std::endl;
-
-    if (par().test_axial)
-    {
-        envGetTmp(PropagatorField, psi);
-        envGetTmp(LatticeComplex, PP);
-        envGetTmp(LatticeComplex, axial_defect);
-        envGetTmp(LatticeComplex, PJ5q);
-        std::vector<TComplex> axial_buf;
-
-        // Compute <P|D_mu A_mu>, D is backwards derivative.
-        axial_defect = zero;
-        for (unsigned int mu = 0; mu < Nd; ++mu)
-        {
-            act.ContractConservedCurrent(q, q, tmp, Current::Axial, mu);
-            tmp -= Cshift(tmp, mu, -1);
-            axial_defect += trace(g5*tmp);
-        }
-
-        // Get <P|J5q> for 5D (zero for 4D) and <P|P>.
-        PJ5q = zero;
-        if (Ls_ > 1)
-        {
-            // <P|P>
-            ExtractSlice(tmp, q, 0, 0);
-            psi  = 0.5 * (tmp - g5*tmp);
-            ExtractSlice(tmp, q, Ls_ - 1, 0);
-            psi += 0.5 * (tmp + g5*tmp);
-            PP = trace(adj(psi)*psi);
-
-            // <P|5Jq>
-            ExtractSlice(tmp, q, Ls_/2 - 1, 0);
-            psi  = 0.5 * (tmp + g5*tmp);
-            ExtractSlice(tmp, q, Ls_/2, 0);
-            psi += 0.5 * (tmp - g5*tmp);
-            PJ5q = trace(adj(psi)*psi);
-        }
-        else
-        {
-            PP = trace(adj(q)*q);
-        }
-
-        // Test ward identity <P|D_mu A_mu> = 2m<P|P> + 2<P|J5q>
-        LOG(Message) << "|D_mu A_mu|^2 = " << norm2(axial_defect) << std::endl;
-        LOG(Message) << "|PP|^2        = " << norm2(PP) << std::endl;
-        LOG(Message) << "|PJ5q|^2      = " << norm2(PJ5q) << std::endl;
-        LOG(Message) << "Axial Ward Identity defect Delta_mu A_mu = "
-                     << norm2(axial_defect) << std::endl;
-    
-        // Axial defect by timeslice.
-        axial_defect -= 2.*(par().mass*PP + PJ5q);
-        LOG(Message) << "Check Axial defect by timeslice" << std::endl;
-        sliceSum(axial_defect, axial_buf, Tp);
-        for (int t = 0; t < axial_buf.size(); ++t)
-        {
-            LOG(Message) << "t = " << t << ": " 
-                         << TensorRemove(axial_buf[t]) << std::endl;
-        }
-    }
-}
-
-END_MODULE_NAMESPACE
-
-END_HADRONS_NAMESPACE
-
-#endif // Hadrons_WardIdentity_hpp_

Hadrons/Modules/MContraction/WeakEye3pt.cc

@@ -2,9 +2,9 @@
 
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: Hadrons/Modules/MUtilities/TestSeqGamma.cc
+Source file: Hadrons/Modules/MContraction/WeakEye3pt.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 
@@ -25,11 +25,10 @@ with this program; if not, write to the Free Software Foundation, Inc.,
 See the full license in the file "LICENSE" in the top level distribution directory
 *************************************************************************************/
 /*  END LEGAL */
-#include <Hadrons/Modules/MUtilities/TestSeqGamma.hpp>
+#include <Hadrons/Modules/MContraction/WeakEye3pt.hpp>
 
 using namespace Grid;
 using namespace Hadrons;
-using namespace MUtilities;
-
-template class Grid::Hadrons::MUtilities::TTestSeqGamma<FIMPL>;
+using namespace MContraction;
 
+template class Grid::Hadrons::MContraction::TWeakEye3pt<FIMPL>;

Hadrons/Modules/MContraction/WeakEye3pt.hpp

@@ -0,0 +1,200 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MContraction/WeakEye3pt.hpp
+
+Copyright (C) 2015-2019
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+Author: Lanny91 <andrew.lawson@gmail.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_MContraction_WeakEye3pt_hpp_
+#define Hadrons_MContraction_WeakEye3pt_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/*
+ * Weak Hamiltonian meson 3-pt diagrams, eye topologies.
+ * 
+ * Schematics:       loop                 |                  
+ *                  /-<-¬                 |                             
+ *                 /     \                |            qbl     G     qbr
+ *                 \     /                |        /----<------*------<----¬         
+ *            qbl   \   /    qbr          |       /          /-*-¬          \
+ *       /-----<-----* *-----<----¬       |      /          /  G  \          \
+ *  gIn *            G G           * gOut | gIn *           \     /  loop    * gOut
+ *       \                        /       |      \           \->-/          /   
+ *        \                      /        |       \                        /       
+ *         \---------->---------/         |        \----------->----------/        
+ *                   qs                   |                   qs                  
+ *                                        |
+ *                one trace               |                two traces
+ * 
+ * one trace : tr(qbr*gOut*qs*adj(gIn)*g5*adj(qbl)*g5*G*loop*G*qbr*gOut)
+ * two traces: tr(qbr*gOut*qs*adj(gIn)*g5*adj(qbl)*g5*G)*tr(loop*G)
+ * 
+ */
+
+BEGIN_MODULE_NAMESPACE(MContraction)
+
+class WeakEye3ptPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(WeakEye3ptPar,
+                                    std::string,    qBarLeft,
+                                    std::string,    qBarRight,
+                                    std::string,    qSpectator,
+                                    std::string,    loop,
+                                    unsigned int,   tOut,
+                                    Gamma::Algebra, gammaIn,
+                                    Gamma::Algebra, gammaOut,
+                                    std::string,    output);
+};
+
+template <typename FImpl>
+class TWeakEye3pt: public Module<WeakEye3ptPar>
+{
+public:
+    FERM_TYPE_ALIASES(FImpl,);
+    class Metadata: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Metadata,
+                                        Gamma::Algebra, in,
+                                        Gamma::Algebra, out,
+                                        Gamma::Algebra, op,
+                                        unsigned int,   trace);
+    };
+    typedef Correlator<Metadata> Result;
+public:
+    // constructor
+    TWeakEye3pt(const std::string name);
+    // destructor
+    virtual ~TWeakEye3pt(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(WeakEye3pt, TWeakEye3pt<FIMPL>, MContraction);
+
+/******************************************************************************
+ *                        TWeakEye3pt implementation                          *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TWeakEye3pt<FImpl>::TWeakEye3pt(const std::string name)
+: Module<WeakEye3ptPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TWeakEye3pt<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().qBarLeft, par().qBarRight, 
+                                   par().qSpectator};
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TWeakEye3pt<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TWeakEye3pt<FImpl>::setup(void)
+{
+    envTmpLat(ComplexField, "corr");
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TWeakEye3pt<FImpl>::execute(void)
+{
+    LOG(Message) << "Computing mesonic weak 3pt contractions, eye topologies" << std::endl;
+    LOG(Message) << "gIn : " << par().gammaIn << std::endl;
+    LOG(Message) << "gOut: " << par().gammaIn << std::endl;
+    LOG(Message) << "tOut: " << par().tOut << std::endl;
+    LOG(Message) << "qbl : " << par().qBarLeft << std::endl;
+    LOG(Message) << "qbr : " << par().qBarRight << std::endl;
+    LOG(Message) << "qs  : " << par().qSpectator << std::endl;
+    LOG(Message) << "loop: " << par().loop << std::endl;
+
+    std::vector<Result> result;
+    Result              r;
+    auto                &qbl  = envGet(PropagatorField, par().qBarLeft);
+    auto                &qbr  = envGet(PropagatorField, par().qBarRight);
+    auto                &loop = envGet(PropagatorField, par().loop);
+    auto                &qs   = envGet(SlicedPropagator, par().qSpectator);
+    auto                qst   = qs[par().tOut];
+    Gamma               gIn(par().gammaIn), gOut(par().gammaOut);
+    Gamma               g5(Gamma::Algebra::Gamma5);
+
+    envGetTmp(ComplexField, corr);
+    r.info.in  = par().gammaIn;
+    r.info.out = par().gammaOut;
+    for (auto &G: Gamma::gall)
+    {
+        SlicedComplex buf;
+
+        r.info.op = G.g;
+        // one trace
+        corr = trace(qbr*gOut*qst*adj(gIn)*g5*adj(qbl)*g5*G*loop*G*qbr*gOut);
+        sliceSum(corr, buf, Tp);
+        r.corr.clear();
+        for (unsigned int t = 0; t < buf.size(); ++t)
+        {
+            r.corr.push_back(TensorRemove(buf[t]));
+        }
+        r.info.trace = 1;
+        result.push_back(r);
+        // two traces
+        corr = trace(qbr*gOut*qst*adj(gIn)*g5*adj(qbl)*g5*G)*trace(loop*G);
+        sliceSum(corr, buf, Tp);
+        r.corr.clear();
+        for (unsigned int t = 0; t < buf.size(); ++t)
+        {
+            r.corr.push_back(TensorRemove(buf[t]));
+        }
+        r.info.trace = 2;
+        result.push_back(r);
+    }
+    saveResult(par().output, "weakEye3pt", result);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MContraction_WeakEye3pt_hpp_

Hadrons/Modules/MContraction/WeakHamiltonian.hpp

@@ -1,118 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: Hadrons/Modules/MContraction/WeakHamiltonian.hpp
-
-Copyright (C) 2015-2018
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-Author: Lanny91 <andrew.lawson@gmail.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_MContraction_WeakHamiltonian_hpp_
-#define Hadrons_MContraction_WeakHamiltonian_hpp_
-
-#include <Hadrons/Global.hpp>
-#include <Hadrons/Module.hpp>
-#include <Hadrons/ModuleFactory.hpp>
-
-BEGIN_HADRONS_NAMESPACE
-
-/******************************************************************************
- *                         WeakHamiltonian                                    *
- ******************************************************************************/
-BEGIN_MODULE_NAMESPACE(MContraction)
-
-/*******************************************************************************
- * Utilities for contractions involving the Weak Hamiltonian.
- ******************************************************************************/
-//// Sum and store correlator.
-#define MAKE_DIAG(exp, buf, res, n)\
-sliceSum(exp, buf, Tp);\
-res.name = (n);\
-res.corr.resize(buf.size());\
-for (unsigned int t = 0; t < buf.size(); ++t)\
-{\
-    res.corr[t] = TensorRemove(buf[t]);\
-}
-
-//// Contraction of mu index: use 'mu' variable in exp.
-#define SUM_MU(buf,exp)\
-buf = zero;\
-for (unsigned int mu = 0; mu < ndim; ++mu)\
-{\
-    buf += exp;\
-}
-
-enum 
-{
-  i_V = 0,
-  i_A = 1,
-  n_i = 2
-};
-
-class WeakHamiltonianPar: Serializable
-{
-public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(WeakHamiltonianPar,
-                                    std::string, q1,
-                                    std::string, q2,
-                                    std::string, q3,
-                                    std::string, q4,
-                                    unsigned int, tSnk,
-                                    std::string, output);
-};
-
-#define MAKE_WEAK_MODULE(modname)\
-class T##modname: public Module<WeakHamiltonianPar>\
-{\
-public:\
-    FERM_TYPE_ALIASES(FIMPL,)\
-    class Result: Serializable\
-    {\
-    public:\
-        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,\
-                                        std::string, name,\
-                                        std::vector<Complex>, corr);\
-    };\
-public:\
-    /* constructor */ \
-    T##modname(const std::string name);\
-    /* destructor */ \
-    virtual ~T##modname(void) {};\
-    /* dependency relation */ \
-    virtual std::vector<std::string> getInput(void);\
-    virtual std::vector<std::string> getOutput(void);\
-public:\
-    std::vector<std::string> VA_label = {"V", "A"};\
-protected:\
-    /* setup */ \
-    virtual void setup(void);\
-    /* execution */ \
-    virtual void execute(void);\
-};\
-MODULE_REGISTER(modname, T##modname, MContraction);
-
-END_MODULE_NAMESPACE
-
-END_HADRONS_NAMESPACE
-
-#endif // Hadrons_MContraction_WeakHamiltonian_hpp_

Hadrons/Modules/MContraction/WeakHamiltonianEye.cc

@@ -1,151 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: Hadrons/Modules/MContraction/WeakHamiltonianEye.cc
-
-Copyright (C) 2015-2018
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-Author: Lanny91 <andrew.lawson@gmail.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 <Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp>
-
-using namespace Grid;
-using namespace Hadrons;
-using namespace MContraction;
-
-/*
- * Weak Hamiltonian current-current contractions, Eye-type.
- * 
- * These contractions are generated by the Q1 and Q2 operators in the physical
- * basis (see e.g. Fig 3 of arXiv:1507.03094).
- * 
- * Schematics:        q4                 |                  
- *                  /-<-¬                |                             
- *                 /     \               |             q2           q3
- *                 \     /               |        /----<------*------<----¬                        
- *            q2    \   /    q3          |       /          /-*-¬          \
- *       /-----<-----* *-----<----¬      |      /          /     \          \
- *    i *            H_W           * f   |   i *           \     /  q4      * f
- *       \                        /      |      \           \->-/          /   
- *        \                      /       |       \                        /       
- *         \---------->---------/        |        \----------->----------/        
- *                   q1                  |                   q1                  
- *                                       |
- *                Saucer (S)             |                  Eye (E)
- * 
- * S: trace(q3*g5*q1*adj(q2)*g5*gL[mu][p_1]*q4*gL[mu][p_2])
- * E: trace(q3*g5*q1*adj(q2)*g5*gL[mu][p_1])*trace(q4*gL[mu][p_2])
- * 
- * Note q1 must be sink smeared.
- */
-
-/******************************************************************************
- *                  TWeakHamiltonianEye implementation                        *
- ******************************************************************************/
-// constructor /////////////////////////////////////////////////////////////////
-TWeakHamiltonianEye::TWeakHamiltonianEye(const std::string name)
-: Module<WeakHamiltonianPar>(name)
-{}
-
-// dependencies/products ///////////////////////////////////////////////////////
-std::vector<std::string> TWeakHamiltonianEye::getInput(void)
-{
-    std::vector<std::string> in = {par().q1, par().q2, par().q3, par().q4};
-    
-    return in;
-}
-
-std::vector<std::string> TWeakHamiltonianEye::getOutput(void)
-{
-    std::vector<std::string> out = {};
-    
-    return out;
-}
-
-// setup ///////////////////////////////////////////////////////////////////////
-void TWeakHamiltonianEye::setup(void)
-{
-    unsigned int ndim = env().getNd();
-
-    envTmpLat(LatticeComplex,  "expbuf");
-    envTmpLat(PropagatorField, "tmp1");
-    envTmpLat(LatticeComplex,  "tmp2");
-    envTmp(std::vector<PropagatorField>, "S_body", 1, ndim, PropagatorField(env().getGrid()));
-    envTmp(std::vector<PropagatorField>, "S_loop", 1, ndim, PropagatorField(env().getGrid()));
-    envTmp(std::vector<LatticeComplex>,  "E_body", 1, ndim, LatticeComplex(env().getGrid()));
-    envTmp(std::vector<LatticeComplex>,  "E_loop", 1, ndim, LatticeComplex(env().getGrid()));
-}
-
-// execution ///////////////////////////////////////////////////////////////////
-void TWeakHamiltonianEye::execute(void)
-{
-    LOG(Message) << "Computing Weak Hamiltonian (Eye type) contractions '" 
-                 << getName() << "' using quarks '" << par().q1 << "', '" 
-                 << par().q2 << ", '" << par().q3 << "' and '" << par().q4 
-                 << "'." << std::endl;
-
-    auto                   &q1 = envGet(SlicedPropagator, par().q1);
-    auto                   &q2 = envGet(PropagatorField, par().q2);
-    auto                   &q3 = envGet(PropagatorField, par().q3);
-    auto                   &q4 = envGet(PropagatorField, par().q4);
-    Gamma                  g5  = Gamma(Gamma::Algebra::Gamma5);
-    std::vector<TComplex>  corrbuf;
-    std::vector<Result>    result(n_eye_diag);
-    unsigned int ndim    = env().getNd();
-
-    envGetTmp(LatticeComplex,               expbuf); 
-    envGetTmp(PropagatorField,              tmp1);
-    envGetTmp(LatticeComplex,               tmp2);
-    envGetTmp(std::vector<PropagatorField>, S_body);
-    envGetTmp(std::vector<PropagatorField>, S_loop);
-    envGetTmp(std::vector<LatticeComplex>,  E_body);
-    envGetTmp(std::vector<LatticeComplex>,  E_loop);
-
-    // Get sink timeslice of q1.
-    SitePropagator q1Snk = q1[par().tSnk];
-
-    // Setup for S-type contractions.
-    for (int mu = 0; mu < ndim; ++mu)
-    {
-        S_body[mu] = MAKE_SE_BODY(q1Snk, q2, q3, GammaL(Gamma::gmu[mu]));
-        S_loop[mu] = MAKE_SE_LOOP(q4, GammaL(Gamma::gmu[mu]));
-    }
-
-    // Perform S-type contractions.    
-    SUM_MU(expbuf, trace(S_body[mu]*S_loop[mu]))
-    MAKE_DIAG(expbuf, corrbuf, result[S_diag], "HW_S")
-
-    // Recycle sub-expressions for E-type contractions.
-    for (unsigned int mu = 0; mu < ndim; ++mu)
-    {
-        E_body[mu] = trace(S_body[mu]);
-        E_loop[mu] = trace(S_loop[mu]);
-    }
-
-    // Perform E-type contractions.
-    SUM_MU(expbuf, E_body[mu]*E_loop[mu])
-    MAKE_DIAG(expbuf, corrbuf, result[E_diag], "HW_E")
-
-    // IO
-    saveResult(par().output, "HW_Eye", result);
-}

Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp

@@ -1,59 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp
-
-Copyright (C) 2015-2018
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-Author: Lanny91 <andrew.lawson@gmail.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_MContraction_WeakHamiltonianEye_hpp_
-#define Hadrons_MContraction_WeakHamiltonianEye_hpp_
-
-#include <Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
-
-BEGIN_HADRONS_NAMESPACE
-
-/******************************************************************************
- *                         WeakHamiltonianEye                                 *
- ******************************************************************************/
-BEGIN_MODULE_NAMESPACE(MContraction)
-
-enum
-{
-    S_diag = 0,
-    E_diag = 1,
-    n_eye_diag = 2
-};
-
-// Saucer and Eye subdiagram contractions.
-#define MAKE_SE_BODY(Q_1, Q_2, Q_3, gamma) (Q_3*g5*Q_1*adj(Q_2)*g5*gamma)
-#define MAKE_SE_LOOP(Q_loop, gamma) (Q_loop*gamma)
-
-MAKE_WEAK_MODULE(WeakHamiltonianEye)
-
-END_MODULE_NAMESPACE
-
-END_HADRONS_NAMESPACE
-
-#endif // Hadrons_MContraction_WeakHamiltonianEye_hpp_

Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc

@@ -1,148 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc
-
-Copyright (C) 2015-2018
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-Author: Lanny91 <andrew.lawson@gmail.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 <Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp>
-
-using namespace Grid;
-using namespace Hadrons;
-using namespace MContraction;
-
-/*
- * Weak Hamiltonian current-current contractions, Non-Eye-type.
- * 
- * These contractions are generated by the Q1 and Q2 operators in the physical
- * basis (see e.g. Fig 3 of arXiv:1507.03094).
- * 
- * Schematic:     
- *            q2             q3          |           q2              q3
- *          /--<--¬       /--<--¬        |        /--<--¬         /--<--¬       
- *         /       \     /       \       |       /       \       /       \      
- *        /         \   /         \      |      /         \     /         \     
- *       /           \ /           \     |     /           \   /           \    
- *    i *             * H_W         *  f |  i *             * * H_W         * f 
- *      \             *             |    |     \           /   \           /
- *       \           / \           /     |      \         /     \         /    
- *        \         /   \         /      |       \       /       \       /  
- *         \       /     \       /       |        \-->--/         \-->--/      
- *          \-->--/       \-->--/        |          q1               q4 
- *            q1             q4          |
- *                Connected (C)          |                 Wing (W)
- *
- * C: trace(q1*adj(q2)*g5*gL[mu]*q3*adj(q4)*g5*gL[mu])
- * W: trace(q1*adj(q2)*g5*gL[mu])*trace(q3*adj(q4)*g5*gL[mu])
- * 
- */
-
-/******************************************************************************
- *                  TWeakHamiltonianNonEye implementation                     *
- ******************************************************************************/
-// constructor /////////////////////////////////////////////////////////////////
-TWeakHamiltonianNonEye::TWeakHamiltonianNonEye(const std::string name)
-: Module<WeakHamiltonianPar>(name)
-{}
-
-// dependencies/products ///////////////////////////////////////////////////////
-std::vector<std::string> TWeakHamiltonianNonEye::getInput(void)
-{
-    std::vector<std::string> in = {par().q1, par().q2, par().q3, par().q4};
-    
-    return in;
-}
-
-std::vector<std::string> TWeakHamiltonianNonEye::getOutput(void)
-{
-    std::vector<std::string> out = {};
-    
-    return out;
-}
-
-// setup ///////////////////////////////////////////////////////////////////////
-void TWeakHamiltonianNonEye::setup(void)
-{
-    unsigned int ndim = env().getNd();
-
-    envTmpLat(LatticeComplex,  "expbuf");
-    envTmpLat(PropagatorField, "tmp1");
-    envTmpLat(LatticeComplex,  "tmp2");
-    envTmp(std::vector<PropagatorField>, "C_i_side_loop", 1, ndim, PropagatorField(env().getGrid()));
-    envTmp(std::vector<PropagatorField>, "C_f_side_loop", 1, ndim, PropagatorField(env().getGrid()));
-    envTmp(std::vector<LatticeComplex>,  "W_i_side_loop", 1, ndim, LatticeComplex(env().getGrid()));
-    envTmp(std::vector<LatticeComplex>,  "W_f_side_loop", 1, ndim, LatticeComplex(env().getGrid()));
-}
-
-// execution ///////////////////////////////////////////////////////////////////
-void TWeakHamiltonianNonEye::execute(void)
-{
-    LOG(Message) << "Computing Weak Hamiltonian (Non-Eye type) contractions '" 
-                 << getName() << "' using quarks '" << par().q1 << "', '" 
-                 << par().q2 << ", '" << par().q3 << "' and '" << par().q4 
-                 << "'." << std::endl;
-    
-    auto                  &q1 = envGet(PropagatorField, par().q1);
-    auto                  &q2 = envGet(PropagatorField, par().q2);
-    auto                  &q3 = envGet(PropagatorField, par().q3);
-    auto                  &q4 = envGet(PropagatorField, par().q4);
-    Gamma                 g5  = Gamma(Gamma::Algebra::Gamma5);
-    std::vector<TComplex> corrbuf;
-    std::vector<Result>   result(n_noneye_diag); 
-    unsigned int          ndim = env().getNd();
-
-    envGetTmp(LatticeComplex,               expbuf); 
-    envGetTmp(PropagatorField,              tmp1);
-    envGetTmp(LatticeComplex,               tmp2);
-    envGetTmp(std::vector<PropagatorField>, C_i_side_loop);
-    envGetTmp(std::vector<PropagatorField>, C_f_side_loop);
-    envGetTmp(std::vector<LatticeComplex>,  W_i_side_loop);
-    envGetTmp(std::vector<LatticeComplex>,  W_f_side_loop);
-
-    // Setup for C-type contractions.
-    for (int mu = 0; mu < ndim; ++mu)
-    {
-        C_i_side_loop[mu] = MAKE_CW_SUBDIAG(q1, q2, GammaL(Gamma::gmu[mu]));
-        C_f_side_loop[mu] = MAKE_CW_SUBDIAG(q3, q4, GammaL(Gamma::gmu[mu]));
-    }
-
-    // Perform C-type contractions.    
-    SUM_MU(expbuf, trace(C_i_side_loop[mu]*C_f_side_loop[mu]))
-    MAKE_DIAG(expbuf, corrbuf, result[C_diag], "HW_C")
-
-    // Recycle sub-expressions for W-type contractions.
-    for (unsigned int mu = 0; mu < ndim; ++mu)
-    {
-        W_i_side_loop[mu] = trace(C_i_side_loop[mu]);
-        W_f_side_loop[mu] = trace(C_f_side_loop[mu]);
-    }
-
-    // Perform W-type contractions.
-    SUM_MU(expbuf, W_i_side_loop[mu]*W_f_side_loop[mu])
-    MAKE_DIAG(expbuf, corrbuf, result[W_diag], "HW_W")
-
-    // IO
-    saveResult(par().output, "HW_NonEye", result);
-}

Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp

@@ -1,58 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp
-
-Copyright (C) 2015-2018
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-Author: Lanny91 <andrew.lawson@gmail.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_MContraction_WeakHamiltonianNonEye_hpp_
-#define Hadrons_MContraction_WeakHamiltonianNonEye_hpp_
-
-#include <Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
-
-BEGIN_HADRONS_NAMESPACE
-
-/******************************************************************************
- *                         WeakHamiltonianNonEye                              *
- ******************************************************************************/
-BEGIN_MODULE_NAMESPACE(MContraction)
-
-enum
-{
-    W_diag = 0,
-    C_diag = 1,
-    n_noneye_diag = 2
-};
-
-// Wing and Connected subdiagram contractions
-#define MAKE_CW_SUBDIAG(Q_1, Q_2, gamma) (Q_1*adj(Q_2)*g5*gamma)
-
-MAKE_WEAK_MODULE(WeakHamiltonianNonEye)
-
-END_MODULE_NAMESPACE
-
-END_HADRONS_NAMESPACE
-
-#endif // Hadrons_MContraction_WeakHamiltonianNonEye_hpp_

Hadrons/Modules/MContraction/WeakNeutral4ptDisc.cc

@@ -1,142 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: Hadrons/Modules/MContraction/WeakNeutral4ptDisc.cc
-
-Copyright (C) 2015-2018
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-Author: Lanny91 <andrew.lawson@gmail.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 <Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp>
-
-using namespace Grid;
-using namespace Hadrons;
-using namespace MContraction;
-
-/*
- * Weak Hamiltonian + current contractions, disconnected topology for neutral 
- * mesons.
- * 
- * These contractions are generated by operators Q_1,...,10 of the dS=1 Weak
- * Hamiltonian in the physical basis and an additional current J (see e.g. 
- * Fig 11 of arXiv:1507.03094).
- * 
- * Schematic:
- *                        
- *           q2          q4             q3
- *       /--<--¬     /---<--¬       /---<--¬
- *     /         \ /         \     /        \
- *  i *           * H_W      |  J *          * f
- *     \         / \         /     \        /
- *      \--->---/   \-------/       \------/
- *          q1 
- * 
- * options
- * - q1: input propagator 1 (string)
- * - q2: input propagator 2 (string)
- * - q3: input propagator 3 (string), assumed to be sequential propagator 
- * - q4: input propagator 4 (string), assumed to be a loop
- * 
- * type 1: trace(q1*adj(q2)*g5*gL[mu])*trace(loop*gL[mu])*trace(q3*g5)
- * type 2: trace(q1*adj(q2)*g5*gL[mu]*loop*gL[mu])*trace(q3*g5)
- */
-
-/*******************************************************************************
- *                  TWeakNeutral4ptDisc implementation                         *
- ******************************************************************************/
-// constructor /////////////////////////////////////////////////////////////////
-TWeakNeutral4ptDisc::TWeakNeutral4ptDisc(const std::string name)
-: Module<WeakHamiltonianPar>(name)
-{}
-
-// dependencies/products ///////////////////////////////////////////////////////
-std::vector<std::string> TWeakNeutral4ptDisc::getInput(void)
-{
-    std::vector<std::string> in = {par().q1, par().q2, par().q3, par().q4};
-    
-    return in;
-}
-
-std::vector<std::string> TWeakNeutral4ptDisc::getOutput(void)
-{
-    std::vector<std::string> out = {};
-    
-    return out;
-}
-
-// setup ///////////////////////////////////////////////////////////////////////
-void TWeakNeutral4ptDisc::setup(void)
-{
-    unsigned int ndim = env().getNd();
-
-    envTmpLat(LatticeComplex,  "expbuf");
-    envTmpLat(PropagatorField, "tmp");
-    envTmpLat(LatticeComplex,  "curr");
-    envTmp(std::vector<PropagatorField>, "meson", 1, ndim, PropagatorField(env().getGrid()));
-    envTmp(std::vector<PropagatorField>, "loop", 1, ndim,  PropagatorField(env().getGrid()));
-}
-
-// execution ///////////////////////////////////////////////////////////////////
-void TWeakNeutral4ptDisc::execute(void)
-{
-    LOG(Message) << "Computing Weak Hamiltonian neutral disconnected contractions '" 
-                 << getName() << "' using quarks '" << par().q1 << "', '" 
-                 << par().q2 << ", '" << par().q3 << "' and '" << par().q4 
-                 << "'." << std::endl;
-
-    auto                  &q1 = envGet(PropagatorField, par().q1);
-    auto                  &q2 = envGet(PropagatorField, par().q2);
-    auto                  &q3 = envGet(PropagatorField, par().q3);
-    auto                  &q4 = envGet(PropagatorField, par().q4);
-    Gamma                 g5  = Gamma(Gamma::Algebra::Gamma5);
-    std::vector<TComplex> corrbuf;
-    std::vector<Result>   result(n_neut_disc_diag);
-    unsigned int          ndim = env().getNd();
-
-    envGetTmp(LatticeComplex,               expbuf); 
-    envGetTmp(PropagatorField,              tmp);
-    envGetTmp(LatticeComplex,               curr);
-    envGetTmp(std::vector<PropagatorField>, meson);
-    envGetTmp(std::vector<PropagatorField>, loop);
-
-    // Setup for type 1 contractions.
-    for (int mu = 0; mu < ndim; ++mu)
-    {
-        meson[mu] = MAKE_DISC_MESON(q1, q2, GammaL(Gamma::gmu[mu]));
-        loop[mu] = MAKE_DISC_LOOP(q4, GammaL(Gamma::gmu[mu]));
-    }
-    curr = MAKE_DISC_CURR(q3, GammaL(Gamma::Algebra::Gamma5));
-
-    // Perform type 1 contractions.    
-    SUM_MU(expbuf, trace(meson[mu]*loop[mu]))
-    expbuf *= curr;
-    MAKE_DIAG(expbuf, corrbuf, result[neut_disc_1_diag], "HW_disc0_1")
-
-    // Perform type 2 contractions.
-    SUM_MU(expbuf, trace(meson[mu])*trace(loop[mu]))
-    expbuf *= curr;
-    MAKE_DIAG(expbuf, corrbuf, result[neut_disc_2_diag], "HW_disc0_2")
-
-    // IO
-    saveResult(par().output, "HW_disc0", result);
-}

Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp

@@ -1,60 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp
-
-Copyright (C) 2015-2018
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-Author: Lanny91 <andrew.lawson@gmail.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_MContraction_WeakNeutral4ptDisc_hpp_
-#define Hadrons_MContraction_WeakNeutral4ptDisc_hpp_
-
-#include <Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
-
-BEGIN_HADRONS_NAMESPACE
-
-/******************************************************************************
- *                         WeakNeutral4ptDisc                                 *
- ******************************************************************************/
-BEGIN_MODULE_NAMESPACE(MContraction)
-
-enum
-{
-    neut_disc_1_diag = 0,
-    neut_disc_2_diag = 1,
-    n_neut_disc_diag = 2
-};
-
-// Neutral 4pt disconnected subdiagram contractions.
-#define MAKE_DISC_MESON(Q_1, Q_2, gamma) (Q_1*adj(Q_2)*g5*gamma)
-#define MAKE_DISC_LOOP(Q_LOOP, gamma) (Q_LOOP*gamma)
-#define MAKE_DISC_CURR(Q_c, gamma) (trace(Q_c*gamma))
-
-MAKE_WEAK_MODULE(WeakNeutral4ptDisc)
-
-END_MODULE_NAMESPACE
-
-END_HADRONS_NAMESPACE
-
-#endif // Hadrons_MContraction_WeakNeutral4ptDisc_hpp_

Hadrons/Modules/MContraction/WeakNonEye3pt.cc

@@ -2,9 +2,9 @@
 
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: Hadrons/Modules/MLoop/NoiseLoop.cc
+Source file: Hadrons/Modules/MContraction/WeakNonEye3pt.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 
@@ -25,11 +25,10 @@ with this program; if not, write to the Free Software Foundation, Inc.,
 See the full license in the file "LICENSE" in the top level distribution directory
 *************************************************************************************/
 /*  END LEGAL */
-#include <Hadrons/Modules/MLoop/NoiseLoop.hpp>
+#include <Hadrons/Modules/MContraction/WeakNonEye3pt.hpp>
 
 using namespace Grid;
 using namespace Hadrons;
-using namespace MLoop;
-
-template class Grid::Hadrons::MLoop::TNoiseLoop<FIMPL>;
+using namespace MContraction;
 
+template class Grid::Hadrons::MContraction::TWeakNonEye3pt<FIMPL>;

Hadrons/Modules/MContraction/WeakNonEye3pt.hpp

@@ -0,0 +1,198 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MContraction/WeakNonEye3pt.hpp
+
+Copyright (C) 2015-2019
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+Author: Lanny91 <andrew.lawson@gmail.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_MContraction_WeakNonEye3pt_hpp_
+#define Hadrons_MContraction_WeakNonEye3pt_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/*
+ * Weak Hamiltonian meson 3-pt diagrams, non-eye topologies.
+ * 
+ * Schematic:     
+ *            qbl           qbr            |            qbl             qbr
+ *          /--<--¬       /--<--¬          |          /--<--¬         /--<--¬       
+ *         /       \     /       \         |         /       \       /       \      
+ *        /         \   /         \        |        /         \     /         \     
+ *       /           \ /           \       |       /           \   /           \    
+ *  gIn *             * G           * gOut |  gIn *           G * * G           * gOut
+ *      \             * G           |      |       \           /   \           /
+ *       \           / \           /       |        \         /     \         /    
+ *        \         /   \         /        |         \       /       \       /  
+ *         \       /     \       /         |          \-->--/         \-->--/      
+ *          \-->--/       \-->--/          |            ql              qr 
+ *            ql            qr             |
+ *               one trace                 |              two traces
+ *
+ * one trace : tr(ql*adj(gIn)*g5*adj(qbl)*g5*G*qbr*gOut*g5*adj(qr)*g5*G)
+ * two traces: tr(ql*adj(gIn)*g5*adj(qbl)*g5*G)*tr(qbr*gOut*g5*adj(qr)*g5*G)
+ * 
+ */
+
+BEGIN_MODULE_NAMESPACE(MContraction)
+
+class WeakNonEye3ptPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(WeakNonEye3ptPar,
+                                    std::string,    qLeft,
+                                    std::string,    qBarLeft,
+                                    std::string,    qRight,
+                                    std::string,    qBarRight,
+                                    Gamma::Algebra, gammaIn,
+                                    Gamma::Algebra, gammaOut,
+                                    std::string,    output);
+};
+
+template <typename FImpl>
+class TWeakNonEye3pt: public Module<WeakNonEye3ptPar>
+{
+public:
+    FERM_TYPE_ALIASES(FImpl,);
+    class Metadata: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Metadata,
+                                        Gamma::Algebra, in,
+                                        Gamma::Algebra, out,
+                                        Gamma::Algebra, op,
+                                        unsigned int,   trace);
+    };
+    typedef Correlator<Metadata> Result;
+public:
+    // constructor
+    TWeakNonEye3pt(const std::string name);
+    // destructor
+    virtual ~TWeakNonEye3pt(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(WeakNonEye3pt, TWeakNonEye3pt<FIMPL>, MContraction);
+
+/******************************************************************************
+ *                 TWeakNonEye3pt implementation                             *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TWeakNonEye3pt<FImpl>::TWeakNonEye3pt(const std::string name)
+: Module<WeakNonEye3ptPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TWeakNonEye3pt<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().qLeft, par().qBarLeft, 
+                                   par().qRight, par().qBarRight};
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TWeakNonEye3pt<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TWeakNonEye3pt<FImpl>::setup(void)
+{
+    envTmpLat(ComplexField, "corr");
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TWeakNonEye3pt<FImpl>::execute(void)
+{
+    LOG(Message) << "Computing mesonic weak 3pt contractions, non-eye topologies" << std::endl;
+    LOG(Message) << "gIn : " << par().gammaIn << std::endl;
+    LOG(Message) << "gOut: " << par().gammaIn << std::endl;
+    LOG(Message) << "ql  : " << par().qLeft << std::endl;
+    LOG(Message) << "qbl : " << par().qBarLeft << std::endl;
+    LOG(Message) << "qr  : " << par().qRight << std::endl;
+    LOG(Message) << "qbr : " << par().qBarRight << std::endl;
+
+    std::vector<Result> result;
+    Result              r;
+    auto                &ql  = envGet(PropagatorField, par().qLeft);
+    auto                &qbl = envGet(PropagatorField, par().qBarLeft);
+    auto                &qr  = envGet(PropagatorField, par().qRight);
+    auto                &qbr = envGet(PropagatorField, par().qBarRight);
+    Gamma               gIn(par().gammaIn), gOut(par().gammaOut);
+    Gamma               g5(Gamma::Algebra::Gamma5);
+
+    envGetTmp(ComplexField, corr);
+    r.info.in  = par().gammaIn;
+    r.info.out = par().gammaOut;
+    for (auto &G: Gamma::gall)
+    {
+        SlicedComplex buf;
+
+        r.info.op = G.g;
+        // one trace
+        corr = trace(ql*adj(gIn)*g5*adj(qbl)*g5*G*qbr*gOut*g5*adj(qr)*g5*G);
+        sliceSum(corr, buf, Tp);
+        r.corr.clear();
+        for (unsigned int t = 0; t < buf.size(); ++t)
+        {
+            r.corr.push_back(TensorRemove(buf[t]));
+        }
+        r.info.trace = 1;
+        result.push_back(r);
+        // two traces
+        corr = trace(ql*adj(gIn)*g5*adj(qbl)*g5*G)*trace(qbr*gOut*g5*adj(qr)*g5*G);
+        sliceSum(corr, buf, Tp);
+        r.corr.clear();
+        for (unsigned int t = 0; t < buf.size(); ++t)
+        {
+            r.corr.push_back(TensorRemove(buf[t]));
+        }
+        r.info.trace = 2;
+        result.push_back(r);
+    }
+    saveResult(par().output, "weakNonEye3pt", result);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MContraction_WeakNonEye3pt_hpp_

Hadrons/Modules/MFermion/FreeProp.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MFermion/FreeProp.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Vera Guelpers <vmg1n14@soton.ac.uk>

Hadrons/Modules/MFermion/FreeProp.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MFermion/FreeProp.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Vera Guelpers <vmg1n14@soton.ac.uk>

Hadrons/Modules/MFermion/GaugeProp.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MFermion/GaugeProp.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MFermion/GaugeProp.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MFermion/GaugeProp.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Guido Cossu <guido.cossu@ed.ac.uk>

Hadrons/Modules/MGauge/Electrify.cc

@@ -4,9 +4,11 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MGauge/Electrify.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
+Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Vera Guelpers <Vera.Guelpers@ed.ac.uk>
+Author: Vera Guelpers <vmg1n14@soton.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

Hadrons/Modules/MGauge/Electrify.hpp

@@ -4,9 +4,11 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MGauge/Electrify.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
+Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Vera Guelpers <Vera.Guelpers@ed.ac.uk>
+Author: Vera Guelpers <vmg1n14@soton.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

Hadrons/Modules/MGauge/FundtoHirep.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MGauge/FundtoHirep.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Guido Cossu <guido.cossu@ed.ac.uk>

Hadrons/Modules/MGauge/FundtoHirep.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MGauge/FundtoHirep.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: pretidav <david.preti@csic.es>

Hadrons/Modules/MGauge/GaugeFix.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MGauge/GaugeFix.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>

Hadrons/Modules/MGauge/GaugeFix.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MGauge/GaugeFix.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>

Hadrons/Modules/MGauge/Random.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MGauge/Random.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MGauge/Random.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MGauge/Random.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MGauge/StochEm.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MGauge/StochEm.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: James Harrison <jch1g10@soton.ac.uk>

Hadrons/Modules/MGauge/StochEm.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MGauge/StochEm.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: James Harrison <jch1g10@soton.ac.uk>

Hadrons/Modules/MGauge/StoutSmearing.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MGauge/StoutSmearing.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MGauge/StoutSmearing.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MGauge/StoutSmearing.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MGauge/Unit.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MGauge/Unit.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MGauge/Unit.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MGauge/Unit.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MGauge/UnitEm.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MGauge/UnitEm.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: James Harrison <jch1g10@soton.ac.uk>

Hadrons/Modules/MGauge/UnitEm.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MGauge/UnitEm.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: James Harrison <jch1g10@soton.ac.uk>

Hadrons/Modules/MIO/LoadA2AVectors.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MIO/LoadA2AVectors.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MIO/LoadA2AVectors.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MIO/LoadA2AVectors.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MIO/LoadBinary.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MIO/LoadBinary.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MIO/LoadBinary.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MIO/LoadBinary.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MIO/LoadCoarseEigenPack.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MIO/LoadCoarseEigenPack.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 
@@ -32,4 +32,6 @@ using namespace Hadrons;
 using namespace MIO;
 
 template class Grid::Hadrons::MIO::TLoadCoarseEigenPack<CoarseFermionEigenPack<FIMPL,HADRONS_DEFAULT_LANCZOS_NBASIS>>;
-
+#ifdef GRID_DEFAULT_PRECISION_DOUBLE
+template class Grid::Hadrons::MIO::TLoadCoarseEigenPack<CoarseFermionEigenPack<FIMPL,HADRONS_DEFAULT_LANCZOS_NBASIS, FIMPLF>>;
+#endif

Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 
@@ -56,7 +56,11 @@ template <typename Pack>
 class TLoadCoarseEigenPack: public Module<LoadCoarseEigenPackPar>
 {
 public:
-    typedef CoarseEigenPack<typename Pack::Field, typename Pack::CoarseField> BasePack;
+    typedef typename Pack::Field                Field;
+    typedef typename Pack::FieldIo              FieldIo;
+    typedef typename Pack::CoarseField          CoarseField;
+    typedef typename Pack::CoarseFieldIo        CoarseFieldIo;
+    typedef CoarseEigenPack<Field, CoarseField, FieldIo, CoarseFieldIo> BasePack;
     template <typename vtype> 
     using iImplScalar = iScalar<iScalar<iScalar<vtype>>>;
     typedef iImplScalar<typename Pack::Field::vector_type> SiteComplex;
@@ -74,7 +78,12 @@ public:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_TMP(LoadCoarseFermionEigenPack, ARG(TLoadCoarseEigenPack<CoarseFermionEigenPack<FIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>>), MIO);
+MODULE_REGISTER_TMP(LoadCoarseFermionEigenPack, 
+                    ARG(TLoadCoarseEigenPack<CoarseFermionEigenPack<FIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>>), MIO);
+#ifdef GRID_DEFAULT_PRECISION_DOUBLE
+MODULE_REGISTER_TMP(LoadCoarseFermionEigenPackIo32, 
+                    ARG(TLoadCoarseEigenPack<CoarseFermionEigenPack<FIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS, FIMPLF>>), MIO);
+#endif
 
 /******************************************************************************
  *                 TLoadCoarseEigenPack implementation                             *
@@ -106,18 +115,27 @@ std::vector<std::string> TLoadCoarseEigenPack<Pack>::getOutput(void)
 template <typename Pack>
 void TLoadCoarseEigenPack<Pack>::setup(void)
 {
-    env().createGrid(par().Ls);
-    env().createCoarseGrid(par().blockSize, par().Ls);
+    GridBase     *gridIo = nullptr, *gridCoarseIo = nullptr;
+
+    if (typeHash<Field>() != typeHash<FieldIo>())
+    {
+        gridIo = envGetRbGrid(FieldIo, par().Ls);
+    }
+    if (typeHash<CoarseField>() != typeHash<CoarseFieldIo>())
+    {
+        gridCoarseIo = envGetCoarseGrid(CoarseFieldIo, par().blockSize, par().Ls);
+    }
     envCreateDerived(BasePack, Pack, getName(), par().Ls, par().sizeFine,
-                     par().sizeCoarse, env().getRbGrid(par().Ls), 
-                     env().getCoarseGrid(par().blockSize, par().Ls));
+                     par().sizeCoarse, envGetRbGrid(Field, par().Ls), 
+                     envGetCoarseGrid(CoarseField, par().blockSize, par().Ls),
+                     gridIo, gridCoarseIo);
 }
 
 // execution ///////////////////////////////////////////////////////////////////
 template <typename Pack>
 void TLoadCoarseEigenPack<Pack>::execute(void)
 {
-    auto                 cg     = env().getCoarseGrid(par().blockSize, par().Ls);
+    auto                 cg     = envGetCoarseGrid(CoarseField, par().blockSize, par().Ls);
     auto                 &epack = envGetDerived(BasePack, Pack, getName());
     Lattice<SiteComplex> dummy(cg);
 

Hadrons/Modules/MIO/LoadCosmHol.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MIO/LoadCosmHol.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MIO/LoadCosmHol.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MIO/LoadCosmHol.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MIO/LoadEigenPack.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MIO/LoadEigenPack.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MIO/LoadEigenPack.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MIO/LoadEigenPack.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MIO/LoadNersc.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MIO/LoadNersc.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MIO/LoadNersc.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MIO/LoadNersc.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MNPR/Amputate.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MNPR/Amputate.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>

Hadrons/Modules/MNPR/Amputate.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MNPR/Amputate.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Julia Kettle J.R.Kettle-2@sms.ed.ac.uk

Hadrons/Modules/MNPR/Bilinear.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MNPR/Bilinear.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>

Hadrons/Modules/MNPR/Bilinear.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MNPR/Bilinear.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Julia Kettle J.R.Kettle-2@sms.ed.ac.uk

Hadrons/Modules/MNPR/FourQuark.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MNPR/FourQuark.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>

Hadrons/Modules/MNPR/FourQuark.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MNPR/FourQuark.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Julia Kettle J.R.Kettle-2@sms.ed.ac.uk

Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.cc

@@ -4,10 +4,11 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Vera Guelpers <Vera.Guelpers@ed.ac.uk>
+Author: Vera Guelpers <vmg1n14@soton.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

Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.hpp

@@ -4,10 +4,11 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Vera Guelpers <Vera.Guelpers@ed.ac.uk>
+Author: Vera Guelpers <vmg1n14@soton.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

Hadrons/Modules/MNoise/TimeDilutedSpinColorDiagonal.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MNoise/TimeDilutedSpinColorDiagonal.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MNoise/TimeDilutedSpinColorDiagonal.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MNoise/TimeDilutedSpinColorDiagonal.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MScalar/ChargedProp.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MScalar/ChargedProp.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: James Harrison <jch1g10@soton.ac.uk>

Hadrons/Modules/MScalar/ChargedProp.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MScalar/ChargedProp.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: James Harrison <jch1g10@soton.ac.uk>

Hadrons/Modules/MScalar/FreeProp.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MScalar/FreeProp.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MScalar/FreeProp.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MScalar/FreeProp.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MScalar/Scalar.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MScalar/Scalar.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MScalar/ScalarVP.cc

@@ -1,564 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: 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
-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/MScalar/ChargedProp.hpp>
-#include <Hadrons/Modules/MScalar/ScalarVP.hpp>
-#include <Hadrons/Modules/MScalar/Scalar.hpp>
-
-using namespace Grid;
-using namespace Hadrons;
-using namespace MScalar;
-
-/*
- * Scalar QED vacuum polarisation up to O(alpha)
- *
- * Conserved vector 2-point function diagram notation:
- *        _______
- *       /       \
- * U_nu *         * U_mu
- *       \_______/
- *
- *                (   adj(S(a\hat{nu}|x)) U_mu(x) S(0|x+a\hat{mu}) U_nu(0)    )
- *          = 2 Re(                             -                             )
- *                ( adj(S(a\hat{nu}|x+a\hat{mu})) adj(U_mu(x)) S(0|x) U_nu(0) )
- *  
- *
- *            _______
- *           /       \
- * free = 1 *         * 1
- *           \_______/
- *
- *
- *
- *             _______
- *            /       \
- * S = iA_nu *         * iA_mu
- *            \_______/
- *
- *
- *         Delta_1
- *         ___*___
- *        /       \
- * X = 1 *         * 1
- *        \___*___/
- *         Delta_1
- *
- *          Delta_1                     Delta_1
- *          ___*___                     ___*___
- *         /       \                   /       \
- *      1 *         * iA_mu  +  iA_nu *         * 1
- *         \_______/                   \_______/
- * 4C =        _______                     _______
- *            /       \                   /       \
- *      +  1 *         * iA_mu  +  iA_nu *         * 1
- *            \___*___/                   \___*___/
- *             Delta_1                     Delta_1
- *
- *     Delta_1   Delta_1
- *          _*___*_             _______
- *         /       \           /       \
- * 2E = 1 *         * 1  +  1 *         * 1
- *         \_______/           \_*___*_/
- *                         Delta_1   Delta_1
- *
- *          Delta_2
- *          ___*___             _______
- *         /       \           /       \
- * 2T = 1 *         * 1  +  1 *         * 1
- *         \_______/           \___*___/
- *                              Delta_2
- *
- *
- *                    _______
- *                   /       \
- * srcT = -A_nu^2/2 *         * 1
- *                   \_______/
- *
- *
- *
- *            _______
- *           /       \
- * snkT = 1 *         * -A_mu^2/2
- *           \_______/
- *
- * Full VP to O(alpha) = free + q^2*(S+X+4C+2E+2T+srcT+snkT)
- */
-
-/******************************************************************************
-*                  TScalarVP implementation                             *
-******************************************************************************/
-// constructor /////////////////////////////////////////////////////////////////
-TScalarVP::TScalarVP(const std::string name)
-: Module<ScalarVPPar>(name)
-{}
-
-// dependencies/products ///////////////////////////////////////////////////////
-std::vector<std::string> TScalarVP::getInput(void)
-{
-    prop0Name_ = par().scalarProp + "_0";
-    propQName_ = par().scalarProp + "_Q";
-    propSunName_ = par().scalarProp + "_Sun";
-    propTadName_ = par().scalarProp + "_Tad";
-
-	std::vector<std::string> in = {par().emField, prop0Name_, propQName_,
-                                   propSunName_, propTadName_};
-    
-    return in;
-}
-
-std::vector<std::string> TScalarVP::getOutput(void)
-{
-    std::vector<std::string> out;
-    
-    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
-    {
-        // out.push_back(getName() + "_propQ_" + std::to_string(mu));
-
-        for (unsigned int nu = 0; nu < env().getNd(); ++nu)
-        {
-            out.push_back(getName() + "_" + std::to_string(mu)
-                          + "_" + std::to_string(nu));
-        }
-    }
-
-    return out;
-}
-
-// setup ///////////////////////////////////////////////////////////////////////
-void TScalarVP::setup(void)
-{
-	freeMomPropName_ = FREEMOMPROP(static_cast<TChargedProp *>(vm().getModule(par().scalarProp))->par().mass);
-	GFSrcName_ = par().scalarProp + "_DinvSrc";
-    fftName_   = par().scalarProp + "_fft";
-	phaseName_.clear();
-	muPropQName_.clear();
-    vpTensorName_.clear();
-    momPhaseName_.clear();
-	for (unsigned int mu = 0; mu < env().getNd(); ++mu)
-    {
-        phaseName_.push_back("_shiftphase_" + std::to_string(mu));
-        muPropQName_.push_back(getName() + "_propQ_" + std::to_string(mu));
-
-        std::vector<std::string> vpTensorName_mu;
-        for (unsigned int nu = 0; nu < env().getNd(); ++nu)
-        {
-            vpTensorName_mu.push_back(getName() + "_" + std::to_string(mu)
-                                      + "_" + std::to_string(nu));
-        }
-        vpTensorName_.push_back(vpTensorName_mu);
-    }
-    if (!par().output.empty())
-    {
-        for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
-        {
-            momPhaseName_.push_back("_momentumphase_" + std::to_string(i_p));
-        }
-    }
-
-    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
-	{
-	    envCreateLat(ScalarField, muPropQName_[mu]);
-
-        for (unsigned int nu = 0; nu < env().getNd(); ++nu)
-        {
-            envCreateLat(ScalarField, vpTensorName_[mu][nu]);
-        }
-	}
-    if (!par().output.empty())
-    {
-        momPhasesDone_ = env().hasCreatedObject(momPhaseName_[0]);
-        for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
-        {
-            envCacheLat(ScalarField, momPhaseName_[i_p]);
-        }
-    }
-    envTmpLat(ScalarField, "buf");
-    envTmpLat(ScalarField, "result");
-    envTmpLat(ScalarField, "Amu");
-    envTmpLat(ScalarField, "Usnk");
-    envTmpLat(ScalarField, "tmpProp");
-}
-
-// execution ///////////////////////////////////////////////////////////////////
-void TScalarVP::execute(void)
-{
-    // CACHING ANALYTIC EXPRESSIONS
-    makeCaches();
-
-    Complex ci(0.0,1.0);
-    Real    q        = static_cast<TChargedProp *>(vm().getModule(par().scalarProp))->par().charge;
-    auto    &prop0   = envGet(ScalarField, prop0Name_);
-    auto    &propQ   = envGet(ScalarField, propQName_);
-    auto    &propSun = envGet(ScalarField, propSunName_);
-    auto    &propTad = envGet(ScalarField, propTadName_);
-    auto    &GFSrc   = envGet(ScalarField, GFSrcName_);
-    auto    &G       = envGet(ScalarField, freeMomPropName_);
-    auto    &fft     = envGet(FFT, fftName_);
-    phase_.clear();
-    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
-    {
-        auto &phmu = envGet(ScalarField, phaseName_[mu]);
-        phase_.push_back(&phmu);
-    }
-    
-    // PROPAGATORS FROM SHIFTED SOURCES
-    LOG(Message) << "Computing O(q) charged scalar propagators..."
-                 << std::endl;
-    std::vector<ScalarField *> muPropQ;
-    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
-    {
-        auto &propmu = envGet(ScalarField, muPropQName_[mu]);
-
-        // -G*momD1*G*F*tau_mu*Src (momD1 = F*D1*Finv)
-        propmu = adj(*phase_[mu])*GFSrc;
-        momD1(propmu, fft);
-        propmu = -G*propmu;
-        fft.FFT_all_dim(propmu, propmu, FFT::backward);
-
-        muPropQ.push_back(&propmu);
-    }
-
-    // CONTRACTIONS
-    auto        &A = envGet(EmField, par().emField);
-    envGetTmp(ScalarField, buf);
-    envGetTmp(ScalarField, result);
-    envGetTmp(ScalarField, Amu);
-    envGetTmp(ScalarField, Usnk);
-    envGetTmp(ScalarField, tmpProp);
-    TComplex    Anu0, Usrc;
-    std::vector<int> coor0 = {0, 0, 0, 0};
-    std::vector<std::vector<ScalarField *> > vpTensor;
-    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
-    {
-        std::vector<ScalarField *> vpTensor_mu;
-        for (unsigned int nu = 0; nu < env().getNd(); ++nu)
-        {
-            auto &vpmunu = envGet(ScalarField, vpTensorName_[mu][nu]);
-            vpTensor_mu.push_back(&vpmunu);
-        }
-        vpTensor.push_back(vpTensor_mu);
-    }
-
-    // Prepare output data structure if necessary
-    Result outputData;
-    if (!par().output.empty())
-    {
-        outputData.projection.resize(par().outputMom.size());
-        outputData.lattice_size = env().getGrid()->_fdimensions;
-        outputData.mass = static_cast<TChargedProp *>(vm().getModule(par().scalarProp))->par().mass;
-        outputData.charge = q;
-        for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
-        {
-            outputData.projection[i_p].momentum = strToVec<int>(par().outputMom[i_p]);
-            outputData.projection[i_p].pi.resize(env().getNd());
-            outputData.projection[i_p].pi_free.resize(env().getNd());
-            outputData.projection[i_p].pi_2E.resize(env().getNd());
-            outputData.projection[i_p].pi_2T.resize(env().getNd());
-            outputData.projection[i_p].pi_S.resize(env().getNd());
-            outputData.projection[i_p].pi_4C.resize(env().getNd());
-            outputData.projection[i_p].pi_X.resize(env().getNd());
-            outputData.projection[i_p].pi_srcT.resize(env().getNd());
-            outputData.projection[i_p].pi_snkT.resize(env().getNd());
-            for (unsigned int nu = 0; nu < env().getNd(); ++nu)
-            {
-                outputData.projection[i_p].pi[nu].resize(env().getNd());
-                outputData.projection[i_p].pi_free[nu].resize(env().getNd());
-                outputData.projection[i_p].pi_2E[nu].resize(env().getNd());
-                outputData.projection[i_p].pi_2T[nu].resize(env().getNd());
-                outputData.projection[i_p].pi_S[nu].resize(env().getNd());
-                outputData.projection[i_p].pi_4C[nu].resize(env().getNd());
-                outputData.projection[i_p].pi_X[nu].resize(env().getNd());
-                outputData.projection[i_p].pi_srcT[nu].resize(env().getNd());
-                outputData.projection[i_p].pi_snkT[nu].resize(env().getNd());
-            }
-        }
-    }
-
-    // Do contractions
-    for (unsigned int nu = 0; nu < env().getNd(); ++nu)
-    {
-        peekSite(Anu0, peekLorentz(A, nu), coor0);
-
-        for (unsigned int mu = 0; mu < env().getNd(); ++mu)
-        {
-            LOG(Message) << "Computing Pi[" << mu << "][" << nu << "]..."
-                         << std::endl;
-            Amu = peekLorentz(A, mu);
-
-            // free
-            tmpProp = Cshift(prop0, nu, -1);     // S_0(0|x-a\hat{\nu})
-                                                 // = S_0(a\hat{\nu}|x)
-            Usrc    = Complex(1.0,0.0);
-            vpContraction(result, prop0, tmpProp, Usrc, mu);
-            *vpTensor[mu][nu] = result;
-            // Do momentum projections if necessary
-            if (!par().output.empty())
-            {
-                for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
-                {
-                    project(outputData.projection[i_p].pi_free[mu][nu], result,
-                            i_p);
-                }
-            }
-            tmpProp = result; // Just using tmpProp as a temporary ScalarField
-                              // here (buf is modified by calls to writeVP())
-
-            // srcT
-            result = tmpProp * (-0.5)*Anu0*Anu0;
-            *vpTensor[mu][nu] += q*q*result;
-            // Do momentum projections if necessary
-            if (!par().output.empty())
-            {
-                for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
-                {
-                    project(outputData.projection[i_p].pi_srcT[mu][nu], result,
-                            i_p);
-                }
-            }
-
-            // snkT
-            result = tmpProp * (-0.5)*Amu*Amu;
-            *vpTensor[mu][nu] += q*q*result;
-            // Do momentum projections if necessary
-            if (!par().output.empty())
-            {
-                for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
-                {
-                    project(outputData.projection[i_p].pi_snkT[mu][nu], result,
-                            i_p);
-                }
-            }
-
-            // S
-            tmpProp = Cshift(prop0, nu, -1);     // S_0(a\hat{\nu}|x)
-            Usrc    = ci*Anu0;
-            Usnk    = ci*Amu;
-            vpContraction(result, prop0, tmpProp, Usrc, Usnk, mu);
-            *vpTensor[mu][nu] += q*q*result;
-            // Do momentum projections if necessary
-            if (!par().output.empty())
-            {
-                for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
-                {
-                    project(outputData.projection[i_p].pi_S[mu][nu], result,
-                            i_p);
-                }
-            }
-
-            // 4C
-            tmpProp = Cshift(prop0, nu, -1);     // S_0(a\hat{\nu}|x)
-            Usrc    = Complex(1.0,0.0);
-            Usnk    = ci*Amu;
-            vpContraction(result, propQ, tmpProp, Usrc, Usnk, mu);
-            Usrc    = ci*Anu0;
-            vpContraction(buf, propQ, tmpProp, Usrc, mu);
-            result += buf;
-            vpContraction(buf, prop0, *muPropQ[nu], Usrc, mu);
-            result += buf;
-            Usrc = Complex(1.0,0.0);
-            Usnk = ci*Amu;
-            vpContraction(buf, prop0, *muPropQ[nu], Usrc, Usnk, mu);
-            result += buf;
-            *vpTensor[mu][nu] += q*q*result;
-            // Do momentum projections if necessary
-            if (!par().output.empty())
-            {
-                for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
-                {
-                    project(outputData.projection[i_p].pi_4C[mu][nu], result,
-                            i_p);
-                }
-            }
-
-            // X
-            Usrc = Complex(1.0,0.0);
-            vpContraction(result, propQ, *muPropQ[nu], Usrc, mu);
-            *vpTensor[mu][nu] += q*q*result;
-            // Do momentum projections if necessary
-            if (!par().output.empty())
-            {
-                for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
-                {
-                    project(outputData.projection[i_p].pi_X[mu][nu], result,
-                            i_p);
-                }
-            }
-
-            // 2E
-            tmpProp = Cshift(prop0, nu, -1);     // S_0(a\hat{\nu}|x)
-            Usrc    = Complex(1.0,0.0);
-            vpContraction(result, propSun, tmpProp, Usrc, mu);
-            tmpProp = Cshift(propSun, nu, -1);     // S_\Sigma(0|x-a\hat{\nu})
-                               //(Note: <S(0|x-a\hat{\nu})> = <S(a\hat{\nu}|x)>)
-            vpContraction(buf, prop0, tmpProp, Usrc, mu);
-            result += buf;
-            *vpTensor[mu][nu] += q*q*result;
-            // Do momentum projections if necessary
-            if (!par().output.empty())
-            {
-                for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
-                {
-                    project(outputData.projection[i_p].pi_2E[mu][nu], result,
-                            i_p);
-                }
-            }
-
-            // 2T
-            tmpProp = Cshift(prop0, nu, -1);     // S_0(a\hat{\nu}|x)
-            Usrc    = Complex(1.0,0.0);
-            vpContraction(result, propTad, tmpProp, Usrc, mu);
-            tmpProp = Cshift(propTad, nu, -1);     // S_T(0|x-a\hat{\nu})
-            vpContraction(buf, prop0, tmpProp, Usrc, mu);
-            result += buf;
-            *vpTensor[mu][nu] += q*q*result;
-            // Do momentum projections if necessary
-            if (!par().output.empty())
-            {
-                for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
-                {
-                    project(outputData.projection[i_p].pi_2T[mu][nu], result,
-                            i_p);
-                }
-            }
-
-            // Do momentum projections of full VP if necessary
-            if (!par().output.empty())
-            {
-                for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
-                {
-                    project(outputData.projection[i_p].pi[mu][nu],
-                            *vpTensor[mu][nu], i_p);
-                }
-            }
-        }
-    }
-
-    // OUTPUT IF NECESSARY
-    if (!par().output.empty())
-    {
-        LOG(Message) << "Saving momentum-projected HVP to '"
-                     << RESULT_FILE_NAME(par().output, vm().getTrajectory()) << "'..."
-                     << std::endl;
-        saveResult(par().output, "HVP", outputData);
-    }
-}
-
-void TScalarVP::makeCaches(void)
-{
-    envGetTmp(ScalarField, buf);
-
-    if ( (!par().output.empty()) && (!momPhasesDone_) )
-    {
-        LOG(Message) << "Caching phases for momentum projections..."
-                     << std::endl;
-        std::vector<int> &l = env().getGrid()->_fdimensions;
-        Complex          ci(0.0,1.0);
-
-        // Calculate phase factors
-        for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
-        {
-            std::vector<int> mom = strToVec<int>(par().outputMom[i_p]);
-            auto &momph_ip = envGet(ScalarField, momPhaseName_[i_p]);
-            momph_ip = zero;
-            for (unsigned int j = 0; j < env().getNd()-1; ++j)
-            {
-                Real twoPiL = M_PI*2./l[j];
-                LatticeCoordinate(buf, j);
-                buf = mom[j]*twoPiL*buf;
-                momph_ip = momph_ip + buf;
-            }
-            momph_ip = exp(-ci*momph_ip);
-            momPhase_.push_back(&momph_ip);
-        }
-    }
-}
-
-void TScalarVP::vpContraction(ScalarField &vp,
-                   ScalarField &prop_0_x, ScalarField &prop_nu_x,
-                   TComplex u_src, ScalarField &u_snk, int mu)
-{
-    // Note: this function assumes a point source is used.
-    vp = adj(prop_nu_x) * u_snk * Cshift(prop_0_x, mu, 1) * u_src;
-    vp -= Cshift(adj(prop_nu_x), mu, 1) * adj(u_snk) * prop_0_x * u_src;
-    vp = 2.0*real(vp);
-}
-
-void TScalarVP::vpContraction(ScalarField &vp,
-                   ScalarField &prop_0_x, ScalarField &prop_nu_x,
-                   TComplex u_src, int mu)
-{
-    // Note: this function assumes a point source is used.
-    vp = adj(prop_nu_x) * Cshift(prop_0_x, mu, 1) * u_src;
-    vp -= Cshift(adj(prop_nu_x), mu, 1) * prop_0_x * u_src;
-    vp = 2.0*real(vp);
-}
-
-void TScalarVP::project(std::vector<Complex> &projection, const ScalarField &vp, int i_p)
-{
-    std::vector<TComplex>   vecBuf;
-    envGetTmp(ScalarField, buf);
-
-    buf = vp*(*momPhase_[i_p]);
-    sliceSum(buf, vecBuf, Tp);
-    projection.resize(vecBuf.size());
-    for (unsigned int t = 0; t < vecBuf.size(); ++t)
-    {
-        projection[t] = TensorRemove(vecBuf[t]);
-    }
-}
-
-void TScalarVP::momD1(ScalarField &s, FFT &fft)
-{
-    auto        &A = envGet(EmField, par().emField);
-    Complex     ci(0.0,1.0);
-
-    envGetTmp(ScalarField, buf);
-    envGetTmp(ScalarField, result);
-    envGetTmp(ScalarField, Amu);
-
-    result = zero;
-    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
-    {
-        Amu = peekLorentz(A, mu);
-        buf = (*phase_[mu])*s;
-        fft.FFT_all_dim(buf, buf, FFT::backward);
-        buf = Amu*buf;
-        fft.FFT_all_dim(buf, buf, FFT::forward);
-        result = result - ci*buf;
-    }
-    fft.FFT_all_dim(s, s, FFT::backward);
-    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
-    {
-        Amu = peekLorentz(A, mu);
-        buf = Amu*s;
-        fft.FFT_all_dim(buf, buf, FFT::forward);
-        result = result + ci*adj(*phase_[mu])*buf;
-    }
-
-    s = result;
-}

Hadrons/Modules/MScalar/ScalarVP.hpp

@@ -1,129 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: 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
-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_MScalar_ScalarVP_hpp_
-#define Hadrons_MScalar_ScalarVP_hpp_
-
-#include <Hadrons/Global.hpp>
-#include <Hadrons/Module.hpp>
-#include <Hadrons/ModuleFactory.hpp>
-
-BEGIN_HADRONS_NAMESPACE
-
-/******************************************************************************
- *                         Scalar vacuum polarisation                         *
- ******************************************************************************/
-BEGIN_MODULE_NAMESPACE(MScalar)
-
-class ScalarVPPar: Serializable
-{
-public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(ScalarVPPar,
-                                    std::string, emField,
-                                    std::string, scalarProp,
-                                    std::string, output,
-                                    std::vector<std::string>, outputMom);
-};
-
-class TScalarVP: public Module<ScalarVPPar>
-{
-public:
-    BASIC_TYPE_ALIASES(SIMPL,);
-    typedef PhotonR::GaugeField     EmField;
-    typedef PhotonR::GaugeLinkField EmComp;
-    class Result: Serializable
-    {
-    public:
-        class Projection: Serializable
-        {
-        public:
-            GRID_SERIALIZABLE_CLASS_MEMBERS(Projection,
-                                            std::vector<int>,     momentum,
-                                            std::vector<std::vector<std::vector<Complex>>>, pi,
-                                            std::vector<std::vector<std::vector<Complex>>>, pi_free,
-                                            std::vector<std::vector<std::vector<Complex>>>, pi_2E,
-                                            std::vector<std::vector<std::vector<Complex>>>, pi_2T,
-                                            std::vector<std::vector<std::vector<Complex>>>, pi_S,
-                                            std::vector<std::vector<std::vector<Complex>>>, pi_4C,
-                                            std::vector<std::vector<std::vector<Complex>>>, pi_X,
-                                            std::vector<std::vector<std::vector<Complex>>>, pi_srcT,
-                                            std::vector<std::vector<std::vector<Complex>>>, pi_snkT);
-        };
-        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
-                                        std::vector<int>,        lattice_size,
-                                        double,                  mass,
-                                        double,                  charge,
-                                        std::vector<Projection>, projection);
-    };
-public:
-    // constructor
-    TScalarVP(const std::string name);
-    // destructor
-    virtual ~TScalarVP(void) {};
-    // dependency relation
-    virtual std::vector<std::string> getInput(void);
-    virtual std::vector<std::string> getOutput(void);
-protected:
-    // setup
-    virtual void setup(void);
-    // execution
-    virtual void execute(void);
-private:
-    void makeCaches(void);
-    // conserved vector two-point contraction
-    void vpContraction(ScalarField &vp,
-                       ScalarField &prop_0_x, ScalarField &prop_nu_x,
-                       TComplex u_src, ScalarField &u_snk, int mu);
-    // conserved vector two-point contraction with unit gauge link at sink
-    void vpContraction(ScalarField &vp,
-                       ScalarField &prop_0_x, ScalarField &prop_nu_x,
-                       TComplex u_src, int mu);
-    // write momentum-projected vacuum polarisation to file(s)
-    void project(std::vector<Complex> &projection, const ScalarField &vp,
-                 int i_p);
-    // momentum-space Delta_1 insertion
-    void momD1(ScalarField &s, FFT &fft);
-private:
-    bool                                        momPhasesDone_;
-    std::string                                 freeMomPropName_, GFSrcName_,
-                                                prop0Name_, propQName_,
-                                                propSunName_, propTadName_,
-                                                fftName_;
-    std::vector<std::string>                    phaseName_, muPropQName_,
-                                                momPhaseName_;
-    std::vector<std::vector<std::string> >      vpTensorName_;
-    std::vector<ScalarField *>                  phase_, momPhase_;
-};
-
-MODULE_REGISTER(ScalarVP, TScalarVP, MScalar);
-
-END_MODULE_NAMESPACE
-
-END_HADRONS_NAMESPACE
-
-#endif // Hadrons_MScalar_ScalarVP_hpp_

Hadrons/Modules/MScalar/VPCounterTerms.cc

@@ -1,260 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: 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
-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/MScalar/VPCounterTerms.hpp>
-#include <Hadrons/Modules/MScalar/Scalar.hpp>
-
-using namespace Grid;
-using namespace Hadrons;
-using namespace MScalar;
-
-/******************************************************************************
-*                  TVPCounterTerms implementation                             *
-******************************************************************************/
-// constructor /////////////////////////////////////////////////////////////////
-TVPCounterTerms::TVPCounterTerms(const std::string name)
-: Module<VPCounterTermsPar>(name)
-{}
-
-// dependencies/products ///////////////////////////////////////////////////////
-std::vector<std::string> TVPCounterTerms::getInput(void)
-{
-    std::vector<std::string> in = {par().source};
-    
-    return in;
-}
-
-std::vector<std::string> TVPCounterTerms::getOutput(void)
-{
-    std::vector<std::string> out;
-    
-    return out;
-}
-
-// setup ///////////////////////////////////////////////////////////////////////
-void TVPCounterTerms::setup(void)
-{
-	freeMomPropName_ = FREEMOMPROP(par().mass);
-    phaseName_.clear();
-    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
-    {
-        phaseName_.push_back("_shiftphase_" + std::to_string(mu));
-    }
-    GFSrcName_ = getName() + "_DinvSrc";
-    phatsqName_ = getName() + "_pHatSquared";
-    prop0Name_ = getName() + "_freeProp";
-    twoscalarName_ = getName() + "_2scalarProp";
-    psquaredName_ = getName() + "_psquaredProp";
-    if (!par().output.empty())
-    {
-        for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
-        {
-            momPhaseName_.push_back("_momentumphase_" + std::to_string(i_p));
-        }
-    }
-
-    envCreateLat(ScalarField, freeMomPropName_);
-    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
-    {
-        envCreateLat(ScalarField, phaseName_[mu]);
-    }
-    envCreateLat(ScalarField, phatsqName_);
-    envCreateLat(ScalarField, GFSrcName_);
-    envCreateLat(ScalarField, prop0Name_);
-    envCreateLat(ScalarField, twoscalarName_);
-    envCreateLat(ScalarField, psquaredName_);
-    if (!par().output.empty())
-    {
-        for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
-        {
-            envCacheLat(ScalarField, momPhaseName_[i_p]);
-        }
-    }
-    envTmpLat(ScalarField, "buf");
-    envTmpLat(ScalarField, "tmp_vp");
-    envTmpLat(ScalarField, "vpPhase");
-}
-
-// execution ///////////////////////////////////////////////////////////////////
-void TVPCounterTerms::execute(void)
-{
-	auto &source = envGet(ScalarField, par().source);
-    Complex     ci(0.0,1.0);
-    FFT         fft(env().getGrid());
-    envGetTmp(ScalarField, buf);
-    envGetTmp(ScalarField, tmp_vp);
-    
-    // Momentum-space free scalar propagator
-    auto &G = envGet(ScalarField, freeMomPropName_);
-    SIMPL::MomentumSpacePropagator(G, par().mass);
-
-    // Phases and hat{p}^2
-    auto &phatsq = envGet(ScalarField, phatsqName_);
-    std::vector<int> &l = env().getGrid()->_fdimensions;
-    
-    LOG(Message) << "Calculating shift phases..." << std::endl;
-    phatsq = zero;
-    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
-    {
-        Real    twoPiL = M_PI*2./l[mu];
-        auto &phmu  = envGet(ScalarField, phaseName_[mu]);
-
-        LatticeCoordinate(buf, mu);
-        phmu = exp(ci*twoPiL*buf);
-        phase_.push_back(&phmu);
-        buf = 2.*sin(.5*twoPiL*buf);
-		phatsq = phatsq + buf*buf;
-    }
-
-    // G*F*src
-    auto &GFSrc       = envGet(ScalarField, GFSrcName_);
-    fft.FFT_all_dim(GFSrc, source, FFT::forward);
-    GFSrc = G*GFSrc;
-
-    // Position-space free scalar propagator
-    auto &prop0       = envGet(ScalarField, prop0Name_);
-    prop0 = GFSrc;
-    fft.FFT_all_dim(prop0, prop0, FFT::backward);
-
-    // Propagators for counter-terms
-    auto &twoscalarProp        = envGet(ScalarField, twoscalarName_);
-    auto &psquaredProp         = envGet(ScalarField, psquaredName_);
-
-    twoscalarProp = G*GFSrc;
-    fft.FFT_all_dim(twoscalarProp, twoscalarProp, FFT::backward);
-
-    psquaredProp = G*phatsq*GFSrc;
-    fft.FFT_all_dim(psquaredProp, psquaredProp, FFT::backward);
-
-    // Prepare output data structure if necessary
-    Result outputData;
-    if (!par().output.empty())
-    {
-        outputData.projection.resize(par().outputMom.size());
-        outputData.lattice_size = env().getGrid()->_fdimensions;
-        outputData.mass = par().mass;
-        for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
-        {
-            outputData.projection[i_p].momentum = strToVec<int>(par().outputMom[i_p]);
-            outputData.projection[i_p].twoScalar.resize(env().getNd());
-            outputData.projection[i_p].threeScalar.resize(env().getNd());
-            outputData.projection[i_p].pSquaredInsertion.resize(env().getNd());
-            for (unsigned int nu = 0; nu < env().getNd(); ++nu)
-            {
-                outputData.projection[i_p].twoScalar[nu].resize(env().getNd());
-                outputData.projection[i_p].threeScalar[nu].resize(env().getNd());
-                outputData.projection[i_p].pSquaredInsertion[nu].resize(env().getNd());
-            }
-            // Calculate phase factors
-            auto &momph_ip = envGet(ScalarField, momPhaseName_[i_p]);
-            momph_ip = zero;
-            for (unsigned int j = 0; j < env().getNd()-1; ++j)
-            {
-                Real twoPiL = M_PI*2./l[j];
-                LatticeCoordinate(buf, j);
-                buf = outputData.projection[i_p].momentum[j]*twoPiL*buf;
-                momph_ip = momph_ip + buf;
-            }
-            momph_ip = exp(-ci*momph_ip);
-            momPhase_.push_back(&momph_ip);
-        }
-    }
-
-    // Contractions
-    for (unsigned int nu = 0; nu < env().getNd(); ++nu)
-    {
-    	buf = adj(Cshift(prop0, nu, -1));
-        for (unsigned int mu = 0; mu < env().getNd(); ++mu)
-        {
-            // Two-scalar loop
-            tmp_vp = buf * Cshift(prop0, mu, 1);
-            tmp_vp -= Cshift(buf, mu, 1) * prop0;
-            tmp_vp = 2.0*real(tmp_vp);
-            // Output if necessary
-            if (!par().output.empty())
-            {
-                for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
-                {
-                    project(outputData.projection[i_p].twoScalar[mu][nu],
-                            tmp_vp, i_p);
-                }
-            }
-
-        	// Three-scalar loop (no vertex)
-    		tmp_vp = buf * Cshift(twoscalarProp, mu, 1);
-            tmp_vp -= Cshift(buf, mu, 1) * twoscalarProp;
-            tmp_vp = 2.0*real(tmp_vp);
-            // Output if necessary
-            if (!par().output.empty())
-            {
-                for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
-                {
-                    project(outputData.projection[i_p].threeScalar[mu][nu],
-                            tmp_vp, i_p);
-                }
-            }
-
-            // Three-scalar loop (hat{p}^2 insertion)
-    		tmp_vp = buf * Cshift(psquaredProp, mu, 1);
-            tmp_vp -= Cshift(buf, mu, 1) * psquaredProp;
-            tmp_vp = 2.0*real(tmp_vp);
-            // Output if necessary
-            if (!par().output.empty())
-            {
-                for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p)
-                {
-                    project(outputData.projection[i_p].pSquaredInsertion[mu][nu],
-                            tmp_vp, i_p);
-                }
-            }
-        }
-    }
-
-    // OUTPUT IF NECESSARY
-    if (!par().output.empty())
-    {
-        LOG(Message) << "Saving momentum-projected correlators to '"
-                     << RESULT_FILE_NAME(par().output, vm().getTrajectory()) << "'..."
-                     << std::endl;
-        saveResult(par().output, "scalar_loops", outputData);
-    }
-}
-
-void TVPCounterTerms::project(std::vector<Complex> &projection, const ScalarField &vp, int i_p)
-{
-    std::vector<TComplex>   vecBuf;
-    envGetTmp(ScalarField, vpPhase);
-
-    vpPhase = vp*(*momPhase_[i_p]);
-    sliceSum(vpPhase, vecBuf, Tp);
-    projection.resize(vecBuf.size());
-    for (unsigned int t = 0; t < vecBuf.size(); ++t)
-    {
-        projection[t] = TensorRemove(vecBuf[t]);
-    }
-}

Hadrons/Modules/MScalar/VPCounterTerms.hpp

@@ -1,103 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: 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
-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_MScalar_VPCounterTerms_hpp_
-#define Hadrons_MScalar_VPCounterTerms_hpp_
-
-#include <Hadrons/Global.hpp>
-#include <Hadrons/Module.hpp>
-#include <Hadrons/ModuleFactory.hpp>
-
-BEGIN_HADRONS_NAMESPACE
-
-/******************************************************************************
- *                         VPCounterTerms                                 *
- ******************************************************************************/
-BEGIN_MODULE_NAMESPACE(MScalar)
-
-class VPCounterTermsPar: Serializable
-{
-public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(VPCounterTermsPar,
-                                    std::string, source,
-                                    double,      mass,
-                                    std::string, output,
-                                    std::vector<std::string>, outputMom);
-};
-
-class TVPCounterTerms: public Module<VPCounterTermsPar>
-{
-public:
-    BASIC_TYPE_ALIASES(SIMPL,);
-    class Result: Serializable
-    {
-    public:
-        class Projection: Serializable
-        {
-        public:
-            GRID_SERIALIZABLE_CLASS_MEMBERS(Projection,
-                                            std::vector<int>,     momentum,
-                                            std::vector<std::vector<std::vector<Complex>>>, twoScalar,
-                                            std::vector<std::vector<std::vector<Complex>>>, threeScalar,
-                                            std::vector<std::vector<std::vector<Complex>>>, pSquaredInsertion);
-        };
-        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
-                                        std::vector<int>,        lattice_size,
-                                        double,                  mass,
-                                        std::vector<Projection>, projection);
-    };
-public:
-    // constructor
-    TVPCounterTerms(const std::string name);
-    // destructor
-    virtual ~TVPCounterTerms(void) {};
-    // dependency relation
-    virtual std::vector<std::string> getInput(void);
-    virtual std::vector<std::string> getOutput(void);
-protected:
-    // setup
-    virtual void setup(void);
-    // execution
-    virtual void execute(void);
-private:
-    void project(std::vector<Complex> &projection, const ScalarField &vp, int i_p);
-private:
-    std::string                freeMomPropName_, GFSrcName_, phatsqName_, prop0Name_,
-                               twoscalarName_, twoscalarVertexName_,
-                               psquaredName_, psquaredVertexName_;
-    std::vector<std::string>   phaseName_, momPhaseName_;
-    std::vector<ScalarField *> phase_, momPhase_;
-};
-
-MODULE_REGISTER(VPCounterTerms, TVPCounterTerms, MScalar);
-
-END_MODULE_NAMESPACE
-
-END_HADRONS_NAMESPACE
-
-#endif // Hadrons_MScalar_VPCounterTerms_hpp_

Hadrons/Modules/MScalarSUN/Div.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MScalarSUN/Div.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MScalarSUN/Div.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MScalarSUN/Div.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MScalarSUN/EMT.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MScalarSUN/EMT.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MScalarSUN/EMT.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MScalarSUN/EMT.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MScalarSUN/Grad.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MScalarSUN/Grad.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MScalarSUN/Grad.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MScalarSUN/Grad.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MScalarSUN/ShiftProbe.cc

@@ -1,39 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: Hadrons/Modules/MScalarSUN/ShiftProbe.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 <Hadrons/Modules/MScalarSUN/ShiftProbe.hpp>
-
-using namespace Grid;
-using namespace Hadrons;
-using namespace MScalarSUN;
-
-template class Grid::Hadrons::MScalarSUN::TShiftProbe<ScalarNxNAdjImplR<2>>;
-template class Grid::Hadrons::MScalarSUN::TShiftProbe<ScalarNxNAdjImplR<3>>;
-template class Grid::Hadrons::MScalarSUN::TShiftProbe<ScalarNxNAdjImplR<4>>;
-template class Grid::Hadrons::MScalarSUN::TShiftProbe<ScalarNxNAdjImplR<5>>;
-template class Grid::Hadrons::MScalarSUN::TShiftProbe<ScalarNxNAdjImplR<6>>;
-

Hadrons/Modules/MScalarSUN/ShiftProbe.hpp

@@ -1,177 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: Hadrons/Modules/MScalarSUN/ShiftProbe.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_MScalarSUN_ShiftProbe_hpp_
-#define Hadrons_MScalarSUN_ShiftProbe_hpp_
-
-#include <Hadrons/Global.hpp>
-#include <Hadrons/Module.hpp>
-#include <Hadrons/ModuleFactory.hpp>
-#include <Hadrons/Modules/MScalarSUN/Utils.hpp>
-
-BEGIN_HADRONS_NAMESPACE
-
-/******************************************************************************
- *         Ward identity phi^n probe with fields at different positions       *
- ******************************************************************************/
-BEGIN_MODULE_NAMESPACE(MScalarSUN)
-
-typedef std::pair<int, int> ShiftPair;
-
-class ShiftProbePar: Serializable
-{
-public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(ShiftProbePar,
-                                    std::string, field,
-                                    std::string, shifts,
-                                    std::string, output);
-};
-
-class ShiftProbeResult: Serializable
-{
-public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(ShiftProbeResult,
-                                    std::string, shifts,
-                                    Complex,     value);
-};
-
-template <typename SImpl>
-class TShiftProbe: public Module<ShiftProbePar>
-{
-public:
-    typedef typename SImpl::Field                          Field;
-    typedef typename SImpl::ComplexField                   ComplexField;
-public:
-    // constructor
-    TShiftProbe(const std::string name);
-    // destructor
-    virtual ~TShiftProbe(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(ShiftProbeSU2, TShiftProbe<ScalarNxNAdjImplR<2>>, MScalarSUN);
-MODULE_REGISTER_TMP(ShiftProbeSU3, TShiftProbe<ScalarNxNAdjImplR<3>>, MScalarSUN);
-MODULE_REGISTER_TMP(ShiftProbeSU4, TShiftProbe<ScalarNxNAdjImplR<4>>, MScalarSUN);
-MODULE_REGISTER_TMP(ShiftProbeSU5, TShiftProbe<ScalarNxNAdjImplR<5>>, MScalarSUN);
-MODULE_REGISTER_TMP(ShiftProbeSU6, TShiftProbe<ScalarNxNAdjImplR<6>>, MScalarSUN);
-
-/******************************************************************************
- *                        TShiftProbe implementation                          *
- ******************************************************************************/
-// constructor /////////////////////////////////////////////////////////////////
-template <typename SImpl>
-TShiftProbe<SImpl>::TShiftProbe(const std::string name)
-: Module<ShiftProbePar>(name)
-{}
-
-// dependencies/products ///////////////////////////////////////////////////////
-template <typename SImpl>
-std::vector<std::string> TShiftProbe<SImpl>::getInput(void)
-{
-    std::vector<std::string> in = {par().field};
-    
-    return in;
-}
-
-template <typename SImpl>
-std::vector<std::string> TShiftProbe<SImpl>::getOutput(void)
-{
-    std::vector<std::string> out = {getName()};
-    
-    return out;
-}
-
-// setup ///////////////////////////////////////////////////////////////////////
-template <typename SImpl>
-void TShiftProbe<SImpl>::setup(void)
-{
-    envTmpLat(Field, "acc");
-    envCreateLat(ComplexField, getName());
-}
-
-// execution ///////////////////////////////////////////////////////////////////
-template <typename SImpl>
-void TShiftProbe<SImpl>::execute(void)
-{
-    LOG(Message) << "Creating shift probe for shifts " << par().shifts
-                 << std::endl;
-
-    std::vector<ShiftPair> shift;
-    double                 sign;
-    auto                   &phi   = envGet(Field, par().field);
-    auto                   &probe = envGet(ComplexField, getName());
-
-    shift = strToVec<ShiftPair>(par().shifts);
-    if (shift.size() % 2 != 0)
-    {
-        HADRONS_ERROR(Size, "the number of shifts is odd");
-    }
-    sign = (shift.size() % 4 == 0) ? 1. : -1.;
-    for (auto &s: shift)
-    {
-        if (s.first >= env().getNd())
-        {
-            HADRONS_ERROR(Size, "dimension to large for shift <" 
-                               + std::to_string(s.first) + " " 
-                               + std::to_string(s.second) + ">" );
-        }
-    }
-    envGetTmp(Field, acc);
-    acc = 1.;
-    for (unsigned int i = 0; i < shift.size(); ++i)
-    {
-        if (shift[i].second == 0)
-        {
-            acc *= phi;
-        }
-        else
-        {
-            acc *= Cshift(phi, shift[i].first, shift[i].second);
-        }
-    }
-    probe = sign*trace(acc);
-    if (!par().output.empty())
-    {
-        ShiftProbeResult r;
-
-        r.shifts = par().shifts;
-        r.value  = TensorRemove(sum(probe));
-        saveResult(par().output, "probe", r);
-    }
-}
-
-END_MODULE_NAMESPACE
-
-END_HADRONS_NAMESPACE
-
-#endif // Hadrons_MScalarSUN_ShiftProbe_hpp_

Hadrons/Modules/MScalarSUN/StochFreeField.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MScalarSUN/StochFreeField.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MScalarSUN/StochFreeField.hpp

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MScalarSUN/StochFreeField.hpp
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 

Hadrons/Modules/MScalarSUN/TrKinetic.cc

@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: Hadrons/Modules/MScalarSUN/TrKinetic.cc
 
-Copyright (C) 2015-2018
+Copyright (C) 2015-2019
 
 Author: Antonin Portelli <antonin.portelli@me.com>