Bridge++  Ver. 2.0.2
corr4pt_4spinor.cpp
Go to the documentation of this file.
1 
14 #include "corr4pt_4spinor.h"
15 
16 const std::string Corr4pt_4spinor::class_name = "Corr4pt_4spinor";
17 
18 //====================================================================
19 void Corr4pt_4spinor::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 
32 
33 //====================================================================
34 void Corr4pt_4spinor::get_parameters(Parameters& params) const
35 {
36  params.set_string("filename_output", m_filename_output);
37  params.set_string("verbose_level", vout.get_verbose_level(m_vl));
38 }
39 
40 
41 //====================================================================
42 double Corr4pt_4spinor::meson_all(const std::vector<Field_F>& sq1,
43  const std::vector<Field_F>& sq2,
44  const std::vector<Field_F>& sq3,
45  const std::vector<Field_F>& sq4)
46 {
47  const int Lt = CommonParameters::Lt();
48 
49  std::vector<dcomplex> corr_global(Lt);
50  GammaMatrix gm_sink_12, gm_sink_34, gm_src_12, gm_src_34;
51 
52  std::ostream& log_file_previous = vout.getStream();
53  std::ofstream log_file;
54 
55  if (m_filename_output != "stdout") {
56  log_file.open(m_filename_output.c_str(), std::ios::app);
57  vout.init(log_file);
58  }
59 
60 
61  gm_src_12 = m_gmset->get_GM(m_gmset->GAMMA5);
62  gm_src_34 = m_gmset->get_GM(m_gmset->GAMMA5);
63  gm_sink_12 = m_gmset->get_GM(m_gmset->GAMMA5);
64  gm_sink_34 = m_gmset->get_GM(m_gmset->GAMMA5);
65  vout.general(m_vl, "PS <-- PS 4pt_correlator:\n");
66  meson_correlator(corr_global, gm_sink_12, gm_sink_34, gm_src_12, gm_src_34, sq1, sq2, sq3, sq4);
67  double result = real(corr_global[0]);
68 
69 #if 0
70  gm_src_12 = m_gmset->get_GM(m_gmset->GAMMA1);
71  gm_src_34 = m_gmset->get_GM(m_gmset->GAMMA1);
72  gm_sink_12 = m_gmset->get_GM(m_gmset->GAMMA1);
73  gm_sink_34 = m_gmset->get_GM(m_gmset->GAMMA1);
74  vout.general(m_vl, "V1 <-- V1 4pt_correlator:\n");
75  meson_correlator(corr_global, gm_sink_12, gm_sink_34, gm_src_12, gm_src_34, sq1, sq2, sq3, sq4);
76 
77  gm_src_12 = m_gmset->get_GM(m_gmset->GAMMA2);
78  gm_src_34 = m_gmset->get_GM(m_gmset->GAMMA2);
79  gm_sink_12 = m_gmset->get_GM(m_gmset->GAMMA2);
80  gm_sink_34 = m_gmset->get_GM(m_gmset->GAMMA2);
81  vout.general(m_vl, "V2 <-- V2 4pt_correlator:\n");
82  meson_correlator(corr_global, gm_sink_12, gm_sink_34, gm_src_12, gm_src_34, sq1, sq2, sq3, sq4);
83 
84  gm_src_12 = m_gmset->get_GM(m_gmset->GAMMA3);
85  gm_src_34 = m_gmset->get_GM(m_gmset->GAMMA3);
86  gm_sink_12 = m_gmset->get_GM(m_gmset->GAMMA3);
87  gm_sink_34 = m_gmset->get_GM(m_gmset->GAMMA3);
88  vout.general(m_vl, "V3 <-- V3 4pt_correlator:\n");
89  meson_correlator(corr_global, gm_sink_12, gm_sink_34, gm_src_12, gm_src_34, sq1, sq2, sq3, sq4);
90 
91  gm_src_12 = m_gmset->get_GM(m_gmset->GAMMA5);
92  gm_src_34 = m_gmset->get_GM(m_gmset->GAMMA5);
93  gm_sink_12 = m_gmset->get_GM(m_gmset->GAMMA45);
94  gm_sink_34 = m_gmset->get_GM(m_gmset->GAMMA45);
95  vout.general(m_vl, "A4 <-- PS 4pt_correlator:\n");
96  meson_correlator(corr_global, gm_sink_12, gm_sink_34, gm_src_12, gm_src_34, sq1, sq2, sq3, sq4);
97 
98  gm_src_12 = m_gmset->get_GM(m_gmset->GAMMA54);
99  gm_src_34 = m_gmset->get_GM(m_gmset->GAMMA54);
100  gm_sink_12 = m_gmset->get_GM(m_gmset->GAMMA5);
101  gm_sink_34 = m_gmset->get_GM(m_gmset->GAMMA5);
102  vout.general(m_vl, "PS <-- A4 4pt_correlator:\n");
103  meson_correlator(corr_global, gm_sink_12, gm_sink_34, gm_src_12, gm_src_34, sq1, sq2, sq3, sq4);
104 
105  gm_src_12 = m_gmset->get_GM(m_gmset->UNITY);
106  gm_src_34 = m_gmset->get_GM(m_gmset->UNITY);
107  gm_sink_12 = m_gmset->get_GM(m_gmset->UNITY);
108  gm_sink_34 = m_gmset->get_GM(m_gmset->UNITY);
109  vout.general(m_vl, "S <-- S 4pt_correlator:\n");
110  meson_correlator(corr_global, gm_sink_12, gm_sink_34, gm_src_12, gm_src_34, sq1, sq2, sq3, sq4);
111 
112  gm_src_12 = m_gmset->get_GM(m_gmset->GAMMA51);
113  gm_src_34 = m_gmset->get_GM(m_gmset->GAMMA51);
114  gm_sink_12 = m_gmset->get_GM(m_gmset->GAMMA51);
115  gm_sink_34 = m_gmset->get_GM(m_gmset->GAMMA51);
116  vout.general(m_vl, "GAMMA51 <-- GAMMA51 4pt_correlator:\n");
117  meson_correlator(corr_global, gm_sink_12, gm_sink_34, gm_src_12, gm_src_34, sq1, sq2, sq3, sq4);
118 
119  gm_src_12 = m_gmset->get_GM(m_gmset->GAMMA52);
120  gm_src_34 = m_gmset->get_GM(m_gmset->GAMMA52);
121  gm_sink_12 = m_gmset->get_GM(m_gmset->GAMMA52);
122  gm_sink_34 = m_gmset->get_GM(m_gmset->GAMMA52);
123  vout.general(m_vl, "GAMMA52 <-- GAMMA52 4pt_correlator:\n");
124  meson_correlator(corr_global, gm_sink_12, gm_sink_34, gm_src_12, gm_src_34, sq1, sq2, sq3, sq4);
125 
126  gm_src_12 = m_gmset->get_GM(m_gmset->GAMMA53);
127  gm_src_34 = m_gmset->get_GM(m_gmset->GAMMA53);
128  gm_sink_12 = m_gmset->get_GM(m_gmset->GAMMA53);
129  gm_sink_34 = m_gmset->get_GM(m_gmset->GAMMA53);
130  vout.general(m_vl, "GAMMA53 <-- GAMMA53 4pt_correlator:\n");
131  meson_correlator(corr_global, gm_sink_12, gm_sink_34, gm_src_12, gm_src_34, sq1, sq2, sq3, sq4);
132 
133  gm_src_12 = m_gmset->get_GM(m_gmset->SIGMA12);
134  gm_src_34 = m_gmset->get_GM(m_gmset->SIGMA12);
135  gm_sink_12 = m_gmset->get_GM(m_gmset->SIGMA12);
136  gm_sink_34 = m_gmset->get_GM(m_gmset->SIGMA12);
137  vout.general(m_vl, "SIGMA12 <-- SIGMA12 4pt_correlator:\n");
138  meson_correlator(corr_global, gm_sink_12, gm_sink_34, gm_src_12, gm_src_34, sq1, sq2, sq3, sq4);
139 
140  gm_src_12 = m_gmset->get_GM(m_gmset->SIGMA23);
141  gm_src_34 = m_gmset->get_GM(m_gmset->SIGMA23);
142  gm_sink_12 = m_gmset->get_GM(m_gmset->SIGMA23);
143  gm_sink_34 = m_gmset->get_GM(m_gmset->SIGMA23);
144  vout.general(m_vl, "SIGMA23 <-- SIGMA23 4pt_correlator:\n");
145  meson_correlator(corr_global, gm_sink_12, gm_sink_34, gm_src_12, gm_src_34, sq1, sq2, sq3, sq4);
146 
147  gm_src_12 = m_gmset->get_GM(m_gmset->SIGMA31);
148  gm_src_34 = m_gmset->get_GM(m_gmset->SIGMA31);
149  gm_sink_12 = m_gmset->get_GM(m_gmset->SIGMA31);
150  gm_sink_34 = m_gmset->get_GM(m_gmset->SIGMA31);
151  vout.general(m_vl, "SIGMA31 <-- SIGMA31 4pt_correlator:\n");
152  meson_correlator(corr_global, gm_sink_12, gm_sink_34, gm_src_12, gm_src_34, sq1, sq2, sq3, sq4);
153 #endif
154 
155  if (m_filename_output != "stdout") {
156  log_file.close();
157  vout.init(log_file_previous);
158  }
159 
160  return result;
161 }
162 
163 
164 //====================================================================
165 void Corr4pt_4spinor::meson_correlator(std::vector<dcomplex>& corr_global,
166  const GammaMatrix& gm_sink_12,
167  const GammaMatrix& gm_sink_34,
168  const GammaMatrix& gm_src_12,
169  const GammaMatrix& gm_src_34,
170  const std::vector<Field_F>& sq1,
171  const std::vector<Field_F>& sq2,
172  const std::vector<Field_F>& sq3,
173  const std::vector<Field_F>& sq4)
174 {
175  const int Nc = CommonParameters::Nc();
176  const int Nd = CommonParameters::Nd();
177  const int Nx = CommonParameters::Nx();
178  const int Ny = CommonParameters::Ny();
179  const int Nz = CommonParameters::Nz();
180  const int Nt = CommonParameters::Nt();
181  const int Nvol = CommonParameters::Nvol();
182 
183  const int Lt = CommonParameters::Lt();
184 
185  assert(corr_global.size() == Lt);
186 
187  //- M_{ij}(x,t) := \bar{q_i}^c \Gamma_ij q_j^c
188  // Prop_ij(t) := q_i(t_sink) \bar{q_j}(t_src)
189  //
190  //- corr_direct = tr(Prop_11(t) Prop_22^{\dagger}(t)) tr(Prop_33(t) Prop_44^{\dagger}(t))
191  // +tr(Prop_13(t) Prop_24^{\dagger}(t)) tr(Prop_31(t) Prop_42^{\dagger}(t))
192  // x x' x x'
193  // t_sink 1 2 3 4 1 2 3 4 1 2
194  // | | | | / / / / / /
195  // ^ v ^ v + v ^ v ^ v
196  // | | | | / / / /
197  // t_src 1 2 3 4 1 2 3 4
198  // c0 c1 c0 c1
199  //
200  //- corr_cross = Prop_11(t)_{c0,c2} Prop_24^{\dagger}(t)_{c2,c1} Prop_33(t)_{c1,c3} Prop_42^{\dagger}(t)_{c3,c0}
201  // +Prop_31(t)_{c0,c3} Prop_44^{\dagger}(t)_{c3,c1} Prop_13(t)_{c1,c2} Prop_22^{\dagger}(t)_{c2,c0}
202  // = Prop_1124(t)_{c0,c1} Prop_3342(t)_{c1,c0}
203  // +Prop_3144(t)_{c0,c1} Prop_1322(t)_{c1,c0}
204  // x x' x x'
205  // t_sink 1 2 3 4 1 2 3 4
206  // | \|/ \|/ |
207  // ^ * + * v
208  // | /|\ /|\ |
209  // t_src 1 2 3 4 1 2 3 4
210  // c0 c1 c0 c1
211 
212  const GammaMatrix gm5 = m_gmset->get_GM(m_gmset->GAMMA5);
213 
214  const GammaMatrix gm_gm5_src_12 = gm_src_12.mult(gm5);
215  const GammaMatrix gm_gm5_src_34 = gm_src_34.mult(gm5);
216 
217  const GammaMatrix gm5_gm_sink_12 = gm5.mult(gm_sink_12);
218  const GammaMatrix gm5_gm_sink_34 = gm5.mult(gm_sink_34);
219 
220  //- corr_direct = tr(Prop_11(t) Prop_22^{\dagger}(t)) tr(Prop_33(t) Prop_44^{\dagger}(t))
221  // +tr(Prop_13(t) Prop_24^{\dagger}(t)) tr(Prop_31(t) Prop_42^{\dagger}(t))
222 
223  std::vector<dcomplex> corr_direct_1122_global(Lt);
224  std::vector<dcomplex> corr_direct_3344_global(Lt);
225  std::vector<dcomplex> corr_direct_1324_global(Lt);
226  std::vector<dcomplex> corr_direct_3142_global(Lt);
227 
228  corr_direct(corr_direct_1122_global,
229  gm5_gm_sink_12, gm_gm5_src_12,
230  sq1, sq2);
231  corr_direct(corr_direct_3344_global,
232  gm5_gm_sink_34, gm_gm5_src_34,
233  sq3, sq4);
234  corr_direct(corr_direct_1324_global,
235  gm5_gm_sink_12, gm_gm5_src_34,
236  sq3, sq4);
237  corr_direct(corr_direct_3142_global,
238  gm5_gm_sink_34, gm_gm5_src_12,
239  sq1, sq2);
240 
241  std::vector<dcomplex> corr_direct_global(Lt);
242  for (int t_global = 0; t_global < Lt; ++t_global) {
243  corr_direct_global[t_global]
244  = corr_direct_1122_global[t_global] * corr_direct_3344_global[t_global]
245  + corr_direct_1324_global[t_global] * corr_direct_3142_global[t_global];
246  }
247 
248  //- corr_cross = Prop_11(t)_{c0,c2} Prop_24^{\dagger}(t)_{c2,c1} Prop_33(t)_{c1,c3} Prop_42^{\dagger}(t)_{c3,c0}
249  // +Prop_31(t)_{c0,c3} Prop_44^{\dagger}(t)_{c3,c1} Prop_13(t)_{c1,c2} Prop_22^{\dagger}(t)_{c2,c0}
250  // = Prop_1124(t)_{c0,c1} Prop_3342(t)_{c1,c0}
251  // +Prop_3144(t)_{c0,c1} Prop_1322(t)_{c1,c0}
252 
253  typedef std::vector<dcomplex> CorrSet;
254  std::vector<CorrSet> corr_cross_1124_global(Lt);
255  std::vector<CorrSet> corr_cross_3342_global(Lt);
256  std::vector<CorrSet> corr_cross_1322_global(Lt);
257  std::vector<CorrSet> corr_cross_3144_global(Lt);
258  for (int t_global = 0; t_global < Lt; ++t_global) {
259  corr_cross_1124_global[t_global].resize(Nc * Nd * Nc * Nd);
260  corr_cross_3342_global[t_global].resize(Nc * Nd * Nc * Nd);
261  corr_cross_1322_global[t_global].resize(Nc * Nd * Nc * Nd);
262  corr_cross_3144_global[t_global].resize(Nc * Nd * Nc * Nd);
263  }
264 
265  corr_cross_sink(corr_cross_1124_global,
266  gm5_gm_sink_12, gm_gm5_src_12,
267  sq1, sq4);
268  corr_cross_sink(corr_cross_3342_global,
269  gm5_gm_sink_34, gm_gm5_src_34,
270  sq3, sq2);
271  corr_cross_sink(corr_cross_3144_global,
272  gm5_gm_sink_34, gm_gm5_src_12,
273  sq1, sq4);
274  corr_cross_sink(corr_cross_1322_global,
275  gm5_gm_sink_12, gm_gm5_src_34,
276  sq3, sq2);
277 
278  //- sum up at src
279  std::vector<dcomplex> corr_cross_global(Lt, cmplx(0.0, 0.0));
280  for (int t_global = 0; t_global < Lt; ++t_global) {
281  for (int c0 = 0; c0 < Nc; ++c0) {
282  for (int d0 = 0; d0 < Nd; ++d0) {
283  for (int c1 = 0; c1 < Nc; ++c1) {
284  for (int d1 = 0; d1 < Nd; ++d1) {
285  int i_cd2_01 = c0 + Nc * d0 + Nc * Nd * c1 + Nc * Nd * Nc * d1;
286  int i_cd2_10 = c1 + Nc * d1 + Nc * Nd * c0 + Nc * Nd * Nc * d0;
287 
288  dcomplex corr_1124 = corr_cross_1124_global[t_global][i_cd2_01];
289  dcomplex corr_3342 = corr_cross_3342_global[t_global][i_cd2_10];
290 
291  corr_cross_global[t_global] += corr_1124 * corr_3342;
292 
293  dcomplex corr_3144 = corr_cross_3144_global[t_global][i_cd2_01];
294  dcomplex corr_1322 = corr_cross_1322_global[t_global][i_cd2_10];
295 
296  corr_cross_global[t_global] += corr_3144 * corr_1322;
297  }
298  }
299  }
300  }
301  }
302 
303 
304  //- write outputs
305  for (int t_global = 0; t_global < Lt; ++t_global) {
306  corr_global[t_global] = corr_direct_global[t_global] - corr_cross_global[t_global];
307 
308  vout.general(m_vl, " %4d %20.12e %20.12e\n",
309  t_global, real(corr_global[t_global]), imag(corr_global[t_global]));
310  }
311  vout.general(m_vl, "\n");
312 
313 
314  vout.detailed(m_vl, " 4pt_correlator_direct:\n");
315  for (int t_global = 0; t_global < Lt; ++t_global) {
316  vout.detailed(m_vl, " %4d %20.12e %20.12e\n",
317  t_global, real(corr_direct_global[t_global]), imag(corr_direct_global[t_global]));
318  }
319  vout.detailed(m_vl, "\n");
320 
321  vout.detailed(m_vl, " 4pt_correlator_cross:\n");
322  for (int t_global = 0; t_global < Lt; ++t_global) {
323  vout.detailed(m_vl, " %4d %20.12e %20.12e\n",
324  t_global, real(corr_cross_global[t_global]), imag(corr_cross_global[t_global]));
325  }
326  vout.detailed(m_vl, "\n");
327 }
328 
329 
330 //====================================================================
331 double Corr4pt_4spinor::meson_momentum_all(const std::vector<Field_F>& sq1,
332  const std::vector<Field_F>& sq2,
333  const std::vector<Field_F>& sq3,
334  const std::vector<Field_F>& sq4,
335  const std::vector<int>& source_position)
336 {
337  const int Ndim = CommonParameters::Ndim();
338  const int Lt = CommonParameters::Lt();
339 
340  std::vector<dcomplex> corr_global(Lt);
341  GammaMatrix gm_sink_12, gm_sink_34, gm_src_12, gm_src_34;
342 
343  const int N_momentum = 10;
344 
345  typedef std::vector<int> MomentumSet;
346  std::vector<MomentumSet> momentum_sink(N_momentum);
347  for (int i_momentum = 0; i_momentum < N_momentum; i_momentum++) {
348  momentum_sink[i_momentum].resize(Ndim - 1);
349  }
350 
351  //- momentum_sink[0] = (1,0,0)
352  int i_momentum = 0;
353  momentum_sink[i_momentum][0] = 1;
354  momentum_sink[i_momentum][1] = 0;
355  momentum_sink[i_momentum][2] = 0;
356 
357  //- momentum_sink[1] = (0,1,0)
358  i_momentum = 1;
359  momentum_sink[i_momentum][0] = 0;
360  momentum_sink[i_momentum][1] = 1;
361  momentum_sink[i_momentum][2] = 0;
362 
363  //- momentum_sink[2] = (0,0,1)
364  i_momentum = 2;
365  momentum_sink[i_momentum][0] = 0;
366  momentum_sink[i_momentum][1] = 0;
367  momentum_sink[i_momentum][2] = 1;
368 
369  //- momentum_sink[3] = (1,1,0)
370  i_momentum = 3;
371  momentum_sink[i_momentum][0] = 1;
372  momentum_sink[i_momentum][1] = 1;
373  momentum_sink[i_momentum][2] = 0;
374 
375  //- momentum_sink[4] = (0,1,1)
376  i_momentum = 4;
377  momentum_sink[i_momentum][0] = 0;
378  momentum_sink[i_momentum][1] = 1;
379  momentum_sink[i_momentum][2] = 1;
380 
381  //- momentum_sink[5] = (1,0,1)
382  i_momentum = 5;
383  momentum_sink[i_momentum][0] = 1;
384  momentum_sink[i_momentum][1] = 0;
385  momentum_sink[i_momentum][2] = 1;
386 
387  //- momentum_sink[6] = (1,1,1)
388  i_momentum = 6;
389  momentum_sink[i_momentum][0] = 1;
390  momentum_sink[i_momentum][1] = 1;
391  momentum_sink[i_momentum][2] = 1;
392 
393  //- momentum_sink[7] = (2,0,0)
394  i_momentum = 7;
395  momentum_sink[i_momentum][0] = 2;
396  momentum_sink[i_momentum][1] = 0;
397  momentum_sink[i_momentum][2] = 0;
398 
399  //- momentum_sink[8] = (0,2,0)
400  i_momentum = 8;
401  momentum_sink[i_momentum][0] = 0;
402  momentum_sink[i_momentum][1] = 2;
403  momentum_sink[i_momentum][2] = 0;
404 
405  //- momentum_sink[9] = (0,0,2)
406  i_momentum = 9;
407  momentum_sink[i_momentum][0] = 0;
408  momentum_sink[i_momentum][1] = 0;
409  momentum_sink[i_momentum][2] = 2;
410 
411 
412  std::ostream& log_file_previous = vout.getStream();
413  std::ofstream log_file;
414 
415  if (m_filename_output != "stdout") {
416  log_file.open(m_filename_output.c_str(), std::ios::app);
417  vout.init(log_file);
418  }
419 
420 
421  gm_src_12 = m_gmset->get_GM(m_gmset->GAMMA5);
422  gm_src_34 = m_gmset->get_GM(m_gmset->GAMMA5);
423  gm_sink_12 = m_gmset->get_GM(m_gmset->GAMMA5);
424  gm_sink_34 = m_gmset->get_GM(m_gmset->GAMMA5);
425  for (int i_momentum = 0; i_momentum < N_momentum; i_momentum++) {
426  vout.general(m_vl, "PS_momentum(%d %d %d) <-- PS correlator:\n",
427  momentum_sink[i_momentum][0],
428  momentum_sink[i_momentum][1],
429  momentum_sink[i_momentum][2]);
430  meson_momentum_correlator(corr_global, momentum_sink[i_momentum],
431  gm_sink_12, gm_sink_34, gm_src_12, gm_src_34,
432  sq1, sq2, sq3, sq4, source_position);
433  }
434 
435  gm_src_12 = m_gmset->get_GM(m_gmset->GAMMA1);
436  gm_src_34 = m_gmset->get_GM(m_gmset->GAMMA1);
437  gm_sink_12 = m_gmset->get_GM(m_gmset->GAMMA1);
438  gm_sink_34 = m_gmset->get_GM(m_gmset->GAMMA1);
439  for (int i_momentum = 0; i_momentum < N_momentum; i_momentum++) {
440  vout.general(m_vl, "V1_momentum(%d %d %d) <-- V1 correlator:\n",
441  momentum_sink[i_momentum][0],
442  momentum_sink[i_momentum][1],
443  momentum_sink[i_momentum][2]);
444  meson_momentum_correlator(corr_global, momentum_sink[i_momentum],
445  gm_sink_12, gm_sink_34, gm_src_12, gm_src_34,
446  sq1, sq2, sq3, sq4, source_position);
447  }
448 
449  gm_src_12 = m_gmset->get_GM(m_gmset->GAMMA2);
450  gm_src_34 = m_gmset->get_GM(m_gmset->GAMMA2);
451  gm_sink_12 = m_gmset->get_GM(m_gmset->GAMMA2);
452  gm_sink_34 = m_gmset->get_GM(m_gmset->GAMMA2);
453  for (int i_momentum = 0; i_momentum < N_momentum; i_momentum++) {
454  vout.general(m_vl, "V2_momentum(%d %d %d) <-- V2 correlator:\n",
455  momentum_sink[i_momentum][0],
456  momentum_sink[i_momentum][1],
457  momentum_sink[i_momentum][2]);
458  meson_momentum_correlator(corr_global, momentum_sink[i_momentum],
459  gm_sink_12, gm_sink_34, gm_src_12, gm_src_34,
460  sq1, sq2, sq3, sq4, source_position);
461  }
462 
463  gm_src_12 = m_gmset->get_GM(m_gmset->GAMMA3);
464  gm_src_34 = m_gmset->get_GM(m_gmset->GAMMA3);
465  gm_sink_12 = m_gmset->get_GM(m_gmset->GAMMA3);
466  gm_sink_34 = m_gmset->get_GM(m_gmset->GAMMA3);
467  for (int i_momentum = 0; i_momentum < N_momentum; i_momentum++) {
468  vout.general(m_vl, "V3_momentum(%d %d %d) <-- V3 correlator:\n",
469  momentum_sink[i_momentum][0],
470  momentum_sink[i_momentum][1],
471  momentum_sink[i_momentum][2]);
472  meson_momentum_correlator(corr_global, momentum_sink[i_momentum],
473  gm_sink_12, gm_sink_34, gm_src_12, gm_src_34,
474  sq1, sq2, sq3, sq4, source_position);
475  }
476 
477 
478  if (m_filename_output != "stdout") {
479  log_file.close();
480  vout.init(log_file_previous);
481  }
482 
483  return EXIT_SUCCESS;
484 }
485 
486 
487 //====================================================================
488 void Corr4pt_4spinor::meson_momentum_correlator(std::vector<dcomplex>& corr_global,
489  const std::vector<int>& momentum_sink,
490  const GammaMatrix& gm_sink_12,
491  const GammaMatrix& gm_sink_34,
492  const GammaMatrix& gm_src_12,
493  const GammaMatrix& gm_src_34,
494  const std::vector<Field_F>& sq1,
495  const std::vector<Field_F>& sq2,
496  const std::vector<Field_F>& sq3,
497  const std::vector<Field_F>& sq4,
498  const std::vector<int>& source_position)
499 {
500  const int Nc = CommonParameters::Nc();
501  const int Nd = CommonParameters::Nd();
502  const int Lt = CommonParameters::Lt();
503  const int Nt = CommonParameters::Nt();
504 
505  assert(corr_global.size() == Lt);
506 
507 #if 0
508  GammaMatrix gm_gm5_src, gm5_gm_sink, gm5;
509  gm5 = m_gmset->get_GM(m_gmset->GAMMA5);
510  gm_gm5_src = gm_src.mult(gm5);
511  gm5_gm_sink = gm5.mult(gm_sink);
512 
513  std::vector<dcomplex> corr_local(Nt);
514 
515  for (int c0 = 0; c0 < Nc; ++c0) {
516  for (int d0 = 0; d0 < Nd; ++d0) {
517  int d1 = gm_gm5_src.index(d0);
518 
519  for (int t = 0; t < Nt; ++t) {
520  dcomplex corr_t;
521 
522  contract_at_t(corr_t, momentum_sink, gm5_gm_sink, source_position,
523  sq1[c0 + Nc * d0], sq2[c0 + Nc * d1], t);
524 
525  corr_local[t] += gm_gm5_src.value(d0) * corr_t;
526  }
527  }
528  }
529 
530  global_corr_t(corr_global, corr_local);
531 
532  for (int t = 0; t < corr_global.size(); ++t) {
533  vout.general(m_vl, " %4d %20.12e %20.12e\n",
534  t, real(corr_global[t]), imag(corr_global[t]));
535  }
536 #endif
537 }
538 
539 
540 //====================================================================
541 void Corr4pt_4spinor::corr_direct(std::vector<dcomplex>& corr_direct_global,
542  const GammaMatrix& gm5_gm_sink,
543  const GammaMatrix& gm_gm5_src,
544  const std::vector<Field_F>& sq1,
545  const std::vector<Field_F>& sq2)
546 {
547  const int Nc = CommonParameters::Nc();
548  const int Nd = CommonParameters::Nd();
549  const int Nt = CommonParameters::Nt();
550 
551  std::vector<dcomplex> corr_local(Nt, cmplx(0.0, 0.0));
552 
553  for (int c0 = 0; c0 < Nc; ++c0) {
554  for (int d0 = 0; d0 < Nd; ++d0) {
555  int d1 = gm_gm5_src.index(d0);
556 
557  for (int t = 0; t < Nt; ++t) {
558  dcomplex corr_t;
559 
560  contract_at_t(corr_t, gm5_gm_sink,
561  sq1[c0 + Nc * d0], sq2[c0 + Nc * d1], t);
562 
563  corr_local[t] += gm_gm5_src.value(d0) * corr_t;
564  }
565  }
566  }
567  global_corr_t(corr_direct_global, corr_local);
568 }
569 
570 
571 //====================================================================
572 void Corr4pt_4spinor::corr_cross_sink(std::vector<std::vector<dcomplex> >& corr_cross_global,
573  const GammaMatrix& gm5_gm_sink,
574  const GammaMatrix& gm_gm5_src,
575  const std::vector<Field_F>& sq1,
576  const std::vector<Field_F>& sq2)
577 {
578  const int Nc = CommonParameters::Nc();
579  const int Nd = CommonParameters::Nd();
580  const int Nt = CommonParameters::Nt();
581  const int Lt = CommonParameters::Lt();
582 
583  typedef std::vector<dcomplex> CorrSet;
584  std::vector<CorrSet> corr_local_idx(Nt);
585  for (int t = 0; t < Nt; ++t) {
586  corr_local_idx[t].resize(Nc * Nd * Nc * Nd);
587  }
588  for (int t = 0; t < Nt; ++t) {
589  for (int i_cd2 = 0; i_cd2 < Nc * Nd * Nc * Nd; ++i_cd2) {
590  corr_local_idx[t][i_cd2] = cmplx(0.0, 0.0);
591  }
592  }
593 
594  for (int c0 = 0; c0 < Nc; ++c0) {
595  for (int d0 = 0; d0 < Nd; ++d0) {
596  for (int c1 = 0; c1 < Nc; ++c1) {
597  for (int d1 = 0; d1 < Nd; ++d1) {
598  int d51 = gm_gm5_src.index(d1);
599 
600  for (int t = 0; t < Nt; ++t) {
601  dcomplex corr_t;
602 
603  contract_at_t(corr_t, gm5_gm_sink,
604  sq1[c0 + Nc * d0], sq2[c1 + Nc * d51], t);
605 
606  int i_cd2 = c0 + Nc * d0 + Nc * Nd * c1 + Nc * Nd * Nc * d51;
607 
608  corr_local_idx[t][i_cd2] = gm_gm5_src.value(d1) * corr_t;
609  }
610  }
611  }
612  }
613  }
614 
615 
616  for (int i_cd2 = 0; i_cd2 < Nc * Nd * Nc * Nd; ++i_cd2) {
617  std::vector<dcomplex> corr_local(Nt);
618  std::vector<dcomplex> corr_global(Lt);
619 
620  for (int t = 0; t < Nt; ++t) {
621  corr_local[t] = corr_local_idx[t][i_cd2];
622  }
623 
624  global_corr_t(corr_global, corr_local);
625 
626  for (int t_global = 0; t_global < Lt; ++t_global) {
627  corr_cross_global[t_global][i_cd2] = corr_global[t_global];
628  }
629  }
630 }
631 
632 
633 //====================================================================
634 void Corr4pt_4spinor::global_corr_t(std::vector<dcomplex>& corr_global,
635  const std::vector<dcomplex>& corr_local)
636 {
637  const int Lt = CommonParameters::Lt();
638  const int Nt = CommonParameters::Nt();
639 
640  assert(corr_global.size() == Lt);
641  assert(corr_local.size() == Nt);
642 
643  const int ipe_t = Communicator::ipe(3);
644 
645  std::vector<dcomplex> corr_tmp(Lt, cmplx(0.0, 0.0));
646 
647  for (int t = 0; t < Nt; ++t) {
648  int t_global = t + ipe_t * Nt;
649  corr_tmp[t_global] = corr_local[t];
650  }
651 
652  for (int t_global = 0; t_global < Lt; ++t_global) {
653  double cr_r = Communicator::reduce_sum(real(corr_tmp[t_global]));
654  double cr_i = Communicator::reduce_sum(imag(corr_tmp[t_global]));
655 
656  corr_global[t_global] = cmplx(cr_r, cr_i);
657  }
658 }
659 
660 
661 //====================================================================
662 //============================================================END=====
GammaMatrix::mult
GammaMatrix mult(GammaMatrix) const
Definition: gammaMatrix.cpp:39
GammaMatrixSet::GAMMA51
@ GAMMA51
Definition: gammaMatrixSet.h:49
CommonParameters::Ny
static int Ny()
Definition: commonParameters.h:106
CommonParameters::Nz
static int Nz()
Definition: commonParameters.h:107
Corr4pt_4spinor::global_corr_t
void global_corr_t(std::vector< dcomplex > &corr_global, const std::vector< dcomplex > &corr_local)
transform node-local correlator in t to global.
Definition: corr4pt_4spinor.cpp:634
GammaMatrixSet::GAMMA5
@ GAMMA5
Definition: gammaMatrixSet.h:48
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
GammaMatrixSet::GAMMA1
@ GAMMA1
Definition: gammaMatrixSet.h:48
GammaMatrix::index
int index(int row) const
Definition: gammaMatrix.h:83
CommonParameters::Ndim
static int Ndim()
Definition: commonParameters.h:117
Parameters
Class for parameters.
Definition: parameters.h:46
GammaMatrix
Gamma Matrix class.
Definition: gammaMatrix.h:44
GammaMatrixSet::UNITY
@ UNITY
Definition: gammaMatrixSet.h:48
Bridge::BridgeIO::detailed
void detailed(const char *format,...)
Definition: bridgeIO.cpp:219
CommonParameters::Nvol
static int Nvol()
Definition: commonParameters.h:109
Corr4pt_4spinor::corr_cross_sink
void corr_cross_sink(std::vector< std::vector< dcomplex > > &corr_cross_global, const GammaMatrix &gm5_gm_sink, const GammaMatrix &gm_gm5_src, const std::vector< Field_F > &sq1, const std::vector< Field_F > &sq2)
Definition: corr4pt_4spinor.cpp:572
Corr4pt_4spinor::m_gmset
GammaMatrixSet * m_gmset
Definition: corr4pt_4spinor.h:43
GammaMatrixSet::GAMMA3
@ GAMMA3
Definition: gammaMatrixSet.h:48
GammaMatrixSet::GAMMA53
@ GAMMA53
Definition: gammaMatrixSet.h:49
GammaMatrixSet::GAMMA54
@ GAMMA54
Definition: gammaMatrixSet.h:49
Corr4pt_4spinor::corr_direct
void corr_direct(std::vector< dcomplex > &corr_direct_global, const GammaMatrix &gm5_gm_sink, const GammaMatrix &gm_gm5_src, const std::vector< Field_F > &sq1, const std::vector< Field_F > &sq2)
totally antisymmetric tensor: index.
Definition: corr4pt_4spinor.cpp:541
CommonParameters::Nx
static int Nx()
Definition: commonParameters.h:105
Corr4pt_4spinor::set_parameters
void set_parameters(const Parameters &params)
Definition: corr4pt_4spinor.cpp:19
Corr4pt_4spinor::get_parameters
void get_parameters(Parameters &params) const
Definition: corr4pt_4spinor.cpp:34
Corr4pt_4spinor::meson_momentum_all
double meson_momentum_all(const std::vector< Field_F > &sq1, const std::vector< Field_F > &sq2, const std::vector< Field_F > &sq3, const std::vector< Field_F > &sq4, const std::vector< int > &source_position)
Definition: corr4pt_4spinor.cpp:331
CommonParameters::Lt
static int Lt()
Definition: commonParameters.h:94
Communicator::reduce_sum
static int reduce_sum(int count, dcomplex *recv_buf, dcomplex *send_buf, int pattern=0)
make a global sum of an array of dcomplex over the communicator. pattern specifies the dimensions to ...
Definition: communicator.cpp:263
CommonParameters::Nc
static int Nc()
Definition: commonParameters.h:115
GammaMatrix::value
dcomplex value(int row) const
Definition: gammaMatrix.h:88
contract_at_t
void contract_at_t(dcomplex &corr, const GammaMatrix &gm_sink, const Field_F &v1, const Field_F &v2, const int time)
Contraction of hadron for 4-spinor fermion.
Definition: contract_4spinor.cpp:37
CommonParameters::Nt
static int Nt()
Definition: commonParameters.h:108
Corr4pt_4spinor::meson_momentum_correlator
void meson_momentum_correlator(std::vector< dcomplex > &corr_global, const std::vector< int > &momentum_sink, const GammaMatrix &gm_sink_12, const GammaMatrix &gm_sink_34, const GammaMatrix &gm_src_21, const GammaMatrix &gm_src_43, const std::vector< Field_F > &sq1, const std::vector< Field_F > &sq2, const std::vector< Field_F > &sq3, const std::vector< Field_F > &sq4, const std::vector< int > &source_position)
Definition: corr4pt_4spinor.cpp:488
Corr4pt_4spinor::meson_all
double meson_all(const std::vector< Field_F > &sq1, const std::vector< Field_F > &sq2, const std::vector< Field_F > &sq3, const std::vector< Field_F > &sq4)
Definition: corr4pt_4spinor.cpp:42
GammaMatrixSet::GAMMA45
@ GAMMA45
Definition: gammaMatrixSet.h:50
Corr4pt_4spinor::class_name
static const std::string class_name
Definition: corr4pt_4spinor.h:35
CommonParameters::Nd
static int Nd()
Definition: commonParameters.h:116
GammaMatrixSet::GAMMA52
@ GAMMA52
Definition: gammaMatrixSet.h:49
Bridge::BridgeIO::set_verbose_level
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:133
Corr4pt_4spinor::m_filename_output
std::string m_filename_output
Definition: corr4pt_4spinor.h:41
Corr4pt_4spinor::meson_correlator
void meson_correlator(std::vector< dcomplex > &corr_global, const GammaMatrix &gm_sink_12, const GammaMatrix &gm_sink_34, const GammaMatrix &gm_src_21, const GammaMatrix &gm_src_43, const std::vector< Field_F > &sq1, const std::vector< Field_F > &sq2, const std::vector< Field_F > &sq3, const std::vector< Field_F > &sq4)
Definition: corr4pt_4spinor.cpp:165
GammaMatrixSet::SIGMA23
@ SIGMA23
Definition: gammaMatrixSet.h:51
Communicator::ipe
static int ipe(const int dir)
logical coordinate of current proc.
Definition: communicator.cpp:105
Parameters::fetch_string
int fetch_string(const string &key, string &value) const
Definition: parameters.cpp:378
GammaMatrixSet::SIGMA12
@ SIGMA12
Definition: gammaMatrixSet.h:51
GammaMatrixSet::get_GM
GammaMatrix get_GM(GMspecies spec)
Definition: gammaMatrixSet.h:76
Parameters::get_string
string get_string(const string &key) const
Definition: parameters.cpp:221
GammaMatrixSet::GAMMA2
@ GAMMA2
Definition: gammaMatrixSet.h:48
GammaMatrixSet::SIGMA31
@ SIGMA31
Definition: gammaMatrixSet.h:51
Bridge::BridgeIO::getStream
std::ostream & getStream()
Definition: bridgeIO.cpp:396
Bridge::BridgeIO::general
void general(const char *format,...)
Definition: bridgeIO.cpp:200
corr4pt_4spinor.h
Bridge::vout
BridgeIO vout
Definition: bridgeIO.cpp:512
Corr4pt_4spinor::m_vl
Bridge::VerboseLevel m_vl
Definition: corr4pt_4spinor.h:38
Bridge::BridgeIO::get_verbose_level
static std::string get_verbose_level(const VerboseLevel vl)
Definition: bridgeIO.cpp:154