diff --git a/sfcf.jl b/sfcf.jl
index 5044562..b2d1e70 100644
--- a/sfcf.jl
+++ b/sfcf.jl
@@ -4,6 +4,8 @@ using TimerOutputs
 using ArgParse
 using CUDA
 using InteractiveUtils
+using BDIO
+using MD5
 
 include("./src/io.jl")
 include("./src/meas.jl")
@@ -26,10 +28,10 @@ end
 
 @timeit "Correlators contractions" compute_correlators()
 
+@timeit "Correlators saving" save_correlators()
 
 print("\n\n")
 print_timer(log_file)
 
 flush(log_file)
 close(log_file)
-
diff --git a/src/io.jl b/src/io.jl
index f507c41..ca4f1b6 100644
--- a/src/io.jl
+++ b/src/io.jl
@@ -1,8 +1,10 @@
 using LatticeGPU
 using TOML
-using TimerOutputs
 using ArgParse
 using CUDA
+using BDIO
+using InteractiveUtils
+using MD5
 
 
 """
@@ -51,10 +53,14 @@ function load_gauge_field()
     println(log_file,"\nReading gauge field with LGPU native convention from: ", parsed_args["c"], "...")
         U,_ = read_cnfg(parsed_args["c"])
     else
-    println(log_file,"\nReading gauge field with CERN convention from: ", parsed_args["c"], "...\n")
+    println(log_file,"\nReading gauge field with CERN convention from: ", parsed_args["c"], "...")
         U = read_cnfg_cern(parsed_args["c"],lp)
     end
 
+    f = open(parsed_args["c"],"r")
+    println(log_file,"MD5 checksum of gauge config: ",bytes2hex(md5(f),"\n"))
+    close(f)
+
     Csw!(dws, U, gp, lp)
     return U
 end
@@ -176,3 +182,54 @@ function write_log()
     return nothing
 end
 
+
+function save_correlators()
+
+    ihdr = [convert(Int32,1708683512)]
+    fname = "./output"*params["Run"]["name"]*".bdio"
+
+    if isfile(fname)
+        fb = BDIO_open(fname, "a")
+        println(log_file,"Appending output to "*fname)
+    else
+        fb = BDIO_open(fname, "w", "BDIO output from sfcf.jl")
+        println(log_file,"Creating new BDIO output file "*fname)
+
+        BDIO_start_record!(fb, BDIO_BIN_GENERIC, 14)
+        BDIO_write!(fb, ihdr)
+        BDIO_write_hash!(fb)
+
+        BDIO_start_record!(fb, BDIO_BIN_GENERIC, 1)
+        BDIO_write!(fb, [convert(Int32, 4)])
+        BDIO_write!(fb, [convert(Int32, lp.iL[i]) for i in 1:4])
+        BDIO_write!(fb, [convert(Int32, lp.ntw[i]) for i in 1:6])
+        BDIO_write!(fb, [dpar.m0, dpar.csw, dpar.tm, dpar.ct])
+        BDIO_write!(fb, [dpar.th[i] for i in 1:4])
+
+        BDIO_write_hash!(fb)
+    end
+
+
+    BDIO_start_record!(fb, BDIO_BIN_GENERIC, 8)
+    BDIO_write!(fb,basename(parsed_args["c"]))
+    fg = open(parsed_args["c"],"r")
+    BDIO_write!(fb,bytes2hex(md5(fg)))
+    close(fg)
+
+    BDIO_write!(fb,fP)
+    BDIO_write!(fb,fA)
+    BDIO_write!(fb,[f1])
+    BDIO_write!(fb,gP)
+    BDIO_write!(fb,gA)
+    BDIO_write!(fb,kV)
+    BDIO_write!(fb,lV)
+    BDIO_write!(fb,[k1])
+    BDIO_write!(fb,kT)
+    BDIO_write!(fb,lT)
+
+    BDIO_write_hash!(fb)
+
+    BDIO_close!(fb)
+
+    return nothing
+end