CHROMA
baryon_s.cc
Go to the documentation of this file.
1 /* + */
2 /* $Id: baryon_s.cc,v 3.1 2007-02-22 21:11:49 bjoo Exp $ ($Date: 2007-02-22 21:11:49 $) */
3 
4 /* This routine is specific to staggered fermions! */
5 
6 /* Construct baryon propagator and Write in (pseudo) NAMELIST format */
7 
8 /* quark_propagator -- quark propagator ( Read ) */
9 /* t_source -- cartesian coordinates of the source ( Read ) */
10 /* j_decay -- direction of the exponential decay ( Read ) */
11 /* bc_spec -- boundary condition for spectroscopy ( Read ) */
12 
13 /* ____ */
14 /* \ */
15 /* b(t) = > < b(t_source, 0) b(t + t_source, x) > */
16 /* / */
17 /* ---- */
18 /* x */
19 
20 /* The colour components are contracted with the totally antisymmetric */
21 /* 'tensor' eps(a,b,c) = antisym_tensor(a,b,c). */
22 
23 /*
24  Run through reorder by mcneile
25  Also move output to outside this file.
26 
27 */
28 
29 #include "chromabase.h"
30 
31 namespace Chroma {
32 
33 
34 //! Function object used for constructing the time-slice set
35  class TimeSliceFunc : public SetFunc
36  {
37  public:
38  TimeSliceFunc(int dir): dir_decay(dir) {}
39 
40  int operator() (const multi1d<int>& coordinate) const {return coordinate[dir_decay];}
41  int numSubsets() const {return Layout::lattSize()[dir_decay];}
42 
43  int dir_decay;
44 
45  private:
46  TimeSliceFunc() {} // hide default constructor
47  };
48 
49 
50 
51  void baryon_s(LatticeStaggeredPropagator & quark_propagator_in,
52  multi1d<Complex> & barprop,
53  multi1d<int> & t_source,
54  int j_decay, int bc_spec)
55  {
56  LatticeColorMatrix quark_propagator;
57 
58  quark_propagator =
59  peekSpin(quark_propagator_in,0,0);
60 
61 
62  LatticeComplex b_prop;
63  LatticeComplex uu_quark;
64 
65  LatticeComplex quark_cpt_A ;
66  LatticeComplex quark_cpt_B ;
67 
68  Complex cdummy;
69  int t0;
70  int t;
71  int t_eff;
72  int c;
73  int cc;
74  int tmp;
75 
76  /* Loop counters */
77  int ci_1;
78  int ci_2;
79  int ci_3;
80  int cf_1;
81  int cf_2;
82  int cf_3;
83 
84 
85 
86  /*# 3 dimensional totally antisymmetric tensor */
87 
88  int antisym_tensor[3][3][3] ;
89  for(int ii=0 ; ii < 3 ; ++ii)
90  for(int jj=0 ; jj < 3 ; ++jj)
91  for(int kk=0 ; kk < 3 ; ++kk)
92  antisym_tensor[ii][jj][kk] = 0 ;
93 
94  antisym_tensor[2][1][0] = 1;
95  antisym_tensor[1][2][0] = -1;
96  antisym_tensor[2][0][1] = -1;
97  antisym_tensor[0][2][1] = 1;
98  antisym_tensor[1][0][2] = 1;
99  antisym_tensor[0][1][2] = -1;
100 
101 
102 
103 
104  // --- pasted from mesons_w ---
105  // Create the time-slice set
106  Set timeslice;
107  timeslice.make(TimeSliceFunc(j_decay));
108 
109  // Length of lattice in j_decay direction
110  int length = timeslice.numSubsets();
111 
112 
113  // not sure why this had to be Dcomplex
114  multi1d<DComplex> hsum(length);
115 
116  t0 = t_source[j_decay]; /* Note j_decay = 0 is not permitted! */
117 
118 
119  /* Set baryon propagator to zero */
120  b_prop = 0;
121 
122  ci_3 = 2;
123  ci_1 = 0;
124  ci_2 = 1;
125  c = 1;
126 
127 
128  /* Sum over second sink u-quark colour */
129  for(cf_3 = 0;cf_3 < ( Nc); ++cf_3 )
130  {
131 
132  /* Sum over sink d-u diquark colour */
133  for(cf_1 = 0;cf_1 < ( Nc); ++cf_1 )
134  for(cf_2 = 0;cf_2 < ( Nc); ++cf_2 )
135  {
136  tmp = antisym_tensor[cf_3][cf_2][cf_1];
137  if(tmp != 0) /* otherwise no contribution! */
138  {
139  cc = c * tmp;
140 
141 
142  /* Build the u-u "diquark" [only direct term (23-->23) is needed!] */
143 
144  quark_cpt_A = peekColor(quark_propagator,cf_2,ci_2 );
145  quark_cpt_B = peekColor(quark_propagator,cf_3,ci_3);
146 
147  uu_quark = quark_cpt_A * quark_cpt_B ;
148 
149  /* Finally tie the u-u "diquark" to the d-quark (1) to form the baryon */
150 
151  quark_cpt_A = peekColor(quark_propagator,cf_1,ci_1);
152  switch(cc)
153  {
154  case +1:
155  b_prop += quark_cpt_A * uu_quark;
156  break;
157  case -1:
158  b_prop -= quark_cpt_A * uu_quark;
159  break;
160  }
161 
162 
163  } /* end if cc (sink diquark colour) */
164  } /* end sum sink diquark colour (cf_1, cf_2) */
165  } /* end sum second sink u-quark colour (cf_3) */
166 
167  /* Project on zero momentum: Do a slice-wise sum. */
168  hsum = sumMulti(b_prop, timeslice);
169 
170  for(t = 0;t < ( length); ++t )
171  {
172  t_eff = (t - t0 + length) % length;
173 
174  if ( bc_spec < 0 && (t_eff+t0) >= length)
175  {
176  barprop[t_eff] = -hsum[t] ;
177  }
178  else
179  barprop[t_eff] = hsum[t];
180  }
181 
182 
183 
184  }
185 
186 
187  void baryon_local_s(LatticeStaggeredPropagator & quark_propagator_in,
188  multi1d<Complex> & barprop,
189  multi1d<int> & t_source,
190  int j_decay, int bc_spec)
191  {
192  LatticeColorMatrix quark_propagator;
193 
194  quark_propagator =
195  peekSpin(quark_propagator_in,0,0);
196 
197 
198  LatticeComplex b_prop;
199  LatticeComplex b_prop_local;
200 
201  LatticeComplex uu_quark;
202 
203  LatticeComplex quark_cpt_A ;
204  LatticeComplex quark_cpt_B ;
205 
206  Complex cdummy;
207  int t0;
208  int t;
209  int t_eff;
210  int c;
211  int cc;
212  int tmp;
213 
214  /* Loop counters */
215  int ci_1;
216  int ci_2;
217  int ci_3;
218  int cf_1;
219  int cf_2;
220  int cf_3;
221 
222 
223 
224  /*# 3 dimensional totally antisymmetric tensor */
225 
226  int antisym_tensor[3][3][3] ;
227  for(int ii=0 ; ii < 3 ; ++ii)
228  for(int jj=0 ; jj < 3 ; ++jj)
229  for(int kk=0 ; kk < 3 ; ++kk)
230  antisym_tensor[ii][jj][kk] = 0 ;
231 
232  antisym_tensor[2][1][0] = 1;
233  antisym_tensor[1][2][0] = -1;
234  antisym_tensor[2][0][1] = -1;
235  antisym_tensor[0][2][1] = 1;
236  antisym_tensor[1][0][2] = 1;
237  antisym_tensor[0][1][2] = -1;
238 
239 
240 
241 
242  // --- pasted from mesons_w ---
243  // Create the time-slice set
244  Set timeslice;
245  timeslice.make(TimeSliceFunc(j_decay));
246 
247  // Length of lattice in j_decay direction
248  int length = timeslice.numSubsets();
249 
250 
251  // not sure why this had to be Dcomplex
252  multi1d<DComplex> hsum(length);
253 
254  t0 = t_source[j_decay]; /* Note j_decay = 0 is not permitted! */
255 
256 
257  /* Set baryon propagator to zero */
258  b_prop = 0;
259 
260  ci_3 = 2;
261  ci_1 = 0;
262  ci_2 = 1;
263  c = 1;
264 
265 
266  /* Sum over second sink u-quark colour */
267  for(cf_3 = 0;cf_3 < ( Nc); ++cf_3 )
268  {
269 
270  /* Sum over sink d-u diquark colour */
271  for(cf_1 = 0;cf_1 < ( Nc); ++cf_1 )
272  for(cf_2 = 0;cf_2 < ( Nc); ++cf_2 )
273  {
274  tmp = antisym_tensor[cf_3][cf_2][cf_1];
275  if(tmp != 0) /* otherwise no contribution! */
276  {
277  cc = c * tmp;
278 
279 
280  /* Build the u-u "diquark" [only direct term (23-->23) is needed!] */
281 
282  quark_cpt_A = peekColor(quark_propagator,cf_2,ci_2 );
283  quark_cpt_B = peekColor(quark_propagator,cf_3,ci_3);
284 
285  uu_quark = quark_cpt_A * quark_cpt_B ;
286 
287  /* Finally tie the u-u "diquark" to the d-quark (1) to form the baryon */
288 
289  quark_cpt_A = peekColor(quark_propagator,cf_1,ci_1);
290  switch(cc)
291  {
292  case +1:
293  b_prop += quark_cpt_A * uu_quark;
294  break;
295  case -1:
296  b_prop -= quark_cpt_A * uu_quark;
297  break;
298  }
299 
300 
301  } /* end if cc (sink diquark colour) */
302  } /* end sum sink diquark colour (cf_1, cf_2) */
303  } /* end sum second sink u-quark colour (cf_3) */
304 
305 
306  // now only sum over part of the hypercube
307  LatticeBoolean ltest(true) ;
308 
309  for(int m = 0; m < Nd; ++m)
310  if ( m != j_decay )
311  ltest &= ( (Layout::latticeCoordinate(m) % 2) == 0);
312 
313  b_prop_local = where(ltest, b_prop, LatticeComplex(zero));
314 
315  /* Project on zero momentum: Do a slice-wise sum. */
316  hsum = sumMulti(b_prop_local, timeslice);
317 
318  for(t = 0;t < ( length); ++t )
319  {
320  t_eff = (t - t0 + length) % length;
321 
322  if ( bc_spec < 0 && (t_eff+t0) >= length)
323  {
324  barprop[t_eff] = -hsum[t] ;
325  }
326  else
327  barprop[t_eff] = hsum[t];
328  }
329 
330 
331 
332  }
333 
334 
335  void lower_contract(LatticeComplex & b_prop,
336  const LatticeColorMatrix & quark_propagator_a,
337  const LatticeColorMatrix & quark_propagator_b,
338  const LatticeColorMatrix & quark_propagator_c,
339  int cf_1, int ci_1,
340  int cf_2, int ci_2,
341  int cf_3, int ci_3,
342  int sgn_epsilon,
343  int wick_sign
344  ) ;
345 
346 
347 
348  //
349  // Version needed for smearing
350  //
351  //
352  //
353 
354 
355  void baryon_s(
356  LatticeStaggeredPropagator & quark_propagator_in_a,
357  LatticeStaggeredPropagator & quark_propagator_in_b,
358  LatticeStaggeredPropagator & quark_propagator_in_c,
359  multi1d<Complex> & barprop,
360  multi1d<int> & t_source,
361  int j_decay, int bc_spec)
362  {
363  LatticeColorMatrix quark_propagator_a;
364  LatticeColorMatrix quark_propagator_b;
365  LatticeColorMatrix quark_propagator_c;
366 
367  quark_propagator_a =
368  peekSpin(quark_propagator_in_a,0,0);
369  quark_propagator_b =
370  peekSpin(quark_propagator_in_b,0,0);
371  quark_propagator_c =
372  peekSpin(quark_propagator_in_c,0,0);
373 
374 
375  LatticeComplex b_prop;
376  LatticeComplex b_prop_local;
377 
378  Complex cdummy;
379  int t0;
380  int t;
381  int t_eff;
382  int c;
383  int cc;
384  int tmp;
385 
386  /* Loop counters */
387  int ci_1;
388  int ci_2;
389  int ci_3;
390  int cf_1;
391  int cf_2;
392  int cf_3;
393 
394  StopWatch swatch;
395  swatch.start();
396 
397 
398  /*# 3 dimensional totally antisymmetric tensor */
399 
400  int antisym_tensor[3][3][3] ;
401  for(int ii=0 ; ii < 3 ; ++ii)
402  for(int jj=0 ; jj < 3 ; ++jj)
403  for(int kk=0 ; kk < 3 ; ++kk)
404  antisym_tensor[ii][jj][kk] = 0 ;
405 
406  antisym_tensor[2][1][0] = 1;
407  antisym_tensor[1][2][0] = -1;
408  antisym_tensor[2][0][1] = -1;
409  antisym_tensor[0][2][1] = 1;
410  antisym_tensor[1][0][2] = 1;
411  antisym_tensor[0][1][2] = -1;
412 
413 
414 
415 
416  // --- pasted from mesons_w ---
417  // Create the time-slice set
418  Set timeslice;
419  timeslice.make(TimeSliceFunc(j_decay));
420 
421  // Length of lattice in j_decay direction
422  int length = timeslice.numSubsets();
423 
424 
425  // not sure why this had to be Dcomplex
426  multi1d<DComplex> hsum(length);
427 
428  t0 = t_source[j_decay]; /* Note j_decay = 0 is not permitted! */
429 
430 
431  /* Set baryon propagator to zero */
432  b_prop = 0;
433 
434  ci_3 = 2;
435  ci_1 = 0;
436  ci_2 = 1;
437  c = 1;
438 
439  int total = 0 ;
440  for(ci_3= 0;ci_3 < Nc; ++ci_3)
441  for(ci_1= 0;ci_1 < Nc; ++ci_1)
442  for(ci_2= 0;ci_2 < Nc; ++ci_2)
443  {
444 
445  /* Sum over second sink u-quark colour */
446  for(cf_3 = 0;cf_3 < ( Nc); ++cf_3 )
447  {
448  for(cf_1 = 0;cf_1 < ( Nc); ++cf_1 )
449  for(cf_2 = 0;cf_2 < ( Nc); ++cf_2 )
450  {
451  tmp =
452  antisym_tensor[cf_3][cf_2][cf_1] *
453  antisym_tensor[ci_3][ci_2][ci_1];
454 
455  if(tmp != 0) /* otherwise no contribution! */
456  {
457  ++total ;
458 
459  int sgn_epsilon = c * tmp;
460  int wick_sign = 1 ;
461 
462  lower_contract(b_prop,
463  quark_propagator_a,
464  quark_propagator_b,
465  quark_propagator_c,
466  cf_1, ci_1,
467  cf_2, ci_2,
468  cf_3, ci_3,
469  sgn_epsilon,1);
470 
471  lower_contract(b_prop,
472  quark_propagator_a,
473  quark_propagator_b,
474  quark_propagator_c,
475  cf_1, ci_1,
476  cf_2, ci_3,
477  cf_3, ci_2,
478  sgn_epsilon,-1);
479 
480 
481  lower_contract(b_prop,
482  quark_propagator_a,
483  quark_propagator_b,
484  quark_propagator_c,
485  cf_1, ci_2,
486  cf_2, ci_3,
487  cf_3, ci_1,
488  sgn_epsilon,1);
489 
490  lower_contract(b_prop,
491  quark_propagator_a,
492  quark_propagator_b,
493  quark_propagator_c,
494  cf_1, ci_2,
495  cf_2, ci_1,
496  cf_3, ci_3,
497  sgn_epsilon,-1);
498 
499  lower_contract(b_prop,
500  quark_propagator_a,
501  quark_propagator_b,
502  quark_propagator_c,
503  cf_1, ci_3,
504  cf_2, ci_1,
505  cf_3, ci_2,
506  sgn_epsilon,1);
507 
508  lower_contract(b_prop,
509  quark_propagator_a,
510  quark_propagator_b,
511  quark_propagator_c,
512  cf_1, ci_3,
513  cf_2, ci_2,
514  cf_3, ci_1,
515  sgn_epsilon,-1);
516 
517 
518  } /* end if cc (sink diquark colour) */
519  } /* end sum sink diquark colour (cf_1, cf_2) */
520  } /* end sum second sink u-quark colour (cf_3) */
521 
522  } // end loop over ci_1, ci_2, ci_3
523 
524  std::cout << "Number of Wick contractions " << total << std::endl ;
525 
526  // only sum over part of the origin of each hypercube
527  LatticeBoolean ltest(true) ;
528 
529  for(int m = 0; m < Nd; ++m)
530  if ( m != j_decay )
531  ltest &= ( (Layout::latticeCoordinate(m) % 2) == 0);
532 
533  b_prop_local = where(ltest, b_prop, LatticeComplex(zero));
534 
535  /* Project on zero momentum: Do a slice-wise sum. */
536  hsum = sumMulti(b_prop_local, timeslice);
537 
538  for(t = 0;t < ( length); ++t )
539  {
540  t_eff = (t - t0 + length) % length;
541 
542  if ( bc_spec < 0 && (t_eff+t0) >= length)
543  {
544  barprop[t_eff] = -hsum[t] ;
545  }
546  else
547  barprop[t_eff] = hsum[t];
548  }
549 
550 
551 
552  swatch.stop();
553  double time_in_sec = swatch.getTimeInSeconds();
554 
555  QDPIO::cout << "Baryon_s took " << time_in_sec
556  << " secs" << std::endl;
557 
558 
559 
560  }
561 
562 
563 
564  //
565  // Some utility routines for the baryon code
566  // These may need to be moved someplace else soon.
567  //
568 
569  //
570  // This is meant to be a slow but safe
571  // routine for the Wick contractions.
572  //
573  //
574 
575  void lower_contract(LatticeComplex & b_prop,
576  const LatticeColorMatrix & quark_propagator_a,
577  const LatticeColorMatrix & quark_propagator_b,
578  const LatticeColorMatrix & quark_propagator_c,
579  int cf_1, int ci_1,
580  int cf_2, int ci_2,
581  int cf_3, int ci_3,
582  int sgn_epsilon,
583  int wick_sign
584  )
585  {
586 
587  LatticeComplex uu_quark;
588 
589  LatticeComplex quark_cpt_A ;
590  LatticeComplex quark_cpt_B ;
591 
592 
593  quark_cpt_A = peekColor(quark_propagator_a,cf_2,ci_2 );
594  quark_cpt_B = peekColor(quark_propagator_b,cf_3,ci_3);
595 
596  uu_quark = quark_cpt_A * quark_cpt_B ;
597 
598  /* Finally tie the u-u "diquark" to the d-quark (1) to form the baryon */
599 
600  quark_cpt_A = peekColor(quark_propagator_c,cf_1,ci_1);
601 
602  int cc = sgn_epsilon * wick_sign ;
603 
604  switch(cc)
605  {
606  case +1:
607  b_prop += quark_cpt_A * uu_quark;
608  break;
609  case -1:
610  b_prop -= quark_cpt_A * uu_quark;
611  break;
612  }
613 
614 
615 
616 
617  }
618 
619 
620 } // end namespace Chroma
chromabase.h
Primary include file for CHROMA library code.
Chroma::TimeSliceFunc
Function object used for constructing the time-slice set.
Definition: barQll_w.h:95
Chroma::TimeSliceFunc::numSubsets
int numSubsets() const
Definition: baryon_s.cc:41
Chroma::TimeSliceFunc::operator()
int operator()(const multi1d< int > &coordinate) const
Definition: barQll_w.h:99
Chroma::TimeSliceFunc::TimeSliceFunc
TimeSliceFunc(int dir)
Definition: baryon_s.cc:38
TimeSliceFunc
Function object used for constructing the time-slice set.
Definition: t_follana_io_s.cc:15
bc_spec
int bc_spec
Definition: inline_barspec_db_w.cc:486
t0
int t0
Definition: inline_barspec_db_w.cc:484
t_source
int t_source
Definition: inline_prop_and_matelem_distillation2_w.cc:196
m
static int m[4]
Definition: make_seeds.cc:16
j_decay
int j_decay
Definition: meslate.cc:22
t
int t
Definition: meslate.cc:37
Nd
Nd
Definition: meslate.cc:74
Chroma
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
Chroma::tmp
LatticeFermion tmp
Definition: mespbg5p_w.cc:36
Chroma::baryon_local_s
void baryon_local_s(LatticeStaggeredPropagator &quark_propagator_in, multi1d< Complex > &barprop, multi1d< int > &t_source, int j_decay, int bc_spec)
Definition: baryon_s.cc:187
Chroma::baryon_s
void baryon_s(LatticeStaggeredPropagator &quark_propagator_in, multi1d< Complex > &barprop, multi1d< int > &t_source, int j_decay, int bc_spec)
Definition: baryon_s.cc:51
Chroma::c
Double c
Definition: invbicg.cc:108
Chroma::lower_contract
void lower_contract(LatticeComplex &b_prop, const LatticeColorMatrix &quark_propagator_a, const LatticeColorMatrix &quark_propagator_b, const LatticeColorMatrix &quark_propagator_c, int cf_1, int ci_1, int cf_2, int ci_2, int cf_3, int ci_3, int sgn_epsilon, int wick_sign)
Definition: baryon_s.cc:575
Chroma::zero
Double zero
Definition: invbicg.cc:106
Generated by doxygen 1.9.1