diff --git a/Nwpw/nwpwlib/D3dB/Pneb.cpp b/Nwpw/nwpwlib/D3dB/Pneb.cpp
index d3561ee7..e67515a7 100644
--- a/Nwpw/nwpwlib/D3dB/Pneb.cpp
+++ b/Nwpw/nwpwlib/D3dB/Pneb.cpp
@@ -977,12 +977,12 @@ void Pneb::ggm_sym_Multiply(double *psi1, double *psi2, double *hml)
       }
       d3db::parall->Vector_SumAll(1,ne[0]*ne[0]+ne[1]*ne[1],hml);
    }
-   if (d3db::parall->is_master())
-      std::cout << "hml= [" << Efmt(15,10) 
-                            << hml[0] << " " << hml[4] << " " << hml[8]  << " " << hml[12]  << std::endl 
-                << "      " << hml[1] << " " << hml[5] << " " << hml[9]  << " " << hml[13]  << std::endl
-                << "      " << hml[2] << " " << hml[6] << " " << hml[10] << " " << hml[14]  << std::endl
-                << "      " << hml[3] << " " << hml[7] << " " << hml[11] << " " << hml[15]  << "]"<< std::endl << std::endl;
+   //if (d3db::parall->is_master())
+   //   std::cout << "hml= [" << Efmt(15,10) 
+   //                         << hml[0] << " " << hml[4] << " " << hml[8]  << " " << hml[12]  << std::endl 
+   //             << "      " << hml[1] << " " << hml[5] << " " << hml[9]  << " " << hml[13]  << std::endl
+   //             << "      " << hml[2] << " " << hml[6] << " " << hml[10] << " " << hml[14]  << std::endl
+   //             << "      " << hml[3] << " " << hml[7] << " " << hml[11] << " " << hml[15]  << "]"<< std::endl << std::endl;
 }
 
 
@@ -2784,7 +2784,6 @@ void Pneb::g_project_out_filled_above(double *psi, const int ms, const int k, do
      int indx = 2*PGrid::npack(1)*ka + ishift;
      double w = -PGrid::cc_pack_dot(1,psi+indx,Horb);
      PGrid::cc_pack_daxpy(1,w,psi+indx,Horb);
-     std::cout << "project_above k=" << k << " ka=" << ka << " w=" << w << std::endl;
    }
 }
 
diff --git a/Nwpw/pspw/lib/molecule/Molecule.cpp b/Nwpw/pspw/lib/molecule/Molecule.cpp
index 56452761..916ed8fd 100644
--- a/Nwpw/pspw/lib/molecule/Molecule.cpp
+++ b/Nwpw/pspw/lib/molecule/Molecule.cpp
@@ -315,12 +315,10 @@ double Molecule::psi_KS_update_orb(const int ms, const int k, const int maxit_or
       e0 = mygrid->cc_pack_dot(1,orb,g);
 
       e0 = -e0;
-      std::cout << "   -- sub it =" << it << " e0=" << e0 << std::endl;
 
 
       double percent_error = 0.0;
       if(error0>1.0e-11) percent_error = std::abs(e0-eold)/error0;
-      std::cout << "   -- percent_error=" << percent_error << std::endl;
 
       precondition = (std::abs(e0-eold)>(sp*maxerror));
 
@@ -329,7 +327,6 @@ double Molecule::psi_KS_update_orb(const int ms, const int k, const int maxit_or
       mygrid->cc_pack_copy(1,g,r1);
       mygrid->cc_pack_daxpy(1,(e0),orb,r1);
 
-      std::cout << "   -- precondition=" << precondition << std::endl;
       //preconditioning 
       if (precondition)
       {
@@ -339,7 +336,6 @@ double Molecule::psi_KS_update_orb(const int ms, const int k, const int maxit_or
 
       //determine conjuagate direction ***
       double lmbda_r1 = mygrid->cc_pack_dot(1,r1,r1);
-      std::cout << "   -- lmbda_r1=" << lmbda_r1 << " lmbda_r0=" << lmbda_r0 << " lmbda_r1/lmbda_r0=" << lmbda_r1/lmbda_r0 << std::endl;
 
       mygrid->cc_pack_copy(1,r1,t);
 
@@ -423,9 +419,6 @@ double Molecule::psi_KS_update(const int maxit_orb, const double maxerror,
       int ishift = ms*neq[0]*2*mygrid->PGrid::npack(1);
       for (auto i=neq[ms]-1; i>=0; --i)
       {
-         double total_energy0 = this->energy();
-         std::cout << "i0=" << i << " TOTAL ENERGY=" << total_energy0 << std::endl;
-
          int indx = 2*mygrid->PGrid::npack(1)*i + ishift;
          double *orb = psi + indx;
 
@@ -433,25 +426,14 @@ double Molecule::psi_KS_update(const int maxit_orb, const double maxerror,
          //mygrid->g_project_out_filled_below(psi1, ms, i, orb);
          mygrid->g_project_out_filled_above(psi1, ms, i, orb);
 
-         double total_energy1 = this->energy();
-         std::cout << "i1=" << i << " TOTAL ENERGY=" << total_energy1 << std::endl;
-
          // normalize
          double norm = mygrid->cc_pack_dot(1,orb,orb);
          norm = 1.0/std::sqrt(norm);
          mygrid->c_pack_SMul(1,norm,orb);
 
-         double total_energy2 = this->energy();
-         std::cout << "i2=" << i << " TOTAL ENERGY=" << total_energy2 << std::endl;
-
          double e0 = psi_KS_update_orb(ms, i, maxit_orb, maxerror, perror, vall, orb,
                                        error_out, coutput);
          esum += e0;
-
-         std::cout << "i=" << i << " e0="<< e0<< std::endl;
-
-         double total_energy3 = this->energy();
-         std::cout << "i3=" << i << " TOTAL ENERGY=" << total_energy3 << std::endl;
       }
    }
 
diff --git a/Nwpw/pspw/lib/molecule/Molecule.hpp b/Nwpw/pspw/lib/molecule/Molecule.hpp
index 7d054f74..d1098059 100644
--- a/Nwpw/pspw/lib/molecule/Molecule.hpp
+++ b/Nwpw/pspw/lib/molecule/Molecule.hpp
@@ -174,7 +174,6 @@ class Molecule {
    {
       myelectron->run(psi1, rho1, dng1, rho1_all);
       myelectron->gen_energies_en(psi1, rho1, dng1, rho1_all, E, en);
-      std::cout << "A E[0] =" << E[0] << std::endl;
       
       /*  ion-ion energy */
       if (myelectron->is_periodic())
@@ -183,7 +182,6 @@ class Molecule {
         E[4] = myion->ion_ion_energy();
       
       E[0] += E[4];
-      std::cout << "B E[0] =" << E[0] << std::endl;
 
       /* get contraints energies */
       if (myion->has_ion_bond_constraints()) {