Insertion of photon field in seqential conserved current

extras/Hadrons/Modules.hpp

@@ -28,7 +28,6 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #include <Grid/Hadrons/Modules/MContraction/WardIdentity.hpp>
 #include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp>
 #include <Grid/Hadrons/Modules/MFermion/GaugeProp.hpp>
-#include <Grid/Hadrons/Modules/MSource/SeqConservedSummed.hpp>
 #include <Grid/Hadrons/Modules/MSource/SeqGamma.hpp>
 #include <Grid/Hadrons/Modules/MSource/Point.hpp>
 #include <Grid/Hadrons/Modules/MSource/Wall.hpp>

extras/Hadrons/Modules/MSource/SeqConserved.hpp

@@ -2,12 +2,13 @@
 
 Grid physics library, www.github.com/paboyle/Grid 
 
-Source file: extras/Hadrons/Modules/MSource/SeqConserved.hpp
+Source file: extras/Hadrons/Modules/MContraction/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
@@ -34,13 +35,17 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #include <Grid/Hadrons/Module.hpp>
 #include <Grid/Hadrons/ModuleFactory.hpp>
 
+#include <Grid/qcd/action/gauge/GaugeImplTypes.h>
+
 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 +53,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 +74,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 +85,8 @@ class TSeqConserved: public Module<SeqConservedPar>
 {
 public:
     FERM_TYPE_ALIASES(FImpl,);
+public:
+    typedef PhotonR::GaugeField     EmField;
 public:
     // constructor
     TSeqConserved(const std::string name);
@@ -93,6 +104,7 @@ protected:
 
 MODULE_REGISTER_NS(SeqConserved, TSeqConserved<FIMPL>, MSource);
 
+
 /******************************************************************************
  *                      TSeqConserved implementation                          *
  ******************************************************************************/
@@ -107,7 +119,8 @@ 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;
 }
 
@@ -116,7 +129,7 @@ std::vector<std::string> TSeqConserved<FImpl>::getOutput(void)
 {
     std::vector<std::string> out = {getName()};
     
-    return out;
+   return out;
 }
 
 // setup ///////////////////////////////////////////////////////////////////////
@@ -125,6 +138,10 @@ void TSeqConserved<FImpl>::setup(void)
 {
     auto Ls_ = env().getObjectLs(par().action);
     envCreateLat(PropagatorField, getName(), Ls_);
+    envTmpLat(PropagatorField, "src_tmp");
+    envTmpLat(LatticeComplex, "mom_phase");
+    envTmpLat(LatticeComplex, "coor");
+    envTmpLat(LatticeComplex, "latt_compl");
 }
 
 // execution ///////////////////////////////////////////////////////////////////
@@ -134,27 +151,71 @@ 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 << 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;
     }
     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, mom_phase);
+    envGetTmp(LatticeComplex, coor);
+    envGetTmp(LatticeComplex, latt_compl);
 
+    src = zero;
+
+    //exp(ipx)
     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);
+    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);
+    LOG(Message) << "Inserting momentum " << 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_

extras/Hadrons/Modules/MSource/SeqConservedSummed.hpp

@@ -1,169 +0,0 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: extras/Hadrons/Modules/MContraction/SeqConservedSummed.hpp
-
-Copyright (C) 2015-2018
-
-Author: Antonin Portelli <antonin.portelli@me.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
-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_MSource_SeqConservedSummed_hpp_
-#define Hadrons_MSource_SeqConservedSummed_hpp_
-
-#include <Grid/Hadrons/Global.hpp>
-#include <Grid/Hadrons/Module.hpp>
-#include <Grid/Hadrons/ModuleFactory.hpp>
-
-BEGIN_HADRONS_NAMESPACE
-
-/*
- 
- Sequential source summed over the Lorentz index of current
- -----------------------------
- * src_x = sum_mu   q_x * theta(x_3 - tA) * theta(tB - x_3) * J_mu * exp(i x.mom)
- 
- * options:
- - q: input propagator (string)
- - action: fermion action used for propagator q (string)
- - tA: begin timeslice (integer)
- - tB: end timesilce (integer)
- - curr_type: type of conserved current to insert (Current)
- - mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0.")
- 
- */
-
-/******************************************************************************
- *                              SeqConservedSummed                                  *
- ******************************************************************************/
-BEGIN_MODULE_NAMESPACE(MSource)
-
-class SeqConservedSummedPar: Serializable
-{
-public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(SeqConservedSummedPar,
-                                    std::string,  q,
-                                    std::string,  action,
-                                    unsigned int, tA,
-                                    unsigned int, tB,
-                                    Current,      curr_type,
-                                    std::string,  mom);
-};
-
-template <typename FImpl>
-class TSeqConservedSummed: public Module<SeqConservedSummedPar>
-{
-public:
-    FERM_TYPE_ALIASES(FImpl,);
-public:
-    // constructor
-    TSeqConservedSummed(const std::string name);
-    // destructor
-    virtual ~TSeqConservedSummed(void) = default;
-    // 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_NS(SeqConservedSummed, TSeqConservedSummed<FIMPL>, MSource);
-
-/******************************************************************************
- *                      TSeqConservedSummed implementation                          *
- ******************************************************************************/
-// constructor /////////////////////////////////////////////////////////////////
-template <typename FImpl>
-TSeqConservedSummed<FImpl>::TSeqConservedSummed(const std::string name)
-: Module<SeqConservedSummedPar>(name)
-{}
-
-// dependencies/products ///////////////////////////////////////////////////////
-template <typename FImpl>
-std::vector<std::string> TSeqConservedSummed<FImpl>::getInput(void)
-{
-    std::vector<std::string> in = {par().q, par().action};
-    
-    return in;
-}
-
-template <typename FImpl>
-std::vector<std::string> TSeqConservedSummed<FImpl>::getOutput(void)
-{
-    std::vector<std::string> out = {getName()};
-    
-    return out;
-}
-
-// setup ///////////////////////////////////////////////////////////////////////
-template <typename FImpl>
-void TSeqConservedSummed<FImpl>::setup(void)
-{
-    auto Ls_ = env().getObjectLs(par().action);
-    envCreateLat(PropagatorField, getName(), Ls_);
-    envTmpLat(PropagatorField, "src_tmp");
-}
-
-// execution ///////////////////////////////////////////////////////////////////
-template <typename FImpl>
-void TSeqConservedSummed<FImpl>::execute(void)
-{
-    if (par().tA == par().tB)
-    {
-        LOG(Message) << "Generating sequential source with conserved "
-                     << par().curr_type << " current insertion summed over mu at " 
-		     << "t = " << par().tA << std::endl;
-    }
-    else
-    {
-        LOG(Message) << "Generating sequential source with conserved "
-                     << par().curr_type << " current insertion summed over mu for " 
-                     << par().tA << " <= t <= " 
-                     << par().tB << 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);
-
-    std::vector<Real> mom = strToVec<Real>(par().mom);
-    src = zero;
-    for(int mu=0;mu<=3;mu++)
-    {
-    	mat.SeqConservedCurrent(q, src_tmp, par().curr_type, mu, 
-                            mom, par().tA, par().tB);
-	src += src_tmp;
-
-    }	
-
-}
-
-END_MODULE_NAMESPACE
-
-END_HADRONS_NAMESPACE
-
-#endif // Hadrons_SeqConservedSummed_hpp_

extras/Hadrons/modules.inc

@@ -20,7 +20,6 @@ modules_hpp =\
   Modules/MContraction/WardIdentity.hpp \
   Modules/MContraction/WeakHamiltonianEye.hpp \
   Modules/MFermion/GaugeProp.hpp \
-  Modules/MSource/SeqConservedSummed.hpp \
   Modules/MSource/SeqGamma.hpp \
   Modules/MSource/Point.hpp \
   Modules/MSource/Wall.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,
+				       Lattice<iSinglet<Simd>> &lattice_cmplx)=0;
     };
 
   }

lib/qcd/action/fermion/ImprovedStaggeredFermion.cc

@@ -407,17 +407,19 @@ void ImprovedStaggeredFermion<Impl>::ContractConservedCurrent(PropagatorField &q
 }
 
 template <class Impl>
-void ImprovedStaggeredFermion<Impl>::SeqConservedCurrent(PropagatorField &q_in,
+void ImprovedStaggeredFermion<Impl>::SeqConservedCurrent(PropagatorField &q_in, 
-                                                         PropagatorField &q_out,
+                                              PropagatorField &q_out,
-                                                         Current curr_type,
+                                              Current curr_type,
-                                                         unsigned int mu, 
+                                              unsigned int mu,
-                                                         std::vector<Real> mom,
+                                              unsigned int tmin, 
-                                                         unsigned int tmin,
+                                              unsigned int tmax,
-                                                         unsigned int tmax)
+					      Lattice<iSinglet<Simd>> &lattice_cmplx)
 {
     assert(0);
+
 }
 
+
 FermOpStaggeredTemplateInstantiate(ImprovedStaggeredFermion);
 
   //AdjointFermOpTemplateInstantiate(ImprovedStaggeredFermion);

lib/qcd/action/fermion/ImprovedStaggeredFermion.h

@@ -166,13 +166,13 @@ class ImprovedStaggeredFermion : public StaggeredKernels<Impl>, public ImprovedS
                                 PropagatorField &q_out,
                                 Current curr_type,
                                 unsigned int mu);
-  void SeqConservedCurrent(PropagatorField &q_in,
+  void SeqConservedCurrent(PropagatorField &q_in, 
                            PropagatorField &q_out,
-                           Current curr_type,
+                           Current curr_type, 
-                           unsigned int mu, 
+                           unsigned int mu,
-                           std::vector<Real> mom,
+                           unsigned int tmin, 
-                           unsigned int tmin,
+                           unsigned int tmax,
-                           unsigned int tmax);
+			   Lattice<iSinglet<Simd>> &lattice_cmplx);
 };
 
 typedef ImprovedStaggeredFermion<StaggeredImplF> ImprovedStaggeredFermionF;

lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc

@@ -419,15 +419,16 @@ void ImprovedStaggeredFermion5D<Impl>::ContractConservedCurrent(PropagatorField
 }
 
 template <class Impl>
-void ImprovedStaggeredFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in,
+void ImprovedStaggeredFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in, 
-                                                          PropagatorField &q_out,
+                                              PropagatorField &q_out,
-                                                          Current curr_type,
+                                              Current curr_type,
-                                                          unsigned int mu, 
+                                              unsigned int mu,
-                                                          std::vector<Real> mom,
+                                              unsigned int tmin, 
-                                                          unsigned int tmin,
+                                              unsigned int tmax,
-                                                          unsigned int tmax)
+					      Lattice<iSinglet<Simd>> &lattice_cmplx)
 {
     assert(0);
+
 }
 
 FermOpStaggeredTemplateInstantiate(ImprovedStaggeredFermion5D);

lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h

@@ -178,13 +178,13 @@ namespace QCD {
                                   PropagatorField &q_out,
                                   Current curr_type,
                                   unsigned int mu);
-    void SeqConservedCurrent(PropagatorField &q_in,
+    void SeqConservedCurrent(PropagatorField &q_in, 
                              PropagatorField &q_out,
-                             Current curr_type,
+                             Current curr_type, 
-                             unsigned int mu, 
+                             unsigned int mu,
-                             std::vector<Real> mom,
+                             unsigned int tmin, 
-                             unsigned int tmin,
+                             unsigned int tmax,
-                             unsigned int tmax);
+                 	     Lattice<iSinglet<Simd>> &lattice_cmplx);
   };
 
 }}

lib/qcd/action/fermion/WilsonFermion.cc

@@ -381,40 +381,100 @@ 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)
+//{
+//    conformable(_grid, q_in._grid);
+//    conformable(_grid, q_out._grid);
+//    Lattice<iSinglet<Simd>> ph(_grid), coor(_grid);
+//    Complex 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((Real)(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;
+//    tmp = ph*q_in;
+//    tmpBwd = Cshift(tmp, mu, -1);
+//
+//    parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU)
+//    {
+//        // Compute the sequential conserved current insertion only if our simd
+//        // object contains a timeslice we need.
+//        vInteger t_mask   = ((coords._odata[sU] >= tmin) &&
+//                             (coords._odata[sU] <= tmax));
+//        Integer timeSlices = Reduce(t_mask);
+//
+//        if (timeSlices > 0)
+//        {
+//            Kernels::SeqConservedCurrentSiteFwd(tmpFwd._odata[sU], 
+//                                                q_out._odata[sU], 
+//                                                Umu, sU, mu, t_mask);
+//        }
+//
+//        // Repeat for backward direction.
+//        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);
+//
+//        if (timeSlices > 0)
+//        {
+//            Kernels::SeqConservedCurrentSiteBwd(tmpBwd._odata[sU], 
+//                                                q_out._odata[sU], 
+//                                                Umu, sU, mu, t_mask);
+//        }
+//    }
+//}
+
 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,
+					      Lattice<iSinglet<Simd>> &lattice_cmplx)
 {
     conformable(_grid, q_in._grid);
     conformable(_grid, q_out._grid);
-    Lattice<iSinglet<Simd>> ph(_grid), coor(_grid);
     Complex 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((Real)(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)
@@ -449,6 +509,8 @@ void WilsonFermion<Impl>::SeqConservedCurrent(PropagatorField &q_in,
                                                 Umu, sU, mu, t_mask);
         }
     }
+
+
 }
 
 FermOpTemplateInstantiate(WilsonFermion);

lib/qcd/action/fermion/WilsonFermion.h

@@ -155,13 +155,13 @@ class WilsonFermion : public WilsonKernels<Impl>, public WilsonFermionStatic {
                                 PropagatorField &q_out,
                                 Current curr_type,
                                 unsigned int mu);
-  void SeqConservedCurrent(PropagatorField &q_in,
+    void SeqConservedCurrent(PropagatorField &q_in, 
-                           PropagatorField &q_out,
+                             PropagatorField &q_out,
-                           Current curr_type,
+                             Current curr_type, 
-                           unsigned int mu, 
+                             unsigned int mu,
-                           std::vector<Real> mom,
+                             unsigned int tmin, 
-                           unsigned int tmin,
+                             unsigned int tmax,
-                           unsigned int tmax);
+			     Lattice<iSinglet<Simd>> &lattice_cmplx);
 };
 
 typedef WilsonFermion<WilsonImplF> WilsonFermionF;

lib/qcd/action/fermion/WilsonFermion5D.cc

@@ -779,18 +779,110 @@ 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)
+//{
+//    conformable(q_in._grid, FermionGrid());
+//    conformable(q_in._grid, q_out._grid);
+//    Lattice<iSinglet<Simd>> ph(FermionGrid()), coor(FermionGrid());
+//    PropagatorField tmpFwd(FermionGrid()), tmpBwd(FermionGrid()),
+//                    tmp(FermionGrid());
+//    Complex 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((Real)(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.
+//    tmp = Cshift(q_in, mu + 1, 1);
+//    tmpFwd = tmp*ph;
+//    tmp = ph*q_in;
+//    tmpBwd = Cshift(tmp, mu + 1, -1);
+//
+//    parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU)
+//    {
+//        // Compute the sequential conserved current insertion only if our simd
+//        // object contains a timeslice we need.
+//        vInteger t_mask   = ((coords._odata[sU] >= tmin) &&
+//                             (coords._odata[sU] <= tmax));
+//        Integer timeSlices = Reduce(t_mask);
+//
+//        if (timeSlices > 0)
+//        {
+//            unsigned int sF = sU * LLs;
+//            for (unsigned int s = 0; s < LLs; ++s)
+//            {
+//                bool axial_sign = ((curr_type == Current::Axial) && (s < (LLs / 2)));
+//		bool tadpole_sign = (curr_type == Current::Tadpole);
+//		bool switch_sgn = tadpole_sign || axial_sign;
+//
+//                Kernels::SeqConservedCurrentSiteFwd(tmpFwd._odata[sF], 
+//                                                    q_out._odata[sF], Umu, sU,
+//                                                    mu, t_mask, switch_sgn);
+//                ++sF;
+//            }
+//        }
+//
+//        // Repeat for backward direction.
+//        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);
+//
+//        if (timeSlices > 0)
+//        {
+//            unsigned int sF = sU * LLs;
+//            for (unsigned int s = 0; s < LLs; ++s)
+//            {
+//                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;
+//            }
+//        }
+//    }
+//}
+
+
+
+
 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,
+						Lattice<iSinglet<Simd>> &lattice_cmplx)
 {
     conformable(q_in._grid, FermionGrid());
     conformable(q_in._grid, q_out._grid);
-    Lattice<iSinglet<Simd>> ph(FermionGrid()), coor(FermionGrid());
     PropagatorField tmpFwd(FermionGrid()), tmpBwd(FermionGrid()),
                     tmp(FermionGrid());
     Complex i(0.0, 1.0);
@@ -798,25 +890,26 @@ void WilsonFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in,
     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((Real)(2*M_PI)*i*ph);
-
     q_out = zero;
     LatticeInteger coords(_FourDimGrid);
     LatticeCoordinate(coords, Tp);
 
+
+    //QED: photon field is 4dim, but need a 5dim object to multiply to
+    //     DWF PropagatorField
+    Lattice<iSinglet<Simd>> lattice_cmplx_5d(FermionGrid());
+    for (unsigned int s = 0; s < LLs; ++s)
+    {   
+         InsertSlice(lattice_cmplx,lattice_cmplx_5d, s, 0);
+    }
+
+
+
     // Need q(x + mu, s) and q(x - mu, s). 5D lattice so shift 4D coordinate mu
     // by one.
     tmp = Cshift(q_in, mu + 1, 1);
-    tmpFwd = tmp*ph;
+    tmpFwd = tmp*lattice_cmplx_5d;
-    tmp = ph*q_in;
+    tmp = lattice_cmplx_5d*q_in;
     tmpBwd = Cshift(tmp, mu + 1, -1);
 
     parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU)
@@ -868,6 +961,10 @@ void WilsonFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in,
     }
 }
 
+
+
+
+
 FermOpTemplateInstantiate(WilsonFermion5D);
 GparityFermOpTemplateInstantiate(WilsonFermion5D);
   

lib/qcd/action/fermion/WilsonFermion5D.h

@@ -222,13 +222,13 @@ namespace QCD {
                                   PropagatorField &q_out,
                                   Current curr_type, 
                                   unsigned int mu);
-    void SeqConservedCurrent(PropagatorField &q_in,
+    void SeqConservedCurrent(PropagatorField &q_in, 
                              PropagatorField &q_out,
-                             Current curr_type,
+                             Current curr_type, 
                              unsigned int mu,
-                             std::vector<Real> mom,
+                             unsigned int tmin, 
-                             unsigned int tmin,
+                             unsigned int tmax,
-                             unsigned int tmax);
+			     Lattice<iSinglet<Simd>> &lattice_cmplx);
   };
 
 }}

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);
     }