From 6e440ad9c13948352d339f75daff90f48c08eb6f Mon Sep 17 00:00:00 2001 From: "Peter Doak (epd)" Date: Tue, 18 Mar 2025 16:49:04 -0400 Subject: [PATCH 1/8] kcontainer now encapsulated and provides k elements at working prec --- src/Particle/LongRange/EwaldHandler2D.cpp | 10 +- src/Particle/LongRange/EwaldHandler3D.cpp | 14 +- src/Particle/LongRange/EwaldHandler3D.h | 16 +- .../LongRange/EwaldHandlerQuasi2D.cpp | 10 +- src/Particle/LongRange/KContainer.cpp | 134 ++++++++++-- src/Particle/LongRange/KContainer.h | 128 +++++++++--- src/Particle/LongRange/LRBreakup.h | 13 +- src/Particle/LongRange/LRHandlerBase.h | 14 +- src/Particle/LongRange/LRHandlerSRCoulomb.h | 42 ++-- src/Particle/LongRange/LRHandlerTemp.h | 8 +- src/Particle/LongRange/LRRPABFeeHandlerTemp.h | 12 +- src/Particle/LongRange/LRRPAHandlerTemp.h | 8 +- src/Particle/LongRange/StructFact.cpp | 12 +- src/Particle/LongRange/StructFact.h | 2 +- .../LongRange/tests/test_StructFact.cpp | 12 +- src/Particle/LongRange/tests/test_ewald3d.cpp | 6 +- .../LongRange/tests/test_kcontainer.cpp | 195 ++++++++++++++---- .../LongRange/tests/test_lrhandler.cpp | 7 +- src/Particle/LongRange/tests/test_srcoul.cpp | 6 +- src/Particle/LongRange/tests/test_temp.cpp | 3 +- src/QMCHamiltonians/CoulombPBCAA.cpp | 10 +- src/QMCHamiltonians/CoulombPBCAB.cpp | 10 +- src/QMCHamiltonians/SkAllEstimator.cpp | 10 +- src/QMCHamiltonians/SkEstimator.cpp | 8 +- src/QMCHamiltonians/SkPot.cpp | 6 +- src/QMCHamiltonians/SkPot.h | 4 +- src/QMCHamiltonians/StaticStructureFactor.cpp | 4 +- src/QMCHamiltonians/StressPBC.cpp | 4 +- .../tests/test_SkAllEstimator.cpp | 4 +- src/QMCTools/QMCFiniteSize/QMCFiniteSize.cpp | 33 +-- .../ElectronGas/FreeOrbitalBuilder.cpp | 8 +- .../Fermion/BackflowBuilder.cpp | 6 +- .../Fermion/Backflow_ee_kSpace.h | 2 +- src/QMCWaveFunctions/Jastrow/J2KECorrection.h | 2 +- src/QMCWaveFunctions/Jastrow/RPAJastrow.cpp | 4 +- .../Jastrow/RadialJastrowBuilder.cpp | 4 +- src/Utilities/for_testing/checkVector.hpp | 27 +++ .../tests/for_testing/test_checkVector.cpp | 1 + 38 files changed, 559 insertions(+), 240 deletions(-) diff --git a/src/Particle/LongRange/EwaldHandler2D.cpp b/src/Particle/LongRange/EwaldHandler2D.cpp index 667b1f33af..f630f3a166 100644 --- a/src/Particle/LongRange/EwaldHandler2D.cpp +++ b/src/Particle/LongRange/EwaldHandler2D.cpp @@ -34,16 +34,18 @@ void EwaldHandler2D::fillFk(const KContainer& KList) const mRealType kalpha = 1.0 / (2.0*alpha); mRealType kmag, uk; - Fk.resize(KList.kpts_cart.size()); - MaxKshell = KList.kshell.size() - 1; + Fk.resize(KList.getKptsCartWorking().size()); + const auto& kshell = KList.getKShell(); + MaxKshell = kshell.size() - 1; Fk_symm.resize(MaxKshell); + const auto& ksq = KList.getKSQWorking(); for (int ks = 0, ki = 0; ks < Fk_symm.size(); ks++) { - kmag = std::sqrt(KList.ksq[ki]); + kmag = std::sqrt(ksq[ki]); uk = knorm * erfc(kalpha*kmag)/kmag; Fk_symm[ks] = uk; - while (ki < KList.kshell[ks + 1] && ki < Fk.size()) + while (ki < kshell[ks + 1] && ki < Fk.size()) Fk[ki++] = uk; } } diff --git a/src/Particle/LongRange/EwaldHandler3D.cpp b/src/Particle/LongRange/EwaldHandler3D.cpp index d5282321d0..e4f767a578 100644 --- a/src/Particle/LongRange/EwaldHandler3D.cpp +++ b/src/Particle/LongRange/EwaldHandler3D.cpp @@ -53,9 +53,10 @@ EwaldHandler3D::EwaldHandler3D(const EwaldHandler3D& aLR, ParticleSet& ref) void EwaldHandler3D::fillFk(const KContainer& KList) { - Fk.resize(KList.kpts_cart.size()); - Fkg.resize(KList.kpts_cart.size()); - const std::vector& kshell(KList.kshell); + const auto& kpts_cart = KList.getKptsCartWorking(); + Fk.resize(kpts_cart.size()); + Fkg.resize(kpts_cart.size()); + const std::vector& kshell(KList.getKShell()); MaxKshell = kshell.size() - 1; @@ -63,12 +64,13 @@ void EwaldHandler3D::fillFk(const KContainer& KList) kMag.resize(MaxKshell); mRealType kgauss = 1.0 / (4 * Sigma * Sigma); mRealType knorm = 4 * M_PI / Volume; + const auto ksq = KList.getKSQWorking(); for (int ks = 0, ki = 0; ks < Fk_symm.size(); ks++) { - mRealType t2e = KList.ksq[ki] * kgauss; - mRealType uk = knorm * std::exp(-t2e) / KList.ksq[ki]; + mRealType t2e = ksq[ki] * kgauss; + mRealType uk = knorm * std::exp(-t2e) / ksq[ki]; Fk_symm[ks] = uk; - while (ki < KList.kshell[ks + 1] && ki < Fk.size()) + while (ki < kshell[ks + 1] && ki < Fk.size()) Fk[ki++] = uk; } diff --git a/src/Particle/LongRange/EwaldHandler3D.h b/src/Particle/LongRange/EwaldHandler3D.h index 4c49e6fa4b..2b8fce566e 100644 --- a/src/Particle/LongRange/EwaldHandler3D.h +++ b/src/Particle/LongRange/EwaldHandler3D.h @@ -103,24 +103,26 @@ class EwaldHandler3D : public LRHandlerBase void fillYkgstrain(const KContainer& KList) { - Fkgstrain.resize(KList.kpts_cart.size()); - const std::vector& kshell(KList.kshell); + Fkgstrain.resize(KList.getKptsCartWorking().size()); + const std::vector& kshell(KList.getKShell()); MaxKshell = kshell.size() - 1; + const auto& ksq = KList.getKSQWorking(); for (int ks = 0, ki = 0; ks < MaxKshell; ks++) { - mRealType uk = evalYkgstrain(std::sqrt(KList.ksq[ki])); - while (ki < KList.kshell[ks + 1] && ki < Fkgstrain.size()) + mRealType uk = evalYkgstrain(std::sqrt(ksq[ki])); + while (ki < kshell[ks + 1] && ki < Fkgstrain.size()) Fkgstrain[ki++] = uk; } } void filldFk_dk(const KContainer& KList) { - dFk_dstrain.resize(KList.kpts_cart.size()); - + const auto& kpts_cart = KList.getKptsCartWorking(); + dFk_dstrain.resize(kpts_cart.size()); + const auto& ksq = KList.getKSQWorking(); for (int ki = 0; ki < dFk_dstrain.size(); ki++) { - dFk_dstrain[ki] = evaluateLR_dstrain(KList.kpts_cart[ki], std::sqrt(KList.ksq[ki])); + dFk_dstrain[ki] = evaluateLR_dstrain(kpts_cart[ki], std::sqrt(ksq[ki])); } } diff --git a/src/Particle/LongRange/EwaldHandlerQuasi2D.cpp b/src/Particle/LongRange/EwaldHandlerQuasi2D.cpp index 599c5a853a..0d234d5b2e 100644 --- a/src/Particle/LongRange/EwaldHandlerQuasi2D.cpp +++ b/src/Particle/LongRange/EwaldHandlerQuasi2D.cpp @@ -35,19 +35,19 @@ void EwaldHandlerQuasi2D::fillFk(const KContainer& KList) { const mRealType knorm = M_PI / area; mRealType kmag, uk; - - Fk.resize(KList.kpts_cart.size()); - MaxKshell = KList.kshell.size() - 1; + const auto& kpts_cart = KList.getKptsCartWorking(); + Fk.resize(kpts_cart.size()); + MaxKshell = KList.getKShell().size() - 1; Fk_symm.resize(MaxKshell); kmags.resize(MaxKshell); for (int ks = 0, ki = 0; ks < Fk_symm.size(); ks++) { - kmag = std::sqrt(KList.ksq[ki]); + kmag = std::sqrt(KList.getKSQWorking()[ki]); kmags[ks] = kmag; // store k magnitutes uk = knorm/kmag; Fk_symm[ks] = uk; - while (ki < KList.kshell[ks + 1] && ki < Fk.size()) + while (ki < KList.getKShell()[ks + 1] && ki < Fk.size()) Fk[ki++] = uk; } } diff --git a/src/Particle/LongRange/KContainer.cpp b/src/Particle/LongRange/KContainer.cpp index 72d4c8bd17..f7f8d4fe57 100644 --- a/src/Particle/LongRange/KContainer.cpp +++ b/src/Particle/LongRange/KContainer.cpp @@ -2,11 +2,12 @@ // This file is distributed under the University of Illinois/NCSA Open Source License. // See LICENSE file in top directory for details. // -// Copyright (c) 2016 Jeongnim Kim and QMCPACK developers. +// Copyright (c) 2025 QMCPACK developers. // // File developed by: Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign // Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign // Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory +// Peter W. Doak, doakpw@ornl.gov, Oak Ridge National Laboratory // // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign ////////////////////////////////////////////////////////////////////////////////////// @@ -21,11 +22,39 @@ namespace qmcplusplus { -void KContainer::updateKLists(const ParticleLayout& lattice, - RealType kc, - unsigned ndim, - const PosType& twist, - bool useSphere) + +template +const std::vector::AppPosition>& KContainerT::getKptsCartWorking() const +{ + if constexpr (std::is_same_v) + return kpts_cart_; + else + return kpts_cart_working_; +} + +template +const std::vector& KContainerT::getKSQWorking() const +{ + if constexpr (std::is_same::value) + return ksq_; + else + return ksq_working_; +} + +template +int KContainerT::getMinusK(int k) const +{ + assert(k < minusk.size()); + return minusk[k]; +} + +template +template +void KContainerT::updateKLists(const CrystalLattice& lattice, + FullPrecReal kc, + unsigned ndim, + const Position& twist, + bool useSphere) { kcutoff = kc; if (kcutoff <= 0.0) @@ -37,11 +66,13 @@ void KContainer::updateKLists(const ParticleLayout& lattice, app_log() << " KContainer initialised with cutoff " << kcutoff << std::endl; app_log() << " # of K-shell = " << kshell.size() << std::endl; - app_log() << " # of K points = " << kpts.size() << std::endl; + app_log() << " # of K points = " << kpts_.size() << std::endl; app_log() << std::endl; } -void KContainer::findApproxMMax(const ParticleLayout& lattice, unsigned ndim) +template +template +void KContainerT::findApproxMMax(const CrystalLattice& lattice, unsigned ndim) { //Estimate the size of the parallelpiped that encompasses a sphere of kcutoff. //mmax is stored as integer translations of the reciprocal cell vectors. @@ -102,19 +133,23 @@ void KContainer::findApproxMMax(const ParticleLayout& lattice, unsigned ndim) mmax[1] = 0; } -void KContainer::BuildKLists(const ParticleLayout& lattice, const PosType& twist, bool useSphere) +template +template +void KContainerT::BuildKLists(const CrystalLattice& lattice, + const Position& twist, + bool useSphere) { TinyVector TempActualMax; TinyVector kvec; - TinyVector kvec_cart; - RealType modk2; + TinyVector kvec_cart; + FullPrecReal modk2; std::vector> kpts_tmp; - std::vector kpts_cart_tmp; - std::vector ksq_tmp; + std::vector kpts_cart_tmp; + std::vector ksq_tmp; // reserve the space for memory efficiency if (useSphere) { - const RealType kcut2 = kcutoff * kcutoff; + const FullPrecReal kcut2 = kcutoff * kcutoff; //Loop over guesses for valid k-points. for (int i = -mmax[0]; i <= mmax[0]; i++) { @@ -208,10 +243,10 @@ void KContainer::BuildKLists(const ParticleLayout& lattice, const PosType& twist } } std::map*>::iterator it(kpts_sorted.begin()); - kpts.resize(numk); - kpts_cart.resize(numk); + kpts_.resize(numk); + kpts_cart_.resize(numk); kpts_cart_soa_.resize(numk); - ksq.resize(numk); + ksq_.resize(numk); kshell.resize(kpts_sorted.size() + 1, 0); int ok = 0, ish = 0; while (it != kpts_sorted.end()) @@ -220,10 +255,10 @@ void KContainer::BuildKLists(const ParticleLayout& lattice, const PosType& twist while (vit != (*it).second->end()) { int ik = (*vit); - kpts[ok] = kpts_tmp[ik]; - kpts_cart[ok] = kpts_cart_tmp[ik]; + kpts_[ok] = kpts_tmp[ik]; + kpts_cart_[ok] = kpts_cart_tmp[ik]; kpts_cart_soa_(ok) = kpts_cart_tmp[ik]; - ksq[ok] = ksq_tmp[ik]; + ksq_[ok] = ksq_tmp[ik]; ++vit; ++ok; } @@ -232,6 +267,11 @@ void KContainer::BuildKLists(const ParticleLayout& lattice, const PosType& twist ++ish; } kpts_cart_soa_.updateTo(); + if constexpr (!std::is_same::value) + { + std::copy(kpts_cart_.begin(), kpts_cart_.end(), std::back_inserter(kpts_cart_working_)); + std::copy(ksq_.begin(), ksq_.end(), std::back_inserter(ksq_working_)); + } it = kpts_sorted.begin(); std::map*>::iterator e_it(kpts_sorted.end()); while (it != e_it) @@ -262,13 +302,63 @@ void KContainer::BuildKLists(const ParticleLayout& lattice, const PosType& twist std::map hashToIndex; for (int ki = 0; ki < numk; ki++) { - hashToIndex[getHashOfVec(kpts[ki], numk)] = ki; + hashToIndex[getHashOfVec(kpts_[ki], numk)] = ki; } // Use the map to find the index of -k from the index of k for (int ki = 0; ki < numk; ki++) { - minusk[ki] = hashToIndex[getHashOfVec(-1 * kpts[ki], numk)]; + minusk[ki] = hashToIndex[getHashOfVec(-1 * kpts_[ki], numk)]; } } +#ifdef MIXED_PRECISION +template class KContainerT; +template void KContainerT::updateKLists(const CrystalLattice& lattice, + FullPrecReal kc, + unsigned ndim, + const Position& twist = Position(), + bool useSphere = true); +template void KContainerT::updateKLists(const CrystalLattice& lattice, + FullPrecReal kc, + unsigned ndim, + const Position& twist = Position(), + bool useSphere = true); +template void KContainerT::findApproxMMax(const CrystalLattice& lattice, unsigned ndim); +template void KContainerT::findApproxMMax(const CrystalLattice& lattice, + unsigned ndim); + + +template class KContainerT; +template void KContainerT::updateKLists(const CrystalLattice& lattice, + FullPrecReal kc, + unsigned ndim, + const Position& twist = Position(), + bool useSphere = true); +template void KContainerT::updateKLists(const CrystalLattice& lattice, + FullPrecReal kc, + unsigned ndim, + const Position& twist = Position(), + bool useSphere = true); +template void KContainerT::findApproxMMax(const CrystalLattice& lattice, + unsigned ndim); +template void KContainerT::findApproxMMax(const CrystalLattice& lattice, + unsigned ndim); +#else +template class KContainerT; +template void KContainerT::updateKLists(const CrystalLattice& lattice, + FullPrecReal kc, + unsigned ndim, + const Position& twist = Position(), + bool useSphere = true); +template void KContainerT::updateKLists(const CrystalLattice& lattice, + FullPrecReal kc, + unsigned ndim, + const Position& twist = Position(), + bool useSphere = true); +template void KContainerT::findApproxMMax(const CrystalLattice& lattice, + unsigned ndim); +template void KContainerT::findApproxMMax(const CrystalLattice& lattice, + unsigned ndim); + +#endif } // namespace qmcplusplus diff --git a/src/Particle/LongRange/KContainer.h b/src/Particle/LongRange/KContainer.h index 3f1f0fad63..650f7baa69 100644 --- a/src/Particle/LongRange/KContainer.h +++ b/src/Particle/LongRange/KContainer.h @@ -2,11 +2,12 @@ // This file is distributed under the University of Illinois/NCSA Open Source License. // See LICENSE file in top directory for details. // -// Copyright (c) 2016 Jeongnim Kim and QMCPACK developers. +// Copyright (c) 2025 QMCPACK developers. // // File developed by: Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign // Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign // Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory +// Peter W. Doak, doakpw@ornl.gov, Oak Ridge National Laboratory // // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign ////////////////////////////////////////////////////////////////////////////////////// @@ -26,16 +27,52 @@ namespace qmcplusplus * It generates a set of k-points that are unit-translations of the reciprocal-space * cell. K-points are generated within a spherical cutoff set by the supercell */ -class KContainer : public QMCTraits +template +class KContainerT { +public: + using Real = REAL; + using FullPrecReal = QMCTraits::FullPrecRealType; + using PositionFull = QMCTraits::QTFull::PosType; + static constexpr int DIM = OHMMS_DIM; + using AppPosition = typename QMCTraits::PosType; + using Position = typename QMCTypes::PosType; + private: /// The cutoff up to which k-vectors are generated. - RealType kcutoff; + FullPrecReal kcutoff; public: //Typedef for the lattice-type - using ParticleLayout = Lattice; + using ParticleLayout = CrystalLattice; + + const auto& get_kpts_cart_soa() const { return kpts_cart_soa_; } + + const std::vector>& getKpts() const { return kpts_; } + const std::vector& getKptsCartWorking() const; + + const std::vector& getKShell() const { return kshell; }; + + const std::vector& getKSQWorking() const; + + int getMinusK(int k) const; + + int getNumK() const { return numk; } + + /** update k-vectors + * @param sc supercell + * @param kc cutoff radius in the K + * @param twist shifts the center of the grid of k-vectors + * @param useSphere if true, use the |K| + */ + template + void updateKLists(const CrystalLattice& lattice, + FullPrecReal kc, + unsigned ndim, + const Position& twist = Position(), + bool useSphere = true); +private: ///number of k-points int numk; @@ -47,13 +84,16 @@ class KContainer : public QMCTraits /** K-vector in reduced coordinates */ - std::vector> kpts; + std::vector> kpts_; /** K-vector in Cartesian coordinates */ - std::vector kpts_cart; + std::vector kpts_cart_; + std::vector kpts_cart_working_; /** squre of kpts in Cartesian coordniates */ - std::vector ksq; + std::vector ksq_; + std::vector ksq_working_; + /** Given a k index, return index to -k */ std::vector minusk; @@ -67,33 +107,73 @@ class KContainer : public QMCTraits */ //std::map*> kpts_sorted; - /** update k-vectors - * @param sc supercell - * @param kc cutoff radius in the K - * @param twist shifts the center of the grid of k-vectors - * @param useSphere if true, use the |K| - */ - void updateKLists(const ParticleLayout& lattice, - RealType kc, - unsigned ndim, - const PosType& twist = PosType(), - bool useSphere = true); - - const auto& get_kpts_cart_soa() const { return kpts_cart_soa_; } - private: /** compute approximate parallelpiped that surrounds kc * @param lattice supercell */ - void findApproxMMax(const ParticleLayout& lattice, unsigned ndim); + template + void findApproxMMax(const CrystalLattice& lattice, unsigned ndim); /** construct the container for k-vectors */ - void BuildKLists(const ParticleLayout& lattice, const PosType& twist, bool useSphere); + template + void BuildKLists(const CrystalLattice& lattice, const Position& twist, bool useSphere); /** K-vector in Cartesian coordinates in SoA layout */ - VectorSoaContainer> kpts_cart_soa_; + VectorSoaContainer> kpts_cart_soa_; }; +using KContainer = KContainerT; + +#ifdef MIXED_PRECISION +extern template class KContainerT; +extern template void KContainerT::updateKLists(const CrystalLattice& lattice, + FullPrecReal kc, + unsigned ndim, + const Position& twist = Position(), + bool useSphere = true); +extern template void KContainerT::updateKLists(const CrystalLattice& lattice, + FullPrecReal kc, + unsigned ndim, + const Position& twist = Position(), + bool useSphere = true); +extern template void KContainerT::findApproxMMax(const CrystalLattice& lattice, + unsigned ndim); +extern template void KContainerT::findApproxMMax(const CrystalLattice& lattice, + unsigned ndim); + +extern template class KContainerT; +extern template void KContainerT::updateKLists(const CrystalLattice& lattice, + FullPrecReal kc, + unsigned ndim, + const Position& twist = Position(), + bool useSphere = true); +extern template void KContainerT::updateKLists(const CrystalLattice& lattice, + FullPrecReal kc, + unsigned ndim, + const Position& twist = Position(), + bool useSphere = true); +extern template void KContainerT::findApproxMMax(const CrystalLattice& lattice, + unsigned ndim); +extern template void KContainerT::findApproxMMax(const CrystalLattice& lattice, + unsigned ndim); +#else +extern template class KContainerT; +extern template void KContainerT::updateKLists(const CrystalLattice& lattice, + FullPrecReal kc, + unsigned ndim, + const Position& twist = Position(), + bool useSphere = true); +extern template void KContainerT::updateKLists(const CrystalLattice& lattice, + FullPrecReal kc, + unsigned ndim, + const Position& twist = Position(), + bool useSphere = true); +extern template void KContainerT::findApproxMMax(const CrystalLattice& lattice, + unsigned ndim); +extern template void KContainerT::findApproxMMax(const CrystalLattice& lattice, + unsigned ndim); +#endif + } // namespace qmcplusplus #endif diff --git a/src/Particle/LongRange/LRBreakup.h b/src/Particle/LongRange/LRBreakup.h index 98c5c2c6c0..e9b6bb906c 100644 --- a/src/Particle/LongRange/LRBreakup.h +++ b/src/Particle/LongRange/LRBreakup.h @@ -334,23 +334,24 @@ int LRBreakup::SetupKVecs(mRealType kc, mRealType kcont, mRealType { //Add low |k| ( < kcont) k-points with exact degeneracy KContainer kexact; - kexact.updateKLists(Basis.get_Lattice(), kcont, Basis.get_Lattice().ndim); + const auto& basis_lattice = Basis.get_Lattice(); + kexact.updateKLists>::Real>(basis_lattice, kcont, basis_lattice.ndim); bool findK = true; mRealType kc2 = kc * kc; //use at least one shell size_t ks = 0; - kc2 = std::max(kc2, static_cast(kexact.ksq[kexact.kshell[ks]])); + kc2 = std::max(kc2, static_cast(kexact.getKSQWorking()[kexact.getKShell()[ks]])); while (findK) { - if (kexact.ksq[kexact.kshell[ks]] > kc2) + if (kexact.getKSQWorking()[kexact.getKShell()[ks]] > kc2) findK = false; else ks++; } size_t maxkshell = ks; - size_t numk = kexact.numk - kexact.kshell[ks]; - for (; ks < kexact.kshell.size() - 1; ks++) - AddKToList(std::sqrt(kexact.ksq[kexact.kshell[ks]]), kexact.kshell[ks + 1] - kexact.kshell[ks]); + size_t numk = kexact.getNumK() - kexact.getKShell()[ks]; + for (; ks < kexact.getKShell().size() - 1; ks++) + AddKToList(std::sqrt(kexact.getKSQWorking()[kexact.getKShell()[ks]]), kexact.getKShell()[ks + 1] - kexact.getKShell()[ks]); ////Add these vectors to the internal list //int numk=0; //mRealType modk2; diff --git a/src/Particle/LongRange/LRHandlerBase.h b/src/Particle/LongRange/LRHandlerBase.h index 40fb15456c..a5518dc650 100644 --- a/src/Particle/LongRange/LRHandlerBase.h +++ b/src/Particle/LongRange/LRHandlerBase.h @@ -136,7 +136,7 @@ struct LRHandlerBase const Matrix& e2ikrA_i = A.getSK().eikr_i; const pRealType* rhokB_r = B.getSK().rhok_r[specB]; const pRealType* rhokB_i = B.getSK().rhok_i[specB]; - const std::vector& kpts = A.getSimulationCell().getKLists().kpts_cart; + const std::vector& kpts = A.getSimulationCell().getKLists().getKptsCartWorking(); for (int ki = 0; ki < Fk.size(); ki++) { PosType k = kpts[ki]; @@ -224,23 +224,23 @@ struct DummyLRHandler : public LRHandlerBase mRealType norm = 4.0 * M_PI / ref.getLattice().Volume; mRealType kcsq = LR_kc * LR_kc; auto& KList(ref.getSimulationCell().getKLists()); - int maxshell = KList.kshell.size() - 1; - const auto& kk(KList.ksq); + int maxshell = KList.getKShell().size() - 1; + const auto& kk(KList.getKSQWorking()); int ksh = 0, ik = 0; while (ksh < maxshell) { if (kk[ik] > kcsq) break; //exit - ik = KList.kshell[++ksh]; + ik = KList.getKShell()[++ksh]; } MaxKshell = ksh; Fk_symm.resize(MaxKshell); - Fk.resize(KList.kpts_cart.size()); + Fk.resize(KList.getKptsCartWorking().size()); for (ksh = 0, ik = 0; ksh < MaxKshell; ksh++, ik++) { - mRealType v = norm * myFunc(kk[KList.kshell[ksh]]); //rpa=u0/kk[ik]; + mRealType v = norm * myFunc(kk[KList.getKShell()[ksh]]); //rpa=u0/kk[ik]; Fk_symm[ksh] = v; - for (; ik < KList.kshell[ksh + 1]; ik++) + for (; ik < KList.getKShell()[ksh + 1]; ik++) Fk[ik] = v; } } diff --git a/src/Particle/LongRange/LRHandlerSRCoulomb.h b/src/Particle/LongRange/LRHandlerSRCoulomb.h index 46c224fc0d..f34e70a4d4 100644 --- a/src/Particle/LongRange/LRHandlerSRCoulomb.h +++ b/src/Particle/LongRange/LRHandlerSRCoulomb.h @@ -399,40 +399,43 @@ class LRHandlerSRCoulomb : public LRHandlerBase void fillYk(KContainer& KList) { - Fk.resize(KList.kpts_cart.size()); - const std::vector& kshell(KList.kshell); + Fk.resize(KList.getKptsCartWorking().size()); + const std::vector& kshell(KList.getKShell()); if (MaxKshell >= kshell.size()) MaxKshell = kshell.size() - 1; Fk_symm.resize(MaxKshell); + const auto& ksq = KList.getKSQWorking(); for (int ks = 0, ki = 0; ks < Fk_symm.size(); ks++) { - mRealType uk = evalYk(std::sqrt(KList.ksq[ki])); + mRealType uk = evalYk(std::sqrt(ksq[ki])); Fk_symm[ks] = uk; - while (ki < KList.kshell[ks + 1] && ki < Fk.size()) + while (ki < kshell[ks + 1] && ki < Fk.size()) Fk[ki++] = uk; } - //for(int ki=0; ki& kshell(KList.kshell); + Fk.resize(kpts_cart.size()); + const std::vector& kshell(KList.getKShell()); if (MaxKshell >= kshell.size()) MaxKshell = kshell.size() - 1; + const auto& ksq = KList.getKSQWorking(); for (int ks = 0, ki = 0; ks < MaxKshell; ks++) { - mRealType uk = evalYkg(std::sqrt(KList.ksq[ki])); + mRealType uk = evalYkg(std::sqrt(ksq[ki])); - while (ki < KList.kshell[ks + 1] && ki < Fkg.size()) + while (ki < kshell[ks + 1] && ki < Fkg.size()) Fkg[ki++] = uk; } //Have to set this, because evaluate and evaluateGrad for LR piece uses @@ -444,24 +447,25 @@ class LRHandlerSRCoulomb : public LRHandlerBase void fillYkgstrain(KContainer& KList) { APP_ABORT("Stresses not supported yet\n"); - Fkgstrain.resize(KList.kpts_cart.size()); - const std::vector& kshell(KList.kshell); + Fkgstrain.resize(KList.getKptsCartWorking().size()); + const std::vector& kshell(KList.getKShell()); if (MaxKshell >= kshell.size()) MaxKshell = kshell.size() - 1; + const auto& ksq = KList.getKSQWorking(); for (int ks = 0, ki = 0; ks < MaxKshell; ks++) { - mRealType uk = evalYkgstrain(std::sqrt(KList.ksq[ki])); - while (ki < KList.kshell[ks + 1] && ki < Fkgstrain.size()) + mRealType uk = evalYkgstrain(std::sqrt(ksq[ki])); + while (ki < kshell[ks + 1] && ki < Fkgstrain.size()) Fkgstrain[ki++] = uk; } } void filldFk_dk(KContainer& KList) { - APP_ABORT("Stresses not supported yet\n"); - dFk_dstrain.resize(KList.kpts_cart.size()); + throw std::runtime_error("Stresses not supported yet\n"); + // dFk_dstrain.resize(KList.getKptsCartWorking().size()); - for (int ki = 0; ki < dFk_dstrain.size(); ki++) - dFk_dstrain[ki] = evaluateLR_dstrain(KList.kpts_cart[ki], std::sqrt(KList.ksq[ki])); + // for (int ki = 0; ki < dFk_dstrain.size(); ki++) + // dFk_dstrain[ki] = evaluateLR_dstrain(KList.getKptsCartWorking()[ki], std::sqrt(KList.getKSQWorking()[ki])); } }; } // namespace qmcplusplus diff --git a/src/Particle/LongRange/LRHandlerTemp.h b/src/Particle/LongRange/LRHandlerTemp.h index 1cfe699ad7..aba37b2881 100644 --- a/src/Particle/LongRange/LRHandlerTemp.h +++ b/src/Particle/LongRange/LRHandlerTemp.h @@ -267,16 +267,16 @@ class LRHandlerTemp : public LRHandlerBase void fillFk(const KContainer& KList) { - Fk.resize(KList.kpts_cart.size()); - const std::vector& kshell(KList.kshell); + Fk.resize(KList.getKptsCartWorking().size()); + const std::vector& kshell(KList.getKShell()); if (MaxKshell >= kshell.size()) MaxKshell = kshell.size() - 1; Fk_symm.resize(MaxKshell); for (int ks = 0, ki = 0; ks < Fk_symm.size(); ks++) { - mRealType uk = evalFk(std::sqrt(KList.ksq[ki])); + mRealType uk = evalFk(std::sqrt(KList.getKSQWorking()[ki])); Fk_symm[ks] = uk; - while (ki < KList.kshell[ks + 1] && ki < Fk.size()) + while (ki < kshell[ks + 1] && ki < Fk.size()) Fk[ki++] = uk; } //for(int ki=0; ki& kshell(KList.kshell); + Fk.resize(KList.getKptsCartWorking().size()); + const std::vector& kshell(KList.getKShell()); if (MaxKshell >= kshell.size()) MaxKshell = kshell.size() - 1; Fk_symm.resize(MaxKshell); // std::cout<<"Filling FK :"<& kshell(KList.kshell); + Fk.resize(KList.getKptsCartWorking().size()); + const std::vector& kshell(KList.getKShell()); if (MaxKshell >= kshell.size()) MaxKshell = kshell.size() - 1; Fk_symm.resize(MaxKshell); // std::cout<<"Filling FK :"<& sk_list const size_t nw = p_list.size(); const size_t num_species = p_leader.groups(); const auto& kpts_cart = sk_leader.k_lists_.get_kpts_cart_soa(); - const size_t nk = sk_leader.k_lists_.numk; + const size_t nk = sk_leader.k_lists_.getNumK(); const size_t nk_padded = kpts_cart.capacity(); auto& coordinates_leader = static_cast(p_leader.getCoordinates()); @@ -139,11 +139,11 @@ void StructFact::computeRhok(const ParticleSet& P) { const size_t num_ptcls = P.getTotalNum(); const size_t num_species = P.groups(); - const size_t nk = k_lists_.numk; + const size_t nk = k_lists_.getNumK(); resize(nk, num_species, num_ptcls); - rhok_r = 0.0; rhok_i = 0.0; + const auto& kpts_cart = k_lists_.getKptsCartWorking(); if (StorePerParticle) { // save per particle and species value @@ -157,7 +157,7 @@ void StructFact::computeRhok(const ParticleSet& P) #pragma omp simd for (int ki = 0; ki < nk; ki++) { - qmcplusplus::sincos(dot(k_lists_.kpts_cart[ki], pos), &eikr_i_ptr[ki], &eikr_r_ptr[ki]); + qmcplusplus::sincos(dot(kpts_cart[ki], pos), &eikr_i_ptr[ki], &eikr_r_ptr[ki]); rhok_r_ptr[ki] += eikr_r_ptr[ki]; rhok_i_ptr[ki] += eikr_i_ptr[ki]; } @@ -176,7 +176,7 @@ void StructFact::computeRhok(const ParticleSet& P) for (int ki = 0; ki < nk; ki++) { RealType s, c; - qmcplusplus::sincos(dot(k_lists_.kpts_cart[ki], pos), &s, &c); + qmcplusplus::sincos(dot(kpts_cart[ki], pos), &s, &c); rhok_r_ptr[ki] += c; rhok_i_ptr[ki] += s; } @@ -191,7 +191,7 @@ void StructFact::computeRhok(const ParticleSet& P) const size_t offset = ib * kblock_size; const size_t this_block_size = std::min(kblock_size, nk - offset); for (int ki = 0; ki < this_block_size; ki++) - phiV[ki] = dot(k_lists_.kpts_cart[ki + offset], pos); + phiV[ki] = dot(kpts_cart[ki + offset], pos); eval_e2iphi(this_block_size, phiV, eikr_r_temp, eikr_i_temp); for (int ki = 0; ki < this_block_size; ki++) { diff --git a/src/Particle/LongRange/StructFact.h b/src/Particle/LongRange/StructFact.h index 615f096a71..ade58fbd12 100644 --- a/src/Particle/LongRange/StructFact.h +++ b/src/Particle/LongRange/StructFact.h @@ -21,11 +21,11 @@ #include #include #include +#include "KContainer.h" namespace qmcplusplus { class ParticleSet; -class KContainer; struct SKMultiWalkerMem; /** @ingroup longrange diff --git a/src/Particle/LongRange/tests/test_StructFact.cpp b/src/Particle/LongRange/tests/test_StructFact.cpp index a7a53c7d02..b5f93c57f0 100644 --- a/src/Particle/LongRange/tests/test_StructFact.cpp +++ b/src/Particle/LongRange/tests/test_StructFact.cpp @@ -48,22 +48,22 @@ TEST_CASE("StructFact", "[lrhandler]") ref.R[2] = {0.3, 4.0, 1.4}; ref.R[3] = {3.2, 4.7, 0.7}; - REQUIRE(simulation_cell.getKLists().numk == 263786); + REQUIRE(simulation_cell.getKLists().getNumK() == 263786); StructFact sk(ref.getLRBox(), simulation_cell.getKLists()); sk.updateAllPart(ref); CHECK(sk.rhok_r.rows() == ref.groups()); CHECK(sk.rhok_i.rows() == ref.groups()); - CHECK(sk.rhok_r.cols() == simulation_cell.getKLists().numk); - CHECK(sk.rhok_i.cols() == simulation_cell.getKLists().numk); + CHECK(sk.rhok_r.cols() == simulation_cell.getKLists().getNumK()); + CHECK(sk.rhok_i.cols() == simulation_cell.getKLists().getNumK()); std::vector> rhok_sum_ref{-125.80618630936, 68.199075127271}; for (int i = 0; i < ref.groups(); i++) { - std::complex rhok_sum, rhok_even_sum; - for (int ik = 0; ik < simulation_cell.getKLists().numk; ik++) - rhok_sum += std::complex(sk.rhok_r[i][ik], sk.rhok_i[i][ik]); + std::complex rhok_sum, rhok_even_sum; + for (int ik = 0; ik < simulation_cell.getKLists().getNumK(); ik++) + rhok_sum += std::complex(sk.rhok_r[i][ik], sk.rhok_i[i][ik]); //std::cout << std::setprecision(14) << rhok_sum << std::endl; CHECK(ComplexApprox(rhok_sum).epsilon(5e-5) == rhok_sum_ref[i]); diff --git a/src/Particle/LongRange/tests/test_ewald3d.cpp b/src/Particle/LongRange/tests/test_ewald3d.cpp index 13b9de25ee..92d807406a 100644 --- a/src/Particle/LongRange/tests/test_ewald3d.cpp +++ b/src/Particle/LongRange/tests/test_ewald3d.cpp @@ -46,8 +46,7 @@ TEST_CASE("ewald3d", "[lrhandler]") CHECK(handler.Sigma == Approx(std::sqrt(lattice.LR_kc / (2.0 * lattice.LR_rc)))); std::cout << "handler.MaxKshell is " << handler.MaxKshell << std::endl; - CHECK( (std::is_same::value ? - handler.MaxKshell == 78 : handler.MaxKshell >= 117 && handler.MaxKshell <= 128 )); + CHECK(handler.MaxKshell == 78); CHECK(handler.LR_rc == Approx(2.5)); CHECK(handler.LR_kc == Approx(12)); @@ -95,8 +94,7 @@ TEST_CASE("ewald3d df", "[lrhandler]") CHECK(handler.Sigma == Approx(std::sqrt(lattice.LR_kc / (2.0 * lattice.LR_rc)))); std::cout << "handler.MaxKshell is " << handler.MaxKshell << std::endl; - CHECK( (std::is_same::value ? - handler.MaxKshell == 78 : handler.MaxKshell >= 117 && handler.MaxKshell <= 128 )); + CHECK(handler.MaxKshell == 78); CHECK(handler.LR_rc == Approx(2.5)); CHECK(handler.LR_kc == Approx(12)); diff --git a/src/Particle/LongRange/tests/test_kcontainer.cpp b/src/Particle/LongRange/tests/test_kcontainer.cpp index b15922e500..8057fe022d 100644 --- a/src/Particle/LongRange/tests/test_kcontainer.cpp +++ b/src/Particle/LongRange/tests/test_kcontainer.cpp @@ -2,9 +2,10 @@ // This file is distributed under the University of Illinois/NCSA Open Source License. // See LICENSE file in top directory for details. // -// Copyright (c) 2022 QMCPACK developers. +// Copyright (c) 2025 QMCPACK developers. // // File developed by: Yubo "Paul" Yang, yubo.paul.yang@gmail.com, CCQ @ Flatiron +// Peter W. Doak, doakpw@ornl.gov, Oak Ridge National Lab // // File created by: Yubo "Paul" Yang, yubo.paul.yang@gmail.com, CCQ @ Flatiron ////////////////////////////////////////////////////////////////////////////////////// @@ -12,6 +13,10 @@ #include "Particle/LongRange/KContainer.h" #include "OhmmsPETE/TinyVector.h" +#include "Particle/Lattice/CrystalLattice.h" +#include "Particle/tests/MinimalParticlePool.h" +#include "Utilities/for_testing/checkVector.hpp" +#include "Utilities/for_testing/NativeInitializerPrint.hpp" namespace qmcplusplus { @@ -27,7 +32,7 @@ TEST_CASE("kcontainer at gamma in 3D", "[longrange]") const std::vector kcs = {blat, std::sqrt(2) * blat, std::sqrt(3) * blat}; const std::vector nks = {6, 18, 26}; - Lattice lattice; + CrystalLattice lattice; lattice.R.diagonal(1.0); lattice.set(lattice.R); // compute Rv and Gv from R @@ -35,10 +40,10 @@ TEST_CASE("kcontainer at gamma in 3D", "[longrange]") for (int ik = 0; ik < kcs.size(); ik++) { const double kc = kcs[ik] + 1e-6; - klists.updateKLists(lattice, kc, ndim); - CHECK(klists.kpts.size() == nks[ik]); + klists.updateKLists>::Scalar_t>(lattice, kc, ndim); + CHECK(klists.getKpts().size() == nks[ik]); } - const int mxk = klists.kpts.size(); + const int mxk = klists.getKpts().size(); int gvecs[26][3] = {{-1, 0, 0}, {0, -1, 0}, {0, 0, -1}, {0, 0, 1}, {0, 1, 0}, {1, 0, 0}, {-1, -1, 0}, {-1, 0, -1}, {-1, 0, 1}, {-1, 1, 0}, {0, -1, -1}, {0, -1, 1}, {0, 1, -1}, {0, 1, 1}, {1, -1, 0}, {1, 0, -1}, {1, 0, 1}, {1, 1, 0}, {-1, -1, -1}, {-1, -1, 1}, {-1, 1, -1}, @@ -48,8 +53,8 @@ TEST_CASE("kcontainer at gamma in 3D", "[longrange]") { for (int ldim = 0; ldim < ndim; ldim++) { - CHECK(klists.kpts[ik][ldim] == gvecs[ik][ldim]); - CHECK(klists.kpts[ik][ldim] * blat == Approx(klists.kpts_cart[ik][ldim])); + CHECK(klists.getKpts()[ik][ldim] == gvecs[ik][ldim]); + CHECK(klists.getKpts()[ik][ldim] * blat == Approx(klists.getKptsCartWorking()[ik][ldim])); } } } @@ -63,7 +68,7 @@ TEST_CASE("kcontainer at twist in 3D", "[longrange]") // twist one shell of kvectors const double kc = blat + 1e-6; - Lattice lattice; + CrystalLattice lattice; lattice.R.diagonal(1.0); lattice.set(lattice.R); // compute Rv and Gv from R @@ -71,17 +76,17 @@ TEST_CASE("kcontainer at twist in 3D", "[longrange]") PosType twist; twist[0] = 0.1; - klists.updateKLists(lattice, kc, ndim, twist); - CHECK(klists.kpts.size() == 1); - CHECK(klists.kpts_cart[0][0] == Approx(blat * (twist[0] - 1))); + klists.updateKLists>::Scalar_t>(lattice, kc, ndim, twist); + CHECK(klists.getKpts().size() == 1); + CHECK(klists.getKptsCartWorking()[0][0] == Approx(blat * (twist[0] - 1))); twist = {-0.5, 0, 0.5}; - klists.updateKLists(lattice, kc, ndim, twist); + klists.updateKLists>::Scalar_t>(lattice, kc, ndim, twist); int gvecs[3][3] = {{0, 0, -1}, {1, 0, -1}, {1, 0, 0}}; - CHECK(klists.kpts.size() == 3); - for (int ik = 0; ik < klists.kpts.size(); ik++) + CHECK(klists.getKpts().size() == 3); + for (int ik = 0; ik < klists.getKpts().size(); ik++) for (int ldim = 0; ldim < ndim; ldim++) - CHECK(klists.kpts_cart[ik][ldim] == Approx(blat * (twist[ldim] + gvecs[ik][ldim]))); + CHECK(klists.getKptsCartWorking()[ik][ldim] == Approx(blat * (twist[ldim] + gvecs[ik][ldim]))); } TEST_CASE("kcontainer at gamma in 2D", "[longrange]") @@ -94,7 +99,7 @@ TEST_CASE("kcontainer at gamma in 2D", "[longrange]") const std::vector kcs = {blat, std::sqrt(2) * blat, 2 * blat}; const std::vector nks = {4, 8, 12}; - Lattice lattice; + CrystalLattice lattice; lattice.R.diagonal(1.0); lattice.set(lattice.R); // compute Rv and Gv from R @@ -102,10 +107,10 @@ TEST_CASE("kcontainer at gamma in 2D", "[longrange]") for (int ik = 0; ik < kcs.size(); ik++) { const double kc = kcs[ik] + 1e-6; - klists.updateKLists(lattice, kc, ndim); - CHECK(klists.kpts.size() == nks[ik]); + klists.updateKLists>::Scalar_t>(lattice, kc, ndim); + CHECK(klists.getKpts().size() == nks[ik]); } - const int mxk = klists.kpts.size(); + const int mxk = klists.getKpts().size(); int gvecs[12][3] = { {-1, 0, 0}, {0, -1, 0}, {0, 1, 0}, {1, 0, 0}, {-1, -1, 0}, {-1, 1, 0}, {1, -1, 0}, {1, 1, 0}, {-2, 0, 0}, {0, -2, 0}, {0, 2, 0}, {2, 0, 0}, @@ -115,8 +120,8 @@ TEST_CASE("kcontainer at gamma in 2D", "[longrange]") { for (int ldim = 0; ldim < ndim; ldim++) { - CHECK(klists.kpts[ik][ldim] == gvecs[ik][ldim]); - CHECK(klists.kpts[ik][ldim] * blat == Approx(klists.kpts_cart[ik][ldim])); + CHECK(klists.getKpts()[ik][ldim] == gvecs[ik][ldim]); + CHECK(klists.getKpts()[ik][ldim] * blat == Approx(klists.getKptsCartWorking()[ik][ldim])); } } } @@ -130,7 +135,7 @@ TEST_CASE("kcontainer at twist in 2D", "[longrange]") // twist one shell of kvectors const double kc = blat + 1e-6; - Lattice lattice; + CrystalLattice lattice; lattice.R.diagonal(1.0); lattice.set(lattice.R); // compute Rv and Gv from R @@ -138,29 +143,139 @@ TEST_CASE("kcontainer at twist in 2D", "[longrange]") PosType twist; twist[0] = 0.1; - klists.updateKLists(lattice, kc, ndim, twist); - CHECK(klists.kpts.size() == 1); - CHECK(klists.kpts_cart[0][0] == Approx(blat * (twist[0] - 1))); + klists.updateKLists>::Scalar_t>(lattice, kc, ndim, twist); + CHECK(klists.getKpts().size() == 1); + CHECK(klists.getKptsCartWorking()[0][0] == Approx(blat * (twist[0] - 1))); twist[1] = 0.1; - klists.updateKLists(lattice, kc, ndim, twist); - CHECK(klists.kpts.size() == 2); - CHECK(klists.kpts_cart[0][0] == Approx(blat * (twist[0] - 1))); - CHECK(klists.kpts_cart[0][1] == Approx(blat * twist[1])); - CHECK(klists.kpts_cart[1][0] == Approx(blat * (twist[0]))); - CHECK(klists.kpts_cart[1][1] == Approx(blat * (twist[1] - 1))); + klists.updateKLists>::Scalar_t>(lattice, kc, ndim, twist); + CHECK(klists.getKpts().size() == 2); + CHECK(klists.getKptsCartWorking()[0][0] == Approx(blat * (twist[0] - 1))); + CHECK(klists.getKptsCartWorking()[0][1] == Approx(blat * twist[1])); + CHECK(klists.getKptsCartWorking()[1][0] == Approx(blat * (twist[0]))); + CHECK(klists.getKptsCartWorking()[1][1] == Approx(blat * (twist[1] - 1))); twist = {-0.5, 0.5, 0}; - klists.updateKLists(lattice, kc, ndim, twist); - CHECK(klists.kpts.size() == 3); + klists.updateKLists>::Scalar_t>(lattice, kc, ndim, twist); + CHECK(klists.getKpts().size() == 3); //for (int ik=0;ik<3;ik++) - // app_log() << klists.kpts_cart[ik] << std::endl; - CHECK(klists.kpts_cart[0][0] == Approx(blat * (twist[0] - 0))); - CHECK(klists.kpts_cart[0][1] == Approx(blat * (twist[1] - 1))); - CHECK(klists.kpts_cart[1][0] == Approx(blat * (twist[0] + 1))); - CHECK(klists.kpts_cart[1][1] == Approx(blat * (twist[1] - 1))); - CHECK(klists.kpts_cart[2][0] == Approx(blat * (twist[0] + 1))); - CHECK(klists.kpts_cart[2][1] == Approx(blat * twist[1])); + // app_log() << klists.getKptsCartWorking()[ik] << std::endl; + CHECK(klists.getKptsCartWorking()[0][0] == Approx(blat * (twist[0] - 0))); + CHECK(klists.getKptsCartWorking()[0][1] == Approx(blat * (twist[1] - 1))); + CHECK(klists.getKptsCartWorking()[1][0] == Approx(blat * (twist[0] + 1))); + CHECK(klists.getKptsCartWorking()[1][1] == Approx(blat * (twist[1] - 1))); + CHECK(klists.getKptsCartWorking()[2][0] == Approx(blat * (twist[0] + 1))); + CHECK(klists.getKptsCartWorking()[2][1] == Approx(blat * twist[1])); } +TEST_CASE("kcontainer for diamond", "[longrange]") +{ + // Communicate* comm = OHMMS::Controller; + // ParticleSetPool particle_pool{MinimalParticlePool::make_diamondC_1x1x1(comm)}; + // ParticleSet pset_ions{*(particle_pool.getParticleSet("ion"))}; + + int ndim = 3; + + // auto kcs = pset_ions.getSimulationCell().getKLists().numk; + // std::cout << "kcs:" << kcs << '\n'; + + using Real = QMCTraits::RealType; + using FullPrecReal = QMCTraits::FullPrecRealType; + + CrystalLattice lattice; + lattice.BoxBConds = {true, true, true}; + Tensor lattice_mat{3.37316115, 3.37316115, 0.00000000, 0.00000000, 3.37316115, + 3.37316115, 3.37316115, 0.00000000, 3.37316115}; + lattice.set(lattice_mat); + lattice.LR_dim_cutoff = 15; + SimulationCell cell(lattice); + + const KContainerT& klists = cell.getKLists(); + std::remove_cv_t> kpoint_lists = { + {-1, -1, -1}, {-1, 0, 0}, {0, -1, 0}, {0, 0, -1}, {0, 0, 1}, {0, 1, 0}, {1, 0, 0}, {1, 1, 1}, + {-1, -1, 0}, {-1, 0, -1}, {0, -1, -1}, {0, 1, 1}, {1, 0, 1}, {1, 1, 0}, {-2, -1, -1}, {-1, -2, -1}, + {-1, -1, -2}, {-1, 0, 1}, {-1, 1, 0}, {0, -1, 1}, {0, 1, -1}, {1, -1, 0}, {1, 0, -1}, {1, 1, 2}, + {1, 2, 1}, {2, 1, 1}, {-2, -2, -1}, {-2, -1, -2}, {-2, -1, 0}, {-2, 0, -1}, {-1, -2, -2}, {-1, -2, 0}, + {-1, -1, 1}, {-1, 0, -2}, {-1, 1, -1}, {-1, 1, 1}, {0, -2, -1}, {0, -1, -2}, {0, 1, 2}, {0, 2, 1}, + {1, -1, -1}, {1, -1, 1}, {1, 0, 2}, {1, 1, -1}, {1, 2, 0}, {1, 2, 2}, {2, 0, 1}, {2, 1, 0}, + {2, 1, 2}, {2, 2, 1}, {-2, -2, -2}, {-2, 0, 0}, {0, -2, 0}, {0, 0, -2}, {0, 0, 2}, {0, 2, 0}, + {2, 0, 0}, {2, 2, 2}, {-2, -2, 0}, {-2, 0, -2}, {0, -2, -2}, {0, 2, 2}, {2, 0, 2}, {2, 2, 0}, + {-3, -2, -2}, {-3, -1, -1}, {-2, -3, -2}, {-2, -2, -3}, {-2, 0, 1}, {-2, 1, 0}, {-1, -3, -1}, {-1, -1, -3}, + {-1, 0, 2}, {-1, 2, 0}, {0, -2, 1}, {0, -1, 2}, {0, 1, -2}, {0, 2, -1}, {1, -2, 0}, {1, 0, -2}, + {1, 1, 3}, {1, 3, 1}, {2, -1, 0}, {2, 0, -1}, {2, 2, 3}, {2, 3, 2}, {3, 1, 1}, {3, 2, 2}, + {-3, -2, -1}, {-3, -1, -2}, {-2, -3, -1}, {-2, -1, -3}, {-2, -1, 1}, {-2, 1, -1}, {-1, -3, -2}, {-1, -2, -3}, + {-1, -2, 1}, {-1, 1, -2}, {-1, 1, 2}, {-1, 2, 1}, {1, -2, -1}, {1, -1, -2}, {1, -1, 2}, {1, 2, -1}, + {1, 2, 3}, {1, 3, 2}, {2, -1, 1}, {2, 1, -1}, {2, 1, 3}, {2, 3, 1}, {3, 1, 2}, {3, 2, 1}, + {-3, -3, -2}, {-3, -2, -3}, {-3, -1, 0}, {-3, 0, -1}, {-2, -3, -3}, {-2, 1, 1}, {-1, -3, 0}, {-1, -1, 2}, + {-1, 0, -3}, {-1, 2, -1}, {0, -3, -1}, {0, -1, -3}, {0, 1, 3}, {0, 3, 1}, {1, -2, 1}, {1, 0, 3}, + {1, 1, -2}, {1, 3, 0}, {2, -1, -1}, {2, 3, 3}, {3, 0, 1}, {3, 1, 0}, {3, 2, 3}, {3, 3, 2}, + {-3, -3, -3}, {-3, -3, -1}, {-3, -2, 0}, {-3, -1, -3}, {-3, 0, -2}, {-3, 0, 0}, {-2, -3, 0}, {-2, -2, 1}, + {-2, 0, -3}, {-2, 1, -2}, {-1, -3, -3}, {-1, 2, 2}, {0, -3, -2}, {0, -3, 0}, {0, -2, -3}, {0, 0, -3}, + {0, 0, 3}, {0, 2, 3}, {0, 3, 0}, {0, 3, 2}, {1, -2, -2}, {1, 3, 3}, {2, -1, 2}, {2, 0, 3}, + {2, 2, -1}, {2, 3, 0}, {3, 0, 0}, {3, 0, 2}, {3, 1, 3}, {3, 2, 0}, {3, 3, 1}, {3, 3, 3}, + {-4, -2, -2}, {-2, -4, -2}, {-2, -2, -4}, {-2, 0, 2}, {-2, 2, 0}, {0, -2, 2}, {0, 2, -2}, {2, -2, 0}, + {2, 0, -2}, {2, 2, 4}, {2, 4, 2}, {4, 2, 2}, {-4, -3, -2}, {-4, -2, -3}, {-4, -2, -1}, {-4, -1, -2}, + {-3, -4, -2}, {-3, -2, -4}, {-3, -1, 1}, {-3, 1, -1}, {-2, -4, -3}, {-2, -4, -1}, {-2, -3, -4}, {-2, -1, -4}, + {-2, -1, 2}, {-2, 1, 2}, {-2, 2, -1}, {-2, 2, 1}, {-1, -4, -2}, {-1, -3, 1}, {-1, -2, -4}, {-1, -2, 2}, + {-1, 1, -3}, {-1, 1, 3}, {-1, 2, -2}, {-1, 3, 1}, {1, -3, -1}, {1, -2, 2}, {1, -1, -3}, {1, -1, 3}, + {1, 2, -2}, {1, 2, 4}, {1, 3, -1}, {1, 4, 2}, {2, -2, -1}, {2, -2, 1}, {2, -1, -2}, {2, 1, -2}, + {2, 1, 4}, {2, 3, 4}, {2, 4, 1}, {2, 4, 3}, {3, -1, 1}, {3, 1, -1}, {3, 2, 4}, {3, 4, 2}, + {4, 1, 2}, {4, 2, 1}, {4, 2, 3}, {4, 3, 2}, {-4, -3, -3}, {-4, -1, -1}, {-3, -4, -3}, {-3, -3, -4}, + {-3, -3, 0}, {-3, 0, -3}, {-3, 0, 1}, {-3, 1, 0}, {-1, -4, -1}, {-1, -1, -4}, {-1, 0, 3}, {-1, 3, 0}, + {0, -3, -3}, {0, -3, 1}, {0, -1, 3}, {0, 1, -3}, {0, 3, -1}, {0, 3, 3}, {1, -3, 0}, {1, 0, -3}, + {1, 1, 4}, {1, 4, 1}, {3, -1, 0}, {3, 0, -1}, {3, 0, 3}, {3, 3, 0}, {3, 3, 4}, {3, 4, 3}, + {4, 1, 1}, {4, 3, 3}, {-4, -3, -1}, {-4, -1, -3}, {-3, -4, -1}, {-3, -2, 1}, {-3, -1, -4}, {-3, 1, -2}, + {-2, -3, 1}, {-2, 1, -3}, {-1, -4, -3}, {-1, -3, -4}, {-1, 2, 3}, {-1, 3, 2}, {1, -3, -2}, {1, -2, -3}, + {1, 3, 4}, {1, 4, 3}, {2, -1, 3}, {2, 3, -1}, {3, -1, 2}, {3, 1, 4}, {3, 2, -1}, {3, 4, 1}, + {4, 1, 3}, {4, 3, 1}, {-4, -4, -3}, {-4, -3, -4}, {-4, -1, 0}, {-4, 0, -1}, {-3, -4, -4}, {-3, 1, 1}, + {-1, -4, 0}, {-1, -1, 3}, {-1, 0, -4}, {-1, 3, -1}, {0, -4, -1}, {0, -1, -4}, {0, 1, 4}, {0, 4, 1}, + {1, -3, 1}, {1, 0, 4}, {1, 1, -3}, {1, 4, 0}, {3, -1, -1}, {3, 4, 4}, {4, 0, 1}, {4, 1, 0}, + {4, 3, 4}, {4, 4, 3}, {-4, -4, -2}, {-4, -2, -4}, {-4, -2, 0}, {-4, 0, -2}, {-2, -4, -4}, {-2, -4, 0}, + {-2, -2, 2}, {-2, 0, -4}, {-2, 2, -2}, {-2, 2, 2}, {0, -4, -2}, {0, -2, -4}, {0, 2, 4}, {0, 4, 2}, + {2, -2, -2}, {2, -2, 2}, {2, 0, 4}, {2, 2, -2}, {2, 4, 0}, {2, 4, 4}, {4, 0, 2}, {4, 2, 0}, + {4, 2, 4}, {4, 4, 2}, {-4, -4, -4}, {-4, 0, 0}, {0, -4, 0}, {0, 0, -4}, {0, 0, 4}, {0, 4, 0}, + {4, 0, 0}, {4, 4, 4}, {-5, -3, -3}, {-5, -2, -2}, {-4, -4, -1}, {-4, -3, 0}, {-4, -1, -4}, {-4, 0, -3}, + {-3, -5, -3}, {-3, -4, 0}, {-3, -3, -5}, {-3, -3, 1}, {-3, 0, -4}, {-3, 0, 2}, {-3, 1, -3}, {-3, 2, 0}, + {-2, -5, -2}, {-2, -2, -5}, {-2, 0, 3}, {-2, 3, 0}, {-1, -4, -4}, {-1, 3, 3}, {0, -4, -3}, {0, -3, -4}, + {0, -3, 2}, {0, -2, 3}, {0, 2, -3}, {0, 3, -2}, {0, 3, 4}, {0, 4, 3}, {1, -3, -3}, {1, 4, 4}, + {2, -3, 0}, {2, 0, -3}, {2, 2, 5}, {2, 5, 2}, {3, -2, 0}, {3, -1, 3}, {3, 0, -2}, {3, 0, 4}, + {3, 3, -1}, {3, 3, 5}, {3, 4, 0}, {3, 5, 3}, {4, 0, 3}, {4, 1, 4}, {4, 3, 0}, {4, 4, 1}, + {5, 2, 2}, {5, 3, 3}, {-5, -3, -2}, {-5, -2, -3}, {-3, -5, -2}, {-3, -2, -5}, {-3, -1, 2}, {-3, 2, -1}, + {-2, -5, -3}, {-2, -3, -5}, {-2, 1, 3}, {-2, 3, 1}, {-1, -3, 2}, {-1, 2, -3}, {1, -2, 3}, {1, 3, -2}, + {2, -3, -1}, {2, -1, -3}, {2, 3, 5}, {2, 5, 3}, {3, -2, 1}, {3, 1, -2}, {3, 2, 5}, {3, 5, 2}, + {5, 2, 3}, {5, 3, 2}, {-5, -4, -3}, {-5, -3, -4}, {-5, -2, -1}, {-5, -1, -2}, {-4, -5, -3}, {-4, -3, -5}, + {-4, -1, 1}, {-4, 1, -1}, {-3, -5, -4}, {-3, -4, -5}, {-3, 1, 2}, {-3, 2, 1}, {-2, -5, -1}, {-2, -1, -5}, + {-2, -1, 3}, {-2, 3, -1}, {-1, -5, -2}, {-1, -4, 1}, {-1, -2, -5}, {-1, -2, 3}, {-1, 1, -4}, {-1, 1, 4}, + {-1, 3, -2}, {-1, 4, 1}, {1, -4, -1}, {1, -3, 2}, {1, -1, -4}, {1, -1, 4}, {1, 2, -3}, {1, 2, 5}, + {1, 4, -1}, {1, 5, 2}, {2, -3, 1}, {2, 1, -3}, {2, 1, 5}, {2, 5, 1}, {3, -2, -1}, {3, -1, -2}, + {3, 4, 5}, {3, 5, 4}, {4, -1, 1}, {4, 1, -1}, {4, 3, 5}, {4, 5, 3}, {5, 1, 2}, {5, 2, 1}, + {5, 3, 4}, {5, 4, 3}, {-5, -4, -4}, {-5, -4, -2}, {-5, -3, -1}, {-5, -2, -4}, {-5, -1, -3}, {-5, -1, -1}, + {-4, -5, -4}, {-4, -5, -2}, {-4, -4, -5}, {-4, -2, -5}, {-4, -2, 1}, {-4, 0, 1}, {-4, 1, -2}, {-4, 1, 0}, + {-3, -5, -1}, {-3, -2, 2}, {-3, -1, -5}, {-3, 2, -2}, {-2, -5, -4}, {-2, -4, -5}, {-2, -4, 1}, {-2, -3, 2}, + {-2, 1, -4}, {-2, 2, -3}, {-2, 2, 3}, {-2, 3, 2}, {-1, -5, -3}, {-1, -5, -1}, {-1, -3, -5}, {-1, -1, -5}, + {-1, 0, 4}, {-1, 2, 4}, {-1, 4, 0}, {-1, 4, 2}, {0, -4, 1}, {0, -1, 4}, {0, 1, -4}, {0, 4, -1}, + {1, -4, -2}, {1, -4, 0}, {1, -2, -4}, {1, 0, -4}, {1, 1, 5}, {1, 3, 5}, {1, 5, 1}, {1, 5, 3}, + {2, -3, -2}, {2, -2, -3}, {2, -2, 3}, {2, -1, 4}, {2, 3, -2}, {2, 4, -1}, {2, 4, 5}, {2, 5, 4}, + {3, -2, 2}, {3, 1, 5}, {3, 2, -2}, {3, 5, 1}, {4, -1, 0}, {4, -1, 2}, {4, 0, -1}, {4, 2, -1}, + {4, 2, 5}, {4, 4, 5}, {4, 5, 2}, {4, 5, 4}, {5, 1, 1}, {5, 1, 3}, {5, 2, 4}, {5, 3, 1}, + {5, 4, 2}, {5, 4, 4}, {-4, -4, 0}, {-4, 0, -4}, {0, -4, -4}, {0, 4, 4}, {4, 0, 4}, {4, 4, 0}, + {-5, -5, -3}, {-5, -3, -5}, {-5, -2, 0}, {-5, 0, -2}, {-3, -5, -5}, {-3, 2, 2}, {-2, -5, 0}, {-2, -2, 3}, + {-2, 0, -5}, {-2, 3, -2}, {0, -5, -2}, {0, -2, -5}, {0, 2, 5}, {0, 5, 2}, {2, -3, 2}, {2, 0, 5}, + {2, 2, -3}, {2, 5, 0}, {3, -2, -2}, {3, 5, 5}, {5, 0, 2}, {5, 2, 0}, {5, 3, 5}, {5, 5, 3}, + {-5, -5, -4}, {-5, -4, -5}, {-5, -4, -1}, {-5, -1, -4}, {-5, -1, 0}, {-5, 0, -1}, {-4, -5, -5}, {-4, -5, -1}, + {-4, -3, 1}, {-4, -1, -5}, {-4, 1, -3}, {-4, 1, 1}, {-3, -4, 1}, {-3, 1, -4}, {-1, -5, -4}, {-1, -5, 0}, + {-1, -4, -5}, {-1, -1, 4}, {-1, 0, -5}, {-1, 3, 4}, {-1, 4, -1}, {-1, 4, 3}, {0, -5, -1}, {0, -1, -5}, + {0, 1, 5}, {0, 5, 1}, {1, -4, -3}, {1, -4, 1}, {1, -3, -4}, {1, 0, 5}, {1, 1, -4}, {1, 4, 5}, + {1, 5, 0}, {1, 5, 4}, {3, -1, 4}, {3, 4, -1}, {4, -1, -1}, {4, -1, 3}, {4, 1, 5}, {4, 3, -1}, + {4, 5, 1}, {4, 5, 5}, {5, 0, 1}, {5, 1, 0}, {5, 1, 4}, {5, 4, 1}, {5, 4, 5}, {5, 5, 4}, + }; + double tolerance = 0.1; + { + INFO("Checking kpoint_lists"); + auto check = checkVector(klists.getKpts(), kpoint_lists, true); + CHECKED_ELSE(check.result) { FAIL(check.result_message); } + } +} + + } // namespace qmcplusplus diff --git a/src/Particle/LongRange/tests/test_lrhandler.cpp b/src/Particle/LongRange/tests/test_lrhandler.cpp index 5fe47def87..b74ddef148 100644 --- a/src/Particle/LongRange/tests/test_lrhandler.cpp +++ b/src/Particle/LongRange/tests/test_lrhandler.cpp @@ -48,8 +48,7 @@ TEST_CASE("dummy", "[lrhandler]") handler.initBreakup(ref); std::cout << "handler.MaxKshell is " << handler.MaxKshell << std::endl; - CHECK( (std::is_same::value ? - handler.MaxKshell == 78 : handler.MaxKshell >= 117 && handler.MaxKshell <= 128 )); + CHECK( handler.MaxKshell == 78); CHECK(handler.LR_kc == Approx(12)); CHECK(handler.LR_rc == Approx(0)); @@ -61,8 +60,8 @@ TEST_CASE("dummy", "[lrhandler]") // the full Coulomb potential should be retained in kspace for (int ish = 0; ish < handler.MaxKshell; ish++) { - int ik = ref.getSimulationCell().getKLists().kshell[ish]; - double k2 = ref.getSimulationCell().getKLists().ksq[ik]; + int ik = ref.getSimulationCell().getKLists().getKShell()[ish]; + double k2 = ref.getSimulationCell().getKLists().getKSQWorking()[ik]; double fk_expect = fk(k2); CHECK(handler.Fk_symm[ish] == Approx(norm * fk_expect)); } diff --git a/src/Particle/LongRange/tests/test_srcoul.cpp b/src/Particle/LongRange/tests/test_srcoul.cpp index 5f44773cda..6be2f080dd 100644 --- a/src/Particle/LongRange/tests/test_srcoul.cpp +++ b/src/Particle/LongRange/tests/test_srcoul.cpp @@ -54,8 +54,7 @@ TEST_CASE("srcoul", "[lrhandler]") handler.initBreakup(ref); std::cout << "handler.MaxKshell is " << handler.MaxKshell << std::endl; - CHECK( (std::is_same::value ? - handler.MaxKshell == 78 : handler.MaxKshell >= 117 && handler.MaxKshell <= 128 )); + CHECK( handler.MaxKshell == 78); CHECK(Approx(handler.LR_rc) == 2.5); CHECK(Approx(handler.LR_kc) == 12); @@ -99,8 +98,7 @@ TEST_CASE("srcoul df", "[lrhandler]") handler.initBreakup(ref); std::cout << "handler.MaxKshell is " << handler.MaxKshell << std::endl; - CHECK( (std::is_same::value ? - handler.MaxKshell == 78 : handler.MaxKshell >= 117 && handler.MaxKshell <= 128 )); + CHECK( handler.MaxKshell == 78); CHECK(Approx(handler.LR_rc) == 2.5); CHECK(Approx(handler.LR_kc) == 12); diff --git a/src/Particle/LongRange/tests/test_temp.cpp b/src/Particle/LongRange/tests/test_temp.cpp index 039197370a..1e7d7c8e6a 100644 --- a/src/Particle/LongRange/tests/test_temp.cpp +++ b/src/Particle/LongRange/tests/test_temp.cpp @@ -55,8 +55,7 @@ TEST_CASE("temp3d", "[lrhandler]") handler.initBreakup(ref); std::cout << "handler.MaxKshell is " << handler.MaxKshell << std::endl; - CHECK( (std::is_same::value ? - handler.MaxKshell == 78 : handler.MaxKshell >= 117 && handler.MaxKshell <= 128 )); + CHECK( handler.MaxKshell == 78); CHECK(Approx(handler.LR_rc) == 2.5); CHECK(Approx(handler.LR_kc) == 12); diff --git a/src/QMCHamiltonians/CoulombPBCAA.cpp b/src/QMCHamiltonians/CoulombPBCAA.cpp index 165cbb8402..1d8b0f7605 100644 --- a/src/QMCHamiltonians/CoulombPBCAA.cpp +++ b/src/QMCHamiltonians/CoulombPBCAA.cpp @@ -306,7 +306,7 @@ void CoulombPBCAA::mw_evaluatePerParticle(const RefVectorWithLeaderevaluate(pset.getSimulationCell().getKLists().kshell, PtclRhoK.rhok_r[s], PtclRhoK.rhok_i[s], + cpbcaa.AA->evaluate(pset.getSimulationCell().getKLists().getKShell(), PtclRhoK.rhok_r[s], PtclRhoK.rhok_i[s], PtclRhoK.eikr_r[i], PtclRhoK.eikr_i[i]); v_sample[i] += v1; Vlr += v1; @@ -420,7 +420,7 @@ CoulombPBCAA::Return_t CoulombPBCAA::evaluate_sp(ParticleSet& P) v1 = 0.0; for (int s = 0; s < NumSpecies; ++s) v1 += z * Zspec[s] * - AA->evaluate(P.getSimulationCell().getKLists().kshell, PtclRhoK.rhok_r[s], PtclRhoK.rhok_i[s], + AA->evaluate(P.getSimulationCell().getKLists().getKShell(), PtclRhoK.rhok_r[s], PtclRhoK.rhok_i[s], PtclRhoK.eikr_r[i], PtclRhoK.eikr_i[i]); V_samp(i) += v1; Vlr += v1; @@ -623,7 +623,7 @@ CoulombPBCAA::Return_t CoulombPBCAA::evalConsts(bool report) } // perform long-range Madelung sum const StructFact& PtclRhoK(Ps.getSK()); - v1 = AA->evaluate_slab(0, Ps.getSimulationCell().getKLists().kshell, PtclRhoK.eikr_r[0], PtclRhoK.eikr_i[0], + v1 = AA->evaluate_slab(0, Ps.getSimulationCell().getKLists().getKShell(), PtclRhoK.eikr_r[0], PtclRhoK.eikr_i[0], PtclRhoK.eikr_r[0], PtclRhoK.eikr_i[0]); if (report) app_log() << " LR Madelung = " << v1 << std::endl; @@ -794,7 +794,7 @@ CoulombPBCAA::Return_t CoulombPBCAA::evalLR(const ParticleSet& P) const { const RealType z = std::abs(dr[jat][slab_dir]); u += Zat[jat] * - AA->evaluate_slab(z, P.getSimulationCell().getKLists().kshell, PtclRhoK.eikr_r[iat], PtclRhoK.eikr_i[iat], + AA->evaluate_slab(z, P.getSimulationCell().getKLists().getKShell(), PtclRhoK.eikr_r[iat], PtclRhoK.eikr_i[iat], PtclRhoK.eikr_r[jat], PtclRhoK.eikr_i[jat]); } res += Zat[iat] * u; @@ -807,7 +807,7 @@ CoulombPBCAA::Return_t CoulombPBCAA::evalLR(const ParticleSet& P) const mRealType Z1 = Zspec[spec1]; for (int spec2 = spec1; spec2 < NumSpecies; spec2++) { - mRealType temp = AA->evaluate(P.getSimulationCell().getKLists().kshell, PtclRhoK.rhok_r[spec1], + mRealType temp = AA->evaluate(P.getSimulationCell().getKLists().getKShell(), PtclRhoK.rhok_r[spec1], PtclRhoK.rhok_i[spec1], PtclRhoK.rhok_r[spec2], PtclRhoK.rhok_i[spec2]); if (spec2 == spec1) temp *= 0.5; diff --git a/src/QMCHamiltonians/CoulombPBCAB.cpp b/src/QMCHamiltonians/CoulombPBCAB.cpp index dcddac7814..5c2cab788f 100644 --- a/src/QMCHamiltonians/CoulombPBCAB.cpp +++ b/src/QMCHamiltonians/CoulombPBCAB.cpp @@ -201,7 +201,7 @@ CoulombPBCAB::Return_t CoulombPBCAB::evaluate_sp(ParticleSet& P) v1 = 0.0; for (int s = 0; s < NumSpeciesA; s++) v1 += Zspec[s] * q * - AB->evaluate(pset_ions_.getSimulationCell().getKLists().kshell, RhoKA.rhok_r[s], RhoKA.rhok_i[s], + AB->evaluate(pset_ions_.getSimulationCell().getKLists().getKShell(), RhoKA.rhok_r[s], RhoKA.rhok_i[s], RhoKB.eikr_r[i], RhoKB.eikr_i[i]); Ve_samp(i) += v1; Vlr += v1; @@ -212,7 +212,7 @@ CoulombPBCAB::Return_t CoulombPBCAB::evaluate_sp(ParticleSet& P) v1 = 0.0; for (int s = 0; s < NumSpeciesB; s++) v1 += Qspec[s] * q * - AB->evaluate(P.getSimulationCell().getKLists().kshell, RhoKB.rhok_r[s], RhoKB.rhok_i[s], RhoKA.eikr_r[i], + AB->evaluate(P.getSimulationCell().getKLists().getKShell(), RhoKB.rhok_r[s], RhoKB.rhok_i[s], RhoKA.eikr_r[i], RhoKA.eikr_i[i]); Vi_samp(i) += v1; Vlr += v1; @@ -341,7 +341,7 @@ void CoulombPBCAB::mw_evaluatePerParticle(const RefVectorWithLeaderevaluate(pset_source.getSimulationCell().getKLists().kshell, RhoKA.rhok_r[s], + cpbcab.AB->evaluate(pset_source.getSimulationCell().getKLists().getKShell(), RhoKA.rhok_r[s], RhoKA.rhok_i[s], RhoKB.eikr_r[i], RhoKB.eikr_i[i]); ve_sample[i] += v1; Vlr += v1; @@ -353,7 +353,7 @@ void CoulombPBCAB::mw_evaluatePerParticle(const RefVectorWithLeaderevaluate(pset.getSimulationCell().getKLists().kshell, RhoKB.rhok_r[s], RhoKB.rhok_i[s], + cpbcab.AB->evaluate(pset.getSimulationCell().getKLists().getKShell(), RhoKB.rhok_r[s], RhoKB.rhok_i[s], RhoKA.eikr_r[i], RhoKA.eikr_i[i]); vi_sample[i] += v1; Vlr += v1; @@ -469,7 +469,7 @@ CoulombPBCAB::Return_t CoulombPBCAB::evalLR(ParticleSet& P) mRealType esum = 0.0; for (int j = 0; j < NumSpeciesB; j++) esum += Qspec[j] * - AB->evaluate(pset_ions_.getSimulationCell().getKLists().kshell, RhoKA.rhok_r[i], RhoKA.rhok_i[i], + AB->evaluate(pset_ions_.getSimulationCell().getKLists().getKShell(), RhoKA.rhok_r[i], RhoKA.rhok_i[i], RhoKB.rhok_r[j], RhoKB.rhok_i[j]); res += Zspec[i] * esum; } diff --git a/src/QMCHamiltonians/SkAllEstimator.cpp b/src/QMCHamiltonians/SkAllEstimator.cpp index 66cf16eb4f..da2dc7f3cf 100644 --- a/src/QMCHamiltonians/SkAllEstimator.cpp +++ b/src/QMCHamiltonians/SkAllEstimator.cpp @@ -31,9 +31,9 @@ SkAllEstimator::SkAllEstimator(ParticleSet& source, ParticleSet& target) NumIonSpecies = ions->getSpeciesSet().getTotalNum(); update_mode_.set(COLLECTABLE, 1); - NumK = source.getSimulationCell().getKLists().numk; + NumK = source.getSimulationCell().getKLists().getNumK(); OneOverN = 1.0 / static_cast(source.getTotalNum()); - Kshell = source.getSimulationCell().getKLists().kshell; + Kshell = source.getSimulationCell().getKLists().getKShell(); MaxKshell = Kshell.size() - 1; RhokTot_r.resize(NumK); @@ -46,7 +46,7 @@ SkAllEstimator::SkAllEstimator(ParticleSet& source, ParticleSet& target) OneOverDnk.resize(MaxKshell); for (int ks = 0; ks < MaxKshell; ks++) { - Kmag[ks] = std::sqrt(source.getSimulationCell().getKLists().ksq[Kshell[ks]]); + Kmag[ks] = std::sqrt(source.getSimulationCell().getKLists().getKSQWorking()[Kshell[ks]]); OneOverDnk[ks] = 1.0 / static_cast(Kshell[ks + 1] - Kshell[ks]); } hdf5_out = false; @@ -71,7 +71,7 @@ void SkAllEstimator::evaluateIonIon() for (int k = 0; k < NumK; k++) { - PosType kvec = ions->getSimulationCell().getKLists().kpts_cart[k]; + PosType kvec = ions->getSimulationCell().getKLists().getKptsCartWorking()[k]; filebuffer << kvec; for (int i = 0; i < NumIonSpecies; i++) @@ -230,7 +230,7 @@ void SkAllEstimator::registerCollectables(std::vector& h5desc, hdf_path hdf_name{name_}; h5desc.emplace_back(hdf_name / "kpoints"); auto& ohKPoints = h5desc.back(); - ohKPoints.addProperty(const_cast&>(ions->getSimulationCell().getKLists().kpts_cart), "value", + ohKPoints.addProperty(const_cast&>(ions->getSimulationCell().getKLists().getKptsCartWorking()), "value", file); // Add electron-electron S(k) diff --git a/src/QMCHamiltonians/SkEstimator.cpp b/src/QMCHamiltonians/SkEstimator.cpp index 18a5543696..a5f97f4e89 100644 --- a/src/QMCHamiltonians/SkEstimator.cpp +++ b/src/QMCHamiltonians/SkEstimator.cpp @@ -25,9 +25,9 @@ SkEstimator::SkEstimator(ParticleSet& source) sourcePtcl = &source; update_mode_.set(COLLECTABLE, 1); NumSpecies = source.getSpeciesSet().getTotalNum(); - NumK = source.getSimulationCell().getKLists().numk; + NumK = source.getSimulationCell().getKLists().getNumK(); OneOverN = 1.0 / static_cast(source.getTotalNum()); - Kshell = source.getSimulationCell().getKLists().kshell; + Kshell = source.getSimulationCell().getKLists().getKShell(); MaxKshell = Kshell.size() - 1; RhokTot_r.resize(NumK); RhokTot_i.resize(NumK); @@ -36,7 +36,7 @@ SkEstimator::SkEstimator(ParticleSet& source) OneOverDnk.resize(MaxKshell); for (int ks = 0; ks < MaxKshell; ks++) { - Kmag[ks] = std::sqrt(source.getSimulationCell().getKLists().ksq[Kshell[ks]]); + Kmag[ks] = std::sqrt(source.getSimulationCell().getKLists().getKSQWorking()[Kshell[ks]]); OneOverDnk[ks] = 1.0 / static_cast(Kshell[ks + 1] - Kshell[ks]); } hdf5_out = true; @@ -131,7 +131,7 @@ void SkEstimator::registerCollectables(std::vector& h5desc, hd hid_t k_set = H5Dcreate(file.getFileID(), kpath.c_str(), H5T_NATIVE_DOUBLE, k_space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); hid_t mem_space = H5Screate_simple(2, kdims, NULL); - auto* ptr = &(sourcePtcl->getSimulationCell().getKLists().kpts_cart[0][0]); + auto* ptr = &(sourcePtcl->getSimulationCell().getKLists().getKptsCartWorking()[0][0]); herr_t ret = H5Dwrite(k_set, H5T_NATIVE_DOUBLE, mem_space, k_space, H5P_DEFAULT, ptr); H5Dclose(k_set); H5Sclose(mem_space); diff --git a/src/QMCHamiltonians/SkPot.cpp b/src/QMCHamiltonians/SkPot.cpp index 5cc6055290..741962f141 100644 --- a/src/QMCHamiltonians/SkPot.cpp +++ b/src/QMCHamiltonians/SkPot.cpp @@ -22,9 +22,9 @@ SkPot::SkPot(ParticleSet& source) { sourcePtcl = &source; NumSpecies = source.getSpeciesSet().getTotalNum(); - NumK = source.getSimulationCell().getKLists().numk; + NumK = source.getSimulationCell().getKLists().getNumK(); OneOverN = 1.0 / static_cast(source.getTotalNum()); - Kshell = source.getSimulationCell().getKLists().kshell; + Kshell = source.getSimulationCell().getKLists().getKShell(); MaxKshell = Kshell.size() - 1; RhokTot.resize(NumK); Fk.resize(NumK); @@ -32,7 +32,7 @@ SkPot::SkPot(ParticleSet& source) OneOverDnk.resize(MaxKshell); for (int ks = 0; ks < MaxKshell; ks++) { - Kmag[ks] = std::sqrt(source.getSimulationCell().getKLists().ksq[Kshell[ks]]); + Kmag[ks] = std::sqrt(source.getSimulationCell().getKLists().getKSQWorking()[Kshell[ks]]); OneOverDnk[ks] = 1.0 / static_cast(Kshell[ks + 1] - Kshell[ks]); } } diff --git a/src/QMCHamiltonians/SkPot.h b/src/QMCHamiltonians/SkPot.h index 49bcf41ecd..dd67fc9ee7 100644 --- a/src/QMCHamiltonians/SkPot.h +++ b/src/QMCHamiltonians/SkPot.h @@ -42,8 +42,8 @@ class SkPot : public OperatorBase { for (int ki = 0; ki < NumK; ki++) { - RealType k = dot(sourcePtcl->getSimulationCell().getKLists().kpts_cart[ki], - sourcePtcl->getSimulationCell().getKLists().kpts_cart[ki]); + RealType k = dot(sourcePtcl->getSimulationCell().getKLists().getKptsCartWorking()[ki], + sourcePtcl->getSimulationCell().getKLists().getKptsCartWorking()[ki]); k = std::sqrt(k) - K_0; Fk[ki] = OneOverN * V_0 * std::exp(-k * k); // app_log()<& { h5desc.emplace_back(hdf_path(name_) / "kpoints"); auto& oh = h5desc.back(); - oh.addProperty(const_cast&>(Pinit.getSimulationCell().getKLists().kpts_cart), "value", file); + oh.addProperty(const_cast&>(Pinit.getSimulationCell().getKLists().getKptsCartWorking()), "value", file); std::vector ng(2); ng[0] = 2; diff --git a/src/QMCHamiltonians/StressPBC.cpp b/src/QMCHamiltonians/StressPBC.cpp index cd814eb071..286e6bddce 100644 --- a/src/QMCHamiltonians/StressPBC.cpp +++ b/src/QMCHamiltonians/StressPBC.cpp @@ -102,7 +102,7 @@ SymTensor StressPBC::evaluateLR_AB(ParticleSet& esum = 0.0; for (int j = 0; j < NumSpeciesB; j++) esum += Qspec[j] * - AA->evaluateStress(P.getSimulationCell().getKLists().kshell, RhoKA.rhok_r[i], RhoKA.rhok_i[i], + AA->evaluateStress(P.getSimulationCell().getKLists().getKShell(), RhoKA.rhok_r[i], RhoKA.rhok_i[i], RhoKB.rhok_r[j], RhoKB.rhok_i[j]); res += Zspec[i] * esum; } @@ -173,7 +173,7 @@ SymTensor StressPBC::evaluateLR_AA(ParticleSet& for (int spec2 = spec1; spec2 < NumSpecies; spec2++) { SymTensor temp = - AA->evaluateStress(P.getSimulationCell().getKLists().kshell, PtclRhoK.rhok_r[spec1], PtclRhoK.rhok_i[spec1], + AA->evaluateStress(P.getSimulationCell().getKLists().getKShell(), PtclRhoK.rhok_r[spec1], PtclRhoK.rhok_i[spec1], PtclRhoK.rhok_r[spec2], PtclRhoK.rhok_i[spec2]); if (spec2 == spec1) temp *= 0.5; diff --git a/src/QMCHamiltonians/tests/test_SkAllEstimator.cpp b/src/QMCHamiltonians/tests/test_SkAllEstimator.cpp index 7dc78c39d3..ae5d9389b2 100644 --- a/src/QMCHamiltonians/tests/test_SkAllEstimator.cpp +++ b/src/QMCHamiltonians/tests/test_SkAllEstimator.cpp @@ -208,7 +208,7 @@ TEST_CASE("SkAll", "[hamiltonian]") // In order to compare to analytic result, need the list // of k-vectors in cartesian coordinates. // Luckily, ParticleSet stores that in SK->getKLists().kpts_cart - int nkpts = elec->getSimulationCell().getKLists().numk; + int nkpts = elec->getSimulationCell().getKLists().getNumK(); std::cout << "\n"; std::cout << "SkAll results:\n"; std::cout << std::fixed; @@ -229,7 +229,7 @@ TEST_CASE("SkAll", "[hamiltonian]") std::cout << std::setprecision(5); for (int k = 0; k < nkpts; k++) { - auto kvec = elec->getSimulationCell().getKLists().kpts_cart[k]; + auto kvec = elec->getSimulationCell().getKLists().getKptsCartWorking()[k]; RealType kx = kvec[0]; RealType ky = kvec[1]; RealType kz = kvec[2]; diff --git a/src/QMCTools/QMCFiniteSize/QMCFiniteSize.cpp b/src/QMCTools/QMCFiniteSize/QMCFiniteSize.cpp index 65b1270ebe..e1afd75187 100644 --- a/src/QMCTools/QMCFiniteSize/QMCFiniteSize.cpp +++ b/src/QMCTools/QMCFiniteSize/QMCFiniteSize.cpp @@ -144,11 +144,11 @@ void QMCFiniteSize::initBreakup() app_log() << "=========================================================\n"; app_log() << " Initializing Long Range Breakup (Esler) \n"; app_log() << "=========================================================\n"; - P = ptclPool.getParticleSet("e"); - AA = LRCoulombSingleton::getHandler(*P); - myRcut = AA->get_rc(); + P = ptclPool.getParticleSet("e"); + AA = LRCoulombSingleton::getHandler(*P); + myRcut = AA->get_rc(); auto myGrid = LinearGrid(); - int ng = P->getLattice().num_ewald_grid_points; + int ng = P->getLattice().num_ewald_grid_points; myGrid.set(0, myRcut, ng); if (rVs == nullptr) { @@ -362,9 +362,9 @@ void QMCFiniteSize::initialize() rs = std::pow(3.0 / (4 * M_PI) * Vol / RealType(Ne), 1.0 / 3.0); rho = RealType(Ne) / Vol; Klist = P->getSimulationCell().getKLists(); - kpts = Klist.kpts; //These are in reduced coordinates. - //Easier to spline, but will have to convert - //for real space integration. + kpts = Klist.getKpts(); //These are in reduced coordinates. + //Easier to spline, but will have to convert + //for real space integration. if (!skparser->has_grid()) skparser->set_grid(kpts); @@ -375,14 +375,14 @@ void QMCFiniteSize::initialize() void QMCFiniteSize::printSkRawSphAvg(const std::vector& sk) { - std::vector vsk_1d(Klist.kshell.size()); + std::vector vsk_1d(Klist.getKShell().size()); // Average within each shell - for (int ks = 0; ks < Klist.kshell.size() - 1; ks++) + for (int ks = 0; ks < Klist.getKShell().size() - 1; ks++) { RealType u = 0; RealType n = 0; - for (int ki = Klist.kshell[ks]; ki < Klist.kshell[ks + 1]; ki++) + for (int ki = Klist.getKShell()[ks]; ki < Klist.getKShell()[ks + 1]; ki++) { u += sk[ki]; n++; @@ -401,13 +401,14 @@ void QMCFiniteSize::printSkRawSphAvg(const std::vector& sk) app_log() << "\nSpherically averaged raw S(k):\n"; app_log() << std::setw(12) << "k" << std::setw(12) << "S(k)" << std::setw(12) << "vk" << "\n"; - for (int ks = 0; ks < Klist.kshell.size() - 1; ks++) + for (int ks = 0; ks < Klist.getKShell().size() - 1; ks++) { - app_log() << std::setw(12) << std::setprecision(8) << std::sqrt(Klist.ksq[Klist.kshell[ks]]) << std::setw(12) - << std::setprecision(8) << vsk_1d[ks] << std::setw(12) << std::setprecision(8) << AA->Fk_symm[ks] << '\n'; + app_log() << std::setw(12) << std::setprecision(8) << std::sqrt(Klist.getKSQWorking()[Klist.getKShell()[ks]]) + << std::setw(12) << std::setprecision(8) << vsk_1d[ks] << std::setw(12) << std::setprecision(8) + << AA->Fk_symm[ks] << '\n'; } - if (vsk_1d[Klist.kshell.size() - 2] < 0.99) + if (vsk_1d[Klist.getKShell().size() - 2] < 0.99) { app_log() << "####################################################################\n"; app_log() << "WARNING: The S(k) in the largest kshell is less than 0.99\n"; @@ -460,7 +461,7 @@ void QMCFiniteSize::printSkSplineSphAvg(UBspline_3d_d* spline) QMCFiniteSize::RealType QMCFiniteSize::calcPotentialDiscrete(std::vector sk) { //This is the \frac{1}{Omega} \sum_{\mathbf{k}} \frac{v_k}{2} S(\mathbf{k}) term. - return 0.5 * AA->evaluate_w_sk(Klist.kshell, sk.data()); + return 0.5 * AA->evaluate_w_sk(Klist.getKShell(), sk.data()); } QMCFiniteSize::RealType QMCFiniteSize::calcPotentialInt(std::vector sk) @@ -468,7 +469,7 @@ QMCFiniteSize::RealType QMCFiniteSize::calcPotentialInt(std::vector sk auto spline = std::unique_ptr{getSkSpline(sk), destroy_Bspline}; RealType kmax = AA->get_kc(); - IndexType ngrid = 2 * Klist.kshell.size() - 1; //make a lager kmesh + IndexType ngrid = 2 * Klist.getKShell().size() - 1; //make a lager kmesh std::vector unigrid1d, k2vksk; RealType dk = kmax / ngrid; diff --git a/src/QMCWaveFunctions/ElectronGas/FreeOrbitalBuilder.cpp b/src/QMCWaveFunctions/ElectronGas/FreeOrbitalBuilder.cpp index 5861dc9d0f..06f4b2f12d 100644 --- a/src/QMCWaveFunctions/ElectronGas/FreeOrbitalBuilder.cpp +++ b/src/QMCWaveFunctions/ElectronGas/FreeOrbitalBuilder.cpp @@ -62,20 +62,20 @@ std::unique_ptr FreeOrbitalBuilder::createSPOSetFromXML(xmlNodePtr cur) #ifdef QMC_COMPLEX for (int ik = 1; ik < npw; ik++) { - kpts[ik] = klists.kpts_cart[ik - 1]; + kpts[ik] = klists.getKptsCartWorking()[ik - 1]; } #else - const int nktot = klists.kpts.size(); + const int nktot = klists.getKpts().size(); std::vector mkidx(npw, 0); int ik = 1; for (int jk = 0; jk < nktot; jk++) { // check if -k is already chosen - const int jmk = klists.minusk[jk]; + const int jmk = klists.getMinusK(jk); if (in_list(jk, mkidx)) continue; // if not, then add this kpoint - kpts[ik] = klists.kpts_cart[jk]; + kpts[ik] = klists.getKptsCartWorking()[jk]; mkidx[ik] = jmk; // keep track of its minus ik++; if (ik >= npw) diff --git a/src/QMCWaveFunctions/Fermion/BackflowBuilder.cpp b/src/QMCWaveFunctions/Fermion/BackflowBuilder.cpp index a8b8a39e99..c01d641049 100644 --- a/src/QMCWaveFunctions/Fermion/BackflowBuilder.cpp +++ b/src/QMCWaveFunctions/Fermion/BackflowBuilder.cpp @@ -403,8 +403,8 @@ std::unique_ptr BackflowBuilder::addRPA(xmlNodePtr cur) Rs = 100.0; } } - int indx = targetPtcl.getSimulationCell().getKLists().ksq.size() - 1; - RealType Kc_max = std::pow(targetPtcl.getSimulationCell().getKLists().ksq[indx], 0.5); + int indx = targetPtcl.getSimulationCell().getKLists().getKSQWorking().size() - 1; + RealType Kc_max = std::pow(targetPtcl.getSimulationCell().getKLists().getKSQWorking()[indx], 0.5); if (Kc < 0) { Kc = 2.0 * std::pow(2.25 * M_PI, 1.0 / 3.0) / tlen; @@ -563,7 +563,7 @@ void BackflowBuilder::makeLongRange_twoBody(xmlNodePtr cur, Backflow_ee_kSpace* { fout << std::sqrt(targetPtcl.getSimulationCell() .getKLists() - .ksq[targetPtcl.getSimulationCell().getKLists().kshell[i]]) + .getKSQWorking()[targetPtcl.getSimulationCell().getKLists().getKShell()[i]]) << " " << yk[i] << std::endl; } fout.close(); diff --git a/src/QMCWaveFunctions/Fermion/Backflow_ee_kSpace.h b/src/QMCWaveFunctions/Fermion/Backflow_ee_kSpace.h index 5178de3f5c..35f45dd4de 100644 --- a/src/QMCWaveFunctions/Fermion/Backflow_ee_kSpace.h +++ b/src/QMCWaveFunctions/Fermion/Backflow_ee_kSpace.h @@ -64,7 +64,7 @@ class Backflow_ee_kSpace : public BackflowFunctionBase { NumKShells = yk.size(); Fk = yk; - NumKVecs = P.getSimulationCell().getKLists().kshell[NumKShells + 1]; + NumKVecs = P.getSimulationCell().getKLists().getKShell()[NumKShells + 1]; Rhok.resize(NumKVecs); if (Optimize) numParams = NumKShells; diff --git a/src/QMCWaveFunctions/Jastrow/J2KECorrection.h b/src/QMCWaveFunctions/Jastrow/J2KECorrection.h index 82d89d519d..a2ff24e209 100644 --- a/src/QMCWaveFunctions/Jastrow/J2KECorrection.h +++ b/src/QMCWaveFunctions/Jastrow/J2KECorrection.h @@ -46,7 +46,7 @@ class J2KECorrection num_elec_in_groups_.push_back(targetPtcl.last(i) - targetPtcl.first(i)); if (SK_enabled) - G0mag = std::sqrt(targetPtcl.getSimulationCell().getKLists().ksq[0]); + G0mag = std::sqrt(targetPtcl.getSimulationCell().getKLists().getKSQWorking()[0]); } RT computeKEcorr() diff --git a/src/QMCWaveFunctions/Jastrow/RPAJastrow.cpp b/src/QMCWaveFunctions/Jastrow/RPAJastrow.cpp index f18c00c28b..3df82f28fa 100644 --- a/src/QMCWaveFunctions/Jastrow/RPAJastrow.cpp +++ b/src/QMCWaveFunctions/Jastrow/RPAJastrow.cpp @@ -96,8 +96,8 @@ void RPAJastrow::buildOrbital(const std::string& name, Rs = 100.0; } } - int indx = targetPtcl.getSimulationCell().getKLists().ksq.size() - 1; - double Kc_max = std::pow(targetPtcl.getSimulationCell().getKLists().ksq[indx], 0.5); + int indx = targetPtcl.getSimulationCell().getKLists().getKSQWorking().size() - 1; + double Kc_max = std::pow(targetPtcl.getSimulationCell().getKLists().getKSQWorking()[indx], 0.5); if (Kc < 0) { Kc = 2.0 * std::pow(2.25 * M_PI, 1.0 / 3.0) / tlen; diff --git a/src/QMCWaveFunctions/Jastrow/RadialJastrowBuilder.cpp b/src/QMCWaveFunctions/Jastrow/RadialJastrowBuilder.cpp index 7748688d83..e006feb6a2 100644 --- a/src/QMCWaveFunctions/Jastrow/RadialJastrowBuilder.cpp +++ b/src/QMCWaveFunctions/Jastrow/RadialJastrowBuilder.cpp @@ -267,9 +267,9 @@ void RadialJastrowBuilder::computeJ2uk(const std::vector& functors throw std::runtime_error("Error generating filename"); fout = fopen(fname.data(), "w"); } - for (int iG = 0; iG < targetPtcl.getSimulationCell().getKLists().ksq.size(); iG++) + for (int iG = 0; iG < targetPtcl.getSimulationCell().getKLists().getKSQWorking().size(); iG++) { - RealType Gmag = std::sqrt(targetPtcl.getSimulationCell().getKLists().ksq[iG]); + RealType Gmag = std::sqrt(targetPtcl.getSimulationCell().getKLists().getKSQWorking()[iG]); RealType sum = 0.0; RealType uk = 0.0; for (int i = 0; i < targetPtcl.groups(); i++) diff --git a/src/Utilities/for_testing/checkVector.hpp b/src/Utilities/for_testing/checkVector.hpp index 4e7df62ab8..91c9ef4f88 100644 --- a/src/Utilities/for_testing/checkVector.hpp +++ b/src/Utilities/for_testing/checkVector.hpp @@ -64,6 +64,33 @@ struct CheckVectorResult * \param[in] eps - add a tolerance for Catch Approx checks. Default to same as in Approx. * The semantics of the return value are discussed above. */ +template +CheckVectorResult checkVector(const std::vector>& a_vec, + const std::vector>& b_vec, + const bool check_all = false) +{ + // This allows use to check a padded b matrix with a nonpadded a + if (a_vec.size() > b_vec.size()) + return {false, "b_vec is too small for a_vec to be a checkable segment"}; + std::stringstream result_msg; + auto vectorElementError = [&result_msg](int i, auto& a_vec, auto& b_vec) { + result_msg << "checkVector found bad element at " << i << " " << a_vec[i] << " != " << b_vec[i] << '\n'; + }; + bool all_elements_match = true; + for (int i = 0; i < a_vec.size(); i++) + { + bool equality = a_vec[i] == b_vec[i]; + if (!equality) + { + vectorElementError(i, a_vec, b_vec); + all_elements_match = false; + if (!check_all) + return {false, result_msg.str()}; + } + } + return {all_elements_match, result_msg.str()}; +} + template CheckVectorResult checkVector(const M1& a_vec, const M2& b_vec, diff --git a/src/Utilities/tests/for_testing/test_checkVector.cpp b/src/Utilities/tests/for_testing/test_checkVector.cpp index e37f6eb559..704a657ff3 100644 --- a/src/Utilities/tests/for_testing/test_checkVector.cpp +++ b/src/Utilities/tests/for_testing/test_checkVector.cpp @@ -10,6 +10,7 @@ ////////////////////////////////////////////////////////////////////////////////////// #include "catch.hpp" +#include "OhmmsPETE/TinyVector.h" #include "checkVector.hpp" #include #include From ccb67560b2f13fbbaf5ecdb328d359e12f457469 Mon Sep 17 00:00:00 2001 From: "Peter Doak (epd)" Date: Thu, 20 Mar 2025 16:53:40 -0400 Subject: [PATCH 2/8] update heg mixed reference values now closer to full prec --- tests/heg/heg_14_gamma/CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/heg/heg_14_gamma/CMakeLists.txt b/tests/heg/heg_14_gamma/CMakeLists.txt index 3672d5f95e..ab0e4e83fe 100644 --- a/tests/heg/heg_14_gamma/CMakeLists.txt +++ b/tests/heg/heg_14_gamma/CMakeLists.txt @@ -116,9 +116,9 @@ qmc_run_and_check( # batch deterministic if(QMC_MIXED_PRECISION) - list(APPEND DET_HEG14GSSJ_SCALARS "totenergy" "-1.24693322 0.00002") + list(APPEND DET_HEG14GSSJ_SCALARS "totenergy" "-1.24699382 0.00002") list(APPEND DET_HEG14GSSJ_SCALARS "kinetic" "0.77668442 0.000005") - list(APPEND DET_HEG14GSSJ_SCALARS "potential" "-2.02361764 0.00002") + list(APPEND DET_HEG14GSSJ_SCALARS "potential" "-2.02367775 0.00002") else() list(APPEND DET_HEG14GSSJ_SCALARS "totenergy" "-1.24699388 0.000001") list(APPEND DET_HEG14GSSJ_SCALARS "kinetic" "0.77668380 0.000001") From ed773a392aaed75f3d892953ee9bb22348e82e93 Mon Sep 17 00:00:00 2001 From: "Peter Doak (epd)" Date: Fri, 21 Mar 2025 11:38:13 -0400 Subject: [PATCH 3/8] documentation and formatting consistent full/mixed prec KContainer --- src/Particle/LongRange/KContainer.cpp | 4 +++ src/Particle/LongRange/KContainer.h | 36 +++++++++++++++------------ 2 files changed, 24 insertions(+), 16 deletions(-) diff --git a/src/Particle/LongRange/KContainer.cpp b/src/Particle/LongRange/KContainer.cpp index f7f8d4fe57..1950d62fa3 100644 --- a/src/Particle/LongRange/KContainer.cpp +++ b/src/Particle/LongRange/KContainer.cpp @@ -26,6 +26,7 @@ namespace qmcplusplus template const std::vector::AppPosition>& KContainerT::getKptsCartWorking() const { + // This is an `if constexpr` so it should not cost a branch at runtime. if constexpr (std::is_same_v) return kpts_cart_; else @@ -35,6 +36,7 @@ const std::vector::AppPosition>& KContainerT::g template const std::vector& KContainerT::getKSQWorking() const { + // This is an `if constexpr` so it should not cost a branch at runtime. if constexpr (std::is_same::value) return ksq_; else @@ -269,6 +271,8 @@ void KContainerT::BuildKLists(const CrystalLattice::value) { + // This copy implicity does the precision reduction. + // the working vectors are not used or initialized for full precision builds. std::copy(kpts_cart_.begin(), kpts_cart_.end(), std::back_inserter(kpts_cart_working_)); std::copy(ksq_.begin(), ksq_.end(), std::back_inserter(ksq_working_)); } diff --git a/src/Particle/LongRange/KContainer.h b/src/Particle/LongRange/KContainer.h index 650f7baa69..5325110de6 100644 --- a/src/Particle/LongRange/KContainer.h +++ b/src/Particle/LongRange/KContainer.h @@ -47,13 +47,22 @@ class KContainerT using ParticleLayout = CrystalLattice; const auto& get_kpts_cart_soa() const { return kpts_cart_soa_; } - const std::vector>& getKpts() const { return kpts_; } + + /** @ingroup Working Precision accessors + * @brief return reference kpt values in the working precision of build + * these encapsulate the cached reduced precision or pass through of the + * full precision values. + * @{ + */ + ///get cartesian kpt representation in the working precision of the application. const std::vector& getKptsCartWorking() const; + ///get ksqr in the working precision of the application. + const std::vector& getKSQWorking() const; + /// @} const std::vector& getKShell() const { return kshell; }; - const std::vector& getKSQWorking() const; int getMinusK(int k) const; @@ -82,17 +91,20 @@ class KContainerT */ TinyVector mmax; - /** K-vector in reduced coordinates + /** k-vector in reduced coordinates */ std::vector> kpts_; - /** K-vector in Cartesian coordinates - */ + /// k-vectors in Cartesian coordinates std::vector kpts_cart_; - std::vector kpts_cart_working_; - /** squre of kpts in Cartesian coordniates - */ + /// k squared at full precision std::vector ksq_; + /** @ingroup Cached Working Precision values + * @brief only used or initialized for mixed precision + * @{ + */ + std::vector kpts_cart_working_; std::vector ksq_working_; + /// @} /** Given a k index, return index to -k */ @@ -100,14 +112,6 @@ class KContainerT /** kpts which belong to the ith-shell [kshell[i], kshell[i+1]) */ std::vector kshell; - /** k points sorted by the |k| excluding |k|=0 - * - * The first for |k| - * The second for a map to the full index. The size of the second is the degeneracy. - */ - //std::map*> kpts_sorted; - -private: /** compute approximate parallelpiped that surrounds kc * @param lattice supercell */ From 17f0313e50ad7adfa0429c167b836347c56745f7 Mon Sep 17 00:00:00 2001 From: "Peter Doak (epd)" Date: Fri, 21 Mar 2025 11:50:34 -0400 Subject: [PATCH 4/8] Revert change from Lattice, always test both precisions KContainer --- .../LongRange/tests/test_kcontainer.cpp | 53 +++++++++++-------- 1 file changed, 30 insertions(+), 23 deletions(-) diff --git a/src/Particle/LongRange/tests/test_kcontainer.cpp b/src/Particle/LongRange/tests/test_kcontainer.cpp index 8057fe022d..c011192477 100644 --- a/src/Particle/LongRange/tests/test_kcontainer.cpp +++ b/src/Particle/LongRange/tests/test_kcontainer.cpp @@ -32,7 +32,7 @@ TEST_CASE("kcontainer at gamma in 3D", "[longrange]") const std::vector kcs = {blat, std::sqrt(2) * blat, std::sqrt(3) * blat}; const std::vector nks = {6, 18, 26}; - CrystalLattice lattice; + Lattice lattice; lattice.R.diagonal(1.0); lattice.set(lattice.R); // compute Rv and Gv from R @@ -40,7 +40,8 @@ TEST_CASE("kcontainer at gamma in 3D", "[longrange]") for (int ik = 0; ik < kcs.size(); ik++) { const double kc = kcs[ik] + 1e-6; - klists.updateKLists>::Scalar_t>(lattice, kc, ndim); + klists.updateKLists>::Scalar_t>(lattice, kc, + ndim); CHECK(klists.getKpts().size() == nks[ik]); } const int mxk = klists.getKpts().size(); @@ -68,7 +69,7 @@ TEST_CASE("kcontainer at twist in 3D", "[longrange]") // twist one shell of kvectors const double kc = blat + 1e-6; - CrystalLattice lattice; + Lattice lattice; lattice.R.diagonal(1.0); lattice.set(lattice.R); // compute Rv and Gv from R @@ -76,12 +77,14 @@ TEST_CASE("kcontainer at twist in 3D", "[longrange]") PosType twist; twist[0] = 0.1; - klists.updateKLists>::Scalar_t>(lattice, kc, ndim, twist); + klists.updateKLists>::Scalar_t>(lattice, kc, + ndim, twist); CHECK(klists.getKpts().size() == 1); CHECK(klists.getKptsCartWorking()[0][0] == Approx(blat * (twist[0] - 1))); twist = {-0.5, 0, 0.5}; - klists.updateKLists>::Scalar_t>(lattice, kc, ndim, twist); + klists.updateKLists>::Scalar_t>(lattice, kc, + ndim, twist); int gvecs[3][3] = {{0, 0, -1}, {1, 0, -1}, {1, 0, 0}}; CHECK(klists.getKpts().size() == 3); for (int ik = 0; ik < klists.getKpts().size(); ik++) @@ -99,7 +102,7 @@ TEST_CASE("kcontainer at gamma in 2D", "[longrange]") const std::vector kcs = {blat, std::sqrt(2) * blat, 2 * blat}; const std::vector nks = {4, 8, 12}; - CrystalLattice lattice; + Lattice lattice; lattice.R.diagonal(1.0); lattice.set(lattice.R); // compute Rv and Gv from R @@ -107,7 +110,8 @@ TEST_CASE("kcontainer at gamma in 2D", "[longrange]") for (int ik = 0; ik < kcs.size(); ik++) { const double kc = kcs[ik] + 1e-6; - klists.updateKLists>::Scalar_t>(lattice, kc, ndim); + klists.updateKLists>::Scalar_t>(lattice, kc, + ndim); CHECK(klists.getKpts().size() == nks[ik]); } const int mxk = klists.getKpts().size(); @@ -135,7 +139,7 @@ TEST_CASE("kcontainer at twist in 2D", "[longrange]") // twist one shell of kvectors const double kc = blat + 1e-6; - CrystalLattice lattice; + Lattice lattice; lattice.R.diagonal(1.0); lattice.set(lattice.R); // compute Rv and Gv from R @@ -143,12 +147,14 @@ TEST_CASE("kcontainer at twist in 2D", "[longrange]") PosType twist; twist[0] = 0.1; - klists.updateKLists>::Scalar_t>(lattice, kc, ndim, twist); + klists.updateKLists>::Scalar_t>(lattice, kc, + ndim, twist); CHECK(klists.getKpts().size() == 1); CHECK(klists.getKptsCartWorking()[0][0] == Approx(blat * (twist[0] - 1))); twist[1] = 0.1; - klists.updateKLists>::Scalar_t>(lattice, kc, ndim, twist); + klists.updateKLists>::Scalar_t>(lattice, kc, + ndim, twist); CHECK(klists.getKpts().size() == 2); CHECK(klists.getKptsCartWorking()[0][0] == Approx(blat * (twist[0] - 1))); CHECK(klists.getKptsCartWorking()[0][1] == Approx(blat * twist[1])); @@ -156,7 +162,8 @@ TEST_CASE("kcontainer at twist in 2D", "[longrange]") CHECK(klists.getKptsCartWorking()[1][1] == Approx(blat * (twist[1] - 1))); twist = {-0.5, 0.5, 0}; - klists.updateKLists>::Scalar_t>(lattice, kc, ndim, twist); + klists.updateKLists>::Scalar_t>(lattice, kc, + ndim, twist); CHECK(klists.getKpts().size() == 3); //for (int ik=0;ik<3;ik++) // app_log() << klists.getKptsCartWorking()[ik] << std::endl; @@ -170,19 +177,12 @@ TEST_CASE("kcontainer at twist in 2D", "[longrange]") TEST_CASE("kcontainer for diamond", "[longrange]") { - // Communicate* comm = OHMMS::Controller; - // ParticleSetPool particle_pool{MinimalParticlePool::make_diamondC_1x1x1(comm)}; - // ParticleSet pset_ions{*(particle_pool.getParticleSet("ion"))}; - int ndim = 3; - // auto kcs = pset_ions.getSimulationCell().getKLists().numk; - // std::cout << "kcs:" << kcs << '\n'; - - using Real = QMCTraits::RealType; + using Real = QMCTraits::RealType; using FullPrecReal = QMCTraits::FullPrecRealType; - - CrystalLattice lattice; + + Lattice lattice; lattice.BoxBConds = {true, true, true}; Tensor lattice_mat{3.37316115, 3.37316115, 0.00000000, 0.00000000, 3.37316115, 3.37316115, 3.37316115, 0.00000000, 3.37316115}; @@ -190,7 +190,8 @@ TEST_CASE("kcontainer for diamond", "[longrange]") lattice.LR_dim_cutoff = 15; SimulationCell cell(lattice); - const KContainerT& klists = cell.getKLists(); + // In this test we check that both full and reduced precision KContainers agree + const KContainerT& klists = cell.getKLists(); std::remove_cv_t> kpoint_lists = { {-1, -1, -1}, {-1, 0, 0}, {0, -1, 0}, {0, 0, -1}, {0, 0, 1}, {0, 1, 0}, {1, 0, 0}, {1, 1, 1}, {-1, -1, 0}, {-1, 0, -1}, {0, -1, -1}, {0, 1, 1}, {1, 0, 1}, {1, 1, 0}, {-2, -1, -1}, {-1, -2, -1}, @@ -269,12 +270,18 @@ TEST_CASE("kcontainer for diamond", "[longrange]") {1, 5, 0}, {1, 5, 4}, {3, -1, 4}, {3, 4, -1}, {4, -1, -1}, {4, -1, 3}, {4, 1, 5}, {4, 3, -1}, {4, 5, 1}, {4, 5, 5}, {5, 0, 1}, {5, 1, 0}, {5, 1, 4}, {5, 4, 1}, {5, 4, 5}, {5, 5, 4}, }; - double tolerance = 0.1; { INFO("Checking kpoint_lists"); auto check = checkVector(klists.getKpts(), kpoint_lists, true); CHECKED_ELSE(check.result) { FAIL(check.result_message); } } + + const KContainerT& klists_full = cell.getKLists(); + { + INFO("Checking kpoint_lists"); + auto check = checkVector(klists_full.getKpts(), kpoint_lists, true); + CHECKED_ELSE(check.result) { FAIL(check.result_message); } + } } From 1505dad2fa0f454ba4e7f0e0df6a83b976bcf0b8 Mon Sep 17 00:00:00 2001 From: "Peter Doak (epd)" Date: Fri, 21 Mar 2025 15:46:58 -0400 Subject: [PATCH 5/8] follow through the simplication of Lattice always beging full prec --- src/Particle/LongRange/KContainer.cpp | 56 ++---------------- src/Particle/LongRange/KContainer.h | 55 +----------------- src/Particle/LongRange/LRBreakup.h | 2 +- .../LongRange/tests/test_kcontainer.cpp | 24 +++----- src/Particle/SimulationCell.cpp | 58 ++++++++++--------- src/Particle/SimulationCell.h | 21 +++---- 6 files changed, 60 insertions(+), 156 deletions(-) diff --git a/src/Particle/LongRange/KContainer.cpp b/src/Particle/LongRange/KContainer.cpp index 1950d62fa3..bebb43fdb7 100644 --- a/src/Particle/LongRange/KContainer.cpp +++ b/src/Particle/LongRange/KContainer.cpp @@ -24,7 +24,8 @@ namespace qmcplusplus { template -const std::vector::AppPosition>& KContainerT::getKptsCartWorking() const +const std::vector::AppPosition>& KContainerT< + REAL>::getKptsCartWorking() const { // This is an `if constexpr` so it should not cost a branch at runtime. if constexpr (std::is_same_v) @@ -51,8 +52,7 @@ int KContainerT::getMinusK(int k) const } template -template -void KContainerT::updateKLists(const CrystalLattice& lattice, +void KContainerT::updateKLists(const Lattice& lattice, FullPrecReal kc, unsigned ndim, const Position& twist, @@ -73,8 +73,7 @@ void KContainerT::updateKLists(const CrystalLattice -template -void KContainerT::findApproxMMax(const CrystalLattice& lattice, unsigned ndim) +void KContainerT::findApproxMMax(const Lattice& lattice, unsigned ndim) { //Estimate the size of the parallelpiped that encompasses a sphere of kcutoff. //mmax is stored as integer translations of the reciprocal cell vectors. @@ -136,8 +135,7 @@ void KContainerT::findApproxMMax(const CrystalLattice -template -void KContainerT::BuildKLists(const CrystalLattice& lattice, +void KContainerT::BuildKLists(const Lattice& lattice, const Position& twist, bool useSphere) { @@ -317,52 +315,8 @@ void KContainerT::BuildKLists(const CrystalLattice; -template void KContainerT::updateKLists(const CrystalLattice& lattice, - FullPrecReal kc, - unsigned ndim, - const Position& twist = Position(), - bool useSphere = true); -template void KContainerT::updateKLists(const CrystalLattice& lattice, - FullPrecReal kc, - unsigned ndim, - const Position& twist = Position(), - bool useSphere = true); -template void KContainerT::findApproxMMax(const CrystalLattice& lattice, unsigned ndim); -template void KContainerT::findApproxMMax(const CrystalLattice& lattice, - unsigned ndim); - - template class KContainerT; -template void KContainerT::updateKLists(const CrystalLattice& lattice, - FullPrecReal kc, - unsigned ndim, - const Position& twist = Position(), - bool useSphere = true); -template void KContainerT::updateKLists(const CrystalLattice& lattice, - FullPrecReal kc, - unsigned ndim, - const Position& twist = Position(), - bool useSphere = true); -template void KContainerT::findApproxMMax(const CrystalLattice& lattice, - unsigned ndim); -template void KContainerT::findApproxMMax(const CrystalLattice& lattice, - unsigned ndim); #else template class KContainerT; -template void KContainerT::updateKLists(const CrystalLattice& lattice, - FullPrecReal kc, - unsigned ndim, - const Position& twist = Position(), - bool useSphere = true); -template void KContainerT::updateKLists(const CrystalLattice& lattice, - FullPrecReal kc, - unsigned ndim, - const Position& twist = Position(), - bool useSphere = true); -template void KContainerT::findApproxMMax(const CrystalLattice& lattice, - unsigned ndim); -template void KContainerT::findApproxMMax(const CrystalLattice& lattice, - unsigned ndim); - #endif } // namespace qmcplusplus diff --git a/src/Particle/LongRange/KContainer.h b/src/Particle/LongRange/KContainer.h index 5325110de6..a78631ab75 100644 --- a/src/Particle/LongRange/KContainer.h +++ b/src/Particle/LongRange/KContainer.h @@ -43,9 +43,6 @@ class KContainerT FullPrecReal kcutoff; public: - //Typedef for the lattice-type - using ParticleLayout = CrystalLattice; - const auto& get_kpts_cart_soa() const { return kpts_cart_soa_; } const std::vector>& getKpts() const { return kpts_; } @@ -74,8 +71,7 @@ class KContainerT * @param twist shifts the center of the grid of k-vectors * @param useSphere if true, use the |K| */ - template - void updateKLists(const CrystalLattice& lattice, + void updateKLists(const Lattice& lattice, FullPrecReal kc, unsigned ndim, const Position& twist = Position(), @@ -115,11 +111,9 @@ class KContainerT /** compute approximate parallelpiped that surrounds kc * @param lattice supercell */ - template - void findApproxMMax(const CrystalLattice& lattice, unsigned ndim); + void findApproxMMax(const Lattice& lattice, unsigned ndim); /** construct the container for k-vectors */ - template - void BuildKLists(const CrystalLattice& lattice, const Position& twist, bool useSphere); + void BuildKLists(const Lattice& lattice, const Position& twist, bool useSphere); /** K-vector in Cartesian coordinates in SoA layout */ @@ -130,52 +124,9 @@ using KContainer = KContainerT; #ifdef MIXED_PRECISION extern template class KContainerT; -extern template void KContainerT::updateKLists(const CrystalLattice& lattice, - FullPrecReal kc, - unsigned ndim, - const Position& twist = Position(), - bool useSphere = true); -extern template void KContainerT::updateKLists(const CrystalLattice& lattice, - FullPrecReal kc, - unsigned ndim, - const Position& twist = Position(), - bool useSphere = true); -extern template void KContainerT::findApproxMMax(const CrystalLattice& lattice, - unsigned ndim); -extern template void KContainerT::findApproxMMax(const CrystalLattice& lattice, - unsigned ndim); - extern template class KContainerT; -extern template void KContainerT::updateKLists(const CrystalLattice& lattice, - FullPrecReal kc, - unsigned ndim, - const Position& twist = Position(), - bool useSphere = true); -extern template void KContainerT::updateKLists(const CrystalLattice& lattice, - FullPrecReal kc, - unsigned ndim, - const Position& twist = Position(), - bool useSphere = true); -extern template void KContainerT::findApproxMMax(const CrystalLattice& lattice, - unsigned ndim); -extern template void KContainerT::findApproxMMax(const CrystalLattice& lattice, - unsigned ndim); #else extern template class KContainerT; -extern template void KContainerT::updateKLists(const CrystalLattice& lattice, - FullPrecReal kc, - unsigned ndim, - const Position& twist = Position(), - bool useSphere = true); -extern template void KContainerT::updateKLists(const CrystalLattice& lattice, - FullPrecReal kc, - unsigned ndim, - const Position& twist = Position(), - bool useSphere = true); -extern template void KContainerT::findApproxMMax(const CrystalLattice& lattice, - unsigned ndim); -extern template void KContainerT::findApproxMMax(const CrystalLattice& lattice, - unsigned ndim); #endif } // namespace qmcplusplus diff --git a/src/Particle/LongRange/LRBreakup.h b/src/Particle/LongRange/LRBreakup.h index e9b6bb906c..51505db7d0 100644 --- a/src/Particle/LongRange/LRBreakup.h +++ b/src/Particle/LongRange/LRBreakup.h @@ -335,7 +335,7 @@ int LRBreakup::SetupKVecs(mRealType kc, mRealType kcont, mRealType //Add low |k| ( < kcont) k-points with exact degeneracy KContainer kexact; const auto& basis_lattice = Basis.get_Lattice(); - kexact.updateKLists>::Real>(basis_lattice, kcont, basis_lattice.ndim); + kexact.updateKLists(basis_lattice, kcont, basis_lattice.ndim); bool findK = true; mRealType kc2 = kc * kc; //use at least one shell diff --git a/src/Particle/LongRange/tests/test_kcontainer.cpp b/src/Particle/LongRange/tests/test_kcontainer.cpp index c011192477..2b31eedafc 100644 --- a/src/Particle/LongRange/tests/test_kcontainer.cpp +++ b/src/Particle/LongRange/tests/test_kcontainer.cpp @@ -40,7 +40,7 @@ TEST_CASE("kcontainer at gamma in 3D", "[longrange]") for (int ik = 0; ik < kcs.size(); ik++) { const double kc = kcs[ik] + 1e-6; - klists.updateKLists>::Scalar_t>(lattice, kc, + klists.updateKLists(lattice, kc, ndim); CHECK(klists.getKpts().size() == nks[ik]); } @@ -77,13 +77,13 @@ TEST_CASE("kcontainer at twist in 3D", "[longrange]") PosType twist; twist[0] = 0.1; - klists.updateKLists>::Scalar_t>(lattice, kc, + klists.updateKLists(lattice, kc, ndim, twist); CHECK(klists.getKpts().size() == 1); CHECK(klists.getKptsCartWorking()[0][0] == Approx(blat * (twist[0] - 1))); twist = {-0.5, 0, 0.5}; - klists.updateKLists>::Scalar_t>(lattice, kc, + klists.updateKLists(lattice, kc, ndim, twist); int gvecs[3][3] = {{0, 0, -1}, {1, 0, -1}, {1, 0, 0}}; CHECK(klists.getKpts().size() == 3); @@ -110,7 +110,7 @@ TEST_CASE("kcontainer at gamma in 2D", "[longrange]") for (int ik = 0; ik < kcs.size(); ik++) { const double kc = kcs[ik] + 1e-6; - klists.updateKLists>::Scalar_t>(lattice, kc, + klists.updateKLists(lattice, kc, ndim); CHECK(klists.getKpts().size() == nks[ik]); } @@ -147,13 +147,13 @@ TEST_CASE("kcontainer at twist in 2D", "[longrange]") PosType twist; twist[0] = 0.1; - klists.updateKLists>::Scalar_t>(lattice, kc, + klists.updateKLists(lattice, kc, ndim, twist); CHECK(klists.getKpts().size() == 1); CHECK(klists.getKptsCartWorking()[0][0] == Approx(blat * (twist[0] - 1))); twist[1] = 0.1; - klists.updateKLists>::Scalar_t>(lattice, kc, + klists.updateKLists(lattice, kc, ndim, twist); CHECK(klists.getKpts().size() == 2); CHECK(klists.getKptsCartWorking()[0][0] == Approx(blat * (twist[0] - 1))); @@ -162,7 +162,7 @@ TEST_CASE("kcontainer at twist in 2D", "[longrange]") CHECK(klists.getKptsCartWorking()[1][1] == Approx(blat * (twist[1] - 1))); twist = {-0.5, 0.5, 0}; - klists.updateKLists>::Scalar_t>(lattice, kc, + klists.updateKLists(lattice, kc, ndim, twist); CHECK(klists.getKpts().size() == 3); //for (int ik=0;ik<3;ik++) @@ -190,8 +190,7 @@ TEST_CASE("kcontainer for diamond", "[longrange]") lattice.LR_dim_cutoff = 15; SimulationCell cell(lattice); - // In this test we check that both full and reduced precision KContainers agree - const KContainerT& klists = cell.getKLists(); + const KContainer& klists = cell.getKLists(); std::remove_cv_t> kpoint_lists = { {-1, -1, -1}, {-1, 0, 0}, {0, -1, 0}, {0, 0, -1}, {0, 0, 1}, {0, 1, 0}, {1, 0, 0}, {1, 1, 1}, {-1, -1, 0}, {-1, 0, -1}, {0, -1, -1}, {0, 1, 1}, {1, 0, 1}, {1, 1, 0}, {-2, -1, -1}, {-1, -2, -1}, @@ -275,13 +274,6 @@ TEST_CASE("kcontainer for diamond", "[longrange]") auto check = checkVector(klists.getKpts(), kpoint_lists, true); CHECKED_ELSE(check.result) { FAIL(check.result_message); } } - - const KContainerT& klists_full = cell.getKLists(); - { - INFO("Checking kpoint_lists"); - auto check = checkVector(klists_full.getKpts(), kpoint_lists, true); - CHECKED_ELSE(check.result) { FAIL(check.result_message); } - } } diff --git a/src/Particle/SimulationCell.cpp b/src/Particle/SimulationCell.cpp index cbd9b5f053..365ad5ea7c 100644 --- a/src/Particle/SimulationCell.cpp +++ b/src/Particle/SimulationCell.cpp @@ -2,64 +2,70 @@ // 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. +// Copyright (c) 2025 QMCPACK developers. // // File developed by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory +// Peter W. Doak, doakpw@ornl.gov, Oak Ridge National Laboratory // // File created by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory ////////////////////////////////////////////////////////////////////////////////////// #include "SimulationCell.h" +#include "Lattice/CrystalLattice.h" +#include namespace qmcplusplus { SimulationCell::SimulationCell() = default; -SimulationCell::SimulationCell(const Lattice& lattice) - : lattice_(lattice) -{ - resetLRBox(); -} +template +constexpr auto always_false = false; + +SimulationCell::SimulationCell(const Lattice& lattice) : lattice_(lattice) { resetLRBox(); } void SimulationCell::resetLRBox() { - if (lattice_.SuperCellEnum != SUPERCELL_OPEN) + if (lattice_.getSuperCellEnum() != SUPERCELL_OPEN) { - lattice_.SetLRCutoffs(lattice_.Rv); - LRBox_ = lattice_; + lattice_.SetLRCutoffs(lattice_.getRv()); + lrbox_ = lattice_; bool changed = false; - if (lattice_.SuperCellEnum == SUPERCELL_SLAB && lattice_.VacuumScale != 1.0) + if (lattice_.getSuperCellEnum() == SUPERCELL_SLAB && lattice_.getVacuumScale() != 1.0) { - LRBox_.R(2, 0) *= lattice_.VacuumScale; - LRBox_.R(2, 1) *= lattice_.VacuumScale; - LRBox_.R(2, 2) *= lattice_.VacuumScale; + lrbox_.getR()(2, 0) *= lattice_.getVacuumScale(); + lrbox_.getR()(2, 1) *= lattice_.getVacuumScale(); + lrbox_.getR()(2, 2) *= lattice_.getVacuumScale(); changed = true; } - else if (lattice_.SuperCellEnum == SUPERCELL_WIRE && lattice_.VacuumScale != 1.0) + else if (lattice_.getSuperCellEnum() == SUPERCELL_WIRE && lattice_.getVacuumScale() != 1.0) { - LRBox_.R(1, 0) *= lattice_.VacuumScale; - LRBox_.R(1, 1) *= lattice_.VacuumScale; - LRBox_.R(1, 2) *= lattice_.VacuumScale; - LRBox_.R(2, 0) *= lattice_.VacuumScale; - LRBox_.R(2, 1) *= lattice_.VacuumScale; - LRBox_.R(2, 2) *= lattice_.VacuumScale; + lrbox_.getR()(1, 0) *= lattice_.getVacuumScale(); + lrbox_.getR()(1, 1) *= lattice_.getVacuumScale(); + lrbox_.getR()(1, 2) *= lattice_.getVacuumScale(); + lrbox_.getR()(2, 0) *= lattice_.getVacuumScale(); + lrbox_.getR()(2, 1) *= lattice_.getVacuumScale(); + lrbox_.getR()(2, 2) *= lattice_.getVacuumScale(); changed = true; } - LRBox_.reset(); - LRBox_.SetLRCutoffs(LRBox_.Rv); - LRBox_.printCutoffs(app_log()); + lrbox_.reset(); + lrbox_.SetLRCutoffs(lrbox_.getRv()); + lrbox_.printCutoffs(app_log()); if (changed) { app_summary() << " Simulation box changed by vacuum supercell conditions" << std::endl; app_log() << "--------------------------------------- " << std::endl; - LRBox_.print(app_log()); + lrbox_.print(app_log()); app_log() << "--------------------------------------- " << std::endl; } + lrbox_.reset(); + lrbox_.SetLRCutoffs(lrbox_.getRv()); + lrbox_.printCutoffs(app_log()); - k_lists_.updateKLists(LRBox_, LRBox_.LR_kc, LRBox_.ndim); + k_lists_.updateKLists(lrbox_, lrbox_.LR_kc, lrbox_.ndim); } } -} + +} // namespace qmcplusplus diff --git a/src/Particle/SimulationCell.h b/src/Particle/SimulationCell.h index 52cd2a2951..b1615dcf96 100644 --- a/src/Particle/SimulationCell.h +++ b/src/Particle/SimulationCell.h @@ -2,9 +2,10 @@ // 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. +// Copyright (c) 2025 QMCPACK developers. // // File developed by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory +// Peter W. Doak, doakpw@ornl.gov, Oak Ridge National Laboratory // // File created by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory ////////////////////////////////////////////////////////////////////////////////////// @@ -23,31 +24,31 @@ class ParticleSetPool; class SimulationCell { public: - + using FullPrecReal = QMCTraits::FullPrecRealType; SimulationCell(); SimulationCell(const Lattice& lattice); const Lattice& getLattice() const { return lattice_; } - const Lattice& getPrimLattice() const { return primative_lattice_; } - const Lattice& getLRBox() const { return LRBox_; } - - void resetLRBox(); - - /// access k_lists_ read only + const Lattice& getPrimLattice() const { return primitive_lattice_; } + const Lattice& getLRBox() const { return lrbox_; } const KContainer& getKLists() const { return k_lists_; } + Lattice& getLattice() { return lattice_; } + + void resetLRBox(); private: ///simulation cell lattice Lattice lattice_; ///Primative cell lattice - Lattice primative_lattice_; + Lattice primitive_lattice_; ///long-range box - Lattice LRBox_; + Lattice lrbox_; /// K-Vector List. KContainer k_lists_; friend class ParticleSetPool; }; + } // namespace qmcplusplus #endif From 02e0d5cf997c37d16af885de272f8145d11f0ec9 Mon Sep 17 00:00:00 2001 From: "Peter Doak (epd)" Date: Mon, 24 Mar 2025 15:55:32 -0400 Subject: [PATCH 6/8] return 5382 fix --- src/Particle/LongRange/tests/test_StructFact.cpp | 3 ++- tests/solids/diamondC_1x1x1_pp/CMakeLists.txt | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/Particle/LongRange/tests/test_StructFact.cpp b/src/Particle/LongRange/tests/test_StructFact.cpp index b5f93c57f0..3c91f67f17 100644 --- a/src/Particle/LongRange/tests/test_StructFact.cpp +++ b/src/Particle/LongRange/tests/test_StructFact.cpp @@ -61,7 +61,8 @@ TEST_CASE("StructFact", "[lrhandler]") for (int i = 0; i < ref.groups(); i++) { - std::complex rhok_sum, rhok_even_sum; + std::complex rhok_sum; + std::complex rhok_even_sum; for (int ik = 0; ik < simulation_cell.getKLists().getNumK(); ik++) rhok_sum += std::complex(sk.rhok_r[i][ik], sk.rhok_i[i][ik]); diff --git a/tests/solids/diamondC_1x1x1_pp/CMakeLists.txt b/tests/solids/diamondC_1x1x1_pp/CMakeLists.txt index 611b18fcdb..64ee45dcff 100644 --- a/tests/solids/diamondC_1x1x1_pp/CMakeLists.txt +++ b/tests/solids/diamondC_1x1x1_pp/CMakeLists.txt @@ -1502,7 +1502,7 @@ if(QMC_MIXED_PRECISION) list(APPEND DET_DIAMOND_DMC_MULTI2_SCALARS "ionion" "-12.77566507 0.000001") list(APPEND DET_DIAMOND_DMC_MULTI2_SCALARS "localecp" "-15.41387689 0.00001471") list(APPEND DET_DIAMOND_DMC_MULTI2_SCALARS "nonlocalecp" "6.10588419 0.00001448") - list(APPEND DET_DIAMOND_DMC_MULTI2_SCALARS "mpc" "-2.61905471 0.00000179") + list(APPEND DET_DIAMOND_DMC_MULTI2_SCALARS "mpc" "-2.61905731 0.00000179") list(APPEND DET_DIAMOND_DMC_MULTI2_SCALARS "samples" "9 0.0") # DMC cycle 3 list(APPEND DET_DIAMOND_DMC_MULTI3_SCALARS "totenergy" "-8.93222965 0.00001861") From 5d16623eb18519edf092fb700c6cf6450ef0870b Mon Sep 17 00:00:00 2001 From: "Peter Doak (epd)" Date: Mon, 24 Mar 2025 17:44:20 -0400 Subject: [PATCH 7/8] CI failure resulted in hard look at SimpleParser --- src/Utilities/SimpleParser.cpp | 4 ++- src/Utilities/SimpleParser.h | 4 +++ src/Utilities/tests/test_parser.cpp | 53 +++++++++++++++++++++++++++-- 3 files changed, 57 insertions(+), 4 deletions(-) diff --git a/src/Utilities/SimpleParser.cpp b/src/Utilities/SimpleParser.cpp index 567905b61a..c84bfa50e2 100644 --- a/src/Utilities/SimpleParser.cpp +++ b/src/Utilities/SimpleParser.cpp @@ -27,7 +27,9 @@ char* readLine(char* s, int max, std::istream& fp) { - char ch; + // not initializing this risks in undefined behavior in the case of empty input since + // anything including \n or ; could be in this. + char ch = '\0'; int i = 0; while (fp.get(ch) && !(ch == '\n' || ch == ';')) { diff --git a/src/Utilities/SimpleParser.h b/src/Utilities/SimpleParser.h index 49cf2dbf03..bb3323127d 100644 --- a/src/Utilities/SimpleParser.h +++ b/src/Utilities/SimpleParser.h @@ -26,6 +26,10 @@ #include #include +/** Legacy readLine function, this readline eats the \n and treats ; as \n + * weird side effect if fp is empty and max > 0 then s[0] set to '\0' + * plus some corner cases with line continuation character and \n or ; following + */ char* readLine(char* s, int max, std::istream& fp); int getwords(std::vector& slist, std::istream& fp, int dummy = 0, const std::string& extra_tokens = ""); diff --git a/src/Utilities/tests/test_parser.cpp b/src/Utilities/tests/test_parser.cpp index 8220127005..f972e463a7 100644 --- a/src/Utilities/tests/test_parser.cpp +++ b/src/Utilities/tests/test_parser.cpp @@ -109,7 +109,6 @@ TEST_CASE("readLine", "[utilities]") {"12345678901234567890extra", {"1234567890123456789"}}, // assuming bufLen=20 below }; - for (auto& tc : tlist) { SECTION(string("Parsing string: ") + tc.input) @@ -123,15 +122,63 @@ TEST_CASE("readLine", "[utilities]") REQUIRE(buf == tc.output[i]); if (i == tc.output.size() - 1) { - REQUIRE(out == NULL); + REQUIRE(out == nullptr); } else { - REQUIRE(out != NULL); + REQUIRE(out != nullptr); } } } } + + SECTION("empty input legacy side effect") + { + char buf[1]; + buf[0] = '\n'; + std::istringstream input; + char* out = readLine(buf, 1, input); + CHECK(out == nullptr); + CHECK(buf[0] == '\0'); + } + + SECTION("line continuation character corner cases") + { + { + char buf[1]; + buf[0] = '\n'; + std::istringstream input{"\\"}; + char* out = readLine(buf, 1, input); + CHECK(out == nullptr); + CHECK(buf[0] == '\0'); + } + { + char buf[2]; + buf[0] = '\n'; + std::istringstream input{"\\\n"}; + char* out = readLine(buf, 2, input); + CHECK(out == buf); + CHECK(buf[0] == '\0'); + } + { + char buf[2]; + buf[0] = '\n'; + std::istringstream input{"\\;"}; + char* out = readLine(buf, 2, input); + CHECK(out == buf); + CHECK(buf[0] == '\\'); + CHECK(buf[1] == '\0'); + } + } + SECTION("semicolon corner case") + { + char buf[2]; + buf[0] = '\n'; + std::istringstream input{";"}; + char* out = readLine(buf, 2, input); + CHECK(out == buf); + CHECK(buf[0] == '\0'); + } } From 0c0d197162e91b8806acb300a7121c0fba6fb20f Mon Sep 17 00:00:00 2001 From: "Peter Doak (epd)" Date: Tue, 1 Apr 2025 17:37:32 -0400 Subject: [PATCH 8/8] remove no longer needed include --- src/Particle/LongRange/tests/test_kcontainer.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/Particle/LongRange/tests/test_kcontainer.cpp b/src/Particle/LongRange/tests/test_kcontainer.cpp index 2b31eedafc..a585f159ec 100644 --- a/src/Particle/LongRange/tests/test_kcontainer.cpp +++ b/src/Particle/LongRange/tests/test_kcontainer.cpp @@ -12,7 +12,6 @@ #include "catch.hpp" #include "Particle/LongRange/KContainer.h" -#include "OhmmsPETE/TinyVector.h" #include "Particle/Lattice/CrystalLattice.h" #include "Particle/tests/MinimalParticlePool.h" #include "Utilities/for_testing/checkVector.hpp"