staples_eo.cpp

Go to the documentation of this file.

#include "staples_eo.h"

const std::string Staples_eo::class_name = "Staples_eo";

//====================================================================
double Staples_eo::plaquette(const Field_G *U)
{
  return (plaq_s(U) + plaq_t(U)) / 2;
}


//====================================================================
double Staples_eo::plaq_s(const Field_G *U)
{
  double         plaq = 0.0;
  static Field_G staple;

  staple = upper(U, 0, 1);
  // staple = lower(U,0,1);
  for (int site = 0; site < Nvol; site++) {
    plaq += ReTr(U->mat(site, 0) * staple.mat_dag(site));   // P_xy
  }

  staple = upper(U, 1, 2);
  // staple = lower(U,1,2);
  for (int site = 0; site < Nvol; site++) {
    plaq += ReTr(U->mat(site, 1) * staple.mat_dag(site));   // P_yz
  }

  staple = upper(U, 2, 0);
  // staple = lower(U,2,0);
  for (int site = 0; site < Nvol; site++) {
    plaq += ReTr(U->mat(site, 2) * staple.mat_dag(site));   // P_zx
  }

  plaq = comm->reduce_sum(plaq);

  return plaq / (Lvol * Nc * 3.0);
}


//====================================================================
double Staples_eo::plaq_t(const Field_G *U)
{
  double         plaq = 0.0;
  static Field_G staple;

  for (int nu = 0; nu < Ndim - 1; nu++) {
    staple = lower(U, 3, nu);
    //    staple = upper(U,3,nu);
    for (int site = 0; site < Nvol; site++) {
      plaq += ReTr(U->mat(site, 3) * staple.mat_dag(site));   // P_zx
    }
  }

  plaq = comm->reduce_sum(plaq);

  return plaq / (Lvol * Nc * 3.0);
}


//====================================================================
void Staples_eo::staple(Field_G& W, const Field_G *U, const int mu)
{
  W = 0.0;
  for (int nu = 0; nu < Ndim; nu++) {
    if (nu != mu) {
      W += upper(U, mu, nu);
      W += lower(U, mu, nu);
    }
  }
}


//====================================================================
Field_G Staples_eo::upper(const Field_G *U, const int mu, const int nu)
{
  // (1)  mu (2)
  //    +-->--+
  // nu |     |
  //   i+     +

  Field_G c;

  Umu.setpart_ex(0, *U, mu);
  Unu.setpart_ex(0, *U, nu);

  shift.backward(v, Unu, mu);
  shift.backward(c, Umu, nu);

  mult_Field_Gnd(w, 0, c, 0, v, 0);
  mult_Field_Gnn(c, 0, Unu, 0, w, 0);

  return c;
}


//====================================================================
Field_G Staples_eo::lower(const Field_G *U, const int mu, const int nu)
{
  //    +     +
  // nu |     |
  //   i+-->--+
  //  (1)  mu (2)

  Field_G c;

  Umu.setpart_ex(0, *U, mu);
  Unu.setpart_ex(0, *U, nu);

  shift.backward(w, Unu, mu);
  mult_Field_Gnn(v, 0, Umu, 0, w, 0);
  //mult_Field_Gnn(v, 0,*U,mu,w,0);
  mult_Field_Gdn(w, 0, Unu, 0, v, 0);
  shift.forward(c, w, nu);

  return c;
}


//====================================================================
//============================================================END=====