Added routines for the computation of the MD force

This commit is contained in:
Alberto Ramos 2021-10-06 17:35:52 +02:00
parent 97e46a538c
commit 55762f12d3
2 changed files with 58 additions and 5 deletions

View file

@ -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)

52
src/Scalar/ScalarForce.jl Normal file
View file

@ -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