Bridge++  Version 1.4.4
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
contract_4spinor.cpp
Go to the documentation of this file.
1 
14 #include "Measurements/Fermion/contract_4spinor.h"
15 #include "Field/index_lex.h"
16 
17 #include <cassert>
18 
19 #if defined USE_GROUP_SU3
20 #define C1 0
21 #define C2 1
22 #define C3 2
23 #endif
24 
25 //====================================================================
26 void contract_at_t(dcomplex& corr,
27  const GammaMatrix& gm_sink,
28  const Field_F& v1, const Field_F& v2,
29  const int time)
30 {
31  int Nc = CommonParameters::Nc();
32  int Nd = CommonParameters::Nd();
33  int Nx = CommonParameters::Nx();
34  int Ny = CommonParameters::Ny();
35  int Nz = CommonParameters::Nz();
36  int Nt = CommonParameters::Nt();
37  int Nvol = CommonParameters::Nvol();
38 
39  assert(Nvol == v1.nvol());
40  assert(Nvol == v2.nvol());
41  assert(time < Nt);
42 
43  Index_lex index;
44  std::vector<int> gm_index(Nd);
45  std::vector<double> corr_r(Nd), corr_i(Nd);
46 
47  for (int i = 0; i < Nd; ++i) {
48  gm_index[i] = gm_sink.index(i);
49  corr_r[i] = 0.0;
50  corr_i[i] = 0.0;
51  }
52 
53  for (int z = 0; z < Nz; ++z) {
54  for (int y = 0; y < Ny; ++y) {
55  for (int x = 0; x < Nx; ++x) {
56  int site = index.site(x, y, z, time);
57 
58  for (int s0 = 0; s0 < Nd; ++s0) {
59  int s1 = gm_index[s0];
60 
61  for (int c1 = 0; c1 < Nc; ++c1) {
62  corr_r[s0] += v1.cmp_r(c1, s1, site) * v2.cmp_r(c1, s0, site)
63  + v1.cmp_i(c1, s1, site) * v2.cmp_i(c1, s0, site);
64 
65  corr_i[s0] += -v1.cmp_r(c1, s1, site) * v2.cmp_i(c1, s0, site)
66  + v1.cmp_i(c1, s1, site) * v2.cmp_r(c1, s0, site);
67  }
68  }
69  }
70  }
71  }
72 
73  corr = cmplx(0.0, 0.0);
74  for (int s0 = 0; s0 < Nd; ++s0) {
75  corr += gm_sink.value(s0) * cmplx(corr_r[s0], corr_i[s0]);
76  }
77 }
78 
79 
80 //====================================================================
81 void contract_at_t(dcomplex& corr,
82  const std::vector<int>& momentum_sink,
83  const GammaMatrix& gm_sink,
84  const std::vector<int>& source_position,
85  const Field_F& v1, const Field_F& v2,
86  const int time)
87 {
88  const int Nc = CommonParameters::Nc();
89  const int Nd = CommonParameters::Nd();
90  const int Nx = CommonParameters::Nx();
91  const int Ny = CommonParameters::Ny();
92  const int Nz = CommonParameters::Nz();
93  const int Nt = CommonParameters::Nt();
94  const int Nvol = CommonParameters::Nvol();
95  const int Ndim = CommonParameters::Ndim();
96  const int Lx = CommonParameters::Lx();
97  const int Ly = CommonParameters::Ly();
98  const int Lz = CommonParameters::Lz();
99 
100  assert(Nvol == v1.nvol());
101  assert(Nvol == v2.nvol());
102  assert(time < Nt);
103  assert(momentum_sink.size() == Ndim - 1);
104 
105  Index_lex index;
106  std::vector<int> gm_index(Nd);
107  std::vector<double> corr_r(Nd), corr_i(Nd);
108 
109  for (int i = 0; i < Nd; ++i) {
110  gm_index[i] = gm_sink.index(i);
111  corr_r[i] = 0.0;
112  corr_i[i] = 0.0;
113  }
114 
115  static const double PI = 4.0 * atan(1.0);
116  std::vector<double> p_unit(Ndim - 1);
117  p_unit[0] = (2.0 * PI / Lx) * momentum_sink[0];
118  p_unit[1] = (2.0 * PI / Ly) * momentum_sink[1];
119  p_unit[2] = (2.0 * PI / Lz) * momentum_sink[2];
120 
121  std::vector<int> ipe(Ndim - 1);
122  ipe[0] = Communicator::ipe(0);
123  ipe[1] = Communicator::ipe(1);
124  ipe[2] = Communicator::ipe(2);
125 
126  for (int z = 0; z < Nz; ++z) {
127  for (int y = 0; y < Ny; ++y) {
128  for (int x = 0; x < Nx; ++x) {
129  int site = index.site(x, y, z, time);
130 
131  int x_global = x + ipe[0] * Nx;
132  int y_global = y + ipe[1] * Ny;
133  int z_global = z + ipe[2] * Nz;
134 
135  double p_x = p_unit[0] * (x_global - source_position[0]);
136  double p_y = p_unit[1] * (y_global - source_position[1]);
137  double p_z = p_unit[2] * (z_global - source_position[2]);
138 
139  double cos_p_xyz = cos(p_x + p_y + p_z);
140  double sin_p_xyz = sin(p_x + p_y + p_z);
141 
142  for (int s0 = 0; s0 < Nd; ++s0) {
143  int s1 = gm_index[s0];
144 
145  double v1_v2_r = 0.0;
146  double v1_v2_i = 0.0;
147 
148  for (int c1 = 0; c1 < Nc; ++c1) {
149  v1_v2_r += v1.cmp_r(c1, s1, site) * v2.cmp_r(c1, s0, site)
150  + v1.cmp_i(c1, s1, site) * v2.cmp_i(c1, s0, site);
151 
152  v1_v2_i += -v1.cmp_r(c1, s1, site) * v2.cmp_i(c1, s0, site)
153  + v1.cmp_i(c1, s1, site) * v2.cmp_r(c1, s0, site);
154  }
155 
156  //- corr[s0] += v2^dagger * v1 * exp(i * p_i * x_i)
157  corr_r[s0] += v1_v2_r * cos_p_xyz - v1_v2_i * sin_p_xyz;
158  corr_i[s0] += v1_v2_r * sin_p_xyz + v1_v2_i * cos_p_xyz;
159  }
160  }
161  }
162  }
163 
164  corr = cmplx(0.0, 0.0);
165  for (int s0 = 0; s0 < Nd; ++s0) {
166  corr += gm_sink.value(s0) * cmplx(corr_r[s0], corr_i[s0]);
167  }
168 }
169 
170 
171 //====================================================================
172 void contract_at_t(dcomplex& corr,
173  const GammaMatrix& gm_sink,
174  const int i_alpha,
175  const Field_F& v1, const Field_F& v2, const Field_F& v3,
176  const int time)
177 {
178 #if defined USE_GROUP_SU3
179  int Nc = CommonParameters::Nc();
180  int Nd = CommonParameters::Nd();
181  int Nx = CommonParameters::Nx();
182  int Ny = CommonParameters::Ny();
183  int Nz = CommonParameters::Nz();
184  int Nt = CommonParameters::Nt();
185  int Nvol = CommonParameters::Nvol();
186 
187  assert(Nc == 3);
188  assert(Nvol == v1.nvol());
189  assert(Nvol == v2.nvol());
190  assert(Nvol == v3.nvol());
191  assert(time < Nt);
192 
193  Index_lex index;
194  std::vector<int> gm_index(Nd);
195  std::vector<double> c_r(Nd), c_i(Nd);
196 
197  for (int i = 0; i < Nd; ++i) {
198  gm_index[i] = gm_sink.index(i);
199  c_r[i] = 0.0;
200  c_i[i] = 0.0;
201  }
202 
203  for (int z = 0; z < Nz; ++z) {
204  for (int y = 0; y < Ny; ++y) {
205  for (int x = 0; x < Nx; ++x) {
206  int site = index.site(x, y, z, time);
207 
208  for (int d1 = 0; d1 < Nd; ++d1) {
209  int d2 = gm_index[d1];
210  int d3 = i_alpha;
211 
212  c_r[d1] += (v1.cmp_r(C1, d1, site) * v2.cmp_r(C2, d2, site)
213  - v1.cmp_i(C1, d1, site) * v2.cmp_i(C2, d2, site)) * v3.cmp_r(C3, d3, site)
214  - (v1.cmp_r(C1, d1, site) * v2.cmp_i(C2, d2, site)
215  + v1.cmp_i(C1, d1, site) * v2.cmp_r(C2, d2, site)) * v3.cmp_i(C3, d3, site);
216  c_i[d1] += (v1.cmp_r(C1, d1, site) * v2.cmp_r(C2, d2, site)
217  - v1.cmp_i(C1, d1, site) * v2.cmp_i(C2, d2, site)) * v3.cmp_i(C3, d3, site)
218  + (v1.cmp_r(C1, d1, site) * v2.cmp_i(C2, d2, site)
219  + v1.cmp_i(C1, d1, site) * v2.cmp_r(C2, d2, site)) * v3.cmp_r(C3, d3, site);
220 
221  c_r[d1] += (v1.cmp_r(C2, d1, site) * v2.cmp_r(C3, d2, site)
222  - v1.cmp_i(C2, d1, site) * v2.cmp_i(C3, d2, site)) * v3.cmp_r(C1, d3, site)
223  - (v1.cmp_r(C2, d1, site) * v2.cmp_i(C3, d2, site)
224  + v1.cmp_i(C2, d1, site) * v2.cmp_r(C3, d2, site)) * v3.cmp_i(C1, d3, site);
225  c_i[d1] += (v1.cmp_r(C2, d1, site) * v2.cmp_r(C3, d2, site)
226  - v1.cmp_i(C2, d1, site) * v2.cmp_i(C3, d2, site)) * v3.cmp_i(C1, d3, site)
227  + (v1.cmp_r(C2, d1, site) * v2.cmp_i(C3, d2, site)
228  + v1.cmp_i(C2, d1, site) * v2.cmp_r(C3, d2, site)) * v3.cmp_r(C1, d3, site);
229 
230  c_r[d1] += (v1.cmp_r(C3, d1, site) * v2.cmp_r(C1, d2, site)
231  - v1.cmp_i(C3, d1, site) * v2.cmp_i(C1, d2, site)) * v3.cmp_r(C2, d3, site)
232  - (v1.cmp_r(C3, d1, site) * v2.cmp_i(C1, d2, site)
233  + v1.cmp_i(C3, d1, site) * v2.cmp_r(C1, d2, site)) * v3.cmp_i(C2, d3, site);
234  c_i[d1] += (v1.cmp_r(C3, d1, site) * v2.cmp_r(C1, d2, site)
235  - v1.cmp_i(C3, d1, site) * v2.cmp_i(C1, d2, site)) * v3.cmp_i(C2, d3, site)
236  + (v1.cmp_r(C3, d1, site) * v2.cmp_i(C1, d2, site)
237  + v1.cmp_i(C3, d1, site) * v2.cmp_r(C1, d2, site)) * v3.cmp_r(C2, d3, site);
238 
239  c_r[d1] -= (v1.cmp_r(C3, d1, site) * v2.cmp_r(C2, d2, site)
240  - v1.cmp_i(C3, d1, site) * v2.cmp_i(C2, d2, site)) * v3.cmp_r(C1, d3, site)
241  - (v1.cmp_r(C3, d1, site) * v2.cmp_i(C2, d2, site)
242  + v1.cmp_i(C3, d1, site) * v2.cmp_r(C2, d2, site)) * v3.cmp_i(C1, d3, site);
243  c_i[d1] -= (v1.cmp_r(C3, d1, site) * v2.cmp_r(C2, d2, site)
244  - v1.cmp_i(C3, d1, site) * v2.cmp_i(C2, d2, site)) * v3.cmp_i(C1, d3, site)
245  + (v1.cmp_r(C3, d1, site) * v2.cmp_i(C2, d2, site)
246  + v1.cmp_i(C3, d1, site) * v2.cmp_r(C2, d2, site)) * v3.cmp_r(C1, d3, site);
247 
248  c_r[d1] -= (v1.cmp_r(C2, d1, site) * v2.cmp_r(C1, d2, site)
249  - v1.cmp_i(C2, d1, site) * v2.cmp_i(C1, d2, site)) * v3.cmp_r(C3, d3, site)
250  - (v1.cmp_r(C2, d1, site) * v2.cmp_i(C1, d2, site)
251  + v1.cmp_i(C2, d1, site) * v2.cmp_r(C1, d2, site)) * v3.cmp_i(C3, d3, site);
252  c_i[d1] -= (v1.cmp_r(C2, d1, site) * v2.cmp_r(C1, d2, site)
253  - v1.cmp_i(C2, d1, site) * v2.cmp_i(C1, d2, site)) * v3.cmp_i(C3, d3, site)
254  + (v1.cmp_r(C2, d1, site) * v2.cmp_i(C1, d2, site)
255  + v1.cmp_i(C2, d1, site) * v2.cmp_r(C1, d2, site)) * v3.cmp_r(C3, d3, site);
256 
257  c_r[d1] -= (v1.cmp_r(C1, d1, site) * v2.cmp_r(C3, d2, site)
258  - v1.cmp_i(C1, d1, site) * v2.cmp_i(C3, d2, site)) * v3.cmp_r(C2, d3, site)
259  - (v1.cmp_r(C1, d1, site) * v2.cmp_i(C3, d2, site)
260  + v1.cmp_i(C1, d1, site) * v2.cmp_r(C3, d2, site)) * v3.cmp_i(C2, d3, site);
261  c_i[d1] -= (v1.cmp_r(C1, d1, site) * v2.cmp_r(C3, d2, site)
262  - v1.cmp_i(C1, d1, site) * v2.cmp_i(C3, d2, site)) * v3.cmp_i(C2, d3, site)
263  + (v1.cmp_r(C1, d1, site) * v2.cmp_i(C3, d2, site)
264  + v1.cmp_i(C1, d1, site) * v2.cmp_r(C3, d2, site)) * v3.cmp_r(C2, d3, site);
265  }
266  }
267  }
268  }
269 
270  corr = cmplx(0.0, 0.0);
271  for (int s0 = 0; s0 < Nd; ++s0) {
272  corr += gm_sink.value(s0) * cmplx(c_r[s0], c_i[s0]);
273  }
274 #endif // (USE_GROUP_SU3)
275 }
276 
277 
278 //====================================================================
279 //============================================================END=====
Field_F::cmp_i
double cmp_i(const int cc, const int s, const int site, const int e=0) const
Definition: field_F.h:101
Index_lex::site
int site(const int &x, const int &y, const int &z, const int &t) const
Definition: index_lex.h:53
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:34
Field::nvol
int nvol() const
Definition: field.h:116
Communicator::ipe
static int ipe(const int dir)
logical coordinate of current proc.
Definition: communicator.cpp:97
Field_F
Wilson-type fermion field.
Definition: field_F.h:37
GammaMatrix
Gamma Matrix class.
Definition: gammaMatrix.h:44
index_lex.h
Index_lex
Lexical site index.
Definition: index_lex.h:34
contract_4spinor.h
GammaMatrix::value
dcomplex value(int row) const
Definition: gammaMatrix.h:88
GammaMatrix::index
int index(int row) const
Definition: gammaMatrix.h:83
Field_F::cmp_r
double cmp_r(const int cc, const int s, const int site, const int e=0) const
Definition: field_F.h:95