Thermostats and Ensemble Control in Molecular Dynamics

Molecular dynamics integrates Hamilton's equations for $N$ interacting particles; because the flow conserves energy and preserves phase-space volume, a long trajectory explores a constant-energy hypersurface and its time averages approximate microcanonical ( $NVE$ ) ensemble averages. This is owned by the molecular-dynamics article. But most experiments are performed at fixed temperature and pressure, not fixed energy and volume. Ensemble control is the art of modifying the equations of motion so that the time-averaged sampling reproduces the canonical $NVT$ or isothermal–isobaric $NPT$ measure while disturbing the underlying physics as little as possible — and the central subtlety, which this article owns, is that getting the average temperature right does not guarantee getting the distribution right. Two philosophies coexist. Deterministic (extended-Lagrangian) thermostats — pioneered by Nosé and recast by Hoover — enlarge phase space with a small number of reservoir degrees of freedom whose equations of motion are engineered so that the marginal distribution of the physical variables is canonical. They are time-reversible and preserve smooth dynamical correlation functions, but they can fail to be ergodic. Stochastic thermostats — Andersen collisions, the Langevin equation, and global velocity rescaling — inject randomness that guarantees ergodicity and a strictly correct stationary distribution via detailed balance and the fluctuation–dissipation relation, at the price of replacing deterministic trajectories with a Markov process that damps long-time correlations. The Langevin thermostat's noise–friction lock $\sigma^2=2\gamma m k_BT$ is the second-kind fluctuation–dissipation relation owned by the langevin-equation and fluctuation-dissipation-theorem articles, recalled here by name. The barostat extends the same idea to the volume: treating $V$ (Andersen) or the full cell matrix (Parrinello–Rahman) as a piston with inertia $W$ , driven by the difference between the instantaneous internal pressure and the external set pressure, one samples constant-pressure ensembles and can observe diffusionless structural transformations. The theoretical demand throughout is that the modified flow possess the intended distribution as an invariant measure and, for trustworthy averages, be ergodic — a distinction dramatized by the notorious non-ergodicity of the bare Nosé–Hoover harmonic oscillator, which satisfies equipartition on average while sampling the wrong distribution. This article owns that owned result: correct average temperature does not imply correct ensemble.

Key takeaways

Plain Newtonian MD conserves energy and samples the microcanonical ( $NVE$ ) ensemble; to reach $NVT$ or $NPT$ one couples the system to fictitious or stochastic reservoirs that exchange energy (a thermostat) and volume (a barostat).
The Nosé–Hoover thermostat adds a dynamical friction variable $\xi$ with feedback $\dot\xi\propto(\sum p^2/m-g k_BT)$ ; the deterministic extended flow has a conserved quantity and, when ergodic, generates the canonical distribution in the physical subspace.
Stochastic thermostats — Andersen (velocity randomization), Langevin (friction + noise, locked by fluctuation–dissipation), and Bussi–Donadio–Parrinello velocity rescaling — enforce the canonical ensemble and are provably ergodic, at the cost of perturbing dynamical correlations.
Barostats (Andersen, Parrinello–Rahman, MTK) promote the volume (or cell matrix) to a dynamical variable driven by the internal–external pressure imbalance, yielding $NPT$ and enabling structural phase transitions.
A correct thermostat must reproduce the right distribution — equipartition $\langle\tfrac12 m v^2\rangle=\tfrac12 k_BT$ AND, ideally, ergodicity; the bare Nosé–Hoover scheme is famously non-ergodic for a single harmonic oscillator, cured by Nosé–Hoover chains.
Key idea: augment Hamilton's equations with reservoir variables (deterministic feedback) or with friction-plus-noise (stochastic) so that the canonical/NPT measure is a stationary — and ideally unique — solution of the phase-space continuity equation; the physical observables are then the marginals.

The core insight

Augment Hamilton's equations with reservoir variables (deterministic feedback) or friction-plus-noise (stochastic) chosen so that the canonical/ $NPT$ measure is a stationary — and ideally unique — solution of the phase-space continuity equation; the physical observables are recovered as marginals. Correct average temperature is necessary but not sufficient — the full distribution, hence ergodicity, is what matters.

Where this lives today

Constant-temperature production runs of proteins in water at 300 K use either a Nosé–Hoover chain (for accurate energy-conserving dynamics and transport properties) or the Bussi–Donadio–Parrinello velocity-rescaling thermostat (for robust, ergodic equilibration that barely disturbs diffusion constants); the once-popular Berendsen thermostat is now confined to initial equilibration because it samples no correct ensemble. An MTK or Parrinello–Rahman barostat lets the simulation cell change shape, enabling diffusionless solid–solid (martensitic) transformations and accurate equations of state. Because deterministic thermostats preserve smooth trajectories, they are preferred when extracting Green–Kubo transport coefficients (diffusion, viscosity, thermal conductivity), while stochastic thermostats with small friction are acceptable if the coupling stays weak versus the dynamical timescale. Replica-exchange (parallel tempering) MD runs many thermostatted copies across a temperature ladder and swaps with a Metropolis criterion — correct canonical sampling at each rung, guaranteed by a rigorous thermostat, is what makes the exchange obey detailed balance.

The physical picture

Imagine the simulated system bolted to a heat bath through a single rubber band whose stiffness is the thermostat. If the system runs hot (kinetic energy above $\tfrac{g}{2}k_BT$ ) the feedback variable grows and applies a drag that bleeds energy out; if it runs cold, the same variable turns negative and pumps energy in. The Nosé–Hoover $\xi$ is exactly this self-adjusting friction, with its own inertia $Q$ so the exchange is smooth and time-reversible rather than instantaneous. A stochastic thermostat is instead like randomly bumping particles: Andersen occasionally replaces a particle's velocity with a fresh Maxwell–Boltzmann draw, while Langevin continuously kicks each particle with white noise as friction removes the resulting drift. The single most important lesson — the owned thesis — is that equipartition is necessary but not sufficient. The bare Nosé–Hoover oscillator can satisfy $\langle p^2/m\rangle=k_BT$ while its trajectory is confined to a thin torus: the average temperature is right, but the distribution $f(q,p)$ is wrong. This is not a numerical artifact but a genuine failure of ergodicity — a single feedback variable cannot supply enough mixing for a stiff system. Adding a chain of thermostats (each thermostatting the previous) breaks the residual conservation laws and restores the full Gaussian. A kinetic-temperature histogram makes the failure visible: bare $NVE$ or non-ergodic Nosé–Hoover gives too-narrow fluctuations, while a proper canonical thermostat reproduces the broad canonical distribution. Detailed balance and fluctuation–dissipation are what guarantee stochastic thermostats sample correctly. For the Langevin thermostat the noise strength and friction are locked, $\sigma^2=2\gamma m k_BT$ — the second-kind fluctuation–dissipation relation owned by the langevin-equation article — so the bumps and the drag conspire to leave the canonical density unchanged, and the non-degenerate noise makes the process ergodic where deterministic schemes get stuck. The barostat applies the same logic to volume: the box becomes a piston pushed outward by the internal pressure and inward by $P_{\rm ext}$ , its response time (not the precise piston mass $W$ ) being the physical knob.

Figure

Kinetic-temperature histogram: NVE / non-ergodic NH too narrow, canonical thermostat is correct — the punchline

Figure

Two thermostat philosophies: Nosé–Hoover feedback friction ξ vs Langevin friction+noise

The core equations

Bare Newtonian MD conserves $H$ and phase-space volume, so it samples the microcanonical ensemble — the starting point that everything modifies.

\dot q_i=\frac{p_i}{m_i},\qquad \dot p_i=-\frac{\partial U}{\partial q_i} \qquad (1)

where: $q_i,p_i$ — position and momentum of degree of freedom $i$ ; $m_i$ — mass; $U(q)$ — potential. The Nosé–Hoover thermostat adds a feedback friction $\xi$ that drives the time-averaged kinetic energy to its equipartition value.

\dot q_i=\frac{p_i}{m_i},\quad \dot p_i=-\frac{\partial U}{\partial q_i}-\xi p_i,\quad \dot\xi=\frac{1}{Q}\!\left(\sum_i\frac{p_i^2}{m_i}-g\,k_BT\right) \qquad (2)

where: $\xi$ — thermostat friction variable; $Q$ — thermostat mass (reservoir inertia, units energy·time²); $g$ — number of physical degrees of freedom; $k_BT$ — target thermal energy. The extended system conserves a quantity $\mathcal H'$ — not a true Hamiltonian in these variables, but a stringent numerical check analogous to energy conservation in $NVE$ .

\mathcal H'=\sum_i\frac{p_i^2}{2m_i}+U(q)+\frac{Q\,\xi^2}{2}+g\,k_BT\,\eta,\qquad \dot\eta=\xi \qquad (3)

where: $\eta$ — thermostat coordinate; $\mathcal H'$ — conserved extended "energy." The Langevin thermostat replaces feedback with friction plus white noise; the two are locked by the fluctuation–dissipation relation so the canonical density is stationary.

\dot p_i=-\frac{\partial U}{\partial q_i}-\gamma\,p_i+\sqrt{2\gamma m_i k_BT}\;\eta_i(t),\quad \langle\eta_i(t)\eta_j(t')\rangle=\delta_{ij}\delta(t-t') \qquad (4)

where: $\gamma$ — friction coefficient (inverse coupling time); $\eta_i(t)$ — unit-variance Gaussian white noise; $\sqrt{2\gamma m_i k_BT}$ — noise amplitude fixed by fluctuation–dissipation (owned by the langevin-equation article). The Andersen thermostat instead thermalizes by stochastic collisions, each drawing a fresh Maxwellian velocity — trivially detailed-balanced against the canonical density.

p'_{i,\alpha}\sim\mathcal N(0,m_i k_BT)\ \text{at Poisson rate}\ \nu \qquad (5)

where: $\nu$ — stochastic collision frequency per particle; $\mathcal N(0,m_ik_BT)$ — Maxwell–Boltzmann velocity distribution at $T$ . The isotropic Andersen/MTK barostat promotes the volume to a piston driven by the pressure imbalance, giving (with a thermostat) the $NPT$ ensemble.

W\ddot\epsilon=3V(P_{\rm int}-P_{\rm ext}),\qquad \epsilon=\tfrac13\ln(V/V_0) \qquad (6)

where: $\epsilon$ — logarithmic strain (volume) variable; $W$ — barostat (piston) mass; $P_{\rm int}$ — instantaneous virial pressure; $P_{\rm ext}$ — external set pressure.

Derivation ACanonical invariance of the Nosé–Hoover flow

The goal here is to prove that the extended distribution $f\propto\exp[-\beta(H+Q\xi^2/2)]$ is a stationary solution of the (compressible) Liouville equation generated by the Nosé–Hoover equations, so that its physical $(q,p)$ marginal is exactly canonical. The subtle point is that the Nosé–Hoover flow is non-Hamiltonian — its phase-space divergence does not vanish — so we must use the continuity equation with the compressibility term, and the whole argument hinges on an exact cancellation between that term and the convective derivative of the trial density.

Step 1. Compute the phase-space compressibility. For a non-Hamiltonian flow $\dot{\mathbf x}=\mathbf v$ the density obeys $\partial_t f=-\mathbf v\!\cdot\!\nabla f-f(\nabla\!\cdot\!\mathbf v)$ . Only the friction term contributes to the divergence:

\nabla\!\cdot\!\mathbf v=\sum_i\partial_{p_i}(-\xi p_i)=-g\,\xi,

so the flow is compressible. Why $-g\xi$ : the friction $-\xi p_i$ acts on all $g$ momentum components, each contributing $-\xi$ .

Step 2. Convective derivative of the trial density $f\propto\exp[-\beta(H+Q\xi^2/2)]$ (independent of $\eta$ ):

\mathbf v\!\cdot\!\nabla\ln f=-\beta\Big[\sum_i(\dot q_i\partial_{q_i}H+\dot p_i\partial_{p_i}H)+Q\xi\dot\xi\Big].

Step 3. Evaluate the bracket. The Hamiltonian part gives $\sum_i\tfrac{p_i}{m_i}\partial_{q_i}U+(-\partial_{q_i}U-\xi p_i)\tfrac{p_i}{m_i}=-\xi\sum_i\tfrac{p_i^2}{m_i}$ , and $Q\xi\dot\xi=\xi(\sum_i p_i^2/m_i-g k_BT)$ . Adding them, the $\xi\sum p^2/m$ terms cancel, leaving $-g k_BT\,\xi$ . Hence $\mathbf v\!\cdot\!\nabla\ln f=-\beta(-g k_BT\,\xi)=g\xi$ , so $\mathbf v\!\cdot\!\nabla f=g\xi f$ .

Step 4. Combine with the compressibility term:

\partial_t f=-\mathbf v\!\cdot\!\nabla f-f(\nabla\!\cdot\!\mathbf v)=-g\xi f-(-g\xi f)=0.

Why this is the whole result: the trial density is stationary — the cancellation $g\xi f-g\xi f=0$ is exact. Marginalizing over the Gaussian $\xi$ and the flat $\eta$ returns $f(q,p)\propto e^{-\beta H}$ , the canonical ensemble, provided the flow is ergodic. Verified: the Nosé–Hoover conserved quantity $\mathcal H'$ is constant along a numerically integrated SHO+thermostat trajectory, confirming the flow, in verify_aq-sm-thermostats-md-ensembles.py.

Go deeper · Derivation B · 1 minLangevin thermostat: canonical stationary density fixes the noise amplitude

Derivation BLangevin thermostat: canonical stationary density fixes the noise amplitude

This derivation shows that the Langevin thermostat samples the canonical ensemble only when the noise strength and friction obey the fluctuation–dissipation relation $\sigma^2=2\gamma m k_BT$ — and that this same non-degenerate noise makes the process ergodic, curing the Nosé–Hoover pathology. We insert the canonical density into the Klein–Kramers (Fokker–Planck) equation of the Langevin dynamics and demand stationarity; the requirement pins the noise amplitude uniquely. The general fluctuation–dissipation content is owned by the langevin-equation and fluctuation-dissipation-theorem articles; here we use the result to fix the thermostat.

Step 1. Write the Klein–Kramers equation for one degree of freedom, $\dot q=p/m$ , $\dot p=-\partial_q U-\gamma p+\sigma\eta$ :

\partial_t f=-\frac{p}{m}\partial_q f+\partial_p\!\big[(\partial_q U+\gamma p)f\big]+\frac{\sigma^2}{2}\partial_p^2 f.

Step 2. Insert the canonical ansatz $f_{\rm st}\propto\exp[-\beta(p^2/2m+U)]$ . The streaming (Hamiltonian) terms $-\tfrac{p}{m}\partial_q f+\partial_p(\partial_q U\,f)$ cancel identically because $f_{\rm st}$ depends on $q,p$ only through $H$ . Why: a function of $H$ is annihilated by the Liouville streaming operator.

Step 3. Balance friction against diffusion. Using $\partial_p f=-\beta\tfrac{p}{m}f$ and $\partial_p^2 f=(\beta^2 p^2/m^2-\beta/m)f$ , the remaining terms are

\gamma f+\gamma p\,\partial_p f+\tfrac{\sigma^2}{2}\partial_p^2 f=f\Big[\gamma-\gamma\tfrac{\beta p^2}{m}+\tfrac{\sigma^2}{2}\big(\tfrac{\beta^2 p^2}{m^2}-\tfrac{\beta}{m}\big)\Big].

Step 4. Require this to vanish for all $p$ . The $p^2$ coefficient gives $-\tfrac{\gamma\beta}{m}+\tfrac{\sigma^2\beta^2}{2m^2}=0\Rightarrow\boxed{\sigma^2=2\gamma m k_BT}$ , and the constant term vanishes for the same $\sigma^2$ . Why this is the fluctuation–dissipation relation: the noise variance is locked to the friction and temperature; mismatch them and the stationary temperature differs from the target. Because the noise is non-degenerate, the process is ergodic — it samples the full canonical density regardless of system stiffness, unlike bare Nosé–Hoover. Verified: the Langevin thermostat on a single harmonic degree of freedom samples $\langle\tfrac12 m v^2\rangle=\tfrac12 k_BT$ with $\sigma^2=2\gamma m k_BT$ in verify_aq-sm-thermostats-md-ensembles.py.

Limiting cases

Case	Result	Meaning
Coupling $\to0$	$Q\to\infty,\ \gamma\to0$	The reservoir decouples; dynamics revert to microcanonical $NVE$ . Energy drifts only from integration error.
Strong Langevin friction	$\gamma\gg\omega_{\rm char}$	Overdamped (Smoluchowski) limit: inertia irrelevant, dynamics diffusive, vibrational spectra lost — but sampling stays canonical.
Andersen collision rate	$\nu\to0$ or $\nu\to\infty$	Too small → poor thermalization (near- $NVE$ ); too large → velocity correlations destroyed each step (unphysical). Optimal $\nu\sim$ one collision per correlation time.
Bare Nosé–Hoover on SHO	$H=\tfrac{p^2}{2m}+\tfrac12 m\omega^2 q^2$	Famously NON-ergodic — phase space splits into invariant tori. A Nosé–Hoover chain or stochastic thermostat is required.
Barostat mass extremes	$W\to\infty\ (NVT)$ , $W\to0$ (unstable)	Large piston mass freezes the volume (recovers $NVT$ ); vanishing mass gives runaway oscillations.

Worked examples

Example 1 — Langevin thermostat samples the canonical variance. Integrate a single harmonic degree of freedom ( $m=1$ , $\omega=1$ , $k_BT=1$ ) under the Langevin thermostat with $\gamma=1$ and noise $\sigma=\sqrt{2\gamma m k_BT}$ , and verify the long-run kinetic energy.

Step 1. Noise amplitude. The fluctuation–dissipation lock gives $\sigma^2=2\gamma m k_BT=2(1)(1)(1)=2$ , so $\sigma=\sqrt2=1.4142$ .

Step 2. Canonical target. Equipartition predicts $\langle\tfrac12 m v^2\rangle=\tfrac12 k_BT=\mathbf{0.5}$ per degree of freedom, i.e. $\langle v^2\rangle=1$ .

Step 3. Simulation. Integrating the Langevin SDE for $\sim10^6$ steps (Euler–Maruyama, $\Delta t=0.01$ ) and time-averaging gives $\langle\tfrac12 v^2\rangle\to0.5$ within Monte-Carlo error, and the velocity histogram matches the Maxwell–Boltzmann Gaussian of variance $k_BT/m=1$ . The thermostat reproduces the exact canonical kinetic energy $\tfrac12 k_BT$ per degree of freedom — the correct distribution, not merely the correct mean.

Example 2 — Nosé–Hoover conserved quantity stays constant. For the same SHO with a Nosé–Hoover thermostat ( $Q=1$ , $g=1$ , $k_BT=1$ ), check that the extended "energy" $\mathcal H'$ is conserved along the trajectory.

Step 1. Initial value. With $q_0=1$ , $p_0=0$ , $\xi_0=0$ , $\eta_0=0$ : $\mathcal H'_0=\tfrac12 p_0^2+\tfrac12\omega^2 q_0^2+\tfrac12 Q\xi_0^2+g k_BT\,\eta_0=0+0.5+0+0=\mathbf{0.5}$ .

Step 2. Time evolution. Integrating (2) with a time-reversible scheme, $\mathcal H'$ should be conserved to integrator accuracy. After $10^5$ steps the drift $|\mathcal H'(t)-\mathcal H'_0|$ stays below $\sim10^{-3}$ , oscillating rather than growing — the same bounded behaviour as energy in a symplectic $NVE$ run.

Step 3. Interpretation. A conserved $\mathcal H'$ is the thermostat's analogue of energy conservation and the primary numerical health check; but — the owned caveat — even a perfectly conserved $\mathcal H'$ does not guarantee the physical marginal is canonical, because the bare scheme may be non-ergodic.

Origins · 1 minmolecular dynamics and how this was discovered

molecular dynamics

method, thermostats developed 1980–1994 · The extension of molecular dynamics to controlled ensembles was built by a sequence of key ideas: Andersen's stochastic-collision thermostat and constant-pressure MD (1980); Nosé's extended-Lagrangian formulation of canonical dynamics (1984), recast by Hoover into the friction-variable form now universally used (1985); the Parrinello–Rahman variable-cell barostat (1981); the Martyna–Klein–Tuckerman Nosé–Hoover chains that repaired the ergodicity failure (1992); and the Bussi–Donadio–Parrinello stochastic velocity rescaling (2007). Together they turned MD from a microcanonical curiosity into a rigorous sampler of the ensembles experiments actually probe.

History

1980

— Andersen introduces the stochastic-collision thermostat and constant-pressure MD

1981

— Parrinello and Rahman formulate the variable-cell (shape-changing) barostat

1984

— Nosé derives canonical MD from an extended Lagrangian

1985

— Hoover recasts Nosé's scheme into the friction-variable equations of motion

1992

— Martyna, Klein and Tuckerman introduce Nosé–Hoover chains, restoring ergodicity

2007

— Bussi, Donadio and Parrinello introduce stochastic velocity rescaling

Going further

References. Frenkel & Smit, Understanding Molecular Simulation, 2nd ed., Ch. 6 · Tuckerman, Statistical Mechanics: Theory and Molecular Simulation, Ch. 4–5 · Nosé, J. Chem. Phys. 81, 511 (1984) · Hoover, Phys. Rev. A 31, 1695 (1985) · Andersen, J. Chem. Phys. 72, 2384 (1980) · Martyna, Klein & Tuckerman, J. Chem. Phys. 97, 2635 (1992) · Bussi, Donadio & Parrinello, J. Chem. Phys. 126, 014101 (2007) · Parrinello & Rahman, J. Appl. Phys. 52, 7182 (1981).

Provenance. A verifier should: (1) confirm the Nosé–Hoover stationarity cancellation $g\xi f-g\xi f=0$ and that $\mathcal H'$ is conserved along a numerically integrated trajectory; (2) reproduce the canonical solution of the Klein–Kramers equation fixing $\sigma^2=2\gamma m k_BT$ , and show the Langevin thermostat samples $\langle\tfrac12 m v^2\rangle=\tfrac12 k_BT$ ; (3) note the non-ergodicity of the bare Nosé–Hoover harmonic oscillator. All expressions are analytic and textbook-standard. The fluctuation–dissipation relation and Langevin physics are owned by the langevin-equation and fluctuation-dissipation-theorem articles; MD integration by the molecular-dynamics article; the canonical ensemble by SM2.

Builds on: Molecular Dynamics Simulation · Langevin Equation and Stochastic Differential Equations · Fluctuation–Dissipation Theorem · Fokker–Planck Equation · Detailed Balance and Microscopic Reversibility. Leads to: Monte Carlo Method · Kinetic Theory and the Boltzmann Equation · Metropolis–Hastings Algorithm.

Why isn't velocity rescaling a valid thermostat?

Simple velocity rescaling (multiply all velocities by $\sqrt{T_{\rm target}/T_{\rm inst}}$ each step) forces the kinetic energy to a fixed value, which is exactly wrong: the canonical ensemble requires the kinetic energy to fluctuate with variance $\tfrac{N_f}{2}(k_BT)^2$ . Suppressing those fluctuations samples no well-defined ensemble and drains energy into center-of-mass motion (the "flying ice cube"). The Bussi–Donadio–Parrinello scheme fixes this by rescaling stochastically toward a target drawn from the correct canonical kinetic-energy distribution, restoring the right fluctuations while keeping the method simple.

What exactly goes wrong with the bare Nosé–Hoover oscillator?

For a single harmonic oscillator the Nosé–Hoover dynamics live in the 3D space $(q,p,\xi)$ (with $\eta$ decoupled). Numerically, for typical $Q$ a large fraction of initial conditions lie on smooth invariant 2-tori: the trajectory never reaches the high- $|\xi|$ regions that correspond to the tails of the canonical density. Equipartition $\langle p^2\rangle=k_BT$ can still hold by symmetry, but the distribution $f(q,p)$ is wrong. The single feedback variable simply cannot supply enough mixing — the effective dynamics are quasi-integrable. This is the textbook demonstration that a correct average is not a correct ensemble.

How does a Nosé–Hoover chain fix the non-ergodicity?

Instead of one thermostat variable, you couple $\xi_1$ to a second thermostat $\xi_2$ , and that to $\xi_3$ , and so on: $\dot\xi_1=\tfrac1{Q_1}(p^2-k_BT)-\xi_1\xi_2$ , $\dot\xi_2=\tfrac1{Q_2}(Q_1\xi_1^2-k_BT)-\xi_2\xi_3,\dots$ The chain destroys the spurious conservation laws that confine the bare scheme to tori, restoring ergodicity; a chain of length $3$ – $5$ suffices in practice (Martyna–Klein–Tuckerman 1992). One verifies, exactly as for the single thermostat, that $f\propto\exp[-\beta(H+\sum_k Q_k\xi_k^2/2)]$ is invariant, so the physical marginal remains canonical.

Deterministic or stochastic thermostat — how do I choose?

It depends on what you measure. Deterministic thermostats (Nosé–Hoover chains) preserve smooth, time-reversible trajectories, so they are preferred for dynamical quantities — Green–Kubo transport coefficients, spectra, correlation functions — where you do not want noise contaminating the dynamics. Stochastic thermostats (Langevin, Andersen, BDP) are provably ergodic and give strictly correct static sampling, ideal for equilibration and free-energy estimates, but their noise damps long-time correlations and must be kept weak (small $\gamma$ ) if you also want dynamics. A common workflow: equilibrate stochastically, then run production with a Nosé–Hoover chain.

What does the barostat actually do to the equations of motion?

It promotes the volume (or the full cell matrix) to a dynamical variable with its own inertia $W$ , driven by the difference between the instantaneous internal (virial) pressure and the external set pressure: $W\ddot\epsilon=3V(P_{\rm int}-P_{\rm ext})$ . The particle coordinates are simultaneously scaled with the box so that fractional positions are preserved. When $P_{\rm int}>P_{\rm ext}$ the box expands; equilibrium is reached at $\langle P_{\rm int}\rangle=P_{\rm ext}$ . Combined with a thermostat this yields the $NPT$ ensemble. The Parrinello–Rahman variant lets the box change shape as well as size, which is essential for capturing anisotropic stress and solid–solid structural transitions.

Why does the Andersen thermostat trivially satisfy detailed balance?

Because each Andersen "collision" replaces a particle's velocity with a fresh draw from the equilibrium Maxwell–Boltzmann distribution $\mathcal N(0,m k_BT)$ — the very distribution the canonical ensemble assigns to that velocity. Proposing a move directly from the target distribution means the acceptance is effectively unity and the canonical density is left invariant. Between collisions the dynamics are ordinary energy-conserving Newtonian motion, and the collisions are Poisson events at rate $\nu$ . The trade-off is that each collision erases that particle's velocity memory, so too-frequent collisions destroy the physical dynamics — the collision rate must be tuned to about one per correlation time.