From 52a856b4a8a9d1fc4b76f8d10e2e638aac267fed Mon Sep 17 00:00:00 2001 From: Vera Guelpers Date: Thu, 3 May 2018 12:33:20 +0100 Subject: [PATCH] FreeProp module for Hadrons --- extras/Hadrons/Modules.hpp | 1 + extras/Hadrons/Modules/MFermion/FreeProp.cc | 36 +++ extras/Hadrons/Modules/MFermion/FreeProp.hpp | 189 ++++++++++++++ extras/Hadrons/modules.inc | 2 + tests/hadrons/Test_free_prop.cc | 245 +++++++++++++++++++ 5 files changed, 473 insertions(+) create mode 100644 extras/Hadrons/Modules/MFermion/FreeProp.cc create mode 100644 extras/Hadrons/Modules/MFermion/FreeProp.hpp create mode 100644 tests/hadrons/Test_free_prop.cc diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index fc536393..2c356cfc 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -37,6 +37,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include #include #include #include diff --git a/extras/Hadrons/Modules/MFermion/FreeProp.cc b/extras/Hadrons/Modules/MFermion/FreeProp.cc new file mode 100644 index 00000000..b8eb0529 --- /dev/null +++ b/extras/Hadrons/Modules/MFermion/FreeProp.cc @@ -0,0 +1,36 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MFermion/FreeProp.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 +*************************************************************************************/ +/* END LEGAL */ +#include + +using namespace Grid; +using namespace Hadrons; +using namespace MFermion; + +template class Grid::Hadrons::MFermion::TFreeProp; + diff --git a/extras/Hadrons/Modules/MFermion/FreeProp.hpp b/extras/Hadrons/Modules/MFermion/FreeProp.hpp new file mode 100644 index 00000000..b1038ddd --- /dev/null +++ b/extras/Hadrons/Modules/MFermion/FreeProp.hpp @@ -0,0 +1,189 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MFermion/FreeProp.hpp + +Copyright (C) 2015-2018 + +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 +*************************************************************************************/ +/* END LEGAL */ + + +#ifndef Hadrons_MFermion_FreeProp_hpp_ +#define Hadrons_MFermion_FreeProp_hpp_ + +#include +#include +#include + +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * FreeProp * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MFermion) + +class FreePropPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(FreePropPar, + std::string, source, + std::string, action, + double, mass, + std::string, twist); +}; + +template +class TFreeProp: public Module +{ +public: + FGS_TYPE_ALIASES(FImpl,); +public: + // constructor + TFreeProp(const std::string name); + // destructor + virtual ~TFreeProp(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: + unsigned int Ls_; +}; + +MODULE_REGISTER_TMP(FreeProp, TFreeProp, MFermion); + +/****************************************************************************** + * TFreeProp implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TFreeProp::TFreeProp(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TFreeProp::getInput(void) +{ + std::vector in = {par().source, par().action}; + + return in; +} + +template +std::vector TFreeProp::getOutput(void) +{ + std::vector out = {getName(), getName() + "_5d"}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TFreeProp::setup(void) +{ + Ls_ = env().getObjectLs(par().action); + envCreateLat(PropagatorField, getName()); + envTmpLat(FermionField, "source", Ls_); + envTmpLat(FermionField, "sol", Ls_); + envTmpLat(FermionField, "tmp"); + if (Ls_ > 1) + { + envCreateLat(PropagatorField, getName() + "_5d", Ls_); + } +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TFreeProp::execute(void) +{ + LOG(Message) << "Computing free fermion propagator '" << getName() << "'" + << std::endl; + + std::string propName = (Ls_ == 1) ? getName() : (getName() + "_5d"); + auto &prop = envGet(PropagatorField, propName); + auto &fullSrc = envGet(PropagatorField, par().source); + auto &mat = envGet(FMat, par().action); + RealD mass = par().mass; + + envGetTmp(FermionField, source); + envGetTmp(FermionField, sol); + envGetTmp(FermionField, tmp); + LOG(Message) << "Calculating a free Propagator with mass " << mass + << " using the action '" << par().action + << "' on source '" << par().source << "'" << std::endl; + for (unsigned int s = 0; s < Ns; ++s) + for (unsigned int c = 0; c < FImpl::Dimension; ++c) + { + LOG(Message) << "Calculation for spin= " << s << ", color= " << c + << std::endl; + // source conversion for 4D sources + if (!env().isObject5d(par().source)) + { + if (Ls_ == 1) + { + PropToFerm(source, fullSrc, s, c); + } + else + { + PropToFerm(tmp, fullSrc, s, c); + make_5D(tmp, source, Ls_); + } + } + // source conversion for 5D sources + else + { + if (Ls_ != env().getObjectLs(par().source)) + { + HADRONS_ERROR(Size, "Ls mismatch between quark action and source"); + } + else + { + PropToFerm(source, fullSrc, s, c); + } + } + sol = zero; + std::vector twist = strToVec(par().twist); + if(twist.size() != Nd) HADRONS_ERROR(Size, "number of twist angles does not match number of dimensions"); + mat.FreePropagator(source,sol,mass,twist); + FermToProp(prop, sol, s, c); + // create 4D propagators from 5D one if necessary + if (Ls_ > 1) + { + PropagatorField &p4d = envGet(PropagatorField, getName()); + make_4D(sol, tmp, Ls_); + FermToProp(p4d, tmp, s, c); + } + } +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MFermion_FreeProp_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index ad3a8727..69463f9d 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -7,6 +7,7 @@ modules_cc =\ Modules/MContraction/WardIdentity.cc \ Modules/MContraction/DiscLoop.cc \ Modules/MContraction/Gamma3pt.cc \ + Modules/MFermion/FreeProp.cc \ Modules/MFermion/GaugeProp.cc \ Modules/MSource/Point.cc \ Modules/MSource/Wall.cc \ @@ -54,6 +55,7 @@ modules_hpp =\ Modules/MContraction/Gamma3pt.hpp \ Modules/MContraction/WardIdentity.hpp \ Modules/MContraction/WeakHamiltonianEye.hpp \ + Modules/MFermion/FreeProp.hpp \ Modules/MFermion/GaugeProp.hpp \ Modules/MSource/SeqGamma.hpp \ Modules/MSource/Point.hpp \ diff --git a/tests/hadrons/Test_free_prop.cc b/tests/hadrons/Test_free_prop.cc new file mode 100644 index 00000000..a1a5aadd --- /dev/null +++ b/tests/hadrons/Test_free_prop.cc @@ -0,0 +1,245 @@ +/******************************************************************************* + Grid physics library, www.github.com/paboyle/Grid + + Source file: tests/hadrons/Test_free_prop.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 +#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 }; + std::vector lepton_flavour = {"mu"}; + std::vector lepton_mass = {.2}; + + 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"); + // unit gauge field for lepton + application.createModule("free_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"; + + + //Propagators from FFT and Feynman rules + for (unsigned int i = 0; i < lepton_mass.size(); ++i) + { + //DWF actions + MAction::DWF::Par actionPar_lep; + actionPar_lep.gauge = "free_gauge"; + actionPar_lep.Ls = 8; + actionPar_lep.M5 = 1.8; + actionPar_lep.mass = lepton_mass[i]; + actionPar_lep.boundary = boundary; + application.createModule("free_DWF_" + lepton_flavour[i], actionPar_lep); + + //DWF free propagators + MFermion::FreeProp::Par freePar; + freePar.source = "pt"; + freePar.action = "free_DWF_" + lepton_flavour[i]; + freePar.twist = "0 0 0 0.5"; + freePar.mass = lepton_mass[i]; + application.createModule("Lpt_" + lepton_flavour[i], + freePar); + + //Wilson actions + MAction::Wilson::Par actionPar_lep_W; + actionPar_lep_W.gauge = "free_gauge"; + actionPar_lep_W.mass = lepton_mass[i]; + actionPar_lep_W.boundary = boundary; + application.createModule("free_W_" + lepton_flavour[i], actionPar_lep_W); + + //Wilson free propagators + MFermion::FreeProp::Par freePar_W; + freePar_W.source = "pt"; + freePar_W.action = "free_W_" + lepton_flavour[i]; + freePar_W.twist = "0 0 0 0.5"; + freePar_W.mass = lepton_mass[i]; + application.createModule("W_Lpt_" + lepton_flavour[i], + freePar_W); + + + } + + + + //Propagators from inversion + for (unsigned int i = 0; i < flavour.size(); ++i) + { + //DWF 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; + solverPar.maxIteration = 10000; + application.createModule("CG_" + flavour[i], + solverPar); + + //DWF propagators + MFermion::GaugeProp::Par quarkPar; + quarkPar.solver = "CG_" + flavour[i]; + quarkPar.source = "pt"; + application.createModule("Qpt_" + flavour[i], + quarkPar); + + + + //Wilson actions + MAction::Wilson::Par actionPar_W; + actionPar_W.gauge = "gauge"; + actionPar_W.mass = mass[i]; + actionPar_W.boundary = boundary; + application.createModule("W_" + flavour[i], actionPar_W); + + + // solvers + MSolver::RBPrecCG::Par solverPar_W; + solverPar_W.action = "W_" + flavour[i]; + solverPar_W.residual = 1.0e-8; + solverPar_W.maxIteration = 10000; + application.createModule("W_CG_" + flavour[i], + solverPar_W); + + //Wilson propagators + MFermion::GaugeProp::Par quarkPar_W; + quarkPar_W.solver = "W_CG_" + flavour[i]; + quarkPar_W.source = "pt"; + application.createModule("W_Qpt_" + flavour[i], + quarkPar_W); + + } + + + //2pt contraction for Propagators from FFT and Feynman rules + for (unsigned int i = 0; i < lepton_flavour.size(); ++i) + for (unsigned int j = i; j < lepton_flavour.size(); ++j) + { + //2pt function contraction DWF + MContraction::Meson::Par freemesPar; + freemesPar.output = "2pt_free/DWF_L_pt_" + lepton_flavour[i] + lepton_flavour[j]; + freemesPar.q1 = "Lpt_" + lepton_flavour[i]; + freemesPar.q2 = "Lpt_" + lepton_flavour[j]; + freemesPar.gammas = "(Gamma5 Gamma5)"; + freemesPar.sink = "sink"; + application.createModule("meson_L_pt_" + + lepton_flavour[i] + lepton_flavour[j], + freemesPar); + + //2pt function contraction Wilson + MContraction::Meson::Par freemesPar_W; + freemesPar_W.output = "2pt_free/W_L_pt_" + lepton_flavour[i] + lepton_flavour[j]; + freemesPar_W.q1 = "W_Lpt_" + lepton_flavour[i]; + freemesPar_W.q2 = "W_Lpt_" + lepton_flavour[j]; + freemesPar_W.gammas = "(Gamma5 Gamma5)"; + freemesPar_W.sink = "sink"; + application.createModule("W_meson_L_pt_" + + lepton_flavour[i] + lepton_flavour[j], + freemesPar_W); + + } + + //2pt contraction for Propagators from inverion + for (unsigned int i = 0; i < flavour.size(); ++i) + for (unsigned int j = i; j < flavour.size(); ++j) + { + //2pt function contraction DWF + MContraction::Meson::Par mesPar; + mesPar.output = "2pt_free/DWF_pt_" + flavour[i] + flavour[j]; + mesPar.q1 = "Qpt_" + flavour[i]; + mesPar.q2 = "Qpt_" + flavour[j]; + mesPar.gammas = "(Gamma5 Gamma5)"; + mesPar.sink = "sink"; + application.createModule("meson_pt_" + + flavour[i] + flavour[j], + mesPar); + + + //2pt function contraction Wilson + MContraction::Meson::Par mesPar_W; + mesPar_W.output = "2pt_free/W_pt_" + flavour[i] + flavour[j]; + mesPar_W.q1 = "W_Qpt_" + flavour[i]; + mesPar_W.q2 = "W_Qpt_" + flavour[j]; + mesPar_W.gammas = "(Gamma5 Gamma5)"; + mesPar_W.sink = "sink"; + application.createModule("W_meson_pt_" + + flavour[i] + flavour[j], + mesPar_W); + + } + + + + // execution + application.saveParameterFile("free_prop.xml"); + application.run(); + + // epilogue + LOG(Message) << "Grid is finalizing now" << std::endl; + Grid_finalize(); + + return EXIT_SUCCESS; +}