Implement evaluateDerivRatios in J1 and J2

src/Particle/DistanceTable.h

@@ -197,23 +197,6 @@ class DistanceTable
       dt_list[iw].finalizePbyP(p_list[iw]);
   }
 
-  /** build a compact list of a neighbor for the iat source
-   * @param iat source particle id
-   * @param rcut cutoff radius
-   * @param jid compressed index
-   * @param dist compressed distance
-   * @param displ compressed displacement
-   * @return number of target particles within rcut
-   */
-  virtual size_t get_neighbors(int iat,
-                               RealType rcut,
-                               int* restrict jid,
-                               RealType* restrict dist,
-                               PosType* restrict displ) const
-  {
-    return 0;
-  }
-
   /** find the first nearest neighbor
    * @param iat source particle id
    * @param r distance

src/Particle/SoaDistanceTableAB.h

@@ -102,28 +102,6 @@ struct SoaDistanceTableAB : public DTD_BConds<T, D, SC>, public DistanceTableAB
       std::copy_n(temp_dr_.data(idim), num_sources_, displacements_[iat].data(idim));
   }
 
-  size_t get_neighbors(int iat,
-                       RealType rcut,
-                       int* restrict jid,
-                       RealType* restrict dist,
-                       PosType* restrict displ) const override
-  {
-    constexpr T cminus(-1);
-    size_t nn = 0;
-    for (int jat = 0; jat < num_targets_; ++jat)
-    {
-      const RealType rij = distances_[jat][iat];
-      if (rij < rcut)
-      { //make the compact list
-        jid[nn]   = jat;
-        dist[nn]  = rij;
-        displ[nn] = cminus * displacements_[jat][iat];
-        nn++;
-      }
-    }
-    return nn;
-  }
-
   int get_first_neighbor(IndexType iat, RealType& r, PosType& dr, bool newpos) const override
   {
     RealType min_dist = std::numeric_limits<RealType>::max();
@@ -160,21 +138,6 @@ struct SoaDistanceTableAB : public DTD_BConds<T, D, SC>, public DistanceTableAB
     return index;
   }
 
-  size_t get_neighbors(int iat, RealType rcut, RealType* restrict dist) const
-  {
-    size_t nn = 0;
-    for (int jat = 0; jat < num_targets_; ++jat)
-    {
-      const RealType rij = distances_[jat][iat];
-      if (rij < rcut)
-      { //make the compact list
-        dist[nn] = rij;
-        nn++;
-      }
-    }
-    return nn;
-  }
-
 private:
   /// timer for evaluate()
   NewTimer& evaluate_timer_;

src/Particle/SoaDistanceTableABOMPTarget.h

@@ -377,28 +377,6 @@ class SoaDistanceTableABOMPTarget : public DTD_BConds<T, D, SC>, public Distance
       std::copy_n(temp_dr_.data(idim), num_sources_, displacements_[iat].data(idim));
   }
 
-  size_t get_neighbors(int iat,
-                       RealType rcut,
-                       int* restrict jid,
-                       RealType* restrict dist,
-                       PosType* restrict displ) const override
-  {
-    constexpr T cminus(-1);
-    size_t nn = 0;
-    for (int jat = 0; jat < num_targets_; ++jat)
-    {
-      const RealType rij = distances_[jat][iat];
-      if (rij < rcut)
-      { //make the compact list
-        jid[nn]   = jat;
-        dist[nn]  = rij;
-        displ[nn] = cminus * displacements_[jat][iat];
-        nn++;
-      }
-    }
-    return nn;
-  }
-
   int get_first_neighbor(IndexType iat, RealType& r, PosType& dr, bool newpos) const override
   {
     RealType min_dist = std::numeric_limits<RealType>::max();
@@ -435,21 +413,6 @@ class SoaDistanceTableABOMPTarget : public DTD_BConds<T, D, SC>, public Distance
     return index;
   }
 
-  size_t get_neighbors(int iat, RealType rcut, RealType* restrict dist) const
-  {
-    size_t nn = 0;
-    for (int jat = 0; jat < num_targets_; ++jat)
-    {
-      const RealType rij = distances_[jat][iat];
-      if (rij < rcut)
-      { //make the compact list
-        dist[nn] = rij;
-        nn++;
-      }
-    }
-    return nn;
-  }
-
 private:
   /// timer for offload portion
   NewTimer& offload_timer_;

src/QMCHamiltonians/NonLocalECPComponent.cpp

@@ -75,9 +75,7 @@ void NonLocalECPComponent::resize_warrays(int n, int m, int l)
   vrad.resize(m);
   dvrad.resize(m);
   vgrad.resize(m);
-  wvec.resize(m);
-  Amat.resize(n * m);
-  dAmat.resize(n * m);
+  wvec.resize(n);
   lpol.resize(l + 1, 1.0);
   dlpol.resize(l + 1, 0.0);
   rrotsgrid_m.resize(n);

src/QMCHamiltonians/NonLocalECPComponent.h

@@ -61,7 +61,7 @@ class NonLocalECPComponent : public QMCTraits
   ///weight of the spherical grid
   std::vector<RealType> sgridweight_m;
   ///Working arrays
-  std::vector<ValueType> wvec, Amat, dAmat;
+  std::vector<ValueType> wvec;
 
   //Position delta for virtual moves.
   std::vector<PosType> deltaV;

src/QMCHamiltonians/NonLocalECPotential.deriv.cpp

@@ -63,34 +63,43 @@ NonLocalECPComponent::RealType NonLocalECPComponent::evaluateValueAndDerivatives
                                                                                  const std::vector<ValueType>& dlogpsi,
                                                                                  std::vector<ValueType>& dhpsioverpsi)
 {
-  dratio.resize(optvars.num_active_vars, nknot);
+  const size_t num_vars = optvars.num_active_vars;
+  dratio.resize(nknot, num_vars);
   dlogpsi_vp.resize(dlogpsi.size());
 
-  ValueType pairpot;
-  ParticleSet::ParticlePos deltarV(nknot);
+  deltaV.resize(nknot);
 
   //displacements wrt W.R[iel]
   for (int j = 0; j < nknot; j++)
-    deltarV[j] = r * rrotsgrid_m[j] - dr;
+    deltaV[j] = r * rrotsgrid_m[j] - dr;
 
-  for (int j = 0; j < nknot; j++)
+  if (VP)
   {
-    PosType pos_now = W.R[iel];
-    W.makeMove(iel, deltarV[j]);
-    psiratio[j] = psi.calcRatio(W, iel);
-    psi.acceptMove(W, iel);
-    W.acceptMove(iel);
-
-    //use existing methods
-    std::fill(dlogpsi_vp.begin(), dlogpsi_vp.end(), 0.0);
-    psi.evaluateDerivativesWF(W, optvars, dlogpsi_vp);
-    for (int v = 0; v < dlogpsi_vp.size(); ++v)
-      dratio(v, j) = (dlogpsi_vp[v] - dlogpsi[v]);
-
-    W.makeMove(iel, -deltarV[j]);
-    psi.calcRatio(W, iel);
-    psi.acceptMove(W, iel);
-    W.acceptMove(iel);
+    // Compute ratios with VP
+    VP->makeMoves(iel, W.R[iel], deltaV, true, iat);
+    psi.evaluateDerivRatios(*VP, optvars, psiratio, dratio);
+  }
+  else
+  {
+    for (int j = 0; j < nknot; j++)
+    {
+      PosType pos_now = W.R[iel];
+      W.makeMove(iel, deltaV[j]);
+      psiratio[j] = psi.calcRatio(W, iel);
+      psi.acceptMove(W, iel);
+      W.acceptMove(iel);
+
+      //use existing methods
+      std::fill(dlogpsi_vp.begin(), dlogpsi_vp.end(), 0.0);
+      psi.evaluateDerivativesWF(W, optvars, dlogpsi_vp);
+      for (int v = 0; v < dlogpsi_vp.size(); ++v)
+        dratio(j, v) = dlogpsi_vp[v] - dlogpsi[v];
+
+      W.makeMove(iel, -deltaV[j]);
+      psi.calcRatio(W, iel);
+      psi.acceptMove(W, iel);
+      W.acceptMove(iel);
+    }
   }
 
   for (int j = 0; j < nknot; ++j)
@@ -99,6 +108,7 @@ NonLocalECPComponent::RealType NonLocalECPComponent::evaluateValueAndDerivatives
   for (int ip = 0; ip < nchannel; ip++)
     vrad[ip] = nlpp_m[ip]->splint(r) * wgt_angpp_m[ip];
 
+  RealType pairpot(0);
   const RealType rinv = RealType(1) / r;
   // Compute spherical harmonics on grid
   for (int j = 0, jl = 0; j < nknot; j++)
@@ -112,33 +122,18 @@ NonLocalECPComponent::RealType NonLocalECPComponent::evaluateValueAndDerivatives
       lpol[l + 1] = (Lfactor1[l] * zz * lpol[l] - l * lpolprev) * Lfactor2[l];
       lpolprev    = lpol[l];
     }
-    for (int l = 0; l < nchannel; l++, jl++)
-      Amat[jl] = lpol[angpp_m[l]];
-  }
-  if (nchannel == 1)
-  {
-    pairpot = vrad[0] * BLAS::dot(nknot, &Amat[0], &psiratio[0]);
-    for (int v = 0; v < dhpsioverpsi.size(); ++v)
-    {
-      for (int j = 0; j < nknot; ++j)
-        dratio(v, j) = psiratio[j] * dratio(v, j);
-      dhpsioverpsi[v] += vrad[0] * BLAS::dot(nknot, &Amat[0], dratio[v]);
-    }
-  }
-  else
-  {
-    BLAS::gemv(nknot, nchannel, &Amat[0], &psiratio[0], &wvec[0]);
-    pairpot = BLAS::dot(nchannel, &vrad[0], &wvec[0]);
-    for (int v = 0; v < dhpsioverpsi.size(); ++v)
-    {
-      for (int j = 0; j < nknot; ++j)
-        dratio(v, j) = psiratio[j] * dratio(v, j);
-      BLAS::gemv(nknot, nchannel, &Amat[0], dratio[v], &wvec[0]);
-      dhpsioverpsi[v] += BLAS::dot(nchannel, &vrad[0], &wvec[0]);
-    }
+
+    RealType lsum = 0.0;
+    for (int l = 0; l < nchannel; l++)
+      lsum += vrad[l] * lpol[angpp_m[l]];
+
+    wvec[j] = lsum * psiratio[j];
+    pairpot += std::real(wvec[j]);
   }
 
-  return std::real(pairpot);
+  BLAS::gemv('N', num_vars, nknot, 1.0, dratio.data(), num_vars, wvec.data(), 1, 1.0, dhpsioverpsi.data(), 1);
+
+  return pairpot;
 }
 
 } // namespace qmcplusplus

src/QMCHamiltonians/tests/test_ecp.cpp

@@ -315,10 +315,10 @@ TEST_CASE("Evaluate_ecp", "[hamiltonian]")
   };
 
   {
-    test_evaluateOne();
     nlpp->initVirtualParticle(elec);
     test_evaluateOne();
     nlpp->deleteVirtualParticle();
+    test_evaluateOne();
   }
 
   opt_variables_type optvars;
@@ -367,10 +367,10 @@ TEST_CASE("Evaluate_ecp", "[hamiltonian]")
   };
 
   {
-    test_evaluateValueAndDerivatives();
     nlpp->initVirtualParticle(elec);
     test_evaluateValueAndDerivatives();
     nlpp->deleteVirtualParticle();
+    test_evaluateValueAndDerivatives();
   }
 
   double Value2(0.0);

src/QMCWaveFunctions/Jastrow/BsplineFunctor.h

@@ -674,7 +674,7 @@ struct BsplineFunctor : public OptimizableFunctorBase
     return true;
   }
 
-  inline bool evaluateDerivatives(real_type r, std::vector<real_type>& derivs)
+  inline bool evaluateDerivatives(real_type r, std::vector<real_type>& derivs) override
   {
     if (r >= cutoff_radius)
       return false;

src/QMCWaveFunctions/Jastrow/FakeFunctor.h

@@ -125,10 +125,8 @@ struct FakeFunctor : public OptimizableFunctorBase
   }
 
 
-  // inline bool evaluateDerivatives(real_type r, std::vector<real_type>& derivs)
-
   /// compute derivatives with respect to variational parameters
-  inline bool evaluateDerivatives(real_type r, std::vector<real_type>& derivs) { return true; }
+  inline bool evaluateDerivatives(real_type r, std::vector<real_type>& derivs) override { return true; }
 
   bool put(xmlNodePtr cur) override { return true; }