Added module for integrators

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

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

src/YM/YM.jl

@@ -16,6 +16,7 @@ using CUDA, Random, StructArrays
 using ..Space
 using ..Groups
 using ..Fields
+using ..MD
 
 import Base.show
 

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,39 +77,29 @@ 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
+    ee = int.eps*gp.beta/gp.ng
-    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
     force_gauge(ymws, U, gp.c0, lp)
-    for i in 1:ns
+    mom .= mom .+ (int.r[1]*ee) .* ymws.frc1
-        mom .= mom .+ (r1*ee) .* ymws.frc1
+    for i in 1:int.ns
-        U .= expm.(U, mom, eps*r2)
+        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)
-        mom .= mom .+ (r3*ee) .* ymws.frc1
+            if (i < int.ns) && (k == 1)
-        U .= expm.(U, mom, eps*r4)
+                mom .= mom .+ (2*int.r[k]*ee) .* ymws.frc1
-
+            else
-        force_gauge(ymws, U, gp.c0, lp)
+                mom .= mom .+ (int.r[k]*ee) .* ymws.frc1
-        mom .= mom .+ (r5*ee) .* ymws.frc1
+            end
-        U .= expm.(U, mom, eps*r6)
+            k += off
-
+        end
-        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
     end
     
     return nothing