About the computation of finite temperatureensemble averages of hybrid quantum-classicalsystems with Molecular Dynamics

Molecular or condensed matter systems are often well approximated by hybrid quantum-classical models: the electrons retain their quantum character, whereas the ions are considered to be classical particles. We discuss various alternative approaches for the computation of equilibrium (canonical) ensemble averages for observables of these hybrid quantum-classical systems through the use of molecular dynamics (MD), i.e. by performing dynamics in the presence of a thermostat and computing time averages over the trajectories. Often, in classical or ab initio MD, the temperature of the electrons is ignored and they are assumed to remain at the instantaneous ground state given by each ionic configuration during the evolution. Here, however, we discuss the general case that considers both classical and quantum subsystems at finite temperature canonical equilibrium. Inspired by a recent formal derivation for the canonical ensemble for quantum classical hybrids, we discuss previous approaches found in the literature, and provide some new formulas.


Introduction
Molecular Dynamics (MD) 1 is conventionally considered to be the theoretical description of molecular or condensed matter systems that assumes the nuclei to be classical particles.
Therefore, these move according to Newton's equations, in the presence of their mutual interaction and of a force that somehow approximates the electron influence. In its traditional formulation (classical MD), the forces are parameterised in some analytical expressions that have been carefully developed over the years, and the numerical problem amounts to the propagation of a purely classical Hamiltonian system with a predefined potential function.
The so-called ab initio or first-principles MD 2 substitutes those analytical force definitions by the on the fly calculation of the quantum electronic structure problem, that provides the forces on the ions due to the electrons in a more precise -yet more costly -manner.
Still, a first principles MD simulation also consists in the integration of a purely classical problem, even though one needs to use quantum mechanics to obtain the forces at each time step. In both classical and first principles MD, the electrons are usually assumed to remain in their ground state, adiabatically adapting to the ions as they move. Therefore, both the classical MD and the (Born-Oppenheimer ground-state) first-principles MD are not, strictly speaking, hybrid quantum-classical dynamics. These methods have been widely employed for both equilibrium and out-of-equilibrium problems.
In the equilibrium case, when studying a system in, e.g., the canonical ensemble, one is normally not interested in the particular trajectory followed by a microstate, but on the ensemble average values of a given property. The MD simulations are then used to compute the multi-dimensional integrals that define those averages, by substituting them with timeaverages over dynamical trajectories of the system -typically, coupled to a thermostat. 3 But in any case, despite the finite temperature, normally the electrons are assumed to be frozen in the ground state. This may be a very good approximation if the electronic excited state energies are far higher than the thermal energies at that temperature. Yet in many circumstances those excited states cannot be ignored.
Out of equilibrium, this may in fact happen very frequently. For example, when dealing with photo-chemistry, that naturally involves electronic excitations. This situation calls for a non-adiabatic extension of the previous MD concept, that allows for "live" electrons: the dynamics must be that of a true hybrid quantum-classical model, in which both classical and quantum particles evolve simultaneously through a set of coupled equations. Two prototypical examples of truly hybrid dynamics are Ehrenfest equations 4 and surface hopping 5 (in this latter case, the electronic motion is stochastic and consists of "jumps" between the adiabatic eigenstates).
In equilibrium, one may also need to lift the approximation of ground state electrons, if the temperature is high enough that the thermal population of the excited states is not negligible. This is the situation discussed in this article. In molecular physics this may happen rarely, but in condensed matter, the metallic or near metallic systems naturally call for a computation of ensemble averages that acknowledges the non-zero population of electronic excited states at non-zero temperature, even if low. In any case, this situation begs the questions: what are the canonical ensemble averages for observables of hybrid systems, and can one compute them using MD? The standard procedure used within classical or adiabatic first principles MD is no longer directly applicable if one is simultaneously propagating nuclei and electrons. The idea of assuming ergodicity and attaching a thermostat to the dynamics, a concept designed for purely classical systems, is dubious at best. This issue has been addressed by performing MD propagating only the classical nuclei, but somehow incorporating the electronic temperature in the definition of the forces, instead of deriving them by merely assuming the electrons to be in the ground state. One possible route, based on the use of density-functional theory, was theorised by Alavi et al. 6 Their method was based on the use of density-functional theory 7 (DFT) to solve for the electronic structure problem. In particular, on the finite temperature extension of DFT (FTDFT), 8,9 that substitutes the ground-state energy functional by a free energy functional. Then, to perform first principles MD at finite temperature, the forces used to propagate the classical ions are given by the gradient of this free energy functional. In essence, this is the same idea that underlies the approach given in Refs. 10,11, except that in this latter case the formulation is general, and not tied to DFT. The fact that the electronic free energy -considered at fixed nuclear configurations -can be viewed as an effective classical Hamiltonian from which the hybrid quantum-classical partition function can be computed was already found by Zwanzig. 12 DFT is in fact the most common electronic structure method for the purpose of performing ab initio MD. The inclusion of electronic temperature effects is therefore usually managed with some form of FTDFT. In practice, this procedure consists of using a Fermi-Dirac distribution for the population of the Kohn-Sham orbitals that constitute the fictitious auxiliary non-interacting system employed to substitute the true interacting many-electron problem. The resulting density is used to compute the ionic forces, in lieu of the groundstate density. The procedure should be completed with the use of temperature-dependent exchange-and-correlation functionals, but this is often ignored, as the development of these functionals has proved to be very difficult.
In this work, we discuss this and other possible routes to obtain the rigorous canonical ensemble averages through the use of thermosthatted MD. The goal is to establish a clear theoretical link between the definition of hybrid ensemble averages and the manners that one can use to compute them using some form of MD with a thermostat. The basic idea consists of generating an ensemble with some form of MD, even if the generated ensemble is wrong, and then using a reweighting formula to compute the right averages. From the analysis, it emerges that, in fact, various possibilities exist. The dynamics for the classical particles moving on the free energy surface will be shown to be a very particular case of this general class. The relative efficiency of the various options may depend on the particular system and choice of electronic structure method. The idea of performing wrong or ficititious dynamics, and then correcting with some reweighting procedure, has been used in the past in the field of MD mostly with the objective of accelerating rare events -see for example Refs. 13,14. In this work we extend the same idea to the simulation of hybrid quantum-classical systems.
In Section 2 we present the expressions for the ensemble averages that constitute the target of the current work. Section 3 discusses the possibility of computing these using a hybrid quantum-classical non-adiabatic MD such as Ehrenfest dynamics. Although the naïve computation of time-averages over thermostatted Ehrenfest dynamics trajectories lead to wrong results (as already noted in earlier works [15][16][17] ), we discuss ways to correct this issue. Section 4 discusses approaches based on classical-only MD propagations.

2
The canonical ensemble of hybrid quantum-classical systems We start by recalling which are the ensemble averages that we are addressing in this work.
The first step should be to clarify the mathematical description of a hybrid model, an issue that is not at all obvious, as demostrated by the various proposals that have been put forward, and by the discussions about their internal consistency. [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35] We will however assume the following very broad assumptions. The classical part is described by a set of position Q ∈ R n and momentum P ∈ R n variables, that we will hereafter collectively group as ξ = (Q, P ). Normally, they correspond to N particles, such that n = 3N in three dimensions. The quantum part is described by a complex Hilbert space H.
The first term is the kinetic energy of the classical particles, whereas the second part iŝ where the first term is the kinetic electronic operator, the second term is the electron-nucleus interaction potential, and the last (purely classical) term is the nucleus-nucleus interaction potential.
Ensembles of hybrid quantum-classical systems can be described 34-36 by ξ-dependent density matrices,ρ(ξ), normalized as: These fully characterize the ensemble, i.e. they permit to obtain the probabilities associated to any measurement. For example, the probability associated to finding the classical subsystem at ξ, and measuring a for observableÂ(ξ), is given by Tr [ρ(ξ)π a (ξ)], whereπ a (ξ) is the projector associated to the eigenvalue a ofÂ(ξ). Or, the probability density associated to the classical subsystem, regardless of the quantum part, is given by F C (ξ) = Trρ(ξ). Likewise, given any observableÂ, the ensemble average is given by: One route to the definition of equilibrium ensembles is the principle of maximization of entropy (MaxEnt). Recently, we argued 36 that the proper definition of entropy for a hybrid quantum-classical system must be: Likewise, we also showed that the maximization of this entropy, subject to the constraint of a given value for the average energy, leads to the hybrid canonical ensemble: Z HC (β) = dµ(ξ)Tre −βĤ(ξ) .
Therefore, the canonical ensemble average of any observable is: The computation of these averages is challenging. First, depending on the model and quantum level of theory used, the calculation of the traces, that in principle require all excited states, can be problematic. But, more importantly, the integral over the classical phase space is difficult because of its very large dimensionality (6N for N classical particles in 3D).
This latter problem is of course akin to the one encountered when studying purely classical systems. It is therefore natural to ask whether it is possible to circumvent it by doing some form of MD.
3 The failure of Ehrenfest dynamics, and some ways to correct it One possibility that immediately comes to mind is the use of a hybrid quantum-classical MD, such as Ehrenfest's, and attaching a thermostat in order to simulate the presence of a bath that would permit to generate the canonical ensemble along the trajectory. In other words, replicating the procedure invented for "standard" MD, but using a hybrid dynamics that requires the explicit propagation of the electrons.
It was soon realized, however, that this procedure leads to wrong ensemble averages. [15][16][17] In the following, we will reexamine this fact in the light of the Hamiltonian character of Ehrenfest dynamics. This analysis will help to understand the correction procedures that in fact permit to use this dynamics to obtain the true ensemble averages.

Fast recap of Hamiltonian dynamics
Let us first recap the basics of Hamiltonian theory. A system can be characterized by providing a phase space M of even dimension 2n. The Poisson bracket is an operation defined over functions in this phase space (the observables), which in the canonical coordinates (q, p) ∈ M reads: The dynamics is determined by the definition of a Hamiltonian function H: the equations of motion for the coordinates q i or p i of any state in or equivalently, as they are more often encountered, in the form of Hamilton's equations: If there is no certainty about the system state, instead of a single point, one must use a probability distribution ρ(q, p, t) defined over the phase space, also known as an ensemble.
This distribution may change in time, according to Liouville's equation: The entropy of any ensemble can be computed as (hereafter, we will group all variables q, p as y): The maximization of this entropy over all possible ensembles subject to the constraint of a given Hamiltonian ensemble average or energy, H ρ = dµ(y)H(y)ρ(y), leads to the canonical ensemble: 37 Z CC (β) = dµ(y) e −βH(y) .
Here, β = 1 k B T is inversely proportional to the temperature T , "CC" stands for "classical canonical", Z CC (β) is the partition function, and the integrals extend over all phase space. This is an equilibrium ensemble, lacking the time-dependence because it is station- The averages, for any observable A, over this canonical ensemble are then given by: The obvious numerical difficulty of computing these very high-dimensional integrals can then be circumvented by integrating a single dynamical trajectory, and using the ergodic hypothesis to identify a time average with the phase space integral: where y β (t) is a trajectory obtained by solving the equations of the motion, modified with a thermostat, i.e.:ẏ Here, we have symbolically added to the Poisson bracket {·, ·} a thermostat X β (t) (it may represent Langevin's stochastic term, 3,38 a Nose-Hoover chain, 39 etc.)

Schrödinger dynamics as a Hamiltonian system
The theory summarized in subsection 3.1 can be applied to any Hamiltonian dynamicsfor example, to Schrödinger's equation, which despite its quantum character, is a "classical" Hamiltonian system from a mathematical perspective. We summarize this fact here -for the mathematical conditions and functional spaces (both finite and infinite dimensional) on which this formalism can be applied, see [ 40,41]; in [ 42], one can follow the standard approach that is summarized here.
Indeed, Schrödinger's equation ( = 1 is assumed throughout this paper), is easy to rewrite as a set of Hamiltonian equations. First, one expands the wavefunction in an orthonormal basis {ϕ i } i , and rewrites Schrödinger's equation for the coefficients c i = ϕ i |ψ : where H ij are the elements of the Hamiltonian matrix H ij = ϕ i |Ĥ|ϕ j . We define a set of "position" and "momenta" variables by taking the real and imaginary parts of this coefficients, i.e. Hamiltonian's equations, defining H(q, p) as the expectation value ofĤ for the wavefunction determined by the (q, p) coefficients: H(q, p) = ψ(q, p)|Ĥ|ψ(q, p) . Note that these "position" and "momemtum" variables should not be given any particular physical meaning, forcing any analogy with classical mechanics.
The Poisson bracket can then be defined in the usual way [Eq. (9)]. We will use the notation {·, ·} Q for this Poisson bracket defined in this new "quantum" phase space, M Q , defined by the variables (q, p). The system dynamics then reduces to: for any function f defined in M Q . For the particular case of the coordinate functions q and p, one obtains the Hamilton equations above. Taking into account the dependence of ψ(q, p), they are entirely equivalent to Schrödinger's equation (see 41,42 for details).
Gibbs canonical ensemble associated with the expectation value of the energy H(q, p), Eq. (14), is stationary under the dynamics, and in principle one could attach one of the typical classical-type thermostats to Schrödinger equation, and produce this ensemble through a trajectory. If one did that, however, the resulting ensemble averages would be wrong, because Gibbs ensemble is not the true quantum canonical ensemble, which is defined through a density matrix as:ρ This is the density matrix that maximizes the von Neumann entropy, which is the real entropy of a quantum system, and not Eq. (13), from which the classical canonical ensemble is derived. It is obvious that the fact that the "thermostatted" Schrödinger dynamics does not produce the correct thermal averages is not a defect of Schrödinger equation, but results of the erroneous application of a technique invented for classical systems to quantum ones.

Ehrenfest dynamics as a Hamiltonian system
Ehrenfest dynamics, usually introduced as a partial classical limit of the full-quantum dynamics 4 , constitutes also a Hamiltonian system, as shown for example in [ 4,43,44]. We summarize here this fact.
We consider a hybrid quantum-classical system, as defined in Section 2. The classical variables ξ = (Q, P ) define a classical phase space M C , whereas the quantum variables η = (q, p), associated to a wavefunction ψ(η) as explained above, define a quantum phase space M Q . We may put these together and define a full, hybrid phase space: Furthermore, we may define a Poisson bracket for functions defined on this full hybrid space by adding the two classical and quantum brackets, defined over the ξ and η variables, respectively: Thus, the classical bracket derivates only with respect to the classical coordinates (Q, P ), and the quantum bracket with respect to the quantum ones (q, p). It is a well known result of Poisson geometry that the addition of two brackets in a Cartesian product results in a bracket that fulfills all the necessary properties.
Finally, given any hybrid observableÂ(ξ), we may define a real function over M as: If, in particular, we consider the Hamiltonian operatorĤ(ξ), which is dependent on the classical degrees of freedom ξ, the hybrid dynamics is generated by the function H(η, ξ) = ψ(η)|Ĥ(ξ)|ψ(η) and the Poisson bracket: where in the last equality of both lines the classical and quantum nature of ξ and η respectively has been invoked to make use of {ξ, f } Q = {η, f } C = 0 ∀f ∈ C ∞ (M).
One may further expand those equations: for the hybrid Poisson bracket acting on the classical variables ξ i , one has: If one then considers the cases ξ i = Q k or ξ i = P k separately, one arrives to Newton's-like On the other hand, the dynamical equation of the the quantum variables, is exactly the same as Eq. (23), but with a dependence of the Hamiltonian operator on the classical variables. Therefore, this equation, as it was shown in the previous section for the quantum-only case, is equivalent to Schrödinger's equation for ψ(η) -although one must maintain that parametric dependence of the Hamiltonian operator on the classical variables ξ = Q, P : Eqs. (33)(34) and (36) (1) and (2) for a set of electrons and nuclei, one gets: Allured by the Hamiltonian character of this set of equations, one may be tempted to consider the Gibbs equilibrium ensemble, Eq. (14), to be the hybrid canonical one. In terms of the quantum-classical variables, it reads: It is also a stationary ensemble in the hybrid case. The averages over this ensemble would be the ones obtained if one attaches a thermostat tuned to temperature T = 1 k B β to Ehrenfest dynamics, propagates a trajectory (η β (t), ξ β (t)), and computes the time averages: For example, one practical way to proceed is to use Langevin's dynamics (although there are various other thermostat definitions that have been invented over the years), that essentially consists in substituting the equation for the force (38) by: where η I (t) are stochastic Gaussian processes that must verify: The α, β indices run over the three spatial dimensions; see for example Ref. 3  The underlying reason behind the difference of the two ensembles is that, in classical systems, all points in the phase space are mutually exclusive, and are given a Boltzmann weight in the canonical ensemble. The MD procedure (and in particular the thermostats) was designed to produce a phase space visitation consistent with this. In the hybrid phase space, due to the quantum character of one of its parts, not all distinct points are mutually exclusive events, 36 and therefore the ensemble targetted by the thermostats, determined by a Boltzmann weight over the phase space, does not match the HC ensemble.

Corrected averages for Ehrenfest dynamics
Nevertheless, Eq. (42) can be useful, as we will show now. The thermostatted Ehrenfest dynamics does sample the phase space, and it generates an ensemble, even if wrong. One may then apply a reweighting procedure -essentially, modifying the averaging in the time integral -and obtain the correct hybrid ensemble averages. This can be done in fact in several ways.
The first thing to notice is that Eq. (42) holds for any function g(η, ξ) on M, not only on the ones that result of a hybrid observable as A(η, ξ) = η|Â(ξ)|η . Then one may ask the question: for any hybrid observableÂ, can one find a function gÂ(η, ξ), such that: If so, one could then perform the dynamics and use Eq. (42) with gÂ in order to obtain gÂ CC (β), and therefore the true hybrid ensemble average Â HC (β).
The answer is positive, and there is not only one, but many possible functions that can be used. In the following, we consider two examples: 1. Equation (46) holds if gÂ is defined as: .
On the other hand, we know that for the identity operator, Î HC (β) = 1, and therefore: Thus, we may compute µ(β) as 1/ gÎ CC (β), and gÎ CC (β) can be obtained from a dynamics propagation, i.e.: Summarizing, a final formula that permits to compute the hybrid ensemble averages is: Therefore, the procedure consists of performing a thermostatted Ehrenfest dynamics, and computing the previous time integrals over the obtained trajectory (η β (t), ξ β (t)).
One obvious difficulty lies in the computation of the traces over the quantum Hilbert space, whose difficulty depends on the level of theory used to deal with the quantum electronic problem.

Equation (46) also holds if gÂ is defined as:
where η α (ξ) are the adiabatic states: and The difficulty due to the computation of the λ(β) factor can be solved in a similar way to the method used in the previous case, leading to the following final formula: .

(57)
This formula avoids the need to compute all the electronic excited states, necessary for the traces present in Eq. (52). In exchange, it contains a probably worse numerical difficulty: the presence of the delta functions. The interpretation of these is the following: during the trajectories, one should not count in the average the state that is being visited, unless the trajectory passes by an eigenstate of the Hamiltonian (a state of the adiabatic basis). In other words, apart from the normalization factor given by the denominator, this formula is a modification of the straigthforward average given in Eq. (42), that discards all states except for the adiabatic eigenstates.
That correction is easy to understand intuitively. Let us first rewrite the hybrid canonical ensemble density matrix,ρ in terms of its spectral decomposition for each ξ: where E α (ξ) are the eigenvalues, andη α (ξ) the projectors on the eigenspaces of the HamiltonianĤ(ξ) (we assume, for simplicity, that there is no degeneration; otherwise one would just need to use an orthogonal basis for each degenerate subspace). Now, this expression can be written in terms of a (generalized) probability distribution function in M, as: This distribution determinesρ HC (ξ), since: ρ HC (ξ) = dµ(η)ρ HC (η, ξ) |η η| η|η . In this section, we show how this problem can be circumvented by making use of dynamics that do not explicitly propagate the electrons, such as ground-state Born-Oppenheimer MD -including the necessary correction to account for the hot electrons -, or the dynamics based on the electronic free energy surface that has already been used in the past. In this way, we frame these approaches into the theoretical setup described above.
Let us suppose that we perform a MD for the classical particles, based on a Hamiltonian function H(ξ) (to be specified below). In this case, the dynamics is not hybrid: the propagation equations involve only the classical particles, moving under the influence of H(ξ).
The ergodic assumption, if it holds, permits to compute: for any function g(ξ). Notice that the classical canonical ensemble that we are using now refers to the classical degrees of freedom ξ only, as opposed to the one used in the previous section, that included the quantum ones.
In the same manner as we did in the previous section, one may wonder the following: for a given hybrid observableÂ(ξ), does there exist some function gÂ(ξ) such that Once again, the answer is affirmative, and in more ways than one. One obvious possibility, analogous to the first one used for Ehrenfest dynamics, is: As it happened in the previous section, the calculation of the normalization factor µ(β) does not require of the explicit computation of the partition functions (that may be impractical), but may result from the MD propagation itself, using the identity Î HC (β) = 1. Using this fact and the same procedure shown in the previous section, one arrives to the final formula: .

(66)
This formula is very similar to Eq. (52). However, the trajectory ξ β (t) to be used here must be obtained through a thermostatted classical-only MD determined by a Hamiltonian function H(ξ), in contrast to the hybrid quantum-classical Ehrenfest dynamics used in the previous section. The Hamiltonian function H(ξ) is in fact arbitrary, although a bad choice for this object could lead to a very bad convergence with respect to the total propagation time t f -since using the ergodic hypothesis requires an accurate sampling of phase space.
Two options that immediately come to mind are: 1. Using the electronic free energy: This actually permits to simplify Eq. (66) into a very appealing form: where at each classical phase space point in the trajectory one must compute the quantum ensemble average: Eq. (68) reminds of the usual MD ergodic averaging formula, just substituting the observable value by the thermal quantum average. In fact, if the observable that one is interested in is purely classical,Â(ξ) = A(ξ)Î, the formula is identical: Therefore, for purely classical observables, if one uses the electronic free energy instead of the ground-state adiabatic energy as the Hamiltonian driving the ionic movement, the resulting MD provides the hybrid canonical averages using the "standard" ergodic average. If the observable is itself hybrid, one must compute at each point during the trajectory the quantum thermal average.
This propagation of the classical variables following the electronic free energy surface underlies the scheme put forward by Alavi et al., 6 although in that work the procedure is tightly tied to the use of FTDFT as the scheme that handles the electronic structure problem (computation of the free energy, and of its gradients). The same concept was also suggested by some of the current authors in Refs. 10,11.
2. Using the ground-state Born-Oppenheimer energy: In this case we would just need to do the usual ground-state Born-Oppenheimer MD, which has the advantage of being a very well known and tested technique, for which plenty of codes and tools exist. In order to obtain the hybrid ensemble averages that do not ignore the electronic temperature, however, one must use the averaging formula (66), which for this case can be transformed into: .
where α runs over all the adiabatic eigenstates, and are the electronic excitations.
On top of the usual ground-state Born-Oppenheimer MD, the added difficulty here would be the computation of these excitations, which may be more or less demanding depending on the level of theory used to model the many-electron problem.
Note that if the observableÂ(ξ) is actually a classical observable A(ξ)Î, this scheme can also be rewritten as: .
Here, we also write the formula in terms of the free energy. Computationally, the difference with respect to the previous approach given in formula (70) is that one does not need the gradients of the free energy, necessary in the previous approach for the computation of the forces in the dynamics.
All previous formulas have assumed that the thermostat is fixed to the target temperature. However, the dynamics can be performed at a different temperature (a technique that has been used in MD to probe larger regions of configuration space in less simulation time), as long as the reweighting corrects for this. Take, for example, formula (66), that we repeat here for convenience, although we now use two different temperatures β and β : .
In this formula the temperature dependence is twofold: 1. The temperature used to define the thermostat, which appears in the β labeling the trajectories ξ β (t). This should be equal to the temperature used for defining the re-weighting factors, i.e. the inverse Boltzmann weight e β H(ξ β (t)) .
2. The "target" temperature, that is the one that should be used in the exponent of the un-normalized hybrid canonical ensemble density matrix, appearing inside the trace, These two temperatures can be different, and formula (75) still holds. In practical applications, this would permit to obtain hybrid canonical ensemble averages at different temperatures β by computing a single thermostatted trajectory at a fixed "ergodic temperature" β .
This can also be done when using Ehrenfest dynamics, and formulas (52) and (57)  one would then need to compute the free energy at those two temperatures). Note also that, in practice, this procedure cannot be indefinitely extended to any temperature range, since the ergodic visitation will not be effective unless the two temperatures are similar.
Summarizing, the previous formulas permit to use well known MD techniques and obtain canonical averages that correctly account for the electronic temperature. Looking, for example, at Eq. (74), the procedure entails two steps: 1. One first performs a standard first principles MD simulation using, for example, the common technique based on ground-state DFT.
2. Then, either on the fly as the trajectory is being generated, or later in a post-processing procedure, one computes the electronic free-energy at the trajectory points, using the finite-temperature DFT extension. With such information, one can use Eq. (74) to correct the time averages that, without this averaging method, would fail to converge to the real canonical ensemble.
This scheme can be applied on top of trajectories obtained previously, had they been saved.
Note that, for the reasons explained above, one may recycle trajectories obtained with ground-state BOMD at some (nuclear only) temperature, to compute ensemble averages at various different global temperatures. This may be an advantage over the procedure implied by Eq. (70) (MD with forces computed on the electronic free-energy surface), as in that case one trajectory must be generated at each temperature. Finally, of course DFT need not be the method to be used for the computation of the forces -and the finite temperature DFT need not be the procedure to obtain the free energy, as one may use for example TDDFT to compute the electronic excitations and apply formula (72).

Numerical example
We finish with a numerical demonstration of the validity of the formulas given above, using a simple model and the very last of the presented schemes: "standard" ground-state Born-Oppheheimer MD with a correction formula. Thus, we consider a simple dimer model, using the internuclear distance Q as the only classical position variable (being P the corresponding momentum), and the subspace generated by the two lowest electronic adiabatic states as the quantum space. Futhermore, we consider that these two adiabatic states correspond to Morse potentials. Hence, the Hamiltonian operator ruling the hybrid dynamics can be written, in the basis of its eigenstates, as: where m is the dimer reduced mass, and the Morse potentials are given by: The parameters defining the Morse potential for the i-th energy level V i (Q) have an easy interpretation: ∆ i is a global shift that sets the value of the curve at its minimum; q i the position at that minimum (V i (Q = q i ) = ∆ i ); D i defines how quickly the potential ascends for Q > q i , and also determines the value of the gap between the minimum (∆ i ) and the big Q limit of the potential: lim Q→∞ V i (Q) = ∆ i + D i . Lastly, b i defines how narrow the well is, how sharply it grows when Q → 0, and also how rapidly it reaches the plateau for Q > q i .
The vibrational frequency associated to each potential well is given by ω i = 2b 2 D m . Using ground-state Born-Oppenheimer MD means that the classical degrees of freedom follow the Hamiltonian system defined by the function: The system is coupled to a Langevin thermostat, at the temperature given by T = 1 k B β . This dynamics provides an ergodic curve over the classical phase space with a visitation weight given by the Boltzmann factor e −βE 0 (Q,P ) . Using the corrected averaging procedure defined by Eq. (72), we can obtain the hybrid canonical ensemble averages. Although this formula is valid for any observable, we chose to compute the average value of the length of the dimer, a purely classical observable:Â(Q, P ) = QÎ.  We stress that, in the procedure presented above, the computation of the ergodic trajectory is an indirect way to perform phase space integrals over the classical phase space.
In principle, any infinite (t f → ∞) ergodic trajectory could be used as a basis to apply the corrected averaging procedure, as long as the implicit distribution over the phase space that results of the dynamics is compensated in the time averages: the use of the ground-state potential energy surface to generate the dynamics is one of the possible many choices. This only holds if the trajectory provides a dense enough visitation of the phase space.
However, in practice, the simulations provide only finite-time trajectories and, therefore, the visitation of phase space is not dense. The effectiveness of a given dynamics will depend on what regions of phase space it probes more frequently. Some thermostatted MD trajectories will be more cost-effective than others, if they visit more frequently the regions of the phase space that are relevant to the target distribution.
In our example, the target distribution is the hybrid canonical ensemble, and one should choose a dynamics that is likely to force the system to spend time on its high probability regions. This is not the only factor to consider, however. For example, it is likely that Ehrenfest dynamics fulfills this condition, but the cost of propagating Ehrenfest equations is high, due to the need to propagate the electrons. Likewise, it may happen that using the free energy as the driving Hamiltonian is costly due to the requirement of computing its gradients with respect to the classical degrees of freedom in order to obtain the forces. A dynamics that requires longer times t f to achieve the convergence of the time average can be however computationally cheaper if the cost of performing the propagation itself is lower. We consider that ground-state MD can be a good compromise, specially at low temperatures, but the analysis strongly depends on the particular model, the electronic structure method, etc.

Conclusions
We have examined the problem of computing the canonical ensemble averages through MD calculations, for hybrid quantum classical systems (typically, quantum electrons and classical nuclei in molecular or condensed matter physics and chemistry). If the temperature is high enough so that the electronic excited states cannot be ignored, performing ground state Born Oppenheimer MD and computing the ergodic averages on the generated trajectories does not yield the correct ensemble averages.
The fact that one cannot assume that the electrons are inert, adiabatically adapting to the ground state, naturally seems to demand for a truly hybrid dynamics, such as Ehrenfest's.
However, the addition of a thermostat to these equations, followed by the computation of the resulting observable time-averages, does not produce the right averages either. The quantum character of part of the electrons cannot be handled by the standard MD + thermostat procedure, designed to produce essentially classical equilibrium ensembles.
Nevertheless, performing a thermostatted dynamics, be it Ehrenfest, or a purely classical one such as ground state MD, does generate a trajectory (i.e. an ensemble) in phase space that can be reweighted in order to obtain the true hybrid averages. This amounts to correcting the time averaging formulas. It has been the purpose of this work to examine here the various options that exist, setting them in a common language. The procedure can be followed using Ehrenfest dynamics, which is a hybrid dynamics, but can also be followed using classical-only dynamics driven by, for example, the ground-state Born-Oppenheimer Hamiltonian. Likewise, this framework does include, as a particular case, the possibility of performing the nuclear dynamics on the Hamiltonian that results of considering the electronic free energy. In this case, if one is interested in computing averages of purely classical observables, the time averages do not require correction, as the factors cancel out. The suitability of any of these procedures over the others depends on the particular model and level of theory used to handle the electronic structure problem.