Bridge++  Version 1.5.4
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
fieldStrength.cpp
Go to the documentation of this file.
1 
14 #include "fieldStrength.h"
15 
16 const std::string FieldStrength::class_name = "FieldStrengh";
17 
18 //====================================================================
19 void FieldStrength::construct_Fmunu_1x1(Field_G& Fmunu_1x1,
20  const int mu, const int nu, const Field_G& U)
21 {
22  const int Nvol = CommonParameters::Nvol();
23 
24  //- building blocks
25  // (1) mu (2)
26  // +->-+
27  // nu | |
28  // i+ +
29  Field_G Cup;
30 
31  m_staple.upper(Cup, U, mu, nu);
32 
33  // (1) mu (2)
34  // i+ +
35  // nu | |
36  // +->-+
37  Field_G Cdn;
38  m_staple.lower(Cdn, U, mu, nu);
39 
40  Field_G Umu;
41  Umu.setpart_ex(0, U, mu);
42 
43  //---------------------------------
44  //- right clover leaf
45  // (2) +-<-+
46  // nu | |
47  // i+->-+
48  // mu (1)
49  mult_Field_Gnd(Fmunu_1x1, 0, Umu, 0, Cup, 0);
50 
51  // (1) i+-<-+
52  // nu | |
53  // +->-+
54  // mu (2)
55  multadd_Field_Gnd(Fmunu_1x1, 0, Cdn, 0, Umu, 0, 1.0);
56  //---------------------------------
57 
58  //---------------------------------
59  //- left clover leaf
60  // mu (2)
61  // +-<-+ (1)
62  // | | nu mu (1)
63  // +->-+i + +-<-+i (2)
64  // | | nu
65  // +->-+
66  Field_G v;
67  mult_Field_Gdn(v, 0, Cup, 0, Umu, 0);
68  multadd_Field_Gdn(v, 0, Umu, 0, Cdn, 0, 1.0);
69 
70  // NB. shift.forward(mu)
71  // mu (2) mu (2)
72  // +-<-+ (1) +-<-+ (1)
73  // | | nu -> | | nu
74  // i+->-+ +->-+i
75  Field_G v2;
76  m_shift.forward(v2, v, mu);
77 
78  axpy(Fmunu_1x1, 1.0, v2);
79  //---------------------------------
80 
81  ah_Field_G(Fmunu_1x1, 0);
82 
83  scal(Fmunu_1x1, -0.25);
84  Fmunu_1x1.xI();
85 }
86 
87 
88 //====================================================================
89 void FieldStrength::construct_Fmunu_1x2(Field_G& Fmunu_1x2,
90  const int mu, const int nu, const Field_G& U)
91 {
92  const int Nvol = CommonParameters::Nvol();
93 
94 
95  //-- building blocks
96  // mu (2)
97  // (1) +->-+
98  // nu | |
99  // i+ +
100  Field_G Cup1;
101 
102  m_staple.upper(Cup1, U, mu, nu);
103 
104  // +-<-+ (2)
105  // | nu
106  // i+->-+
107  // mu (1)
108  Field_G Cup2;
109  m_staple.upper(Cup2, U, nu, mu);
110 
111  // (1) i+ +
112  // nu | |
113  // +->-+
114  // mu (2)
115  Field_G Cdn1;
116  m_staple.lower(Cdn1, U, mu, nu);
117 
118  // (2) +->-+
119  // nu |
120  // +-<-+i
121  // mu (1)
122  Field_G Cdn2;
123  m_staple.lower(Cdn2, U, nu, mu);
124 
125  Field_G Umu;
126  Umu.setpart_ex(0, U, mu);
127 
128  Field_G Unu;
129  Unu.setpart_ex(0, U, nu);
130 
131  // +->-+
132  //
133  // i+
134  Field_G Umu_nu;
135  m_shift.backward(Umu_nu, Umu, nu);
136 
137  // +
138  // |
139  // i+ +
140  Field_G Unu_mu;
141  m_shift.backward(Unu_mu, Unu, mu);
142 
143  //---------------------------------
144  //-- 2x1 part
145 
146  //- upper-right(2x1)
147  // +---<---+ + +-<-+ <---+
148  // nu | | = | + + | +
149  // i+--->---+ i+ i+ i+ >---+ i+->-+
150  Field_G c;
151  m_shift.backward(c, Cup2, mu);
152 
153  Field_G v;
154  mult_Field_Gdd(v, 0, Umu_nu, 0, Unu, 0);
155 
156  Field_G w;
157  mult_Field_Gnn(w, 0, c, 0, v, 0);
158 
159  Field_G rect;
160  mult_Field_Gnn(rect, 0, Umu, 0, w, 0);
161 
162  Fmunu_1x2 = rect;
163 
164  //- lower-right(2x1)
165  // (1) +---<---+ (1) i+---<---+
166  // nu | | -> nu | |
167  // i+--->---+ +--->---+
168  // mu (2) mu (2)
169  mult_Field_Gnd(v, 0, c, 0, Umu_nu, 0);
170  mult_Field_Gnn(w, 0, Umu, 0, v, 0);
171 
172  mult_Field_Gdn(rect, 0, Unu, 0, w, 0);
173 
174  m_shift.forward(v, rect, nu);
175 
176  axpy(Fmunu_1x2, 1.0, v);
177 
178  //- upper-left(2x1)
179  // mu (2)
180  // +---<---+ (1) +-<-+ +-<-+ +
181  // | | nu = + | + + |
182  // +---i---+ i+->-+ +->-+i +i +i +
183  //
184  // NB. shift.forward(mu)
185  // mu (2) mu (2)
186  // +---<---+ (1) +---<---+ (1)
187  // | | nu -> | | nu
188  // +---i---+ +--->---+i
189  mult_Field_Gdn(v, 0, Cdn2, 0, Umu, 0);
190  mult_Field_Gdn(w, 0, Umu_nu, 0, v, 0);
191 
192  mult_Field_Gnn(rect, 0, Unu_mu, 0, w, 0);
193 
194  m_shift.forward(v, rect, mu);
195 
196  axpy(Fmunu_1x2, 1.0, v);
197 
198  //- lower-left(2x1)
199  // mu (1)
200  // +---<---+ (2) + +-<-+ +-<-+
201  // | | nu = | + + | +
202  // +---i---+ +i + i+->-+ +->-+i +i
203  //
204  // NB. shift.forward(mu+nu)
205  // mu (1) mu (1) mu (1)
206  // +---<---+ (2) +---<---+ (2) +---<---+i (2)
207  // | | nu -> | | nu -> | | nu
208  // +---i---+ +--->---+i +--->---+
209  mult_Field_Gnn(v, 0, Umu, 0, Unu_mu, 0);
210  mult_Field_Gdn(w, 0, Cdn2, 0, v, 0);
211 
212  mult_Field_Gdn(rect, 0, Umu_nu, 0, w, 0);
213 
214  m_shift.forward(w, rect, mu);
215  m_shift.forward(v, w, nu);
216 
217  axpy(Fmunu_1x2, 1.0, v);
218  //---------------------------------
219 
220 
221  //---------------------------------
222  //-- 1x2 part
223 
224  //- upper-right(1x2)
225  // +-<-+ +-<-+
226  // | | | |
227  // (2) + + = + + + + + + +
228  // nu | | | |
229  // i+->-+ i+ i+ i+ + i+->-+
230  // mu (1)
231  m_shift.backward(c, Cup1, nu);
232 
233  mult_Field_Gdd(v, 0, c, 0, Unu, 0);
234  mult_Field_Gnn(w, 0, Unu_mu, 0, v, 0);
235 
236  mult_Field_Gnn(rect, 0, Umu, 0, w, 0);
237 
238  axpy(Fmunu_1x2, 1.0, rect);
239 
240  //- lower-right(1x2)
241  // mu (2)
242  // (1) +-<-+ +-<-+ + +
243  // nu | | | |
244  // i+ + = i+ + i+ + + i+ + + i+
245  // | | | |
246  // +->-+ +->-+
247  //
248  // NB. shift.forward(nu)
249  // mu (2) mu (2)
250  // (1) +-<-+ (1) i+-<-+
251  // nu | | nu | |
252  // i+ + -> + +
253  // | | | |
254  // +->-+ +->-+
255  mult_Field_Gnd(v, 0, Unu_mu, 0, Umu_nu, 0);
256  mult_Field_Gnn(w, 0, Cdn1, 0, v, 0);
257 
258  mult_Field_Gdn(rect, 0, Unu, 0, w, 0);
259 
260  m_shift.forward(v, rect, nu);
261 
262  axpy(Fmunu_1x2, 1.0, v);
263 
264  //- upper-left(1x2)
265  // +-<-+ +-<-+
266  // | | | |
267  // + + (1) = + + + + + + +
268  // | | nu | |
269  // i+->-+ i+->-+ i+ i+ i+ +
270  // mu (2)
271  //
272  // NB. shift.forward(mu)
273  // +-<-+ +-<-+
274  // | | | |
275  // + + (1) -> + + (1)
276  // | | nu | | nu
277  // i+->-+ +->-+i
278  // mu (2) mu (2)
279  mult_Field_Gdn(v, 0, Unu, 0, Umu, 0);
280  mult_Field_Gdn(w, 0, c, 0, v, 0);
281 
282  mult_Field_Gnn(rect, 0, Unu_mu, 0, w, 0);
283 
284  m_shift.forward(v, rect, mu);
285 
286  axpy(Fmunu_1x2, 1.0, v);
287 
288  //- lower-left(1x2)
289  // mu (1)
290  // (2) +-<-+ + + +-<-+
291  // nu | | | |
292  // i+ + = i+ + + i+ + + i+ + i+
293  // | | | |
294  // +->-+ +->-+
295  //
296  // NB. shift.forward(mu+nu)
297  // mu (1) mu (1)
298  // (2) +-<-+ (2) +-<-+i
299  // nu | | nu | |
300  // i+ + -> + +
301  // | | | |
302  // +->-+ +->-+
303  mult_Field_Gnn(v, 0, Cdn1, 0, Unu_mu, 0);
304  mult_Field_Gdn(w, 0, Unu, 0, v, 0);
305 
306  mult_Field_Gdn(rect, 0, Umu_nu, 0, w, 0);
307 
308  m_shift.forward(w, rect, mu);
309  m_shift.forward(v, w, nu);
310 
311  axpy(Fmunu_1x2, 1.0, v);
312  //---------------------------------
313 
314  ah_Field_G(Fmunu_1x2, 0);
315 
316  //- normalization = 1/8
317  // i.e. overall factor = 1/4, average of 1x2 and 2x1 = 1/2
318  scal(Fmunu_1x2, -0.125);
319  Fmunu_1x2.xI();
320 }
321 
322 
323 //====================================================================
324 void FieldStrength::construct_Fmunu_1x1_traceless(Field_G& Fmunu_1x1,
325  const int mu, const int nu, const Field_G& U)
326 {
327  const int Nvol = CommonParameters::Nvol();
328 
329  //- building blocks
330  // (1) mu (2)
331  // +->-+
332  // nu | |
333  // i+ +
334  Field_G Cup;
335 
336  m_staple.upper(Cup, U, mu, nu);
337 
338  // (1) mu (2)
339  // i+ +
340  // nu | |
341  // +->-+
342  Field_G Cdn;
343  m_staple.lower(Cdn, U, mu, nu);
344 
345  Field_G Umu;
346  Umu.setpart_ex(0, U, mu);
347 
348  //---------------------------------
349  //- right clover leaf
350  // (2) +-<-+
351  // nu | |
352  // i+->-+
353  // mu (1)
354  mult_Field_Gnd(Fmunu_1x1, 0, Umu, 0, Cup, 0);
355 
356  // (1) i+-<-+
357  // nu | |
358  // +->-+
359  // mu (2)
360  multadd_Field_Gnd(Fmunu_1x1, 0, Cdn, 0, Umu, 0, 1.0);
361  //---------------------------------
362 
363  //---------------------------------
364  //- left clover leaf
365  // mu (2)
366  // +-<-+ (1)
367  // | | nu mu (1)
368  // +->-+i + +-<-+i (2)
369  // | | nu
370  // +->-+
371  Field_G v;
372  mult_Field_Gdn(v, 0, Cup, 0, Umu, 0);
373  multadd_Field_Gdn(v, 0, Umu, 0, Cdn, 0, 1.0);
374 
375  // NB. shift.forward(mu)
376  // mu (2) mu (2)
377  // +-<-+ (1) +-<-+ (1)
378  // | | nu -> | | nu
379  // i+->-+ +->-+i
380  Field_G v2;
381  m_shift.forward(v2, v, mu);
382 
383  axpy(Fmunu_1x1, 1.0, v2);
384  //---------------------------------
385 
386  // ah_Field_G(Fmunu_1x1, 0);
387  at_Field_G(Fmunu_1x1, 0);
388 
389  scal(Fmunu_1x1, -0.25);
390  Fmunu_1x1.xI();
391 }
392 
393 
394 //====================================================================
395 void FieldStrength::construct_Fmunu_1x2_traceless(Field_G& Fmunu_1x2,
396  const int mu, const int nu, const Field_G& U)
397 {
398  const int Nvol = CommonParameters::Nvol();
399 
400  //-- building blocks
401  // mu (2)
402  // (1) +->-+
403  // nu | |
404  // i+ +
405  Field_G Cup1;
406 
407  m_staple.upper(Cup1, U, mu, nu);
408 
409  // +-<-+ (2)
410  // | nu
411  // i+->-+
412  // mu (1)
413  Field_G Cup2;
414  m_staple.upper(Cup2, U, nu, mu);
415 
416  // (1) i+ +
417  // nu | |
418  // +->-+
419  // mu (2)
420  Field_G Cdn1;
421  m_staple.lower(Cdn1, U, mu, nu);
422 
423  // (2) +->-+
424  // nu |
425  // +-<-+i
426  // mu (1)
427  Field_G Cdn2;
428  m_staple.lower(Cdn2, U, nu, mu);
429 
430  Field_G Umu;
431  Umu.setpart_ex(0, U, mu);
432 
433  Field_G Unu;
434  Unu.setpart_ex(0, U, nu);
435 
436  // +->-+
437  //
438  // i+
439  Field_G Umu_nu;
440  m_shift.backward(Umu_nu, Umu, nu);
441 
442  // +
443  // |
444  // i+ +
445  Field_G Unu_mu;
446  m_shift.backward(Unu_mu, Unu, mu);
447 
448  //---------------------------------
449  //-- 2x1 part
450 
451  //- upper-right(2x1)
452  // +---<---+ + +-<-+ <---+
453  // nu | | = | + + | +
454  // i+--->---+ i+ i+ i+ >---+ i+->-+
455  Field_G c;
456  m_shift.backward(c, Cup2, mu);
457 
458  Field_G v;
459  mult_Field_Gdd(v, 0, Umu_nu, 0, Unu, 0);
460 
461  Field_G w;
462  mult_Field_Gnn(w, 0, c, 0, v, 0);
463 
464  Field_G rect;
465  mult_Field_Gnn(rect, 0, Umu, 0, w, 0);
466 
467  Fmunu_1x2 = rect;
468 
469  //- lower-right(2x1)
470  // (1) +---<---+ (1) i+---<---+
471  // nu | | -> nu | |
472  // i+--->---+ +--->---+
473  // mu (2) mu (2)
474  mult_Field_Gnd(v, 0, c, 0, Umu_nu, 0);
475  mult_Field_Gnn(w, 0, Umu, 0, v, 0);
476 
477  mult_Field_Gdn(rect, 0, Unu, 0, w, 0);
478 
479  m_shift.forward(v, rect, nu);
480 
481  axpy(Fmunu_1x2, 1.0, v);
482 
483  //- upper-left(2x1)
484  // mu (2)
485  // +---<---+ (1) +-<-+ +-<-+ +
486  // | | nu = + | + + |
487  // +---i---+ i+->-+ +->-+i +i +i +
488  //
489  // NB. shift.forward(mu)
490  // mu (2) mu (2)
491  // +---<---+ (1) +---<---+ (1)
492  // | | nu -> | | nu
493  // +---i---+ +--->---+i
494  mult_Field_Gdn(v, 0, Cdn2, 0, Umu, 0);
495  mult_Field_Gdn(w, 0, Umu_nu, 0, v, 0);
496 
497  mult_Field_Gnn(rect, 0, Unu_mu, 0, w, 0);
498 
499  m_shift.forward(v, rect, mu);
500 
501  axpy(Fmunu_1x2, 1.0, v);
502 
503  //- lower-left(2x1)
504  // mu (1)
505  // +---<---+ (2) + +-<-+ +-<-+
506  // | | nu = | + + | +
507  // +---i---+ +i + i+->-+ +->-+i +i
508  //
509  // NB. shift.forward(mu+nu)
510  // mu (1) mu (1) mu (1)
511  // +---<---+ (2) +---<---+ (2) +---<---+i (2)
512  // | | nu -> | | nu -> | | nu
513  // +---i---+ +--->---+i +--->---+
514  mult_Field_Gnn(v, 0, Umu, 0, Unu_mu, 0);
515  mult_Field_Gdn(w, 0, Cdn2, 0, v, 0);
516 
517  mult_Field_Gdn(rect, 0, Umu_nu, 0, w, 0);
518 
519  m_shift.forward(w, rect, mu);
520  m_shift.forward(v, w, nu);
521 
522  axpy(Fmunu_1x2, 1.0, v);
523  //---------------------------------
524 
525 
526  //---------------------------------
527  //-- 1x2 part
528 
529  //- upper-right(1x2)
530  // +-<-+ +-<-+
531  // | | | |
532  // (2) + + = + + + + + + +
533  // nu | | | |
534  // i+->-+ i+ i+ i+ + i+->-+
535  // mu (1)
536  m_shift.backward(c, Cup1, nu);
537 
538  mult_Field_Gdd(v, 0, c, 0, Unu, 0);
539  mult_Field_Gnn(w, 0, Unu_mu, 0, v, 0);
540 
541  mult_Field_Gnn(rect, 0, Umu, 0, w, 0);
542 
543  axpy(Fmunu_1x2, 1.0, rect);
544 
545  //- lower-right(1x2)
546  // mu (2)
547  // (1) +-<-+ +-<-+ + +
548  // nu | | | |
549  // i+ + = i+ + i+ + + i+ + + i+
550  // | | | |
551  // +->-+ +->-+
552  //
553  // NB. shift.forward(nu)
554  // mu (2) mu (2)
555  // (1) +-<-+ (1) i+-<-+
556  // nu | | nu | |
557  // i+ + -> + +
558  // | | | |
559  // +->-+ +->-+
560  mult_Field_Gnd(v, 0, Unu_mu, 0, Umu_nu, 0);
561  mult_Field_Gnn(w, 0, Cdn1, 0, v, 0);
562 
563  mult_Field_Gdn(rect, 0, Unu, 0, w, 0);
564 
565  m_shift.forward(v, rect, nu);
566 
567  axpy(Fmunu_1x2, 1.0, v);
568 
569  //- upper-left(1x2)
570  // +-<-+ +-<-+
571  // | | | |
572  // + + (1) = + + + + + + +
573  // | | nu | |
574  // i+->-+ i+->-+ i+ i+ i+ +
575  // mu (2)
576  //
577  // NB. shift.forward(mu)
578  // +-<-+ +-<-+
579  // | | | |
580  // + + (1) -> + + (1)
581  // | | nu | | nu
582  // i+->-+ +->-+i
583  // mu (2) mu (2)
584  mult_Field_Gdn(v, 0, Unu, 0, Umu, 0);
585  mult_Field_Gdn(w, 0, c, 0, v, 0);
586 
587  mult_Field_Gnn(rect, 0, Unu_mu, 0, w, 0);
588 
589  m_shift.forward(v, rect, mu);
590 
591  axpy(Fmunu_1x2, 1.0, v);
592 
593  //- lower-left(1x2)
594  // mu (1)
595  // (2) +-<-+ + + +-<-+
596  // nu | | | |
597  // i+ + = i+ + + i+ + + i+ + i+
598  // | | | |
599  // +->-+ +->-+
600  //
601  // NB. shift.forward(mu+nu)
602  // mu (1) mu (1)
603  // (2) +-<-+ (2) +-<-+i
604  // nu | | nu | |
605  // i+ + -> + +
606  // | | | |
607  // +->-+ +->-+
608  mult_Field_Gnn(v, 0, Cdn1, 0, Unu_mu, 0);
609  mult_Field_Gdn(w, 0, Unu, 0, v, 0);
610 
611  mult_Field_Gdn(rect, 0, Umu_nu, 0, w, 0);
612 
613  m_shift.forward(w, rect, mu);
614  m_shift.forward(v, w, nu);
615 
616  axpy(Fmunu_1x2, 1.0, v);
617  //---------------------------------
618 
619  // ah_Field_G(Fmunu_1x2, 0);
620  at_Field_G(Fmunu_1x2, 0);
621 
622  //- normalization = 1/8
623  // i.e. overall factor = 1/4, average of 1x2 and 2x1 = 1/2
624  scal(Fmunu_1x2, -0.125);
625  Fmunu_1x2.xI();
626 }
627 
628 
629 //====================================================================
630 void FieldStrength::construct_Fmunu_plaq_traceless(Field_G& Fmunu_plaq,
631  const int mu, const int nu, const Field_G& U)
632 {
633  const int Nvol = CommonParameters::Nvol();
634 
635  //- building blocks
636  // (1) mu (2)
637  // +->-+
638  // nu | |
639  // i+ +
640  Field_G Cup;
641 
642  m_staple.upper(Cup, U, mu, nu);
643 
644  Field_G Umu;
645  Umu.setpart_ex(0, U, mu);
646 
647  //- right clover leaf
648  // (2) +-<-+
649  // nu | |
650  // i+->-+
651  // mu (1)
652  mult_Field_Gnd(Fmunu_plaq, 0, Umu, 0, Cup, 0);
653 
654  // ah_Field_G(Fmunu_plaq, 0);
655  at_Field_G(Fmunu_plaq, 0);
656 
657  Fmunu_plaq.xI();
658 }
659 
660 
661 //====================================================================
662 double FieldStrength::contract(const Field_G& Fmunu_1, const Field_G& Fmunu_2)
663 {
664  const int Nvol = CommonParameters::Nvol();
665 
666  double FF = 0.0;
667 
668  for (int site = 0; site < Nvol; ++site) {
669  FF += ReTr(Fmunu_1.mat(site) * Fmunu_2.mat(site));
670  }
671  const double result = Communicator::reduce_sum(FF);
672 
673  return result;
674 }
675 
676 
677 //====================================================================
678 void FieldStrength::contract_at_t(std::vector<double>& corr_global,
679  const Field_G& Fmunu_1, const Field_G& Fmunu_2)
680 {
681  const int Nx = CommonParameters::Nx();
682  const int Ny = CommonParameters::Ny();
683  const int Nz = CommonParameters::Nz();
684  const int Nt = CommonParameters::Nt();
685  const int Lt = CommonParameters::Lt();
686 
687  assert(corr_global.size() == Lt);
688 
689  std::vector<double> corr_local(Nt, 0.0);
690  for (int t = 0; t < Nt; ++t) {
691  for (int z = 0; z < Nz; ++z) {
692  for (int y = 0; y < Ny; ++y) {
693  for (int x = 0; x < Nx; ++x) {
694  int site = x + Nx * (y + Ny * (z + Nz * t));
695 
696  corr_local[t] += ReTr(Fmunu_1.mat(site) * Fmunu_2.mat(site));
697  }
698  }
699  }
700  }
701  global_corr_t(corr_global, corr_local);
702 }
703 
704 
705 //====================================================================
706 void FieldStrength::contract_at_t(std::vector<double>& corr_global,
707  const Field_G& Fmunu_1,
708  const Field_G& Fmunu_2,
709  const std::vector<int>& momentum_sink,
710  const std::vector<int>& source_position)
711 {
712  int Nx = CommonParameters::Nx();
713  int Ny = CommonParameters::Ny();
714  int Nz = CommonParameters::Nz();
715  int Nt = CommonParameters::Nt();
716  const int Lx = CommonParameters::Lx();
717  const int Ly = CommonParameters::Ly();
718  const int Lz = CommonParameters::Lz();
719  const int Lt = CommonParameters::Lt();
720  const int ND = CommonParameters::Nd();
721 
722  assert(corr_global.size() == Lt);
723  assert(momentum_sink.size() == ND - 1);
724  assert(source_position.size() == ND);
725 
726  static const double PI = 4.0 * atan(1.0);
727  std::vector<double> p_unit(ND - 1);
728  p_unit[0] = (2.0 * PI / Lx) * momentum_sink[0];
729  p_unit[1] = (2.0 * PI / Ly) * momentum_sink[1];
730  p_unit[2] = (2.0 * PI / Lz) * momentum_sink[2];
731 
732  std::vector<int> ipe(ND - 1);
733  ipe[0] = Communicator::ipe(0);
734  ipe[1] = Communicator::ipe(1);
735  ipe[2] = Communicator::ipe(2);
736 
737  std::vector<double> corr_local(Nt, 0);
738  for (int t = 0; t < Nt; ++t) {
739  for (int z = 0; z < Nz; ++z) {
740  for (int y = 0; y < Ny; ++y) {
741  for (int x = 0; x < Nx; ++x) {
742  int site = x + Nx * (y + Ny * (z + Nz * t));
743  int x_global = x + ipe[0] * Nx;
744  int y_global = y + ipe[1] * Ny;
745  int z_global = z + ipe[2] * Nz;
746  double p_x = p_unit[0] * (x_global - source_position[0]);
747  double p_y = p_unit[1] * (y_global - source_position[1]);
748  double p_z = p_unit[2] * (z_global - source_position[2]);
749  double cos_p_xyz = cos(p_x + p_y + p_z);
750 
751  double FF = ReTr(Fmunu_1.mat(site) * Fmunu_2.mat(site));
752  corr_local[t] += FF * cos_p_xyz;
753  }
754  }
755  }
756  }
757  global_corr_t(corr_global, corr_local);
758 }
759 
760 
761 //====================================================================
762 void FieldStrength::contract_at_x(std::vector<double>& corr_global,
763  const Field_G& Fmunu_1,
764  const Field_G& Fmunu_2)
765 {
766  const int Nx = CommonParameters::Nx();
767  const int Ny = CommonParameters::Ny();
768  const int Nz = CommonParameters::Nz();
769  const int Nt = CommonParameters::Nt();
770  const int Lx = CommonParameters::Lx();
771 
772  assert(corr_global.size() == Lx);
773 
774  std::vector<double> corr_local(Nx, 0.0);
775  for (int x = 0; x < Nx; ++x) {
776  for (int t = 0; t < Nt; ++t) {
777  for (int z = 0; z < Nz; ++z) {
778  for (int y = 0; y < Ny; ++y) {
779  int site = x + Nx * (y + Ny * (z + Nz * t));
780 
781  corr_local[x] += ReTr(Fmunu_1.mat(site) * Fmunu_2.mat(site));
782  }
783  }
784  }
785  }
786  global_corr_x(corr_global, corr_local);
787 }
788 
789 
790 //====================================================================
791 void FieldStrength::contract_at_x(std::vector<double>& corr_global,
792  const Field_G& Fmunu_1,
793  const Field_G& Fmunu_2,
794  const std::vector<int>& momentum_sink,
795  const std::vector<int>& source_position)
796 {
797  const int Nx = CommonParameters::Nx();
798  const int Ny = CommonParameters::Ny();
799  const int Nz = CommonParameters::Nz();
800  const int Nt = CommonParameters::Nt();
801  const int Lx = CommonParameters::Lx();
802  const int Ly = CommonParameters::Ly();
803  const int Lz = CommonParameters::Lz();
804  const int Lt = CommonParameters::Lt();
805  const int ND = CommonParameters::Nd();
806 
807  assert(corr_global.size() == Lx);
808  assert(momentum_sink.size() == ND - 1);
809  assert(source_position.size() == ND);
810 
811  static const double PI = 4.0 * atan(1.0);
812  std::vector<double> p_unit(ND - 1);
813  p_unit[0] = (2.0 * PI / Ly) * momentum_sink[0];
814  p_unit[1] = (2.0 * PI / Lz) * momentum_sink[1];
815  p_unit[2] = (2.0 * PI / Lt) * momentum_sink[2];
816 
817  std::vector<int> ipe(ND);
818  ipe[0] = Communicator::ipe(0);
819  ipe[1] = Communicator::ipe(1);
820  ipe[2] = Communicator::ipe(2);
821  ipe[3] = Communicator::ipe(3);
822 
823  std::vector<double> corr_local(Nx, 0);
824  for (int x = 0; x < Nx; ++x) {
825  for (int t = 0; t < Nt; ++t) {
826  for (int z = 0; z < Nz; ++z) {
827  for (int y = 0; y < Ny; ++y) {
828  int site = x + Nx * (y + Ny * (z + Nz * t));
829  int y_global = y + ipe[1] * Ny;
830  int z_global = z + ipe[2] * Nz;
831  int t_global = t + ipe[3] * Nt;
832  double p_y = p_unit[0] * (y_global - source_position[1]);
833  double p_z = p_unit[1] * (z_global - source_position[2]);
834  double p_t = p_unit[2] * (t_global - source_position[3]);
835  double cos_p_yzt = cos(p_t + p_y + p_z);
836  double FF = ReTr(Fmunu_1.mat(site) * Fmunu_2.mat(site));
837  corr_local[x] += FF * cos_p_yzt;
838  }
839  }
840  }
841  }
842  global_corr_x(corr_global, corr_local);
843 }
844 
845 
846 //====================================================================
847 void FieldStrength::contract_at_y(std::vector<double>& corr_global,
848  const Field_G& Fmunu_1,
849  const Field_G& Fmunu_2)
850 {
851  const int Nx = CommonParameters::Nx();
852  const int Ny = CommonParameters::Ny();
853  const int Nz = CommonParameters::Nz();
854  const int Nt = CommonParameters::Nt();
855  const int Ly = CommonParameters::Ly();
856 
857  assert(corr_global.size() == Ly);
858 
859  std::vector<double> corr_local(Ny, 0.0);
860  for (int y = 0; y < Ny; ++y) {
861  for (int t = 0; t < Nt; ++t) {
862  for (int z = 0; z < Nz; ++z) {
863  for (int x = 0; x < Nx; ++x) {
864  int site = x + Nx * (y + Ny * (z + Nz * t));
865 
866  corr_local[y] += ReTr(Fmunu_1.mat(site) * Fmunu_2.mat(site));
867  }
868  }
869  }
870  }
871  global_corr_y(corr_global, corr_local);
872 }
873 
874 
875 //====================================================================
876 void FieldStrength::contract_at_y(std::vector<double>& corr_global,
877  const Field_G& Fmunu_1,
878  const Field_G& Fmunu_2,
879  const std::vector<int>& momentum_sink,
880  const std::vector<int>& source_position)
881 {
882  const int Nx = CommonParameters::Nx();
883  const int Ny = CommonParameters::Ny();
884  const int Nz = CommonParameters::Nz();
885  const int Nt = CommonParameters::Nt();
886  const int Lx = CommonParameters::Lx();
887  const int Ly = CommonParameters::Ly();
888  const int Lz = CommonParameters::Lz();
889  const int Lt = CommonParameters::Lt();
890  const int ND = CommonParameters::Nd();
891 
892  assert(corr_global.size() == Ly);
893  assert(momentum_sink.size() == ND - 1);
894  assert(source_position.size() == ND);
895 
896  static const double PI = 4.0 * atan(1.0);
897  std::vector<double> p_unit(ND - 1);
898  p_unit[0] = (2.0 * PI / Ly) * momentum_sink[0];
899  p_unit[1] = (2.0 * PI / Lz) * momentum_sink[1];
900  p_unit[2] = (2.0 * PI / Lt) * momentum_sink[2];
901 
902  std::vector<int> ipe(ND);
903  ipe[0] = Communicator::ipe(0);
904  ipe[1] = Communicator::ipe(1);
905  ipe[2] = Communicator::ipe(2);
906  ipe[3] = Communicator::ipe(3);
907 
908  std::vector<double> corr_local(Ny, 0);
909  for (int y = 0; y < Ny; ++y) {
910  for (int t = 0; t < Nt; ++t) {
911  for (int z = 0; z < Nz; ++z) {
912  for (int x = 0; x < Nx; ++x) {
913  int site = x + Nx * (y + Ny * (z + Nz * t));
914  int x_global = x + ipe[0] * Nx;
915  int z_global = z + ipe[2] * Nz;
916  int t_global = t + ipe[3] * Nt;
917  double p_x = p_unit[0] * (x_global - source_position[0]);
918  double p_z = p_unit[1] * (z_global - source_position[2]);
919  double p_t = p_unit[2] * (t_global - source_position[3]);
920  double cos_p_xzt = cos(p_t + p_x + p_z);
921  double FF = ReTr(Fmunu_1.mat(site) * Fmunu_2.mat(site));
922  corr_local[y] += FF * cos_p_xzt;
923  }
924  }
925  }
926  }
927  global_corr_y(corr_global, corr_local);
928 }
929 
930 
931 //====================================================================
932 void FieldStrength::contract_at_z(std::vector<double>& corr_global,
933  const Field_G& Fmunu_1,
934  const Field_G& Fmunu_2)
935 {
936  const int Nx = CommonParameters::Nx();
937  const int Ny = CommonParameters::Ny();
938  const int Nz = CommonParameters::Nz();
939  const int Nt = CommonParameters::Nt();
940  const int Lz = CommonParameters::Lz();
941 
942  assert(corr_global.size() == Lz);
943 
944  std::vector<double> corr_local(Nz, 0.0);
945  for (int z = 0; z < Nz; ++z) {
946  for (int t = 0; t < Nt; ++t) {
947  for (int y = 0; y < Ny; ++y) {
948  for (int x = 0; x < Nx; ++x) {
949  int site = x + Nx * (y + Ny * (z + Nz * t));
950 
951  corr_local[z] += ReTr(Fmunu_1.mat(site) * Fmunu_2.mat(site));
952  }
953  }
954  }
955  }
956  global_corr_z(corr_global, corr_local);
957 }
958 
959 
960 //====================================================================
961 void FieldStrength::contract_at_z(std::vector<double>& corr_global,
962  const Field_G& Fmunu_1,
963  const Field_G& Fmunu_2,
964  const std::vector<int>& momentum_sink,
965  const std::vector<int>& source_position)
966 {
967  const int Nx = CommonParameters::Nx();
968  const int Ny = CommonParameters::Ny();
969  const int Nz = CommonParameters::Nz();
970  const int Nt = CommonParameters::Nt();
971  const int Lx = CommonParameters::Lx();
972  const int Ly = CommonParameters::Ly();
973  const int Lz = CommonParameters::Lz();
974  const int Lt = CommonParameters::Lt();
975  const int ND = CommonParameters::Nd();
976 
977  assert(corr_global.size() == Lz);
978  assert(momentum_sink.size() == ND - 1);
979  assert(source_position.size() == ND);
980 
981  static const double PI = 4.0 * atan(1.0);
982  std::vector<double> p_unit(ND - 1);
983  p_unit[0] = (2.0 * PI / Ly) * momentum_sink[0];
984  p_unit[1] = (2.0 * PI / Lz) * momentum_sink[1];
985  p_unit[2] = (2.0 * PI / Lt) * momentum_sink[2];
986 
987  std::vector<int> ipe(ND);
988  ipe[0] = Communicator::ipe(0);
989  ipe[1] = Communicator::ipe(1);
990  ipe[2] = Communicator::ipe(2);
991  ipe[3] = Communicator::ipe(3);
992 
993  std::vector<double> corr_local(Nz, 0);
994  for (int z = 0; z < Nz; ++z) {
995  for (int t = 0; t < Nt; ++t) {
996  for (int y = 0; y < Ny; ++y) {
997  for (int x = 0; x < Nx; ++x) {
998  int site = x + Nx * (y + Ny * (z + Nz * t));
999  int x_global = x + ipe[0] * Nx;
1000  int y_global = y + ipe[1] * Ny;
1001  int t_global = t + ipe[3] * Nt;
1002  double p_x = p_unit[0] * (x_global - source_position[0]);
1003  double p_y = p_unit[1] * (y_global - source_position[1]);
1004  double p_t = p_unit[2] * (t_global - source_position[3]);
1005  double cos_p_xyt = cos(p_t + p_x + p_y);
1006  double FF = ReTr(Fmunu_1.mat(site) * Fmunu_2.mat(site));
1007  corr_local[z] += FF * cos_p_xyt;
1008  }
1009  }
1010  }
1011  }
1012  global_corr_z(corr_global, corr_local);
1013 }
1014 
1015 
1016 //====================================================================
1017 void FieldStrength::global_corr_x(std::vector<double>& corr_global,
1018  const std::vector<double>& corr_local)
1019 {
1020  const int Lx = CommonParameters::Lx();
1021  const int Nx = CommonParameters::Nx();
1022 
1023  assert(corr_global.size() == Lx);
1024  assert(corr_local.size() == Nx);
1025 
1026  const int ipe_x = Communicator::ipe(0);
1027 
1028  std::vector<double> corr_tmp(Lx, 0.0);
1029 
1030  for (int x = 0; x < Nx; ++x) {
1031  int x_global = x + ipe_x * Nx;
1032  corr_tmp[x_global] = corr_local[x];
1033  }
1034 
1035  for (int x_global = 0; x_global < Lx; ++x_global) {
1036  double cr_r = Communicator::reduce_sum(corr_tmp[x_global]);
1037  corr_global[x_global] = cr_r;
1038  }
1039 }
1040 
1041 
1042 //====================================================================
1043 void FieldStrength::global_corr_y(std::vector<double>& corr_global,
1044  const std::vector<double>& corr_local)
1045 {
1046  const int Ly = CommonParameters::Ly();
1047  const int Ny = CommonParameters::Ny();
1048 
1049  assert(corr_global.size() == Ly);
1050  assert(corr_local.size() == Ny);
1051 
1052  const int ipe_y = Communicator::ipe(1);
1053 
1054  std::vector<double> corr_tmp(Ly, 0.0);
1055 
1056  for (int y = 0; y < Ny; ++y) {
1057  int y_global = y + ipe_y * Ny;
1058  corr_tmp[y_global] = corr_local[y];
1059  }
1060 
1061  for (int y_global = 0; y_global < Ly; ++y_global) {
1062  double cr_r = Communicator::reduce_sum(corr_tmp[y_global]);
1063  corr_global[y_global] = cr_r;
1064  }
1065 }
1066 
1067 
1068 //====================================================================
1069 void FieldStrength::global_corr_z(std::vector<double>& corr_global,
1070  const std::vector<double>& corr_local)
1071 {
1072  const int Lz = CommonParameters::Lz();
1073  const int Nz = CommonParameters::Nz();
1074 
1075  assert(corr_global.size() == Lz);
1076  assert(corr_local.size() == Nz);
1077 
1078  const int ipe_z = Communicator::ipe(2);
1079 
1080  std::vector<double> corr_tmp(Lz, 0.0);
1081 
1082  for (int z = 0; z < Nz; ++z) {
1083  int z_global = z + ipe_z * Nz;
1084  corr_tmp[z_global] = corr_local[z];
1085  }
1086 
1087  for (int z_global = 0; z_global < Lz; ++z_global) {
1088  double cr_r = Communicator::reduce_sum(corr_tmp[z_global]);
1089  corr_global[z_global] = cr_r;
1090  }
1091 }
1092 
1093 
1094 //====================================================================
1095 void FieldStrength::global_corr_t(std::vector<double>& corr_global,
1096  const std::vector<double>& corr_local)
1097 {
1098  const int Lt = CommonParameters::Lt();
1099  const int Nt = CommonParameters::Nt();
1100 
1101  assert(corr_global.size() == Lt);
1102  assert(corr_local.size() == Nt);
1103 
1104  const int ipe_t = Communicator::ipe(3);
1105 
1106  std::vector<double> corr_tmp(Lt, 0.0);
1107 
1108  for (int t = 0; t < Nt; ++t) {
1109  int t_global = t + ipe_t * Nt;
1110  corr_tmp[t_global] = corr_local[t];
1111  }
1112 
1113  for (int t_global = 0; t_global < Lt; ++t_global) {
1114  double cr_r = Communicator::reduce_sum(corr_tmp[t_global]);
1115  corr_global[t_global] = cr_r;
1116  }
1117 }
1118 
1119 
1120 //====================================================================
1121 //============================================================END=====
FieldStrength::contract_at_y
void contract_at_y(std::vector< double > &corr_global, const Field_G &Fmunu_1, const Field_G &Fmunu_2)
Calculate corr[y]= . Intended to be used to calculate correlations function of the topological charge...
Definition: fieldStrength.cpp:847
scal
void scal(Field &x, const double a)
scal(x, a): x = a * x
Definition: field.cpp:433
FieldStrength::m_staple
Staple_lex m_staple
Definition: fieldStrength.h:52
FieldStrength::construct_Fmunu_1x2
void construct_Fmunu_1x2(Field_G &Fmunu, const int mu, const int nu, const Field_G &U)
Constructs the anti-Hermitian field strength with eight 1x2 rectangular clover leaves.
Definition: fieldStrength.cpp:89
FieldStrength::m_shift
ShiftField_lex m_shift
Definition: fieldStrength.h:51
ah_Field_G
void ah_Field_G(Field_G &W, const int ex)
Definition: field_G_imp.cpp:501
FieldStrength::contract_at_z
void contract_at_z(std::vector< double > &corr_global, const Field_G &Fmunu_1, const Field_G &Fmunu_2)
Calculate corr[z]= . Intended to be used to calculate correlations function of the topological charge...
Definition: fieldStrength.cpp:932
multadd_Field_Gnd
void multadd_Field_Gnd(Field_G &W, const int ex, const Field_G &U1, const int ex1, const Field_G &U2, const int ex2, const double ff)
Definition: field_G_imp.cpp:365
ND
#define ND
Definition: field_F_imp_SU2-inc.h:5
at_Field_G
void at_Field_G(Field_G &W, const int ex)
Definition: field_G_imp.cpp:455
mult_Field_Gdn
void mult_Field_Gdn(Field_G &W, const int ex, const Field_G &U1, const int ex1, const Field_G &U2, const int ex2)
Definition: field_G_imp.cpp:149
FieldStrength::class_name
static const std::string class_name
Definition: fieldStrength.h:45
FieldStrength::global_corr_t
void global_corr_t(std::vector< double > &corr_global, const std::vector< double > &corr_local)
transform node-local correlator in t into global.
Definition: fieldStrength.cpp:1095
Communicator::ipe
static int ipe(const int dir)
logical coordinate of current proc.
Definition: communicator.cpp:105
FieldStrength::construct_Fmunu_1x1_traceless
void construct_Fmunu_1x1_traceless(Field_G &Fmunu, const int mu, const int nu, const Field_G &U)
Constructs the anti-Hermitian traceless field strength with four 1x1 plquette clover leaves...
Definition: fieldStrength.cpp:324
Staple_lex::lower
void lower(Field_G &, const Field_G &, const int mu, const int nu)
constructs lower staple in mu-nu plane.
Definition: staple_lex.cpp:177
Field_G
SU(N) gauge field.
Definition: field_G.h:38
FieldStrength::contract_at_t
void contract_at_t(std::vector< double > &corr_global, const Field_G &Fmunu_1, const Field_G &Fmunu_2)
Calculate corr[t]= . Intended to be used to calculate correlations function of the topological charge...
Definition: fieldStrength.cpp:678
FieldStrength::global_corr_y
void global_corr_y(std::vector< double > &corr_global, const std::vector< double > &corr_local)
transform node-local correlator in y into global.
Definition: fieldStrength.cpp:1043
FieldStrength::construct_Fmunu_1x1
void construct_Fmunu_1x1(Field_G &Fmunu, const int mu, const int nu, const Field_G &U)
Constructs the anti-Hermitian field strength with four 1x1 plquette clover leaves.
Definition: fieldStrength.cpp:19
mult_Field_Gdd
void mult_Field_Gdd(Field_G &W, const int ex, const Field_G &U1, const int ex1, const Field_G &U2, const int ex2)
Definition: field_G_imp.cpp:233
ShiftField_lex::backward
void backward(Field &, const Field &, const int mu)
Definition: shiftField_lex_imp.cpp:20
FieldStrength::construct_Fmunu_plaq_traceless
void construct_Fmunu_plaq_traceless(Field_G &Fmunu, const int mu, const int nu, const Field_G &U)
Constructs the anti-Hermitian traceless field strength with an imaginary part of the plaquette...
Definition: fieldStrength.cpp:630
Staple_lex::upper
void upper(Field_G &, const Field_G &, const int mu, const int nu)
constructs upper staple in mu-nu plane.
Definition: staple_lex.cpp:151
axpy
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
Definition: field.cpp:319
mult_Field_Gnn
void mult_Field_Gnn(Field_G &W, const int ex, const Field_G &U1, const int ex1, const Field_G &U2, const int ex2)
Definition: field_G_imp.cpp:107
Communicator::reduce_sum
static int reduce_sum(int count, double *recv_buf, double *send_buf, int pattern=0)
make a global sum of an array of double over the communicator. pattern specifies the dimensions to be...
Definition: communicator.cpp:244
FieldStrength::global_corr_z
void global_corr_z(std::vector< double > &corr_global, const std::vector< double > &corr_local)
transform node-local correlator in z into global.
Definition: fieldStrength.cpp:1069
FieldStrength::construct_Fmunu_1x2_traceless
void construct_Fmunu_1x2_traceless(Field_G &Fmunu, const int mu, const int nu, const Field_G &U)
Constructs the anti-Hermitian traceless field strength with eight 1x2 rectangular clover leaves...
Definition: fieldStrength.cpp:395
FieldStrength::global_corr_x
void global_corr_x(std::vector< double > &corr_global, const std::vector< double > &corr_local)
transform node-local correlator in x into global.
Definition: fieldStrength.cpp:1017
Field::setpart_ex
void setpart_ex(int ex, const Field &w, int exw)
Definition: field.h:197
FieldStrength::contract
double contract(const Field_G &Fmunu_1, const Field_G &Fmunu_2)
Calculate and returns its value. Intended to be used for the topological charge and energy momentum ...
Definition: fieldStrength.cpp:662
FieldStrength::contract_at_x
void contract_at_x(std::vector< double > &corr_global, const Field_G &Fmunu_1, const Field_G &Fmunu_2)
Calculate corr[x]= . Intended to be used to calculate correlations function of the topological charge...
Definition: fieldStrength.cpp:762
Field_G::mat
Mat_SU_N mat(const int site, const int mn=0) const
Definition: field_G.h:114
SU_N::ReTr
double ReTr(const Mat_SU_N &m)
Definition: mat_SU_N.h:488
multadd_Field_Gdn
void multadd_Field_Gdn(Field_G &W, const int ex, const Field_G &U1, const int ex1, const Field_G &U2, const int ex2, const double ff)
Definition: field_G_imp.cpp:320
ShiftField_lex::forward
void forward(Field &, const Field &, const int mu)
Definition: shiftField_lex_imp.cpp:41
mult_Field_Gnd
void mult_Field_Gnd(Field_G &W, const int ex, const Field_G &U1, const int ex1, const Field_G &U2, const int ex2)
Definition: field_G_imp.cpp:191
Field_G::xI
void xI()
Definition: field_G.h:184