integrator.cpp

Go to the documentation of this file.

#include "hmc_General.h"

//====================================================================
void HMC_General::set_parameters(double Estep, int Nmdc,
                                 int Mtpl_test,
                                 int Naction0, int Naction1, int Naction2,
                                 int Nmd1, int Nmd2)
{
  vout.general(m_vl, "HMC (3-level leapfrog) setup:\n");

  vout.general(m_vl, "  Number of actions: %4d\n", m_action.size());

  m_Estep     = Estep;
  m_Nmdc      = Nmdc;
  m_Mtpl_test = Mtpl_test;

  m_Naction0 = Naction0;
  m_Naction1 = Naction1;
  m_Naction2 = Naction2;
  m_Nmd1     = Nmd1;
  m_Nmd2     = Nmd2;

  assert(m_action.size() == m_Naction0 + m_Naction1 + m_Naction2);

  m_Nprec = 8;

  vout.general(m_vl, "  Estep    = %10.6f\n", m_Estep);
  vout.general(m_vl, "  Nmdc     = %4d\n", m_Nmdc);
  vout.general(m_vl, "  Naction0 = %4d\n", m_Naction0);
  vout.general(m_vl, "  Naction1 = %4d\n", m_Naction1);
  vout.general(m_vl, "  Naction2 = %4d\n", m_Naction2);
  vout.general(m_vl, "  Nmd1     = %4d\n", m_Nmd1);
  vout.general(m_vl, "  Nmd2     = %4d\n", m_Nmd2);
  vout.general(m_vl, "  Nprec    = %4d\n", m_Nprec);
}


//====================================================================
void HMC_General::update(Field_G& conf_org)
{
  vout.general(m_vl, "HMC (single leapfrog) start.\n");

  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();

  int Nc   = CommonParameters::Nc();
  int Lvol = CommonParameters::Lvol();

  Field_G iP(Nvol, Nex);

  // Langevin step
  double H_tot0 = langevin(iP, U);
  vout.general(m_vl, "  H_tot(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
  integrate(iP, U);

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

  double plaq1 = m_staple->plaquette(&U);
  vout.general(m_vl, "  plaq(fin)      = %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 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 = langevin_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]->langevin(m_rand);
  }

  double H_tot = H_iP + H_actions;
  return H_tot;
}


//====================================================================
double HMC_General::langevin_P(Field_G& iP)
{
  int Nin  = iP.nin();
  int Nvol = iP.nvol();
  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);

  m_rand->gauss_lex_global(Hrand);

  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.
      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.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();
  iP2 = 0.5 * iP2 * iP2;

  return iP2;
}


//====================================================================
double HMC_General::calc_Hamiltonian(Field_G& iP, Field_G& U)
{
  vout.general(m_vl, "Hamiltonian calculation.\n");

  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;

  // 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;
}


//====================================================================
void HMC_General::integrate(Field_G& iP, Field_G& U)
{
  vout.general(m_vl, "Integration level-0 start.\n");

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

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

  Field force(Nin, Nvol, Nex), force1(Nin, Nvol, Nex);

  // Initial half step of update of iP
  if (m_Nmdc > 0) {
    int imdc = 0;
    vout.general(m_vl, "imdc = %d\n", imdc);
    force = 0.0;
    for (int i = 0; i < m_Naction0; ++i) {
      force1 = m_action[i]->force();
      force += esteph * force1;
    }
    iP += (Field_G)force;
  }

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

    integrate1(iP, U);

    estep2 = estep;
    if (imdc == m_Nmdc) estep2 = esteph;

    force = 0.0;
    for (int i = 0; i < m_Naction0; ++i) {
      force1 = m_action[i]->force();
      force += estep2 * force1;
    }
    iP += (Field_G)force;
  }  // here imdc loop ends

  vout.general(m_vl, "Integration level-0 finished.\n");
}


//====================================================================
void HMC_General::integrate1(Field_G& iP, Field_G& U)
{
  vout.general(m_vl, "Integration level-1 start.\n");

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

  double estep  = m_Estep / ((double)m_Nmd1);
  double esteph = 0.5 * estep;
  double estep2;

  Field force(Nin, Nvol, Nex), force1(Nin, Nvol, Nex);

  // Initial half step of update of h
  if (m_Nmd1 > 0) {
    int imd1 = 0;
    vout.general(m_vl, "imd1 = %d\n", imd1);
    force = 0.0;
    for (int i = 0; i < m_Naction1; ++i) {
      int i2 = i + m_Naction0;
      force1 = m_action[i2]->force();
      force += esteph * force1;
    }
    iP += (Field_G)force;
  }

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

    integrate2(iP, U);

    estep2 = estep;
    if (imd1 == m_Nmd1) estep2 = esteph;

    force = 0.0;
    for (int i = 0; i < m_Naction1; ++i) {
      int i2 = i + m_Naction0;
      force1 = m_action[i2]->force();
      force += estep2 * force1;
    }
    iP += (Field_G)force;
  }  // here imdc loop ends

  vout.general(m_vl, "Integration level-1 finished.\n");
}


//====================================================================
void HMC_General::integrate2(Field_G& iP, Field_G& U)
{
  vout.general(m_vl, "Integration level-2 start.\n");

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

  double estep  = m_Estep / ((double)(m_Nmd1 * m_Nmd2));
  double esteph = 0.5 * estep;
  double estep2;

  Field force(Nin, Nvol, Nex), force1(Nin, Nvol, Nex);

  // Initial half step of update of h
  if (m_Nmd2 > 0) {
    int imd2 = 0;
    vout.general(m_vl, "imd2 = %d\n", imd2);
    force = 0.0;
    for (int i = 0; i < m_Naction2; ++i) {
      int i2 = i + m_Naction0 + m_Naction1;
      force1 = m_action[i2]->force();
      force += esteph * force1;
    }
    iP += (Field_G)force;
  }

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

    update_U(estep, iP, U);

    estep2 = estep;
    if (imd2 == m_Nmd2) estep2 = esteph;

    force = 0.0;
    for (int i = 0; i < m_Naction2; ++i) {
      int i2 = i + m_Naction0 + m_Naction1;
      force1 = m_action[i2]->force();
      force += estep2 * force1;
    }
    iP += (Field_G)force;
  }  // here imdc loop ends

  vout.general(m_vl, "Integration level-2 finished.\n");
}


//====================================================================
void HMC_General::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);
        u2  = h1 * u1;
        u2 *= exf;
        u1  = u2;
        u1 += u0;
      }

      u1.reunit();
      U.set_mat(site, ex, u1);
    }
  }

  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();
  }
}


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