PathIntegralMonteCarlo Class Reference

The main driver class for the entire path integral monte carlo program. More...

#include <pimc.h>

Public Member Functions
	PathIntegralMonteCarlo (boost::ptr_vector< Path > &, MTRand &, boost::ptr_vector< move_vector > &, boost::ptr_vector< estimator_vector > &, const bool, uint32 binSize=100)
	Constructor. More...
	~PathIntegralMonteCarlo ()
	Destructor.
void	equilStepDiagonal ()
	Diagonal Equilibration. More...
bool	equilStepRelaxmu ()
	Relax the chemical potential to find a target number of particles. More...
bool	equilStepRelaxC0 ()
	Relax the worm constant C0 until we have found the desired diagonal fraction. More...
void	equilStep (const uint32, const bool, const bool)
	Equilibration. More...
void	step ()
	PIMC step. More...
void	finalOutput ()
	Output simulation statistics to disk. More...
void	saveState (const int finalSave=0)
	Save the state of the simulation to disk, including all path, worm, move and estimator data.
void	outputPDB ()
	Output the worldline configuration to disk using PDB , suitable for plotting using vmd. More...
void	printWormState ()
string	printHistogram ()
	Output a histogram to the terminal. More...

Data Fields
int	numStoredBins
	Number of stored estimators.
int	numDiagonal
	Number of consecutive diagonal configs.
int	numConfig
	Number of configurations;.
int	numCoMAttempted
	Number of Center of Mass moves.
int	prevNumCoMAttempted
	Number of Center of Mass moves attempted.
int	numCoMAccepted
	Number of equil CoM moves accepted.
int	numDisplaceAttempted
	Number of equil Displace moves.
int	numDisplaceAccepted
	Number of equil Displace moves accepted.
int	numMuAttempted
	Number of moves between mu adjustments.
int	numNAttempted
	The number of particle measurements.
int	numStepsAttempted
	Number of steps for relaxing C0.

Detailed Description

The main driver class for the entire path integral monte carlo program.

Holds the path, action, move and estimator objects along with methods which actually perform the monte carlo sampling procedure.

Definition at line 37 of file pimc.h.

Constructor & Destructor Documentation

PathIntegralMonteCarlo()

PathIntegralMonteCarlo::PathIntegralMonteCarlo	(	boost::ptr_vector< Path > &	_pathPtrVec,
		MTRand &	_random,
		boost::ptr_vector< move_vector > &	_movePtrVec,
		boost::ptr_vector< estimator_vector > &	_estimatorPtrVec,
		const bool	_startWithState,
		uint32	_binSize = 100
	)

Constructor.

Here we initialize all data structures, moves and estimators that will be required with peforming a path integral quantum monte carlo simulation. The initialization depends on whether or not we are restarting, or starting from a user supplied state.

Definition at line 27 of file pimc.cpp.

                                                    :
    random(_random),
    binSize(_binSize),
    Npaths(_pathPtrVec.size()),
    pathPtrVec(_pathPtrVec),
    path(pathPtrVec.front()),
    movePtrVec(_movePtrVec),
    move(movePtrVec.front()),
    estimatorPtrVec(_estimatorPtrVec),
    estimator(estimatorPtrVec.front())
{
    /* Are we starting from a saved state? */
    startWithState = _startWithState;
    
    /* Initialize stateStrings */
    stateStrings.resize(Npaths);
    
    /* We keep simulating until we have stored a certain number of measurements */
    numStoredBins = 0;
    
    /* Initialize the number of sweeps and updates per sweep*/
    numSteps = 0;
 
    /* This is required in case we have zero particles in the simulation */
    numUpdates = int(ceil(max(constants()->initialNumParticles(),1)*constants()->numTimeSlices() / 
                (1.0*constants()->Mbar())));
    if (numUpdates < 100)
        numUpdates = 100;
    
    /* Calculate the number of sweeps to make sure we touch every bead */
    numImagTimeSweeps = int(ceil((1.0*constants()->numTimeSlices()/(1.0*constants()->Mbar()))));
    
    /* These counters are used in the equilibration process */
    numConfig = 0;
    numDiagonal = 0;
    numCoMAttempted = 200;
    prevNumCoMAttempted = 200;
    numCoMAccepted = 0;
    numDisplaceAttempted = 200;
    numDisplaceAccepted = 0;
    numNAttempted = 0;
 
    numStepsAttempted = 2000*numUpdates;
    numMuAttempted = 2000*numUpdates;
 
    relaxmuMessage = false;
    relaxC0Message = false;
    equilMessage = false;
    equilODMessage = false;
 
    /* The target diagonal fraction (hard coded) */
    targetDiagFrac = 0.70;
 
    /* Do we need to shift C0 up or down? */
    sgnDiagFrac = 0;
    shiftC0 = 0.25;
 
    /* Have we found the desired C0 value for the diagonal fraction */
    foundC0 = false;
 
    /* The target/initial number of particles */
    N0 = constants()->initialNumParticles();
 
    /* Number probability distribution used when relaxing chemical potential.
     * We assume the maximum possible number is 2000 particles */
    PN.resize(2000);
    PN = 0;
    foundmu = false;
    muFactor = 1.0;
    sgnAveN = 0;
 
    /* Used for optimization of μ search */
    bestPN = 0;
    bestmu = constants()->mu();
    bestDiffAveN = 5000;
    
    /* Intialize the config number to zero */
    configNumber = 0;
    
    /* Determine the cumulative attempt move probabilities, and the indices of
     * various diagonal moves */
    double cumDiagProb = 0.0;
    double cumOffDiagProb = 0.0;
    string moveName;
 
    for (auto movePtr = move.begin(); movePtr != move.end(); ++movePtr) {
 
        /* Get the namne of the move and check if it is the generic diagonal
         * move */
        moveName = movePtr->getName();
        if ((moveName == "bisection") || (moveName == "staging"))
            moveName = "diagonal";
 
        /* Accumulate the diagonal moves */
        if ( (movePtr->operateOnConfig == DIAGONAL) || movePtr->operateOnConfig == ANY ) {
            cumDiagProb += constants()->attemptProb(moveName);
            attemptDiagProb.push_back(cumDiagProb);
        }
        else 
            attemptDiagProb.push_back(cumDiagProb);
 
        /* Accumulate the off-diagonal moves */
        if ( (movePtr->operateOnConfig == OFFDIAGONAL) || movePtr->operateOnConfig == ANY) {
            cumOffDiagProb += constants()->attemptProb(moveName);
            attemptOffDiagProb.push_back(cumOffDiagProb);
        }
        else
            attemptOffDiagProb.push_back(cumOffDiagProb);
 
        /* Find the indices of moves in case we need them */
        moveIndex[moveName] = std::distance(move.begin(), movePtr);
    }
    
    /* Make sure the move cumulative probability arrays add up to 1 */
    attemptDiagProb.back() = 1.0 + EPS;
    attemptOffDiagProb.back() = 1.0 + EPS;
    
    /* If we are restarting, or loading a state from disk, do so now. */
    if (startWithState || constants()->restart())
        loadState();
 
    /* Setup all the estimators for measurement i/o */
    for (auto &&estPtr : estimatorPtrVec) 
        for (auto &est : estPtr)
            est.prepare();
 
    /* Make a list of estimator names for the 0th estimator */
    for (auto estimatorPtr = estimator.begin(); estimatorPtr != estimator.end(); ++estimatorPtr) 
        estimatorIndex[estimatorPtr->getName()] = std::distance(estimator.begin(), estimatorPtr);
}

Here is the call graph for this function:

Member Function Documentation

equilStep()

void PathIntegralMonteCarlo::equilStep	(	const uint32	iStep,
		const bool	relaxC0,
		const bool	relaxmu
	)

Equilibration.

The equilibration method, where we perform fully diagonal moves 20% of the time, finding an optimal, COM step, then if it is desired, we attempt to find an optimal value of μ and C0.

Definition at line 561 of file pimc.cpp.

                                                                                                 {
 
    /* How far are we into the equilibration? */
    double equilFrac = (1.0*iStep) / (1.0*constants()->numEqSteps());
 
    /* We begin by attempting to perform a purely diagonal equilibration */
    if (equilFrac < 0.1 && !startWithState)  {
        if (iStep == 0) {
            numRemainingSteps = constants()->numEqSteps()/10;
            if (numRemainingSteps == 0)
                numRemainingSteps = 1;
            barStepSize = numRemainingSteps/44;
            if (!barStepSize)
                barStepSize = 1;
            cout << "[" << std::flush;
            numBars = 0;
        }
        equilStepDiagonal();
 
        /* Output a progress bar */
        if ((iStep > 0) && (iStep % barStepSize == 0) && (numBars < 44)) {
            cout << "▇" << std::flush;
            numBars++;
        }
        
        if (iStep == numRemainingSteps-1) {
            for (int nb = numBars; nb < 44; nb++)
                cout << "▇" << std::flush;
            cout << "] \t. Diagonal Pre-Equilibration." << endl << std::flush;
        }
 
    }
    else if ((equilFrac < 0.2) && (relaxmu || relaxC0)) {
 
        if (!equilODMessage) {
            equilODMessage = true;
            cout << "[";
            numRemainingSteps = int(0.2*constants()->numEqSteps())-iStep;
            barStepSize = numRemainingSteps/44;
            if (!barStepSize)
                barStepSize = 1;
            numBars = 0;
        }
 
        /* Output a progress bar */
        if ((iStep % barStepSize == 0) && (numBars < 44)) {
            cout << "▇" << std::flush;
            numBars++;
        }
 
        if (iStep == constants()->numEqSteps()/5-1) {
            for (int nb = numBars; nb < 44; nb++)
                cout << "▇" << std::flush;
            cout << "] \t. Off-Diagonal Pre-Equilibration." << endl << std::flush;
        }
 
        /* If we are performing any parameter relaxtion, we do some additional 
         * off-diagonal pre-equilibration */
        for (int n = 0; n < numUpdates; n++) {
 
            /* Generate random number and run through all moves */
            double x = random.rand();
            string mName;
            for(uint32 p=0;p<Npaths;p++){
                mName = update(x,n,p);
            }
        }
    }
    else {
        /* Perform possible parameter relaxation schemes */
        if (relaxmu && !foundmu)
            foundmu = equilStepRelaxmu();
        else if (relaxC0 && !foundC0) {
            foundC0 = equilStepRelaxC0();
        }
        /* Otherwise (or after converged) equilibrate */
        else {
            if (!equilMessage) {
                equilMessage = true;
                cout << format("[PIMCID: %s] - Equilibration Stage.") % constants()->id() << endl << "[";
                numRemainingSteps = constants()->numEqSteps()-iStep;
                barStepSize = numRemainingSteps/44;
                if (!barStepSize)
                    barStepSize = 1;
                numBars = 0;
            }
 
            /* Output a progress bar */
            if ((iStep % barStepSize == 0) && (numBars < 44)) {
                cout << "▇" << std::flush;
                numBars++;
            }
 
            /* This is the regular equilibration portion after we have the desired 
             * value of C0 and μ. */
            for (int n = 0; n < numUpdates; n++) {
 
                /* Generate random number and run through all moves */
                double x = random.rand();
                string mName;
                for(uint32 p=0;p<Npaths;p++){
                    mName = update(x,n,p);
                }
            } // numUpdates
        }
 
    } // equilfrac > 0.2
 
    /* Save a state every binsize equilibrium steps provided we are diagonal*/
    if ( path.worm.isConfigDiagonal && (iStep > 0) && (iStep % binSize) == 0) 
        saveState();
 
    if ((iStep == constants()->numEqSteps()-1) && equilMessage){
        for (int nb = numBars; nb < 44; nb++)
            cout << "▇" << std::flush;
        cout << "]" << endl;
    }
}

Here is the call graph for this function:

equilStepDiagonal()

void PathIntegralMonteCarlo::equilStepDiagonal

(

)

Diagonal Equilibration.

The diagonal-only portion of the equilibration.

Definition at line 198 of file pimc.cpp.

                                               {
 
    for (int n = 0; n < constants()->initialNumParticles(); n++) {
            
        double x = random.rand();
            
        /* Here we do the diagonal pre-equilibration moves, and allow for 
         * optimization of simulation parameters */
        if (x < constants()->attemptProb("center of mass")) {
 
            /* index for the center of mass */
            int index = moveIndex["center of mass"];
 
            /* perform the CoM move */
            for (auto & cmove : movePtrVec)
                cmove.at(index).attemptMove();
 
            /* We check how many CoM moves we have tried.  Every 200 moves, we see if we need
             * to adjust comDelta, provided we are in the pre-equilibration diagonal state. */
            if ( (move.at(index).getNumAttempted() > 0) 
                    && (move.at(index).getNumAttempted() > prevNumCoMAttempted)
                    && (move.at(index).getNumAttempted() % numCoMAttempted == 0) 
                    && (constants()->comDelta() < 0.5*blitz::min(path.boxPtr->side)) ) {
 
                numCoMAccepted  = move.at(index).getNumAccepted() - numCoMAccepted;
                double CoMRatio = 1.0*numCoMAccepted / numCoMAttempted;
                if (CoMRatio < 0.2)
                    constants()->shiftCoMDelta(-0.6);
                else if (CoMRatio < 0.3)
                    constants()->shiftCoMDelta(-0.4);
                else if (CoMRatio < 0.4)
                    constants()->shiftCoMDelta(-0.2);
                else if (CoMRatio > 0.6)
                    constants()->shiftCoMDelta(0.2);
                else if (CoMRatio > 0.7)
                    constants()->shiftCoMDelta(0.4);
                else if (CoMRatio > 0.8)
                    constants()->shiftCoMDelta(0.6);
 
                /* Reset the counters */
                numCoMAccepted = move.at(index).getNumAccepted();
                prevNumCoMAttempted = move.at(index).getNumAttempted();
            } // CoM Delta Shift
 
        } // Center of mass move
        /* Now try a displace move */
        else if (x < constants()->attemptProb("center of mass") + constants()->attemptProb("displace")) {
 
            int index = moveIndex["displace"];
            for (auto & cmove : movePtrVec)
                cmove.at(index).attemptMove();
 
            /* We check how many displace moves we have tried.  Every numDisplaceAttempted moves, we see if we need
             * to adjust delta, provided we are in the pre-equilibration diagonal state. */
            if ( (move.at(index).getNumAttempted() > 0) 
                    && (move.at(index).getNumAttempted() % numDisplaceAttempted == 0) ) {
 
                numDisplaceAccepted = move.at(index).getNumAccepted() - numDisplaceAccepted;
                double displaceRatio = 1.0*numDisplaceAccepted / numDisplaceAttempted;
                if (displaceRatio < 0.2)
                    constants()->shiftDisplaceDelta(-0.6);
                else if (displaceRatio < 0.3)
                    constants()->shiftDisplaceDelta(-0.4);
                else if (displaceRatio < 0.4)
                    constants()->shiftDisplaceDelta(-0.2);
                else if (displaceRatio > 0.6)
                    constants()->shiftDisplaceDelta(0.2);
                else if (displaceRatio > 0.7)
                    constants()->shiftDisplaceDelta(0.4);
                else if (displaceRatio > 0.8)
                    constants()->shiftDisplaceDelta(0.6);
 
                //cout << "delta: " << constants()->delta() << " " << displaceRatio << endl;
                /* Reset the counters */
                numDisplaceAccepted = move.at(index).getNumAccepted();
            }
        } // displace move
        else {
            /* Attemp a diagonal path update*/
            for (int sweep = 0; sweep < numImagTimeSweeps; sweep++){
                for (auto &cmove : movePtrVec)
                    cmove.at(moveIndex["diagonal"]).attemptMove();
            }
        }
 
    } // initNumParticles
}

Here is the call graph for this function:

Here is the caller graph for this function:

equilStepRelaxC0()

bool PathIntegralMonteCarlo::equilStepRelaxC0

(

)

Relax the worm constant C0 until we have found the desired diagonal fraction.

This portion of the equilibration attempts to find the appropriate value of the worm constant C0 until we have arrived at a desired diagonal fraction fixed to be 75% here.

Definition at line 412 of file pimc.cpp.

                                              {
 
    /* Print a message when starting the relaxation */
    if (!relaxC0Message) {
        relaxC0Message = true;
        cout << format("[PIMCID: %s] - Relax Worm Constant.\n") % constants()->id() << endl;
    }
 
    for (int n = 0; n < numUpdates; n++) {
 
        /* Generate random number and run through all moves */
        double x = random.rand();
        string mName;
        for(uint32 p=0; p<Npaths; p++) {
            mName = update(x,n,p);
        }
 
        /* We accumulate the number of diagonal configurations */
        numConfig++;
        if (path.worm.isConfigDiagonal) 
            ++numDiagonal;
 
        /* Every numStepsAttempted steps, we check the diagonal fraction and update the
         * worm constant. */
        if ( numConfig == numStepsAttempted) {
 
            double diagFrac = 1.0*numDiagonal / (1.0*numConfig);
 
            if ( (diagFrac > (targetDiagFrac-0.05)) && (diagFrac <= (targetDiagFrac+0.05)) ) {
                cout << format("\nConverged on C0 = %8.5f\n\n") % constants()->C0();
                /* for (int i = 0; i < diagFracVals.size(); i++) */ 
                /*     cout << format("%12.5f%12.5f\n") % C0Vals[i] % diagFracVals[i]; */
                return true;
            }
            else {
                cout << format("%4.2f\t%8.5f\t") % diagFrac % constants()->C0();
 
                /* Store the values of C0 and the diagonal fraciton */
                C0Vals.push_back(constants()->C0());
                diagFracVals.push_back(diagFrac);
 
                /* Perform an iterative linaer regression and suggest a new
                 * value */
                constants()->setC0(linearRegressionC0());
 
                cout << format("%8.5f\t%5d\t%8.6f\n") % constants()->C0() 
                    % path.getTrueNumParticles() 
                    % (1.0*path.getTrueNumParticles()/path.boxPtr->volume);
 
                /* Reset the counters */
                numDiagonal = 0;
                numConfig = 0;
 
            } // Haven't converged yet
 
        } // numConfig == numStepsAttempted
 
    } // numUpdates
 
    return false;
}

Here is the call graph for this function:

Here is the caller graph for this function:

equilStepRelaxmu()

bool PathIntegralMonteCarlo::equilStepRelaxmu

(

)

Relax the chemical potential to find a target number of particles.

This portion of the equilibration attempts to find the appropriate chemical potential to converge onto a target number of particles.

Definition at line 292 of file pimc.cpp.

                                              {
 
    /* Print a message when starting the relaxation */
    if (!relaxmuMessage) {
        relaxmuMessage = true;
        cout << format("[PIMCID: %s] - Relax Chemical Potential.") % constants()->id() << endl;
    }
 
    double shiftmu = 1.0;
    double mu;
 
    for (int n = 0; n < numUpdates; n++) {
 
        /* Generate random number and run through all moves */
        double x = random.rand();
        string mName;
        for(uint32 p=0; p<Npaths; p++) {
            mName = update(x,n,p);
        }
 
        /* Relax the chemical potential to get a desired average number of
         * particles */
        int cN = round(1.0*path.worm.getNumBeadsOn()/constants()->numTimeSlices());
        int curSize = PN.size();
        if (cN >= curSize) {
            PN.resizeAndPreserve(cN+1);
            for (int n = curSize; n < cN; n++)
                PN(n) = 0;
        }
 
        PN(cN)++;
        numNAttempted += 1;
 
        if ( (numNAttempted == numMuAttempted) && (N0 > 0) ) {
 
            /* Output the distribution to the terminal*/
            cout << printHistogram();
 
            /* Compute the peak loation and average number of particles */
            int peakN = blitz::maxIndex(PN)[0];
            firstIndex i;
            int aveN = round(blitz::sum(i*PN)/blitz::sum(PN));
 
            /* If we have shifted the peak to the desired value, exit */
            if (peakN == N0) {
                cout << format("Converged on μ = %8.5f\n\n") % constants()->mu();
                return true;
            }
            else {
                string method;
 
                /* We make sure the peak is in a window 10 particles away from
                 * the target. */
                bool peakInWindow = ((peakN >= N0-5) && (peakN <= N0+5));
 
                if ((PN(N0) > 0) &&  (PN(N0+1) > 0) && (PN(N0-1) > 0) && peakInWindow) {
                    method = "Equal";
                    mu = constants()->mu() - 0.5*constants()->T()*log(1.0*PN(N0+1)/PN(N0-1));
                }
                else if ((PN(N0) > 0) &&  (PN(N0+1) > 0) && peakInWindow) {
                    method = "Down";
                    mu = constants()->mu() - constants()->T()*log(1.0*PN(N0+1)/PN(N0));
                }
                else if ((PN(N0) > 0) &&  (PN(N0-1) > 0) && peakInWindow) {
                    method = "Up";
                    mu = constants()->mu() - constants()->T()*log(1.0*PN(N0)/PN(N0-1));
                }
                else {
                    if (!inWindow) {
                        cout << "Reseting μ to previous best value." << endl;
                        mu = bestmu;
                    }
                    else {
                        if (peakN > N0) {
                            mu = constants()->mu() - shiftmu*muFactor;
                            method = "Down";
                        }
                        else {
                            mu = constants()->mu() + shiftmu*muFactor;
                            method = "Up";
                        }
 
                        /* Determine if we need to alter the muFactor. */
                        int sgn = ((N0-aveN) > 0) - ((N0-aveN) < 0);
                        if (sgnAveN != 0) {
                            if (sgn == sgnAveN)
                                muFactor *= 1.618;
                            else
                                muFactor /= 1.61803399;
                        }
                        sgnAveN = sgn;
                    }
                }
 
                cout << format("Shifting %5s:%10.2f%10.2f%10d%10d%10d") % method % 
                    constants()->mu() % mu % PN(N0-1) % PN(N0) % PN(N0+1);
 
                constants()->setmu(mu);
 
                cout << format("%12d%12d%12d\n") % peakN % aveN % path.getTrueNumParticles();
 
                numNAttempted = 0;
                PN = 0;
 
            } // haven't moved peak yet
 
        } // numNAttempted == numMuAttempted
    } // numupdates
 
    return false;
}

Here is the call graph for this function:

Here is the caller graph for this function:

finalOutput()

void PathIntegralMonteCarlo::finalOutput

(

)

Output simulation statistics to disk.

We perform this once at the very end of the simulation. It saves the details of accetance probabilities as well as estimators to a file.

Definition at line 752 of file pimc.cpp.

                                         {
 
    /* Output the acceptance data to the log file */
    communicate()->file("log")->stream() << endl;
    communicate()->file("log")->stream() << endl;
 
    communicate()->file("log")->stream() << "---------- Begin Acceptance Data ---------------" << endl;
    communicate()->file("log")->stream() << endl;
    communicate()->file("log")->stream() << format("%-29s\t:\t%7.5f\n") % "Total Rate" 
        % move.front().getTotAcceptanceRatio();
    communicate()->file("log")->stream() << endl;
 
    /* Ouptut all the move acceptance information to disk */
    string moveName;
    for (auto &cmove : move){
 
        moveName = cmove.getName();
 
        /* We only output levels for those moves which have a variable size */
        if (cmove.variableLength) {
            for (int n = 0; n <= constants()->b(); n++) {
                communicate()->file("log")->stream() << format("%-12s Level %-10d\t:\t%7.5f\t(%d/%d)\n") 
                    % moveName % n % cmove.getAcceptanceRatioLevel(n) 
                    % cmove.numAcceptedLevel(n) % cmove.numAttemptedLevel(n);
            }
        }
        communicate()->file("log")->stream() << format("%-29s\t:\t%7.5f\t(%d/%d)\n") % moveName
            % cmove.getAcceptanceRatio() % cmove.numAccepted % cmove.numAttempted;
        communicate()->file("log")->stream() << endl;
    
    }
    communicate()->file("log")->stream() << "---------- End Acceptance Data -----------------" << endl;
 
    communicate()->file("log")->stream() << endl;
    communicate()->file("log")->stream() << endl;
 
    /* Output the estimator statistics to the log file */
    communicate()->file("log")->stream() << "---------- Begin Estimator Data ----------------" << endl;
    communicate()->file("log")->stream() << endl;
    for (auto &cestimator : estimator) {
        communicate()->file("log")->stream() << format("%-33s\t:\t%16d\t%16d\n") % cestimator.getName()
            % cestimator.getNumSampled() % cestimator.getTotNumAccumulated();
    
    }
    communicate()->file("log")->stream() << endl;
    communicate()->file("log")->stream() << "---------- End Estimator Data ------------------" << endl;
}

Here is the call graph for this function:

outputPDB()

void PathIntegralMonteCarlo::outputPDB

(

)

Output the worldline configuration to disk using PDB , suitable for plotting using vmd.

We must post-process the final pdb file and split it up due to connectivity changes.

See also

For the PDB specification: http://www.wwpdb.org/documentation/format32/v3.2.html

Definition at line 1161 of file pimc.cpp.

                                       {
 
    numParticles = path.getNumParticles();
    int numTimeSlices = path.numTimeSlices;
 
    configNumber++;
 
    /* We go through all beads, and find the start and end bead for each
     * worldline, adding them to an array */
    Array <beadLocator,1> startBead,endBead;
    startBead.resize(numParticles);
    endBead.resize(numParticles);
 
    /* We sort the output by the number of beads in a worldline */
    Array <int,1> wlLength(numParticles);
    wlLength = 0;
 
    /* This is the index-beadNumber mapping array */
    Array <int,2> beadNum(numTimeSlices,numParticles);
    beadNum = 0;
 
    int numWorldLines = 0;
 
    /* Get the list of beads that are active in the simulation */
    Array <bool,2> doBead(numTimeSlices,numParticles);      
    doBead = cast<bool>(path.worm.getBeads());
 
    /* We go through each particle/worldline */
    int nwl = 0;
    int beadNumber = 0;
    for (int n = 0; n < numParticles; n++) {
 
        /* The initial bead to be moved */
        startBead(nwl) = 0,n;
 
        /* We make sure we don't try to touch the same worldline twice */
        if (doBead(startBead(nwl))) {
 
            /* Check to see if the start bead is on a worm.  If it is, we start
             * at the worm tail and end at its head. */
            if (path.worm.foundBead(path,startBead(nwl))) {
                startBead(nwl) = path.worm.tail;
                endBead(nwl)   = path.next(path.worm.head);
            }
            /* Otherwise, we loop around until we find the initial bead */
            else
                endBead(nwl) = startBead(nwl);
 
            /* Mark the beads as touched and increment the number of worldlines */
            beadLocator beadIndex;
            beadIndex = startBead(nwl);
            int length = 1;
            do {
                doBead(beadIndex) = false;
                beadNum(beadIndex) = beadNumber;
                beadNumber++;
                length++;
                beadIndex = path.next(beadIndex);
            } while (!all(beadIndex==endBead(nwl)));
 
            /* We label each trajectory by the number of particles it contains.
             * a worm is always given label 0 */
            if ((length % numTimeSlices) == 0)
                wlLength(nwl) = length/numTimeSlices;
            else 
                wlLength(nwl) = 0;
 
            nwl++;
        } // doBead
 
    } // n
    numWorldLines = nwl;
 
    /* Output the PDB header */
    communicate()->file("wl")->stream() << format("REMARK [CONFIG %04d]\n") % configNumber;
 
    /* Output the unit cell information.  It is always cubic.  Everything is scaled by
     * an overall factor for better visualization. */
 
    double scale = 10.0;                
    int i;
    communicate()->file("wl")->stream() << format("%-6s") % "CRYST1";
    for (i = 0; i < NDIM; i++) 
        communicate()->file("wl")->stream() << format("%9.3f") % (scale*path.boxPtr->side[i]);
    while (i < 3) {
        communicate()->file("wl")->stream() << format("%9.3f") % 1.0;
        i++;
    }
    communicate()->file("wl")->stream() << format("%7.2f%7.2f%7.2f %-11s%4d\n") % 90.0 % 90.0 % 90.0 % "P 1" % 1;
 
    /* We output the atom block */
    beadLocator beadIndex;
    for (int n = 0; n < numWorldLines;  n++) {
        beadIndex = startBead(n);
        do {
            /* We give the zero-time-slice bead a special name */
            if (beadIndex[0] == 0) {
                communicate()->file("wl")->stream() << format("%-6s%5d %-4s %03d %9s") % "ATOM" 
                    % beadNum(beadIndex) % "H0" % wlLength(n) % " ";
            }
            else {
                communicate()->file("wl")->stream() << format("%-6s%5d %-4s %03d %9s") % "ATOM" 
                    % beadNum(beadIndex) % "HE" % wlLength(n) % " ";
            }
 
            /* Output the coordinates in 3D */
            for (i = 0; i < NDIM; i++) {
                communicate()->file("wl")->stream() << format("%8.3f") % (scale*path(beadIndex)[i]);
            }
            while (i < 3) {
                communicate()->file("wl")->stream() << format("%8.3f") % 0.0;
                i++;
            }
            communicate()->file("wl")->stream() << format("%14s\n") % "HE";
 
            beadIndex = path.next(beadIndex);
        } while (!all(beadIndex==endBead(n)));
    }
    communicate()->file("wl")->stream() <<("TER\n");
 
    /* Now output the connect block */
    for (int n = 0; n < numWorldLines;  n++) {
        beadIndex = startBead(n);
        do {
            communicate()->file("wl")->stream() << format("%-6s%5d") % "CONECT" % beadNum(beadIndex);
            beadLocator prevIndex,nextIndex;
            prevIndex = path.prev(beadIndex);
            nextIndex = path.next(beadIndex);
            
            /* We make sure that we don't connect beads linked by periodic bondary 
             * conditions */
 
            /* Check the previous bead */
            if (path.worm.beadOn(prevIndex)) {
                dVec sep;
                sep = path(beadIndex) - path(prevIndex);
                if (dot(sep,sep) < path.boxPtr->rcut2)
                    communicate()->file("wl")->stream() << format("%5d") % beadNum(prevIndex);
            }
 
            /* Now the next bead */
            if (path.worm.beadOn(nextIndex)) {
                dVec sep;
                sep = path(nextIndex) - path(beadIndex);
                if (dot(sep,sep) < path.boxPtr->rcut2)
                    communicate()->file("wl")->stream() << format("%5d") % beadNum(nextIndex);
            }
            communicate()->file("wl")->stream() << endl;
 
            beadIndex = path.next(beadIndex);
        } while (!all(beadIndex==endBead(n)));
    }
    communicate()->file("wl")->stream() <<("END\n");
 
    /* Free up memory */
    startBead.free();
    endBead.free();
    wlLength.free();
    beadNum.free();
    doBead.free();
}

Here is the call graph for this function:

printHistogram()

string PathIntegralMonteCarlo::printHistogram

(

)

Output a histogram to the terminal.

Useful for visualizing distributions.

Definition at line 1361 of file pimc.cpp.

                                              {
 
    string histogramDisplay = "\n";
 
    /* Setup the window of the histogram to be analyzed.*/
    int start,end,window;
    window = 5;
 
    start = N0-window;
    if (start < 0) 
        start = 0;
    end = N0 + window;
 
    /* We make sure we haven't moved past the end of the probability
     * distribution array */
    if (end > PN.size()-1)
        end = PN.size()-1;
 
    double factor = 50.0/blitz::max(PN);
    int sumPN = 0;
    for (int n = start; n <= end; n++) {
        int numStars = int(PN(n)*factor);
        histogramDisplay +=  str(format("%03d|") % n);
        for (int i = 0; i < numStars; i++)
            histogramDisplay +=  "▇";
        histogramDisplay += "\n";
        sumPN += PN(n);
    }
    histogramDisplay += "\n";
 
    /* if (sumPN > bestPN) { */
    /*     bestPN = sumPN; */
        /* bestmu = constants()->mu(); */
        /* cout << "setting mu" << endl; */
    /* } */
 
    firstIndex i;
    double aveN = (blitz::sum(i*PN)/blitz::sum(PN));
 
    double diffAveN = abs(N0-aveN);
 
    if (diffAveN < bestDiffAveN) {
        bestDiffAveN = diffAveN;
        bestmu = constants()->mu();
    }
 
    inWindow = (sumPN > 0);
    if (bestPN == 0) 
        inWindow = true;
 
    return histogramDisplay;
}

Here is the call graph for this function:

Here is the caller graph for this function:

step()

void PathIntegralMonteCarlo::step

(

)

PIMC step.

This method performs the metropolis sampling, which is actually a complicated multi-step operation which consists of various types of moves.
They can in general occur at different frequencies. We also measure all estimators.

Definition at line 688 of file pimc.cpp.

                                  {
 
    string moveName;
 
    /* perform updates on each set of paths */
    for (uint32 pIdx=0; pIdx<Npaths; pIdx++) {
 
        /* We run through all moves, making sure that we could have touched each bead at least once */
        for (int n = 0; n < numUpdates ; n++)  {
            moveName = update(random.rand(),n,pIdx);
        }
 
        /* Perform all measurements */
        for (auto& est : estimatorPtrVec[pIdx])
            est.sample();
        
        /* Every binSize measurements, we output averages to disk and record the
         * state of the simulation on disk.  */
        if (estimatorPtrVec[pIdx].size() > 0){
            if (estimatorPtrVec[pIdx].front().getNumAccumulated() >= binSize) {
 
                for (auto& est : estimatorPtrVec[pIdx]) {
                    /* cout << est.getNumAccumulated() << endl; */
                    if (est.getNumAccumulated() >= binSize)
                        est.output();
                }
                if(Npaths==1) 
                    saveState();
                if (pIdx == 0)
                    ++numStoredBins;
            }
        }
    }
        
    if(estimatorPtrVec.size() > Npaths) {
 
        /* Multi-Path estimators are at the end of the estimator vectors */
        for (auto& est : estimatorPtrVec.back()) 
            est.sample();
 
        /* Every binSize measurements, we output averages to disk and record the
         * state of the simulation on disk.  */
        if (estimatorPtrVec.back().front().getNumAccumulated() >= binSize) {
 
            for (auto& est : estimatorPtrVec.back()) 
                if (est.getNumAccumulated() >= binSize) 
                    est.output();
 
            /* Save to disk or store a state file */
            saveState();
 
            if (estimator.size() == 0)
                ++numStoredBins;
        }
    }
}

Here is the call graph for this function:

---

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

• include/pimc.h
• src/pimc.cpp