Workflow
Defining binding rules
A BindingRules object holds a list of particle species and the bonds allowed between their binding sites. A bond is written [species_i site_i species_j site_j], meaning site site_i of species species_i may bind site site_j of species species_j.
julia> using Roly
julia> bonds = [1 3 2 3;
2 2 3 2;
2 1 4 1;
3 1 4 1];
julia> rules = BindingRules(bonds, UnitTriangle)
2d BindingRules[n=4, k=4]Pass a single species if all blocks have the same shape, otherwise a vector with one species per index used in bonds.
In 2D Roly ships UnitTriangle, UnitSquare, UnitHexagon and PolygonParticleSpecies for regular polygons, and PatchyDisk for a disk with sites on its rim. In 3D it ships UnitTetrahedron, UnitCube, UnitOctahedron, UnitDodecahedron, UnitIcosahedron, UnitPyramid, UnitPrism, UnitAntiprism and PolyhedronParticleSpecies for any convex polyhedron with one site per face, and PatchySphere for a sphere whose patches inherit a polyhedron's rotation group. See Custom particle species to define your own.
Colors and symmetry
You describe a species by coloring its binding sites, and the bond table refers to those colors. Which sites the particle cannot tell apart follows from the coloring and the geometry, and Roly derives it.
julia> symmetrynumber(PolyhedronParticleSpecies(Cube())) # every face distinct
1
julia> symmetrynumber(PolyhedronParticleSpecies(Cube(); colors=fill(1, 6))) # all faces alike
24
julia> caps = [abs(Roly.facenormal(Cube(), i)[3]) > 0.5 ? 2 : 1 for i in 1:6];
julia> symmetrynumber(PolyhedronParticleSpecies(Cube(); colors=caps)) # caps apart from sides
8Faces come in no particular order, so the third example picks the caps by their normals. The answer 8 is D_4: telling two opposite faces apart leaves the 4-fold axis through them and the 2-fold axes across it.
Two more keywords say what a bond at a face means, rather than which bonds exist. locking decides whether a site holds its partner in the orientation its frame names, and twists turns a site about its normal to pick which orientation that is. See Orientation and twists.
Symmetry groups
Bodies are built by name: Cube, Prism, Antiprism and the rest all return a Polyhedron. rotationgroup lists the rotations a body has, and grouporder counts those of a named group.
julia> length(Roly.rotationgroup(Cube()))
24
julia> grouporder(Octahedral())
24
julia> grouporder(Dihedral(5))
10The named groups are Cyclic, Dihedral, Tetrahedral, Octahedral and Icosahedral, the rotation-only point groups a rigid body can have.
Sketching rules interactively
ruleeditor builds a bond table geometrically: place blocks on a lattice and it reads the rules off every pair of touching sites.
rules = ruleeditor(UnitSquare) # also works for UnitTriangle and UnitHexagonArrow keys move the cursor, Enter places, Space erases, r and R rotate, digits 1 to 9 switch species, q accepts.
┌ Editor ──────────┐┌ Particles ──────────────────────────┐
│── Species ── ││ ┌─────┐ ┌─────┐ │
│▶ ■ species 1 ││ │ │ │ │ │
│ ■ species 2 ││ │ → │ │ ← │ │
│ ││ │ │ │ │ │
│── Keys ────── ││ └─────┘ └─────┘ │
│arrows move ││ ┌─────┐ ┌─────┐ │
│enter place ││ │ │ │ │ │
│space erase ││ │ ↑ │ │ ↑ │ │
│r / R rotate ││ │ │ │ │ │
│1-9 species ││ └─────┘ └─────┘ │
│c clear ││ │
│q accept ││ │
└──────────────────┘└─────────────────────────────────────┘Pass output=:bonds or output=:matrix for a copy-pasteable result instead of a BindingRules:
bonds = ruleeditor(UnitSquare; output=:bonds) # n×4 integer matrix
rules = BindingRules(bonds, UnitSquare) # reproduces the same rulesEnumerating polyforms
polyenum walks every polyform the rules allow.
julia> result = polyenum(rules; maxsize=20, maxstrs=100_000);
julia> result.nstructures
16
julia> result.largest_size
5
julia> result.status
Finished::RSStatus = 0Cap maxsize (particles per polyform) or maxstrs (total polyforms) when the rules allow unbounded growth. status says why the run stopped: Finished, MaxDepthReached, MaxVerticesReached or BreakTriggered.
Storing polyforms
polygen returns the polyforms in a Vector, sorted by size.
julia> polys = polygen(rules; maxsize=20);
julia> length(polys)
16Counting polyforms
countpolyforms counts without storing anything, switching to an unbiased sampled estimate when exact enumeration gets too expensive.
julia> c = countpolyforms(rules);
julia> c.n
16.0
julia> c.exact
true
julia> c.uncertainty
0.0It requires an explicit maxsize when the rules allow polyforms of unbounded size.
Applying constraints
polyenum takes a callback that runs at each polyform, receiving it and its size and returning one of three signals:
ACCEPTcounts the polyform and keeps exploring it,REJECTskips the polyform and everything grown from it,BREAKstops the enumeration.
julia> constraint(s, _) = composition(s)[4] <= 1 ? ACCEPT : REJECT;
julia> polyenum(constraint, rules).nstructures
14REJECT prunes a whole subtree, so the constraint must be monotone. If a polyform violates it, everything grown from it must violate it too.
Visualizing polyforms
Load any Makie backend to activate the plotting extension.
using GLMakie # or CairoMakie, WGLMakie, ...
render(polys[end]) # display a single polyformrender picks a 2D or 3D axis to match the polyform. Use GLMakie or WGLMakie for 3D, since CairoMakie sorts primitives instead of depth-testing them and shows seams where faces meet. 2D output is fine in any backend.
A species renders too, which is the quickest way to see how its faces are colored and which way its sites face.
render(PolyhedronParticleSpecies(Prism(3); colors=[1, 2, 2, 2, 1]))
render(UnitCube; bindingrules=rules) # sites no bond can use are drawn inertpolyformplot!(ax, poly) draws onto an existing Makie axis.