diff --git a/src/Particle/CMakeLists.txt b/src/Particle/CMakeLists.txt index b7055795f4..4f47940a8f 100644 --- a/src/Particle/CMakeLists.txt +++ b/src/Particle/CMakeLists.txt @@ -26,7 +26,6 @@ set(PARTICLE MCWalkerConfigurationT.cpp WalkerConfigurations.cpp SpeciesSet.cpp - SampleStack.cpp SampleStackT.cpp createDistanceTableAA.cpp createDistanceTableAB.cpp diff --git a/src/Particle/Reptile.h b/src/Particle/Reptile.h index 6e7f7cd7cd..d32143a4ca 100644 --- a/src/Particle/Reptile.h +++ b/src/Particle/Reptile.h @@ -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; - - std::vector Action; - std::vector 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()<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; } // namespace qmcplusplus #endif diff --git a/src/Particle/SampleStack.cpp b/src/Particle/SampleStack.cpp deleted file mode 100644 index c2720b9b8c..0000000000 --- a/src/Particle/SampleStack.cpp +++ /dev/null @@ -1,60 +0,0 @@ -////////////////////////////////////////////////////////////////////////////////////// -// This file is distributed under the University of Illinois/NCSA Open Source License. -// See LICENSE file in top directory for details. -// -// Copyright (c) 2020 QMCPACK developers. -// -// File developed by: Mark Dewing, mdewing@anl.gov, Argonne National Laboratory -// -// File created by: Mark Dewing, mdewing@anl.gov, Argonne National Laboratory -////////////////////////////////////////////////////////////////////////////////////// - - -#include "SampleStack.h" -#include "Utilities/IteratorUtility.h" - -namespace qmcplusplus -{ - -/** allocate the SampleStack - * @param n number of samples per rank - * @param num_ranks number of ranks. Used to set global number of samples. - */ -void SampleStack::setMaxSamples(size_t n, size_t num_ranks) -{ - max_samples_ = n; - global_num_samples_ = n * num_ranks; - current_sample_count_ = std::min(current_sample_count_, max_samples_); - sample_vector_.resize(n, MCSample(0)); -} - -const MCSample& SampleStack::getSample(size_t i) const { return sample_vector_[i]; } - -void SampleStack::appendSample(MCSample&& sample) -{ - // Ignore samples in excess of the expected number of samples - if (current_sample_count_ < max_samples_) - { - sample_vector_[current_sample_count_] = std::move(sample); - current_sample_count_++; - } -} - -/** load a single sample from SampleStack - */ -void SampleStack::loadSample(ParticleSet& pset, size_t iw) const -{ - pset.R = sample_vector_[iw].R; - pset.spins = sample_vector_[iw].spins; -} - -void SampleStack::clearEnsemble() -{ - sample_vector_.clear(); - current_sample_count_ = 0; -} - -void SampleStack::resetSampleCount() { current_sample_count_ = 0; } - - -} // namespace qmcplusplus diff --git a/src/Particle/SampleStack.h b/src/Particle/SampleStack.h index 3614f53558..9565ceac02 100644 --- a/src/Particle/SampleStack.h +++ b/src/Particle/SampleStack.h @@ -18,50 +18,12 @@ #ifndef QMCPLUSPLUS_SAMPLE_STACK_H #define QMCPLUSPLUS_SAMPLE_STACK_H -#include -#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 sample_vector_; -}; - +using SampleStack = SampleStackT; } // namespace qmcplusplus #endif diff --git a/src/QMCDrivers/MCPopulation.h b/src/QMCDrivers/MCPopulation.h index 1d876030be..fd38c242d0 100644 --- a/src/QMCDrivers/MCPopulation.h +++ b/src/QMCDrivers/MCPopulation.h @@ -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 { diff --git a/src/QMCDrivers/VMC/VMCBatched.h b/src/QMCDrivers/VMC/VMCBatched.h index 589e6ee6da..c404510167 100644 --- a/src/QMCDrivers/VMC/VMCBatched.h +++ b/src/QMCDrivers/VMC/VMCBatched.h @@ -18,7 +18,7 @@ #include "QMCDrivers/MCPopulation.h" #include "QMCDrivers/ContextForSteps.h" #include "QMCDrivers/GreenFunctionModifiers/DriftModifierBase.h" - +#include "Particle/SampleStack.h" #include "Utilities/Timer.h" namespace qmcplusplus diff --git a/src/QMCDrivers/VMC/VMCFactoryNew.h b/src/QMCDrivers/VMC/VMCFactoryNew.h index 0dc19ff9a8..e1247a5d5a 100644 --- a/src/QMCDrivers/VMC/VMCFactoryNew.h +++ b/src/QMCDrivers/VMC/VMCFactoryNew.h @@ -19,6 +19,7 @@ #include "QMCWaveFunctions/WaveFunctionPool.h" #include "Message/Communicate.h" #include "Particle/ParticleSetPool.h" +#include "Particle/SampleStack.h" namespace qmcplusplus { diff --git a/src/QMCDrivers/WFOpt/QMCFixedSampleLinearOptimize.cpp b/src/QMCDrivers/WFOpt/QMCFixedSampleLinearOptimize.cpp index d7e0dae882..e45b2cebff 100644 --- a/src/QMCDrivers/WFOpt/QMCFixedSampleLinearOptimize.cpp +++ b/src/QMCDrivers/WFOpt/QMCFixedSampleLinearOptimize.cpp @@ -28,6 +28,7 @@ #include "CPU/Blasf.h" #include "Numerics/MatrixOperators.h" #include "Message/UniformCommunicateError.h" +#include "Particle/SampleStack.h" #include #ifdef HAVE_LMY_ENGINE #include "formic/utils/matrix.h" diff --git a/src/QMCDrivers/WFOpt/QMCFixedSampleLinearOptimizeBatched.cpp b/src/QMCDrivers/WFOpt/QMCFixedSampleLinearOptimizeBatched.cpp index 7aaac1e60b..fa4c35012c 100644 --- a/src/QMCDrivers/WFOpt/QMCFixedSampleLinearOptimizeBatched.cpp +++ b/src/QMCDrivers/WFOpt/QMCFixedSampleLinearOptimizeBatched.cpp @@ -30,6 +30,8 @@ #include "Numerics/MatrixOperators.h" #include "EstimatorInputDelegates.h" #include "Message/UniformCommunicateError.h" +#include "Particle/SampleStack.h" + #include #ifdef HAVE_LMY_ENGINE #include "formic/utils/matrix.h" diff --git a/src/QMCDrivers/WFOpt/QMCWFOptFactoryNew.h b/src/QMCDrivers/WFOpt/QMCWFOptFactoryNew.h index d78ffb5ff4..010ef554ad 100644 --- a/src/QMCDrivers/WFOpt/QMCWFOptFactoryNew.h +++ b/src/QMCDrivers/WFOpt/QMCWFOptFactoryNew.h @@ -13,6 +13,7 @@ #define QMCPLUSPLUS_WFOPTFACTORYNEW_H #include "QMCDrivers/QMCDriverInterface.h" +#include "Particle/SampleStack.h" class Communicate; @@ -22,7 +23,6 @@ class MCPopulation; class WaveFunctionPool; class QMCHamiltonian; class TrialWaveFunction; -class SampleStack; class QMCFixedSampleLinearOptimizeBatched; class ProjectData;