fopr_TMWilson.cpp

Go to the documentation of this file.

#include "fopr_TMWilson.h"

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

using std::valarray;

//- parameter entries
namespace {
  void append_entry(Parameters& param)
  {
    param.Register_string("gamma_matrix_type", "NULL");

    param.Register_double("hopping_parameter", 0.0);
    param.Register_double("twisted_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.TMWilson", append_entry);
#endif
}
//- end

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

//====================================================================
void Fopr_TMWilson::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, tw_mass;
  valarray<int> bc;

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

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


  set_parameters(kappa, tw_mass, bc);
}


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

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

  //- store values
  m_kappa   = kappa;
  m_tw_mass = tw_mass;

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

  //- propagate parameters
  m_fopr_w->set_parameters(m_kappa, m_boundary);
}


//====================================================================
void Fopr_TMWilson::D(Field& v, const Field& f)
{
  Field_F w2(f.nvol(), f.nex());

  m_fopr_w->D(v, f);
  m_fopr_w->mult_gm5(w2, f);
  w2.xI();
  double tm_2K = 2.0 * m_kappa * m_tw_mass;
  v += tm_2K * (Field)w2;
}


//====================================================================
void Fopr_TMWilson::Ddag(Field& v, const Field& f)
{
  Field_F w2(f.nvol(), f.nex());

  m_fopr_w->Ddag(v, f);
  m_fopr_w->mult_gm5(w2, f);
  w2.xI();
  double tm_2K = 2.0 * m_kappa * m_tw_mass;
  v -= tm_2K * (Field)w2;
}


//====================================================================
void Fopr_TMWilson::H(Field& v, const Field& f)
{
  Field_F w2(f.nvol(), f.nex());

  m_fopr_w->H(v, f);
  w2.Ix(f);
  double tm_2K = 2.0 * m_kappa * m_tw_mass;
  v += tm_2K * (Field)w2;
}


//====================================================================
void Fopr_TMWilson::Hdag(Field& v, const Field& f)
{
  Field_F w2(f.nvol(), f.nex());

  m_fopr_w->H(v, f);
  w2.Ix(f);
  double tm_2K = 2.0 * m_kappa * m_tw_mass;
  v -= tm_2K * (Field)w2;
}


//====================================================================
void Fopr_TMWilson::DdagD(Field& v, const Field& f)
{
  m_fopr_w->DdagD(v, f);
  double tm_2K = 2.0 * m_kappa * m_tw_mass;
  v += (tm_2K * tm_2K) * f;
}


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