mw_ APIs for SOECPotential and SOECPComponent

src/Particle/VirtualParticleSet.cpp

@@ -37,7 +37,8 @@ struct VPMultiWalkerMem : public Resource
   Resource* makeClone() const override { return new VPMultiWalkerMem(*this); }
 };
 
-VirtualParticleSet::VirtualParticleSet(const ParticleSet& p, int nptcl, size_t dt_count_limit) : ParticleSet(p.getSimulationCell())
+VirtualParticleSet::VirtualParticleSet(const ParticleSet& p, int nptcl, size_t dt_count_limit)
+    : ParticleSet(p.getSimulationCell())
 {
   setName("virtual");
 
@@ -194,4 +195,54 @@ void VirtualParticleSet::mw_makeMoves(const RefVectorWithLeader<VirtualParticleS
   ParticleSet::mw_update(p_list);
 }
 
+void VirtualParticleSet::mw_makeMovesWithSpin(const RefVectorWithLeader<VirtualParticleSet>& vp_list,
+                                              const RefVectorWithLeader<ParticleSet>& refp_list,
+                                              const RefVector<const std::vector<PosType>>& deltaV_list,
+                                              const RefVector<const std::vector<RealType>>& deltaS_list,
+                                              const RefVector<const NLPPJob<RealType>>& joblist,
+                                              bool sphere)
+{
+  auto& vp_leader = vp_list.getLeader();
+  if (!vp_leader.isSpinor())
+    throw std::runtime_error(
+        "VirtualParticleSet::mw_makeMovesWithSpin should not be called if particle sets aren't spionor types");
+  vp_leader.onSphere = sphere;
+  vp_leader.refPS    = refp_list.getLeader();
+
+  const size_t nVPs = countVPs(vp_list);
+  auto& mw_refPctls = vp_leader.getMultiWalkerRefPctls();
+  mw_refPctls.resize(nVPs);
+
+  RefVectorWithLeader<ParticleSet> p_list(vp_leader);
+  p_list.reserve(vp_list.size());
+
+  size_t ivp = 0;
+  for (int iw = 0; iw < vp_list.size(); iw++)
+  {
+    VirtualParticleSet& vp(vp_list[iw]);
+    const std::vector<PosType>& deltaV(deltaV_list[iw]);
+    const std::vector<RealType>& deltaS(deltaS_list[iw]);
+    const NLPPJob<RealType>& job(joblist[iw]);
+
+    vp.onSphere      = sphere;
+    vp.refPS         = refp_list[iw];
+    vp.refPtcl       = job.electron_id;
+    vp.refSourcePtcl = job.ion_id;
+    assert(vp.R.size() == deltaV.size());
+    assert(vp.spins.size() == deltaS.size());
+    assert(vp.R.size() == vp.spins.size());
+    for (size_t k = 0; k < vp.R.size(); k++, ivp++)
+    {
+      vp.R[k]          = refp_list[iw].R[vp.refPtcl] + deltaV[k];
+      vp.spins[k]      = refp_list[iw].spins[vp.refPtcl] + deltaS[k];
+      mw_refPctls[ivp] = vp.refPtcl;
+    }
+    p_list.push_back(vp);
+  }
+  assert(ivp == nVPs);
+
+  mw_refPctls.updateTo();
+  ParticleSet::mw_update(p_list);
+}
+
 } // namespace qmcplusplus

src/Particle/VirtualParticleSet.h

@@ -88,7 +88,11 @@ class VirtualParticleSet : public ParticleSet
      * @param sphere set true if VP are on a sphere around the reference source particle
      * @param iat reference source particle
      */
-  void makeMoves(const ParticleSet& refp, int jel, const std::vector<PosType>& deltaV, bool sphere = false, int iat = -1);
+  void makeMoves(const ParticleSet& refp,
+                 int jel,
+                 const std::vector<PosType>& deltaV,
+                 bool sphere = false,
+                 int iat     = -1);
 
   /** move virtual particles to new postions and update distance tables
      * @param refp reference particle set
@@ -111,6 +115,13 @@ class VirtualParticleSet : public ParticleSet
                            const RefVector<const NLPPJob<RealType>>& joblist,
                            bool sphere);
 
+  static void mw_makeMovesWithSpin(const RefVectorWithLeader<VirtualParticleSet>& vp_list,
+                                   const RefVectorWithLeader<ParticleSet>& p_list,
+                                   const RefVector<const std::vector<PosType>>& deltaV_list,
+                                   const RefVector<const std::vector<RealType>>& deltaS_list,
+                                   const RefVector<const NLPPJob<RealType>>& joblist,
+                                   bool sphere);
+
   static RefVectorWithLeader<ParticleSet> RefVectorWithLeaderParticleSet(
       const RefVectorWithLeader<VirtualParticleSet>& vp_list)
   {

src/QMCHamiltonians/SOECPComponent.cpp

@@ -15,6 +15,7 @@
 #include "SOECPComponent.h"
 #include "Numerics/Ylm.h"
 #include "CPU/BLAS.hpp"
+#include "NLPPJob.h"
 
 namespace qmcplusplus
 {
@@ -141,10 +142,6 @@ SOECPComponent::RealType SOECPComponent::evaluateOne(ParticleSet& W,
                                                      RealType r,
                                                      const PosType& dr)
 {
-
-  for (int ip = 0; ip < nchannel_; ip++)
-    vrad_[ip] = sopp_m_[ip]->splint(r);
-
   RealType sold = W.spins[iel];
   buildTotalQuadrature(r, dr, sold);
 
@@ -162,6 +159,11 @@ SOECPComponent::RealType SOECPComponent::evaluateOne(ParticleSet& W,
       Psi.resetPhaseDiff();
     }
 
+  return calculateProjector(r, dr, sold);
+}
+
+SOECPComponent::RealType SOECPComponent::calculateProjector(RealType r, const PosType& dr, RealType sold)
+{
   ComplexType pairpot;
   for (int iq = 0; iq < total_knots_; iq++)
   {
@@ -190,6 +192,81 @@ SOECPComponent::RealType SOECPComponent::evaluateOne(ParticleSet& W,
   return std::real(pairpot);
 }
 
+void SOECPComponent::mw_evaluateOne(const RefVectorWithLeader<SOECPComponent>& soecp_component_list,
+                                    const RefVectorWithLeader<ParticleSet>& p_list,
+                                    const RefVectorWithLeader<TrialWaveFunction>& psi_list,
+                                    const RefVector<const NLPPJob<RealType>>& joblist,
+                                    std::vector<RealType>& pairpots,
+                                    ResourceCollection& collection)
+{
+  auto& soecp_component_leader = soecp_component_list.getLeader();
+  if (soecp_component_leader.VP_)
+  {
+    // Compute ratios with VP
+    RefVectorWithLeader<VirtualParticleSet> vp_list(*soecp_component_leader.VP_);
+    RefVectorWithLeader<const VirtualParticleSet> const_vp_list(*soecp_component_leader.VP_);
+    RefVector<const std::vector<PosType>> deltaV_list;
+    RefVector<const std::vector<RealType>> deltaS_list;
+    RefVector<std::vector<ValueType>> psiratios_list;
+    vp_list.reserve(soecp_component_list.size());
+    const_vp_list.reserve(soecp_component_list.size());
+    deltaV_list.reserve(soecp_component_list.size());
+    deltaS_list.reserve(soecp_component_list.size());
+    psiratios_list.reserve(soecp_component_list.size());
+
+    for (size_t i = 0; i < soecp_component_list.size(); i++)
+    {
+      SOECPComponent& component(soecp_component_list[i]);
+      const NLPPJob<RealType>& job = joblist[i];
+      const RealType sold          = p_list[i].spins[job.electron_id];
+
+      component.buildTotalQuadrature(job.ion_elec_dist, job.ion_elec_displ, sold);
+
+      vp_list.push_back(*component.VP_);
+      const_vp_list.push_back(*component.VP_);
+      deltaV_list.push_back(component.deltaV_);
+      deltaS_list.push_back(component.deltaS_);
+      psiratios_list.push_back(component.psiratio_);
+    }
+
+    ResourceCollectionTeamLock<VirtualParticleSet> vp_res_lock(collection, vp_list);
+
+    VirtualParticleSet::mw_makeMovesWithSpin(vp_list, p_list, deltaV_list, deltaS_list, joblist, true);
+
+    TrialWaveFunction::mw_evaluateRatios(psi_list, const_vp_list, psiratios_list);
+  }
+  else
+  {
+    // Compute ratios without VP. Slow
+    for (size_t i = 0; i < p_list.size(); i++)
+    {
+      SOECPComponent& component(soecp_component_list[i]);
+      ParticleSet& W(p_list[i]);
+      TrialWaveFunction& psi(psi_list[i]);
+      const NLPPJob<RealType>& job = joblist[i];
+
+      const RealType sold = W.spins[job.electron_id];
+      component.buildTotalQuadrature(job.ion_elec_dist, job.ion_elec_displ, sold);
+
+      for (int j = 0; j < component.total_knots_; j++)
+      {
+        W.makeMoveWithSpin(job.electron_id, component.deltaV_[j], component.deltaS_[j]);
+        component.psiratio_[j] = psi.calcRatio(W, job.electron_id);
+        W.rejectMove(job.electron_id);
+        psi.resetPhaseDiff();
+      }
+    }
+  }
+
+  for (size_t i = 0; i < p_list.size(); i++)
+  {
+    SOECPComponent& component(soecp_component_list[i]);
+    const NLPPJob<RealType>& job = joblist[i];
+    const RealType sold          = p_list[i].spins[job.electron_id];
+    pairpots[i]                  = component.calculateProjector(job.ion_elec_dist, job.ion_elec_displ, sold);
+  }
+}
+
 SOECPComponent::RealType SOECPComponent::evaluateValueAndDerivatives(ParticleSet& W,
                                                                      int iat,
                                                                      TrialWaveFunction& Psi,
@@ -317,6 +394,10 @@ void SOECPComponent::buildTotalQuadrature(const RealType r, const PosType& dr, c
   addSpatialQuadrature(sknot_, r, dr, smax - sold, RealType(1. / 3.) * dS);
 
   assert(count == total_knots_);
+
+  //also set the radial function
+  for (int ip = 0; ip < nchannel_; ip++)
+    vrad_[ip] = sopp_m_[ip]->splint(r);
 }
 
 void SOECPComponent::rotateQuadratureGrid(const TensorType& rmat)

src/QMCHamiltonians/SOECPComponent.h

@@ -16,13 +16,19 @@
 #include "QMCHamiltonians/OperatorBase.h"
 #include "QMCHamiltonians/RandomRotationMatrix.h"
 #include "QMCWaveFunctions/TrialWaveFunction.h"
+#include <ResourceCollection.h>
 #include "Numerics/OneDimGridBase.h"
 #include "Numerics/OneDimGridFunctor.h"
 #include "Numerics/OneDimLinearSpline.h"
 #include "Numerics/OneDimCubicSpline.h"
 
 namespace qmcplusplus
 {
+
+namespace testing
+{
+class TestSOECPotential;
+}
 /** class SOECPComponent
  **  brief Computes the nonlocal spin-orbit interaction \f$\Delta V_SO(r) |ljm_j><ljm_j|\f$.
  **  details This computes the nonlocal spin-orbit interaction between a single ion species and 
@@ -72,7 +78,7 @@ class SOECPComponent : public QMCTraits
   Vector<ValueType> dlogpsi_vp_;
   VirtualParticleSet* VP_;
 
-  //This builds the full quadrature grid for the Simpsons rule used for spin integrals as well as 
+  //This builds the full quadrature grid for the Simpsons rule used for spin integrals as well as
   //the spatial quadrature. In this function, it specifies the deltaS_ and deltaV_ for all the quadrature points and sets the interal weights
   //in spin_quad_weights
   //If there are s0,s1,...sN spin integral points and q0,q1,...qM spatial quadrature points, the order is
@@ -108,6 +114,16 @@ class SOECPComponent : public QMCTraits
    */
   RealType evaluateOne(ParticleSet& W, int iat, TrialWaveFunction& Psi, int iel, RealType r, const PosType& dr);
 
+  RealType calculateProjector(RealType r, const PosType& dr, RealType sold);
+
+
+  static void mw_evaluateOne(const RefVectorWithLeader<SOECPComponent>& soecp_component_list,
+                             const RefVectorWithLeader<ParticleSet>& p_list,
+                             const RefVectorWithLeader<TrialWaveFunction>& psi_list,
+                             const RefVector<const NLPPJob<RealType>>& joblist,
+                             std::vector<RealType>& pairpots,
+                             ResourceCollection& collection);
+
   RealType evaluateValueAndDerivatives(ParticleSet& P,
                                        int iat,
                                        TrialWaveFunction& psi,
@@ -130,8 +146,12 @@ class SOECPComponent : public QMCTraits
   inline int getNknot() const { return nknot_; }
   inline int getSknot() const { return sknot_; }
 
+  const VirtualParticleSet* getVP() const { return VP_; };
+
   friend struct ECPComponentBuilder;
   friend void copyGridUnrotatedForTest(SOECPComponent& nlpp);
+
+  friend class testing::TestSOECPotential;
 };
 
 } // namespace qmcplusplus