CHROMA
group_baryon_run.cc
Go to the documentation of this file.
1 /*
2  * \brief Group Theory based Baryon Operator Construction with All2All
3  */
4 #include "chroma.h"
5 #include "chromabase.h"
11 using namespace Chroma;
12 using namespace GroupBaryonOperatorEnv;
13 
14  void write( XMLWriter& xml, const std::string& path,
16  {
17  push( xml, path );
18  write( xml, "BaryonOpElement", input );
19  pop( xml );
20  }
21  void write( XMLWriter& xml, const std::string& path,
23  {
24  push( xml, path );
25  write( xml, "BaryonOpIndex", input );
26  pop( xml );
27  }
28  void write( XMLWriter& xml, const std::string& path,
30  {
31  push( xml, path );
32  write( xml, "BaryonOpIns", input );
33  pop( xml );
34  }
35  void write( XMLWriter& xml, const std::string& path,
37  {
38  push( xml, path );
39  write( xml, "gauge_id", input.gauge_id );
40  //write( xml, "BaryonOp", input.Bops );
41  pop( xml );
42  }
43  //! BaryonOperator header writer
44  void write( XMLWriter& xml, const std::string& path, const BaryonOperator_t& param )
45  {
46  if( path != "." ) push( xml, path );
47  int version = 1;
48  write( xml, "version", version );
49  write( xml, "mom2_max", param.mom2_max );
50  write( xml, "j_decay", param.j_decay );
51  write( xml, "seed_l", param.seed_l );
52  write( xml, "seed_m", param.seed_m );
53  write( xml, "seed_r", param.seed_r );
54  write( xml, "perms", param.quarkOrderings );
55  if( path != "." ) pop( xml );
56  }
57  bool linkageHack(void)
58  {
59  bool foo = true;
60  // Inline Measurements
62  return foo;
63  }
64 //
65 //! Make a group theory based baryon operator with all-to-all quark
66 //! propagators using the dilution method of TrinLat (2004)
67 //
68 /*! \defgroup GroupBaryonOperator
69  * \ingroup main
70  *
71  * Main program for making a group baryon operator
72  */
73 int main( int argc, char *argv[] )
74 {
75  // ==================================
76  // Put the machine into a known state
77  // ==================================
78  Chroma::initialize( &argc, &argv );
79  START_CODE();
80  QDPIO::cout << "GroupBaryonAll2All: Group Baryon Operators with All2All" << std::endl;
81  QDPIO::cout << "Linkage = " << linkageHack() << std::endl;
82  StopWatch swatch;
83  // ===============
84  // Read input data
85  // ===============
86  // Instantiate xml reader
87  XMLReader xml_in( Chroma::getXMLInputFileName() );
88  const std::string path = "/GroupBaryonAll2All";
89  Params params( xml_in, path );
90  // ==================
91  // Setup lattice size
92  // ==================
93  Layout::setLattSize( params.nrow );
94  Layout::create();
95  // ============
96  // XML business
97  // ============
98  XMLFileWriter& xml_out = Chroma::getXMLOutputInstance(); // XMLDAT by default
99  push( xml_out, "GroupBaryonAll2AllRun" );
100  proginfo( xml_out ); // Print out basic program info
101  write( xml_out, "InputCopy", xml_in ); // save a copy of the input
102  xml_out.flush();
103  // ============
104  // Startup gauge
105  // =============
106  XMLReader gauge_file_xml, gauge_xml;
107  gaugeStartup(gauge_file_xml, gauge_xml, params.gaugestuff.u, params.gaugestuff.cfg);
108  // test the config
111  QDPIO::cout<< "Gauge field measurements " << w_plaq <<" "<< s_plaq <<" "<< t_plaq <<" "<< link <<" "<<std::endl;
112 
113  // Smear the gauge field if needed
114  swatch.start();
115  std::istringstream xml_l( params.gaugestuff.link_smearing.xml );
116  XMLReader linktop( xml_l );
117  const std::string link_path = "/LinkSmearing";
118  QDPIO::cout << "Link smearing type = " << params.gaugestuff.link_smearing.id << std::endl;
120  linkSmearing( TheLinkSmearingFactory::Instance().createObject(
122  linktop,
123  link_path ) );
124  ( *linkSmearing ) ( params.gaugestuff.u );
125  swatch.stop();
126  QDPIO::cout << "Gauge links smeared : time = " << swatch.getTimeInSeconds() << " secs" << std::endl;
128  QDPIO::cout<< "gauge meas " << w_plaq <<" "<< s_plaq <<" "<< t_plaq <<" "<< link <<" "<<std::endl;
129  // ====================
130  // qqq's and baryonop's
131  // ====================
132  multi1d<GroupBaryonQQQ> AQQQ; // annihilation
133  multi1d<GroupBaryonQQQ> CQQQ; // creation
134  multi1d<GroupBaryonOp> AOp;
135  multi1d<GroupBaryonOp> COp;
136  // ====================================
137  // Read in the text file produced by
138  // Gen_Input_For_All2All_Baryons2.pl
139  // called "Chroma_Input.txt"
140  // ====================================
141  ReadTextInput( params, AOp, COp, AQQQ, CQQQ );
142  // ===============
143  // Noise Solutions
144  // ===============
145  multi1d<Params::QuarkSourceSolutions_t> quarks( params.qprop.solns.size() );
146  try
147  {
148  swatch.start();
149  for(int n=0; n < quarks.size(); ++n)
150  {
151  QDPIO::cout << "Attempt to read solutions for source number=" << n << std::endl;
152  quarks[n].dilutions.resize( params.qprop.solns[n].soln_file_names.size() );
153  QDPIO::cout << "dilutions.size= " << quarks[n].dilutions.size() << std::endl;
154  for(int i=0; i < quarks[n].dilutions.size(); ++i)
155  {
156  XMLReader file_xml, record_xml;
157  QDPFileReader from(file_xml, params.qprop.solns[n].soln_file_names[i], QDPIO_SERIAL);
158  read(from, record_xml, quarks[n].dilutions[i].soln);
159  close(from);
160  // need the seed
161  read(record_xml, "/Propagator/PropSource/Source/ran_seed", quarks[n].seed);
162  params.dilution[ n ].ran_seed = quarks[ n ].seed;
163  // read headers
164  read(record_xml, "/Propagator/PropSource", quarks[n].dilutions[i].source_header);
165  read(record_xml, "/Propagator/ForwardProp", quarks[n].dilutions[i].prop_header);
166  }
167  }
168  swatch.stop();
169  }
170  catch (const std::string& e)
171  {
172  QDPIO::cerr << "Error extracting headers: " << e << std::endl;
173  QDP_abort(1);
174  }
175  QDPIO::cout << "Read SOLUTION vectors from files : time = " << swatch.getTimeInSeconds() << " secs" << std::endl;
176  // ====================================
177  //
178  // Main Part of Program
179  //
180  // Run through the QQQ combinations and
181  // multiply by appropriate coefficients
182  // to make the various operators
183  //
184  // ====================================
185  multi1d<LatticeComplex> vresult(1);
186  LatticeComplex result;
187 #ifdef REDUCETOTIMEDILUTION
188  multi1d<Params::QuarkSourceSolutions_t> quarks2( params.qprop.solns.size() );
189 #endif
190  // don't average over equivalent momenta
191  SftMom phases(params.mom2_max, true, params.j_decay);
192  multi2d<DComplex> elem( params.Nmomenta, params.nrow[3] );
193 
194 #ifdef MAKE_SINK_OPERATORS
195 
196 #ifdef REDUCETOTIMEDILUTION
197  // =========================================
198  //
199  // Change the dilution scheme to a lower one
200  //
201  // =========================================
202  // From TfSxCfGf ==> TfSxCxGx hard-wired for test
203  try
204  {
205  swatch.start();
206  for(int n=0; n < quarks.size(); ++n)
207  {
208  QDPIO::cout << "Reducing spin-colour dilution levels for quark " << n << std::endl;
209  int indx, indx2;
210  //
211  // Fix these for now ... temporary fix
212  //
213  int NTimeDilution = params.nrow[3];
214  int NSpinDilution = 4;
215  int NColourDilution = 3;
216  int NSpaceDilution = 1;
217  quarks2[n].dilutions.resize( quarks[n].dilutions.size()/(params.NdilReduce) );
218  quarks2[n].seed = quarks[n].seed;
219  for(int i=0; i < (NTimeDilution); ++i)
220  {
221  // fix this
222  quarks2[n].dilutions[i].source_header = quarks[n].dilutions[i*(params.NdilReduce)].source_header;
223  quarks2[n].dilutions[i].prop_header = quarks[n].dilutions[i*(params.NdilReduce)].prop_header;
224  }
225  // compute the indices needed assuming they are
226  // in the order as below, time--colour--spin--space
227  // index = x + Nx * ( s + Ns * ( c + Nc * ( t ) ) )
228  for(int t=0; t < NTimeDilution; ++t)
229  {
230  for(int x=0; x < NSpaceDilution; ++x)
231  {
232  indx2 = x + NSpaceDilution * ( t );
233  LatticeFermion Q=zero;
234  for(int c=0; c < NColourDilution; ++c)
235  {
236  for(int s=0; s < NSpinDilution; ++s)
237  {
238  indx = x + NSpaceDilution * ( s + NSpinDilution * ( c + NColourDilution * ( t ) ) );
239  //indx = x + NSpaceDilution * s + NSpaceDilution * NSpinDilution * c + NSpaceDilution * NSpinDilution * NColourDilution * t;
240  Q += quarks[n].dilutions[ indx ].soln;
241  }
242  }
243  quarks2[n].dilutions[ indx2 ].soln = Q;
244  }
245  }
246  }
247  } // try
248  catch ( const std::string& e )
249  {
250  QDPIO::cerr << ": Caught Exception removing COLOUR and SPIN-dilution: " << e << std::endl;
251  QDP_abort( 1 );
252  }
253  swatch.stop();
254  QDPIO::cout << "COLOUR and SPIN dilution removed: time= " << swatch.getTimeInSeconds() << " secs" << std::endl;
255 #endif
256  try
257  { //
258  // Annihilation Operator (plus) first (only one term)
259  //
260  QDPIO::cout << "Making the SINK operators" << std::endl;
261  int ordering = 0;
262  for(int n=0; n < params.Noperators; ++n)
263  {
264 #ifdef REDUCETOTIMEDILUTION
265  AOp[ n ].baryonoperator.seed_l = quarks2[0].seed;
266  AOp[ n ].baryonoperator.seed_m = quarks2[1].seed;
267  AOp[ n ].baryonoperator.seed_r = quarks2[2].seed;
268 #else
269  AOp[ n ].baryonoperator.seed_l = quarks[0].seed;
270  AOp[ n ].baryonoperator.seed_m = quarks[1].seed;
271  AOp[ n ].baryonoperator.seed_r = quarks[2].seed;
272 #endif
273  }
274  swatch.start();
275  // Sink smear all the quarks up-front
276  std::istringstream xml_s( params.sink_smearing.sink.xml );
277  XMLReader sinktop( xml_s );
279  sinkSmearing(
281  params.sink_smearing.sink.id, sinktop,
283 #ifdef REDUCETOTIMEDILUTION
284  for(int n=0; n < quarks2.size(); ++n)
285  {
286  for(int i=0; i < (params.NH[ordering][n]/params.NdilReduce); ++i)
287  {
288  ( *sinkSmearing ) ( quarks2[ n ].dilutions[ i ].soln );
289  }
290  }
291 #else
292  for(int n=0; n < quarks.size(); ++n)
293  {
294  for(int i=0; i < (params.NH[ordering][n]/params.NdilReduce); ++i)
295  {
296  ( *sinkSmearing ) ( quarks[ n ].dilutions[ i ].soln );
297  }
298  }
299 #endif
300  swatch.stop();
301  QDPIO::cout << "SINK smearings done: time= " << swatch.getTimeInSeconds() << " secs" << std::endl;
302 
303  swatch.start();
304  for(int qqq=0; qqq < params.NQQQs; ++qqq)
305  {
306  for(int i=0; i < (params.NH[ordering][0]/params.NdilReduce); ++i)
307  for(int j=0; j < (params.NH[ordering][1]/params.NdilReduce); ++j)
308  for(int k=0; k < (params.NH[ordering][2]/params.NdilReduce); ++k)
309  {
310 #ifdef REDUCETOTIMEDILUTION
311  vresult = AQQQ[ qqq ](
312  quarks2[ 0 ].dilutions[ i ].soln,
313  quarks2[ 1 ].dilutions[ j ].soln,
314  quarks2[ 2 ].dilutions[ k ].soln,
315  PLUS
316  );
317 #else
318  vresult = AQQQ[ qqq ](
319  quarks[ 0 ].dilutions[ i ].soln,
320  quarks[ 1 ].dilutions[ j ].soln,
321  quarks[ 2 ].dilutions[ k ].soln,
322  PLUS
323  );
324 #endif
325  elem = phases.sft( vresult[ 0 ] );
326  for(int b=0; b < AQQQ[ qqq ].NBaryonOps; ++b)
327  {
328  for(int p=0; p < params.Nmomenta; ++p)
329  {
330  for(int t=0; t < params.nrow[ 3 ]; ++t)
331  {
332  AQQQ[ qqq ].baryon[ b ]->baryonoperator.orderings[ 0 ].op( i, j, k ).ind[ 0 ].elem( p, t )
333  += ( AQQQ[ qqq ].coef[ b ] * elem( p, t ) );
334  } // t
335  } // p
336  } // Bop_index b
337  } // i,j,k
338  } // qqq
339  swatch.stop();
340  } // try
341  catch ( const std::string& e )
342  {
343  QDPIO::cerr << ": Caught Exception creating sink baryon operator: " << e << std::endl;
344  QDP_abort( 1 );
345  }
346  QDPIO::cout << "SINK operators done: time= " << swatch.getTimeInSeconds() << " secs" << std::endl;
347  // Save the operators
348  // ONLY SciDAC output format is supported!
349  swatch.start();
350  for(int b=0; b < params.Noperators; ++b)
351  {
352  XMLBufferWriter file_xml;
353  push( file_xml, params.Names[ b ] );
354  file_xml << params.param.baryon_operator;
355  write( file_xml, "Config_info", gauge_xml );
356  pop( file_xml );
357  std::string output=params.Names[ b ]+"_snk";
358 
359  QDPFileWriter to( file_xml, output,
360  QDPIO_SINGLEFILE, QDPIO_SERIAL, QDPIO_OPEN );
361  /*
362  // Write the scalar data
363  {
364  XMLBufferWriter record_xml;
365  write( record_xml, "SourceBaryonOperator", COp[ b ].baryonoperator );
366  write( to, record_xml, COp[ b ].baryonoperator.serialize() );
367  }
368  */
369  // Write the scalar data
370  {
371  XMLBufferWriter record_xml;
372  write( record_xml, "SinkBaryonOperator", AOp[ b ].baryonoperator );
373  write( to, record_xml, AOp[ b ].baryonoperator.serialize() );
374  }
375  close( to );
376  }
377  swatch.stop();
378  QDPIO::cout << "Operators written: time= " << swatch.getTimeInSeconds() << " secs" << std::endl;
379 #endif // end annihilation operators
380 
381 #ifdef MAKE_SOURCE_OPERATORS
382  // =============
383  // Noise Sources ... re-generate from the rng seeds
384  // =============
385  QDPIO::cout << "Re-generating the noise SOURCES" << std::endl;
386  try
387  {
388  swatch.start();
389  int N;
390  for(int n=0; n < quarks.size(); ++n)
391  {
392  bool first = true;
393  LatticeFermion quark_noise; // noisy source on entire lattice
394  for(int i=0; i < quarks[ n ].dilutions.size(); ++i)
395  {
396  std::istringstream xml_s( quarks[ n ].dilutions[ i ].source_header.source.xml );
397  XMLReader sourcetop( xml_s );
398  //QDPIO::cout << "Source = " << quarks[n].dilutions[i].source_header.source.id << std::endl;
399  if ( quarks[ n ].dilutions[ i ].source_header.source.id != DiluteZNQuarkSourceConstEnv::name )
400  {
401  QDPIO::cerr << "Expected source_type = " << DiluteZNQuarkSourceConstEnv::name << std::endl;
402  QDP_abort( 1 );
403  }
404  //QDPIO::cout << "Quark num= " << n << " dilution num= " << i << std::endl;
405  DiluteZNQuarkSourceConstEnv::Params srcParams( sourcetop,
406  quarks[ n ].dilutions[ i ].source_header.source.path );
408  if ( first )
409  {
410  first = false;
411  quarks[ 0 ].j_decay = srcParams.j_decay;
412  // Grab N
413  N = srcParams.N;
414  // Set the seed to desired value
415  quarks[ n ].seed = srcParams.ran_seed;
416  QDP::RNG::setrn( quarks[ n ].seed );
417  // Create the noisy quark source on the entire lattice
418  zN_src( quark_noise, N );
419  }
420  // The seeds must always agree - here the seed is the unique id of the source
421  if ( toBool( srcParams.ran_seed != quarks[ n ].seed ) )
422  {
423  QDPIO::cerr << "quark source=" << n << " dilution=" << i << " seed does not match" << std::endl;
424  QDP_abort( 1 );
425  }
426  // The N's must always agree
427  if ( toBool( srcParams.N != N ) )
428  {
429  QDPIO::cerr << "quark source=" << n << " dilution=" << i << " N does not match" << std::endl;
430  QDP_abort( 1 );
431  }
432  // Use a trick here, create the source and subtract it from the global noisy
433  // Check at the end that the global noisy is zero everywhere.
434  // NOTE: the seed will be set every call
435  quarks[ n ].dilutions[ i ].source = srcConst( params.gaugestuff.u );
436  quark_noise -= quarks[ n ].dilutions[ i ].source;
437  } // end for i
438  Double dcnt = norm2( quark_noise );
439  if ( toDouble( dcnt ) != 0.0 ) // problematic - seems to work with unnormalized sources
440  {
441  QDPIO::cerr << "Noise not saturated by all potential solutions: dcnt=" << dcnt << std::endl;
442  QDP_abort( 1 );
443  }
444  } // end for n
445  swatch.stop();
446  } // end try
447  catch ( const std::string& e )
448  {
449  QDPIO::cerr << ": Caught Exception creating source: " << e << std::endl;
450  QDP_abort( 1 );
451  }
452  QDPIO::cout << "SOURCE vectors reconstructed from Seeds : time = " << swatch.getTimeInSeconds() << " secs" << std::endl;
453 
454 #ifdef REDUCETOTIMEDILUTION
455  // =========================================
456  //
457  // Change the dilution scheme to a lower one
458  //
459  // =========================================
460  // From TfSxCfGf ==> TfSxCxGx hard-wired for test
461  try
462  {
463  swatch.start();
464  for(int n=0; n < quarks.size(); ++n)
465  {
466  QDPIO::cout << "REDUCING spin-colour dilution levels for quark " << n << std::endl;
467  int indx, indx2;
468  //
469  // Fix these for now ... temporary fix
470  //
471  int NTimeDilution = params.nrow[3];
472  int NSpinDilution = 4;
473  int NColourDilution = 3;
474  int NSpaceDilution = 1;
475  quarks2[n].dilutions.resize( quarks[n].dilutions.size()/(params.NdilReduce) );
476  quarks2[n].seed = quarks[n].seed;
477  for(int i=0; i < (NTimeDilution); ++i)
478  {
479  // fix this
480  quarks2[n].dilutions[i].source_header = quarks[n].dilutions[i*(params.NdilReduce)].source_header;
481  quarks2[n].dilutions[i].prop_header = quarks[n].dilutions[i*(params.NdilReduce)].prop_header;
482  }
483  // compute the indices needed assuming they are
484  // in the order as below, time--colour--spin--space
485  // index = x + Nx * ( s + Ns * ( c + Nc * ( t ) ) )
486  for(int t=0; t < NTimeDilution; ++t)
487  {
488  for(int x=0; x < NSpaceDilution; ++x)
489  {
490  indx2 = x + NSpaceDilution * ( t );
491  LatticeFermion Q=zero;
492  for(int c=0; c < NColourDilution; ++c)
493  {
494  for(int s=0; s < NSpinDilution; ++s)
495  {
496  indx = x + NSpaceDilution * ( s + NSpinDilution * ( c + NColourDilution * ( t ) ) );
497  Q += quarks[n].dilutions[ indx ].source;
498  }
499  }
500  quarks2[n].dilutions[ indx2 ].source = Q;
501  }
502  }
503  }
504  } // try
505  catch ( const std::string& e )
506  {
507  QDPIO::cerr << ": Caught Exception removing COLOUR and SPIN-dilution: " << e << std::endl;
508  QDP_abort( 1 );
509  }
510  swatch.stop();
511  QDPIO::cout << "COLOUR and SPIN DILUTION REMOVED: time= " << swatch.getTimeInSeconds() << " secs" << std::endl;
512 #endif
513  //
514  // setting the seeds for output (baryon_operator file)
515  // should also put the mass of the quarks in here as well
516  //
517  for(int n=0; n < params.Noperators; ++n)
518  {
519 #ifdef REDUCETOTIMEDILUTION
520  COp[ n ].baryonoperator.seed_l = quarks2[0].seed;
521  COp[ n ].baryonoperator.seed_m = quarks2[1].seed;
522  COp[ n ].baryonoperator.seed_r = quarks2[2].seed;
523 #else
524  COp[ n ].baryonoperator.seed_l = quarks[0].seed;
525  COp[ n ].baryonoperator.seed_m = quarks[1].seed;
526  COp[ n ].baryonoperator.seed_r = quarks[2].seed;
527 #endif
528  }
529  try
530  { //
531  // Creation Operator (MINUS)
532  //
533  for(int ordering=0; ordering < params.NsrcOrderings; ++ordering)
534  {
535  QDPIO::cout << "SOURCE Operator: ordering = " << ordering << std::endl;
536  swatch.start();
537  Real signs;
538  int not_used;
539  const int L=0,M=1,R=2;
540  //
541  // t = quarks[n].dilutions[i].source_header.t_source = timeslice of source
542  //
543  int qL = DilSwap( ordering, 0,1,2, L );
544  int qM = DilSwap( ordering, 0,1,2, M );
545  int qR = DilSwap( ordering, 0,1,2, R );
546  #define FIRSTTYPE
547  switch ( ordering )
548  {
549  case 0: // [012]
550  signs = 2.0; break;
551  case 1: // [210]
552  #ifdef FIRSTTYPE
553  signs = 2.0; break;
554  #else
555  signs = -2.0; break;
556  #endif
557  case 2: // [021]
558  #ifdef FIRSTTYPE
559  signs = -1.0; break;
560  #else
561  signs = 1.0; break;
562  #endif
563  case 3: // [102]
564  #ifdef FIRSTTYPE
565  signs = -1.0; break;
566  #else
567  signs = 1.0; break;
568  #endif
569  case 4: // [120]
570  signs = -1.0; break;
571  case 5: // [201]
572  signs = -1.0; break;
573  }
574  for(int qqq=0; qqq < params.NQQQs; ++qqq)
575  {
576  for(int iL=0; iL < (params.NH[ordering][L]/params.NdilReduce); ++iL)
577  for(int iM=0; iM < (params.NH[ordering][M]/params.NdilReduce); ++iM)
578  for(int iR=0; iR < (params.NH[ordering][R]/params.NdilReduce); ++iR)
579  {
580 #ifdef REDUCETOTIMEDILUTION
581  if( (quarks2[qL].dilutions[iL].source_header.t_source
582  == quarks2[qM].dilutions[iM].source_header.t_source)
583  &&(quarks2[qM].dilutions[iM].source_header.t_source
584  == quarks2[qR].dilutions[iR].source_header.t_source) )
585 #else
586  if( (quarks[qL].dilutions[iL].source_header.t_source
587  == quarks[qM].dilutions[iM].source_header.t_source)
588  &&(quarks[qM].dilutions[iM].source_header.t_source
589  == quarks[qR].dilutions[iR].source_header.t_source) )
590 #endif
591  {
592  int i0 = DilSwapInv( ordering, iL, iM, iR, 0 );
593  int i1 = DilSwapInv( ordering, iL, iM, iR, 1 );
594  int i2 = DilSwapInv( ordering, iL, iM, iR, 2 );
595  for(int t=0; t < phases.getSet().numSubsets(); ++t)
596  {
597  result[ phases.getSet()[t] ]
598 #ifdef REDUCETOTIMEDILUTION
599  = CQQQ[ qqq ](
600  quarks2[ qL ].dilutions[ iL ].source,
601  quarks2[ qM ].dilutions[ iM ].source,
602  quarks2[ qR ].dilutions[ iR ].source,
603  not_used,
604  MINUS
605  );
606 #else
607  = CQQQ[ qqq ](
608  quarks[ qL ].dilutions[ iL ].source,
609  quarks[ qM ].dilutions[ iM ].source,
610  quarks[ qR ].dilutions[ iR ].source,
611  not_used,
612  MINUS
613  );
614 #endif
615  }
616  elem = phases.sft( result );
617  for(int b=0; b < CQQQ[ qqq ].NBaryonOps; ++b)
618  {
619  for(int p=0; p < params.Nmomenta; ++p)
620  {
621  for(int t=0; t < params.nrow[ 3 ]; ++t)
622  {
623  CQQQ[ qqq ].baryon[ b ]->baryonoperator.orderings[ 0 ].op( i0, i1, i2 ).ind[ 0 ].elem( p, t )
624  += ( CQQQ[ qqq ].coef[ b ] * signs * elem( p, t ) );
625  } // t
626  } // p
627  } // Bop_index
628  } // if t_sources are all the same
629  } // iL,iM,iR
630  } // qqq
631  swatch.stop();
632  QDPIO::cout << "SOURCE operator " << ordering << " done: time= " << swatch.getTimeInSeconds() << " secs" << std::endl;
633  } // ordering
634  } // try
635  catch ( const std::string& e )
636  {
637  QDPIO::cerr << ": Caught Exception creating source baryon operator: " << e << std::endl;
638  QDP_abort( 1 );
639  }
640  QDPIO::cout << "All SOURCE operators done" << std::endl;
641 
642  // Save the operators
643  // ONLY SciDAC output format is supported!
644  swatch.start();
645  for(int b=0; b < params.Noperators; ++b)
646  {
647  XMLBufferWriter file_xml;
648  push( file_xml, params.Names[ b ] );
649  file_xml << params.param.baryon_operator;
650  write( file_xml, "Config_info", gauge_xml );
651  pop( file_xml );
652  std::string output=params.Names[ b ]+"_src";
653 
654  QDPFileWriter to( file_xml, output,
655  QDPIO_SINGLEFILE, QDPIO_SERIAL, QDPIO_OPEN );
656  // Write the scalar data
657  {
658  XMLBufferWriter record_xml;
659  write( record_xml, "SourceBaryonOperator", COp[ b ].baryonoperator );
660  write( to, record_xml, COp[ b ].baryonoperator.serialize() );
661  }
662  /*
663  // Write the scalar data
664  {
665  XMLBufferWriter record_xml;
666  write( record_xml, "SinkBaryonOperator", AOp[ b ].baryonoperator );
667  write( to, record_xml, AOp[ b ].baryonoperator.serialize() );
668  }
669  */
670  close( to );
671  }
672  swatch.stop();
673  QDPIO::cout << "Operators written: time= " << swatch.getTimeInSeconds() << " secs" << std::endl;
674 #endif // end of creation operators
675 
676  QDPIO::cout<<"end of group baryon operator calculation"<<std::endl;
677 
678  // Close the namelist output file XMLDAT
679  pop( xml_out ); // GroupBaryonAll2All
680  END_CODE();
681  // Time to bolt
683  exit( 0 );
684 }
Primary include file for CHROMA in application codes.
Primary include file for CHROMA library code.
Random complex Z(N) sources using dilution.
Class for counted reference semantics.
Definition: handle.h:33
Fourier transform phase factor support.
Definition: sftmom.h:35
multi2d< DComplex > sft(const LatticeComplex &cf) const
Do a sumMulti(cf*phases,getSet())
Definition: sftmom.cc:524
const Set & getSet() const
The set to be used in sumMulti.
Definition: sftmom.h:57
static T & Instance()
Definition: singleton.h:432
Construct group baryon operators.
int main(int argc, char *argv[])
void gaugeStartup(XMLReader &gauge_file_xml, XMLReader &gauge_xml, multi1d< LatticeColorMatrix > &u, Cfg_t &cfg)
Initialize the gauge fields.
void read(XMLReader &xml, const std::string &path, GroupBaryonOperatorEnv::Params::Qprop_t::Solutions_t &input)
Reader.
void proginfo(XMLWriter &xml)
Print out basic information about this program.
Definition: proginfo.cc:24
void zN_src(LatticeFermion &a, int N)
Volume source of Z(N) noise.
Definition: zN_src.cc:41
ForwardProp_t prop_header
Params params
unsigned j
Definition: ldumul_w.cc:35
unsigned n
Definition: ldumul_w.cc:36
void setrn(int iseed[4])
Definition: make_seeds.cc:38
int x
Definition: meslate.cc:34
int t
Definition: meslate.cc:37
const std::string name
Name to be used.
int DilSwapInv(int ord, int i, int j, int k, int which)
void ReadTextInput(Params &params, multi1d< GroupBaryonOp > &AB, multi1d< GroupBaryonOp > &CB, multi1d< GroupBaryonQQQ > &AQQQ, multi1d< GroupBaryonQQQ > &CQQQ)
Reader.
void write(XMLWriter &xml, const std::string &path, const GroupBaryonOperatorEnv::Params &param)
int DilSwap(int ord, int i, int j, int k, int which)
bool registerAll()
Register all the factories.
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
LinOpSysSolverMGProtoClover::Q Q
push(xml_out,"Condensates")
Double c
Definition: invbicg.cc:108
int i
Definition: pbg5p_w.cc:55
void initialize(int *argc, char ***argv)
Chroma initialisation routine.
Definition: chroma_init.cc:114
@ MINUS
Definition: chromabase.h:45
@ PLUS
Definition: chromabase.h:45
void finalize(void)
Chroma finalization routine.
Definition: chroma_init.cc:308
bool linkageHack(void)
Definition: const_hmc.cc:660
std::string getXMLInputFileName()
Get input file name.
Definition: chroma_init.cc:88
pop(xml_out)
void MesPlq(const multi1d< LatticeColorMatrixF3 > &u, multi2d< Double > &plane_plaq, Double &link)
Definition: mesplq.cc:83
START_CODE()
Complex b
Definition: invbicg.cc:96
XMLFileWriter & getXMLOutputInstance()
Get xml output instance.
Definition: chroma_init.cc:359
Double zero
Definition: invbicg.cc:106
int k
Definition: invbicg.cc:119
multi1d< LatticeFermion > s(Ncb)
FloatingPoint< double > Double
Definition: gtest.h:7351
::std::string string
Definition: gtest.h:1979
Double link
Definition: pade_trln_w.cc:146
Double t_plaq
Definition: pade_trln_w.cc:145
Double w_plaq
Definition: pade_trln_w.cc:143
Double s_plaq
Definition: pade_trln_w.cc:144
All make sink constructors.
Factory for producing quark prop sinks.
All source smearing.
Factory for producing quark smearing objects.
Random complex Z(N) sources using dilution.
struct Chroma::GroupBaryonOperatorEnv::Params::Qprop_t qprop
struct Chroma::GroupBaryonOperatorEnv::Params::param_t param
struct Chroma::GroupBaryonOperatorEnv::Params::gaugestuff_t gaugestuff