diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index 4b420517..502c14af 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -22,12 +22,15 @@ #include #include #include +#include #include #include #include #include -#include #include +#include +#include +#include #include #include #include diff --git a/extras/Hadrons/Modules/MGauge/StochEm.cc b/extras/Hadrons/Modules/MGauge/StochEm.cc index 7b07557e..1727706c 100644 --- a/extras/Hadrons/Modules/MGauge/StochEm.cc +++ b/extras/Hadrons/Modules/MGauge/StochEm.cc @@ -7,6 +7,7 @@ Source file: extras/Hadrons/Modules/MGauge/StochEm.cc Copyright (C) 2015-2018 Author: Antonin Portelli +Author: James Harrison Author: Vera Guelpers This program is free software; you can redistribute it and/or modify @@ -58,12 +59,8 @@ 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; - } + weightDone_ = env().hasCreatedObject("_" + getName() + "_weight"); + envCacheLat(EmComp, "_" + getName() + "_weight"); envCreateLat(EmField, getName()); } @@ -72,13 +69,14 @@ void TStochEm::execute(void) { LOG(Message) << "Generating stochastic EM potential..." << std::endl; - PhotonR photon(par().gauge, par().zmScheme); + std::vector improvements = strToVec(par().improvement); + PhotonR photon(par().gauge, par().zmScheme, improvements, par().G0_qedInf); auto &a = envGet(EmField, getName()); auto &w = envGet(EmComp, "_" + getName() + "_weight"); - if (create_weight) + if (!weightDone_) { - LOG(Message) << "Caching stochatic EM potential weight (gauge: " + LOG(Message) << "Caching stochastic EM potential weight (gauge: " << par().gauge << ", zero-mode scheme: " << par().zmScheme << ")..." << std::endl; photon.StochasticWeight(w); diff --git a/extras/Hadrons/Modules/MGauge/StochEm.hpp b/extras/Hadrons/Modules/MGauge/StochEm.hpp index a933afab..5eb9f1e9 100644 --- a/extras/Hadrons/Modules/MGauge/StochEm.hpp +++ b/extras/Hadrons/Modules/MGauge/StochEm.hpp @@ -7,6 +7,7 @@ Source file: extras/Hadrons/Modules/MGauge/StochEm.hpp Copyright (C) 2015-2018 Author: Antonin Portelli +Author: James Harrison Author: Vera Guelpers This program is free software; you can redistribute it and/or modify @@ -45,7 +46,9 @@ class StochEmPar: Serializable public: GRID_SERIALIZABLE_CLASS_MEMBERS(StochEmPar, PhotonR::Gauge, gauge, - PhotonR::ZmScheme, zmScheme); + PhotonR::ZmScheme, zmScheme, + std::string, improvement, + Real, G0_qedInf); }; class TStochEm: public Module @@ -61,13 +64,13 @@ public: // dependency relation virtual std::vector getInput(void); virtual std::vector getOutput(void); -private: - bool create_weight; protected: // setup virtual void setup(void); // execution virtual void execute(void); +private: + bool weightDone_; }; MODULE_REGISTER(StochEm, TStochEm, MGauge); diff --git a/extras/Hadrons/Modules/MGauge/UnitEm.cc b/extras/Hadrons/Modules/MGauge/UnitEm.cc new file mode 100644 index 00000000..6260f10a --- /dev/null +++ b/extras/Hadrons/Modules/MGauge/UnitEm.cc @@ -0,0 +1,69 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MGauge/StochEm.cc + +Copyright (C) 2015 +Copyright (C) 2016 + +Author: James Harrison + +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 + +using namespace Grid; +using namespace Hadrons; +using namespace MGauge; + +/****************************************************************************** +* TStochEm implementation * +******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +TUnitEm::TUnitEm(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +std::vector TUnitEm::getInput(void) +{ + return std::vector(); +} + +std::vector TUnitEm::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +void TUnitEm::setup(void) +{ + envCreateLat(EmField, getName()); +} + +// execution /////////////////////////////////////////////////////////////////// +void TUnitEm::execute(void) +{ + PhotonR photon(0, 0); // Just chose arbitrary input values here + auto &a = envGet(EmField, getName()); + LOG(Message) << "Generating unit EM potential..." << std::endl; + photon.UnitField(a); +} diff --git a/extras/Hadrons/Modules/MGauge/UnitEm.hpp b/extras/Hadrons/Modules/MGauge/UnitEm.hpp new file mode 100644 index 00000000..07ee5bba --- /dev/null +++ b/extras/Hadrons/Modules/MGauge/UnitEm.hpp @@ -0,0 +1,69 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MGauge/StochEm.hpp + +Copyright (C) 2015 +Copyright (C) 2016 + +Author: James Harrison + +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_MGauge_UnitEm_hpp_ +#define Hadrons_MGauge_UnitEm_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * StochEm * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MGauge) + +class TUnitEm: public Module +{ +public: + typedef PhotonR::GaugeField EmField; + typedef PhotonR::GaugeLinkField EmComp; +public: + // constructor + TUnitEm(const std::string name); + // destructor + virtual ~TUnitEm(void) {}; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); +protected: + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER(UnitEm, TUnitEm, MGauge); + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MGauge_UnitEm_hpp_ diff --git a/extras/Hadrons/Modules/MScalar/ChargedProp.cc b/extras/Hadrons/Modules/MScalar/ChargedProp.cc index 7d59ac3e..88c3dffd 100644 --- a/extras/Hadrons/Modules/MScalar/ChargedProp.cc +++ b/extras/Hadrons/Modules/MScalar/ChargedProp.cc @@ -51,7 +51,8 @@ std::vector TChargedProp::getInput(void) std::vector TChargedProp::getOutput(void) { - std::vector out = {getName()}; + std::vector out = {getName(), getName()+"_0", getName()+"_Q", + getName()+"_Sun", getName()+"_Tad"}; return out; } @@ -66,18 +67,27 @@ void TChargedProp::setup(void) phaseName_.push_back("_shiftphase_" + std::to_string(mu)); } GFSrcName_ = getName() + "_DinvSrc"; + prop0Name_ = getName() + "_0"; + propQName_ = getName() + "_Q"; + propSunName_ = getName() + "_Sun"; + propTadName_ = getName() + "_Tad"; fftName_ = getName() + "_fft"; freeMomPropDone_ = env().hasCreatedObject(freeMomPropName_); GFSrcDone_ = env().hasCreatedObject(GFSrcName_); phasesDone_ = env().hasCreatedObject(phaseName_[0]); + prop0Done_ = env().hasCreatedObject(prop0Name_); envCacheLat(ScalarField, freeMomPropName_); for (unsigned int mu = 0; mu < env().getNd(); ++mu) { envCacheLat(ScalarField, phaseName_[mu]); } envCacheLat(ScalarField, GFSrcName_); + envCacheLat(ScalarField, prop0Name_); envCreateLat(ScalarField, getName()); + envCreateLat(ScalarField, propQName_); + envCreateLat(ScalarField, propSunName_); + envCreateLat(ScalarField, propTadName_); envTmpLat(ScalarField, "buf"); envTmpLat(ScalarField, "result"); envTmpLat(ScalarField, "Amu"); @@ -95,79 +105,125 @@ void TChargedProp::execute(void) << " (mass= " << par().mass << ", charge= " << par().charge << ")..." << std::endl; - auto &prop = envGet(ScalarField, getName()); - auto &GFSrc = envGet(ScalarField, GFSrcName_); - auto &G = envGet(ScalarField, freeMomPropName_); - auto &fft = envGet(FFT, fftName_); - double q = par().charge; - envGetTmp(ScalarField, result); - envGetTmp(ScalarField, buf); + auto &prop = envGet(ScalarField, getName()); + 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_); + double q = par().charge; + envGetTmp(ScalarField, buf); - // G*F*Src - prop = GFSrc; + // -G*momD1*G*F*Src (momD1 = F*D1*Finv) + propQ = GFSrc; + momD1(propQ, fft); + propQ = -G*propQ; + propSun = -propQ; + fft.FFT_dim(propQ, propQ, env().getNd()-1, FFT::backward); - // - q*G*momD1*G*F*Src (momD1 = F*D1*Finv) - buf = GFSrc; - momD1(buf, fft); - buf = G*buf; - prop = prop - q*buf; + // G*momD1*G*momD1*G*F*Src (here buf = G*momD1*G*F*Src) + momD1(propSun, fft); + propSun = G*propSun; + fft.FFT_dim(propSun, propSun, env().getNd()-1, FFT::backward); - // + q^2*G*momD1*G*momD1*G*F*Src (here buf = G*momD1*G*F*Src) - momD1(buf, fft); - prop = prop + q*q*G*buf; - - // - q^2*G*momD2*G*F*Src (momD2 = F*D2*Finv) - buf = GFSrc; - momD2(buf, fft); - prop = prop - q*q*G*buf; - - // final FT - fft.FFT_all_dim(prop, prop, FFT::backward); + // -G*momD2*G*F*Src (momD2 = F*D2*Finv) + propTad = GFSrc; + momD2(propTad, fft); + propTad = -G*propTad; + fft.FFT_dim(propTad, propTad, env().getNd()-1, FFT::backward); + // full charged scalar propagator + fft.FFT_dim(buf, GFSrc, env().getNd()-1, FFT::backward); + prop = buf + q*propQ + q*q*propSun + q*q*propTad; + // OUTPUT IF NECESSARY if (!par().output.empty()) { - std::string filename = par().output + "." + - std::to_string(vm().getTrajectory()); - - LOG(Message) << "Saving zero-momentum projection to '" - << filename << "'..." << std::endl; - - std::vector vecBuf; - std::vector result; - - sliceSum(prop, vecBuf, Tp); - result.resize(vecBuf.size()); - for (unsigned int t = 0; t < vecBuf.size(); ++t) + Result result; + TComplex site; + std::vector siteCoor; + + LOG(Message) << "Saving momentum-projected propagator to '" + << RESULT_FILE_NAME(par().output) << "'..." + << std::endl; + result.projection.resize(par().outputMom.size()); + result.lattice_size = env().getGrid()->_fdimensions; + result.mass = par().mass; + result.charge = q; + siteCoor.resize(env().getNd()); + for (unsigned int i_p = 0; i_p < par().outputMom.size(); ++i_p) { - result[t] = TensorRemove(vecBuf[t]); + result.projection[i_p].momentum = strToVec(par().outputMom[i_p]); + + LOG(Message) << "Calculating (" << par().outputMom[i_p] + << ") momentum projection" << std::endl; + + result.projection[i_p].corr_0.resize(env().getGrid()->_fdimensions[env().getNd()-1]); + result.projection[i_p].corr.resize(env().getGrid()->_fdimensions[env().getNd()-1]); + result.projection[i_p].corr_Q.resize(env().getGrid()->_fdimensions[env().getNd()-1]); + result.projection[i_p].corr_Sun.resize(env().getGrid()->_fdimensions[env().getNd()-1]); + result.projection[i_p].corr_Tad.resize(env().getGrid()->_fdimensions[env().getNd()-1]); + + for (unsigned int j = 0; j < env().getNd()-1; ++j) + { + siteCoor[j] = result.projection[i_p].momentum[j]; + } + + for (unsigned int t = 0; t < result.projection[i_p].corr.size(); ++t) + { + siteCoor[env().getNd()-1] = t; + peekSite(site, prop, siteCoor); + result.projection[i_p].corr[t]=TensorRemove(site); + peekSite(site, buf, siteCoor); + result.projection[i_p].corr_0[t]=TensorRemove(site); + peekSite(site, propQ, siteCoor); + result.projection[i_p].corr_Q[t]=TensorRemove(site); + peekSite(site, propSun, siteCoor); + result.projection[i_p].corr_Sun[t]=TensorRemove(site); + peekSite(site, propTad, siteCoor); + result.projection[i_p].corr_Tad[t]=TensorRemove(site); + } } - saveResult(par().output, "charge", q); saveResult(par().output, "prop", result); } + + std::vector mask(env().getNd(),1); + mask[env().getNd()-1] = 0; + fft.FFT_dim_mask(prop, prop, mask, FFT::backward); + fft.FFT_dim_mask(propQ, propQ, mask, FFT::backward); + fft.FFT_dim_mask(propSun, propSun, mask, FFT::backward); + fft.FFT_dim_mask(propTad, propTad, mask, FFT::backward); } void TChargedProp::makeCaches(void) { auto &freeMomProp = envGet(ScalarField, freeMomPropName_); auto &GFSrc = envGet(ScalarField, GFSrcName_); + auto &prop0 = envGet(ScalarField, prop0Name_); auto &fft = envGet(FFT, fftName_); if (!freeMomPropDone_) { - LOG(Message) << "Caching momentum space free scalar propagator" + LOG(Message) << "Caching momentum-space free scalar propagator" << " (mass= " << par().mass << ")..." << std::endl; SIMPL::MomentumSpacePropagator(freeMomProp, par().mass); } if (!GFSrcDone_) { - FFT fft(env().getGrid()); auto &source = envGet(ScalarField, par().source); - + LOG(Message) << "Caching G*F*src..." << std::endl; fft.FFT_all_dim(GFSrc, source, FFT::forward); GFSrc = freeMomProp*GFSrc; } + if (!prop0Done_) + { + LOG(Message) << "Caching position-space free scalar propagator..." + << std::endl; + fft.FFT_all_dim(prop0, GFSrc, FFT::backward); + } if (!phasesDone_) { std::vector &l = env().getGrid()->_fdimensions; @@ -184,6 +240,14 @@ void TChargedProp::makeCaches(void) phase_.push_back(&phmu); } } + else + { + phase_.clear(); + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + { + phase_.push_back(env().getObject(phaseName_[mu])); + } + } } void TChargedProp::momD1(ScalarField &s, FFT &fft) diff --git a/extras/Hadrons/Modules/MScalar/ChargedProp.hpp b/extras/Hadrons/Modules/MScalar/ChargedProp.hpp index 549d6154..3ca7940f 100644 --- a/extras/Hadrons/Modules/MScalar/ChargedProp.hpp +++ b/extras/Hadrons/Modules/MScalar/ChargedProp.hpp @@ -7,6 +7,7 @@ Source file: extras/Hadrons/Modules/MScalar/ChargedProp.hpp Copyright (C) 2015-2018 Author: Antonin Portelli +Author: James Harrison 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 @@ -47,7 +48,8 @@ public: std::string, source, double, mass, double, charge, - std::string, output); + std::string, output, + std::vector, outputMom); }; class TChargedProp: public Module @@ -56,6 +58,26 @@ public: SCALAR_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, momentum, + std::vector, corr, + std::vector, corr_0, + std::vector, corr_Q, + std::vector, corr_Sun, + std::vector, corr_Tad); + }; + GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + std::vector, lattice_size, + double, mass, + double, charge, + std::vector, projection); + }; public: // constructor TChargedProp(const std::string name); @@ -74,8 +96,10 @@ private: void momD1(ScalarField &s, FFT &fft); void momD2(ScalarField &s, FFT &fft); private: - bool freeMomPropDone_, GFSrcDone_, phasesDone_; - std::string freeMomPropName_, GFSrcName_, fftName_; + bool freeMomPropDone_, GFSrcDone_, prop0Done_, + phasesDone_; + std::string freeMomPropName_, GFSrcName_, prop0Name_, + propQName_, propSunName_, propTadName_, fftName_; std::vector phaseName_; std::vector phase_; }; diff --git a/extras/Hadrons/Modules/MScalar/ScalarVP.cc b/extras/Hadrons/Modules/MScalar/ScalarVP.cc new file mode 100644 index 00000000..fdb9586e --- /dev/null +++ b/extras/Hadrons/Modules/MScalar/ScalarVP.cc @@ -0,0 +1,536 @@ +#include +#include +#include + +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(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +std::vector TScalarVP::getInput(void) +{ + prop0Name_ = par().scalarProp + "_0"; + propQName_ = par().scalarProp + "_Q"; + propSunName_ = par().scalarProp + "_Sun"; + propTadName_ = par().scalarProp + "_Tad"; + + std::vector in = {par().emField, prop0Name_, propQName_, + propSunName_, propTadName_}; + + return in; +} + +std::vector TScalarVP::getOutput(void) +{ + std::vector 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(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 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(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 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 coor0 = {0, 0, 0, 0}; + std::vector > vpTensor; + for (unsigned int mu = 0; mu < env().getNd(); ++mu) + { + std::vector 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(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(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: = ) + 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) << "'..." + << 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 &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 mom = strToVec(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 &projection, const ScalarField &vp, int i_p) +{ + std::vector 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; +} diff --git a/extras/Hadrons/Modules/MScalar/ScalarVP.hpp b/extras/Hadrons/Modules/MScalar/ScalarVP.hpp new file mode 100644 index 00000000..741325d1 --- /dev/null +++ b/extras/Hadrons/Modules/MScalar/ScalarVP.hpp @@ -0,0 +1,101 @@ +#ifndef Hadrons_MScalar_ScalarVP_hpp_ +#define Hadrons_MScalar_ScalarVP_hpp_ + +#include +#include +#include + +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, outputMom); +}; + +class TScalarVP: public Module +{ +public: + SCALAR_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, momentum, + std::vector>>, pi, + std::vector>>, pi_free, + std::vector>>, pi_2E, + std::vector>>, pi_2T, + std::vector>>, pi_S, + std::vector>>, pi_4C, + std::vector>>, pi_X, + std::vector>>, pi_srcT, + std::vector>>, pi_snkT); + }; + GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + std::vector, lattice_size, + double, mass, + double, charge, + std::vector, projection); + }; +public: + // constructor + TScalarVP(const std::string name); + // destructor + virtual ~TScalarVP(void) {}; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector 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 &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 phaseName_, muPropQName_, + momPhaseName_; + std::vector > vpTensorName_; + std::vector phase_, momPhase_; +}; + +MODULE_REGISTER(ScalarVP, TScalarVP, MScalar); + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MScalar_ScalarVP_hpp_ diff --git a/extras/Hadrons/Modules/MScalar/VPCounterTerms.cc b/extras/Hadrons/Modules/MScalar/VPCounterTerms.cc new file mode 100644 index 00000000..b3393679 --- /dev/null +++ b/extras/Hadrons/Modules/MScalar/VPCounterTerms.cc @@ -0,0 +1,232 @@ +#include +#include + +using namespace Grid; +using namespace Hadrons; +using namespace MScalar; + +/****************************************************************************** +* TVPCounterTerms implementation * +******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +TVPCounterTerms::TVPCounterTerms(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +std::vector TVPCounterTerms::getInput(void) +{ + std::vector in = {par().source}; + + return in; +} + +std::vector TVPCounterTerms::getOutput(void) +{ + std::vector 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 &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(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) << "'..." + << std::endl; + saveResult(par().output, "scalar_loops", outputData); + } +} + +void TVPCounterTerms::project(std::vector &projection, const ScalarField &vp, int i_p) +{ + std::vector 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]); + } +} diff --git a/extras/Hadrons/Modules/MScalar/VPCounterTerms.hpp b/extras/Hadrons/Modules/MScalar/VPCounterTerms.hpp new file mode 100644 index 00000000..99e1b1ce --- /dev/null +++ b/extras/Hadrons/Modules/MScalar/VPCounterTerms.hpp @@ -0,0 +1,75 @@ +#ifndef Hadrons_MScalar_VPCounterTerms_hpp_ +#define Hadrons_MScalar_VPCounterTerms_hpp_ + +#include +#include +#include + +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, outputMom); +}; + +class TVPCounterTerms: public Module +{ +public: + SCALAR_TYPE_ALIASES(SIMPL,); + class Result: Serializable + { + public: + class Projection: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Projection, + std::vector, momentum, + std::vector>>, twoScalar, + std::vector>>, threeScalar, + std::vector>>, pSquaredInsertion); + }; + GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + std::vector, lattice_size, + double, mass, + std::vector, projection); + }; +public: + // constructor + TVPCounterTerms(const std::string name); + // destructor + virtual ~TVPCounterTerms(void) {}; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); +protected: + // setup + virtual void setup(void); + // execution + virtual void execute(void); +private: + void project(std::vector &projection, const ScalarField &vp, int i_p); +private: + std::string freeMomPropName_, GFSrcName_, phatsqName_, prop0Name_, + twoscalarName_, twoscalarVertexName_, + psquaredName_, psquaredVertexName_; + std::vector phaseName_, momPhaseName_; + std::vector phase_, momPhase_; +}; + +MODULE_REGISTER(VPCounterTerms, TVPCounterTerms, MScalar); + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MScalar_VPCounterTerms_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index 57dc8d20..0e44176b 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -19,6 +19,7 @@ modules_cc =\ Modules/MSolver/RBPrecCG.cc \ Modules/MSolver/LocalCoherenceLanczos.cc \ Modules/MGauge/Unit.cc \ + Modules/MGauge/UnitEm.cc \ Modules/MGauge/StochEm.cc \ Modules/MGauge/Random.cc \ Modules/MGauge/FundtoHirep.cc \ @@ -27,6 +28,8 @@ modules_cc =\ Modules/MLoop/NoiseLoop.cc \ Modules/MScalar/FreeProp.cc \ Modules/MScalar/ChargedProp.cc \ + Modules/MScalar/ScalarVP.cc \ + Modules/MScalar/VPCounterTerms.cc \ Modules/MAction/Wilson.cc \ Modules/MAction/ZMobiusDWF.cc \ Modules/MAction/WilsonClover.cc \ @@ -69,16 +72,19 @@ modules_hpp =\ Modules/MSink/Point.hpp \ Modules/MSolver/LocalCoherenceLanczos.hpp \ Modules/MSolver/RBPrecCG.hpp \ - Modules/MGauge/Unit.hpp \ Modules/MGauge/Random.hpp \ Modules/MGauge/FundtoHirep.hpp \ Modules/MGauge/StochEm.hpp \ + Modules/MGauge/Unit.hpp \ + Modules/MGauge/UnitEm.hpp \ Modules/MUtilities/TestSeqGamma.hpp \ Modules/MUtilities/TestSeqConserved.hpp \ Modules/MLoop/NoiseLoop.hpp \ + Modules/MScalar/ChargedProp.hpp \ Modules/MScalar/FreeProp.hpp \ Modules/MScalar/Scalar.hpp \ - Modules/MScalar/ChargedProp.hpp \ + Modules/MScalar/ScalarVP.hpp \ + Modules/MScalar/VPCounterTerms.hpp \ Modules/MAction/DWF.hpp \ Modules/MAction/Wilson.hpp \ Modules/MAction/WilsonClover.hpp \ diff --git a/lib/lattice/Lattice_rng.h b/lib/lattice/Lattice_rng.h index edf9dd23..654c5dd5 100644 --- a/lib/lattice/Lattice_rng.h +++ b/lib/lattice/Lattice_rng.h @@ -158,10 +158,19 @@ namespace Grid { // tens of seconds per trajectory so this is clean in all reasonable cases, // and margin of safety is orders of magnitude. // We could hack Sitmo to skip in the higher order words of state if necessary + // + // Replace with 2^30 ; avoid problem on large volumes + // ///////////////////////////////////////////////////////////////////////////////////// // uint64_t skip = site+1; // Old init Skipped then drew. Checked compat with faster init + const int shift = 30; + uint64_t skip = site; - skip = skip<<40; + + skip = skip<> shift)==site); // check for overflow + eng.discard(skip); // std::cout << " Engine " < improvements); + Photon(Gauge gauge, ZmScheme zmScheme, Real G0); + Photon(Gauge gauge, ZmScheme zmScheme, std::vector improvements, Real G0); virtual ~Photon(void) = default; void FreePropagator(const GaugeField &in, GaugeField &out); void MomentumSpacePropagator(const GaugeField &in, GaugeField &out); @@ -68,21 +71,44 @@ namespace QCD{ void StochasticField(GaugeField &out, GridParallelRNG &rng); void StochasticField(GaugeField &out, GridParallelRNG &rng, const GaugeLinkField &weight); + void UnitField(GaugeField &out); private: + void infVolPropagator(GaugeLinkField &out); void invKHatSquared(GaugeLinkField &out); void zmSub(GaugeLinkField &out); private: Gauge gauge_; ZmScheme zmScheme_; + std::vector improvement_; + Real G0_; }; typedef Photon PhotonR; template Photon::Photon(Gauge gauge, ZmScheme zmScheme) - : gauge_(gauge), zmScheme_(zmScheme) + : gauge_(gauge), zmScheme_(zmScheme), improvement_(std::vector()), + G0_(0.15493339023106021408483720810737508876916113364521) {} - + + template + Photon::Photon(Gauge gauge, ZmScheme zmScheme, + std::vector improvements) + : gauge_(gauge), zmScheme_(zmScheme), improvement_(improvements), + G0_(0.15493339023106021408483720810737508876916113364521) + {} + + template + Photon::Photon(Gauge gauge, ZmScheme zmScheme, Real G0) + : gauge_(gauge), zmScheme_(zmScheme), improvement_(std::vector()), G0_(G0) + {} + + template + Photon::Photon(Gauge gauge, ZmScheme zmScheme, + std::vector improvements, Real G0) + : gauge_(gauge), zmScheme_(zmScheme), improvement_(improvements), G0_(G0) + {} + template void Photon::FreePropagator (const GaugeField &in,GaugeField &out) { @@ -95,6 +121,34 @@ namespace QCD{ MomentumSpacePropagator(prop_k,in_k); theFFT.FFT_all_dim(out,prop_k,FFT::backward); } + + template + void Photon::infVolPropagator(GaugeLinkField &out) + { + auto *grid = dynamic_cast(out._grid); + LatticeReal xmu(grid); + GaugeLinkField one(grid); + const unsigned int nd = grid->_ndimension; + std::vector &l = grid->_fdimensions; + std::vector x0(nd,0); + TComplex Tone = Complex(1.0,0.0); + TComplex Tzero = Complex(G0_,0.0); + FFT fft(grid); + + one = Complex(1.0,0.0); + out = zero; + for(int mu = 0; mu < nd; mu++) + { + LatticeCoordinate(xmu,mu); + Real lo2 = l[mu]/2.0; + xmu = where(xmu < lo2, xmu, xmu-double(l[mu])); + out = out + toComplex(4*M_PI*M_PI*xmu*xmu); + } + pokeSite(Tone, out, x0); + out = one/out; + pokeSite(Tzero, out, x0); + fft.FFT_all_dim(out, out, FFT::forward); + } template void Photon::invKHatSquared(GaugeLinkField &out) @@ -127,6 +181,7 @@ namespace QCD{ { GridBase *grid = out._grid; const unsigned int nd = grid->_ndimension; + std::vector &l = grid->_fdimensions; switch (zmScheme_) { @@ -148,28 +203,49 @@ namespace QCD{ for(int d = 0; d < grid->_ndimension - 1; d++) { LatticeCoordinate(coor,d); + coor = where(coor < Integer(l[d]/2), coor, coor-Integer(l[d])); spNrm = spNrm + coor*coor; } out = where(spNrm == Integer(0), 0.*out, out); - - break; + + // IR improvement + for(int i = 0; i < improvement_.size(); i++) + { + Real f = sqrt(improvement_[i]+1); + out = where(spNrm == Integer(i+1), f*out, out); + } } default: break; } } - + template void Photon::MomentumSpacePropagator(const GaugeField &in, GaugeField &out) { - GridBase *grid = out._grid; - LatticeComplex k2Inv(grid); + GridBase *grid = out._grid; + LatticeComplex momProp(grid); - invKHatSquared(k2Inv); - zmSub(k2Inv); + switch (zmScheme_) + { + case ZmScheme::qedTL: + case ZmScheme::qedL: + { + invKHatSquared(momProp); + zmSub(momProp); + break; + } + case ZmScheme::qedInf: + { + infVolPropagator(momProp); + break; + } + default: + break; + } - out = in*k2Inv; + out = in*momProp; } template @@ -179,14 +255,30 @@ namespace QCD{ const unsigned int nd = grid->_ndimension; std::vector latt_size = grid->_fdimensions; - Integer vol = 1; - for(int d = 0; d < nd; d++) + switch (zmScheme_) { - vol = vol * latt_size[d]; + case ZmScheme::qedTL: + case ZmScheme::qedL: + { + Integer vol = 1; + for(int d = 0; d < nd; d++) + { + vol = vol * latt_size[d]; + } + invKHatSquared(weight); + weight = sqrt(vol)*sqrt(weight); + zmSub(weight); + break; + } + case ZmScheme::qedInf: + { + infVolPropagator(weight); + weight = sqrt(real(weight)); + break; + } + default: + break; } - invKHatSquared(weight); - weight = sqrt(vol)*sqrt(weight); - zmSub(weight); } template @@ -209,13 +301,52 @@ namespace QCD{ GaugeField aTilde(grid); FFT fft(grid); + switch (zmScheme_) + { + case ZmScheme::qedTL: + case ZmScheme::qedL: + { + for(int mu = 0; mu < nd; mu++) + { + gaussian(rng, r); + r = weight*r; + pokeLorentz(aTilde, r, mu); + } + break; + } + case ZmScheme::qedInf: + { + Complex shift(1., 1.); // This needs to be a GaugeLink element? + for(int mu = 0; mu < nd; mu++) + { + bernoulli(rng, r); + r = weight*(2.*r - shift); + pokeLorentz(aTilde, r, mu); + } + break; + } + default: + break; + } + + fft.FFT_all_dim(out, aTilde, FFT::backward); + + out = real(out); + } + + template + void Photon::UnitField(GaugeField &out) + { + auto *grid = dynamic_cast(out._grid); + const unsigned int nd = grid->_ndimension; + GaugeLinkField r(grid); + + r = Complex(1.0,0.0); + for(int mu = 0; mu < nd; mu++) { - gaussian(rng, r); - r = weight*r; - pokeLorentz(aTilde, r, mu); + pokeLorentz(out, r, mu); } - fft.FFT_all_dim(out, aTilde, FFT::backward); out = real(out); }