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