Bridge++  Ver. 2.0.2
energyMomentumTensor.cpp
Go to the documentation of this file.
1 
14 #include "energyMomentumTensor.h"
15 
16 const std::string EnergyMomentumTensor::class_name = "EnergyMomentumTensor";
17 
18 //====================================================================
19 void EnergyMomentumTensor::set_parameters(const Parameters& params)
20 {
21  m_filename_output = params.get_string("filename_output");
22  if (m_filename_output.empty()) {
23  m_filename_output = "stdout";
24  }
25 
26  std::string vlevel;
27  if (!params.fetch_string("verbose_level", vlevel)) {
28  m_vl = vout.set_verbose_level(vlevel);
29  }
30 
31  //- fetch and check input parameters
32  double c_plaq, c_rect;
33  int max_mom;
34 
35  int err = 0;
36  err += params.fetch_double("c_plaq", c_plaq);
37  err += params.fetch_double("c_rect", c_rect);
38  err += params.fetch_int("max_momentum", max_mom);
39 
40  if (err) {
41  vout.crucial(m_vl, "Error at %s: input parameter not found.\n", class_name.c_str());
42  exit(EXIT_FAILURE);
43  }
44 
45 
46  set_parameters(c_plaq, c_rect, max_mom);
47 }
48 
49 
50 //====================================================================
51 void EnergyMomentumTensor::get_parameters(Parameters& params) const
52 {
53  params.set_double("c_plaq", m_c_plaq);
54  params.set_double("c_rect", m_c_rect);
55  params.set_int("max_momentum", m_max_mom);
56 
57  params.set_string("filename_output", m_filename_output);
58  params.set_string("verbose_level", vout.get_verbose_level(m_vl));
59 }
60 
61 
62 //====================================================================
63 void EnergyMomentumTensor::set_parameters(const double c_plaq, const double c_rect, const int max_mom)
64 {
65  //- print input parameters
66  vout.general(m_vl, "Topological Charge measurement:\n");
67  vout.general(m_vl, " c_plaq = %12.6f\n", c_plaq);
68  vout.general(m_vl, " c_rect = %12.6f\n", c_rect);
69  vout.general(m_vl, " max_momentum = %d\n", max_mom);
70 
71  //- range check
72  // NB. beta,c_plaq,c_rect == 0 is allowed.
73 
74  //- store values
75  m_c_plaq = c_plaq;
76  m_c_rect = c_rect;
77  m_max_mom = max_mom;
78 }
79 
80 
81 //====================================================================
82 double EnergyMomentumTensor::measure_EMT(const double t_flow)
83 {
84  assert(m_flag_field_set == 1);
85 
86  const int Ndim = CommonParameters::Ndim();
87  const int Nvol = CommonParameters::Nvol();
88  const int NPE = CommonParameters::NPE();
89 
90  static const double l_c_rect = 1.0 / 8.0;
91  m_c_rect = -1.0 / 12.0;
92  m_c_plaq = 1.0 - 8 * m_c_rect;
93 
94  //- output
95  std::ostream& log_file_previous = vout.getStream();
96  std::ofstream log_file;
97 
98  if (m_filename_output != "stdout") {
99  log_file.open(m_filename_output.c_str(), std::ios::app);
100  vout.init(log_file);
101  }
102 
103  double O2_1x1 = 0.0;
104  double O2_1x2 = 0.0;
105  double O2_plaq = 0.0;
106  for (int mu = 0; mu < 6; ++mu) {
107  O2_1x1 += m_field_strength.contract(m_Fmunu_1x1[mu], m_Fmunu_1x1[mu]);
108  O2_1x2 += m_field_strength.contract(m_Fmunu_1x2[mu], m_Fmunu_1x2[mu]);
109  O2_plaq += m_field_strength.contract(m_Fmunu_plaq[mu], m_Fmunu_plaq[mu]);
110  }
111  O2_1x1 *= 4.0 / Nvol / NPE;
112  vout.general(m_vl, " O2_clover_plaq = %.8f %.16e\n", t_flow, O2_1x1);
113  O2_1x2 *= 8.0 / Nvol / NPE;
114  const double O2_imp = (m_c_plaq * O2_1x1 + m_c_rect * O2_1x2);
115  vout.general(m_vl, " O2_clover_imp = %.8f %.16e\n", t_flow, O2_imp);
116  const double O2_rect = l_c_rect * O2_1x2;
117  vout.general(m_vl, " O2_clover_rect = %.8f %.16e\n", t_flow, O2_rect);
118  O2_plaq *= 4.0 / Nvol / NPE;
119  vout.general(m_vl, " O2_plaq = %.8f %.16e\n", t_flow, O2_plaq);
120 
121  for (int mu = 0; mu < Ndim; ++mu) {
122  for (int nu = 0; nu < Ndim; ++nu) {
123  double O1_1x1 = 0.0;
124  double O1_1x2 = 0.0;
125  double O1_plaq = 0.0;
126  for (int rho = 0; rho < Ndim; rho++) {
127  if ((mu != rho) && (nu != rho)) {
128  int ind1 = index_munu2i(mu, rho);
129  int ind2 = index_munu2i(nu, rho);
130  int fac1 = factor(mu, rho);
131  int fac2 = factor(nu, rho);
132 
133  double scr = m_field_strength.contract(m_Fmunu_1x1[ind1], m_Fmunu_1x1[ind2]);
134  scr *= fac1 * fac2;
135  O1_1x1 += scr;
136 
137  scr = m_field_strength.contract(m_Fmunu_1x2[ind1], m_Fmunu_1x2[ind2]);
138  scr *= fac1 * fac2;
139  O1_1x2 += scr;
140 
141  scr = m_field_strength.contract(m_Fmunu_plaq[ind1], m_Fmunu_plaq[ind2]);
142  scr *= fac1 * fac2;
143  O1_plaq += scr;
144  }
145  }
146  O1_1x1 *= 2.0 / Nvol / NPE;
147  O1_1x2 *= 4.0 / Nvol / NPE;
148  O1_plaq *= 2.0 / Nvol / NPE;
149  vout.general(m_vl, " O1_clover_plaq = %d %d %.8f %.16e\n", mu, nu, t_flow, O1_1x1);
150  double O1_rect = l_c_rect * O1_1x2;
151  vout.general(m_vl, " O1_clover_1x2 = %d %d %.8f %.16e\n", mu, nu, t_flow, O1_1x2);
152  vout.general(m_vl, " O1_clover_rect = %d %d %.8f %.16e\n", mu, nu, t_flow, O1_rect);
153  double O1_imp = (m_c_plaq * O1_1x1 + m_c_rect * O1_1x2);
154  vout.general(m_vl, " O1_clover_imp = %d %d %.8f %.16e\n", mu, nu, t_flow, O1_imp);
155  vout.general(m_vl, " O1_plaq = %d %d %.8f %.16e\n", mu, nu, t_flow, O1_plaq);
156  }
157  }
158 
159  if (m_filename_output != "stdout") {
160  log_file.close();
161  vout.init(log_file_previous);
162  }
163  return O2_imp;
164 }
165 
166 
167 //====================================================================
168 double EnergyMomentumTensor::measure_EMT_at_t(const double t_flow)
169 {
170  assert(m_flag_field_set == 1);
171 
172  const int Lt = CommonParameters::Lt();
173 
174  const int Nvol = CommonParameters::Nvol();
175  const int NPE = CommonParameters::NPE();
176 
177  const int Ndim = CommonParameters::Ndim();
178 
179  const double normalization = double(Lt) / Nvol / NPE;
180  double result = 0.0; // dummy initialization
181 
182  static const double l_c_rect = 1.0 / 8.0;
183  m_c_rect = -1.0 / 12.0;
184  m_c_plaq = 1.0 - 8 * m_c_rect;
185 
186  //- output
187  std::ostream& log_file_previous = vout.getStream();
188  std::ofstream log_file;
189 
190  if (m_filename_output != "stdout") {
191  log_file.open(m_filename_output.c_str(), std::ios::app);
192  vout.init(log_file);
193  }
194 
195  for (int mu = 0; mu < Ndim; ++mu) {
196  for (int nu = 0; nu < Ndim; ++nu) {
197  std::vector<double> corr_plaq(Lt, 0);
198  std::vector<double> corr_1x1(Lt, 0);
199  std::vector<double> corr_1x2(Lt, 0);
200  std::vector<double> corr_scr(Lt);
201  for (int rho = 0; rho < Ndim; rho++) {
202  if ((mu != rho) && (nu != rho)) {
203  int ind1 = index_munu2i(mu, rho);
204  int ind2 = index_munu2i(nu, rho);
205  int fac1 = factor(mu, rho);
206  int fac2 = factor(nu, rho);
207  m_field_strength.contract_at_t(corr_scr, m_Fmunu_plaq[ind1], m_Fmunu_plaq[ind2]);
208  for (int t = 0; t < Lt; ++t) {
209  corr_plaq[t] += (fac1 * fac2 * corr_scr[t]);
210  }
211  m_field_strength.contract_at_t(corr_scr, m_Fmunu_1x1[ind1], m_Fmunu_1x1[ind2]);
212  for (int t = 0; t < Lt; ++t) {
213  corr_1x1[t] += (fac1 * fac2 * corr_scr[t]);
214  }
215  m_field_strength.contract_at_t(corr_scr, m_Fmunu_1x2[ind1], m_Fmunu_1x2[ind2]);
216  for (int t = 0; t < Lt; ++t) {
217  corr_1x2[t] += (fac1 * fac2 * corr_scr[t]);
218  }
219  }
220  }
221  //double O1_plaq = 0.0;
222  //double O1_1x1 = 0.0;
223  //double O1_1x2 = 0.0;
224  for (int t = 0; t < Lt; ++t) {
225  corr_plaq[t] *= 2 * normalization;
226  //O1_plaq += corr_plaq[t];
227  corr_1x1[t] *= 2 * normalization;
228  //O1_1x1 += corr_1x1[t];
229  corr_1x2[t] *= 4 * normalization;
230  //O1_1x2 += corr_1x2[t];
231 
232  vout.general(m_vl, " O1_clover_plaq_t = %d %d %.8f %d %.16e\n", mu, nu, t_flow, t, corr_1x1[t]);
233  double O1_rect = l_c_rect * corr_1x2[t];
234  vout.general(m_vl, " O1_clover_rect_t = %d %d %.8f %d %.16e\n", mu, nu, t_flow, t, O1_rect);
235  double O1_imp = (m_c_plaq * corr_1x1[t] + m_c_rect * corr_1x2[t]);
236  vout.general(m_vl, " O1_clover_imp_t = %d %d %.8f %d %.16e\n", mu, nu, t_flow, t, O1_imp);
237  result = O1_imp;
238  vout.general(m_vl, " O1_plaq_t = %d %d %.8f %d %.16e\n", mu, nu, t_flow, t, corr_plaq[t]);
239  }
240 
241  /*
242  double O1_imp = m_c_plaq * O1_1x1 + m_c_rect * O1_1x2;
243  double O1_rect = l_c_rect * O1_1x2;
244  vout.general(m_vl, " O1_clover_plaq_sumt = %d %d %.8f %.16e\n", mu, nu, t_flow, O1_1x1/Lt);
245  //vout.general(m_vl, " O1_clover_1x2_sumt = %d %d %.8f %.16e\n", mu, nu, t_flow, O1_1x2/Lt);
246  vout.general(m_vl, " O1_clover_rect_sumt = %d %d %.8f %.16e\n", mu, nu, t_flow, O1_rect/Lt);
247  vout.general(m_vl, " O1_clover_imp_sumt = %d %d %.8f %.16e\n", mu, nu, t_flow, O1_imp/Lt);
248  vout.general(m_vl, " O1_plaq_sumt = %d %d %.8f %.16e\n", mu, nu, t_flow, O1_plaq/Lt);
249  */
250  }
251  }
252 
253  if (m_filename_output != "stdout") {
254  log_file.close();
255  vout.init(log_file_previous);
256  }
257  return result;
258 }
259 
260 
261 //====================================================================
262 double EnergyMomentumTensor::measure_EMT_at_t_FT(const double t_flow)
263 {
264  assert(m_flag_field_set == 1);
265 
266  const int Lt = CommonParameters::Lt();
267 
268  const int Nvol = CommonParameters::Nvol();
269  const int NPE = CommonParameters::NPE();
270 
271  const int Ndim = CommonParameters::Ndim();
272 
273  const double normalization = double(Lt) / Nvol / NPE;
274  double result = 0.0; // dummy initialization
275 
276  static const double l_c_rect = 1.0 / 8.0;
277  m_c_rect = -1.0 / 12.0;
278  m_c_plaq = 1.0 - 8 * m_c_rect;
279 
280  const int Np = (2 * m_max_mom + 1);
281  vector<int> source_position(4, 0);
282  vector<int> momentum_sink(3);
283 
284  //- output
285  std::ostream& log_file_previous = vout.getStream();
286  std::ofstream log_file;
287 
288  if (m_filename_output != "stdout") {
289  log_file.open(m_filename_output.c_str(), std::ios::app);
290  vout.init(log_file);
291  }
292 
293  for (int ipx = 0; ipx < m_max_mom + 1; ipx++) {
294  for (int ipy = 0; ipy < Np; ipy++) {
295  for (int ipz = 0; ipz < Np; ipz++) {
296  momentum_sink[0] = ipx;
297  momentum_sink[1] = ipy - m_max_mom;
298  momentum_sink[2] = ipz - m_max_mom;
299  for (int mu = 0; mu < Ndim; ++mu) {
300  for (int nu = 0; nu < Ndim; ++nu) {
301  std::vector<double> corr_plaq(Lt, 0);
302  std::vector<double> corr_1x1(Lt, 0);
303  std::vector<double> corr_1x2(Lt, 0);
304  std::vector<double> corr_scr(Lt);
305  for (int rho = 0; rho < Ndim; rho++) {
306  if ((mu != rho) && (nu != rho)) {
307  int ind1 = index_munu2i(mu, rho);
308  int ind2 = index_munu2i(nu, rho);
309  int fac1 = factor(mu, rho);
310  int fac2 = factor(nu, rho);
311  m_field_strength.contract_at_t(corr_scr, m_Fmunu_plaq[ind1], m_Fmunu_plaq[ind2], momentum_sink, source_position);
312  for (int t = 0; t < Lt; ++t) {
313  corr_plaq[t] += (fac1 * fac2 * corr_scr[t]);
314  }
315  m_field_strength.contract_at_t(corr_scr, m_Fmunu_1x1[ind1], m_Fmunu_1x1[ind2], momentum_sink, source_position);
316  for (int t = 0; t < Lt; ++t) {
317  corr_1x1[t] += (fac1 * fac2 * corr_scr[t]);
318  }
319  m_field_strength.contract_at_t(corr_scr, m_Fmunu_1x2[ind1], m_Fmunu_1x2[ind2], momentum_sink, source_position);
320  for (int t = 0; t < Lt; ++t) {
321  corr_1x2[t] += (fac1 * fac2 * corr_scr[t]);
322  }
323  }
324  }
325  //double O1_plaq = 0.0;
326  // double O1_1x1 = 0.0;
327  // double O1_1x2 = 0.0;
328  for (int t = 0; t < Lt; ++t) {
329  corr_plaq[t] *= 2 * normalization;
330  //O1_plaq += corr_plaq[t];
331  corr_1x1[t] *= 2 * normalization;
332  //O1_1x1 += corr_1x1[t];
333  corr_1x2[t] *= 4 * normalization;
334  //O1_1x2 += corr_1x2[t];
335  vout.general(m_vl, " O1_clover_plaq_t_FT = %d %d %d %d %d %.8f %d %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, t, corr_1x1[t]);
336  double O1_rect = l_c_rect * corr_1x2[t];
337  vout.general(m_vl, " O1_clover_rect_t_FT = %d %d %d %d %d %.8f %d %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, t, O1_rect);
338  double O1_imp = (m_c_plaq * corr_1x1[t] + m_c_rect * corr_1x2[t]);
339  vout.general(m_vl, " O1_clover_imp_t_FT = %d %d %d %d %d %.8f %d %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, t, O1_imp);
340  result = O1_imp;
341  vout.general(m_vl, " O1_plaq_t_FT = %d %d %d %d %d %.8f %d %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, t, corr_plaq[t]);
342  }
343 
344  /*
345  double O1_imp = m_c_plaq * O1_1x1 + m_c_rect * O1_1x2;
346  double O1_rect = l_c_rect * O1_1x2;
347  vout.general(m_vl, " O1_clover_plaq_sumt_FT = %d %d %d %d %d %.8f %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, O1_1x1/Lt);
348  //vout.general(m_vl, " O1_clover_1x2_sumt_FT = %d %d %.8f %.16e\n", mu, nu, t_flow, O1_1x2/Lt);
349  vout.general(m_vl, " O1_clover_rect_sumt_FT = %d %d %d %d %d %.8f %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, O1_rect/Lt);
350  vout.general(m_vl, " O1_clover_imp_sumt_FT = %d %d %d %d %d %.8f %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, O1_imp/Lt);
351  vout.general(m_vl, " O1_plaq_sumt_FT = %d %d %d %d %d %.8f %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, O1_plaq/Lt);
352  */
353  }
354  }
355  }
356  }
357  }
358 
359  if (m_filename_output != "stdout") {
360  log_file.close();
361  vout.init(log_file_previous);
362  }
363  return result;
364 }
365 
366 
367 //====================================================================
368 double EnergyMomentumTensor::measure_EMT_at_x(const double t_flow)
369 {
370  assert(m_flag_field_set == 1);
371 
372  const int Lx = CommonParameters::Lx();
373 
374  const int Nvol = CommonParameters::Nvol();
375  const int NPE = CommonParameters::NPE();
376 
377  const int Ndim = CommonParameters::Ndim();
378 
379  const double normalization = double(Lx) / Nvol / NPE;
380  double result = 0.0; // dummy initialization
381 
382  static const double l_c_rect = 1.0 / 8.0;
383  m_c_rect = -1.0 / 12.0;
384  m_c_plaq = 1.0 - 8 * m_c_rect;
385 
386  std::ostream& log_file_previous = vout.getStream();
387  std::ofstream log_file;
388 
389  if (m_filename_output != "stdout") {
390  log_file.open(m_filename_output.c_str(), std::ios::app);
391  vout.init(log_file);
392  }
393 
394  for (int mu = 0; mu < Ndim; ++mu) {
395  for (int nu = 0; nu < Ndim; ++nu) {
396  std::vector<double> corr_plaq(Lx, 0);
397  std::vector<double> corr_1x1(Lx, 0);
398  std::vector<double> corr_1x2(Lx, 0);
399  std::vector<double> corr_scr(Lx);
400  for (int rho = 0; rho < Ndim; rho++) {
401  if ((mu != rho) && (nu != rho)) {
402  int ind1 = index_munu2i(mu, rho);
403  int ind2 = index_munu2i(nu, rho);
404  int fac1 = factor(mu, rho);
405  int fac2 = factor(nu, rho);
406  m_field_strength.contract_at_x(corr_scr, m_Fmunu_plaq[ind1], m_Fmunu_plaq[ind2]);
407  for (int x = 0; x < Lx; ++x) {
408  corr_plaq[x] += (fac1 * fac2 * corr_scr[x]);
409  }
410  m_field_strength.contract_at_x(corr_scr, m_Fmunu_1x1[ind1], m_Fmunu_1x1[ind2]);
411  for (int x = 0; x < Lx; ++x) {
412  corr_1x1[x] += (fac1 * fac2 * corr_scr[x]);
413  }
414  m_field_strength.contract_at_x(corr_scr, m_Fmunu_1x2[ind1], m_Fmunu_1x2[ind2]);
415  for (int x = 0; x < Lx; ++x) {
416  corr_1x2[x] += (fac1 * fac2 * corr_scr[x]);
417  }
418  }
419  }
420  //double O1_plaq = 0.0;
421  //double O1_1x1 = 0.0;
422  //double O1_1x2 = 0.0;
423  for (int x = 0; x < Lx; ++x) {
424  corr_plaq[x] *= 2 * normalization;
425  //O1_plaq += corr_plaq[x];
426  corr_1x1[x] *= 2 * normalization;
427  //O1_1x1 += corr_1x1[x];
428  corr_1x2[x] *= 4 * normalization;
429  //O1_1x2 += corr_1x2[x];
430 
431  vout.general(m_vl, " O1_clover_plaq_x = %d %d %.8f %d %.16e\n", mu, nu, t_flow, x, corr_1x1[x]);
432  double O1_rect = l_c_rect * corr_1x2[x];
433  vout.general(m_vl, " O1_clover_rect_x = %d %d %.8f %d %.16e\n", mu, nu, t_flow, x, O1_rect);
434  double O1_imp = (m_c_plaq * corr_1x1[x] + m_c_rect * corr_1x2[x]);
435  vout.general(m_vl, " O1_clover_imp_x = %d %d %.8f %d %.16e\n", mu, nu, t_flow, x, O1_imp);
436  result = O1_imp;
437  vout.general(m_vl, " O1_plaq_x = %d %d %.8f %d %.16e\n", mu, nu, t_flow, x, corr_plaq[x]);
438  }
439 
440  /*
441  double O1_imp = m_c_plaq * O1_1x1 + m_c_rect * O1_1x2;
442  double O1_rect = l_c_rect * O1_1x2;
443  vout.general(m_vl, " O1_clover_plaq_sumx = %d %d %.8f %.16e\n", mu, nu, t_flow, O1_1x1/Lx);
444  vout.general(m_vl, " O1_clover_rect_sumx = %d %d %.8f %.16e\n", mu, nu, t_flow, O1_rect/Lx);
445  vout.general(m_vl, " O1_clover_imp_sumx = %d %d %.8f %.16e\n", mu, nu, t_flow, O1_imp/Lx);
446  vout.general(m_vl, " O1_plaq_sumx = %d %d %.8f %.16e\n", mu, nu, t_flow, O1_plaq/Lx);
447  */
448  }
449  }
450 
451  if (m_filename_output != "stdout") {
452  log_file.close();
453  vout.init(log_file_previous);
454  }
455  return result;
456 }
457 
458 
459 //====================================================================
460 double EnergyMomentumTensor::measure_EMT_at_x_FT(const double t_flow)
461 {
462  assert(m_flag_field_set == 1);
463 
464  const int Lx = CommonParameters::Lx();
465 
466  const int Nvol = CommonParameters::Nvol();
467  const int NPE = CommonParameters::NPE();
468 
469  const int Ndim = CommonParameters::Ndim();
470 
471  const double normalization = double(Lx) / Nvol / NPE;
472  double result = 0.0; // dummy initialization
473 
474  static const double l_c_rect = 1.0 / 8.0;
475  m_c_rect = -1.0 / 12.0;
476  m_c_plaq = 1.0 - 8 * m_c_rect;
477 
478  const int Np = (2 * m_max_mom + 1);
479  vector<int> source_position(4, 0);
480  vector<int> momentum_sink(3);
481 
482  //- output
483  std::ostream& log_file_previous = vout.getStream();
484  std::ofstream log_file;
485 
486  if (m_filename_output != "stdout") {
487  log_file.open(m_filename_output.c_str(), std::ios::app);
488  vout.init(log_file);
489  }
490 
491  for (int ipx = 0; ipx < m_max_mom + 1; ipx++) {
492  for (int ipy = 0; ipy < Np; ipy++) {
493  for (int ipz = 0; ipz < Np; ipz++) {
494  momentum_sink[0] = ipx;
495  momentum_sink[1] = ipy - m_max_mom;
496  momentum_sink[2] = ipz - m_max_mom;
497  for (int mu = 0; mu < Ndim; ++mu) {
498  for (int nu = 0; nu < Ndim; ++nu) {
499  std::vector<double> corr_plaq(Lx, 0);
500  std::vector<double> corr_1x1(Lx, 0);
501  std::vector<double> corr_1x2(Lx, 0);
502  std::vector<double> corr_scr(Lx);
503  for (int rho = 0; rho < Ndim; rho++) {
504  if ((mu != rho) && (nu != rho)) {
505  int ind1 = index_munu2i(mu, rho);
506  int ind2 = index_munu2i(nu, rho);
507  int fac1 = factor(mu, rho);
508  int fac2 = factor(nu, rho);
509  m_field_strength.contract_at_x(corr_scr, m_Fmunu_plaq[ind1], m_Fmunu_plaq[ind2], momentum_sink, source_position);
510  for (int t = 0; t < Lx; ++t) {
511  corr_plaq[t] += (fac1 * fac2 * corr_scr[t]);
512  }
513  m_field_strength.contract_at_x(corr_scr, m_Fmunu_1x1[ind1], m_Fmunu_1x1[ind2], momentum_sink, source_position);
514  for (int t = 0; t < Lx; ++t) {
515  corr_1x1[t] += (fac1 * fac2 * corr_scr[t]);
516  }
517  m_field_strength.contract_at_x(corr_scr, m_Fmunu_1x2[ind1], m_Fmunu_1x2[ind2], momentum_sink, source_position);
518  for (int t = 0; t < Lx; ++t) {
519  corr_1x2[t] += (fac1 * fac2 * corr_scr[t]);
520  }
521  }
522  }
523  //double O1_plaq = 0.0;
524  //double O1_1x1 = 0.0;
525  //double O1_1x2 = 0.0;
526  for (int t = 0; t < Lx; ++t) {
527  corr_plaq[t] *= 2 * normalization;
528  //O1_plaq += corr_plaq[t];
529  corr_1x1[t] *= 2 * normalization;
530  //O1_1x1 += corr_1x1[t];
531  corr_1x2[t] *= 4 * normalization;
532  //O1_1x2 += corr_1x2[t];
533  vout.general(m_vl, " O1_clover_plaq_x_FT = %d %d %d %d %d %.8f %d %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, t, corr_1x1[t]);
534  double O1_rect = l_c_rect * corr_1x2[t];
535  vout.general(m_vl, " O1_clover_rect_x_FT = %d %d %d %d %d %.8f %d %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, t, O1_rect);
536  double O1_imp = (m_c_plaq * corr_1x1[t] + m_c_rect * corr_1x2[t]);
537  vout.general(m_vl, " O1_clover_imp_x_FT = %d %d %d %d %d %.8f %d %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, t, O1_imp);
538  result = O1_imp;
539  vout.general(m_vl, " O1_plaq_x_FT = %d %d %d %d %d %.8f %d %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, t, corr_plaq[t]);
540  }
541 
542  /*
543  double O1_imp = m_c_plaq * O1_1x1 + m_c_rect * O1_1x2;
544  double O1_rect = l_c_rect * O1_1x2;
545  vout.general(m_vl, " O1_clover_plaq_sumx_FT = %d %d %d %d %d %.8f %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, O1_1x1/Lx);
546  //vout.general(m_vl, " O1_clover_1x2_sumx_FT = %d %d %.8f %.16e\n", mu, nu, t_flow, O1_1x2/Lx);
547  vout.general(m_vl, " O1_clover_rect_sumx_FT = %d %d %d %d %d %.8f %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, O1_rect/Lx);
548  vout.general(m_vl, " O1_clover_imp_sumx_FT = %d %d %d %d %d %.8f %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, O1_imp/Lx);
549  vout.general(m_vl, " O1_plaq_sumx_FT = %d %d %d %d %d %.8f %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, O1_plaq/Lx);
550  */
551  }
552  }
553  }
554  }
555  }
556 
557  if (m_filename_output != "stdout") {
558  log_file.close();
559  vout.init(log_file_previous);
560  }
561  return result;
562 }
563 
564 
565 //====================================================================
566 double EnergyMomentumTensor::measure_EMT_at_y(const double t_flow)
567 {
568  assert(m_flag_field_set == 1);
569 
570  const int Ly = CommonParameters::Ly();
571 
572  const int Nvol = CommonParameters::Nvol();
573  const int NPE = CommonParameters::NPE();
574 
575  const int Ndim = CommonParameters::Ndim();
576 
577  const double normalization = double(Ly) / Nvol / NPE;
578  double result = 0.0; // dummy initialization
579 
580  static const double l_c_rect = 1.0 / 8.0;
581  m_c_rect = -1.0 / 12.0;
582  m_c_plaq = 1.0 - 8 * m_c_rect;
583 
584  //- output
585  std::ostream& log_file_previous = vout.getStream();
586  std::ofstream log_file;
587 
588  if (m_filename_output != "stdout") {
589  log_file.open(m_filename_output.c_str(), std::ios::app);
590  vout.init(log_file);
591  }
592 
593  for (int mu = 0; mu < Ndim; ++mu) {
594  for (int nu = 0; nu < Ndim; ++nu) {
595  std::vector<double> corr_plaq(Ly, 0);
596  std::vector<double> corr_1x1(Ly, 0);
597  std::vector<double> corr_1x2(Ly, 0);
598  std::vector<double> corr_scr(Ly);
599  for (int rho = 0; rho < Ndim; rho++) {
600  if ((mu != rho) && (nu != rho)) {
601  int ind1 = index_munu2i(mu, rho);
602  int ind2 = index_munu2i(nu, rho);
603  int fac1 = factor(mu, rho);
604  int fac2 = factor(nu, rho);
605  m_field_strength.contract_at_y(corr_scr, m_Fmunu_plaq[ind1], m_Fmunu_plaq[ind2]);
606  for (int y = 0; y < Ly; ++y) {
607  corr_plaq[y] += (fac1 * fac2 * corr_scr[y]);
608  }
609  m_field_strength.contract_at_y(corr_scr, m_Fmunu_1x1[ind1], m_Fmunu_1x1[ind2]);
610  for (int y = 0; y < Ly; ++y) {
611  corr_1x1[y] += (fac1 * fac2 * corr_scr[y]);
612  }
613  m_field_strength.contract_at_y(corr_scr, m_Fmunu_1x2[ind1], m_Fmunu_1x2[ind2]);
614  for (int y = 0; y < Ly; ++y) {
615  corr_1x2[y] += (fac1 * fac2 * corr_scr[y]);
616  }
617  }
618  }
619  //double O1_plaq = 0.0;
620  //double O1_1x1 = 0.0;
621  //double O1_1x2 = 0.0;
622  for (int y = 0; y < Ly; ++y) {
623  corr_plaq[y] *= 2 * normalization;
624  //O1_plaq += corr_plaq[y];
625  corr_1x1[y] *= 2 * normalization;
626  //O1_1x1 += corr_1x1[y];
627  corr_1x2[y] *= 4 * normalization;
628  //O1_1x2 += corr_1x2[y];
629 
630  vout.general(m_vl, " O1_clover_plaq_y = %d %d %.8f %d %.16e\n", mu, nu, t_flow, y, corr_1x1[y]);
631  double O1_rect = l_c_rect * corr_1x2[y];
632  vout.general(m_vl, " O1_clover_rect_y = %d %d %.8f %d %.16e\n", mu, nu, t_flow, y, O1_rect);
633  double O1_imp = (m_c_plaq * corr_1x1[y] + m_c_rect * corr_1x2[y]);
634  vout.general(m_vl, " O1_clover_imp_y = %d %d %.8f %d %.16e\n", mu, nu, t_flow, y, O1_imp);
635  result = O1_imp;
636  vout.general(m_vl, " O1_plaq_y = %d %d %.8f %d %.16e\n", mu, nu, t_flow, y, corr_plaq[y]);
637  }
638 
639  /*
640  double O1_imp = m_c_plaq * O1_1x1 + m_c_rect * O1_1x2;
641  double O1_rect = l_c_rect * O1_1x2;
642  vout.general(m_vl, " O1_clover_plaq_sumy = %d %d %.8f %.16e\n", mu, nu, t_flow, O1_1x1/Ly);
643  vout.general(m_vl, " O1_clover_rect_sumy = %d %d %.8f %.16e\n", mu, nu, t_flow, O1_rect/Ly);
644  vout.general(m_vl, " O1_clover_imp_sumy = %d %d %.8f %.16e\n", mu, nu, t_flow, O1_imp/Ly);
645  vout.general(m_vl, " O1_plaq_sumy = %d %d %.8f %.16e\n", mu, nu, t_flow, O1_plaq/Ly);
646  */
647  }
648  }
649 
650  if (m_filename_output != "stdout") {
651  log_file.close();
652  vout.init(log_file_previous);
653  }
654  return result;
655 }
656 
657 
658 //====================================================================
659 double EnergyMomentumTensor::measure_EMT_at_y_FT(const double t_flow)
660 {
661  assert(m_flag_field_set == 1);
662 
663  const int Ly = CommonParameters::Ly();
664 
665  const int Nvol = CommonParameters::Nvol();
666  const int NPE = CommonParameters::NPE();
667 
668  const int Ndim = CommonParameters::Ndim();
669 
670  const double normalization = double(Ly) / Nvol / NPE;
671  double result = 0.0; // dummy initialization
672 
673  static const double l_c_rect = 1.0 / 8.0;
674  m_c_rect = -1.0 / 12.0;
675  m_c_plaq = 1.0 - 8 * m_c_rect;
676 
677  const int Np = (2 * m_max_mom + 1);
678  vector<int> source_position(4, 0);
679  vector<int> momentum_sink(3);
680 
681  //- output
682  std::ostream& log_file_previous = vout.getStream();
683  std::ofstream log_file;
684 
685  if (m_filename_output != "stdout") {
686  log_file.open(m_filename_output.c_str(), std::ios::app);
687  vout.init(log_file);
688  }
689 
690  for (int ipx = 0; ipx < m_max_mom + 1; ipx++) {
691  for (int ipy = 0; ipy < Np; ipy++) {
692  for (int ipz = 0; ipz < Np; ipz++) {
693  momentum_sink[0] = ipx;
694  momentum_sink[1] = ipy - m_max_mom;
695  momentum_sink[2] = ipz - m_max_mom;
696  for (int mu = 0; mu < Ndim; ++mu) {
697  for (int nu = 0; nu < Ndim; ++nu) {
698  std::vector<double> corr_plaq(Ly, 0);
699  std::vector<double> corr_1x1(Ly, 0);
700  std::vector<double> corr_1x2(Ly, 0);
701  std::vector<double> corr_scr(Ly);
702  for (int rho = 0; rho < Ndim; rho++) {
703  if ((mu != rho) && (nu != rho)) {
704  int ind1 = index_munu2i(mu, rho);
705  int ind2 = index_munu2i(nu, rho);
706  int fac1 = factor(mu, rho);
707  int fac2 = factor(nu, rho);
708  m_field_strength.contract_at_y(corr_scr, m_Fmunu_plaq[ind1], m_Fmunu_plaq[ind2], momentum_sink, source_position);
709  for (int t = 0; t < Ly; ++t) {
710  corr_plaq[t] += (fac1 * fac2 * corr_scr[t]);
711  }
712  m_field_strength.contract_at_y(corr_scr, m_Fmunu_1x1[ind1], m_Fmunu_1x1[ind2], momentum_sink, source_position);
713  for (int t = 0; t < Ly; ++t) {
714  corr_1x1[t] += (fac1 * fac2 * corr_scr[t]);
715  }
716  m_field_strength.contract_at_y(corr_scr, m_Fmunu_1x2[ind1], m_Fmunu_1x2[ind2], momentum_sink, source_position);
717  for (int t = 0; t < Ly; ++t) {
718  corr_1x2[t] += (fac1 * fac2 * corr_scr[t]);
719  }
720  }
721  }
722  //double O1_plaq = 0.0;
723  //double O1_1x1 = 0.0;
724  //double O1_1x2 = 0.0;
725  for (int t = 0; t < Ly; ++t) {
726  corr_plaq[t] *= 2 * normalization;
727  //O1_plaq += corr_plaq[t];
728  corr_1x1[t] *= 2 * normalization;
729  //O1_1x1 += corr_1x1[t];
730  corr_1x2[t] *= 4 * normalization;
731  //O1_1x2 += corr_1x2[t];
732  vout.general(m_vl, " O1_clover_plaq_y_FT = %d %d %d %d %d %.8f %d %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, t, corr_1x1[t]);
733  double O1_rect = l_c_rect * corr_1x2[t];
734  vout.general(m_vl, " O1_clover_rect_y_FT = %d %d %d %d %d %.8f %d %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, t, O1_rect);
735  double O1_imp = (m_c_plaq * corr_1x1[t] + m_c_rect * corr_1x2[t]);
736  vout.general(m_vl, " O1_clover_imp_y_FT = %d %d %d %d %d %.8f %d %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, t, O1_imp);
737  result = O1_imp;
738  vout.general(m_vl, " O1_plaq_y_FT = %d %d %d %d %d %.8f %d %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, t, corr_plaq[t]);
739  }
740 
741  /*
742  double O1_imp = m_c_plaq * O1_1x1 + m_c_rect * O1_1x2;
743  double O1_rect = l_c_rect * O1_1x2;
744  vout.general(m_vl, " O1_clover_plaq_sumy_FT = %d %d %d %d %d %.8f %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, O1_1x1/Ly);
745  //vout.general(m_vl, " O1_clover_1x2_sumy_FT = %d %d %.8f %.16e\n", mu, nu, t_flow, O1_1x2/Ly);
746  vout.general(m_vl, " O1_clover_rect_sumy_FT = %d %d %d %d %d %.8f %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, O1_rect/Ly);
747  vout.general(m_vl, " O1_clover_imp_sumy_FT = %d %d %d %d %d %.8f %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, O1_imp/Ly);
748  vout.general(m_vl, " O1_plaq_sumy_FT = %d %d %d %d %d %.8f %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, O1_plaq/Ly);
749  */
750  }
751  }
752  }
753  }
754  }
755 
756  if (m_filename_output != "stdout") {
757  log_file.close();
758  vout.init(log_file_previous);
759  }
760  return result;
761 }
762 
763 
764 //====================================================================
765 double EnergyMomentumTensor::measure_EMT_at_z(const double t_flow)
766 {
767  assert(m_flag_field_set == 1);
768 
769  const int Lz = CommonParameters::Lz();
770 
771  const int Nvol = CommonParameters::Nvol();
772  const int NPE = CommonParameters::NPE();
773 
774  const int Ndim = CommonParameters::Ndim();
775 
776  const double normalization = double(Lz) / Nvol / NPE;
777  double result = 0.0; // dummy initialization
778 
779  static const double l_c_rect = 1.0 / 8.0;
780  m_c_rect = -1.0 / 12.0;
781  m_c_plaq = 1.0 - 8 * m_c_rect;
782 
783  //- output
784  std::ostream& log_file_previous = vout.getStream();
785  std::ofstream log_file;
786 
787  if (m_filename_output != "stdout") {
788  log_file.open(m_filename_output.c_str(), std::ios::app);
789  vout.init(log_file);
790  }
791 
792  for (int mu = 0; mu < Ndim; ++mu) {
793  for (int nu = 0; nu < Ndim; ++nu) {
794  std::vector<double> corr_plaq(Lz, 0);
795  std::vector<double> corr_1x1(Lz, 0);
796  std::vector<double> corr_1x2(Lz, 0);
797  std::vector<double> corr_scr(Lz);
798  for (int rho = 0; rho < Ndim; rho++) {
799  if ((mu != rho) && (nu != rho)) {
800  int ind1 = index_munu2i(mu, rho);
801  int ind2 = index_munu2i(nu, rho);
802  int fac1 = factor(mu, rho);
803  int fac2 = factor(nu, rho);
804  m_field_strength.contract_at_z(corr_scr, m_Fmunu_plaq[ind1], m_Fmunu_plaq[ind2]);
805  for (int z = 0; z < Lz; ++z) {
806  corr_plaq[z] += (fac1 * fac2 * corr_scr[z]);
807  }
808  m_field_strength.contract_at_z(corr_scr, m_Fmunu_1x1[ind1], m_Fmunu_1x1[ind2]);
809  for (int z = 0; z < Lz; ++z) {
810  corr_1x1[z] += (fac1 * fac2 * corr_scr[z]);
811  }
812  m_field_strength.contract_at_z(corr_scr, m_Fmunu_1x2[ind1], m_Fmunu_1x2[ind2]);
813  for (int z = 0; z < Lz; ++z) {
814  corr_1x2[z] += (fac1 * fac2 * corr_scr[z]);
815  }
816  }
817  }
818  //double O1_plaq = 0.0;
819  //double O1_1x1 = 0.0;
820  //double O1_1x2 = 0.0;
821  for (int z = 0; z < Lz; ++z) {
822  corr_plaq[z] *= 2 * normalization;
823  //O1_plaq += corr_plaq[z];
824  corr_1x1[z] *= 2 * normalization;
825  //O1_1x1 += corr_1x1[z];
826  corr_1x2[z] *= 4 * normalization;
827  //O1_1x2 += corr_1x2[z];
828 
829  vout.general(m_vl, " O1_clover_plaq_z = %d %d %.8f %d %.16e\n", mu, nu, t_flow, z, corr_1x1[z]);
830  double O1_rect = l_c_rect * corr_1x2[z];
831  vout.general(m_vl, " O1_clover_rect_z = %d %d %.8f %d %.16e\n", mu, nu, t_flow, z, O1_rect);
832  double O1_imp = (m_c_plaq * corr_1x1[z] + m_c_rect * corr_1x2[z]);
833  vout.general(m_vl, " O1_clover_imp_z = %d %d %.8f %d %.16e\n", mu, nu, t_flow, z, O1_imp);
834  result = O1_imp;
835  vout.general(m_vl, " O1_plaq_z = %d %d %.8f %d %.16e\n", mu, nu, t_flow, z, corr_plaq[z]);
836  }
837 
838  /*
839  double O1_imp = m_c_plaq * O1_1x1 + m_c_rect * O1_1x2;
840  double O1_rect = l_c_rect * O1_1x2;
841  vout.general(m_vl, " O1_clover_plaq_sumz = %d %d %.8f %.16e\n", mu, nu, t_flow, O1_1x1/Lz);
842  vout.general(m_vl, " O1_clover_rect_sumz = %d %d %.8f %.16e\n", mu, nu, t_flow, O1_rect/Lz);
843  vout.general(m_vl, " O1_clover_imp_sumz = %d %d %.8f %.16e\n", mu, nu, t_flow, O1_imp/Lz);
844  vout.general(m_vl, " O1_plaq_sumz = %d %d %.8f %.16e\n", mu, nu, t_flow, O1_plaq/Lz);
845  */
846  }
847  }
848 
849  if (m_filename_output != "stdout") {
850  log_file.close();
851  vout.init(log_file_previous);
852  }
853  return result;
854 }
855 
856 
857 //====================================================================
858 double EnergyMomentumTensor::measure_EMT_at_z_FT(const double t_flow)
859 {
860  assert(m_flag_field_set == 1);
861 
862  const int Lz = CommonParameters::Lz();
863 
864  const int Nvol = CommonParameters::Nvol();
865  const int NPE = CommonParameters::NPE();
866 
867  const int Ndim = CommonParameters::Ndim();
868 
869  const double normalization = double(Lz) / Nvol / NPE;
870  double result = 0.0; // dummy initialization
871 
872  static const double l_c_rect = 1.0 / 8.0;
873  m_c_rect = -1.0 / 12.0;
874  m_c_plaq = 1.0 - 8 * m_c_rect;
875 
876  const int Np = (2 * m_max_mom + 1);
877  vector<int> source_position(4, 0);
878  vector<int> momentum_sink(3);
879 
880  //- output
881  std::ostream& log_file_previous = vout.getStream();
882  std::ofstream log_file;
883 
884  if (m_filename_output != "stdout") {
885  log_file.open(m_filename_output.c_str(), std::ios::app);
886  vout.init(log_file);
887  }
888 
889  for (int ipx = 0; ipx < m_max_mom + 1; ipx++) {
890  for (int ipy = 0; ipy < Np; ipy++) {
891  for (int ipz = 0; ipz < Np; ipz++) {
892  momentum_sink[0] = ipx;
893  momentum_sink[1] = ipy - m_max_mom;
894  momentum_sink[2] = ipz - m_max_mom;
895  for (int mu = 0; mu < Ndim; ++mu) {
896  for (int nu = 0; nu < Ndim; ++nu) {
897  std::vector<double> corr_plaq(Lz, 0);
898  std::vector<double> corr_1x1(Lz, 0);
899  std::vector<double> corr_1x2(Lz, 0);
900  std::vector<double> corr_scr(Lz);
901  for (int rho = 0; rho < Ndim; rho++) {
902  if ((mu != rho) && (nu != rho)) {
903  int ind1 = index_munu2i(mu, rho);
904  int ind2 = index_munu2i(nu, rho);
905  int fac1 = factor(mu, rho);
906  int fac2 = factor(nu, rho);
907  m_field_strength.contract_at_z(corr_scr, m_Fmunu_plaq[ind1], m_Fmunu_plaq[ind2], momentum_sink, source_position);
908  for (int t = 0; t < Lz; ++t) {
909  corr_plaq[t] += (fac1 * fac2 * corr_scr[t]);
910  }
911  m_field_strength.contract_at_z(corr_scr, m_Fmunu_1x1[ind1], m_Fmunu_1x1[ind2], momentum_sink, source_position);
912  for (int t = 0; t < Lz; ++t) {
913  corr_1x1[t] += (fac1 * fac2 * corr_scr[t]);
914  }
915  m_field_strength.contract_at_z(corr_scr, m_Fmunu_1x2[ind1], m_Fmunu_1x2[ind2], momentum_sink, source_position);
916  for (int t = 0; t < Lz; ++t) {
917  corr_1x2[t] += (fac1 * fac2 * corr_scr[t]);
918  }
919  }
920  }
921  //double O1_plaq = 0.0;
922  //double O1_1x1 = 0.0;
923  //double O1_1x2 = 0.0;
924  for (int t = 0; t < Lz; ++t) {
925  corr_plaq[t] *= 2 * normalization;
926  //O1_plaq += corr_plaq[t];
927  corr_1x1[t] *= 2 * normalization;
928  //O1_1x1 += corr_1x1[t];
929  corr_1x2[t] *= 4 * normalization;
930  //O1_1x2 += corr_1x2[t];
931  vout.general(m_vl, " O1_clover_plaq_z_FT = %d %d %d %d %d %.8f %d %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, t, corr_1x1[t]);
932  double O1_rect = l_c_rect * corr_1x2[t];
933  vout.general(m_vl, " O1_clover_rect_z_FT = %d %d %d %d %d %.8f %d %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, t, O1_rect);
934  double O1_imp = (m_c_plaq * corr_1x1[t] + m_c_rect * corr_1x2[t]);
935  vout.general(m_vl, " O1_clover_imp_z_FT = %d %d %d %d %d %.8f %d %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, t, O1_imp);
936  result = O1_imp;
937  vout.general(m_vl, " O1_plaq_z_FT = %d %d %d %d %d %.8f %d %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, t, corr_plaq[t]);
938  }
939 
940  /*
941  double O1_imp = m_c_plaq * O1_1x1 + m_c_rect * O1_1x2;
942  double O1_rect = l_c_rect * O1_1x2;
943  vout.general(m_vl, " O1_clover_plaq_sumz_FT = %d %d %d %d %d %.8f %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, O1_1x1/Lz);
944  //vout.general(m_vl, " O1_clover_1x2_sumz_FT = %d %d %.8f %.16e\n", mu, nu, t_flow, O1_1x2/Lz);
945  vout.general(m_vl, " O1_clover_rect_sumz_FT = %d %d %d %d %d %.8f %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, O1_rect/Lz);
946  vout.general(m_vl, " O1_clover_imp_sumz_FT = %d %d %d %d %d %.8f %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, O1_imp/Lz);
947  vout.general(m_vl, " O1_plaq_sumz_FT = %d %d %d %d %d %.8f %.16e\n", momentum_sink[0], momentum_sink[1], momentum_sink[2], mu, nu, t_flow, O1_plaq/Lz);
948  */
949  }
950  }
951  }
952  }
953  }
954 
955  if (m_filename_output != "stdout") {
956  log_file.close();
957  vout.init(log_file_previous);
958  }
959  return result;
960 }
961 
962 
963 //====================================================================
964 void EnergyMomentumTensor::set_field_strength(const Field_G& U)
965 {
966  const int Ndim = CommonParameters::Ndim();
967 
968  // NB. #(mu,nu)=6 i.e. (1,2),(1,3),(1,4),(2,3),(2,4),(3,4)
969  m_Fmunu_plaq.resize(6);
970  m_Fmunu_1x1.resize(6);
971  m_Fmunu_1x2.resize(6);
972  int i_munu = 0;
973  for (int mu = 0; mu < Ndim; ++mu) {
974  for (int nu = mu + 1; nu < Ndim; ++nu) {
975  m_field_strength.construct_Fmunu_plaq_traceless(m_Fmunu_plaq[i_munu], mu, nu, U);
976  m_field_strength.construct_Fmunu_1x1_traceless(m_Fmunu_1x1[i_munu], mu, nu, U);
977  m_field_strength.construct_Fmunu_1x2_traceless(m_Fmunu_1x2[i_munu], mu, nu, U);
978  ++i_munu;
979  }
980  }
981  m_flag_field_set = 1;
982 }
983 
984 
985 //====================================================================
986 int EnergyMomentumTensor::factor(const int mu, const int nu)
987 {
988  if (mu < nu) {
989  return 1;
990  } else if (mu > nu) {
991  return -1;
992  } else {
993  return 0;
994  }
995 }
996 
997 
998 //====================================================================
999 int EnergyMomentumTensor::index_munu2i(int mu, int nu)
1000 {
1001  const int Ndim = CommonParameters::Ndim();
1002 
1003  assert(mu < Ndim);
1004  assert(nu < Ndim);
1005 
1006  if (mu > nu) {
1007  int scr = mu;
1008 
1009  mu = nu;
1010  nu = scr;
1011  }
1012  if ((mu == 0) && (nu == 1)) return 0;
1013  else if ((mu == 0) && (nu == 2)) return 1;
1014  else if ((mu == 0) && (nu == 3)) return 2;
1015  else if ((mu == 1) && (nu == 2)) return 3;
1016  else if ((mu == 1) && (nu == 3)) return 4;
1017  else if ((mu == 2) && (nu == 3)) return 5;
1018  else return 0;
1019 }
1020 
1021 
1022 //====================================================================
1023 //============================================================END=====
EnergyMomentumTensor::m_vl
Bridge::VerboseLevel m_vl
Definition: energyMomentumTensor.h:45
energyMomentumTensor.h
EnergyMomentumTensor::measure_EMT_at_z
double measure_EMT_at_z(const double t_flow)
Measure energy momentum tensor density in z direction and print out the result using an argument t_f...
Definition: energyMomentumTensor.cpp:765
EnergyMomentumTensor::m_field_strength
FieldStrength m_field_strength
Definition: energyMomentumTensor.h:55
Bridge::BridgeIO::init
void init(const std::string &filename)
Definition: bridgeIO.cpp:53
Parameters::set_string
void set_string(const string &key, const string &value)
Definition: parameters.cpp:39
EnergyMomentumTensor::measure_EMT_at_x_FT
double measure_EMT_at_x_FT(const double t_flow)
Definition: energyMomentumTensor.cpp:460
EnergyMomentumTensor::m_Fmunu_plaq
std::vector< Field_G > m_Fmunu_plaq
Definition: energyMomentumTensor.h:57
CommonParameters::Ndim
static int Ndim()
Definition: commonParameters.h:117
EnergyMomentumTensor::measure_EMT_at_t
double measure_EMT_at_t(const double t_flow)
Measure energy momentum tensor density as a function of time and print out the result using an argum...
Definition: energyMomentumTensor.cpp:168
Parameters
Class for parameters.
Definition: parameters.h:46
Parameters::set_double
void set_double(const string &key, const double value)
Definition: parameters.cpp:33
EnergyMomentumTensor::index_munu2i
int index_munu2i(int mu, int nu)
Returns array number [0-5] of vector m_Fmunu correponding to (mu,nu).
Definition: energyMomentumTensor.cpp:999
EnergyMomentumTensor::m_Fmunu_1x1
std::vector< Field_G > m_Fmunu_1x1
Definition: energyMomentumTensor.h:58
CommonParameters::Ly
static int Ly()
Definition: commonParameters.h:92
CommonParameters::Nvol
static int Nvol()
Definition: commonParameters.h:109
EnergyMomentumTensor::m_c_plaq
double m_c_plaq
Definition: energyMomentumTensor.h:50
EnergyMomentumTensor::set_field_strength
void set_field_strength(const Field_G &U)
Construct the anti-Hermitian traceless field strength by the flowed link U. Should be called before ...
Definition: energyMomentumTensor.cpp:964
FieldStrength::construct_Fmunu_1x2_traceless
void construct_Fmunu_1x2_traceless(Field_G &Fmunu, const int mu, const int nu, const Field_G &U)
Constructs the anti-Hermitian traceless field strength with eight 1x2 rectangular clover leaves.
Definition: fieldStrength.cpp:395
FieldStrength::contract
double contract(const Field_G &Fmunu_1, const Field_G &Fmunu_2)
Calculate and returns its value. Intended to be used for the topological charge and energy momentum ...
Definition: fieldStrength.cpp:662
FieldStrength::contract_at_x
void contract_at_x(std::vector< double > &corr_global, const Field_G &Fmunu_1, const Field_G &Fmunu_2)
Calculate corr[x]= . Intended to be used to calculate correlations function of the topological charge...
Definition: fieldStrength.cpp:762
EnergyMomentumTensor::measure_EMT_at_y_FT
double measure_EMT_at_y_FT(const double t_flow)
Definition: energyMomentumTensor.cpp:659
EnergyMomentumTensor::measure_EMT_at_y
double measure_EMT_at_y(const double t_flow)
Measure energy momentum tensor density in y direction and print out the result using an argument t_f...
Definition: energyMomentumTensor.cpp:566
FieldStrength::construct_Fmunu_plaq_traceless
void construct_Fmunu_plaq_traceless(Field_G &Fmunu, const int mu, const int nu, const Field_G &U)
Constructs the anti-Hermitian traceless field strength with an imaginary part of the plaquette.
Definition: fieldStrength.cpp:630
FieldStrength::contract_at_z
void contract_at_z(std::vector< double > &corr_global, const Field_G &Fmunu_1, const Field_G &Fmunu_2)
Calculate corr[z]= . Intended to be used to calculate correlations function of the topological charge...
Definition: fieldStrength.cpp:932
CommonParameters::Lt
static int Lt()
Definition: commonParameters.h:94
CommonParameters::Lx
static int Lx()
Definition: commonParameters.h:91
CommonParameters::Lz
static int Lz()
Definition: commonParameters.h:93
FieldStrength::contract_at_y
void contract_at_y(std::vector< double > &corr_global, const Field_G &Fmunu_1, const Field_G &Fmunu_2)
Calculate corr[y]= . Intended to be used to calculate correlations function of the topological charge...
Definition: fieldStrength.cpp:847
EnergyMomentumTensor::set_parameters
void set_parameters(const Parameters &params)
setting parameters.
Definition: energyMomentumTensor.cpp:19
EnergyMomentumTensor::m_Fmunu_1x2
std::vector< Field_G > m_Fmunu_1x2
Definition: energyMomentumTensor.h:59
EnergyMomentumTensor::m_max_mom
int m_max_mom
maximum of momentum for Fourier transformation: p_x=[0,max_mom], p_y=[-max_mom,max_mom],...
Definition: energyMomentumTensor.h:53
CommonParameters::NPE
static int NPE()
Definition: commonParameters.h:101
EnergyMomentumTensor::m_c_rect
double m_c_rect
Definition: energyMomentumTensor.h:51
EnergyMomentumTensor::measure_EMT_at_x
double measure_EMT_at_x(const double t_flow)
Measure energy momentum tensor density in x direction and print out the result using an argument t_f...
Definition: energyMomentumTensor.cpp:368
Bridge::BridgeIO::set_verbose_level
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:133
FieldStrength::contract_at_t
void contract_at_t(std::vector< double > &corr_global, const Field_G &Fmunu_1, const Field_G &Fmunu_2)
Calculate corr[t]= . Intended to be used to calculate correlations function of the topological charge...
Definition: fieldStrength.cpp:678
Parameters::set_int
void set_int(const string &key, const int value)
Definition: parameters.cpp:36
EnergyMomentumTensor::measure_EMT_at_t_FT
double measure_EMT_at_t_FT(const double t_flow)
Definition: energyMomentumTensor.cpp:262
FieldStrength::construct_Fmunu_1x1_traceless
void construct_Fmunu_1x1_traceless(Field_G &Fmunu, const int mu, const int nu, const Field_G &U)
Constructs the anti-Hermitian traceless field strength with four 1x1 plquette clover leaves.
Definition: fieldStrength.cpp:324
Parameters::fetch_string
int fetch_string(const string &key, string &value) const
Definition: parameters.cpp:378
Parameters::fetch_double
int fetch_double(const string &key, double &value) const
Definition: parameters.cpp:327
EnergyMomentumTensor::class_name
static const std::string class_name
Definition: energyMomentumTensor.h:42
EnergyMomentumTensor::factor
int factor(const int mu, const int nu)
Returns +1 for mu<nu, -1 for mu>nu and 0 otherwise.
Definition: energyMomentumTensor.cpp:986
Parameters::get_string
string get_string(const string &key) const
Definition: parameters.cpp:221
EnergyMomentumTensor::m_filename_output
std::string m_filename_output
Definition: energyMomentumTensor.h:48
Bridge::BridgeIO::crucial
void crucial(const char *format,...)
Definition: bridgeIO.cpp:180
EnergyMomentumTensor::measure_EMT
double measure_EMT(const double t_flow)
Definition: energyMomentumTensor.cpp:82
Field_G
SU(N) gauge field.
Definition: field_G.h:38
Bridge::BridgeIO::getStream
std::ostream & getStream()
Definition: bridgeIO.cpp:396
EnergyMomentumTensor::m_flag_field_set
int m_flag_field_set
Definition: energyMomentumTensor.h:60
EnergyMomentumTensor::measure_EMT_at_z_FT
double measure_EMT_at_z_FT(const double t_flow)
Definition: energyMomentumTensor.cpp:858
Parameters::fetch_int
int fetch_int(const string &key, int &value) const
Definition: parameters.cpp:346
Bridge::BridgeIO::general
void general(const char *format,...)
Definition: bridgeIO.cpp:200
EnergyMomentumTensor::get_parameters
void get_parameters(Parameters &params) const
Definition: energyMomentumTensor.cpp:51
Bridge::vout
BridgeIO vout
Definition: bridgeIO.cpp:512
Bridge::BridgeIO::get_verbose_level
static std::string get_verbose_level(const VerboseLevel vl)
Definition: bridgeIO.cpp:154