diff --git a/src/Estimators/CMakeLists.txt b/src/Estimators/CMakeLists.txt index 66efbc0be6..89ffef4289 100644 --- a/src/Estimators/CMakeLists.txt +++ b/src/Estimators/CMakeLists.txt @@ -40,7 +40,8 @@ set(QMCEST_SRC PerParticleHamiltonianLogger.cpp ReferencePointsInput.cpp NEReferencePoints.cpp - NESpaceGrid.cpp) + NESpaceGrid.cpp + EnergyDensityEstimator.cpp) #################################### # create libqmcestimators diff --git a/src/Estimators/EnergyDensityEstimator.cpp b/src/Estimators/EnergyDensityEstimator.cpp new file mode 100644 index 0000000000..3b5209795d --- /dev/null +++ b/src/Estimators/EnergyDensityEstimator.cpp @@ -0,0 +1,420 @@ +///////////////////////////////////////////////////////////////////////////////////// +// 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 +// Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory +// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory +// Peter W. Doak, doakpw@ornl.gov, Oak Ridge National Laboratory +// +// Some code refactored from: QMCHamiltonian/EnergyDensityEstimator.cpp +////////////////////////////////////////////////////////////////////////////////////// + +#include "EnergyDensityEstimator.h" + +namespace qmcplusplus +{ + +struct PosCharge +{ + NEEnergyDensityEstimator::ParticlePos r_ptcls; + std::vector z_ptcls; +}; + +auto NEEnergyDensityEstimator::extractIonPositionsAndCharge(const ParticleSet& pset) +{ + const SpeciesSet& species(pset.getSpeciesSet()); + int charge_index = species.findAttribute("charge"); + int nspecies = species.TotalNum; + int nps = pset.getTotalNum(); + std::vector z_spec; + z_spec.resize(nspecies); + std::vector z_ptcls; + z_ptcls.resize(nps); + for (int spec = 0; spec < nspecies; spec++) + z_spec[spec] = species(charge_index, spec); + for (int i = 0; i < nps; i++) + z_ptcls[i] = z_spec[pset.GroupID[i]]; + ParticlePos r_ptcls; + r_ptcls.resize(pset.R.size()); + for (int i = 0; i < pset.R.size(); i++) + r_ptcls[i] = pset.R[i]; + if (pset.getLattice().SuperCellEnum != SUPERCELL_OPEN) + pset.applyMinimumImage(r_ptcls); + return PosCharge{r_ptcls, z_ptcls}; +} + +NEEnergyDensityEstimator::NEEnergyDensityEstimator(const EnergyDensityInput& input, + const PSPool& pset_pool, + DataLocality data_locality) + : OperatorEstBase(data_locality), input_(input), pset_dynamic_(getParticleSet(pset_pool, input.get_dynamic())) +{ + requires_listener_ = true; + my_name_ = "NEEnergyDensityEstimator"; + n_particles_ = pset_dynamic_.getTotalNum(); + if (!(input_.get_static().empty())) + { + pset_static_.emplace(getParticleSet(pset_pool, input.get_static())); + if (!input.get_ion_points()) + { + n_ions_ = pset_static_->getTotalNum(); + n_particles_ += n_ions_; + ed_ion_values_.resize(n_ions_, N_EDVALS); + r_ion_work_.resize(n_ions_, OHMMS_DIM); + } + } + + constructToReferencePoints(pset_dynamic_, pset_static_); + + bool periodic = pset_dynamic_.getLattice().SuperCellEnum != SUPERCELL_OPEN; + + spacegrid_inputs_ = input.get_space_grid_inputs(); + spacegrids_.reserve(spacegrid_inputs_.size()); + for (int ig = 0; ig < spacegrid_inputs_.size(); ++ig) + { + if (pset_static_) + { + auto [r_ptcls, z_ptcls] = extractIonPositionsAndCharge(*pset_static_); + spacegrids_.emplace_back(std::make_unique>(spacegrid_inputs_[ig], ref_points_->get_points(), r_ptcls, + z_ptcls, pset_dynamic_.getTotalNum(), N_EDVALS, periodic)); + } + else + spacegrids_.emplace_back( + std::make_unique>(spacegrid_inputs_[ig], ref_points_->get_points(), N_EDVALS, periodic)); + } +#ifndef NDEBUG + std::cout << "Instantiated " << spacegrids_.size() << " spacegrids\n"; +#endif +} + +NEEnergyDensityEstimator::NEEnergyDensityEstimator(const NEEnergyDensityEstimator& ede, const DataLocality dl) + : OperatorEstBase(dl), + input_(ede.input_), + pset_dynamic_(ede.pset_dynamic_), + pset_static_(ede.pset_static_), + n_particles_(ede.n_particles_), + n_ions_(ede.n_ions_), + spacegrid_inputs_(ede.spacegrid_inputs_) +{ + requires_listener_ = true; + data_locality_ = dl; + + constructToReferencePoints(pset_dynamic_, pset_static_); + + for (const auto& space_grid : ede.spacegrids_) + spacegrids_.emplace_back(std::make_unique>(*space_grid)); +} + +void NEEnergyDensityEstimator::constructToReferencePoints(ParticleSet& pset_dynamic, + const std::optional& pset_static) +{ + // Bringing a bunch of adhoc setup from legacy. + // removed redundant turnOnPerParticleSK this is handled via the EstimatorManagerNew if listeners are detected through + // CoulombPBCAA{AB} which are the actual operators that need it. + // Although it is possible that our copies of the particle sets will need per particle structure factors I don't think + // they do. + + RefVector pset_refs; + + if (pset_static_) + { + dtable_index_ = pset_dynamic_.addTable(pset_static_.value()); + pset_refs.push_back(pset_static_.value()); + } + + r_work_.resize(n_particles_); + ed_values_.resize(n_particles_, N_EDVALS); + + // right now this is only the case when pset_static_ && input.get_ion_points_ are true + if (n_ions_ > 0) + { + ed_ion_values_.resize(n_ions_, N_EDVALS); + r_ion_work_.resize(n_ions_, OHMMS_DIM); + for (int i = 0; i < n_ions_; ++i) + for (int d = 0; d < OHMMS_DIM; ++d) + r_ion_work_(i, d) = pset_static_->R[i][d]; + } + particles_outside_.resize(n_particles_, true); + ref_points_ = std::make_unique(input_.get_ref_points_input(), pset_dynamic_, pset_refs); +} + +NEEnergyDensityEstimator::~NEEnergyDensityEstimator(){}; + +void NEEnergyDensityEstimator::registerListeners(QMCHamiltonian& ham_leader) +{ + ListenerVector kinetic_listener("kinetic", getListener(kinetic_values_)); + QMCHamiltonian::mw_registerKineticListener(ham_leader, kinetic_listener); + ListenerVector potential_listener("potential", getListener(local_pot_values_)); + QMCHamiltonian::mw_registerLocalPotentialListener(ham_leader, potential_listener); + ListenerVector ion_potential_listener("potential", getListener(local_ion_pot_values_)); + QMCHamiltonian::mw_registerLocalIonPotentialListener(ham_leader, ion_potential_listener); +} + +/** This function collects the per particle energies. + * right now these are indentified by a string for each type. This could be optimized but + * could also be an insiginificant cost versus the frequently large number of values handled. + * The values themselves are a vector of size particle_num. + */ +ListenerVector::ReportingFunction NEEnergyDensityEstimator::getListener( + CrowdEnergyValues& local_values) +{ + return [&local_values](const int walker_index, const std::string& name, const Vector& inputV) { + if (walker_index >= local_values[name].size()) + local_values[name].resize(walker_index + 1); + local_values[name][walker_index] = inputV; + }; +} + +const ParticleSet& NEEnergyDensityEstimator::getParticleSet(const PSPool& psetpool, const std::string& psname) const +{ + auto pset_iter(psetpool.find(psname)); + if (pset_iter == psetpool.end()) + { + throw UniformCommunicateError("Particle set pool does not contain \"" + psname + + "\" so NEEnergyDensityEstimator::get_particleset fails!"); + + } + return *(pset_iter->second.get()); +} + +std::size_t NEEnergyDensityEstimator::getFullDataSize() const +{ + std::size_t size = OperatorEstBase::get_data().size(); + for (const UPtr>& grid : spacegrids_) + size += grid->getDataVector().size(); + size += 1; // for nsamples; + return size; +} + +void NEEnergyDensityEstimator::packData(PooledData& buffer) +{ + OperatorEstBase::packData(buffer); + for(auto& grid: spacegrids_) + buffer.add(grid->getDataVector().begin(),grid->getDataVector().end()); + buffer.add(nsamples_); +} + +void NEEnergyDensityEstimator::unpackData(PooledData& buffer) +{ + OperatorEstBase::unpackData(buffer); + for(auto& grid: spacegrids_) + buffer.get(grid->getDataVector().begin(),grid->getDataVector().end()); + buffer.get(nsamples_); +} + +// void NEEnergyDensityEstimator::unpackData(std::vector& operator_receive_buffer) +// { + +// } + +void NEEnergyDensityEstimator::accumulate(const RefVector& walkers, + const RefVector& psets, + const RefVector& wfns, + const RefVector& hams, + RandomBase& rng) +{ + // Variable population is possible during DMC. + reduced_local_kinetic_values_.resize(walkers.size()); + reduced_local_pot_values_.resize(walkers.size()); + reduced_local_ion_pot_values_.resize(walkers.size()); + + // Depending on Hamiltonian setup one or more of these values could be absent. + if (!local_pot_values_.empty()) + combinePerParticleEnergies(local_pot_values_, reduced_local_pot_values_); + if (!kinetic_values_.empty()) + combinePerParticleEnergies(kinetic_values_, reduced_local_kinetic_values_); + if (pset_static_) + { + if (!local_ion_pot_values_.empty()) + combinePerParticleEnergies(local_ion_pot_values_, reduced_local_ion_pot_values_); + } + for (int iw = 0; iw < walkers.size(); ++iw) + { + walkers_weight_ += walkers[iw].get().Weight; + evaluate(psets[iw], walkers[iw], iw); + } +} + +void NEEnergyDensityEstimator::evaluate(ParticleSet& pset, const MCPWalker& walker, const int walker_index) +{ + //Collect positions from ParticleSets + int p_count = 0; + { + const ParticlePos& Rs = pset.R; + for (int i = 0; i < Rs.size(); i++) + { + r_work_[p_count] = Rs[i]; + p_count++; + } + } + if (pset_static_ && !input_.get_ion_points()) + { + const ParticlePos& Rs = pset_static_->R; + for (int i = 0; i < Rs.size(); i++) + { + r_work_[p_count] = Rs[i]; + p_count++; + } + } + + if (pset.getLattice().SuperCellEnum != SUPERCELL_OPEN) + pset.applyMinimumImage(r_work_); + //Convert information accumulated in ParticleSets into EnergyDensity quantities + Real weight = walker.Weight; + p_count = 0; + { + const auto& Ts = reduced_local_kinetic_values_[walker_index]; + const auto& Vd = reduced_local_pot_values_[walker_index]; + for (int i = 0; i < pset.getTotalNum(); i++) + { + ed_values_(p_count, W) = weight; + ed_values_(p_count, T) = weight * Ts[i]; + ed_values_(p_count, V) = weight * Vd[i]; + p_count++; + } + } + if (pset_static_) + { + const ParticleSet& Ps = *pset_static_; + const auto& Vs = reduced_local_ion_pot_values_[walker_index]; + if (!input_.get_ion_points()) + for (int i = 0; i < Ps.getTotalNum(); i++) + { + ed_values_(p_count, W) = weight; + ed_values_(p_count, T) = 0.0; + if (Vs.size() > i) + ed_values_(p_count, V) = weight * Vs[i]; + p_count++; + } + else + for (int i = 0; i < Ps.getTotalNum(); i++) + { + ed_ion_values_(i, W) = weight; + ed_ion_values_(i, T) = 0.0; + ed_ion_values_(i, V) = weight * Vs[i]; + } + } + //Accumulate energy density in spacegrids + //const auto& dtab(pset.getDistTableAB(dtable_index_)); + fill(particles_outside_.begin(), particles_outside_.end(), true); + for (int i = 0; i < spacegrids_.size(); i++) + { + NESpaceGrid& sg = *(spacegrids_[i]); + sg.accumulate(r_work_, ed_values_, particles_outside_); //, dtab); + } + + //Accumulate energy density of particles outside any spacegrid + int bi, v; + const int bimax = outside_buffer_offset + N_EDVALS; + for (int p = 0; p < particles_outside_.size(); p++) + { + if (particles_outside_[p]) + { + for (bi = outside_buffer_offset, v = 0; bi < bimax; bi++, v++) + { + data_[bi] += ed_values_(p, v); + } + } + } + if (input_.get_ion_points()) + { + // Accumulate energy density for ions at a point field + bi = outside_buffer_offset + N_EDVALS; + for (int i = 0; i < n_ions_; i++) + for (v = 0; v < N_EDVALS; v++, bi++) + { + data_[bi] += ed_ion_values_(i, v); + } + } + nsamples_++; +} + +void NEEnergyDensityEstimator::collect(const RefVector& type_erased_operator_estimators) +{ + int num_crowds = type_erased_operator_estimators.size(); + for (int ig = 0; ig < spacegrids_.size(); ++ig) + { + RefVector> crowd_grids; + crowd_grids.reserve(num_crowds); + for (OperatorEstBase& crowd_oeb : type_erased_operator_estimators) + { + NEEnergyDensityEstimator& crowd_ede = dynamic_cast(crowd_oeb); + NESpaceGrid& grid_ref = *(crowd_ede.spacegrids_[ig]); + crowd_grids.push_back(grid_ref); + //nsamples += grid_ref.samples; + } + NESpaceGrid::collect(*(spacegrids_[ig]), crowd_grids); + + } + OperatorEstBase::collect(type_erased_operator_estimators); +} + +void NEEnergyDensityEstimator::registerOperatorEstimator(hdf_archive& file) +{ + hdf_path hdf_name{my_name_}; + + hdf_path path_variables = hdf_name / std::string_view("variables"); + file.push(path_variables, true); + file.write(const_cast(n_particles_), "nparticles"); + auto nspacegrids = input_.get_space_grid_inputs().size(); + file.write(nspacegrids, "nspacegrids"); + file.write(const_cast(nsamples_), "nsamples"); + + if (input_.get_ion_points()) + { + file.write(const_cast(n_ions_), "nions"); + file.write(const_cast&>(r_ion_work_), "ion_positions"); + } + file.pop(); + + file.push(hdf_name); + + ref_points_->write(file); + + h5desc_.emplace_back(hdf_name / "outside"); + auto& ohOutside = h5desc_.back(); + std::vector ng(1); + ng[0] = N_EDVALS; + ohOutside.set_dimensions(ng, outside_buffer_offset); + for (int i = 0; i < spacegrids_.size(); i++) + spacegrids_[i]->registerGrid(file, i); + + file.pop(); +} + +void NEEnergyDensityEstimator::write(hdf_archive& file) +{ + hdf_path hdf_name{my_name_}; + file.push(hdf_name); + file.write(particles_outside_, "outside"); + for (int i = 0; i < spacegrids_.size(); i++) + spacegrids_[i]->write(file); + if (input_.get_ion_points()) + file.write(ed_ion_values_, "ions"); + file.pop(); +} + +std::unique_ptr NEEnergyDensityEstimator::spawnCrowdClone() const +{ + auto spawn_data_locality = data_locality_; + auto data_size = this->data_.size(); + UPtr spawn(std::make_unique(*this, spawn_data_locality)); + spawn->get_data().resize(data_size); + return spawn; +} + +RefVector> NEEnergyDensityEstimator::getSpaceGrids() { return convertUPtrToRefVector(spacegrids_); } + +void NEEnergyDensityEstimator::startBlock(int steps) {} + +// RefVector>& NEEnergyDensityEstimator::getExtraData() { +// RefVector> refs; +// for(auto& space_grid : space_grids_) { +// refs.push_back(getDataVector()); +// } +// return refs; +// } + +} // namespace qmcplusplus diff --git a/src/Estimators/EnergyDensityEstimator.h b/src/Estimators/EnergyDensityEstimator.h new file mode 100644 index 0000000000..9f309071d2 --- /dev/null +++ b/src/Estimators/EnergyDensityEstimator.h @@ -0,0 +1,202 @@ +////////////////////////////////////////////////////////////////////////////////////// +// 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: Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory +// Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory +// Peter W. Doak, doakpw@ornl.gov, Oak Ridge National Laboratory +// +// Some code refactored from: QMCHamiltonian/EnergyDensityEstimator.h +////////////////////////////////////////////////////////////////////////////////////// + +#ifndef QMCPLUSPLUS_ESTIMATOR_ENERGY_DENSITY_H +#define QMCPLUSPLUS_ESTIMATOR_ENERGY_DENSITY_H + +#include "OperatorEstBase.h" +#include "OhmmsPETE/OhmmsMatrix.h" +#include "NEReferencePoints.h" +#include "NESpaceGrid.h" +#include "EnergyDensityInput.h" +#include + +namespace qmcplusplus +{ +class NEEnergyDensityEstimator : public OperatorEstBase +{ +public: + using Real = QMCTraits::RealType; + using FullPrecReal = QMCTraits::FullPrecRealType; + using Point = typename NEReferencePoints::Point; + using Points = typename NEReferencePoints::Points; + using POLT = PtclOnLatticeTraits; + using ParticlePos = POLT::ParticlePos; + using PSPool = typename ParticleSetPool::PoolType; + + /** Constructor + * + * \param[in] input immutable parameters from user input + * \param[in] pset_pool global particle set pool object to allow arbitrarily named static pset to be used. (Necessary?) + * \param[in] dl datalocality tag for alternate memory management implementations. Currently EDE only implements + * crowd locality which duplicates the full estimator data object per crowd + 1 at ranks scope. + */ + NEEnergyDensityEstimator(const EnergyDensityInput& input, const PSPool& PSP, DataLocality dl = DataLocality::crowd); + NEEnergyDensityEstimator(const NEEnergyDensityEstimator& ede, const DataLocality dl); + +private: + /** shared construction code + * stateful setup from legacy. + * removed redundant turnOnPerParticleSK this is handled via the EstimatorManagerNew if listeners are detected through + * CoulombPBCAA{AB} which are the actual operators that need it. + * Although it is possible that our copies of the particle sets will need per particle structure factors I don't think + * they do. I think it is needed for the actual particle sets involved in walking. + */ + void constructToReferencePoints(ParticleSet& pset_dynamic, const std::optional& pset_static); + +public: + /** @ingroup OperatorEstBase Overrides + * should NESpaceGrid have anything to do with OperatorEstBase API? + * @{ + */ + ~NEEnergyDensityEstimator() override; + + /** Register listeners for energy values + */ + void registerListeners(QMCHamiltonian& ham_leader) override; + + ListenerVector::ReportingFunction getListener(); + + /** accumulate 1 or more walkers of EnergyDensity samples + */ + void accumulate(const RefVector& walkers, + const RefVector& psets, + const RefVector& wfns, + const RefVector& hams, + RandomBase& rng) override; + + void evaluate(ParticleSet& pset, const MCPWalker& walker, const int walker_index); + + /** this allows the EstimatorManagerNew to reduce without needing to know the details + * of SpaceGrid's data. + * + * can use base class default until crowd level MomentumDistribution + * estimators don't have a copy of the density grid. + */ + void collect(const RefVector& type_erased_operator_estimators) override; + + std::size_t getFullDataSize() const override; + + void packData(PooledData& buffer) override; + void unpackData(PooledData& buffer) override; + + UPtr spawnCrowdClone() const override; + + /** start block entry point + */ + void startBlock(int steps) override; + + /// @} + + /** return lambda function to register as listener + * the purpose of this function is to factor out the production of the lambda for unit testing + * \param[out] values + */ + ListenerVector::ReportingFunction getListener(CrowdEnergyValues& values); + + const ParticleSet& getParticleSet(const PSPool& psetpool, const std::string& psname) const; + + void registerOperatorEstimator(hdf_archive& file) override; + + void write(hdf_archive& file); + + RefVector> getSpaceGrids(); + + RefVector>& getExtraData(); +private: + auto extractIonPositionsAndCharge(const ParticleSet& pset); + + EnergyDensityInput input_; + + CrowdEnergyValues kinetic_values_; + CrowdEnergyValues local_pot_values_; + CrowdEnergyValues local_ion_pot_values_; + + /// working space for reduced_values_; + std::vector> reduced_local_pot_values_; + std::vector> reduced_local_kinetic_values_; + std::vector> reduced_local_ion_pot_values_; + + ParticleSet pset_dynamic_; + std::optional pset_static_; + int dtable_index_ = -1; + + int n_particles_; + int n_ions_; + + UPtr ref_points_; + + /// unboxed SpaceGridInputs for child and clone objects to refrence + std::vector spacegrid_inputs_; + + /** EnergyDenstity quantities + * legacy style enum into vector not ideal. + */ + enum + { + W = 0, + T, + V, + N_EDVALS + }; + + /** @ingroup outside particles + * for the EnergyDensity of particles falling outside any spacegrid + * @{ */ + /** in legacy the starting index into the collectibles buffer. + * Maintained to use more legacy code without modificaiton in the short term. + * In the long term its possible the entire way the grid data is structured in memory should be redesigned. + * + * Right now we put the ed_values and ed_ion_values_ into the base class data_. + */ + int outside_buffer_offset{0}; + std::vector particles_outside_; + + /// @} + + /// spacegrids are used to find which cell domain contains the Energy information of particles + UPtrVector> spacegrids_; + + /** @ingroup working variables reused for each walker + * Do not use for persistent state. + * @{ */ + + /// particle positions includes ion positions if pset_stat_ and !input_.get_ion_points() + ParticlePos r_work_; + /// ion positions if pset_static_ && input.get_ion_points_ are true + Matrix r_ion_work_; + + /** preallocated value matrices + * It seems like these should be local variables in evaluate, + * but this is how it was done in legacy so just ported for now. + * depending on these as state doesn't seem to be done now and shouldn't be. + */ + Matrix ed_values_; + /// if input_.get_ion_points and pset_static_ then ion values end up here + Matrix ed_ion_values_; + /// @} + + //number of samples accumulated + int nsamples_; + + DataLocality data_locality_{DataLocality::crowd}; + + // //needed (temporarily) for chempot + // //ParticleSet should carry Zptcl so it doesn't have + // // to be computed everywhere from species + // std::vector Zptcl; + // ParticlePos Rptcl; +}; +} // namespace qmcplusplus + +#endif diff --git a/src/Estimators/OperatorEstBase.cpp b/src/Estimators/OperatorEstBase.cpp index 5b4674655f..dcc96e8bee 100644 --- a/src/Estimators/OperatorEstBase.cpp +++ b/src/Estimators/OperatorEstBase.cpp @@ -61,6 +61,14 @@ void OperatorEstBase::write(hdf_archive& file) file.pop(); } +void OperatorEstBase::packData(PooledData& buffer) { + buffer.add(data_.begin(),data_.end()); +} + +void OperatorEstBase::unpackData(PooledData& buffer) { + buffer.get(data_.begin(), data_.end()); +} + void OperatorEstBase::zero() { if (data_locality_ == DataLocality::rank || data_locality_ == DataLocality::crowd) diff --git a/src/Estimators/OperatorEstBase.h b/src/Estimators/OperatorEstBase.h index 7198cfbd5c..2d924dd0ab 100644 --- a/src/Estimators/OperatorEstBase.h +++ b/src/Estimators/OperatorEstBase.h @@ -91,9 +91,33 @@ class OperatorEstBase virtual void startBlock(int steps) = 0; + const std::vector& get_data() const { return data_; } std::vector& get_data() { return data_; } + virtual std::size_t getFullDataSize() const { return data_.size(); } + /** @ingroup Functions to add or remove estimator data from PooledData + * @brief used for MPI reduction. + * These are only used on the rank estimator owned by EstimatorManagerNew. + * The rank EstimatorManagerNew owns the buffer. + * It is not intended to store the state of the estimator. + * The packing and unpacking functions must follow the same sequence of adds or gets + * as PooledData is a stateful sequence of bytes with an internal position cursor. + * @{ + */ + + /** Packs data from native container types in a subtype of Operator est base + * to buffer of type Real for reduction over MPI. + * I.e. writes to pooled data. + */ + virtual void packData(PooledData& buffer); + /** Unpacks data from mpi buffer of type Real into native container types + * after a reduction over MPI. + * i.e. reads from pooled data. + */ + virtual void unpackData(PooledData& buffer); + ///@} + /*** create and tie OperatorEstimator's observable_helper hdf5 wrapper to stat.h5 file * @param gid hdf5 group to which the observables belong * diff --git a/src/Estimators/PerParticleHamiltonianLogger.cpp b/src/Estimators/PerParticleHamiltonianLogger.cpp index 6ca92e82f0..155e21fa61 100644 --- a/src/Estimators/PerParticleHamiltonianLogger.cpp +++ b/src/Estimators/PerParticleHamiltonianLogger.cpp @@ -17,7 +17,7 @@ namespace qmcplusplus { PerParticleHamiltonianLogger::PerParticleHamiltonianLogger(PerParticleHamiltonianLoggerInput&& input, int rank) - : OperatorEstBase(DataLocality::crowd), rank_estimator_(nullptr), input_(input), rank_(rank) + : OperatorEstBase(DataLocality::crowd), input_(input), rank_(rank) { requires_listener_ = true; my_name_ = "PerParticleHamiltonianLogger"; @@ -26,8 +26,8 @@ PerParticleHamiltonianLogger::PerParticleHamiltonianLogger(PerParticleHamiltonia rank_fstream_.open(filename, std::ios::out); } -PerParticleHamiltonianLogger::PerParticleHamiltonianLogger(const PerParticleHamiltonianLogger& pphl, DataLocality dl) - : OperatorEstBase(dl), rank_estimator_(const_cast(&pphl)), input_(pphl.input_) +PerParticleHamiltonianLogger::PerParticleHamiltonianLogger(PerParticleHamiltonianLogger& pphl, DataLocality dl) + : OperatorEstBase(dl), rank_estimator_(makeOptionalRef(pphl)), input_(pphl.input_) { requires_listener_ = true; my_name_ = pphl.name_; @@ -69,18 +69,27 @@ void PerParticleHamiltonianLogger::accumulate(const RefVector& walker walker_ids_.clear(); for (MCPWalker& walker : walkers) walker_ids_.push_back(walker.getWalkerID()); - rank_estimator_->write(values_, walker_ids_); + rank_estimator_->get().write(values_, walker_ids_); // \todo some per crowd reduction. // clear log values } +PerParticleHamiltonianLogger::Real PerParticleHamiltonianLogger::sumOverAll() const +{ + Real sum{0}; + for (auto& [component, values] : values_) + for (auto& a_vector : values) + sum += std::accumulate(a_vector.begin(), a_vector.end(), 0.0); + return sum; +} + std::unique_ptr PerParticleHamiltonianLogger::spawnCrowdClone() const { std::size_t data_size = data_.size(); auto spawn_data_locality = data_locality_; - auto spawn = std::make_unique(*this, spawn_data_locality); + auto spawn = std::make_unique(const_cast(*this), spawn_data_locality); spawn->get_data().resize(data_size, 0.0); return spawn; } @@ -104,6 +113,7 @@ void PerParticleHamiltonianLogger::registerListeners(QMCHamiltonian& ham_leader) { ListenerVector listener(name_, getLogger()); QMCHamiltonian::mw_registerLocalEnergyListener(ham_leader, listener); + QMCHamiltonian::mw_registerLocalIonPotentialListener(ham_leader, listener); } void PerParticleHamiltonianLogger::startBlock(int steps) diff --git a/src/Estimators/PerParticleHamiltonianLogger.h b/src/Estimators/PerParticleHamiltonianLogger.h index 54bf76132a..43f435f71a 100644 --- a/src/Estimators/PerParticleHamiltonianLogger.h +++ b/src/Estimators/PerParticleHamiltonianLogger.h @@ -21,19 +21,18 @@ #include #include "OhmmsPETE/OhmmsVector.h" #include "QMCHamiltonians/Listener.hpp" - +#include "type_traits/OptionalRef.hpp" namespace qmcplusplus { class PerParticleHamiltonianLogger : public OperatorEstBase { public: - using Real = QMCTraits::RealType; + using Real = QMCTraits::RealType; using CrowdLogValues = std::unordered_map>>; - PerParticleHamiltonianLogger(PerParticleHamiltonianLoggerInput&& input, int rank); - PerParticleHamiltonianLogger(const PerParticleHamiltonianLogger& other, DataLocality data_locality); + PerParticleHamiltonianLogger(PerParticleHamiltonianLogger& other, DataLocality data_locality); void accumulate(const RefVector& walkers, const RefVector& psets, @@ -55,20 +54,28 @@ class PerParticleHamiltonianLogger : public OperatorEstBase void write(CrowdLogValues& values, const std::vector& walkers_ids); + /** This function is supplied for testing it sums over all values currently stored by the + * per particle logger. This is useful in unit testing some of the other estimators + * that make use of Listeners. + */ + Real sumOverAll() const; + int get_block() { return block_; } private: bool crowd_clone = false; - PerParticleHamiltonianLogger * const rank_estimator_; + const OptionalRef rank_estimator_; PerParticleHamiltonianLoggerInput input_; int rank_; CrowdLogValues values_; std::vector walker_ids_; const std::string name_{"PerParticleHamiltonianLogger"}; + /// rank owned fstream std::fstream rank_fstream_; + /// use it to prevent race when writing per-crowd. std::mutex write_lock; int block_ = 0; }; - -} + +} // namespace qmcplusplus #endif diff --git a/src/Estimators/tests/CMakeLists.txt b/src/Estimators/tests/CMakeLists.txt index 8f51c78f00..0f9857dc8c 100644 --- a/src/Estimators/tests/CMakeLists.txt +++ b/src/Estimators/tests/CMakeLists.txt @@ -54,7 +54,10 @@ set(SRCS test_EstimatorManagerCrowd.cpp test_SpaceGridInput.cpp test_NESpaceGrid.cpp + EnergyDensityTest.cpp + MockGoldWalkerElements.cpp test_EnergyDensityInput.cpp + test_EnergyDensityEstimator.cpp ) add_executable(${UTEST_EXE} ${SRCS}) diff --git a/src/Estimators/tests/EnergyDensityTest.cpp b/src/Estimators/tests/EnergyDensityTest.cpp new file mode 100644 index 0000000000..e0f1f1a9e5 --- /dev/null +++ b/src/Estimators/tests/EnergyDensityTest.cpp @@ -0,0 +1,96 @@ +////////////////////////////////////////////////////////////////////////////////////// +// This file is distributed under the University of Illinois/NCSA Open Source License. +// See LICENSE file in top directory for details. +// +// Copyright (c) 2024 QMCPACK developers. +// +// File developed by: Peter Doak, doakpw@ornl.gov, Oak Ridge National Laboratory +// +// File created by: Peter Doak, doakpw@ornl.gov, Oak Ridge National Laboratory +////////////////////////////////////////////////////////////////////////////////////// + +#include "EnergyDensityTest.h" +#include "Utilities/ProjectData.h" +#include "Utilities/ResourceCollection.h" +#include "Utilities/StdRandom.h" +#include "GenerateRandomParticleSets.h" +#include "ValidEnergyDensityInput.h" + +namespace qmcplusplus +{ +namespace testing +{ + +EnergyDensityTest::EnergyDensityTest(Communicate* comm, int num_walkers, bool generate_test_data) : test_project_("test", ProjectData::DriverVersion::BATCH), gold_elem_(testing::makeGoldWalkerElementsWithEE(comm, test_project_.getRuntimeOptions())), walkers_(num_walkers, MCPWalker(gold_elem_.pset_elec.getTotalNum())) +{ + // This has to be first because before the golden Hamiltonian is copied + // it needs to know if there are per particle listeners. + gold_elem_.ham.informOperatorsOfListener(); + + psets_ = testing::generateRandomParticleSets(gold_elem_.pset_elec, gold_elem_.pset_ions, deterministic_rs_, num_walkers, + generate_test_data); + + auto pset_refs(makeRefVector(psets_)); + auto& trial_wavefunction = gold_elem_.twf; + + Libxml2Document doc; + using Input = testing::EnergyDensityInputs; + Input input; + bool okay = doc.parseFromString(Input::getXml(Input::valid::CELL)); + xmlNodePtr node = doc.getRoot(); + + UPtr edein; + edein = std::make_unique(node); + NEEnergyDensityEstimator e_den_est(*edein, gold_elem_.particle_pool.getPool()); + + for (int iw = 0; iw < num_walkers; ++iw) + { + twfs_.emplace_back(gold_elem_.twf.makeClone(psets_[iw])); + hams_.emplace_back(gold_elem_.ham.makeClone(psets_[iw], *twfs_.back())); + } + + auto walker_refs = makeRefVector(walkers_); + + RefVector ham_refs = convertUPtrToRefVector(hams_); + + hams_[0]->createResource(ham_res_); + psets_[0].createResource(pset_res_); + twfs_[0]->createResource(wfc_res_); +} + +NEEnergyDensityEstimator& EnergyDensityTest::getEnergyDensityEstimator() +{ + if (!eden_est_) + { + Libxml2Document doc; + using Input = testing::EnergyDensityInputs; + bool okay = doc.parseFromString(Input::getXml(Input::valid::CELL)); + xmlNodePtr node = doc.getRoot(); + edein_ = std::make_unique(node); + eden_est_ = std::make_unique(*edein_, gold_elem_.particle_pool.getPool()); + } + return *eden_est_; +} + +RefVectorWithLeader EnergyDensityTest::getWalkerList() { + return {walkers_[0], makeRefVector(walkers_)}; } +RefVectorWithLeader EnergyDensityTest::getPSetList() { + return {psets_[0], makeRefVector(psets_)}; } +RefVectorWithLeader EnergyDensityTest::getHamList() { + return {*hams_[0], convertUPtrToRefVector(hams_)}; } +RefVectorWithLeader EnergyDensityTest::getTwfList() { + return {*twfs_[0], convertUPtrToRefVector(twfs_)}; } + +ResourceCollection& EnergyDensityTest::getPSetRes() { + return pset_res_; +} +ResourceCollection& EnergyDensityTest::getHamRes() { + return ham_res_; +} +ResourceCollection& EnergyDensityTest::getTwfRes() { + return wfc_res_; +} + +} // namespace testing +} // namespace qmcplusplus + diff --git a/src/Estimators/tests/EnergyDensityTest.h b/src/Estimators/tests/EnergyDensityTest.h new file mode 100644 index 0000000000..11dfe5c23a --- /dev/null +++ b/src/Estimators/tests/EnergyDensityTest.h @@ -0,0 +1,114 @@ +////////////////////////////////////////////////////////////////////////////////////// +// This file is distributed under the University of Illinois/NCSA Open Source License. +// See LICENSE file in top directory for details. +// +// Copyright (c) 2024 QMCPACK developers. +// +// File developed by: Peter Doak, doakpw@ornl.gov, Oak Ridge National Laboratory +// +// File created by: Peter Doak, doakpw@ornl.gov, Oak Ridge National Laboratory +////////////////////////////////////////////////////////////////////////////////////// + +#ifndef QMCPLUSPLUS_ENERGYDENSITYTEST_H +#define QMCPLUSPLUS_ENERGYDENSITYTEST_H + +#include "EnergyDensityEstimator.h" +#include "ResourceCollection.h" +#include "Utilities/ProjectData.h" +#include "MockGoldWalkerElements.h" + +namespace qmcplusplus +{ +namespace testing +{ + +/** This class is here to allow sharing of the code to stand up an evironment to test + * the energy density estimator. + */ +class EnergyDensityTest +{ +public: + using MCPWalker = typename OperatorEstBase::MCPWalker; + + EnergyDensityTest(Communicate* comm, int num_walkers, bool generate_test_data = false); + + RefVector getPsetRefs() { return makeRefVector(psets_); } + NEEnergyDensityEstimator& getEnergyDensityEstimator(); + + MockGoldWalkerElements& getGoldElements() { return gold_elem_; } + + RefVectorWithLeader getWalkerList(); + RefVectorWithLeader getPSetList(); + RefVectorWithLeader getHamList(); + RefVectorWithLeader getTwfList(); + + ResourceCollection& getPSetRes(); + ResourceCollection& getHamRes(); + ResourceCollection& getTwfRes(); + +private: + ProjectData test_project_; + MockGoldWalkerElements gold_elem_; + + UPtrVector hams_; + UPtrVector twfs_; + std::vector psets_; + std::vector walkers_; + + ResourceCollection pset_res_{"test_pset_res"}; + ResourceCollection ham_res_{"test_ham_res"}; + ResourceCollection wfc_res_{"test_wfc_res"}; + + UPtr eden_est_; + UPtr edein_; + + /// Canned test data so rng differences don't cause fails. + std::vector deterministic_rs_ = { + { + {0.515677886, 0.9486072745, -1.17423246}, + {-0.3166678423, 1.376550506, 1.445290031}, + {1.96071365, 2.47265689, 1.051449486}, + {0.745853269, 0.5551359072, 4.080774681}, + {-0.3515016103, -0.5192222523, 0.9941510909}, + {-0.8354426872, 0.7071638258, -0.3409843552}, + {0.4386044751, 1.237378731, 2.331874152}, + {2.125850717, 0.3221067321, 0.5825731561}, + }, + { + {-0.4633736785, 0.06318772224, -0.8580153742}, + {-1.174926354, -0.6276503679, 0.07458759314}, + {1.327618206, 2.085829379, 1.415749862}, + {0.9114727103, 0.1789183931, -0.08135540251}, + {-2.267908723, 0.802928773, 0.9522812957}, + {1.502715257, -1.84493529, 0.2805620469}, + {3.168934617, 0.1348337978, 1.371092768}, + {0.8310229518, 1.070827168, 1.18016733}, + }, + { + {-0.04847732172, -1.201739871, -1.700527771}, + {0.1589259538, -0.3096047065, -2.066626415}, + {2.255976232, 1.629132391, -0.8024446773}, + {2.534792993, 3.121092901, 1.609780703}, + {-0.2892376071, -0.152022511, -2.727613712}, + {0.2477154804, 0.5039232765, 2.995702733}, + {3.679345099, 3.037770313, 2.808899306}, + {0.6418578532, 1.935944544, 1.515637954}, + }, + { + {0.91126951, 0.0234699242, 1.442297821}, + {-0.9240061217, -0.1014997844, 0.9081020061}, + {1.887794866, 2.210192703, 2.209118551}, + {2.758945014, -1.21140421, 1.3337907}, + {0.376540703, 0.3485486555, 0.9056881595}, + {-0.3512770187, -0.4056820917, -2.068499576}, + {0.5358460986, 2.720153363, 1.41999706}, + {2.284020089, 1.173071915, 1.044597715}, + }, + }; + +}; + +} // namespace testing +} // namespace qmcplusplus + +#endif diff --git a/src/Estimators/tests/GenerateRandomParticleSets.cpp b/src/Estimators/tests/GenerateRandomParticleSets.cpp index f74c26145b..bb5f85b497 100644 --- a/src/Estimators/tests/GenerateRandomParticleSets.cpp +++ b/src/Estimators/tests/GenerateRandomParticleSets.cpp @@ -22,49 +22,39 @@ namespace qmcplusplus { namespace testing { -template std::vector generateRandomParticleSets(ParticleSet& pset_target, ParticleSet& pset_source, std::vector& deterministic_rs, - int num_psets) + int num_psets, + bool generate_test_data) { int nwalkers = num_psets; std::vector psets(num_psets, pset_target); - if constexpr (GEN_TEST_DATA) + if (generate_test_data) { std::cout << "Initialize OneBodyDensityMatrices::accumulate psets with:\n{"; std::vector psets; for (int iw = 0; iw < nwalkers; ++iw) { - //psets.emplace_back(pset_target); + psets.emplace_back(pset_target); psets.back().randomizeFromSource(pset_source); std::cout << "{"; for (auto r : psets.back().R) std::cout << NativePrint(r) << ","; std::cout << "},\n"; + psets.back().update(); } std::cout << "}\n"; } else { - for (int iw = 0; iw < nwalkers; ++iw) + for (int iw = 0; iw < nwalkers; ++iw) { psets[iw].R = deterministic_rs[iw]; + psets[iw].update(); + } } return psets; } -template -std::vector generateRandomParticleSets(ParticleSet& pset_target, - ParticleSet& pset_source, - std::vector& deterministic_rs, - int num_psets); - -template -std::vector generateRandomParticleSets(ParticleSet& pset_target, - ParticleSet& pset_source, - std::vector& deterministic_rs, - int num_psets); - - } // namespace testing } // namespace qmcplusplus diff --git a/src/Estimators/tests/GenerateRandomParticleSets.h b/src/Estimators/tests/GenerateRandomParticleSets.h index a20ced1d73..f253d56217 100644 --- a/src/Estimators/tests/GenerateRandomParticleSets.h +++ b/src/Estimators/tests/GenerateRandomParticleSets.h @@ -28,22 +28,11 @@ namespace testing /** * This function sets particle set positions from a set of Rs or writes out a set of positions for test reproducibility. */ -template std::vector generateRandomParticleSets(ParticleSet& pset_target, ParticleSet& pset_source, std::vector& deterministic_rs, - int num_psets); - -extern template std::vector generateRandomParticleSets( - ParticleSet& pset_target, - ParticleSet& pset_source, - std::vector& deterministic_rs, - int num_psets); -extern template std::vector generateRandomParticleSets( - ParticleSet& pset_target, - ParticleSet& pset_source, - std::vector& deterministic_rs, - int num_psets); + int num_psets, + bool generate_test_data = false); } // namespace testing } // namespace qmcplusplus diff --git a/src/Estimators/tests/MockGoldWalkerElements.cpp b/src/Estimators/tests/MockGoldWalkerElements.cpp new file mode 100644 index 0000000000..e09aff07eb --- /dev/null +++ b/src/Estimators/tests/MockGoldWalkerElements.cpp @@ -0,0 +1,50 @@ +////////////////////////////////////////////////////////////////////////////////////// +// 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: Peter Doak, doakpw@ornl.gov, Oak Ridge National Lab +// +// File created by: Peter Doak, doakpw@ornl.gov, Oak Ridge National Lab +////////////////////////////////////////////////////////////////////////////////////// + +#include "MockGoldWalkerElements.h" + +namespace qmcplusplus +{ +namespace testing +{ +MockGoldWalkerElements::MockGoldWalkerElements(Communicate* comm, + RuntimeOptions& runtime_opt, + WaveFunctionPoolFactoryFunc wavefunction_pool_fac_func, + HamPoolFactoryFunc ham_pool_fac_func) + : particle_pool(MinimalParticlePool::make_diamondC_1x1x1(comm)), + pset_elec(*(particle_pool.getParticleSet("e"))), + pset_ions(*(particle_pool.getParticleSet("ion"))), + wavefunction_pool(wavefunction_pool_fac_func(runtime_opt, comm, particle_pool)), + hamiltonian_pool(ham_pool_fac_func(comm, particle_pool, wavefunction_pool)), + twf(*(wavefunction_pool.getPrimary())), + ham(*(hamiltonian_pool.getPrimary())) +{} + +MockGoldWalkerElements makeGoldWalkerElementsWithEE(Communicate* comm, RuntimeOptions runtime_opt) +{ + using namespace std::placeholders; + MockGoldWalkerElements::WaveFunctionPoolFactoryFunc wfp_diamondC = + std::bind(&MinimalWaveFunctionPool::make_diamondC_1x1x1, _1, _2, _3); + MockGoldWalkerElements::HamPoolFactoryFunc hamp_ee = std::bind(MinimalHamiltonianPool::make_hamWithEE, _1, _2, _3); + return MockGoldWalkerElements(comm, runtime_opt, wfp_diamondC, hamp_ee); +} // namespace qmcplusplus + +MockGoldWalkerElements makeGoldWalkerElementsWithEEEI(Communicate* comm, RuntimeOptions runtime_opt) +{ + using namespace std::placeholders; + MockGoldWalkerElements::WaveFunctionPoolFactoryFunc wfp_diamondC = + std::bind(&MinimalWaveFunctionPool::make_diamondC_1x1x1, _1, _2, _3); + MockGoldWalkerElements::HamPoolFactoryFunc hamp_ee = std::bind(MinimalHamiltonianPool::makeHamWithEEEI, _1, _2, _3); + return MockGoldWalkerElements(comm, runtime_opt, wfp_diamondC, hamp_ee); +} + +} // namespace testing +} // namespace qmcplusplus diff --git a/src/Estimators/tests/MockGoldWalkerElements.h b/src/Estimators/tests/MockGoldWalkerElements.h new file mode 100644 index 0000000000..cbff360cb1 --- /dev/null +++ b/src/Estimators/tests/MockGoldWalkerElements.h @@ -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: Peter Doak, doakpw@ornl.gov, Oak Ridge National Lab +// +// File created by: Peter Doak, doakpw@ornl.gov, Oak Ridge National Lab +////////////////////////////////////////////////////////////////////////////////////// + +/** \file + * Some ParticleSet functions use the global Random so we need some helper functions to + * avoid interminant test state when multiple tests are run from a single test program. + */ + +#ifndef QMCPLUSPLUS_MOCK_WALKER_ELEMENTS_FOR_ESTIMATOR_TEST_H +#define QMCPLUSPLUS_MOCK_WALKER_ELEMENTS_FOR_ESTIMATOR_TEST_H + +#include "Particle/tests/MinimalParticlePool.h" +#include "QMCWaveFunctions/tests/MinimalWaveFunctionPool.h" +#include "QMCHamiltonians/tests/MinimalHamiltonianPool.h" +#include "Utilities/RuntimeOptions.h" +#include "Message/Communicate.h" +#include + +namespace qmcplusplus +{ +namespace testing +{ +class MockGoldWalkerElements +{ +public: + using WaveFunctionPoolFactoryFunc = + std::function; + using HamPoolFactoryFunc = + std::function; + MockGoldWalkerElements(Communicate* comm, + RuntimeOptions& runtime_opt, + WaveFunctionPoolFactoryFunc wfp_func, + HamPoolFactoryFunc ham_pool_fac_func); + + ParticleSetPool particle_pool; + ParticleSet& pset_elec; + ParticleSet& pset_ions; + WaveFunctionPool wavefunction_pool; + HamiltonianPool hamiltonian_pool; + TrialWaveFunction& twf; + QMCHamiltonian& ham; +}; + +MockGoldWalkerElements makeGoldWalkerElementsWithEE(Communicate*, RuntimeOptions run_time_opt); +MockGoldWalkerElements makeGoldWalkerElementsWithEEEI(Communicate*, RuntimeOptions run_time_opt); +} // namespace testing +} // namespace qmcplusplus + +#endif diff --git a/src/Estimators/tests/test_EnergyDensityEstimator.cpp b/src/Estimators/tests/test_EnergyDensityEstimator.cpp new file mode 100644 index 0000000000..3c19414600 --- /dev/null +++ b/src/Estimators/tests/test_EnergyDensityEstimator.cpp @@ -0,0 +1,150 @@ +////////////////////////////////////////////////////////////////////////////////////// +// This file is distributed under the University of Illinois/NCSA Open Source License. +// See LICENSE file in top directory for details. +// +// Copyright (c) 2024 QMCPACK developers. +// +// File developed by: Peter Doak, doakpw@ornl.gov, Oak Ridge National Lab +// +// File created by: Peter Doak, doakpw@ornl.gov, Oak Ridge National Lab +////////////////////////////////////////////////////////////////////////////////////// + + +#include "catch.hpp" + +#include "EnergyDensityEstimator.h" +#include +#include "EstimatorTesting.h" +#include "EnergyDensityTest.h" +#include "PerParticleHamiltonianLogger.h" +#include "ValidEnergyDensityInput.h" +#include "Particle/Walker.h" +#include "Particle/tests/MinimalParticlePool.h" +#include "QMCWaveFunctions/tests/MinimalWaveFunctionPool.h" +#include "QMCHamiltonians/tests/MinimalHamiltonianPool.h" +#include "Utilities/for_testing/NativeInitializerPrint.hpp" + +constexpr bool generate_test_data = false; + +namespace qmcplusplus +{ + +TEST_CASE("NEEnergyDensityEstimator::Constructor", "[estimators]") +{ + Communicate* comm; + comm = OHMMS::Controller; + + ParticleSetPool particle_pool{MinimalParticlePool::make_diamondC_1x1x1(comm)}; + + ParticleSet pset_elec{*(particle_pool.getParticleSet("e"))}; + ParticleSet pset_ions{*(particle_pool.getParticleSet("ion"))}; + + pset_elec.R = + ParticleSet::ParticlePos{{1.451870349, 1.381521229, 1.165202269}, {1.244515371, 1.382273176, 1.21105285}, + {0.000459944, 1.329603408, 1.265030556}, {0.748660329, 1.63420622, 1.393637791}, + {0.033228526, 1.391869137, 0.654413566}, {1.114198787, 1.654334594, 0.231075822}, + {1.657151589, 0.883870516, 1.201243939}, {0.97317591, 1.245644974, 0.284564732}}; + + Libxml2Document doc; + using Input = testing::EnergyDensityInputs; + bool okay = doc.parseFromString(Input::getXml(Input::valid::CELL)); + xmlNodePtr node = doc.getRoot(); + EnergyDensityInput edein{node}; + { + NEEnergyDensityEstimator e_den_est(edein, particle_pool.getPool()); + } +} + +TEST_CASE("NEEnergyDensityEstimator::spawnCrowdClone", "[estimators]") +{ + Communicate* comm; + comm = OHMMS::Controller; + + ParticleSetPool particle_pool{MinimalParticlePool::make_diamondC_1x1x1(comm)}; + + ParticleSet pset_elec{*(particle_pool.getParticleSet("e"))}; + ParticleSet pset_ions{*(particle_pool.getParticleSet("ion"))}; + + pset_elec.R = + ParticleSet::ParticlePos{{1.451870349, 1.381521229, 1.165202269}, {1.244515371, 1.382273176, 1.21105285}, + {0.000459944, 1.329603408, 1.265030556}, {0.748660329, 1.63420622, 1.393637791}, + {0.033228526, 1.391869137, 0.654413566}, {1.114198787, 1.654334594, 0.231075822}, + {1.657151589, 0.883870516, 1.201243939}, {0.97317591, 1.245644974, 0.284564732}}; + + Libxml2Document doc; + using Input = testing::EnergyDensityInputs; + bool okay = doc.parseFromString(Input::getXml(Input::valid::CELL)); + xmlNodePtr node = doc.getRoot(); + EnergyDensityInput edein{node}; + { + NEEnergyDensityEstimator original_e_den_est(edein, particle_pool.getPool()); + + auto clone = original_e_den_est.spawnCrowdClone(); + REQUIRE(clone != nullptr); + REQUIRE(clone.get() != &original_e_den_est); + REQUIRE(dynamic_cast(clone.get()) != nullptr); + } +} + +TEST_CASE("NEEnergyDensityEstimator::AccumulateIntegration", "[estimators]") +{ + Communicate* comm = OHMMS::Controller; + + testing::EnergyDensityTest eden_test(comm, 4 /*num_walkers*/, generate_test_data); + + auto ham_list = eden_test.getHamList(); + auto& ham_leader = ham_list.getLeader(); + auto ham_lock = ResourceCollectionTeamLock(eden_test.getHamRes(), ham_list); + eden_test.getEnergyDensityEstimator().registerListeners(ham_leader); + + PerParticleHamiltonianLogger pph_logger({}, 0); + pph_logger.registerListeners(ham_leader); + + auto pset_list = eden_test.getPSetList(); + auto pset_lock = ResourceCollectionTeamLock(eden_test.getPSetRes(), pset_list); + + pset_list[0].L[0] = 1.0; + pset_list[1].L[1] = 1.0; + pset_list[2].L[2] = 1.0; + pset_list[3].L[3] = 1.0; + + ParticleSet::mw_update(pset_list); + + auto twf_list = eden_test.getTwfList(); + auto twf_lock = ResourceCollectionTeamLock(eden_test.getTwfRes(), twf_list); + TrialWaveFunction::mw_evaluateLog(twf_list, pset_list); + QMCHamiltonian::mw_evaluate(ham_list, twf_list, pset_list); + + hdf_archive hd; + std::string test_file{"ede_test.hdf"}; + bool okay = hd.create(test_file); + REQUIRE(okay); + std::vector h5desc; + + auto& e_den_est = eden_test.getEnergyDensityEstimator(); + e_den_est.registerOperatorEstimator(hd); + + StdRandom rng; + rng.init(101); + + e_den_est.accumulate(eden_test.getWalkerList(), pset_list, twf_list, ham_list, rng); + auto spacegrids = e_den_est.getSpaceGrids(); + + decltype(spacegrids)::value_type::type& grid = spacegrids[0]; + + double summed_grid = 0; + // grid memory layout is (W)eight (T) Kinetic (V) potential + for (int i = 0; i < 16000; i++) + summed_grid += *(grid.getDataVector().begin() + i * 3 + 2) + *(grid.getDataVector().begin() + i * 3 + 1); + + auto expected_sum = pph_logger.sumOverAll(); + //Here we check the sum of logged energies against the total energy in the grid. + CHECK(summed_grid == Approx(expected_sum)); + + e_den_est.write(hd); + std::cout << "wrote success\n"; +} + +TEST_CASE("NEEnergyDensityEstimator::Collect", "[estimators]") {} + +} // namespace qmcplusplus diff --git a/src/Estimators/tests/test_OneBodyDensityMatrices.cpp b/src/Estimators/tests/test_OneBodyDensityMatrices.cpp index b637679596..cd3a90ef36 100644 --- a/src/Estimators/tests/test_OneBodyDensityMatrices.cpp +++ b/src/Estimators/tests/test_OneBodyDensityMatrices.cpp @@ -287,7 +287,7 @@ TEST_CASE("OneBodyDensityMatrices::generateSamples", "[estimators]") auto& species_set = pset_target.getSpeciesSet(); auto& spo_map = wavefunction_pool.getWaveFunction("wavefunction")->getSPOMap(); - auto samplingCaseRunner = [&pset_target, &species_set, &spo_map](Input::valid test_case) { + auto samplingCaseRunner = [&input, &pset_target, &species_set, &spo_map](Input::valid test_case) { Libxml2Document doc; bool okay = doc.parseFromString(Input::getXml(test_case)); if (!okay) @@ -435,7 +435,7 @@ TEST_CASE("OneBodyDensityMatrices::accumulate", "[estimators]") {2.535993099, 1.637133598, 3.689830303}, }}; std::vector psets = - testing::generateRandomParticleSets(pset_target, pset_source, deterministic_rs, nwalkers); + testing::generateRandomParticleSets(pset_target, pset_source, deterministic_rs, nwalkers, generate_test_data); auto& trial_wavefunction = *(wavefunction_pool.getPrimary()); std::vector> twfcs(nwalkers);