mod gle

module gle

Generalized Langevin equation colored-noise thermostat (optimal sampling). Generalized Langevin equation (GLE) colored-noise thermostat, in the Markovian extended-phase-space form of Ceriotti, Bussi and Parrinello.

White-noise Langevin dynamics has a single friction gamma, so it can critically damp only one frequency; across a spectrum of curvatures (an ill-conditioned objective) most modes are far from critical and decorrelate slowly. The GLE replaces the white noise by a colored noise generated from ns auxiliary momenta, whose drift matrix A shapes a frequency-dependent friction. A well-chosen A makes the sampling efficiency nearly flat across a target frequency band – “optimal sampling” – so conditioning is handled by the noise spectrum exactly as the 1/sqrt(D) scale handles dimension.

The per-step propagator follows i-PI’s ThermoGLE: T = exp(-dt A), S S^T = k_B (C - T C T^T), and s <- T s + S xi`, with xi` standard normal and ``s[0] the physical (mass-scaled) momentum. For canonical sampling C = temperature * I, which makes the auxiliary process leave the Boltzmann momentum distribution invariant for ANY valid drift A – the “a la carte” decoupling of correctness (set by C) from acceleration (set by A). matrix_exp is a Taylor scaling-and-squaring and the noise factor is an LDL^T square root, so no linear-algebra backend is needed.

Variables

const OPTIMAL_SAMPLING_NS: usize

Number of auxiliary momenta in the fitted optimal-sampling drift.

Functions

fn ldl_sqrt(m: &Array2<f64>) -> Array2<f64>

Stabilised LDL^T square root: returns lower-triangular S with S S^T = M for a symmetric positive-semidefinite M, zeroing any negative pivot. This is i-PI’s stab_cholesky; it needs no eigensolver and is numerically robust for the GLE noise covariance, which can be singular at small dt.

fn matrix_exp(m: &Array2<f64>, n_taylor: usize, n_square: usize) -> Array2<f64>

Matrix exponential by Taylor series with scaling and squaring.

exp(M) = (exp(M / 2^k))^{2^k} with the inner exponential a truncated Taylor series evaluated by Horner. Matches i-PI’s matrix_exp; for the small GLE drift matrices (dimension ns + 1, a handful of rows) it is exact to machine precision with n_taylor = 20, n_square = 10.

fn optimal_sampling_drift(omega0: f64) -> Array2<f64>

Build the fitted optimal-sampling drift, scaled to characteristic frequency omega0 (the band becomes [omega0, 100 omega0]).

The drift has units of frequency, so the whole fitted reference scales linearly: A(omega0) = omega0 * A_ref. Each bath is a damped-oscillator block with eigenvalues -Gamma +/- i Omega skew-coupled to the physical momentum; the oscillator and coupling terms are skew, so A + A^T is diagonal and non-negative and the fluctuation-dissipation relation A C + C A^T = B B^T holds for C = I (canonical sampling). Spreading the baths across the band flattens the friction – hence the sampling efficiency – the colored-noise analogue of critical damping over a whole spectrum rather than at one frequency.

Structs and Unions

struct GleThermostat

A GLE thermostat with precomputed drift and noise factors for one timestep.

ns: usize

Number of auxiliary momenta (matrix dimension is ns + 1).

Implementations

impl GleThermostat

Functions

fn canonical(a: &Array2<f64>, dt: f64, temperature: f64, kb: f64) -> Self

Canonical thermostat: C = temperature * I, leaving the Boltzmann momentum law invariant for any valid drift a.

fn new(a: &Array2<f64>, c: &Array2<f64>, dt: f64, kb: f64) -> Self

Build from an explicit drift a and covariance c at timestep dt (Boltzmann constant kb, in reduced units kb = 1). a and c are (ns+1) x (ns+1).

fn sample_stationary<R: Rng + ?Sized>(&self, c: &Array2<f64>, dim: usize, kb: f64, rng: &mut R) -> Array2<f64>

Draw an auxiliary state from its stationary covariance C (mass-scaled physical momentum in row 0), used to initialise a chain at equilibrium.

fn step<R: Rng + ?Sized>(&self, s: &mut ArrayViewMut2<f64>, rng: &mut R)

Advance the auxiliary state one timestep in place. s is (ns+1) x dim; row 0 holds the physical (mass-scaled) momentum of each coordinate, rows 1..=ns the auxiliary momenta. Applies s <- T s + S xi.