FreeOrbital SPOSet

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
 
   <wavefunction name="psi0" target="e">
-    <determinantset type="electron-gas" shell="7" shell2="7" randomize="true">
+  <sposet_builder type="free">
+    <sposet name="spo-ud" size="37" twist="0 0 0"/>
+  </sposet_builder>
+  <determinantset>
+    <slaterdeterminant>
+      <determinant id="updet" sposet="spo-ud"/>
+      <determinant id="dndet" sposet="spo-ud"/>
+    </slaterdeterminant>
   </determinantset>
   <jastrow name="J2" type="Two-Body" function="Bspline" print="no">
     <correlation speciesA="u" speciesB="u" size="8" cusp="0">

src/Particle/LongRange/KContainer.cpp

@@ -21,15 +21,19 @@
 
 namespace qmcplusplus
 {
-void KContainer::updateKLists(const ParticleLayout& lattice, RealType kc, unsigned ndim, bool useSphere)
+void KContainer::updateKLists(const ParticleLayout& lattice,
+                              RealType kc,
+                              unsigned ndim,
+                              const PosType& twist,
+                              bool useSphere)
 {
   kcutoff = kc;
   if (kcutoff <= 0.0)
   {
     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 +102,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<int, DIM + 1> TempActualMax;
   TinyVector<int, DIM> kvec;
@@ -125,7 +129,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)

src/Particle/LongRange/KContainer.h

@@ -72,16 +72,21 @@ 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,
+                    const PosType& twist = PosType(),
+                    bool useSphere       = true);
 
   const auto& get_kpts_cart_soa() const { return kpts_cart_soa_; }
+
 private:
   /** compute approximate parallelpiped that surrounds kc
    * @param lattice supercell
    */
   void findApproxMMax(const ParticleLayout& lattice, unsigned ndim);
   /** 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
    */

src/Particle/LongRange/StructFact.h

@@ -22,7 +22,6 @@
 
 namespace qmcplusplus
 {
-
 class KContainer;
 /** @ingroup longrange
  *\brief Calculates the structure-factor for a particle set
@@ -71,6 +70,9 @@ class StructFact : public QMCTraits
   /// accessor of StorePerParticle
   bool isStorePerParticle() const { return StorePerParticle; }
 
+  /// accessor of k_lists_
+  const KContainer getKLists() const { return k_lists_; }
+
 private:
   /// Compute all rhok elements from the start
   void computeRhok(const ParticleSet& P);

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)

src/Particle/LongRange/tests/test_kcontainer.cpp

@@ -0,0 +1,194 @@
+#include "catch.hpp"
+
+#include "Particle/LongRange/KContainer.h"
+#include "OhmmsPETE/TinyVector.h"
+
+namespace qmcplusplus
+{
+using PosType = TinyVector<double, 3>;
+
+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<double> kcs = {blat, std::sqrt(2)*blat, std::sqrt(3)*blat};
+  const std::vector<int> nks = {6, 18, 26};
+
+  CrystalLattice<OHMMS_PRECISION, OHMMS_DIM> lattice;
+  lattice.R.diagonal(1.0);
+  lattice.set(lattice.R); // compute Rv and Gv from R
+
+  KContainer klists;
+  for (int ik=0;ik<kcs.size();ik++)
+  {
+    const double kc = kcs[ik]+1e-6;
+    klists.updateKLists(lattice, kc, ndim);
+    CHECK(klists.kpts.size() == nks[ik]);
+  }
+  const int mxk = klists.kpts.size();
+  int gvecs[26][3] = {
+    {-1,  0,  0},
+    { 0, -1,  0},
+    { 0,  0, -1},
+    { 0,  0,  1},
+    { 0,  1,  0},
+    { 1,  0,  0},
+    {-1, -1,  0},
+    {-1,  0, -1},
+    {-1,  0,  1},
+    {-1,  1,  0},
+    { 0, -1, -1},
+    { 0, -1,  1},
+    { 0,  1, -1},
+    { 0,  1,  1},
+    { 1, -1,  0},
+    { 1,  0, -1},
+    { 1,  0,  1},
+    { 1,  1,  0},
+    {-1, -1, -1},
+    {-1, -1,  1},
+    {-1,  1, -1},
+    {-1,  1,  1},
+    { 1, -1, -1},
+    { 1, -1,  1},
+    { 1,  1, -1},
+    { 1,  1,  1}
+  };
+
+  for (int ik=0;ik<mxk;ik++)
+  {
+    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]);
+    }
+  }
+}
+
+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<OHMMS_PRECISION, OHMMS_DIM> 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<klists.kpts.size();ik++)
+    for (int ldim=0;ldim<ndim;ldim++)
+      CHECK(klists.kpts_cart[ik][ldim] == Approx(blat*(twist[ldim]+gvecs[ik][ldim])));
+}
+
+TEST_CASE("kcontainer at gamma in 2D", "[longrange]")
+{
+  const int ndim = 2;
+  const double alat = 1.0;
+  const double blat = 2*M_PI/alat;
+
+  // check first 3 shells of kvectors
+  const std::vector<double> kcs = {blat, std::sqrt(2)*blat, 2*blat};
+  const std::vector<int> nks = {4, 8, 12};
+
+  CrystalLattice<OHMMS_PRECISION, OHMMS_DIM> lattice;
+  lattice.R.diagonal(1.0);
+  lattice.set(lattice.R); // compute Rv and Gv from R
+
+  KContainer klists;
+  for (int ik=0;ik<kcs.size();ik++)
+  {
+    const double kc = kcs[ik]+1e-6;
+    klists.updateKLists(lattice, kc, ndim);
+    CHECK(klists.kpts.size() == nks[ik]);
+  }
+  const int mxk = klists.kpts.size();
+  int gvecs[12][3] = {
+    {-1,  0,  0},
+    { 0, -1,  0},
+    { 0,  1,  0},
+    { 1,  0,  0},
+    {-1, -1,  0},
+    {-1,  1,  0},
+    { 1, -1,  0},
+    { 1,  1,  0},
+    {-2,  0,  0},
+    { 0, -2,  0},
+    { 0,  2,  0},
+    { 2,  0,  0},
+  };
+
+  for (int ik=0;ik<mxk;ik++)
+  {
+    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]);
+    }
+  }
+}
+
+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<OHMMS_PRECISION, OHMMS_DIM> 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

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/FreeOrbital.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<EGOSet>(kup, k2up);
-  auto spo_dn = std::make_unique<EGOSet>(kdn, k2dn);
+  auto spo_up = std::make_unique<FreeOrbital>(kup);
+  auto spo_dn = std::make_unique<FreeOrbital>(kdn);
 
   auto spinor_set = std::make_unique<SpinorSet>();
   spinor_set->set_spos(std::move(spo_up), std::move(spo_dn));

src/QMCWaveFunctions/CMakeLists.txt

@@ -63,12 +63,7 @@ set(JASTROW_SRCS
 set(JASTROW_OMPTARGET_SRCS
     Jastrow/J2OMPTarget.cpp
     Jastrow/BsplineFunctor.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)
+set(FERMION_SRCS ${FERMION_SRCS} ElectronGas/FreeOrbital.cpp ElectronGas/FreeOrbitalBuilder.cpp)
 
 # wavefunctions only availbale to 3-dim problems
 if(OHMMS_DIM MATCHES 3)