Merge pull request #5173 from camelto2/fast_soc_evaluateValueAndDeriv…

…ative

Allow exact spin integration for SOECP alongside wave function optimization

src/QMCHamiltonians/SOECPComponent.cpp

@@ -185,6 +185,8 @@ SOECPComponent::RealType SOECPComponent::evaluateOne(ParticleSet& W,
 
 SOECPComponent::RealType SOECPComponent::calculateProjector(RealType r, const PosType& dr, RealType sold)
 {
+#ifdef QMC_COMPLEX
+  wvec_.resize(total_knots_); //contribution from each quarature point
   ComplexType pairpot;
   for (int iq = 0; iq < total_knots_; iq++)
   {
@@ -208,9 +210,13 @@ SOECPComponent::RealType SOECPComponent::calculateProjector(RealType r, const Po
       }
       lsum += vrad_[il] * msums;
     }
-    pairpot += psiratio_[iq] * lsum * spin_quad_weights_[iq];
+    wvec_[iq] = psiratio_[iq] * lsum * spin_quad_weights_[iq];
+    pairpot += wvec_[iq];
   }
   return std::real(pairpot);
+#else
+  throw std::runtime_error("SOECPComponent::calculateProjector only implemented in complex build.");
+#endif
 }
 
 void SOECPComponent::setupExactSpinProjector(RealType r, const PosType& dr, RealType sold)
@@ -429,7 +435,6 @@ SOECPComponent::RealType SOECPComponent::evaluateValueAndDerivatives(ParticleSet
   const size_t num_vars = optvars.num_active_vars;
   dratio_.resize(total_knots_, num_vars);
   dlogpsi_vp_.resize(dlogpsi.size());
-  wvec_.resize(total_knots_);
 
   RealType sold = W.spins[iel];
   buildTotalQuadrature(r, dr, sold);
@@ -462,38 +467,48 @@ SOECPComponent::RealType SOECPComponent::evaluateValueAndDerivatives(ParticleSet
       W.acceptMove(iel);
     }
 
-  ComplexType pairpot;
-  for (int iq = 0; iq < total_knots_; iq++)
-  {
-    ComplexType lsum;
-    for (int il = 0; il < nchannel_; il++)
-    {
-      int l = il + 1;
-      ComplexType msums;
-      for (int m1 = -l; m1 <= l; m1++)
-      {
-        ComplexType Y = sphericalHarmonic(l, m1, dr);
-        for (int m2 = -l; m2 <= l; m2++)
-        {
-          ComplexType ldots;
-          for (int id = 0; id < 3; id++)
-            ldots += lmMatrixElements(l, m1, m2, id) * sMatrixElements(W.spins[iel], W.spins[iel] + deltaS_[iq], id);
-          ComplexType cY = std::conj(sphericalHarmonic(l, m2, rrotsgrid_m_[iq % nknot_]));
-          msums += Y * cY * ldots;
-        }
-      }
-      lsum += sopp_m_[il]->splint(r) * msums;
-    }
-    wvec_[iq] = lsum * psiratio_[iq] * spin_quad_weights_[iq];
-    pairpot += wvec_[iq];
-  }
-
+  ComplexType pairpot = calculateProjector(r, dr, sold);
+  //calculateProjector stores quad points in wvec_
   BLAS::gemv('N', num_vars, total_knots_, 1.0, dratio_.data(), num_vars, wvec_.data(), 1, 1.0, dhpsioverpsi.data(), 1);
 
   return std::real(pairpot);
 #endif
 }
 
+SOECPComponent::RealType SOECPComponent::evaluateValueAndDerivativesExactSpinIntegration(
+    ParticleSet& W,
+    int iat,
+    TrialWaveFunction& Psi,
+    int iel,
+    RealType r,
+    const PosType& dr,
+    const opt_variables_type& optvars,
+    const Vector<ValueType>& dlogpsi,
+    Vector<ValueType>& dhpsioverpsi)
+{
+#ifndef QMC_COMPLEX
+  throw std::runtime_error("SOECPComponent::evaluateValueAndDerivatives should not be called in real build\n");
+#else
+  const size_t num_vars = optvars.num_active_vars;
+  dratio_.resize(total_knots_, num_vars);
+  dlogpsi_vp_.resize(dlogpsi.size());
+
+  RealType sold = W.spins[iel];
+  for (int j = 0; j < nknot_; j++)
+    deltaV_[j] = r * rrotsgrid_m_[j] - dr;
+  vp_->makeMoves(W, iel, deltaV_, true, iat);
+  setupExactSpinProjector(r, dr, sold);
+  Psi.evaluateSpinorDerivRatios(*vp_, spinor_multiplier_, optvars, psiratio_, dratio_);
+
+  BLAS::gemv('N', num_vars, total_knots_, 1.0, dratio_.data(), num_vars, psiratio_.data(), 1, 1.0, dhpsioverpsi.data(), 1);
+
+  ComplexType pairpot;
+  for (size_t iq = 0; iq < nknot_; iq++)
+    pairpot += psiratio_[iq];
+  return std::real(pairpot);
+#endif
+}
+
 void SOECPComponent::buildTotalQuadrature(const RealType r, const PosType& dr, const RealType sold)
 {
   int count = 0;

src/QMCHamiltonians/SOECPComponent.h

@@ -153,6 +153,16 @@ class SOECPComponent : public QMCTraits
                                        const Vector<ValueType>& dlogpsi,
                                        Vector<ValueType>& dhpsioverpsi);
 
+  RealType evaluateValueAndDerivativesExactSpinIntegration(ParticleSet& P,
+                                                           int iat,
+                                                           TrialWaveFunction& psi,
+                                                           int iel,
+                                                           RealType r,
+                                                           const PosType& dr,
+                                                           const opt_variables_type& optvars,
+                                                           const Vector<ValueType>& dlogpsi,
+                                                           Vector<ValueType>& dhpsioverpsi);
+
   void print(std::ostream& os);
 
   void initVirtualParticle(const ParticleSet& qp);

src/QMCHamiltonians/tests/test_ecp.cpp

@@ -546,26 +546,30 @@ TEST_CASE("Evaluate_soecp", "[hamiltonian]")
   //Need to set up temporary data for this configuration in trial wavefunction.  Needed for ratios.
   auto logpsi = psi.evaluateLog(elec);
 
-  auto test_evaluateOne = [&]() {
+  auto test_evaluateOne = [&](bool exact) {
     RealType Value1(0.0);
     for (int jel = 0; jel < elec.getTotalNum(); jel++)
     {
       const auto& dist  = myTable.getDistRow(jel);
       const auto& displ = myTable.getDisplRow(jel);
       for (int iat = 0; iat < ions.getTotalNum(); iat++)
         if (sopp != nullptr && dist[iat] < sopp->getRmax())
-          Value1 += sopp->evaluateOne(elec, iat, psi, jel, dist[iat], RealType(-1) * displ[iat]);
+          if (exact)
+            Value1 += sopp->evaluateOneExactSpinIntegration(elec, iat, psi, jel, dist[iat], RealType(-1) * displ[iat]);
+          else
+            Value1 += sopp->evaluateOne(elec, iat, psi, jel, dist[iat], RealType(-1) * displ[iat]);
     }
     REQUIRE(Value1 == Approx(-3.530511241));
   };
 
   {
     //test with VPs
     sopp->initVirtualParticle(elec);
-    test_evaluateOne();
+    test_evaluateOne(false);
+    test_evaluateOne(true);
     sopp->deleteVirtualParticle();
     //test without VPs
-    test_evaluateOne();
+    test_evaluateOne(false);
   }
 
   //Check evaluateValueAndDerivatives
@@ -591,7 +595,7 @@ TEST_CASE("Evaluate_soecp", "[hamiltonian]")
                                              -8.169955304e-14};
 
 
-  auto test_evaluateValueAndDerivatives = [&]() {
+  auto test_evaluateValueAndDerivatives = [&](bool exact) {
     dlogpsi.resize(NumOptimizables, ValueType(0));
     dhpsioverpsi.resize(NumOptimizables, ValueType(0));
     psi.evaluateDerivatives(elec, optvars, dlogpsi, dhpsioverpsi);
@@ -609,8 +613,12 @@ TEST_CASE("Evaluate_soecp", "[hamiltonian]")
       const auto& displ = myTable.getDisplRow(jel);
       for (int iat = 0; iat < ions.getTotalNum(); iat++)
         if (sopp != nullptr && dist[iat] < sopp->getRmax())
-          Value1 += sopp->evaluateValueAndDerivatives(elec, iat, psi, jel, dist[iat], -displ[iat], optvars, dlogpsi,
-                                                      dhpsioverpsi);
+          if (exact)
+            Value1 += sopp->evaluateValueAndDerivativesExactSpinIntegration(elec, iat, psi, jel, dist[iat], -displ[iat],
+                                                                            optvars, dlogpsi, dhpsioverpsi);
+          else
+            Value1 += sopp->evaluateValueAndDerivatives(elec, iat, psi, jel, dist[iat], -displ[iat], optvars, dlogpsi,
+                                                        dhpsioverpsi);
     }
     REQUIRE(Value1 == Approx(-3.530511241).epsilon(2.e-5));
 
@@ -620,9 +628,10 @@ TEST_CASE("Evaluate_soecp", "[hamiltonian]")
 
   {
     sopp->initVirtualParticle(elec);
-    test_evaluateValueAndDerivatives();
+    test_evaluateValueAndDerivatives(false);
+    test_evaluateValueAndDerivatives(true);
     sopp->deleteVirtualParticle();
-    test_evaluateValueAndDerivatives();
+    test_evaluateValueAndDerivatives(false);
   }
 }
 #endif

src/QMCWaveFunctions/Fermion/DiracDeterminant.cpp

@@ -448,7 +448,9 @@ void DiracDeterminant<DU_TYPE>::evaluateRatios(const VirtualParticleSet& VP, std
 }
 
 template<typename DU_TYPE>
-void DiracDeterminant<DU_TYPE>::evaluateSpinorRatios(const VirtualParticleSet& VP, const std::pair<ValueVector, ValueVector>& spinor_multiplier, std::vector<ValueType>& ratios)
+void DiracDeterminant<DU_TYPE>::evaluateSpinorRatios(const VirtualParticleSet& VP,
+                                                     const std::pair<ValueVector, ValueVector>& spinor_multiplier,
+                                                     std::vector<ValueType>& ratios)
 {
   {
     ScopedTimer local_timer(RatioTimer);
@@ -521,6 +523,19 @@ void DiracDeterminant<DU_TYPE>::evaluateDerivRatios(const VirtualParticleSet& VP
   Phi->evaluateDerivRatios(VP, optvars, psiV, invRow, ratios, dratios, FirstIndex, LastIndex);
 }
 
+template<typename DU_TYPE>
+void DiracDeterminant<DU_TYPE>::evaluateSpinorDerivRatios(const VirtualParticleSet& VP,
+                                                          const std::pair<ValueVector, ValueVector>& spinor_multiplier,
+                                                          const opt_variables_type& optvars,
+                                                          std::vector<ValueType>& ratios,
+                                                          Matrix<ValueType>& dratios)
+{
+  const int WorkingIndex = VP.refPtcl - FirstIndex;
+  assert(WorkingIndex >= 0);
+  std::copy_n(psiM[WorkingIndex], invRow.size(), invRow.data());
+  Phi->evaluateSpinorDerivRatios(VP, spinor_multiplier, optvars, psiV, invRow, ratios, dratios, FirstIndex, LastIndex);
+}
+
 template<typename DU_TYPE>
 void DiracDeterminant<DU_TYPE>::evaluateRatiosAlltoOne(ParticleSet& P, std::vector<ValueType>& ratios)
 {

src/QMCWaveFunctions/Fermion/DiracDeterminant.h

@@ -106,7 +106,9 @@ class DiracDeterminant : public DiracDeterminantBase
    */
   void evaluateRatios(const VirtualParticleSet& VP, std::vector<ValueType>& ratios) override;
 
-  void evaluateSpinorRatios(const VirtualParticleSet& VP, const std::pair<ValueVector, ValueVector>& spinor_multipler, std::vector<ValueType>& ratios) override;
+  void evaluateSpinorRatios(const VirtualParticleSet& VP,
+                            const std::pair<ValueVector, ValueVector>& spinor_multipler,
+                            std::vector<ValueType>& ratios) override;
 
   void mw_evaluateRatios(const RefVectorWithLeader<WaveFunctionComponent>& wfc_list,
                          const RefVectorWithLeader<const VirtualParticleSet>& vp_list,
@@ -117,6 +119,12 @@ class DiracDeterminant : public DiracDeterminantBase
                            std::vector<ValueType>& ratios,
                            Matrix<ValueType>& dratios) override;
 
+  void evaluateSpinorDerivRatios(const VirtualParticleSet& VP,
+                                 const std::pair<ValueVector, ValueVector>& spinor_multipler,
+                                 const opt_variables_type& optvars,
+                                 std::vector<ValueType>& ratios,
+                                 Matrix<ValueType>& dratios) override;
+
   PsiValue ratioGrad(ParticleSet& P, int iat, GradType& grad_iat) override;
 
   PsiValue ratioGradWithSpin(ParticleSet& P, int iat, GradType& grad_iat, ComplexType& spingrad) final;

src/QMCWaveFunctions/Fermion/DiracDeterminantBatched.cpp

@@ -671,7 +671,7 @@ void DiracDeterminantBatched<PL, VT, FPVT>::mw_evaluateGL(const RefVectorWithLea
 
     for (int iw = 0; iw < nw; iw++)
     {
-      auto& det = wfc_list.getCastedElement<DiracDeterminantBatched<PL, VT, FPVT>>(iw);
+      auto& det      = wfc_list.getCastedElement<DiracDeterminantBatched<PL, VT, FPVT>>(iw);
       det.UpdateMode = ORB_WALKER;
       det.computeGL(G_list[iw], L_list[iw]);
     }
@@ -806,7 +806,7 @@ void DiracDeterminantBatched<PL, VT, FPVT>::evaluateRatios(const VirtualParticle
 template<PlatformKind PL, typename VT, typename FPVT>
 void DiracDeterminantBatched<PL, VT, FPVT>::evaluateSpinorRatios(
     const VirtualParticleSet& VP,
-    const std::pair<ValueVector, ValueVector>& spinor_multipler,
+    const std::pair<ValueVector, ValueVector>& spinor_multiplier,
     std::vector<Value>& ratios)
 {
   {
@@ -816,7 +816,7 @@ void DiracDeterminantBatched<PL, VT, FPVT>::evaluateSpinorRatios(
   }
   {
     ScopedTimer local_timer(SPOVTimer);
-    Phi->evaluateDetSpinorRatios(VP, psiV_host_view, spinor_multipler, d2psiV_host_view, ratios);
+    Phi->evaluateDetSpinorRatios(VP, psiV_host_view, spinor_multiplier, d2psiV_host_view, ratios);
   }
 }
 
@@ -871,6 +871,20 @@ void DiracDeterminantBatched<PL, VT, FPVT>::evaluateDerivRatios(const VirtualPar
   Phi->evaluateDerivRatios(VP, optvars, psiV_host_view, d2psiV_host_view, ratios, dratios, FirstIndex, LastIndex);
 }
 
+template<PlatformKind PL, typename VT, typename FPVT>
+void DiracDeterminantBatched<PL, VT, FPVT>::evaluateSpinorDerivRatios(
+    const VirtualParticleSet& VP,
+    const std::pair<ValueVector, ValueVector>& spinor_multiplier,
+    const opt_variables_type& optvars,
+    std::vector<ValueType>& ratios,
+    Matrix<ValueType>& dratios)
+{
+  const int WorkingIndex = VP.refPtcl - FirstIndex;
+  assert(WorkingIndex >= 0);
+  std::copy_n(psiMinv_[WorkingIndex], d2psiV.size(), d2psiV.data());
+  Phi->evaluateSpinorDerivRatios(VP, spinor_multiplier, optvars, psiV_host_view, d2psiV_host_view, ratios, dratios, FirstIndex, LastIndex);
+}
+
 template<PlatformKind PL, typename VT, typename FPVT>
 void DiracDeterminantBatched<PL, VT, FPVT>::evaluateRatiosAlltoOne(ParticleSet& P, std::vector<Value>& ratios)
 {

src/QMCWaveFunctions/Fermion/DiracDeterminantBatched.h

@@ -139,6 +139,12 @@ class DiracDeterminantBatched : public DiracDeterminantBase
                            std::vector<ValueType>& ratios,
                            Matrix<ValueType>& dratios) override;
 
+  void evaluateSpinorDerivRatios(const VirtualParticleSet& VP,
+                            const std::pair<ValueVector, ValueVector>& spinor_multiplier,
+                           const opt_variables_type& optvars,
+                           std::vector<ValueType>& ratios,
+                           Matrix<ValueType>& dratios) override;
+
   PsiValue ratioGrad(ParticleSet& P, int iat, Grad& grad_iat) override;
 
   void mw_ratioGrad(const RefVectorWithLeader<WaveFunctionComponent>& wfc_list,

src/QMCWaveFunctions/Fermion/SlaterDet.cpp

@@ -99,7 +99,7 @@ void SlaterDet::mw_evalGradWithSpin(const RefVectorWithLeader<WaveFunctionCompon
                                     const RefVectorWithLeader<ParticleSet>& p_list,
                                     int iat,
                                     std::vector<GradType>& grad_now,
-                                    std::vector<ComplexType>& spingrad_now) const 
+                                    std::vector<ComplexType>& spingrad_now) const
 {
   const int det_id = getDetID(iat);
   Dets[det_id]->mw_evalGradWithSpin(extract_DetRef_list(wfc_list, det_id), p_list, iat, grad_now, spingrad_now);
@@ -119,6 +119,15 @@ void SlaterDet::evaluateDerivRatios(const VirtualParticleSet& VP,
   return Dets[getDetID(VP.refPtcl)]->evaluateDerivRatios(VP, optvars, ratios, dratios);
 }
 
+void SlaterDet::evaluateSpinorDerivRatios(const VirtualParticleSet& VP,
+                                          const std::pair<ValueVector, ValueVector>& spinor_multiplier,
+                                          const opt_variables_type& optvars,
+                                          std::vector<ValueType>& ratios,
+                                          Matrix<ValueType>& dratios)
+{
+  return Dets[getDetID(VP.refPtcl)]->evaluateSpinorDerivRatios(VP, spinor_multiplier, optvars, ratios, dratios);
+}
+
 SlaterDet::LogValue SlaterDet::evaluateLog(const ParticleSet& P,
                                            ParticleSet::ParticleGradient& G,
                                            ParticleSet::ParticleLaplacian& L)

src/QMCWaveFunctions/Fermion/SlaterDet.h

@@ -103,7 +103,9 @@ class SlaterDet : public WaveFunctionComponent
     return Dets[getDetID(VP.refPtcl)]->evaluateRatios(VP, ratios);
   }
 
-  inline void evaluateSpinorRatios(const VirtualParticleSet& VP, const std::pair<ValueVector, ValueVector>& spinor_multiplier, std::vector<ValueType>& ratios) override
+  inline void evaluateSpinorRatios(const VirtualParticleSet& VP,
+                                   const std::pair<ValueVector, ValueVector>& spinor_multiplier,
+                                   std::vector<ValueType>& ratios) override
   {
     return Dets[getDetID(VP.refPtcl)]->evaluateSpinorRatios(VP, spinor_multiplier, ratios);
   }
@@ -113,6 +115,12 @@ class SlaterDet : public WaveFunctionComponent
                            std::vector<ValueType>& ratios,
                            Matrix<ValueType>& dratios) override;
 
+  void evaluateSpinorDerivRatios(const VirtualParticleSet& VP,
+                                 const std::pair<ValueVector, ValueVector>& spinor_multiplier,
+                                 const opt_variables_type& optvars,
+                                 std::vector<ValueType>& ratios,
+                                 Matrix<ValueType>& dratios) override;
+
   inline void mw_evaluateRatios(const RefVectorWithLeader<WaveFunctionComponent>& wfc_list,
                                 const RefVectorWithLeader<const VirtualParticleSet>& vp_list,
                                 std::vector<std::vector<ValueType>>& ratios) const override

src/QMCWaveFunctions/SPOSet.cpp

@@ -303,6 +303,25 @@ void SPOSet::evaluateDerivRatios(const VirtualParticleSet& VP,
                            "must be overloaded when the SPOSet is optimizable.");
 }
 
+void SPOSet::evaluateSpinorDerivRatios(const VirtualParticleSet& VP,
+                                       const std::pair<ValueVector, ValueVector>& spinor_multiplier,
+                                       const opt_variables_type& optvars,
+                                       ValueVector& psi,
+                                       const ValueVector& psiinv,
+                                       std::vector<ValueType>& ratios,
+                                       Matrix<ValueType>& dratios,
+                                       int FirstIndex,
+                                       int LastIndex)
+{
+  // Match the fallback in WaveFunctionComponent that evaluates just the ratios
+  evaluateDetSpinorRatios(VP, psi, spinor_multiplier, psiinv, ratios);
+
+  if (isOptimizable())
+    throw std::logic_error("Bug!! " + getClassName() +
+                           "::evaluateSpinorDerivRatios "
+                           "must be overloaded when the SPOSet is optimizable.");
+}
+
 
 /** Evaluate the derivative of the optimized orbitals with respect to the parameters
    *  this is used only for MSD, to be refined for better serving both single and multi SD