Path Integral Quantum Monte Carlo
Public Member Functions
LJHourGlassPotential Class Reference

Computes the value of the external wall potential for an hour-glass shaped cavity. More...

#include <potential.h>

+ Inheritance diagram for LJHourGlassPotential:
+ Collaboration diagram for LJHourGlassPotential:

Public Member Functions

 LJHourGlassPotential (const double, const double, const double)
 Constructor. More...
 
 ~LJHourGlassPotential ()
 Destructor.
 
double V (const dVec &)
 The integrated LJ hour glass potential. More...
 
dVec gradV (const dVec &pos)
 The gradient of the potential. More...
 
double grad2V (const dVec &pos)
 Grad^2 of the potential.
 
Array< dVec, 1 > initialConfig (const Container *, MTRand &, const int)
 Initial configuration corresponding to the LJ cylinder potential. More...
 
Array< dVec, 1 > initialConfig1 (const Container *, MTRand &, const int)
 Return a set of initial positions inside the hourglass. More...
 
- Public Member Functions inherited from PotentialBase
 PotentialBase ()
 Constructor.
 
virtual ~PotentialBase ()
 Destructor.
 
virtual double V (const dVec &, const dVec &)
 The effective potential for the pair product approximation.
 
virtual double dVdlambda (const dVec &, const dVec &)
 The derivative of the effective potential with respect to lambda and tau.
 
virtual double dVdtau (const dVec &, const dVec &)
 
void output (const double)
 A debug method that output's the potential to a supplied separation. More...
 
virtual Array< double, 1 > getExcLen ()
 Array to hold data elements. More...
 

Additional Inherited Members

- Data Fields inherited from PotentialBase
double tailV
 Tail correction factor.
 
- Protected Member Functions inherited from PotentialBase
double deltaSeparation (double sep1, double sep2) const
 Return the minimum image difference for 1D separations.
 

Detailed Description

Computes the value of the external wall potential for an hour-glass shaped cavity.

Definition at line 622 of file potential.h.

Constructor & Destructor Documentation

◆ LJHourGlassPotential()

LJHourGlassPotential::LJHourGlassPotential ( const double  radius,
const double  deltaRadius,
const double  deltaWidth 
)

Constructor.

We create a Lennard-Jones hourglass, which for now, directly evaluates the potential based on breaking up the pore into finite size pieces, each with their own radius.

See also
C. Chakravarty J. Phys. Chem. B, 101, 1878 (1997).
Parameters
radiusThe radius of the cylinder
deltaRadiusR(z=0) = radius - deltaRadius
deltaWidthThe width of the hourglass constriction

Definition at line 1273 of file potential.cpp.

1274  : PotentialBase()
1275 {
1276  /* The radius of the tube */
1277  R = radius;
1278 
1279  /* The modification in radius: R(z=0) = R - dR */
1280  dR = deltaRadius;
1281 
1282  /* the length of the pore */
1283  L = constants()->L();
1284 
1285  /* The inverse length scale over which the radius changes */
1286  invd = 1.0/deltaWidth;
1287  R0 = 1.0/tanh(0.5*L*invd);
1288 
1289  /* The density of nitrogen in silicon nitride */
1290  density = 0.078; // atoms / angstrom^3
1291 
1292  /* These are the values we have historically used for amorphous silicon
1293  * nitride. */
1294  epsilon = 10.22; // Kelvin
1295  sigma = 2.628; // angstroms
1296 }
double L() const
Get maximum side length.
Definition: constants.h:52
PotentialBase()
Constructor.
Definition: potential.cpp:25
ConstantParameters * constants()
Global public access to the constants.
Definition: constants.h:201
+ Here is the call graph for this function:

Member Function Documentation

◆ gradV()

dVec LJHourGlassPotential::gradV ( const dVec pos)
inlinevirtual

The gradient of the potential.

Use the prrimitive approx.

Reimplemented from PotentialBase.

Definition at line 632 of file potential.h.

632  {
633  return (0.0*pos);
634  }

◆ initialConfig()

Array< dVec, 1 > LJHourGlassPotential::initialConfig ( const Container boxPtr,
MTRand &  random,
const int  numParticles 
)
virtual

Initial configuration corresponding to the LJ cylinder potential.

Return a set of initial positions inside the hourglass.

Parameters
boxPtrA pointer to the simulation cell
randomThe random number generator
numParticlesThe initial number of particles
Returns
An array of classical particle positions

Reimplemented from PotentialBase.

Definition at line 1429 of file potential.cpp.

1430  {
1431 
1432  /* The particle configuration */
1433  Array<dVec,1> initialPos(numParticles);
1434  initialPos = 0.0;
1435 
1436  dVec pos;
1437  pos = 0.0;
1438 
1439  /* We randomly place the particles inside the cylinder taking acount of the
1440  * pinched radius.
1441  * @see http://mathworld.wolfram.com/DiskPointPicking.html
1442  */
1443  for (int n = 0; n < numParticles; n++) {
1444 
1445  /* Uniform position along the pore */
1446  pos[NDIM-1] = L*(-0.5 + random.rand());
1447 
1448  /* Uniform position in a disk of z-dependent radius*/
1449  double theta = 2.0*M_PI*random.rand();
1450  double r = (Rtanh(pos[NDIM-1]) - sigma)*sqrt(random.rand());
1451 
1452  pos[0] = r*cos(theta);
1453  pos[1] = r*sin(theta);
1454 
1455  boxPtr->putInside(pos);
1456 
1457  initialPos(n) = pos;
1458  }
1459 
1460  return initialPos;
1461 }
virtual void putInside(dVec &) const =0
Place a vector inside the simulation cell.
#define NDIM
Number of spatial dimnsions.
Definition: common.h:71
TinyVector< double, NDIM > dVec
A NDIM-vector of type double.
Definition: common.h:111
+ Here is the call graph for this function:

◆ initialConfig1()

Array< dVec, 1 > LJHourGlassPotential::initialConfig1 ( const Container boxPtr,
MTRand &  random,
const int  numParticles 
)

Return a set of initial positions inside the hourglass.

Parameters
boxPtrA pointer to the simulation cell
randomThe random number generator
numParticlesThe initial number of particles
Returns
An array of classical particle positions

Definition at line 1354 of file potential.cpp.

1355  {
1356 
1357  /* The particle configuration */
1358  Array<dVec,1> initialPos(numParticles);
1359  initialPos = 0.0;
1360 
1361  /* We shift the radius inward to account for the excluded volume from the
1362  * hard wall. This represents the largest prism that can be put
1363  * inside a cylinder. */
1364  dVec lside;
1365  lside[0] = lside[1] = sqrt(2.0)*(R-dR-sigma);
1366  lside[2] = boxPtr->side[NDIM-1];
1367 
1368  /* Make sure our radius isn't too small */
1369  PIMC_ASSERT(lside[0]>0.0);
1370 
1371  /* Get the linear size per particle */
1372  double initSide = pow((1.0*numParticles/product(lside)),-1.0/(1.0*NDIM));
1373 
1374  /* We determine the number of initial grid boxes there are in
1375  * in each dimension and compute their size */
1376  int totNumGridBoxes = 1;
1377  iVec numNNGrid;
1378  dVec sizeNNGrid;
1379 
1380  for (int i = 0; i < NDIM; i++) {
1381  numNNGrid[i] = static_cast<int>(ceil((lside[i] / initSide) - EPS));
1382 
1383  /* Make sure we have at least one grid box */
1384  if (numNNGrid[i] < 1)
1385  numNNGrid[i] = 1;
1386 
1387  /* Compute the actual size of the grid */
1388  sizeNNGrid[i] = lside[i] / (1.0 * numNNGrid[i]);
1389 
1390  /* Determine the total number of grid boxes */
1391  totNumGridBoxes *= numNNGrid[i];
1392  }
1393 
1394  /* Now, we place the particles at the middle of each box */
1395  PIMC_ASSERT(totNumGridBoxes>=numParticles);
1396  dVec pos;
1397  for (int n = 0; n < totNumGridBoxes; n++) {
1398 
1399  iVec gridIndex;
1400  for (int i = 0; i < NDIM; i++) {
1401  int scale = 1;
1402  for (int j = i+1; j < NDIM; j++)
1403  scale *= numNNGrid[j];
1404  gridIndex[i] = (n/scale) % numNNGrid[i];
1405  }
1406 
1407  for (int i = 0; i < NDIM; i++)
1408  pos[i] = (gridIndex[i]+0.5)*sizeNNGrid[i] - 0.5*lside[i];
1409 
1410  boxPtr->putInside(pos);
1411 
1412  if (n < numParticles)
1413  initialPos(n) = pos;
1414  else
1415  break;
1416  }
1417 
1418  return initialPos;
1419 }
dVec side
The linear dimensions of the box.
Definition: container.h:31
#define EPS
A small number.
Definition: common.h:94
#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
+ Here is the call graph for this function:

◆ V()

double LJHourGlassPotential::V ( const dVec r)
virtual

The integrated LJ hour glass potential.

Return the actual value of the LJ Cylinder potential, for a distance r from the surface of the wall.

Checked and working with Mathematica (hourglass_potential_test.nb) on 2014-09-04.

Parameters
rthe position of a particle
Returns
the confinement potential

Reimplemented from PotentialBase.

Definition at line 1314 of file potential.cpp.

1314  {
1315 
1316  double x = 0.0;
1317  for (int i = 0; i < NDIM-1; i++)
1318  x += r[i]*r[i];
1319 
1320  double Rz = Rtanh(r[NDIM-1]);
1321 
1322  x = sqrt(x)/Rz;
1323 
1324  if (x >= 1.0)
1325  x = 1.0 - EPS;
1326 
1327  double x2 = x*x;
1328  double x4 = x2*x2;
1329  double x6 = x2*x4;
1330  double x8 = x4*x4;
1331  double f1 = 1.0 / (1.0 - x2);
1332  double sigoR3 = pow(sigma/Rz,3.0);
1333  double sigoR9 = sigoR3*sigoR3*sigoR3;
1334 
1335  double Kx2 = boost::math::ellint_1(x);
1336  double Ex2 = boost::math::ellint_2(x);
1337 
1338  double v9 = (1.0*pow(f1,9.0)/(240.0)) * (
1339  (1091.0 + 11156*x2 + 16434*x4 + 4052*x6 + 35*x8)*Ex2 -
1340  8.0*(1.0 - x2)*(1.0 + 7*x2)*(97.0 + 134*x2 + 25*x4)*Kx2);
1341  double v3 = 2.0*pow(f1,3.0) * ((7.0 + x2)*Ex2 - 4.0*(1.0-x2)*Kx2);
1342 
1343  return ((M_PI*epsilon*sigma*sigma*sigma*density/3.0)*(sigoR9*v9 - sigoR3*v3));
1344 }

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