Statistical Mechanical Ensembles Research Articles

The imprint of conservation laws on correlated particle production

The study of event-by-event fluctuations of net-baryon number in a subspace of full phase space is a promising direction for deciphering the structure of strongly interacting matter created in collisions of relativistic heavy nuclei. Such fluctuations are generally suppressed by exact baryon number conservation. Moreover, the suppression is stronger if baryon number is conserved locally. In this report we present a conceptually new approach to quantify correlations in rapidity space between baryon-antibaryon, baryon-baryon, and antibaryon-antibaryon pairs and demonstrate their impact on net-baryon number fluctuations. For the special case of Gaussian rapidity distributions, we use the Cholesky factorization of the covariance matrix, while the general case is introduced by exploiting the well-known Metropolis and Simulated Annealing methods. The approach is based on the use of the canonical ensemble of statistical mechanics for baryon number and can be applied to study correlations between baryons as well as strange and/or charm hadrons. It can also be applied to describe relativistic nuclear collisions leading to the production of multi-particle final states. One application of our method is the search for formation of proton clusters at low collision energies emerging as a harbinger of the anticipated first-order chiral phase transition. In a first step, the results obtained are compared to the recent measurements from the CERN ALICE collaboration. Such investigations are key to explore the phase diagram of strongly interacting matter and baryon production mechanisms at energy scales from several GeV to several TeV.

Methods for Classical-Mechanical Molecular Simulation in Chemistry: Achievements, Limitations, Perspectives.

More than a half century ago it became feasible to simulate, using classical-mechanical equations of motion, the dynamics of molecular systems on a computer. Since then classical-physical molecular simulation has become an integral part of chemical research. It is widely applied in a variety of branches of chemistry and has significantly contributed to the development of chemical knowledge. It offers understanding and interpretation of experimental results, semiquantitative predictions for measurable and nonmeasurable properties of substances, and allows the calculation of properties of molecular systems under conditions that are experimentally inaccessible. Yet, molecular simulation is built on a number of assumptions, approximations, and simplifications which limit its range of applicability and its accuracy. These concern the potential-energy function used, adequate sampling of the vast statistical-mechanical configurational space of a molecular system and the methods used to compute particular properties of chemical systems from statistical-mechanical ensembles. During the past half century various methodological ideas to improve the efficiency and accuracy of classical-physical molecular simulation have been proposed, investigated, evaluated, implemented in general simulation software or were abandoned. The latter because of fundamental flaws or, while being physically sound, computational inefficiency. Some of these methodological ideas are briefly reviewed and the most effective methods are highlighted. Limitations of classical-physical simulation are discussed and perspectives are sketched.

Statistical Mechanics Approaches for Studying Temperature and Rate Effects in Multistable Systems

Systems with a multistable energy landscape are widespread in physics, biophysics, technology, and materials science. They are strongly influenced by thermal fluctuations and external mechanical actions that can be applied at different rates, moving the system from equilibrium to non-equilibrium regimes. In this paper, we focus on a simple system involving a single breaking phenomenon to describe the various theoretical approaches used to study these problems. To begin with, we propose the exact solution at thermodynamic equilibrium based on the calculation of the partition function without approximations. We then introduce the technique of spin variables, which is able to simplify the treatment even for systems with a large number of coordinates. We then analyze the energy balance of the system to better understand its underlying physics. Finally, we introduce a technique based on transition state theory useful for studying the non-equilibrium dynamical regimes of these systems. This method is appropriate for the evaluation of rate effects and hysteresis loops. These approaches are developed for both the Helmholtz ensemble (prescribed extension) and the Gibbs ensemble (applied force) of statistical mechanics. The symmetry and duality of these two ensembles is discussed in depth. While these techniques are used here for a simple system with theoretical purposes, they can be applied to complex systems of interest for several physical, biophysical, and technological applications.

Statistical field theory of polarizable polymer chains with nonlocal dipolar interactions.

The electromechanical response of polymeric soft matter to applied electric fields is of fundamental scientific interest as well as relevant to technologies for sensing and actuation. Several existing theoretical and numerical approaches for polarizable polymers subject to a combined applied electric field and stretch are based on discrete monomer models. In these models, accounting for the interactions between the induced dipoles on monomers is challenging due to the nonlocality of these interactions. On the other hand, the framework of statistical field theory provides a continuous description of polymer chains that potentially enables a tractable way to account for these interactions. However, prior formulations using this framework have been restricted to the case of weak anisotropy of the monomer polarizability. This paper formulates a general approach based in the framework of statistical field theory to account for the nonlocal nature of the dipolar interactions without any restrictions on the anisotropy or nonlinearity of the polarizability of the monomer. The approach is based on three key elements: (1) the statistical field theory framework, in which the discrete monomers are regularized to a continuous dipole distribution, (2) a replacement of the nonlocal dipole-dipole interactions by the local electrostatics partial differential equationwith the continuous dipole distribution as the forcing, and (3) the use of a completely general relation between the polarization and the local electric field. Rather than treat the dipole-dipole interactions directly, the continuous description in the field theory enables the computationally tractable nonlocal-to-local transformation. Further, it enables the use of a realistic statistical-mechanical ensemble wherein the average far-field applied electric field is prescribed, rather than prescribing the applied field at every point in the polymer domain. The model is applied, using the finite element method, to study the electromechanical response of a polymer chain in the ensemble with fixed far-field applied electric field and fixed chain stretch. The nonlocal dipolar interactions are found to increase, over the case where dipole-dipole interactions are neglected, the magnitudes of the polarization and electric field by orders of magnitude as well as significantly change their spatial distributions. Next, the effect of the relative orientation between the applied field and the chain on the local electric field and polarization is studied. The model predicts that the elastic response of the polymer chain is linear, consistent with the Gaussian approximation, and largely unchanged by the orientation of the applied electric field, though the polarization and local electric field distributions are significantly impacted.

Statistical structural inference from edge weights using a mixture of gamma distributions

Abstract The inference of reliable and meaningful connectivity information from weights representing the affinity between nodes in a graph is an outstanding problem in network science. Usually, this is achieved by simply thresholding the edge weights to distinguish true links from false ones and to obtain a sparse set of connections. Tools developed in statistical mechanics have provided particularly effective ways to locate the optimal threshold so as to preserve the statistical properties of the network structure. Thermodynamic analogies together with statistical mechanical ensembles have been proven to be useful in analysing edge-weighted networks. To extend this work, in this article, we use a statistical mechanical model to describe the probability distribution for edge weights. This models the distribution of edge weights using a mixture of Gamma distributions. Using a two-component Gamma mixture model with components describing the edge and non-edge weight distributions, we use the Expectation–Maximization algorithm to estimate the corresponding Gamma distribution parameters and mixing proportions. This gives the optimal threshold to convert weighted networks to sets of binary-valued connections. Numerical analysis shows that it provides a new way to describe the edge weight probability. Furthermore, using a physical analogy in which the weights are the energies of molecules in a solid, the probability density function for nodes is identical to the degree distribution resulting from a uniform weight on edges. This provides an alternative way to study the degree distribution with the nodal probability function in unweighted networks. We observe a phase transition in the low-temperature region, corresponding to a structural transition caused by applying the threshold. Experimental results on real-world weighted and unweighted networks reveal an improved performance for inferring binary edge connections from edge weights.

Stable recursive auxiliary field quantum Monte Carlo algorithm in the canonical ensemble: Applications to thermometry and the Hubbard model.

Many experimentally accessible, finite-sized interacting quantum systems are most appropriately described by the canonical ensemble of statistical mechanics. Conventional numerical simulation methods either approximate them as being coupled to a particle bath or use projective algorithms which may suffer from nonoptimal scaling with system size or large algorithmic prefactors. In this paper, we introduce a highly stable, recursive auxiliary field quantum Monte Carlo approach that can directly simulate systems in the canonical ensemble. We apply the method to the fermion Hubbard model in one and two spatial dimensions in a regime known to exhibit a significant "sign" problem and find improved performance over existing approaches including rapid convergence to ground-state expectation values. The effects of excitations above the ground state are quantified using an estimator-agnostic approach including studying the temperature dependence of the purity and overlap fidelity of the canonical and grand canonical density matrices. As an important application, we show that thermometry approaches often exploited in ultracold atoms that employ an analysis of the velocity distribution in the grand canonical ensemble may be subject to errors leading to an underestimation of extracted temperatures with respect to the Fermi temperature.

H-Theorem in an Isolated Quantum Harmonic Oscillator.

We consider the H-theorem in an isolated quantum harmonic oscillator through the time-dependent Schrödinger equation. The effect of potential in producing entropy is investigated in detail, and we found that including a barrier potential into a harmonic trap would lead to the thermalization of the system, while a harmonic trap alone would not thermalize the system. During thermalization, Shannon entropy increases, which shows that a microscopic quantum system still obeys the macroscopic thermodynamics law. Meanwhile, initial coherent mechanical energy transforms to incoherent thermal energy during thermalization, which exhibiting the decoherence of an oscillating wave packet featured by a large decreasing of autocorrelation length. When reaching thermal equilibrium, the wave packet comes to a halt, with the density distributions both in position and momentum spaces well-fitted by a microcanonical ensemble of statistical mechanics.

Physics-constrained Bayesian inference of state functions in classical density-functional theory.

We develop a novel data-driven approach to the inverse problem of classical statistical mechanics: Given the experimental data on the collective motion of a classical many-body system, how does one characterize the free energy landscape of that system? By combining non-parametric Bayesian inference with physically motivated constraints, we develop an efficient learning algorithm that automates the construction of approximate free-energy functionals. In contrast to optimization-based machine learning approaches, which seek to minimize a cost function, the central idea of the proposed Bayesian inference is to propagate a set of prior assumptions through the model, derived from physical principles. The experimental data are used to probabilistically weigh the possible model predictions. This naturally leads to humanly interpretable algorithms with full uncertainty quantification of predictions. In our case, the output of the learning algorithm is a probability distribution over a family of free energy functionals, consistent with the observed particle data. We find that surprisingly small data samples contain sufficient information for inferring highly accurate analytic expressions of the underlying free-energy functionals, making our algorithm highly data efficient. In particular, we consider classical particle systems with excluded volume interactions, which are ubiquitous in nature, while being highly challenging in terms of free energy modeling. We validate our approach on the paradigmatic case of one-dimensional fluid and develop inference algorithms for the canonical and grand-canonical statistical-mechanical ensembles. Extensions to higher dimensional systems are conceptually straightforward, while standard coarse-graining techniques allow one to easily incorporate attractive interactions.

Chaitin’s Omega and an algorithmic phase transition

We consider the statistical mechanical ensemble of bit string histories that are computed by a universal Turing machine. The role of the energy is played by the program size. We show that this ensemble has a first-order phase transition at a critical temperature, at which the partition function equals Chaitin’s halting probability Ω. This phase transition has curious properties: the free energy is continuous near the critical temperature, but almost jumps: it converges more slowly to its finite critical value than any computable function. At the critical temperature, the average size of the bit strings diverges. We define a non-universal Turing machine that approximates this behavior of the partition function in a computable way by a super-logarithmic singularity, and discuss its thermodynamic properties. We also discuss analogies and differences between Chaitin’s Omega and the partition function of a quantum mechanical particle, and with quantum Turing machines. For universal Turing machines, we conjecture that the ensemble of bit string histories at the critical temperature has a continuum formulation in terms of a string theory.

Critical Role of Thermal Fluctuations for CO Binding on Electrocatalytic Metal Surfaces.

This work considers the evaluation of density functional theory (DFT) when comparing against experimental observations of CO binding trends on the strong binding Pt(111) and intermediate binding Cu(111) and for weak binding Ag(111) and Au(111) surfaces important in electrocatalysis. By introducing thermal fluctuations using appropriate statistical mechanical NVT and NPT ensembles, we find that the RPBE and B97M-rV DFT functionals yield qualitatively better metal surface strain trends and CO enthalpies of binding for Cu(111) and Pt(111) than found at 0 K, thereby correcting the overbinding by 0.2 to 0.3 eV to yield better agreement with the enthalpies determined from experiment. The importance of dispersion effects are manifest for the weak CO binding Ag(111) and Au(111) surfaces at finite temperatures in which the RPBE functional does not bind CO at all, while the B97M-rV functional shows that the CO–metal interactions are a mixture of chemisorbed and physisorbed species with binding enthalpies that are within ∼0.05 eV of experiment. Across all M(111) surfaces, we show that the B97M-rV functional consistently predicts the correct atop site preference for all metals due to thermally induced surface distortions that preferentially favor the undercoordinated site. This study demonstrates the need to fully account for finite temperature fluctuations to make contact with the binding enthalpies from surface science experiments and electrocatalysis applications.

Systematic derivation of hydrodynamic equations for viscoelastic phase separation

We present a detailed derivation of a simple hydrodynamic two-fluid model, which aims at the description of the phase separation of non-entangled polymer solutions, where viscoelastic effects play a role. It is directly based upon the coarse-graining of a well-defined molecular model, such that all degrees of freedom have a clear and unambiguous molecular interpretation. The considerations are based upon a free-energy functional, and the dynamics is split into a conservative and a dissipative part, where the latter satisfies the Onsager relations and the second law of thermodynamics. The model is therefore fully consistent with both equilibrium and non-equilibrium thermodynamics. The derivation proceeds in two steps: firstly, we derive an extended model comprising two scalar and four vector fields, such that inertial dynamics of the macromolecules and of the relative motion of the two fluids is taken into account. In the second step, we eliminate these inertial contributions and, as a replacement, introduce phenomenological dissipative terms, which can be modeled easily by taking into account the principles of non-equilibrium thermodynamics. The final simplified model comprises the momentum conservation equation, which includes both interfacial and elastic stresses, a convection–diffusion equation where interfacial and elastic contributions occur as well, and a suitably convected relaxation equation for the end-to-end vector field. In contrast to the traditional two-scale description that is used to derive rheological equations of motion, we here treat the hydrodynamic and the macromolecular degrees of freedom on the same basis. Nevertheless, the resulting model is fairly similar, though not fully identical, to models that have been discussed previously. Notably, we find a rheological constitutive equation that differs from the standard Oldroyd-B model. Within the framework of kinetic theory, this difference may be traced back to a different underlying statistical-mechanical ensemble that is used for averaging the stress. To what extent the model is able to reproduce the full phenomenology of viscoelastic phase separation is presently an open question, which shall be investigated in the future.

A statistical mechanics-based investigation of a solute atmosphere interacting with a crack: Solute redistribution and image effects

Although continuum frameworks have been developed in the past for studying how material defects behave in solid solutions, recent studies have demonstrated that the self-stress field of the solutes which enter the chemical potential in the prevalent continuum models must be excluded, and instead the image stresses included, in a correct treatment. Here, working within the framework of Gibbs canonical ensemble of statistical mechanics, and accounting for the image stresses of the solutes arising due to the presence of the boundaries, we examine the interaction between a solute atmosphere and a crack in the presence of an external loading. Relying on a plane strain model, the solute atmosphere is treated as a distribution of misfitting solute rods, infinitely long in the direction normal to the plane of analysis. Since solutes interact through their image stress fields and can move around in the solid, Monte Carlo simulations, fully accounting for both the energy of interaction between the solute rods and their interaction with the external loading, are performed to explore various configurations of the solute rods in the host, and to compute the temperature-dependent ensemble averages of the enthalpy, equilibrium solute distribution, and the energy release rate (ERR) of crack growth. It is shown that solute redistribution induced by a tensile loading of the solid results in an increase in the ERR. Our simulations further interestingly reveal that even in the limit of large temperatures where a nearly uniform solute distribution prevails, the interaction between the collective image stress field of the solutes and the stress field of the loading results in an added contribution to the ERR which, at a fixed average concentration and a fixed loading-induced stress intensity factor, does not fade away with increasing solid dimensions, and further strengthens both with the solute concentration and with the external loading.

Parallel Combination of Inner Capacitance and Ionic Capacitance, Apparently Inconsistent with Stern’s Model

A double layer capacitance (DLC) has mainly been brought about in the Helmholtz layer rather than in the diffuse layer, as was demonstrated with the invariance of DLC to salt concentration, c, less than 0.5 M (M = mol dm−3). The DLC measured here increased with concentrations of KCl and HCl solutions as high as 1 M at a platinum electrode by the ac impedance method. It was represented as a sum of the Helmholtz capacitance and the ionic one which had 0.7 power of the concentrations. The simple addition implies that the Helmholtz contribution and the ionic one should be represented by a parallel combination rather than a series one such as in the Stern model. The disagreement of the experimental values of the DLC with the Gouy–Chapman theory at high concentrations has been conventionally attributed to the effects of packing of ions over their sizes. In this paper, a model of avoiding the packing was introduced, in which ions were distributed in the direction normal to the electrode in the balance of electric motive force and the thermal energy, keeping the uniform distribution on a plane projected to the electrode. The energy balance was taken by using the grand canonical ensemble in statistical mechanics. The ionic contribution had a linear relation with the applied voltage rather than exponential dependence. When a series combination was applied to the Helmholtz capacitance and the ionic one under the condition of difference between the locally anionic DLC and the cationic one, we obtained approximately a parallel combination of the two capacitances because either the anionic or the cationic DLC works predominantly.

Emergence of a thermal equilibrium in a subsystem of a pure ground state by quantum entanglement

By numerically exact calculations of spin-1/2 antiferromagnetic Heisenberg models on small clusters, we demonstrate that quantum entanglement between subsystems $A$ and $B$ in a pure ground state of a whole system $A+B$ can induce thermal equilibrium in subsystem $A$. Here, the whole system is bipartitoned with the entanglement cut that covers the entire volume of subsystem $A$. Temperature ${\cal T}_{A}$ of subsystem $A$ is not a parameter but can be determined from the entanglement von Neumann entropy ${\cal S}_{A}$ and the total energy ${\cal E}_{A}$ of subsystem $A$ calculated for the ground state of the whole system. We show that temperature ${\cal T}_{A}$ can be derived by minimizing the relative entropy for the reduced density matrix operator of subsystem $A$ and the Gibbs state (i.e., thermodynamic density matrix operator) of subsystem $A$ with respect to the coupling strength between subsystems $A$ and $B$. Temperature ${\cal T}_{A}$ is essentially identical to the thermodynamic temperature, for which the entropy and the internal energy evaluated using the canonical ensemble in statistical mechanics for the isolated subsystem $A$ agree numerically with the entanglement entropy ${\cal S}_{A}$ and the total energy ${\cal E}_{A}$ of subsystem $A$.Fidelity calculations ascertain that the reduced density matrix operator of subsystem $A$ for the pure but entangled ground state of the whole system $A+B$ matches, within a maximally $1.5\%$ error in the finite size clusters studied, the thermodynamic density matrix operator of subsystem $A$ at temperature ${\cal T}_{A}$. We argue that quantum fluctuation in an entangled pure state can mimic thermal fluctuation in a subsystem. We also provide two simple but nontrivial analytical examples of free bosons and free fermions for which these statements are exact. We furthermore discuss implications and possible applications of our finding.

Alchemical molecular dynamics for inverse design

We present a molecular dynamics (MD) implementation of an extended statistical mechanical ensemble that includes ‘alchemical’ degrees of freedom describing particle attributes as thermodynamic variables. We demonstrate the use of this alchemical MD method in inverse design simulations of particles interacting via the Oscillating Pair Potential (OPP) and the Lennard–Jones–Gauss potential (LJG) – two general, previously studied models for which phase diagrams are known. We show that alchemical MD can quickly and efficiently optimise pair potentials for target structures within a specified design space in the low-temperature regime, where internal energy adequately represents the features of the alchemical free energy landscape. We show that alchemical MD can be also used to inversely design pair potentials to achieve target materials properties (here, bulk modulus) directly, without explicit knowledge of the structure–property relationship. Alchemical MD can easily be generalised and applied to any target materials properties or structures and used with any differentiable interaction potential.

Squeezed ensemble for systems with first-order phase transitions

All ensembles of statistical mechanics are equivalent in the sense that they give the equivalent thermodynamic functions in the thermodynamic limit. However, when investigating microscopic structures in the first-order phase transition region, one must choose an appropriate statistical ensemble. The appropriate choice is particularly important when one investigates finite systems, for which even the equivalence of ensembles does not hold. We propose a class of statistical ensembles, which always give the correct equilibrium state even in the first-order phase transition region. We derive various formulas for this class of ensembles, including the one by which temperature is obtained directly from energy without knowing entropy. Moreover, these ensembles are convenient for practical calculations because of good analytic properties. We also derive formulas which relate statistical-mechanical quantities of different ensembles, including the conventional ones, for finite systems. The formulas are useful for obtaining results with smaller finite-size effects, and for improving the computational efficiency. The advantages of the squeezed ensembles are confirmed by applying them to the Heisenberg model and the frustrated Ising model.

Thermodynamic derivation of classical density functional theory

Classical density functional theory has evolved into a major branch of statistical and condensed matter physics. The fundamental equation of the equilibrium theory is where is the thermal Helmholtz free energy of the system as a functional of its non-uniform local number density , is the external potential, μ is the uniform chemical potential, and denotes an isothermal functional derivative. This equation implicitly determines , and is ordinarily derived from the grand canonical ensemble (GCE) of statistical mechanics. Here we show that it can also be simply derived from thermodynamics alone, and is therefore not inherently statistical in character or specific to the GCE. This derivation further shows that the isothermal functional derivative of At is identically equal to the isentropic functional derivative of the internal energy Et, which is not immediately obvious in the statistical theory. This development shows that certain aspects of density functional theory are essentially macroscopic rather than statistical in character, thereby clarifying the physical content of the theory and making it more accessible to students and nonspecialists.

A methodology to calculate small-angle scattering profiles of macromolecular solutions from molecular simulations in the grand-canonical ensemble.

The theoretical framework to evaluate small-angle scattering (SAS) profiles for multi-component macromolecular solutions is re-examined from the standpoint of molecular simulations in the grand-canonical ensemble, where the chemical potentials of all species in solution are fixed. This statistical mechanical ensemble resembles more closely scattering experiments, capturing concentration fluctuations that arise from the exchange of molecules between the scattering volume and the bulk solution. The resulting grand-canonical expression relates scattering intensities to the different intra- and intermolecular pair distribution functions, as well as to the distribution of molecular concentrations on the scattering volume. This formulation represents a generalized expression that encompasses most of the existing methods to evaluate SAS profiles from molecular simulations. The grand-canonical SAS methodology is probed for a series of different implicit-solvent, homogeneous systems at conditions ranging from dilute to concentrated. These systems consist of spherical colloids, dumbbell particles, and highly flexible polymer chains. Comparison of the resulting SAS curves against classical methodologies based on either theoretical approaches or canonical simulations (i.e., at a fixed number of molecules) shows equivalence between the different scattering intensities so long as interactions between molecules are net repulsive or weakly attractive. On the other hand, for strongly attractive interactions, grand-canonical SAS profiles deviate in the low- and intermediate-q range from those calculated in a canonical ensemble. Such differences are due to the distribution of molecules becoming asymmetric, which yields a higher contribution from configurations with molecular concentrations larger than the nominal value. Additionally, for flexible systems, explicit discrimination between intra- and inter-molecular SAS contributions permits the implementation of model-free, structural analysis such as Guinier's plots at high molecular concentrations, beyond what the traditional limits are for such analysis.

Isomerization kinetics of flexible molecules in the gas phase: Atomistic versus coarse-grained sampling.

Isomerization kinetics of molecules in the gas phase naturally falls on the microcanonical ensemble of statistical mechanics, which for small systems might significantly differ from the more traditional canonical ensemble. In this work, we explore the examples of cis-trans isomerization in butane and bibenzyl and to what extent the fully atomistic rate constants in isolated molecules can be reproduced by coarse-graining the system into a lower dimensional potential of mean force (PMF) along a reaction coordinate of interest, the orthogonal degrees of freedom acting as a canonical bath in a Langevin description. Time independent microcanonical rate constants can be properly defined from appropriate state residence time correlation functions; however, the resulting rate constants acquire some time dependence upon canonical averaging of initial conditions. Stationary rate constants are recovered once the molecule is placed into a real condensed environment pertaining to the canonical ensemble. The effective one-dimensional kinetics along the PMF, based on appropriately chosen inertia and damping parameters, quantitatively reproduces the atomistic rate constants at short times but deviates systematically over long times owing to the neglect of some couplings between the system and the bath that are all intrinsically present in the atomistic treatment. In bibenzyl, where stronger temperature effects are noted than in butane, the effective Langevin dynamics along the PMF still performs well at short times, indicating the potential interest of this extremely simplified approach for sampling high-dimensional energy surfaces and evaluating reaction rate constants.

Dynamic principle for ensemble control tools.

Dynamical equations describing physical systems in contact with a thermal bath are commonly extended by mathematical tools called "thermostats." These tools are designed for sampling ensembles in statistical mechanics. Here we propose a dynamic principle underlying a range of thermostats which is derived using fundamental laws of statistical physics and ensures invariance of the canonical measure. The principle covers both stochastic and deterministic thermostat schemes. Our method has a clear advantage over a range of proposed and widely used thermostat schemes that are based on formal mathematical reasoning. Following the derivation of the proposed principle, we show its generality and illustrate its applications including design of temperature control tools that differ from the Nosé-Hoover-Langevin scheme.