PWOribitalSetT and PWBasisT

Add FullRealType in SPOSet and RotatedSPOs

Move generic definition after specialization

add implicit implementations

Fix some errors

initial commit of templated PWOribitSetT that compiles

cleanup

templateitze PWBasis as well, as is dependancy

remove inaccurate comment

remove polluted commit

src/QMCWaveFunctions/CMakeLists.txt

@@ -111,9 +111,9 @@ if(OHMMS_DIM MATCHES 3)
   endif(HAVE_EINSPLINE)
 
   # plane wave SPO
-  set(FERMION_SRCS ${FERMION_SRCS} PlaneWave/PWBasis.cpp PlaneWave/PWParameterSet.cpp PlaneWave/PWOrbitalBuilder.cpp)
+  set(FERMION_SRCS ${FERMION_SRCS} PlaneWave/PWBasis.cpp PlaneWave/PWBasisT.cpp PlaneWave/PWParameterSet.cpp PlaneWave/PWOrbitalBuilder.cpp)
   if(QMC_COMPLEX)
-    set(FERMION_SRCS ${FERMION_SRCS} PlaneWave/PWOrbitalSet.cpp)
+    set(FERMION_SRCS ${FERMION_SRCS} PlaneWave/PWOrbitalSet.cpp PlaneWave/PWOrbitalSetT.cpp)
   else()
     set(FERMION_SRCS ${FERMION_SRCS} PlaneWave/PWRealOrbitalSet.cpp PlaneWave/PWRealOrbitalSetT.cpp)
   endif(QMC_COMPLEX)

src/QMCWaveFunctions/PlaneWave/PWBasisT.cpp

@@ -0,0 +1,197 @@
+//////////////////////////////////////////////////////////////////////////////////////
+// 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, [email protected], University of Illinois at Urbana-Champaign
+//                    Jeremy McMinnis, [email protected], University of Illinois at Urbana-Champaign
+//                    Mark A. Berrill, [email protected], Oak Ridge National Laboratory
+//                    Mark Dewing, [email protected], University of Illinois at Urbana-Champaign
+//
+// File created by: Jeongnim Kim, [email protected], University of Illinois at Urbana-Champaign
+//////////////////////////////////////////////////////////////////////////////////////
+
+
+/** @file PWBasisT.cpp
+ * @brief Definition of member functions of Plane-wave basis set
+ */
+#include "PWBasisT.h"
+
+namespace qmcplusplus
+{
+template<class T>
+int PWBasisT<T>::readbasis(hdf_archive& h5basisgroup,
+                       RealType ecutoff,
+                       const ParticleLayout& lat,
+                       const std::string& pwname,
+                       const std::string& pwmultname,
+                       bool resizeContainer)
+{
+  ///make a local copy
+  Lattice = lat;
+  ecut    = ecutoff;
+  app_log() << "  PWBasisT<T>::" << pwmultname << " is found " << std::endl;
+  h5basisgroup.read(gvecs, "/electrons/kpoint_0/gvectors");
+  NumPlaneWaves = std::max(gvecs.size(), kplusgvecs_cart.size());
+  if (NumPlaneWaves == 0)
+  {
+    app_error() << "  PWBasisT<T>::readbasis Basis is missing. Abort " << std::endl;
+    abort(); //FIX_ABORT
+  }
+  if (kplusgvecs_cart.empty())
+  {
+    kplusgvecs_cart.resize(NumPlaneWaves);
+    for (int i = 0; i < NumPlaneWaves; i++)
+      kplusgvecs_cart[i] = Lattice.k_cart(gvecs[i]);
+  }
+  //app_log() << "  Gx Gy Gz " << std::endl;
+  //for(int i=0; i<kplusgvecs_cart.size(); i++)
+  //{
+  //  app_log() << kplusgvecs_cart[i] << std::endl;
+  //}
+  //Now remove elements outside Ecut. At the same time, fill k+G and |k+G| lists.
+  //Also keep track of the original index ordering (using indexmap[]) so that
+  //orbital coefficients can be ordered and trimmed for ecut in the same way.
+  //support older parser
+  if (resizeContainer)
+    reset();
+  //std::copy(gvecs.begin(),gvecs.end(),std::ostream_iterator<GIndex_t>(std::cout,"\n"));
+  return NumPlaneWaves;
+}
+
+template<class T>
+void PWBasisT<T>::setTwistAngle(const PosType& tang)
+{
+  PosType dang   = twist - tang;
+  bool sameTwist = dot(dang, dang) < std::numeric_limits<RealType>::epsilon();
+  if (maxmaxg && sameTwist)
+    return;
+  twist = tang;
+  reset();
+}
+
+template<class T>
+void PWBasisT<T>::reset()
+{
+  trimforecut();
+  //logC.resize(3,2*maxmaxg+1);
+  Z.resize(NumPlaneWaves, 2 + DIM);
+  Zv.resize(NumPlaneWaves);
+  phi.resize(NumPlaneWaves);
+}
+
+/** Remove basis elements if kinetic energy > ecut.
+ *
+ * Keep and indexmap so we know how to match coefficients on read.
+ */
+template<class T>
+void PWBasisT<T>::trimforecut()
+{
+  //Convert the twist angle to Cartesian coordinates.
+  twist_cart = Lattice.k_cart(twist);
+  inputmap.resize(NumPlaneWaves);
+  app_log() << "  PWBasisT<T>::TwistAngle (unit) =" << twist << std::endl;
+  app_log() << "  PWBasisT<T>::TwistAngle (cart) =" << twist_cart << std::endl;
+  app_log() << "  PWBasisT<T>::trimforecut NumPlaneWaves (before) =" << NumPlaneWaves << std::endl;
+  std::vector<GIndex_t> gvecCopy(gvecs);
+  std::vector<PosType> gcartCopy(kplusgvecs_cart);
+  gvecs.clear();
+  kplusgvecs_cart.clear();
+  minusModKplusG2.reserve(NumPlaneWaves);
+  //  RealType kcutoff2 = 2.0*ecut; //std::sqrt(2.0*ecut);
+  int ngIn = NumPlaneWaves;
+  for (int ig = 0, newig = 0; ig < ngIn; ig++)
+  {
+    //PosType tempvec = Lattice.k_cart(gvecCopy[ig]+twist);
+    PosType tempvec = gcartCopy[ig] + twist_cart;
+    RealType mod2   = dot(tempvec, tempvec);
+
+    // Keep all the g-vectors
+    // The cutoff energy is not stored in the HDF file now.
+    // Is truncating the gvectors to a spherical shell necessary?
+    if (true)
+    {
+      gvecs.push_back(gvecCopy[ig]);
+      kplusgvecs_cart.push_back(tempvec);
+      minusModKplusG2.push_back(-mod2);
+      //Remember which position in the HDF5 file this came from...for coefficients
+      inputmap[ig] = newig++;
+    }
+#if 0
+    if(mod2<=kcutoff2)
+    {
+      gvecs.push_back(gvecCopy[ig]);
+      kplusgvecs_cart.push_back(tempvec);
+      minusModKplusG2.push_back(-mod2);
+      //Remember which position in the HDF5 file this came from...for coefficients
+      inputmap[ig] = newig++;
+    }
+    else
+    {
+      inputmap[ig] = -1; //Temporary value...need to know final NumPlaneWaves.
+      NumPlaneWaves--;
+    }
+#endif
+  }
+#if defined(PWBasisT_USE_RECURSIVE)
+  //Store the maximum number of translations, within ecut, of any reciprocal cell vector.
+  for (int ig = 0; ig < NumPlaneWaves; ig++)
+    for (int i = 0; i < OHMMS_DIM; i++)
+      if (std::abs(gvecs[ig][i]) > maxg[i])
+        maxg[i] = std::abs(gvecs[ig][i]);
+  gvecs_shifted.resize(NumPlaneWaves);
+  for (int ig = 0; ig < NumPlaneWaves; ig++)
+    gvecs_shifted[ig] = gvecs[ig] + maxg;
+  maxmaxg = std::max(maxg[0], std::max(maxg[1], maxg[2]));
+  //changes the order???? ok
+  C.resize(3, 2 * maxmaxg + 2);
+#else
+  maxmaxg = 1;
+#endif
+  //    //make a copy of input to gvecCopy
+  ////    for(int ig=0, newig=0; ig<ngIn; ig++) {
+  //      //Check size of this g-vector
+  //      PosType tempvec = Lattice.k_cart(gvecCopy[ig]+twist);
+  //      RealType mod2 = dot(tempvec,tempvec);
+  //      if(mod2<=kcutoff2){ //Keep this element
+  //        gvecs.push_back(gvecCopy[ig]);
+  //        kplusgvecs_cart.push_back(tempvec);
+  //        minusModKplusG2.push_back(-mod2);
+  //        //Remember which position in the HDF5 file this came from...for coefficients
+  //        inputmap[ig] = newig++;
+  ////#if !defined(QMC_COMPLEX)
+  ////        //Build the negative vector. See comment at declaration (above) for details.
+  ////        if(gvecCopy[ig][0] < 0)
+  ////          negative.push_back(0);
+  ////        else if(gvecCopy[ig][0] > 0)
+  ////          negative.push_back(1);
+  ////        else { //gx == 0, test gy
+  ////          if(gvecCopy[ig][1] < 0)
+  ////            negative.push_back(0);
+  ////          else if(gvecCopy[ig][1] > 0)
+  ////            negative.push_back(1);
+  ////          else { //gx == gy == 0; test gz. If gz==0 also, take negative=1 (arbitrary)
+  ////            if(gvecCopy[ig][2] < 0)
+  ////              negative.push_back(0);
+  ////            else
+  ////              negative.push_back(1);
+  ////          }
+  ////        }
+  ////#endif
+  //      } else {
+  //        inputmap[ig] = -1; //Temporary value...need to know final NumPlaneWaves.
+  //        NumPlaneWaves--;
+  //      }
+  //    }
+  //Finalize the basis. Fix temporary values of inputmap.
+  //for(int ig=0; ig<inputmap.size(); ig++)
+  //  if(inputmap[ig] == -1)
+  //    inputmap[ig] = NumPlaneWaves; //For dumping coefficients of PWs>ecut
+  app_log() << "                       NumPlaneWaves (after)  =" << NumPlaneWaves << std::endl;
+}
+// template class PWBasisT<double>;
+// template class PWBasisT<float>;
+template class PWBasisT<std::complex<double>>;
+template class PWBasisT<std::complex<float>>;
+} // namespace qmcplusplus