Working version. Not consistent fA/gA yet

input/sfcf.in

@@ -1,21 +1,22 @@
 [Run]
-user    = "fperez"
-name    = "SFtest_Setup10"
+user    = "Lattice user"
+name    = "Template run"
 
 [Space]
 size       = [8,8,8,8]
 blocks     = [4,4,4,4]
-phi0   = [-1.047197551196598, 0.0]               # Boundary fields for SF
-phiT   = [-3.141592653589793, 1.047197551196598]
-cG     = 1.0  # Correction for boundaries (SF only)
+phi0   = [0.0, 0.0]               # Boundary fields for SF
+phiT   = [0.0, 0.0]
+cG     = 9.6553053029661895e-01  # Correction for boundaries (SF only)
 
 [Fermion]
-beta   = 6.0   #Used only if cG/csw/ct are not input
-kappa  = 0.1
+beta   = 8.4044000000000008e+00   #Used only if cG/csw/ct are not input
+kappa  = 1.3247670000000000e-01
 theta  = 0.5
-csw    = 1.0
-ct     = 0.9
+csw    = 1.2394261961623441e+00
+ct     = 9.8718528389891003e-01
 
 [Solver]
-tolerance = 1.0e-10
+tolerance = 1.0e-13
 maxiter   = 1000
+

output/randfield.fA

@@ -0,0 +1,8 @@
+1    0.0 + 0.0im
+2    -31.04305231567254 + 2.1252056646610396e-18im
+3    -1.4671112666898436 - 9.031488800098977e-20im
+4    -0.06537188111892875 + 2.716461131624143e-21im
+5    -0.002933795526834015 - 2.365736551901659e-22im
+6    -0.00013043629270008468 + 1.8495241508802853e-23im
+7    -5.875876629837171e-6 + 4.476186102841274e-25im
+8    -2.726581392340069e-7 + 6.563322731438818e-27im

output/randfield.fP

@@ -0,0 +1,8 @@
+1    0.0 + 0.0im
+2    73.67395934956765 + 0.0im
+3    3.3887421125734036 + 0.0im
+4    0.15144928799903828 + 0.0im
+5    0.006828696539887742 + 0.0im
+6    0.0003023150684989386 + 0.0im
+7    1.361849306449979e-5 + 0.0im
+8    6.375793198531621e-7 + 0.0im

output/randfield.gA

@@ -0,0 +1,8 @@
+1    0.0 + 0.0im
+2    -2.7310053780937123e-7 - 4.550149701314796e-27im
+3    -5.8143401218682364e-6 + 1.6206358436860467e-25im
+4    -0.00013056821426116682 - 5.932435955653745e-24im
+5    -0.0029261383446798125 + 7.222320723353623e-23im
+6    -0.06619259853464403 - 3.22931311140702e-21im
+7    -1.4692707713700914 + 2.5710096927128185e-19im
+8    -30.89429497223169 - 3.0001906991747318e-18im

output/randfield.gP

@@ -0,0 +1,8 @@
+1    0.0 + 0.0im
+2    6.435062273546668e-7 + 0.0im
+3    1.3621688998858976e-5 + 0.0im
+4    0.00030680086560982034 + 0.0im
+5    0.006857460528564787 + 0.0im
+6    0.1539180922119846 + 0.0im
+7    3.3981870666381786 + 0.0im
+8    73.17751168291873 + 0.0im

output/sfqcd_L8_b8.4044_k0.1324767_r0_id10n1.fA

@@ -0,0 +1,8 @@
+1    0.0 + 0.0im
+2    -0.6030372212882809 + 8.523269031746397e-20im
+3    -0.6038618620887128 + 1.4463087574367178e-19im
+4    -0.5926760024257145 - 8.533856943587076e-20im
+5    -0.601666224177839 + 4.129285617864714e-21im
+6    -0.584098047171397 - 5.108667463127499e-20im
+7    -0.5976874848953928 + 3.4430565816907224e-20im
+8    -0.6031406468183864 + 2.1996386849010112e-20im

output/sfqcd_L8_b8.4044_k0.1324767_r0_id10n1.fP

@@ -0,0 +1,8 @@
+1    0.0 + 0.0im
+2    4.073707933719405 + 0.0im
+3    3.263441992848985 + 0.0im
+4    2.599329122980939 + 0.0im
+5    2.0726204194655655 + 0.0im
+6    1.7097535104448742 + 0.0im
+7    1.4965323741017251 + 0.0im
+8    1.3535414135779962 + 0.0im

output/sfqcd_L8_b8.4044_k0.1324767_r0_id10n1.gA

@@ -0,0 +1,8 @@
+1    0.0 + 0.0im
+2    -0.6084153970265829 + 5.61159327555974e-21im
+3    -0.5874030271133971 + 4.200754022789296e-20im
+4    -0.5990867578116633 + 1.8767073737603092e-20im
+5    -0.5901710309969539 + 6.858319894799663e-20im
+6    -0.6045239706245472 - 7.273895434546304e-20im
+7    -0.6002224832169679 - 1.4463087574367178e-19im
+8    -0.5874210789132396 + 4.6375053862172944e-20im

output/sfqcd_L8_b8.4044_k0.1324767_r0_id10n1.gP

@@ -0,0 +1,8 @@
+1    0.0 + 0.0im
+2    1.3469498829880007 + 0.0im
+3    1.4652378537401343 + 0.0im
+4    1.6960689316800241 + 0.0im
+5    1.9966909928964056 + 0.0im
+6    2.496926368571815 + 0.0im
+7    3.188758530279535 + 0.0im
+8    3.983946740235404 + 0.0im

sfcf.jl

@@ -18,6 +18,10 @@ function parse_commandline()
         required = true
         arg_type = String
 
+        "--cern"
+        help = "Config written with the export_cnfg_cern() convention"
+        action = :store_true
+
     end
 
     return parse_args(s)
@@ -70,7 +74,81 @@ dws = DiracWorkspace(SU3fund,Float64,lp);
 ymws = YMworkspace(SU3,Float64,lp);
 
 println("Reading gauge field from: ", parsed_args["c"], "...")
-U,_ = read_cnfg(parsed_args["c"])
+
+
+
+
+if !parsed_args["cern"]
+    U,_ = read_cnfg(parsed_args["c"])
+else 
+    U = vector_field(SU3{Float64}, lp);
+    file = open(parsed_args["c"])
+    for t in 1:lp.iL[4]
+        for i in 1:lp.iL[1]
+            for j in 1:lp.iL[2]
+                for k in 1:lp.iL[3]
+
+                    for d in [4,1,2,3]
+
+                        f,r = point_index(CartesianIndex((i,j,k,t)),lp)
+
+                        #a11 !!
+
+                        re11 = read(file,Float64)
+                        co11 = read(file,Float64)
+                        
+                        #a12 !!
+                        
+                        re12 = read(file,Float64)
+                        co12 = read(file,Float64)
+                        
+                        #a13 !!
+
+                        re13 = read(file,Float64)
+                        co13 = read(file,Float64)
+
+                        #a21 !!
+
+                        re21 = read(file,Float64)
+                        co21 = read(file,Float64)
+
+                        #a22 !!
+
+                        re22 = read(file,Float64)
+                        co22 = read(file,Float64)
+
+                        #a23 !!
+
+                        re23 = read(file,Float64)
+                        co23 = read(file,Float64)
+
+                        #a31
+
+                        re31 = read(file,Float64)
+                        co31 = read(file,Float64)
+                        
+                        #a32
+                        
+                        re32 = read(file,Float64)
+                        co32 = read(file,Float64)
+                        
+                        #a33
+
+                        re33 = read(file,Float64)
+                        co33 = read(file,Float64)
+
+                        CUDA.@allowscalar (U[f,d,r] = SU3{Float64}(re11 + im*co11, re12 + im*co12, re13 + im*co13,
+		                                                           re21 + im*co21, re22 + im*co22, re23 + im*co23))
+                    end
+                end
+            end
+        end
+    end
+    length(read(file)) == (prod(lp.iL[1:3])*4*8*9*2) ? nothing : error("File not fully read")
+    close(file)
+end
+
+
 Csw!(dws, U, gp, lp)
 println("--------------------------------------------------")
 
@@ -154,7 +232,7 @@ println("--------------------------------------------------")
 println("Computing correlators")
 
 println("Computing fP...")
-file = open(basename(parsed_args["i"])*".fP","w+")
+file = open("./output/"*basename(parsed_args["c"])*".fP","w+")
 fP = Vector{Complex{Float64}}(undef,lp.iL[4])
 for t in 1:lp.iL[4]
     fP[t] = 0.0
@@ -162,50 +240,58 @@ for t in 1:lp.iL[4]
         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])
     end end end
-    println(file,t,"    ",fP[t])
+    println(file,t,"    ",fP[t]/prod(lp.iL[1:3]))
 end
 flush(file)
 close(file)
 
 println("Computing fA...")
-file = open(basename(parsed_args["i"])*".fA","w+")
+file = open("./output/"*basename(parsed_args["c"])*".fA","w+")
 fA = Vector{Complex{Float64}}(undef,lp.iL[4])
 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],A13[b,r]) + dot(A21[b,r],A23[b,r]) + dot(A31[b,r],A33[b,r]) + dot(A12[b,r],A14[b,r]) + dot(A22[b,r],A24[b,r]) + dot(A32[b,r],A34[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]))
+
     end end end
-    println(file,t,"    ",fA[t])
+    println(file,t,"    ",fA[t]/prod(lp.iL[1:3]))
 end
 flush(file)
 close(file)
 
 println("Computing gP...")
-file = open(basename(parsed_args["i"])*".gP","w+")
+file = open("./output/"*basename(parsed_args["c"])*".gP","w+")
 gP = Vector{Complex{Float64}}(undef,lp.iL[4])
 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])
-    println(file,t,"    ",gP[t])
     end end end
+    println(file,t,"    ",gP[t]/prod(lp.iL[1:3]))
 end
 flush(file)
 close(file)
 
 
 println("Computing gA...")
-file = open(basename(parsed_args["i"])*".gA","w+")
+file = open("./output/"*basename(parsed_args["c"])*".gA","w+")
 gA = Vector{Complex{Float64}}(undef,lp.iL[4])
 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],B13[b,r]) + dot(B21[b,r],B23[b,r]) + dot(B31[b,r],B33[b,r]) + dot(B12[b,r],B14[b,r]) + dot(B22[b,r],B24[b,r]) + dot(B32[b,r],B34[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]))
     end end end
-    println(file,t,"    ",gA[t])
+    println(file,t,"    ",gA[t]/prod(lp.iL[1:3]))
 end
 flush(file)
 close(file)
+
+
+
+print("\n\n")
+print_timer()