staples_SF.cpp

Go to the documentation of this file.

#include "staples_SF.h"

#ifdef USE_PARAMETERS_FACTORY
#include "parameters_factory.h"
#endif


//- parameter entries
namespace {
  void append_entry(Parameters& param)
  {
    param.Register_double_vector("phi", std::vector<double>());
    param.Register_double_vector("phipr", std::vector<double>());
    param.Register_double_vector("pomega", std::vector<double>());

    param.Register_string("verbose_level", "NULL");
  }


#ifdef USE_PARAMETERS_FACTORY
  bool init_param = ParametersFactory::Register("Staples_SF", append_entry);
#endif
}
//- end

//- parameters class
Parameters_Staples_SF::Parameters_Staples_SF() { append_entry(*this); }
//- end

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

//********************************************************************
void Staples_SF::set_parameters(const Parameters& params)
{
  const string str_vlevel = params.get_string("verbose_level");

  m_vl = vout.set_verbose_level(str_vlevel);

  //- fetch and check input parameters
  std::vector<double> phi, phipr, pomega;

  int err = 0;
  err += params.fetch_double_vector("phi", phi);
  err += params.fetch_double_vector("phipr", phipr);
  err += params.fetch_double_vector("pomega", pomega);

  if (err) {
    vout.crucial(m_vl, "%s: fetch error, input parameter not found.\n", class_name.c_str());
    exit(EXIT_FAILURE);
  }


  set_parameters(phi, phipr, pomega);  // call std::vector version
}


//====================================================================
void Staples_SF::set_parameters(std::vector<double>& phi, std::vector<double>& phipr,
                                std::vector<double>& pomega)
{
  //- print input parameters
  vout.general(m_vl, "Parameters of %s:\n", class_name.c_str());
  vout.general(m_vl, "  phi1    = %12.6f\n", phi[0]);
  vout.general(m_vl, "  phi2    = %12.6f\n", phi[1]);
  vout.general(m_vl, "  phi3    = %12.6f\n", phi[2]);
  vout.general(m_vl, "  phipr1  = %12.6f\n", phipr[0]);
  vout.general(m_vl, "  phipr2  = %12.6f\n", phipr[1]);
  vout.general(m_vl, "  phipr3  = %12.6f\n", phipr[2]);
  vout.general(m_vl, "  pomega1 = %12.6f\n", pomega[0]);
  vout.general(m_vl, "  pomega2 = %12.6f\n", pomega[1]);
  vout.general(m_vl, "  pomega3 = %12.6f\n", pomega[2]);

  //- range check
  // NB. phi,phipr,pomega == 0 is allowed.
  assert(phi.size() == 3);
  assert(phipr.size() == 3);
  assert(pomega.size() == 3);

  //- store values
  set_parameters(&phi[0], &phipr[0], &pomega[0]);  // call double[] version
}


//====================================================================
void Staples_SF::set_parameters(double *phi, double *phipr)
{
  double c0r, c0i, c1r, c1i, c2r, c2i;

  c0r = cos(phi[0] / Lx);
  c0i = sin(phi[0] / Lx);
  c1r = cos(phi[1] / Lx);
  c1i = sin(phi[1] / Lx);
  c2r = cos(phi[2] / Lx);
  c2i = sin(phi[2] / Lx);

  wk.zero();
  wk.set(0, 0, c0r, c0i);
  wk.set(1, 1, c1r, c1i);
  wk.set(2, 2, c2r, c2i);

  c0r = cos(phipr[0] / Lx);
  c0i = sin(phipr[0] / Lx);
  c1r = cos(phipr[1] / Lx);
  c1i = sin(phipr[1] / Lx);
  c2r = cos(phipr[2] / Lx);
  c2i = sin(phipr[2] / Lx);

  wkpr.zero();
  wkpr.set(0, 0, c0r, c0i);
  wkpr.set(1, 1, c1r, c1i);
  wkpr.set(2, 2, c2r, c2i);

  iomega0.zero();
  iomega0.set(0, 0, 0.0, 1.0);
  iomega0.set(1, 1, 0.0, -0.5);
  iomega0.set(2, 2, 0.0, -0.5);

  initialized = 1;
}


//====================================================================

void Staples_SF::set_parameters(const double *phi, const double *phipr, const double *pomega)
{
  double c0r, c0i, c1r, c1i, c2r, c2i;

  c0r = cos(phi[0] / Lx);
  c0i = sin(phi[0] / Lx);
  c1r = cos(phi[1] / Lx);
  c1i = sin(phi[1] / Lx);
  c2r = cos(phi[2] / Lx);
  c2i = sin(phi[2] / Lx);
  wk.zero();
  wk.set(0, 0, c0r, c0i);
  wk.set(1, 1, c1r, c1i);
  wk.set(2, 2, c2r, c2i);

  c0r = cos(phipr[0] / Lx);
  c0i = sin(phipr[0] / Lx);
  c1r = cos(phipr[1] / Lx);
  c1i = sin(phipr[1] / Lx);
  c2r = cos(phipr[2] / Lx);
  c2i = sin(phipr[2] / Lx);
  wkpr.zero();
  wkpr.set(0, 0, c0r, c0i);
  wkpr.set(1, 1, c1r, c1i);
  wkpr.set(2, 2, c2r, c2i);

  iomega0.zero();
  iomega0.set(0, 0, 0.0, pomega[0]);
  iomega0.set(1, 1, 0.0, pomega[1]);
  iomega0.set(2, 2, 0.0, pomega[2]);

  initialized = 1;
}


//====================================================================

double Staples_SF::sf_coupling_plaq(const Field_G& U, double ct)
{
  if (!initialized) {
    vout.crucial(m_vl, "%s: Paramer is not initialized.\n", class_name.c_str());
    exit(EXIT_FAILURE);
  }

  double            plaq   = 0.0;
  double            plaqt0 = 0.0;
  double            plaqtT = 0.0;
  double            scr;
  Mat_SU_N          up(Nc);
  static Field_G_SF staple;

  int site;
  for (int nu = 0; nu < Ndim - 1; nu++) {
    upper(staple, U, 3, nu);

    for (int z = 0; z < Nz; z++) {
      for (int y = 0; y < Ny; y++) {
        for (int x = 0; x < Nx; x++) {
          int t;
          // boundary
          if (Communicator::ipe(3) == 0) {
            t    = 0;
            site = index.site(x, y, z, t);
            up   = staple.mat(site) * U.mat_dag(site, 3);
            scr  = ReTr(iomega0 * up);

            /*
            up.unit();
            up *= staple.mat(site);
            up *= U->mat_dag(site,3);
            up *= iomega0;
            scr = ReTr( up );
            */
            plaq   -= scr;
            plaqt0 += scr;
          }
          // boundary
          if (Communicator::ipe(3) == NPEt - 1) {
            t    = Nt - 1;
            site = index.site(x, y, z, t);
            up   = staple.mat_dag(site) * U.mat(site, 3);
            scr  = ReTr(iomega0 * up);

            /*
            up.unit();
            up *= staple.mat_dag(site);
            up *= U->mat(site,3);
            up *= iomega0;
            scr = ReTr( up );
            */
            plaq   += scr;
            plaqtT += scr;
          }
        }
      }
    }
  }

  plaq    = Communicator::reduce_sum(plaq);
  plaqt0  = Communicator::reduce_sum(plaqt0);
  plaqtT  = Communicator::reduce_sum(plaqtT);
  plaq   *= ct / Lx;
  plaqt0 *= ct / Lx;
  plaqtT *= ct / Lx;

  vout.general(m_vl, "SF_delSg_plaq, from 0, from T = %.8f %.8f %.8f\n", plaq, plaqt0, plaqtT);

  return plaq;
}


//====================================================================

double Staples_SF::sf_coupling_rect(const Field_G& m_U, double ctr)
{
  if (!initialized) {
    vout.crucial(m_vl, "%s: Paramer is not initialized.\n", class_name.c_str());
    exit(EXIT_FAILURE);
  }

  double rect;
  double rect01 = 0.0;
  double rect02 = 0.0;
  double rect03 = 0.0;
  double rectt1 = 0.0;
  double rectt2 = 0.0;
  double rectt3 = 0.0;
  // double scr;

  Field_G_SF Cup1(Nvol, 1), Cup2(Nvol, 1);
  Field_G_SF Cdn1(Nvol, 1), Cdn2(Nvol, 1);
  Field_G    Umu(Nvol, 1), Unu(Nvol, 1);
  Field_G    v(Nvol, 1), c(Nvol, 1);

  Mat_SU_N wmat(Nc), cmat(Nc);

  int site;
  int nu = 3;
  for (int mu = 0; mu < Ndim - 1; mu++) {
    // rect01
    //      <---<---+
    //      |       |
    //  t=0 x--->---+
    //  omega0

    // rect02
    //      <---<---+
    //      |       |
    //  t=0 +---x---+
    //       omega0

    // rectt1
    // omega0
    // t=Nt x--->---+
    //      |       |
    //      +---<---+

    // rectt2
    //       omega0
    // t=Nt +---x---+
    //      |       |
    //      <---<---+
    upper(Cup2, m_U, nu, mu);

    Umu.setpart_ex(0, m_U, mu);
    Unu.setpart_ex(0, m_U, nu);
    shift.backward(v, Cup2, mu);
    shift.backward(c, Umu, nu);
    for (int z = 0; z < Nz; z++) {
      for (int y = 0; y < Ny; y++) {
        for (int x = 0; x < Nx; x++) {
          int t;
          t = 0;
          if (Communicator::ipe(3) == 0) {
            site    = index.site(x, y, z, t);
            wmat    = wk * v.mat(site);
            cmat    = wmat * c.mat_dag(site);
            wmat    = cmat * Unu.mat_dag(site);
            rect01 += ReTr(iomega0 * wmat);
            wmat    = iomega0 * v.mat(site);
            cmat    = wmat * c.mat_dag(site);
            wmat    = cmat * Unu.mat_dag(site);
            rect02 += ReTr(wk * wmat);
          }
          t = Nt - 1;
          if (Communicator::ipe(3) == NPEt - 1) {
            site    = index.site(x, y, z, t);
            cmat    = iomega0 * wkpr;
            wmat    = cmat * v.mat_dag(site);
            cmat    = Unu.mat(site) * wmat;
            rectt1 += ReTr(cmat * m_U.mat_dag(site, mu));
            cmat    = iomega0 * v.mat_dag(site);
            wmat    = wkpr * cmat;
            cmat    = Unu.mat(site) * wmat;
            rectt2 += ReTr(cmat * m_U.mat_dag(site, mu));
          }
        }
      }
    }

    // rect03
    //      +---+
    //      |   |
    //      v   ^
    //      |   |
    //  t=0 x--->
    //   omega0
    upper(Cup1, m_U, mu, nu);

    shift.backward(v, Unu, mu);
    shift.backward(c, Cup1, nu);
    for (int z = 0; z < Nz; z++) {
      for (int y = 0; y < Ny; y++) {
        for (int x = 0; x < Nx; x++) {
          int t;
          t = 0;
          if (Communicator::ipe(3) == 0) {
            site    = index.site(x, y, z, t);
            wmat    = c.mat(site) * v.mat_dag(site);
            cmat    = Unu.mat(site) * wmat;
            wmat    = wk * cmat.dag();
            rect03 += ReTr(iomega0 * wmat);
          }
        }
      }
    }

    // rectt3
    //   omega0
    // t=Nt x--->
    //      |   |
    //      ^   v
    //      |   |
    //      +---+
    lower(Cdn1, m_U, mu, nu);

    shift.backward(v, Unu, mu);
    for (int z = 0; z < Nz; z++) {
      for (int y = 0; y < Ny; y++) {
        for (int x = 0; x < Nx; x++) {
          int t;
          t    = Nt - 1;
          site = index.site(x, y, z, t);
          if (Communicator::ipe(3) == NPEt - 1) {
            wmat    = iomega0 * wkpr;
            cmat    = wmat * v.mat_dag(site);
            wmat    = cmat * Cdn1.mat_dag(site);
            rectt3 += ReTr(Unu.mat(site) * wmat);
          }
        }
      }
    }
  }

  rect01  = Communicator::reduce_sum(rect01);
  rect02  = Communicator::reduce_sum(rect02);
  rect03  = Communicator::reduce_sum(rect03);
  rectt1  = Communicator::reduce_sum(rectt1);
  rectt2  = Communicator::reduce_sum(rectt2);
  rectt3  = Communicator::reduce_sum(rectt3);
  rect01 *= ctr / Lx;
  rect02 *= ctr / Lx;
  rect03 /= Lx;
  rectt1 *= ctr / Lx;
  rectt2 *= ctr / Lx;
  rectt3 /= Lx;
  rect    = -rect01 - rect02 - rect03 + rectt1 + rectt2 + rectt3;

  vout.general(m_vl, "SF_delSg_rect, at 01, 02, 03, at T1, T2, T3 = %.8f %.8f %.8f %.8f %.8f %.8f %.8f\n",
               rect, rect01, rect02, rect03, rectt1, rectt2, rectt3);

  return rect;
}


//====================================================================

double Staples_SF::plaquette(const Field_G& U)
{
  if (!initialized) {
    vout.crucial(m_vl, "%s: Paramer is not initialized.\n", class_name.c_str());
    exit(EXIT_FAILURE);
  }
  //  return (plaq_s(U) +plaq_t(U))/2;
  return(plaq_s(U) + plaq_t(U));
}


//====================================================================

double Staples_SF::plaquette_ct(const Field_G& U, double ct)
{
  if (!initialized) {
    vout.crucial(m_vl, "%s: Paramer is not initialized.\n", class_name.c_str());
    exit(EXIT_FAILURE);
  }
  return(plaq_s(U) + plaq_t_ct(U, ct));
}


//====================================================================

double Staples_SF::plaq_s(const Field_G& U)
{
  double            plaq = 0.0;
  static Field_G_SF staple;

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

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

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

  plaq = Communicator::reduce_sum(plaq);

  return plaq;
}


//====================================================================

double Staples_SF::plaq_t(const Field_G& U)
{
  if (!initialized) {
    vout.crucial(m_vl, "%s: Paramer is not initialized.\n", class_name.c_str());
    exit(EXIT_FAILURE);
  }

  double            plaq = 0.0;
  static Field_G_SF staple;

  for (int nu = 0; nu < Ndim - 1; nu++) {
    lower(staple, 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_tk
    }
  }

  plaq = Communicator::reduce_sum(plaq);

  return plaq;
}


//====================================================================

double Staples_SF::plaq_t_ct(const Field_G& U, double ct)
{
  if (!initialized) {
    vout.crucial(m_vl, "%s: Paramer is not initialized.\n", class_name.c_str());
    exit(EXIT_FAILURE);
  }

  double            plaq = 0.0;
  static Field_G_SF staple;

  for (int nu = 0; nu < Ndim - 1; nu++) {
    lower(staple, U, 3, nu);
    // If the node is at the boundary the temporal plaquette is multiplied with ct.
    if (Communicator::ipe(3) == 0) {
      staple.mult_ct_boundary(0, ct);
    }
    if (Communicator::ipe(3) == NPEt - 1) {
      staple.mult_ct_boundary(Nt - 1, ct);
    }
    for (int site = 0; site < Nvol; site++) {
      plaq += ReTr(U.mat(site, 3) * staple.mat_dag(site));   // P_tk
    }
  }

  plaq = Communicator::reduce_sum(plaq);

  return plaq;
}


//====================================================================

void Staples_SF::staple(Field_G& W, const Field_G& U, const int mu)
{
  if (!initialized) {
    vout.crucial(m_vl, "%s: Paramer is not initialized.\n", class_name.c_str());
    exit(EXIT_FAILURE);
  }

  W.set(0.0);
  Field_G_SF c_tmp;
  for (int nu = 0; nu < Ndim; nu++) {
    if (nu != mu) {
      upper(c_tmp, U, mu, nu);
      axpy(W, 1.0, c_tmp);
      lower(c_tmp, U, mu, nu);
      axpy(W, 1.0, c_tmp);
    }
  }
}


//====================================================================

void Staples_SF::staple_ct(Field_G& W, const Field_G& U, const int mu, double ct)
{
  W.set(0.0);
  Field_G_SF staple_upper;
  Field_G_SF staple_lower;

  for (int nu = 0; nu < Ndim; nu++) {
    if (nu != mu) {
      upper(staple_upper, U, mu, nu);
      lower(staple_lower, U, mu, nu);
      if (Communicator::ipe(3) == 0) {
        if (mu == 3) {
          staple_upper.mult_ct_boundary(0, ct);
          staple_lower.mult_ct_boundary(0, ct);
        }
        if (nu == 3) {
          staple_lower.mult_ct_boundary(1, ct);
        }
      }
      if (Communicator::ipe(3) == NPEt - 1) {
        if (mu == 3) {
          staple_upper.mult_ct_boundary(Nt - 1, ct);
          staple_lower.mult_ct_boundary(Nt - 1, ct);
        }
        if (nu == 3) {
          staple_upper.mult_ct_boundary(Nt - 1, ct);
        }
      }
      axpy(W, 1.0, staple_upper);
      axpy(W, 1.0, staple_lower);
    }
  }
}


//====================================================================

void Staples_SF::upper(Field_G_SF& c, const Field_G& U, const int mu, const int nu)
{
  if (!initialized) {
    vout.crucial(m_vl, "%s: Paramer is not initialized.\n", class_name.c_str());
    exit(EXIT_FAILURE);
  }
  // (1)  mu (2)
  //    +-->--+
  // nu |     |
  //   i+     +

  Umu.setpart_ex(0, U, mu);
  Unu.setpart_ex(0, U, nu);
  if (mu != 3) Umu.set_boundary_wk(wk);
  if (nu != 3) Unu.set_boundary_wk(wk);

  shift.backward(v, Unu, mu);
  shift.backward(c, Umu, nu);
  if (mu == 3) v.set_boundary_wkpr(wkpr);
  if (nu == 3) c.set_boundary_wkpr(wkpr);

  mult_Field_Gnd(w, 0, c, 0, v, 0);
  mult_Field_Gnn(c, 0, Unu, 0, w, 0);
  if (mu != 3) c.set_boundary_zero();
}


//====================================================================

void Staples_SF::lower(Field_G_SF& c, const Field_G& U, const int mu, const int nu)
{
  if (!initialized) {
    vout.crucial(m_vl, "%s: Paramer is not initialized.\n", class_name.c_str());
    exit(EXIT_FAILURE);
  }
  //    +     +
  // nu |     |
  //   i+-->--+
  //  (1)  mu (2)

  Umu.setpart_ex(0, U, mu);
  Unu.setpart_ex(0, U, nu);
  if (mu != 3) Umu.set_boundary_wk(wk);
  if (nu != 3) Unu.set_boundary_wk(wk);

  shift.backward(w, Unu, mu);
  if (mu == 3) w.set_boundary_wkpr(wkpr);

  mult_Field_Gnn(v, 0, Umu, 0, w, 0);
  mult_Field_Gdn(w, 0, Unu, 0, v, 0);

  shift.forward(c, w, nu);
  if (mu != 3) c.set_boundary_zero();
}


//====================================================================

void Staples_SF::print_plaquette(const Field_G& U)
{
  double plaq  = plaquette(U);
  double plaq2 = plaq + 3 * 3 * Lx * Ly * Lz;

  vout.general(m_vl, "plaq_SF without boundary spatial plaq = %.8f\n",
               plaq / (3 * Lx * Ly * Lz * (6 * Lt - 3)));
  vout.general(m_vl, "plaq_SF with boundary spatial plaq = %.8f\n",
               plaq2 / (3 * 6 * Lx * Ly * Lz * Lt));
}


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