25 bool init = Projection::Factory::Register(
"Stout_SU3", create_object);
36 const string str_vlevel = params.
get_string(
"verbose_level");
76 int Nvol = Uorg.
nvol();
78 assert(Cst.
nex() == Nex);
79 assert(Cst.
nvol() == Nvol);
80 assert(U.
nex() == Nex);
81 assert(U.
nvol() == Nvol);
83 int NinG = Uorg.
nin();
92 for (
int mu = 0; mu < Nex; ++mu) {
93 for (
int site = 0; site < Nvol; ++site) {
94 Uorg.
mat(ut, site, mu);
95 Cst.
mat(ct, site, mu);
98 iQ2.mult_nn(iQ1, iQ1);
101 double norm = iQ1.norm2();
102 if (norm > 1.0e-10) {
106 for (
int cc = 0; cc <
NC *
NC; ++cc) {
107 dcomplex qt = f0 * cmplx(iQ0.r(cc), iQ0.i(cc))
108 + f1 * cmplx(iQ1.i(cc), -iQ1.r(cc))
109 - f2 * cmplx(iQ2.r(cc), iQ2.i(cc));
110 e_iQ.
set_r(cc, real(qt));
111 e_iQ.
set_i(cc, imag(qt));
138 int Nvol = iQ.
nvol();
148 for (
int mu = 0; mu < Nex; ++mu) {
149 for (
int site = 0; site < Nvol; ++site) {
150 iQ1 = iQ.
mat(site, mu);
154 double norm = iQ1.
norm2();
155 if (norm > 1.0e-10) {
159 for (
int cc = 0; cc <
NC *
NC; ++cc) {
160 dcomplex qt = f0 * cmplx(iQ0.r(cc), iQ0.i(cc))
161 + f1 * cmplx(iQ1.i(cc), -iQ1.r(cc))
162 - f2 * cmplx(iQ2.r(cc), iQ2.i(cc));
163 e_iQ.
set_ri(cc, site, mu, real(qt), imag(qt));
190 double alpha,
const Field_G& Sigmap,
203 assert(Xi.
nvol() == Nvol);
204 assert(iTheta.
nvol() == Nvol);
205 assert(Sigmap.
nvol() == Nvol);
206 assert(Cst.
nvol() == Nvol);
207 assert(Uorg.
nvol() == Nvol);
208 assert(iTheta.
nex() == Nex);
209 assert(Sigmap.
nex() == Nex);
210 assert(Cst.
nex() == Nex);
211 assert(Uorg.
nex() == Nex);
224 for (
int mu = 0; mu < Nex; ++mu) {
225 for (
int site = 0; site < Nvol; ++site) {
227 Cst.
mat(C_tmp, site, mu);
229 Uorg.
mat(U_tmp, site, mu);
231 Sigmap.
mat(Sigmap_tmp, site, mu);
236 iQ2.mult_nn(iQ1, iQ1);
237 iQ3.mult_nn(iQ1, iQ2);
240 double norm = iQ1.norm2();
241 if (norm > 1.0e-10) {
245 for (
int cc = 0; cc <
NC *
NC; ++cc) {
246 dcomplex qt = f0 * cmplx(iQ0.r(cc), iQ0.i(cc))
247 + f1 * cmplx(iQ1.i(cc), -iQ1.r(cc))
248 - f2 * cmplx(iQ2.r(cc), iQ2.i(cc));
249 e_iQ.
set(cc, real(qt), imag(qt));
256 double cos_w = cos(w);
258 dcomplex emiu = cmplx(cos(u), -sin(u));
259 dcomplex e2iu = cmplx(cos(2.0 * u), sin(2.0 * u));
261 dcomplex r01 = cmplx(2.0 * u, 2.0 * (u2 - w2)) * e2iu
262 + emiu * cmplx(16.0 * u * cos_w + 2.0 * u * (3.0 * u2 + w2) * xi0,
263 -8.0 * u2 * cos_w + 2.0 * (9.0 * u2 + w2) * xi0);
265 dcomplex r11 = cmplx(2.0, 4.0 * u) * e2iu
266 + emiu * cmplx(-2.0 * cos_w + (3.0 * u2 - w2) * xi0,
267 2.0 * u * cos_w + 6.0 * u * xi0);
269 dcomplex r21 = cmplx(0.0, 2.0) * e2iu
270 + emiu * cmplx(-3.0 * u * xi0, cos_w - 3.0 * xi0);
272 dcomplex r02 = cmplx(-2.0, 0.0) * e2iu
273 + emiu * cmplx(-8.0 * u2 * xi0,
274 2.0 * u * (cos_w + xi0 + 3.0 * u2 * xi1));
276 dcomplex r12 = emiu * cmplx(2.0 * u * xi0,
277 -cos_w - xi0 + 3.0 * u2 * xi1);
279 dcomplex r22 = emiu * cmplx(xi0, -3.0 * u * xi1);
281 double fden = 1.0 / (2 * (9.0 * u2 - w2) * (9.0 * u2 - w2));
283 dcomplex b10 = cmplx(2.0 * u, 0.0) * r01 + cmplx(3.0 * u2 - w2, 0.0) * r02
284 - cmplx(30.0 * u2 + 2.0 * w2, 0.0) * f0;
285 dcomplex b11 = cmplx(2.0 * u, 0.0) * r11 + cmplx(3.0 * u2 - w2, 0.0) * r12
286 - cmplx(30.0 * u2 + 2.0 * w2, 0.0) * f1;
287 dcomplex b12 = cmplx(2.0 * u, 0.0) * r21 + cmplx(3.0 * u2 - w2, 0.0) * r22
288 - cmplx(30.0 * u2 + 2.0 * w2, 0.0) * f2;
290 dcomplex b20 = r01 - cmplx(3.0 * u, 0.0) * r02 - cmplx(24.0 * u, 0.0) * f0;
291 dcomplex b21 = r11 - cmplx(3.0 * u, 0.0) * r12 - cmplx(24.0 * u, 0.0) * f1;
292 dcomplex b22 = r21 - cmplx(3.0 * u, 0.0) * r22 - cmplx(24.0 * u, 0.0) * f2;
294 b10 *= cmplx(fden, 0.0);
295 b11 *= cmplx(fden, 0.0);
296 b12 *= cmplx(fden, 0.0);
297 b20 *= cmplx(fden, 0.0);
298 b21 *= cmplx(fden, 0.0);
299 b22 *= cmplx(fden, 0.0);
301 for (
int cc = 0; cc < NC *
NC; ++cc) {
302 dcomplex qt1 = b10 * cmplx(iQ0.r(cc), iQ0.i(cc))
303 + b11 * cmplx(iQ1.i(cc), -iQ1.r(cc))
304 - b12 * cmplx(iQ2.r(cc), iQ2.i(cc));
305 B1.set(cc, real(qt1), imag(qt1));
307 dcomplex qt2 = b20 * cmplx(iQ0.r(cc), iQ0.i(cc))
308 + b21 * cmplx(iQ1.i(cc), -iQ1.r(cc))
309 - b22 * cmplx(iQ2.r(cc), iQ2.i(cc));
310 B2.
set(cc, real(qt2), imag(qt2));
313 USigmap.mult_nn(U_tmp, Sigmap_tmp);
315 tmp1.mult_nn(USigmap, B1);
318 dcomplex tr1 = cmplx(tmp1.r(0) + tmp1.r(4) + tmp1.r(8),
319 tmp1.i(0) + tmp1.i(4) + tmp1.i(8));
320 dcomplex tr2 = cmplx(tmp2.
r(0) + tmp2.
r(4) + tmp2.
r(8),
321 tmp2.
i(0) + tmp2.
i(4) + tmp2.
i(8));
323 iQUS.mult_nn(iQ1, USigmap);
324 iUSQ.mult_nn(USigmap, iQ1);
326 for (
int cc = 0; cc < NC *
NC; ++cc) {
327 dcomplex qt = tr1 * cmplx(iQ1.i(cc), -iQ1.r(cc))
328 - tr2 * cmplx(iQ2.r(cc), iQ2.i(cc))
329 + f1 * cmplx(USigmap.r(cc), USigmap.i(cc))
330 + f2 * cmplx(iQUS.i(cc), -iQUS.r(cc))
331 + f2 * cmplx(iUSQ.i(cc), -iUSQ.r(cc));
332 iGamma.
set(cc, -imag(qt), real(qt));
342 iTheta_tmp.
mult_nn(iGamma, U_tmp);
344 iTheta.
set_mat(site, mu, iTheta_tmp);
346 Xi_tmp.mult_nn(Sigmap_tmp, e_iQ);
347 Xi_tmp.multadd_dn(C_tmp, iGamma);
367 const double& u,
const double& w)
372 double cos_w = cos(w);
374 double cos_u = cos(u);
375 double sin_u = sin(u);
377 dcomplex emiu = cmplx(cos_u, -sin_u);
378 dcomplex e2iu = cmplx(cos_u * cos_u - sin_u * sin_u, 2.0 * sin_u * cos_u);
380 dcomplex h0 = e2iu * cmplx(u2 - w2, 0.0)
381 + emiu * cmplx(8.0 * u2 * cos_w, 2.0 * u * (3.0 * u2 + w2) * xi0);
382 dcomplex h1 = cmplx(2 * u, 0.0) * e2iu
383 - emiu * cmplx(2.0 * u * cos_w, -(3.0 * u2 - w2) * xi0);
384 dcomplex h2 = e2iu - emiu * cmplx(cos_w, 3.0 * u * xi0);
386 double fden = 1.0 / (9.0 * u2 - w2);
399 double c0 = -(iQ3.
i(0, 0) + iQ3.
i(1, 1) + iQ3.
i(2, 2)) / 3.0;
400 double c1 = -0.5 * (iQ2.
r(0, 0) + iQ2.
r(1, 1) + iQ2.
r(2, 2));
401 double c13r = sqrt(c1 / 3.0);
402 double c0max = 2.0 * c13r * c13r * c13r;
404 double theta = acos(c0 / c0max);
407 u = c13r * cos(theta / 3.0);
408 w = sqrt(c1) * sin(theta / 3.0);
428 static double c0 = -1.0 / 3.0;
429 static double c1 = 1.0 / 30.0;
430 static double c2 = -1.0 / 840.0;
431 static double c3 = 1.0 / 45360.0;
432 static double c4 = -1.0 / 3991680.0;
434 return c0 + w2 * (c1 + w2 * (c2 + w2 * (c3 + w2 * c4)));
436 return (w * cos(w) - sin(w)) / (w * w * w);
448 int Nvol = iQ.
nvol();
454 for (
int ex = 0; ex < Nex; ++ex) {
455 for (
int site = 0; site < Nvol; ++site) {
458 h1 = iQ.
mat(site, ex);
460 for (
int iprec = 0; iprec < Nprec; ++iprec) {
461 double exf = 1.0 / (Nprec - iprec);
void general(const char *format,...)
void force_recursive(Field_G &Xi, Field_G &iTheta, double alpha, const Field_G &Sigmap, const Field_G &C, const Field_G &U)
determination of fields for force calculation
void set_uw(double &u, double &w, const Mat_SU_N &iQ1, const Mat_SU_N &iQ2)
Mat_SU_N & at()
antihermitian traceless
void set_i(int c, const double &im)
double func_xi0(double w)
void project(Field_G &U, double alpha, const Field_G &C, const Field_G &Uorg)
projection U = P[alpha, C, Uorg]
void set_r(int c, const double &re)
double func_xi1(double w)
void mult_nd(const Mat_SU_N &u1, const Mat_SU_N &u2)
void set_fj(dcomplex &f0, dcomplex &f1, dcomplex &f2, const double &u, const double &w)
base class for projection operator into gauge group.
void set_parameters(const Parameters ¶m)
void exp_iQ(Field_G &e_iQ, const Field_G &iQ)
void exp_iQ_bf(Field_G &e_iQ, const Field_G &iQ)
static const std::string class_name
static double get_time()
obtain a wall-clock time.
void set(int c, double re, const double &im)
string get_string(const string &key) const
void set_mat(const int site, const int mn, const Mat_SU_N &U)
Stout(exponential)-type projection to SU(N) gauge group.
Mat_SU_N mat(const int site, const int mn=0) const
void mult_nn(const Mat_SU_N &u1, const Mat_SU_N &u2)
static VerboseLevel set_verbose_level(const std::string &str)
Bridge::VerboseLevel m_vl
void set_ri(const int cc, const int site, const int mn, const double re, const double im)