diff --git a/Nwpw/band/lib/cpsi/cpsi.cpp b/Nwpw/band/lib/cpsi/cpsi.cpp
index 0108f7ce..e372b5a2 100644
--- a/Nwpw/band/lib/cpsi/cpsi.cpp
+++ b/Nwpw/band/lib/cpsi/cpsi.cpp
@@ -17,31 +17,28 @@ namespace pwdft {
 
 /*****************************************************
  *                                                   *
- *               wvfnc_expander_convert              *
+ *              cwvfnc_expander_convert              *
  *                                                   *
  *****************************************************/
 
 static void cwvfnc_expander_convert(int ngrid[], double *psi1, int dngrid[], double *psi2) 
 {
    int indx, dindx, i2, j2, k2;
-   int nfft3d = (ngrid[0]) * ngrid[1] * ngrid[2];
-   int dnfft3d = (dngrid[0]) * dngrid[1] * dngrid[2];
-   int n2ft3d = 2 * nfft3d;
-   int dn2ft3d = 2 * dnfft3d;
-   int inc2 = (ngrid[0]);
-   int dinc2 = (dngrid[0]);
-   int inc3 = (ngrid[0]) * ngrid[1];
-   int dinc3 = (dngrid[0]) * dngrid[1];
+   int nfft3d  =  ngrid[0] *  ngrid[1] *  ngrid[2];
+   int dnfft3d = dngrid[0] * dngrid[1] * dngrid[2];
+   int n2ft3d  = 2*nfft3d;
+   int dn2ft3d = 2*dnfft3d;
+   int inc2  = ngrid[0];
+   int dinc2 = dngrid[0];
+   int inc3  = ngrid[0]*ngrid[1];
+   int dinc3 = dngrid[0]*dngrid[1];
  
    int n1 = ngrid[0];
    int n2 = ngrid[1];
    int n3 = ngrid[2];
-   if (n1 > dngrid[0])
-     n1 = dngrid[0];
-   if (n2 > dngrid[1])
-     n2 = dngrid[1];
-   if (n3 > dngrid[2])
-     n3 = dngrid[2];
+   if (n1 > dngrid[0]) n1 = dngrid[0];
+   if (n2 > dngrid[1]) n2 = dngrid[1];
+   if (n3 > dngrid[2]) n3 = dngrid[2];
  
    int idiff = dngrid[0] - ngrid[0];
    int jdiff = dngrid[1] - ngrid[1];
@@ -49,60 +46,68 @@ static void cwvfnc_expander_convert(int ngrid[], double *psi1, int dngrid[], dou
    bool ireverse = (idiff < 0);
    bool jreverse = (jdiff < 0);
    bool kreverse = (kdiff < 0);
-   if (jreverse)
-     jdiff = -jdiff;
-   if (kreverse)
-     kdiff = -kdiff;
- 
-   std::memset(psi2, 0, dn2ft3d * sizeof(double));
-   for (auto k = 0; k < n3; ++k)
-     for (auto j = 0; j < n2; ++j)
-       for (auto i = 0; i < (n1); ++i) {
-         indx = i;
-         dindx = i;
- 
-         if (k < (n3 / 2))
-           k2 = k;
-         else
-           k2 = kdiff + k;
- 
-         if (j < (n2 / 2))
-           j2 = j;
-         else
-           j2 = jdiff + j;
-
-         if (i < (n1 / 2))
-           i2 = i;
-         else
-           i2 = idiff + i;
-
-         if (ireverse) {
-           indx = indx + i2;
-           dindx = dindx + i;
-         } else {
-           indx = indx + i;
-           dindx = dindx + i2;
-         }
- 
-         if (jreverse) {
-           indx = indx + j2 * inc2;
-           dindx = dindx + j * dinc2;
-         } else {
-           indx = indx + j * inc2;
-           dindx = dindx + j2 * dinc2;
-         }
- 
-         if (kreverse) {
-           indx = indx + k2 * inc3;
-           dindx = dindx + k * dinc3;
-         } else {
-           indx = indx + k * inc3;
-           dindx = dindx + k2 * dinc3;
-         }
- 
-         psi2[2 * dindx] = psi1[2 * indx];
-         psi2[2 * dindx + 1] = psi1[2 * indx + 1];
-       }
+   if (jreverse) jdiff = -jdiff;
+   if (kreverse) kdiff = -kdiff;
+ 
+   std::memset(psi2, 0, dn2ft3d*sizeof(double));
+   for (auto k=0; k<n3; ++k)
+   for (auto j=0; j<n2; ++j)
+   for (auto i=0; i<n1; ++i) 
+   {
+      indx  = 0;
+      dindx = 0;
+ 
+      if (k < (n3/2))
+        k2 = k;
+      else
+        k2 = kdiff + k;
+ 
+      if (j < (n2/2))
+        j2 = j;
+      else
+        j2 = jdiff + j;
+
+      if (i < (n1/2))
+        i2 = i;
+      else
+        i2 = idiff + i;
+
+      if (ireverse) 
+      {
+         indx = indx   + i2;
+         dindx = dindx + i;
+      } 
+      else 
+      {
+         indx = indx   + i;
+         dindx = dindx + i2;
+      }
+ 
+      if (jreverse) 
+      {
+         indx = indx   + j2 * inc2;
+         dindx = dindx + j  * dinc2;
+      } 
+      else 
+      {
+         indx = indx   + j  * inc2;
+         dindx = dindx + j2 * dinc2;
+      }
+ 
+      if (kreverse) 
+      {
+         indx = indx   + k2 * inc3;
+         dindx = dindx + k  * dinc3;
+      } 
+      else 
+      {
+         indx = indx   + k * inc3;
+         dindx = dindx + k2 * dinc3;
+      }
+ 
+      psi2[2*dindx]   = psi1[2*indx];
+      psi2[2*dindx+1] = psi1[2*indx + 1];
+   }
 }
 
 /*****************************************************
@@ -154,14 +159,15 @@ static void cwvfnc_expander(Cneb *mycneb, char *filename, std::ostream &coutput)
      iwrite(6, &nbrillouin, 1);
      iwrite(6, &occupation, 1);
  
-     int n2ft3d = (nfft[0] + 2) * nfft[1] * nfft[2];
-     int dn2ft3d = (dnfft[0] + 2) * dnfft[1] * dnfft[2];
+     int n2ft3d  = 2* nfft[0] *  nfft[1] *  nfft[2];
+     int dn2ft3d = 2*dnfft[0] * dnfft[1] * dnfft[2];
      double *psi1 = new double[n2ft3d];
      double *psi2 = new double[dn2ft3d];
+     for (auto nb=0; nb<nbrillouin; ++nb)
      for (auto ms=0; ms<ispin; ++ms)
      for (auto n=0; n<ne[ms]; ++n) 
      {
-        if (lprint) coutput << " converting .... psi:" << n + 1 << " spin:" << ms + 1 << std::endl;
+        if (lprint) coutput << " converting .... cpsi: nb=" << nb+1 << " n=" << n+1 << " spin=" << ms + 1 << std::endl;
         dread(4, psi1, n2ft3d);
         cwvfnc_expander_convert(nfft, psi1, dnfft, psi2);
         dwrite(6, psi2, dn2ft3d);
diff --git a/Nwpw/nwpw.cpp b/Nwpw/nwpw.cpp
index adbb44f6..f1de6f37 100644
--- a/Nwpw/nwpw.cpp
+++ b/Nwpw/nwpw.cpp
@@ -863,7 +863,8 @@ int main(int argc, char *argv[]) {
      std::cout << "First task=" << task << std::endl << std::endl;
 
   // Initialize wavefunction
-  if ((task<10) && (parse_initialize_wvfnc(rtdbstr, true))) 
+  //if ((task<10) && (parse_initialize_wvfnc(rtdbstr, true))) 
+  if (parse_initialize_wvfnc(rtdbstr, true)) 
   {
      bool wvfnc_initialize = parse_initialize_wvfnc(rtdbstr, false);
      if (wvfnc_initialize)
@@ -891,12 +892,24 @@ int main(int argc, char *argv[]) {
               dum_rtdbstr = parse_initialize_wvfnc_set(dum_rtdbstr, true);
               wvfnc_initialize = false;
            }
-           if (oprint)
-              std::cout << std::endl
-                  << "Running staged energy optimization - lowlevel_rtdbstr = "
-                  << dum_rtdbstr << std::endl
-                  << std::endl;
-           ierr += pspw_minimizer(MPI_COMM_WORLD, dum_rtdbstr, std::cout);
+           if (task<10)
+           {
+              if (oprint)
+                 std::cout << std::endl
+                           << "Running staged pspw energy optimization - lowlevel_rtdbstr = "
+                           << dum_rtdbstr << std::endl
+                           << std::endl;
+              ierr += pwdft::pspw_minimizer(MPI_COMM_WORLD, dum_rtdbstr, std::cout);
+           }
+           else
+           {
+              if (oprint)
+                 std::cout << std::endl
+                           << "Running staged band energy optimization - lowlevel_rtdbstr = "
+                           << dum_rtdbstr << std::endl
+                           << std::endl;
+              ierr += pwdft::band_minimizer(MPI_COMM_WORLD, dum_rtdbstr, std::cout);
+           }
         }
      }
   }
diff --git a/Nwpw/pspw/lib/psi/psi.cpp b/Nwpw/pspw/lib/psi/psi.cpp
index 1391c7fe..0e0d9e94 100644
--- a/Nwpw/pspw/lib/psi/psi.cpp
+++ b/Nwpw/pspw/lib/psi/psi.cpp
@@ -29,73 +29,80 @@ namespace pwdft {
  *                                                   *
  *****************************************************/
 
-static void wvfnc_expander_convert(int ngrid[], double *psi1, int dngrid[],
-                                   double *psi2) {
-  int indx, dindx, j2, k2;
-  int nfft3d = (ngrid[0] / 2 + 1) * ngrid[1] * ngrid[2];
-  int dnfft3d = (dngrid[0] / 2 + 1) * dngrid[1] * dngrid[2];
-  int n2ft3d = 2 * nfft3d;
-  int dn2ft3d = 2 * dnfft3d;
-  int inc2 = (ngrid[0] / 2 + 1);
-  int dinc2 = (dngrid[0] / 2 + 1);
-  int inc3 = (ngrid[0] / 2 + 1) * ngrid[1];
-  int dinc3 = (dngrid[0] / 2 + 1) * dngrid[1];
-
-  int n1 = ngrid[0];
-  int n2 = ngrid[1];
-  int n3 = ngrid[2];
-  if (n1 > dngrid[0])
-    n1 = dngrid[0];
-  if (n2 > dngrid[1])
-    n2 = dngrid[1];
-  if (n3 > dngrid[2])
-    n3 = dngrid[2];
-
-  int jdiff = dngrid[1] - ngrid[1];
-  int kdiff = dngrid[2] - ngrid[2];
-  bool jreverse = (jdiff < 0);
-  bool kreverse = (kdiff < 0);
-  if (jreverse)
-    jdiff = -jdiff;
-  if (kreverse)
-    kdiff = -kdiff;
-
-  std::memset(psi2, 0, dn2ft3d * sizeof(double));
-  for (auto k = 0; k < n3; ++k)
-    for (auto j = 0; j < n2; ++j)
-      for (auto i = 0; i < (n1 / 2 + 1); ++i) {
-        indx = i;
-        dindx = i;
-
-        if (k < (n3 / 2))
-          k2 = k;
-        else
-          k2 = kdiff + k;
-
-        if (j < (n2 / 2))
-          j2 = j;
-        else
-          j2 = jdiff + j;
-
-        if (jreverse) {
-          indx = indx + j2 * inc2;
-          dindx = dindx + j * dinc2;
-        } else {
-          indx = indx + j * inc2;
-          dindx = dindx + j2 * dinc2;
-        }
-
-        if (kreverse) {
-          indx = indx + k2 * inc3;
-          dindx = dindx + k * dinc3;
-        } else {
-          indx = indx + k * inc3;
-          dindx = dindx + k2 * dinc3;
-        }
-
-        psi2[2 * dindx] = psi1[2 * indx];
-        psi2[2 * dindx + 1] = psi1[2 * indx + 1];
+static void wvfnc_expander_convert(int ngrid[], double *psi1, int dngrid[], double *psi2) 
+{
+   int indx, dindx, j2, k2;
+   int nfft3d = (ngrid[0] / 2 + 1) * ngrid[1] * ngrid[2];
+   int dnfft3d = (dngrid[0] / 2 + 1) * dngrid[1] * dngrid[2];
+   int n2ft3d = 2 * nfft3d;
+   int dn2ft3d = 2 * dnfft3d;
+   int inc2 = (ngrid[0] / 2 + 1);
+   int dinc2 = (dngrid[0] / 2 + 1);
+   int inc3 = (ngrid[0] / 2 + 1) * ngrid[1];
+   int dinc3 = (dngrid[0] / 2 + 1) * dngrid[1];
+ 
+   int n1 = ngrid[0];
+   int n2 = ngrid[1];
+   int n3 = ngrid[2];
+   if (n1 > dngrid[0])
+      n1 = dngrid[0];
+   if (n2 > dngrid[1])
+      n2 = dngrid[1];
+   if (n3 > dngrid[2])
+      n3 = dngrid[2];
+ 
+   int jdiff = dngrid[1] - ngrid[1];
+   int kdiff = dngrid[2] - ngrid[2];
+   bool jreverse = (jdiff < 0);
+   bool kreverse = (kdiff < 0);
+   if (jreverse)
+      jdiff = -jdiff;
+   if (kreverse)
+      kdiff = -kdiff;
+
+   std::memset(psi2, 0, dn2ft3d * sizeof(double));
+   for (auto k = 0; k < n3; ++k)
+   for (auto j = 0; j < n2; ++j)
+   for (auto i = 0; i < (n1/2+1); ++i) 
+   {
+      indx = i;
+      dindx = i;
+
+      if (k < (n3/2))
+         k2 = k;
+      else
+         k2 = kdiff + k;
+
+      if (j < (n2/2))
+         j2 = j;
+      else
+         j2 = jdiff + j;
+
+      if (jreverse) 
+      {
+         indx = indx   + j2*inc2;
+         dindx = dindx + j*dinc2;
+      } 
+      else 
+      {
+         indx = indx   + j*inc2;
+         dindx = dindx + j2*dinc2;
       }
+
+      if (kreverse) 
+      {
+         indx = indx   + k2*inc3;
+         dindx = dindx + k*dinc3;
+      } 
+      else 
+      {
+         indx = indx   + k*inc3;
+         dindx = dindx + k2*dinc3;
+      }
+
+      psi2[2*dindx]   = psi1[2*indx];
+      psi2[2*dindx+1] = psi1[2*indx+1];
+   }
 }
 
 /*****************************************************
@@ -115,67 +122,66 @@ static void wvfnc_expander(Pneb *mypneb, char *filename, std::ostream &coutput)
    strcpy(tmpfilename, filename);
    strcat(tmpfilename, ".wvfnc_expander");
  
-   if (myparall->is_master()) {
-     openfile(4, filename, "r");
-     iread(4, &version, 1);
-     iread(4, nfft, 3);
-     dread(4, unita, 9);
-     iread(4, &ispin, 1);
-     iread(4, ne, 2);
-     iread(4, &occupation, 1);
- 
-     dnfft[0] = mypneb->nx;
-     dnfft[1] = mypneb->ny;
-     dnfft[2] = mypneb->nz;
-     dunita[0] = mypneb->lattice->unita1d(0);
-     dunita[1] = mypneb->lattice->unita1d(1);
-     dunita[2] = mypneb->lattice->unita1d(2);
-     dunita[3] = mypneb->lattice->unita1d(3);
-     dunita[4] = mypneb->lattice->unita1d(4);
-     dunita[5] = mypneb->lattice->unita1d(5);
-     dunita[6] = mypneb->lattice->unita1d(6);
-     dunita[7] = mypneb->lattice->unita1d(7);
-     dunita[8] = mypneb->lattice->unita1d(8);
- 
-     openfile(6, tmpfilename, "w");
-     iwrite(6, &version, 1);
-     iwrite(6, dnfft, 3);
-     dwrite(6, dunita, 9);
-     iwrite(6, &ispin, 1);
-     iwrite(6, ne, 2);
-     iwrite(6, &occupation, 1);
- 
-     int n2ft3d = (nfft[0] + 2) * nfft[1] * nfft[2];
-     int dn2ft3d = (dnfft[0] + 2) * dnfft[1] * dnfft[2];
-     double *psi1 = new double[n2ft3d];
-     double *psi2 = new double[dn2ft3d];
-     for (auto ms = 0; ms < ispin; ++ms)
-       for (auto n = 0; n < ne[ms]; ++n) {
-         if (lprint)
-           coutput << " converting .... psi:" << n + 1 << " spin:" << ms + 1
-                   << std::endl;
+   if (myparall->is_master())
+   {
+      openfile(4, filename, "r");
+      iread(4, &version, 1);
+      iread(4, nfft, 3);
+      dread(4, unita, 9);
+      iread(4, &ispin, 1);
+      iread(4, ne, 2);
+      iread(4, &occupation, 1);
+     
+      dnfft[0] = mypneb->nx;
+      dnfft[1] = mypneb->ny;
+      dnfft[2] = mypneb->nz;
+      dunita[0] = mypneb->lattice->unita1d(0);
+      dunita[1] = mypneb->lattice->unita1d(1);
+      dunita[2] = mypneb->lattice->unita1d(2);
+      dunita[3] = mypneb->lattice->unita1d(3);
+      dunita[4] = mypneb->lattice->unita1d(4);
+      dunita[5] = mypneb->lattice->unita1d(5);
+      dunita[6] = mypneb->lattice->unita1d(6);
+      dunita[7] = mypneb->lattice->unita1d(7);
+      dunita[8] = mypneb->lattice->unita1d(8);
+     
+      openfile(6, tmpfilename, "w");
+      iwrite(6, &version, 1);
+      iwrite(6, dnfft, 3);
+      dwrite(6, dunita, 9);
+      iwrite(6, &ispin, 1);
+      iwrite(6, ne, 2);
+      iwrite(6, &occupation, 1);
+     
+      int n2ft3d = (nfft[0] + 2) * nfft[1] * nfft[2];
+      int dn2ft3d = (dnfft[0] + 2) * dnfft[1] * dnfft[2];
+      double *psi1 = new double[n2ft3d];
+      double *psi2 = new double[dn2ft3d];
+      for (auto ms = 0; ms < ispin; ++ms)
+      for (auto n = 0;  n  < ne[ms]; ++n) 
+      {
+         if (lprint) coutput << " converting .... psi:" << n + 1 << " spin:" << ms + 1 << std::endl;
          dread(4, psi1, n2ft3d);
          wvfnc_expander_convert(nfft, psi1, dnfft, psi2);
          dwrite(6, psi2, dn2ft3d);
-       }
-     if (lprint)
-       coutput << std::endl;
-     if (occupation > 0) {
-       double *occ1 = new double[ne[0] + ne[1]];
-       dread(4, occ1, (ne[0] + ne[1]));
-       dwrite(6, occ1, (ne[0] + ne[1]));
-       delete[] occ1;
-     }
- 
-     delete[] psi1;
-     delete[] psi2;
- 
-     closefile(4);
-     closefile(6);
- 
-     /* copy and delete the tmpfilenane file */
-     std::rename(tmpfilename, filename);
-     // std::remove(tmpfilename);
+      }
+      if (lprint) coutput << std::endl;
+      if (occupation > 0) {
+        double *occ1 = new double[ne[0] + ne[1]];
+        dread(4, occ1, (ne[0] + ne[1]));
+        dwrite(6, occ1, (ne[0] + ne[1]));
+        delete[] occ1;
+      }
+     
+      delete[] psi1;
+      delete[] psi2;
+     
+      closefile(4);
+      closefile(6);
+     
+      /* copy and delete the tmpfilenane file */
+      std::rename(tmpfilename, filename);
+      // std::remove(tmpfilename);
    }
 }
 
@@ -186,22 +192,24 @@ static void wvfnc_expander(Pneb *mypneb, char *filename, std::ostream &coutput)
  *****************************************************/
 
 void psi_get_header(Parallel *myparall, int *version, int nfft[],
-                    double unita[], int *ispin, int ne[], char *filename) {
-  if (myparall->is_master()) {
-    // char *fname = control_input_movecs_filename();
-    openfile(4, filename, "r");
-    iread(4, version, 1);
-    iread(4, nfft, 3);
-    dread(4, unita, 9);
-    iread(4, ispin, 1);
-    iread(4, ne, 2);
-    closefile(4);
-  }
-  myparall->Brdcst_iValue(0, 0, version);
-  myparall->Brdcst_iValues(0, 0, 3, nfft);
-  myparall->Brdcst_Values(0, 0, 9, unita);
-  myparall->Brdcst_iValue(0, 0, ispin);
-  myparall->Brdcst_iValues(0, 0, 2, ne);
+                    double unita[], int *ispin, int ne[], char *filename) 
+{
+   if (myparall->is_master()) 
+   {
+      // char *fname = control_input_movecs_filename();
+      openfile(4, filename, "r");
+      iread(4, version, 1);
+      iread(4, nfft, 3);
+      dread(4, unita, 9);
+      iread(4, ispin, 1);
+      iread(4, ne, 2);
+      closefile(4);
+   }
+   myparall->Brdcst_iValue(0, 0, version);
+   myparall->Brdcst_iValues(0, 0, 3, nfft);
+   myparall->Brdcst_Values(0, 0, 9, unita);
+   myparall->Brdcst_iValue(0, 0, ispin);
+   myparall->Brdcst_iValues(0, 0, 2, ne);
 }
 
 /*****************************************************
@@ -209,25 +217,25 @@ void psi_get_header(Parallel *myparall, int *version, int nfft[],
  *                psi_check_convert                  *
  *                                                   *
  *****************************************************/
-static bool psi_check_convert(Pneb *mypneb, char *filename,
-                              std::ostream &coutput) {
-  Parallel *myparall = mypneb->d3db::parall;
-  int version0, ispin0, nfft0[3], ne0[2];
-  double unita0[9];
-  bool converted = false;
-
-  psi_get_header(myparall, &version0, nfft0, unita0, &ispin0, ne0, filename);
-  if ((nfft0[0] != mypneb->nx) || (nfft0[1] != mypneb->ny) || (nfft0[2] != mypneb->nz)) 
-  {
-     if (myparall->base_stdio_print)
-        coutput << " psi grids are being converted: " << std::endl
-                << " -----------------------------: " << std::endl;
-    
-     wvfnc_expander(mypneb, filename, coutput);
-     converted = true;
-  }
-
-  return converted;
+static bool psi_check_convert(Pneb *mypneb, char *filename, std::ostream &coutput) 
+{
+   Parallel *myparall = mypneb->d3db::parall;
+   int version0, ispin0, nfft0[3], ne0[2];
+   double unita0[9];
+   bool converted = false;
+ 
+   psi_get_header(myparall, &version0, nfft0, unita0, &ispin0, ne0, filename);
+   if ((nfft0[0] != mypneb->nx) || (nfft0[1] != mypneb->ny) || (nfft0[2] != mypneb->nz)) 
+   {
+      if (myparall->base_stdio_print)
+         coutput << " psi grids are being converted: " << std::endl
+                 << " -----------------------------: " << std::endl;
+     
+      wvfnc_expander(mypneb, filename, coutput);
+      converted = true;
+   }
+ 
+   return converted;
 }
 
 /*****************************************************
@@ -242,33 +250,34 @@ static bool psi_check_convert(Pneb *mypneb, char *filename,
 
 */
 void psi_read0(Pneb *mypneb, int *version, int nfft[], double unita[],
-               int *ispin, int ne[], double *psi, char *filename) {
-  int occupation;
-
-  Parallel *myparall = mypneb->d3db::parall;
-
-  if (myparall->is_master()) {
-    openfile(4, filename, "r");
-    iread(4, version, 1);
-    iread(4, nfft, 3);
-    dread(4, unita, 9);
-    iread(4, ispin, 1);
-    iread(4, ne, 2);
-    iread(4, &occupation, 1);
-  }
-  myparall->Brdcst_iValue(0, 0, version);
-  myparall->Brdcst_iValues(0, 0, 3, nfft);
-  myparall->Brdcst_Values(0, 0, 9, unita);
-  myparall->Brdcst_iValue(0, 0, ispin);
-  myparall->Brdcst_iValues(0, 0, 2, ne);
-
-  /* reads in c format and automatically packs the result to g format */
-  //mypneb->g_read(4,ispin,psi);
-  mypneb->g_read_ne(4,ne,psi);
-  
-
-  if (myparall->is_master())
-    closefile(4);
+               int *ispin, int ne[], double *psi, char *filename) 
+{
+   int occupation;
+ 
+   Parallel *myparall = mypneb->d3db::parall;
+ 
+   if (myparall->is_master()) {
+     openfile(4, filename, "r");
+     iread(4, version, 1);
+     iread(4, nfft, 3);
+     dread(4, unita, 9);
+     iread(4, ispin, 1);
+     iread(4, ne, 2);
+     iread(4, &occupation, 1);
+   }
+   myparall->Brdcst_iValue(0, 0, version);
+   myparall->Brdcst_iValues(0, 0, 3, nfft);
+   myparall->Brdcst_Values(0, 0, 9, unita);
+   myparall->Brdcst_iValue(0, 0, ispin);
+   myparall->Brdcst_iValues(0, 0, 2, ne);
+ 
+   /* reads in c format and automatically packs the result to g format */
+   //mypneb->g_read(4,ispin,psi);
+   mypneb->g_read_ne(4,ne,psi);
+   
+ 
+   if (myparall->is_master())
+     closefile(4);
 }
 
 /*****************************************************
@@ -386,17 +395,18 @@ void psi_write(Pneb *mypneb, int *version, int nfft[], double unita[],
  *                psi_filefind                       *
  *                                                   *
  *****************************************************/
-bool psi_filefind(Pneb *mypneb, char *filename) {
-  int ifound = 0;
-  Parallel *myparall = mypneb->d3db::parall;
-
-  if (myparall->is_master()) {
-    ifound = cfileexists(filename);
-  }
-
-  myparall->Brdcst_iValue(0, 0, &ifound);
-
-  return (ifound > 0);
+bool psi_filefind(Pneb *mypneb, char *filename) 
+{
+   int ifound = 0;
+   Parallel *myparall = mypneb->d3db::parall;
+ 
+   if (myparall->is_master()) {
+     ifound = cfileexists(filename);
+   }
+ 
+   myparall->Brdcst_iValue(0, 0, &ifound);
+ 
+   return (ifound > 0);
 }
 
 /*
diff --git a/QA/BAND-B/band9d.nw b/QA/BAND-B/band9d.nw
index 94207d20..a2fa1969 100644
--- a/QA/BAND-B/band9d.nw
+++ b/QA/BAND-B/band9d.nw
@@ -36,7 +36,7 @@ nwpw
     tolerances 1.0e-9 1.0e-9
   end
   loop 10 1000
-  tolerances 1.0e-8 1.0e-8
+  #tolerances 1.0e-8 1.0e-8
   #vectors input  band-pbe-B2-mylabel.movecs output junk.elc
   vectors input  junk.elc output junk.elc