diff --git a/src/Particle/VirtualParticleSet.cpp b/src/Particle/VirtualParticleSet.cpp index b3b619dc68..c78b4e331b 100644 --- a/src/Particle/VirtualParticleSet.cpp +++ b/src/Particle/VirtualParticleSet.cpp @@ -37,7 +37,8 @@ struct VPMultiWalkerMem : public Resource Resource* makeClone() const override { return new VPMultiWalkerMem(*this); } }; -VirtualParticleSet::VirtualParticleSet(const ParticleSet& p, int nptcl, size_t dt_count_limit) : ParticleSet(p.getSimulationCell()) +VirtualParticleSet::VirtualParticleSet(const ParticleSet& p, int nptcl, size_t dt_count_limit) + : ParticleSet(p.getSimulationCell()) { setName("virtual"); @@ -194,4 +195,54 @@ void VirtualParticleSet::mw_makeMoves(const RefVectorWithLeader& vp_list, + const RefVectorWithLeader& refp_list, + const RefVector>& deltaV_list, + const RefVector>& deltaS_list, + const RefVector>& joblist, + bool sphere) +{ + auto& vp_leader = vp_list.getLeader(); + if (!vp_leader.isSpinor()) + throw std::runtime_error( + "VirtualParticleSet::mw_makeMovesWithSpin should not be called if particle sets aren't spionor types"); + vp_leader.onSphere = sphere; + vp_leader.refPS = refp_list.getLeader(); + + const size_t nVPs = countVPs(vp_list); + auto& mw_refPctls = vp_leader.getMultiWalkerRefPctls(); + mw_refPctls.resize(nVPs); + + RefVectorWithLeader p_list(vp_leader); + p_list.reserve(vp_list.size()); + + size_t ivp = 0; + for (int iw = 0; iw < vp_list.size(); iw++) + { + VirtualParticleSet& vp(vp_list[iw]); + const std::vector& deltaV(deltaV_list[iw]); + const std::vector& deltaS(deltaS_list[iw]); + const NLPPJob& job(joblist[iw]); + + vp.onSphere = sphere; + vp.refPS = refp_list[iw]; + vp.refPtcl = job.electron_id; + vp.refSourcePtcl = job.ion_id; + assert(vp.R.size() == deltaV.size()); + assert(vp.spins.size() == deltaS.size()); + assert(vp.R.size() == vp.spins.size()); + for (size_t k = 0; k < vp.R.size(); k++, ivp++) + { + vp.R[k] = refp_list[iw].R[vp.refPtcl] + deltaV[k]; + vp.spins[k] = refp_list[iw].spins[vp.refPtcl] + deltaS[k]; + mw_refPctls[ivp] = vp.refPtcl; + } + p_list.push_back(vp); + } + assert(ivp == nVPs); + + mw_refPctls.updateTo(); + ParticleSet::mw_update(p_list); +} + } // namespace qmcplusplus diff --git a/src/Particle/VirtualParticleSet.h b/src/Particle/VirtualParticleSet.h index 6d40d3c824..ae406eacc6 100644 --- a/src/Particle/VirtualParticleSet.h +++ b/src/Particle/VirtualParticleSet.h @@ -88,7 +88,11 @@ class VirtualParticleSet : public ParticleSet * @param sphere set true if VP are on a sphere around the reference source particle * @param iat reference source particle */ - void makeMoves(const ParticleSet& refp, int jel, const std::vector& deltaV, bool sphere = false, int iat = -1); + void makeMoves(const ParticleSet& refp, + int jel, + const std::vector& deltaV, + bool sphere = false, + int iat = -1); /** move virtual particles to new postions and update distance tables * @param refp reference particle set @@ -111,6 +115,13 @@ class VirtualParticleSet : public ParticleSet const RefVector>& joblist, bool sphere); + static void mw_makeMovesWithSpin(const RefVectorWithLeader& vp_list, + const RefVectorWithLeader& p_list, + const RefVector>& deltaV_list, + const RefVector>& deltaS_list, + const RefVector>& joblist, + bool sphere); + static RefVectorWithLeader RefVectorWithLeaderParticleSet( const RefVectorWithLeader& vp_list) { diff --git a/src/QMCHamiltonians/ECPComponentBuilder.2.cpp b/src/QMCHamiltonians/ECPComponentBuilder.2.cpp index 6ffa2b03b3..ac1c9dc48c 100644 --- a/src/QMCHamiltonians/ECPComponentBuilder.2.cpp +++ b/src/QMCHamiltonians/ECPComponentBuilder.2.cpp @@ -135,7 +135,7 @@ void ECPComponentBuilder::buildSemiLocalAndLocal(std::vector& semiPt // we may not know which one is local yet. std::vector angList; - std::vector angListSO; //For spin-orbit, if it exists + std::vector angListSO; //For spin-orbit, if it exists std::vector vpsPtr; std::vector vpsoPtr; //For spin-orbit, if it exists. Lmax = -1; @@ -364,8 +364,8 @@ void ECPComponentBuilder::buildSO(const std::vector& angList, app->spline(); pp_so->add(angList[l], app); } - NumSO = angList.size(); - pp_so->Rmax_ = rmax; + NumSO = angList.size(); + pp_so->setRmax(rmax); } bool ECPComponentBuilder::parseCasino(const std::string& fname, xmlNodePtr cur) diff --git a/src/QMCHamiltonians/SOECPComponent.cpp b/src/QMCHamiltonians/SOECPComponent.cpp index 8f82f586ea..408b250b6a 100644 --- a/src/QMCHamiltonians/SOECPComponent.cpp +++ b/src/QMCHamiltonians/SOECPComponent.cpp @@ -15,36 +15,37 @@ #include "SOECPComponent.h" #include "Numerics/Ylm.h" #include "CPU/BLAS.hpp" +#include "NLPPJob.h" namespace qmcplusplus { SOECPComponent::SOECPComponent() - : lmax_(0), nchannel_(0), nknot_(0), sknot_(0), total_knots_(0), Rmax_(-1), VP_(nullptr) + : lmax_(0), nchannel_(0), nknot_(0), sknot_(0), total_knots_(0), rmax_(-1), vp_(nullptr) {} SOECPComponent::~SOECPComponent() { for (int i = 0; i < sopp_m_.size(); i++) delete sopp_m_[i]; - if (VP_) - delete VP_; + if (vp_) + delete vp_; } void SOECPComponent::print(std::ostream& os) {} void SOECPComponent::initVirtualParticle(const ParticleSet& qp) { - assert(VP_ == nullptr); + assert(vp_ == nullptr); outputManager.pause(); - VP_ = new VirtualParticleSet(qp, total_knots_); + vp_ = new VirtualParticleSet(qp, total_knots_); outputManager.resume(); } void SOECPComponent::deleteVirtualParticle() { - if (VP_) - delete VP_; - VP_ = nullptr; + if (vp_) + delete vp_; + vp_ = nullptr; } void SOECPComponent::add(int l, RadialPotentialType* pp) @@ -58,8 +59,8 @@ SOECPComponent* SOECPComponent::makeClone(const ParticleSet& qp) SOECPComponent* myclone = new SOECPComponent(*this); for (int i = 0; i < sopp_m_.size(); i++) myclone->sopp_m_[i] = sopp_m_[i]->makeClone(); - if (VP_) - myclone->VP_ = new VirtualParticleSet(qp, total_knots_); + if (vp_) + myclone->vp_ = new VirtualParticleSet(qp, total_knots_); return myclone; } @@ -141,17 +142,13 @@ SOECPComponent::RealType SOECPComponent::evaluateOne(ParticleSet& W, RealType r, const PosType& dr) { - - for (int ip = 0; ip < nchannel_; ip++) - vrad_[ip] = sopp_m_[ip]->splint(r); - RealType sold = W.spins[iel]; buildTotalQuadrature(r, dr, sold); - if (VP_) + if (vp_) { - VP_->makeMovesWithSpin(W, iel, deltaV_, deltaS_, true, iat); - Psi.evaluateRatios(*VP_, psiratio_); + vp_->makeMovesWithSpin(W, iel, deltaV_, deltaS_, true, iat); + Psi.evaluateRatios(*vp_, psiratio_); } else for (int iq = 0; iq < total_knots_; iq++) @@ -162,6 +159,11 @@ SOECPComponent::RealType SOECPComponent::evaluateOne(ParticleSet& W, Psi.resetPhaseDiff(); } + return calculateProjector(r, dr, sold); +} + +SOECPComponent::RealType SOECPComponent::calculateProjector(RealType r, const PosType& dr, RealType sold) +{ ComplexType pairpot; for (int iq = 0; iq < total_knots_; iq++) { @@ -190,6 +192,81 @@ SOECPComponent::RealType SOECPComponent::evaluateOne(ParticleSet& W, return std::real(pairpot); } +void SOECPComponent::mw_evaluateOne(const RefVectorWithLeader& soecp_component_list, + const RefVectorWithLeader& p_list, + const RefVectorWithLeader& psi_list, + const RefVector>& joblist, + std::vector& pairpots, + ResourceCollection& collection) +{ + auto& soecp_component_leader = soecp_component_list.getLeader(); + if (soecp_component_leader.vp_) + { + // Compute ratios with VP + RefVectorWithLeader vp_list(*soecp_component_leader.vp_); + RefVectorWithLeader const_vp_list(*soecp_component_leader.vp_); + RefVector> deltaV_list; + RefVector> deltaS_list; + RefVector> psiratios_list; + vp_list.reserve(soecp_component_list.size()); + const_vp_list.reserve(soecp_component_list.size()); + deltaV_list.reserve(soecp_component_list.size()); + deltaS_list.reserve(soecp_component_list.size()); + psiratios_list.reserve(soecp_component_list.size()); + + for (size_t i = 0; i < soecp_component_list.size(); i++) + { + SOECPComponent& component(soecp_component_list[i]); + const NLPPJob& job = joblist[i]; + const RealType sold = p_list[i].spins[job.electron_id]; + + component.buildTotalQuadrature(job.ion_elec_dist, job.ion_elec_displ, sold); + + vp_list.push_back(*component.vp_); + const_vp_list.push_back(*component.vp_); + deltaV_list.push_back(component.deltaV_); + deltaS_list.push_back(component.deltaS_); + psiratios_list.push_back(component.psiratio_); + } + + ResourceCollectionTeamLock vp_res_lock(collection, vp_list); + + VirtualParticleSet::mw_makeMovesWithSpin(vp_list, p_list, deltaV_list, deltaS_list, joblist, true); + + TrialWaveFunction::mw_evaluateRatios(psi_list, const_vp_list, psiratios_list); + } + else + { + // Compute ratios without VP. Slow + for (size_t i = 0; i < p_list.size(); i++) + { + SOECPComponent& component(soecp_component_list[i]); + ParticleSet& W(p_list[i]); + TrialWaveFunction& psi(psi_list[i]); + const NLPPJob& job = joblist[i]; + + const RealType sold = W.spins[job.electron_id]; + component.buildTotalQuadrature(job.ion_elec_dist, job.ion_elec_displ, sold); + + for (int j = 0; j < component.total_knots_; j++) + { + W.makeMoveWithSpin(job.electron_id, component.deltaV_[j], component.deltaS_[j]); + component.psiratio_[j] = psi.calcRatio(W, job.electron_id); + W.rejectMove(job.electron_id); + psi.resetPhaseDiff(); + } + } + } + + for (size_t i = 0; i < p_list.size(); i++) + { + SOECPComponent& component(soecp_component_list[i]); + const NLPPJob& job = joblist[i]; + const RealType sold = p_list[i].spins[job.electron_id]; + pairpots[i] = component.calculateProjector(job.ion_elec_dist, job.ion_elec_displ, sold); + } +} + SOECPComponent::RealType SOECPComponent::evaluateValueAndDerivatives(ParticleSet& W, int iat, TrialWaveFunction& Psi, @@ -214,11 +291,11 @@ SOECPComponent::RealType SOECPComponent::evaluateValueAndDerivatives(ParticleSet //Now we have all the spin and spatial quadrature points acculated to use in evaluation //Now we need to obtain dlogpsi and dlogpsi_vp - if (VP_) + if (vp_) { // Compute ratios with VP - VP_->makeMovesWithSpin(W, iel, deltaV_, deltaS_, true, iat); - Psi.evaluateDerivRatios(*VP_, optvars, psiratio_, dratio_); + vp_->makeMovesWithSpin(W, iel, deltaV_, deltaS_, true, iat); + Psi.evaluateDerivRatios(*vp_, optvars, psiratio_, dratio_); } else for (int iq = 0; iq < total_knots_; iq++) @@ -317,6 +394,10 @@ void SOECPComponent::buildTotalQuadrature(const RealType r, const PosType& dr, c addSpatialQuadrature(sknot_, r, dr, smax - sold, RealType(1. / 3.) * dS); assert(count == total_knots_); + + //also set the radial function + for (int ip = 0; ip < nchannel_; ip++) + vrad_[ip] = sopp_m_[ip]->splint(r); } void SOECPComponent::rotateQuadratureGrid(const TensorType& rmat) diff --git a/src/QMCHamiltonians/SOECPComponent.h b/src/QMCHamiltonians/SOECPComponent.h index a1ee8796a7..66ed70e305 100644 --- a/src/QMCHamiltonians/SOECPComponent.h +++ b/src/QMCHamiltonians/SOECPComponent.h @@ -16,6 +16,7 @@ #include "QMCHamiltonians/OperatorBase.h" #include "QMCHamiltonians/RandomRotationMatrix.h" #include "QMCWaveFunctions/TrialWaveFunction.h" +#include #include "Numerics/OneDimGridBase.h" #include "Numerics/OneDimGridFunctor.h" #include "Numerics/OneDimLinearSpline.h" @@ -23,6 +24,11 @@ namespace qmcplusplus { + +namespace testing +{ +class TestSOECPotential; +} /** class SOECPComponent ** brief Computes the nonlocal spin-orbit interaction \f$\Delta V_SO(r) |ljm_j> angpp_m_; ///Non-Local part of the pseudo-potential @@ -70,9 +76,9 @@ class SOECPComponent : public QMCTraits //scratch spaces used by evaluateValueAndDerivative Matrix dratio_; Vector dlogpsi_vp_; - VirtualParticleSet* VP_; + VirtualParticleSet* vp_; - //This builds the full quadrature grid for the Simpsons rule used for spin integrals as well as + //This builds the full quadrature grid for the Simpsons rule used for spin integrals as well as //the spatial quadrature. In this function, it specifies the deltaS_ and deltaV_ for all the quadrature points and sets the interal weights //in spin_quad_weights //If there are s0,s1,...sN spin integral points and q0,q1,...qM spatial quadrature points, the order is @@ -108,6 +114,16 @@ class SOECPComponent : public QMCTraits */ RealType evaluateOne(ParticleSet& W, int iat, TrialWaveFunction& Psi, int iel, RealType r, const PosType& dr); + RealType calculateProjector(RealType r, const PosType& dr, RealType sold); + + + static void mw_evaluateOne(const RefVectorWithLeader& soecp_component_list, + const RefVectorWithLeader& p_list, + const RefVectorWithLeader& psi_list, + const RefVector>& joblist, + std::vector& pairpots, + ResourceCollection& collection); + RealType evaluateValueAndDerivatives(ParticleSet& P, int iat, TrialWaveFunction& psi, @@ -123,15 +139,19 @@ class SOECPComponent : public QMCTraits void initVirtualParticle(const ParticleSet& qp); void deleteVirtualParticle(); - inline void setRmax(RealType rmax) { Rmax_ = rmax; } - inline RealType getRmax() const { return Rmax_; } - inline void setLmax(int Lmax) { lmax_ = Lmax; } + inline void setRmax(RealType rmax) { rmax_ = rmax; } + inline RealType getRmax() const { return rmax_; } + inline void setLmax(int lmax) { lmax_ = lmax; } inline int getLmax() const { return lmax_; } inline int getNknot() const { return nknot_; } inline int getSknot() const { return sknot_; } + const VirtualParticleSet* getVP() const { return vp_; }; + friend struct ECPComponentBuilder; friend void copyGridUnrotatedForTest(SOECPComponent& nlpp); + + friend class testing::TestSOECPotential; }; } // namespace qmcplusplus diff --git a/src/QMCHamiltonians/SOECPotential.cpp b/src/QMCHamiltonians/SOECPotential.cpp index e680e8968a..3adaf2c60c 100644 --- a/src/QMCHamiltonians/SOECPotential.cpp +++ b/src/QMCHamiltonians/SOECPotential.cpp @@ -12,55 +12,89 @@ #include "Particle/DistanceTable.h" #include "SOECPotential.h" #include "Utilities/IteratorUtility.h" +#include "NLPPJob.h" +#include +#include namespace qmcplusplus { + +struct SOECPotential::SOECPotentialMultiWalkerResource : public Resource +{ + SOECPotentialMultiWalkerResource() : Resource("SOECPotential") {} + + Resource* makeClone() const override; + + ResourceCollection collection{"SOPPcollection"}; + + //crowds worth of per particle soecp values + Matrix ve_samples; + Matrix vi_samples; +}; + /** constructor *\param ionic positions *\param els electronic poitions *\param psi Trial wave function */ SOECPotential::SOECPotential(ParticleSet& ions, ParticleSet& els, TrialWaveFunction& psi) - : myRNG(nullptr), IonConfig(ions), Psi(psi), Peln(els), ElecNeighborIons(els), IonNeighborElecs(ions) + : my_rng_(nullptr), ion_config_(ions), psi_(psi), peln_(els), elec_neighbor_ions_(els), ion_neighbor_elecs_(ions) { setEnergyDomain(POTENTIAL); twoBodyQuantumDomain(ions, els); - myTableIndex = els.addTable(ions); - NumIons = ions.getTotalNum(); - PP.resize(NumIons, nullptr); - PPset.resize(IonConfig.getSpeciesSet().getTotalNum()); + my_table_index_ = els.addTable(ions); + num_ions_ = ions.getTotalNum(); + pp_.resize(num_ions_, nullptr); + ppset_.resize(ion_config_.getSpeciesSet().getTotalNum()); + sopp_jobs_.resize(els.groups()); + for (size_t ig = 0; ig < els.groups(); ig++) + sopp_jobs_[ig].reserve(2 * els.groupsize(ig)); } +SOECPotential::~SOECPotential() = default; + void SOECPotential::resetTargetParticleSet(ParticleSet& P) {} SOECPotential::Return_t SOECPotential::evaluate(ParticleSet& P) +{ + evaluateImpl(P, false); + return value_; +} + +SOECPotential::Return_t SOECPotential::evaluateDeterministic(ParticleSet& P) +{ + evaluateImpl(P, true); + return value_; +} + +void SOECPotential::evaluateImpl(ParticleSet& P, bool keep_grid) { value_ = 0.0; - for (int ipp = 0; ipp < PPset.size(); ipp++) - if (PPset[ipp]) - PPset[ipp]->rotateQuadratureGrid(generateRandomRotationMatrix(*myRNG)); + if (!keep_grid) + for (int ipp = 0; ipp < ppset_.size(); ipp++) + if (ppset_[ipp]) + ppset_[ipp]->rotateQuadratureGrid(generateRandomRotationMatrix(*my_rng_)); - const auto& myTable = P.getDistTableAB(myTableIndex); - for (int iat = 0; iat < NumIons; iat++) - IonNeighborElecs.getNeighborList(iat).clear(); + const auto& ble = P.getDistTableAB(my_table_index_); + for (int iat = 0; iat < num_ions_; iat++) + ion_neighbor_elecs_.getNeighborList(iat).clear(); for (int jel = 0; jel < P.getTotalNum(); jel++) - ElecNeighborIons.getNeighborList(jel).clear(); + elec_neighbor_ions_.getNeighborList(jel).clear(); for (int jel = 0; jel < P.getTotalNum(); jel++) { - const auto& dist = myTable.getDistRow(jel); - const auto& displ = myTable.getDisplRow(jel); - std::vector& NeighborIons = ElecNeighborIons.getNeighborList(jel); - for (int iat = 0; iat < NumIons; iat++) - if (PP[iat] != nullptr && dist[iat] < PP[iat]->getRmax()) + const auto& dist = ble.getDistRow(jel); + const auto& displ = ble.getDisplRow(jel); + std::vector& NeighborIons = elec_neighbor_ions_.getNeighborList(jel); + for (int iat = 0; iat < num_ions_; iat++) + if (pp_[iat] != nullptr && dist[iat] < pp_[iat]->getRmax()) { - RealType pairpot = PP[iat]->evaluateOne(P, iat, Psi, jel, dist[iat], -displ[iat]); + RealType pairpot = pp_[iat]->evaluateOne(P, iat, psi_, jel, dist[iat], -displ[iat]); value_ += pairpot; NeighborIons.push_back(iat); - IonNeighborElecs.getNeighborList(iat).push_back(jel); + ion_neighbor_elecs_.getNeighborList(iat).push_back(jel); } } - return value_; } SOECPotential::Return_t SOECPotential::evaluateValueAndDerivatives(ParticleSet& P, @@ -69,38 +103,237 @@ SOECPotential::Return_t SOECPotential::evaluateValueAndDerivatives(ParticleSet& Vector& dhpsioverpsi) { value_ = 0.0; - for (int ipp = 0; ipp < PPset.size(); ipp++) - if (PPset[ipp]) - PPset[ipp]->rotateQuadratureGrid(generateRandomRotationMatrix(*myRNG)); + for (int ipp = 0; ipp < ppset_.size(); ipp++) + if (ppset_[ipp]) + ppset_[ipp]->rotateQuadratureGrid(generateRandomRotationMatrix(*my_rng_)); - const auto& myTable = P.getDistTableAB(myTableIndex); + const auto& ble = P.getDistTableAB(my_table_index_); for (int jel = 0; jel < P.getTotalNum(); jel++) { - const auto& dist = myTable.getDistRow(jel); - const auto& displ = myTable.getDisplRow(jel); - for (int iat = 0; iat < NumIons; iat++) - if (PP[iat] != nullptr && dist[iat] < PP[iat]->getRmax()) - value_ += PP[iat]->evaluateValueAndDerivatives(P, iat, Psi, jel, dist[iat], -displ[iat], optvars, dlogpsi, - dhpsioverpsi); + const auto& dist = ble.getDistRow(jel); + const auto& displ = ble.getDisplRow(jel); + for (int iat = 0; iat < num_ions_; iat++) + if (pp_[iat] != nullptr && dist[iat] < pp_[iat]->getRmax()) + value_ += pp_[iat]->evaluateValueAndDerivatives(P, iat, psi_, jel, dist[iat], -displ[iat], optvars, dlogpsi, + dhpsioverpsi); } return value_; } +void SOECPotential::mw_evaluate(const RefVectorWithLeader& o_list, + const RefVectorWithLeader& wf_list, + const RefVectorWithLeader& p_list) const +{ + mw_evaluateImpl(o_list, wf_list, p_list, std::nullopt); +} + +void SOECPotential::mw_evaluatePerParticle(const RefVectorWithLeader& o_list, + const RefVectorWithLeader& wf_list, + const RefVectorWithLeader& p_list, + const std::vector>& listeners, + const std::vector>& listeners_ions) const +{ + std::optional> l_opt(std::in_place, listeners, listeners_ions); + mw_evaluateImpl(o_list, wf_list, p_list, l_opt); +} + +void SOECPotential::mw_evaluateImpl(const RefVectorWithLeader& o_list, + const RefVectorWithLeader& wf_list, + const RefVectorWithLeader& p_list, + const std::optional> listeners, + bool keep_grid) +{ + auto& O_leader = o_list.getCastedLeader(); + ParticleSet& pset_leader = p_list.getLeader(); + const size_t nw = o_list.size(); + + for (size_t iw = 0; iw < nw; iw++) + { + auto& O = o_list.getCastedElement(iw); + const ParticleSet& P(p_list[iw]); + + if (!keep_grid) + for (size_t ipp = 0; ipp < O.ppset_.size(); ipp++) + if (O.ppset_[ipp]) + O.ppset_[ipp]->rotateQuadratureGrid(generateRandomRotationMatrix(*O.my_rng_)); + + //loop over all the ions + const auto& ble = P.getDistTableAB(O.my_table_index_); + //clear elec and ion neighbor lists + for (size_t iat = 0; iat < O.num_ions_; iat++) + O.ion_neighbor_elecs_.getNeighborList(iat).clear(); + for (size_t jel = 0; jel < P.getTotalNum(); jel++) + O.elec_neighbor_ions_.getNeighborList(jel).clear(); + + for (size_t ig = 0; ig < P.groups(); ig++) // loop over species + { + auto& joblist = O.sopp_jobs_[ig]; + joblist.clear(); + + for (size_t jel = P.first(ig); jel < P.last(ig); jel++) + { + const auto& dist = ble.getDistRow(jel); + const auto& displ = ble.getDisplRow(jel); + std::vector& NeighborIons = O.elec_neighbor_ions_.getNeighborList(jel); + for (size_t iat = 0; iat < O.num_ions_; iat++) + if (O.pp_[iat] != nullptr && dist[iat] < O.pp_[iat]->getRmax()) + { + NeighborIons.push_back(iat); + O.ion_neighbor_elecs_.getNeighborList(iat).push_back(jel); + joblist.emplace_back(iat, jel, dist[iat], -displ[iat]); + } + } + } + O.value_ = 0.0; + } + + if (listeners) + { + auto& ve_samples = O_leader.mw_res_->ve_samples; + auto& vi_samples = O_leader.mw_res_->vi_samples; + ve_samples.resize(nw, pset_leader.getTotalNum()); + vi_samples.resize(nw, O_leader.ion_config_.getTotalNum()); + } + + auto pp_component = std::find_if(O_leader.ppset_.begin(), O_leader.ppset_.end(), [](auto& ptr) { return bool(ptr); }); + assert(pp_component != std::end(O_leader.ppset_)); + + RefVector soecp_potential_list; + RefVectorWithLeader soecp_component_list(**pp_component); + RefVectorWithLeader pset_list(pset_leader); + RefVectorWithLeader psi_list(O_leader.psi_); + assert(&O_leader.psi_ == &wf_list.getLeader()); + for (size_t iw = 0; iw < nw; iw++) + assert(&o_list.getCastedElement(iw).psi_ == &wf_list[iw]); + + RefVector> batch_list; + std::vector pairpots(nw); + + soecp_potential_list.reserve(nw); + soecp_component_list.reserve(nw); + pset_list.reserve(nw); + psi_list.reserve(nw); + batch_list.reserve(nw); + + for (size_t ig = 0; ig < pset_leader.groups(); ig++) + { + TrialWaveFunction::mw_prepareGroup(wf_list, p_list, ig); + + size_t max_num_jobs = 0; + for (size_t iw = 0; iw < nw; iw++) + { + const auto& O = o_list.getCastedElement(iw); + max_num_jobs = std::max(max_num_jobs, O.sopp_jobs_[ig].size()); + } + + for (size_t jobid = 0; jobid < max_num_jobs; jobid++) + { + soecp_potential_list.clear(); + soecp_component_list.clear(); + pset_list.clear(); + psi_list.clear(); + batch_list.clear(); + for (size_t iw = 0; iw < nw; iw++) + { + auto& O = o_list.getCastedElement(iw); + if (jobid < O.sopp_jobs_[ig].size()) + { + const auto& job = O.sopp_jobs_[ig][jobid]; + soecp_potential_list.push_back(O); + soecp_component_list.push_back(*O.pp_[job.ion_id]); + pset_list.push_back(p_list[iw]); + psi_list.push_back(wf_list[iw]); + batch_list.push_back(job); + } + } + + SOECPComponent::mw_evaluateOne(soecp_component_list, pset_list, psi_list, batch_list, pairpots, + O_leader.mw_res_->collection); + + for (size_t j = 0; j < soecp_potential_list.size(); j++) + { + soecp_potential_list[j].get().value_ += pairpots[j]; + + if (listeners) + { + auto& ve_samples = O_leader.mw_res_->ve_samples; + auto& vi_samples = O_leader.mw_res_->vi_samples; + int iw = j; + ve_samples(iw, batch_list[j].get().electron_id) += pairpots[j]; + vi_samples(iw, batch_list[j].get().ion_id) += pairpots[j]; + } + } + } + } + + if (listeners) + { + auto& ve_samples = O_leader.mw_res_->ve_samples; + auto& vi_samples = O_leader.mw_res_->vi_samples; + int num_electrons = pset_leader.getTotalNum(); + for (int iw = 0; iw < nw; iw++) + { + Vector ve_sample(ve_samples.begin(iw), num_electrons); + Vector vi_sample(vi_samples.begin(iw), O_leader.num_ions_); + for (const ListenerVector& listener : listeners->electron_values) + listener.report(iw, O_leader.getName(), ve_sample); + for (const ListenerVector& listener : listeners->ion_values) + listener.report(iw, O_leader.getName(), vi_sample); + } + ve_samples = 0.0; + vi_samples = 0.0; + } +} + std::unique_ptr SOECPotential::makeClone(ParticleSet& qp, TrialWaveFunction& psi) { - std::unique_ptr myclone = std::make_unique(IonConfig, qp, psi); - for (int ig = 0; ig < PPset.size(); ++ig) - if (PPset[ig]) - myclone->addComponent(ig, std::unique_ptr(PPset[ig]->makeClone(qp))); + std::unique_ptr myclone = std::make_unique(ion_config_, qp, psi); + for (int ig = 0; ig < ppset_.size(); ++ig) + if (ppset_[ig]) + myclone->addComponent(ig, std::unique_ptr(ppset_[ig]->makeClone(qp))); return myclone; } void SOECPotential::addComponent(int groupID, std::unique_ptr&& ppot) { - for (int iat = 0; iat < PP.size(); iat++) - if (IonConfig.GroupID[iat] == groupID) - PP[iat] = ppot.get(); - PPset[groupID] = std::move(ppot); + for (int iat = 0; iat < pp_.size(); iat++) + if (ion_config_.GroupID[iat] == groupID) + pp_[iat] = ppot.get(); + ppset_[groupID] = std::move(ppot); +} + +void SOECPotential::createResource(ResourceCollection& collection) const +{ + auto new_res = std::make_unique(); + for (int ig = 0; ig < ppset_.size(); ig++) + if (ppset_[ig]->getVP()) + { + ppset_[ig]->getVP()->createResource(new_res->collection); + break; + } + auto resource_index = collection.addResource(std::move(new_res)); +} + +void SOECPotential::acquireResource(ResourceCollection& collection, + const RefVectorWithLeader& o_list) const +{ + auto& O_leader = o_list.getCastedLeader(); + auto res_ptr = dynamic_cast(collection.lendResource().release()); + if (!res_ptr) + throw std::runtime_error("SOECPotential::acquireResource dynamic_cast failed"); + O_leader.mw_res_.reset(res_ptr); +} + +void SOECPotential::releaseResource(ResourceCollection& collection, + const RefVectorWithLeader& o_list) const +{ + auto& O_leader = o_list.getCastedLeader(); + collection.takebackResource(std::move(O_leader.mw_res_)); +} + +Resource* SOECPotential::SOECPotentialMultiWalkerResource::makeClone() const +{ + return new SOECPotentialMultiWalkerResource(*this); } } // namespace qmcplusplus diff --git a/src/QMCHamiltonians/SOECPotential.h b/src/QMCHamiltonians/SOECPotential.h index b6d5801dae..60471cf246 100644 --- a/src/QMCHamiltonians/SOECPotential.h +++ b/src/QMCHamiltonians/SOECPotential.h @@ -17,27 +17,50 @@ namespace qmcplusplus { +template +struct NLPPJob; + +namespace testing +{ +class TestSOECPotential; +} + class SOECPotential : public OperatorBase { + struct SOECPotentialMultiWalkerResource; + using Real = QMCTraits::RealType; + public: SOECPotential(ParticleSet& ions, ParticleSet& els, TrialWaveFunction& psi); + ~SOECPotential() override; bool dependsOnWaveFunction() const override { return true; } std::string getClassName() const override { return "SOECPotential"; } void resetTargetParticleSet(ParticleSet& P) override; Return_t evaluate(ParticleSet& P) override; + Return_t evaluateDeterministic(ParticleSet& P) override; Return_t evaluateValueAndDerivatives(ParticleSet& P, const opt_variables_type& optvars, const Vector& dlogpsi, Vector& dhpsioverpsi) override; + void mw_evaluate(const RefVectorWithLeader& o_list, + const RefVectorWithLeader& wf_list, + const RefVectorWithLeader& p_list) const override; + + void mw_evaluatePerParticle(const RefVectorWithLeader& o_list, + const RefVectorWithLeader& wf_list, + const RefVectorWithLeader& p_list, + const std::vector>& listeners, + const std::vector>& listeners_ions) const override; + bool put(xmlNodePtr cur) override { return true; } bool get(std::ostream& os) const override { - os << "SOECPotential: " << IonConfig.getName(); + os << "SOECPotential: " << ion_config_.getName(); return true; } @@ -45,27 +68,51 @@ class SOECPotential : public OperatorBase void addComponent(int groupID, std::unique_ptr&& pp); - void setRandomGenerator(RandomGenerator* rng) override { myRNG = rng; } + void setRandomGenerator(RandomGenerator* rng) override { my_rng_ = rng; } + + //initialize shared resource and hand to collection + void createResource(ResourceCollection& collection) const override; + + //acquire a shared resource + void acquireResource(ResourceCollection& collection, const RefVectorWithLeader& o_list) const override; + + //return a shared resource to a collection + void releaseResource(ResourceCollection& collection, const RefVectorWithLeader& o_list) const override; protected: - RandomGenerator* myRNG; - std::vector PP; - std::vector> PPset; - ParticleSet& IonConfig; - TrialWaveFunction& Psi; + RandomGenerator* my_rng_; + std::vector pp_; + std::vector> ppset_; + ParticleSet& ion_config_; + TrialWaveFunction& psi_; + static void mw_evaluateImpl(const RefVectorWithLeader& o_list, + const RefVectorWithLeader& wf_list, + const RefVectorWithLeader& p_list, + std::optional> listeners, + bool keep_grid = false); private: ///number of ions - int NumIons; + int num_ions_; ///index of distance table for ion-el pair - int myTableIndex; + int my_table_index_; ///reference to the electrons - ParticleSet& Peln; + ParticleSet& peln_; ///neighborlist of electrons - NeighborLists ElecNeighborIons; + NeighborLists elec_neighbor_ions_; ///neighborlist of ions - NeighborLists IonNeighborElecs; + NeighborLists ion_neighbor_elecs_; + //job list for evaluation + std::vector>> sopp_jobs_; + //multi walker resource + std::unique_ptr mw_res_; + + void evaluateImpl(ParticleSet& elec, bool keep_grid = false); + + //for testing + friend class testing::TestSOECPotential; }; + } // namespace qmcplusplus #endif diff --git a/src/QMCHamiltonians/tests/CMakeLists.txt b/src/QMCHamiltonians/tests/CMakeLists.txt index 84fd6ed391..6506cdb662 100644 --- a/src/QMCHamiltonians/tests/CMakeLists.txt +++ b/src/QMCHamiltonians/tests/CMakeLists.txt @@ -32,6 +32,8 @@ set(HAM_SRCS if(NOT QMC_COMPLEX) set(HAM_SRCS ${HAM_SRCS} test_RotatedSPOs_NLPP.cpp) +else() + set(HAM_SRCS ${HAM_SRCS} test_SOECPotential.cpp) endif() set(FORCE_SRCS ${FORCE_SRCS} test_ion_derivs.cpp) diff --git a/src/QMCHamiltonians/tests/test_SOECPotential.cpp b/src/QMCHamiltonians/tests/test_SOECPotential.cpp new file mode 100644 index 0000000000..69359eccb8 --- /dev/null +++ b/src/QMCHamiltonians/tests/test_SOECPotential.cpp @@ -0,0 +1,264 @@ +////////////////////////////////////////////////////////////////////////////////////// +// 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: Cody A. Melton, cmelton@sandia.gov, Sandia National Laboratories +// +// File created by: Cody A. Melton, cmelton@sandia.gov, Sandia Nationaln Laboratories +////////////////////////////////////////////////////////////////////////////////////// + +#include "catch.hpp" + +#include "Configuration.h" +#include "QMCHamiltonians/SOECPotential.h" +#include "QMCWaveFunctions/ElectronGas/FreeOrbital.h" +#include "QMCWaveFunctions/SpinorSet.h" +#include "QMCWaveFunctions/Fermion/SlaterDet.h" +#include "QMCWaveFunctions/Fermion/DiracDeterminantBatched.h" +#include "QMCWaveFunctions/Jastrow/BsplineFunctor.h" +#include "QMCWaveFunctions/Jastrow/RadialJastrowBuilder.h" +#include "QMCHamiltonians/ECPComponentBuilder.h" +#include "TestListenerFunction.h" + +namespace qmcplusplus +{ +namespace testing +{ +class TestSOECPotential +{ + using Real = QMCTraits::RealType; + +public: + static void copyGridUnrotatedForTest(SOECPotential& so_ecp) + { + so_ecp.ppset_[0]->rrotsgrid_m_ = so_ecp.ppset_[0]->sgridxyz_m_; + } + static bool didGridChange(SOECPotential& so_ecp) + { + return so_ecp.ppset_[0]->rrotsgrid_m_ != so_ecp.ppset_[0]->sgridxyz_m_; + } + static void addVPs(const RefVectorWithLeader& o_list, const RefVectorWithLeader& p_list) + { + for (size_t iw = 0; iw < o_list.size(); iw++) + { + auto& sopp = o_list.getCastedElement(iw); + auto& pset = p_list[iw]; + for (auto& uptr_comp : sopp.ppset_) + uptr_comp.get()->initVirtualParticle(pset); + } + } + static void mw_evaluateImpl(SOECPotential& so_ecp, + const RefVectorWithLeader& o_list, + const RefVectorWithLeader& twf_list, + const RefVectorWithLeader& p_list, + const std::optional> listener_opt, + bool keep_grid) + { + so_ecp.mw_evaluateImpl(o_list, twf_list, p_list, listener_opt, keep_grid); + } +}; +} // namespace testing + +void doSOECPotentialTest(bool use_VPs) +{ + using Real = QMCTraits::RealType; + using FullPrecReal = QMCTraits::FullPrecRealType; + using Value = QMCTraits::ValueType; + using Pos = QMCTraits::PosType; + using testing::getParticularListener; + + Communicate* c = OHMMS::Controller; + + //Cell definition: + + CrystalLattice lattice; + lattice.BoxBConds = false; // periodic + lattice.R.diagonal(20); + lattice.LR_dim_cutoff = 15; + lattice.reset(); + + const SimulationCell simulation_cell(lattice); + auto ions_uptr = std::make_unique(simulation_cell); + auto elec_uptr = std::make_unique(simulation_cell); + ParticleSet& ions(*ions_uptr); + ParticleSet& elec(*elec_uptr); + ions.setName("ion0"); + ions.create({1}); + ions.R[0] = {0.0, 0.0, 0.0}; + + SpeciesSet& ion_species = ions.getSpeciesSet(); + int pIdx = ion_species.addSpecies("H"); + int pChargeIdx = ion_species.addAttribute("charge"); + int iatnumber = ion_species.addAttribute("atomic_number"); + ion_species(pChargeIdx, pIdx) = 0; + ion_species(iatnumber, pIdx) = 1; + ions.createSK(); + + elec.setName("e"); + elec.setSpinor(true); + elec.create({2}); + elec.R[0] = {0.138, -0.24, 0.216}; + elec.R[1] = {-0.216, 0.24, -0.138}; + elec.spins = {0.2, 0.51}; + + SpeciesSet& tspecies = elec.getSpeciesSet(); + int upIdx = tspecies.addSpecies("u"); + int chargeIdx = tspecies.addAttribute("charge"); + int massIdx = tspecies.addAttribute("mass"); + tspecies(chargeIdx, upIdx) = -1; + tspecies(massIdx, upIdx) = 1.0; + + elec.createSK(); + + ions.resetGroups(); + elec.resetGroups(); + + ParticleSet elec2(elec); + elec2.update(); + + RefVector ptcls{elec, elec2}; + RefVectorWithLeader p_list(elec, ptcls); + ResourceCollection pset_res("test_pset_res"); + ResourceCollectionTeamLock mw_pset_lock(pset_res, p_list); + + TrialWaveFunction psi; + + std::vector kup, kdn; + std::vector k2up, k2dn; + QMCTraits::IndexType nelec = elec.getTotalNum(); + REQUIRE(nelec == 2); + + kup.resize(nelec); + kup[0] = Pos(1, 1, 1); + kup[1] = Pos(2, 2, 2); + + kdn.resize(nelec); + kdn[0] = Pos(2, 2, 2); + kdn[1] = Pos(1, 1, 1); + + auto spo_up = std::make_unique("free_orb_up", kup); + auto spo_dn = std::make_unique("free_orb_dn", kdn); + + auto spinor_set = std::make_unique("free_orb_spinor"); + spinor_set->set_spos(std::move(spo_up), std::move(spo_dn)); + QMCTraits::IndexType norb = spinor_set->getOrbitalSetSize(); + REQUIRE(norb == 2); + + auto dd = std::make_unique>(std::move(spinor_set), 0, nelec); + std::vector> dirac_dets; + dirac_dets.push_back(std::move(dd)); + auto sd = std::make_unique(elec, std::move(dirac_dets)); + psi.addComponent(std::move(sd)); + + //Add the two body jastrow, parameters from test_J2_bspline + //adding jastrow will allow for adding WF parameter derivatives since FreeOrbital doesn't + //support that + const char* particles = R"( + + + 0.02904699284 -0.1004179 -0.1752703883 -0.2232576505 -0.2728029201 + + + + )"; + Libxml2Document doc; + bool okay = doc.parseFromString(particles); + REQUIRE(okay); + xmlNodePtr root = doc.getRoot(); + xmlNodePtr jas2 = xmlFirstElementChild(root); + RadialJastrowBuilder jastrow(c, elec); + psi.addComponent(jastrow.buildComponent(jas2)); + + std::unique_ptr psi_clone = psi.makeClone(elec2); + RefVectorWithLeader twf_list(psi, {psi, *psi_clone}); + ResourceCollection twf_res("test_twf_res"); + psi.createResource(twf_res); + ResourceCollectionTeamLock mw_twf_lock(twf_res, twf_list); + + //Now we set up the SO ECP component. + SOECPotential so_ecp(ions, elec, psi); + ECPComponentBuilder ecp_comp_builder("test_read_soecp", c); + okay = ecp_comp_builder.read_pp_file("so_ecp_test.xml"); + REQUIRE(okay); + UPtr so_ecp_comp = std::move(ecp_comp_builder.pp_so); + so_ecp.addComponent(0, std::move(so_ecp_comp)); + UPtr so_ecp2_ptr = so_ecp.makeClone(elec2, *psi_clone); + auto& so_ecp2 = dynamic_cast(*so_ecp2_ptr); + + StdRandom rng(10101); + StdRandom rng2(10201); + so_ecp.setRandomGenerator(&rng); + so_ecp2.setRandomGenerator(&rng2); + + RefVector so_ecps{so_ecp, so_ecp2}; + RefVectorWithLeader o_list(so_ecp, so_ecps); + if (use_VPs) + testing::TestSOECPotential::addVPs(o_list, p_list); + ResourceCollection so_ecp_res("test_so_ecp_res"); + so_ecp.createResource(so_ecp_res); + ResourceCollectionTeamLock so_ecp_lock(so_ecp_res, o_list); + + testing::TestSOECPotential::copyGridUnrotatedForTest(so_ecp); + testing::TestSOECPotential::copyGridUnrotatedForTest(so_ecp2); + + CHECK(!testing::TestSOECPotential::didGridChange(so_ecp)); + CHECK(!testing::TestSOECPotential::didGridChange(so_ecp2)); + + int num_walkers = 2; + int max_values = 10; + Matrix local_pots(num_walkers, max_values); + Matrix local_pots2(num_walkers, max_values); + std::vector> listeners; + listeners.emplace_back("sopotential", getParticularListener(local_pots)); + listeners.emplace_back("sopotential", getParticularListener(local_pots2)); + + Matrix ion_pots(num_walkers, max_values); + Matrix ion_pots2(num_walkers, max_values); + std::vector> ion_listeners; + ion_listeners.emplace_back("sopotential", getParticularListener(ion_pots)); + ion_listeners.emplace_back("sopotential", getParticularListener(ion_pots2)); + + ParticleSet::mw_update(p_list); + TrialWaveFunction::mw_evaluateLog(twf_list, p_list); + + ListenerOption listener_opt{listeners, ion_listeners}; + testing::TestSOECPotential::mw_evaluateImpl(so_ecp, o_list, twf_list, p_list, listener_opt, true); + + //use single walker API to get reference value + auto value = o_list[0].evaluateDeterministic(p_list[0]); + CHECK(std::accumulate(local_pots.begin(), local_pots.begin() + local_pots.cols(), 0.0) == Approx(value)); + CHECK(std::accumulate(local_pots2.begin(), local_pots2.begin() + local_pots2.cols(), 0.0) == Approx(value)); + CHECK(std::accumulate(ion_pots.begin(), ion_pots.begin() + ion_pots.cols(), 0.0) == Approx(value)); + CHECK(std::accumulate(ion_pots2.begin(), ion_pots2.begin() + ion_pots2.cols(), 0.0) == Approx(value)); + + CHECK(!testing::TestSOECPotential::didGridChange(so_ecp)); + CHECK(!testing::TestSOECPotential::didGridChange(so_ecp2)); + + elec.R[0] = {0.05, 0.0, -0.05}; + elec.update(); + testing::TestSOECPotential::mw_evaluateImpl(so_ecp, o_list, twf_list, p_list, listener_opt, true); + + CHECK(!testing::TestSOECPotential::didGridChange(so_ecp)); + CHECK(!testing::TestSOECPotential::didGridChange(so_ecp2)); + auto value2 = o_list[0].evaluateDeterministic(elec); + + CHECK(std::accumulate(local_pots.begin(), local_pots.begin() + local_pots.cols(), 0.0) == Approx(value2)); + // check the second walker which will be unchanged. + CHECK(std::accumulate(local_pots2[1], local_pots2[1] + local_pots2.cols(), 0.0) == Approx(value)); + + //also check whether or not reference value from single_walker API is actually correct + //this value comes directly from the reference code soecp_eval_reference.cpp + CHECK(value == Approx(-3.530511241)); +} + +TEST_CASE("SOECPotential", "[hamiltonian]") +{ + //do test using VPs. This uses mw_ APIs for TWF in the mw_ SOECP APIs + doSOECPotentialTest(true); + //do test without VPs. This uses legacy APIs for TWF in the mw_ SOECP APIs + doSOECPotentialTest(false); +} + +} // namespace qmcplusplus