portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-09 23:45:36 +00:00
Code Releases Activity

Merge pull request #214 from guelpers/feature/kl2QEDseq

Browse Source 
Kl2 contraction with sequential propagators
This commit is contained in:
Antonin Portelli 2019-05-23 20:31:41 +01:00 committed by GitHub
commit c7c0a1065f
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
11 changed files with 848 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
Grid/qcd/action/fermion/DomainWallFermion.h
View File

@ -43,7 +43,7 @@ namespace Grid {
INHERIT_IMPL_TYPES(Impl); INHERIT_IMPL_TYPES(Impl);
public: public:
void FreePropagator(const FermionField &in,FermionField &out,RealD mass, std::vector<double> twist, bool fiveD) { void FreePropagator(const FermionField &in,FermionField &out,RealD mass,std::vector<Complex> boundary, std::vector<double> twist, bool fiveD) {
FermionField in_k(in._grid); FermionField in_k(in._grid);
FermionField prop_k(in._grid); FermionField prop_k(in._grid);
@ -55,15 +55,20 @@ namespace Grid {
FermionField in_buf(in._grid); in_buf = zero; FermionField in_buf(in._grid); in_buf = zero;
Scalar ci(0.0,1.0); Scalar ci(0.0,1.0);
assert(twist.size() == Nd);//check that twist is Nd assert(twist.size() == Nd);//check that twist is Nd
assert(boundary.size() == Nd);//check that boundary conditions is Nd
int shift = 0; int shift = 0;
if(fiveD) shift = 1; if(fiveD) shift = 1;
for(unsigned int nu = 0; nu < Nd; nu++) for(unsigned int nu = 0; nu < Nd; nu++)
{ {
// Shift coordinate lattice index by 1 to account for 5th dimension. // Shift coordinate lattice index by 1 to account for 5th dimension.
LatticeCoordinate(coor, nu + shift); LatticeCoordinate(coor, nu + shift);
ph = ph + twist[nu]*coor*((1./(in._grid->_fdimensions[nu+shift]))); double boundary_phase = ::acos(real(boundary[nu]));
ph = ph + boundary_phase*coor*((1./(in._grid->_fdimensions[nu+shift])));
//momenta for propagator shifted by twist+boundary
twist[nu] = twist[nu] + boundary_phase/((2.0*M_PI));
} }
in_buf = exp(Scalar(2.0*M_PI)*ci*ph*(-1.0))*in; in_buf = exp(ci*ph*(-1.0))*in;
if(fiveD){//FFT only on temporal and spatial dimensions if(fiveD){//FFT only on temporal and spatial dimensions
std::vector<int> mask(Nd+1,1); mask[0] = 0; std::vector<int> mask(Nd+1,1); mask[0] = 0;
@ -76,25 +81,28 @@ namespace Grid {
this->MomentumSpacePropagatorHt(prop_k,in_k,mass,twist); this->MomentumSpacePropagatorHt(prop_k,in_k,mass,twist);
theFFT.FFT_all_dim(out,prop_k,FFT::backward); theFFT.FFT_all_dim(out,prop_k,FFT::backward);
} }
//phase for boundary condition //phase for boundary condition
out = out * exp(Scalar(2.0*M_PI)*ci*ph); out = out * exp(ci*ph);
}; };
virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass,std::vector<double> twist) { virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass,std::vector<Complex> boundary,std::vector<double> twist) {
bool fiveD = true; //5d propagator by default bool fiveD = true; //5d propagator by default
FreePropagator(in,out,mass,twist,fiveD); FreePropagator(in,out,mass,boundary,twist,fiveD);
}; };
virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass, bool fiveD) { virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass, bool fiveD) {
std::vector<double> twist(Nd,0.0); //default: periodic boundarys in all directions std::vector<double> twist(Nd,0.0); //default: periodic boundarys in all directions
FreePropagator(in,out,mass,twist,fiveD); std::vector<Complex> boundary;
for(int i=0;i<Nd;i++) boundary.push_back(1);//default: periodic boundary conditions
FreePropagator(in,out,mass,boundary,twist,fiveD);
}; };
virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass) { virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass) {
bool fiveD = true; //5d propagator by default bool fiveD = true; //5d propagator by default
std::vector<double> twist(Nd,0.0); //default: periodic boundarys in all directions std::vector<double> twist(Nd,0.0); //default: twist angle 0
FreePropagator(in,out,mass,twist,fiveD); std::vector<Complex> boundary;
for(int i=0;i<Nd;i++) boundary.push_back(1); //default: periodic boundary conditions
FreePropagator(in,out,mass,boundary,twist,fiveD);
}; };
virtual void Instantiatable(void) {}; virtual void Instantiatable(void) {};

17
Grid/qcd/action/fermion/FermionOperator.h
View File

@ -96,7 +96,7 @@ namespace Grid {
virtual void MomentumSpacePropagator(FermionField &out,const FermionField &in,RealD _m,std::vector<double> twist) { assert(0);}; virtual void MomentumSpacePropagator(FermionField &out,const FermionField &in,RealD _m,std::vector<double> twist) { assert(0);};
virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass,std::vector<double> twist) { virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass,std::vector<Complex> boundary,std::vector<double> twist) {
FFT theFFT((GridCartesian *) in._grid); FFT theFFT((GridCartesian *) in._grid);
FermionField in_k(in._grid); FermionField in_k(in._grid);
@ -108,24 +108,31 @@ namespace Grid {
FermionField in_buf(in._grid); in_buf = zero; FermionField in_buf(in._grid); in_buf = zero;
Scalar ci(0.0,1.0); Scalar ci(0.0,1.0);
assert(twist.size() == Nd);//check that twist is Nd assert(twist.size() == Nd);//check that twist is Nd
assert(boundary.size() == Nd);//check that boundary conditions is Nd
for(unsigned int nu = 0; nu < Nd; nu++) for(unsigned int nu = 0; nu < Nd; nu++)
{ {
LatticeCoordinate(coor, nu); LatticeCoordinate(coor, nu);
ph = ph + twist[nu]*coor*((1./(in._grid->_fdimensions[nu]))); double boundary_phase = ::acos(real(boundary[nu]));
ph = ph + boundary_phase*coor*((1./(in._grid->_fdimensions[nu])));
//momenta for propagator shifted by twist+boundary
twist[nu] = twist[nu] + boundary_phase/((2.0*M_PI));
} }
in_buf = exp(Scalar(-2.0*M_PI)*ci*ph)*in; in_buf = exp(ci*ph*(-1.0))*in;
theFFT.FFT_all_dim(in_k,in_buf,FFT::forward); theFFT.FFT_all_dim(in_k,in_buf,FFT::forward);
this->MomentumSpacePropagator(prop_k,in_k,mass,twist); this->MomentumSpacePropagator(prop_k,in_k,mass,twist);
theFFT.FFT_all_dim(out,prop_k,FFT::backward); theFFT.FFT_all_dim(out,prop_k,FFT::backward);
//phase for boundary condition //phase for boundary condition
out = out * exp(Scalar(2.0*M_PI)*ci*ph); out = out * exp(ci*ph);
}; };
virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass) { virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass) {
std::vector<Complex> boundary;
for(int i=0;i<Nd;i++) boundary.push_back(1);//default: periodic boundary conditions
std::vector<double> twist(Nd,0.0); //default: periodic boundarys in all directions std::vector<double> twist(Nd,0.0); //default: periodic boundarys in all directions
FreePropagator(in,out,mass,twist); FreePropagator(in,out,mass,boundary,twist);
}; };
/////////////////////////////////////////////// ///////////////////////////////////////////////

3
Hadrons/Modules.hpp
View File

@ -6,10 +6,13 @@
#include <Hadrons/Modules/MContraction/Meson.hpp> #include <Hadrons/Modules/MContraction/Meson.hpp>
#include <Hadrons/Modules/MContraction/DiscLoop.hpp> #include <Hadrons/Modules/MContraction/DiscLoop.hpp>
#include <Hadrons/Modules/MContraction/Gamma3pt.hpp> #include <Hadrons/Modules/MContraction/Gamma3pt.hpp>
#include <Hadrons/Modules/MContraction/WeakMesonDecayKl2.hpp>
#include <Hadrons/Modules/MContraction/WeakNonEye3pt.hpp> #include <Hadrons/Modules/MContraction/WeakNonEye3pt.hpp>
#include <Hadrons/Modules/MFermion/FreeProp.hpp> #include <Hadrons/Modules/MFermion/FreeProp.hpp>
#include <Hadrons/Modules/MFermion/GaugeProp.hpp> #include <Hadrons/Modules/MFermion/GaugeProp.hpp>
#include <Hadrons/Modules/MFermion/EMLepton.hpp>
#include <Hadrons/Modules/MSource/SeqGamma.hpp> #include <Hadrons/Modules/MSource/SeqGamma.hpp>
#include <Hadrons/Modules/MSource/SeqAslash.hpp>
#include <Hadrons/Modules/MSource/Point.hpp> #include <Hadrons/Modules/MSource/Point.hpp>
#include <Hadrons/Modules/MSource/Wall.hpp> #include <Hadrons/Modules/MSource/Wall.hpp>
#include <Hadrons/Modules/MSource/Z2.hpp> #include <Hadrons/Modules/MSource/Z2.hpp>

36
Hadrons/Modules/MContraction/WeakMesonDecayKl2.cc Normal file
View File

@ -0,0 +1,36 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Hadrons/Modules/MContraction/WeakMesonDecayKl2.cc
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
Author: Vera Guelpers <Vera.Guelpers@ed.ac.uk>
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 <Hadrons/Modules/MContraction/WeakMesonDecayKl2.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MContraction;
template class Grid::Hadrons::MContraction::TWeakMesonDecayKl2<FIMPL>;

210
Hadrons/Modules/MContraction/WeakMesonDecayKl2.hpp Normal file
View File

@ -0,0 +1,210 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Hadrons/Modules/MContraction/WeakMesonDecayKl2.hpp
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
Author: Vera Guelpers <Vera.Guelpers@ed.ac.uk>
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_MContraction_WeakMesonDecayKl2_hpp_
#define Hadrons_MContraction_WeakMesonDecayKl2_hpp_
#include <Hadrons/Global.hpp>
#include <Hadrons/Module.hpp>
#include <Hadrons/ModuleFactory.hpp>
BEGIN_HADRONS_NAMESPACE
/*
* Kl2 contraction
* -----------------------------
*
* contraction for Kl2 decay, including the lepton
*
* trace(q1*adj(q2)*g5*gL[mu]) * (gL[mu] * lepton)_{a,b}
*
* with open spinor indices (a,b) for the lepton part
*
* q1 lepton
* /------------\ /------------
* / \ /
* / \H_W/
* g_5 * * *
* \ /
* \ /
* \____________/
* q2
*
* * options:
* - q1: input propagator 1 (string)
* - q2: input propagator 2 (string)
* - lepton: input lepton (string)
*/
/******************************************************************************
* TWeakMesonDecayKl2 *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MContraction)
class WeakMesonDecayKl2Par: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(WeakMesonDecayKl2Par,
std::string, q1,
std::string, q2,
std::string, lepton,
std::string, output);
};
template <typename FImpl>
class TWeakMesonDecayKl2: public Module<WeakMesonDecayKl2Par>
{
public:
FERM_TYPE_ALIASES(FImpl,);
class Metadata: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(Metadata,
int, spinidx1,
int, spinidx2);
};
typedef Correlator<Metadata> Result;
public:
// constructor
TWeakMesonDecayKl2(const std::string name);
// destructor
virtual ~TWeakMesonDecayKl2(void) {};
// dependencies/products
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
protected:
// execution
virtual void setup(void);
// execution
virtual void execute(void);
};
MODULE_REGISTER_TMP(WeakMesonDecayKl2, TWeakMesonDecayKl2<FIMPL>, MContraction);
/******************************************************************************
* TWeakMesonDecayKl2 implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TWeakMesonDecayKl2<FImpl>::TWeakMesonDecayKl2(const std::string name)
: Module<WeakMesonDecayKl2Par>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TWeakMesonDecayKl2<FImpl>::getInput(void)
{
std::vector<std::string> input = {par().q1, par().q2, par().lepton};
return input;
}
template <typename FImpl>
std::vector<std::string> TWeakMesonDecayKl2<FImpl>::getOutput(void)
{
std::vector<std::string> output = {};
return output;
}
// setup ////////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TWeakMesonDecayKl2<FImpl>::setup(void)
{
envTmpLat(LatticeComplex, "c");
envTmpLat(PropagatorField, "prop_buf");
envCreateLat(PropagatorField, getName());
envTmpLat(LatticeComplex, "buf");
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TWeakMesonDecayKl2<FImpl>::execute(void)
{
LOG(Message) << "Computing QED Kl2 contractions '" << getName() << "' using"
<< " quarks '" << par().q1 << "' and '" << par().q2 << "' and"
<< "lepton '" << par().lepton << "'" << std::endl;
auto &res = envGet(PropagatorField, getName()); res = zero;
Gamma g5(Gamma::Algebra::Gamma5);
int nt = env().getDim(Tp);
auto &q1 = envGet(PropagatorField, par().q1);
auto &q2 = envGet(PropagatorField, par().q2);
auto &lepton = envGet(PropagatorField, par().lepton);
envGetTmp(LatticeComplex, buf);
std::vector<TComplex> res_summed;
envGetTmp(LatticeComplex, c);
envGetTmp(PropagatorField, prop_buf);
std::vector<Result> result;
Result r;
for (unsigned int mu = 0; mu < 4; ++mu)
{
c = zero;
//hadronic part: trace(q1*adj(q2)*g5*gL[mu])
c = trace(q1*adj(q2)*g5*GammaL(Gamma::gmu[mu]));
prop_buf = 1.;
//multiply lepton part
res += c * prop_buf * GammaL(Gamma::gmu[mu]) * lepton;
}
//loop over spinor index of lepton part
unsigned int i = 0;
for (unsigned int s1 = 0; s1 < Ns ; ++s1)
for (unsigned int s2 = 0; s2 < Ns ; ++s2)
{
buf = peekColour(peekSpin(res,s1,s2),0,0);
sliceSum(buf, res_summed, Tp);
r.corr.clear();
for (unsigned int t = 0; t < nt; ++t)
{
r.corr.push_back(TensorRemove(res_summed[t]));
}
r.info.spinidx1 = s1;
r.info.spinidx2 = s2;
result.push_back(r);
i+=1;
}
saveResult(par().output, "weakdecay", result);
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MContraction_WeakMesonDecayKl2_hpp_

35
Hadrons/Modules/MFermion/EMLepton.cc Normal file
View File

@ -0,0 +1,35 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Hadrons/Modules/MFermion/EMLepton.cc
Copyright (C) 2015-2019
Author: Vera Guelpers <Vera.Guelpers@ed.ac.uk>
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 <Hadrons/Modules/MFermion/EMLepton.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MFermion;
template class Grid::Hadrons::MFermion::TEMLepton<FIMPL>;

295
Hadrons/Modules/MFermion/EMLepton.hpp Normal file
View File

@ -0,0 +1,295 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Hadrons/Modules/MFermion/EMLepton.hpp
Copyright (C) 2015-2019
Author: Vera Guelpers <Vera.Guelpers@ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#ifndef Hadrons_MFermion_EMLepton_hpp_
#define Hadrons_MFermion_EMLepton_hpp_
#include <Hadrons/Global.hpp>
#include <Hadrons/Module.hpp>
#include <Hadrons/ModuleFactory.hpp>
BEGIN_HADRONS_NAMESPACE
/*******************************************************************************
*
* Calculates a free lepton propagator with a sequential insertion of
* i*\gamma_mu A_mu with a photon field A_mu
*
* L(x) = \sum_y S(x,y) i*\gamma_mu*A_mu S(y,xl) \delta_{(tl-x0),dt}
*
* with xl = (0,0,0,tl)
*
* In addition outputs the propagator without photon vertex
*
* L^{free}(x) = S(x,xl) \delta_{(tl-x0),dt}
*
*
* options:
* - action: fermion action used for propagator (string)
* - emField: photon field A_mu (string)
* - mass: input mass for the lepton propagator
* - twist: twisted boundary for lepton propagator, e.g. "0.0 0.0 0.0 0.5"
* - deltat: source-sink separation
*
*******************************************************************************/
/******************************************************************************
* EMLepton *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MFermion)
class EMLeptonPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(EMLeptonPar,
std::string, action,
std::string, emField,
double, mass,
std::string , boundary,
std::string, twist,
unsigned int, deltat);
};
template <typename FImpl>
class TEMLepton: public Module<EMLeptonPar>
{
public:
FERM_TYPE_ALIASES(FImpl,);
public:
typedef PhotonR::GaugeField EmField;
public:
// constructor
TEMLepton(const std::string name);
// destructor
virtual ~TEMLepton(void) {};
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
protected:
// setup
virtual void setup(void);
// execution
virtual void execute(void);
private:
unsigned int Ls_;
};
MODULE_REGISTER_TMP(EMLepton, TEMLepton<FIMPL>, MFermion);
/******************************************************************************
* TEMLepton implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TEMLepton<FImpl>::TEMLepton(const std::string name)
: Module<EMLeptonPar>(name)
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TEMLepton<FImpl>::getInput(void)
{
std::vector<std::string> in = {par().action, par().emField};
return in;
}
template <typename FImpl>
std::vector<std::string> TEMLepton<FImpl>::getOutput(void)
{
std::vector<std::string> out = {getName(), getName() + "_free"};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TEMLepton<FImpl>::setup(void)
{
Ls_ = env().getObjectLs(par().action);
envCreateLat(PropagatorField, getName());
envCreateLat(PropagatorField, getName() + "_free");
envTmpLat(FermionField, "source", Ls_);
envTmpLat(FermionField, "sol", Ls_);
envTmpLat(FermionField, "tmp");
envTmpLat(PropagatorField, "sourcetmp");
envTmpLat(PropagatorField, "proptmp");
envTmpLat(PropagatorField, "freetmp");
envTmp(Lattice<iScalar<vInteger>>, "tlat",1, envGetGrid(LatticeComplex));
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TEMLepton<FImpl>::execute(void)
{
LOG(Message) << "Computing free fermion propagator '" << getName() << "'"
<< std::endl;
auto &mat = envGet(FMat, par().action);
RealD mass = par().mass;
Complex ci(0.0,1.0);
PropagatorField &Aslashlep = envGet(PropagatorField, getName());
PropagatorField &lep = envGet(PropagatorField, getName() + "_free");
envGetTmp(FermionField, source);
envGetTmp(FermionField, sol);
envGetTmp(FermionField, tmp);
LOG(Message) << "Calculating a lepton Propagator with sequential Aslash insertion with lepton mass "
<< mass << " using the action '" << par().action
<< "' for fixed source-sink separation of " << par().deltat << std::endl;
envGetTmp(Lattice<iScalar<vInteger>>, tlat);
LatticeCoordinate(tlat, Tp);
std::vector<double> twist = strToVec<double>(par().twist);
if(twist.size() != Nd)
{
HADRONS_ERROR(Size, "number of twist angles does not match number of dimensions");
}
std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
if(boundary.size() != Nd)
{
HADRONS_ERROR(Size, "number of boundary conditions does not match number of dimensions");
}
auto &stoch_photon = envGet(EmField, par().emField);
unsigned int nt = env().getDim(Tp);
envGetTmp(PropagatorField, proptmp);
envGetTmp(PropagatorField, freetmp);
envGetTmp(PropagatorField, sourcetmp);
std::vector<int> position;
SitePropagator id;
id = 1.;
unsigned int tl=0;
//point source at (0,0,0,tl)
position.clear();
for(int tt=0;tt<Nd-1;tt++) position.push_back(0);
position.push_back(tl);
sourcetmp = zero;
pokeSite(id, sourcetmp, position);
//free propagator from pt source
for (unsigned int s = 0; s < Ns; ++s)
{
LOG(Message) << "Calculation for spin= " << s << std::endl;
if (Ls_ == 1)
{
PropToFerm<FImpl>(source, sourcetmp, s, 0);
}
else
{
PropToFerm<FImpl>(tmp, sourcetmp, s, 0);
// 5D source if action is 5d
mat.ImportPhysicalFermionSource(tmp, source);
}
sol = zero;
mat.FreePropagator(source,sol,mass,boundary,twist);
if (Ls_ == 1)
{
FermToProp<FImpl>(freetmp, sol, s, 0);
}
// create 4D propagators from 5D one if necessary
if (Ls_ > 1)
{
mat.ExportPhysicalFermionSolution(sol, tmp);
FermToProp<FImpl>(freetmp, tmp, s, 0);
}
}
for(tl=0;tl<nt;tl++){
//shift free propagator to different source positions
//account for possible anti-periodic boundary in time
proptmp = Cshift(freetmp,Tp, -tl);
proptmp = where( tlat < tl, boundary[Tp]*proptmp, proptmp);
// free propagator for fixed source-sink separation
lep = where(tlat == (tl-par().deltat+nt)%nt, proptmp, lep);
// i*A_mu*gamma_mu
sourcetmp = zero;
for(unsigned int mu=0;mu<=3;mu++)
{
Gamma gmu(Gamma::gmu[mu]);
sourcetmp += ci * PeekIndex<LorentzIndex>(stoch_photon, mu) * (gmu * proptmp );
}
proptmp = zero;
//sequential propagator from i*Aslash*S
LOG(Message) << "Sequential propagator for t= " << tl << std::endl;
for (unsigned int s = 0; s < Ns; ++s)
{
LOG(Message) << "Calculation for spin= " << s << std::endl;
if (Ls_ == 1)
{
PropToFerm<FImpl>(source, sourcetmp, s, 0);
}
else
{
PropToFerm<FImpl>(tmp, sourcetmp, s, 0);
// 5D source if action is 5d
mat.ImportPhysicalFermionSource(tmp, source);
}
sol = zero;
mat.FreePropagator(source,sol,mass,boundary,twist);
if (Ls_ == 1)
{
FermToProp<FImpl>(proptmp, sol, s, 0);
}
// create 4D propagators from 5D one if necessary
if (Ls_ > 1)
{
mat.ExportPhysicalFermionSolution(sol, tmp);
FermToProp<FImpl>(proptmp, tmp, s, 0);
}
}
// keep the result for the desired delta t
Aslashlep = where(tlat == (tl-par().deltat+nt)%nt, proptmp, Aslashlep);
}
//account for possible anti-periodic boundary in time
Aslashlep = where( tlat >= nt-par().deltat, boundary[Tp]*Aslashlep, Aslashlep);
lep = where( tlat >= nt-par().deltat, boundary[Tp]*lep, lep);
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MFermion_EMLepton_hpp_

13
Hadrons/Modules/MFermion/FreeProp.hpp
View File

@ -49,6 +49,7 @@ public:
std::string, source, std::string, source,
std::string, action, std::string, action,
double, mass, double, mass,
std::string , boundary,
std::string, twist); std::string, twist);
}; };
@ -168,8 +169,16 @@ void TFreeProp<FImpl>::execute(void)
} }
sol = zero; sol = zero;
std::vector<double> twist = strToVec<double>(par().twist); std::vector<double> twist = strToVec<double>(par().twist);
if(twist.size() != Nd) HADRONS_ERROR(Size, "number of twist angles does not match number of dimensions"); if(twist.size() != Nd)
mat.FreePropagator(source,sol,mass,twist); {
HADRONS_ERROR(Size, "number of twist angles does not match number of dimensions");
}
std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
if(boundary.size() != Nd)
{
HADRONS_ERROR(Size, "number of boundary conditions does not match number of dimensions");
}
mat.FreePropagator(source,sol,mass,boundary,twist);
FermToProp<FImpl>(prop, sol, s, c); FermToProp<FImpl>(prop, sol, s, c);
// create 4D propagators from 5D one if necessary // create 4D propagators from 5D one if necessary
if (Ls_ > 1) if (Ls_ > 1)

36
Hadrons/Modules/MSource/SeqAslash.cc Normal file
View File

@ -0,0 +1,36 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Hadrons/Modules/MSource/SeqAslash.cc
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
Author: Vera Guelpers <Vera.Guelpers@ed.ac.uk>
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 <Hadrons/Modules/MSource/SeqAslash.hpp>
using namespace Grid;
using namespace Hadrons;
using namespace MSource;
template class Grid::Hadrons::MSource::TSeqAslash<FIMPL>;

186
Hadrons/Modules/MSource/SeqAslash.hpp Normal file
View File

@ -0,0 +1,186 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: Hadrons/Modules/MSource/SeqAslash.hpp
Copyright (C) 2015-2018
Author: Antonin Portelli <antonin.portelli@me.com>
Author: Lanny91 <andrew.lawson@gmail.com>
Author: Vera Guelpers <Vera.Guelpers@ed.ac.uk>
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_MSource_SeqAslash_hpp_
#define Hadrons_MSource_SeqAslash_hpp_
#include <Hadrons/Global.hpp>
#include <Hadrons/Module.hpp>
#include <Hadrons/ModuleFactory.hpp>
BEGIN_HADRONS_NAMESPACE
/*
Sequential source
-----------------------------
* src_x = q_x * theta(x_3 - tA) * theta(tB - x_3) * i * A_mu g_mu * exp(i x.mom)
* options:
- q: input propagator (string)
- tA: begin timeslice (integer)
- tB: end timesilce (integer)
- emField: input photon field (string)
- mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0.")
*/
/******************************************************************************
* SeqAslash *
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MSource)
class SeqAslashPar: Serializable
{
public:
GRID_SERIALIZABLE_CLASS_MEMBERS(SeqAslashPar,
std::string, q,
unsigned int, tA,
unsigned int, tB,
std::string, emField,
std::string, mom);
};
template <typename FImpl>
class TSeqAslash: public Module<SeqAslashPar>
{
public:
FERM_TYPE_ALIASES(FImpl,);
public:
typedef PhotonR::GaugeField EmField;
public:
// constructor
TSeqAslash(const std::string name);
// destructor
virtual ~TSeqAslash(void) {};
// dependency relation
virtual std::vector<std::string> getInput(void);
virtual std::vector<std::string> getOutput(void);
protected:
// setup
virtual void setup(void);
// execution
virtual void execute(void);
private:
bool hasPhase_{false};
std::string momphName_, tName_;
};
MODULE_REGISTER_TMP(SeqAslash, TSeqAslash<FIMPL>, MSource);
/******************************************************************************
* TSeqAslash implementation *
******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
template <typename FImpl>
TSeqAslash<FImpl>::TSeqAslash(const std::string name)
: Module<SeqAslashPar>(name)
, momphName_ (name + "_momph")
, tName_ (name + "_t")
{}
// dependencies/products ///////////////////////////////////////////////////////
template <typename FImpl>
std::vector<std::string> TSeqAslash<FImpl>::getInput(void)
{
std::vector<std::string> in = {par().q,par().emField};
return in;
}
template <typename FImpl>
std::vector<std::string> TSeqAslash<FImpl>::getOutput(void)
{
std::vector<std::string> out = {getName()};
return out;
}
// setup ///////////////////////////////////////////////////////////////////////
template <typename FImpl>
void TSeqAslash<FImpl>::setup(void)
{
envCreateLat(PropagatorField, getName());
envCache(Lattice<iScalar<vInteger>>, tName_, 1, envGetGrid(LatticeComplex));
envCacheLat(LatticeComplex, momphName_);
envTmpLat(LatticeComplex, "coor");
}
// execution ///////////////////////////////////////////////////////////////////
template <typename FImpl>
void TSeqAslash<FImpl>::execute(void)
{
if (par().tA == par().tB)
{
LOG(Message) << "Generating Aslash sequential source at t= " << par().tA
<< " using the photon field " << par().emField << std::endl;
}
else
{
LOG(Message) << "Generating Aslash sequential source for "
<< par().tA << " <= t <= " << par().tB
<< " using the photon field " << par().emField << std::endl;
}
auto &src = envGet(PropagatorField, getName()); src=zero;
auto &q = envGet(PropagatorField, par().q);
auto &ph = envGet(LatticeComplex, momphName_);
auto &t = envGet(Lattice<iScalar<vInteger>>, tName_);
if (!hasPhase_)
{
Complex i(0.0,1.0);
std::vector<Real> p;
envGetTmp(LatticeComplex, coor);
p = strToVec<Real>(par().mom);
ph = zero;
for(unsigned int mu = 0; mu < env().getNd(); mu++)
{
LatticeCoordinate(coor, mu);
ph = ph + (p[mu]/env().getDim(mu))*coor;
}
ph = exp((Real)(2*M_PI)*i*ph);
LatticeCoordinate(t, Tp);
hasPhase_ = true;
}
auto &stoch_photon = envGet(EmField, par().emField);
Complex ci(0.0,1.0);
for(unsigned int mu=0;mu<=3;mu++)
{
Gamma gmu(Gamma::gmu[mu]);
src = src + where((t >= par().tA) and (t <= par().tB), ci * PeekIndex<LorentzIndex>(stoch_photon, mu) *ph*(gmu*q), 0.*q);
}
}
END_MODULE_NAMESPACE
END_HADRONS_NAMESPACE
#endif // Hadrons_MSource_SeqAslash_hpp_

6
Hadrons/modules.inc
View File

@ -1,6 +1,7 @@
modules_cc =\ modules_cc =\
Modules/MContraction/Baryon.cc \ Modules/MContraction/Baryon.cc \
Modules/MContraction/Meson.cc \ Modules/MContraction/Meson.cc \
Modules/MContraction/WeakMesonDecayKl2.cc \
Modules/MContraction/WeakEye3pt.cc \ Modules/MContraction/WeakEye3pt.cc \
Modules/MContraction/A2ALoop.cc \ Modules/MContraction/A2ALoop.cc \
Modules/MContraction/WeakNonEye3pt.cc \ Modules/MContraction/WeakNonEye3pt.cc \
@ -10,11 +11,13 @@ modules_cc =\
Modules/MContraction/Gamma3pt.cc \ Modules/MContraction/Gamma3pt.cc \
Modules/MFermion/FreeProp.cc \ Modules/MFermion/FreeProp.cc \
Modules/MFermion/GaugeProp.cc \ Modules/MFermion/GaugeProp.cc \
Modules/MFermion/EMLepton.cc \
Modules/MSource/Momentum.cc \ Modules/MSource/Momentum.cc \
Modules/MSource/Point.cc \ Modules/MSource/Point.cc \
Modules/MSource/Wall.cc \ Modules/MSource/Wall.cc \
Modules/MSource/SeqConserved.cc \ Modules/MSource/SeqConserved.cc \
Modules/MSource/SeqGamma.cc \ Modules/MSource/SeqGamma.cc \
Modules/MSource/SeqAslash.cc \
Modules/MSource/Z2.cc \ Modules/MSource/Z2.cc \
Modules/MSink/Point.cc \ Modules/MSink/Point.cc \
Modules/MSink/Smear.cc \ Modules/MSink/Smear.cc \
@ -72,10 +75,13 @@ modules_hpp =\
Modules/MContraction/Meson.hpp \ Modules/MContraction/Meson.hpp \
Modules/MContraction/DiscLoop.hpp \ Modules/MContraction/DiscLoop.hpp \
Modules/MContraction/Gamma3pt.hpp \ Modules/MContraction/Gamma3pt.hpp \
Modules/MContraction/WeakMesonDecayKl2.hpp \
Modules/MContraction/WeakNonEye3pt.hpp \ Modules/MContraction/WeakNonEye3pt.hpp \
Modules/MFermion/FreeProp.hpp \ Modules/MFermion/FreeProp.hpp \
Modules/MFermion/GaugeProp.hpp \ Modules/MFermion/GaugeProp.hpp \
Modules/MFermion/EMLepton.hpp \
Modules/MSource/SeqGamma.hpp \ Modules/MSource/SeqGamma.hpp \
Modules/MSource/SeqAslash.hpp \
Modules/MSource/Point.hpp \ Modules/MSource/Point.hpp \
Modules/MSource/Wall.hpp \ Modules/MSource/Wall.hpp \
Modules/MSource/Z2.hpp \ Modules/MSource/Z2.hpp \