diff --git a/extras/Hadrons/Modules/MGauge/StochEm.cc b/extras/Hadrons/Modules/MGauge/StochEm.cc index 21b7f626..cfbd28ef 100644 --- a/extras/Hadrons/Modules/MGauge/StochEm.cc +++ b/extras/Hadrons/Modules/MGauge/StochEm.cc @@ -57,9 +57,11 @@ std::vector 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: " diff --git a/extras/Hadrons/Modules/MGauge/StochEm.hpp b/extras/Hadrons/Modules/MGauge/StochEm.hpp index 87b70880..05162443 100644 --- a/extras/Hadrons/Modules/MGauge/StochEm.hpp +++ b/extras/Hadrons/Modules/MGauge/StochEm.hpp @@ -60,6 +60,8 @@ public: // dependency relation virtual std::vector getInput(void); virtual std::vector getOutput(void); +private: + bool create_weight; protected: // setup virtual void setup(void); diff --git a/extras/Hadrons/Modules/MSource/SeqConserved.hpp b/extras/Hadrons/Modules/MSource/SeqConserved.hpp index ee8d8d56..5c4a320b 100644 --- a/extras/Hadrons/Modules/MSource/SeqConserved.hpp +++ b/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 -Author: Lanny91 +Author: Andrew Lawson +Author: Vera Guelpers 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, - std::string, mom); + unsigned int, mu_min, + unsigned int, mu_max, + std::string, mom, + std::string, photon); }; template @@ -76,6 +83,8 @@ class TSeqConserved: public Module { 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, MSource); + /****************************************************************************** * TSeqConserved implementation * ******************************************************************************/ @@ -100,6 +113,7 @@ MODULE_REGISTER_NS(SeqConserved, TSeqConserved, MSource); template TSeqConserved::TSeqConserved(const std::string name) : Module(name) +, SeqmomphName_ (name + "_Seqmomph") {} // dependencies/products /////////////////////////////////////////////////////// @@ -107,7 +121,8 @@ template std::vector TSeqConserved::getInput(void) { std::vector in = {par().q, par().action}; - + if (!par().photon.empty()) in.push_back(par().photon); + return in; } @@ -116,7 +131,7 @@ std::vector TSeqConserved::getOutput(void) { std::vector out = {getName()}; - return out; + return out; } // setup /////////////////////////////////////////////////////////////////////// @@ -125,6 +140,10 @@ void TSeqConserved::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::execute(void) if (par().tA == par().tB) { LOG(Message) << "Generating sequential source with conserved " - << par().curr_type << " current insertion (mu = " - << par().mu << ") at " << "t = " << par().tA << std::endl; + << par().curr_type << " current at " + << "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().mu << ") for " << par().tA << " <= t <= " - << par().tB << std::endl; + << par().curr_type << " current for " + << par().tA << " <= t <= " + << 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); - std::vector mom = strToVec(par().mom); - mat.SeqConservedCurrent(q, src, par().curr_type, par().mu, - mom, par().tA, par().tB); + src = zero; + + //exp(ipx) + auto &mom_phase = envGet(LatticeComplex, SeqmomphName_); + if (!SeqhasPhase_) + { + std::vector mom = strToVec(par().mom); + mom_phase = zero; + 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(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(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_ diff --git a/lib/qcd/action/fermion/FermionOperator.h b/lib/qcd/action/fermion/FermionOperator.h index 1d395d53..9eb74baa 100644 --- a/lib/qcd/action/fermion/FermionOperator.h +++ b/lib/qcd/action/fermion/FermionOperator.h @@ -125,9 +125,9 @@ namespace Grid { PropagatorField &q_out, Current curr_type, unsigned int mu, - std::vector mom, unsigned int tmin, - unsigned int tmax)=0; + unsigned int tmax, + ComplexField &lattice_cmplx)=0; }; } diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc b/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc index 5005b1a2..ed07fb5b 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion.cc @@ -407,17 +407,19 @@ void ImprovedStaggeredFermion::ContractConservedCurrent(PropagatorField &q } template -void ImprovedStaggeredFermion::SeqConservedCurrent(PropagatorField &q_in, - PropagatorField &q_out, - Current curr_type, - unsigned int mu, - std::vector mom, - unsigned int tmin, - unsigned int tmax) +void ImprovedStaggeredFermion::SeqConservedCurrent(PropagatorField &q_in, + PropagatorField &q_out, + Current curr_type, + unsigned int mu, + unsigned int tmin, + unsigned int tmax, + ComplexField &lattice_cmplx) { assert(0); + } + FermOpStaggeredTemplateInstantiate(ImprovedStaggeredFermion); //AdjointFermOpTemplateInstantiate(ImprovedStaggeredFermion); diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion.h b/lib/qcd/action/fermion/ImprovedStaggeredFermion.h index 9d5270c6..d2426ef7 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion.h +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion.h @@ -166,13 +166,13 @@ class ImprovedStaggeredFermion : public StaggeredKernels, 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, - unsigned int mu, - std::vector mom, - unsigned int tmin, - unsigned int tmax); + Current curr_type, + unsigned int mu, + unsigned int tmin, + unsigned int tmax, + ComplexField &lattice_cmplx); }; typedef ImprovedStaggeredFermion ImprovedStaggeredFermionF; diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc index 7308ff30..3510b8f1 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc +++ b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc @@ -419,15 +419,16 @@ void ImprovedStaggeredFermion5D::ContractConservedCurrent(PropagatorField } template -void ImprovedStaggeredFermion5D::SeqConservedCurrent(PropagatorField &q_in, - PropagatorField &q_out, - Current curr_type, - unsigned int mu, - std::vector mom, - unsigned int tmin, - unsigned int tmax) +void ImprovedStaggeredFermion5D::SeqConservedCurrent(PropagatorField &q_in, + PropagatorField &q_out, + Current curr_type, + unsigned int mu, + unsigned int tmin, + unsigned int tmax, + ComplexField &lattice_cmplx) { assert(0); + } FermOpStaggeredTemplateInstantiate(ImprovedStaggeredFermion5D); diff --git a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h b/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h index 22fb9e7d..e9cf2013 100644 --- a/lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h +++ b/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, - unsigned int mu, - std::vector mom, - unsigned int tmin, - unsigned int tmax); + Current curr_type, + unsigned int mu, + unsigned int tmin, + unsigned int tmax, + ComplexField &lattice_cmplx); }; }} diff --git a/lib/qcd/action/fermion/WilsonFermion.cc b/lib/qcd/action/fermion/WilsonFermion.cc index dfaa6758..2d9cf22d 100644 --- a/lib/qcd/action/fermion/WilsonFermion.cc +++ b/lib/qcd/action/fermion/WilsonFermion.cc @@ -407,40 +407,30 @@ void WilsonFermion::ContractConservedCurrent(PropagatorField &q_in_1, } } + template void WilsonFermion::SeqConservedCurrent(PropagatorField &q_in, PropagatorField &q_out, Current curr_type, unsigned int mu, - std::vector mom, unsigned int tmin, - unsigned int tmax) + unsigned int tmax, + ComplexField &lattice_cmplx) { conformable(_grid, q_in._grid); conformable(_grid, q_out._grid); - Lattice> 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; - tmp = ph*q_in; + tmpFwd = tmp*lattice_cmplx; + 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::SeqConservedCurrent(PropagatorField &q_in, Umu, sU, mu, t_mask); } } + + } FermOpTemplateInstantiate(WilsonFermion); diff --git a/lib/qcd/action/fermion/WilsonFermion.h b/lib/qcd/action/fermion/WilsonFermion.h index c0db827e..ea25ed7f 100644 --- a/lib/qcd/action/fermion/WilsonFermion.h +++ b/lib/qcd/action/fermion/WilsonFermion.h @@ -174,13 +174,13 @@ class WilsonFermion : public WilsonKernels, public WilsonFermionStatic { PropagatorField &q_out, Current curr_type, unsigned int mu); - void SeqConservedCurrent(PropagatorField &q_in, - PropagatorField &q_out, - Current curr_type, - unsigned int mu, - std::vector mom, - unsigned int tmin, - unsigned int tmax); + void SeqConservedCurrent(PropagatorField &q_in, + PropagatorField &q_out, + Current curr_type, + unsigned int mu, + unsigned int tmin, + unsigned int tmax, + ComplexField &lattice_cmplx); }; typedef WilsonFermion WilsonFermionF; diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index 3e58fed6..6f82aad2 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -779,92 +779,89 @@ void WilsonFermion5D::ContractConservedCurrent(PropagatorField &q_in_1, } + template void WilsonFermion5D::SeqConservedCurrent(PropagatorField &q_in, PropagatorField &q_out, Current curr_type, unsigned int mu, - std::vector 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> ph(FermionGrid()), coor(FermionGrid()); - PropagatorField tmpFwd(FermionGrid()), tmpBwd(FermionGrid()), - tmp(FermionGrid()); - ComplexD i(0.0, 1.0); + PropagatorField tmp(GaugeGrid()),tmp2(GaugeGrid()); 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. - 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) + for (unsigned int s = 0; s < LLs; ++s) { - // 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); + bool axial_sign = ((curr_type == Current::Axial) && (s < (LLs / 2))); + bool tadpole_sign = (curr_type == Current::Tadpole); + bool switch_sgn = tadpole_sign || axial_sign; - if (timeSlices > 0) - { - unsigned int sF = sU * LLs; - for (unsigned int s = 0; s < LLs; ++s) + + //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) + { + // 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) { - 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); - ++sF; + unsigned int sF = sU * LLs + s; + Kernels::SeqConservedCurrentSiteFwd(tmp2._odata[sU], + q_out._odata[sF], Umu, sU, + mu, t_mask, switch_sgn); } + } - // Repeat for backward direction. - t_mask = ((coords._odata[sU] >= (tmin + tshift)) && - (coords._odata[sU] <= (tmax + 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) - //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 )); + parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU) + { + vInteger t_mask = ((coords._odata[sU] >= (tmin + tshift)) && + (coords._odata[sU] <= (tmax + tshift))); - timeSlices = Reduce(t_mask); + //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 )); - if (timeSlices > 0) - { - unsigned int sF = sU * LLs; - for (unsigned int s = 0; s < LLs; ++s) + Integer timeSlices = Reduce(t_mask); + + if (timeSlices > 0) { - 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; + unsigned int sF = sU * LLs + s; + Kernels::SeqConservedCurrentSiteBwd(tmp2._odata[sU], + q_out._odata[sF], Umu, sU, + mu, t_mask, axial_sign); } - } + } } } + + + + FermOpTemplateInstantiate(WilsonFermion5D); GparityFermOpTemplateInstantiate(WilsonFermion5D); diff --git a/lib/qcd/action/fermion/WilsonFermion5D.h b/lib/qcd/action/fermion/WilsonFermion5D.h index d66f4a1d..21da4c31 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.h +++ b/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 mom, - unsigned int tmin, - unsigned int tmax); + unsigned int tmin, + unsigned int tmax, + ComplexField &lattice_cmplx); }; }} diff --git a/tests/hadrons/Test_QED.cc b/tests/hadrons/Test_QED.cc new file mode 100644 index 00000000..44dd0578 --- /dev/null +++ b/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 + Author: Vera Guelpers + + 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 + +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 flavour = {"h"}; //{"l", "s", "c1", "c2", "c3"}; + std::vector 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("gauge"); + // pt source + MSource::Point::Par ptPar; + ptPar.position = "0 0 0 0"; + application.createModule("pt", ptPar); + // sink + MSink::Point::Par sinkPar; + sinkPar.mom = "0 0 0"; + application.createModule("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("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("DWF_" + flavour[i], actionPar); + + + // solvers + MSolver::RBPrecCG::Par solverPar; + solverPar.action = "DWF_" + flavour[i]; + solverPar.residual = 1.0e-8; + application.createModule("CG_" + flavour[i], + solverPar); + + // propagators + MFermion::GaugeProp::Par quarkPar; + quarkPar.solver = "CG_" + flavour[i]; + quarkPar.source = "pt"; + application.createModule("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("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("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("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("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("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("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 = ""; + mesPar.sink = "sink"; + application.createModule("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 = ""; + mesPar_seq_T.sink = "sink"; + application.createModule("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 = ""; + mesPar_seq_E.sink = "sink"; + application.createModule("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 = ""; + mesPar_seq_S.sink = "sink"; + application.createModule("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; +} diff --git a/tests/hadrons/Test_hadrons.hpp b/tests/hadrons/Test_hadrons.hpp index 0265f5a6..6188a580 100644 --- a/tests/hadrons/Test_hadrons.hpp +++ b/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(srcName, seqPar); }