diff --git a/Nwpw/pspw/lib/molecule/Molecule.cpp b/Nwpw/pspw/lib/molecule/Molecule.cpp
index 9089da11..55e1d7a2 100644
--- a/Nwpw/pspw/lib/molecule/Molecule.cpp
+++ b/Nwpw/pspw/lib/molecule/Molecule.cpp
@@ -144,7 +144,7 @@ void Molecule::psi_minimize(double *vall, std::ostream &coutput)
             bool continue_inner_loop = true;
             for (int l=0; continue_inner_loop && l<=(1+(l2-1)*3); ++l)
             {
-               eorb0 = psi_KS_update(ms,k,2,tole,0.001,vall,orb,&error_out,coutput);
+               eorb0 = psi_KS_update_orb(ms,k,2,tole,0.001,vall,orb,&error_out,coutput);
                if (error_out <= tole)
                   continue_inner_loop = false; // Exit the inner loop
             }
@@ -235,7 +235,7 @@ void Molecule::psi_get_gradient(double *orb, double *vall, double *Horb)
 
 /********************************************
  *                                          *
- *         Molecule::psi_KS_update          *
+ *         Molecule::psi_KS_update_orb      *
  *                                          *
  ********************************************/
 //    This routine performs a KS update on orbital i
@@ -283,7 +283,7 @@ void Molecule::psi_get_gradient(double *orb, double *vall, double *Horb)
  * - `mygrid::cc_pack_daxpy()`
  * - `myelectron::get_myke()->ke_precondition()`
  */
-double Molecule::psi_KS_update(const int ms, const int k, const int maxit_orb, const double maxerror,
+double Molecule::psi_KS_update_orb(const int ms, const int k, const int maxit_orb, const double maxerror,
                                         const double perror, double *vall, double *orb,
                                         double *error_out, std::ostream &coutput)
 {        
@@ -403,6 +403,40 @@ double Molecule::psi_KS_update(const int ms, const int k, const int maxit_orb, c
    return e0;
 }
 
+/********************************************
+ *                                          *
+ *         Molecule::psi_KS_update          *
+ *                                          *
+ ********************************************/
+double Molecule::psi_KS_update(const int maxit_orb, const double maxerror,
+                               const double perror, double *vall, 
+                               const int ispin, const int *neq, double *psi,
+                               double *error_out, std::ostream &coutput)
+
+{        
+
+   for (auto ms=0; ms<ispin; ++ms)
+   {
+      int ishift = ms*neq[0]*2*mygrid->PGrid::npack(1);
+      for (auto i=0; i<neq[ms]; ++i)
+      {
+         int indx = 2*mygrid->PGrid::npack(1)*i + ishift;
+         double *orb = psi + indx;
+
+         // orthogonalize to lower orbitals
+         mygrid->g_project_out_filled_below(psi1, ms, i, orb);
+
+         // normalize
+         double norm = mygrid->cc_pack_dot(1,orb,orb);
+         norm = 1.0/std::sqrt(norm);
+         mygrid->c_pack_SMul(1,norm,orb);
+
+         double e0 = psi_KS_update_orb(ms, i, maxit_orb, maxerror, perror, vall, orb,
+                                       error_out, coutput);
+      }
+   }
+}        
+
 /********************************************
  *                                          *
  *        psi_linesearch_update             *
diff --git a/Nwpw/pspw/lib/molecule/Molecule.hpp b/Nwpw/pspw/lib/molecule/Molecule.hpp
index 8abe3053..62976a84 100644
--- a/Nwpw/pspw/lib/molecule/Molecule.hpp
+++ b/Nwpw/pspw/lib/molecule/Molecule.hpp
@@ -121,7 +121,13 @@ class Molecule {
 
    void psi_minimize(double *, std::ostream &);
    void psi_get_gradient(double *, double *, double *);
-   double psi_KS_update(const int, const int, const int, const double, const double, double *, double *, double *, std::ostream &);
+   double psi_KS_update_orb(const int, const int, const int, const double, const double, double *, double *, double *, std::ostream &);
+
+
+   double psi_KS_update(const int, const double, const double, double *, 
+                        const int, const int *, double *,
+                        double *, std::ostream &);
+
 
    void psi_linesearch_update(double, double, double *, double *, double *, double *);
    void psi_sort();
diff --git a/Nwpw/pspw/minimizer/cgsd.hpp b/Nwpw/pspw/minimizer/cgsd.hpp
index 7731a519..b054063e 100644
--- a/Nwpw/pspw/minimizer/cgsd.hpp
+++ b/Nwpw/pspw/minimizer/cgsd.hpp
@@ -30,7 +30,7 @@ extern double cgsd_bybminimize(Molecule &, Geodesic *, double *, double *,
 extern double cgsd_bybminimize2(Molecule &, Geodesic *, double *, double *,
                               double *, int, int,int,
                               nwpw_scf_mixing &,
-                              double, double);
+                              double, double, std::ostream &);
 
 } // namespace pwdft
 #endif
diff --git a/Nwpw/pspw/minimizer/cgsd_bybminimize2.cpp b/Nwpw/pspw/minimizer/cgsd_bybminimize2.cpp
index 27d4d4ff..35f48b6c 100644
--- a/Nwpw/pspw/minimizer/cgsd_bybminimize2.cpp
+++ b/Nwpw/pspw/minimizer/cgsd_bybminimize2.cpp
@@ -31,7 +31,7 @@ double cgsd_bybminimize2(Molecule &mymolecule, Geodesic *mygeodesic, double *E,
                         double *deltae, double *deltac, int current_iteration,
                         int ks_it_in, int ks_it_out, 
                         nwpw_scf_mixing &scfmix,
-                        double tole, double tolc) 
+                        double tole, double tolc, std::ostream &coutput) 
 {
    bool done = false;
    double tmin = 0.0;
@@ -40,11 +40,19 @@ double cgsd_bybminimize2(Molecule &mymolecule, Geodesic *mygeodesic, double *E,
    double sum0, sum1, scale, total_energy;
    double dE, max_sigma, min_sigma;
    double Eold, dEold, Enew;
-   double tmin0, deltae0;
+   double tmin0, deltae0, perror,error_out;
  
    Pneb *mygrid = mymolecule.mygrid;
    mygeodesic_ptr = mygeodesic;
-   double *vall_out;
+
+   double ks_deltae = tole;
+   int ispin = mygrid->ispin;
+   int *neq = mygrid->neq;
+   double *vall = mygrid->r_nalloc(ispin);
+   mymolecule.gen_vall();
+   mymolecule.get_vall(vall);
+
+
 
    // ion-ion energy 
    double eion = mymolecule.eion();
@@ -89,6 +97,9 @@ double cgsd_bybminimize2(Molecule &mymolecule, Geodesic *mygeodesic, double *E,
    //                       scf_algorithm,scf_alpha,diis_histories,
    //                       mygrid->ispin,mygrid->n2ft3d,vall_out);
 
+   double e0 = mymolecule.psi_KS_update(ks_it_in,ks_deltae,perror,vall,
+                                        ispin, neq, mymolecule.psi1,
+                                        deltae, coutput);
 
  
    /* iniitialize blocked cg */
diff --git a/Nwpw/pspw/minimizer/cgsd_energy.cpp b/Nwpw/pspw/minimizer/cgsd_energy.cpp
index c61e5e04..f40cdc04 100644
--- a/Nwpw/pspw/minimizer/cgsd_energy.cpp
+++ b/Nwpw/pspw/minimizer/cgsd_energy.cpp
@@ -71,7 +71,7 @@ double cgsd_energy(Control2 &control, Molecule &mymolecule, bool doprint, std::o
    double kerker_g0 = control.kerker_g0();
    int diis_histories = control.diis_histories();
    int scf_algorithm = control.scf_algorithm();
- 
+
    double dt = control.time_step();
    double dte = dt/sqrt(control.fake_mass());
  
@@ -103,8 +103,8 @@ double cgsd_energy(Control2 &control, Molecule &mymolecule, bool doprint, std::o
    }
  
    // if (minimizer > 1) pspw_Grsm_list_start()
-   if ((minimizer == 5) || (minimizer == 8))
-      it_out = 1;
+   //if ((minimizer == 5) || (minimizer == 8))
+   //   it_out = 1;
  
    Geodesic12 mygeodesic12(minimizer, &mymolecule, control);
  
@@ -378,8 +378,7 @@ double cgsd_energy(Control2 &control, Molecule &mymolecule, bool doprint, std::o
                              scf_algorithm,scf_alpha,diis_histories,
                              mygrid->ispin,mygrid->n2ft3d,mymolecule.rho1);
 
-
-      while ((icount < it_out*it_in) && (!converged))
+      while ((icount < (it_out*it_in)) && (!converged))
       {
          ++icount;
          if (stalled)
@@ -407,7 +406,7 @@ double cgsd_energy(Control2 &control, Molecule &mymolecule, bool doprint, std::o
 
          total_energy = cgsd_bybminimize2(mymolecule,mygeodesic12.mygeodesic1,E,&deltae,
                                          &deltac,bfgscount,ks_it_in,ks_it_out,
-                                         scfmix,tole,tolc);
+                                         scfmix,tole,tolc,coutput);
         ++bfgscount;
         if (oprint) 
           coutput << Ifmt(10) << icount
@@ -416,9 +415,10 @@ double cgsd_energy(Control2 &control, Molecule &mymolecule, bool doprint, std::o
                   << Efmt(16,6) << deltac << std::endl;
         if ((std::fabs(deltae) > fabs(deltae_old)) ||
             (std::fabs(deltae) > 1.0e-2) || (deltae > 0.0))
-          stalled = true;
+           stalled = true;
         else
-          stalled = false;
+           stalled = false;
+
         converged = (std::fabs(deltae) < tole) && (deltac < tolc);
       }    
    } else if (minimizer == 9) {