Bridge++  Ver. 1.3.x
contract_4spinor.cpp
Go to the documentation of this file.
1 
14 #include "contract_4spinor.h"
15 #include "index_lex.h"
16 
17 #include <cassert>
18 
19 #define NC 3
20 #define C1 0
21 #define C2 1
22 #define C3 2
23 
24 //====================================================================
25 void contract_at_t(dcomplex& corr, const GammaMatrix& gm,
26  const Field_F& v1, const Field_F& v2, int time)
27 {
28  int Nc = CommonParameters::Nc();
29  int Nd = CommonParameters::Nd();
30  int Nx = CommonParameters::Nx();
31  int Ny = CommonParameters::Ny();
32  int Nz = CommonParameters::Nz();
33  int Nt = CommonParameters::Nt();
34  int Nvol = CommonParameters::Nvol();
35 
36  assert(Nvol == v1.nvol());
37  assert(Nvol == v2.nvol());
38  assert(time < Nt);
39 
40  Index_lex index;
41  std::vector<int> gmindex(Nd);
42  std::vector<double> corr_r(Nd), corr_i(Nd);
43 
44  for (int i = 0; i < Nd; ++i) {
45  gmindex[i] = gm.index(i);
46  corr_r[i] = 0.0;
47  corr_i[i] = 0.0;
48  }
49 
50  for (int z = 0; z < Nz; ++z) {
51  for (int y = 0; y < Ny; ++y) {
52  for (int x = 0; x < Nx; ++x) {
53  int site = index.site(x, y, z, time);
54 
55  for (int s0 = 0; s0 < Nd; ++s0) {
56  int s1 = gmindex[s0];
57 
58  for (int c1 = 0; c1 < Nc; ++c1) {
59  corr_r[s0] += v1.cmp_r(c1, s1, site) * v2.cmp_r(c1, s0, site)
60  + v1.cmp_i(c1, s1, site) * v2.cmp_i(c1, s0, site);
61 
62  corr_i[s0] += v1.cmp_r(c1, s1, site) * v2.cmp_i(c1, s0, site)
63  - v1.cmp_i(c1, s1, site) * v2.cmp_r(c1, s0, site);
64  }
65  }
66  }
67  }
68  }
69 
70  corr = cmplx(0.0, 0.0);
71  for (int s0 = 0; s0 < Nd; ++s0) {
72  corr += gm.value(s0) * cmplx(corr_r[s0], corr_i[s0]);
73  }
74 }
75 
76 
77 //====================================================================
78 void contract_at_t(dcomplex& corr, const GammaMatrix& gm, int d3,
79  const Field_F& v1, const Field_F& v2, const Field_F& v3,
80  int time)
81 {
82  int Nc = CommonParameters::Nc();
83  int Nd = CommonParameters::Nd();
84  int Nx = CommonParameters::Nx();
85  int Ny = CommonParameters::Ny();
86  int Nz = CommonParameters::Nz();
87  int Nt = CommonParameters::Nt();
88  int Nvol = CommonParameters::Nvol();
89 
90  assert(Nc == 3);
91  assert(Nvol == v1.nvol());
92  assert(Nvol == v2.nvol());
93  assert(Nvol == v3.nvol());
94  assert(time < Nt);
95 
96  Index_lex index;
97  std::vector<int> gmindex(Nd);
98  std::vector<double> cr(Nd), ci(Nd);
99 
100  for (int i = 0; i < Nd; ++i) {
101  gmindex[i] = gm.index(i);
102  cr[i] = 0.0;
103  ci[i] = 0.0;
104  }
105 
106  for (int z = 0; z < Nz; ++z) {
107  for (int y = 0; y < Ny; ++y) {
108  for (int x = 0; x < Nx; ++x) {
109  int site = index.site(x, y, z, time);
110 
111  for (int d1 = 0; d1 < Nd; ++d1) {
112  int d2 = gmindex[d1];
113 
114  cr[d1] += (v1.cmp_r(C1, d1, site) * v2.cmp_r(C2, d2, site)
115  - v1.cmp_i(C1, d1, site) * v2.cmp_i(C2, d2, site)) * v3.cmp_r(C3, d3, site)
116  - (v1.cmp_r(C1, d1, site) * v2.cmp_i(C2, d2, site)
117  + v1.cmp_i(C1, d1, site) * v2.cmp_r(C2, d2, site)) * v3.cmp_i(C3, d3, site);
118  ci[d1] += (v1.cmp_r(C1, d1, site) * v2.cmp_r(C2, d2, site)
119  - v1.cmp_i(C1, d1, site) * v2.cmp_i(C2, d2, site)) * v3.cmp_i(C3, d3, site)
120  + (v1.cmp_r(C1, d1, site) * v2.cmp_i(C2, d2, site)
121  + v1.cmp_i(C1, d1, site) * v2.cmp_r(C2, d2, site)) * v3.cmp_r(C3, d3, site);
122 
123  cr[d1] += (v1.cmp_r(C2, d1, site) * v2.cmp_r(C3, d2, site)
124  - v1.cmp_i(C2, d1, site) * v2.cmp_i(C3, d2, site)) * v3.cmp_r(C1, d3, site)
125  - (v1.cmp_r(C2, d1, site) * v2.cmp_i(C3, d2, site)
126  + v1.cmp_i(C2, d1, site) * v2.cmp_r(C3, d2, site)) * v3.cmp_i(C1, d3, site);
127  ci[d1] += (v1.cmp_r(C2, d1, site) * v2.cmp_r(C3, d2, site)
128  - v1.cmp_i(C2, d1, site) * v2.cmp_i(C3, d2, site)) * v3.cmp_i(C1, d3, site)
129  + (v1.cmp_r(C2, d1, site) * v2.cmp_i(C3, d2, site)
130  + v1.cmp_i(C2, d1, site) * v2.cmp_r(C3, d2, site)) * v3.cmp_r(C1, d3, site);
131 
132  cr[d1] += (v1.cmp_r(C3, d1, site) * v2.cmp_r(C1, d2, site)
133  - v1.cmp_i(C3, d1, site) * v2.cmp_i(C1, d2, site)) * v3.cmp_r(C2, d3, site)
134  - (v1.cmp_r(C3, d1, site) * v2.cmp_i(C1, d2, site)
135  + v1.cmp_i(C3, d1, site) * v2.cmp_r(C1, d2, site)) * v3.cmp_i(C2, d3, site);
136  ci[d1] += (v1.cmp_r(C3, d1, site) * v2.cmp_r(C1, d2, site)
137  - v1.cmp_i(C3, d1, site) * v2.cmp_i(C1, d2, site)) * v3.cmp_i(C2, d3, site)
138  + (v1.cmp_r(C3, d1, site) * v2.cmp_i(C1, d2, site)
139  + v1.cmp_i(C3, d1, site) * v2.cmp_r(C1, d2, site)) * v3.cmp_r(C2, d3, site);
140 
141  cr[d1] -= (v1.cmp_r(C3, d1, site) * v2.cmp_r(C2, d2, site)
142  - v1.cmp_i(C3, d1, site) * v2.cmp_i(C2, d2, site)) * v3.cmp_r(C1, d3, site)
143  - (v1.cmp_r(C3, d1, site) * v2.cmp_i(C2, d2, site)
144  + v1.cmp_i(C3, d1, site) * v2.cmp_r(C2, d2, site)) * v3.cmp_i(C1, d3, site);
145  ci[d1] -= (v1.cmp_r(C3, d1, site) * v2.cmp_r(C2, d2, site)
146  - v1.cmp_i(C3, d1, site) * v2.cmp_i(C2, d2, site)) * v3.cmp_i(C1, d3, site)
147  + (v1.cmp_r(C3, d1, site) * v2.cmp_i(C2, d2, site)
148  + v1.cmp_i(C3, d1, site) * v2.cmp_r(C2, d2, site)) * v3.cmp_r(C1, d3, site);
149 
150  cr[d1] -= (v1.cmp_r(C2, d1, site) * v2.cmp_r(C1, d2, site)
151  - v1.cmp_i(C2, d1, site) * v2.cmp_i(C1, d2, site)) * v3.cmp_r(C3, d3, site)
152  - (v1.cmp_r(C2, d1, site) * v2.cmp_i(C1, d2, site)
153  + v1.cmp_i(C2, d1, site) * v2.cmp_r(C1, d2, site)) * v3.cmp_i(C3, d3, site);
154  ci[d1] -= (v1.cmp_r(C2, d1, site) * v2.cmp_r(C1, d2, site)
155  - v1.cmp_i(C2, d1, site) * v2.cmp_i(C1, d2, site)) * v3.cmp_i(C3, d3, site)
156  + (v1.cmp_r(C2, d1, site) * v2.cmp_i(C1, d2, site)
157  + v1.cmp_i(C2, d1, site) * v2.cmp_r(C1, d2, site)) * v3.cmp_r(C3, d3, site);
158 
159  cr[d1] -= (v1.cmp_r(C1, d1, site) * v2.cmp_r(C3, d2, site)
160  - v1.cmp_i(C1, d1, site) * v2.cmp_i(C3, d2, site)) * v3.cmp_r(C2, d3, site)
161  - (v1.cmp_r(C1, d1, site) * v2.cmp_i(C3, d2, site)
162  + v1.cmp_i(C1, d1, site) * v2.cmp_r(C3, d2, site)) * v3.cmp_i(C2, d3, site);
163  ci[d1] -= (v1.cmp_r(C1, d1, site) * v2.cmp_r(C3, d2, site)
164  - v1.cmp_i(C1, d1, site) * v2.cmp_i(C3, d2, site)) * v3.cmp_i(C2, d3, site)
165  + (v1.cmp_r(C1, d1, site) * v2.cmp_i(C3, d2, site)
166  + v1.cmp_i(C1, d1, site) * v2.cmp_r(C3, d2, site)) * v3.cmp_r(C2, d3, site);
167  }
168  }
169  }
170  }
171 
172  corr = cmplx(0.0, 0.0);
173  for (int s0 = 0; s0 < Nd; ++s0) {
174  corr += gm.value(s0) * cmplx(cr[s0], ci[s0]);
175  }
176 }
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
C3
#define C3
Definition: contract_4spinor.cpp:22
Field::nvol
int nvol() const
Definition: field.h:116
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
C1
#define C1
Definition: contract_4spinor.cpp:20
contract_at_t
void contract_at_t(dcomplex &corr, const GammaMatrix &gm, const Field_F &v1, const Field_F &v2, int time)
contraction for meson at a given time t.
Definition: contract_4spinor.cpp:30
contract_4spinor.h
GammaMatrix::value
dcomplex value(int row) const
Definition: gammaMatrix.h:88
C2
#define C2
Definition: contract_4spinor.cpp:21
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