1
0
mirror of https://github.com/paboyle/Grid.git synced 2024-11-10 07:55:35 +00:00

Adding some wilson loop support

This commit is contained in:
Azusa Yamaguchi 2015-06-05 10:02:36 +01:00
parent 8bd9fb4427
commit 58cdcbb5e4
4 changed files with 196 additions and 0 deletions

View File

@ -0,0 +1,20 @@
#ifndef QCD_GAUGE_ACTION_BASE
#define QCD_GAUGE_ACTION_BASE
namespace Grid {
namespace QCD{
template<class GaugeField>
class GaugeActionBase { // derive this from TermInAction?
public:
virtual RealD S(const GaugeField &U) = 0; // evaluate the action
virtual void deriv(const GaugeField &U,GaugeField & dSdU ) = 0; // evaluate the action derivative
virtual void staple(const GaugeField &U,GaugeField & dSdU ) = 0; // evaluate the action derivative
virtual void refresh(const GaugeField & ) {};
// Boundary conditions?
// Heatbath?
virtual ~GaugeActionBase() {};
};
}}
#endif

View File

@ -0,0 +1,24 @@
#ifndef QCD_WILSON_GAUGE_ACTION_H
#define QCD_WILSON_GAUGE_ACTION_H
////////////////////////////////////////////////////////////////////////
// Wilson Gauge Action .. should I template the Nc etc..
////////////////////////////////////////////////////////////////////////
template<class GaugeField,class MatrixField>
class WilsonGaugeAction : public GaugeActionBase<GaugeField> {
public:
virtual RealD S(const GaugeField &U) {
return WilsonLoops<MatrixField,GaugeField>::sumPlaquette(U);
};
virtual RealD deriv(GaugeField &U,GaugeField & dSdU ) {
WilsonLoops<MatrixField,GaugeField>::Staple(dSdU,U,mu);
};
virtual void staple(const MatrixField &stap,GaugeField & U,int mu ) {
WilsonLoops<MatrixField,GaugeField>::Staple(stap,U,mu);
};
};
#endif
#endif

View File

@ -0,0 +1,22 @@
#ifndef QCD_UTILS_COVARIANT_CSHIFT_H
#define QCD_UTILS_COVARIANT_CSHIFT_H
namespace Grid {
namespace QCD {
////////////////////////////////////////////////////////////////////////
// Low performance implementation of CovariantCshift API
////////////////////////////////////////////////////////////////////////
template<class covariant,class gauge> Lattice<covariant> CovShiftForward(const Lattice<gauge> &Link,
int mu,
const Lattice<covariant> &field)
{
return Link*Cshift(field,mu,1);// moves towards negative mu
}
template<class covariant,class gauge> Lattice<covariant> CovShiftBackward(const Lattice<gauge> &Link,
int mu,
const Lattice<covariant> &field)
{
Lattice<covariant> tmp(field._grid);
tmp = adj(Link)*field;
return Cshift(tmp,mu,-1);// moves towards positive mu
}
}}

130
lib/qcd/utils/WilsonLoops.h Normal file
View File

@ -0,0 +1,130 @@
#ifndef QCD_UTILS_WILSON_LOOPS_H
#define QCD_UTILS_WILSON_LOOPS_H
namespace Grid {
namespace QCD {
// Common wilson loop observables
template<class GaugeMat,class GaugeLorentz>
class WilsonLoops {
public:
//////////////////////////////////////////////////
// directed plaquette oriented in mu,nu plane
//////////////////////////////////////////////////
static void dirPlaquette(GaugeMat &plaq,const std::vector<GaugeMat> &U, const int mu, const int nu)
{
plaq=CovShiftForward(U[mu],mu,U[nu])*adj(CovShiftForward(U[nu],nu,U[mu]));
}
//////////////////////////////////////////////////
// trace of directed plaquette oriented in mu,nu plane
//////////////////////////////////////////////////
static void traceDirPlaquette(LatticeComplex &plaq, const std::vector<GaugeMat> &U, const int mu, const int nu)
{
GaugeMat sp(U[0]._grid);
dirPlaquette(sp,U,mu,nu);
plaq=trace(sp);
}
//////////////////////////////////////////////////
// sum over all planes of plaquette
//////////////////////////////////////////////////
static void sitePlaquette(LatticeComplex &Plaq,const std::vector<GaugeMat> &U)
{
LatticeComplex sitePlaq(U[0]._grid);
for(int mu=1;mu<Nd;mu++){
for(int nu=0;nu<mu;nu++){
traceDirPlaquette(sitePlaq,U,mu,nu);
Plaq = Plaq + sitePlaq;
}
}
}
//////////////////////////////////////////////////
// sum over all x,y,z,t and over all planes of plaquette
//////////////////////////////////////////////////
static RealD sumPlaquette(const GaugeLorentz &Umu){
std::vector<GaugeMat> U(4,Umu._grid);
for(int mu=0;mu<Nd;mu++){
U[mu] = peekIndex<LorentzIndex>(Umu,mu);
}
LatticeComplex Plaq(Umu._grid);
sitePlaquette(Plaq,U);
TComplex Tp = sum(Plaq);
Complex p = TensorRemove(Tp);
return p.real();
}
//////////////////////////////////////////////////
// average over all x,y,z,t and over all planes of plaquette
//////////////////////////////////////////////////
static RealD avgPlaquette(const GaugeLorentz &Umu){
RealD sumplaq = sumPlaquette(Umu);
double vol = Umu._grid->gSites();
double faces = (1.0*Nd*(Nd-1))/2.0;
return sumplaq/vol/faces/Nc; // Nd , Nc dependent... FIXME
}
//////////////////////////////////////////////////
// the sum over all staples on each site
//////////////////////////////////////////////////
static void Staple(GaugeMat &staple,GaugeLorentz &Umu,int mu){
std::vector<GaugeMat> U(4,Umu._grid);
for(int d=0;d<Nd;d++){
U[d] = peekIndex<LorentzIndex>(Umu,d);
}
staple = zero;
GaugeMat tmp(U[0]._grid);
for(int nu=0;nu<Nd;nu++){
if(nu != mu) {
// mu
// ^
// |__ nu
// __ __
// | |
// __| = __| *
//
staple += CovShiftForward(U[nu],nu,U[mu])*Cshift(adj(U[nu]),mu,+1);
//
// __ __
// | |
// |__ = | * __
//
//
tmp = CovShiftForward (U[mu],mu,U[nu]);
staple+= CovShiftBackward(U[nu],nu,tmp);
}
}
}
//////////////////////////////////////////////////////
// Similar to above for rectangle is required
//////////////////////////////////////////////////////
/*
void siteRectangle(GaugeMat &plaq,const std::vector<GaugeMat> &U, const int mu, const int nu){
RealD avgRectangle(const std::vector<GaugeMat> &U){}
RealD avgRectangle(const std::vector<GaugeMat> &U, const int mu, const int nu){}
void traceRectangle(LatticeComplex &plaq,const std::vector<GaugeMat> &U, const int mu, const int nu){}
void siteRectangle(GaugeMat &plaq,const std::vector<GaugeMat> &U, const int mu, const int nu){}
*/
};
}}
#endif