Path Integral Quantum Monte Carlo
Public Member Functions | Data Fields | Protected Member Functions | Protected Attributes | Static Protected Attributes | Friends
MoveBase Class Referenceabstract

The base class that all moves will be derived from. More...

#include <move.h>

+ Inheritance diagram for MoveBase:
+ Collaboration diagram for MoveBase:

Public Member Functions

 MoveBase (Path &, ActionBase *, MTRand &, ensemble _operateOnConfig=ANY, bool _varLength=false)
 Move naming conventions: More...
 
virtual ~MoveBase ()
 Destructor.
 
virtual string getName ()
 return the move name
 
double getAcceptanceRatio ()
 Get the acceptance ratio.
 
double getTotAcceptanceRatio ()
 Get the total acceptance ratio.
 
double getAcceptanceRatioLevel (int n)
 Get the acceptance ratio by level.
 
int getNumAttempted ()
 Get the number of moves attempted.
 
int getNumAccepted ()
 Get the number of moves accepted.
 
int getNumAttemptedLevel (int n)
 Get the number of moves attempted by level.
 
int getNumAcceptedLevel (int n)
 Get the number of moves accepted by level.
 
virtual bool attemptMove ()=0
 Attempt the move (will be overloaded).
 
void resetTotAccept ()
 Reset the total accepted counter.
 
void resetAccept ()
 Reset the number accepted counter.
 

Data Fields

ensemble operateOnConfig
 What configurations do we operate on?
 
bool variableLength
 Does the move have a variable length?
 
string name1
 

Protected Member Functions

virtual void keepMove ()
 Keep the move. More...
 
virtual void undoMove ()=0
 undo the move
 
dVec newStagingPosition (const beadLocator &, const beadLocator &, const int, const int)
 Returns a new staging position which will exactly sample the kinetic action. More...
 
dVec newStagingPosition (const beadLocator &, const beadLocator &, const int, const int, iVec &)
 Returns a new staging position which will exactly sample the kinetic action in different winding sectors. More...
 
iVec sampleWindingSector (const beadLocator &, const beadLocator &, const int, double &)
 Obtain a winding sector for a stage-like move. More...
 
iVec getWindingNumber (const beadLocator &, const beadLocator &)
 Find the winding number for a path between two beads. More...
 
dVec newFreeParticlePosition (const beadLocator &)
 Generates a new position, which exactly samples the free particle density matrix. More...
 
dVec newBisectionPosition (const beadLocator &, const int)
 Returns a new bisection position which will exactly sample the kinetic action. More...
 
void printMoveState (string)
 
void checkMove (int, double)
 

Protected Attributes

Pathpath
 A reference to the paths.
 
ActionBaseactionPtr
 A base pointer to the action.
 
MTRand & random
 A reference to the RNG.
 
bool success
 Did we sucessfully perform a move?
 
uint32 numAccepted
 The number of accepted moves.
 
uint32 numAttempted
 The number of attempted moves.
 
int numToMove
 The number of particles moved.
 
int numLevels
 
Array< uint32, 1 > numAcceptedLevel
 The number of moves accepted at each level.
 
Array< uint32, 1 > numAttemptedLevel
 The number of moves attempted at each level.
 
Array< dVec, 1 > originalPos
 The original particle positions.
 
Array< dVec, 1 > newPos
 New particle positions.
 
vector< iVecwinding
 The winding vectors

 
vector< int > windingSector
 Used to index different winding sectors.
 
vector< double > cumrho0
 Used for tower-sampling winding sectors.
 
int maxWind
 The largest winding number.
 
int numWind
 The total number of winding vectors.
 
double oldAction
 The original potential action.
 
double newAction
 The new potential action.
 
double deltaAction
 The action difference.
 
double sqrt2LambdaTau
 sqrt(2 * Lambda * tau)
 
double sqrtLambdaTau
 sqrt(Lambda * tau)
 
beadLocator nBeadIndex
 Neighbor bead index.
 
dVec neighborPos
 Staging neighbor position.
 
dVec newRanPos
 Staing random position.
 
double newK
 
double oldK
 The old and new kinetic action.
 
double newV
 
double oldV
 The old and new potential action.
 

Static Protected Attributes

static uint32 totAccepted = 0
 The total number of moves accepted.
 
static uint32 totAttempted = 0
 The total number of moves attempted.
 

Friends

class PathIntegralMonteCarlo
 

Detailed Description

The base class that all moves will be derived from.

There will be a bunch of different types of moves, those which move individual particles, the center of mass of an entire worldine loop etc. They will all share the basic functionality defined here.

Definition at line 30 of file move.h.

Constructor & Destructor Documentation

◆ MoveBase()

MoveBase::MoveBase ( Path _path,
ActionBase _actionPtr,
MTRand &  _random,
ensemble  _operateOnConfig = ANY,
bool  _varLength = false 
)

Move naming conventions:

1) be as descriptive as possible 2) spaces are fine Constructor.

Definition at line 60 of file move.cpp.

61  :
62  operateOnConfig(_operateOnConfig),
63  variableLength(_varLength),
64  path(_path),
65  actionPtr(_actionPtr),
66  random(_random) {
67 
68  /* Initialize private data to zero */
70  success = false;
71 
72  /* Setup the move attempt/accept arrays */
73  int b = int (ceil(log(1.0*constants()->Mbar()) / log(2.0)-EPS));
74  numAcceptedLevel.resize(b+1);
75  numAttemptedLevel.resize(b+1);
76 
77  numAcceptedLevel = 0;
79 
80  sqrtLambdaTau = sqrt(constants()->lambda() * constants()->tau());
81  sqrt2LambdaTau = sqrt(2.0)*sqrtLambdaTau;
82 
83  /* Setup the free density matrix arrays for sampling different
84  * winding sectors. We will sample w = -maxWind ... maxWind */
85  maxWind = constants()->maxWind();
86  numWind = ipow(2*maxWind + 1,NDIM);
87 
88  /* initialize the cumulative probability distribution */
89  cumrho0.resize(numWind);
90 
91  /* For each integer labelling a winding sector, we construct the winding
92  * vector and append to a matrix */
93  iVec wind;
94  for (int n = 0; n < numWind; n++ ) {
95  for (int i = 0; i < NDIM; i++) {
96  int scale = 1;
97  for (int j = i+1; j < NDIM; j++)
98  scale *= (2*maxWind + 1);
99  wind[i] = (n/scale) % (2*maxWind + 1);
100 
101  /* Wrap into the appropriate winding sector */
102  wind[i] -= (wind[i] > maxWind)*(2*maxWind + 1);
103  }
104 
105  /* Adjust for any non-periodic boundary conditions */
106  for (int i = 0; i < NDIM; i++) {
107  if (!path.boxPtr->periodic[i])
108  wind[i] = 0;
109  }
110 
111  /* Store the winding number */
112  winding.push_back(wind);
113  }
114 
115  /* Now we would like to sort the winding number array for the effecient
116  * calculation of maximal probabilities. We sort based on the winding
117  * sector. */
118  std::stable_sort(winding.begin(), winding.end(), [](const iVec& w1, const iVec& w2) {
119  return blitz::max(abs(w1)) < blitz::max(abs(w2));
120  /* return (blitz::dot(w1,w1) < blitz::dot(w2,w2)); */
121  });
122 
123  /* Now we determine the indices of the different winding sectors. These
124  * are used for optimization purposes during tower sampling */
125  for (int n = 0; n < numWind-1; n++) {
126  /* if (abs(dot(winding(n),winding(n)) - dot(winding(n+1),winding(n+1))) > EPS) */
127  if ( abs( max(abs(winding[n+1])) - max(abs(winding[n])) ) > EPS)
128  windingSector.push_back(n);
129  }
130  /* Add the last index */
131  if (windingSector.back() != numWind-1)
132  windingSector.push_back(numWind-1);
133 }
int maxWind()
Get the maximum winding number sampled.
Definition: constants.h:101
TinyVector< unsigned int, NDIM > periodic
Determines which dimensions have periodic bc.
Definition: container.h:29
double sqrtLambdaTau
sqrt(Lambda * tau)
Definition: move.h:112
bool variableLength
Does the move have a variable length?
Definition: move.h:38
uint32 numAttempted
The number of attempted moves.
Definition: move.h:87
vector< iVec > winding
The winding vectors
Definition: move.h:100
Path & path
A reference to the paths.
Definition: move.h:80
int numWind
The total number of winding vectors.
Definition: move.h:105
ensemble operateOnConfig
What configurations do we operate on?
Definition: move.h:37
MTRand & random
A reference to the RNG.
Definition: move.h:82
int maxWind
The largest winding number.
Definition: move.h:104
ActionBase * actionPtr
A base pointer to the action.
Definition: move.h:81
int numToMove
The number of particles moved.
Definition: move.h:88
uint32 numAccepted
The number of accepted moves.
Definition: move.h:86
vector< int > windingSector
Used to index different winding sectors.
Definition: move.h:101
Array< uint32, 1 > numAttemptedLevel
The number of moves attempted at each level.
Definition: move.h:95
bool success
Did we sucessfully perform a move?
Definition: move.h:84
double sqrt2LambdaTau
sqrt(2 * Lambda * tau)
Definition: move.h:111
Array< uint32, 1 > numAcceptedLevel
The number of moves accepted at each level.
Definition: move.h:94
vector< double > cumrho0
Used for tower-sampling winding sectors.
Definition: move.h:102
const Container * boxPtr
A constant reference to the container class.
Definition: path.h:43
Ttype & max(Ttype &x, Ttype &y)
Maximum of two inputs.
Definition: common.h:146
#define NDIM
Number of spatial dimnsions.
Definition: common.h:71
int ipow(int base, int power)
Return the integer value of a number raised to a power.
Definition: common.h:136
#define EPS
A small number.
Definition: common.h:94
TinyVector< int, NDIM > iVec
A NDIM-vector of type integer.
Definition: common.h:114
ConstantParameters * constants()
Global public access to the constants.
Definition: constants.h:201
+ Here is the call graph for this function:

Member Function Documentation

◆ getWindingNumber()

iVec MoveBase::getWindingNumber ( const beadLocator startBead,
const beadLocator endBead 
)
protected

Find the winding number for a path between two beads.

Find the winding number for a given path between two beads.

By following a trajectory between two beads, accumulate the winding number by tracking when periodic boundary conditions are invoked.

Parameters
startBeadThe index of the start of the stage
endBeadThe index of the final bead in the stage
Returns
A integer NDIM-vector which holds the winding vector of the path.

Definition at line 409 of file move.cpp.

409  {
410 
411  iVec wind;
412  wind = 0;
413  beadLocator beadIndex;
414  beadIndex = startBead;
415  dVec vel;
416  do {
417  /* Get the vector separation */
418  vel = path(path.next(beadIndex)) - path(beadIndex);
419 
420  for (int i = 0; i < NDIM; i++) {
421 
422  /* Only worry about the winding number for PBC */
423  if (path.boxPtr->periodic[i]) {
424  if (vel[i] < -0.5*path.boxPtr->side[i])
425  ++wind[i];
426  if (vel[i] >= 0.5*path.boxPtr->side[i])
427  --wind[i];
428  }
429  }
430 
431  beadIndex = path.next(beadIndex);
432  } while (!all(beadIndex==endBead));
433 
434  return wind;
435 }
dVec side
The linear dimensions of the box.
Definition: container.h:31
beadLocator & next(int slice, int ptcl)
Move one link forward in imaginary time.
Definition: path.h:86
TinyVector< int, 2 > beadLocator
time-slice,bead-number world line index
Definition: common.h:117
TinyVector< double, NDIM > dVec
A NDIM-vector of type double.
Definition: common.h:111
+ Here is the call graph for this function:

◆ keepMove()

DEBUG void MoveBase::keepMove ( )
protectedvirtual

Keep the move.

For all move types, if the move is accepted, we simply increment our accept counters.

Definition at line 248 of file move.cpp.

248  {
249  numAccepted++;
250  totAccepted++;
251 
252  /* Restore the shift level for the time step to 1 */
253  actionPtr->setShift(1);
254 
255  success = true;
256 }
void setShift(int _shift)
The public method that sets the tau scaling factor.
Definition: action.h:89
static uint32 totAccepted
The total number of moves accepted.
Definition: move.h:91
+ Here is the call graph for this function:
+ Here is the caller graph for this function:

◆ newBisectionPosition()

dVec MoveBase::newBisectionPosition ( const beadLocator beadIndex,
const int  lshift 
)
protected

Returns a new bisection position which will exactly sample the kinetic action.

Parameters
beadIndexThe bead that we want a new position for
lshiftThe number of time slices between beads moved at this level
Returns
A NDIM-vector which holds a new random position.

Definition at line 469 of file move.cpp.

470  {
471 
472  /* The size of the move */
473  double delta = sqrtLambdaTau*sqrt(1.0*lshift);
474 
475  /* We first get the index and position of the 'previous' neighbor bead */
476  nBeadIndex = path.prev(beadIndex,lshift);
477 
478  /* We now get the midpoint between the previous and next beads */
479  newRanPos = path.getSeparation(path.next(beadIndex,lshift),nBeadIndex);
480  newRanPos *= 0.5;
482 
483  /* This is the gausian distributed random kick around that midpoint */
484  for (int i = 0; i < NDIM; i++)
485  newRanPos[i] = random.randNorm(newRanPos[i],delta);
486 
487  /* Put in PBC and return */
489  return newRanPos;
490 }
virtual void putInside(dVec &) const =0
Place a vector inside the simulation cell.
beadLocator nBeadIndex
Neighbor bead index.
Definition: move.h:114
dVec newRanPos
Staing random position.
Definition: move.h:116
dVec getSeparation(const beadLocator &, const beadLocator &) const
Return the separation vector between two particles in the same timeslice.
Definition: path.h:173
beadLocator & prev(int slice, int ptcl)
Move one link backward in imaginary time.
Definition: path.h:95
+ Here is the call graph for this function:
+ Here is the caller graph for this function:

◆ newFreeParticlePosition()

dVec MoveBase::newFreeParticlePosition ( const beadLocator neighborIndex)
protected

Generates a new position, which exactly samples the free particle density matrix.

Compute a position which is selected from a guassian distribution with a mean at a neighboring position, and variance equal to 2 * lambda * tau.

Parameters
neighborIndexthe beadLocator for a neighboring bead
Returns
A randomly generated position which exactly samples 1/2 the kinetic action.

Definition at line 448 of file move.cpp.

448  {
449 
450  PIMC_ASSERT(path.worm.beadOn(neighborIndex));
451 
452  /* The Gaussian distributed random position */
453  for (int i = 0; i < NDIM; i++)
454  newRanPos[i] = random.randNorm(path(neighborIndex)[i],sqrt2LambdaTau);
455 
457 
458  return newRanPos;
459 }
Worm worm
Details on the worm.
Definition: path.h:44
int beadOn(int, int) const
Safely get a bead (int indexed)
Definition: worm.h:148
#define PIMC_ASSERT(X)
Rename assert method.
Definition: common.h:64
+ Here is the call graph for this function:
+ Here is the caller graph for this function:

◆ newStagingPosition() [1/2]

dVec MoveBase::newStagingPosition ( const beadLocator neighborIndex,
const beadLocator endIndex,
const int  stageLength,
const int  k 
)
protected

Returns a new staging position which will exactly sample the kinetic action.

Parameters
neighborIndexThe index of the bead to be updated's neighbor
endIndexThe index of the final bead in the stage
stageLengthThe length of the stage
kThe position along the stage
Returns
A NDIM-vector which holds a new random position.

Definition at line 268 of file move.cpp.

269  {
270 
271  PIMC_ASSERT(path.worm.beadOn(neighborIndex));
272 
273  /* The rescaled value of lambda used for staging */
274  double f1 = 1.0 * (stageLength - k - 1);
275  double f2 = 1.0 / (1.0*(stageLength - k));
276  double sqrtLambdaKTau = sqrt2LambdaTau * sqrt(f1 * f2);
277 
278  /* We find the new 'midpoint' position which exactly samples the kinetic
279  * density matrix */
280  neighborPos = path(neighborIndex);
281  newRanPos = path(endIndex) - neighborPos;
283  newRanPos *= f2;
285 
286  /* This is the random kick around that midpoint */
287  for (int i = 0; i < NDIM; i++)
288  newRanPos[i] = random.randNorm(newRanPos[i],sqrtLambdaKTau);
289 
291 
292  return newRanPos;
293 }
void putInBC(dVec &r) const
Place a vector in boundary conditions.
Definition: container.h:50
dVec neighborPos
Staging neighbor position.
Definition: move.h:115
+ Here is the call graph for this function:
+ Here is the caller graph for this function:

◆ newStagingPosition() [2/2]

dVec MoveBase::newStagingPosition ( const beadLocator neighborIndex,
const beadLocator endIndex,
const int  stageLength,
const int  k,
iVec wind 
)
protected

Returns a new staging position which will exactly sample the kinetic action in different winding sectors.

Parameters
neighborIndexThe index of the bead to be updated's neighbor
endIndexThe index of the final bead in the stage
stageLengthThe length of the stage
kThe position along the stage
Returns
A NDIM-vector which holds a new random position.

Definition at line 305 of file move.cpp.

306  {
307 
308  PIMC_ASSERT(path.worm.beadOn(neighborIndex));
309 
310  /* The rescaled value of lambda used for staging */
311  double f1 = 1.0 * (stageLength - k - 1);
312  double f2 = 1.0 / (1.0*(stageLength - k));
313  double sqrtLambdaKTau = sqrt2LambdaTau * sqrt(f1 * f2);
314 
315  /* We find the new 'midpoint' position which exactly samples the kinetic
316  * density matrix */
317  neighborPos = path(neighborIndex);
318  newRanPos = (path(endIndex)+path.boxPtr->side*wind)-neighborPos;
319  newRanPos *= f2;
321 
322  /* This is the random kick around that midpoint */
323  for (int i = 0; i < NDIM; i++)
324  newRanPos[i] = random.randNorm(newRanPos[i],sqrtLambdaKTau);
325 
326  /* Make sure we choose the correct winding trajectory */
327  for (int i = 0; i < NDIM; i++) {
328  while (newRanPos[i] < -0.5*path.boxPtr->side[i])
329  {
330  wind[i]++;
331  newRanPos[i] += path.boxPtr->side[i];
332  }
333  while (newRanPos[i] >= 0.5*path.boxPtr->side[i])
334  {
335  wind[i]--;
336  newRanPos[i] -= path.boxPtr->side[i];
337  }
338  }
339 
341  assert(newRanPos[0] < 0.5*path.boxPtr->side[0] || newRanPos[0] >= -0.5*path.boxPtr->side[0]);
342 
343  return newRanPos;
344 }
+ Here is the call graph for this function:

◆ sampleWindingSector()

iVec MoveBase::sampleWindingSector ( const beadLocator startBead,
const beadLocator endBead,
const int  stageLength,
double &  totalrho0 
)
protected

Obtain a winding sector for a stage-like move.

Determine the winding sector to sample for a stage like move.

Perform tower sampling to select a winding sector from the free particle kinetic density matrix.

Parameters
startBeadThe index of the start of the stage
endBeadThe index of the final bead in the stage
stageLengthThe length of the stage
Returns
A integer NDIM-vector which holds the winding vector to be sampled.

Definition at line 357 of file move.cpp.

358  {
359 
360  /* For now, we hard-code the tolerance at 0.001% */
361  /* double tolerance = 1.0E-6; */
362 
363  /* Get the w = 0 sector separation */
364  dVec vel,velW;
365  vel = path(endBead) - path(startBead);
366 
367  /* If we aren't sampling winding sectors we always
368  * get the min sector. */
369  /* path.boxPtr->putInBC(vecl); */
370 
371  /* Initialize the probability and cumulative probabilities */
372  /* fill(cumrho0.begin(), cumrho0.end(), 0); */
373  cumrho0[0] = actionPtr->rho0(vel,stageLength);
374  totalrho0 = cumrho0[0];
375 
376  /* Sample the free density matrix for different winding sectors */
377  for (int n = 1; n < numWind; n++) {
378  velW = vel + winding.at(n)*path.boxPtr->side;
379  double crho0 = actionPtr->rho0(velW,stageLength);
380  totalrho0 += crho0;
381  cumrho0[n] = cumrho0[n-1] + crho0;
382  }
383 
384  /* Normalize the cumulative probability */
385  for (auto& crho0 : cumrho0)
386  crho0 /= totalrho0;
387 
388  /* for (uint32 n = 0; n < cumrho0.size(); ++n) */
389  /* cumrho0[n] /= totalrho0; */
390 
391  /* Perform tower sampling to select the winding vector */
392  int index;
393  index = std::lower_bound(cumrho0.begin(),cumrho0.end(),random.rand())
394  - cumrho0.begin();
395 
396  return winding.at(index);
397 }
double rho0(const dVec &, const dVec &, int)
The free-particle density matrix.
Definition: action.cpp:83
+ Here is the call graph for this function:
+ Here is the caller graph for this function:

The documentation for this class was generated from the following files: