hmc_General.cpp

Go to the documentation of this file.

 
#include "hmc_General.h"
 
const std::string HMC_General::class_name = "HMC_General";
 
//====================================================================
HMC_General::HMC_General(const ActionList& action_list,
                         std::vector<Director *> director,
                         Integrator *integrator,
                         RandomNumbers *rand)
  : m_vl(CommonParameters::Vlevel())
{
  ActionSet action_set = action_list.get_actions();
 
  m_action.resize(action_set.size());
  for (int i = 0; i < action_set.size(); ++i) {
    m_action[i] = action_set[i];
  }
  m_director.resize(director.size());
  for (int i = 0; i < director.size(); ++i) {
    m_director[i] = director[i];
  }
  m_integrator        = integrator;
  m_rand              = rand;
  m_staple            = new Staple_lex;
  m_Metropolis_test   = false;
  m_Langevin_P        = new Langevin_Momentum(m_rand);
  m_trajectory_length = 0.0;
}
 
 
//====================================================================
//  with params
HMC_General::HMC_General(const ActionList& action_list,
                         std::vector<Director *> director,
                         Integrator *integrator,
                         RandomNumbers *rand,
                         const Parameters& params)
  : m_vl(CommonParameters::Vlevel())
{
  ActionSet action_set = action_list.get_actions();
 
  m_action.resize(action_set.size());
  for (int i = 0; i < action_set.size(); ++i) {
    m_action[i] = action_set[i];
  }
  m_director.resize(director.size());
  for (int i = 0; i < director.size(); ++i) {
    m_director[i] = director[i];
  }
  m_integrator        = integrator;
  m_rand              = rand;
  m_staple            = new Staple_lex;
  m_Metropolis_test   = false;
  m_Langevin_P        = new Langevin_Momentum(m_rand);
  m_trajectory_length = 0.0;
 
  set_parameters(params);
}
 
 
//====================================================================
HMC_General::HMC_General(const ActionList& action_list,
                         Integrator *integrator,
                         RandomNumbers *rand)
  : m_vl(CommonParameters::Vlevel())
{
  ActionSet action_set = action_list.get_actions();
 
  m_action.resize(action_set.size());
  for (int i = 0; i < action_set.size(); ++i) {
    m_action[i] = action_set[i];
  }
 
  m_director.resize(0);
 
  m_integrator        = integrator;
  m_rand              = rand;
  m_staple            = new Staple_lex;
  m_Metropolis_test   = false;
  m_Langevin_P        = new Langevin_Momentum(m_rand);
  m_trajectory_length = 0.0;
}
 
 
//====================================================================
//  with params
HMC_General::HMC_General(const ActionList& action_list,
                         Integrator *integrator,
                         RandomNumbers *rand,
                         const Parameters& params)
  : m_vl(CommonParameters::Vlevel())
{
  ActionSet action_set = action_list.get_actions();
 
  m_action.resize(action_set.size());
  for (int i = 0; i < action_set.size(); ++i) {
    m_action[i] = action_set[i];
  }
 
  m_director.resize(0);
 
  m_integrator        = integrator;
  m_rand              = rand;
  m_staple            = new Staple_lex;
  m_Metropolis_test   = false;
  m_Langevin_P        = new Langevin_Momentum(m_rand);
  m_trajectory_length = 0.0;
 
  set_parameters(params);
}
 
 
//====================================================================
HMC_General::~HMC_General()
{
  delete m_Langevin_P;
  delete m_staple;
}
 
 
//====================================================================
void HMC_General::set_parameters(const Parameters& params)
{
  std::string vlevel;
  if (!params.fetch_string("verbose_level", vlevel)) {
    m_vl = vout.set_verbose_level(vlevel);
  }
 
  //- fetch and check input parameters
  double traj_length;
  bool   Metropolis_test;
 
  int err = 0;
  err += params.fetch_double("trajectory_length", traj_length);
  err += params.fetch_bool("Metropolis_test", Metropolis_test);
 
  if (err) {
    vout.crucial(m_vl, "Error at %s: input parameter not found.\n",
                 class_name.c_str());
    exit(EXIT_FAILURE);
  }
 
  set_parameters(traj_length, Metropolis_test);
}
 
 
//====================================================================
void HMC_General::get_parameters(Parameters& params) const
{
  params.set_double("trajectory_length", m_trajectory_length);
  params.set_bool("Metropolis_test", m_Metropolis_test);
 
  params.set_string("verbose_level", vout.get_verbose_level(m_vl));
}
 
 
//====================================================================
void HMC_General::set_parameters(const double trajectory_length, const int Metropolis_test)
{
  vout.crucial(m_vl, "%s: warning: integer variable for Metroplis_test is obsolete. use boolean parameter.\n", class_name.c_str());
 
  return set_parameters(trajectory_length, (Metropolis_test == 0) ? false : true);
}
 
 
//====================================================================
void HMC_General::set_parameters(const double trajectory_length, const bool Metropolis_test)
{
  //- print input parameters
  vout.general(m_vl, "%s:\n", class_name.c_str());
  vout.general(m_vl, "  Number of actions: %4d\n", m_action.size());
  vout.general(m_vl, "  traj_length     = %8.6f\n", trajectory_length);
  vout.general(m_vl, "  Metropolis_test = %s\n", Metropolis_test ? "true" : "false");
 
  //- range check
  // NB. Metropolis_test == 0 is allowed.
 
  //- store values
  m_trajectory_length = trajectory_length;
  m_Metropolis_test   = Metropolis_test;
}
 
 
//====================================================================
double HMC_General::update(Field_G& Uorg)
{
  const int Nc = CommonParameters::Nc();
 
  Field_G U(Uorg);
 
  for (int i = 0; i < m_action.size(); ++i) {
 
    m_action[i]->set_config(&U);
  }
 
  const int Nin  = U.nin();
  const int Nvol = U.nvol();
  const int Nex  = U.nex();
 
  Field_G iP(Nvol, Nex);
 
  vout.general(m_vl, "\n");
  vout.general(m_vl, "HMC (general) start.\n");
  vout.general(m_vl, "\n");
 
  // Langevin step
  const double H_total0 = langevin(iP, U);
  vout.general(m_vl, "H_total(init)  = %20.10f\n", H_total0);
 
  const double plaq0 = m_staple->plaquette(U);
  vout.general(m_vl, "plaq(init)     = %20.12f\n", plaq0);
  vout.general(m_vl, "\n");
 
  // initial Hamiltonian
  //    H_total0 = calc_Hamiltonian(iP,U);
 
  // molecular dynamical integration
  m_integrator->evolve(m_trajectory_length, iP, U);
  vout.general(m_vl, "\n");
 
  // trial Hamiltonian
  const double H_total1 = calc_Hamiltonian(iP, U);
  vout.general(m_vl, "H_total(trial) = %20.10f\n", H_total1);
 
  const 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");
 
  const double diff_H           = H_total1 - H_total0;
  const double exp_minus_diff_H = exp(-diff_H);
 
  double rand = -1.0;
  if (m_Metropolis_test) {
    rand = m_rand->get();
  }
  vout.general(m_vl, "  H(diff)       = %14.8f\n", diff_H);
  vout.general(m_vl, "  exp(-diff_H)  = %14.8f\n", exp_minus_diff_H);
  vout.general(m_vl, "  Random number = %14.8f\n", rand);
 
  if (rand <= exp_minus_diff_H) {  // accepted
    vout.general(m_vl, "  Metropolis Accepted\n");
    Uorg = U;
  } else {  // rejected
    vout.general(m_vl, "  Metropolis Rejected\n");
  }
 
  const double plaq_final = m_staple->plaquette(Uorg);
  vout.general(m_vl, "plaq(final)    = %20.12f\n", plaq_final);
 
  return plaq_final;
}
 
 
//====================================================================
double HMC_General::langevin(Field_G& iP, const Field_G& U)
{
  const int Nc   = CommonParameters::Nc();
  const int Ndim = CommonParameters::Ndim();
  const int NcA  = Nc * Nc - 1;
 
  const int Nvol = CommonParameters::Nvol();
  const int NPE  = CommonParameters::NPE();
 
  vout.general(m_vl, "Langevin step:\n");
 
  // discard caches
  m_integrator->invalidate_cache();
 
  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);
 
  double Fgauge = 1.0/(double(NcA * Nvol * Ndim) * double(NPE));
  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 * Fgauge);
 
  vout.increase_indent();
 
  double H_actions = 0.0;
  for (int i = 0; i < m_action.size(); ++i) {
    H_actions += m_action[i]->langevin(m_rand);
  }
 
  vout.decrease_indent();
 
  double H_total = H_iP + H_actions;
  return H_total;
}
 
 
//====================================================================
double HMC_General::calc_Hamiltonian(const Field_G& iP, const Field_G& U)
{
  const int Nin  = U.nin();
  const int Nvol = U.nvol();
  const int Nex  = U.nex();
 
  const int Nc  = CommonParameters::Nc();
  const int Nd  = CommonParameters::Nd();
  const int NcA = Nc * Nc - 1;
 
  const int NPE = CommonParameters::NPE();
 
  vout.general(m_vl, "Hamiltonian calculation:\n");
 
  // kinetic term
  double H_iP   = calcH_P(iP);
  double Fgauge = 1.0/(double(NcA * Nvol * Nex) * double(NPE));
  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 * Fgauge);
 
  vout.increase_indent();
 
  double H_actions = 0.0;
  for (int i = 0; i < m_action.size(); ++i) {
    H_actions += m_action[i]->calcH();
  }
 
  vout.decrease_indent();
 
  double H_total = H_iP + H_actions;
  return H_total;
}
 
 
//====================================================================
double HMC_General::calcH_P(const Field_G& iP)
{
  const double hn   = iP.norm();
  const double H_iP = 0.5 * hn * hn;
 
  return H_iP;
}
 
 
//============================================================END=====