Path Integral Quantum Monte Carlo
pdrive.cpp
Go to the documentation of this file.
1 
8 #include "common.h"
9 #include "constants.h"
10 #include "container.h"
11 #include "path.h"
12 #include "potential.h"
13 #include "action.h"
14 #include "wavefunction.h"
15 #include "pimc.h"
16 #include "lookuptable.h"
17 #include "communicator.h"
18 #include "setup.h"
19 #include "cmc.h"
20 #include "move.h"
21 
30 int main (int argc, char *argv[]) {
31 
32  /* Get initial time */
33  time_t start_time = time(NULL);
34  time_t current_time; //current time
35  bool wallClockReached = false;
36 
37  uint32 seed = 139853; // The seed for the random number generator
38 
39  Setup setup;
40 
41  /* Attempt to parse the command line options */
42  try {
43  setup.getOptions(argc,argv);
44  }
45  catch(exception& ex) {
46  cerr << "error: " << ex.what() << "\n";
47  return 1;
48  }
49  catch(...) {
50  cerr << "Exception of unknown type!\n";
51  }
52 
53  /* Parse the setup options and possibly exit */
54  if (setup.parseOptions())
55  return 1;
56 
57  /* The global random number generator, we add the process number to the seed (for
58  * use in parallel simulations.*/
59  seed = setup.seed(seed);
60  MTRand random(seed);
61 
62  /* Get the simulation box */
63  Container *boxPtr = setup.cell();
64 
65  /* Create the worldlines */
66  if (setup.worldlines())
67  return 1;
68 
69  /* Setup the simulation constants */
70  setup.setConstants();
71 
72  /* Setup the simulation communicator */
73  setup.communicator();
74 
75  /* Get number of paths to use */
76  int Npaths = constants()->Npaths();
77 
78  /* Create and initialize the Nearest Neighbor Lookup Table */
79  boost::ptr_vector<LookupTable> lookupPtrVec;
80  for(int i=0; i<Npaths; i++){
81  lookupPtrVec.push_back(
82  new LookupTable(boxPtr,constants()->numTimeSlices(),
83  constants()->initialNumParticles()));
84  }
85 
86  /* Create and initialize the potential pointers */
87  PotentialBase *interactionPotentialPtr = setup.interactionPotential(boxPtr);
88  PotentialBase *externalPotentialPtr = setup.externalPotential(boxPtr);
89  if ((constants()->extPotentialType() == "graphenelut3dtobinary") ||
90  (constants()->extPotentialType() == "graphenelut3dtotext") ||
91  (constants()->extPotentialType() == "graphenelut3dgenerate") ) {
92  return 99;
93  }
94 
95  /* Get the initial conditions associated with the external potential */
96  /* Must use the copy constructor as we return a copy */
97  Array<dVec,1> initialPos =
98  externalPotentialPtr->initialConfig(boxPtr,random,constants()->initialNumParticles());
99 
100  /* Perform a classical canonical pre-equilibration to obtain a suitable
101  * initial state if N > 0*/
102  if (!constants()->restart() && (constants()->initialNumParticles() > 0) ) {
103  ClassicalMonteCarlo CMC(externalPotentialPtr,interactionPotentialPtr,random,boxPtr,
104  initialPos);
105  CMC.run(constants()->numEqSteps(),0);
106  }
107 
108  /* Setup the path data variable */
109  boost::ptr_vector<Path> pathPtrVec;
110  for(int i=0; i<Npaths; i++){
111  pathPtrVec.push_back(
112  new Path(boxPtr,lookupPtrVec[i],constants()->numTimeSlices(),
113  initialPos,constants()->numBroken()));
114  }
115 
116  /* The Trial Wave Function (constant for pimc) */
117  WaveFunctionBase *waveFunctionPtr = setup.waveFunction(pathPtrVec.front(),lookupPtrVec.front());
118 
119  /* Setup the action */
120  boost::ptr_vector<ActionBase> actionPtrVec;
121  for(int i=0; i<Npaths; i++){
122  actionPtrVec.push_back(
123  setup.action(pathPtrVec[i],lookupPtrVec[i],externalPotentialPtr,
124  interactionPotentialPtr,waveFunctionPtr) );
125  }
126 
127  /* The list of Monte Carlo updates (moves) that will be performed */
128  boost::ptr_vector< boost::ptr_vector<MoveBase> > movesPtrVec;
129  for(int i=0; i<Npaths;i++){
130  movesPtrVec.push_back(
131  setup.moves(pathPtrVec[i],&actionPtrVec[i],random));
132  }
133 
134  /* The list of estimators that will be performed */
135  boost::ptr_vector< boost::ptr_vector<EstimatorBase> > estimatorsPtrVec;
136  for(int i=0; i<Npaths;i++){
137  estimatorsPtrVec.push_back(
138  setup.estimators(pathPtrVec[i],&actionPtrVec[i],random));
139 
140  /* Add labels to estimator output files for multiple paths */
141  if(i > 0) {
142  for(uint32 j = 0; j < estimatorsPtrVec.back().size(); j++)
143  estimatorsPtrVec.back().at(j).appendLabel(str(format("%d") % (i+1)));
144  }
145  }
146 
147  /* Setup the multi-path estimators */
148  if(Npaths>1){
149  estimatorsPtrVec.push_back(setup.estimators(pathPtrVec,actionPtrVec,random));
150  }
151 
152  /* Setup the pimc object */
153  PathIntegralMonteCarlo pimc(pathPtrVec,random,movesPtrVec,estimatorsPtrVec,
154  !setup.params["start_with_state"].as<string>().empty(),
155  setup.params["bin_size"].as<int>());
156 
157  /* A silly banner */
158  if (PIGS)
159  cout << endl
160  << " _____ _____ _____ _____" << endl
161  << "| __ \\ |_ _| / ____| / ____|" << endl
162  << "| |__) | | | | | __ | (___" << endl
163  << "| ___/ | | | | |_ | \\___ \\" << endl
164  << "| | _| |_ | |__| | ____) |" << endl
165  << "|_| |_____| \\_____| |_____/" << endl
166  << endl;
167  else
168  cout << endl
169  << " _____ _____ __ __ _____" << endl
170  << " | __ \\ |_ _| | \\/ | / ____|" << endl
171  << " | |__) | | | | \\ / | | | " << endl
172  << " | ___/ | | | |\\/| | | | " << endl
173  << " | | _| |_ | | | | | |____ " << endl
174  << " |_| |_____| |_| |_| \\_____|" << endl
175  << endl;
176 
177  /* If this is a fresh run, we equilibrate and output simulation parameters to disk */
178  if (!constants()->restart()) {
179 
180  /* Equilibrate */
181  cout << format("[PIMCID: %s] - Pre-Equilibration Stage.") % constants()->id() << endl;
182  for (uint32 n = 0; n < constants()->numEqSteps(); n++)
183  pimc.equilStep(n,setup.params("relax"),setup.params("relaxmu"));
184 
185  /* If we have relaxed the chemical potential, need to update file names
186  * in the grand canonical ensemble */
187  if (!setup.params("canonical") && setup.params("relaxmu")) {
189  }
190 
191  /* Output simulation details/parameters */
192  setup.outputOptions(argc,argv,seed,boxPtr,lookupPtrVec.front().getNumNNGrid());
193  }
194 
195  cout << format("[PIMCID: %s] - Measurement Stage.") % constants()->id() << endl;
196 
197  /* Sample */
198  int oldNumStored = 0;
199  int outNum = 0;
200  int numOutput = setup.params["output_config"].as<int>();
201  uint32 n = 0;
202  do {
203  pimc.step();
204  if (pimc.numStoredBins > oldNumStored) {
205  oldNumStored = pimc.numStoredBins;
206  cout << format("[PIMCID: %s] - Bin #%5d stored to disk.") % constants()->id()
207  % oldNumStored << endl;
208  }
209  n++;
210 
211  /* Output configurations to disk */
212  if ((numOutput > 0) && ((n % numOutput) == 0)) {
213  pathPtrVec.front().outputConfig(outNum);
214  outNum++;
215  }
216 
217  /* Check if we've reached the wall clock limit*/
218  if(constants()->wallClockOn()){
219  current_time = time(NULL);
220  if ( uint32(current_time) > (uint32(start_time) + constants()->wallClock()) ){
221  wallClockReached = true;
222  break;
223  }
224  }
225  } while (pimc.numStoredBins < setup.params["number_bins_stored"].as<int>());
226  if (wallClockReached)
227  cout << format("[PIMCID: %s] - Wall clock limit reached.") % constants()->id() << endl;
228  else
229  cout << format("[PIMCID: %s] - Measurement complete.") % constants()->id() << endl;
230 
231  /* Output Results */
232  if (!constants()->saveStateFiles())
233  pimc.saveState(1);
234  pimc.finalOutput();
235 
236  /* Free up memory */
237  delete interactionPotentialPtr;
238  delete externalPotentialPtr;
239  delete boxPtr;
240  delete waveFunctionPtr;
241 
242  initialPos.free();
243 
244  return 1;
245 }
Action class definitions.
Pre-equilibration via classical Monte Carlo.
Definition: cmc.h:26
void run(uint32, bool)
Perform the Monte Carlo equilibration.
Definition: cmc.cpp:88
void updateNames()
Update the data name and rename any existing files.
string id()
Get simulation UUID.
Definition: constants.h:102
uint32 numEqSteps()
Get the number of equilibration steps.
Definition: constants.h:103
The base class which holds details on the generalized box that our system will be simulated inside of...
Definition: container.h:24
The particle (bead) lookup table.
Definition: lookuptable.h:29
The main driver class for the entire path integral monte carlo program.
Definition: pimc.h:37
void saveState(const int finalSave=0)
Save the state of the simulation to disk, including all path, worm, move and estimator data.
Definition: pimc.cpp:804
void equilStep(const uint32, const bool, const bool)
Equilibration.
Definition: pimc.cpp:561
void finalOutput()
Output simulation statistics to disk.
Definition: pimc.cpp:752
void step()
PIMC step.
Definition: pimc.cpp:688
int numStoredBins
Number of stored estimators.
Definition: pimc.h:64
The space-time trajectories.
Definition: path.h:29
The base class from which all specific potentials are derived from.
Definition: potential.h:32
virtual Array< dVec, 1 > initialConfig(const Container *, MTRand &, const int)
Default Initial configuration of particles.
Definition: potential.cpp:41
Setup the simulation.
Definition: setup.h:255
ActionBase * action(const Path &, LookupTable &, PotentialBase *, PotentialBase *, WaveFunctionBase *)
Setup the action.
Definition: setup.cpp:1168
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
Holds a base class that all trial wave function classes will be derived from.
Definition: wavefunction.h:26
ClassicalMonteCarlo class definition.
Global common header with shared dependencies and methods.
unsigned long uint32
Unsigned integer type, at least 32 bits.
Definition: common.h:105
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.
LookupTable class definition.
Move class definitions.
Path class definition.
int main(int argc, char *argv[])
Main driver.
Definition: pdrive.cpp:30
PathIntegralMonteCarlo class definition.
All possible potential classes.
Parses command line input and sets up the details of the simulation.
Action class definitions.