fopr_Domainwall.cpp

Go to the documentation of this file.

#include "fopr_Domainwall.h"

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

using std::valarray;
using std::string;

//- parameter entries
namespace {
  void append_entry(Parameters& param)
  {
    param.Register_string("gamma_matrix_type", "NULL");
    param.Register_double("quark_mass", 0.0);
    param.Register_double("domain_wall_height", 0.0);
    param.Register_int("extent_of_5th_dimension", 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.Domainwall", append_entry);
#endif
}
//- end

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

//====================================================================
void Fopr_Domainwall::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
  string        str_gmset_type;
  double        mq, M0;
  int           Ns;
  valarray<int> bc;

  int err          = 0;
  int err_optional = 0;
  err_optional += params.fetch_string("gamma_matrix_type", str_gmset_type);
  err          += params.fetch_double("quark_mass", mq);
  err          += params.fetch_double("domain_wall_height", M0);
  err          += params.fetch_int("extent_of_5th_dimension", Ns);
  err          += params.fetch_int_vector("boundary_condition", bc);

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


  set_parameters(mq, M0, Ns, bc);
}


//====================================================================
void Fopr_Domainwall::set_parameters(const double mq, const double M0,
                                     const int Ns, const valarray<int> bc)
{
  int Ndim = CommonParameters::Ndim();

  //- print input parameters
  vout.general(m_vl, "Parameters of domain-wall fermion operator:\n");
  vout.general(m_vl, "  mq   = %8.4f\n", mq);
  vout.general(m_vl, "  M0   = %8.4f\n", M0);
  vout.general(m_vl, "  Ns   = %4d\n", Ns);
  for (int mu = 0; mu < Ndim; ++mu) {
    vout.general(m_vl, "  boundary[%d] = %2d\n", mu, bc[mu]);
  }

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

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

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

  //- store values
  m_mq = mq;
  m_M0 = M0;
  m_Ns = Ns;

  m_boundary.resize(Ndim);
  for (int mu = 0; mu < Ndim; ++mu) {
    m_boundary[mu] = bc[mu];
  }

  //- propagate parameters
  double kappa = 1.0 / (8.0 - 2.0 * m_M0);
  m_fopr_w->set_parameters(kappa, m_boundary);
}


//====================================================================
const Field Fopr_Domainwall::DdagD(const Field& w)
{
  return H(H(w));
}


//====================================================================
const Field Fopr_Domainwall::mult_chproj_4d(const Field& w, int ipm)
{
  Field v(w);

  if (ipm == 1) {
    v += m_fopr_w->mult_gm5(w);
  } else if (ipm == -1) {
    v -= m_fopr_w->mult_gm5(w);
  } else {
    vout.crucial(m_vl, "Fopr_Domainwall: illegal chirality = %d.\n", ipm);
    abort();
  }

  v *= 0.5;

  return v;
}


//====================================================================
const Field Fopr_Domainwall::Ddag(const Field& w)
{
  return H(mult_gm5R(w));
}


//====================================================================
const Field Fopr_Domainwall::H(const Field& w)
{
  return mult_gm5R(D(w));
}


//====================================================================
const Field Fopr_Domainwall::mult_gm5R(const Field& w)
{
  int Nc   = CommonParameters::Nc();
  int Nd   = CommonParameters::Nd();
  int Nvol = CommonParameters::Nvol();
  int NinF = 2 * Nc * Nd;

  assert(w.nin() == NinF);
  assert(w.nvol() == Nvol);
  assert(w.nex() == m_Ns);

  Field v(NinF, Nvol, m_Ns);
  Field w4(NinF, Nvol, 1), v4(NinF, Nvol, 1);

  for (int is = 0; is < m_Ns; ++is) {
    w4.setpart_ex(0, w, is);
    v4 = m_fopr_w->mult_gm5(w4);
    v.setpart_ex(m_Ns - 1 - is, v4, 0);
  }

  return v;
}


//====================================================================
const Field Fopr_Domainwall::D(const Field& w)
{
  int Nc   = CommonParameters::Nc();
  int Nd   = CommonParameters::Nd();
  int Nvol = CommonParameters::Nvol();
  int NinF = 2 * Nc * Nd;

  assert(w.nin() == NinF);
  assert(w.nvol() == Nvol);
  assert(w.nex() == m_Ns);

  Field v(NinF, Nvol, m_Ns);
  Field w4(NinF, Nvol, 1), w4_up(NinF, Nvol, 1), w4_dn(NinF, Nvol, 1);
  Field v4(NinF, Nvol, 1);

  for (int is = 0; is < m_Ns; ++is) {
    w4.setpart_ex(0, w, is);
    m_fopr_w->set_mode("D");
    v4 = m_fopr_w->mult(w4);
    m_fopr_w->fopr_normalize(v4);

    int    is_up = (is + 1) % m_Ns;
    double Fup   = -0.5;
    if (is == m_Ns - 1) Fup = m_mq / 2.0;
    w4_up.setpart_ex(0, w, is_up);
    w4_up -= m_fopr_w->mult_gm5(w4_up);
    w4.addpart_ex(0, w4_up, 0, Fup);

    int    is_dn = (is - 1 + m_Ns) % m_Ns;
    double Fdn   = -0.5;
    if (is == 0) Fdn = m_mq / 2.0;
    w4_dn.setpart_ex(0, w, is_dn);
    w4_dn += m_fopr_w->mult_gm5(w4_dn);
    w4.addpart_ex(0, w4_dn, 0, Fdn);

    v.setpart_ex(is, v4, 0);
    v.addpart_ex(is, w4, 0);
  }

  return v;
}


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