diff --git a/src/main/test_kapp.jl b/src/main/test_kapp.jl
new file mode 100644
index 0000000..c5071af
--- /dev/null
+++ b/src/main/test_kapp.jl
@@ -0,0 +1,115 @@
+using CUDA, Logging, StructArrays, Random, TimerOutputs, ADerrors
+
+CUDA.allowscalar(true)
+import Pkg
+Pkg.activate("/lhome/ific/a/alramos/s.images/julia/workspace/LatticeGPU")
+#Pkg.activate("/home/alberto/code/julia/LatticeGPU")
+using LatticeGPU
+
+lp = SpaceParm{4}((16,16,16,16), (4,4,4,4))
+gp = GaugeParm(2.25, 1.0, (0.0,0.0), 2)
+
+NSC = 1
+println("Space  Parameters: ", lp)
+println("Gauge  Parameters: ", gp)
+GRP  = SU2
+ALG  = SU2alg
+SCL  = SU2fund
+PREC = Float64
+println("Precision:         ", PREC)
+
+println("Allocating YM workspace")
+ymws = YMworkspace(GRP, PREC, lp)
+println("Allocating Scalar workspace")
+sws  = ScalarWorkspace(PREC, NSC, lp)
+
+# Seed RNG
+println("Seeding CURAND...")
+Random.seed!(CURAND.default_rng(), 1234)
+Random.seed!(1234)
+
+# Main program
+println("Allocating gauge field")
+U = vector_field(GRP{PREC}, lp)
+fill!(U, one(GRP{PREC}))
+println("Allocating scalar field")
+Phi = nscalar_field(SCL{PREC}, NSC, lp)
+fill!(Phi, zero(SCL{PREC}))
+
+dt  = 0.15
+ns  = 12
+
+nth    = 500
+nms    = 5000
+nsteps = 30
+h      = 0.6/nsteps
+pl     = Array{Float64, 2}(undef, nms+nth, nsteps)
+rho    = Array{Float64, 2}(undef, nms+nth, nsteps)
+dh     = Array{Float64, 2}(undef, nms+nth, nsteps)
+acc    = Array{Bool, 2}(undef, nms+nth, nsteps)
+
+
+for i in 1:nsteps
+    sp = ScalarParm((h*(i-1),), (0.5,))
+    println("## Simulating Scalar parameters: ")
+    println(sp)
+    
+    k = 0
+    for j in 1:nth
+        k = k + 1
+        dh[k,i], acc[k,i] = HMC!(U,Phi, dt,ns,lp, gp, sp, ymws, sws)
+        pl[k,i]  = plaquette(U,lp, gp, ymws)
+        rho[k,i] = CUDA.mapreduce(norm2, +, Phi)/prod(lp.iL)
+        
+        println("  THM $j/$nth (kappa: $(sp.kap[1])):   ", acc[k,i], " ",
+                dh[k,i], " ", pl[k,i], " ", rho[k,i])
+    end
+    println(" ")
+    
+    for j in 1:nms
+        k = k + 1
+        dh[k,i], acc[k,i] = HMC!(U,Phi, dt,ns,lp, gp, sp, ymws, sws)
+        pl[k,i]  = plaquette(U,lp, gp, ymws)
+        rho[k,i] = CUDA.mapreduce(norm2, +, Phi)/prod(lp.iL)
+        
+        println("  MSM $j/$nth (kappa: $(sp.kap[1])):   ", acc[k,i], " ",
+                dh[k,i], " ", pl[k,i], " ", rho[k,i])
+    end
+
+    println("\n\n")
+end
+
+println("## Summary of results")
+for i in 1:nsteps
+    kap = h*(i-1)
+    println(" # kappa: ", kap)
+
+    println("   - Acceptance: ", count(acc[nth+1:end, i])/nms)
+
+    uw = uwreal(exp.(.-dh[nth+1:end, i]), "Foo")
+    try 
+        uwerr(uw)
+        println("   - I am one:   ", uw, "    (tauint: ", taui(uw,"Foo"), ")")
+    catch e
+        println("No Error ")
+        println("   - I am one:   ", value(uw))
+    end
+        
+    uw = uwreal(pl[nth+1:end, i], "Foo")
+    uwerr(uw)
+    println("   - Plaquette:  ", uw, "    (tauint: ", taui(uw,"Foo"), ")")
+
+    uw = uwreal(rho[nth+1:end, i], "Foo")
+    uwerr(uw)
+    println("   - Field mod:  ", uw, "    (tauint: ", taui(uw,"Foo"), ")")
+
+    println("")
+end
+
+
+println("## Timming results")
+print_timer(linechars = :ascii)
+
+
+
+