diff --git a/src/qball/CurrentDensity.cc b/src/qball/CurrentDensity.cc
index 704eb793..b1ab0341 100644
--- a/src/qball/CurrentDensity.cc
+++ b/src/qball/CurrentDensity.cc
@@ -36,7 +36,7 @@ CurrentDensity::CurrentDensity(const Sample& s, const Wavefunction & wf):
   
 }
 
-void CurrentDensity::update_current(EnergyFunctional & energy_functional, const Wavefunction & dwf){
+void CurrentDensity::update_current(EnergyFunctional & energy_functional, const Wavefunction & dwf,bool output){
 
   Wavefunction rwf(wf_);
   Wavefunction rdwf(wf_);
@@ -100,10 +100,11 @@ void CurrentDensity::update_current(EnergyFunctional & energy_functional, const
   // TODO: Reduce total current over spin
   assert(wf_.nspin() == 1);
 	 
-  if ( wf_.context().onpe0() ){
-    std::cout << "  total_electronic_current:\t" << std::fixed << std::setw( 20 ) << std::setprecision( 12 ) << total_current[0] << '\t' << total_current[1] << '\t' << total_current[2] << std::endl;
+  if (output){
+    if ( wf_.context().onpe0() ){
+        std::cout << "  total_electronic_current:\t" << std::fixed << std::setw( 20 ) << std::setprecision( 12 ) << total_current[0] << '\t' << total_current[1] << '\t' << total_current[2] << std::endl;
+    }
   }
-
 }
 
 void CurrentDensity::plot(const Sample * s, const std::string & filename){
diff --git a/src/qball/CurrentDensity.h b/src/qball/CurrentDensity.h
index a55f5009..0e2125f0 100644
--- a/src/qball/CurrentDensity.h
+++ b/src/qball/CurrentDensity.h
@@ -25,6 +25,8 @@
 #include "Wavefunction.h"
 #include "EnergyFunctional.h"
 #include "Sample.h"
+#include <config.h>
+#include <iomanip>
 
 class CurrentDensity : private ChargeDensity
 {
@@ -41,7 +43,7 @@ class CurrentDensity : private ChargeDensity
   ~CurrentDensity (){
   }
 
-  void update_current(EnergyFunctional & energy, const Wavefunction & dwf);
+  void update_current(EnergyFunctional & energy, const Wavefunction & dwf, bool output=true);
 
   void plot(const Sample *, const std::string &);
   void plot_vtk(const Sample *, const std::string &);
diff --git a/src/qball/EhrenSampleStepper.cc b/src/qball/EhrenSampleStepper.cc
index b0ba5272..9e19f176 100644
--- a/src/qball/EhrenSampleStepper.cc
+++ b/src/qball/EhrenSampleStepper.cc
@@ -277,6 +277,13 @@ void EhrenSampleStepper::step(int niter)
   QB_Pstart(14,scfloop);
 #endif
   tmap["total_niter"].start();
+  //yyf: Update vector field and current from where it stops
+  if(ef_.vp) {
+    ef_.vp->propagate(s_.ctrl.tddt*(s_.ctrl.mditer-1), s_.ctrl.tddt);
+    ef_.vector_potential_changed(compute_stress);
+    currd_.update_current(ef_, dwf,false);
+  }
+
   for ( int iter = 0; iter < niter; iter++ )
   {
 
@@ -357,7 +364,7 @@ void EhrenSampleStepper::step(int niter)
 
     tmap["current"].start();
     currd_.update_current(ef_, dwf);
-    tmap["current"].start();
+    tmap["current"].stop();
 
     if(ef_.vp && oncoutpe){
       std::cout << "<!-- vector_potential: " << ef_.vp->value() << " -->\n";
@@ -477,7 +484,7 @@ void EhrenSampleStepper::step(int niter)
     wf_stepper->preupdate();
     tmap["preupdate"].stop();
 
-    if(ef_.vp) ef_.vp->propagate(s_.ctrl.tddt*(iter + 1), s_.ctrl.tddt);
+    if(ef_.vp) ef_.vp->propagate(s_.ctrl.tddt*s_.ctrl.mditer, s_.ctrl.tddt);
     
     tmap["ionic"].start();
     if ( atoms_move )