Cfl split. Phat correlator.

src/io.jl

@@ -123,7 +123,7 @@ function save_data()
         BDIO_write!(fb, [dpar.m0, dpar.csw, params["Fermion"]["Cfl"]])
         BDIO_write!(fb, [dpar.th[i] for i in 1:4])
         BDIO_write!(fb, [convert(Int32,length(flow_times))])
-        BDIO_write!(fb, flow_times)
+        length(flow_times) > 0 ? BDIO_write!(fb, flow_times) : nothing
         BDIO_write!(fb, [convert(Int32,N_noise)])
 
         BDIO_write_hash!(fb)
@@ -141,6 +141,7 @@ function save_data()
             BDIO_write!(fb,ap_corr_t[:,noi,fl])
             BDIO_write!(fb,[Quark_cond[noi,fl]])
             BDIO_write!(fb,[Quark_cond_cfl[noi,fl]])
+            BDIO_write!(fb,[phat_t[noi,fl]])
         end
     end
 
@@ -152,8 +153,6 @@ function save_data()
 end
 
 
-
-
 """
     function load_gauge_field()
 

src/mc.jl

@@ -23,10 +23,12 @@ function Gauge_update()
     println(log_file,"Performing "*string(params["HMC"]["dmeas"])*" HMC steps...")
     flush(log_file)
 
+    global dH = 0.0
+    global acc = true
+
     for _ in 1:params["HMC"]["dmeas"]
         dH,acc = HMC!(U,intsch,lp,gp,ymws)
-        acc ? Nmc_acc += 1 : println(log_file,"Rejected HCM step")
-        acc ? Nmc_acc += 1  : nothing
+        acc ? (global Nmc_acc += 1) : println(log_file,"Rejected HCM step")
     end
     global MCid = MCid + 1
     println(log_file,"Updated to Gauge configuration n"*string(MCid))

src/meas.jl

@@ -21,6 +21,7 @@ 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 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));
     return nothing
@@ -34,8 +35,8 @@ Stores the two point function in the variables 'pp_corr_t0', 'ap_corr_t0', 'pp_c
 """
 function two_pt()
 
-        println(log_file,"Measuring 2pt...")
-        flush(log_file)
+    println(log_file,"Measuring 2pt...")
+    flush(log_file)
     for noi in 1:N_noise
 
         propagator!(psi,U, dpar, dws, lp,params["Measurements"]["maxiter"], params["Measurements"]["tolerance"],1)
@@ -92,15 +93,16 @@ end
 """
 function one_pt()
 
-Stores the one point function in the variable 'Quark_cond'.
+Stores the one point function in the variable 'Quark_cond' and 'Quark_cond_cfl'.
+The final result is 'Quark_cond' + C_fl*'Quark_cond_cfl'
 """
 function one_pt()
 
     copyto!(U,U_CPU)
 
     for fl in 1:length(flow_times)
-            println(log_file,"Measuring 1pt at t = "*string(flow_times[fl])*"...")
-            flush(log_file)
+        println(log_file,"Measuring 1pt at t = "*string(flow_times[fl])*"...")
+        flush(log_file)
 
         for noi in 1:N_noise
 
@@ -131,3 +133,42 @@ function one_pt()
     return nothing
 end
 
+
+
+
+
+"""
+function two_pt_rnd()
+
+Stores the random field evolved in flow time in the variable 'phat_t'
+"""
+function two_pt_rnd()
+
+    copyto!(U,U_CPU)
+
+    println(log_file,"Measuring PhatPt correlator...")
+    flush(log_file)
+    for noi in 1:N_noise
+
+        pfrandomize!(psi,lp)
+
+        delta_t = 0.0
+        neps = int.eps_ini
+        for fl in 1:length(flow_times)
+            delta_t = flow_times[fl] - delta_t
+            _,epslist = flw_adapt(U, psi, int, delta_t, neps, gp, dpar, lp, ymws, dws)
+            neps = epslist[end]
+
+            phat_t[noi,fl] = sum(norm2.(psi))
+        end
+
+        copyto!(U,U_CPU)
+    end
+
+    phat_t .= phat_t./prod(lp.iL)
+
+    println(log_file,"Finished measuring PhatPt")
+    flush(log_file)
+
+    return nothing
+end

src/utils.jl

@@ -31,12 +31,14 @@ function read_ff(name::String)
     TvecR = Vector{Float64}(undef,iL[4]);
     TvecC = Vector{ComplexF64}(undef,iL[4]);
     CN = Vector{ComplexF64}(undef,1);
+    RN = Vector{Float64}(undef,1);
 
     pp_corr = Array{Float64}(undef,1,iL[4],N_noise);
     ap_corr = Array{ComplexF64}(undef,1,iL[4],N_noise);
     pp_corr_t = Array{Float64}(undef,1,iL[4],N_noise,nflow);
     ap_corr_t = Array{ComplexF64}(undef,1,iL[4],N_noise,nflow);
     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);
 
     pp_corr_mc = Array{Float64}(undef,0,iL[4],N_noise);
@@ -44,6 +46,7 @@ function read_ff(name::String)
     pp_corr_t_mc = Array{Float64}(undef,0,iL[4],N_noise,nflow);
     ap_corr_t_mc = Array{ComplexF64}(undef,0,iL[4],N_noise,nflow);
     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);
 
     while BDIO_get_uinfo(file) != 8
@@ -70,6 +73,8 @@ function read_ff(name::String)
                 Sigma[1,noi,fl]= CN[1]
                 BDIO_read(file,CN)
                 Sigma_cfl[1,noi,fl]= CN[1]
+                BDIO_read(file,RN)
+                Phat_t[1,noi,fl]= RN[1]
             end
         end
 
@@ -79,6 +84,7 @@ function read_ff(name::String)
         ap_corr_t_mc=[ap_corr_t_mc;ap_corr_t]
         Sigma_mc=[Sigma_mc;Sigma]
         Sigma_cfl_mc=[Sigma_cfl_mc;Sigma_cfl]
+        Phat_t_mc=[Phat_t_mc;Phat_t]
 
         BDIO_seek!(file)
         BDIO_seek!(file)
@@ -86,7 +92,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
+    return pp_corr_mc,ap_corr_mc,pp_corr_t_mc,ap_corr_t_mc,Sigma_mc,Sigma_cfl_mc,Phat_t_mc
 end
 
 
@@ -97,7 +103,7 @@ 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 = read_ff(file)
+    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]
@@ -113,5 +119,6 @@ function uwff(file::String)
     global Sigma = [uwreal(real.(sum(Sigma_mc,dims = 2)[:,1,i]./Ns),runame) for i in 1:Nfl]
     global Sigma_cfl = [uwreal(real.(sum(Sigma_cfl_mc,dims = 2)[:,1,i]./Ns),runame) for i in 1:Nfl]
 
+    global Phat_t = [uwreal(sum(Phat_t_mc,dims = 2)[:,1,i]./Ns,runame) for i in 1:Nfl]
     return nothing
 end