From 860f725d9933a335ac2b6fad9562544e065fbf87 Mon Sep 17 00:00:00 2001
From: Alberto Ramos <alberto.ramos@ific.uv.es>
Date: Wed, 3 Nov 2021 15:08:51 +0100
Subject: [PATCH] Added reunitarization after trajectory

---
 src/Groups/GroupSU2.jl |  6 +++---
 src/Groups/GroupSU3.jl | 16 +++++++++++++++-
 src/Groups/Groups.jl   |  2 +-
 src/LatticeGPU.jl      |  2 +-
 src/YM/YMhmc.jl        |  2 ++
 5 files changed, 22 insertions(+), 6 deletions(-)

diff --git a/src/Groups/GroupSU2.jl b/src/Groups/GroupSU2.jl
index 3d98b0a..e45c65a 100644
--- a/src/Groups/GroupSU2.jl
+++ b/src/Groups/GroupSU2.jl
@@ -24,11 +24,11 @@ tr(g::SU2{T})            where T <: AbstractFloat = complex(2.0*real(g.t1), 0.0)
 dev_one(g::SU2{T})       where T <: AbstractFloat = sqrt(( abs2(g.t1 - one(T)) + abs2(g.t2))/2)
 
 """
-    function normalize(a::T) where {T <: Group}
+    function unitarize(a::T) where {T <: Group}
 
-Return a normalized element of the group.
+Return a unitarized element of the group.
 """
-function normalize(a::SU2{T}) where T <: AbstractFloat
+function unitarize(a::SU2{T}) where T <: AbstractFloat
     dr = sqrt(abs2(a.t1) + abs2(a.t2))
     if (dr == 0.0)
         return SU2(0.0)
diff --git a/src/Groups/GroupSU3.jl b/src/Groups/GroupSU3.jl
index 9ac8830..bf7d2f8 100644
--- a/src/Groups/GroupSU3.jl
+++ b/src/Groups/GroupSU3.jl
@@ -14,8 +14,22 @@ dag(a::SU3{T})           where T <: AbstractFloat = inverse(a)
 tr(a::SU3{T})            where T <: AbstractFloat = a.u11+a.u22+conj(a.u11*a.u22 - a.u12*a.u21)
 dev_one(g::SU3{T}) where T <: AbstractFloat = sqrt(( abs2(g.u11 - one(T)) + abs2(g.u12) + abs2(g.u13) + abs2(g.u21) + abs2(g.u22 - one(T)) + abs2(g.u23) )/6)
 
+function unitarize(g::SU3{T}) where T <: AbstractFloat
 
-
+    dv = sqrt(abs2(g.u11)+abs2(g.u12)+abs2(g.u13))
+    gu11 = g.u11/dv
+    gu12 = g.u12/dv
+    gu13 = g.u13/dv
+    
+    z    = g.u21*conj(gu11) + g.u22*conj(gu12) + g.u23*conj(gu13)
+    gu21 = g.u21 - z*gu11
+    gu22 = g.u22 - z*gu12
+    gu23 = g.u23 - z*gu13
+    dv = sqrt(abs2(gu21)+abs2(gu22)+abs2(gu23))
+    
+    return SU3{T}(gu11, gu12, gu13, gu21/dv, gu22/dv, gu23/dv)
+end
+    
 function Base.:*(a::SU3{T},b::SU3{T}) where T <: AbstractFloat
 
     bu31 = conj(b.u12*b.u23 - b.u13*b.u22)
diff --git a/src/Groups/Groups.jl b/src/Groups/Groups.jl
index 044d0e6..856ddce 100644
--- a/src/Groups/Groups.jl
+++ b/src/Groups/Groups.jl
@@ -47,7 +47,7 @@ include("GroupU1.jl")
 export U1, U1alg
 
 
-export dot, expm, exp, dag, normalize, inverse, tr, projalg, norm, norm2, isgroup, alg2mat, dev_one
+export dot, expm, exp, dag, unitarize, inverse, tr, projalg, norm, norm2, isgroup, alg2mat, dev_one
 
 
 end # module
diff --git a/src/LatticeGPU.jl b/src/LatticeGPU.jl
index f8184ab..ad2229d 100644
--- a/src/LatticeGPU.jl
+++ b/src/LatticeGPU.jl
@@ -17,7 +17,7 @@ include("Groups/Groups.jl")
 using .Groups
 export Group, Algebra
 export SU2, SU2alg, SU3, SU3alg, M3x3, M2x2, U1, U1alg
-export dot, expm, exp, dag, normalize, inverse, tr, projalg, norm, norm2, isgroup, alg2mat, dev_one
+export dot, expm, exp, dag, unitarize, inverse, tr, projalg, norm, norm2, isgroup, alg2mat, dev_one
 
 include("Space/Space.jl")
 
diff --git a/src/YM/YMhmc.jl b/src/YM/YMhmc.jl
index 1a3dd8d..8d3fb45 100644
--- a/src/YM/YMhmc.jl
+++ b/src/YM/YMhmc.jl
@@ -84,6 +84,8 @@ function HMC!(U, int::IntrScheme, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspac
                 acc = false
             end
         end
+
+        U .= unitarize.(U)
         
     end
     return dh, acc