Added import for CERN format. Corrected bug in Qtop

src/LatticeGPU.jl

@@ -45,6 +45,6 @@ export Eoft_clover, Eoft_plaq, Qtop
 export FlowIntr, wfl_euler, zfl_euler, wfl_rk2, zfl_rk2, wfl_rk3, zfl_rk3
 export flw, flw_adapt
 export sfcoupling, bndfield, setbndfield
-export import_lex64
+export import_lex64, import_cern64
 
 end # module

src/Space/Space.jl

@@ -317,7 +317,22 @@ end
 
     return (b,r)
 end
-    
+
+"""
+    function point_color(p::NTuple{2,Int64}, lp::SpaceParm)
+
+Returns the sum of the cartesian coordinates of the point p=(b,r).
+"""
+@inline function point_color(p::NTuple{2,Int64}, lp::SpaceParm)
+
+    s = cnt(p[1], p[2], 1, lp)
+    for i in 2:lp.ndim
+        s = s + cnt(p[1], p[2], i, lp)
+    end
+
+    return s
+end
+
 export up, dw, updw, point_index, point_coord, point_time
 export BC_PERIODIC, BC_OPEN, BC_SF_AFWB, BC_SF_ORBI
 

src/YM/YM.jl

@@ -147,6 +147,6 @@ include("YMsf.jl")
 export sfcoupling, bndfield, setbndfield
 
 include("YMio.jl")
-export import_lex64
+export import_lex64, import_cern64
 
 end

src/YM/YMflow.jl

@@ -296,21 +296,21 @@ function Qtop(Qslc, U, gp::GaugeParm, lp::SpaceParm{4,M,B,D}, ymws::YMworkspace)
         
         fill!(ymws.rm, zero(eltype(ymws.rm)))
         CUDA.@sync begin
-            CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_field_tensor!(ymws.frc1, ymws.frc2, U, gp.Ubnd, 1,6, ztw[1], ztw[6], lp)
+            CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_field_tensor!(ymws.frc1, ymws.frc2, U, gp.Ubnd, 1,5, ztw[1], ztw[5], lp)
         end
         CUDA.@sync begin
             CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_add_qd!(ymws.rm, -, ymws.frc1, ymws.frc2, U, lp)
         end
     
         CUDA.@sync begin
-            CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_field_tensor!(ymws.frc1, ymws.frc2, U, gp.Ubnd, 2,5, ztw[2], ztw[5], lp)
+            CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_field_tensor!(ymws.frc1, ymws.frc2, U, gp.Ubnd, 2,4, ztw[2], ztw[4], lp)
         end
         CUDA.@sync begin
             CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_add_qd!(ymws.rm, +, ymws.frc1, ymws.frc2, U, lp)
         end
     
         CUDA.@sync begin
-            CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_field_tensor!(ymws.frc1, ymws.frc2, U, gp.Ubnd, 3,4, ztw[3], ztw[4], lp)
+            CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_field_tensor!(ymws.frc1, ymws.frc2, U, gp.Ubnd, 3,6, ztw[3], ztw[6], lp)
         end
         CUDA.@sync begin
             CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_add_qd!(ymws.rm, -, ymws.frc1, ymws.frc2, U, lp)
@@ -491,10 +491,10 @@ function krnl_field_tensor!(frc1::AbstractArray{TA}, frc2, U::AbstractArray{T},
         frc2[bd,3,rd]   = zero(TA)
         frc2[bu2,4,ru2] = projalg(l2*l1)
     else
-        frc2[b,1,r]     = projalg(ztw1, Ush[b,1]*l1/Ush[b,2])
+        frc2[b,1,r]     = projalg(ztw2, Ush[b,1]*l1/Ush[b,2])
-        frc2[bu1,2,ru1] = projalg(ztw1, l1*l2)
+        frc2[bu1,2,ru1] = projalg(ztw2, l1*l2)
-        frc2[bd,3,rd]   = projalg(ztw1, gt2\(l2*gt1))
+        frc2[bd,3,rd]   = projalg(ztw2, gt2\(l2*gt1))
-        frc2[bu2,4,ru2] = projalg(ztw1, l2*l1)
+        frc2[bu2,4,ru2] = projalg(ztw2, l2*l1)
     end
 
     return nothing

src/YM/YMio.jl

@@ -45,3 +45,47 @@ function import_lex64(fname, lp::SpaceParm)
 
     return CuArray(Ucpu), Ubnd
 end
+
+"""
+    function import_cern64(fname::String, ibc, lp::SpaceParm)
+
+import a double precision configuration in CERN format. 
+"""
+function import_cern64(fname, ibc, lp::SpaceParm; log=true)
+
+    fp = open(fname, "r")
+    iL = Vector{Int32}(undef, 4)
+    read!(fp, iL)
+    avgpl = Vector{Float64}(undef, 1)
+    read!(fp, avgpl)
+    if log
+        println("# [import_cern64] Read from conf file: ", iL, " (plaq: ", avgpl, ")")
+    end
+    
+    dtr = [4,1,2,3]
+    assign(V, ic) = SU3{Float64}(V[1,ic],V[2,ic],V[3,ic],V[4,ic],V[5,ic],V[6,ic])
+    
+    Ucpu = Array{SU3{Float64}, 3}(undef, lp.bsz, lp.ndim, lp.rsz)
+    V = Array{ComplexF64, 2}(undef, 9, 2)
+    for i4 in 1:lp.iL[4]
+        for i1 in 1:lp.iL[1]
+            for i2 in 1:lp.iL[2]
+                for i3 in 1:lp.iL[3]
+                    if (mod(i1+i2+i3+i4-4, 2) == 1)
+                        b, r = point_index(CartesianIndex(i1,i2,i3,i4), lp)
+                        for id in 1:lp.ndim
+                            read!(fp, V)
+                            Ucpu[b,dtr[id],r] = assign(V, 1)
+
+                            bd, rd = dw((b,r), dtr[id], lp)
+                            Ucpu[bd,dtr[id],rd] = assign(V, 2)
+                        end
+                    end
+                end
+            end
+        end
+    end
+    close(fp)
+
+    return CuArray(Ucpu)
+end