From 0ae90b253437a7bf427f482eed9698cbc358986d Mon Sep 17 00:00:00 2001
From: Alberto Ramos <alberto.ramos@ific.uv.es>
Date: Sat, 9 Oct 2021 00:40:04 +0200
Subject: [PATCH] Added module for integrators

---
 src/LatticeGPU.jl |  5 ++++
 src/MD/MD.jl      | 35 ++++++++++++++++++++++++++++
 src/YM/YM.jl      |  1 +
 src/YM/YMhmc.jl   | 59 ++++++++++++++++++++---------------------------
 4 files changed, 66 insertions(+), 34 deletions(-)
 create mode 100644 src/MD/MD.jl

diff --git a/src/LatticeGPU.jl b/src/LatticeGPU.jl
index 337e5bb..c920257 100644
--- a/src/LatticeGPU.jl
+++ b/src/LatticeGPU.jl
@@ -29,6 +29,11 @@ include("Fields/Fields.jl")
 using .Fields
 export vector_field, scalar_field, nscalar_field
 
+include("MD/MD.jl")
+using .MD
+export IntrScheme
+export omf4, leapfrog
+
 include("YM/YM.jl")
 
 using .YM
diff --git a/src/MD/MD.jl b/src/MD/MD.jl
new file mode 100644
index 0000000..a67d43b
--- /dev/null
+++ b/src/MD/MD.jl
@@ -0,0 +1,35 @@
+###
+### "THE BEER-WARE LICENSE":
+### Alberto Ramos wrote this file. As long as you retain this 
+### notice you can do whatever you want with this stuff. If we meet some 
+### day, and you think this stuff is worth it, you can buy me a beer in 
+### return. <alberto.ramos@cern.ch>
+###
+### file:    MD.jl
+### created: Fri Oct  8 21:53:14 2021
+###                               
+
+module MD
+
+# Dalla Brida / Luscher coefficients of
+# fourth order integrator
+const r1 =  0.08398315262876693
+const r2 =  0.25397851084105950
+const r3 =  0.68223653357190910
+const r4 = -0.03230286765269967
+const r5 =  0.5-r1-r3
+const r6 =  1.0-2.0*(r2+r4)
+
+struct IntrScheme{N, T}
+    r::NTuple{N, T}
+    eps::T
+    ns::Int64    
+end
+
+omf4(::Type{T}, eps, ns) where T = IntrScheme{6,T}((r1,r2,r3,r4,r5,r6), eps, ns)
+leapfrog(::Type{T}, eps, ns) where T = IntrScheme{2,T}((0.5,1.0,0.5), eps, ns)
+
+export IntrScheme, omf4, leapfrog
+
+
+end
diff --git a/src/YM/YM.jl b/src/YM/YM.jl
index f257b70..4fecea5 100644
--- a/src/YM/YM.jl
+++ b/src/YM/YM.jl
@@ -16,6 +16,7 @@ using CUDA, Random, StructArrays
 using ..Space
 using ..Groups
 using ..Fields
+using ..MD
 
 import Base.show
 
diff --git a/src/YM/YMhmc.jl b/src/YM/YMhmc.jl
index 55bfbbc..07579e2 100644
--- a/src/YM/YMhmc.jl
+++ b/src/YM/YMhmc.jl
@@ -48,14 +48,15 @@ function hamiltonian(mom, U, lp, gp, ymws)
     return K+V
 end
 
-function HMC!(U, eps, ns, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace; noacc=false)
+function HMC!(U, eps, ns, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace{T}; noacc=false) where T
 
+    int = omf4(T, eps, ns)
     ymws.U1 .= U
 
     randomize!(ymws.mom, lp, ymws)
     hini = hamiltonian(ymws.mom, U, lp, gp, ymws)
 
-    OMF4!(ymws.mom, U, eps, ns, lp, gp, ymws)
+    MD!(ymws.mom, U, int, lp, gp, ymws)
     
     dh   = hamiltonian(ymws.mom, U, lp, gp, ymws) - hini
     pacc = exp(-dh)
@@ -76,40 +77,30 @@ function HMC!(U, eps, ns, lp::SpaceParm, gp::GaugeParm, ymws::YMworkspace; noacc
     return dh, acc
 end
 
-function OMF4!(mom, U, eps, ns, lp::SpaceParm, gp::GaugeParm{T}, ymws::YMworkspace{T}) where T <: AbstractFloat
+function MD!(mom, U, int::IntrScheme{NI, T}, lp::SpaceParm, gp::GaugeParm{T}, ymws::YMworkspace{T}) where {NI, T <: AbstractFloat}
 
-    r1::T =  0.08398315262876693
-    r2::T =  0.2539785108410595
-    r3::T =  0.6822365335719091
-    r4::T = -0.03230286765269967
-    r5::T =  0.5-r1-r3
-    r6::T =  1.0-2.0*(r2+r4)
-
-    ee = eps*gp.beta/gp.ng
+    ee = int.eps*gp.beta/gp.ng
     force_gauge(ymws, U, gp.c0, lp)
-    for i in 1:ns
-        mom .= mom .+ (r1*ee) .* ymws.frc1
-        U .= expm.(U, mom, eps*r2)
-    
-        force_gauge(ymws, U, gp.c0, lp)
-        mom .= mom .+ (r3*ee) .* ymws.frc1
-        U .= expm.(U, mom, eps*r4)
-
-        force_gauge(ymws, U, gp.c0, lp)
-        mom .= mom .+ (r5*ee) .* ymws.frc1
-        U .= expm.(U, mom, eps*r6)
-
-        force_gauge(ymws, U, gp.c0, lp)
-        mom .= mom .+ (r5*ee) .* ymws.frc1
-        U .= expm.(U, mom, eps*r4)
-
-        force_gauge(ymws, U, gp.c0, lp)
-        mom .= mom .+ (r3*ee) .* ymws.frc1
-        U .= expm.(U, mom, eps*r2)
-
-        force_gauge(ymws, U, gp.c0, lp)
-        mom .= mom .+ (r1*ee) .* ymws.frc1
+    mom .= mom .+ (int.r[1]*ee) .* ymws.frc1
+    for i in 1:int.ns
+        k   = 2
+        off = 1
+        for j in 1:NI-1
+            U .= expm.(U, mom, int.eps*int.r[k])
+            if k == NI
+                off = -1
+            end
+            k += off
+            
+            force_gauge(ymws, U, gp.c0, lp)
+            if (i < int.ns) && (k == 1)
+                mom .= mom .+ (2*int.r[k]*ee) .* ymws.frc1
+            else
+                mom .= mom .+ (int.r[k]*ee) .* ymws.frc1
+            end
+            k += off
+        end
     end
-
+    
     return nothing
 end