Bridge++  Ver. 2.0.4
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  const bool use_outputfile = (m_filename_output != "stdout");
53  if (use_outputfile) {
54  int rank_io = 0;
55  vout.init(m_filename_output, rank_io, std::ios::app);
56  }
57 
58  gm_src_12 = m_gmset->get_GM(m_gmset->GAMMA5);
59  gm_src_34 = m_gmset->get_GM(m_gmset->GAMMA5);
60  gm_sink_12 = m_gmset->get_GM(m_gmset->GAMMA5);
61  gm_sink_34 = m_gmset->get_GM(m_gmset->GAMMA5);
62  vout.general(m_vl, "PS <-- PS 4pt_correlator:\n");
63  meson_correlator(corr_global, gm_sink_12, gm_sink_34, gm_src_12, gm_src_34, sq1, sq2, sq3, sq4);
64  double result = real(corr_global[0]);
65 
66 #if 0
67  gm_src_12 = m_gmset->get_GM(m_gmset->GAMMA1);
68  gm_src_34 = m_gmset->get_GM(m_gmset->GAMMA1);
69  gm_sink_12 = m_gmset->get_GM(m_gmset->GAMMA1);
70  gm_sink_34 = m_gmset->get_GM(m_gmset->GAMMA1);
71  vout.general(m_vl, "V1 <-- V1 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->GAMMA2);
75  gm_src_34 = m_gmset->get_GM(m_gmset->GAMMA2);
76  gm_sink_12 = m_gmset->get_GM(m_gmset->GAMMA2);
77  gm_sink_34 = m_gmset->get_GM(m_gmset->GAMMA2);
78  vout.general(m_vl, "V2 <-- V2 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->GAMMA3);
82  gm_src_34 = m_gmset->get_GM(m_gmset->GAMMA3);
83  gm_sink_12 = m_gmset->get_GM(m_gmset->GAMMA3);
84  gm_sink_34 = m_gmset->get_GM(m_gmset->GAMMA3);
85  vout.general(m_vl, "V3 <-- V3 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->GAMMA5);
89  gm_src_34 = m_gmset->get_GM(m_gmset->GAMMA5);
90  gm_sink_12 = m_gmset->get_GM(m_gmset->GAMMA45);
91  gm_sink_34 = m_gmset->get_GM(m_gmset->GAMMA45);
92  vout.general(m_vl, "A4 <-- PS 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->GAMMA54);
96  gm_src_34 = m_gmset->get_GM(m_gmset->GAMMA54);
97  gm_sink_12 = m_gmset->get_GM(m_gmset->GAMMA5);
98  gm_sink_34 = m_gmset->get_GM(m_gmset->GAMMA5);
99  vout.general(m_vl, "PS <-- A4 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->UNITY);
103  gm_src_34 = m_gmset->get_GM(m_gmset->UNITY);
104  gm_sink_12 = m_gmset->get_GM(m_gmset->UNITY);
105  gm_sink_34 = m_gmset->get_GM(m_gmset->UNITY);
106  vout.general(m_vl, "S <-- S 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->GAMMA51);
110  gm_src_34 = m_gmset->get_GM(m_gmset->GAMMA51);
111  gm_sink_12 = m_gmset->get_GM(m_gmset->GAMMA51);
112  gm_sink_34 = m_gmset->get_GM(m_gmset->GAMMA51);
113  vout.general(m_vl, "GAMMA51 <-- GAMMA51 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->GAMMA52);
117  gm_src_34 = m_gmset->get_GM(m_gmset->GAMMA52);
118  gm_sink_12 = m_gmset->get_GM(m_gmset->GAMMA52);
119  gm_sink_34 = m_gmset->get_GM(m_gmset->GAMMA52);
120  vout.general(m_vl, "GAMMA52 <-- GAMMA52 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->GAMMA53);
124  gm_src_34 = m_gmset->get_GM(m_gmset->GAMMA53);
125  gm_sink_12 = m_gmset->get_GM(m_gmset->GAMMA53);
126  gm_sink_34 = m_gmset->get_GM(m_gmset->GAMMA53);
127  vout.general(m_vl, "GAMMA53 <-- GAMMA53 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->SIGMA12);
131  gm_src_34 = m_gmset->get_GM(m_gmset->SIGMA12);
132  gm_sink_12 = m_gmset->get_GM(m_gmset->SIGMA12);
133  gm_sink_34 = m_gmset->get_GM(m_gmset->SIGMA12);
134  vout.general(m_vl, "SIGMA12 <-- SIGMA12 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->SIGMA23);
138  gm_src_34 = m_gmset->get_GM(m_gmset->SIGMA23);
139  gm_sink_12 = m_gmset->get_GM(m_gmset->SIGMA23);
140  gm_sink_34 = m_gmset->get_GM(m_gmset->SIGMA23);
141  vout.general(m_vl, "SIGMA23 <-- SIGMA23 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 
144  gm_src_12 = m_gmset->get_GM(m_gmset->SIGMA31);
145  gm_src_34 = m_gmset->get_GM(m_gmset->SIGMA31);
146  gm_sink_12 = m_gmset->get_GM(m_gmset->SIGMA31);
147  gm_sink_34 = m_gmset->get_GM(m_gmset->SIGMA31);
148  vout.general(m_vl, "SIGMA31 <-- SIGMA31 4pt_correlator:\n");
149  meson_correlator(corr_global, gm_sink_12, gm_sink_34, gm_src_12, gm_src_34, sq1, sq2, sq3, sq4);
150 #endif
151 
152  if (use_outputfile) {
153  vout.unset();
154  }
155 
156  return result;
157 }
158 
159 
160 //====================================================================
161 void Corr4pt_4spinor::meson_correlator(std::vector<dcomplex>& corr_global,
162  const GammaMatrix& gm_sink_12,
163  const GammaMatrix& gm_sink_34,
164  const GammaMatrix& gm_src_12,
165  const GammaMatrix& gm_src_34,
166  const std::vector<Field_F>& sq1,
167  const std::vector<Field_F>& sq2,
168  const std::vector<Field_F>& sq3,
169  const std::vector<Field_F>& sq4)
170 {
171  const int Nc = CommonParameters::Nc();
172  const int Nd = CommonParameters::Nd();
173  const int Nx = CommonParameters::Nx();
174  const int Ny = CommonParameters::Ny();
175  const int Nz = CommonParameters::Nz();
176  const int Nt = CommonParameters::Nt();
177  const int Nvol = CommonParameters::Nvol();
178 
179  const int Lt = CommonParameters::Lt();
180 
181  assert(corr_global.size() == Lt);
182 
183  //- M_{ij}(x,t) := \bar{q_i}^c \Gamma_ij q_j^c
184  // Prop_ij(t) := q_i(t_sink) \bar{q_j}(t_src)
185  //
186  //- corr_direct = tr(Prop_11(t) Prop_22^{\dagger}(t)) tr(Prop_33(t) Prop_44^{\dagger}(t))
187  // +tr(Prop_13(t) Prop_24^{\dagger}(t)) tr(Prop_31(t) Prop_42^{\dagger}(t))
188  // x x' x x'
189  // t_sink 1 2 3 4 1 2 3 4 1 2
190  // | | | | / / / / / /
191  // ^ v ^ v + v ^ v ^ v
192  // | | | | / / / /
193  // t_src 1 2 3 4 1 2 3 4
194  // c0 c1 c0 c1
195  //
196  //- corr_cross = Prop_11(t)_{c0,c2} Prop_24^{\dagger}(t)_{c2,c1} Prop_33(t)_{c1,c3} Prop_42^{\dagger}(t)_{c3,c0}
197  // +Prop_31(t)_{c0,c3} Prop_44^{\dagger}(t)_{c3,c1} Prop_13(t)_{c1,c2} Prop_22^{\dagger}(t)_{c2,c0}
198  // = Prop_1124(t)_{c0,c1} Prop_3342(t)_{c1,c0}
199  // +Prop_3144(t)_{c0,c1} Prop_1322(t)_{c1,c0}
200  // x x' x x'
201  // t_sink 1 2 3 4 1 2 3 4
202  // | \|/ \|/ |
203  // ^ * + * v
204  // | /|\ /|\ |
205  // t_src 1 2 3 4 1 2 3 4
206  // c0 c1 c0 c1
207 
208  const GammaMatrix gm5 = m_gmset->get_GM(m_gmset->GAMMA5);
209 
210  const GammaMatrix gm_gm5_src_12 = gm_src_12.mult(gm5);
211  const GammaMatrix gm_gm5_src_34 = gm_src_34.mult(gm5);
212 
213  const GammaMatrix gm5_gm_sink_12 = gm5.mult(gm_sink_12);
214  const GammaMatrix gm5_gm_sink_34 = gm5.mult(gm_sink_34);
215 
216  //- corr_direct = tr(Prop_11(t) Prop_22^{\dagger}(t)) tr(Prop_33(t) Prop_44^{\dagger}(t))
217  // +tr(Prop_13(t) Prop_24^{\dagger}(t)) tr(Prop_31(t) Prop_42^{\dagger}(t))
218 
219  std::vector<dcomplex> corr_direct_1122_global(Lt);
220  std::vector<dcomplex> corr_direct_3344_global(Lt);
221  std::vector<dcomplex> corr_direct_1324_global(Lt);
222  std::vector<dcomplex> corr_direct_3142_global(Lt);
223 
224  corr_direct(corr_direct_1122_global,
225  gm5_gm_sink_12, gm_gm5_src_12,
226  sq1, sq2);
227  corr_direct(corr_direct_3344_global,
228  gm5_gm_sink_34, gm_gm5_src_34,
229  sq3, sq4);
230  corr_direct(corr_direct_1324_global,
231  gm5_gm_sink_12, gm_gm5_src_34,
232  sq3, sq4);
233  corr_direct(corr_direct_3142_global,
234  gm5_gm_sink_34, gm_gm5_src_12,
235  sq1, sq2);
236 
237  std::vector<dcomplex> corr_direct_global(Lt);
238  for (int t_global = 0; t_global < Lt; ++t_global) {
239  corr_direct_global[t_global]
240  = corr_direct_1122_global[t_global] * corr_direct_3344_global[t_global]
241  + corr_direct_1324_global[t_global] * corr_direct_3142_global[t_global];
242  }
243 
244  //- corr_cross = Prop_11(t)_{c0,c2} Prop_24^{\dagger}(t)_{c2,c1} Prop_33(t)_{c1,c3} Prop_42^{\dagger}(t)_{c3,c0}
245  // +Prop_31(t)_{c0,c3} Prop_44^{\dagger}(t)_{c3,c1} Prop_13(t)_{c1,c2} Prop_22^{\dagger}(t)_{c2,c0}
246  // = Prop_1124(t)_{c0,c1} Prop_3342(t)_{c1,c0}
247  // +Prop_3144(t)_{c0,c1} Prop_1322(t)_{c1,c0}
248 
249  typedef std::vector<dcomplex> CorrSet;
250  std::vector<CorrSet> corr_cross_1124_global(Lt);
251  std::vector<CorrSet> corr_cross_3342_global(Lt);
252  std::vector<CorrSet> corr_cross_1322_global(Lt);
253  std::vector<CorrSet> corr_cross_3144_global(Lt);
254  for (int t_global = 0; t_global < Lt; ++t_global) {
255  corr_cross_1124_global[t_global].resize(Nc * Nd * Nc * Nd);
256  corr_cross_3342_global[t_global].resize(Nc * Nd * Nc * Nd);
257  corr_cross_1322_global[t_global].resize(Nc * Nd * Nc * Nd);
258  corr_cross_3144_global[t_global].resize(Nc * Nd * Nc * Nd);
259  }
260 
261  corr_cross_sink(corr_cross_1124_global,
262  gm5_gm_sink_12, gm_gm5_src_12,
263  sq1, sq4);
264  corr_cross_sink(corr_cross_3342_global,
265  gm5_gm_sink_34, gm_gm5_src_34,
266  sq3, sq2);
267  corr_cross_sink(corr_cross_3144_global,
268  gm5_gm_sink_34, gm_gm5_src_12,
269  sq1, sq4);
270  corr_cross_sink(corr_cross_1322_global,
271  gm5_gm_sink_12, gm_gm5_src_34,
272  sq3, sq2);
273 
274  //- sum up at src
275  std::vector<dcomplex> corr_cross_global(Lt, cmplx(0.0, 0.0));
276  for (int t_global = 0; t_global < Lt; ++t_global) {
277  for (int c0 = 0; c0 < Nc; ++c0) {
278  for (int d0 = 0; d0 < Nd; ++d0) {
279  for (int c1 = 0; c1 < Nc; ++c1) {
280  for (int d1 = 0; d1 < Nd; ++d1) {
281  int i_cd2_01 = c0 + Nc * d0 + Nc * Nd * c1 + Nc * Nd * Nc * d1;
282  int i_cd2_10 = c1 + Nc * d1 + Nc * Nd * c0 + Nc * Nd * Nc * d0;
283 
284  dcomplex corr_1124 = corr_cross_1124_global[t_global][i_cd2_01];
285  dcomplex corr_3342 = corr_cross_3342_global[t_global][i_cd2_10];
286 
287  corr_cross_global[t_global] += corr_1124 * corr_3342;
288 
289  dcomplex corr_3144 = corr_cross_3144_global[t_global][i_cd2_01];
290  dcomplex corr_1322 = corr_cross_1322_global[t_global][i_cd2_10];
291 
292  corr_cross_global[t_global] += corr_3144 * corr_1322;
293  }
294  }
295  }
296  }
297  }
298 
299 
300  //- write outputs
301  for (int t_global = 0; t_global < Lt; ++t_global) {
302  corr_global[t_global] = corr_direct_global[t_global] - corr_cross_global[t_global];
303 
304  vout.general(m_vl, " %4d %20.12e %20.12e\n",
305  t_global, real(corr_global[t_global]), imag(corr_global[t_global]));
306  }
307  vout.general(m_vl, "\n");
308 
309 
310  vout.detailed(m_vl, " 4pt_correlator_direct:\n");
311  for (int t_global = 0; t_global < Lt; ++t_global) {
312  vout.detailed(m_vl, " %4d %20.12e %20.12e\n",
313  t_global, real(corr_direct_global[t_global]), imag(corr_direct_global[t_global]));
314  }
315  vout.detailed(m_vl, "\n");
316 
317  vout.detailed(m_vl, " 4pt_correlator_cross:\n");
318  for (int t_global = 0; t_global < Lt; ++t_global) {
319  vout.detailed(m_vl, " %4d %20.12e %20.12e\n",
320  t_global, real(corr_cross_global[t_global]), imag(corr_cross_global[t_global]));
321  }
322  vout.detailed(m_vl, "\n");
323 }
324 
325 
326 //====================================================================
327 double Corr4pt_4spinor::meson_momentum_all(const std::vector<Field_F>& sq1,
328  const std::vector<Field_F>& sq2,
329  const std::vector<Field_F>& sq3,
330  const std::vector<Field_F>& sq4,
331  const std::vector<int>& source_position)
332 {
333  const int Ndim = CommonParameters::Ndim();
334  const int Lt = CommonParameters::Lt();
335 
336  std::vector<dcomplex> corr_global(Lt);
337  GammaMatrix gm_sink_12, gm_sink_34, gm_src_12, gm_src_34;
338 
339  const int N_momentum = 10;
340 
341  typedef std::vector<int> MomentumSet;
342  std::vector<MomentumSet> momentum_sink(N_momentum);
343  for (int i_momentum = 0; i_momentum < N_momentum; i_momentum++) {
344  momentum_sink[i_momentum].resize(Ndim - 1);
345  }
346 
347  //- momentum_sink[0] = (1,0,0)
348  int i_momentum = 0;
349  momentum_sink[i_momentum][0] = 1;
350  momentum_sink[i_momentum][1] = 0;
351  momentum_sink[i_momentum][2] = 0;
352 
353  //- momentum_sink[1] = (0,1,0)
354  i_momentum = 1;
355  momentum_sink[i_momentum][0] = 0;
356  momentum_sink[i_momentum][1] = 1;
357  momentum_sink[i_momentum][2] = 0;
358 
359  //- momentum_sink[2] = (0,0,1)
360  i_momentum = 2;
361  momentum_sink[i_momentum][0] = 0;
362  momentum_sink[i_momentum][1] = 0;
363  momentum_sink[i_momentum][2] = 1;
364 
365  //- momentum_sink[3] = (1,1,0)
366  i_momentum = 3;
367  momentum_sink[i_momentum][0] = 1;
368  momentum_sink[i_momentum][1] = 1;
369  momentum_sink[i_momentum][2] = 0;
370 
371  //- momentum_sink[4] = (0,1,1)
372  i_momentum = 4;
373  momentum_sink[i_momentum][0] = 0;
374  momentum_sink[i_momentum][1] = 1;
375  momentum_sink[i_momentum][2] = 1;
376 
377  //- momentum_sink[5] = (1,0,1)
378  i_momentum = 5;
379  momentum_sink[i_momentum][0] = 1;
380  momentum_sink[i_momentum][1] = 0;
381  momentum_sink[i_momentum][2] = 1;
382 
383  //- momentum_sink[6] = (1,1,1)
384  i_momentum = 6;
385  momentum_sink[i_momentum][0] = 1;
386  momentum_sink[i_momentum][1] = 1;
387  momentum_sink[i_momentum][2] = 1;
388 
389  //- momentum_sink[7] = (2,0,0)
390  i_momentum = 7;
391  momentum_sink[i_momentum][0] = 2;
392  momentum_sink[i_momentum][1] = 0;
393  momentum_sink[i_momentum][2] = 0;
394 
395  //- momentum_sink[8] = (0,2,0)
396  i_momentum = 8;
397  momentum_sink[i_momentum][0] = 0;
398  momentum_sink[i_momentum][1] = 2;
399  momentum_sink[i_momentum][2] = 0;
400 
401  //- momentum_sink[9] = (0,0,2)
402  i_momentum = 9;
403  momentum_sink[i_momentum][0] = 0;
404  momentum_sink[i_momentum][1] = 0;
405  momentum_sink[i_momentum][2] = 2;
406 
407  const bool use_outputfile = (m_filename_output != "stdout");
408  if (use_outputfile) {
409  int rank_io = 0;
410  vout.init(m_filename_output, rank_io, std::ios::app);
411  }
412 
413  gm_src_12 = m_gmset->get_GM(m_gmset->GAMMA5);
414  gm_src_34 = m_gmset->get_GM(m_gmset->GAMMA5);
415  gm_sink_12 = m_gmset->get_GM(m_gmset->GAMMA5);
416  gm_sink_34 = m_gmset->get_GM(m_gmset->GAMMA5);
417  for (int i_momentum = 0; i_momentum < N_momentum; i_momentum++) {
418  vout.general(m_vl, "PS_momentum(%d %d %d) <-- PS correlator:\n",
419  momentum_sink[i_momentum][0],
420  momentum_sink[i_momentum][1],
421  momentum_sink[i_momentum][2]);
422  meson_momentum_correlator(corr_global, momentum_sink[i_momentum],
423  gm_sink_12, gm_sink_34, gm_src_12, gm_src_34,
424  sq1, sq2, sq3, sq4, source_position);
425  }
426 
427  gm_src_12 = m_gmset->get_GM(m_gmset->GAMMA1);
428  gm_src_34 = m_gmset->get_GM(m_gmset->GAMMA1);
429  gm_sink_12 = m_gmset->get_GM(m_gmset->GAMMA1);
430  gm_sink_34 = m_gmset->get_GM(m_gmset->GAMMA1);
431  for (int i_momentum = 0; i_momentum < N_momentum; i_momentum++) {
432  vout.general(m_vl, "V1_momentum(%d %d %d) <-- V1 correlator:\n",
433  momentum_sink[i_momentum][0],
434  momentum_sink[i_momentum][1],
435  momentum_sink[i_momentum][2]);
436  meson_momentum_correlator(corr_global, momentum_sink[i_momentum],
437  gm_sink_12, gm_sink_34, gm_src_12, gm_src_34,
438  sq1, sq2, sq3, sq4, source_position);
439  }
440 
441  gm_src_12 = m_gmset->get_GM(m_gmset->GAMMA2);
442  gm_src_34 = m_gmset->get_GM(m_gmset->GAMMA2);
443  gm_sink_12 = m_gmset->get_GM(m_gmset->GAMMA2);
444  gm_sink_34 = m_gmset->get_GM(m_gmset->GAMMA2);
445  for (int i_momentum = 0; i_momentum < N_momentum; i_momentum++) {
446  vout.general(m_vl, "V2_momentum(%d %d %d) <-- V2 correlator:\n",
447  momentum_sink[i_momentum][0],
448  momentum_sink[i_momentum][1],
449  momentum_sink[i_momentum][2]);
450  meson_momentum_correlator(corr_global, momentum_sink[i_momentum],
451  gm_sink_12, gm_sink_34, gm_src_12, gm_src_34,
452  sq1, sq2, sq3, sq4, source_position);
453  }
454 
455  gm_src_12 = m_gmset->get_GM(m_gmset->GAMMA3);
456  gm_src_34 = m_gmset->get_GM(m_gmset->GAMMA3);
457  gm_sink_12 = m_gmset->get_GM(m_gmset->GAMMA3);
458  gm_sink_34 = m_gmset->get_GM(m_gmset->GAMMA3);
459  for (int i_momentum = 0; i_momentum < N_momentum; i_momentum++) {
460  vout.general(m_vl, "V3_momentum(%d %d %d) <-- V3 correlator:\n",
461  momentum_sink[i_momentum][0],
462  momentum_sink[i_momentum][1],
463  momentum_sink[i_momentum][2]);
464  meson_momentum_correlator(corr_global, momentum_sink[i_momentum],
465  gm_sink_12, gm_sink_34, gm_src_12, gm_src_34,
466  sq1, sq2, sq3, sq4, source_position);
467  }
468 
469  if (use_outputfile) {
470  vout.unset();
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=====
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:624
GammaMatrixSet::GAMMA5
@ GAMMA5
Definition: gammaMatrixSet.h:48
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
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:281
Bridge::BridgeIO::unset
void unset()
Definition: bridgeIO.cpp:142
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:562
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:531
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:327
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:478
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:195
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:161
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::general
void general(const char *format,...)
Definition: bridgeIO.cpp:262
corr4pt_4spinor.h
Bridge::vout
BridgeIO vout
Definition: bridgeIO.cpp:569
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:216