diff --git a/src/BOSampleStepper.C b/src/BOSampleStepper.C
index da507c5..559541a 100644
--- a/src/BOSampleStepper.C
+++ b/src/BOSampleStepper.C
@@ -1012,7 +1012,8 @@ void BOSampleStepper::step(int niter)
 
         if ( fractional_occ )
         {
-          wf.update_occ(s_.ctrl.fermi_temp);
+          if ( s_.ctrl.fermi_temp > 0 )
+            wf.update_occ(s_.ctrl.fermi_temp);
 #if 0
           if ( onpe0 )
           {
diff --git a/src/Makefile b/src/Makefile
index a047193..ec499a3 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -548,6 +548,9 @@ Nrowmax.o: Control.h
 Nspin.o: Sample.h AtomSet.h Context.h blacs.h Atom.h D3vector.h UnitCell.h
 Nspin.o: D3tensor.h blas.h ConstraintSet.h ExtForceSet.h Wavefunction.h
 Nspin.o: Control.h
+Occ.o: Sample.h AtomSet.h Context.h blacs.h Atom.h D3vector.h UnitCell.h
+Occ.o: D3tensor.h blas.h ConstraintSet.h ExtForceSet.h Wavefunction.h
+Occ.o: Control.h
 PBEFunctional.o: PBEFunctional.h XCFunctional.h
 PBEFunctional.o: XCFunctional.h
 PSDAWavefunctionStepper.o: PSDAWavefunctionStepper.h WavefunctionStepper.h
@@ -837,8 +840,8 @@ qb.o: CellLock.h CellMass.h ChargeMixCoeff.h ChargeMixNdim.h ChargeMixRcut.h
 qb.o: Debug.h Dspin.h Ecut.h Ecutprec.h Ecuts.h Efield.h ForceTol.h
 qb.o: Polarization.h Emass.h ExtStress.h FermiTemp.h IterCmd.h
 qb.o: IterCmdPeriod.h Dt.h MuRSH.h Nempty.h NetCharge.h Nrowmax.h Nspin.h
-qb.o: RefCell.h ScfTol.h Stress.h StressTol.h Thermostat.h ThTemp.h ThTime.h
-qb.o: ThWidth.h Vext.h ExternalPotential.h WfDiag.h WfDyn.h Xc.h
+qb.o: Occ.h RefCell.h ScfTol.h Stress.h StressTol.h Thermostat.h ThTemp.h
+qb.o: ThTime.h ThWidth.h Vext.h ExternalPotential.h WfDiag.h WfDyn.h Xc.h
 qbox_xmlns.o: qbox_xmlns.h
 release.o: release.h
 sinft.o: sinft.h
diff --git a/src/Occ.h b/src/Occ.h
new file mode 100644
index 0000000..be873cf
--- /dev/null
+++ b/src/Occ.h
@@ -0,0 +1,179 @@
+////////////////////////////////////////////////////////////////////////////////
+//
+// Copyright (c) 2008 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/>.
+//
+////////////////////////////////////////////////////////////////////////////////
+//
+// Occ.h
+//
+////////////////////////////////////////////////////////////////////////////////
+
+#ifndef OCC_H
+#define OCC_H
+
+#include<iostream>
+#include<iomanip>
+#include<sstream>
+#include<stdlib.h>
+
+#include "Sample.h"
+
+class Occ : public Var
+{
+  Sample *s;
+
+  public:
+
+  const char *name ( void ) const { return "occ"; };
+
+  int set ( int argc, char **argv )
+  {
+    if ( (argc!=3) && (argc!=4) )
+    {
+      if ( ui->onpe0() )
+      {
+        cout << " use: set occ [ispin] n f" << endl;
+        cout << " ispin = {1|2}" << endl;
+      }
+      return 1;
+    }
+
+    if ( s->wf.nkp() != 1 )
+    {
+      if ( ui->onpe0() )
+        cout << " set occ: not implemented for multiple k-points" << endl;
+      return 1;
+    }
+
+    if ( argc == 3 )
+    {
+      // set occ n f
+      if ( s->wf.nspin() != 1 )
+      {
+        if ( ui->onpe0() )
+          cout << "nspin must be 1" << endl;
+        return 1;
+      }
+
+      const int n = atoi(argv[1]);
+      if ( (n < 1) || n > s->wf.nst() )
+      {
+        if ( ui->onpe0() )
+          cout << " n must be in [1," << s->wf.nst() << "]" << endl;
+        return 1;
+      }
+
+      const double f = atof(argv[2]);
+      if ( (f < 0.0) || (f > 2.0) )
+      {
+        if ( ui->onpe0() )
+          cout << " f must be in [0,2]" << endl;
+        return 1;
+      }
+
+      Wavefunction& wf = s->wf;
+      SlaterDet& sd = *wf.sd(0,0);
+      // n-1 in next line: states are numbered starting from 1 in
+      // the set occ command
+      sd.set_occ(n-1,f);
+      // recompute total electronic charge
+      double sum = 0.0;
+      for ( int i = 0; i < sd.nst(); i++ )
+        sum += sd.occ(i);
+      if ( ui->onpe0() )
+        cout << " total electronic charge: " << sum << endl;
+      // adjust total number of electrons
+      //wf.set_nel((int)sum);
+    }
+
+    if ( argc == 4 )
+    {
+      // set occ ispin n f
+      if ( s->wf.nspin() != 2 )
+      {
+        if ( ui->onpe0() )
+          cout << "nspin must be 2 when using ispin" << endl;
+        return 1;
+      }
+
+      const int ispin  = atoi(argv[1]);
+      if ( (ispin < 1) || (ispin > 2) )
+      {
+        if ( ui->onpe0() )
+          cout << " ispin must be 1 or 2" << endl;
+        return 1;
+      }
+
+      const int n = atoi(argv[2]);
+      if ( (n < 1) || n > s->wf.nst(ispin-1) )
+      {
+        if ( ui->onpe0() )
+          cout << " n must be in [1," << s->wf.nst(ispin-1) << "]" << endl;
+        return 1;
+      }
+
+      const double f = atof(argv[3]);
+      if ( (f < 0.0) || (f > 1.0) )
+      {
+        if ( ui->onpe0() )
+          cout << " f must be in [0,1] when using ispin" << endl;
+        return 1;
+      }
+
+      Wavefunction& wf = s->wf;
+      // ispin-1 in next line: spins are numbered starting from 1 in
+      // the set occ command
+      SlaterDet& sd = *wf.sd(ispin-1,0);
+      // n-1 in next line: states are numbered starting from 1 in
+      // the set occ command
+      sd.set_occ(n-1,f);
+      // recompute total electronic charge
+      double sum = 0.0;
+      for ( int isp = 0; isp < wf.nspin(); isp++ )
+        for ( int i = 0; i < wf.nst(isp); i++ )
+          sum += wf.sd(isp,0)->occ(i);
+      if ( ui->onpe0() )
+        cout << " total electronic charge: " << sum << endl;
+      // adjust total number of electrons
+      //wf.set_nel((int)sum);
+    }
+    return 0;
+  }
+
+  string print (void) const
+  {
+    ostringstream st;
+    st.setf(ios::right,ios::adjustfield);
+    st.setf(ios::fixed,ios::floatfield);
+
+    const Wavefunction& wf = s->wf;
+    st << " occupation numbers" << endl;
+    for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
+    {
+      const SlaterDet& sd = *wf.sd(ispin,0);
+      const int nst = sd.nst();
+      for ( int n = 0; n < nst; n++ )
+      {
+        st << setw(7) << setprecision(4) << sd.occ(n);
+        if ( ( n%10 ) == 9 ) st << endl;
+      }
+      if ( nst % 10 != 0 )
+        st << endl;
+      if ( (wf.nspin() == 2) && (ispin == 0) )
+        st << endl;
+    }
+    st << " total electronic charge: " << wf.nel();
+    return st.str();
+ }
+
+  Occ(Sample *sample) : s(sample) {};
+};
+#endif
diff --git a/src/SlaterDet.h b/src/SlaterDet.h
index 7d0e39e..aa2a581 100644
--- a/src/SlaterDet.h
+++ b/src/SlaterDet.h
@@ -89,6 +89,7 @@ class SlaterDet
   const double* occ_ptr(int i) const { return &occ_[i]; }
   void set_occ(std::vector<double>& occ)
     { assert(occ_.size()==occ.size()); occ_ = occ; }
+  void set_occ(int i, double f) { occ_[i] = f; }
   void set_eig(std::vector<double>& eig)
     { assert(eig_.size()==eig.size()); eig_ = eig; }
   void set_eig(std::valarray<double>& eig)
diff --git a/src/Wavefunction.C b/src/Wavefunction.C
index 47633c4..13b9c88 100644
--- a/src/Wavefunction.C
+++ b/src/Wavefunction.C
@@ -42,15 +42,12 @@ nempty_(0), nspin_(1), deltaspin_(0), ecut_(0.0), nrowmax_(32)
 
 ////////////////////////////////////////////////////////////////////////////////
 Wavefunction::Wavefunction(const Wavefunction& wf) : ctxt_(wf.ctxt_),
-nel_(wf.nel_), nempty_(wf.nempty_), nspin_(wf.nspin_),
+nel_(wf.nel_), nempty_(wf.nempty_), nspin_(wf.nspin_), nst_(wf.nst_),
 deltaspin_(wf.deltaspin_), nrowmax_(wf.nrowmax_),
 cell_(wf.cell_), refcell_(wf.refcell_),
 ecut_(wf.ecut_), weight_(wf.weight_), kpoint_(wf.kpoint_)
 {
   // Create a Wavefunction using the dimensions of the argument
-
-  compute_nst();
-
   // Next lines: do special allocation of contexts to ensure that
   // contexts are same as those of wf
   spincontext_ = new Context(*wf.spincontext_);
@@ -58,8 +55,7 @@ ecut_(wf.ecut_), weight_(wf.weight_), kpoint_(wf.kpoint_)
   sdcontext_ = new Context(*wf.sdcontext_);
 
   allocate();
-
-  resize(cell_,refcell_,ecut_);
+  resize();
   init();
 }
 
@@ -147,6 +143,34 @@ double Wavefunction::entropy(void) const
 }
 
 ////////////////////////////////////////////////////////////////////////////////
+void Wavefunction::resize(void)
+{
+  try
+  {
+    // resize all SlaterDets using current cell_, refcell_, ecut_, nst_[ispin]
+    for ( int ispin = 0; ispin < nspin_; ispin++ )
+    {
+      for ( int ikp = 0; ikp < kpoint_.size(); ikp++ )
+      {
+        sd_[ispin][ikp]->resize(cell_,refcell_,ecut_,nst_[ispin]);
+      }
+    }
+  }
+  catch ( const SlaterDetException& sdex )
+  {
+    cout << " Wavefunction: SlaterDetException during resize: " << endl
+         << sdex.msg << endl;
+    // no resize took place
+    return;
+  }
+  catch ( bad_alloc )
+  {
+    cout << " Wavefunction: insufficient memory for resize operation" << endl;
+    return;
+  }
+}
+
+////////////////////////////////////////////////////////////////////////////////
 void Wavefunction::resize(const UnitCell& cell, const UnitCell& refcell,
   double ecut)
 {
@@ -176,7 +200,6 @@ void Wavefunction::resize(const UnitCell& cell, const UnitCell& refcell,
     cout << " Wavefunction: insufficient memory for resize operation" << endl;
     return;
   }
-
 }
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -252,7 +275,7 @@ void Wavefunction::set_nel(int nel)
 
   nel_ = nel;
   compute_nst();
-  resize(cell_,refcell_,ecut_);
+  resize();
 }
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -267,7 +290,7 @@ void Wavefunction::set_nempty(int nempty)
   nempty_ = nempty;
   compute_nst();
   update_occ(0.0);
-  resize(cell_,refcell_,ecut_);
+  resize();
 }
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -281,7 +304,7 @@ void Wavefunction::set_nspin(int nspin)
 
   compute_nst();
   allocate();
-  resize(cell_,refcell_,ecut_);
+  resize();
   init();
   update_occ(0.0);
 }
@@ -306,7 +329,7 @@ void Wavefunction::set_deltaspin(int deltaspin)
   deltaspin_ = deltaspin;
   compute_nst();
   allocate();
-  resize(cell_,refcell_,ecut_);
+  resize();
   init();
   update_occ(0.0);
 }
@@ -345,7 +368,7 @@ void Wavefunction::set_nrowmax(int n)
   create_contexts();
   compute_nst();
   allocate();
-  resize(cell_,refcell_,ecut_);
+  resize();
   init();
 }
 
@@ -421,7 +444,7 @@ void Wavefunction::del_kpoint(D3vector kpoint)
   kpoint_.erase(pk);
   weight_.erase(pw);
   allocate();
-  resize(cell_,refcell_,ecut_);
+  resize();
   init();
   update_occ(0.0);
 }
diff --git a/src/qb.C b/src/qb.C
index 596ebc2..66774a7 100644
--- a/src/qb.C
+++ b/src/qb.C
@@ -108,6 +108,7 @@ using namespace std;
 #include "NetCharge.h"
 #include "Nrowmax.h"
 #include "Nspin.h"
+#include "Occ.h"
 #include "RefCell.h"
 #include "ScfTol.h"
 #include "Stress.h"
@@ -301,6 +302,7 @@ int main(int argc, char **argv, char **envp)
   ui.addVar(new NetCharge(s));
   ui.addVar(new Nrowmax(s));
   ui.addVar(new Nspin(s));
+  ui.addVar(new Occ(s));
   ui.addVar(new Dspin(s));
   ui.addVar(new RefCell(s));
   ui.addVar(new ScfTol(s));