Merge branch 'feature/hadrons' of github.com:paboyle/Grid into feature/hadrons

extras/Hadrons/Modules/MGauge/StochEm.cc

@@ -57,9 +57,11 @@ std::vector<std::string> TStochEm::getOutput(void)
 // setup ///////////////////////////////////////////////////////////////////////
 void TStochEm::setup(void)
 {
+    create_weight = false; 
     if (!env().hasCreatedObject("_" + getName() + "_weight"))
     {
         envCacheLat(EmComp, "_" + getName() + "_weight");
+        create_weight = true;
     }
     envCreateLat(EmField, getName());
 }
@@ -67,13 +69,13 @@ void TStochEm::setup(void)
 // execution ///////////////////////////////////////////////////////////////////
 void TStochEm::execute(void)
 {
-    LOG(Message) << "Generating stochatic EM potential..." << std::endl;
+    LOG(Message) << "Generating stochastic EM potential..." << std::endl;
 
     PhotonR photon(par().gauge, par().zmScheme);
     auto    &a = envGet(EmField, getName());
     auto    &w = envGet(EmComp, "_" + getName() + "_weight");
     
-    if (!env().hasCreatedObject("_" + getName() + "_weight"))
+    if (create_weight)
     {
         LOG(Message) << "Caching stochatic EM potential weight (gauge: "
                      << par().gauge << ", zero-mode scheme: "

extras/Hadrons/Modules/MGauge/StochEm.hpp

@@ -60,6 +60,8 @@ public:
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
+private:
+    bool create_weight;
 protected:
     // setup
     virtual void setup(void);

extras/Hadrons/Modules/MSource/SeqConserved.hpp

@@ -7,7 +7,8 @@ Source file: extras/Hadrons/Modules/MSource/SeqConserved.hpp
 Copyright (C) 2015-2018
 
 Author: Antonin Portelli <antonin.portelli@me.com>
-Author: Lanny91 <andrew.lawson@gmail.com>
+Author: Andrew Lawson    <andrew.lawson1991@gmail.com>
+Author: Vera Guelpers    <v.m.guelpers@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
@@ -38,9 +39,11 @@ BEGIN_HADRONS_NAMESPACE
 
 /*
  
- Sequential source
+ Sequential source with insertion of conserved current. 
+ Additionally optional insertion of a photon field A_\mu(x).
  -----------------------------
- * src_x = q_x * theta(x_3 - tA) * theta(tB - x_3) * J_mu * exp(i x.mom)
+ * src_x = sum_{mu=mu_min}^{mu_max} 
+     q_x * theta(x_3 - tA) * theta(tB - x_3) * J_mu * exp(i x.mom) (* A_\mu(x))
  
  * options:
  - q: input propagator (string)
@@ -48,8 +51,10 @@ BEGIN_HADRONS_NAMESPACE
  - tA: begin timeslice (integer)
  - tB: end timesilce (integer)
  - curr_type: type of conserved current to insert (Current)
- - mu: Lorentz index of current to insert (integer)
+ - mu_min: begin Lorentz Index (integer)
+ - mu_max: end Lorentz Index (integer)
  - mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0.")
+ - photon: optional photon field (string)
  
  */
 
@@ -67,8 +72,10 @@ public:
                                     unsigned int, tA,
                                     unsigned int, tB,
                                     Current,      curr_type,
-                                    unsigned int, mu,
+                                    unsigned int, mu_min,
-                                    std::string,  mom);
+                                    unsigned int, mu_max,
+                                    std::string,  mom,
+                                    std::string,  photon);
 };
 
 template <typename FImpl>
@@ -76,6 +83,8 @@ class TSeqConserved: public Module<SeqConservedPar>
 {
 public:
     FERM_TYPE_ALIASES(FImpl,);
+public:
+    typedef PhotonR::GaugeField     EmField;
 public:
     // constructor
     TSeqConserved(const std::string name);
@@ -89,10 +98,14 @@ protected:
     virtual void setup(void);
     // execution
     virtual void execute(void);
+private:
+    bool        SeqhasPhase_{false}; 
+    std::string SeqmomphName_;
 };
 
 MODULE_REGISTER_NS(SeqConserved, TSeqConserved<FIMPL>, MSource);
 
+
 /******************************************************************************
  *                      TSeqConserved implementation                          *
  ******************************************************************************/
@@ -100,6 +113,7 @@ MODULE_REGISTER_NS(SeqConserved, TSeqConserved<FIMPL>, MSource);
 template <typename FImpl>
 TSeqConserved<FImpl>::TSeqConserved(const std::string name)
 : Module<SeqConservedPar>(name)
+, SeqmomphName_ (name + "_Seqmomph")
 {}
 
 // dependencies/products ///////////////////////////////////////////////////////
@@ -107,6 +121,7 @@ template <typename FImpl>
 std::vector<std::string> TSeqConserved<FImpl>::getInput(void)
 {
     std::vector<std::string> in = {par().q, par().action};
+    if (!par().photon.empty()) in.push_back(par().photon);
         
     return in;
 }
@@ -125,6 +140,10 @@ void TSeqConserved<FImpl>::setup(void)
 {
     auto Ls_ = env().getObjectLs(par().action);
     envCreateLat(PropagatorField, getName(), Ls_);
+    envTmpLat(PropagatorField, "src_tmp");
+    envCacheLat(LatticeComplex, SeqmomphName_);
+    envTmpLat(LatticeComplex, "coor");
+    envTmpLat(LatticeComplex, "latt_compl");
 }
 
 // execution ///////////////////////////////////////////////////////////////////
@@ -134,27 +153,79 @@ void TSeqConserved<FImpl>::execute(void)
     if (par().tA == par().tB)
     {
         LOG(Message) << "Generating sequential source with conserved "
-                     << par().curr_type << " current insertion (mu = " 
+                     << par().curr_type << " current at " 
-                     << par().mu << ") at " << "t = " << par().tA << std::endl;
+		     << "t = " << par().tA << " summed over the indices " 
+		     << par().mu_min << " <= mu <= " << par().mu_max 
+		     << std::endl;
     }
     else
     {
         LOG(Message) << "Generating sequential source with conserved "
-                     << par().curr_type << " current insertion (mu = " 
+                     << par().curr_type << " current for " 
-                     << par().mu << ") for " << par().tA << " <= t <= " 
+                     << par().tA << " <= t <= " 
-                     << par().tB << std::endl;
+                     << par().tB << " summed over the indices " 
+		     << par().mu_min << " <= mu <= " << par().mu_max
+	             << std::endl;
     }
     auto &src = envGet(PropagatorField, getName());
+    envGetTmp(PropagatorField, src_tmp);
+    src_tmp = src;
     auto &q   = envGet(PropagatorField, par().q);
     auto &mat = envGet(FMat, par().action);
+    envGetTmp(LatticeComplex, latt_compl);
 
+    src = zero;
+
+    //exp(ipx)
+    auto &mom_phase = envGet(LatticeComplex, SeqmomphName_);
+    if (!SeqhasPhase_)
+    {    
         std::vector<Real> mom = strToVec<Real>(par().mom);
-    mat.SeqConservedCurrent(q, src, par().curr_type, par().mu, 
+        mom_phase = zero;
-                            mom, par().tA, par().tB);
+        Complex           i(0.0,1.0);
+        envGetTmp(LatticeComplex, coor);
+        for(unsigned int mu = 0; mu < env().getNd(); mu++)
+        {
+            LatticeCoordinate(coor, mu);
+            mom_phase = mom_phase + (mom[mu]/env().getGrid()->_fdimensions[mu])*coor;
         }
+        mom_phase = exp((Real)(2*M_PI)*i*mom_phase);
+        SeqhasPhase_ = true;
+    }
+    LOG(Message) << "Inserting momentum " << strToVec<Real>(par().mom) << std::endl;
+
+
+
+    if (!par().photon.empty())    	
+    {
+	 LOG(Message) << "Inserting the stochastic photon field " << par().photon << std::endl;
+    }
+
+    for(unsigned int mu=par().mu_min;mu<=par().mu_max;mu++)
+    {
+        if (!par().photon.empty())    	
+        {
+	    //Get the stochastic photon field, if required
+            auto &stoch_photon = envGet(EmField,  par().photon);
+    	    latt_compl =  PeekIndex<LorentzIndex>(stoch_photon, mu) * mom_phase;
+        }
+        else
+        {
+            latt_compl = mom_phase;
+        } 
+
+    	mat.SeqConservedCurrent(q, src_tmp, par().curr_type, mu, 
+                             par().tA, par().tB, latt_compl);
+	src += src_tmp;
+
+    }	
+
+ 
+}
+
 
 END_MODULE_NAMESPACE
 
 END_HADRONS_NAMESPACE
 
-#endif // Hadrons_SeqConserved_hpp_
+#endif // Hadrons_MSource_SeqConserved_hpp_

lib/qcd/action/fermion/FermionOperator.h

@@ -125,9 +125,9 @@ namespace Grid {
                                        PropagatorField &q_out,
                                        Current curr_type,
                                        unsigned int mu,
-                                       std::vector<Real> mom,
                                        unsigned int tmin, 
-                                       unsigned int tmax)=0;
+                                       unsigned int tmax,
+				       ComplexField &lattice_cmplx)=0;
     };
 
   }

lib/qcd/action/fermion/ImprovedStaggeredFermion.cc

@@ -411,13 +411,15 @@ void ImprovedStaggeredFermion<Impl>::SeqConservedCurrent(PropagatorField &q_in,
                                               PropagatorField &q_out,
                                               Current curr_type,
                                               unsigned int mu,
-                                                         std::vector<Real> mom,
                                               unsigned int tmin, 
-                                                         unsigned int tmax)
+                                              unsigned int tmax,
+					      ComplexField &lattice_cmplx)
 {
     assert(0);
+
 }
 
+
 FermOpStaggeredTemplateInstantiate(ImprovedStaggeredFermion);
 
   //AdjointFermOpTemplateInstantiate(ImprovedStaggeredFermion);

lib/qcd/action/fermion/ImprovedStaggeredFermion.h

@@ -170,9 +170,9 @@ class ImprovedStaggeredFermion : public StaggeredKernels<Impl>, public ImprovedS
                            PropagatorField &q_out,
                            Current curr_type, 
                            unsigned int mu,
-                           std::vector<Real> mom,
                            unsigned int tmin, 
-                           unsigned int tmax);
+                           unsigned int tmax,
+			   ComplexField &lattice_cmplx);
 };
 
 typedef ImprovedStaggeredFermion<StaggeredImplF> ImprovedStaggeredFermionF;

lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc

@@ -423,11 +423,12 @@ void ImprovedStaggeredFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in
                                               PropagatorField &q_out,
                                               Current curr_type,
                                               unsigned int mu,
-                                                          std::vector<Real> mom,
                                               unsigned int tmin, 
-                                                          unsigned int tmax)
+                                              unsigned int tmax,
+					      ComplexField &lattice_cmplx)
 {
     assert(0);
+
 }
 
 FermOpStaggeredTemplateInstantiate(ImprovedStaggeredFermion5D);

lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h

@@ -182,9 +182,9 @@ namespace QCD {
                              PropagatorField &q_out,
                              Current curr_type, 
                              unsigned int mu,
-                             std::vector<Real> mom,
                              unsigned int tmin, 
-                             unsigned int tmax);
+                             unsigned int tmax,
+                 	     ComplexField &lattice_cmplx);
   };
 
 }}

lib/qcd/action/fermion/WilsonFermion.cc

@@ -407,40 +407,30 @@ void WilsonFermion<Impl>::ContractConservedCurrent(PropagatorField &q_in_1,
     }
 }
 
+
 template <class Impl>
 void WilsonFermion<Impl>::SeqConservedCurrent(PropagatorField &q_in, 
                                               PropagatorField &q_out,
                                               Current curr_type,
                                               unsigned int mu,
-                                              std::vector<Real> mom,
                                               unsigned int tmin, 
-                                              unsigned int tmax)
+                                              unsigned int tmax,
+					      ComplexField &lattice_cmplx)
 {
     conformable(_grid, q_in._grid);
     conformable(_grid, q_out._grid);
-    Lattice<iSinglet<Simd>> ph(_grid), coor(_grid);
-    ComplexD i(0.0,1.0);
     PropagatorField tmpFwd(_grid), tmpBwd(_grid), tmp(_grid);
     unsigned int tshift = (mu == Tp) ? 1 : 0;
     unsigned int LLt    = GridDefaultLatt()[Tp];
 
-    // Momentum projection
-    ph = zero;
-    for(unsigned int mu = 0; mu < Nd - 1; mu++)
-    {
-        LatticeCoordinate(coor, mu);
-        ph = ph + mom[mu]*coor*((1./(_grid->_fdimensions[mu])));
-    }
-    ph = exp((RealD)(2*M_PI)*i*ph);
-
     q_out = zero;
     LatticeInteger coords(_grid);
     LatticeCoordinate(coords, Tp);
 
     // Need q(x + mu) and q(x - mu).
     tmp = Cshift(q_in, mu, 1);
-    tmpFwd = tmp*ph;
+    tmpFwd = tmp*lattice_cmplx;
-    tmp = ph*q_in;
+    tmp = lattice_cmplx*q_in;
     tmpBwd = Cshift(tmp, mu, -1);
 
     parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU)
@@ -475,6 +465,8 @@ void WilsonFermion<Impl>::SeqConservedCurrent(PropagatorField &q_in,
                                                 Umu, sU, mu, t_mask);
         }
     }
+
+
 }
 
 FermOpTemplateInstantiate(WilsonFermion);

lib/qcd/action/fermion/WilsonFermion.h

@@ -178,9 +178,9 @@ class WilsonFermion : public WilsonKernels<Impl>, public WilsonFermionStatic {
                              PropagatorField &q_out,
                              Current curr_type, 
                              unsigned int mu,
-                           std::vector<Real> mom,
                              unsigned int tmin, 
-                           unsigned int tmax);
+                             unsigned int tmax,
+			     ComplexField &lattice_cmplx);
 };
 
 typedef WilsonFermion<WilsonImplF> WilsonFermionF;

lib/qcd/action/fermion/WilsonFermion5D.cc

@@ -779,45 +779,39 @@ void WilsonFermion5D<Impl>::ContractConservedCurrent(PropagatorField &q_in_1,
 }
 
 
+
 template <class Impl>
 void WilsonFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in, 
                                                 PropagatorField &q_out,
                                                 Current curr_type, 
                                                 unsigned int mu,
-                                                std::vector<Real> mom,
                                                 unsigned int tmin, 
-                                                unsigned int tmax)
+                                                unsigned int tmax,
+						ComplexField &lattice_cmplx)
 {
     conformable(q_in._grid, FermionGrid());
     conformable(q_in._grid, q_out._grid);
-    Lattice<iSinglet<Simd>> ph(FermionGrid()), coor(FermionGrid());
+    PropagatorField tmp(GaugeGrid()),tmp2(GaugeGrid());
-    PropagatorField tmpFwd(FermionGrid()), tmpBwd(FermionGrid()),
-                    tmp(FermionGrid());
-    ComplexD i(0.0, 1.0);
     unsigned int tshift = (mu == Tp) ? 1 : 0;
     unsigned int LLs = q_in._grid->_rdimensions[0];
     unsigned int LLt    = GridDefaultLatt()[Tp];
 
-    // Momentum projection.
-    ph = zero;
-    for(unsigned int nu = 0; nu < Nd - 1; nu++)
-    {
-        // Shift coordinate lattice index by 1 to account for 5th dimension.
-        LatticeCoordinate(coor, nu + 1);
-        ph = ph + mom[nu]*coor*((1./(_FourDimGrid->_fdimensions[nu])));
-    }
-    ph = exp((RealD)(2*M_PI)*i*ph);
-
     q_out = zero;
     LatticeInteger coords(_FourDimGrid);
     LatticeCoordinate(coords, Tp);
 
-    // Need q(x + mu, s) and q(x - mu, s). 5D lattice so shift 4D coordinate mu
+
-    // by one.
+    for (unsigned int s = 0; s < LLs; ++s)
-    tmp = Cshift(q_in, mu + 1, 1);
+    {
-    tmpFwd = tmp*ph;
+        bool axial_sign = ((curr_type == Current::Axial) && (s < (LLs / 2)));
-    tmp = ph*q_in;
+	bool tadpole_sign = (curr_type == Current::Tadpole);
-    tmpBwd = Cshift(tmp, mu + 1, -1);
+	bool switch_sgn = tadpole_sign || axial_sign;
+
+
+        //forward direction: Need q(x + mu, s)*A(x)
+        ExtractSlice(tmp2, q_in, s, 0);  //q(x,s) 
+        tmp = Cshift(tmp2, mu, 1);	 //q(x+mu,s)
+        tmp2 = tmp*lattice_cmplx;	 //q(x+mu,s)*A(x)	
 
     	parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU)
     	{
@@ -829,42 +823,45 @@ void WilsonFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in,
 
             if (timeSlices > 0)
             {
-            unsigned int sF = sU * LLs;
+		unsigned int sF = sU * LLs + s;
-            for (unsigned int s = 0; s < LLs; ++s)
+                Kernels::SeqConservedCurrentSiteFwd(tmp2._odata[sU], 
-            {
-                bool axial_sign = ((curr_type == Current::Axial) && (s < (LLs / 2)));
-                Kernels::SeqConservedCurrentSiteFwd(tmpFwd._odata[sF], 
                                               q_out._odata[sF], Umu, sU,
-                                                    mu, t_mask, axial_sign);
+                                              mu, t_mask, switch_sgn);
-                ++sF;
-            }
             }
 
-        // Repeat for backward direction.
+        }
-        t_mask     = ((coords._odata[sU] >= (tmin + tshift)) && 
+
+        //backward direction: Need q(x - mu, s)*A(x-mu)
+        ExtractSlice(tmp2, q_in, s, 0);  //q(x,s)
+        tmp = lattice_cmplx*tmp2;	 //q(x,s)*A(x)
+        tmp2 = Cshift(tmp, mu, -1);	 //q(x-mu,s)*A(x-mu,s)
+
+    	parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU)
+    	{
+            vInteger  t_mask     = ((coords._odata[sU] >= (tmin + tshift)) && 
                    	  	    (coords._odata[sU] <= (tmax + tshift)));
 
 	    //if tmax = LLt-1 (last timeslice) include timeslice 0 if the time is shifted (mu=3)	
 	    unsigned int t0 = 0;
 	    if((tmax==LLt-1) && (tshift==1)) t_mask = (t_mask || (coords._odata[sU] == t0 ));
 
-        timeSlices = Reduce(t_mask);
+            Integer timeSlices = Reduce(t_mask);
 
             if (timeSlices > 0)
             {
-            unsigned int sF = sU * LLs;
+		unsigned int sF = sU * LLs + s; 
-            for (unsigned int s = 0; s < LLs; ++s)
+        	Kernels::SeqConservedCurrentSiteBwd(tmp2._odata[sU], 
-            {
-                bool axial_sign = ((curr_type == Current::Axial) && (s < (LLs / 2)));
-                Kernels::SeqConservedCurrentSiteBwd(tmpBwd._odata[sF], 
                                              q_out._odata[sF], Umu, sU,
                                              mu, t_mask, axial_sign);
-                ++sF;
             }
 	}
     }
 }
 
+
+
+
+
 FermOpTemplateInstantiate(WilsonFermion5D);
 GparityFermOpTemplateInstantiate(WilsonFermion5D);
   

lib/qcd/action/fermion/WilsonFermion5D.h

@@ -226,9 +226,9 @@ namespace QCD {
                              PropagatorField &q_out,
                              Current curr_type, 
                              unsigned int mu,
-                             std::vector<Real> mom,
                              unsigned int tmin, 
-                             unsigned int tmax);
+                             unsigned int tmax,
+			     ComplexField &lattice_cmplx);
   };
 
 }}

tests/hadrons/Test_QED.cc

@@ -0,0 +1,260 @@
+/*******************************************************************************
+ Grid physics library, www.github.com/paboyle/Grid
+ 
+ Source file: tests/hadrons/Test_hadrons_spectrum.cc
+ 
+ Copyright (C) 2015-2018
+ 
+ Author: Antonin Portelli <antonin.portelli@me.com>
+ Author: Vera Guelpers    <v.m.guelpers@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.
+ *******************************************************************************/
+
+#include <Grid/Hadrons/Application.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+
+int main(int argc, char *argv[])
+{
+    // initialization //////////////////////////////////////////////////////////
+    Grid_init(&argc, &argv);
+    HadronsLogError.Active(GridLogError.isActive());
+    HadronsLogWarning.Active(GridLogWarning.isActive());
+    HadronsLogMessage.Active(GridLogMessage.isActive());
+    HadronsLogIterative.Active(GridLogIterative.isActive());
+    HadronsLogDebug.Active(GridLogDebug.isActive());
+    LOG(Message) << "Grid initialized" << std::endl;
+
+
+    // run setup ///////////////////////////////////////////////////////////////
+    Application              application;
+    std::vector<std::string> flavour = {"h"}; //{"l", "s", "c1", "c2", "c3"};
+    std::vector<double>      mass    = {.2}; //{.01, .04, .2  , .25 , .3  };
+
+    unsigned int  nt    = GridDefaultLatt()[Tp];
+    
+    // global parameters
+    Application::GlobalPar globalPar;
+    globalPar.trajCounter.start = 1500;
+    globalPar.trajCounter.end   = 1520;
+    globalPar.trajCounter.step  = 20;
+    globalPar.seed              = "1 2 3 4";
+    application.setPar(globalPar);
+    // gauge field
+    application.createModule<MGauge::Unit>("gauge");
+    // pt source
+    MSource::Point::Par ptPar;
+    ptPar.position = "0 0 0 0";
+    application.createModule<MSource::Point>("pt", ptPar);
+    // sink
+    MSink::Point::Par sinkPar;
+    sinkPar.mom = "0 0 0";
+    application.createModule<MSink::ScalarPoint>("sink", sinkPar);
+    
+    // set fermion boundary conditions to be periodic space, antiperiodic time.
+    std::string boundary = "1 1 1 -1";
+
+    //stochastic photon field
+    MGauge::StochEm::Par photonPar;
+    photonPar.gauge = PhotonR::Gauge::feynman;
+    photonPar.zmScheme = PhotonR::ZmScheme::qedL;
+    application.createModule<MGauge::StochEm>("ph_field", photonPar);
+
+
+
+    for (unsigned int i = 0; i < flavour.size(); ++i)
+    {
+        // actions
+        MAction::DWF::Par actionPar;
+        actionPar.gauge = "gauge";
+        actionPar.Ls    = 8;
+        actionPar.M5    = 1.8;
+        actionPar.mass  = mass[i];
+        actionPar.boundary = boundary;
+        application.createModule<MAction::DWF>("DWF_" + flavour[i], actionPar);
+
+        
+        // solvers
+        MSolver::RBPrecCG::Par solverPar;
+        solverPar.action   = "DWF_" + flavour[i];
+        solverPar.residual = 1.0e-8;
+        application.createModule<MSolver::RBPrecCG>("CG_" + flavour[i],
+                                                    solverPar);
+        
+        // propagators
+        MFermion::GaugeProp::Par quarkPar;
+        quarkPar.solver = "CG_" + flavour[i];
+        quarkPar.source = "pt";
+        application.createModule<MFermion::GaugeProp>("Qpt_" + flavour[i],
+							 quarkPar);
+
+
+	//seq sources with tadpole insertion
+        MSource::SeqConserved::Par seqPar_T;
+        seqPar_T.q         = "Qpt_" + flavour[i] + "_5d";
+        seqPar_T.action    = "DWF_" + flavour[i];
+        seqPar_T.tA        = 0;
+        seqPar_T.tB        = nt-1;
+        seqPar_T.curr_type = Current::Tadpole;
+	seqPar_T.mu_min	   = 0;
+	seqPar_T.mu_max	   = 3;
+        seqPar_T.mom       = "0. 0. 0. 0.";
+        application.createModule<MSource::SeqConserved>("Qpt_" + flavour[i]
+							 + "_seq_T", seqPar_T);
+        // seq propagator with tadpole insertion
+        MFermion::GaugeProp::Par quarkPar_seq_T;
+        quarkPar_seq_T.solver = "CG_" + flavour[i];
+        quarkPar_seq_T.source = "Qpt_" + flavour[i] + "_seq_T";
+        application.createModule<MFermion::GaugeProp>("Qpt_" + flavour[i]
+							+ "_seq_T" + flavour[i],
+							quarkPar_seq_T);
+
+
+
+	//seq sources with conserved vector and photon insertion
+        MSource::SeqConserved::Par seqPar_V;
+        seqPar_V.q         = "Qpt_" + flavour[i] + "_5d";
+        seqPar_V.action    = "DWF_" + flavour[i];
+        seqPar_V.tA        = 0;
+        seqPar_V.tB        = nt-1;
+        seqPar_V.curr_type = Current::Vector;
+	seqPar_V.mu_min	   = 0;
+	seqPar_V.mu_max	   = 3;
+        seqPar_V.mom       = "0. 0. 0. 0.";
+	seqPar_V.photon	   = "ph_field";
+        application.createModule<MSource::SeqConserved>("Qpt_" + flavour[i] 
+						    + "_seq_V_ph", seqPar_V);
+        // seq propagator with conserved vector and photon insertion
+        MFermion::GaugeProp::Par quarkPar_seq_V;
+        quarkPar_seq_V.solver = "CG_" + flavour[i];
+        quarkPar_seq_V.source = "Qpt_" + flavour[i] + "_seq_V_ph";
+        application.createModule<MFermion::GaugeProp>("Qpt_" + flavour[i] 
+						+ "_seq_V_ph_" + flavour[i], 
+							quarkPar_seq_V);
+
+
+
+	//double seq sources with conserved vector and photon insertion
+	//(for self energy)
+        MSource::SeqConserved::Par seqPar_VV;
+        seqPar_VV.q         = "Qpt_" + flavour[i] + "_seq_V_ph_" 
+				+ flavour[i] + "_5d";
+        seqPar_VV.action    = "DWF_" + flavour[i];
+        seqPar_VV.tA        = 0;
+        seqPar_VV.tB        = nt-1;
+        seqPar_VV.curr_type = Current::Vector;
+	seqPar_VV.mu_min    = 0;
+	seqPar_VV.mu_max    = 3;
+        seqPar_VV.mom       = "0. 0. 0. 0.";
+	seqPar_VV.photon    = "ph_field";
+        application.createModule<MSource::SeqConserved>("Qpt_" + flavour[i] 
+						+ "_seq_V_ph" + flavour[i] 
+						+ "_seq_V_ph", seqPar_VV);
+        //double seq propagator with conserved vector and photon insertion
+        MFermion::GaugeProp::Par quarkPar_seq_VV;
+        quarkPar_seq_VV.solver = "CG_" + flavour[i];
+        quarkPar_seq_VV.source = "Qpt_" + flavour[i] + "_seq_V_ph" 
+						+ flavour[i] + "_seq_V_ph";
+        application.createModule<MFermion::GaugeProp>("Qpt_" + flavour[i] 
+						+ "_seq_V_ph_" + flavour[i] 
+						+ "_seq_V_ph_" + flavour[i], 
+							quarkPar_seq_VV);
+
+
+
+    }
+    for (unsigned int i = 0; i < flavour.size(); ++i)
+    for (unsigned int j = i; j < flavour.size(); ++j)
+    {
+        //2pt function contraction
+        MContraction::Meson::Par mesPar;
+        mesPar.output  = "QED/pt_" + flavour[i] + flavour[j];
+        mesPar.q1      = "Qpt_" + flavour[i];
+        mesPar.q2      = "Qpt_" + flavour[j];
+        mesPar.gammas  = "<Gamma5 Gamma5>";
+        mesPar.sink    = "sink";
+        application.createModule<MContraction::Meson>("meson_pt_"
+                                                      + flavour[i] + flavour[j],
+                                                      mesPar);
+
+
+
+        //tadpole contraction
+        MContraction::Meson::Par mesPar_seq_T;
+        mesPar_seq_T.output  = "QED/tadpole_pt_" + flavour[i] + "_T_" 
+					+ flavour[i] + "__" + flavour[j];
+        mesPar_seq_T.q1      = "Qpt_" + flavour[i] + "_seq_T" + flavour[i];
+        mesPar_seq_T.q2      = "Qpt_" + flavour[j];
+        mesPar_seq_T.gammas  = "<Gamma5 Gamma5>";
+        mesPar_seq_T.sink    = "sink";
+        application.createModule<MContraction::Meson>("meson_tadpole_pt_" +
+                                                      flavour[i] + "_seq_T" 
+						      + flavour[i] + flavour[j],
+                                                      mesPar_seq_T);
+
+
+
+        //photon exchange contraction
+        MContraction::Meson::Par mesPar_seq_E;
+        mesPar_seq_E.output  = "QED/exchange_pt_" + flavour[i] + "_V_ph_" 
+				+ flavour[i] + "__" + flavour[j] + "_V_ph_"
+				+ flavour[j];
+        mesPar_seq_E.q1      = "Qpt_" + flavour[i] + "_seq_V_ph_" + flavour[i];
+        mesPar_seq_E.q2      = "Qpt_" + flavour[j] + "_seq_V_ph_" + flavour[j];
+        mesPar_seq_E.gammas  = "<Gamma5 Gamma5>";
+        mesPar_seq_E.sink    = "sink";
+        application.createModule<MContraction::Meson>("meson_exchange_pt_" 
+					+ flavour[i] + "_seq_V_ph_" + flavour[i] 
+					+ flavour[j] + "_seq_V_ph_" + flavour[j],
+                                                      mesPar_seq_E);
+
+
+
+        //self energy contraction
+        MContraction::Meson::Par mesPar_seq_S;
+        mesPar_seq_S.output  = "QED/selfenergy_pt_" + flavour[i] + "_V_ph_" 
+				+ flavour[i] + "_V_ph_" + flavour[i] + "__" 
+				+  flavour[j];
+        mesPar_seq_S.q1      = "Qpt_" + flavour[i] + "_seq_V_ph_" + flavour[i] 
+				+ "_seq_V_ph_" + flavour[i];
+        mesPar_seq_S.q2      = "Qpt_" + flavour[j];
+        mesPar_seq_S.gammas  = "<Gamma5 Gamma5>";
+        mesPar_seq_S.sink    = "sink";
+        application.createModule<MContraction::Meson>("meson_selfenergy_pt_" 
+						    + flavour[i] + "_seq_V_ph_" 
+						    + flavour[i] + "_seq_V_ph_" 
+						    + flavour[i] + flavour[j],
+                                                       mesPar_seq_S);
+
+    }
+
+
+
+    
+    // execution
+    application.saveParameterFile("QED.xml");
+    application.run();
+    
+    // epilogue
+    LOG(Message) << "Grid is finalizing now" << std::endl;
+    Grid_finalize();
+    
+    return EXIT_SUCCESS;
+}

tests/hadrons/Test_hadrons.hpp

@@ -263,7 +263,8 @@ inline void makeConservedSequentialSource(Application &application,
         seqPar.tA        = tS;
         seqPar.tB        = tS;
         seqPar.curr_type = curr;
-        seqPar.mu        = mu;
+        seqPar.mu_min    = mu;
+        seqPar.mu_min    = mu;
         seqPar.mom       = mom;
         application.createModule<MSource::SeqConserved>(srcName, seqPar);
     }