mod dos

module dos

A posterior over the density of minima, and acceptance by entropy rather than by energy. A posterior over the density of minima, and acceptance by entropy.

Measured negative as an acceptance rule. Flat-histogram acceptance built on this estimator was refuted against paired controls (numbers and mechanism in docs/derivations/), and no acceptance rule here is part of the recommended configuration. The joint multi-sweep estimator itself stands, with its tests, as a density-of-states tool.

Why the acceptance rule is the trap

Basin hopping samples minima from pi(m) ∝ exp(-E~_m / T) on the quenched surface. Marginalised onto energy that is

p(E~) ∝ g(E~) exp(-E~ / T)

where g counts minima per unit quenched energy. On a multi-funnel landscape g is not flat and not slowly varying: the icosahedral funnel of a 38-point Lennard-Jones cluster holds exponentially more minima than the face-centred-cubic funnel that contains the global minimum. At a temperature loose enough to keep the chain moving, g dominates and the chain sits where the states are. At a temperature tight enough for energy to beat multiplicity, the chain stops moving. The obstruction is a property of the target distribution and no proposal fixes it, which is why better moves – Hamiltonian dynamics, tempering, surrogate-accelerated acceptance, escape along soft modes – all measure at chance on this landscape.

Sampling flat in energy instead

Give each state weight 1 / g(E~) rather than exp(-E~ / T). The marginal over energy becomes g(E~) / g(E~), a constant, so the chain spends as much time at each quenched energy as at any other, and the rare deep energies get the same share as the abundant shallow ones. With S = ln g the acceptance for a symmetric proposal is

a = min(1, exp(-[S(E~_new) - S(E~_old)]))

Metropolis with entropy where the energy was. This is the multicanonical construction of Berg and Neuhaus (doi:10.1016/0370-2693(91)90256-P) and the flat-histogram target of Wang and Landau (doi:10.1103/PhysRevLett.86.2050), applied to the density of minima rather than of configurations, which is the object Bogdan, Wales and Calvo estimate for cluster thermodynamics (doi:10.1063/1.2148958).

It is also the coarse-graining the landscape actually admits. Under this target two minima at the same quenched energy are the same state, so the state space collapses from an exponential number of minima onto a single energy axis, and no structural descriptor is needed to say which states are equivalent.

What is Bayesian about it

S is unknown and has to be learned from the run. Wang and Landau learn it by a histogram with a multiplicative schedule f -> sqrt(f), whose error saturates rather than converging (Belardinelli and Pereyra, doi:10.1103/PhysRevE.75.046701). Here S carries a posterior instead:

• Likelihood. Under a frozen weight w = exp(-S^) the chain’s stationary distribution over bins is p_k ∝ exp(S_k - S^_k), so a visit histogram n_k observes the residual directly: n_k ~ Poisson(exp(S_k - S^_k + c)).

• Prior. A second-difference Gaussian random walk on S, the intrinsic smoothing prior of Poisson density estimation. It is proper on second differences, improper on the two degrees of freedom it deliberately does not constrain – level and slope – which is what lets the posterior extrapolate S linearly past the deepest energy yet seen, with variance that grows as it goes.

• Exploration. The sampling weight is a draw from the posterior rather than its mean, so a bin the run has little evidence about is over-weighted in proportion to how little is known, and gets visited. Thompson sampling on the weight function.

The unvisited bin is where this separates from Wang-Landau in kind rather than in degree: a histogram method has nothing to say about a bin with no counts, and the prior gives it a value with an honest error bar.

Validity

Changing the weight function invalidates detailed balance, so the weight is frozen for a fixed sweep and refreshed between sweeps. Each sweep is then an exact Markov chain for its own target and the sequence is the standard stochastic-approximation scheme, not an adaptive chain whose invariance has to be argued.

Variables

const BINS: usize

Default number of energy bins.

Functions

fn cut_from(seen: &[f64], quantile: f64) -> Option<(f64, f64)>

The cut and width implied by a sweep’s visited energies.

The cut is a quantile of what the chain saw and the width is a spread of the same sample, so both descend as the run descends and neither is a tuned constant. A run that has stopped improving keeps the cut where it is, which is the right behaviour: the schedule follows progress rather than the clock.

Structs and Unions

struct CutWeight

A flat target below a cut and a Boltzmann one above it.

Flat sampling across a whole energy range is the wrong target for a search. It buys barrier crossing by giving every energy an equal share of the run, and most of the range holds nothing worth the share: measured on 38 points, flat over the visited range solved 3 seeds in 24 where the plain Metropolis rule solved 15, because the budget went to the high-energy sea.

Restricting the flat region fixes that, and the restriction cannot be a wall. Crossing between funnels requires rising in quenched energy, so forbidding the rise rebuilds the trap the flat target was adopted to dissolve. What works is flat below and suppressed above:

-ln pi(E~) = S(E~) + max(0, E~ - c) / w

Below c the multiplicity of a funnel does not decide anything, so a barrier inside the reachable region is invisible. Above c the cost grows linearly, so the chain can still climb w-worth to cross and does not walk off to the top. The cut descends with the chain, which is the annealing: the schedule is read off the run’s own visited energies rather than supplied as a cooling curve.

weight: Weight

The entropy curve.

cut: f64

Energy above which the target stops being flat.

width: f64

Width of the suppression above the cut, in energy.

Implementations

impl CutWeight

Functions

fn accept_prob(&self, e_old: f64, e_new: f64, bias_delta: f64) -> f64

Acceptance probability for a move between two quenched energies, with an additive term for any external bias the caller carries.

fn cost(&self, e: f64) -> f64

The negative log target at an energy.

struct DensityOfStates

A posterior over S = ln g on a binned energy axis.

refreshes: usize

Refreshes performed.

clamped: usize

Samples that fell outside the window and were clamped into an end bin.

Implementations

impl DensityOfStates

Functions

fn bins(&self) -> usize

Number of bins.

fn centre(&self, k: usize) -> f64

Centre energy of bin k.

fn draw<R: Rng + ?Sized>(&mut self, rng: &mut R) -> Weight

A weight function drawn from the posterior, frozen for one sweep.

Thompson sampling on S: one draw per bin, correlated through the same standard normal along the axis so the drawn curve is smooth rather than a bin-independent rattle, since an unsmooth weight would make the acceptance ratio between neighbouring bins meaningless.

fn entropy(&self, e: f64) -> (f64, f64)

Posterior mean and standard deviation of S at an energy.

Inside the window this reads the bin. Below the window it extrapolates linearly from the slope at the low end, with the standard deviation growing linearly in the distance extrapolated, which is the behaviour the second-difference prior implies: level and slope are unconstrained by it, so a step past the evidence keeps the slope and loses confidence at a constant rate.

fn inside(&self, e: f64) -> bool

Whether an energy sits inside the binned window.

fn mean_weight(&mut self) -> Weight

The posterior mean as a weight function, for runs that want no exploration term.

fn new(lo: f64, hi: f64, bins: usize) -> Self

A posterior over [lo, hi] split into bins bins, flat and uncertain before any evidence arrives.

The initial state is deliberately S = 0 everywhere with a wide spread: flat S makes the first sweep plain random-walk acceptance, which is the least committed thing to do before the run has seen an energy.

fn observe(&mut self, e: f64)

Records a visit at quenched energy e.

fn pending(&self) -> usize

Samples pending since the last refresh.

fn refresh(&mut self) -> bool

Folds the pending histogram into the posterior and clears it.

Returns false when there is nothing to learn from.

fn temperature(&self, e: f64) -> (f64, f64)

Posterior mean and standard deviation of the statistical temperature at an energy.

The statistical temperature is the reciprocal slope of the entropy, T_S = (dS/dE~)^-1, and it is the temperature at which a chain standing at that energy is critically mobile: the Metropolis ratio for a typical move is order one, neither frozen nor free. Kim, Straub and Keyes use it to drive dynamics directly (doi:10.1103/PhysRevLett.97.050601).

Its behaviour is what a search wants and is the opposite of a cooling curve. At a funnel floor few new minima appear per unit energy, so the slope is small, the temperature is high, and the chain climbs out. In the high-energy sea minima are dense, the slope is large, the temperature is low, and the chain descends. It is read off the run’s own entropy rather than scheduled against the clock, so a chain that has stopped finding new minima at its current depth heats itself.

The slope is taken by central difference on the fitted mean, which is smooth by construction under the second-difference prior, and its spread follows from the two bins it is taken across.

fn top_seen(&self) -> usize

Highest bin index the run has ever visited.

struct EnergyBias

A well-tempered bias deposited in quenched energy.

The per-basin bias this crate carries fills the basin the chain stands in. Measured at 38 points, that is not the shape of the trap: of 17 failing seeds, 12 end at exactly -173.252378 and 4 at -173.134317, the floor of the icosahedral funnel. A funnel holds exponentially many basins, so filling them one at a time cannot fill it, and coarsening the basin metric does not help because a single length in coordinate space cannot tell a funnel’s variants from a different funnel: at a radius of 0.7 the run registers 365 basins and solves 55 of 72, at 2.0 it registers 7 and solves 24.

Energy separates what the length cannot. The funnel floor is an energy the chain returns to, so depositing there fills the funnel rather than one of its basins, and the coordinate costs nothing because every hop computes it.

Deposits are well tempered (Barducci, Bussi and Parrinello, doi:10.1103/PhysRevLett.100.020603): height falls as exp(-V/((gamma-1) T)) where the bias already stands, so the sum converges instead of growing without bound, and the sampled distribution is the well-tempered ensemble of Bonomi and Parrinello (doi:10.1103/PhysRevLett.104.190601) rather than a flat one. This is the part the flat-histogram acceptance got wrong: it forced every energy to an equal share, where tempering only broadens what the chain already samples.

w0: f64

Initial deposit height.

gamma: f64

Tempering factor. One is no bias; large is untempered metadynamics.

sigma_bins: f64

Width of a deposit, in bins.

deposits: usize

Deposits made.

Implementations

impl EnergyBias

Variables

const FILL_DEPOSITS: f64

The number of deposits that fill a well one standard deviation deep.

The one free number in this construction, and it is dimensionless: every other scale is taken from the run’s own quenched-energy distribution, so nothing here carries units of a particular system’s energy and nothing is set per system.

Functions

fn at(&self, e: f64) -> f64

The bias at an energy, held flat outside the binned window so a chain that leaves the range is neither pushed back nor pulled out.

fn delta(&self, e_old: f64, e_new: f64, temp: f64) -> f64

The bias difference a move carries, in units of temperature, ready to be added to a Metropolis exponent.

fn deposit(&mut self, e: f64, temp: f64)

Deposits at an energy, at the well-tempered height for the bias already standing there.

fn from_sample(seen: &[f64], temp: f64, bins: usize) -> Option<Self>

A bias whose scales come from a sample of quenched energies.

The tempering factor is set by (gamma - 1) T = sigma, so the bias broadens the energy distribution by about its own width, and the deposit height by w0 = sigma / FILL_DEPOSITS, so a well one standard deviation deep fills in a fixed number of deposits whatever the system’s energy scale. The window spans the sample padded by a standard deviation on each side.

Returns None when the sample is too small or degenerate to set a scale, which is the honest answer rather than a default.

fn new(lo: f64, hi: f64, bins: usize, w0: f64, gamma: f64) -> Self

A bias over [lo, hi], empty until something is deposited.

fn peak(&self) -> f64

Largest bias standing anywhere, for reporting how filled the range is.

struct Weight

A frozen weight function: the entropy curve one sweep accepts against.

Implementations

impl Weight

Functions

fn accept_prob(&self, e_old: f64, e_new: f64) -> f64

Acceptance probability for a move between two quenched energies.

Metropolis against 1 / g: a move to a rarer energy is always taken, and a move to a more abundant one is taken with the ratio of multiplicities. Both directions in energy are treated alike, which is what makes the sampled histogram flat rather than funnel-weighted.

fn at(&self, e: f64) -> f64

The entropy this weight assigns to an energy, extrapolated linearly outside the binned window.