API reference

AssemblySystem

Wraps a set of Roly.jl BindingRules together with an EntropyModel used to compute structure partition functions. Structure enumeration, the composition matrix, and partition functions are computed lazily and cached on first use.

AssemblySystem(bindingrules::BindingRules, entropymodel::EntropyModel; maxsize=Inf, verbose=true)

Construct an AssemblySystem from bindingrules and an EntropyModel used to compute structure partition functions. Only structures of at most maxsize particles are considered.

compositionmatrix(asys::AssemblySystem; enum_kwargs...)

Return the composition matrix of asys, computing and caching it on first use. Unless the structures themselves are already cached, they are not retained.

countstructures(asys::AssemblySystem; kwargs...)

Count how many structures asys would enumerate, without enumerating them. Cheap enough to run before committing to structures or partitionfunctions on a large system.

Returns a StructureCount, exact when the enumeration fits in exact_budget and a statistical estimate otherwise — as opposed to nstructures, which returns the plain number of structures actually enumerated and pays for the enumeration to do it.

entropydimension(asys::AssemblySystem)

The dimension in which the entropy model of asys scores its structures: 3 whenever the model embeds, and dimension(asys) otherwise.

iscomplete(asys::AssemblySystem)

Whether asys covers every structure its rules allow, i.e. whether maxsize cut the enumeration short. missing if nothing has been enumerated yet.

nstructures(asys::AssemblySystem; enum_kwargs...)

Return the number of enumerated structures of asys, without retaining the structures themselves.

partitionfunctions(asys::AssemblySystem; enum_kwargs...)

Return the per-structure partition functions of asys, computed via entropy using asys.bondpotential, computing and caching them on first use.

structures(asys::AssemblySystem; enum_kwargs...)

Return the enumerated structures of asys, computing and caching them on first use.

BondPotential

Supertype of the bond potentials. See RigidSpringPotential.

RigidSpringPotential(σ=1; k=1)

Harmonic springs of total stiffness k between the patch clouds of two bonded binding sites.

Two dimensional, with the patches spread uniformly along a segment of length σ across the site, which has second moment σ²/12.

RigidSpringPotential

A bond modelled as harmonic springs between two rigid clouds of patches, one carried by each of the two bonded binding sites,

V = -ε + Σ_α kα/2 ‖pα₁(q₁) - pα₂(q₂)‖² ,

with the patch α of site i sitting at pαᵢ(qᵢ) = xᵢ + R(ψᵢ) aαᵢ. Call it as potential(pose1, pose2, color1, color2) on the lab poses of the two binding sites; the binding energy ε is not included, as it is carried separately by the reaction network.

Only the total spring constant k = Σ kα and the patch (cross-)covariances Sᵢⱼ = k⁻¹ Σ kα aαᵢ aαⱼᵀ enter, so those are what is stored. Since the poses are those of the binding sites rather than of the particles, and bonded sites coincide, every cloud is taken to be centred on its own site: patch positions are used only through their spread aα - ā.

See bondvolume for the configurational integral of a single such bond.

RigidSpringPotential(A::AbstractMatrix; k=1)

Harmonic springs of total stiffness k between the patch clouds of two bonded binding sites.

A is an explicit d x n cloud of patch positions in the binding site's frame, carried by every site; k is either a number or a vector of the per-patch stiffnesses.

RigidSpringPotential(As::AbstractVector{<:AbstractMatrix}; k=1)

Harmonic springs of total stiffness k between the patch clouds of two bonded binding sites.

One d x n cloud of patch positions per binding site color, As[c] in the frame of color c, paired column by column; k is either a number or a vector of the per-patch stiffnesses. The patches of two colors are in general not compatible, unless their clouds are related by the rotation that brings the two sites face to face. If incompatible, the bond carries a strainenergy.

RigidSpringPotential(σx, σy; k=1)

Harmonic springs of total stiffness k between the patch clouds of two bonded binding sites.

Three dimensional, with the patches spread uniformly over a σx x σy sheet across the site, which has second moments σx²/12 and σy²/12.

bondenergy(potential, pose1, pose2, color1, color2)

Energy of one bond, given the lab poses and colors of the two binding sites that form it.

bondvolume(potential, color1=1, color2=1)
logbondvolume(potential, color1=1, color2=1)

The configurational integral ω(S, k) = exp(-∑ε) ∫ exp(-V(g)) dg of a single bond, over the relative pose g ∈ SE(d) of the two binding sites that form it, excluding the binding energy.

checkpotential(potential, rules::BindingRules, d::Integer)

Check that potential can describe every bond rules allows, in d dimensions, and throw if it cannot.

contactexcess(potential, color1=1, color2=1)

How far above its own minimum a bond sits when its two binding sites are placed as Roly bonds them, k [Σ_μ σ̃_μ(S₁₂) - tr(S₁₂ R⋆ᵀ)]. Zero exactly when the potential is minimised at that pose.

Distinct from strainenergy, which asks only whether the springs can relax somewhere. A potential is free to relax elsewhere and bondvolume will still be right, since it integrates over all poses; the Laplace solvers expand about the bonded pose and so reject it.

bondvolume(potential, color1=1, color2=1)
logbondvolume(potential, color1=1, color2=1)

The configurational integral ω(S, k) = exp(-∑ε) ∫ exp(-V(g)) dg of a single bond, over the relative pose g ∈ SE(d) of the two binding sites that form it, excluding the binding energy.

strainenergy(potential, color1=1, color2=1)

Energy k [tr(S₁₁ + S₂₂)/2 - Σ_μ σ̃_μ(S₁₂)] left in the springs of a bond that cannot relax them all at once, in the stiff-spring limit. Non-negative, and zero exactly when the two patch clouds are compatible.

Note the sign against the binding energies εs, which are bond strengths: a bond contributes exp(ε - strainenergy) to a density, so strain has to be subtracted from ε to fold it in, not added.

COMLaplace()

Entropy solver that expands the bond potential to second order about the bonded configuration, working in the frame of the structure's centre of mass.

EntropyModel(potential, solver=COMLaplace(); embed3d=false)
EntropyModel(solver=TreeLike(); embed3d=false)

Everything needed to turn a structure into a partition function: the bond potential, the solver used to evaluate it, and whether to embed the structure in 3d. The second form is for solvers that take no potential.

EntropySolver

Supertype of the methods for turning a structure and a bond potential into a partition function. See TreeLike, MeanField, TetheredLaplace and COMLaplace.

MeanField()

Entropy solver that is exact at mean field level for a RigidSpringPotential, where every particle sees its neighbours' patches at their thermally averaged positions rather than at their instantaneous ones.

Needs an unstressed structure, so that those averaged positions are the ones the particles are glued to; each particle then contributes ω(kᵢ, ⟨Sᵢ⟩) built from the spread of its own patches, and the structure's topology enters only through how many patches each particle carries.

Being variational it can only undercount. A dimer it gets exactly right, since tethering one of two particles leaves a single free one and the product ansatz is then no ansatz at all.

TetheredLaplace()

Entropy solver that expands the bond potential to second order about the bonded configuration, holding one particle fixed to remove the global translations and rotations.

TreeLike(omega=1)

Entropy solver that treats every structure as a tree, charging n - 1 bond volumes on top of the structure's orientational volume. A structure that really is a tree has exactly n - 1 bonds and the answer is exact. For non-trees, the bond volumes are averaged geometrically over the bonds. With a bond potential the bondvolumes are computed from bondvolume, with no potential each bond contributes omega.

bondcolors(poly::Polyform)

The binding site colors (c1, c2) of every bond of poly, in the order Roly.bonds gives them.

entropy(model::EntropyModel, poly::Polyform)

Partition function of the structure poly under model.

entropydimension(m::EntropyModel, x)

The dimension in which m scores the structures of x, a BindingRules or a Polyform: 3 whenever m embeds, and the dimension of x itself otherwise.

hessian(bond_potential, poly; embed3d=false)

Hessian of map_potential at the bonded configuration, of shape (dtot*nparticles, dtot*nparticles).

map_potential(bond_potential, poly::Polyform; embed3d=false)

Build the function that gives the total bond energy of poly as a function of a dtot x nparticles matrix of displacements from the bonded configuration.

_initialpotentials(ϕs, εs, N, B, Ωs, μrange, T)

Initial guess for the chemical potentials of the species in μrange: the ideal-gas value μ_i = log(ϕ_i) - log(Ω_i), obtained by putting every particle of species i in its own monomer.

Binding only adds to ϕ_i at fixed μ_i, so this is an upper bound on the root, and it is the one place the concentrations enter. They routinely span several decades across species, and a guess that ignores them can leave the solver a very long way from the root: a guess set by the binding energies alone is the same number for every species, however far apart their concentrations are.

Species with no monomer among the structures fall back to -1.1 * median(abs.(εs)).

chemicalpotentials(asys::AssemblySystem, ϕs, εs)

Compute the chemical potentials of each particle species used in asys as a function of particle concentrations (ϕs) and binding energies (εs).

alg selects the nonlinear solver and abstol/reltol its tolerances; any further keyword arguments are passed on to solve.

densities(asys::AssemblySystem, ϕs, εs)

Compute the equilibrium number densities of the structures of asys as a function of particle concentrations (ϕs) and binding energies (εs).

densities(asys::AssemblySystem, ξ)

Compute the equilibrium number densities of the structures of asys as a function of ξ, a vector containing chemical potentials and binding energies.

density_jacobian(asys::AssemblySystem, ϕs, εs)
density_jacobian(ϕs, εs, M, Ωs)

Derivative of the structure number densities with respect to particle concentrations (ϕs) and binding energies (εs).

density_jacobian(asys::AssemblySystem, ξ)
density_jacobian(ξ, M, Ωs)

Derivative of the structure number densities with respect to ξ, a vector containing chemical potentials and binding energies.

logdensities(asys::AssemblySystem, ϕs, εs)

Compute the log equilibrium number densities of the structures of asys as a function of particle concentrations (ϕs) and binding energies (εs).

logdensities(asys::AssemblySystem, ξ)

Compute the log equilibrium number densities of the structures of asys as a function of ξ, a vector containing chemical potentials and binding energies.

logparticledensities(asys::AssemblySystem, ϕs, εs)

Compute the log total equilibrium number densities of each particle species used in asys as a function of particle concentrations (ϕs) and binding energies (εs).

logparticledensities(asys::AssemblySystem, ξ)

Compute the log total equilibrium number densities of each particle species used in asys as a function of ξ, a vector containing chemical potentials and binding energies.

logyields(asys::AssemblySystem, ϕs, εs)

Compute the log equilibrium yields of the structures of asys as a function of particle concentrations (ϕs) and binding energies (εs).

logyields(asys::AssemblySystem, ξ)

Compute the log equilibrium yields of the structures of asys as a function of ξ, a vector containing chemical potentials and binding energies.

particledensities(asys::AssemblySystem, ϕs, εs)

Compute the total equilibrium number densities of each particle species used in asys as a function of particle concentrations (ϕs) and binding energies (εs).

particledensities(asys::AssemblySystem, ξ)

Compute the total equilibrium number densities of each particle species used in asys as a function of ξ, a vector containing chemical potentials and binding energies.

particledensity_jacobian(asys::AssemblySystem, ξ)
particledensity_jacobian(ξ, M, Ωs; ns=_nspecies(M))

Derivative of the total particle densities with respect to ξ, a vector containing chemical potentials and binding energies.

topotentials(asys::AssemblySystem, ϕs, εs)
topotentials(ϕs, εs, M, Ωs)

Convert particle concentrations (ϕs) and binding energies (εs) into the ξ vector, consisting of chemical potentials and binding energies.

toyields(densities)

Normalize a list of number densities into yields. If densities has multiple axes, the normalization is carried out over dims.

yield_jacobian(asys::AssemblySystem, ϕs, εs)
yield_jacobian(ϕs, εs, M, Ωs)

Derivative of the structure yields with respect to particle concentrations (ϕs) and binding energies (εs).

yield_jacobian(asys::AssemblySystem, ξ)
yield_jacobian(ξ, M, Ωs)

Derivative of the structure yields with respect to ξ, a vector containing chemical potentials and binding energies.

yields(asys::AssemblySystem, ϕs, εs)

Compute the equilibrium yields of the structures of asys as a function of particle concentrations (ϕs) and binding energies (εs).

yields(asys::AssemblySystem, ξ)

Compute the equilibrium yields of the structures of asys as a function of ξ, a vector containing chemical potentials and binding energies.

Reaction

A single aggregation/fragmentation reaction i + j ⟷ k, as handed to a rate kernel.

A lightweightview into a ReactionNetwork; use indices, reactants, product, bondcounts, cut and halves to read it.

ReactionNetwork(asys::AssemblySystem; maxbonds=Inf, fwdkernel=Returns(1.0), bwdkernel=fwdkernel)

The reactions i + j ⟷ k among the structures of asys that break at most maxbonds bonds at a time, with a forward and a backward rate for each.

Each kernel takes a Reaction and returns a rate, read back with fwdrates and bwdrates. active marks the reactions that ended up with a nonzero rate; the others cannot contribute and may be skipped.

Iterating a ReactionNetwork yields one Reaction per entry. Use rate! to apply different kernels without repeating the cut enumeration, which is the expensive part.

assemblysystem(net::ReactionNetwork)

The AssemblySystem whose structures net reacts.

bondcounts(rxn::Reaction)

How many bonds of each type are formed, equivalently broken, as a vector of length nbonds. Obtained from the composition of the product less those of its fragments; the particle counts cancel, since particles are conserved.

reactions(net::ReactionNetwork)
fwdrates(net::ReactionNetwork)
bwdrates(net::ReactionNetwork)

The (i, j, k) index triples of every reaction, and their forward and backward rates. Positions line up with net[r], with iterating net, and with each other.

These are the network's own storage rather than copies, so they allocate nothing but should be treated as read-only; use rate! to change rates.

cut(rxn::Reaction)

The edges broken to separate the product into its two fragments.

reactions(net::ReactionNetwork)
fwdrates(net::ReactionNetwork)
bwdrates(net::ReactionNetwork)

The (i, j, k) index triples of every reaction, and their forward and backward rates. Positions line up with net[r], with iterating net, and with each other.

These are the network's own storage rather than copies, so they allocate nothing but should be treated as read-only; use rate! to change rates.

generate_cuts(g; maxbonds)

Every way of splitting g into exactly two connected pieces by breaking at most maxbonds bonds. Returns the broken edges and the vertices of each piece.

halves(rxn::Reaction)

The vertices of the product belonging to each fragment.

indices(rxn::Reaction)

(i, j, k): the positions of the two reactants and the product in the structure list.

isdetailedbalanced(net::ReactionNetwork)

Whether every active reaction has equal forward and backward rates.

nreactions(net::ReactionNetwork)

The number of reactions in net, active or not.

product(rxn::Reaction)

The structure that forms, carrying the geometry its fragments have once bound.

rate!(net::ReactionNetwork; fwdkernel=Returns(1.0), bwdkernel=fwdkernel)

Set the rates of the reactions of net in place. The kernels should take a Reaction and return a finite, non-negative rate; anything else is rejected here rather than downstream, where a bad rate is mostly silent. A negative one yields a plausible-looking stability matrix and a DomainError from the kinetics, and a NaN reports itself as a detailed-balance violation, since NaN != NaN.

reactants(rxn::Reaction)

The two fragments, as the canonical structures of their isomorphism class. These describe the same fragments as product restricted to halves, but in their free rather than bound poses.

reactions(net::ReactionNetwork)
fwdrates(net::ReactionNetwork)
bwdrates(net::ReactionNetwork)

The (i, j, k) index triples of every reaction, and their forward and backward rates. Positions line up with net[r], with iterating net, and with each other.

These are the network's own storage rather than copies, so they allocate nothing but should be treated as read-only; use rate! to change rates.

simulate_kinetics(net::ReactionNetwork, ϕs, εs; T, kwargs...)

As above, but as a function of particle concentrations (ϕs) and binding energies (εs).

simulate_kinetics(net::ReactionNetwork, ξ; T, [initial_densities, saveat])

Simulate the assembly kinetics over the reactions of net as a function of ξ, a vector containing chemical potentials and binding energies. Unless initial_densities is specified, the initial state consists of monomers at the concentrations implied by ξ.

Keyword arguments:

• T: total simulation time.
• initial_densities (optional): initial densities of all structures at t=0.
• saveat (optional): time points at which the simulation state should be stored.
• sparsejac (optional): store the Jacobian as a sparse matrix.
• alg (optional): the ODE algorithm, Rodas5P() by default.

Remaining keyword arguments go to solve.

Returns a vector of time points and a matrix of structure number densities of shape (nstructures, ntimepoints).

DensityBasis()

Coordinates are structure number densities, so S = ∂(dρ/dt)/∂ρ as it stands.

StabilityBasis

The basis a stability matrix is expressed in. DensityBasis is the plain one, in structure number densities; SymmetricBasis is the similarity transform that makes it symmetric.

SymmetricBasis()

Coordinates are rescaled by D = diagm(sqrt.(ρeq)), giving inv(D) * S * D. Detailed balance makes that symmetric, so its spectrum is real, and it is the same spectrum as in DensityBasis since the two differ by a similarity transform.

correlationtime(net::ReactionNetwork, ξ; nconserved=nspecies(assemblysystem(net)), scale=1)
correlationtime(S::AbstractMatrix; nconserved)

The slowest relaxation time of the linearised kinetics, -1/λ for the least negative eigenvalue of the stability matrix that is not one of its zero modes.

• nconserved: how many zero eigenvalues to skip, one per conserved particle count.
• rtol: eigenvalues within rtol * maximum(abs, S) of zero count as zero modes.

Returns Inf when the selected mode is itself a zero mode, meaning the kinetics has more conserved quantities than nconserved says.

The network method builds S in the SymmetricBasis, which makes the spectrum real; that spectrum is shared with DensityBasis, so the basis only affects accuracy. It inherits the isdetailedbalanced requirement from stabilitymatrix.

stabilityjacobian(net::ReactionNetwork, ξ, [basis]; scale=1)
stabilityjacobian!(J, net::ReactionNetwork, ξ, [basis]; scale=1)

Derivative of stabilitymatrix with respect to ξ, of shape nstructures x nstructures x length(ξ). Requires isdetailedbalanced, as above.

stabilityjacobian(net::ReactionNetwork, ξ, [basis]; scale=1)
stabilityjacobian!(J, net::ReactionNetwork, ξ, [basis]; scale=1)

Derivative of stabilitymatrix with respect to ξ, of shape nstructures x nstructures x length(ξ). Requires isdetailedbalanced, as above.

stabilitymatrix(net::ReactionNetwork, ξ, [basis]; scale=1)
stabilitymatrix!(S, net::ReactionNetwork, ξ, [basis]; scale=1)

The nstructures x nstructures stability matrix of net about the equilibrium at ξ, expressed in basis (default DensityBasis). scale multiplies every rate, and so the whole matrix.

Requires isdetailedbalanced: the equilibrium densities are only a steady state of the kinetics when each reaction's forward and backward rates agree, so there is nothing to linearise about otherwise.

stabilitymatrix(net::ReactionNetwork, ξ, [basis]; scale=1)
stabilitymatrix!(S, net::ReactionNetwork, ξ, [basis]; scale=1)

The nstructures x nstructures stability matrix of net about the equilibrium at ξ, expressed in basis (default DensityBasis). scale multiplies every rate, and so the whole matrix.

Requires isdetailedbalanced: the equilibrium densities are only a steady state of the kinetics when each reaction's forward and backward rates agree, so there is nothing to linearise about otherwise.

capfraction(δ)

Fraction of a sphere covered by a cap of half-angle δ.

centerofdiffusion(xs, rs)

The point about which the translation-rotation coupling of the cluster is symmetric.

diffusionconstants(xs, rs; cod=nothing)
diffusionconstants(p::Polyform; cod=nothing)

The orientationally averaged translational and rotational diffusion constants, (Dlin, Drot).

diffusiontensor(xs, rs; cod=nothing)

The 6 x 6 diffusion tensor of the cluster about its center of diffusion.

frictiontensor(xs, rs)

The 6 x 6 grand friction tensor of spheres of radii rs at positions xs, held rigidly together. Blocks are translation-translation, translation-rotation and rotation-rotation, about the origin.

orientedbindingrate(rxn::Reaction; siteradius=0.5)

Diffusion-limited association rate when binding also requires the two partners to be correctly oriented. Takes relevant parameters from the geometry reacting particles. siteradius determines the size of the binding patches.

Usable directly as a kernel for rate!.

orientedbindingrate(; D, Dr1, δ1, Dr2, δ2, R)
orientedbindingrate(rxn::Reaction; siteradius=0.5)

Diffusion-limited association rate when binding also requires the two partners to be correctly oriented: reactive caps of half-angle δ1, δ2 on spheres of separation R, with relative translational diffusion D and rotational diffusion Dr1, Dr2.

Reduces to smoluchowskirate as the caps grow to cover their spheres.

rotneprager(Δx; ri=1, rj=1)

Rotne-Prager-Yamakawa pair mobility of two spheres separated by Δx.

Overlapping spheres take Yamakawa's near-field branch. The far-field form is not positive definite there: it returns pair mobilities exceeding the self mobility, which makes the grand mobility matrix indefinite and diffusion constants come out negative. The two branches agree at contact.

smoluchowskirate(rxn::Reaction)

Diffusion-limited association rate of two spheres that react on contact whatever their orientation, 4π D R. Usable directly as a kernel for rate!.

smoluchowskirate(; D, R)

Diffusion-limited association rate of two spheres that react on contact whatever their orientation, 4π D R.

lineardesign(M, idxs; optimizer, kwargs...)
lineardesign(asys::AssemblySystem, idxs; optimizer, kwargs...)

Find parameters that make the target structures idxs the only ones at zero excess free energy.

• idxs: index or indices of the target structures
• optimizer: any MathOptInterface optimizer, e.g. Clarabel.Optimizer. Configure solver tolerances on the optimizer itself, with MOI.OptimizerWithAttributes
• preprocess = true: drop structures built from particles or bonds the targets do not use
• refine = true: if the targets are not designable on their own, optimize for the smallest designable set containing them
• atol = 1e-6: tolerance for deciding that a structure sits at zero excess free energy
• silent = true: suppress solver output
• infval = 100: finite stand-in for infinite parameters

Returns (ξ, residual). Requires Convex to be loaded.

maxyielddesign(M, idxs; omegas, maxdensity, optimizer, kwargs...)
maxyielddesign(asys::AssemblySystem, idxs; maxdensity, optimizer, kwargs...)

Maximize the yield of the target structures idxs within a bond energy budget.

• idxs: index or indices of the target structures
• maxdensity: cap on the total particle density, in the units omegas is expressed in
• relative_yields = nothing: densities to hold the targets at, relative to each other, one entry per target; nothing holds them equal
• omegas: entropic partition function of every structure
• energy_budget = 1: cap on the bond energies, applied to energy_measure of them
• energy_measure = mean: what the budget bounds, e.g. mean or maximum
• uniform_energy = false: hold every bond type at the same energy
• target_stoichiometry = false: keep the particle densities in the targets' stoichiometric ratio
• optimizer: any MathOptInterface optimizer, e.g. Clarabel.Optimizer
• preprocess, silent, infval: as in lineardesign

The budget is a cap, not a target: bond energy also favours the larger off-target structures, so the solution can come in under it.

target_stoichiometry is a cap too, since holding the ratios exactly is not a convex constraint. It caps each species at the targets' ratio, which leaves under-supply allowed; that is expected to bind rather than to relax, and a warning reports the solves where it does not.

Returns (ξ, residual), where residual is log Σ_{j∉idxs} ρ_j/ρ_i. Requires Convex to be loaded.

minenergydesign(M, idxs; minyield, omegas, maxdensity, optimizer, kwargs...)
minenergydesign(asys::AssemblySystem, idxs; minyield, maxdensity, optimizer, kwargs...)

Find the weakest bonds that still reach a yield of minyield for the target structures idxs.

The inverse of maxyielddesign: the yield enters as a constraint and the bond energies become the objective.

• minyield: yield the targets must reach, as a fraction of all structures present
• energy_measure = mean: what is minimized, e.g. mean or maximum of the bond energies
• uniform_energy = false: hold every bond type at the same energy
• target_stoichiometry = false: as in maxyielddesign, and a cap there too
• other arguments as in maxyielddesign, energy_budget excepted: here the energy is the objective, not a constraint

Returns (ξ, residual), where residual is the achieved energy_measure. Requires Convex to be loaded.