fopr_CloverTerm_eo.cpp

Go to the documentation of this file.

#include "fopr_CloverTerm_eo.h"

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

using std::valarray;

//====================================================================
//- parameter entry
namespace {
  void append_entry(Parameters& param)
  {
    param.Register_double("hopping_parameter", 0.0);
    param.Register_double("clover_coefficient", 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.CloverTerm_eo", append_entry);
#endif
}
//- end

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

//====================================================================
void Fopr_CloverTerm_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        kappa, cSW;
  valarray<int> bc;

  int err = 0;
  err += params.fetch_double("hopping_parameter", kappa);
  err += params.fetch_double("clover_coefficient", cSW);
  err += params.fetch_int_vector("boundary_condition", bc);

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


  set_parameters(kappa, cSW, bc);
}


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

  //- range check
  // NB. kappa,cSW == 0 is allowed.
  assert(bc.size() == m_Ndim);

  //- store values
  m_kappa = kappa;
  m_cSW   = cSW;
  assert(bc.size() == m_Ndim);
  for (int mu = 0; mu < m_Ndim; ++mu) {
    m_boundary[mu] = bc[mu];
  }
}


//====================================================================
void Fopr_CloverTerm_eo::set_config(Field *Ueo)
{
  m_Ueo = (Field_G *)Ueo;
  set_csw();
}


//====================================================================
void Fopr_CloverTerm_eo::init(std::string repr)
{
  m_Ueo = 0;

  m_repr = repr;

  m_boundary.resize(m_Ndim);
  m_GM.resize(m_Ndim + 1);
  m_SG.resize(m_Ndim * m_Ndim);

  GammaMatrixSet *gmset = GammaMatrixSet::New(m_repr);

  m_GM[0] = gmset->get_GM(gmset->GAMMA1);
  m_GM[1] = gmset->get_GM(gmset->GAMMA2);
  m_GM[2] = gmset->get_GM(gmset->GAMMA3);
  m_GM[3] = gmset->get_GM(gmset->GAMMA4);
  m_GM[4] = gmset->get_GM(gmset->GAMMA5);

  m_SG[sg_index(0, 1)] = gmset->get_GM(gmset->SIGMA12);
  m_SG[sg_index(1, 2)] = gmset->get_GM(gmset->SIGMA23);
  m_SG[sg_index(2, 0)] = gmset->get_GM(gmset->SIGMA31);
  m_SG[sg_index(3, 0)] = gmset->get_GM(gmset->SIGMA41);
  m_SG[sg_index(3, 1)] = gmset->get_GM(gmset->SIGMA42);
  m_SG[sg_index(3, 2)] = gmset->get_GM(gmset->SIGMA43);

  m_SG[sg_index(1, 0)] = m_SG[sg_index(0, 1)].mult(-1);
  m_SG[sg_index(2, 1)] = m_SG[sg_index(1, 2)].mult(-1);
  m_SG[sg_index(0, 2)] = m_SG[sg_index(2, 0)].mult(-1);
  m_SG[sg_index(0, 3)] = m_SG[sg_index(3, 0)].mult(-1);
  m_SG[sg_index(1, 3)] = m_SG[sg_index(3, 1)].mult(-1);
  m_SG[sg_index(2, 3)] = m_SG[sg_index(3, 2)].mult(-1);

  m_SG[sg_index(0, 0)] = gmset->get_GM(gmset->UNITY);
  m_SG[sg_index(1, 1)] = gmset->get_GM(gmset->UNITY);
  m_SG[sg_index(2, 2)] = gmset->get_GM(gmset->UNITY);
  m_SG[sg_index(3, 3)] = gmset->get_GM(gmset->UNITY);
  // these 4 gamma matrices are actually not used.

  delete gmset;
}


//====================================================================
std::vector<double> Fopr_CloverTerm_eo::csmatrix(const int& site)
{
  std::vector<double> matrix(m_Nc * m_Nc * m_Nd * m_Nd * 2);

  for (int ispin = 0; ispin < m_Nd / 2; ++ispin) {
    for (int icolor = 0; icolor < m_Nc; ++icolor) {
      int ics = icolor + ispin * m_Nc;
      for (int jspin = 0; jspin < m_Nd; ++jspin) {
        int js2 = (jspin + m_Nd / 2) % m_Nd;
        for (int jcolor = 0; jcolor < m_Nc; ++jcolor) {
          int cs1 = jcolor + m_Nc * (jspin + m_Nd * ics);
          int cs2 = jcolor + m_Nc * (jspin + m_Nd * (ics + m_Nc * m_Nd / 2));
          int cc  = jcolor + icolor * m_Nc;
          int ss1 = jspin + ispin * m_Nd;
          int ss2 = js2 + ispin * m_Nd;

          matrix[2 * cs1]     = m_T.cmp_r(cc, site, ss1);
          matrix[2 * cs1 + 1] = m_T.cmp_i(cc, site, ss1);

          matrix[2 * cs2]     = m_T.cmp_r(cc, site, ss2);
          matrix[2 * cs2 + 1] = m_T.cmp_i(cc, site, ss2);
        }
      }
    }
  }

  return matrix;
}


//====================================================================
void Fopr_CloverTerm_eo::D(Field& v, const Field& f, const int ieo)
{
  // Field_F& vf = (Field_F&) v;
  // Field_F& ff = (Field_F&) f;

  Field_F vf(v.nvol(), v.nex());
  Field_F ff(f.nvol(), f.nex());

  ff = (Field_F)f;

  int n_ex = f.nex();

  Vec_SU_N u_vec, d_vec;
  for (int iex = 0; iex < n_ex; ++iex) {
    for (int isite = 0; isite < m_Nvol2; ++isite) {
      for (int ispin = 0; ispin < m_Nd / 2; ++ispin) {
        u_vec.zero();
        d_vec.zero();
        for (int jspin = 0; jspin < m_Nd; ++jspin) {
          int u_spin = jspin + ispin * m_Nd;
          u_vec += m_T.mat(idx.site(isite, ieo), u_spin)
                   * ff.vec(jspin, isite, iex);
          int d_spin = (jspin + m_Nd / 2) % m_Nd + ispin * m_Nd;
          d_vec += m_T.mat(idx.site(isite, ieo), d_spin)
                   * ff.vec(jspin, isite, iex);
        }
        vf.set_vec(ispin, isite, iex, u_vec);
        vf.set_vec(ispin + m_Nd / 2, isite, iex, d_vec);
      }
    }
  }

  v = (Field)vf;
}


//====================================================================
void Fopr_CloverTerm_eo::mult_isigma(Field_F& v, const Field_F& w,
                                     const int mu, const int nu)
{
  assert(mu != nu);
  v.mult_iGM(m_SG[sg_index(mu, nu)], w);
}


//====================================================================
void Fopr_CloverTerm_eo::set_csw()
{
  m_T = 0.0;

  //- sigma23
  Field_G F;
  set_fieldstrength(F, 1, 2);
  F.xI();
  m_T.addpart_ex(1, F, 0);
  m_T.addpart_ex(4, F, 0);

  //- sigma31
  set_fieldstrength(F, 2, 0);
  m_T.addpart_ex(1, F, 0);
  m_T.addpart_ex(4, -F, 0);

  //- sigma12
  set_fieldstrength(F, 0, 1);
  F.xI();
  m_T.addpart_ex(0, F, 0);
  m_T.addpart_ex(5, -F, 0);

  //- sigma41
  set_fieldstrength(F, 3, 0);
  F.xI();
  F *= -1.0;
  m_T.addpart_ex(3, F, 0);
  m_T.addpart_ex(6, F, 0);

  //- sigma42
  set_fieldstrength(F, 3, 1);
  m_T.addpart_ex(6, F, 0);
  m_T.addpart_ex(3, -F, 0);

  //- sigma43
  set_fieldstrength(F, 3, 2);
  F.xI();
  m_T.addpart_ex(7, F, 0);
  m_T.addpart_ex(2, -F, 0);

  m_T *= -m_kappa * m_cSW;

  //- add unit color matrix for diag part of dirac matrix;
  for (int ispin = 0; ispin < m_Nd / 2; ++ispin) {
    int spin_diag = ispin + m_Nd * ispin;
    for (int isite = 0; isite < m_Nvol; ++isite) {
      for (int icolor = 0; icolor < m_Nc; ++icolor) {
        int cc2 = 2 * (icolor * m_Nc + icolor);
        m_T.add(cc2, isite, spin_diag, 1.0);
      }
    }
  }
}


//====================================================================
void Fopr_CloverTerm_eo::set_fieldstrength(Field_G& Fst,
                                           const int mu, const int nu)
{
  Staples_eo staple;

  Field_G Cup(m_Nvol, 1), Cdn(m_Nvol, 1);
  Field_G Umu(m_Nvol, 1);
  Field_G w(m_Nvol, 1), v(m_Nvol, 1), v2(m_Nvol, 1);

  Cup = staple.upper(m_Ueo, mu, nu);
  Cdn = staple.lower(m_Ueo, mu, nu);
  Umu.setpart_ex(0, *m_Ueo, mu);

  for (int site = 0; site < m_Nvol; ++site) {
    w.set_mat(site, 0, Umu.mat(site) * Cup.mat_dag(site));
  }

  for (int site = 0; site < m_Nvol; ++site) {
    v2.set_mat(site, 0, Umu.mat(site) * Cdn.mat_dag(site));
  }

  w -= v2;

  for (int site = 0; site < m_Nvol; ++site) {
    v.set_mat(site, 0, Cup.mat_dag(site) * Umu.mat(site));
  }

  for (int site = 0; site < m_Nvol; ++site) {
    v2.set_mat(site, 0, Cdn.mat_dag(site) * Umu.mat(site));
  }

  v -= v2;

  m_shift_eo.forward(v2, v, mu);

  w += v2;

  for (int site = 0; site < m_Nvol; ++site) {
    Fst.set_mat(site, 0, w.mat(site).ah());
  }

  Fst *= 0.25;
}


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