MoveBase Class Reference

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
Path &	path
	A reference to the paths.
ActionBase *	actionPtr
	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< iVec >	winding
	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.

                                                    :
    operateOnConfig(_operateOnConfig),
    variableLength(_varLength),
    path(_path), 
    actionPtr(_actionPtr), 
    random(_random) {
 
    /* Initialize private data to zero */
    numAccepted = numAttempted = numToMove = 0;
    success = false;
 
    /* Setup the move attempt/accept arrays */
    int b  = int (ceil(log(1.0*constants()->Mbar()) / log(2.0)-EPS));
    numAcceptedLevel.resize(b+1);
    numAttemptedLevel.resize(b+1);
 
    numAcceptedLevel = 0;
    numAttemptedLevel = 0;
 
    sqrtLambdaTau = sqrt(constants()->lambda() * constants()->tau());
    sqrt2LambdaTau = sqrt(2.0)*sqrtLambdaTau;
    
    /* Setup the free density matrix arrays for sampling different
     * winding sectors.  We will sample w = -maxWind ... maxWind */
    maxWind = constants()->maxWind();
    numWind = ipow(2*maxWind + 1,NDIM);
 
    /* initialize the cumulative probability distribution */
    cumrho0.resize(numWind);
 
    /* For each integer labelling a winding sector, we construct the winding
     * vector and append to a matrix */
    iVec wind;
    for (int n = 0; n < numWind; n++ ) {
        for (int i = 0; i < NDIM; i++) {
            int scale = 1;
            for (int j = i+1; j < NDIM; j++) 
                scale *= (2*maxWind + 1);
            wind[i] = (n/scale) % (2*maxWind + 1);
 
            /* Wrap into the appropriate winding sector */
            wind[i] -= (wind[i] > maxWind)*(2*maxWind + 1);
        }
 
        /* Adjust for any non-periodic boundary conditions */
        for (int i = 0; i < NDIM; i++) {
            if (!path.boxPtr->periodic[i])
                wind[i] = 0;
        }
 
        /* Store the winding number */
        winding.push_back(wind);
    }
 
    /* Now we would like to sort the winding number array for the effecient
     * calculation of maximal probabilities. We sort based on the winding
     * sector. */
    std::stable_sort(winding.begin(), winding.end(), [](const iVec& w1, const iVec& w2) {
        return blitz::max(abs(w1)) < blitz::max(abs(w2));
        /* return (blitz::dot(w1,w1) < blitz::dot(w2,w2)); */
    });
 
    /* Now we determine the indices of the different winding sectors.  These
     * are used for optimization purposes during tower sampling */
    for (int n = 0; n < numWind-1; n++) {
        /* if (abs(dot(winding(n),winding(n)) - dot(winding(n+1),winding(n+1))) > EPS) */
        if ( abs( max(abs(winding[n+1])) - max(abs(winding[n])) ) > EPS)
            windingSector.push_back(n);
    }
    /* Add the last index */
    if (windingSector.back() != numWind-1)
        windingSector.push_back(numWind-1);
}

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

startBead	The index of the start of the stage
endBead	The 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.

                                                                                        { 
 
    iVec wind;
    wind = 0;
    beadLocator beadIndex;
    beadIndex = startBead;
    dVec vel;
    do {
        /* Get the vector separation */
        vel = path(path.next(beadIndex)) - path(beadIndex);
 
        for (int i = 0; i < NDIM; i++) {
 
            /* Only worry about the winding number for PBC */
            if (path.boxPtr->periodic[i]) {
                if (vel[i] < -0.5*path.boxPtr->side[i])
                    ++wind[i];
                if (vel[i] >= 0.5*path.boxPtr->side[i])
                    --wind[i];
            }
        }
 
        beadIndex = path.next(beadIndex);
    } while (!all(beadIndex==endBead));
 
    return wind;
}

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.

                        {
    numAccepted++;
    totAccepted++;
 
    /* Restore the shift level for the time step to 1 */
    actionPtr->setShift(1);
 
    success = true;
}

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

beadIndex	The bead that we want a new position for
lshift	The 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.

                          { 
 
    /* The size of the move */
    double delta = sqrtLambdaTau*sqrt(1.0*lshift);
 
    /* We first get the index and position of the 'previous' neighbor bead */
    nBeadIndex = path.prev(beadIndex,lshift);
 
    /* We now get the midpoint between the previous and next beads */
    newRanPos = path.getSeparation(path.next(beadIndex,lshift),nBeadIndex);
    newRanPos *= 0.5;
    newRanPos += path(nBeadIndex);
 
    /* This is the gausian distributed random kick around that midpoint */
    for (int i = 0; i < NDIM; i++)
        newRanPos[i] = random.randNorm(newRanPos[i],delta);
 
    /* Put in PBC and return */
    path.boxPtr->putInside(newRanPos);
    return newRanPos;
}

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

neighborIndex

the 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.

                                                                       {
 
    PIMC_ASSERT(path.worm.beadOn(neighborIndex));
 
    /* The Gaussian distributed random position */
    for (int i = 0; i < NDIM; i++)
        newRanPos[i] = random.randNorm(path(neighborIndex)[i],sqrt2LambdaTau);
 
    path.boxPtr->putInside(newRanPos);
 
    return newRanPos;
}

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

neighborIndex	The index of the bead to be updated's neighbor
endIndex	The index of the final bead in the stage
stageLength	The length of the stage
k	The position along the stage

Returns

A NDIM-vector which holds a new random position.

Definition at line 268 of file move.cpp.

                                            {
 
    PIMC_ASSERT(path.worm.beadOn(neighborIndex));
 
    /* The rescaled value of lambda used for staging */
    double f1 = 1.0 * (stageLength - k - 1);
    double f2 = 1.0 / (1.0*(stageLength - k));
    double sqrtLambdaKTau = sqrt2LambdaTau * sqrt(f1 * f2);
 
    /* We find the new 'midpoint' position which exactly samples the kinetic 
     * density matrix */
    neighborPos = path(neighborIndex);
    newRanPos = path(endIndex) - neighborPos;
    path.boxPtr->putInBC(newRanPos);
    newRanPos *= f2;
    newRanPos += neighborPos;
 
    /* This is the random kick around that midpoint */
    for (int i = 0; i < NDIM; i++)
        newRanPos[i] = random.randNorm(newRanPos[i],sqrtLambdaKTau);
 
    path.boxPtr->putInside(newRanPos);
 
    return newRanPos;
}

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

neighborIndex	The index of the bead to be updated's neighbor
endIndex	The index of the final bead in the stage
stageLength	The length of the stage
k	The position along the stage

Returns

A NDIM-vector which holds a new random position.

Definition at line 305 of file move.cpp.

                                                        {
    
    PIMC_ASSERT(path.worm.beadOn(neighborIndex));
    
    /* The rescaled value of lambda used for staging */
    double f1 = 1.0 * (stageLength - k - 1);
    double f2 = 1.0 / (1.0*(stageLength - k));
    double sqrtLambdaKTau = sqrt2LambdaTau * sqrt(f1 * f2);
    
    /* We find the new 'midpoint' position which exactly samples the kinetic 
     * density matrix */
    neighborPos = path(neighborIndex);
    newRanPos = (path(endIndex)+path.boxPtr->side*wind)-neighborPos;
    newRanPos *= f2;
    newRanPos += neighborPos;
    
    /* This is the random kick around that midpoint */
    for (int i = 0; i < NDIM; i++)
        newRanPos[i] = random.randNorm(newRanPos[i],sqrtLambdaKTau);
    
    /* Make sure we choose the correct winding trajectory */
    for (int i = 0; i < NDIM; i++) {
        while (newRanPos[i] < -0.5*path.boxPtr->side[i])
        {
            wind[i]++;
            newRanPos[i] += path.boxPtr->side[i];
        }
        while (newRanPos[i] >= 0.5*path.boxPtr->side[i])
        {
            wind[i]--;
            newRanPos[i] -= path.boxPtr->side[i];
        }
    }
 
    path.boxPtr->putInside(newRanPos);
    assert(newRanPos[0] < 0.5*path.boxPtr->side[0] || newRanPos[0] >= -0.5*path.boxPtr->side[0]);
    
    return newRanPos;
}

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

startBead	The index of the start of the stage
endBead	The index of the final bead in the stage
stageLength	The 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.

                                                  {
 
    /* For now, we hard-code the tolerance at 0.001% */
    /* double tolerance = 1.0E-6; */
 
    /* Get the w = 0 sector separation */ 
    dVec vel,velW;
    vel = path(endBead) - path(startBead);
 
    /* If we aren't sampling winding sectors we always
     * get the min sector. */
    /* path.boxPtr->putInBC(vecl); */
 
    /* Initialize the probability and cumulative probabilities */
    /* fill(cumrho0.begin(), cumrho0.end(), 0); */
    cumrho0[0] = actionPtr->rho0(vel,stageLength);
    totalrho0 = cumrho0[0];
 
    /* Sample the free density matrix for different winding sectors */
    for (int n = 1; n < numWind; n++) {
        velW = vel + winding.at(n)*path.boxPtr->side;
        double crho0 = actionPtr->rho0(velW,stageLength);
        totalrho0 += crho0;
        cumrho0[n] = cumrho0[n-1] + crho0;
    }
 
    /* Normalize the cumulative probability */
    for (auto& crho0 : cumrho0)
        crho0 /= totalrho0;
 
    /* for (uint32 n = 0; n < cumrho0.size(); ++n) */
    /*    cumrho0[n] /= totalrho0; */
 
    /* Perform tower sampling to select the winding vector */
    int index;
    index = std::lower_bound(cumrho0.begin(),cumrho0.end(),random.rand())
            - cumrho0.begin();
 
    return winding.at(index);
}

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:

• include/move.h
• src/move.cpp