mdlj — wafer-scale short-range molecular dynamics

Classical MD: N point particles in a periodic cubic box, a short-ranged spherically cut-off pair potential, velocity-Verlet integration. The box is cut into nx × ny × nz cells, one cell per PE, and mapped onto a 2D patch of the wafer by folding the third dimension into the two physical directions. Each timestep every PE exchanges a halo of neighbour-cell atoms, builds a pair list over a half shell, computes forces, reduces partial forces back to the owning PEs, integrates, and migrates atoms that crossed a cell boundary.

Despite the name, the pair potential is the purely repulsive U(r) = 0.1·r^-12 (ljforce() in mdlj.w), as description.txt says. The initial condition comes from a LAMMPS trajectory (lj2.5-traj.dump, 108 atoms); the kernel cuts the potential off at rcut = 2.5.

---

Algorithm

Pair potential and force

ljforce() works on the flat interaction list, with rinv[m] = 1/r_ij:

u = 0.1 · r^-12                      (eval; one term per pair)
F = 0.1 · 12 · r^-14 · (r_i − r_j)   (df[m] = 1.2·r^-14; dxij[m][k] *= df[m])

Hard cutoff at rcut (d2 < rcut2, tested in detect_platelet()); the potential is not shifted at the cutoff. eval is scaled by 0.1 and not by 1/2, which is right because the pair list holds each pair once.

1/r comes from rsqrt_rinv(): a bit-trick seed (magic = 0x5f37 minus the right-shifted high word of r²) plus three Newton iterations y ← y(1.5 − 0.5·x·y²), each three ∀j loops over the pair list.

Time integration

Velocity Verlet, split over two loops so the step count is not limited by fp32 counter resolution:

for iter in [0, niter):               # outer, fp32
  for inner in [0, ninner]:           # inner, fp32
    halo exchange, pair list, forces, force reduction
    dtv = (inner == 0 or inner == ninner) ? 0.5·dt : dt
    v += F · invmass[type] · dtv
    if inner < ninner:  x += v · dt
    atom migration

The half kick at inner == ninner and the one at inner == 0 of the next outer iteration happen at the same positions, so together they form the full dt kick — exactly velocity Verlet over ninner timesteps per outer iteration. make-inputbin.py ships niter = 10, ninner = 1, nsample = 1: 10 timesteps, 2 force evaluations each.

Two known gaps, from description.txt: net momentum is never re-zeroed, and the neighbour list is rebuilt every timestep (rebuild/rebuildnext are compile-time 1).

Decomposition and the NZX/NZY fold

nz = nzx·nzy, so the program occupies nxtot × nytot = (nx·nzx) × (ny·nzy) tiles. Two mappings compose:

1. Serpentine folding of each logical axis. jtox/xtoj (in layoutindex.py and 3dmd_route3sdr.cc) map PE column group pix0 to logical cell pix0/2 when pix0 is even and nx-1-(pix0-1)/2 when odd. Logical cells 0 and nx-1 therefore land on physically adjacent groups, so the periodic logical ring becomes a local physical path with no long wrap wire. hix = pix0 % 2 says which of the two interleaved lanes a tile is on; hiy likewise for rows.
2. Folding z into x and y. Within each group of nzx consecutive PE columns, izx selects one of nzx z-slices (reversed when hix = 1); rows carry izy. (izx, izy) becomes a single z cell index via a Hamiltonian cycle over the nzx × nzy block, emitted as zidx<nzx>x<nzy>.dat. For nzx=4, nzy=2:

izy=0:   0   7   6   5
izy=1:   1   2   3   4         (closed cycle 0→1→2→…→7→0)

chain_jtox(nx,nzx,pix) = nzx·jtox(nx,pix0) + izx (in 3dmd_route3sdr.cc) is a tile's position along the folded chain running through a physical row; the Python chain_jtox in layoutindex.py returns the tuple (jtox(nx,pix0), izx, hix) instead. One logical x-cell step is +nzx along the chain, i.e. 2·nzx physical hops (the tiles in between belong to the other lane and only forward). With nx=8, nzx=4:

pix   : 0  1  2  3 | 4  5  6  7 | 8  9 10 11 | ...
pix0  :     0      |     1      |     2      |
hix   : 0  0  0  0 | 1  1  1  1 | 0  0  0  0 |
x     : 0  0  0  0 | 7  7  7  7 | 1  1  1  1 |
izx   : 0  1  2  3 | 3  2  1  0 | 0  1  2  3 |
chain : 0  1  2  3 |31 30 29 28 | 4  5  6  7 |

so chain index 0 (column 0) talks to chain index 4 (column 8) — same izx, logical x one higher.

Halo exchange, the "platelet", force reduction

range = 3 cells each way. Communication is a shift-register relay: at step r a PE forwards the block it received at step r-1 and receives the block r+1 cells away.

• getcount() relays atom counts — three rounds on each of x+, x−, y+, y−, filling {xm,xp,ym,yp}count[1..3] (slot 0 is the PE's own nlocal). Then zpcount[0] is set to the whole x-extent (Σ xmcount[0..3] + Σ xpcount[1..3] = local + both x halos) and three more rounds relay that along z−.
• exchange() relays the 4-word position+id records over the same five directions, using those counts and the derived *first offsets.

Afterwards pos holds:

[0,     nlocal)   owned atoms
[ymptr, y0ptr)    atoms from cells y-1 .. y-3
[y0ptr, ypptr)    second copy of the owned atoms
[ypptr, yend)     atoms from cells y+1 .. y+3
[xmptr, x0ptr)    atoms from cells x-1 .. x-3
[x0ptr, xpptr)    third copy of the owned atoms
[xpptr, xend)     atoms from cells x+1 .. x+3
[zpptr, zend)     the whole x-line (x-halo + local + x+halo) of cells z+1 .. z+3

detect_platelet() builds the pair list over a half shell (the "platelet"), so each pair is found on exactly one PE:

• i ∈ [ymptr, yend) (whole y-line) × j ∈ [xpptr, zend) (x+ halo, z+ halos);
• i ∈ [y0ptr, ypptr) (local copy) × j ∈ (i, yend) (local–local pairs and the y+ halo).

Minimum-image wrapping is applied inline per component with the copysign(convmagic, ·) round-to-nearest trick (convmagic = 2^23). remap_positions() recovers each atom's type from the sign of the id in pos[i][3] (type 1 stored negative, type 2 positive — see make-inputbin.py); klist = itype + jtype becomes a 0/1/2 pair-type index and atomtype of the owned atoms becomes the 0/1 index into invmass[].

Because partial forces were accumulated onto halo copies, reduce() runs the relay backwards with accumulating receives (<+-) on z+, x+, x−, y+, y−. The two extra local copies are then folded in: ff1[j] += ff1[3·x0ptr+j] + ff1[3·y0ptr+j].

Atom migration

After integration, for k = 0,1,2 each PE minimum-images pos[j][k] about center[k] and tests it against cell[k][0..1]. Departing atoms are packed as 8-word records (4 position words + 3 velocity words, padded to two pos rows), counts are exchanged with the two nearest neighbours in that direction, then the records are exchanged and appended. One hop, same colours as the halo relay.

Fabric topology, colours, routing

Ten colours carry simulation traffic; mdlj.paint binds the sockets:

colour	socket	meaning
0 / 1	s_xm	x− low lane (hix=0) / high lane (hix=1)
2 / 3	s_xp	x+ low / high lane
4 / 5	s_ym	y− low lane (hiy=0) / high lane
6 / 7	s_yp	y+ low / high lane
8 / 9	s_zm / s_zp	z− / z+ along the Hamiltonian cycle
12–15	ping	global barrier network (the socket is colour 12)

Paint binds only the low-lane colours; at startup main rewrites the input- and output-queue colour fields (set_iq_color/set_oq_color over the registers at 0x7bc0/0x7b80) to +1 where hix (0,2 → 1,3) or hiy (4,6 → 5,7), so each tile participates in its own lane's colour and silently forwards the other lane's.

The router/switch configuration is not written by Paint. 3dmd_route3sdr generates it and the Makefile misets it under the magic symbols router (order y x r) and switch (order y x r n), covering 24 colours over the (nxtot+2) × (nytot+2) array including the moat ring. For colours 0–7, layout2Dpbc() gives every tile a four-position switch state machine:

pos 0:  upstream -> downstream   (forward)
pos 1:  core     -> downstream   (inject)    "tail"
pos 2:  core     -> core         (loopback)
pos 3:  upstream -> core         (deliver)   "head"

The initial position follows the tile's chain index modulo nzx+1 (nzy+1 for y): tail at ≡ 0, head at ≡ nzx, forwarders between, plus a fixup that makes a trailing partial segment all heads. Every data burst is followed by the two control wavelets xycmds ((adv,adv,rst) then (adv)) or, for z, zcmds ((adv,rst) then (adv)), which march the state machine so that one send + control round equals one logical-cell shift. mdlj.w documents the sequence beside the xycmds/zcmds declarations (* marks the no-CE bit), and the comment at the end of chain_xtoj() in 3dmd_route3sdr.cc gives the switch positions it steps through.

The periodic wrap of the serpentine is closed by the moat, the one-tile ring around the core rectangle (built in mdlj.paint, name('moat), no program on it). Moat tiles are RtoRswap/LtoLswap/UtoUswap/DtoDswap: receive from the inward side, transmit straight back, and swap the colour to the partner of the pair. A wavelet running off the end of the low lane turns around and returns as the high lane — the fold turnaround.

Colours 0–7 are emitted with the teardown bit (1<<15) set; main clears it with clear_teardown() on all eight before any communication.

Ping network (global barrier)

pingorigin.w runs on a single tile, rect(0 0 1 1) (the corner moat tile at Paint (0, 0)), named ping. Colours 12–15 form a comb across the padded array, hard-coded in 3dmd_route3sdr.cc as pingmat[3][3] — one pattern per corner/edge/interior tile class. Every nsample outer iterations each core sends ping[] <- control(adv[0]); and then receives twice (junk <- ping[];), while the origin sends a probe plus control(adv[0]), waits for it to come back (which only happens once every core is at the barrier), sets frame[0], then sends a second value plus control(rst[0]) to release everyone. The optional external_driver handshake (spin on cont[0] until a host writes it) is compiled out (external_driver = 0).

Phases of one inner iteration

 1. getcount()          relay counts: 3 rounds × (x±, y±), then 3 × z−
 2. pointer setup       ymptr/y0ptr/ypptr/yend/xmptr/x0ptr/xpptr/xend/zpptr/zend
 3. copy local atoms    into the x0 and y0 slots
 4. exchange()          relay 4-word position+id records over x±, y±, z−
 5. remap_positions()   id sign -> atomtype
 6. detect_platelet()   half-shell pair list: dxij, dist2, ilist, jlist, klist
 7. rsqrt_rinv()        rinv = 1/sqrt(dist2), vectorised, 3 Newton steps
 8. zeroforce()         ff1[0 .. 3·zend) = 0
 9. ljforce()           energy, df = 1.2·r^-14, scatter ±F onto i and j
10. reduce()            accumulate partial forces back on z+, x±, y±
11. fold local copies   ff1[j] += ff1[3·x0ptr+j] + ff1[3·y0ptr+j]
12. integrate           half/full kick, drift
13. migrate             per-axis ownership test, count exchange, record exchange

Every nsample outer iterations (when isample == nsample) the PE computes keval = 0.5 Σ v², calls recordtime() (the three 16-bit words of the timestamp counter at 0x7b30, plus eval/keval, into the 128-deep ring buffers tstamp, evec, kevec) and joins the ping barrier. recordtime() also runs once before the loop, so a run records 1 + ⌊niter / nsample⌋ frames.

---

Source files

File	Role
bench.paint	what paint is invoked on: a copy of mdlj.paint
mdlj.paint	floorplan: core rectangle, moat ring, ping origin, socket/colour bindings
pinggrid.paint	Paint functions that paint a ping-colour grid; mdlj.paint includes it but does not call them
mdlj.w	the per-PE MD program
params.w	paramstruct (int16 per-PE params) and fparamstruct (fp32 per-PE params)
pingorigin.w	single-tile barrier driver on the corner moat tile
3dmd_route3sdr.cc	C++ generator for the router/switch tables and the z index map (SDR only)
3dmd_route3sdr	the generator, compiled; the Makefile runs it and has no rule to rebuild it from the .cc. The experimental SDK ships it compiled from the .cc
Makefile	PAINTVARS, PAINTFLAGS, the miset/miget wiring, the check steps
test.rc	one row per configuration
make-inputbin.py	reads the LAMMPS dump, bins atoms into cells, writes params/fparams/posin/velin .bin
check-exchange.py	correctness oracle for the halo bookkeeping
ecalc.py	energy/timing post-processing and the pass/fail energy check
layoutindex.py	Python jtox/xtoj/chain_jtox/chain_xtoj; its chain_jtox returns a tuple
readdump.py, writedump.py	LAMMPS dump reader / writer helpers
posbin2dump.py	writes per-PE atom data as a LAMMPS dump (manual tool; see Correctness checking)
frame.sh	a host loop that pauses a running system, saves a frame and resumes it; it names a map file and symbols (ping.iter, core00.h) this app does not have, and the Makefile does not run it
lj2.5-traj.dump	LAMMPS trajectory, 108 atoms, box ≈ 10.2599³; readdump takes the first frame (timestep 25000)
data1.dat	reference energies, 10 rows × (potential, kinetic, total)
copyright.txt	the header at the top of most source files

mdlj.w. After the import of params and a few globals, the file lays out a hand-written memory map: fpmem[1500] and its integer alias xpmem[3000] form one arena, and every array (countblock/firstblock with their five per-direction views, vel, pos, ff, dxij, rinv, tmp1, tmp2) is placed with an explicit address(&previous + size) so that varend ≤ memend — a comment says to verify that in bench.map. The extern aliases (pos0, vel1, ff1, countblock, firstblock, rinv, dxij) exist so the Makefile can miset/miget them by name. Then come register helpers (set_iq_color, get_iq_color, set_oq_color, get_oq_color, clear_teardown), an addressinit() whose body is commented out, and initparams(). The communication layer is getcount(), exchange(), reduce(). Then come remap_positions() (id → type), compute_dist2() (never called), detect_platelet() (the scalar half-shell pair-list build), rsqrt_rinv(), zeroforce(), ljforce(), and timerinit()/recordtime(). main does the colour fixup and teardown clearing, then runs the nested iteration loop. It also fills ptrs[22] with diagnostics: nlocal, the five halo counts, the eight region pointers, nint, the candidate-pair counters icount1/icount2, the inner-iteration count, and running maxima of nlocal, zend and nint (slots 18/19/20 — the right way to size NLOCMAX, NMAX, NINTMAX).

params.w. Layout of the per-PE blobs that miset writes: paramstruct = {nx ny nzx nzy, ix iy iz izx izy hix hiy, nlocal nsample} and fparamstruct = {dt rcut rscale rhoscale, niter ninner, invmass[2], box[3], boxinv[3], cell[3][2], center[3]}. rscale/rhoscale are 0 and unused. cell[k][0..1] are the cell bounds relative to the cell centre and center[k] the cell centre in box coordinates; both feed the migration test and the minimum-image code.

mdlj.paint. Declares colours 0–9 plus pingcol = color(12), sizes the core as rect(1 1 nxtot nytot), places code('mdlj [...]) on it with the socket bindings and the parameters nlocmax/nintmax/nmax (named core), builds the four-sided moat (no code, routing only), places pingorigin on rect(0 0 1 1) (named ping), and unions the three into all. It includes pinggrid.paint, but the call to its ping-grid function is commented out. It contains no paint(tile() c [...]) routing for the simulation colours — all of that arrives via miset.

pingorigin.w. It reads the core's niter/nsample from its own params/fparams copies (which the Makefile misets from params.bin.origin/fparams.bin.origin) and drives the probe/release handshake.

3dmd_route3sdr.cc. Run as ./3dmd_route3sdr NX NY NZX NZY, it writes zidx<nzx>x<nzy>.dat (Hamiltonian cycle index per (izx,izy)), rconfig<nxtot>x<nytot>.bin (uint16[nytot+2][nxtot+2][24] router words) and swconfig<nxtot>x<nytot>.bin (uint16[...][24][4] switch words). Its pieces: jtox/xtoj/chain_jtox/chain_xtoj (the fold); routerreg/router (bit encoders for the router word, two switch words and the position/status word); makeZloop (Hamiltonian cycle on the nzx × nzy block, recursing with a transpose when na is odd); makeZchain (walks the cycle to assign indices, then converts it to a red/black one-hop marching chain — "first 2k sends to 2k+1, then 2k+1 sends to 2k+2"); layout2Dpbc (the colour 0–7 chains and forwarders); lookup/zlookup (replicate the z pattern into every block, mirrored by xflip/yflip per hix/hiy); and main, which also emits the moat swap routers and the pingmat table. It prints arrow diagrams of the x-row, y-column and z-loop routing to stdout, which land in test-<tag>/make-log — the quickest sanity check on a new geometry.

---

Parameters

Set per line in test.rc; the first seven reach Paint as bare globals via PAINTVARS.

Name	Controls
NLOCMAX	maximum atoms owned by one PE. Sizes vel[NLOCMAX][3] and the pos0/vel1 miset/miget windows. Must exceed the peak of ptrs[18].
NMAX	maximum atoms resident on one PE (owned + two local copies + all halos). Sizes pos[NMAX][4] and ff[NMAX][3]. Must exceed the peak of ptrs[19] (zend).
NINTMAX	maximum pairs in the interaction list. Sizes dxij, rinv, tmp1, tmp2, dist2, df, ilist, jlist, klist. detect_platelet() clamps rather than overflows, so an undersized value silently drops pairs — watch ptrs[20].
NX, NY	logical cells along x and y, i.e. PE column/row groups.
NZX, NZY	the factorisation nz = NZX·NZY of the logical z extent into the two physical directions. Also sets the marching-chain period (NZX+1, NZY+1).
ARCH	sdr (WSE-3), passed to Paint as arch=. Only sdr works here: the route generator is SDR-specific and the config-register addresses in mdlj.w (0x7c00, 0x7bc0, 0x7b80, 0x7b30) are marked SDR.

Derived: NXTOT = NX·NZX, NYTOT = NY·NZY, PE count NXTOT·NYTOT plus the moat ring. The shipped configuration NX=8 NY=8 NZX=4 NZY=2 is 8×8×8 = 512 cells on a 32×16 patch — hence the tag t512.

Physics parameters are not in test.rc; they are hard-coded near the top of make-inputbin.py: dt = 0.001, rcut = 2.5, invmass = [1, 1], niter = 10, ninner = 1, nsample = 1, dump file ../lj2.5-traj.dump. Per description.txt, to change them cheaply delete test-<tag>/bench.img but keep the directory and rebuild: the compiled program is kept, and only the initial conditions are regenerated before the simulation reruns.

---

Building and running

The experimental SDK installs this example beside the tutorial, as projects/mdlj, with the route generator already compiled. Run every row of its test.rc with tungsten-tutorials -f projects/mdlj, or work by hand in a copy of the directory, set up as for the tutorials, where <sdk> is the SDK's directory:

cp -r <sdk>/share/tutorial/projects/mdlj . && cd mdlj
make TEST=t512                 # build, simulate and check
make TEST='*'                  # every row of test.rc
make TEST=t512 GOAL=bench.img  # build the image only, no simulation
make TEST=t512 GOAL=debugx     # also trace, and unpack the debug viewer's files
make clean                     # remove all test-*/ directories

test.rc has 1 row, tagged both t512 and pr_smoke. $(SIMFABRIC) runs the image on the Cerebras SDK simulator. MODE (default sim) only names the outputs (out-$(MODE).bin); sim is the only mode.

Output lands in test-<tag>/: make-log (full build + simulation log, including the route generator's arrow diagrams), bench.map, bench.obj, bench.img, bench-sim.log, bench-sim.core, sim_stats.json, and the mdlj_artefacts/mdlj1.yeti and pingorigin_artefacts/pingorigin2.yeti listings. make's own stdout is one PASS/FAIL line per configuration.

Three build steps are specific to this app. generate_bins runs the route generator as ../3dmd_route3sdr NX NY NZX NZY, then python3 ../make-inputbin.py NX NY NZX NZY NLOCMAX params.bin posin.bin velin.bin. The Makefile does not compile 3dmd_route3sdr.cc: after changing it, rebuild 3dmd_route3sdr yourself in the project directory with a C++ compiler, g++ -O2 -o 3dmd_route3sdr 3dmd_route3sdr.cc. A copy of the project that has only the .cc needs the same command before its first run. bench.img is bench.obj patched by miset: core.params/core.fparams from params.bin/fparams.bin, ping.params/ping.fparams from the .origin single-tile variants, router/switch from the generated route tables, and core.pos0/core.vel1 from posin.bin/velin.bin — with the order= strings scraped out of bench.map by grep/sed. out-sim.bin runs miget against bench-sim.core for pos0→posout.bin, vel1→velout.bin, plus ptrs, firstblock, countblock, rinv, df, dist2, tstamp, ilist, jlist, dxij, evec, kevec and ff1.

Correctness checking

Two independent checks run in checkbins:

• check-exchange.py rebuilds the halo bookkeeping in numpy. It un-folds the physical 2D arrays back to a logical nz × ny × nx grid (remap23, using zidx and chain_jtox), recomputes every {xm,xp,ym,yp,zp}count with np.roll over the periodic cell grid, derives the *first offsets and all eight region pointers, and compares them element-by-element against what the device computed. Any mismatch raises. The device's values come from its last rebuild and the expected ones from the final nlocal; they agree because the last inner step does no drift, so its migration moves no atom. This validates the communication pattern, not the physics.
• ecalc.py is the physics oracle. It hard-codes NLOCMAX = 20, the ring-buffer depth 128 and 10 frames, and takes the core count from ptrs.bin. It sums per-PE potential (evec) and kinetic (kevec) energy over all PEs for each recorded frame after the first, turns the three 16-bit timer words into 48-bit counts and prints PE 0's cycles since the first timestamp, and compares (potential, kinetic, total) against the 10 reference rows of data1.dat with max|Δ| / (1 + max|ref|) ≤ 1e-5. It also prints the kinetic energy of the final velocities in velout.bin (keend), but does not check it. Non-finite results raise explicitly.

posbin2dump.py NLOCMAX fparams.bin PARAMS POS VEL FORCES out.dump writes per-PE atoms as one id-sorted LAMMPS dump, recovering type from the sign of the stored id. It is not run by the Makefile, and it does not take the Makefile's outputs as they are: it reads each PE's atom count from the nlocal field of 13-word paramstruct records (the layout of params.bin, which holds the counts from before the run), not from ptrs.bin, and it needs FORCES to hold NLOCMAX atoms per PE, while ff1.bin holds NMAX.

References

• Parameters
• Array attributes
• Hardware register access
• Parallel blocks
• Accumulate assignment
• control()
• Raw wavelets
• Paint code attachment

Run it

The experimental SDK installs this project beside the tutorial, as projects/mdlj. Run it like a tutorial: tungsten-tutorials -f projects/mdlj.