pimc.cpp

Go to the documentation of this file.

 
#include "pimc.h"
#include "path.h"
#include "lookuptable.h"
#include "move.h"
 
// ---------------------------------------------------------------------------
// ---------------------------------------------------------------------------
// PATH INTEGRAL MONTE CARLO CLASS -------------------------------------------
// ---------------------------------------------------------------------------
// ---------------------------------------------------------------------------
 
/**************************************************************************/
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) :
    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);
}
 
/**************************************************************************/
PathIntegralMonteCarlo::~PathIntegralMonteCarlo () {
    PN.free();
}
 
/**************************************************************************/
string PathIntegralMonteCarlo::update(const double x, const int sweep, const int pathIdx=0) { 
 
    success = false;
    string moveName = "NONE";
    int index;
 
    /* Determine the index of the move to be performed */
    if (pathPtrVec[pathIdx].worm.isConfigDiagonal)
        index = std::lower_bound(attemptDiagProb.begin(),attemptDiagProb.end(),x)
            - attemptDiagProb.begin();
    else 
        index = std::lower_bound(attemptOffDiagProb.begin(),attemptOffDiagProb.end(),x)
            - attemptOffDiagProb.begin();
 
    /* Perform the move */
    moveName = movePtrVec[pathIdx].at(index).getName();
    success = movePtrVec[pathIdx].at(index).attemptMove();
 
    return moveName;
}
 
 
/**************************************************************************/
void PathIntegralMonteCarlo::equilStepDiagonal() {
 
    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
}
 
/**************************************************************************/
bool PathIntegralMonteCarlo::equilStepRelaxmu() {
 
    /* 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;
}
 
/**************************************************************************/
bool PathIntegralMonteCarlo::equilStepRelaxC0() {
 
    /* 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;
}
 
/**************************************************************************/
double PathIntegralMonteCarlo::linearRegressionC0() {
 
    auto num = diagFracVals.size();
    double diagFrac = diagFracVals.back();
 
    if (num == 1) {
 
        sgnDiagFrac = ((targetDiagFrac-diagFrac) > 0) - ((targetDiagFrac-diagFrac) < 0);
 
        if (diagFrac < targetDiagFrac) {
            /* double C0new = C0Vals.back() - shiftC0; */
            /* return (C0new > 0.0) ? C0new : 1.0E-4; */
            return 0.4*C0Vals.back();
        }
        else
            return 1.6*C0Vals.back();
            /* return C0Vals.back() + shiftC0; */
    }
    else {
 
        /* We first find the closest values to targetC0 */
        /* double closestC0; */
        /* double minDiff = 1000; */
        /* for (int i = 0; i < num; i++) { */
        /*     double cDiff = abs(diagFracVals[i]-targetC0); */
        /*     if (cDiff < minDiff) { */
        /*         closestC0 = C0Vals[i]; */
        /*         minDiff = cDiff; */
        /*     } */
        /* } */
 
        double sumX,sumY,sumX2,sumXY,sum1;
 
        sum1 = sumX = sumY = sumXY = sumX2 = 0.0;
 
        for (int i = 0; i < num; i++) {
            sumX += C0Vals[i];
            sumX2 += C0Vals[i]*C0Vals[i];
            sumY += diagFracVals[i];
            sumXY += C0Vals[i]*diagFracVals[i];
            sum1 += 1.0;
        }
 
        double denom = sum1*sumX2 - sumX*sumX;
            double a0 = (sumY*sumX2 - sumX*sumXY)/denom;
        double a1 = (sum1*sumXY - sumY*sumX)/denom;
 
        double C0guess = (targetDiagFrac-a0)/a1;
 
        if ((C0guess < 1E4) && ( C0guess > 1E-4) && (a1 < 0.0))
            return C0guess;
        else
        {
 
            /* Determine if we need to alter the muFactor. */
            
            /* This is commented out for now, return to see if this can be
             * improved. */
            /* int sgn = ((targetDiagFrac-diagFrac) > 0) - ((targetDiagFrac-diagFrac) < 0); */
            /* if (sgn == sgnDiagFrac) */
            /*     shiftC0 *= 1.618; */
            /* else */
            /*     shiftC0 /= 1.61803399; */
            /* sgnDiagFrac = sgn; */
 
            if (diagFrac < targetDiagFrac) {
                /* double C0new = C0Vals.back() - shiftC0; */
                /* return (C0new > 0.0) ? C0new : 1.0E-4; */
                return (1.0/1.618)*C0Vals.back();
            }
            else
                return 1.61803399*C0Vals.back();
                /* return C0Vals.back() + shiftC0; */
        }
    }
}
 
/**************************************************************************/
void PathIntegralMonteCarlo::equilStep(const uint32 iStep, const bool relaxC0, const bool relaxmu) {
 
    /* 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;
    }
}
 
/**************************************************************************/
void PathIntegralMonteCarlo::step() {
 
    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;
        }
    }
}
 
/**************************************************************************/
void PathIntegralMonteCarlo::finalOutput() {
 
    /* 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;
}
 
/**************************************************************************/
void PathIntegralMonteCarlo::saveState(const int finalSave) {
 
    stringstream stateStrStrm;
    
    /* We only update the string during the simulation */
    if (!finalSave) {
 
        for(uint32 pIdx=0; pIdx<Npaths; pIdx++){
 
            stateStrStrm.str("");
 
            /* leftpack and update the lookup table arrays */
            pathPtrVec[pIdx].leftPack();
            pathPtrVec[pIdx].lookup.updateGrid(path);
 
            /* We First write the current total number of world lines */
            stateStrStrm << pathPtrVec[pIdx].getNumParticles() << endl;
 
            /* Now write the total acceptance information for all moves */
            stateStrStrm << format("%16d\t%16d\n")
                % movePtrVec[pIdx].front().totAccepted % movePtrVec[pIdx].front().totAttempted;
 
            /* Now record the individual move acceptance information,
             * first for the diagonal, then off-diagonal*/
            for (const auto &cmove : movePtrVec[pIdx])
                stateStrStrm << format("%16d\t%16d\n") % cmove.numAccepted % cmove.numAttempted;
 
            /* Output the estimator sampling information */
            for (const auto &cestimator : estimatorPtrVec[pIdx])
                stateStrStrm << format("%16d\t%16d\n") % cestimator.getTotNumAccumulated()
                    % cestimator.getNumSampled();
 
            /* Now we output the actual path and worldline data */
            stateStrStrm << setprecision(16) << pathPtrVec[pIdx].beads << endl;
            stateStrStrm << pathPtrVec[pIdx].nextLink << endl;
            stateStrStrm << pathPtrVec[pIdx].prevLink << endl;
 
            /* Output the worm data */
            stateStrStrm << pathPtrVec[pIdx].worm.beads << endl;
 
            /* Save the state of the random number generator */
            uint32 randomState[random.SAVE];
            random.save(randomState);
            for (int i = 0; i < random.SAVE; i++)
                stateStrStrm << randomState[i] << " ";
            stateStrStrm << endl;
 
            /* store the state string */
            stateStrings[pIdx] = stateStrStrm.str();
 
        } // for pidx
    }
 
    /* Save the state file to disk */
    string stateFileName;
    if (constants()->saveStateFiles() || finalSave) {
 
        for(uint32 pIdx=0; pIdx<Npaths; pIdx++) {
 
            stateFileName = "state";
            if(pIdx > 0)
                stateFileName += str(format("%d") % (pIdx+1));
 
            /* Prepare the state file for writing */
            communicate()->file(stateFileName.c_str())->reset();
 
            /* Write the stateString to disk */
            communicate()->file(stateFileName.c_str())->stream() << stateStrings[pIdx];
 
            /* Rename and copy the file. */
            communicate()->file(stateFileName.c_str())->rename();
        }
    }
}
 
/**************************************************************************/
void PathIntegralMonteCarlo::loadClassicalState(Array <dVec,2> &tempBeads,
        Array <unsigned int, 2> &tempWormBeads, int numWorldLines) {
 
    /* We go through each active worldline and create a new classical
     * configuration */
    int ptcl;
    beadLocator beadIndex;
    
    for(uint32 pIdx=0; pIdx<Npaths; pIdx++){
        ptcl = 0;
        for (int n = 0; n < numWorldLines; n++) {
            beadIndex = 0,n;
 
            /* Check if the bead is on */
            if (tempWormBeads(beadIndex)) {
                
                /* Assign the classical configuration */
                pathPtrVec[pIdx].beads(Range::all(),ptcl) = tempBeads(beadIndex);
                pathPtrVec[pIdx].worm.beads(Range::all(),ptcl) = 1;
                ptcl++;
            }
        }
    }
}
 
/**************************************************************************/
void PathIntegralMonteCarlo::loadQuantumState(Array <dVec,2> &tempBeads, 
        Array <beadLocator,2> &tempNextLink, Array <beadLocator,2> &tempPrevLink,
        int numTimeSlices, int tempNumWorldLines) {
 
    /* Prevent double counting of worldlines */
    Array <bool, 1> doBead(tempNumWorldLines);
 
    beadLocator startBead,beadIndex;
    beadLocator newStartBead;
    int ptcl;
    int slice;
    
    for(uint32 pIdx=0; pIdx<Npaths; pIdx++) {
 
        newStartBead = 0,0;
        ptcl = 0;
        slice = 0;
        doBead = true;
        
        /* Now we iterate through each worldline exactly once */
        for (int n = 0; n < tempNumWorldLines; n++) {
 
            /* The initial bead to be moved */
            startBead = 0,n;
 
            /* We make sure we don't try to touch the same worldline twice */
            if (doBead(n)) {
 
                beadIndex = startBead;
 
                /* The world line length, we simply advance until we have looped back on 
                 * ourselves. */
                slice = 0;
                newStartBead = slice,ptcl;
                do {
                    /* We turn off any zero-slice beads we have touched */
                    if (beadIndex[0]==0)
                        doBead(beadIndex[1]) = false;
 
                    pathPtrVec[pIdx].beads(slice % numTimeSlices,ptcl) = tempBeads(beadIndex);
                    pathPtrVec[pIdx].worm.beads(slice % numTimeSlices,ptcl) = 1;
 
                    beadIndex = tempNextLink(beadIndex);
                    ++slice;
 
                    /* Do a forward reconnection, provided we are not at the
                     * last bead */
                    if ( ((slice % numTimeSlices) == 0) && !all(beadIndex==startBead)) {
                        pathPtrVec[pIdx].nextLink(numTimeSlices-1,ptcl) = 0,ptcl+1;
                        pathPtrVec[pIdx].prevLink(0,ptcl+1) = numTimeSlices-1,ptcl;
                        ++ptcl;
                    }
                } while (!all(beadIndex==startBead));
 
                /* Now we have to add the remaining beads and perform the final
                 * reconnection */
                for (int tslice = (slice % numTimeSlices); tslice < numTimeSlices; tslice++) {
                    pathPtrVec[pIdx].beads(tslice,ptcl) = tempBeads(beadIndex);
                    pathPtrVec[pIdx].worm.beads(tslice,ptcl) = 1;
                }
                pathPtrVec[pIdx].nextLink(numTimeSlices-1,ptcl) = newStartBead;
                pathPtrVec[pIdx].prevLink(newStartBead) = numTimeSlices-1,ptcl;
                ++ptcl;
 
            } // doBead
        } // n
    } // pIdx
 
    /* Free local memory */
    doBead.free();
 
}
 
/**************************************************************************/
void PathIntegralMonteCarlo::loadState() {
 
    string tempString;
    
    string fileInitStr = "init";
    
    for( uint32 pIdx=0; pIdx<Npaths; pIdx++){
        
        if(pIdx>0){
            fileInitStr = str(format("init%d") % (pIdx+1));
        }
 
        /* Reset the total acceptance information */
        movePtrVec[pIdx].front().resetTotAccept();
 
        /* Reset all the individual move acceptance information */
        for (auto &cmove : movePtrVec[pIdx])
            cmove.resetAccept();
 
        /* Reset estimator sampling information */
        for (auto &cestimator : estimatorPtrVec[pIdx])
            cestimator.restart(0,0);
 
        /* We first read the former total number of world lines */
        int numWorldLines;
        int numTimeSlices = pathPtrVec[pIdx].numTimeSlices;
        communicate()->file(fileInitStr)->stream() >> numWorldLines;
 
        /* Now we skip through the input file until we find the beads matrix.  This
         * is signalled by the appearance of an open bracket "(" */
        while (!communicate()->file(fileInitStr)->stream().eof()) {
            if (communicate()->file(fileInitStr)->stream().peek() != '(') 
                getline (communicate()->file(fileInitStr)->stream(), tempString);
            else
                break;
        }
 
        /* Now we resize all path data members and read them from the init state file */
        pathPtrVec[pIdx].beads.resize(numTimeSlices,numWorldLines);
        pathPtrVec[pIdx].nextLink.resize(numTimeSlices,numWorldLines);
        pathPtrVec[pIdx].prevLink.resize(numTimeSlices,numWorldLines);
        pathPtrVec[pIdx].worm.beads.resize(numTimeSlices,numWorldLines);
 
        /* A temporary container for the beads array */
        Array <dVec,2> tempBeads;
 
        /* Get the worldline configuration */
        communicate()->file(fileInitStr)->stream() >> tempBeads;
 
        /* The temporary number of time slices */
        int tempNumTimeSlices = tempBeads.rows();
 
        if (tempNumTimeSlices == numTimeSlices) {
        
            /* Copy over the beads array */
            pathPtrVec[pIdx].beads = tempBeads;
 
            /* Get the link arrays */
            communicate()->file(fileInitStr)->stream() >> pathPtrVec[pIdx].nextLink;
            communicate()->file(fileInitStr)->stream() >> pathPtrVec[pIdx].prevLink;
 
            /* Repeat for the worm file */
            communicate()->file(fileInitStr)->stream() >> pathPtrVec[pIdx].worm.beads;
 
        } // locBeads.rows() == numTimeSlices
        else {
 
            /* Initialize the links */
            firstIndex i1;
            secondIndex i2;
            pathPtrVec[pIdx].prevLink[0] = i1-1;
            pathPtrVec[pIdx].prevLink[1] = i2;
            pathPtrVec[pIdx].nextLink[0] = i1+1;
            pathPtrVec[pIdx].nextLink[1] = i2;
        
            /* Here we implement the initial periodic boundary conditions in 
             * imaginary time */
            pathPtrVec[pIdx].prevLink(0,Range::all())[0] = numTimeSlices-1;
            pathPtrVec[pIdx].nextLink(numTimeSlices-1,Range::all())[0] = 0;
 
            /* Reset the worm.beads array */
            pathPtrVec[pIdx].worm.beads = 0;
 
            /* Temporary containers for the links and worm beads */
            Array <beadLocator,2> tempNextLink;
            Array <beadLocator,2> tempPrevLink;
            Array <unsigned int,2> tempWormBeads;
 
            /* Get the link arrays and worm file */
            communicate()->file(fileInitStr)->stream() >> tempNextLink;
            communicate()->file(fileInitStr)->stream() >> tempPrevLink;
            communicate()->file(fileInitStr)->stream() >> tempWormBeads;
 
            /* Load a classical (all time slice positions equal) from the input
             * file */
            loadClassicalState(tempBeads,tempWormBeads, numWorldLines);
 
            /* Load a quantum initial state from a file */
            //if (tempNumTimeSlices < numTimeSlices) {
            //    loadQuantumState(tempBeads,tempNextLink,tempPrevLink,
            //            numTimeSlices,int(sum(tempWormBeads)/tempNumTimeSlices));
            //}
 
            /* Now we make sure all empty beads are unlinked */
            beadLocator beadIndex;
            for (int slice = 0; slice < numTimeSlices; slice++) {
                for (int ptcl = 0; ptcl < numWorldLines; ptcl++) {
                    beadIndex = slice,ptcl;
                    if (!pathPtrVec[pIdx].worm.beads(beadIndex)) {
                        pathPtrVec[pIdx].nextLink(beadIndex) = XXX;
                        pathPtrVec[pIdx].prevLink(beadIndex) = XXX;
                    }
                }
            }
 
            /* Free local memory */
            tempPrevLink.free();
            tempNextLink.free();
            tempWormBeads.free();
 
        } // locBeads.rows() != numTimeSlices
 
        /* Load the state of the random number generator, only if we are restarting 
         * the simulation */
        if (constants()->restart()) {
            uint32 randomState[random.SAVE];
            for (int i = 0; i < random.SAVE; i++) 
                communicate()->file(fileInitStr)->stream() >> randomState[i];
            random.load(randomState);
        }
 
        /* Reset the number of on beads */
        pathPtrVec[pIdx].worm.resetNumBeadsOn();
 
        /* Go through all beads, and make sure they fit inside the simulation cell.
         * At the same time, determine how many active beads there are per slice */
        beadLocator beadIndex;
        pathPtrVec[pIdx].numBeadsAtSlice = 0;
        for (beadIndex[0] = 0; beadIndex[0] < numTimeSlices; ++beadIndex[0]) {
            for (beadIndex[1] = 0; beadIndex[1] < numWorldLines; ++beadIndex[1]) {
                pathPtrVec[pIdx].boxPtr->putInside(path(beadIndex));
                if (pathPtrVec[pIdx].worm.beadOn(beadIndex))
                    ++pathPtrVec[pIdx].numBeadsAtSlice(beadIndex[0]);
            } // particles
        } // slices
 
        /* Reset the worm parameters */
        pathPtrVec[pIdx].worm.isConfigDiagonal = true;
        pathPtrVec[pIdx].worm.reset();
 
        /* Resize and update the lookup table arrays */
        pathPtrVec[pIdx].lookup.resizeList(numWorldLines);
        pathPtrVec[pIdx].lookup.updateGrid(path);
 
        /* Reset the broken/closed worldline vectors */
        pathPtrVec[pIdx].resetBrokenClosedVecs();
        
        /* Close the file */
        communicate()->file(fileInitStr)->close();
        
        /* Free up memory */
        tempBeads.free();
        
    }
}
 
/**************************************************************************/
void PathIntegralMonteCarlo::outputPDB() {
 
    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();
}
 
/**************************************************************************/
void PathIntegralMonteCarlo::printWormState() {
 
    /* We make a list of all the beads contained in the worm */
    Array <beadLocator,1> wormBeads;    // Used for debugging
    wormBeads.resize(path.worm.length+1);
    wormBeads = XXX;
 
    /* Output the worldline configuration */
    communicate()->file("debug")->stream() << " (" << path.getTrueNumParticles() << ")" << endl;
    communicate()->file("debug")->stream() << "head " << path.worm.head[0] << " " << path.worm.head[1]
        << " tail " << path.worm.tail[0] << " " << path.worm.tail[1]
        << " length " << path.worm.length 
        << " gap " << path.worm.gap << endl;
 
    if (!path.worm.isConfigDiagonal) {
        beadLocator beadIndex;
        beadIndex = path.worm.tail;
        int n = 0;
        do {
            wormBeads(n) = beadIndex;
            beadIndex = path.next(beadIndex);
            ++n;
        } while(!all(beadIndex==path.next(path.worm.head)));
    }
 
    path.printWormConfig(wormBeads);
    path.printLinks<fstream>(communicate()->file("debug")->stream());
    wormBeads.free();
}
 
/**************************************************************************/
string PathIntegralMonteCarlo::printHistogram() {
 
    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;
}