CHROMA
t_fermion_loop_w.cc
Go to the documentation of this file.
1 /*! \file
2  * \brief Main code for generation of disconnected
3  * loops
4  *
5  * This version is for Wilson fermions.
6  * This code should work for su3 or su4.
7  *
8  * I don't fully understand the conventions for kappa values
9  * and masses. I will figure this out when I compare against
10  * the UKQCD code.
11  */
12 
13 #include <iostream>
14 #include <cstdio>
15 
16 #define MAIN
17 
18 // Include everything...
19 #include "chroma.h"
20 
21 
22 using namespace Chroma;
23 
24 
25 // copied from t_ritz.cc
27 
28 
29 void loops(const LatticeFermion &q_source,
30  const LatticeFermion &psi,
31  int length,
32  XMLWriter& xml_gamma,
33  const std::string& xml_tag) ;
34 
35 void z2_src(LatticeFermion& a) ;
36 
37 void z2_src(LatticeFermion& a, int slice, int mu) ;
38 
39 
40 /*
41  * Input
42  */
43 
44 
45 // Parameters which must be determined from the XML input
46 // and written to the XML output
47 struct Param_t
48 {
50  Real Mass; // Staggered mass
51  Real u0; // Tadpole Factor
52 
53  GaugeStartType cfg_type; // storage order for stored gauge configuration
54  PropType prop_type; // storage order for stored propagator
55 
56  SysSolverCGParams invParam;
57 
58  Real GFAccu, OrPara; // Gauge fixing tolerance and over-relaxation param
59  int GFMax; // Maximum gauge fixing iterations
60  int number_sample ; // number of z2_noise samples
61 
62  multi1d<int> nrow;
63  multi1d<int> boundary;
64  multi1d<int> t_srce;
65 };
66 
67 
68 struct Prop_t
69 {
70  std::string source_file;
71  std::string prop_file;
72 };
73 
74 struct Propagator_input_t
75 {
77  Param_t param;
78  Cfg_t cfg;
79  Prop_t prop;
80 };
81 
82 
83 //
84 void read(XMLReader& xml, const std::string& path, Prop_t& input)
85 {
86  XMLReader inputtop(xml, path);
87 
88 // read(inputtop, "source_file", input.source_file);
89  read(inputtop, "prop_file", input.prop_file);
90 }
91 
92 
93 
94 // Reader for input parameters
95 void read(XMLReader& xml, const std::string& path, Propagator_input_t& input)
96 {
97  XMLReader inputtop(xml, path);
98 
99 
100  // First, read the input parameter version. Then, if this version
101  // includes 'Nc' and 'Nd', verify they agree with values compiled
102  // into QDP++
103 
104  // Read in the IO_version
105  try
106  {
107  read(inputtop, "IO_version/version", input.io_version.version);
108  }
109  catch (const std::string& e)
110  {
111  QDPIO::cerr << "Error reading data: " << e << std::endl;
112  throw;
113  }
114 
115 
116  // Currently, in the supported IO versions, there is only a small difference
117  // in the inputs. So, to make code simpler, extract the common bits
118 
119  // Read the uncommon bits first
120  try
121  {
122  XMLReader paramtop(inputtop, "param"); // push into 'param' group
123 
124  switch (input.io_version.version)
125  {
126  /**************************************************************************/
127  case 1 :
128  /**************************************************************************/
129  break;
130 
131  default :
132  /**************************************************************************/
133 
134  QDPIO::cerr << "Input parameter version " << input.io_version.version << " unsupported." << std::endl;
135  QDP_abort(1);
136  }
137  }
138  catch (const std::string& e)
139  {
140  QDPIO::cerr << "Error reading data: " << e << std::endl;
141  throw;
142  }
143 
144 
145  // Read the common bits
146  try
147  {
148  XMLReader paramtop(inputtop, "param"); // push into 'param' group
149 
150  {
151  std::string ferm_type_str;
152  read(paramtop, "FermTypeP", ferm_type_str);
153  if (ferm_type_str == "WILSON") {
155  }
156  }
157 
158  // GTF NOTE: I'm going to switch on FermTypeP here because I want
159  // to leave open the option of treating masses differently.
160  switch (input.param.FermTypeP) {
161  case FERM_TYPE_WILSON :
162 
163  QDPIO::cout << "Compute fermion loops for Wilson fermions" << std::endl;
164 
165  read(paramtop, "Mass", input.param.Mass);
166  read(paramtop, "u0" , input.param.u0);
167  read(paramtop, "number_sample" , input.param.number_sample);
168 
169 #if 0
170  for (int i=0; i < input.param.numKappa; ++i) {
171  if (toBool(input.param.Kappa[i] < 0.0)) {
172  QDPIO::cerr << "Unreasonable value for Kappa." << std::endl;
173  QDPIO::cerr << " Kappa[" << i << "] = " << input.param.Kappa[i] << std::endl;
174  QDP_abort(1);
175  } else {
176  QDPIO::cout << " Spectroscopy Kappa: " << input.param.Kappa[i] << std::endl;
177  }
178  }
179 #endif
180 
181  break;
182 
183  default :
184  QDP_error_exit("Fermion type not supported\n.");
185  }
186 
187  {
188  std::string cfg_type_str;
189  read(paramtop, "cfg_type", cfg_type_str);
190  if (cfg_type_str == "NERSC") {
192  }
193  else if (cfg_type_str == "HOT") {
194  input.param.cfg_type = HOT_START;
195  } else if (cfg_type_str == "COLD") {
196  input.param.cfg_type = COLD_START;
197  } else {
198  QDP_error_exit("Only know NERSC/HOT/COLD files yet");
199  }
200 
201  }
202 
203  {
204  std::string prop_type_str;
205  read(paramtop, "prop_type", prop_type_str);
206  if (prop_type_str == "SZIN") {
208  } else {
209  QDP_error_exit("Dont know non SZIN files yet");
210  }
211  }
212 
213 // read(paramtop, "invType", input.param.invType);
214 // input.param.invParam.invType = CG_INVERTER; //need to fix this
215  read(paramtop, "RsdCG", input.param.invParam.RsdCG);
216  read(paramtop, "MaxCG", input.param.invParam.MaxCG);
217  // read(paramtop, "GFAccu", input.param.GFAccu);
218  // read(paramtop, "OrPara", input.param.OrPara);
219  // read(paramtop, "GFMax", input.param.GFMax);
220 
221  read(paramtop, "nrow", input.param.nrow);
222  read(paramtop, "boundary", input.param.boundary);
223  read(paramtop, "t_srce", input.param.t_srce);
224  }
225  catch (const std::string& e)
226  {
227  QDPIO::cerr << "Error reading data: " << e << std::endl;
228  throw;
229  }
230 
231 
232  // Read in the gauge configuration file name
233  try
234  {
235  read(inputtop, "Cfg", input.cfg);
236  read(inputtop, "Prop", input.prop);
237  }
238  catch (const std::string& e)
239  {
240  QDPIO::cerr << "Error reading data: " << e << std::endl;
241  throw;
242  }
243 }
244 
245 //! Test fermion loops
246 /*! \defgroup t_fermion_loop Test fermion loops
247  * \ingroup testsmain
248  */
249 
250 int main(int argc, char **argv)
251 {
252  // Put the machine into a known state
253  Chroma::initialize(&argc, &argv);
254 
255  // Input parameter structure
256  Propagator_input_t input;
257 
258  // Instantiate xml reader for DATA
259 // XMLReader xml_in("INPUT_W.xml");
260  XMLReader xml_in(Chroma::getXMLInputFileName());
261 
262  // Read data
263  read(xml_in, "/propagator", input);
264 
265  // Specify lattice size, shape, etc.
266  Layout::setLattSize(input.param.nrow);
267  Layout::create();
268 
269  // Read in the configuration along with relevant information.
270  multi1d<LatticeColorMatrix> u(Nd);
271 
272  XMLReader gauge_xml;
273 
274  QDPIO::cout << "Calculation for SU(" << Nc << ")" << std::endl;
275  switch (input.param.cfg_type)
276  {
277  case FILE_START_NERSC :
278  // su3 specific (at the moment)
279  readArchiv(gauge_xml, u, input.cfg.cfg_file);
280  break;
281  case HOT_START :
282  // create a hot configuration
283  for(int dir = 0 ; dir < Nd ; ++dir)
284  {
285  gaussian(u[dir]);
286  reunit(u[dir]) ;
287  }
288  QDPIO::cout << "Hot/Random configuration created" << std::endl;
289  break;
290  default :
291  QDP_error_exit("Configuration type is unsupported.");
292  }
293 
294  // Check if the gauge field configuration is unitarized
295  unitarityCheck(u);
296 
297  // Instantiate XML writer for XMLDAT
298  XMLFileWriter& xml_out = Chroma::getXMLOutputInstance();
299  push(xml_out, "z2_loops");
300 
301  // Write out the input
302  write(xml_out, "Input", xml_in);
303 
304  // Write out the config header
305  write(xml_out, "Config_info", gauge_xml);
306 
307  // Write out the source header
308  // write(xml_out, "Source_info", source_xml);
309 
310  push(xml_out, "Output_version");
311  write(xml_out, "out_version", 1);
312  pop(xml_out);
313 
314  xml_out.flush();
315 
316 
317  // Calculate some gauge invariant observables just for info.
318  MesPlq(xml_out, "Observables", u);
319  xml_out.flush();
320 
321  //
322  // gauge invariance test
323  //
324 
325  // this parameter will be read from the input file
326  bool do_gauge_transform ;
327  do_gauge_transform = false ;
328  // do_gauge_transform = true ;
329 
330 
331  if( do_gauge_transform )
332  {
333  // gauge transform the gauge fields
334  multi1d<LatticeColorMatrix> u_trans(Nd);
335 
336  // create a random gauge transform
337  LatticeColorMatrix v ;
338 
339  gaussian(v);
340  reunit(v) ;
341 
342  for(int dir = 0 ; dir < Nd ; ++dir)
343  {
344  u_trans[dir] = v*u[dir]*adj(shift(v,FORWARD,dir)) ;
345  u[dir] = u_trans[dir] ;
346  }
347 
348  QDPIO::cout << "Random gauge transform done" << std::endl;
349 
350  } // end of gauge transform
351 
352 
353  int j_decay = Nd-1;
354 
355 
356  // set up the calculation of quark propagators
357  // Typedefs to save typing
358  typedef LatticeFermion T;
359  typedef multi1d<LatticeColorMatrix> P;
360  typedef multi1d<LatticeColorMatrix> Q;
361 
362 
363  // Create a fermion BC. Note, the handle is on an ABSTRACT type.
365 
366  //
367  // Initialize fermion action
368  //
369  UnprecWilsonFermAct S_f(cfs,input.param.Mass);
370 
371  // Set up a state for the current u,
372  Handle< FermState<T,P,Q> > state(S_f.createState(u));
373  GroupXML_t inv_param;
374  {
375  XMLBufferWriter xml_buf;
376  write(xml_buf, "InvertParam", input.param.invParam);
377  XMLReader xml_in(xml_buf);
378  inv_param = readXMLGroup(xml_in, "/InvertParam", "invType");
379  }
380  Handle< SystemSolver<LatticeFermion> > qprop(S_f.qprop(state,inv_param));
381 
382  //
383  // Loop over the source color and spin , creating the source
384  // and calling the relevant propagator routines.
385  //
386  //
387  LatticePropagator quark_propagator;
388  XMLBufferWriter xml_buf;
389  int ncg_had = 0;
390  int n_count;
391 
392 
393  multi2d<DComplex> loops(16, input.param.nrow[3]);
394 
395 
396  LatticeFermion q_source, psi;
397 
398 
399  for(int sample = 0 ; sample < input.param.number_sample ; ++sample)
400  {
401  QDPIO::cout << "Inversion for sample = " << sample << std::endl;
402 
403  q_source = zero ;
404  // gaussian(q_source);
405  z2_src(q_source) ;
406  // z2_src(q_source,0,3) ;
407 
408  // DEBUG write out the source
409  // write(xml_out, "q_source", q_source);
410 
411 
412  // initial guess is zero
413  psi = zero; // note this is ``zero'' and not 0
414 
415  // Compute the propagator for given source color/spin
416  // int n_count;
417 
418  SystemSolverResults_t res = (*qprop)(psi, q_source);
419  ncg_had += res.n_count;
420 
421  push(xml_out,"Qprop");
422  write(xml_out, "Mass" , input.param.Mass);
423  write(xml_out, "RsdCG", input.param.invParam.RsdCG);
424  write(xml_out, "Sample" , sample);
425  write(xml_out, "n_count", res.n_count);
426  pop(xml_out);
427 
428 
429  std::string xml_tag = "dis_loops" ;
430  ::loops(q_source,psi,input.param.nrow[3],xml_out,xml_tag) ;
431 
432  } // numnber of samples
433 
434 
435  // write out the loops
436  push(xml_out,"loop_diagrams");
437 
438  pop(xml_out);
439 
440 
441  pop(xml_out);
442  xml_out.close();
443  xml_in.close();
444 
445  // Time to bolt
447 
448  exit(0);
449 }
Primary include file for CHROMA in application codes.
Class for counted reference semantics.
Definition: handle.h:33
Concrete class for all gauge actions with simple boundary conditions.
Definition: simple_fermbc.h:42
Unpreconditioned Wilson fermion action.
Real u0
EXTERN int FermTypeP
int mu
Definition: cool.cc:24
void read(XMLReader &xml, const std::string &path, AsqtadFermActParams &param)
Read parameters.
void write(XMLWriter &xml, const std::string &path, const AsqtadFermActParams &param)
Writer parameters.
void unitarityCheck(const multi1d< LatticeColorMatrixF3 > &u)
Check the unitarity of color matrix in SU(N)
Definition: unit_check.cc:20
void loops(const LatticeFermion &q_source, const LatticeFermion &psi, int length, XMLWriter &xml_gamma, const std::string &xml_tag)
Fermion loop code.
Definition: loops_w.cc:53
FermType
Fermion type.
PropType
Propagator type.
GroupXML_t readXMLGroup(XMLReader &xml_in, const std::string &path, const std::string &type_name)
Read group and return as a std::string.
@ FERM_TYPE_WILSON
void z2_src(LatticeFermion &a)
Z2-source.
Definition: z2_src.cc:53
std::map< std::string, SinkPropContainer_t > prop
multi1d< int > t_srce
int j_decay
Definition: meslate.cc:22
Nd
Definition: meslate.cc:74
multi1d< LatticeColorMatrix > P
Asqtad Staggered-Dirac operator.
Definition: klein_gord.cc:10
gaussian(aux)
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)
static multi1d< LatticeColorMatrix > u
int n_count
Definition: invbicg.cc:78
LinOpSysSolverMGProtoClover::Q Q
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > const multi1d< Real > & Mass
Definition: pbg5p_w.cc:29
push(xml_out,"Condensates")
LinOpSysSolverMGProtoClover::T T
int i
Definition: pbg5p_w.cc:55
void initialize(int *argc, char ***argv)
Chroma initialisation routine.
Definition: chroma_init.cc:114
void finalize(void)
Chroma finalization routine.
Definition: chroma_init.cc:308
void reunit(LatticeColorMatrixF3 &xa)
Definition: reunit.cc:467
Complex a
Definition: invbicg.cc:95
std::string getXMLInputFileName()
Get input file name.
Definition: chroma_init.cc:88
LatticeFermion psi
Definition: mespbg5p_w.cc:35
pop(xml_out)
void MesPlq(const multi1d< LatticeColorMatrixF3 > &u, multi2d< Double > &plane_plaq, Double &link)
Definition: mesplq.cc:83
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > state
Definition: pbg5p_w.cc:28
XMLFileWriter & getXMLOutputInstance()
Get xml output instance.
Definition: chroma_init.cc:359
const WilsonTypeFermAct< multi1d< LatticeFermion > > & S_f
Definition: pbg5p_w.cc:27
Double zero
Definition: invbicg.cc:106
::std::string string
Definition: gtest.h:1979
#define FORWARD
Definition: primitives.h:82
Gauge configuration structure.
Definition: cfgtype_io.h:16
std::string cfg_file
Definition: cfgtype_io.h:18
Hold group xml and type id.
Params for CG inverter.
Holds return info from SystemSolver call.
Definition: syssolver.h:17
Parameters for running program.
Definition: qpropadd.cc:17
GaugeStartType cfg_type
PropType prop_type
SysSolverCGParams invParam
FermType FermTypeP
multi1d< int > t_srce
multi1d< int > nrow
Definition: qpropadd.cc:18
CfgType cfg_type
multi1d< int > boundary
Propagators.
std::string prop_file
IO_version_t io_version
int main(int argc, char **argv)
GaugeStartType
@ HOT_START
@ COLD_START
@ FILE_START_NERSC