diff --git a/docs/methods.rst b/docs/methods.rst index 44bd5943c2..bb0913f8cf 100644 --- a/docs/methods.rst +++ b/docs/methods.rst @@ -1758,6 +1758,41 @@ Batched ``dmc`` driver (experimental) - ``spin_mass`` Optional parameter to allow the user to change the rate of spin sampling. If spin sampling is on using ``spinor`` == yes in the electron ParticleSet input, the spin mass determines the rate of spin sampling, resulting in an effective spin timestep :math:`\tau_s = \frac{\tau}{\mu_s}`. The algorithm is described in detail in :cite:`Melton2016-1` and :cite:`Melton2016-2`. +- ``warmupsteps``: These are the steps at the beginning of a DMC run in + which the instantaneous population average energy is used to update the trial + energy and updates happen at every step. The aim is to rapidly equilibrate the population while avoiding overly large population fluctuations. + Unlike VMC, these warmupsteps are included in the requested DMC step count. + +.. math:: + + E_\text{trial} = E_\text{pop\_avg}+(\ln \texttt{targetwalkers}-\ln N_\text{pop}) / \texttt{timestep} + +where :math:`E_\text{pop\_avg}` is the local energy average over the walker population at the current step +and :math:`N_\text{pop}` is the current walker population size. +After the warm-up phase, the trial energy is updated as + +.. math:: + + E_\text{trial} = E_\text{ref}+\texttt{feedback}\cdot(\ln\texttt{targetWalkers}-\ln N_\text{pop}) + +where :math:`E_\text{ref}` is the :math:`E_\text{pop\_avg}` average over all the post warm-up steps up to the current step. The update frequency is controlled by ``energyUpdateInterval``. + +- ``energyUpdateInterval``: Post warm-up, the trial energy is updated every + ``energyUpdateInterval`` steps. Default value is 1 (every step). + +- ``refEnergy``: The default reference energy is taken from the VMC run + that precedes the DMC run. This value is updated to the current mean + whenever branching happens. + +- ``feedback``: This variable is used to determine how strong to react + to population fluctuations when doing population control. Default value is 1. See the + equation in ``warmupsteps`` for more details. + +- ``refenergy_update_scheme``: choose a scheme for updating :math:`E_\text{ref}` used in calculating :math:`E_\text{trial}`. + 'unlimited_history' uses unweighted average of :math:`E_\text{pop\_avg}` of all the steps collected post warm-up. + 'limited_history' uses weighted average of :math:`E_\text{pop\_avg}` of the latest at maximum + min(1, int(1.0 / (feedback * tau))) steps collected post warm-up. Default 'unlimited_history'. + .. code-block:: :caption: The following is an example of a minimal DMC section using the batched ``dmc`` driver :name: Listing 48b diff --git a/src/Estimators/SizeLimitedDataQueue.hpp b/src/Estimators/SizeLimitedDataQueue.hpp new file mode 100644 index 0000000000..3d0bc0a519 --- /dev/null +++ b/src/Estimators/SizeLimitedDataQueue.hpp @@ -0,0 +1,83 @@ +////////////////////////////////////////////////////////////////////////////////////// +// This file is distributed under the University of Illinois/NCSA Open Source License. +// See LICENSE file in top directory for details. +// +// Copyright (c) 2023 QMCPACK developers. +// +// File developed by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory +// +// File created by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory +////////////////////////////////////////////////////////////////////////////////////// +// -*- C++ -*- +#ifndef QMCPLUSPLUS_SIZELIMITEDDATAQUEUE_H +#define QMCPLUSPLUS_SIZELIMITEDDATAQUEUE_H + +#include +#include +#include + +namespace qmcplusplus +{ + +/** collect data with a history limit. + * data stored in std::deque> + */ +template +class SizeLimitedDataQueue +{ +public: + struct HistoryElement + { + T weight; + std::array properties; + }; + + using value_type = HistoryElement; + + SizeLimitedDataQueue(size_t size_limit) : size_limit_(size_limit) {} + + /// add a new record + void push(const value_type& val) + { + if (data.size() == size_limit_) + data.pop_front(); + assert(data.size() < size_limit_); + data.push_back(val); + } + + /// add a new record + void push(value_type&& val) + { + if (data.size() == size_limit_) + data.pop_front(); + assert(data.size() < size_limit_); + data.push_back(val); + } + + /// return weighted average + auto weighted_avg() const + { + std::array avg; + std::fill(avg.begin(), avg.end(), T(0)); + T weight_sum = 0; + for (auto& element : data) + { + weight_sum += element.weight; + for (size_t i = 0; i < NUM_FIELDS; i++) + avg[i] += element.properties[i] * element.weight; + } + for (size_t i = 0; i < NUM_FIELDS; i++) + avg[i] /= weight_sum; + return avg; + } + + /// return the number of records + auto size() const { return data.size(); } + +private: + std::deque data; + const size_t size_limit_; +}; + +} // namespace qmcplusplus +#endif diff --git a/src/Estimators/tests/CMakeLists.txt b/src/Estimators/tests/CMakeLists.txt index f4e7bcce25..b19b2370d0 100644 --- a/src/Estimators/tests/CMakeLists.txt +++ b/src/Estimators/tests/CMakeLists.txt @@ -33,6 +33,7 @@ set(SRCS test_InputSection.cpp test_EstimatorManagerInput.cpp test_ScalarEstimatorInputs.cpp + test_SizeLimitedDataQueue.cpp test_MomentumDistribution.cpp test_OneBodyDensityMatricesInput.cpp test_OneBodyDensityMatrices.cpp diff --git a/src/Estimators/tests/test_SizeLimitedDataQueue.cpp b/src/Estimators/tests/test_SizeLimitedDataQueue.cpp new file mode 100644 index 0000000000..88925cff67 --- /dev/null +++ b/src/Estimators/tests/test_SizeLimitedDataQueue.cpp @@ -0,0 +1,49 @@ +////////////////////////////////////////////////////////////////////////////////////// +// This file is distributed under the University of Illinois/NCSA Open Source License. +// See LICENSE file in top directory for details. +// +// Copyright (c) 2023 QMCPACK developers. +// +// File developed by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory +// +// File created by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory +////////////////////////////////////////////////////////////////////////////////////// + +#include "catch.hpp" +#include "SizeLimitedDataQueue.hpp" + +namespace qmcplusplus +{ + +TEST_CASE("SizeLimitedDataQueue", "[estimators]") +{ + SizeLimitedDataQueue weight_and_energy(3); + CHECK(weight_and_energy.size() == 0); + { + weight_and_energy.push({1.0, {2.0}}); + CHECK(weight_and_energy.size() == 1); + auto avg = weight_and_energy.weighted_avg(); + CHECK(Approx(avg[0]) == 2.0); + } + { + weight_and_energy.push({3.0, {1.0}}); + CHECK(weight_and_energy.size() == 2); + auto avg = weight_and_energy.weighted_avg(); + CHECK(Approx(avg[0]) == 1.25); + } + { + SizeLimitedDataQueue::HistoryElement temp{0.5, {3.0}}; + weight_and_energy.push(std::move(temp)); + CHECK(weight_and_energy.size() == 3); + auto avg = weight_and_energy.weighted_avg(); + CHECK(Approx(avg[0]) == 1.444444444); + } + { + weight_and_energy.push({0.5, {3.0}}); + CHECK(weight_and_energy.size() == 3); + auto avg = weight_and_energy.weighted_avg(); + CHECK(Approx(avg[0]) == 1.5); + } +} + +} // namespace qmcplusplus diff --git a/src/QMCDrivers/CMakeLists.txt b/src/QMCDrivers/CMakeLists.txt index 80cda7e478..5e6ef66f66 100644 --- a/src/QMCDrivers/CMakeLists.txt +++ b/src/QMCDrivers/CMakeLists.txt @@ -64,6 +64,7 @@ set(QMCDRIVERS DMC/DMCFactoryNew.cpp DMC/DMCBatched.cpp DMC/DMCDriverInput.cpp + DMC/DMCRefEnergy.cpp DMC/WalkerControlFactory.cpp DMC/WalkerReconfiguration.cpp DMC/WalkerControl.cpp diff --git a/src/QMCDrivers/DMC/DMCBatched.cpp b/src/QMCDrivers/DMC/DMCBatched.cpp index 6935034035..58b2e4c08c 100644 --- a/src/QMCDrivers/DMC/DMCBatched.cpp +++ b/src/QMCDrivers/DMC/DMCBatched.cpp @@ -394,7 +394,12 @@ void DMCBatched::process(xmlNodePtr node) // I'd like to do away with this method in DMCBatched. app_log() << " Creating the branching engine and walker controler" << std::endl; - branch_engine_ = std::make_unique(qmcdriver_input_.get_tau(), population_.get_num_global_walkers()); + const auto refE_update_scheme = dmcdriver_input_.get_refenergy_update_scheme(); + app_log() << " Reference energy is updated using the " + << (refE_update_scheme == DMCRefEnergyScheme::UNLIMITED_HISTORY ? "unlimited_history" : "limited_history") << " scheme" + << std::endl; + branch_engine_ = + std::make_unique(qmcdriver_input_.get_tau(), dmcdriver_input_.get_feedback(), refE_update_scheme); branch_engine_->put(node); walker_controller_ = std::make_unique(myComm, Random, dmcdriver_input_.get_reconfiguration()); diff --git a/src/QMCDrivers/DMC/DMCDriverInput.cpp b/src/QMCDrivers/DMC/DMCDriverInput.cpp index 820645d960..d954069722 100644 --- a/src/QMCDrivers/DMC/DMCDriverInput.cpp +++ b/src/QMCDrivers/DMC/DMCDriverInput.cpp @@ -19,10 +19,13 @@ void DMCDriverInput::readXML(xmlNodePtr node) { ParameterSet parameter_set_; std::string reconfig_str; + std::string refE_update_scheme_str; parameter_set_.add(reconfig_str, "reconfiguration", {"no", "yes", "runwhileincorrect"}); parameter_set_.add(NonLocalMove, "nonlocalmove", {"no", "yes", "v0", "v1", "v3"}); parameter_set_.add(NonLocalMove, "nonlocalmoves", {"no", "yes", "v0", "v1", "v3"}); parameter_set_.add(max_age_, "MaxAge"); + parameter_set_.add(feedback_, "feedback"); + parameter_set_.add(refE_update_scheme_str, "refenergy_update_scheme", {"unlimited_history", "limited_history"}); // from DMC.cpp put(xmlNodePtr) parameter_set_.add(branch_interval_, "branchInterval"); @@ -64,6 +67,11 @@ void DMCDriverInput::readXML(xmlNodePtr node) if (reserve_ < 1.0) throw std::runtime_error("You can only reserve walkers above the target walker count"); + + if (refE_update_scheme_str == "unlimited_history") + refenergy_update_scheme_ = DMCRefEnergyScheme::UNLIMITED_HISTORY; + else + refenergy_update_scheme_ = DMCRefEnergyScheme::LIMITED_HISTORY; } std::ostream& operator<<(std::ostream& o_stream, const DMCDriverInput& dmci) { return o_stream; } diff --git a/src/QMCDrivers/DMC/DMCDriverInput.h b/src/QMCDrivers/DMC/DMCDriverInput.h index f357e6da41..d729d48ead 100644 --- a/src/QMCDrivers/DMC/DMCDriverInput.h +++ b/src/QMCDrivers/DMC/DMCDriverInput.h @@ -14,6 +14,7 @@ #include "Configuration.h" #include "OhmmsData/ParameterSet.h" +#include "DMC/DMCRefEnergyScheme.h" namespace qmcplusplus { @@ -32,6 +33,8 @@ class DMCDriverInput bool get_reconfiguration() const { return reconfiguration_; } IndexType get_max_age() const { return max_age_; } IndexType get_branch_interval() const { return branch_interval_; } + double get_feedback() const { return feedback_; } + DMCRefEnergyScheme get_refenergy_update_scheme() const { return refenergy_update_scheme_; } const std::string& get_non_local_move() const { return NonLocalMove; } double get_alpha() const { return alpha_; } double get_gamma() const { return gamma_; } @@ -46,6 +49,10 @@ class DMCDriverInput */ ///Interval between branching IndexType branch_interval_ = 1; + ///feed back parameter for population control + double feedback_ = 1.0; + ///input std::string to determine reference energy update scheme + DMCRefEnergyScheme refenergy_update_scheme_; ///input std::string to determine kill walkers or not std::string KillWalker; ///input std::string to determine swap walkers among mpi processors diff --git a/src/QMCDrivers/DMC/DMCRefEnergy.cpp b/src/QMCDrivers/DMC/DMCRefEnergy.cpp new file mode 100644 index 0000000000..0ae22692ac --- /dev/null +++ b/src/QMCDrivers/DMC/DMCRefEnergy.cpp @@ -0,0 +1,57 @@ +////////////////////////////////////////////////////////////////////////////////////// +// This file is distributed under the University of Illinois/NCSA Open Source License. +// See LICENSE file in top directory for details. +// +// Copyright (c) 2023 QMCPACK developers. +// +// File developed by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory +// +// File created by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory +////////////////////////////////////////////////////////////////////////////////////// + +#include +#include "DMCRefEnergy.h" + +namespace qmcplusplus +{ + +using FullPrecReal = DMCRefEnergy::FullPrecReal; + +DMCRefEnergy::DMCRefEnergy(DMCRefEnergyScheme scheme, size_t history_limit) + : scheme_(scheme), energy_and_variance_(history_limit) +{} + +std::tuple DMCRefEnergy::getEnergyVariance() const +{ + if (scheme_ == DMCRefEnergyScheme::LIMITED_HISTORY) + { + auto avg = energy_and_variance_.weighted_avg(); + return {avg[ENERGY], avg[VARIANCE]}; + } + else + return {energy_hist_.mean(), variance_hist_.mean()}; +} + +void DMCRefEnergy::pushWeightEnergyVariance(FullPrecReal weight, FullPrecReal ene, FullPrecReal var) +{ + if (scheme_ == DMCRefEnergyScheme::LIMITED_HISTORY) + energy_and_variance_.push({weight, {ene, var}}); + else + { + energy_hist_(ene); + variance_hist_(var); + } +} + +size_t DMCRefEnergy::count() const +{ + if (scheme_ == DMCRefEnergyScheme::LIMITED_HISTORY) + return energy_and_variance_.size(); + else + { + assert(energy_hist_.count() == variance_hist_.count()); + return energy_hist_.count(); + } +} + +} // namespace qmcplusplus diff --git a/src/QMCDrivers/DMC/DMCRefEnergy.h b/src/QMCDrivers/DMC/DMCRefEnergy.h new file mode 100644 index 0000000000..cb122ea508 --- /dev/null +++ b/src/QMCDrivers/DMC/DMCRefEnergy.h @@ -0,0 +1,64 @@ +////////////////////////////////////////////////////////////////////////////////////// +// This file is distributed under the University of Illinois/NCSA Open Source License. +// See LICENSE file in top directory for details. +// +// Copyright (c) 2023 QMCPACK developers. +// +// File developed by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory +// +// File created by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory +////////////////////////////////////////////////////////////////////////////////////// +// -*- C++ -*- +#ifndef QMCPLUSPLUS_DMCREFENERGY_H +#define QMCPLUSPLUS_DMCREFENERGY_H + +#include +#include +#include +#include +#include "DMCRefEnergyScheme.h" + +namespace qmcplusplus +{ +/** Handle updating Eref used for calculating the trial energy. + */ +class DMCRefEnergy +{ +public: + using FullPrecReal = QMCTraits::FullPrecRealType; + + enum DataLayout + { + ENERGY = 0, + VARIANCE, + DATA_SIZE + }; + +private: + /// scheme + DMCRefEnergyScheme scheme_; + + // legacy scheme data + ///a simple accumulator for energy + accumulator_set energy_hist_; + ///a simple accumulator for variance + accumulator_set variance_hist_; + + // limited memory scheme data + SizeLimitedDataQueue energy_and_variance_; + +public: + DMCRefEnergy(DMCRefEnergyScheme scheme, size_t history_limit); + + /// return energy and variance + std::tuple getEnergyVariance() const; + + /// record weight, energy and variance. + void pushWeightEnergyVariance(FullPrecReal weight, FullPrecReal ene, FullPrecReal var); + + /// return record count. + size_t count() const; +}; + +} // namespace qmcplusplus +#endif diff --git a/src/QMCDrivers/DMC/DMCRefEnergyScheme.h b/src/QMCDrivers/DMC/DMCRefEnergyScheme.h new file mode 100644 index 0000000000..caa4d8f339 --- /dev/null +++ b/src/QMCDrivers/DMC/DMCRefEnergyScheme.h @@ -0,0 +1,25 @@ +////////////////////////////////////////////////////////////////////////////////////// +// This file is distributed under the University of Illinois/NCSA Open Source License. +// See LICENSE file in top directory for details. +// +// Copyright (c) 2023 QMCPACK developers. +// +// File developed by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory +// +// File created by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory +////////////////////////////////////////////////////////////////////////////////////// +// -*- C++ -*- +#ifndef QMCPLUSPLUS_DMCREFENERGYSCHEME_H +#define QMCPLUSPLUS_DMCREFENERGYSCHEME_H + +namespace qmcplusplus +{ +/** DMCRefEnergy schemes + */ +enum class DMCRefEnergyScheme +{ + UNLIMITED_HISTORY, + LIMITED_HISTORY +}; +} // namespace qmcplusplus +#endif diff --git a/src/QMCDrivers/SFNBranch.cpp b/src/QMCDrivers/SFNBranch.cpp index 6696827d5f..1860030397 100644 --- a/src/QMCDrivers/SFNBranch.cpp +++ b/src/QMCDrivers/SFNBranch.cpp @@ -31,7 +31,11 @@ enum DUMMYOPT }; -SFNBranch::SFNBranch(RealType tau, int nideal) : WarmUpToDoSteps(0), EtrialUpdateToDoSteps(0), myNode(NULL) +SFNBranch::SFNBranch(RealType tau, RealType feedback, DMCRefEnergyScheme refenergy_update_scheme) + : WarmUpToDoSteps(0), + EtrialUpdateToDoSteps(0), + myNode(NULL), + ref_energy_collector(refenergy_update_scheme, std::max(1, static_cast(1.0 / (feedback * tau)))) { BranchMode.set(B_DMCSTAGE, 0); //warmup stage BranchMode.set(B_POPCONTROL, 1); //use standard DMC @@ -41,7 +45,7 @@ SFNBranch::SFNBranch(RealType tau, int nideal) : WarmUpToDoSteps(0), EtrialUpdat vParam.fill(1.0); vParam[SBVP::TAU] = tau; vParam[SBVP::TAUEFF] = tau; - vParam[SBVP::FEEDBACK] = 1.0; + vParam[SBVP::FEEDBACK] = feedback; vParam[SBVP::FILTERSCALE] = 10; vParam[SBVP::SIGMA_BOUND] = 10; R2Accepted(1.0e-10); @@ -80,8 +84,6 @@ void SFNBranch::registerParameters() m_param.add(vParam[SBVP::TAU], "TimeStep"); //filterscale: sets the filtercutoff to sigma*filterscale m_param.add(vParam[SBVP::FILTERSCALE], "filterscale"); - //feed back parameter for population control - m_param.add(vParam[SBVP::FEEDBACK], "feedback"); m_param.add(vParam[SBVP::SIGMA_BOUND], "sigmaBound"); //turn on/off effective tau onl for time-step error comparisons m_param.add(sParam[USETAUOPT], "useBareTau"); @@ -106,8 +108,7 @@ int SFNBranch::initParam(const MCPopulation& population, vParam[SBVP::SIGMA2] = var; vParam[SBVP::TAUEFF] = vParam[SBVP::TAU] * R2Accepted.result() / R2Proposed.result(); /// FIXME, magic number 50 - setBranchCutoff(vParam[SBVP::SIGMA2], vParam[SBVP::SIGMA_BOUND], 50, - population.get_golden_electrons().getTotalNum()); + setBranchCutoff(vParam[SBVP::SIGMA2], vParam[SBVP::SIGMA_BOUND], 50, population.get_golden_electrons().getTotalNum()); int nwtot_now = population.get_num_global_walkers(); if (iParam[B_TARGETWALKERS] == 0) @@ -149,12 +150,13 @@ void SFNBranch::updateParamAfterPopControl(int pop_int, const MCDataType namespace qmcplusplus @@ -148,7 +149,7 @@ class SFNBranch : public QMCTraits using VParamType = VParams; ///Constructor - SFNBranch(RealType tau, int nideal); + SFNBranch(RealType tau, RealType feedback, DMCRefEnergyScheme); ///copy constructor SFNBranch(const SFNBranch& abranch) = delete; @@ -305,10 +306,8 @@ class SFNBranch : public QMCTraits int EtrialUpdateToDoSteps; ///save xml element xmlNodePtr myNode; - ///a simple accumulator for energy - accumulator_set EnergyHist; - ///a simple accumulator for variance - accumulator_set VarianceHist; + /// collect energy and variance history + DMCRefEnergy ref_energy_collector; ///a simple accumulator for energy accumulator_set R2Accepted; ///a simple accumulator for energy diff --git a/src/QMCDrivers/tests/test_SFNBranch.cpp b/src/QMCDrivers/tests/test_SFNBranch.cpp index 292ec54e81..d2ccdc1d5d 100644 --- a/src/QMCDrivers/tests/test_SFNBranch.cpp +++ b/src/QMCDrivers/tests/test_SFNBranch.cpp @@ -49,7 +49,7 @@ class SetupSFNBranch walkers[0].get().R[0] = 1.0; walkers[1].get().R[0] = 0.5; - auto sfnb = std::make_unique(tau_, num_global_walkers_); + auto sfnb = std::make_unique(tau_, 1.0, DMCRefEnergyScheme::LIMITED_HISTORY); createMyNode(*sfnb, valid_dmc_input_sections[valid_dmc_input_dmc_batch_index]);