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