Bridge++  Ver. 1.1.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
action_G_Rectangle.cpp
Go to the documentation of this file.
1 
14 #include "action_G_Rectangle.h"
15 
16 #ifdef USE_PARAMETERS_FACTORY
17 #include "parameters_factory.h"
18 #endif
19 
20 //- parameter entries
21 namespace {
22  void append_entry(Parameters& param)
23  {
24  param.Register_double("beta", 0.0);
25  param.Register_double("c_plaq", 0.0);
26  param.Register_double("c_rect", 0.0);
27 
28  param.Register_string("verbose_level", "NULL");
29  }
30 
31 
32 #ifdef USE_PARAMETERS_FACTORY
33  bool init_param = ParametersFactory::Register("Action.G_Rectangle", append_entry);
34 #endif
35 }
36 //- end
37 
38 //- parameters class
39 Parameters_Action_G_Rectangle::Parameters_Action_G_Rectangle() { append_entry(*this); }
40 //- end
41 
42 //====================================================================
43 void Action_G_Rectangle::set_parameters(const Parameters& params)
44 {
45  const string str_vlevel = params.get_string("verbose_level");
46 
47  m_vl = vout.set_verbose_level(str_vlevel);
48 
49  //- fetch and check input parameters
50  double beta, c_plaq, c_rect;
51 
52  int err = 0;
53  err += params.fetch_double("beta", beta);
54  err += params.fetch_double("c_plaq", c_plaq);
55  err += params.fetch_double("c_rect", c_rect);
56 
57  if (err) {
58  vout.crucial(m_vl, "Action_G_Rectangle: fetch error, input parameter not found.\n");
59  abort();
60  }
61 
62 
63  set_parameters(beta, c_plaq, c_rect);
64 }
65 
66 
67 //====================================================================
68 void Action_G_Rectangle::set_parameters(double beta,
69  double c_plaq, double c_rect)
70 {
71  //- print input parameters
72  vout.general(m_vl, "Action_G_Rectangle:\n");
73  vout.general(m_vl, " beta = %12.6f\n", beta);
74  vout.general(m_vl, " c_plaq = %12.6f\n", c_plaq);
75  vout.general(m_vl, " c_rect = %12.6f\n", c_rect);
76 
77  //- range check
78  // NB. beta,c_plaq,c_rect == 0 is allowed.
79 
80  //- store values
81  m_beta = beta;
82  m_c_plaq = c_plaq;
83  m_c_rect = c_rect;
84 
85  //- post-process
86  int Nc = CommonParameters::Nc();
87  int Nvol = CommonParameters::Nvol();
88  int Ndim = CommonParameters::Ndim();
89  int NinG = 2 * Nc * Nc;
90  m_force.reset(NinG, Nvol, Ndim);
91 
92  m_status_linkv = 0;
93 }
94 
95 
96 //====================================================================
97 double Action_G_Rectangle::langevin(RandomNumbers *rand)
98 {
99  m_status_linkv = 0;
100 
101  double H_U = calcH();
102 
103  return H_U;
104 }
105 
106 
107 //====================================================================
108 double Action_G_Rectangle::calcH()
109 {
110  int Nc = CommonParameters::Nc();
111  int Ndim = CommonParameters::Ndim();
112  int Nvol = CommonParameters::Nvol();
113  int Lvol = CommonParameters::Lvol();
114 
115  int Ndim2 = Ndim * (Ndim - 1) / 2;
116  int size_U = Lvol * Ndim2;
117 
118  Field_G Cup1(Nvol, 1), Cup2(Nvol, 1);
119  Field_G Cdn1(Nvol, 1), Cdn2(Nvol, 1);
120  Field_G Umu(Nvol, 1), Unu(Nvol, 1);
121  Field_G v(Nvol, 1), w(Nvol, 1), c(Nvol, 1);
122 
123  double plaqF = 0.0;
124  double rectF = 0.0;
125 
126  vout.general(m_vl, " Action_G_Rectangle: %s\n", m_label.c_str());
127 
128  for (int mu = 0; mu < Ndim; ++mu) {
129  for (int nu = mu + 1; nu < Ndim; ++nu) {
130  Cup1 = m_staple.upper(*m_U, mu, nu);
131  Cup2 = m_staple.upper(*m_U, nu, mu);
132 
133  // plaquette term
134  for (int site = 0; site < Nvol; ++site) {
135  plaqF += ReTr(m_U->mat(site, mu) * Cup1.mat_dag(site));
136  }
137 
138  // rectangular terms
139 
140  // +---+---+
141  // | | term
142  // x <---+
143 
144  Umu.setpart_ex(0, *m_U, mu);
145  Unu.setpart_ex(0, *m_U, nu);
146 
147  m_shift.backward(v, Cup2, mu);
148  m_shift.backward(c, Umu, nu);
149 
150  w.mult_Field_Gnd(0, c, 0, v, 0);
151  c.mult_Field_Gnn(0, Unu, 0, w, 0);
152  for (int site = 0; site < Nvol; ++site) {
153  rectF += ReTr(m_U->mat(site, mu) * c.mat_dag(site));
154  }
155 
156  // +---+
157  // | |
158  // + + term
159  // | |
160  // x v
161 
162  m_shift.backward(v, Unu, mu);
163  m_shift.backward(c, Cup1, nu);
164 
165  w.mult_Field_Gnd(0, c, 0, v, 0);
166  c.mult_Field_Gnn(0, Unu, 0, w, 0);
167  for (int site = 0; site < Nvol; ++site) {
168  rectF += ReTr(m_U->mat(site, mu) * c.mat_dag(site));
169  }
170  }
171  }
172 
173  plaqF = Communicator::reduce_sum(plaqF);
174  rectF = Communicator::reduce_sum(rectF);
175 
176  double plaq = plaqF / Nc;
177  vout.general(m_vl, " Plaquette = %18.8f\n", plaq / size_U);
178 
179  double H_U = m_c_plaq * (Ndim2 * Lvol - plaqF / Nc)
180  + m_c_rect * (Ndim2 * Lvol * 2 - rectF / Nc);
181 
182  H_U = m_beta * H_U;
183 
184  vout.general(m_vl, " H_Grectangle = %18.8f\n", H_U);
185  vout.general(m_vl, " H_G/dof = %18.8f\n", H_U / size_U);
186 
187  return H_U;
188 }
189 
190 
191 //====================================================================
192 const Field Action_G_Rectangle::force()
193 {
194  if (m_status_linkv == 0) {
195  int Nc = CommonParameters::Nc();
196  int Nvol = CommonParameters::Nvol();
197  int Ndim = CommonParameters::Ndim();
198 
199  assert(m_U->nin() == Nc * Nc * 2);
200  assert(m_U->nvol() == Nvol);
201  assert(m_U->nex() == Ndim);
202 
203  double betaNc = m_beta / Nc;
204  Mat_SU_N ut(Nc);
205 
206  Field_G force(Nvol, Ndim);
207  Field_G force1(Nvol, 1);
208 
209  Field_G Cup1(Nvol, 1), Cup2(Nvol, 1);
210  Field_G Cdn1(Nvol, 1), Cdn2(Nvol, 1);
211  Field_G Umu(Nvol, 1), Unu(Nvol, 1);
212  Field_G v(Nvol, 1), w(Nvol, 1), c(Nvol, 1);
213 
214  vout.general(m_vl, " Action_G_Rectangle: %s\n", m_label.c_str());
215 
216  for (int mu = 0; mu < Ndim; ++mu) {
217  force1 = 0.0;
218  for (int nu = 0; nu < Ndim; ++nu) {
219  if (nu == mu) continue;
220 
221  Cup1 = m_staple.upper(*m_U, mu, nu);
222  Cup2 = m_staple.upper(*m_U, nu, mu);
223  Cdn1 = m_staple.lower(*m_U, mu, nu);
224  Cdn2 = m_staple.lower(*m_U, nu, mu);
225 
226  // plaquette term
227  force1.addpart_ex(0, Cup1, 0, m_c_plaq);
228  force1.addpart_ex(0, Cdn1, 0, m_c_plaq);
229 
230  // rectangular term
231 
232  Umu.setpart_ex(0, *m_U, mu);
233  Unu.setpart_ex(0, *m_U, nu);
234 
235  // +---+---+
236  // | | term
237  // x <---+
238 
239  m_shift.backward(v, Cup2, mu);
240  m_shift.backward(c, Umu, nu);
241 
242  w.mult_Field_Gnd(0, c, 0, v, 0);
243  c.mult_Field_Gnn(0, Unu, 0, w, 0);
244 
245  force1.addpart_ex(0, c, 0, m_c_rect);
246 
247  // +---+
248  // | |
249  // + + term
250  // | |
251  // x v
252 
253  m_shift.backward(v, Unu, mu);
254  m_shift.backward(c, Cup1, nu);
255 
256  w.mult_Field_Gnd(0, c, 0, v, 0);
257  c.mult_Field_Gnn(0, Unu, 0, w, 0);
258 
259  force1.addpart_ex(0, c, 0, m_c_rect);
260 
261  // +---+---+
262  // | | term
263  // +---x v
264 
265  m_shift.backward(v, Unu, mu);
266  m_shift.backward(c, Umu, nu);
267 
268  w.mult_Field_Gnd(0, c, 0, v, 0);
269  c.mult_Field_Gnn(0, Cdn2, 0, w, 0);
270 
271  force1.addpart_ex(0, c, 0, m_c_rect);
272 
273  // x <---+
274  // | | term
275  // +---+---+
276 
277  m_shift.backward(v, Cup2, mu);
278 
279  w.mult_Field_Gnn(0, Umu, 0, v, 0);
280  v.mult_Field_Gdn(0, Unu, 0, w, 0);
281 
282  m_shift.forward(c, v, nu);
283 
284  force1.addpart_ex(0, c, 0, m_c_rect);
285 
286  // x ^
287  // | |
288  // + + term
289  // | |
290  // +---+
291 
292  m_shift.backward(v, Unu, mu);
293 
294  w.mult_Field_Gnn(0, Cdn1, 0, v, 0);
295  v.mult_Field_Gdn(0, Unu, 0, w, 0);
296 
297  m_shift.forward(c, v, nu);
298 
299  force1.addpart_ex(0, c, 0, m_c_rect);
300 
301  // +---x ^
302  // | | term
303  // +---+---+
304 
305  m_shift.backward(v, Unu, mu);
306 
307  w.mult_Field_Gnn(0, Umu, 0, v, 0);
308  v.mult_Field_Gdn(0, Cdn2, 0, w, 0);
309 
310  m_shift.forward(c, v, nu);
311 
312  force1.addpart_ex(0, c, 0, m_c_rect);
313  }
314 
315  force.mult_Field_Gnd(mu, *m_U, mu, force1, 0);
316  force.at_Field_G(mu);
317  }
318 
319  force *= -betaNc;
320 
321  m_force = (Field)force;
322  ++m_status_linkv;
323 
324  double Fave, Fmax, Fdev;
325  m_force.stat(Fave, Fmax, Fdev);
326  vout.general(m_vl, " Frectangle_ave = %12.6f Frectangle_max = %12.6f Frectangle_dev = %12.6f\n",
327  Fave, Fmax, Fdev);
328 
329  return m_force;
330  } else {
331  vout.general(m_vl, " Action_G_Rectangle returns previous force.\n");
332  return m_force;
333  }
334 }
335 
336 
337 //====================================================================
338 //============================================================END=====