Batched driver move abstraction

docs/methods.rst

@@ -302,6 +302,9 @@ The following is an example of VMC section storing configurations (walker sample
   +--------------------------------+--------------+-------------------------+-------------+-----------------------------------------------+
   | ``debug_checks``               | text         | see additional info     | dep.        | Turn on/off additional recompute and checks   |
   +--------------------------------+--------------+-------------------------+-------------+-----------------------------------------------+
+  | ``spin_mass``                  | real         | :math:`\geq 0`          | 1.0         | Effective mass for spin sampling              |
+  +--------------------------------+--------------+-------------------------+-------------+-----------------------------------------------+
+
 
 Additional information:
 
@@ -357,6 +360,9 @@ Additional information:
 
 - ``debug_checks`` valid values are 'no', 'all', 'checkGL_after_load', 'checkGL_after_moves', 'checkGL_after_tmove'. If the build type is `debug`, the default value is 'all'. Otherwise, the default value is 'no'.
 
+- ``spin_mass`` Optional parameter to allow the user to change the rate of spin sampling. If spin sampling is on using ``spinor`` == yes in the electron ParticleSet input,  the spin mass determines the rate
+  of spin sampling, resulting in an effective spin timestep :math:`\tau_s = \frac{\tau}{\mu_s}`. The algorithm is described in detail in :cite:`Melton2016-1` and :cite:`Melton2016-2`.
+
 An example VMC section for a simple ``vmc_batch`` run:
 
 ::
@@ -1550,6 +1556,9 @@ Combining VMC and DMC in a single run (wavefunction optimization can be combined
   +--------------------------------+--------------+-------------------------+-------------+-----------------------------------------------+
   | ``debug_checks``               | text         | see additional info     | dep.        | Turn on/off additional recompute and checks   |
   +--------------------------------+--------------+-------------------------+-------------+-----------------------------------------------+
+  | ``spin_mass``                  | real         | :math:`\geq 0`          | 1.0         | Effective mass for spin sampling              |
+  +--------------------------------+--------------+-------------------------+-------------+-----------------------------------------------+
+
 
 - ``crowds`` The number of crowds that the walkers are subdivided into on each MPI rank. If not provided, it is set equal to the number of OpenMP threads.
 
@@ -1566,6 +1575,9 @@ Combining VMC and DMC in a single run (wavefunction optimization can be combined
 
 - ``debug_checks`` valid values are 'no', 'all', 'checkGL_after_load', 'checkGL_after_moves', 'checkGL_after_tmove'. If the build type is `debug`, the default value is 'all'. Otherwise, the default value is 'no'.
 
+- ``spin_mass`` Optional parameter to allow the user to change the rate of spin sampling. If spin sampling is on using ``spinor`` == yes in the electron ParticleSet input,  the spin mass determines the rate
+  of spin sampling, resulting in an effective spin timestep :math:`\tau_s = \frac{\tau}{\mu_s}`. The algorithm is described in detail in :cite:`Melton2016-1` and :cite:`Melton2016-2`.
+
 .. code-block::
   :caption: The following is an example of a minimal DMC section using the ``dmc_batch`` driver
   :name: Listing 48b

src/Particle/MCCoords.cpp

@@ -2,7 +2,7 @@
 // This file is distributed under the University of Illinois/NCSA Open Source License.
 // See LICENSE file in top directory for details.
 //
-// Copyright (c) 2022 developers.
+// Copyright (c) 2022 QMCPACK developers.
 //
 // File developed by: Peter Doak, [email protected], Oak Ridge National Laboratory
 //
@@ -13,12 +13,13 @@
 
 namespace qmcplusplus
 {
-template<CoordsType MCT>
-void MCCoords<MCT>::resize(const std::size_t size)
+void MCCoords<CoordsType::POS>::resize(const std::size_t size) { positions.resize(size); }
+
+void MCCoords<CoordsType::POS_SPIN>::resize(const std::size_t size)
 {
   positions.resize(size);
+  spins.resize(size);
 }
-
 template struct MCCoords<CoordsType::POS>;
 template struct MCCoords<CoordsType::POS_SPIN>;
 } // namespace qmcplusplus

src/Particle/MCCoords.hpp

@@ -2,7 +2,7 @@
 // This file is distributed under the University of Illinois/NCSA Open Source License.
 // See LICENSE file in top directory for details.
 //
-// Copyright (c) 2022 developers.
+// Copyright (c) 2022 QMCPACK developers.
 //
 // File developed by: Peter Doak, [email protected], Oak Ridge National Laboratory
 //
@@ -19,75 +19,32 @@
 
 namespace qmcplusplus
 {
-
 enum class CoordsType
 {
   POS,
   POS_SPIN
 };
 
-template<CoordsType MCT = CoordsType::POS>
-struct MCCoords
+template<CoordsType MCT>
+struct MCCoords;
+
+template<>
+struct MCCoords<CoordsType::POS>
 {
-  // This cleans up some other code.
   void resize(const std::size_t size);
+
   std::vector<QMCTraits::PosType> positions;
 };
 
 template<>
 struct MCCoords<CoordsType::POS_SPIN>
 {
-  // This cleans up some other code.
-  void resize(const std::size_t size)
-  {
-    positions.resize(size);
-    spins.resize(size);
-  }
+  void resize(const std::size_t size);
+
   std::vector<QMCTraits::PosType> positions;
   std::vector<QMCTraits::FullPrecRealType> spins;
 };
 
-/** Object to encapsulate appropriate tau derived values
- *  for a particular MCCoords specialization
- */
-template<typename Real, CoordsType CT = CoordsType::POS>
-struct Taus
-{
-  Real tauovermass;
-  Real oneover2tau;
-  Real sqrttau;
-  Taus(Real tau, Real grp_inv_mass)
-  {
-    tauovermass = tau * grp_inv_mass;
-    oneover2tau = 0.5 / (tauovermass);
-    sqrttau     = std::sqrt(tauovermass);
-  }
-};
-
-template<typename Real>
-struct Taus<Real, CoordsType::POS_SPIN> : public Taus<Real, CoordsType::POS>
-{
-  using Base = Taus<Real, CoordsType::POS>;
-  Real spin_tauovermass;
-  Real spin_oneover2tau;
-  Real spin_sqrttau;
-  Taus(Real tau, Real grp_inv_mass, Real spin_mass) : Base(tau, grp_inv_mass)
-  {
-    spin_tauovermass = Base::tauovermass / spin_mass;
-    spin_oneover2tau = 0.5 / (spin_tauovermass);
-    spin_sqrttau     = std::sqrt(spin_tauovermass);
-  }
-};
-
-/** Factory function for Taus based on MCCoordsTypes
- *  Note as in previous code value of tau derived values is not full precision.
- */
-template<CoordsType CT, typename... ARGS>
-auto makeTaus(MCCoords<CT>& mc_coords, const ARGS&... args)
-{
-  return Taus<QMCTraits::RealType, CT>(args...);
-}
-
 extern template struct MCCoords<CoordsType::POS>;
 extern template struct MCCoords<CoordsType::POS_SPIN>;
 } // namespace qmcplusplus

src/Particle/PSdispatcher.cpp

@@ -42,36 +42,29 @@ void PSdispatcher::flex_update(const RefVectorWithLeader<ParticleSet>& p_list, b
       pset.update(skipSK);
 }
 
+template<CoordsType CT>
 void PSdispatcher::flex_makeMove(const RefVectorWithLeader<ParticleSet>& p_list,
                                  int iat,
-                                 const std::vector<SingleParticlePos>& displs) const
+                                 const MCCoords<CT>& displs) const
 {
   if (use_batch_)
     ParticleSet::mw_makeMove(p_list, iat, displs);
   else
     for (size_t iw = 0; iw < p_list.size(); iw++)
-      p_list[iw].makeMove(iat, displs[iw]);
-}
-
-void PSdispatcher::flex_makeMoveWithSpin(const RefVectorWithLeader<ParticleSet>& p_list,
-                                         int iat,
-                                         const std::vector<SingleParticlePos>& displs,
-                                         const std::vector<Scalar_t>& sdispls) const
-{
-  if (use_batch_)
-    ParticleSet::mw_makeMoveWithSpin(p_list, iat, displs, sdispls);
-  else
-    for (size_t iw = 0; iw < p_list.size(); iw++)
-      p_list[iw].makeMoveWithSpin(iat, displs[iw], sdispls[iw]);
+      if constexpr (CT == CoordsType::POS_SPIN)
+        p_list[iw].makeMoveWithSpin(iat, displs.positions[iw], displs.spins[iw]);
+      else
+        p_list[iw].makeMove(iat, displs.positions[iw]);
 }
 
+template<CoordsType CT>
 void PSdispatcher::flex_accept_rejectMove(const RefVectorWithLeader<ParticleSet>& p_list,
                                           int iat,
                                           const std::vector<bool>& isAccepted,
                                           bool forward_mode) const
 {
   if (use_batch_)
-    ParticleSet::mw_accept_rejectMove(p_list, iat, isAccepted, forward_mode);
+    ParticleSet::mw_accept_rejectMove<CT>(p_list, iat, isAccepted, forward_mode);
   else
     for (size_t iw = 0; iw < p_list.size(); iw++)
       p_list[iw].accept_rejectMove(iat, isAccepted[iw], forward_mode);
@@ -96,4 +89,18 @@ void PSdispatcher::flex_saveWalker(const RefVectorWithLeader<ParticleSet>& p_lis
       p_list[iw].saveWalker(walkers[iw]);
 }
 
+template void PSdispatcher::flex_makeMove<CoordsType::POS>(const RefVectorWithLeader<ParticleSet>& p_list,
+                                                           int iat,
+                                                           const MCCoords<CoordsType::POS>& displs) const;
+template void PSdispatcher::flex_makeMove<CoordsType::POS_SPIN>(const RefVectorWithLeader<ParticleSet>& p_list,
+                                                                int iat,
+                                                                const MCCoords<CoordsType::POS_SPIN>& displs) const;
+template void PSdispatcher::flex_accept_rejectMove<CoordsType::POS>(const RefVectorWithLeader<ParticleSet>& p_list,
+                                                                    int iat,
+                                                                    const std::vector<bool>& isAccepted,
+                                                                    bool forward_mode) const;
+template void PSdispatcher::flex_accept_rejectMove<CoordsType::POS_SPIN>(const RefVectorWithLeader<ParticleSet>& p_list,
+                                                                         int iat,
+                                                                         const std::vector<bool>& isAccepted,
+                                                                         bool forward_mode) const;
 } // namespace qmcplusplus

src/Particle/PSdispatcher.h

@@ -37,15 +37,10 @@ class PSdispatcher
 
   void flex_update(const RefVectorWithLeader<ParticleSet>& p_list, bool skipSK = false) const;
 
-  void flex_makeMove(const RefVectorWithLeader<ParticleSet>& p_list,
-                     int iat,
-                     const std::vector<SingleParticlePos>& displs) const;
-
-  void flex_makeMoveWithSpin(const RefVectorWithLeader<ParticleSet>& p_list,
-                             int iat,
-                             const std::vector<SingleParticlePos>& displs,
-                             const std::vector<Scalar_t>& sdispls) const;
+  template<CoordsType CT>
+  void flex_makeMove(const RefVectorWithLeader<ParticleSet>& p_list, int iat, const MCCoords<CT>& displs) const;
 
+  template<CoordsType CT>
   void flex_accept_rejectMove(const RefVectorWithLeader<ParticleSet>& p_list,
                               int iat,
                               const std::vector<bool>& isAccepted,

src/Particle/ParticleBase/RandomSeqGenerator.h

@@ -79,14 +79,6 @@ inline void makeGaussRandomWithEngine(std::vector<TinyVector<T, D>>& a, RG& rng)
   assignGaussRand(&(a[0][0]), a.size() * D, rng);
 }
 
-template<CoordsType CT, class RG>
-inline void makeGaussRandomWithEngine(MCCoords<CT>& a, RG& rng)
-{
-  makeGaussRandomWithEngine(a.positions, rng);
-  if constexpr (std::is_same<MCCoords<CT>, MCCoords<CoordsType::POS_SPIN>>::value)
-    makeGaussRandomWithEngine(a.spins, rng);
-}
-
 template<typename T, class RG>
 inline void makeGaussRandomWithEngine(std::vector<T>& a, RG& rng)
 {
@@ -101,5 +93,13 @@ inline void makeGaussRandomWithEngine(ParticleAttrib<T>& a, RG& rng)
   assignGaussRand(&(a[0]), a.size(), rng);
 }
 
+template<CoordsType CT, class RG>
+inline void makeGaussRandomWithEngine(MCCoords<CT>& a, RG& rng)
+{
+  makeGaussRandomWithEngine(a.positions, rng);
+  if constexpr (CT == CoordsType::POS_SPIN)
+    makeGaussRandomWithEngine(a.spins, rng);
+}
+
 } // namespace qmcplusplus
 #endif