From 7335258963ff3deccc27c0c9eeb1235ee630110e Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Tue, 28 Jun 2022 14:07:54 -0400 Subject: [PATCH 01/34] add free particle orbitals --- src/Particle/LongRange/StructFact.h | 3 + src/QMCWaveFunctions/CMakeLists.txt | 1 + .../ElectronGas/FreeParticle.cpp | 277 ++++++++++++++++++ .../ElectronGas/FreeParticle.h | 68 +++++ .../ElectronGas/FreeParticleBuilder.cpp | 47 +++ .../ElectronGas/FreeParticleBuilder.h | 19 ++ src/QMCWaveFunctions/SPOSetBuilderFactory.cpp | 7 +- 7 files changed, 419 insertions(+), 3 deletions(-) create mode 100644 src/QMCWaveFunctions/ElectronGas/FreeParticle.cpp create mode 100644 src/QMCWaveFunctions/ElectronGas/FreeParticle.h create mode 100644 src/QMCWaveFunctions/ElectronGas/FreeParticleBuilder.cpp create mode 100644 src/QMCWaveFunctions/ElectronGas/FreeParticleBuilder.h diff --git a/src/Particle/LongRange/StructFact.h b/src/Particle/LongRange/StructFact.h index 81278ec581..85da4a6308 100644 --- a/src/Particle/LongRange/StructFact.h +++ b/src/Particle/LongRange/StructFact.h @@ -71,6 +71,9 @@ class StructFact : public QMCTraits /// accessor of StorePerParticle bool isStorePerParticle() const { return StorePerParticle; } + /// accessor of k_lists_ + KContainer getKLists() const { return k_lists_; } + private: /// Compute all rhok elements from the start void computeRhok(const ParticleSet& P); diff --git a/src/QMCWaveFunctions/CMakeLists.txt b/src/QMCWaveFunctions/CMakeLists.txt index 9b79e0a7d4..7d4a617db8 100644 --- a/src/QMCWaveFunctions/CMakeLists.txt +++ b/src/QMCWaveFunctions/CMakeLists.txt @@ -63,6 +63,7 @@ set(JASTROW_SRCS set(JASTROW_OMPTARGET_SRCS Jastrow/J2OMPTarget.cpp Jastrow/BsplineFunctor.cpp) +set(FERMION_SRCS ${FERMION_SRCS} ElectronGas/FreeParticle.cpp ElectronGas/FreeParticleBuilder.cpp) if(QMC_COMPLEX) set(FERMION_SRCS ${FERMION_SRCS} ElectronGas/ElectronGasComplexOrbitalBuilder.cpp) else(QMC_COMPLEX) diff --git a/src/QMCWaveFunctions/ElectronGas/FreeParticle.cpp b/src/QMCWaveFunctions/ElectronGas/FreeParticle.cpp new file mode 100644 index 0000000000..24ef635efa --- /dev/null +++ b/src/QMCWaveFunctions/ElectronGas/FreeParticle.cpp @@ -0,0 +1,277 @@ +#include "FreeParticle.h" + +namespace qmcplusplus +{ + +FreeParticle::FreeParticle(const std::vector& kpts_cart) + : K(kpts_cart), +#ifdef QMC_COMPLEX + mink(0), // first k at twist may not be 0 +#else + mink(1), // treat k=0 as special case +#endif + maxk(kpts_cart.size()) +{ +#ifdef QMC_COMPLEX + OrbitalSetSize = maxk; +#else + OrbitalSetSize = 2*maxk-1; // k=0 has no (cos, sin) split +#endif + mK2.resize(maxk); + for (int ik=0; ik& kpts_cart); + ~FreeParticle(); + + // phi[i][j] is phi_j(r_i), i.e. electron i in orbital j + // i \in [first, last) + void evaluate_notranspose( + const ParticleSet& P, + int first, + int last, + ValueMatrix& phi, + GradMatrix& dphi, + ValueMatrix& d2phi) override; + + // plug r_i into all orbitals + void evaluateVGL( + const ParticleSet& P, + int i, + ValueVector& pvec, + GradVector& dpvec, + ValueVector& d2pvec + ) override; + void evaluateValue(const ParticleSet& P, int iat, ValueVector& pvec) override; + + // hessian matrix is needed by backflow + void evaluate_notranspose( + const ParticleSet& P, + int first, + int last, + ValueMatrix& phi, + GradMatrix& dphi, + HessMatrix& d2phi_mat) override; + + // derivative of hessian is needed to optimize backflow + void evaluate_notranspose( + const ParticleSet& P, + int first, + int last, + ValueMatrix& phi, + GradMatrix& dphi, + HessMatrix& d2phi_mat, + GGGMatrix& d3phi_mat) override; + + // ---- begin required overrides + std::unique_ptr makeClone() const {return std::make_unique(*this);} + void resetParameters(const opt_variables_type& optVariables) override {} //called by BFTrans} + void setOrbitalSetSize(int norbs) override {APP_ABORT("not implemented")}; + // required overrides end ---- + void report(const std::string& pad) const override; +private: + const std::vector K; // K vectors + const int mink; // minimum k index + const int maxk; // maximum number of K vectors + std::vector mK2; // minus K^2 +}; + +} // qmcplusplus +#endif diff --git a/src/QMCWaveFunctions/ElectronGas/FreeParticleBuilder.cpp b/src/QMCWaveFunctions/ElectronGas/FreeParticleBuilder.cpp new file mode 100644 index 0000000000..b162ce516c --- /dev/null +++ b/src/QMCWaveFunctions/ElectronGas/FreeParticleBuilder.cpp @@ -0,0 +1,47 @@ +#include "FreeParticleBuilder.h" +#include "OhmmsData/AttributeSet.h" +#include "LongRange/StructFact.h" +#include "FreeParticle.h" + +namespace qmcplusplus +{ + +FreeParticleBuilder::FreeParticleBuilder(ParticleSet& els, Communicate* comm, xmlNodePtr cur) + : SPOSetBuilder("PW", comm), + targetPtcl(els) +{} + +std::unique_ptr FreeParticleBuilder::createSPOSetFromXML(xmlNodePtr cur) +{ + int npw, norb; + PosType twist(0.0); + OhmmsAttributeSet attrib; + attrib.add(norb, "size"); + attrib.add(twist, "twist"); + attrib.put(cur); + +#ifdef QMC_COMPLEX + npw = norb; +#else + npw = std::ceil((norb+1.0)/2); +#endif + targetPtcl.setTwist(twist); + PosType tvec = targetPtcl.getLattice().k_cart(twist); + app_log() << "twist fraction = " << twist << std::endl; + app_log() << "twist cartesian = " << tvec << std::endl; + + // extract npw k-points from container + // kpts_cart is sorted by magnitude + std::vector kpts; + kpts.resize(npw); + kpts[0] = tvec; + for (int ik=1;ik(kpts); + sposet->report(" "); + return sposet; +} + +} // qmcplusplus diff --git a/src/QMCWaveFunctions/ElectronGas/FreeParticleBuilder.h b/src/QMCWaveFunctions/ElectronGas/FreeParticleBuilder.h new file mode 100644 index 0000000000..4dbfac3224 --- /dev/null +++ b/src/QMCWaveFunctions/ElectronGas/FreeParticleBuilder.h @@ -0,0 +1,19 @@ +#ifndef QMCPLUSPLUS_FREE_PARTICLE_BUILDER_H +#define QMCPLUSPLUS_FREE_PARTICLE_BUILDER_H + +#include "QMCWaveFunctions/SPOSetBuilder.h" + +namespace qmcplusplus +{ +class FreeParticleBuilder : public SPOSetBuilder +{ +public: + FreeParticleBuilder(ParticleSet& els, Communicate* comm, xmlNodePtr cur); + ~FreeParticleBuilder(){} + + std::unique_ptr createSPOSetFromXML(xmlNodePtr cur); +private: + ParticleSet& targetPtcl; +}; +} // qmcplusplus +#endif diff --git a/src/QMCWaveFunctions/SPOSetBuilderFactory.cpp b/src/QMCWaveFunctions/SPOSetBuilderFactory.cpp index 7b5428ced0..20e7bc816c 100644 --- a/src/QMCWaveFunctions/SPOSetBuilderFactory.cpp +++ b/src/QMCWaveFunctions/SPOSetBuilderFactory.cpp @@ -21,6 +21,7 @@ #include "QMCWaveFunctions/ElectronGas/ElectronGasOrbitalBuilder.h" #include "QMCWaveFunctions/HarmonicOscillator/SHOSetBuilder.h" #include "ModernStringUtils.hpp" +#include "QMCWaveFunctions/ElectronGas/FreeParticleBuilder.h" #if OHMMS_DIM == 3 #include "QMCWaveFunctions/LCAO/LCAOrbitalBuilder.h" @@ -95,10 +96,10 @@ std::unique_ptr SPOSetBuilderFactory::createSPOSetBuilder(xmlNode app_log() << "Composite SPO set with existing SPOSets." << std::endl; bb = std::make_unique(myComm, *this); } - else if (type == "jellium" || type == "heg") + else if (type == "jellium" || type == "heg" || type == "free") { - app_log() << "Electron gas SPO set" << std::endl; - bb = std::make_unique(targetPtcl, myComm, rootNode); + app_log() << "Free-particle SPO set" << std::endl; + bb = std::make_unique(targetPtcl, myComm, rootNode); } else if (type == "sho") { From e58fb6493bb693c4532c690e45accfe8141fd301 Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Tue, 28 Jun 2022 15:01:30 -0400 Subject: [PATCH 02/34] update heg tests --- tests/heg/heg_14_gamma/det_heg-short-SJB_opt.xml | 9 ++++++++- tests/heg/heg_14_gamma/heg.ni.wfs.xml | 10 +++++++++- tests/heg/heg_14_gamma/heg.sj.wfs.xml | 10 +++++++++- tests/heg/heg_14_gamma/heg.sjb.wfs.xml | 9 ++++++++- tests/heg/heg_54_J2rpa/heg.rpa.rs5.xml | 11 ++++++++++- 5 files changed, 44 insertions(+), 5 deletions(-) diff --git a/tests/heg/heg_14_gamma/det_heg-short-SJB_opt.xml b/tests/heg/heg_14_gamma/det_heg-short-SJB_opt.xml index 5bf607cce2..8fd090e047 100644 --- a/tests/heg/heg_14_gamma/det_heg-short-SJB_opt.xml +++ b/tests/heg/heg_14_gamma/det_heg-short-SJB_opt.xml @@ -9,7 +9,14 @@ Sample qmc run for Slater-Jastrow-Backflow HEG. - + + + + + + + + diff --git a/tests/heg/heg_14_gamma/heg.ni.wfs.xml b/tests/heg/heg_14_gamma/heg.ni.wfs.xml index d5754ecb9f..87ce64ae73 100644 --- a/tests/heg/heg_14_gamma/heg.ni.wfs.xml +++ b/tests/heg/heg_14_gamma/heg.ni.wfs.xml @@ -1,6 +1,14 @@ - + + + + + + + + + diff --git a/tests/heg/heg_14_gamma/heg.sj.wfs.xml b/tests/heg/heg_14_gamma/heg.sj.wfs.xml index a061d60059..eea423f6d1 100644 --- a/tests/heg/heg_14_gamma/heg.sj.wfs.xml +++ b/tests/heg/heg_14_gamma/heg.sj.wfs.xml @@ -1,7 +1,15 @@ - + + + + + + + + + 1.082858193 0.6653279375 0.4358910287 0.2243616172 0.1102948764 diff --git a/tests/heg/heg_14_gamma/heg.sjb.wfs.xml b/tests/heg/heg_14_gamma/heg.sjb.wfs.xml index f2c1d11570..f720bbacda 100644 --- a/tests/heg/heg_14_gamma/heg.sjb.wfs.xml +++ b/tests/heg/heg_14_gamma/heg.sjb.wfs.xml @@ -1,7 +1,14 @@ - + + + + + + + + diff --git a/tests/heg/heg_54_J2rpa/heg.rpa.rs5.xml b/tests/heg/heg_54_J2rpa/heg.rpa.rs5.xml index 8ad670491c..438e59b00f 100755 --- a/tests/heg/heg_54_J2rpa/heg.rpa.rs5.xml +++ b/tests/heg/heg_54_J2rpa/heg.rpa.rs5.xml @@ -23,7 +23,16 @@ xsi:noNamespaceSchemaLocation="http://www.mcc.uiuc.edu/qmc/schema/molecu.xsd"> - + + + + + + + + + + From 9e4a4d74a8a0e6c984e4ee808d9cd8e5dac70351 Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Tue, 28 Jun 2022 15:03:43 -0400 Subject: [PATCH 03/34] transfer old header also K -> kvecs --- .../ElectronGas/FreeParticle.cpp | 84 +++++++++---------- .../ElectronGas/FreeParticle.h | 22 ++++- 2 files changed, 61 insertions(+), 45 deletions(-) diff --git a/src/QMCWaveFunctions/ElectronGas/FreeParticle.cpp b/src/QMCWaveFunctions/ElectronGas/FreeParticle.cpp index 24ef635efa..71bd668dd6 100644 --- a/src/QMCWaveFunctions/ElectronGas/FreeParticle.cpp +++ b/src/QMCWaveFunctions/ElectronGas/FreeParticle.cpp @@ -4,7 +4,7 @@ namespace qmcplusplus { FreeParticle::FreeParticle(const std::vector& kpts_cart) - : K(kpts_cart), + : kvecs(kpts_cart), #ifdef QMC_COMPLEX mink(0), // first k at twist may not be 0 #else @@ -17,10 +17,10 @@ FreeParticle::FreeParticle(const std::vector& kpts_cart) #else OrbitalSetSize = 2*maxk-1; // k=0 has no (cos, sin) split #endif - mK2.resize(maxk); + k2neg.resize(maxk); for (int ik=0; ik K; // K vectors + const std::vector kvecs; // kvecs vectors const int mink; // minimum k index - const int maxk; // maximum number of K vectors - std::vector mK2; // minus K^2 + const int maxk; // maximum number of kvecs vectors + std::vector k2neg; // minus kvecs^2 }; } // qmcplusplus From 28bb03497d49c505cac10fe97879fdbb41d513ed Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Tue, 28 Jun 2022 15:23:04 -0400 Subject: [PATCH 04/34] remove ElectronGasOrbital and its Complex --- src/QMCHamiltonians/tests/test_ecp.cpp | 18 +- src/QMCWaveFunctions/CMakeLists.txt | 6 - .../ElectronGasComplexOrbitalBuilder.cpp | 127 ------- .../ElectronGasComplexOrbitalBuilder.h | 139 ------- .../ElectronGas/ElectronGasOrbitalBuilder.cpp | 172 --------- .../ElectronGas/ElectronGasOrbitalBuilder.h | 346 ------------------ src/QMCWaveFunctions/SPOSetBuilderFactory.cpp | 1 - src/QMCWaveFunctions/WaveFunctionFactory.cpp | 16 +- .../tests/test_DiracDeterminant.cpp | 23 +- tests/heg/heg_14_gamma/heg-long-HF.xml | 2 +- tests/heg/heg_14_gamma/heg-short-HF.xml | 2 +- tests/heg/heg_14_gamma/heg.hf.wfs.xml | 6 - tests/heg/heg_54_J2rpa/det_heg.rpa.rs5.xml | 10 +- tests/heg/heg_54_J2rpa/heg.rpa.rs5.xml | 0 14 files changed, 22 insertions(+), 846 deletions(-) delete mode 100644 src/QMCWaveFunctions/ElectronGas/ElectronGasComplexOrbitalBuilder.cpp delete mode 100644 src/QMCWaveFunctions/ElectronGas/ElectronGasComplexOrbitalBuilder.h delete mode 100644 src/QMCWaveFunctions/ElectronGas/ElectronGasOrbitalBuilder.cpp delete mode 100644 src/QMCWaveFunctions/ElectronGas/ElectronGasOrbitalBuilder.h delete mode 100644 tests/heg/heg_14_gamma/heg.hf.wfs.xml mode change 100755 => 100644 tests/heg/heg_54_J2rpa/det_heg.rpa.rs5.xml mode change 100755 => 100644 tests/heg/heg_54_J2rpa/heg.rpa.rs5.xml diff --git a/src/QMCHamiltonians/tests/test_ecp.cpp b/src/QMCHamiltonians/tests/test_ecp.cpp index 6e9e89489c..569afda532 100644 --- a/src/QMCHamiltonians/tests/test_ecp.cpp +++ b/src/QMCHamiltonians/tests/test_ecp.cpp @@ -35,9 +35,8 @@ #include "Particle/ParticleSet.h" #include "LongRange/EwaldHandler3D.h" -#ifdef QMC_COMPLEX //This is for the spinor test. -#include "QMCWaveFunctions/ElectronGas/ElectronGasComplexOrbitalBuilder.h" -#endif +//This is for the spinor test. +#include "QMCWaveFunctions/ElectronGas/FreeParticle.h" namespace qmcplusplus { @@ -501,20 +500,11 @@ TEST_CASE("Evaluate_soecp", "[hamiltonian]") kup.resize(nelec); kup[0] = PosType(1, 1, 1); - k2up.resize(nelec); - //For some goofy reason, EGOSet needs to be initialized with: - //1.) A k-vector list (fine). - //2.) A list of -|k|^2. To save on expensive - sign multiplication apparently. - k2up[0] = -dot(kup[0], kup[0]); - kdn.resize(nelec); kdn[0] = PosType(2, 2, 2); - k2dn.resize(nelec); - k2dn[0] = -dot(kdn[0], kdn[0]); - - auto spo_up = std::make_unique(kup, k2up); - auto spo_dn = std::make_unique(kdn, k2dn); + auto spo_up = std::make_unique(kup); + auto spo_dn = std::make_unique(kdn); auto spinor_set = std::make_unique(); spinor_set->set_spos(std::move(spo_up), std::move(spo_dn)); diff --git a/src/QMCWaveFunctions/CMakeLists.txt b/src/QMCWaveFunctions/CMakeLists.txt index 7d4a617db8..8024e8aea0 100644 --- a/src/QMCWaveFunctions/CMakeLists.txt +++ b/src/QMCWaveFunctions/CMakeLists.txt @@ -64,12 +64,6 @@ set(JASTROW_OMPTARGET_SRCS Jastrow/J2OMPTarget.cpp Jastrow/BsplineFunctor.cpp) set(FERMION_SRCS ${FERMION_SRCS} ElectronGas/FreeParticle.cpp ElectronGas/FreeParticleBuilder.cpp) -if(QMC_COMPLEX) - set(FERMION_SRCS ${FERMION_SRCS} ElectronGas/ElectronGasComplexOrbitalBuilder.cpp) -else(QMC_COMPLEX) - set(FERMION_SRCS ${FERMION_SRCS} ElectronGas/ElectronGasOrbitalBuilder.cpp) - -endif(QMC_COMPLEX) # wavefunctions only availbale to 3-dim problems if(OHMMS_DIM MATCHES 3) diff --git a/src/QMCWaveFunctions/ElectronGas/ElectronGasComplexOrbitalBuilder.cpp b/src/QMCWaveFunctions/ElectronGas/ElectronGasComplexOrbitalBuilder.cpp deleted file mode 100644 index 1fb29682ef..0000000000 --- a/src/QMCWaveFunctions/ElectronGas/ElectronGasComplexOrbitalBuilder.cpp +++ /dev/null @@ -1,127 +0,0 @@ -////////////////////////////////////////////////////////////////////////////////////// -// 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: Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign -// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign -// Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign -// Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory -// 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 "ElectronGasComplexOrbitalBuilder.h" -#include "QMCWaveFunctions/Fermion/SlaterDet.h" -#include "QMCWaveFunctions/Fermion/DiracDeterminant.h" -#include "OhmmsData/AttributeSet.h" - -namespace qmcplusplus -{ -/** constructor for EGOSet - * @param norb number of orbitals for the EGOSet - * @param k list of unique k points in Cartesian coordinate excluding gamma - * @param k2 k2[i]=dot(k[i],k[i]) - */ -EGOSet::EGOSet(const std::vector& k, const std::vector& k2) : K(k), mK2(k2) -{ - KptMax = k.size(); - OrbitalSetSize = k.size(); - className = "EGOSet"; - //assign_energies(); -} - -EGOSet::EGOSet(const std::vector& k, const std::vector& k2, const std::vector& d) - : K(k), mK2(k2) -{ - KptMax = k.size(); - OrbitalSetSize = k.size(); - className = "EGOSet"; - //assign_energies(); - //assign_degeneracies(d); -} - -ElectronGasComplexOrbitalBuilder::ElectronGasComplexOrbitalBuilder(Communicate* comm, ParticleSet& els) - : WaveFunctionComponentBuilder(comm, els) -{} - - -std::unique_ptr ElectronGasComplexOrbitalBuilder::buildComponent(xmlNodePtr cur) -{ - int nc = 0; - PosType twist(0.0); - OhmmsAttributeSet aAttrib; - aAttrib.add(nc, "shell"); - aAttrib.add(twist, "twist"); - aAttrib.put(cur); - //using Det_t = DiracDeterminant ; - //using SlaterDeterminant_t = SlaterDeterminant; - using Det_t = DiracDeterminant<>; - using SlaterDeterminant_t = SlaterDet; - int nat = targetPtcl.getTotalNum(); - int nup = nat / 2; - HEGGrid egGrid(targetPtcl.getLattice()); - if (nc == 0) - nc = egGrid.getShellIndex(nup); - egGrid.createGrid(nc, nup, twist); - targetPtcl.setTwist(twist); - std::vector> dets; - //create up determinant - dets.push_back(std::make_unique(std::make_unique(egGrid.kpt, egGrid.mk2), 0, nup)); - //create down determinant - dets.push_back(std::make_unique(std::make_unique(egGrid.kpt, egGrid.mk2), nup, nup + nup)); - //create a Slater determinant - return std::make_unique(targetPtcl, std::move(dets)); -} - -ElectronGasSPOBuilder::ElectronGasSPOBuilder(ParticleSet& p, Communicate* comm, xmlNodePtr cur) - : SPOSetBuilder("ElectronGas", comm), has_twist(false), unique_twist(-1.0), egGrid(p.getLattice()), spo_node(NULL) -{ - ClassName = "ElectronGasSPOBuilder"; -} - -std::unique_ptr ElectronGasSPOBuilder::createSPOSetFromXML(xmlNodePtr cur) -{ - app_log() << "ElectronGasSPOBuilder::createSPOSet " << std::endl; - int nc = 0; - int ns = 0; - PosType twist(0.0); - std::string spo_name("heg"); - OhmmsAttributeSet aAttrib; - aAttrib.add(ns, "size"); - aAttrib.add(twist, "twist"); - aAttrib.add(spo_name, "name"); - aAttrib.add(spo_name, "id"); - aAttrib.put(cur); - if (has_twist) - twist = unique_twist; - else - { - unique_twist = twist; - has_twist = true; - for (int d = 0; d < OHMMS_DIM; ++d) - has_twist &= (unique_twist[d] + 1.0) > 1e-6; - } - if (ns > 0) - nc = egGrid.getShellIndex(ns); - if (nc < 0) - { - app_error() << " HEG Invalid Shell." << std::endl; - APP_ABORT("ElectronGasSPOBuilder::put"); - } - egGrid.createGrid(nc, ns, twist); - return std::make_unique(egGrid.kpt, egGrid.mk2, egGrid.deg); -} - - -std::unique_ptr ElectronGasSPOBuilder::createSPOSetFromIndices(indices_t& indices) -{ - egGrid.createGrid(indices); - return std::make_unique(egGrid.kpt, egGrid.mk2, egGrid.deg); -} - - -} // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/ElectronGas/ElectronGasComplexOrbitalBuilder.h b/src/QMCWaveFunctions/ElectronGas/ElectronGasComplexOrbitalBuilder.h deleted file mode 100644 index 443c97a6df..0000000000 --- a/src/QMCWaveFunctions/ElectronGas/ElectronGasComplexOrbitalBuilder.h +++ /dev/null @@ -1,139 +0,0 @@ -////////////////////////////////////////////////////////////////////////////////////// -// 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: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign -// Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign -// Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory -// 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_ELECTRONGAS_COMPLEXORBITALS_H -#define QMCPLUSPLUS_ELECTRONGAS_COMPLEXORBITALS_H - -#include "QMCWaveFunctions/WaveFunctionComponentBuilder.h" -#include "QMCWaveFunctions/SPOSet.h" -#include "QMCWaveFunctions/SPOSetBuilder.h" -#include "QMCWaveFunctions/ElectronGas/HEGGrid.h" -#include "CPU/math.hpp" - - -namespace qmcplusplus -{ -struct EGOSet : public SPOSet -{ - int KptMax; - std::vector K; - std::vector mK2; - - EGOSet(const std::vector& k, const std::vector& k2); - EGOSet(const std::vector& k, const std::vector& k2, const std::vector& d); - - std::unique_ptr makeClone() const override { return std::make_unique(*this); } - - void resetParameters(const opt_variables_type& optVariables) override {} - void setOrbitalSetSize(int norbs) override {} - - inline void evaluateValue(const ParticleSet& P, int iat, ValueVector& psi) override - { - const PosType& r = P.activeR(iat); - RealType sinkr, coskr; - for (int ik = 0; ik < KptMax; ik++) - { - qmcplusplus::sincos(dot(K[ik], r), &sinkr, &coskr); - psi[ik] = ValueType(coskr, sinkr); - } - } - - /** generic inline function to handle a row - * @param P current ParticleSet - * @param iat active particle - * @param psi value row - * @param dpsi gradient row - * @param d2psi laplacian row - */ - inline void evaluateVGL(const ParticleSet& P, - int iat, - ValueVector& psi, - GradVector& dpsi, - ValueVector& d2psi) override - { - const PosType& r = P.activeR(iat); - RealType sinkr, coskr; - for (int ik = 0; ik < KptMax; ik++) - { - qmcplusplus::sincos(dot(K[ik], r), &sinkr, &coskr); - psi[ik] = ValueType(coskr, sinkr); - dpsi[ik] = ValueType(-sinkr, coskr) * K[ik]; - d2psi[ik] = ValueType(mK2[ik] * coskr, mK2[ik] * sinkr); - } - } - - void evaluate_notranspose(const ParticleSet& P, - int first, - int last, - ValueMatrix& logdet, - GradMatrix& dlogdet, - ValueMatrix& d2logdet) override - { - for (int iat = first, i = 0; iat < last; ++iat, ++i) - { - ValueVector v(logdet[i], OrbitalSetSize); - GradVector g(dlogdet[i], OrbitalSetSize); - ValueVector l(d2logdet[i], OrbitalSetSize); - evaluateVGL(P, iat, v, g, l); - } - } - void evaluate_notranspose(const ParticleSet& P, - int first, - int last, - ValueMatrix& logdet, - GradMatrix& dlogdet, - HessMatrix& grad_grad_logdet, - GGGMatrix& grad_grad_grad_logdet) override - { - APP_ABORT( - "Incomplete implementation EGOSet::evaluate(P,first,last,lodget,dlodet,grad_grad_logdet,grad_grad_grad_logdet"); - } -}; - - -/** OrbitalBuilder for Slater determinants of electron-gas -*/ -class ElectronGasComplexOrbitalBuilder : public WaveFunctionComponentBuilder -{ -public: - ///constructor - ElectronGasComplexOrbitalBuilder(Communicate* comm, ParticleSet& els); - - ///implement vritual function - std::unique_ptr buildComponent(xmlNodePtr cur) override; -}; - -/** OrbitalBuilder for Slater determinants of electron-gas -*/ -class ElectronGasSPOBuilder : public SPOSetBuilder -{ -protected: - bool has_twist; - PosType unique_twist; - HEGGrid egGrid; - xmlNodePtr spo_node; - -public: - ///constructor - ElectronGasSPOBuilder(ParticleSet& p, Communicate* comm, xmlNodePtr cur); - - /** initialize the Antisymmetric wave function for electrons - *@param cur the current xml node - */ - std::unique_ptr createSPOSetFromXML(xmlNodePtr cur) override; - std::unique_ptr createSPOSetFromIndices(indices_t& indices); -}; -} // namespace qmcplusplus -#endif diff --git a/src/QMCWaveFunctions/ElectronGas/ElectronGasOrbitalBuilder.cpp b/src/QMCWaveFunctions/ElectronGas/ElectronGasOrbitalBuilder.cpp deleted file mode 100644 index 3f15d1e4f0..0000000000 --- a/src/QMCWaveFunctions/ElectronGas/ElectronGasOrbitalBuilder.cpp +++ /dev/null @@ -1,172 +0,0 @@ -////////////////////////////////////////////////////////////////////////////////////// -// 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: 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 -// Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign -// Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory -// 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 "ElectronGasOrbitalBuilder.h" -#include "QMCWaveFunctions/Fermion/SlaterDet.h" -#include "OhmmsData/AttributeSet.h" -#include "QMCWaveFunctions/Fermion/BackflowBuilder.h" -#include "QMCWaveFunctions/Fermion/BackflowTransformation.h" -#include "QMCWaveFunctions/Fermion/SlaterDetWithBackflow.h" -#include "QMCWaveFunctions/Fermion/DiracDeterminant.h" -#include "QMCWaveFunctions/Fermion/DiracDeterminantWithBackflow.h" - -namespace qmcplusplus -{ -/** constructor for EGOSet - * @param norb number of orbitals for the EGOSet - * @param k list of unique k points in Cartesian coordinate excluding gamma - * @param k2 k2[i]=dot(k[i],k[i]) - */ -RealEGOSet::RealEGOSet(const std::vector& k, const std::vector& k2) : K(k), mK2(k2) -{ - KptMax = k.size(); - OrbitalSetSize = 2 * k.size() + 1; - className = "EGOSet"; -} - -ElectronGasOrbitalBuilder::ElectronGasOrbitalBuilder(Communicate* comm, ParticleSet& els) - : WaveFunctionComponentBuilder(comm, els), UseBackflow(false) -{} - -std::unique_ptr ElectronGasOrbitalBuilder::buildComponent(xmlNodePtr cur) -{ - int nc(0), nc2(-2); - ValueType bosonic_eps(-999999); - ValueType rntype(0); - PosType twist(0.0); - OhmmsAttributeSet aAttrib; - aAttrib.add(nc, "shell"); - aAttrib.add(nc2, "shell2"); - aAttrib.add(bosonic_eps, "eps"); - aAttrib.add(rntype, "primary"); - aAttrib.add(twist, "twist"); - aAttrib.put(cur); - if (nc2 == -2) - nc2 = nc; - xmlNodePtr curRoot = cur; - xmlNodePtr BFNode(NULL); - cur = curRoot->children; - while (cur != NULL) //check the basis set - { - std::string cname(getNodeName(cur)); - if (cname == backflow_tag) - { - // FIX FIX FIX !!! - UseBackflow = true; - if (BFNode == NULL) - BFNode = cur; - } - cur = cur->next; - } - HEGGrid egGrid(targetPtcl.getLattice()); - int nat = targetPtcl.getTotalNum(); - if (nc == 0) - nc = nc2 = egGrid.getShellIndex(nat / 2); - int nup = egGrid.getNumberOfKpoints(nc); - int ndn(0); - if (nc2 > -1) - ndn = egGrid.getNumberOfKpoints(nc2); - if (nc < 0) - { - app_error() << " HEG Invalid Shell." << std::endl; - APP_ABORT("ElectronGasOrbitalBuilder::put"); - } - if (nat != (nup + ndn)) - { - app_error() << " The number of particles " << nup << "/" << ndn << " does not match to the shell." << std::endl; - app_error() << " Suggested values for the number of particles " << std::endl; - app_error() << " " << 2 * egGrid.getNumberOfKpoints(nc) << " for shell " << nc << std::endl; - app_error() << " " << 2 * egGrid.getNumberOfKpoints(nc - 1) << " for shell " << nc - 1 << std::endl; - APP_ABORT("ElectronGasOrbitalBuilder::put"); - return nullptr; - } - - //create a E(lectron)G(as)O(rbital)Set - int nkpts = (nup - 1) / 2; - egGrid.createGrid(nc, nkpts); - auto psiu = std::make_unique(egGrid.kpt, egGrid.mk2); - std::unique_ptr psid; - if (ndn > 0) - { - if (nup != ndn) - { - int nkpts2 = (ndn - 1) / 2; - HEGGrid egGrid2(targetPtcl.getLattice()); - egGrid2.createGrid(nc2, nkpts2); - } - psid = std::make_unique(egGrid.kpt, egGrid.mk2); - } - - //create a Slater determinant - if (UseBackflow) - { - app_log() << "Creating Backflow transformation in ElectronGasOrbitalBuilder::put(xmlNodePtr cur).\n"; - PSetMap dummy; - BackflowBuilder bfbuilder(targetPtcl, dummy); - auto BFTrans = bfbuilder.buildBackflowTransformation(BFNode); - std::vector> dets; - //create up determinant - dets.push_back(std::make_unique(std::move(psiu), *BFTrans, 0, nup)); - //create down determinant - if (ndn > 0) - dets.push_back(std::make_unique(std::move(psid), *BFTrans, nup, nup + ndn)); - auto sdet = std::make_unique(targetPtcl, std::move(dets), std::move(BFTrans)); - return sdet; - } - else - { - std::vector> dets; - //create up determinant - dets.push_back(std::make_unique>(std::move(psiu), 0, nup)); - //create down determinant - if (ndn > 0) - dets.push_back(std::make_unique>(std::move(psid), nup, nup + ndn)); - return std::make_unique(targetPtcl, std::move(dets)); - } -} - -ElectronGasSPOBuilder::ElectronGasSPOBuilder(ParticleSet& p, Communicate* comm, xmlNodePtr cur) - : SPOSetBuilder("ElectronGas", comm), egGrid(p.getLattice()) -{ - ClassName = "ElectronGasSPOBuilder"; -} - -std::unique_ptr ElectronGasSPOBuilder::createSPOSetFromXML(xmlNodePtr cur) -{ - int nc = 0; - int ns = 0; - PosType twist(0.0); - std::string spo_name("heg"); - OhmmsAttributeSet aAttrib; - aAttrib.add(ns, "size"); - aAttrib.add(twist, "twist"); - aAttrib.add(spo_name, "name"); - aAttrib.add(spo_name, "id"); - aAttrib.put(cur); - - if (ns > 0) - nc = egGrid.getShellFromStates(ns); - if (nc < 0) - { - app_error() << " HEG Invalid Shell." << std::endl; - APP_ABORT("ElectronGasOrbitalBuilder::put"); - } - egGrid.createGrid(nc, (ns - 1) / 2); - return std::make_unique(egGrid.kpt, egGrid.mk2); -} - -} // namespace qmcplusplus diff --git a/src/QMCWaveFunctions/ElectronGas/ElectronGasOrbitalBuilder.h b/src/QMCWaveFunctions/ElectronGas/ElectronGasOrbitalBuilder.h deleted file mode 100644 index 815ea964e7..0000000000 --- a/src/QMCWaveFunctions/ElectronGas/ElectronGasOrbitalBuilder.h +++ /dev/null @@ -1,346 +0,0 @@ -////////////////////////////////////////////////////////////////////////////////////// -// 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 -// Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign -// Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory -// 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_ELECTRONGAS_ORBITALS_H -#define QMCPLUSPLUS_ELECTRONGAS_ORBITALS_H - -#include "QMCWaveFunctions/WaveFunctionComponentBuilder.h" -#include "QMCWaveFunctions/SPOSet.h" -#include "CPU/math.hpp" - -#include "QMCWaveFunctions/SPOSetBuilder.h" -#include "QMCWaveFunctions/ElectronGas/HEGGrid.h" - -#if defined(QMC_COMPLEX) -#include "QMCWaveFunctions/ElectronGas/ElectronGasComplexOrbitalBuilder.h" -#else /** declare real HEG orbitals **/ -namespace qmcplusplus -{ -//forward declaration -class BackflowTransformation; - -struct RealEGOSet : public SPOSet -{ - int KptMax; - RealType kdotr; - std::vector K; - std::vector mK2; - - RealEGOSet(const std::vector& k, const std::vector& k2); - - void resetParameters(const opt_variables_type& optVariables) override {} - void setOrbitalSetSize(int norbs) override {} - - std::unique_ptr makeClone() const override { return std::make_unique(*this); } - - PosType get_k(int i) override - { - //Only used in the same_k part of the optimizable SPO set. we allow optimization to k points in the same direction - if (i > 0) - { - int ik = (i - 1) / 2; - // int even=(i-1)%2; - PosType k_tmp = K[ik]; - k_tmp *= 1.0 / std::sqrt(-mK2[ik]); - // if(even) - // return -1*k_tmp; - // else - return k_tmp; - } - else - return PosType(); - }; - - inline ValueType f(const PosType& pos, int i) - { - if (i > 0) - { - int ik = (i - 1) / 2; - int even = (i - 1) % 2; - kdotr = dot(K[ik], pos); - if (even) - return std::cos(kdotr); - else - return std::sin(kdotr); - } - else - return 1.0; - } - - void evaluateValue(const ParticleSet& P, int iat, ValueVector& psi) override - { - const PosType& r = P.activeR(iat); - RealType sinkr, coskr; - psi[0] = 1.0; - for (int ik = 0, j = 1; ik < KptMax; ik++) - { - qmcplusplus::sincos(dot(K[ik], r), &sinkr, &coskr); - psi[j++] = coskr; - psi[j++] = sinkr; - } - } - - /** generic inline function to handle a row - * @param P current ParticleSet - * @param iat active particle - * @param psi value row - * @param dpsi gradient row - * @param d2psi laplacian row - */ - inline void evaluateVGL(const ParticleSet& P, - int iat, - ValueVector& psi, - GradVector& dpsi, - ValueVector& d2psi) override - { - psi[0] = 1.0; - dpsi[0] = 0.0; - d2psi[0] = 0.0; - const PosType& r = P.activeR(iat); - RealType coskr, sinkr; - for (int ik = 0, j1 = 1; ik < KptMax; ik++, j1 += 2) - { - int j2 = j1 + 1; - qmcplusplus::sincos(dot(K[ik], r), &sinkr, &coskr); - psi[j1] = coskr; - psi[j2] = sinkr; - dpsi[j1] = -sinkr * K[ik]; - dpsi[j2] = coskr * K[ik]; - d2psi[j1] = mK2[ik] * coskr; - d2psi[j2] = mK2[ik] * sinkr; - } - } - - /** generic inline function to handle a row - * @param P current ParticleSet - * @param iat active particle - * @param psi value row - * @param dpsi gradient row - * @param hess hessian row - */ - inline void evaluateVGH(const ParticleSet& P, int iat, ValueVector& psi, GradVector& dpsi, HessVector& hess) override - { - psi[0] = 1.0; - dpsi[0] = 0.0; - hess[0] = 0.0; - const PosType& r = P.activeR(iat); - RealType coskr, sinkr; - for (int ik = 0, j1 = 1; ik < KptMax; ik++, j1 += 2) - { - int j2 = j1 + 1; - qmcplusplus::sincos(dot(K[ik], r), &sinkr, &coskr); - psi[j1] = coskr; - psi[j2] = sinkr; - dpsi[j1] = -sinkr * K[ik]; - dpsi[j2] = coskr * K[ik]; - for (int la = 0; la < OHMMS_DIM; la++) - { - (hess[j1])(la, la) = -coskr * (K[ik])[la] * (K[ik])[la]; - (hess[j2])(la, la) = -sinkr * (K[ik])[la] * (K[ik])[la]; - for (int lb = la + 1; lb < OHMMS_DIM; lb++) - { - (hess[j1])(la, lb) = -coskr * (K[ik])[la] * (K[ik])[lb]; - (hess[j2])(la, lb) = -sinkr * (K[ik])[la] * (K[ik])[lb]; - (hess[j1])(lb, la) = (hess[j1])(la, lb); - (hess[j2])(lb, la) = (hess[j2])(la, lb); - } - } - } - } - - void evaluate_notranspose(const ParticleSet& P, - int first, - int last, - ValueMatrix& logdet, - GradMatrix& dlogdet, - ValueMatrix& d2logdet) override - { - for (int iat = first, i = 0; iat < last; ++iat, ++i) - { - ValueVector v(logdet[i], OrbitalSetSize); - GradVector g(dlogdet[i], OrbitalSetSize); - ValueVector l(d2logdet[i], OrbitalSetSize); - evaluateVGL(P, iat, v, g, l); - } - } - - void evaluate_notranspose(const ParticleSet& P, - int first, - int last, - ValueMatrix& logdet, - GradMatrix& dlogdet, - HessMatrix& grad_grad_logdet) override - { - for (int i = 0, iat = first; iat < last; i++, iat++) - { - ValueType* psi = logdet[i]; - GradType* dpsi = dlogdet[i]; - HessType* hess = grad_grad_logdet[i]; - psi[0] = 1.0; - dpsi[0] = 0.0; - hess[0] = 0.0; - RealType coskr, sinkr; - for (int ik = 0, j1 = 1; ik < KptMax; ik++, j1 += 2) - { - int j2 = j1 + 1; - qmcplusplus::sincos(dot(K[ik], P.R[iat]), &sinkr, &coskr); - psi[j1] = coskr; - psi[j2] = sinkr; - dpsi[j1] = -sinkr * K[ik]; - dpsi[j2] = coskr * K[ik]; - for (int la = 0; la < OHMMS_DIM; la++) - { - (hess[j1])(la, la) = -coskr * (K[ik])[la] * (K[ik])[la]; - (hess[j2])(la, la) = -sinkr * (K[ik])[la] * (K[ik])[la]; - for (int lb = la + 1; lb < OHMMS_DIM; lb++) - { - (hess[j1])(la, lb) = -coskr * (K[ik])[la] * (K[ik])[lb]; - (hess[j2])(la, lb) = -sinkr * (K[ik])[la] * (K[ik])[lb]; - (hess[j1])(lb, la) = (hess[j1])(la, lb); - (hess[j2])(lb, la) = (hess[j2])(la, lb); - } - } - } - } - } - - void evaluate_notranspose(const ParticleSet& P, - int first, - int last, - ValueMatrix& logdet, - GradMatrix& dlogdet, - HessMatrix& grad_grad_logdet, - GGGMatrix& grad_grad_grad_logdet) override - { - for (int i = 0, iat = first; iat < last; i++, iat++) - { - ValueType* psi = logdet[i]; - GradType* dpsi = dlogdet[i]; - HessType* hess = grad_grad_logdet[i]; - GGGType* ggg = grad_grad_grad_logdet[i]; - psi[0] = 1.0; - dpsi[0] = 0.0; - hess[0] = 0.0; - ggg[0] = 0.0; - RealType coskr, sinkr; - for (int ik = 0, j1 = 1; ik < KptMax; ik++, j1 += 2) - { - int j2 = j1 + 1; - qmcplusplus::sincos(dot(K[ik], P.R[iat]), &sinkr, &coskr); - psi[j1] = coskr; - psi[j2] = sinkr; - dpsi[j1] = -sinkr * K[ik]; - dpsi[j2] = coskr * K[ik]; - for (int la = 0; la < OHMMS_DIM; la++) - { - (hess[j1])(la, la) = -coskr * (K[ik])[la] * (K[ik])[la]; - (hess[j2])(la, la) = -sinkr * (K[ik])[la] * (K[ik])[la]; - ((ggg[j1])[la])(la, la) = sinkr * (K[ik])[la] * (K[ik])[la] * (K[ik])[la]; - ((ggg[j2])[la])(la, la) = -coskr * (K[ik])[la] * (K[ik])[la] * (K[ik])[la]; - for (int lb = la + 1; lb < OHMMS_DIM; lb++) - { - (hess[j1])(la, lb) = -coskr * (K[ik])[la] * (K[ik])[lb]; - (hess[j2])(la, lb) = -sinkr * (K[ik])[la] * (K[ik])[lb]; - (hess[j1])(lb, la) = (hess[j1])(la, lb); - (hess[j2])(lb, la) = (hess[j2])(la, lb); - ((ggg[j1])[la])(lb, la) = sinkr * (K[ik])[la] * (K[ik])[lb] * (K[ik])[la]; - ((ggg[j2])[la])(lb, la) = -coskr * (K[ik])[la] * (K[ik])[lb] * (K[ik])[la]; - ((ggg[j1])[la])(la, lb) = ((ggg[j1])[la])(lb, la); - ((ggg[j2])[la])(la, lb) = ((ggg[j2])[la])(lb, la); - ((ggg[j1])[lb])(la, la) = ((ggg[j1])[la])(lb, la); - ((ggg[j2])[lb])(la, la) = ((ggg[j2])[la])(lb, la); - ((ggg[j1])[la])(lb, lb) = sinkr * (K[ik])[la] * (K[ik])[lb] * (K[ik])[lb]; - ((ggg[j2])[la])(lb, lb) = -coskr * (K[ik])[la] * (K[ik])[lb] * (K[ik])[lb]; - ((ggg[j1])[lb])(la, lb) = ((ggg[j1])[la])(lb, lb); - ((ggg[j2])[lb])(la, lb) = ((ggg[j2])[la])(lb, lb); - ((ggg[j1])[lb])(lb, la) = ((ggg[j1])[la])(lb, lb); - ((ggg[j2])[lb])(lb, la) = ((ggg[j2])[la])(lb, lb); - for (int lc = lb + 1; lc < OHMMS_DIM; lc++) - { - ((ggg[j1])[la])(lb, lc) = sinkr * (K[ik])[la] * (K[ik])[lb] * (K[ik])[lc]; - ((ggg[j2])[la])(lb, lc) = -coskr * (K[ik])[la] * (K[ik])[lb] * (K[ik])[lc]; - ((ggg[j1])[la])(lc, lb) = ((ggg[j1])[la])(lb, lc); - ((ggg[j2])[la])(lc, lb) = ((ggg[j2])[la])(lb, lc); - ((ggg[j1])[lb])(la, lc) = ((ggg[j1])[la])(lb, lc); - ((ggg[j2])[lb])(la, lc) = ((ggg[j2])[la])(lb, lc); - ((ggg[j1])[lb])(lc, la) = ((ggg[j1])[la])(lb, lc); - ((ggg[j2])[lb])(lc, la) = ((ggg[j2])[la])(lb, lc); - ((ggg[j1])[lc])(la, lb) = ((ggg[j1])[la])(lb, lc); - ((ggg[j2])[lc])(la, lb) = ((ggg[j2])[la])(lb, lc); - ((ggg[j1])[lc])(lb, la) = ((ggg[j1])[la])(lb, lc); - ((ggg[j2])[lc])(lb, la) = ((ggg[j2])[la])(lb, lc); - } - } - } - //#if OHMMS_DIM ==3 - // for(int la=0; la<3; la++) { - // for(int lb=0; lb<3; lb++) { - // for(int lc=0; lc<3; lc++) { - // ( (ggg[j1])[la] )(lb,lc) = sinkr*(K[ik])[la]*(K[ik])[lb]*(K[ik])[lc]; - // ( (ggg[j2])[la] )(lb,lc) = -coskr*(K[ik])[la]*(K[ik])[lb]*(K[ik])[lc]; - // } - // } - // } - //#elif OHMMS_DIM ==2 - // for(int la=0; la<2; la++) { - // for(int lb=0; lb<2; lb++) { - // for(int lc=0; lc<2; lc++) { - // ( (ggg[j1])[la] )(lb,lc) = sinkr*(K[ik])[la]*(K[ik])[lb]*(K[ik])[lc]; - // ( (ggg[j2])[la] )(lb,lc) = -coskr*(K[ik])[la]*(K[ik])[lb]*(K[ik])[lc]; - // } - // } - // } - //#endif - } - } - } -}; - -/** OrbitalBuilder for Slater determinants of electron-gas -*/ -class ElectronGasOrbitalBuilder : public WaveFunctionComponentBuilder -{ -public: - ///constructor - ElectronGasOrbitalBuilder(Communicate* comm, ParticleSet& els); - - ///implement vritual function - std::unique_ptr buildComponent(xmlNodePtr cur) override; - - bool UseBackflow; -}; - -/** OrbitalBuilder for Slater determinants of electron-gas -*/ -class ElectronGasSPOBuilder : public SPOSetBuilder -{ -protected: - HEGGrid egGrid; - xmlNodePtr spo_node; - -public: - ///constructor - ElectronGasSPOBuilder(ParticleSet& p, Communicate* comm, xmlNodePtr cur); - - /** initialize the Antisymmetric wave function for electrons - *@param cur the current xml node - */ - std::unique_ptr createSPOSetFromXML(xmlNodePtr cur) override; -}; - -} // namespace qmcplusplus -#endif /** QMC_COMPLEX **/ -#endif diff --git a/src/QMCWaveFunctions/SPOSetBuilderFactory.cpp b/src/QMCWaveFunctions/SPOSetBuilderFactory.cpp index 20e7bc816c..c8fbed5e75 100644 --- a/src/QMCWaveFunctions/SPOSetBuilderFactory.cpp +++ b/src/QMCWaveFunctions/SPOSetBuilderFactory.cpp @@ -18,7 +18,6 @@ #include "SPOSetBuilderFactory.h" #include "QMCWaveFunctions/SPOSetScanner.h" -#include "QMCWaveFunctions/ElectronGas/ElectronGasOrbitalBuilder.h" #include "QMCWaveFunctions/HarmonicOscillator/SHOSetBuilder.h" #include "ModernStringUtils.hpp" #include "QMCWaveFunctions/ElectronGas/FreeParticleBuilder.h" diff --git a/src/QMCWaveFunctions/WaveFunctionFactory.cpp b/src/QMCWaveFunctions/WaveFunctionFactory.cpp index 8af5d8144a..1cddae31af 100644 --- a/src/QMCWaveFunctions/WaveFunctionFactory.cpp +++ b/src/QMCWaveFunctions/WaveFunctionFactory.cpp @@ -24,13 +24,6 @@ #include "QMCWaveFunctions/Fermion/SlaterDetBuilder.h" #include "QMCWaveFunctions/LatticeGaussianProductBuilder.h" #include "QMCWaveFunctions/ExampleHeBuilder.h" - -#if defined(QMC_COMPLEX) -#include "QMCWaveFunctions/ElectronGas/ElectronGasComplexOrbitalBuilder.h" -#else -#include "QMCWaveFunctions/ElectronGas/ElectronGasOrbitalBuilder.h" -#endif - #include "QMCWaveFunctions/PlaneWave/PWOrbitalBuilder.h" #if OHMMS_DIM == 3 && !defined(QMC_COMPLEX) #include "QMCWaveFunctions/AGPDeterminantBuilder.h" @@ -189,11 +182,10 @@ bool WaveFunctionFactory::addFermionTerm(TrialWaveFunction& psi, SPOSetBuilderFa std::unique_ptr detbuilder; if (orbtype == "electron-gas") { -#if defined(QMC_COMPLEX) - detbuilder = std::make_unique(myComm, targetPtcl); -#else - detbuilder = std::make_unique(myComm, targetPtcl); -#endif + std::ostringstream msg; + msg << "electron-gas in determinantset is deprecated"; + msg << " please use \"free\" orbitals in sposet_builder" << std::endl; + throw std::runtime_error(msg.str()); } else if (orbtype == "PWBasis" || orbtype == "PW" || orbtype == "pw") { diff --git a/src/QMCWaveFunctions/tests/test_DiracDeterminant.cpp b/src/QMCWaveFunctions/tests/test_DiracDeterminant.cpp index 1d2e28cb41..f686215035 100644 --- a/src/QMCWaveFunctions/tests/test_DiracDeterminant.cpp +++ b/src/QMCWaveFunctions/tests/test_DiracDeterminant.cpp @@ -18,11 +18,7 @@ #include "QMCWaveFunctions/WaveFunctionComponent.h" #include "QMCWaveFunctions/Fermion/DiracDeterminant.h" #include "QMCWaveFunctions/tests/FakeSPO.h" - -#ifdef QMC_COMPLEX //This is for the spinor test. -#include "QMCWaveFunctions/ElectronGas/ElectronGasComplexOrbitalBuilder.h" -#endif - +#include "QMCWaveFunctions/ElectronGas/FreeParticle.h" #include "QMCWaveFunctions/SpinorSet.h" #include @@ -539,26 +535,13 @@ void test_DiracDeterminant_spinor_update(const DetMatInvertor inverter_kind) kup[1] = PosType(0.1, 0.2, 0.3); kup[2] = PosType(0.4, 0.5, 0.6); - k2up.resize(nelec); - //For some goofy reason, EGOSet needs to be initialized with: - //1.) A k-vector list (fine). - //2.) A list of -|k|^2. To save on expensive - sign multiplication apparently. - k2up[0] = -dot(kup[0], kup[0]); - k2up[1] = -dot(kup[1], kup[1]); - k2up[2] = -dot(kup[2], kup[2]); - kdn.resize(nelec); kdn[0] = PosType(0, 0, 0); kdn[1] = PosType(-0.1, 0.2, -0.3); kdn[2] = PosType(0.4, -0.5, 0.6); - k2dn.resize(nelec); - k2dn[0] = -dot(kdn[0], kdn[0]); - k2dn[1] = -dot(kdn[1], kdn[1]); - k2dn[2] = -dot(kdn[2], kdn[2]); - - auto spo_up = std::make_unique(kup, k2up); - auto spo_dn = std::make_unique(kdn, k2dn); + auto spo_up = std::make_unique(kup); + auto spo_dn = std::make_unique(kdn); auto spinor_set = std::make_unique(); spinor_set->set_spos(std::move(spo_up), std::move(spo_dn)); diff --git a/tests/heg/heg_14_gamma/heg-long-HF.xml b/tests/heg/heg_14_gamma/heg-long-HF.xml index 8347d2090e..b9dc3655a7 100644 --- a/tests/heg/heg_14_gamma/heg-long-HF.xml +++ b/tests/heg/heg_14_gamma/heg-long-HF.xml @@ -6,7 +6,7 @@ Sample qmc run for Hartree-Fock HEG. - + 200 diff --git a/tests/heg/heg_14_gamma/heg-short-HF.xml b/tests/heg/heg_14_gamma/heg-short-HF.xml index 0da42c68fb..5e6112045a 100644 --- a/tests/heg/heg_14_gamma/heg-short-HF.xml +++ b/tests/heg/heg_14_gamma/heg-short-HF.xml @@ -6,7 +6,7 @@ Sample qmc run for Hartree-Fock HEG. - + 200 diff --git a/tests/heg/heg_14_gamma/heg.hf.wfs.xml b/tests/heg/heg_14_gamma/heg.hf.wfs.xml deleted file mode 100644 index d5754ecb9f..0000000000 --- a/tests/heg/heg_14_gamma/heg.hf.wfs.xml +++ /dev/null @@ -1,6 +0,0 @@ - - - - - - diff --git a/tests/heg/heg_54_J2rpa/det_heg.rpa.rs5.xml b/tests/heg/heg_54_J2rpa/det_heg.rpa.rs5.xml old mode 100755 new mode 100644 index 4e3115727f..8f4c4ab469 --- a/tests/heg/heg_54_J2rpa/det_heg.rpa.rs5.xml +++ b/tests/heg/heg_54_J2rpa/det_heg.rpa.rs5.xml @@ -23,7 +23,15 @@ xsi:noNamespaceSchemaLocation="http://www.mcc.uiuc.edu/qmc/schema/molecu.xsd"> - + + + + + + + + + diff --git a/tests/heg/heg_54_J2rpa/heg.rpa.rs5.xml b/tests/heg/heg_54_J2rpa/heg.rpa.rs5.xml old mode 100755 new mode 100644 From b137c6bab168d5d65b2cacf1d957bf1930ae3292 Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Tue, 28 Jun 2022 17:53:07 -0400 Subject: [PATCH 05/34] occupy half sphere of kvecs in real code --- .../ElectronGas/FreeParticleBuilder.cpp | 50 +++++++++++++++---- .../ElectronGas/FreeParticleBuilder.h | 1 + 2 files changed, 42 insertions(+), 9 deletions(-) diff --git a/src/QMCWaveFunctions/ElectronGas/FreeParticleBuilder.cpp b/src/QMCWaveFunctions/ElectronGas/FreeParticleBuilder.cpp index b162ce516c..338378e878 100644 --- a/src/QMCWaveFunctions/ElectronGas/FreeParticleBuilder.cpp +++ b/src/QMCWaveFunctions/ElectronGas/FreeParticleBuilder.cpp @@ -1,7 +1,8 @@ -#include "FreeParticleBuilder.h" #include "OhmmsData/AttributeSet.h" #include "LongRange/StructFact.h" -#include "FreeParticle.h" +#include "LongRange/KContainer.h" +#include "QMCWaveFunctions/ElectronGas/FreeParticle.h" +#include "QMCWaveFunctions/ElectronGas/FreeParticleBuilder.h" namespace qmcplusplus { @@ -20,28 +21,59 @@ std::unique_ptr FreeParticleBuilder::createSPOSetFromXML(xmlNodePtr cur) attrib.add(twist, "twist"); attrib.put(cur); + PosType tvec = targetPtcl.getLattice().k_cart(twist); #ifdef QMC_COMPLEX npw = norb; -#else - npw = std::ceil((norb+1.0)/2); -#endif targetPtcl.setTwist(twist); - PosType tvec = targetPtcl.getLattice().k_cart(twist); app_log() << "twist fraction = " << twist << std::endl; app_log() << "twist cartesian = " << tvec << std::endl; +#else + npw = std::ceil((norb+1.0)/2); + for (int ldim=0;ldim 1e-16) throw "no twist for real orbitals"; + } +#endif // extract npw k-points from container // kpts_cart is sorted by magnitude - std::vector kpts; - kpts.resize(npw); + // k=0 is not in kpts_cart + std::vector kpts(npw); kpts[0] = tvec; + const KContainer klists = targetPtcl.getSK().getKLists(); +#ifdef QMC_COMPLEX for (int ik=1;ik mkidx(npw, 0); + int ik = 1; + for (int jk=0;jk= npw) break; + } +#endif auto sposet = std::make_unique(kpts); sposet->report(" "); return sposet; } +bool FreeParticleBuilder::in_list(const int j, const std::vector l) +{ + for (int i=0;i createSPOSetFromXML(xmlNodePtr cur); private: ParticleSet& targetPtcl; + bool in_list(const int j, const std::vector l); }; } // qmcplusplus #endif From 9af84e4ca084ae6e1fbf3f200018e713a959b589 Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Tue, 28 Jun 2022 18:36:32 -0400 Subject: [PATCH 06/34] add twist to KContainer --- src/Particle/LongRange/KContainer.cpp | 8 ++++---- src/Particle/LongRange/KContainer.h | 4 ++-- src/Particle/LongRange/LRBreakup.h | 3 ++- src/Particle/SimulationCell.cpp | 3 ++- 4 files changed, 10 insertions(+), 8 deletions(-) diff --git a/src/Particle/LongRange/KContainer.cpp b/src/Particle/LongRange/KContainer.cpp index b13c5e05e5..4aa5c02aec 100644 --- a/src/Particle/LongRange/KContainer.cpp +++ b/src/Particle/LongRange/KContainer.cpp @@ -21,7 +21,7 @@ namespace qmcplusplus { -void KContainer::updateKLists(const ParticleLayout& lattice, RealType kc, unsigned ndim, bool useSphere) +void KContainer::updateKLists(const ParticleLayout& lattice, RealType kc, unsigned ndim, PosType twist, bool useSphere) { kcutoff = kc; if (kcutoff <= 0.0) @@ -29,7 +29,7 @@ void KContainer::updateKLists(const ParticleLayout& lattice, RealType kc, unsign APP_ABORT(" Illegal cutoff for KContainer"); } findApproxMMax(lattice, ndim); - BuildKLists(lattice, useSphere); + BuildKLists(lattice, twist, useSphere); app_log() << " KContainer initialised with cutoff " << kcutoff << std::endl; app_log() << " # of K-shell = " << kshell.size() << std::endl; @@ -98,7 +98,7 @@ void KContainer::findApproxMMax(const ParticleLayout& lattice, unsigned ndim) mmax[1] = 0; } -void KContainer::BuildKLists(const ParticleLayout& lattice, bool useSphere) +void KContainer::BuildKLists(const ParticleLayout& lattice, PosType twist, bool useSphere) { TinyVector TempActualMax; TinyVector kvec; @@ -125,7 +125,7 @@ void KContainer::BuildKLists(const ParticleLayout& lattice, bool useSphere) if (i == 0 && j == 0 && k == 0) continue; //Convert kvec to Cartesian - kvec_cart = lattice.k_cart(kvec); + kvec_cart = lattice.k_cart(kvec+twist); //Find modk modk2 = dot(kvec_cart, kvec_cart); if (modk2 > kcut2) diff --git a/src/Particle/LongRange/KContainer.h b/src/Particle/LongRange/KContainer.h index 2496ccc0f2..9ad6ac19a8 100644 --- a/src/Particle/LongRange/KContainer.h +++ b/src/Particle/LongRange/KContainer.h @@ -72,7 +72,7 @@ class KContainer : public QMCTraits * @param kc cutoff radius in the K * @param useSphere if true, use the |K| */ - void updateKLists(const ParticleLayout& lattice, RealType kc, unsigned ndim, bool useSphere = true); + void updateKLists(const ParticleLayout& lattice, RealType kc, unsigned ndim, PosType twist, bool useSphere = true); const auto& get_kpts_cart_soa() const { return kpts_cart_soa_; } private: @@ -81,7 +81,7 @@ class KContainer : public QMCTraits */ void findApproxMMax(const ParticleLayout& lattice, unsigned ndim); /** construct the container for k-vectors */ - void BuildKLists(const ParticleLayout& lattice, bool useSphere); + void BuildKLists(const ParticleLayout& lattice, PosType twist, bool useSphere); /** K-vector in Cartesian coordinates in SoA layout */ diff --git a/src/Particle/LongRange/LRBreakup.h b/src/Particle/LongRange/LRBreakup.h index 38233d0b18..2f49c82a23 100644 --- a/src/Particle/LongRange/LRBreakup.h +++ b/src/Particle/LongRange/LRBreakup.h @@ -334,7 +334,8 @@ int LRBreakup::SetupKVecs(mRealType kc, mRealType kcont, mRealType { //Add low |k| ( < kcont) k-points with exact degeneracy KContainer kexact; - kexact.updateKLists(Basis.get_Lattice(), kcont, Basis.get_Lattice().ndim); + PosType twist(0); + kexact.updateKLists(Basis.get_Lattice(), kcont, Basis.get_Lattice().ndim, twist); bool findK = true; mRealType kc2 = kc * kc; //use at least one shell diff --git a/src/Particle/SimulationCell.cpp b/src/Particle/SimulationCell.cpp index cbd9b5f053..b90805b61f 100644 --- a/src/Particle/SimulationCell.cpp +++ b/src/Particle/SimulationCell.cpp @@ -59,7 +59,8 @@ void SimulationCell::resetLRBox() app_log() << "--------------------------------------- " << std::endl; } - k_lists_.updateKLists(LRBox_, LRBox_.LR_kc, LRBox_.ndim); + QMCTraits::PosType twist(0); + k_lists_.updateKLists(LRBox_, LRBox_.LR_kc, LRBox_.ndim, twist); } } } From 406b5f14e1ba51d9b61e1722f4446274cd73d3c9 Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Tue, 28 Jun 2022 19:08:37 -0400 Subject: [PATCH 07/34] use twist in free-particle orbitals --- .../ElectronGas/FreeParticleBuilder.cpp | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/src/QMCWaveFunctions/ElectronGas/FreeParticleBuilder.cpp b/src/QMCWaveFunctions/ElectronGas/FreeParticleBuilder.cpp index 338378e878..5d5c3ca747 100644 --- a/src/QMCWaveFunctions/ElectronGas/FreeParticleBuilder.cpp +++ b/src/QMCWaveFunctions/ElectronGas/FreeParticleBuilder.cpp @@ -21,7 +21,9 @@ std::unique_ptr FreeParticleBuilder::createSPOSetFromXML(xmlNodePtr cur) attrib.add(twist, "twist"); attrib.put(cur); - PosType tvec = targetPtcl.getLattice().k_cart(twist); + auto lattice = targetPtcl.getLattice(); + + PosType tvec = lattice.k_cart(twist); #ifdef QMC_COMPLEX npw = norb; targetPtcl.setTwist(twist); @@ -31,20 +33,24 @@ std::unique_ptr FreeParticleBuilder::createSPOSetFromXML(xmlNodePtr cur) npw = std::ceil((norb+1.0)/2); for (int ldim=0;ldim 1e-16) throw "no twist for real orbitals"; + 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 - // k=0 is not in kpts_cart 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; - const KContainer klists = targetPtcl.getSK().getKLists(); #ifdef QMC_COMPLEX for (int ik=1;ik Date: Tue, 28 Jun 2022 19:25:29 -0400 Subject: [PATCH 08/34] cut down HEGGrid --- src/QMCWaveFunctions/ElectronGas/HEGGrid.h | 199 --------------------- 1 file changed, 199 deletions(-) diff --git a/src/QMCWaveFunctions/ElectronGas/HEGGrid.h b/src/QMCWaveFunctions/ElectronGas/HEGGrid.h index a4657e7d60..7414a10c98 100644 --- a/src/QMCWaveFunctions/ElectronGas/HEGGrid.h +++ b/src/QMCWaveFunctions/ElectronGas/HEGGrid.h @@ -69,20 +69,12 @@ struct HEGGrid HEGGrid(const PL_t& lat) : Lattice(lat) {} - - ~HEGGrid() = default; /** return the estimated number of grid in each direction */ inline int getNC(int nup) const { return static_cast(std::pow(static_cast(nup), 1.0 / 3.0)) / 2 + 1; } - /** return the estimated number of grid in each direction (upper bound) */ - inline int get_nc(int nstates) const - { - return static_cast(std::pow(static_cast(nstates), 1.0 / 3.0) * .7) + 1; - } - //return the number of k-points upto nsh-shell inline int getNumberOfKpoints(int nsh) const { @@ -92,15 +84,6 @@ struct HEGGrid return -1; } - - inline int getShellFromStates(int nst) - { - for (int i = 0; i < n_within_shell.size(); i++) - if (n_within_shell[i] == nst) - return i; - return -1; - } - //return the shell index for nkpt k-points inline int getShellIndex(int nkpt) const { @@ -117,188 +100,6 @@ struct HEGGrid */ inline T getCellLength(int nptcl, T rs_in) const { return std::pow(4.0 * M_PI * nptcl / 3.0, 1.0 / 3.0) * rs_in; } - void sortGrid(int nc) - { - int first_ix2, first_ix3; - for (int ix1 = 0; ix1 <= nc; ix1++) - { - if (ix1 == 0) - first_ix2 = 0; - else - first_ix2 = -nc; - for (int ix2 = first_ix2; ix2 <= nc; ix2++) - { - if (ix1 == 0 && ix2 == 0) - first_ix3 = 1; - else - first_ix3 = -nc; - for (int ix3 = first_ix3; ix3 <= nc; ix3++) - { - int ih = ix1 * ix1 + ix2 * ix2 + ix3 * ix3; - if (auto it = rs.find(ih); it == rs.end()) - rs[ih] = {PosType(ix1, ix2, ix3)}; - else - it->second.push_back(PosType(ix1, ix2, ix3)); - } - } - } - } - - void createGrid(int nc, int nkpts) - { - if (rs.empty()) - sortGrid(nc); - NumKptsHalf = nkpts; - kpt.resize(nkpts); - mk2.resize(nkpts); - int ikpt = 0; - MaxKsq = 0.0; - auto rs_it = rs.begin(); - auto rs_end = rs.end(); - while (ikpt < nkpts && rs_it != rs_end) - { - auto ns_it = rs_it->second.begin(); - auto ns_end = rs_it->second.end(); - T minus_ksq = -Lattice.ksq(*ns_it); - while (ikpt < nkpts && ns_it != ns_end) - { - kpt[ikpt] = Lattice.k_cart(*ns_it); - mk2[ikpt] = minus_ksq; - ++ikpt; - ++ns_it; - } - ++rs_it; - } - MaxKsq = Lattice.ksq(rs_it->second.front()); - app_log() << "List of kpoints (half-sphere) " << std::endl; - for (int ik = 0; ik < kpt.size(); ik++) - { - app_log() << ik << " " << kpt[ik] << " " << mk2[ik] << std::endl; - } - } - - - void clear_kpoints() { kpoints_grid.reset(); } - - - void create_kpoints(int nc, const PosType& tw, T tol = 1e-6) - { - if (!kpoints_grid.has_value()) - kpoints_grid = kpoints_t(); - else if (nc <= nctmp) - return; - nctmp = nc; - kpoints_t& kpoints = *kpoints_grid; - app_log() << " resizing kpoint grid" << std::endl; - app_log() << " current size = " << kpoints.size() << std::endl; - // make space for the kpoint grid - int nkpoints = pow(2 * (nc + 1) + 1, 3); - kpoints.resize(nkpoints); - app_log() << " cubic size = " << kpoints.size() << std::endl; - typename kpoints_t::iterator kptmp, kp = kpoints.begin(), kp_end = kpoints.end(); - // make the kpoint grid - T k2max = std::numeric_limits::max(); - for (int i0 = -nc - 1; i0 <= nc + 1; ++i0) - for (int i1 = -nc - 1; i1 <= nc + 1; ++i1) - for (int i2 = -nc - 1; i2 <= nc + 1; ++i2) - { - PosType k(i0 + tw[0], i1 + tw[1], i2 + tw[2]); - kp->k = Lattice.k_cart(k); - kp->k2 = Lattice.ksq(k); - if (std::abs(i0) == (nc + 1) || std::abs(i1) == (nc + 1) || std::abs(i2) == (nc + 1)) - k2max = std::min(k2max, kp->k2); - ++kp; - } - // sort kpoints by magnitude - sort(kpoints.begin(), kpoints.end(), kpdata_comp); - // eliminate kpoints outside of inscribing sphere - int nkp = 0; - kp = kpoints.begin(); - while (kp != kp_end && kp->k2 < k2max + 1e-3) - { - nkp++; - ++kp; - } - kpoints.resize(nkp); - app_log() << " new spherical size = " << kpoints.size() << std::endl; - kp_end = kpoints.end(); - // count degeneracies - kp = kpoints.begin(); - while (kp != kp_end) - { - T k2 = kp->k2; - kptmp = kp; - int g = 1; - ++kptmp; - // look ahead to count - while (kptmp != kp_end && std::abs(kptmp->k2 - k2) < tol) - { - g++; - ++kptmp; - } - kp->g = g; - // run over degenerate states to assign - for (int n = 0; n < g - 1; ++n) - (++kp)->g = g; - ++kp; - } - //app_log()<<"create_kpoints"<< std::endl; - //app_log()<<" nkpoints = "<k2<<" "<g<<" "<k<< std::endl; - //APP_ABORT("end create_kpoints"); - } - - - void createGrid(int nc, int nkpts, const PosType& twistAngle) - { - twist = twistAngle; - create_kpoints(nc, twistAngle); - kpoints_t& kpoints = *kpoints_grid; - if (nkpts > kpoints.size()) - APP_ABORT("HEGGrid::createGrid requested more kpoints than created"); - kpt.resize(nkpts); - mk2.resize(nkpts); - deg.resize(nkpts); - for (int i = 0; i < nkpts; ++i) - { - const kpdata_t& kp = kpoints[i]; - kpt[i] = kp.k; - mk2[i] = -kp.k2; - deg[i] = kp.g; - } - app_log() << "List of kpoints with twist = " << twistAngle << std::endl; - for (int ik = 0; ik < kpt.size(); ik++) - app_log() << ik << " " << kpt[ik] << " " << -mk2[ik] << std::endl; - } - - - void createGrid(const std::vector& states, T tol = 1e-6) { createGrid(states, twist, tol); } - - - void createGrid(const std::vector& states, const PosType& twistAngle, T tol = 1e-6) - { - int smax = 0; - for (int i = 0; i < states.size(); ++i) - smax = std::max(smax, states[i]); - smax++; - create_kpoints(get_nc(smax), twistAngle, tol); - kpoints_t& kpoints = *kpoints_grid; - if (smax > kpoints.size()) - APP_ABORT("HEGGrid::createGrid(states) requested more kpoints than created"); - int nkpts = states.size(); - kpt.resize(nkpts); - mk2.resize(nkpts); - deg.resize(nkpts); - for (int i = 0; i < states.size(); ++i) - { - const kpdata_t& kp = kpoints[states[i]]; - kpt[i] = kp.k; - mk2[i] = -kp.k2; - deg[i] = kp.g; - } - } }; } // namespace qmcplusplus From a6c02e04ad87fa48575ee132ea53e00cc66f8dbf Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Tue, 28 Jun 2022 19:28:09 -0400 Subject: [PATCH 09/34] cut HEG some more --- src/QMCWaveFunctions/ElectronGas/HEGGrid.h | 37 +--------------------- 1 file changed, 1 insertion(+), 36 deletions(-) diff --git a/src/QMCWaveFunctions/ElectronGas/HEGGrid.h b/src/QMCWaveFunctions/ElectronGas/HEGGrid.h index 7414a10c98..3c8659d036 100644 --- a/src/QMCWaveFunctions/ElectronGas/HEGGrid.h +++ b/src/QMCWaveFunctions/ElectronGas/HEGGrid.h @@ -8,6 +8,7 @@ // Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign // Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory // Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory +// Yubo "Paul" Yang, yubo.paul.yang@gmail.com, CCQ @ Flatiron // // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign ////////////////////////////////////////////////////////////////////////////////////// @@ -17,56 +18,20 @@ #define QMCPLUSPLUS_HEGGRID_H #include "Lattice/CrystalLattice.h" -#include -#include -#include namespace qmcplusplus { -template -struct kpdata -{ - TinyVector k; - T k2; - int g; -}; - - -template -bool kpdata_comp(const kpdata& left, const kpdata& right) -{ - return left.k2 < right.k2; -} //three-d specialization template struct HEGGrid { using PL_t = CrystalLattice; - using PosType = typename PL_t::SingleParticlePos; - using RealType = typename PL_t::Scalar_t; - ///number of kpoints of a half sphere excluding gamma - int NumKptsHalf; - ///maxmim ksq - T MaxKsq; const PL_t& Lattice; - std::map> rs; - std::vector kpt; - std::vector mk2; - std::vector deg; static constexpr std::array n_within_shell{{1, 7, 19, 27, 33, 57, 81, 93, 123, 147, 171, 179, 203, 251, 257, 305, 341, 365, 389, 437, 461, 485, 515, 587, 619, 691, 739, 751, 799, 847, 895}}; - PosType twist{0.0}; - - - using kpdata_t = kpdata; - using kpoints_t = std::vector; - - std::optional kpoints_grid; - int nctmp{-1}; - HEGGrid(const PL_t& lat) : Lattice(lat) {} ~HEGGrid() = default; From c6d95e61b4e8da2b703fb5f46a2cc5ba06e01fef Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Tue, 28 Jun 2022 19:58:01 -0400 Subject: [PATCH 10/34] improve backflow tests --- tests/heg/heg_14_gamma/CMakeLists.txt | 2 +- tests/heg/heg_14_gamma/heg-long-SJB.xml | 4 ++-- tests/heg/heg_14_gamma/heg-short-SJB.xml | 4 ++-- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/heg/heg_14_gamma/CMakeLists.txt b/tests/heg/heg_14_gamma/CMakeLists.txt index a87997fac3..9d091b1db3 100644 --- a/tests/heg/heg_14_gamma/CMakeLists.txt +++ b/tests/heg/heg_14_gamma/CMakeLists.txt @@ -231,7 +231,7 @@ if(NOT QMC_CUDA) HEG14GLSJ_DMC_SCALARS # DMC ) - if(NOT QMC_COMPLEX) + if(QMC_COMPLEX) # # HEG14G - Slater-Jastrow-Backflow VMC # Run 16x1, 4x4, 1x16 to test SJB cloning diff --git a/tests/heg/heg_14_gamma/heg-long-SJB.xml b/tests/heg/heg_14_gamma/heg-long-SJB.xml index 5a1a10af54..c7434c8797 100644 --- a/tests/heg/heg_14_gamma/heg-long-SJB.xml +++ b/tests/heg/heg_14_gamma/heg-long-SJB.xml @@ -8,10 +8,10 @@ Sample qmc run for Slater-Jastrow-Backflow HEG. - + 200 800 400 -5.0 +1.0 diff --git a/tests/heg/heg_14_gamma/heg-short-SJB.xml b/tests/heg/heg_14_gamma/heg-short-SJB.xml index 74c1a440e9..b230de32db 100644 --- a/tests/heg/heg_14_gamma/heg-short-SJB.xml +++ b/tests/heg/heg_14_gamma/heg-short-SJB.xml @@ -8,10 +8,10 @@ Sample qmc run for Slater-Jastrow-Backflow HEG. - + 200 800 40 -5.0 +1.0 From fa6aacdf43c6b3e437e3236db7622969b50f404d Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Tue, 28 Jun 2022 20:01:47 -0400 Subject: [PATCH 11/34] mark makeClone as override --- src/QMCWaveFunctions/ElectronGas/FreeParticle.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/QMCWaveFunctions/ElectronGas/FreeParticle.h b/src/QMCWaveFunctions/ElectronGas/FreeParticle.h index 27b2b6345a..a8bbd5f3bb 100644 --- a/src/QMCWaveFunctions/ElectronGas/FreeParticle.h +++ b/src/QMCWaveFunctions/ElectronGas/FreeParticle.h @@ -68,7 +68,7 @@ class FreeParticle : public SPOSet GGGMatrix& d3phi_mat) override; // ---- begin required overrides - std::unique_ptr makeClone() const {return std::make_unique(*this);} + std::unique_ptr makeClone() const override {return std::make_unique(*this);} void resetParameters(const opt_variables_type& optVariables) override {} //called by BFTrans} void setOrbitalSetSize(int norbs) override {APP_ABORT("not implemented")}; // required overrides end ---- From d83e24ecd453b30144248ab2db7e011f59a94251 Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Tue, 28 Jun 2022 20:13:58 -0400 Subject: [PATCH 12/34] document free SPOSet --- docs/intro_wavefunction.rst | 18 +++++++++++------- tests/heg/heg_14_gamma/CMakeLists.txt | 2 -- 2 files changed, 11 insertions(+), 9 deletions(-) diff --git a/docs/intro_wavefunction.rst b/docs/intro_wavefunction.rst index b9a59a3950..924088b6dd 100644 --- a/docs/intro_wavefunction.rst +++ b/docs/intro_wavefunction.rst @@ -706,14 +706,11 @@ Plane-wave basis sets Homogeneous electron gas ~~~~~~~~~~~~~~~~~~~~~~~~ -The interacting Fermi liquid has its own special ``determinantset`` for filling up a -Fermi surface. The shell number can be specified separately for both spin-up and spin-down. -This determines how many electrons to include of each time; only closed shells are currently -implemented. The shells are filled according to the rules of a square box; if other lattice -vectors are used, the electrons might not fill up a complete shell. +The interacting Fermi liquid can be created using a determinant of free-particle orbitals. +The lowest-energy plane-wave states compatible with the boundary condition are occupied. This following example can also be used for Helium simulations by specifying the -proper pair interaction in the Hamiltonian section. +proper pair interaction in the Hamiltonian section and using a bosonic wavefunction. .. code-block:: :caption: 2D Fermi liquid example: particle specification @@ -740,7 +737,14 @@ proper pair interaction in the Hamiltonian section. :name: Listing 9 - + + + + + + + + diff --git a/tests/heg/heg_14_gamma/CMakeLists.txt b/tests/heg/heg_14_gamma/CMakeLists.txt index 9d091b1db3..faa99a8549 100644 --- a/tests/heg/heg_14_gamma/CMakeLists.txt +++ b/tests/heg/heg_14_gamma/CMakeLists.txt @@ -231,7 +231,6 @@ if(NOT QMC_CUDA) HEG14GLSJ_DMC_SCALARS # DMC ) - if(QMC_COMPLEX) # # HEG14G - Slater-Jastrow-Backflow VMC # Run 16x1, 4x4, 1x16 to test SJB cloning @@ -325,7 +324,6 @@ if(NOT QMC_CUDA) list(APPEND DET_HEG14GSSJB_SCALARS "totenergy" "-0.92845399 0.000001") list(APPEND DET_HEG14GSSJB_SCALARS "kinetic" "0.76941751 0.000001") list(APPEND DET_HEG14GSSJB_SCALARS "potential" "-1.69787150 0.000001") - endif() qmc_run_and_check( deterministic-heg_14_gamma-sjb From a7f2c180a465e91300a64786eb5e55209934a1ee Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Tue, 28 Jun 2022 20:15:06 -0400 Subject: [PATCH 13/34] override createSPOSetFromXML --- src/QMCWaveFunctions/ElectronGas/FreeParticleBuilder.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/QMCWaveFunctions/ElectronGas/FreeParticleBuilder.h b/src/QMCWaveFunctions/ElectronGas/FreeParticleBuilder.h index 5d5bc834c3..a3db8ce3f5 100644 --- a/src/QMCWaveFunctions/ElectronGas/FreeParticleBuilder.h +++ b/src/QMCWaveFunctions/ElectronGas/FreeParticleBuilder.h @@ -11,7 +11,7 @@ class FreeParticleBuilder : public SPOSetBuilder FreeParticleBuilder(ParticleSet& els, Communicate* comm, xmlNodePtr cur); ~FreeParticleBuilder(){} - std::unique_ptr createSPOSetFromXML(xmlNodePtr cur); + std::unique_ptr createSPOSetFromXML(xmlNodePtr cur) override; private: ParticleSet& targetPtcl; bool in_list(const int j, const std::vector l); From b1ec05ed17b6109b6d32aaabc91dd874ddfa97e3 Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Tue, 28 Jun 2022 20:31:39 -0400 Subject: [PATCH 14/34] fix cmake for heg test --- tests/heg/heg_14_gamma/CMakeLists.txt | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/tests/heg/heg_14_gamma/CMakeLists.txt b/tests/heg/heg_14_gamma/CMakeLists.txt index faa99a8549..5af9b64558 100644 --- a/tests/heg/heg_14_gamma/CMakeLists.txt +++ b/tests/heg/heg_14_gamma/CMakeLists.txt @@ -324,6 +324,7 @@ if(NOT QMC_CUDA) list(APPEND DET_HEG14GSSJB_SCALARS "totenergy" "-0.92845399 0.000001") list(APPEND DET_HEG14GSSJB_SCALARS "kinetic" "0.76941751 0.000001") list(APPEND DET_HEG14GSSJB_SCALARS "potential" "-1.69787150 0.000001") + endif() qmc_run_and_check( deterministic-heg_14_gamma-sjb @@ -367,10 +368,6 @@ if(NOT QMC_CUDA) ) endif() - - else() - message(VERBOSE "Skipping HEG backflow tests because they are not supported by complex build (QMC_COMPLEX=1)") - endif() else() message(VERBOSE "Skipping all HEG tests because they are not supported by CUDA build (QMC_CUDA=1)") endif() From 53318ece9df5968aa378901c4a534bb40cd206eb Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Tue, 28 Jun 2022 20:36:05 -0400 Subject: [PATCH 15/34] remove twist test in real code --- tests/heg/heg_54_J2rpa/CMakeLists.txt | 29 ++++++++++++++------------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/tests/heg/heg_54_J2rpa/CMakeLists.txt b/tests/heg/heg_54_J2rpa/CMakeLists.txt index fbd0c3db8b..111b1c4376 100644 --- a/tests/heg/heg_54_J2rpa/CMakeLists.txt +++ b/tests/heg/heg_54_J2rpa/CMakeLists.txt @@ -35,20 +35,6 @@ if(NOT QMC_CUDA) list(APPEND DET_HEG54J2RPA_SCALARS "potential" "-7.18586631 0.000001") list(APPEND DET_HEG54J2RPA_SCALARS "eeenergy" "-7.18586631 0.000001") endif() - else() - if(QMC_MIXED_PRECISION) - list(APPEND DET_HEG54J2RPA_SCALARS "totenergy" "-3.92025567 0.000005") - list(APPEND DET_HEG54J2RPA_SCALARS "kinetic" "3.63462568 0.000005") - list(APPEND DET_HEG54J2RPA_SCALARS "potential" "-7.55488135 0.000005") - list(APPEND DET_HEG54J2RPA_SCALARS "eeenergy" "-7.55488135 0.000005") - else() - list(APPEND DET_HEG54J2RPA_SCALARS "totenergy" "-3.92025567 0.000001") - list(APPEND DET_HEG54J2RPA_SCALARS "kinetic" "3.63462568 0.000001") - list(APPEND DET_HEG54J2RPA_SCALARS "potential" "-7.55488135 0.000001") - list(APPEND DET_HEG54J2RPA_SCALARS "eeenergy" "-7.55488135 0.000001") - endif() - endif() - qmc_run_and_check( deterministic-heg_54_J2rpa "${qmcpack_SOURCE_DIR}/tests/heg/heg_54_J2rpa" @@ -59,6 +45,21 @@ if(NOT QMC_CUDA) TRUE 1 DET_HEG54J2RPA_SCALARS) + else() + ## real code should not be ran with a twist + ##if(QMC_MIXED_PRECISION) + ## list(APPEND DET_HEG54J2RPA_SCALARS "totenergy" "-3.92025567 0.000005") + ## list(APPEND DET_HEG54J2RPA_SCALARS "kinetic" "3.63462568 0.000005") + ## list(APPEND DET_HEG54J2RPA_SCALARS "potential" "-7.55488135 0.000005") + ## list(APPEND DET_HEG54J2RPA_SCALARS "eeenergy" "-7.55488135 0.000005") + ##else() + ## list(APPEND DET_HEG54J2RPA_SCALARS "totenergy" "-3.92025567 0.000001") + ## list(APPEND DET_HEG54J2RPA_SCALARS "kinetic" "3.63462568 0.000001") + ## list(APPEND DET_HEG54J2RPA_SCALARS "potential" "-7.55488135 0.000001") + ## list(APPEND DET_HEG54J2RPA_SCALARS "eeenergy" "-7.55488135 0.000001") + ##endif() + endif() + else() message(VERBOSE "Skipping HEG J2 RPA test because they are not supported by CUDA build (QMC_CUDA=1)") From e78c3ce685fc26ed0e60a7a90388f0e5fd4ef59d Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Thu, 30 Jun 2022 11:19:32 -0400 Subject: [PATCH 16/34] add batch driver tests --- tests/heg/heg_14_gamma/CMakeLists.txt | 37 +++++++++++ .../heg_14_gamma/batch/deter-heg-SJ-batch.xml | 60 +++++++++++++++++ .../batch/deter-heg-SJB-batch.xml | 66 +++++++++++++++++++ .../batch/short-heg-SJ-batched.xml | 58 ++++++++++++++++ 4 files changed, 221 insertions(+) create mode 100644 tests/heg/heg_14_gamma/batch/deter-heg-SJ-batch.xml create mode 100644 tests/heg/heg_14_gamma/batch/deter-heg-SJB-batch.xml create mode 100644 tests/heg/heg_14_gamma/batch/short-heg-SJ-batched.xml diff --git a/tests/heg/heg_14_gamma/CMakeLists.txt b/tests/heg/heg_14_gamma/CMakeLists.txt index 5af9b64558..daf076fe43 100644 --- a/tests/heg/heg_14_gamma/CMakeLists.txt +++ b/tests/heg/heg_14_gamma/CMakeLists.txt @@ -89,6 +89,7 @@ if(NOT QMC_CUDA) list(APPEND HEG14GSSJ_NEW_SCALARS "totenergy" "-1.073286 0.000271") # total energy list(APPEND HEG14GSSJ_NEW_SCALARS "variance" " 0.024564 0.000342") # energy variance + # legacy driver qmc_run_and_check( short-heg_14_gamma-sj_new "${qmcpack_SOURCE_DIR}/tests/heg/heg_14_gamma" @@ -101,6 +102,42 @@ if(NOT QMC_CUDA) HEG14GSSJ_NEW_SCALARS # VMC ) + # batch driver + qmc_run_and_check( + short-heg_14_gamma-sj_batch + "${qmcpack_SOURCE_DIR}/tests/heg/heg_14_gamma/batch" + short_heg_SJ_batch + short-heg-SJ-batch.xml + 1 + 16 + TRUE + 0 + HEG14GSSJ_NEW_SCALARS # VMC + ) + + # batch deterministic + if(QMC_MIXED_PRECISION) + list(APPEND DET_HEG14GSSJ_SCALARS "totenergy" "-1.24699388 0.00002") + list(APPEND DET_HEG14GSSJ_SCALARS "kinetic" "0.77668380 0.000005") + list(APPEND DET_HEG14GSSJ_SCALARS "potential" "-2.02367769 0.00002") + else() + list(APPEND DET_HEG14GSSJ_SCALARS "totenergy" "-1.24699388 0.000001") + list(APPEND DET_HEG14GSSJ_SCALARS "kinetic" "0.77668380 0.000001") + list(APPEND DET_HEG14GSSJ_SCALARS "potential" "-2.02367769 0.000001") + endif() + + qmc_run_and_check( + deterministic-heg_14_gamma-sj-batch + "${qmcpack_SOURCE_DIR}/tests/heg/heg_14_gamma/batch" + deter_heg_SJ_batch + deter-heg-SJ-batch.xml + 1 + 1 + TRUE + 0 + DET_HEG14GSSJ_SCALARS # VMC + ) + # # HEG14G - Slater-Jastrow DMC # diff --git a/tests/heg/heg_14_gamma/batch/deter-heg-SJ-batch.xml b/tests/heg/heg_14_gamma/batch/deter-heg-SJ-batch.xml new file mode 100644 index 0000000000..9d36a6d9d9 --- /dev/null +++ b/tests/heg/heg_14_gamma/batch/deter-heg-SJ-batch.xml @@ -0,0 +1,60 @@ + + + + batch + + + + + + p p p + + 6 + 5.0 + 14 + + + + -1 + 1.0 + + + -1 + 1.0 + + + + + + + + + + + + + + + + 1.082858193 0.6653279375 0.4358910287 0.2243616172 0.1102948764 + + + + + 1.696171854 1.047722154 0.6275148566 0.3175982878 0.1446706214 + + + + + + + + + + + 2 + 3 + 3 + 5.0 + + diff --git a/tests/heg/heg_14_gamma/batch/deter-heg-SJB-batch.xml b/tests/heg/heg_14_gamma/batch/deter-heg-SJB-batch.xml new file mode 100644 index 0000000000..312aa59422 --- /dev/null +++ b/tests/heg/heg_14_gamma/batch/deter-heg-SJB-batch.xml @@ -0,0 +1,66 @@ + + + + batch + + + + + + p p p + + 6 + 5.0 + 14 + + + + -1 + 1.0 + + + -1 + 1.0 + + + + + + + + + + + + + + + 0.01062989927 0.06173674368 0.08246051331 0.01993157184 0.02842145713 + + + 0.5824836964 0.272210709 0.1644887754 0.07508022611 0.03887818308 + + + + + + + 1.070711605 0.8007940734 0.5536996638 0.2534227771 0.1367602245 + + + 1.841049409 1.15635041 0.6935822234 0.3482003193 0.1525416364 + + + + + + + + + + 2 + 2 + 3 + 5.0 + + diff --git a/tests/heg/heg_14_gamma/batch/short-heg-SJ-batched.xml b/tests/heg/heg_14_gamma/batch/short-heg-SJ-batched.xml new file mode 100644 index 0000000000..7ff31840f6 --- /dev/null +++ b/tests/heg/heg_14_gamma/batch/short-heg-SJ-batched.xml @@ -0,0 +1,58 @@ + + + + batch + + + + + p p p + + 6 + 5.0 + 14 + + + + -1 + 1.0 + + + -1 + 1.0 + + + + + + + + + + + + + + + + 1.082858193 0.6653279375 0.4358910287 0.2243616172 0.1102948764 + + + + + 1.696171854 1.047722154 0.6275148566 0.3175982878 0.1446706214 + + + + + + + + + + 200 + 800 + 40 + 5.0 + + From 5ad23b157bc08a5fdf44b9ae247301527329eb0d Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Thu, 30 Jun 2022 14:52:52 -0400 Subject: [PATCH 17/34] rename FreeParticle to FreeOrbital --- src/QMCHamiltonians/tests/test_ecp.cpp | 6 ++--- src/QMCWaveFunctions/CMakeLists.txt | 2 +- .../{FreeParticle.cpp => FreeOrbital.cpp} | 22 +++++++++---------- .../{FreeParticle.h => FreeOrbital.h} | 12 +++++----- ...icleBuilder.cpp => FreeOrbitalBuilder.cpp} | 12 +++++----- ...ParticleBuilder.h => FreeOrbitalBuilder.h} | 10 ++++----- src/QMCWaveFunctions/SPOSetBuilderFactory.cpp | 4 ++-- .../tests/test_DiracDeterminant.cpp | 6 ++--- 8 files changed, 37 insertions(+), 37 deletions(-) rename src/QMCWaveFunctions/ElectronGas/{FreeParticle.cpp => FreeOrbital.cpp} (94%) rename src/QMCWaveFunctions/ElectronGas/{FreeParticle.h => FreeOrbital.h} (92%) rename src/QMCWaveFunctions/ElectronGas/{FreeParticleBuilder.cpp => FreeOrbitalBuilder.cpp} (80%) rename src/QMCWaveFunctions/ElectronGas/{FreeParticleBuilder.h => FreeOrbitalBuilder.h} (53%) diff --git a/src/QMCHamiltonians/tests/test_ecp.cpp b/src/QMCHamiltonians/tests/test_ecp.cpp index 569afda532..0432077194 100644 --- a/src/QMCHamiltonians/tests/test_ecp.cpp +++ b/src/QMCHamiltonians/tests/test_ecp.cpp @@ -36,7 +36,7 @@ #include "LongRange/EwaldHandler3D.h" //This is for the spinor test. -#include "QMCWaveFunctions/ElectronGas/FreeParticle.h" +#include "QMCWaveFunctions/ElectronGas/FreeOrbital.h" namespace qmcplusplus { @@ -503,8 +503,8 @@ TEST_CASE("Evaluate_soecp", "[hamiltonian]") kdn.resize(nelec); kdn[0] = PosType(2, 2, 2); - auto spo_up = std::make_unique(kup); - auto spo_dn = std::make_unique(kdn); + auto spo_up = std::make_unique(kup); + auto spo_dn = std::make_unique(kdn); auto spinor_set = std::make_unique(); spinor_set->set_spos(std::move(spo_up), std::move(spo_dn)); diff --git a/src/QMCWaveFunctions/CMakeLists.txt b/src/QMCWaveFunctions/CMakeLists.txt index 8024e8aea0..3bf398e335 100644 --- a/src/QMCWaveFunctions/CMakeLists.txt +++ b/src/QMCWaveFunctions/CMakeLists.txt @@ -63,7 +63,7 @@ set(JASTROW_SRCS set(JASTROW_OMPTARGET_SRCS Jastrow/J2OMPTarget.cpp Jastrow/BsplineFunctor.cpp) -set(FERMION_SRCS ${FERMION_SRCS} ElectronGas/FreeParticle.cpp ElectronGas/FreeParticleBuilder.cpp) +set(FERMION_SRCS ${FERMION_SRCS} ElectronGas/FreeOrbital.cpp ElectronGas/FreeOrbitalBuilder.cpp) # wavefunctions only availbale to 3-dim problems if(OHMMS_DIM MATCHES 3) diff --git a/src/QMCWaveFunctions/ElectronGas/FreeParticle.cpp b/src/QMCWaveFunctions/ElectronGas/FreeOrbital.cpp similarity index 94% rename from src/QMCWaveFunctions/ElectronGas/FreeParticle.cpp rename to src/QMCWaveFunctions/ElectronGas/FreeOrbital.cpp index 71bd668dd6..e06ca0eea2 100644 --- a/src/QMCWaveFunctions/ElectronGas/FreeParticle.cpp +++ b/src/QMCWaveFunctions/ElectronGas/FreeOrbital.cpp @@ -1,9 +1,9 @@ -#include "FreeParticle.h" +#include "FreeOrbital.h" namespace qmcplusplus { -FreeParticle::FreeParticle(const std::vector& kpts_cart) +FreeOrbital::FreeOrbital(const std::vector& kpts_cart) : kvecs(kpts_cart), #ifdef QMC_COMPLEX mink(0), // first k at twist may not be 0 @@ -24,9 +24,9 @@ FreeParticle::FreeParticle(const std::vector& kpts_cart) } } -FreeParticle::~FreeParticle(){} +FreeOrbital::~FreeOrbital(){} -void FreeParticle::evaluateVGL( +void FreeOrbital::evaluateVGL( const ParticleSet& P, int iat, ValueVector& pvec, @@ -60,7 +60,7 @@ void FreeParticle::evaluateVGL( #endif } -void FreeParticle::evaluateValue( +void FreeOrbital::evaluateValue( const ParticleSet& P, int iat, ValueVector& pvec) @@ -84,7 +84,7 @@ void FreeParticle::evaluateValue( #endif } -void FreeParticle::evaluate_notranspose( +void FreeOrbital::evaluate_notranspose( const ParticleSet& P, int first, int last, @@ -101,7 +101,7 @@ void FreeParticle::evaluate_notranspose( } } -void FreeParticle::evaluate_notranspose( +void FreeOrbital::evaluate_notranspose( const ParticleSet& P, int first, int last, @@ -165,7 +165,7 @@ void FreeParticle::evaluate_notranspose( } } -void FreeParticle::evaluate_notranspose( +void FreeOrbital::evaluate_notranspose( const ParticleSet& P, int first, int last, @@ -264,14 +264,14 @@ void FreeParticle::evaluate_notranspose( } } -void FreeParticle::report(const std::string& pad) const +void FreeOrbital::report(const std::string& pad) const { - app_log() << pad << "FreeParticle report" << std::endl; + app_log() << pad << "FreeOrbital report" << std::endl; for (int ik=0;ik& kpts_cart); - ~FreeParticle(); + FreeOrbital(const std::vector& kpts_cart); + ~FreeOrbital(); // phi[i][j] is phi_j(r_i), i.e. electron i in orbital j // i \in [first, last) @@ -68,7 +68,7 @@ class FreeParticle : public SPOSet GGGMatrix& d3phi_mat) override; // ---- begin required overrides - std::unique_ptr makeClone() const override {return std::make_unique(*this);} + std::unique_ptr makeClone() const override {return std::make_unique(*this);} void resetParameters(const opt_variables_type& optVariables) override {} //called by BFTrans} void setOrbitalSetSize(int norbs) override {APP_ABORT("not implemented")}; // required overrides end ---- diff --git a/src/QMCWaveFunctions/ElectronGas/FreeParticleBuilder.cpp b/src/QMCWaveFunctions/ElectronGas/FreeOrbitalBuilder.cpp similarity index 80% rename from src/QMCWaveFunctions/ElectronGas/FreeParticleBuilder.cpp rename to src/QMCWaveFunctions/ElectronGas/FreeOrbitalBuilder.cpp index 5d5c3ca747..32ee034b5a 100644 --- a/src/QMCWaveFunctions/ElectronGas/FreeParticleBuilder.cpp +++ b/src/QMCWaveFunctions/ElectronGas/FreeOrbitalBuilder.cpp @@ -1,18 +1,18 @@ #include "OhmmsData/AttributeSet.h" #include "LongRange/StructFact.h" #include "LongRange/KContainer.h" -#include "QMCWaveFunctions/ElectronGas/FreeParticle.h" -#include "QMCWaveFunctions/ElectronGas/FreeParticleBuilder.h" +#include "QMCWaveFunctions/ElectronGas/FreeOrbital.h" +#include "QMCWaveFunctions/ElectronGas/FreeOrbitalBuilder.h" namespace qmcplusplus { -FreeParticleBuilder::FreeParticleBuilder(ParticleSet& els, Communicate* comm, xmlNodePtr cur) +FreeOrbitalBuilder::FreeOrbitalBuilder(ParticleSet& els, Communicate* comm, xmlNodePtr cur) : SPOSetBuilder("PW", comm), targetPtcl(els) {} -std::unique_ptr FreeParticleBuilder::createSPOSetFromXML(xmlNodePtr cur) +std::unique_ptr FreeOrbitalBuilder::createSPOSetFromXML(xmlNodePtr cur) { int npw, norb; PosType twist(0.0); @@ -68,12 +68,12 @@ std::unique_ptr FreeParticleBuilder::createSPOSetFromXML(xmlNodePtr cur) if (ik >= npw) break; } #endif - auto sposet = std::make_unique(kpts); + auto sposet = std::make_unique(kpts); sposet->report(" "); return sposet; } -bool FreeParticleBuilder::in_list(const int j, const std::vector l) +bool FreeOrbitalBuilder::in_list(const int j, const std::vector l) { for (int i=0;i createSPOSetFromXML(xmlNodePtr cur) override; private: diff --git a/src/QMCWaveFunctions/SPOSetBuilderFactory.cpp b/src/QMCWaveFunctions/SPOSetBuilderFactory.cpp index c8fbed5e75..55a16ca2be 100644 --- a/src/QMCWaveFunctions/SPOSetBuilderFactory.cpp +++ b/src/QMCWaveFunctions/SPOSetBuilderFactory.cpp @@ -20,7 +20,7 @@ #include "QMCWaveFunctions/SPOSetScanner.h" #include "QMCWaveFunctions/HarmonicOscillator/SHOSetBuilder.h" #include "ModernStringUtils.hpp" -#include "QMCWaveFunctions/ElectronGas/FreeParticleBuilder.h" +#include "QMCWaveFunctions/ElectronGas/FreeOrbitalBuilder.h" #if OHMMS_DIM == 3 #include "QMCWaveFunctions/LCAO/LCAOrbitalBuilder.h" @@ -98,7 +98,7 @@ std::unique_ptr SPOSetBuilderFactory::createSPOSetBuilder(xmlNode else if (type == "jellium" || type == "heg" || type == "free") { app_log() << "Free-particle SPO set" << std::endl; - bb = std::make_unique(targetPtcl, myComm, rootNode); + bb = std::make_unique(targetPtcl, myComm, rootNode); } else if (type == "sho") { diff --git a/src/QMCWaveFunctions/tests/test_DiracDeterminant.cpp b/src/QMCWaveFunctions/tests/test_DiracDeterminant.cpp index f686215035..1e2bef713b 100644 --- a/src/QMCWaveFunctions/tests/test_DiracDeterminant.cpp +++ b/src/QMCWaveFunctions/tests/test_DiracDeterminant.cpp @@ -18,7 +18,7 @@ #include "QMCWaveFunctions/WaveFunctionComponent.h" #include "QMCWaveFunctions/Fermion/DiracDeterminant.h" #include "QMCWaveFunctions/tests/FakeSPO.h" -#include "QMCWaveFunctions/ElectronGas/FreeParticle.h" +#include "QMCWaveFunctions/ElectronGas/FreeOrbital.h" #include "QMCWaveFunctions/SpinorSet.h" #include @@ -540,8 +540,8 @@ void test_DiracDeterminant_spinor_update(const DetMatInvertor inverter_kind) kdn[1] = PosType(-0.1, 0.2, -0.3); kdn[2] = PosType(0.4, -0.5, 0.6); - auto spo_up = std::make_unique(kup); - auto spo_dn = std::make_unique(kdn); + auto spo_up = std::make_unique(kup); + auto spo_dn = std::make_unique(kdn); auto spinor_set = std::make_unique(); spinor_set->set_spos(std::move(spo_up), std::move(spo_dn)); From 4a1fb172eb7fdd13094b84be8d31fe78d25eb09f Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Thu, 30 Jun 2022 14:58:25 -0400 Subject: [PATCH 18/34] rename test file --- .../batch/{short-heg-SJ-batched.xml => short-heg-SJ-batch.xml} | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) rename tests/heg/heg_14_gamma/batch/{short-heg-SJ-batched.xml => short-heg-SJ-batch.xml} (98%) diff --git a/tests/heg/heg_14_gamma/batch/short-heg-SJ-batched.xml b/tests/heg/heg_14_gamma/batch/short-heg-SJ-batch.xml similarity index 98% rename from tests/heg/heg_14_gamma/batch/short-heg-SJ-batched.xml rename to tests/heg/heg_14_gamma/batch/short-heg-SJ-batch.xml index 7ff31840f6..2e313a1ebc 100644 --- a/tests/heg/heg_14_gamma/batch/short-heg-SJ-batched.xml +++ b/tests/heg/heg_14_gamma/batch/short-heg-SJ-batch.xml @@ -1,6 +1,6 @@ - + batch From 557f73bd3c23c3edbad4503cac7b254ed2d402b1 Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Sun, 3 Jul 2022 13:58:28 -0400 Subject: [PATCH 19/34] default twist --- src/Particle/LongRange/KContainer.cpp | 2 +- src/Particle/LongRange/KContainer.h | 2 +- src/Particle/LongRange/LRBreakup.h | 3 +-- src/Particle/SimulationCell.cpp | 3 +-- 4 files changed, 4 insertions(+), 6 deletions(-) diff --git a/src/Particle/LongRange/KContainer.cpp b/src/Particle/LongRange/KContainer.cpp index 4aa5c02aec..32e7f7bf0b 100644 --- a/src/Particle/LongRange/KContainer.cpp +++ b/src/Particle/LongRange/KContainer.cpp @@ -21,7 +21,7 @@ namespace qmcplusplus { -void KContainer::updateKLists(const ParticleLayout& lattice, RealType kc, unsigned ndim, PosType twist, bool useSphere) +void KContainer::updateKLists(const ParticleLayout& lattice, RealType kc, unsigned ndim, const PosType& twist, bool useSphere) { kcutoff = kc; if (kcutoff <= 0.0) diff --git a/src/Particle/LongRange/KContainer.h b/src/Particle/LongRange/KContainer.h index 9ad6ac19a8..49d32ae7ad 100644 --- a/src/Particle/LongRange/KContainer.h +++ b/src/Particle/LongRange/KContainer.h @@ -72,7 +72,7 @@ class KContainer : public QMCTraits * @param kc cutoff radius in the K * @param useSphere if true, use the |K| */ - void updateKLists(const ParticleLayout& lattice, RealType kc, unsigned ndim, PosType twist, bool useSphere = true); + void updateKLists(const ParticleLayout& lattice, RealType kc, unsigned ndim, const PosType& twist=PosType(), bool useSphere = true); const auto& get_kpts_cart_soa() const { return kpts_cart_soa_; } private: diff --git a/src/Particle/LongRange/LRBreakup.h b/src/Particle/LongRange/LRBreakup.h index 2f49c82a23..38233d0b18 100644 --- a/src/Particle/LongRange/LRBreakup.h +++ b/src/Particle/LongRange/LRBreakup.h @@ -334,8 +334,7 @@ int LRBreakup::SetupKVecs(mRealType kc, mRealType kcont, mRealType { //Add low |k| ( < kcont) k-points with exact degeneracy KContainer kexact; - PosType twist(0); - kexact.updateKLists(Basis.get_Lattice(), kcont, Basis.get_Lattice().ndim, twist); + kexact.updateKLists(Basis.get_Lattice(), kcont, Basis.get_Lattice().ndim); bool findK = true; mRealType kc2 = kc * kc; //use at least one shell diff --git a/src/Particle/SimulationCell.cpp b/src/Particle/SimulationCell.cpp index b90805b61f..cbd9b5f053 100644 --- a/src/Particle/SimulationCell.cpp +++ b/src/Particle/SimulationCell.cpp @@ -59,8 +59,7 @@ void SimulationCell::resetLRBox() app_log() << "--------------------------------------- " << std::endl; } - QMCTraits::PosType twist(0); - k_lists_.updateKLists(LRBox_, LRBox_.LR_kc, LRBox_.ndim, twist); + k_lists_.updateKLists(LRBox_, LRBox_.LR_kc, LRBox_.ndim); } } } From 9845e2bec8c393ddcc94586bc34a1a705dbcc367 Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Sun, 3 Jul 2022 13:59:04 -0400 Subject: [PATCH 20/34] throw runtime_error --- src/QMCWaveFunctions/ElectronGas/FreeOrbital.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/QMCWaveFunctions/ElectronGas/FreeOrbital.h b/src/QMCWaveFunctions/ElectronGas/FreeOrbital.h index a5c903ecd0..d61ee3a2fb 100644 --- a/src/QMCWaveFunctions/ElectronGas/FreeOrbital.h +++ b/src/QMCWaveFunctions/ElectronGas/FreeOrbital.h @@ -67,12 +67,12 @@ class FreeOrbital : public SPOSet HessMatrix& d2phi_mat, GGGMatrix& d3phi_mat) override; + void report(const std::string& pad) const override; // ---- begin required overrides std::unique_ptr makeClone() const override {return std::make_unique(*this);} void resetParameters(const opt_variables_type& optVariables) override {} //called by BFTrans} - void setOrbitalSetSize(int norbs) override {APP_ABORT("not implemented")}; + void setOrbitalSetSize(int norbs) override {throw std::runtime_error("not implemented");} // required overrides end ---- - void report(const std::string& pad) const override; private: const std::vector kvecs; // kvecs vectors const int mink; // minimum k index From cdb7b0e4ffe1c8c715729be0b3b3587a3863f483 Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Sun, 3 Jul 2022 15:01:21 -0400 Subject: [PATCH 21/34] start KContainer tests --- .../LongRange/tests/test_kcontainer.cpp | 69 +++++++++++++++++++ 1 file changed, 69 insertions(+) create mode 100644 src/Particle/LongRange/tests/test_kcontainer.cpp diff --git a/src/Particle/LongRange/tests/test_kcontainer.cpp b/src/Particle/LongRange/tests/test_kcontainer.cpp new file mode 100644 index 0000000000..152b7ae07d --- /dev/null +++ b/src/Particle/LongRange/tests/test_kcontainer.cpp @@ -0,0 +1,69 @@ +#include "catch.hpp" + +#include "Particle/LongRange/KContainer.h" + +namespace qmcplusplus +{ + +TEST_CASE("kcontainer at gamma", "[longrange]") +{ + const int ndim = 3; + const double alat = 1.0; + const double blat = 2*M_PI/alat; + const std::vector kcs = {blat, std::sqrt(2)*blat, std::sqrt(3)*blat}; + const std::vector nks = {6, 18, 26}; + + CrystalLattice lattice; + lattice.BoxBConds = true; + lattice.R.diagonal(1.0); + lattice.set(lattice.R); // compute Rv and Gv from R + + KContainer klists; + // check first 3 shells of kvectors + for (int ik=0;ik Date: Sun, 3 Jul 2022 15:07:42 -0400 Subject: [PATCH 22/34] kvecs in 2D --- .../LongRange/tests/test_kcontainer.cpp | 53 +++++++++++++++++-- 1 file changed, 50 insertions(+), 3 deletions(-) diff --git a/src/Particle/LongRange/tests/test_kcontainer.cpp b/src/Particle/LongRange/tests/test_kcontainer.cpp index 152b7ae07d..6330312cad 100644 --- a/src/Particle/LongRange/tests/test_kcontainer.cpp +++ b/src/Particle/LongRange/tests/test_kcontainer.cpp @@ -5,21 +5,21 @@ namespace qmcplusplus { -TEST_CASE("kcontainer at gamma", "[longrange]") +TEST_CASE("kcontainer at gamma in 3D", "[longrange]") { const int ndim = 3; const double alat = 1.0; const double blat = 2*M_PI/alat; + + // check first 3 shells of kvectors const std::vector kcs = {blat, std::sqrt(2)*blat, std::sqrt(3)*blat}; const std::vector nks = {6, 18, 26}; CrystalLattice lattice; - lattice.BoxBConds = true; lattice.R.diagonal(1.0); lattice.set(lattice.R); // compute Rv and Gv from R KContainer klists; - // check first 3 shells of kvectors for (int ik=0;ik kcs = {blat, std::sqrt(2)*blat, 2*blat}; + const std::vector nks = {4, 8, 12}; + + CrystalLattice lattice; + lattice.R.diagonal(1.0); + lattice.set(lattice.R); // compute Rv and Gv from R + + KContainer klists; + for (int ik=0;ik Date: Sun, 3 Jul 2022 15:26:58 -0400 Subject: [PATCH 23/34] const klist --- src/Particle/LongRange/StructFact.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Particle/LongRange/StructFact.h b/src/Particle/LongRange/StructFact.h index 85da4a6308..9ac72cd5e1 100644 --- a/src/Particle/LongRange/StructFact.h +++ b/src/Particle/LongRange/StructFact.h @@ -72,7 +72,7 @@ class StructFact : public QMCTraits bool isStorePerParticle() const { return StorePerParticle; } /// accessor of k_lists_ - KContainer getKLists() const { return k_lists_; } + const KContainer getKLists() const { return k_lists_; } private: /// Compute all rhok elements from the start From 4fdd530bd8ad88d445a9db22a69b4a6587453929 Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Sun, 3 Jul 2022 15:35:07 -0400 Subject: [PATCH 24/34] twist 2D --- .../LongRange/tests/test_kcontainer.cpp | 48 ++++++++++++++++++- 1 file changed, 46 insertions(+), 2 deletions(-) diff --git a/src/Particle/LongRange/tests/test_kcontainer.cpp b/src/Particle/LongRange/tests/test_kcontainer.cpp index 6330312cad..051461bd57 100644 --- a/src/Particle/LongRange/tests/test_kcontainer.cpp +++ b/src/Particle/LongRange/tests/test_kcontainer.cpp @@ -1,9 +1,11 @@ #include "catch.hpp" #include "Particle/LongRange/KContainer.h" +#include "OhmmsPETE/TinyVector.h" namespace qmcplusplus { +using PosType = TinyVector; TEST_CASE("kcontainer at gamma in 3D", "[longrange]") { @@ -27,7 +29,7 @@ TEST_CASE("kcontainer at gamma in 3D", "[longrange]") CHECK(klists.kpts.size() == nks[ik]); } const int mxk = klists.kpts.size(); - int gvecs[mxk][OHMMS_DIM] = { + int gvecs[26][3] = { {-1, 0, 0}, { 0, -1, 0}, { 0, 0, -1}, @@ -88,7 +90,7 @@ TEST_CASE("kcontainer at gamma in 2D", "[longrange]") CHECK(klists.kpts.size() == nks[ik]); } const int mxk = klists.kpts.size(); - int gvecs[mxk][OHMMS_DIM] = { + int gvecs[12][3] = { {-1, 0, 0}, { 0, -1, 0}, { 0, 1, 0}, @@ -113,4 +115,46 @@ TEST_CASE("kcontainer at gamma in 2D", "[longrange]") } } +TEST_CASE("kcontainer at twist in 2D", "[longrange]") +{ + const int ndim = 2; + const double alat = 1.0; + const double blat = 2*M_PI/alat; + + // twist one shell of kvectors + const double kc = blat+1e-6; + + CrystalLattice lattice; + lattice.R.diagonal(1.0); + lattice.set(lattice.R); // compute Rv and Gv from R + + KContainer klists; + + PosType twist; + twist[0] = 0.1; + klists.updateKLists(lattice, kc, ndim, twist); + CHECK(klists.kpts.size() == 1); + CHECK(klists.kpts_cart[0][0] == Approx(blat*(twist[0]-1))); + + twist[1] = 0.1; + klists.updateKLists(lattice, kc, ndim, twist); + CHECK(klists.kpts.size() == 2); + CHECK(klists.kpts_cart[0][0] == Approx(blat*(twist[0]-1))); + CHECK(klists.kpts_cart[0][1] == Approx(blat*twist[1])); + CHECK(klists.kpts_cart[1][0] == Approx(blat*(twist[0]))); + CHECK(klists.kpts_cart[1][1] == Approx(blat*(twist[1]-1))); + + twist = {-0.5, 0.5, 0}; + klists.updateKLists(lattice, kc, ndim, twist); + CHECK(klists.kpts.size() == 3); + //for (int ik=0;ik<3;ik++) + // app_log() << klists.kpts_cart[ik] << std::endl; + CHECK(klists.kpts_cart[0][0] == Approx(blat*(twist[0]-0))); + CHECK(klists.kpts_cart[0][1] == Approx(blat*(twist[1]-1))); + CHECK(klists.kpts_cart[1][0] == Approx(blat*(twist[0]+1))); + CHECK(klists.kpts_cart[1][1] == Approx(blat*(twist[1]-1))); + CHECK(klists.kpts_cart[2][0] == Approx(blat*(twist[0]+1))); + CHECK(klists.kpts_cart[2][1] == Approx(blat*twist[1])); +} + } // qmcplusplus From 1229433a557743b89c45f740cbb099003fae0ce4 Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Sun, 3 Jul 2022 15:45:46 -0400 Subject: [PATCH 25/34] twist in 3D --- .../LongRange/tests/test_kcontainer.cpp | 34 +++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/src/Particle/LongRange/tests/test_kcontainer.cpp b/src/Particle/LongRange/tests/test_kcontainer.cpp index 051461bd57..aa0f46dd51 100644 --- a/src/Particle/LongRange/tests/test_kcontainer.cpp +++ b/src/Particle/LongRange/tests/test_kcontainer.cpp @@ -68,6 +68,40 @@ TEST_CASE("kcontainer at gamma in 3D", "[longrange]") } } +TEST_CASE("kcontainer at twist in 3D", "[longrange]") +{ + const int ndim = 3; + const double alat = 1.0; + const double blat = 2*M_PI/alat; + + // twist one shell of kvectors + const double kc = blat+1e-6; + + CrystalLattice lattice; + lattice.R.diagonal(1.0); + lattice.set(lattice.R); // compute Rv and Gv from R + + KContainer klists; + + PosType twist; + twist[0] = 0.1; + klists.updateKLists(lattice, kc, ndim, twist); + CHECK(klists.kpts.size() == 1); + CHECK(klists.kpts_cart[0][0] == Approx(blat*(twist[0]-1))); + + twist = {-0.5, 0, 0.5}; + klists.updateKLists(lattice, kc, ndim, twist); + int gvecs[3][3] = { + { 0, 0, -1}, + { 1, 0, -1}, + { 1, 0, 0} + }; + CHECK(klists.kpts.size() == 3); + for (int ik=0;ik Date: Sun, 3 Jul 2022 15:46:12 -0400 Subject: [PATCH 26/34] run KContainer tests --- src/Particle/LongRange/tests/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Particle/LongRange/tests/CMakeLists.txt b/src/Particle/LongRange/tests/CMakeLists.txt index 2899a2355a..cc7aac79e9 100644 --- a/src/Particle/LongRange/tests/CMakeLists.txt +++ b/src/Particle/LongRange/tests/CMakeLists.txt @@ -13,7 +13,7 @@ set(SRC_DIR longrange) set(UTEST_EXE test_${SRC_DIR}) set(UTEST_NAME deterministic-unit_test_${SRC_DIR}) -add_executable(${UTEST_EXE} test_lrhandler.cpp test_ewald3d.cpp test_temp.cpp test_srcoul.cpp test_StructFact.cpp) +add_executable(${UTEST_EXE} test_lrhandler.cpp test_ewald3d.cpp test_temp.cpp test_srcoul.cpp test_StructFact.cpp test_kcontainer.cpp) target_link_libraries(${UTEST_EXE} catch_main qmcparticle) if(USE_OBJECT_TARGET) target_link_libraries(${UTEST_EXE} qmcutil qmcparticle_omptarget) From 199bb8520a8a4c9ff641133b8e12ef7f4dfd9d2d Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Sun, 3 Jul 2022 15:49:45 -0400 Subject: [PATCH 27/34] clang format --- src/Particle/LongRange/KContainer.cpp | 8 +- src/Particle/LongRange/KContainer.h | 7 +- src/Particle/LongRange/StructFact.h | 1 - .../ElectronGas/FreeOrbital.cpp | 162 ++++++++---------- .../ElectronGas/FreeOrbital.h | 62 +++---- .../ElectronGas/FreeOrbitalBuilder.cpp | 31 ++-- .../ElectronGas/FreeOrbitalBuilder.h | 5 +- src/QMCWaveFunctions/ElectronGas/HEGGrid.h | 4 +- 8 files changed, 133 insertions(+), 147 deletions(-) diff --git a/src/Particle/LongRange/KContainer.cpp b/src/Particle/LongRange/KContainer.cpp index 32e7f7bf0b..738c3f73db 100644 --- a/src/Particle/LongRange/KContainer.cpp +++ b/src/Particle/LongRange/KContainer.cpp @@ -21,7 +21,11 @@ namespace qmcplusplus { -void KContainer::updateKLists(const ParticleLayout& lattice, RealType kc, unsigned ndim, const PosType& twist, bool useSphere) +void KContainer::updateKLists(const ParticleLayout& lattice, + RealType kc, + unsigned ndim, + const PosType& twist, + bool useSphere) { kcutoff = kc; if (kcutoff <= 0.0) @@ -125,7 +129,7 @@ void KContainer::BuildKLists(const ParticleLayout& lattice, PosType twist, bool if (i == 0 && j == 0 && k == 0) continue; //Convert kvec to Cartesian - kvec_cart = lattice.k_cart(kvec+twist); + kvec_cart = lattice.k_cart(kvec + twist); //Find modk modk2 = dot(kvec_cart, kvec_cart); if (modk2 > kcut2) diff --git a/src/Particle/LongRange/KContainer.h b/src/Particle/LongRange/KContainer.h index 49d32ae7ad..d30dcc9687 100644 --- a/src/Particle/LongRange/KContainer.h +++ b/src/Particle/LongRange/KContainer.h @@ -72,9 +72,14 @@ class KContainer : public QMCTraits * @param kc cutoff radius in the K * @param useSphere if true, use the |K| */ - void updateKLists(const ParticleLayout& lattice, RealType kc, unsigned ndim, const PosType& twist=PosType(), bool useSphere = true); + void updateKLists(const ParticleLayout& lattice, + RealType kc, + unsigned ndim, + const PosType& twist = PosType(), + bool useSphere = true); const auto& get_kpts_cart_soa() const { return kpts_cart_soa_; } + private: /** compute approximate parallelpiped that surrounds kc * @param lattice supercell diff --git a/src/Particle/LongRange/StructFact.h b/src/Particle/LongRange/StructFact.h index 9ac72cd5e1..c5430b02a8 100644 --- a/src/Particle/LongRange/StructFact.h +++ b/src/Particle/LongRange/StructFact.h @@ -22,7 +22,6 @@ namespace qmcplusplus { - class KContainer; /** @ingroup longrange *\brief Calculates the structure-factor for a particle set diff --git a/src/QMCWaveFunctions/ElectronGas/FreeOrbital.cpp b/src/QMCWaveFunctions/ElectronGas/FreeOrbital.cpp index e06ca0eea2..d04370adb5 100644 --- a/src/QMCWaveFunctions/ElectronGas/FreeOrbital.cpp +++ b/src/QMCWaveFunctions/ElectronGas/FreeOrbital.cpp @@ -2,40 +2,34 @@ namespace qmcplusplus { - FreeOrbital::FreeOrbital(const std::vector& kpts_cart) - : kvecs(kpts_cart), + : kvecs(kpts_cart), #ifdef QMC_COMPLEX - mink(0), // first k at twist may not be 0 + mink(0), // first k at twist may not be 0 #else - mink(1), // treat k=0 as special case + mink(1), // treat k=0 as special case #endif - maxk(kpts_cart.size()) + maxk(kpts_cart.size()) { #ifdef QMC_COMPLEX OrbitalSetSize = maxk; #else - OrbitalSetSize = 2*maxk-1; // k=0 has no (cos, sin) split + OrbitalSetSize = 2 * maxk - 1; // k=0 has no (cos, sin) split #endif k2neg.resize(maxk); - for (int ik=0; ik makeClone() const override {return std::make_unique(*this);} + std::unique_ptr makeClone() const override { return std::make_unique(*this); } void resetParameters(const opt_variables_type& optVariables) override {} //called by BFTrans} - void setOrbitalSetSize(int norbs) override {throw std::runtime_error("not implemented");} + void setOrbitalSetSize(int norbs) override { throw std::runtime_error("not implemented"); } // required overrides end ---- private: const std::vector kvecs; // kvecs vectors - const int mink; // minimum k index - const int maxk; // maximum number of kvecs vectors - std::vector k2neg; // minus kvecs^2 + const int mink; // minimum k index + const int maxk; // maximum number of kvecs vectors + std::vector k2neg; // minus kvecs^2 }; -} // qmcplusplus +} // namespace qmcplusplus #endif diff --git a/src/QMCWaveFunctions/ElectronGas/FreeOrbitalBuilder.cpp b/src/QMCWaveFunctions/ElectronGas/FreeOrbitalBuilder.cpp index 32ee034b5a..f2eb944592 100644 --- a/src/QMCWaveFunctions/ElectronGas/FreeOrbitalBuilder.cpp +++ b/src/QMCWaveFunctions/ElectronGas/FreeOrbitalBuilder.cpp @@ -6,10 +6,8 @@ namespace qmcplusplus { - FreeOrbitalBuilder::FreeOrbitalBuilder(ParticleSet& els, Communicate* comm, xmlNodePtr cur) - : SPOSetBuilder("PW", comm), - targetPtcl(els) + : SPOSetBuilder("PW", comm), targetPtcl(els) {} std::unique_ptr FreeOrbitalBuilder::createSPOSetFromXML(xmlNodePtr cur) @@ -30,10 +28,10 @@ std::unique_ptr FreeOrbitalBuilder::createSPOSetFromXML(xmlNodePtr cur) app_log() << "twist fraction = " << twist << std::endl; app_log() << "twist cartesian = " << tvec << std::endl; #else - npw = std::ceil((norb+1.0)/2); - for (int ldim=0;ldim 1e-16) + if (std::abs(twist[ldim]) > 1e-16) throw std::runtime_error("no twist for real orbitals"); } #endif @@ -48,24 +46,26 @@ std::unique_ptr FreeOrbitalBuilder::createSPOSetFromXML(xmlNodePtr cur) // k0 is not in kpts_cart kpts[0] = tvec; #ifdef QMC_COMPLEX - for (int ik=1;ik mkidx(npw, 0); int ik = 1; - for (int jk=0;jk= npw) break; + if (ik >= npw) + break; } #endif auto sposet = std::make_unique(kpts); @@ -75,11 +75,12 @@ std::unique_ptr FreeOrbitalBuilder::createSPOSetFromXML(xmlNodePtr cur) bool FreeOrbitalBuilder::in_list(const int j, const std::vector l) { - for (int i=0;i createSPOSetFromXML(xmlNodePtr cur) override; + private: ParticleSet& targetPtcl; bool in_list(const int j, const std::vector l); }; -} // qmcplusplus +} // namespace qmcplusplus #endif diff --git a/src/QMCWaveFunctions/ElectronGas/HEGGrid.h b/src/QMCWaveFunctions/ElectronGas/HEGGrid.h index 3c8659d036..1e1045064b 100644 --- a/src/QMCWaveFunctions/ElectronGas/HEGGrid.h +++ b/src/QMCWaveFunctions/ElectronGas/HEGGrid.h @@ -21,12 +21,11 @@ namespace qmcplusplus { - //three-d specialization template struct HEGGrid { - using PL_t = CrystalLattice; + using PL_t = CrystalLattice; const PL_t& Lattice; static constexpr std::array n_within_shell{{1, 7, 19, 27, 33, 57, 81, 93, 123, 147, 171, @@ -64,7 +63,6 @@ struct HEGGrid * @param rs_in rs */ inline T getCellLength(int nptcl, T rs_in) const { return std::pow(4.0 * M_PI * nptcl / 3.0, 1.0 / 3.0) * rs_in; } - }; } // namespace qmcplusplus From bbde8ff5ad54220e3eb93e8ec2f1ad8bc3ab980f Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Tue, 5 Jul 2022 08:37:45 -0400 Subject: [PATCH 28/34] getter return const reference --- src/Particle/LongRange/KContainer.cpp | 2 +- src/Particle/LongRange/KContainer.h | 3 ++- src/Particle/LongRange/StructFact.h | 2 +- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/Particle/LongRange/KContainer.cpp b/src/Particle/LongRange/KContainer.cpp index 738c3f73db..72d4c8bd17 100644 --- a/src/Particle/LongRange/KContainer.cpp +++ b/src/Particle/LongRange/KContainer.cpp @@ -102,7 +102,7 @@ void KContainer::findApproxMMax(const ParticleLayout& lattice, unsigned ndim) mmax[1] = 0; } -void KContainer::BuildKLists(const ParticleLayout& lattice, PosType twist, bool useSphere) +void KContainer::BuildKLists(const ParticleLayout& lattice, const PosType& twist, bool useSphere) { TinyVector TempActualMax; TinyVector kvec; diff --git a/src/Particle/LongRange/KContainer.h b/src/Particle/LongRange/KContainer.h index d30dcc9687..eee91affc7 100644 --- a/src/Particle/LongRange/KContainer.h +++ b/src/Particle/LongRange/KContainer.h @@ -70,6 +70,7 @@ class KContainer : public QMCTraits /** update k-vectors * @param sc supercell * @param kc cutoff radius in the K + * @param twist shifts the center of the grid of k-vectors * @param useSphere if true, use the |K| */ void updateKLists(const ParticleLayout& lattice, @@ -86,7 +87,7 @@ class KContainer : public QMCTraits */ void findApproxMMax(const ParticleLayout& lattice, unsigned ndim); /** construct the container for k-vectors */ - void BuildKLists(const ParticleLayout& lattice, PosType twist, bool useSphere); + void BuildKLists(const ParticleLayout& lattice, const PosType& twist, bool useSphere); /** K-vector in Cartesian coordinates in SoA layout */ diff --git a/src/Particle/LongRange/StructFact.h b/src/Particle/LongRange/StructFact.h index c5430b02a8..ad1de84d53 100644 --- a/src/Particle/LongRange/StructFact.h +++ b/src/Particle/LongRange/StructFact.h @@ -71,7 +71,7 @@ class StructFact : public QMCTraits bool isStorePerParticle() const { return StorePerParticle; } /// accessor of k_lists_ - const KContainer getKLists() const { return k_lists_; } + const KContainer& getKLists() const { return k_lists_; } private: /// Compute all rhok elements from the start From 98ad0101aafae755d36ee799c5783292d86eaa83 Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Tue, 5 Jul 2022 08:37:57 -0400 Subject: [PATCH 29/34] fix missing Approx --- src/Particle/LongRange/tests/test_kcontainer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Particle/LongRange/tests/test_kcontainer.cpp b/src/Particle/LongRange/tests/test_kcontainer.cpp index aa0f46dd51..2c3be38785 100644 --- a/src/Particle/LongRange/tests/test_kcontainer.cpp +++ b/src/Particle/LongRange/tests/test_kcontainer.cpp @@ -63,7 +63,7 @@ TEST_CASE("kcontainer at gamma in 3D", "[longrange]") for (int ldim=0;ldim Date: Tue, 5 Jul 2022 09:40:06 -0400 Subject: [PATCH 30/34] fewer parenthesis --- .../ElectronGas/FreeOrbital.cpp | 90 +++++++++---------- 1 file changed, 44 insertions(+), 46 deletions(-) diff --git a/src/QMCWaveFunctions/ElectronGas/FreeOrbital.cpp b/src/QMCWaveFunctions/ElectronGas/FreeOrbital.cpp index d04370adb5..6e320b6314 100644 --- a/src/QMCWaveFunctions/ElectronGas/FreeOrbital.cpp +++ b/src/QMCWaveFunctions/ElectronGas/FreeOrbital.cpp @@ -18,9 +18,7 @@ FreeOrbital::FreeOrbital(const std::vector& kpts_cart) #endif k2neg.resize(maxk); for (int ik = 0; ik < maxk; ik++) - { k2neg[ik] = -dot(kvecs[ik], kvecs[ik]); - } } FreeOrbital::~FreeOrbital() {} @@ -118,11 +116,11 @@ void FreeOrbital::evaluate_notranspose(const ParticleSet& P, dp[ik] = ValueType(-sinkr, coskr) * kvecs[ik]; for (int la = 0; la < OHMMS_DIM; la++) { - (hess[ik])(la, la) = -phi_of_r * (kvecs[ik])[la] * (kvecs[ik])[la]; + hess[ik](la, la) = -phi_of_r * (kvecs[ik])[la] * (kvecs[ik])[la]; for (int lb = la + 1; lb < OHMMS_DIM; lb++) { - (hess[ik])(la, lb) = -phi_of_r * (kvecs[ik])[la] * (kvecs[ik])[lb]; - (hess[ik])(lb, la) = (hess[ik])(la, lb); + hess[ik](la, lb) = -phi_of_r * (kvecs[ik])[la] * (kvecs[ik])[lb]; + hess[ik](lb, la) = hess[ik](la, lb); } } #else @@ -134,14 +132,14 @@ void FreeOrbital::evaluate_notranspose(const ParticleSet& P, dp[j2] = coskr * kvecs[ik]; for (int la = 0; la < OHMMS_DIM; la++) { - (hess[j1])(la, la) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[la]; - (hess[j2])(la, la) = -sinkr * (kvecs[ik])[la] * (kvecs[ik])[la]; + hess[j1](la, la) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[la]; + hess[j2](la, la) = -sinkr * (kvecs[ik])[la] * (kvecs[ik])[la]; for (int lb = la + 1; lb < OHMMS_DIM; lb++) { - (hess[j1])(la, lb) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[lb]; - (hess[j2])(la, lb) = -sinkr * (kvecs[ik])[la] * (kvecs[ik])[lb]; - (hess[j1])(lb, la) = (hess[j1])(la, lb); - (hess[j2])(lb, la) = (hess[j2])(la, lb); + hess[j1](la, lb) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[lb]; + hess[j2](la, lb) = -sinkr * (kvecs[ik])[la] * (kvecs[ik])[lb]; + hess[j1](lb, la) = hess[j1](la, lb); + hess[j2](lb, la) = hess[j2](la, lb); } } #endif @@ -182,11 +180,11 @@ void FreeOrbital::evaluate_notranspose(const ParticleSet& P, dp[ik] = compi * phi_of_r * kvecs[ik]; for (int la = 0; la < OHMMS_DIM; la++) { - (hess[ik])(la, la) = -phi_of_r * (kvecs[ik])[la] * (kvecs[ik])[la]; + hess[ik](la, la) = -phi_of_r * (kvecs[ik])[la] * (kvecs[ik])[la]; for (int lb = la + 1; lb < OHMMS_DIM; lb++) { - (hess[ik])(la, lb) = -phi_of_r * (kvecs[ik])[la] * (kvecs[ik])[lb]; - (hess[ik])(lb, la) = (hess[ik])(la, lb); + hess[ik](la, lb) = -phi_of_r * (kvecs[ik])[la] * (kvecs[ik])[lb]; + hess[ik](lb, la) = hess[ik](la, lb); } } for (int la = 0; la < OHMMS_DIM; la++) @@ -202,42 +200,42 @@ void FreeOrbital::evaluate_notranspose(const ParticleSet& P, dp[j2] = coskr * kvecs[ik]; for (int la = 0; la < OHMMS_DIM; la++) { - (hess[j1])(la, la) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[la]; - (hess[j2])(la, la) = -sinkr * (kvecs[ik])[la] * (kvecs[ik])[la]; - ((ggg[j1])[la])(la, la) = sinkr * (kvecs[ik])[la] * (kvecs[ik])[la] * (kvecs[ik])[la]; - ((ggg[j2])[la])(la, la) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[la] * (kvecs[ik])[la]; + hess[j1](la, la) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[la]; + hess[j2](la, la) = -sinkr * (kvecs[ik])[la] * (kvecs[ik])[la]; + ggg[j1][la](la, la) = sinkr * (kvecs[ik])[la] * (kvecs[ik])[la] * (kvecs[ik])[la]; + ggg[j2][la](la, la) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[la] * (kvecs[ik])[la]; for (int lb = la + 1; lb < OHMMS_DIM; lb++) { - (hess[j1])(la, lb) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[lb]; - (hess[j2])(la, lb) = -sinkr * (kvecs[ik])[la] * (kvecs[ik])[lb]; - (hess[j1])(lb, la) = (hess[j1])(la, lb); - (hess[j2])(lb, la) = (hess[j2])(la, lb); - ((ggg[j1])[la])(lb, la) = sinkr * (kvecs[ik])[la] * (kvecs[ik])[lb] * (kvecs[ik])[la]; - ((ggg[j2])[la])(lb, la) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[lb] * (kvecs[ik])[la]; - ((ggg[j1])[la])(la, lb) = ((ggg[j1])[la])(lb, la); - ((ggg[j2])[la])(la, lb) = ((ggg[j2])[la])(lb, la); - ((ggg[j1])[lb])(la, la) = ((ggg[j1])[la])(lb, la); - ((ggg[j2])[lb])(la, la) = ((ggg[j2])[la])(lb, la); - ((ggg[j1])[la])(lb, lb) = sinkr * (kvecs[ik])[la] * (kvecs[ik])[lb] * (kvecs[ik])[lb]; - ((ggg[j2])[la])(lb, lb) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[lb] * (kvecs[ik])[lb]; - ((ggg[j1])[lb])(la, lb) = ((ggg[j1])[la])(lb, lb); - ((ggg[j2])[lb])(la, lb) = ((ggg[j2])[la])(lb, lb); - ((ggg[j1])[lb])(lb, la) = ((ggg[j1])[la])(lb, lb); - ((ggg[j2])[lb])(lb, la) = ((ggg[j2])[la])(lb, lb); + hess[j1](la, lb) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[lb]; + hess[j2](la, lb) = -sinkr * (kvecs[ik])[la] * (kvecs[ik])[lb]; + hess[j1](lb, la) = hess[j1](la, lb); + hess[j2](lb, la) = hess[j2](la, lb); + ggg[j1][la](lb, la) = sinkr * (kvecs[ik])[la] * (kvecs[ik])[lb] * (kvecs[ik])[la]; + ggg[j2][la](lb, la) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[lb] * (kvecs[ik])[la]; + ggg[j1][la](la, lb) = ggg[j1][la](lb, la); + ggg[j2][la](la, lb) = ggg[j2][la](lb, la); + ggg[j1][lb](la, la) = ggg[j1][la](lb, la); + ggg[j2][lb](la, la) = ggg[j2][la](lb, la); + ggg[j1][la](lb, lb) = sinkr * (kvecs[ik])[la] * (kvecs[ik])[lb] * (kvecs[ik])[lb]; + ggg[j2][la](lb, lb) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[lb] * (kvecs[ik])[lb]; + ggg[j1][lb](la, lb) = ggg[j1][la](lb, lb); + ggg[j2][lb](la, lb) = ggg[j2][la](lb, lb); + ggg[j1][lb](lb, la) = ggg[j1][la](lb, lb); + ggg[j2][lb](lb, la) = ggg[j2][la](lb, lb); for (int lc = lb + 1; lc < OHMMS_DIM; lc++) { - ((ggg[j1])[la])(lb, lc) = sinkr * (kvecs[ik])[la] * (kvecs[ik])[lb] * (kvecs[ik])[lc]; - ((ggg[j2])[la])(lb, lc) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[lb] * (kvecs[ik])[lc]; - ((ggg[j1])[la])(lc, lb) = ((ggg[j1])[la])(lb, lc); - ((ggg[j2])[la])(lc, lb) = ((ggg[j2])[la])(lb, lc); - ((ggg[j1])[lb])(la, lc) = ((ggg[j1])[la])(lb, lc); - ((ggg[j2])[lb])(la, lc) = ((ggg[j2])[la])(lb, lc); - ((ggg[j1])[lb])(lc, la) = ((ggg[j1])[la])(lb, lc); - ((ggg[j2])[lb])(lc, la) = ((ggg[j2])[la])(lb, lc); - ((ggg[j1])[lc])(la, lb) = ((ggg[j1])[la])(lb, lc); - ((ggg[j2])[lc])(la, lb) = ((ggg[j2])[la])(lb, lc); - ((ggg[j1])[lc])(lb, la) = ((ggg[j1])[la])(lb, lc); - ((ggg[j2])[lc])(lb, la) = ((ggg[j2])[la])(lb, lc); + ggg[j1][la](lb, lc) = sinkr * (kvecs[ik])[la] * (kvecs[ik])[lb] * (kvecs[ik])[lc]; + ggg[j2][la](lb, lc) = -coskr * (kvecs[ik])[la] * (kvecs[ik])[lb] * (kvecs[ik])[lc]; + ggg[j1][la](lc, lb) = ggg[j1][la](lb, lc); + ggg[j2][la](lc, lb) = ggg[j2][la](lb, lc); + ggg[j1][lb](la, lc) = ggg[j1][la](lb, lc); + ggg[j2][lb](la, lc) = ggg[j2][la](lb, lc); + ggg[j1][lb](lc, la) = ggg[j1][la](lb, lc); + ggg[j2][lb](lc, la) = ggg[j2][la](lb, lc); + ggg[j1][lc](la, lb) = ggg[j1][la](lb, lc); + ggg[j2][lc](la, lb) = ggg[j2][la](lb, lc); + ggg[j1][lc](lb, la) = ggg[j1][la](lb, lc); + ggg[j2][lc](lb, la) = ggg[j2][la](lb, lc); } } } From ad592eb08cea04a3b17826699280d4924ffa29df Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Tue, 5 Jul 2022 09:42:03 -0400 Subject: [PATCH 31/34] add license header --- src/Particle/LongRange/tests/test_kcontainer.cpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/Particle/LongRange/tests/test_kcontainer.cpp b/src/Particle/LongRange/tests/test_kcontainer.cpp index 2c3be38785..93bc36cbfd 100644 --- a/src/Particle/LongRange/tests/test_kcontainer.cpp +++ b/src/Particle/LongRange/tests/test_kcontainer.cpp @@ -1,3 +1,13 @@ +////////////////////////////////////////////////////////////////////////////////////// +// This file is distributed under the University of Illinois/NCSA Open Source License. +// See LICENSE file in top directory for details. +// +// Copyright (c) 2022 QMCPACK developers. +// +// File developed by: Yubo "Paul" Yang, yubo.paul.yang@gmail.com, CCQ @ Flatiron +// +// File created by: Yubo "Paul" Yang, yubo.paul.yang@gmail.com, CCQ @ Flatiron +////////////////////////////////////////////////////////////////////////////////////// #include "catch.hpp" #include "Particle/LongRange/KContainer.h" From bf73840c55fa73e81affba3048ab2abf6c50b26f Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Tue, 5 Jul 2022 09:43:27 -0400 Subject: [PATCH 32/34] Ewald2D license headers --- src/QMCHamiltonians/tests/test_ewald2d.cpp | 10 ++++++++++ src/QMCHamiltonians/tests/test_ewald_quasi2d.cpp | 10 ++++++++++ 2 files changed, 20 insertions(+) diff --git a/src/QMCHamiltonians/tests/test_ewald2d.cpp b/src/QMCHamiltonians/tests/test_ewald2d.cpp index 652c25f95d..c5908ee642 100644 --- a/src/QMCHamiltonians/tests/test_ewald2d.cpp +++ b/src/QMCHamiltonians/tests/test_ewald2d.cpp @@ -1,3 +1,13 @@ +////////////////////////////////////////////////////////////////////////////////////// +// This file is distributed under the University of Illinois/NCSA Open Source License. +// See LICENSE file in top directory for details. +// +// Copyright (c) 2022 QMCPACK developers. +// +// File developed by: Yubo "Paul" Yang, yubo.paul.yang@gmail.com, CCQ @ Flatiron +// +// File created by: Yubo "Paul" Yang, yubo.paul.yang@gmail.com, CCQ @ Flatiron +////////////////////////////////////////////////////////////////////////////////////// #include "catch.hpp" #include "Particle/ParticleSet.h" diff --git a/src/QMCHamiltonians/tests/test_ewald_quasi2d.cpp b/src/QMCHamiltonians/tests/test_ewald_quasi2d.cpp index 8855b9bd4f..7d9d6848ad 100644 --- a/src/QMCHamiltonians/tests/test_ewald_quasi2d.cpp +++ b/src/QMCHamiltonians/tests/test_ewald_quasi2d.cpp @@ -1,3 +1,13 @@ +////////////////////////////////////////////////////////////////////////////////////// +// This file is distributed under the University of Illinois/NCSA Open Source License. +// See LICENSE file in top directory for details. +// +// Copyright (c) 2022 QMCPACK developers. +// +// File developed by: Yubo "Paul" Yang, yubo.paul.yang@gmail.com, CCQ @ Flatiron +// +// File created by: Yubo "Paul" Yang, yubo.paul.yang@gmail.com, CCQ @ Flatiron +////////////////////////////////////////////////////////////////////////////////////// #include "catch.hpp" #include "Particle/ParticleSet.h" From fc38389eb8a626ef90e2ce6d92f898b38a632c3f Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Tue, 5 Jul 2022 09:44:56 -0400 Subject: [PATCH 33/34] clang format --- .../LongRange/tests/test_kcontainer.cpp | 136 +++++++----------- 1 file changed, 49 insertions(+), 87 deletions(-) diff --git a/src/Particle/LongRange/tests/test_kcontainer.cpp b/src/Particle/LongRange/tests/test_kcontainer.cpp index 93bc36cbfd..337c00dd66 100644 --- a/src/Particle/LongRange/tests/test_kcontainer.cpp +++ b/src/Particle/LongRange/tests/test_kcontainer.cpp @@ -19,73 +19,49 @@ using PosType = TinyVector; TEST_CASE("kcontainer at gamma in 3D", "[longrange]") { - const int ndim = 3; + const int ndim = 3; const double alat = 1.0; - const double blat = 2*M_PI/alat; + const double blat = 2 * M_PI / alat; // check first 3 shells of kvectors - const std::vector kcs = {blat, std::sqrt(2)*blat, std::sqrt(3)*blat}; - const std::vector nks = {6, 18, 26}; + const std::vector kcs = {blat, std::sqrt(2) * blat, std::sqrt(3) * blat}; + const std::vector nks = {6, 18, 26}; CrystalLattice lattice; lattice.R.diagonal(1.0); lattice.set(lattice.R); // compute Rv and Gv from R KContainer klists; - for (int ik=0;ik lattice; lattice.R.diagonal(1.0); @@ -97,76 +73,62 @@ TEST_CASE("kcontainer at twist in 3D", "[longrange]") twist[0] = 0.1; klists.updateKLists(lattice, kc, ndim, twist); CHECK(klists.kpts.size() == 1); - CHECK(klists.kpts_cart[0][0] == Approx(blat*(twist[0]-1))); + CHECK(klists.kpts_cart[0][0] == Approx(blat * (twist[0] - 1))); twist = {-0.5, 0, 0.5}; klists.updateKLists(lattice, kc, ndim, twist); - int gvecs[3][3] = { - { 0, 0, -1}, - { 1, 0, -1}, - { 1, 0, 0} - }; + int gvecs[3][3] = {{0, 0, -1}, {1, 0, -1}, {1, 0, 0}}; CHECK(klists.kpts.size() == 3); - for (int ik=0;ik kcs = {blat, std::sqrt(2)*blat, 2*blat}; - const std::vector nks = {4, 8, 12}; + const std::vector kcs = {blat, std::sqrt(2) * blat, 2 * blat}; + const std::vector nks = {4, 8, 12}; CrystalLattice lattice; lattice.R.diagonal(1.0); lattice.set(lattice.R); // compute Rv and Gv from R KContainer klists; - for (int ik=0;ik lattice; lattice.R.diagonal(1.0); @@ -178,27 +140,27 @@ TEST_CASE("kcontainer at twist in 2D", "[longrange]") twist[0] = 0.1; klists.updateKLists(lattice, kc, ndim, twist); CHECK(klists.kpts.size() == 1); - CHECK(klists.kpts_cart[0][0] == Approx(blat*(twist[0]-1))); + CHECK(klists.kpts_cart[0][0] == Approx(blat * (twist[0] - 1))); twist[1] = 0.1; klists.updateKLists(lattice, kc, ndim, twist); CHECK(klists.kpts.size() == 2); - CHECK(klists.kpts_cart[0][0] == Approx(blat*(twist[0]-1))); - CHECK(klists.kpts_cart[0][1] == Approx(blat*twist[1])); - CHECK(klists.kpts_cart[1][0] == Approx(blat*(twist[0]))); - CHECK(klists.kpts_cart[1][1] == Approx(blat*(twist[1]-1))); + CHECK(klists.kpts_cart[0][0] == Approx(blat * (twist[0] - 1))); + CHECK(klists.kpts_cart[0][1] == Approx(blat * twist[1])); + CHECK(klists.kpts_cart[1][0] == Approx(blat * (twist[0]))); + CHECK(klists.kpts_cart[1][1] == Approx(blat * (twist[1] - 1))); twist = {-0.5, 0.5, 0}; klists.updateKLists(lattice, kc, ndim, twist); CHECK(klists.kpts.size() == 3); //for (int ik=0;ik<3;ik++) // app_log() << klists.kpts_cart[ik] << std::endl; - CHECK(klists.kpts_cart[0][0] == Approx(blat*(twist[0]-0))); - CHECK(klists.kpts_cart[0][1] == Approx(blat*(twist[1]-1))); - CHECK(klists.kpts_cart[1][0] == Approx(blat*(twist[0]+1))); - CHECK(klists.kpts_cart[1][1] == Approx(blat*(twist[1]-1))); - CHECK(klists.kpts_cart[2][0] == Approx(blat*(twist[0]+1))); - CHECK(klists.kpts_cart[2][1] == Approx(blat*twist[1])); + CHECK(klists.kpts_cart[0][0] == Approx(blat * (twist[0] - 0))); + CHECK(klists.kpts_cart[0][1] == Approx(blat * (twist[1] - 1))); + CHECK(klists.kpts_cart[1][0] == Approx(blat * (twist[0] + 1))); + CHECK(klists.kpts_cart[1][1] == Approx(blat * (twist[1] - 1))); + CHECK(klists.kpts_cart[2][0] == Approx(blat * (twist[0] + 1))); + CHECK(klists.kpts_cart[2][1] == Approx(blat * twist[1])); } -} // qmcplusplus +} // namespace qmcplusplus From f4ed10d38e3dfba612ad344c8735ab4062003326 Mon Sep 17 00:00:00 2001 From: "Yubo \"Paul\" Yang" Date: Tue, 5 Jul 2022 10:11:30 -0400 Subject: [PATCH 34/34] one more missed Approx --- src/Particle/LongRange/tests/test_kcontainer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Particle/LongRange/tests/test_kcontainer.cpp b/src/Particle/LongRange/tests/test_kcontainer.cpp index 337c00dd66..b4e7f80597 100644 --- a/src/Particle/LongRange/tests/test_kcontainer.cpp +++ b/src/Particle/LongRange/tests/test_kcontainer.cpp @@ -116,7 +116,7 @@ TEST_CASE("kcontainer at gamma in 2D", "[longrange]") for (int ldim = 0; ldim < ndim; ldim++) { CHECK(klists.kpts[ik][ldim] == gvecs[ik][ldim]); - CHECK(klists.kpts[ik][ldim] * blat == klists.kpts_cart[ik][ldim]); + CHECK(klists.kpts[ik][ldim] * blat == Approx(klists.kpts_cart[ik][ldim])); } } }