diff --git a/src/QMCWaveFunctions/Determinant.h b/src/QMCWaveFunctions/Determinant.h
index 9f8b547a4..862057fa5 100644
--- a/src/QMCWaveFunctions/Determinant.h
+++ b/src/QMCWaveFunctions/Determinant.h
@@ -178,7 +178,8 @@ void checkDiff(const MT1 &a, const MT2 &b, const std::string &tag)
 
 struct DiracDeterminant : public WaveFunctionComponentBase
 {
-  DiracDeterminant(int nels, RandomGenerator<RealType> RNG)
+  DiracDeterminant(int nels, RandomGenerator<RealType> RNG, int First=0)
+  : FirstIndex(First)
   {
     psiMinv.resize(nels, nels);
     psiV.resize(nels);
@@ -222,7 +223,7 @@ struct DiracDeterminant : public WaveFunctionComponentBase
 
   GradType evalGrad(ParticleSet &P, int iat) {}
 
-  ValueType ratioGrad(ParticleSet &P, int iat, GradType &grad) {}
+  ValueType ratioGrad(ParticleSet &P, int iat, GradType &grad) { return ratio(P, iat); }
 
   void evaluateGL(ParticleSet &P, ParticleSet::ParticleGradient_t &G,
                   ParticleSet::ParticleLaplacian_t &L, bool fromscratch = false) {}
@@ -246,7 +247,7 @@ struct DiracDeterminant : public WaveFunctionComponentBase
     constexpr double czero(0);
     for (int j = 0; j < nels; ++j)
       psiV[j] = myRandom() - shift;
-    curRatio = inner_product_n(psiV.data(), psiMinv[iel], nels, czero);
+    curRatio = inner_product_n(psiV.data(), psiMinv[iel-FirstIndex], nels, czero);
     return curRatio;
   }
 
@@ -254,8 +255,8 @@ struct DiracDeterminant : public WaveFunctionComponentBase
   inline void acceptMove(ParticleSet &P, int iel)
   {
     const int nels = psiV.size();
-    updateRow(psiMinv.data(), psiV.data(), nels, nels, iel, curRatio);
-    std::copy_n(psiV.data(), nels, psiMsave[iel]);
+    updateRow(psiMinv.data(), psiV.data(), nels, nels, iel-FirstIndex, curRatio);
+    std::copy_n(psiV.data(), nels, psiMsave[iel-FirstIndex]);
   }
 
   /** accessor functions for checking */
@@ -269,6 +270,8 @@ struct DiracDeterminant : public WaveFunctionComponentBase
   double curRatio;
   /// workspace size
   int LWork;
+  /// initial particle index
+  const int FirstIndex;
   /// inverse matrix to be update
   Matrix<RealType> psiMinv;
   /// a SPO set for the row update
diff --git a/src/QMCWaveFunctions/DeterminantRef.h b/src/QMCWaveFunctions/DeterminantRef.h
index e1531eb54..a766e5485 100644
--- a/src/QMCWaveFunctions/DeterminantRef.h
+++ b/src/QMCWaveFunctions/DeterminantRef.h
@@ -182,7 +182,8 @@ struct DiracDeterminantRef : public qmcplusplus::WaveFunctionComponentBase
 {
   using ParticleSet = qmcplusplus::ParticleSet;
 
-  DiracDeterminantRef(int nels, qmcplusplus::RandomGenerator<RealType> RNG)
+  DiracDeterminantRef(int nels, qmcplusplus::RandomGenerator<RealType> RNG, int First=0)
+  : FirstIndex(First)
   {
     psiMinv.resize(nels, nels);
     psiV.resize(nels);
@@ -225,7 +226,7 @@ struct DiracDeterminantRef : public qmcplusplus::WaveFunctionComponentBase
     // FIXME do we want remainder of evaluateLog?
   }
   GradType evalGrad(ParticleSet &P, int iat) {}
-  ValueType ratioGrad(ParticleSet &P, int iat, GradType &grad) {}
+  ValueType ratioGrad(ParticleSet &P, int iat, GradType &grad) { return ratio(P, iat); }
   void evaluateGL(ParticleSet &P, ParticleSet::ParticleGradient_t &G,
                   ParticleSet::ParticleLaplacian_t &L, bool fromscratch = false) {}
 
@@ -248,7 +249,7 @@ struct DiracDeterminantRef : public qmcplusplus::WaveFunctionComponentBase
     constexpr double czero(0);
     for (int j = 0; j < nels; ++j)
       psiV[j] = myRandom() - shift;
-    curRatio = inner_product_n(psiV.data(), psiMinv[iel], nels, czero);
+    curRatio = inner_product_n(psiV.data(), psiMinv[iel-FirstIndex], nels, czero);
     return curRatio;
   }
 
@@ -256,8 +257,8 @@ struct DiracDeterminantRef : public qmcplusplus::WaveFunctionComponentBase
   inline void acceptMove(ParticleSet &P, int iel)
   {
     const int nels = psiV.size();
-    updateRow(psiMinv.data(), psiV.data(), nels, nels, iel, curRatio);
-    std::copy_n(psiV.data(), nels, psiMsave[iel]);
+    updateRow(psiMinv.data(), psiV.data(), nels, nels, iel-FirstIndex, curRatio);
+    std::copy_n(psiV.data(), nels, psiMsave[iel-FirstIndex]);
   }
 
   /** accessor functions for checking */
@@ -271,6 +272,8 @@ struct DiracDeterminantRef : public qmcplusplus::WaveFunctionComponentBase
   double curRatio;
   /// workspace size
   int LWork;
+  /// initial particle index
+  const int FirstIndex;
   /// inverse matrix to be update
   qmcplusplus::Matrix<RealType> psiMinv;
   /// a SPO set for the row update
diff --git a/src/QMCWaveFunctions/WaveFunction.cpp b/src/QMCWaveFunctions/WaveFunction.cpp
index d4f602e27..cadf628c7 100644
--- a/src/QMCWaveFunctions/WaveFunction.cpp
+++ b/src/QMCWaveFunctions/WaveFunction.cpp
@@ -33,7 +33,9 @@ WaveFunction::WaveFunction(ParticleSet &ions, ParticleSet &els,
   d_ie = DistanceTable::add(ions, els, DT_SOA);
 
   // determinant component
-  Det = new DetType(els.getTotalNum(), RNG);
+  nelup = els.getTotalNum()/2;
+  Det_up = new DetType(nelup, RNG, 0);
+  Det_dn = new DetType(els.getTotalNum()-nelup, RNG, nelup);
 
   // J1 component
   J1 = new J1OrbType(ions, els);
@@ -48,7 +50,14 @@ WaveFunction::WaveFunction(ParticleSet &ions, ParticleSet &els,
   buildJeeI(*J3, els.Lattice.WignerSeitzRadius);
 }
 
-WaveFunction::~WaveFunction() { delete J2; }
+WaveFunction::~WaveFunction()
+{
+  delete Det_up;
+  delete Det_dn;
+  delete J1;
+  delete J2;
+  delete J3;
+}
 
 void WaveFunction::evaluateLog(ParticleSet &P)
 {
@@ -57,7 +66,8 @@ void WaveFunction::evaluateLog(ParticleSet &P)
   {
     P.G       = czero;
     P.L       = czero;
-    LogValue  = Det->evaluateLog(P, P.G, P.L);
+    LogValue  = Det_up->evaluateLog(P, P.G, P.L);
+    LogValue += Det_dn->evaluateLog(P, P.G, P.L);
     LogValue += J1->evaluateLog(P, P.G, P.L);
     LogValue += J2->evaluateLog(P, P.G, P.L);
     LogValue += J3->evaluateLog(P, P.G, P.L);
@@ -67,26 +77,35 @@ void WaveFunction::evaluateLog(ParticleSet &P)
 
 WaveFunctionBase::posT WaveFunction::evalGrad(ParticleSet &P, int iat)
 {
-  return Det->evalGrad(P, iat) + J1->evalGrad(P, iat) + J2->evalGrad(P, iat) +
-         J3->evalGrad(P, iat);
+  return ( iat<nelup ? Det_up->evalGrad(P, iat) : Det_dn->evalGrad(P, iat) )
+         + J1->evalGrad(P, iat)
+         + J2->evalGrad(P, iat)
+         + J3->evalGrad(P, iat);
 }
 
 WaveFunctionBase::valT WaveFunction::ratioGrad(ParticleSet &P, int iat,
                                                posT &grad)
 {
-  return Det->ratioGrad(P, iat, grad) + J1->ratioGrad(P, iat, grad) +
-         J2->ratioGrad(P, iat, grad) + J3->ratioGrad(P, iat, grad);
+  return ( iat<nelup ? Det_up->ratioGrad(P, iat, grad) : Det_dn->ratioGrad(P, iat, grad) )
+         * J1->ratioGrad(P, iat, grad)
+         * J2->ratioGrad(P, iat, grad)
+         * J3->ratioGrad(P, iat, grad);
 }
 
 WaveFunctionBase::valT WaveFunction::ratio(ParticleSet &P, int iat)
 {
-  return Det->ratio(P, iat) * J1->ratio(P, iat) * J2->ratio(P, iat) *
-         J3->ratio(P, iat);
+  return ( iat<nelup ? Det_up->ratio(P, iat) : Det_dn->ratio(P, iat) )
+         * J1->ratio(P, iat)
+         * J2->ratio(P, iat)
+         * J3->ratio(P, iat);
 }
 
 void WaveFunction::acceptMove(ParticleSet &P, int iat)
 {
-  Det->acceptMove(P, iat);
+  if(iat<nelup)
+    Det_up->acceptMove(P, iat);
+  else
+    Det_dn->acceptMove(P, iat);
   J1->acceptMove(P, iat);
   J2->acceptMove(P, iat);
   J3->acceptMove(P, iat);
@@ -99,7 +118,8 @@ void WaveFunction::evaluateGL(ParticleSet &P)
   constexpr valT czero(0);
   P.G = czero;
   P.L = czero;
-  Det->evaluateGL(P, P.G, P.L);
+  Det_up->evaluateGL(P, P.G, P.L);
+  Det_dn->evaluateGL(P, P.G, P.L);
   J1->evaluateGL(P, P.G, P.L);
   J2->evaluateGL(P, P.G, P.L);
   J3->evaluateGL(P, P.G, P.L);
diff --git a/src/QMCWaveFunctions/WaveFunction.h b/src/QMCWaveFunctions/WaveFunction.h
index f8209c902..1c9a25fba 100644
--- a/src/QMCWaveFunctions/WaveFunction.h
+++ b/src/QMCWaveFunctions/WaveFunction.h
@@ -73,7 +73,9 @@ struct WaveFunction : public WaveFunctionBase
   J1OrbType *J1;
   J2OrbType *J2;
   J3OrbType *J3;
-  DetType *Det;
+  int nelup;
+  DetType *Det_up;
+  DetType *Det_dn;
 
   WaveFunction(ParticleSet &ions, ParticleSet &els,
                RandomGenerator<RealType> RNG);
@@ -106,7 +108,9 @@ struct WaveFunctionRef : public qmcplusplus::WaveFunctionBase
   J1OrbType *J1;
   J2OrbType *J2;
   J3OrbType *J3;
-  DetType *Det;
+  int nelup;
+  DetType *Det_up;
+  DetType *Det_dn;
   PooledData<valT> Buffer;
 
   WaveFunctionRef(ParticleSet &ions, ParticleSet &els,
diff --git a/src/QMCWaveFunctions/WaveFunctionRef.cpp b/src/QMCWaveFunctions/WaveFunctionRef.cpp
index a0e9e9f41..4e2a878f3 100644
--- a/src/QMCWaveFunctions/WaveFunctionRef.cpp
+++ b/src/QMCWaveFunctions/WaveFunctionRef.cpp
@@ -39,7 +39,9 @@ WaveFunctionRef::WaveFunctionRef(ParticleSet &ions, ParticleSet &els,
   d_ie = DistanceTable::add(ions, els, DT_SOA);
 
   // determinant component
-  Det = new DetType(els.getTotalNum(), RNG);
+  nelup = els.getTotalNum()/2;
+  Det_up = new DetType(nelup, RNG, 0);
+  Det_dn = new DetType(els.getTotalNum()-nelup, RNG, nelup);
 
   // J1 component
   J1 = new J1OrbType(ions, els);
@@ -54,7 +56,14 @@ WaveFunctionRef::WaveFunctionRef(ParticleSet &ions, ParticleSet &els,
   buildJeeI(*J3, els.Lattice.WignerSeitzRadius);
 }
 
-WaveFunctionRef::~WaveFunctionRef() { delete J2; }
+WaveFunctionRef::~WaveFunctionRef()
+{
+  delete Det_up;
+  delete Det_dn;
+  delete J1;
+  delete J2;
+  delete J3;
+}
 
 void WaveFunctionRef::evaluateLog(ParticleSet &P)
 {
@@ -63,7 +72,8 @@ void WaveFunctionRef::evaluateLog(ParticleSet &P)
   {
     P.G       = czero;
     P.L       = czero;
-    LogValue  = Det->evaluateLog(P, P.G, P.L);
+    LogValue  = Det_up->evaluateLog(P, P.G, P.L);
+    LogValue += Det_dn->evaluateLog(P, P.G, P.L);
     LogValue += J1->evaluateLog(P, P.G, P.L);
     LogValue += J2->evaluateLog(P, P.G, P.L);
     LogValue += J3->evaluateLog(P, P.G, P.L);
@@ -73,26 +83,35 @@ void WaveFunctionRef::evaluateLog(ParticleSet &P)
 
 WaveFunctionBase::posT WaveFunctionRef::evalGrad(ParticleSet &P, int iat)
 {
-  return Det->evalGrad(P, iat) + J1->evalGrad(P, iat) + J2->evalGrad(P, iat) +
-         J3->evalGrad(P, iat);
+  return ( iat<nelup ? Det_up->evalGrad(P, iat) : Det_dn->evalGrad(P, iat) )
+         + J1->evalGrad(P, iat)
+         + J2->evalGrad(P, iat)
+         + J3->evalGrad(P, iat);
 }
 
 WaveFunctionBase::valT WaveFunctionRef::ratioGrad(ParticleSet &P, int iat,
                                                   posT &grad)
 {
-  return Det->ratioGrad(P, iat, grad) + J1->ratioGrad(P, iat, grad) +
-         J2->ratioGrad(P, iat, grad) + J3->ratioGrad(P, iat, grad);
+  return ( iat<nelup ? Det_up->ratioGrad(P, iat, grad) : Det_dn->ratioGrad(P, iat, grad) )
+         * J1->ratioGrad(P, iat, grad)
+         * J2->ratioGrad(P, iat, grad)
+         * J3->ratioGrad(P, iat, grad);
 }
 
 WaveFunctionBase::valT WaveFunctionRef::ratio(ParticleSet &P, int iat)
 {
-  return Det->ratio(P, iat) * J1->ratio(P, iat) * J2->ratio(P, iat) *
-         J3->ratio(P, iat);
+  return ( iat<nelup ? Det_up->ratio(P, iat) : Det_dn->ratio(P, iat) )
+         * J1->ratio(P, iat)
+         * J2->ratio(P, iat)
+         * J3->ratio(P, iat);
 }
 
 void WaveFunctionRef::acceptMove(ParticleSet &P, int iat)
 {
-  Det->acceptMove(P, iat);
+  if(iat<nelup)
+    Det_up->acceptMove(P, iat);
+  else
+    Det_dn->acceptMove(P, iat);
   J1->acceptMove(P, iat);
   J2->acceptMove(P, iat);
   J3->acceptMove(P, iat);
@@ -105,7 +124,8 @@ void WaveFunctionRef::evaluateGL(ParticleSet &P)
   constexpr valT czero(0);
   P.G = czero;
   P.L = czero;
-  Det->evaluateGL(P, P.G, P.L);
+  Det_up->evaluateGL(P, P.G, P.L);
+  Det_dn->evaluateGL(P, P.G, P.L);
   J1->evaluateGL(P, P.G, P.L);
   J2->evaluateGL(P, P.G, P.L);
   J3->evaluateGL(P, P.G, P.L);