mod known_basin

module known_basin

Invert known wells so rgmin does not walk a Leave back into them. Invert basins already occupied so rgmin will not walk back into them.

Henkelman and Jónsson, J. Chem. Phys. 1999, 111, 7010 <https://doi.org/10.1063/1.480097>: the dimer replaces (F) by (F-2(Fcdot P)P) so a first-order stepper walks up the lowest mode to a saddle. Occupancy Leave applies the same Householder in the DECAF packing map, not in Cartesian radius. (mu) is the mean of per-center crate::soap::local_nu3_z (SOAP plus ACE (nu=3)) at crate::catalog::packing::PACKING_SPEC. For the nearest archived packing (mu_k), (hat u_varphi=(mu-mu_k)/|mu-mu_k|) and (P=J_mu^{mathsf T}hat u_varphi) is the Cartesian pullback through the analytic stacked Jacobian: the same increment-on-every-center lift as crate::soap::kick_packing_nu3. Then (gleftarrow g-2(gcdothat P)hat P) when (gcdot P>0). Leftover SOAP (p_i-mu) collapses across the paper funnels; this map does not.

When the packing mean is unavailable (fewer than a closed-shell neighbour cloud), the same Householder falls back to the COM-free Cartesian radius of each well. That fallback is not a packing walk.

The Householder on (g) is not conservative. The line search is run on the matching PES (E+V), a Gaussian hill on each known well in the same map. After a transformed quench that changes DECAF family, a raw-(E) polish sits on a true minimum of that well.

Variables

const LEAVE_BARRIER_FLOOR: f64

Barrier the first rung aims at, in units of the well depth per atom.

The ladder is walked in energy, not in length, because a barrier is what a Leave has to clear and the step that clears it depends on the curvature of the structure it starts from. The unit is the depth per atom of the minimum being left, which every cluster potential supplies and into which no morphology enters.

const LEAVE_BARRIER_GROWTH: f64

Ratio between the barriers successive rungs aim at.

const LEAVE_RIDGE_RUNGS: usize

Packing-map steps one Leave may accumulate along one direction.

The ceiling at LEAVE_WALK_CLIMB times the depth per atom is what ends the walk. This count lets that ceiling fire first at the first rung’s barrier.

const LEAVE_RUNGS: usize

Rungs a single Leave walks before it reports a refusal.

const LEAVE_RUNG_BACKTRACKS: usize

Bracketing and bisection steps the rung line search may take.

const LEAVE_RUNG_EXTRA: usize

Rungs walked past the first escape before the best of them is taken.

The first rung that leaves is not the best one to leave by. Measured from the LJ75 icosahedral minimum, the first escape lands anywhere between 5 and 25 (varepsilon) above the floor, and the spread across neighbouring rungs of one direction is most of that range.

const LEAVE_RUNG_RMSD: f64

Fallback Leave size when no curvature is available, in Cartesian RMSD.

A length is the wrong thing to fix. At a minimum the gradient vanishes, so a displacement of root-mean-square size (delta) over (N) atoms reaches (lambda Ndelta^2/2) under curvature (lambda), and the step that reaches a barrier (Delta) is rung_rmsd. Measured by Lanczos on the sealed LJ75 icosahedral minimum, (lambda_{min}=12.83), so the 7.48 (varepsilon) ico-Marks barrier of Wales and Doye is a step of (0.125) and the energy it actually costs is (7.11). This constant is nearly three times that step, and reach grows as (delta^2), so along the same mode it spends 58.9 (varepsilon): eight times the barrier, which is a melt and not a crossing. It survives only for callers that cannot afford a curvature pass.

Wales, D. J.; Doye, J. P. K. J. Phys. Chem. A 1997, 101, 5111. <https://doi.org/10.1021/jp970984n>

const LEAVE_WALK_CLIMB: f64

Largest climb in raw (E), in units of the well depth per atom.

The span rule alone accepts a step that puts two atoms on top of each other, because overlapping atoms move (mu) further than any rearrangement does. Wales and Doye put the LJ75 ico-Marks barriers at 8.69 and 7.48 (varepsilon) against a well depth of 5.28 (varepsilon) per atom, so four times the depth per atom clears the ridge the walk is meant to cross and refuses the ones it is not.

const LEAVE_WALK_DESCENTS: usize

Consecutive falls in raw (E) that name the far side of a ridge.

const LEAVE_WALK_SPAN: f64

Largest Cartesian RMSD one armed walk may cover in total.

A cooperative rearrangement moves atoms by about a nearest-neighbour separation. Past that the walk is not crossing a ridge between packings, it is pulling the cluster apart, and the ridge rule cannot end a climb that never turns over.

const LEAVE_WALK_STEP: f64

Largest Cartesian RMSD one armed step may take.

The hill width is the size of the well being left; the walk across it is not one step of that size. A trust radius equal to the width lets an accepted uphill step carry atoms through each other, and the walk then climbs for the whole iteration budget and reports a structure at (10^{11}varepsilon) that no quench recovers.

Functions

fn arm_leave(origin: ArrayView1<f64>, sigma_rmsd: f64, references: &[Vec<f64>])

Arm the transformed quench for one occupancy Leave.

origin is the live well this extra is leaving, and references are the packings already on file, as coordinates. Every one of them contributes a mode the Householder inverts, so the walk is pushed away from all of them and not only out of the well it starts in. The leftover-SOAP archive is a list of cloud means, not structures, so it cannot supply these: a Leave that reads it sees vectors of the wrong length and inverts nothing but its own well.

fn arm_leave_free(origin: ArrayView1<f64>, sigma_rmsd: f64, references: &[crate::catalog::PackingReference], temperature: f64, delta_t: f64)

Arm the Leave against the free-energy depth of the known packings.

The hill on well (k) is its depth along the packing coordinate, and a depth in potential energy is the wrong quantity to fill. What holds a chain is (F_k=E_k-TS_k), and on a cluster the configurational entropy of a packing is the log of how many ways the run reaches it, which is the arrival count the reference cloud now carries. On LJ75 that difference decides the answer: the Marks decahedron is 1.210082 (varepsilon) below the icosahedral floor, but the icosahedral shelf carries hundreds of minima against one narrow decahedral well, so at the run temperature (0.8) an entropy ratio of a thousand is already (5.5,varepsilon) in the other direction. A chain that stays icosahedral is at the free-energy minimum, and no deposit shaped like a potential energy can say otherwise (Mandelshtam, Frantsuzov, Calvo, J. Phys. Chem. A 2006, 110, 5326).

delta_t is the well-tempered scale of Barducci, Bussi and Parrinello (Phys. Rev. Lett. 2008, 100, 020603): a hill on a well that already carries (V_k) is scaled by (e^{-V_k/Delta T}), so the pile converges to (-Delta T/(T+Delta T)) times the free energy instead of growing without bound. Zero leaves every hill at full height, which is the non-tempered case of Laio and Parrinello (Proc. Natl. Acad. Sci. 2002, 99, 12562).

fn direction_curvature<G>(x: ArrayView1<f64>, direction: ArrayView1<f64>, epsilon: f64, mut grad: G) -> Option<f64>
where
    G: FnMut(ArrayView1<f64>) -> Option<Array1<f64>>

Curvature of the potential along a unit Cartesian direction.

(lambda_u = hat u^{mathsf T} H hat u) by central difference of the gradient, two evaluations. Reported for the record; the rung is sized by leave_packing_rung_to, which measures the energy instead of trusting this.

fn disarm()

Drop the transform. Later rgmin calls see the raw PES.

fn effective(x: ArrayView1<f64>, energy: f64, grad: Array1<f64>) -> (f64, Array1<f64>)

((E,g)) seen by rgmin: identity when unarmed, (E+V) and Householder-(g+nabla V) when a Leave is in flight.

fn is_armed() -> bool

Whether a Leave quench is currently transformed.

fn leave_packing_ladder<Q, E>(x: ArrayView1<f64>, cover_index: usize, references: &[Vec<f64>], species: Option<&[u32]>, mobile: Option<&[usize]>, depth_per_atom: f64, rungs: usize, quench: Q, mut energy: E) -> Option<(f64, Array1<f64>, usize)>
where
    Q: FnMut(ArrayView1<f64>) -> (f64, Array1<f64>),
    E: FnMut(ArrayView1<f64>) -> Option<f64>

Walk the Leave ladder until the quench installs a packing.

The ladder is walked in barrier, not in length. Rung k aims at rung_barrier and takes the step rung_rmsd that reaches it at the measured curvature of the structure being left, so the same ladder is the right size on a stiff cluster and on a soft one and needs no length constant. curvature is the softest non-rigid eigenvalue from a Lanczos pass (crate::curvature::curvature_features) and depth_per_atom is the well depth the ladder is scaled in; a caller with neither falls back to LEAVE_RUNG_RMSD.

quench is the caller’s relaxation, and the caller arms arm_leave around it so the walk is not pulled back into the wells on file. The rung that lands outside every packing on file is the Leave. None means the ladder is spent, which is a refusal to report, not a reason to fall back on a leftover-SOAP hole: a hole of the occupied packing quenches into the occupied packing.

fn leave_packing_ridge<R>(origin: ArrayView1<f64>, cover_index: usize, references: &[Vec<f64>], species: Option<&[u32]>, mobile: Option<&[usize]>, depth_per_atom: f64, relax_steps: usize, eval: R) -> Option<(f64, Array1<f64>, usize)>
where
    R: FnMut(ArrayView1<f64>, usize) -> (f64, Array1<f64>)

Follow one packing-map covering direction in barrier-sized steps.

Each step is leave_packing_rung_to from the current point, not from the well. The accept quench is raw (E): an invert-armed landing is not a packing. Quench starts only after the unquenched rise has reached rung_barrier at rung 2, which on LJ75 is past the Wales–Doye ico–Marks barrier scale. Measured covering ladders from the sealed icosahedral minimum leave to novel packings at (-371) to (-391); this walk is the same increment, accumulated far enough that a quench can sit on the far side.

fn leave_packing_rung(x: ArrayView1<f64>, cover_index: usize, rmsd: f64, references: &[Vec<f64>], species: Option<&[u32]>, mobile: Option<&[usize]>) -> Array1<f64>

One rung of the Leave ladder: a packing-map step of Cartesian size rmsd along the covering direction cover_index, pointed away from the packings on file.

fn leave_packing_rung_to<E>(x: ArrayView1<f64>, cover_index: usize, barrier: f64, references: &[Vec<f64>], species: Option<&[u32]>, mobile: Option<&[usize]>, energy: E) -> Array1<f64>
where
    E: FnMut(ArrayView1<f64>) -> Option<f64>

Rung sized so the measured energy rise is at most barrier.

A root-mean-square cap is the wrong bound for this direction. The packing pullback is not spread over the cluster: it concentrates on the few centres whose environment the increment changes, so an RMSD of 0.35 over 75 atoms can be one atom moving 1.7 sigma, straight through its neighbours. A crushed cluster has enormous negative transverse curvature, (-r^{-1},mathrm{d}V/mathrm{d}r) with the pair force deep in the repulsive wall, so a min-mode climb started there reports (lambda<0) and a flipped force on its first step and calls the ridge crossed. Measured on the sealed LJ75 icosahedral minimum: every one of twelve covering starts declared a crossing after one step at curvatures between (-7times10^5) and (-10^{13}), and the quench from there returned to the floor it started on.

rung_rmsd supplies the first trial length from the harmonic identity, and the length is then halved until the potential itself agrees that the step costs no more than the barrier. That bound cannot be met by a crushed structure, so the climb always starts on the landscape.

fn leave_packing_rung_to_dir<E>(x: ArrayView1<f64>, direction: &[f64], barrier: f64, species: Option<&[u32]>, mobile: Option<&[usize]>, mut energy: E) -> Array1<f64>
where
    E: FnMut(ArrayView1<f64>) -> Option<f64>

leave_packing_rung_to along a packed feature direction already in hand.

fn leave_packing_starts<Q>(x: ArrayView1<f64>, starts: &[Array1<f64>], references: &[Vec<f64>], mut quench: Q) -> Option<(f64, Array1<f64>, usize)>
where
    Q: FnMut(ArrayView1<f64>) -> (f64, Array1<f64>)

Walk a prepared ladder of starts until one quenches outside every packing on file, and keep the lowest-energy escape.

Split from leave_packing_ladder because building the rungs needs a gradient and walking them needs a quench, and a caller whose budget ledger backs both cannot hold the two closures at once.

fn leave_packing_toward<R>(origin: ArrayView1<f64>, target: ArrayView1<f64>, references: &[Vec<f64>], species: Option<&[u32]>, mobile: Option<&[usize]>, depth_per_atom: f64, relax_steps: usize, eval: R) -> Option<(f64, Array1<f64>, usize)>
where
    R: FnMut(ArrayView1<f64>, usize) -> (f64, Array1<f64>)

Walk the packing map from origin toward target’s packing mean.

Same energy-capped steps as leave_packing_ridge, but the direction is (widehat{mu_{mathrm{target}}-mu}) recomputed at the current point rather than a covering index. A covering of the high-dimensional packing sphere does not hit this vector.

fn lift() -> Option<(f64, f64)>

Hill amplitude (A) and width (sigma_varphi) of the armed transform, once the first transformed evaluation has set them.

The deposited potential is (Asum_k e^{-r_k^2/2sigma_varphi^2}), so what decides whether a wider ensemble can lift a chain out of its funnel is (A) against the barrier it has to clear: on LJ75 the icosahedral–Marks saddles sit 8.69 and 7.48 (varepsilon) above the shelf (Doye, Wales, Berry, J. Chem. Phys. 1995, 103, 4234). Reported so a run can be read against that number rather than against whether it happened to escape.

fn packing_cover_direction(x: ArrayView1<f64>, cover_index: usize, references: &[Vec<f64>]) -> Option<Vec<f64>>

Unit packing-map direction: covering point cover_index, orthogonal to (mu), signed away from the nearest packing on file.

Plasencia Gutiérrez, M.; Argáez, C.; Jónsson, H. J. Chem. Theory Comput. 2017, 13 (1), 125-134. <https://doi.org/10.1021/acs.jctc.5b01216>

fn packing_direction_between(from: ArrayView1<f64>, toward: ArrayView1<f64>) -> Option<Vec<f64>>

Unit packing-map direction from from toward toward.

(widehat{mu_{mathrm{to}}-mu_{mathrm{from}}}) with the breath along (mu_{mathrm{from}}) removed. This is the increment that walks one known packing onto another in the same map the covering Leave uses; a covering of (S^{d-1}) at (dsim 10^{2}) does not hit it by chance.

fn propose_leave_covers<E, R>(origin: ArrayView1<f64>, references: &[Vec<f64>], depth_per_atom: f64, species: Option<&[u32]>, n_try: usize, rng: &mut R, mut energy: E) -> Vec<(usize, Vec<f64>)>
where
    E: FnMut(ArrayView1<f64>) -> Option<f64>,
    R: rand::Rng + ?Sized

Unquenched Leave starts scored by packing histogram, for FunnelModel EI.

Each probe is one first-rung step or the fivefold residual. The quench is paid only for the cover EI (or Thompson) then selects.

fn rung_barrier(depth_per_atom: f64, rung: usize) -> f64

Barrier the rung at index rung aims at, from the well depth per atom.

fn rung_rmsd(curvature: f64, atoms: usize, barrier: f64) -> Option<f64>

Root-mean-square step whose harmonic reach is barrier.

(delta=sqrt{2Delta/(lambda N)}), which is Hop.rung_reaches_barrier in proofs/lean/Hop/LeavePacking.lean. None when the curvature, the count or the barrier is not positive, so a caller with no measurement falls back rather than inventing a length.

fn span(x: ArrayView1<f64>) -> f64

Distance from x to the nearest armed well.

Packing L2 (min_k|mu-mu_k|) when the wells carry a (nu=3) mean, otherwise COM-free RMSD. Zero when unarmed. Occupancy Leave keeps the invert walk when this rises.

fn step_rgmin<F>(opt: &mut crate::methods::warm_lbfgs::WarmLbfgs, x0: ArrayView1<f64>, max_iter: usize, mut fg: F) -> (f64, Array1<f64>)
where
    F: FnMut(ArrayView1<f64>) -> Option<(f64, Array1<f64>)>

rgmin step on the transformed surface: two-loop direction, accept a step that increases the span from the known wells. Span is packing L2 (min_k|mu-mu_k|) when the wells carry a (nu=3) mean, otherwise COM-free RMSD. Raw (E) may rise; that is the dimer walk away from an occupied packing.

The walk ends at the ridge, the way an activation-relaxation walk does: raw (E) climbs while the invert holds the descent component that points back at the wells on file, and the first steps where it falls again are the far side. Quenching from there is what decides which packing the walk landed in.

Two stopping rules that do not work, both measured from the LJ75 icosahedral minimum. Crossing the DECAF grain stops the walk on a distorted geometry that has not crossed anything: the quench that follows returns to (-396.282249), the floor it started on, every time. Having no rule at all lets the walk climb for the whole iteration budget and report structures between (10^4) and (10^{11}varepsilon), which is atoms sitting on each other. LEAVE_WALK_SPAN bounds the second case, since a raw energy that rises monotonically never turns over.

fn with_disarmed<T>(body: impl FnOnce) -> T) -> T

Run body on the raw surface, then restore the armed transform.

The polish that decides where a walk landed has to run on (E), not on (E+V): a hill the Leave put there is not part of the landscape, and a minimum of the sum is not a minimum of the potential.

fn with_hill_only<T>(body: impl FnOnce) -> T) -> T

Suppress only the Householder while body runs, keeping the hill.

The reflection is a poor-man’s min-mode: it inverts the force along one direction chosen in advance, which is a guess about where the exit is. A min-mode search does not need that guess, but it does need a force that is the gradient of something, and a reflected force is not.

The hill is a different object and is wanted: filling the occupied well by (A) lowers the saddle out of it by the same amount, so a dimer on (E+V) looks for a barrier the deposit has already paid down. That is the pairing metadynamics and saddle search are usually put in, and it is not available while the two are welded together.