CHROMA
wilslp.cc
Go to the documentation of this file.
1 /*! \file
2  * \brief Calculate Wilson loops
3  */
4 
5 #include "chromabase.h"
6 #include "meas/glue/wilslp.h"
7 #include "meas/gfix/axgauge.h"
8 
9 namespace Chroma
10 {
11 
12  //! Calculate Wilson loops
13  /*!
14  * \ingroup glue
15  *
16  * Calculates, depending on option, (1) "space-like" planar
17  * Wilson loops in the directions perpendicular to j_decay
18  * that have equal length, (2) "time-like" planar Wilson loops
19  * with time direction j_decay and space directions the
20  * perpendicular ones that have equal length and (3) off-axis
21  * "time-like" Wilson loops along 3 paricular paths in the
22  * space directions that have equal length.
23  *
24  * \param u gauge field (Read)
25  * \param j_decay decay direction (Read)
26  * \param t_dir time direction (Read)
27  * \param kind binary-combined YES/NO [1/0] of the three options (Read)
28  * e.g. kind = 2 gives planar t-like, kind=6 is
29  * planar + off-axis: sqrt(2), sqrt(5), sqrt(3)
30  */
31 
32  void wilslp(const multi1d<LatticeColorMatrix>& u,
33  int j_decay, int t_dir, int kind,
34  XMLWriter& xml, const std::string& xml_group)
35  {
36  START_CODE();
37 
38  multi1d<int> nrow(Nd);
39  nrow = Layout::lattSize();
40 
41  multi1d<LatticeColorMatrix> ug(Nd);
42  multi1d<int> space_dir(Nd);
43  LatticeColorMatrix u_t;
44  LatticeColorMatrix u_tmp;
45  LatticeColorMatrix u_diag;
46  LatticeColorMatrix u_space;
47  LatticeColorMatrix tmp_2;
48  LatticeColorMatrix tmp_3;
49  LatticeReal wl_trace;
50  LatticeInteger t_coord;
51  LatticeBoolean btmp;
52  Double dummy;
53  Real ftmp;
54 
55  int lengthr;
56  int lengtht;
57  int lsizet;
58  int lsizer;
59  int nspace;
60  int mu;
61  int nu;
62  int rho;
63  int r;
64  int t;
65  int tt;
66  int nr;
67  int i;
68  int j;
69  int r_off;
70 
71  lsizet = nrow[j_decay];
72  lengtht = lsizet / 2;
73 
74  /* Determine 'lattice shape' */
75  space_dir = 1;
76  nr = 1;
77  for(mu = 0;mu <= ( Nd-2); ++mu )
78  if ( mu != j_decay && mu != t_dir )
79  {
80  lsizer = nrow[mu];
81  for(nu = mu+1;nu < ( Nd); ++nu )
82  if ( nu != j_decay && nu != t_dir )
83  {
84  r = nrow[nu];
85  if( r == lsizer )
86  space_dir[mu] += nr;
87  }
88  }
89 
90  nspace = space_dir[0];
91  nu = 0;
92  for(mu = 1;mu <= ( Nd-2); ++mu )
93  if ( mu != j_decay && mu != t_dir )
94  {
95  nr = space_dir[mu];
96  if( nr > nspace )
97  {
98  nspace = nr;
99  nu = mu;
100  }
101  }
102 
103  lsizer = nrow[nu];
104  nr = 0;
105  space_dir[nr] = nu;
106  for(mu = nu+1;mu < ( Nd); ++mu )
107  if ( mu != j_decay && mu != t_dir )
108  {
109  r = nrow[mu];
110  if( r == lsizer )
111  {
112  nr = nr + 1;
113  space_dir[nr] = mu;
114  }
115  }
116 
117  QDPIO::cout << "nspace = " << nspace
118  << " space_dir = "
119  << space_dir[0] << " "
120  << space_dir[1] << " "
121  << space_dir[2] << " "
122  << space_dir[3] << std::endl;
123 
124  if( nr != (nspace-1) )
125  {
126  QDPIO::cerr << __func__
127  << ": Trouble with space dimensions for Wilson loops: nr = " << nr
128  << " nspace = " << nspace
129  << " nu = " << nu << std::endl;
130  QDP_abort(1);
131  }
132 
133 // lengthr = lsizer/2; // Old style with measurement on L/2
134  lengthr = lsizer; // Now use all L
135 
136  /* Compute "space-like" planar Wilson loops, if desired */
137  if ( (kind & 1) != 0 )
138  {
139  if( nspace < 2 )
140  for(mu = 0;mu < ( Nd); ++mu )
141  QDP_error_exit("Wrong lattice size for space-like Wilson loops: ",
142  mu, nrow[mu]);
143 
144  multi2d<Double> wils_loop1(lengthr, lengthr);
145  wils_loop1 = 0; /* initialize the Wilson loops */
146 
147  for(j = 1;j < ( nspace); ++j )
148  {
149  nu = space_dir[j];
150 
151  /*+ */
152  /* Fix to axial gauge in nu-direction */
153  /*- */
154  ug = u;
155  axGauge (ug, nu);
156  t_coord = Layout::latticeCoordinate(nu);
157 
158  for(i = 0;i < ( j); ++i )
159  {
160  mu = space_dir[i];
161 
162  for(r = 0;r < ( lengthr); ++r )
163  {
164  /* Gather t-links (i.e. nu-links) from r-direction (i.e. mu) */
165  /* and make "r-links" */
166 
167  if ( r == 0 )
168  {
169  u_t = shift(ug[nu], FORWARD, mu);
170  u_space = ug[mu];
171  }
172  else
173  {
174  tmp_2 = shift(u_t, FORWARD, mu);
175  u_t = tmp_2;
176 
177  tmp_2 = shift(u_space, FORWARD, mu);
178  u_space = ug[mu] * tmp_2;
179  }
180 
181  for(t = 0;t < ( lengthr); ++t )
182  {
183  /* Gather r-link from t-direction (i.e. nu) */
184 
185  if ( t == 0 )
186  {
187  u_tmp = shift(u_space, FORWARD, nu);
188  }
189  else
190  {
191  tmp_2 = shift(u_tmp, FORWARD, nu);
192  u_tmp = tmp_2;
193  }
194 
195  tt = lsizer - t - 1;
196  btmp = t_coord < tt;
197 
198  tmp_2 = ug[nu] * u_tmp;
199  tmp_3 = tmp_2 * adj(u_t);
200  copymask(tmp_3, btmp, u_tmp);
201  tmp_2 = tmp_3 * adj(u_space);
202  wl_trace = real(trace(tmp_2));
203  wils_loop1[r][t] += sum(wl_trace);
204  } /* end t loop */
205  } /* end r loop */
206  } /* end i loop (for mu) */
207  } /* end j loop (for nu) */
208 
209  dummy = 2.0 / double (Layout::vol()*Nc*nspace*(nspace-1)) ;
210 
211  push(xml, "wils_loop1"); // XML tag for wils_wloop1
212  write(xml, "lengthr", lengthr);
213  push(xml, "wloop1");
214 
215  multi1d<Double> wloop1(lengthr);
216 
217  for(r = 0; r < lengthr; ++r)
218  {
219  for(t = 0; t < lengthr; ++t)
220  {
221  wils_loop1[t][r] = wils_loop1[t][r] * dummy;
222  wloop1[t] = wils_loop1[t][r];
223  }
224  push(xml, "elem");
225  write(xml, "r", r);
226  write(xml, "loop", wloop1); // write out wils_wloop1
227  pop(xml); // elem
228  } // end for r
229 
230  pop(xml); // XML end tag for wloop1
231  pop(xml); // XML end tag for wils_wloop1
232  QDPIO::cout << "wils_loop1 data written to .xml file " << std::endl;
233 
234  } /* end of option "space-like planar Wilson loops" */
235 
236  /* Fix to axial gauge */
237  ug = u;
238  axGauge (ug, j_decay);
239  t_coord = Layout::latticeCoordinate(j_decay);
240 
241  /* Compute "time-like" planar Wilson loops, if desired */
242  if ( (kind & 2) != 0 )
243  {
244  multi2d<Double> wils_loop2(lengtht, lengthr);
245  wils_loop2 = 0; /* initialize the Wilson loops */
246 
247  QDPIO::cout << "computing time-like Wilson loops" << std::endl;
248 
249  for(i = 0;i < ( nspace); ++i )
250  {
251  mu = space_dir[i];
252  for(r = 0;r < ( lengthr); ++r )
253  {
254  /* Gather t-links from space-direction (i.e. mu) */
255  /* and make "space-links" */
256 
257  if ( r == 0 )
258  {
259  u_t = shift(ug[j_decay], FORWARD, mu);
260  u_space = ug[mu];
261  }
262  else
263  {
264  tmp_2 = shift(u_t, FORWARD, mu);
265  u_t = tmp_2;
266 
267  tmp_2 = shift(u_space, FORWARD, mu);
268  tmp_3 = ug[mu] * tmp_2;
269  u_space = tmp_3;
270  }
271 
272  for(t = 0;t < ( lengtht); ++t )
273  {
274  /* Gather space-link from t-direction (i.e. j_decay) */
275 
276  if ( t == 0 )
277  {
278  u_tmp = shift(u_space, FORWARD, j_decay);
279  }
280  else
281  {
282  tmp_2 = shift(u_tmp, FORWARD, j_decay);
283  u_tmp = tmp_2;
284  }
285  tt = lsizet - t - 1;
286  btmp = t_coord < tt;
287 
288  tmp_2 = ug[j_decay] * u_tmp;
289  tmp_3 = tmp_2 * adj(u_t);
290  copymask(tmp_3, btmp, u_tmp);
291  tmp_2 = tmp_3 * adj(u_space);
292  wl_trace = real(trace(tmp_2));
293  wils_loop2[t][r] += sum(wl_trace);
294 
295  } /* end t loop */
296  } /* end r loop */
297  } /* end i loop (for mu) */
298 
299  dummy = 1.0 / double (Layout::vol()*Nc*nspace) ;
300 
301  push(xml,"wils_loop2"); // XML tag for wils_wloop2
302  write(xml, "lengthr", lengthr);
303  write(xml, "lengtht", lengtht);
304  push(xml, "wloop2");
305 
306  multi1d<Double> wloop2(lengtht);
307 
308  for(r = 0; r < lengthr; ++r)
309  {
310  for(t = 0; t < lengtht; ++t)
311  {
312  wils_loop2[t][r] *= dummy;
313  wloop2[t] = wils_loop2[t][r];
314  }
315  push(xml, "elem");
316  write(xml, "r", r);
317  write(xml, "loop", wloop2); // write out wils_loop2
318  pop(xml); // elem
319  } // end for r
320 
321  pop(xml); // XML end tag for loop2
322  pop(xml); // XML end tag for wils_loop2
323  QDPIO::cout << "wils_loop2 data written to .xml file " << std::endl;
324  } /* end on option "time-like planar Wilson loops" */
325 
326  /* Compute "time-like" off-axis Wilson loops, if desired */
327  if ( (kind & 4) != 0 )
328  {
329  int length;
330  if( nspace < 2 )
331  {
332  for(mu = 0;mu < ( Nd); ++mu )
333  QDP_error_exit("Wrong lattice size for off-axis Wilson loops: ",
334  mu, nrow[mu]);
335  }
336  else if ( nspace == 2 )
337  {
338  length = lengthr + lengthr/2;
339  }
340  else
341  length = 2 * lengthr + lengthr/2;
342 
343  QDPIO::cout << "wils_loop3: lengtht=" << lengtht << " length=" << length
344  << " nspace=" << nspace << " j_decay=" << j_decay
345  << " t_dir=" << t_dir << std::endl;
346  multi2d<Double> wils_loop3(lengtht, length);
347  wils_loop3 = 0; /* initialize the non-planar Wilson loops */
348 
349  for(j = 1;j < ( nspace); ++j )
350  {
351  nu = space_dir[j];
352 
353  for(i = 0;i < ( j); ++i )
354  {
355  mu = space_dir[i];
356 
357  /*+ */
358  /* Do off-axis "sqrt(2)" loops in (mu,nu)-plane */
359  /*- */
360 
361  /* Make the "corner link" in the (mu,nu) direction */
362  ftmp = 1.0 / 2.0 ;
363  tmp_2 = shift(ug[nu], FORWARD, mu);
364  u_diag = ug[mu] * tmp_2;
365  tmp_2 = shift(ug[mu], FORWARD, nu);
366  u_diag += ug[nu] * tmp_2;
367  u_diag = u_diag * ftmp;
368 
369  for(r = 0;r < ( lengthr); ++r )
370  {
371  /* Gather t-links from space-directions (i.e. mu and nu) */
372  /* and make "space-links" */
373 
374  if ( r == 0 )
375  {
376  tmp_2 = shift(ug[j_decay], FORWARD, mu);
377  u_t = shift(tmp_2, FORWARD, nu);
378  u_space = u_diag;
379  }
380  else
381  {
382  tmp_2 = shift(u_t, FORWARD, mu);
383  u_t = shift(tmp_2, FORWARD, nu);
384  tmp_2 = shift(u_space, FORWARD, mu);
385  tmp_3 = shift(tmp_2, FORWARD, nu);
386  u_space = u_diag * tmp_3;
387  }
388 
389  for(t = 0;t < ( lengtht); ++t )
390  {
391  /* Gather space-link from t-direction (i.e. j_decay) */
392  if ( t == 0 )
393  {
394  u_tmp = shift(u_space, FORWARD, j_decay);
395  }
396  else
397  {
398  tmp_2 = shift(u_tmp, FORWARD, j_decay);
399  u_tmp = tmp_2;
400  }
401 
402  tt = lsizet - t - 1;
403  btmp = t_coord < tt;
404  tmp_2 = ug[j_decay] * u_tmp;
405  tmp_3 = tmp_2 * adj(u_t);
406  copymask(tmp_3, btmp, u_tmp);
407  tmp_2 = tmp_3 * adj(u_space);
408  wl_trace = real(trace(tmp_2));
409  wils_loop3[t][r] += sum(wl_trace);
410  } /* end t loop */
411  } /* end r loop */
412 
413  /*+ */
414  /* Do off-axis "sqrt(2)" loops in (mu,-nu)-plane */
415  /*- */
416 
417  /* Make the "corner link" in the (mu,-nu) direction */
418  ftmp = 1.0 / 2.0;
419  tmp_2 = adj(ug[nu]) * ug[mu];
420  u_diag = shift(tmp_2, BACKWARD, nu);
421  tmp_2 = shift(ug[nu], FORWARD, mu);
422  tmp_3 = shift(tmp_2, BACKWARD, nu);
423  u_diag += ug[mu] * adj(tmp_3);
424  u_diag = u_diag * ftmp;
425 
426  for(r = 0;r < ( lengthr); ++r )
427  {
428  /* Gather t-links from space-directions (i.e. mu and -nu) */
429  /* and make "space-links" */
430 
431  if ( r == 0 )
432  {
433  tmp_2 = shift(ug[j_decay], FORWARD, mu);
434  u_t = shift(tmp_2, BACKWARD, nu);
435  u_space = u_diag;
436  }
437  else
438  {
439  tmp_2 = shift(u_t, FORWARD, mu);
440  u_t = shift(tmp_2, BACKWARD, nu);
441  tmp_2 = shift(u_space, FORWARD, mu);
442  tmp_3 = shift(tmp_2, BACKWARD, nu);
443  u_space = u_diag * tmp_3;
444  }
445 
446  for(t = 0;t < ( lengtht); ++t )
447  {
448  /* Gather space-link from t-direction (i.e. j_decay) */
449  if ( t == 0 )
450  {
451  u_tmp = shift(u_space, FORWARD, j_decay);
452  }
453  else
454  {
455  tmp_2 = shift(u_tmp, FORWARD, j_decay);
456  u_tmp = tmp_2;
457  }
458 
459  tt = lsizet - t - 1;
460  btmp = t_coord < tt;
461  tmp_2 = ug[j_decay] * u_tmp;
462  tmp_3 = tmp_2 * adj(u_t);
463  copymask(tmp_3, btmp, u_tmp);
464  tmp_2 = tmp_3 * adj(u_space);
465  wl_trace = real(trace(tmp_2));
466  wils_loop3[t][r] += sum(wl_trace);
467 
468  } /* end t loop */
469  } /* end r loop */
470 
471  /*+ */
472  /* Do off-axis "sqrt(5)" loops in (2*mu,nu)-plane */
473  /*- */
474  r_off = lengthr;
475 
476  /* Make the "corner link" in the (mu,nu) direction */
477  ftmp = 1.0 / 2.0;
478 
479  tmp_2 = shift(ug[nu], FORWARD, mu);
480  tmp_3 = ug[mu] * tmp_2;
481  tmp_2 = shift(tmp_3, FORWARD, mu);
482  u_diag = ug[mu] * tmp_2;
483 
484  tmp_2 = shift(ug[mu], FORWARD, mu);
485  tmp_3 = ug[mu] * tmp_2;
486  tmp_2 = shift(tmp_3, FORWARD, nu);
487  u_diag += ug[nu] * tmp_2;
488 
489  u_diag = u_diag * ftmp;
490 
491  for(r = 0;r < ( lengthr/2); ++r )
492  {
493  /* Gather t-links from space-directions (i.e. 2*mu and nu) */
494  /* and make "space-links" */
495 
496  if ( r == 0 )
497  {
498  tmp_2 = shift(ug[j_decay], FORWARD, mu);
499  tmp_3 = shift(tmp_2, FORWARD, mu);
500  u_t = shift(tmp_3, FORWARD, nu);
501  u_space = u_diag;
502  }
503  else
504  {
505  tmp_2 = shift(u_t, FORWARD, mu);
506  tmp_3 = shift(tmp_2, FORWARD, mu);
507  u_tmp = shift(tmp_3, FORWARD, nu);
508  u_t = u_tmp;
509 
510  tmp_2 = shift(u_space, FORWARD, mu);
511  tmp_3 = shift(tmp_2, FORWARD, mu);
512  tmp_2 = shift(tmp_3, FORWARD, nu);
513  u_tmp = u_diag * tmp_2;
514  u_space = u_tmp;
515  }
516 
517  for(t = 0;t < ( lengtht); ++t )
518  {
519  /* Gather space-link from t-direction (i.e. j_decay) */
520  if ( t == 0 )
521  {
522  u_tmp = shift(u_space, FORWARD, j_decay);
523  }
524  else
525  {
526  tmp_2 = shift(u_tmp, FORWARD, j_decay);
527  u_tmp = tmp_2;
528  }
529 
530  tt = lsizet - t - 1;
531  btmp = t_coord < tt;
532  tmp_2 = ug[j_decay] * u_tmp;
533  tmp_3 = tmp_2 * adj(u_t);
534  copymask(tmp_3, btmp, u_tmp);
535  tmp_2 = tmp_3 * adj(u_space);
536  wl_trace = real(trace(tmp_2));
537  wils_loop3[t][r_off+r] += sum(wl_trace);
538 
539  } /* end t loop */
540  } /* end r loop */
541 
542  /*+ */
543  /* Do off-axis "sqrt(5)" loops in (2*mu,-nu)-plane */
544  /*- */
545 
546  /* Make the "corner link" in the (2*mu,-nu) direction */
547  ftmp = 1.0 / 2.0;
548 
549  tmp_2 = shift(ug[mu], FORWARD, mu);
550  tmp_3 = ug[mu] * tmp_2;
551  tmp_2 = adj(ug[nu]) * tmp_3;
552  u_diag = shift(tmp_2, BACKWARD, nu);
553 
554  tmp_2 = shift(ug[nu], FORWARD, mu);
555  tmp_3 = shift(tmp_2, BACKWARD, nu);
556  tmp_2 = ug[mu] * adj(tmp_3);
557  tmp_3 = shift(tmp_2, FORWARD, mu);
558  u_diag += ug[mu] * tmp_3;
559 
560  u_diag = u_diag * ftmp;
561 
562  for(r = 0;r < ( lengthr/2); ++r )
563  {
564  /* Gather t-links from space-directions (i.e. 2*mu and -nu) */
565  /* and make "space-links" */
566 
567  if ( r == 0 )
568  {
569  tmp_2 = shift(ug[j_decay], FORWARD, mu);
570  tmp_3 = shift(tmp_2, FORWARD, mu);
571  u_t = shift(tmp_3, BACKWARD, nu);
572  u_space = u_diag;
573  }
574  else
575  {
576  tmp_2 = shift(u_t, FORWARD, mu);
577  tmp_3 = shift(tmp_2, FORWARD, mu);
578  u_tmp = shift(tmp_3, BACKWARD, nu);
579  u_t = u_tmp;
580 
581  tmp_2 = shift(u_space, FORWARD, mu);
582  tmp_3 = shift(tmp_2, FORWARD, mu);
583  tmp_2 = shift(tmp_3, BACKWARD, nu);
584  u_tmp = u_diag * tmp_2;
585  u_space = u_tmp;
586  }
587 
588  for(t = 0;t < ( lengtht); ++t )
589  {
590  /* Gather space-link from t-direction (i.e. j_decay) */
591  if ( t == 0 )
592  {
593  u_tmp = shift(u_space, FORWARD, j_decay);
594  }
595  else
596  {
597  tmp_2 = shift(u_tmp, FORWARD, j_decay);
598  u_tmp = tmp_2;
599  }
600 
601 
602  tt = lsizet - t - 1;
603  btmp = t_coord < tt;
604 
605  tmp_2 = ug[j_decay] * u_tmp;
606  tmp_3 = tmp_2 * adj(u_t);
607  copymask(tmp_3, btmp, u_tmp);
608  tmp_2 = tmp_3 * adj(u_space);
609  wl_trace = real(trace(tmp_2));
610  wils_loop3[t][r_off+r] += sum(wl_trace);
611  } /* end t loop */
612  } /* end r loop */
613 
614 
615  /*+ */
616  /* Do off-axis "sqrt(5)" loops in (2*nu,mu)-plane */
617  /*- */
618 
619  /* Make the "corner link" in the (2*nu,mu) direction */
620  ftmp = 1.0 / 2.0;
621 
622  tmp_2 = shift(ug[mu], FORWARD, nu);
623  tmp_3 = ug[nu] * tmp_2;
624  tmp_2 = shift(tmp_3, FORWARD, nu);
625  u_diag = ug[nu] * tmp_2;
626 
627  tmp_2 = shift(ug[nu], FORWARD, nu);
628  tmp_3 = ug[nu] * tmp_2;
629  tmp_2 = shift(tmp_3, FORWARD, mu);
630  u_diag += ug[mu] * tmp_2;
631 
632  u_diag = u_diag * ftmp;
633 
634  for(r = 0;r < ( lengthr/2); ++r )
635  {
636  /* Gather t-links from space-directions (i.e. 2*nu and mu) */
637  /* and make "space-links" */
638 
639  if ( r == 0 )
640  {
641  tmp_2 = shift(ug[j_decay], FORWARD, nu);
642  tmp_3 = shift(tmp_2, FORWARD, nu);
643  u_t = shift(tmp_3, FORWARD, mu);
644  u_space = u_diag;
645  }
646  else
647  {
648  tmp_2 = shift(u_t, FORWARD, nu);
649  tmp_3 = shift(tmp_2, FORWARD, nu);
650  u_tmp = shift(tmp_3, FORWARD, mu);
651  u_t = u_tmp;
652 
653  tmp_2 = shift(u_space, FORWARD, nu);
654  tmp_3 = shift(tmp_2, FORWARD, nu);
655  tmp_2 = shift(tmp_3, FORWARD, mu);
656  u_tmp = u_diag * tmp_2;
657  u_space = u_tmp;
658  }
659 
660  for(t = 0;t < ( lengtht); ++t )
661  {
662  /* Gather space-link from t-direction (i.e. j_decay) */
663  if ( t == 0 )
664  {
665  u_tmp = shift(u_space, FORWARD, j_decay);
666  }
667  else
668  {
669  tmp_2 = shift(u_tmp, FORWARD, j_decay);
670  u_tmp = tmp_2;
671  }
672 
673  tt = lsizet - t - 1;
674  btmp = t_coord < tt;
675 
676  tmp_2 = ug[j_decay] * u_tmp;
677  tmp_3 = tmp_2 * adj(u_t);
678  copymask(tmp_3, btmp, u_tmp);
679  tmp_2 = tmp_3 * adj(u_space);
680  wl_trace = real(trace(tmp_2));
681  wils_loop3[t][r_off+r] += sum(wl_trace);
682  } /* end t loop */
683  } /* end r loop */
684 
685 
686  /*+ */
687  /* Do off-axis "sqrt(5)" loops in (2*nu,-mu)-plane */
688  /*- */
689 
690  /* Make the "corner link" in the (2*nu,-mu) direction */
691  ftmp = 1.0 / 2.0;
692  tmp_2 = shift(ug[nu], FORWARD, nu);
693  tmp_3 = ug[nu] * tmp_2;
694  tmp_2 = adj(ug[mu]) * tmp_3;
695  u_diag = shift(tmp_2, BACKWARD, mu);
696 
697  tmp_2 = shift(ug[mu], FORWARD, nu);
698  tmp_3 = shift(tmp_2, BACKWARD, mu);
699  tmp_2 = ug[nu] * adj(tmp_3);
700  tmp_3 = shift(tmp_2, FORWARD, nu);
701  u_diag += ug[nu] * tmp_3;
702 
703  u_diag = u_diag * ftmp;
704 
705  for(r = 0;r < ( lengthr/2); ++r )
706  {
707  /* Gather t-links from space-directions (i.e. 2*nu and -mu) */
708  /* and make "space-links" */
709 
710  if ( r == 0 )
711  {
712  tmp_2 = shift(ug[j_decay], FORWARD, nu);
713  tmp_3 = shift(tmp_2, FORWARD, nu);
714  u_t = shift(tmp_3, BACKWARD, mu);
715  u_space = u_diag;
716  }
717  else
718  {
719  tmp_2 = shift(u_t, FORWARD, nu);
720  tmp_3 = shift(tmp_2, FORWARD, nu);
721  u_tmp = shift(tmp_3, BACKWARD, mu);
722  u_t = u_tmp;
723 
724  tmp_2 = shift(u_space, FORWARD, nu);
725  tmp_3 = shift(tmp_2, FORWARD, nu);
726  tmp_2 = shift(tmp_3, BACKWARD, mu);
727  u_tmp = u_diag * tmp_2;
728  u_space = u_tmp;
729  }
730 
731  for(t = 0;t < ( lengtht); ++t )
732  {
733  /* Gather space-link from t-direction (i.e. j_decay) */
734  if ( t == 0 )
735  {
736  u_tmp = shift(u_space, FORWARD, j_decay);
737  }
738  else
739  {
740  tmp_2 = shift(u_tmp, FORWARD, j_decay);
741  u_tmp = tmp_2;
742  }
743 
744  tt = lsizet - t - 1;
745  btmp = t_coord < tt;
746 
747  tmp_2 = ug[j_decay] * u_tmp;
748  tmp_3 = tmp_2 * adj(u_t);
749  copymask(tmp_3, btmp, u_tmp);
750  tmp_2 = tmp_3 * adj(u_space);
751  wl_trace = real(trace(tmp_2));
752  wils_loop3[t][r_off+r] += sum(wl_trace);
753 
754  } /* end t loop */
755  } /* end r loop */
756  } /* end i loop (for mu) */
757  } /* end j loop (for nu) */
758 
759  dummy = 1.0 / double (Layout::vol()*Nc*nspace*(nspace-1)) ;
760  for(t = 0;t < ( lengtht); ++t )
761  for(r = 0;r < ( lengthr); ++r )
762  wils_loop3[t][r] = wils_loop3[t][r] * dummy;
763 
764  dummy = 1.0 / double (Layout::vol()*Nc*2*nspace*(nspace-1)) ;
765  for(t = 0;t < ( lengtht); ++t )
766  for(r = r_off;r < ( r_off+lengthr/2); ++r )
767  wils_loop3[t][r] = wils_loop3[t][r] * dummy;
768 
769  if ( nspace > 2 )
770  {
771  int k;
772  /*+ */
773  /* Do off-axis "sqrt(3)" loops */
774  /*- */
775  r_off = lengthr + lengthr/2;
776 
777  for(k = 2;k < ( nspace); ++k )
778  {
779  rho = space_dir[k];
780 
781  for(j = 1;j < ( k); ++j )
782  {
783  nu = space_dir[j];
784 
785  for(i = 0;i < ( j); ++i )
786  {
787  mu = space_dir[i];
788 
789  /*+ */
790  /* Do off-axis "sqrt(3)" loops in (mu,nu,rho)-plane */
791  /*- */
792 
793  /* Make the "corner link" in the (mu,nu,rho) direction */
794  ftmp = 1.0 / 6.0;
795  tmp_2 = shift(ug[rho], FORWARD, nu);
796  u_tmp = ug[nu] * tmp_2;
797 
798  tmp_2 = shift(ug[nu], FORWARD, rho);
799  u_tmp += ug[rho] * tmp_2;
800 
801  tmp_2 = shift(u_tmp, FORWARD, mu);
802  u_diag = ug[mu] * tmp_2;
803 
804  tmp_2 = shift(ug[rho], FORWARD, mu);
805  u_tmp = ug[mu] * tmp_2;
806 
807  tmp_2 = shift(ug[mu], FORWARD, rho);
808  u_tmp += ug[rho] * tmp_2;
809 
810  tmp_2 = shift(u_tmp, FORWARD, nu);
811  u_diag += ug[nu] * tmp_2;
812 
813  tmp_2 = shift(ug[nu], FORWARD, mu);
814  u_tmp = ug[mu] * tmp_2;
815 
816  tmp_2 = shift(ug[mu], FORWARD, nu);
817  u_tmp += ug[nu] * tmp_2;
818 
819  tmp_2 = shift(u_tmp, FORWARD, rho);
820  u_diag += ug[rho] * tmp_2;
821 
822  u_diag = u_diag * ftmp;
823 
824  for(r = 0;r < ( lengthr); ++r )
825  {
826  /* Gather t-links from space-directions (i.e. mu, nu and rho) */
827  /* and make "space-links" */
828 
829  if ( r == 0 )
830  {
831  tmp_2 = shift(ug[j_decay], FORWARD, mu);
832  tmp_3 = shift(tmp_2, FORWARD, nu);
833  u_t = shift(tmp_3, FORWARD, rho);
834 
835  u_space = u_diag;
836  }
837  else
838  {
839  tmp_2 = shift(u_t, FORWARD, mu);
840  tmp_3 = shift(tmp_2, FORWARD, nu);
841  u_tmp = shift(tmp_3, FORWARD, rho);
842  u_t = u_tmp;
843 
844  tmp_2 = shift(u_space, FORWARD, mu);
845  tmp_3 = shift(tmp_2, FORWARD, nu);
846  tmp_2 = shift(tmp_3, FORWARD, rho);
847  u_tmp = u_diag * tmp_2;
848  u_space = u_tmp;
849  }
850 
851  for(t = 0;t < ( lengtht); ++t )
852  {
853  /* Gather space-link from t-direction (i.e. j_decay) */
854  if ( t == 0 )
855  {
856  u_tmp = shift(u_space, FORWARD, j_decay);
857  }
858  else
859  {
860  tmp_2 = shift(u_tmp, FORWARD, j_decay);
861  u_tmp = tmp_2;
862  }
863 
864  tt = lsizet - t - 1;
865  btmp = t_coord < tt;
866 
867  tmp_2 = ug[j_decay] * u_tmp;
868  tmp_3 = tmp_2 * adj(u_t);
869  copymask(tmp_3, btmp, u_tmp);
870  tmp_2 = tmp_3 * adj(u_space);
871  wl_trace = real(trace(tmp_2));
872  wils_loop3[t][r_off+r] += sum(wl_trace);
873  } /* end t loop */
874  } /* end r loop */
875 
876  /*+ */
877  /* Do off-axis "sqrt(3)" loops in (mu,nu,-rho)-plane */
878  /*- */
879 
880  /* Make the "corner link" in the (mu,nu,-rho) direction */
881  ftmp = 1.0 / 6.0;
882  tmp_2 = shift(ug[nu], FORWARD, mu);
883  u_tmp = ug[mu] * tmp_2;
884 
885  tmp_2 = shift(ug[mu], FORWARD, nu);
886  u_tmp += ug[nu] * tmp_2;
887 
888  tmp_2 = adj(ug[rho]) * u_tmp;
889  u_diag = shift(tmp_2, BACKWARD, rho);
890 
891  tmp_2 = adj(ug[rho]) * ug[nu];
892  u_tmp = shift(tmp_2, BACKWARD, rho);
893 
894  tmp_2 = shift(ug[rho], FORWARD, nu);
895  tmp_3 = shift(tmp_2, BACKWARD, rho);
896  u_tmp += ug[nu] * adj(tmp_3);
897 
898  tmp_2 = shift(u_tmp, FORWARD, mu);
899  u_diag += ug[mu] * tmp_2;
900 
901  tmp_2 = adj(ug[rho]) * ug[mu];
902  u_tmp = shift(tmp_2, BACKWARD, rho);
903 
904  tmp_2 = shift(ug[rho], FORWARD, mu);
905  tmp_3 = shift(tmp_2, BACKWARD, rho);
906  u_tmp += ug[mu] * adj(tmp_3);
907 
908  tmp_2 = shift(u_tmp, FORWARD, nu);
909  u_diag += ug[nu] * tmp_2;
910 
911  u_diag = u_diag * ftmp;
912 
913  for(r = 0;r < ( lengthr); ++r )
914  {
915  /* Gather t-links from space-directions (i.e. mu, nu and -rho) */
916  /* and make "space-links" */
917 
918  if ( r == 0 )
919  {
920  tmp_2 = shift(ug[j_decay], FORWARD, mu);
921  tmp_3 = shift(tmp_2, FORWARD, nu);
922  u_t = shift(tmp_3, BACKWARD, rho);
923  u_space = u_diag;
924  }
925  else
926  {
927  tmp_2 = shift(u_t, FORWARD, mu);
928  tmp_3 = shift(tmp_2, FORWARD, nu);
929  u_tmp = shift(tmp_3, BACKWARD, rho);
930  u_t = u_tmp;
931 
932  tmp_2 = shift(u_space, FORWARD, mu);
933  tmp_3 = shift(tmp_2, FORWARD, nu);
934  tmp_2 = shift(tmp_3, BACKWARD, rho);
935  u_tmp = u_diag * tmp_2;
936  u_space = u_tmp;
937  }
938 
939  for(t = 0;t < ( lengtht); ++t )
940  {
941  /* Gather space-link from t-direction (i.e. j_decay) */
942  if ( t == 0 )
943  {
944  u_tmp = shift(u_space, FORWARD, j_decay);
945  }
946  else
947  {
948  tmp_2 = shift(u_tmp, FORWARD, j_decay);
949  u_tmp = tmp_2;
950  }
951 
952  tt = lsizet - t - 1;
953  btmp = t_coord < tt;
954 
955  tmp_2 = ug[j_decay] * u_tmp;
956  tmp_3 = tmp_2 * adj(u_t);
957  copymask(tmp_3, btmp, u_tmp);
958  tmp_2 = tmp_3 * adj(u_space);
959  wl_trace = real(trace(tmp_2));
960  wils_loop3[t][r_off+r] += sum(wl_trace);
961 
962  } /* end t loop */
963  } /* end r loop */
964 
965  /*+ */
966  /* Do off-axis "sqrt(3)" loops in (mu,-nu,rho)-plane */
967  /*- */
968 
969  /* Make the "corner link" in the (mu,-nu,rho) direction */
970  ftmp = 1.0 / 6.0;
971  tmp_2 = shift(ug[rho], FORWARD, mu);
972  u_tmp = ug[mu] * tmp_2;
973 
974  tmp_2 = shift(ug[mu], FORWARD, rho);
975  u_tmp += ug[rho] * tmp_2;
976 
977  tmp_2 = adj(ug[nu]) * u_tmp;
978  u_diag = shift(tmp_2, BACKWARD, nu);
979 
980  tmp_2 = adj(ug[nu]) * ug[rho];
981  u_tmp = shift(tmp_2, BACKWARD, nu);
982 
983  tmp_2 = shift(ug[nu], FORWARD, rho);
984  tmp_3 = shift(tmp_2, BACKWARD, nu);
985  u_tmp += ug[rho] * adj(tmp_3);
986 
987  tmp_2 = shift(u_tmp, FORWARD, mu);
988  u_diag += ug[mu] * tmp_2;
989 
990  tmp_2 = adj(ug[nu]) * ug[mu];
991  u_tmp = shift(tmp_2, BACKWARD, nu);
992 
993  tmp_2 = shift(ug[nu], FORWARD, mu);
994  tmp_3 = shift(tmp_2, BACKWARD, nu);
995  u_tmp += ug[mu] * adj(tmp_3);
996 
997  tmp_2 = shift(u_tmp, FORWARD, rho);
998  u_diag += ug[rho] * tmp_2;
999 
1000  u_diag = u_diag * ftmp;
1001 
1002  for(r = 0;r < ( lengthr); ++r )
1003  {
1004  /* Gather t-links from space-directions (i.e. mu, -nu and rho) */
1005  /* and make "space-links" */
1006 
1007  if ( r == 0 )
1008  {
1009  tmp_2 = shift(ug[j_decay], FORWARD, mu);
1010  tmp_3 = shift(tmp_2, BACKWARD, nu);
1011  u_t = shift(tmp_3, FORWARD, rho);
1012  u_space = u_diag;
1013  }
1014  else
1015  {
1016  tmp_2 = shift(u_t, FORWARD, mu);
1017  tmp_3 = shift(tmp_2, BACKWARD, nu);
1018  u_tmp = shift(tmp_3, FORWARD, rho);
1019  u_t = u_tmp;
1020 
1021  tmp_2 = shift(u_space, FORWARD, mu);
1022  tmp_3 = shift(tmp_2, BACKWARD, nu);
1023  tmp_2 = shift(tmp_3, FORWARD, rho);
1024  u_tmp = u_diag * tmp_2;
1025  u_space = u_tmp;
1026  }
1027 
1028  for(t = 0;t < ( lengtht); ++t )
1029  {
1030  /* Gather space-link from t-direction (i.e. j_decay) */
1031  if ( t == 0 )
1032  {
1033  u_tmp = shift(u_space, FORWARD, j_decay);
1034  }
1035  else
1036  {
1037  tmp_2 = shift(u_tmp, FORWARD, j_decay);
1038  u_tmp = tmp_2;
1039  }
1040 
1041  tt = lsizet - t - 1;
1042  btmp = t_coord < tt;
1043 
1044  tmp_2 = ug[j_decay] * u_tmp;
1045  tmp_3 = tmp_2 * adj(u_t);
1046  copymask(tmp_3, btmp, u_tmp);
1047  tmp_2 = tmp_3 * adj(u_space);
1048  wl_trace = real(trace(tmp_2));
1049  wils_loop3[t][r_off+r] += sum(wl_trace);
1050 
1051  } /* end t loop */
1052  } /* end r loop */
1053 
1054  /*+ */
1055  /* Do off-axis "sqrt(3)" loops in (mu,-nu,-rho)-plane */
1056  /*- */
1057 
1058  /* Make the "corner link" in the (mu,-nu,-rho) direction */
1059  ftmp = 1.0 / 6.0;
1060 
1061  tmp_2 = adj(ug[nu]) * ug[mu];
1062  u_tmp = shift(tmp_2, BACKWARD, nu);
1063 
1064  tmp_2 = shift(ug[nu], FORWARD, mu);
1065  tmp_3 = shift(tmp_2, BACKWARD, nu);
1066  u_tmp += ug[mu] * adj(tmp_3);
1067 
1068  tmp_2 = adj(ug[rho]) * u_tmp;
1069  u_diag = shift(tmp_2, BACKWARD, rho);
1070 
1071  tmp_2 = adj(ug[rho]) * ug[mu];
1072  u_tmp = shift(tmp_2, BACKWARD, rho);
1073 
1074  tmp_2 = shift(ug[rho], FORWARD, mu);
1075  tmp_3 = shift(tmp_2, BACKWARD, rho);
1076  u_tmp += ug[mu] * adj(tmp_3);
1077 
1078  tmp_2 = adj(ug[nu]) * u_tmp;
1079  tmp_3 = shift(tmp_2, BACKWARD, nu);
1080  u_diag += tmp_3;
1081 
1082  tmp_2 = shift(ug[rho], BACKWARD, rho);
1083  tmp_3 = adj(ug[nu]) * adj(tmp_2);
1084  u_tmp = shift(tmp_3, BACKWARD, nu);
1085 
1086  tmp_2 = shift(ug[nu], BACKWARD, nu);
1087  tmp_3 = adj(ug[rho]) * adj(tmp_2);
1088  tmp_2 = shift(tmp_3, BACKWARD, rho);
1089  u_tmp += tmp_2;
1090 
1091  tmp_2 = shift(u_tmp, FORWARD, mu);
1092  u_diag += ug[mu] * tmp_2;
1093 
1094  u_diag = u_diag * ftmp;
1095 
1096  for(r = 0;r < ( lengthr); ++r )
1097  {
1098  /* Gather t-links from space-directions (i.e. mu, -nu and -rho) */
1099  /* and make "space-links" */
1100 
1101  if ( r == 0 )
1102  {
1103  tmp_2 = shift(ug[j_decay], FORWARD, mu);
1104  tmp_3 = shift(tmp_2, BACKWARD, nu);
1105  u_t = shift(tmp_3, BACKWARD, rho);
1106  u_space = u_diag;
1107  }
1108  else
1109  {
1110  tmp_2 = shift(u_t, FORWARD, mu);
1111  tmp_3 = shift(tmp_2, BACKWARD, nu);
1112  u_tmp = shift(tmp_3, BACKWARD, rho);
1113  u_t = u_tmp;
1114 
1115  tmp_2 = shift(u_space, FORWARD, mu);
1116  tmp_3 = shift(tmp_2, BACKWARD, nu);
1117  tmp_2 = shift(tmp_3, BACKWARD, rho);
1118  u_tmp = u_diag * tmp_2;
1119  u_space = u_tmp;
1120  }
1121 
1122  for(t = 0;t < ( lengtht); ++t )
1123  {
1124  /* Gather space-link from t-direction (i.e. j_decay) */
1125  if ( t == 0 )
1126  {
1127  u_tmp = shift(u_space, FORWARD, j_decay);
1128  }
1129  else
1130  {
1131  tmp_2 = shift(u_tmp, FORWARD, j_decay);
1132  u_tmp = tmp_2;
1133  }
1134 
1135  tt = lsizet - t - 1;
1136  btmp = t_coord < tt;
1137 
1138  tmp_2 = ug[j_decay] * u_tmp;
1139  tmp_3 = tmp_2 * adj(u_t);
1140  copymask(tmp_3, btmp, u_tmp);
1141  tmp_2 = tmp_3 * adj(u_space);
1142  wl_trace = real(trace(tmp_2));
1143  wils_loop3[t][r_off+r] += sum(wl_trace);
1144 
1145  } /* end t loop */
1146  } /* end r loop */
1147  } /* end i loop (for mu) */
1148  } /* end j loop (for nu) */
1149  } /* end k loop (for rho) */
1150 
1151 
1152 
1153 
1154  dummy = 3.0 / double (Layout::vol()*Nc*2*nspace*(nspace-1)*(nspace-2)) ;
1155  for(t = 0;t < ( lengtht); ++t )
1156  for(r = r_off;r < ( r_off+lengthr); ++r )
1157  wils_loop3[t][r] = wils_loop3[t][r] * dummy;
1158 
1159  } /* end if ( nspace > 2 ) */
1160 
1161  push(xml,"wils_loop3"); // XML tag for wils_wloop3
1162  write(xml, "lengthr", lengthr);
1163  write(xml, "lengtht", lengtht);
1164  write(xml, "length", length);
1165  push(xml, "wloop3");
1166 
1167  multi1d<Double> wloop3(wils_loop3.size2());
1168 
1169  for(r = 0; r < wils_loop3.size1(); ++r)
1170  {
1171  for(t = 0; t < wils_loop3.size2(); ++t)
1172  {
1173  wloop3[t] = wils_loop3[t][r];
1174  }
1175  push(xml, "elem");
1176  write(xml, "r", r);
1177  write(xml, "loop", wloop3); // write out wils_wloop3
1178  pop(xml); // elem
1179  } // end for r
1180 
1181  pop(xml); // XML end tag for wloop3
1182  pop(xml); // XML end tag for wils_wloop3
1183  QDPIO::cout << "wils_loop3 data written to .xml file " << std::endl;
1184  } /* end of option "off-axis Wilson loops" */
1185 
1186  QDPIO::cout << "All wils_loop data written to .xml file " << std::endl;
1187 
1188  END_CODE();
1189  } // end of wilslp
1190 
1191 } // end namespace Chroma
1192 
axgauge.h
Axial gauge fixing.
chromabase.h
Primary include file for CHROMA library code.
mu
int mu
Definition: cool.cc:24
nu
int nu
Definition: cool.cc:25
Chroma::write
void write(XMLWriter &xml, const std::string &path, const AsqtadFermActParams &param)
Writer parameters.
Definition: asqtad_fermact_params_s.cc:40
Chroma::axGauge
void axGauge(multi1d< LatticeColorMatrix > &ug, int decay_dir)
Axial gauge fixing.
Definition: axgauge.cc:27
Chroma::wilslp
void wilslp(const multi1d< LatticeColorMatrix > &u, int j_decay, int t_dir, int kind, XMLWriter &xml, const std::string &xml_group)
Calculate Wilson loops.
Definition: wilslp.cc:32
j
unsigned j
Definition: ldumul_w.cc:35
tmp_2
LatticeColorMatrix tmp_2
Definition: meslate.cc:51
j_decay
int j_decay
Definition: meslate.cc:22
tmp_3
LatticeColorMatrix tmp_3
Definition: meslate.cc:52
t
int t
Definition: meslate.cc:37
Nd
Nd
Definition: meslate.cc:74
Chroma
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
Chroma::QDP_error_exit
QDP_error_exit("too many BiCG iterations", n_count, rsd_sq, cp, c, re_rvr, im_rvr, re_a, im_a, re_b, im_b)
Chroma::u
static multi1d< LatticeColorMatrix > u
Definition: syssolver_linop_qop_mg_w.cc:39
Chroma::push
push(xml_out,"Condensates")
Chroma::i
int i
Definition: pbg5p_w.cc:55
Chroma::pop
pop(xml_out)
Chroma::START_CODE
START_CODE()
Chroma::k
int k
Definition: invbicg.cc:119
testing::internal::Double
FloatingPoint< double > Double
Definition: gtest.h:7351
testing::internal::string
::std::string string
Definition: gtest.h:1979
FORWARD
#define FORWARD
Definition: primitives.h:82
BACKWARD
#define BACKWARD
Definition: primitives.h:83
copymask
copymask(lcoord, lbit, ltmp_1, REPLACE)
sum
Double sum
Definition: qtopcor.cc:37
dummy
Real dummy
Definition: qtopcor.cc:36
wilslp.h
Calculate Wilson loops.
Generated by doxygen 1.9.1