Architecture¶
Developer-oriented view of the crate. Read concepts
first for the abstractions this chapter assumes.
This page covers four things, in order:
- Module map — where everything lives
- Data flow — how values travel from user input to packed frame
- Algorithms — pseudo-code for the three nested loops
- Hot path — what one objective evaluation does
- Invariants and conventions — load-bearing rules
Module map¶
src/
├── lib.rs the public surface (one path per item) + rustdoc chapters
├── pipeline/ the lifecycle body — the only one in the crate
│ ├── mod.rs Pipeline + the five-part run() (validate → space/state →
│ │ chain check → run each stage → assemble)
│ ├── engine.rs StageFactory + PackEngine traits, EngineSetup
│ ├── combinators.rs Repeat (Until) + Guarded (OnViolation) — combinators
│ │ are Stage impls, so mod.rs needs no branch for either
│ └── bracket.rs the callback bracket: adopted-callback tagging, open_bracket, close_bracket
├── settings.rs PackSettings + LogSpec ┐ shared by every run — no
├── pack_space.rs density / pbc / cell resolution │ lifecycle, never import
│ + restraint broadcast │ pipeline/
├── state.rs State + the verbatim Placements ┘
│ a run hands the next one
├── target.rs Target — molecule type + per-molecule restraints + fixed_from
├── template.rs coord_rows — Frame::coords as [x, y, z] rows (Å); rotatable_bonds
│ — the one unclassed-bond policy (crate-root leaf over molrs)
├── restraint/ AtomRestraint trait; region.rs (RegionRestraint over molrs Region),
│ cell.rs (CellRestraint), collective/, crate-private geometric/ (.inp kernels)
├── callback.rs Callback trait + LogLevel + 4 built-in callbacks
├── objective.rs Objective trait + its PackSystem impl over the crate-private
│ compute_f / compute_g / compute_fg
├── eval.rs EvalMode / EvalOutput
├── outcome.rs StageOutcome
├── system/ PackSystem = single owner of mutable packing state
│ ├── pack_system.rs
│ ├── pack_state.rs PackState — system + Placed + RigidView; evaluate_unscaled
│ ├── rigid_view.rs RigidView — the 6·ntotmol COM + Euler placement vector
│ ├── grid.rs install_resolved_cell / coverage_radmax (diameter)
│ ├── geometry.rs GeometryKey — evaluation-cache identity
│ ├── build.rs build_system — system + CSR restraint pool from Targets
│ └── work_buffers.rs scratch arrays (xcart, gxcar, …)
├── pack/ rigid-body family (crate-private)
│ ├── gencan/ GencanPack (gencan_pack.rs), GencanStage, pgencan / gencan / linesearch
│ ├── initial.rs initial random placement
│ ├── restmol.rs single-molecule constraint pre-fit
│ └── movebad.rs worst-molecule perturbation heuristic
├── stage.rs Stage seam — the interface every packing algorithm implements
├── invariant.rs Layers (L0–L5 repair-cost ladder) + Invariant trait +
│ Violation + RestraintsSatisfied — consumed by combinators.rs::Guarded
├── grow/ chain-growth path, peer of the GENCAN path
│ ├── cbmc_grow.rs CbmcGrow — the chain-growth engine (honest verdicts)
│ ├── lattice/ LatticeStage — diamond-lattice SAW for melt density
│ │ (lattice_grow.rs LatticeGrow engine, saw.rs walk,
│ │ decorate.rs backbone on sites, config.rs leaf)
│ ├── config.rs GrowConfig (leaf — no target/engine imports)
│ ├── prior.rs TorsionPrior / AnglePrior + C∞ calibration
│ ├── internal.rs template bond graph → internal-coordinate tree
│ ├── field.rs OverlapField — cell-listed hard-core / soft-shell probe
│ ├── driver.rs GrowStage round loop (seeding, retraction, softening)
│ ├── moves.rs propose / commit / retract / relax primitives
│ └── error.rs GrowError, including lattice-only variants
├── optimizer/ in-loop conformation optimizers (outside the pack family)
├── euler.rs Euler angles ↔ rotation matrices (shared; crate-private)
├── random.rs molpack's own RNG streams (kept for bit parity)
├── assemble.rs packed coords + targets → topology-complete Frame
├── script/ .inp parser + lowering to Targets
└── bin/molpack/ CLI front-end (cli feature)
Dependency direction¶
lib.rs
│
pipeline/ (the lifecycle — depends on everything below)
│
┌──────────┬────────┬───┴─────┬───────────┬──────────┐
▼ ▼ ▼ ▼ ▼ ▼
settings target pack/ grow/ callback
pack_space │ (gencan, (lattice
state │ initial, is a peer
│ │ movebad) inside)
▼ ▼ │ │
└───────────► system/PackSystem (+ grid)
│
▼
objective.rs ← hot path
│
└── eval.rs (EvalMode / EvalOutput)
target / template / restraint / region are pure data — no driver imports.
template.rs owns coord_rows (Å) and rotatable_bonds (an unclassed
bond is a rotatable single bond — for growth and the torsion optimizer alike).
Bond graphs are molrs::core::Topology;
molpack does not ship a parallel Topology type.
pipeline/ is the only module that imports everything else; settings,
pack_space and state import nothing from pipeline/ — the arrow points
one way, pipeline/ reads them, never the reverse.
objective is the narrow waist through which all per-atom work flows.
The chain-growth path enters at the same level as pack: a preset's
stages() (StageFactory) builds one Stage from
the stage seam, Pipeline::run drives it, and
grow/ (the GrowStage) consumes the same PackState and is judged by the
same objective — it never calls the GENCAN internals.
A stage is one interchangeable packing algorithm behind four methods:
name, requires, guarantees, and run. The middle two declare the
shape of the state the stage needs on the way in and promises on the way
out, as a Placed marker — Placed::None (nothing placed yet) or
Placed::All (every free molecule has a placement). run receives a
PackState: the run's PackSystem, plus that marker, plus the rigid
placement vector RigidView (three centre-of-mass and three Euler values
per free molecule). run returns Result<StageOutcome, PackError>: a stage
that cannot do its job fails with a named error, never by disguising failure
as converged = false. On the Ok side, StageOutcome is only what
the stage alone knows — whether it met its own convergence criterion, and
how many times it had to relax a constructive guarantee. The run's
fdist / frest verdict is read off the state afterwards, so no
algorithm grades its own paper.
A Pipeline composes several stages behind one
call — Pipeline::new().with_stage(CbmcGrow::new(prior)).with_stage(GencanPack::new()).run(..)
runs chain growth, then rigid-body push-off, in one lifecycle, continuing
from the first stage's placements rather than re-placing from scratch. Still
out of scope: parallel or branching stage graphs (v1 is a linear sequence
only). For two genuinely independent packs, run each to completion and hand
the earlier result to the next as a fixed obstacle via
Target::fixed_from(&result.frame).
Data flow¶
USER INPUTS ─→ Target / PackEngine builders
Frame, count, restraints, (GencanPack | CbmcGrow | Pipeline)
callbacks, tolerance, seed
─→ PackEngine::run() (one line per preset)
─→ Pipeline::run() the lifecycle body
a. broadcast global → per-target restraints
b. snapshot every Target
c. build PackSystem, wrap into PackState
WorkBuffers (xcart, gxcar, scratch)
d. flatten restraints → CSR pool
e. per stage: Stage::run(state, targets, …)
GENCAN: install grid, then initial
placement — or continue from a seed
PER-ITERATION ─→ evaluate(x, mode, &mut g)
(inside a stage's own loop — → expand_molecules: x → xcart
GencanStage for the rigid path) → restraint penalties per atom
reads f / g via → cell list + pair penalties
&mut dyn Objective → project gradient back: gxcar → g
returns f_total, fdist, frest
OUTPUT ─→ State
frame, plus converged, fdist,
frest, degraded
Three rules govern this flow:
PackSystemowns mutable state. GENCAN, movebad, callbacks, and the phase driver all take&mut PackSystem(writers) or&PackSystem(observers). No other module owns mutable state across iterations.Arc<dyn Restraint>for polymorphic storage. Cheap clone (refcount bump) into the per-atom CSR pool. The hot path does one virtual call per restraint per atom.- GENCAN is decoupled.
gencan/pgencantakes&mut dyn Objective, not&mut PackSystem. Synthetic objectives (Rosenbrock, Booth, Beale) exercise the optimizer in isolation.
Coordinate layout¶
The optimizer variable vector x packs centers of mass and Euler angles:
Cartesian atom positions xcart: Vec<[F; 3]> of length ntotat are
expanded each evaluation:
where i is molecule type, m is copy index, a is atom index.
Algorithms¶
One lifecycle loop over stages, then three nested loops inside the rigid-body stage.
The lifecycle: Pipeline::run() (one call)¶
fn run(targets, max_loops):
validate inputs (non-empty, valid PBC, atoms > 0)
broadcast settings.global_restraints → each target's molecule_restraints
resolve packing space; build PackSystem, wrap into PackState
check the stage chain (empty list / bad order / a preset's non-default settings)
callbacks.on_start
for stage in stages: // one stage for a preset's own run
state.invalidate_geometry_cache()
callbacks.on_stage_start
outcome := stage.run(state, targets, budget, callbacks)? // named PackError
// skips on_stage_end + on_finish
state.set_placed(stage.guarantees().placed)
callbacks.on_stage_end
if callbacks.should_stop(): break
rebuild xcart from the final rigid view; callbacks.on_finish
assemble Frame into State (+ converged / fdist / frest / degraded / intra)
Every preset's PackEngine::run is one line —
Pipeline::single(self).run(targets, max_loops) — so GencanPack::run(),
CbmcGrow::run() and LatticeGrow::run() all resolve to the loop above. It
lives once, in src/pipeline/mod.rs, never duplicated per engine.
Outer: GencanStage::run() (one stage)¶
fn run(state, targets, budget, callbacks):
if state already placed, or this stage carries a seed:
install the box + cell grid
if this stage carries a seed:
inject the seed's placements verbatim; state.placed := All
push_off := state.placed == All
if push_off:
write_xcart(state) // continue from what is there
else:
run init_passes of restmol(): // geometric pre-fit, no pair kernel
for each free target type:
place molecules randomly inside their restraints
relax restraint penalties only
callbacks.on_initialized
for phase in 0 ..= ntype:
if phase < ntype:
comptype[i] := (i == phase) // PER-TYPE pre-compaction
else:
comptype[i] := true // ALL-TYPES main phase
report := run_phase(phase, max_loops, …)
if report.error_phase: break
return Ok(StageOutcome::new(converged, 0))
The preamble — box/grid install, seed injection, the initial()-vs-push-off
choice — and the phase loop both live in GencanStage::run
(src/pack/gencan/solver.rs); nothing above the stage boundary decides when
initial() (and the movebad heuristic it configures) runs. CbmcGrow's
GrowStage (src/grow/driver.rs) is a peer stage under the same lifecycle,
with the phase loop above replaced by its own growth round loop.
Why per-type pre-compaction first: if every type optimizes simultaneously from a random start, cross-type interference traps the solver in shallow minima. Compacting one type at a time inside its own restraint volume gives the all-types phase a much better seed.
Middle: run_phase (one phase)¶
fn run_phase(phase_id, max_loops):
callbacks.on_phase_start(phase_info)
radscale := discale // start with inflated radii (default 1.1)
// Quick-exit: if the unscaled objective is already below precision,
// skip the whole phase.
if evaluate_unscaled(sys, x).below(precision): return Converged
for loop_idx in 0 .. max_loops:
result := run_iteration(loop_idx, radscale, optimizer_bindings)
radscale := decay(radscale) // → 1.0 over the phase
callbacks.on_step(step_info, sys)
if result.converged: return Converged
if callbacks.should_stop(): return EarlyStop
return MaxLoops
radscale starts at discale (1.1) and decays toward 1.0 over the
phase. This soft-starts the pair penalty: the optimizer first sees
slightly oversized atoms (easier to push apart) and tightens to true
tolerance as the phase progresses.
Inner: run_iteration (one outer step)¶
fn run_iteration(loop_idx, radscale, optimizer_bindings):
// 1. Movebad — relocate the K worst molecules.
if movebad enabled:
identify atoms with largest restraint + pair penalty
perturb their COM/Euler within init_box_half_size
// 2. In-loop optimizers — all-type phase only, so that
// COM/Euler indexing covers every molecule.
for binding in optimizer_bindings:
assemble a Frame per selection (moving copies + frozen neighbours)
binding.optimizer.minimize(&mut frame)
map the displacement back into each copy's reference conformer
revert the copy if the packing objective got worse
// 3. GENCAN — bound-constrained quasi-Newton solve.
pgencan(x, &mut sys, params, precision)
// Internally: tn_linesearch → CG inner solve → SPG fallback,
// each step calls sys.evaluate(x, mode, g).
// 4. Convergence check on the unscaled objective.
f_unscaled := evaluate_unscaled(sys, x)
fimp := percentage improvement vs previous loop
converged := fdist < precision AND frest < precision
return { converged, fimp, fdist, frest }
GENCAN itself runs three nested solvers:
pgencan: project x onto bounds, then call gencan
gencan: truncated-Newton outer; calls tn_linesearch
tn_ls: conjugate-gradient line search; SPG fallback if CG stalls
Each leaf step calls sys.evaluate(x, mode, &mut g) — the hot path.
Hot path: objective evaluation¶
Objective::evaluate on PackSystem is invoked O(10³–10⁴) times per run().
Performance lives here.
evaluate(x, mode, g) dispatches by mode:
FOnly → compute_f
GradientOnly → compute_g
FAndGradient → compute_fg
compute_fg is the canonical path — it does five steps:
1. expand_molecules(x):
for each molecule type t, copy m, atom a:
xcart[icart] := com_t,m + R(eul_t,m) · ref_coords[t, a]
2. accumulate_constraint_value_and_gradient (per atom icart):
range := iratom_offsets[icart] .. iratom_offsets[icart + 1]
for &irest in iratom_indices[range]:
f += sys.restraints[irest].fg(xcart[icart], scale, scale2,
&mut grad_xcart[icart])
// Linear penalties consume `scale`; quadratic consume `scale2`.
3. insert_atom_in_cell (per atom):
linked-list bucket atoms into cells
cell side ≈ 2 × max_radius × radscale
4. accumulate_pair_fg (or _parallel under rayon):
for each non-empty cell c:
for each neighbor cell c′ in 13-cell stencil:
for each (i ∈ c, j ∈ c′):
d := pbc_distance(xi, xj)
σ := (rᵢ + rⱼ) · radscale
if d < σ:
penalty := (σ − d)²
grad_xcart[i] += d penalty / d xi
grad_xcart[j] += d penalty / d xj
5. project_cartesian_gradient:
for each molecule m, atom a:
g_com[m] += grad_xcart[icart]
g_euler[m] += Jᵀ(eul_m, ref_a) · grad_xcart[icart]
// J = ∂xcart/∂eul, derived once per molecule from R(eul).
Cost breakdown: steps 1–3 are O(N_atoms); step 4 is
O(N_atoms × neighbor_avg) ≈ O(N_atoms × 32) and dominates wall time on
realistic workloads. Step 4 is the rayon parallelization point
(accumulate_pair_fg_parallel), reducing into per-atom gradient slots
via AtomicU64 (since Cell<f64> is not Sync).
The Arc<dyn Restraint> virtual call in step 2 measured at +0.22% e2e
versus the prior monomorphic dispatch — a negligible cost for the
flexibility of user-defined restraints.
Invariants and conventions¶
Gradient accumulation. Restraint::fg accumulates the true
gradient (∂penalty/∂x) into g with +=. Optimizer negates for descent.
Multiple restraints may touch one atom, so never overwrite.
Two-scale contract. Linear penalties (Packmol kinds ⅔/6/7/10/11)
consume scale; quadratic penalties (kinds ⅘/8/9/12/13/14/15) consume
scale2. Each impl Restraint picks one internally.
Rotation convention. R_new = δR · R_old (LEFT multiplication).
Single-atom tests cannot detect LEFT/RIGHT bugs — always test with
≥ 2 atoms.
Coordinate layout. GENCAN's x is [com₀..n, eul₀..n] of length
6·ntotmol. Cartesian atom positions xcart are Vec<[F; 3]> of length
ntotat.
Thread safety. All trait objects are Send + Sync. Interior
mutability inside parallel reductions uses AtomicU64 with
f64::to_bits / f64::from_bits — Cell<f64> is not Sync.
Scope equivalence.
There is no separate global-restraint storage path. The broadcast at
Pipeline::run() entry is the implementation.
Restraint vs Constraint. Packmol implements all 15 "constraints" as
soft penalties. Naming reflects mechanism, not user intent → Restraint.
Direction-3 extension pattern. Every extension trait follows the
same shape: public trait, N concrete pub struct impls, user types
impl Trait identically. No Builtin* / Native* wrappers in the
public API.
init1 short-circuit. Set during the initial geometric pre-fit.
Skips the pair kernel — the restraint-only objective is enough to get
atoms into their regions before pair conflicts matter.
Cheatsheet¶
| Question | Where to look |
|---|---|
| How is one restraint's penalty computed for one atom? | restraint/*::f / *::fg |
Where does with_global_restraint broadcast? |
pack_space.rs::broadcast_global_restraints |
| Where is the per-atom CSR pool built? | system/build.rs::build_system (CSR build loop) |
How are x ↔ Cartesian coords expanded? |
objective.rs::expand_molecules, euler.rs::eulerrmat |
| Where is the pair-overlap kernel? | objective.rs::accumulate_pair_fg_parallel |
| What does the initial pre-fit do? | pack/initial.rs::initial, pack/restmol.rs::restmol |
| How is precision-based termination tested? | pack/gencan/search.rs::converged (Packmol's packmolprecision) |
What does movebad do? |
pack/movebad.rs::movebad |
| How is torsion MC wired in? | optimizer/torsion_mc.rs::TorsionMcOptimizer::minimize, called from optimizer/mod.rs::run_optimizer_bindings |
| Where does periodic boundary wrap apply? | system/pack_system.rs::pbc_distance |