diff --git a/src/QMCWaveFunctions/EinsplineSet.cpp b/src/QMCWaveFunctions/EinsplineSet.cpp
index b0cd052d3..f1e3f1df1 100644
--- a/src/QMCWaveFunctions/EinsplineSet.cpp
+++ b/src/QMCWaveFunctions/EinsplineSet.cpp
@@ -288,11 +288,6 @@ namespace qmcplusplus {
   EinsplineSetExtended<StorageType>::setOrbitalSetSize(int norbs)
   {
     OrbitalSetSize = norbs;
-    eikr.resize(norbs);
-    phase.resize(norbs);
-    StorageValueVector.resize(norbs);
-    StorageGradVector.resize(norbs);
-    StorageHessVector.resize(norbs);
   }
   
   template<typename StorageType> void
@@ -304,14 +299,23 @@ namespace qmcplusplus {
     for (int i=0; i<OHMMS_DIM; i++)
       ru[i] -= std::floor (ru[i]);
     EinsplineMultiEval (MultiSpline, ru, StorageValueVector);
-    computePhaseFactors(r);
-    for (int i=0; i<psi.size(); i++) {
+    //computePhaseFactors(r);
+    int N = StorageValueVector.size();
+    int psiIndex = 0;
+    for (int i=0; i<N; i++) {
       PosType k = kPoints[i];
       double s,c;
       double phase = -dot(r, k);
       sincos (phase, &s, &c);
       complex<double> e_mikr (c,s);
-      convert (e_mikr*StorageValueVector[i], psi[i]);
+      complex<double> psi_val = e_mikr*StorageValueVector[i];
+      psi[psiIndex] = real(psi_val);
+      psiIndex++;
+      if (MakeTwoCopies[i]) {
+	psi[psiIndex] = imag(psi_val);
+	psiIndex++;
+      }
+      //convert (e_mikr*StorageValueVector[i], psi[i]);
     }
   }
 
@@ -325,7 +329,7 @@ namespace qmcplusplus {
     for (int i=0; i<OHMMS_DIM; i++)
       ru[i] -= std::floor (ru[i]);
     EinsplineMultiEval (MultiSpline, ru, StorageValueVector);
-    computePhaseFactors(r);
+    //computePhaseFactors(r);
     for (int i=0; i<psi.size(); i++) {
       PosType k = kPoints[i];
       double s,c;
@@ -362,7 +366,7 @@ namespace qmcplusplus {
       ru[i] -= std::floor (ru[i]);
     EinsplineMultiEval (MultiSpline, ru, StorageValueVector,
 			StorageGradVector, StorageHessVector);
-    computePhaseFactors(r);
+    //computePhaseFactors(r);
     complex<double> eye (0.0, 1.0);
     for (int j=0; j<psi.size(); j++) {
       complex<double> u, laplu;
@@ -397,7 +401,7 @@ namespace qmcplusplus {
       ru[i] -= std::floor (ru[i]);
     EinsplineMultiEval (MultiSpline, ru, StorageValueVector,
 			StorageGradVector, StorageHessVector);
-    computePhaseFactors(r);
+    //computePhaseFactors(r);
     complex<double> eye (0.0, 1.0);
     for (int j=0; j<psi.size(); j++) {
       complex<double> u, laplu;
@@ -451,9 +455,11 @@ namespace qmcplusplus {
 	ru[n] -= std::floor (ru[n]);
       EinsplineMultiEval (MultiSpline, ru, StorageValueVector,
 			  StorageGradVector, StorageHessVector);
-      computePhaseFactors(r);
+      //computePhaseFactors(r);
       complex<double> eye (0.0, 1.0);
-      for (int j=0; j<OrbitalSetSize; j++) {
+      int N = StorageValueVector.size();
+      int psiIndex = 0;
+      for (int j=0; j<N; j++) {
 	complex<double> u, laplu;
 	TinyVector<complex<double>, OHMMS_DIM> gradu;
 	u = StorageValueVector[j];
@@ -467,14 +473,35 @@ namespace qmcplusplus {
 	double phase = -dot(r, k);
 	sincos (phase, &s, &c);
 	complex<double> e_mikr (c,s);
-	psi(j,i) = real(e_mikr*u);
+
+	complex<double> psi_val = e_mikr*u;
+	TinyVector<complex<double>,OHMMS_DIM> psi_grad =
+	  e_mikr*(-eye*u*ck + gradu);
+	complex<double> psi_lapl = e_mikr*(-dot(k,k)*u - 2.0*eye*dot(ck,gradu) + laplu);
+	
+	psi(psiIndex,i) = real(psi_val);
 	for (int n=0; n<3; n++)
-	  dpsi(i,j) = real(e_mikr*(-eye*u*ck + gradu));
-	d2psi(i,j) = real(e_mikr*(-dot(k,k)*u - 2.0*eye*dot(ck,gradu) + laplu));
+	  dpsi(i,psiIndex)[n] = real(psi_grad[n]);
+	d2psi(i,psiIndex) = real(psi_lapl);
+	psiIndex++;
+
+	if (MakeTwoCopies[j]) {
+	  psi(psiIndex,i) = imag(psi_val);
+	  for (int n=0; n<3; n++)
+	    dpsi(i,psiIndex)[n] = imag(psi_grad[n]);
+	  d2psi(i,psiIndex) = imag(psi_lapl);
+	  psiIndex++;
+	}
+
+	// psi(j,i) = real(e_mikr*u);
+	// for (int n=0; n<3; n++)
+	//   dpsi(i,j)[n] = real(e_mikr*(-eye*u*ck + gradu));
+	// d2psi(i,j) = real(e_mikr*(-dot(k,k)*u - 2.0*eye*dot(ck,gradu) + laplu));
 	//convert(e_mikr * u, psi(j,i));
 	//convertVec(e_mikr*(-eye*u*ck + gradu), dpsi(i,j));
 	//convert(e_mikr*(-dot(k,k)*u - 2.0*eye*dot(ck,gradu) + laplu), d2psi(i,j));
       } 
+      //cerr << "VGH psiIndex = " << psiIndex << endl;
     }
   }
   
@@ -492,7 +519,7 @@ namespace qmcplusplus {
 	ru[n] -= std::floor (ru[n]);
       EinsplineMultiEval (MultiSpline, ru, StorageValueVector,
 			  StorageGradVector, StorageHessVector);
-      computePhaseFactors(r);
+      //computePhaseFactors(r);
       complex<double> eye (0.0, 1.0);
       for (int j=0; j<OrbitalSetSize; j++) {
 	complex<double> u, laplu;
diff --git a/src/QMCWaveFunctions/EinsplineSet.h b/src/QMCWaveFunctions/EinsplineSet.h
index bf659fff7..0b426e92e 100644
--- a/src/QMCWaveFunctions/EinsplineSet.h
+++ b/src/QMCWaveFunctions/EinsplineSet.h
@@ -179,7 +179,7 @@ namespace qmcplusplus {
     StorageGradVector_t  StorageGradVector;
     StorageHessVector_t  StorageHessVector;
     // True if we should unpack this orbital into two copies
-    Vector<bool>         MakeTwoCopies;
+    vector<bool>         MakeTwoCopies;
     // k-points for each orbital
     Vector<TinyVector<double,OHMMS_DIM> > kPoints;
 
@@ -221,12 +221,12 @@ namespace qmcplusplus {
   inline void EinsplineSetExtended<StorageType>::computePhaseFactors(TinyVector<double,OHMMS_DIM> r)
   {
 #ifdef HAVE_MKL
-    for (int i=0; i<OrbitalSetSize; i++) 
+    for (int i=0; i<kPoints.size(); i++) 
       phase[i] = -dot(r, kPoints[i]);
     vzCIS(OrbitalSetSize, phase, (double*)eikr.data());
 #else
     double s, c;
-    for (int i=0; i<OrbitalSetSize; i++) {
+    for (int i=0; i<kPoints.size(); i++) {
       phase[i] = -dot(r, kPoints[i]);
       sincos (phase[i], &s, &c);
       eikr[i] = complex<double>(c,s);
diff --git a/src/QMCWaveFunctions/EinsplineSetBuilder.cpp b/src/QMCWaveFunctions/EinsplineSetBuilder.cpp
index 4138d2a58..270c6fdea 100644
--- a/src/QMCWaveFunctions/EinsplineSetBuilder.cpp
+++ b/src/QMCWaveFunctions/EinsplineSetBuilder.cpp
@@ -464,8 +464,10 @@ namespace qmcplusplus {
     for (int i=0; i<superSets[TwistNum].size(); i++)
       IncludeTwists.push_back(superSets[TwistNum][i]);
 
-    vector<int> copyTwists;
     // Now, find out which twists are distinct
+    DistinctTwists.clear();
+#ifndef QMC_COMPLEX
+    vector<int> copyTwists;
     for (int i=0; i<IncludeTwists.size(); i++) {
       int ti        = IncludeTwists[i];
       PosType twist_i = TwistAngles[ti];
@@ -475,11 +477,8 @@ namespace qmcplusplus {
 	PosType twist_j = TwistAngles[tj];
 	PosType sum  = twist_i + twist_j;
 	PosType diff = twist_i - twist_j;
-	if (TwistPair (twist_i, twist_j)) {
-	  cerr << "twist " << twist_i << " is paired with twist "
-	       << twist_j << endl;
+	if (TwistPair (twist_i, twist_j)) 
 	  distinct = false;
-	}
       }
       if (distinct)
 	DistinctTwists.push_back(ti);
@@ -501,6 +500,15 @@ namespace qmcplusplus {
       fprintf (stderr, "Using %d copies of twist angle [%6.3f, %6.3f, %6.3f]\n",
 	       MakeTwoCopies[i] ? 2 : 1, twist_i[0], twist_i[1], twist_i[2]);
     }
+#else
+    DistinctTwists.resize(IncludeTwists.size());
+    MakeTwoCopies.resize(IncludeTwists.size());
+    for (int i=0; i<IncludeTwists.size(); i++) {
+      DistinctTwists[i] = IncludeTwists[i];
+      MakeTwoCopies[i] = false;
+    }
+
+#endif
   }
   
   // This function analyzes the twist vectors to see if they lay on a
@@ -624,12 +632,13 @@ namespace qmcplusplus {
       eigenstatesGroup = "/eigenstates";
 
     SortBands.clear();
-    for (int ti=0; ti<IncludeTwists.size(); ti++) {
-      int tindex = IncludeTwists[ti];
+    for (int ti=0; ti<DistinctTwists.size(); ti++) {
+      int tindex = DistinctTwists[ti];
       for (int bi=0; bi<NumBands; bi++) {
 	BandInfo band;
 	band.TwistIndex = tindex;
 	band.BandIndex  = bi;
+	band.MakeTwoCopies = MakeTwoCopies[ti];
 	
 	// Read eigenenergy from file
 	ostringstream ePath, sPath;
@@ -664,6 +673,18 @@ namespace qmcplusplus {
       cerr << "Sorting the bands now:\n";
       sort (SortBands.begin(), SortBands.end());
     }
+    
+    int orbIndex = 0;
+    int numOrbs = 0;
+    while (numOrbs < OrbitalSet->getOrbitalSetSize()) {
+      if (SortBands[orbIndex].MakeTwoCopies)
+	numOrbs += 2;
+      else
+	numOrbs++;
+      orbIndex++;
+    }
+    NumDistinctOrbitals = orbIndex;
+    cerr << "We will read " << NumDistinctOrbitals << " distinct orbitals.\n";
   }
 
 
@@ -749,13 +770,21 @@ namespace qmcplusplus {
   EinsplineSetBuilder::ReadBands 
   (int spin, EinsplineSetExtended<double>* orbitalSet)
   {
-    int N = orbitalSet->getOrbitalSetSize();
+    //int N = orbitalSet->getOrbitalSetSize();
+    int N = NumDistinctOrbitals;
+    orbitalSet->kPoints.resize(N);
+    orbitalSet->MakeTwoCopies.resize(N);
+    orbitalSet->StorageValueVector.resize(N);
+    orbitalSet->StorageGradVector.resize(N);
+    orbitalSet->StorageHessVector.resize(N);
+    orbitalSet->phase.resize(N);
+    orbitalSet->eikr.resize(N);
     // Read in k-points
-    orbitalSet->kPoints.resize(orbitalSet->getOrbitalSetSize());
     for (int iorb=0; iorb<N; iorb++) {
       int ti = SortBands[iorb].TwistIndex;
       PosType twist  = TwistAngles[ti];
       orbitalSet->kPoints[iorb] = orbitalSet->PrimLattice.k_cart(twist);
+      orbitalSet->MakeTwoCopies[iorb] = SortBands[iorb].MakeTwoCopies;
     }
     
     // First, check to see if we have already read this in
@@ -876,13 +905,22 @@ namespace qmcplusplus {
   EinsplineSetBuilder::ReadBands 
   (int spin, EinsplineSetExtended<complex<double> >* orbitalSet)
   {
-    int N = orbitalSet->getOrbitalSetSize();
+    // int N = orbitalSet->getOrbitalSetSize();
+    int N = NumDistinctOrbitals;
+    cerr << "NumDistinctOrbitals = " << N << endl;
+    orbitalSet->kPoints.resize(N);
+    orbitalSet->MakeTwoCopies.resize(N);
+    orbitalSet->StorageValueVector.resize(N);
+    orbitalSet->StorageGradVector.resize(N);
+    orbitalSet->StorageHessVector.resize(N);
+    orbitalSet->phase.resize(N);
+    orbitalSet->eikr.resize(N);
     // Read in k-points
-    orbitalSet->kPoints.resize(orbitalSet->getOrbitalSetSize());
     for (int iorb=0; iorb<N; iorb++) {
       int ti = SortBands[iorb].TwistIndex;
       PosType twist  = TwistAngles[ti];
       orbitalSet->kPoints[iorb] = orbitalSet->PrimLattice.k_cart(twist);
+      orbitalSet->MakeTwoCopies[iorb] = SortBands[iorb].MakeTwoCopies;
     }
     
     // First, check to see if we have already read this in
diff --git a/src/QMCWaveFunctions/EinsplineSetBuilder.h b/src/QMCWaveFunctions/EinsplineSetBuilder.h
index 51dce8f12..18a83135b 100644
--- a/src/QMCWaveFunctions/EinsplineSetBuilder.h
+++ b/src/QMCWaveFunctions/EinsplineSetBuilder.h
@@ -80,6 +80,9 @@ namespace qmcplusplus {
   struct BandInfo {
     int TwistIndex, BandIndex, Spin;
     double Energy;
+    // This is true if we should make distinct copies 
+    // represeninting a +k, -k pair
+    bool MakeTwoCopies;
     inline bool operator<(BandInfo other) const
     { return Energy < other.Energy; }
   };
@@ -142,6 +145,7 @@ namespace qmcplusplus {
     // clone 
     vector<TinyVector<int,OHMMS_DIM> > UseTwists;
     vector<int> IncludeTwists, DistinctTwists;
+    int NumDistinctOrbitals;
     // This is true if the corresponding twist in DistinctTwists should
     // should be used to generate two distinct orbitals from the real and 
     // imaginary parts.