diff --git a/src/QMCWaveFunctions/Fermion/MultiDiracDeterminant.2.cpp b/src/QMCWaveFunctions/Fermion/MultiDiracDeterminant.2.cpp index 6b8c4dd0e..81fb7a8de 100644 --- a/src/QMCWaveFunctions/Fermion/MultiDiracDeterminant.2.cpp +++ b/src/QMCWaveFunctions/Fermion/MultiDiracDeterminant.2.cpp @@ -114,7 +114,7 @@ void MultiDiracDeterminant::BuildDotProductsAndCalculateRatios(int ref, { const ValueType det0 = ratios(ref, iat)[dx]; BuildDotProductsAndCalculateRatios_impl(ref, det0, WorkSpace.data(), psiinv, psi, dotProducts, data, pairs, sign); - for (size_t count = 0; count < NumDets; ++count) + for (size_t count = 0; count < getNumDets(); ++count) ratios(count, iat)[dx] = WorkSpace[count]; #if 0 ValueType det0 = ratios(ref,iat)[dx]; @@ -160,7 +160,7 @@ void MultiDiracDeterminant::BuildDotProductsAndCalculateRatios(int ref, const ValueType det0 = ratios(ref, iat); BuildDotProductsAndCalculateRatios_impl(ref, det0, WorkSpace.data(), psiinv, psi, dotProducts, data, pairs, sign); //splatt - for (size_t count = 0; count < NumDets; ++count) + for (size_t count = 0; count < getNumDets(); ++count) ratios(count, iat) = WorkSpace[count]; #if 0 ValueType det0 = ratios(ref,iat); diff --git a/src/QMCWaveFunctions/Fermion/MultiDiracDeterminant.cpp b/src/QMCWaveFunctions/Fermion/MultiDiracDeterminant.cpp index 8c42c5b88..d3713c270 100644 --- a/src/QMCWaveFunctions/Fermion/MultiDiracDeterminant.cpp +++ b/src/QMCWaveFunctions/Fermion/MultiDiracDeterminant.cpp @@ -30,25 +30,18 @@ namespace qmcplusplus { -/** set the index of the first particle in the determinant and reset the size of the determinant - *@param first index of first particle - *@param nel number of particles in the determinant - */ -void MultiDiracDeterminant::set(int first, int nel) +void MultiDiracDeterminant::set(int first, int nel, int ref_det_id) { - APP_ABORT(ClassName + "set(int first, int nel) is disabled. \n"); + assert(ciConfigList); + assert(ciConfigList->size() > 0); + + FirstIndex = first; + ReferenceDeterminant = ref_det_id; + resize(nel); + createDetData((*ciConfigList)[ReferenceDeterminant], *detData, *uniquePairs, *DetSigns); } -void MultiDiracDeterminant::set(int first, int nel, int norb) -{ - FirstIndex = first; - DetCalculator.resize(nel); - resize(nel, norb); - // mmorales; remove later - // testDets(); -} - -void MultiDiracDeterminant::createDetData(ci_configuration2& ref, +void MultiDiracDeterminant::createDetData(const ci_configuration2& ref, std::vector& data, std::vector>& pairs, std::vector& sign) @@ -222,9 +215,6 @@ void MultiDiracDeterminant::copyFromBuffer(ParticleSet& P, WFBufferType& buf) buf.get(FirstAddressOfGrads, LastAddressOfGrads); buf.get(lapls.first_address(), lapls.last_address()); // only used with ORB_PBYP_ALL, - //psiM_temp = psiM; - //dpsiM_temp = dpsiM; - //d2psiM_temp = d2psiM; psiMinv_temp = psiMinv; int n1 = psiM.extent(0); int n2 = psiM.extent(1); @@ -316,24 +306,20 @@ MultiDiracDeterminant::MultiDiracDeterminant(const MultiDiracDeterminant& s) readMatGradTimer(*timer_manager.createTimer(ClassName + "::readMatGrad")), buildTableGradTimer(*timer_manager.createTimer(ClassName + "::buildTableGrad")), ExtraStuffTimer(*timer_manager.createTimer(ClassName + "::ExtraStuff")), + Phi(s.Phi->makeClone()), + NumOrbitals(Phi->getOrbitalSetSize()), NP(0), FirstIndex(s.FirstIndex), - ciConfigList(nullptr) + ciConfigList(s.ciConfigList), + ReferenceDeterminant(s.ReferenceDeterminant), + detData(s.detData), + uniquePairs(s.uniquePairs), + DetSigns(s.DetSigns) { - IsCloned = true; - - ReferenceDeterminant = s.ReferenceDeterminant; - ciConfigList = s.ciConfigList; - NumDets = s.NumDets; - detData = s.detData; - uniquePairs = s.uniquePairs; - DetSigns = s.DetSigns; - Optimizable = s.Optimizable; + Optimizable = s.Optimizable; registerTimers(); - Phi.reset(s.Phi->makeClone()); - this->resize(s.NumPtcls, s.NumOrbitals); - this->DetCalculator.resize(s.NumPtcls); + resize(s.NumPtcls); } SPOSetPtr MultiDiracDeterminant::clonePhi() const { return Phi->makeClone(); } @@ -361,48 +347,24 @@ MultiDiracDeterminant::MultiDiracDeterminant(std::unique_ptr&& spos, int readMatGradTimer(*timer_manager.createTimer(ClassName + "::readMatGrad")), buildTableGradTimer(*timer_manager.createTimer(ClassName + "::buildTableGrad")), ExtraStuffTimer(*timer_manager.createTimer(ClassName + "::ExtraStuff")), + Phi(std::move(spos)), + NumOrbitals(Phi->getOrbitalSetSize()), NP(0), FirstIndex(first), - Phi(std::move(spos)), - ciConfigList(nullptr), ReferenceDeterminant(0) { (Phi->isOptimizable() == true) ? Optimizable = true : Optimizable = false; - IsCloned = false; - - ciConfigList = new std::vector; - detData = new std::vector; - uniquePairs = new std::vector>; - DetSigns = new std::vector; + ciConfigList = std::make_shared>(); + detData = std::make_shared>(); + uniquePairs = std::make_shared>>(); + DetSigns = std::make_shared>(); registerTimers(); } ///default destructor -MultiDiracDeterminant::~MultiDiracDeterminant() {} - -MultiDiracDeterminant& MultiDiracDeterminant::operator=(const MultiDiracDeterminant& s) -{ - if (this == &s) - return *this; - - NP = 0; - IsCloned = true; - ReferenceDeterminant = s.ReferenceDeterminant; - ciConfigList = s.ciConfigList; - NumDets = s.NumDets; - - detData = s.detData; - uniquePairs = s.uniquePairs; - FirstIndex = s.FirstIndex; - DetSigns = s.DetSigns; - - resize(s.NumPtcls, s.NumOrbitals); - this->DetCalculator.resize(s.NumPtcls); - - return *this; -} +MultiDiracDeterminant::~MultiDiracDeterminant() = default; void MultiDiracDeterminant::registerData(ParticleSet& P, WFBufferType& buf) { @@ -412,7 +374,7 @@ void MultiDiracDeterminant::registerData(ParticleSet& P, WFBufferType& buf) //int norb = cols(); NP = P.getTotalNum(); FirstAddressOfGrads = &(grads(0, 0)[0]); - LastAddressOfGrads = FirstAddressOfGrads + NumPtcls * DIM * NumDets; + LastAddressOfGrads = FirstAddressOfGrads + NumPtcls * DIM * getNumDets(); FirstAddressOfdpsiM = &(dpsiM(0, 0)[0]); //(*dpsiM.begin())[0]); LastAddressOfdpsiM = FirstAddressOfdpsiM + NumPtcls * NumOrbitals * DIM; } @@ -428,40 +390,27 @@ void MultiDiracDeterminant::registerData(ParticleSet& P, WFBufferType& buf) } -void MultiDiracDeterminant::setDetInfo(int ref, std::vector* list) -{ - ReferenceDeterminant = ref; - ciConfigList = list; - NumDets = list->size(); -} - ///reset the size: with the number of particles and number of orbtials /// morb is the total number of orbitals, including virtual -void MultiDiracDeterminant::resize(int nel, int morb) +void MultiDiracDeterminant::resize(int nel) { - if (nel <= 0 || morb <= 0) + if (nel <= 0) { APP_ABORT(" ERROR: MultiDiracDeterminant::resize arguments equal to zero. \n"); } - if (NumDets == 0 || NumDets != ciConfigList->size()) - { - APP_ABORT(" ERROR: MultiDiracDeterminant::resize problems with NumDets. \n"); - } - NumPtcls = nel; - NumOrbitals = morb; - LastIndex = FirstIndex + nel; + const int NumDets = getNumDets(); + + NumPtcls = nel; + LastIndex = FirstIndex + nel; psiV_temp.resize(nel); psiV.resize(NumOrbitals); dpsiV.resize(NumOrbitals); d2psiV.resize(NumOrbitals); - psiM.resize(nel, morb); - dpsiM.resize(nel, morb); - d2psiM.resize(nel, morb); - //psiM_temp.resize(nel,morb); - //dpsiM_temp.resize(nel,morb); - //d2psiM_temp.resize(nel,morb); - TpsiM.resize(morb, nel); + psiM.resize(nel, NumOrbitals); + dpsiM.resize(nel, NumOrbitals); + d2psiM.resize(nel, NumOrbitals); + TpsiM.resize(NumOrbitals, nel); psiMinv.resize(nel, nel); dpsiMinv.resize(nel, nel); psiMinv_temp.resize(nel, nel); @@ -476,14 +425,8 @@ void MultiDiracDeterminant::resize(int nel, int morb) new_grads.resize(NumDets, nel); lapls.resize(NumDets, nel); new_lapls.resize(NumDets, nel); - dotProducts.resize(morb, morb); - - //if(ciConfigList==nullptr) - //{ - // APP_ABORT("ciConfigList was not properly initialized.\n"); - //} - if (!IsCloned) - createDetData((*ciConfigList)[ReferenceDeterminant], *detData, *uniquePairs, *DetSigns); + dotProducts.resize(NumOrbitals, NumOrbitals); + DetCalculator.resize(nel); } void MultiDiracDeterminant::registerTimers() diff --git a/src/QMCWaveFunctions/Fermion/MultiDiracDeterminant.h b/src/QMCWaveFunctions/Fermion/MultiDiracDeterminant.h index ecca9b63a..54419ff59 100644 --- a/src/QMCWaveFunctions/Fermion/MultiDiracDeterminant.h +++ b/src/QMCWaveFunctions/Fermion/MultiDiracDeterminant.h @@ -69,7 +69,7 @@ public: */ MultiDiracDeterminant(const MultiDiracDeterminant& s); - MultiDiracDeterminant& operator=(const MultiDiracDeterminant& s); + MultiDiracDeterminant& operator=(const MultiDiracDeterminant& s) = delete; /** return a clone of Phi */ @@ -77,18 +77,14 @@ public: SPOSetPtr getPhi() { return Phi.get(); }; - inline IndexType rows() const { return NumPtcls; } - - inline IndexType cols() const { return NumOrbitals; } - /** set the index of the first particle in the determinant and reset the size of the determinant *@param first index of first particle *@param nel number of particles in the determinant - *@param norb total number of orbitals (including unoccupied) + *@param ref_det_id id of the reference determinant + * + * Note: ciConfigList should have been populated when calling this function */ - void set(int first, int nel, int norb); - - void set(int first, int nel); + void set(int first, int nel, int ref_det_id); void setBF(BackflowTransformation* bf) {} @@ -166,9 +162,6 @@ public: inline void reportStatus(std::ostream& os) {} - ///reset the size: with the number of particles and number of orbtials - virtual void resize(int nel, int morb); - void registerData(ParticleSet& P, WFBufferType& buf); LogValueType updateBuffer(ParticleSet& P, WFBufferType& buf, bool fromscratch = false); @@ -228,7 +221,7 @@ public: // create necessary structures used in the evaluation of the determinants // this works with confgList, which shouldn't change during a simulation - void createDetData(ci_configuration2& ref, + void createDetData(const ci_configuration2& ref, std::vector& data, std::vector>& pairs, std::vector& sign); @@ -370,8 +363,6 @@ public: } */ - void setDetInfo(int ref, std::vector* list); - /** evaluate the value of all the unique determinants with one electron moved. Used by the table method *@param P particle set which provides the positions *@param iat the index of the moved electron @@ -389,43 +380,49 @@ public: // full evaluation of all the structures from scratch, used in evaluateLog for example void evaluateForWalkerMove(const ParticleSet& P, bool fromScratch = true); + // accessors + inline int getNumDets() const { return ciConfigList->size(); } + inline int getNumPtcls() const { return NumPtcls; } + inline int getFirstIndex() const { return FirstIndex; } + inline std::vector& getCIConfigList() { return *ciConfigList; } + + /// store determinants (old and new). FIXME: move to private + ValueVector_t detValues, new_detValues; + GradMatrix_t grads, new_grads; + ValueMatrix_t lapls, new_lapls; + +private: + ///reset the size: with the number of particles + void resize(int nel); + + ///a set of single-particle orbitals used to fill in the values of the matrix + const std::unique_ptr Phi; + ///number of single-particle orbitals which belong to this Dirac determinant + const int NumOrbitals; ///total number of particles int NP; - ///number of single-particle orbitals which belong to this Dirac determinant - int NumOrbitals; ///number of particles which belong to this Dirac determinant int NumPtcls; ///index of the first particle with respect to the particle set int FirstIndex; ///index of the last particle with respect to the particle set int LastIndex; - ///a set of single-particle orbitals used to fill in the values of the matrix - std::unique_ptr Phi; - /// number of determinants handled by this object - int NumDets; - ///bool to cleanup - bool IsCloned; ///use shared_ptr - std::vector* ciConfigList; + std::shared_ptr> ciConfigList; // the reference determinant never changes, so there is no need to store it. // if its value is zero, then use a data from backup, but always use this one // by default int ReferenceDeterminant; - /// store determinants (old and new) - ValueVector_t detValues, new_detValues; - GradMatrix_t grads, new_grads; - ValueMatrix_t lapls, new_lapls; - /// psiM(i,j) \f$= \psi_j({\bf r}_i)\f$ /// TpsiM(i,j) \f$= psiM(j,i) \f$ - ValueMatrix_t psiM, psiM_temp, TpsiM, psiMinv, psiMinv_temp; + ValueMatrix_t psiM, TpsiM, psiMinv, psiMinv_temp; /// dpsiM(i,j) \f$= \nabla_i \psi_j({\bf r}_i)\f$ - GradMatrix_t dpsiM, dpsiM_temp; + GradMatrix_t dpsiM; // temporaty storage ValueMatrix_t dpsiMinv; /// d2psiM(i,j) \f$= \nabla_i^2 \psi_j({\bf r}_i)\f$ - ValueMatrix_t d2psiM, d2psiM_temp; + ValueMatrix_t d2psiM; /// value of single-particle orbital for particle-by-particle update ValueVector_t psiV, psiV_temp; @@ -452,9 +449,9 @@ public: * -i1,i2,...,in : occupied orbital to be replaced (these must be numbers from 0:Nptcl-1) * -a1,a2,...,an : excited states that replace the orbitals (these can be anything) */ - std::vector* detData; - std::vector>* uniquePairs; - std::vector* DetSigns; + std::shared_ptr> detData; + std::shared_ptr>> uniquePairs; + std::shared_ptr> DetSigns; MultiDiracDeterminantCalculator DetCalculator; }; diff --git a/src/QMCWaveFunctions/Fermion/MultiSlaterDeterminant.cpp b/src/QMCWaveFunctions/Fermion/MultiSlaterDeterminant.cpp index caafed177..3d9eccccf 100644 --- a/src/QMCWaveFunctions/Fermion/MultiSlaterDeterminant.cpp +++ b/src/QMCWaveFunctions/Fermion/MultiSlaterDeterminant.cpp @@ -72,35 +72,15 @@ WaveFunctionComponentPtr MultiSlaterDeterminant::makeClone(ParticleSet& tqp) con clone->CSFexpansion = CSFexpansion; clone->DetsPerCSF = DetsPerCSF; } - // SPOSetProxyForMSD* spo = clone->spo_up; - // spo->occup.resize(uniqueConfg_up.size(),clone->nels_up); for (int i = 0; i < dets_up.size(); i++) { - // int nq=0; - // // configuration& ci = uniqueConfg_up[i]; - // for(int k=0; koccup(i,nq++) = k; - // } - // } - SingleDet_t* adet = new SingleDet_t(std::static_pointer_cast(clone->spo_up), 0); - adet->set(clone->FirstIndex_up, clone->nels_up); - clone->dets_up.push_back(adet); + clone->dets_up.push_back(std::make_unique(std::static_pointer_cast(clone->spo_up), 0)); + clone->dets_up.back()->set(clone->FirstIndex_up, clone->nels_up); } - // spo = clone->spo_dn; - // spo->occup.resize(uniqueConfg_dn.size(),clone->nels_dn); for (int i = 0; i < dets_dn.size(); i++) { - // int nq=0; - // // configuration& ci = uniqueConfg_dn[i]; - // for(int k=0; koccup(i,nq++) = k; - // } - // } - SingleDet_t* adet = new SingleDet_t(std::static_pointer_cast(clone->spo_dn), 0); - adet->set(clone->FirstIndex_dn, clone->nels_dn); - clone->dets_dn.push_back(adet); + clone->dets_dn.emplace_back(std::make_unique(std::static_pointer_cast(clone->spo_dn), 0)); + clone->dets_dn.back()->set(clone->FirstIndex_dn, clone->nels_dn); } clone->Optimizable = Optimizable; clone->C = C; diff --git a/src/QMCWaveFunctions/Fermion/MultiSlaterDeterminant.h b/src/QMCWaveFunctions/Fermion/MultiSlaterDeterminant.h index 5598d99d0..d0d24a454 100644 --- a/src/QMCWaveFunctions/Fermion/MultiSlaterDeterminant.h +++ b/src/QMCWaveFunctions/Fermion/MultiSlaterDeterminant.h @@ -54,8 +54,6 @@ public: NewTimer &RatioTimer, &RatioGradTimer, &RatioAllTimer, &UpdateTimer, &EvaluateTimer; NewTimer &Ratio1Timer, &Ratio1GradTimer, &Ratio1AllTimer, &AccRejTimer, &evalOrbTimer; - typedef DiracDeterminantBase* DiracDeterminantBasePtr; - typedef SPOSet* SPOSetPtr; typedef OrbitalSetTraits::IndexVector_t IndexVector_t; typedef OrbitalSetTraits::ValueVector_t ValueVector_t; typedef OrbitalSetTraits::GradVector_t GradVector_t; @@ -129,8 +127,8 @@ public: std::shared_ptr spo_up; std::shared_ptr spo_dn; - std::vector dets_up; - std::vector dets_dn; + std::vector> dets_up; + std::vector> dets_dn; // map determinant in linear combination to unique det list std::vector C2node_up; diff --git a/src/QMCWaveFunctions/Fermion/MultiSlaterDeterminantFast.cpp b/src/QMCWaveFunctions/Fermion/MultiSlaterDeterminantFast.cpp index 8e38a87e5..80e3b0c6d 100644 --- a/src/QMCWaveFunctions/Fermion/MultiSlaterDeterminantFast.cpp +++ b/src/QMCWaveFunctions/Fermion/MultiSlaterDeterminantFast.cpp @@ -174,12 +174,12 @@ WaveFunctionComponent::PsiValueType MultiSlaterDeterminantFast::evaluate_vgl_imp for (size_t ig = 0; ig < Dets.size(); ig++) precomputeC_otherDs(P, ig); - for (size_t i = 0; i < Dets[0]->NumDets; ++i) + for (size_t i = 0; i < Dets[0]->getNumDets(); ++i) psi += C_otherDs[0][i] * Dets[0]->detValues[i]; for (size_t id = 0; id < Dets.size(); id++) - for (size_t i = 0; i < Dets[id]->NumDets; ++i) - for (int k = 0, n = Dets[id]->FirstIndex; k < Dets[id]->NumPtcls; k++, n++) + for (size_t i = 0; i < Dets[id]->getNumDets(); ++i) + for (int k = 0, n = Dets[id]->getFirstIndex(); k < Dets[id]->getNumPtcls(); k++, n++) { g_tmp[n] += C_otherDs[id][i] * Dets[id]->grads(i, k); l_tmp[n] += C_otherDs[id][i] * Dets[id]->lapls(i, k); @@ -229,10 +229,10 @@ WaveFunctionComponent::PsiValueType MultiSlaterDeterminantFast::evalGrad_impl(Pa const GradMatrix_t& grads = (newpos) ? Dets[det_id]->new_grads : Dets[det_id]->grads; const ValueType* restrict detValues0 = (newpos) ? Dets[det_id]->new_detValues.data() : Dets[det_id]->detValues.data(); - const size_t noffset = Dets[det_id]->FirstIndex; + const size_t noffset = Dets[det_id]->getFirstIndex(); PsiValueType psi(0); - for (size_t i = 0; i < Dets[det_id]->NumDets; i++) + for (size_t i = 0; i < Dets[det_id]->getNumDets(); i++) { psi += detValues0[i] * C_otherDs[det_id][i]; g_at += C_otherDs[det_id][i] * grads(i, iat - noffset); @@ -257,7 +257,7 @@ WaveFunctionComponent::PsiValueType MultiSlaterDeterminantFast::evalGrad_impl_no const size_t* restrict det0 = (*C2node)[det_id].data(); const ValueType* restrict cptr = C->data(); const size_t nc = C->size(); - const size_t noffset = Dets[det_id]->FirstIndex; + const size_t noffset = Dets[det_id]->getFirstIndex(); PsiValueType psi(0); for (size_t i = 0; i < nc; ++i) { @@ -324,7 +324,7 @@ WaveFunctionComponent::PsiValueType MultiSlaterDeterminantFast::ratio_impl(Parti // psi=Det_Coeff[i]*Det_Value[unique_det_up]*Det_Value[unique_det_dn]*Det_Value[unique_det_AnyOtherType] // Since only one electron group is moved at the time, identified by det_id, We precompute: // C_otherDs[det_id][i]=Det_Coeff[i]*Det_Value[unique_det_dn]*Det_Value[unique_det_AnyOtherType] - for (size_t i = 0; i < Dets[det_id]->NumDets; i++) + for (size_t i = 0; i < Dets[det_id]->getNumDets(); i++) psi += detValues0[i] * C_otherDs[det_id][i]; return psi; @@ -385,7 +385,7 @@ void MultiSlaterDeterminantFast::evaluateRatios(const VirtualParticleSet& VP, st PsiValueType psiNew(0); if (use_pre_computing_) - for (size_t i = 0; i < Dets[det_id]->NumDets; i++) + for (size_t i = 0; i < Dets[det_id]->getNumDets(); i++) psiNew += detValues0[i] * C_otherDs[det_id][i]; else { @@ -591,7 +591,7 @@ void MultiSlaterDeterminantFast::evaluateDerivatives(ParticleSet& P, for (size_t i = 0; i < laplSum[id].size(); i++) { laplSum[id][i] = 0.0; - for (size_t k = 0; k < Dets[id]->NumPtcls; k++) + for (size_t k = 0; k < Dets[id]->getNumPtcls(); k++) laplSum[id][i] += Dets[id]->lapls[i][k]; } } @@ -599,12 +599,12 @@ void MultiSlaterDeterminantFast::evaluateDerivatives(ParticleSet& P, ValueType lapl_sum = 0.0; myG_temp = 0.0; for (size_t id = 0; id < Dets.size(); id++) - for (size_t i = 0; i < Dets[id]->NumDets; i++) + for (size_t i = 0; i < Dets[id]->getNumDets(); i++) { // assume C_otherDs prepared by evaluateLog already ValueType tmp = C_otherDs[id][i] * psiinv; lapl_sum += tmp * laplSum[id][i]; - for (size_t k = 0, j = Dets[id]->FirstIndex; k < Dets[id]->NumPtcls; k++, j++) + for (size_t k = 0, j = Dets[id]->getFirstIndex(); k < Dets[id]->getNumPtcls(); k++, j++) myG_temp[j] += tmp * Dets[id]->grads(i, k); } @@ -646,7 +646,7 @@ void MultiSlaterDeterminantFast::evaluateDerivatives(ParticleSet& P, tmp *= detValues_otherspin[otherspinC]; } q0 += tmp * laplSum[id][spinC]; - for (size_t l = 0, j = Dets[id]->FirstIndex; l < Dets[id]->NumPtcls; l++, j++) + for (size_t l = 0, j = Dets[id]->getFirstIndex(); l < Dets[id]->getNumPtcls(); l++, j++) v[id] += tmp * static_cast(dot(P.G[j], grads_spin(spinC, l)) - dot(myG_temp[j], grads_spin(spinC, l))); } @@ -681,7 +681,7 @@ void MultiSlaterDeterminantFast::evaluateDerivatives(ParticleSet& P, tmp *= Dets[other_id]->detValues[otherspinC]; } q0 += tmp * laplSum[id][spinC]; - for (size_t l = 0, j = Dets[id]->FirstIndex; l < Dets[id]->NumPtcls; l++, j++) + for (size_t l = 0, j = Dets[id]->getFirstIndex(); l < Dets[id]->getNumPtcls(); l++, j++) v[id] += tmp * static_cast(dot(P.G[j], grads_spin(spinC, l)) - dot(myG_temp[j], grads_spin(spinC, l))); } @@ -833,7 +833,7 @@ void MultiSlaterDeterminantFast::precomputeC_otherDs(const ParticleSet& P, int i // C_otherDs(2, :) stores C x D_up x D_dn ScopedTimer local_timer(PrepareGroupTimer); - C_otherDs[ig].resize(Dets[ig]->NumDets); + C_otherDs[ig].resize(Dets[ig]->getNumDets()); std::fill(C_otherDs[ig].begin(), C_otherDs[ig].end(), ValueType(0)); for (size_t i = 0; i < C->size(); i++) { diff --git a/src/QMCWaveFunctions/Fermion/MultiSlaterDeterminantWithBackflow.cpp b/src/QMCWaveFunctions/Fermion/MultiSlaterDeterminantWithBackflow.cpp index bcd93d3dc..bdac9ebb2 100644 --- a/src/QMCWaveFunctions/Fermion/MultiSlaterDeterminantWithBackflow.cpp +++ b/src/QMCWaveFunctions/Fermion/MultiSlaterDeterminantWithBackflow.cpp @@ -52,17 +52,15 @@ WaveFunctionComponentPtr MultiSlaterDeterminantWithBackflow::makeClone(ParticleS } for (int i = 0; i < dets_up.size(); i++) { - DiracDeterminantWithBackflow* dclne = - (DiracDeterminantWithBackflow*)dets_up[i]->makeCopy(std::static_pointer_cast(clone->spo_up)); - dclne->BFTrans = tr; - clone->dets_up.push_back(dclne); + clone->dets_up.emplace_back(dets_up[i]->makeCopy(std::static_pointer_cast(clone->spo_up))); + auto& dclne = dynamic_cast(*clone->dets_up.back()); + dclne.BFTrans = tr; } for (int i = 0; i < dets_dn.size(); i++) { - DiracDeterminantWithBackflow* dclne = - (DiracDeterminantWithBackflow*)dets_dn[i]->makeCopy(std::static_pointer_cast(clone->spo_dn)); - dclne->BFTrans = tr; - clone->dets_dn.push_back(dclne); + clone->dets_dn.emplace_back(dets_dn[i]->makeCopy(std::static_pointer_cast(clone->spo_dn))); + auto& dclne = dynamic_cast(*clone->dets_dn.back()); + dclne.BFTrans = tr; } clone->Optimizable = Optimizable; clone->C = C; diff --git a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp index 5c08286ca..584ed33ca 100644 --- a/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp +++ b/src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp @@ -594,10 +594,7 @@ bool SlaterDetBuilder::createMSDFast(std::vectorReferenceDeterminant = 0; // for now - Dets[grp]->NumDets = uniqueConfgs[grp].size(); - std::vector& list = *Dets[grp]->ciConfigList; + std::vector& list = Dets[grp]->getCIConfigList(); list.resize(uniqueConfgs[grp].size()); for (int i = 0; i < list.size(); i++) { @@ -612,7 +609,8 @@ bool SlaterDetBuilder::createMSDFast(std::vectorset(targetPtcl.first(grp), nptcls[grp], Dets[grp]->Phi->getOrbitalSetSize()); + // you should choose the det with highest weight for reference. for now choosing 0 + Dets[grp]->set(targetPtcl.first(grp), nptcls[grp], 0); } if (CSFcoeff.size() == 1) @@ -725,9 +723,11 @@ bool SlaterDetBuilder::createMSD(MultiSlaterDeterminant* multiSD, xmlNodePtr cur multiSD->C2node_up = C2nodes[0]; multiSD->C2node_dn = C2nodes[1]; multiSD->resize(uniqueConfgs[0].size(), uniqueConfgs[1].size()); + // alpha dets { auto& spo = multiSD->spo_up; spo->occup.resize(uniqueConfgs[0].size(), multiSD->nels_up); + multiSD->dets_up.reserve(uniqueConfgs[0].size()); for (int i = 0; i < uniqueConfgs[0].size(); i++) { int nq = 0; @@ -749,12 +749,14 @@ bool SlaterDetBuilder::createMSD(MultiSlaterDeterminant* multiSD, xmlNodePtr cur adet = new DiracDeterminant<>(std::static_pointer_cast(spo), 0); } adet->set(multiSD->FirstIndex_up, multiSD->nels_up); - multiSD->dets_up.push_back(adet); + multiSD->dets_up.emplace_back(adet); } } + // beta dets { auto& spo = multiSD->spo_dn; spo->occup.resize(uniqueConfgs[1].size(), multiSD->nels_dn); + multiSD->dets_dn.reserve(uniqueConfgs[1].size()); for (int i = 0; i < uniqueConfgs[1].size(); i++) { int nq = 0; @@ -776,7 +778,7 @@ bool SlaterDetBuilder::createMSD(MultiSlaterDeterminant* multiSD, xmlNodePtr cur adet = new DiracDeterminant<>(std::static_pointer_cast(spo), 0); } adet->set(multiSD->FirstIndex_dn, multiSD->nels_dn); - multiSD->dets_dn.push_back(adet); + multiSD->dets_dn.emplace_back(adet); } } if (multiSD->CSFcoeff.size() == 1 || multiSD->C.size() == 1) diff --git a/src/QMCWaveFunctions/Fermion/ci_configuration2.h b/src/QMCWaveFunctions/Fermion/ci_configuration2.h index cebbe0643..6db1751dd 100644 --- a/src/QMCWaveFunctions/Fermion/ci_configuration2.h +++ b/src/QMCWaveFunctions/Fermion/ci_configuration2.h @@ -18,6 +18,7 @@ #include "CPU/SIMD/simd.hpp" #include #include +#include namespace qmcplusplus { @@ -37,16 +38,11 @@ struct ci_configuration2 inline bool operator==(const ci_configuration2& c) const { if (occup.size() != c.occup.size()) - { - APP_ABORT("ci_configuration2::operator==() - ci_configuration2s are not compatible."); - } + throw std::runtime_error("ci_configuration2::operator==() - ci_configuration2s are not compatible."); + for (int i = 0; i < occup.size(); i++) - { if (occup[i] != c.occup[i]) - { return false; - } - } return true; } @@ -57,15 +53,12 @@ struct ci_configuration2 size_t calculateNumOfExcitations(const ci_configuration2& c) const { if (occup.size() != c.occup.size()) - { - APP_ABORT("ci_configuration2::operator==() - ci_configuration2s are not compatible."); - } + throw std::runtime_error("ci_configuration2::operator==() - ci_configuration2s are not compatible."); + size_t n = 0; for (size_t i = 0; i < occup.size(); i++) - { if (std::find(c.occup.begin(), c.occup.end(), occup[i]) == c.occup.end()) n++; - } return n; } @@ -86,30 +79,27 @@ struct ci_configuration2 std::vector& uno) const { if (occup.size() != c.occup.size()) - { - APP_ABORT("ci_configuration2::operator==() - ci_configuration2s are not compatible."); - } + throw std::runtime_error("ci_configuration2::operator==() - ci_configuration2s are not compatible."); + n = 0; for (size_t i = 0; i < occup.size(); i++) - { if (std::find(c.occup.begin(), c.occup.end(), occup[i]) == c.occup.end()) { pos[n] = i; ocp[n++] = occup[i]; } - } + if (n == 0) return 1.0; + size_t cnt = 0; for (size_t i = 0; i < c.occup.size(); i++) - { if (std::find(occup.begin(), occup.end(), c.occup[i]) == occup.end()) uno[cnt++] = c.occup[i]; - } + if (cnt != n) - { - APP_ABORT(" Error #1 in ci_configuration2::calculateExcitations() \n"); - } + throw std::runtime_error(" Error #1 in ci_configuration2::calculateExcitations() \n"); + double res = 1.0; // this is needed because ci coefficients are given wrt standard ordering, // but by defining the determinant through excitations from a reference might change @@ -118,9 +108,9 @@ struct ci_configuration2 auto ref0(occup); for (size_t i = 0; i < n; i++) ref0[pos[i]] = uno[i]; + for (size_t i = 0; i < ref0.size(); i++) for (size_t k = i + 1; k < ref0.size(); k++) - { if (ref0[i] > ref0[k]) { size_t q = ref0[i]; @@ -129,10 +119,8 @@ struct ci_configuration2 res *= -1.0; } else if (ref0[i] == ref0[k]) - { - APP_ABORT(" Error #2 in ci_configuration2::calculateExcitations() \n"); - } - } + throw std::runtime_error(" Error #2 in ci_configuration2::calculateExcitations() \n"); + return res; } }; diff --git a/src/QMCWaveFunctions/tests/CMakeLists.txt b/src/QMCWaveFunctions/tests/CMakeLists.txt index 76d54cfda..abdb0018b 100644 --- a/src/QMCWaveFunctions/tests/CMakeLists.txt +++ b/src/QMCWaveFunctions/tests/CMakeLists.txt @@ -61,7 +61,7 @@ SET(JASTROW_SRC test_bspline_jastrow.cpp test_counting_jastrow.cpp test_polynomi test_rpa_jastrow.cpp test_user_jastrow.cpp test_kspace_jastrow.cpp test_pade_jastrow.cpp test_short_range_cusp_jastrow.cpp test_DiffTwoBodyJastrowOrbital.cpp) SET(DETERMINANT_SRC FakeSPO.cpp makeRngSpdMatrix.cpp test_DiracDeterminantBatched.cpp test_DiracDeterminantBatched.cpp - test_multi_dirac_determinant.cpp test_dirac_matrix.cpp + test_multi_dirac_determinant.cpp test_dirac_matrix.cpp test_ci_configuration.cpp test_multi_slater_determinant.cpp) IF(ENABLE_CUDA) diff --git a/src/QMCWaveFunctions/tests/test_ci_configuration.cpp b/src/QMCWaveFunctions/tests/test_ci_configuration.cpp new file mode 100644 index 000000000..4226dea11 --- /dev/null +++ b/src/QMCWaveFunctions/tests/test_ci_configuration.cpp @@ -0,0 +1,75 @@ +////////////////////////////////////////////////////////////////////////////////////// +// This file is distributed under the University of Illinois/NCSA Open Source License. +// See LICENSE file in top directory for details. +// +// Copyright (c) 2021 QMCPACK developers. +// +// File developed by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory +// +// File created by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory +////////////////////////////////////////////////////////////////////////////////////// + + +#include "catch.hpp" + +#include "Fermion/ci_configuration2.h" + +namespace qmcplusplus +{ +TEST_CASE("ci_configuration2", "[wavefunction]") +{ + const int n_states = 6; + std::vector ref{0, 1, 2, 3}; + ci_configuration2 ref_state(ref); + + size_t n_excited; + double sign = 0.0; + std::vector pos(n_states); + std::vector occupied(n_states); + std::vector unoccupied(n_states); + + std::vector ext1{0, 1, 4, 5}; + ci_configuration2 ext1_state(ext1); + + n_excited = ext1_state.calculateNumOfExcitations(ref_state); + REQUIRE(n_excited == 2); + + sign = ext1_state.calculateExcitations(ref_state, n_excited, pos, occupied, unoccupied); + + REQUIRE(n_excited == 2); + CHECK(sign == 1.0); + CHECK(pos[0] == 2); + CHECK(pos[1] == 3); + CHECK(occupied[0] == 4); + CHECK(occupied[1] == 5); + CHECK(unoccupied[0] == 2); + CHECK(unoccupied[1] == 3); + + std::vector ext2{0, 1, 2, 6}; + ci_configuration2 ext2_state(ext2); + + n_excited = ext2_state.calculateNumOfExcitations(ref_state); + REQUIRE(n_excited == 1); + + sign = ext2_state.calculateExcitations(ref_state, n_excited, pos, occupied, unoccupied); + + REQUIRE(n_excited == 1); + CHECK(sign == 1.0); + CHECK(pos[0] == 3); + CHECK(occupied[0] == 6); + CHECK(unoccupied[0] == 3); + + std::vector ref2{0, 1}; + ci_configuration2 ref2_state(ref2); + std::vector ext3{1, 6}; + ci_configuration2 ext3_state(ext3); + + sign = ext3_state.calculateExcitations(ref2_state, n_excited, pos, occupied, unoccupied); + REQUIRE(n_excited == 1); + CHECK(sign == -1.0); + + sign = ref2_state.calculateExcitations(ext3_state, n_excited, pos, occupied, unoccupied); + REQUIRE(n_excited == 1); + CHECK(sign == -1.0); +} +} // namespace qmcplusplus