diff --git a/src/BOSampleStepper.C b/src/BOSampleStepper.C
index ff24986..d79129c 100644
--- a/src/BOSampleStepper.C
+++ b/src/BOSampleStepper.C
@@ -740,8 +740,10 @@ void BOSampleStepper::step(int niter)
         mixer.restart();
 
       double ehart, ehart_m;
+      bool scf_converged = false;
+      int itscf = 0;
 
-      for ( int itscf = 0; itscf < nitscf_; itscf++ )
+      while ( !scf_converged && itscf < nitscf_ )
       {
         if ( nite_ > 0 && onpe0 )
           cout << "  BOSampleStepper: start scf iteration" << endl;
@@ -852,13 +854,17 @@ void BOSampleStepper::step(int niter)
         // if ( onpe0 && nite_ > 0 )
         //   cout << " delta_ehart = " << delta_ehart << endl;
         int ite = 0;
-        double etotal_int, etotal_int_m;
-        double eigenvalue_sum, eigenvalue_sum_m;
+        double etotal_int, etotal_int_m = 0.0;
+        double etotal, etotal_m = 0.0;
+
+        double eigenvalue_sum, eigenvalue_sum_m = 0.0;
         // if nite == 0: do 1 iteration, no screening in charge mixing
         // if nite > 0: do nite iterations, use screening in charge mixing
+        //
+        double energy;
         while ( !nonscf_converged && ite < max(nite_,1) )
         {
-          double energy = ef_.energy(true,dwf,false,fion,false,sigma_eks);
+          energy = ef_.energy(true,dwf,false,fion,false,sigma_eks);
           double enthalpy = energy;
 
           if ( ite > 0 )
@@ -930,6 +936,10 @@ void BOSampleStepper::step(int niter)
           ite++;
         }
 
+        if ( itscf > 0 )
+          etotal_m = etotal;
+        etotal = energy;
+
         // if ( onpe0 && nite_ > 0 && ite >= nite_ )
         //   cout << " BOSampleStepper::step: nscf loop not converged after "
         //        << nite_ << " iterations" << endl;
@@ -987,7 +997,11 @@ void BOSampleStepper::step(int niter)
 
         if ( nite_ > 0 && onpe0 )
           cout << "  BOSampleStepper: end scf iteration" << endl;
-      } // for itscf
+
+        double delta_etotal = fabs(etotal - etotal_m);
+        scf_converged |= (delta_etotal < s_.ctrl.scf_tol);
+        itscf++;
+      } // while scf
 
       if ( compute_mlwf || compute_mlwfc )
       {
diff --git a/src/Control.h b/src/Control.h
index e4a5dc4..aa6edec 100644
--- a/src/Control.h
+++ b/src/Control.h
@@ -66,5 +66,7 @@ struct Control
 
   int blHF[3];
   double btHF;
+
+  double scf_tol;
 };
 #endif
diff --git a/src/Makefile b/src/Makefile
index 4c4ef3f..7c815c9 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -533,7 +533,8 @@ qb.o: SlaterDet.h Basis.h Matrix.h AtomsDyn.h BlHF.h BtHF.h Cell.h CellDyn.h
 qb.o: CellLock.h CellMass.h ChargeMixCoeff.h ChargeMixNdim.h ChargeMixRcut.h
 qb.o: Debug.h Dspin.h Ecut.h Ecutprec.h Ecuts.h Emass.h ExtStress.h
 qb.o: FermiTemp.h Dt.h Nempty.h NetCharge.h Nrowmax.h Nspin.h RefCell.h
-qb.o: Stress.h Thermostat.h ThTemp.h ThTime.h ThWidth.h WfDiag.h WfDyn.h Xc.h
+qb.o: ScfTol.h Stress.h Thermostat.h ThTemp.h ThTime.h ThWidth.h WfDiag.h
+qb.o: WfDyn.h Xc.h
 qbox_xmlns.o: qbox_xmlns.h
 QuitCmd.o: UserInterface.h Sample.h AtomSet.h Context.h Atom.h D3vector.h
 QuitCmd.o: UnitCell.h ConstraintSet.h ExtForceSet.h Wavefunction.h Control.h
@@ -588,6 +589,8 @@ SaveCmd.o: D3vector.h UnitCell.h ConstraintSet.h ExtForceSet.h Wavefunction.h
 SaveCmd.o: Control.h SampleWriter.h isodate.h release.h
 SaveCmd.o: UserInterface.h Sample.h AtomSet.h Context.h Atom.h D3vector.h
 SaveCmd.o: UnitCell.h ConstraintSet.h ExtForceSet.h Wavefunction.h Control.h
+ScfTol.o: Sample.h AtomSet.h Context.h Atom.h D3vector.h UnitCell.h
+ScfTol.o: ConstraintSet.h ExtForceSet.h Wavefunction.h Control.h
 SDAIonicStepper.o: SDAIonicStepper.h IonicStepper.h Sample.h AtomSet.h
 SDAIonicStepper.o: Context.h Atom.h D3vector.h UnitCell.h ConstraintSet.h
 SDAIonicStepper.o: ExtForceSet.h Wavefunction.h Control.h Species.h
@@ -659,8 +662,8 @@ testEnergyFunctional.o: UnitCell.h ConstraintSet.h ExtForceSet.h
 testEnergyFunctional.o: Wavefunction.h Control.h EnergyFunctional.h
 testEnergyFunctional.o: StructureFactor.h Timer.h
 test_fftw.o: Timer.h readTSC.h
-testFourierTransform.o: Context.h Basis.h D3vector.h UnitCell.h
-testFourierTransform.o: FourierTransform.h Timer.h
+testFourierTransform.o: Basis.h D3vector.h UnitCell.h FourierTransform.h
+testFourierTransform.o: Timer.h
 testjacobi.o: Timer.h Context.h Matrix.h jacobi.h
 testjade.o: Timer.h Context.h Matrix.h jade.h
 testLDAFunctional.o: LDAFunctional.h XCFunctional.h
diff --git a/src/ScfTol.h b/src/ScfTol.h
new file mode 100644
index 0000000..520e7c2
--- /dev/null
+++ b/src/ScfTol.h
@@ -0,0 +1,69 @@
+////////////////////////////////////////////////////////////////////////////////
+//
+// Copyright (c) 2014 The Regents of the University of California
+//
+// This file is part of Qbox
+//
+// Qbox is distributed under the terms of the GNU General Public License
+// as published by the Free Software Foundation, either version 2 of
+// the License, or (at your option) any later version.
+// See the file COPYING in the root directory of this distribution
+// or <http://www.gnu.org/licenses/>.
+//
+////////////////////////////////////////////////////////////////////////////////
+//
+// ScfTol.h
+//
+////////////////////////////////////////////////////////////////////////////////
+#ifndef SCFTOL_H
+#define SCFTOL_H
+
+#include<iostream>
+#include<iomanip>
+#include<sstream>
+#include<stdlib.h>
+
+#include "Sample.h"
+
+class ScfTol : public Var
+{
+  Sample *s;
+
+  public:
+
+  const char *name ( void ) const { return "scf_tol"; };
+
+  int set ( int argc, char **argv )
+  {
+    if ( argc != 2 )
+    {
+      if ( ui->onpe0() )
+      cout << " scf_tol takes only one value" << endl;
+      return 1;
+    }
+
+    double v = atof(argv[1]);
+    if ( v <= 0.0 )
+    {
+      if ( ui->onpe0() )
+        cout << " scf_tol must be non-negative" << endl;
+      return 1;
+    }
+
+    s->ctrl.scf_tol = v;
+    return 0;
+  }
+
+  string print (void) const
+  {
+     ostringstream st;
+     st.setf(ios::left,ios::adjustfield);
+     st << setw(10) << name() << " = ";
+     st.setf(ios::right,ios::adjustfield);
+     st << setw(10) << s->ctrl.scf_tol;
+     return st.str();
+  }
+
+  ScfTol(Sample *sample) : s(sample) { s->ctrl.scf_tol = 0.0; }
+};
+#endif
diff --git a/src/qb.C b/src/qb.C
index 970f260..62fdba6 100644
--- a/src/qb.C
+++ b/src/qb.C
@@ -101,6 +101,7 @@ using namespace std;
 #include "Nrowmax.h"
 #include "Nspin.h"
 #include "RefCell.h"
+#include "ScfTol.h"
 #include "Stress.h"
 #include "Thermostat.h"
 #include "ThTemp.h"
@@ -310,6 +311,7 @@ int main(int argc, char **argv, char **envp)
   ui.addVar(new Nspin(s));
   ui.addVar(new Dspin(s));
   ui.addVar(new RefCell(s));
+  ui.addVar(new ScfTol(s));
   ui.addVar(new Stress(s));
   ui.addVar(new Thermostat(s));
   ui.addVar(new ThTemp(s));