Advanced Sampling Methods for Quantum Free Energy Calculations

Advanced Sampling Methods for

Quantum Free Energy Calculations

By

ROBERTOVERDEL ARANDA

Supervisor: Dr. Sara Bonella

A dissertation submitted to the University of Amsterdam, the University La Sapienza of Rome, and ENS de Lyon & UCB Lyon 1 to obtained the triple diploma Erasmus Mundus AtoSiM 2.0 Master of Science.

A

BSTRACT

I

n this thesis, we address the problem of formulating an accurate and efficient scheme to compute quantum free energies. The inclusion of quantum effects in free energy calculations is important in systems that, presenting the problem of thermally activated events, also contain light atoms and/or are at low temperatures. Therefore, this topic is crucial in several branches of physics, chemistry and biochemistry. This work is divided in two parts. In the first part, a detailed review of two advanced computational methods is made. The first approach is the Quantum Single Sweep method [Bonella and Vuilleumier, (in press)], which is based on the path integral formalism. Though very precise in the reconstruction of quantum free energy profiles, this scheme is computationally expensive. The second technique that is revised here is the Quantum Thermal Bath [Dammak et al., Phys. Rev. Lett., 103:190601, (2009)], an approach that uses a generalized Langevin equation to effectively include quantum effects in the form a random force and a dissipative term. This technique has the advantage that is less costly than path-integral-based methods; however, it only gives exact results for harmonic systems. Furthermore, the method’s applicability for free energy calculations has never been tested systematically.

The second part of the present project studies precisely a route to combine the former two methodologies in a way that keeps the best elements of them, with the goal to generate an algorithm that is less computationally demanding than the quantum single sweep method but without the problems of the quantum thermal bath. While carrying out this program, we found that a direct use of the quantum thermal bath in the dynamical equations of the quantum single sweep does not yield the correct free energy landscape. Instead, the classical result is obtained for one-dimensional systems. These problems are carefully analysed and explained in this thesis. A second outcome of the present work is a novel scheme that is proposed to overcome the issues mentioned above. This new method combines some steps of the quantum single sweep with a hybrid Monte Carlo calculation in which the trial moves are proposed by evolving a molecular dynamics trajectory based on the quantum thermal bath. Although the performance of this novel technique is not numerically tested here, we present its theoretical formulation and the main points to push forward in its implementation, a task that we believe would be very interesting to follow, since to our knowledge it has not yet been done elsewhere.

A

CKNOWLEDGEMENTS AND

D

EDICATION

I

would like to thank my supervisor, Dr. Sara Bonella, for her advice and patience during the realization of this thesis. I learned a lot working with Dr. Bonella and I greatly appreciate her clarity and enthusiasm during our discussions. I also benefited from the meetings with some members of the group of Low-dimensional Oxides at INSP in Paris. Therefore, I wish to acknowledge Dr. Philippe Depondt, Dr. Fabio Finocchi, and Dr. Simon Huppert, for the useful discussions and remarks on my work and for providing me with the code to perform the quantum thermal bath calculations.

I am very grateful to Alexandra Guilleminot at ENS de Lyon, and Guy Byron at CECAM-EPFL, for their help in dealing with all the administrative formalities required for the beginning of my internship at CECAM-EPFL.

I also want to thank my AtoSiM classmates for their friendship and support, not only in the last five months but in the last two years. Thanks also to Przemysław Juda, my office (and hiking) mate in Lausanne, for the interesting discussions during the coffee breaks.

I acknowledge the European Commission and EPFL for awarding me an Erasmus Mundus grant and a scholarship within the Swiss-European Mobility Programme, respectively.

Finalmente, quiero agradecer de forma infinita a mi familia por todo el amor y el apoyo que siempre me han brindado, ingredientes imprescindibles para haber llegado hasta aquí.

T

ABLE OF

C

ONTENTS

Page

1 Introduction 1

1.1 The Classical Free Energy . . . 1

1.2 The Problem of Rare Events . . . 4

1.3 Nuclear Quantum Effects . . . 6

2 The Quantum Single Sweep Method 9 2.1 The Path Integral Representation of the Quantum Free Energy . . . 9

2.2 Temperature Accelerated Molecular Dynamics . . . 14

2.3 Free Energy Reconstruction . . . 18

2.4 Comments on the Accuracy and Efficiency of the QSS Method . . . 19

2.5 Other Methods to Compute Quantum Free Energies . . . 24

3 The Quantum Thermal Bath 27 3.1 Theoretical Aspects of QTB . . . 27

3.2 Applications of QTB in 1D Systems . . . 32

3.3 Drawbacks and Limitations . . . 35

3.4 A Further Approach: QTB with Path Integrals. . . 36

4 A New Hybrid Methodology: TAMD with QTB-MC 39 4.1 Direct Coupling of the TAMD System to a QTB . . . 40

4.2 Analysis of the Results Obtained using the Direct TAMD-QTB Methodology . . . 41

4.3 A New Approach: TAMD with QTB-MC . . . 45

5 Conclusions and Perspectives 51

A Appendix A: Exact Results for the Quantum Harmonic Oscillator in 1D 53

C

H A P T E R

1

I

NTRODUCTION

T

his chapter presents the necessary background to understand the problem that is tackled in this thesis, namely the problem of computing free energies in systems in which nuclear quantum effects are important. In Section 1.1, we revise the concept of classical free energy. In particular, a definition of the free energy with respect to some collective variables that characterize the progress of a reaction of interest is presented. Section1.2explains why standard molecular simulation techniques cannot be simply applied to compute free energies. This is known as the problem of rare events, and some sampling schemes that allow for an efficient computation of classical free energies are briefly summarized here. Finally, in Section1.3the importance of nuclear quantum effects is underlined, giving the motivation to study quantum free energies in which advanced sampling methodologies have to be employed as well, in order to overcome the rare event issue.

1.1

The Classical Free Energy

The free energy is a key concept in statistical mechanics and related fields. In general, the free energy is defined in the different thermodynamic ensembles as the logarithm the corresponding partition function multiplied by the factor −kBT, where kB is Boltzmann’s constant1 and T is

the temperature of the system. For instance, in the canonical ensemble, we have the following relation that gives the statistical mechanical definition of the Helmholtz free energy:

(1.1) _{F = −k}BT ln Q(N, V , T),

where Q(N, V , T) is the canonical partition function and it is defined as

1_{In S.I. units k}

(1.2) _{Q(N, V , T) = C} Z

dqdp e−H(q,p)/kBT_,

where C is a normalization constant, and H(q, p) is the Hamiltonian of the system of interest that contains N particles and occupies a volume V . q denotes the set of all particle positions and p is the vector containing the associated momenta. The point q belongs to the so-called configuration space, which is defined by the containing volume V and it is denoted D(V ) ⊂ RdN for a d-dimensional system.2The integral in (1.2) is carried out over the phase space, which is the space of all possible microscopic configurations (q, p). Equation (1.1) is crucial, for it serves as a generator of other thermodynamic quantities via differentiation [1,2].

The free energy is also the key quantity for calculation of processes that involve thermally activated events [1–4]. In an activated event the system under consideration must overcome a free energy barrier for the process to go forward. The starting and ending points of the transformation involved in this sort of processes are in most cases long lived states of the system. Such states are sometimes called metastable. They are in fact associated with the minima of the free energy surface, and an activated event can therefore be described as a transition between two of these minima (see Figure1.1).

The free energy surface also gives information about the relative stability of the long lived states of the system. This can be obtained by computing the free energy difference of the corre-sponding system configurations. Moreover, the location and height of the barriers separating the minima of the free energy landscape define transition states and transition rates involved in an activated process. In elementary chemical reaction theory, the transition rate k of a given reaction follows Arrhenius equation [5]:

(1.3) _{k = Ae}−∆F/kBT_,

where A is the so-called ‘frequency factor’ and∆F, which is normally referred to as the ‘activation energy’, is the height of the barrier that separates ‘products’ and ‘reactants’ in the free energy landscape. Relation (1.3) will helps us to understand the problem behind free energy calculations, as discussed in the next section.

Apart from the case of low-dimensional model systems, the free energy landscape of more interesting systems found in applications in physics, chemistry and biology, is in general a very complex mathematical object that typically contains a huge amount of saddle points (transition states) and local minima (metastable states), as illustrated in Figure1.2. However, most of these metastable basins are irrelevant when studying a given reaction or transition. Moreover, as in any other sampling problem, we would like to spend as much time as possible sampling the

2_{In the following, it will be implicitly assumed that all the integrals done with respect to the coordinates are}

1.1. THE CLASSICAL FREE ENERGY

Figure 1.1: Schematic representation of a free energy surface of a molecular system. The free energy is a function of a reaction coordinate (collective variable) that indexes the reaction of interest. Two metastable states labelled ‘A’ and ‘B’ are separated by a free energy barrier defining a transition state. (Reproduced from [6].)

rather relevant regions that give the main contribution in the calculation of ensemble averages. This can be somewhat achieved by using a set of generalized coordinates which characterize the progress of the reaction under study. Such generalized coordinates are called collective variables.3 Some examples are interatomic distances, dihedral angles, dipole moments, etc.

Let us analyse the utility of working in collective variables. To this end, let us consider a set ofνcollective variables denoted byθ_{(q) = {}θ1(q), . . . ,θν(q)}, where we have assumed that the collective variables are functions only of the coordinates q. One can now define a free energy associated withθ_{(q). If z = {z}1, . . . , zν}is a given realization ofθ(q), then the associated free energy is given by [1,7]

(1.4) _{F(z) = −kB}T ln P(z),

with P(z) being the marginal probability of θ, defined as

(1.5) _{P(z) =} C Q(N, V , T) Z dqdp e−H(q,p)/kBT ν Y α=1δ (_{θα(q) − zα}),

where C is the same constant that appears in Eqn. (1.2), and_δdenotes the Dirac delta function. Equation (1.5) can be further simplified if the Hamiltonian of the system is separable, i.e., if it can be written as H(q, p) = K(p) +U(q), where K(p) is the kinetic energy of the system and U(q) is the potential energy, since then the integral with respect to p as well as the prefactor C drop out, and (1.5) reduces to

Figure 1.2: Example of a complex free energy landscape with a large number of critical points (metastable states, transitions states, etc.). (Reproduced from [8].)

(1.6) _{P(z) =} 1 Z(N, V , T) Z dq e−U(q)/kBT ν Y α=1δ (_{θα(q) − zα}),

where Z(N, V , T) =R dq e−U(q)/kBT _{is the configurational partition function.}

The free energy defined in Eqn. (1.4)4 encodes the statistical weights of the states connected through the transitions that are indexed by the collective variables. Equation (1.4) also constitutes the starting point for several techniques to compute classical free energies as discussed in Section

1.2. Lastly, note that the discussion presented here has been limited to classical systems. This is for an illustrative purpose and we shall define a quantum version of (1.4) in Chapter2, where a method to compute quantum free energies is presented.

1.2

The Problem of Rare Events

Equation (1.4) can be used to devise a direct approach to reconstruct F(z) in a computer simulation. Indeed, one only needs to evolve a trajectory of the system x(t) = (q(t), p(t)), using any kind of dynamics provided it is ergodic with respect to the canonical ensemble5 [4, 9], e.g. Langevin dynamics (1.7)    ˙ p(t) _{= −}γ_{p(t) −}∂U(q)_{∂q +}p2mγ η(t), ˙ q(t) = p(t)_m ,

4_{This free energy is sometimes called the Landau free energy.}

5_{In this thesis, we shall restrict ourselves to the canonical ensemble. However, most of the concepts and techniques}

1.2. THE PROBLEM OF RARE EVENTS

where m is the mass of the particles,6_γis the friction coefficient, and_η(t) is a Gaussian white noise with mean zero and covariance satisfying:

(1.8) _〈η(t)_η(t0_{)〉 = k}BTδ(t − t0).

Note that the amplitude of the random forceη(t) is tuned such that the equipartition theorem is verified [10].

From the evolution of (1.7) one can build a histogram of the values of_θ(q(t)) visited during the course of the simulation. Such histogram is used to approximate P(z), from which F(z) is obtained via (1.4). In spite of its simplicity, this rather naive algorithm is not useful in practice whenever the system under study has a rough free energy landscape with deep minima where it can get trapped. To understand the breakdown of this approach, let us consider a system with a free energy surface such as the one depicted in Figure1.2, and let us call∆F the height of the highest relevant barrier. According to Eqn. (1.3), the time needed to overcome such a barrier is of the order of

(1.9) _{t ∼ O (e}∆F/kBT_).

In order to make an accurate sampling of F(z) the timescale used in our simulations must be much larger than (1.9). However, if ∆F is large compared to kBT this timescale becomes

prohibitively long. This is the problem of rare events.

There exists a vast number of methods to deal with this issue for classical systems. Here some of the most popular schemes are summarized (for a more detailed description of these and some other techniques the reader is referred to [1–4]).

• Thermodynamic perturbation: In this method the free energy difference between two states A and B, is expressed as a standard canonical average, namely

(1.10) ∆F_{AB= −k}BT ln

D

e−(UB−UA)/kBTE

A,

where the average on the right-hand side is taken with respect to the equilibrium distri-bution of system A. Additionally, when working in collective variables, usual sampling techniques relying on histogram methods can also be employed [4].

• Umbrella sampling: The basic idea of this method is that in order to get an accurate estimate of the free energy difference between two states or systems A and B, the relevant regions of configuration space that are accessible to both systems must be sampled simultaneously.

This is achieved by replacing the Boltzmann factor that is used in a standard Monte Carlo (MC) simulation by a nonnegative, but otherwise arbitrary, weight functionπ(q). The key equation of this method is

(1.11) e−∆FAB₌D_e−(UB−UA)/kBTE A= D e−UB/kBT_/πE π D e−UA/kBT_/π E π .

• Adiabatic molecular dynamics: It accelerates the sampling in the collective variables space by making use of an artificially high temperature in the equations of motion. The dynamical equations are then reformulated in a new coordinate system such that the collective variables themselves are explicit coordinates. As the name of the method suggests, this technique relies on the assumption of adiabatic separation between the collective variables and the physical degrees of freedom. This crucial assumption allows to obtain the correct free energy profile, even though the thermodynamics of the system is, strictly speaking, somehow affected by the fact that there are two temperatures.

• Metadynamics: It belongs to the class of methods that use a Lagrangian formulation in an extended phase space. In this case, the extended system includes the set of collective variables. This is done in order to obtain a sampler in the collective variables space. The sampling is then accelerated by introducing a bias in the form of a history-dependent potential that allows to keep the system out of the regions of phase space that have already been visited. Such a potential is constructed as a sum of Gaussians centered along the trajectory followed by the set of collective variables.

In Section 2.2 another method to efficiently explore and reconstruct classical free energies is discussed in more detail. This method bears a lot of similarities with adiabatic molecular dynamics and metadynamics, as it will be shown the aforementioned section.

1.3

Nuclear Quantum Effects

In classical molecular simulations the nuclei of the atoms that form a molecular system are regarded as classical particles in the presence of either an effective potential, which is obtained through quantum mechanical calculations, or an empirical potential fitted so as to agree with experimental data or some other quantum mechanical computations [11]. There is, nevertheless, a broad variety of systems in which nuclear quantum effects (NQEs) are expected to play an important role in the physical description of the system.

Some of the NQEs that might be important for the proper description of a molecular system are the zero-point energy (ZPE), tunneling, and exchange effects. The ZPE is perhaps one of first quantum effects studied in an undergraduate course of quantum mechanics. In fact, this concept

1.3. NUCLEAR QUANTUM EFFECTS

first appears in the solution of very simple one-dimensional problems such as the quantum harmonic oscillator and the square-well potential [1,12,13]. ZPE essentially refers to the ground state energy, that is, the lowest possible energy of a bound state potential like the ones mentioned above. It is basically a consequence of Heisenberg’s uncertainty principle: If a particle is at bottom of the potential, and therefore very localized in space, the uncertainty relation tells us that it will acquire a finite momentum, which in turn leads to a minimum finite kinetic energy. For the simple systems mentioned before, it is even possible to show that an accurate evaluation of the momentum uncertainty of a localized particle yields a ZPE that is exactly equal to the ground state energy (see, for example, Problem 4.2 in [13]). Thus, ZPE is the minimum motion required for a particle that is localized or confined in space, and it will always be higher than the minimum value of the potential energy, contrary to what happens in classical mechanics in which the particle can be in the minimum of the potential with zero kinetic energy.

Tunneling is encountered in the analysis of systems with a potential barrier. This is why this quantum effect is so relevant for free energy calculations of quantum systems. Tunneling is a purely quantum mechanical effect that has no classical analog. It refers to the finite probability that a quantum particle has to penetrate in or pass through a finite potential barrier, being this a region with an energy larger than the particle energy. This NQE has far-reaching applications in several branches of modern physics such as nuclear and semiconductor physics [13].

Exchange effects are also purely quantum mechanical and come from the symmetry of the many-body wave function. If we imagine a system composed by identical particles, that is, particles that are exact replicas of each other in the sense that there is no experiment capable of detecting any fundamental difference among them, from a classical point of view one is able to tell the difference between physical situations related by the exchange of any pair of these identical particles. In the classical world this is possible because identical particles can be distinguished by following their nonidentical trajectories. On the other hand, in quantum theory the concept of a continuous trajectory for a particle makes no sense, and hence, such physical situations must be treated as one. The consequences of this restriction leads to the spin statistics theorem, which basically states that particles in nature can be either in a symmetric state or an antisymmetric state.7The particles belonging to first group are called bosons, whereas the particles in the second category are fermions. Spin statistics has allowed us to discover very fascinating effects such as the Bose-Einstein condensation [1], in bosonic systems, or to talk about the Fermi surface that is a key concept in solid state physics [14].

In general, NQEs must be taken into account whenever we are treating systems that contain light nuclei and/or at low temperatures [1,15,16]. Accordingly, hydrogen-containing systems have become a major point of interest of studies concerning NQEs in the recent years (see, for example, [1,16–18]). Of particular interest, specially from a biological perspective, are the investigations of the role played by NQEs in determining the static and dynamic properties of

7_{The quantum states are said to be symmetric or antisymmetric depending on whether the many-body wave}

water [19,20]. However, NQEs studies are not limited to hydrogen-containing systems. Other systems studied in the context of NQEs are liquid helium (for this particular system ZPE is crucial to explain what prevents helium from freezing at very low temperatures [13]), lithium and related compounds, and carbon-based crystals, just to mention a few examples (see [21,22] and references therein).

This set of examples shows us the need to develop methods to compute quantum free energies. Such methods must be such that they can deal with the problem of thermal activation discussed in the last section, while at the same time accounting for NQEs. In the next chapter, we first define a quantum version of the free energy (1.4) and then present a state-of-the-art method for accurately sampling quantum free energies.

C

H A P T E R

2

T

HE

Q

UANTUM

S

INGLE

S

WEEP

M

ETHOD

I

n spite of all the work that there exists to overcome the problem of rare events in classical systems (see Section1.2), there is much less work regarding quantum mechanical systems, i.e., systems in which ZPE, tunneling, and exchange effects must be rigorously included. In this chapter, a method belonging to the state-of-the-art of quantum free energy calculations is presented: The quantum single sweep (QSS) method [23]. Section2.1gives a detailed mathemati-cal formulation of the starting point of the QSS method, namely, the path integral representation of the quantum free energy. QSS is an extension of a classical method to sample and reconstruct free energies, which is the single sweep method [24]. Hence, in Sections2.2and2.3the main ingredients of the latter method are explained, but in a context that adapts the derivations to the quantum case. Section2.4finishes the discussion on QSS by analysing the accuracy and efficiency of this method. Finally, some other methodologies to calculate quantum free energies are briefly summarized in Section2.5.

2.1

The Path Integral Representation of the Quantum Free

Energy

As most of the existing methods to compute quantum free energies, QSS is based on the Feynman path integral formulation of quantum statistical mechanics [25]. There are two main reasons for the use of this formalism. In the first place, the Feynman path integral gives an intuitive physical picture of quantum processes, which can be adapted to the classical limit in a natural fashion [1,21]. Secondly, the path integral formulation of equilibrium statistical mechanics is ideally suited for simulation. Indeed, the so-called classical isomorphism [26] allows to map the quantum system onto an analogous classical polymer system, meaning that one can make use

of all the tools already developed for simulating such kind of classical systems. This analogy is further explained later on in this section.

Let us start by defining the quantum free energy of a system of N distinguishable1 particles. This is done in analogy with the classical definition given in (1.4). To this end, we shall assume that the collective variables are functions of the coordinate operators of the system ˆx, i.e.,

θ( ˆx) = {θ1( ˆx), . . . ,θν( ˆx)}. We have (2.1) _{F(z) = −}β−1ln Tr " e−β ˆH Q δ(θ( ˆx) − ˆIz) # ,

whereβ_{= 1/(k}BT), ˆH is the Hamiltonian of the quantum system, Q = Tr[e−β ˆH] is the quantum

canonical partition function,δ(θ_{( ˆx) − ˆIz) ≡}Qν

α=1δ(θα( ˆx) − ˆIzα), ˆI is the identity operator, and Tr denotes the trace of an operator. The quantity ˆ_ρ( ˆ_{H) := e}−β ˆH_{/Q is called the canonical density}

matrix or quantum Boltzmann density.

Despite the similar form of Eqn. (2.1) and its classical counterpart (see Eqn. (1.4) and (1.5)), there are three major differences between these two expressions. First of all, Eqn. (2.1) involves the use of operators (indicated by a hat over the corresponding quantities) that are associated with observables (Hamiltonian and collective variables). Secondly, the expression for the quantum free energy is more abstract than that of the classical free energy, since the form of the operators involved in the former quantity is not completely specified until we select a representation in a given basis. Finally, the expectation values in the quantum case are expressed as traces, which only take the form of an integral if a continuous representation is chosen.

Since the trace of an operator does not depend on the choice of the basis set,2 we can use any representation we like to manipulate Eqn. (2.1). Here we will use the position representation, which is a continuous one and turns out to be very convenient to address the problem at hand. Thus, using the coordinate basis, we get

(2.2) _{F(z) = −β}−1ln Q−1 Z

dx_{1 〈x1|e}−β ˆH_δ(_{θ( ˆx) − ˆIz)|x1〉.}

In the equation above, and in what follows, we use one dimensional notation for the coordinates of the system. Nonetheless, the generalization to more dimensions is straightforward, at least in the case of distinguishable particles. Now let us multiply by the unit operator ˆ_{I =}R dx|x〉〈x|, and use the fact that the collective variables are diagonal in the coordinate basis, i.e., 〈x|δ(θ_{( ˆx) − ˆIz)|x〉 =}

δ(θ_{(x) − z). We have}

1_{The assumption of distinguishability will be discussed at the end of the this section.} 2_{Proof : Tr[ ˆA}0_{] = Tr[ ˆU ˆA ˆU}†

] = Tr[ ˆU†U ˆA] = Tr[ ˆA], where ˆˆ U represents a unitary transformation and we use the cyclic property of the trace in the second step.

2.1. THE PATH INTEGRAL REPRESENTATION OF THE QUANTUM FREE ENERGY F(z) = −β−1_{ln Q}−1Z _dx 1dx2 〈x1|e−β ˆH|x2〉〈x2|δ(θ( ˆx) − ˆIz)|x1〉 = −β−1ln Q−1 Z dx1dx2 〈x1|e−β ˆH|x2〉δ(θ(x1) − z)δ(x1− x2) = −β−1ln Q−1 Z dx1 〈x1|e−β ˆH|x1〉δ(θ(x1) − z). (2.3)

In the second line we used the fact that the position vectors are orthonormal, i.e., 〈x1|x2〉 =

δ(x1− x2). Now our problem is to evaluate the matrix element appearing in the last line of the

former equation. In general, it is very difficult —if not impossible— to determine the effect of e−β ˆH_{upon the position vector |x〉. Indeed, the complexity of this problem is equivalent to that of}

solving the Schrödinger equation. The problem arises from the fact that the Hamiltonian ˆH is the sum of a kinetic energy term ˆK and a potential energy part ˆU, and these two operators do not generally commute with each other, i.e., [ ˆK , ˆ_{U] 6= 0. Consequently}

(2.4) e−β ˆH_{= e}−β( ˆK + ˆU)_{6= e}−β ˆKe−β ˆU.

If this separation would be true matters would simplify considerably. However, we can overcome this difficulty by using the same trick that is employed when studying the action of the classical time evolution operator upon a phase space vector (see, for instance, Section 3.10 of [1]). The trick consists of using the Trotter theorem, which states the following:

(2.5) e−β( ˆK + ˆU)_{= lim}

P→∞

h

e−β ˆU/2Pe−β ˆK /Pe−β ˆU/2PiP.

This is the first step in deriving the path integral representation of the quantum free energy. The next step is to insert P − 1 realizations of the identity operator into the matrix element that we wish to compute, that is,

〈x1|e−β ˆH|x1〉 = lim P→∞〈x1| h e−β ˆU/2Pe−β ˆK /Pe−β ˆU/2PiP|x1〉 = lim P→∞ Z dx2· · · dxP 〈x1|e−β ˆU/2Pe−β ˆK /Pe−β ˆU/2P|x2〉 (2.6)

× 〈x2|e−β ˆU/2Pe−β ˆK /Pe−β ˆU/2P|x3〉 · · · 〈xP|e−β ˆU/2Pe−β ˆK /Pe−β ˆU/2P|x1〉.

Note that in the second line of (2.6) there are now P matrix elements to be evaluated. The main advantage of using the Trotter theorem is that each of these matrix elements can be computed in a closed form. Indeed, let us consider one of these matrix elements

where we used the fact that ˆU depends only on the positions of the particles, and hence it is diagonal in the position representation. On the other hand, the kinetic energy operator ˆ_{K = ˆp}2/2m, m being the mass of the particles, is not diagonal in the coordinate basis set, but in the momentum representation. Hence, it is convenient to switch to the latter representation, namely

〈xi|e−β ˆK /P|xi+1〉 = Z dp dp0 _{〈xi|p〉〈p|e}−β ˆp2/2mP_|p0_〉〈p0_|xi+1〉 = 1 2π~ Z dp ei p(xi−xi+1)/~_e−βp2/2mP_, (2.8)

where we have used that 〈x|p〉 = (2π~)−1/2ei p·x/~ _{and 〈p|p}0_{〉 =δ(p − p}0_{), with~ being Planck’s}

constant3divided by 2_π. The integral in (2.8) can be performed analytically as it takes the form of a Gaussian integral upon completing the square in the argument of the exponential, yielding the following result

(2.9) 〈xi|e−β ˆK /P|xi+1〉 = µ _mP 2πβ~2 ¶1/2 exp · − mP 2β~2(xi− xi+1) 2¸_.

Combining Eqn. (2.9), (2.7) and (2.6) with (2.3) leads to the following equation for the free energy:

(2.10) _{F(z) = −}β−1lnQ−1 Z dx1· · · dxP exp " − P X i=1 mP 2_β~2(xi− xi+1) 2 −β P P X i=1 U(xi) # δ(θ(x1) − z),

whereQ is the path integral representation of the canonical partition function Q (except for a multiplicative factor) (2.11) Q = Z dx1· · · dxP exp " − P X i=1 mP 2β~2(xi− xi+1) 2 −β P P X i=1 U(xi) # .

where x_{P+1= x}1in Eqn. (2.10) and (2.11). Also, note that these equations have been written for

finite P, and only become exact for P → ∞. In fact, the leading error in (2.6) for finite P is of the orderO((_β/P)3) [1].

Next, let us define an effective potential as follows

(2.12) U_eff(x₁, . . . , x_P;_{β) :=} "_P X i=1 1 2mω 2 P(xi+1− xi)2+ 1 P P X i=1 U(x_i) # , withωP= p

P /β~and x_{P+1= x}1. Using this definition, we arrive at the very convenient final form

of the quantum free energy in the path integral representation, namely

2.1. THE PATH INTEGRAL REPRESENTATION OF THE QUANTUM FREE ENERGY

(2.13) _{F(z) = −β}−1lnQ−1 Z

dx_{1· · · dxP} e−βUeff(x1,...,xP;β)δ₍θ_(x

1) − z).

It is important to stress that Eq. (2.13_{) becomes exact in the limit P → ∞. For finite P, the integral}

in the expression above is completely equivalent to the configuration integral of a closed ring polymer (the fact that the polymer is a closed ring is due to the condition x_P+1= x1). The ‘beads’ of

the polymer, which would be described by the coordinates x_i, are linked together through nearest neighbour harmonic couplings in the effective potential (2.12), and each one of them interacts with the P-th fraction of the potential U, which describes an external and/or interaction potential. Such harmonic terms represent the quantum kinetic energy of the system, as they come from the evaluation of the kinetic part (see Eqn. (2.9)). This is the so-called classical isomorphism mentioned at the beginning of this section.

When looking at Eqn. (2.13), it seems that the coordinate x₁ plays a special role, for it is the only coordinate that appears in the delta function that constrains the collective variables at the value z. However, this is just a consequence of our choice in the initial representation used to write Eqn. (2.2) and there is really nothing particularly special about such a coordinate. In fact, one can substitute the delta function in (2.13) by its average with respect to the coordinates of the beads xi, i.e.,PP_i=1δ(θ(xi)− z)/P. This is possible since the path integral is invariant under a cyclic

reassignment of all the coordinates xi→ xi+M, with M being an integer, and taking into account

the wraparound condition x_{P+M= x}M [1]. Nevertheless, it turns out that making a distinction

between the coordinate x1and the rest of the coordinates is more convenient from a numerical

point of view than treating all the coordinates xion equal footing in the QSS method. This will

be discussed in Section2.2.

Let us now discuss one of the starting assumptions in the developments presented here, namely, the fact that the quantum particles are regarded as distinguishable. This assumption, that in the case of an ideal gas of quantum particles yields the Boltzmann distribution (e.g. see Section 11.3 in [1]), consists of neglecting the spin statistics, or in other words, the exchange effects4that were discussed in Section1.3. Exchange effects involve long-range correlations of delocalized wave functions [1], and therefore, they can usually be ignored for the systems that we are interested in, such as protons, which have a highly localized wave function, at least at not too low temperatures. This localization can be quantified, for example, by comparing the thermal wavelength5 with the typical (average) positions of the particles in the system. Therefore, the path integral treatments that include the spin statistics are not discussed in the present work, however, they can be found in the literature (see, for example, [27] and [28]).

4_{In the context of path integrals, exchange effects involve endpoint conditions that ‘tie’ the paths of the different}

particles together at the endpoints.

5_{The thermal wavelength is defined asλ}

th≡ h/p2πmkBT . This quantity allows us to determine the conditions

under which quantum effects become relevant. For instance, if the particle concentration N/V of an ideal gas is such that (N/V )1/3≤ λththen the gas will obey a quantum statistics. Otherwise, the Boltzmann distribution serves as a

Figure 2.1: Schematic illustration of two quantum Boltzmann particles represented as closed ring polymers. The beads of each of these polymers are linked together through harmonic couplings, whereas the interaction between the quantum particles occurs via the interaction of beads with the same index only. In this figure, P = 8. (Reproduced from [1].)

Having justified the focus on Boltzmann particles, let us write down explicitly the path integral that appears in (2.10) and (2.11) for N of these particles in d dimensions

(2.14) Z N Y i=1 dr(1)_i · · · dr(P)_i exp ½ −β P X j=1 · N X i=1 1 2miω 2 P ³ r( j)_i − r( j+1)_i ´2+1 PU ³ r( j)₁ , . . . , r( j)_N´ ¸¾ .

To avoid confusion, let us remark that each quantum particle is represented by a closed ring polymer. Therefore, in the path integral representation a system of N quantum particles is mapped onto a system of N classical cyclic polymers. An interesting aspect in this case is that the interaction between the quantum particles are now represented as interactions between beads with the same index j of the corresponding ring polymers (see the potential energy part in (2.14)). This is schematically represented in Figure2.1. In the following, we stick to the one dimensional case for notational convenience.

2.2

Temperature Accelerated Molecular Dynamics

Equation (2.13) can be use as starting point to formulate different methodologies to compute quantum free energies. In effect, since (2.13) can be associated with the free energy of an isomorphic classical system, one can apply basically any method to compute classical free energies (see Section1.2) to explore and reconstruct this free energy. The QSS method, introduced by Bonella and Vuilleumier [23] (see also the Appendix in [16]), uses for this purpose the single sweep method developed originally by Maragliano and Vanden-Eijnden [24]. The single sweep method is a two-step methodology. The first step, which is called temperature accelerated molecular dynamics (TAMD) [9], rapidly explores the important regions of the free energy landscape associated with a set of collective variables, and allows to get local estimates of the gradient of

2.2. TEMPERATURE ACCELERATED MOLECULAR DYNAMICS

the free energy function in such regions. The second step carries out the reconstruction of the free energy globally, by means of a variational procedure. In this section, the details of the exploration step are revised, whereas the reconstruction step is explained in Section2.3.

The basic idea of TAMD is to introduce an extended system where the collective variables are treated as dynamical ones; next, the sampling is accelerated by using an artificially high temperature for the collective variables. The extended system used in TAMD to explore (2.13) is described by the following equations of motion:

           m ¨x1= −∇x1Ueff(x;β) −κ ν X α=1 (_θα(x1) − zα)∇x1θα(x1) + Therm(β), m ¨x_{i= −∇xi}U_eff(x;_{β) + Therm(β), i = 2,..., P,} (2.15) M ¨z_α=κ(θ_α(x1) − zα) + Therm( ¯β), α= 1, . . . ,ν.

Here x denotes a vector containing all the coordinates of the polymer beads, i.e., x = {x1, . . . , xP}.

Therm(β) represents a thermostat at the temperature specified byβ. The nature of this ther-mostat is not important as long as we make sure that the canonical distribution is sampled [9]. Yet, further clarification is needed regarding the use of a different inverse temperature ¯βfor the auxiliary variables z. We shall touch upon this delicate point later on in this section. Note also that a different ‘mass’6 M is assigned to the auxiliary variables.

The dynamical system (2.15) corresponds to a potential energy, which not only includes the effective potential Ueff(x) given by Eqn. (2.12), but also a harmonic coupling betweenθ(x1) and

the auxiliary variables z, that is

(2.16) U_{κ(x, z) = U}eff(x) + 1 2κ ν X α=1 (θ_α(x1) − zα)2.

where the dimension of the parameterκis such that the right-end side of (2.16) has dimensions of an energy.

There are three parameters in the system (2.15) to be adjusted, namely ¯_β,_κ, and M. A right choice of ¯βallows the auxiliary variables to efficiently explore the interesting regions in collective variable space of the free energy landscape. An appropriate choice ofκand M is what actually makes this method correct from a physical viewpoint. Let us start by looking at the artificial mass M. We consider the limit behaviour of (2.15_{) when M À m.}7In this case, adiabatic separation between the physical variables x and the auxiliary variables z is induced. To be more precise x(t) will evolve at fixed z, on a timescale much shorter than the typical timescale of z(t). If this

6_{In fact, this mass-like parameter determines the timescale on which z evolves dynamically.}

7_{Actually one might have to choose properly the value of some of the thermostats parameters too. For instance, if}

Langevin thermostats are used in the TAMD system, which include a random force and a friction term (see Eqn.(1.7)), the friction coefficient associated with z, say ¯γ, must be much larger than the one in thermostat of the x variables, i.e.,

¯

condition is true then the physical variables thermalize on the conditional probability of x at given z:

(2.17) ρ_{κ(x|z) =} e

−βUκ(x,z)

R dx e−βUκ(x,z).

On the other hand, z(t) will evolve feeling the average value of x(t). This to say that z(t) will effectively evolve according to an equation that is the average of the third equation in (2.15) with respect toρ_{κ(x|z) [}9]. The only term affected by the average is the first one on the right-hand side of the aforementioned equation, that is

(2.18) _{− ∇}zFκ(z) = −κ Z

dx (z −θ(x1))ρκ(x|z), where the quantity F_κ(z) is defined as

(2.19) F_{κ(z) := −β}−1ln Z_κ−1 Z dx e−βUeff(x)_e−βκ2 Pν α=1(θα(x1)−zα)2_, with Z_κ=R dxdz e−βUκ(x,z)_{. F}

κ(z) is a ‘smoothed’ version of the target free energy F(z) given in (2.13). In fact, (2.13) is formally recovered from (2.19) in the limitβκ→ ∞. This is true as the Dirac delta function can be approximated by a Gaussian [29]

(2.20) δ_{(x) = lim} σ→0 1 p 2_πσ2e −x2/2σ2 .

In the present case_σ2_{= 1/(βκ}). From this point we can understand the role played by the coupling constantκ, namely, it controls how smooth F_κ(z) is with respect to F(z). In practical terms, we can adjust the value of 1/(βκ) depending on the level of resolution we want for F(z): Its square root should be smaller than the desired scale of resolution. This approximation is the reason why we took a harmonic coupling between the collective variables and the auxiliary ones in the first place.

Regarding the approximation of the delta function by a Gaussian, there is a comment that can be added. In Section2.1it was mentioned that in the expression for the free energy (2.13) we constrain the collective variables only through the coordinate x1instead of treating all the xi

coordinates on equal footing, and that this is done this way because of numerical convenience. The reason for this is that if we consider the average of the delta function,PP

i=1δ(θ(xi) − z)/P,

which is a sum of several delta functions, we cannot approximate such an average with a single Gaussian, unless the logarithm is taken first. However, this might yield possible singular forces that may be difficult to include in our simulations.

We have then proven that in the limit situations M À m andβκÀ 1, z(t) evolves approxi-mately according to the equation

2.2. TEMPERATURE ACCELERATED MOLECULAR DYNAMICS

(2.21) _{M ¨z = −∇zF(z) + Therm( ¯β}).

Let us stress that F(z) in (2.21) is the target free energy at the physical inverse temperatureβ. Now the equilibrium probability density of z(t), solution of the TAMD system (2.15), sayρ_κ,M,ξ(z), whereξrepresents all thermostat’s parameters that must be adjusted to guarantee adiabatic separation, is approximately the same as that of z(t), solution of (2.21), namely

(2.22) _{ρκ,M,ξ(z) ≈} e

− ¯βF(z)

R dz e− ¯βF(z).

The result above is what makes TAMD a useful tool to explore free energy landscapes in collective variables space. Indeed, all one needs to do is to monitor the evolution of z(t) in (2.15), build

ρκ,M,ξ(z) for example by histogramming, and then F(z) is obtained via (2.22). This approach is

efficient because (2.22) holds for whatever value of ¯_β, hence we can rise the artificial temperature ¯

T, or analogously decrease ¯β, so that ¯β∆F = O(1), in which case the timescale for the auxiliary variables to overcome the free energy barrier∆F is of the order O(e1) (see Eqn. (1.9)), accessible to numerical simulation. Moreover, it can be shown that this approximation can systematically be improved by increasing the values of M andβκ, as the error is of the order ofO(m/M) + O(βk) [9].

Let us finish this section by discussing two delicate points related to the fact that we use different temperatures for the physical and the auxiliary variables. In principle this might seem as a recipe for disaster, for the thermodynamics of the system of interest, i.e., the one described by the variables x at temperature T, will be generally affected by this trick. Indeed, it was mentioned that the nature of the thermostat used in (2.15) is irrelevant as long as it samples the canonical distribution. Now, that such distribution is sampled in this method is just a conjecture, since strictly speaking the system (2.15) is not at equilibrium due to the fact that it involves two temperatures. This conjecture is, nonetheless, expected to be true when adiabatic separation between x and z holds. Loosely speaking, if z are adiabatically separated from the physical degrees of freedom, the heat exchange between them occurs very slowly, and therefore, both x and z are, to a good approximation, canonical atβand ¯β(see equations (2.17) and (2.22)8), respectively. A different way to argue why a scheme such as TAMD, which introduces an artificially high temperature and makes use of adiabatic separation, gives the correct free energy can be adapted from the work of Rosso et al. [30], who developed a method that bears a lot of similarities with TAMD and that was mentioned in Section1.2, namely adiabatic free energy dynamics.

Finally, one may question the validity of adiabatic separation in the context of TAMD, as another crucial ingredient of the method consists of taking the limit_{κ→ ∞, this parameter being}

the one that couples_θ(x) and z. Unfortunately, to the best of our knowledge, there is no rigorous proof showing that both conditions can always be reconciled. One should, instead, try different values of the parameters in preliminary simulations to make sure that adiabatic separation holds while having a large value ofκ. Adiabatic separation can be monitored in a simulation by checking that 〈12m ˙x2〉 =

P

2kBT and 〈 1

2M ˙z2〉 =ν2kBT, to a given degree of precision.¯

2.3

Free Energy Reconstruction

In the previous section, we explained the way in which TAMD efficiently explores the free energy landscape in collective variable space. We also mentioned that we could even use a TAMD simulation to reconstruct the free energy profile by histogramming the collective variables chosen for a given problem. This approach for reconstructing the free energy profile is, however, not so efficient if the number of collective variablesνis large, for the number of bins of the histogram scales exponentially withν. In this section, we present a more efficient way to perform this task.

This reconstruction step consists in turn of three steps. First, we need to select points along a z-trajectory generated in a TAMD simulation.9We shall refer to this points as centers and denote them z_{k, with k = 1,..., K. The centers can be dropped along z(t) as follows: start by choosing} z1= z(t = 0), next place a new center zk when the trajectory z(t) has reached a point that is

separated from all previous centers at least by a distance dz. The parameter dz controls the

density of the covering by the centers: The smaller dz, the more centers are deposited and the

higher the precision of the reconstructed free energy profile.

In the second step, we compute estimates of the mean force locally at each of these centers z_k. The mean force is defined as minus the gradient of the free energy profile,10and can therefore be estimated from Eqn. (2.18). Invoking the ergodic hypothesis, we have

(2.23) fk= lim_τ→∞κ

τ

Z τ

0

dt (θ( ¯x1(t)) − zk),

where fk denotes the estimate of the mean force at the center zk. ¯x1(t) is obtained from the

evolution of (2.15_{) at fixed z, namely z ≡ zk}. The calculations of the different f_k are independent from each other, and can therefore be made in parallel. The error made in this estimates is also due to the finiteness ofκ. In order to ensure the convergence of the mean force calculations, one might need to take the value of_κin (2.23) even larger than the one used in the TAMD simulation [24].

Lastly, the free energy profile is globally interpolated based on the local estimates of the mean force. To this end, the free energy is expressed as a linear combination of Gaussians centered on the zk, that is

9_{It is assumed that such a trajectory has visited at least once all the relevant regions of the free energy landscape.} 10_{For this reason, F(z) is also called the potential of mean force.}

2.4. COMMENTS ON THE ACCURACY AND EFFICIENCY OF THE QSS METHOD (2.24) _{F(z) =}˜ K X k=1 ake− (z−zk)2 2σ2 _,

where {ak}andσ> 0 are adjustable parameters, which are optimized by minimizing the so-called

objective function (2.25) E(a,σ_{) =} K X k=1 εk ¯ ¯ ¯ ¯ f_k+ µ∂_F(z)˜ ∂z ¶ z=zk ¯ ¯ ¯ ¯ 2 ,

which measures the discrepancy between minus the gradient of ˜F(z) at zkand the corresponding

estimate of the mean force. Following Poma et al. [16],_{εkmust be set equal to 1/|fk|}2. The detailed procedure to minimize E(a,σ) can be consulted in [24,31].

To conclude this section, let us mention that this reconstruction scheme resembles one of the first techniques developed to compute free energy differences, namely thermodynamic integration [32], in the sense that both methods use local information of the mean force to reconstruct the free energy globally. For one-dimensional systems, it is in fact advisable to use the latter method as it is much simpler than the single sweep. On the other hand, it is more advantageous to employ single sweep in more than one dimension with centers lying on an irregular grid [31].

2.4

Comments on the Accuracy and Efficiency of the QSS

Method

To finish the discussion of the QSS method, let us analyse its accuracy and efficiency. To this end, the method is tested in two one-dimensional model systems. The first model is the harmonic oscillator (2.26) UHO(x) = 1 2mω 2 0x2

withω0 being the frequency of the oscillator, and the second model is the double-well potential

(2.27) UDW(x) = U0 ·µ_x a ¶2 − 1 ¸2 ,

where U0is the height of the barrier of the potential and a gives the location of its minima. These

model systems are very important because they are perhaps the simplest systems in which NQEs such as ZPE and tunneling are present, as mentioned in Section1.3. Let us stress the fact that the harmonic oscillator does not have energy barriers and as such, when using the QSS to study this system we will not be testing the capability of the method to explore rugged free energy

profiles. Nevertheless, it is a good test for the numerical implementation of the method, since we can compare the outcome of simulations with analytical results (see AppendixA). Besides, as we shall see in Chapter4, this systems shows quite well some differences between classical and quantum computations (in particular, see Figure4.1).

In order to study how well the QSS method performs on these model systems, we take the simplest possible choice of the collective variable, that is, the identityθ(x1) = x1. The reconstructed

quantum free energy of the harmonic oscillator, obtained using the methodology exposed in this chapter, is shown in Figure2.2_{for two temperatures T = 30 K and T = 300 K. The frequency}

of the oscillator isω0= 100 THz. For the simulation of the system (2.15), a classical Langevin

thermostat (see Eqn. (1.7)) was used both for the polymer bead coordinates and the collective variable. The friction coefficient_γassociated to the physical degrees of freedom is set equal to 0.1 THz. The coupling constantκfor the first part of the calculation is set equal to 0.5 in atomic units (a.u.). The mass of the oscillator was taken equal to the mass of a proton mp= 1.67262 × 10−27

kg, whereas adiabatic separation between the polymer bead coordinates and z was ensured with an artificial ‘mass’ M = 105× mpand an artificial temperature ¯T = 5 × 105 K for the auxiliary

variable.

The reconstructed free energies (red dashed lines in Figure 2.2) are compared with the exact quantum free energy (blue solid lines). The expression for the quantum free energy of the harmonic oscillator can be derived in a closed form as shown in AppendixA(see Eqn. (A.8)). We can observe that in both cases (T = 30 K and T = 300 K) the QSS method allows to build a free energy profile that is in perfect agreement with the theoretical result. In more precise terms, the relative error of the results shown in Figure2.2is below 0.3% in both cases. The degree of precision is, of course, controlled by the bead number P: The more beads are included in the polymer integrals, the more accurate the results are from the point of view of quantum mechanics. In the present example the accuracy mentioned above was reached using 1350 beads at T = 30 K, and only 125 beads at T = 300 K. The bead number, thus, can also serve as an indicator of how quantum the system is. Indeed, the former system has a thermal wavelength larger than that of the latter.

The examples considered here are really quantum as we are treating them as isolated systems. For example, in the work of Poma et al. [16_{], NQEs of protons embedded in a crystal at T = 380 K}

are computed using the present method, and only 16 beads were enough to give a reasonable accuracy. This is because of many-body effects (the interaction of the proton with the crystal) and the higher temperature used in this work, which somewhat wash away the ‘quantum-ness’ of the system.

The mean forces used in the reconstruction of the free energy profiles shown in Figure2.2are displayed in Figure2.3as a function of the center location. The red points in the figure represent the estimates obtained using Eqn. (2.23), while the blue solid line is a polynomial interpolation of these points. In fact, for the present cases a direct integration of (minus) this interpolation

2.4. COMMENTS ON THE ACCURACY AND EFFICIENCY OF THE QSS METHOD

Figure 2.2: Quantum free energy of a harmonic oscillator in one dimension at two temperatures (a) T = 30 K and (b) T = 300 K. The blue solid curves are the exact (analytic) result, whereas the red dashed lines are the profiles obtained with the QSS method. The bead numbers used in the simulations are (a) P = 1350 and (b) P = 125.

Figure 2.3: Mean forces used in the reconstruction of the free energy profiles shown in Figure2.2. (a) T = 30 K. (b) T = 300 K.

allows us to recover the free energy profile. This is actually the procedure followed to generate the curves in Figure2.2, as well as in Figure2.4for the double-well potential. In all cases, the offset of the reconstructed free energy profile was set such that the minima of this curve coincide with the minima of the exact free energy plot.

As shown in Figure2.3, 11 centers were drawn to perform the free energy calculation for the harmonic oscillator. For each point the evolution of (2.15) at fixed z was follow for 500000 fs, with a relaxation time of 100000 fs. The time-step used here was setδ_{t = 0.1 fs. Note however that} neither the simulation length nor the time-step have to be the same for all centers. As a rule of thumb, the optimal values of these parameters are determined by monitoring the convergence of the average (2.23) to a certain limit value. In this example, nevertheless, the aforementioned values work well for all centers. Moreover, in the case of the harmonic oscillator, it is also possible to derive an analytic expression for the estimates of the mean force, which depends parametrically

on_κ. This is expression is also given in AppendixA(see Eqn. (A.14)). Such an expression allows us to assess the precision of the estimates obtained using Eqn. (2.23) for different values ofκ. Here we found that withκ= 400 a.u. it is possible to keep the relative error below 0.3 % as well. Similar simulations were carried out for a particle of mass mpin the presence of a double-well

potential with parameters U_{0= 0.4 eV and a = 0.5 Å. The results are displayed in Figure}2.4(free energy profiles) and Figure2.5(mean forces). The expression of the quantum free energy of a double-well potential is not known analytically. Therefore, the comparison here is made with the free energy profile obtained by numerically solving the time-independent Schrödinger equation an populating the resulting eigenstates with the proper Boltzmann factor. To show how this is done, let us recast Eqn. (2.3_{) into the basis of the eigenvectors |n〉 of the Hamiltonian of the}

system. We also denote the eigenvalues and eigenfunctions of ˆH as_²nand_ψn(x),11respectively. Introducing a realization of the identity ˆ_{I =}P

n|n〉〈n|, we have F(z) = −β−1_{ln Q}−1X n Z dx 〈x|e−β ˆH|n〉〈n|x〉δ(_{θ(x) − z)} = −β−1ln Q−1X n Z dx 〈x|e−β²n |n〉〈n|x〉δ(θ_{(x) − z)} = −β−1ln Q−1X n e−β²n Z dx 〈x|n〉〈n|x〉δ(θ_{(x) − z)} = −β−1ln Q−1X n e−β²n Z dx |ψn(x)|2δ(θ(x) − z), (2.28)

where in the second line we used that e−β ˆH

|n〉 = e−β²n_{|n〉, and in the last line we expressed the}

inner products of the vectors |x〉 and |n〉 in terms of the eigenfunctionsψn(x). Ifθ(x) = x, as in

the present case, we can do the integral above directly. Moreover, if we express Q in the basis of eigenvectors of ˆH as well, we obtain

(2.29) _{F(z) = −}β−1ln · P ne−β²n|ψn(z)|2 P ne−β²n ¸ ,

which is the expression used to get the exact quantum free energy in Figure2.4. For the numerical calculations, only the first 50 eigenvalues and eigenfunctions of ˆH were considered.

Figure2.4shows us that QSS also allows us to reconstruct the quantum free energy of the double-well potential with a very good precision, although in this case the relative error in some parts of the free energy landscape is slightly bigger than in the case of the harmonic oscillator. This means that a few more beads have to be included in this calculation.

The QSS method has also been tested in more challenging systems [16,23], and it has proven to be useful in including NQEs. For example, in the work of Poma et al. [16], it was found that

11_{The eigenfunctionψ}

n(x) is obtained by taking the inner product of |n〉 with an eigenvector of the coordinate

2.4. COMMENTS ON THE ACCURACY AND EFFICIENCY OF THE QSS METHOD

Figure 2.4: Quantum free energy of a double-well potential in one dimension at two temperatures (a) T = 30 K and (b) T = 300 K. The blue solid curves are the exact results, obtained by numerically solving the Schrödinger equation. The red dashed lines are the profiles obtained with the QSS method. The bead numbers used in the simulations are (a) P = 1350 and (b) P = 125.

Figure 2.5: Mean forces used in the reconstruction of the free energy profiles shown in Figure2.4. (a) T = 30 K. (b) T = 300 K.

the inclusion of quantum effects yields a decrease in the free energy barrier of the process of hydrogen vacancy diffusion in Na3AlH6, improving the agreement with neutron scattering

experiments. The quantum free energy obtained in this work also shows some manifestations of tunnelling effects as explained in Figure2.6. Despite the fact that this calculation was performed in combination with a Density Functional Theory (DFT) based evolution of the potential, and the number of nuclear degrees of freedom per simulation box is 3792 (79 atoms in 3 dimension, with 16 beads per atom), these simulations can be done using computational resources that are available nowadays. This illustrates the reasonable numerical cost of the method when employed to study a realistic condensed-phase system.

A very positive feature of QSS (as well as of the ‘classical’ single sweep method) is the fact that, the evaluation of the mean forces (2.23), which is the bulk of the calculations as in any molecular dynamics (MD) scheme, can readily be parallelized as it is composed of independent simulations.

Figure 2.6: Classical (black curve) and quantum (grey curve) free energies for the process of H vacancy diffusion in Na₃AlH₆. The quantum free energy was calculated via the quantum single sweep method. The shift of the minima in the quantum case towards the origin with respect to the classical case, is often regarded as a manifestation of tunnelling effects. (Reproduced from [16].)

It has also been argued that the efficiency of the method should be kept when a larger number of collective variables is considered [23]. In spite of these factors, methods employing the path integral formalism are in general computationally expensive, tending to be one or two orders of magnitude more costly than treating the nuclei classically. This is due to the fact the number of degrees of freedom in the original description of the system is multiplied by the bead number P, that has to be large enough to give accurate results.

2.5

Other Methods to Compute Quantum Free Energies

For the sake of completeness, alternative methods to calculate quantum free energies are briefly reviewed in this section. The list shown here is by no means exhaustive; it nonetheless gives an idea of the alternative methods that exist to carry out this sort of calculations. It is noteworthy that most of these methods are formulated in the same spirit as the QSS scheme, that is, by recasting a method for computing classical free energies into the language of path integrals presented in this chapter.

• Wigner-Kirkwood expansion: A first route to include quantum corrections to the free energy is the Wigner-Kirkwood expansion, which essentially consists of an expansion of the quantum Boltzmann density in powers of Planck’s constant_~. Using this expansion it is possible to show [33] that the Helmholtz free energy, to second order in_~, is given by

2.5. OTHER METHODS TO COMPUTE QUANTUM FREE ENERGIES (2.30) _{F ≈ Fcl+}β 2 ~2 24 N X i=1 1 mi〈(∇i U)2_〉cl

where the subindex ‘cl’ stands for ‘classical’, ∇iindicates differentiation with respect to the

Cartesian coordinates of the i-th particle, and U is the total potential energy of the system. • Non-Markovian path integral Monte Carlo: Proposed by Crespo et al. [34], this scheme is

implemented as a history-dependent MC using a path integral formulation. The method was inspired by the Wang-Landau sampling (e.g. see Section 7.6 of [1]) and the metadynamics method (see Section1.2). In the original work by Crespo and coworkers the free energy of the two-dimensional quantum Ising model was computed, showing an accuracy comparable with the best methods to study such kind of systems, such as schemes based on the Wang-Landau method. The advantage of this approach over the latter ones is that it can easily be applied to study off-lattice quantum problems.

• Nonequilibrium path integrals: This method was introduced in 2008 by van Zon et al. [35]. This is an extension to the quantum realm of classical nonequilibrium methods based on the Jarzynski equation [36]

(2.31) _〈e−βW_{〉 = e}−β∆FAB_,

where the brackets denote an ensemble average over measurements of W, the irreversible work made on the system for bringing it from state A to state B. Therefore, the left-hand side of (2.31) is an out-of-equilibrium quantity. On the contrary, the right-hand side involves the equilibrium free energy difference∆F_AB. The quantum method is built by expressing nonequilibrium relations such as (2.31) in the language of continuous path integrals. • Ab-initio path integral metadynamics: This is a hybrid methodology that builds upon a

technique for sampling classical free energies, namely, metadynamics, and the path integral approach explained in this chapter, with the particularity that the interatomic interactions, i.e., the potential U(x), are computed on-the-fly from electronic structure calculations. This is why the name of the method starts with the term ab-initio. The method was developed very recently by Ivanov et al. [37] and it has brought some light on elucidating the reaction mechanism behind double proton transfer in the formic acid dimer. As pointed out in the previous section, the idea of coupling path integrals with ab-initio calculations has also been employed within the context of QSS, see [16].

C

H A P T E R

3

T

HE

Q

UANTUM

T

HERMAL

B

ATH

Q

uantum thermal baths (QTBs) are stochastic models that approximately include NQEs through an effective noise and a dissipative force, which are adjusted so as to have the power spectral density given by the quantum fluctuation-dissipation theorem [22,38]. In this approach, the description of the system of interest is made in terms of a reduced set of dynamical variables with the remaining degrees of freedom being taken into account by the effective noise mentioned above. Therefore, the computational cost of QTB methods is relatively low, in particular when compared with path integral techniques [39]. In this chapter, we study a particular kind of QTB that was introduced by Dammak et al. [38]. The theoretical aspects of this approach are detailed in Section3.1. Next, the method is tested using two model systems in one dimension: The harmonic oscillator and the double-well potential. While QTB is in principle exact for the former system, it just leads to approximate, yet good qualitative results for the latter. These applications are presented in Section3.2. Some words of caution are given in Section3.3, where the drawbacks and limitations of the QTB approximation are discussed. Finally, in Section

3.4an approach proposed recently by Brieuc et al. [15] combining path integrals with QTB is analysed. This latter part makes the connection between the present chapter and the former one.

3.1

Theoretical Aspects of QTB

The QTB technique for MD simulation was introduced in 2009 by Dammak et al. [38] and Ceriotti et al. [40]. The basic idea in both papers is to use a thermostat such as the one given by Eqn. (1.7), i.e., a Langevin thermostat but with a coloured noise in place of the white noise_η(t).1The main difference between these two approaches is that in [40] a non-Markovian friction term is used