Ewald cleanup

src/Particle/Lattice/CrystalLattice.h

@@ -241,6 +241,7 @@ struct CrystalLattice : public LRBreakupParameters<T, D>
     Base::LR_kc         = rhs.LR_kc;
     Base::LR_rc         = rhs.LR_rc;
     Base::ewaldAlpha    = rhs.ewaldAlpha;
+    Base::ewaldGrid     = rhs.ewaldGrid;
     Base::nlat          = rhs.nlat;
     Base::ndim          = rhs.ndim;
     Base::dgate         = rhs.dgate;

src/Particle/Lattice/LRBreakupParameters.h

@@ -33,6 +33,7 @@ class LRBreakupParameters<T, 3>
   T LR_tol;
   T ewaldAlpha;
   unsigned nlat;
+  unsigned ewaldGrid;
   ///number of strictly enforced periodic spatial dimensions
   /// ewald_strict2d sets ndim=2, otherwise ndim=3
   unsigned ndim;
@@ -42,7 +43,7 @@ class LRBreakupParameters<T, 3>
 
   ///default constructor
   LRBreakupParameters() : LR_dim_cutoff(15.0), LR_rc(1e6), LR_kc(0.0), LR_tol(3e-4),
-                          ewaldAlpha(-1.0), nlat(0), ndim(3), dgate(-1.0), mimg(100) {}
+                          ewaldAlpha(-1.0), nlat(0), ewaldGrid(1001), ndim(3), dgate(-1.0), mimg(100) {}
 
   ///Set LR_rc = radius of smallest sphere inside box and kc=dim/rc
   void SetLRCutoffs(const TinyVector<TinyVector<T, 3>, 3>& a)

src/Particle/LongRange/EwaldHandler2D.cpp

@@ -22,7 +22,12 @@ EwaldHandler2D::EwaldHandler2D(ParticleSet& ref, mRealType kc_in)
   LR_rc = ref.getLattice().LR_rc; // CoulombPBC needs get_rc() to createSpline4RbyVs
   LR_kc = ref.getLattice().LR_kc; // get_kc() is used in QMCFiniteSize
   alpha = ref.getLattice().ewaldAlpha;
-  if (alpha < 0) alpha = std::sqrt(LR_kc/2.0/LR_rc);
+  if (alpha < 0)
+  {
+    alpha = std::sqrt(LR_kc/2.0/LR_rc);
+    app_log() << "init alpha = " << alpha << "\n";
+    alpha = guess_ewald_alpha(LR_rc);
+  }
   area = ref.getLattice().Volume/ref.getLattice().R(2,2);
   // report
   app_log() << "    alpha = " << alpha << " area = " << area << std::endl;

src/Particle/LongRange/EwaldHandler2D.h

@@ -56,7 +56,6 @@ class EwaldHandler2D : public LRHandlerBase
   void resetTargetParticleSet(ParticleSet& ref) override {}
   // overrides end
 private:
-  mRealType alpha;
   mRealType area;
 };
 } // qmcplusplus

src/Particle/LongRange/EwaldHandlerQuasi2D.cpp

@@ -25,7 +25,12 @@ EwaldHandlerQuasi2D::EwaldHandlerQuasi2D(ParticleSet& ref, mRealType kc_in)
   LR_rc = ref.getLattice().LR_rc; // CoulombPBC needs get_rc() to createSpline4RbyVs
   LR_kc = ref.getLattice().LR_kc; // get_kc() is used in QMCFiniteSize
   alpha = ref.getLattice().ewaldAlpha;
-  if (alpha < 0) alpha = std::sqrt(LR_kc/2.0/LR_rc);
+  if (alpha < 0)
+  {
+    alpha = std::sqrt(LR_kc/2.0/LR_rc);
+    app_log() << "init alpha = " << alpha << "\n";
+    alpha = guess_ewald_alpha(LR_rc);
+  }
   area = ref.getLattice().Volume/ref.getLattice().R(2,2);
   // report
   app_log() << "    alpha = " << alpha << " area = " << area << std::endl;

src/Particle/LongRange/EwaldHandlerQuasi2D.h

@@ -64,7 +64,6 @@ class EwaldHandlerQuasi2D : public LRHandlerBase
   void resetTargetParticleSet(ParticleSet& ref) override {}
   // overrides end
 private:
-  mRealType alpha;
   mRealType area;
   ///store |k|
   std::vector<mRealType> kmags;

src/Particle/LongRange/EwaldScreen2D.cpp

@@ -28,6 +28,8 @@ EwaldScreen2D::EwaldScreen2D(ParticleSet& ref, mRealType kc_in)
   {
     const mRealType alpha0 = std::sqrt(LR_kc/2.0/LR_rc);
     alpha = std::max(0.61093226575644/dgate, alpha0);
+    app_log() << "init alpha = " << alpha << "\n";
+    alpha = guess_ewald_alpha(LR_rc);
   }
   area = ref.getLattice().Volume/ref.getLattice().R(2,2);
   // report

src/Particle/LongRange/EwaldScreen2D.h

@@ -51,7 +51,6 @@ class EwaldScreen2D : public LRHandlerBase
   void resetTargetParticleSet(ParticleSet& ref) override {}
   // overrides end
 private:
-  mRealType alpha;
   mRealType area;
   mRealType dgate;
 };

src/Particle/LongRange/LRCoulombSingleton.cpp

@@ -153,13 +153,6 @@ std::unique_ptr<OneDimCubicSpline<T>> createSpline4RbyVs_temp(const LRHandlerBas
                                                               const T vtol)
 {
   using func_type = OneDimCubicSpline<T>;
-  std::unique_ptr<LinearGrid<T>> agrid_local;
-  if (agrid == nullptr)
-  {
-    agrid_local = std::make_unique<LinearGrid<T>>();
-    agrid_local->set(0.0, rcut, 1001);
-    agrid = agrid_local.get();
-  }
   const int ng = agrid->size();
   std::vector<T> v(ng);
   // check last grid point is close to zero
@@ -202,6 +195,7 @@ std::unique_ptr<OneDimCubicSpline<T>> createSpline4RbyVs_temp(const LRHandlerBas
     std::ostringstream msg;
     msg << "LRCoulombSingleton::createSpline4RbyVs_temp. Short-range potential\n";
     msg << "  interpolation chi^2 = " << chi2 << "\n";
+    msg << "  increase ewald_grid. \n";
     throw std::runtime_error(msg.str());
   }
   return V0;
@@ -213,13 +207,6 @@ std::unique_ptr<OneDimCubicSpline<T>> createSpline4RbyVsDeriv_temp(const LRHandl
                                                                    const LinearGrid<T>* agrid)
 {
   using func_type = OneDimCubicSpline<T>;
-  std::unique_ptr<LinearGrid<T>> agrid_local;
-  if (agrid == nullptr)
-  {
-    agrid_local = std::make_unique<LinearGrid<T>>();
-    agrid_local->set(0.0, rcut, 1001);
-    agrid = agrid_local.get();
-  }
   int ng = agrid->size();
   std::vector<T> v(ng);
   T r = (*agrid)[0];

src/Particle/LongRange/LRHandlerBase.h

@@ -37,6 +37,8 @@ struct LRHandlerBase
   mRealType LR_kc;
   /// Maximum r cutoff
   mRealType LR_rc;
+  /// Ewald screening parameter
+  mRealType alpha;
   ///evaluate long-range piece (Fk is not all zero)
   bool llr;
   ///Fourier component for all the k-point
@@ -62,7 +64,7 @@ struct LRHandlerBase
 
 
   //constructor
-  explicit LRHandlerBase(mRealType kc) : MaxKshell(0), LR_kc(kc), LR_rc(0), llr(true), ClassName("LRHandlerBase") {}
+  explicit LRHandlerBase(mRealType kc) : MaxKshell(0), LR_kc(kc), LR_rc(0), alpha(LR_kc), llr(true), ClassName("LRHandlerBase") {}
 
   // virtual destructor
   virtual ~LRHandlerBase() = default;
@@ -72,6 +74,52 @@ struct LRHandlerBase
   //return k cutoff
   inline mRealType get_kc() const { return LR_kc; }
 
+  inline mRealType guess_ewald_alpha(mRealType rcut)
+  {
+    const int miter = 1000;
+    const mRealType tol = 1e-10;
+    const mRealType edge = 0.9;
+    mRealType alpha_max = 1e6/rcut;
+    mRealType alpha_min = 1e-6/rcut;
+    mRealType vr = evaluate(rcut, 1./rcut);
+    mRealType vr_edge = evaluate(edge*rcut, 1./rcut/edge);
+    bool ltoo_small = vr > tol;
+    bool ltoo_large = vr_edge < tol;
+    int niter = 0;
+    app_log() << "guessing alpha" << std::endl;
+    while ((niter < miter) and (ltoo_small or ltoo_large))
+    {
+      app_log() << niter << " amin = " << alpha_min << " a = " << alpha << " amax = " << alpha_max << "\n";
+      if (alpha_max-alpha_min < 2*tol) break;
+      // adjust
+      if (ltoo_small)
+      {
+        alpha = (alpha+alpha_max)/2.0;
+      }
+      else if (ltoo_large)
+      {
+        alpha = (alpha+alpha_min)/2.0;
+      }
+      else
+      {
+        break;
+      }
+      // re-evaluate
+      vr = evaluate(rcut, 1./rcut);
+      vr_edge = evaluate(rcut/edge, 1./rcut/edge);
+      ltoo_small = vr > tol;
+      ltoo_large = vr_edge < tol;
+      if (ltoo_small)
+        alpha_min = alpha;
+      else if (ltoo_large)
+        alpha_max = alpha;
+      else
+        break;
+      niter++;
+    }
+    return alpha;
+  }
+
   inline mRealType evaluate_w_sk(const std::vector<int>& kshell, const pRealType* restrict sk) const
   {
     mRealType vk = 0.0;

src/Particle/LongRange/Screen2D.cpp

@@ -24,9 +24,32 @@ Screen2D::Screen2D(ParticleSet& ref, mRealType kc_in)
   LR_rc = ref.getLattice().LR_rc; // CoulombPBC needs get_rc() to createSpline4RbyVs
   LR_kc = ref.getLattice().LR_kc; // get_kc() is used in QMCFiniteSize
   area = ref.getLattice().Volume/ref.getLattice().R(2,2);
+  int nlat = ref.getLattice().nlat;
+  mRealType reach = LR_rc + nlat*2.0*LR_rc;
+  // find necessary rcut
+  mRealType new_rcut = findRcut();
+  LR_rc = std::max(new_rcut, LR_rc);
   // report
   app_log() << "    dgate = " << dgate << std::endl;
   app_log() << "    mimg  = " << mimg << std::endl;
+  app_log() << "     rcut = " << LR_rc << std::endl;
+  app_log() << "     nlat = " << nlat << std::endl;
+  app_log() << "    reach = " << reach << std::endl;
+  bool ltoo_small = (new_rcut > reach);
+  if (ltoo_small)
+  {
+    const int mlat = 10;
+    while (nlat < mlat)
+    {
+      nlat++;
+      reach = LR_rc + nlat*2.0*LR_rc;
+      if (reach >= new_rcut) break;
+    }
+    std::ostringstream msg;
+    msg << "nlat = " << ref.getLattice().nlat << " is too small"
+        << "  try increase to " << nlat << "\n";
+    throw std::runtime_error(msg.str());
+  }
   fillFk(ref.getSimulationCell().getKLists());
 }
 
@@ -66,4 +89,39 @@ Screen2D::mRealType Screen2D::evaluate(mRealType r, mRealType rinv) const
   return pre*vsr;
 }
 
+Screen2D::mRealType Screen2D::findRcut()
+{
+  const int miter = 1000;
+  const mRealType edge = 0.9;
+  const mRealType tol = 1e-10;
+  mRealType r = LR_rc;
+  mRealType rmax = 10*LR_rc;
+  mRealType rmin = 0.0;
+  int iter = 0;
+  mRealType vtail = evaluate(r, 1./r);
+  mRealType vedge = evaluate(edge*r, 1./r/edge);
+  app_log() << "  Screen2D::findRcut\n";
+  while ((iter < miter) and ((vedge<tol) or (vtail>tol)))
+  {
+    app_log() << iter << " " << rmin << " " << r << " " << rmax << "\n";
+    if ((rmax-rmin) < 2*tol) break;
+    if (vtail>tol)
+      r = (r+rmax)/2.0;
+    else if (vedge<tol)
+      r = (r+rmin)/2.0;
+    else
+      break;
+    vtail = evaluate(r, 1./r);
+    vedge = evaluate(edge*r, 1./r/edge);
+    if (vtail>tol)
+      rmin = r;
+    else if (vedge<tol)
+      rmax = r;
+    else
+      break;
+    iter++;
+  }
+  return r;
+}
+
 } // qmcplusplus

src/Particle/LongRange/Screen2D.h

@@ -27,6 +27,7 @@ class Screen2D : public LRHandlerBase
   // short-range part
   mRealType evaluate(mRealType r, mRealType rinv) const override;
   inline mRealType evaluateLR_r0() const override {return 0.0;}
+  mRealType findRcut();
 
   // long-range part
   inline mRealType evaluateLR(mRealType r) const override {return 0.0;}

src/Particle/ParticleIO/LatticeIO.cpp

@@ -120,6 +120,10 @@ bool LatticeParser::put(xmlNodePtr cur)
       {
         putContent(ref_.ewaldAlpha, cur);
       }
+      else if (aname == "ewald_grid")
+      {
+        putContent(ref_.ewaldGrid, cur);
+      }
       else if (aname == "nlat")
       {
         putContent(ref_.nlat, cur);

src/QMCHamiltonians/CoulombPBCAA.cpp

@@ -495,10 +495,18 @@ void CoulombPBCAA::initBreakup(ParticleSet& P)
   //AA->initBreakup(*PtclRef);
   myConst = evalConsts();
   myRcut  = AA->get_rc();
-  myRcut  = myRcut + 2*myRcut*nlat;
+
+  if (myGrid == nullptr)
+  {
+    myGrid = std::make_shared<LinearGrid<RealType>>();
+    int ng = P.getLattice().ewaldGrid;
+    app_log() << "    CoulombPBCAB::initBreakup\n  Setting a linear grid=[0," << myRcut << ") number of grid points =" << ng
+              << std::endl;
+    myGrid->set(0, myRcut, ng);
+  }
 
   if (rVs == nullptr)
-    rVs = LRCoulombSingleton::createSpline4RbyVs(AA.get(), myRcut);
+    rVs = LRCoulombSingleton::createSpline4RbyVs(AA.get(), myRcut, myGrid.get());
 
   rVs_offload = std::make_shared<const OffloadSpline>(*rVs);
 
@@ -507,7 +515,7 @@ void CoulombPBCAA::initBreakup(ParticleSet& P)
     dAA = LRCoulombSingleton::getDerivHandler(P);
     if (rVsforce == nullptr)
     {
-      rVsforce = LRCoulombSingleton::createSpline4RbyVs(dAA.get(), myRcut);
+      rVsforce = LRCoulombSingleton::createSpline4RbyVs(dAA.get(), myRcut, myGrid.get());
     }
   }
 }

src/QMCHamiltonians/CoulombPBCAA.h

@@ -62,7 +62,10 @@ struct CoulombPBCAA : public OperatorBase, public ForceBase
   int MemberAttribIndx;
   int NumCenters;
   Return_t myConst;
+  ///cutoff radius of the short-range part
   RealType myRcut;
+  ///radial grid
+  std::shared_ptr<GridType> myGrid;
   std::string PtclRefName;
 
   std::vector<RealType> Zat, Zspec;