portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-11-04 05:54:32 +00:00
Code Releases Activity

Merge branch 'develop' into feature/clover

Browse Source
This commit is contained in:
Guido Cossu
2017-08-04 12:19:57 +01:00
parent 62a64d9108 175f393f9d
commit fde71c3c52
296 changed files with 34674 additions and 9453 deletions
Download Patch File Download Diff File Expand all files Collapse all files

197
lib/qcd/utils/CovariantLaplacian.h Normal file
View File

@@ -0,0 +1,197 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/action/scalar/CovariantLaplacian.h
Copyright (C) 2016
Author: Guido Cossu <guido.cossu@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 COVARIANT_LAPLACIAN_H
#define COVARIANT_LAPLACIAN_H
namespace Grid {
namespace QCD {
struct LaplacianParams : Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(LaplacianParams,
RealD, lo,
RealD, hi,
int, MaxIter,
RealD, tolerance,
int, degree,
int, precision);
// constructor
LaplacianParams(RealD lo = 0.0,
RealD hi = 1.0,
int maxit = 1000,
RealD tol = 1.0e-8,
int degree = 10,
int precision = 64)
: lo(lo),
hi(hi),
MaxIter(maxit),
tolerance(tol),
degree(degree),
precision(precision){};
};
////////////////////////////////////////////////////////////
// Laplacian operator L on adjoint fields
//
// phi: adjoint field
// L: D_mu^dag D_mu
//
// L phi(x) = Sum_mu [ U_mu(x)phi(x+mu)U_mu(x)^dag +
// U_mu(x-mu)^dag phi(x-mu)U_mu(x-mu)
// -2phi(x)]
//
// Operator designed to be encapsulated by
// an HermitianLinearOperator<.. , ..>
////////////////////////////////////////////////////////////
template <class Impl>
class LaplacianAdjointField: public Metric<typename Impl::Field> {
OperatorFunction<typename Impl::Field> &Solver;
LaplacianParams param;
MultiShiftFunction PowerHalf;
MultiShiftFunction PowerInvHalf;
public:
INHERIT_GIMPL_TYPES(Impl);
LaplacianAdjointField(GridBase* grid, OperatorFunction<GaugeField>& S, LaplacianParams& p, const RealD k = 1.0)
: U(Nd, grid), Solver(S), param(p), kappa(k){
AlgRemez remez(param.lo,param.hi,param.precision);
std::cout<<GridLogMessage << "Generating degree "<<param.degree<<" for x^(1/2)"<<std::endl;
remez.generateApprox(param.degree,1,2);
PowerHalf.Init(remez,param.tolerance,false);
PowerInvHalf.Init(remez,param.tolerance,true);
};
void Mdir(const GaugeField&, GaugeField&, int, int){ assert(0);}
void Mdiag(const GaugeField&, GaugeField&){ assert(0);}
void ImportGauge(const GaugeField& _U) {
for (int mu = 0; mu < Nd; mu++) {
U[mu] = PeekIndex<LorentzIndex>(_U, mu);
}
}
void M(const GaugeField& in, GaugeField& out) {
// in is an antihermitian matrix
// test
//GaugeField herm = in + adj(in);
//std::cout << "AHermiticity: " << norm2(herm) << std::endl;
GaugeLinkField tmp(in._grid);
GaugeLinkField tmp2(in._grid);
GaugeLinkField sum(in._grid);
for (int nu = 0; nu < Nd; nu++) {
sum = zero;
GaugeLinkField in_nu = PeekIndex<LorentzIndex>(in, nu);
GaugeLinkField out_nu(out._grid);
for (int mu = 0; mu < Nd; mu++) {
tmp = U[mu] * Cshift(in_nu, mu, +1) * adj(U[mu]);
tmp2 = adj(U[mu]) * in_nu * U[mu];
sum += tmp + Cshift(tmp2, mu, -1) - 2.0 * in_nu;
}
out_nu = (1.0 - kappa) * in_nu - kappa / (double(4 * Nd)) * sum;
PokeIndex<LorentzIndex>(out, out_nu, nu);
}
}
void MDeriv(const GaugeField& in, GaugeField& der) {
// in is anti-hermitian
RealD factor = -kappa / (double(4 * Nd));
for (int mu = 0; mu < Nd; mu++){
GaugeLinkField der_mu(der._grid);
der_mu = zero;
for (int nu = 0; nu < Nd; nu++){
GaugeLinkField in_nu = PeekIndex<LorentzIndex>(in, nu);
der_mu += U[mu] * Cshift(in_nu, mu, 1) * adj(U[mu]) * in_nu;
}
// the minus sign comes by using the in_nu instead of the
// adjoint in the last multiplication
PokeIndex<LorentzIndex>(der, -2.0 * factor * der_mu, mu);
}
}
// separating this temporarily
void MDeriv(const GaugeField& left, const GaugeField& right,
GaugeField& der) {
// in is anti-hermitian
RealD factor = -kappa / (double(4 * Nd));
for (int mu = 0; mu < Nd; mu++) {
GaugeLinkField der_mu(der._grid);
der_mu = zero;
for (int nu = 0; nu < Nd; nu++) {
GaugeLinkField left_nu = PeekIndex<LorentzIndex>(left, nu);
GaugeLinkField right_nu = PeekIndex<LorentzIndex>(right, nu);
der_mu += U[mu] * Cshift(left_nu, mu, 1) * adj(U[mu]) * right_nu;
der_mu += U[mu] * Cshift(right_nu, mu, 1) * adj(U[mu]) * left_nu;
}
PokeIndex<LorentzIndex>(der, -factor * der_mu, mu);
}
}
void Minv(const GaugeField& in, GaugeField& inverted){
HermitianLinearOperator<LaplacianAdjointField<Impl>,GaugeField> HermOp(*this);
Solver(HermOp, in, inverted);
}
void MSquareRoot(GaugeField& P){
GaugeField Gp(P._grid);
HermitianLinearOperator<LaplacianAdjointField<Impl>,GaugeField> HermOp(*this);
ConjugateGradientMultiShift<GaugeField> msCG(param.MaxIter,PowerHalf);
msCG(HermOp,P,Gp);
P = Gp;
}
void MInvSquareRoot(GaugeField& P){
GaugeField Gp(P._grid);
HermitianLinearOperator<LaplacianAdjointField<Impl>,GaugeField> HermOp(*this);
ConjugateGradientMultiShift<GaugeField> msCG(param.MaxIter,PowerInvHalf);
msCG(HermOp,P,Gp);
P = Gp;
}
private:
RealD kappa;
std::vector<GaugeLinkField> U;
};
}
}
#endif

190
lib/qcd/utils/GaugeFix.h Normal file
View File

@@ -0,0 +1,190 @@
/*************************************************************************************
grid` physics library, www.github.com/paboyle/Grid
Copyright (C) 2015
Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
Author: Peter Boyle <paboyle@ph.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 <Grid/Grid.h>
namespace Grid {
namespace QCD {
template <class Gimpl>
class FourierAcceleratedGaugeFixer : public Gimpl {
public:
INHERIT_GIMPL_TYPES(Gimpl);
typedef typename Gimpl::GaugeLinkField GaugeMat;
typedef typename Gimpl::GaugeField GaugeLorentz;
static void GaugeLinkToLieAlgebraField(const std::vector<GaugeMat> &U,std::vector<GaugeMat> &A) {
for(int mu=0;mu<Nd;mu++){
Complex cmi(0.0,-1.0);
A[mu] = Ta(U[mu]) * cmi;
}
}
static void DmuAmu(const std::vector<GaugeMat> &A,GaugeMat &dmuAmu) {
dmuAmu=zero;
for(int mu=0;mu<Nd;mu++){
dmuAmu = dmuAmu + A[mu] - Cshift(A[mu],mu,-1);
}
}
static void SteepestDescentGaugeFix(GaugeLorentz &Umu,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false) {
GridBase *grid = Umu._grid;
Real org_plaq =WilsonLoops<Gimpl>::avgPlaquette(Umu);
Real org_link_trace=WilsonLoops<Gimpl>::linkTrace(Umu);
Real old_trace = org_link_trace;
Real trG;
std::vector<GaugeMat> U(Nd,grid);
GaugeMat dmuAmu(grid);
for(int i=0;i<maxiter;i++){
for(int mu=0;mu<Nd;mu++) U[mu]= PeekIndex<LorentzIndex>(Umu,mu);
if ( Fourier==false ) {
trG = SteepestDescentStep(U,alpha,dmuAmu);
} else {
trG = FourierAccelSteepestDescentStep(U,alpha,dmuAmu);
}
for(int mu=0;mu<Nd;mu++) PokeIndex<LorentzIndex>(Umu,U[mu],mu);
// Monitor progress and convergence test
// infrequently to minimise cost overhead
if ( i %20 == 0 ) {
Real plaq =WilsonLoops<Gimpl>::avgPlaquette(Umu);
Real link_trace=WilsonLoops<Gimpl>::linkTrace(Umu);
if (Fourier)
std::cout << GridLogMessage << "Fourier Iteration "<<i<< " plaq= "<<plaq<< " dmuAmu " << norm2(dmuAmu)<< std::endl;
else
std::cout << GridLogMessage << " Iteration "<<i<< " plaq= "<<plaq<< " dmuAmu " << norm2(dmuAmu)<< std::endl;
Real Phi = 1.0 - old_trace / link_trace ;
Real Omega= 1.0 - trG;
std::cout << GridLogMessage << " Iteration "<<i<< " Phi= "<<Phi<< " Omega= " << Omega<< " trG " << trG <<std::endl;
if ( (Omega < Omega_tol) && ( ::fabs(Phi) < Phi_tol) ) {
std::cout << GridLogMessage << "Converged ! "<<std::endl;
return;
}
old_trace = link_trace;
}
}
};
static Real SteepestDescentStep(std::vector<GaugeMat> &U,Real & alpha, GaugeMat & dmuAmu) {
GridBase *grid = U[0]._grid;
std::vector<GaugeMat> A(Nd,grid);
GaugeMat g(grid);
GaugeLinkToLieAlgebraField(U,A);
ExpiAlphaDmuAmu(A,g,alpha,dmuAmu);
Real vol = grid->gSites();
Real trG = TensorRemove(sum(trace(g))).real()/vol/Nc;
SU<Nc>::GaugeTransform(U,g);
return trG;
}
static Real FourierAccelSteepestDescentStep(std::vector<GaugeMat> &U,Real & alpha, GaugeMat & dmuAmu) {
GridBase *grid = U[0]._grid;
Real vol = grid->gSites();
FFT theFFT((GridCartesian *)grid);
LatticeComplex Fp(grid);
LatticeComplex psq(grid); psq=zero;
LatticeComplex pmu(grid);
LatticeComplex one(grid); one = Complex(1.0,0.0);
GaugeMat g(grid);
GaugeMat dmuAmu_p(grid);
std::vector<GaugeMat> A(Nd,grid);
GaugeLinkToLieAlgebraField(U,A);
DmuAmu(A,dmuAmu);
theFFT.FFT_all_dim(dmuAmu_p,dmuAmu,FFT::forward);
//////////////////////////////////
// Work out Fp = psq_max/ psq...
//////////////////////////////////
std::vector<int> latt_size = grid->GlobalDimensions();
std::vector<int> coor(grid->_ndimension,0);
for(int mu=0;mu<Nd;mu++) {
Real TwoPiL = M_PI * 2.0/ latt_size[mu];
LatticeCoordinate(pmu,mu);
pmu = TwoPiL * pmu ;
psq = psq + 4.0*sin(pmu*0.5)*sin(pmu*0.5);
}
Complex psqMax(16.0);
Fp = psqMax*one/psq;
/*
static int once;
if ( once == 0 ) {
std::cout << " Fp " << Fp <<std::endl;
once ++;
}*/
pokeSite(TComplex(1.0),Fp,coor);
dmuAmu_p = dmuAmu_p * Fp;
theFFT.FFT_all_dim(dmuAmu,dmuAmu_p,FFT::backward);
GaugeMat ciadmam(grid);
Complex cialpha(0.0,-alpha);
ciadmam = dmuAmu*cialpha;
SU<Nc>::taExp(ciadmam,g);
Real trG = TensorRemove(sum(trace(g))).real()/vol/Nc;
SU<Nc>::GaugeTransform(U,g);
return trG;
}
static void ExpiAlphaDmuAmu(const std::vector<GaugeMat> &A,GaugeMat &g,Real & alpha, GaugeMat &dmuAmu) {
GridBase *grid = g._grid;
Complex cialpha(0.0,-alpha);
GaugeMat ciadmam(grid);
DmuAmu(A,dmuAmu);
ciadmam = dmuAmu*cialpha;
SU<Nc>::taExp(ciadmam,g);
}
};
}
}

226
lib/qcd/utils/Metric.h Normal file
View File

@@ -0,0 +1,226 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/hmc/integrators/Integrator.h
Copyright (C) 2015
Author: Guido Cossu <guido.cossu@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 METRIC_H
#define METRIC_H
namespace Grid{
namespace QCD{
template <typename Field>
class Metric{
public:
virtual void ImportGauge(const Field&) = 0;
virtual void M(const Field&, Field&) = 0;
virtual void Minv(const Field&, Field&) = 0;
virtual void MSquareRoot(Field&) = 0;
virtual void MInvSquareRoot(Field&) = 0;
virtual void MDeriv(const Field&, Field&) = 0;
virtual void MDeriv(const Field&, const Field&, Field&) = 0;
};
// Need a trivial operator
template <typename Field>
class TrivialMetric : public Metric<Field>{
public:
virtual void ImportGauge(const Field&){};
virtual void M(const Field& in, Field& out){
out = in;
}
virtual void Minv(const Field& in, Field& out){
out = in;
}
virtual void MSquareRoot(Field& P){
// do nothing
}
virtual void MInvSquareRoot(Field& P){
// do nothing
}
virtual void MDeriv(const Field& in, Field& out){
out = zero;
}
virtual void MDeriv(const Field& left, const Field& right, Field& out){
out = zero;
}
};
///////////////////////////////
// Generalised momenta
///////////////////////////////
template <typename Implementation>
class GeneralisedMomenta{
public:
typedef typename Implementation::Field MomentaField; //for readability
typedef typename Implementation::GaugeLinkField MomentaLinkField; //for readability
Metric<MomentaField>& M;
MomentaField Mom;
// Auxiliary fields
// not hard coded but inherit the type from the metric
// created Nd new fields
// hide these in the metric?
//typedef Lattice<iVector<iScalar<iMatrix<vComplex, Nc> >, Nd/2 > > AuxiliaryMomentaType;
MomentaField AuxMom;
MomentaField AuxField;
GeneralisedMomenta(GridBase* grid, Metric<MomentaField>& M): M(M), Mom(grid), AuxMom(grid), AuxField(grid){}
// Correct
void MomentaDistribution(GridParallelRNG& pRNG){
// Generate a distribution for
// P^dag G P
// where G = M^-1
// Generate gaussian momenta
Implementation::generate_momenta(Mom, pRNG);
// Modify the distribution with the metric
M.MSquareRoot(Mom);
if (1) {
// Auxiliary momenta
// do nothing if trivial, so hide in the metric
MomentaField AuxMomTemp(Mom._grid);
Implementation::generate_momenta(AuxMom, pRNG);
Implementation::generate_momenta(AuxField, pRNG);
// Modify the distribution with the metric
// Aux^dag M Aux
M.MInvSquareRoot(AuxMom); // AuxMom = M^{-1/2} AuxMomTemp
}
}
// Correct
RealD MomentaAction(){
MomentaField inv(Mom._grid);
inv = zero;
M.Minv(Mom, inv);
LatticeComplex Hloc(Mom._grid);
Hloc = zero;
for (int mu = 0; mu < Nd; mu++) {
// This is not very general
// hide in the metric
auto Mom_mu = PeekIndex<LorentzIndex>(Mom, mu);
auto inv_mu = PeekIndex<LorentzIndex>(inv, mu);
Hloc += trace(Mom_mu * inv_mu);
}
if (1) {
// Auxiliary Fields
// hide in the metric
M.M(AuxMom, inv);
for (int mu = 0; mu < Nd; mu++) {
// This is not very general
// hide in the operators
auto inv_mu = PeekIndex<LorentzIndex>(inv, mu);
auto am_mu = PeekIndex<LorentzIndex>(AuxMom, mu);
auto af_mu = PeekIndex<LorentzIndex>(AuxField, mu);
Hloc += trace(am_mu * inv_mu);// p M p
Hloc += trace(af_mu * af_mu);
}
}
Complex Hsum = sum(Hloc);
return Hsum.real();
}
// Correct
void DerivativeU(MomentaField& in, MomentaField& der){
// Compute the derivative of the kinetic term
// with respect to the gauge field
MomentaField MDer(in._grid);
MomentaField X(in._grid);
X = zero;
M.Minv(in, X); // X = G in
M.MDeriv(X, MDer); // MDer = U * dS/dU
der = Implementation::projectForce(MDer); // Ta if gauge fields
}
void AuxiliaryFieldsDerivative(MomentaField& der){
der = zero;
if (1){
// Auxiliary fields
MomentaField der_temp(der._grid);
MomentaField X(der._grid);
X=zero;
//M.M(AuxMom, X); // X = M Aux
// Two derivative terms
// the Mderiv need separation of left and right terms
M.MDeriv(AuxMom, der);
// this one should not be necessary (identical to the previous one)
//M.MDeriv(X, AuxMom, der_temp); der += der_temp;
der = -1.0*Implementation::projectForce(der);
}
}
void DerivativeP(MomentaField& der){
der = zero;
M.Minv(Mom, der);
// is the projection necessary here?
// no for fields in the algebra
der = Implementation::projectForce(der);
}
void update_auxiliary_momenta(RealD ep){
if(1){
AuxMom -= ep * AuxField;
}
}
void update_auxiliary_fields(RealD ep){
if (1) {
MomentaField tmp(AuxMom._grid);
MomentaField tmp2(AuxMom._grid);
M.M(AuxMom, tmp);
// M.M(tmp, tmp2);
AuxField += ep * tmp; // M^2 AuxMom
// factor of 2?
}
}
};
}
}
#endif //METRIC_H

6
lib/qcd/utils/SUn.h
View File

@@ -170,6 +170,7 @@ class SU {
ta()()(i2, i1) = 1.0;
ta = ta * 0.5;
}
template <class cplx>
static void generatorSigmaX(int su2Index, iSUnMatrix<cplx> &ta) {
ta = zero;
@@ -194,6 +195,8 @@ class SU {
ta = ta * nrm;
}
////////////////////////////////////////////////////////////////////////
// Map a su2 subgroup number to the pair of rows that are non zero
////////////////////////////////////////////////////////////////////////
@@ -713,8 +716,7 @@ template<typename GaugeField,typename GaugeMat>
for (int a = 0; a < AdjointDimension; a++) {
generator(a, Ta);
auto tmp = - 2.0 * (trace(timesI(Ta) * in)) * scale;// 2.0 for the normalization of the trace in the fundamental rep
pokeColour(h_out, tmp, a);
pokeColour(h_out, - 2.0 * (trace(timesI(Ta) * in)) * scale, a);
}
}

96
lib/qcd/utils/ScalarObjs.h Normal file
View File

@@ -0,0 +1,96 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./lib/qcd/utils/WilsonLoops.h
Copyright (C) 2015
Author: neo <cossu@post.kek.jp>
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 SCALAR_OBJS_H
#define SCALAR_OBJS_H
namespace Grid {
// FIXME drop the QCD namespace in Nd
// Scalar field obs
template <class Impl>
class ScalarObs {
public:
//////////////////////////////////////////////////
// squared field
//////////////////////////////////////////////////
static void phisquared(typename Impl::Field &fsq,
const typename Impl::Field &f) {
fsq = f * f;
}
//////////////////////////////////////////////////
// phi^4 interaction term
//////////////////////////////////////////////////
static void phifourth(typename Impl::Field &fsq,
const typename Impl::Field &f) {
fsq = f * f * f * f;
}
//////////////////////////////////////////////////
// phi(x)phi(x+mu)
//////////////////////////////////////////////////
static void phider(typename Impl::Field &fsq,
const typename Impl::Field &f) {
fsq = Cshift(f, 0, -1) * f;
for (int mu = 1; mu < QCD::Nd; mu++) fsq += Cshift(f, mu, -1) * f;
}
//////////////////////////////////////////////////
// Vol sum of the previous obs.
//////////////////////////////////////////////////
static RealD sumphider(const typename Impl::Field &f) {
typename Impl::Field tmp(f._grid);
tmp = Cshift(f, 0, -1) * f;
for (int mu = 1; mu < QCD::Nd; mu++) {
tmp += Cshift(f, mu, -1) * f;
}
return -sum(trace(tmp));
}
static RealD sumphisquared(const typename Impl::Field &f) {
typename Impl::Field tmp(f._grid);
tmp = f * f;
return sum(trace(tmp));
}
static RealD sumphifourth(const typename Impl::Field &f) {
typename Impl::Field tmp(f._grid);
phifourth(tmp, f);
return sum(trace(tmp));
}
};
}
#endif

6
lib/qcd/utils/Utils.h
View File

@@ -3,7 +3,13 @@
#include <Grid/qcd/utils/SpaceTimeGrid.h>
#include <Grid/qcd/utils/LinalgUtils.h>
#include <Grid/qcd/utils/CovariantCshift.h>
// Scalar field
#include <Grid/qcd/utils/ScalarObjs.h>
// Include representations
#include <Grid/qcd/utils/SUn.h>
#include <Grid/qcd/utils/SUnAdjoint.h>
#include <Grid/qcd/utils/SUnTwoIndex.h>
#endif

127
lib/qcd/utils/WilsonLoops.h
View File

@@ -54,14 +54,26 @@ public:
// resolution throughout the usage in this file, and rather defeats the
// purpose of deriving
// from Gimpl.
/*
plaq = Gimpl::CovShiftBackward(
U[mu], mu, Gimpl::CovShiftBackward(
U[nu], nu, Gimpl::CovShiftForward(U[mu], mu, U[nu])));
*/
// _
//|< _|
plaq = Gimpl::CovShiftForward(U[mu],mu,
Gimpl::CovShiftForward(U[nu],nu,
Gimpl::CovShiftBackward(U[mu],mu,
Gimpl::CovShiftIdentityBackward(U[nu], nu))));
}
//////////////////////////////////////////////////
// trace of directed plaquette oriented in mu,nu plane
//////////////////////////////////////////////////
static void traceDirPlaquette(LatticeComplex &plaq,
static void traceDirPlaquette(ComplexField &plaq,
const std::vector<GaugeMat> &U, const int mu,
const int nu) {
GaugeMat sp(U[0]._grid);
@@ -71,9 +83,9 @@ public:
//////////////////////////////////////////////////
// sum over all planes of plaquette
//////////////////////////////////////////////////
static void sitePlaquette(LatticeComplex &Plaq,
static void sitePlaquette(ComplexField &Plaq,
const std::vector<GaugeMat> &U) {
LatticeComplex sitePlaq(U[0]._grid);
ComplexField sitePlaq(U[0]._grid);
Plaq = zero;
for (int mu = 1; mu < Nd; mu++) {
for (int nu = 0; nu < mu; nu++) {
@@ -86,20 +98,21 @@ public:
// 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);
std::vector<GaugeMat> U(Nd, Umu._grid);
// inefficient here
for (int mu = 0; mu < Nd; mu++) {
U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
}
LatticeComplex Plaq(Umu._grid);
ComplexField Plaq(Umu._grid);
sitePlaquette(Plaq, U);
TComplex Tp = sum(Plaq);
Complex p = TensorRemove(Tp);
auto Tp = sum(Plaq);
auto p = TensorRemove(Tp);
return p.real();
}
//////////////////////////////////////////////////
// average over all x,y,z,t and over all planes of plaquette
//////////////////////////////////////////////////
@@ -114,17 +127,17 @@ public:
// average over traced single links
//////////////////////////////////////////////////
static RealD linkTrace(const GaugeLorentz &Umu) {
std::vector<GaugeMat> U(4, Umu._grid);
std::vector<GaugeMat> U(Nd, Umu._grid);
LatticeComplex Tr(Umu._grid);
ComplexField Tr(Umu._grid);
Tr = zero;
for (int mu = 0; mu < Nd; mu++) {
U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
Tr = Tr + trace(U[mu]);
}
TComplex Tp = sum(Tr);
Complex p = TensorRemove(Tp);
auto Tp = sum(Tr);
auto p = TensorRemove(Tp);
double vol = Umu._grid->gSites();
@@ -139,7 +152,7 @@ public:
GridBase *grid = Umu._grid;
std::vector<GaugeMat> U(4, grid);
std::vector<GaugeMat> U(Nd, grid);
for (int d = 0; d < Nd; d++) {
U[d] = PeekIndex<LorentzIndex>(Umu, d);
}
@@ -175,6 +188,32 @@ public:
}
}
// For the force term
static void StapleMult(GaugeMat &staple, const GaugeLorentz &Umu, int mu) {
GridBase *grid = Umu._grid;
std::vector<GaugeMat> U(Nd, grid);
for (int d = 0; d < Nd; d++) {
// this operation is taking too much time
U[d] = PeekIndex<LorentzIndex>(Umu, d);
}
staple = zero;
GaugeMat tmp1(grid);
GaugeMat tmp2(grid);
for (int nu = 0; nu < Nd; nu++) {
if (nu != mu) {
// this is ~10% faster than the Staple
tmp1 = Cshift(U[nu], mu, 1);
tmp2 = Cshift(U[mu], nu, 1);
staple += tmp1* adj(U[nu]*tmp2);
tmp2 = adj(U[mu]*tmp1)*U[nu];
staple += Cshift(tmp2, nu, -1);
}
}
staple = U[mu]*staple;
}
//////////////////////////////////////////////////
// the sum over all staples on each site
//////////////////////////////////////////////////
@@ -187,7 +226,6 @@ public:
U[d] = PeekIndex<LorentzIndex>(Umu, d);
}
staple = zero;
GaugeMat tmp(grid);
for (int nu = 0; nu < Nd; nu++) {
@@ -201,7 +239,7 @@ public:
// |
// __|
//
staple += Gimpl::ShiftStaple(
Gimpl::CovShiftForward(
U[nu], nu,
@@ -214,10 +252,10 @@ public:
// |__
//
//
staple += Gimpl::ShiftStaple(
Gimpl::CovShiftBackward(U[nu], nu,
Gimpl::CovShiftBackward(U[mu], mu, U[nu])),
mu);
Gimpl::CovShiftBackward(U[mu], mu, U[nu])), mu);
}
}
}
@@ -227,15 +265,12 @@ public:
//////////////////////////////////////////////////
static void StapleUpper(GaugeMat &staple, const GaugeLorentz &Umu, int mu,
int nu) {
staple = zero;
if (nu != mu) {
GridBase *grid = Umu._grid;
std::vector<GaugeMat> U(4, grid);
std::vector<GaugeMat> U(Nd, grid);
for (int d = 0; d < Nd; d++) {
U[d] = PeekIndex<LorentzIndex>(Umu, d);
U[d] = PeekIndex<LorentzIndex>(Umu, d);// some redundant copies
}
// mu
@@ -247,7 +282,7 @@ public:
// __|
//
staple += Gimpl::ShiftStaple(
staple = Gimpl::ShiftStaple(
Gimpl::CovShiftForward(
U[nu], nu,
Gimpl::CovShiftBackward(
@@ -282,6 +317,7 @@ public:
Gimpl::CovShiftBackward(U[nu], nu,
Gimpl::CovShiftBackward(U[mu], mu, U[nu])),
mu);
}
}
@@ -296,14 +332,35 @@ public:
// +--<--+ +--<--+
GaugeMat Vup(Umu._grid), Vdn(Umu._grid);
StapleUpper(Vup, Umu, mu, nu);// coalesce these two (up low)
StapleUpper(Vup, Umu, mu, nu);
StapleLower(Vdn, Umu, mu, nu);
GaugeMat v = adj(Vup) - adj(Vdn);
GaugeMat v = Vup - Vdn;
GaugeMat u = PeekIndex<LorentzIndex>(Umu, mu); // some redundant copies
GaugeMat vu = v*u;
FS = 0.25*Ta(u*v + Cshift(vu, mu, +1));
FS = 0.25*Ta(u*v + Cshift(vu, mu, -1));
}
static Real TopologicalCharge(GaugeLorentz &U){
// 4d topological charge
assert(Nd==4);
// Bx = -iF(y,z), By = -iF(z,y), Bz = -iF(x,y)
GaugeMat Bx(U._grid), By(U._grid), Bz(U._grid);
FieldStrength(Bx, U, Ydir, Zdir);
FieldStrength(By, U, Zdir, Xdir);
FieldStrength(Bz, U, Xdir, Ydir);
// Ex = -iF(t,x), Ey = -iF(t,y), Ez = -iF(t,z)
GaugeMat Ex(U._grid), Ey(U._grid), Ez(U._grid);
FieldStrength(Ex, U, Tdir, Xdir);
FieldStrength(Ey, U, Tdir, Ydir);
FieldStrength(Ez, U, Tdir, Zdir);
double coeff = 8.0/(32.0*M_PI*M_PI);
ComplexField qfield = coeff*trace(Bx*Ex + By*Ey + Bz*Ez);
auto Tq = sum(qfield);
return TensorRemove(Tq).real();
}
//////////////////////////////////////////////////////
@@ -321,16 +378,16 @@ public:
adj(Gimpl::CovShiftForward(
U[nu], nu, Gimpl::CovShiftForward(U[nu], nu, U[mu])));
}
static void traceDirRectangle(LatticeComplex &rect,
static void traceDirRectangle(ComplexField &rect,
const std::vector<GaugeMat> &U, const int mu,
const int nu) {
GaugeMat sp(U[0]._grid);
dirRectangle(sp, U, mu, nu);
rect = trace(sp);
}
static void siteRectangle(LatticeComplex &Rect,
static void siteRectangle(ComplexField &Rect,
const std::vector<GaugeMat> &U) {
LatticeComplex siteRect(U[0]._grid);
ComplexField siteRect(U[0]._grid);
Rect = zero;
for (int mu = 1; mu < Nd; mu++) {
for (int nu = 0; nu < mu; nu++) {
@@ -350,12 +407,12 @@ public:
U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
}
LatticeComplex Rect(Umu._grid);
ComplexField Rect(Umu._grid);
siteRectangle(Rect, U);
TComplex Tp = sum(Rect);
Complex p = TensorRemove(Tp);
auto Tp = sum(Rect);
auto p = TensorRemove(Tp);
return p.real();
}
//////////////////////////////////////////////////
@@ -425,8 +482,8 @@ public:
// |___ ___|
//
// tmp= Staple2x1* Cshift(U[mu],mu,-2);
// Stap+= Cshift(tmp,mu,1) ;
// tmp= Staple2x1* Cshift(U[mu],mu,-2);
// Stap+= Cshift(tmp,mu,1) ;
Stap += Cshift(Staple2x1, mu, 1) * Cshift(U[mu], mu, -1);
;