Path Integral Quantum Monte Carlo
setup.cpp
Go to the documentation of this file.
1 
9 #include "setup.h"
10 #include "container.h"
11 #include "constants.h"
12 #include "communicator.h"
13 #include "potential.h"
14 #include "wavefunction.h"
15 #include "action.h"
16 #include "move.h"
17 #include "estimator.h"
18 
19 /**************************************************************************/
25 string getList(const vector<string> &options) {
26 
27  ostringstream optionList;
28  std::copy(options.begin(),options.end()-1,
29  std::ostream_iterator<string>(optionList, ", "));
30  optionList << options.back();
31  return optionList.str();
32 }
33 
34 /**************************************************************************/
37 ostream& operator<<(ostream& os, const vector<string>& vals)
38 {
39  for (auto const& val :vals) {
40  if (val != vals.back())
41  os << val << ", ";
42  else
43  os << val;
44  }
45  return os;
46 }
47 
48 // ---------------------------------------------------------------------------
49 // ---------------------------------------------------------------------------
50 // PARAMETER CLASS -----------------------------------------------------------
51 // ---------------------------------------------------------------------------
52 // ---------------------------------------------------------------------------
53 
54 /**************************************************************************/
61 vector<string> Parameters::split(const string &s, char delim) {
62  vector<string> elems;
63  stringstream ss(s);
64  string item;
65  while (getline(ss, item, delim))
66  elems.push_back(item);
67  return elems;
68 }
69 
70 /**************************************************************************/
76 
77  for(auto const& par : params ) {
78  if (!par.second.empty()) {
79  string key = par.first;
80  cout << key << ":\t";
81  if (type.at(key) == typeid(bool))
82  cout << params.count(key);
83  else
84  extract[par.first](par.second);
85  cout << endl;
86  }
87  }
88 }
89 
90 /**************************************************************************/
97 void Parameters::setupCommandLine(boost::ptr_map<string,po::options_description> &_options) {
98 
99  /* We iterate through the parameters map */
100  for (auto & par : params) {
101 
102  /* get the key */
103  string key = par.first;
104 
105  /* Get the correct option class */
106  po::options_description &option = _options[pClass[key]];
107 
108  /* construct a command line label out of the long and short name */
109  string label = key;
110  if (shortName[key] != "")
111  label += "," + shortName[key];
112 
113  /* insert the option for each possible data type */
114  if (type.at(key) == typeid(bool))
115  option.add_options()(label.c_str(),helpMessage[key].c_str());
116  /* option.add_options()(label.c_str(),po::bool_switch()->default_value(false),helpMessage[key].c_str()); */
117  else if (type.at(key) == typeid(double))
118  option.add_options()(label.c_str(),initValue<double>(key),helpMessage[key].c_str());
119  else if (type.at(key) == typeid(int))
120  option.add_options()(label.c_str(),initValue<int>(key),helpMessage[key].c_str());
121  else if (type.at(key) == typeid(uint32))
122  option.add_options()(label.c_str(),initValue<uint32>(key),helpMessage[key].c_str());
123  else if (type.at(key) == typeid(string))
124  option.add_options()(label.c_str(),initValue<string>(key),helpMessage[key].c_str());
125  else if (type.at(key) == typeid(vector<string>)) {
126 
127  if (state[key] == DEFAULTED) {
128  /* a temporary copy of the vector of strings */
129  vector<string> ops = par.second.as<vector<string>>();
130  option.add_options()(label.c_str(), po::value<vector<string>>()->
131  default_value(ops, getList(ops))->composing(), helpMessage[key].c_str());
132  }
133  else
134  option.add_options()(label.c_str(),
135  po::value<vector<string>>()->composing(), helpMessage[key].c_str());
136  }
137  else
138  cerr << "insertOption Failed to find a valid type.";
139  }
140 }
141 
142 /**************************************************************************/
151 void Parameters::update(int argc, char *argv[], po::options_description &cmdLineOptions) {
152 
153  /* create a temporary parameter map */
154  po::variables_map cmdparams;
155 
156  /* Get the values from the command line */
157  po::store(po::parse_command_line(argc, argv, cmdLineOptions), cmdparams);
158  po::notify(cmdparams);
159 
160  /* Now we go through our local parameter map and adjust values, keeping in
161  * mind that we store boolean values */
162  for (auto & par : params) {
163 
164  /* get the key and value */
165  string key = par.first;
166  auto &val = par.second;
167 
168  /* If the parameter exists, copy it over */
169  if (cmdparams.count(key)) {
170 
171  val = cmdparams[key];
172 
173  /* update the state */
174  if (cmdparams[key].empty())
175  state[key] = UNSET;
176  else if (cmdparams[key].defaulted())
177  state[key] = DEFAULTED;
178  else
179  state[key] = SET;
180  }
181  }
182 
183  /* Now we load a potential xml file from disk. It is optional so we use a
184  * try/catch. */
185  if (!params["param_file"].empty()) {
186 
187  try {
188  pt::ptree xmlParams;
189  pt::read_xml(params["param_file"].as<string>(),xmlParams);
190  xmlParams = xmlParams.get_child("pimc");
191 
192  /* Go through our local parameter map and determine what needs to be
193  * set from the xml file */
194  for (auto & par : params) {
195 
196  /* get the key */
197  string key = par.first;
198 
199  /* the parameter class string */
200  string paramClass = pClass[key] + "_parameters";
201 
202  /* Get a clean xml node */
203  pt::ptree xml;
204  if (getNode(xmlParams,xml,paramClass)) {
205 
206  if (type.at(key) == typeid(bool)){
207  set<bool>(key,xml);
208  }
209  else if (type.at(key) == typeid(double)) {
210  set<double>(key,xml);
211  }
212  else if (type.at(key) == typeid(int)) {
213  set<int>(key,xml);
214  }
215  else if (type.at(key) == typeid(uint32)) {
216  set<uint32>(key,xml);
217  }
218  else if (type.at(key) == typeid(string)) {
219  set<string>(key,xml);
220  }
221  else if (type.at(key) == typeid(vector<string>)) {
222 
223  string pluralKey = key + "s";
224  if (xml.count(pluralKey) != 0) {
225 
226  /* a temporary copy of the vector of strings */
227  vector<string> ops;
228 
229  /* get the options */
230  for (auto const &opName : xml.get_child(pluralKey))
231  ops.push_back(opName.second.data());
232 
233  /* update the parameter */
234  set<vector<string>>(key,ops);
235  state[key] = SET;
236  }
237  }
238  else
239  cerr << "xmlInsertOption Failed to find a valid type.";
240  } //if getNode
241  } //for par
242 
243  }
244  catch(...) {
245  cerr << "Cannot read parameter file: " << params["param_file"].as<string>() << endl;
246  }
247  }
248 }
249 
250 // ---------------------------------------------------------------------------
251 // ---------------------------------------------------------------------------
252 // SETUP CLASS ---------------------------------------------------------------
253 // ---------------------------------------------------------------------------
254 // ---------------------------------------------------------------------------
255 //
256 /**************************************************************************/
263  params(),
264  cmdLineOptions("Command Line Options")
265 {
266  /* Initialize the option class names */
267  optionClassNames = {"simulation","physical","cell","algorithm","potential",
268  "measurement"};
269 
270  /* Define the allowed interaction potential names */
271  interactionPotentialName = {"aziz", "szalewicz", "delta", "lorentzian", "sutherland",
272  "hard_sphere", "hard_rod", "free", "delta1D", "harmonic", "dipole"};
273  interactionNames = getList(interactionPotentialName);
274 
275  /* Define the allowed external potential names */
276  externalPotentialName = {"free", "harmonic", "osc_tube", "lj_tube", "plated_lj_tube",
277  "hard_tube", "hg_tube", "fixed_aziz", "gasp_prim", "fixed_lj", "graphene", "graphenelut",
278  "graphenelut3d", "graphenelut3dgenerate", "graphenelut3dtobinary", "graphenelut3dtotext"};
279  externalNames = getList(externalPotentialName);
280 
281  /* Define the allowed action names */
282  actionName = {"primitive", "li_broughton", "gsf", "pair_product"};
283  actionNames = getList(actionName);
284 
285  /* Define the allowed trial wave function names */
286  waveFunctionName = {"constant", "sech", "jastrow", "lieb", "sutherland"};
287  waveFunctionNames = getList(waveFunctionName);
288 
289  /* Define the allowed random number generator names */
290  randomGeneratorName = {"boost_mt19937","std_mt19937", "pimc_mt19937"};
291  randomGeneratorNames = getList(randomGeneratorName);
292 
293  /* Get the allowed estimator names */
294  estimatorName = estimatorFactory()->getNames();
295  vector<string> multiEstimatorName = multiEstimatorFactory()->getNames();
296  estimatorName.insert(estimatorName.end(), multiEstimatorName.begin(),
297  multiEstimatorName.end());
298  estimatorNames = getList(estimatorName);
299 
300  /* Get the allowed move names */
301  moveName = moveFactory()->getNames();
302  moveNames = getList(moveName);
303 }
304 
305 /**************************************************************************/
309 void Setup::initParameters() {
310 
311  /* Instantiate all the option classes */
312  for(string oClass : optionClassNames) {
313  string optionClassName = oClass;
314  optionClassName[0] = toupper(optionClassName[0]);
315  optionClassName += " Options";
316  optionClasses.insert(oClass,new po::options_description(optionClassName));
317  }
318 
319  /* Initialize the simulation options */
320  string oClass = "simulation";
321  params.add<bool>("help,h","produce help message",oClass);
322  params.add<bool>("version","output repo version",oClass);
323  params.add<bool>("validate","validate command line or xml options",oClass);
324  params.add<bool>("dimension","output currently compiled dimension",oClass);
325  params.add<int>("output_config,o","number of output configurations",oClass,0);
326  params.add<uint32>("process,p","process or cpu number",oClass,0);
327  params.add<string>("restart,R","restart running simulation with PIMCID",oClass);
328  params.add<double>("wall_clock,W","set wall clock limit in hours",oClass);
329  params.add<string>("start_with_state,s", "start simulation with a supplied state file.",oClass,"");
330  params.add<bool>("no_save_state","Only save a state file at the end of a simulation",oClass);
331  params.add<bool>("estimator_list","Output a list of estimators in xml format.",oClass);
332  params.add<bool>("update_list","Output a list of updates in xml format.",oClass);
333  params.add<string>("label","a label to append to all estimator files.",oClass,"");
334  params.add<string>("rng,G",str(format("random number generator type:\n%s") % randomGeneratorNames).c_str(),oClass,"pimc_mt19937");
335  params.add<string>("param_file","a valid path to the parameters input xml file.",oClass);
336 
337  /* Initialize the cell options */
338  oClass = "cell";
339  params.add<string>("geometry,b","simulation cell type [prism,cylinder]",oClass,"prism");
340  params.add<double>("size,L","linear system size [angstroms]",oClass);
341  params.add<double>("Lx","system size in x-direction [angstroms]",oClass);
342  params.add<double>("Ly","system size in y-direction [angstroms]",oClass);
343  params.add<double>("Lz","system size in z-direction [angstroms]",oClass);
344  params.add<double>("radius,r","tube radius [angstroms]",oClass);
345 
346  /* Initialize the potential options */
347  oClass = "potential";
348  params.add<string>("interaction,I",str(format("interaction potential type:\n%s") % interactionNames).c_str(),oClass,"aziz");
349  params.add<string>("external,X",str(format("external potential type:\n%s") % externalNames).c_str(),oClass,"free");
350  params.add<double>("scattering_length,a","scattering length [angstroms]",oClass,1.0);
351  params.add<double>("delta_width","delta function potential width",oClass,1.0E-3);
352  params.add<double>("delta_strength,c","delta function potential integrated strength",oClass,10.0);
353  params.add<double>("interaction_strength,g","interaction parameter",oClass,1.0);
354  params.add<double>("omega","harmonic interaction potential frequency",oClass,1.0);
355  params.add<double>("lj_sigma","Lennard-Jones hard-core radius [angstroms]",oClass,2.74);
356  params.add<double>("lj_epsilon","Lennard-Jones energy scale [kelvin]",oClass,16.2463);
357  params.add<double>("lj_width","Radial with of LJ plated cylinder material [angstroms]",oClass);
358  params.add<double>("lj_density","Density LJ plated cylinder material [angstroms^(-3)]",oClass);
359  params.add<double>("hourglass_radius","differential radius for hourglass potential [angstroms]",oClass,0.0);
360  params.add<double>("hourglass_width","constriction width for hourglass potential [angstroms]",oClass,0.0);
361  params.add<string>("fixed,f","input file name for fixed atomic positions.",oClass,"");
362  params.add<double>("potential_cutoff,l","interaction potential cutoff length [angstroms]",oClass);
363  params.add<double>("empty_width_y,y","how much space (in y-) around Gasparini barrier",oClass);
364  params.add<double>("empty_width_z,z","how much space (in z-) around Gasparini barrier",oClass);
365 
366  /* These are graphene potential options */
367  params.add<double>("strain","strain of graphene lattice in y-direction",oClass,0.00);
368  params.add<double>("poisson","Poisson's ratio for graphene",oClass,0.165);
369  params.add<double>("carbon_carbon_dist,A","Carbon-Carbon distance for graphene",oClass,1.42);
370  params.add<string>("graphenelut3d_file_prefix","GrapheneLUT3D file prefix <prefix>serialized.{dat|txt}",oClass,"");
371 
372  /* Initialize the physical options */
373  oClass = "physical";
374  params.add<bool>("canonical","perform a canonical simulation",oClass);
375  params.add<double>("mass,m","particle mass [amu]",oClass,4.0030);
376  params.add<double>("density,n",str(format("initial density [angstroms^(-%d)]") % NDIM).c_str(),oClass);
377  params.add<int>("number_particles,N","number of particles",oClass);
378  params.add<double>("temperature,T","temperature [kelvin]",oClass);
379  params.add<string>("wavefunction",str(format("trial wave function type:\n%s")
380  % waveFunctionNames).c_str(),oClass,"constant");
381  params.add<double>("R_LL_wfn","length scale of the lieb liniger wave function",oClass,0.0);
382  params.add<double>("k_LL_wfn","wave number of the lieb liniger wave function",oClass,0.0);
383  params.add<double>("end_factor","end bead potential action multiplicatave factor",oClass,1.0);
384  params.add<double>("chemical_potential,u","chemical potential [kelvin]",oClass,0.0);
385  params.add<int>("number_paths","number of paths",oClass,1);
386 
387  /* Initialize the algorithm options */
388  oClass = "algorithm";
389  params.add<bool>("relax","perform a worm constant relaxation",oClass);
390  params.add<bool>("relaxmu", "perform a chemical potential relaxation to target a fixed density",oClass);
391  params.add<int>("number_time_slices,P","number of time slices",oClass);
392  params.add<int>("window","set particle number window",oClass);
393  params.add<double>("gaussian_window_width", "set gaussian ensemble weight",oClass);
394  params.add<double>("imaginary_time_step,t", "imaginary time step [kelvin^(-1)]",oClass);
395  params.add<double>("imaginary_time_length","total path length in imaginary time [kelvin^(-1)]",oClass);
396  params.add<double>("worm_constant,C", "worm acceptance constant",oClass,1.0);
397  params.add<double>("com_delta,D", "center of mass update radius[angstroms]",oClass);
398  params.add<double>("displace_delta,d", "displace update radius [angstroms]",oClass);
399  params.add<int>("update_length,M", "non-local update length, (Mbar), -1 sets maximal value",oClass);
400  params.add<string>("action", str(format("action type:\n%s") % actionNames).c_str(),oClass,"gsf");
401  params.add<bool>("full_updates", "perform full staging updates",oClass);
402  params.add<bool>("var_updates", "perform variable length diagonal updates",oClass);
403  params.add<bool>("staging", "perform staging instead of bisection for the diagonal update",oClass);
404  params.add<int>("max_winding", "the maximum winding number to sample",oClass,1);
405 
406  /* Initialize the measurement options */
407  oClass = "measurement";
408  params.add<uint32>("number_eq_steps,E", "number of equilibration steps",oClass,1);
409  params.add<int>("number_bins_stored,S", "number of estimator bins stored",oClass,1);
410  params.add<int>("bin_size", "number of updates per bin",oClass,100);
411  params.add<double>("estimator_radius,w", "maximum radius for cylinder estimators",oClass,2.0);
412  params.add<int>("virial_window,V", "centroid virial energy estimator window",oClass,5);
413  params.add<int>("number_broken", "number of broken world-lines",oClass,0);
414  params.add<double>("spatial_subregion", "define a spatial subregion",oClass);
415 
416  /* The updates, measurement defaults, and ensemble can depend on PIGS vs PIMC */
417  vector<string> estimatorsToMeasure;
418  vector<string> movesToPerform;
419  if (PIGS) {
420  params.set<bool>("canonical",true);
421  estimatorsToMeasure = {EnergyEstimator::name, TimeEstimator::name};
422  movesToPerform = {CenterOfMassMove::name, StagingMove::name, EndStagingMove::name,
423  DisplaceMove::name};
424  }
425  else {
426  estimatorsToMeasure = {EnergyEstimator::name, NumberParticlesEstimator::name,
427  TimeEstimator::name, DiagonalFractionEstimator::name};
428 
429  movesToPerform = {CenterOfMassMove::name, BisectionMove::name, OpenMove::name,
430  CloseMove::name, InsertMove::name, RemoveMove::name, AdvanceHeadMove::name,
431  RecedeHeadMove::name, AdvanceTailMove::name, RecedeTailMove::name, SwapHeadMove::name,
432  SwapTailMove::name};
433 
434  }
435  params.add<vector<string>>("update", str(format("Monte Carlo updates to be performed:\n%s")
436  % moveNames).c_str(), "algorithm",movesToPerform);
437  params.add<vector<string>>("estimator,e", str(format("estimators to be measured:\n%s")
438  % estimatorNames).c_str(),"measurement",estimatorsToMeasure);
439 
440 }
441 
442 /**************************************************************************/
448 string Setup::getXMLOptionList(const vector<string> &options, const string tag) {
449 
450  string optionList = "";
451  for(auto name : options)
452  optionList += str(format("<%s>\n %s\n<\\%s>\n") % tag % name % tag);
453  return optionList;
454 }
455 
456 /**************************************************************************/
464 void Setup::getOptions(int argc, char *argv[])
465 {
466  /* Initialize all possible options */
467  initParameters();
468 
469  /* Insert all options to the appropriate optionClass map for reading
470  * from the command line. */
471  params.setupCommandLine(optionClasses);
472 
473  /* Add all classes to the complete set */
474  for(string oClass : optionClassNames)
475  cmdLineOptions.add(optionClasses[oClass]);
476 
477  /* Update the parameter map from the commamand line and XML */
478  params.update(argc,argv,cmdLineOptions);
479 }
480 
481 /**************************************************************************/
489 
490  /* Do we need help? */
491  if (params("help")) {
492  cout << cmdLineOptions << endl;
493  return true;
494  }
495 
496  /* Output the dimension the code was compiled with then exit */
497  if (params("dimension")) {
498  cout << endl << format("Code was compiled for a %d-dimensional system.") % NDIM
499  << endl << endl;
500  return true;
501  }
502 
503  /* Output the git version that the code was compiled with */
504  if (params("version")) {
505  cout << endl << format("Code was compiled with repo version %s.") % REPO_VERSION
506  << endl << endl;
507  return true;
508  }
509 
510  /* Print a list of allowed estimators in xml format */
511  if (params("estimator_list")) {
512  cout << endl << "The following estimators can be included in an xml input file." << endl;
513  cout << getXMLOptionList(estimatorName, "estimator") << endl << endl;
514  return true;
515  }
516 
517  /* Have we defined a temperature for a PIMC simulation?*/
518  if (!params("temperature") && !PIGS ) {
519  cerr << endl << "ERROR: No temperature defined!" << endl << endl;
520  cerr << "Action: specify temperature (T)" << endl;
521  return true;
522  }
523 
524  /* Have we mistakingly defined one for a PIGS simulation? */
525  if (params("temperature") && PIGS) {
526  cerr << endl << "ERROR: Did you mean to define a non-zero temperature for PIGS?" << endl << endl;
527  cerr << "Action: remove temperature (T)" << endl;
528  return true;
529  }
530 
531  /* Have we physically defined a simulation cell? */
532  dVec side;
533  definedCell = false;
534  if (params("size")) {
535  definedCell = true;
536  side = params["size"].as<double>();
537  params.set<dVec>("side",side);
538 
539  }
540  else if ((params("Lx") + params("Ly") + params("Lz")) == NDIM) {
541  definedCell = true;
542 
543  /* The side labels */
544  vector<string> sides = {"Lx","Ly","Lz"};
545  for (int i = 0; i < NDIM; i++)
546  side[i] = params[sides[i]].as<double>();
547  params.set<dVec>("side",side);
548  }
549 
550  /* Make sure we have defined enough options to create the simulation cell */
551  if (!( (params("density") && params("number_particles")) ||
552  (definedCell && params("number_particles")) ||
553  (definedCell && params("density")) ) ) {
554  cerr << endl << "ERROR: Cannot create the simulation cell!" << endl << endl;
555  cerr << "Action: define a valid simulation cell." << endl;
556  cerr << "Need: [number_particles (N) AND density (n)] OR " << endl;
557  cerr << " [number_particles (N) AND size (L) or Lx,Ly,Lz] OR" << endl;
558  cerr << " [size (L) or Lx,Ly,Lz AND density (n)]" << endl << endl;
559  cerr << optionClasses["cell"] << endl;
560  return true;
561  }
562 
563  /* Make sure we have selected a valid cell type */
564  if (!( (params["geometry"].as<string>() == "cylinder") ||
565  (params["geometry"].as<string>() == "prism") ))
566  {
567  cerr << endl << "ERROR: Invalid simulation cell type." << endl << endl;
568  cerr << "Action: change cell (b) to one of:" << endl
569  << "\t[prism,cylinder]" << endl;
570  return true;
571  }
572 
573  /* Can we create the worldlines? */
574  if (!((params("imaginary_time_length") && params("number_time_slices")) ||
575  (params("imaginary_time_length") && params("imaginary_time_step")) ||
576  (params("number_time_slices") && params("imaginary_time_step")) ||
577  (params("temperature") && (params("number_time_slices") ||
578  (params("imaginary_time_step")))) ) ) {
579  cerr << endl << "ERROR: Cannot create imaginary time paths!" << endl << endl;
580  cerr << "Action: define imaginary_time_length AND imaginary_time_step (t) OR" << endl;
581  cerr << " define imaginary_time_length AND number_time_steps (P) OR" << endl;
582  cerr << " define number_time_slices (P) AND imaginary_time_step (t)" << endl << endl;
583  cerr << " define temperature (T) AND (number_time_slices (P) OR imaginary_time_step (t))" << endl << endl;
584  cerr << optionClasses["algorithm"] << endl;
585  return true;
586  }
587 
588  /* Make sure we have selected a valid interaction potential */
589  if (std::find(interactionPotentialName.begin(), interactionPotentialName.end(),
590  params["interaction"].as<string>()) == interactionPotentialName.end()) {
591  cerr << endl << "ERROR: Invalid interaction potential!" << endl << endl;
592  cerr << "Action: set interaction (I) to one of:" << endl
593  << "\t[" << interactionNames << "]" << endl;
594  return true;
595  }
596 
597  /* Make sure we have selected a valid external potential */
598  if (std::find(externalPotentialName.begin(), externalPotentialName.end(),
599  params["external"].as<string>()) == externalPotentialName.end()) {
600  cerr << endl << "ERROR: Invalid external potential!" << endl << endl;
601  cerr << "Action: set external (X) must to one of:" << endl
602  << "\t[" << externalNames << "]" << endl;
603  return true;
604  }
605 
606  /* Make sure we have selected a valid action type */
607  if (std::find(actionName.begin(), actionName.end(),
608  params["action"].as<string>()) == actionName.end()) {
609  cerr << endl << "ERROR: Invalid action!" << endl << endl;
610  cerr << "Action: set action to one of:" << endl
611  << "\t[" << actionNames << "]" << endl;
612  return true;
613  }
614 
615  /* Make sure we have selected a valid trial wave fucntion */
616  if (std::find(waveFunctionName.begin(), waveFunctionName.end(),
617  params["wavefunction"].as<string>()) == waveFunctionName.end()) {
618  cerr << endl << "ERROR: Invalid trial wave function!" << endl << endl;
619  cerr << "Action: set wave function to one of :" << endl
620  << "\t[" << waveFunctionNames << "]" << endl;
621  return true;
622  }
623 
624  /* Make sure we use the pair product approximation for discontinuous
625  * potentials */
626  if (((params["interaction"].as<string>() == "hard_sphere") ||
627  (params["interaction"].as<string>() == "hard_rod") ||
628  (params["interaction"].as<string>() == "delta1D") ) &&
629  (params["action"].as<string>() != "pair_product") ) {
630  cout << endl;
631  cerr << "ERROR: Need to use the pair product approximation with discontinuous potentials!";
632  cout << endl << endl;
633  cerr << "Action: change the action to pair_product." << endl;
634  return 1;
635  }
636 
637  /* We can only use the cylinder potentials for a 3D system */
638  if ((params["external"].as<string>().find("tube") != string::npos) && (NDIM != 3)) {
639  cerr << endl << "ERROR: Can only use tube potentials for a 3D system!" << endl << endl;
640  cerr << "Action: change the potential or recompile with ndim=3." << endl;
641  return 1;
642  }
643 
644  /* Need to specify a radius for the tube potentials */
645  if ( (params["external"].as<string>().find("tube") != string::npos) && (!params("radius")) ) {
646  cerr << endl << "ERROR: Incomplete specification for external potential!" << endl << endl;
647  cerr << "Action: specfiy a radius (r) for the tube potentials." << endl;
648  return 1;
649  }
650 
651  /* Need to specify a y- barrier width scale factor for Gasparini potential */
652  if ( (params["external"].as<string>().find("gasp_prim") != string::npos) &&
653  (!params("empty_width_y")) ) {
654  cerr << endl << "ERROR: Incomplete specification for external potential!" << endl << endl;
655  cerr << "Action: specify a y- scale factor (y) for the Gasparini potential." << endl;
656  return 1;
657  }
658 
659  /* Need to specify a z- barrier width scale factor for Gasparini potential */
660  if ( (params["external"].as<string>().find("gasp_prim") != string::npos) &&
661  (!params("empty_width_z")) ) {
662  cerr << endl << "ERROR: Incomplete specification for external potential!" << endl << endl;
663  cerr << "Action: specify a z- scale factor (z) for the Gasparini potential." << endl;
664  return 1;
665  }
666 
667  /* We can only use the hard sphere potential in a 3D system */
668  if ((params["interaction"].as<string>().find("hard_sphere") != string::npos) && (NDIM != 3)) {
669  cerr << endl << "ERROR: Can only use hard sphere potentials for a 3D system!" << endl << endl;
670  cerr << "Action: change the potential or recompile with ndim=3." << endl;
671  return 1;
672  }
673 
674  /* We can only use the hard rod potential in a 1D system */
675  if ((params["interaction"].as<string>().find("hard_rod") != string::npos) && (NDIM != 1)) {
676  cerr << endl << "ERROR: Can only use hard rod potentials for a 1D system!" << endl << endl;
677  cerr << "Action: change the potential or recompile with ndim=1." << endl;
678  return 1;
679  }
680 
681  /* We can only use the delta1D potential in a 1D system */
682  if ((params["interaction"].as<string>().find("delta1D") != string::npos) && (NDIM != 1)) {
683  cerr << endl << "ERROR: Can only use delta1D potentials for a 1D system!" << endl << endl;
684  cerr << "Action: change the potential or recompile with ndim=1." << endl;
685  return 1;
686  }
687 
688  /* If a list of estimators has been supplied, we need to verify */
689  for (string name : params["estimator"].as<vector<string>>()){
690  if (std::find(estimatorName.begin(), estimatorName.end(), name) == estimatorName.end()) {
691  cerr << "ERROR: Tried to measure a non-existent estimator: " << name << endl;
692  cerr << "Action: set estimator to one of:" << endl
693  << "\t[" << estimatorNames<< "]" << endl;
694  return true;
695  }
696  }
697 
698  /* Make sure we don't measure any pigs estimators if we have T > 0 */
699  /* if (!PIGS) { */
700  /* for (string name : params["estimator"].as<vector<string>>()) { */
701  /* bool foundPIGSEstimator = (name.find("pigs") != string::npos); */
702  /* if (foundPIGSEstimator) { */
703  /* cerr << "ERROR: Tried to measure a PIGS estimator when T > 0: " << name << endl; */
704  /* cerr << "Action: remove " << name << " estimator." << endl; */
705  /* return true; */
706  /* } */
707  /* } */
708  /* } */
709  /* /1* Make sure all our estimators are pigs estimators *1/ */
710  /* else { */
711  /* for (string name : params["estimator"].as<vector<string> >()) { */
712  /* bool foundPIGSEstimator = (name.find("pigs") != string::npos) */
713  /* || (name.find("time") != string::npos); */
714 
715  /* if (!foundPIGSEstimator) { */
716  /* cerr << "ERROR: Tried to measure a non-PIGS estimator when T = 0: " << name << endl; */
717  /* cerr << "Action: remove " << name << " estimator." << endl; */
718  /* return true; */
719  /* } */
720  /* } */
721  /* } */
722 
723  /* Only measure cylinder estimators when we have a cylinder cell */
724  if (params["geometry"].as<string>() == "prism") {
725  for (string name : params["estimator"].as<vector<string> >()) {
726  bool foundCylinderEstimator = (name.find("cylinder") != string::npos);
727  if (foundCylinderEstimator) {
728  cerr << "ERROR: Tried to measure a cylinder estimator in a prism: " << name << endl;
729  cerr << "Action: remove " << name << " estimator." << endl;
730  return true;
731  }
732  }
733  }
734 
735  /* Validate the move list */
736  for (string name : params["update"].as<vector<string>>()){
737  if (std::find(moveName.begin(), moveName.end(), name) == moveName.end()) {
738  cerr << "ERROR: Tried to perform a non-existent move: " << name << endl;
739  cerr << "Action: set move to one of:" << endl
740  << "\t[" << moveNames << "]" << endl;
741  return true;
742  }
743  }
744 
745  if (params("validate")) {
746  cerr << "SUCCESS: All command line and/or xml options have been verified." << endl;
747  cerr << "Action: remove --validate flag to proceed with simulation." << endl;
748  return true;
749  }
750 
751  return false;
752 }
753 
754 /**************************************************************************/
761 uint32 Setup::seed (const uint32 startSeed) {
762  return startSeed + params["process"].as<uint32>();
763 }
764 
765 /**************************************************************************/
772 
773  /* Initialize the local box pointer */
774  Container *boxPtr = NULL;
775 
776  /* Setup a cylindrical simulation cell */
777  if (params["geometry"].as<string>() == "cylinder") {
778 
779  double radius = params["radius"].as<double>();
780 
781  /* We need to adjust the maximal possible radius for a hour glass
782  * potential that bows outwards */
783  if ( (params["external"].as<string>() == "hg_tube") &&
784  (params["hourglass_radius"].as<double>() > 0) )
785  radius += params["hourglass_radius"].as<double>();
786 
787 
788  if (definedCell && params("number_particles"))
789  boxPtr = new Cylinder(radius,params["side"].as<dVec>()[NDIM-1]);
790  else if (definedCell && params("density")) {
791  boxPtr = new Cylinder(radius,params["side"].as<dVec>()[NDIM-1]);
792  params.set<int>("number_particles", int(boxPtr->volume*params["density"].as<double>()));
793  }
794  else
795  boxPtr = new Cylinder(params["density"].as<double>(),
796  radius,params["number_particles"].as<int>());
797  }
798  /* Setup a hyperprism */
799  else if (params["geometry"].as<string>() == "prism") {
800 
801  /* we determine if we are using a non-periodic cell for
802  * the graphene potential */
803  iVec periodic;
804  periodic = 1;
805  if (params["external"].as<string>().find("graphene") != string::npos)
806  periodic[NDIM-1] = 0;
807 
808  if (definedCell && params("number_particles"))
809  boxPtr = new Prism(params["side"].as<dVec>(),periodic);
810  else if (definedCell && params("density")) {
811  boxPtr = new Prism(params["side"].as<dVec>(),periodic);
812  params.set<int>("number_particles", int(boxPtr->volume*params["density"].as<double>()));
813  }
814  else
815  boxPtr = new Prism(params["density"].as<double>(),params["number_particles"].as<int>());
816  }
817 
818  /* Add the volume to the parameter map */
819  params.set<double>("volume",boxPtr->volume);
820  params.set<dVec>("side",boxPtr->side);
821 
822  return boxPtr;
823 }
824 
825  /* if ((params["external"].as<string>().find("tube") != string::npos) && (NDIM != 3)) { */
826 /**************************************************************************/
834 
835  int numTimeSlices,numDeltaTau;
836  double tau,imaginaryTimeLength;
837  bool pathBreak = (params["number_broken"].as<int>() > 0 ||
838  params("spatial_subregion"));
839 
840  /* !!NB!! We separate the worldline building into two sections, one for PIGS and
841  * one for PIMC. This can probably be cleaned up in the future. */
842  if (!PIGS) {
843 
844  /* Set the imaginary time length*/
845  if (!params("imaginary_time_length"))
846  params.set<double>("imaginary_time_length",1.0/params["temperature"].as<double>());
847 
848  /* We determine if we have fixed P or tau. We require that the number of time
849  * slices is always even. */
850  if (!params("number_time_slices")) {
851  tau = params["imaginary_time_step"].as<double>();
852  numTimeSlices = static_cast<int>(1.0/(params["temperature"].as<double>() * tau) + EPS);
853  if ((numTimeSlices % 2) != 0)
854  numTimeSlices--;
855  params.set<int>("number_time_slices",numTimeSlices);
856  }
857  else {
858  numTimeSlices = params["number_time_slices"].as<int>();
859  if ((numTimeSlices % 2) != 0)
860  numTimeSlices--;
861  tau = 1.0/(params["temperature"].as<double>() * numTimeSlices);
862  params.set<int>("number_time_slices",numTimeSlices);
863  params.set<double>("imaginary_time_step",tau);
864  }
865  }
866  /* This is for PIGS (T=0) simulations */
867  else {
868  /* Fixing the total imaginary time and the number of time slices */
869  if ( params("imaginary_time_length") && params("number_time_slices") ) {
870 
871  /* We make sure we have an odd number of time slices */
872  numTimeSlices = params["number_time_slices"].as<int>();
873  if ((numTimeSlices % 2) == 0)
874  numTimeSlices++;
875  /* Make sure the center of the path is an odd (even) slice for a
876  closed (open) path calucation */
877  numDeltaTau = (numTimeSlices-1)/2;
878  if ( ((pathBreak)&&(numDeltaTau%2==0)) || ((!pathBreak)&&(numDeltaTau%2==1)) )
879  numTimeSlices += 2;
880 
881  tau = params["imaginary_time_length"].as<double>() / (numTimeSlices-1);
882  params.set<int>("number_time_slices",numTimeSlices);
883  params.set<double>("imaginary_time_step",tau);
884  }
885  /* Fixing the total imaginary time and the time step. We need to make sure
886  * we can get an odd number of slices */
887  else if ( params("imaginary_time_length") && params("imaginary_time_step") ) {
888  tau = params["imaginary_time_step"].as<double>();
889  numTimeSlices = static_cast<int>((params["imaginary_time_length"].as<double>() / tau) + EPS)+1;
890 
891  /* We make sure we have an odd number of time slices */
892  if ((numTimeSlices % 2) == 0)
893  numTimeSlices++;
894 
895  /* Make sure the center of the path is an odd (even) slice for a
896  closed (open) path calucation */
897  numDeltaTau = (numTimeSlices-1)/2;
898  if ( ((pathBreak)&&(numDeltaTau%2==0)) || ((!pathBreak)&&(numDeltaTau%2==1)) )
899  numTimeSlices += 2;
900 
901  imaginaryTimeLength = (numTimeSlices-1)*tau;
902  params.set<double>("imaginary_time_length",imaginaryTimeLength);
903 
904  params.set<int>("number_time_slices",numTimeSlices);
905  }
906  /* Fixing the number of time steps and the size of the time step. */
907  else {
908  /* We make sure we have an odd number of time slices */
909  numTimeSlices = params["number_time_slices"].as<int>();
910  if ((numTimeSlices % 2) == 0)
911  numTimeSlices++;
912  /* Make sure the center of the path is an odd (even) slice for a
913  closed (open) path calucation */
914  numDeltaTau = (numTimeSlices-1)/2;
915  if ( ((pathBreak)&&(numDeltaTau%2==0)) || ((!pathBreak)&&(numDeltaTau%2==1)) )
916  numTimeSlices += 2;
917 
918  /* Set the updated number of time slices and imaginary time length */
919  params.set<int>("number_time_slices",numTimeSlices);
920  imaginaryTimeLength = (numTimeSlices-1)*params["imaginary_time_step"].as<double>();
921  params.set<double>("imaginary_time_length",imaginaryTimeLength);
922  }
923 
924  /* We set the effective temperature to 1.0/imagTimeLength */
925  params.set<double>("temperature",1.0/params["imaginary_time_length"].as<double>());
926  }
927 
928  /* If we haven't fixed Mbar, do so now. We use the value defined in
929  * PRE 74, 036701 (2006). We make sure it is always >= 8 */
930  int Mbar = 0;
931  if (!params("update_length")) {
932  /* Compute the average inter-particle separation */
933  double ell = pow((1.0*params["number_particles"].as<int>() /
934  params["volume"].as<double>()),-1.0/(1.0*NDIM));
935  double lambda = 24.24 / params["mass"].as<double>();
936  Mbar = int(ell*ell/(16*lambda*params["imaginary_time_step"].as<double>()));
937  if (Mbar < 2)
938  Mbar = 2;
939  else if (Mbar > numTimeSlices)
940  Mbar = numTimeSlices/2;
941  params.set<int>("update_length",Mbar);
942  }
943  else
944  {
945  /* Any negative integer sets the maximum Mbar, useful for free particles */
946  if (params["update_length"].as<int>() < 0)
947  Mbar = numTimeSlices - 2;
948  else
949  Mbar = params["update_length"].as<int>();
950  }
951 
952  /* Now we make sure it is even and not too large*/
953  if (Mbar % 2)
954  Mbar++;
955 
956  params.set<int>("update_length",Mbar);
957 
958  if (Mbar > numTimeSlices) {
959  cerr << Mbar << " " << numTimeSlices << endl;
960  cerr << endl << "ERROR: Update length > number time slices!" << endl << endl;
961  cerr << "Action: Increase number_time_slices (P) OR" << endl;
962  cerr << " Decrease update_length (M) OR" << endl;
963  cerr << " Decrease imaginary_time_step (t) OR" << endl;
964  cerr << " Increase imaginary_time_length" << endl;
965  return true;
966  }
967 
968  return false;
969 }
970 
971 /**************************************************************************/
977 
978  /* At present, we need to make sure that if a pair_product action has been
979  * selected, that we turn off the cuttoff by making it the size of the box */
980  if (params["action"].as<string>() == "pair_product" ||
981  !params("potential_cutoff") )
982  params.set<double>("potential_cutoff",params["side"].as<dVec>()[NDIM-1]);
983 
984  /* Set the required constants */
986 
987  /* If we have specified either the center of mass or displace shift on the
988  * command line, we update their values. */
989  if (params("com_delta"))
990  constants()->setCoMDelta(params["com_delta"].as<double>());
991  else {
992  params.set<double>("com_delta",constants()->comDelta());
993  }
994  if (params("displace_delta"))
995  constants()->setDisplaceDelta(params["displace_delta"].as<double>());
996  else {
997  params.set<double>("displace_delta",constants()->displaceDelta());
998  }
999 
1000 }
1001 
1002 /**************************************************************************/
1014 
1015  communicate()->init(params["imaginary_time_step"].as<double>(),
1016  (params["output_config"].as<int>() > 0),params["start_with_state"].as<string>(),
1017  params["fixed"].as<string>());
1018 }
1019 
1020 /*************************************************************************/
1027 
1028  PotentialBase *interactionPotentialPtr = NULL;
1029  if (constants()->intPotentialType() == "free")
1030  interactionPotentialPtr = new FreePotential();
1031  else if (constants()->intPotentialType() == "delta")
1032  interactionPotentialPtr = new DeltaPotential(params["delta_width"].as<double>(),
1033  params["delta_strength"].as<double>());
1034  else if (constants()->intPotentialType() == "sutherland")
1035  interactionPotentialPtr = new SutherlandPotential(params["interaction_strength"].as<double>());
1036  else if (constants()->intPotentialType() == "hard_sphere")
1037  interactionPotentialPtr = new HardSpherePotential(params["scattering_length"].as<double>());
1038  else if (constants()->intPotentialType() == "hard_rod")
1039  interactionPotentialPtr = new HardRodPotential(params["scattering_length"].as<double>());
1040  else if (constants()->intPotentialType() == "delta1D")
1041  interactionPotentialPtr = new Delta1DPotential(params["delta_strength"].as<double>());
1042  else if (constants()->intPotentialType() == "lorentzian")
1043  interactionPotentialPtr = new LorentzianPotential(params["delta_width"].as<double>(),
1044  params["delta_strength"].as<double>());
1045  else if (constants()->intPotentialType() == "aziz")
1046  interactionPotentialPtr = new AzizPotential(boxPtr);
1047  else if (constants()->intPotentialType() == "szalewicz")
1048  interactionPotentialPtr = new SzalewiczPotential(boxPtr);
1049  else if (constants()->intPotentialType() == "harmonic")
1050  interactionPotentialPtr = new HarmonicPotential(params["omega"].as<double>());
1051  else if (constants()->intPotentialType() == "dipole")
1052  interactionPotentialPtr = new DipolePotential();
1053 
1054  return interactionPotentialPtr;
1055 }
1056 
1057 /*************************************************************************/
1064 
1065  PotentialBase *externalPotentialPtr = NULL;
1066 
1067  if (constants()->extPotentialType() == "harmonic")
1068  externalPotentialPtr = new HarmonicPotential(params["omega"].as<double>());
1069  else if (constants()->extPotentialType() == "free")
1070  externalPotentialPtr = new FreePotential();
1071  else if (constants()->extPotentialType() == "osc_tube")
1072  externalPotentialPtr = new HarmonicCylinderPotential(params["radius"].as<double>());
1073  else if (constants()->extPotentialType() == "hg_tube")
1074  externalPotentialPtr = new LJHourGlassPotential(params["radius"].as<double>(),
1075  params["hourglass_radius"].as<double>(), params["hourglass_width"].as<double>());
1076  else if (constants()->extPotentialType() == "plated_lj_tube")
1077  externalPotentialPtr = new PlatedLJCylinderPotential(params["radius"].as<double>(),
1078  params["lj_width"].as<double>(), params["lj_sigma"].as<double>(),
1079  params["lj_epsilon"].as<double>(), params["lj_density"].as<double>());
1080  else if (constants()->extPotentialType() == "lj_tube")
1081  externalPotentialPtr = new LJCylinderPotential(params["radius"].as<double>());
1082  else if (constants()->extPotentialType() == "hard_tube")
1083  externalPotentialPtr = new HardCylinderPotential(params["radius"].as<double>());
1084  else if (constants()->extPotentialType() == "single_well")
1085  externalPotentialPtr = new SingleWellPotential();
1086  else if (constants()->extPotentialType() == "fixed_aziz")
1087  externalPotentialPtr = new FixedAzizPotential(boxPtr);
1088  else if (constants()->extPotentialType() == "gasp_prim")
1089  externalPotentialPtr = new Gasparini_1_Potential(params["empty_width_z"].as<double>(),
1090  params["empty_width_y"].as<double>(),boxPtr);
1091  else if (constants()->extPotentialType() == "graphene")
1092  externalPotentialPtr = new GraphenePotential(params["strain"].as<double>(),
1093  params["poisson"].as<double>(),
1094  params["carbon_carbon_dist"].as<double>(),
1095  params["lj_sigma"].as<double>(),
1096  params["lj_epsilon"].as<double>()
1097  );
1098  else if (constants()->extPotentialType() == "graphenelut")
1099  externalPotentialPtr = new GrapheneLUTPotential(params["strain"].as<double>(),
1100  params["poisson"].as<double>(),
1101  params["carbon_carbon_dist"].as<double>(),
1102  params["lj_sigma"].as<double>(),
1103  params["lj_epsilon"].as<double>(),
1104  boxPtr
1105  );
1106  else if (constants()->extPotentialType() == "fixed_lj")
1107  externalPotentialPtr = new FixedPositionLJPotential(params["lj_sigma"].as<double>(),
1108  params["lj_epsilon"].as<double>(),boxPtr);
1109  else if (constants()->extPotentialType() == "graphenelut3d")
1110  externalPotentialPtr = new GrapheneLUT3DPotential(
1111  params["graphenelut3d_file_prefix"].as<string>(),
1112  boxPtr
1113  );
1114  else if (constants()->extPotentialType() == "graphenelut3dgenerate")
1115  externalPotentialPtr = new GrapheneLUT3DPotentialGenerate(
1116  params["strain"].as<double>(),
1117  params["poisson"].as<double>(),
1118  params["carbon_carbon_dist"].as<double>(),
1119  params["lj_sigma"].as<double>(),
1120  params["lj_epsilon"].as<double>(),
1121  boxPtr
1122  );
1123  else if (constants()->extPotentialType() == "graphenelut3dtobinary")
1124  externalPotentialPtr = new GrapheneLUT3DPotentialToBinary(
1125  params["graphenelut3d_file_prefix"].as<string>(),
1126  boxPtr
1127  );
1128  else if (constants()->extPotentialType() == "graphenelut3dtotext")
1129  externalPotentialPtr = new GrapheneLUT3DPotentialToText(
1130  params["graphenelut3d_file_prefix"].as<string>(),
1131  boxPtr
1132  );
1133 
1134  return externalPotentialPtr;
1135 }
1136 
1137 /*************************************************************************/
1144 
1145  WaveFunctionBase *waveFunctionPtr = NULL;
1146 
1147  if (constants()->waveFunctionType() == "constant")
1148  waveFunctionPtr = new WaveFunctionBase(path,lookup);
1149  else if (constants()->waveFunctionType() == "sech")
1150  waveFunctionPtr = new SechWaveFunction(path,lookup);
1151  else if (constants()->waveFunctionType() == "jastrow")
1152  waveFunctionPtr = new JastrowWaveFunction(path,lookup);
1153  else if (constants()->waveFunctionType() == "lieb")
1154  waveFunctionPtr = new LiebLinigerWaveFunction(path,lookup);
1155  else if (constants()->waveFunctionType() == "sutherland")
1156  waveFunctionPtr = new SutherlandWaveFunction(path,lookup,
1157  params["interaction_strength"].as<double>());
1158 
1159  return waveFunctionPtr;
1160 }
1161 
1162 /*************************************************************************/
1168 ActionBase * Setup::action(const Path &path, LookupTable &lookup,
1169  PotentialBase * externalPotentialPtr,
1170  PotentialBase * interactionPotentialPtr,
1171  WaveFunctionBase * waveFunctionPtr) {
1172 
1173  ActionBase *actionPtr = NULL;
1174 
1175  /* Non Local Actions */
1176  if (constants()->actionType() == "pair_product") {
1177  actionPtr = new NonLocalAction(path,lookup,externalPotentialPtr,
1178  interactionPotentialPtr,waveFunctionPtr,false,constants()->actionType());
1179  }
1180  /* Local Actions */
1181  else {
1182 
1183  /* The factors needed for local actions. */
1184  TinyVector <double,2> VFactor;
1185  TinyVector <double,2> gradVFactor;
1186 
1187  VFactor = 1.0;
1188  gradVFactor = 0.0;
1189  int period = 1;
1190 
1191  if (constants()->actionType() == "gsf") {
1192 
1193  /* There are different pre-factors for even/odd slices, we use the
1194  * Prokof'ev version. I have taken alpha = 0 here. */
1195  VFactor[0] = 2.0/3.0;
1196  VFactor[1] = 4.0/3.0;
1197 
1198  double alpha = 0.0;
1199  gradVFactor[0] = 2.0*alpha/9.0;
1200  gradVFactor[1] = 2.0*(1.0-alpha)/9.0;
1201  period = 2;
1202  }
1203  else if (constants()->actionType() == "li_broughton") {
1204  gradVFactor = 1.0 / 12.0;
1205  period = 2;
1206  }
1207 
1208  /* Do we want to use full staging moves? */
1209  bool local = !params("full_updates");
1210 
1211  actionPtr = new LocalAction(path,lookup,externalPotentialPtr,
1212  interactionPotentialPtr,waveFunctionPtr,VFactor,gradVFactor,local,
1213  constants()->actionType(),constants()->endFactor(),period);
1214  }
1215 
1216  return actionPtr;
1217 }
1218 
1219 /*************************************************************************/
1227 boost::ptr_vector<MoveBase>* Setup::moves(Path &path, ActionBase *actionPtr,
1228  MTRand &random) {
1229 
1230  /* Adapt the move list based on defined parameters */
1231  vector<string> updates = params["update"].as<vector<string>>();
1232  if (PIGS) {
1233 
1234  /* If we have broken paths */
1235  if (params["number_broken"].as<int>() > 0) {
1236 
1237  /* Add the swap break move */
1238  if (std::find(updates.begin(), updates.end(), SwapBreakMove::name)
1239  == updates.end())
1240  updates.push_back(SwapBreakMove::name);
1241 
1242  params.set<vector<string>>("update",updates);
1243  constants()->setAttemptProb("diagonal",0.5);
1244  constants()->setAttemptProb("swap break",0.1);
1245 
1246  } else if (params("spatial_subregion")) {
1247 
1248  /* Add the mid staging move */
1249  if (std::find(updates.begin(), updates.end(), MidStagingMove::name)
1250  == updates.end())
1251  updates.push_back(MidStagingMove::name);
1252 
1253  params.set<vector<string>>("update",updates);
1254  constants()->setAttemptProb("diagonal",0.5);
1255  constants()->setAttemptProb("mid staging",0.1);
1256  }
1257  }
1258  else {
1259 
1260  /* We potentially replace bisection with staging */
1261  if ( params("full_updates") || params("staging") ||
1262  (params["action"].as<string>() == "pair_product") ||
1263  (params["max_winding"].as<int>() > 1) )
1264  {
1265 
1266  /* Replace bisection with staging */
1267  auto it = std::find(updates.begin(), updates.end(), BisectionMove::name);
1268  if (it == updates.end())
1269  updates.push_back(StagingMove::name);
1270  else
1271  updates.at(std::distance(updates.begin(),it)) = StagingMove::name;
1272 
1273  params.set<vector<string>>("update",updates);
1274  }
1275  }
1276 
1277  /* Create a new list of moves that will be returned */
1278  boost::ptr_vector<MoveBase>* movePtr = new boost::ptr_vector<MoveBase>();
1279 
1280  /* Instatiate the moves */
1281  for (auto& name : params["update"].as<vector<string>>())
1282  movePtr->push_back(moveFactory()->Create(name,path,actionPtr,random));
1283 
1284  return movePtr;
1285 }
1286 
1287 /*************************************************************************/
1296 boost::ptr_vector<EstimatorBase> * Setup::estimators(Path &path,
1297  ActionBase *actionPtr, MTRand &random) {
1298 
1299  /* Create the list of estimator pointers */
1300  boost::ptr_vector<EstimatorBase>* estimatorPtr = new boost::ptr_vector<EstimatorBase>();
1301 
1302  /* Instatiate the single path estimators */
1303  for (auto& name : params["estimator"].as<vector<string>>()) {
1304 
1305  if (name.find("multi") == string::npos)
1306  estimatorPtr->push_back(estimatorFactory()->Create(name,path,
1307  actionPtr,random,params["estimator_radius"].as<double>()));
1308  }
1309 
1310  /* We determine where a line break is needed for all estimators writing to
1311  * a common estimator file */
1312  for (const auto & common : {"estimator","cyl_estimator"}) {
1313  auto ePtr = find_if(estimatorPtr->rbegin(), estimatorPtr->rend(),
1314  [common](EstimatorBase &e) { return e.getLabel() == common; });
1315  if (ePtr != estimatorPtr->rend())
1316  ePtr->addEndLine();
1317  }
1318 
1319  return estimatorPtr;
1320 }
1321 
1322 /*************************************************************************/
1329 boost::ptr_vector<EstimatorBase> * Setup::estimators(
1330  boost::ptr_vector<Path> &pathPtrVec,
1331  boost::ptr_vector<ActionBase> &actionPtrVec, MTRand &random) {
1332 
1333  /* Create the list of estimator pointers */
1334  boost::ptr_vector<EstimatorBase>* multiEstimatorPtr = new boost::ptr_vector<EstimatorBase>();
1335 
1336  /* Instatiate the multi path estimators */
1337  for (auto& name : params["estimator"].as<vector<string>>()) {
1338 
1339  if (name.find("multi") != string::npos)
1340  multiEstimatorPtr->push_back(multiEstimatorFactory()->Create(name,pathPtrVec[0],
1341  pathPtrVec[1],&actionPtrVec[0],&actionPtrVec[1],random,
1342  params["estimator_radius"].as<double>()));
1343  }
1344 
1345  return multiEstimatorPtr;
1346 }
1347 
1348 /*************************************************************************/
1354 void Setup::cleanCommandLineOptions(int argc, char*argv[], vector<string> &cmdArg,
1355  vector<string> &cmdSep, vector<string> &cmdVal) {
1356 
1357  /* Get the program name */
1358  cmdArg.push_back(argv[0]);
1359  cmdSep.push_back("");
1360  cmdVal.push_back("");
1361 
1362  for (int n = 1; n < argc; n++) {
1363 
1364  string carg(argv[n]);
1365 
1366  /* Determine if we have a long or short option */
1367  auto numDashes = std::count(carg.begin(), carg.end(), '-');
1368 
1369  /* Short option */
1370  if (numDashes == 1) {
1371  if (carg.length() <= 2) {
1372  cmdArg.push_back(carg);
1373  /* make sure we aren't at the end of the list */
1374  if ((n+1) < argc) {
1375  cmdSep.push_back(" ");
1376  cmdVal.push_back(string(argv[n+1]));
1377  }
1378  else {
1379  cmdSep.push_back("");
1380  cmdVal.push_back("");
1381  }
1382  ++n;
1383  }
1384  else {
1385  cmdArg.push_back(carg.substr(0,2));
1386  cmdSep.push_back(" ");
1387  cmdVal.push_back(carg.substr(2,carg.length()-2));
1388  }
1389  }
1390  /* Long option */
1391  else if (numDashes == 2) {
1392 
1393  /* Do we have an equal sign? */
1394  auto posEqual = carg.find("=");
1395  if (posEqual != string::npos) {
1396  cmdArg.push_back(carg.substr(0,posEqual));
1397  cmdSep.push_back("=");
1398  string nextArg(carg.substr(posEqual+1));
1399 
1400  /* Do we have a space in our value? */
1401  if (nextArg.find(" ") != string::npos)
1402  cmdVal.push_back(string("\"") + nextArg + string("\""));
1403  else
1404  cmdVal.push_back(nextArg);
1405  }
1406  else {
1407  cmdArg.push_back(carg);
1408 
1409  /* make sure we aren't at the end of the list */
1410  if ((n+1) < argc) {
1411  string nextArg(argv[n+1]);
1412 
1413  /* is this option a boolean switch? */
1414  if (nextArg.find('-') == string::npos) {
1415 
1416  cmdSep.push_back(" ");
1417 
1418  /* Do we have a space in our value? */
1419  if (nextArg.find(" ") != string::npos)
1420  cmdVal.push_back(string("\"") + nextArg + string("\""));
1421  else
1422  cmdVal.push_back(nextArg);
1423  ++n;
1424  }
1425  else {
1426  cmdSep.push_back("");
1427  cmdVal.push_back("");
1428  }
1429  }
1430  else {
1431  cmdSep.push_back("");
1432  cmdVal.push_back("");
1433  }
1434  }
1435  }
1436  }
1437 }
1438 
1439 /*************************************************************************/
1451 void Setup::outputOptions(int argc, char *argv[], const uint32 _seed,
1452  const Container *boxPtr, const iVec &nnGrid) {
1453 
1454  /* Pre-process the command line options to make them suitable for log
1455  * output */
1456  vector<string> cmdArg, cmdSep, cmdVal;
1457  cleanCommandLineOptions(argc,argv,cmdArg,cmdSep,cmdVal);
1458 
1459  communicate()->file("log")->stream() << endl << "# ";
1460 
1461  /* Construct the command that would be required to restart the simulation */
1462  bool outputC0 = false;
1463  bool outputmu = false;
1464  bool outputD = false;
1465  bool outputd = false;
1466  bool outputEstimator = false;
1467  bool outputUpdate = false;
1468 
1469  for (uint32 n = 0; n < cmdArg.size(); ++n) {
1470 
1471  auto arg = cmdArg[n];
1472 
1473  if ( (arg == "-s") || (arg == "--start_with_state") ) {
1474  /* do nothing */
1475  }
1476  else if ( (arg == "-C") || (arg == "worm_constant") ) {
1477  communicate()->file("log")->stream() << format("-C %21.15e ") % constants()->C0();
1478  outputC0 = true;
1479  }
1480  else if ( (arg == "-u") || (arg == "chemical_potential") ) {
1481  communicate()->file("log")->stream() << format("-u %21.15e ") % constants()->mu();
1482  outputmu = true;
1483  }
1484  else if ((arg == "-p") || (arg == "--process")) {
1485  communicate()->file("log")->stream() << format("-p %03d ") % params["process"].as<uint32>();
1486  }
1487  else if ((arg == "-D") || (arg == "--com_delta")) {
1488  communicate()->file("log")->stream() << format("-D %21.15e ") % constants()->comDelta();
1489  outputD = true;
1490  }
1491  else if ((arg == "-d") || (arg == "--displace_delta")) {
1492  communicate()->file("log")->stream() << format("-d %21.15e ") % constants()->displaceDelta();
1493  outputd = true;
1494  }
1495  else if ((arg == "--estimator") || (arg == "-e")) {
1496  outputEstimator = true;
1497  }
1498  else if (arg == "--update") {
1499  outputUpdate = false;
1500  }
1501  else
1502  communicate()->file("log")->stream() << cmdArg[n] << cmdSep[n] << cmdVal[n] << " ";
1503  }
1504 
1505  /* Output the restart flag */
1506  communicate()->file("log")->stream() << format("-R %s ") % constants()->id();
1507 
1508  /* If we haven't specified the worm constant, output it now */
1509  if (!outputC0)
1510  communicate()->file("log")->stream() << format("-C %21.15e ") % constants()->C0();
1511 
1512  /* If we haven't specified the chemical potential , output it now */
1513  if (!outputmu)
1514  communicate()->file("log")->stream() << format("-u %21.15e ") % constants()->mu();
1515 
1516  /* If we haven't specified the center of mass Delta, output it now */
1517  if (!outputD)
1518  communicate()->file("log")->stream() << format("-D %21.15e ") % constants()->comDelta();
1519 
1520  /* If we haven't specified the displace delta, output it now */
1521  if (PIGS && !outputd)
1522  communicate()->file("log")->stream() << format("-d %21.15e ") % constants()->displaceDelta();
1523 
1524  /* If we specified estimators, add them to the restart string */
1525  if (outputEstimator) {
1526  for (const auto &estName : params["estimator"].as<vector<string>>()) {
1527  string wrapper("");
1528  if (estName.find(" ") != string::npos)
1529  wrapper = "\"";
1530  communicate()->file("log")->stream() << "--estimator=" << wrapper << estName << wrapper << " ";
1531  }
1532  }
1533 
1534  /* If we specified updates, add them to the restart string */
1535  if (outputUpdate) {
1536  for (const auto &mvName : params["updates"].as<vector<string>>()) {
1537  string wrapper("");
1538  if (mvName.find(" ") != string::npos)
1539  wrapper = "\"";
1540  communicate()->file("log")->stream() << "--update=" << wrapper << mvName << wrapper << " ";
1541  }
1542  }
1543 
1544  communicate()->file("log")->stream() << endl << endl;
1545  communicate()->file("log")->stream() << "---------- Begin Simulation Parameters ----------" << endl;
1546  communicate()->file("log")->stream() << endl;
1547 
1548  /* record the full command line string */
1549  communicate()->file("log")->stream() << format("%-24s\t:\t") % "Command String";
1550 
1551  for (uint32 n = 0; n < cmdArg.size(); n++)
1552  communicate()->file("log")->stream() << cmdArg[n] << cmdSep[n] << cmdVal[n] << " ";
1553  communicate()->file("log")->stream() << endl;
1554 
1555  if (constants()->canonical())
1556  communicate()->file("log")->stream() << format("%-24s\t:\t%s\n") % "Ensemble" % "canonical";
1557  else
1558  communicate()->file("log")->stream() << format("%-24s\t:\t%s\n") % "Ensemble" % "grand canonical";
1559 
1560  if (PIGS)
1561  communicate()->file("log")->stream() << format("%-24s\t:\t%s\n") % "Simulation Type" % "PIGS";
1562  else
1563  communicate()->file("log")->stream() << format("%-24s\t:\t%s\n") % "Simulation Type" % "PIMC";
1564  communicate()->file("log")->stream() << format("%-24s\t:\t%s\n") % "Action Type" % params["action"].as<string>();
1565  communicate()->file("log")->stream() << format("%-24s\t:\t%s\n") % "Number of paths" % params["number_paths"].as<int>();
1566  communicate()->file("log")->stream() << format("%-24s\t:\t%s\n") % "Interaction Potential" %
1567  params["interaction"].as<string>();
1568 
1569  /* Ouptut a possible delta function width and strength */
1570  if ( (params["interaction"].as<string>().find("delta") != string::npos) ||
1571  (params["interaction"].as<string>().find("lorentzian") != string::npos) ) {
1572  communicate()->file("log")->stream() << format("%-24s\t:\t%8.3e\n") % "Delta Width"
1573  % params["delta_width"].as<double>();
1574  communicate()->file("log")->stream() << format("%-24s\t:\t%-7.2f\n") % "Delta Strength"
1575  % params["delta_strength"].as<double>();
1576  }
1577 
1578  /* Ouptut a possible sutherland model interaction strength*/
1579  if (params["interaction"].as<string>().find("sutherland") != string::npos) {
1580  communicate()->file("log")->stream() << format("%-24s\t:\t%8.3e\n") % "Interaction Strength"
1581  % params["interaction_strength"].as<double>();
1582  }
1583 
1584  /* Output harmonic interaction frequecy */
1585  if ( (params["interaction"].as<string>().find("harmonic") != string::npos) ) {
1586  communicate()->file("log")->stream() << format("%-24s\t:\t%-7.2f\n") % "Harmonic Int. Freq."
1587  % params["omega"].as<double>();
1588  }
1589 
1590  /* Output a possible scattering length */
1591  if ( (params["interaction"].as<string>().find("hard_sphere") != string::npos) ||
1592  (params["interaction"].as<string>().find("hard_rod") != string::npos) ) {
1593  communicate()->file("log")->stream() << format("%-24s\t:\t%-7.2f\n") % "Scattering Length"
1594  % params["scattering_length"].as<double>();
1595  }
1596  communicate()->file("log")->stream() << format("%-24s\t:\t%s\n")
1597  % "External Potential" % params["external"].as<string>();
1598 
1599  /* output a possible hourglass radius and width for the hourglass potential. */
1600  if (params["external"].as<string>().find("hg_tube") != string::npos) {
1601  communicate()->file("log")->stream() << format("%-24s\t:\t%-7.2f\n")
1602  % "HourGlass Radius" % params["hourglass_radius"].as<double>();
1603  communicate()->file("log")->stream() << format("%-24s\t:\t%-7.2f\n")
1604  % "HourGlass Width" % params["hourglass_width"].as<double>();
1605  }
1606 
1607  /* output a possible carbon-carbon distance, strain value and LJ
1608  * parameters for the graphene potential */
1609  if (params["external"].as<string>().find("graphene") != string::npos) {
1610  communicate()->file("log")->stream() << format("%-24s\t:\t%-7.2f\n")
1611  % "Carbon Carbon Distance" % params["carbon_carbon_dist"].as<double>();
1612  communicate()->file("log")->stream() << format("%-24s\t:\t%-7.2f\n")
1613  % "Graphene Strain %" % params["strain"].as<double>();
1614  communicate()->file("log")->stream() << format("%-24s\t:\t%-7.2f\n")
1615  % "Graphene Poission Ratio %" % params["poisson"].as<double>();
1616  communicate()->file("log")->stream() << format("%-24s\t:\t%-7.2f\n")
1617  % "Graphene-Carbon LJ Sigma" % params["lj_sigma"].as<double>();
1618  communicate()->file("log")->stream() << format("%-24s\t:\t%-7.2f\n")
1619  % "Graphene-Carbon LJ Epsilon" % params["lj_epsilon"].as<double>();
1620  }
1621 
1622  /* output possible paramters of the plated LJ cylinder potential */
1623  if (params["external"].as<string>().find("plated") != string::npos) {
1624  communicate()->file("log")->stream() << format("%-24s\t:\t%-12.5e\n")
1625  % "Plating Radial Width" % params["lj_width"].as<double>();
1626  communicate()->file("log")->stream() << format("%-24s\t:\t%-12.5e\n")
1627  % "Plating LJ Sigma" % params["lj_sigma"].as<double>();
1628  communicate()->file("log")->stream() << format("%-24s\t:\t%-12.5e\n")
1629  % "Plating LJ Epsilon" % params["lj_epsilon"].as<double>();
1630  communicate()->file("log")->stream() << format("%-24s\t:\t%-12.5e\n")
1631  % "Plating LJ Density" % params["lj_density"].as<double>();
1632  }
1633 
1634  if (PIGS) {
1635  communicate()->file("log")->stream() << format("%-24s\t:\t%s\n")
1636  % "Wavefunction Type" % params["wavefunction"].as<string>();
1637  communicate()->file("log")->stream() <<
1638  format("%-24s\t:\t%7.5f\n") % "End Factor" % params["end_factor"].as<double>();
1639  /* Output possible wave function parameters */
1640  if ( (params["wavefunction"].as<string>().find("lieb") != string::npos) ) {
1641  communicate()->file("log")->stream() << format("%-24s\t:\t%-7.2f\n") % "Wavefunction length scale"
1642  % params["R_LL_wfn"].as<double>();
1643  communicate()->file("log")->stream() << format("%-24s\t:\t%-7.4f\n") % "Wavefunction wave number"
1644  % params["k_LL_wfn"].as<double>();
1645  }
1646  }
1647  communicate()->file("log")->stream() <<
1648  format("%-24s\t:\t%7.5f\n") % "Temperature" % params["temperature"].as<double>();
1649  communicate()->file("log")->stream() <<
1650  format("%-24s\t:\t%7.5f\n") % "Chemical Potential" % constants()->mu();
1651  communicate()->file("log")->stream() <<
1652  format("%-24s\t:\t%7.5f\n") % "Particle Mass" % params["mass"].as<double>();
1653  communicate()->file("log")->stream() <<
1654  format("%-24s\t:\t%d\n") % "Number Time Slices" % constants()->numTimeSlices();
1655  communicate()->file("log")->stream() <<
1656  format("%-24s\t:\t%7.5f\n") % "Specified Imaginary Time Step" % params["imaginary_time_step"].as<double>();
1657  communicate()->file("log")->stream() <<
1658  format("%-24s\t:\t%7.5f\n") % "Imaginary Time Step" % constants()->tau();
1659  communicate()->file("log")->stream() <<
1660  format("%-24s\t:\t%7.5f\n") % "Imaginary Time Length" % constants()->imagTimeLength();
1661  communicate()->file("log")->stream() <<
1662  format("%-24s\t:\t%d\n") % "Initial Number Particles" % params["number_particles"].as<int>();
1663  communicate()->file("log")->stream() <<
1664  format("%-24s\t:\t%7.5f\n") % "Initial Density"
1665  % (1.0*params["number_particles"].as<int>()/boxPtr->volume);
1666  communicate()->file("log")->stream() <<
1667  format("%-24s\t:\t%d\n") % "Num. Broken World-lines" % params["number_broken"].as<int>();
1668  if ( constants()->spatialSubregionOn()){
1669  communicate()->file("log")->stream() <<
1670  format("%-24s\t:\t%d\n") % "Spatial Subregion" % params["spatial_subregion"].as<double>();
1671  }
1672  communicate()->file("log")->stream() <<
1673  format("%-24s\t:\t%s\n") % "Container Type" % boxPtr->name;
1674  communicate()->file("log")->stream() << format("%-24s\t:\t") % "Container Dimensions";
1675  for (int i = 0; i < NDIM; i++) {
1676  communicate()->file("log")->stream() << format("%7.5f") % boxPtr->side[i];
1677  if (i < (NDIM-1))
1678  communicate()->file("log")->stream() << " x ";
1679  else
1680  communicate()->file("log")->stream() << endl;
1681  }
1682  communicate()->file("log")->stream() << format("%-24s\t:\t%7.5f\n") % "Container Volume"
1683  % boxPtr->volume;
1684  communicate()->file("log")->stream() << format("%-24s\t:\t") % "Lookup Table";
1685  for (int i = 0; i < NDIM; i++) {
1686  communicate()->file("log")->stream() << format("%d") % nnGrid[i];
1687  if (i < (NDIM-1))
1688  communicate()->file("log")->stream() << " x ";
1689  else
1690  communicate()->file("log")->stream() << endl;
1691  }
1692  communicate()->file("log")->stream() <<
1693  format("%-24s\t:\t%d\n") % "Maximum Winding Sector" % params["max_winding"].as<int>();
1694  communicate()->file("log")->stream() << format("%-24s\t:\t%7.5f\n") % "Initial Worm Constant" %
1695  params["worm_constant"].as<double>();
1696  communicate()->file("log")->stream() << format("%-24s\t:\t%7.5f\n") % "Worm Constant" % constants()->C0();
1697  communicate()->file("log")->stream() << format("%-24s\t:\t%7.5f\n") % "Initial CoM Delta" % params["com_delta"].as<double>();
1698  communicate()->file("log")->stream() << format("%-24s\t:\t%7.5f\n") % "CoM Delta" % constants()->comDelta();
1699  if (PIGS) {
1700  communicate()->file("log")->stream() << format("%-24s\t:\t%7.5f\n") % "Initial Displace Delta" % params["displace_delta"].as<double>();
1701  communicate()->file("log")->stream() << format("%-24s\t:\t%7.5f\n") % "Displace Delta" % constants()->displaceDelta();
1702  }
1703  communicate()->file("log")->stream() <<
1704  format("%-24s\t:\t%d\n") % "Bisection Parameter" % constants()->b();
1705  communicate()->file("log")->stream() <<
1706  format("%-24s\t:\t%d\n") % "Update Length" % constants()->Mbar();
1707  communicate()->file("log")->stream() <<
1708  format("%-24s\t:\t%7.5f\n") % "Potential Cutoff Length" % params["potential_cutoff"].as<double>();
1709  communicate()->file("log")->stream() <<
1710  format("%-24s\t:\t%d\n") % "Bin Size" % params["bin_size"].as<int>();
1711  communicate()->file("log")->stream() <<
1712  format("%-24s\t:\t%d\n") % "Number EQ Steps" % params["number_eq_steps"].as<uint32>();
1713  communicate()->file("log")->stream() <<
1714  format("%-24s\t:\t%d\n") % "Number Bins Stored" % params["number_bins_stored"].as<int>();
1715  communicate()->file("log")->stream() << format("%-24s\t:\t%d\n") % "Random Number Seed" % _seed;
1716  if (params["virial_window"].as<int>() == 0){
1717  communicate()->file("log")->stream() <<
1718  format("%-24s\t:\t%d\n") % "Virial Window" % constants()->numTimeSlices();
1719  }
1720  else{
1721  communicate()->file("log")->stream() <<
1722  format("%-24s\t:\t%d\n") % "Virial Window" % params["virial_window"].as<int>();
1723  }
1724  communicate()->file("log")->stream() << endl;
1725  communicate()->file("log")->stream() << "---------- End Simulation Parameters ------------" << endl;
1726 }
Action class definitions.
Holds a base class that all action classes will be derived from.
Definition: action.h:29
Computes the value of the semi-empircal Aziz potential that is known to be accurate for He-4.
Definition: potential.h:678
File * file(string type)
Get method returning file object.
Definition: communicator.h:85
void init(double, bool, string, string)
Initialize the output files.
int numTimeSlices()
Get number of time slices.
Definition: constants.h:99
double mu() const
Get chemical potential.
Definition: constants.h:43
void setCoMDelta(double _comDelta)
Set the CoM move size.
Definition: constants.h:122
double displaceDelta() const
Get center of mass shift.
Definition: constants.h:54
double imagTimeLength() const
Get the extent in imaginary time.
Definition: constants.h:42
int Mbar()
Get Mbar.
Definition: constants.h:97
string id()
Get simulation UUID.
Definition: constants.h:102
void initConstants(po::variables_map &)
Initialize all constants from command line, XMl and defaults.
Definition: constants.cpp:32
double C0() const
Get worm factor C0.
Definition: constants.h:49
double tau() const
Get imaginary time step.
Definition: constants.h:44
int b()
Get bisection level.
Definition: constants.h:98
double comDelta() const
Get center of mass shift.
Definition: constants.h:53
string actionType() const
Get wave action type.
Definition: constants.h:113
void setDisplaceDelta(double _displaceDelta)
Set the displace move size.
Definition: constants.h:123
The base class which holds details on the generalized box that our system will be simulated inside of...
Definition: container.h:24
double volume
The volume of the container in A^3.
Definition: container.h:35
string name
The name of the container.
Definition: container.h:39
dVec side
The linear dimensions of the box.
Definition: container.h:31
A three dimensional cylinder with fixed boundary conditions in the x and y directions and periodic bo...
Definition: container.h:143
Computes the effective potential from the exact two-body density matrix for delta interactions in 1D.
Definition: potential.h:1060
Computes the potential energy for delta function interaction potential, approximated here as the limi...
Definition: potential.h:239
Computes the potential energy for polarized electrical dipoles with strength D in reduced units where...
Definition: potential.h:360
The base class that all estimator classes will be derived from.
Definition: estimator.h:28
Computes the potential energy resulting from a series of fixed helium atoms that are not updated and ...
Definition: potential.h:924
Returns Lennard-Jones potential between adatoms and fixed postions in FILENAME.
Definition: potential.h:961
Free potential.
Definition: potential.h:123
Computes potential energy for Gasparini potential.
Definition: potential.h:990
FIXME Returns van der Waals' potential between a helium adatom and a graphene sheet using summation i...
Definition: potential.h:1422
FIXME Returns van der Waals' potential between a helium adatom and a graphene sheet using summation i...
Definition: potential.h:1624
FIXME Returns van der Waals' potential between a helium adatom and a graphene sheet using summation i...
Definition: potential.h:1650
Returns van der Waals' potential between a helium adatom and a graphene sheet using summation in reci...
Definition: potential.h:1296
Returns van der Waals' potential between a helium adatom and a graphene sheet using summation in reci...
Definition: potential.h:1222
The smooth non-corregated version of the helium-carbon nanotube potential.
Definition: potential.h:1174
Computes the value of the external wall potential for a hard-walled cylindrical cavity.
Definition: potential.h:405
Computes the effective potential from the exact two-body density matrix for hard rods in 1D.
Definition: potential.h:1102
Computes the effective potential from the exact two-body density matrix for hard spheres in 3D.
Definition: potential.h:1032
Computes the potential energy for an external harmonic potential with axial symmetry.
Definition: potential.h:202
Computes the potential energy for an external harmonic potential.
Definition: potential.h:145
Implementation of a Jastrow trial wave function suitable for He.
Definition: wavefunction.h:78
Computes the value of the external wall potential for a cylindrical cavity.
Definition: potential.h:535
Computes the value of the external wall potential for an hour-glass shaped cavity.
Definition: potential.h:622
Implementation of a Jastrow trial wave function suitable for He.
Definition: wavefunction.h:107
A base class to be inherited by actions that are local in imaginary time.
Definition: action.h:150
The particle (bead) lookup table.
Definition: lookuptable.h:29
Computes the potential energy for delta function interaction potential, approximated here as the limi...
Definition: potential.h:275
A base class to be inherited by actions that are non-local in imaginary time.
Definition: action.h:257
void print()
print out the parameter map
Definition: setup.cpp:75
void set(const string &, const Ttype)
Set a parameter from a value.
Definition: setup.h:192
void update(int, char *[], po::options_description &)
Update parameters from command line.
Definition: setup.cpp:151
void add(string, string, string)
Add a parameter to the map without a default value.
Definition: setup.h:125
void setupCommandLine(boost::ptr_map< string, po::options_description > &)
Insert options into the options_description data structure to be read from the command line.
Definition: setup.cpp:97
The space-time trajectories.
Definition: path.h:29
Computes the value of the external wall potential for a plated cylindrical cavity.
Definition: potential.h:440
The base class from which all specific potentials are derived from.
Definition: potential.h:32
A NDIM-dimensional hyperprism with periodic boundary conditions.
Definition: container.h:106
Implementation of the Psi_T = sech(a*x) trial wave function suitable for the simple harmonic osscilat...
Definition: wavefunction.h:58
ActionBase * action(const Path &, LookupTable &, PotentialBase *, PotentialBase *, WaveFunctionBase *)
Setup the action.
Definition: setup.cpp:1168
Setup()
Setup the program_options variables.
Definition: setup.cpp:262
uint32 seed(const uint32)
Return the random seed.
Definition: setup.cpp:761
bool parseOptions()
Parse the command line options for obvious errors and return values.
Definition: setup.cpp:488
PotentialBase * interactionPotential(const Container *)
Setup the interaction potential.
Definition: setup.cpp:1026
WaveFunctionBase * waveFunction(const Path &, LookupTable &)
Setup the trial wave function.
Definition: setup.cpp:1143
void setConstants()
Setup the simulation constants.
Definition: setup.cpp:976
void outputOptions(int, char *[], const uint32, const Container *, const iVec &)
Output the simulation parameters to a log file.
Definition: setup.cpp:1451
void communicator()
Setup the communicator.
Definition: setup.cpp:1013
boost::ptr_vector< EstimatorBase > * estimators(Path &, ActionBase *, MTRand &)
Create a list of estimators to be measured.
Definition: setup.cpp:1296
void getOptions(int, char *[])
Define all command line options and get them from the command line.
Definition: setup.cpp:464
boost::ptr_vector< MoveBase > * moves(Path &, ActionBase *, MTRand &)
Define the Monte Carlo updates that will be performed.
Definition: setup.cpp:1227
Parameters params
All simulation parameters.
Definition: setup.h:302
Container * cell()
Setup the simulation cell.
Definition: setup.cpp:771
bool worldlines()
Setup the worldlines.
Definition: setup.cpp:833
PotentialBase * externalPotential(const Container *)
Setup the external potential.
Definition: setup.cpp:1063
Computes the potential energy for an external single well potential.
Definition: potential.h:175
Computes the potential energy for the periodic Sutherland model which approximates long-range 1/r^2 i...
Definition: potential.h:313
Implementation of the Sutherland model exact wavefunction.
Definition: wavefunction.h:136
Computes the value of the semi-empircal Szalewicz potential that is known to be accurate for He-4.
Definition: potential.h:774
Holds a base class that all trial wave function classes will be derived from.
Definition: wavefunction.h:26
#define NDIM
Number of spatial dimnsions.
Definition: common.h:71
unsigned long uint32
Unsigned integer type, at least 32 bits.
Definition: common.h:105
#define EPS
A small number.
Definition: common.h:94
TinyVector< double, NDIM > dVec
A NDIM-vector of type double.
Definition: common.h:111
TinyVector< int, NDIM > iVec
A NDIM-vector of type integer.
Definition: common.h:114
Class definitions for all file input/output.
Communicator * communicate()
Global public access to the communcator singleton.
Definition: communicator.h:121
ConstantParameters class definition.
ConstantParameters * constants()
Global public access to the constants.
Definition: constants.h:201
The simulation cell.
Estimator class definitions.
Move class definitions.
All possible potential classes.
string getList(const vector< string > &options)
Create a comma separated list from a vector of strings.
Definition: setup.cpp:25
ostream & operator<<(ostream &os, const vector< string > &vals)
Overload << to print a vector of strings to the terminal.
Definition: setup.cpp:37
Parses command line input and sets up the details of the simulation.
Action class definitions.