diff --git a/src/bench/benchStagProp.nim b/src/bench/benchStagProp.nim
index 9580ffc..2ac4f8f 100644
--- a/src/bench/benchStagProp.nim
+++ b/src/bench/benchStagProp.nim
@@ -18,8 +18,11 @@ var v1 = lo.ColorVector()
 var v2 = lo.ColorVector()
 var r = lo.ColorVector()
 var rs = newRNGField(RngMilc6, lo, intParam("seed", 987654321).uint64)
+let warm0 = 0.29 + 500.0/lo.physVol.float
+var warm = floatParam("warm", warm0)
 threads:
-  g.random rs
+  #g.random rs
+  g.warm warm, rs
   g.setBC
   g.stagPhase
   v1 := 0
@@ -67,6 +70,7 @@ for i in 0..3:
   g3[2*i] = g[i]
   g3[2*i+1] = lo.ColorMatrix()
   g3[2*i+1].randomSU rs
+  g3[2*i+1] *= 0.1
 var s3 = newStag3(g3)
 #s3.D(v2, v1, m)
 s3.solve(v2, v1, mass, sp)
diff --git a/src/comms/halo.nim b/src/comms/halo.nim
new file mode 100644
index 0000000..ba2cc95
--- /dev/null
+++ b/src/comms/halo.nim
@@ -0,0 +1,441 @@
+import qex, base/hyper, comms/gather
+getOptimPragmas()
+
+type
+  HaloLayout*[L] = ref object
+    lo*: L  # layout
+    outerExt*: seq[int32]  # extended outer lattice size
+    offset*: seq[int32]  # offset of outer lattice in extended outer
+    lex*: seq[int32]  # outerExt lex index of extended sites
+    index*: seq[int32]  # index of extended site for given lex index
+    neighborFwd*: seq[seq[int32]]  # fwd neighbor for extended outer lattice
+    neighborBck*: seq[seq[int32]]  # bck neighbor for extended outer lattice
+    nOut*: int  # sites in outer lattice
+    nExt*: int  # sites in extended outer lattice
+  Halo*[L,F,T] = ref object
+    layout*: HaloLayout[L]
+    field*: F
+    halo*: alignedMem[T]
+    nOut*: int  # sites in outer lattice
+    nExt*: int  # sites in extended outer lattice
+  HaloMap*[L] = ref object
+    layout*: HaloLayout[L]
+    offsets*: seq[seq[int32]]  # offsets of outer sites needed
+    gather*: GatherMap  # gather map to fetch needed halo sites
+    gatherRev*: GatherMap  # gather map to send halo sites back
+  GatherHalo*[F,T;Rev:static bool] = object
+    src: F
+    dest: ptr UncheckedArray[T]
+    elemsize: int
+    vlen: int
+
+proc norm2*(h: Halo): auto =
+  var s = h.halo[0].norm2
+  for i in 1..<h.halo.len:
+    s += h.halo[i].norm2
+  var t = s.simdReduce
+  h.field.l.threadRankSum(t)
+  t += h.field.norm2
+  t
+
+proc sum*(h: Halo): auto =
+  var s = h.halo[0]
+  for i in 1..<h.halo.len:
+    s += h.halo[i]
+  var t = s.simdSum
+  h.field.l.threadRankSum(t)
+  t += h.field.sum
+  t
+
+proc inside[T,U,V:openArray[auto]](x: T, lo: U, hi: V): bool {.inline.} =
+  result = true
+  for i in 0..<x.len:
+    if x[i] < lo[i] or x[i] >= hi[i]:
+      result = false
+      break
+
+proc init[T](s: var seq[T], x: openArray[SomeNumber]) =
+  s.newSeq(x.len)
+  for i in 0..<x.len: s[i] = T x[i]
+
+# TODO: order halo by layers
+proc makeHaloLayout*[L:Layout](lo: L, fwdOffset,bckOffset: openarray[SomeInteger]): HaloLayout[L] =
+  tic("makeHaloLayout")
+  let outerHigh = lo.outerGeom + bckOffset
+  let outerExt = outerHigh + fwdOffset
+  let nExt = outerExt.product
+  let nOut = lo.nSitesOuter
+  var lex = newSeq[int32](nExt)
+  var index = newSeq[int32](nExt)
+  let nd = lo.nDim
+  var x = newSeq[int32](nd)
+  var k = int32 nOut
+  for i in 0..<nExt:  # loop over lex indices
+    x.lexCoord(i, outerExt)
+    if x.inside(bckOffset, outerHigh):
+      var y = x - bckOffset
+      let loc = lo.localIndex(y) div L.V
+      lex[loc] = int32 i
+      index[i] = int32 loc
+    else:
+      lex[k] = int32 i
+      index[i] = int32 k
+      inc k
+  toc("lex/index")
+  result.new
+  result.lo = lo
+  result.outerExt.init outerExt
+  result.offset.init bckOffset
+  result.lex = lex
+  result.index = index
+  result.nOut = nOut
+  result.nExt = nExt
+  result.neighborFwd = newSeq[seq[int32]](nd)
+  result.neighborBck = newSeq[seq[int32]](nd)
+  for mu in 0..<nd:
+    result.neighborFwd[mu] = newSeq[int32](nExt)
+    result.neighborBck[mu] = newSeq[int32](nExt)
+  for i in 0..<nExt:  # loop over lex indices
+    let ii = index[i]
+    x.lexCoord(i, outerExt)
+    for mu in 0..<nd:
+      if x[mu] == 0:
+        result.neighborBck[mu][ii] = -1
+      x[mu] += 1
+      if x[mu] < outerExt[mu]:
+        let j = x.lexIndex(outerExt)
+        let ij = index[j]
+        result.neighborFwd[mu][ii] = ij
+        result.neighborBck[mu][ij] = ii
+      else:
+        #echo x
+        result.neighborFwd[mu][ii] = -1
+      x[mu] -= 1
+  toc("neighbors")
+
+proc makeHaloMap*[L](hl: HaloLayout[L], c: Comm, offsets: seq[seq[int32]]): HaloMap[L] =
+  tic("makeHaloMap")
+  result.new
+  result.layout = hl
+  result.offsets = offsets
+  let outerHigh = hl.offset + hl.lo.outerGeom
+  let nd = hl.lo.nDim
+  var x = newSeq[int32](nd)
+  var rl = newSeq[RecvList](0)
+  var sl = newSeq[SendList](0)
+  let vlen = L.V
+  var vecOffset = newSeq[seq[int32]](vlen)
+  for k in 0 ..< vlen:
+    vecOffset[k] = newSeq[int32](nd)
+    let li0 = hl.offset.lexIndex(hl.outerExt)
+    let i0 = hl.index[li0]
+    let ik = vlen*i0 + k
+    hl.lo.coord(x, ik)
+    vecOffset[k] = x - hl.offset
+  for i in hl.nOut ..< hl.nExt:  # loop over halo sites
+    let li = hl.lex[i]
+    x.lexCoord(li, hl.outerExt)
+    var keep = false
+    for j in 0 ..< offsets.len:
+      let y = x - offsets[j]
+      if y.inside(hl.offset, outerHigh):
+        keep = true
+        break
+    if keep:
+      for k in 0 ..< vlen:
+        var y = x + vecOffset[k]
+        for l in 0..<nd:
+          if y[l] < 0: y[l] += int32 hl.lo.physGeom[l]
+          if y[l] >= hl.lo.physGeom[l]: y[l] -= int32 hl.lo.physGeom[l]
+        let ri = hl.lo.rankIndex(y)
+        let didx = vlen*(i-hl.nOut)+k
+        rl.add RecvList(didx: int32 didx, srank: int32 ri.rank, sidx: int32 ri.index)
+        sl.add SendList(sidx: int32 didx, drank: int32 ri.rank, didx: int32 ri.index)
+  toc("RecvList")
+  result.gather = c.makeGatherMap(rl)
+  result.gatherRev = c.makeGatherMap(sl)
+  toc("makeGatherMap")
+proc makeHaloMap*[L,N](hl: HaloLayout[L], c: Comm, offsets: seq[array[N,SomeInteger]]): HaloMap[L] =
+  var o = newSeq[seq[int32]](offsets.len)
+  for i in 0..<offsets.len:
+    let t = offsets[i]
+    o[i] = @t
+  makeHaloMap(hl, c, o)
+
+# T is vectorized type
+proc makeHalo*[L,F,T](hl: HaloLayout[L], f: F, t: typedesc[T]): Halo[L,F,T] =
+  result.new
+  result.layout = hl
+  result.field = f
+  result.nOut = hl.nOut
+  result.nExt = hl.nExt
+  let nhalo = result.nExt - result.nOut
+  result.halo.newU(nhalo)
+
+template makeHalo*[L,F](hl: HaloLayout[L], f: F): auto =
+  makeHalo(hl, f, eval(F.type[0]))
+
+template copy*[F,T;Rev:static bool](gh: GatherHalo[F,T,Rev], d: pointer, s: SomeInteger) =
+  type E = eval(index(type T, type asSimd(0)))
+  let p = cast[ptr E](d)
+  when Rev:
+    let o = s div gh.vlen
+    let i = s mod gh.vlen
+    p[] := gh.dest[o][asSimd(i)]
+  else:
+    p[] := gh.src{s}
+template copy*[F,T;Rev:static bool](gh: GatherHalo[F,T,Rev], d: SomeInteger, s: pointer) =
+  type E = eval(index(type T,type asSimd(0)))
+  let p = cast[ptr E](s)
+  when Rev:
+    #threadCritical:
+    gh.src{d} += p[]
+  else:
+    let o = d div gh.vlen
+    let i = d mod gh.vlen
+    gh.dest[o][asSimd(i)] = p[]
+template copy*[F,T;Rev:static bool](gh: GatherHalo[F,T,Rev], d: SomeInteger, s: SomeInteger) =
+  when Rev:
+    let o = s div gh.vlen
+    let i = s mod gh.vlen
+    #threadCritical:
+    gh.src{d} += gh.dest[o][asSimd(i)]
+  else:
+    let o = d div gh.vlen
+    let i = d mod gh.vlen
+    gh.dest[o][asSimd(i)] = gh.src{s}
+#[
+template copy*[F,T](gh: GatherHalo[F,T], dv: SomeInteger, di: array,
+                    s: array, n: SomeInteger) =
+  echo n, " ", dv, " ", di, " ", s
+  #for i in 0..<n:
+  #  gh.dest[dv][asSimd(di[i])] = gh.src{s[i]}
+  let t = gh.dest[dv]
+  for i in 0..<n:
+    t[asSimd(di[i])] = gh.src{s[i]}
+  gh.dest[dv] = t
+]#
+
+proc update*[L,F,T](h: Halo[L,F,T], hm: HaloMap[L], c: Comm) =
+  tic("Halo update")
+  let elemSize = sizeof(T) div L.V
+  var gh: GatherHalo[F,T,false]
+  gh.src = h.field
+  gh.dest = cast[ptr UncheckedArray[T]](unsafeAddr h.halo[0])
+  gh.elemsize = elemSize
+  gh.vlen = L.V
+  c.gather(hm.gather, gh)
+  toc("gather")
+
+proc updateRev*[L,F,T](h: Halo[L,F,T], hm: HaloMap[L], c: Comm) =
+  tic("Halo update")
+  let elemSize = sizeof(T) div L.V
+  var gh: GatherHalo[F,T,true]
+  gh.src = h.field
+  gh.dest = cast[ptr UncheckedArray[T]](unsafeAddr h.halo[0])
+  gh.elemsize = elemSize
+  gh.vlen = L.V
+  #let g = hm.gather.reverse
+  c.gather(hm.gatherRev, gh)
+  toc("gather")
+
+#[
+proc `[]`*(h: Halo, i: SomeInteger): auto {.inline,noInit.} =
+  let k = i - h.nOut
+  if k < 0:
+    h.field[i]
+  else:
+    h.halo[k]
+
+proc `[]`*(h: var Halo, i: SomeInteger): var auto {.inline.} =
+  var t {.noInit.}: ptr type(h.field[i])
+  let k = i - h.nOut
+  if k < 0:
+    t = addr h.field[i]
+  else:
+    t = addr h.halo[k]
+  t[]
+template `[]`*(h: Halo, i: SomeInteger): auto =
+  var t {.noInit.}: ptr type(h.field[i])
+  let k = i - h.nOut
+  if k < 0:
+    t = addr h.field[i]
+  else:
+    t = addr h.halo[k]
+  t[]
+]#
+
+proc indexPtr*[L,F,T](h: Halo[L,F,T], i: SomeInteger): ptr T {.alwaysInline.} =
+  let k = i - h.nOut
+  result = if k<0: addr h.field[i] else: addr h.halo[k]
+template `[]`*(h: Halo, i: SomeInteger): auto = indexPtr(h,i)[]
+
+proc `[]=`*(h: Halo, i: SomeInteger, x: auto) {.alwaysInline.} =
+  let k = i - h.nOut
+  if k < 0:
+    h.field[i] = x
+  else:
+    h.halo[k] := x
+
+proc neighbor*(h: Halo, i: SomeInteger, mu: SomeInteger, fb: SomeInteger): int32 =
+  if fb > 0:
+    h.map.neighborFwd[mu][i]
+  else:
+    h.map.neighborBck[mu][i]
+
+when isMainModule:
+  qexInit()
+  tic("main")
+  var defaultLat = @[4,4,4,4]
+  defaultSetup()
+  let nd = lo.nDim
+  var seed = 987654321'u
+  #var rng = newRngField(lo, RngMilc6, seed)
+  var rng = newRngField(lo, MRG32k3a, seed)
+  g.gaussian rng
+  toc "gaussian"
+  #var r0 = lo.Real()
+  #var cv1 = lo.ColorVector()
+  #var cv2 = lo.ColorVector()
+  echo 6.0 * g.plaq
+  toc "plaq"
+  #cv0.gaussian rng
+  resetTimers()
+
+  proc testPlaq =
+    tic("testPlaq")
+    #let hl = lo.makeHaloLayout([1,1,1,1],[1,1,1,1])
+    let hl = lo.makeHaloLayout([1,1,1,1],[0,0,0,0])
+    toc "makeHaloLayout"
+    type HM = HaloMap[type lo]
+    let comm = getDefaultComm()
+    var hm = newSeq[HM](nd)
+    for mu in 0..<nd:
+      var offsets = newSeq[seq[int32]](0)
+      for nu in 0..<nd:
+        if nu == mu: continue
+        var t = newSeq[int32](nd)
+        t[nu] = 1
+        offsets.add t
+      hm[mu] = hl.makeHaloMap(comm, offsets)
+    toc "makeHaloMap"
+    type H = type makeHalo(hl, g[0])
+    var h = newSeq[H](nd)
+    for d in 0..<nd:
+      h[d] = makeHalo(hl, g[d])
+    toc "makeHalo"
+    for d in 0..<nd:
+      h[d].update hm[d], comm
+    toc "update"
+    type PT = evalType(g[0][0][0,0].re)
+    var pl = newSeq[float](6)
+    threads:
+      var spl: array[6,PT]
+      for i in g[0]:
+        var k = 0
+        for mu in 1..<4:
+          let n0 = hl.neighborFwd[mu][i]
+          for nu in 0..<mu:
+            let n1 = hl.neighborFwd[nu][i]
+            let a = g[mu][i] * h[nu][n0]
+            let b = g[nu][i] * h[mu][n1]
+            let p = redot(a, b)
+            spl[k] += p
+            inc k
+      for k in 0..<6:
+        let tpl = simdReduce spl[k]
+        threadCritical:
+          pl[k] += tpl[k]
+      #var tpl: array[6,float]
+      #for k in 0..<6: tpl[k] = simdReduce spl[k]
+      #threadSum tpl
+      #threadSingle:
+      #    for k in 0..<6: pl[k] = tpl[k]
+    rankSum pl
+    toc "pl"
+    let vf = 1.0/(g[0][0].nRows*lo.physVol)
+    echo pl * vf
+    echo 6.0 * g.plaq
+  testPlaq()
+
+  proc testStaple =
+    var st = lo.newGauge
+    var gf = lo.newGauge
+    for mu in 0..<nd: st[mu] := 0
+    tic("testStaple")
+    let hl = lo.makeHaloLayout([1,1,1,1],[1,1,1,1])
+    toc "makeHaloLayout"
+    type HM = HaloMap[type lo]
+    let comm = getDefaultComm()
+    var hm = newSeq[HM](nd)
+    for mu in 0..<nd:
+      var offsets = newSeq[seq[int32]](0)
+      for nu in 0..<nd:
+        var t = newSeq[int32](nd)
+        t[nu] = -1
+        offsets.add @t
+        if nu == mu: continue
+        t[nu] = 1
+        offsets.add @t
+        t[mu] = -1
+        offsets.add @t
+      hm[mu] = hl.makeHaloMap(comm, offsets)
+      #echo offsets
+    toc "makeHaloMap"
+    type H = type makeHalo(hl, g[0])
+    var h = newSeq[H](nd)
+    for mu in 0..<nd:
+      h[mu] = makeHalo(hl, g[mu])
+    toc "makeHalo"
+    for mu in 0..<nd:
+      h[mu].update hm[mu], comm
+    toc "update"
+    var cp = 1.0/3.0
+    threads:
+      for i in g[0]:
+        for mu in 1..<4:
+          let fmu = hl.neighborFwd[mu][i]
+          let bmu = hl.neighborBck[mu][i]
+          for nu in 0..<mu:
+            let fnu = hl.neighborFwd[nu][i]
+            let bnu = hl.neighborBck[nu][i]
+            let fmubnu = hl.neighborBck[nu][fmu]
+            let fnubmu = hl.neighborFwd[nu][bmu]
+            let t0 = h[mu][fnu] * h[nu][fmu].adj
+            let t1 = g[nu][i] * t0
+            st[mu][i] += cp*t1
+            let t2 = g[mu][i] * t0.adj
+            st[nu][i] += cp*t2
+            let t3 = h[mu][bnu] * h[nu][fmubnu]
+            let t4 = h[nu][bnu].adj * t3
+            st[mu][i] += cp*t4
+            let t5 = h[nu][bmu] * h[mu][fnubmu]
+            let t6 = h[mu][bmu].adj * t5
+            st[nu][i] += cp*t6
+    toc "st"
+    gf.gaugeForce(g)
+    #gf.gaugeForce2(g)
+    #gf.gaugeForce3(g)
+    toc "gaugeForce"
+    for mu in 0..<nd:
+      for e in st[mu]:
+        let s = g[mu][e] * st[mu][e].adj
+        st[mu][e].projectTAH s
+      echo gf[mu].norm2, " ", st[mu].norm2, " ", (gf[mu]-st[mu]).norm2
+  testStaple()
+
+  echoTimers()
+  qexFinalize()
+
+
+# u[mu=0] -> [0,1,0,0],[0,0,1,0],[0,0,0,1],
+#            [-1,0,0,0],[0,-1,0,0],[0,0,-1,0],[0,0,0,-1],
+#            [-1,1,0,0],[-1,0,1,0],[-1,0,0,1]
+
+# [1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1],
+# [-1,0,0,0],[0,-1,0,0],[0,0,-1,0],[0,0,0,-1],
+# [-1,1,0,0],[-1,0,1,0],[-1,0,0,1],
+# [1,-1,0,0],[0,-1,1,0],[0,-1,0,1],
+# [1,0,-1,0],[0,1,-1,0],[0,0,-1,1],
+# [1,0,0,-1],[0,1,0,-1],[0,0,1,-1]
+
diff --git a/src/experimental/pgag.nim b/src/experimental/pgag.nim
new file mode 100644
index 0000000..02c2b03
--- /dev/null
+++ b/src/experimental/pgag.nim
@@ -0,0 +1,785 @@
+import qex, gauge, physics/[qcdTypes,stagSolve]
+import times, macros, algorithm
+import hmc/metropolis
+import observables/sources
+import quda/qudaWrapper
+import hmc/agradOps
+
+qexinit()
+
+let
+  lat = intSeqParam("lat", @[8,8,8,8])
+  beta = floatParam("beta", 5.6)
+  tau0 = floatParam("tau", 0.04)
+  lam0 = floatParam("lam", 0.19)
+  sig0 = floatParam("sig", 0.3)
+  rho0 = floatParam("rho", 0.2)
+  theta0 = floatParam("theta", 0.2)
+  ups0 = floatParam("ups", 0.2)
+  upit = intParam("upit", 1)
+  gsteps0 = intParam("gsteps", 32)
+  nwarm = intParam("nwarm", 10)
+  trajs = intParam("trajs", 10)
+  seed0 = defaultComm.broadcast(int(1000*epochTime()))
+  seed = uint64 intParam("seed", seed0)
+  infn = stringParam("infn", "")
+  outfn = stringParam("outfn", "")
+var
+  tau = tau0
+  lam = lam0
+  sig = sig0
+  rho = rho0
+  ups = ups0
+  theta = theta0
+  gsteps = gsteps0
+
+macro echoparam(x: typed): untyped =
+  let n = x.repr
+  result = quote do:
+    echo `n`, ": ", `x`
+
+echoparam(beta)
+echoparam(tau)
+echoparam(lam)
+echoparam(sig)
+echoparam(rho)
+echoparam(theta)
+echoparam(ups)
+echoparam(upit)
+echoparam(gsteps)
+echoparam(nwarm)
+echoparam(trajs)
+echoparam(seed)
+echoparam(infn)
+echoparam(outfn)
+
+let
+  gc = GaugeActionCoeffs(plaq: beta, adjplaq: 0)
+  lo = lat.newLayout
+  vol = lo.physVol
+
+var r = lo.newRNGField(RngMilc6, seed)
+var R:RngMilc6  # global RNG
+R.seed(seed, 987654321)
+
+var g = lo.newgauge
+type
+  Gauge = typeof(g)
+  GaugeV = AgVar[Gauge]
+template newGaugeV(c: AgTape, x: Gauge): auto = newGaugeFV(c, x)
+case infn
+of "":
+  g.random r
+  #g.unit
+of "hot":
+  g.random r
+of "cold":
+  g.unit
+else:
+  echo "Loading gauge field from file: ", infn
+  let err = g.loadGauge infn
+
+echo "plaq: ", 6.0*g.plaq
+echo "gaugeAction2: ", g.gaugeAction2 gc
+echo "actionA: ", gc.actionA g
+
+var
+  p = lo.newgauge
+  f = lo.newgauge
+  g0 = lo.newgauge
+
+proc norm2*(x: seq): float =
+  for i in 0..<x.len:
+    result += x[i].norm2
+    #echo result
+  #result /= (x.len * x[0].l.physVol).float
+
+proc norm2v*(x: seq): float =
+  for i in 0..<x.len:
+    result += x[i].norm2
+    #echo result
+  result /= (x.len * x[0].l.physVol).float
+
+proc norm2subtract*(x: Field, y: float): float =
+  var s: evalType(norm2(toDouble(x[0])))
+  for i in x:
+    s += x[i].toDouble.norm2 - y
+  result = s.simdReduce
+  x.l.threadRankSum(result)
+
+proc norm2subtract*(x: seq, y: float): float =
+  for i in 0..<x.len:
+    result += x[i].norm2subtract(y)
+
+proc forceX(c: GaugeActionCoeffs, g: auto, f: auto) =
+  when defined(qudaDir):
+    #var ff = f.newOneOf
+    #gc.forceA(g, ff)
+    gc.qudaGaugeForce(g, f)
+    #echo "QEX:  ", ff[0].norm2
+    #echo "QUDA: ", f[0].norm2
+    #echo ff[0][0]
+    #echo f[0][0]
+    #echo (f[0]-ff[0]).norm2
+  else:
+    gc.forceA(g, f)
+
+var gutime = 0.0
+proc mdt(t: float) =
+  tic()
+  threads:
+    for mu in 0..<g.len:
+      for s in g[mu]:
+        g[mu][s] := exp(t*p[mu][s])*g[mu][s]
+  gutime += getElapsedTime()
+  toc("mdt")
+var gf2s = 0.0
+var ngf2 = 0
+var gftime = 0.0
+proc mdv(t: float) =
+  tic()
+  gc.forceX(g, f)
+  let gf2 = f.norm2v
+  gf2s += gf2
+  inc ngf2
+  threads:
+    for mu in 0..<f.len:
+      p[mu] -= t*f[mu]
+  gftime += getElapsedTime()
+  toc("mdv")
+
+proc mdvx(t: float) =
+  tic()
+  gc.forceX(g, f)
+  let gf2 = f.norm2v
+  gf2s += gf2
+  inc ngf2
+  threads:
+    for mu in 0..<f.len:
+      p[mu] -= t*f[mu]
+  gftime += getElapsedTime()
+  toc("mdv")
+
+var
+  gauges = newSeq[Gauge](0)
+  moms = newSeq[Gauge](0)
+  gaugevs = newSeq[GaugeV](0)
+  momvs = newSeq[GaugeV](0)
+  params = newSeq[FloatV](0)
+  ptemps = newSeq[FloatV](0)
+var tape = newAgTape()
+proc pushParam(x: float): FloatV =
+  params.add tape.newFloatV(x)
+  result = params[^1]
+
+gauges.add g0
+gaugevs.add tape.newGaugeV(g0)
+var p0 = lo.newGauge
+moms.add p0
+momvs.add tape.newGaugeV(p0)
+#params.add tape.newFloatV(tau)
+#let vtau = params[0]
+let vtau = pushParam(tau)
+
+proc pushTemp =
+  ptemps.add tape.newFloatV
+proc pushGauge =
+  gauges.add lo.newGauge
+  gaugevs.add tape.newGaugeV(gauges[^1])
+proc pushMom =
+  moms.add lo.newGauge
+  momvs.add tape.newGaugeV(moms[^1])
+
+proc addT(veps: FloatV) =
+  pushGauge()
+  exp(gaugevs[^1], veps, momvs[^1])
+  pushGauge()
+  mul(gaugevs[^1], gaugevs[^2], gaugevs[^3])
+
+proc addGf(g: GaugeV) =
+  pushMom()
+  gc.gderiv(momvs[^1], g)
+  pushMom()
+  mulna(momvs[^1], g, momvs[^2])
+  pushMom()
+  projtah(momvs[^1], momvs[^2])
+
+proc addGx(veps: FloatV, p,g: GaugeV) =
+  addGf(g)
+  pushMom()
+  xpay(momvs[^1], p, veps, momvs[^2])
+
+proc addG(veps: FloatV) =
+  addGx(veps, momvs[^1], gaugevs[^1])
+
+proc addGF(va, vb: FloatV) =
+  let p = momvs[^1]
+  addGf(gaugevs[^1])
+  pushMom()
+  exp(momvs[^1], vb, momvs[^2])
+  pushMom()
+  mul(momvs[^1], momvs[^2], gaugevs[^1])
+  addGx(va, p, momvs[^1])
+
+proc setupMDx =
+  pushTemp()
+  let vtau2 = ptemps[^1]
+  mul(vtau2, 0.5, vtau)
+  addT(vtau2)
+  addG(vtau)
+  addT(vtau2)
+
+proc setupMDa =
+  pushTemp()
+  let vtau2 = ptemps[^1]
+  mul(vtau2, 0.5, vtau)
+  pushTemp()
+  let vtau4 = ptemps[^1]
+  mul(vtau4, 0.25, vtau)
+  addT(vtau4)
+  addG(vtau2)
+  addT(vtau2)
+  addG(vtau2)
+  addT(vtau4)
+
+proc `-`(x: FloatV, y: FloatV): FloatV =
+  pushTemp()
+  result = ptemps[^1]
+  sub(result, x, y)
+proc `*`(x: float, y: FloatV): FloatV =
+  pushTemp()
+  result = ptemps[^1]
+  mul(result, x, y)
+proc `*`(x: FloatV, y: FloatV): FloatV =
+  pushTemp()
+  result = ptemps[^1]
+  mul(result, x, y)
+
+proc setupMD2 =
+  let lam = pushParam(lam0)
+  let lamt = lam * vtau
+  let sig = vtau - 2 * lamt
+  let tau2 = 0.5 * vtau
+  addT(lamt)
+  addG(tau2)
+  addT(sig)
+  addG(tau2)
+  addT(lamt)
+
+proc setupMD2p =
+  let lam = pushParam(lam0)
+  let lamt = lam * vtau
+  let sig = vtau - 2 * lamt
+  let tau2 = 0.5 * vtau
+  addG(lamt)
+  addT(tau2)
+  addG(sig)
+  addT(tau2)
+  addG(lamt)
+
+proc setupMD2g =
+  let lam = pushParam(lam0)
+  let lamt = lam * vtau
+  let sig = vtau - 2 * lamt
+  let tau2 = 0.5 * vtau
+  let rho = pushParam(rho0)
+  let rhot = rho * vtau * vtau * vtau
+  addG(lamt)
+  addT(tau2)
+  addGF(sig, rhot)
+  addT(tau2)
+  addG(lamt)
+
+proc setupMD3 =
+  let lam = pushParam(lam0)
+  let lamt = lam * vtau
+  let sig = pushParam(sig0)
+  let sigt = sig * vtau
+  let alp = 0.5*vtau - lamt
+  let bet = vtau - 2*sigt
+  addT(lamt)
+  addG(sigt)
+  addT(alp)
+  addG(bet)
+  addT(alp)
+  addG(sigt)
+  addT(lamt)
+
+proc setupMD3g =
+  let lam = pushParam(lam0)
+  let lamt = lam * vtau
+  let sig = pushParam(sig0)
+  let sigt = sig * vtau
+  let alp = 0.5*vtau - lamt
+  let bet = vtau - 2*sigt
+  let rho = pushParam(rho0)
+  let rhot = rho * vtau * vtau * vtau
+  addT(lamt)
+  addG(sigt)
+  addT(alp)
+  addGF(bet, rhot)
+  addT(alp)
+  addG(sigt)
+  addT(lamt)
+
+proc setupMD3p =
+  let lam = pushParam(lam0)
+  let lamt = lam * vtau
+  let sig = pushParam(sig0)
+  let sigt = sig * vtau
+  let alp = 0.5*vtau - lamt
+  let bet = vtau - 2*sigt
+  addG(lamt)
+  addT(sigt)
+  addG(alp)
+  addT(bet)
+  addG(alp)
+  addT(sigt)
+  addG(lamt)
+
+# 2x OMF2: 0.1 0.25 0.3
+proc setupMD4 =
+  let lam = pushParam(lam0)
+  let lamt = lam * vtau
+  let sig = pushParam(sig0)
+  let sigt = sig * vtau
+  let rho = pushParam(rho0)
+  let rhot = rho * vtau
+  let alp = 0.5*vtau - sigt
+  let bet = vtau - 2*lamt - 2*rhot
+  addT(lamt)
+  addG(sigt)
+  addT(rhot)
+  addG(alp)
+  addT(bet)
+  addG(alp)
+  addT(rhot)
+  addG(sigt)
+  addT(lamt)
+
+proc setupMD4g =
+  let lam = pushParam(lam0)
+  let lamt = lam * vtau
+  let sig = pushParam(sig0)
+  let sigt = sig * vtau
+  let rho = pushParam(rho0)
+  let rhot = rho * vtau
+  let alp = 0.5*vtau - sigt
+  let bet = vtau - 2*lamt - 2*rhot
+  let theta = pushParam(theta0)
+  let thetat = theta * vtau * vtau * vtau
+  addG(lamt)
+  addT(sigt)
+  addG(rhot)
+  addT(alp)
+  addGF(bet, thetat)
+  addT(alp)
+  addG(rhot)
+  addT(sigt)
+  addG(lamt)
+
+proc setupMD5 =
+  let lam = pushParam(lam0)
+  let lamt = lam * vtau
+  let sig = pushParam(sig0)
+  let sigt = sig * vtau
+  let rho = pushParam(rho0)
+  let rhot = rho * vtau
+  let theta = pushParam(theta0)
+  let thetat = theta * vtau
+  let alp = 0.5*vtau - lamt - rhot
+  let bet = vtau - 2*sigt - 2*thetat
+  addT(lamt)
+  addG(sigt)
+  addT(rhot)
+  addG(thetat)
+  addT(alp)
+  addG(bet)
+  addT(alp)
+  addG(thetat)
+  addT(rhot)
+  addG(sigt)
+  addT(lamt)
+
+# 6x LF: 0.08 0.17 0.17 0.17 0.17
+# 3x OMF2: 0.06 0.17 0.2 0.17 0.13
+proc setupMD6 =
+  let lam = pushParam(lam0)
+  let lamt = lam * vtau
+  let sig = pushParam(sig0)
+  let sigt = sig * vtau
+  let rho = pushParam(rho0)
+  let rhot = rho * vtau
+  let theta = pushParam(theta0)
+  let thetat = theta * vtau
+  let ups = pushParam(ups0)
+  let upst = ups * vtau
+  let alp = 0.5*vtau - sigt - thetat
+  let bet = vtau - 2*lamt - 2*rhot - 2*upst
+  addT(lamt)
+  addG(sigt)
+  addT(rhot)
+  addG(thetat)
+  addT(upst)
+  addG(alp)
+  addT(bet)
+  addG(alp)
+  addT(upst)
+  addG(thetat)
+  addT(rhot)
+  addG(sigt)
+  addT(lamt)
+
+proc update =
+  tape.run
+  for mu in 0..<g.len:
+    g[mu] := gauges[^1][mu]
+    p[mu] := moms[^1][mu]
+
+var
+  p2xv = tape.newFloatV
+  p2v = tape.newFloatV
+  gav = tape.newFloatV
+  hv = tape.newFloatV
+proc addAction(p: GaugeV, g: GaugeV) =
+  norm2subtract(p2xv, p, 8.0)
+  mul(p2v, 0.5, p2xv)
+  gaction(gc, gav, g)
+  add(hv, p2v, gav)
+
+proc setupAction =
+  tape.addTrack
+  addAction(momvs[0], gaugevs[0])
+  tape.addTrack
+  addAction(momvs[^1], gaugevs[^1])
+  tape.setTrack 0
+  #echo tape
+
+type
+  Met = ref object of MetropolisRoot
+    state: int
+
+proc init(m: var Met) =
+  m.new
+  var r = MetropolisRoot m
+  init(r)
+  m.verbosity = 1
+
+proc start*(m: var Met) =
+  tic()
+  m.state = 0
+  threads:
+    for i in 0..<g.len:
+      g0[i] := g[i]
+    p.randomTAH r
+    for i in 0..<p.len:
+      p0[i] := p[i]
+    threadBarrier()
+    #psi.gaussian r
+    #threadBarrier()
+    #g.rephase
+    #threadBarrier()
+    #stag.D(phi, psi, mass)
+    #threadBarrier()
+    #phi.odd := 0
+    #g.rephase
+  toc("init p, phi")
+
+proc getH*(m: Met): float =
+  tic()
+  var p2 = 0.0
+  var f2 = 0.0
+  #threads:
+  #  g.rephase
+  #stag.solve(psi, phi, mass, spa)
+  #toc("fa solve")
+  #echo "psi e: ", psi.even.norm2
+  #echo "psi o: ", psi.odd.norm2
+  threads:
+    #g.rephase
+    var p2t = 0.0
+    for i in 0..<p.len:
+      #p2t += p[i].norm2
+      p2t += p[i].norm2subtract(8.0)
+    threadMaster: p2 = p2t
+    #var psi2 = psi.norm2
+    #var psi2 = psi.norm2subtract(3.0)
+    #threadMaster: f2 = psi2
+  let
+    ga0 = gc.actionA g
+    #fa0 = 0.5*f2 - (1.5*vol).float
+    #t0 = 0.5*p2 - (16*vol).float
+    fa0 = 0.5*f2
+    t0 = 0.5*p2
+    h0 = ga0 + fa0 + t0
+  result = h0
+  toc("end getH")
+  if m.state == 0:
+    inc m.state
+    echo "Begin H: ",h0,"  Sg: ",ga0,"  Sf: ",fa0,"  T: ",t0
+    tape.setTrack 1
+    tape.run
+    echo &"      H: {hv.obj}  Sg: {gav.obj}  T: {p2v.obj}"
+    tape.setTrack 0
+  else:
+    echo "End H: ",h0,"  Sg: ",ga0,"  Sf: ",fa0,"  T: ",t0
+    tape.setTrack 2
+    tape.run
+    echo &"      H: {hv.obj}  Sg: {gav.obj}  T: {p2v.obj}"
+    tape.setTrack 0
+
+proc getCost(m: Met): seq[float] =
+  result.newSeq(0)
+  var cost = vtau.obj * m.avgPAccept
+  result.add cost
+  for i in 0..<params.len:
+    var costg = (if i==0: 1.0 else: 0.0)
+    if m.hNew > m.hOld:
+      costg = m.pAccept * (costg - vtau.obj * params[i].grad)
+    result.add costg
+
+proc checkGrad(m: Met) =
+  #let tg = vtau.grad
+  #let f0 = vtau.obj * exp(m.hOld-m.hNew)
+  let eps = 1e-7
+  var gs = newSeq[float](0)
+  for i in 0..<params.len: gs.add params[i].grad
+  for i in 0..<params.len:
+    let t = params[i].obj
+    params[i].obj += eps
+    update()
+    let h = m.getH
+    echo "params[",i,"] grad: ", gs[i]
+    echo "params[",i,"] diff: ", (h-m.hNew)/eps
+    params[i].obj = t
+  update()
+
+# w = sum_i w_i
+# sum = sum_i w_i x_i
+# sum2 = sum_i w_i x_i^2
+# X = sum/w
+# var = sum_i w_i (x_i - X)^2 = sum_i w_i x_i^2 - w X^2
+type DecayStat = object
+  n: int
+  fac: float
+  w: float
+  sum: float
+  sum2: float
+proc newDecayStat(f: float): DecayStat =
+  result.fac = f
+proc push(d: var DecayStat, x: float) =
+  inc d.n
+  let f = d.fac
+  let o = 1.0 - f
+  d.w = o*d.w + f
+  d.sum = o*d.sum + f*x
+  d.sum2 = o*d.sum2 + f*x*x
+proc mean(d: DecayStat): float =
+  result = d.sum / d.w
+proc variance(d: DecayStat): float =
+  let m = d.sum / d.w
+  let m2 = d.sum2 / d.w
+  result = sqrt(m2 - m*m)
+#var decay = 0.1
+var decay = 0.2
+var cgstat = newSeq[RunningStat](0)
+#var cgstat = newSeq[DecayStat](0)
+
+proc getgrad(m: Met) =
+  hv.grad = 1.0
+  tape.setTrack 2
+  tape.grad
+  tape.setTrack 0
+  #gc.gaugeDeriv2(gaugevs[^1].obj, gaugevs[^1].grad)
+  #for mu in 0..<g.len:
+  #  momvs[^1].grad[mu] := momvs[^1].obj[mu]
+  tape.grad
+  #let tg = vtau.grad
+  let cst = getCost(m)
+  echo "cost:  ", cst[0]
+  if cgstat.len < params.len: cgstat.setLen(params.len)
+  #if cgstat.len < params.len:
+  #  for i in 0..<params.len: cgstat.add newDecayStat(decay)
+  for i in 0..<cgstat.len:
+    cgstat[i].push cst[i+1]
+    let m = cgstat[i].mean
+    let v = cgstat[i].variance
+    #echo "costg: ", cgstat[i].mean, " ", cgstat[i].standardDeviationS
+    #echo "costg: ", cgstat[i].mean, " ", cst[i+1]
+    echo &"costg: {m/v:8.6f} {m:10.6f} {cst[i+1]}"
+
+proc updateParams(rate: float) =
+  let n = cgstat.len
+  var s = 0.0
+  var p = newSeq[float](n)
+  var g = newSeq[float](n)
+  for i in 0..<n:
+    let m = cgstat[i].mean
+    let v = cgstat[i].variance
+    let t = m / v
+    g[i] = t
+    s += t*t
+    p[i] = params[i].obj
+  echo "Params: ", p
+  s = rate*sqrt(n/s)
+  for i in 0..<n:
+    #let m = cgstat[i].mean
+    let d = s*g[i]
+    params[i].obj += d
+    p[i] = params[i].obj
+  echo "Params: ", p
+  for i in 0..<n:
+    clear cgstat[i]
+
+proc finish*(m: var Met) =
+  discard
+  #for mu in 0..<g.len:
+  #  g[mu] := g0[mu]
+  #  p[mu] := p0[mu]
+  #discard m.getH
+  #mdt(0.5*tau)
+  #echo "mdt: ", g.norm2, "  ", vg2.obj.norm2
+  #mdv(tau)
+  #echo "mdv: ", p.norm2, "  ", vp2.obj.norm2
+  #mdt(0.5*tau)
+  #echo "mdt: ", g.norm2, "  ", vg6.obj.norm2
+  # pacc = max(1, exp(hOld-hNew))
+  #if m.hNew <= m.hold:
+  #gc.gaugeDeriv2(gaugevs[^1].obj, gaugevs[^1].grad)
+  #for mu in 0..<g.len:
+  #  momvs[^1].grad[mu] := momvs[^1].obj[mu]
+  #tape.grad
+  ##let tg = vtau.grad
+  #let cst = getCost(m)
+  #echo "cost:  ", cst[0]
+  ##if cgstat.len < params.len: cgstat.setLen(params.len)
+  #if cgstat.len < params.len:
+  #  for i in 0..<params.len: cgstat.add newDecayStat(decay)
+  #for i in 0..<cgstat.len:
+  #  cgstat[i].push cst[i+1]
+  #  #echo "costg: ", cgstat[i].mean, " ", cgstat[i].standardDeviationS
+  #  echo "costg: ", cgstat[i].mean, " ", cst[i+1]
+  #checkGrad(m)
+
+proc disp(m: Met) =
+  let p = g.plaq
+  echo "plaq: ", 6.0*p
+  echo "tplaq: ", p.sum
+  echo "gauge force: ", sqrt(gf2s/ngf2)
+  gf2s = 0.0
+  ngf2 = 0
+  echo "rmsDeltaH: ", sqrt(m.avgDeltaH2)
+  echo "avgPAccept: ", m.avgPAccept
+
+proc accept*(m: var Met) =
+  disp(m)
+
+proc reject*(m: var Met) =
+  threads:
+    for i in 0..<g.len:
+      g[i] := g0[i]
+  disp(m)
+
+var alwaysAccept = false
+proc globalRand*(m: Met): float =
+  if not alwaysAccept:
+    result = R.uniform
+
+proc generate*(m: var Met) =
+  update()
+
+var m: Met
+m.init
+
+var src = lo.ColorVector
+var prop = lo.ColorVector
+let nt = lat[3]
+let nt2 = (nt+2) div 2
+#var picorr = newSeq[float](lat[3])
+var picorr = newSeq[RunningStat](nt)
+var picorrf = newSeq[RunningStat](nt2)
+var sp = newSolverParams()
+sp.r2req = 1e-20
+block:
+  var v: evalType(src{0})
+  v := 0
+  v[0] = 1
+  threads:
+    src.wallSource(0, v)
+  #echo src.norm2slice(3)
+
+proc measure =
+  #threads:
+  #  g.rephase
+  #stag.solve(prop, src, mass, sp)
+  #threads:
+  #  g.rephase
+  #picorr += prop.norm2slice(3)
+  #echo picorr
+  #let pic = prop.norm2slice(3)
+  #for i in 0..<nt:
+  #  picorr[i].push pic[i]
+  #  let k = min(i, nt-i)
+  #  picorrf[k].push pic[i]
+  #for i in 0..<nt:
+  #  echo i, " ", picorr[i].mean, " ", picorr[i].standardDeviationS
+  #for i in 0..<nt2:
+  #  let k = nt2 - 1 - i
+  #  echo i, " ", picorr[k].mean, " ", picorr[k].standardDeviationS
+  discard
+
+#setupMD2()
+#setupMD2p()
+setupMD2g()
+#setupMD3()
+#setupMD3g()
+#setupMD4()
+#setupMD4g()
+#setupMD6()
+
+setupAction()
+
+if nwarm > 0:
+  echo "Starting warmups"
+  #gsteps = 60
+  #fsteps = 40
+  #setupMDx()
+  alwaysAccept = true
+  for n in 1..nwarm:
+    m.update
+  m.clearStats
+
+echo "Starting HMC"
+#gsteps = gsteps0
+#fsteps = fsteps0
+#setupMD5()
+alwaysAccept = false
+gutime = 0.0
+gftime = 0.0
+#fftime = 0.0
+block:
+  tic()
+  for n in 1..trajs:
+    echo "Starting trajectory: ", n
+    tic()
+    m.update
+    getGrad(m)
+    if upit > 0:
+      if n mod upit == 0:
+        updateParams(0.001)
+    let tup = getElapsedTime()
+    measure()
+    let ttot = getElapsedTime()
+    echo "End trajectory update: ", tup, "  measure: ", ttot-tup, "  total: ", ttot
+  let et = getElapsedTime()
+  toc()
+  echo "HMC time: ", et
+  #let at = gutime + gftime + fftime
+  #echo &"gu: {gutime}  gf: {gftime}  ff: {fftime}  ot: {et-at}  tt: {et}"
+
+if outfn != "":
+  echo "Saving gauge field to file: ", outfn
+  let err = g.saveGauge outfn
+
+#echoTimers()
+qexfinalize()
diff --git a/src/gauge/hypsmear2.nim b/src/gauge/hypsmear2.nim
new file mode 100644
index 0000000..2dda26d
--- /dev/null
+++ b/src/gauge/hypsmear2.nim
@@ -0,0 +1,544 @@
+import base, layout, gauge, hypsmear, comms/halo, bitops
+getOptimPragmas()
+
+#nflop = 61632.0
+# l1[mu,nu] : g[mu](0,nu,-nu) g[nu](0,mu)(0,-nu)
+# l2[mu,nu] : g[mu](0) l1[mu,b](a,-a) l1[a,b](0,mu)(0,-a)
+# fl[mu] : g[mu](0) l2[mu,nu](nu,-nu) l2[nu,mu](0,mu)(0,-nu)
+# l2[mu,nu] : g[mu](0)+(0,b,-b)(a,-a)
+#             g[a](0,mu)(0,-a)(0,b,-b)
+#             g[b](0,mu)(0,-b)(0,a,-a)
+# g[0]: (0,+-1,+-1,+-1) (-1,1,+-1,+-1) (-1,+-1,1,+-1) (-1,+-1,+-1,1)
+
+type
+  HypTemps*[L,F,T] = object
+    gf: array[4,F]
+    hl: HaloLayout[L]
+    hm: HaloMap[L]
+    l1x: FieldArray[L.V,T]
+    l2x: FieldArray[L.V,T]
+    l1: FieldArray[L.V,T]
+    l2: FieldArray[L.V,T]
+    flx: FieldArray[L.V,T]
+    hgf: array[4,Halo[L,F,T]]
+    h1x: array[4,array[4,Halo[L,F,T]]]
+    h2x: array[4,array[4,Halo[L,F,T]]]
+    h1: array[4,array[4,Halo[L,F,T]]]
+    h2: array[4,array[4,Halo[L,F,T]]]
+    info: PerfInfo
+
+template proj(x,y: auto) = x.projectU y
+proc proj[T](x: T): T {.noInit,alwaysInline.} = result.proj x
+template projDeriv(r: auto, x: auto, c: auto) = r.projectUDeriv(x, c)
+template projDeriv(r: auto, u, x: auto, c: auto) = r.projectUDeriv(u, x, c)
+
+proc init*[L,F,T,G](ht: var HypTemps[L,F,T], gf: G, comm: Comm) =
+  tic "HypTemps init"
+  static: doAssert(type(gf[0]) is F)
+  static: doAssert(type(gf[0].l) is L)
+  static: doAssert(type(gf[0][0]) is T)
+  let lo = gf[0].l
+  doAssert(lo.nDim == 4)
+  ht.hl = lo.makeHaloLayout([1,1,1,1],[1,1,1,1])
+  toc "makeHaloLayout"
+  ht.l1x = newFieldArray2(lo,F,[4,4],mu!=nu)
+  ht.l2x = newFieldArray2(lo,F,[4,4],mu!=nu)
+  ht.l1 = newFieldArray2(lo,F,[4,4],mu!=nu)
+  ht.l2 = newFieldArray2(lo,F,[4,4],mu!=nu)
+  ht.flx = newFieldArray(lo,F,4)
+  for mu in 0..<4:
+    ht.gf[mu] = gf[mu]
+    ht.hgf[mu] = makeHalo(ht.hl, gf[mu])
+    for nu in 0..<4:
+      if nu == mu: continue
+      ht.h1x[mu][nu] = makeHalo(ht.hl, ht.l1x[mu,nu])
+      ht.h2x[mu][nu] = makeHalo(ht.hl, ht.l2x[mu,nu])
+      ht.h1[mu][nu] = makeHalo(ht.hl, ht.l1[mu,nu])
+      ht.h2[mu][nu] = makeHalo(ht.hl, ht.l2[mu,nu])
+  toc "makeHalo"
+  const hmtype = 0
+  when hmtype == 0:
+    var offsets = newSeq[array[4,int32]](0)
+    for k in 0..15:  # use all corners for simplicity
+      var t = [1'i32,1,1,1]
+      for i in 0..<4:
+        if k.testBit i: t[i] = -1
+      offsets.add t
+    ht.hm = ht.hl.makeHaloMap(comm, offsets)
+  elif hmtype == 1:
+    var offsets = newSeq[array[4,int32]](0)
+    for k in 0..15:
+      var t = [1'i32,1,1,1]
+      for i in 0..<4:
+        if k.testBit i: t[i] = -1
+      if k != 0 and k != 15: offsets.add t
+      for mu in 0..<4:
+        if t[mu] == 1:
+          t[mu] = 0
+          offsets.add t
+          t[mu] = 1
+    ht.hm = hl.makeHaloMap(comm, offsets)
+  toc "makeHaloMap"
+proc newHypTemps*[G](gf: G): auto =
+  type
+    L = type gf[0].l
+    F = type gf[0]
+    T = type gf[0][0]
+  var ht: HypTemps[L,F,T]
+  let comm = getDefaultComm()
+  ht.init(gf, comm)
+  ht
+
+proc symstaple(s: var auto, a: float, nu0,mu0,nu1,nu2,mu1,nu3: auto) =
+  let t0 = mu0 * nu1.adj
+  let t1 = nu0 * t0
+  s += a*t1
+  let t2 = mu1 * nu3
+  let t3 = nu2.adj * t2
+  s += a*t3
+
+proc symstaple(s: var auto, a: float, x,y: auto, i,fnu,fmu,bnu,fmubnu: SomeInteger,
+               p: static bool = false) =
+  when p:
+    template prj(m: auto): auto = proj(m)
+  else:
+    template prj(m: auto): auto = m
+  let xi = prj x[i]
+  let yfnu = prj y[fnu]
+  let xfmu = prj x[fmu]
+  let xbnu = prj x[bnu]
+  let ybnu = prj y[bnu]
+  let xfmubnu = prj x[fmubnu]
+  var t = xi * yfnu * xfmu.adj
+  t += xbnu.adj * ybnu * xfmubnu
+  s += a*t
+template symstaplep(s: var auto, a: float, x,y: auto, i,fnu,fmu,bnu,fmubnu: SomeInteger) =
+  symstaple(s, a, x,y, i,fnu,fmu,bnu,fmubnu, true)
+
+proc symderiv(s: var auto, x,y,cx,cy: auto, i,fnu,fmu,bnu,fmubnu: SomeInteger,
+              p: static bool = false): int =
+  mixin projectUflops
+  when p:
+    template prj(m: auto): auto = proj(m)
+  else:
+    template prj(m: auto): auto = m
+  if fnu>=0 and fmu>=0:
+    let xi = prj x[i]
+    let yfnu = prj y[fnu]
+    let xfmu = prj x[fmu]
+    s += cx[i] * yfnu * xfmu.adj
+    s += xi * cy[fnu] * xfmu.adj
+    s += xi * yfnu * cx[fmu].adj
+    result += 3*projectUflops(3) + 3*(18+2*3*66)
+  if bnu>=0 and fmubnu>=0:
+    let xbnu = prj x[bnu]
+    let ybnu = prj y[bnu]
+    let xfmubnu = prj x[fmubnu]
+    s += cx[bnu].adj * ybnu * xfmubnu
+    s += xbnu.adj * cy[bnu] * xfmubnu
+    s += xbnu.adj * ybnu * cx[fmubnu]
+    result += 3*projectUflops(3) + 3*(18+2*3*66)
+template symderivp(s: var auto, x,y,cx,cy: auto, i,fnu,fmu,bnu,fmubnu: SomeInteger): int =
+  symderiv(s, x,y,cx,cy, i,fnu,fmu,bnu,fmubnu, true)
+
+proc symderiv3(s: var auto, x,y,cx,cy: auto, i,fnu,fmu,bnu,fmubnu,nOut: SomeInteger,
+               p: static bool = false): int =
+  mixin projectUflops
+  when p:
+    template prj(m: auto): auto = proj(m)
+  else:
+    template prj(m: auto): auto = m
+  if fnu>=0 and fmu>=0:
+    let xi = prj x[i]
+    let yfnu = prj y[fnu]
+    let xfmu = prj x[fmu]
+    result += 3*projectUflops(3)
+    if i<nOut:
+      s += cx[i] * yfnu * xfmu.adj
+      result += 18+2*3*66
+    if fnu<nOut:
+      s += xi * cy[fnu] * xfmu.adj
+      result += 18+2*3*66
+    if fmu<nOUt:
+      s += xi * yfnu * cx[fmu].adj
+      result += 18+2*3*66
+  if bnu>=0 and fmubnu>=0:
+    let xbnu = prj x[bnu]
+    let ybnu = prj y[bnu]
+    result += 2*projectUflops(3)
+    if bnu<nOut:
+      let xfmubnu = prj x[fmubnu]
+      s += cx[bnu].adj * ybnu * xfmubnu
+      s += xbnu.adj * cy[bnu] * xfmubnu
+      result += projectUflops(3) + 2*(18+2*3*66)
+    if fmubnu<nOut:
+      s += xbnu.adj * ybnu * cx[fmubnu]
+      result += 18+2*3*66
+template symderiv3p(s: var auto, x,y,cx,cy: auto, i,fnu,fmu,bnu,fmubnu,nOut: SomeInteger): int =
+  symderiv3(s, x,y,cx,cy, i,fnu,fmu,bnu,fmubnu,nOut, true)
+
+proc symderiv2(s: var auto, x,y,cx,cy: auto, i,fnu,fmu,bnu,fmubnu: SomeInteger): int =
+  template prj(m: auto): auto = proj(m)
+  if fnu>=0 and fmu>=0:
+    let xi = prj x[i]
+    let yfnu = prj y[fnu]
+    let xfmu = prj x[fmu]
+    s += cx[i] * yfnu * xfmu.adj
+    s += xi * cy[fnu] * xfmu.adj
+    s += xi * yfnu * cx[fmu].adj
+    inc result, 3
+  if bnu>=0 and fmubnu>=0:
+    let xbnu = prj x[bnu]
+    let ybnu = prj y[bnu]
+    let xfmubnu = prj x[fmubnu]
+    s += cx[bnu].adj * ybnu * xfmubnu
+    s += xbnu.adj * cy[bnu] * xfmubnu
+    s += xbnu.adj * ybnu * cx[fmubnu]
+    inc result, 3
+
+#[
+# x: side links, y: middle links
+# snu: c xmu ynu'
+# snu: y xmu cnu'
+# snu: c' x-mu ynu-mu
+# snu: y' x-mu cnu-mu
+# smu: x cnu xmu'
+# smu: x-nu' c-nu xmu-nu
+# x, xmu, x-mu, x-nu, xmu-nu; y, ynu, ynu-mu; c, cnu, c-nu
+proc symderiv(f: var auto, a: float, nu0,mu0,nu1,nu2,mu1,nu3,c,cnu: auto) =
+  let t0 = mu0 * nu1.adj
+  let t1 = nu0 * t0
+  f += a*t1
+  let t2 = mu1 * nu3
+  let t3 = nu2.adj * t2
+  f += a*t3
+]#
+
+proc smear*[L,F,T,G](ht: HypTemps[L,F,T], coef: HypCoefs, fl: G) =
+  tic "HypTemps smear"
+  static: doAssert(type(fl[0]) is F)
+  static: doAssert(type(fl[0].l) is L)
+  static: doAssert(type(fl[0][0]) is T)
+  let lo = fl[0].l
+  doAssert(lo.nDim == 4)
+  let comm = getDefaultComm()
+  for mu in 0..<4:
+    ht.hgf[mu].halo := 0
+    ht.hgf[mu].update ht.hm, comm
+  toc "update"
+  const V = L.V
+  var
+    gf = ht.gf
+    hl = ht.hl
+    l1x = ht.l1x
+    l2x = ht.l2x
+    flx = ht.flx
+    hgf = ht.hgf
+    h1x = ht.h1x
+    h2x = ht.h2x
+    h1 = ht.h1
+    h2 = ht.h2
+  let
+    alp1 = coef.alpha1 / 2.0
+    alp2 = coef.alpha2 / 4.0
+    alp3 = coef.alpha3 / 6.0
+    ma1 = 1 - coef.alpha1
+    ma2 = 1 - coef.alpha2
+    ma3 = 1 - coef.alpha3
+  let nhalo = hl.nExt
+  threads:
+    let nc = gf[0][0].nrows
+    let staplesFlops = float((4*(6*nc+2*(nc-1))+4*2)*nc*nc*V)
+    let siteFlops1 = staplesFlops + float(2*nc*nc*V)
+    let siteFlops2 = 2*staplesFlops+float((2*nc*nc+12*projectUflops(nc))*V)
+    let siteFlops3 = 3*staplesFlops+float((2*nc*nc+19*projectUflops(nc))*V)
+    var flops = 0.0
+    # h1x[mu][nu] mu,nu,-nu,mu-nu
+    tfor i, 0..<nhalo:
+      for mu in 0..<4:
+        #for nu in 0..<4:
+        #  if nu == mu: continue
+        #  h1x[mu][nu][i] = 0
+        let fmu = hl.neighborFwd[mu][i]
+        if fmu<0: continue
+        for nu in 0..<4:
+          if nu == mu: continue
+          let fnu = hl.neighborFwd[nu][i]
+          if fnu<0: continue
+          let bnu = hl.neighborBck[nu][i]
+          if bnu<0: continue
+          let fmubnu = hl.neighborBck[nu][fmu]
+          if fmubnu<0: continue
+          h1x[mu][nu][i] = ma1 * hgf[mu][i]
+          symstaple(h1x[mu][nu][i], alp1, hgf[nu],hgf[mu], i,fnu,fmu,bnu,fmubnu)
+          h1[mu][nu][i].proj h1x[mu][nu][i]
+          flops += siteFlops1
+    #threadBarrier()
+    flops.threadSum
+    toc("1",flops=flops)
+    flops = 0
+    # h2x[mu][nu]: h1x[a][b]: 0,mu,-a,mu-a  h1x[mu][b]: a,-a
+    tfor i, 0..<nhalo:
+      for mu in 0..<4:
+        #for nu in 0..<4:
+        #  if nu == mu: continue
+        #  h2x[mu][nu][i] = 0
+        let fmu = hl.neighborFwd[mu][i]
+        if fmu<0: continue
+        for nu in 0..<4:
+          if nu == mu: continue
+          h2x[mu][nu][i] = ma2 * hgf[mu][i]
+          for a in 0..<4:
+            if a == mu or a == nu: continue
+            let b = 1+2+3-mu-nu-a
+            let fa = hl.neighborFwd[a][i]
+            if fa<0: continue
+            let ba = hl.neighborBck[a][i]
+            if ba<0: continue
+            let fmuba = hl.neighborBck[a][fmu]
+            if fmuba<0: continue
+            #symstaplep(h2x[mu][nu][i], alp2, h1x[a][b],h1x[mu][b], i,fa,fmu,ba,fmuba)
+            symstaple(h2x[mu][nu][i], alp2, h1[a][b],h1[mu][b], i,fa,fmu,ba,fmuba)
+          h2[mu][nu][i].proj h2x[mu][nu][i]
+          flops += siteFlops2
+    #threadBarrier()
+    flops.threadSum
+    toc("2",flops=flops)
+    flops = 0
+    # flx[mu]: h2x[nu][mu]: 0,mu,-nu,mu-nu  h2x[mu][nu]: nu,-nu
+    for i in lo:
+      for mu in 0..<4:
+        flx[mu][i] = ma3 * gf[mu][i]
+        let fmu = hl.neighborFwd[mu][i]
+        for nu in 0..<4:
+          if nu == mu: continue
+          let fnu = hl.neighborFwd[nu][i]
+          let bnu = hl.neighborBck[nu][i]
+          let fmubnu = hl.neighborBck[nu][fmu]
+          #symstaplep(flx[mu][i], alp3, h2x[nu][mu],h2x[mu][nu], i,fnu,fmu,bnu,fmubnu)
+          symstaple(flx[mu][i], alp3, h2[nu][mu],h2[mu][nu], i,fnu,fmu,bnu,fmubnu)
+        fl[mu][i].proj flx[mu][i]
+        flops += siteFlops3
+    #threadBarrier()
+    flops.threadSum
+    toc("3",flops=flops)
+  toc("threads end")
+
+proc smear2*[G](coef: HypCoefs, gf: G, fl: G, info: var PerfInfo) =
+  tic "smear2"
+  type
+    F = type gf[0]
+    L = type gf[0].l
+    T = type gf[0][0]
+  let comm = getDefaultComm()
+  var ht: HypTemps[L,F,T]
+  ht.init gf, comm
+  toc "init"
+  ht.smear coef, fl
+  toc "smear"
+
+proc force*[L,F,T,G,C](ht: HypTemps[L,F,T], coef: HypCoefs, f: G, chain: C) =
+  tic "HypTemps force"
+  static: doAssert(type(f[0]) is F)
+  static: doAssert(type(f[0].l) is L)
+  static: doAssert(type(f[0][0]) is T)
+  let lo = f[0].l
+  doAssert(lo.nDim == 4)
+  let comm = getDefaultComm()
+  const V = L.V
+  var
+    gf = ht.gf
+    hl = ht.hl
+    l1x = ht.l1x
+    l2x = ht.l2x
+    flx = ht.flx
+    hgf = ht.hgf
+    h1x = ht.h1x
+    h2x = ht.h2x
+    h1 = ht.h1
+    h2 = ht.h2
+    fc = newFieldArray(lo,F,4)
+    fl1 = newFieldArray2(lo,F,[4,4],mu!=nu)
+    fl2 = newFieldArray2(lo,F,[4,4],mu!=nu)
+    hf: array[4,Halo[L,F,T]]
+    hfc: array[4,Halo[L,F,T]]
+    hl1: array[4,array[4,Halo[L,F,T]]]
+    hl2: array[4,array[4,Halo[L,F,T]]]
+  let
+    alp1 = coef.alpha1 / 2.0
+    alp2 = coef.alpha2 / 4.0
+    alp3 = coef.alpha3 / 6.0
+    ma1 = 1 - coef.alpha1
+    ma2 = 1 - coef.alpha2
+    ma3 = 1 - coef.alpha3
+  for mu in 0..<4:
+    hf[mu] = makeHalo(ht.hl, f[mu])
+    hfc[mu] = makeHalo(ht.hl, fc[mu])
+    for nu in 0..<4:
+      if nu == mu: continue
+      hl1[mu][nu] = makeHalo(ht.hl, fl1[mu,nu])
+      hl2[mu][nu] = makeHalo(ht.hl, fl2[mu,nu])
+  let nhalo = hl.nExt
+  let nOut = hl.nOut
+  threads:
+    let nc = gf[0][0].nrows
+    var flops = 0.0
+    for mu in 0..<4:
+      for i in lo:
+        fc[mu][i].projDeriv(flx[mu][i], chain[mu][i])
+        f[mu][i] = ma3 * fc[mu][i]
+        fc[mu][i] *= alp3
+      tfor i, nOut..<nhalo:
+        hf[mu][i] = 0
+        hfc[mu][i] = 0
+    flops.threadSum
+    toc "fc"
+    flops = 0
+    tfor i, 0..<nhalo:
+      for mu in 0..<4:
+        let fmu = hl.neighborFwd[mu][i]
+        for nu in 0..<4:
+          if nu==mu: continue
+          hl2[mu][nu][i] = 0
+          let fnu = hl.neighborFwd[nu][i]
+          let bnu = hl.neighborBck[nu][i]
+          let fmubnu = if fmu<0: -1 else: hl.neighborBck[nu][fmu]
+          #let flp = symderiv3p(hl2[mu][nu][i], h2x[nu][mu],h2x[mu][nu],hfc[nu],hfc[mu], i,fnu,fmu,bnu,fmubnu,nOut)
+          let flp = symderiv3(hl2[mu][nu][i], h2[nu][mu],h2[mu][nu],hfc[nu],hfc[mu], i,fnu,fmu,bnu,fmubnu,nOut)
+          if flp > 0:
+            hl2[mu][nu][i].projDeriv(h2x[mu][nu][i], hl2[mu][nu][i])
+            hf[mu][i] += ma2 * hl2[mu][nu][i]
+            hl2[mu][nu][i] *= alp2
+            flops += V*(flp+54+2*projectUflops(nc)) # guess for projDeriv
+    flops.threadSum
+    toc "3", flops=flops
+    flops = 0
+    tfor i, 0..<nhalo:
+      for mu in 0..<4:
+        let fmu = hl.neighborFwd[mu][i]
+        var flp = 0
+        for nu in 0..<4:
+          if nu == mu: continue
+          hl1[mu][nu][i] = 0
+          #var t {.noInit.}: evalType(hl1[mu][nu][i])
+          #t := 0
+          for a in 0..<4:
+            if a == mu or a == nu: continue
+            let b = 1+2+3-mu-nu-a
+            let fa = hl.neighborFwd[a][i]
+            let ba = hl.neighborBck[a][i]
+            let fmuba = if fmu<0: -1 else: hl.neighborBck[a][fmu]
+            #flp += symderivp(hl1[mu][nu][i], h1x[a][nu],h1x[mu][nu],hl2[a][b],hl2[mu][b], i,fa,fmu,ba,fmuba)
+            flp += symderiv(hl1[mu][nu][i], h1[a][nu],h1[mu][nu],hl2[a][b],hl2[mu][b], i,fa,fmu,ba,fmuba)
+            #flp += symderiv(t, h1[a][nu],h1[mu][nu],hl2[a][b],hl2[mu][b], i,fa,fmu,ba,fmuba)
+          if flp > 0:
+            hl1[mu][nu][i].projDeriv(h1x[mu][nu][i], hl1[mu][nu][i])
+            hf[mu][i] += ma1 * hl1[mu][nu][i]
+            hl1[mu][nu][i] *= alp1
+            #t.projDeriv(h1x[mu][nu][i], t)
+            #hf[mu][i] += ma1 * t
+            #hl1[mu][nu][i] = alp1 * t
+            flops += V*(flp+54+2*projectUflops(nc)) # guess for projDeriv
+    flops.threadSum
+    toc "2", flops=flops
+    flops = 0
+    tfor i, 0..<nhalo:
+      for mu in 0..<4:
+        let fmu = hl.neighborFwd[mu][i]
+        for nu in 0..<4:
+          if nu == mu: continue
+          let fnu = hl.neighborFwd[nu][i]
+          let bnu = hl.neighborBck[nu][i]
+          let fmubnu = if fmu<0: -1 else: hl.neighborBck[nu][fmu]
+          let flp = symderiv(hf[mu][i], hgf[nu],hgf[mu],hl1[nu][mu],hl1[mu][nu], i,fnu,fmu,bnu,fmubnu)
+          flops += V*flp
+    flops.threadSum
+    toc "1", flops=flops
+  for mu in 0..<4:
+    hf[mu].updateRev ht.hm, comm
+  toc("end")
+
+when isMainModule:
+  import qex
+  qexInit()
+  tic("main")
+  var defaultLat = @[4,4,4,4]
+  defaultSetup()
+  var seed = 987654321'u
+  #var rng = newRngField(lo, RngMilc6, seed)
+  var rng = newRngField(lo, MRG32k3a, seed)
+  g.gaussian rng
+  toc "gaussian"
+  echo 6.0 * g.plaq
+  toc "plaq"
+
+  var info: PerfInfo
+  var coef: HypCoefs
+  coef.alpha1 = 0.4
+  coef.alpha2 = 0.5
+  coef.alpha3 = 0.5
+  echo coef
+
+  var fl = lo.newGauge()
+  var fl2 = lo.newGauge()
+
+  proc testSmear =
+    resetTimers()
+    tic "start"
+    coef.smear(g, fl, info)
+    toc "smear"
+    coef.smear2(g, fl2, info)
+    toc "smear2"
+    echoTimers()
+    echo fl.plaq
+    echo fl2.plaq
+  #testSmear()
+
+  proc testForce =
+    var f = lo.newGauge()
+    var f2 = lo.newGauge()
+    var c = lo.newGauge()
+    c.gaussian rng
+    block:
+      let fn = coef.smearGetForce(g, fl, info)
+      fn(f, c)
+    resetTimers()
+    block:
+      tic()
+      let fn = coef.smearGetForce(g, fl, info)
+      toc "smear"
+      fn(f, c)
+      toc "force"
+    block:
+      tic()
+      let ht = newHypTemps(g)
+      toc "newHypTemps"
+      freezeTimers()
+      ht.smear(coef, fl2)
+      ht.force(coef, f2, c)
+      thawTimers()
+      block:
+        tic()
+        ht.smear(coef, fl2)
+        toc "HTsmear"
+        ht.force(coef, f2, c)
+        toc "HTforce"
+    echoTimers()
+    echo fl.plaq
+    echo fl2.plaq
+    echo f.plaq
+    echo f2.plaq
+    var d = fl[0].newOneOf
+    var d2fl = newSeq[float](fl.len)
+    var d2f = newSeq[float](f.len)
+    for mu in 0..<fl.len:
+      d := fl2[mu] - fl[mu]
+      d2fl[mu] = sqrt(d.norm2/fl[mu].norm2)
+      d := f2[mu] - f[mu]
+      d2f[mu] = sqrt(d.norm2/f[mu].norm2)
+    echo "error smear: ", d2fl
+    echo "error force: ", d2f
+
+  testForce()
+
+  qexFinalize()
diff --git a/src/gauge/staples.nim b/src/gauge/staples.nim
index fb53619..4940e70 100644
--- a/src/gauge/staples.nim
+++ b/src/gauge/staples.nim
@@ -22,6 +22,54 @@ proc startCornerShifts*[T](u0: openArray[T]): auto =
           s[mu][nu].startSB(u[mu][ix])
   return s
 
+proc startPlaqShifts*[T](u0: openArray[T]): auto =
+  var u = @u0
+  var sf:seq[seq[ShiftB[type(u[0][0])]]]
+  var sb:seq[seq[ShiftB[type(u[0][0])]]]
+  let nd = u.len
+  sf.newSeq(nd)
+  sb.newSeq(nd)
+  for mu in 0..<nd:
+    sf[mu].newSeq(nd)
+    sb[mu].newSeq(nd)
+    for nu in 0..<nd:
+      sb[mu][nu].initShiftB(u[mu], nu, -1, "all")
+      if mu!=nu:
+        sf[mu][nu].initShiftB(u[mu], nu, 1, "all")
+  threads:
+    for mu in 0..<nd:
+      sb[mu][mu].startSB(u[mu][ix])
+    for mu in 0..<nd:
+      for nu in 0..<nd:
+        if mu!=nu:
+          sf[mu][nu].startSB(u[mu][ix])
+          sb[mu][nu].startSB(u[mu][ix])
+  return (sf,sb)
+
+proc finishPlaqShifts*[T](u0,b0: openArray[T], sf,sb: seq): auto =
+  var u = @u0
+  var b = @b0
+  var sc:seq[seq[ShiftB[type(u[0][0])]]]
+  let nd = u.len
+  sc.newSeq(nd)
+  for mu in 0..<nd:
+    sc[mu].newSeq(nd)
+    for nu in 0..<nd:
+      if mu!=nu:
+        sc[mu][nu].initShiftB(u[mu], nu, 1, "all")
+  threads:
+    for mu in 0..<nd:
+      for ir in b[mu]:
+        localSB(sb[mu][mu], ir, assign(b[mu][ir], it), u[mu][ix])
+    for mu in 0..<nd:
+      boundarySB(sb[mu][mu], assign(b[mu][ir], it))
+    threadBarrier()
+    for mu in 0..<nd:
+      for nu in 0..<nd:
+        if mu!=nu:
+          sc[mu][nu].startSB(b[mu][ix])
+  return (sc)
+
 proc startStapleShifts*[T](u: openArray[T]): auto =
   var s:seq[seq[seq[ShiftB[type(u[0][0][0])]]]]
   let nd = u.len