From 47bc03b2939c6e580b377bed3700ad5cb1e0e4d8 Mon Sep 17 00:00:00 2001 From: Philip Fackler Date: Fri, 25 Aug 2023 12:22:16 -0400 Subject: [PATCH] Implement SPOSetBuilderFactoryT and most required builders --- src/QMCWaveFunctions/CMakeLists.txt | 26 +- src/QMCWaveFunctions/CompositeSPOSetT.cpp | 196 +-- src/QMCWaveFunctions/CompositeSPOSetT.h | 180 +-- .../ElectronGas/FreeOrbitalBuilderT.cpp | 108 ++ .../ElectronGas/FreeOrbitalBuilderT.h | 25 + .../LCAO/CuspCorrectionConstructionT.cpp | 801 ++++++++++++ .../LCAO/CuspCorrectionConstructionT.h | 313 +++++ src/QMCWaveFunctions/LCAO/CuspCorrectionT.h | 114 ++ .../LCAO/LCAOrbitalBuilderT.cpp | 1088 +++++++++++++++++ .../LCAO/LCAOrbitalBuilderT.h | 141 +++ src/QMCWaveFunctions/LCAO/LCAOrbitalSetT.h | 660 +++++----- .../LCAO/LCAOrbitalSetWithCorrectionT.h | 3 +- .../LCAO/MultiFunctorAdapter.h | 2 + .../LCAO/SoaCuspCorrectionT.h | 4 +- .../LCAO/SoaLocalizedBasisSet.cpp | 2 +- .../LCAO/SoaLocalizedBasisSet.h | 6 +- .../SPOSetBuilderFactoryT.cpp | 261 ++++ src/QMCWaveFunctions/SPOSetBuilderFactoryT.h | 88 ++ src/QMCWaveFunctions/SPOSetBuilderT.h | 1 + src/QMCWaveFunctions/SPOSetScannerT.h | 216 ++++ 20 files changed, 3755 insertions(+), 480 deletions(-) create mode 100644 src/QMCWaveFunctions/ElectronGas/FreeOrbitalBuilderT.cpp create mode 100644 src/QMCWaveFunctions/ElectronGas/FreeOrbitalBuilderT.h create mode 100644 src/QMCWaveFunctions/LCAO/CuspCorrectionConstructionT.cpp create mode 100644 src/QMCWaveFunctions/LCAO/CuspCorrectionConstructionT.h create mode 100644 src/QMCWaveFunctions/LCAO/CuspCorrectionT.h create mode 100644 src/QMCWaveFunctions/LCAO/LCAOrbitalBuilderT.cpp create mode 100644 src/QMCWaveFunctions/LCAO/LCAOrbitalBuilderT.h create mode 100644 src/QMCWaveFunctions/SPOSetBuilderFactoryT.cpp create mode 100644 src/QMCWaveFunctions/SPOSetBuilderFactoryT.h create mode 100644 src/QMCWaveFunctions/SPOSetScannerT.h diff --git a/src/QMCWaveFunctions/CMakeLists.txt b/src/QMCWaveFunctions/CMakeLists.txt index b1f86e739b..40e8e5a17d 100644 --- a/src/QMCWaveFunctions/CMakeLists.txt +++ b/src/QMCWaveFunctions/CMakeLists.txt @@ -67,21 +67,36 @@ set(JASTROW_SRCS set(JASTROW_OMPTARGET_SRCS Jastrow/TwoBodyJastrow.cpp Jastrow/BsplineFunctor.cpp) -set(FERMION_SRCS ${FERMION_SRCS} ElectronGas/FreeOrbital.cpp ElectronGas/FreeOrbitalT.cpp ElectronGas/FreeOrbitalBuilder.cpp) +set(FERMION_SRCS ${FERMION_SRCS} + ElectronGas/FreeOrbital.cpp + ElectronGas/FreeOrbitalT.cpp + ElectronGas/FreeOrbitalBuilder.cpp + ElectronGas/FreeOrbitalBuilderT.cpp) # wavefunctions only availbale to 3-dim problems if(OHMMS_DIM MATCHES 3) set(JASTROW_SRCS ${JASTROW_SRCS} Jastrow/eeI_JastrowBuilder.cpp Jastrow/CountingJastrowBuilder.cpp) - set(FERMION_SRCS ${FERMION_SRCS} LCAO/LCAOrbitalSet.cpp LCAO/LCAOrbitalSetT.cpp LCAO/LCAOrbitalBuilder.cpp LCAO/MultiQuinticSpline1D.cpp - LCAO/AOBasisBuilder.cpp LCAO/SoaLocalizedBasisSet.cpp) + set(FERMION_SRCS ${FERMION_SRCS} + LCAO/LCAOrbitalSet.cpp + LCAO/LCAOrbitalSetT.cpp + LCAO/LCAOrbitalBuilder.cpp + LCAO/LCAOrbitalBuilderT.cpp + LCAO/MultiQuinticSpline1D.cpp + LCAO/AOBasisBuilder.cpp + LCAO/SoaLocalizedBasisSet.cpp) if(QMC_COMPLEX) set(FERMION_SRCS ${FERMION_SRCS} LCAO/LCAOSpinorBuilder.cpp) else(QMC_COMPLEX) #LCAO cusp correction is not ready for complex - set(FERMION_SRCS ${FERMION_SRCS} LCAO/LCAOrbitalSetWithCorrection.cpp LCAO/LCAOrbitalSetWithCorrectionT.cpp - LCAO/CuspCorrectionConstruction.cpp LCAO/SoaCuspCorrection.cpp LCAO/SoaCuspCorrectionT.cpp) + set(FERMION_SRCS ${FERMION_SRCS} + LCAO/LCAOrbitalSetWithCorrection.cpp + LCAO/LCAOrbitalSetWithCorrectionT.cpp + LCAO/CuspCorrectionConstruction.cpp + LCAO/CuspCorrectionConstructionT.cpp + LCAO/SoaCuspCorrection.cpp + LCAO/SoaCuspCorrectionT.cpp) endif(QMC_COMPLEX) if(HAVE_EINSPLINE) @@ -134,6 +149,7 @@ set(FERMION_SRCS Fermion/DiracDeterminantWithBackflow.cpp Fermion/SlaterDetWithBackflow.cpp SPOSetBuilderFactory.cpp + SPOSetBuilderFactoryT.cpp TrialWaveFunction.cpp TWFdispatcher.cpp TWFFastDerivWrapper.cpp diff --git a/src/QMCWaveFunctions/CompositeSPOSetT.cpp b/src/QMCWaveFunctions/CompositeSPOSetT.cpp index 1d635e8a41..1a0c574b5b 100644 --- a/src/QMCWaveFunctions/CompositeSPOSetT.cpp +++ b/src/QMCWaveFunctions/CompositeSPOSetT.cpp @@ -14,10 +14,9 @@ // at Urbana-Champaign ////////////////////////////////////////////////////////////////////////////////////// -#include "CompositeSPOSetT.h" +#include "QMCWaveFunctions/CompositeSPOSetT.h" #include "OhmmsData/AttributeSet.h" -#include "QMCWaveFunctions/SPOSetBuilderFactory.h" #include "Utilities/IteratorUtility.h" #include @@ -38,27 +37,27 @@ template inline void insert_columns(const MAT1& small, MAT2& big, int offset_c) { - const int c = small.cols(); - for (int i = 0; i < small.rows(); ++i) - std::copy(small[i], small[i] + c, big[i] + offset_c); + const int c = small.cols(); + for (int i = 0; i < small.rows(); ++i) + std::copy(small[i], small[i] + c, big[i] + offset_c); } } // namespace MatrixOperators template CompositeSPOSetT::CompositeSPOSetT(const std::string& my_name) : - SPOSetT(my_name) + SPOSetT(my_name) { - this->OrbitalSetSize = 0; - component_offsets.reserve(4); + this->OrbitalSetSize = 0; + component_offsets.reserve(4); } template CompositeSPOSetT::CompositeSPOSetT(const CompositeSPOSetT& other) : - SPOSetT(other) + SPOSetT(other) { - for (auto& element : other.components) { - this->add(element->makeClone()); - } + for (auto& element : other.components) { + this->add(element->makeClone()); + } } template @@ -68,120 +67,120 @@ template void CompositeSPOSetT::add(std::unique_ptr> component) { - if (components.empty()) - component_offsets.push_back(0); // add 0 + if (components.empty()) + component_offsets.push_back(0); // add 0 - int norbs = component->size(); - components.push_back(std::move(component)); - component_values.emplace_back(norbs); - component_gradients.emplace_back(norbs); - component_laplacians.emplace_back(norbs); + int norbs = component->size(); + components.push_back(std::move(component)); + component_values.emplace_back(norbs); + component_gradients.emplace_back(norbs); + component_laplacians.emplace_back(norbs); - this->OrbitalSetSize += norbs; - component_offsets.push_back(this->OrbitalSetSize); + this->OrbitalSetSize += norbs; + component_offsets.push_back(this->OrbitalSetSize); } template void CompositeSPOSetT::report() { - app_log() << "CompositeSPOSetT" << std::endl; - app_log() << " ncomponents = " << components.size() << std::endl; - app_log() << " components" << std::endl; - for (int i = 0; i < components.size(); ++i) { - app_log() << " " << i << std::endl; - components[i]->basic_report(" "); - } + app_log() << "CompositeSPOSetT" << std::endl; + app_log() << " ncomponents = " << components.size() << std::endl; + app_log() << " components" << std::endl; + for (int i = 0; i < components.size(); ++i) { + app_log() << " " << i << std::endl; + components[i]->basic_report(" "); + } } template std::unique_ptr> CompositeSPOSetT::makeClone() const { - return std::make_unique>(*this); + return std::make_unique>(*this); } template void CompositeSPOSetT::evaluateValue( - const ParticleSet& P, int iat, ValueVector& psi) + const ParticleSet& P, int iat, ValueVector& psi) { - int n = 0; - for (int c = 0; c < components.size(); ++c) { - SPOSetT& component = *components[c]; - ValueVector& values = component_values[c]; - component.evaluateValue(P, iat, values); - std::copy(values.begin(), values.end(), psi.begin() + n); - n += component.size(); - } + int n = 0; + for (int c = 0; c < components.size(); ++c) { + SPOSetT& component = *components[c]; + ValueVector& values = component_values[c]; + component.evaluateValue(P, iat, values); + std::copy(values.begin(), values.end(), psi.begin() + n); + n += component.size(); + } } template void CompositeSPOSetT::evaluateVGL(const ParticleSet& P, int iat, - ValueVector& psi, GradVector& dpsi, ValueVector& d2psi) + ValueVector& psi, GradVector& dpsi, ValueVector& d2psi) { - int n = 0; - for (int c = 0; c < components.size(); ++c) { - SPOSetT& component = *components[c]; - ValueVector& values = component_values[c]; - GradVector& gradients = component_gradients[c]; - ValueVector& laplacians = component_laplacians[c]; - component.evaluateVGL(P, iat, values, gradients, laplacians); - std::copy(values.begin(), values.end(), psi.begin() + n); - std::copy(gradients.begin(), gradients.end(), dpsi.begin() + n); - std::copy(laplacians.begin(), laplacians.end(), d2psi.begin() + n); - n += component.size(); - } + int n = 0; + for (int c = 0; c < components.size(); ++c) { + SPOSetT& component = *components[c]; + ValueVector& values = component_values[c]; + GradVector& gradients = component_gradients[c]; + ValueVector& laplacians = component_laplacians[c]; + component.evaluateVGL(P, iat, values, gradients, laplacians); + std::copy(values.begin(), values.end(), psi.begin() + n); + std::copy(gradients.begin(), gradients.end(), dpsi.begin() + n); + std::copy(laplacians.begin(), laplacians.end(), d2psi.begin() + n); + n += component.size(); + } } template void CompositeSPOSetT::evaluate_notranspose(const ParticleSet& P, int first, - int last, ValueMatrix& logdet, GradMatrix& dlogdet, ValueMatrix& d2logdet) + int last, ValueMatrix& logdet, GradMatrix& dlogdet, ValueMatrix& d2logdet) { - const int nat = last - first; - for (int c = 0; c < components.size(); ++c) { - int norb = components[c]->size(); - ValueMatrix v(nat, norb); - GradMatrix g(nat, norb); - ValueMatrix l(nat, norb); - components[c]->evaluate_notranspose(P, first, last, v, g, l); - int n = component_offsets[c]; - MatrixOperators::insert_columns(v, logdet, n); - MatrixOperators::insert_columns(g, dlogdet, n); - MatrixOperators::insert_columns(l, d2logdet, n); - } + const int nat = last - first; + for (int c = 0; c < components.size(); ++c) { + int norb = components[c]->size(); + ValueMatrix v(nat, norb); + GradMatrix g(nat, norb); + ValueMatrix l(nat, norb); + components[c]->evaluate_notranspose(P, first, last, v, g, l); + int n = component_offsets[c]; + MatrixOperators::insert_columns(v, logdet, n); + MatrixOperators::insert_columns(g, dlogdet, n); + MatrixOperators::insert_columns(l, d2logdet, n); + } } template void CompositeSPOSetT::evaluate_notranspose(const ParticleSet& P, int first, - int last, ValueMatrix& logdet, GradMatrix& dlogdet, - HessMatrix& grad_grad_logdet) + int last, ValueMatrix& logdet, GradMatrix& dlogdet, + HessMatrix& grad_grad_logdet) { - const int nat = last - first; - for (int c = 0; c < components.size(); ++c) { - int norb = components[c]->size(); - ValueMatrix v(nat, norb); - GradMatrix g(nat, norb); - HessMatrix h(nat, norb); - components[c]->evaluate_notranspose(P, first, last, v, g, h); - int n = component_offsets[c]; - MatrixOperators::insert_columns(v, logdet, n); - MatrixOperators::insert_columns(g, dlogdet, n); - MatrixOperators::insert_columns(h, grad_grad_logdet, n); - } + const int nat = last - first; + for (int c = 0; c < components.size(); ++c) { + int norb = components[c]->size(); + ValueMatrix v(nat, norb); + GradMatrix g(nat, norb); + HessMatrix h(nat, norb); + components[c]->evaluate_notranspose(P, first, last, v, g, h); + int n = component_offsets[c]; + MatrixOperators::insert_columns(v, logdet, n); + MatrixOperators::insert_columns(g, dlogdet, n); + MatrixOperators::insert_columns(h, grad_grad_logdet, n); + } } template void CompositeSPOSetT::evaluate_notranspose(const ParticleSet& P, int first, - int last, ValueMatrix& logdet, GradMatrix& dlogdet, - HessMatrix& grad_grad_logdet, GGGMatrix& grad_grad_grad_logdet) + int last, ValueMatrix& logdet, GradMatrix& dlogdet, + HessMatrix& grad_grad_logdet, GGGMatrix& grad_grad_grad_logdet) { - not_implemented( - "evaluate_notranspose(P,first,last,logdet,dlogdet,ddlogdet,dddlogdet)"); + not_implemented( + "evaluate_notranspose(P,first,last,logdet,dlogdet,ddlogdet,dddlogdet)"); } // Class concrete types from ValueType @@ -190,4 +189,37 @@ template class CompositeSPOSetT; template class CompositeSPOSetT>; template class CompositeSPOSetT>; +template +std::unique_ptr> +CompositeSPOSetBuilderT::createSPOSetFromXML(xmlNodePtr cur) +{ + std::vector spolist; + putContent(spolist, cur); + if (spolist.empty()) { + return nullptr; + } + + auto spo_now = std::make_unique>( + getXMLAttributeValue(cur, "name")); + for (int i = 0; i < spolist.size(); ++i) { + const SPOSetT* spo = sposet_builder_factory_.getSPOSet(spolist[i]); + if (spo) + spo_now->add(spo->makeClone()); + } + return (spo_now->size()) ? std::unique_ptr>{std::move(spo_now)} : + nullptr; +} + +template +std::unique_ptr> +CompositeSPOSetBuilderT::createSPOSet(xmlNodePtr cur, SPOSetInputInfo& input) +{ + return createSPOSetFromXML(cur); +} + +template class CompositeSPOSetBuilderT; +template class CompositeSPOSetBuilderT; +template class CompositeSPOSetBuilderT>; +template class CompositeSPOSetBuilderT>; + } // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/CompositeSPOSetT.h b/src/QMCWaveFunctions/CompositeSPOSetT.h index c8d156ac0c..db7344bff9 100644 --- a/src/QMCWaveFunctions/CompositeSPOSetT.h +++ b/src/QMCWaveFunctions/CompositeSPOSetT.h @@ -18,8 +18,8 @@ #define QMCPLUSPLUS_COMPOSITE_SPOSETT_H #include "QMCWaveFunctions/BasisSetBase.h" -#include "QMCWaveFunctions/SPOSetBuilder.h" -#include "QMCWaveFunctions/SPOSetBuilderFactory.h" +#include "QMCWaveFunctions/SPOSetBuilderFactoryT.h" +#include "QMCWaveFunctions/SPOSetBuilderT.h" #include "QMCWaveFunctions/SPOSetT.h" namespace qmcplusplus @@ -28,83 +28,105 @@ template class CompositeSPOSetT : public SPOSetT { public: - using ValueVector = typename SPOSetT::ValueVector; - using GradVector = typename SPOSetT::GradVector; - using ValueMatrix = typename SPOSetT::ValueMatrix; - using GradMatrix = typename SPOSetT::GradMatrix; - using HessMatrix = typename SPOSetT::HessMatrix; - using GGGMatrix = typename SPOSetT::GGGMatrix; - - /// component SPOSets - std::vector>> components; - /// temporary storage for values - std::vector component_values; - /// temporary storage for gradients - std::vector component_gradients; - /// temporary storage for laplacians - std::vector component_laplacians; - /// store the precomputed offsets - std::vector component_offsets; - - CompositeSPOSetT(const std::string& my_name); - /** - * @TODO: do we want template copy constructor - * (i.e., copy from other with different type argument)? - */ - CompositeSPOSetT(const CompositeSPOSetT& other); - ~CompositeSPOSetT() override; - - std::string - getClassName() const override - { - return "CompositeSPOSetT"; - } - - /// add a sposet component to this composite sposet - void - add(std::unique_ptr> component); - - /// print out component info - void - report(); - - // SPOSet interface methods - /// size is determined by component sposets and nothing else - inline void - setOrbitalSetSize(int norbs) override - { - } - - std::unique_ptr> - makeClone() const override; - - void - evaluateValue(const ParticleSet& P, int iat, ValueVector& psi) override; - - void - evaluateVGL(const ParticleSet& P, int iat, ValueVector& psi, - GradVector& dpsi, ValueVector& d2psi) override; - - /// unimplemented functions call this to abort - inline void - not_implemented(const std::string& method) - { - APP_ABORT("CompositeSPOSetT::" + method + " has not been implemented"); - } - - // methods to be implemented in the future (possibly) - void - evaluate_notranspose(const ParticleSet& P, int first, int last, - ValueMatrix& logdet, GradMatrix& dlogdet, - ValueMatrix& d2logdet) override; - void - evaluate_notranspose(const ParticleSet& P, int first, int last, - ValueMatrix& logdet, GradMatrix& dlogdet, - HessMatrix& ddlogdet) override; - void - evaluate_notranspose(const ParticleSet& P, int first, int last, - ValueMatrix& logdet, GradMatrix& dlogdet, HessMatrix& ddlogdet, - GGGMatrix& dddlogdet) override; + using ValueVector = typename SPOSetT::ValueVector; + using GradVector = typename SPOSetT::GradVector; + using ValueMatrix = typename SPOSetT::ValueMatrix; + using GradMatrix = typename SPOSetT::GradMatrix; + using HessMatrix = typename SPOSetT::HessMatrix; + using GGGMatrix = typename SPOSetT::GGGMatrix; + + /// component SPOSets + std::vector>> components; + /// temporary storage for values + std::vector component_values; + /// temporary storage for gradients + std::vector component_gradients; + /// temporary storage for laplacians + std::vector component_laplacians; + /// store the precomputed offsets + std::vector component_offsets; + + CompositeSPOSetT(const std::string& my_name); + /** + * @TODO: do we want template copy constructor + * (i.e., copy from other with different type argument)? + */ + CompositeSPOSetT(const CompositeSPOSetT& other); + ~CompositeSPOSetT() override; + + std::string + getClassName() const override + { + return "CompositeSPOSetT"; + } + + /// add a sposet component to this composite sposet + void + add(std::unique_ptr> component); + + /// print out component info + void + report(); + + // SPOSet interface methods + /// size is determined by component sposets and nothing else + inline void + setOrbitalSetSize(int norbs) override + { + } + + std::unique_ptr> + makeClone() const override; + + void + evaluateValue(const ParticleSet& P, int iat, ValueVector& psi) override; + + void + evaluateVGL(const ParticleSet& P, int iat, ValueVector& psi, + GradVector& dpsi, ValueVector& d2psi) override; + + /// unimplemented functions call this to abort + inline void + not_implemented(const std::string& method) + { + APP_ABORT("CompositeSPOSetT::" + method + " has not been implemented"); + } + + // methods to be implemented in the future (possibly) + void + evaluate_notranspose(const ParticleSet& P, int first, int last, + ValueMatrix& logdet, GradMatrix& dlogdet, + ValueMatrix& d2logdet) override; + void + evaluate_notranspose(const ParticleSet& P, int first, int last, + ValueMatrix& logdet, GradMatrix& dlogdet, + HessMatrix& ddlogdet) override; + void + evaluate_notranspose(const ParticleSet& P, int first, int last, + ValueMatrix& logdet, GradMatrix& dlogdet, HessMatrix& ddlogdet, + GGGMatrix& dddlogdet) override; +}; + +template +class CompositeSPOSetBuilderT : public SPOSetBuilderT +{ +public: + CompositeSPOSetBuilderT( + Communicate* comm, const SPOSetBuilderFactoryT& factory) : + SPOSetBuilderT("Composite", comm), + sposet_builder_factory_(factory) + { + } + + // SPOSetBuilder interface + std::unique_ptr> + createSPOSetFromXML(xmlNodePtr cur) override; + + std::unique_ptr> + createSPOSet(xmlNodePtr cur, SPOSetInputInfo& input) override; + + /// reference to the sposet_builder_factory + const SPOSetBuilderFactoryT& sposet_builder_factory_; }; } // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/ElectronGas/FreeOrbitalBuilderT.cpp b/src/QMCWaveFunctions/ElectronGas/FreeOrbitalBuilderT.cpp new file mode 100644 index 0000000000..4257021557 --- /dev/null +++ b/src/QMCWaveFunctions/ElectronGas/FreeOrbitalBuilderT.cpp @@ -0,0 +1,108 @@ +#include "OhmmsData/AttributeSet.h" +#include "LongRange/StructFact.h" +#include "LongRange/KContainer.h" +#include "QMCWaveFunctions/ElectronGas/FreeOrbitalT.h" +#include "QMCWaveFunctions/ElectronGas/FreeOrbitalBuilderT.h" + +namespace qmcplusplus +{ +template +FreeOrbitalBuilderT::FreeOrbitalBuilderT(ParticleSet& els, Communicate* comm, xmlNodePtr cur) + : SPOSetBuilderT("PW", comm), targetPtcl(els) +{} + +template +std::unique_ptr> FreeOrbitalBuilderT::createSPOSetFromXML(xmlNodePtr cur) +{ + int norb = -1; + std::string spo_object_name; + PosType twist(0.0); + OhmmsAttributeSet attrib; + attrib.add(norb, "size"); + attrib.add(twist, "twist"); + attrib.add(spo_object_name, "name"); + attrib.put(cur); + + if (norb < 0) + throw std::runtime_error("free orbital SPO set require the \"size\" input"); + + auto lattice = targetPtcl.getLattice(); + + PosType tvec = lattice.k_cart(twist); +#ifdef QMC_COMPLEX + const int npw = norb; + targetPtcl.setTwist(twist); + app_log() << "twist fraction = " << twist << std::endl; + app_log() << "twist cartesian = " << tvec << std::endl; +#else + const int npw = std::ceil((norb + 1.0) / 2); + if (2 * npw - 1 != norb) + { + std::ostringstream msg; + msg << "norb = " << norb << " npw = " << npw; + msg << " cannot be ran in real PWs (sin, cos)" << std::endl; + msg << "either use complex build or change the size of SPO set" << std::endl; + msg << "ideally, set size to a closed shell of PWs." << std::endl; + throw std::runtime_error(msg.str()); + } + for (int ldim = 0; ldim < twist.size(); ldim++) + { + if (std::abs(twist[ldim]) > 1e-16) + throw std::runtime_error("no twist for real orbitals"); + } +#endif + + // extract npw k-points from container + // kpts_cart is sorted by magnitude + std::vector kpts(npw); + KContainer klists; + RealType kcut = lattice.LR_kc; // to-do: reduce kcut to >~ kf + klists.updateKLists(lattice, kcut, lattice.ndim, twist); + + // k0 is not in kpts_cart + kpts[0] = tvec; +#ifdef QMC_COMPLEX + for (int ik = 1; ik < npw; ik++) + { + kpts[ik] = klists.kpts_cart[ik - 1]; + } +#else + const int nktot = klists.kpts.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]; + if (in_list(jk, mkidx)) + continue; + // if not, then add this kpoint + kpts[ik] = klists.kpts_cart[jk]; + mkidx[ik] = jmk; // keep track of its minus + ik++; + if (ik >= npw) + break; + } +#endif + auto sposet = std::make_unique>(spo_object_name, kpts); + sposet->report(" "); + return sposet; +} + +template +bool FreeOrbitalBuilderT::in_list(const int j, const std::vector l) +{ + for (int i = 0; i < l.size(); i++) + { + if (j == l[i]) + return true; + } + return false; +} + +template class FreeOrbitalBuilderT; +template class FreeOrbitalBuilderT; +template class FreeOrbitalBuilderT>; +template class FreeOrbitalBuilderT>; + +} // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/ElectronGas/FreeOrbitalBuilderT.h b/src/QMCWaveFunctions/ElectronGas/FreeOrbitalBuilderT.h new file mode 100644 index 0000000000..dcd69fd4b8 --- /dev/null +++ b/src/QMCWaveFunctions/ElectronGas/FreeOrbitalBuilderT.h @@ -0,0 +1,25 @@ +#ifndef QMCPLUSPLUS_FREE_ORBITAL_BUILDERT_H +#define QMCPLUSPLUS_FREE_ORBITAL_BUILDERT_H + +#include "QMCWaveFunctions/SPOSetBuilderT.h" + +namespace qmcplusplus +{ +template +class FreeOrbitalBuilderT : public SPOSetBuilderT +{ +public: + using RealType = typename SPOSetBuilderT::RealType; + using PosType = typename SPOSetBuilderT::PosType; + + FreeOrbitalBuilderT(ParticleSet& els, Communicate* comm, xmlNodePtr cur); + ~FreeOrbitalBuilderT() {} + + std::unique_ptr> createSPOSetFromXML(xmlNodePtr cur) override; + +private: + ParticleSet& targetPtcl; + bool in_list(const int j, const std::vector l); +}; +} // namespace qmcplusplus +#endif diff --git a/src/QMCWaveFunctions/LCAO/CuspCorrectionConstructionT.cpp b/src/QMCWaveFunctions/LCAO/CuspCorrectionConstructionT.cpp new file mode 100644 index 0000000000..1178491533 --- /dev/null +++ b/src/QMCWaveFunctions/LCAO/CuspCorrectionConstructionT.cpp @@ -0,0 +1,801 @@ +////////////////////////////////////////////////////////////////////////////////////// +// This file is distributed under the University of Illinois/NCSA Open Source +// License. See LICENSE file in top directory for details. +// +// Copyright (c) 2023 QMCPACK developers. +// +// File developed by: Mark Dewing, mdewing@anl.gov, Argonne National Laboratory +// Ye Luo, yeluo@anl.gov, Argonne National Laboratory +// Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore +// National Laboratory Jeremy McMinnis, jmcminis@gmail.com, +// University of Illinois at Urbana-Champaign Mark A. +// Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory +// +// File created by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory +////////////////////////////////////////////////////////////////////////////////////// + +#include "CuspCorrectionConstructionT.h" + +#include "Message/Communicate.h" +#include "MultiQuinticSpline1D.h" +#include "Numerics/MinimizeOneDim.h" +#include "OhmmsData/AttributeSet.h" +#include "SoaAtomicBasisSet.h" +#include "SoaCuspCorrectionBasisSet.h" +#include "SoaLocalizedBasisSet.h" +#include "Utilities/FairDivide.h" + +namespace qmcplusplus +{ +template +void +CuspCorrectionConstructionT::splitPhiEta(int center, + const std::vector& corrCenter, LCAOrbitalSetT& Phi, + LCAOrbitalSetT& Eta) +{ + std::vector is_s_orbital(Phi.myBasisSet->BasisSetSize, false); + std::vector correct_this_center(corrCenter.size(), false); + correct_this_center[center] = corrCenter[center]; + + Phi.myBasisSet->queryOrbitalsForSType(correct_this_center, is_s_orbital); + + int nOrbs = Phi.getOrbitalSetSize(); + int bss = Phi.getBasisSetSize(); + + for (int i = 0; i < bss; i++) { + if (is_s_orbital[i]) { + auto& cref(*(Eta.C)); + for (int k = 0; k < nOrbs; k++) + cref(k, i) = 0.0; // Eta->C(k,i) = 0.0; + } + else { + auto& cref(*(Phi.C)); + for (int k = 0; k < nOrbs; k++) + cref(k, i) = 0.0; // Phi->C(k,i) = 0.0; + } + } +} + +template +void +CuspCorrectionConstructionT::removeSTypeOrbitals( + const std::vector& corrCenter, LCAOrbitalSetT& Phi) +{ + std::vector is_s_orbital(Phi.myBasisSet->BasisSetSize, false); + + Phi.myBasisSet->queryOrbitalsForSType(corrCenter, is_s_orbital); + + int nOrbs = Phi.getOrbitalSetSize(); + int bss = Phi.getBasisSetSize(); + + for (int i = 0; i < bss; i++) { + if (is_s_orbital[i]) { + auto& cref(*(Phi.C)); + for (int k = 0; k < nOrbs; k++) + cref(k, i) = 0.0; + } + } +} + +// Will be the corrected value for r < rc and the original wavefunction for r > +// rc +template +void +CuspCorrectionConstructionT::computeRadialPhiBar(ParticleSet* targetP, + ParticleSet* sourceP, int curOrb_, int curCenter_, SPOSetT* Phi, + Vector& xgrid, Vector& rad_orb, + const CuspCorrectionParametersT& data) +{ + OneMolecularOrbitalT phiMO(targetP, sourceP, Phi); + phiMO.changeOrbital(curCenter_, curOrb_); + CuspCorrectionT cusp(data); + + for (int i = 0; i < xgrid.size(); i++) { + rad_orb[i] = phiBar(cusp, xgrid[i], phiMO); + } +} + +// Get the ideal local energy at one point +// Eq. 17 in the paper. Coefficients are taken from the paper. +template +typename CuspCorrectionConstructionT::RealType +CuspCorrectionConstructionT::getOneIdealLocalEnergy( + RealType r, RealType Z, RealType beta0) +{ + RealType beta[7] = { + 3.25819, -15.0126, 33.7308, -42.8705, 31.2276, -12.1316, 1.94692}; + RealType idealEL = beta0; + RealType r1 = r * r; + for (int i = 0; i < 7; i++) { + idealEL += beta[i] * r1; + r1 *= r; + } + return idealEL * Z * Z; +} + +// Get the ideal local energy for a vector of positions +template +void +CuspCorrectionConstructionT::getIdealLocalEnergy(const ValueVector& pos, + RealType Z, RealType Rc, RealType ELorigAtRc, ValueVector& ELideal) +{ + // assert(pos.size() == ELideal.size() + RealType beta0 = 0.0; + RealType tmp = getOneIdealLocalEnergy(Rc, Z, beta0); + beta0 = (ELorigAtRc - tmp) / (Z * Z); + for (int i = 0; i < pos.size(); i++) { + ELideal[i] = getOneIdealLocalEnergy(pos[i], Z, beta0); + } +} + +// Evaluate constraints. Equations 9-13 in the paper. +template +void +CuspCorrectionConstructionT::evalX(RealType valRc, GradType gradRc, + ValueType lapRc, RealType Rc, RealType Z, RealType C, RealType valAtZero, + RealType eta0, TinyVector& X) +{ + X[0] = std::log(std::abs(valRc - C)); + X[1] = gradRc[0] / (valRc - C); + X[2] = (lapRc - 2.0 * gradRc[0] / Rc) / (valRc - C); + X[3] = -Z * (valAtZero + eta0) / (valAtZero - C); + X[4] = std::log(std::abs(valAtZero - C)); +} + +// Compute polynomial coefficients from constraints. Eq. 14 in the paper. +template +void +CuspCorrectionConstructionT::X2alpha(const TinyVector& X, + RealType Rc, TinyVector& alpha) +{ + RealType RcInv = 1.0 / Rc, RcInv2 = RcInv * RcInv; + alpha[0] = X[4]; + alpha[1] = X[3]; + alpha[2] = 6.0 * X[0] * RcInv2 - 3.0 * X[1] * RcInv + X[2] * 0.5 - + 3.0 * X[3] * RcInv - 6.0 * X[4] * RcInv2 - 0.5 * X[1] * X[1]; + alpha[3] = -8.0 * X[0] * RcInv2 * RcInv + 5.0 * X[1] * RcInv2 - + X[2] * RcInv + 3.0 * X[3] * RcInv2 + 8.0 * X[4] * RcInv2 * RcInv + + X[1] * X[1] * RcInv; + alpha[4] = 3.0 * X[0] * RcInv2 * RcInv2 - 2.0 * X[1] * RcInv2 * RcInv + + 0.5 * X[2] * RcInv2 - X[3] * RcInv2 * RcInv - + 3.0 * X[4] * RcInv2 * RcInv2 - 0.5 * X[1] * X[1] * RcInv2; +} + +// Eq. 16 in the paper. +template +typename CuspCorrectionConstructionT::RealType +CuspCorrectionConstructionT::getZeff( + RealType Z, RealType etaAtZero, RealType phiBarAtZero) +{ + return Z * (1.0 + etaAtZero / phiBarAtZero); +} + +template +typename CuspCorrectionConstructionT::RealType +CuspCorrectionConstructionT::phiBar( + const CuspCorrectionT& cusp, RealType r, OneMolecularOrbitalT& phiMO) +{ + if (r <= cusp.cparam.Rc) + return cusp.cparam.C + cusp.Rr(r); + else + return phiMO.phi(r); +} + +// Compute the effective one-electron local energy at a vector of points. +// Eq. 15 in the paper for r < Rc. Normal local energy for R > Rc. +template +void +CuspCorrectionConstructionT::getCurrentLocalEnergy(const ValueVector& pos, + RealType Zeff, RealType Rc, RealType originalELatRc, + CuspCorrectionT& cusp, OneMolecularOrbitalT& phiMO, + ValueVector& ELcurr) +{ + // assert(pos.size() == ELcurr.size()); + ValueType val; + GradType grad; + ValueType lap; + phiMO.phi_vgl(Rc, val, grad, lap); + RealType dE = originalELatRc - (-0.5 * lap / val - Zeff / Rc); + for (int i = 0; i < pos.size(); i++) { + RealType r = pos[i]; + // prevent NaN's if phiBar is zero + RealType offset = 1e-12; + if (r <= Rc) { + RealType dp = cusp.dpr(r); + ELcurr[i] = -0.5 * cusp.Rr(r) * + (2.0 * dp / r + cusp.d2pr(r) + dp * dp) / + (offset + phiBar(cusp, r, phiMO)) - + Zeff / r + dE; + } + else { + phiMO.phi_vgl(pos[i], val, grad, lap); + ELcurr[i] = -0.5 * lap / val - Zeff / r + dE; + } + } +} + +// Return value is local energy at Rc +template +typename CuspCorrectionConstructionT::RealType +CuspCorrectionConstructionT::getOriginalLocalEnergy(const ValueVector& pos, + RealType Zeff, RealType Rc, OneMolecularOrbitalT& phiMO, + ValueVector& ELorig) +{ + // assert(pos.size() == ELorig.size()); + + ValueType val; + GradType grad; + ValueType lap; + for (int i = 0; i < pos.size(); i++) { + RealType r = pos[i]; + phiMO.phi_vgl(r, val, grad, lap); + ELorig[i] = -0.5 * lap / val - Zeff / r; + } + + phiMO.phi_vgl(Rc, val, grad, lap); + return -0.5 * lap / val - Zeff / Rc; +} + +// Sum of squares difference between the current local energy and the ideal +// local energy. +// This is the objective function to minimize. +template +typename CuspCorrectionConstructionT::RealType +CuspCorrectionConstructionT::getELchi2( + const ValueVector& ELcurr, const ValueVector& ELideal) +{ + assert(ELcurr.size() == ELideal.size()); + + RealType chi2 = 0.0; + for (int i = 0; i < ELcurr.size(); i++) { + RealType diff = ELcurr[i] - ELideal[i]; + chi2 += diff * diff; + } + return chi2; +} + +// Compute the chi squared distance given a value for phi at zero. +template +typename CuspCorrectionConstructionT::RealType +CuspCorrectionConstructionT::evaluateForPhi0Body(RealType phi0, + ValueVector& pos, ValueVector& ELcurr, ValueVector& ELideal, + CuspCorrectionT& cusp, OneMolecularOrbitalT& phiMO, + ValGradLap phiAtRc, RealType etaAtZero, RealType ELorigAtRc, RealType Z) +{ + cusp.cparam.sg = phi0 > 0.0 ? 1.0 : -1.0; + cusp.cparam.C = (phiAtRc.val * phi0 < 0.0) ? 1.5 * phiAtRc.val : 0.0; + TinyVector X; + evalX(phiAtRc.val, phiAtRc.grad, phiAtRc.lap, cusp.cparam.Rc, Z, + cusp.cparam.C, phi0, etaAtZero, X); + X2alpha(X, cusp.cparam.Rc, cusp.cparam.alpha); + RealType Zeff = getZeff(Z, etaAtZero, phiBar(cusp, 0.0, phiMO)); + getCurrentLocalEnergy( + pos, Zeff, cusp.cparam.Rc, ELorigAtRc, cusp, phiMO, ELcurr); + RealType chi2 = getELchi2(ELcurr, ELideal); + return chi2; +} + +// Optimize free parameter (value of phi at zero) to minimize distance to ideal +// local energy. Output is return value and parameter values are in cusp.cparam +template +typename CuspCorrectionConstructionT::RealType +CuspCorrectionConstructionT::minimizeForPhiAtZero(CuspCorrectionT& cusp, + OneMolecularOrbitalT& phiMO, RealType Z, RealType eta0, ValueVector& pos, + ValueVector& ELcurr, ValueVector& ELideal, RealType start_phi0) +{ + ValGradLap vglAtRc; + ValueVector tmp_pos(0); + ValueVector ELorig(0); + RealType Zeff = getZeff(Z, eta0, phiBar(cusp, 0.0, phiMO)); + + RealType ELorigAtRc = + getOriginalLocalEnergy(tmp_pos, Zeff, cusp.cparam.Rc, phiMO, ELorig); + getIdealLocalEnergy(pos, Z, cusp.cparam.Rc, ELorigAtRc, ELideal); + phiMO.phi_vgl(cusp.cparam.Rc, vglAtRc.val, vglAtRc.grad, vglAtRc.lap); + + Bracket_min_t bracket(start_phi0, 0.0, 0.0, false); + try { + bracket = bracket_minimum( + [&](RealType x) -> RealType { + return evaluateForPhi0Body(x, pos, ELcurr, ELideal, cusp, phiMO, + vglAtRc, eta0, ELorigAtRc, Z); + }, + start_phi0); + } + catch (const std::runtime_error& e) { + APP_ABORT("Bracketing minimum failed for finding phi0. \n"); + } + + auto min_res = find_minimum( + [&](RealType x) -> RealType { + return evaluateForPhi0Body(x, pos, ELcurr, ELideal, cusp, phiMO, + vglAtRc, eta0, ELorigAtRc, Z); + }, + bracket); + + start_phi0 = min_res.first; + + return min_res.second; +} + +// Optimize the cutoff radius. There is an inner loop optimizing for phi0 for +// each value of Rc. Elcurr and ELideal are expected to have the correct size on +// input (same size as pos) Output is parameter values in cusp.cparam +template +void +CuspCorrectionConstructionT::minimizeForRc(CuspCorrectionT& cusp, + OneMolecularOrbitalT& phiMO, RealType Z, RealType Rc_init, + RealType Rc_max, RealType eta0, ValueVector& pos, ValueVector& ELcurr, + ValueVector& ELideal) +{ + Bracket_min_t bracket(Rc_init, 0.0, 0.0, false); + RealType start_phi0 = phiMO.phi(0.0); + try { + bracket = bracket_minimum( + [&](RealType x) -> RealType { + cusp.cparam.Rc = x; + return minimizeForPhiAtZero( + cusp, phiMO, Z, eta0, pos, ELcurr, ELideal, start_phi0); + }, + Rc_init, Rc_max); + } + catch (const std::runtime_error& e) { + APP_ABORT("Bracketing minimum failed for finding rc. \n"); + } + + if (bracket.success) { + auto min_res = find_minimum( + [&](RealType x) -> RealType { + cusp.cparam.Rc = x; + return minimizeForPhiAtZero( + cusp, phiMO, Z, eta0, pos, ELcurr, ELideal, start_phi0); + }, + bracket); + } + else { + cusp.cparam.Rc = bracket.a; + minimizeForPhiAtZero( + cusp, phiMO, Z, eta0, pos, ELcurr, ELideal, start_phi0); + } +} + +// Modifies orbital set lcwc +template +void +CuspCorrectionConstructionT::applyCuspCorrection( + const Matrix>& info, ParticleSet& targetPtcl, + ParticleSet& sourcePtcl, LCAOrbitalSetT& lcao, + SoaCuspCorrectionT& cusp, const std::string& id) +{ + const int num_centers = info.rows(); + const int orbital_set_size = info.cols(); + using RealType = typename SPOSetT::RealType; + + NewTimer& cuspApplyTimer = createGlobalTimer( + "CuspCorrectionConstruction::applyCuspCorrection", timer_level_medium); + + ScopedTimer cuspApplyTimerWrapper(cuspApplyTimer); + + LCAOrbitalSetT phi("phi", + std::unique_ptr::basis_type>( + lcao.myBasisSet->makeClone())); + phi.setOrbitalSetSize(lcao.getOrbitalSetSize()); + + LCAOrbitalSetT eta("eta", + std::unique_ptr::basis_type>( + lcao.myBasisSet->makeClone())); + eta.setOrbitalSetSize(lcao.getOrbitalSetSize()); + + std::vector corrCenter(num_centers, "true"); + + // What's this grid's lifespan? Why on the heap? + auto radial_grid = std::make_unique>(); + radial_grid->set(0.000001, 100.0, 1001); + + Vector xgrid; + Vector rad_orb; + xgrid.resize(radial_grid->size()); + rad_orb.resize(radial_grid->size()); + for (int ig = 0; ig < radial_grid->size(); ig++) { + xgrid[ig] = radial_grid->r(ig); + } + + for (int ic = 0; ic < num_centers; ic++) { + *eta.C = *lcao.C; + *phi.C = *lcao.C; + + splitPhiEta(ic, corrCenter, phi, eta); + + // loop over MO index - cot must be an array (of len MO size) + // the loop is inside cot - in the multiqunitic + auto cot = std::make_unique>(); + cot->initializeRadialSet(*radial_grid, orbital_set_size); + // How is this useful? + // cot->ID.resize(orbital_set_size); + // for (int mo_idx = 0; mo_idx < orbital_set_size; mo_idx++) { + // cot->ID[mo_idx] = mo_idx; + // } + + for (int mo_idx = 0; mo_idx < orbital_set_size; mo_idx++) { + computeRadialPhiBar(&targetPtcl, &sourcePtcl, mo_idx, ic, &phi, + xgrid, rad_orb, info(ic, mo_idx)); + RealType yprime_i = (rad_orb[1] - rad_orb[0]) / + (radial_grid->r(1) - radial_grid->r(0)); + OneDimQuinticSpline radial_spline( + radial_grid->makeClone(), rad_orb); + radial_spline.spline(0, yprime_i, rad_orb.size() - 1, 0.0); + cot->addSpline(mo_idx, radial_spline); + + if (outputManager.isDebugActive()) { + // For testing against AoS output + // Output phiBar to soaOrbs.downdet.C0.MO0 + int nElms = 500; + RealType dx = info(ic, mo_idx).Rc * 1.2 / nElms; + Vector pos; + Vector output_orb; + pos.resize(nElms); + output_orb.resize(nElms); + for (int i = 0; i < nElms; i++) { + pos[i] = (i + 1.0) * dx; + } + computeRadialPhiBar(&targetPtcl, &sourcePtcl, mo_idx, ic, &phi, + pos, output_orb, info(ic, mo_idx)); + std::string filename = "soaOrbs." + id + ".C" + + std::to_string(ic) + ".MO" + std::to_string(mo_idx); + std::cout << "Writing to " << filename << std::endl; + std::ofstream out(filename.c_str()); + out << "# r phiBar(r)" << std::endl; + for (int i = 0; i < nElms; i++) { + out << pos[i] << " " << output_orb[i] << std::endl; + } + out.close(); + } + } + cusp.add(ic, std::move(cot)); + } + removeSTypeOrbitals(corrCenter, lcao); +} + +template +void +CuspCorrectionConstructionT::generateCuspInfo( + Matrix>& info, const ParticleSet& targetPtcl, + const ParticleSet& sourcePtcl, const LCAOrbitalSetT& lcao, + const std::string& id, Communicate& Comm) +{ + const int num_centers = info.rows(); + const int orbital_set_size = info.cols(); + using RealType = typename SPOSetT::RealType; + using ValueVector = typename SPOSetT::ValueVector; + + NewTimer& cuspCreateTimer = createGlobalTimer( + "CuspCorrectionConstruction::createCuspParameters", timer_level_medium); + NewTimer& splitPhiEtaTimer = createGlobalTimer( + "CuspCorrectionConstruction::splitPhiEta", timer_level_fine); + NewTimer& computeTimer = createGlobalTimer( + "CuspCorrectionConstruction::computeCorrection", timer_level_fine); + + ScopedTimer createCuspTimerWrapper(cuspCreateTimer); + + LCAOrbitalSetT phi("phi", + std::unique_ptr::basis_type>( + lcao.myBasisSet->makeClone())); + phi.setOrbitalSetSize(lcao.getOrbitalSetSize()); + + LCAOrbitalSetT eta("eta", + std::unique_ptr::basis_type>( + lcao.myBasisSet->makeClone())); + eta.setOrbitalSetSize(lcao.getOrbitalSetSize()); + + std::vector corrCenter(num_centers, "true"); + + using GridType = OneDimGridBase; + int npts = 500; + + // Parallelize correction of MO's across MPI ranks + std::vector offset; + FairDivideLow(orbital_set_size, Comm.size(), offset); + + int start_mo = offset[Comm.rank()]; + int end_mo = offset[Comm.rank() + 1]; + app_log() + << " Number of molecular orbitals to compute correction on this rank: " + << end_mo - start_mo << std::endl; + +// Specify dynamic scheduling explicitly for load balancing. Each iteration +// should take enough time that scheduling overhead is not an issue. +#pragma omp parallel for schedule(dynamic) collapse(2) + for (int center_idx = 0; center_idx < num_centers; center_idx++) { + for (int mo_idx = start_mo; mo_idx < end_mo; mo_idx++) { + ParticleSet localTargetPtcl(targetPtcl); + ParticleSet localSourcePtcl(sourcePtcl); + + LCAOrbitalSetT local_phi("local_phi", + std::unique_ptr::basis_type>( + phi.myBasisSet->makeClone())); + local_phi.setOrbitalSetSize(phi.getOrbitalSetSize()); + + LCAOrbitalSetT local_eta("local_eta", + std::unique_ptr::basis_type>( + eta.myBasisSet->makeClone())); + local_eta.setOrbitalSetSize(eta.getOrbitalSetSize()); + +#pragma omp critical + app_log() << " Working on MO: " << mo_idx + << " Center: " << center_idx << std::endl; + + { + ScopedTimer local_timer(splitPhiEtaTimer); + + *local_eta.C = *lcao.C; + *local_phi.C = *lcao.C; + splitPhiEta(center_idx, corrCenter, local_phi, local_eta); + } + + bool corrO = false; + auto& cref(*(local_phi.C)); + for (int ip = 0; ip < cref.cols(); ip++) { + if (std::abs(cref(mo_idx, ip)) > 0) { + corrO = true; + break; + } + } + + if (corrO) { + OneMolecularOrbitalT etaMO( + &localTargetPtcl, &localSourcePtcl, &local_eta); + etaMO.changeOrbital(center_idx, mo_idx); + + OneMolecularOrbitalT phiMO( + &localTargetPtcl, &localSourcePtcl, &local_phi); + phiMO.changeOrbital(center_idx, mo_idx); + + SpeciesSet& tspecies(localSourcePtcl.getSpeciesSet()); + int iz = tspecies.addAttribute("charge"); + RealType Z = tspecies(iz, localSourcePtcl.GroupID[center_idx]); + + RealType Rc_max = 0.2; + RealType rc = 0.1; + + RealType dx = rc * 1.2 / npts; + ValueVector pos(npts); + ValueVector ELideal(npts); + ValueVector ELcurr(npts); + for (int i = 0; i < npts; i++) { + pos[i] = (i + 1.0) * dx; + } + + RealType eta0 = etaMO.phi(0.0); + ValueVector ELorig(npts); + CuspCorrectionT cusp(info(center_idx, mo_idx)); + { + ScopedTimer local_timer(computeTimer); + minimizeForRc( + cusp, phiMO, Z, rc, Rc_max, eta0, pos, ELcurr, ELideal); + } + // Update shared object. Each iteration accesses a different + // element and this is an array (no bookkeeping data to update), + // so no synchronization is necessary. + info(center_idx, mo_idx) = cusp.cparam; + } + } + } + + for (int root = 0; root < Comm.size(); root++) { + int start_mo = offset[root]; + int end_mo = offset[root + 1]; + for (int mo_idx = start_mo; mo_idx < end_mo; mo_idx++) { + for (int center_idx = 0; center_idx < num_centers; center_idx++) { + broadcastCuspInfo(info(center_idx, mo_idx), Comm, root); + } + } + } +} + +template +void +CuspCorrectionConstructionT::broadcastCuspInfo( + CuspCorrectionParametersT& param, Communicate& Comm, int root) +{ +#ifdef HAVE_MPI + std::vector buffer(9); + buffer[0] = param.Rc; + buffer[1] = param.C; + buffer[2] = param.sg; + buffer[3] = param.alpha[0]; + buffer[4] = param.alpha[1]; + buffer[5] = param.alpha[2]; + buffer[6] = param.alpha[3]; + buffer[7] = param.alpha[4]; + buffer[8] = param.redo; + + Comm.comm.broadcast(buffer.begin(), buffer.end(), root); + + param.Rc = buffer[0]; + param.C = buffer[1]; + param.sg = buffer[2]; + param.alpha[0] = buffer[3]; + param.alpha[1] = buffer[4]; + param.alpha[2] = buffer[5]; + param.alpha[3] = buffer[6]; + param.alpha[4] = buffer[7]; + param.redo = buffer[8] == 0.0 ? 0 : 1; +#endif +} + +template +bool +CuspCorrectionConstructionT::readCuspInfo(const std::string& cuspInfoFile, + const std::string& objectName, int OrbitalSetSize, + Matrix>& info) +{ + bool success = true; + int ncenter = info.rows(); + app_log() << "Reading cusp info from : " << cuspInfoFile << std::endl; + Libxml2Document adoc; + if (!adoc.parse(cuspInfoFile)) { + app_log() << "Could not find precomputed cusp data for spo set: " + << objectName << std::endl; + app_log() << "Recalculating data.\n"; + return false; + } + xmlNodePtr head = adoc.getRoot(); + head = head->children; + xmlNodePtr cur = NULL, ctr; + while (head != NULL) { + std::string cname(getNodeName(head)); + if (cname == "sposet") { + std::string name; + OhmmsAttributeSet spoAttrib; + spoAttrib.add(name, "name"); + spoAttrib.put(head); + if (name == objectName) { + cur = head; + break; + } + } + head = head->next; + } + if (cur == NULL) { + app_log() << "Could not find precomputed cusp data for spo set: " + << objectName << std::endl; + app_log() << "Recalculating data.\n"; + return false; + } + else { + app_log() << "Found precomputed cusp data for spo set: " << objectName + << std::endl; + } + cur = cur->children; + while (cur != NULL) { + std::string cname(getNodeName(cur)); + if (cname == "center") { + int num = -1; + OhmmsAttributeSet Attrib; + Attrib.add(num, "num"); + Attrib.put(cur); + if (num < 0 || num >= ncenter) { + APP_ABORT("Error with cusp info xml block. incorrect center " + "number. \n"); + } + ctr = cur->children; + while (ctr != NULL) { + std::string cname(getNodeName(ctr)); + if (cname == "orbital") { + int orb = -1; + OhmmsAttributeSet orbAttrib; + QMCTraits::RealType a1(0.0), a2, a3, a4, a5, a6, a7, a8, a9; + orbAttrib.add(orb, "num"); + orbAttrib.add(a1, "redo"); + orbAttrib.add(a2, "C"); + orbAttrib.add(a3, "sg"); + orbAttrib.add(a4, "rc"); + orbAttrib.add(a5, "a1"); + orbAttrib.add(a6, "a2"); + orbAttrib.add(a7, "a3"); + orbAttrib.add(a8, "a4"); + orbAttrib.add(a9, "a5"); + orbAttrib.put(ctr); + if (orb < OrbitalSetSize) { + info(num, orb).redo = a1; + info(num, orb).C = a2; + info(num, orb).sg = a3; + info(num, orb).Rc = a4; + info(num, orb).alpha[0] = a5; + info(num, orb).alpha[1] = a6; + info(num, orb).alpha[2] = a7; + info(num, orb).alpha[3] = a8; + info(num, orb).alpha[4] = a9; + } + } + ctr = ctr->next; + } + } + cur = cur->next; + } + return success; +} + +template +void +CuspCorrectionConstructionT::saveCusp(const std::string& filename, + const Matrix>& info, const std::string& id) +{ + const int num_centers = info.rows(); + const int orbital_set_size = info.cols(); + xmlDocPtr doc = xmlNewDoc((const xmlChar*)"1.0"); + xmlNodePtr cuspRoot = xmlNewNode(NULL, BAD_CAST "qmcsystem"); + xmlNodePtr spo = xmlNewNode(NULL, (const xmlChar*)"sposet"); + xmlNewProp(spo, (const xmlChar*)"name", (const xmlChar*)id.c_str()); + xmlAddChild(cuspRoot, spo); + xmlDocSetRootElement(doc, cuspRoot); + + for (int center_idx = 0; center_idx < num_centers; center_idx++) { + xmlNodePtr ctr = xmlNewNode(NULL, (const xmlChar*)"center"); + std::ostringstream num; + num << center_idx; + xmlNewProp( + ctr, (const xmlChar*)"num", (const xmlChar*)num.str().c_str()); + + for (int mo_idx = 0; mo_idx < orbital_set_size; mo_idx++) { + std::ostringstream num0, C, sg, rc, a1, a2, a3, a4, a5; + xmlNodePtr orb = xmlNewNode(NULL, (const xmlChar*)"orbital"); + num0 << mo_idx; + xmlNewProp( + orb, (const xmlChar*)"num", (const xmlChar*)num0.str().c_str()); + + C.setf(std::ios::scientific, std::ios::floatfield); + C.precision(14); + C << info(center_idx, mo_idx).C; + sg.setf(std::ios::scientific, std::ios::floatfield); + sg.precision(14); + sg << info(center_idx, mo_idx).sg; + rc.setf(std::ios::scientific, std::ios::floatfield); + rc.precision(14); + rc << info(center_idx, mo_idx).Rc; + a1.setf(std::ios::scientific, std::ios::floatfield); + a1.precision(14); + a1 << info(center_idx, mo_idx).alpha[0]; + a2.setf(std::ios::scientific, std::ios::floatfield); + a2.precision(14); + a2 << info(center_idx, mo_idx).alpha[1]; + a3.setf(std::ios::scientific, std::ios::floatfield); + a3.precision(14); + a3 << info(center_idx, mo_idx).alpha[2]; + a4.setf(std::ios::scientific, std::ios::floatfield); + a4.precision(14); + a4 << info(center_idx, mo_idx).alpha[3]; + a5.setf(std::ios::scientific, std::ios::floatfield); + a5.precision(14); + a5 << info(center_idx, mo_idx).alpha[4]; + xmlNewProp( + orb, (const xmlChar*)"C", (const xmlChar*)C.str().c_str()); + xmlNewProp( + orb, (const xmlChar*)"sg", (const xmlChar*)sg.str().c_str()); + xmlNewProp( + orb, (const xmlChar*)"rc", (const xmlChar*)rc.str().c_str()); + xmlNewProp( + orb, (const xmlChar*)"a1", (const xmlChar*)a1.str().c_str()); + xmlNewProp( + orb, (const xmlChar*)"a2", (const xmlChar*)a2.str().c_str()); + xmlNewProp( + orb, (const xmlChar*)"a3", (const xmlChar*)a3.str().c_str()); + xmlNewProp( + orb, (const xmlChar*)"a4", (const xmlChar*)a4.str().c_str()); + xmlNewProp( + orb, (const xmlChar*)"a5", (const xmlChar*)a5.str().c_str()); + xmlAddChild(ctr, orb); + } + xmlAddChild(spo, ctr); + } + + app_log() << "Saving resulting cusp Info xml block to: " << filename + << std::endl; + xmlSaveFormatFile(filename.c_str(), doc, 1); + xmlFreeDoc(doc); +} + +template class CuspCorrectionConstructionT; +template class CuspCorrectionConstructionT; + +} // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/LCAO/CuspCorrectionConstructionT.h b/src/QMCWaveFunctions/LCAO/CuspCorrectionConstructionT.h new file mode 100644 index 0000000000..300443c4a0 --- /dev/null +++ b/src/QMCWaveFunctions/LCAO/CuspCorrectionConstructionT.h @@ -0,0 +1,313 @@ +////////////////////////////////////////////////////////////////////////////////////// +// 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. +// +// File developed by: Mark Dewing, mdewing@anl.gov, Argonne National Laboratory +// Ye Luo, yeluo@anl.gov, Argonne National Laboratory +// +// File created by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory +////////////////////////////////////////////////////////////////////////////////////// + +#ifndef QMCPLUSPLUS_CUSP_CORRECTION_CONSTRUCTORT_H +#define QMCPLUSPLUS_CUSP_CORRECTION_CONSTRUCTORT_H + +#include "CuspCorrectionT.h" +#include "LCAOrbitalSetT.h" +#include "QMCWaveFunctions/OrbitalSetTraits.h" +#include "SoaCuspCorrectionT.h" + +class Communicate; +namespace qmcplusplus +{ + +class ParticleSet; + +template +class OneMolecularOrbitalT +{ +public: + using RealType = typename OrbitalSetTraits::RealType; + using ValueType = typename OrbitalSetTraits::ValueType; + using GradType = typename OrbitalSetTraits::GradType; + using ValueVector = typename OrbitalSetTraits::ValueVector; + using GradVector = typename OrbitalSetTraits::GradVector; + using SPOSetPtr = SPOSetT*; + + ValueType + phi(RealType r) + { + TinyVector dr = 0; + dr[0] = r; + + targetPtcl->R[0] = sourcePtcl->R[curCenter]; + targetPtcl->makeMove(0, dr); + Psi1->evaluateValue(*targetPtcl, 0, val1); + + return val1[curOrb]; + } + + void + phi_vgl(RealType r, RealType& val, GradType& grad, RealType& lap) + { + TinyVector dr = 0; + dr[0] = r; + + targetPtcl->R[0] = sourcePtcl->R[curCenter]; + targetPtcl->makeMove(0, dr); + Psi1->evaluateVGL(*targetPtcl, 0, val1, grad1, lap1); + + val = val1[curOrb]; + grad = grad1[curOrb]; + lap = lap1[curOrb]; + } + + OneMolecularOrbitalT( + ParticleSet* targetP, ParticleSet* sourceP, SPOSetPtr Phi) : + targetPtcl(targetP), + sourcePtcl(sourceP), + curOrb(0), + curCenter(0) + { + Psi1 = Phi; + int norb = Psi1->getOrbitalSetSize(); + val1.resize(norb); + grad1.resize(norb); + lap1.resize(norb); + } + + void + changeOrbital(int centerIdx, int orbIdx) + { + curCenter = centerIdx; + curOrb = orbIdx; + } + +private: + /// Temporary storage for real wavefunction values + ValueVector val1; + GradVector grad1; + ValueVector lap1; + + /// target ParticleSet + ParticleSet* targetPtcl; + /// source ParticleSet + ParticleSet* sourcePtcl; + + /// Index of orbital + int curOrb; + + /// Index of atomic center + int curCenter; + + SPOSetPtr Psi1; +}; + +template +class CuspCorrectionConstructionT +{ +public: + using RealType = typename OrbitalSetTraits::RealType; + using ValueType = typename OrbitalSetTraits::ValueType; + using ValueVector = typename OrbitalSetTraits::ValueVector; + using GradType = typename OrbitalSetTraits::GradType; + using GradVector = typename OrbitalSetTraits::GradVector; + + struct ValGradLap + { + ValueType val; + GradType grad; + ValueType lap; + }; + + /// Divide molecular orbital into atomic S-orbitals on this center (phi), + /// and everything else (eta). + static void + splitPhiEta(int center, const std::vector& corrCenter, + LCAOrbitalSetT& phi, LCAOrbitalSetT& eta); + + /// Remove S atomic orbitals from all molecular orbitals on all centers. + static void + removeSTypeOrbitals( + const std::vector& corrCenter, LCAOrbitalSetT& Phi); + + /// Compute the radial part of the corrected wavefunction + static void + computeRadialPhiBar(ParticleSet* targetP, ParticleSet* sourceP, int curOrb_, + int curCenter_, SPOSetT* Phi, Vector& xgrid, + Vector& rad_orb, const CuspCorrectionParametersT& data); + + /** Ideal local energy at one point + * @param r input radial distance + * @param Z nuclear charge + * @param beta0 adjustable parameter to make energy continuous at Rc + */ + static RealType + getOneIdealLocalEnergy(RealType r, RealType Z, RealType beta0); + + /** Ideal local energy at a vector of points + * @param pos input vector of radial distances + * @param Z nuclear charge + * @param Rc cutoff radius where the correction meets the actual orbital + * @param ELorigAtRc local energy at Rc. beta0 is adjusted to make energy + * continuous at Rc + * @param ELideal - output the ideal local energy at pos values + */ + static void + getIdealLocalEnergy(const ValueVector& pos, RealType Z, RealType Rc, + RealType ELorigAtRc, ValueVector& ELideal); + + /** Evaluate various orbital quantities that enter as constraints on the + * correction + * @param valRc orbital value at Rc + * @param gradRc orbital gradient at Rc + * @param lapRc orbital laplacian at Rc + * @param Rc cutoff radius + * @param Z nuclear charge + * @param C offset to keep correction to a single sign + * @param valAtZero orbital value at zero + * @param eta0 value of non-corrected pieces of the orbital at zero + * @param X output + */ + static void + evalX(RealType valRc, GradType gradRc, ValueType lapRc, RealType Rc, + RealType Z, RealType C, RealType valAtZero, RealType eta0, + TinyVector& X); + + /** Convert constraints to polynomial parameters + * @param X input from evalX + * @param Rc cutoff radius + * @param alpha output the polynomial parameters for the correction + */ + static void + X2alpha(const TinyVector& X, RealType Rc, + TinyVector& alpha); + + /** Effective nuclear charge to keep effective local energy finite at zero + * @param Z nuclear charge + * @param etaAtZero value of non-S orbitals at this center + * @param phiBarAtZero value of corrected orbital at zero + */ + static RealType + getZeff(RealType Z, RealType etaAtZero, RealType phiBarAtZero); + + static RealType + phiBar(const CuspCorrectionT& cusp, RealType r, + OneMolecularOrbitalT& phiMO); + + /** Compute effective local energy at vector of points + * @param pos input vector of radial distances + * @param Zeff effective charge from getZeff + * @param Rc cutoff radius + * @param originalELatRc Local energy at the center from the uncorrected + * orbital + * @param cusp cusp correction parameters + * @param phiMO uncorrected orbital (S-orbitals on this center only) + * @param ELcurr output local energy at each distance in pos + */ + static void + getCurrentLocalEnergy(const ValueVector& pos, RealType Zeff, RealType Rc, + RealType originalELatRc, CuspCorrectionT& cusp, + OneMolecularOrbitalT& phiMO, ValueVector& ELcurr); + + /** Local energy from uncorrected orbital + * @param pos input vector of radial distances + * @param Zeff nuclear charge + * @param Rc cutoff radius + * @param phiMO uncorrected orbital (S-orbitals on this center only) + * @param ELorig output local energy at each distance in pos + * + * Return is value of local energy at zero. This is the value needed for + * subsequent computations. The routine can be called with an empty vector + * of positions to get just this value. + */ + static RealType + getOriginalLocalEnergy(const ValueVector& pos, RealType Zeff, RealType Rc, + OneMolecularOrbitalT& phiMO, ValueVector& Elorig); + + /** Sum of squares difference between the current and ideal local energies + * This is the objective function to be minimized. + * @param Elcurr current local energy + * @param Elideal ideal local energy + */ + static RealType + getELchi2(const ValueVector& ELcurr, const ValueVector& ELideal); + + /** Minimize chi2 with respect to phi at zero for a fixed Rc + * @param cusp correction parameters + * @param phiMO uncorrected orbital (S-orbitals on this center only) + * @param Z nuclear charge + * @param eta0 value at zero for parts of the orbital that don't require + * correction - the non-S-orbitals on this center and all orbitals on other + * centers + * @param pos vector of radial positions + * @param Elcurr storage for current local energy + * @param Elideal storage for ideal local energy + */ + static RealType + minimizeForPhiAtZero(CuspCorrectionT& cusp, + OneMolecularOrbitalT& phiMO, RealType Z, RealType eta0, + ValueVector& pos, ValueVector& ELcurr, ValueVector& ELideal, + RealType start_phi0); + + /** Minimize chi2 with respect to Rc and phi at zero. + * @param cusp correction parameters + * @param phiMO uncorrected orbital (S-orbitals on this center only) + * @param Z nuclear charge + * @param Rc_init initial value for Rc + * @param Rc_max maximum value for Rc + * @param eta0 value at zero for parts of the orbital that don't require + * correction - the non-S-orbitals on this center and all orbitals on other + * centers + * @param pos vector of radial positions + * @param Elcurr storage for current local energy + * @param Elideal storage for ideal local energy + * + * Output is parameter values in cusp.cparam + */ + static void + minimizeForRc(CuspCorrectionT& cusp, OneMolecularOrbitalT& phiMO, + RealType Z, RealType Rc_init, RealType Rc_max, RealType eta0, + ValueVector& pos, ValueVector& ELcurr, ValueVector& ELideal); + + // Modifies orbital set lcwc + static void + applyCuspCorrection(const Matrix>& info, + ParticleSet& targetPtcl, ParticleSet& sourcePtcl, + LCAOrbitalSetT& lcao, SoaCuspCorrectionT& cusp, + const std::string& id); + + static void + generateCuspInfo(Matrix>& info, + const ParticleSet& targetPtcl, const ParticleSet& sourcePtcl, + const LCAOrbitalSetT& lcao, const std::string& id, + Communicate& Comm); + + /// Broadcast cusp correction parameters + static void + broadcastCuspInfo( + CuspCorrectionParametersT& param, Communicate& Comm, int root); + + /// Read cusp correction parameters from XML file + static bool + readCuspInfo(const std::string& cuspInfoFile, const std::string& objectName, + int OrbitalSetSize, Matrix>& info); + + /// save cusp correction info to a file. + static void + saveCusp(const std::string& filename, + const Matrix>& info, + const std::string& id); + +private: + static RealType + evaluateForPhi0Body(RealType phi0, ValueVector& pos, ValueVector& ELcurr, + ValueVector& ELideal, CuspCorrectionT& cusp, + OneMolecularOrbitalT& phiMO, ValGradLap phiAtRc, RealType etaAtZero, + RealType ELorigAtRc, RealType Z); +}; + +} // namespace qmcplusplus + +#endif diff --git a/src/QMCWaveFunctions/LCAO/CuspCorrectionT.h b/src/QMCWaveFunctions/LCAO/CuspCorrectionT.h new file mode 100644 index 0000000000..18fa1ed531 --- /dev/null +++ b/src/QMCWaveFunctions/LCAO/CuspCorrectionT.h @@ -0,0 +1,114 @@ +////////////////////////////////////////////////////////////////////////////////////// +// This file is distributed under the University of Illinois/NCSA Open Source +// License. See LICENSE file in top directory for details. +// +// Copyright (c) 2018 Jeongnim Kim and QMCPACK developers. +// +// File developed by: Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore +// National Laboratory +// Jeremy McMinnis, jmcminis@gmail.com, University of +// Illinois at Urbana-Champaign Mark A. Berrill, +// berrillma@ornl.gov, Oak Ridge National Laboratory Mark +// Dewing, mdewing@anl.gov, Argonne National Laboratory +// +// File created by: Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore +// National Laboratory +////////////////////////////////////////////////////////////////////////////////////// + +/** @file CuspCorrection.h + * @brief Corrections to electron-nucleus cusp for all-electron molecular + * calculations. + */ + +#ifndef QMCPLUSPLUS_CUSPCORRECTIONT_H +#define QMCPLUSPLUS_CUSPCORRECTIONT_H + +#include "Configuration.h" +#include "OrbitalSetTraits.h" + +#include + +namespace qmcplusplus +{ +/** + * @brief Cusp correction parameters + * + * From "Scheme for adding electron-nuclear cusps to Gaussian orbitals" Ma, + * Towler, Drummond, and Needs JCP 122, 224322 (2005) + * + * Equations 7 and 8 in the paper define the correction. These are the + * parameters in those equations. + */ + +template +struct CuspCorrectionParametersT +{ + using ValueType = typename OrbitalSetTraits::ValueType; + using RealType = typename OrbitalSetTraits::RealType; + + /// The cutoff radius + RealType Rc; + + /// A shift to keep correction to a single sign + RealType C; + + /// The sign of the wavefunction at the nucleus + RealType sg; + + /// The coefficients of the polynomial \f$p(r)\f$ in Eq 8 + TinyVector alpha; + + /// Flag to indicate the correction should be recalculated + int redo; + + CuspCorrectionParametersT() : Rc(0.0), C(0.0), sg(1.0), alpha(0.0), redo(0) + { + } +}; + +/// Formulas for applying the cusp correction + +template +class CuspCorrectionT +{ + using RealType = typename OrbitalSetTraits::RealType; + +public: + inline RealType + Rr(RealType r) const + { + return cparam.sg * std::exp(pr(r)); + } + + inline RealType + pr(RealType r) const + { + auto& alpha = cparam.alpha; + return alpha[0] + alpha[1] * r + alpha[2] * r * r + + alpha[3] * r * r * r + alpha[4] * r * r * r * r; + } + + inline RealType + dpr(RealType r) const + { + auto& alpha = cparam.alpha; + return alpha[1] + 2.0 * alpha[2] * r + 3.0 * alpha[3] * r * r + + 4.0 * alpha[4] * r * r * r; + } + + inline RealType + d2pr(RealType r) const + { + auto& alpha = cparam.alpha; + return 2.0 * alpha[2] + 6.0 * alpha[3] * r + 12.0 * alpha[4] * r * r; + } + + CuspCorrectionT(const CuspCorrectionParametersT& param) : cparam(param) + { + } + + CuspCorrectionParametersT cparam; +}; +} // namespace qmcplusplus + +#endif diff --git a/src/QMCWaveFunctions/LCAO/LCAOrbitalBuilderT.cpp b/src/QMCWaveFunctions/LCAO/LCAOrbitalBuilderT.cpp new file mode 100644 index 0000000000..4a9de5f847 --- /dev/null +++ b/src/QMCWaveFunctions/LCAO/LCAOrbitalBuilderT.cpp @@ -0,0 +1,1088 @@ +////////////////////////////////////////////////////////////////////////////////////// +// 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. +// +// File developed by: Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore +// National Laboratory +// Jeremy McMinnis, jmcminis@gmail.com, University of +// Illinois at Urbana-Champaign Jaron T. Krogel, +// krogeljt@ornl.gov, Oak Ridge National Laboratory Jeongnim +// Kim, jeongnim.kim@gmail.com, University of Illinois at +// Urbana-Champaign Ye Luo, yeluo@anl.gov, Argonne National +// Laboratory Mark A. Berrill, berrillma@ornl.gov, Oak Ridge +// National Laboratory +// +// File created by: Jeongnim Kim, jeongnim.kim@intel.com, Intel Corp. +////////////////////////////////////////////////////////////////////////////////////// + +#include "LCAOrbitalBuilderT.h" + +#include "AOBasisBuilder.h" +#include "LCAOrbitalSetT.h" +#include "MultiFunctorAdapter.h" +#include "MultiQuinticSpline1D.h" +#include "Numerics/SoaCartesianTensor.h" +#include "Numerics/SoaSphericalTensor.h" +#include "OhmmsData/AttributeSet.h" +#include "QMCWaveFunctions/SPOSetT.h" +#include "SoaAtomicBasisSet.h" +#include "SoaLocalizedBasisSet.h" +#if !defined(QMC_COMPLEX) +#include "CuspCorrectionConstructionT.h" +#include "LCAOrbitalSetWithCorrectionT.h" +#endif +#include "CPU/math.hpp" +#include "Message/CommOperators.h" +#include "Utilities/ProgressReportEngine.h" +#include "hdf/hdf_archive.h" + +#include + +namespace qmcplusplus +{ +/** traits for a localized basis set; used by createBasisSet + * + * T radial function value type + * ORBT orbital value type, can be complex + * ROT {0=numuerica;, 1=gto; 2=sto} + * SH {0=cartesian, 1=spherical} + * If too confusing, inroduce enumeration. + */ +template +struct ao_traits +{ +}; + +/** specialization for numerical-cartesian AO */ +template +struct ao_traits +{ + using radial_type = MultiQuinticSpline1D; + using angular_type = SoaCartesianTensor; + using ao_type = SoaAtomicBasisSet; + using basis_type = SoaLocalizedBasisSet; +}; + +/** specialization for numerical-spherical AO */ +template +struct ao_traits +{ + using radial_type = MultiQuinticSpline1D; + using angular_type = SoaSphericalTensor; + using ao_type = SoaAtomicBasisSet; + using basis_type = SoaLocalizedBasisSet; +}; + +/** specialization for GTO-cartesian AO */ +template +struct ao_traits +{ + using radial_type = MultiFunctorAdapter>; + using angular_type = SoaCartesianTensor; + using ao_type = SoaAtomicBasisSet; + using basis_type = SoaLocalizedBasisSet; +}; + +/** specialization for GTO-cartesian AO */ +template +struct ao_traits +{ + using radial_type = MultiFunctorAdapter>; + using angular_type = SoaSphericalTensor; + using ao_type = SoaAtomicBasisSet; + using basis_type = SoaLocalizedBasisSet; +}; + +/** specialization for STO-spherical AO */ +template +struct ao_traits +{ + using radial_type = MultiFunctorAdapter>; + using angular_type = SoaSphericalTensor; + using ao_type = SoaAtomicBasisSet; + using basis_type = SoaLocalizedBasisSet; +}; + +inline bool +is_same(const xmlChar* a, const char* b) +{ + return !strcmp((const char*)a, b); +} + +template +LCAOrbitalBuilderT::LCAOrbitalBuilderT( + ParticleSet& els, ParticleSet& ions, Communicate* comm, xmlNodePtr cur) : + SPOSetBuilderT("LCAO", comm), + targetPtcl(els), + sourcePtcl(ions), + h5_path(""), + SuperTwist(0.0), + doCuspCorrection(false) +{ + this->ClassName = "LCAOrbitalBuilder"; + ReportEngine PRE(this->ClassName, "createBasisSet"); + + std::string cuspC("no"); // cusp correction + OhmmsAttributeSet aAttrib; + aAttrib.add(cuspC, "cuspCorrection"); + aAttrib.add(h5_path, "href"); + aAttrib.add(PBCImages, "PBCimages"); + aAttrib.add(SuperTwist, "twist"); + aAttrib.put(cur); + + if (cuspC == "yes") + doCuspCorrection = true; + // Evaluate the Phase factor. Equals 1 for OBC. + EvalPeriodicImagePhaseFactors(SuperTwist, PeriodicImagePhaseFactors); + + // no need to wait but load the basis set + processChildren( + cur, [&](const std::string& cname, const xmlNodePtr element) { + if (cname == "basisset") { + std::string basisset_name_input( + getXMLAttributeValue(element, "name")); + std::string basisset_name(basisset_name_input.empty() ? + "LCAOBSet" : + basisset_name_input); + if (basisset_map_.find(basisset_name) != basisset_map_.end()) { + std::ostringstream err_msg; + err_msg << "Cannot create basisset " << basisset_name + << " which already exists." << std::endl; + throw std::runtime_error(err_msg.str()); + } + if (h5_path != "") + basisset_map_[basisset_name] = loadBasisSetFromH5(element); + else + basisset_map_[basisset_name] = + loadBasisSetFromXML(element, cur); + } + }); + + // deprecated h5 basis set handling when basisset element is missing + if (basisset_map_.size() == 0 && h5_path != "") { + app_warning() + << "!!!!!!! Deprecated input style: missing basisset element. " + << "LCAO needs an explicit basisset XML element. " + << "Fallback on loading an implicit one." << std::endl; + basisset_map_["LCAOBSet"] = loadBasisSetFromH5(cur); + } + + if (basisset_map_.size() == 0) + throw std::runtime_error("No basisset found in the XML input!"); +} + +template +LCAOrbitalBuilderT::~LCAOrbitalBuilderT() +{ + // properly cleanup +} + +template +int +LCAOrbitalBuilderT::determineRadialOrbType(xmlNodePtr cur) const +{ + std::string keyOpt; + std::string transformOpt; + OhmmsAttributeSet aAttrib; + aAttrib.add(keyOpt, "keyword"); + aAttrib.add(keyOpt, "key"); + aAttrib.add(transformOpt, "transform"); + aAttrib.put(cur); + + int radialOrbType = -1; + if (transformOpt == "yes" || keyOpt == "NMO") + radialOrbType = 0; + else { + if (keyOpt == "GTO") + radialOrbType = 1; + if (keyOpt == "STO") + radialOrbType = 2; + } + return radialOrbType; +} + +template +std::unique_ptr::BasisSet_t> +LCAOrbitalBuilderT::loadBasisSetFromXML(xmlNodePtr cur, xmlNodePtr parent) +{ + ReportEngine PRE(this->ClassName, "loadBasisSetFromXML(xmlNodePtr)"); + int ylm = -1; + { + xmlNodePtr cur1 = cur->xmlChildrenNode; + while (cur1 != NULL && ylm < 0) { + if (is_same(cur1->name, "atomicBasisSet")) { + std::string sph; + OhmmsAttributeSet att; + att.add(sph, "angular"); + att.put(cur1); + ylm = (sph == "cartesian") ? 0 : 1; + } + cur1 = cur1->next; + } + } + + if (ylm < 0) + PRE.error("Missing angular attribute of atomicBasisSet.", true); + + int radialOrbType = determineRadialOrbType(cur); + if (radialOrbType < 0) { + app_warning() << "Radial orbital type cannot be determined based on " + "the attributes of basisset line. " + << "Trying the parent element." << std::endl; + radialOrbType = determineRadialOrbType(parent); + } + + if (radialOrbType < 0) + PRE.error("Unknown radial function for LCAO orbitals. Specify " + "keyword=\"NMO/GTO/STO\" .", + true); + + BasisSet_t* myBasisSet = nullptr; + /** process atomicBasisSet per ion species */ + switch (radialOrbType) { + case (0): // numerical + app_log() << " LCAO: SoaAtomicBasisSet" + << std::endl; + if (ylm) + myBasisSet = createBasisSet<0, 1>(cur); + else + myBasisSet = createBasisSet<0, 0>(cur); + break; + case (1): // gto + app_log() << " LCAO: SoaAtomicBasisSet" + << std::endl; + if (ylm) + myBasisSet = createBasisSet<1, 1>(cur); + else + myBasisSet = createBasisSet<1, 0>(cur); + break; + case (2): // sto + app_log() << " LCAO: SoaAtomicBasisSet" + << std::endl; + myBasisSet = createBasisSet<2, 1>(cur); + break; + default: + PRE.error("Cannot construct SoaAtomicBasisSet.", true); + break; + } + + return std::unique_ptr(myBasisSet); +} + +template +std::unique_ptr::BasisSet_t> +LCAOrbitalBuilderT::loadBasisSetFromH5(xmlNodePtr parent) +{ + ReportEngine PRE(this->ClassName, "loadBasisSetFromH5()"); + + hdf_archive hin(this->myComm); + int ylm = -1; + if (this->myComm->rank() == 0) { + if (!hin.open(h5_path, H5F_ACC_RDONLY)) + PRE.error("Could not open H5 file", true); + + hin.push("basisset", false); + + std::string sph; + std::string ElemID0 = "atomicBasisSet0"; + + hin.push(ElemID0.c_str(), false); + + if (!hin.readEntry(sph, "angular")) + PRE.error("Could not find name of basisset group in H5; Probably " + "Corrupt H5 file", + true); + ylm = (sph == "cartesian") ? 0 : 1; + hin.close(); + } + + this->myComm->bcast(ylm); + if (ylm < 0) + PRE.error("Missing angular attribute of atomicBasisSet.", true); + + int radialOrbType = determineRadialOrbType(parent); + if (radialOrbType < 0) + PRE.error("Unknown radial function for LCAO orbitals. Specify " + "keyword=\"NMO/GTO/STO\" .", + true); + + BasisSet_t* myBasisSet = nullptr; + /** process atomicBasisSet per ion species */ + switch (radialOrbType) { + case (0): // numerical + app_log() << " LCAO: SoaAtomicBasisSet" + << std::endl; + if (ylm) + myBasisSet = createBasisSetH5<0, 1>(); + else + myBasisSet = createBasisSetH5<0, 0>(); + break; + case (1): // gto + app_log() << " LCAO: SoaAtomicBasisSet" + << std::endl; + if (ylm) + myBasisSet = createBasisSetH5<1, 1>(); + else + myBasisSet = createBasisSetH5<1, 0>(); + break; + case (2): // sto + app_log() << " LCAO: SoaAtomicBasisSet" + << std::endl; + myBasisSet = createBasisSetH5<2, 1>(); + break; + default: + PRE.error("Cannot construct SoaAtomicBasisSet.", true); + break; + } + return std::unique_ptr(myBasisSet); +} + +template +template +typename LCAOrbitalBuilderT::BasisSet_t* +LCAOrbitalBuilderT::createBasisSet(xmlNodePtr cur) +{ + ReportEngine PRE(this->ClassName, "createBasisSet(xmlNodePtr)"); + + using ao_type = typename ao_traits::ao_type; + using basis_type = typename ao_traits::basis_type; + + basis_type* mBasisSet = new basis_type(sourcePtcl, targetPtcl); + + // list of built centers + std::vector ao_built_centers; + + /** process atomicBasisSet per ion species */ + cur = cur->xmlChildrenNode; + while (cur != NULL) // loop over unique ioons + { + std::string cname((const char*)(cur->name)); + + if (cname == "atomicBasisSet") { + std::string elementType; + std::string sph; + OhmmsAttributeSet att; + att.add(elementType, "elementType"); + att.put(cur); + + if (elementType.empty()) + PRE.error( + "Missing elementType attribute of atomicBasisSet.", true); + + auto it = std::find( + ao_built_centers.begin(), ao_built_centers.end(), elementType); + if (it == ao_built_centers.end()) { + AOBasisBuilder any(elementType, this->myComm); + any.put(cur); + auto aoBasis = any.createAOSet(cur); + if (aoBasis) { + // add the new atomic basis to the basis set + int activeCenter = + sourcePtcl.getSpeciesSet().findSpecies(elementType); + mBasisSet->add(activeCenter, std::move(aoBasis)); + } + ao_built_centers.push_back(elementType); + } + } + cur = cur->next; + } // done with basis set + mBasisSet->setBasisSetSize(-1); + mBasisSet->setPBCParams(PBCImages, SuperTwist, PeriodicImagePhaseFactors); + return mBasisSet; +} + +template +template +typename LCAOrbitalBuilderT::BasisSet_t* +LCAOrbitalBuilderT::createBasisSetH5() +{ + ReportEngine PRE(this->ClassName, "createBasisSetH5(xmlNodePtr)"); + + using ao_type = typename ao_traits::ao_type; + using basis_type = typename ao_traits::basis_type; + + basis_type* mBasisSet = new basis_type(sourcePtcl, targetPtcl); + + // list of built centers + std::vector ao_built_centers; + + int Nb_Elements(0); + std::string basiset_name; + + /** process atomicBasisSet per ion species */ + app_log() << "Reading BasisSet from HDF5 file:" << h5_path << std::endl; + + hdf_archive hin(this->myComm); + if (this->myComm->rank() == 0) { + if (!hin.open(h5_path, H5F_ACC_RDONLY)) + PRE.error("Could not open H5 file", true); + + hin.push("basisset", false); + + hin.read(Nb_Elements, "NbElements"); + } + + this->myComm->bcast(Nb_Elements); + if (Nb_Elements < 1) + PRE.error("Missing elementType attribute of atomicBasisSet.", true); + + for (int i = 0; i < Nb_Elements; i++) { + std::string elementType, dataset; + std::stringstream tempElem; + std::string ElemID0 = "atomicBasisSet", ElemType; + tempElem << ElemID0 << i; + ElemType = tempElem.str(); + + if (this->myComm->rank() == 0) { + hin.push(ElemType.c_str(), false); + + if (!hin.readEntry(basiset_name, "name")) + PRE.error("Could not find name of basisset group in H5; " + "Probably Corrupt H5 file", + true); + if (!hin.readEntry(elementType, "elementType")) + PRE.error("Could not read elementType in H5; Probably Corrupt " + "H5 file", + true); + } + this->myComm->bcast(basiset_name); + this->myComm->bcast(elementType); + + auto it = std::find( + ao_built_centers.begin(), ao_built_centers.end(), elementType); + if (it == ao_built_centers.end()) { + AOBasisBuilder any(elementType, this->myComm); + any.putH5(hin); + auto aoBasis = any.createAOSetH5(hin); + if (aoBasis) { + // add the new atomic basis to the basis set + int activeCenter = + sourcePtcl.getSpeciesSet().findSpecies(elementType); + mBasisSet->add(activeCenter, std::move(aoBasis)); + } + ao_built_centers.push_back(elementType); + } + + if (this->myComm->rank() == 0) + hin.pop(); + } + + if (this->myComm->rank() == 0) { + hin.pop(); + hin.close(); + } + mBasisSet->setBasisSetSize(-1); + mBasisSet->setPBCParams(PBCImages, SuperTwist, PeriodicImagePhaseFactors); + return mBasisSet; +} + +template +std::unique_ptr> +LCAOrbitalBuilderT::createSPOSetFromXML(xmlNodePtr cur) +{ + ReportEngine PRE(this->ClassName, "createSPO(xmlNodePtr)"); + std::string spo_name(""), cusp_file(""), optimize("no"); + std::string basisset_name("LCAOBSet"); + OhmmsAttributeSet spoAttrib; + spoAttrib.add(spo_name, "name"); + spoAttrib.add(spo_name, "id"); + spoAttrib.add(cusp_file, "cuspInfo"); + spoAttrib.add(basisset_name, "basisset"); + spoAttrib.put(cur); + + std::unique_ptr myBasisSet; + if (basisset_map_.find(basisset_name) == basisset_map_.end()) + this->myComm->barrier_and_abort( + "basisset \"" + basisset_name + "\" cannot be found\n"); + else + myBasisSet.reset(basisset_map_[basisset_name]->makeClone()); + + std::unique_ptr> sposet; + if (doCuspCorrection) { +#if defined(QMC_COMPLEX) + this->myComm->barrier_and_abort( + "LCAOrbitalBuilder::createSPOSetFromXML cusp correction is not " + "supported on complex LCAO."); +#else + app_summary() << " Using cusp correction." << std::endl; + auto lcwc = std::make_unique>( + spo_name, sourcePtcl, targetPtcl, std::move(myBasisSet)); + loadMO(lcwc->lcao, cur); + lcwc->setOrbitalSetSize(lcwc->lcao.getOrbitalSetSize()); + sposet = std::move(lcwc); +#endif + } + else { + auto lcos = std::make_unique>( + spo_name, std::move(myBasisSet)); + loadMO(*lcos, cur); + sposet = std::move(lcos); + } + +#if !defined(QMC_COMPLEX) + if (doCuspCorrection) { + // Create a temporary particle set to use for cusp initialization. + // The particle coordinates left at the end are unsuitable for further + // computations. The coordinates get set to nuclear positions, which + // leads to zero e-N distance, which causes a NaN in SoaAtomicBasisSet.h + // This problem only appears when the electron positions are specified + // in the input. The random particle placement step executes after this + // part of the code, overwriting the leftover positions from the cusp + // initialization. + ParticleSet tmp_targetPtcl(targetPtcl); + + const int num_centers = sourcePtcl.getTotalNum(); + auto& lcwc = dynamic_cast&>(*sposet); + + const int orbital_set_size = lcwc.getOrbitalSetSize(); + Matrix> info( + num_centers, orbital_set_size); + + // set a default file name if not given + if (cusp_file.empty()) + cusp_file = spo_name + ".cuspInfo.xml"; + + bool file_exists( + this->myComm->rank() == 0 && std::ifstream(cusp_file).good()); + this->myComm->bcast(file_exists); + app_log() << " Cusp correction file " << cusp_file + << (file_exists ? " exits." : " doesn't exist.") << std::endl; + + // validate file if it exists + if (file_exists) { + bool valid = 0; + if (this->myComm->rank() == 0) + valid = CuspCorrectionConstructionT::readCuspInfo( + cusp_file, spo_name, orbital_set_size, info); + this->myComm->bcast(valid); + if (!valid) + this->myComm->barrier_and_abort( + "Invalid cusp correction file " + cusp_file); +#ifdef HAVE_MPI + for (int orb_idx = 0; orb_idx < orbital_set_size; orb_idx++) + for (int center_idx = 0; center_idx < num_centers; center_idx++) + broadcastCuspInfo( + info(center_idx, orb_idx), *this->myComm, 0); +#endif + } + else { + CuspCorrectionConstructionT::generateCuspInfo(info, + tmp_targetPtcl, sourcePtcl, lcwc.lcao, spo_name, *this->myComm); + if (this->myComm->rank() == 0) + CuspCorrectionConstructionT::saveCusp( + cusp_file, info, spo_name); + } + + CuspCorrectionConstructionT::applyCuspCorrection( + info, tmp_targetPtcl, sourcePtcl, lcwc.lcao, lcwc.cusp, spo_name); + } +#endif + + return sposet; +} + +/** Parse the xml file for information on the Dirac determinants. + *@param cur the current xmlNode + */ +template +bool +LCAOrbitalBuilderT::loadMO(LCAOrbitalSetT& spo, xmlNodePtr cur) +{ +#undef FunctionName +#define FunctionName \ + printf("Calling FunctionName from %s\n", __FUNCTION__); \ + FunctionNameReal + // Check if HDF5 present + ReportEngine PRE("LCAOrbitalBuilder", "put(xmlNodePtr)"); + + // initialize the number of orbital by the basis set size + int norb = spo.getBasisSetSize(); + std::string debugc("no"); + double orbital_mix_magnitude = 0.0; + bool PBC = false; + OhmmsAttributeSet aAttrib; + aAttrib.add(norb, "orbitals"); + aAttrib.add(norb, "size"); + aAttrib.add(debugc, "debug"); + aAttrib.add(orbital_mix_magnitude, "orbital_mix_magnitude"); + aAttrib.put(cur); + xmlNodePtr occ_ptr = NULL; + xmlNodePtr coeff_ptr = NULL; + cur = cur->xmlChildrenNode; + while (cur != NULL) { + std::string cname((const char*)(cur->name)); + if (cname == "occupation") { + occ_ptr = cur; + } + else if (cname.find("coeff") < cname.size() || cname == "parameter" || + cname == "Var") { + coeff_ptr = cur; + } + cur = cur->next; + } + if (coeff_ptr == NULL) { + app_log() << " Using Identity for the LCOrbitalSet " << std::endl; + return true; + } + spo.setOrbitalSetSize(norb); + bool success = putOccupation(spo, occ_ptr); + if (h5_path == "") + success = putFromXML(spo, coeff_ptr); + else { + hdf_archive hin(this->myComm); + + if (this->myComm->rank() == 0) { + if (!hin.open(h5_path, H5F_ACC_RDONLY)) + APP_ABORT("LCAOrbitalBuilder::putFromH5 missing or incorrect " + "path to H5 file."); + + try { + hin.push("PBC", false); + PBC = true; + } + catch (const std::exception& e) { + app_debug() << e.what() << std::endl; + PBC = false; + } + + if (PBC) + hin.read(PBC, "PBC"); + + hin.close(); + } + this->myComm->bcast(PBC); + if (PBC) + success = putPBCFromH5(spo, coeff_ptr); + else + success = putFromH5(spo, coeff_ptr); + } + + // Ye: used to construct cusp correction + // bool success2 = transformSPOSet(); + if (debugc == "yes") { + app_log() << " Single-particle orbital coefficients dims=" + << spo.C->rows() << " x " << spo.C->cols() << std::endl; + app_log() << *spo.C << std::endl; + } + + return success; +} + +template +bool +LCAOrbitalBuilderT::putFromXML(LCAOrbitalSetT& spo, xmlNodePtr coeff_ptr) +{ + int norbs = 0; + OhmmsAttributeSet aAttrib; + aAttrib.add(norbs, "size"); + aAttrib.add(norbs, "orbitals"); + aAttrib.put(coeff_ptr); + if (norbs < spo.getOrbitalSetSize()) { + return false; + APP_ABORT("LCAOrbitalBuilder::putFromXML missing or incorrect size"); + } + if (norbs) { + std::vector Ctemp; + int BasisSetSize = spo.getBasisSetSize(); + Ctemp.resize(norbs * BasisSetSize); + putContent(Ctemp, coeff_ptr); + int n = 0, i = 0; + typename std::vector::iterator cit(Ctemp.begin()); + while (i < spo.getOrbitalSetSize()) { + if (Occ[n] > std::numeric_limits::epsilon()) { + std::copy(cit, cit + BasisSetSize, (*spo.C)[i]); + i++; + } + n++; + cit += BasisSetSize; + } + } + return true; +} + +/** read data from a hdf5 file + * @param norb number of orbitals to be initialized + * @param coeff_ptr xmlnode for coefficients + */ +template +bool +LCAOrbitalBuilderT::putFromH5(LCAOrbitalSetT& spo, xmlNodePtr coeff_ptr) +{ + int neigs = spo.getBasisSetSize(); + int setVal = -1; + OhmmsAttributeSet aAttrib; + aAttrib.add(setVal, "spindataset"); + aAttrib.add(neigs, "size"); + aAttrib.add(neigs, "orbitals"); + aAttrib.put(coeff_ptr); + hdf_archive hin(this->myComm); + if (this->myComm->rank() == 0) { + if (!hin.open(h5_path, H5F_ACC_RDONLY)) + APP_ABORT("LCAOrbitalBuilder::putFromH5 missing or incorrect path " + "to H5 file."); + + Matrix Ctemp; + std::array name; + + // This is to make sure of Backward compatibility with previous tags. + int name_len = std::snprintf( + name.data(), name.size(), "%s%d", "/Super_Twist/eigenset_", setVal); + if (name_len < 0) + throw std::runtime_error("Error generating name"); + std::string setname(name.data(), name_len); + if (!hin.readEntry(Ctemp, setname)) { + name_len = std::snprintf( + name.data(), name.size(), "%s%d", "/KPTS_0/eigenset_", setVal); + if (name_len < 0) + throw std::runtime_error("Error generating name"); + setname = std::string(name.data(), name_len); + hin.read(Ctemp, setname); + } + hin.close(); + + if (Ctemp.cols() != spo.getBasisSetSize()) { + std::ostringstream err_msg; + err_msg << "Basis set size " << spo.getBasisSetSize() + << " mismatched the number of MO coefficients columns " + << Ctemp.cols() << " from h5." << std::endl; + this->myComm->barrier_and_abort(err_msg.str()); + } + + int norbs = spo.getOrbitalSetSize(); + if (Ctemp.rows() < norbs) { + std::ostringstream err_msg; + err_msg << "Need " << norbs + << " orbitals. Insufficient rows of MO coefficients " + << Ctemp.rows() << " from h5." << std::endl; + this->myComm->barrier_and_abort(err_msg.str()); + } + + int n = 0, i = 0; + while (i < norbs) { + if (Occ[n] > 0.0) { + std::copy(Ctemp[n], Ctemp[n + 1], (*spo.C)[i]); + i++; + } + n++; + } + } + this->myComm->bcast(spo.C->data(), spo.C->size()); + return true; +} + +/** read data from a hdf5 file + * @param norb number of orbitals to be initialized + * @param coeff_ptr xmlnode for coefficients + */ +template +bool +LCAOrbitalBuilderT::putPBCFromH5( + LCAOrbitalSetT& spo, xmlNodePtr coeff_ptr) +{ + ReportEngine PRE("LCAOrbitalBuilder", "LCAOrbitalBuilder::putPBCFromH5"); + int norbs = spo.getOrbitalSetSize(); + int neigs = spo.getBasisSetSize(); + int setVal = -1; + bool IsComplex = false; + bool MultiDet = false; + PosType SuperTwist(0.0); + PosType SuperTwistH5(0.0); + OhmmsAttributeSet aAttrib; + aAttrib.add(setVal, "spindataset"); + aAttrib.add(neigs, "size"); + aAttrib.add(neigs, "orbitals"); + aAttrib.put(coeff_ptr); + hdf_archive hin(this->myComm); + + xmlNodePtr curtemp = coeff_ptr; + + std::string xmlTag("determinantset"); + std::string MSDTag("sposet"); + std::string SDTag("determinant"); + std::string EndTag("qmcsystem"); + std::string curname; + + do { + std::stringstream ss; + curtemp = curtemp->parent; + ss << curtemp->name; + ss >> curname; + if (curname == MSDTag) + MultiDet = true; /// Used to know if running an MSD calculation - + /// needed for order of Orbitals. + if (curname == SDTag) + MultiDet = false; + + } while ((xmlTag != curname) && (curname != EndTag)); + if (curname == EndTag) { + APP_ABORT("Could not find in wf file the \"sposet\" or \"determinant\" " + "tags. Please verify input or contact developers"); + } + + aAttrib.add(SuperTwist, "twist"); + aAttrib.put(curtemp); + + if (this->myComm->rank() == 0) { + if (!hin.open(h5_path, H5F_ACC_RDONLY)) + APP_ABORT("LCAOrbitalBuilder::putFromH5 missing or incorrect path " + "to H5 file."); + hin.push("parameters"); + hin.read(IsComplex, "IsComplex"); + hin.pop(); + + std::string setname("/Super_Twist/Coord"); + hin.read(SuperTwistH5, setname); + if (std::abs(SuperTwistH5[0] - SuperTwist[0]) >= 1e-6 || + std::abs(SuperTwistH5[1] - SuperTwist[1]) >= 1e-6 || + std::abs(SuperTwistH5[2] - SuperTwist[2]) >= 1e-6) { + app_log() << "Super Twist in XML : " << SuperTwist[0] + << " In H5:" << SuperTwistH5[0] << std::endl; + app_log() << " " << SuperTwist[1] + << " " << SuperTwistH5[1] << std::endl; + app_log() << " " << SuperTwist[2] + << " " << SuperTwistH5[2] << std::endl; + app_log() << "Diff in Coord x :" + << std::abs(SuperTwistH5[0] - SuperTwist[0]) << std::endl; + app_log() << " y :" + << std::abs(SuperTwistH5[1] - SuperTwist[1]) << std::endl; + app_log() << " z :" + << std::abs(SuperTwistH5[2] - SuperTwist[2]) << std::endl; + APP_ABORT("Requested Super Twist in XML and Super Twist in HDF5 do " + "not Match!!! Aborting."); + } + // SuperTwist=SuperTwistH5; + Matrix Ctemp; + LoadFullCoefsFromH5(hin, setVal, SuperTwist, Ctemp, MultiDet); + + int n = 0, i = 0; + while (i < norbs) { + if (Occ[n] > 0.0) { + std::copy(Ctemp[n], Ctemp[n + 1], (*spo.C)[i]); + i++; + } + n++; + } + + hin.close(); + } +#ifdef HAVE_MPI + this->myComm->comm.broadcast_n(spo.C->data(), spo.C->size()); +#endif + return true; +} + +template +bool +LCAOrbitalBuilderT::putOccupation(LCAOrbitalSetT& spo, xmlNodePtr occ_ptr) +{ + // die?? + if (spo.getBasisSetSize() == 0) { + APP_ABORT( + "LCAOrbitalBuilder::putOccupation detected ZERO BasisSetSize"); + return false; + } + Occ.resize(std::max(spo.getBasisSetSize(), spo.getOrbitalSetSize())); + Occ = 0.0; + for (int i = 0; i < spo.getOrbitalSetSize(); i++) + Occ[i] = 1.0; + std::vector occ_in; + std::string occ_mode("table"); + if (occ_ptr == NULL) { + occ_mode = "ground"; + } + else { + const std::string o(getXMLAttributeValue(occ_ptr, "mode")); + if (!o.empty()) + occ_mode = o; + } + // Do nothing if mode == ground + if (occ_mode == "excited") { + putContent(occ_in, occ_ptr); + for (int k = 0; k < occ_in.size(); k++) { + if (occ_in[k] < 0) // remove this, -1 is to adjust the base + Occ[-occ_in[k] - 1] = 0.0; + else + Occ[occ_in[k] - 1] = 1.0; + } + } + else if (occ_mode == "table") { + putContent(Occ, occ_ptr); + } + return true; +} + +template +void +LCAOrbitalBuilderT::readRealMatrixFromH5( + hdf_archive& hin, const std::string& setname, Matrix& Creal) const +{ + hin.read(Creal, setname); +} + +template +void +LCAOrbitalBuilderT::LoadFullCoefsFromH5(hdf_archive& hin, int setVal, + PosType& SuperTwist, Matrix>& Ctemp, bool MultiDet) +{ + Matrix Creal; + Matrix Ccmplx; + + std::array name; + int name_len{0}; + /// When running Single Determinant calculations, MO coeff loaded based on + /// occupation and lowest eingenvalue. However, for solids with + /// multideterminants, orbitals are order by kpoints; first all MOs for + /// kpoint 1, then 2 etc + /// The multideterminants occupation is specified in the input/HDF5 and + /// theefore as long as there is consistency between the order in which we + /// read the orbitals and the occupation, we are safe. In the case of + /// Multideterminants generated by pyscf and Quantum Package, They are + /// stored in the same order as generated for quantum package and one + /// should use the orbitals labelled eigenset_unsorted. + + if (MultiDet == false) + name_len = std::snprintf( + name.data(), name.size(), "%s%d", "/Super_Twist/eigenset_", setVal); + else + name_len = std::snprintf(name.data(), name.size(), "%s%d", + "/Super_Twist/eigenset_unsorted_", setVal); + if (name_len < 0) + throw std::runtime_error("Error generating name"); + + std::string setname(name.data(), name_len); + readRealMatrixFromH5(hin, setname, Creal); + + bool IsComplex = true; + hin.read(IsComplex, "/parameters/IsComplex"); + if (IsComplex == false) { + Ccmplx.resize(Creal.rows(), Creal.cols()); + Ccmplx = 0.0; + } + else { + setname += "_imag"; + readRealMatrixFromH5(hin, setname, Ccmplx); + } + + Ctemp.resize(Creal.rows(), Creal.cols()); + for (int i = 0; i < Ctemp.rows(); i++) + for (int j = 0; j < Ctemp.cols(); j++) + Ctemp[i][j] = std::complex(Creal[i][j], Ccmplx[i][j]); +} + +template +void +LCAOrbitalBuilderT::LoadFullCoefsFromH5(hdf_archive& hin, int setVal, + PosType& SuperTwist, Matrix& Creal, bool MultiDet) +{ + bool IsComplex = false; + hin.read(IsComplex, "/parameters/IsComplex"); + if (IsComplex && + (std::abs(SuperTwist[0]) >= 1e-6 || std::abs(SuperTwist[1]) >= 1e-6 || + std::abs(SuperTwist[2]) >= 1e-6)) { + std::string setname( + "This Wavefunction is Complex and you are using the real version " + "of QMCPACK. " + "Please re-run this job with the Complex build of QMCPACK."); + APP_ABORT(setname.c_str()); + } + + std::array name; + int name_len{0}; + bool PBC = false; + hin.read(PBC, "/PBC/PBC"); + if (MultiDet && PBC) + name_len = std::snprintf(name.data(), name.size(), "%s%d", + "/Super_Twist/eigenset_unsorted_", setVal); + else + name_len = std::snprintf( + name.data(), name.size(), "%s%d", "/Super_Twist/eigenset_", setVal); + if (name_len < 0) + throw std::runtime_error("Error generating name"); + + readRealMatrixFromH5(hin, std::string(name.data(), name_len), Creal); +} + +/// Periodic Image Phase Factors computation to be determined +template +void +LCAOrbitalBuilderT::EvalPeriodicImagePhaseFactors( + PosType SuperTwist, std::vector& LocPeriodicImagePhaseFactors) +{ + const int NbImages = + (PBCImages[0] + 1) * (PBCImages[1] + 1) * (PBCImages[2] + 1); + LocPeriodicImagePhaseFactors.resize(NbImages); + for (size_t i = 0; i < NbImages; i++) + LocPeriodicImagePhaseFactors[i] = 1.0; +} + +template +void +LCAOrbitalBuilderT::EvalPeriodicImagePhaseFactors(PosType SuperTwist, + std::vector>& LocPeriodicImagePhaseFactors) +{ + // Allow computation to continue with no HDF file if the system has open + // boundary conditions. The complex build is usually only used with open BC + // for testing. + bool usesOpenBC = + PBCImages[0] == 0 && PBCImages[1] == 0 && PBCImages[2] == 0; + + /// Exp(ik.g) where i is imaginary, k is the supertwist and g is the + /// translation vector PBCImage. + if (h5_path != "" && !usesOpenBC) { + hdf_archive hin(this->myComm); + if (this->myComm->rank() == 0) { + if (!hin.open(h5_path, H5F_ACC_RDONLY)) + APP_ABORT("Could not open H5 file"); + + hin.push("Cell", false); + + hin.read(Lattice, "LatticeVectors"); + hin.close(); + } + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + this->myComm->bcast(Lattice(i, j)); + } + else if (!usesOpenBC) { + APP_ABORT("Attempting to run PBC LCAO with no HDF5 support. Behaviour " + "is unknown. Safer to exit"); + } + + int phase_idx = 0; + int TransX, TransY, TransZ; + RealType phase; + + for (int i = 0; i <= PBCImages[0]; i++) // loop Translation over X + { + TransX = ((i % 2) * 2 - 1) * ((i + 1) / 2); + for (int j = 0; j <= PBCImages[1]; j++) // loop Translation over Y + { + TransY = ((j % 2) * 2 - 1) * ((j + 1) / 2); + for (int k = 0; k <= PBCImages[2]; k++) // loop Translation over Z + { + TransZ = ((k % 2) * 2 - 1) * ((k + 1) / 2); + RealType s, c; + PosType Val; + Val[0] = TransX * Lattice(0, 0) + TransY * Lattice(1, 0) + + TransZ * Lattice(2, 0); + Val[1] = TransX * Lattice(0, 1) + TransY * Lattice(1, 1) + + TransZ * Lattice(2, 1); + Val[2] = TransX * Lattice(0, 2) + TransY * Lattice(1, 2) + + TransZ * Lattice(2, 2); + + phase = dot(SuperTwist, Val); + qmcplusplus::sincos(phase, &s, &c); + + LocPeriodicImagePhaseFactors.emplace_back(c, s); + } + } + } +} + +template class LCAOrbitalBuilderT; +template class LCAOrbitalBuilderT; +template class LCAOrbitalBuilderT>; +template class LCAOrbitalBuilderT>; +} // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/LCAO/LCAOrbitalBuilderT.h b/src/QMCWaveFunctions/LCAO/LCAOrbitalBuilderT.h new file mode 100644 index 0000000000..a746326df7 --- /dev/null +++ b/src/QMCWaveFunctions/LCAO/LCAOrbitalBuilderT.h @@ -0,0 +1,141 @@ +////////////////////////////////////////////////////////////////////////////////////// +// 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. +// +// File developed by: Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore +// National Laboratory +// Jeremy McMinnis, jmcminis@gmail.com, University of +// Illinois at Urbana-Champaign Jaron T. Krogel, +// krogeljt@ornl.gov, Oak Ridge National Laboratory Jeongnim +// Kim, jeongnim.kim@gmail.com, University of Illinois at +// Urbana-Champaign Ye Luo, yeluo@anl.gov, Argonne National +// Laboratory Mark A. Berrill, berrillma@ornl.gov, Oak Ridge +// National Laboratory +// +// File created by: Jeongnim Kim, jeongnim.kim@intel.com, Intel Corp. +////////////////////////////////////////////////////////////////////////////////////// + +#ifndef QMCPLUSPLUS_SOA_LCAO_ORBITAL_BUILDERT_H +#define QMCPLUSPLUS_SOA_LCAO_ORBITAL_BUILDERT_H + +#include "QMCWaveFunctions/BasisSetBase.h" +#include "QMCWaveFunctions/LCAO/LCAOrbitalSetT.h" +#include "QMCWaveFunctions/SPOSetBuilderT.h" + +#include + +namespace qmcplusplus +{ +/** SPOSetBuilder using new LCAOrbitalSet and Soa versions + * + * Reimplement MolecularSPOSetBuilder + * - support both CartesianTensor and SphericalTensor + */ +template +class LCAOrbitalBuilderT : public SPOSetBuilderT +{ +public: + using BasisSet_t = typename LCAOrbitalSetT::basis_type; + using RealType = typename LCAOrbitalSetT::RealType; + using PosType = typename LCAOrbitalSetT::PosType; + + /** constructor + * \param els reference to the electrons + * \param ions reference to the ions + */ + LCAOrbitalBuilderT( + ParticleSet& els, ParticleSet& ions, Communicate* comm, xmlNodePtr cur); + ~LCAOrbitalBuilderT() override; + std::unique_ptr> + createSPOSetFromXML(xmlNodePtr cur) override; + +protected: + /// target ParticleSet + ParticleSet& targetPtcl; + /// source ParticleSet + ParticleSet& sourcePtcl; + /// localized basis set map + std::map> basisset_map_; + /// if true, add cusp correction to orbitals + bool cuspCorr; + /// Path to HDF5 Wavefunction + std::string h5_path; + /// Number of periodic Images for Orbital evaluation + TinyVector PBCImages; + /// Coordinates Super Twist + PosType SuperTwist; + /// Periodic Image Phase Factors. Correspond to the phase from the + /// PBCImages. Computed only once. + std::vector PeriodicImagePhaseFactors; + /// Store Lattice parameters from HDF5 to use in PeriodicImagePhaseFactors + Tensor Lattice; + + /// Enable cusp correction + bool doCuspCorrection; + + /** create basis set + * + * Use ao_traits to match (ROT)x(SH) combo + */ + template + BasisSet_t* + createBasisSet(xmlNodePtr cur); + template + BasisSet_t* + createBasisSetH5(); + + // The following items were previously in SPOSet + /// occupation number + Vector Occ; + bool + loadMO(LCAOrbitalSetT& spo, xmlNodePtr cur); + bool + putOccupation(LCAOrbitalSetT& spo, xmlNodePtr occ_ptr); + bool + putFromXML(LCAOrbitalSetT& spo, xmlNodePtr coeff_ptr); + bool + putFromH5(LCAOrbitalSetT& spo, xmlNodePtr coeff_ptr); + bool + putPBCFromH5(LCAOrbitalSetT& spo, xmlNodePtr coeff_ptr); + // the dimensions of Ctemp are determined by the dataset on file + void + LoadFullCoefsFromH5(hdf_archive& hin, int setVal, PosType& SuperTwist, + Matrix>& Ctemp, bool MultiDet); + // the dimensions of Creal are determined by the dataset on file + void + LoadFullCoefsFromH5(hdf_archive& hin, int setVal, PosType& SuperTwist, + Matrix& Creal, bool Multidet); + void + EvalPeriodicImagePhaseFactors(PosType SuperTwist, + std::vector& LocPeriodicImagePhaseFactors); + void + EvalPeriodicImagePhaseFactors(PosType SuperTwist, + std::vector>& LocPeriodicImagePhaseFactors); + /** read matrix from h5 file + * \param[in] hin: hdf5 arhive to be read from + * \param setname: where to read from in hdf5 archive + * \param[out] Creal: matrix read from h5 + * + * added in header to allow use from derived class LCAOSpinorBuilder as well + */ + void + readRealMatrixFromH5(hdf_archive& hin, const std::string& setname, + Matrix& Creal) const; + +private: + /// load a basis set from XML input + std::unique_ptr + loadBasisSetFromXML(xmlNodePtr cur, xmlNodePtr parent); + /// load a basis set from h5 file + std::unique_ptr + loadBasisSetFromH5(xmlNodePtr parent); + /// determine radial orbital type based on "keyword" and "transform" + /// attributes + int + determineRadialOrbType(xmlNodePtr cur) const; +}; + +} // namespace qmcplusplus +#endif diff --git a/src/QMCWaveFunctions/LCAO/LCAOrbitalSetT.h b/src/QMCWaveFunctions/LCAO/LCAOrbitalSetT.h index 6df0013bd5..974add33b6 100644 --- a/src/QMCWaveFunctions/LCAO/LCAOrbitalSetT.h +++ b/src/QMCWaveFunctions/LCAO/LCAOrbitalSetT.h @@ -1,6 +1,6 @@ ////////////////////////////////////////////////////////////////////////////////////// -// This file is distributed under the University of Illinois/NCSA Open Source License. -// See LICENSE file in top directory for details. +// This file is distributed under the University of Illinois/NCSA Open Source +// License. See LICENSE file in top directory for details. // // Copyright (c) 2016 Jeongnim Kim and QMCPACK developers. // @@ -9,328 +9,374 @@ // File created by: Jeongnim Kim, jeongnim.kim@intel.com, Intel Corp. ////////////////////////////////////////////////////////////////////////////////////// - #ifndef QMCPLUSPLUS_SOA_LINEARCOMIBINATIONORBITALSETT_H #define QMCPLUSPLUS_SOA_LINEARCOMIBINATIONORBITALSETT_H -#include -#include "QMCWaveFunctions/SPOSetT.h" +#include "Numerics/DeterminantOperators.h" +#include "Numerics/MatrixOperators.h" #include "QMCWaveFunctions/BasisSetBase.h" +#include "QMCWaveFunctions/SPOSetT.h" -#include "Numerics/MatrixOperators.h" -#include "Numerics/DeterminantOperators.h" +#include namespace qmcplusplus { -/** class to handle linear combinations of basis orbitals used to evaluate the Dirac determinants. - * - * SoA verson of LCOrtbitalSet - * Localized basis set is always real - */ -template +/** class to handle linear combinations of basis orbitals used to evaluate the + * Dirac determinants. + * + * SoA verson of LCOrtbitalSet + * Localized basis set is always real + */ +template class LCAOrbitalSetT : public SPOSetT { public: - using basis_type = SoaBasisSetBase; - using vgl_type = typename basis_type::vgl_type; - using vgh_type = typename basis_type::vgh_type; - using vghgh_type = typename basis_type::vghgh_type; - - using IndexType = typename SPOSetT::IndexType; - using RealType = typename SPOSetT::RealType; - using ComplexType = typename SPOSetT::ComplexType; - using ValueVector = typename SPOSetT::ValueVector; - using ValueMatrix = typename SPOSetT::ValueMatrix; - using GradVector = typename SPOSetT::GradVector; - using GradMatrix = typename SPOSetT::GradMatrix; - using HessMatrix = typename SPOSetT::HessMatrix; - using PosType = typename SPOSetT::PosType; - using HessVector = typename SPOSetT::HessVector; - using GGGMatrix = typename SPOSetT::GGGMatrix; - using GGGVector = typename SPOSetT::GGGVector; - using GradType = typename SPOSetT::GradType; - using OffloadMWVGLArray = Array>; // [VGL, walker, Orbs] - using OffloadMWVArray = Array>; // [walker, Orbs] - - ///pointer to the basis set - std::unique_ptr myBasisSet; - /// pointer to matrix containing the coefficients - std::shared_ptr C; - - /** constructor + using basis_type = SoaBasisSetBase; + using vgl_type = typename basis_type::vgl_type; + using vgh_type = typename basis_type::vgh_type; + using vghgh_type = typename basis_type::vghgh_type; + + using IndexType = typename SPOSetT::IndexType; + using RealType = typename SPOSetT::RealType; + using ComplexType = typename SPOSetT::ComplexType; + using ValueVector = typename SPOSetT::ValueVector; + using ValueMatrix = typename SPOSetT::ValueMatrix; + using GradVector = typename SPOSetT::GradVector; + using GradMatrix = typename SPOSetT::GradMatrix; + using HessMatrix = typename SPOSetT::HessMatrix; + using PosType = typename SPOSetT::PosType; + using HessVector = typename SPOSetT::HessVector; + using GGGMatrix = typename SPOSetT::GGGMatrix; + using GGGVector = typename SPOSetT::GGGVector; + using GradType = typename SPOSetT::GradType; + using OffloadMWVGLArray = typename basis_type::OffloadMWVGLArray; + using OffloadMWVArray = typename basis_type::OffloadMWVArray; + + /// pointer to the basis set + std::unique_ptr myBasisSet; + /// pointer to matrix containing the coefficients + std::shared_ptr C; + + /** constructor * @param bs pointer to the BasisSet */ - LCAOrbitalSetT(const std::string& my_name, std::unique_ptr&& bs); - - LCAOrbitalSetT(const LCAOrbitalSetT& in); - - std::string getClassName() const final { return "LCAOrbitalSetT"; } - - bool isRotationSupported() const final { return true; } - - bool hasIonDerivs() const final { return true; } - - std::unique_ptr> makeClone() const final; - - void storeParamsBeforeRotation() final { C_copy = std::make_shared(*C); } - - void applyRotation(const ValueMatrix& rot_mat, bool use_stored_copy) final; - - /** set the OrbitalSetSize and Identity=false and initialize internal storages - */ - void setOrbitalSetSize(int norbs) final; - - /** return the size of the basis set - */ - int getBasisSetSize() const { return (myBasisSet == nullptr) ? 0 : myBasisSet->getBasisSetSize(); } - - bool isIdentity() const { return Identity; }; - - /** check consistency between Identity and C - * - */ - void checkObject() const final; - - void evaluateValue(const ParticleSet& P, int iat, ValueVector& psi) final; - - void evaluateVGL(const ParticleSet& P, int iat, ValueVector& psi, GradVector& dpsi, ValueVector& d2psi) final; - - void mw_evaluateValue(const RefVectorWithLeader>& spo_list, - const RefVectorWithLeader& P_list, - int iat, - const RefVector& psi_v_list) const final; - - void mw_evaluateVGL(const RefVectorWithLeader>& spo_list, - const RefVectorWithLeader& P_list, - int iat, - const RefVector& psi_v_list, - const RefVector& dpsi_v_list, - const RefVector& d2psi_v_list) const final; - - void mw_evaluateDetRatios(const RefVectorWithLeader>& spo_list, - const RefVectorWithLeader& vp_list, - const RefVector& psi_list, - const std::vector& invRow_ptr_list, - std::vector>& ratios_list) const final; - - void evaluateDetRatios(const VirtualParticleSet& VP, - ValueVector& psi, - const ValueVector& psiinv, - std::vector& ratios) final; - - void mw_evaluateVGLandDetRatioGrads(const RefVectorWithLeader>& spo_list, - const RefVectorWithLeader& P_list, - int iat, - const std::vector& invRow_ptr_list, - OffloadMWVGLArray& phi_vgl_v, - std::vector& ratios, - std::vector& grads) const final; - - void evaluateVGH(const ParticleSet& P, - int iat, - ValueVector& psi, - GradVector& dpsi, - HessVector& grad_grad_psi) final; - - void evaluateVGHGH(const ParticleSet& P, - int iat, - ValueVector& psi, - GradVector& dpsi, - HessVector& grad_grad_psi, - GGGVector& grad_grad_grad_psi) final; - - void evaluate_notranspose(const ParticleSet& P, - int first, - int last, - ValueMatrix& logdet, - GradMatrix& dlogdet, - ValueMatrix& d2logdet) final; - - void evaluate_notranspose(const ParticleSet& P, - int first, - int last, - ValueMatrix& logdet, - GradMatrix& dlogdet, - HessMatrix& grad_grad_logdet) final; - - void evaluate_notranspose(const ParticleSet& P, - int first, - int last, - ValueMatrix& logdet, - GradMatrix& dlogdet, - HessMatrix& grad_grad_logdet, - GGGMatrix& grad_grad_grad_logdet) final; - - //NOTE: The data types get complicated here, so here's an overview of the - // data types associated with ionic derivatives, and how to get their data. - // - //NOTE: These data structures hold the data for one particular ion, and so the ID is implicit. - // It's up to the user to keep track of which ion these derivatives refer to. - // - // 1.) GradMatrix grad_phi: Holds the ionic derivatives of each SPO for each electron. - // Example: grad_phi[iel][iorb][idim]. iel -- electron index. - // iorb -- orbital index. - // idim -- cartesian index of ionic derivative. - // X=0, Y=1, Z=2. - // - // 2.) HessMatrix grad_grad_phi: Holds the ionic derivatives of the electron gradient components - // for each SPO and each electron. - // Example: grad_grad_phi[iel][iorb](idim,edim) iel -- electron index. - // iorb -- orbital index. - // idim -- ionic derivative's cartesian index. - // X=0, Y=1, Z=2 - // edim -- electron derivative's cartesian index. - // x=0, y=1, z=2. - // - // 3.) GradMatrix grad_lapl_phi: Holds the ionic derivatives of the electron laplacian for each SPO and each electron. - // Example: grad_lapl_phi[iel][iorb][idim]. iel -- electron index. - // iorb -- orbital index. - // idim -- cartesian index of ionic derivative. - // X=0, Y=1, Z=2. - - /** - * \brief Calculate ion derivatives of SPO's. - * - * @param P Electron particle set. - * @param first index of first electron - * @@param last index of last electron - * @param source Ion particle set. - * @param iat_src Index of ion. - * @param gradphi Container storing ion gradients for all particles and all orbitals. - */ - void evaluateGradSource(const ParticleSet& P, - int first, - int last, - const ParticleSet& source, - int iat_src, - GradMatrix& grad_phi) final; - - /** - * \brief Calculate ion derivatives of SPO's, their gradients, and their laplacians. - * - * @param P Electron particle set. - * @param first index of first electron. - * @@param last index of last electron - * @param source Ion particle set. - * @param iat_src Index of ion. - * @param grad_phi Container storing ion gradients for all particles and all orbitals. - * @param grad_grad_phi Container storing ion gradients of electron gradients for all particles and all orbitals. - * @param grad_lapl_phi Container storing ion gradients of SPO laplacians for all particles and all orbitals. - */ - void evaluateGradSource(const ParticleSet& P, - int first, - int last, - const ParticleSet& source, - int iat_src, - GradMatrix& grad_phi, - HessMatrix& grad_grad_phi, - GradMatrix& grad_lapl_phi) final; - - void evaluateGradSourceRow(const ParticleSet& P, - int iel, - const ParticleSet& source, - int iat_src, - GradVector& grad_phi) final; - - void createResource(ResourceCollection& collection) const final; - void acquireResource(ResourceCollection& collection, const RefVectorWithLeader>& spo_list) const final; - void releaseResource(ResourceCollection& collection, const RefVectorWithLeader>& spo_list) const final; + LCAOrbitalSetT( + const std::string& my_name, std::unique_ptr&& bs); + + LCAOrbitalSetT(const LCAOrbitalSetT& in); + + std::string + getClassName() const final + { + return "LCAOrbitalSetT"; + } + + bool + isRotationSupported() const final + { + return true; + } + + bool + hasIonDerivs() const final + { + return true; + } + + std::unique_ptr> + makeClone() const final; + + void + storeParamsBeforeRotation() final + { + C_copy = std::make_shared(*C); + } + + void + applyRotation(const ValueMatrix& rot_mat, bool use_stored_copy) final; + + /** set the OrbitalSetSize and Identity=false and initialize internal + * storages + */ + void + setOrbitalSetSize(int norbs) final; + + /** return the size of the basis set + */ + int + getBasisSetSize() const + { + return (myBasisSet == nullptr) ? 0 : myBasisSet->getBasisSetSize(); + } + + bool + isIdentity() const + { + return Identity; + }; + + /** check consistency between Identity and C + * + */ + void + checkObject() const final; + + void + evaluateValue(const ParticleSet& P, int iat, ValueVector& psi) final; + + void + evaluateVGL(const ParticleSet& P, int iat, ValueVector& psi, + GradVector& dpsi, ValueVector& d2psi) final; + + void + mw_evaluateValue(const RefVectorWithLeader>& spo_list, + const RefVectorWithLeader& P_list, int iat, + const RefVector& psi_v_list) const final; + + void + mw_evaluateVGL(const RefVectorWithLeader>& spo_list, + const RefVectorWithLeader& P_list, int iat, + const RefVector& psi_v_list, + const RefVector& dpsi_v_list, + const RefVector& d2psi_v_list) const final; + + void + mw_evaluateDetRatios(const RefVectorWithLeader>& spo_list, + const RefVectorWithLeader& vp_list, + const RefVector& psi_list, + const std::vector& invRow_ptr_list, + std::vector>& ratios_list) const final; + + void + evaluateDetRatios(const VirtualParticleSet& VP, ValueVector& psi, + const ValueVector& psiinv, std::vector& ratios) final; + + void + mw_evaluateVGLandDetRatioGrads( + const RefVectorWithLeader>& spo_list, + const RefVectorWithLeader& P_list, int iat, + const std::vector& invRow_ptr_list, + OffloadMWVGLArray& phi_vgl_v, std::vector& ratios, + std::vector& grads) const final; + + void + evaluateVGH(const ParticleSet& P, int iat, ValueVector& psi, + GradVector& dpsi, HessVector& grad_grad_psi) final; + + void + evaluateVGHGH(const ParticleSet& P, int iat, ValueVector& psi, + GradVector& dpsi, HessVector& grad_grad_psi, + GGGVector& grad_grad_grad_psi) final; + + void + evaluate_notranspose(const ParticleSet& P, int first, int last, + ValueMatrix& logdet, GradMatrix& dlogdet, ValueMatrix& d2logdet) final; + + void + evaluate_notranspose(const ParticleSet& P, int first, int last, + ValueMatrix& logdet, GradMatrix& dlogdet, + HessMatrix& grad_grad_logdet) final; + + void + evaluate_notranspose(const ParticleSet& P, int first, int last, + ValueMatrix& logdet, GradMatrix& dlogdet, HessMatrix& grad_grad_logdet, + GGGMatrix& grad_grad_grad_logdet) final; + + // NOTE: The data types get complicated here, so here's an overview of the + // data types associated with ionic derivatives, and how to get their + // data. + // + // NOTE: These data structures hold the data for one particular ion, and so + // the ID is implicit. + // It's up to the user to keep track of which ion these derivatives + // refer to. + // + // 1.) GradMatrix grad_phi: Holds the ionic derivatives of each SPO for + // each electron. + // Example: grad_phi[iel][iorb][idim]. iel -- electron index. + // iorb -- orbital index. + // idim -- cartesian index + // of ionic derivative. + // X=0, Y=1, Z=2. + // + // 2.) HessMatrix grad_grad_phi: Holds the ionic derivatives of the + // electron gradient components + // for each SPO and each electron. + // Example: grad_grad_phi[iel][iorb](idim,edim) iel -- + // electron index. + // iorb -- + // orbital index. + // idim -- ionic + // derivative's + // cartesian + // index. + // X=0, Y=1, + // Z=2 + // edim -- + // electron + // derivative's + // cartesian + // index. + // x=0, y=1, + // z=2. + // + // 3.) GradMatrix grad_lapl_phi: Holds the ionic derivatives of the + // electron laplacian for each SPO and each electron. + // Example: grad_lapl_phi[iel][iorb][idim]. iel -- electron + // index. + // iorb -- orbital + // index. idim -- + // cartesian index of + // ionic derivative. + // X=0, Y=1, Z=2. + + /** + * \brief Calculate ion derivatives of SPO's. + * + * @param P Electron particle set. + * @param first index of first electron + * @@param last index of last electron + * @param source Ion particle set. + * @param iat_src Index of ion. + * @param gradphi Container storing ion gradients for all particles and all + * orbitals. + */ + void + evaluateGradSource(const ParticleSet& P, int first, int last, + const ParticleSet& source, int iat_src, GradMatrix& grad_phi) final; + + /** + * \brief Calculate ion derivatives of SPO's, their gradients, and their + * laplacians. + * + * @param P Electron particle set. + * @param first index of first electron. + * @@param last index of last electron + * @param source Ion particle set. + * @param iat_src Index of ion. + * @param grad_phi Container storing ion gradients for all particles and + * all orbitals. + * @param grad_grad_phi Container storing ion gradients of electron + * gradients for all particles and all orbitals. + * @param grad_lapl_phi Container storing ion gradients of SPO laplacians + * for all particles and all orbitals. + */ + void + evaluateGradSource(const ParticleSet& P, int first, int last, + const ParticleSet& source, int iat_src, GradMatrix& grad_phi, + HessMatrix& grad_grad_phi, GradMatrix& grad_lapl_phi) final; + + void + evaluateGradSourceRow(const ParticleSet& P, int iel, + const ParticleSet& source, int iat_src, GradVector& grad_phi) final; + + void + createResource(ResourceCollection& collection) const final; + void + acquireResource(ResourceCollection& collection, + const RefVectorWithLeader>& spo_list) const final; + void + releaseResource(ResourceCollection& collection, + const RefVectorWithLeader>& spo_list) const final; protected: - ///number of Single-particle orbitals - const IndexType BasisSetSize; - /// a copy of the original C before orbital rotation is applied; - std::shared_ptr C_copy; - - ///true if C is an identity matrix - bool Identity; - ///Temp(BasisSetSize) : Row index=V,Gx,Gy,Gz,L - vgl_type Temp; - ///Tempv(OrbitalSetSize) Tempv=C*Temp - vgl_type Tempv; - - ///These are temporary VectorSoAContainers to hold value, gradient, and hessian for - ///all basis or SPO functions evaluated at a given point. - ///Nbasis x [1(value)+3(gradient)+6(hessian)] - vgh_type Temph; - ///Norbitals x [1(value)+3(gradient)+6(hessian)] - vgh_type Temphv; - - ///These are temporary VectorSoAContainers to hold value, gradient, hessian, and - /// gradient hessian for all basis or SPO functions evaluated at a given point. - ///Nbasis x [1(value)+3(gradient)+6(hessian)+10(grad_hessian)] - vghgh_type Tempgh; - ///Nbasis x [1(value)+3(gradient)+6(hessian)+10(grad_hessian)] - vghgh_type Tempghv; + /// number of Single-particle orbitals + const IndexType BasisSetSize; + /// a copy of the original C before orbital rotation is applied; + std::shared_ptr C_copy; + + /// true if C is an identity matrix + bool Identity; + /// Temp(BasisSetSize) : Row index=V,Gx,Gy,Gz,L + vgl_type Temp; + /// Tempv(OrbitalSetSize) Tempv=C*Temp + vgl_type Tempv; + + /// These are temporary VectorSoAContainers to hold value, gradient, and + /// hessian for all basis or SPO functions evaluated at a given point. + /// Nbasis x [1(value)+3(gradient)+6(hessian)] + vgh_type Temph; + /// Norbitals x [1(value)+3(gradient)+6(hessian)] + vgh_type Temphv; + + /// These are temporary VectorSoAContainers to hold value, gradient, + /// hessian, and + /// gradient hessian for all basis or SPO functions evaluated at a given + /// point. + /// Nbasis x [1(value)+3(gradient)+6(hessian)+10(grad_hessian)] + vghgh_type Tempgh; + /// Nbasis x [1(value)+3(gradient)+6(hessian)+10(grad_hessian)] + vghgh_type Tempghv; private: - ///helper functions to handle Identity - void evaluate_vgl_impl(const vgl_type& temp, ValueVector& psi, GradVector& dpsi, ValueVector& d2psi) const; - - void evaluate_vgl_impl(const vgl_type& temp, - int i, - ValueMatrix& logdet, - GradMatrix& dlogdet, - ValueMatrix& d2logdet) const; - ///These two functions unpack the data in vgh_type temp object into wavefunction friendly data structures. - - - ///This unpacks temp into vectors psi, dpsi, and d2psi. - void evaluate_vgh_impl(const vgh_type& temp, ValueVector& psi, GradVector& dpsi, HessVector& d2psi) const; - - ///Unpacks temp into the ith row (or electron index) of logdet, dlogdet, dhlogdet. - void evaluate_vgh_impl(const vgh_type& temp, - int i, - ValueMatrix& logdet, - GradMatrix& dlogdet, - HessMatrix& dhlogdet) const; - ///Unpacks data in vghgh_type temp object into wavefunction friendly data structures for value, gradient, hessian - ///and gradient hessian. - void evaluate_vghgh_impl(const vghgh_type& temp, - ValueVector& psi, - GradVector& dpsi, - HessVector& d2psi, - GGGVector& dghpsi) const; - - void evaluate_vghgh_impl(const vghgh_type& temp, - int i, - ValueMatrix& logdet, - GradMatrix& dlogdet, - HessMatrix& dhlogdet, - GGGMatrix& dghlogdet) const; - - - ///Unpacks data in vgl object and calculates/places ionic gradient result into dlogdet. - void evaluate_ionderiv_v_impl(const vgl_type& temp, int i, GradMatrix& dlogdet) const; - - ///Unpacks data in vgl object and calculates/places ionic gradient of value, - /// electron gradient, and electron laplacian result into dlogdet, dglogdet, and dllogdet respectively. - void evaluate_ionderiv_vgl_impl(const vghgh_type& temp, - int i, - GradMatrix& dlogdet, - HessMatrix& dglogdet, - GradMatrix& dllogdet) const; - - ///Unpacks data in vgl object and calculates/places ionic gradient of a single row (phi_j(r)) into dlogdet. - void evaluate_ionderiv_v_row_impl(const vgl_type& temp, GradVector& dlogdet) const; - - void mw_evaluateVGLImplGEMM(const RefVectorWithLeader>& spo_list, - const RefVectorWithLeader& P_list, - int iat, - OffloadMWVGLArray& phi_vgl_v) const; - - /// packed walker GEMM implementation - void mw_evaluateValueImplGEMM(const RefVectorWithLeader>& spo_list, - const RefVectorWithLeader& P_list, - int iat, - OffloadMWVArray& phi_v) const; - - struct LCAOMultiWalkerMem; - ResourceHandle mw_mem_handle_; - /// timer for basis set - NewTimer& basis_timer_; - /// timer for MO - NewTimer& mo_timer_; + /// helper functions to handle Identity + void + evaluate_vgl_impl(const vgl_type& temp, ValueVector& psi, GradVector& dpsi, + ValueVector& d2psi) const; + + void + evaluate_vgl_impl(const vgl_type& temp, int i, ValueMatrix& logdet, + GradMatrix& dlogdet, ValueMatrix& d2logdet) const; + /// These two functions unpack the data in vgh_type temp object into + /// wavefunction friendly data structures. + + /// This unpacks temp into vectors psi, dpsi, and d2psi. + void + evaluate_vgh_impl(const vgh_type& temp, ValueVector& psi, GradVector& dpsi, + HessVector& d2psi) const; + + /// Unpacks temp into the ith row (or electron index) of logdet, dlogdet, + /// dhlogdet. + void + evaluate_vgh_impl(const vgh_type& temp, int i, ValueMatrix& logdet, + GradMatrix& dlogdet, HessMatrix& dhlogdet) const; + /// Unpacks data in vghgh_type temp object into wavefunction friendly data + /// structures for value, gradient, hessian and gradient hessian. + void + evaluate_vghgh_impl(const vghgh_type& temp, ValueVector& psi, + GradVector& dpsi, HessVector& d2psi, GGGVector& dghpsi) const; + + void + evaluate_vghgh_impl(const vghgh_type& temp, int i, ValueMatrix& logdet, + GradMatrix& dlogdet, HessMatrix& dhlogdet, GGGMatrix& dghlogdet) const; + + /// Unpacks data in vgl object and calculates/places ionic gradient result + /// into dlogdet. + void + evaluate_ionderiv_v_impl( + const vgl_type& temp, int i, GradMatrix& dlogdet) const; + + /// Unpacks data in vgl object and calculates/places ionic gradient of + /// value, + /// electron gradient, and electron laplacian result into dlogdet, + /// dglogdet, and dllogdet respectively. + void + evaluate_ionderiv_vgl_impl(const vghgh_type& temp, int i, + GradMatrix& dlogdet, HessMatrix& dglogdet, GradMatrix& dllogdet) const; + + /// Unpacks data in vgl object and calculates/places ionic gradient of a + /// single row (phi_j(r)) into dlogdet. + void + evaluate_ionderiv_v_row_impl( + const vgl_type& temp, GradVector& dlogdet) const; + + void + mw_evaluateVGLImplGEMM(const RefVectorWithLeader>& spo_list, + const RefVectorWithLeader& P_list, int iat, + OffloadMWVGLArray& phi_vgl_v) const; + + /// packed walker GEMM implementation + void + mw_evaluateValueImplGEMM(const RefVectorWithLeader>& spo_list, + const RefVectorWithLeader& P_list, int iat, + OffloadMWVArray& phi_v) const; + + struct LCAOMultiWalkerMem; + ResourceHandle mw_mem_handle_; + /// timer for basis set + NewTimer& basis_timer_; + /// timer for MO + NewTimer& mo_timer_; }; } // namespace qmcplusplus #endif diff --git a/src/QMCWaveFunctions/LCAO/LCAOrbitalSetWithCorrectionT.h b/src/QMCWaveFunctions/LCAO/LCAOrbitalSetWithCorrectionT.h index b1fc69cf6e..30c3f188e6 100644 --- a/src/QMCWaveFunctions/LCAO/LCAOrbitalSetWithCorrectionT.h +++ b/src/QMCWaveFunctions/LCAO/LCAOrbitalSetWithCorrectionT.h @@ -64,7 +64,8 @@ class LCAOrbitalSetWithCorrectionT : public SPOSetT GradMatrix& dlogdet, ValueMatrix& d2logdet) final; - friend class LCAOrbitalBuilder; + template + friend class LCAOrbitalBuilderT; private: LCAOrbitalSetT lcao; diff --git a/src/QMCWaveFunctions/LCAO/MultiFunctorAdapter.h b/src/QMCWaveFunctions/LCAO/MultiFunctorAdapter.h index 74df4c0b59..110491c006 100644 --- a/src/QMCWaveFunctions/LCAO/MultiFunctorAdapter.h +++ b/src/QMCWaveFunctions/LCAO/MultiFunctorAdapter.h @@ -19,6 +19,8 @@ #include "Message/MPIObjectBase.h" #include "ModernStringUtils.hpp" #include "hdf/hdf_archive.h" +#include "LCAO/MultiQuinticSpline1D.h" +#include "LCAO/SoaAtomicBasisSet.h" namespace qmcplusplus { diff --git a/src/QMCWaveFunctions/LCAO/SoaCuspCorrectionT.h b/src/QMCWaveFunctions/LCAO/SoaCuspCorrectionT.h index f20bfa5730..dca3912f90 100644 --- a/src/QMCWaveFunctions/LCAO/SoaCuspCorrectionT.h +++ b/src/QMCWaveFunctions/LCAO/SoaCuspCorrectionT.h @@ -12,8 +12,8 @@ /** @file SoaCuspCorrectionT.h */ -#ifndef QMCPLUSPLUS_SOA_CUSPCORRECTION_H -#define QMCPLUSPLUS_SOA_CUSPCORRECTION_H +#ifndef QMCPLUSPLUS_SOA_CUSPCORRECTIONT_H +#define QMCPLUSPLUS_SOA_CUSPCORRECTIONT_H #include "Configuration.h" #include "QMCWaveFunctions/SPOSetT.h" diff --git a/src/QMCWaveFunctions/LCAO/SoaLocalizedBasisSet.cpp b/src/QMCWaveFunctions/LCAO/SoaLocalizedBasisSet.cpp index 6f4afa9d53..eca9c70729 100644 --- a/src/QMCWaveFunctions/LCAO/SoaLocalizedBasisSet.cpp +++ b/src/QMCWaveFunctions/LCAO/SoaLocalizedBasisSet.cpp @@ -52,7 +52,7 @@ SoaLocalizedBasisSet::SoaLocalizedBasisSet(const SoaLocalizedBasisSet template void SoaLocalizedBasisSet::setPBCParams(const TinyVector& PBCImages, const TinyVector Sup_Twist, - const std::vector& phase_factor) + const std::vector& phase_factor) { for (int i = 0; i < LOBasisSet.size(); ++i) LOBasisSet[i]->setPBCParams(PBCImages, Sup_Twist, phase_factor); diff --git a/src/QMCWaveFunctions/LCAO/SoaLocalizedBasisSet.h b/src/QMCWaveFunctions/LCAO/SoaLocalizedBasisSet.h index 191ccdc8f4..4a9b8bf8f8 100644 --- a/src/QMCWaveFunctions/LCAO/SoaLocalizedBasisSet.h +++ b/src/QMCWaveFunctions/LCAO/SoaLocalizedBasisSet.h @@ -45,8 +45,8 @@ class SoaLocalizedBasisSet : public SoaBasisSetBase using vgh_type = typename BaseType::vgh_type; using vghgh_type = typename BaseType::vghgh_type; using PosType = typename ParticleSet::PosType; - using OffloadMWVGLArray = Array>; // [VGL, walker, Orbs] - using OffloadMWVArray = Array>; // [walker, Orbs] + using OffloadMWVGLArray = typename BaseType::OffloadMWVGLArray; + using OffloadMWVArray = typename BaseType::OffloadMWVArray; using BaseType::BasisSetSize; @@ -91,7 +91,7 @@ class SoaLocalizedBasisSet : public SoaBasisSetBase */ void setPBCParams(const TinyVector& PBCImages, const TinyVector Sup_Twist, - const std::vector& phase_factor); + const std::vector& phase_factor); /** set BasisSetSize and allocate mVGL container */ diff --git a/src/QMCWaveFunctions/SPOSetBuilderFactoryT.cpp b/src/QMCWaveFunctions/SPOSetBuilderFactoryT.cpp new file mode 100644 index 0000000000..284eaa7efe --- /dev/null +++ b/src/QMCWaveFunctions/SPOSetBuilderFactoryT.cpp @@ -0,0 +1,261 @@ +////////////////////////////////////////////////////////////////////////////////////// +// This file is distributed under the University of Illinois/NCSA Open Source +// License. See LICENSE file in top directory for details. +// +// Copyright (c) 2020 QMCPACK developers. +// +// File developed by: Bryan Clark, bclark@Princeton.edu, Princeton University +// Ken Esler, kpesler@gmail.com, University of Illinois at +// Urbana-Champaign Miguel Morales, moralessilva2@llnl.gov, +// Lawrence Livermore National Laboratory Jeremy McMinnis, +// jmcminis@gmail.com, University of Illinois at +// Urbana-Champaign Jaron T. Krogel, krogeljt@ornl.gov, Oak +// Ridge National Laboratory Jeongnim Kim, +// jeongnim.kim@gmail.com, University of Illinois at +// Urbana-Champaign Mark A. Berrill, berrillma@ornl.gov, Oak +// Ridge National Laboratory +// +// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois +// at Urbana-Champaign +////////////////////////////////////////////////////////////////////////////////////// + +#include "SPOSetBuilderFactoryT.h" + +#include "ModernStringUtils.hpp" +#include "QMCWaveFunctions/ElectronGas/FreeOrbitalBuilderT.h" +#include "QMCWaveFunctions/HarmonicOscillator/SHOSetBuilderT.h" +#include "QMCWaveFunctions/SPOSetScannerT.h" +#if OHMMS_DIM == 3 +#include "QMCWaveFunctions/LCAO/LCAOrbitalBuilderT.h" + +#if defined(QMC_COMPLEX) +#include "QMCWaveFunctions/EinsplineSpinorSetBuilder.h" +#include "QMCWaveFunctions/LCAO/LCAOSpinorBuilder.h" +#endif + +#if defined(HAVE_EINSPLINE) +#include "QMCWaveFunctions/EinsplineSetBuilder.h" +#endif +#endif +#include "Message/MPIObjectBase.h" +#include "OhmmsData/AttributeSet.h" +#include "QMCWaveFunctions/CompositeSPOSetT.h" +#include "Utilities/IteratorUtility.h" +#include "Utilities/ProgressReportEngine.h" + +namespace qmcplusplus +{ + +template +const SPOSetT* +SPOSetBuilderFactoryT::getSPOSet(const std::string& name) const +{ + if (auto spoit = sposets.find(name); spoit == sposets.end()) { + // keep this commented until legacy input styles are moved. + // In legacy input styles, this look up may fail and need to build + // SPOSetT on the fly. + return nullptr; + } + else + return spoit->second.get(); +} + +/** constructor + * \param els reference to the electrons + * \param psi reference to the wavefunction + * \param ions reference to the ions + */ +template +SPOSetBuilderFactoryT::SPOSetBuilderFactoryT( + Communicate* comm, ParticleSet& els, const PSetMap& psets) : + MPIObjectBase(comm), + targetPtcl(els), + ptclPool(psets) +{ + ClassName = "SPOSetBuilderFactoryT"; +} + +template +SPOSetBuilderFactoryT::~SPOSetBuilderFactoryT() +{ + DEBUG_MEMORY("SPOSetBuilderFactoryT::~SPOSetBuilderFactoryT"); +} + +template +std::unique_ptr> +SPOSetBuilderFactoryT::createSPOSetBuilder(xmlNodePtr rootNode) +{ + ReportEngine PRE(ClassName, "createSPOSetBuilder"); + std::string sourceOpt("ion0"); + std::string type(""); + std::string name(""); + OhmmsAttributeSet aAttrib; + aAttrib.add(sourceOpt, "source"); + aAttrib.add(type, "type"); + aAttrib.add(name, "name"); + + if (rootNode != NULL) + aAttrib.put(rootNode); + + std::string type_in = type; + type = lowerCase(type); + + // when name is missing, type becomes the input + if (name.empty()) + name = type_in; + + std::unique_ptr> bb; + + if (type == "composite") { + app_log() << "Composite SPO set with existing SPOSets." << std::endl; + bb = std::make_unique>(myComm, *this); + } + else if (type == "jellium" || type == "heg" || type == "free") { + app_log() << "Free-particle SPO set" << std::endl; + bb = std::make_unique>( + targetPtcl, myComm, rootNode); + } + else if (type == "sho") { + app_log() << "Harmonic Oscillator SPO set" << std::endl; + bb = std::make_unique>(targetPtcl, myComm); + } +#if OHMMS_DIM == 3 + else if (type.find("spline") < type.size()) { + if (targetPtcl.isSpinor()) { +#ifdef QMC_COMPLEX + app_log() << "Einspline Spinor Set\n"; + // FIXME + // bb = std::make_unique(targetPtcl, + // ptclPool, myComm, rootNode); +#else + PRE.error("Use of einspline spinors requires QMC_COMPLEX=1. " + "Rebuild with this option"); +#endif + } + else { +#if defined(HAVE_EINSPLINE) + PRE << "EinsplineSetBuilder: using libeinspline for B-spline " + "orbitals.\n"; + // FIXME + // bb = std::make_unique(targetPtcl, ptclPool, + // myComm, rootNode); +#else + PRE.error("Einspline is missing for B-spline orbitals", true); +#endif + } + } + else if (type == "molecularorbital" || type == "mo") { + ParticleSet* ions = nullptr; + // initialize with the source tag + auto pit(ptclPool.find(sourceOpt)); + if (pit == ptclPool.end()) + PRE.error("Missing basisset/@source.", true); + else + ions = pit->second.get(); + if (targetPtcl.isSpinor()) +#ifdef QMC_COMPLEX + bb = std::make_unique( + targetPtcl, *ions, myComm, rootNode); +#else + PRE.error("Use of lcao spinors requires QMC_COMPLEX=1. Rebuild " + "with this option"); +#endif + else + bb = std::make_unique>( + targetPtcl, *ions, myComm, rootNode); + } +#endif // OHMMS_DIM==3 + PRE.flush(); + + if (!bb) + myComm->barrier_and_abort("SPOSetBuilderFactoryT::createSPOSetBuilder " + "SPOSetBuilderT creation failed."); + + app_log() << " Created SPOSetT builder named '" << name << "' of type " + << type << std::endl; + return bb; +} + +template +void +SPOSetBuilderFactoryT::buildSPOSetCollection(xmlNodePtr cur) +{ + std::string collection_name; + std::string collection_type; + OhmmsAttributeSet attrib; + attrib.add(collection_name, "name"); + attrib.add(collection_type, "type"); + attrib.put(cur); + + // use collection_type as collection_name if collection_name is not given + if (collection_name.empty()) + collection_name = collection_type; + + app_summary() << std::endl; + app_summary() << " Single particle orbitals (SPO) collection" + << std::endl; + app_summary() << " -----------------------------------------" + << std::endl; + app_summary() << " Name: " << collection_name + << " Type input: " << collection_type << std::endl; + app_summary() << std::endl; + + // create the SPOSetT builder + auto bb = createSPOSetBuilder(cur); + + // going through a list of sposet entries + int nsposets = 0; + processChildren( + cur, [&](const std::string& cname, const xmlNodePtr element) { + if (cname == "sposet") { + addSPOSet( + std::unique_ptr>(bb->createSPOSet(element))); + nsposets++; + } + if (cname == "rotated_sposet") { + addSPOSet(std::unique_ptr>( + bb->createRotatedSPOSet(element))); + nsposets++; + } + }); + + if (nsposets == 0) + myComm->barrier_and_abort( + "SPOSetBuilderFactoryT::buildSPOSetCollection no " + "elements found"); + + // going through a list of spo_scanner entries + processChildren( + cur, [&](const std::string& cname, const xmlNodePtr element) { + if (cname == "spo_scanner") + if (myComm->rank() == 0) { + SPOSetScannerT ascanner(sposets, targetPtcl, ptclPool); + ascanner.put(element); + } + }); +} + +template +void +SPOSetBuilderFactoryT::addSPOSet(std::unique_ptr> spo) +{ + if (spo->getName().empty()) + myComm->barrier_and_abort( + "sposet created in sposet_collection must have a name!"); + + if (sposets.find(spo->getName()) != sposets.end()) + myComm->barrier_and_abort("The name of each sposet must be unique! '" + + spo->getName() + "' exists."); + else + sposets.emplace(spo->getName(), std::move(spo)); +} + +template +std::string SPOSetBuilderFactoryT::basisset_tag = "basisset"; + +template class SPOSetBuilderFactoryT; +template class SPOSetBuilderFactoryT; +template class SPOSetBuilderFactoryT>; +template class SPOSetBuilderFactoryT>; + +} // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/SPOSetBuilderFactoryT.h b/src/QMCWaveFunctions/SPOSetBuilderFactoryT.h new file mode 100644 index 0000000000..ce1e9b89da --- /dev/null +++ b/src/QMCWaveFunctions/SPOSetBuilderFactoryT.h @@ -0,0 +1,88 @@ +////////////////////////////////////////////////////////////////////////////////////// +// 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. +// +// File developed by: Jeremy McMinnis, jmcminis@gmail.com, University of +// Illinois at Urbana-Champaign +// Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National +// Laboratory Jeongnim Kim, jeongnim.kim@gmail.com, +// University of Illinois at Urbana-Champaign Mark A. +// Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory +// +// File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois +// at Urbana-Champaign +////////////////////////////////////////////////////////////////////////////////////// + +#ifndef QMCPLUSPLUS_BASISSETFACTORYT_H +#define QMCPLUSPLUS_BASISSETFACTORYT_H + +#include "QMCWaveFunctions/SPOSetBuilderT.h" +#include "QMCWaveFunctions/WaveFunctionComponentBuilder.h" +#include "type_traits/template_types.hpp" + +namespace qmcplusplus +{ +template +class SPOSetBuilderFactoryT : public MPIObjectBase +{ +public: + using SPOMap = typename SPOSetT::SPOMap; + using PSetMap = std::map>; + + /** constructor + * \param comm communicator + * \param els reference to the electrons + * \param ions reference to the ions + */ + SPOSetBuilderFactoryT( + Communicate* comm, ParticleSet& els, const PSetMap& psets); + + ~SPOSetBuilderFactoryT(); + + std::unique_ptr> + createSPOSetBuilder(xmlNodePtr rootNode); + + /** returns a named sposet from the pool + * only use in serial portion of execution + * ie during initialization prior to threaded code + */ + const SPOSetT* + getSPOSet(const std::string& name) const; + + void + buildSPOSetCollection(xmlNodePtr cur); + + bool + empty() const + { + return sposets.empty(); + } + + /** add an SPOSet to sposets map. + * This is only used to handle legacy SPOSet input styles without using + * sposet_collection + */ + void addSPOSet(std::unique_ptr>); + + SPOMap&& + exportSPOSets() + { + return std::move(sposets); + } + +private: + /// reference to the target particle + ParticleSet& targetPtcl; + + /// reference to the particle pool + const PSetMap& ptclPool; + + /// list of all sposets created by the builders of this factory + SPOMap sposets; + + static std::string basisset_tag; +}; +} // namespace qmcplusplus +#endif diff --git a/src/QMCWaveFunctions/SPOSetBuilderT.h b/src/QMCWaveFunctions/SPOSetBuilderT.h index 060451a94d..8bb3071df6 100644 --- a/src/QMCWaveFunctions/SPOSetBuilderT.h +++ b/src/QMCWaveFunctions/SPOSetBuilderT.h @@ -48,6 +48,7 @@ template class SPOSetBuilderT : public QMCTraits, public MPIObjectBase { public: + using PosType = typename SPOSetT::PosType; using RealType = typename SPOSetT::RealType; using indices_t = std::vector; using energies_t = std::vector; diff --git a/src/QMCWaveFunctions/SPOSetScannerT.h b/src/QMCWaveFunctions/SPOSetScannerT.h new file mode 100644 index 0000000000..9a3bb418a1 --- /dev/null +++ b/src/QMCWaveFunctions/SPOSetScannerT.h @@ -0,0 +1,216 @@ +////////////////////////////////////////////////////////////////////////////////////// +// 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. +// +// File developed by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory +// +// File created by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory +////////////////////////////////////////////////////////////////////////////////////// + + +#ifndef QMCPLUSPLUS_SPOSET_SCANNERT_H +#define QMCPLUSPLUS_SPOSET_SCANNERT_H + +#include "Particle/ParticleSet.h" +#include "QMCWaveFunctions/OrbitalSetTraits.h" +#include "QMCWaveFunctions/SPOSetT.h" +#include "OhmmsData/AttributeSet.h" + +namespace qmcplusplus +{ +/** a scanner for all the SPO sets. + */ +template +class SPOSetScannerT +{ +public: + using PtclPool = std::map>; + using SPOSetMap = typename SPOSetT::SPOMap; + using RealType = typename SPOSetT::RealType; + using ValueVector = typename SPOSetT::ValueVector; + using GradVector = typename SPOSetT::GradVector; + using HessVector = typename SPOSetT::HessVector; + + RealType myfabs(RealType s) { return std::fabs(s); } + template + std::complex myfabs(std::complex& s) + { + return std::complex(myfabs(s.real()), myfabs(s.imag())); + } + template + TinyVector myfabs(TinyVector& s) + { + return TinyVector(myfabs(s[0]), myfabs(s[1]), myfabs(s[2])); + } + + const SPOSetMap& sposets; + ParticleSet& target; + const PtclPool& ptcl_pool_; + ParticleSet* ions; + + // construction/destruction + SPOSetScannerT(const SPOSetMap& sposets_in, ParticleSet& targetPtcl, const PtclPool& psets) + : sposets(sposets_in), target(targetPtcl), ptcl_pool_(psets), ions(0){}; + //~SPOSetScannerT(){}; + + // processing scanning + void put(xmlNodePtr cur) + { + app_log() << "Entering the SPO set scanner!" << std::endl; + // check in the source particle set and search for it in the pool. + std::string sourcePtcl("ion0"); + OhmmsAttributeSet aAttrib; + aAttrib.add(sourcePtcl, "source"); + aAttrib.put(cur); + auto pit(ptcl_pool_.find(sourcePtcl)); + if (pit == ptcl_pool_.end()) + app_log() << "Source particle set not found. Can not be used as reference point." << std::endl; + else + ions = pit->second.get(); + + // scanning the SPO sets + xmlNodePtr cur_save = cur; + for (const auto& [name, sposet] : sposets) + { + app_log() << " Processing SPO " << sposet->getName() << std::endl; + // scanning the paths + cur = cur_save->children; + while (cur != NULL) + { + std::string trace_name("no name"); + OhmmsAttributeSet aAttrib; + aAttrib.add(trace_name, "name"); + aAttrib.put(cur); + std::string cname(getNodeName(cur)); + std::string prefix(sposet->getName() + "_" + cname + "_" + trace_name); + if (cname == "path") + { + app_log() << " Scanning a " << cname << " called " << trace_name << " and writing to " + << prefix + "_v/g/l/report.dat" << std::endl; + auto spo = sposet->makeClone(); + scan_path(cur, *spo, prefix); + } + else + { + if (cname != "text" && cname != "comment") + app_log() << " Unknown type of scanning " << cname << std::endl; + } + cur = cur->next; + } + } + app_log() << "Exiting the SPO set scanner!" << std::endl << std::endl; + } + + // scanning a path + void scan_path(xmlNodePtr cur, SPOSetT& sposet, std::string prefix) + { + std::string file_name; + file_name = prefix + "_v.dat"; + std::ofstream output_v(file_name.c_str()); + file_name = prefix + "_g.dat"; + std::ofstream output_g(file_name.c_str()); + file_name = prefix + "_l.dat"; + std::ofstream output_l(file_name.c_str()); + file_name = prefix + "_report.dat"; + std::ofstream output_report(file_name.c_str()); + + int nknots(2); + int from_atom(-1); + int to_atom(-1); + TinyVector from_pos(0.0, 0.0, 0.0); + TinyVector to_pos(0.0, 0.0, 0.0); + + OhmmsAttributeSet aAttrib; + aAttrib.add(nknots, "nknots"); + aAttrib.add(from_atom, "from_atom"); + aAttrib.add(to_atom, "to_atom"); + aAttrib.add(from_pos, "from_pos"); + aAttrib.add(to_pos, "to_pos"); + aAttrib.put(cur); + + // sanity check + if (nknots < 2) + nknots = 2; + // check out the reference atom coordinates + if (ions) + { + if (from_atom >= 0 && from_atom < ions->R.size()) + from_pos = ions->R[from_atom]; + if (to_atom >= 0 && to_atom < ions->R.size()) + to_pos = ions->R[to_atom]; + } + + // prepare a fake particle set + ValueVector SPO_v, SPO_l, SPO_v_avg, SPO_l_avg; + GradVector SPO_g, SPO_g_avg; + int OrbitalSize(sposet.size()); + SPO_v.resize(OrbitalSize); + SPO_g.resize(OrbitalSize); + SPO_l.resize(OrbitalSize); + SPO_v_avg.resize(OrbitalSize); + SPO_g_avg.resize(OrbitalSize); + SPO_l_avg.resize(OrbitalSize); + SPO_v_avg = 0.0; + SPO_g_avg = 0.0; + SPO_l_avg = 0.0; + double Delta = 1.0 / (nknots - 1); + int elec_count = target.R.size(); + auto R_saved = target.R; + ParticleSet::SingleParticlePos zero_pos(0.0, 0.0, 0.0); + for (int icount = 0, ind = 0; icount < nknots; icount++, ind++) + { + if (ind == elec_count) + ind = 0; + target.R[ind][0] = (to_pos[0] - from_pos[0]) * Delta * icount + from_pos[0]; + target.R[ind][1] = (to_pos[1] - from_pos[1]) * Delta * icount + from_pos[1]; + target.R[ind][2] = (to_pos[2] - from_pos[2]) * Delta * icount + from_pos[2]; + target.makeMove(ind, zero_pos); + sposet.evaluateVGL(target, ind, SPO_v, SPO_g, SPO_l); + std::ostringstream o; + o << "x_y_z " << std::fixed << std::setprecision(7) << target.R[ind][0] << " " << target.R[ind][1] << " " + << target.R[ind][2]; + output_v << o.str() << " : " << std::scientific << std::setprecision(12); + output_g << o.str() << " : " << std::scientific << std::setprecision(12); + output_l << o.str() << " : " << std::scientific << std::setprecision(12); + for (int iorb = 0; iorb < OrbitalSize; iorb++) + { + SPO_v_avg[iorb] += myfabs(SPO_v[iorb]); + SPO_g_avg[iorb] += myfabs(SPO_g[iorb]); + SPO_l_avg[iorb] += myfabs(SPO_l[iorb]); + output_v << SPO_v[iorb] << " "; + output_g << SPO_g[iorb][0] << " " << SPO_g[iorb][1] << " " << SPO_g[iorb][2] << " "; + output_l << SPO_l[iorb] << " "; + } + output_v << std::endl; + output_g << std::endl; + output_l << std::endl; + } + // restore the whole target. + target.R = R_saved; + target.update(); +#ifdef QMC_COMPLEX + output_report << "# Report: Orb Value_avg I/R Gradients_avg Laplacian_avg" << std::endl; +#else + output_report << "# Report: Orb Value_avg Gradients_avg Laplacian_avg" << std::endl; +#endif + for (int iorb = 0; iorb < OrbitalSize; iorb++) + output_report << "\t" << iorb << " " << std::scientific + << SPO_v_avg[iorb] * static_cast(1.0 / nknots) << " " +#ifdef QMC_COMPLEX + << SPO_v_avg[iorb].imag() / SPO_v_avg[iorb].real() << " " +#endif + << SPO_g_avg[iorb][0] * static_cast(1.0 / nknots) << " " + << SPO_g_avg[iorb][1] * static_cast(1.0 / nknots) << " " + << SPO_g_avg[iorb][2] * static_cast(1.0 / nknots) << " " + << SPO_l_avg[iorb] * static_cast(1.0 / nknots) << std::fixed << std::endl; + output_v.close(); + output_g.close(); + output_l.close(); + output_report.close(); + } +}; +} // namespace qmcplusplus + +#endif