Path Integral Quantum Monte Carlo
container.cpp
Go to the documentation of this file.
1 
8 #include "container.h"
9 #include "constants.h"
10 
11 // ---------------------------------------------------------------------------
12 // ---------------------------------------------------------------------------
13 // CONTAINER CLASS -----------------------------------------------------------
14 // ---------------------------------------------------------------------------
15 // ---------------------------------------------------------------------------
16 
17 /**************************************************************************/
21  side = 0.0;
22  sideInv = 0.0;
23  sideInv2 = 0.0;
24  pSide = 0.0;
25  periodic = 1;
26  maxSep = 0.0;
27  volume = 0.0;
28  rcut2 = 0.0;
29  name = "";
30  fullyPeriodic = true;
31 
32  /* Determine the number of grid boxes in the lookup table */
33  numGrid = 1;
34  for (int i = 0; i < NDIM; i++)
35  numGrid *= NGRIDSEP;
36 }
37 
38 
39 /**************************************************************************/
43 }
44 
45 /**************************************************************************/
51 double Container::gridRadius2(const int n) const {
52  iVec _gridIndex;
53  for (int i = 0; i < NDIM; i++) {
54  int scale = 1;
55  for (int j = i+1; j < NDIM; j++)
56  scale *= NGRIDSEP;
57  _gridIndex[i] = (n/scale) % NGRIDSEP;
58  PIMC_ASSERT(_gridIndex[i]<NGRIDSEP);
59  }
60 
61  double r2 = 0.0;
62  for (int i = 0; i < 2; i++) {
63  double ri = -0.5*side[i] + (_gridIndex[i] + 0.5)*gridSize[i];
64  r2 += ri*ri;
65  }
66  return r2;
67 }
68 
69 // ---------------------------------------------------------------------------
70 // ---------------------------------------------------------------------------
71 // PRISM CLASS ---------------------------------------------------------------
72 // ---------------------------------------------------------------------------
73 // ---------------------------------------------------------------------------
74 
75 /**************************************************************************/
82 Prism::Prism(double density, int numParticles) {
83 
84  /* Here we can only create a cube */
85  /* Setup the cube size in each dimension */
86  side = pow(1.0*numParticles / density, 1.0/(1.0*NDIM));
87  sideInv = 1.0/side;
88  sideInv2 = 2.0*sideInv;
89 
90  rcut2 = 0.25*side[NDIM-1]*side[NDIM-1];
91 
92  /* The hyper cube has periodic boundary conditions */
93  periodic = 1;
94  pSide = side;
95  fullyPeriodic = true;
96 
97  /* Compute the maximum possible separation possible inside the box */
98  maxSep = sqrt(dot(side/(periodic + 1.0),side/(periodic + 1.0)));
99 
100  /* Calculate the volume of the cube */
101  volume = product(side);
102 
103  /* The grid size for the lookup table */
105 
106  name = "Prism";
107 }
108 
109 /**************************************************************************/
115 Prism::Prism(const dVec &_side, const iVec &_periodic) {
116 
117  /* Setup the cube size in each dimension */
118  side = _side;
119  sideInv = 1.0/side;
120  sideInv2 = 2.0*sideInv;
121 
122  rcut2 = 0.25*side[NDIM-1]*side[NDIM-1];
123 
124  /* Setup the periodic boundary conditions */
125  periodic = _periodic;
126  pSide = periodic*side;
127 
128  /* are there any non-periodic boundary conditions? */
129  fullyPeriodic = all(periodic==1);
130 
131  /* Compute the maximum possible separation possible inside the box */
132  maxSep = sqrt(dot(side/(periodic + 1.0),side/(periodic + 1.0)));
133 
134  /* Calculate the volume of the cube */
135  volume = product(side);
136 
137  /* The grid size for the lookup table */
139 
140  name = "Prism";
141 }
142 
143 /**************************************************************************/
147 }
148 
149 /**************************************************************************/
152 dVec Prism::randPosition(MTRand &random) const {
153  dVec randPos;
154  for (int i = 0; i < NDIM; i++)
155  randPos[i] = side[i]*(-0.5 + random.randExc());
156  return randPos;
157 }
158 
159 /**************************************************************************/
162 dVec Prism::randUpdate (MTRand &random, const dVec &pos) const {
163  dVec randPos;
164  randPos = pos;
165  for (int i = 0; i < NDIM; i++)
166  randPos[i] += 4.0*constants()->displaceDelta()*(-0.5 + random.rand());
167  putInside(randPos);
168 
169  /* Make sure we don't move a particle outside the box */
170  if (!fullyPeriodic) {
171  for (int i = 0; i < NDIM; i++) {
172  if (!periodic[i]) {
173  if (randPos[i] >= 0.5*side[i])
174  randPos[i] = 0.5*side[i] - 2*EPS;
175  if (randPos[i] < -0.5*side[i])
176  randPos[i] = -0.5*side[i] + 2*EPS;
177  }
178  }
179  }
180  return randPos;
181 }
182 
183 /**************************************************************************/
192 int Prism::gridIndex(const dVec &pos) const {
193 
194  int gNumber = 0;
195  for (int i = 0; i < NDIM; i++) {
196  int scale = 1;
197  for (int j = i+1; j < NDIM; j++)
198  scale *= NGRIDSEP;
199  gNumber += scale *
200  static_cast<int>(abs( pos[i] + 0.5*side[i] - EPS ) / (gridSize[i] + EPS));
201  }
202  PIMC_ASSERT(gNumber<numGrid);
203  return gNumber;
204 }
205 
206 /**************************************************************************/
212 double Prism::gridBoxVolume(const int n) const {
213  return blitz::product(gridSize);
214 }
215 
216 // ---------------------------------------------------------------------------
217 // ---------------------------------------------------------------------------
218 // CYLINDER CLASS ------------------------------------------------------------
219 // ---------------------------------------------------------------------------
220 // ---------------------------------------------------------------------------
221 
222 /**************************************************************************/
229 Cylinder::Cylinder(const double _rho, const double radius, const int numParticles) {
230 
231  /* We can only make a cylinder in 3 dimensions */
232  if (NDIM != 3) {
233  cerr << "You can only create a cylinder in 3 dimensions, change NDIM!"
234  << endl;
235  exit(EXIT_FAILURE);
236  }
237  else {
238  /* We can compute the linear length from the density, number of particles
239  * and radius */
240  double L = (1.0*numParticles)/(_rho * M_PI * radius * radius);
241 
242  /* We check to make sure that our aspect ratio is at least 2:1 */
243  if (L < 2.0*radius)
244  cerr << "L:r is smaller than 2:1!" << endl;
245 
246  /* Setup the prism size in each of the three dimensions which the
247  * cylinder will be inscribed insde of. We make it 2R X 2R X L */
248  side[0] = side[1] = 2.0 * radius;
249  side[2] = L;
250 
251  rcut2 = 0.25*side[2]*side[2];
252 
253  sideInv = 1.0/side;
254  sideInv2 = 2.0*sideInv;
255 
256  /* setup which dimensions have periodic boundary conditions, 1.0 for yes
257  * and 0.0 for no. */
258  periodic[0] = periodic[1] = 0;
259  periodic[2] = 1;
260 
261  pSide = periodic*side;
262 
263  /* Compute the maximum possible separation possible inside the box */
264  maxSep = sqrt(dot(side/(periodic + 1.0),side/(periodic + 1.0)));
265 
266  /* Compute the cylinder volume. We use the radius here instead of the actual
267  * side. This is the 'active' volume */
268  volume = M_PI*radius*radius*L;
269 
270  name = "Cylinder";
271 
272  /* The grid size for the lookup table */
274 
275  } // end else
276 }
277 
278 /**************************************************************************/
286 Cylinder::Cylinder(const double radius, const double L) {
287 
288  /* We can only make a cylinder in 3 dimensions */
289  if (NDIM != 3) {
290  cerr << "You can only create a cylinder in 3 dimensions, change NDIM!"
291  << endl;
292  exit(EXIT_FAILURE);
293  }
294  else {
295  /* We check to make sure that our aspect ratio is at least 2:1 */
296  if (L < 2.0*radius)
297  cerr << "L:r is smaller than 2:1!" << endl;
298 
299  /* Setup the prism size in each of the three dimensions which the
300  * cylinder will be inscribed insde of. We make it 2R X 2R X L
301  * to account for */
302  side[0] = side[1] = 2.0 * radius;
303  side[2] = L;
304 
305  rcut2 = 0.25*side[2]*side[2];
306 
307  sideInv = 1.0/side;
308  sideInv2 = 2.0*sideInv;
309 
310  /* setup which dimensions have periodic boundary conditions, 1 for yes
311  * and 0 for no. */
312  periodic[0] = periodic[1] = 0;
313  periodic[2] = 1;
314 
315  /* Compute the maximum possible separation possible inside the box */
316  maxSep = sqrt(dot(side/(periodic + 1.0),side/(periodic + 1.0)));
317 
318  pSide = periodic*side;
319 
320  /* Compute the cylinder volume. We use the radius here instead of the actual
321  * side. This is the 'active' volume */
322  volume = M_PI*radius*radius*L;
323 
324  name = "Cylinder";
325 
326  /* The grid size for the lookup table */
328 
329  } // end else
330 }
331 
332 /**************************************************************************/
336 }
337 
338 /**************************************************************************/
342 dVec Cylinder::randPosition(MTRand &random) const {
343 
344  dVec randPos;
345  double r = 0.5*side[0]*sqrt(random.randExc());
346  double phi = random.randExc(2.0*M_PI);
347  randPos[0] = r * cos(phi);
348  randPos[1] = r * sin(phi);
349  randPos[2] = side[2]*(-0.5 + random.randExc());
350 
351  return randPos;
352 }
353 
354 /**************************************************************************/
357 dVec Cylinder::randUpdate (MTRand &random, const dVec &pos) const {
358  dVec randPos;
359  if (random.rand() < 0.5)
360  randPos = randUpdateJumpShell(random,pos);
361  else
362  randPos = randUpdateSmall(random,pos);
363  return randPos;
364 }
365 
366 /**************************************************************************/
369 dVec Cylinder::randUpdateJumpShell (MTRand &random, const dVec &pos) const {
370  dVec randPos;
371  randPos = pos;
372  double theta = atan2(pos[1],pos[0]);
373  double oldr = sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
374  double newr;
375 
376  if (random.rand() > 0.5)
377  newr = oldr + (1.00 + 2.50*random.rand());
378  else
379  newr = oldr - (1.00 + 2.50*random.rand());
380 
381  if (newr < 0.0)
382  newr *= -1.0;
383 
384  double ranTheta = M_PI*(-0.05 + 0.1*random.rand());
385  randPos[0] = newr*cos(theta + ranTheta);
386  randPos[1] = newr*sin(theta + ranTheta);
387  randPos[2] += constants()->displaceDelta()*(-0.5 + random.rand());
388 
389  putInside(randPos);
390  return randPos;
391 }
392 
393 // THIS IS WHAT WAS HERE BEFORE
394 /**************************************************************************/
397 dVec Cylinder::randUpdateSmall (MTRand &random, const dVec &pos) const {
398  dVec randPos;
399  randPos = pos;
400  for (int i = 0; i < NDIM; i++)
401  randPos[i] += constants()->displaceDelta()*(-0.5 + random.rand());
402  putInside(randPos);
403  return randPos;
404 }
405 /**************************************************************************/
408 void Cylinder::putInside(dVec &r) const {
409 
410  /* For the x and y axis, any outside position is placed near the boundary */
411  for (int i = 0; i < NDIM-1; i++) {
412  if (r[i] >= 0.5*side[i])
413  r[i] = 0.5*side[i] - EPS;
414  if (r[i] < -0.5*side[i])
415  r[i] = -0.5*side[i] + EPS;
416  }
417 
418  /* Now place the z-coordinate in PBC */
419  r[2] -= (r[2] >= 0.5*side[2])*side[2];
420  r[2] += (r[2] < -0.5*side[2])*side[2];
421 }
422 
423 /**************************************************************************/
432 int Cylinder::gridIndex(const dVec &pos) const {
433 
434  // /* Get the r and theta components */
435  // double r = sqrt(pos[0]*pos[0]+ pos[1]*pos[1]);
436  // double theta = atan2(pos[1],pos[0]);
437  // theta += (theta < 0.0)*2.0*M_PI;
438 
439  // /* Get the 3d vector index */
440  // iVec grid;
441  // grid[0] = static_cast<int>(abs(r-EPS)/(gridSize[0]+EPS));
442  // grid[1] = static_cast<int>(abs(theta-EPS)/(gridSize[1]+EPS));
443  // grid[2] = static_cast<int>(abs(pos[2] + 0.5*side[2] - EPS )/(gridSize[2] + EPS));
444 
445  // /* return the flattened index */
446  // int gNumber = grid[0]*NGRIDSEP*NGRIDSEP + grid[1]*NGRIDSEP + grid[2];
447 
448  // PIMC_ASSERT(gNumber<numGrid);
449  // return gNumber;
450 
451  int gNumber = 0;
452  for (int i = 0; i < NDIM; i++) {
453  int scale = 1;
454  for (int j = i+1; j < NDIM; j++)
455  scale *= NGRIDSEP;
456  gNumber += scale *
457  static_cast<int>(abs( pos[i] + 0.5*side[i] - EPS ) / (gridSize[i] + EPS));
458  }
459  PIMC_ASSERT(gNumber<numGrid);
460 
461  return gNumber;
462 }
463 
464 /**************************************************************************/
470 double Cylinder::gridBoxVolume(const int n) const {
471 
472  /* Get the first grid box index */
473  return blitz::product(gridSize);
474 }
475 
double displaceDelta() const
Get center of mass shift.
Definition: constants.h:54
dVec gridSize
The grid size in each dimension.
Definition: container.h:44
Container()
Initialize all variables.
Definition: container.cpp:20
double volume
The volume of the container in A^3.
Definition: container.h:35
int numGrid
The number of grid boxes for the position grid.
Definition: container.h:41
dVec sideInv
The inverse box dimensions.
Definition: container.h:32
dVec pSide
Periodic * side.
Definition: container.h:97
TinyVector< unsigned int, NDIM > periodic
Determines which dimensions have periodic bc.
Definition: container.h:29
dVec sideInv2
2 times the inverse box dimensions
Definition: container.h:33
double gridRadius2(const int) const
The radius of a grid box.
Definition: container.cpp:51
double rcut2
The smallest separation squared.
Definition: container.h:36
string name
The name of the container.
Definition: container.h:39
bool fullyPeriodic
Is the prism fully periodic?
Definition: container.h:42
dVec side
The linear dimensions of the box.
Definition: container.h:31
virtual ~Container()
Empty destructor.
Definition: container.cpp:42
double maxSep
The maximum possible separation for 2 beads on the same timeslice.
Definition: container.h:37
~Cylinder()
Destructor.
Definition: container.cpp:335
dVec randUpdateJumpShell(MTRand &, const dVec &) const
Return a random position close to the supplied one.
Definition: container.cpp:369
dVec randUpdateSmall(MTRand &, const dVec &) const
Return a random position close to the supplied one.
Definition: container.cpp:397
Cylinder(const double, const double, const int)
Create a cylinder given density, radius and number of particles.
Definition: container.cpp:229
dVec randPosition(MTRand &) const
Return a random position inside the cylinder.
Definition: container.cpp:342
double gridBoxVolume(const int) const
Given a grid index, return the hyper volume of the associated grid box.
Definition: container.cpp:470
dVec randUpdate(MTRand &, const dVec &) const
Return a random position close to the supplied one.
Definition: container.cpp:357
void putInside(dVec &r) const
Make sure that a suplied vector is put inside the cylinder.
Definition: container.cpp:408
int gridIndex(const dVec &) const
Given a particle position, return a single integer which maps to a unique grid position.
Definition: container.cpp:432
double gridBoxVolume(const int) const
Given a grid index, return the hyper volume of the associated grid box.
Definition: container.cpp:212
~Prism()
Empty destructor.
Definition: container.cpp:146
dVec randPosition(MTRand &) const
Return a random position inside the cube.
Definition: container.cpp:152
Prism(const double, const int)
Create a NDIM-dimensional hyperprism given density and number of particles.
Definition: container.cpp:82
int gridIndex(const dVec &) const
Given a particle position, return a single integer which maps to a unique grid position.
Definition: container.cpp:192
void putInside(dVec &r) const
For PBC, this is identical to putInBC.
Definition: container.h:113
dVec randUpdate(MTRand &, const dVec &) const
Return a random position close to the supplied one.
Definition: container.cpp:162
#define NDIM
Number of spatial dimnsions.
Definition: common.h:71
#define EPS
A small number.
Definition: common.h:94
#define NGRIDSEP
Spatial separations to be used in each dimension of the particle position grid.
Definition: common.h:93
TinyVector< double, NDIM > dVec
A NDIM-vector of type double.
Definition: common.h:111
#define PIMC_ASSERT(X)
Rename assert method.
Definition: common.h:64
TinyVector< int, NDIM > iVec
A NDIM-vector of type integer.
Definition: common.h:114
ConstantParameters class definition.
ConstantParameters * constants()
Global public access to the constants.
Definition: constants.h:201
The simulation cell.