mirror of
https://github.com/paboyle/Grid.git
synced 2025-04-04 11:15:55 +01:00
Merge branch 'feature/hadrons' into develop
This commit is contained in:
commit
c9e9e8061d
@ -162,7 +162,8 @@ void Application::saveParameterFile(const std::string parameterFileName)
|
||||
sizeString((size)*locVol_) << " (" << sizeString(size) << "/site)"
|
||||
|
||||
#define DEFINE_MEMPEAK \
|
||||
auto memPeak = [this](const std::vector<unsigned int> &program)\
|
||||
GeneticScheduler<unsigned int>::ObjFunc memPeak = \
|
||||
[this](const std::vector<unsigned int> &program)\
|
||||
{\
|
||||
unsigned int memPeak;\
|
||||
bool msg;\
|
||||
|
@ -145,6 +145,15 @@ std::string typeName(void)
|
||||
return typeName(typeIdPt<T>());
|
||||
}
|
||||
|
||||
// default writers/readers
|
||||
#ifdef HAVE_HDF5
|
||||
typedef Hdf5Reader CorrReader;
|
||||
typedef Hdf5Writer CorrWriter;
|
||||
#else
|
||||
typedef XmlReader CorrReader;
|
||||
typedef XmlWriter CorrWriter;
|
||||
#endif
|
||||
|
||||
END_HADRONS_NAMESPACE
|
||||
|
||||
#endif // Hadrons_Global_hpp_
|
||||
|
@ -29,12 +29,20 @@ See the full license in the file "LICENSE" in the top level distribution directo
|
||||
#include <Grid/Hadrons/Modules/MAction/DWF.hpp>
|
||||
#include <Grid/Hadrons/Modules/MAction/Wilson.hpp>
|
||||
#include <Grid/Hadrons/Modules/MContraction/Baryon.hpp>
|
||||
#include <Grid/Hadrons/Modules/MContraction/DiscLoop.hpp>
|
||||
#include <Grid/Hadrons/Modules/MContraction/Gamma3pt.hpp>
|
||||
#include <Grid/Hadrons/Modules/MContraction/Meson.hpp>
|
||||
#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
|
||||
#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp>
|
||||
#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp>
|
||||
#include <Grid/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp>
|
||||
#include <Grid/Hadrons/Modules/MGauge/Load.hpp>
|
||||
#include <Grid/Hadrons/Modules/MGauge/Random.hpp>
|
||||
#include <Grid/Hadrons/Modules/MGauge/Unit.hpp>
|
||||
#include <Grid/Hadrons/Modules/MLoop/NoiseLoop.hpp>
|
||||
#include <Grid/Hadrons/Modules/MSolver/RBPrecCG.hpp>
|
||||
#include <Grid/Hadrons/Modules/MSource/Point.hpp>
|
||||
#include <Grid/Hadrons/Modules/MSource/SeqGamma.hpp>
|
||||
#include <Grid/Hadrons/Modules/MSource/Wall.hpp>
|
||||
#include <Grid/Hadrons/Modules/MSource/Z2.hpp>
|
||||
#include <Grid/Hadrons/Modules/Quark.hpp>
|
||||
|
@ -112,7 +112,7 @@ void TBaryon<FImpl1, FImpl2, FImpl3>::execute(void)
|
||||
<< " quarks '" << par().q1 << "', '" << par().q2 << "', and '"
|
||||
<< par().q3 << "'" << std::endl;
|
||||
|
||||
XmlWriter writer(par().output);
|
||||
CorrWriter writer(par().output);
|
||||
PropagatorField1 &q1 = *env().template getObject<PropagatorField1>(par().q1);
|
||||
PropagatorField2 &q2 = *env().template getObject<PropagatorField2>(par().q2);
|
||||
PropagatorField3 &q3 = *env().template getObject<PropagatorField3>(par().q2);
|
||||
|
144
extras/Hadrons/Modules/MContraction/DiscLoop.hpp
Normal file
144
extras/Hadrons/Modules/MContraction/DiscLoop.hpp
Normal file
@ -0,0 +1,144 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: extras/Hadrons/Modules/MContraction/DiscLoop.hpp
|
||||
|
||||
Copyright (C) 2017
|
||||
|
||||
Author: Andrew Lawson <andrew.lawson1991@gmail.com>
|
||||
|
||||
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_DiscLoop_hpp_
|
||||
#define Hadrons_DiscLoop_hpp_
|
||||
|
||||
#include <Grid/Hadrons/Global.hpp>
|
||||
#include <Grid/Hadrons/Module.hpp>
|
||||
#include <Grid/Hadrons/ModuleFactory.hpp>
|
||||
|
||||
BEGIN_HADRONS_NAMESPACE
|
||||
|
||||
/******************************************************************************
|
||||
* DiscLoop *
|
||||
******************************************************************************/
|
||||
BEGIN_MODULE_NAMESPACE(MContraction)
|
||||
|
||||
class DiscLoopPar: Serializable
|
||||
{
|
||||
public:
|
||||
GRID_SERIALIZABLE_CLASS_MEMBERS(DiscLoopPar,
|
||||
std::string, q_loop,
|
||||
Gamma::Algebra, gamma,
|
||||
std::string, output);
|
||||
};
|
||||
|
||||
template <typename FImpl>
|
||||
class TDiscLoop: public Module<DiscLoopPar>
|
||||
{
|
||||
TYPE_ALIASES(FImpl,);
|
||||
class Result: Serializable
|
||||
{
|
||||
public:
|
||||
GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
|
||||
Gamma::Algebra, gamma,
|
||||
std::vector<Complex>, corr);
|
||||
};
|
||||
public:
|
||||
// constructor
|
||||
TDiscLoop(const std::string name);
|
||||
// destructor
|
||||
virtual ~TDiscLoop(void) = default;
|
||||
// dependency relation
|
||||
virtual std::vector<std::string> getInput(void);
|
||||
virtual std::vector<std::string> getOutput(void);
|
||||
// setup
|
||||
virtual void setup(void);
|
||||
// execution
|
||||
virtual void execute(void);
|
||||
};
|
||||
|
||||
MODULE_REGISTER_NS(DiscLoop, TDiscLoop<FIMPL>, MContraction);
|
||||
|
||||
/******************************************************************************
|
||||
* TDiscLoop implementation *
|
||||
******************************************************************************/
|
||||
// constructor /////////////////////////////////////////////////////////////////
|
||||
template <typename FImpl>
|
||||
TDiscLoop<FImpl>::TDiscLoop(const std::string name)
|
||||
: Module<DiscLoopPar>(name)
|
||||
{}
|
||||
|
||||
// dependencies/products ///////////////////////////////////////////////////////
|
||||
template <typename FImpl>
|
||||
std::vector<std::string> TDiscLoop<FImpl>::getInput(void)
|
||||
{
|
||||
std::vector<std::string> in = {par().q_loop};
|
||||
|
||||
return in;
|
||||
}
|
||||
|
||||
template <typename FImpl>
|
||||
std::vector<std::string> TDiscLoop<FImpl>::getOutput(void)
|
||||
{
|
||||
std::vector<std::string> out = {getName()};
|
||||
|
||||
return out;
|
||||
}
|
||||
|
||||
// setup ///////////////////////////////////////////////////////////////////////
|
||||
template <typename FImpl>
|
||||
void TDiscLoop<FImpl>::setup(void)
|
||||
{
|
||||
|
||||
}
|
||||
|
||||
// execution ///////////////////////////////////////////////////////////////////
|
||||
template <typename FImpl>
|
||||
void TDiscLoop<FImpl>::execute(void)
|
||||
{
|
||||
LOG(Message) << "Computing disconnected loop contraction '" << getName()
|
||||
<< "' using '" << par().q_loop << "' with " << par().gamma
|
||||
<< " insertion." << std::endl;
|
||||
|
||||
CorrWriter writer(par().output);
|
||||
PropagatorField &q_loop = *env().template getObject<PropagatorField>(par().q_loop);
|
||||
LatticeComplex c(env().getGrid());
|
||||
Gamma gamma(par().gamma);
|
||||
std::vector<TComplex> buf;
|
||||
Result result;
|
||||
|
||||
c = trace(gamma*q_loop);
|
||||
sliceSum(c, buf, Tp);
|
||||
|
||||
result.gamma = par().gamma;
|
||||
result.corr.resize(buf.size());
|
||||
for (unsigned int t = 0; t < buf.size(); ++t)
|
||||
{
|
||||
result.corr[t] = TensorRemove(buf[t]);
|
||||
}
|
||||
|
||||
write(writer, "disc", result);
|
||||
}
|
||||
|
||||
END_MODULE_NAMESPACE
|
||||
|
||||
END_HADRONS_NAMESPACE
|
||||
|
||||
#endif // Hadrons_DiscLoop_hpp_
|
170
extras/Hadrons/Modules/MContraction/Gamma3pt.hpp
Normal file
170
extras/Hadrons/Modules/MContraction/Gamma3pt.hpp
Normal file
@ -0,0 +1,170 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: extras/Hadrons/Modules/MContraction/Gamma3pt.hpp
|
||||
|
||||
Copyright (C) 2017
|
||||
|
||||
Author: Andrew Lawson <andrew.lawson1991@gmail.com>
|
||||
|
||||
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_Gamma3pt_hpp_
|
||||
#define Hadrons_Gamma3pt_hpp_
|
||||
|
||||
#include <Grid/Hadrons/Global.hpp>
|
||||
#include <Grid/Hadrons/Module.hpp>
|
||||
#include <Grid/Hadrons/ModuleFactory.hpp>
|
||||
|
||||
BEGIN_HADRONS_NAMESPACE
|
||||
|
||||
/*
|
||||
* 3pt contraction with gamma matrix insertion.
|
||||
*
|
||||
* Schematic:
|
||||
*
|
||||
* q2 q3
|
||||
* /----<------*------<----¬
|
||||
* / gamma \
|
||||
* / \
|
||||
* i * * f
|
||||
* \ /
|
||||
* \ /
|
||||
* \----------->----------/
|
||||
* q1
|
||||
*
|
||||
* trace(g5*q1*adj(q2)*g5*gamma*q3)
|
||||
*/
|
||||
|
||||
/******************************************************************************
|
||||
* Gamma3pt *
|
||||
******************************************************************************/
|
||||
BEGIN_MODULE_NAMESPACE(MContraction)
|
||||
|
||||
class Gamma3ptPar: Serializable
|
||||
{
|
||||
public:
|
||||
GRID_SERIALIZABLE_CLASS_MEMBERS(Gamma3ptPar,
|
||||
std::string, q1,
|
||||
std::string, q2,
|
||||
std::string, q3,
|
||||
Gamma::Algebra, gamma,
|
||||
std::string, output);
|
||||
};
|
||||
|
||||
template <typename FImpl1, typename FImpl2, typename FImpl3>
|
||||
class TGamma3pt: public Module<Gamma3ptPar>
|
||||
{
|
||||
TYPE_ALIASES(FImpl1, 1);
|
||||
TYPE_ALIASES(FImpl2, 2);
|
||||
TYPE_ALIASES(FImpl3, 3);
|
||||
class Result: Serializable
|
||||
{
|
||||
public:
|
||||
GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
|
||||
Gamma::Algebra, gamma,
|
||||
std::vector<Complex>, corr);
|
||||
};
|
||||
public:
|
||||
// constructor
|
||||
TGamma3pt(const std::string name);
|
||||
// destructor
|
||||
virtual ~TGamma3pt(void) = default;
|
||||
// dependency relation
|
||||
virtual std::vector<std::string> getInput(void);
|
||||
virtual std::vector<std::string> getOutput(void);
|
||||
// setup
|
||||
virtual void setup(void);
|
||||
// execution
|
||||
virtual void execute(void);
|
||||
};
|
||||
|
||||
MODULE_REGISTER_NS(Gamma3pt, ARG(TGamma3pt<FIMPL, FIMPL, FIMPL>), MContraction);
|
||||
|
||||
/******************************************************************************
|
||||
* TGamma3pt implementation *
|
||||
******************************************************************************/
|
||||
// constructor /////////////////////////////////////////////////////////////////
|
||||
template <typename FImpl1, typename FImpl2, typename FImpl3>
|
||||
TGamma3pt<FImpl1, FImpl2, FImpl3>::TGamma3pt(const std::string name)
|
||||
: Module<Gamma3ptPar>(name)
|
||||
{}
|
||||
|
||||
// dependencies/products ///////////////////////////////////////////////////////
|
||||
template <typename FImpl1, typename FImpl2, typename FImpl3>
|
||||
std::vector<std::string> TGamma3pt<FImpl1, FImpl2, FImpl3>::getInput(void)
|
||||
{
|
||||
std::vector<std::string> in = {par().q1, par().q2, par().q3};
|
||||
|
||||
return in;
|
||||
}
|
||||
|
||||
template <typename FImpl1, typename FImpl2, typename FImpl3>
|
||||
std::vector<std::string> TGamma3pt<FImpl1, FImpl2, FImpl3>::getOutput(void)
|
||||
{
|
||||
std::vector<std::string> out = {getName()};
|
||||
|
||||
return out;
|
||||
}
|
||||
|
||||
// setup ///////////////////////////////////////////////////////////////////////
|
||||
template <typename FImpl1, typename FImpl2, typename FImpl3>
|
||||
void TGamma3pt<FImpl1, FImpl2, FImpl3>::setup(void)
|
||||
{
|
||||
|
||||
}
|
||||
|
||||
// execution ///////////////////////////////////////////////////////////////////
|
||||
template <typename FImpl1, typename FImpl2, typename FImpl3>
|
||||
void TGamma3pt<FImpl1, FImpl2, FImpl3>::execute(void)
|
||||
{
|
||||
LOG(Message) << "Computing 3pt contractions '" << getName() << "' using"
|
||||
<< " quarks '" << par().q1 << "', '" << par().q2 << "' and '"
|
||||
<< par().q3 << "', with " << par().gamma << " insertion."
|
||||
<< std::endl;
|
||||
|
||||
CorrWriter writer(par().output);
|
||||
PropagatorField1 &q1 = *env().template getObject<PropagatorField1>(par().q1);
|
||||
PropagatorField2 &q2 = *env().template getObject<PropagatorField2>(par().q2);
|
||||
PropagatorField3 &q3 = *env().template getObject<PropagatorField3>(par().q3);
|
||||
LatticeComplex c(env().getGrid());
|
||||
Gamma g5(Gamma::Algebra::Gamma5);
|
||||
Gamma gamma(par().gamma);
|
||||
std::vector<TComplex> buf;
|
||||
Result result;
|
||||
|
||||
c = trace(g5*q1*adj(q2)*(g5*gamma)*q3);
|
||||
sliceSum(c, buf, Tp);
|
||||
|
||||
result.gamma = par().gamma;
|
||||
result.corr.resize(buf.size());
|
||||
for (unsigned int t = 0; t < buf.size(); ++t)
|
||||
{
|
||||
result.corr[t] = TensorRemove(buf[t]);
|
||||
}
|
||||
|
||||
write(writer, "gamma3pt", result);
|
||||
}
|
||||
|
||||
END_MODULE_NAMESPACE
|
||||
|
||||
END_HADRONS_NAMESPACE
|
||||
|
||||
#endif // Hadrons_Gamma3pt_hpp_
|
@ -6,8 +6,10 @@ Source file: extras/Hadrons/Modules/MContraction/Meson.hpp
|
||||
|
||||
Copyright (C) 2015
|
||||
Copyright (C) 2016
|
||||
Copyright (C) 2017
|
||||
|
||||
Author: Antonin Portelli <antonin.portelli@me.com>
|
||||
Andrew Lawson <andrew.lawson1991@gmail.com>
|
||||
|
||||
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
|
||||
@ -36,20 +38,39 @@ See the full license in the file "LICENSE" in the top level distribution directo
|
||||
|
||||
BEGIN_HADRONS_NAMESPACE
|
||||
|
||||
/*
|
||||
|
||||
Meson contractions
|
||||
-----------------------------
|
||||
|
||||
* options:
|
||||
- q1: input propagator 1 (string)
|
||||
- q2: input propagator 2 (string)
|
||||
- gammas: gamma products to insert at sink & source, pairs of gamma matrices
|
||||
(space-separated strings) in angled brackets (i.e. <g_sink g_src>),
|
||||
in a sequence (e.g. "<Gamma5 Gamma5><Gamma5 GammaT>").
|
||||
|
||||
Special values: "all" - perform all possible contractions.
|
||||
- mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0."),
|
||||
given as multiples of (2*pi) / L.
|
||||
*/
|
||||
|
||||
/******************************************************************************
|
||||
* TMeson *
|
||||
******************************************************************************/
|
||||
BEGIN_MODULE_NAMESPACE(MContraction)
|
||||
|
||||
typedef std::pair<Gamma::Algebra, Gamma::Algebra> GammaPair;
|
||||
|
||||
class MesonPar: Serializable
|
||||
{
|
||||
public:
|
||||
GRID_SERIALIZABLE_CLASS_MEMBERS(MesonPar,
|
||||
std::string, q1,
|
||||
std::string, q2,
|
||||
std::string, output,
|
||||
Gamma::Algebra, gammaSource,
|
||||
Gamma::Algebra, gammaSink);
|
||||
std::string, q1,
|
||||
std::string, q2,
|
||||
std::string, gammas,
|
||||
std::string, mom,
|
||||
std::string, output);
|
||||
};
|
||||
|
||||
template <typename FImpl1, typename FImpl2>
|
||||
@ -61,7 +82,10 @@ public:
|
||||
class Result: Serializable
|
||||
{
|
||||
public:
|
||||
GRID_SERIALIZABLE_CLASS_MEMBERS(Result, std::vector<Complex>, corr);
|
||||
GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
|
||||
Gamma::Algebra, gamma_snk,
|
||||
Gamma::Algebra, gamma_src,
|
||||
std::vector<Complex>, corr);
|
||||
};
|
||||
public:
|
||||
// constructor
|
||||
@ -71,6 +95,7 @@ public:
|
||||
// dependencies/products
|
||||
virtual std::vector<std::string> getInput(void);
|
||||
virtual std::vector<std::string> getOutput(void);
|
||||
virtual void parseGammaString(std::vector<GammaPair> &gammaList);
|
||||
// execution
|
||||
virtual void execute(void);
|
||||
};
|
||||
@ -103,6 +128,32 @@ std::vector<std::string> TMeson<FImpl1, FImpl2>::getOutput(void)
|
||||
return output;
|
||||
}
|
||||
|
||||
template <typename FImpl1, typename FImpl2>
|
||||
void TMeson<FImpl1, FImpl2>::parseGammaString(std::vector<GammaPair> &gammaList)
|
||||
{
|
||||
// Determine gamma matrices to insert at source/sink.
|
||||
if (par().gammas.compare("all") == 0)
|
||||
{
|
||||
// Do all contractions.
|
||||
unsigned int n_gam = Ns * Ns;
|
||||
gammaList.resize(n_gam*n_gam);
|
||||
for (unsigned int i = 1; i < Gamma::nGamma; i += 2)
|
||||
{
|
||||
for (unsigned int j = 1; j < Gamma::nGamma; j += 2)
|
||||
{
|
||||
gammaList.push_back(std::make_pair((Gamma::Algebra)i,
|
||||
(Gamma::Algebra)j));
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
// Parse individual contractions from input string.
|
||||
gammaList = strToVec<GammaPair>(par().gammas);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// execution ///////////////////////////////////////////////////////////////////
|
||||
template <typename FImpl1, typename FImpl2>
|
||||
void TMeson<FImpl1, FImpl2>::execute(void)
|
||||
@ -111,21 +162,44 @@ void TMeson<FImpl1, FImpl2>::execute(void)
|
||||
<< " quarks '" << par().q1 << "' and '" << par().q2 << "'"
|
||||
<< std::endl;
|
||||
|
||||
XmlWriter writer(par().output);
|
||||
PropagatorField1 &q1 = *env().template getObject<PropagatorField1>(par().q1);
|
||||
PropagatorField2 &q2 = *env().template getObject<PropagatorField2>(par().q2);
|
||||
LatticeComplex c(env().getGrid());
|
||||
Gamma gSrc(par().gammaSource), gSnk(par().gammaSink);
|
||||
Gamma g5(Gamma::Algebra::Gamma5);
|
||||
std::vector<TComplex> buf;
|
||||
Result result;
|
||||
|
||||
c = trace(gSnk*q1*adj(gSrc)*g5*adj(q2)*g5);
|
||||
sliceSum(c, buf, Tp);
|
||||
result.corr.resize(buf.size());
|
||||
for (unsigned int t = 0; t < buf.size(); ++t)
|
||||
CorrWriter writer(par().output);
|
||||
PropagatorField1 &q1 = *env().template getObject<PropagatorField1>(par().q1);
|
||||
PropagatorField2 &q2 = *env().template getObject<PropagatorField2>(par().q2);
|
||||
LatticeComplex c(env().getGrid());
|
||||
Gamma g5(Gamma::Algebra::Gamma5);
|
||||
std::vector<GammaPair> gammaList;
|
||||
std::vector<TComplex> buf;
|
||||
std::vector<Result> result;
|
||||
std::vector<Real> p;
|
||||
|
||||
p = strToVec<Real>(par().mom);
|
||||
LatticeComplex ph(env().getGrid()), coor(env().getGrid());
|
||||
Complex i(0.0,1.0);
|
||||
ph = zero;
|
||||
for(unsigned int mu = 0; mu < env().getNd(); mu++)
|
||||
{
|
||||
result.corr[t] = TensorRemove(buf[t]);
|
||||
LatticeCoordinate(coor, mu);
|
||||
ph = ph + p[mu]*coor*((1./(env().getGrid()->_fdimensions[mu])));
|
||||
}
|
||||
ph = exp((Real)(2*M_PI)*i*ph);
|
||||
|
||||
parseGammaString(gammaList);
|
||||
|
||||
result.resize(gammaList.size());
|
||||
for (unsigned int i = 0; i < result.size(); ++i)
|
||||
{
|
||||
Gamma gSnk(gammaList[i].first);
|
||||
Gamma gSrc(gammaList[i].second);
|
||||
c = trace((g5*gSnk)*q1*(adj(gSrc)*g5)*adj(q2))*ph;
|
||||
sliceSum(c, buf, Tp);
|
||||
|
||||
result[i].gamma_snk = gammaList[i].first;
|
||||
result[i].gamma_src = gammaList[i].second;
|
||||
result[i].corr.resize(buf.size());
|
||||
for (unsigned int t = 0; t < buf.size(); ++t)
|
||||
{
|
||||
result[i].corr[t] = TensorRemove(buf[t]);
|
||||
}
|
||||
}
|
||||
write(writer, "meson", result);
|
||||
}
|
||||
|
114
extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp
Normal file
114
extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp
Normal file
@ -0,0 +1,114 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp
|
||||
|
||||
Copyright (C) 2017
|
||||
|
||||
Author: Andrew Lawson <andrew.lawson1991@gmail.com>
|
||||
|
||||
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_WeakHamiltonian_hpp_
|
||||
#define Hadrons_WeakHamiltonian_hpp_
|
||||
|
||||
#include <Grid/Hadrons/Global.hpp>
|
||||
#include <Grid/Hadrons/Module.hpp>
|
||||
#include <Grid/Hadrons/ModuleFactory.hpp>
|
||||
|
||||
BEGIN_HADRONS_NAMESPACE
|
||||
|
||||
/******************************************************************************
|
||||
* WeakHamiltonian *
|
||||
******************************************************************************/
|
||||
BEGIN_MODULE_NAMESPACE(MContraction)
|
||||
|
||||
/*******************************************************************************
|
||||
* Utilities for contractions involving the Weak Hamiltonian.
|
||||
******************************************************************************/
|
||||
//// Sum and store correlator.
|
||||
#define MAKE_DIAG(exp, buf, res, n)\
|
||||
sliceSum(exp, buf, Tp);\
|
||||
res.name = (n);\
|
||||
res.corr.resize(buf.size());\
|
||||
for (unsigned int t = 0; t < buf.size(); ++t)\
|
||||
{\
|
||||
res.corr[t] = TensorRemove(buf[t]);\
|
||||
}
|
||||
|
||||
//// Contraction of mu index: use 'mu' variable in exp.
|
||||
#define SUM_MU(buf,exp)\
|
||||
buf = zero;\
|
||||
for (unsigned int mu = 0; mu < ndim; ++mu)\
|
||||
{\
|
||||
buf += exp;\
|
||||
}
|
||||
|
||||
enum
|
||||
{
|
||||
i_V = 0,
|
||||
i_A = 1,
|
||||
n_i = 2
|
||||
};
|
||||
|
||||
class WeakHamiltonianPar: Serializable
|
||||
{
|
||||
public:
|
||||
GRID_SERIALIZABLE_CLASS_MEMBERS(WeakHamiltonianPar,
|
||||
std::string, q1,
|
||||
std::string, q2,
|
||||
std::string, q3,
|
||||
std::string, q4,
|
||||
std::string, output);
|
||||
};
|
||||
|
||||
#define MAKE_WEAK_MODULE(modname)\
|
||||
class T##modname: public Module<WeakHamiltonianPar>\
|
||||
{\
|
||||
public:\
|
||||
TYPE_ALIASES(FIMPL,)\
|
||||
class Result: Serializable\
|
||||
{\
|
||||
public:\
|
||||
GRID_SERIALIZABLE_CLASS_MEMBERS(Result,\
|
||||
std::string, name,\
|
||||
std::vector<Complex>, corr);\
|
||||
};\
|
||||
public:\
|
||||
/* constructor */ \
|
||||
T##modname(const std::string name);\
|
||||
/* destructor */ \
|
||||
virtual ~T##modname(void) = default;\
|
||||
/* dependency relation */ \
|
||||
virtual std::vector<std::string> getInput(void);\
|
||||
virtual std::vector<std::string> getOutput(void);\
|
||||
/* setup */ \
|
||||
virtual void setup(void);\
|
||||
/* execution */ \
|
||||
virtual void execute(void);\
|
||||
std::vector<std::string> VA_label = {"V", "A"};\
|
||||
};\
|
||||
MODULE_REGISTER_NS(modname, T##modname, MContraction);
|
||||
|
||||
END_MODULE_NAMESPACE
|
||||
|
||||
END_HADRONS_NAMESPACE
|
||||
|
||||
#endif // Hadrons_WeakHamiltonian_hpp_
|
137
extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc
Normal file
137
extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc
Normal file
@ -0,0 +1,137 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc
|
||||
|
||||
Copyright (C) 2017
|
||||
|
||||
Author: Andrew Lawson <andrew.lawson1991@gmail.com>
|
||||
|
||||
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 <Grid/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp>
|
||||
|
||||
using namespace Grid;
|
||||
using namespace Hadrons;
|
||||
using namespace MContraction;
|
||||
|
||||
/*
|
||||
* Weak Hamiltonian current-current contractions, Eye-type.
|
||||
*
|
||||
* These contractions are generated by the Q1 and Q2 operators in the physical
|
||||
* basis (see e.g. Fig 3 of arXiv:1507.03094).
|
||||
*
|
||||
* Schematics: q4 |
|
||||
* /-<-¬ |
|
||||
* / \ | q2 q3
|
||||
* \ / | /----<------*------<----¬
|
||||
* q2 \ / q3 | / /-*-¬ \
|
||||
* /-----<-----* *-----<----¬ | / / \ \
|
||||
* i * H_W * f | i * \ / q4 * f
|
||||
* \ / | \ \->-/ /
|
||||
* \ / | \ /
|
||||
* \---------->---------/ | \----------->----------/
|
||||
* q1 | q1
|
||||
* |
|
||||
* Saucer (S) | Eye (E)
|
||||
*
|
||||
* S: trace(q3*g5*q1*adj(q2)*g5*gL[mu][p_1]*q4*gL[mu][p_2])
|
||||
* E: trace(q3*g5*q1*adj(q2)*g5*gL[mu][p_1])*trace(q4*gL[mu][p_2])
|
||||
*/
|
||||
|
||||
/******************************************************************************
|
||||
* TWeakHamiltonianEye implementation *
|
||||
******************************************************************************/
|
||||
// constructor /////////////////////////////////////////////////////////////////
|
||||
TWeakHamiltonianEye::TWeakHamiltonianEye(const std::string name)
|
||||
: Module<WeakHamiltonianPar>(name)
|
||||
{}
|
||||
|
||||
// dependencies/products ///////////////////////////////////////////////////////
|
||||
std::vector<std::string> TWeakHamiltonianEye::getInput(void)
|
||||
{
|
||||
std::vector<std::string> in = {par().q1, par().q2, par().q3, par().q4};
|
||||
|
||||
return in;
|
||||
}
|
||||
|
||||
std::vector<std::string> TWeakHamiltonianEye::getOutput(void)
|
||||
{
|
||||
std::vector<std::string> out = {getName()};
|
||||
|
||||
return out;
|
||||
}
|
||||
|
||||
// setup ///////////////////////////////////////////////////////////////////////
|
||||
void TWeakHamiltonianEye::setup(void)
|
||||
{
|
||||
|
||||
}
|
||||
|
||||
// execution ///////////////////////////////////////////////////////////////////
|
||||
void TWeakHamiltonianEye::execute(void)
|
||||
{
|
||||
LOG(Message) << "Computing Weak Hamiltonian (Eye type) contractions '"
|
||||
<< getName() << "' using quarks '" << par().q1 << "', '"
|
||||
<< par().q2 << ", '" << par().q3 << "' and '" << par().q4
|
||||
<< "'." << std::endl;
|
||||
|
||||
CorrWriter writer(par().output);
|
||||
PropagatorField &q1 = *env().template getObject<PropagatorField>(par().q1);
|
||||
PropagatorField &q2 = *env().template getObject<PropagatorField>(par().q2);
|
||||
PropagatorField &q3 = *env().template getObject<PropagatorField>(par().q3);
|
||||
PropagatorField &q4 = *env().template getObject<PropagatorField>(par().q4);
|
||||
Gamma g5 = Gamma(Gamma::Algebra::Gamma5);
|
||||
LatticeComplex expbuf(env().getGrid());
|
||||
std::vector<TComplex> corrbuf;
|
||||
std::vector<Result> result(n_eye_diag);
|
||||
unsigned int ndim = env().getNd();
|
||||
|
||||
PropagatorField tmp1(env().getGrid());
|
||||
LatticeComplex tmp2(env().getGrid());
|
||||
std::vector<PropagatorField> S_body(ndim, tmp1);
|
||||
std::vector<PropagatorField> S_loop(ndim, tmp1);
|
||||
std::vector<LatticeComplex> E_body(ndim, tmp2);
|
||||
std::vector<LatticeComplex> E_loop(ndim, tmp2);
|
||||
|
||||
// Setup for S-type contractions.
|
||||
for (int mu = 0; mu < ndim; ++mu)
|
||||
{
|
||||
S_body[mu] = MAKE_SE_BODY(q1, q2, q3, GammaL(Gamma::gmu[mu]));
|
||||
S_loop[mu] = MAKE_SE_LOOP(q4, GammaL(Gamma::gmu[mu]));
|
||||
}
|
||||
|
||||
// Perform S-type contractions.
|
||||
SUM_MU(expbuf, trace(S_body[mu]*S_loop[mu]))
|
||||
MAKE_DIAG(expbuf, corrbuf, result[S_diag], "HW_S")
|
||||
|
||||
// Recycle sub-expressions for E-type contractions.
|
||||
for (unsigned int mu = 0; mu < ndim; ++mu)
|
||||
{
|
||||
E_body[mu] = trace(S_body[mu]);
|
||||
E_loop[mu] = trace(S_loop[mu]);
|
||||
}
|
||||
|
||||
// Perform E-type contractions.
|
||||
SUM_MU(expbuf, E_body[mu]*E_loop[mu])
|
||||
MAKE_DIAG(expbuf, corrbuf, result[E_diag], "HW_E")
|
||||
|
||||
write(writer, "HW_Eye", result);
|
||||
}
|
58
extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp
Normal file
58
extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp
Normal file
@ -0,0 +1,58 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp
|
||||
|
||||
Copyright (C) 2017
|
||||
|
||||
Author: Andrew Lawson <andrew.lawson1991@gmail.com>
|
||||
|
||||
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_WeakHamiltonianEye_hpp_
|
||||
#define Hadrons_WeakHamiltonianEye_hpp_
|
||||
|
||||
#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
|
||||
|
||||
BEGIN_HADRONS_NAMESPACE
|
||||
|
||||
/******************************************************************************
|
||||
* WeakHamiltonianEye *
|
||||
******************************************************************************/
|
||||
BEGIN_MODULE_NAMESPACE(MContraction)
|
||||
|
||||
enum
|
||||
{
|
||||
S_diag = 0,
|
||||
E_diag = 1,
|
||||
n_eye_diag = 2
|
||||
};
|
||||
|
||||
// Saucer and Eye subdiagram contractions.
|
||||
#define MAKE_SE_BODY(Q_1, Q_2, Q_3, gamma) (Q_3*g5*Q_1*adj(Q_2)*g5*gamma)
|
||||
#define MAKE_SE_LOOP(Q_loop, gamma) (Q_loop*gamma)
|
||||
|
||||
MAKE_WEAK_MODULE(WeakHamiltonianEye)
|
||||
|
||||
END_MODULE_NAMESPACE
|
||||
|
||||
END_HADRONS_NAMESPACE
|
||||
|
||||
#endif // Hadrons_WeakHamiltonianEye_hpp_
|
139
extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc
Normal file
139
extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc
Normal file
@ -0,0 +1,139 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc
|
||||
|
||||
Copyright (C) 2017
|
||||
|
||||
Author: Andrew Lawson <andrew.lawson1991@gmail.com>
|
||||
|
||||
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 <Grid/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp>
|
||||
|
||||
using namespace Grid;
|
||||
using namespace Hadrons;
|
||||
using namespace MContraction;
|
||||
|
||||
/*
|
||||
* Weak Hamiltonian current-current contractions, Non-Eye-type.
|
||||
*
|
||||
* These contractions are generated by the Q1 and Q2 operators in the physical
|
||||
* basis (see e.g. Fig 3 of arXiv:1507.03094).
|
||||
*
|
||||
* Schematic:
|
||||
* q2 q3 | q2 q3
|
||||
* /--<--¬ /--<--¬ | /--<--¬ /--<--¬
|
||||
* / \ / \ | / \ / \
|
||||
* / \ / \ | / \ / \
|
||||
* / \ / \ | / \ / \
|
||||
* i * * H_W * f | i * * * H_W * f
|
||||
* \ * | | \ / \ /
|
||||
* \ / \ / | \ / \ /
|
||||
* \ / \ / | \ / \ /
|
||||
* \ / \ / | \-->--/ \-->--/
|
||||
* \-->--/ \-->--/ | q1 q4
|
||||
* q1 q4 |
|
||||
* Connected (C) | Wing (W)
|
||||
*
|
||||
* C: trace(q1*adj(q2)*g5*gL[mu]*q3*adj(q4)*g5*gL[mu])
|
||||
* W: trace(q1*adj(q2)*g5*gL[mu])*trace(q3*adj(q4)*g5*gL[mu])
|
||||
*
|
||||
*/
|
||||
|
||||
/******************************************************************************
|
||||
* TWeakHamiltonianNonEye implementation *
|
||||
******************************************************************************/
|
||||
// constructor /////////////////////////////////////////////////////////////////
|
||||
TWeakHamiltonianNonEye::TWeakHamiltonianNonEye(const std::string name)
|
||||
: Module<WeakHamiltonianPar>(name)
|
||||
{}
|
||||
|
||||
// dependencies/products ///////////////////////////////////////////////////////
|
||||
std::vector<std::string> TWeakHamiltonianNonEye::getInput(void)
|
||||
{
|
||||
std::vector<std::string> in = {par().q1, par().q2, par().q3, par().q4};
|
||||
|
||||
return in;
|
||||
}
|
||||
|
||||
std::vector<std::string> TWeakHamiltonianNonEye::getOutput(void)
|
||||
{
|
||||
std::vector<std::string> out = {getName()};
|
||||
|
||||
return out;
|
||||
}
|
||||
|
||||
// setup ///////////////////////////////////////////////////////////////////////
|
||||
void TWeakHamiltonianNonEye::setup(void)
|
||||
{
|
||||
|
||||
}
|
||||
|
||||
// execution ///////////////////////////////////////////////////////////////////
|
||||
void TWeakHamiltonianNonEye::execute(void)
|
||||
{
|
||||
LOG(Message) << "Computing Weak Hamiltonian (Non-Eye type) contractions '"
|
||||
<< getName() << "' using quarks '" << par().q1 << "', '"
|
||||
<< par().q2 << ", '" << par().q3 << "' and '" << par().q4
|
||||
<< "'." << std::endl;
|
||||
|
||||
CorrWriter writer(par().output);
|
||||
PropagatorField &q1 = *env().template getObject<PropagatorField>(par().q1);
|
||||
PropagatorField &q2 = *env().template getObject<PropagatorField>(par().q2);
|
||||
PropagatorField &q3 = *env().template getObject<PropagatorField>(par().q3);
|
||||
PropagatorField &q4 = *env().template getObject<PropagatorField>(par().q4);
|
||||
Gamma g5 = Gamma(Gamma::Algebra::Gamma5);
|
||||
LatticeComplex expbuf(env().getGrid());
|
||||
std::vector<TComplex> corrbuf;
|
||||
std::vector<Result> result(n_noneye_diag);
|
||||
unsigned int ndim = env().getNd();
|
||||
|
||||
PropagatorField tmp1(env().getGrid());
|
||||
LatticeComplex tmp2(env().getGrid());
|
||||
std::vector<PropagatorField> C_i_side_loop(ndim, tmp1);
|
||||
std::vector<PropagatorField> C_f_side_loop(ndim, tmp1);
|
||||
std::vector<LatticeComplex> W_i_side_loop(ndim, tmp2);
|
||||
std::vector<LatticeComplex> W_f_side_loop(ndim, tmp2);
|
||||
|
||||
// Setup for C-type contractions.
|
||||
for (int mu = 0; mu < ndim; ++mu)
|
||||
{
|
||||
C_i_side_loop[mu] = MAKE_CW_SUBDIAG(q1, q2, GammaL(Gamma::gmu[mu]));
|
||||
C_f_side_loop[mu] = MAKE_CW_SUBDIAG(q3, q4, GammaL(Gamma::gmu[mu]));
|
||||
}
|
||||
|
||||
// Perform C-type contractions.
|
||||
SUM_MU(expbuf, trace(C_i_side_loop[mu]*C_f_side_loop[mu]))
|
||||
MAKE_DIAG(expbuf, corrbuf, result[C_diag], "HW_C")
|
||||
|
||||
// Recycle sub-expressions for W-type contractions.
|
||||
for (unsigned int mu = 0; mu < ndim; ++mu)
|
||||
{
|
||||
W_i_side_loop[mu] = trace(C_i_side_loop[mu]);
|
||||
W_f_side_loop[mu] = trace(C_f_side_loop[mu]);
|
||||
}
|
||||
|
||||
// Perform W-type contractions.
|
||||
SUM_MU(expbuf, W_i_side_loop[mu]*W_f_side_loop[mu])
|
||||
MAKE_DIAG(expbuf, corrbuf, result[W_diag], "HW_W")
|
||||
|
||||
write(writer, "HW_NonEye", result);
|
||||
}
|
@ -0,0 +1,57 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp
|
||||
|
||||
Copyright (C) 2017
|
||||
|
||||
Author: Andrew Lawson <andrew.lawson1991@gmail.com>
|
||||
|
||||
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_WeakHamiltonianNonEye_hpp_
|
||||
#define Hadrons_WeakHamiltonianNonEye_hpp_
|
||||
|
||||
#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
|
||||
|
||||
BEGIN_HADRONS_NAMESPACE
|
||||
|
||||
/******************************************************************************
|
||||
* WeakHamiltonianNonEye *
|
||||
******************************************************************************/
|
||||
BEGIN_MODULE_NAMESPACE(MContraction)
|
||||
|
||||
enum
|
||||
{
|
||||
W_diag = 0,
|
||||
C_diag = 1,
|
||||
n_noneye_diag = 2
|
||||
};
|
||||
|
||||
// Wing and Connected subdiagram contractions
|
||||
#define MAKE_CW_SUBDIAG(Q_1, Q_2, gamma) (Q_1*adj(Q_2)*g5*gamma)
|
||||
|
||||
MAKE_WEAK_MODULE(WeakHamiltonianNonEye)
|
||||
|
||||
END_MODULE_NAMESPACE
|
||||
|
||||
END_HADRONS_NAMESPACE
|
||||
|
||||
#endif // Hadrons_WeakHamiltonianNonEye_hpp_
|
135
extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.cc
Normal file
135
extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.cc
Normal file
@ -0,0 +1,135 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.cc
|
||||
|
||||
Copyright (C) 2017
|
||||
|
||||
Author: Andrew Lawson <andrew.lawson1991@gmail.com>
|
||||
|
||||
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 <Grid/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp>
|
||||
|
||||
using namespace Grid;
|
||||
using namespace Hadrons;
|
||||
using namespace MContraction;
|
||||
|
||||
/*
|
||||
* Weak Hamiltonian + current contractions, disconnected topology for neutral
|
||||
* mesons.
|
||||
*
|
||||
* These contractions are generated by operators Q_1,...,10 of the dS=1 Weak
|
||||
* Hamiltonian in the physical basis and an additional current J (see e.g.
|
||||
* Fig 11 of arXiv:1507.03094).
|
||||
*
|
||||
* Schematic:
|
||||
*
|
||||
* q2 q4 q3
|
||||
* /--<--¬ /---<--¬ /---<--¬
|
||||
* / \ / \ / \
|
||||
* i * * H_W | J * * f
|
||||
* \ / \ / \ /
|
||||
* \--->---/ \-------/ \------/
|
||||
* q1
|
||||
*
|
||||
* options
|
||||
* - q1: input propagator 1 (string)
|
||||
* - q2: input propagator 2 (string)
|
||||
* - q3: input propagator 3 (string), assumed to be sequential propagator
|
||||
* - q4: input propagator 4 (string), assumed to be a loop
|
||||
*
|
||||
* type 1: trace(q1*adj(q2)*g5*gL[mu])*trace(loop*gL[mu])*trace(q3*g5)
|
||||
* type 2: trace(q1*adj(q2)*g5*gL[mu]*loop*gL[mu])*trace(q3*g5)
|
||||
*/
|
||||
|
||||
/*******************************************************************************
|
||||
* TWeakNeutral4ptDisc implementation *
|
||||
******************************************************************************/
|
||||
// constructor /////////////////////////////////////////////////////////////////
|
||||
TWeakNeutral4ptDisc::TWeakNeutral4ptDisc(const std::string name)
|
||||
: Module<WeakHamiltonianPar>(name)
|
||||
{}
|
||||
|
||||
// dependencies/products ///////////////////////////////////////////////////////
|
||||
std::vector<std::string> TWeakNeutral4ptDisc::getInput(void)
|
||||
{
|
||||
std::vector<std::string> in = {par().q1, par().q2, par().q3, par().q4};
|
||||
|
||||
return in;
|
||||
}
|
||||
|
||||
std::vector<std::string> TWeakNeutral4ptDisc::getOutput(void)
|
||||
{
|
||||
std::vector<std::string> out = {getName()};
|
||||
|
||||
return out;
|
||||
}
|
||||
|
||||
// setup ///////////////////////////////////////////////////////////////////////
|
||||
void TWeakNeutral4ptDisc::setup(void)
|
||||
{
|
||||
|
||||
}
|
||||
|
||||
// execution ///////////////////////////////////////////////////////////////////
|
||||
void TWeakNeutral4ptDisc::execute(void)
|
||||
{
|
||||
LOG(Message) << "Computing Weak Hamiltonian neutral disconnected contractions '"
|
||||
<< getName() << "' using quarks '" << par().q1 << "', '"
|
||||
<< par().q2 << ", '" << par().q3 << "' and '" << par().q4
|
||||
<< "'." << std::endl;
|
||||
|
||||
CorrWriter writer(par().output);
|
||||
PropagatorField &q1 = *env().template getObject<PropagatorField>(par().q1);
|
||||
PropagatorField &q2 = *env().template getObject<PropagatorField>(par().q2);
|
||||
PropagatorField &q3 = *env().template getObject<PropagatorField>(par().q3);
|
||||
PropagatorField &q4 = *env().template getObject<PropagatorField>(par().q4);
|
||||
Gamma g5 = Gamma(Gamma::Algebra::Gamma5);
|
||||
LatticeComplex expbuf(env().getGrid());
|
||||
std::vector<TComplex> corrbuf;
|
||||
std::vector<Result> result(n_neut_disc_diag);
|
||||
unsigned int ndim = env().getNd();
|
||||
|
||||
PropagatorField tmp(env().getGrid());
|
||||
std::vector<PropagatorField> meson(ndim, tmp);
|
||||
std::vector<PropagatorField> loop(ndim, tmp);
|
||||
LatticeComplex curr(env().getGrid());
|
||||
|
||||
// Setup for type 1 contractions.
|
||||
for (int mu = 0; mu < ndim; ++mu)
|
||||
{
|
||||
meson[mu] = MAKE_DISC_MESON(q1, q2, GammaL(Gamma::gmu[mu]));
|
||||
loop[mu] = MAKE_DISC_LOOP(q4, GammaL(Gamma::gmu[mu]));
|
||||
}
|
||||
curr = MAKE_DISC_CURR(q3, GammaL(Gamma::Algebra::Gamma5));
|
||||
|
||||
// Perform type 1 contractions.
|
||||
SUM_MU(expbuf, trace(meson[mu]*loop[mu]))
|
||||
expbuf *= curr;
|
||||
MAKE_DIAG(expbuf, corrbuf, result[neut_disc_1_diag], "HW_disc0_1")
|
||||
|
||||
// Perform type 2 contractions.
|
||||
SUM_MU(expbuf, trace(meson[mu])*trace(loop[mu]))
|
||||
expbuf *= curr;
|
||||
MAKE_DIAG(expbuf, corrbuf, result[neut_disc_2_diag], "HW_disc0_2")
|
||||
|
||||
write(writer, "HW_disc0", result);
|
||||
}
|
59
extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp
Normal file
59
extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp
Normal file
@ -0,0 +1,59 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp
|
||||
|
||||
Copyright (C) 2017
|
||||
|
||||
Author: Andrew Lawson <andrew.lawson1991@gmail.com>
|
||||
|
||||
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_WeakNeutral4ptDisc_hpp_
|
||||
#define Hadrons_WeakNeutral4ptDisc_hpp_
|
||||
|
||||
#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
|
||||
|
||||
BEGIN_HADRONS_NAMESPACE
|
||||
|
||||
/******************************************************************************
|
||||
* WeakNeutral4ptDisc *
|
||||
******************************************************************************/
|
||||
BEGIN_MODULE_NAMESPACE(MContraction)
|
||||
|
||||
enum
|
||||
{
|
||||
neut_disc_1_diag = 0,
|
||||
neut_disc_2_diag = 1,
|
||||
n_neut_disc_diag = 2
|
||||
};
|
||||
|
||||
// Neutral 4pt disconnected subdiagram contractions.
|
||||
#define MAKE_DISC_MESON(Q_1, Q_2, gamma) (Q_1*adj(Q_2)*g5*gamma)
|
||||
#define MAKE_DISC_LOOP(Q_LOOP, gamma) (Q_LOOP*gamma)
|
||||
#define MAKE_DISC_CURR(Q_c, gamma) (trace(Q_c*gamma))
|
||||
|
||||
MAKE_WEAK_MODULE(WeakNeutral4ptDisc)
|
||||
|
||||
END_MODULE_NAMESPACE
|
||||
|
||||
END_HADRONS_NAMESPACE
|
||||
|
||||
#endif // Hadrons_WeakNeutral4ptDisc_hpp_
|
132
extras/Hadrons/Modules/MLoop/NoiseLoop.hpp
Normal file
132
extras/Hadrons/Modules/MLoop/NoiseLoop.hpp
Normal file
@ -0,0 +1,132 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: extras/Hadrons/Modules/MLoop/NoiseLoop.hpp
|
||||
|
||||
Copyright (C) 2016
|
||||
|
||||
Author: Andrew Lawson <andrew.lawson1991@gmail.com>
|
||||
|
||||
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_NoiseLoop_hpp_
|
||||
#define Hadrons_NoiseLoop_hpp_
|
||||
|
||||
#include <Grid/Hadrons/Global.hpp>
|
||||
#include <Grid/Hadrons/Module.hpp>
|
||||
#include <Grid/Hadrons/ModuleFactory.hpp>
|
||||
|
||||
BEGIN_HADRONS_NAMESPACE
|
||||
|
||||
/*
|
||||
|
||||
Noise loop propagator
|
||||
-----------------------------
|
||||
* loop_x = q_x * adj(eta_x)
|
||||
|
||||
* options:
|
||||
- q = Result of inversion on noise source.
|
||||
- eta = noise source.
|
||||
|
||||
*/
|
||||
|
||||
|
||||
/******************************************************************************
|
||||
* NoiseLoop *
|
||||
******************************************************************************/
|
||||
BEGIN_MODULE_NAMESPACE(MLoop)
|
||||
|
||||
class NoiseLoopPar: Serializable
|
||||
{
|
||||
public:
|
||||
GRID_SERIALIZABLE_CLASS_MEMBERS(NoiseLoopPar,
|
||||
std::string, q,
|
||||
std::string, eta);
|
||||
};
|
||||
|
||||
template <typename FImpl>
|
||||
class TNoiseLoop: public Module<NoiseLoopPar>
|
||||
{
|
||||
public:
|
||||
TYPE_ALIASES(FImpl,);
|
||||
public:
|
||||
// constructor
|
||||
TNoiseLoop(const std::string name);
|
||||
// destructor
|
||||
virtual ~TNoiseLoop(void) = default;
|
||||
// dependency relation
|
||||
virtual std::vector<std::string> getInput(void);
|
||||
virtual std::vector<std::string> getOutput(void);
|
||||
// setup
|
||||
virtual void setup(void);
|
||||
// execution
|
||||
virtual void execute(void);
|
||||
};
|
||||
|
||||
MODULE_REGISTER_NS(NoiseLoop, TNoiseLoop<FIMPL>, MLoop);
|
||||
|
||||
/******************************************************************************
|
||||
* TNoiseLoop implementation *
|
||||
******************************************************************************/
|
||||
// constructor /////////////////////////////////////////////////////////////////
|
||||
template <typename FImpl>
|
||||
TNoiseLoop<FImpl>::TNoiseLoop(const std::string name)
|
||||
: Module<NoiseLoopPar>(name)
|
||||
{}
|
||||
|
||||
// dependencies/products ///////////////////////////////////////////////////////
|
||||
template <typename FImpl>
|
||||
std::vector<std::string> TNoiseLoop<FImpl>::getInput(void)
|
||||
{
|
||||
std::vector<std::string> in = {par().q, par().eta};
|
||||
|
||||
return in;
|
||||
}
|
||||
|
||||
template <typename FImpl>
|
||||
std::vector<std::string> TNoiseLoop<FImpl>::getOutput(void)
|
||||
{
|
||||
std::vector<std::string> out = {getName()};
|
||||
|
||||
return out;
|
||||
}
|
||||
|
||||
// setup ///////////////////////////////////////////////////////////////////////
|
||||
template <typename FImpl>
|
||||
void TNoiseLoop<FImpl>::setup(void)
|
||||
{
|
||||
env().template registerLattice<PropagatorField>(getName());
|
||||
}
|
||||
|
||||
// execution ///////////////////////////////////////////////////////////////////
|
||||
template <typename FImpl>
|
||||
void TNoiseLoop<FImpl>::execute(void)
|
||||
{
|
||||
PropagatorField &loop = *env().template createLattice<PropagatorField>(getName());
|
||||
PropagatorField &q = *env().template getObject<PropagatorField>(par().q);
|
||||
PropagatorField &eta = *env().template getObject<PropagatorField>(par().eta);
|
||||
loop = q*adj(eta);
|
||||
}
|
||||
|
||||
END_MODULE_NAMESPACE
|
||||
|
||||
END_HADRONS_NAMESPACE
|
||||
|
||||
#endif // Hadrons_NoiseLoop_hpp_
|
@ -6,6 +6,7 @@ Source file: extras/Hadrons/Modules/MSource/SeqGamma.hpp
|
||||
|
||||
Copyright (C) 2015
|
||||
Copyright (C) 2016
|
||||
Copyright (C) 2017
|
||||
|
||||
Author: Antonin Portelli <antonin.portelli@me.com>
|
||||
|
||||
@ -149,9 +150,9 @@ void TSeqGamma<FImpl>::execute(void)
|
||||
for(unsigned int mu = 0; mu < env().getNd(); mu++)
|
||||
{
|
||||
LatticeCoordinate(coor, mu);
|
||||
ph = ph + p[mu]*coor;
|
||||
ph = ph + p[mu]*coor*((1./(env().getGrid()->_fdimensions[mu])));
|
||||
}
|
||||
ph = exp(i*ph);
|
||||
ph = exp((Real)(2*M_PI)*i*ph);
|
||||
LatticeCoordinate(t, Tp);
|
||||
src = where((t >= par().tA) and (t <= par().tB), ph*(g*q), 0.*q);
|
||||
}
|
||||
|
147
extras/Hadrons/Modules/MSource/Wall.hpp
Normal file
147
extras/Hadrons/Modules/MSource/Wall.hpp
Normal file
@ -0,0 +1,147 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: extras/Hadrons/Modules/MSource/Wall.hpp
|
||||
|
||||
Copyright (C) 2017
|
||||
|
||||
Author: Andrew Lawson <andrew.lawson1991@gmail.com>
|
||||
|
||||
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_WallSource_hpp_
|
||||
#define Hadrons_WallSource_hpp_
|
||||
|
||||
#include <Grid/Hadrons/Global.hpp>
|
||||
#include <Grid/Hadrons/Module.hpp>
|
||||
#include <Grid/Hadrons/ModuleFactory.hpp>
|
||||
|
||||
BEGIN_HADRONS_NAMESPACE
|
||||
|
||||
/*
|
||||
|
||||
Wall source
|
||||
-----------------------------
|
||||
* src_x = delta(x_3 - tW) * exp(i x.mom)
|
||||
|
||||
* options:
|
||||
- tW: source timeslice (integer)
|
||||
- mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0.")
|
||||
|
||||
*/
|
||||
|
||||
/******************************************************************************
|
||||
* Wall *
|
||||
******************************************************************************/
|
||||
BEGIN_MODULE_NAMESPACE(MSource)
|
||||
|
||||
class WallPar: Serializable
|
||||
{
|
||||
public:
|
||||
GRID_SERIALIZABLE_CLASS_MEMBERS(WallPar,
|
||||
unsigned int, tW,
|
||||
std::string, mom);
|
||||
};
|
||||
|
||||
template <typename FImpl>
|
||||
class TWall: public Module<WallPar>
|
||||
{
|
||||
public:
|
||||
TYPE_ALIASES(FImpl,);
|
||||
public:
|
||||
// constructor
|
||||
TWall(const std::string name);
|
||||
// destructor
|
||||
virtual ~TWall(void) = default;
|
||||
// dependency relation
|
||||
virtual std::vector<std::string> getInput(void);
|
||||
virtual std::vector<std::string> getOutput(void);
|
||||
// setup
|
||||
virtual void setup(void);
|
||||
// execution
|
||||
virtual void execute(void);
|
||||
};
|
||||
|
||||
MODULE_REGISTER_NS(Wall, TWall<FIMPL>, MSource);
|
||||
|
||||
/******************************************************************************
|
||||
* TWall implementation *
|
||||
******************************************************************************/
|
||||
// constructor /////////////////////////////////////////////////////////////////
|
||||
template <typename FImpl>
|
||||
TWall<FImpl>::TWall(const std::string name)
|
||||
: Module<WallPar>(name)
|
||||
{}
|
||||
|
||||
// dependencies/products ///////////////////////////////////////////////////////
|
||||
template <typename FImpl>
|
||||
std::vector<std::string> TWall<FImpl>::getInput(void)
|
||||
{
|
||||
std::vector<std::string> in;
|
||||
|
||||
return in;
|
||||
}
|
||||
|
||||
template <typename FImpl>
|
||||
std::vector<std::string> TWall<FImpl>::getOutput(void)
|
||||
{
|
||||
std::vector<std::string> out = {getName()};
|
||||
|
||||
return out;
|
||||
}
|
||||
|
||||
// setup ///////////////////////////////////////////////////////////////////////
|
||||
template <typename FImpl>
|
||||
void TWall<FImpl>::setup(void)
|
||||
{
|
||||
env().template registerLattice<PropagatorField>(getName());
|
||||
}
|
||||
|
||||
// execution ///////////////////////////////////////////////////////////////////
|
||||
template <typename FImpl>
|
||||
void TWall<FImpl>::execute(void)
|
||||
{
|
||||
LOG(Message) << "Generating wall source at t = " << par().tW
|
||||
<< " with momentum " << par().mom << std::endl;
|
||||
|
||||
PropagatorField &src = *env().template createLattice<PropagatorField>(getName());
|
||||
Lattice<iScalar<vInteger>> t(env().getGrid());
|
||||
LatticeComplex ph(env().getGrid()), coor(env().getGrid());
|
||||
std::vector<Real> p;
|
||||
Complex i(0.0,1.0);
|
||||
|
||||
p = strToVec<Real>(par().mom);
|
||||
ph = zero;
|
||||
for(unsigned int mu = 0; mu < Nd; mu++)
|
||||
{
|
||||
LatticeCoordinate(coor, mu);
|
||||
ph = ph + p[mu]*coor*((1./(env().getGrid()->_fdimensions[mu])));
|
||||
}
|
||||
ph = exp((Real)(2*M_PI)*i*ph);
|
||||
LatticeCoordinate(t, Tp);
|
||||
src = 1.;
|
||||
src = where((t == par().tW), src*ph, 0.*src);
|
||||
}
|
||||
|
||||
END_MODULE_NAMESPACE
|
||||
|
||||
END_HADRONS_NAMESPACE
|
||||
|
||||
#endif // Hadrons_WallSource_hpp_
|
@ -173,7 +173,7 @@ void TQuark<FImpl>::execute(void)
|
||||
*env().template getObject<PropagatorField>(getName());
|
||||
|
||||
axpby_ssp_pminus(sol, 0., sol, 1., sol, 0, 0);
|
||||
axpby_ssp_pplus(sol, 0., sol, 1., sol, 0, Ls_-1);
|
||||
axpby_ssp_pplus(sol, 1., sol, 1., sol, 0, Ls_-1);
|
||||
ExtractSlice(tmp, sol, 0, 0);
|
||||
FermToProp(p4d, tmp, s, c);
|
||||
}
|
||||
|
@ -1,4 +1,7 @@
|
||||
modules_cc =\
|
||||
Modules/MContraction/WeakHamiltonianEye.cc \
|
||||
Modules/MContraction/WeakHamiltonianNonEye.cc \
|
||||
Modules/MContraction/WeakNeutral4ptDisc.cc \
|
||||
Modules/MGauge/Load.cc \
|
||||
Modules/MGauge/Random.cc \
|
||||
Modules/MGauge/Unit.cc
|
||||
@ -7,13 +10,21 @@ modules_hpp =\
|
||||
Modules/MAction/DWF.hpp \
|
||||
Modules/MAction/Wilson.hpp \
|
||||
Modules/MContraction/Baryon.hpp \
|
||||
Modules/MContraction/DiscLoop.hpp \
|
||||
Modules/MContraction/Gamma3pt.hpp \
|
||||
Modules/MContraction/Meson.hpp \
|
||||
Modules/MContraction/WeakHamiltonian.hpp \
|
||||
Modules/MContraction/WeakHamiltonianEye.hpp \
|
||||
Modules/MContraction/WeakHamiltonianNonEye.hpp \
|
||||
Modules/MContraction/WeakNeutral4ptDisc.hpp \
|
||||
Modules/MGauge/Load.hpp \
|
||||
Modules/MGauge/Random.hpp \
|
||||
Modules/MGauge/Unit.hpp \
|
||||
Modules/MLoop/NoiseLoop.hpp \
|
||||
Modules/MSolver/RBPrecCG.hpp \
|
||||
Modules/MSource/Point.hpp \
|
||||
Modules/MSource/SeqGamma.hpp \
|
||||
Modules/MSource/Wall.hpp \
|
||||
Modules/MSource/Z2.hpp \
|
||||
Modules/Quark.hpp
|
||||
|
||||
|
@ -138,6 +138,54 @@ namespace Grid {
|
||||
unsigned int dimInd_{0};
|
||||
};
|
||||
|
||||
// Pair IO utilities /////////////////////////////////////////////////////////
|
||||
// helper function to parse input in the format "<obj1 obj2>"
|
||||
template <typename T1, typename T2>
|
||||
inline std::istream & operator>>(std::istream &is, std::pair<T1, T2> &buf)
|
||||
{
|
||||
T1 buf1;
|
||||
T2 buf2;
|
||||
char c;
|
||||
|
||||
// Search for "pair" delimiters.
|
||||
do
|
||||
{
|
||||
is.get(c);
|
||||
} while (c != '<' && !is.eof());
|
||||
if (c == '<')
|
||||
{
|
||||
int start = is.tellg();
|
||||
do
|
||||
{
|
||||
is.get(c);
|
||||
} while (c != '>' && !is.eof());
|
||||
if (c == '>')
|
||||
{
|
||||
int end = is.tellg();
|
||||
int psize = end - start - 1;
|
||||
|
||||
// Only read data between pair limiters.
|
||||
is.seekg(start);
|
||||
std::string tmpstr(psize, ' ');
|
||||
is.read(&tmpstr[0], psize);
|
||||
std::istringstream temp(tmpstr);
|
||||
temp >> buf1 >> buf2;
|
||||
buf = std::make_pair(buf1, buf2);
|
||||
is.seekg(end);
|
||||
}
|
||||
}
|
||||
is.peek();
|
||||
return is;
|
||||
}
|
||||
|
||||
// output to streams for pairs
|
||||
template <class T1, class T2>
|
||||
inline std::ostream & operator<<(std::ostream &os, const std::pair<T1, T2> &p)
|
||||
{
|
||||
os << "<" << p.first << " " << p.second << ">";
|
||||
return os;
|
||||
}
|
||||
|
||||
// Abstract writer/reader classes ////////////////////////////////////////////
|
||||
// static polymorphism implemented using CRTP idiom
|
||||
class Serializable;
|
||||
|
@ -113,6 +113,7 @@ int main(int argc,char **argv)
|
||||
// test serializable class writing
|
||||
myclass obj(1234); // non-trivial constructor
|
||||
std::vector<myclass> vec;
|
||||
std::pair<myenum, myenum> pair;
|
||||
|
||||
std::cout << "-- serialisable class writing to 'bother.xml'..." << std::endl;
|
||||
write(WR,"obj",obj);
|
||||
@ -120,6 +121,8 @@ int main(int argc,char **argv)
|
||||
vec.push_back(myclass(1234));
|
||||
vec.push_back(myclass(5678));
|
||||
vec.push_back(myclass(3838));
|
||||
pair = std::make_pair(myenum::red, myenum::blue);
|
||||
|
||||
write(WR, "objvec", vec);
|
||||
std::cout << "-- serialisable class writing to std::cout:" << std::endl;
|
||||
std::cout << obj << std::endl;
|
||||
@ -127,21 +130,30 @@ int main(int argc,char **argv)
|
||||
std::cout << "vec[0] == obj: " << ((vec[0] == obj) ? "true" : "false") << std::endl;
|
||||
std::cout << "vec[1] == obj: " << ((vec[1] == obj) ? "true" : "false") << std::endl;
|
||||
|
||||
write(WR, "objpair", pair);
|
||||
std::cout << "-- pair writing to std::cout:" << std::endl;
|
||||
std::cout << pair << std::endl;
|
||||
|
||||
// read tests
|
||||
std::cout << "\n==== IO self-consistency tests" << std::endl;
|
||||
//// XML
|
||||
ioTest<XmlWriter, XmlReader>("iotest.xml", obj, "XML (object) ");
|
||||
ioTest<XmlWriter, XmlReader>("iotest.xml", vec, "XML (vector of objects)");
|
||||
ioTest<XmlWriter, XmlReader>("iotest.xml", pair, "XML (pair of objects)");
|
||||
//// binary
|
||||
ioTest<BinaryWriter, BinaryReader>("iotest.bin", obj, "binary (object) ");
|
||||
ioTest<BinaryWriter, BinaryReader>("iotest.bin", vec, "binary (vector of objects)");
|
||||
ioTest<BinaryWriter, BinaryReader>("iotest.bin", pair, "binary (pair of objects)");
|
||||
//// text
|
||||
ioTest<TextWriter, TextReader>("iotest.dat", obj, "text (object) ");
|
||||
ioTest<TextWriter, TextReader>("iotest.dat", vec, "text (vector of objects)");
|
||||
ioTest<TextWriter, TextReader>("iotest.dat", pair, "text (pair of objects)");
|
||||
//// HDF5
|
||||
#undef HAVE_HDF5
|
||||
#ifdef HAVE_HDF5
|
||||
ioTest<Hdf5Writer, Hdf5Reader>("iotest.h5", obj, "HDF5 (object) ");
|
||||
ioTest<Hdf5Writer, Hdf5Reader>("iotest.h5", vec, "HDF5 (vector of objects)");
|
||||
ioTest<Hdf5Writer, Hdf5Reader>("iotest.h5", pair, "HDF5 (pair of objects)");
|
||||
#endif
|
||||
|
||||
std::cout << "\n==== vector flattening/reconstruction" << std::endl;
|
||||
|
368
tests/hadrons/Test_hadrons.hpp
Normal file
368
tests/hadrons/Test_hadrons.hpp
Normal file
@ -0,0 +1,368 @@
|
||||
/*******************************************************************************
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: tests/hadrons/Test_hadrons.hpp
|
||||
|
||||
Copyright (C) 2017
|
||||
|
||||
Author: Andrew Lawson <andrew.lawson1991@gmail.com>
|
||||
|
||||
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 <Grid/Hadrons/Application.hpp>
|
||||
|
||||
using namespace Grid;
|
||||
using namespace Hadrons;
|
||||
|
||||
/*******************************************************************************
|
||||
* Macros to reduce code duplication.
|
||||
******************************************************************************/
|
||||
// Useful definitions
|
||||
#define ZERO_MOM "0. 0. 0. 0."
|
||||
#define INIT_INDEX(s, n) (std::string(s) + "_" + std::to_string(n))
|
||||
#define ADD_INDEX(s, n) (s + "_" + std::to_string(n))
|
||||
#define LABEL_3PT(s, t1, t2) ADD_INDEX(INIT_INDEX(s, t1), t2)
|
||||
#define LABEL_4PT(s, t1, t2, t3) ADD_INDEX(ADD_INDEX(INIT_INDEX(s, t1), t2), t3)
|
||||
#define LABEL_4PT_NOISE(s, t1, t2, t3, nn) ADD_INDEX(ADD_INDEX(ADD_INDEX(INIT_INDEX(s, t1), t2), t3), nn)
|
||||
|
||||
// Wall source/sink macros
|
||||
#define NAME_3MOM_WALL_SOURCE(t, mom) ("wall_" + std::to_string(t) + "_" + mom)
|
||||
#define NAME_WALL_SOURCE(t) NAME_3MOM_WALL_SOURCE(t, ZERO_MOM)
|
||||
#define NAME_POINT_SOURCE(pos) ("point_" + pos)
|
||||
|
||||
#define MAKE_3MOM_WALL_PROP(tW, mom, propName, solver)\
|
||||
{\
|
||||
std::string srcName = NAME_3MOM_WALL_SOURCE(tW, mom);\
|
||||
makeWallSource(application, srcName, tW, mom);\
|
||||
makePropagator(application, propName, srcName, solver);\
|
||||
}
|
||||
|
||||
#define MAKE_WALL_PROP(tW, propName, solver)\
|
||||
MAKE_3MOM_WALL_PROP(tW, ZERO_MOM, propName, solver)
|
||||
|
||||
// Sequential source macros
|
||||
#define MAKE_SEQUENTIAL_PROP(tS, qSrc, mom, propName, solver)\
|
||||
{\
|
||||
std::string srcName = ADD_INDEX(qSrc + "_seq", tS);\
|
||||
makeSequentialSource(application, srcName, qSrc, tS, mom);\
|
||||
makePropagator(application, propName, srcName, solver);\
|
||||
}
|
||||
|
||||
// Point source macros
|
||||
#define MAKE_POINT_PROP(pos, propName, solver)\
|
||||
{\
|
||||
std::string srcName = NAME_POINT_SOURCE(pos);\
|
||||
makePointSource(application, srcName, pos);\
|
||||
makePropagator(application, propName, srcName, solver);\
|
||||
}
|
||||
|
||||
/*******************************************************************************
|
||||
* Functions for propagator construction.
|
||||
******************************************************************************/
|
||||
|
||||
/*******************************************************************************
|
||||
* Name: makePointSource
|
||||
* Purpose: Construct point source and add to application module.
|
||||
* Parameters: application - main application that stores modules.
|
||||
* srcName - name of source module to create.
|
||||
* pos - Position of point source.
|
||||
* Returns: None.
|
||||
******************************************************************************/
|
||||
inline void makePointSource(Application &application, std::string srcName,
|
||||
std::string pos)
|
||||
{
|
||||
// If the source already exists, don't make the module again.
|
||||
if (!(Environment::getInstance().hasModule(srcName)))
|
||||
{
|
||||
MSource::Point::Par pointPar;
|
||||
pointPar.position = pos;
|
||||
application.createModule<MSource::Point>(srcName, pointPar);
|
||||
}
|
||||
}
|
||||
|
||||
/*******************************************************************************
|
||||
* Name: makeSequentialSource
|
||||
* Purpose: Construct sequential source and add to application module.
|
||||
* Parameters: application - main application that stores modules.
|
||||
* srcName - name of source module to create.
|
||||
* qSrc - Input quark for sequential inversion.
|
||||
* tS - sequential source timeslice.
|
||||
* mom - momentum insertion (default is zero).
|
||||
* Returns: None.
|
||||
******************************************************************************/
|
||||
inline void makeSequentialSource(Application &application, std::string srcName,
|
||||
std::string qSrc, unsigned int tS,
|
||||
std::string mom = ZERO_MOM)
|
||||
{
|
||||
// If the source already exists, don't make the module again.
|
||||
if (!(Environment::getInstance().hasModule(srcName)))
|
||||
{
|
||||
MSource::SeqGamma::Par seqPar;
|
||||
seqPar.q = qSrc;
|
||||
seqPar.tA = tS;
|
||||
seqPar.tB = tS;
|
||||
seqPar.mom = mom;
|
||||
seqPar.gamma = Gamma::Algebra::GammaT;
|
||||
application.createModule<MSource::SeqGamma>(srcName, seqPar);
|
||||
}
|
||||
}
|
||||
|
||||
/*******************************************************************************
|
||||
* Name: makeWallSource
|
||||
* Purpose: Construct wall source and add to application module.
|
||||
* Parameters: application - main application that stores modules.
|
||||
* srcName - name of source module to create.
|
||||
* tW - wall source timeslice.
|
||||
* mom - momentum insertion (default is zero).
|
||||
* Returns: None.
|
||||
******************************************************************************/
|
||||
inline void makeWallSource(Application &application, std::string srcName,
|
||||
unsigned int tW, std::string mom = ZERO_MOM)
|
||||
{
|
||||
// If the source already exists, don't make the module again.
|
||||
if (!(Environment::getInstance().hasModule(srcName)))
|
||||
{
|
||||
MSource::Wall::Par wallPar;
|
||||
wallPar.tW = tW;
|
||||
wallPar.mom = mom;
|
||||
application.createModule<MSource::Wall>(srcName, wallPar);
|
||||
}
|
||||
}
|
||||
|
||||
/*******************************************************************************
|
||||
* Name: makeWallSink
|
||||
* Purpose: Wall sink smearing of a propagator.
|
||||
* Parameters: application - main application that stores modules.
|
||||
* propName - name of input propagator.
|
||||
* wallName - name of smeared propagator.
|
||||
* mom - momentum insertion (default is zero).
|
||||
* Returns: None.
|
||||
******************************************************************************/
|
||||
inline void makeWallSink(Application &application, std::string propName,
|
||||
std::string wallName, std::string mom = ZERO_MOM)
|
||||
{
|
||||
// If the propagator has already been smeared, don't smear it again.
|
||||
// Temporarily removed, strategy for sink smearing likely to change.
|
||||
/*if (!(Environment::getInstance().hasModule(wallName)))
|
||||
{
|
||||
MSink::Wall::Par wallPar;
|
||||
wallPar.q = propName;
|
||||
wallPar.mom = mom;
|
||||
application.createModule<MSink::Wall>(wallName, wallPar);
|
||||
}*/
|
||||
}
|
||||
|
||||
/*******************************************************************************
|
||||
* Name: makePropagator
|
||||
* Purpose: Construct source and propagator then add to application module.
|
||||
* Parameters: application - main application that stores modules.
|
||||
* propName - name of propagator module to create.
|
||||
* srcName - name of source module to use.
|
||||
* solver - solver to use (default is CG).
|
||||
* Returns: None.
|
||||
******************************************************************************/
|
||||
inline void makePropagator(Application &application, std::string &propName,
|
||||
std::string &srcName, std::string &solver)
|
||||
{
|
||||
// If the propagator already exists, don't make the module again.
|
||||
if (!(Environment::getInstance().hasModule(propName)))
|
||||
{
|
||||
Quark::Par quarkPar;
|
||||
quarkPar.source = srcName;
|
||||
quarkPar.solver = solver;
|
||||
application.createModule<Quark>(propName, quarkPar);
|
||||
}
|
||||
}
|
||||
|
||||
/*******************************************************************************
|
||||
* Name: makeLoop
|
||||
* Purpose: Use noise source and inversion result to make loop propagator, then
|
||||
* add to application module.
|
||||
* Parameters: application - main application that stores modules.
|
||||
* propName - name of propagator module to create.
|
||||
* srcName - name of noise source module to use.
|
||||
* resName - name of inversion result on given noise source.
|
||||
* Returns: None.
|
||||
******************************************************************************/
|
||||
inline void makeLoop(Application &application, std::string &propName,
|
||||
std::string &srcName, std::string &resName)
|
||||
{
|
||||
// If the loop propagator already exists, don't make the module again.
|
||||
if (!(Environment::getInstance().hasModule(propName)))
|
||||
{
|
||||
MLoop::NoiseLoop::Par loopPar;
|
||||
loopPar.q = resName;
|
||||
loopPar.eta = srcName;
|
||||
application.createModule<MLoop::NoiseLoop>(propName, loopPar);
|
||||
}
|
||||
}
|
||||
|
||||
/*******************************************************************************
|
||||
* Contraction module creation.
|
||||
******************************************************************************/
|
||||
|
||||
/*******************************************************************************
|
||||
* Name: mesonContraction
|
||||
* Purpose: Create meson contraction module and add to application module.
|
||||
* Parameters: application - main application that stores modules.
|
||||
* npt - specify n-point correlator (for labelling).
|
||||
* q1 - quark propagator 1.
|
||||
* q2 - quark propagator 2.
|
||||
* label - unique label to construct module name.
|
||||
* mom - momentum to project (default is zero)
|
||||
* gammas - gamma insertions at source and sink.
|
||||
* Returns: None.
|
||||
******************************************************************************/
|
||||
inline void mesonContraction(Application &application, unsigned int npt,
|
||||
std::string &q1, std::string &q2,
|
||||
std::string &label,
|
||||
std::string mom = ZERO_MOM,
|
||||
std::string gammas = "<Gamma5 Gamma5>")
|
||||
{
|
||||
std::string modName = std::to_string(npt) + "pt_" + label;
|
||||
if (!(Environment::getInstance().hasModule(modName)))
|
||||
{
|
||||
MContraction::Meson::Par mesPar;
|
||||
mesPar.output = std::to_string(npt) + "pt/" + label;
|
||||
mesPar.q1 = q1;
|
||||
mesPar.q2 = q2;
|
||||
mesPar.mom = mom;
|
||||
mesPar.gammas = gammas;
|
||||
application.createModule<MContraction::Meson>(modName, mesPar);
|
||||
}
|
||||
}
|
||||
|
||||
/*******************************************************************************
|
||||
* Name: gamma3ptContraction
|
||||
* Purpose: Create gamma3pt contraction module and add to application module.
|
||||
* Parameters: application - main application that stores modules.
|
||||
* npt - specify n-point correlator (for labelling).
|
||||
* q1 - quark propagator 1.
|
||||
* q2 - quark propagator 2.
|
||||
* q3 - quark propagator 3.
|
||||
* label - unique label to construct module name.
|
||||
* gamma - gamma insertions between q2 and q3.
|
||||
* Returns: None.
|
||||
******************************************************************************/
|
||||
inline void gamma3ptContraction(Application &application, unsigned int npt,
|
||||
std::string &q1, std::string &q2,
|
||||
std::string &q3, std::string &label,
|
||||
Gamma::Algebra gamma = Gamma::Algebra::Identity)
|
||||
{
|
||||
std::string modName = std::to_string(npt) + "pt_" + label;
|
||||
if (!(Environment::getInstance().hasModule(modName)))
|
||||
{
|
||||
MContraction::Gamma3pt::Par gamma3ptPar;
|
||||
gamma3ptPar.output = std::to_string(npt) + "pt/" + label;
|
||||
gamma3ptPar.q1 = q1;
|
||||
gamma3ptPar.q2 = q2;
|
||||
gamma3ptPar.q3 = q3;
|
||||
gamma3ptPar.gamma = gamma;
|
||||
application.createModule<MContraction::Gamma3pt>(modName, gamma3ptPar);
|
||||
}
|
||||
}
|
||||
|
||||
/*******************************************************************************
|
||||
* Name: weakContraction[Eye,NonEye]
|
||||
* Purpose: Create Weak Hamiltonian contraction module for Eye/NonEye topology
|
||||
* and add to application module.
|
||||
* Parameters: application - main application that stores modules.
|
||||
* npt - specify n-point correlator (for labelling).
|
||||
* q1 - quark propagator 1.
|
||||
* q2 - quark propagator 2.
|
||||
* q3 - quark propagator 3.
|
||||
* q4 - quark propagator 4.
|
||||
* label - unique label to construct module name.
|
||||
* Returns: None.
|
||||
******************************************************************************/
|
||||
#define HW_CONTRACTION(top) \
|
||||
inline void weakContraction##top(Application &application, unsigned int npt,\
|
||||
std::string &q1, std::string &q2, \
|
||||
std::string &q3, std::string &q4, \
|
||||
std::string &label)\
|
||||
{\
|
||||
std::string modName = std::to_string(npt) + "pt_" + label;\
|
||||
if (!(Environment::getInstance().hasModule(modName)))\
|
||||
{\
|
||||
MContraction::WeakHamiltonian##top::Par weakPar;\
|
||||
weakPar.output = std::to_string(npt) + "pt/" + label;\
|
||||
weakPar.q1 = q1;\
|
||||
weakPar.q2 = q2;\
|
||||
weakPar.q3 = q3;\
|
||||
weakPar.q4 = q4;\
|
||||
application.createModule<MContraction::WeakHamiltonian##top>(modName, weakPar);\
|
||||
}\
|
||||
}
|
||||
HW_CONTRACTION(Eye) // weakContractionEye
|
||||
HW_CONTRACTION(NonEye) // weakContractionNonEye
|
||||
|
||||
/*******************************************************************************
|
||||
* Name: disc0Contraction
|
||||
* Purpose: Create contraction module for 4pt Weak Hamiltonian + current
|
||||
* disconnected topology for neutral mesons and add to application
|
||||
* module.
|
||||
* Parameters: application - main application that stores modules.
|
||||
* q1 - quark propagator 1.
|
||||
* q2 - quark propagator 2.
|
||||
* q3 - quark propagator 3.
|
||||
* q4 - quark propagator 4.
|
||||
* label - unique label to construct module name.
|
||||
* Returns: None.
|
||||
******************************************************************************/
|
||||
inline void disc0Contraction(Application &application,
|
||||
std::string &q1, std::string &q2,
|
||||
std::string &q3, std::string &q4,
|
||||
std::string &label)
|
||||
{
|
||||
std::string modName = "4pt_" + label;
|
||||
if (!(Environment::getInstance().hasModule(modName)))
|
||||
{
|
||||
MContraction::WeakNeutral4ptDisc::Par disc0Par;
|
||||
disc0Par.output = "4pt/" + label;
|
||||
disc0Par.q1 = q1;
|
||||
disc0Par.q2 = q2;
|
||||
disc0Par.q3 = q3;
|
||||
disc0Par.q4 = q4;
|
||||
application.createModule<MContraction::WeakNeutral4ptDisc>(modName, disc0Par);
|
||||
}
|
||||
}
|
||||
|
||||
/*******************************************************************************
|
||||
* Name: discLoopContraction
|
||||
* Purpose: Create contraction module for disconnected loop and add to
|
||||
* application module.
|
||||
* Parameters: application - main application that stores modules.
|
||||
* q_loop - loop quark propagator.
|
||||
* modName - unique module name.
|
||||
* gamma - gamma matrix to use in contraction.
|
||||
* Returns: None.
|
||||
******************************************************************************/
|
||||
inline void discLoopContraction(Application &application,
|
||||
std::string &q_loop, std::string &modName,
|
||||
Gamma::Algebra gamma = Gamma::Algebra::Identity)
|
||||
{
|
||||
if (!(Environment::getInstance().hasModule(modName)))
|
||||
{
|
||||
MContraction::DiscLoop::Par discPar;
|
||||
discPar.output = "disc/" + modName;
|
||||
discPar.q_loop = q_loop;
|
||||
discPar.gamma = gamma;
|
||||
application.createModule<MContraction::DiscLoop>(modName, discPar);
|
||||
}
|
||||
}
|
@ -30,14 +30,6 @@
|
||||
using namespace Grid;
|
||||
using namespace Hadrons;
|
||||
|
||||
static Gamma::Algebra gmu[4] =
|
||||
{
|
||||
Gamma::Algebra::GammaX,
|
||||
Gamma::Algebra::GammaY,
|
||||
Gamma::Algebra::GammaZ,
|
||||
Gamma::Algebra::GammaT
|
||||
};
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
// initialization //////////////////////////////////////////////////////////
|
||||
@ -110,7 +102,7 @@ int main(int argc, char *argv[])
|
||||
seqName.push_back(std::vector<std::string>(Nd));
|
||||
for (unsigned int mu = 0; mu < Nd; ++mu)
|
||||
{
|
||||
seqPar.gamma = gmu[mu];
|
||||
seqPar.gamma = 0x1 << mu;
|
||||
seqName[i][mu] = "G" + std::to_string(seqPar.gamma)
|
||||
+ "_" + std::to_string(seqPar.tA) + "-"
|
||||
+ qName[i];
|
||||
@ -135,11 +127,11 @@ int main(int argc, char *argv[])
|
||||
for (unsigned int i = 0; i < flavour.size(); ++i)
|
||||
for (unsigned int j = i; j < flavour.size(); ++j)
|
||||
{
|
||||
mesPar.output = "mesons/Z2_" + flavour[i] + flavour[j];
|
||||
mesPar.q1 = qName[i];
|
||||
mesPar.q2 = qName[j];
|
||||
mesPar.gammaSource = Gamma::Algebra::Gamma5;
|
||||
mesPar.gammaSink = Gamma::Algebra::Gamma5;
|
||||
mesPar.output = "mesons/Z2_" + flavour[i] + flavour[j];
|
||||
mesPar.q1 = qName[i];
|
||||
mesPar.q2 = qName[j];
|
||||
mesPar.gammas = "all";
|
||||
mesPar.mom = "0. 0. 0. 0.";
|
||||
application.createModule<MContraction::Meson>("meson_Z2_"
|
||||
+ std::to_string(t)
|
||||
+ "_"
|
||||
@ -157,6 +149,8 @@ int main(int argc, char *argv[])
|
||||
+ std::to_string(mu);
|
||||
mesPar.q1 = qName[i];
|
||||
mesPar.q2 = seqName[j][mu];
|
||||
mesPar.gammas = "all";
|
||||
mesPar.mom = "0. 0. 0. 0.";
|
||||
application.createModule<MContraction::Meson>("3pt_Z2_"
|
||||
+ std::to_string(t)
|
||||
+ "_"
|
||||
|
337
tests/hadrons/Test_hadrons_rarekaon.cc
Normal file
337
tests/hadrons/Test_hadrons_rarekaon.cc
Normal file
@ -0,0 +1,337 @@
|
||||
/*******************************************************************************
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: tests/hadrons/Test_hadrons_rarekaon.cc
|
||||
|
||||
Copyright (C) 2017
|
||||
|
||||
Author: Andrew Lawson <andrew.lawson1991@gmail.com>
|
||||
|
||||
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 "Test_hadrons.hpp"
|
||||
|
||||
using namespace Grid;
|
||||
using namespace Hadrons;
|
||||
|
||||
enum quarks
|
||||
{
|
||||
light = 0,
|
||||
strange = 1,
|
||||
charm = 2
|
||||
};
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
// parse command line //////////////////////////////////////////////////////
|
||||
std::string configStem;
|
||||
|
||||
if (argc < 2)
|
||||
{
|
||||
std::cerr << "usage: " << argv[0] << " <configuration filestem> [Grid options]";
|
||||
std::cerr << std::endl;
|
||||
std::exit(EXIT_FAILURE);
|
||||
}
|
||||
configStem = argv[1];
|
||||
|
||||
// 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<double> mass = {.01, .04, .2};
|
||||
std::vector<std::string> flavour = {"l", "s", "c"};
|
||||
std::vector<std::string> solvers = {"CG_l", "CG_s", "CG_c"};
|
||||
std::string kmom = "0. 0. 0. 0.";
|
||||
std::string pmom = "1. 0. 0. 0.";
|
||||
std::string qmom = "-1. 0. 0. 0.";
|
||||
std::string mqmom = "1. 0. 0. 0.";
|
||||
std::vector<unsigned int> tKs = {0};
|
||||
unsigned int dt_pi = 16;
|
||||
std::vector<unsigned int> tJs = {8};
|
||||
unsigned int n_noise = 1;
|
||||
unsigned int nt = 32;
|
||||
bool do_disconnected(false);
|
||||
|
||||
// Global parameters.
|
||||
Application::GlobalPar globalPar;
|
||||
globalPar.trajCounter.start = 1500;
|
||||
globalPar.trajCounter.end = 1520;
|
||||
globalPar.trajCounter.step = 20;
|
||||
globalPar.seed = "1 2 3 4";
|
||||
globalPar.genetic.maxGen = 1000;
|
||||
globalPar.genetic.maxCstGen = 200;
|
||||
globalPar.genetic.popSize = 20;
|
||||
globalPar.genetic.mutationRate = .1;
|
||||
application.setPar(globalPar);
|
||||
|
||||
// gauge field
|
||||
if (configStem == "None")
|
||||
{
|
||||
application.createModule<MGauge::Unit>("gauge");
|
||||
}
|
||||
else
|
||||
{
|
||||
MGauge::Load::Par gaugePar;
|
||||
gaugePar.file = configStem;
|
||||
application.createModule<MGauge::Load>("gauge", gaugePar);
|
||||
}
|
||||
for (unsigned int i = 0; i < flavour.size(); ++i)
|
||||
{
|
||||
// actions
|
||||
MAction::DWF::Par actionPar;
|
||||
actionPar.gauge = "gauge";
|
||||
actionPar.Ls = 16;
|
||||
actionPar.M5 = 1.8;
|
||||
actionPar.mass = mass[i];
|
||||
application.createModule<MAction::DWF>("DWF_" + flavour[i], actionPar);
|
||||
|
||||
// solvers
|
||||
// RBPrecCG -> CG
|
||||
MSolver::RBPrecCG::Par solverPar;
|
||||
solverPar.action = "DWF_" + flavour[i];
|
||||
solverPar.residual = 1.0e-8;
|
||||
application.createModule<MSolver::RBPrecCG>(solvers[i],
|
||||
solverPar);
|
||||
}
|
||||
|
||||
// Create noise propagators for loops.
|
||||
std::vector<std::string> noiseSrcs;
|
||||
std::vector<std::vector<std::string>> noiseRes;
|
||||
std::vector<std::vector<std::string>> noiseProps;
|
||||
if (n_noise > 0)
|
||||
{
|
||||
MSource::Z2::Par noisePar;
|
||||
noisePar.tA = 0;
|
||||
noisePar.tB = nt - 1;
|
||||
std::string loop_stem = "loop_";
|
||||
|
||||
noiseRes.resize(flavour.size());
|
||||
noiseProps.resize(flavour.size());
|
||||
for (unsigned int nn = 0; nn < n_noise; ++nn)
|
||||
{
|
||||
std::string eta = INIT_INDEX("noise", nn);
|
||||
application.createModule<MSource::Z2>(eta, noisePar);
|
||||
noiseSrcs.push_back(eta);
|
||||
|
||||
for (unsigned int f = 0; f < flavour.size(); ++f)
|
||||
{
|
||||
std::string loop_prop = INIT_INDEX(loop_stem + flavour[f], nn);
|
||||
std::string loop_res = loop_prop + "_res";
|
||||
makePropagator(application, loop_res, eta, solvers[f]);
|
||||
makeLoop(application, loop_prop, eta, loop_res);
|
||||
noiseRes[f].push_back(loop_res);
|
||||
noiseProps[f].push_back(loop_prop);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Translate rare kaon decay across specified timeslices.
|
||||
for (unsigned int i = 0; i < tKs.size(); ++i)
|
||||
{
|
||||
// Zero-momentum wall source propagators for kaon and pion.
|
||||
unsigned int tK = tKs[i];
|
||||
unsigned int tpi = (tK + dt_pi) % nt;
|
||||
std::string q_Kl_0 = INIT_INDEX("Q_l_0", tK);
|
||||
std::string q_pil_0 = INIT_INDEX("Q_l_0", tpi);
|
||||
MAKE_WALL_PROP(tK, q_Kl_0, solvers[light]);
|
||||
MAKE_WALL_PROP(tpi, q_pil_0, solvers[light]);
|
||||
|
||||
// Wall sources for kaon and pion with momentum insertion. If either
|
||||
// p or k are zero, or p = k, re-use the existing name to avoid
|
||||
// duplicating a propagator.
|
||||
std::string q_Ks_k = INIT_INDEX("Q_Ks_k", tK);
|
||||
std::string q_Ks_p = INIT_INDEX((kmom == pmom) ? "Q_Ks_k" : "Q_Ks_p", tK);
|
||||
std::string q_pil_k = INIT_INDEX((kmom == ZERO_MOM) ? "Q_l_0" : "Q_l_k", tpi);
|
||||
std::string q_pil_p = INIT_INDEX((pmom == kmom) ? q_pil_k : ((pmom == ZERO_MOM) ? "Q_l_0" : "Q_l_p"), tpi);
|
||||
MAKE_3MOM_WALL_PROP(tK, kmom, q_Ks_k, solvers[strange]);
|
||||
MAKE_3MOM_WALL_PROP(tK, pmom, q_Ks_p, solvers[strange]);
|
||||
MAKE_3MOM_WALL_PROP(tpi, kmom, q_pil_k, solvers[light]);
|
||||
MAKE_3MOM_WALL_PROP(tpi, pmom, q_pil_p, solvers[light]);
|
||||
|
||||
/***********************************************************************
|
||||
* CONTRACTIONS: pi and K 2pt contractions with mom = p, k.
|
||||
**********************************************************************/
|
||||
// Wall-Point
|
||||
std::string PW_K_k = INIT_INDEX("PW_K_k", tK);
|
||||
std::string PW_K_p = INIT_INDEX("PW_K_p", tK);
|
||||
std::string PW_pi_k = INIT_INDEX("PW_pi_k", tpi);
|
||||
std::string PW_pi_p = INIT_INDEX("PW_pi_p", tpi);
|
||||
mesonContraction(application, 2, q_Kl_0, q_Ks_k, PW_K_k, kmom);
|
||||
mesonContraction(application, 2, q_Kl_0, q_Ks_p, PW_K_p, pmom);
|
||||
mesonContraction(application, 2, q_pil_k, q_pil_0, PW_pi_k, kmom);
|
||||
mesonContraction(application, 2, q_pil_p, q_pil_0, PW_pi_p, pmom);
|
||||
// Wall-Wall, to be done - requires modification of meson module.
|
||||
|
||||
/***********************************************************************
|
||||
* CONTRACTIONS: 3pt Weak Hamiltonian, C & W (non-Eye type) classes.
|
||||
**********************************************************************/
|
||||
std::string HW_CW_k = LABEL_3PT("HW_CW_k", tK, tpi);
|
||||
std::string HW_CW_p = LABEL_3PT("HW_CW_p", tK, tpi);
|
||||
weakContractionNonEye(application, 3, q_Kl_0, q_Ks_k, q_pil_k, q_pil_0, HW_CW_k);
|
||||
weakContractionNonEye(application, 3, q_Kl_0, q_Ks_p, q_pil_p, q_pil_0, HW_CW_p);
|
||||
|
||||
/***********************************************************************
|
||||
* CONTRACTIONS: 3pt sd insertion.
|
||||
**********************************************************************/
|
||||
// Note: eventually will use wall sink smeared q_Kl_0 instead.
|
||||
std::string sd_k = LABEL_3PT("sd_k", tK, tpi);
|
||||
std::string sd_p = LABEL_3PT("sd_p", tK, tpi);
|
||||
gamma3ptContraction(application, 3, q_Kl_0, q_Ks_k, q_pil_k, sd_k);
|
||||
gamma3ptContraction(application, 3, q_Kl_0, q_Ks_p, q_pil_p, sd_p);
|
||||
|
||||
for (unsigned int nn = 0; nn < n_noise; ++nn)
|
||||
{
|
||||
/*******************************************************************
|
||||
* CONTRACTIONS: 3pt Weak Hamiltonian, S and E (Eye type) classes.
|
||||
******************************************************************/
|
||||
// Note: eventually will use wall sink smeared q_Kl_0 instead.
|
||||
for (unsigned int f = 0; f < flavour.size(); ++f)
|
||||
{
|
||||
if ((f != strange) || do_disconnected)
|
||||
{
|
||||
std::string HW_SE_k = LABEL_3PT("HW_SE_k_" + flavour[f], tK, tpi);
|
||||
std::string HW_SE_p = LABEL_3PT("HW_SE_p_" + flavour[f], tK, tpi);
|
||||
std::string loop_q = noiseProps[f][nn];
|
||||
weakContractionEye(application, 3, q_Kl_0, q_Ks_k, q_pil_k, loop_q, HW_CW_k);
|
||||
weakContractionEye(application, 3, q_Kl_0, q_Ks_p, q_pil_p, loop_q, HW_CW_p);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Perform separate contractions for each t_J position.
|
||||
for (unsigned int j = 0; j < tJs.size(); ++j)
|
||||
{
|
||||
// Sequential sources for current insertions. Local for now,
|
||||
// gamma_0 only.
|
||||
unsigned int tJ = (tJs[j] + tK) % nt;
|
||||
MSource::SeqGamma::Par seqPar;
|
||||
std::string q_KlCl_q = LABEL_3PT("Q_KlCl_q", tK, tJ);
|
||||
std::string q_KsCs_mq = LABEL_3PT("Q_KsCs_mq", tK, tJ);
|
||||
std::string q_pilCl_q = LABEL_3PT("Q_pilCl_q", tpi, tJ);
|
||||
std::string q_pilCl_mq = LABEL_3PT("Q_pilCl_mq", tpi, tJ);
|
||||
MAKE_SEQUENTIAL_PROP(tJ, q_Kl_0, qmom, q_KlCl_q, solvers[light]);
|
||||
MAKE_SEQUENTIAL_PROP(tJ, q_Ks_k, mqmom, q_KsCs_mq, solvers[strange]);
|
||||
MAKE_SEQUENTIAL_PROP(tJ, q_pil_p, qmom, q_pilCl_q, solvers[light]);
|
||||
MAKE_SEQUENTIAL_PROP(tJ, q_pil_0, mqmom, q_pilCl_mq, solvers[light]);
|
||||
|
||||
/*******************************************************************
|
||||
* CONTRACTIONS: pi and K 3pt contractions with current insertion.
|
||||
******************************************************************/
|
||||
// Wall-Point
|
||||
std::string C_PW_Kl = LABEL_3PT("C_PW_Kl", tK, tJ);
|
||||
std::string C_PW_Ksb = LABEL_3PT("C_PW_Ksb", tK, tJ);
|
||||
std::string C_PW_pilb = LABEL_3PT("C_PW_pilb", tK, tJ);
|
||||
std::string C_PW_pil = LABEL_3PT("C_PW_pil", tK, tJ);
|
||||
mesonContraction(application, 3, q_KlCl_q, q_Ks_k, C_PW_Kl, pmom);
|
||||
mesonContraction(application, 3, q_Kl_0, q_KsCs_mq, C_PW_Ksb, pmom);
|
||||
mesonContraction(application, 3, q_pil_0, q_pilCl_q, C_PW_pilb, kmom);
|
||||
mesonContraction(application, 3, q_pilCl_mq, q_pil_p, C_PW_pil, kmom);
|
||||
// Wall-Wall, to be done.
|
||||
|
||||
/*******************************************************************
|
||||
* CONTRACTIONS: 4pt contractions, C & W classes.
|
||||
******************************************************************/
|
||||
std::string CW_Kl = LABEL_4PT("CW_Kl", tK, tJ, tpi);
|
||||
std::string CW_Ksb = LABEL_4PT("CW_Ksb", tK, tJ, tpi);
|
||||
std::string CW_pilb = LABEL_4PT("CW_pilb", tK, tJ, tpi);
|
||||
std::string CW_pil = LABEL_4PT("CW_pil", tK, tJ, tpi);
|
||||
weakContractionNonEye(application, 4, q_KlCl_q, q_Ks_k, q_pil_p, q_pil_0, CW_Kl);
|
||||
weakContractionNonEye(application, 4, q_Kl_0, q_KsCs_mq, q_pil_p, q_pil_0, CW_Ksb);
|
||||
weakContractionNonEye(application, 4, q_Kl_0, q_Ks_k, q_pilCl_q, q_pil_0, CW_pilb);
|
||||
weakContractionNonEye(application, 4, q_Kl_0, q_Ks_k, q_pil_p, q_pilCl_mq, CW_pil);
|
||||
|
||||
/*******************************************************************
|
||||
* CONTRACTIONS: 4pt contractions, sd insertions.
|
||||
******************************************************************/
|
||||
// Note: eventually will use wall sink smeared q_Kl_0/q_KlCl_q instead.
|
||||
std::string sd_Kl = LABEL_4PT("sd_Kl", tK, tJ, tpi);
|
||||
std::string sd_Ksb = LABEL_4PT("sd_Ksb", tK, tJ, tpi);
|
||||
std::string sd_pilb = LABEL_4PT("sd_pilb", tK, tJ, tpi);
|
||||
gamma3ptContraction(application, 4, q_KlCl_q, q_Ks_k, q_pil_p, sd_Kl);
|
||||
gamma3ptContraction(application, 4, q_Kl_0, q_KsCs_mq, q_pil_p, sd_Ksb);
|
||||
gamma3ptContraction(application, 4, q_Kl_0, q_Ks_k, q_pilCl_q, sd_pilb);
|
||||
|
||||
// Sequential sources for each noise propagator.
|
||||
for (unsigned int nn = 0; nn < n_noise; ++nn)
|
||||
{
|
||||
std::string loop_stem = "loop_";
|
||||
|
||||
// Contraction required for each quark flavour - alternatively
|
||||
// drop the strange loop if not performing disconnected
|
||||
// contractions or neglecting H_W operators Q_3 -> Q_10.
|
||||
for (unsigned int f = 0; f < flavour.size(); ++f)
|
||||
{
|
||||
if ((f != strange) || do_disconnected)
|
||||
{
|
||||
std::string eta = noiseSrcs[nn];
|
||||
std::string loop_q = noiseProps[f][nn];
|
||||
std::string loop_qCq = LABEL_3PT(loop_stem + flavour[f], tJ, nn);
|
||||
std::string loop_qCq_res = loop_qCq + "_res";
|
||||
MAKE_SEQUENTIAL_PROP(tJ, noiseRes[f][nn], qmom,
|
||||
loop_qCq_res, solvers[f]);
|
||||
makeLoop(application, loop_qCq, eta, loop_qCq_res);
|
||||
|
||||
/*******************************************************
|
||||
* CONTRACTIONS: 4pt contractions, S & E classes.
|
||||
******************************************************/
|
||||
// Note: eventually will use wall sink smeared q_Kl_0/q_KlCl_q instead.
|
||||
std::string SE_Kl = LABEL_4PT_NOISE("SE_Kl", tK, tJ, tpi, nn);
|
||||
std::string SE_Ksb = LABEL_4PT_NOISE("SE_Ksb", tK, tJ, tpi, nn);
|
||||
std::string SE_pilb = LABEL_4PT_NOISE("SE_pilb", tK, tJ, tpi, nn);
|
||||
std::string SE_loop = LABEL_4PT_NOISE("SE_loop", tK, tJ, tpi, nn);
|
||||
weakContractionEye(application, 4, q_KlCl_q, q_Ks_k, q_pil_p, loop_q, SE_Kl);
|
||||
weakContractionEye(application, 4, q_Kl_0, q_KsCs_mq, q_pil_p, loop_q, SE_Ksb);
|
||||
weakContractionEye(application, 4, q_Kl_0, q_Ks_k, q_pilCl_q, loop_q, SE_pilb);
|
||||
weakContractionEye(application, 4, q_Kl_0, q_Ks_k, q_pil_p, loop_qCq, SE_loop);
|
||||
|
||||
/*******************************************************
|
||||
* CONTRACTIONS: 4pt contractions, pi0 disconnected
|
||||
* loop.
|
||||
******************************************************/
|
||||
std::string disc0 = LABEL_4PT_NOISE("disc0", tK, tJ, tpi, nn);
|
||||
disc0Contraction(application, q_Kl_0, q_Ks_k, q_pilCl_q, loop_q, disc0);
|
||||
|
||||
/*******************************************************
|
||||
* CONTRACTIONS: Disconnected loop.
|
||||
******************************************************/
|
||||
std::string discLoop = "disc_" + loop_qCq;
|
||||
discLoopContraction(application, loop_qCq, discLoop);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
// execution
|
||||
std::string par_file_name = "rarekaon_000_100_tK0_tpi16_tJ8_noloop_mc0.2.xml";
|
||||
application.saveParameterFile(par_file_name);
|
||||
application.run();
|
||||
|
||||
// epilogue
|
||||
LOG(Message) << "Grid is finalizing now" << std::endl;
|
||||
Grid_finalize();
|
||||
|
||||
return EXIT_SUCCESS;
|
||||
}
|
@ -96,12 +96,16 @@ int main(int argc, char *argv[])
|
||||
mesPar.output = "mesons/pt_" + flavour[i] + flavour[j];
|
||||
mesPar.q1 = "Qpt_" + flavour[i];
|
||||
mesPar.q2 = "Qpt_" + flavour[j];
|
||||
mesPar.gammas = "all";
|
||||
mesPar.mom = "0. 0. 0. 0.";
|
||||
application.createModule<MContraction::Meson>("meson_pt_"
|
||||
+ flavour[i] + flavour[j],
|
||||
mesPar);
|
||||
mesPar.output = "mesons/Z2_" + flavour[i] + flavour[j];
|
||||
mesPar.q1 = "QZ2_" + flavour[i];
|
||||
mesPar.q2 = "QZ2_" + flavour[j];
|
||||
mesPar.gammas = "all";
|
||||
mesPar.mom = "0. 0. 0. 0.";
|
||||
application.createModule<MContraction::Meson>("meson_Z2_"
|
||||
+ flavour[i] + flavour[j],
|
||||
mesPar);
|
||||
|
Loading…
x
Reference in New Issue
Block a user