force_F_Domainwall.cpp

Go to the documentation of this file.

#include "force_F_Domainwall.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_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("Force.F_Domainwall", append_entry);
#endif
}
//- end

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

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

  int err = 0;
  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, "Force_F_Domainwall: fetch error, input parameter not found.\n");
    abort();
  }


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


//====================================================================
void Force_F_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 force:\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, "Force_F_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];
  }

  //- post-process

  //- Domain-wall operator
  m_fopr_dw->set_parameters(m_mq, m_M0, m_Ns, m_boundary);

  //- Wilson opr/force
  double kappa = 1.0 / (8.0 - 2.0 * m_M0);
  m_fopr_w->set_parameters(kappa, m_boundary);
  m_force_w->set_parameters(kappa, m_boundary);
}


//====================================================================
Field Force_F_Domainwall::force_core(const Field& eta)
{
  int Nc   = CommonParameters::Nc();
  int Nvol = CommonParameters::Nvol();
  int Ndim = CommonParameters::Ndim();
  int NinG = 2 * Nc * Nc;

  Field_G  force(Nvol, Ndim), force1(Nvol, Ndim);
  Mat_SU_N ut(Nc);

  force1 = force_udiv(eta);

  for (int mu = 0; mu < Ndim; ++mu) {
    for (int site = 0; site < Nvol; ++site) {
      ut = m_U->mat(site, mu) * force1.mat(site, mu);
      ut.at();
      ut *= -2.0;
      force.set_mat(site, mu, ut);
    }
  }

  return (Field)force;
}


//====================================================================
Field Force_F_Domainwall::force_udiv(const Field& eta)
{
  int Nc   = CommonParameters::Nc();
  int Nd   = CommonParameters::Nd();
  int Nvol = CommonParameters::Nvol();
  int Ndim = CommonParameters::Ndim();
  int NinG = 2 * Nc * Nc;
  int NinF = 2 * Nc * Nd;

  Field force(NinG, Nvol, Ndim), force2(NinG, Nvol, Ndim);
  Field zeta(NinF, Nvol, m_Ns);

  zeta = m_fopr_dw->H(eta);

  force  = force_udiv1(zeta, eta);
  force2 = force_udiv1(eta, zeta);
  force += force2;

  return (Field)force;
}


//====================================================================
Field Force_F_Domainwall::force_udiv1(const Field& zeta,
                                      const Field& eta)
{
  int Nc   = CommonParameters::Nc();
  int Nd   = CommonParameters::Nd();
  int Nvol = CommonParameters::Nvol();
  int Ndim = CommonParameters::Ndim();
  int NinG = 2 * Nc * Nc;
  int NinF = 2 * Nc * Nd;

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

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

  Field eta4(NinF, Nvol, 1), zeta4(NinF, Nvol, 1);

  Field force(NinG, Nvol, Ndim);

  force = 0.0;

  for (int is = 0; is < m_Ns; ++is) {
    zeta4.setpart_ex(0, zeta, m_Ns - 1 - is);
    eta4.setpart_ex(0, eta, is);
    m_fopr_w->fopr_normalize(eta4);

    force += m_force_w->force_udiv1((Field_F)zeta4, (Field_F)eta4);
  }

  return force;
}


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