diff --git a/src/io.jl b/src/io.jl
index 0039bb10280c5aa7c2f02f9ac77000946cb189e2..c21841763fbdda4097e1ee036c6084757259eee0 100644
--- a/src/io.jl
+++ b/src/io.jl
@@ -144,6 +144,7 @@ function save_data()
             BDIO_write!(fb,[phat_t[noi,fl]])
         end
     end
+    BDIO_write!(fb,Eoft)
 
     BDIO_write_hash!(fb)
 
diff --git a/src/mc.jl b/src/mc.jl
index 2e9b62882e975311281ac11bac39cbf54d417716..96bcf06510fc26f03293b13f92466b890548c827 100644
--- a/src/mc.jl
+++ b/src/mc.jl
@@ -33,7 +33,7 @@ function Gauge_update()
     global MCid = MCid + 1
     println(log_file,"Updated to Gauge configuration n"*string(MCid))
     println(log_file,"dH   = ",dH)
-    println(log_file,"acc  = ",Nmc_acc/(params["HMC"]["dmeas"] * Mcid))
+    println(log_file,"acc  = ",Nmc_acc/(params["HMC"]["dmeas"] * MCid))
     println(log_file,"Plaq = ", plaquette(U, lp, gp, ymws))
     println(log_file,"Qtop = ", Qtop(U, gp, lp, ymws))
 
diff --git a/src/meas.jl b/src/meas.jl
index 466bcbb25e4b947920cb67f1dcf615396616623a..3f63dcccee5ab91c8decc4892cc6c758d4a00b64 100644
--- a/src/meas.jl
+++ b/src/meas.jl
@@ -21,6 +21,8 @@ function load_fields()
     global pp_corr_t = fill(zero(Float64),(lp.iL[4],N_noise,length(flow_times)));
     global ap_corr_t = fill(zero(ComplexF64),(lp.iL[4],N_noise,length(flow_times)));
 
+    global Eoft = Array{Complex{Float64}}(undef,1+length(flow_times));
+
     global phat_t = Array{Float64}(undef,N_noise,length(flow_times));
     global Quark_cond = Array{Complex{Float64}}(undef,N_noise,length(flow_times));
     global Quark_cond_cfl = Array{Complex{Float64}}(undef,N_noise,length(flow_times));
@@ -44,6 +46,8 @@ function two_pt()
         pp_density .= Array(norm2.(psi))
         ap_density .= Array(dot.(psi,dmul.(Gamma{4},psi)))
 
+        noi == 1 ? Eoft[1] = Eoft_plaq(U, gp, lp, ymws) : nothing
+
         pp_corr_t0[:,noi] .= zero(Float64)
         ap_corr_t0[:,noi] .= zero(ComplexF64)
         for t in 1:lp.iL[4]
@@ -61,6 +65,8 @@ function two_pt()
             _,epslist = flw_adapt(U, psi, int, delta_t, neps, gp, dpar, lp, ymws, dws)
             neps = epslist[end]
 
+            noi == 1 ? Eoft[fl+1] = Eoft_plaq(U, gp, lp, ymws) : nothing
+
             pp_density .= Array(norm2.(psi))
             ap_density .= Array(dot.(psi,dmul.(Gamma{4},psi)))
             pp_corr_t[:,noi,fl] .= zero(Float64)
diff --git a/src/utils.jl b/src/utils.jl
index 590bfa6d93bbb20b9fba8f8472929d13478891a5..e9c0210cd897cbe4bbaaf5a5daa7834703e14ccb 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -40,6 +40,7 @@ function read_ff(name::String)
     Sigma = Array{ComplexF64}(undef,1,N_noise,nflow);
     Phat_t = Array{Float64}(undef,1,N_noise,nflow);
     Sigma_cfl = Array{ComplexF64}(undef,1,N_noise,nflow);
+    Eoft = Array{Complex{Float64}}(undef,1+nflow);
 
     pp_corr_mc = Array{Float64}(undef,0,iL[4],N_noise);
     ap_corr_mc = Array{ComplexF64}(undef,0,iL[4],N_noise);
@@ -48,6 +49,7 @@ function read_ff(name::String)
     Sigma_mc = Array{ComplexF64}(undef,0,N_noise,nflow);
     Phat_t_mc = Array{Float64}(undef,0,N_noise,nflow);
     Sigma_cfl_mc = Array{ComplexF64}(undef,0,N_noise,nflow);
+    Eoft_mc = Array{Complex{Float64}}(undef,0,1+nflow);
 
     while BDIO_get_uinfo(file) != 8
         BDIO_seek!(file)
@@ -77,6 +79,8 @@ function read_ff(name::String)
                 Phat_t[1,noi,fl]= RN[1]
             end
         end
+        BDIO_read(file,Eoft)
+        Eoft_mc=[Eoft_mc;Eoft]
 
         pp_corr_mc=[pp_corr_mc;pp_corr]
         ap_corr_mc=[ap_corr_mc;ap_corr]
@@ -92,7 +96,7 @@ function read_ff(name::String)
 
     BDIO_close!(file)
 
-    return pp_corr_mc,ap_corr_mc,pp_corr_t_mc,ap_corr_t_mc,Sigma_mc,Sigma_cfl_mc,Phat_t_mc
+    return Eoft_mc,pp_corr_mc,ap_corr_mc,pp_corr_t_mc,ap_corr_t_mc,Sigma_mc,Sigma_cfl_mc,Phat_t_mc
 end
 
 
@@ -103,13 +107,15 @@ Loads in the variables 'pp_corr','ap_corr','pp_corr_t','ap_corr_t' and 'Sigma' t
 'uwreal' quantities. The indices are, when corresponding, eucliden time and flow times
 """
 function uwff(file::String)
-    pp_corr_mc,ap_corr_mc,pp_corr_t_mc,ap_corr_t_mc,Sigma_mc,Phat_t_mc = read_ff(file)
+    Eoft_mc,pp_corr_mc,ap_corr_mc,pp_corr_t_mc,ap_corr_t_mc,Sigma_mc,Phat_t_mc = read_ff(file)
     runame = String(split(split(file,"/")[end],".")[1])
 
     T = size(pp_corr_mc)[2]
     Ns = size(pp_corr_mc)[3]
     Nfl = size(pp_corr_t_mc)[4]
 
+    global Eoft = [uwreal(real.(Eoft_mc[:,i]),runame) for i in 1:T+1]
+
     global pp_corr = [uwreal(sum(pp_corr_mc,dims = 3)[:,i,1]./Ns,runame) for i in 1:T]
     global ap_corr = [uwreal(real.(sum(ap_corr_mc,dims = 3)[:,i,1]./Ns),runame) for i in 1:T]