molcrafts-molpack 0.1.0

Faithful Rust port of Packmol molecular packing (split off from molrs)
Documentation
# Core Concepts

This chapter defines each abstraction in the crate in one place.
Cross-link to the types for full API details.

## AtomRestraint

An [`AtomRestraint`](crate::AtomRestraint) is a **soft penalty** applied per
atom: `f(x, scale, scale2) -> F` and `fg(x, scale, scale2, g) -> F`.
It contributes to the packing objective and — in all current
implementations — derives from Packmol's `comprest.f90` / `gwalls.f90`.

```text
pub trait AtomRestraint: Send + Sync + std::fmt::Debug {
    fn f (&self, x: &[F; 3], scale: F, scale2: F) -> F;
    fn fg(&self, x: &[F; 3], scale: F, scale2: F, g: &mut [F; 3]) -> F;
    fn is_parallel_safe(&self) -> bool { true }
    fn name(&self) -> &'static str { std::any::type_name::<Self>() }
}
```

The crate ships 14 concrete `*Restraint` structs (one per Packmol
`kind` 2..=15), each holding its own semantically-named geometric
fields. User types `impl Restraint` sit in the same type slot —
there is no `Builtin*` wrapper in the public API. See
[`extending`](crate::extending) for a tutorial.

### Gradient convention

`fg` accumulates the TRUE gradient (∂penalty/∂x) INTO `g` with `+=`.
Do not overwrite: many restraints may touch the same atom. The
optimizer negates for descent.

### Two-scale contract

Packmol convention (mirrored in the port):

- Linear penalties — kinds 2, 3, 6, 7, 10, 11 (box / cube / plane) —
  consume `scale`.
- Quadratic penalties — kinds 4, 5, 8, 9, 12, 13, 14, 15 (sphere /
  ellipsoid / cylinder / gaussian) — consume `scale2`.

Each `impl Restraint` picks one internally. User-defined restraints
may ignore both knobs and use their own stiffness coefficient as an
instance field.

## Region

A [`Region`](crate::Region) is a **geometric predicate** with a signed
distance function:

```text
pub trait Region: Send + Sync + std::fmt::Debug {
    fn contains(&self, x: &[F; 3]) -> bool;
    fn signed_distance(&self, x: &[F; 3]) -> F;
    fn signed_distance_grad(&self, x: &[F; 3]) -> [F; 3] { /* default FD */ }
    fn bounding_box(&self) -> Option<Aabb> { None }
}
```

Regions compose via the zero-cost combinators
[`And`](crate::And) / [`Or`](crate::Or) / [`Not`](crate::Not), with
analytic chain-rule gradients (max / min / negate). The
[`RegionExt`](crate::RegionExt) trait gives every `Region` ergonomic
`.and(...)` / `.or(...)` / `.not()` methods.

Any `Region` lifts to a `Restraint` via
[`RegionRestraint<R>`](crate::RegionRestraint):

```text
penalty(x) = scale2 * max(0, signed_distance(x))²
```

Use `Region` when you want compositional geometry (intersection /
union / complement). Use `Restraint` directly when you want a specific
penalty shape (linear vs quadratic, custom stiffness, multi-atom).

## Relaxer

A [`Relaxer`](crate::Relaxer) modifies a target's **reference geometry**
between outer optimizer calls. Use cases: torsion-MC sampling for
flexible chains, local MD relaxation, gradient descent on bond-angle
targets.

Two-part design — builder + runner:

- `Relaxer::spawn(&self, frame, ref_coords) -> Box<dyn RelaxerRunner>` is
  called once at `pack()` entry.
- `RelaxerRunner::on_iter(&mut self, coords, f_current, evaluate, rng)`
  runs between movebad and GENCAN each outer iteration; returns
  `Some(new_coords)` on accept, `None` on reject.

Relaxers require `count == 1` because all copies share the same
reference coords.

Built-in: [`TorsionMcRelaxer`](crate::TorsionMcRelaxer) (Metropolis
torsion sampling with self-avoidance).

## Handler

A [`Handler`](crate::Handler) is an observer invoked at well-defined
lifecycle points:

```text
pub trait Handler: Send {
    fn on_start        (&mut self, ntotat, ntotmol)       {}
    fn on_initialized  (&mut self, sys: &PackContext)     {}
    fn on_step         (&mut self, info: &StepInfo, sys);   // required
    fn on_phase_start  (&mut self, info: &PhaseInfo)      {}
    fn on_phase_end    (&mut self, info, report: &PhaseReport) {}
    fn on_inner_iter   (&mut self, iter, f, sys)          {}
    fn on_finish       (&mut self, sys: &PackContext)     {}
    fn should_stop     (&self) -> bool                    { false }
}
```

Observer contract: `sys` is always `&PackContext`, never `&mut`.
Handlers cannot modify packer state — use a `Relaxer` if you need to.

Built-ins: [`NullHandler`](crate::NullHandler),
[`ProgressHandler`](crate::ProgressHandler),
[`EarlyStopHandler`](crate::EarlyStopHandler),
[`XYZHandler`](crate::XYZHandler).

## Objective

The [`Objective`](crate::objective::Objective) trait abstracts over
what GENCAN sees. `PackContext` implements it; synthetic test
objectives (Rosenbrock / Booth / Beale) can implement it to exercise
the optimizer in isolation.

```text
pub trait Objective {
    fn evaluate(&mut self, x: &[F], mode: EvalMode, g: Option<&mut [F]>) -> EvalOutput;
    fn fdist(&self) -> F;
    fn frest(&self) -> F;
    fn ncf(&self) -> u32;
    fn ncg(&self) -> u32;
    fn reset_eval_counters(&mut self);
    fn bounds(&self, l: &mut [F], u: &mut [F]);
}
```

GENCAN (`pgencan`, `gencan`, `tn_ls`, `spg`, `cg`) takes `&mut dyn
Objective` rather than `&mut PackContext` — the optimizer is
decoupled from the packing state.

## Target

A [`Target`](crate::Target) describes one molecule type:

- Input coordinates + centered reference coordinates.
- Van der Waals radii, element symbols, copy count, name.
- Its attached restraints (per-target + per-atom-subset).
- Its attached relaxers.
- Optional fixed placement (Euler + translation).
- Optional Euler-angle bounds (`with_rotation_bound(Axis, Angle, Angle)`).

Targets are snapshotted at `pack()` entry — mutating a `Target` after
passing it to the packer has no effect.

## Molpack

[`Molpack`](crate::Molpack) is the builder facade:

```text
Molpack::new()
    .with_log_level(...)
    .with_handler(...)
    .with_global_restraint(...)  // broadcast to every target
    .with_periodic_box(min, max) // or via periodic InsideBoxRestraint
    .pack(&[targets], max_loops)
```

Every tuning knob (`with_tolerance`, `with_precision`,
`with_inner_iterations`, `with_seed`, `with_avoid_overlap`, …) has a
Packmol-matching default, so `Molpack::new().pack(&targets, max_loops)`
is a complete call. You only set a knob to *change* its default — e.g.
`with_avoid_overlap(false)` to let solvent seed inside a fixed solute
(on by default), or `with_seed(n)` to pick a different RNG stream (the
default seed is Packmol's `1_234_567`).

Every setter consumes and returns `self`. `pack` takes `&mut self`
(handlers are invoked through it).

## PackContext

[`PackContext`](crate::PackContext) is the single owner of mutable
packing state — coordinates, cell lists, restraint pool, rotation
buffers, counters. All optimizer / movebad / handler code paths take
`&mut PackContext` (for writers) or `&PackContext` (for observers).

Structure (`molpack/src/context/`):

- `ModelData` — topology and inputs (immutable after init).
- `RuntimeState` — mutable per-iteration state (x, coor, radius).
- `WorkBuffers` — scratch arrays (xcart, gxcar, radiuswork).

Users rarely touch `PackContext` directly — it's passed through
handlers and relaxers. Power users implementing a custom `Objective`
against synthetic test problems will interact with it.

## Scope equivalence law

```text
molpack.with_global_restraint(r)
    ≡  for t in targets { t.with_restraint(r.clone()) }
```

There is no separate "global-restraint" storage path in `PackContext`.
The broadcast happens inside `pack()`; each target receives an
`Arc::clone` of every global restraint (refcount bump, not a deep
copy).

Per-atom-subset scope is a method-argument pair
`(indices: &[usize], restraint: impl Restraint)`, not a wrapper
struct. There is no `AtomRestraint` public type.

## Restraint versus Constraint

- **Restraint** = soft penalty (violable; pays energy cost).
- **Constraint** = hard constraint (must satisfy; SHAKE / RATTLE /
  LINCS / Lagrange-multiplier mechanisms).

Packmol implements all 15 "constraints" as soft penalties
(`scale * max(0, d)` or `scale2 * max(0, d)²`). Honest naming ⇒
`Restraint`. This crate does not currently define a `Constraint`
trait; adding hard constraints is future work.

## Direction-3 extension pattern

Every extension trait in this crate follows the same shape:

1. Public trait: `pub trait X`.
2. N concrete `pub struct` types that `impl X`, each holding its own
   semantically-named fields.
3. User types `impl X` identically — zero type-level distinction from
   built-ins.

Forbidden in the public API:

- `Builtin*` / `Native*` / `Packmol*` prefixed wrapper types.
- Tagged-union enums that package N built-ins as a single exposed
  type.
- Builder pattern (`X::new().add(...).add(...)`).
- Injection of composition operators into the main trait
  (composition lives on separate traits, e.g. `Region` vs `Restraint`).
- Wrapper types for per-atom-subset scope (that's a method-argument
  pair, not a type).

If you need internal AoS performance structures (e.g. tagged unions
for hot-path match dispatch), they go behind `pub(crate)` and opt-in
via a crate-private hook — invisible to users.