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