diff --git a/Hadrons/Modules/MNPR/ConservedBilinear.cpp b/Hadrons/Modules/MNPR/ConservedBilinear.cpp new file mode 100644 index 00000000..f9e695b5 --- /dev/null +++ b/Hadrons/Modules/MNPR/ConservedBilinear.cpp @@ -0,0 +1,33 @@ +/* + * ConservedBilinear.cpp, part of Hadrons (https://github.com/aportelli/Hadrons) + * + * Copyright (C) 2015 - 2022 + * + * Author: Antonin Portelli + * + * Hadrons 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. + * + * Hadrons 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 Hadrons. If not, see . + * + * 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 MNPR; + +template class Grid::Hadrons::MNPR::TConservedBilinear; + diff --git a/Hadrons/Modules/MNPR/ConservedBilinear.hpp b/Hadrons/Modules/MNPR/ConservedBilinear.hpp new file mode 100644 index 00000000..094d1a26 --- /dev/null +++ b/Hadrons/Modules/MNPR/ConservedBilinear.hpp @@ -0,0 +1,220 @@ +/* + * ConservedBilinear.hpp, part of Hadrons (https://github.com/aportelli/Hadrons) + * + * Copyright (C) 2015 - 2022 + * + * Author: Antonin Portelli + * Author: Julia Kettle J.R.Kettle-2@sms.ed.ac.uk + * Author: Peter Boyle + * Author: Fabian Joswig + * Author: Felix Erben + * + * Hadrons 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. + * + * Hadrons 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 Hadrons. If not, see . + * + * See the full license in the file "LICENSE" in the top level distribution + * directory. + */ + +/* END LEGAL */ + +#ifndef Hadrons_ConservedBilinear_hpp_ +#define Hadrons_ConservedBilinear_hpp_ + +#include +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * TConservedBilinear * + Performs bilinear contractions of the type tr[g5*adj(qOut')*g5*G*qIn'] + Suitable for non exceptional momenta in Rome-Southampton NPR where G is + the non-local conserved vector current which depends on the action. + +Returns a spin-colour matrix, with indices si,sj, ci,cj +**************************************************************************/ +BEGIN_MODULE_NAMESPACE(MNPR) + +class ConservedBilinearPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(ConservedBilinearPar, + std::string, action, + std::string, qIn, + std::string, qOut, + std::string, pIn, + std::string, pOut, + std::string, output); +}; + +template +class TConservedBilinear: public Module +{ +public: + FERM_TYPE_ALIASES(FImpl,); + class Metadata: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Metadata, + Gamma::Algebra, gamma, + std::string, pIn, + std::string, pOut); + }; + typedef Correlator Result; +public: + // constructor + TConservedBilinear(const std::string name); + // destructor + virtual ~TConservedBilinear(void) {}; + // dependencies/products + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + virtual std::vector getOutputFiles(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +private: +}; + +MODULE_REGISTER_TMP(ConservedBilinear, ARG(TConservedBilinear), MNPR); + +/****************************************************************************** + * TConservedBilinear implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TConservedBilinear::TConservedBilinear(const std::string name) +: Module(name) +{} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TConservedBilinear::setup(void) +{ + LOG(Message) << "Running setup for ConservedBilinear" << std::endl; + + // The propagator can be 4d or 5d, but must match the action + const unsigned int ActionLs_{ env().getObjectLs(par().action) }; + if (ActionLs_ != 1) + { + std::string sError{ "ConservedBilinear currently only implemented for 4d actions."}; + HADRONS_ERROR(Size, sError); + } + + envTmpLat(PropagatorField, "qIn_phased"); + envTmpLat(PropagatorField, "qOut_phased"); + envTmpLat(PropagatorField, "lret"); + envTmpLat(PropagatorField, "dummy_phys_source"); + envTmpLat(ComplexField, "pDotXIn"); + envTmpLat(ComplexField, "pDotXOut"); + envTmpLat(ComplexField, "xMu"); +} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TConservedBilinear::getInput(void) +{ + std::vector input = {par().action, par().qIn, par().qOut}; + + return input; +} + +template +std::vector TConservedBilinear::getOutput(void) +{ + std::vector out = {}; + + return out; +} + +template +std::vector TConservedBilinear::getOutputFiles(void) +{ + std::vector output = {resultFilename(par().output)}; + + return output; +} + +template +void TConservedBilinear::execute(void) +{ + + LOG(Message) << "Computing conserved bilinear contractions '" << getName() << "' using" + << " propagators '" << par().qIn << "' and '" << par().qOut << "'" + << "with action " << par().action << std::endl; + + auto &act = envGet(FMat, par().action); + + // Propagators + auto &qIn = envGet(PropagatorField, par().qIn); + auto &qOut = envGet(PropagatorField, par().qOut); + envGetTmp(PropagatorField, qIn_phased); + envGetTmp(PropagatorField, qOut_phased); + envGetTmp(PropagatorField, lret); + envGetTmp(PropagatorField, dummy_phys_source); + envGetTmp(ComplexField, pDotXIn); + envGetTmp(ComplexField, pDotXOut); + envGetTmp(ComplexField, xMu); + + // momentum on legs + std::vector pIn = strToVec(par().pIn), + pOut = strToVec(par().pOut); + Coordinate latt_size = GridDefaultLatt(); + Complex Ci(0.0,1.0); + std::vector result; + Result r; + + Real volume = 1.0; + for (int mu = 0; mu < Nd; mu++) { + volume *= latt_size[mu]; + } + + NPRUtils::dot(pDotXIn,pIn); + qIn_phased = qIn * exp(-Ci * pDotXIn); + NPRUtils::dot(pDotXOut,pOut); + qOut_phased = qOut * exp(-Ci * pDotXOut); + + + r.info.pIn = par().pIn; + r.info.pOut = par().pOut; + + Gamma::Algebra Gmu [] = { + Gamma::Algebra::GammaX, + Gamma::Algebra::GammaY, + Gamma::Algebra::GammaZ, + Gamma::Algebra::GammaT, + }; + + for (int mu = 0; mu < Nd; mu++) + { + r.info.gamma = Gmu[mu]; + act.ContractConservedCurrent(qOut_phased, qIn_phased, lret, dummy_phys_source, Current::Vector, mu); + r.corr.push_back( (1.0 / volume) * sum_large(lret) ); + result.push_back(r); + r.corr.erase(r.corr.begin()); + } + + ////////////////////////////////////////////////// + saveResult(par().output, "ConservedBilinear", result); + LOG(Message) << "Complete. Writing results to " << par().output << std::endl; +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_ConservedBilinear_hpp_ \ No newline at end of file diff --git a/tests/Make.inc b/tests/Make.inc index e353baaf..a2041d57 100644 --- a/tests/Make.inc +++ b/tests/Make.inc @@ -1,8 +1,11 @@ -tests-local: Test_QED_Local Test_QED Test_diskvector Test_exact_distil Test_free_prop Test_hadrons_meson_3pt Test_hadrons_spectrum Test_sigma_to_nucleon Test_stoch_distil Test_xi_to_sigma -EXTRA_PROGRAMS = Test_QED_Local Test_QED Test_diskvector Test_exact_distil Test_free_prop Test_hadrons_meson_3pt Test_hadrons_spectrum Test_sigma_to_nucleon Test_stoch_distil Test_xi_to_sigma +tests-local: Test_conserved_bilinear Test_QED_Local Test_QED Test_diskvector Test_exact_distil Test_free_prop Test_hadrons_meson_3pt Test_hadrons_spectrum Test_sigma_to_nucleon Test_stoch_distil Test_xi_to_sigma +EXTRA_PROGRAMS = Test_conserved_bilinear Test_QED_Local Test_QED Test_diskvector Test_exact_distil Test_free_prop Test_hadrons_meson_3pt Test_hadrons_spectrum Test_sigma_to_nucleon Test_stoch_distil Test_xi_to_sigma -Test_QED_Local_SOURCES=Test_QED.cc -Test_QED_LDADD=-lHadrons -lGrid +Test_conserved_bilinear_SOURCES=Test_conserved_bilinear.cc +Test_conserved_bilinear_LDADD=-lHadrons -lGrid + +Test_QED_Local_SOURCES=Test_QED_Local.cc +Test_QED_Local_LDADD=-lHadrons -lGrid Test_QED_SOURCES=Test_QED.cc Test_QED_LDADD=-lHadrons -lGrid diff --git a/tests/Makefile.am b/tests/Makefile.am index b96e8d24..c9c47829 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -1,7 +1,8 @@ AM_LDFLAGS += -L$(top_builddir)/Hadrons AM_CXXFLAGS += -I$(top_srcdir) - + EXTRA_PROGRAMS = \ + Test_conserved_bilinear \ Test_QED_Local \ Test_QED \ Test_database \ @@ -21,6 +22,9 @@ CLEANFILES = $(EXTRA_PROGRAMS) tests-local: $(EXTRA_PROGRAMS) +Test_conserved_bilinear_SOURCES=Test_conserved_bilinear.cpp +Test_conserved_bilinear_LDADD=-lHadrons -lGrid + Test_QED_Local_SOURCES=Test_QED_Local.cpp Test_QED_Local_LDADD=-lHadrons -lGrid diff --git a/tests/Test_conserved_bilinear.cpp b/tests/Test_conserved_bilinear.cpp new file mode 100644 index 00000000..8d3eed71 --- /dev/null +++ b/tests/Test_conserved_bilinear.cpp @@ -0,0 +1,137 @@ +/* + * Test_hadrons_meson_3pt.cpp, part of Hadrons (https://github.com/aportelli/Hadrons) + * + * Copyright (C) 2015 - 2022 + * + * Author: Antonin Portelli + * Author: Fabian Joswig + * + * Hadrons 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. + * + * Hadrons 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 Hadrons. If not, see . + * + * See the full license in the file "LICENSE" in the top level distribution + * directory. + */ + +/* END LEGAL */ + +#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; + + // global parameters + Application::GlobalPar globalPar; + globalPar.trajCounter.start = 1500; + globalPar.trajCounter.end = 1520; + globalPar.trajCounter.step = 20; + globalPar.runId = "test_conserved_bilinear"; + globalPar.genetic.maxGen = 1000; + globalPar.genetic.maxCstGen = 200; + globalPar.genetic.popSize = 20; + globalPar.genetic.mutationRate = .1; + application.setPar(globalPar); + + // gauge field + application.createModule("gauge"); + + // set fermion boundary conditions to be periodic space, antiperiodic time. + std::string boundary = "1 1 1 1"; + std::string twist = "0. 0. 0. 0."; + + // actions + MAction::WilsonExpClover::Par actionPar; + actionPar.gauge = "gauge"; + actionPar.mass = 0.01; + actionPar.csw_r = 1.12; + actionPar.csw_t = 1.12; + actionPar.cF = 1.0; + actionPar.boundary = boundary; + actionPar.twist = twist; + application.createModule("eWC", actionPar); + + // solvers + MSolver::RBPrecCG::Par solverPar; + solverPar.action = "eWC"; + solverPar.residual = 1.0e-15; + solverPar.maxIteration = 10000; + application.createModule("CG", solverPar); + + // momentum source1 + MSource::Momentum::Par momentumPar1; + momentumPar1.mom = "2 2 0 0"; + application.createModule("mom1", momentumPar1); + + // momentum source2 + MSource::Momentum::Par momentumPar2; + momentumPar2.mom = "2 0 2 0"; + application.createModule("mom2", momentumPar2); + + // propagator1 + MFermion::GaugeProp::Par quarkPar1; + quarkPar1.solver = "CG"; + quarkPar1.source = "mom1"; + application.createModule("prop1", quarkPar1); + + // propagator2 + MFermion::GaugeProp::Par quarkPar2; + quarkPar2.solver = "CG"; + quarkPar2.source = "mom2"; + application.createModule("prop2", quarkPar2); + + // ExternalLeg1 + MNPR::ExternalLeg::Par externalLegPar1; + externalLegPar1.qIn = "prop1"; + externalLegPar1.pIn = "2 2 0 0"; + application.createModule("leg1", externalLegPar1); + + // ExternalLeg2 + MNPR::ExternalLeg::Par externalLegPar2; + externalLegPar2.qIn = "prop2"; + externalLegPar2.pIn = "2 0 2 0"; + application.createModule("leg2", externalLegPar2); + + // ConservedBilinear + MNPR::ConservedBilinear::Par ConservedBilinearPar; + ConservedBilinearPar.action = "eWC"; + ConservedBilinearPar.qIn = "prop1"; + ConservedBilinearPar.pIn = "2 2 0 0"; + ConservedBilinearPar.qOut = "prop2"; + ConservedBilinearPar.pOut = "2 0 2 0"; + application.createModule("ConservedBilinear", ConservedBilinearPar); + + // execution + application.saveParameterFile("conserved_bilinear.xml"); + application.run(); + + // epilogue + LOG(Message) << "Grid is finalizing now" << std::endl; + Grid_finalize(); + + return EXIT_SUCCESS; +}