About the computation of finite temperature ensemble averages of hybrid quantum-classical systems with molecular dynamics

Let us first recap the basics of Hamiltonian theory. A system can be characterized by providing a phase space $\mathcal{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 $\left(q,p\right)\in \mathcal{M}$ reads:

$$\begin{equation}\left\{A,B\right\}=\sum\limits _{i=1}^{n}\left[\frac{\partial A}{\partial {q}_{i}}\frac{\partial B}{\partial {p}_{i}}-\frac{\partial A}{\partial {p}_{i}}\frac{\partial B}{\partial {q}_{i}}\right].\end{equation} \tag{ 9 }$$

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 $\mathcal{M}$ are ${\dot {q}}_{i}=\left\{{q}_{i},H\right\}$ and ${\dot {p}}_{i}=\left\{{p}_{i},H\right\}$, or equivalently, as they are more often encountered, in the form of Hamilton's equations:

$$\begin{equation}{\dot {q}}_{i}=\frac{\partial H}{\partial {p}_{i}},\end{equation} \tag{ 10 }$$

$$\begin{equation}{\dot {p}}_{i}=-\frac{\partial H}{\partial {q}_{i}}.\end{equation} \tag{ 11 }$$

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:

$$\begin{equation}\frac{\partial \rho }{\partial t}=\left\{H,\rho \right\}.\end{equation} \tag{ 12 }$$

The entropy of any ensemble can be computed as (hereafter, we will group all variables q, p as y):

$$\begin{equation}S\left[\rho \right]=-{k}_{\text{B}}\int \mathrm{d}\mu \left(y\right)\rho \left(y\right)\mathrm{log}\enspace \rho \left(y\right).\end{equation} \tag{ 13 }$$

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 [41]:

$$\begin{equation}{\rho }_{\text{CC}}\left(y\right)=\frac{1}{{Z}_{\text{CC}\;}\left(\beta \right)}{\text{e}}^{-\beta H\left(y\right)},\end{equation} \tag{ 14 }$$

$$\begin{equation}{Z}_{\text{CC}}\left(\beta \right)=\int \mathrm{d}\mu \left(y\right){\text{e}}^{-\beta H\left(y\right)}.\end{equation} \tag{ 15 }$$

Here, $\beta =\frac{1}{{k}_{\text{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 stationary: {H, ρ_CC} = 0.

The averages, for any observable A, over this canonical ensemble are then given by:

$$\begin{equation}\langle A{\rangle }_{\text{CC}}\left(\beta \right)=\frac{1}{{Z}_{\text{CC}}\left(\beta \right)}\int \mathrm{d}\mu \left(y\right){\text{e}}^{-\beta H\left(y\right)}A\left(y\right).\end{equation} \tag{ 16 }$$

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:

$$\begin{equation}\langle A{\rangle }_{\text{CC}}\left(\beta \right)=\underset{{t}_{f}\to \infty }{\mathrm{lim}}\enspace \frac{1}{{t}_{f}}{\int }_{0}^{{t}_{f}} \mathrm{d}t\enspace A\left({y}^{\beta }\left(t\right)\right),\end{equation} \tag{ 17 }$$

where y^β (t) is a trajectory obtained by solving the equations of the motion, modified with a thermostat, i.e.:

$$\begin{equation}{\dot {y}}_{i}^{\beta }=\left\{{y}_{i}^{\beta },H\right\}+{X}_{i}^{\beta }\left(t\right).\end{equation} \tag{ 18 }$$

Here, we have symbolically added to the Poisson bracket {⋅, ⋅} a thermostat X_β (t) (it may represent Langevin's stochastic term [29, 45], a Nose–Hoover chain [30], etc.)

The theory summarized in subsection 3.1 can be applied to any Hamiltonian dynamics—for 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 [16, 26]; in [24], one can follow the standard approach that is summarized here.

Indeed, Schrödinger's equation (ℏ = 1 is assumed throughout this paper),

$$\begin{equation}\text{i}\frac{\mathrm{d}}{\mathrm{d}t}\vert \psi \left(t\right)\rangle =\hat{H}\vert \psi \left(t\right)\rangle ,\end{equation} \tag{ 19 }$$

is easy to rewrite as a set of Hamiltonian equations. First, one expands the wavefunction in an orthonormal basis ${\left\{{\varphi }_{i}\right\}}_{i}$, and rewrites Schrödinger's equation for the coefficients c_i = ⟨φ_i |ψ⟩:

$$\begin{equation}\text{i}\frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{matrix}\hfill {c}_{1}\hfill \\ \hfill {c}_{2}\hfill \\ \hfill {\vdots}\hfill \\ \hfill {c}_{n}\hfill \end{matrix}\right)=\left(\begin{matrix}\hfill {H}_{11}\hfill & \hfill {H}_{12}\hfill & \hfill \dots \hfill & \hfill {H}_{1n}\hfill \\ \hfill {H}_{21}\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill {H}_{2n}\hfill \\ \hfill {\vdots}\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill {\vdots}\hfill \\ \hfill {H}_{n1}\hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill {H}_{nn}\hfill \end{matrix}\right)\left(\begin{matrix}\hfill {c}_{1}\hfill \\ \hfill {c}_{2}\hfill \\ \hfill {\vdots}\hfill \\ \hfill {c}_{n}\hfill \end{matrix}\right),\end{equation} \tag{ 20 }$$

where H_ij are the elements of the Hamiltonian matrix ${H}_{ij}=\langle {\varphi }_{i}\vert \hat{H}\vert {\varphi }_{j}\rangle$. We define a set of 'position' and 'momenta' variables by taking the real and imaginary parts of this coefficients, respectively: $c=\frac{1}{\sqrt{2}}\left(q+\text{i}p\right)$. One may then show [24, 26] that equation (19) is equivalent to

$$\begin{equation}{\dot {q}}_{i}=\frac{\partial H}{\partial {p}_{i}},\end{equation} \tag{ 21 }$$

$$\begin{equation}{\dot {p}}_{i}=-\frac{\partial H}{\partial {q}_{i}},\end{equation} \tag{ 22 }$$

i.e. Hamiltonian's equations, defining H(q, p) as the expectation value of $\hat{H}$ for the wavefunction determined by the (q, p) coefficients: $H\left(q,p\right)=\langle \psi \left(q,p\right)\vert \hat{H}\vert \psi \left(q,p\right)\rangle$. Note that these 'position' and 'momentum' 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 (equation (9)). We will use the notation {⋅,⋅}_Q for this Poisson bracket defined in this new 'quantum' phase space, ${\mathcal{M}}_{Q}$, defined by the variables (q, p). The system dynamics then reduces to:

$$\begin{equation}\dot {f}={\left\{\enspace f,H\right\}}_{Q},\end{equation} \tag{ 23 }$$

for any function f defined in ${\mathcal{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 [24, 26] for details).

Gibbs canonical ensemble associated with the expectation value of the energy H(q, p), equation (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:

$$\begin{equation}{\hat{\rho }}_{\text{QC}}=\frac{1}{{Z}_{\text{QC}}\left(\beta \right)}{\text{e}}^{-\beta \hat{H}},\end{equation} \tag{ 24 }$$

$$\begin{equation}{Z}_{\text{QC}}\left(\beta \right)=\mathrm{Tr}\enspace {\text{e}}^{-\beta \hat{H}}.\end{equation} \tag{ 25 }$$

This is the density matrix that maximizes the von Neumann entropy,

$$\begin{equation}S\left[\hat{\rho }\right]=-{k}_{\text{B}}\enspace \mathrm{Tr}\left(\hat{\rho }\enspace \mathrm{log}\enspace \hat{\rho }\right),\end{equation} \tag{ 26 }$$

which is the real entropy of a quantum system, and not equation (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, usually introduced as a partial classical limit of the full-quantum dynamics [11], constitutes also a Hamiltonian system, as shown for example in [6, 7, 11]. 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 ${\mathcal{M}}_{\text{C}}$, whereas the quantum variables η = (q, p), associated to a wavefunction ψ(η) as explained above, define a quantum phase space ${\mathcal{M}}_{Q}$. We may put these together and define a full, hybrid phase space:

$$\begin{equation}\mathcal{M}={\mathcal{M}}_{\text{C}}{\times}{\mathcal{M}}_{Q}.\end{equation} \tag{ 27 }$$

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:

$$\begin{align}{\left\{A,B\right\}}_{\text{H}}& ={\left\{A,B\right\}}_{\text{C}}+{\left\{A,B\right\}}_{Q}\\ & =\sum\limits _{i}\left({\partial }_{{Q}_{i}}A{\partial }_{{P}_{i}}B-{\partial }_{{P}_{i}}A{\partial }_{{Q}_{i}}B\right)+\sum\limits _{i}\left({\partial }_{{q}_{i}}A{\partial }_{{p}_{i}}B-{\partial }_{{p}_{i}}A{\partial }_{{q}_{i}}B\right).\end{align} \tag{ 28 }$$

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 $\hat{A}\left(\xi \right)$, we may define a real function over $\mathcal{M}$ as:

$$\begin{equation}A\left(\eta ,\xi \right)=\langle \psi \left(\eta \right)\vert \hat{A}\left(\xi \right)\vert \psi \left(\eta \right)\rangle .\end{equation} \tag{ 29 }$$

If, in particular, we consider the Hamiltonian operator $\hat{H}\left(\xi \right)$, which is dependent on the classical degrees of freedom ξ, the hybrid dynamics is generated by the function $H\left(\eta ,\xi \right)=\langle \psi \left(\eta \right)\vert \hat{H}\left(\xi \right)\vert \psi \left(\eta \right)\rangle$ and the Poisson bracket:

$$\begin{equation}{\dot {\xi }}_{i}={\left\{{\xi }_{i},H\right\}}_{\text{H}}={\left\{{\xi }_{i},H\right\}}_{\text{C}},\end{equation} \tag{ 30 }$$

$$\begin{equation}{\dot {\eta }}_{a}={\left\{{\eta }_{a},H\right\}}_{\text{H}}={\left\{{\eta }_{a},H\right\}}_{\mathit{\text{Q}}},\end{equation} \tag{ 31 }$$

where in the last equality of both lines the classical and quantum nature of ξ and η respectively has been invoked to make use of ${\left\{\xi ,f\right\}}_{Q}={\left\{\eta ,f\right\}}_{\text{C}}=0\forall \enspace f\in {C}^{\infty }\left(\mathcal{M}\right)$.

One may further expand those equations: for the hybrid Poisson bracket acting on the classical variables ξ_i , one has:

$$\begin{equation}{\dot {\xi }}_{i}={\left\{{\xi }_{i},H\right\}}_{\text{C}}=\sum\limits _{j}\left(\left({\partial }_{{Q}_{j}}{\xi }_{i}\right)\langle \psi \vert {\partial }_{{P}_{j}}\hat{H}\left(\xi \right)\vert \psi \rangle -\left({\partial }_{{P}_{j}}{\xi }_{i}\right)\langle \psi \vert {\partial }_{{Q}_{j}}\hat{H}\left(\xi \right)\vert \psi \rangle \right).\end{equation} \tag{ 32 }$$

If one then considers the cases ξ_i = Q_k or ξ_i = P_k separately, one arrives to Newton's-like equations for Q_i , P_i :

$$\begin{equation}{\dot {Q}}_{k}=\langle \psi \vert \frac{\partial \hat{H}}{\partial {P}_{k}}\vert \psi \rangle ,\end{equation} \tag{ 33 }$$

$$\begin{equation}{\dot {P}}_{k}=-\langle \psi \vert \frac{\partial \hat{H}}{\partial {Q}_{k}}\vert \psi \rangle .\end{equation} \tag{ 34 }$$

On the other hand, the dynamical equation of the quantum variables,

$$\begin{equation}{\dot {\eta }}_{a}={\left\{{\eta }_{a},H\left(\xi ,\eta \right)\right\}}_{Q},\end{equation} \tag{ 35 }$$

is exactly the same as equation (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:

$$\begin{equation}\frac{\mathrm{d}}{\mathrm{d}t}\vert \psi \rangle =-\text{i}\hat{H}\left(Q,P\right)\vert \psi \rangle .\end{equation} \tag{ 36 }$$

Equations (33), (34) and (36) are Ehrenfest's equations for a hybrid model. From our previous analysis of the classical and the quantum case, it becomes clear that they compose a Hamiltonian system with the Poisson bracket defined as in equation (28). Perhaps the form given in equations (33), (34) and (36) is still not the most recognizable; if one uses the Hamiltonian defined in equations (1) and (2) for a set of electrons and nuclei, one gets:

$$\begin{equation}{\dot { \overrightarrow {Q}}}_{I}=\frac{{ \overrightarrow {P}}_{I}}{{M}_{I}}\end{equation} \tag{ 37 }$$

$$\begin{equation}{\dot { \overrightarrow {P}}}_{I}=-\langle \psi \vert { \overrightarrow {\nabla }}_{I}{\hat{H}}_{\mathrm{e}}\left(Q\right)\vert \psi \rangle ,\end{equation} \tag{ 38 }$$

$$\begin{equation}\frac{\mathrm{d}}{\mathrm{d}t}\vert \psi \rangle =-\text{i}\hat{H}\left(Q,P\right)\vert \psi \rangle .\end{equation} \tag{ 39 }$$

Allured by the Hamiltonian character of this set of equations, one may be tempted to consider the Gibbs equilibrium ensemble, equation (14), to be the hybrid canonical one. In terms of the quantum-classical variables, it reads:

$$\begin{equation}{\rho }_{\text{CC}}\left(\eta ,\xi \right)=\frac{1}{{Z}_{\text{CC}\;}\left(\beta \right)}{\text{e}}^{-\beta H\left(\eta ,\xi \right)},\end{equation} \tag{ 40 }$$

$$\begin{equation}{Z}_{\text{CC}}\left(\beta \right)=\int \mathrm{d}\mu \left(\xi \right)\mathrm{d}\mu \left(\eta \right)\enspace {\text{e}}^{-\beta H\left(\eta ,\xi \right)}.\end{equation} \tag{ 41 }$$

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=\frac{1}{{k}_{\text{B}}\beta }$ to Ehrenfest dynamics, propagates a trajectory (η^β (t), ξ^β (t)), and computes the time averages:

$$\begin{align}\langle A{\rangle }_{\text{CC}}\left(\beta \right)& =\underset{{t}_{f}\to \infty }{\mathrm{lim}}\enspace \frac{1}{{t}_{f}}{\int }_{0}^{{t}_{f}} \mathrm{d}t\enspace A\left({\eta }^{\beta }\left(t\right),{\xi }^{\beta }\left(t\right)\right)\\ & =\frac{1}{{Z}_{\text{CC}\;}\left(\beta \right)}\int \mathrm{d}\mu \left(\xi \right)\mathrm{d}\mu \left(\eta \right)A\left(\eta ,\xi \right){\text{e}}^{-\beta H\left(\eta ,\xi \right)}.\end{align} \tag{ 42 }$$

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:

$$\begin{equation}{\dot { \overrightarrow {P}}}_{I}=-\langle \psi \vert { \overrightarrow {\nabla }}_{I}{\hat{H}}_{\mathrm{e}}\left(Q\right)\vert \psi \rangle -\beta \gamma \frac{{ \overrightarrow {P}}_{I}}{{M}_{I}}+{ \overrightarrow {\eta }}_{I}\left(t\right),\end{equation} \tag{ 43 }$$

where ${ \overrightarrow {\eta }}_{I}\left(t\right)$ are stochastic Gaussian processes that must verify:

$$\begin{equation}\langle { \overrightarrow {\eta }}_{I}\left(t\right)\rangle =0,\end{equation} \tag{ 44 }$$

$$\begin{equation}\langle {\eta }_{I\alpha }\left(t\right){\eta }_{J\beta }\left({t}^{\prime }\right)\rangle =2\gamma {\delta }_{IJ}{\delta }_{\alpha \beta }\delta \left(t,{t}^{\prime }\right).\end{equation} \tag{ 45 }$$

The α, β indices run over the three spatial dimensions; see for example reference [45] for details. In any case, the values thus obtained are not the ensemble average values that one would wish to obtain, given above in equation (8), hence the previously documented numerical failure of this approach—see for example references [32, 34, 35]. It should be noted, however, that this fact by itself should not be considered a failure of Ehrenfest dynamics—inasmuch as the same fact noted above for the quantum case cannot be considered a failure of Schrödinger equation. It results, once again, of the erroneous application of a technique invented for purely classical systems to hybrid ones, that contain some quantum variables.

A correct statistical-mechanical definition of a system requires the definition of a sample space, i.e., a basis of mutually exclusive events, which can be unequivocally characterized by the results of an experiment (see [5, 8, 14, 21]). The underlying reason behind the difference between classical and hybrid (or quantum) ensembles is that, in classical systems, all points in the phase space are mutually exclusive, whereas they are not in the quantum and hybrid case. Indeed, in classical mechanics, two different points of the phase space represent statistical events which are mutually exclusive: if the system is at state ξ, one cannot measure it to be at state ξ'. One may then define an entropy function using the classical phase space as sample space (i.e. a space that contains mutually exclusive events), and its maximization, subject to a fixed ensemble energy, leads to the classical Gibbs canonical ensemble.

In contrast, in quantum mechanics, two different states of the Hilbert space are not always mutually exclusive: if the system is at state ϕ, the probability of measuring it to be at ϕ' is not zero, unless the two states are orthogonal, ⟨ϕ|ϕ'⟩ = 0. Therefore, using the standard theory of probability, we should consider only orthogonal states to represent events (see sections 4.5 and 4.6 of [21]). This consideration led von Neumann and others to the introduction of the concept of density matrix in order to define a statistical ensemble—and to the definition of an entropy based on these objects, whose maximization results in the quantum canonical ensemble.

This precaution about the exclusivity of events must also be taken in the hybrid case—a fact that is discussed at more length in reference [5]. The points of the hybrid phase space are not mutually exclusive unless their quantum parts are orthogonal, or their classical parts are different. It should be noticed, however, that the MD procedure (and in particular the thermostats) produces a visitation of the phase space that disregards this fact. Thus, a MD trajectory visits (assuming ergodicity) all points in a hybrid phase space, assigning to each one a Boltzmann weight, and therefore producing a purely classical canonical ensemble that does not take into account the non-exclusivity of hybrid events.

Nevertheless, equation (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 equation (42) holds for any function g(η, ξ) on $\mathcal{M}$, not only on the ones that result of a hybrid observable as $A\left(\eta ,\xi \right)=\langle \eta \vert \hat{A}\left(\xi \right)\vert \eta \rangle$. Then one may ask the question: for any hybrid observable $\hat{A}$, can one find a function ${g}_{\hat{A}}\left(\eta ,\xi \right)$, such that:

$$\begin{equation}\langle {\hat{A}\rangle }_{\text{HC}}\left(\beta \right)=\langle {g/A\hat{A}\rangle }_{\text{CC}}\left(\beta \right)\enspace ?\end{equation} \tag{ 46 }$$

If so, one could then perform the dynamics and use equation (42) with ${g}_{\hat{A}}$ in order to obtain $\langle {g}_{\hat{A}}{\rangle }_{\text{CC}}\left(\beta \right)$, and therefore the true hybrid ensemble average $\langle \hat{A}{\rangle }_{\text{HC}}\left(\beta \right)$.

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:

• (a)  
  Equation (46) holds if ${g}_{\hat{A}}$ is defined as:

  $$\begin{equation}{g}_{\hat{A}}\left(\eta ,\xi \right)=\mu \left(\beta \right){\text{e}}^{\beta H\left(\eta ,\xi \right)}\enspace \mathrm{Tr}\left[{\text{e}}^{-\beta \hat{H}\left(\xi \right)}\hat{A}\left(\xi \right)\right],\quad \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\end{equation} \tag{ 47 }$$

  $$\begin{equation}\mu \left(\beta \right)=\frac{{Z}_{\text{CC}}\left(\beta \right)}{\left(\int \mathrm{d}\mu \left(\eta \right)\right){Z}_{\text{HC}}\left(\beta \right)}.\end{equation} \tag{ 48 }$$

  The computation of the normalization factor μ(β) may seem problematic, but it can be obtained from the dynamical trajectory, in the following way: for each ${g}_{\hat{A}}$, we define an 'unnormalized' function

  $$\begin{equation}{\tilde {g}}_{\hat{A}}\left(\eta ,\xi \right)=\frac{{g}_{\hat{A}}\left(\eta ,\xi \right)}{\mu \left(\beta \right)}={\text{e}}^{\beta H\left(\eta ,\xi \right)}\enspace \mathrm{Tr}\left[{\text{e}}^{-\beta \hat{H}\left(\xi \right)}\hat{A}\left(\xi \right)\right],\end{equation} \tag{ 49 }$$

  such that $\langle {g}_{\hat{A}}{\rangle }_{\text{CC}}\left(\beta \right)=\mu \left(\beta \right)\langle {\tilde {g}}_{\hat{A}}{\rangle }_{\text{CC}}\left(\beta \right)$.On the other hand, we know that for the identity operator, $\langle \hat{I}{\rangle }_{\text{HC}}\left(\beta \right)=1$, and therefore:

  $$\begin{equation}{\langle {g}_{\hat{I}}\rangle }_{\text{CC}}\left(\beta \right)=\mu \left(\beta \right){\langle {\tilde {g}}_{\hat{I}}\rangle }_{\text{CC}}\left(\beta \right)={\langle \hat{I}\rangle }_{\text{HC}}\left(\beta \right)=1.\end{equation} \tag{ 50 }$$

  Thus, we may compute μ(β) as $1/{\langle {\tilde {g}}_{\hat{I}}\rangle }_{\text{CC}}\left(\beta \right)$, and ${\langle {\tilde {g}}_{\hat{I}}\rangle }_{\text{CC}}\left(\beta \right)$ can be obtained from a dynamics propagation, i.e.:

  $$\begin{equation}\frac{1}{\mu \left(\beta \right)}=\underset{{t}_{f}\to \infty }{\mathrm{lim}}\enspace \frac{1}{{t}_{f}}{\int }_{0}^{{t}_{f}} \mathrm{d}t\enspace {g}_{\hat{I}}\left({\eta }^{\beta }\left(t\right),{\xi }^{\beta }\left(t\right)\right).\end{equation} \tag{ 51 }$$

  Summarizing, a final formula that permits to compute the hybrid ensemble averages is:

  $$\begin{equation}{\langle \hat{A}\rangle }_{\text{HC}}\left(\beta \right)=\frac{\underset{{t}_{f}\to \infty \enspace }{\mathrm{lim}\enspace }\underset{0}{\overset{{t}_{f}}{\int }} \mathrm{d}t\enspace {\text{e}}^{\beta H\left({\eta }^{\beta }\left(t\right),{\xi }^{\beta }\left(t\right)\right)}\enspace \mathrm{Tr}\left[{\text{e}}^{-\beta \hat{H}\left({\xi }^{\beta }\left(t\right)\right)}\hat{A}\left({\xi }^{\beta }\left(t\right)\right)\right]}{\underset{{t}_{f}\to \infty \enspace }{\mathrm{lim}\enspace }\underset{0}{\overset{{t}_{f}}{\int }} \mathrm{d}t\enspace {\text{e}}^{\beta H\left({\eta }^{\beta }\left(t\right),{\xi }^{\beta }\left(t\right)\right)}\enspace \mathrm{Tr}\left[{\text{e}}^{-\beta \hat{H}\left({\xi }^{\beta }\left(t\right)\right)}\right]}.\end{equation} \tag{ 52 }$$

  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.
• (b)  
  Equation (46) also holds if ${g}_{\hat{A}}$ is defined as:

  $$\begin{equation}{g}_{\hat{A}}\left(\eta ,\xi \right)=\lambda \left(\beta \right)\sum\limits _{\alpha }\delta \left(\eta -{\eta }_{\alpha }\left(\xi \right)\right){A}_{\alpha \alpha }\left(\xi \right),\end{equation} \tag{ 53 }$$

  $$\begin{equation}\lambda \left(\beta \right)=\frac{{Z}_{\text{CC}}\left(\beta \right)}{{Z}_{\text{HC}}\left(\beta \right)},\end{equation} \tag{ 54 }$$

  where η_α (ξ) are the adiabatic states:

  $$\begin{equation}\hat{H}\left(\xi \right)\vert {\eta }_{\alpha }\left(\xi \right)\rangle ={E}_{\alpha }\left(\xi \right)\vert {\eta }_{\alpha }\left(\xi \right)\rangle ,\end{equation} \tag{ 55 }$$

  and

  $$\begin{equation}{A}_{\alpha \alpha }\left(\xi \right)=\langle {\eta }_{\alpha }\left(\xi \right)\vert \hat{A}\left(\xi \right)\vert {\eta }_{\alpha }\left(\xi \right)\rangle .\end{equation} \tag{ 56 }$$

  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:

  $$\begin{equation}{\langle \hat{A}\rangle }_{\text{HC}}\left(\beta \right)=\frac{\underset{{t}_{f}\to \infty \enspace }{\mathrm{lim}\enspace }\underset{0}{\overset{{t}_{f}}{\int }} \mathrm{d}t\enspace \sum\limits _{\alpha }\delta \left({\eta }^{\beta }\left(t\right)-{\eta }_{\alpha }\left({\xi }^{\beta }\left(t\right)\right)\right){A}_{\alpha \alpha }\left({\xi }^{\beta }\left(t\right)\right)}{\underset{{t}_{f}\to \infty \enspace }{\mathrm{lim}\enspace }\underset{0}{\overset{{t}_{f}}{\int }} \mathrm{d}t\enspace \sum\limits _{\alpha }\delta \left({\eta }^{\beta }\left(t\right)-{\eta }_{\alpha }\left({\xi }^{\beta }\left(t\right)\right)\right)}.\end{equation} \tag{ 57 }$$

  This formula avoids the need to compute all the electronic excited states, necessary for the traces present in equation (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 straightforward average given in equation (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,

  $$\begin{equation}{\hat{\rho }}_{\text{HC}}\left(\xi \right)=\frac{1}{{Z}_{\text{HC}}\left(\beta \right)}{\text{e}}^{-\beta \hat{H}\left(\xi \right)},\end{equation} \tag{ 58 }$$

  in terms of its spectral decomposition for each ξ:

  $$\begin{equation}{\hat{\rho }}_{\text{HC}}\left(\xi \right)=\frac{1}{{Z}_{\text{HC}}\left(\beta \right)}\sum\limits _{\alpha }{\text{e}}^{-\beta {E}_{\alpha }\left(\xi \right)}{\hat{\eta }}_{\alpha }\left(\xi \right),\end{equation} \tag{ 59 }$$

  where E_α (ξ) are the eigenvalues, and ${\hat{\eta }}_{\alpha }\left(\xi \right)$ the projectors on the eigenspaces of the Hamiltonian $\hat{H}\left(\xi \right)$ (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 $\mathcal{M}$, as:

  $$\begin{equation}{\rho }_{\text{HC}}\left(\eta ,\xi \right)=\frac{1}{{Z}_{\text{HC}}\left(\beta \right)}\sum\limits _{\alpha }\delta \left(\eta -{\eta }_{\alpha }\left(\xi \right)\right){\text{e}}^{-\beta {E}_{\alpha }\left(\xi \right)}.\end{equation} \tag{ 60 }$$

  This distribution determines ${\hat{\rho }}_{\text{HC}}\left(\xi \right)$, since:

  $$\begin{equation}{\hat{\rho }}_{\text{HC}}\left(\xi \right)=\int \mathrm{d}\mu \left(\eta \right){\rho }_{\text{HC}}\left(\eta ,\xi \right)\frac{\vert \eta \rangle \langle \eta \vert }{\langle \eta \vert \eta \rangle }.\end{equation} \tag{ 61 }$$

  By comparing equation (60) with equation (40), it becomes clear that the error that this latter equation does is counting all possible states, whereas the true hybrid ensemble only counts the states in the adiabatic basis. Of course, one could choose a different basis, but the point is that the 'classical' Gibbs distribution (40) overcounts the quantum states. For a deeper discussion on this issue, we refer the reader to reference [5].In order to implement this procedure numerically, one should of course use some finite representation of the delta functions, giving them a non-zero width. It is unclear, however, that this would lead to an efficient scheme, since the propagation would probably have to be very long in order to obtain an accurate sampling of the quantum states.

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–Oppenheimer 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. Furthermore, 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:

$$\begin{equation}\hat{H}\left(Q,P\right)=\frac{{P}^{2}}{2m}\hat{\mathbb{I}}+\left(\begin{matrix}{V}_{0}\left(Q\right)& 0\\ 0& {V}_{1}\left(Q\right)\end{matrix}\right),\end{equation} \tag{ 76 }$$

where m is the dimer reduced mass, and the Morse potentials are given by:

$$\begin{equation}{V}_{i}\left(Q\right)={D}_{i}{\left(1-{\text{e}}^{-{b}_{i}\left(Q-{q}_{i}\right)}\right)}^{2}+{{\Delta}}_{i}.\end{equation} \tag{ 77 }$$

The parameters defining the Morse potential for the ith 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 ${\omega }_{i}=\sqrt{\frac{2{b}^{2}D}{m}}$. Figure 1 depicts these potential energy curves ⁶ .

Zoom In Zoom Out Reset image size

Using ground-state Born–Oppenheimer MD means that the classical degrees of freedom follow the Hamiltonian system defined by the function:

$$\begin{equation}\mathcal{H}\left(Q,P\right)={E}_{0}\left(Q,P\right)=\frac{{P}^{2}}{2m}+{V}_{0}\left(Q\right).\end{equation} \tag{ 78 }$$

The system is coupled to a Langevin thermostat, at the temperature given by $T=\frac{1}{{k}_{\text{B}}\beta }$. This dynamics provides an ergodic curve over the classical phase space with a visitation weight given by the Boltzmann factor ${\text{e}}^{-\beta {E}_{0}\left(Q,P\right)}$. Using the corrected averaging procedure defined by equation (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: $\hat{A}\left(Q,P\right)=Q\hat{I}$.

Figure 2 shows the results. As the model is particularly simple, we can display both the exact values (i.e. the hybrid canonical ensemble averages computed by performing the direct integration in phase space, using equation (8)), and the values produced by using the time-averages over the dynamics. In this latter case, we display both the corrected averages, that result from formula (72), and the ergodic averages using ground-state MD, which corresponds to the purely-classical canonical ensemble average. It is clear how the proposed reweighting formula yields the correct numbers.

Zoom In Zoom Out Reset image size

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.

Finally, the blue crosses in figure 2 (marked as MD T') are the results obtained by performing a single trajectory at a fixed temperature T', and then correcting appropriately for each different temperature. We have chosen T' ∼ 1.35ω₀/k_B, corresponding to the sixth data point in the figure, as it is in the middle of the temperature range that we want to probe. The accuracy seems to be slightly better than the results obtained by performing different trajectories at different temperatures. This should not be taken as a general rule: we have performed convergence analysis with both methods and we cannot conclude that the use of the different-temperature method always leads to more precise averages for equal total propagation time. Note, however, than the single-trajectory method may be more advantageous as the computation time saved by not computing the multiple trajectories at different temperatures can be used to converge better the single trajectory with respect to its total propagation time.