CHROMA
plaq_plus_two_plaq_gaugeact.cc
Go to the documentation of this file.
1 /*! \file
2  * \brief Plaquette plus two-plaquette (plaquette squared) gauge action
3  */
4 
5 #include "chromabase.h"
10 #include "io/aniso_io.h"
11 
12 namespace Chroma
13 {
14 
15  namespace PlaqPlusTwoPlaqGaugeActEnv
16  {
17  namespace
18  {
20  multi1d<LatticeColorMatrix> >* createGaugeAct(XMLReader& xml,
21  const std::string& path)
22  {
23  return new GaugeAct(CreateGaugeStateEnv::reader(xml, path),
24  Params(xml, path));
25  }
26 
27  const std::string name = "PLAQ_PLUS_TWO_PLAQ_GAUGEACT";
28 
29  //! Local registration flag
30  static bool registered = false;
31  }
32 
33  //! Register all the factories
34  bool registerAll()
35  {
36  bool success = true;
37  if (! registered)
38  {
39  success &= TheGaugeActFactory::Instance().registerObject(name, createGaugeAct);
40  registered = true;
41  }
42  return success;
43  }
44 
45 
46  Params::Params(XMLReader& xml_in, const std::string& path)
47  {
48  XMLReader paramtop(xml_in, path);
49 
50  try
51  {
52  read(paramtop, "beta_F", beta_F);
53  read(paramtop, "beta_A", beta_A);
54  }
55  catch( const std::string& e ) {
56  QDPIO::cerr << "Error reading XML: " << e << std::endl;
57  QDP_abort(1);
58  }
59  }
60 
61 
62  //! Compute the action
64  {
65  // Action at the site level
66  multi2d<LatticeReal> plaq_site;
67  this->siteAction(plaq_site, state);
68 
69  // Total action
70  // Fundamental part
71  Double act_F = zero;
72 
73  for(int mu=1; mu < Nd; ++mu)
74  {
75  for(int nu=0; nu < mu; ++nu)
76  {
77  // Sum over plaquettes
78  act_F += sum(plaq_site[mu][nu]);
79  }
80  }
81 
82  // TwoPlaq part
83  Double act_A = zero;
84 
85  for(int mu=1; mu < Nd; ++mu)
86  {
87  for(int nu=0; nu < mu; ++nu)
88  {
89  // Sum over plaquettes
90  act_A += sum(plaq_site[mu][nu] * plaq_site[mu][nu]);
91  }
92  }
93 
94  // Normalize
95  Real act = -param.beta_F * act_F - Real(0.5)*param.beta_A * act_A;
96 
97  return act;
98  }
99 
100 
101  //! Compute the site-level action
102  void GaugeAct::siteAction(multi2d<LatticeReal>& site_act, const Handle< GaugeState<P,Q> >& state) const
103  {
104  START_CODE();
105 
106  // Initialize
107  site_act.resize(Nd,Nd);
108  site_act = zero;
109 
110  // Handle< const GaugeState<P,Q> > u_bc(createState(u));
111  // Apply boundaries
112  const multi1d<LatticeColorMatrix>& u = state->getLinks();
113 
114  // Compute the average plaquettes
115  for(int mu=1; mu < Nd; ++mu)
116  {
117  for(int nu=0; nu < mu; ++nu)
118  {
119  /* tmp_0 = u(x+mu,nu)*u_dag(x+nu,mu) */
120  /* tmp_1 = tmp_0*u_dag(x,nu)=u(x+mu,nu)*u_dag(x+nu,mu)*u_dag(x,nu) */
121  /* wplaq_tmp = tr(u(x,mu)*tmp_1=u(x,mu)*u(x+mu,nu)*u_dag(x+nu,mu)*u_dag(x,nu)) */
122  site_act[mu][nu] += real(trace(u[mu]*shift(u[nu],FORWARD,mu)*adj(shift(u[mu],FORWARD,nu))*adj(u[nu])));
123 
124  // Account for normalization
125  site_act[mu][nu] *= Real(1) / Real(Nc);
126 
127  // Keep a copy
128  site_act[nu][mu] = site_act[mu][nu];
129  }
130  }
131 
132 
133  END_CODE();
134  }
135 
136 
137  //! Compute staple
138  /*! Default version. Derived class should override this if needed. */
139  void GaugeAct::staple(LatticeColorMatrix& result,
140  const Handle< GaugeState<P,Q> >& state,
141  int mu, int cb) const
142  {
143  QDPIO::cerr << __func__ << ": staple not possible\n";
144  QDP_abort(1);
145  }
146 
147 
148  //! Compute dS/dU
149  void GaugeAct::deriv(multi1d<LatticeColorMatrix>& ds_u,
150  const Handle< GaugeState<P,Q> >& state) const
151  {
152  START_CODE();
153 
154  // Derivate at the site level
155  multi1d<LatticeColorMatrix> deriv_fun;
156  this->derivPlaqFun(deriv_fun, state);
157 
158  // Derivate at the site level
159  multi1d<LatticeColorMatrix> deriv_two;
160  this->derivPlaqTwo(deriv_two, state);
161 
162  // Total derivative
163  // Fold in normalization. The (1/2) in the two_plaq is removed in the derivative via the chain rule
164  ds_u.resize(Nd);
165 
166  for(int mu = 0; mu < Nd; mu++)
167  {
168  ds_u[mu] = (-param.beta_F) * deriv_fun[mu];
169  ds_u[mu] += (-param.beta_A) * deriv_two[mu];
170  }
171 
172  END_CODE();
173  }
174 
175 
176  //! Compute dS/dU
177  void GaugeAct::derivPlaqFun(multi1d<LatticeColorMatrix>& ds_u,
178  const Handle< GaugeState<P,Q> >& state) const
179  {
180  START_CODE();
181 
182  LatticeColorMatrix tmp_0;
183  LatticeColorMatrix tmp_1;
184  LatticeColorMatrix tmp_2;
185 
186  const multi1d<LatticeColorMatrix>& u = state->getLinks();
187 
188  ds_u.resize(Nd);
189  ds_u =zero;
190 
191  for(int mu=0; mu < Nd; mu++)
192  {
193  LatticeColorMatrix G;
194  G = zero;
195 
196  for(int nu = 0; nu < Nd; nu++)
197  {
198  if (mu == nu) continue;
199 
200  LatticeColorMatrix tmp_1 = shift(u[nu], FORWARD, mu);
201  LatticeColorMatrix tmp_2 = shift(u[mu], FORWARD, nu);
202 
203  LatticeColorMatrix up_staple = tmp_1*adj(tmp_2)*adj(u[nu]);
204  LatticeColorMatrix down_staple = adj(tmp_1)*adj(u[mu])*u[nu];
205 
206  G += up_staple + shift(down_staple, BACKWARD, nu);
207  }
208 
209  ds_u[mu] = u[mu]*G;
210  }
211 
212  // It is 1/(4Nc) to account for normalisation relevant to fermions
213  // in the taproj, which is a factor of 2 different from the
214  // one used here.
215  for(int mu=0; mu < Nd; mu++)
216  {
217  ds_u[mu] *= Real(1)/(Real(2*Nc));
218  }
219 
220  // Zero the force on any fixed boundaries
221  getGaugeBC().zero(ds_u);
222 
223  END_CODE();
224  }
225 
226 
227 
228  //! Compute dS/dU
229  void GaugeAct::derivPlaqTwo(multi1d<LatticeColorMatrix>& ds_u,
230  const Handle< GaugeState<P,Q> >& state) const
231  {
232  START_CODE();
233 
234  // Action at the site level
235  multi2d<LatticeReal> plaq_site;
236  this->siteAction(plaq_site, state);
237 
238  // Derivative part
239  LatticeColorMatrix tmp_0;
240  LatticeColorMatrix tmp_1;
241  LatticeColorMatrix tmp_2;
242 
243  const multi1d<LatticeColorMatrix>& u = state->getLinks();
244 
245  ds_u.resize(Nd);
246  ds_u =zero;
247 
248  for(int mu=0; mu < Nd; mu++)
249  {
250  LatticeColorMatrix G;
251  G = zero;
252 
253  for(int nu = 0; nu < Nd; nu++)
254  {
255  if (mu == nu) continue;
256 
257  LatticeColorMatrix tmp_1 = shift(u[nu], FORWARD, mu);
258  LatticeColorMatrix tmp_2 = shift(u[mu], FORWARD, nu);
259 
260  LatticeColorMatrix up_staple = tmp_1*adj(tmp_2)*adj(u[nu]);
261  LatticeColorMatrix down_staple = adj(tmp_1)*adj(u[mu])*u[nu];
262 
263  G += up_staple * plaq_site[mu][nu];
264  G += shift(down_staple * plaq_site[mu][nu], BACKWARD, nu);
265  }
266 
267  ds_u[mu] = u[mu]*G;
268  }
269 
270  // It is 1/(4Nc) to account for normalisation relevant to fermions
271  // in the taproj, which is a factor of 2 different from the
272  // one used here.
273  for(int mu=0; mu < Nd; mu++)
274  {
275  ds_u[mu] *= Real(1)/(Real(2*Nc));
276  }
277 
278  // Zero the force on any fixed boundaries
279  getGaugeBC().zero(ds_u);
280 
281  END_CODE();
282  }
283 
284  }
285 
286 }
Anisotropy parameters.
Primary include file for CHROMA library code.
Abstract base class for gauge actions.
Definition: gaugeact.h:25
virtual const GaugeBC< multi1d< LatticeColorMatrix >, multi1d< LatticeColorMatrix > > & getGaugeBC() const
Return the gauge BC object for this action.
Definition: gaugeact.h:46
virtual void zero(P &ds_u) const =0
Zero some gauge-like field in place on the masked links.
Support class for fermion actions and linear operators.
Definition: state.h:74
Class for counted reference semantics.
Definition: handle.h:33
void derivPlaqTwo(multi1d< LatticeColorMatrix > &ds_u, const Handle< GaugeState< P, Q > > &state) const
Compute dS/dU.
void siteAction(multi2d< LatticeReal > &site_act, const Handle< GaugeState< P, Q > > &state) const
Compute the site-level action.
Double S(const Handle< GaugeState< P, Q > > &state) const
Compute the actions.
void deriv(multi1d< LatticeColorMatrix > &result, const Handle< GaugeState< P, Q > > &state) const
Compute dS/dU.
void derivPlaqFun(multi1d< LatticeColorMatrix > &ds_u, const Handle< GaugeState< P, Q > > &state) const
Compute dS/dU.
void staple(LatticeColorMatrix &result, const Handle< GaugeState< P, Q > > &state, int mu, int cb) const
Compute staple.
static T & Instance()
Definition: singleton.h:432
int GaugeAct
int mu
Definition: cool.cc:24
int nu
Definition: cool.cc:25
All gauge create-state method.
Gauge create state factory.
Fermion action factories.
void read(XMLReader &xml, const std::string &path, AsqtadFermActParams &param)
Read parameters.
LatticeColorMatrix tmp_2
Definition: meslate.cc:51
LatticeColorMatrix tmp_1
Definition: meslate.cc:50
Nd
Definition: meslate.cc:74
GaugeAction< multi1d< LatticeColorMatrix >, multi1d< LatticeColorMatrix > > * createGaugeAct(XMLReader &xml, const std::string &path)
static bool registered
Local registration flag.
const std::string name
Name to be used.
Handle< CreateGaugeState< multi1d< LatticeColorMatrix >, multi1d< LatticeColorMatrix > > > reader(XMLReader &xml_in, const std::string &path)
Helper function for the CreateGaugeState readers.
bool registerAll()
Register all the factories.
multi1d< LatticeColorMatrix > G
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
static multi1d< LatticeColorMatrix > u
START_CODE()
int cb
Definition: invbicg.cc:120
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > state
Definition: pbg5p_w.cc:28
Double zero
Definition: invbicg.cc:106
FloatingPoint< double > Double
Definition: gtest.h:7351
::std::string string
Definition: gtest.h:1979
Plaquette plus two-plaquette (plaquette squared) gauge action.
#define FORWARD
Definition: primitives.h:82
#define BACKWARD
Definition: primitives.h:83
Double sum
Definition: qtopcor.cc:37