Simulation of Nuclear Quantum Effects in Condensed Matter Systems via Quantum Baths

Abstract

This paper reviews methods that aim at simulating nuclear quantum effects (NQEs) using generalized thermal baths. Generalized (or quantum) baths simulate statistical quantum features, and in particular zero-point energy effects, through non-Markovian stochastic dynamics. They make use of generalized Langevin Equations (GLEs), in which the quantum Bose–Einstein energy distribution is enforced by tuning the random and friction forces, while the system degrees of freedom remain classical. Although these baths have been formally justified only for harmonic oscillators, they perform well for several systems, while keeping the cost of the simulations comparable to the classical ones. We review the formal properties and main characteristics of classical and quantum GLEs, in relation with the fluctuation–dissipation theorems. Then, we describe the quantum thermostat and quantum thermal bath, the two generalized baths currently most used, providing several examples of applications for condensed matter systems, including the calculation of vibrational spectra. The most important drawback of these methods, zero-point energy leakage, is discussed in detail with the help of model systems, and a recently proposed scheme to monitor and mitigate or eliminate it—the adaptive quantum thermal bath—is summarised. This approach considerably extends the domain of application of generalized baths, leading, for instance, to the successful simulation of liquid water, where a subtle interplay of NQEs is at play. The paper concludes by overviewing further development opportunities and open challenges of generalized baths.

1. Introduction

This paper reviews methods that use generalized thermal baths, that we shall indicate generically as quantum baths. Other expressions, such as generalized or quasiclassical baths, are used in the literature to designate a family of methods revolving around various flavors of dynamics aiming at mimicking quantum properties via generalized Langevin equations. We adopt here the generic expression quantum bath, which is clear enough to specify the field without entering into the technical details separating these schemes to approximate nuclear quantum effects (NQEs) and, in particular, atomic zero-point energy (ZPE).

It is well-recognised that NQEs cannot be neglected for systems at low temperatures or at high pressure, where they can have, for example, a strong influence on phase transitions [1,2,3]. Nuclear quantum effects, however, also often come into play at ambient conditions, typically when the thermal energy of chemical bonds in molecules is comparable or smaller than the associated vibrational ZPE or when characteristic interatomic distances in the system are comparable to the de Broglie wavelength of the nuclei. These conditions are met, in particular, in systems containing light atoms and most notably hydrogen [4,5], which is a ubiquitous constituent of inorganic compounds and a basic element in biological systems. Thus, nuclear quantum effects influence, for example, the stability of crystal polymorphs of pharmaceutical interest, the spectroscopy of ice and water, or enzymatic reactions in living organisms [6,7,8,9]. NQEs are also required to reveal the effect of isotopic substitutions on equilibrium properties, as a classical description of the nuclei predicts identical statistical properties for all isotopes. The exact simulation of these effects is, however, a challenging problem due to the large (sometimes unachievable) computational cost of available methods. Time-independent statistical properties of quantum nuclei at thermal equilibrium can be obtained using the so-called classical isomorphism of path integrals (see Section 6), in which a quantum particle is represented as a collection of beads interacting via an effective classical-like potential. While this approach guarantees theoretical convergence to the exact quantum result in the limit of an infinite number of beads, it necessitates in practice to simulate a number of degrees of freedom that can become prohibitive, in particular when first-principles methods are used for the interatomic interactions. The situation is even more problematic for time-dependent properties for which exact simulation methods scale exponentially with the number of degrees of freedom. Different approaches have been proposed, such as in [10,11], each with their own strengths and weaknesses; but, in spite of several interesting results, no clear, practical and general reference method has emerged.

In the following, we focus on the family of approximate methods for simulating quantum nuclear effects that rely on the use of quantum baths. Although these methods can substantially differ from one another in their practical implementation, they all aim at introducing quantum statistical properties into classical trajectory-based dynamics. This is achieved via a stochastic evolution that takes the form of a generalized Langevin Equation (GLE), in which the quantum Bose–Einstein energy distribution is enforced instead of equipartition, by means of an appropriate tuning of the random and friction forces. The main advantage of these approaches is that the number of degrees of freedom in the system and the deterministic part of the propagation remain classical, thus enabling them to mimic quantum properties at a numerical cost directly comparable to that of standard Langevin dynamics. Contrary to standard Langevin, however, quantum baths exhibit a complex non-Markovian dynamics whose formal properties are not obvious. In particular, while these methods deliver remarkable performances in a growing number of applications, they cannot be formally justified except for harmonic systems. This intriguing scenario motivates our approach for this review. In Section 2, we start by stating some fundamental results in the statistical mechanics of classical and quantum systems, the fluctuation–dissipation theorems (FDT). We then introduce the standard and generalized Langevin equations and set the stage for the following developments. This is done, in particular, by recalling how the GLE can be derived starting from the physical picture of a generic system (not necessarily harmonic) bilinearly coupled to an environment modelled as a set of harmonic oscillators. We shall consider both classical and quantum versions of this system-bath model and describe how the quantum version can be used, via an appropriate quasiclassical limit, to motivate the quantum bath algorithms on which we focus in this review. While more general models for the bath are possible (including, in particular, anharmonicity) these are not particularly relevant for this work and will therefore not be discussed, nor included in examples. Section 3 focuses on using the quantum FDT to build GLEs that combine classical deterministic evolution with quantum stochastic behaviour to mimic NQEs. In particular, we introduce two approaches, the Quantum Thermal Bath [12] (QTB) and the Quantum Thermostat [13,14] (QT) that have recently attracted considerable interest. Section 4 aims at understanding more in depth the capabilities of the quantum thermal bath via the detailed analysis of simple low-dimensional systems that enable a comparison with exact quantum results. In particular, the most important pathology of this type of methods, namely the ZPE leakage, is illustrated on an elementary example. Although this section focuses on the QTB, which is formally simpler and has fewer parameters than the QT, the general observations apply to both frameworks. In Section 5, we review the capabilities of quantum bath methods by summarizing some significant applications to realistic models of condensed matter systems. Remarkably, these techniques prove useful not only in the simulation of static equilibrium properties of quantum nuclei, but also in the calculation of vibrational spectra [15], including subtle anharmonic features that elude other trajectory-based approximations to quantum dynamics [16]. In Section 6, we complete the overview of important applications of thermal baths by showing how they can be used to improve the convergence of calculations based on the path integral formalism. We then introduce in Section 7 a recent development that, by exploiting the quantum fluctuation–dissipation theorem, enables to monitor and efficiently compensate ZPE leakage. This approach, known as the adaptive quantum thermal bath [17] (adQTB) extends considerably the domain of application of the QTB, leading, for example, to the successful simulation of liquid water [18]. The review ends with a discussion of remaining limitations, open questions and possible future developments of quantum baths as effective and accurate tools for the simulation of condensed phase systems.

We conclude this Introduction by noting that GLEs have a long and rich history, with applications in a very broad range of domains that we do not account for here. The literature on quantum baths is also rapidly increasing. While we have tried to provide appropriate references throughout the text, we have chosen to present the material focusing on results and derivations that seemed more directly related to the methods at the core of this work, intended as a relatively self-contained introduction to a still-developing set of exciting new algorithms for growing classes of quantum problems in physics and chemistry.

2. The Langevin Equation as a Thermal Bath

Originally introduced as a description of the Brownian motion, the Langevin equation is widely used in molecular dynamics simulations as a practical way to implement the canonical ensemble sampling. Starting from the microcanonical ensemble for a closed system, the canonical distribution is obtained by considering a small portion of the system that can exchange energy with its environment, which is referred to as the bath. In order to avoid the explicit representation of the bath degrees of freedom, the Langevin equation describes the bath via the combination of a random force and a friction mechanism (see Figure 1), that both act on the system degrees of freedom.

Since any stable physical system at equilibrium fulfills the fluctuation–dissipation theorem (FDT), we start by recalling this pivotal result in statistical mechanics, which will be used all along this review. We then introduce the Langevin equation and its non-Markovian generalization and recall how they can be derived from an explicit harmonic model for the bath. We end this section by reviewing some instances of the generalized Langevin approach in the literature, pointing out more specifically some developments that are relevant in the perspective of the implementation of quantum bath methods to approximate NQEs, as presented in Section 3.

2.1. Linear Response and Fluctuation–Dissipation Theorems

Linear response theory provides useful relations characterizing dynamical properties under external perturbations via averages over the equilibrium probability of the system, e.g., thermal probability density. We first introduce the complex admittance or susceptibility $\chi \left(\omega \right)$ that characterizes the linear response of a system subject to a perturbative force $F\left(t\right)=\mathrm{Re}\left[F_0{e}^{i\omega t}\right]$ along a given direction x. At first order in the perturbation, the change in the velocity along x induced by the perturbation (here and in the remainder of the section, we use one-dimensional notations; the generalization to higher dimension is straightforward) reads:

$$\Delta v\left(t\right)=\mathrm{Re}\left[F_0\chi \left(\omega \right)e^{i\omega t}\right]$$

This relation can be considered as a definition for $\chi \left(\omega \right)$. For a classical system at thermal equilibrium at temperature T, it can be shown that the following fluctuation–dissipation theorem (FDT) holds [19,20]:

$$C_{vv}\left(\omega \right)=2k_{B}T\mathrm{Re}\left[\chi \left(\omega \right)\right]$$

(1)

The FDT as derived by Kubo in [19] is actually more general than Equation (1) and can be used to relate the cross-correlation spectrum $C_{AB}\left(\omega \right)$ of two arbitrary observables A and B to the corresponding linear response function. In our case, Equation (1) relates the real part of the susceptibility (that is dissipative) to the equilibrium velocity autocorrelation spectrum (that characterizes fluctuations):

$$C_{vv}\left(\omega \right)={\int }_{-\infty }^{\infty }dt\langle v\left(0\right)v\left(t\right)\rangle {\mathrm{e}}^{-i\omega t}$$

(2)

The FDT characterizes the frequency distribution of energy at equilibrium, in its classical Formulation (1), it expresses the equipartition of energy in the system. However, a slightly different FDT can be derived for quantum systems [19,21]:

$$C_{vv}^s\left(\omega \right)=2\theta (\omega ,T)\mathrm{Re}\left[\chi \left(\omega \right)\right]\mathrm{with}\theta (\omega ,T)=\frac{\hslash \omega }{2}coth\left(\frac{\hslash \omega }{2k_{B}T}\right)$$

(3)

Here the classical thermal energy $k_{B}T$ has been replaced by the quantum Bose–Einstein distribution $\theta (\omega ,T)$, and the autocorrelation spectrum that is used is now the Fourier transform of the symmetrized correlation function:

$$C_{vv}^s\left(\omega \right)={\int }_{-\infty }^{\infty }dt\mathrm{Tr}\left[\frac{1}{2}{\rho }_{eq}\left\{v\left(0\right)v\left(t\right)+v\left(t\right)v\left(0\right)\right\}\right]{\mathrm{e}}^{-i\omega t}$$

with ${\rho }_{eq}$ the equilibrium density operator. $v\left(t\right)$ now designates the time-dependent velocity operator in the Heisenberg picture. The definition of the linear susceptibility is also modified in the quantum case and is related to the change induced by the perturbative force on the density operator $\rho$. The fluctuation–dissipation theorem plays a key role in the generalized Langevin equations used to sample the classical canonical ensemble, as we will develop in the following sections. More recently, the quantum version of this theorem has also been employed to design the trajectory-based approximations for the simulation of NQEs at the core of this review.

2.2. The Langevin Equation in Classical Molecular Dynamics Simulations

In its basic form, the Langevin equation of motion writes:

$$m\ddot{x}=-\frac{dV}{dx}-m\gamma \dot{x}+R\left(t\right)$$

(4)

where $R\left(t\right)$ is a random force with zero mean (usually taken from a Gaussian distribution) and characterized by a white noise correlation function:

$$\langle R\left(t\right)R(t+\tau )\rangle =2\gamma m{k}_{B}T\delta \left(\tau \right)$$

(5)

The random and friction forces model the effect of a thermal bath (see Figure 1) and ensure the canonical statistics, including the energy equipartition theorem and the fluctuation–dissipation theorem as stated in Section 2.1. The constant $\gamma$ is the strength of the coupling of the thermal bath with the simulated system. Its inverse $1/\gamma$ provides an indication of the relaxation time of the system, that is, the time it takes for the random force and the dissipation to equilibrate.

Used as a thermostat in molecular dynamics simulations, the Langevin Equation (4) provides, for any $\gamma$, an exact sampling of the canonical Boltzmann distribution. More precisely, it can be proved that the probability distribution converges exponentially to the invariant one starting from a wide range of different initial conditions [22], thus avoiding temperature oscillations or other spurious effects that other thermostats, such as Nosé-Hoover, can generate. Dynamical properties, on the other hand, are affected by the coupling with the bath and depend on the value of $\gamma$. In particular, the damping term generally broadens the peaks corresponding to modes in the vibrational spectra over a typical width $\gamma$.

2.3. Generalized Langevin Equation (GLE)

The standard (Markovian) Langevin Equation (4) can be generalized in the following manner:

$$m\ddot{x}=-\frac{dV}{dx}-{\int }_{-\infty }^{t}dsK(t-s)\dot{x}\left(s\right)+R\left(t\right)$$

(6)

This equation is non-Markovian since the friction force does not only depend on the velocity $\dot{x}$ at time t, but is rather expressed as an integral over past values of $\dot{x}$ through the memory kernel $K\left(\tau \right)$, which determines the actual dependence. In order to enforce the canonical ensemble statistics and the equipartition of energy, the random force should then be related to $K\left(\tau \right)$ by the following relation:

$$\langle R\left(t\right)R(t+\tau )\rangle =k_{B}TK\left(\tau \right)$$

(7)

which reduces to (5) in the Markovian limit $K\left(\tau \right)=2m\gamma \delta \left(\tau \right)$. As the bath is supposed to be at equilibrium, the time correlation function $\langle R\left(t\right)R(t+\tau )\rangle$ does not depend on the arbitrary time origin t. Defining the Fourier transformed quantities (here, we modify the memory kernel, which is in principle defined for non-negative times, by symmetrization as $K(-\tau )=K\left(\tau \right)$):

$$\begin{array}{ccc}\hfill \tilde{K}\left(\omega \right)& =& {\int }_{-\infty }^{\infty }d\tau K\left(\tau \right){\mathrm{e}}^{-i\omega \tau }\hfill \\ \hfill C_{RR}\left(\omega \right)& =& {\int }_{-\infty }^{\infty }d\tau \langle R\left(t\right)R(t+\tau )\rangle {\mathrm{e}}^{-i\omega \tau },\hfill \end{array}$$

yields the following relation:

$$C_{RR}\left(\omega \right)=k_{B}T\tilde{K}\left(\omega \right).$$

(8)

Equation (8)—or its time-domain version, (7)—follows from the equilibrium between the system and the thermal bath, and it is generally referred to as a fluctuation–dissipation theorem (FDT), similarly to Equation (1). The latter, however, is a general result of linear response theory that involves only intrinsic observables of the system (its velocity autocorrelation and its linear response function) whereas Equations (7) and (8), characterize the system–bath interactions via the friction and random forces in the particular context of the generalized Langevin dynamics. As pointed out by Kubo, Toda and Hatsishume [20], Equation (1) can be regarded as more fundamental, and following these authors, we will refer to it as the first-kind FDT, while Equations (7) and (8) (that can be derived from it as a corollary when the potential $V\left(x\right)$ is harmonic) will be designated as second-kind FDT.

2.4. Coupling to a Harmonic Bath

The generalized Langevin equation (6) appears in a variety of contexts. In particular, a slightly different form of GLE can be derived under very general assumptions using projection operator techniques [20,23], an approach that has recently lead to interesting developments [24,25]. However, the damping and random forces obtained in this formalism are generally not explicit. In this section, we focus on an alternative, more practical and explicit approach: we derive the GLE from the description of a system bilinearly coupled to a bath of independent harmonic oscillators. In passing, we note that, in principle, the bath can be constructed in different ways, for instance by including non-harmonic oscillators. However, the bath of harmonic oscillators is at the same time easy-to-treat analytically and can describe all features that are relevant for the physical system. The extended Hamiltonian for the bath and the physical system, which can be harmonic or not, (the latter case being the most physically relevant one) is then:

$$H=\frac{p^2}{2m}+V\left(x\right)+\sum _j\left[\frac{p_j^2}{2m_j}+\frac{1}{2}{m}_j{\omega }_j^2{(q_j-x)}^2\right]$$

(9)

where x and p are the system coordinate and momentum, while $q_j$, $p_j$ are an ensemble of harmonic oscillators that constitute the thermal bath. The coupling to the system arises from the fact that the equilibrium position of the bath oscillators is taken at $q_j-x$ and therefore changes in time with the variations of x. It is then easily shown that the equation of motion for the system coordinate can be recast into a GLE of the form of (6), where the memory kernel $K\left(\tau \right)$ is defined as:

$$\begin{array}{c}K\left(\tau \right)=\sum _j{m}_j{\omega }_j^{2}cos\left({\omega }_j\tau \right),\\ \mathrm{or}\mathrm{equivalently},\tilde{K}\left(\omega \right)=\pi \sum _j{m}_j{\omega }_j^2\left[\delta (\omega -{\omega }_j)+\delta (\omega +{\omega }_j)\right]\end{array}$$

(10)

The damping force is therefore generally non-Markovian with a kernel that depends on the spectral density of the bath oscillators. In particular, in the limit of an infinite number of bath oscillators with a spectral density proportional to $1/{\omega }^2$, $\tilde{K}\left(\omega \right)$ becomes constant and the Markovian limit is recovered. The random force $R\left(t\right)$ is related to the initial configuration of the bath at an initial time $t_0$:

$$R\left(t\right)=-K(t-t_0)x\left(t_0\right)+\sum _j\left[m_j{\omega }_j^2{q}_j\left(t_0\right)cos\left[{\omega }_j(t-t_0)\right]+{\omega }_j{p}_j\left(t_0\right)sin\left[{\omega }_j(t-t_0)\right]\right]$$

(11)

If $t_0$ is in a distant past, the first term vanishes. Moreover, if the bath oscillators are assumed to be initially (i.e., at $t_0$) in thermal equilibrium, the random force autocorrelation function can be shown to obey the second-kind FDT of Equations (7) and (8).

Interestingly, if the system–bath Hamiltonian (9) is treated quantum-mechanically, a GLE formally identical to Equation (6) can be derived for the time-dependent operator x in the Heisenberg picture. The expression of the memory kernel $K\left(\tau \right)$ and the random force $R\left(t\right)$ are the same as in the classical case (only $R\left(t\right)$ is now an operator), but the second-kind FDT becomes (compare to (8)):

$$C_{RR}^s\left(\omega \right)={\int }_{-\infty }^{\infty }dt\frac{1}{2}〈R\left(t\right)R(t+\tau )+R(t+\tau )R\left(t\right)〉{\mathrm{e}}^{-i\omega \tau }=\theta (\omega ,T)\tilde{K}\left(\omega \right)$$

(12)

where, similarly to what was reported for the first-kind FDT (3), the symmetrised correlation function is used and the classical average thermal energy $k_{B}T$ is replaced by its quantum counterpart $\theta (\omega ,T)$.

In the developments above, we followed mostly the work and formalism of Ford and coworkers in [26,27]. In their studies, the authors point out that the quite simple harmonic bath model is more general than it appears. First, they show that, for the GLE (6) to be physically meaningful, the memory kernel $\tilde{K}\left(\omega \right)$ should be positive for all frequencies. (Ford and coauthors actually consider more generally the Laplace transform:

$$\mu \left(z\right)={\int }_0^{\infty }dt{\mathrm{e}}^{-iz\tau }K\left(\tau \right),$$

that is related to the Fourier transform of the (symmetrized) memory kernel by $\tilde{K}\left(\omega \right)=2\mathrm{Re}\left[\mu (\omega +i{0}^{+})\right]$. They show that $\mu \left(z\right)$ is analytical in the upper half-plane (z is a complex variable), and belongs to the class of the so-called positive real functions that possess several specific mathematical properties that the authors then use in their argumentation on the generality of the harmonic bath model [27]). Now, according to Equation (10), with an appropriate choice of the frequencies ${\omega }_j$ (potentially in infinite number), the harmonic bath can be tailored to produce any arbitrary positive memory kernel $\tilde{K}\left(\omega \right)$. Therefore, it can be regarded as a very generic prototype for the GLE (6) and the authors proceed by pointing out the relation between the harmonic bath of Equation (9) and some of the pre-existing models. Interestingly, in a former work [28], Cortes et al. had shown that the simple and general second-kind FDT in (12) might not be valid when considering forms of the system-bath coupling more complex than the bilinear interaction. In that case, and for a quantum mechanical treatment (in the classical case the FDT remains of a simple form), the relation between the dissipation and the fluctuations depends on the details of the actual system–bath coupling. In particular, Cortes et al. derive an expression for the FDT in the case of a quadratic interaction.

2.5. Quasiclassical Limit of Harmonic Bath

The Hamiltonian (9) or closely related harmonic bath models are common starting points to derive GLE also when looking at quasiclassical approximations for quantum dynamical properties. In general, these developments start from a fully quantum dynamical description of a generic system coupled to the harmonic bath and then approximate the dynamics in some flavor of the classical limit. The specific form of the results depends both on the procedure enforced to impose the classical limit (e.g., high temperature, formal limit for $\hslash \to 0$, stationary phase approximation of the propagator…), and on the adopted representation of quantum mechanics. For example, starting from the Wigner representation of quantum averages, nuclear quantum effects can be approximated using generalized Langevin equations where the friction and random forces are kept Markovian but become position-dependent [29,30]. Furthermore, one can take the classical limit for all degrees of freedom (semiclassical methods) or exploit knowledge of the exact solution of the harmonic part of the problem to derive a mixed quantum–classical scheme in which the bath is described quantum mechanically, while the dynamics of the system, including the coupling to the harmonic bath, is treated classically (quasiclassical or mixed quantum–classical limit). Mixed quantum–classical limits, however, are delicate. For instance, Bader and Berne [31] considered the—apparently simple—case of a harmonic system linearly coupled to a harmonic bath to investigate the problem of vibrational energy relaxation. They recovered the well-known result that, for this kind of system, quantum and classical dynamics lead to the same time-evolution for averages and showed that the quantum relaxation time can be evaluated from a purely classical calculation. However, they also proved that—when a mixed quantum–classical dynamics is imposed on the system—the result for the relaxation time is incorrect and different depending on whether the bath is treated quantum mechanically and the system classically or vice-versa.

A careful analysis of the mixed quantum–classical limit that nicely sets the stage for the methods presented in the next chapter was performed by Schmid [32]. Schmid focused on the time evolution of the density matrix of a generic quantum system linearly coupled to a bath of quantum harmonic oscillators (as in (9)) as written in the path integral formulation of quantum mechanics (see also Section 6). In the coordinate representation, the density matrix determines the probability to find the system at a given position at time $t=0$ and at another assigned position at time t. Within the path integral formalism, the density matrix can be computed as a functional integral over the set of all possible paths connecting the initial and final positions of the system and bath in the assigned time t. The harmonic nature of the bath enables an explicit expression for the marginal density matrix of the system to be obtained by integrating over the paths of the bath. This still leaves an integral over the system’s paths to perform. When the coupling with the bath is sufficiently strong, quantum coherence effects are strongly suppressed. The integral can then be approximated to argue that the relevant paths for the system become stochastic trajectories of the form:

$$m\ddot{x}=-\frac{dV}{dx}-m\gamma \dot{x}+R\left(t\right)$$

(13)

where $R\left(t\right)$ is a Gaussian stochastic process such that, in Fourier space,

$$C_{RR}\left(\omega \right)=2m\gamma \theta (\omega ,T)$$

(14)

The result above is obtained by introducing a non-trivial (and somewhat heuristic) probabilistic interpretation of the path integral representation of the density matrix and enforcing a quasiclassical approximation based on a truncation of an expansion of the potential that is not always easy to justify. Although Equations (13) and (14) seem formally identical to the quantum generalized Langevin Equations (6) and (12) in the particular case of a Markovian friction kernel $K\left(\tau \right)=2m\gamma \delta \left(\tau \right)$, a fundamental difference has been introduced with the quasiclassical approximation: whereas in Equation (6), x referred to the quantum position operator (time-dependent in the Heisenberg picture), in Equation (13), x is now the position of a classical-like system, following a stochastic trajectory. This result thus combines a fully classical evolution of the system (given by the deterministic force $-dV/dx$) with a colored thermostat that injects quantum properties, and in particular zero point energy, in the dynamics. This early approach bears remarkable similarities with the two more recent methods detailed in the next section.

3. Enforcing Quantum Statistics via Generalized Bath

In this section, we focus on two approaches, the Quantum Thermal Bath (QTB) and the Quantum Thermostat (QT) that were developed independently in 2009 by two different groups. Though these methods differ in their practical implementation, they both prescribe a form of GLE that combines a classical evolution for the system and the coupling to a quantum bath designed to induce relevant quantum properties (e.g., zero-point energy) in the otherwise classical evolution of the system.

3.1. The Quantum Thermal Bath

The quantum thermal bath was proposed by Dammak and coworkers [12] to approximate NQEs in general molecular dynamics simulations via a generalized Langevin equation. They use a Markovian friction kernel and a colored random force $R\left(t\right)$ with an autocorrelation spectrum which is tailored according to the energy distribution $\theta (\omega ,T)$ of the quantum harmonic oscillator. The equations of the resulting dynamics are identical to those derived from the quasiclassical approximation by Schmid [32]. For one-dimensional systems, Equations (13) and (14) hold; their generalization to multiple atomic degrees of freedom is straightforward. However, contrary to Schmid, who aimed at modelling the dynamics of a system coupled to an actual physical bath, in the QTB, the Langevin equation is simply used as a thermostat in order to enforce the quantum statistical distribution of energy (including zero-point energy effects). The expression of the random force autocorrelation, Equation (14), is universal, i.e., system-independent, so that the method does not require a priori knowledge of the system or of its detailed spectral density. Heuristically, this form of GLE tends to thermalize the vibration modes of the system with the average energy distribution of quantum harmonic oscillators in equilibrium at temperature T, hence the lexical origin of the quantum thermal bath.

Plots of $\theta (\omega ,T)$ and its temperature derivative, which enters the definition of key quantities such as the heat capacity, are provided in Figure 2 for three typical frequencies $\omega$. One can easily check that:

$$\theta (\omega ,T)=\frac{\hslash \omega }{2}\coth \left(\frac{\hslash \omega }{2k_{B}T}\right)=\frac{\hslash \omega }{2}\left[\frac{2k_{B}T}{\hslash \omega }+\frac{\hslash \omega }{6k_{B}T}-o{\left(\frac{\hslash \omega }{2k_{B}T}\right)}^3\right]\mathrm{for}\frac{\hslash \omega }{2k_{B}T}\ll 1$$

(15)

Therefore in the classical high-temperature limit, i.e., when $k_{B}T\gg \hslash \omega$ for all relevant frequencies, $\theta (\omega ,T)$ approaches its classical value $k_{B}T$, and substituting in Equation (14) shows that the QTB reduces to a classical Langevin thermostat. However, the typical temperature where the classical limit is reached can be rather high (see Figure 2). Even more delicate are quantities such as the heat capacity, i.e., the temperature derivative of the average energy that for a harmonic oscillator is given by $d\theta (\omega ,T)/dT$ which shows a very slow convergence to the classical value.

Apart from a slight dependence on the friction parameter $\gamma$ that is further analyzed in Section 4.1, the probability distributions obtained in QTB simulations are exact for harmonic systems, by construction of the method. For more general cases, it can provide a good approximation to zero-point energy effects in realistic systems, under some conditions that are further developed in Section 4, Section 5 and Section 7.

3.2. The Quantum Thermostat

Independently from the work by Dammak and coworkers, Ceriotti, Bussi and Parrinello also proposed to approximate NQEs via a non-Markovian GLE. In two papers in 2009 [13,14], they introduced a general framework to implement and propagate GLEs of the form of (6) using auxiliary variables. The authors aim at thermalizing different degrees of freedom with distinct characteristic frequencies at different effective temperatures. This tool can be used for different purposes: in particular, the authors exploit it to build a quantum thermostat (QT) and thermalize vibrational modes with an effective temperature that includes a zero-point motion contribution [14], similarly to the QTB. They also propose to apply this method in the context of Car-Parrinello molecular dynamics [33], where the ionic degrees of freedom should evolve at a given (physical) temperature while the faster electronic degrees of freedom should oscillate around the ground state that corresponds to the instantaneous ionic configuration, at much smaller temperatures [13]. The authors successfully applied the colored bath thermostat on a model polarizable system and verified that the fast degrees of freedom (i.e., the motion of the centers of mass of the shells with respect to their nuclei) are correctly thermalized. More generally, they showed that the colored bath thermostat can be adapted to specific cases once the density of vibrational states is known and that it allows a faster thermalization than the Nosé–Hoover chains [13,22].

Starting from the Orstein–Uhlenbeck theory of stochastic processes, Ceriotti and coworkers show how a non-Markovian GLE can be implemented via a Markovian Langevin dynamics for an extended set of variables of which two are the physical position and momentum (for a one-dimensional physical system) and N are additional auxiliary momenta [34]. In this way, the colored thermostat GLE can be efficiently applied by integrating the Markovian dynamics in the $(N+2)$-dimensional space via usual molecular dynamics techniques. This dynamics is governed by the following equations:

$$\left[\begin{array}{c}\dot{x}\\ \dot{p}\\ \dot{\mathbf{s}}\end{array}\right]=-\left[\begin{array}{ccc}0& -1/m& \mathbf{0}\\ V^{\prime }\left(x\right)& a_{pp}& {\mathbf{a}}_{\mathbf{p}}^{\mathbf{T}}\\ \mathbf{0}& {\overline{\mathbf{a}}}_{\mathbf{p}}& \mathbf{A}\end{array}\right]\left[\begin{array}{c}x\\ p\\ \mathbf{s}\end{array}\right]+\left[\begin{array}{ccc}0& 0& \mathbf{0}\\ 0& b_{pp}& {\mathbf{b}}_{\mathbf{p}}^{\mathbf{T}}\\ \mathbf{0}& {\overline{\mathbf{b}}}_{\mathbf{p}}& \mathbf{B}\end{array}\right]\left[\begin{array}{c}0\\ r_p\\ \mathbf{r}\end{array}\right]$$

(16)

where $\mathbf{s}$ and $\dot{\mathbf{s}}$ are the N-component vectors of the auxiliary momenta and their time derivatives, $\mathbf{A}$ and $\mathbf{B}$ are $(N\times N)$ matrices of adjustable parameters that characterize respectively the friction and the random force memory kernels, together with the N-component parameter vectors ${\mathbf{a}}_{\mathbf{p}}$, ${\overline{\mathbf{a}}}_{\mathbf{p}}$ and ${\mathbf{b}}_{\mathbf{p}}$, ${\overline{\mathbf{b}}}_{\mathbf{p}}$. $V^{\prime }\left(x\right)$ is the spatial derivative of the potential energy and the $(N+1)$-component vector $(r_p,\mathbf{r})$ contains normalized white noise variables: $\langle r_i\left(t\right)r_j\left(t^{\prime }\right)\rangle ={\delta }_{ij}\delta (t-t^{\prime })$. The dynamics in the extended $(x,p,\mathbf{s})$-space is therefore Markovian and it can be solved analytically for a harmonic potential $V\left(x\right)=m{\omega }^2{x}^2/2$. In more general non-linear cases, a numerical algorithm based on the symmetric Trotter decomposition of the Liouvillian can be used, where the linear drift and the random force terms in $(\dot{p},\dot{\mathbf{s}})$ are integrated exactly according to the theory of Ornstein–Uhlenbeck processes, whereas the evolution of x and the non-linear part of $\dot{p}$ are integrated via a velocity Verlet scheme [34].

As described in detail in appendix in [34], by integrating away the auxiliary momenta $\mathbf{s}$, Equation (16) reduces to a GLE in the form of (6) for the physical degree of freedom x, with a memory kernel:

$$K\left(\tau \right)=2a_{pp}\delta \left(\tau \right)-{{\mathbf{a}}_{\mathbf{p}}}^T{e}^{-\left|\tau \right|\mathbf{A}}{\overline{\mathbf{a}}}_{\mathbf{p}}$$

In Appendix B, we illustrate the method by examining two simple cases for $N=1$ and $N=2$, which correspond to an exponentially decaying memory kernel and to a damped oscillatory kernel, respectively, and we explicitly write the associated $\mathbf{A}$ and $\mathbf{B}$ matrices. Generally, the kernels that can build via the extended variable formalism include (although not exclusively) functions of the form $K\left(\tau \right)=\mathrm{Re}\left[c_k{e}^{-{\alpha }_k\left|\tau \right|+i{\omega }_{k}t}\right]$ with arbitrary parameters $c_k$, ${\alpha }_k$ and ${\omega }_k$. Appendix A of ref. [34] provides a detailed mathematical formalism.

The classical FDT (8) is enforced when

$${\mathbf{B}}_p{\mathbf{B}}_p^T=m{k}_{B}T({\mathbf{A}}_p+{\mathbf{A}}_p^T),$$

(17)

where ${\mathbf{B}}_p$ is the $(N+1)\times (N+1)$ matrix that comprises $b_{pp}$, ${\mathbf{b}}_p$, ${\overline{\mathbf{b}}}_p$ and $\mathbf{B}$ (and similarly for ${\mathbf{A}}_p$). However, the noise matrix ${\mathbf{B}}_p$ can also be chosen not to fulfill the classical FDT, so that the GLE can be arbitrarily tuned in order to thermalize the different frequencies of the system with well-chosen effective temperatures.

While in [13], the focus was on the thermalization of classical degrees of freedom at different temperatures, in the following paper [14] the same authors applied the colored baths strategy in order to include quantum corrections to the classical dynamics of ions in systems that are weakly anharmonic, such as diamond and ordinary ice $Ih$. In order to design this quantum thermostat (QT), they considered the particular case of the harmonic oscillator $V\left(x\right)=\frac{1}{2}m{\omega }^2{x}^2$ for which the stochastic equations of motion (16) can be solved analytically and the exact kinetic and potential energies of the physical system explicitly computed. The parameters of the GLE are then tuned in order for the average position and momentum distributions to be as close as possible to that of a quantum harmonic oscillator that is Gaussian, as in the classical case, but with a width that is a function of the oscillation period $\omega$:

$$\frac{1}{m}\langle p^2\rangle =m{\omega }^2\langle x^2\rangle =\frac{\hslash \omega }{2}coth\left(\frac{\hslash \omega }{2k_{B}T}\right)=\theta (\omega ,T)$$

(18)

By properly adjusting $K\left(\tau \right)$ and the noise correlation function, one can tune the average fluctuations in position and in momentum so as to reproduce the quantum–mechanical behavior in Equation (18) over a wide range of vibrational frequencies. In practice, this is done via a fitting procedure of the ${\mathbf{A}}_{\mathbf{p}}$ and ${\mathbf{B}}_{\mathbf{p}}$ matrices that have been introduced in Equation (16). This method is general and versatile and offers a wide choice for the friction and random force kernels; nonetheless, it requires a fine tuning of the characteristic matrices, which might depend on the system under consideration. The QT method has been made available in the open-source i-PI suite [35], which has enabled its wider use [36,37,38].

4. Critical Analysis of Quantum Baths: Model Systems

In Section 3.1 and Section 3.2, the generalized baths (QTB and QT) are introduced as tools to impose the quantum statistics to a system that generally follows classical laws of motion. However, the question of the trustworthiness of quantum baths methods should be addressed. In this section, we provide a systematic analysis of the performances of the QTB on 1D and 2D model systems for which exact solutions are available analytically or can be obtained for comparison from a direct numerical resolution of Schrödinger’s equation (another example with a 1D double well is given in Appendix C). Although we focus here on the QTB formalism, similar results could be obtained with the quantum thermostat of [14].

4.1. Harmonic Systems

A pedagogical introduction to the QTB has been given by Barrat and Rodney, with several useful remarks, in particular, concerning the application of the QTB to a harmonic system, $V\left(x\right)=\frac{1}{2}m{\omega }_0^2{x}^2$ [39]. In that case, as the force is linear with position, an explicit expression can be obtained for the Fourier transform of the position and the velocity and hence for their respective correlation spectra:

$$C_{xx}\left(\omega \right)=\frac{2\gamma \theta (\omega ,T)/m}{{({\omega }^2-{\omega }_0^2)}^2+{\omega }^2{\gamma }^2}\mathrm{and}{C}_{vv}\left(\omega \right)=\frac{2\gamma {\omega }^2\theta (\omega ,T)/m}{{({\omega }^2-{\omega }_0^2)}^2+{\omega }^2{\gamma }^2}$$

(19)

These spectra correspond to Lorentzian peaks centered at ${\omega }_0$, and with a width $\gamma$. They are good approximations to the quantum symmetrized correlation spectra $C_{xx}^s\left(\omega \right)$ and $C_{vv}^s\left(\omega \right)$ apart from a spectral broadening with respect to the Dirac $\delta —$function expected for the 1D quantum system, which is only recovered in the $\gamma \to 0$ limit. The integrals of these two spectra provide, respectively, the average potential and kinetic energies as:

$$\begin{array}{ccc}\hfill \langle E_{pot}\rangle & =& \int \frac{d\omega }{2\pi }\frac{\gamma {\omega }_0^2}{{({\omega }^2-{\omega }_0^2)}^2+{\omega }^2{\gamma }^2}\theta (\omega ,T)\hfill \end{array}$$

(20)

$$\begin{array}{ccc}\hfill \langle E_{kin}\rangle & =& \int \frac{d\omega }{2\pi }\frac{\gamma {\omega }^2}{{({\omega }^2-{\omega }_0^2)}^2+{\omega }^2{\gamma }^2}\theta (\omega ,T)\hfill \end{array}$$

(21)

As the zero-point energy term in $\theta (\omega ,T)$ is linear for large $\omega$, the integral in (21) diverges logarithmically. Barrat and Rodney therefore suggested the introduction of a cutoff frequency ${\omega }_{max}$ in the random force spectrum, that should be larger than the highest frequency of the physical system. (Note that a finite integration time step $\delta t$ implicitly imposes a maximal frequency $\delta t^{-1}/2$. However, the use of a noise power spectrum without an explicit high-frequency cutoff ${\omega }_{max}$ poses various numerical problems). In Appendix A, the practical noise generation as well as the choice of the algorithms to integrate the Langevin equation is briefly discussed. With this precaution, in the limit of small values of $\gamma$, the Lorentzian factors in Equations (20) and (21) tend to Dirac $\delta —$functions and both the kinetic and the potential energy equal $\theta ({\omega }_0,T)/2$ as expected for the quantum harmonic oscillator (for nonzero $\gamma$, the average energies are slightly different from $\theta ({\omega }_0,T)/2$, although this difference can be corrected using spectral deconvolution techniques [18,40]). The position and the momentum probability distributions obtained in QTB simulations of the harmonic oscillator are Gaussian, with widths fixed by Equations (20) and (21), therefore the method provides exact estimates (at least in the $\gamma \to 0$ limit) for the quantum average of any static observable depending only on position or momentum, including zero-point energy effects. As mentioned above, it also yields a good approximation to some dynamical properties such as position and velocity autocorrelation spectra. This distinction is important in the context of nuclear quantum effects simulations as path-integral techniques can provide exact references for static properties (though at an elevated computational cost, see Section 6), but only approximate methods are available to simulate the quantum dynamical properties of complex systems with large numbers of atoms.

The paragraphs below examine the effect of anharmonicity on the accuracy of QTB simulations. However, even at the harmonic level, a remark should be made about static observables with a joint position and momentum dependence such as the total energy fluctuations:

$$\begin{array}{ccc}\hfill \Delta E^2=\langle E^2\rangle -{\langle E\rangle }^2=\langle {(E_{kin}+E_{pot})}^2\rangle -{\langle E_{kin}+E_{pot}\rangle }^2& =& {\left[\frac{\hslash {\omega }_0}{2}sinh{\left(\frac{\hslash {\omega }_0}{2k_{B}T}\right)}^{-1}\right]}^2\hfill \\ & =& k_B{T}^2\frac{d\theta ({\omega }_0,T)}{dT}\hfill \end{array}$$

(22)

for the quantum harmonic oscillator at temperature T. In contrast, in QTB simulations, the energy fluctuations are equal to $\theta {({\omega }_0,T)}^2$, as the method is equivalent in the harmonic case to a classical Langevin dynamics at the effective temperature $\theta ({\omega }_0,T)/k_B$ (at least in the small $\gamma$ limit). Indeed, even though the QTB describes correctly the fluctuations of the kinetic and potential energies separately, it cannot capture the intrinsically quantum correlation between the two observables, which leads to a systematic overestimation of $\Delta E^2$, in particular at low temperatures $\frac{\hslash {\omega }_0}{2k_{B}T}\ll 1$. As a consequence, in QTB simulations, quantities such as the heat capacity should be evaluated by differentiation of the average energy with respect to temperature, as they cannot be directly related to energy fluctuations as it is usual in classical MD. The overestimation of the energy fluctuations may also have indirect consequences on the estimation of some dynamical observables in more complex systems (in particular barrier crossing rates).

4.2. Morse Potential in One Dimension

The anharmonic Morse potential (Figure 3) can be written as:

$$V\left(x\right)=V_0{\mathrm{e}}^{-\frac{x-x_0}{d}}\left({\mathrm{e}}^{-\frac{x-x_0}{d}}-2\right)$$

(23)

where $V_0$ is the depth of the well, $x_0$ the position of the minimum and d the ‘width’ of the well which in practice controls the anharmonicity thereof.

4.2.1. A ‘Mildly’ Anharmonic Example

Simulations were carried out with the following parameters: $V_0=3.0$ eV, $x_0=1.0$ Å, $d=0.4$ Å, $T=300$ K, $\gamma =0.1$ THz and cutoff frequency $f_{\max }=500$ THz. The mass of the particle is that of a proton. These parameters yield a harmonic frequency of 3185 cm${}^{-1}$, of the same order of magnitude as that of OH stretching modes in various systems. The zero-point energy over $V_0$ ratio is 0.065, so this example can be considered as mildly anharmonic, given that, at this low temperature, only the ground state is significantly populated.

The quantum analysis and the QTB both yield a mean postition $\langle x\rangle =1.02$ Å, while the classical simulation yields 1.00 Å. The quantum total energy (actually the zero-point energy at this temperature) turns out to be 194 meV while the QTB gives 197 meV (for the small friction coefficient $\gamma$ used in this case, the kinetic energy overestimation effect mentioned in Section 4.1 is minor and the dependence on the frequency cutoff is negligible) and the classical only 26 meV which, of course, does not take into account the large ZPE. The quantum distribution is shown in green in Figure 3 while the classical distribution (yellow) in much sharper. The distribution generated by the QTB (blue) is adequately broadened, however, its shape is slightly modified: the maximum is at the the classical position, while its tail reaches out further.

Figure 4 shows the corresponding spectra. The QTB spectrum is broadened and red-shifted with respect to the classical one as the anharmonic part of the potential is explored further, but red-shift of the actual quantum transition (between the ground state and first excited states) is even larger; therefore, the QTB captures the correct trend while slightly underestimating the quantum effect. The inset shows the overtone resonance at approximately twice the main peak frequency. Similar shifts are observed between the classical, QTB and exact quantum results as for the main resonance. The detailed perturbation analysis performed in [16] shows that an analogy can be drawn between the QTB and linearized semiclassical initial value representation (LSC-IVR) methods, that combine quantum phase space sampling with classical molecular dynamics for the evaluation of time correlation functions [41,42]. Indeed, in QTB simulations, the short-term dynamics is essentially classical, whereas the coupling with the bath introduces quantum statistics in the phase space sampling. It can be shown that the overtone intensity is underestimated by approximately a factor of two due to the absence of position–momentum correlation in the QTB phase space sampling [43,44]. This error on the overtone intensity can be detected (and potentially corrected) using the first-kind fluctuation–dissipation theorem criterion defined in Section 7, and it remains limited compared to that of path integral methods (see Section 6), which yield a much lower overtone intensity, similar to that of classical simulations [16].

4.2.2. A ‘Strongly’ Anharmonic Example

Now the following parameters are used: $V_0=1.7$ eV, $x_0=1.0$ Å, $d=0.2$ Å and $T=300$ K. This yields a harmonic frequency at 3185 cm${}^{-1}$ and ground state energy over depth ratio of 0.75, that is significantly more anharmonic. Another qualitative criterion for anharmonicity can be found in how the exact quantum distribution obtained from solving Schrödinger’s equation extends beyond the potential inflexion point: beyond this point, the restoring force decreases instead of increasing, and therefore contributes strongly to anharmonic behavior. In this example, the exact distribution (green line in Figure 5) does significantly extend beyond that point ($x\approx 1.14$ Å), while in the previous (Figure 3) it does not.

Figure 5 shows the results. The features already noticed for the mildly anharmonic case are enhanced: the maximum still remains at the classical position, while the tail extends even further towards larger distances, requiring the introduction of $x_{\max }=2.6$ Å to prevent the particle from escaping the well: in this strongly anharmonic situation, the calculated QTB mean position therefore will depend on how $x_{\max }$ is chosen and ends up being unreliable (${\langle x\rangle }_{\mathrm{QTB}}=1.27$ Å, versus the exact 1.04 Å). The total energy is also over-estimated (exact quantum: 184.4 meV, QTB: 250.4 meV, classical: 26.3 meV), mainly because the increased probability at large distances x leads to an over-estimation of the average potential energy. The spectra show qualitatively the same features as above, but are much broader as expected from the increased anharmonicity.

In these simple 1D examples, the QTB can be considered to a certain extent as a classical dynamics simulation at a higher effective temperature: $T_{\mathrm{eff}}\simeq 1800$ K in the above example. Indeed, Figure 6 shows that the classical distribution, obtained at $T=T_{\mathrm{eff}}$, is quite similar to that of the QTB at $T=300$ K. This should not mean that it suffices to estimate the effective temperature and carry out a classical simulation to describe NQEs correctly! In realistic systems with many degrees of freedom, each vibration mode is characterized by a different frequency and therefore a different effective temperature, so that performing classical simulations at an averaged $T_{\mathrm{eff}}$ would yield erroneous results [45]. The main advantage of the QTB and QT approaches is that the effective temperature imposed to each mode naturally corresponds to their frequency by construction of the generalized bath (including approximate anharmonic shifts, as shown above), without requiring prior knowledge on the system nor harmonic frequency calculations.

4.3. Several Degrees of Freedom and Zero-Point Energy Leakage

While in one-dimensional systems, the Quantum Thermal Bath introduces a significant correction to the classical Langevin simulations by taking into account the quantum zero-point energy, the issue of energy transfer between modes due to anharmonicity is irrelevant for lack of several modes to transfer between. However, most simulations deal with systems that involve a large number of degrees of freedom and such transfers are quite likely to occur. The unwanted consequence thereof is Zero-Point Energy Leakage (ZPEL), first illustrated with a simple two-dimensional model and then described from a more general viewpoint: a diagnosis and a cure are presented in Section 7.

4.3.1. A 2D Model for an O–H Vibration

We extend the above analysis to the following 2D model:

$$V(x,y)=V_0{\mathrm{e}}^{-\frac{r(x,y)-r_0}{d}}\left({\mathrm{e}}^{-\frac{r(x,y)-r_0}{d}}-2\right)+\frac{1}{2}k\theta {(x,y)}^2$$

(24)

$r(x,y)$ and $\theta (x,y)$ represent the polar coordinates associated with the Cartesian coordinate $(x,y)$. This model represents an O–H bond in the direction x (associated to the ‘mildly’ anharmonic Morse potential of Section 4.2 for the stretching vibration) which allows bending vibrations (with a harmonic potential for the angle $\theta$ with force constant $k=9.48{\mathrm{eV}.\mathrm{rad}}^{-2}$). The mass for both degrees of freedom is set to that of a proton. Figure 7 shows the quantum, classical and QTB probability distributions at 300 K. As for the 1D models, the QTB captures very well the strong broadening caused by zero-point energy effects, although the details of the shape of the QTB distribution display slight discrepancies with respect to the quantum reference.

4.3.2. Zero-Point Energy Leakage

Table 1. Average energies for the 2D model of O–H bond at 300 K (energy values in meV). Quantum results are obtained from a PIMD simulation with 96 beads.

	$\langle {\mathit{E}}_{\mathbf{kin}}\rangle$	$\langle {\mathit{E}}_{\mathbf{stretch}}\rangle$	$\langle {\mathit{E}}_{\mathbf{bend}}\rangle$	$\langle {\mathit{E}}_{\mathbf{tot}}\rangle$
Classical	26	13	13	52
QTB ($\gamma =1$ THz)	143	87	64	294
QTB ($\gamma =10$ THz)	146	90	60	296
QTB ($\gamma =20$ THz)	149	93	56	298
Quantum	143	99	49	291

shows different energy contributions calculated exactly, with the QTB and with classical MD. We easily notice that the QTB provides a massive energy correction to the classical MD and that the total energy is in good agreement with the exact one. The total kinetic energy in QTB is slightly overestimated due to spectral broadening effects for finite $\gamma$ coefficients (as presented in Section 3.1). The most notable energy errors are found in the stretching and bending contributions to the potential energy: the former is underestimated while the latter is overestimated. This discrepancy is typical of the zero-point energy leakage (ZPEL). Note that the impact of the ZPEL is reduced when the coupling coefficient $\gamma$ is increased.

Zero-point energy leakage is a well-known issue, initially described in the context of LSC-IVR simulations, in which the system is first prepared with the appropriate quantum energy distribution, after which short classical simulations are carried out to assess the dynamics. As the dynamics are indeed classical, the equipartition theorem will take over and eventually destroy the quantum distribution [46,47,48]. This is clearly connected with anharmonicity that allows vibrational modes to exchange energy, resulting, in practice, in a flow of energy from high-frequency modes to low-frequency modes. As the effective temperature can be very high, the outcome might be dramatic, as the melting down of the system in solid-state simulations. In the case of LSC-IVR simulations, the problem can be addressed by carrying out short enough classical dynamics, so that the leakage does not significantly alter the energy distribution and thus the dynamical results.

It is not so with the quantum baths methods [34,49] as the quantum energy distribution is generated on the fly through the colored Langevin thermostat. The result is therefore a compromise between how much energy is being pumped into the system by the colored noise and how fast this energy is flowing to low frequency vibrational modes finally to be dissipated by the friction forces. A systematic endeavour [50] to measure ZPEL in QTB simulations as a function of anharmonicity and the coupling coefficient $\gamma$ in simple model systems such as non-linearly coupled harmonic oscillators showed that ZPEL indeed increases with anharmonicity and that increasing the coupling coefficient $\gamma$ significantly reduces the effects of ZPEL. Similar observations were made regarding the QT in the early papers on this method [14,34] and strongly coupled generalized baths were designed in order to mitigate the effect of ZPEL, an approach comparable to the increase of the coupling coefficient $\gamma$ of the QTB. However, the spectral broadening induced by large coupling coefficients can hinder the use of quantum baths for the computation of vibration spectra. Even more importantly, large system–bath couplings do not completely suppress ZPEL, which can remain significant in strongly anharmonic systems such as liquid water where it has massive consequences, as examined in Section 7.3. Some authors also proposed to modify the colored noise memory kernel in the QTB method in order to compensate for the ZPEL and restore the correct energy distribution between the different modes [49,51]. However, without an appropriate criterion to detect and quantify the ZPEL in general cases, these attempts remained somewhat ad hoc and could not be generalized. Such a criterion was finally derived in [17], and is the basis for the adaptive QTB method presented in Section 7.

5. Applications of Quantum Baths to Realistic Systems

Quantum baths have been used to introduce the quantum statistics in the simulation of several systems, and yield valuable results. In the following, we discuss some selected applications of these methods to simulate the properties of condensed matter systems, by splitting the discussion between structural properties and time-dependent observables. These examples illustrate the usefulness of quantum bath approaches and highlight their most serious pathology, i.e., zero-point energy leakage, which motivates the detailed analysis of ZPEL and the recent work to mitigate its effects presented in Section 7. The studies reviewed in this section include both QTB and QT simulations, though fewer in number for the latter method. We note that the GLE framework developed in the context of the QT also found a wide application in its extension to path integral molecular dynamics (briefly presented in Section 6).

5.1. Structural and Thermodynamic Properties

An important quantity for crystals is the Debye temperature ${\mathsf{\Theta }}_D$, below which the specific heat and other thermodynamic properties diverge from their classical behavior [52]. In particular, crystals with ${\mathsf{\Theta }}_D>300$ K can show significant nuclear quantum effects at room temperature.

5.1.1. MgO

Dammak and coworkers [12] showed that the QTB reproduces the experimental trends of the lattice parameter and the heat capacity of the MgO crystal as a function of temperature. While recovering the classical limit for temperatures close to or above the Debye temperature (${\mathsf{\Theta }}_D\simeq 940$ K in MgO), the specific heat follows the expected quantum behavior at low temperatures ($C_V\to 0$ as $T\to 0$), in excellent agreement with experimental data.

5.1.2. Diamond

Diamond has a particularly high Debye temperature ${\mathsf{\Theta }}_D\simeq 2000$ K. Ceriotti and coworkers [14] applied the quantum thermostat (QT) to model diamond crystals, where interatomic interactions were described with a semi-empirical potential, and compared their results to existing path integral (PI) simulations. The internal energy $E_{int}\left(T\right)$ and lattice parameter $a\left(T\right)$ as provided by the QT follow the PI trends with temperature down to $T\sim 0.06{\mathsf{\Theta }}_D$. Below this temperature, the quantum thermostat failed to counterbalance the strong phonon–phonon coupling, which results in zero-point energy leakage from high to low frequencies and a too-short phonon lifetime. Therefore, $E_{int}\left(T\right)$ and $a\left(T\right)$ sensitively diverge from the PI behavior at low temperatures.

Isotope Effects: LiH vs. LiD

Simulations that rely on classical statistical mechanics cannot account for isotope effects straightforwardly, because statistical averages are independent of the nuclear mass. A marked isotope effect takes place in LiH/LiD crystals; the Debye temperatures, as obtained by fitting the thermal conductivity of LiH and LiD crystals, are 851 K and 638 K, respectively (see [53] and references therein). One of the first QTB simulations in conjunction with the first-principle description of interatomic forces via the density functional theory [53] explained why the lattice parameter of LiH is significantly larger than that of LiD from very low temperatures up to 600 K, as experimentally observed. The volume expansion coefficient is also different in this T range between the two isotopes. The elastic properties of LiH and LiD with increasing pressure are markedly distinct, which reflects in two different equations of state for the two isotopes as measured by X-ray diffraction experiments in the 0–94 GPa pressure range [54]. Remarkably, the quasi-harmonic approximation fails in reproducing the isotope shift in pressure $\Delta P$, which is defined as the difference in pressure between LiH and LiD at fixed volume (see Figure 8). The QTB simulations follow the experimental trend in $\Delta P$, while the quasi-harmonic approximation deviates from it, especially when P increases. This has been interpreted as a consequence of the importance of anharmonic contributions to the thermal properties of those crystals.

5.1.3. Energetic and Structural Properties of Atomic and Molecular Clusters

Hernández-Roja and coworkers applied the QTB to study rare-gas clusters, modelled via Lennard–Jones interatomic potentials [49]. They found that the method generally reproduces energetic properties relatively well, reinforcing earlier results on cationic neon clusters [55]. However, close to the liquid-like to solid-like (l-s) transition, the QTB performed unevenly against PI-based simulations: the l-s transition in Ar clusters is shifted to low temperatures and the Ne clusters show liquid-like behavior at all temperatures, even when PI simulations would predict a solid-like behavior. As pointed out by the authors, this failure is due to ZPEL. They also encountered similar difficulties in the simulation of water octamers with the QTB, where the strong leakage of zero-point energy stemming from high-frequency intra-molecular modes completely distorts the cluster structure and causes spurious melting, even down to low temperatures [49].

5.1.4. Proton Momentum Distribution

In the original quantum thermostat paper, the authors applied the newly introduced method to compute the proton momentum distribution of ordinary ice Ih [14]. NQEs strongly broaden this distribution with respect classical simulations and shift its maximum to larger momentum values. The QT was able to capture these effects and yielded close agreement with experimental measurements. In a later study, the same method was applied to lithium imide Li₂NH, using density functional theory to model atomic interactions, and the results were compared to inelastic neutron scattering measurements [56]. In this system, the N-H bond is strongly anharmonic, as evidenced by the divergence between molecular dynamics and quasi-harmonic results for the radial distribution of the N-H groups. Still, the QT method captures NQEs efficiently in this context and allows the recovery of a fine agreement with the experimental proton momentum distribution.

5.2. Structural Properties: Hydrogen Bonds

Hydrogen bond symmetrization under pressure provides a stringent test on the capability of quantum baths to reproduce nuclear quantum effects at extreme conditions that are relevant for planetary physics.

5.2.1. High-Pressure VII-X Transition in Ice

In the crystal phase VII of ice under high pressure, hydrogen atoms are connected to their nearest neighbour oxygen via a covalent bond and to their next-nearest neighbour via a hydrogen bond: O−H ⋯ O. As pressure increases, the oxygen atoms become closer and the proton tends to delocalize over two equivalent positions (resonating O−H⋯O and O⋯ H−O configurations) [57]. Beyond a critical pressure $P_c$, the asymmetry vanishes: there is no more a (longer) hydrogen bond and a (shorter) ionic–covalent bond; the proton forms two equivalent bonds to the neighboring O atoms. This proton-centered phase (ice X) is cubic with $Pn\overline{3}m$ symmetry (space group 224). The order parameter x for this transition can be taken as the difference of the proton distances from the two neighboring O atoms [57]. By first-principle constrained minimization at fixed x and excluding all nuclear quantum effects, an expression can be obtained fitting the proton effective potential energy [1]:

$$V(x;P)=a{x}^4+b(P-P_0)x^2+\frac{b^2{(P-P_0)}^2}{4a}\mathrm{for}P\le P_0$$

(25)

with $a=7.2$ eV/Å${}^4$, $b=0.04$ eV/(GPa Å${}^2$), and $P_0=100$ GPa. Under increasing pressure P, the O-O distance shrinks; the potential wells come closer and the barrier lowers (see Figure 9). In the classical framework, the proton centering occurs when the barrier disappears, so the critical pressure is $P_c=P_0$; in the quantum framework, the proton is delocalized, so that proton centering occurs at much lower pressure, when the zero-point energy equals the barrier height so that $P_c\simeq 70\mathrm{GPa}<P_0$. Experimentally, the VII to X phase transition in ice at room temperature under pressure takes place at approximately 65 GPa, while classical simulations, that do not take NQEs into account, mostly predict a transition pressure between 90 and 100 GPa. In QTB simulations at different pressures [1], the onset of hydrogen-bond symmetrization occurs at $P\simeq$ 65–70 GPa, in agreement with previous results obtained via a path-integral based method [57].

5.2.2. AlOOH

Under increasing pressure, the $\delta$ phase of AlOOH crystal undergoes a phase transition from a $P{2}_{1}nm$ structure with asymmetric and disordered O−H⋯O bonds to a stiffer $Pnnm$ phase with symmetric hydrogen bonds, which should be stable within the pressure and temperature ranges typical for the Earth’s mantle. In QTB simulations at pressures as low as 10 GPa, the mean proton position equals the midpoint of the O-O distance [58], while only at much larger pressures (around 30 GPa) does the $Pnnm$ phase with symmetric hydrogen bonds become stable according to the classical viewpoint [59]. However, the proton centering occurs when the effective potential that is seen by the proton is still asymmetric. The lack of the full inversion symmetry impacts the proton distribution, which is skewed as provided by the QTB even beyond the critical pressure [58]. At the transition the O-H stretching modes soften considerably and fade out around 10 GPa in the $P{2}_{1}nm$ structure, when thermal and nuclear quantum effects are taken into account in the simulations. Later X-ray diffraction experiments [60] also confirm the transition from a hydrogen-off-centered to a hydrogen-centered phase above 10 GPa.

Overall, the QTB provides a satisfying description of hydrogen bonds under increasing pressure, either where a center of inversion exists (as in ice X) or does not (as in $\delta$-AlOOH).

5.3. Dynamical Properties

Although initially devised to introduce the Bose–Einstein statistics for the evaluation of equilibrium properties, the quantum baths have been also used to evaluate dynamical properties. The following results show that the performances of quantum baths are rather system- and property-dependent, although in most cases the introduction of NQEs via the quantum baths improves the description of vibrational properties, even in very anharmonic systems. Although the coupling with the bath tends to distort vibrational spectra, particularly in the context of the QT method, this effects can be efficiently corrected a posteriori using spectral deconvolution techniques [18,40]. Moreover, quantum baths have been employed to study intrinsically nonequilibrium properties, such as thermal conductivity in crystals (see Section 5.3.3).

5.3.1. Phonon Spectra in Pure and Salty Ices under Pressure

As observed in Section 4, QTB simulations provide direct access to vibrational spectra, accounting consistently (though approximately) for anharmonic and quantum effects. This property has been exploited in different studies to confront QTB results with other, more established simulation methods [15,18], or with experimental data. In particular, the phonon spectrum of ice was extracted from QTB simulations at various pressures, computed from the Fourier transform of the velocity–velocity autocorrelation function [1].

These QTB results are shown in Figure 10, and compared with the infra-red absorption and Raman scattering data. The QTB reproduces the experimental data well in the whole pressure range, including the softening of the O-H stretching mode close to the VII-X transition. This softening is an essential feature of the transition that can be related to the response of the system: in the classical picture, the squared frequency of the soft phonon mode is proportional to the inverse susceptibility ${\chi }^{-1}$ and, in the simple one-dimensional double-well model, $\chi$ is maximum at the transition pressure $P_c$. This relates the O-H stretching mode softening with the enhanced fluctuations of the dipole moment that typically occurs at the transition. Interestingly, the O-H stretching mode softening occurs at larger pressures when a small quantity (∼1%) of LiCl or NaCl salt is incorporated into the ice crystal, an effect that is well captured by the QTB [61].

Contrary to previous beliefs, this effect is not due to steric hindrance, but to the a-symmetrization of the OHO triplet because of the long-range dipolar field generated by the dissociated Li${}^{+}$ Cl${}^{-}$ or Na${}^{+}$ Cl${}^{-}$ salts. Accordingly, the critical pressure for the VII-X transition in salty ice under pressure is upshifted and much closer to its classical value than in pure ice. This counterintuitive phenomenon (salty ice is “more classical” than pure ice) shows the subtleties of nuclear quantum effects in hydrogen-bonded systems.

5.3.2. Highly-Excited Molecules

Calvo and coworkers [62] studied the response of naphtalene to a strong laser pulse at varying frequencies by using classical molecular dynamics or accounting for nuclear quantum effects via path-integral molecular dynamics or the QTB. First, they observed that the spectrum of the transferred energy from the laser to the molecule is much distinct in the classical and the quantum frames; the path-integral and the quantum bath spectra both yield absorption at lower frequencies than the classical molecular dynamics. From the trajectories, and using a pre-computed rate constant, they predicted the emission of a hydrogen atom. Moreover, for the dissociation probability spectrum, QTB and path-integral yield consistent descriptions, which both differ from the harmonic picture.

5.3.3. From Equilibrium to Nonequilibrium

In this section, we first introduce the approach adopted by Dhar and Roy [63] to compute the thermal conductivity within a fully quantum formalism in the case of a one-dimensional harmonic chain of atoms. We then shortly describe how this method was generalized to anharmonic models by Wang [64], by introducing a quasiclassical approximation in which the thermal conductivity is evaluated through molecular dynamics simulations with a quantum heat bath (QHB). This approach has enabled the study of two realistic systems [65] and it presents striking formal similarities with the QTB (and QT) methods. We conclude this section by briefly reviewing a few studies that exploit the QTB formalism to evaluate similar nonequilibrium properties.

We briefly recall the connection between equilibrium properties and time-dependent propagators via atomic Green’s functions, which is the subject of a vast amount of literature [66]. Here, we follow the formalism adopted by Dhar and Roy for a harmonic system [63]. They considered a one-dimensional harmonic chain of N particles of mass M, with displacement $u_j$ ($j=1,\dots ,N$) with respect to their equilibrium position. The chain is linearly coupled at its left (L) and right (R) sides with two semi-infinite harmonic lattices that act as heat reservoirs at two different temperatures $T_L$ and $T_R$. Substituting the formal solution of the quantum propagation for these reservoirs, one can eliminate the bath degrees of freedom and obtain the quantum Langevin equation for the displacement vector: $u^C=(u_1,\dots ,u_N)$:

$${\ddot{u}}^C=-K{u}^C-{\int }_{t_0}^{t}ds{\mathsf{\Sigma }}^{+}(t,s)u^C\left(s\right)+{\xi }_L\left(t\right)+{\xi }_R\left(t\right)$$

(26)

where K is the $N\times N$ coupling matrix describing the harmonic forces within the chain, and ${\mathsf{\Sigma }}^{+}={\mathsf{\Sigma }}_L^{+}+{\mathsf{\Sigma }}_R^{+}$ is the self-energy matrix due to the combined interaction with the left and right reservoirs. Since the baths are assumed to be harmonic, ${\mathsf{\Sigma }}^{+}(t,s)$ can be expressed analytically in the frequency domain. The noise vectors ${\xi }_L\left(t\right)$ and ${\xi }_R\left(t\right)$ can then be related to the self-energies by assuming that the reservoirs are at thermal equilibrium and therefore obey Bose–Einstein statistics, which yields the following fluctuation–dissipation relation for the symmetrized correlation function in the frequency domain:

$$\frac{1}{2}\langle \tilde{\xi }\left(\omega \right){\tilde{\xi }}^T\left({\omega }^{\prime }\right)+\tilde{\xi }\left({\omega }^{\prime }\right){\tilde{\xi }}^T\left(\omega \right)\rangle =\delta (\omega +{\omega }^{\prime })\frac{\hslash \Gamma \left(\omega \right)}{2\pi }\coth \left(\frac{\hslash \omega }{2k_{B}T}\right)$$

(27)

where $\Gamma \left(\omega \right)=\mathrm{Im}\left[{\mathsf{\Sigma }}^{+}\left(\omega \right)\right]$ and $\tilde{\xi }\left(\omega \right)$ refer either to the left or the right lead, with the appropriate temperature $T=T_L,T_R$. Here the bath self-energy provides the connection between the equilibrium correlation function of the bath and the time propagator of the system that is needed to evaluate transport and nonequilibrium properties. The quantum Langevin Equation (26) is formally very similar to the GLE (6), that was also derived under the assumption of bilinear coupling with an harmonic bath and extended by Dhar and Roy to a nonequilibrium case with two baths at different temperatures.

The atomistic Green’s function for the chain can then be expressed as:

$$G^{\pm }\left(\omega \right)={\left[-M{\omega }^2+K+{\mathsf{\Sigma }}_L^{\pm }\left(\omega \right)+{\mathsf{\Sigma }}_R^{\pm }\left(\omega \right)\right]}^{-1}$$

(28)

In this expression, the harmonic frequencies of the isolated chain are renormalized by the self-energies of the interaction with the right and left bath leads. Provided ${\mathsf{\Sigma }}_L^{\pm }\left(\omega \right)$ and ${\mathsf{\Sigma }}_R^{\pm }\left(\omega \right)$ are known (the advanced and retarded self-energies ${\mathsf{\Sigma }}^{\pm }$ are related to the interaction at the interfaces between the system sites and the bath sites; therefore, the self-energies depend on the bath mode distribution and on their geometry—for a thorough account of the atomistic Green’s function method, see, for instance [67]), an expression is then derived for the thermal conductivity, by relating the energy flux though the chain and the temperature difference $T_L-T_R$ between the two reservoirs. In this fully quantum treatment, which is made possible by the fact that the system is harmonic, the thermal conductivity tends to zero at low temperatures, contrary to the classical description of heat transport.

5.3.4. Thermal Conductivity of Insulating Crystals

Wang aimed at simulating heat transport in realistic nanostructures from the ballistic regime at low temperatures to the diffusive regime at high T, via molecular dynamics [64]. He noted that the kinetic theory for phonons provides the thermal conductivity as $\kappa =\frac{1}{3}cvl$, where c is the heat capacity, v the sound velocity and l the mean free path. Therefore, classical molecular dynamics cannot provide a reliable estimation of heat transport in the low-temperature regime as c goes to a non-zero constant for $T\to 0$ in the classical framework, which clearly contradicts the experimental facts.

Wang then followed the Dhar and Roy’s approach [63], that he extended to anharmonic systems. In order to make the problem tractable in that case, Wang introduces a quasiclassical approximation similar to that of Schmid [32], and ended up with a generalized Langevin equation. It is formally identical to (26), apart from the fact that the harmonic force term $-K{u}^C$ is replaced by a general nonlinear force $F\left(u^C\right)$. However, the interpretation of this equation is completely modified in the quasiclassical approximation as, instead of quantum position operators, the vector $u^C\left(t\right)$ now describes the classical-like stochastic trajectories of the atoms that can be obtained by integrating the GLE using molecular dynamics techniques. Wang named such a method quantum molecular dynamics, but in the context of the present review, we refer to this approach as quantum heat bath (QHB): he showed that the QHB allows the recovery of the exact quantum thermal conductivity in the case of harmonic chains, and provides consistent results for both the high- and low-temperature regimes for anharmonic one-dimensional chains [64].

Later on, Wang and coworkers [65] applied the QHB to analyze thermal transport in a two-dimensional realistic model of graphene nanoribbon. They show that the QHB results are more consistent than other schemes that use a posteriori quantum corrections to otherwise classical dynamics. For instance, some authors run classical MD simulations at a “quantum-equivalent classical temperature” $T_{MD}=\int d\omega D\left(\omega \right)\theta (\omega ,T)/k_B$, where $D\left(\omega \right)$ is the phonon density of states [68]. Such a posteriori correction schemes are nevertheless not fully consistent, even for the one-dimensional harmonic chain, as they differ from the exact result by a numerical factor [65]. Interestingly, Wang and coworkers also provide a diagrammatic perturbative expansion of the QHB results for weakly anharmonic systems and analyze their discrepancies with respect to a fully quantum-mechanical nonequilibrium Green functions approach. From this analysis, they argue that QHB is exact in leading order in the anharmonic perturbation for the phonon lifetime.

The thermal conductivity of solid argon, modelled via Lennard–Jones interatomic potentials, was also investigated using nonequilibrium QTB molecular dynamics by Bedoya-Martinez et al. [51]. In these calculations, the simulation cell is divided into slabs, two of which are thermalized at different temperatures. The ratio between the imposed temperature gradient and the resulting heat flux then yields the thermal conductivity. As was observed for other materials [12,69], the QTB greatly improves the estimation of the heat capacity with respect to classical simulations. Nonetheless, these authors observed that, at low temperatures (where NQEs are expected to be most relevant), the experimental trends for the thermal conductivity were much better predicted when a classical Langevin equation was used as a thermostat than with the QTB. This result, which contrasts the outcomes by Wang and coworkers [64,65], was related by Bedoya-Martinez et al. to the underestimation of the phonon lifetimes. After the analysis of the simulations, they attributed the too-short phonon lifetime to the presence of a strong ZPE leakage, which tends to equalize the effective temperature for all phonon modes irrespective of their frequency, and partly to a more fundamental limitation of the QTB method, that would be intrinsically unable to capture this quantity accurately (see also our analysis in Section 4).

5.3.5. Shock Compression

Qi and Reed [69] proposed atomistic simulations of shock-compressed materials that incorporates quantum nuclear effects on the fly. The authors combined the MultiScale Shock Technique to reproduce the kinetics of liquid CH₄, described through semiempirical interatomic potentials, with the QTB as a bath. The heat capacity as well as the temperature and pressure versus density curves are significantly improved by the introduction of NQEs in comparison with experimental data.

6. Path Integrals and Colored Bath

Quantum baths are also used to reduce the numerical cost of the reference method to include NQEs in simulations of realistic systems: path-integral molecular dynamics (PIMD). It can be shown, at least for time-independent properties, that PIMD systematically converges to the exact quantum result irrespective of the form of the interaction potential. This convergence can, however, come at a considerable numerical cost that sometimes hinders ambitious applications. The successful use of generalized baths to facilitate the convergence of this approach is therefore an additional motivation for their interest, regardless of the theoretical limitations discussed in previous sections. In this chapter, a very brief introduction to PIMD is provided, followed by a description of two approaches based on generalised baths that effectively improve the efficiency of path integral calculations.

6.1. Path Integral Molecular Dynamics

PIMD implements numerically a discrete version of Richard Feynman’s path-integral formalism applied in imaginary time to the thermal density matrix. In practice, there is a correspondence between the quantum equilibrium distribution and that of a classical equivalent system in which each particle is replaced by a set of P replicas or “beads”, subject to the following potential:

$$V_P(x_1,\cdots ,x_P)=\frac{mP}{2{\hslash }^2{\beta }^2}\sum _{\ell =1}^P{(x_{\ell }-x_{\ell +1})}^2+\frac{1}{P}\sum _{\ell =1}^{P}V\left(x_{\ell }\right)$$

(29)

The first term in the right-hand side is the bead–bead interaction that arises from the kinetic energy operator, while the second term V is the physical potential to which the original particle is submitted, e.g., interactions with other particles [70]. Standard molecular dynamics methods can then be used to obtain trajectories and sample a probability distribution of the classical equivalent system, which effectively includes NQEs. The main advantage of this method is that it converges towards the exact quantum distribution as the number of beads P is increased to infinity. Its computational cost is, however, an important drawback, especially at low temperature where the number of beads needed to reach convergence becomes very large. This issue can be mitigated by combining path integrals with quantum baths (see Section 6.2 and Section 6.3 below). Another more fundamental limitation concerns dynamical observables such as vibrational spectra: PIMD is based on an imaginary-time path integral mapping of the quantum equilibrium density to that of the classical equivalent system described by Equation (29), but this mapping does not capture the quantum dynamics (even in the large P limit). PIMD is nontheless at the basis of different approximations to the quantum dynamics among which, in particular, centroid MD and ring-polymer MD as well as other more computationally demanding methods [10,16,30].

6.2. Combining Path Integrals with Langevin Equation

In 2011, Ceriotti and coauthors proposed to use generalized Langevin thermostats in the context of PIMD simulations [71] (this combination was first denoted PI+GLE). In a standard PIMD simulation, the equivalent potential $V_P$ of (29) is sampled using molecular dynamics with a classical thermostat (e.g., a white noise Langevin thermostat) and the sampled distribution tends to the exact quantum one for large P. In PI+GLE, on the contrary, the equivalent system is attached to a GLE thermostat designed in such a way that, for a harmonic oscillator, the sampled distribution equals the quantum one for any number of beads P. For a single bead, this strategy reduces to the quantum thermostat approach described in Section 3.2, while for large P, it tends to the standard PIMD. By construction, the method is exact for harmonic systems, and, for general anharmonic systems, the number of beads required to converge is significantly reduced with respect to standard PIMD. For example, for liquid water at ambient temperature, 6 beads yield similar accuracy as 32 in standard PIMD [71]. A further refinement of the method was then introduced under the name PIGLET [72], in which the GLE thermostat is applied in the ring-polymer normal mode basis in order to speed up the convergence of the quantum kinetic energy estimator. The PIGLET approach considerably reduces the computational burden of PIMD while retaining the systematic convergence with increasing number of beads. This approach enables accurate simulations at ultra-low temperatures feasible [73]. Although path-integral generalized Langevin equation methods have been limited to the study of static properties, Kapil and coworkers recently introduced a post-processing scheme that estimates the dynamical perturbation introduced by the generalized bath and recover time-correlation functions [74]. The PIGLET availability in the i-PI simulation platform [35] has enabled the wide use of this appealing compromise [75,76,77].

6.3. The Quantum Thermal Bath and Path-Integral Molecular Dynamics

Following the work by Ceriotti et al., Brieuc and coauthors proposed to combine PIMD with the QTB Formalism [78]. In this approach, denoted PIQTB, the path-integral equivalent system of Equation (29) is thermalized via a Langevin equation with a Markovian friction and a colored random force (as in the QTB). As for the PIGLET method, PIQTB tends to standard white noise Langevin PIMD for large number of beads while it reduces to QTB in the single-bead case. For intermediate numbers of beads, the memory kernel of the random force is chosen in such a way as to enforce a modified second kind FDT, in order to ensure sampling of the exact quantum distribution when applied to harmonic oscillators. For anharmonic systems, PIQTB offers similar convergence speed-up as PIGLET down to very low temperatures [79]. In particular, combined with other recent developments, PIQTB has enabled accurate simulations of chemical systems in superfluid helium nanodroplet down to ultra-low temperatures on the order of 1 K, whereas the computational cost of converging such calculations using standard PIMD techniques would be prohibitive [80].

7. Overcoming the Leakage via the Adaptive Quantum Thermal Bath (adQTB)

Although quantum baths are very attractive for their simplicity and efficiency, the ZPEL is a major limitation, which needs to be cured. In this section we show how the first-kind quantum FDT can be used as a quantitative diagnostic tool for the ZPEL and to design an adaptive method in which it is systematically compensated for.

7.1. Quantitative Assessment of the ZPEL

As mentioned in Section 4.3 and Section 5, ZPEL and its consequences were observed in several studies with either the QTB or the QT methods. The ZPEL can be mitigated by increasing the damping coefficient $\gamma$, or equivalently in the QT formalism, by designing more strongly coupled GLE thermostats [34]. However, strong system–bath couplings also affect the accuracy of the GLE methods (even for harmonic systems, the QTB and the QT are only rigorously exact in the weak coupling limit), and distort the vibrational spectra. Furthermore, even for large $\gamma$ values, the ZPEL is not completely suppressed and it can still have non-negligible effects. Therefore, a first step to design an effective cure is to provide a reliable diagnostic tool. Such a diagnosis is derived from the first kind fluctuation–dissipation theorem (see Section 2.1) that, in the context of the QTB, can be written as:

$$2\theta (\omega ,T)\mathrm{Re}\left[\chi \left(\omega \right)\right]=C_{vv}\left(\omega \right)$$

(30)

where $\theta (\omega ,T)$ is the quantum energy distribution, $\chi \left(\omega \right)$ is the complex admittance or susceptibility that characterizes the response of the system to a perturbation, and $C_{vv}\left(\omega \right)$ is the Fourier transform of the velocity–velocity correlation function. By using the quantum first-kind FDT (Equation (3)), the symmetrized autocorrelation spectrum $C_{vv}^s\left(\omega \right)$ (replaced here by $C_{vv}\left(\omega \right)$), can be evaluated from the Fourier transform of the QTB trajectories as it would be in a classical framework, consistently with the quasiclassical approximation at the root of the method.

As introduced in Section 2.1, the first-kind FDT is a very general relation that holds for a quantum system at thermal equilibrium and indicates that the vibration power density is distributed according to the quantum energy distribution $\theta (\omega ,T)$. The first-kind FDT is more fundamental than the second-kind FDT as used in the original formulation of the QTB, which relates the strength of the friction and random forces in the generalized Langevin equation. Indeed, the QTB would thermalize the vibrational modes with the appropriate average energy if they were non interacting, but in general anharmonic systems, the modes are coupled. This generates an energy flow from high to low frequencies, even though the system-bath coupling obeys the second-kind FDT. This flow creates an energy unbalance that can be detected through the violation of the first-kind FDT (30).

To exploit this relation, it was shown that, in a Langevin dynamics, the linear susceptibility can be estimated from the velocity–random-force correlation [17]:

$$\mathrm{Re}\left[\chi \left(\omega \right)\right]=\frac{\mathrm{Re}\left[C_{vR}\left(\omega \right)\right]}{C_{RR}\left(\omega \right)}=\frac{\mathrm{Re}\left[C_{vR}\left(\omega \right)\right]}{2m\gamma \theta (\omega ,T)}$$

(31)

The first-kind FDT can then be rewritten as (here we use one-dimensional notations; generally, the FDT should be satisfied for each degree of freedom and for each frequency independently):

$$\mathrm{Re}\left[C_{vR}\left(\omega \right)\right]=m\gamma C_{vv}\left(\omega \right)$$

(32)

where $C_{vR}\left(\omega \right)$ is the Fourier transform of the velocity–random-force correlation function. All terms in Equation (32) are now easily computed during a simulation run, and deviations from it are a clear indication that ZPEL occurs. Indeed, the ratio between the two sides of the equation can be used to define an effective frequency-dependent temperature $T_{\mathrm{eff}}\left(\omega \right)$, that, if the first kind FDT holds, should follow the quantum distribution of energy:

$$T_{\mathrm{eff}}\left(\omega \right)=\frac{\mathrm{Re}\left[C_{vR}\left(\omega \right)\right]}{m\gamma C_{vv}\left(\omega \right)}{k}_B^{-1}\theta (\omega ,T)=k_B^{-1}\theta (\omega ,T)\mathrm{iff}\mathrm{FDT}\mathrm{holds}$$

(33)

The effect of ZPEL on this effective temperature is represented in Figure 11 in the case of liquid water. Even with a relatively large value of $\gamma$, the effective temperature deviates from the first kind FDT reference (shown as a black line), which is a clear evidence of ZPEL. In particular, at low frequencies, the effective temperature is overestimated by more than 50 K for hydrogen atoms (around 30 K for oxygen). This is a consequence of the leakage of the large zero-point energy, which is significant in the high-frequency bending and stretching molecular modes towards low frequencies, and it can lead to massive errors as it strongly distorts the water’s molecular structure (see also Figure 12).

7.2. A Cure through Strict Enforcement of the First Kind FDT

To design a cure for ZPEL, an adaptive procedure was proposed in order to compensate the energy unbalance measured via the first-kind FDT [17]. To that end, the second-kind FDT is broken and two different frequency-dependent coupling coefficients are introduced: ${\gamma }_r\left(\omega \right)$, that replaces $\gamma$ in the random force power spectrum, Equation (14), and ${\gamma }_f\left(\omega \right)$ that characterizes the memory kernel of the friction force that is rendered non-Markovian (it is related to the memory kernel in Equation (6), by $\tilde{K}\left(\omega \right)=2m\mathrm{Re}[{\gamma }_f\left(\omega \right)]$). The first-kind FDT is then written as [17]:

$$\mathrm{Re}\left[C_{vR}\left(\omega \right)\right]=m{\gamma }_r\left(\omega \right)C_{vv}\left(\omega \right)$$

(34)

In the adaptive QTB (adQTB) method, the frequency-dependent coefficients ${\gamma }_r$ and ${\gamma }_f$ are adjusted on the fly during the simulation through a relaxation process that nullifies the difference between the two sides of Equation (34). This means that the random and the friction forces adapt by modifying the energy balance for each frequency to compensate for ZPEL and enforce the proper quantum energy distribution. In practice, only one of the two coefficients ${\gamma }_f$ or ${\gamma }_r$ has to be adapted. Both possibilities were tested in [17] and provide similar accuracy in compensating the ZPEL, although the variant where ${\gamma }_r$ is adapted while ${\gamma }_f$ is kept fixed is much simpler to apply, as it does not require the implementation of a non-Markovian friction force.

The adaptive procedure was successfully checked on several model systems and on Ne₁₃ clusters where it was able to stop the spurious melting of the solid clusters that ZPEL causes in QTB simulations [17,49]. Plots of ${\gamma }_r\left(\omega \right)$ on Figure 13 clearly show how the bath adapts, reducing the energy input from the bath at low frequencies, while increasing it at high frequencies.

7.3. Application to Hydrogen-Bonded Systems: Liquid Water

In a very recent study, liquid water was simulated using adQTB [18] in its adaptive random force version (i.e., adapting ${\gamma }_r\left(\omega \right)$ while keeping the friction Markovian with a constant coefficient $\gamma$). Different observables such as average energies and radial distributions were compared to PIMD simulations and shown to closely follow this reference technique (see Figure 12). Constant-pressure simulations were also performed to evaluate the density and the vaporization enthalpy as a function of temperature, in good agreement with the path integral calculations, at a computational cost comparable with that of classical simulations.

Water is a particularly challenging problem for the (ad)QTB method. Indeed, it was shown in the literature that in hydrogen-bonded systems, two opposite NQEs are in competition. On the one side, the stretching zero-point motion tends to strengthen hydrogen bonding, while on the other, the bending zero-point energy weakens it [82,83]. In liquid water, these two effects almost cancel each other, so that overall NQEs remain limited. Capturing this subtle balance is therefore essential to reproduce the properties of water accurately. In the standard QTB method, the molecular structure of water is strongly washed out by the ZPEL, as illustrated in Figure 12, but the adaptive procedure restores the correct distribution of energy, including zero-point contributions, and hence recovers radial distribution functions very close to the PIMD ones.

Interestingly, this study also confirms the potential interest of the adQTB for the inclusion of NQEs in the estimation of vibrational spectra: the infrared absorption spectrum of water was computed and compared to thermostated ring-polymer molecular dynamics (TRPMD), one of the most common approximations derived from the path-integral formalism for the calculation of time-dependent properties. The adQTB provides comparable accuracy to TRPMD, and it even captures better the high-frequency part of the infrared spectrum, where overtone and combination bands appear, whose amplitude the TRPMD method is unable to reproduce [16,84,85].

8. Quantumness from Baths: State of the Art and Open Issues

Quantum baths rely on a kind of blend of classical and quantum features: the system degrees of freedom follow a classical dynamics in which quantum effects are introduced via a generalized Langevin equation. Langevin dynamics turns out to be remarkably versatile and enables various choices for the bath characteristics and the strength of the coupling to the system. In quantum baths methods, it is used to reproduce zero-point energy (ZPE) effects by thermalizing the modes of the system at a well-chosen frequency-dependent effective temperature. These approaches are particularly suited for the simulation of the statistical properties of systems that are at the borderline between the classical and the quantum worlds. This is typically the case of light nuclei in molecules and condensed matter, such as H, He, or Li, or even heavier nuclei at low temperatures. In these cases, ZPEs can dramatically change the energetic landscape of the system, especially for highly anharmonic systems and shallow potential energy profiles; at low temperatures (much smaller than the ZPE) some fluctuations can be frozen, which impacts the computed response functions, such as the heat capacity; when energy barriers are present between several distinct configurations, the quantum description can greatly differ from the classical one, whenever the extension of the nuclear wavefunction is comparable to the width of the energy barrier. Isotope effects are another important example of NQEs, as purely classical averages are mass-independent. Therefore, the variation under isotope substitution of the thermodynamical properties (geometric isotope effect) or the reaction rates (kinetic isotope effect) are often clear experimental markers of the importance of quantum mechanics. In addition, within the realm of condensed matter physics, simulations involve a large number of degrees of freedom subject to non trivial potentials that are often expensive to calculate. When NQEs cannot be neglected—in many cases even at room temperature and pressure—methods to account for them must therefore be adapted to these specific constraints. Quantum baths attempt to use a classical formalism in which the energy distribution is altered to follow the Bose–Einstein distribution via a stochastic process included into a generalized Langevin equation. Apart from some technical subtleties that come along with the use of the generalized Langevin formalism, these methods do not involve increasing the computational load substantially, which make them highly desirable for the simulation of condensed matter systems, with thousands of atoms and more, that could be unfeasible otherwise.

The idea of enforcing the correct quantum energy distribution in realistic molecular dynamics simulations through a generalized Langevin equation is rather recent. As for many other approaches [86], this quantum–classical blend is clearly an approximation, which must be assessed.

An important shortcoming of these approaches, which is common also to some other mixed quantum–classical methods, is the zero-point energy leakage, that is, the trend to transfer a significant part of the excess ZPE from high to low frequencies. In this respect, the use of a generalized fluctuation–dissipation theorem (FDT) to monitor and regulate energy fluxes between the bath and the distinct degrees of freedom of the system has provided an elegant and easily implementable solution to this problem. The last-generation adaptive QTBs enforce the first-kind FDT to every frequency separately rather than the second-kind FDT, as used in previous quantum baths methods, and, by keeping the ZPE leakage below a certain threshold via an effective built-in feedback mechanism, is affirming as a quite safe mixed quantum–classical approximate scheme to appreciate NQEs in many systems. Moreover, this solution to the ZPE leakage made it possible to distinguish the intrinsic inaccuracies of the QTB method from those that are induced by the leakage, thus opening the way to other possible improvements of quantum baths. We note that, although the adaptive procedure has so far only been explored in the context of the QTB, in principle it could be extended to the QT formalism as well.

We stress that, at variance with their classical counterpart, the equilibrium distributions obtained from quantum baths cannot be derived analytically apart from a few cases, such as the harmonic oscillators. We illustrated through a few examples the performances of quantum baths (focusing more particularly on the QTB) both on equilibrium, time-independent, averages, and time-dependent correlation functions. Quite a wide experience has thus been gained on various systems, differing by the kind of interatomic bonding and the thermodynamic conditions. Rather generally, the QTB yields the maxima of the particle distributions at the minima of the potential energy surface, as do the classical ones. This is a minor shortcoming for mildly anharmonic cases, but might be a serious limitation when dealing with strongly anharmonic systems, characterized by double wells or when shallow potential energy surfaces are combined with large ZPE contributions; in such cases, a comparison with exact probability distributions as obtained via path-integral-based simulations could be helpful. However, in many realistic problems, and even close to phase transitions (for which the potential typically presents double well features), quantum baths much improve the classical description, capturing the correct trends and yielding quantitatively reliable results for a number of observables.

It turns out the quantum baths provide also encouraging results when dynamical processes are under scrutiny. They combine approximate sampling of the quantum thermal distribution with an essentially classical time propagation. The vibrational spectra, which can be easily obtained and analyzed as in the purely classical dynamics, can provide a relevant description of vibrational properties, including anharmonic properties, overtones and resonances, although a built-in criterion for attaining the convergence with respect to fully quantum simulations has not emerged yet. Many questions remain open in that respect, and quantum transport phenomena represent an interesting challenge at the time of writing.

Another recent use of quantum baths relies on their coupling to other methods. In this respect, the combination of quantum baths with the path-integral formalism often provides a good compromise: while path-integral-based methods converge towards exact equilibrium states when the number of beads is increased, they are costly. When the beads are thermalized or coupled to quantum baths, their number can be significantly reduced. This enables the simulations of molecules and condensed matter systems at very low temperatures, which are otherwise out of reach within pure path-integral molecular dynamics, by keeping the number of beads relatively small.

To conclude, quantum baths have known an increasing number of applications for the simulation of material properties since their conception. In this review, we analyzed the difficulties in their use as well as their shortcomings; at the same time, we prospected new solutions for making them more reliable and physically sound. We anticipate that other improvements will be built, so as to make generalized (quantum) baths an increasingly reliable simulation tool for nuclear quantum effects in material science.