All the correlation functions added.

sfcf.jl

@@ -33,51 +33,15 @@ println("--------------------------------------------------")
 println("Reading input file from:", parsed_args["i"], "...")
 params = TOML.parsefile(parsed_args["i"])
 
-cG = try
-    params["Space"]["cG"]
-catch
-    try
-        println("Using default value for cG")
-        1.0 - 0.089*6/params["Fermion"]["beta"]
-    catch
-        error("Either cG of beta has to be defined in the input file")
-    end
-end
-
-ct = try
-    params["Fermion"]["ct"]
-catch
-    try
-        error("ct not defined")
-        println("Using default value for ct")
-    catch
-        error("Either ct of beta has to be defined in the input file")
-    end
-end
-
-csw = try
-    params["Fermion"]["csw"]
-catch
-    try
-        error("csw not defined")
-        println("Using default value for csw")
-    catch
-        error("Either csw of beta has to be defined in the input file")
-    end
-end
-
 lp = SpaceParm{4}(tuple(params["Space"]["size"]...), tuple(params["Space"]["blocks"]...),BC_SF_ORBI, (0,0,0,0,0,0))
-gp = GaugeParm{Float64}(SU3{Float64},params["Fermion"]["beta"],1.0,(cG,0.0),params["Space"]["phiT"],lp.iL);
+gp = GaugeParm{Float64}(SU3{Float64},params["Fermion"]["beta"],1.0,(params["Space"]["cG"],0.0),params["Space"]["phiT"],lp.iL);
-dpar = DiracParam{Float64}(SU3fund,(1/(2*params["Fermion"]["kappa"])) - 4,csw,ntuple(i -> exp((i!=4)*im*params["Fermion"]["theta"]/lp.iL[i]),4),0.0,ct);
+dpar = DiracParam{Float64}(SU3fund,(1/(2*params["Fermion"]["kappa"])) - 4,params["Fermion"]["csw"],ntuple(i -> exp((i!=4)*im*params["Fermion"]["theta"]/lp.iL[i]),4),0.0,params["Fermion"]["ct"]);
 
 dws = DiracWorkspace(SU3fund,Float64,lp);
 ymws = YMworkspace(SU3,Float64,lp);
 
 println("Reading gauge field from: ", parsed_args["c"], "...")
 
-
-
-
 if !parsed_args["cern"]
     U,_ = read_cnfg(parsed_args["c"])
 else 
@@ -168,44 +132,43 @@ println("--------------------------------------------------")
 
 psi = scalar_field(Spinor{4,SU3fund{Float64}},lp)
 
-f1 = 0.0
 
 println("Computing propagators from t=0 to the bulk")
 niter = bndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 1, 1)
-println("CG converged for s=1 and c=1 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
+println("CG converged for c=1 and s=1 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
 A11 = Array(psi)
-f1 += norm2(bndtobnd(psi, U, dpar, dws, lp))
+C11 = bndtobnd(psi, U, dpar, dws, lp)
 
 niter = bndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 2, 1)
 println("CG converged for c=2 and s=1 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
 A21 = Array(psi)
-f1 += norm2(bndtobnd(psi, U, dpar, dws, lp))
+C21 = bndtobnd(psi, U, dpar, dws, lp)
 
 niter = bndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 3, 1)
 println("CG converged for c=3 and s=1 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
 A31 = Array(psi)
-f1 += norm2(bndtobnd(psi, U, dpar, dws, lp))
+C31 = bndtobnd(psi, U, dpar, dws, lp)
 
 niter = bndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 1, 2)
 println("CG converged for c=1 and s=2 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
 A12 = Array(psi)
-f1 += norm2(bndtobnd(psi, U, dpar, dws, lp))
+C12 = bndtobnd(psi, U, dpar, dws, lp)
 
 niter = bndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 2, 2)
 println("CG converged for c=2 and s=2 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
 A22 = Array(psi)
-f1 += norm2(bndtobnd(psi, U, dpar, dws, lp))
+C22 = bndtobnd(psi, U, dpar, dws, lp)
 
 niter = bndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 3, 2)
 println("CG converged for c=3 and s=2 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
 A32 = Array(psi)
-f1 += norm2(bndtobnd(psi, U, dpar, dws, lp))
+C32 = bndtobnd(psi, U, dpar, dws, lp)
 
 
 
 println("Computing propagators from t=T to the bulk")
 niter = Tbndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 1, 1)
-println("CG converged for s=1 and c=1 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
+println("CG converged for c=1 and s=1 after ",niter," iterations with absolute residue ", CUDA.mapreduce(x -> norm(x), +, dws.sr))
 B11 = Array(psi)
 
 niter = Tbndpropagator!(psi, U, dpar, dws, lp, params["Solver"]["maxiter"], params["Solver"]["tolerance"], 2, 1)
@@ -231,6 +194,28 @@ B32 = Array(psi)
 println("--------------------------------------------------")
 println("Computing correlators")
 
+
+
+println("Computing f1...")
+file = open("./output/"*basename(parsed_args["c"])*".f1","w+")
+f1 = (norm2(C11) + norm2(C21) + norm2(C31) + norm2(C12) + norm2(C22) + norm2(C32))
+println(file,"    ",f1)
+flush(file)
+close(file)
+
+println("Computing k1...")
+file = open("./output/"*basename(parsed_args["c"])*".k1","w+")
+k1 = (dot(C11,dmul(Gamma{6},imm(C12))) + dot(C11,dmul(Gamma{7},-C12)) + dot(C11,dmul(Gamma{8}, imm(C11)))
+     +dot(C12,dmul(Gamma{6},imm(C11))) + dot(C12,dmul(Gamma{7}, C11)) + dot(C12,dmul(Gamma{8},mimm(C12)))
+     +dot(C21,dmul(Gamma{6},imm(C22))) + dot(C21,dmul(Gamma{7},-C22)) + dot(C21,dmul(Gamma{8}, imm(C21)))
+     +dot(C22,dmul(Gamma{6},imm(C21))) + dot(C22,dmul(Gamma{7}, C21)) + dot(C22,dmul(Gamma{8},mimm(C22)))
+     +dot(C31,dmul(Gamma{6},imm(C32))) + dot(C31,dmul(Gamma{7},-C32)) + dot(C31,dmul(Gamma{8}, imm(C31)))
+     +dot(C32,dmul(Gamma{6},imm(C31))) + dot(C32,dmul(Gamma{7}, C31)) + dot(C32,dmul(Gamma{8},mimm(C32))) )/3
+
+println(file,"    ",k1)
+flush(file)
+close(file)
+
 println("Computing fP...")
 file = open("./output/"*basename(parsed_args["c"])*".fP","w+")
 fP = Vector{Complex{Float64}}(undef,lp.iL[4])
@@ -238,9 +223,9 @@ for t in 1:lp.iL[4]
     fP[t] = 0.0
     for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
         b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
-        fP[t] += norm2(A11[b,r]) + norm2(A21[b,r]) + norm2(A31[b,r]) + norm2(A12[b,r]) + norm2(A22[b,r]) + norm2(A32[b,r])
+        fP[t] += (norm2(A11[b,r]) + norm2(A21[b,r]) + norm2(A31[b,r]) + norm2(A12[b,r]) + norm2(A22[b,r]) + norm2(A32[b,r]))/prod(lp.iL[1:3])
     end end end
-    println(file,t,"    ",fP[t]/prod(lp.iL[1:3]))
+    println(file,t,"    ",fP[t])
 end
 flush(file)
 close(file)
@@ -252,11 +237,11 @@ for t in 1:lp.iL[4]
     fA[t] = 0.0
     for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
         b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
-        fA[t] -= dot(A11[b,r],dmul(Gamma{4},A11[b,r])) + dot(A21[b,r],dmul(Gamma{4},A21[b,r])) + dot(A31[b,r],dmul(Gamma{4},A31[b,r]))
+        fA[t] -= (dot(A11[b,r],dmul(Gamma{4},A11[b,r])) + dot(A21[b,r],dmul(Gamma{4},A21[b,r])) + dot(A31[b,r],dmul(Gamma{4},A31[b,r]))
-        + dot(A12[b,r],dmul(Gamma{4},A12[b,r])) + dot(A22[b,r],dmul(Gamma{4},A22[b,r])) + dot(A32[b,r],dmul(Gamma{4},A32[b,r]))
+        + dot(A12[b,r],dmul(Gamma{4},A12[b,r])) + dot(A22[b,r],dmul(Gamma{4},A22[b,r])) + dot(A32[b,r],dmul(Gamma{4},A32[b,r])))/prod(lp.iL[1:3])
 
     end end end
-    println(file,t,"    ",fA[t]/prod(lp.iL[1:3]))
+    println(file,t,"    ",fA[t])
 end
 flush(file)
 close(file)
@@ -268,9 +253,9 @@ for t in 1:lp.iL[4]
     gP[t] = 0.0
     for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
         b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
-        gP[t] += norm2(B11[b,r]) + norm2(B21[b,r]) + norm2(B31[b,r]) + norm2(B12[b,r]) + norm2(B22[b,r]) + norm2(B32[b,r])
+        gP[t] += (norm2(B11[b,r]) + norm2(B21[b,r]) + norm2(B31[b,r]) + norm2(B12[b,r]) + norm2(B22[b,r]) + norm2(B32[b,r]))/prod(lp.iL[1:3])
     end end end
-    println(file,t,"    ",gP[t]/prod(lp.iL[1:3]))
+    println(file,t,"    ",gP[t])
 end
 flush(file)
 close(file)
@@ -283,15 +268,97 @@ for t in 1:lp.iL[4]
     gA[t] = 0.0
     for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
         b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
-        gA[t] += dot(B11[b,r],dmul(Gamma{4},B11[b,r])) + dot(B21[b,r],dmul(Gamma{4},B21[b,r])) + dot(B31[b,r],dmul(Gamma{4},B31[b,r]))
+        gA[t] += (dot(B11[b,r],dmul(Gamma{4},B11[b,r])) + dot(B21[b,r],dmul(Gamma{4},B21[b,r])) + dot(B31[b,r],dmul(Gamma{4},B31[b,r]))
-        + dot(B12[b,r],dmul(Gamma{4},B12[b,r])) + dot(B22[b,r],dmul(Gamma{4},B22[b,r])) + dot(B32[b,r],dmul(Gamma{4},B32[b,r]))
+        + dot(B12[b,r],dmul(Gamma{4},B12[b,r])) + dot(B22[b,r],dmul(Gamma{4},B22[b,r])) + dot(B32[b,r],dmul(Gamma{4},B32[b,r])))/prod(lp.iL[1:3])
     end end end
-    println(file,t,"    ",gA[t]/prod(lp.iL[1:3]))
+    println(file,t,"    ",gA[t])
+end
+flush(file)
+close(file)
+
+#### print f1
+
+println("Computing kV...")
+file = open("./output/"*basename(parsed_args["c"])*".kV","w+")
+kv = Vector{Complex{Float64}}(undef,lp.iL[4])
+for t in 1:lp.iL[4]
+    kv[t] = 0.0
+    for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
+        b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
+        kv[t] += ( dot(A11[b,r],dmul(Gamma{6},imm(A12[b,r]))) + dot(A11[b,r],dmul(Gamma{7},-A12[b,r])) + dot(A11[b,r],dmul(Gamma{8}, imm(A11[b,r])))
+                  +dot(A12[b,r],dmul(Gamma{6},imm(A11[b,r]))) + dot(A12[b,r],dmul(Gamma{7}, A11[b,r])) + dot(A12[b,r],dmul(Gamma{8},mimm(A12[b,r])))
+                  +dot(A21[b,r],dmul(Gamma{6},imm(A22[b,r]))) + dot(A21[b,r],dmul(Gamma{7},-A22[b,r])) + dot(A21[b,r],dmul(Gamma{8}, imm(A21[b,r])))
+                  +dot(A22[b,r],dmul(Gamma{6},imm(A21[b,r]))) + dot(A22[b,r],dmul(Gamma{7}, A21[b,r])) + dot(A22[b,r],dmul(Gamma{8},mimm(A22[b,r])))
+                  +dot(A31[b,r],dmul(Gamma{6},imm(A32[b,r]))) + dot(A31[b,r],dmul(Gamma{7},-A32[b,r])) + dot(A31[b,r],dmul(Gamma{8}, imm(A31[b,r])))
+                  +dot(A32[b,r],dmul(Gamma{6},imm(A31[b,r]))) + dot(A32[b,r],dmul(Gamma{7}, A31[b,r])) + dot(A32[b,r],dmul(Gamma{8},mimm(A32[b,r]))) )/(3*prod(lp.iL[1:3]))
+    end end end
+    println(file,t,"    ",kv[t])
+end
+flush(file)
+close(file)
+
+
+println("Computing lV...")
+file = open("./output/"*basename(parsed_args["c"])*".lV","w+")
+lv = Vector{Complex{Float64}}(undef,lp.iL[4])
+for t in 1:lp.iL[4]
+    lv[t] = 0.0
+    for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
+        b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
+        lv[t] += ( dot(B11[b,r],dmul(Gamma{6},mimm(B12[b,r]))) + dot(B11[b,r],dmul(Gamma{7}, B12[b,r])) + dot(B11[b,r],dmul(Gamma{8},mimm(B11[b,r])))
+                  +dot(B12[b,r],dmul(Gamma{6},mimm(B11[b,r]))) + dot(B12[b,r],dmul(Gamma{7},-B11[b,r])) + dot(B12[b,r],dmul(Gamma{8}, imm(B12[b,r])))
+                  +dot(B21[b,r],dmul(Gamma{6},mimm(B22[b,r]))) + dot(B21[b,r],dmul(Gamma{7}, B22[b,r])) + dot(B21[b,r],dmul(Gamma{8},mimm(B21[b,r])))
+                  +dot(B22[b,r],dmul(Gamma{6},mimm(B21[b,r]))) + dot(B22[b,r],dmul(Gamma{7},-B21[b,r])) + dot(B22[b,r],dmul(Gamma{8}, imm(B22[b,r])))
+                  +dot(B31[b,r],dmul(Gamma{6},mimm(B32[b,r]))) + dot(B31[b,r],dmul(Gamma{7}, B32[b,r])) + dot(B31[b,r],dmul(Gamma{8},mimm(B31[b,r])))
+                  +dot(B32[b,r],dmul(Gamma{6},mimm(B31[b,r]))) + dot(B32[b,r],dmul(Gamma{7},-B31[b,r])) + dot(B32[b,r],dmul(Gamma{8}, imm(B32[b,r]))) )/(3*prod(lp.iL[1:3]))
+    end end end
+    println(file,t,"    ",lv[t])
 end
 flush(file)
 close(file)
 
 
 
+println("Computing kT...")
+file = open("./output/"*basename(parsed_args["c"])*".kT","w+")
+kt = Vector{Complex{Float64}}(undef,lp.iL[4])
+for t in 1:lp.iL[4]
+    kt[t] = 0.0
+    for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
+        b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
+        kt[t] +=im*(dot(A11[b,r],dmul(Gamma{10},imm(dmul(Gamma{5},A12[b,r])))) + dot(A11[b,r],dmul(Gamma{11},dmul(Gamma{5},-A12[b,r]))) + dot(A11[b,r],dmul(Gamma{12},dmul(Gamma{5}, imm(A11[b,r]))))
+                   +dot(A12[b,r],dmul(Gamma{10},imm(dmul(Gamma{5},A11[b,r])))) + dot(A12[b,r],dmul(Gamma{11},dmul(Gamma{5}, A11[b,r]))) + dot(A12[b,r],dmul(Gamma{12},dmul(Gamma{5},mimm(A12[b,r]))))
+                   +dot(A21[b,r],dmul(Gamma{10},imm(dmul(Gamma{5},A22[b,r])))) + dot(A21[b,r],dmul(Gamma{11},dmul(Gamma{5},-A22[b,r]))) + dot(A21[b,r],dmul(Gamma{12},dmul(Gamma{5}, imm(A21[b,r]))))
+                   +dot(A22[b,r],dmul(Gamma{10},imm(dmul(Gamma{5},A21[b,r])))) + dot(A22[b,r],dmul(Gamma{11},dmul(Gamma{5}, A21[b,r]))) + dot(A22[b,r],dmul(Gamma{12},dmul(Gamma{5},mimm(A22[b,r]))))
+                   +dot(A31[b,r],dmul(Gamma{10},imm(dmul(Gamma{5},A32[b,r])))) + dot(A31[b,r],dmul(Gamma{11},dmul(Gamma{5},-A32[b,r]))) + dot(A31[b,r],dmul(Gamma{12},dmul(Gamma{5}, imm(A31[b,r]))))
+                   +dot(A32[b,r],dmul(Gamma{10},imm(dmul(Gamma{5},A31[b,r])))) + dot(A32[b,r],dmul(Gamma{11},dmul(Gamma{5}, A31[b,r]))) + dot(A32[b,r],dmul(Gamma{12},dmul(Gamma{5},mimm(A32[b,r])))) )/(3*prod(lp.iL[1:3]))
+    end end end
+    println(file,t,"    ",kt[t])
+end
+flush(file)
+close(file)
+
+
+
+println("Computing lT...")
+file = open("./output/"*basename(parsed_args["c"])*".lT","w+")
+lt = Vector{Complex{Float64}}(undef,lp.iL[4])
+for t in 1:lp.iL[4]
+    lt[t] = 0.0
+    for i in 1:lp.iL[1] for j in 1:lp.iL[2] for k in 1:lp.iL[3]
+        b,r = point_index(CartesianIndex{lp.ndim}((i,j,k,t)),lp)
+        lt[t] +=-im*(dot(B11[b,r],dmul(Gamma{10},mimm(dmul(Gamma{5},B12[b,r])))) + dot(B11[b,r],dmul(Gamma{11},dmul(Gamma{5}, B12[b,r]))) + dot(B11[b,r],dmul(Gamma{12},dmul(Gamma{5},mimm(B11[b,r]))))
+                    +dot(B12[b,r],dmul(Gamma{10},mimm(dmul(Gamma{5},B11[b,r])))) + dot(B12[b,r],dmul(Gamma{11},dmul(Gamma{5},-B11[b,r]))) + dot(B12[b,r],dmul(Gamma{12},dmul(Gamma{5}, imm(B12[b,r]))))
+                    +dot(B21[b,r],dmul(Gamma{10},mimm(dmul(Gamma{5},B22[b,r])))) + dot(B21[b,r],dmul(Gamma{11},dmul(Gamma{5}, B22[b,r]))) + dot(B21[b,r],dmul(Gamma{12},dmul(Gamma{5},mimm(B21[b,r]))))
+                    +dot(B22[b,r],dmul(Gamma{10},mimm(dmul(Gamma{5},B21[b,r])))) + dot(B22[b,r],dmul(Gamma{11},dmul(Gamma{5},-B21[b,r]))) + dot(B22[b,r],dmul(Gamma{12},dmul(Gamma{5}, imm(B22[b,r]))))
+                    +dot(B31[b,r],dmul(Gamma{10},mimm(dmul(Gamma{5},B32[b,r])))) + dot(B31[b,r],dmul(Gamma{11},dmul(Gamma{5}, B32[b,r]))) + dot(B31[b,r],dmul(Gamma{12},dmul(Gamma{5},mimm(B31[b,r]))))
+                    +dot(B32[b,r],dmul(Gamma{10},mimm(dmul(Gamma{5},B31[b,r])))) + dot(B32[b,r],dmul(Gamma{11},dmul(Gamma{5},-B31[b,r]))) + dot(B32[b,r],dmul(Gamma{12},dmul(Gamma{5}, imm(B32[b,r])))) )/(3*prod(lp.iL[1:3]))
+    end end end
+    println(file,t,"    ",lt[t])
+end
+flush(file)
+close(file)
+
+
 print("\n\n")
 print_timer()