Merge branch 'develop' into develop

.gitignore

@@ -20,3 +20,7 @@ tests/solids/NiO_a4_e48_pp/NiO-fcc-supertwist111-supershift000-S1.h5
 
 # Visual Studio code settings
 .vscode/
+
+# Doxygen
+doxygen/output/
+doxygen/qmcpack_all.tag

docs/methods.rst

@@ -151,8 +151,6 @@ Variational Monte Carlo
   +--------------------------------+--------------+-------------------------+-------------+-----------------------------------------------+
   | ``blocks_between_recompute``   | integer      | :math:`\geq 0`          | dep.        | Wavefunction recompute frequency              |
   +--------------------------------+--------------+-------------------------+-------------+-----------------------------------------------+
-  | ``spinMoves``                  | text         | yes,no                  | no          | Whether or not to sample the electron spins   |
-  +--------------------------------+--------------+-------------------------+-------------+-----------------------------------------------+
   | ``spinMass``                   | real         | :math:`> 0`             | 1.0         | Effective mass for spin sampling              |
   +--------------------------------+--------------+-------------------------+-------------+-----------------------------------------------+
 
@@ -228,12 +226,8 @@ Additional information:
   recompute) by default when not using mixed precision. Recomputing
   introduces a performance penalty dependent on system size.
 
-- ``spinMoves`` Determines whether or not the spin variables are sampled following
+- ``spinMass`` 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
-  :cite:`Melton2016-1` and :cite:`Melton2016-2`. If a relativistic calculation is desired using pseudopotentials,
+  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`.
-  spin variable sampling is required.
-
-- ``spinMass`` If spin sampling is on using ``spinMoves`` == yes, the spin mass determines the rate
-  of spin sampling, resulting in an effective spin timestep :math:`\tau_s = \frac{\tau}{\mu_s}`.
 
 An example VMC section for a simple VMC run:
 
@@ -1189,8 +1183,6 @@ parameters:
   +--------------------------------+--------------+-------------------------+-------------+-----------------------------------------------+
   | ``blocks_between_recompute``   | integer      | :math:`\geq 0`          | dep.        | Wavefunction recompute frequency              |
   +--------------------------------+--------------+-------------------------+-------------+-----------------------------------------------+
-  | ``spinMoves``                  | text         | yes,no                  | no          | Whether or not to sample the electron spins   |
-  +--------------------------------+--------------+-------------------------+-------------+-----------------------------------------------+
   | ``spinMass``                   | real         | :math:`> 0`             | 1.0         | Effective mass for spin sampling              |
   +--------------------------------+--------------+-------------------------+-------------+-----------------------------------------------+
 
@@ -1280,12 +1272,9 @@ Additional information:
    elapsed, the program will finalize the simulation even if all blocks
    are not completed.
 
--  ``spinMoves`` Determines whether or not the spin variables are sampled following :cite:`Melton2016-1` 
+-  ``spinMass`` This is an 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 
-   and :cite:`Melton2016-2`. If a relativistic calculation is desired using pseudopotentials, spin variable sampling is required.
-
--  ``spinMass`` If spin sampling is on using ``spinMoves`` == yes, the spin mass determines the rate 
    of spin sampling, resulting in an effective spin timestep :math:`\tau_s = \frac{\tau}{\mu_s}` where 
-   :math:`\tau` is the normal spatial timestep and :math:`\mu_s` is the value of the spin mass.
+   :math:`\tau` is the normal spatial timestep and :math:`\mu_s` is the value of the spin mass. The algorithm is described in detail in :cite:`Melton2016-1` and :cite:`Melton2016-2`.
 
 
 -  ``energyUpdateInterval``: The default is to update the trial energy

docs/simulationcell.rst

@@ -160,17 +160,19 @@ Input specification
 
 Attribute:
 
-  +----------------------------------------+----------+----------------------+---------+------------------------------+
+  +----------------------------------------+----------+----------------------+---------+-------------------------------+
   | Name                                   | Datatype | Values               | Default | Description                   |
-  +========================================+==========+======================+=========+==============================+
+  +========================================+==========+======================+=========+===============================+
   | ``name/id``                            | Text     | *Any*                | e       | Name of particle set          |
-  +----------------------------------------+----------+----------------------+---------+------------------------------+
+  +----------------------------------------+----------+----------------------+---------+-------------------------------+
   | ``size``:math:`^o`                     | Integer  | *Any*                | 0       | Number of particles in set    |
-  +----------------------------------------+----------+----------------------+---------+------------------------------+
+  +----------------------------------------+----------+----------------------+---------+-------------------------------+
   | ``random``\ :math:`^o`                 | Text     | Yes/no               | No      | Randomize starting positions  |
-  +----------------------------------------+----------+----------------------+---------+------------------------------+
+  +----------------------------------------+----------+----------------------+---------+-------------------------------+
   | ``randomsrc``/``randomsrc``:math:`^o`  | Text     | ``particleset.name`` | *None*  | Particle set to randomize     |
-  +----------------------------------------+----------+----------------------+---------+------------------------------+
+  +----------------------------------------+----------+----------------------+---------+-------------------------------+
+  | ``spinor``:math:`^o`                   | Text     | Yes/no               | No      | particleset treated as spinor |
+  +----------------------------------------+----------+----------------------+---------+-------------------------------+
 
 Detailed attribute description
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -251,6 +253,11 @@ Optional particleset attributes
    | Specify source particle set around which to randomize the initial
      positions of this particle set.
 
+-  | ``spinor``
+   | Sets an internal flag that the particleset (usually for electrons) is
+     a spinor object. This is used in the wavefunction builders and QMC drivers
+     to determiane if spin sampling will be used
+
 Required name attributes
 ^^^^^^^^^^^^^^^^^^^^^^^^
 

docs/spin_orbit.rst

@@ -58,7 +58,7 @@ where we now utilize determinants of spinors, as opposed to the usual product of
   <?xml version="1.0"?>
   <qmcsystem>
     <wavefunction name="psi0" target="e">
-      <sposet_builder name="spo_builder" type="spinorbspline" href="eshdf.h5" tilematrix="100010001" twistnum="0" source="ion0" size="10">
+      <sposet_builder name="spo_builder" type="bspline" href="eshdf.h5" tilematrix="100010001" twistnum="0" source="ion0" size="10">
         <sposet type="bspline" name="myspo" size="10">
           <occupation mode="ground"/>
         </sposet>
@@ -95,7 +95,7 @@ We also make a small modification in the particleset specification:
   :caption: specification for the electron particle when performing spin-orbit calculations
   :name: slisting2
 
-  <particleset name="e" random="yes" randomsrc="ion0">
+  <particleset name="e" random="yes" randomsrc="ion0" spinor="yes">
      <group name="u" size="10" mass="1.0">
         <parameter name="charge"              >    -1                    </parameter>
         <parameter name="mass"                >    1.0                   </parameter>
@@ -103,6 +103,8 @@ We also make a small modification in the particleset specification:
   </particleset>
 
 Note that we only provide a single electron group to represent all electrons  in the system, as opposed to the usual separation of up and down electrons. 
+The additional keyword ``spinor=yes`` is the *only* required keyword for spinors.
+This will be used internally to determine which movers to use in QMC drivers (e.g. VMCUpdatePbyP vs SOVMCUpdatePbyP) and which SPOSets to use in the trial wave function (spinors vs. normal orbitals)
 
 *note*: In the current implementation, spinor wavefunctions are only 
 supported at the single determinant level. Multideterminant spinor wave functions will be supported in a future release. 
@@ -133,8 +135,7 @@ Green's function for the spatial evolution and the *spin kinetic energy*
 operator introduces a Green's function for the spin variables. 
 Note that this includes a contribution from the *spin drift* :math:`\mathbf{v}_{\mathbf{S}}(\mathbf{S}) =  \nabla_{\mathbf{S}} \ln \Psi_T(\mathbf{S})`.
 
-In both the VMC and DMC methods, the spin sampling is controlled by two input
+In both the VMC and DMC methods, there are no required changes to a typical input
-parameters in the ``xml`` blocks. 
 
 .. code-block::
 
@@ -143,14 +144,17 @@ parameters in the ``xml`` blocks.
     <parameter name="blocks"   >    50   </parameter>
     <parameter name="walkers"  >    10   </parameter>
     <parameter name="timestep" >  0.01   </parameter>
-    <parameter name="spinMoves">   yes   </parameter>
-    <parameter name="spinMass" >   1.0   </parameter>
   </qmc>
 
-The ``spinMoves`` flag turns on the spin sampling, which is off by default. 
+Whether or not spin moves are used is determined internally by the ``spinor=yes`` flag in particleset.
-The ``spinMass`` flag sets the :math:`\mu_s` parameter used in the 
+
-particle updates, and effectively controls the rate of sampling for the spin 
+By default, the spin mass :math:`\mu_s` (which controls the rate of spin sampling relative to the spatial sampling) is set to 1.0. 
-coordinates relative to the spatial coordinates. 
+This can be changed by adding an additional parameter to the QMC input
+
+.. code-block::
+
+ <parameter name="spinMass" > 0.25 </parameter>
+
 A larger/smaller spin mass corresponds to slower/faster spin sampling relative to the spatial coordinates.
 
 Spin-Orbit Effective Core Potentials

external_codes/catch/complex_approx.hpp

@@ -7,6 +7,7 @@
 // Copy and modify the Approx class to handle complex numbers
 
 namespace Catch {
+namespace Detail {
 class ComplexApprox
 {
     std::complex<double> m_value;
@@ -87,20 +88,22 @@ public:
        return os;
     }
 };
+}
   
 template<>
-struct StringMaker<ComplexApprox> {
+struct StringMaker<Catch::Detail::ComplexApprox> {
-  static std::string convert(ComplexApprox const &value);
+  static std::string convert(Catch::Detail::ComplexApprox const &value);
 };
 
 #ifdef CATCH_IMPL
-std::string StringMaker<ComplexApprox>::convert(ComplexApprox const& value)
+std::string StringMaker<Catch::Detail::ComplexApprox>::convert(Catch::Detail::ComplexApprox const& value)
 {
   return value.toString();
 }
 #endif
+
 }
 
-using Catch::ComplexApprox;
+using Catch::Detail::ComplexApprox;
 
 #endif

src/Configuration.h

@@ -101,9 +101,9 @@ struct PtclOnLatticeTraits
 //  Check if we are compiling with Catch defined.  Could use other symbols if needed.
 #ifdef TEST_CASE
 #ifdef QMC_COMPLEX
-typedef ComplexApprox ValueApprox;
+  using ValueApprox = Catch::Detail::ComplexApprox;
 #else
-typedef Approx ValueApprox;
+  using ValueApprox = Catch::Detail::Approx;
 #endif
 #endif
 

src/Estimators/tests/test_local_energy_est.cpp

@@ -48,9 +48,9 @@ TEST_CASE("LocalEnergy", "[estimators]")
   QMCHamiltonian H;
   LocalEnergyEstimator le_est(H, false);
 
-  LocalEnergyEstimator* le_est2 = dynamic_cast<LocalEnergyEstimator*>(le_est.clone());
+  std::unique_ptr<LocalEnergyEstimator> le_est2{dynamic_cast<LocalEnergyEstimator*>(le_est.clone())};
-  REQUIRE(le_est2 != NULL);
+  REQUIRE(le_est2 != nullptr);
-  REQUIRE(le_est2 != &le_est);
+  REQUIRE(le_est2.get() != &le_est);
 
   MCWalkerConfiguration W;
   W.setName("electrons");

src/Estimators/tests/test_trace_manager.cpp

@@ -33,8 +33,6 @@ TEST_CASE("TraceManager", "[estimators]")
 {
   Communicate* c = OHMMS::Controller;
 
-  Array<TraceInt, 1>* int_sample;
-
   TraceManager tm(c);
 
   tm.put(NULL, true, "test");
@@ -47,7 +45,7 @@ TEST_CASE("TraceManager", "[estimators]")
 
   tm.request.determine_stream_write();
 
-  int_sample       = tm.checkout_int<1>("scalar1");
+  auto int_sample  = std::unique_ptr<Array<TraceInt, 1>>{tm.checkout_int<1>("scalar1")};
   (*int_sample)(0) = 2;
 
   tm.update_status();
@@ -85,88 +83,71 @@ TEST_CASE("TraceManager check_trace_build", "[estimators]")
   TraceBuffer<int> tbi;
   TraceBuffer<TraceReal> tbr;
 
-  Array<int, 1>* ai1;
+  auto ai1 = std::unique_ptr<Array<int, 1>>{tssi.checkout_array<1>(domain, name1, shape1)};
-  Array<int, 2>* ai2;
+  ai1.reset(tssi.checkout_array<1>(P, name, shape1));
-  Array<int, 3>* ai3;
+  auto ai2 = std::unique_ptr<Array<int, 2>>{tssi.checkout_array<2>(domain, name2, shape2)};
-  Array<int, 4>* ai4;
+  ai2.reset(tssi.checkout_array<2>(P, name2, shape2));
-  ai1 = tssi.checkout_array<1>(domain, name1, shape1);
+  auto ai3 = std::unique_ptr<Array<int, 3>>{tssi.checkout_array<3>(domain, name3, shape3)};
-  ai1 = tssi.checkout_array<1>(P, name, shape1);
+  ai3.reset(tssi.checkout_array<3>(P, name3, shape3));
-  ai2 = tssi.checkout_array<2>(domain, name2, shape2);
+  auto ai4 = std::unique_ptr<Array<int, 4>>{tssi.checkout_array<4>(domain, name4, shape4)};
-  ai2 = tssi.checkout_array<2>(P, name2, shape2);
+  ai4.reset(tssi.checkout_array<4>(P, name4, shape4));
-  ai3 = tssi.checkout_array<3>(domain, name3, shape3);
-  ai3 = tssi.checkout_array<3>(P, name3, shape3);
-  ai4 = tssi.checkout_array<4>(domain, name4, shape4);
-  ai4 = tssi.checkout_array<4>(P, name4, shape4);
 
-  Array<TraceReal, 1>* ar1;
+  auto ar1 = std::unique_ptr<Array<TraceReal, 1>>{tssr.checkout_array<1>(domain, name1, shape1)};
-  Array<TraceReal, 2>* ar2;
+  ar1.reset(tssr.checkout_array<1>(P, name1, shape1));
-  Array<TraceReal, 3>* ar3;
+  auto ar2 = std::unique_ptr<Array<TraceReal, 2>>{tssr.checkout_array<2>(domain, name2, shape2)};
-  Array<TraceReal, 4>* ar4;
+  ar2.reset(tssr.checkout_array<2>(P, name2, shape2));
-  ar1 = tssr.checkout_array<1>(domain, name1, shape1);
+  auto ar3 = std::unique_ptr<Array<TraceReal, 3>>{tssr.checkout_array<3>(domain, name3, shape3)};
-  ar1 = tssr.checkout_array<1>(P, name1, shape1);
+  ar3.reset(tssr.checkout_array<3>(P, name3, shape3));
-  ar2 = tssr.checkout_array<2>(domain, name2, shape2);
+  auto ar4 = std::unique_ptr<Array<TraceReal, 4>>{tssr.checkout_array<4>(domain, name4, shape4)};
-  ar2 = tssr.checkout_array<2>(P, name2, shape2);
+  ar4.reset(tssr.checkout_array<4>(P, name4, shape4));
-  ar3 = tssr.checkout_array<3>(domain, name3, shape3);
-  ar3 = tssr.checkout_array<3>(P, name3, shape3);
-  ar4 = tssr.checkout_array<4>(domain, name4, shape4);
-  ar4 = tssr.checkout_array<4>(P, name4, shape4);
 
-
+  auto ac1 = std::unique_ptr<Array<TraceComp, 1>>{tssc.checkout_array<1>(domain, name1, shape1)};
-  Array<TraceComp, 1>* ac1;
+  ac1.reset(tssc.checkout_array<1>(P, name1, shape1));
-  Array<TraceComp, 2>* ac2;
+  auto ac2 = std::unique_ptr<Array<TraceComp, 2>>{tssc.checkout_array<2>(domain, name2, shape2)};
-  Array<TraceComp, 3>* ac3;
+  ac2.reset(tssc.checkout_array<2>(P, name2, shape2));
-  Array<TraceComp, 4>* ac4;
+  auto ac3 = std::unique_ptr<Array<TraceComp, 3>>{tssc.checkout_array<3>(domain, name3, shape3)};
-  ac1 = tssc.checkout_array<1>(domain, name1, shape1);
+  ac3.reset(tssc.checkout_array<3>(P, name3, shape3));
-  ac1 = tssc.checkout_array<1>(P, name1, shape1);
+  auto ac4 = std::unique_ptr<Array<TraceComp, 4>>{tssc.checkout_array<4>(domain, name4, shape4)};
-  ac2 = tssc.checkout_array<2>(domain, name2, shape2);
+  ac4.reset(tssc.checkout_array<4>(P, name4, shape4));
-  ac2 = tssc.checkout_array<2>(P, name2, shape2);
-  ac3 = tssc.checkout_array<3>(domain, name3, shape3);
-  ac3 = tssc.checkout_array<3>(P, name3, shape3);
-  ac4 = tssc.checkout_array<4>(domain, name4, shape4);
-  ac4 = tssc.checkout_array<4>(P, name4, shape4);
 
   TraceManager tm;
-  Array<long, 1>* al1;
+  auto al1 = std::unique_ptr<Array<long, 1>>{tm.checkout_int<1>(name1)};
-  Array<long, 2>* al2;
+  al1.reset(tm.checkout_int<1>(domain, name1, 10));
-  Array<long, 3>* al3;
+  al1.reset(tm.checkout_int<1>(name1, P));
-  Array<long, 4>* al4;
+  auto al2 = std::unique_ptr<Array<long, 2>>{tm.checkout_int<2>(name2, 5, 6)};
-  al1 = tm.checkout_int<1>(name1);
+  al2.reset(tm.checkout_int<2>(domain, name2, 10, 11));
-  al1 = tm.checkout_int<1>(domain, name1, 10);
+  al2.reset(tm.checkout_int<2>(name2, P, 11));
-  al1 = tm.checkout_int<1>(name1, P);
+  auto al3 = std::unique_ptr<Array<long, 3>>{tm.checkout_int<3>(name3, 5, 6, 7)};
-  al2 = tm.checkout_int<2>(name2, 5, 6);
+  al3.reset(tm.checkout_int<3>(domain, name3, 10, 11, 12));
-  al2 = tm.checkout_int<2>(domain, name2, 10, 11);
+  al3.reset(tm.checkout_int<3>(name3, P, 11, 12));
-  al2 = tm.checkout_int<2>(name2, P, 11);
+  auto al4 = std::unique_ptr<Array<long, 4>>{tm.checkout_int<4>(name4, 5, 6, 7, 8)};
-  al3 = tm.checkout_int<3>(name3, 5, 6, 7);
+  al4.reset(tm.checkout_int<4>(domain, name4, 10, 11, 12, 13));
-  al3 = tm.checkout_int<3>(domain, name3, 10, 11, 12);
+  al4.reset(tm.checkout_int<4>(name4, P, 11, 12, 13));
-  al3 = tm.checkout_int<3>(name3, P, 11, 12);
+  ar1.reset(tm.checkout_real<1>(name1));
-  al4 = tm.checkout_int<4>(name4, 5, 6, 7, 8);
+  ar1.reset(tm.checkout_real<1>(domain, name1, 10));
-  al4 = tm.checkout_int<4>(domain, name4, 10, 11, 12, 13);
+  ar1.reset(tm.checkout_real<1>(name1, P));
-  al4 = tm.checkout_int<4>(name4, P, 11, 12, 13);
+  ar2.reset(tm.checkout_real<2>(name2, 5, 6));
-  ar1 = tm.checkout_real<1>(name1);
+  ar2.reset(tm.checkout_real<2>(domain, name2, 10, 11));
-  ar1 = tm.checkout_real<1>(domain, name1, 10);
+  ar2.reset(tm.checkout_real<2>(name2, P, 11));
-  ar1 = tm.checkout_real<1>(name1, P);
+  ar3.reset(tm.checkout_real<3>(name3, 5, 6, 7));
-  ar2 = tm.checkout_real<2>(name2, 5, 6);
+  ar3.reset(tm.checkout_real<3>(domain, name3, 10, 11, 12));
-  ar2 = tm.checkout_real<2>(domain, name2, 10, 11);
+  ar3.reset(tm.checkout_real<3>(name3, P, 11, 12));
-  ar2 = tm.checkout_real<2>(name2, P, 11);
+  ar4.reset(tm.checkout_real<4>(name4, 5, 6, 7, 8));
-  ar3 = tm.checkout_real<3>(name3, 5, 6, 7);
+  ar4.reset(tm.checkout_real<4>(domain, name4, 10, 11, 12, 13));
-  ar3 = tm.checkout_real<3>(domain, name3, 10, 11, 12);
+  ar4.reset(tm.checkout_real<4>(name4, P, 11, 12, 13));
-  ar3 = tm.checkout_real<3>(name3, P, 11, 12);
-  ar4 = tm.checkout_real<4>(name4, 5, 6, 7, 8);
-  ar4 = tm.checkout_real<4>(domain, name4, 10, 11, 12, 13);
-  ar4 = tm.checkout_real<4>(name4, P, 11, 12, 13);
 
-  ac1 = tm.checkout_complex<1>(name1);
+  ac1.reset(tm.checkout_complex<1>(name1));
-  ac1 = tm.checkout_complex<1>(domain, name1, 10);
+  ac1.reset(tm.checkout_complex<1>(domain, name1, 10));
-  ac1 = tm.checkout_complex<1>(name1, P);
+  ac1.reset(tm.checkout_complex<1>(name1, P));
-  ac2 = tm.checkout_complex<2>(name2, 5, 6);
+  ac2.reset(tm.checkout_complex<2>(name2, 5, 6));
-  ac2 = tm.checkout_complex<2>(domain, name2, 10, 11);
+  ac2.reset(tm.checkout_complex<2>(domain, name2, 10, 11));
-  ac2 = tm.checkout_complex<2>(name2, P, 11);
+  ac2.reset(tm.checkout_complex<2>(name2, P, 11));
-  ac3 = tm.checkout_complex<3>(name3, 5, 6, 7);
+  ac3.reset(tm.checkout_complex<3>(name3, 5, 6, 7));
-  ac3 = tm.checkout_complex<3>(domain, name3, 10, 11, 12);
+  ac3.reset(tm.checkout_complex<3>(domain, name3, 10, 11, 12));
-  ac3 = tm.checkout_complex<3>(name3, P, 11, 12);
+  ac3.reset(tm.checkout_complex<3>(name3, P, 11, 12));
-  ac4 = tm.checkout_complex<4>(name4, 5, 6, 7, 8);
+  ac4.reset(tm.checkout_complex<4>(name4, 5, 6, 7, 8));
-  ac4 = tm.checkout_complex<4>(domain, name4, 10, 11, 12, 13);
+  ac4.reset(tm.checkout_complex<4>(domain, name4, 10, 11, 12, 13));
-  ac4 = tm.checkout_complex<4>(name4, P, 11, 12, 13);
+  ac4.reset(tm.checkout_complex<4>(name4, P, 11, 12, 13));
 }
 
 } // namespace qmcplusplus

src/Particle/LongRange/LRHandlerBase.h

@@ -69,10 +69,10 @@ struct LRHandlerBase
   //return k cutoff
   inline mRealType get_kc() const { return LR_kc; }
 
-  /** evaluate \f$\sum_k F_{k} \rho^1_{-{\bf k} \rho^2_{\bf k}\f$
+  /** evaluate \f$\sum_k F_{k} \rho^1_{-{\bf k}} \rho^2_{\bf k}\f$
    * @param kshell degeneracies of the vectors
-   * @param rk1 starting address of \f$\rho^1_{{\bf k}\f$
+   * @param rk1 starting address of \f$\rho^1_{{\bf k}}\f$
-   * @param rk2 starting address of \f$\rho^2_{{\bf k}\f$
+   * @param rk2 starting address of \f$\rho^2_{{\bf k}}\f$
    *
    * Valid for the strictly ordered k and \f$F_{k}\f$.
    */
@@ -169,11 +169,11 @@ struct LRHandlerBase
     return vk;
   }
 
-  /** evaluate \f$\sum_k F_{k} \rho^1_{-{\bf k} \rho^2_{\bf k}\f$
+  /** evaluate \f$\sum_k F_{k} \rho^1_{-{\bf k}} \rho^2_{\bf k}\f$
-   * and \f$\sum_k F_{k} \rho^1_{-{\bf k} \rho^2_{\bf k}\f$
+   * and \f$\sum_k F_{k} \rho^1_{-{\bf k}} \rho^2_{\bf k}\f$
    * @param kshell degeneracies of the vectors
-   * @param rk1 starting address of \f$\rho^1_{{\bf k}\f$
+   * @param rk1 starting address of \f$\rho^1_{{\bf k}}\f$
-   * @param rk2 starting address of \f$\rho^2_{{\bf k}\f$
+   * @param rk2 starting address of \f$\rho^2_{{\bf k}}\f$
    *
    * Valid for the strictly ordered k and \f$F_{k}\f$.
    */

src/Particle/LongRange/LRRPABFeeHandlerTemp.h

@@ -132,10 +132,10 @@ struct LRRPABFeeHandlerTemp : public LRHandlerBase
     //       for(int n=0; n<coefs.size(); n++) v -= coefs[n]*Basis.h(n,r);
   }
 
-  /** evaluate \f$\sum_k F_{k} \rho^1_{-{\bf k} \rho^2_{\bf k}\f$
+  /** evaluate \f$\sum_k F_{k} \rho^1_{-{\bf k}} \rho^2_{\bf k}\f$
    * @param kshell degeneracies of the vectors
-   * @param rk1 starting address of \f$\rho^1_{{\bf k}\f$
+   * @param rk1 starting address of \f$\rho^1_{{\bf k}}\f$
-   * @param rk2 starting address of \f$\rho^2_{{\bf k}\f$
+   * @param rk2 starting address of \f$\rho^2_{{\bf k}}\f$
    *
    * Valid for the strictly ordered k and \f$F_{k}\f$.
    */

src/Particle/LongRange/LRRPAHandlerTemp.h

@@ -130,10 +130,10 @@ struct LRRPAHandlerTemp : public LRHandlerBase
     //       for(int n=0; n<coefs.size(); n++) v -= coefs[n]*Basis.h(n,r);
   }
 
-  /** evaluate \f$\sum_k F_{k} \rho^1_{-{\bf k} \rho^2_{\bf k}\f$
+  /** evaluate \f$\sum_k F_{k} \rho^1_{-{\bf k}} \rho^2_{\bf k}\f$
    * @param kshell degeneracies of the vectors
-   * @param rk1 starting address of \f$\rho^1_{{\bf k}\f$
+   * @param rk1 starting address of \f$\rho^1_{{\bf k}}\f$
-   * @param rk2 starting address of \f$\rho^2_{{\bf k}\f$
+   * @param rk2 starting address of \f$\rho^2_{{\bf k}}\f$
    *
    * Valid for the strictly ordered k and \f$F_{k}\f$.
    */

src/Particle/ParticleSet.cpp

@@ -68,6 +68,7 @@ ParticleSet::ParticleSet(const DynamicCoordinateKind kind)
       IsGrouped(true),
       SameMass(true),
       ThreadID(0),
+      is_spinor_(false),
       activePtcl(-1),
       Properties(0, 0, 1, WP::MAXPROPERTIES),
       myTwist(0.0),
@@ -83,6 +84,7 @@ ParticleSet::ParticleSet(const ParticleSet& p)
     : IsGrouped(p.IsGrouped),
       SameMass(true),
       ThreadID(0),
+      is_spinor_(false),
       activePtcl(-1),
       mySpecies(p.getSpeciesSet()),
       Properties(p.Properties),

src/Particle/ParticleSet.h

@@ -104,6 +104,8 @@ public:
   bool SameMass;
   ///threa id
   Index_t ThreadID;
+  ///true is a dynamic spin calculation
+  bool is_spinor_;
   /** the index of the active particle during particle-by-particle moves
    *
    * when a single particle move is proposed, the particle id is assigned to activePtcl
@@ -562,6 +564,7 @@ public:
     spins            = ptclin.spins;
     ID               = ptclin.ID;
     GroupID          = ptclin.GroupID;
+    is_spinor_       = ptclin.is_spinor_;
     if (ptclin.SubPtcl.size())
     {
       SubPtcl.resize(ptclin.SubPtcl.size());

src/Particle/ParticleSetPool.cpp

@@ -157,6 +157,7 @@ bool ParticleSetPool::put(xmlNodePtr cur)
   std::string randomR("no");
   std::string randomsrc;
   std::string useGPU;
+  std::string spinor;
   OhmmsAttributeSet pAttrib;
   pAttrib.add(id, "id");
   pAttrib.add(id, "name");
@@ -164,6 +165,7 @@ bool ParticleSetPool::put(xmlNodePtr cur)
   pAttrib.add(randomR, "random");
   pAttrib.add(randomsrc, "randomsrc");
   pAttrib.add(randomsrc, "random_source");
+  pAttrib.add(spinor, "spinor", {"no", "yes"});
 #if defined(ENABLE_OFFLOAD)
   pAttrib.add(useGPU, "gpu", {"yes", "no"});
 #endif
@@ -206,6 +208,7 @@ bool ParticleSetPool::put(xmlNodePtr cur)
       randomize_nodes.push_back(anode);
     }
     pTemp->setName(id);
+    pTemp->is_spinor_ = spinor == "yes";
     app_summary() << "  Particle set size: " << pTemp->getTotalNum() << std::endl;
     app_summary() << std::endl;
     return success;

src/Particle/tests/test_particle_pool.cpp

@@ -91,7 +91,7 @@ TEST_CASE("ParticleSetPool random", "[qmcapp]")
     0.1 0.2 0.3 \
   </attrib> \
 </particleset> \
-<particleset name=\"elec\" random=\"yes\" randomsrc=\"ion0\"> \
+<particleset name=\"elec\" random=\"yes\" randomsrc=\"ion0\" spinor=\"yes\"> \
    <group name=\"u\" size=\"4\"> \
       <parameter name=\"charge\">-1</parameter> \
    </group> \
@@ -110,12 +110,15 @@ TEST_CASE("ParticleSetPool random", "[qmcapp]")
 
   ParticleSet* ions = pp.getParticleSet("ion0");
   REQUIRE(ions != NULL);
+  REQUIRE(!ions->is_spinor_);
 
   ParticleSet* elec = pp.getParticleSet("elec");
   REQUIRE(ions != NULL);
+  REQUIRE(elec->is_spinor_);
   REQUIRE(elec->R.size() == 4);
   REQUIRE(elec->spins.size() == 4);
 
+
   // should do something
   pp.randomize();
 

src/QMCDrivers/DMC/DMC.cpp

@@ -70,6 +70,7 @@ void DMC::resetUpdateEngines()
               "reconfiguration to \"runwhileincorrect\" instead of \"yes\" to restore consistent behaviour.")
   makeClones(W, Psi, H);
   Timer init_timer;
+  bool spinor = false;
   if (Movers.empty())
   {
     W.loadEnsemble(wClones);
@@ -108,8 +109,6 @@ void DMC::resetUpdateEngines()
         o << "\n  Walkers are killed when a node crossing is detected";
       else
         o << "\n  DMC moves are rejected when a node crossing is detected";
-      if (SpinMoves == "yes")
-        o << "\n  Spins treated as dynamic variable with SpinMass: " << SpinMass;
       app_log() << o.str() << std::endl;
     }
 #pragma omp parallel for
@@ -126,8 +125,9 @@ void DMC::resetUpdateEngines()
       Rng[ip] = new RandomGenerator_t(*RandomNumberControl::Children[ip]);
       hClones[ip]->setRandomGenerator(Rng[ip]);
 #endif
-      if (SpinMoves == "yes")
+      if (W.is_spinor_)
       {
+        spinor = true;
         if (qmc_driver_mode[QMC_UPDATE_MODE])
         {
           Movers[ip] = new SODMCUpdatePbyPWithRejectionFast(*wClones[ip], *psiClones[ip], *hClones[ip], *Rng[ip]);
@@ -177,6 +177,10 @@ void DMC::resetUpdateEngines()
     }
   }
 #endif
+
+  if (spinor)
+    app_log() << "   Spins treated as dynamic variable with SpinMass: " << SpinMass << std::endl;
+
   branchEngine->checkParameters(W);
   int mxage = mover_MaxAge;
   if (fixW)

src/QMCDrivers/QMCDriver.cpp

@@ -117,9 +117,7 @@ QMCDriver::QMCDriver(MCWalkerConfiguration& w,
 
   m_param.add(nTargetPopulation, "samples");
 
-  SpinMoves = "no";
   SpinMass = 1.0;
-  m_param.add(SpinMoves, "SpinMoves");
   m_param.add(SpinMass, "SpinMass");
 
   Tau = 0.1;
@@ -520,14 +518,6 @@ bool QMCDriver::putQMCInfo(xmlNodePtr cur)
   if (!AppendRun)
     CurrentStep = 0;
 
-  tolower(SpinMoves);
-  if (SpinMoves != "yes" && SpinMoves != "no")
-    myComm->barrier_and_abort("SpinMoves must be yes/no!\n");
-#if defined(QMC_CUDA)
-  if (SpinMoves == "yes")
-    myComm->barrier_and_abort("Spin moves are not supported in legacy CUDA build.");
-#endif
-
   //if walkers are initialized via <mcwalkerset/>, use the existing one
   if (qmc_common.qmc_counter || qmc_common.is_restart)
   {

src/QMCDrivers/QMCDriver.h

@@ -331,8 +331,7 @@ protected:
   ///temporary storage for random displacement
   ParticleSet::ParticlePos_t deltaR;
 
-  ///turn on spin moves
+  ///spin mass for spinor calcs
-  std::string SpinMoves;
   RealType SpinMass;
 
   bool putQMCInfo(xmlNodePtr cur);

src/QMCDrivers/VMC/VMC.cpp

@@ -156,6 +156,7 @@ void VMC::resetRun()
   app_log() << std::endl;
 
   bool movers_created = false;
+  bool spinors        = false;
   if (Movers.empty())
   {
     movers_created = true;
@@ -179,9 +180,9 @@ void VMC::resetRun()
       Rng[ip] = new RandomGenerator_t(*(RandomNumberControl::Children[ip]));
 #endif
       hClones[ip]->setRandomGenerator(Rng[ip]);
-
+      if (W.is_spinor_)
-      if (SpinMoves == "yes")
       {
+        spinors = true;
         if (qmc_driver_mode[QMC_UPDATE_MODE])
         {
           Movers[ip] = new SOVMCUpdatePbyP(*wClones[ip], *psiClones[ip], *hClones[ip], *Rng[ip]);
@@ -239,7 +240,7 @@ void VMC::resetRun()
       Movers[i]->UseDrift = false;
   }
 
-  if (SpinMoves == "yes")
+  if (spinors)
   {
     app_log() << "  Spins treated as dynamic variable with SpinMass: " << SpinMass << std::endl;
     for (int i = 0; i < Movers.size(); i++)
@@ -391,7 +392,6 @@ bool VMC::put(xmlNodePtr q)
   app_log() << "  target samples = " << nTargetPopulation << std::endl;
   app_log() << "  walkers/mpi    = " << W.getActiveWalkers() << std::endl << std::endl;
   app_log() << "  stepsbetweensamples = " << nStepsBetweenSamples << std::endl;
-  app_log() << "  SpinMoves      = " << SpinMoves << std::endl;
 
   m_param.get(app_log());
 

src/QMCDrivers/tests/test_dmc_driver.cpp

@@ -150,6 +150,7 @@ TEST_CASE("SODMC", "[drivers][dmc]")
   elec.R[0][1]    = 0.0;
   elec.R[0][2]    = 0.0;
   elec.spins[0]   = 0.0;
+  elec.is_spinor_ = true;
   elec.createWalkers(1);
 
   SpeciesSet& tspecies       = elec.getSpeciesSet();
@@ -185,7 +186,6 @@ TEST_CASE("SODMC", "[drivers][dmc]")
    <parameter name=\"steps\">1</parameter> \
    <parameter name=\"blocks\">1</parameter> \
    <parameter name=\"timestep\">0.1</parameter> \
-   <parameter name=\"SpinMoves\">yes</parameter> \
    <parameter name=\"SpinMass\">0.25</parameter> \
   </qmc> \
   ";

src/QMCDrivers/tests/test_vmc_driver.cpp

@@ -135,6 +135,7 @@ TEST_CASE("SOVMC", "[drivers][vmc]")
   elec.create(agroup);
   elec.R[0]       = {1.0, 0.0, 0.0};
   elec.spins[0]   = 0.0;
+  elec.is_spinor_ = true;
   elec.createWalkers(1);
 
   SpeciesSet& tspecies       = elec.getSpeciesSet();
@@ -170,7 +171,6 @@ TEST_CASE("SOVMC", "[drivers][vmc]")
    <parameter name=\"blocks\">1</parameter> \
    <parameter name=\"timestep\">0.1</parameter> \
    <parameter name=\"usedrift\">no</parameter> \
-   <parameter name=\"SpinMoves\">yes</parameter> \
    <parameter name=\"SpinMass\">0.25</parameter> \
   </qmc> \
   ";
@@ -221,6 +221,7 @@ TEST_CASE("SOVMC-alle", "[drivers][vmc]")
   elec.create(agroup);
   elec.R[0]       = {1.0, 0.0, 0.0};
   elec.spins[0]   = 0.0;
+  elec.is_spinor_ = true;
   elec.createWalkers(1);
 
   SpeciesSet& tspecies       = elec.getSpeciesSet();
@@ -256,7 +257,6 @@ TEST_CASE("SOVMC-alle", "[drivers][vmc]")
    <parameter name=\"blocks\">1</parameter> \
    <parameter name=\"timestep\">0.1</parameter> \
    <parameter name=\"usedrift\">no</parameter> \
-   <parameter name=\"SpinMoves\">yes</parameter> \
    <parameter name=\"SpinMass\">0.25</parameter> \
   </qmc> \
   ";

src/QMCWaveFunctions/DiffWaveFunctionComponent.h

@@ -27,9 +27,9 @@ namespace qmcplusplus
  * @brief Base base class for derivatives of WaveFunctionComponent
  *
  * Derived classes implement the differentiate function which evaluates
- * - \f$\fraction{\partial \log\Psi}{\partial \alpha}\f$
+ * - \f${\partial \log \Psi}/{\partial \alpha}\f$
- * - \f$\nabla \fraction{\partial \log\Psi}{\partial \alpha}\f$
+ * - \f$\nabla {\partial \log \Psi}/{\partial \alpha}\f$
- * - \f$\nabla^2 \fraction{\partial \log\Psi}{\partial \alpha}\f$
+ * - \f$\nabla^2 {\partial \log \Psi}/{\partial \alpha}\f$
  * Each object handles one or more parameters during the optimiation.
  * The data type of refOrbital, std::vector<WaveFunctionComponent*> is intended for the cases
  * when a common variable is used by several WaveFunctionComponent class, e.g.,

src/QMCWaveFunctions/Fermion/MultiDiracDeterminant.2.cpp

@@ -114,7 +114,7 @@ void MultiDiracDeterminant::BuildDotProductsAndCalculateRatios(int ref,
 {
   const ValueType det0 = ratios(ref, iat)[dx];
   BuildDotProductsAndCalculateRatios_impl(ref, det0, WorkSpace.data(), psiinv, psi, dotProducts, data, pairs, sign);
-  for (size_t count = 0; count < NumDets; ++count)
+  for (size_t count = 0; count < getNumDets(); ++count)
     ratios(count, iat)[dx] = WorkSpace[count];
 #if 0
     ValueType det0 = ratios(ref,iat)[dx];
@@ -160,7 +160,7 @@ void MultiDiracDeterminant::BuildDotProductsAndCalculateRatios(int ref,
   const ValueType det0 = ratios(ref, iat);
   BuildDotProductsAndCalculateRatios_impl(ref, det0, WorkSpace.data(), psiinv, psi, dotProducts, data, pairs, sign);
   //splatt
-  for (size_t count = 0; count < NumDets; ++count)
+  for (size_t count = 0; count < getNumDets(); ++count)
     ratios(count, iat) = WorkSpace[count];
 #if 0
     ValueType det0 = ratios(ref,iat);
@@ -230,15 +230,8 @@ void MultiDiracDeterminant::evaluateDetsForPtclMove(const ParticleSet& P, int ia
     for (size_t i = 0; i < NumOrbitals; i++)
       TpsiM(i, WorkingIndex) = psiV[i];
     ExtraStuffTimer.stop();
-    BuildDotProductsAndCalculateRatios(ReferenceDeterminant,
+    BuildDotProductsAndCalculateRatios(ReferenceDeterminant, WorkingIndex, new_detValues, psiMinv_temp, TpsiM,
-                                       WorkingIndex,
+                                       dotProducts, *detData, *uniquePairs, *DetSigns);
-                                       new_detValues,
-                                       psiMinv_temp,
-                                       TpsiM,
-                                       dotProducts,
-                                       *detData,
-                                       *uniquePairs,
-                                       *DetSigns);
     // check comment above
     for (size_t i = 0; i < NumOrbitals; i++)
       TpsiM(i, WorkingIndex) = psiM(WorkingIndex, i);
@@ -289,15 +282,8 @@ void MultiDiracDeterminant::evaluateDetsAndGradsForPtclMove(const ParticleSet& P
     for (size_t i = 0; i < NumOrbitals; i++)
       TpsiM(i, WorkingIndex) = psiV[i];
     ExtraStuffTimer.stop();
-    BuildDotProductsAndCalculateRatios(ReferenceDeterminant,
+    BuildDotProductsAndCalculateRatios(ReferenceDeterminant, WorkingIndex, new_detValues, psiMinv_temp, TpsiM,
-                                       WorkingIndex,
+                                       dotProducts, *detData, *uniquePairs, *DetSigns);
-                                       new_detValues,
-                                       psiMinv_temp,
-                                       TpsiM,
-                                       dotProducts,
-                                       *detData,
-                                       *uniquePairs,
-                                       *DetSigns);
     for (size_t idim = 0; idim < OHMMS_DIM; idim++)
     {
       ExtraStuffTimer.start();
@@ -310,16 +296,8 @@ void MultiDiracDeterminant::evaluateDetsAndGradsForPtclMove(const ParticleSet& P
       for (size_t i = 0; i < NumOrbitals; i++)
         TpsiM(i, WorkingIndex) = dpsiV[i][idim];
       ExtraStuffTimer.stop();
-      BuildDotProductsAndCalculateRatios(ReferenceDeterminant,
+      BuildDotProductsAndCalculateRatios(ReferenceDeterminant, WorkingIndex, new_grads, dpsiMinv, TpsiM, dotProducts,
-                                         WorkingIndex,
+                                         *detData, *uniquePairs, *DetSigns, idim);
-                                         new_grads,
-                                         dpsiMinv,
-                                         TpsiM,
-                                         dotProducts,
-                                         *detData,
-                                         *uniquePairs,
-                                         *DetSigns,
-                                         idim);
     }
     // check comment above
     for (int i = 0; i < NumOrbitals; i++)
@@ -364,16 +342,8 @@ void MultiDiracDeterminant::evaluateGrads(ParticleSet& P, int iat)
       InverseUpdateByColumn(dpsiMinv, psiV_temp, workV1, workV2, WorkingIndex, ratioG);
       for (size_t i = 0; i < NumOrbitals; i++)
         TpsiM(i, WorkingIndex) = dpsiM(WorkingIndex, i)[idim];
-      BuildDotProductsAndCalculateRatios(ReferenceDeterminant,
+      BuildDotProductsAndCalculateRatios(ReferenceDeterminant, WorkingIndex, grads, dpsiMinv, TpsiM, dotProducts,
-                                         WorkingIndex,
+                                         *detData, *uniquePairs, *DetSigns, idim);
-                                         grads,
-                                         dpsiMinv,
-                                         TpsiM,
-                                         dotProducts,
-                                         *detData,
-                                         *uniquePairs,
-                                         *DetSigns,
-                                         idim);
     }
     // check comment above
     for (size_t i = 0; i < NumOrbitals; i++)
@@ -428,15 +398,8 @@ void MultiDiracDeterminant::evaluateAllForPtclMove(const ParticleSet& P, int iat
     InverseUpdateByColumn(psiMinv_temp, psiV_temp, workV1, workV2, WorkingIndex, ratioRef);
     for (size_t i = 0; i < NumOrbitals; i++)
       TpsiM(i, WorkingIndex) = psiV[i];
-    BuildDotProductsAndCalculateRatios(ReferenceDeterminant,
+    BuildDotProductsAndCalculateRatios(ReferenceDeterminant, WorkingIndex, new_detValues, psiMinv_temp, TpsiM,
-                                       WorkingIndex,
+                                       dotProducts, *detData, *uniquePairs, *DetSigns);
-                                       new_detValues,
-                                       psiMinv_temp,
-                                       TpsiM,
-                                       dotProducts,
-                                       *detData,
-                                       *uniquePairs,
-                                       *DetSigns);
     for (size_t jat = 0; jat < NumPtcls; jat++)
     {
       it = confgList[ReferenceDeterminant].occup.begin();
@@ -461,16 +424,8 @@ void MultiDiracDeterminant::evaluateAllForPtclMove(const ParticleSet& P, int iat
           InverseUpdateByColumn(dpsiMinv, psiV_temp, workV1, workV2, jat, gradRatio[idim]);
           for (size_t i = 0; i < NumOrbitals; i++)
             TpsiM(i, jat) = dpsiV[i][idim];
-          BuildDotProductsAndCalculateRatios(ReferenceDeterminant,
+          BuildDotProductsAndCalculateRatios(ReferenceDeterminant, jat, new_grads, dpsiMinv, TpsiM, dotProducts,
-                                             jat,
+                                             *detData, *uniquePairs, *DetSigns, idim);
-                                             new_grads,
-                                             dpsiMinv,
-                                             TpsiM,
-                                             dotProducts,
-                                             *detData,
-                                             *uniquePairs,
-                                             *DetSigns,
-                                             idim);
         }
         dpsiMinv = psiMinv_temp;
         it       = confgList[ReferenceDeterminant].occup.begin();
@@ -479,15 +434,8 @@ void MultiDiracDeterminant::evaluateAllForPtclMove(const ParticleSet& P, int iat
         InverseUpdateByColumn(dpsiMinv, psiV_temp, workV1, workV2, jat, ratioLapl);
         for (size_t i = 0; i < NumOrbitals; i++)
           TpsiM(i, jat) = d2psiV[i];
-        BuildDotProductsAndCalculateRatios(ReferenceDeterminant,
+        BuildDotProductsAndCalculateRatios(ReferenceDeterminant, jat, new_lapls, dpsiMinv, TpsiM, dotProducts, *detData,
-                                           jat,
+                                           *uniquePairs, *DetSigns);
-                                           new_lapls,
-                                           dpsiMinv,
-                                           TpsiM,
-                                           dotProducts,
-                                           *detData,
-                                           *uniquePairs,
-                                           *DetSigns);
         for (size_t i = 0; i < NumOrbitals; i++)
           TpsiM(i, jat) = psiV[i];
       }
@@ -510,16 +458,8 @@ void MultiDiracDeterminant::evaluateAllForPtclMove(const ParticleSet& P, int iat
           InverseUpdateByColumn(dpsiMinv, psiV_temp, workV1, workV2, jat, gradRatio[idim]);
           for (size_t i = 0; i < NumOrbitals; i++)
             TpsiM(i, jat) = dpsiM(jat, i)[idim];
-          BuildDotProductsAndCalculateRatios(ReferenceDeterminant,
+          BuildDotProductsAndCalculateRatios(ReferenceDeterminant, jat, new_grads, dpsiMinv, TpsiM, dotProducts,
-                                             jat,
+                                             *detData, *uniquePairs, *DetSigns, idim);
-                                             new_grads,
-                                             dpsiMinv,
-                                             TpsiM,
-                                             dotProducts,
-                                             *detData,
-                                             *uniquePairs,
-                                             *DetSigns,
-                                             idim);
         }
         dpsiMinv = psiMinv_temp;
         it       = confgList[ReferenceDeterminant].occup.begin();
@@ -528,15 +468,8 @@ void MultiDiracDeterminant::evaluateAllForPtclMove(const ParticleSet& P, int iat
         InverseUpdateByColumn(dpsiMinv, psiV_temp, workV1, workV2, jat, ratioLapl);
         for (size_t i = 0; i < NumOrbitals; i++)
           TpsiM(i, jat) = d2psiM(jat, i);
-        BuildDotProductsAndCalculateRatios(ReferenceDeterminant,
+        BuildDotProductsAndCalculateRatios(ReferenceDeterminant, jat, new_lapls, dpsiMinv, TpsiM, dotProducts, *detData,
-                                           jat,
+                                           *uniquePairs, *DetSigns);
-                                           new_lapls,
-                                           dpsiMinv,
-                                           TpsiM,
-                                           dotProducts,
-                                           *detData,
-                                           *uniquePairs,
-                                           *DetSigns);
         for (size_t i = 0; i < NumOrbitals; i++)
           TpsiM(i, jat) = psiM(jat, i);
       }

src/QMCWaveFunctions/Fermion/MultiDiracDeterminant.cpp

@@ -30,25 +30,18 @@
 
 namespace qmcplusplus
 {
-/** set the index of the first particle in the determinant and reset the size of the determinant
+void MultiDiracDeterminant::set(int first, int nel, int ref_det_id)
- *@param first index of first particle
- *@param nel number of particles in the determinant
- */
-void MultiDiracDeterminant::set(int first, int nel)
 {
-  APP_ABORT(ClassName + "set(int first, int nel) is disabled. \n");
+  assert(ciConfigList);
-}
+  assert(ciConfigList->size() > 0);
 
-void MultiDiracDeterminant::set(int first, int nel, int norb)
-{
   FirstIndex           = first;
-  DetCalculator.resize(nel);
+  ReferenceDeterminant = ref_det_id;
-  resize(nel, norb);
+  resize(nel);
-  // mmorales; remove later
+  createDetData((*ciConfigList)[ReferenceDeterminant], *detData, *uniquePairs, *DetSigns);
-  //    testDets();
 }
 
-void MultiDiracDeterminant::createDetData(ci_configuration2& ref,
+void MultiDiracDeterminant::createDetData(const ci_configuration2& ref,
                                           std::vector<int>& data,
                                           std::vector<std::pair<int, int>>& pairs,
                                           std::vector<RealType>& sign)
@@ -222,9 +215,6 @@ void MultiDiracDeterminant::copyFromBuffer(ParticleSet& P, WFBufferType& buf)
   buf.get(FirstAddressOfGrads, LastAddressOfGrads);
   buf.get(lapls.first_address(), lapls.last_address());
   // only used with ORB_PBYP_ALL,
-  //psiM_temp = psiM;
-  //dpsiM_temp = dpsiM;
-  //d2psiM_temp = d2psiM;
   psiMinv_temp = psiMinv;
   int n1       = psiM.extent(0);
   int n2       = psiM.extent(1);
@@ -316,24 +306,20 @@ MultiDiracDeterminant::MultiDiracDeterminant(const MultiDiracDeterminant& s)
       readMatGradTimer(*timer_manager.createTimer(ClassName + "::readMatGrad")),
       buildTableGradTimer(*timer_manager.createTimer(ClassName + "::buildTableGrad")),
       ExtraStuffTimer(*timer_manager.createTimer(ClassName + "::ExtraStuff")),
+      Phi(s.Phi->makeClone()),
+      NumOrbitals(Phi->getOrbitalSetSize()),
       NP(0),
       FirstIndex(s.FirstIndex),
-      ciConfigList(nullptr)
+      ciConfigList(s.ciConfigList),
+      ReferenceDeterminant(s.ReferenceDeterminant),
+      detData(s.detData),
+      uniquePairs(s.uniquePairs),
+      DetSigns(s.DetSigns)
 {
-  IsCloned = true;
-
-  ReferenceDeterminant = s.ReferenceDeterminant;
-  ciConfigList         = s.ciConfigList;
-  NumDets              = s.NumDets;
-  detData              = s.detData;
-  uniquePairs          = s.uniquePairs;
-  DetSigns             = s.DetSigns;
   Optimizable = s.Optimizable;
 
   registerTimers();
-  Phi.reset(s.Phi->makeClone());
+  resize(s.NumPtcls);
-  this->resize(s.NumPtcls, s.NumOrbitals);
-  this->DetCalculator.resize(s.NumPtcls);
 }
 
 SPOSetPtr MultiDiracDeterminant::clonePhi() const { return Phi->makeClone(); }
@@ -361,48 +347,24 @@ MultiDiracDeterminant::MultiDiracDeterminant(std::unique_ptr<SPOSet>&& spos, int
       readMatGradTimer(*timer_manager.createTimer(ClassName + "::readMatGrad")),
       buildTableGradTimer(*timer_manager.createTimer(ClassName + "::buildTableGrad")),
       ExtraStuffTimer(*timer_manager.createTimer(ClassName + "::ExtraStuff")),
+      Phi(std::move(spos)),
+      NumOrbitals(Phi->getOrbitalSetSize()),
       NP(0),
       FirstIndex(first),
-      Phi(std::move(spos)),
-      ciConfigList(nullptr),
       ReferenceDeterminant(0)
 {
   (Phi->isOptimizable() == true) ? Optimizable = true : Optimizable = false;
 
-  IsCloned = false;
+  ciConfigList = std::make_shared<std::vector<ci_configuration2>>();
-
+  detData      = std::make_shared<std::vector<int>>();
-  ciConfigList = new std::vector<ci_configuration2>;
+  uniquePairs  = std::make_shared<std::vector<std::pair<int, int>>>();
-  detData      = new std::vector<int>;
+  DetSigns     = std::make_shared<std::vector<RealType>>();
-  uniquePairs  = new std::vector<std::pair<int, int>>;
-  DetSigns     = new std::vector<RealType>;
 
   registerTimers();
 }
 
 ///default destructor
-MultiDiracDeterminant::~MultiDiracDeterminant() {}
+MultiDiracDeterminant::~MultiDiracDeterminant() = default;
-
-MultiDiracDeterminant& MultiDiracDeterminant::operator=(const MultiDiracDeterminant& s)
-{
-  if (this == &s)
-    return *this;
-
-  NP                   = 0;
-  IsCloned             = true;
-  ReferenceDeterminant = s.ReferenceDeterminant;
-  ciConfigList         = s.ciConfigList;
-  NumDets              = s.NumDets;
-
-  detData     = s.detData;
-  uniquePairs = s.uniquePairs;
-  FirstIndex  = s.FirstIndex;
-  DetSigns    = s.DetSigns;
-
-  resize(s.NumPtcls, s.NumOrbitals);
-  this->DetCalculator.resize(s.NumPtcls);
-
-  return *this;
-}
 
 void MultiDiracDeterminant::registerData(ParticleSet& P, WFBufferType& buf)
 {
@@ -412,7 +374,7 @@ void MultiDiracDeterminant::registerData(ParticleSet& P, WFBufferType& buf)
     //int norb = cols();
     NP                  = P.getTotalNum();
     FirstAddressOfGrads = &(grads(0, 0)[0]);
-    LastAddressOfGrads  = FirstAddressOfGrads + NumPtcls * DIM * NumDets;
+    LastAddressOfGrads  = FirstAddressOfGrads + NumPtcls * DIM * getNumDets();
     FirstAddressOfdpsiM = &(dpsiM(0, 0)[0]); //(*dpsiM.begin())[0]);
     LastAddressOfdpsiM  = FirstAddressOfdpsiM + NumPtcls * NumOrbitals * DIM;
   }
@@ -428,40 +390,27 @@ void MultiDiracDeterminant::registerData(ParticleSet& P, WFBufferType& buf)
 }
 
 
-void MultiDiracDeterminant::setDetInfo(int ref, std::vector<ci_configuration2>* list)
-{
-  ReferenceDeterminant = ref;
-  ciConfigList         = list;
-  NumDets              = list->size();
-}
-
 ///reset the size: with the number of particles and number of orbtials
 /// morb is the total number of orbitals, including virtual
-void MultiDiracDeterminant::resize(int nel, int morb)
+void MultiDiracDeterminant::resize(int nel)
 {
-  if (nel <= 0 || morb <= 0)
+  if (nel <= 0)
   {
     APP_ABORT(" ERROR: MultiDiracDeterminant::resize arguments equal to zero. \n");
   }
-  if (NumDets == 0 || NumDets != ciConfigList->size())
+
-  {
+  const int NumDets = getNumDets();
-    APP_ABORT(" ERROR: MultiDiracDeterminant::resize problems with NumDets. \n");
-  }
 
   NumPtcls  = nel;
-  NumOrbitals = morb;
   LastIndex = FirstIndex + nel;
   psiV_temp.resize(nel);
   psiV.resize(NumOrbitals);
   dpsiV.resize(NumOrbitals);
   d2psiV.resize(NumOrbitals);
-  psiM.resize(nel, morb);
+  psiM.resize(nel, NumOrbitals);
-  dpsiM.resize(nel, morb);
+  dpsiM.resize(nel, NumOrbitals);
-  d2psiM.resize(nel, morb);
+  d2psiM.resize(nel, NumOrbitals);
-  //psiM_temp.resize(nel,morb);
+  TpsiM.resize(NumOrbitals, nel);
-  //dpsiM_temp.resize(nel,morb);
-  //d2psiM_temp.resize(nel,morb);
-  TpsiM.resize(morb, nel);
   psiMinv.resize(nel, nel);
   dpsiMinv.resize(nel, nel);
   psiMinv_temp.resize(nel, nel);
@@ -476,14 +425,8 @@ void MultiDiracDeterminant::resize(int nel, int morb)
   new_grads.resize(NumDets, nel);
   lapls.resize(NumDets, nel);
   new_lapls.resize(NumDets, nel);
-  dotProducts.resize(morb, morb);
+  dotProducts.resize(NumOrbitals, NumOrbitals);
-
+  DetCalculator.resize(nel);
-  //if(ciConfigList==nullptr)
-  //{
-  //  APP_ABORT("ciConfigList was not properly initialized.\n");
-  //}
-  if (!IsCloned)
-    createDetData((*ciConfigList)[ReferenceDeterminant], *detData, *uniquePairs, *DetSigns);
 }
 
 void MultiDiracDeterminant::registerTimers()

src/QMCWaveFunctions/Fermion/MultiDiracDeterminant.h

@@ -69,7 +69,7 @@ public:
    */
   MultiDiracDeterminant(const MultiDiracDeterminant& s);
 
-  MultiDiracDeterminant& operator=(const MultiDiracDeterminant& s);
+  MultiDiracDeterminant& operator=(const MultiDiracDeterminant& s) = delete;
 
   /** return a clone of Phi
    */
@@ -77,18 +77,14 @@ public:
 
   SPOSetPtr getPhi() { return Phi.get(); };
 
-  inline IndexType rows() const { return NumPtcls; }
-
-  inline IndexType cols() const { return NumOrbitals; }
-
   /** set the index of the first particle in the determinant and reset the size of the determinant
    *@param first index of first particle
    *@param nel number of particles in the determinant
-   *@param norb total number of orbitals (including unoccupied)
+   *@param ref_det_id id of the reference determinant
+   *
+   * Note: ciConfigList should have been populated when calling this function
    */
-  void set(int first, int nel, int norb);
+  void set(int first, int nel, int ref_det_id);
-
-  void set(int first, int nel);
 
   void setBF(BackflowTransformation* bf) {}
 
@@ -159,29 +155,13 @@ public:
     const ValueMatrix_t& Minv_up      = psiMinv;
     const ValueMatrix_t& Minv_dn      = pseudo_dn.psiMinv;
 
-    Phi->evaluateDerivativesWF(P,
+    Phi->evaluateDerivativesWF(P, optvars, dlogpsi, psiCurrent, Coeff, C2node_up, C2node_dn, detValues_up, detValues_dn,
-                                optvars, 
+                               M_up, M_dn, Minv_up, Minv_dn, *detData, lookup_tbl);
-                                dlogpsi, 
-                                psiCurrent, 
-                                Coeff, 
-                                C2node_up, 
-                                C2node_dn, 
-                                detValues_up,
-                                detValues_dn, 
-                                M_up, 
-                                M_dn, 
-                                Minv_up, 
-                                Minv_dn, 
-                                *detData, 
-                                lookup_tbl);
   }
 
 
   inline void reportStatus(std::ostream& os) {}
 
-  ///reset the size: with the number of particles and number of orbtials
-  virtual void resize(int nel, int morb);
-
   void registerData(ParticleSet& P, WFBufferType& buf);
 
   LogValueType updateBuffer(ParticleSet& P, WFBufferType& buf, bool fromscratch = false);
@@ -220,7 +200,9 @@ public:
     return PsiValueType();
   }
 
-  LogValueType evaluateLog(const ParticleSet& P, ParticleSet::ParticleGradient_t& G, ParticleSet::ParticleLaplacian_t& L)
+  LogValueType evaluateLog(const ParticleSet& P,
+                           ParticleSet::ParticleGradient_t& G,
+                           ParticleSet::ParticleLaplacian_t& L)
   {
     APP_ABORT("  MultiDiracDeterminant: This should not be called. \n");
     return 0.0;
@@ -239,7 +221,7 @@ public:
 
   // create necessary structures used in the evaluation of the determinants
   // this works with confgList, which shouldn't change during a simulation
-  void createDetData(ci_configuration2& ref,
+  void createDetData(const ci_configuration2& ref,
                      std::vector<int>& data,
                      std::vector<std::pair<int, int>>& pairs,
                      std::vector<RealType>& sign);
@@ -381,8 +363,6 @@ public:
       }
   */
 
-  void setDetInfo(int ref, std::vector<ci_configuration2>* list);
-
   /** evaluate the value of all the unique determinants with one electron moved. Used by the table method
    *@param P particle set which provides the positions
    *@param iat the index of the moved electron
@@ -400,43 +380,49 @@ public:
   // full evaluation of all the structures from scratch, used in evaluateLog for example
   void evaluateForWalkerMove(const ParticleSet& P, bool fromScratch = true);
 
+  // accessors
+  inline int getNumDets() const { return ciConfigList->size(); }
+  inline int getNumPtcls() const { return NumPtcls; }
+  inline int getFirstIndex() const { return FirstIndex; }
+  inline std::vector<ci_configuration2>& getCIConfigList() { return *ciConfigList; }
+
+  /// store determinants (old and new). FIXME: move to private
+  ValueVector_t detValues, new_detValues;
+  GradMatrix_t grads, new_grads;
+  ValueMatrix_t lapls, new_lapls;
+
+private:
+  ///reset the size: with the number of particles
+  void resize(int nel);
+
+  ///a set of single-particle orbitals used to fill in the  values of the matrix
+  const std::unique_ptr<SPOSet> Phi;
+  ///number of single-particle orbitals which belong to this Dirac determinant
+  const int NumOrbitals;
   ///total number of particles
   int NP;
-  ///number of single-particle orbitals which belong to this Dirac determinant
-  int NumOrbitals;
   ///number of particles which belong to this Dirac determinant
   int NumPtcls;
   ///index of the first particle with respect to the particle set
   int FirstIndex;
   ///index of the last particle with respect to the particle set
   int LastIndex;
-  ///a set of single-particle orbitals used to fill in the  values of the matrix
-  std::unique_ptr<SPOSet> Phi;
-  /// number of determinants handled by this object
-  int NumDets;
-  ///bool to cleanup
-  bool IsCloned;
   ///use shared_ptr
-  std::vector<ci_configuration2>* ciConfigList;
+  std::shared_ptr<std::vector<ci_configuration2>> ciConfigList;
   // the reference determinant never changes, so there is no need to store it.
   // if its value is zero, then use a data from backup, but always use this one
   // by default
   int ReferenceDeterminant;
 
-  /// store determinants (old and new)
-  ValueVector_t detValues, new_detValues;
-  GradMatrix_t grads, new_grads;
-  ValueMatrix_t lapls, new_lapls;
-
   /// psiM(i,j) \f$= \psi_j({\bf r}_i)\f$
   /// TpsiM(i,j) \f$= psiM(j,i) \f$
-  ValueMatrix_t psiM, psiM_temp, TpsiM, psiMinv, psiMinv_temp;
+  ValueMatrix_t psiM, TpsiM, psiMinv, psiMinv_temp;
   /// dpsiM(i,j) \f$= \nabla_i \psi_j({\bf r}_i)\f$
-  GradMatrix_t dpsiM, dpsiM_temp;
+  GradMatrix_t dpsiM;
   // temporaty storage
   ValueMatrix_t dpsiMinv;
   /// d2psiM(i,j) \f$= \nabla_i^2 \psi_j({\bf r}_i)\f$
-  ValueMatrix_t d2psiM, d2psiM_temp;
+  ValueMatrix_t d2psiM;
 
   /// value of single-particle orbital for particle-by-particle update
   ValueVector_t psiV, psiV_temp;
@@ -463,9 +449,9 @@ public:
    *     -i1,i2,...,in : occupied orbital to be replaced (these must be numbers from 0:Nptcl-1)
    *     -a1,a2,...,an : excited states that replace the orbitals (these can be anything)
    */
-  std::vector<int>* detData;
+  std::shared_ptr<std::vector<int>> detData;
-  std::vector<std::pair<int, int>>* uniquePairs;
+  std::shared_ptr<std::vector<std::pair<int, int>>> uniquePairs;
-  std::vector<RealType>* DetSigns;
+  std::shared_ptr<std::vector<RealType>> DetSigns;
   MultiDiracDeterminantCalculator<ValueType> DetCalculator;
 };
 

src/QMCWaveFunctions/Fermion/MultiSlaterDeterminant.cpp

@@ -56,8 +56,10 @@ MultiSlaterDeterminant::MultiSlaterDeterminant(ParticleSet& targetPtcl,
 WaveFunctionComponentPtr MultiSlaterDeterminant::makeClone(ParticleSet& tqp) const
 {
   typedef DiracDeterminant<> SingleDet_t;
-  auto spo_up_C = std::make_unique<SPOSetProxyForMSD>(std::unique_ptr<SPOSet>(spo_up->refPhi->makeClone()), FirstIndex_up, LastIndex_up);
+  auto spo_up_C   = std::make_unique<SPOSetProxyForMSD>(std::unique_ptr<SPOSet>(spo_up->refPhi->makeClone()),
-  auto spo_dn_C = std::make_unique<SPOSetProxyForMSD>(std::unique_ptr<SPOSet>(spo_dn->refPhi->makeClone()), FirstIndex_dn, LastIndex_dn);
+                                                      FirstIndex_up, LastIndex_up);
+  auto spo_dn_C   = std::make_unique<SPOSetProxyForMSD>(std::unique_ptr<SPOSet>(spo_dn->refPhi->makeClone()),
+                                                      FirstIndex_dn, LastIndex_dn);
   spo_up_C->occup = spo_up->occup;
   spo_dn_C->occup = spo_dn->occup;
   MultiSlaterDeterminant* clone = new MultiSlaterDeterminant(tqp, std::move(spo_up_C), std::move(spo_dn_C));
@@ -70,35 +72,15 @@ WaveFunctionComponentPtr MultiSlaterDeterminant::makeClone(ParticleSet& tqp) con
     clone->CSFexpansion = CSFexpansion;
     clone->DetsPerCSF   = DetsPerCSF;
   }
-  //     SPOSetProxyForMSD* spo = clone->spo_up;
-  //     spo->occup.resize(uniqueConfg_up.size(),clone->nels_up);
   for (int i = 0; i < dets_up.size(); i++)
   {
-    //       int nq=0;
+    clone->dets_up.push_back(std::make_unique<SingleDet_t>(std::static_pointer_cast<SPOSet>(clone->spo_up), 0));
-    // //       configuration& ci = uniqueConfg_up[i];
+    clone->dets_up.back()->set(clone->FirstIndex_up, clone->nels_up);
-    //       for(int k=0; k<uniqueConfg_up[i].occup.size(); k++) {
-    //         if(uniqueConfg_up[i].occup[k]) {
-    //           spo->occup(i,nq++) = k;
-    //         }
-    //       }
-    SingleDet_t* adet = new SingleDet_t(std::static_pointer_cast<SPOSet>(clone->spo_up), 0);
-    adet->set(clone->FirstIndex_up, clone->nels_up);
-    clone->dets_up.push_back(adet);
   }
-  //     spo = clone->spo_dn;
-  //     spo->occup.resize(uniqueConfg_dn.size(),clone->nels_dn);
   for (int i = 0; i < dets_dn.size(); i++)
   {
-    //       int nq=0;
+    clone->dets_dn.emplace_back(std::make_unique<SingleDet_t>(std::static_pointer_cast<SPOSet>(clone->spo_dn), 0));
-    // //       configuration& ci = uniqueConfg_dn[i];
+    clone->dets_dn.back()->set(clone->FirstIndex_dn, clone->nels_dn);
-    //       for(int k=0; k<uniqueConfg_dn[i].occup.size(); k++) {
-    //         if(uniqueConfg_dn[i].occup[k]) {
-    //           spo->occup(i,nq++) = k;
-    //         }
-    //       }
-    SingleDet_t* adet = new SingleDet_t(std::static_pointer_cast<SPOSet>(clone->spo_dn), 0);
-    adet->set(clone->FirstIndex_dn, clone->nels_dn);
-    clone->dets_dn.push_back(adet);
   }
   clone->Optimizable = Optimizable;
   clone->C           = C;

src/QMCWaveFunctions/Fermion/MultiSlaterDeterminant.h

@@ -54,8 +54,6 @@ public:
   NewTimer &RatioTimer, &RatioGradTimer, &RatioAllTimer, &UpdateTimer, &EvaluateTimer;
   NewTimer &Ratio1Timer, &Ratio1GradTimer, &Ratio1AllTimer, &AccRejTimer, &evalOrbTimer;
 
-  typedef DiracDeterminantBase* DiracDeterminantBasePtr;
-  typedef SPOSet* SPOSetPtr;
   typedef OrbitalSetTraits<ValueType>::IndexVector_t IndexVector_t;
   typedef OrbitalSetTraits<ValueType>::ValueVector_t ValueVector_t;
   typedef OrbitalSetTraits<ValueType>::GradVector_t GradVector_t;
@@ -85,7 +83,9 @@ public:
   ///set BF pointers
   virtual void setBF(BackflowTransformation* BFTrans) {}
 
-  virtual ValueType evaluate(const ParticleSet& P, ParticleSet::ParticleGradient_t& G, ParticleSet::ParticleLaplacian_t& L);
+  virtual ValueType evaluate(const ParticleSet& P,
+                             ParticleSet::ParticleGradient_t& G,
+                             ParticleSet::ParticleLaplacian_t& L);
 
   virtual LogValueType evaluateLog(const ParticleSet& P,
                                    ParticleSet::ParticleGradient_t& G,
@@ -127,8 +127,8 @@ public:
   std::shared_ptr<SPOSetProxyForMSD> spo_up;
   std::shared_ptr<SPOSetProxyForMSD> spo_dn;
 
-  std::vector<DiracDeterminantBasePtr> dets_up;
+  std::vector<std::unique_ptr<DiracDeterminantBase>> dets_up;
-  std::vector<DiracDeterminantBasePtr> dets_dn;
+  std::vector<std::unique_ptr<DiracDeterminantBase>> dets_dn;
 
   // map determinant in linear combination to unique det list
   std::vector<size_t> C2node_up;

src/QMCWaveFunctions/Fermion/MultiSlaterDeterminantFast.cpp

@@ -174,12 +174,12 @@ WaveFunctionComponent::PsiValueType MultiSlaterDeterminantFast::evaluate_vgl_imp
   for (size_t ig = 0; ig < Dets.size(); ig++)
     precomputeC_otherDs(P, ig);
 
-  for (size_t i = 0; i < Dets[0]->NumDets; ++i)
+  for (size_t i = 0; i < Dets[0]->getNumDets(); ++i)
     psi += C_otherDs[0][i] * Dets[0]->detValues[i];
 
   for (size_t id = 0; id < Dets.size(); id++)
-    for (size_t i = 0; i < Dets[id]->NumDets; ++i)
+    for (size_t i = 0; i < Dets[id]->getNumDets(); ++i)
-      for (int k = 0, n = Dets[id]->FirstIndex; k < Dets[id]->NumPtcls; k++, n++)
+      for (int k = 0, n = Dets[id]->getFirstIndex(); k < Dets[id]->getNumPtcls(); k++, n++)
       {
         g_tmp[n] += C_otherDs[id][i] * Dets[id]->grads(i, k);
         l_tmp[n] += C_otherDs[id][i] * Dets[id]->lapls(i, k);
@@ -229,10 +229,10 @@ WaveFunctionComponent::PsiValueType MultiSlaterDeterminantFast::evalGrad_impl(Pa
 
   const GradMatrix_t& grads            = (newpos) ? Dets[det_id]->new_grads : Dets[det_id]->grads;
   const ValueType* restrict detValues0 = (newpos) ? Dets[det_id]->new_detValues.data() : Dets[det_id]->detValues.data();
-  const size_t noffset                 = Dets[det_id]->FirstIndex;
+  const size_t noffset                 = Dets[det_id]->getFirstIndex();
 
   PsiValueType psi(0);
-  for (size_t i = 0; i < Dets[det_id]->NumDets; i++)
+  for (size_t i = 0; i < Dets[det_id]->getNumDets(); i++)
   {
     psi += detValues0[i] * C_otherDs[det_id][i];
     g_at += C_otherDs[det_id][i] * grads(i, iat - noffset);
@@ -257,7 +257,7 @@ WaveFunctionComponent::PsiValueType MultiSlaterDeterminantFast::evalGrad_impl_no
   const size_t* restrict det0          = (*C2node)[det_id].data();
   const ValueType* restrict cptr       = C->data();
   const size_t nc                      = C->size();
-  const size_t noffset                 = Dets[det_id]->FirstIndex;
+  const size_t noffset                 = Dets[det_id]->getFirstIndex();
   PsiValueType psi(0);
   for (size_t i = 0; i < nc; ++i)
   {
@@ -324,7 +324,7 @@ WaveFunctionComponent::PsiValueType MultiSlaterDeterminantFast::ratio_impl(Parti
   // psi=Det_Coeff[i]*Det_Value[unique_det_up]*Det_Value[unique_det_dn]*Det_Value[unique_det_AnyOtherType]
   // Since only one electron group is moved at the time, identified by det_id, We precompute:
   // C_otherDs[det_id][i]=Det_Coeff[i]*Det_Value[unique_det_dn]*Det_Value[unique_det_AnyOtherType]
-  for (size_t i = 0; i < Dets[det_id]->NumDets; i++)
+  for (size_t i = 0; i < Dets[det_id]->getNumDets(); i++)
     psi += detValues0[i] * C_otherDs[det_id][i];
 
   return psi;
@@ -385,7 +385,7 @@ void MultiSlaterDeterminantFast::evaluateRatios(const VirtualParticleSet& VP, st
 
     PsiValueType psiNew(0);
     if (use_pre_computing_)
-      for (size_t i = 0; i < Dets[det_id]->NumDets; i++)
+      for (size_t i = 0; i < Dets[det_id]->getNumDets(); i++)
         psiNew += detValues0[i] * C_otherDs[det_id][i];
     else
     {
@@ -591,7 +591,7 @@ void MultiSlaterDeterminantFast::evaluateDerivatives(ParticleSet& P,
         for (size_t i = 0; i < laplSum[id].size(); i++)
         {
           laplSum[id][i] = 0.0;
-          for (size_t k = 0; k < Dets[id]->NumPtcls; k++)
+          for (size_t k = 0; k < Dets[id]->getNumPtcls(); k++)
             laplSum[id][i] += Dets[id]->lapls[i][k];
         }
       }
@@ -599,12 +599,12 @@ void MultiSlaterDeterminantFast::evaluateDerivatives(ParticleSet& P,
       ValueType lapl_sum = 0.0;
       myG_temp           = 0.0;
       for (size_t id = 0; id < Dets.size(); id++)
-        for (size_t i = 0; i < Dets[id]->NumDets; i++)
+        for (size_t i = 0; i < Dets[id]->getNumDets(); i++)
         {
           // assume C_otherDs prepared by evaluateLog already
           ValueType tmp = C_otherDs[id][i] * psiinv;
           lapl_sum += tmp * laplSum[id][i];
-          for (size_t k = 0, j = Dets[id]->FirstIndex; k < Dets[id]->NumPtcls; k++, j++)
+          for (size_t k = 0, j = Dets[id]->getFirstIndex(); k < Dets[id]->getNumPtcls(); k++, j++)
             myG_temp[j] += tmp * Dets[id]->grads(i, k);
         }
 
@@ -646,7 +646,7 @@ void MultiSlaterDeterminantFast::evaluateDerivatives(ParticleSet& P,
                 tmp *= detValues_otherspin[otherspinC];
               }
               q0 += tmp * laplSum[id][spinC];
-              for (size_t l = 0, j = Dets[id]->FirstIndex; l < Dets[id]->NumPtcls; l++, j++)
+              for (size_t l = 0, j = Dets[id]->getFirstIndex(); l < Dets[id]->getNumPtcls(); l++, j++)
                 v[id] += tmp *
                     static_cast<ValueType>(dot(P.G[j], grads_spin(spinC, l)) - dot(myG_temp[j], grads_spin(spinC, l)));
             }
@@ -681,7 +681,7 @@ void MultiSlaterDeterminantFast::evaluateDerivatives(ParticleSet& P,
               tmp *= Dets[other_id]->detValues[otherspinC];
             }
             q0 += tmp * laplSum[id][spinC];
-            for (size_t l = 0, j = Dets[id]->FirstIndex; l < Dets[id]->NumPtcls; l++, j++)
+            for (size_t l = 0, j = Dets[id]->getFirstIndex(); l < Dets[id]->getNumPtcls(); l++, j++)
               v[id] += tmp *
                   static_cast<ValueType>(dot(P.G[j], grads_spin(spinC, l)) - dot(myG_temp[j], grads_spin(spinC, l)));
           }
@@ -833,7 +833,7 @@ void MultiSlaterDeterminantFast::precomputeC_otherDs(const ParticleSet& P, int i
   // C_otherDs(2, :) stores C x D_up x D_dn
 
   ScopedTimer local_timer(PrepareGroupTimer);
-  C_otherDs[ig].resize(Dets[ig]->NumDets);
+  C_otherDs[ig].resize(Dets[ig]->getNumDets());
   std::fill(C_otherDs[ig].begin(), C_otherDs[ig].end(), ValueType(0));
   for (size_t i = 0; i < C->size(); i++)
   {

src/QMCWaveFunctions/Fermion/MultiSlaterDeterminantWithBackflow.cpp

@@ -22,7 +22,8 @@ MultiSlaterDeterminantWithBackflow::MultiSlaterDeterminantWithBackflow(ParticleS
                                                                        std::unique_ptr<SPOSetProxyForMSD>&& upspo,
                                                                        std::unique_ptr<SPOSetProxyForMSD>&& dnspo,
                                                                        BackflowTransformation* BF)
-    : MultiSlaterDeterminant(targetPtcl, std::move(upspo), std::move(dnspo), "MultiSlaterDeterminantWithBackflow"), BFTrans(BF)
+    : MultiSlaterDeterminant(targetPtcl, std::move(upspo), std::move(dnspo), "MultiSlaterDeterminantWithBackflow"),
+      BFTrans(BF)
 {
   Optimizable  = false;
   is_fermionic = true;
@@ -32,11 +33,14 @@ WaveFunctionComponentPtr MultiSlaterDeterminantWithBackflow::makeClone(ParticleS
 {
   // mmorales: the proxy classes read from the particle set inside BFTrans
   BackflowTransformation* tr = BFTrans->makeClone(tqp);
-  auto spo_up_C = std::make_unique<SPOSetProxyForMSD>(std::unique_ptr<SPOSet>(spo_up->refPhi->makeClone()), FirstIndex_up, LastIndex_up);
+  auto spo_up_C              = std::make_unique<SPOSetProxyForMSD>(std::unique_ptr<SPOSet>(spo_up->refPhi->makeClone()),
-  auto spo_dn_C = std::make_unique<SPOSetProxyForMSD>(std::unique_ptr<SPOSet>(spo_dn->refPhi->makeClone()), FirstIndex_dn, LastIndex_dn);
+                                                      FirstIndex_up, LastIndex_up);
+  auto spo_dn_C              = std::make_unique<SPOSetProxyForMSD>(std::unique_ptr<SPOSet>(spo_dn->refPhi->makeClone()),
+                                                      FirstIndex_dn, LastIndex_dn);
   spo_up_C->occup            = spo_up->occup;
   spo_dn_C->occup            = spo_dn->occup;
-  MultiSlaterDeterminantWithBackflow* clone = new MultiSlaterDeterminantWithBackflow(tqp, std::move(spo_up_C), std::move(spo_dn_C), tr);
+  MultiSlaterDeterminantWithBackflow* clone =
+      new MultiSlaterDeterminantWithBackflow(tqp, std::move(spo_up_C), std::move(spo_dn_C), tr);
   clone->C2node_up = C2node_up;
   clone->C2node_dn = C2node_dn;
   clone->resize(dets_up.size(), dets_dn.size());
@@ -48,15 +52,15 @@ WaveFunctionComponentPtr MultiSlaterDeterminantWithBackflow::makeClone(ParticleS
   }
   for (int i = 0; i < dets_up.size(); i++)
   {
-    DiracDeterminantWithBackflow* dclne = (DiracDeterminantWithBackflow*)dets_up[i]->makeCopy(std::static_pointer_cast<SPOSet>(clone->spo_up));
+    clone->dets_up.emplace_back(dets_up[i]->makeCopy(std::static_pointer_cast<SPOSet>(clone->spo_up)));
-    dclne->BFTrans                      = tr;
+    auto& dclne = dynamic_cast<DiracDeterminantWithBackflow&>(*clone->dets_up.back());
-    clone->dets_up.push_back(dclne);
+    dclne.BFTrans = tr;
   }
   for (int i = 0; i < dets_dn.size(); i++)
   {
-    DiracDeterminantWithBackflow* dclne = (DiracDeterminantWithBackflow*)dets_dn[i]->makeCopy(std::static_pointer_cast<SPOSet>(clone->spo_dn));
+    clone->dets_dn.emplace_back(dets_dn[i]->makeCopy(std::static_pointer_cast<SPOSet>(clone->spo_dn)));
-    dclne->BFTrans                      = tr;
+    auto& dclne = dynamic_cast<DiracDeterminantWithBackflow&>(*clone->dets_dn.back());
-    clone->dets_dn.push_back(dclne);
+    dclne.BFTrans = tr;
   }
   clone->Optimizable = Optimizable;
   clone->C           = C;

src/QMCWaveFunctions/Fermion/SlaterDetBuilder.cpp

@@ -594,10 +594,7 @@ bool SlaterDetBuilder::createMSDFast(std::vector<std::unique_ptr<MultiDiracDeter
 
   for (int grp = 0; grp < nGroups; grp++)
   {
-    // you should choose the det with highest weight for reference
+    std::vector<ci_configuration2>& list = Dets[grp]->getCIConfigList();
-    Dets[grp]->ReferenceDeterminant      = 0; // for now
-    Dets[grp]->NumDets                   = uniqueConfgs[grp].size();
-    std::vector<ci_configuration2>& list = *Dets[grp]->ciConfigList;
     list.resize(uniqueConfgs[grp].size());
     for (int i = 0; i < list.size(); i++)
     {
@@ -612,7 +609,8 @@ bool SlaterDetBuilder::createMSDFast(std::vector<std::unique_ptr<MultiDiracDeter
                   << grp << ", problems with ci configuration list. \n");
       }
     }
-    Dets[grp]->set(targetPtcl.first(grp), nptcls[grp], Dets[grp]->Phi->getOrbitalSetSize());
+    // you should choose the det with highest weight for reference. for now choosing 0
+    Dets[grp]->set(targetPtcl.first(grp), nptcls[grp], 0);
   }
 
   if (CSFcoeff.size() == 1)
@@ -725,9 +723,11 @@ bool SlaterDetBuilder::createMSD(MultiSlaterDeterminant* multiSD, xmlNodePtr cur
   multiSD->C2node_up = C2nodes[0];
   multiSD->C2node_dn = C2nodes[1];
   multiSD->resize(uniqueConfgs[0].size(), uniqueConfgs[1].size());
+  // alpha dets
   {
     auto& spo = multiSD->spo_up;
     spo->occup.resize(uniqueConfgs[0].size(), multiSD->nels_up);
+    multiSD->dets_up.reserve(uniqueConfgs[0].size());
     for (int i = 0; i < uniqueConfgs[0].size(); i++)
     {
       int nq               = 0;
@@ -749,12 +749,14 @@ bool SlaterDetBuilder::createMSD(MultiSlaterDeterminant* multiSD, xmlNodePtr cur
         adet = new DiracDeterminant<>(std::static_pointer_cast<SPOSet>(spo), 0);
       }
       adet->set(multiSD->FirstIndex_up, multiSD->nels_up);
-      multiSD->dets_up.push_back(adet);
+      multiSD->dets_up.emplace_back(adet);
     }
   }
+  // beta dets
   {
     auto& spo = multiSD->spo_dn;
     spo->occup.resize(uniqueConfgs[1].size(), multiSD->nels_dn);
+    multiSD->dets_dn.reserve(uniqueConfgs[1].size());
     for (int i = 0; i < uniqueConfgs[1].size(); i++)
     {
       int nq               = 0;
@@ -776,7 +778,7 @@ bool SlaterDetBuilder::createMSD(MultiSlaterDeterminant* multiSD, xmlNodePtr cur
         adet = new DiracDeterminant<>(std::static_pointer_cast<SPOSet>(spo), 0);
       }
       adet->set(multiSD->FirstIndex_dn, multiSD->nels_dn);
-      multiSD->dets_dn.push_back(adet);
+      multiSD->dets_dn.emplace_back(adet);
     }
   }
   if (multiSD->CSFcoeff.size() == 1 || multiSD->C.size() == 1)

src/QMCWaveFunctions/Fermion/ci_configuration2.h

@@ -18,6 +18,7 @@
 #include "CPU/SIMD/simd.hpp"
 #include <algorithm>
 #include <iostream>
+#include <stdexcept>
 
 namespace qmcplusplus
 {
@@ -37,16 +38,11 @@ struct ci_configuration2
   inline bool operator==(const ci_configuration2& c) const
   {
     if (occup.size() != c.occup.size())
-    {
+      throw std::runtime_error("ci_configuration2::operator==() - ci_configuration2s are not compatible.");
-      APP_ABORT("ci_configuration2::operator==() - ci_configuration2s are not compatible.");
+
-    }
     for (int i = 0; i < occup.size(); i++)
-    {
       if (occup[i] != c.occup[i])
-      {
         return false;
-      }
-    }
     return true;
   }
 
@@ -57,15 +53,12 @@ struct ci_configuration2
   size_t calculateNumOfExcitations(const ci_configuration2& c) const
   {
     if (occup.size() != c.occup.size())
-    {
+      throw std::runtime_error("ci_configuration2::operator==() - ci_configuration2s are not compatible.");
-      APP_ABORT("ci_configuration2::operator==() - ci_configuration2s are not compatible.");
+
-    }
     size_t n = 0;
     for (size_t i = 0; i < occup.size(); i++)
-    {
       if (std::find(c.occup.begin(), c.occup.end(), occup[i]) == c.occup.end())
         n++;
-    }
     return n;
   }
 
@@ -86,30 +79,27 @@ struct ci_configuration2
                               std::vector<size_t>& uno) const
   {
     if (occup.size() != c.occup.size())
-    {
+      throw std::runtime_error("ci_configuration2::operator==() - ci_configuration2s are not compatible.");
-      APP_ABORT("ci_configuration2::operator==() - ci_configuration2s are not compatible.");
+
-    }
     n = 0;
     for (size_t i = 0; i < occup.size(); i++)
-    {
       if (std::find(c.occup.begin(), c.occup.end(), occup[i]) == c.occup.end())
       {
         pos[n]   = i;
         ocp[n++] = occup[i];
       }
-    }
+
     if (n == 0)
       return 1.0;
+
     size_t cnt = 0;
     for (size_t i = 0; i < c.occup.size(); i++)
-    {
       if (std::find(occup.begin(), occup.end(), c.occup[i]) == occup.end())
         uno[cnt++] = c.occup[i];
-    }
+
     if (cnt != n)
-    {
+      throw std::runtime_error(" Error #1 in ci_configuration2::calculateExcitations() \n");
-      APP_ABORT(" Error #1 in ci_configuration2::calculateExcitations() \n");
+
-    }
     double res = 1.0;
     // this is needed because ci coefficients are given wrt standard ordering,
     // but by defining the determinant through excitations from a reference might change
@@ -118,9 +108,9 @@ struct ci_configuration2
     auto ref0(occup);
     for (size_t i = 0; i < n; i++)
       ref0[pos[i]] = uno[i];
+
     for (size_t i = 0; i < ref0.size(); i++)
       for (size_t k = i + 1; k < ref0.size(); k++)
-      {
         if (ref0[i] > ref0[k])
         {
           size_t q = ref0[i];
@@ -129,10 +119,8 @@ struct ci_configuration2
           res *= -1.0;
         }
         else if (ref0[i] == ref0[k])
-        {
+          throw std::runtime_error(" Error #2 in ci_configuration2::calculateExcitations() \n");
-          APP_ABORT(" Error #2 in ci_configuration2::calculateExcitations() \n");
+
-        }
-      }
     return res;
   }
 };

src/QMCWaveFunctions/LCAO/LCAOSpinorBuilder.cpp

@@ -47,8 +47,10 @@ SPOSet* LCAOSpinorBuilder::createSPOSetFromXML(xmlNodePtr cur)
   if (optimize == "yes")
     app_log() << "  SPOSet " << spo_name << " is optimizable\n";
 
-  std::unique_ptr<LCAOrbitalSet> upspo = std::make_unique<LCAOrbitalSet>(std::unique_ptr<BasisSet_t>(myBasisSet->makeClone()), optimize == "yes");
+  std::unique_ptr<LCAOrbitalSet> upspo =
-  std::unique_ptr<LCAOrbitalSet> dnspo = std::make_unique<LCAOrbitalSet>(std::unique_ptr<BasisSet_t>(myBasisSet->makeClone()), optimize == "yes");
+      std::make_unique<LCAOrbitalSet>(std::unique_ptr<BasisSet_t>(myBasisSet->makeClone()), optimize == "yes");
+  std::unique_ptr<LCAOrbitalSet> dnspo =
+      std::make_unique<LCAOrbitalSet>(std::unique_ptr<BasisSet_t>(myBasisSet->makeClone()), optimize == "yes");
 
   loadMO(*upspo, *dnspo, cur);
 

src/QMCWaveFunctions/LCAO/LCAOrbitalSet.h

@@ -271,7 +271,6 @@ private:
                                   GradMatrix_t& dlogdet,
                                   HessMatrix_t& dglogdet,
                                   GradMatrix_t& dllogdet) const;
-
 };
 } // namespace qmcplusplus
 #endif

src/QMCWaveFunctions/LCAO/LCAOrbitalSetWithCorrection.cpp

@@ -29,10 +29,7 @@ void LCAOrbitalSetWithCorrection::setOrbitalSetSize(int norbs)
 }
 
 
-SPOSet* LCAOrbitalSetWithCorrection::makeClone() const
+SPOSet* LCAOrbitalSetWithCorrection::makeClone() const { return new LCAOrbitalSetWithCorrection(*this); }
-{
-  return new LCAOrbitalSetWithCorrection(*this);
-}
 
 void LCAOrbitalSetWithCorrection::evaluateValue(const ParticleSet& P, int iat, ValueVector_t& psi)
 {

src/QMCWaveFunctions/SPOSetBuilderFactory.cpp

@@ -41,7 +41,6 @@
 
 namespace qmcplusplus
 {
-
 SPOSet* SPOSetBuilderFactory::getSPOSet(const std::string& name) const
 {
   int nfound  = 0;
@@ -164,7 +163,9 @@ SPOSetBuilder* SPOSetBuilderFactory::createSPOSetBuilder(xmlNodePtr rootNode)
     bb = new SHOSetBuilder(targetPtcl, myComm);
   }
 #if OHMMS_DIM == 3
-  else if (type == "spinorbspline")
+  else if (type.find("spline") < type.size())
+  {
+    if (targetPtcl.is_spinor_)
     {
 #ifdef QMC_COMPLEX
       app_log() << "Einspline Spinor Set\n";
@@ -173,7 +174,7 @@ SPOSetBuilder* SPOSetBuilderFactory::createSPOSetBuilder(xmlNodePtr rootNode)
       PRE.error("Use of einspline spinors requires QMC_COMPLEX=1.  Rebuild with this option");
 #endif
     }
-  else if (type.find("spline") < type.size())
+    else
     {
 #if defined(HAVE_EINSPLINE)
       PRE << "EinsplineSetBuilder:  using libeinspline for B-spline orbitals.\n";
@@ -182,6 +183,7 @@ SPOSetBuilder* SPOSetBuilderFactory::createSPOSetBuilder(xmlNodePtr rootNode)
       PRE.error("Einspline is missing for B-spline orbitals", true);
 #endif
     }
+  }
   else if (type == "molecularorbital" || type == "mo")
   {
     ParticleSet* ions = 0;
@@ -191,22 +193,14 @@ SPOSetBuilder* SPOSetBuilderFactory::createSPOSetBuilder(xmlNodePtr rootNode)
       PRE.error("Missing basisset/@source.", true);
     else
       ions = (*pit).second;
-    bb = new LCAOrbitalBuilder(targetPtcl, *ions, myComm, rootNode);
+    if (targetPtcl.is_spinor_)
-  }
-  else if (type == "molecularspinor")
-  {
 #ifdef QMC_COMPLEX
-    ParticleSet* ions = 0;
-    //initialize with the source tag
-    PtclPoolType::iterator pit(ptclPool.find(sourceOpt));
-    if (pit == ptclPool.end())
-      PRE.error("Missing basisset/@source.", true);
-    else
-      ions = (*pit).second;
       bb = new LCAOSpinorBuilder(targetPtcl, *ions, myComm, rootNode);
 #else
       PRE.error("Use of lcao spinors requires QMC_COMPLEX=1.  Rebuild with this option");
 #endif
+    else
+      bb = new LCAOrbitalBuilder(targetPtcl, *ions, myComm, rootNode);
   }
 #endif //OHMMS_DIM==3
   PRE.flush();

src/QMCWaveFunctions/WaveFunctionComponent.h

@@ -492,7 +492,7 @@ struct WaveFunctionComponent : public QMCTraits
     }
   }
 
-  /** Calculates the derivatives of \f$ \grad(\textrm{log}(\psif)) \f$ with respect to
+  /** Calculates the derivatives of \f$ \nabla \textnormal{log} \psi_f \f$ with respect to
       the optimizable parameters, and the dot product of this is then
       performed with the passed-in G_in gradient vector. This object is then
       returned as dgradlogpsi.

src/QMCWaveFunctions/tests/CMakeLists.txt

@@ -61,7 +61,7 @@ SET(JASTROW_SRC test_bspline_jastrow.cpp test_counting_jastrow.cpp test_polynomi
                 test_rpa_jastrow.cpp test_user_jastrow.cpp test_kspace_jastrow.cpp test_pade_jastrow.cpp
                 test_short_range_cusp_jastrow.cpp test_DiffTwoBodyJastrowOrbital.cpp)
 SET(DETERMINANT_SRC FakeSPO.cpp makeRngSpdMatrix.cpp test_DiracDeterminantBatched.cpp test_DiracDeterminantBatched.cpp
-                    test_multi_dirac_determinant.cpp test_dirac_matrix.cpp
+                    test_multi_dirac_determinant.cpp test_dirac_matrix.cpp test_ci_configuration.cpp
                     test_multi_slater_determinant.cpp)
 
 IF(ENABLE_CUDA)

src/QMCWaveFunctions/tests/test_DiracDeterminant.cpp

@@ -463,6 +463,7 @@ TEST_CASE("DiracDeterminant_spinor_update", "[wavefunction][fermion]")
   elec_.spins[0]   = 0.0;
   elec_.spins[1]   = 0.2;
   elec_.spins[2]   = 0.4;
+  elec_.is_spinor_ = true;
 
   // O2 test example from pwscf non-collinear calculation.
   elec_.Lattice.R(0, 0) = 5.10509515;

src/QMCWaveFunctions/tests/test_MO_spinor.cpp

@@ -17,7 +17,7 @@
 #include "Particle/ParticleSet.h"
 #include "Particle/ParticleSetPool.h"
 #include "Particle/DistanceTableData.h"
-#include "QMCWaveFunctions/LCAO/LCAOSpinorBuilder.h"
+#include "QMCWaveFunctions/SPOSetBuilderFactory.h"
 
 namespace qmcplusplus
 {
@@ -53,6 +53,7 @@ void test_lcao_spinor()
   elec_.R[0][1]    = -0.3;
   elec_.R[0][2]    = 1.7;
   elec_.spins[0]   = 0.6;
+  elec_.is_spinor_ = true;
 
   SpeciesSet& tspecies       = elec_.getSpeciesSet();
   int upIdx                  = tspecies.addSpecies("u");
@@ -68,7 +69,7 @@ void test_lcao_spinor()
   ptcl.addParticleSet(std::move(ions_uptr));
 
   const char* particles = "<tmp> \
-   <sposet_builder name=\"spinorbuilder\" type=\"molecularspinor\" href=\"lcao_spinor.h5\" source=\"ion\" precision=\"float\"> \
+   <sposet_builder name=\"spinorbuilder\" type=\"molecularorbital\" href=\"lcao_spinor.h5\" source=\"ion\" precision=\"float\"> \
      <basisset transform=\"yes\"/> \
      <sposet name=\"myspo\" size=\"1\"/> \
    </sposet_builder> \
@@ -82,13 +83,14 @@ void test_lcao_spinor()
   xmlNodePtr root = doc.getRoot();
 
   xmlNodePtr bnode = xmlFirstElementChild(root);
-  LCAOSpinorBuilder bb(elec_, ions_, c, bnode);
+  SPOSetBuilderFactory fac(c, elec_, ptcl.getPool());
+  SPOSetBuilder* bb = fac.createSPOSetBuilder(bnode);
 
   // only pick up the last sposet
   std::unique_ptr<SPOSet> spo;
   processChildren(bnode, [&](const std::string& cname, const xmlNodePtr element) {
     if (cname == "sposet")
-      spo.reset(bb.createSPOSetFromXML(element));
+      spo.reset(bb->createSPOSet(element));
   });
   REQUIRE(spo);
 
@@ -232,6 +234,7 @@ void test_lcao_spinor_excited()
   elec_.R[0][1]    = -0.3;
   elec_.R[0][2]    = 1.7;
   elec_.spins[0]   = 0.6;
+  elec_.is_spinor_ = true;
 
   SpeciesSet& tspecies       = elec_.getSpeciesSet();
   int upIdx                  = tspecies.addSpecies("u");
@@ -247,7 +250,7 @@ void test_lcao_spinor_excited()
   ptcl.addParticleSet(std::move(ions_uptr));
 
   const char* particles = "<tmp> \
-   <sposet_builder name=\"spinorbuilder\" type=\"molecularspinor\" href=\"lcao_spinor.h5\" source=\"ion\" precision=\"float\"> \
+   <sposet_builder name=\"spinorbuilder\" type=\"molecularorbital\" href=\"lcao_spinor.h5\" source=\"ion\" precision=\"float\"> \
      <basisset name=\"myset\" transform=\"yes\"/> \
      <sposet name=\"myspo\" basisset=\"myset\" size=\"1\"> \
        <occupation mode=\"excited\"> \
@@ -265,13 +268,14 @@ void test_lcao_spinor_excited()
   xmlNodePtr root = doc.getRoot();
 
   xmlNodePtr bnode = xmlFirstElementChild(root);
-  LCAOSpinorBuilder bb(elec_, ions_, c, bnode);
+  SPOSetBuilderFactory fac(c, elec_, ptcl.getPool());
+  SPOSetBuilder* bb = fac.createSPOSetBuilder(bnode);
 
   // only pick up the last sposet
   std::unique_ptr<SPOSet> spo;
   processChildren(bnode, [&](const std::string& cname, const xmlNodePtr element) {
     if (cname == "sposet")
-      spo.reset(bb.createSPOSetFromXML(element));
+      spo.reset(bb->createSPOSet(element));
   });
   REQUIRE(spo);
 

src/QMCWaveFunctions/tests/test_ci_configuration.cpp

@@ -0,0 +1,75 @@
+//////////////////////////////////////////////////////////////////////////////////////
+// This file is distributed under the University of Illinois/NCSA Open Source License.
+// See LICENSE file in top directory for details.
+//
+// Copyright (c) 2021 QMCPACK developers.
+//
+// File developed by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory
+//
+// File created by: Ye Luo, yeluo@anl.gov, Argonne National Laboratory
+//////////////////////////////////////////////////////////////////////////////////////
+
+
+#include "catch.hpp"
+
+#include "Fermion/ci_configuration2.h"
+
+namespace qmcplusplus
+{
+TEST_CASE("ci_configuration2", "[wavefunction]")
+{
+  const int n_states = 6;
+  std::vector<size_t> ref{0, 1, 2, 3};
+  ci_configuration2 ref_state(ref);
+
+  size_t n_excited;
+  double sign = 0.0;
+  std::vector<size_t> pos(n_states);
+  std::vector<size_t> occupied(n_states);
+  std::vector<size_t> unoccupied(n_states);
+
+  std::vector<size_t> ext1{0, 1, 4, 5};
+  ci_configuration2 ext1_state(ext1);
+
+  n_excited = ext1_state.calculateNumOfExcitations(ref_state);
+  REQUIRE(n_excited == 2);
+
+  sign = ext1_state.calculateExcitations(ref_state, n_excited, pos, occupied, unoccupied);
+
+  REQUIRE(n_excited == 2);
+  CHECK(sign == 1.0);
+  CHECK(pos[0] == 2);
+  CHECK(pos[1] == 3);
+  CHECK(occupied[0] == 4);
+  CHECK(occupied[1] == 5);
+  CHECK(unoccupied[0] == 2);
+  CHECK(unoccupied[1] == 3);
+
+  std::vector<size_t> ext2{0, 1, 2, 6};
+  ci_configuration2 ext2_state(ext2);
+
+  n_excited = ext2_state.calculateNumOfExcitations(ref_state);
+  REQUIRE(n_excited == 1);
+
+  sign = ext2_state.calculateExcitations(ref_state, n_excited, pos, occupied, unoccupied);
+
+  REQUIRE(n_excited == 1);
+  CHECK(sign == 1.0);
+  CHECK(pos[0] == 3);
+  CHECK(occupied[0] == 6);
+  CHECK(unoccupied[0] == 3);
+
+  std::vector<size_t> ref2{0, 1};
+  ci_configuration2 ref2_state(ref2);
+  std::vector<size_t> ext3{1, 6};
+  ci_configuration2 ext3_state(ext3);
+
+  sign = ext3_state.calculateExcitations(ref2_state, n_excited, pos, occupied, unoccupied);
+  REQUIRE(n_excited == 1);
+  CHECK(sign == -1.0);
+
+  sign = ref2_state.calculateExcitations(ext3_state, n_excited, pos, occupied, unoccupied);
+  REQUIRE(n_excited == 1);
+  CHECK(sign == -1.0);
+}
+} // namespace qmcplusplus

src/QMCWaveFunctions/tests/test_einset_spinor.cpp

@@ -17,8 +17,7 @@
 #include "Particle/ParticleSet.h"
 #include "Particle/ParticleSetPool.h"
 #include "QMCWaveFunctions/WaveFunctionComponent.h"
-#include "QMCWaveFunctions/EinsplineSetBuilder.h"
+#include "QMCWaveFunctions/SPOSetBuilderFactory.h"
-#include "QMCWaveFunctions/EinsplineSpinorSetBuilder.h"
 
 #include <stdio.h>
 #include <string>
@@ -71,6 +70,7 @@ TEST_CASE("Einspline SpinorSet from HDF", "[wavefunction]")
   elec_.spins[0]   = 0.0;
   elec_.spins[1]   = 0.2;
   elec_.spins[2]   = 0.4;
+  elec_.is_spinor_ = true;
 
   // O2 test example from pwscf non-collinear calculation.
   elec_.Lattice.R(0, 0) = 5.10509515;
@@ -97,7 +97,7 @@ TEST_CASE("Einspline SpinorSet from HDF", "[wavefunction]")
 
 
   const char* particles = "<tmp> \
-   <sposet_builder name=\"A\" type=\"spinorbspline\" href=\"o2_45deg_spins.pwscf.h5\" tilematrix=\"1 0 0 0 1 0 0 0 1\" twistnum=\"0\" source=\"ion\" size=\"3\" precision=\"float\"> \
+   <sposet_builder name=\"A\" type=\"einspline\" href=\"o2_45deg_spins.pwscf.h5\" tilematrix=\"1 0 0 0 1 0 0 0 1\" twistnum=\"0\" source=\"ion\" size=\"3\" precision=\"float\"> \
      <sposet name=\"myspo\" size=\"3\"> \
        <occupation mode=\"ground\"/> \
      </sposet> \
@@ -113,8 +113,10 @@ TEST_CASE("Einspline SpinorSet from HDF", "[wavefunction]")
 
   xmlNodePtr ein1 = xmlFirstElementChild(root);
 
-  EinsplineSpinorSetBuilder einSet(elec_, ptcl.getPool(), c, ein1);
+  SPOSetBuilderFactory fac(c, elec_, ptcl.getPool());
-  std::unique_ptr<SPOSet> spo(einSet.createSPOSetFromXML(ein1));
+  SPOSetBuilder* builder = fac.createSPOSetBuilder(ein1);
+
+  std::unique_ptr<SPOSet> spo(builder->createSPOSet(ein1));
   REQUIRE(spo);
 
   SPOSet::ValueMatrix_t psiM(elec_.R.size(), spo->getOrbitalSetSize());

src/einspline/tests/test_3d.cpp

@@ -44,20 +44,20 @@ TEST_CASE("double_3d_natural","[einspline]")
   BCtype_d bc;
   bc.lCode = NATURAL;
   bc.rCode = NATURAL;
-  UBspline_3d_d* s = create_UBspline_3d_d(grid, grid, grid, bc, bc, bc, data);
+  auto s   = std::unique_ptr<UBspline_3d_d, void (*)(void*)>{create_UBspline_3d_d(grid, grid, grid, bc, bc, bc, data),
+                                                           destroy_Bspline};
 
   REQUIRE(s);
 
   double val;
-  eval_UBspline_3d_d(s, 1.0, 1.0, 1.0, &val);
+  eval_UBspline_3d_d(s.get(), 1.0, 1.0, 1.0, &val);
   REQUIRE(val == Approx(2.0));
   
   //This puts the value right below the limit
   double pos=10.0-5*std::numeric_limits<double>::epsilon();
-  eval_UBspline_3d_d(s, pos, pos, pos, &val);
+  eval_UBspline_3d_d(s.get(), pos, pos, pos, &val);
   REQUIRE(val == Approx(9.0));
 
-  destroy_Bspline(s);
 }
 
 TEST_CASE("double_3d_periodic","[einspline]")
@@ -86,44 +86,48 @@ TEST_CASE("double_3d_periodic","[einspline]")
   bc.lCode = PERIODIC;
   bc.rCode = PERIODIC;
 
-  UBspline_3d_d* s = create_UBspline_3d_d(grid, grid, grid, bc, bc, bc, data);
+  auto s = std::unique_ptr<UBspline_3d_d, void (*)(void*)>{create_UBspline_3d_d(grid, grid, grid, bc, bc, bc, data),
+                                                           destroy_Bspline};
   REQUIRE(s);
 
   double val;
-  eval_UBspline_3d_d(s, 0.0, 0.0, 0.0, &val);
+  eval_UBspline_3d_d(s.get(), 0.0, 0.0, 0.0, &val);
   REQUIRE(val == Approx(0.0));
 
   double pos=0.99999999;
-  eval_UBspline_3d_d(s, pos, pos, pos, &val);
+  eval_UBspline_3d_d(s.get(), pos, pos, pos, &val);
   REQUIRE(val == Approx(0.0));
 
-  eval_UBspline_3d_d(s, delta, delta, delta, &val);
+  eval_UBspline_3d_d(s.get(), delta, delta, delta, &val);
   REQUIRE(val == Approx(data[N*N + N + 1]));
 
-  multi_UBspline_3d_d* m = create_multi_UBspline_3d_d(grid, grid, grid, bc, bc, bc, 1);
+  auto m =
+      std::unique_ptr<multi_UBspline_3d_d, void (*)(void*)>{create_multi_UBspline_3d_d(grid, grid, grid, bc, bc, bc, 1),
+                                                            destroy_Bspline};
   REQUIRE(m);
 
-  set_multi_UBspline_3d_d(m, 0, data);
+  set_multi_UBspline_3d_d(m.get(), 0, data);
 
-  eval_multi_UBspline_3d_d(m, delta, delta, delta, &val);
+  eval_multi_UBspline_3d_d(m.get(), delta, delta, delta, &val);
   REQUIRE(val == Approx(data[N*N + N + 1]));
 
-
   BCtype_s bc_s;
   bc_s.lCode = PERIODIC;
   bc_s.rCode = PERIODIC;
 
-  multi_UBspline_3d_s* ms = create_multi_UBspline_3d_s(grid, grid, grid, bc_s, bc_s, bc_s, 1);
+  auto ms = std::unique_ptr<multi_UBspline_3d_s, void (*)(void*)>{create_multi_UBspline_3d_s(grid, grid, grid, bc_s,
+                                                                                             bc_s, bc_s, 1),
+                                                                  destroy_Bspline};
   REQUIRE(ms);
-  set_multi_UBspline_3d_s_d(ms, 0, data);
+  set_multi_UBspline_3d_s_d(ms.get(), 0, data);
 
   float fval;
-  eval_multi_UBspline_3d_s(ms, delta, delta, delta, &fval);
+  eval_multi_UBspline_3d_s(ms.get(), delta, delta, delta, &fval);
   REQUIRE(fval == Approx(data[N*N + N + 1]));
 
   float grads[3];
   float hess[9];
-  eval_multi_UBspline_3d_s_vgh(ms, 0.1, 0.2, 0.0, &fval, grads, hess);
+  eval_multi_UBspline_3d_s_vgh(ms.get(), 0.1, 0.2, 0.0, &fval, grads, hess);
 
   // See miniqmc-python/splines/test_3d.py for values
   REQUIRE(grads[0] == Approx(5.11104213833));

src/einspline/tests/test_one.cpp

@@ -38,26 +38,26 @@ TEST_CASE("double_1d_natural", "[einspline]")
   BCtype_d xBC;
   xBC.lCode = NATURAL;
   xBC.rCode = NATURAL;
-  UBspline_1d_d* s = create_UBspline_1d_d(x_grid, xBC, data.data());
+  auto s =
+      std::unique_ptr<UBspline_1d_d, void (*)(void*)>{create_UBspline_1d_d(x_grid, xBC, data.data()), destroy_Bspline};
 
   REQUIRE(s);
 
   double val;
 
-  eval_UBspline_1d_d(s, 1.0, &val);
+  eval_UBspline_1d_d(s.get(), 1.0, &val);
   REQUIRE(val == Approx(2.0));
 
-  eval_UBspline_1d_d(s, 9.9999999, &val);
+  eval_UBspline_1d_d(s.get(), 9.9999999, &val);
   REQUIRE(val == Approx(3.0));
 
   // This should assert
-  // eval_UBspline_1d_d(s, 10.0, &val);
+  // eval_UBspline_1d_d(s.get(), 10.0, &val);
   // REQUIRE(val == Approx(3.0));
 
-  eval_UBspline_1d_d(s, 5.5, &val);
+  eval_UBspline_1d_d(s.get(), 5.5, &val);
   REQUIRE(val == Approx(2.5));
 
-  destroy_Bspline(s);
 
   // three point case
   x_grid.start = 1.0;
@@ -69,21 +69,18 @@ TEST_CASE("double_1d_natural", "[einspline]")
 
   xBC.lCode = NATURAL;
   xBC.rCode = NATURAL;
-  s = create_UBspline_1d_d(x_grid, xBC, data.data());
+  s.reset(create_UBspline_1d_d(x_grid, xBC, data.data()));
 
   REQUIRE(s);
 
-  eval_UBspline_1d_d(s, 1.0, &val);
+  eval_UBspline_1d_d(s.get(), 1.0, &val);
   REQUIRE(val == Approx(2.0));
 
-  eval_UBspline_1d_d(s, 9.9999999, &val);
+  eval_UBspline_1d_d(s.get(), 9.9999999, &val);
   REQUIRE(val == Approx(3.0));
 
-  eval_UBspline_1d_d(s, 5.5, &val);
+  eval_UBspline_1d_d(s.get(), 5.5, &val);
   REQUIRE(val == Approx(2.7));
-
-  destroy_Bspline(s);
-
 }
 
 
@@ -101,13 +98,14 @@ TEST_CASE("double_1d_multi", "[einspline]")
   BCtype_d xBC;
   xBC.lCode = NATURAL;
   xBC.rCode = NATURAL;
-  multi_UBspline_1d_d* s = create_multi_UBspline_1d_d(x_grid, xBC, 1);
+  auto s    = std::unique_ptr<multi_UBspline_1d_d, void (*)(void*)>{create_multi_UBspline_1d_d(x_grid, xBC, 1),
+                                                                 destroy_Bspline};
   REQUIRE(s);
 
-  set_multi_UBspline_1d_d(s, 0, data);
+  set_multi_UBspline_1d_d(s.get(), 0, data);
 
   double val;
-  eval_multi_UBspline_1d_d(s, 1.0, &val);
+  eval_multi_UBspline_1d_d(s.get(), 1.0, &val);
   REQUIRE(val == Approx(2.0));
 }
 
@@ -133,15 +131,15 @@ TEST_CASE("double_1d_periodic", "[einspline]")
   bc.lCode = PERIODIC;
   bc.rCode = PERIODIC;
 
-  UBspline_1d_d* s = create_UBspline_1d_d(x_grid, bc, data);
+  auto s = std::unique_ptr<UBspline_1d_d, void (*)(void*)>{create_UBspline_1d_d(x_grid, bc, data), destroy_Bspline};
 
   REQUIRE(s);
 
   double val;
-  eval_UBspline_1d_d(s, 0.0, &val);
+  eval_UBspline_1d_d(s.get(), 0.0, &val);
   REQUIRE(val == Approx(0.0));
 
-  eval_UBspline_1d_d(s, delta, &val);
+  eval_UBspline_1d_d(s.get(), delta, &val);
   REQUIRE(val == Approx(data[1]));
 
   double micro_delta = delta / 4.0;
@@ -152,16 +150,14 @@ TEST_CASE("double_1d_periodic", "[einspline]")
     double x = micro_delta * i;
     micro_data[i] = sin(tpi*x);
   }
-  eval_UBspline_1d_d(s, micro_delta * 3, &val);
+  eval_UBspline_1d_d(s.get(), micro_delta * 3, &val);
   REQUIRE(val == Approx(micro_data[3]).epsilon(0.001));
 
-  eval_UBspline_1d_d(s, micro_delta * 17, &val);
+  eval_UBspline_1d_d(s.get(), micro_delta * 17, &val);
   REQUIRE(val == Approx(micro_data[17]).epsilon(0.001));
 
-  eval_UBspline_1d_d(s, micro_delta * 31, &val);
+  eval_UBspline_1d_d(s.get(), micro_delta * 31, &val);
   REQUIRE(val == Approx(micro_data[31]).epsilon(0.001));
-
-  destroy_Bspline(s);
 }
 
 #ifdef QMC_CUDA
@@ -182,23 +178,24 @@ TEST_CASE("multi_cuda_wrapper", "[einspline]")
   BCtype_s xBC;
   xBC.lCode = NATURAL;
   xBC.rCode = NATURAL;
-  multi_UBspline_1d_s* s = create_multi_UBspline_1d_s(x_grid, xBC, 1);
+  auto s    = std::unique_ptr<multi_UBspline_1d_s, void (*)(void*)>{create_multi_UBspline_1d_s(x_grid, xBC, 1),
+                                                                 destroy_Bspline};
   REQUIRE(s);
-  set_multi_UBspline_1d_s(s, 0, data);
+  set_multi_UBspline_1d_s(s.get(), 0, data);
 
   float pos[1];
   pos[0] = 0.0;
 
   // Check the CPU value
   float cpu_val[1];
-  eval_multi_UBspline_1d_s(s, pos[0], cpu_val);
+  eval_multi_UBspline_1d_s(s.get(), pos[0], cpu_val);
   REQUIRE(cpu_val[0] == 2.0);
 
   //this would assert in debug and is an illegal value for pos[0]
   //pos[0] = 11.0;
   pos[0] = 9.99999999;
   // Check the CPU value
-  eval_multi_UBspline_1d_s(s, pos[0], cpu_val);
+  eval_multi_UBspline_1d_s(s.get(), pos[0], cpu_val);
   //std::cout << std::setprecision(24) << cpu_val[0] << " == " << 3.0000f << '\n';
   //With power 9/ clang 8  3.0000f is 3 but cpu_val[0] = 3.0000002384185791015625
   REQUIRE(cpu_val[0] == Approx(3.0000f).epsilon(0.00000025));
@@ -206,11 +203,11 @@ TEST_CASE("multi_cuda_wrapper", "[einspline]")
   // Check the GPU value
   pos[0] = 0.0;
   float vals_output[1];
-  test_multi(s, pos, vals_output);
+  test_multi(s.get(), pos, vals_output);
   REQUIRE(vals_output[0] == 2.0);
 
   pos[0] = 9.99999999;
-  test_multi(s, pos, vals_output);
+  test_multi(s.get(), pos, vals_output);
   REQUIRE(vals_output[0] == 3.0);
 }
 

tests/solids/monoO_noncollinear_1x1x1_pp/qmc_short.in.xml

@@ -16,7 +16,7 @@
          </parameter>
          <parameter name="LR_dim_cutoff"       >    15                 </parameter>
       </simulationcell>
-      <particleset name="e">
+      <particleset name="e" spinor="yes">
          <group name="u" size="12" mass="1.0">
             <parameter name="charge"              >    -1                    </parameter>
             <parameter name="mass"                >    1.0                   </parameter>
@@ -63,7 +63,7 @@
          </group>
       </particleset>
       <wavefunction name="psi0" target="e">
-   <sposet_builder name="A" type="spinorbspline" href="eshdf.h5" tilematrix="1 0 0 0 1 0 0 0 1" twistnum="0" source="ion0" size="12"> 
+   <sposet_builder name="A" type="einspline" href="o2_45deg_spins.pwscf.h5" tilematrix="1 0 0 0 1 0 0 0 1" twistnum="0" source="ion0" size="12"> 
      <sposet name="myspo" size="12"> 
        <occupation mode="ground"/> 
      </sposet>