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/ECPComponentBuilder.2.cpp

@@ -135,7 +135,7 @@ void ECPComponentBuilder::buildSemiLocalAndLocal(std::vector<xmlNodePtr>& semiPt
   // we may not know which one is local yet.
 
   std::vector<int> angList;
-  std::vector<int> angListSO;      //For spin-orbit, if it exists
+  std::vector<int> angListSO; //For spin-orbit, if it exists
   std::vector<xmlNodePtr> vpsPtr;
   std::vector<xmlNodePtr> vpsoPtr; //For spin-orbit, if it exists.
   Lmax   = -1;
@@ -364,8 +364,8 @@ void ECPComponentBuilder::buildSO(const std::vector<int>& angList,
     app->spline();
     pp_so->add(angList[l], app);
   }
-  NumSO        = angList.size();
-  pp_so->Rmax_ = rmax;
+  NumSO = angList.size();
+  pp_so->setRmax(rmax);
 }
 
 bool ECPComponentBuilder::parseCasino(const std::string& fname, xmlNodePtr cur)

src/QMCHamiltonians/SOECPComponent.cpp

@@ -15,36 +15,37 @@
 #include "SOECPComponent.h"
 #include "Numerics/Ylm.h"
 #include "CPU/BLAS.hpp"
+#include "NLPPJob.h"
 
 namespace qmcplusplus
 {
 SOECPComponent::SOECPComponent()
-    : lmax_(0), nchannel_(0), nknot_(0), sknot_(0), total_knots_(0), Rmax_(-1), VP_(nullptr)
+    : lmax_(0), nchannel_(0), nknot_(0), sknot_(0), total_knots_(0), rmax_(-1), vp_(nullptr)
 {}
 
 SOECPComponent::~SOECPComponent()
 {
   for (int i = 0; i < sopp_m_.size(); i++)
     delete sopp_m_[i];
-  if (VP_)
-    delete VP_;
+  if (vp_)
+    delete vp_;
 }
 
 void SOECPComponent::print(std::ostream& os) {}
 
 void SOECPComponent::initVirtualParticle(const ParticleSet& qp)
 {
-  assert(VP_ == nullptr);
+  assert(vp_ == nullptr);
   outputManager.pause();
-  VP_ = new VirtualParticleSet(qp, total_knots_);
+  vp_ = new VirtualParticleSet(qp, total_knots_);
   outputManager.resume();
 }
 
 void SOECPComponent::deleteVirtualParticle()
 {
-  if (VP_)
-    delete VP_;
-  VP_ = nullptr;
+  if (vp_)
+    delete vp_;
+  vp_ = nullptr;
 }
 
 void SOECPComponent::add(int l, RadialPotentialType* pp)
@@ -58,8 +59,8 @@ SOECPComponent* SOECPComponent::makeClone(const ParticleSet& qp)
   SOECPComponent* myclone = new SOECPComponent(*this);
   for (int i = 0; i < sopp_m_.size(); i++)
     myclone->sopp_m_[i] = sopp_m_[i]->makeClone();
-  if (VP_)
-    myclone->VP_ = new VirtualParticleSet(qp, total_knots_);
+  if (vp_)
+    myclone->vp_ = new VirtualParticleSet(qp, total_knots_);
   return myclone;
 }
 
@@ -141,17 +142,13 @@ 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);
 
-  if (VP_)
+  if (vp_)
   {
-    VP_->makeMovesWithSpin(W, iel, deltaV_, deltaS_, true, iat);
-    Psi.evaluateRatios(*VP_, psiratio_);
+    vp_->makeMovesWithSpin(W, iel, deltaV_, deltaS_, true, iat);
+    Psi.evaluateRatios(*vp_, psiratio_);
   }
   else
     for (int iq = 0; iq < total_knots_; iq++)
@@ -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,
@@ -214,11 +291,11 @@ SOECPComponent::RealType SOECPComponent::evaluateValueAndDerivatives(ParticleSet
 
   //Now we have all the spin and spatial quadrature points acculated to use in evaluation
   //Now we need to obtain dlogpsi and dlogpsi_vp
-  if (VP_)
+  if (vp_)
   {
     // Compute ratios with VP
-    VP_->makeMovesWithSpin(W, iel, deltaV_, deltaS_, true, iat);
-    Psi.evaluateDerivRatios(*VP_, optvars, psiratio_, dratio_);
+    vp_->makeMovesWithSpin(W, iel, deltaV_, deltaS_, true, iat);
+    Psi.evaluateDerivRatios(*vp_, optvars, psiratio_, dratio_);
   }
   else
     for (int iq = 0; iq < total_knots_; iq++)
@@ -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 
@@ -46,7 +52,7 @@ class SOECPComponent : public QMCTraits
   int sknot_;
   int total_knots_; //spin + spatial knots
   ///Maximum cutoff the non-local pseudopotential
-  RealType Rmax_;
+  RealType rmax_;
   ///Angular momentum map
   aligned_vector<int> angpp_m_;
   ///Non-Local part of the pseudo-potential
@@ -70,9 +76,9 @@ class SOECPComponent : public QMCTraits
   //scratch spaces used by evaluateValueAndDerivative
   Matrix<ValueType> dratio_;
   Vector<ValueType> dlogpsi_vp_;
-  VirtualParticleSet* 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,
@@ -123,15 +139,19 @@ class SOECPComponent : public QMCTraits
   void initVirtualParticle(const ParticleSet& qp);
   void deleteVirtualParticle();
 
-  inline void setRmax(RealType rmax) { Rmax_ = rmax; }
-  inline RealType getRmax() const { return Rmax_; }
-  inline void setLmax(int Lmax) { lmax_ = Lmax; }
+  inline void setRmax(RealType rmax) { rmax_ = rmax; }
+  inline RealType getRmax() const { return rmax_; }
+  inline void setLmax(int lmax) { lmax_ = lmax; }
   inline int getLmax() const { return lmax_; }
   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