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