From 9ba07b506aa4d8919485692be266c81bb5eefd36 Mon Sep 17 00:00:00 2001
From: Nicolas Lang <nicolas.lang@ific.uv.es>
Date: Mon, 26 May 2025 12:48:04 +0200
Subject: [PATCH] fixed the OBC branch in the Zeuthen term - now results match
 openQCD

---
 src/YM/YM.jl     |  2 +-
 src/YM/YMflow.jl | 30 ++++++++++++++++++++++++++----
 2 files changed, 27 insertions(+), 5 deletions(-)

diff --git a/src/YM/YM.jl b/src/YM/YM.jl
index 5578730..8f862a4 100644
--- a/src/YM/YM.jl
+++ b/src/YM/YM.jl
@@ -165,7 +165,7 @@ include("YMfields.jl")
 export randomize!, zero!, norm2
 
 include("YMact.jl")
-export krnl_plaq!, force_gauge, force_gauge_flw, force_wilson
+export krnl_plaq!, force_gauge, force_gauge_flw, force_wilson, bnd_rescale_flw!
 
 include("YMhmc.jl")
 export gauge_action, hamiltonian, plaquette, HMC!, MD!
diff --git a/src/YM/YMflow.jl b/src/YM/YMflow.jl
index 92c8e2c..62abde1 100644
--- a/src/YM/YMflow.jl
+++ b/src/YM/YMflow.jl
@@ -127,6 +127,22 @@ function add_zth_term(ymws::YMworkspace, U, lp)
     return nothing
 end
 
+
+"""
+    function rescale_bnd(ymws::YMworkspace, lp)
+
+Rescales the time-boundary of ymws.frc1 by a factor of 2.0 (needed for open BC)
+"""
+function rescale_bnd(ymws::YMworkspace, lp)
+
+    CUDA.@sync begin
+        CUDA.@cuda threads=lp.bsz blocks=lp.rsz bnd_rescale_flw!(ymws.frc1,lp::SpaceParm)
+    end
+
+    return nothing
+end
+
+
 function krnl_add_zth!(frc, frc2::AbstractArray{TA}, U::AbstractArray{TG}, lp::SpaceParm{N,M,B,D}) where {TA,TG,N,M,B,D}
 
     @inbounds begin 
@@ -155,12 +171,14 @@ function krnl_add_zth!(frc, frc2::AbstractArray{TA}, U::AbstractArray{TG}, lp::S
                 end
             end
             if OBC
-                if (it > 1) && (it < lp.iL[end])
+		if (id < N)
                     frc2[b,id,r] = (5/6)*frc[b,id,r] + (1/6)*(projalg(Ud\Y*Ud) +
                         projalg(U[b,id,r]*X/U[b,id,r]))
-                elseif ((it == lp.iL[end]) || (it == 1))  && (id < N)
+                elseif (it > 1) && (it < (lp.iL[end]-1))
                     frc2[b,id,r] = (5/6)*frc[b,id,r] + (1/6)*(projalg(Ud\Y*Ud) +
                         projalg(U[b,id,r]*X/U[b,id,r]))
+		elseif (it == 1) || (it == (lp.iL[end]-1))
+		    frc2[b,id,r] = frc[b,id,r]
                 end
             else
                 frc2[b,id,r] = (5/6)*frc[b,id,r] + (1/6)*(projalg(Ud\Y*Ud) +
@@ -171,6 +189,8 @@ function krnl_add_zth!(frc, frc2::AbstractArray{TA}, U::AbstractArray{TG}, lp::S
     return nothing
 end
 
+
+
 """
         function flw(U, int::FlowIntr{NI,T}, ns::Int64, gp::GaugeParm, lp::SpaceParm, ymws::YMworkspace)
 
@@ -204,18 +224,20 @@ flw(U, int::FlowIntr{NI,T}, ns::Int64, gp::GaugeParm, lp::SpaceParm, ymws::YMwor
 function flw(U, int::FlowIntr{NI,T}, ns::Int64, eps, gp::GaugeParm, lp::SpaceParm{N,M,BC_OPEN,D}, ymws::YMworkspace) where {NI,T,N,M,D}
     @timeit "Integrating flow equations" begin
         for i in 1:ns
-            force_gauge_flw(ymws, U, int.c0, 1, gp, lp)
+            force_gauge(ymws, U, int.c0, 1, gp, lp)
             if int.add_zth
                 add_zth_term(ymws::YMworkspace, U, lp)
             end
+            rescale_bnd(ymws, lp)
             ymws.mom .= ymws.frc1
             U .= expm.(U, ymws.mom, 2*eps*int.r)
 
             for k in 1:NI
-                force_gauge_flw(ymws, U, int.c0, 1, gp, lp)
+                force_gauge(ymws, U, int.c0, 1, gp, lp)
                 if int.add_zth
                     add_zth_term(ymws::YMworkspace, U, lp)
                 end
+                rescale_bnd(ymws, lp)
                 ymws.mom .= int.e0[k].*ymws.mom .+ int.e1[k].*ymws.frc1
                 U .= expm.(U, ymws.mom, 2*eps)
             end