diff --git a/sfcf.jl b/sfcf.jl index 8dfc16b..bead746 100644 --- a/sfcf.jl +++ b/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); -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); +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,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])) - + 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])) + 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])))/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])) - + 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])) + 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])))/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()