topologicalCharge.cpp

Go to the documentation of this file.

#include "topologicalCharge.h"

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


//- parameter entries
namespace {
  void append_entry(Parameters& param)
  {
    param.Register_double("c_plaq", 0.0);
    param.Register_double("c_rect", 0.0);

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


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

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

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

//====================================================================
void TopologicalCharge::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 c_plaq, c_rect;

  int err = 0;
  err += params.fetch_double("c_plaq", c_plaq);
  err += params.fetch_double("c_rect", c_rect);

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


  set_parameters(c_plaq, c_rect);
}


//====================================================================
void TopologicalCharge::set_parameters(double c_plaq, double c_rect)
{
  //- print input parameters
  vout.general(m_vl, "Topological Charge measurement:\n");
  vout.general(m_vl, "  c_plaq = %12.6f\n", c_plaq);
  vout.general(m_vl, "  c_rect = %12.6f\n", c_rect);

  //- range check
  // NB. beta,c_plaq,c_rect == 0 is allowed.

  //- store values
  m_c_plaq = c_plaq;
  m_c_rect = c_rect;

  //- post-process
}


//====================================================================
double TopologicalCharge::measure(Field_G& U)
{
  int                 Ndim = CommonParameters::Ndim();
  int                 Nvol = CommonParameters::Nvol();
  static const double PI   = 4.0 * atan(1.0);
  static const double PI2  = PI * PI;

  // NB. #(mu,nu)=6 i.e. (1,2),(1,3),(1,4),(2,3),(2,4),(3,4)
  std::vector<Field_G> Fmunu_1x1(6), Fmunu_1x2(6);


  vout.paranoiac(m_vl, "%s: measurement start.\n", class_name.c_str());

  int i_munu = 0;
  for (int mu = 0; mu < Ndim; ++mu) {
    for (int nu = mu + 1; nu < Ndim; ++nu) {
      construct_Fmunu_1x1(Fmunu_1x1[i_munu], mu, nu, U);
      construct_Fmunu_1x2(Fmunu_1x2[i_munu], mu, nu, U);

      ++i_munu;
    }
  }


  double Q_1x1 = 0.0;

  // (mu,nu, rho,sigma) = (1,2, 3,4)
  Q_1x1 += contract_epsilon_tensor(Fmunu_1x1[0], Fmunu_1x1[5]);

  // (mu,nu, rho,sigma) = (1,3, 2,4)
  Q_1x1 -= contract_epsilon_tensor(Fmunu_1x1[1], Fmunu_1x1[4]);

  // (mu,nu, rho,sigma) = (1,4, 2,3)
  Q_1x1 += contract_epsilon_tensor(Fmunu_1x1[2], Fmunu_1x1[3]);

  // #degeneracy of (mu,nu, rho,sigma) is 8
  Q_1x1 *= 8.0;


  double Q_1x2 = 0.0;
  Q_1x2 += contract_epsilon_tensor(Fmunu_1x2[0], Fmunu_1x2[5]);
  Q_1x2 -= contract_epsilon_tensor(Fmunu_1x2[1], Fmunu_1x2[4]);
  Q_1x2 += contract_epsilon_tensor(Fmunu_1x2[2], Fmunu_1x2[3]);
  Q_1x2 *= 8.0;


  double Q_topo = (m_c_plaq * Q_1x1 + m_c_rect * Q_1x2 / 2.0) / (32.0 * PI2);



  vout.general(m_vl, "  Q_1x1,Q_1x2 = %20.16e %20.16e\n",
               Q_1x1 / (32.0 * PI2),
               Q_1x2 / (2.0 * (32.0 * PI2)));

  vout.paranoiac(m_vl, "%s: measurement finished.\n", class_name.c_str());


  return Q_topo;
}


//====================================================================
void TopologicalCharge::construct_Fmunu_1x1(Field_G& Fmunu_1x1,
                                            const int mu, const int nu, Field_G& U)
{
  int Nvol = CommonParameters::Nvol();

  Staples staple;

  Field_G Cup(Nvol, 1), Cdn(Nvol, 1);
  Field_G Umu(Nvol, 1);
  Field_G v(Nvol, 1), v2(Nvol, 1);


  //- building blocks
  // (1)  mu (2)
  //    +->-+
  // nu |   |
  //   i+   +
  staple.upper(Cup, U, mu, nu);

  // (1)  mu (2)
  //   i+   +
  // nu |   |
  //    +->-+
  staple.lower(Cdn, U, mu, nu);
  Umu.setpart_ex(0, U, mu);


  //- right clover leaf
  // (2)  +-<-+
  //  nu  |   |
  //     i+->-+
  //      mu (1)
  mult_Field_Gnd(Fmunu_1x1, 0, Umu, 0, Cup, 0);

  // (1) i+-<-+
  //  nu  |   |
  //      +->-+
  //      mu (2)
  multadd_Field_Gnd(Fmunu_1x1, 0, Cdn, 0, Umu, 0, 1.0);


  //- left clover leaf
  //   mu (2)
  //   +-<-+ (1)
  //   |   | nu    mu (1)
  //   +->-+i   +  +-<-+i (2)
  //               |   | nu
  //               +->-+
  mult_Field_Gdn(v, 0, Cup, 0, Umu, 0);
  multadd_Field_Gdn(v, 0, Umu, 0, Cdn, 0, 1.0);

  // NB. shift.forward(mu)
  //   mu (2)        mu (2)
  //   +-<-+ (1)     +-<-+ (1)
  //   |   | nu  ->  |   | nu
  //  i+->-+         +->-+i
  m_shift.forward(v2, v, mu);

  axpy(Fmunu_1x1, 1.0, v2);

  ah_Field_G(Fmunu_1x1, 0);

  scal(Fmunu_1x1, -0.25);
  Fmunu_1x1.xI();
}


//====================================================================
void TopologicalCharge::construct_Fmunu_1x2(Field_G& Fmunu_1x2,
                                            const int mu, const int nu, Field_G& U)
{
  int Nvol = CommonParameters::Nvol();

  Staples staple;

  Field_G Cup1(Nvol, 1), Cup2(Nvol, 1);
  Field_G Cdn1(Nvol, 1), Cdn2(Nvol, 1);
  Field_G Umu(Nvol, 1), Unu(Nvol, 1);
  Field_G Umu_nu(Nvol, 1), Unu_mu(Nvol, 1);
  Field_G rect(Nvol, 1);
  Field_G v(Nvol, 1), w(Nvol, 1), c(Nvol, 1);

  //-- building blocks
  //      mu (2)
  // (1)  +->-+
  //  nu  |   |
  //     i+   +
  staple.upper(Cup1, U, mu, nu);

  //      +-<-+ (2)
  //          | nu
  //     i+->-+
  //      mu (1)
  staple.upper(Cup2, U, nu, mu);

  // (1) i+   +
  //  nu  |   |
  //      +->-+
  //      mu (2)
  staple.lower(Cdn1, U, mu, nu);

  // (2)  +->-+
  //  nu  |
  //      +-<-+i
  //      mu (1)
  staple.lower(Cdn2, U, nu, mu);

  Umu.setpart_ex(0, U, mu);
  Unu.setpart_ex(0, U, nu);

  //  +->-+
  //
  // i+
  m_shift.backward(Umu_nu, Umu, nu);

  //      +
  //      |
  // i+   +
  m_shift.backward(Unu_mu, Unu, mu);

  //============================================================
  //-- 2x1 part

  //- upper-right(2x1)
  //    +---<---+      +      +-<-+         <---+
  // nu |       |  =   |  +          +          |  +
  //   i+--->---+     i+     i+         i+  >---+     i+->-+
  m_shift.backward(c, Cup2, mu);

  mult_Field_Gdd(v, 0, Umu_nu, 0, Unu, 0);
  mult_Field_Gnn(w, 0, c, 0, v, 0);

  mult_Field_Gnn(rect, 0, Umu, 0, w, 0);

  Fmunu_1x2 = rect;

  //- lower-right(2x1)
  // (1) +---<---+     (1) i+---<---+
  //  nu |       |  ->  nu  |       |
  //    i+--->---+          +--->---+
  //      mu (2)             mu (2)
  mult_Field_Gnd(v, 0, c, 0, Umu_nu, 0);
  mult_Field_Gnn(w, 0, Umu, 0, v, 0);

  mult_Field_Gdn(rect, 0, Unu, 0, w, 0);

  m_shift.forward(v, rect, nu);

  axpy(Fmunu_1x2, 1.0, v);

  //- upper-left(2x1)
  //      mu (2)
  //  +---<---+ (1)              +-<-+      +-<-+           +
  //  |       | nu  =         +  |       +         +        |
  //  +---i---+       i+->-+     +->-+i     +i         +i   +
  //
  // NB. shift.forward(mu)
  //      mu (2)             mu (2)
  //  +---<---+ (1)      +---<---+ (1)
  //  |       | nu  ->   |       | nu
  //  +---i---+          +--->---+i
  mult_Field_Gdn(v, 0, Cdn2, 0, Umu, 0);
  mult_Field_Gdn(w, 0, Umu_nu, 0, v, 0);

  mult_Field_Gnn(rect, 0, Unu_mu, 0, w, 0);

  m_shift.forward(v, rect, mu);

  axpy(Fmunu_1x2, 1.0, v);

  //- lower-left(2x1)
  //      mu (1)
  //  +---<---+ (2)        +                +-<-+      +-<-+
  //  |       | nu  =      |  +          +  |       +
  //  +---i---+       +i   +     i+->-+     +->-+i     +i
  //
  // NB. shift.forward(mu+nu)
  //      mu (1)             mu (1)            mu (1)
  //  +---<---+ (2)      +---<---+ (2)     +---<---+i (2)
  //  |       | nu  ->   |       | nu  ->  |       |  nu
  //  +---i---+          +--->---+i        +--->---+
  mult_Field_Gnn(v, 0, Umu, 0, Unu_mu, 0);
  mult_Field_Gdn(w, 0, Cdn2, 0, v, 0);

  mult_Field_Gdn(rect, 0, Umu_nu, 0, w, 0);

  m_shift.forward(w, rect, mu);
  m_shift.forward(v, w, nu);

  axpy(Fmunu_1x2, 1.0, v);
  //============================================================


  //============================================================
  //-- 1x2 part

  //- upper-right(1x2)
  //     +-<-+             +-<-+
  //     |   |             |   |
  // (2) +   +  =  +   +   +   +   +      +   +
  //  nu |   |     |                      |
  //    i+->-+    i+      i+         i+   +      i+->-+
  //     mu (1)
  m_shift.backward(c, Cup1, nu);

  mult_Field_Gdd(v, 0, c, 0, Unu, 0);
  mult_Field_Gnn(w, 0, Unu_mu, 0, v, 0);

  mult_Field_Gnn(rect, 0, Umu, 0, w, 0);

  axpy(Fmunu_1x2, 1.0, rect);

  //- lower-right(1x2)
  //      mu (2)
  // (1)  +-<-+      +-<-+           +                     +
  //  nu  |   |                      |                     |
  //     i+   +  =  i+      +   i+   +   +   i+   +   +   i+
  //      |   |                               |   |
  //      +->-+                               +->-+
  //
  // NB. shift.forward(nu)
  //      mu (2)        mu (2)
  // (1)  +-<-+    (1) i+-<-+
  //  nu  |   |     nu  |   |
  //     i+   +  ->     +   +
  //      |   |         |   |
  //      +->-+         +->-+
  mult_Field_Gnd(v, 0, Unu_mu, 0, Umu_nu, 0);
  mult_Field_Gnn(w, 0, Cdn1, 0, v, 0);

  mult_Field_Gdn(rect, 0, Unu, 0, w, 0);

  m_shift.forward(v, rect, nu);

  axpy(Fmunu_1x2, 1.0, v);

  //- upper-left(1x2)
  //  +-<-+                            +-<-+
  //  |   |                            |   |
  //  +   + (1)  =         +   +   +   +   +   +       +
  //  |   | nu                 |                       |
  // i+->-+        i+->-+     i+      i+          i+   +
  //  mu (2)
  //
  // NB. shift.forward(mu)
  //  +-<-+         +-<-+
  //  |   |         |   |
  //  +   + (1) ->  +   + (1)
  //  |   | nu      |   | nu
  // i+->-+         +->-+i
  //  mu (2)        mu (2)
  mult_Field_Gdn(v, 0, Unu, 0, Umu, 0);
  mult_Field_Gdn(w, 0, c, 0, v, 0);

  mult_Field_Gnn(rect, 0, Unu_mu, 0, w, 0);

  m_shift.forward(v, rect, mu);

  axpy(Fmunu_1x2, 1.0, v);

  //- lower-left(1x2)
  //      mu (1)
  // (2)  +-<-+          +                     +        +-<-+
  //  nu  |   |          |                     |
  //     i+   +  =  i+   +   +   i+   +   +   i+   +   i+
  //      |   |                   |   |
  //      +->-+                   +->-+
  //
  // NB. shift.forward(mu+nu)
  //      mu (1)        mu (1)
  // (2)  +-<-+    (2)  +-<-+i
  //  nu  |   |     nu  |   |
  //     i+   +  ->     +   +
  //      |   |         |   |
  //      +->-+         +->-+
  mult_Field_Gnn(v, 0, Cdn1, 0, Unu_mu, 0);
  mult_Field_Gdn(w, 0, Unu, 0, v, 0);

  mult_Field_Gdn(rect, 0, Umu_nu, 0, w, 0);

  m_shift.forward(w, rect, mu);
  m_shift.forward(v, w, nu);

  axpy(Fmunu_1x2, 1.0, v);
  //============================================================

  ah_Field_G(Fmunu_1x2, 0);

  scal(Fmunu_1x2, -0.25);
  Fmunu_1x2.xI();
}


//====================================================================
double TopologicalCharge::contract_epsilon_tensor(Field_G& Fmunu_1, Field_G& Fmunu_2)
{
  int Nvol = CommonParameters::Nvol();

  double Q_topo = 0.0;

  for (int site = 0; site < Nvol; ++site) {
    Q_topo += ReTr(Fmunu_1.mat(site) * Fmunu_2.mat(site));
  }
  double result = Communicator::reduce_sum(Q_topo);

  return result;
}


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