From 42e539c9bdec203bd59a13396c38ff5a4adc8f9c Mon Sep 17 00:00:00 2001
From: Alberto Ramos <alberto.ramos@ific.uv.es>
Date: Wed, 13 Oct 2021 18:54:19 +0200
Subject: [PATCH] Added missing files to working tree

---
 src/Scalar/ScalarFields.jl | 31 +++++++++++++
 src/Scalar/ScalarHMC.jl    | 91 ++++++++++++++++++++++++++++++++++++++
 2 files changed, 122 insertions(+)
 create mode 100644 src/Scalar/ScalarFields.jl
 create mode 100644 src/Scalar/ScalarHMC.jl

diff --git a/src/Scalar/ScalarFields.jl b/src/Scalar/ScalarFields.jl
new file mode 100644
index 0000000..f4a7697
--- /dev/null
+++ b/src/Scalar/ScalarFields.jl
@@ -0,0 +1,31 @@
+###
+### "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:    ScalarFields.jl
+### created: Mon Oct 11 21:18:26 2021
+###                               
+
+function randomize!(f, sp::ScalarParm{NS}, lp::SpaceParm, ymws::YMworkspace) where {NS}
+    
+    m = CUDA.randn(ymws.PRC, lp.bsz, 4, NS, lp.rsz)
+    CUDA.@sync begin
+        CUDA.@cuda threads=lp.bsz blocks=lp.rsz krnl_assign_SU2fund!(f,m,sp,lp)
+    end
+
+    return nothing
+end
+
+function krnl_assign_SU2fund!(f::AbstractArray{T}, m, sp::ScalarParm{NS}, lp::SpaceParm) where {T, NS}
+
+    b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
+    for i in 1:NS
+        f[b,i,r] = SU2fund(complex(m[b,1,i,r], m[b,2,i,r]),
+                           complex(m[b,3,i,r], m[b,4,i,r]))
+    end
+
+    return nothing
+end
diff --git a/src/Scalar/ScalarHMC.jl b/src/Scalar/ScalarHMC.jl
new file mode 100644
index 0000000..964f08b
--- /dev/null
+++ b/src/Scalar/ScalarHMC.jl
@@ -0,0 +1,91 @@
+###
+### "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:    ScalarHMC.jl
+### created: Sun Oct 10 23:40:28 2021
+###                               
+
+
+function hamiltonian(mom, U, pmom, Phi, lp, gp, sp, ymws)
+    
+    SG = gauge_action(U, lp, gp, ymws)
+    SS = scalar_action(U, Phi, lp, sp, ymws)
+    PG = CUDA.mapreduce(norm2, +, mom)/2
+    PS = CUDA.mapreduce(norm2, +, pmom)/2
+
+    println("Hamiltonian: ", SG, " ", SS, " ", PG, " ",PS)
+    
+    return SG+SS+PG+PS
+end
+
+function HMC!(U, Phi, eps, ns, lp::SpaceParm, gp::GaugeParm, sp::ScalarParm, ymws::YMworkspace{T}, sws::ScalarWorkspace; noacc=false) where T
+
+    int = omf4(T, eps, ns)
+    ymws.U1 .= U
+    sws.Phi .= Phi
+    
+    randomize!(ymws.mom, lp, ymws)
+    randomize!(sws.mom, sp, lp, ymws)
+    hini = hamiltonian(ymws.mom, U, sws.mom, Phi, lp, gp, sp, ymws)
+    println(hini)
+    
+    MD!(ymws.mom, U, sws.mom, Phi, int, lp, gp, sp, ymws, sws)
+
+    dh   = hamiltonian(ymws.mom, U, sws.mom, Phi, lp, gp, sp, ymws) - hini
+    println(dh+hini)
+    pacc = exp(-dh)
+
+    acc = true
+    if (noacc)
+        return dh, acc
+    end
+    
+    if (pacc < 1.0)
+        r = rand()
+        if (pacc < r) 
+            U   .= ymws.U1
+            Phi .= sws.Phi
+            acc = false
+        end
+    end
+
+    return dh, acc
+end
+
+function MD!(mom, U, pmom, Phi, int::IntrScheme{NI, T}, lp::SpaceParm, gp::GaugeParm{T}, sp::ScalarParm, ymws::YMworkspace{T}, sws::ScalarWorkspace) where {NI, T <: AbstractFloat}
+
+    YM.force_gauge(ymws, U, gp.c0, lp)
+    force_scalar(ymws, sws, U, Phi, sp, gp, lp)
+    
+    mom  .= mom  .+ (int.r[1]*int.eps) .* ymws.frc1
+    pmom .= pmom .+ (int.r[1]*int.eps) .* sws.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])
+            Phi .= Phi .+ (int.eps*int.r[k]).*pmom
+            if k == NI
+                off = -1
+            end
+            k += off
+            
+            YM.force_gauge(ymws, U, gp.c0, lp)
+            force_scalar(ymws, sws, U, Phi, sp, gp, lp)
+            if (i < int.ns) && (k == 1)
+                mom  .= mom  .+ (2*int.r[k]*int.eps) .* ymws.frc1
+                pmom .= pmom .+ (2*int.r[k]*int.eps) .* sws.frc1
+            else
+                mom  .= mom  .+ (int.r[k]*int.eps) .* ymws.frc1
+                pmom .= pmom .+ (int.r[k]*int.eps) .* sws.frc1
+            end
+            k += off
+        end
+    end
+    
+    return nothing
+end