Bridge++  Ver. 1.2.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
contract_4spinor.cpp
Go to the documentation of this file.
1 #include "contract_4spinor.h"
2 
3 #include <cassert>
4 
5 // This implementation only applies to SU(3) group and Nd=4 case.
6 #define NC 3
7 #define NC2 6
8 #define NDF 18
9 #define ND 4
10 #define NCD 12
11 #define NCD2 24
12 #define C1 0
13 #define C2 2
14 #define C3 4
15 
16 //====================================================================
17 void contract_at_t(dcomplex& corr, const GammaMatrix& gm,
18  const Field_F& v1, const Field_F& v2, int time)
19 {
20  int Nvol = v1.nvol();
21  int Nvol_s = Nvol / CommonParameters::Nt();
22 
23  assert(Nvol == v2.nvol());
24  assert(v1.nex() == 1);
25  assert(v2.nex() == 1);
26 
27  double *w1 = const_cast<Field_F *>(&v1)->ptr(0);
28  double *w2 = const_cast<Field_F *>(&v2)->ptr(0);
29 
30  int id11 = 0;
31  int id12 = NC2;
32  int id13 = 2 * NC2;
33  int id14 = 3 * NC2;
34  int id21 = gm.index(0) * NC2;
35  int id22 = gm.index(1) * NC2;
36  int id23 = gm.index(2) * NC2;
37  int id24 = gm.index(3) * NC2;
38 
39  double cr[ND];
40  double ci[ND];
41  for (int id = 0; id < ND; ++id) {
42  cr[id] = 0.0;
43  ci[id] = 0.0;
44  }
45 
46  for (int ss = 0; ss < Nvol_s; ++ss) {
47  int site = NCD2 * (ss + time * Nvol_s);
48 
49  for (int cc = 0; cc < NC; ++cc) {
50  cr[0] += w1[2 * cc + id21 + site] * w2[2 * cc + id11 + site]
51  + w1[2 * cc + 1 + id21 + site] * w2[2 * cc + 1 + id11 + site];
52  cr[1] += w1[2 * cc + id22 + site] * w2[2 * cc + id12 + site]
53  + w1[2 * cc + 1 + id22 + site] * w2[2 * cc + 1 + id12 + site];
54  cr[2] += w1[2 * cc + id23 + site] * w2[2 * cc + id13 + site]
55  + w1[2 * cc + 1 + id23 + site] * w2[2 * cc + 1 + id13 + site];
56  cr[3] += w1[2 * cc + id24 + site] * w2[2 * cc + id14 + site]
57  + w1[2 * cc + 1 + id24 + site] * w2[2 * cc + 1 + id14 + site];
58 
59  ci[0] += w1[2 * cc + id21 + site] * w2[2 * cc + 1 + id11 + site]
60  - w1[2 * cc + 1 + id21 + site] * w2[2 * cc + id11 + site];
61  ci[1] += w1[2 * cc + id22 + site] * w2[2 * cc + 1 + id12 + site]
62  - w1[2 * cc + 1 + id22 + site] * w2[2 * cc + id12 + site];
63  ci[2] += w1[2 * cc + id23 + site] * w2[2 * cc + 1 + id13 + site]
64  - w1[2 * cc + 1 + id23 + site] * w2[2 * cc + id13 + site];
65  ci[3] += w1[2 * cc + id24 + site] * w2[2 * cc + 1 + id14 + site]
66  - w1[2 * cc + 1 + id24 + site] * w2[2 * cc + id14 + site];
67  }
68  }
69 
70  corr = gm.value(0) * cmplx(cr[0], ci[0])
71  + gm.value(1) * cmplx(cr[1], ci[1])
72  + gm.value(2) * cmplx(cr[2], ci[2])
73  + gm.value(3) * cmplx(cr[3], ci[3]);
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 Nvol = v1.nvol();
83  int Nvol_s = Nvol / CommonParameters::Nt();
84 
85  assert(Nvol == v2.nvol());
86  assert(Nvol == v3.nvol());
87  assert(v1.nex() == 1);
88  assert(v2.nex() == 1);
89  assert(v3.nex() == 1);
90 
91  double *w1 = const_cast<Field_F *>(&v1)->ptr(0);
92  double *w2 = const_cast<Field_F *>(&v2)->ptr(0);
93  double *w3 = const_cast<Field_F *>(&v3)->ptr(0);
94 
95  int s3 = d3 * NC2;
96 
97  int gmd[ND];
98  double cr[ND];
99  double ci[ND];
100 
101  for (int id = 0; id < ND; ++id) {
102  gmd[id] = gm.index(id);
103  cr[id] = 0.0;
104  ci[id] = 0.0;
105  }
106 
107 
108  for (int ss = 0; ss < Nvol_s; ++ss) {
109  int site = NCD2 * (ss + time * Nvol_s);
110 
111  for (int id = 0; id < ND; ++id) {
112  int iw1 = id * NC2 + site;
113  int iw2 = gmd[id] * NC2 + site;
114  int iw3 = s3 + site;
115 
116  cr[id] += (w1[C1 + iw1] * w2[C2 + iw2]
117  - w1[C1 + 1 + iw1] * w2[C2 + 1 + iw2]) * w3[C3 + iw3]
118  - (w1[C1 + iw1] * w2[C2 + 1 + iw2]
119  + w1[C1 + 1 + iw1] * w2[C2 + iw2]) * w3[C3 + 1 + iw3];
120  ci[id] += (w1[C1 + iw1] * w2[C2 + iw2]
121  - w1[C1 + 1 + iw1] * w2[C2 + 1 + iw2]) * w3[C3 + 1 + iw3]
122  + (w1[C1 + iw1] * w2[C2 + 1 + iw2]
123  + w1[C1 + 1 + iw1] * w2[C2 + iw2]) * w3[C3 + iw3];
124 
125  cr[id] += (w1[C2 + iw1] * w2[C3 + iw2]
126  - w1[C2 + 1 + iw1] * w2[C3 + 1 + iw2]) * w3[C1 + iw3]
127  - (w1[C2 + iw1] * w2[C3 + 1 + iw2]
128  + w1[C2 + 1 + iw1] * w2[C3 + iw2]) * w3[C1 + 1 + iw3];
129  ci[id] += (w1[C2 + iw1] * w2[C3 + iw2]
130  - w1[C2 + 1 + iw1] * w2[C3 + 1 + iw2]) * w3[C1 + 1 + iw3]
131  + (w1[C2 + iw1] * w2[C3 + 1 + iw2]
132  + w1[C2 + 1 + iw1] * w2[C3 + iw2]) * w3[C1 + iw3];
133 
134  cr[id] += (w1[C3 + iw1] * w2[C1 + iw2]
135  - w1[C3 + 1 + iw1] * w2[C1 + 1 + iw2]) * w3[C2 + iw3]
136  - (w1[C3 + iw1] * w2[C1 + 1 + iw2]
137  + w1[C3 + 1 + iw1] * w2[C1 + iw2]) * w3[C2 + 1 + iw3];
138  ci[id] += (w1[C3 + iw1] * w2[C1 + iw2]
139  - w1[C3 + 1 + iw1] * w2[C1 + 1 + iw2]) * w3[C2 + 1 + iw3]
140  + (w1[C3 + iw1] * w2[C1 + 1 + iw2]
141  + w1[C3 + 1 + iw1] * w2[C1 + iw2]) * w3[C2 + iw3];
142 
143  cr[id] -= (w1[C3 + iw1] * w2[C2 + iw2]
144  - w1[C3 + 1 + iw1] * w2[C2 + 1 + iw2]) * w3[C1 + iw3]
145  - (w1[C3 + iw1] * w2[C2 + 1 + iw2]
146  + w1[C3 + 1 + iw1] * w2[C2 + iw2]) * w3[C1 + 1 + iw3];
147  ci[id] -= (w1[C3 + iw1] * w2[C2 + iw2]
148  - w1[C3 + 1 + iw1] * w2[C2 + 1 + iw2]) * w3[C1 + 1 + iw3]
149  + (w1[C3 + iw1] * w2[C2 + 1 + iw2]
150  + w1[C3 + 1 + iw1] * w2[C2 + iw2]) * w3[C1 + iw3];
151 
152  cr[id] -= (w1[C2 + iw1] * w2[C1 + iw2]
153  - w1[C2 + 1 + iw1] * w2[C1 + 1 + iw2]) * w3[C3 + iw3]
154  - (w1[C2 + iw1] * w2[C1 + 1 + iw2]
155  + w1[C2 + 1 + iw1] * w2[C1 + iw2]) * w3[C3 + 1 + iw3];
156  ci[id] -= (w1[C2 + iw1] * w2[C1 + iw2]
157  - w1[C2 + 1 + iw1] * w2[C1 + 1 + iw2]) * w3[C3 + 1 + iw3]
158  + (w1[C2 + iw1] * w2[C1 + 1 + iw2]
159  + w1[C2 + 1 + iw1] * w2[C1 + iw2]) * w3[C3 + iw3];
160 
161  cr[id] -= (w1[C1 + iw1] * w2[C3 + iw2]
162  - w1[C1 + 1 + iw1] * w2[C3 + 1 + iw2]) * w3[C2 + iw3]
163  - (w1[C1 + iw1] * w2[C3 + 1 + iw2]
164  + w1[C1 + 1 + iw1] * w2[C3 + iw2]) * w3[C2 + 1 + iw3];
165  ci[id] -= (w1[C1 + iw1] * w2[C3 + iw2]
166  - w1[C1 + 1 + iw1] * w2[C3 + 1 + iw2]) * w3[C2 + 1 + iw3]
167  + (w1[C1 + iw1] * w2[C3 + 1 + iw2]
168  + w1[C1 + 1 + iw1] * w2[C3 + iw2]) * w3[C2 + iw3];
169  }
170  }
171  corr = cmplx(0.0, 0.0);
172  for (int id = 0; id < ND; ++id) {
173  corr += cmplx(cr[id], ci[id]) * gm.value(id);
174  }
175 }
contract_at_t
void contract_at_t(dcomplex &corr, const GammaMatrix &gm, const Field_F &f1, const Field_F &f2, int time)
contraction for meson at a given time t.
Definition: contract_4spinor.cpp:17
NC
#define NC
Definition: contract_4spinor.cpp:6
Field::nvol
int nvol() const
Definition: field.h:101
NCD2
#define NCD2
Definition: contract_4spinor.cpp:11
Field_F
Wilson-type fermion field.
Definition: field_F.h:37
NC2
#define NC2
Definition: contract_4spinor.cpp:7
GammaMatrix
Gamma Matrix class.
Definition: gammaMatrix.h:45
C1
#define C1
Definition: contract_4spinor.cpp:12
Field::nex
int nex() const
Definition: field.h:102
C2
#define C2
Definition: contract_4spinor.cpp:13
ND
#define ND
Definition: contract_4spinor.cpp:9
contract_4spinor.h
C3
#define C3
Definition: contract_4spinor.cpp:14
GammaMatrix::value
dcomplex value(int row) const
Definition: gammaMatrix.h:89
GammaMatrix::index
int index(int row) const
Definition: gammaMatrix.h:84