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

radius	The radius of the cylinder
deltaRadius	R(z=0) = radius - deltaRadius
deltaWidth	The width of the hourglass constriction

Definition at line 1273 of file potential.cpp.

                                                           : PotentialBase()
{
    /* The radius of the tube */
    R = radius;
 
    /* The modification in radius: R(z=0) = R - dR */
    dR = deltaRadius;
 
    /* the length of the pore */
    L = constants()->L();
 
    /* The inverse length scale over which the radius changes */
    invd = 1.0/deltaWidth;
    R0 = 1.0/tanh(0.5*L*invd);
 
    /* The density of nitrogen in silicon nitride */
    density = 0.078; // atoms / angstrom^3
 
    /* These are the values we have historically used for amorphous silicon
     * nitride. */ 
    epsilon = 10.22;    // Kelvin
    sigma   = 2.628;    // angstroms
}

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.

                                    {
            return (0.0*pos);
        }

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

boxPtr	A pointer to the simulation cell
random	The random number generator
numParticles	The initial number of particles

Returns

An array of classical particle positions

Reimplemented from PotentialBase.

Definition at line 1429 of file potential.cpp.

                                                {
 
    /* The particle configuration */
    Array<dVec,1> initialPos(numParticles);
    initialPos = 0.0;
 
    dVec pos;
    pos = 0.0;
 
    /* We randomly place the particles inside the cylinder taking acount of the 
     * pinched radius. 
     * @see http://mathworld.wolfram.com/DiskPointPicking.html
     */
    for (int n = 0; n < numParticles; n++) {
 
            /* Uniform position along the pore */
            pos[NDIM-1] = L*(-0.5 + random.rand());
 
            /* Uniform position in a disk of z-dependent radius*/
            double theta = 2.0*M_PI*random.rand();
            double r = (Rtanh(pos[NDIM-1]) - sigma)*sqrt(random.rand());
 
            pos[0] = r*cos(theta);
            pos[1] = r*sin(theta);
 
            boxPtr->putInside(pos);
 
            initialPos(n) = pos;
    }
 
    return initialPos;
}

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

boxPtr	A pointer to the simulation cell
random	The random number generator
numParticles	The initial number of particles

Returns

An array of classical particle positions

Definition at line 1354 of file potential.cpp.

                                                {
 
    /* The particle configuration */
    Array<dVec,1> initialPos(numParticles);
    initialPos = 0.0;
 
    /* We shift the radius inward to account for the excluded volume from the
     * hard wall.  This represents the largest prism that can be put
     * inside a cylinder. */
    dVec lside;
    lside[0] = lside[1] = sqrt(2.0)*(R-dR-sigma);
    lside[2] = boxPtr->side[NDIM-1];
 
    /* Make sure our radius isn't too small */
    PIMC_ASSERT(lside[0]>0.0);
 
    /* Get the linear size per particle */
    double initSide = pow((1.0*numParticles/product(lside)),-1.0/(1.0*NDIM));
 
    /* We determine the number of initial grid boxes there are in 
     * in each dimension and compute their size */
    int totNumGridBoxes = 1;
    iVec numNNGrid;
    dVec sizeNNGrid;
 
    for (int i = 0; i < NDIM; i++) {
        numNNGrid[i] = static_cast<int>(ceil((lside[i] / initSide) - EPS));
 
        /* Make sure we have at least one grid box */
        if (numNNGrid[i] < 1)
            numNNGrid[i] = 1;
 
        /* Compute the actual size of the grid */
        sizeNNGrid[i] = lside[i] / (1.0 * numNNGrid[i]);
 
        /* Determine the total number of grid boxes */
        totNumGridBoxes *= numNNGrid[i];
    }
 
    /* Now, we place the particles at the middle of each box */
    PIMC_ASSERT(totNumGridBoxes>=numParticles);
    dVec pos;
    for (int n = 0; n < totNumGridBoxes; n++) {
 
        iVec gridIndex;
        for (int i = 0; i < NDIM; i++) {
            int scale = 1;
            for (int j = i+1; j < NDIM; j++) 
                scale *= numNNGrid[j];
            gridIndex[i] = (n/scale) % numNNGrid[i];
        }
 
        for (int i = 0; i < NDIM; i++) 
            pos[i] = (gridIndex[i]+0.5)*sizeNNGrid[i] - 0.5*lside[i];
 
        boxPtr->putInside(pos);
 
        if (n < numParticles)
            initialPos(n) = pos;
        else 
            break;
    }
 
    return initialPos;
}

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

r	the position of a particle

Returns

the confinement potential

Reimplemented from PotentialBase.

Definition at line 1314 of file potential.cpp.

                                            {
 
    double x = 0.0;
    for (int i = 0; i < NDIM-1; i++)
        x += r[i]*r[i];
 
    double Rz = Rtanh(r[NDIM-1]);
 
    x = sqrt(x)/Rz;
 
    if (x >= 1.0)
        x = 1.0 - EPS;
    
    double x2 = x*x;
    double x4 = x2*x2;
    double x6 = x2*x4;
    double x8 = x4*x4;
    double f1 = 1.0 / (1.0 - x2);
    double sigoR3 = pow(sigma/Rz,3.0);
    double sigoR9 = sigoR3*sigoR3*sigoR3;
 
    double Kx2 = boost::math::ellint_1(x);
    double Ex2 = boost::math::ellint_2(x);
 
    double v9 = (1.0*pow(f1,9.0)/(240.0)) * (
            (1091.0 + 11156*x2 + 16434*x4 + 4052*x6 + 35*x8)*Ex2 - 
            8.0*(1.0 - x2)*(1.0 + 7*x2)*(97.0 + 134*x2 + 25*x4)*Kx2);
    double v3 = 2.0*pow(f1,3.0) * ((7.0 + x2)*Ex2 - 4.0*(1.0-x2)*Kx2);
 
    return ((M_PI*epsilon*sigma*sigma*sigma*density/3.0)*(sigoR9*v9 - sigoR3*v3));
}

---

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

• include/potential.h
• src/potential.cpp