Skip to content

Commit

Permalink
checkpoint
Browse files Browse the repository at this point in the history 
Signed-off-by: Steven Hahn <[email protected]>
  • Loading branch information
@quantumsteve
quantumsteve committed Sep 29, 2023
1 parent d65681f commit 538bfcd
Show file tree
Hide file tree
Showing 10 changed files with 12 additions and 354 deletions.
1 change: 0 additions & 1 deletion src/Particle/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,6 @@ set(PARTICLE
MCWalkerConfigurationT.cpp
WalkerConfigurations.cpp
SpeciesSet.cpp
SampleStack.cpp
SampleStackT.cpp
createDistanceTableAA.cpp
createDistanceTableAB.cpp
Expand Down
246 changes: 2 additions & 244 deletions src/Particle/Reptile.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,254 +25,12 @@
#ifndef QMCPLUSPLUS_REPTILE_H
#define QMCPLUSPLUS_REPTILE_H

#include "QMCDrivers/DriftOperators.h"
#include "QMCDrivers/WalkerProperties.h"
#include "Configuration.h"
#include "Walker.h"
#include "Particle/MCWalkerConfiguration.h"
#include "Particle/ReptileT.h"

namespace qmcplusplus
{

class Reptile : public QMCTraits
{
public:
using WP = WalkerProperties::Indexes;
using Walker_t = MCWalkerConfiguration::Walker_t;
using WalkerIter_t = MCWalkerConfiguration::iterator;
using ReptileConfig_t = std::vector<Walker_t::ParticlePos>;

std::vector<IndexType> Action;
std::vector<IndexType> TransProb;

RealType forwardprob;
RealType backwardprob;
RealType forwardaction;
RealType backwardaction;

RealType tau;

MCWalkerConfiguration& w;
WalkerIter_t repstart, repend;
IndexType direction, headindex, nbeads;
Walker_t* prophead;

inline Reptile(MCWalkerConfiguration& W, WalkerIter_t start, WalkerIter_t end)
: w(W),
repstart(start),
repend(end),
direction(1),
headindex(0),
prophead(0) //, r2prop(0.0), r2accept(0.0),tau(0.0)
{
Action.resize(3);
Action[0] = w.addProperty("ActionBackward");
Action[1] = w.addProperty("ActionForward");
Action[2] = w.addProperty("ActionLocal");
TransProb.resize(2);
TransProb[0] = w.addProperty("TransProbBackward");
TransProb[1] = w.addProperty("TransProbForward");

nbeads = repend - repstart;
}

~Reptile() {}

inline IndexType size() { return nbeads; }

inline Walker_t& operator[](IndexType i) { return getWalker(getBeadIndex(i)); }

inline IndexType wrapIndex(IndexType repindex) { return (repindex % nbeads + nbeads) % nbeads; }

inline Walker_t& getWalker(IndexType i)
{
WalkerIter_t bead = repstart + wrapIndex(i);
return **bead;
}

inline IndexType getBeadIndex(IndexType i) { return wrapIndex(headindex + direction * i); }
inline Walker_t& getBead(IndexType i) { return getWalker(getBeadIndex(i)); }
inline Walker_t& getHead() { return getWalker(getBeadIndex(0)); }
inline Walker_t& getTail() { return getWalker(getBeadIndex(nbeads - 1)); }
inline Walker_t& getNext() { return getWalker(getBeadIndex(nbeads - 2)); }
inline Walker_t& getCenter() { return getWalker(getBeadIndex((nbeads - 1) / 2)); }
//inline void setProposedHead(){

inline void flip()
{
// direction*=-1;
// headindex = getBeadIndex(nbeads-1);
headindex = wrapIndex(headindex - direction);
direction *= -1;
}

inline void setDirection(IndexType dir) { direction = dir; }

inline void setBead(Walker_t& walker, IndexType i)
{
IndexType index = getBeadIndex(i);
Walker_t& newbead(getWalker(index));
newbead = walker; //This should be a hard copy
}

inline void setHead(Walker_t& overwrite)
{
//overwrite last element.
headindex = getBeadIndex(nbeads - 1); //sets to position of tail.
Walker_t& newhead(getBead(0));
newhead = overwrite;
}
//This function does two things: 1.) Moves the reptile forward 1 step. 2.) Returns the new head.
inline Walker_t& getNewHead()
{
//overwrite last element.
headindex = getBeadIndex(nbeads - 1); //sets to position of tail.
return getWalker(headindex);
}

void saveAction(Walker_t& walker, IndexType d, RealType val, IndexType nPsi = 0)
{
//IndexType repdirection=circbuffer.get_direction();
IndexType actionindex = 2;
if (direction != 0)
actionindex = (1 - d * direction) / 2;
walker.Properties(nPsi, Action[actionindex]) = val;
}

RealType getDirectionalAction(Walker_t& walker, IndexType d, IndexType nPsi = 0)
{
//IndexType repdirection=circbuffer.get_direction();
IndexType actionindex = 2;
if (d != 0)
actionindex = (1 - direction * d) / 2;

return walker.Properties(nPsi, Action[actionindex]);
}

RealType getLinkAction(Walker_t& new_walker, Walker_t& old_walker, IndexType d, IndexType nPsi = 0)
{
RealType af = getDirectionalAction(old_walker, +1, nPsi);
RealType ab = getDirectionalAction(new_walker, -1, nPsi);
RealType a0 = getDirectionalAction(old_walker, 0, nPsi) + getDirectionalAction(new_walker, 0, nPsi);
return af + ab + a0;
}

void saveTransProb(Walker_t& walker, IndexType d, RealType val, IndexType nPsi = 0)
{
//IndexType repdirection=circbuffer.get_direction();
IndexType transindex = (1 - d * direction) / 2;
walker.Properties(nPsi, TransProb[transindex]) = val;
}

void saveTransProb(ParticleSet& W, IndexType d, RealType val, IndexType nPsi = 0)
{
//IndexType repdirection=circbuffer.get_direction();
IndexType transindex = (1 - d * direction) / 2;
W.Properties(nPsi, TransProb[transindex]) = val;
}
RealType getTransProb(Walker_t& walker, IndexType d, RealType nPsi = 0)
{
//IndexType repdirection=circbuffer.get_direction();
IndexType transindex = (1 - d * direction) / 2;
return walker.Properties(nPsi, TransProb[transindex]);
}
RealType getTransProb(ParticleSet& W, IndexType d, RealType nPsi = 0)
{
//IndexType repdirection=circbuffer.get_direction();
IndexType transindex = (1 - d * direction) / 2;
return W.Properties(nPsi, TransProb[transindex]);
}

inline void printState()
{
app_log() << "********PRINT REPTILE STATE*********\n";
app_log() << "Direction=" << direction << " Headindex=" << headindex << " tail=" << getBeadIndex(nbeads - 1)
<< "\n next=" << getBeadIndex(nbeads - 2) << " nbeads=" << nbeads << std::endl;
app_log() << "BeadIndex\tWrapIndex\tEnergy\tAction[0]\tAction[1]\tAction[2]\t\n";
for (int i = 0; i < nbeads; i++)
{
app_log() << i << "\t" << getBeadIndex(i) << "\t" << getBead(i).Properties(WP::LOCALENERGY) << "\t"
<< getBead(i).Properties(Action[0]) << "\t" << getBead(i).Properties(Action[1]) << "\t"
<< getBead(i).Properties(Action[2]) << "\n";
}
app_log() << "POSITIONS===============:\n";
for (int i = 0; i < nbeads; i++)
{
// app_log()<<i<<"\t1"<<1<<"\t"<<getBead(i).R[0]<<"\n";
// app_log()<<i<<"\t2"<<2<<"\t"<<getBead(i).R[1]<<"\n";
app_log() << "BEAD #" << i << " tau = " << tau * i << std::endl;
app_log() << getBead(i).R << std::endl;
}
app_log() << "GVECS===============:\n";
for (int i = 0; i < nbeads; i++)
{
// app_log()<<i<<"\t1"<<1<<"\t"<<getBead(i).G[0]<<"\n";
// app_log()<<i<<"\t2"<<2<<"\t"<<getBead(i).G[1]<<"\n";
app_log() << "BEAD #" << i << " tau = " << tau * i << std::endl;
app_log() << getBead(i).G << std::endl;
}
app_log() << "************************************\n";
}
inline RealType getTau() { return tau; }
inline void setTau(RealType t) { tau = t; }


//This takes a value of imaginary time "t" and returns a 3N particle position vector, corresponding to a time slice extrapolated
// from the current reptile. If t>length of reptile, then return the last bead. if t<0; return the first bead.
inline Walker_t::ParticlePos linearInterp(RealType t)
{
IndexType nbead =
IndexType(t / tau); //Calculate the lower bound on the timeslice. t is between binnum*Tau and (binnum+1)Tau
RealType beadfrac = t / tau - nbead; //the fractional coordinate between n and n+1 bead
if (nbead <= 0)
{
ParticleSet::ParticlePos result = getHead().R;
return result;
}
else if (nbead >= nbeads - 1)
{
ParticleSet::ParticlePos result = getTail().R;
return result;
}

else
{
Walker_t::ParticlePos dR(getBead(nbead + 1).R), interpR(getBead(nbead).R);
dR = dR - getBead(nbead).R;

interpR = getBead(nbead).R + beadfrac * dR;
return interpR;
}
}
inline ReptileConfig_t getReptileSlicePositions(RealType tau, RealType beta)
{
IndexType nbeads_new = IndexType(beta / tau);
ReptileConfig_t new_reptile_coords(0);

for (IndexType i = 0; i < nbeads_new; i++)
new_reptile_coords.push_back(linearInterp(tau * i));

return new_reptile_coords;
}

inline void setReptileSlicePositions(ReptileConfig_t& rept)
{
if (rept.size() == nbeads)
{
for (int i = 0; i < nbeads; i++)
getBead(i).R = rept[i];
}
else
;
}

inline void setReptileSlicePositions(Walker_t::ParticlePos R)
{
for (int i = 0; i < nbeads; i++)
getBead(i).R = R;
}
};

using Reptile = ReptileT<QMCTraits::ValueType>;

} // namespace qmcplusplus
#endif
60 changes: 0 additions & 60 deletions src/Particle/SampleStack.cpp
View file

This file was deleted.

44 changes: 3 additions & 41 deletions src/Particle/SampleStack.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,50 +18,12 @@
#ifndef QMCPLUSPLUS_SAMPLE_STACK_H
#define QMCPLUSPLUS_SAMPLE_STACK_H

#include <vector>
#include "Particle/ParticleSet.h"
#include "Particle/MCSample.h"
#include "Particle/Walker.h"
#include "Particle/WalkerConfigurations.h"
#include "Configuration.h"
#include "Particle/SampleStackT.h"

namespace qmcplusplus
{
class SampleStack
{
public:
using PropertySetType = QMCTraits::PropertySetType;

size_t getMaxSamples() const { return max_samples_; }

bool empty() const { return sample_vector_.empty(); }

const MCSample& getSample(size_t i) const;

//@{save/load/clear function for optimization
inline size_t getNumSamples() const { return current_sample_count_; }
///set the number of max samples per rank.
void setMaxSamples(size_t n, size_t number_of_ranks = 1);
/// Global number of samples is number of samples per rank * number of ranks
size_t getGlobalNumSamples() const { return global_num_samples_; }
/// load a single sample from SampleStack
void loadSample(ParticleSet& pset, size_t iw) const;

void appendSample(MCSample&& sample);

///clear the ensemble
void clearEnsemble();
//@}
/// Set the sample count to zero but preserve the storage
void resetSampleCount();

private:
size_t max_samples_{10};
size_t current_sample_count_{0};
size_t global_num_samples_{max_samples_};

std::vector<MCSample> sample_vector_;
};

using SampleStack = SampleStackT<QMCTraits::ValueType>;

} // namespace qmcplusplus
#endif
7 changes: 1 addition & 6 deletions src/QMCDrivers/MCPopulation.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,12 +24,7 @@
#include "QMCDrivers/WalkerElementsRef.h"
#include "OhmmsPETE/OhmmsVector.h"
#include "Utilities/FairDivide.h"

// forward declaration
namespace optimize
{
struct VariableSet;
}
#include "QMCWaveFunctions/VariableSet.h"

namespace qmcplusplus
{
Expand Down
Loading

0 comments on commit 538bfcd

Please sign in to comment.