From 4c73f010c13e10e5d7569f375985c8c849fdf8a1 Mon Sep 17 00:00:00 2001 From: Philip Fackler Date: Fri, 29 Sep 2023 11:25:36 -0400 Subject: [PATCH] Implemented WalkerConfigurationsT and dealt with the ensuing mess --- src/Estimators/EstimatorManagerBase.h | 2 +- src/Estimators/EstimatorManagerCrowd.h | 3 +- src/Estimators/OperatorEstBase.h | 3 +- src/Estimators/ScalarEstimatorBase.h | 2 +- .../tests/test_EstimatorManagerCrowd.cpp | 2 +- src/Particle/CMakeLists.txt | 2 +- src/Particle/InitMolecularSystemT.cpp | 2 +- src/Particle/MCWalkerConfigurationT.cpp | 48 +- src/Particle/MCWalkerConfigurationT.h | 14 +- src/Particle/ParticleSetPoolT.cpp | 2 +- src/Particle/VirtualParticleSetT.cpp | 3 + src/Particle/WalkerConfigurations.cpp | 149 --- src/Particle/WalkerConfigurations.h | 163 +--- src/Particle/WalkerConfigurationsT.cpp | 170 ++++ src/Particle/WalkerConfigurationsT.h | 258 ++++++ src/QMCDrivers/MCPopulation.h | 2 +- src/QMCDrivers/WFOpt/QMCCostFunctionBatched.h | 1 + src/QMCDrivers/WalkerElementsRef.h | 42 +- src/QMCDrivers/tests/WalkerConsumer.h | 68 +- src/QMCDrivers/tests/test_Crowd.cpp | 2 +- src/QMCDrivers/tests/test_DMCBatched.cpp | 1 + src/QMCDrivers/tests/test_SFNBranch.cpp | 2 +- .../BsplineFactory/SplineSetReaderT.h | 1 + .../BsplineFactory/createBsplineReaderT.cpp | 10 + src/QMCWaveFunctions/CompositeSPOSetT.cpp | 3 + src/QMCWaveFunctions/EinsplineSetBuilderT.cpp | 5 + .../HarmonicOscillator/SHOSetBuilderT.cpp | 3 + .../HarmonicOscillator/SHOSetT.cpp | 3 + .../LCAO/LCAOrbitalBuilderT.cpp | 3 + src/QMCWaveFunctions/LCAO/LCAOrbitalSetT.cpp | 3 + .../LCAO/SoaLocalizedBasisSetT.cpp | 67 +- .../SPOSetBuilderFactoryT.cpp | 3 + src/QMCWaveFunctions/VariableSetT.cpp | 8 +- src/QMCWaveFunctions/tests/CMakeLists.txt | 7 - .../tests/test_RotatedSPOs.cpp | 869 ------------------ .../tests/test_RotatedSPOsT.cpp | 71 +- 36 files changed, 668 insertions(+), 1329 deletions(-) delete mode 100644 src/Particle/WalkerConfigurations.cpp create mode 100644 src/Particle/WalkerConfigurationsT.cpp create mode 100644 src/Particle/WalkerConfigurationsT.h delete mode 100644 src/QMCWaveFunctions/tests/test_RotatedSPOs.cpp diff --git a/src/Estimators/EstimatorManagerBase.h b/src/Estimators/EstimatorManagerBase.h index be4f912fa2..57f43bc9a1 100644 --- a/src/Estimators/EstimatorManagerBase.h +++ b/src/Estimators/EstimatorManagerBase.h @@ -52,7 +52,7 @@ class EstimatorManagerBase using EstimatorType = ScalarEstimatorBase; using BufferType = std::vector; - using MCPWalker = Walker; + using MCPWalker = MCWalkerConfiguration::Walker_t; ///default constructor EstimatorManagerBase(Communicate* c = 0); diff --git a/src/Estimators/EstimatorManagerCrowd.h b/src/Estimators/EstimatorManagerCrowd.h index e97bc9e8f6..3d845c827b 100644 --- a/src/Estimators/EstimatorManagerCrowd.h +++ b/src/Estimators/EstimatorManagerCrowd.h @@ -22,6 +22,7 @@ #include "Estimators/EstimatorManagerNew.h" #include "Particle/Walker.h" #include "OhmmsPETE/OhmmsVector.h" +#include "Particle/MCWalkerConfiguration.h" namespace qmcplusplus { @@ -38,7 +39,7 @@ class QMCHamiltonian; class EstimatorManagerCrowd { public: - using MCPWalker = Walker; + using MCPWalker = MCWalkerConfiguration::Walker_t; using RealType = EstimatorManagerNew::RealType; using FullPrecRealType = EstimatorManagerNew::FullPrecRealType; diff --git a/src/Estimators/OperatorEstBase.h b/src/Estimators/OperatorEstBase.h index 49244e7b90..489ca4ed29 100644 --- a/src/Estimators/OperatorEstBase.h +++ b/src/Estimators/OperatorEstBase.h @@ -23,6 +23,7 @@ #include "QMCWaveFunctions/OrbitalSetTraits.h" #include "type_traits/DataLocality.h" #include "hdf/hdf_archive.h" +#include "Particle/MCWalkerConfiguration.h" #include namespace qmcplusplus @@ -41,7 +42,7 @@ class OperatorEstBase public: using QMCT = QMCTraits; using FullPrecRealType = QMCT::FullPrecRealType; - using MCPWalker = Walker; + using MCPWalker = MCWalkerConfiguration::Walker_t; using Data = std::vector; diff --git a/src/Estimators/ScalarEstimatorBase.h b/src/Estimators/ScalarEstimatorBase.h index 9ccaae2dc7..c848aee332 100644 --- a/src/Estimators/ScalarEstimatorBase.h +++ b/src/Estimators/ScalarEstimatorBase.h @@ -42,7 +42,7 @@ struct ScalarEstimatorBase using RealType = QMCTraits::FullPrecRealType; using accumulator_type = accumulator_set; using Walker_t = MCWalkerConfiguration::Walker_t; - using MCPWalker = Walker; + using MCPWalker = Walker_t; using WalkerIterator = MCWalkerConfiguration::const_iterator; using RecordListType = RecordNamedProperty; diff --git a/src/Estimators/tests/test_EstimatorManagerCrowd.cpp b/src/Estimators/tests/test_EstimatorManagerCrowd.cpp index b4c0548572..e2eb6fcecc 100644 --- a/src/Estimators/tests/test_EstimatorManagerCrowd.cpp +++ b/src/Estimators/tests/test_EstimatorManagerCrowd.cpp @@ -105,7 +105,7 @@ TEST_CASE("EstimatorManagerCrowd PerParticleHamiltonianLogger integration", "[es EstimatorManagerCrowd emc(emn); - using MCPWalker = Walker; + using MCPWalker = EstimatorManagerCrowd::MCPWalker; std::vector walkers(num_walkers, MCPWalker(pset.getTotalNum())); diff --git a/src/Particle/CMakeLists.txt b/src/Particle/CMakeLists.txt index 4f47940a8f..33e890f882 100644 --- a/src/Particle/CMakeLists.txt +++ b/src/Particle/CMakeLists.txt @@ -24,7 +24,7 @@ set(PARTICLE DynamicCoordinatesT.cpp MCCoordsT.cpp MCWalkerConfigurationT.cpp - WalkerConfigurations.cpp + WalkerConfigurationsT.cpp SpeciesSet.cpp SampleStackT.cpp createDistanceTableAA.cpp diff --git a/src/Particle/InitMolecularSystemT.cpp b/src/Particle/InitMolecularSystemT.cpp index 0fbe90382d..7584dc65ac 100644 --- a/src/Particle/InitMolecularSystemT.cpp +++ b/src/Particle/InitMolecularSystemT.cpp @@ -308,7 +308,7 @@ InitMolecularSystemT::reset() #ifndef QMC_COMPLEX template class InitMolecularSystemT; -template class InitMolecularSystemT; +// template class InitMolecularSystemT; #else template class InitMolecularSystemT>; template class InitMolecularSystemT>; diff --git a/src/Particle/MCWalkerConfigurationT.cpp b/src/Particle/MCWalkerConfigurationT.cpp index 2692bd1b73..753fe89fb7 100644 --- a/src/Particle/MCWalkerConfigurationT.cpp +++ b/src/Particle/MCWalkerConfigurationT.cpp @@ -59,7 +59,7 @@ MCWalkerConfigurationT::MCWalkerConfigurationT( { samples.clearEnsemble(); samples.setMaxSamples(mcw.getMaxSamples()); - setWalkerOffsets(mcw.getWalkerOffsets()); + this->setWalkerOffsets(mcw.getWalkerOffsets()); this->Properties = mcw.Properties; } @@ -70,11 +70,11 @@ template void MCWalkerConfigurationT::createWalkers(int n) { - const int old_nw = getActiveWalkers(); - WalkerConfigurations::createWalkers(n, this->TotalNum); + const int old_nw = this->getActiveWalkers(); + WalkerConfigurationsT::createWalkers(n, this->TotalNum); // no pre-existing walkers, need to initialized based on particleset. if (old_nw == 0) - for (auto& awalker : walker_list_) { + for (auto& awalker : this->walker_list_) { awalker->R = this->R; awalker->spins = this->spins; } @@ -85,16 +85,16 @@ template void MCWalkerConfigurationT::resize(int numWalkers, int numPtcls) { - if (this->TotalNum && walker_list_.size()) + if (this->TotalNum && this->walker_list_.size()) app_warning() << "MCWalkerConfiguration::resize cleans up the walker list." << std::endl; - const int old_nw = getActiveWalkers(); + const int old_nw = this->getActiveWalkers(); ParticleSetT::resize(unsigned(numPtcls)); - WalkerConfigurations::resize(numWalkers, this->TotalNum); + WalkerConfigurationsT::resize(numWalkers, this->TotalNum); // no pre-existing walkers, need to initialized based on particleset. if (old_nw == 0) - for (auto& awalker : walker_list_) { + for (auto& awalker : this->walker_list_) { awalker->R = this->R; awalker->spins = this->spins; } @@ -139,7 +139,7 @@ MCWalkerConfigurationT::resetWalkerProperty(int ncopy) APP_ABORT("Fatal Exception"); } - for (auto& walker : walker_list_) { + for (auto& walker : this->walker_list_) { walker->resizeProperty(ncopy, m); walker->Weight = 1.0; } @@ -153,12 +153,12 @@ MCWalkerConfigurationT::resizeWalkerHistories() // using std::vector > is too costly. int np = this->PropertyHistory.size(); if (np) - for (int iw = 0; iw < walker_list_.size(); ++iw) - walker_list_[iw]->PropertyHistory = this->PropertyHistory; + for (int iw = 0; iw < this->walker_list_.size(); ++iw) + this->walker_list_[iw]->PropertyHistory = this->PropertyHistory; np = this->PHindex.size(); if (np) - for (int iw = 0; iw < walker_list_.size(); ++iw) - walker_list_[iw]->PHindex = this->PHindex; + for (int iw = 0; iw < this->walker_list_.size(); ++iw) + this->walker_list_[iw]->PHindex = this->PHindex; ; } @@ -179,7 +179,7 @@ template void MCWalkerConfigurationT::saveEnsemble() { - saveEnsemble(walker_list_.begin(), walker_list_.end()); + saveEnsemble(this->walker_list_.begin(), this->walker_list_.end()); } /** save the [first,last) walkers to SampleStack @@ -211,14 +211,14 @@ MCWalkerConfigurationT::loadEnsemble() int nsamples = std::min(samples.getMaxSamples(), samples.getNumSamples()); if (samples.empty() || nsamples == 0) return; - Walker_t::PropertyContainer_t prop( + typename Walker_t::PropertyContainer_t prop( 1, this->PropertyList.size(), 1, WP::MAXPROPERTIES); - walker_list_.resize(nsamples); + this->walker_list_.resize(nsamples); for (int i = 0; i < nsamples; ++i) { auto awalker = std::make_unique(this->TotalNum); awalker->Properties.copy(prop); samples.getSample(i).convertToWalker(*awalker); - walker_list_[i] = std::move(awalker); + this->walker_list_[i] = std::move(awalker); } resizeWalkerHistories(); samples.clearEnsemble(); @@ -230,7 +230,7 @@ MCWalkerConfigurationT::dumpEnsemble( std::vector*>& others, HDFWalkerOutput& out, int np, int nBlock) { - WalkerConfigurations wctemp; + WalkerConfigurationsT wctemp; for (auto* mcwc : others) { const auto& astack(mcwc->getSampleStack()); const size_t sample_size = @@ -277,18 +277,18 @@ MCWalkerConfigurationT::loadEnsemble( } int nw_tot = off.back(); if (nw_tot) { - Walker_t::PropertyContainer_t prop( + typename Walker_t::PropertyContainer_t prop( 1, this->PropertyList.size(), 1, WP::MAXPROPERTIES); - while (walker_list_.size()) - pop_back(); - walker_list_.resize(nw_tot); + while (this->walker_list_.size()) + this->pop_back(); + this->walker_list_.resize(nw_tot); for (int i = 0; i < others.size(); ++i) { SampleStackT& astack(others[i]->getSampleStack()); for (int j = 0, iw = off[i]; iw < off[i + 1]; ++j, ++iw) { auto awalker = std::make_unique(this->TotalNum); awalker->Properties.copy(prop); astack.getSample(j).convertToWalker(*awalker); - walker_list_[iw] = std::move(awalker); + this->walker_list_[iw] = std::move(awalker); } if (doclean) others[i]->clearEnsemble(); @@ -307,7 +307,7 @@ MCWalkerConfigurationT::clearEnsemble() #ifndef QMC_COMPLEX template class MCWalkerConfigurationT; -template class MCWalkerConfigurationT; +// template class MCWalkerConfigurationT; #else template class MCWalkerConfigurationT>; template class MCWalkerConfigurationT>; diff --git a/src/Particle/MCWalkerConfigurationT.h b/src/Particle/MCWalkerConfigurationT.h index 49a159e51d..4de261c12a 100644 --- a/src/Particle/MCWalkerConfigurationT.h +++ b/src/Particle/MCWalkerConfigurationT.h @@ -29,7 +29,7 @@ #include "Particle/ParticleSetT.h" #include "Particle/SampleStackT.h" #include "Particle/Walker.h" -#include "Particle/WalkerConfigurations.h" +#include "Particle/WalkerConfigurationsT.h" #include "Utilities/IteratorUtility.h" namespace qmcplusplus @@ -60,7 +60,7 @@ class ReptileT; template class MCWalkerConfigurationT : public ParticleSetT, - public WalkerConfigurations + public WalkerConfigurationsT { public: /**enumeration for update*/ @@ -71,16 +71,16 @@ class MCWalkerConfigurationT : Update_Particle /// move a particle by particle }; - using Walker_t = WalkerConfigurations::Walker_t; + using Walker_t = typename WalkerConfigurationsT::Walker_t; /// container type of the Properties of a Walker - using PropertyContainer_t = Walker_t::PropertyContainer_t; + using PropertyContainer_t = typename Walker_t::PropertyContainer_t; /// container type of Walkers using WalkerList_t = std::vector>; /// FIX: a type alias of iterator for an object should not be for just one /// of many objects it holds. - using iterator = WalkerList_t::iterator; + using iterator = typename WalkerList_t::iterator; /// const_iterator of Walker container - using const_iterator = WalkerList_t::const_iterator; + using const_iterator = typename WalkerList_t::const_iterator; using ReptileList_t = UPtrVector>; @@ -104,7 +104,7 @@ class MCWalkerConfigurationT : resize(int numWalkers, int numPtcls); /// clean up the walker list - using WalkerConfigurations::clear; + using WalkerConfigurationsT::clear; /// resize Walker::PropertyHistory and Walker::PHindex: void resizeWalkerHistories(); diff --git a/src/Particle/ParticleSetPoolT.cpp b/src/Particle/ParticleSetPoolT.cpp index 6e06aa45fc..cea8fecca6 100644 --- a/src/Particle/ParticleSetPoolT.cpp +++ b/src/Particle/ParticleSetPoolT.cpp @@ -272,7 +272,7 @@ ParticleSetPoolT::reset() // explicit instantiations #ifndef QMC_COMPLEX template class ParticleSetPoolT; -template class ParticleSetPoolT; +// template class ParticleSetPoolT; #else template class ParticleSetPoolT>; template class ParticleSetPoolT>; diff --git a/src/Particle/VirtualParticleSetT.cpp b/src/Particle/VirtualParticleSetT.cpp index 1f896405fc..e208c7f8db 100644 --- a/src/Particle/VirtualParticleSetT.cpp +++ b/src/Particle/VirtualParticleSetT.cpp @@ -265,8 +265,11 @@ VirtualParticleSetT::mw_makeMovesWithSpin( ParticleSetT::mw_update(p_list); } +#ifndef QMC_COMPLEX template class VirtualParticleSetT; template class VirtualParticleSetT; +#else template class VirtualParticleSetT>; template class VirtualParticleSetT>; +#endif } // namespace qmcplusplus diff --git a/src/Particle/WalkerConfigurations.cpp b/src/Particle/WalkerConfigurations.cpp deleted file mode 100644 index a3959d1610..0000000000 --- a/src/Particle/WalkerConfigurations.cpp +++ /dev/null @@ -1,149 +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: Jordan E. Vincent, University of Illinois at Urbana-Champaign -// Bryan Clark, bclark@Princeton.edu, Princeton University -// Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign -// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign -// Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign -// Cynthia Gu, zg1@ornl.gov, Oak Ridge National Laboratory -// Ye Luo, yeluo@anl.gov, Argonne National Laboratory -// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory -// -// File created by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory -////////////////////////////////////////////////////////////////////////////////////// - - -#include "WalkerConfigurations.h" -#include -#include "Utilities/IteratorUtility.h" - -namespace qmcplusplus -{ -WalkerConfigurations::WalkerConfigurations() = default; - -///default destructor -WalkerConfigurations::~WalkerConfigurations() { destroyWalkers(walker_list_.begin(), walker_list_.end()); } - - -void WalkerConfigurations::createWalkers(int n, size_t numPtcls) -{ - if (walker_list_.empty()) - { - while (n) - { - walker_list_.push_back(std::make_unique(numPtcls)); - --n; - } - } - else - { - if (walker_list_.size() >= n) - { - int iw = walker_list_.size(); //copy from the back - for (int i = 0; i < n; ++i) - { - walker_list_.push_back(std::make_unique(*walker_list_[--iw])); - } - } - else - { - int nc = n / walker_list_.size(); - int nw0 = walker_list_.size(); - for (int iw = 0; iw < nw0; ++iw) - { - for (int ic = 0; ic < nc; ++ic) - walker_list_.push_back(std::make_unique(*walker_list_[iw])); - } - n -= nc * nw0; - while (n > 0) - { - walker_list_.push_back(std::make_unique(*walker_list_[--nw0])); - --n; - } - } - } -} - - -void WalkerConfigurations::resize(int numWalkers, size_t numPtcls) -{ - int dn = numWalkers - walker_list_.size(); - if (dn > 0) - createWalkers(dn, numPtcls); - if (dn < 0) - { - int nw = -dn; - if (nw < walker_list_.size()) - { - walker_list_.erase(walker_list_.begin(), walker_list_.begin() - dn); - } - } -} - -///returns the next valid iterator -WalkerConfigurations::iterator WalkerConfigurations::destroyWalkers(iterator first, iterator last) -{ - return walker_list_.erase(first, last); -} - -void WalkerConfigurations::createWalkers(iterator first, iterator last) -{ - destroyWalkers(walker_list_.begin(), walker_list_.end()); - while (first != last) - { - walker_list_.push_back(std::make_unique(**first)); - ++first; - } -} - -void WalkerConfigurations::destroyWalkers(int nw) -{ - if (nw > walker_list_.size()) - { - app_warning() << " Cannot remove walkers. Current Walkers = " << walker_list_.size() << std::endl; - return; - } - nw = walker_list_.size() - nw; - int iw = nw; - walker_list_.erase(walker_list_.begin() + nw, walker_list_.end()); -} - -void WalkerConfigurations::copyWalkers(iterator first, iterator last, iterator it) -{ - while (first != last) - { - (*it++)->makeCopy(**first++); - } -} - -/** Make Metropolis move to the walkers and save in a temporary array. - * @param it the iterator of the first walker to work on - * @param tauinv inverse of the time step - * - * R + D + X - */ -void WalkerConfigurations::reset() -{ - for (auto& walker : walker_list_) - { - walker->Weight = 1.0; - walker->Multiplicity = 1.0; - } -} - -void WalkerConfigurations::putConfigurations(Walker_t::RealType* target, QMCTraits::FullPrecRealType* weights) const -{ - for (const auto& walker : walker_list_) - { - std::copy(get_first_address(walker->R), get_last_address(walker->R), target); - target += get_last_address(walker->R) - get_first_address(walker->R); - *weights = walker->Weight; - ++weights; - } -} - -} // namespace qmcplusplus diff --git a/src/Particle/WalkerConfigurations.h b/src/Particle/WalkerConfigurations.h index ae3e96915b..da0fe17853 100644 --- a/src/Particle/WalkerConfigurations.h +++ b/src/Particle/WalkerConfigurations.h @@ -23,169 +23,10 @@ #ifndef QMCPLUSPLUS_WALKERCONFIGURATIONS_H #define QMCPLUSPLUS_WALKERCONFIGURATIONS_H #include "Configuration.h" -#include "Particle/Walker.h" -#include "Utilities/IteratorUtility.h" -#include "Particle/ParticleSetTraits.h" +#include "Particle/WalkerConfigurationsT.h" namespace qmcplusplus { -/** Monte Carlo Data of an ensemble - * - * The quantities are shared by all the nodes in a group - * - NumSamples number of samples - * - Weight total weight of a sample - * - Energy average energy of a sample - * - Variance variance - * - LivingFraction fraction of walkers alive each step. - */ -template -struct MCDataType -{ - T NumSamples; - T RNSamples; - T Weight; - T Energy; - T AlternateEnergy; - T Variance; - T R2Accepted; - T R2Proposed; - T LivingFraction; -}; - -/** A set of light weight walkers that are carried between driver sections and restart - */ -class WalkerConfigurations -{ -public: - /// walker type - using Walker_t = Walker, LatticeParticleTraits>; - using FullPrecRealType = QMCTraits::FullPrecRealType; - ///container type of Walkers - using walker_list__t = std::vector>; - /// FIX: a type alias of iterator for an object should not be for just one of many objects it holds. - using iterator = walker_list__t::iterator; - ///const_iterator of Walker container - using const_iterator = walker_list__t::const_iterator; - - MCDataType EnsembleProperty; - - WalkerConfigurations(); - ~WalkerConfigurations(); - WalkerConfigurations(const WalkerConfigurations&) = delete; - WalkerConfigurations& operator=(const WalkerConfigurations&) = delete; - WalkerConfigurations(WalkerConfigurations&&) = default; - WalkerConfigurations& operator=(WalkerConfigurations&&) = default; - - /** create numWalkers Walkers - * - * Append Walkers to walker_list_. - */ - void createWalkers(int numWalkers, size_t numPtcls); - /** create walkers - * @param first walker iterator - * @param last walker iterator - */ - void createWalkers(iterator first, iterator last); - /** copy walkers - * @param first input walker iterator - * @param last input walker iterator - * @param start first target iterator - * - * No memory allocation is allowed. - */ - void copyWalkers(iterator first, iterator last, iterator start); - - /** destroy Walkers from itstart to itend - *@param first starting iterator of the walkers - *@param last ending iterator of the walkers - */ - iterator destroyWalkers(iterator first, iterator last); - - /** destroy Walkers - *@param nw number of walkers to be destroyed - */ - void destroyWalkers(int nw); - - ///clean up the walker list and make a new list - void resize(int numWalkers, size_t numPtcls); - - ///return the number of active walkers - inline size_t getActiveWalkers() const { return walker_list_.size(); } - ///return the total number of active walkers among a MPI group - inline size_t getGlobalNumWalkers() const { return walker_offsets_.empty() ? 0 : walker_offsets_.back(); } - ///return the total number of active walkers among a MPI group - - inline void setWalkerOffsets(const std::vector& o) { walker_offsets_ = o; } - inline const std::vector& getWalkerOffsets() const { return walker_offsets_; } - - /// return the first iterator - inline iterator begin() { return walker_list_.begin(); } - /// return the last iterator, [begin(), end()) - inline iterator end() { return walker_list_.end(); } - - /// return the first const_iterator - inline const_iterator begin() const { return walker_list_.begin(); } - - /// return the last const_iterator [begin(), end()) - inline const_iterator end() const { return walker_list_.end(); } - /**@}*/ - - /** clear the walker_list_ without destroying them - * - * Provide std::vector::clear interface - */ - inline void clear() { walker_list_.clear(); } - - /** insert elements - * @param it locator where the inserting begins - * @param first starting iterator - * @param last ending iterator - * - * Provide std::vector::insert interface - */ - template - inline void insert(iterator it, INPUT_ITER first, INPUT_ITER last) - { - walker_list_.insert(it, first, last); - } - - /** add Walker_t* at the end - * @param awalker pointer to a walker - * - * Provide std::vector::push_back interface - */ - inline void push_back(std::unique_ptr awalker) { walker_list_.push_back(std::move(awalker)); } - - /** delete the last Walker_t* - * - * Provide std::vector::pop_back interface - */ - inline void pop_back() { walker_list_.pop_back(); } - - inline Walker_t* operator[](int i) { return walker_list_[i].get(); } - - inline const Walker_t* operator[](int i) const { return walker_list_[i].get(); } - - /** reset the Walkers - */ - void reset(); - - ///save the particle positions of all the walkers into target - void putConfigurations(Walker_t::RealType* target, QMCTraits::FullPrecRealType* weights) const; - -protected: - ///a collection of walkers - walker_list__t walker_list_; - -private: - /** starting index of the walkers in a processor group - * - * walker_offsets_[0]=0 and walker_offsets_[walker_offsets_.size()-1]=total number of walkers in a group - * walker_offsets_[processorid+1]-walker_offsets_[processorid] is equal to the number of walkers on a processor, - * i.e., W.getActiveWalkers(). - * walker_offsets_ is added to handle parallel I/O with hdf5 - */ - std::vector walker_offsets_; -}; +using WalkerConfigurations = WalkerConfigurationsT; } // namespace qmcplusplus #endif diff --git a/src/Particle/WalkerConfigurationsT.cpp b/src/Particle/WalkerConfigurationsT.cpp new file mode 100644 index 0000000000..bf5642e2c3 --- /dev/null +++ b/src/Particle/WalkerConfigurationsT.cpp @@ -0,0 +1,170 @@ +////////////////////////////////////////////////////////////////////////////////////// +// 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: Jordan E. Vincent, University of Illinois at +// Urbana-Champaign +// Bryan Clark, bclark@Princeton.edu, Princeton University +// Ken Esler, kpesler@gmail.com, University of Illinois at +// Urbana-Champaign Jeremy McMinnis, jmcminis@gmail.com, +// University of Illinois at Urbana-Champaign Jeongnim Kim, +// jeongnim.kim@gmail.com, University of Illinois at +// Urbana-Champaign Cynthia Gu, zg1@ornl.gov, Oak Ridge +// National Laboratory Ye Luo, yeluo@anl.gov, Argonne +// National Laboratory Mark A. Berrill, berrillma@ornl.gov, +// Oak Ridge National Laboratory +// +// File created by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory +////////////////////////////////////////////////////////////////////////////////////// + +#include "WalkerConfigurationsT.h" + +#include "Utilities/IteratorUtility.h" +#include "Platforms/Host/OutputManager.h" + +#include + +namespace qmcplusplus +{ +template +WalkerConfigurationsT::WalkerConfigurationsT() = default; + +/// default destructor +template +WalkerConfigurationsT::~WalkerConfigurationsT() +{ + destroyWalkers(walker_list_.begin(), walker_list_.end()); +} + +template +void +WalkerConfigurationsT::createWalkers(int n, size_t numPtcls) +{ + if (walker_list_.empty()) { + while (n) { + walker_list_.push_back(std::make_unique(numPtcls)); + --n; + } + } + else { + if (walker_list_.size() >= n) { + int iw = walker_list_.size(); // copy from the back + for (int i = 0; i < n; ++i) { + walker_list_.push_back( + std::make_unique(*walker_list_[--iw])); + } + } + else { + int nc = n / walker_list_.size(); + int nw0 = walker_list_.size(); + for (int iw = 0; iw < nw0; ++iw) { + for (int ic = 0; ic < nc; ++ic) + walker_list_.push_back( + std::make_unique(*walker_list_[iw])); + } + n -= nc * nw0; + while (n > 0) { + walker_list_.push_back( + std::make_unique(*walker_list_[--nw0])); + --n; + } + } + } +} + +template +void +WalkerConfigurationsT::resize(int numWalkers, size_t numPtcls) +{ + int dn = numWalkers - walker_list_.size(); + if (dn > 0) + createWalkers(dn, numPtcls); + if (dn < 0) { + int nw = -dn; + if (nw < walker_list_.size()) { + walker_list_.erase(walker_list_.begin(), walker_list_.begin() - dn); + } + } +} + +/// returns the next valid iterator +template +typename WalkerConfigurationsT::iterator +WalkerConfigurationsT::destroyWalkers(iterator first, iterator last) +{ + return walker_list_.erase(first, last); +} + +template +void +WalkerConfigurationsT::createWalkers(iterator first, iterator last) +{ + destroyWalkers(walker_list_.begin(), walker_list_.end()); + while (first != last) { + walker_list_.push_back(std::make_unique(**first)); + ++first; + } +} + +template +void +WalkerConfigurationsT::destroyWalkers(int nw) +{ + if (nw > walker_list_.size()) { + app_warning() << " Cannot remove walkers. Current Walkers = " + << walker_list_.size() << std::endl; + return; + } + nw = walker_list_.size() - nw; + int iw = nw; + walker_list_.erase(walker_list_.begin() + nw, walker_list_.end()); +} + +template +void +WalkerConfigurationsT::copyWalkers( + iterator first, iterator last, iterator it) +{ + while (first != last) { + (*it++)->makeCopy(**first++); + } +} + +/** Make Metropolis move to the walkers and save in a temporary array. + * @param it the iterator of the first walker to work on + * @param tauinv inverse of the time step + * + * R + D + X + */ +template +void +WalkerConfigurationsT::reset() +{ + for (auto& walker : walker_list_) { + walker->Weight = 1.0; + walker->Multiplicity = 1.0; + } +} + +template +void +WalkerConfigurationsT::putConfigurations( + RealType* target, FullPrecRealType* weights) const +{ + for (const auto& walker : walker_list_) { + std::copy( + get_first_address(walker->R), get_last_address(walker->R), target); + target += get_last_address(walker->R) - get_first_address(walker->R); + *weights = walker->Weight; + ++weights; + } +} + +template class WalkerConfigurationsT; +template class WalkerConfigurationsT; +template class WalkerConfigurationsT>; +template class WalkerConfigurationsT>; + +} // namespace qmcplusplus diff --git a/src/Particle/WalkerConfigurationsT.h b/src/Particle/WalkerConfigurationsT.h new file mode 100644 index 0000000000..7b9cae36d3 --- /dev/null +++ b/src/Particle/WalkerConfigurationsT.h @@ -0,0 +1,258 @@ +////////////////////////////////////////////////////////////////////////////////////// +// 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: Jordan E. Vincent, University of Illinois at +// Urbana-Champaign +// Ken Esler, kpesler@gmail.com, University of Illinois at +// Urbana-Champaign Jeremy McMinnis, jmcminis@gmail.com, +// University of Illinois at Urbana-Champaign Jeongnim Kim, +// jeongnim.kim@gmail.com, University of Illinois at +// Urbana-Champaign Cynthia Gu, zg1@ornl.gov, Oak Ridge +// National Laboratory Raymond Clay III, +// j.k.rofling@gmail.com, Lawrence Livermore National +// Laboratory Ye Luo, yeluo@anl.gov, Argonne National +// Laboratory Mark A. Berrill, berrillma@ornl.gov, Oak Ridge +// National Laboratory +// +// File created by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory +////////////////////////////////////////////////////////////////////////////////////// + +#ifndef QMCPLUSPLUS_WALKERCONFIGURATIONST_H +#define QMCPLUSPLUS_WALKERCONFIGURATIONST_H +#include "Particle/ParticleSetTraits.h" +#include "Particle/Walker.h" +#include "Utilities/IteratorUtility.h" + +namespace qmcplusplus +{ +/** Monte Carlo Data of an ensemble + * + * The quantities are shared by all the nodes in a group + * - NumSamples number of samples + * - Weight total weight of a sample + * - Energy average energy of a sample + * - Variance variance + * - LivingFraction fraction of walkers alive each step. + */ +template +struct MCDataType +{ + T NumSamples; + T RNSamples; + T Weight; + T Energy; + T AlternateEnergy; + T Variance; + T R2Accepted; + T R2Proposed; + T LivingFraction; +}; + +/** A set of light weight walkers that are carried between driver sections and + * restart + */ +template +class WalkerConfigurationsT +{ +public: + /// walker type + using Walker_t = Walker, LatticeParticleTraits>; + using RealType = typename Walker_t::RealType; + using FullPrecRealType = typename ParticleSetTraits::FullPrecRealType; + /// container type of Walkers + using walker_list__t = std::vector>; + /// FIX: a type alias of iterator for an object should not be for just one + /// of many objects it holds. + using iterator = typename walker_list__t::iterator; + /// const_iterator of Walker container + using const_iterator = typename walker_list__t::const_iterator; + + MCDataType EnsembleProperty; + + WalkerConfigurationsT(); + ~WalkerConfigurationsT(); + WalkerConfigurationsT(const WalkerConfigurationsT&) = delete; + WalkerConfigurationsT& + operator=(const WalkerConfigurationsT&) = delete; + WalkerConfigurationsT(WalkerConfigurationsT&&) = default; + WalkerConfigurationsT& + operator=(WalkerConfigurationsT&&) = default; + + /** create numWalkers Walkers + * + * Append Walkers to walker_list_. + */ + void + createWalkers(int numWalkers, size_t numPtcls); + /** create walkers + * @param first walker iterator + * @param last walker iterator + */ + void + createWalkers(iterator first, iterator last); + /** copy walkers + * @param first input walker iterator + * @param last input walker iterator + * @param start first target iterator + * + * No memory allocation is allowed. + */ + void + copyWalkers(iterator first, iterator last, iterator start); + + /** destroy Walkers from itstart to itend + *@param first starting iterator of the walkers + *@param last ending iterator of the walkers + */ + iterator + destroyWalkers(iterator first, iterator last); + + /** destroy Walkers + *@param nw number of walkers to be destroyed + */ + void + destroyWalkers(int nw); + + /// clean up the walker list and make a new list + void + resize(int numWalkers, size_t numPtcls); + + /// return the number of active walkers + inline size_t + getActiveWalkers() const + { + return walker_list_.size(); + } + /// return the total number of active walkers among a MPI group + inline size_t + getGlobalNumWalkers() const + { + return walker_offsets_.empty() ? 0 : walker_offsets_.back(); + } + /// return the total number of active walkers among a MPI group + + inline void + setWalkerOffsets(const std::vector& o) + { + walker_offsets_ = o; + } + inline const std::vector& + getWalkerOffsets() const + { + return walker_offsets_; + } + + /// return the first iterator + inline iterator + begin() + { + return walker_list_.begin(); + } + /// return the last iterator, [begin(), end()) + inline iterator + end() + { + return walker_list_.end(); + } + + /// return the first const_iterator + inline const_iterator + begin() const + { + return walker_list_.begin(); + } + + /// return the last const_iterator [begin(), end()) + inline const_iterator + end() const + { + return walker_list_.end(); + } + /**@}*/ + + /** clear the walker_list_ without destroying them + * + * Provide std::vector::clear interface + */ + inline void + clear() + { + walker_list_.clear(); + } + + /** insert elements + * @param it locator where the inserting begins + * @param first starting iterator + * @param last ending iterator + * + * Provide std::vector::insert interface + */ + template + inline void + insert(iterator it, INPUT_ITER first, INPUT_ITER last) + { + walker_list_.insert(it, first, last); + } + + /** add Walker_t* at the end + * @param awalker pointer to a walker + * + * Provide std::vector::push_back interface + */ + inline void + push_back(std::unique_ptr awalker) + { + walker_list_.push_back(std::move(awalker)); + } + + /** delete the last Walker_t* + * + * Provide std::vector::pop_back interface + */ + inline void + pop_back() + { + walker_list_.pop_back(); + } + + inline Walker_t* + operator[](int i) + { + return walker_list_[i].get(); + } + + inline const Walker_t* + operator[](int i) const + { + return walker_list_[i].get(); + } + + /** reset the Walkers + */ + void + reset(); + + /// save the particle positions of all the walkers into target + void + putConfigurations(RealType* target, FullPrecRealType* weights) const; + +protected: + /// a collection of walkers + walker_list__t walker_list_; + +private: + /** starting index of the walkers in a processor group + * + * walker_offsets_[0]=0 and walker_offsets_[walker_offsets_.size()-1]=total + * number of walkers in a group + * walker_offsets_[processorid+1]-walker_offsets_[processorid] is equal to + * the number of walkers on a processor, i.e., W.getActiveWalkers(). + * walker_offsets_ is added to handle parallel I/O with hdf5 + */ + std::vector walker_offsets_; +}; +} // namespace qmcplusplus +#endif diff --git a/src/QMCDrivers/MCPopulation.h b/src/QMCDrivers/MCPopulation.h index fd38c242d0..b66a86bca0 100644 --- a/src/QMCDrivers/MCPopulation.h +++ b/src/QMCDrivers/MCPopulation.h @@ -33,7 +33,7 @@ class QMCHamiltonian; class MCPopulation { public: - using MCPWalker = Walker; + using MCPWalker = MCWalkerConfiguration::Walker_t; using WFBuffer = MCPWalker::WFBuffer_t; using RealType = QMCTraits::RealType; using Properties = MCPWalker::PropertyContainer_t; diff --git a/src/QMCDrivers/WFOpt/QMCCostFunctionBatched.h b/src/QMCDrivers/WFOpt/QMCCostFunctionBatched.h index b9440a6bfb..305cb180fc 100644 --- a/src/QMCDrivers/WFOpt/QMCCostFunctionBatched.h +++ b/src/QMCDrivers/WFOpt/QMCCostFunctionBatched.h @@ -18,6 +18,7 @@ #include "QMCDrivers/WFOpt/QMCCostFunctionBase.h" #include "QMCDrivers/CloneManager.h" #include "QMCWaveFunctions/OrbitalSetTraits.h" +#include "SampleStack.h" namespace qmcplusplus { diff --git a/src/QMCDrivers/WalkerElementsRef.h b/src/QMCDrivers/WalkerElementsRef.h index 1a11de623c..d5d35a6bcc 100644 --- a/src/QMCDrivers/WalkerElementsRef.h +++ b/src/QMCDrivers/WalkerElementsRef.h @@ -1,6 +1,6 @@ ////////////////////////////////////////////////////////////////////////////////////// -// This file is distributed under the University of Illinois/NCSA Open Source License. -// See LICENSE file in top directory for details. +// 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. // @@ -13,6 +13,7 @@ #define QMCPLUSPLUS_WALKERELEMENTSREF_H #include "Configuration.h" +#include "Particle/ParticleSetTraits.h" #include "Particle/Walker.h" namespace qmcplusplus @@ -22,28 +23,35 @@ class TrialWaveFunction; /** type for returning the walker and its elements from MCPopulation * - * have no expectations for the validity of the references in this structure past - * the context it was returned in. It should not be returned by a call in a + * have no expectations for the validity of the references in this structure + * past the context it was returned in. It should not be returned by a call in a * crowd or threaded context. - * + * * @ye-luo's "fat" walker * - * We need this if we want to "copyFrom" the whole fat walker when it comes off the line - * i.e. mpi. Insuring the "fat" walker is valid at the earliest possible point seems - * less likely to end in tears then just calling copyFrom random other places (hopefully) - * in time, in order to not access an invalid walker element. + * We need this if we want to "copyFrom" the whole fat walker when it comes off + * the line i.e. mpi. Insuring the "fat" walker is valid at the earliest + * possible point seems less likely to end in tears then just calling copyFrom + * random other places (hopefully) in time, in order to not access an invalid + * walker element. */ struct WalkerElementsRef { - /** to allow use of emplace back - */ - WalkerElementsRef(Walker& walker_in, ParticleSet& pset_in, TrialWaveFunction& twf_in) : walker(walker_in), pset(pset_in), twf(twf_in) {} -; - Walker& walker; - ParticleSet& pset; - TrialWaveFunction& twf; + using WalkerType = Walker, + LatticeParticleTraits>; + /** to allow use of emplace back + */ + WalkerElementsRef(WalkerType& walker_in, ParticleSet& pset_in, + TrialWaveFunction& twf_in) : + walker(walker_in), + pset(pset_in), + twf(twf_in){}; + + WalkerType& walker; + ParticleSet& pset; + TrialWaveFunction& twf; }; -} +} // namespace qmcplusplus #endif diff --git a/src/QMCDrivers/tests/WalkerConsumer.h b/src/QMCDrivers/tests/WalkerConsumer.h index 689a2be280..91bc9675bf 100644 --- a/src/QMCDrivers/tests/WalkerConsumer.h +++ b/src/QMCDrivers/tests/WalkerConsumer.h @@ -1,6 +1,6 @@ ////////////////////////////////////////////////////////////////////////////////////// -// This file is distributed under the University of Illinois/NCSA Open Source License. -// See LICENSE file in top directory for details. +// This file is distributed under the University of Illinois/NCSA Open Source +// License. See LICENSE file in top directory for details. // // Copyright (c) 2019 QMCPACK developers. // @@ -12,10 +12,12 @@ #ifndef QMCPLUSPLUS_WALKERCONSUMER_H #define QMCPLUSPLUS_WALKERCONSUMER_H -#include - +#include "Configuration.h" +#include "Particle/ParticleSetTraits.h" #include "Particle/Walker.h" +#include + namespace qmcplusplus { class ResourceCollection; @@ -29,32 +31,38 @@ namespace testing class WalkerConsumer { public: - std::vector>> walkers; - std::vector> walker_elecs_; - std::vector> walker_twfs_; - std::vector> walker_hamiltonians_; - - void initializeResources(const ResourceCollection& twf_resource) {} - - void addWalker(Walker& walker, - ParticleSet& elecs, - TrialWaveFunction& twf, - QMCHamiltonian& hamiltonian) - { - walkers.push_back(walker); - walker_elecs_.push_back(elecs); - walker_twfs_.push_back(twf); - walker_hamiltonians_.push_back(hamiltonian); - } - - void clearWalkers() - { - // We're clearing the refs to the objects not the referred to objects. - walkers.clear(); - walker_elecs_.clear(); - walker_twfs_.clear(); - walker_hamiltonians_.clear(); - } + using WalkerType = Walker, + LatticeParticleTraits>; + + std::vector> walkers; + std::vector> walker_elecs_; + std::vector> walker_twfs_; + std::vector> walker_hamiltonians_; + + void + initializeResources(const ResourceCollection& twf_resource) + { + } + + void + addWalker(WalkerType& walker, ParticleSet& elecs, TrialWaveFunction& twf, + QMCHamiltonian& hamiltonian) + { + walkers.push_back(walker); + walker_elecs_.push_back(elecs); + walker_twfs_.push_back(twf); + walker_hamiltonians_.push_back(hamiltonian); + } + + void + clearWalkers() + { + // We're clearing the refs to the objects not the referred to objects. + walkers.clear(); + walker_elecs_.clear(); + walker_twfs_.clear(); + walker_hamiltonians_.clear(); + } }; } // namespace testing diff --git a/src/QMCDrivers/tests/test_Crowd.cpp b/src/QMCDrivers/tests/test_Crowd.cpp index 96d5d98a1c..ea854406ce 100644 --- a/src/QMCDrivers/tests/test_Crowd.cpp +++ b/src/QMCDrivers/tests/test_Crowd.cpp @@ -29,7 +29,7 @@ namespace testing class CrowdWithWalkers { public: - using MCPWalker = Walker; + using MCPWalker = Crowd::MCPWalker; EstimatorManagerNew em; UPtr crowd_ptr; diff --git a/src/QMCDrivers/tests/test_DMCBatched.cpp b/src/QMCDrivers/tests/test_DMCBatched.cpp index 1a069efb3b..e846141452 100644 --- a/src/QMCDrivers/tests/test_DMCBatched.cpp +++ b/src/QMCDrivers/tests/test_DMCBatched.cpp @@ -20,6 +20,7 @@ #include "Concurrency/Info.hpp" #include "Concurrency/UtilityFunctions.hpp" #include "Platforms/Host/OutputManager.h" +#include "SampleStack.h" namespace qmcplusplus { diff --git a/src/QMCDrivers/tests/test_SFNBranch.cpp b/src/QMCDrivers/tests/test_SFNBranch.cpp index d2ccdc1d5d..5291e7d36f 100644 --- a/src/QMCDrivers/tests/test_SFNBranch.cpp +++ b/src/QMCDrivers/tests/test_SFNBranch.cpp @@ -26,7 +26,7 @@ namespace qmcplusplus { -using MCPWalker = Walker; +using MCPWalker = MCPopulation::MCPWalker; using RealType = double; namespace testing diff --git a/src/QMCWaveFunctions/BsplineFactory/SplineSetReaderT.h b/src/QMCWaveFunctions/BsplineFactory/SplineSetReaderT.h index 95ce5d5590..4d01235b06 100644 --- a/src/QMCWaveFunctions/BsplineFactory/SplineSetReaderT.h +++ b/src/QMCWaveFunctions/BsplineFactory/SplineSetReaderT.h @@ -18,6 +18,7 @@ #ifndef QMCPLUSPLUS_SPLINESET_READERT_H #define QMCPLUSPLUS_SPLINESET_READERT_H +#include "EinsplineSetBuilderT.h" #include "BsplineFactory/BsplineReaderBaseT.h" #include "Utilities/FairDivide.h" #include "Utilities/ProgressReportEngine.h" diff --git a/src/QMCWaveFunctions/BsplineFactory/createBsplineReaderT.cpp b/src/QMCWaveFunctions/BsplineFactory/createBsplineReaderT.cpp index f7b5e17c77..1d3b623267 100644 --- a/src/QMCWaveFunctions/BsplineFactory/createBsplineReaderT.cpp +++ b/src/QMCWaveFunctions/BsplineFactory/createBsplineReaderT.cpp @@ -216,6 +216,7 @@ createBsplineComplexDoubleT( return CreateComplexHelper::createDouble(e, hybrid_rep, useGPU); } +#ifdef QMC_COMPLEX template std::unique_ptr>> createBsplineComplexDoubleT>( EinsplineSetBuilderT>* e, bool hybrid_rep, @@ -226,6 +227,7 @@ createBsplineComplexDoubleT>( EinsplineSetBuilderT>* e, bool hybrid_rep, const std::string& useGPU); +#else template std::unique_ptr> createBsplineComplexDoubleT( EinsplineSetBuilderT* e, bool hybrid_rep, const std::string& useGPU); @@ -233,6 +235,7 @@ createBsplineComplexDoubleT( template std::unique_ptr> createBsplineComplexDoubleT(EinsplineSetBuilderT* e, bool hybrid_rep, const std::string& useGPU); +#endif template std::unique_ptr> @@ -242,6 +245,7 @@ createBsplineComplexSingleT( return CreateComplexHelper::createSingle(e, hybrid_rep, useGPU); } +#ifdef QMC_COMPLEX template std::unique_ptr>> createBsplineComplexSingleT>( EinsplineSetBuilderT>* e, bool hybrid_rep, @@ -252,6 +256,7 @@ createBsplineComplexSingleT>( EinsplineSetBuilderT>* e, bool hybrid_rep, const std::string& useGPU); +#else template std::unique_ptr> createBsplineComplexSingleT( EinsplineSetBuilderT* e, bool hybrid_rep, const std::string& useGPU); @@ -259,6 +264,7 @@ createBsplineComplexSingleT( template std::unique_ptr> createBsplineComplexSingleT(EinsplineSetBuilderT* e, bool hybrid_rep, const std::string& useGPU); +#endif template std::unique_ptr> @@ -286,6 +292,7 @@ createBsplineRealDoubleT( return aReader; } +#ifndef QMC_COMPLEX template std::unique_ptr> createBsplineRealDoubleT( EinsplineSetBuilderT* e, bool hybrid_rep, const std::string& useGPU); @@ -293,6 +300,7 @@ createBsplineRealDoubleT( template std::unique_ptr> createBsplineRealDoubleT(EinsplineSetBuilderT* e, bool hybrid_rep, const std::string& useGPU); +#endif template std::unique_ptr> @@ -320,6 +328,7 @@ createBsplineRealSingleT( return aReader; } +#ifndef QMC_COMPLEX template std::unique_ptr> createBsplineRealSingleT( EinsplineSetBuilderT* e, bool hybrid_rep, const std::string& useGPU); @@ -327,5 +336,6 @@ createBsplineRealSingleT( template std::unique_ptr> createBsplineRealSingleT(EinsplineSetBuilderT* e, bool hybrid_rep, const std::string& useGPU); +#endif } // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/CompositeSPOSetT.cpp b/src/QMCWaveFunctions/CompositeSPOSetT.cpp index 31a3f71399..d39aaaa82e 100644 --- a/src/QMCWaveFunctions/CompositeSPOSetT.cpp +++ b/src/QMCWaveFunctions/CompositeSPOSetT.cpp @@ -217,9 +217,12 @@ CompositeSPOSetBuilderT::createSPOSet(xmlNodePtr cur, SPOSetInputInfo& input) return createSPOSetFromXML(cur); } +#ifndef QMC_COMPLEX template class CompositeSPOSetBuilderT; template class CompositeSPOSetBuilderT; +#else template class CompositeSPOSetBuilderT>; template class CompositeSPOSetBuilderT>; +#endif } // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/EinsplineSetBuilderT.cpp b/src/QMCWaveFunctions/EinsplineSetBuilderT.cpp index 46157f9b28..7d14613afa 100644 --- a/src/QMCWaveFunctions/EinsplineSetBuilderT.cpp +++ b/src/QMCWaveFunctions/EinsplineSetBuilderT.cpp @@ -1598,6 +1598,7 @@ EinsplineSetBuilderT::createBsplineReader( } } +#ifdef QMC_COMPLEX template <> void EinsplineSetBuilderT>::createBsplineReader( @@ -1627,6 +1628,7 @@ EinsplineSetBuilderT>::createBsplineReader( createBsplineComplexDoubleT(this, hybridRep, useGPU); } } +#endif template std::unique_ptr> @@ -1809,9 +1811,12 @@ EinsplineSetBuilderT::ReadGvectors_ESHDF() return hasPsig; } +#ifndef QMC_COMPLEX template class EinsplineSetBuilderT; template class EinsplineSetBuilderT; +#else template class EinsplineSetBuilderT>; template class EinsplineSetBuilderT>; +#endif } // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/HarmonicOscillator/SHOSetBuilderT.cpp b/src/QMCWaveFunctions/HarmonicOscillator/SHOSetBuilderT.cpp index 77ae1eda5a..7c309d5b87 100644 --- a/src/QMCWaveFunctions/HarmonicOscillator/SHOSetBuilderT.cpp +++ b/src/QMCWaveFunctions/HarmonicOscillator/SHOSetBuilderT.cpp @@ -208,9 +208,12 @@ SHOSetBuilderT::report(const std::string& pad) app_log().flush(); } +#ifndef QMC_COMPLEX template class SHOSetBuilderT; template class SHOSetBuilderT; +#else template class SHOSetBuilderT>; template class SHOSetBuilderT>; +#endif } // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/HarmonicOscillator/SHOSetT.cpp b/src/QMCWaveFunctions/HarmonicOscillator/SHOSetT.cpp index b4e55a258d..1286b07393 100644 --- a/src/QMCWaveFunctions/HarmonicOscillator/SHOSetT.cpp +++ b/src/QMCWaveFunctions/HarmonicOscillator/SHOSetT.cpp @@ -553,9 +553,12 @@ SHOSetT::evaluateGradSource(const ParticleSetT& P, int first, int last, } // Class concrete types from ValueType +#ifndef QMC_COMPLEX template class SHOSetT; template class SHOSetT; +#else template class SHOSetT>; template class SHOSetT>; +#endif } // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/LCAO/LCAOrbitalBuilderT.cpp b/src/QMCWaveFunctions/LCAO/LCAOrbitalBuilderT.cpp index 39ea3953ee..635f3e61d9 100644 --- a/src/QMCWaveFunctions/LCAO/LCAOrbitalBuilderT.cpp +++ b/src/QMCWaveFunctions/LCAO/LCAOrbitalBuilderT.cpp @@ -1180,8 +1180,11 @@ LCAOrbitalBuilderT::EvalPeriodicImagePhaseFactors(PosType SuperTwist, } } +#ifndef QMC_COMPLEX template class LCAOrbitalBuilderT; template class LCAOrbitalBuilderT; +#else template class LCAOrbitalBuilderT>; template class LCAOrbitalBuilderT>; +#endif } // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/LCAO/LCAOrbitalSetT.cpp b/src/QMCWaveFunctions/LCAO/LCAOrbitalSetT.cpp index 6abd2d8b22..698c35dc59 100644 --- a/src/QMCWaveFunctions/LCAO/LCAOrbitalSetT.cpp +++ b/src/QMCWaveFunctions/LCAO/LCAOrbitalSetT.cpp @@ -954,9 +954,12 @@ LCAOrbitalSetT::applyRotation( } // Class concrete types from ValueType +#ifndef QMC_COMPLEX template class LCAOrbitalSetT; template class LCAOrbitalSetT; +#else template class LCAOrbitalSetT>; template class LCAOrbitalSetT>; +#endif } // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/LCAO/SoaLocalizedBasisSetT.cpp b/src/QMCWaveFunctions/LCAO/SoaLocalizedBasisSetT.cpp index 7b62735768..a17a845ac7 100644 --- a/src/QMCWaveFunctions/LCAO/SoaLocalizedBasisSetT.cpp +++ b/src/QMCWaveFunctions/LCAO/SoaLocalizedBasisSetT.cpp @@ -365,105 +365,124 @@ SoaLocalizedBasisSetT::add(int icenter, std::unique_ptr aos) // TODO: this should be redone with template template parameters +#ifndef QMC_COMPLEX template class SoaLocalizedBasisSetT< SoaAtomicBasisSetT, SoaCartesianTensor, double>, double>; -template class SoaLocalizedBasisSetT< - SoaAtomicBasisSetT, SoaCartesianTensor, - std::complex>, - std::complex>; template class SoaLocalizedBasisSetT< SoaAtomicBasisSetT, SoaCartesianTensor, float>, float>; +#else +template class SoaLocalizedBasisSetT< + SoaAtomicBasisSetT, SoaCartesianTensor, + std::complex>, + std::complex>; template class SoaLocalizedBasisSetT< SoaAtomicBasisSetT, SoaCartesianTensor, std::complex>, std::complex>; +#endif +#ifndef QMC_COMPLEX template class SoaLocalizedBasisSetT< SoaAtomicBasisSetT, SoaSphericalTensor, double>, double>; -template class SoaLocalizedBasisSetT< - SoaAtomicBasisSetT, SoaSphericalTensor, - std::complex>, - std::complex>; template class SoaLocalizedBasisSetT< SoaAtomicBasisSetT, SoaSphericalTensor, float>, float>; +#else +template class SoaLocalizedBasisSetT< + SoaAtomicBasisSetT, SoaSphericalTensor, + std::complex>, + std::complex>; template class SoaLocalizedBasisSetT< SoaAtomicBasisSetT, SoaSphericalTensor, std::complex>, std::complex>; +#endif +#ifndef QMC_COMPLEX template class SoaLocalizedBasisSetT< SoaAtomicBasisSetT>, SoaCartesianTensor, double>, double>; -template class SoaLocalizedBasisSetT< - SoaAtomicBasisSetT>, - SoaCartesianTensor, std::complex>, - std::complex>; template class SoaLocalizedBasisSetT< SoaAtomicBasisSetT>, SoaCartesianTensor, float>, float>; +#else +template class SoaLocalizedBasisSetT< + SoaAtomicBasisSetT>, + SoaCartesianTensor, std::complex>, + std::complex>; template class SoaLocalizedBasisSetT< SoaAtomicBasisSetT>, SoaCartesianTensor, std::complex>, std::complex>; +#endif +#ifndef QMC_COMPLEX template class SoaLocalizedBasisSetT< SoaAtomicBasisSetT>, SoaSphericalTensor, double>, double>; -template class SoaLocalizedBasisSetT< - SoaAtomicBasisSetT>, - SoaSphericalTensor, std::complex>, - std::complex>; template class SoaLocalizedBasisSetT< SoaAtomicBasisSetT>, SoaSphericalTensor, float>, float>; +#else +template class SoaLocalizedBasisSetT< + SoaAtomicBasisSetT>, + SoaSphericalTensor, std::complex>, + std::complex>; template class SoaLocalizedBasisSetT< SoaAtomicBasisSetT>, SoaSphericalTensor, std::complex>, std::complex>; +#endif +#ifndef QMC_COMPLEX template class SoaLocalizedBasisSetT< SoaAtomicBasisSetT>, SoaCartesianTensor, double>, double>; -template class SoaLocalizedBasisSetT< - SoaAtomicBasisSetT>, - SoaCartesianTensor, std::complex>, - std::complex>; template class SoaLocalizedBasisSetT< SoaAtomicBasisSetT>, SoaCartesianTensor, float>, float>; +#else +template class SoaLocalizedBasisSetT< + SoaAtomicBasisSetT>, + SoaCartesianTensor, std::complex>, + std::complex>; template class SoaLocalizedBasisSetT< SoaAtomicBasisSetT>, SoaCartesianTensor, std::complex>, std::complex>; +#endif +#ifndef QMC_COMPLEX template class SoaLocalizedBasisSetT< SoaAtomicBasisSetT>, SoaSphericalTensor, double>, double>; -template class SoaLocalizedBasisSetT< - SoaAtomicBasisSetT>, - SoaSphericalTensor, std::complex>, - std::complex>; template class SoaLocalizedBasisSetT< SoaAtomicBasisSetT>, SoaSphericalTensor, float>, float>; +#else +template class SoaLocalizedBasisSetT< + SoaAtomicBasisSetT>, + SoaSphericalTensor, std::complex>, + std::complex>; template class SoaLocalizedBasisSetT< SoaAtomicBasisSetT>, SoaSphericalTensor, std::complex>, std::complex>; +#endif + } // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/SPOSetBuilderFactoryT.cpp b/src/QMCWaveFunctions/SPOSetBuilderFactoryT.cpp index 12148277a0..d4cee7fba7 100644 --- a/src/QMCWaveFunctions/SPOSetBuilderFactoryT.cpp +++ b/src/QMCWaveFunctions/SPOSetBuilderFactoryT.cpp @@ -276,8 +276,11 @@ SPOSetBuilderFactoryT::addSPOSet(std::unique_ptr> spo) template std::string SPOSetBuilderFactoryT::basisset_tag = "basisset"; +#ifdef QMC_COMPLEX template class SPOSetBuilderFactoryT>; template class SPOSetBuilderFactoryT>; +#else template class SPOSetBuilderFactoryT; template class SPOSetBuilderFactoryT; +#endif } // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/VariableSetT.cpp b/src/QMCWaveFunctions/VariableSetT.cpp index 064ac26a13..2c49401066 100644 --- a/src/QMCWaveFunctions/VariableSetT.cpp +++ b/src/QMCWaveFunctions/VariableSetT.cpp @@ -338,9 +338,9 @@ VariableSetT::readFromHDF( hin.pop(); } -template struct VariableSetT; -template struct VariableSetT; -template struct VariableSetT>; -template struct VariableSetT>; +template class VariableSetT; +template class VariableSetT; +template class VariableSetT>; +template class VariableSetT>; } // namespace optimize diff --git a/src/QMCWaveFunctions/tests/CMakeLists.txt b/src/QMCWaveFunctions/tests/CMakeLists.txt index df01d4ba35..01db0f9123 100644 --- a/src/QMCWaveFunctions/tests/CMakeLists.txt +++ b/src/QMCWaveFunctions/tests/CMakeLists.txt @@ -145,13 +145,6 @@ add_library(sposets_for_testing FakeSPOT.cpp FakeSPO.cpp ConstantSPOSetT.cpp) target_include_directories(sposets_for_testing PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}") target_link_libraries(sposets_for_testing PUBLIC qmcwfs) -# @TODO: Remove when rotations work for complex stuff -if(NOT QMC_COMPLEX) - if(NOT ENABLE_CUDA) - set(SPOSET_SRC test_RotatedSPOs.cpp ${SPOSET_SRC}) - endif() -endif() - if(ENABLE_CUDA) set(DETERMINANT_SRC ${DETERMINANT_SRC} test_DiracMatrixComputeCUDA.cpp test_cuBLAS_LU.cpp) if(NOT QMC_CUDA2HIP) diff --git a/src/QMCWaveFunctions/tests/test_RotatedSPOs.cpp b/src/QMCWaveFunctions/tests/test_RotatedSPOs.cpp deleted file mode 100644 index da3f10b376..0000000000 --- a/src/QMCWaveFunctions/tests/test_RotatedSPOs.cpp +++ /dev/null @@ -1,869 +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) 2022 QMCPACK developers. -// -// File developed by: Joshua Townsend, jptowns@sandia.gov, Sandia National Laboratories -// -// File created by: Joshua Townsend, jptowns@sandia.gov, Sandia National Laboratories -////////////////////////////////////////////////////////////////////////////////////// - -#include "catch.hpp" - -#include "type_traits/template_types.hpp" -#include "type_traits/ConvertToReal.h" -#include "OhmmsData/Libxml2Doc.h" -#include "OhmmsPETE/OhmmsMatrix.h" -#include "Particle/ParticleSet.h" -#include "Particle/ParticleSetPool.h" -#include "QMCWaveFunctions/WaveFunctionComponent.h" -#include "QMCWaveFunctions/EinsplineSetBuilder.h" -#include "QMCWaveFunctions/RotatedSPOs.h" -#include "checkMatrix.hpp" -#include "FakeSPO.h" -#include -#include "QMCWaveFunctions/VariableSet.h" - -#include -#include -#include - -using std::string; - -namespace qmcplusplus -{ -/* - JPT 04.01.2022: Adapted from test_einset.cpp - Test the spline rotated machinery for SplineR2R (extend to others later). -*/ -TEST_CASE("RotatedSPOs via SplineR2R", "[wavefunction]") -{ - using RealType = QMCTraits::RealType; - - /* - BEGIN Boilerplate stuff to make a simple SPOSet. Copied from test_einset.cpp - */ - - Communicate* c = OHMMS::Controller; - - // We get a "Mismatched supercell lattices" error due to default ctor? - ParticleSet::ParticleLayout lattice; - - // diamondC_1x1x1 - lattice.R = {3.37316115, 3.37316115, 0.0, 0.0, 3.37316115, 3.37316115, 3.37316115, 0.0, 3.37316115}; - - ParticleSetPool ptcl = ParticleSetPool(c); - ptcl.setSimulationCell(lattice); - // LAttice seems fine after this point... - - auto ions_uptr = std::make_unique(ptcl.getSimulationCell()); - auto elec_uptr = std::make_unique(ptcl.getSimulationCell()); - ParticleSet& ions_(*ions_uptr); - ParticleSet& elec_(*elec_uptr); - - ions_.setName("ion"); - ptcl.addParticleSet(std::move(ions_uptr)); - ions_.create({2}); - ions_.R[0] = {0.0, 0.0, 0.0}; - ions_.R[1] = {1.68658058, 1.68658058, 1.68658058}; - elec_.setName("elec"); - ptcl.addParticleSet(std::move(elec_uptr)); - elec_.create({2}); - elec_.R[0] = {0.0, 0.0, 0.0}; - elec_.R[1] = {0.0, 1.0, 0.0}; - SpeciesSet& tspecies = elec_.getSpeciesSet(); - int upIdx = tspecies.addSpecies("u"); - int chargeIdx = tspecies.addAttribute("charge"); - tspecies(chargeIdx, upIdx) = -1; - - //diamondC_1x1x1 - 8 bands available - const char* particles = R"( - - -)"; - - Libxml2Document doc; - bool okay = doc.parseFromString(particles); - REQUIRE(okay); - - xmlNodePtr root = doc.getRoot(); - - xmlNodePtr ein1 = xmlFirstElementChild(root); - - EinsplineSetBuilder einSet(elec_, ptcl.getPool(), c, ein1); - auto spo = einSet.createSPOSetFromXML(ein1); - REQUIRE(spo); - - /* - END Boilerplate stuff. Now we have a SplineR2R wavefunction - ready for rotation. What follows is the actual test. - */ - - // SplineR2R only for the moment, so skip if QMC_COMPLEX is set -#if !defined(QMC_COMPLEX) - - spo->storeParamsBeforeRotation(); - // 1.) Make a RotatedSPOs object so that we can use the rotation routines - auto rot_spo = std::make_unique("one_rotated_set", std::move(spo)); - - // Sanity check for orbs. Expect 2 electrons, 8 orbitals, & 79507 coefs/orb. - const auto orbitalsetsize = rot_spo->getOrbitalSetSize(); - REQUIRE(orbitalsetsize == 8); - - // 2.) Get data for unrotated orbitals. Check that there's no rotation - rot_spo->buildOptVariables(elec_.R.size()); - SPOSet::ValueMatrix psiM_bare(elec_.R.size(), orbitalsetsize); - SPOSet::GradMatrix dpsiM_bare(elec_.R.size(), orbitalsetsize); - SPOSet::ValueMatrix d2psiM_bare(elec_.R.size(), orbitalsetsize); - rot_spo->evaluate_notranspose(elec_, 0, elec_.R.size(), psiM_bare, dpsiM_bare, d2psiM_bare); - - // This stuff checks that no rotation was applied. Copied from test_einset.cpp. - // value - CHECK(std::real(psiM_bare[1][0]) == Approx(-0.8886948824)); - CHECK(std::real(psiM_bare[1][1]) == Approx(1.4194120169)); - // grad - CHECK(std::real(dpsiM_bare[1][0][0]) == Approx(-0.0000183403)); - CHECK(std::real(dpsiM_bare[1][0][1]) == Approx(0.1655139178)); - CHECK(std::real(dpsiM_bare[1][0][2]) == Approx(-0.0000193077)); - CHECK(std::real(dpsiM_bare[1][1][0]) == Approx(-1.3131694794)); - CHECK(std::real(dpsiM_bare[1][1][1]) == Approx(-1.1174004078)); - CHECK(std::real(dpsiM_bare[1][1][2]) == Approx(-0.8462534547)); - // lapl - CHECK(std::real(d2psiM_bare[1][0]) == Approx(1.3313053846)); - CHECK(std::real(d2psiM_bare[1][1]) == Approx(-4.712583065)); - - /* - 3.) Apply a rotation to the orbitals - To do this, construct a params vector and call the - RotatedSPOs::apply_rotation(params) method. That should do the - right thing for this particular spline class. - - For 2 electrons in 8 orbs, we expect 2*(8-2) = 12 params. - */ - const auto rot_size = rot_spo->m_act_rot_inds.size(); - REQUIRE(rot_size == 12); // = Nelec*(Norbs - Nelec) = 2*(8-2) = 12 - std::vector param(rot_size); - for (auto i = 0; i < rot_size; i++) - { - param[i] = 0.01 * static_cast(i); - } - rot_spo->apply_rotation(param, false); // Expect this to call SplineR2R::applyRotation() - - // 4.) Get data for rotated orbitals. - SPOSet::ValueMatrix psiM_rot(elec_.R.size(), orbitalsetsize); - SPOSet::GradMatrix dpsiM_rot(elec_.R.size(), orbitalsetsize); - SPOSet::ValueMatrix d2psiM_rot(elec_.R.size(), orbitalsetsize); - rot_spo->evaluate_notranspose(elec_, 0, elec_.R.size(), psiM_rot, dpsiM_rot, d2psiM_rot); - - /* - Manually encode the unitary transformation. Ugly, but it works. - @TODO: Use the total rotation machinery when it's implemented - - NB: This is truncated to 5 sig-figs, so there is some slop here as - compared to what is done in the splines via apply_rotation(). - So below we reduce the threshold for comparison. This can - probably be ditched once we have a way to grab the actual - rotation matrix... - */ - SPOSet::ValueMatrix rot_mat(orbitalsetsize, orbitalsetsize); - rot_mat[0][0] = 0.99726; - rot_mat[0][1] = -0.00722; - rot_mat[0][2] = 0.00014; - rot_mat[0][3] = -0.00982; - rot_mat[0][4] = -0.01979; - rot_mat[0][5] = -0.02976; - rot_mat[0][6] = -0.03972; - rot_mat[0][7] = -0.04969; - rot_mat[1][0] = -0.00722; - rot_mat[1][1] = 0.97754; - rot_mat[1][2] = -0.05955; - rot_mat[1][3] = -0.06945; - rot_mat[1][4] = -0.07935; - rot_mat[1][5] = -0.08925; - rot_mat[1][6] = -0.09915; - rot_mat[1][7] = -0.10905; - rot_mat[2][0] = -0.00014; - rot_mat[2][1] = 0.05955; - rot_mat[2][2] = 0.99821; - rot_mat[2][3] = -0.00209; - rot_mat[2][4] = -0.00239; - rot_mat[2][5] = -0.00269; - rot_mat[2][6] = -0.00299; - rot_mat[2][7] = -0.00329; - rot_mat[3][0] = 0.00982; - rot_mat[3][1] = 0.06945; - rot_mat[3][2] = -0.00209; - rot_mat[3][3] = 0.99751; - rot_mat[3][4] = -0.00289; - rot_mat[3][5] = -0.00329; - rot_mat[3][6] = -0.00368; - rot_mat[3][7] = -0.00408; - rot_mat[4][0] = 0.01979; - rot_mat[4][1] = 0.07935; - rot_mat[4][2] = -0.00239; - rot_mat[4][3] = -0.00289; - rot_mat[4][4] = 0.99661; - rot_mat[4][5] = -0.00388; - rot_mat[4][6] = -0.00438; - rot_mat[4][7] = -0.00488; - rot_mat[5][0] = 0.02976; - rot_mat[5][1] = 0.08925; - rot_mat[5][2] = -0.00269; - rot_mat[5][3] = -0.00329; - rot_mat[5][4] = -0.00388; - rot_mat[5][5] = 0.99552; - rot_mat[5][6] = -0.00508; - rot_mat[5][7] = -0.00568; - rot_mat[6][0] = 0.03972; - rot_mat[6][1] = 0.09915; - rot_mat[6][2] = -0.00299; - rot_mat[6][3] = -0.00368; - rot_mat[6][4] = -0.00438; - rot_mat[6][5] = -0.00508; - rot_mat[6][6] = 0.99422; - rot_mat[6][7] = -0.00647; - rot_mat[7][0] = 0.04969; - rot_mat[7][1] = 0.10905; - rot_mat[7][2] = -0.00329; - rot_mat[7][3] = -0.00408; - rot_mat[7][4] = -0.00488; - rot_mat[7][5] = -0.00568; - rot_mat[7][6] = -0.00647; - rot_mat[7][7] = 0.99273; - - // Now compute the expected values by hand using the transformation above - double val1 = 0.; - double val2 = 0.; - for (auto i = 0; i < rot_mat.size1(); i++) - { - val1 += psiM_bare[0][i] * rot_mat[i][0]; - val2 += psiM_bare[1][i] * rot_mat[i][0]; - } - - // value - CHECK(std::real(psiM_rot[0][0]) == Approx(val1)); - CHECK(std::real(psiM_rot[1][0]) == Approx(val2)); - - std::vector grad1(3); - std::vector grad2(3); - for (auto j = 0; j < grad1.size(); j++) - { - for (auto i = 0; i < rot_mat.size1(); i++) - { - grad1[j] += dpsiM_bare[0][i][j] * rot_mat[i][0]; - grad2[j] += dpsiM_bare[1][i][j] * rot_mat[i][0]; - } - } - - // grad - CHECK(dpsiM_rot[0][0][0] == Approx(grad1[0]).epsilon(0.0001)); - CHECK(dpsiM_rot[0][0][1] == Approx(grad1[1]).epsilon(0.0001)); - CHECK(dpsiM_rot[0][0][2] == Approx(grad1[2]).epsilon(0.0001)); - CHECK(dpsiM_rot[1][0][0] == Approx(grad2[0]).epsilon(0.0001)); - CHECK(dpsiM_rot[1][0][1] == Approx(grad2[1]).epsilon(0.0001)); - CHECK(dpsiM_rot[1][0][2] == Approx(grad2[2]).epsilon(0.0001)); - - double lap1 = 0.; - double lap2 = 0.; - for (auto i = 0; i < rot_mat.size1(); i++) - { - lap1 += d2psiM_bare[0][i] * rot_mat[i][0]; - lap2 += d2psiM_bare[1][i] * rot_mat[i][0]; - } - - // Lapl - CHECK(std::real(d2psiM_rot[0][0]) == Approx(lap1).epsilon(0.0001)); - CHECK(std::real(d2psiM_rot[1][0]) == Approx(lap2).epsilon(0.0001)); - -#endif -} - -TEST_CASE("RotatedSPOs createRotationIndices", "[wavefunction]") -{ - // No active-active or virtual-virtual rotations - // Only active-virtual - RotatedSPOs::RotationIndices rot_ind; - int nel = 1; - int nmo = 3; - RotatedSPOs::createRotationIndices(nel, nmo, rot_ind); - CHECK(rot_ind.size() == 2); - - // Full rotation contains all rotations - // Size should be number of pairs of orbitals: nmo*(nmo-1)/2 - RotatedSPOs::RotationIndices full_rot_ind; - RotatedSPOs::createRotationIndicesFull(nel, nmo, full_rot_ind); - CHECK(full_rot_ind.size() == 3); - - nel = 2; - RotatedSPOs::RotationIndices rot_ind2; - RotatedSPOs::createRotationIndices(nel, nmo, rot_ind2); - CHECK(rot_ind2.size() == 2); - - RotatedSPOs::RotationIndices full_rot_ind2; - RotatedSPOs::createRotationIndicesFull(nel, nmo, full_rot_ind2); - CHECK(full_rot_ind2.size() == 3); - - nmo = 4; - RotatedSPOs::RotationIndices rot_ind3; - RotatedSPOs::createRotationIndices(nel, nmo, rot_ind3); - CHECK(rot_ind3.size() == 4); - - RotatedSPOs::RotationIndices full_rot_ind3; - RotatedSPOs::createRotationIndicesFull(nel, nmo, full_rot_ind3); - CHECK(full_rot_ind3.size() == 6); -} - -TEST_CASE("RotatedSPOs constructAntiSymmetricMatrix", "[wavefunction]") -{ - using ValueType = SPOSet::ValueType; - using ValueMatrix = SPOSet::ValueMatrix; - - RotatedSPOs::RotationIndices rot_ind; - int nel = 1; - int nmo = 3; - RotatedSPOs::createRotationIndices(nel, nmo, rot_ind); - - ValueMatrix m3(nmo, nmo); - m3 = ValueType(0); - std::vector params = {0.1, 0.2}; - - RotatedSPOs::constructAntiSymmetricMatrix(rot_ind, params, m3); - - // clang-format off - std::vector expected_data = { 0.0, -0.1, -0.2, - 0.1, 0.0, 0.0, - 0.2, 0.0, 0.0 }; - // clang-format on - - ValueMatrix expected_m3(expected_data.data(), 3, 3); - - CheckMatrixResult check_matrix_result = checkMatrix(m3, expected_m3, true); - CHECKED_ELSE(check_matrix_result.result) { FAIL(check_matrix_result.result_message); } - - std::vector params_out(2); - RotatedSPOs::extractParamsFromAntiSymmetricMatrix(rot_ind, m3, params_out); - CHECK(params_out[0] == Approx(0.1)); - CHECK(params_out[1] == Approx(0.2)); -} - -// Expected values of the matrix exponential come from gen_matrix_ops.py -TEST_CASE("RotatedSPOs exponentiate matrix", "[wavefunction]") -{ - using ValueType = SPOSet::ValueType; - using ValueMatrix = SPOSet::ValueMatrix; - - std::vector mat1_data = {0.0}; - SPOSet::ValueMatrix m1(mat1_data.data(), 1, 1); - RotatedSPOs::exponentiate_antisym_matrix(m1); - // Always return 1.0 (the only possible anti-symmetric 1x1 matrix is 0) - CHECK(m1(0, 0) == ValueApprox(1.0)); - - // clang-format off - std::vector mat2_data = { 0.0, -0.1, - 0.1, 0.0 }; - // clang-format on - - SPOSet::ValueMatrix m2(mat2_data.data(), 2, 2); - RotatedSPOs::exponentiate_antisym_matrix(m2); - - // clang-format off - std::vector expected_rot2 = { 0.995004165278026, -0.0998334166468282, - 0.0998334166468282, 0.995004165278026 }; - // clang-format on - - ValueMatrix expected_m2(expected_rot2.data(), 2, 2); - CheckMatrixResult check_matrix_result2 = checkMatrix(m2, expected_m2, true); - CHECKED_ELSE(check_matrix_result2.result) { FAIL(check_matrix_result2.result_message); } - - - // clang-format off - std::vector m3_input_data = { 0.0, -0.3, -0.1, - 0.3, 0.0, -0.2, - 0.1, 0.2, 0.0 }; - - - std::vector expected_rot3 = { 0.950580617906092, -0.302932713402637, -0.0680313164049401, - 0.283164960565074, 0.935754803277919, -0.210191705950743, - 0.127334574917630, 0.180540076694398, 0.975290308953046 }; - - // clang-format on - - ValueMatrix m3(m3_input_data.data(), 3, 3); - ValueMatrix expected_m3(expected_rot3.data(), 3, 3); - - RotatedSPOs::exponentiate_antisym_matrix(m3); - - CheckMatrixResult check_matrix_result3 = checkMatrix(m3, expected_m3, true); - CHECKED_ELSE(check_matrix_result3.result) { FAIL(check_matrix_result3.result_message); } -} - -TEST_CASE("RotatedSPOs log matrix", "[wavefunction]") -{ - using ValueType = SPOSet::ValueType; - using ValueMatrix = SPOSet::ValueMatrix; - - std::vector mat1_data = {1.0}; - SPOSet::ValueMatrix m1(mat1_data.data(), 1, 1); - SPOSet::ValueMatrix out_m1(1, 1); - RotatedSPOs::log_antisym_matrix(m1, out_m1); - // Should always be 1.0 (the only possible anti-symmetric 1x1 matrix is 0) - CHECK(out_m1(0, 0) == ValueApprox(0.0)); - - // clang-format off - std::vector start_rot2 = { 0.995004165278026, -0.0998334166468282, - 0.0998334166468282, 0.995004165278026 }; - - std::vector mat2_data = { 0.0, -0.1, - 0.1, 0.0 }; - // clang-format on - - ValueMatrix rot_m2(start_rot2.data(), 2, 2); - ValueMatrix out_m2(2, 2); - RotatedSPOs::log_antisym_matrix(rot_m2, out_m2); - - SPOSet::ValueMatrix m2(mat2_data.data(), 2, 2); - CheckMatrixResult check_matrix_result2 = checkMatrix(m2, out_m2, true); - CHECKED_ELSE(check_matrix_result2.result) { FAIL(check_matrix_result2.result_message); } - - // clang-format off - std::vector start_rot3 = { 0.950580617906092, -0.302932713402637, -0.0680313164049401, - 0.283164960565074, 0.935754803277919, -0.210191705950743, - 0.127334574917630, 0.180540076694398, 0.975290308953046 }; - - std::vector m3_input_data = { 0.0, -0.3, -0.1, - 0.3, 0.0, -0.2, - 0.1, 0.2, 0.0 }; - // clang-format on - ValueMatrix rot_m3(start_rot3.data(), 3, 3); - ValueMatrix out_m3(3, 3); - RotatedSPOs::log_antisym_matrix(rot_m3, out_m3); - - SPOSet::ValueMatrix m3(m3_input_data.data(), 3, 3); - CheckMatrixResult check_matrix_result3 = checkMatrix(m3, out_m3, true); - CHECKED_ELSE(check_matrix_result3.result) { FAIL(check_matrix_result3.result_message); } -} - -// Test round trip A -> exp(A) -> log(exp(A)) -// The log is multi-valued so this test may fail if the rotation parameters are too large. -// The exponentials will be the same, though -// exp(log(exp(A))) == exp(A) -TEST_CASE("RotatedSPOs exp-log matrix", "[wavefunction]") -{ - using ValueType = SPOSet::ValueType; - using ValueMatrix = SPOSet::ValueMatrix; - - RotatedSPOs::RotationIndices rot_ind; - int nel = 2; - int nmo = 4; - RotatedSPOs::createRotationIndices(nel, nmo, rot_ind); - - ValueMatrix rot_m4(nmo, nmo); - rot_m4 = ValueType(0); - - std::vector params4 = {-1.1, 1.5, 0.2, -0.15}; - - RotatedSPOs::constructAntiSymmetricMatrix(rot_ind, params4, rot_m4); - ValueMatrix orig_rot_m4 = rot_m4; - ValueMatrix out_m4(nmo, nmo); - - RotatedSPOs::exponentiate_antisym_matrix(rot_m4); - - RotatedSPOs::log_antisym_matrix(rot_m4, out_m4); - - CheckMatrixResult check_matrix_result4 = checkMatrix(out_m4, orig_rot_m4, true); - CHECKED_ELSE(check_matrix_result4.result) { FAIL(check_matrix_result4.result_message); } - - std::vector params4out(4); - RotatedSPOs::extractParamsFromAntiSymmetricMatrix(rot_ind, out_m4, params4out); - for (int i = 0; i < params4.size(); i++) - { - CHECK(params4[i] == Approx(params4out[i])); - } -} - -TEST_CASE("RotatedSPOs hcpBe", "[wavefunction]") -{ - using RealType = QMCTraits::RealType; - Communicate* c = OHMMS::Controller; - - ParticleSet::ParticleLayout lattice; - lattice.R = {4.32747284, 0.00000000, 0.00000000, -2.16373642, 3.74770142, - 0.00000000, 0.00000000, 0.00000000, 6.78114995}; - - ParticleSetPool ptcl = ParticleSetPool(c); - ptcl.setSimulationCell(lattice); - auto ions_uptr = std::make_unique(ptcl.getSimulationCell()); - auto elec_uptr = std::make_unique(ptcl.getSimulationCell()); - ParticleSet& ions(*ions_uptr); - ParticleSet& elec(*elec_uptr); - - ions.setName("ion"); - ptcl.addParticleSet(std::move(ions_uptr)); - ions.create({1}); - ions.R[0] = {0.0, 0.0, 0.0}; - - elec.setName("elec"); - ptcl.addParticleSet(std::move(elec_uptr)); - elec.create({1}); - elec.R[0] = {0.0, 0.0, 0.0}; - - SpeciesSet& tspecies = elec.getSpeciesSet(); - int upIdx = tspecies.addSpecies("u"); - int chargeIdx = tspecies.addAttribute("charge"); - tspecies(chargeIdx, upIdx) = -1; - - // Add the attribute save_coefs="yes" to the sposet_builder tag to generate the - // spline file for use in eval_bspline_spo.py - - const char* particles = R"( - - - -)"; - - Libxml2Document doc; - bool okay = doc.parseFromString(particles); - REQUIRE(okay); - - xmlNodePtr root = doc.getRoot(); - - xmlNodePtr sposet_builder = xmlFirstElementChild(root); - xmlNodePtr sposet_ptr = xmlFirstElementChild(sposet_builder); - - EinsplineSetBuilder einSet(elec, ptcl.getPool(), c, sposet_builder); - auto spo = einSet.createSPOSetFromXML(sposet_ptr); - REQUIRE(spo); - - spo->storeParamsBeforeRotation(); - auto rot_spo = std::make_unique("one_rotated_set", std::move(spo)); - - // Sanity check for orbs. Expect 1 electron, 2 orbitals - const auto orbitalsetsize = rot_spo->getOrbitalSetSize(); - REQUIRE(orbitalsetsize == 2); - - rot_spo->buildOptVariables(elec.R.size()); - - SPOSet::ValueMatrix psiM_bare(elec.R.size(), orbitalsetsize); - SPOSet::GradMatrix dpsiM_bare(elec.R.size(), orbitalsetsize); - SPOSet::ValueMatrix d2psiM_bare(elec.R.size(), orbitalsetsize); - rot_spo->evaluate_notranspose(elec, 0, elec.R.size(), psiM_bare, dpsiM_bare, d2psiM_bare); - - // Values generated from eval_bspline_spo.py, the generate_point_values_hcpBe function - CHECK(std::real(psiM_bare[0][0]) == Approx(0.210221765375514)); - CHECK(std::real(psiM_bare[0][1]) == Approx(-2.984345024542937e-06)); - - CHECK(std::real(d2psiM_bare[0][0]) == Approx(5.303848362116568)); - - opt_variables_type opt_vars; - rot_spo->checkInVariablesExclusive(opt_vars); - opt_vars.resetIndex(); - rot_spo->checkOutVariables(opt_vars); - rot_spo->resetParametersExclusive(opt_vars); - - using ValueType = QMCTraits::ValueType; - Vector dlogpsi(1); - Vector dhpsioverpsi(1); - rot_spo->evaluateDerivatives(elec, opt_vars, dlogpsi, dhpsioverpsi, 0, 1); - - CHECK(dlogpsi[0] == ValueApprox(-1.41961753e-05)); - CHECK(dhpsioverpsi[0] == ValueApprox(-0.00060853)); - - std::vector params = {0.1}; - rot_spo->apply_rotation(params, false); - - rot_spo->evaluate_notranspose(elec, 0, elec.R.size(), psiM_bare, dpsiM_bare, d2psiM_bare); - CHECK(std::real(psiM_bare[0][0]) == Approx(0.20917123424337608)); - CHECK(std::real(psiM_bare[0][1]) == Approx(-0.02099012652669549)); - - CHECK(std::real(d2psiM_bare[0][0]) == Approx(5.277362065087747)); - - dlogpsi[0] = 0.0; - dhpsioverpsi[0] = 0.0; - - rot_spo->evaluateDerivatives(elec, opt_vars, dlogpsi, dhpsioverpsi, 0, 1); - CHECK(dlogpsi[0] == ValueApprox(-0.10034901119468914)); - CHECK(dhpsioverpsi[0] == ValueApprox(32.96939041498753)); -} - -// Test construction of delta rotation -TEST_CASE("RotatedSPOs construct delta matrix", "[wavefunction]") -{ - using ValueType = SPOSet::ValueType; - using ValueMatrix = SPOSet::ValueMatrix; - - int nel = 2; - int nmo = 4; - RotatedSPOs::RotationIndices rot_ind; - RotatedSPOs::createRotationIndices(nel, nmo, rot_ind); - RotatedSPOs::RotationIndices full_rot_ind; - RotatedSPOs::createRotationIndicesFull(nel, nmo, full_rot_ind); - // rot_ind size is 4 and full rot_ind size is 6 - - ValueMatrix rot_m4(nmo, nmo); - rot_m4 = ValueType(0); - - // When comparing with gen_matrix_ops.py, be aware of the order of indices - // in full_rot - // rot_ind is (0,2) (0,3) (1,2) (1,3) - // full_rot_ind is (0,2) (0,3) (1,2) (1,3) (0,1) (2,3) - // The extra indices go at the back - std::vector old_params = {1.5, 0.2, -0.15, 0.03, -1.1, 0.05}; - std::vector delta_params = {0.1, 0.3, 0.2, -0.1}; - std::vector new_params(6); - - RotatedSPOs::constructDeltaRotation(delta_params, old_params, rot_ind, full_rot_ind, new_params, rot_m4); - - // clang-format off - std::vector rot_data4 = - { -0.371126931484737, 0.491586564957393, -0.784780958819798, 0.0687480658200083, - -0.373372784561548, 0.66111547793048, 0.610450337985578, 0.225542620014052, - 0.751270334458895, 0.566737323353515, -0.0297901110611425, -0.336918744155143, - 0.398058348785074, 0.00881931472604944, -0.102867783149713, 0.911531672428406 }; - // clang-format on - - ValueMatrix new_rot_m4(rot_data4.data(), 4, 4); - - CheckMatrixResult check_matrix_result4 = checkMatrix(rot_m4, new_rot_m4, true); - CHECKED_ELSE(check_matrix_result4.result) { FAIL(check_matrix_result4.result_message); } - - // Reminder: Ordering! - std::vector expected_new_param = {1.6813965019790489, 0.3623564254653294, -0.05486544454559908, - -0.20574472941408453, -0.9542513302873077, 0.27497788909911774}; - for (int i = 0; i < new_params.size(); i++) - CHECK(new_params[i] == Approx(expected_new_param[i])); - - - // Rotated back to original position - - std::vector new_params2(6); - std::vector reverse_delta_params = {-0.1, -0.3, -0.2, 0.1}; - RotatedSPOs::constructDeltaRotation(reverse_delta_params, new_params, rot_ind, full_rot_ind, new_params2, rot_m4); - for (int i = 0; i < new_params2.size(); i++) - CHECK(new_params2[i] == Approx(old_params[i])); -} - -namespace testing -{ -opt_variables_type& getMyVars(SPOSet& rot) { return rot.myVars; } -opt_variables_type& getMyVarsFull(RotatedSPOs& rot) { return rot.myVarsFull; } -std::vector>& getHistoryParams(RotatedSPOs& rot) { return rot.history_params_; } -} // namespace testing - -// Test using global rotation -TEST_CASE("RotatedSPOs read and write parameters", "[wavefunction]") -{ - auto fake_spo = std::make_unique(); - fake_spo->setOrbitalSetSize(4); - RotatedSPOs rot("fake_rot", std::move(fake_spo)); - int nel = 2; - rot.buildOptVariables(nel); - - optimize::VariableSet vs; - rot.checkInVariablesExclusive(vs); - vs[0] = 0.1; - vs[1] = 0.15; - vs[2] = 0.2; - vs[3] = 0.25; - rot.resetParametersExclusive(vs); - - { - hdf_archive hout; - vs.writeToHDF("rot_vp.h5", hout); - - rot.writeVariationalParameters(hout); - } - - auto fake_spo2 = std::make_unique(); - fake_spo2->setOrbitalSetSize(4); - - RotatedSPOs rot2("fake_rot", std::move(fake_spo2)); - rot2.buildOptVariables(nel); - - optimize::VariableSet vs2; - rot2.checkInVariablesExclusive(vs2); - - hdf_archive hin; - vs2.readFromHDF("rot_vp.h5", hin); - rot2.readVariationalParameters(hin); - - opt_variables_type& var = testing::getMyVars(rot2); - CHECK(var[0] == Approx(vs[0])); - CHECK(var[1] == Approx(vs[1])); - CHECK(var[2] == Approx(vs[2])); - CHECK(var[3] == Approx(vs[3])); - - opt_variables_type& full_var = testing::getMyVarsFull(rot2); - CHECK(full_var[0] == Approx(vs[0])); - CHECK(full_var[1] == Approx(vs[1])); - CHECK(full_var[2] == Approx(vs[2])); - CHECK(full_var[3] == Approx(vs[3])); - CHECK(full_var[4] == Approx(0.0)); - CHECK(full_var[5] == Approx(0.0)); -} - -// Test using history list. -TEST_CASE("RotatedSPOs read and write parameters history", "[wavefunction]") -{ - auto fake_spo = std::make_unique(); - fake_spo->setOrbitalSetSize(4); - RotatedSPOs rot("fake_rot", std::move(fake_spo)); - rot.set_use_global_rotation(false); - int nel = 2; - rot.buildOptVariables(nel); - - optimize::VariableSet vs; - rot.checkInVariablesExclusive(vs); - vs[0] = 0.1; - vs[1] = 0.15; - vs[2] = 0.2; - vs[3] = 0.25; - rot.resetParametersExclusive(vs); - - { - hdf_archive hout; - vs.writeToHDF("rot_vp_hist.h5", hout); - - rot.writeVariationalParameters(hout); - } - - auto fake_spo2 = std::make_unique(); - fake_spo2->setOrbitalSetSize(4); - - RotatedSPOs rot2("fake_rot", std::move(fake_spo2)); - rot2.buildOptVariables(nel); - - optimize::VariableSet vs2; - rot2.checkInVariablesExclusive(vs2); - - hdf_archive hin; - vs2.readFromHDF("rot_vp_hist.h5", hin); - rot2.readVariationalParameters(hin); - - opt_variables_type& var = testing::getMyVars(rot2); - CHECK(var[0] == Approx(vs[0])); - CHECK(var[1] == Approx(vs[1])); - CHECK(var[2] == Approx(vs[2])); - CHECK(var[3] == Approx(vs[3])); - - auto hist = testing::getHistoryParams(rot2); - REQUIRE(hist.size() == 1); - REQUIRE(hist[0].size() == 4); -} - -class DummySPOSetWithoutMW : public SPOSet -{ -public: - DummySPOSetWithoutMW(const std::string& my_name) : SPOSet(my_name) {} - void setOrbitalSetSize(int norbs) override {} - void evaluateValue(const ParticleSet& P, int iat, SPOSet::ValueVector& psi) override - { - assert(psi.size() == 3); - psi[0] = 123; - psi[1] = 456; - psi[2] = 789; - } - void evaluateVGL(const ParticleSet& P, int iat, ValueVector& psi, GradVector& dpsi, ValueVector& d2psi) override {} - void evaluate_notranspose(const ParticleSet& P, - int first, - int last, - ValueMatrix& logdet, - GradMatrix& dlogdet, - ValueMatrix& d2logdet) override - {} - std::string getClassName() const override { return my_name_; } -}; - -class DummySPOSetWithMW : public DummySPOSetWithoutMW -{ -public: - DummySPOSetWithMW(const std::string& my_name) : DummySPOSetWithoutMW(my_name) {} - void mw_evaluateValue(const RefVectorWithLeader& spo_list, - const RefVectorWithLeader& P_list, - int iat, - const RefVector& psi_v_list) const override - { - for (auto& psi : psi_v_list) - { - assert(psi.get().size() == 3); - psi.get()[0] = 321; - psi.get()[1] = 654; - psi.get()[2] = 987; - } - } -}; - -TEST_CASE("RotatedSPOs mw_ APIs", "[wavefunction]") -{ - //checking that mw_ API works in RotatedSPOs and is not defaulting to - //SPOSet default implementation - { - //First check calling the mw_ APIs for RotatedSPOs, for which the - //underlying implementation just calls the underlying SPOSet mw_ API - //In the case that the underlying SPOSet doesn't specialize the mw_ API, - //the underlying SPOSet will fall back to the default SPOSet mw_, which is - //just a loop over the single walker API. - RotatedSPOs rot_spo0("rotated0", std::make_unique("no mw 0")); - RotatedSPOs rot_spo1("rotated1", std::make_unique("no mw 1")); - RefVectorWithLeader spo_list(rot_spo0, {rot_spo0, rot_spo1}); - - ResourceCollection spo_res("test_rot_res"); - rot_spo0.createResource(spo_res); - ResourceCollectionTeamLock mw_sposet_lock(spo_res, spo_list); - - const SimulationCell simulation_cell; - ParticleSet elec0(simulation_cell); - ParticleSet elec1(simulation_cell); - RefVectorWithLeader p_list(elec0, {elec0, elec1}); - - SPOSet::ValueVector psi0(3); - SPOSet::ValueVector psi1(3); - RefVector psi_v_list{psi0, psi1}; - - rot_spo0.mw_evaluateValue(spo_list, p_list, 0, psi_v_list); - for (int iw = 0; iw < spo_list.size(); iw++) - { - CHECK(psi_v_list[iw].get()[0] == Approx(123)); - CHECK(psi_v_list[iw].get()[1] == Approx(456)); - CHECK(psi_v_list[iw].get()[2] == Approx(789)); - } - } - { - //In the case that the underlying SPOSet DOES have mw_ specializations, - //we want to make sure that RotatedSPOs are triggering that appropriately - //This will mean that the underlying SPOSets will do the appropriate offloading - //To check this, DummySPOSetWithMW has an explicit mw_evaluateValue which sets - //different values than what gets set in evaluateValue. By doing this, - //we are ensuring that RotatedSPOs->mw_evaluaeValue is calling the specialization - //in the underlying SPO and not using the default SPOSet implementation which - //loops over single walker APIs (which have different values enforced in - // DummySPOSetWithoutMW - - RotatedSPOs rot_spo0("rotated0", std::make_unique("mw 0")); - RotatedSPOs rot_spo1("rotated1", std::make_unique("mw 1")); - RefVectorWithLeader spo_list(rot_spo0, {rot_spo0, rot_spo1}); - - ResourceCollection spo_res("test_rot_res"); - rot_spo0.createResource(spo_res); - ResourceCollectionTeamLock mw_sposet_lock(spo_res, spo_list); - - const SimulationCell simulation_cell; - ParticleSet elec0(simulation_cell); - ParticleSet elec1(simulation_cell); - RefVectorWithLeader p_list(elec0, {elec0, elec1}); - - SPOSet::ValueVector psi0(3); - SPOSet::ValueVector psi1(3); - RefVector psi_v_list{psi0, psi1}; - - rot_spo0.mw_evaluateValue(spo_list, p_list, 0, psi_v_list); - for (int iw = 0; iw < spo_list.size(); iw++) - { - CHECK(psi_v_list[iw].get()[0] == Approx(321)); - CHECK(psi_v_list[iw].get()[1] == Approx(654)); - CHECK(psi_v_list[iw].get()[2] == Approx(987)); - } - } -} - -} // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/tests/test_RotatedSPOsT.cpp b/src/QMCWaveFunctions/tests/test_RotatedSPOsT.cpp index 5d2a3d0bee..e27ca3524e 100644 --- a/src/QMCWaveFunctions/tests/test_RotatedSPOsT.cpp +++ b/src/QMCWaveFunctions/tests/test_RotatedSPOsT.cpp @@ -28,6 +28,7 @@ #include #include +#include using std::string; @@ -49,10 +50,26 @@ using ValueApprox = typename ValueApproxHelper::Type; namespace testing { -OptVariablesTypeT& getMyVars(SPOSetT& rot) { return rot.myVars; } -OptVariablesTypeT& getMyVars(SPOSetT& rot) { return rot.myVars; } -OptVariablesTypeT& getMyVarsFull(RotatedSPOsT& rot) { return rot.myVarsFull; } -OptVariablesTypeT& getMyVarsFull(RotatedSPOsT& rot) { return rot.myVarsFull; } +OptVariablesTypeT& +getMyVars(SPOSetT& rot) +{ + return rot.myVars; +} +OptVariablesTypeT& +getMyVars(SPOSetT& rot) +{ + return rot.myVars; +} +OptVariablesTypeT& +getMyVarsFull(RotatedSPOsT& rot) +{ + return rot.myVarsFull; +} +OptVariablesTypeT& +getMyVarsFull(RotatedSPOsT& rot) +{ + return rot.myVarsFull; +} std::vector>& getHistoryParams(RotatedSPOsT& rot) { @@ -66,12 +83,14 @@ getHistoryParams(RotatedSPOsT& rot) } } // namespace testing +using TestTypeList = std::tuple; + /* JPT 04.01.2022: Adapted from test_einset.cpp Test the spline rotated machinery for SplineR2R (extend to others later). */ -TEMPLATE_TEST_CASE( - "RotatedSPOs via SplineR2R", "[wavefunction][template]", double, float) +TEMPLATE_LIST_TEST_CASE( + "RotatedSPOs via SplineR2R", "[wavefunction][template]", TestTypeList) { using RealType = typename SPOSetT::RealType; @@ -337,8 +356,8 @@ TEMPLATE_TEST_CASE( #endif } -TEMPLATE_TEST_CASE("RotatedSPOs createRotationIndices", - "[wavefunction][template]", double, float) +TEMPLATE_LIST_TEST_CASE("RotatedSPOs createRotationIndices", + "[wavefunction][template]", TestTypeList) { // No active-active or virtual-virtual rotations // Only active-virtual @@ -373,8 +392,8 @@ TEMPLATE_TEST_CASE("RotatedSPOs createRotationIndices", CHECK(full_rot_ind3.size() == 6); } -TEMPLATE_TEST_CASE("RotatedSPOs constructAntiSymmetricMatrix", - "[wavefunction][template]", double, float) +TEMPLATE_LIST_TEST_CASE("RotatedSPOs constructAntiSymmetricMatrix", + "[wavefunction][template]", TestTypeList) { using ValueType = typename SPOSetT::ValueType; using ValueMatrix = typename SPOSetT::ValueMatrix; @@ -412,8 +431,8 @@ TEMPLATE_TEST_CASE("RotatedSPOs constructAntiSymmetricMatrix", } // Expected values of the matrix exponential come from gen_matrix_ops.py -TEMPLATE_TEST_CASE("RotatedSPOs exponentiate matrix", - "[wavefunction][template]", double, float) +TEMPLATE_LIST_TEST_CASE( + "RotatedSPOs exponentiate matrix", "[wavefunction][template]", TestTypeList) { using ValueType = typename SPOSetT::ValueType; using ValueMatrix = typename SPOSetT::ValueMatrix; @@ -468,8 +487,8 @@ TEMPLATE_TEST_CASE("RotatedSPOs exponentiate matrix", } } -TEMPLATE_TEST_CASE( - "RotatedSPOs log matrix", "[wavefunction][template]", double, float) +TEMPLATE_LIST_TEST_CASE( + "RotatedSPOs log matrix", "[wavefunction][template]", TestTypeList) { using ValueType = typename SPOSetT::ValueType; using ValueMatrix = typename SPOSetT::ValueMatrix; @@ -525,8 +544,8 @@ TEMPLATE_TEST_CASE( // The log is multi-valued so this test may fail if the rotation parameters are // too large. The exponentials will be the same, though // exp(log(exp(A))) == exp(A) -TEMPLATE_TEST_CASE( - "RotatedSPOs exp-log matrix", "[wavefunction][template]", double, float) +TEMPLATE_LIST_TEST_CASE( + "RotatedSPOs exp-log matrix", "[wavefunction][template]", TestTypeList) { using ValueType = typename SPOSetT::ValueType; using ValueMatrix = typename SPOSetT::ValueMatrix; @@ -565,8 +584,8 @@ TEMPLATE_TEST_CASE( } } -TEMPLATE_TEST_CASE( - "RotatedSPOs hcpBe", "[wavefunction][template]", double, float) +TEMPLATE_LIST_TEST_CASE( + "RotatedSPOs hcpBe", "[wavefunction][template]", TestTypeList) { using RealType = typename OrbitalSetTraits::RealType; Communicate* c = OHMMS::Controller; @@ -687,8 +706,8 @@ TEMPLATE_TEST_CASE( } // Test construction of delta rotation -TEMPLATE_TEST_CASE("RotatedSPOs construct delta matrix", - "[wavefunction][template]", double, float) +TEMPLATE_LIST_TEST_CASE("RotatedSPOs construct delta matrix", + "[wavefunction][template]", TestTypeList) { using ValueType = typename SPOSetT::ValueType; using ValueMatrix = typename SPOSetT::ValueMatrix; @@ -751,8 +770,8 @@ TEMPLATE_TEST_CASE("RotatedSPOs construct delta matrix", } // Test using global rotation -TEMPLATE_TEST_CASE("RotatedSPOs read and write parameters", - "[wavefunction][template]", double, float) +TEMPLATE_LIST_TEST_CASE("RotatedSPOs read and write parameters", + "[wavefunction][template]", TestTypeList) { auto fake_spo = std::make_unique>(); fake_spo->setOrbitalSetSize(4); @@ -804,8 +823,8 @@ TEMPLATE_TEST_CASE("RotatedSPOs read and write parameters", } // Test using history list. -TEMPLATE_TEST_CASE("RotatedSPOs read and write parameters history", - "[wavefunction][template]", double, float) +TEMPLATE_LIST_TEST_CASE("RotatedSPOs read and write parameters history", + "[wavefunction][template]", TestTypeList) { auto fake_spo = std::make_unique>(); fake_spo->setOrbitalSetSize(4); @@ -920,8 +939,8 @@ class DummySPOSetWithMWT : public DummySPOSetWithoutMWT } }; -TEMPLATE_TEST_CASE( - "RotatedSPOs mw_ APIs", "[wavefunction][template]", double, float) +TEMPLATE_LIST_TEST_CASE( + "RotatedSPOs mw_ APIs", "[wavefunction][template]", TestTypeList) { // checking that mw_ API works in RotatedSPOs and is not defaulting to // SPOSet default implementation