hmc_Leapfrog.cpp

Go to the documentation of this file.

#include "hmc_Leapfrog.h"

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

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

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


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

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

//====================================================================
void HMC_Leapfrog::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 Estep;
  int    Nmdc, Nprec, Mtpl_test;

  int err = 0;
  err += params.fetch_double("step_size", Estep);
  err += params.fetch_int("number_of_steps", Nmdc);
  err += params.fetch_int("order_of_exp_iP", Nprec);
  err += params.fetch_int("Metropolis_test", Mtpl_test);

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


  set_parameters(Estep, Nmdc, Nprec, Mtpl_test);
}


//====================================================================
void HMC_Leapfrog::set_parameters(double Estep, int Nmdc, int Nprec, int Mtpl_test)
{
  //- print input parameters
  vout.general(m_vl, "HMC (single leapfrog) setup.\n");
  vout.general(m_vl, "  Number of actions: %d\n", m_action.size());
  vout.general(m_vl, "  Estep     = %10.6f\n", Estep);
  vout.general(m_vl, "  Nmdc      = %d\n", Nmdc);
  vout.general(m_vl, "  Nprec     = %d\n", Nprec);
  vout.general(m_vl, "  Mtpl_test = %d\n", Mtpl_test);

  //- range check
  // NB. Estep,Nmdc,Nprec,Mtpl_test == 0 is allowed.

  //- store values
  m_Estep     = Estep;     // step size (SA)
  m_Nmdc      = Nmdc;      // a number of steps (SA)
  m_Nprec     = Nprec;     // order of approximation for e^{iP} (SA)
  m_Mtpl_test = Mtpl_test; // Metropolis test on/off (SA)
}


//====================================================================
double HMC_Leapfrog::update(Field_G& conf_org)
{
  int Nc   = CommonParameters::Nc();   //color (SA)
  int Lvol = CommonParameters::Lvol(); // global volume (SA)

  Field_G U(conf_org);                 // set original gauge conf. to U (SA)

  for (int i = 0; i < m_action.size(); ++i) {
    m_action[i]->set_config(&U); // set gauge conf. for each action (SA)
  }

  int Nin  = U.nin();    // 2x(complex) 3(color) x 3(color) part (SA)
  int Nvol = U.nvol();   // local volume (SA)
  int Nex  = U.nex();    // direction mu (SA)

  Field_G iP(Nvol, Nex); // define momentum iP for gauge fields (SA)


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

  //- Langevin step
  double H_tot0 = langevin(iP, U); // calculate initial hamiltonian (SA)
  // H=p^2/s + S_G + all fermion contributions (SA)
  vout.general(m_vl, "  H_tot(init) = %18.8f\n", H_tot0);

  double plaq0 = m_staple->plaquette(U); // calculate plaquette of old U (SA)
  vout.general(m_vl, "  plaq(init)  = %18.8f\n", plaq0);

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

  //- molecular dynamical integration
  integrate(iP, U); // MD step (SA)

  //- final Hamiltonian
  double H_tot1 = calc_Hamiltonian(iP, U); // calculate final hamiltonian (SA)
  vout.general(m_vl, "  H_tot(fin)  = %18.8f\n", H_tot1);

  double plaq1 = m_staple->plaquette(U); // calculate plaquette of new U (SA)
  vout.general(m_vl, "  plaq(fin)   = %18.8f\n", plaq1);


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

  double diff_H = H_tot1 - H_tot0; // Delta H (SA)
  double Arprob = exp(-diff_H);    // e^{-Delta H} (SA)
  double rn     = -1.0;
  if (m_Mtpl_test != 0) {
    rn = m_rand->get(); //  random number 0 <= rn < 1 (SA)
  }
  vout.general(m_vl, "  H(diff)                = %10.6f\n", diff_H);
  vout.general(m_vl, "  Acceptance probability = %10.6f\n", Arprob);
  vout.general(m_vl, "  Random number          = %10.6f\n", rn);

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

  return m_staple->plaquette(U);
}


//====================================================================
double HMC_Leapfrog::langevin(Field_G& iP, Field_G& U)
{
  int Nc   = CommonParameters::Nc();
  int Lvol = CommonParameters::Lvol(); //global volume (SA)
  int Ndim = CommonParameters::Ndim();
  int NcA  = Nc * Nc - 1;              // Nc^2-1 (SA)

  int size_iP = NcA * Lvol * Ndim;     // global size of momentum iP (SA)

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

  for (int i = 0; i < m_action.size(); ++i) {
    m_action[i]->notify_linkv(); // tell action that gauge field is modified.(SA)
  }

  for (int i = 0; i < m_director.size(); ++i) {
    m_director[i]->notify_linkv(); // tell director that gauge filed is modified.(SA)
  }

  //- kinetic term
  double H_iP = langevin_P(iP); // calculate hamiltonian of momenta iP. (SA)
  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);
    // calculate contribution to hamiltonian from each action. (SA)
  }

  double H_tot = H_iP + H_actions; // total hamiltonian. (SA)
  return H_tot;
}


//====================================================================
double HMC_Leapfrog::langevin_P(Field_G& iP)
{
  int Nin  = iP.nin();
  int Nvol = iP.nvol(); // local volume (SA)
  int Nex  = iP.nex();

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

  // confirm that now gauge group is SU(3)
  assert(NcA == 8);

  Field Hrand(NcA, Nvol, Nex);     // Local random number (SA)

  m_rand->gauss_lex_global(Hrand); // generate gaussian random number Hrand. (SA)

  double sq3r  = 1.0 / sqrt(3.0);
  double sq3r2 = 2.0 * sq3r;

  for (int ex = 0; ex < Nex; ++ex) {
    for (int site = 0; site < Nvol; ++site) {
      //here SU(3) is explicitly assumed. need to generalize.
      // return (jin,site,ex) component(double) of Hrand. (SA)
      double hc1 = Hrand.cmp(0, site, ex);
      double hc2 = Hrand.cmp(1, site, ex);
      double hc3 = Hrand.cmp(2, site, ex);
      double hc4 = Hrand.cmp(3, site, ex);
      double hc5 = Hrand.cmp(4, site, ex);
      double hc6 = Hrand.cmp(5, site, ex);
      double hc7 = Hrand.cmp(6, site, ex);
      double hc8 = Hrand.cmp(7, site, ex);

      /*
      iP = i P^a T^a: T^a: Gellmann matrix
       P=
       | hc3+hc8/sqrt(3), hc1-i*hc2,       hc4-i*hc5     |
       | hc1+i*hc2,      -hc3+hc8/sqrt(3), hc6-i*hc7     |
       | hc4+i*hc5,       hc6+i*hc7,      -2*hc8/sqrt(3) |
      */
      iP.set(0, site, ex, 0.0);
      iP.set(1, site, ex, hc3 + hc8 * sq3r);
      iP.set(2, site, ex, hc2);
      iP.set(3, site, ex, hc1);
      iP.set(4, site, ex, hc5);
      iP.set(5, site, ex, hc4);
      iP.set(6, site, ex, -hc2);
      iP.set(7, site, ex, hc1);
      iP.set(8, site, ex, 0.0);
      iP.set(9, site, ex, -hc3 + hc8 * sq3r);
      iP.set(10, site, ex, hc7);
      iP.set(11, site, ex, hc6);
      iP.set(12, site, ex, -hc5);
      iP.set(13, site, ex, hc4);
      iP.set(14, site, ex, -hc7);
      iP.set(15, site, ex, hc6);
      iP.set(16, site, ex, 0.0);
      iP.set(17, site, ex, -hc8 * sq3r2);
    }
  }

  double iP2 = iP.norm(); //calculate global norm sqrt(|iP2|^2) (SA)
  iP2 = 0.5 * iP2 * iP2;  // |iP|^2/2 (SA)

  return iP2;
}


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

  int Nc   = CommonParameters::Nc();
  int Nd   = CommonParameters::Nd();
  int Lvol = CommonParameters::Lvol(); //Global volume (SA)
  int NcA  = Nc * Nc - 1;

  int size_iP  = NcA * Lvol * Nex; // Global size of iP (SA)
  int size_U   = Lvol * 6.0;       // Global size of U (SA)
  int size_psf = Lvol * Nc * Nd;   // Global size of pseudo-fermion (SA)


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

  // kinetic term
  double H_iP = calcH_P(iP); // calculate hamiltonian for iP (SA)
  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) {
    // calculate contribution to hamiltonian from each action (SA)
    H_actions += m_action[i]->calcH();
  }

  double H_tot = H_iP + H_actions;
  return H_tot;
}


//====================================================================
double HMC_Leapfrog::calcH_P(Field_G& iP)
{
  //calculate hamiltonian for iP: |iP|^2/2 (SA)
  double hn   = iP.norm();
  double H_iP = 0.5 * hn * hn;

  return H_iP;
}


//====================================================================
void HMC_Leapfrog::integrate(Field_G& iP, Field_G& U)
{
  vout.general(m_vl, "Integration.\n");

  double estep  = m_Estep;
  double estep2 = 0.5 * m_Estep;

  // Initial half step of update of h
  if (m_Nmdc > 0) {
    int imdc = 0;
    vout.general(m_vl, "  imdc = %d\n", imdc);
    update_P(estep2, iP, U); // update momentum iP at first half step (SA)
  }

  // Molecular dynamics step
  for (int imdc = 1; imdc < m_Nmdc + 1; imdc++) {
    vout.general(m_vl, "  imdc = %d\n", imdc);

    update_U(estep, iP, U); // update gauge field U at full step (SA)

    if (imdc < m_Nmdc) {
      update_P(estep, iP, U);  // update momentum iP at full step (SA)
    } else {
      update_P(estep2, iP, U); // update momentum iP at last half step (SA)
    }
  }  // here imdc loop ends
}


//====================================================================
void HMC_Leapfrog::update_P(double estep, Field_G& iP, Field_G& U)
{
  int Nin  = U.nin();
  int Nvol = U.nvol();
  int Nex  = U.nex();
  int Nc   = CommonParameters::Nc();

  Field force(Nin, Nvol, Nex);

  force = 0.0;

  Field  force1(Nin, Nvol, Nex);
  double H_actions = 0.0;
  for (int i = 0; i < m_action.size(); ++i) {
    force1 = m_action[i]->force(); // calculate force for each action (SA)
    force += estep * force1;
  }

  iP += (Field_G)force; // iP = iP + step-size * force (SA)
}


//====================================================================
void HMC_Leapfrog::update_U(double estep, Field_G& iP, Field_G& U)
{
  int Nvol = U.nvol();
  int Nex  = U.nex();
  int Nc   = CommonParameters::Nc();

  Mat_SU_N u0(Nc), u1(Nc), u2(Nc);
  Mat_SU_N h1(Nc);

  for (int ex = 0; ex < Nex; ++ex) {
    for (int site = 0; site < Nvol; ++site) {
      u0 = U.mat(site, ex);
      u1 = U.mat(site, ex);
      h1 = iP.mat(site, ex);

      for (int iprec = 0; iprec < m_Nprec; ++iprec) {
        double exf = estep / (m_Nprec - iprec); // step_size/(N-k) (SA)
        u2  = h1 * u1;                          // iP* u1 (SA)
        u2 *= exf;                              // step_size*iP*u1/(N-K)  (SA)
        u1  = u2;
        u1 += u0;                               // U + step_size*iP*u1/(N-K) (SA)
      }
      // u1 =sum_{k=0}^{N-1} (step_size * iP)^k/k!, N=m_Nprec (SA)
      u1.reunit();             // reunitarize u1 (SA)
      U.set_mat(site, ex, u1); // U = u1 (SA)
    }
  }

  for (int i = 0; i < m_action.size(); ++i) {
    m_action[i]->notify_linkv();   // notify all actions about update of U (SA)
  }

  for (int i = 0; i < m_director.size(); ++i) {
    m_director[i]->notify_linkv(); // notify all directors about update of U (SA)
  }
}


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