hmc_General.cpp

Go to the documentation of this file.

#include "hmc_General.h"

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

//- parameter entries
namespace {
  void append_entry(Parameters& param)
  {
    param.Register_int("Metropolis_test", 0);

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


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

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

const std::string HMC_General::class_name = "HMC_General";

//====================================================================
void HMC_General::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
  int Mtpl_test;

  int err = 0;
  err += params.fetch_int("Metropolis_test", Mtpl_test);

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

  set_parameters(Mtpl_test);
}


//====================================================================
void HMC_General::set_parameters(int Mtpl_test)
{
  //- print input parameters
  vout.general(m_vl, "Parameters of %s:\n", class_name.c_str());
  vout.general(m_vl, "  Number of actions: %4d\n", m_action.size());
  vout.general(m_vl, "  Mtlp_test = %4d\n", Mtpl_test);

  //- range check
  // NB. Mtpl_test == 0 is allowed.

  //- store values
  m_Mtpl_test = Mtpl_test;
}


//====================================================================
double HMC_General::update(Field_G& conf_org)
{
  int Nc   = CommonParameters::Nc();
  int Lvol = CommonParameters::Lvol();

  Field_G U(conf_org);

  for (int i = 0; i < m_action.size(); ++i) {
    m_action[i]->set_config(&U);
  }

  int Nin  = U.nin();
  int Nvol = U.nvol();
  int Nex  = U.nex();

  Field_G iP(Nvol, Nex);


  vout.general(m_vl, "HMC (general) start.\n");

  // Langevin step
  double H_tot0 = langevin(iP, U);
  vout.general(m_vl, "  H_total(init)  = %20.10f\n", H_tot0);

  double plaq0 = m_staple->plaquette(U);
  vout.general(m_vl, "  plaq(init)     = %20.12f\n", plaq0);

  // initial Hamiltonian
  //    H_tot0 = calc_Hamiltonian(iP,U);

  // molecular dynamical integration
  m_integrator->evolve(iP, U);

  // trial Hamiltonian
  double H_tot1 = calc_Hamiltonian(iP, U);
  vout.general(m_vl, "  H_total(trial) = %20.10f\n", H_tot1);

  double plaq1 = m_staple->plaquette(U);
  vout.general(m_vl, "  plaq(trial)    = %20.12f\n", plaq1);


  // Metropolis test
  vout.general(m_vl, "Metropolis test.\n");

  double diff_H = H_tot1 - H_tot0;
  double Arprob = exp(-diff_H);
  double rn     = -1.0;
  if (m_Mtpl_test != 0) {
    rn = m_rand->get();
  }
  vout.general(m_vl, "  H(diff)        = %14.8f\n", diff_H);
  vout.general(m_vl, "  Acceptance     = %14.8f\n", Arprob);
  vout.general(m_vl, "  Random number  = %14.8f\n", rn);

  if (rn <= Arprob) {  // accepted
    vout.general(m_vl, "  Accepted\n");
    conf_org = U;
  } else {  // rejected
    vout.general(m_vl, "  Rejected\n");
  }


  double plaq_final = m_staple->plaquette(conf_org);
  vout.general(m_vl, "  plaq(final)    = %20.12f\n", plaq_final);

  return plaq_final;
}


//====================================================================
double HMC_General::langevin(Field_G& iP, Field_G& U)
{
  int Nc   = CommonParameters::Nc();
  int Lvol = CommonParameters::Lvol();
  int Ndim = CommonParameters::Ndim();
  int NcA  = Nc * Nc - 1;

  int size_iP = NcA * Lvol * Ndim;

  vout.general(m_vl, "Langevin step.\n");

  for (int i = 0; i < m_action.size(); ++i) {
    m_action[i]->notify_linkv();
  }

  for (int i = 0; i < m_director.size(); ++i) {
    m_director[i]->notify_linkv();
  }

  // kinetic term
  double H_iP = m_Langevin_P->set_iP(iP);
  vout.general(m_vl, "  Kinetic term:\n");
  vout.general(m_vl, "    H_kin        = %18.8f\n", H_iP);
  vout.general(m_vl, "    H_kin/dof    = %18.8f\n", H_iP / size_iP);


  double H_actions = 0.0;
  for (int i = 0; i < m_action.size(); ++i) {
    H_actions += m_action[i]->langevin(m_rand);
  }

  double H_tot = H_iP + H_actions;
  return H_tot;
}


//====================================================================
double HMC_General::calc_Hamiltonian(Field_G& iP, Field_G& U)
{
  int Nin  = U.nin();
  int Nvol = U.nvol();
  int Nex  = U.nex();

  int Nc   = CommonParameters::Nc();
  int Nd   = CommonParameters::Nd();
  int Lvol = CommonParameters::Lvol();
  int NcA  = Nc * Nc - 1;

  int size_iP  = NcA * Lvol * Nex;
  // int size_U   = Lvol * 6.0;
  // int size_psf = Lvol * Nc * Nd;


  vout.general(m_vl, "Hamiltonian calculation.\n");

  // kinetic term
  double H_iP = calcH_P(iP);
  vout.general(m_vl, "  Kinetic term:\n");
  vout.general(m_vl, "    H_kin        = %18.8f\n", H_iP);
  vout.general(m_vl, "    H_kin/dof    = %18.8f\n", H_iP / size_iP);

  double H_actions = 0.0;
  for (int i = 0; i < m_action.size(); ++i) {
    H_actions += m_action[i]->calcH();
  }

  double H_tot = H_iP + H_actions;
  return H_tot;
}


//====================================================================
double HMC_General::calcH_P(Field_G& iP)
{
  double hn   = iP.norm();
  double H_iP = 0.5 * hn * hn;

  return H_iP;
}


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