From 55762f12d3175cf2e317538cbc6d50c9597a3752 Mon Sep 17 00:00:00 2001
From: Alberto Ramos <alberto.ramos@ific.uv.es>
Date: Wed, 6 Oct 2021 17:35:52 +0200
Subject: [PATCH] Added routines for the computation of the MD force

---
 src/Groups/FundamentalSU2.jl | 11 ++++----
 src/Scalar/ScalarForce.jl    | 52 ++++++++++++++++++++++++++++++++++++
 2 files changed, 58 insertions(+), 5 deletions(-)
 create mode 100644 src/Scalar/ScalarForce.jl

diff --git a/src/Groups/FundamentalSU2.jl b/src/Groups/FundamentalSU2.jl
index e3fb8b8..74de503 100644
--- a/src/Groups/FundamentalSU2.jl
+++ b/src/Groups/FundamentalSU2.jl
@@ -20,10 +20,11 @@ Base.:*(a::SU2{T},b::SU2fund{T}) where T <: AbstractFloat = SU2fund{T}(a.t1*b.t1
 Base.:/(a::SU2fund{T},b::SU2{T}) where T <: AbstractFloat = SU2fund{T}(a.t1*conj(b.t1)+a.t2*conj(b.t2),-a.t1*b.t2+a.t2*b.t1)
 Base.:\(a::SU2{T},b::SU2fund{T}) where T <: AbstractFloat = SU2fund{T}(conj(a.t1)*b.t1+a.t2*conj(b.t2),conj(a.t1)*b.t2-a.t2*conj(b.t1))
 
-"""
-    \(phi::SU2fund{T},u::SU2{T})
+Base.:+(a::SU2fund{T},b::SU2fund{T}) where T <: AbstractFloat = SU2fund{T}(a.t1+b.t1,a.t2+b.t2)
+Base.:-(a::SU2fund{T},b::SU2fund{T}) where T <: AbstractFloat = SU2fund{T}(a.t1-b.t1,a.t2-b.t2)
 
-This operation has to be understood as (phi^+ u) 
-"""
-Base.:\(a::SU2fund{T},b::SU2{T}) where T <: AbstractFloat = SU2fund{T}(conj(a.t1)*b.t1+a.t2*conj(b.t2),conj(a.t1)*b.t2-a.t2*conj(b.t1))
+# Operations with numbers
+Base.:*(a::SU2fund{T},b::Number) where T <: AbstractFloat = SU2fund{T}(b*a.t1,b*a.t2)
+Base.:*(b::Number,a::SU2fund{T}) where T <: AbstractFloat = SU2fund{T}(b*a.t1,b*a.t2)
+Base.:/(a::SU2fund{T},b::Number) where T <: AbstractFloat = SU2fund{T}(a.t1/b,a.t2/b)
 
diff --git a/src/Scalar/ScalarForce.jl b/src/Scalar/ScalarForce.jl
new file mode 100644
index 0000000..e1cbe63
--- /dev/null
+++ b/src/Scalar/ScalarForce.jl
@@ -0,0 +1,52 @@
+###
+### "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:    ScalarForce.jl
+### created: Wed Oct  6 15:39:07 2021
+###                               
+
+function krnl_force_scalar!(fgauge, fscalar, U::AbstractArray{TG}, Phi::AbstractArray{TS}, sp::ScalarParam{NP,T}, lp::SpaceParm{N,M,D}) where {TG,TS,NP,T,N,M,D}
+
+    b, r = CUDA.threadIdx().x, CUDA.blockIdx().x
+
+    Ush = @cuStaticSharedMem(TG, (D,N))
+    Psh = @cuStaticSharedMem(TS, (D,NP))
+    
+    for id in 1:N
+        Ush[b,id] = U[b,id,r]
+    end
+    for i in 1:NP
+        Psh[b,i] = Phi[b,i,r]
+    end
+    sync_threads()
+
+    for n in 1:NP
+        fscalar[b,n,u] = zero(TS)
+        for id in 1:N
+            bu, ru = up((b,r), id, lp)
+            bd, rd = dw((b,r), id, lp)
+
+            if ru == r
+                p1 = Ush[b,id]*Psh[bu,n]
+            else
+                p1 = Ush[b,id]*Phi[bu,n,ru]
+            end
+
+            if rd == r
+                fscalar[b,n,u] += (-2*sp.kap[n])*(p1 + dag(Psh[bd,n])*Ush[bd,id])
+            else
+                fscalar[b,n,u] += (-2*sp.kap[n])*(p1 + dag(Phi[bd,n,rd])*U[bd,id,rd])
+            end
+            
+            fgauge[b,id,r] += (2*sp.kap[n])*projalg(p1*dag(Psh[b,n]))
+        end
+        fscalar[b,n,r] += ( 2 + 2*sp.eta[n]*(dot(Psh[b,n],Psh[b,n])-1) ) * Psh[b,n]
+    end
+
+
+    return nothing
+end