Merge pull request #4933 from camelto2/fast_soc

Fast, efficient implementation of SOC

docs/hamiltonianobservable.rst

@@ -305,31 +305,33 @@ the radial functions :math:`V_{\ell}^{\rm SO}` can be included in the pseudopote
 
 attributes:
 
-  +-----------------------------+--------------+-----------------------+------------------------+--------------------------------------------------+
-  | **Name**                    | **Datatype** | **Values**            | **Default**            | **Description**                                  |
-  +=============================+==============+=======================+========================+==================================================+
-  | ``type``:math:`^r`          | text         | **pseudo**            |                        | Must be pseudo                                   |
-  +-----------------------------+--------------+-----------------------+------------------------+--------------------------------------------------+
-  | ``name/id``:math:`^r`       | text         | *anything*            | PseudoPot              | *No current function*                            |
-  +-----------------------------+--------------+-----------------------+------------------------+--------------------------------------------------+
-  | ``source``:math:`^r`        | text         | ``particleset.name``  | i                      | Ion ``particleset`` name                         |
-  +-----------------------------+--------------+-----------------------+------------------------+--------------------------------------------------+
-  | ``target``:math:`^r`        | text         | ``particleset.name``  | ``hamiltonian.target`` | Electron ``particleset`` name                    |
-  +-----------------------------+--------------+-----------------------+------------------------+--------------------------------------------------+
-  | ``pbc``:math:`^o`           | boolean      | yes/no                | yes*                   | Use Ewald summation                              |
-  +-----------------------------+--------------+-----------------------+------------------------+--------------------------------------------------+
-  | ``forces``                  | boolean      | yes/no                | no                     | *Deprecated*                                     |
-  +-----------------------------+--------------+-----------------------+------------------------+--------------------------------------------------+
-  | ``wavefunction``:math:`^r`  | text         | ``wavefunction.name`` | invalid                | Identify wavefunction                            |
-  +-----------------------------+--------------+-----------------------+------------------------+--------------------------------------------------+
-  | ``format``:math:`^r`        | text         | xml/table             | table                  | Select file format                               |
-  +-----------------------------+--------------+-----------------------+------------------------+--------------------------------------------------+
-  | ``algorithm``:math:`^o`     | text         | batched/non-batched   | batched                | Choose NLPP algorithm                            |
-  +-----------------------------+--------------+-----------------------+------------------------+--------------------------------------------------+
-  | ``DLA``:math:`^o`           | text         | yes/no                | no                     | Use determinant localization approximation       |
-  +-----------------------------+--------------+-----------------------+------------------------+--------------------------------------------------+
-  | ``physicalSO``:math:`^o`    | boolean      | yes/no                | yes                    | Include the SO contribution in the local energy  |
-  +-----------------------------+--------------+-----------------------+------------------------+--------------------------------------------------+
+  +------------------------------+--------------+-----------------------+------------------------+--------------------------------------------------+
+  | **Name**                     | **Datatype** | **Values**            | **Default**            | **Description**                                  |
+  +==============================+==============+=======================+========================+==================================================+
+  | ``type``:math:`^r`           | text         | **pseudo**            |                        | Must be pseudo                                   |
+  +------------------------------+--------------+-----------------------+------------------------+--------------------------------------------------+
+  | ``name/id``:math:`^r`        | text         | *anything*            | PseudoPot              | *No current function*                            |
+  +------------------------------+--------------+-----------------------+------------------------+--------------------------------------------------+
+  | ``source``:math:`^r`         | text         | ``particleset.name``  | i                      | Ion ``particleset`` name                         |
+  +------------------------------+--------------+-----------------------+------------------------+--------------------------------------------------+
+  | ``target``:math:`^r`         | text         | ``particleset.name``  | ``hamiltonian.target`` | Electron ``particleset`` name                    |
+  +------------------------------+--------------+-----------------------+------------------------+--------------------------------------------------+
+  | ``pbc``:math:`^o`            | boolean      | yes/no                | yes*                   | Use Ewald summation                              |
+  +------------------------------+--------------+-----------------------+------------------------+--------------------------------------------------+
+  | ``forces``                   | boolean      | yes/no                | no                     | *Deprecated*                                     |
+  +------------------------------+--------------+-----------------------+------------------------+--------------------------------------------------+
+  | ``wavefunction``:math:`^r`   | text         | ``wavefunction.name`` | invalid                | Identify wavefunction                            |
+  +------------------------------+--------------+-----------------------+------------------------+--------------------------------------------------+
+  | ``format``:math:`^r`         | text         | xml/table             | table                  | Select file format                               |
+  +------------------------------+--------------+-----------------------+------------------------+--------------------------------------------------+
+  | ``algorithm``:math:`^o`      | text         | batched/non-batched   | batched                | Choose NLPP algorithm                            |
+  +------------------------------+--------------+-----------------------+------------------------+--------------------------------------------------+
+  | ``DLA``:math:`^o`            | text         | yes/no                | no                     | Use determinant localization approximation       |
+  +------------------------------+--------------+-----------------------+------------------------+--------------------------------------------------+
+  | ``physicalSO``:math:`^o`     | boolean      | yes/no                | yes                    | Include the SO contribution in the local energy  |
+  +------------------------------+--------------+-----------------------+------------------------+--------------------------------------------------+
+  | ``spin_integrator``:math:`^o`| text         | exact / simpson       | exact                  | Choose which spin integration technique to use   |
+  +------------------------------+--------------+-----------------------+------------------------+--------------------------------------------------+
 
 Additional information:
 
@@ -374,6 +376,12 @@ Additional information:
    ``.xml`` file, this flag allows control over whether the SO contribution
    is included in the local energy. 
 
+-  **spin_integrator** Selects which spin integration technique to use.
+   ``simpson`` uses a numerical integration scheme
+   which can be inefficient but was previously the default. The ``exact`` method exploits 
+   the structure of the Slater-Jastrow wave function in order to analytically 
+   perform the spin integral.
+
 .. code-block::
   :caption: QMCPXML element for pseudopotential electron-ion interaction (psf files).
   :name: Listing 19

src/QMCHamiltonians/ECPComponentBuilder.cpp

@@ -28,13 +28,13 @@
 
 namespace qmcplusplus
 {
-ECPComponentBuilder::ECPComponentBuilder(const std::string& aname, Communicate* c, int nrule, int llocal)
+ECPComponentBuilder::ECPComponentBuilder(const std::string& aname, Communicate* c, int nrule, int llocal, int srule)
     : MPIObjectBase(c),
       NumNonLocal(0),
       Lmax(0),
       Llocal(llocal),
       Nrule(nrule),
-      Srule(8),
+      Srule(srule),
       AtomicNumber(0),
       Zeff(0),
       RcutMax(-1),

src/QMCHamiltonians/ECPComponentBuilder.h

@@ -48,7 +48,10 @@ struct ECPComponentBuilder : public MPIObjectBase, public QMCTraits
   std::unique_ptr<L2RadialPotential> pp_L2;
   std::map<std::string, int> angMon;
 
-  ECPComponentBuilder(const std::string& aname, Communicate* c, int nrule = -1, int llocal = -1);
+  /** constructor
+   * spin grid used for numerical integration. use 0 for exact integration.
+   */
+  ECPComponentBuilder(const std::string& aname, Communicate* c, int nrule = -1, int llocal = -1, int srule = 8);
 
   bool parse(const std::string& fname, xmlNodePtr cur);
   bool put(xmlNodePtr cur);

src/QMCHamiltonians/ECPotentialBuilder.cpp

@@ -65,6 +65,7 @@ bool ECPotentialBuilder::put(xmlNodePtr cur)
   std::string pbc;
   std::string forces;
   std::string physicalSO;
+  std::string spin_integrator;
 
   OhmmsAttributeSet pAttrib;
   pAttrib.add(ecpFormat, "format", {"table", "xml"});
@@ -73,11 +74,14 @@ bool ECPotentialBuilder::put(xmlNodePtr cur)
   pAttrib.add(pbc, "pbc", {"yes", "no"});
   pAttrib.add(forces, "forces", {"no", "yes"});
   pAttrib.add(physicalSO, "physicalSO", {"yes", "no"});
+  pAttrib.add(spin_integrator, "spin_integrator", {"exact", "simpson"});
   pAttrib.put(cur);
 
   bool doForces = (forces == "yes") || (forces == "true");
   if (use_DLA == "yes")
     app_log() << "    Using determinant localization approximation (DLA)" << std::endl;
+
+  use_exact_spin = (spin_integrator == "exact") ? true : false;
   if (ecpFormat == "xml")
   {
     useXmlFormat(cur);
@@ -150,7 +154,7 @@ bool ECPotentialBuilder::put(xmlNodePtr cur)
     else
       APP_ABORT("physicalSO must be set to yes/no. Unknown option given\n");
 
-    std::unique_ptr<SOECPotential> apot = std::make_unique<SOECPotential>(IonConfig, targetPtcl, targetPsi);
+    std::unique_ptr<SOECPotential> apot = std::make_unique<SOECPotential>(IonConfig, targetPtcl, targetPsi, use_exact_spin);
     int nknot_max                       = 0;
     int sknot_max                       = 0;
     for (int i = 0; i < soPot.size(); i++)
@@ -166,10 +170,12 @@ bool ECPotentialBuilder::put(xmlNodePtr cur)
     }
     app_log() << "\n  Using SOECP potential \n"
               << "    Maximum grid on a sphere for SOECPotential: " << nknot_max << std::endl;
-    app_log() << "    Maximum grid for Simpson's rule for spin integral: " << sknot_max << std::endl;
+    if (use_exact_spin)
+      app_log() << "    Using fast SOECP evaluation. Spin integration is exact" << std::endl;
+    else
+      app_log() << "    Maximum grid for Simpson's rule for spin integral: " << sknot_max << std::endl;
     if (NLPP_algo == "batched")
       app_log() << "    Using batched ratio computing in SOECP potential" << std::endl;
-
     if (physicalSO == "yes")
       targetH.addOperator(std::move(apot), "SOECP"); //default is physical operator
     else
@@ -225,7 +231,9 @@ void ECPotentialBuilder::useXmlFormat(xmlNodePtr cur)
       {
         app_log() << std::endl << "  Adding pseudopotential for " << ionName << std::endl;
 
-        ECPComponentBuilder ecp(ionName, myComm, nrule, llocal);
+        //Use simpsons rule for spin integral if not using exact spin integration
+        int srule = use_exact_spin ? 0 : 8;
+        ECPComponentBuilder ecp(ionName, myComm, nrule, llocal, srule);
         if (format == "xml")
         {
           if (href == "none")

src/QMCHamiltonians/ECPotentialBuilder.h

@@ -33,6 +33,7 @@ struct ECPotentialBuilder : public MPIObjectBase, public QMCTraits
   bool hasNonLocalPot;
   bool hasSOPot;
   bool hasL2Pot;
+  bool use_exact_spin;
 
   QMCHamiltonian& targetH;
   ParticleSet& IonConfig;

src/QMCHamiltonians/SOECPComponent.cpp

@@ -71,10 +71,8 @@ void SOECPComponent::resize_warrays(int n, int m, int s)
   nchannel_ = sopp_m_.size();
   nknot_    = sgridxyz_m_.size();
   sknot_    = s;
-  if (sknot_ < 2)
-    throw std::runtime_error("Spin knots must be >= 2\n");
-  if (sknot_ % 2 != 0)
-    throw std::runtime_error("Spin knots uses Simpson's rule. Must have even number of knots");
+  if (sknot_ % 2 != 0 && sknot_ > 0)
+    throw std::runtime_error("Spin integration uses Simpson's rule. Must have even number of knots in this case");
 
   //Need +1 for Simpsons rule to include both end points.
   //sknot here refers to the number of subintervals for integration
@@ -135,6 +133,25 @@ SOECPComponent::ComplexType SOECPComponent::lmMatrixElements(int l, int m1, int
   }
 }
 
+SOECPComponent::ComplexType SOECPComponent::matrixElementDecomposed(int l, int m1, int m2, RealType spin, bool plus)
+{
+  RealType pm = plus ? RealType(1) : RealType(-1);
+  RealType mp = -1 * pm;
+
+
+  ComplexType lx = lmMatrixElements(l, m1, m2, 0);
+  ComplexType ly = lmMatrixElements(l, m1, m2, 1);
+  ComplexType lz = lmMatrixElements(l, m1, m2, 2);
+
+  RealType cos2s = std::cos(2 * spin);
+  RealType sin2s = std::sin(2 * spin);
+  ComplexType phase(cos2s, mp * sin2s);
+  ComplexType eye(0, 1);
+  RealType onehalf = 0.5;
+  ComplexType val  = onehalf * (pm * lz + phase * (lx + pm * eye * ly));
+  return val;
+}
+
 SOECPComponent::RealType SOECPComponent::evaluateOne(ParticleSet& W,
                                                      int iat,
                                                      TrialWaveFunction& Psi,
@@ -192,6 +209,67 @@ SOECPComponent::RealType SOECPComponent::calculateProjector(RealType r, const Po
   return std::real(pairpot);
 }
 
+SOECPComponent::RealType SOECPComponent::evaluateOneExactSpinIntegration(ParticleSet& W,
+                                                                         const int iat,
+                                                                         const TrialWaveFunction& psi,
+                                                                         const int iel,
+                                                                         const RealType r,
+                                                                         const PosType& dr)
+{
+#ifdef QMC_COMPLEX
+  RealType sold = W.spins[iel];
+
+  using ValueVector = SPOSet::ValueVector;
+  std::pair<ValueVector, ValueVector> spinor_multiplier;
+  auto& up_row = spinor_multiplier.first;
+  auto& dn_row = spinor_multiplier.second;
+  up_row.resize(nknot_);
+  dn_row.resize(nknot_);
+
+  //buildQuadraturePointDeltaPositions
+  for (int j = 0; j < nknot_; j++)
+    deltaV_[j] = r * rrotsgrid_m_[j] - dr;
+  vp_->makeMoves(W, iel, deltaV_, true, iat);
+
+  //calculate radial potentials
+  for (int ip = 0; ip < nchannel_; ip++)
+    vrad_[ip] = sopp_m_[ip]->splint(r);
+
+  for (int iq = 0; iq < nknot_; iq++)
+  {
+    up_row[iq] = 0.0;
+    dn_row[iq] = 0.0;
+    for (int il = 0; il < nchannel_; il++)
+    {
+      int l = il + 1;
+      for (int m1 = -l; m1 <= l; m1++)
+      {
+        ComplexType Y = sphericalHarmonic(l, m1, dr);
+        for (int m2 = -l; m2 <= l; m2++)
+        {
+          ComplexType cY    = std::conj(sphericalHarmonic(l, m2, rrotsgrid_m_[iq]));
+          ComplexType so_up = matrixElementDecomposed(l, m1, m2, sold, true);
+          ComplexType so_dn = matrixElementDecomposed(l, m1, m2, sold, false);
+          RealType fourpi   = 4.0 * M_PI;
+          //Note: Need 4pi weight. In Normal NonLocalECP, 1/4Pi generated from transformation to legendre polynomials and gets absorbed into the
+          // quadrature integration. We don't get the 1/4Pi from legendre here, so we need to scale by 4Pi.
+          up_row[iq] += fourpi * vrad_[il] * Y * cY * so_up;
+          dn_row[iq] += fourpi * vrad_[il] * Y * cY * so_dn;
+        }
+      }
+    }
+  }
+
+  psi.evaluateSpinorRatios(*vp_, spinor_multiplier, psiratio_);
+  ComplexType pairpot;
+  for (size_t iq = 0; iq < nknot_; iq++)
+    pairpot += psiratio_[iq] * sgridweight_m_[iq];
+  return std::real(pairpot);
+#else
+  throw std::runtime_error("SOECPComponent::evaluateOneBodyOpMatrixContribution only implemented in complex build");
+#endif
+}
+
 void SOECPComponent::mw_evaluateOne(const RefVectorWithLeader<SOECPComponent>& soecp_component_list,
                                     const RefVectorWithLeader<ParticleSet>& p_list,
                                     const RefVectorWithLeader<TrialWaveFunction>& psi_list,

src/QMCHamiltonians/SOECPComponent.h

@@ -43,6 +43,7 @@ class SOECPComponent : public QMCTraits
   using SpherGridType       = std::vector<PosType>;
   using GridType            = OneDimGridBase<RealType>;
   using RadialPotentialType = OneDimCubicSpline<RealType>;
+  using ValueMatrix         = SPOSet::ValueMatrix;
 
   ///Non Local part: angular momentum, potential and grid
   int lmax_;
@@ -58,9 +59,10 @@ class SOECPComponent : public QMCTraits
   ///Non-Local part of the pseudo-potential
   std::vector<RadialPotentialType*> sopp_m_;
 
-  ComplexType sMatrixElements(RealType s1, RealType s2, int dim);
-  ComplexType lmMatrixElements(int l, int m1, int m2, int dim);
-  int kroneckerDelta(int x, int y);
+  static ComplexType sMatrixElements(RealType s1, RealType s2, int dim);
+  static ComplexType lmMatrixElements(int l, int m1, int m2, int dim);
+  static ComplexType matrixElementDecomposed(int l, int m1, int m2, RealType spin, bool plus = true);
+  static int kroneckerDelta(int x, int y);
 
   std::vector<PosType> deltaV_;
   std::vector<RealType> deltaS_;
@@ -116,6 +118,7 @@ class SOECPComponent : public QMCTraits
 
   RealType calculateProjector(RealType r, const PosType& dr, RealType sold);
 
+  RealType evaluateOneExactSpinIntegration(ParticleSet& W, const int iat, const TrialWaveFunction& psi, const int iel, const RealType r, const PosType& dr);
 
   static void mw_evaluateOne(const RefVectorWithLeader<SOECPComponent>& soecp_component_list,
                              const RefVectorWithLeader<ParticleSet>& p_list,