Introduce USE_MULTIQUINTIC to use MultiQuinticSpline1D

Prepare paper comparing AoS-MO, SoA-MO and SoA-MO+muti

src/QMCWaveFunctions/lcao/AOBasisBuilder.h

@@ -222,8 +222,9 @@ COT* AOBasisBuilder<COT>::createAOSet(xmlNodePtr cur)
     ++it;
   }
 
+  radFuncBuilder.finalize();
   //aos->Rmax can be set small
-  aos->setRmax(0); 
+  //aos->setRmax(0); 
   aos->setBasisSetSize(-1);
   app_log() << "   Maximu Angular Momentum   = " << aos->Ylm.lmax() << std::endl
             << "   Number of Radial functors = " << aos->Rnl.size() << std::endl

src/QMCWaveFunctions/lcao/LCAOrbitalBuilder.h

@@ -18,8 +18,12 @@
 #ifndef QMCPLUSPLUS_LCAO_ORBITAL_BUILDER_H
 #define QMCPLUSPLUS_LCAO_ORBITAL_BUILDER_H
 
+//testing multiquintic
+//#define USE_MULTIQUINTIC
+
 #include "QMCWaveFunctions/BasisSetBase.h"
 #include "QMCWaveFunctions/lcao/NGFunctor.h"
+#include <QMCWaveFunctions/lcao/MultiQuinticSpline1D.h>
 #include "QMCWaveFunctions/lcao/SoaCartesianTensor.h"
 #include "QMCWaveFunctions/lcao/SoaSphericalTensor.h"
 #include "QMCWaveFunctions/lcao/SoaAtomicBasisSet.h"

src/QMCWaveFunctions/lcao/RadialOrbitalSetBuilder.h

@@ -31,6 +31,28 @@
 namespace qmcplusplus
 {
 
+  /** temporary function to compute the cutoff without constructing NGFunctor */
+template<typename Fin, typename T>
+inline double find_cutoff(Fin& in, T rmax)
+{
+  //myInFunc=in;
+  LogGrid<T> agrid;
+  const T delta=0.1;
+  const T eps=1e-5;
+  bool too_small=true;
+  agrid.set(eps,rmax,1001);
+  int i=1000;
+  T r=rmax;
+  while(too_small && i>0)
+  {
+    r=agrid[i--];
+    T x=in.f(r);
+    too_small=(std::abs(x)<eps);
+  }
+  return static_cast<double>(r);
+}
+
+
 /** Build a set of radial orbitals at the origin
  *
  * For a center,
@@ -63,7 +85,7 @@ public:
 
   ///constructor
   RadialOrbitalSetBuilder(xmlNodePtr cur=NULL);
-  ///destructor
+  ///destructor: cleanup gtoTemp, stoTemp
   ~RadialOrbitalSetBuilder();
 
   ///assign a CenteredOrbitalType to work on
@@ -83,11 +105,23 @@ public:
    */
   bool putCommon(xmlNodePtr cur);
 
+  /** This is when the radial orbitals are actually created */
+  void finalize();
+
 private:
+  //only check the cutoff
   void addGaussian(xmlNodePtr cur);
   void addSlater(xmlNodePtr cur);
   void addNumerical(xmlNodePtr cur, const std::string& dsname);
+
   hid_t m_fileid;
+
+  ///safe common cutoff radius
+  double m_rcut_safe;
+  ///store the temporary analytic data 
+  std::vector<GaussianCombo<RealType>*> gtoTemp;
+  std::vector<SlaterCombo<RealType>*> stoTemp;
+
 };
 
   template<typename COT>
@@ -209,6 +243,10 @@ private:
         //}
         //using 0.01 for the time being
       }
+
+      //set zero to use std::max
+      m_rcut_safe=0;
+
       return true;
     }
 
@@ -217,11 +255,6 @@ private:
    * \param nlms a vector containing the quantum numbers \f$(n,l,m,s)\f$
    * \return true is succeeds
    *
-   This function puts the STO on a logarithmic grid and calculates the boundary
-   conditions for the 1D Cubic Spline.  The derivates at the endpoint
-   are assumed to be all zero.  Note: for the radial orbital we use
-   \f[ f(r) = \frac{R(r)}{r^l}, \f] where \f$ R(r) \f$ is the usual
-   radial orbital and \f$ l \f$ is the angular momentum.
    */
   template<typename COT>
     bool
@@ -257,11 +290,13 @@ private:
       {
         addNumerical(cur,dsname);
       }
+#if !defined(USE_MULTIQUINTIC)
       if(lastRnl && m_orbitals->Rnl.size()> lastRnl)
       {
         app_log() << "\tSetting GridManager of " << lastRnl << " radial orbital to false" << std::endl;
         m_orbitals->Rnl[lastRnl]->setGridManager(false);
       }
+#endif
       return true;
     }
 
@@ -269,16 +304,69 @@ private:
   void RadialOrbitalSetBuilder<COT>::addGaussian(xmlNodePtr cur)
   {
     int L= m_nlms[1];
+#if defined(USE_MULTIQUINTIC)
+    GaussianCombo<RealType>* gset=new GaussianCombo<RealType>(L,Normalized);
+    gset->putBasisGroup(cur);
+    auto r0=find_cutoff(*gset,m_orbitals->Grids[0]->rmax());
+    m_rcut_safe=std::max(m_rcut_safe,r0);
+    //add this basisGroup
+    gtoTemp.push_back(gset);
+    m_orbitals->Rnl.push_back(nullptr); //CLEANUP
+    m_orbitals->RnlID.push_back(m_nlms);
+#else
     GaussianCombo<RealType> gset(L,Normalized);
     gset.putBasisGroup(cur);
     GridType* agrid = m_orbitals->Grids[0];
     RadialOrbitalType *radorb = new RadialOrbitalType(agrid);
-    if(m_rcut<0)
-      m_rcut = agrid->rmax();
+    if(m_rcut<0) m_rcut = agrid->rmax();
     Transform2GridFunctor<GaussianCombo<RealType>,RadialOrbitalType> transform(gset, *radorb);
     transform.generate(agrid->rmin(),m_rcut,agrid->size());
     m_orbitals->Rnl.push_back(radorb);
     m_orbitals->RnlID.push_back(m_nlms);
+#endif
+  }
+
+/* Finalize this set using the common grid
+ *
+ * This function puts the STO on a logarithmic grid and calculates the boundary
+ * conditions for the 1D Cubic Spline.  The derivates at the endpoint
+ * are assumed to be all zero.  Note: for the radial orbital we use
+ * \f[ f(r) = \frac{R(r)}{r^l}, \f] where \f$ R(r) \f$ is the usual
+ * radial orbital and \f$ l \f$ is the angular momentum.
+ */
+  template<typename COT>
+  void RadialOrbitalSetBuilder<COT>::finalize()
+  {
+    std::cout << "Going to finalize " << m_rcut_safe << std::endl;
+#if defined(USE_MULTIQUINTIC)
+    GridType* agrid = m_orbitals->Grids[0];
+    agrid->locate(static_cast<RealType>(m_rcut_safe));
+    agrid->set(agrid->rmin(),m_rcut_safe,agrid->Loc); 
+
+    m_orbitals->setRmax(m_rcut_safe);
+    int norbs=gtoTemp.size();
+
+    MultiQuinticSpline1D<RealType>* multiset=new MultiQuinticSpline1D<RealType>;;
+    int npts=agrid->size();
+    multiset->initialize(agrid,norbs);
+
+    if(gtoTemp.size())
+    {
+      RadialOrbitalType radorb(agrid);
+      for(int ib=0; ib<norbs; ++ib)
+      {
+        Transform2GridFunctor<GaussianCombo<RealType>,RadialOrbitalType> transform(*(gtoTemp[ib]), radorb);
+        transform.generate(agrid->rmin(),agrid->rmax(),agrid->size());
+        //m_orbitals->Rnl[ib]=radorb;
+        multiset->add_spline(ib,radorb.myFunc);
+        delete gtoTemp[ib];
+      }
+    }
+    m_orbitals->MultiRnl=multiset;
+#else
+    m_orbitals->setRmax(m_rcut);
+#endif
+
   }
 
   template<typename COT>

src/QMCWaveFunctions/lcao/SoaAtomicBasisSet.h

@@ -41,13 +41,16 @@ namespace qmcplusplus
       aligned_vector<int> LM;
       /**index of the corresponding radial orbital with quantum numbers \f$ (n,l) \f$ */
       aligned_vector<int> NL;
-      ///container for the radial orbitals
-      aligned_vector<ROT*> Rnl;
       ///container for the quantum-numbers
       std::vector<QuantumNumberType> RnlID;
-      ///set of grids
+      ///container for the radial orbitals
+      aligned_vector<ROT*> Rnl;
+#if defined(USE_MULTIQUINTIC)
+      ///Replace Rnl, Grids
+      MultiQuinticSpline1D<value_type> *MultiRnl;
+#endif
+      ///set of grids : keep this until completion
       std::vector<grid_type*> Grids;
-
       ///the constructor
       explicit SoaAtomicBasisSet(int lmax, bool addsignforM=false)
         :Ylm(lmax,addsignforM){}
@@ -58,6 +61,10 @@ namespace qmcplusplus
 
       SoaAtomicBasisSet<ROT,SH>* makeClone() const
       {
+#if defined(USE_MULTIQUINTIC)
+        SoaAtomicBasisSet<ROT,SH>* myclone=new SoaAtomicBasisSet<ROT,SH>(*this);
+        myclone->MultiRnl=MultiRnl->makeClone();
+#else
         grid_type *grid_clone=Grids[0]->makeClone();
         SoaAtomicBasisSet<ROT,SH>* myclone=new SoaAtomicBasisSet<ROT,SH>(*this);
         myclone->Grids[0]=grid_clone;
@@ -65,25 +72,26 @@ namespace qmcplusplus
         {
           myclone->Rnl[i]=new ROT(*Rnl[i],grid_clone,i==0);
         }
+#endif
         return myclone;
       }
 
       void checkInVariables(opt_variables_type& active)
       {
-        for(size_t nl=0; nl<Rnl.size(); nl++)
-          Rnl[nl]->checkInVariables(active);
+        //for(size_t nl=0; nl<Rnl.size(); nl++)
+        //  Rnl[nl]->checkInVariables(active);
       }
 
       void checkOutVariables(const opt_variables_type& active)
       {
-        for(size_t nl=0; nl<Rnl.size(); nl++)
-          Rnl[nl]->checkOutVariables(active);
+        //for(size_t nl=0; nl<Rnl.size(); nl++)
+        //  Rnl[nl]->checkOutVariables(active);
       }
 
       void resetParameters(const opt_variables_type& active)
       {
-        for(size_t nl=0; nl<Rnl.size(); nl++)
-          Rnl[nl]->resetParameters(active);
+        //for(size_t nl=0; nl<Rnl.size(); nl++)
+        //  Rnl[nl]->resetParameters(active);
       }
 
       /** return the number of basis functions
@@ -108,7 +116,11 @@ namespace qmcplusplus
       template<typename T>
       inline void setRmax(T rmax)
       {
-        Rmax=(rmax>0)?rmax:Grids[0]->rmax();
+#if defined(USE_MULTIQUINTIC)
+        Rmax=(rmax>0)? rmax: MultiRnl->rmax();
+#else
+        Rmax=(rmax>0)? rmax:Grids[0]->rmax();
+#endif
       }
 
       /** reset the target ParticleSet
@@ -139,13 +151,19 @@ namespace qmcplusplus
           const T x=-dr[0], y=-dr[1], z=-dr[2];
           Ylm.evaluateVGL(x,y,z);
 
-          const size_t nl_max=Rnl.size();
+          const size_t nl_max=RnlID.size();
           T phi[nl_max]; 
           T dphi[nl_max];
           T d2phi[nl_max];
 
+#if defined(USE_MULTIQUINTIC)
+          MultiRnl->evaluate(r,phi,dphi,d2phi);
+#else
           for(size_t nl=0; nl<nl_max; ++nl)
+          {
             phi[nl]=Rnl[nl]->evaluate(r,dphi[nl],d2phi[nl]);
+          }
+#endif
 
           //V,Gx,Gy,Gz,L
           T* restrict psi   =vgl.data(0)+offset; const T* restrict ylm_v=Ylm[0]; //value
@@ -188,10 +206,14 @@ namespace qmcplusplus
           T ylm_v[Ylm.size()]; 
           Ylm.evaluateV(-dr[0],-dr[1],-dr[2],ylm_v);
 
-          const int nl_max=Rnl.size();
+          const int nl_max=RnlID.size();
           T phi_r[nl_max];
+#if defined(USE_MULTIQUINTIC)
+          MultiRnl->evaluate(r,phi_r);
+#else
           for(size_t nl=0; nl<nl_max; ++nl)
             phi_r[nl]=Rnl[nl]->evaluate(r);
+#endif
 
           for(size_t ib=0; ib<BasisSetSize; ++ib)
           {