fopr_Staggered_eo.cpp

Go to the documentation of this file.

#include "fopr_Staggered_eo.h"

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

using std::valarray;

//- parameter entries
namespace {
  void append_entry(Parameters& param)
  {
    param.Register_double("quark_mass", 0.0);
    param.Register_int_vector("boundary_condition", std::valarray<int>());

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


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

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

//====================================================================
void Fopr_Staggered_eo::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
  double        mq;
  valarray<int> bc;

  int err = 0;
  err += params.fetch_double("quark_mass", mq);
  err += params.fetch_int_vector("boundary_condition", bc);

  if (err) {
    vout.crucial(m_vl, "Fopr_Staggered_eo: fetch error, input parameter not found.\n");
    abort();
  }


  set_parameters(mq, bc);
}


//====================================================================
void Fopr_Staggered_eo::set_parameters(const double mq, const std::valarray<int> bc)
{
  //- print input parameters
  vout.general(m_vl, "Parameters of staggered fermion operator:\n");
  vout.general(m_vl, "  mq   = %8.4f\n", mq);
  for (int mu = 0; mu < m_Ndim; ++mu) {
    vout.general(m_vl, "  boundary[%d] = %2d\n", mu, bc[mu]);
  }

  //- range check
  int err = 0;
  err += ParameterCheck::non_zero(mq);

  if (err) {
    vout.crucial(m_vl, "Fopr_Staggered_eo: parameter range check failed.\n");
    abort();
  }

  assert(bc.size() == m_Ndim);

  //- store values
  m_mq = mq;
  // m_boundary.resize(m_Ndim);  // already resized in the constructor.
  for (int mu = 0; mu < m_Ndim; ++mu) {
    m_boundary[mu] = bc[mu];
  }
}


//====================================================================
void Fopr_Staggered_eo::set_staggered_phase()
{
  int Nt   = CommonParameters::Nt();
  int Nz   = CommonParameters::Nz();
  int Ny   = CommonParameters::Ny();
  int Nx   = CommonParameters::Nx();
  int Nvol = CommonParameters::Nvol();

  Field     phase_lex(1, Nvol, m_Ndim);
  Index_lex idx_lex;

  int ipex = Communicator::ipe(0);
  int ipey = Communicator::ipe(1);
  int ipez = Communicator::ipe(2);
  int ipet = Communicator::ipe(3);

  for (int t = 0; t < Nt; ++t) {
    int t2 = t + ipet * Nt;
    for (int z = 0; z < Nz; ++z) {
      int z2 = z + ipez * Nz;
      for (int y = 0; y < Ny; ++y) {
        int y2 = y + ipey * Ny;
        for (int x = 0; x < Nx; ++x) {
          int x2 = x + ipex * Nx;
          int is = idx_lex.site(x, y, z, t);

          phase_lex.set(0, is, 0, 1.0);
          phase_lex.set(0, is, 1, 1.0);
          phase_lex.set(0, is, 2, 1.0);
          phase_lex.set(0, is, 3, 1.0);

          if (((x2 + 1) % 2) == 1) phase_lex.set(0, is, 1, -1.0);
          if (((x2 + y2 + 2) % 2) == 1) phase_lex.set(0, is, 2, -1.0);
          if (((x2 + y2 + z2 + 3) % 2) == 1) phase_lex.set(0, is, 3, -1.0);
        }
      }
    }
  }

  m_idx.convertField(m_staggered_phase, phase_lex);
}


//====================================================================
void Fopr_Staggered_eo::mult_staggered_phase(Field& v, int mu, int ieo)
{
  int Nvol  = CommonParameters::Nvol();
  int Nvol2 = Nvol / 2;

  int Nin = v.nin();
  int Nex = v.nex();

  assert(v.nvol() == Nvol2);

  for (int ex = 0; ex < Nex; ++ex) {
    for (int site = 0; site < Nvol2; ++site) {
      int site2 = m_idx.site(site, ieo);

      for (int in = 0; in < Nin; ++in) {
        double vt = m_staggered_phase.cmp(0, site2, mu) * v.cmp(in, site, ex);
        v.set(in, site, ex, vt);
      }
    }
  }
}


//====================================================================
void Fopr_Staggered_eo::set_staggered_config(Field_G *Ueo)
{
  int    Nc = CommonParameters::Nc();
  double ur, ui;

  for (int mu = 0; mu < m_Ndim; ++mu) {
    for (int site = 0; site < m_Nvol; ++site) {
      for (int cc = 0; cc < Nc * Nc; ++cc) {
        ur = m_staggered_phase.cmp(0, site, mu) * Ueo->cmp_r(cc, site, mu);
        ui = m_staggered_phase.cmp(0, site, mu) * Ueo->cmp_i(cc, site, mu);

        m_Ueo->set_ri(cc, site, mu, ur, ui);
      }
    }
  }
}


//====================================================================
const Field Fopr_Staggered_eo::Meo(const Field& f, const int ieo)
{
  int NinF = 2 * CommonParameters::Nc();

  assert(f.nin() == NinF);
  assert(f.nvol() == m_Nvol2);
  assert(f.nex() == 1);

  Field_F_1spinor w(m_Nvol2, 1);
  w = 0.0;

  mult_xp(w, (Field_F_1spinor)f, ieo);
  mult_xm(w, (Field_F_1spinor)f, ieo);

  mult_yp(w, (Field_F_1spinor)f, ieo);
  mult_ym(w, (Field_F_1spinor)f, ieo);

  mult_zp(w, (Field_F_1spinor)f, ieo);
  mult_zm(w, (Field_F_1spinor)f, ieo);

  mult_tp(w, (Field_F_1spinor)f, ieo);
  mult_tm(w, (Field_F_1spinor)f, ieo);

  w *= 0.5 / m_mq;

  return (Field)w;
}


//====================================================================
void Fopr_Staggered_eo::mult_xp(Field_F_1spinor& w, const Field_F_1spinor& f,
                                const int ieo)
{
  shift.backward_h(*tr, f, m_boundary[0], 0, ieo);

  for (int ix = 0; ix < m_Nvol2; ++ix) {
    int gx = m_idx.site(ix, ieo);
    w.add_vec(ix, 0, m_Ueo->mat(gx, 0) * tr->vec(ix));
  }
}


//====================================================================
void Fopr_Staggered_eo::mult_xm(Field_F_1spinor& w, const Field_F_1spinor& f,
                                const int ieo)
{
  for (int ix = 0; ix < m_Nvol2; ++ix) {
    int gx = m_idx.site(ix, 1 - ieo);
    tr->set_vec(ix, 0, m_Ueo->mat_dag(gx, 0) * f.vec(ix));
  }

  shift.forward_h(*tr2, *tr, m_boundary[0], 0, ieo);
  w -= *tr2;
}


//====================================================================
void Fopr_Staggered_eo::mult_yp(Field_F_1spinor& w, const Field_F_1spinor& f,
                                const int ieo)
{
  shift.backward_h(*tr, f, m_boundary[1], 1, ieo);

  for (int ix = 0; ix < m_Nvol2; ++ix) {
    int gx = m_idx.site(ix, ieo);
    w.add_vec(ix, 0, m_Ueo->mat(gx, 1) * tr->vec(ix));
  }
}


//====================================================================
void Fopr_Staggered_eo::mult_ym(Field_F_1spinor& w, const Field_F_1spinor& f,
                                const int ieo)
{
  for (int ix = 0; ix < m_Nvol2; ++ix) {
    int gx = m_idx.site(ix, 1 - ieo);
    tr->set_vec(ix, 0, m_Ueo->mat_dag(gx, 1) * f.vec(ix));
  }

  shift.forward_h(*tr2, *tr, m_boundary[1], 1, ieo);

  w -= *tr2;
}


//====================================================================
void Fopr_Staggered_eo::mult_zp(Field_F_1spinor& w, const Field_F_1spinor& f,
                                const int ieo)
{
  shift.backward_h(*tr, f, m_boundary[2], 2, ieo);

  for (int ix = 0; ix < m_Nvol2; ++ix) {
    int gx = m_idx.site(ix, ieo);
    w.add_vec(ix, 0, m_Ueo->mat(gx, 2) * tr->vec(ix));
  }
}


//====================================================================
void Fopr_Staggered_eo::mult_zm(Field_F_1spinor& w, const Field_F_1spinor& f,
                                const int ieo)
{
  for (int ix = 0; ix < m_Nvol2; ++ix) {
    int gx = m_idx.site(ix, 1 - ieo);
    tr->set_vec(ix, 0, m_Ueo->mat_dag(gx, 2) * f.vec(ix));
  }

  shift.forward_h(*tr2, *tr, m_boundary[2], 2, ieo);

  w -= *tr2;
}


//====================================================================
void Fopr_Staggered_eo::mult_tp(Field_F_1spinor& w, const Field_F_1spinor& f,
                                const int ieo)
{
  shift.backward_h(*tr, f, m_boundary[3], 3, ieo);

  for (int ix = 0; ix < m_Nvol2; ++ix) {
    int gx = m_idx.site(ix, ieo);
    w.add_vec(ix, 0, m_Ueo->mat(gx, 3) * tr->vec(ix));
  }
}


//====================================================================
void Fopr_Staggered_eo::mult_tm(Field_F_1spinor& w, const Field_F_1spinor& f,
                                const int ieo)
{
  for (int ix = 0; ix < m_Nvol2; ++ix) {
    int gx = m_idx.site(ix, 1 - ieo);
    tr->set_vec(ix, 0, m_Ueo->mat_dag(gx, 3) * f.vec(ix));
  }

  shift.forward_h(*tr2, *tr, m_boundary[3], 3, ieo);

  w -= *tr2;
}


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