Marc Hayoun

Je Voudrais Ensuite Remercier

Fabio Finocchi

Philippe

Merci À Vous

( Hichem

Fabio Marc

Je Tiens À Remercier

Marie-Pierre Gaigeot

Florent Calvo

Rodolphe Vuilleumier

Grégory Geneste

Hichem Dammak

Rodolphe Vuilleumier D'

Gentien, Charlotte, Camille, Sergei Anastasia Charles

Pavan, Bobo, Cintia, Yousra, Damien, Zied, Mickaël, Allistar, Romain, Tahar, Désirée, Hedia, Halyna Mohamed Yang

Xiaofei

Chen Bertrand

Yi

Yin Ze

Bon

S ) Guilhem

G Brahim

Ingrid Agnès

Gilles

Nicolas

Claire, Thierry etc.. Hubert

Antoine M Louise

Hélène Marie Yan

Je tiens à remercier particulièrement

Keywords:

aurais pas pu réussir cette thèse.

Chapter 1

General introduction 1.1 Nuclear quantum effects

Numerical simulations are now routinely used to investigate the properties of physical systems. In particular, since the pioneer work of Fermi, Pasta, Ulam and Tsingou [1,2], Molecular Dynamics (MD) has become a powerful tool to study condensed matter systems at the atomic scale. MD allows the simulation of the trajectories of the atoms that composed the system by numerically solving their equations of motion. Interatomic forces can be described using different methods such as using force fields or ab initio methods from the electronic structure of the atoms (using first-principle methods such as density-functional theory), while the dynamics of the nuclei is described using classical mechanics (i.e. the atoms follow Newton's equations).

Thus, the quantum nature of the nuclei is not taken into account in standard MD simulations and quantum effects such as quantum fluctuations and delocalization, zero-point energy and quantification or tunnelling effect for example are neglected. However, these nuclear quantum effects can have a major impact on the properties of the system, in particular at low temperatures and/or for systems containing light atoms such as hydrogen. For example, quantum effects are necessary to correctly describe the phase diagram of different systems such as water [3] or ferroelectric crystals like BaTiO 3 or SrTiO 3 [4,5]. The study of isotope effects, observed in various systems, such as lithium hydride [6], ice [7,8] or even biological systems [9,10], also requires the inclusion of nuclear quantum effects. Finally, as a last example, the low temperature behavior of the heat capacity of solids, which vanishes at T=0 K, cannot be reproduced if the quantum statistics of the atoms is not taken into account.

In order to account for nuclear quantum effects in simulations, path integral methods, such as path integral molecular dynamics (PIMD), are often used [11]. PIMD allows the computation of exact quantum results for static equilibrium properties, however dynamical properties are not directly accessible. Moreover, PIMD simulations are computationally demanding, and thus, are generally limited to the study of rather small systems (in particular if the interatomic forces are described using first principle methods). Approximate methods based on a modified Langevin dynamics have been proposed as an alternative to PIMD [12,13]. Among them, is the quantum thermal bath (QTB) method [12]. The method is particularly appealing because of its computational cost that is equivalent to standard MD thus allowing the study of large and complex systems. Although approximate in highly anharmonic cases, the QTB method has been able to give satisfactory results for various anharmonic systems [6,[14][15][16][17][18][19]. Moreover, dynamical properties are directly accessible in QTB-MD simulations. The first part of this thesis is devoted to the study of the QTB method.

Proton conduction in perovskite materials

In a second part of this thesis, we study proton (H + ) conduction in perovskite materials. Proton conduction in solids is a key process for various applications such as fuel cells, electrolysers, sensors, gas separation... [20]. In particular, good proton conductors are needed for application in hydrogen fuel cells or, more precisely, in proton conducting fuel cells (PCFC). The operating principle of such fuel cells is schematically represented in figure 1.1 -left panel. Perovskite materials can exhibit relatively high proton conductivities, and thus have been extensively studied for application in fuel cells, in particular as electrolyte materials in PCFC. These perovskite materials are oxides of formula ABO 3 whose structure is represented in figure 1.1 -right panel. Among these oxides, doped barium zirconate (BaZrO 3 ) combines one of the highest proton conductivity with a good chemical stability which makes it a promising material for future applications as electrolyte material in PCFC [21].

In these oxides, protons are generally introduced by doping the structure on the B site. A fraction of the B site cations is substituted by another cation B' of lower valence. This results in an excess of negative charges which is compensated by the creation of oxygen vacancies (V •• O in Kröger-Vink notation1 ). Then the compound is placed under wet atmosphere in order to incorporate water molecules that will fill some of the vacancies and create protonic defects in the structure (OH • O ) through the following hydration reaction:

H 2 O + V •• O + O × O → 2OH • O (1.1)
After hydration, the incorporated protons (H + ) occupy interstitial positions close to the oxygen sites and can potentially diffuse. The proton generally remains covalently bonded to an oxygen atom creating a protonic defect OH -. Long range migration of the proton in perovskite materials arise from a combination of two motions: transfer and rotation (or reorientation) [21]. During the transfer step, the proton jumps between two neighboring oxygen atoms, while during the rotation step, the proton remains covalently bonded to the nearest oxygen atom and rotates around it. It is generally accepted in the literature that the rotational motion of the proton is fast compared to the transfer step, suggesting that proton transfer is rate limiting in these materials [21]. The proton motion in perovskite materials is shown schematically on figure 1.2.

Quantum effects and proton conduction

In this thesis, we investigate in particular the importance of quantum effects on the diffusion of hydrogen in perovskite materials. Since hydrogen is the lightest element, we expect quantum effects to have a significant impact on its diffusion. The role of quantum effects on hydrogen diffusion has been mainly studied for hydrogen in metals and on metallic surfaces [22][23][24], and it has been demonstrated, that the diffusion exhibits two different regimes [25][26][27]. A classical regime, at high temperature, for which the diffusion coefficient follows the classical Arrhenius law:

D(T ) = D 0 exp(-βE a ) (1.2)
with E a the activation energy associated with the diffusion and a quantum regime, at low temperature, where quantum effects becomes significant and leads to a deviation from the previous Arrhenius law. This behavior can be seen in figure 1.3 for the experimental diffusion coefficient of hydrogen in a niobium crystal. We clearly see that at low temperatures (i.e. for T 250 K) the diffusion coefficient deviates from the high temperature Arrhenius behavior. This is associated to quantum effects, in particular zero-point energy and tunnelling effect, that become non negligible at low temperatures. In this thesis, we study the diffusion of protons (H + ) in barium zirconate (BaZrO 3 ), using PIMD and QTB-MD simulations, in order to investigate if a similar behavior is to be expected for proton diffusion in perovskite materials. • The Empirical Valence Bond model (EVB) to include chemical reactivity for the study of proton diffusion based on the potential proposed by Raiteri et al. [28] • The possibility to run PIMD and QTB-PIMD simulations

• The possibility to run constrained simulations using the SHAKE algorithm [29,30] in order to perform, in particular, free energy computation using thermodynamic integration in the blue moon ensemble [31,32] 1

.5 Construction of the manuscript

The manuscript is made up of six chapters organized in three parts. The first part is introductory and contains a general introduction (chapter 1) and the presentation of the different numerical methods (Standard MD, PIMD and QTB-MD) used during the PhD (chapter 2). In the second part, we first study the QTB method focusing in particular on its major limitation: the zero-point energy leakage problem (chapter 3). Then, we present another way to use the QTB method as a thermostat for PIMD simulations (chapter 4). The last part is devoted to the study of proton conduction in perovskite materials. First, we study the impact of quantum effects on the diffusion of hydrogen in BaZrO 3 (chapter 5). Then, we study the proton diffusion mechanisms in GdBaCo 2 O 5.5 (chapter 6) in which a strong correlation between O and H diffusion is predicted.

Chapter 2

Nuclear quantum effects in molecular dynamics simulations

Introduction

As explained in the general introduction (chapter 1), nuclear quantum effects are generally neglected in simulations. However, the quantum nature of the nuclei can play a major role on the properties of the system in particular at low temperature and/or in systems containing light atoms such as hydrogen. In these cases, quantum effects such as zero-point effects or quantum fluctuation/delocalization need to be somehow included in simulations. Let us already note that quantum effects associated with exchange of identical particles are not taken into account in the following, and thus the fermionic or bosonic nature of the particles is neglected. 1Since direct resolution of Schrödinger equation is only possible for a few degrees of freedom, other methods are needed. A standard way to account for nuclear quantum effects is path integral based methods [1,2] such as Path Integral Molecular Dynamics (PIMD) or Path Integral Monte Carlo (PIMC). These methods are based on Feynman path integral formulation of quantum mechanics [3], and allows to exactly compute quantum static properties even for highly anharmonic systems. However, one important drawback is that dynamical quantities are not directly accessible in the path integral framework [4].

Computing quantum dynamical properties is still a great challenge in computational chemistry and physics nowadays. Over the years, several approximate techniques has been proposed to deal with this issue. In particular, methods like Centroid Molecular Dynamics (CMD) [5,6] and Ring Polymer Molecular Dynamics (RPMD) [7] has been developed to go beyond PIMD in order to compute time correlation functions. These two methods have been able to give satisfactory results on several systems [8][9][10][11], but the computation time required reduces their range of applicability 2 .

Another important class of methods uses semiclassical ideas to adapt classical methods such as molecular dynamics to the quantum case. These semiclassical methods are widely used because they provide a good ratio between computational cost and accuracy. Numerous methods lies into this class [12] including for example the semiclassical initial value representation [13], in which semiclassical formulas are used to express the time evolution propagator as a phase space average over the initial conditions of classical trajectories, the forwardbackward semiclassical dynamics [14], which uses similar semiclassical ideas with a combina-tion of the time evolution propagator and its adjoint to compute time correlation functions, and the Liouville dynamics method [15,16] based in particular on a quantum generalization of the Liouville theorem.

Recently, alternative methods based on a modified Langevin equation have been proposed [17,18]. Among them is the quantum thermal bath (QTB) method which includes nuclear quantum effects in MD simulations through a modified (quantum) Langevin thermostat. Although only exact in the harmonic case [19,20], the QTB method has been able to give satisfactory results for various anharmonic systems [21][22][23][24][25][26][27]. A first advantage of the method is its universal and clear formulation which makes its implementation in an existing MD code quite straightforward. Moreover, the computational cost of QTB-MD simulations is similar to standard MD simulations, and dynamical properties can, in principle, be computed directly.

In this chapter we first present the basics of molecular dynamics and introduce the Langevin thermostat. Then, we present the basics of the QTB and PIMD methods. Finally, we study the behavior of the QTB method on simple one-dimensional systems with increasing degree of anharmonicity and compare the results to PIMD. Since the QTB method is only exact for harmonic systems, we investigate in particular the limitations of the method for strongly anharmonic cases.

Standard molecular dynamics

Molecular dynamics (MD) is a powerful numerical method to investigate the properties of complex systems. The idea is to numerically solve the equations of motions of the atoms or molecules that compose the system. As already explained in the general introduction, the dynamics is classical i.e. the atoms follow Newton equations. These equations of motion can be numerically solved using finite differences methods. The time is discretized with a time step δt and the equations of motions can be integrated using the Verlet algorithm [28]

r i (t + δt) = 2 r i (t) -r i (t -δt) + δt 2 m i f i (t) + O(δt 4 ) (2.1)
with f i (t) the force that applies on the atom i at time t.

In simulations, we are interested in computing the values of some macroscopic quantities, often called observables, which can be compared to experimental values. These observables are ensemble average values A . If the ergodicity hypothesis is valid, these statistical average values can be computed using molecular dynamics

A = lim Ns→∞ δt→0 Ns n=1 A n N s (2.2)
where A n represent the value of A evaluated at time t = nδt and the time duration of the simulated trajectory is N s δt. The typical output of an MD run is the trajectory of the atoms i.e. the sequences { r i (t n )}, { v i (t n )}; n = 1, ..., N s ; i = 1, ..., N from which average values of any observable A({ r i }, { v i }) can be computed using expression (2.2).

For now, we have considered an isolated system so its total energy is conserved. Thus an accurate integration of the equations of motion is supposed to ensure a good energy Ph.D. thesis -F. Brieuc -2016

Standard molecular dynamics

conservation during the simulation. In order to achieve this, the time step δt needs to be small compare to the smallest characteristic time of vibration in the system.

δt << 1 ν max (2.3) with ν max the highest vibration frequency.

One often wants to study a system in the canonical ensemble (N, V, T ) i.e. in contact with a thermal bath. In this case, a thermostat is needed to perform MD simulations at constant temperature. Several types of thermostats exist [29] such as velocitiy rescaling, the Nosé-Hoover thermostat or the Langevin thermostat. In all our simulations we use the Langevin thermostat which samples correctly the canonical ensemble and ensures a very good ergodicity. We describe this thermostat in details in the next section.

The Langevin thermostat

The Langevin thermostat is a well known and widely used thermostat based on the Langevin equation ṗ = fγp + R (2.4) This equation has been proposed in 1908 by P. Langevin [30] to described the Brownian motion of a particle in a fluid (i.e. in the f = 0 case). More generally, the Langevin equation describes the dynamics of a system (here for one particle in 1D) in contact with a thermal bath. The last two terms are the forces associated to the interaction between the system and the bath. The first term -γp is a friction force and the last term R(t) is a stochastic force.

In standard Langevin dynamics the random force R(t) is supposed to have the following properties :

• It is a stationary process.

• Its distribution is Gaussian with zero mean : R(t) = 0 (2.5)

• Its autocorrelation function is

R(t)R(t + τ ) = 2mγk B T δ(τ ) (2.6)
which means that the values of R(t) are uncorrelated.

Since R(t) is a stationary process, it obeys the Wiener-Khinchin theorem which states that the power spectral density I R and the autocorrelation function of R(t) are related by Fourier transform

I R = +∞ -∞ R(t)R(t + τ ) e -iωτ dτ (2.7)
Since the values of the stochastic force are supposed to be uncorrelated (eq. (2.6)), d the power spectral density (PSD) of R(t) is a constant

I R = 2mγk B T (2.8)
thus the stochastic force is a white noise. This expression for the spectral density comes from the classical expression of the fluctuation-dissipation theorem [31] that relates the dissipation in the system (γ) with the thermal fluctuations that come from the stochastic force I R .

With these properties for R(t) it can be shown that the system reaches thermal equilibrium with the proper canonical Boltzmann distribution [32]. In particular, the equipartition theorem holds

p 2 2m = k B T 2 (2.9)
Thus Langevin dynamics can be used as a tool to sample the canonical ensemble in MD simulations. By adding friction and stochastic forces, it is possible to ensure the canonical distribution at a fixed temperature T . In this case, the friction and stochastic forces are only used as a tool to thermalise the system, and the friction coefficient is then a free parameter which has to be correctly chosen.

Numerical aspect of the Langevin thermostat

First, the friction coefficient γ should be chosen small enough so that the forces of the thermostat do not perturb the natural dynamics of the system. Thus the characteristic times in the dynamics should be small compared to the characteristic time of exchange with the thermostat (1/γ). In other words γ << ν c (2.10) where ν c is a characteristic vibrational frequency of the system. Too high values of γ will results in an overdamped dynamics that will primarily affect time correlation functions such as the diffusion coefficient or the vibrational spectrum.

In contrast, the time needed for the system to reach the equilibrium is of the order of a few times 1/γ. Thus the duration of a trajectory should be very high compared to this time

N s δt >> 1 γ (2.11)
where N s is the number of time steps used in the trajectory and δt the time step. A very low value of γ thus means that very long trajectories will be needed.

A good value for the friction coefficient has to be found before launching production runs. In practice, one studies the evolution of a physical quantity O with γ. If the friction coefficient is low enough, the values of O should be independent of γ. Then, one chooses the highest value of γ for which O remains almost constant.

Langevin dynamics can be numerically integrated using the standard Verlet algorithm (eq.

(2.1)). Since the friction force of the Langevin thermostat depends on v(t) one needs to compute the velocity at time t which can be done using the following expression

v(t) = 3x(t) -4x(t -δt) + x(t -2δt) 2δt + O(δt 2 ) (2.
12)

The last thing one needs is to generate the random force R(t) with the correct properties. In particular, the random force should have the correct spectral density I R . From the Wiener-Khinchin theorem (eq. 2.7) one obtains that

R(t)R(t + τ ) = +∞ -∞
I R e iωτ dω 2π (2.13) and with the PSD given in equation (2.8)

R(t)R(t + τ ) = 2mγk B T +∞ -∞
e iωτ dω 2π (2.14) In simulations, the time is discretized with t = nδt leading to the angular frequency step δω = 2π/(N s δt) with ω k = kδω. The inverse Fourier transform takes the standard discrete form

R(nδt)R((n + n ′ )δt) = 2mγk B T 1 N s δt
Ns-1 k=0 e i2πkn ′ /Ns (2.15) Ns (2. 16)

R(nδt)R((n + n ′ )δt) = 2mγk B T δt 1 N s 1 -e i2πn ′ 1 -e i2πn ′ /
R(nδt)R((n + n ′ )δt) = 2mγk B T δt δ n ′ ,0 (2.17)
One sees that after discretization of time the Dirac function δ(τ ) becomes δ n ′ ,0 /δt where δ n ′ ,0 is the Kronecker delta 3 . Thus the correct random force can be obtained by first generating uncorrelated random numbers distributed with the normal distribution (N (0, 1)) and then multiply those numbers by 2mγk B T /δt.

The quantum thermal bath method

The quantum thermal bath (QTB) method has been proposed by Dammak and coworkers in 2009 [17]. The main idea is to modify the standard Langevin thermostat in order to include quantum effects associated with the dynamics of the nuclei in MD simulations. In this section we present the basics of this method.

The dynamics obtained using the Langevin thermostat is classical by construction. In particular, the equipartition theorem holds and thus every harmonic vibrational modes have the same average energy (k B T ). In the quantum case, one expects the average energy of an harmonic vibrational mode to be given by

θ(ω; T ) = ω 1 2 + 1 exp(β ω) -1 (2.18)
which depend on its angular frequency ω. The main idea of the QTB method is to modify the Langevin thermostat to ensure the quantum energy distribution given by equation (2.18).

The properties of the dynamics obtained from Langevin type equations are directly related to the properties of the random force, and in particular to its power spectral density I R . This spectral density is obtained from the fluctuation-dissipation theorem which writes [31,32]:

χ′′ (ω) = ω 2κ(ω; T ) I x (2.19)

Chapter 2. NQE in MD simulations

with κ(ω; T ) the average energy of a vibrational mode of angular frequency ω at temperature T , χ′′ (ω) the imaginary part of χ(ω) and χ(ω) the susceptibility that connects the Fourier transform of the position x(ω) to the Fourier transform of the stochastic force R(ω):

x(ω) = χ(ω) R(ω) (2.20)
Using the definition of the power spectral density 4 we obtain

I x = | χ(ω)| 2 I R (2.21)
and the fluctuation-dissipation theorem can be written in the form

I R = 2κ(ω; T ) ω χ′′ | χ(ω)| 2 (2.22)
The susceptibility χ(ω) is obtained by Fourier transform of the equation of motion (i.e. the Langevin equation -eq. (2.4)). In the case of an harmonic oscillator of angular frequency ω 0 one finds that the susceptibility is given by

χ(ω) = 1 m ω 2 0 -ω 2 -iγω (2.23)
and from equation (2.22), we obtain the following expression for the PSD of the stochastic force

I R = 2mγκ(ω; T ) (2.24)
In the classical version of the fluctuation-dissipation theorem

κ(ω; T ) = k B T (2.25)
and the spectral density is thus a constant given by eq. (2.8) leading to the classical Langevin dynamics of the previous section.

A generalization of the fluctuation-dissipation theorem for quantum systems developed by Callen and Welton in 1951 [33] gives

κ(ω; T ) = θ(ω; T ) = ω 1 2 + 1 exp(β ω) -1 (2.26)
which leads to the following power spectral density for the random force

I R (ω) = 2mγθ(ω; T ) = 2mγ ω 1 2 + 1 exp(β ω) -1 (2.27)
Several authors have proposed to use this expression for the spectral density of the random force to design a Langevin equation valid for quantum systems [34,35]. The quantum thermal bath method apply this idea to MD simulations by modifying the Langevin thermostat: the PSD of the random force is modified according to equation (2.27).

It is worth noting that, first the spectral density of R(t) now depends on ω: the random force is not a white noise anymore and is referred to as a colored noise. Second, the function θ(ω; T ) is the average energy of an harmonic oscillator of angular frequency ω. For high temperatures (i.e. typically for temperature higher than the Debye temperature Θ D5 ) we reach the classical limit and θ(ω, T ) ≈ k B T . We recover the classical expression of the PSD and thus the equipartition of the energy. On the contrary, for T Θ D the function θ(ω; T ) diverges from the classical behavior and quantum effects become important. In particular, in the limit T → 0 the function θ(ω; T ) → ω 2 and the average energy is different than zero because of zero-point energy.

Note also that with the PSD of equation (2.27) the average kinetic energy obtained for a vibrational mode of angular frequency ω is θ(ω; T )/2 which is only exact in the harmonic approximation. Thus the QTB method is expected to give the exact energy distribution only in the harmonic case. This is a first limitation of the quantum thermal bath: unlike PIMD, the quantum thermal bath method becomes approximate for anharmonic systems. Nevertheless, the quantum thermal bath method has been able to give satisfactory results in various anharmonic systems [21][22][23][24][25][26][27]. But the method should still be used with caution for strongly anharmonic systems, in particular because energy is transferred from vibrational modes of high frequencies to vibrational modes of low frequencies leading to the zero-point energy leakage problem. This limitation is studied in details in chapter 3.

Finally, in QTB-MD simulations, the equipartition theorem does not hold, the average kinetic energy of any degree of freedom K includes quantum effects and thus is greater than k B T /2 (or 3k B T /2 in 3D). So the effective temperature T * = 2 K /k B (or 2 K /3k B in 3D) of QTB-MD simulations is generally greater than the target temperature T of the thermostat. The effective temperature T * becomes equal to T only in the classical limit (i.e. at high temperatures).

Another method using colored noise in Langevin type equations has been proposed by Ceriotti and coworkers [18] to include quantum effects in MD simulations. This method is based on a more general equation of motion called the Generalised Langevin Equation (GLE). In this GLE method, the quantum effects are introduced through a frequency dependent friction coefficient whereas, in the QTB case, quantum effects are included through the power spectral density of the random force. The two methods are basically equivalent although, the GLE method requires careful and complex optimization of several parameters in order to recover the quantum fluctuations. Thus the QTB method is simpler to use and to implement in an existing MD code.

Numerical aspect of the quantum thermal bath method

In the QTB method, one needs to generate the random force with the correct power spectral density. Since the stochastic force is not a white noise anymore, the standard way to generate the random force in Langevin thermostat cannot be used. Thus we use a more general procedure [36] to generate a random force with a target spectral density I R . Since in MD simulations the time is discretized t n = nδt we want to generate the values R n = R(t n ). The procedure is first to generate the stochastic force in the Fourier space Rk = R(ω k ) with ω k = kδω and, for a trajectory of total duration N s δt, the angular frequency step δω = 2π/(N s δt). Then the forces R n are obtained using a discrete Fourier transform of Rk . This procedure is described in more details in ref. [37], here we only present the main steps necessary to implement it.

• Generate independent Gaussian random numbers a k and b

k for k = 1, ..., N s /2 -1 • Compute Rk for k = 1, ..., N s /2 -1 using the following expression Rk = N s δt 2 I R (ω k ) (a k + ib k ) (2.28)
• Generate Rk for k = N s /2 + 1, .., N s -1 using the following symmetry property

Rk = R * Ns-k (2.29)
which comes from the fact that R(t n ) has to be a real function.

• Set R0 and RNs/2 to zero

• Compute the values for R n using a discrete Fourier transform

R n = 1 N s δt Ns-1 p=0
Rk e i2πnk/Ns (2.30)

In this procedure, the random forces for the entire time duration of the trajectory are generated before launching the simulation. An alternative way to generate the stochastic forces has been proposed by Barrat and Rodney which allows to generate the forces R(t n ) "on the fly" during the simulation [38]. In this paper, the authors also emphasize the importance of using an angular frequency cut-off ω cut when generating the stochastic forces. So the QTB method finally contains two parameters : the friction coefficient γ and the angular frequency cut-off ω cut . The values for these two parameters have to be carefully chosen. Since the QTB method is based on the Langevin thermostat the friction coefficient γ should be chosen according to the same prescriptions. In particular, γ should be small enough to ensure that the forces associated with the thermostat does not modify the natural dynamics of the system.

Moreover, the angular frequency cut-off ω cut has to be higher than the highest angular frequency observed in the system ω max . However, ω cut can not be chosen arbitrarily high.

It should be chosen of the order of a few times ω max to prevent the inclusion of too high frequencies which could lead to the divergence of the energy. More precisely this divergence behaves as γ ln(ω cut ) [38] and so for a fixed value of ω cut the friction coefficient γ also has to be small enough to avoid any divergence. We found in our simulations that ω cut ≈ 2ω max is an acceptable value which seems to avoid energy divergence for any reasonable values of γ.

Path integral molecular dynamics

Path integral formulation of quantum mechanics has been developed by R.P. Feynman, and is described in details in the book he wrote with A.R. Hibbs in 1965 [39]. In this book, they show that the canonical density matrix of a quantum system can be written in a path integral form which is the basis of path integral molecular dynamics. Starting from this expression for the density matrix, we present the basics of the PIMD method and how quantum static properties can be numerically computed using this method.

The classical isomorphism

For simplicity, we consider the case of one particle of mass m in one dimension (position x) described by the Hamiltonian

Ĥ = K + V = p2 2m + V (x) (2.31)
where K and V are the kinetic and potential energy operators, respectively. The generalization to N particles in three dimensions is straightforward. In the canonical (N, V, T ) ensemble the density operator or density matrix of the system is given by

ρ(β) = e -β Ĥ (2.32)
whose matrix elements in position representation are

ρ(x, x ′ ; β) = x| ρ(β) |x ′ = x| e -β( K+ V ) |x ′ (2.33)
The density matrix can be approximated using the Trotter product formula leading to the following expression

ρ(x, x ′ ; β) ≈ mP 2πβ 2 P/2 dx 1 ... dx P -1 exp - P s=1 mP 2β 2 (x s-1 -x s ) 2 + β P V (x s ) (2.34)
called the discrete imaginary time path integral expression of the density matrix (for a derivation of this formula see the complements of this chapter -section 2.7)

The canonical partition function Z is related to the canonical density matrix through

Z = T r[ρ(β)] = dx x| ρ(β) |x = dx ρ(x, x; β) (2.35)
And the average value of any observable A is given by

A = 1 Z T r[Aρ] (2.36)
Using expression (2.34) for the density matrix with a small change in notation writing x ≡ x P , we can express the canonical partition function as

Z ≈ Z P = mP 2πβ 2 P/2 dx 1 ... dx P exp - P s=1 mP 2β 2 (x s-1 -x s ) 2 + β P V (x s ) (2.37)
with cyclic conditions on the path {x s } i.e. x P = x 1 which comes from the trace operation. Let us recall that expression (2.34) for the density matrix and thus this expression for the canonical partition function are both exact in the limit P → ∞ so

Z = lim P →∞ Z P (2.38)

Chapter 2. NQE in MD simulations

The constant prefactor in the expression of Z P is actually the result of a P -dimensional Gaussian integral and can be rewritten mP 2πβ with m ′ = mP . The value of this prefactor does not affect the results of any average value, and thus the value for the parameter m ′ can be chosen arbitrarily [4]. A natural choice for m ′ is the physical mass m of the particle but other standard (and more complex) choice for m ′ exist [40]. In the following we use m ′ = m in all the PIMD calculations.

Finally, the quantum canonical partition function can be written in the following form

Z ≈ Z P = 1 (2π ) P dp 1 ... dp P dx 1 ... dx P exp (-βH P ) (2.41)
with

H P = P s=1 p 2 s 2m + 1 2 mω 2 P (x s-1 -x s ) 2 + 1 P V (x s ) (2.42)
which is the classical partition function of a system defined by the effective Hamiltonian H P . So the quantum partition function is formally equivalent to the classical partition function of an extended system : this is called the classical isomorphism. This extended system is composed of P particles (called replicas or beads) that are submitted to the external potential V /P . Each replica s interacts with its "nearest neighbour" (s ± 1) via an harmonic coupling of angular frequency ω P = √ P /β . In the case of a many-body system (i = 1, .., N ) the extended system is composed of N × P replicas (or beads). Each particle i is described by an ensemble of P beads (s = 1, .., P ) interacting with each other via the harmonic coupling 1 2 m i ω 2 P s (x i,sx i,s-1 ) 2 . Then each replica s of particle i also interacts with replica s of the other particles j = i through the effective potential V (|x j,sx i,s |)/P . The classical isomorphism is schematically represented in figure 2.1. Let us note that we consider distinguishable particles here. If the particles are indistinguishable, then the effects associated with the exchange of particles needs to be taken into account 6 . In the following, we will only consider the case of distinguishable particles.

For a finite value of P , one can estimate quantum average values for any time-independent observable A by exploiting the classical isomorphism

A ≈ A P = 1 Z P dp 1 ... dp P dx 1 ... dx P A(x, p) exp (-βH P ) (2.43)
These average values can be evaluated using standard molecular dynamics algorithms leading to the Path Integral Molecular Dynamics (PIMD) method 7 . In practice, we run standard MD simulations in the extended system composed of N × P particles to compute average quantities.

V/P V ω P ω P Figure 2.1: Schematic representation of the classical isomorphism between a quantum system and the extended classical system. Here for a system composed of N = 2 particles and for a Trotter number P = 6. Thus, in the extended system each particle is represented by a "ring polymer" composed of P = 6 beads (or replicas). The polymer is closed because of the cyclic conditions (xP +1 = x1) coming from the trace operation in equation (2.36).

Note again that expression (2.41) is formally exact only in the P → ∞ limit which of course is not accessible in simulations. However, the method converges to the exact result when the number of replicas increases. Thus, one can use a finite number of beads P although one has to be sure that P is high enough to ensure a good convergence. In practice, we compute the quantities that we want to study for different values of P in order to find how many replicas are needed to reach the convergence. An example of this convergence is presented in figure 2.2. Once the correct value for P has been found for a calculation at temperature T , the product P × T should be kept constant for every simulations to keep the same level of convergence. This is one important limitation of PIMD : when the temperature investigated becomes very low, the number of beads necessary to converge becomes too high, and the computation time becomes prohibitive.

Note also that here molecular dynamics is only used to explore the phase space and sample the canonical distribution to compute time-independent average values. The dynamics in the extended system cannot be used to compute time-dependent quantities [4]. Indeed, expression (2.34) for the density matrix, and thus expression (2.41) for the partition function comes from a path integral formulation in imaginary time. So the time associated with the dynamics in the extended classical system is not the physical time. For this reason, dynamical quantities such as diffusion coefficent or vibrational spectrum are not directly accessible in the path integral framework. This is probably the major limitation of the PIMD method. Several approximate methods such as ring polymer molecular dynamics [7] or centroid molecular dynamics [5,6] has been designed to solve this problem. These methods have been able to give satisfactory results on several systems however their computational cost is even higher than that of PIMD which reduces their range of applicability.

Finally, it has been shown that ergodicity can sometimes be difficult to reach in PIMD when using an important number of beads [41]. This problem can generally be avoided by using a stochastic thermostat. For this reason, we use the Langevin thermostat to ensure ergodicity and the canonical distribution at the target temperature 8 . Thus the equation of motion of

Chapter 2. NQE in MD simulations

a replica s in our PIMD simulations finally writes

ṗs = f s -mω 2 P (2x s -x s+1 -x s-1 ) -γp s + R s (2.44)
where f s = -(1/P )∂V /∂x s is the external force exerted on replica s and the two last terms are the frictional (friction coefficient γ) and stochastic forces of the Langevin thermostat.

Computation of macroscopic properties

In simulations we are interested to compute macroscopic quantities which are given by average values of the form given in equation (2.36). In the previous section we saw that, in the PIMD method, average values are obtained using equation (2.43). In the following, we will review some standard expressions to compute average values in PIMD. For simplicity in notations we will drop the subscript for PIMD average so that .. ≡ .. P . Thus average values now are computed in the extended phase space of the system with P replicas. The details for the derivation of these expressions are given in the complements of this chaptersection 2.7.

Starting from the usual thermodynamic relation for the average energy

E = - 1
Z P ∂Z P ∂β (2.45)
one obtains the following expression for the average total energy

E = P s=1 p 2 s 2m - 1 2 mω 2 P P s=1 (x s-1 -x s ) 2 + P s=1 1 P V (x s ) (2.46)
The last member of this expression is the average potential energy

U = P s=1 1 P V (x s ) (2.47)
thus the average kinetic energy is given by

K = P s=1 p 2 s 2m - 1 2 mω 2 P P s=1 (x s-1 -x s ) 2 (2.48)
So the kinetic energy can be computed using the following expression called the primitive estimator of the kinetic energy

K prim = P s=1 p 2 s 2m - 1 2 mω 2 P P s=1 (x s-1 -x s ) 2 (2.49)
Since we have equipartition of the energy p 2 s /2m = 1/2β the primitive estimator can also be written

K prim = P 2β - 1 2 mω 2 P P s=1 (x s-1 -x s ) 2 (2.50)
This estimator indeed converges to the exact quantum value of the kinetic energy in the limit P → ∞. But its fluctuations increase with P [4,42] so that it can become inaccurate for high Trotter numbers. Thus another estimator for the average kinetic energy called the virial estimator has been proposed [42]. This estimator is obtained by using a path-integral version of the virial theorem

P 2β - P s=1 1 2 mω 2 P (x s-1 -x s ) 2 = 1 2P P s=1 x s ∂V ∂x s (2.51)
whose demonstration is given in the complements of this chapter -section 2.7. From this, one directly obtains another expression for the kinetic energy called the virial estimator of the kinetic energy

K vir = 1 2P P s=1 x s ∂V ∂x s (2.52)
Finally, this estimator is not translationnally invariant which is a problem for unbounded systems. Thus a generalization valid for both bounded and unbounded systems is generally used [43]

K Cvir = 1 2β + 1 2P P s=1 (x s -x c ) ∂V ∂x s (2.53)
which is sometimes called the centroid-virial estimator of the kinetic energy. With x c the position of the centroid of the ring polymer

x c = 1
P P s=1 x s (2.54)
In figure 2.2 the convergence of those different estimators with the Trotter number P for a simple one-dimensional harmonic oscillator is shown. One clearly see that all the estimators of course give the same average value, and converge to the expected exact quantum value.

All these estimators are for one particle (N = 1) in one dimension. In the more general case of N particles in 3 dimensions the estimators becomes

K prim = 3N P 2β - N i=1 P s=1 1 2 m i ω 2 P ( r i,s-1 -r i,s ) 2 (2.55) K vir = 1 2P N i=1 P s=1 r i,s • ∇ r i,s V (2.56) K Cvir = 3N 2β + 1 2P N i=1 P s=1 ( r i,s -r c ) • ∇ r i,s V (2.57)
One can derive in the same way estimators for other quantities such as the pressure P for example

P = N P βΩ - 2 3Ω N i=1 P s=1 1 2 m i ω 2 P ( r i,s-1 -r i,s ) 2 - 1 3Ω 1 P N i=1 P s=1 r i,s • ∇ r i,s V (2.58) PΩ = 2 3 K prim - 1 3P N i=1 P s=1 r i,s • ∇ r i,s V (2.59)
for an N -body system in 3D and with Ω the volume of the simulation box. Error bars are smaller than the symbols size and thus are not displayed here.

Results on one-dimensional systems

In this section we apply the QTB method on simple one-dimensional systems. Since the quantum thermal bath is only exact in the harmonic case, we are particularly interested to study systems where anharmonicity is important. We first study the case of an harmonic oscillator before switching to the Morse potential with increasing anharmonicity. Finally, we focus on the quartic double-well potential. In each case, we compare the QTB-MD results to PIMD which is known to give the exact results and focus on two important quantities : the average energy and the position density probability (or position distribution).

Harmonic oscillator

Let's first verify that the quantum thermal bath is indeed able to give exact results in the harmonic case. The system is an hydrogen atom in an harmonic potential well

V (x) = 1 2 mω 2 x 2 (2.60)
The temperature T and the angular frequency ω are chosen such that ω/k B T >> 1 more precisely ω/k B T ≈ 13 so the system is mainly in its ground state and T << ω/k B thus quantum effects are expected to be significant. Figure 2.3 shows the probability density of position ρ(x) obtained in the classical case (i.e. with the standard Langevin thermostat), and using QTB-MD and PIMD. As expected the QTB method gives the exact probability density in this case. One can see that the quantum probability density is significantly different than the classical one thus quantum effects indeed are important here. The average energies in reduced units are given in table 2.1. We see that QTB-MD gives the same energy as PIMD and thus gives the exact quantum energy. As expected, the QTB method is exact in the harmonic case.

Morse potential

We now consider the case of an anharmonic potential : the Morse potential

V (x) = D 1 -e -αx 2
(2.61)

The eigenvalues E n for this potential can be expressed in the following form [44]

ε n = E n D = ω D n + 1 2 - 1 λ 2 n + 1 2 2 (2.62)
where ε n are the energies in reduced form and ω = α 2D/m is the angular frequency in the harmonic approximation. The dimensionless parameter λ is defined as

1 λ 2 = 2 α 2 2mD (2.63)
This is a standard potential to model the vibration of diatomic molecules for example. The anharmonicity of the Morse potential can be expressed as the difference between the groundstate energy and its harmonic approximation value, this difference is 1/(4λ 2 ). In this section we will study three different cases with increasing anharmonicity: the first case is weakly anharmonic (1/λ 2 = 0.0015 as for the HCl molecule), in the second case the anharmonicity is increased (1/λ 2 = 0.004) and in the last case the anharmonicity is strong (1/λ 2 = 0.024). The simulations are carried out at a temperature T = 0.02D/k B . Figure 2.4 shows the position probability density for the three cases and table 2.2 gives the average energies in reduced units obtained by QTB-MD and PIMD. Clearly the QTB-MD method becomes approximate when anharmonicity increases. Nevertheless we see that the position distributions obtained with the QTB method are similar to the ones obtained with PIMD. In particular in the weakly anharmonic case the PIMD and QTB-MD results are very close. These results shows that the QTB method is approximate when anharmonicity is strong. This is confirmed by the results on the energy: one can see that the QTB-MD and PIMD energies in the first two cases (1/λ 2 = 0.0015 and 1/λ 2 = 0.004) with low or moderate anharmonicity are very close. But when anharmonicity is strong (1/λ 2 = 0.024) the two methods do not give the same energy anymore.

Quartic double-well

Finally, we study the case of a strongly anharmonic potential: the symmetric quartic doublewell

V (x) = V 0 x a 2 -1 2 (2.64)
The double-well potential appears in many areas of science such as physics, chemistry or biology. It is used in particular to model thermally activated processes such as chemical reactions. The diffusion of point defects in solids is a thermally activated process, and in general, the microscopic mechanisms associated with this diffusion are well described by double well potentials. In the particular case of the hydrogen diffusion in oxide, the proton (H + ) has to overcome a barrier (V 0 ) to jump from one oxygen atom to another. The potential energy associated with this process has the shape of a double well. In the classical case, the proton can only cross the barrier by jumping over, which will only be possible if its energy is sufficient (i.e. if kT V 0 ). On the contrary, in the quantum case, the proton can tunnel through the barrier and thus cross it even if its energy is lower than V 0 . Moreover, because of zero-point energy, the effective barrier height is lower than V 0 and becomes V 0ω/2 in the harmonic approximation. The double-well potential is thus particularly interesting because it allows us to study the behavior of the QTB in cases where tunnelling is possible.

The time-independent Schrödinger equation in the case of the quartic double-well can be written in a dimensionless form

-C d 2 φ dy 2 + (y 2 -1) 2 = εφ (2.65)
where φ(y) are the stationary wavefunctions, y = x/a is the position in a reduced form, ε = E/V 0 is the energy in reduced units and

C = 2 2ma 2 V 0 (2.66)
A numerical resolution of equation (2.65) shows that there is a critical value for the constant C (C 0 = 0.731778) for which the ground state energy is equal to the barrier height V 0 (i.e. ε 0 = 1). Then for C < C 0 , the ground state energy is lower than V 0 (i.e. ε 0 < 1), and for C > C 0 , the ground state energy is higher than the height of the barrier (i.e. ε 0 > 1). In the following we will study three cases. A first one where C = 1 which is greater than C 0 and two cases for which C < C 0 : C = 0.1 and C = 0.3. Figure 2.5 shows the results on the position distribution obtained using PIMD, QTB-MD and Langevin MD.

We clearly see that the quantum thermal bath is approximate in this case. More precisely, the QTB method becomes approximate when the value of the constant C increases and fail to give the correct position distribution for high values of C. Again the quartic double-well potential is interesting here because it allows the investigation of the ability of the quantum thermal bath to describe quantum tunnelling. It is worth noting that there is no reason to expect the QTB method to correctly describe the tunnel effect a priori. Indeed the QTB is only constructed in order to impose the quantum (harmonic) energy distribution while tunnelling is a complex quantum mechanical phenomenon. If the average total energy of the particle E is higher than the barrier height V 0 , the particle will mainly transit from one well to the other through over barrier motion and tunnelling will be negligible. On the contrary, when E < V 0 , tunnelling can become important. In the cases we studied here that corresponds to C = 0.1 and C = 0.3. For C = 0.1, the probability density obtained using the QTB method is similar to that obtained by PIMD. In contrast, for C = 0.3, the QTB method is unable to give the correct position distribution. Tunnelling increases when the barrier height V 0 and/or the inter-well distance 2a decreases which corresponds to an increase of C. Thus quantum tunnelling is especially important for high values of C (with E < V 0 ). So the QTB method gives the correct position distribution in cases where tunnelling effect remains low but becomes approximate when tunnelling is important.

Chapter 2. NQE in MD simulations

C K /V 0 U /V 0 QTB-MD
When the value of C increases, the particle is able to visit more frequently the anharmonic parts of the potential. Since we know that the QTB becomes approximate when anharmonicity increases, it is not very surprising to find that the QTB method fails for high values of C.

Two main features appears when comparing PIMD and QTB position distributions: 1 -the QTB method seems to underestimate the probability of finding the particle at the barrier top and 2 -the QTB method is unable to reproduce the correct position of the two peaks in the probability density : In the exact density the distance between the two peaks decreases when C increases. On the contrary, in the QTB case, the two peaks always corresponds to the bottom of the wells thus the distance between them is constant for any value of C.

Finally, table 2.3 gives the average energy in reduced units obtained using QTB and PIMD in the three cases. We see that, in terms of energy, QTB and PIMD results are comparable for any value of C. In particular, the total energy (K+U ) obtained by the QTB and PIMD methods are basically the same. This can be surprising at first since we found for the highly anharmonic case of the Morse potential that the QTB and PIMD energies were different.

Although, if we analyse the trajectories obtained in this highly anharmonic Morse potential, we find that, in the QTB case, the particle is sometimes able to escape from the well but not on the PIMD case. In contrast, this cannot happen in the double-well case.

Conclusion

When quantum effects are important and needs to be taken into account the PIMD and QTB methods are two possible options. These two methods are different in their formulation and theoretical basis, and they have different advantages and drawbacks. Path integral molecular dynamics gives exact results even for highly anharmonic systems. However, only static quantities can be obtained, and the number of beads required at low temperatures becomes large leading to costly calculations. On the contrary, the computational cost of the QTB method is similar to the computational cost of standard MD simulations. Moreover dynamical quantities are directly accessible. However, we saw that the QTB method is approximate when dealing with highly anharmonic cases. Let us note that, at least in the cases studied here, the QTB-MD results are always closer to the exact quantum results than are the classical MD results. Except in strongly anharmonic systems, the QTB method is expected to give at least qualitatively correct results, and indeed, the method has been able to give satisfactory results in various anharmonic systems [21][22][23][24][25][26][27].

Chapter 2. NQE in MD simulations

Complements of chapter 2 Path integral formulation of the density matrix

We continue to consider the case of one particle of mass m in one dimension (position x) described by the Hamiltonian

Ĥ = K + V = p2 2m + V (x) (2.67)
where K and V are the kinetic and potential energy operators, respectively. The generalization to N particles in three dimensions is straightforward. In the canonical (N, V, T ) ensemble the density operator or density matrix of the system is given by ρ(β) = e -β Ĥ (2.68) whose matrix elements in position representation are

ρ(x, x ′ ; β) = x| ρ(β) |x ′ = x| e -β( K+ V ) |x ′ (2.69)
Now, we would like to split the density matrix in a kinetic and a potential part. Since the kinetic and potential operator do not commute e -β( K+ V ) = e -β K e -β V (2.70) However, we can use the Trotter product formula [45] which allows us 9 to write e -β( K+ V ) = lim which is known as the "primitive" approximation. This procedure is often called Trotter decompositions and P is referred to as the Trotter number. The error associated with this particular approximation is of order P -2 . Other (higher order) Trotter decomposition to go beyond the primitive approximation are sometimes used [46,47].

To simplify notation let us introduce the operator α = e -β P K e -β P V so that the density matrix is ρ(x, x ′ ; β) ≈ x| αP |x ′ (2.74) 9 The Trotter product formula e -τ ( Â+ B) = lim

P →∞ e -τ P Âe -τ P B P
is only valid under several conditions. The two operators Â, B and their sum  + B have to be self-adjoint, moreover the two operators have to be lower bounded. In the case of the kinetic and potential energy operators these conditions are satisfied.

Then, by introducing P -1 times the identity operator 1 = dx |x x|, one obtains the following expression for the canonical density matrix

ρ(x, x ′ ; β) ≈ dx 1 dx 2 ... dx P -1 x| α |x 1 x 1 | α |x 2 ... x P -1 | α |x ′ (2.75)
Now we want to evaluate the term

x i-1 | α |x i = x i-1 | e -β P K e -β P V |x i (2.76)
Since the potential operator V is diagonal in position representation, we can directly write

x i-1 | α |x i = x i-1 | e -β P K |x i e -β P V (x i ) (2.77)
For the kinetic operator, let us introduce the identity operator 1 = dk |k k| in the basis of the eigenvectors |k of the momentum operator

x i-1 | α |x i = dk x i-1 | e -β P K |k k|x i e -β P V (x i ) (2.78)
The kinetic operator is diagonal in momentum space and the eigenvalues are 2 k 2 /2m thus

x i-1 | α |x i = dk x i-1 |k k|x i e -β 2 k 2 /2mP e -β P V (x i ) (2.79)
Now recalling the relation between momentum an position space x|k = 1 √ 2π e ikx (2.80) one finds

x i-1 | α |x i = 1 2π dk e ik(x i-1 -x i ) e -β 2 k 2 /2mP e -β P V (x i ) (2.81)
The integral over k is a Gaussian integral and we obtain after integration

x i-1 | α |x i = mP 2πβ 2 1/2 exp - mP 2β 2 (x i-1 -x i ) 2 - β P V (x i ) (2.82)
Finally, by introducing this expression in equation (2.75), we obtain the following expression for the density matrix

ρ(x, x ′ ; β) ≈ mP 2πβ 2 P/2 dx 1 ... dx P -1 exp - P s=1 mP 2β 2 (x s-1 -x s ) 2 + β P V (x s ) (2.83)
With starting point x 0 fixed at position x 0 = x and ending point x P fixed at x P = x ′ . This is sometimes called a "discrete imaginary time path integral" formulation of the density matrix.

The connection between this expression and path integrals with imaginary time can be understood by first recalling that the time evolution of a quantum system is directly related to the propagator (or evolution operator) Û (t) = e -i Ĥt/ . In particular, if we consider a particle initially prepared in an eigenstate of the position operator |x , the probability amplitude A of finding the particle in the eigenstate |x ′ after the time t is

A = x ′ | Û (t) |x = U (x, x ′ ; t) (2.84)
Then, Feynman path integral formulation of quantum mechanics tells us that this probability amplitude and thus the propagator U (x, x ′ ; t) can be written as a sum over all possible paths to go from x to x ′ in a time t:

U (x, x ′ ; t) = x ′ x Dx(t) e iS[x]/ (2.85) with x ′
x Dx(t) that corresponds to a functional integral with respect to a path x(t) with starting point fixed at position x and ending point fixed at x ′ . Each path contributes by a phase factor exp(iS[x]/ ) with S[x] the classical action:

S[x] = t 0 L(x, ẋ, t) dt = t 0 1 2 m ẋ2 -V (x) dt (2.86)
One can notice the similarity between the density matrix ρ(β) = exp(-β Ĥ) and the propagator Û (t) = exp(-i Ĥt/ ), so that

ρ(β) = Û (-iβ ) (2.87)
thus the density matrix is equivalent to the propagator of the system evolving in imaginary time t ∈ [0, -iβ ]. Since the propagator can be written as a path integral, it is also possible to write the density matrix in a path integral form.

If we go back now to expression (2.83) for the density matrix and take the continuous limit P → ∞ we obtain

ρ(x, x ′ ; β) = lim P →∞ δτ →0 m 2πδτ 2 P/2 dx 1 ... dx P -1 exp -δτ P s=1 m 2 2 x s-1 -x s δτ 2 + V (x s ) (2.88)
where we have used δτ = β/P so that the limit P → ∞ is equivalent to δτ → 0. In this limit, the number of points that connect x to x ′ goes to infinity, and the sequence {x 0 , x 1 , x 2 ...x P -1 , x P } becomes a continuous path x(τ ) from x 0 = x to x P = x ′ . The P -1 integrals now are represented by an integral over the whole function x(τ ) which is called a functional integral. This is formally written using the sign Dx(τ ) so that lim

P →∞ δτ →0 m 2πδτ 2 P/2 dx 1 ... dx P -1 ≡ x ′ x Dx(τ ) (2.89)
Moreover, the argument of the exponential when P → ∞ becomes a standard (Riemann) integral lim

P →∞ δτ →0 δτ P s=1 m 2 2 x s-1 -x s δτ 2 + V (x s ) = β 0 m 2 2 dx dτ 2 + V (x) dτ (2.

90)

Ph.D. thesis -F. Brieuc -2016 Thus the density matrix can be written as follows

ρ(x, x ′ ; β) = x ′ x Dx(τ ) exp - β 0 m 2 2 dx dτ 2 + V (x) dτ (2.91)
Finally, by performing the transformation from τ to imaginary time t = -i τ , one finds the following expression for the density matrix

ρ(x, x ′ ; β) = x ′ x Dx(τ ) exp i -iβ 0 m 2 dx d t 2 -V (x)

PIMD estimators for the kinetic energy

In the following we give more details about the derivation of the standard estimators to compute average values of the kinetic energy in PIMD. We will drop the subscript for PIMD average values so that .. ≡ .. P .

Primitive estimator

Starting from the usual thermodynamic relation for the average energy and using equation (2.43) one can write E as

E = - 1
Z P ∂Z P ∂β (2.94)
which gives using the expression for Z P (eq. (2.41))

E = 1 Z P 1 (2π ) P dp 1 ... dp P dx 1 ... dx P H P + β ∂H P ∂β exp (-βH P ) (2.95) E = H P + β ∂H P ∂β (2.96)
with the expression of H P given in equation (2.42).

H P = P s=1 p 2 s 2m + 1 2 mω 2 P (x s-1 -x s ) 2 + 1 P V (x s )
The derivative of ∂H P /∂β is different from zero since ω P = √ P /β depends on β and we obtain

E = P s=1 p 2 s 2m - 1 2 mω 2 P P s=1 (x s-1 -x s ) 2 + P s=1 1 P V (x s ) (2.97)

Chapter 2. NQE in MD simulations

The last member of this expression is the average potential energy

U = P s=1 1 P V (x s ) (2.98)
thus the average kinetic energy is given by

K = P s=1 p 2 s 2m - 1 2 mω 2 P P s=1 (x s-1 -x s ) 2 (2.99)
The kinetic energy can be computed using the following expression called the primitive estimator

K prim = P s=1 p 2 s 2m - 1 2 mω 2 P P s=1 (x s-1 -x s ) 2 (2.100)
Since we have equipartition of the energy p 2 s /2m = 1/2β the primitive estimator can also be written

K prim = P 2β - 1 2 mω 2 P P s=1 (x s-1 -x s ) 2 (2.101)
As already explained before, this estimator indeed converges to the exact quantum value of the kinetic energy in the limit P → ∞. But its fluctuations increase with P [4,42] so that it can become inaccurate for high Trotter numbers. Thus another estimator for the average kinetic energy called the virial estimator has been proposed [42].

Virial estimators

In PIMD the extended system is described by the effective Hamiltonian H P whose expression is given by (2.42).

H P = P s=1 p 2 s 2m + 1 2 mω 2 P (x s-1 -x s ) 2 + 1 P V (x s )
We define the effective potential x s ∂α P ∂x s + 1

φ P = P s=1 1 2 mω 2 P (x s-1 -x s ) 2 + 1 P V (x s ) = α P + U (2.102) with α P = P s=1 mω 2 P (x s-1 -x s ) 2 /
P P s=1 x s ∂V (x s ) ∂x s (2.109)
And one can easily show that P s=1 x s (∂α P /∂x s ) = 2α P so that 2α P + 1

P P s=1 x s ∂V (x s ) ∂x s = P β (2.110)
So we finally obtain

P 2β - P s=1 1 2 mω 2 P (x s-1 -x s ) 2 = 1 2P P s=1 x s ∂V ∂x s (2.111)
which is a generalization of the virial theorem to the path integral framework (eq. 2.51).

Now from the expression of the primitive estimator (eq. 2.101) and using equation (2.111) we obtain the virial expression of the kinetic energy estimator

K vir = 1 2P P s=1 x s ∂V ∂x s (2.112)
Finally, this estimator is not translationnally invariant which is a problem for unbounded system. Thus a generalization valid for both bounded and unbounded systems is generally used [43]

K Cvir = 1 2β + 1 2P P s=1 (x s -x c ) ∂V ∂x s (2.113)
which is sometimes called the centroid-virial estimator of the kinetic energy, with x c the position of the centroid of the ring polymer. This expression is obtained by developing a similar virial theorem than equation (2.111)

P -1 2β - P s=1 1 2 mω 2 P (x s-1 -x s ) 2 = 1 2P P s=1 (x s -x c ) ∂V ∂x s (2.114)
which can be obtained in the same way as equation (2.111) but starting with calculating

P s=1 (x s -x c ) ∂φ P ∂x s ( 2

Chapter 3

Limitations of the QTB method: zero-point energy leakage

Introduction

In the previous chapter we have seen that, despite being approximate for highly anharmonic systems, the quantum thermal bath method is a good alternative to path integral methods to account for nuclear quantum effects in simulations. Indeed, the method has been able to give satisfactory results in several anharmonic cases [1][2][3][4][5][6][7]. In addition, QTB-MD simulations do not require any additional computational cost compared to standard MD and dynamical properties are directly accessible. However, it has been recently pointed out [8,9] that the QTB method is prone to zero-point energy leakage (ZPEL) in anharmonic systems. ZPEL is a well known problem in methods based on classical trajectories [10][11][12] where part of the energy of high frequency modes are transferred to low frequency ones leading to a wrong energy distribution. In particular, zero-point energy can be transferred leading to some vibrational modes having an energy lower than their zero-point energy. In some cases, this leakage can have dramatic consequences on the computed properties of the system [8,9,13].

In practice, the obtained energy distribution is the result of a competition between the quantum thermal bath and the zero-point energy leakage. On one hand, the QTB gives more energy to high frequency modes than low frequency ones, while on the other hand, energy transfer between the modes tends to homogenize the energy distribution. If the QTB is not able to fully counterbalance the leakage, the obtained energy distribution is intermediate between a completely homogeneous distribution and the expected QTB distribution θ(ω, T ) (eq. 2.18). Zero-point energy leakage occurs in several sytems (such as liquid water, Lennard-Jones systems, diamond crystal...) [8,13,14]. In the particular case of the QTB, it has been first pointed out by Bedoya-Martínez and coworkers in aluminium and argon crystals [8].

In this chapter, we investigate in details the conditions leading to ZPEL and the parameters that influence it. We find, in particular, that increasing the friction coefficient in QTB-MD simulations can significantly reduce ZPEL and even completely remove it. It is worth noting that solutions to the ZPEL problem in QTB-MD simulations has already been proposed [8,15]. In particular, Bedoya-Martńinez et al. proposed to modify the spectral density of the random force in order to counterbalance the leakage. Since the leakage leads to low frequency modes with too much energy and high frequency modes with too low energy, the authors tried to decrease the power spectral density at low frequencies and increase it a high frequencies. This solution has been able to limit the zero-point energy leakage in weakly anharmonic systems. Unfortunately it does not work for strongly anharmonic systems and the procedure is system dependent [8]. Here, we first study a simple model system composed of two one-dimensional harmonic oscillators coupled by anharmonic couplings. This model provides a very clear illustration of the leakage and thus allows the investigation of the parameters driving this phenomenon. Then, we focus on another model system: a one dimensional chain of atoms. This model has been designed to reproduce the basic features of hydrogen bonded systems with realistic orders of magnitude. Here, we focus in particular on the consequences of ZPEL for the computed properties of the system.

This work is the result of a collaboration with Y. Bronstein, P. Depondt, F. Finocchi (INSP, UPMC) and M. Hayoun (LSI, École Polytechnique) and has been accepted for publication in the Journal of Chemical Theory and Computation.

Coupled harmonic oscillators

The first model we study here is composed of two harmonic oscillators coupled by cubic or quartic anharmonic couplings. This system is interesting because it provide a very simple description of vibrational modes coupled by anharmonic coupling terms. The first oscillator is associated with the high frequency modes and the second to low frequency modes (ω 1 > ω 2 ). The Hamiltonian of the system writes

H = 1 2 m ẋ2 1 + 1 2 mω 2 1 x 2 1 + 1 2 m ẋ2 2 + 1 2 mω 2 2 x 2 2 + C 3 (x 1 -x 2 ) 3 + C 4 (x 1 -x 2 ) 4 (3.1)
where x 1 , x 2 are the positions of oscillators 1 and 2 respectively; ω 1 , ω 2 are the angular frequencies of oscillators 1 and 2; m is their mass and C 3 , C 4 are the coupling constants associated with cubic and quartic couplings, respectively. This Hamiltonian can be written in a dimensionless form H = H/ω 1 , so that H = q2

1 2 + q 2 1 2 + q2 2 2 + Ω 2 q 2 2 2 + c 3 (q 1 -q 2 ) 3 + c 4 (q 1 -q 2 ) 4 (3.2)
where the following change of variables has been used:

Ω = ω 2 ω 1 = ν 2 ν 1 q i = x i ξ ξ = mω 1 qi = dq i dt * t * = ω 1 t c 3 = C 3 ξ 3 ω 1 c 4 = C 4 ξ 4 ω 1 (3.3)
with q 1 and q 2 being the reduced position of the two oscillators and Ω the ratio of their frequencies. We study here the influence of the coupling constant c 3 and c 4 on the average energy of the oscillators

ε 1 = q2 1 2 + q 2 1 2 , ε 2 = q2 2 2 + Ω 2 q 2 2 2 . ( 3.4)
The QTB-MD simulations are performed with a friction coefficient γ = 4 × 10 -4 ω 1 , an angular frequency cut-off ω cut = 2ω 1 and a time step δt = 0.05ω -1 1 . Average values are obtained using 30 independent trajectories of 10 compared to the exact results obtained by numerically solving Schrödinger equation1 which is possible here because of the small number of degrees of freedom. The ratio Ω is varied in the range 0.05 -0.8 and the values of the coupling constants c 3 and c 4 lies in the range 0-25×10 -4 and 0-40×10 -4 respectively, so that we cover a large range of coupling energies while remaining in a rather weak coupling regime. The temperature is set to k B T = 0.03 ω 1 in order to ensure that the thermal energy contribution of the two oscillators remains small compared to their zero-point energies. Indeed, the exact calculation shows that the energies of the oscillators are basically equal to their zero-point energies (ε 1 = 0.5 and ε 2 = Ω/2) and are almost independent of the coupling constants c 3 and c 4 for the range of values studied here. Figure 3.1 shows the values of ε 1 and ε 2 obtained by QTB-MD as a function of c 3 and c 4 for two distinct cases: the case of cubic coupling only (c 4 = 0) and Ω = 0.5 (top panel) and the case of quartic coupling only (c 3 = 0) and Ω = 0.25 (bottom panel). As expected, in the uncoupled case (i.e. c 3 = c 4 = 0 ) the QTB method gives the exact results. However, when the coupling constant c 3 or c 4 increases, the QTB-MD energies diverge from the exact results. More precisely, the energy of oscillator 1 is underestimated and the energy of oscillator 2 is overestimated: a part of the energy from the high frequency oscillator has been transferred to the low frequency oscillator. This is a clear illustration of the ZPEL problem in QTB-MD simulations.

This simple model only depends on 3 independent parameters Ω, c 3 and c 4 and we will now investigate their effect on ZPEL. In order to quantify the leakage, we define the following

ζ = ∆ε exact -∆ε QTB ∆ε exact = (ε exact 1 -ε exact 2 ) -(ε QTB 1 -ε QTB 2 ) ε exact 1 -ε exact 2 (3.5)
In this definition, the leakage is maximum for ζ = 1: the system has reached an equipartition of the energy

ε QTB 1 = ε QTB 2
. On the contrary, there is no leakage for ζ = 0 since in this case ε QTB 1,2 = ε exact 1,2 . Figure 3.2 presents the evolution of ζ as a function of Ω in the cubic case for different values of c 3 (top panel) and in the quartic case for different values of c 4 (bottom panel). We see that the zero-point energy leakage is highly sensitive to the value of Ω. Clearly, in the cubic case, ZPEL only occurs for Ω ≈ 0.5. This behavior can be understood by the fact that cubic couplings are responsible for frequency doubling i.e. the creation of additional modes at 2ω. Indeed additional peaks appears in the vibrational spectrum (see figure 3.3 (a)) at 2ω 2 , ω 1ω 2 and ω 1 + ω 2 . The leakage occurs in the particular case of Ω = 0.5 for which there is a clear resonance between the modes 2ω 2 and ω 1 and between the modes ω 1ω 2 and ω 2 . Thus ZPEL seems to be directly related to resonances between different vibrational modes in the system. The quartic case however is more complicated. We expect significant ZPEL for Ω ≈ 1/3 since quartic couplings are responsible for the generation of modes with frequency 3ω. There should be, in particular, a resonance between the modes at 3ω 2 and ω 1 for Ω = 1/3. Significant ZPEL is indeed observed at Ω ≈ 1/3 which confirm the link between the leakage and the resonances. However, we note that important leakage is also observed for Ω < 1/3. arbitrary units For now, all the QTB-MD simulations have been carried out with a fixed value of the friction coefficient γ = 4 × 10 -4 ω 1 . This value has been chosen to be small compared to the lowest angular frequency (ω 2 ) in order to ensure that the thermostat does not perturb the natural dynamics of the system as it is explained in the first chapter. The friction coefficient γ is directly related to the strength of the coupling between the system and the thermostat. Thus one can expect that, if γ increases, the QTB will be able to better counterbalance the leakage. And, for a value of γ that is high enough, the QTB should be able to fully counterbalance the leakage and thus impose the correct energy distribution. In the following, the friction coefficient γ is varied from 4 × 10 -4 ω 1 to 2 × 10 -2 ω 1 and we focus on the cases where ZPEL is important i.e Ω = 0.5 in the cubic case and Ω = 0.25 in the quartic case. Figure 3.4 shows that indeed the leakage strongly depends on the value of the friction coefficient and increasing γ reduces the ZPEL. This is particularly clear in the cubic case for which the increase of γ can completely remove the leakage. Figure 3.3 (c) shows the vibrational spectra obtain with a higher value of γ in the cubic case; we see that the additional modes corresponding to the resonances have disappeared which further confirms the relation between ZPEL and mode resonances. In the quartic case however, increasing γ is not sufficient to completely remove the leakage. Figure 3.3 (d) shows the vibrational spectra in this case and indeed increasing γ is not sufficient to suppress all the additional peaks in the spectrum. Thus resonances are still possible and ZPEL remains.

γ=4x10 -4 ω 1 ω 1 -ω 2 2ω 2 ω 1 +ω 2 (a) 0 0.5 1 1.5 2 γ=4x10 -3 ω 1 (c) 0 0.5 1 1.5 2 arbitrary units ω/ω 1 γ=4x10 -4 ω 1 (b) 0 0.5 1 1.5 2 ω/ω 1 γ=4x10 -3 ω 1 (d)
In conclusion, we have seen that ZPEL is directly related to resonances between vibrational modes. In a realistic system composed of a large number of degrees of freedom these mode resonances are of course very likely to occur and will be almost impossible to avoid in practice. However, we have found that increasing the friction coefficient can largely help to limit the impact of ZPEL in particular in the case of cubic coupling terms. Of course realistic interaction potentials will consist of all the higher-order terms (not just cubic or quartic terms) and thus the question whether the increase of γ will still be helpful in realistic systems arises. For this reason, we study in the next section another model system with realistic interaction potentials (Morse type) between the atoms. Finally, increasing γ raises the issue that usually, in standard Langevin dynamics, the friction coefficient should remain low enough in order not to perturb the dynamics. A high value of γ will for example results in important broadening of the peaks in the vibrational spectrum. In the next section, we focus in particular on the effect of the leakage and of the increase of γ on the properties of the system.

One-dimensional chain of atoms

We now study a one dimensional chain of atoms O-H-O-H-O-H composed of N O = 3 oxygen atoms and N H = 3 hydrogen atoms with periodic boundary conditions. The O-H interaction is described using a Morse type potential proposed by Johannsen [16] that has the following form:

V OH (r) = u 0 a + be a(r-r 0 ) a e -b(r-r 0 ) -1 + b e a(r-r 0 ) -1 -u 0 (3.
V OO (R) = C 0 1 -e -α 0 (R-R 0 ) 2 -C 0 (3.
ν OO = 1 2π 4C 0 α 2 0 m O (3.8)
with m O the mass of the oxygen atom. The value of α 0 is varied so that ν OO lies between 10 and 70 THz. Clearly, this model cannot be used to accurately represent hydrogen-bonded systems such as ice or liquid water for example. However, it exhibits the basic features of those type of systems with realistic orders of magnitude i.e. high-frequency modes related to O-H vibrations (with frequencies around 100 THz) and low frequency modes related to O-O vibrations (with frequencies of several tens of THz).

The potential energy of an hydrogen atom in the O-H-O group is given by V OH (r)+V OH (R-r) which has a double-well potential shape. A normal mode analysis of this system reveals one low frequency mode (of frequency ν 2 ) associated with O-O vibrations and two high frequency modes (of very similar frequencies ν 1 ) that corresponds to O-H stretching modes (see figure 3.5). In analogy with the previous model, the O-H stretching modes can be considered to play the role of the high frequency oscillator while the O-O lattice mode corresponds to the low frequency oscillator. In the following, we are first interested to study the influence of the parameter Ω = ν 2 /ν 1 and the friction coefficient γ on the ZPEL. We want in particular to see if the results obtained for the two oscillators model are confirmed in this system.

The QTB-MD simulations are carried out at T = 600 K with a time step δt = 0.1 fs and average values are obtained using 12 independent trajectories of 3 ns each. The frequency ν 2 is varied through the parameter α 0 while ν OH is kept constant (thus ν 1 remains almost constant). The QTB-MD results are compared to PIMD results obtained with a Trotter number P = 20 which ensures a good convergence in all the cases studied here.

To study the leakage in this case, we choose to compare the kinetic energy of the light atoms (H) which are mainly involved in high frequency modes to the kinetic energy of the heavier atoms (O) which mainly participate in low frequency modes. In order to do this we define the effective temperatures (T H and T O ) of hydrogen and oxygen atoms from their kinetic energies:

k B T H 2 = 1 N H N H i=1 K i , k B T O 2 = 1 N O N O i=1 K i (3.9)
where K i is the average kinetic energy of atom i. In standard MD simulations, the equipartition of the energy results in T H = T O since kinetic energy is equally distributed among all the modes. However, in a quantum system (i.e. for PIMD and QTB-MD simulations), the average kinetic energy is higher for high frequency modes than for low frequency ones and therefore T H > T O . Moreover, the modes have higher kinetic energy than in the classical case (mainly because of their zero-point energy) and thus T H and T O are expected to be higher than T = 600 K. Figure 3.6 shows the values of these effective temperatures during the simulation. We see that the QTB-MD simulations tend to underestimate T H and overestimate T O : this a signature of ZPEL. As expected, the leakage tends to increase the kinetic energy of low frequency modes and decrease the kinetic energy of high frequency ones. Now in order to quantify the leakage we define the following deviation factor

ζ = (T H -T O ) (PIMD) -(T H -T O ) (QTB) (T H -T O ) (PIMD) (3.10)
Similarly to the previous model, ZPEL is maximum for ζ = 1 and there is no leakage for ζ = 0. The evolution of ζ as a function of Ω is represented in figure 3.7. We see that ZPEL is strongly dependent on Ω and that, as in the previous model, significant leakage mostly arises for Ω ≈ 1/2. This is reminiscent of the resonances effect that we have seen on the previous model. However, we also observe important ZPEL for Ω ≈ 0.1: this case corresponds to a highly anharmonic regime where a structural transition occurs and therefore corresponds to a different physical situation than for the other values of Ω. This transition is characterized by a splitting of the O-O distance as it can be seen in figure 3.8 which presents the distribution of distances: two different equilibrium O-O distances are observed for Ω = 0.1 (figure 3.8 (d)). We did not study this transition in details here however we now focus on the cases for which ZPEL is strong i.e. Ω = 0.1 and Ω = 0.5. Figure 3.7 also shows that increasing γ can indeed tremendously help to limit ZPEL in this case: for γ = 20 THz, zero-point energy leakage is lower than 0.2 for any value of Ω, even Ω = 0.1 and Ω = 0.5.

We now study the impact of the ZPEL and of the increase of γ on the properties of the system. First, we focus on the distribution of distances presented in figure 3.8. We see that, in the case of Ω = 0.5 (top panel), the impact of ZPEL on the O-H distance (figure 3.8 (a)) is relatively small. In contrast, the effect of the leakage is more substantial on the distribution of O-O distances (figure 3.8 (b)). In particular, the distribution is too broad because of the excess of kinetic energy coming from the ZPEL. Increasing the friction coefficient γ clearly eliminates the effects of the leakage: for a value of γ = 20 THz the QTB and PIMD distributions coincides. In the case of Ω = 0.1 (bottom panel), we see that the impact of ZPEL on the distribution of distances seems to be more important compared to the Ω = 0.5 case. However, one should be careful with the analysis of this case since it corresponds to a highly anharmonic case. Indeed, the QTB method is known to be approximate in strongly anharmonic cases even if no leakage occurs. The increase of γ here seems to be less efficient since there still is an important difference between the PIMD and the QTB distributions for γ = 20 THz (figure 3.8 (c) and (d)). It is worth noting that for a very high value of γ (200 THz) the QTB and PIMD distributions almost coincide. Thus it seems that the major part of the failure of the QTB in this case is indeed coming from ZPEL and the impact of the leakage can also be removed by increasing γ. However, this value of γ is higher than the typical frequencies in the system, thus dynamical quantities cannot be used in this case. We enter a regime where the QTB-MD simulations are overdamped and thus the thermostat is now only used to sample the phase space and the dynamics of the system is strongly perturbed by the thermostat. As already said, one normally should decrease and not increase γ in Langevin dynamics and thus increasing γ to very high values will have consequences in particular on dynamical properties. Figure 3.9 shows the vibrational spectra obtained from QTB-MD trajectories for two values of γ in the case of Ω = 0.5. For the low value of γ (0.2 THz) ZPEL is present in the system while for the higher value of γ (10 THz) the leakage is almost completely removed. As expected, increasing γ leads to a broadening of the peaks however the positions of these peaks and thus the values of the mode frequencies remains unchanged with the increase of the friction. Thus, even with a rather large value of γ, vibrational spectra computed using QTB-MD contain some useful information.

Practical discussions

We have seen that the ZPEL is related to resonances between vibrational modes. These resonances will generally be unavoidable in realistic systems with many degrees of freedom. However, we have seen that increasing the friction coefficient significantly reduces the leakage. Thus, before launching QTB-MD simulations to compute average values, one should first try to evaluate the impact of ZPEL for the particular system under study. And then choose the value for the friction coefficient γ that allows one to sufficiently reduce the leakage. For systems containing different chemical elements the ratio of their kinetic energies can be used to quantify the leakage. Another more general method is to compute the kinetic energy distribution among the modes and compare it to the expected distribution θ(ω, T ). This is the procedure that Bedoya-Martínez and coworkers have used in their paper (figure 1 and 5 in ref. [8]) to illustrate the ZPEL within QTB-MD simulations (see also figure 3.11 in the complements of this chapter). In order to choose the value of the friction coefficient one can for example monitor the evolution of the kinetic energy distribution with γ. One can also compute different properties of interest using QTB-MD and study their evolution with the increase of the friction coeffcient. When it is possible comparing the results of QTB-MD method with results from PIMD is very useful to assess the impact of ZPEL. The PIMD calculations could be done on a smaller simulation box before launching QTB-MD simulations on larger cells once a good value of γ is found.

High values of γ will have consequences on the vibrational properties such as the vibrational spectrum. At the end, one needs to find the best compromise between the impact of ZPEL and the quality needed for the spectrum. It is worth noting that, vibrational spectra generally are computed in the NVE ensemble using standard (classical) MD. Thermostatted MD simulations are used to generate initial conditions and then the spectrum is computed using NVE trajectories starting from these thermostated initial conditions [3,9]. This procedure allows one to completely remove the impact of the thermostat when computing the spectrum and thus one could try to use it with QTB-MD simulations. However this procedure should not be used to compute spectra using the QTB method. Indeed, one could generate initial conditions using QTB-MD, possibly with a high friction coefficient to limit the effect of ZPEL, and then remove the thermal bath to compute the spectrum. However, once the thermal bath is removed the system follows a classical dynamics and thus will naturally reach an equipartition of the energy in a similar way as for the ZPEL. Thus one should not use NVE trajectories to compute vibrational spectra.

It is worth noting that increasing γ can also have consequences on the energy of the system.

In particular, high values of γ can lead to the divergence of the total energy as already emphasized by Barrat and Rodney [17]. The top panel of figure 3.10 shows the relative deviation of the kinetic and potential parts of the QTB-MD energies from their exact values as a function of γ and for an harmonic oscillator. We clearly see that for high values of γ, the energies computed using QTB-MD diverges from the exact values. In particular, the kinetic part of the energy is very sensitive to the increase of γ. In contrast, the potential part of the energy is less sensitive to the value of the friction coefficient and diverges only for high values of γ (for γ > 0.1ω 0 i.e. γ > 0.6ν 0 ). The average kinetic and potential energies obtained using the QTB method in the case of an harmonic oscillator can be expressed in the following form:

K = K(ω) dω 2π = 1 2 m|v(ω)| 2 dω 2π (3.11) V = V (ω) dω 2π = 1 2 mω 2 0 |x(ω)| 2 dω 2π (3.12)
with K(ω), V (ω) the Fourier transform of the kinetic and potential energy respectively and x(ω), v(ω) the Fourier transform of the position and the velocity which can be computed from the Fourier transform of the equation of motion:

x(ω) = -1 m ω 2 -ω 2 0 + iγω R(ω) (3.13) v(ω) = iω m ω 2 -ω 2 0 + iγω R(ω) (3.14)
Thus we find that the kinetic and potential energies writes:

K = K(ω) dω 2π = γ ω 2 ω 2 -ω 2 0 2 + γ 2 ω 2 θ(ω, T ) dω 2π (3.15) V = V (ω) dω 2π = γ ω 2 0 ω 2 -ω 2 0 2 + γ 2 ω 2 θ(ω, T ) dω 2π (3.16)
K(ω) and V (ω) are plotted for different value of γ on the bottom panel of figure 3.10 we first see that increasing γ leads to a broadening of K(ω) and V (ω). Moreover, for high values of γ, both K(ω) and V (ω) become asymmetric and their maximum is shifted from ω 0 . For the potential part V (ω), the tail for ω < ω 0 is longer than the tail for ω > ω 0 . In addition, the maximum is displaced towards lower values than ω 0 . This behavior explains why the potential energy gives values lower than the exact ones for high γ. In contrast, for the kinetic part K(ω), the tail for ω > ω 0 is longer than the tail for ω < ω 0 and the maximum is displaced towards higher values than ω 0 . Thus explaining why the kinetic energy tends to be overestimated for high values of γ. Moreover, we see that the asymmetry is stronger for K(ω) than for V (ω) which explains why the kinetic energy is more sensitive to the increase of γ than the potential energy. We also see, in figure 3.10, that reducing the angular frequency cut-off ω cut allows to reduce the divergence of the kinetic energy. The effect of the cut-off is to remove a part of the high frequency modes thus compensating for the asymmetry of K(ω) and leading to a decrease of the kinetic energy. Thus, decreasing ω cut limit the divergence of the kinetic energy. In contrast, decreasing ω cut does not impact the divergence of the potential energy since the cut-off only acts on the high frequencies. Let us note that ω cut can not be chosen arbitrarily small and should of course remain higher than all the frequencies in the system. In practice, ω cut should be equal to a few times the highest angular frequency. From figure 3.10 we see that ω cut = 2ω 0 or lower seems to be a good choice which avoid any divergence of the energy up to γ ≈ 0.05ω 0 ≈ 0.3ν 0 . Thus, we generally choose ω cut = 2ω max or lower in our QTB-MD simulations (with ω max being the highest angular frequency in the system). Finally, it is worth noting that the kinetic energy can also be computed using the virial theorem so that:

K vir = - 1 2 N i=1 r i • f i = 1 2 N i=1 r i • ∂V ∂ r i (3.17)
with f i the internal force that applies on the atom i and r i its position. Of course, for one degree of freedom the expression simply becomes K vir = -1/2 f • r . Since this expression is based on the potential energy its behavior with increasing γ is the same as the potential energy (figure 3.10 top panel). For this reason, computing the kinetic energy in this way generally is more robust than using the velocities, in the sense that it is independent of the choice of ω cut . Moreover the virial expression leads to better estimation of the kinetic energy if ω cut has been chosen too high. It is thus preferable to compute the kinetic energy both using the velocities and the virial theorem in QTB-MD simulations and compare these two values which should be equal if there is absolutely no divergence of the energy.

Conclusion

In this chapter we have studied in details the ZPEL problem in QTB-MD simulations. We first found that ZPEL is directly related to resonances between vibrational modes. These resonances will hardly be avoidable in any realistic system. However, we also found that increasing the friction coefficient γ can significantly reduce the leakage and even completely remove it in some cases. This can be understood by recalling that the QTB tries to ensure an energy distribution in which high frequency modes have more energy than low frequency ones while MD tends to homogenize the energy among the modes. The ZPEL comes from the energy transfer between the modes resulting in an intermediate energy distribution between the QTB distribution (θ(ω, T )) and a completely homogeneous one. When the friction coefficient is increased the system becomes more strongly coupled to the thermal bath and thus the QTB is able to better counterbalance the leakage. For a sufficiently high value of γ the QTB ensures the correct energy distribution.

ZPEL can have important consequences on the energy distribution among the modes and on the computed properties such as the distribution of distances. In all the systems studied here, increasing γ have been able to limit the impact of ZPEL for the structural properties.

In some cases, such as in the one-dimensional chain of atoms with Ω = 0.1, a very high value of γ was necessary. In this case, vibrational properties are significantly affected by the friction coefficient and thus can not be used anymore. However, we have seen that, even though increasing γ implies an important broadening of the peak in the vibrational spectrum, it seems that useful informations can still be obtained for rather large values of γ.

In particular, mode frequencies are not dramatically affected by the increase of the friction. We expect the low frequency part of the spectrum to be much more altered by a high value of γ than the high frequency part. Thus, if we focus for example on OH stretching and bending modes in hydrogen-bonded materials, it should be possible to use a relatively high value of the friction coefficient since the frequencies of these modes will usually remain larger than γ.

Increasing γ seems to be a simple and yet effective way to remove or at least limit the impact of zero-point energy leakage in QTB-MD simulations. This has been further confirmed on two different realistic systems: a model crystal of aluminium and a ferroelectric crystal (BaTiO 3 ). On both systems, increasing γ has proved to be an efficient method to limit ZPEL. The results are presented in the complements of this chapter. In strongly anharmonic cases, the value of γ necessary to remove the leakage might be of the same order of magnitude or even higher that the typical frequencies in the system. In this case, the computed vibrational properties becomes unusable and we enter an overdamped regime. However, it seems that overdamped QTB-MD simulations might still be able to provide valuable results for structural properties (see the 1D chain of atom with Ω = 0.1 for example).

60

Chapter 3. Limitations of the QTB method: ZPEL

Complements of chapter 3

The results briefly presented here have been obtained by H. Dammak and M. Hayoun.

Aluminium

The ZPEL problem in QTB-MD simulations have been first pointed out by Bedoya-Martínez and coworkers [8] on a perfect crystal of aluminium modelled by a Lennard-Jones potential

V (r) = 4ε σ r 12 - σ r 6 (3.18)
with ε = 0.125 eV, σ = 2.54 Å and a cutoff of 1.37σ. In their paper, Bedoya-Martínez et al. showed that the kinetic energy distribution obtained by QTB-MD simulations at T = 10 K and using a friction coefficient γ = 1 THz is not the expected QTB distribution θ(ν, T ). Indeed, the leakage leads to an energy flow from high frequency modes to low frequency ones and the QTB is unable to fully counterbalance this leakage. We confirm this results for γ = 0.9 THz as we can see in figure 3.11 (full circles). The energy of low frequency modes is overestimated while the energy of high frequency modes is underestimated because of the leakage. However, increasing γ allows to remove the ZPEL. Indeed, the energy distribution obtained by QTB-MD when γ is increased to 10 THz is very close to the correct QTB distribution θ(ν, T ) (open circles in figure 3.11). The inset of figure 3.11 shows the convergence of the QTB-MD energy distribution with γ, we see that for γ 9 THz the leakage is removed.

Barium Titanate

Barium titanate (BaTiO 3 ) is a ferroelectric material described by a complex energy landscape. More details on this material are given in the next chapter. In summary, BaTiO 3 exhibit a complex sequence of phase transitions with increasing temperature: from a low temperature rhombohedral (R) phase to an orthorhombic (O) phase to a tetragonal (T) phase and finally to a high temperature cubic (C) phase. In the first three phases (R,O and T) the system is ferroelectric i.e. it exhibit a macroscopic electric polarization. Quantum effects in this material has been shown to significantly decrease the transition temperatures by ∼ 30 -50 K [18] and to strongly modify the shape of the pressure-temperature phase diagram [19,20]. The system is modelled here using an effective Hamiltonian proposed by Zhong et al. [21,22] and derived from density functional theory calculations. QTB-MD and PIMD simulations have been performed in the isothermal-isobaric ensemble (using a Langevin barostat [23]) with different values of the friction coefficient γ for QTB-MD simulations and for temperatures ranging from 1 K to 270 K. Figure 3.12 presents the values of the polarization as a function of the temperature. The different phases are visible and transition temperatures are indicated by dashed lines. We see that in the cubic phase the polarization is zero since the phase is not ferroelectric. The QTB method with a small value of γ (0.1 THz) fails to reproduce the correct sequence of phase transitions (see the inset of figure 3.12). More precisely, only tetragonal and cubic phases are predicted (i.e. the system is predicted to be in the tetragonal phase even at T = 1 K). This failure is attributed to ZPEL. Indeed, when γ is increased, the complete sequence of phase transitions is retrieved. Moreover, for γ ≈ 16 THz the effect of ZPEL is suppressed and we see in figure 3.12 that we obtain the correct phases along with transition temperatures that are very similar to those obtained using PIMD, in this case. Chapter 4

Quantum thermal bath for path integral molecular dynamics simulation

The work presented in this chapter has been published in [1]. The structure of the chapter and the majority of the figures it contains are based on this publication

Introduction

Path integral molecular dynamics allows one to account for nuclear quantum effects in simulations, and gives exact results even for strongly anharmonic systems. But computational cost is an important limitation of the method. In particular, at low temperatures, the number of replicas necessary to reach the convergence becomes high, and thus the computation time becomes prohibitive. For this reason, PIMD simulations on large systems are difficult to achieve, in particular using first-principles description of the interatomic forces. Several methods have been designed to reduce the computational cost of standard PIMD. Among them are methods based on higher-order Trotter discretisation [2,3] or ring polymer contraction schemes where the slowly varying (long-range) contributions of the interaction potential are treated with fewer beads than the rapidly varying (short-range) contributions [4] and which have been recently generalized for ab initio simulations [5][6][7].

The quantum thermal bath is a good alternative to PIMD that allows one to include nuclear quantum effects at no additional cost compared to standard MD. Moreover time dependent quantities are directly accessible. However, the method becomes approximate when dealing with strongly anharmonic systems and, more problematic, is prone to zero-point energy leakage. In this section we present another way to use the QTB method. More precisely, we show how the quantum thermal bath can be combined with path integral molecular dynamics in order to decrease the computational cost of PIMD simulations [1]. In this combined QTB-PIMD method, the classical thermostat of standard PIMD simulations is replaced by the quantum thermal bath. A part of the quantum effects is introduced through the QTB leading to a faster convergence with the number of beads than standard PIMD. Fewer replicas are needed to converge so the computation time is decreased. In this chapter we first describe the procedure to combine the quantum thermal bath with PIMD. Then we compare the QTB-PIMD method to standard PIMD in several test cases.

Combining the quantum thermal bath with path integral molecular dynamics

This section describes in details the procedure to combine the QTB with PIMD. This procedure has been first proposed by Ceriotti and coworkers [8] to combine their GLE method with path integral molecular dynamics. Here we show that the combination procedure can been adapted to the QTB case. The main idea is to use the quantum thermal bath as a thermostat for path integral molecular dynamics simulations. In order to achieve this, the QTB method needs to be adapted to the PIMD case. The quantum thermal bath is basically defined by the power spectral density of the random force. Thus the spectral density needs to be adapted to the PIMD case, in particular it will now depend on the number of replicas used during the simulation.

Modified power spectral density

In both PIMD (with a Langevin thermostat) and the combined QTB-PIMD methods the equations of motion are the same. Considering, for simplicity, the case of one particle (N =1) of mass m moving in one dimension, the equation of motion for one bead s of the polymer is given by equation (2.44) and writes ṗs = f smω 2 P (2x sx s+1x s-1 )γp s + R s where f s = -(1/P )∂V /∂x s is the internal force exerted on replica s and the two last terms are the frictional (friction coefficient γ) and stochastic forces of the thermostat, respectively. As explained in the first part (see section 2.3), the random force is defined by its power spectral density I R . This spectral density is obtained from the fluctuation-dissipation theorem which can be written in the following form (equation (2.24))

I R (ω, T ) = 2mγκ(ω, T ) (4.1)
In standard PIMD, the thermostat applied on each bead of the polymer is classical and thus κ(ω, T ) = k B T . In the case of the QTB-PIMD method, the spectral density of the force depends on ω but also on the number of beads P used in the simulation. When P = 1 we are in the case of the bare QTB method, and thus we know that κ(ω, T ) = θ(ω, T ) which is the average energy of a quantum harmonic oscillator. Now we need to determined κ(ω, T ) for P > 1.

For this, let us consider the case of an harmonic oscillator:

V (x) = 1 2 mω 2 x 2 .
The position fluctuation for this system is given by

x 2 = θ(ω, T ) mω 2 = mω 1 2 + 1 exp(β ω) -1 (4.2)
which can be rewritten

x 2 = 2mω coth β ω 2 (4.3)
In PIMD one can transform the coordinates x s (s = 1, ..., P ) into normal modes coordinates q k (k = 0, ..., P -1) of the ring polymer. In the simple harmonic case, these normal modes are independent harmonic oscillators of angular frequency ω k given by

ω 2 k = ω 2 P + 4ω 2 P sin 2 kπ P (4.4)
These two coordinate systems are formally related to each other through (discrete) Fourier transform so that

x s = 1 √ P P -1 k=0
q k e i2πks/P (4.5)

And the position fluctuation of one bead s is given by

x 2 s = 1 P P -1 k=0 q 2 k (4.6)
If the quantum thermal bath (with a spectral density of equation (4.1)) is applied on the normal mode coordinates their average kinetic energy is κ(ω, T )/2 and their position fluctuation is

q 2 k = κ(ω k , T ) mω 2 k (4.7)
Now, if the QTB method is correctly applied to the system, each bead should have the correct quantum position fluctuation 1 (equation (4.3)) so that

x 2 s = 1 P P -1 k=0 q 2 k = 2mω coth β ω 2 (4.8)
Then using equation (4.7) we obtain the following relation

1 P P -1 k=0 κ(ω k , T ) mω 2 k = 2mω coth β ω 2 (4.9)
We now need to determine the function κ(ω, T ) solution of this equation. Defining the following dimensionless quantities

u = β ω 2 h(u) = u coth(u) f (0) P (u) = β P κ 2u β (4.10)
Equation (4.9) can be rewritten in the following dimensionless form

P -1 k=0 u 2 u 2 k f (0) P (u k ) = h(u) (4.11)
with u k the reduced angular frequency which, according to equation (4.4) and (4.10), reads

u 2 k = u 2 P + P sin 2 kπ P (4.12)
Equation (4.11) can be numerically solved using the self-consistent iterative technique proposed in [8] and described in the complements of this chapter (section 4.6). Then the function κ(ω, T ) can be obtained from f where U = s V (x s )/P . In the path integral framework, the kinetic energy is generally computed using the centroid virial estimator (eq. (2.53)) which writes

K = 1 2β + 1 2P P s=1 (x s -x c ) ∂V ∂x s where x c = 1 P P s=1
x s is the position of the centroid. In the harmonic oscillator case, the centroid virial estimator can be rewritten by replacing V (x) by its expression

K = 1 2β + 1 P P s=1 1 2 mω 2 x 2 s - 1 2 mω 2 x 2 c (4.15) K = 1 2β + U - 1 2 mω 2 x 2 c (4.16)
Of course in the harmonic case K = U thus we should have

1 2 mω 2 x 2 c = 1 2β (4.17)
This relation holds in standard PIMD since the dynamics of the beads is classical but not when the quantum thermal bath is applied. Thus if we want to compute the kinetic energy using the centroid virial estimator we need to ensure this relation. In ref. [9] Ceriotti and coworkers proposed to do so by treating the k = 0 normal mode classically. Indeed the k = 0 normal mode is related to the centroid and in particular q 2 0 = P x 2 c . Moreover, its angular frequency is ω 0 = ω/ √ P thus one finds that

1 2 mω 2 x 2 c = 1 2 mω 2 0 q 2 0 (4.18)
So we see that the centroid mode has to follow a classical dynamics if we want expression (4.17) to hold. This can be achieved by applying a classical thermostat to the k = 0 normal mode, using the other k > 0 normal modes to impose the correct quantum fluctuations.

In this case, equation (4.9) is modified as follows

1 P k B T mω 2 0 + 1 P P -1 k=1 κ(ω k , T ) mω 2 k = 2mω coth β ω 2 (4.19)
Using the same dimensionless quantities defined in equations (4.10) we obtain

P -1 k=1 u 2 u 2 k f (1) P (u k ) = h(u) -1 (4.20)
Ph.D. thesis -F. Brieuc -2016 where f (0)

P , f (1)
P are used to differentiate between the two cases. The f (0) P function is associated with the case where all the normal modes are treated in the same way and the f (1) P with the case where the k = 0 normal mode is classically considered. Equation (4.20) can also be solved using the iterative self consistent technique described in the complements of this chapter (section 4.6). Note that all the masses m in the equations are the physical masses and thus the combination requires the use of physical bead masses and not fictitious ones. • Then we generate the random forces R k (t) with a power spectral density given by equation (4.1) :

I R (ω, T ) = 2mγκ(ω, T )
using the same method as for the standard QTB method. Of course, in the f

(1)

P case, the random force associated with the k = 0 normal mode is classical and thus generated with a power spectral density equal to 2mγk B T . Let us note that the random forces R k are the forces to be applied on the normal modes of the ring polymer. In the f (0) P case, one can show that the power spectral density that applies to the beads or to the normal modes are equivalent and thus the random forces R k can be directly applied to the beads. But this is not the case with the f

(1) P function since the centroid mode (k = 0) is not addressed in the same way as the others. In this case one needs to transform the random forces R k that applies to the normal modes to the random forces that applies on the beads, R s . This can be done using the following real transformation

R s = k C s,k R k (4.22)
with the matrix C s,k given by Of course the case k = P/2 is only possible for even values of P and thus is not treated in the case of odd values of P .

C s,k =              1/ √ P
• Finally we apply the quantum thermal bath on the beads using the random forces R s .

Computation of macroscopic properties

In chapter 2 we have seen that several estimators exists to compute macroscopic quantities in PIMD. In particular, the kinetic energy can be computed using different estimators : the primitive estimator (eq. (2.49)), the virial estimator (eq. (2.52)) and the centroid virial estimator (eq. (2.53)) that write

K prim = P s=1 p 2 s 2m - 1 2 mω 2 P P s=1 (x s-1 -x s ) 2 K vir = 1 2P P s=1 x s ∂V ∂x s K Cvir = 1 2β + 1 2P P s=1 (x s -x c ) ∂V ∂x s
Generally, in standard PIMD, the first term of the primitive estimator is rewritten using the equipartition of the energy. But in QTB-PIMD simulations, the dynamics of the beads is not classical so the equipartition theorem no longer holds. Since a part of the quantum fluctuations are included in the momenta we have

p 2 s 2m ≥ 1 2β (4.24)
Thus the general expression of the primitive estimator given in (eq. (2.49)) should be used in QTB-PIMD simulations.

As we have already seen the centroid virial estimator in its standard form (eq. (2.53)) is not directly adapted to QTB-PIMD simulations. To avoid this problem a second formulation (using the f

(1)

P function) has been designed in order to ensure, in particular, that K = U in the harmonic case. In this formulation, the k = 0 normal mode is treated classically. Thus the standard centroid virial estimator cannot be used with the f (0) P formulation of the QTB-PIMD method and should only be used in the f (1) P case. In the purely harmonic case, the normal modes of the ring polymer are uncoupled harmonic oscillators thus the standard centroid virial estimator works perfectly in this case and with the f (1) P formulation. However, when anharmonicity is introduced the normal modes are coupled and thus can exchange energy. In analogy with zero-point energy leakage, an energy flow from the k > 0 modes (with more energy) to the k = 0 mode arise. Thus the centroid kinetic energy can become greater than 1/2β and the centroid virial estimator can become approximate even with the f (1) P function.

If we now go back to the harmonic case, there is another way to ensure that K = U (equation (4.17)). We have seen that the k = 0 normal mode is related to the centroid, in particular 1 2

mω 2 x 2 c = 1 2 mω 2 0 q 2 0
and, since the normal modes are harmonic oscillators in this case, their average potential energy and average kinetic energy are equal thus P case. This modified estimator should also be correct in the f (1) P case and for anharmonic cases. The modified centroid virial estimator reads

K mCvir = P 1 2 m ẋc 2 + 1 2P P s=1 (x s -x c ) ∂V ∂x s (4.27)
Finally, the standard virial estimator 2 (equation (2.52)) is correct in any cases (i.e. f (0)

P or f

(1) P ) and directly leads to K = U in the harmonic case. However, as already explained, it is only valid for bounded system and thus the centroid virial estimator is generally preferred. Now we want to compare the two formulations f (0) P and f

(1) P of the QTB-PIMD method and to choose the adequate kinetic energy estimator between equations (2.49), (2.53) and (4.27). The standard virial estimator is used as a reference here to validate the other estimators. The tests are performed on two cases of the one-dimensional Morse potential already studied in section 2.5 with increasing anharmonicity (1/λ 2 = 0.0015 and 1/λ 2 = 0.024). Figure 4.1 shows the comparison between PIMD and QTB-PIMD. The kinetic energy is computed with the different the estimators. We first notice that, as expected, the QTB-PIMD method clearly allows a significantly faster convergence with the number of beads than standard PIMD.

For both methods, the convergence is slower when anharmonicity is increased. In the weakly anharmonic case (1/λ 2 = 0.0015), the potential energy (Figure 4.1 (a)) remains almost constant with the number of replicas P . Since the anharmonicity is weak, the bare QTB method (for P = 1) already provides a value very close to the exact one. The kinetic energy (Figure 4.1 (b)) depends on the estimator. Note that in the two cases of the Morse potential K prim and K Cvir systematically give the highest and lowest values for the kinetic energy, respectively. As expected, we clearly see that K Cvir fails to give the correct kinetic energy in the more anharmonic case. To conclude, this example shows that both definitions of the f P function results in similar convergence with the number of beads. The primitive, virial and the modified centroid virial estimator can be used in QTB-PIMD simulations. But the standard centroid virial estimator is not adapted to the QTB-PIMD method even with the f (1) P function.

Figure 4.2 shows the influence of the friction coefficient γ on the primitive K prim and the modified centroid virial estimator K mCvir . We first see that K mCvir is clearly less sensitive to γ than K prim . Moreover, K mCvir gives the closest estimation to K vir which serve as a reference here. In particular, when used with the f (1) P function, K mCvir gives exactly the same values as K vir and is particularly insensitive to γ. Thus the best combination here seems to use the modified centroid virial estimator with the f (1) P function. 2 Let us note that, in the general case of a system of volume Ω and composed of N particles in three dimensions, the standard virial estimator can be obtained from the general expression of the pressure P:

PΩ = 2 3 K + 1 3P i,s ri,s • fi,s
In the case of a system for which the pressure is zero, we obtain Kvir = - P . Here the relative deviation from the virial estimator (K -Kvir)/Kvir is plotted. The friction coefficient γ is normalised by the angular frequency ωmin given by the centroid angular frequency in the harmonic approximation (normal mode k = 0) ωmin = ω0 = ω/ √ P . The Morse potential used here is the weakly anharmonic one (1/λ 2 = 0.0015).

Results on simple one-dimensional systems 73

Results on simple one-dimensional systems

Now we first compare the convergence of the QTB-PIMD and PIMD methods on onedimensional systems: a simple harmonic potential and a double well potential.

Harmonic oscillator

As as first example and test case, we consider the same harmonic oscillator as in chapter 2. The system is an hydrogen atom in an harmonic potential well. The temperature T and the angular frequency ω are chosen such that ω/k B T >> 1 more precisely ω/k B T ≈ 13 so the system is mainly in its ground state and T << ω/k B thus quantum effects are expected to be significant. Figure 4.3 shows the kinetic and potential energy convergence with respect to the Trotter number P . Since the bare QTB method is exact in the harmonic case, we see as expected that the QTB-PIMD is already converged for P = 1. When the number of beads increases the QTB-PIMD continue to give the exact results. Note that, in the particular case of the f

Double-well potential

We study here the same double well as in section 2.5. We investigate in particular the position distribution of an hydrogen atom in this double-well. As we already saw in section 2.5, the time-independent Schrödinger equation can be written in a dimensionless form (eq. (2.65)) and we find that the problem is governed by the following dimensionless constant

C = 2 2ma 2 V 0 (4.28)
Ph.D. thesis -F. Brieuc -2016 where V 0 is the height of the barrier and 2a is the distance between the two wells. A numerical resolution of equation (2.65) shows that there is a critical value for the constant C (C 0 = 0.731778) for which the ground state energy is equal to the barrier height V 0 (i.e. ε 0 = 1). Then for C < C 0 the ground state energy is lower than V 0 (i.e. ε 0 < 1) and for C > C 0 the ground state energy is higher than the height of the barrier (i.e. ε 0 > 1). In the following, we study the same three cases as in section 2.5 : a case where C = 1 which is greater than C 0 and two cases for which C < C 0 : C = 0.1 and C = 0.3.

d P = +∞ -∞ (ρ -ρ 0 ) 2 dy +∞ -∞ ρ 2 0 dy (4.29)
which quantifies the difference between the position distribution ρ(y) and the exact one ρ 0 (y) obtained by numerically solving the Schrödinger equation. We first see that, as expected, the QTB-PIMD method converges faster than the standard PIMD method. In the first case (C = 1), the distributions obtained with the QTB-PIMD and PIMD methods converge to the exact one (within an error d P ≈ 2%) at P=16 and P=21, respectively. In this case, the exact distribution shows only one maximum when the bare QTB method gives a distribution with two maxima and clearly fails with an error of d P ≈ 40% . In the other cases, the bare QTB method is already able to give a reasonable approximation to the exact distribution thus the convergence of QTB-PIMD is quite fast: P = 7 for C = 0.3 and P = 4 for C = 0.1 (within an error of d P ≈ 2%). In the latter case, the gain obtained by QTB-PIMD compared to PIMD is not very important (lower than a factor of 2). In constrast, the advantage of the QTB-PIMD is more substantial (around a factor of 3) on the convergence of the total energy (see figure 4.5).

Results on realistic systems

We now compare PIMD and QTB-PIMD on more realistic systems. We focus on a ferroelectric transition in BaTiO 3 and the position distribution of a proton in the proton conductor material BaZrO 3 .

Ferroeletric-paraelectric phase transition in BaTiO 3

These calculations have been carried out by H. Dammak

Barium titanate (BaTiO 3 ) is a prototypical ferroelectric crystal described by a complex energy landscape with multiple wells. It is an ionic crystal with a perovskite structure (ABO 3 ) represented in figure 4.6. Ferroelectricity is the property of certain materials to exhibit a spontaneous electric polarization which can be reversed by the application of an electric field. In BaTiO 3 (BTO), the spontaneous polarization is related to atomic displacements leading to different positions for the barycenters of the positive and negative charges. This results in the creation of an electric dipole moment P in the unit cell. Most ferroelectric materials only exhibit this ferroelectric behaviour on a certain range of temperature. In the particular case of BTO, the material becomes paraelectric (i.e. no spontaneous polarization) when temperatures becomes higher than a critical temperature T c . This temperature is the critical temperature associated with a phase transition between the tetragonal and the cubic phase (see figure 4.6). At temperature higher than T c the material is in the cubic phase and is paraelectric ( P = 0) and for temperature lower than T c the material is in the tetragonal phase and is ferroelectric ( P = 0). BTO exhibits other ferroelectric phases at lower temperatures and undergoes a sequence of phase transitions with increasing temperature: rhombohedral (R) -orthorhombic (O) -tetragonal (T) -cubic (C). Here we only consider the ferroelectric-paraelectric (T-C) phase transition.

Quantum effects in this material have been shown to significantly decrease the transition temperatures by ∼ 30 -50 K [10] and to strongly modify the shape of the pressure-temperature phase diagram [11,12]. The ferroelectric properties of BTO are described here using the effective Hamiltonian of Zhong et al. [13,14] whose parameters have been obtained from first-principles calculations. In this description, the degrees of freedom are the dipole moment (or "local modes") of every unit cell and the components of the strain tensor. This approach has been able to yield an excellent description of the sequence of structural phase transitions that BTO undergoes with temperature [13,14]. The critical temperature (T c ) for the T-C transition (i.e the ferroelectric-paraelectric transition) is equal to 300 K and 260 K when computed using standard MD or PIMD, respectively [15]. Figure 4.7 presents the reduced polarization calculated by PIMD and QTB-PIMD as a function of the temperature around the T-C phase transition. We see that when the temperature becomes higher than a certain temperature the polarization goes to zero which indicates the transition from the tetragonal-ferroelectric phase to the cubic-paraelectric one. In this work, we investigate the convergence of the phase-transition temperature (T c ), as a function of the number of beads, by performing QTB-PIMD and PIMD simulation in the isothermal-isobaric ensemble (using a Langevin thermostat and a Langevin barostat [15]) at zero pressure.

Clearly the QTB-PIMD method provides a faster convergence in this case since the critical temperature T c = 257K obtained with only two beads is very close to the temperature (1) P functions referenced as "(0)" and "(1)" here.

Vertical dashed lines indicates the critical temperatures obtained by QTB-PIMD: 202K(P=1), 257K(P=2) and 259K(P=3). The inset provides the convergence of the polarization in the ferroelectric phase with respect to the number of beads for the three methods.

obtained using PIMD with a Trotter number P = 16 (259K). In addition, the inset in figure 4.7 shows that the convergence of the polarization with respect to the number of beads in the ferroelectric phase is faster when using the f (1) P function within the QTB-PIMD method. We see that the bare QTB method (i.e. QTB-PIMD with P = 1) strongly underestimates the critical temperature by ∼ 55 K. When using QTB-PIMD (even with only P = 2) the failures of the QTB method are fixed. In chapter 3 we have seen that the failure of the QTB method on this system are related to the zero-point energy leakage (ZPEL) problem. Thus the effect of the ZPEL are suppressed by the QTB-PIMD combination.

Proton position distribution in BaZrO 3

We now focus on the proton disorder in another perovskite material. More precisely, we investigate the position distribution of a proton (H + ) in the cubic phase of barium zirconate (BaZrO 3 ). Indeed perovskite-type materials (ABO 3 ) have been shown to exhibit high protonic conductivity [16]. Consequently, these materials are interesting for various technological applications in fuel cells, electrolysers, sensors, gas separators... [17] Among these oxides, doped BaZrO 3 (BZO) exhibits one of the highest conductivities combined with a good chemical stability [16,17] which makes it, in particular, a good potential candidate as electrolyte material for proton conducting fuel cells.

Proton conduction in these perovskite materials has been widely studied and it has been shown that the long-range migration is a combination of transfer and reorientation processes [16]. During the transfer step, the hydrogen atom jumps between the two neighbouring O atoms (Grotthuss mechanism) while, during the reorientation stage, the proton performs a rotation around the nearest O atom. There are two possible rotations: one around the Zr-O-Zr axis and one around an axis orthogonal to the Zr-O-Zr axis, with the H atom remaining in the same O-Zr-O plane. The interatomic forces are obtained using a reactive force field [18] that reproduces the ab initio computed activation energies for the transfer and reorientation processes in BZO [18,19] (see chapter 5). Reactivity for the proton is introduced using the empirical valence bond model (see chapter 5 for more details). No dopant is introduced here so the simulations illustrate the situation far from the dopant atoms. To ensure electrical neutrality in the calculations, a uniform background charge is added to compensate for the proton charge.

The simulations are performed on a 3×3×3 simulation box containing 136 atoms. The threedimensional position distribution of an hydrogen atom in BZO at T = 300 K is computed using PIMD and QTB-PIMD. The distributions in the O-Zr-O plane for increasing values of P are displayed in figure 4.8. The classical distribution obtained by standard MD is also displayed on the top part of the figure as an illustration. Quantum effects are important here since the classical and the quantum distributions are clearly different. In the classical case, the distribution exhibits two peaks thus there is two equilibrium positions for the proton. In contrast, in the quantum case, there is only one broad peak in the distribution thus the proton is able to freely rotate around the nearest oxygen atoms and there is only one equilibrium position for the proton. In addition, proton transfer is rarely observed at T = 300 K. Now we see on figure 4.8 that for QTB-PIMD the two peaks disappear for P = 3 and the distribution is almost converged for P = 4. In contrast, standard PIMD requires a value of at least P = 8 to converge in this case. In order to perform a more accurate comparison of the two methods we compute the divergence factor d P (a generalisation of eq. (4.29) to 3D distributions): the results are presented on figure 4.9. The PIMD calculations are converged within an error of d P = 8% for P = 6, while the QTB-PIMD method requires a smaller number of beads, P = 4. We conclude, as expected, that QTB-PIMD is more efficient than standard PIMD, since fewer replicas are needed for convergence.

Conclusion

We have combined the PIMD method with the QTB approach: the QTB is used as a thermostat for standard PIMD in order to improve the PIMD convergence. Compared to standard PIMD, this combination needs less replicas to converge since a part of the quantum effects is included directly in the dynamics of the beads through the QTB. The gain generally is a factor of 2 or 3 but strongly depends on the system and the physical quantities under study. This combination can also be seen in another way where a small number of replicas in QTB-PIMD is able to correct the failures of the QTB-MD (essentially the zero-point energy leakage). The main advantage of the QTB-PIMD method is its ability to give exact results using less replicas than standard PIMD. Unfortunately, as in PIMD, the major drawback of the new method is that the time-dependent correlation functions are not directly accessible.

The combination with PIMD requires the modification of the power spectral density of the random force within the QTB. This spectral density is proportional to a reduced function (f P ), which can be defined in two ways. In one way (f (0) P ), random forces are applied in the same way to all of the normal modes of the polymer. In the other way (f (1) P ), the centroid mode (k = 0) is addressed with a (classical) Langevin thermostat. In the cases studied here, the f (1) P functions seems to give better convergence of the macroscopic quantities with the number of beads. Considering that some of the quantum fluctuations are included in the momenta through the QTB contribution, a modified centroid-virial estimator of the kinetic energy is proposed.

For future works it would be interesting to investigate the possibility of combining in the same way the QTB with methods such as ring polymer molecular dynamics (RPMD) or centroid molecular dynamics (CMD) to compute time-dependent quantities.

Complements of chapter 4 Derivation of equation (4.6)

As we have seen in this chapter PIMD can be performed on the direct coordinates x s or on the normal modes of the ring polymer q k . During the presentation of the combination QTB-PIMD we have used the following relation (eq. (4.6))

x 2 s = 1 P P -1 k=0 q 2 k
which relates the fluctuations of the beads to the fluctuations of the normal modes. Let us derive this formula. Since direct coordinates and normal modes are formally related using a discrete Fourier transform we have

x 2 s = 1 √ P P -1 k=0
q k e 2iπks/P 2 (4.30)

x 2 s = 1 P q 0 + P -1 k=1
q k e 2iπks/P q 0 + P -1

k ′ =1 q k ′ e 2iπk ′ s/P (4.31)
where we have isolated the first k = 0 term in each sum. Let us now make a small change of variables in the second sum so that k ′′ = Pk ′ . We obtain

x 2 s = 1 P q 0 + P -1 k=1
q k e 2iπks/P q 0 + P -1

k ′′ =1
q P -k ′′ e 2iπ(P -k ′′ )s/P (4.32)

x 2 s = 1 P q 0 + P -1 k=1
q k e 2iπks/P q 0 + P -1

k ′′ =1
q P -k ′′ e -2iπk ′′ s/P (4.33)

where we have used the fact that e 2iπs = 1 ∀s. Now to impose that x s are real we need to impose that q P -k ′′ = q * k ′′ . And we also need to impose q k to be real thus q P -k ′′ = q k ′′ . A real transformation between normal modes and direct coordinates is given in equation (4.23). We then obtain

x 2 s = 1 P P -1 k=0
q k e 2iπks/P P -1

k ′′ =0 q k ′′ e -2iπk ′′ s/P (4.34)
where we have reintroduced the k = 0 terms in the sums and thus we have

x 2 s = 1 P P -1 k=0 P -1 k ′′ =0 q k q k ′′ e 2iπ(k-k ′′ )s/P (4.35)
Finally the normal modes are uncoupled harmonic oscillators (in the harmonic case only) thus q k q k ′′ = q k q k ′′ . Since q k = q k ′′ = 0 we finally obtain

x 2 s = 1 P P -1 k=0 q 2 k (4.36)
which is equation (4.6).

Ph.D. thesis -F. Brieuc -2016

Self-consistent resolution of equations (4.11) and (4.20)

In order to generate the random forces for the QTB-PIMD simulations one needs the function f (0) P (ω) which is solution of (4.11) or the function f

(1) P (ω) which is solution of equation (4.20). In this section we present in details the self-consistent method used to numerically solved these two equations. The self-consistent technique has been proposed by Ceriotti and coworkers [8].

Solving equation(4.11)

We want to solve

P -1 k=0 u 2 u 2 k f (0) P (u k ) = h(u)
We can first isolate the first term and using the fact that u 0 = u/ √ P we obtain

f (0) P u √ P = 1 P h(u) - P -1 k=1 u 2 u 2 k f (0) P (u k ) (4.37)
Before solving, let us rewrite the equation using the function

F P (u) = f (0) P u √ P (4.38)
so that equation (4.37) becomes

F P (u) = 1 P h(u) - P -1 k=1 u 2 u 2 k F P u k √ P (4.39)
we drop the superscript for simplicity in the notation. We first need to choose an initial (guess) solution with good asymptotic behavior :

F (0) P (u) = 1 P h(u/P ) (4.40)
which is the exact solution in the case P = 1. From this initial solution, the equation (4.39) is iteratively solved as

F (i+1) P (u) = α P h(u) - P -1 k=1 u 2 u 2 k F (i) P u k √ P + (1 -α)F (i) P (u) (4.41)
where α is a "weighting parameter" whose value giving the best convergence has been empirically found to be close to 1/P . This method converges to the exact solution F P (u) with a good accuracy after around 30 iterations. The function f In practice, the random forces are generated in the angular frequency space on a range [ω min , ω max ] with ω max being the cut-off frequency for the random force generation (see section 2.3) so ω max = ω cut and ω min is related to the time step δt and the number of steps used in the simulation N s so that ω min = 2π/(N s δt). This angular frequency range then defines the range [u min , u max ] = [β ω min /2, β ω max /2] on which the function F P (u) needs to be determined. Equation (4.41) needs the function F P to be evaluated at u k √ P which reaches a maximum value of u 2 max + P 2 , greater than u max . To overcome this problem, the values of F P for u > u max are linearly extrapolated from the last 20% of the u range. The functions f (0) P obtained for different values of P are presented in figure 4.10 -left panel. We see that the functions approaches 1/P for u → 0 and u/P 3/2 for u → ∞.

Solving equation(4.20)

We now want to solve

P -1 k=1 u 2 u 2 k f (0) P (u k ) = h(u) -1
which can be solved in a very similar way as equation (4.11). This time we can isolate the k = 1 term from the sum and using the fact that u 1 = u 2 P + P sin 2 π P we obtain

u 2 u 2 1 f (1) P (u 1 ) = h(u) -1 - P -1 k=2 u 2 u 2 k f (1) P (u k ) (4.43)
Before solving, let's rewrite the equation using the function

F P (u) = f (1) P (u 1 ) (4.44)
so that equation (4.43) becomes

F P (u) = u 2 1 u 2 h(u) -1 - P -1 k=2 u 2 u 2 k F P P u 2 k -P 2 sin 2 (π/P ) (4.45)
again we drop the superscript for simplicity in the notation. We choose an initial (guess) solution with good asymptotic behavior :

F (0) P (u) = 1 P -1 h u P - 1
P (4.46)
Then equation 4.39 is iteratively solved from this initial solution using

F (i+1) P (u) = α u 2 1 u 2 h(u) -1 - P -1 k=2 u 2 u 2 k F (i) P P u 2 k -P 2 sin 2 (π/P ) + (1-α)F (i) P (u) (4.47)
where α is also giving the best convergence for values close to 1/P . The method converges to the exact solution F P (u) with a good accuracy after around 30 iterations. The function

f (1)
P (u) is obtained from :

f (1)
P (u) = F P P u 2 -P 2 sin 2 (π/P ) (4.48)

Ph.D. thesis -F. Brieuc -2016

Complements of chapter 4 83

Then we generate the stochastic forces from the power spectral density

I R (ω, T ) = 2mγκ(ω, T ) with κ(ω, T ) = P β f (1) P β ω 2 = P β f (1) P (u)
In the same way, the random forces are generated on a range [ω min , ω max ] which defines the range [u min , u max ] = [β ω min /2, β ω max /2] on which the function F P (u) needs to be determined. Equation (4.47) needs the function F P to be evaluated at P u 2 k -P 2 sin 2 (π/P ) which, for the highest u value u max , gives P u 2 max + P 2 1 -sin 2 (π/P ) which is greater or equal than u max . To overcome this problem, the values of F P for u > u max are linearly extrapolated from the last 20% of the u range in the same way as for f (0) P . Another practical problem also arise for the low values of u, from eq.(4.48) one see that that f Quantum effects on the proton conduction in BaZrO 3

Introduction

Proton1 conduction in solids is a key process for various applications. In particular, perovskite type oxides have been shown to exhibit high protonic conductivity [1], and consequently are promising electrolyte materials for proton conducting fuel cells (PCFC). Among the large number of proton conducting oxides, doped barium zirconate (BaZrO 3 ) exhibits one of the highest conductivity coupled with a good chemical stability [2] which makes it one of the most promising candidate for application as electrolyte material for PCFC [1]. Moreover, barium zirconate (BZO) is a perfectly cubic perovskite on a wide range of temperature [3,4]. This highly symmetric structure greatly reduces the number of stable sites for the proton and thus simplify the proton diffusion mechanisms in the material. For all these reasons, BZO is a very good model system to study the quantum effects on the proton conduction in perovskite type oxides.

Since hydrogen is the lightest element of the periodic table, it is expected to exhibit a very strong quantum behaviour compared to the other atoms. In this chapter, we investigate the impact of quantum effects on the protonic diffusion in BZO. As discussed in the general introduction, the importance of quantum effects on the diffusion of hydrogen in solids has been demonstrated in particular for hydrogen in metals [5,6]. It has been shown, both experimentally and numerically, that diffusion in these systems exhibits two different regimes: a classical or semiclassical regime for temperatures higher than a crossover temperature T c where the diffusion is well described by the classical Arrhenius law and a quantum regime for temperatures lower than T c where tunelling becomes important leading to a deviation from the Arrrhenius law [5,6]. A lot of numerical studies have already investigated the diffusion of hydrogen in perovskite oxides in general (see [1,2] for a review) and in barium zirconate in particular [7][8][9][10][11][12][13][14][15][16][17][18]. However, few studies have taken into account the quantum nature of the proton [17][18][19][20]. Among them, Zhang and coworkers [18] have studied proton transport in barium zirconate using ab initio PIMD. Their results show that quantum effects becomes important on the diffusion coefficient for temperatures lower than a crossover temperature T c ≈ 600 K.

In this chapter, we study the diffusion of an interstitial proton (H + ) in barium zirconate (BZO) and investigate the importance of nuclear quantum effects. We are particularly interested to see if the QTB method is able to give valuable results on this type of system and if the obtained results are comparable with those from Zhang and coworkers [18]. First, Section Chapter 5. NQE on the proton conduction in BZO 5.2 gives a description of the elementary processes responsible for the long range diffusion of the proton in BZO. Then, in section 5.3, we give the simulation details and in particular describe the force field used to compute the interatomic interactions. Section 5.4 presents some preliminary results on the structure of BZO and on the diffusion of the proton. Finally, section 5.5 presents the results on the diffusion coefficient as obtained using QTB-MD and the impact of quantum effects on free energy barriers associated to the transfer step obtained using PIMD both for hydrogen and deuterium.

Elementary processes for proton diffusion

As explained in the general introduction (chapter 1), long range migration of protons in perovskites arises from a combination of two types of elementary processes: transfer and reorientation (or rotation) [1,7] (see figure 5.1). During the transfer step (T), the proton jumps between two neighbouring oxygen atoms (Grotthuss mechanism). In contrast, during the reorientation steps, the proton remains covalently bonded to the nearest oxygen atom and performs a rotation around it. In the first rotation (R 1 ), the proton rotates by 90 Let us note that the distribution has been symmetrized using the symmetry of the unit cell in order to obtain the eight equivalent positions since proton transfer is rare at this temperature.

Simulations Details

Force field description

The interatomic interactions are described using an ab initio based force field proposed by Raiteri et al. [13]. This force field reproduces the potential energy barriers associated with the transfer and reorientation mechanisms computed using density functional theory [13,16].

In the following, O 1 represents the nearest oxygen atom to the proton, O 2 the second nearest one (see figure 5.3) and O represents all the other oxygen atoms.

The short range interatomic interactions are described using Born-Mayer potentials of the form

A exp - r ρ (5.1)
with the following parameters: The long range interactions are given by coulombic interactions using nomimal charges: q Ba = +2e, q Zr = +4e and q O = -2e. With e being the elementary electric charge. The Van der Waals interactions between oxygens are neglected in this force field [13].

Buckingham A (eV) ρ (Å)
r 1 r 2 O 1 H 1 2 r 1 r 2 H r 1 r 2 O 2 O 1 O 2

The empirical valence bond model

We need to include chemical reactivity for the hydrogen (i.e. the possibility of bond breaking) to describe the transfer step. This is done through the empirical valence bond (EVB) [21,22] model which has been previously used to simulate proton transfer in various systems [13,[21][22][23][24].

In this model, the potential energy is constructed as a combination of potential energies (V 1 and V 2 ) associated with the proton being in an initial state (reactant state) |1 and in an final state (product state) |2 . This is done in a way reminiscent of the superposition of states in quantum mechanics so that the state vector of the proton would be

|ψ = c 1 |1 + c 2 |2 (5.2)
and the corresponding Hamiltonian could be written in matrix form in the {|1 , |2 } basis set as

Ĥ = V 1 V 12 V 12 V 2 (5.3) with V 1 = 1| Ĥ |1 , V 2 = 2| Ĥ |2 and V 12 = 1| Ĥ |2 = 2| Ĥ |1 .
In the initial state |1 , the proton is covalently bonded to the nearest oxygen atom O 1 while, in the final state |2 , the proton is covalently bonded to the next nearest oxygen atom O 2 (see figure 5.3). V 1 is the total potential energy of the system computed with H covalently bonded to O 1 and V 2 is the total potential energy of the system computed with H covalently bonded to O 2 . The two states are coupled through the coupling term V 12 and the total potential energy of the system V is given by the lowest eigenvalues of the matrix Ĥ:

V = (V 1 + V 2 ) -(V 1 -V 2 ) 2 + 4V 2 12 2 (5.4)
The coupling term V 12 is a function of the reaction coordinate associated with the transfer step ξ( r 1 , r 2 ) = r 2r 1 with r 1 and r 2 being the distance O 1 -H and O 2 -H, respectively (see figure 5.3):

V 12 = λ exp(-α ξ 2 ) (5.5)
with the parameters λ = 0.7998 eV and α = 16 Å -2 whose values have been fitted to reproduce the energy barrier for the transfer step computed using density functional theory [13].

Ph.D. thesis -F. Brieuc -2016 The red curve is the potential energy for the proton in the |1 state and the green is the potential energy for the proton in the |2 state. Finally the blue curve is the coupling between the two states. We see that the EVB construction gives a total energy with a symmetric double well shape (dashed black line).

Simulations Details

93

-2 -1 0 1 2 potential energy ξ (Å) V 1 (ξ) V 2 (ξ) V 12 (ξ) V(ξ)
Finally, the interactions within the O-H-O complex are described differently than the others. Partial charges are used for the hydroxide ion q O 1 = -1.308698e and q H = 0.308698e.

The O 1 -H interaction (covalent bond) consists only of an harmonic interaction with a force constant k = 46.016 eV•Å -2 and an equilibrium distance r 0 = 0.985357 Å (no coulombic or other interactions). The H-O 2 interaction (hydrogen bond) consists of a coulombic interaction using the charges q H = 0.308698e, q O 2 = -2e and a repulsive interaction of the form B/r 2 12 with B = 6.391 eV•Å -12 .

To summarize, the total potential energy of the system for the proton in state |1 is first computed:

V 1 = k 2 (r 1 -r 0 ) 2 + B r 2 12 + 1 4πε 0 q H q O 2 r 2 + U (5.6)
with q O 1 = -1.308698e and q O 2 = -2e. The last term U represents all the other interactions (i.e. coulombic and Buckingham type interactions).

Then the total potential energy of the system with the proton in state |2 is computed:

V 2 = k 2 (r 2 -r 0 ) 2 + B r 1 12 + 1 4πε 0 q H q O 1 r 1 + U ′ (5.7)
The role of O 1 and O 2 is interchanged and in particular the electric charges used in the coulombic interactions now are q O 1 = -2e, q O 2 = -1.308698e. Finally, the total potential energy V is obtained according to equation (5.4). Let us note that the role of O 1 and O 2 is interchanged for all the interactions thus even for the coulombic and the Buckingham ones so that the interactions between Ba-O and O 2 interchanged. Figure 5.4 schematically represent the potential energy associated with the interaction within the O-H-O complex. One see that, as expected, the EVB model gives a potential energy surface with a double well shape along the reaction coordinate for the transfer step ξ.

Choice of the parameters and details of MD simulations

We have performed Langevin MD, QTB-MD and PIMD simulations on a 3× 3× 3 simulation box of cubic BaZrO 3 with one proton (H + ) containing 136 atoms (27 Ba,27 Zr, 81 O and 1 H). Even though doping the structure is necessary to include protons in practice, we did not include any dopant atoms in our calculations. Thus the following results describe the situation in regions far from the dopant. Of course, removing the dopant atom of the study is an important approximation since it is well known that dopant can have a significant impact on the diffusion of the protons [9,12,25]. In particular, the effect of trapping of the proton by the dopant cannot be investigated here. However, the trapping effect is important for high concentration of dopant. At low concentration (lower than ≈ 15%), the effect of the dopant atom on the protonic conduction in Y-doped BaZrO 3 is negligible [1]. Indeed, Raiteri and coworkers have computed the proton diffusion coefficient without dopant and with a small concentration of yttrium and have found very similar values of the diffusion coefficient in both cases [13]. A more complete study should also include dopant atoms. Finally, since no dopant is included, we need to compensate for the proton charge thus a uniform background charge is added to ensure electrical neutrality in the simulation box. Let us also emphasize the fact that we included only one hydrogen atom in our calculations. Thus possible effects coming from interactions between protons are not taken into account in our study. This approximation is reasonable since the proton concentration in these type of compounds is generally rather low (5-20%) [1] so that proton interactions should be negligible at least in first approximation.

Preliminary study and tests

Before carrying out the simulations, we want to check the ability of the force field to correctly describe the system. Even though, Raiteri and coworkers already verify the results obtained during the development of the force field [13], we want in particular to test that our implementation is correct. We first check that some of the basic properties of pure barium zirconate are correctly reproduced, then we include the proton and test that we obtain the correct dynamics. Finally, before including quantum effects we study in details the classical dynamics of the proton and its diffusion.

Structure of barium zirconate

Since we want to compare our results with the results obtained by Raiteri et al. [13], the simulations in this section are carried out with the same size for the simulation box (6×6×6). We focus on the structure of pure BZO and plot the evolution of the relaxed energy (at T = 0K) as a function of the unit cell volume (see figure 5.5). We can then fit the curve with the Murnaghan equation of state [26]:

E(Ω) = E 0 + B 0 Ω B ′ 0 (Ω 0 /Ω) B ′ 0 B ′ 0 -1 + 1 - B 0 Ω 0 B ′ 0 -1 (5.8)
where Ω is the unit cell volume, B 0 and B ′ 0 are the bulk modulus and its pressure derivative at the equilibrium volume Ω 0 and E 0 is the relaxed energy at the equilibrium volume. We obtain, from the fitting, an estimation for the equilibrium cell parameter a 0 = Ω 1/3 0 and for the bulk modulus B 0 . The results are given in table 5.1 and compared to experimental results and the results obtained by Raiteri and coworkers.

We see that the force field reproduces the experimental cell parameter a 0 and that we obtain a value very close to that of Raiteri and coworkers. For the bulk modulus B 0 , we obtain the same value as Raiteri et al. and we see that the force field tends to overestimate the bulk modulus compared to experimental values. This overestimation will not have any impact in the following since all the simulations are carried at zero pressure.

Proton diffusion

We now include a proton in the structure and study in details its dynamics using standard (Langevin) MD. The simulations have been performed with a time step of δt = 0.3 fs and average values have been obtained using 8 independent trajectories of 25×10 6 steps computed after 1×10 6 steps of equilibration. Let us first note that the inclusion of the proton results in a slight increase of the equilibrium cell parameter (∼ 0.01%) whose value is now a 0 ≈ 4.1911Å. Before computing the proton diffusion coefficient we study the elementary processes of the proton diffusion. One can see the two types of elementary processes, rotation and transfer, described in section 5.2. The arrows on the figure indicate a possible diffusion pathway for the proton and we see that a combination of rotation and transfer is necessary for long range diffusion. In order to study these processes more precisely we compute the static potential energy barriers associated with these two processes. Figure 5.7 shows the static potential energy barriers calculated along the reaction coordinates associated with the rotation (R 1 ) and transfer steps. These barriers are obtained by relaxing the system with the reaction coordinate constrained to a certain value. For the transfer step, we choose a standard form for the reaction coordinate ξ( r 1 , r 2 ) = r 2r 1 with r 1 being the O 1 -H distance and r 2 the O 2 -H distance. For the rotation step, we define the reaction coordinate θ as the dihedral angle between the initial O 1 -Zr-O 2 plane and the O 1 -Zr-H plane. The constraints have been included using the SHAKE algorithm [28,29] (see the complements of the chapter for more details). The energy barrier associated with the transfer is higher than the one associated with rotation, indicating that the transfer step might be the limiting step for the diffusion of the proton as it is generally accepted [2]. These energy barrier values (0.38 eV for the transfer and 0.29 eV for the rotation) are higher than the typical values reported in the litterature [9,17,18]. Sundell and coworkers [17] for example report values of around 0.2 eV for the transfer step and around 0.18 eV for the rotation step. However, higher values have also been reported for example Gomez and coworkers [12] report transfer barriers of around 0.4-0.5 eV far from the dopant in Y-doped barium zirconate2 . In their paper, Raiteri and coworkers surprisingly report a static barrier for the transfer of around 0.038 eV (We think that it is a typing error).

Elementary processes of diffusion

Diffusion coefficient

As explained in more details in the complements of this chapter (section 5.7), the diffusion coefficient can be computed in different ways. One standard method is to compute the mean square displacement (MSD), which is related to the diffusion coefficient D through the Einstein relation: lim

t→∞ ( r(t + t 0 ) -r(t 0 )) 2 = 6Dt (5.9)
which, in this form, is only valid for a three dimensional isotropic diffusion. This is the case here since we study barium zirconate in the cubic phase 3 . The quantity ( r(t)r(tt 0 )) 2 is the mean square displacement. In practice, we compute the MSD of the proton and then evaluate the diffusion coefficient from a linear fit in the long time limit. The computed diffusion coefficient of the proton computed from the MSD is compared to the results of Raiteri and coworkers in figure 5.8. Since it is known that the diffusion coefficient can sometimes exhibit a strong dependence on the size of the simulation box [31] we have computed here the diffusion coefficient for two different sizes 3 × 3 × 3 and 6 × 6 × 6. We see that the results are very similar so reducing the size of the box to 3 × 3 × 3 does not have a significant impact on the diffusion coefficient, and thus simulations will now always be done on the 3 × 3 × 3 simulation box. Moreover, the results are very close to those obtained by Raiteri et al. thus validating our implementation of the force field. We obtain, in particular, the same activation energy E a ≈ 0.43 eV but there is a small discrepancy on the prefactor D 0 which seems to be underestimated in our calculations. Raiteri and coworkers have evaluated the diffusion coefficient in a different way based on the jump frequency of the proton. Indeed, since the diffusion of the proton is relatively simple here, the diffusion coefficient can be evaluated by computing the frequencies associated with the rotation (Γ rot ) and the transfer (Γ tr ) of the proton [13]. Figure 5.9 shows the evolution of these two frequencies Γ rot and Γ tr with temperature.

We see that the rotation frequency is significantly higher than the transfer frequency hence the transfer step is indeed the rate limiting step as it is generally accepted in the litterature. This is consistent with the higher static energy barrier obtained for the transfer step. Moreover, from these jump frequencies, we find activation energies of around 0.

E a = 0.425 ± 0.01 eV D 0 = 6.5 10 -4 ± 0.6 10 -4 cm 2 /s E a = 0.437 ± 0.02 eV D 0 = 4.5 10 -4 ± 0.9 10 -4 cm 2 /s E a = 0.429 ± 0.005 eV D 0 = 3.5 10 -4 ± 0.2 10 -4 cm 2 /s Figure 5.8: Diffusion coefficient for the proton in barium zirconate obtained using Langevin MD and computed from the MSD (eq. (5.9)) for a simulation box of size 3 × 3 × 3 and 6 × 6 × 6. The results of Raiteri et al. is also plotted for comparison. We see that the three calculations gives almost the same activation energy Ea ≈ 0.43 eV but there is a small difference on the prefactor which is explained in the following. .9: Jump frequencies for the rotation and transfer steps as obtained from Langevin MD. The two type of frequencies exhibit a linear behavior (Arrhenius behavior) which is expected since they are related to thermally activated processes. Frequencies for rotations are at least one order of magnitude higher than the transfer frequencies on this temperature range thus the transfer step is the limiting step.

0.43 eV for the rotation and the transfer steps respectively, that are comparable with the static energy barriers (T = 0 K) of 0.29 eV for rotation and 0.38 eV for transfer. Moreover, since Γ rot >> Γ tr , the proton can be considered to almost freely rotate around the oxygen atom. Thus the diffusion coefficient is directly related to Γ tr [8,13,32] via the following relation:

D = 1 6 Γ tr d 2 (5.10)
where d is the jump length associated with the transfer. A derivation of this formula is given in the complements of this chapter. Let us note that this relation only holds if the proton diffusion can be well decribed by a three dimensional uncorrelated random walk. In their paper, Raiteri and coworker have used this formula to evaluate the diffusion coefficient with a jump length d equal to the cell parameter a (equation ( 14) in ref. [13]). However, the correct jump length associated with the transfer step here is a/ √ 2 and indeed we see in figure 5.10 that the diffusion coefficient computed using the MSD or the transfer frequencies (eq. 5.10) are equivalent when using d = a/ √ 2. This indicates that, the proton diffusion in this material is well described by an uncorrelated random walk. This difference in the definition of the jump length is responsible for the small discrepancy in the prefactors we found (figure 5.8). We see in figure 5.10 that if we adapt the results of Raiteri et al. using the correct definition of the jump length (d = a/ √ 2) we indeed obtain the same results and the difference in prefactors have now disappeared. The obtained activation energy and prefactor are given in table 5.2 along with some other values from experimental and numerical studies for comparison. We see that the diffusion coefficient obtained using the force field compare quite well with experimental and previous numerical results, in particular the value of the activation energy. Table 5.2: Comparison of the obtained activation energy and prefactor for the proton diffusion coefficient with some experimental and numerical values from litterature. Experimental values have been obtained for Y-doped BaZrO3 with a small dopant concentration (10%) and are taken from [1]. Numerical values have been obtained using the ReaxFF force field and are taken from [10].

Quantum effects on proton diffusion

We now switch to the question of the impact of nuclear quantum effects on the diffusion of the proton. In a first part, we present and analyse the results on the diffusion coefficient as obtained from QTB-MD. Then, in a second part, we study the impact of quantum effects on the free energy barriers associated with the transfer step obtained using PIMD.

Diffusion coefficient from QTB-MD

The QTB-MD simulations have been performed with a time step of δt = 0.3 fs and a friction coefficient of γ = 0.02 THz. The friction coefficient value has been chosen small enough to ensure that the dynamics of the proton is not significantly affected by its value. This is of crucial importance here since we want to study the diffusion of the proton (i.e. we are directly interested in dynamical properties). In practice, we have chosen γ small enough to ensure that the peak in the vibrational density of state associated with vibrations of H is unaffected by the friction. Finally, average values have been obtained using 8 independent trajectories of 25 × 10 6 steps computed after 1 × 10 6 steps of equilibration.

Figure 5.11 shows the diffusion coefficient obtained using Langevin MD and QTB-MD. We see that, as expected, the QTB-MD results deviate from the classical results for temperatures lower than around 800 K. As expected, the impact of the quantum effects is to increase the overall diffusion coefficient.

These results, however, should be regarded with caution: as we have previously seen, a small value of γ as we have used here, could mean that strong zero-point energy leakage (ZPEL) was present in the simulations. In order to verify the presence of ZPEL, we have computed the density of states and the kinetic energy distribution shown in figure 5.12. The presence of ZPEL is clear: the number of high frequency modes in the DOS obtained from QTB-MD is underestimated while the number of low frequency modes is overestimated. The QTB method is not able to counterbalance the leakage and the resulting energy distribution (bottom panel of figure 5.12) is not the expected QTB distribution. We reach an intermediate and almost homogeneous distribution with the small value of γ used here (0.02 THz). As we have seen in chapter 3, increasing γ would reduce this leakage. However, we are interested in dynamical properties here and increasing γ will have an important effect on the computed values of the diffusion coefficient. Thus, increasing γ is not really an option here. It would be useful to define to what extent these results are valuable or not, even with the presence of this strong ZPEL. In order to do this, we compute, in the following section, the free energy barriers associated with the transfer step using both PIMD and standard MD. We want, in particular, to compare the temperature at which quantum effects becomes important on these free energy barriers with the temperature at which the diffusion coefficient starts to deviate from the classical results in QTB-MD simulations (i.e. ≈ 800 K).

Quantum effects on proton diffusion 103

Free energy barriers from PIMD

Like any other dynamical quantities, the diffusion coefficient is not directly accessible in PIMD. However, it can be estimated using transition state theory [18,20,32,33] which provides an approximate expression for the jump rate Γ. As we have seen previously the diffusion coefficient can be related to the jump rate (or jump frequency) through:

D = 1 6 Γd 2 (5.11)
since the diffusion is well described by an uncorrelated random walk in this system. The jump length d = a/ √ 2 and Γ is the total jump frequency thus containing both the rotation frequency (Γ rot ) and the transfer frequency (Γ tr ) through the following relation [8]:

1 Γ = 1 Γ tr + 1 Γ rot (5.12)
Since the reorientation step is very fast compared to the transfer step Γ rot >> Γ tr , the proton can be consider to freely rotate and the diffusion is basically governed by the transfer rate Γ tr leading to equation (5.10). Thus the diffusion coefficient can finally be written as :

D = a 2 12 Γ tr (5.13)
Classical transition state theory allows one to express the jump rate Γ tr as [18,32]:

Γ tr = Γ 0 exp(-β∆F ‡ ) (5.14)
With ∆F ‡ = F (ξ ‡ ) -F (ξ 0 ) the free energy barrier height associated with the reaction: in our case proton transfer. ∆F (ξ) is the free energy profile along the reaction path defined by the reaction coordinate ξ({ r}) (schematically represented in figure 5.13). The saddle point (barrier top) is located at ξ = ξ ‡ while ξ 0 corresponds to the stable site (bottom of the well). The prefactor Γ 0 is generally considered to be independent or weakly dependent on temperature and is often related to an attempt frequency which can be defined as the vibration frequency in the bottom of the well along the reaction coordinate ξ [32].

Here we compute free energy profile ∆F (ξ) along the reaction coordinate for the transfer step ξ( r 1 , r 2 ) = r 2r 1 with r 1 and r 2 the O 1 -H distance and the O 2 -H distance, respectively saddle point (see figure 5.3). This expression for ξ({ r}) is widely used to describe proton transfer [13,18,34,35]. Within this definition the saddle point is located at ξ ‡ = 0. The PIMD simulations are carried out at a constant volume (cell parameter a = 4.975 Å) and at a constant product P × T = 3000 which ensure a good level of convergence of the free energy barriers. The free energy profiles are computed using the thermodynamic integration method which is based on the following relation:

∆F (ξ) = F (ξ) -F (ξ 0 ) = ξ ξ 0 ∂F ∂ξ ξ=ξ ′ dξ ′ (5.15)
In order to obtain the free energy profile, one needs to compute the free energy derivative along the reaction coordinate at fixed values of ξ. This is done using constrained MD simulations (using the SHAKE algorithm [28,29]) with the free energy derivative computed using the blue moon ensemble [36,37] method. More details about free energy computation using this method are available in the complements of this chapter.

In this section, we want to investigate the impact of quantum effects on the free energy barrier for the transfer step. However, transition state theory is a classical theory and thus the expression of Γ tr presented in equation (5.14) is only valid for a classical system. Gillan has proposed a generalisation of transition state theory in the framework of path integrals [38,39] which relates the jump rate Γ to the free energy barrier that the centroid of the ring polymer has to overcome (∆F ‡ c ). In practice, the generalized theory is very similar to the classical transition state theory but with the reaction coordinate ξ replaced by the centroid reaction coordinate ξ c = r 2,cr 1,c 4 . So that the jump rate can be written

Γ tr = Γ c 0 exp(-β∆F ‡ c ) (5.16)
Free energy barrier associated to the centroid ∆F ‡ c can be computed using thermodynamic integration exactly in the same way as in the classical case using equation (5.15) but with the reaction coordinate ξ replaced by the centroid reaction coordinate ξ c . Figure 5.14 present the classical free energy profiles ∆F (ξ) as obtained from constrained standard MD (top panel) and the quantum (centroid) free energy profiles ∆F (ξ c ) as obtained from constrained PIMD (bottom panel) associated with the transfer step. We see first that the free energy barrier, in the classical case, is almost independent of temperature which is consistent with the Arrhenius behavior that we found earlier for the classical diffusion coefficient. Moreover, we find a free energy barrier which is around 0.4 eV consistently with the activation energy of 0.43 eV that we found from the MSD. This further confirmed that the diffusion is mainly governed by the transfer step. As expected, when temperature is lowered, quantum effects become significant and tend to lower the free energy barrier. At 100 K, the classical free energy barrier is ≈ 0.395 eV while the quantum free energy is found to be ≈ 0.274 eV thus quantum effects are responsible for a lowering of 30% of the barrier at this temperature.

Quantum-Classical crossover

Since the dependence of the jump rate with temperature mainly comes from the exponential part of equation (5.14) and (5.16), we have plotted -β∆F ‡ as a function of the inverse temperature in figure 5.15. We see that, for temperatures higher than 300 K, the PIMD and standard MD barrier heights coincides while, for T < 300 K, quantum effects leads to smaller barriers. Thus we expect two different diffusion regimes for the proton: a classical regime at high temperatures, where the diffusion is well described by the Arrhenius law and a quantum regime at low temperatures associated with a deviation from the previous Arrhenius behavior. The transition between this two regimes is around T c = 300 K. This crossover from a classical to a quantum regime has been observed experimentally in the case of hydrogen diffusion in metals [40,41] or on metal surfaces [42,43]. It has also been predicted numerically for various systems including diffusion of hydrogen in metals or on metal surfaces [44][45][46][47][START_REF] Sundell | Quantum motion of hydrogen on Cu(001) using first-principles calculations[END_REF][START_REF] Kimizuka | Effect of temperature on fast hydrogen diffusion in iron: A path-integral quantum dynamics approach[END_REF], hydrogen diffusion in silicon [START_REF] Herrero | Thermally assisted tunneling of hydrogen in silicon: A path-integral Monte Carlo study[END_REF] or proton transfer in molecules [23]. More importantly, this transition between quantum and classical regimes of diffusion is expected to arise at low temperatures for proton diffusion in perovskite oxides [18,20]. Our results are coherent with the results of the path integral study by Zhang et al. which also obtained this quantum-classical crossover on the proton diffusion coefficient in BaZrO 3 . However, Zhang and coworkers found a higher transition temperature T c ≈ 600 K. This is consistent with the fact that, in our calculations the transfer barrier (0.38 eV) is higher than what Zhang and coworkers reports (0.182 eV). From these results, we expect the diffusion coefficient of the proton to deviate from Arrhenius behavior for temperatures lower than around 300 K. It is worth noting that no deviation from the Arrhenius law has been experimentally measured down to 300 K [1] (to the best of our knowledge there are no values of the proton diffusion coefficient in perovskite materials experimentally reported for temperatures below 300 K). In comparison, QTB-MD simulations predict a deviation from Arrhenius behavior for temperatures lower than around 800 -700 K. This discrepancy between QTB-MD and PIMD results is, at least partly, a consequence of zero-point energy leakage. Thus, QTB-MD qualitatively gives the expected behavior i.e. the deviation from the Arrhenius law at low temperature, but significantly overestimates the crossover temperature. In order to get more insights on these two different regimes we plot the probability density of the replicas ρ(ξ) with the centroid fixed at the barrier top (figure 5.16). The spatial extension of the ring polymer gives qualitative information on the diffusion process [18,[START_REF] Mattsson | Quantum Monte Carlo study of surface diffusion[END_REF]. At high temperatures, the distribution of the ring polymer tends to the classical limit ρ(ξ) = δ(ξ-ξ ‡ ). When temperature decreases, the distribution broadens because of quantum fluctuations. We see, in particular, that the shape of the distribution changes with decreasing temperature. This behavior is related to the crossover between the two diffusion regimes. For T > 300 K, the distribution is localised around the saddle point with one clear peak at ξ = ξ ‡ . In this case, the diffusion is mainly governed by overbarrier motions [18,[START_REF] Mattsson | Quantum Monte Carlo study of surface diffusion[END_REF]. In contrast, for temperatures lower than 300 K, the shape of the distribution changes to a bimodal distribution. This is a fingerprint of tunnelling effect [START_REF] Herrero | Thermally assisted tunneling of hydrogen in silicon: A path-integral Monte Carlo study[END_REF] and we enter a regime in which the diffusion is mainly governed by tunnelling events. It is worth noting that tunnelling of hydrogen in perovskite oxides has indeed been experimentally suggested at low temperatures [START_REF] Kuskovsky | Low-temperature dielectric relaxation peaks involving proton tunneling in Ba 1-x Nd x CeO 3[END_REF][START_REF] Cordero | Hydrogen tunneling in the perovskite ionic conductor BaCe 1-x Y x O 3-δ[END_REF].

Isotope effect

Several studies have reported isotope effect on the diffusion coefficient in perovskite oxides [START_REF] Nowick | Isotope effect and proton hopping in hightemperature protonic conductors[END_REF] for hydrogen replaced by deuterium. It is generally found that the diffusion coefficient decreases when hydrogen is replaced by deuterium. The typical shift in the diffusion coefficient can be significant: D H /D D ≈ 1.5 -3 [START_REF] Nowick | Isotope effect and proton hopping in hightemperature protonic conductors[END_REF]. Clearly this effect can come from a change on the free energy barrier ∆F ‡ and/or on the prefactor Γ 0 . The effect on the prefactor is not studied here but we investigate the impact of isotopic substitution on free energy barriers for the transfer step. Figure 5.17 shows the barrier obtained with deuterium or hydrogen as a function of temperature. We see that exchanging H with D tends to increase the barrier for temperature T < 300 K. The crossover between a quantum and a classical regime is also observed in the case of deuterium. The barriers obtained with D are ≈ 0.01 eV higher than the barriers obtained with H. More precisely, the difference ∆F ‡ D -∆F ‡ H varies between 0.008 and 0.014 eV in the temperature range studied here. The isotope effect on the free energy barriers observed here is directly related to quantum effects. In a simple approximation, the isotope effect on the barrier can be related to the difference in zero-point energy between hydrogen and deuterium. If we consider that the frequency at the bottom of the well along the reaction coordinate is close to the OH stretching frequency ν OH , then the free energy difference is approximated by ∆F ‡ D -∆F ‡ H = hν OH (1 -1/ √ 2). Thus, with ν OH ≈ 110 THz which approximately is the value observed in our simulations (see the DOS in figure 5.12), we would expect an isotope effect of around ∆F ‡ D -∆F ‡ H ≈ 0.07 eV. The value we obtained from PIMD simulations is lower (0.01 eV) than this estimation. However this is a very crude approximation: first the frequency at the bottom of the well along ξ is different than the OH frequency. Second, the tunnelling effect, which is mass dependent and, as we have seen, significant at low temperatures, is completely absent of the estimation. Figure 5.18 shows the distribution of the ring polymer with the centroid constrained at the barrier top for both H and D for comparison purposes. First, the system is more localised for deuterium than for hydrogen. This is expected since deuterium is heavier than hydrogen. We also observe the change of the distrbution shape from an unimodal to a bimodal distribution in the case of deuterium. However, in the hydrogen case the change arises around T = 300 K while in the case of deuterium the distribution exhibits two peaks for temperatures lower than 200 K. This is related to the fact that probability for tunnelling events increases when the mass of the tunnelling particle decreases.

Conclusion

In this chapter, we have studied proton conduction in barium zirconate, a potential candidate as electrolyte material for hydrogen fuel cell. In particular, we have focused on the importance of quantum effects on the proton diffusion at intermediate/low temperatures. The interactions are described by a reactive force field based on the empirical valence bond model. In a first part, we have thoroughly tested this force field: we found that the transfer step is rate limiting, as it is generally accepted, and that the computed proton diffusion coefficient compares well with experimental results. In particular, the activation energy is very well reproduced. In a second part, we have focused on the impact of quantum effects. First, we have computed the proton diffusion coefficient using QTB-MD simulations. We found, as expected, that when temperature decreases quantum effects become important and tends to increase the proton diffusion coefficient. More precisely, the QTB-MD results deviates from the standard MD results for temperatures lower than around 800 K. However, since we were interested in the diffusion coefficient here, the simulations have been performed with a very small value of γ in order to ensure that the dynamics of the proton was not perturbed by the thermostat. Thus strong zero-point energy leakage was present and the QTB-MD results have to be regarded with caution. In order to compare the QTB-MD results with PIMD, we have computed the free energy barriers associated with the transfer step using PIMD. These results show that quantum effects become important for temperatures lower than around 300 K and tend to decrease the free energy barrier heights. We thus expect the diffusion to exhibit two diffusion regimes: a classical regime, at high temperature, in which the diffusion is mainly governed by overbarrier motions and the diffusion coefficient follows the Arrhenius law and a quantum regime, at low temperatures, in which tunnelling effect becomes significant leading to a deviation from the Arrhenius behavior. The temperature associated with the crossover between these two regimes is expected to be around 300 K.

The QTB method predicts a similar behavior but significantly overestimates the crossover temperature. From these results, we expect quantum effects to have a significant impact on the proton conduction in BaZrO 3 and in perovskite oxides in general at low temperature. Quantum effects becomes significant for temperatures lower than around 300 K, thus the impact of quantum effects on proton diffusion is probably negligible at the typical working temperature of proton conducting fuel cells (600-900 K). Future works are needed in order to complete this study. First, the free energy barriers associated with the rotation of the proton should be computed. Zhang and coworkers have found that quantum effects have less impact on the rotation free energy barriers than on the transfer ones [18]. This leads to lower barriers for transfer than rotation suggesting that rotation becomes the rate limiting step at low temperatures. It would be interesting to see if we can confirm this result. In a second step, we could try to estimate the prefactor using transition state theory in order to obtain values for the diffusion coefficient.

Complements of chapter 5 Computation of diffusion coefficient

One central quantity in this chapter is the diffusion coefficient D and in particular the proton diffusion coefficient. Here we briefly present some standard methods to compute diffusion coefficients from MD trajectories.

From Mean Square Displacement Probably the most standard way to compute a diffusion coefficient from MD trajectories is to use the well-known Einstein relation (or Einstein-Smoluchowski relation):

lim t→∞ ( r(t + t 0 ) -r(t 0 )) 2 = 6Dt
(5.17)

The quantity ( r(t + t 0 )r(t 0 )) 2 is the mean square displacement. In this form, the relation only holds for an isotropic diffusion i.e. D x = D y = D z = D in 3 dimensions. If the diffusion is not isotropic then the relation needs to be projected on the different directions so that : In practice, the mean square displacement is first computed from the MD trajectories and then the diffusion coefficient is evaluated from a linear fit of the MSD in the long time limit. The computation of the diffusion coefficient from the MSD is completely general and holds for any type of diffusion.

lim t→∞ (x(t + t 0 ) -x(t 0 )) 2 =

From jump frequencies

In general, diffusion in solids can be described as a series of atomic jumps between lattice sites. In this case, the diffusion coefficient can be evaluated from jump frequencies using a random walk model on a lattice. We start from the Einstein relation, writing R(t) = r(t + t 0 )r(t 0 ) the displacement of the atom lim t→∞ R(t) 2 = 6Dt (5.21)

Now if we consider the atom to diffuse via a series of atomic jumps, the displacement is the sum of all the different atomic jumps defined by the vector r i with i = 1, ..., n and n is the number of jumps so that with formally in the t → ∞ limit, the number of jump n → ∞. Now if all the jumps are considered to be independent r i r j = r i r j and if there is no preferred direction r i = r j = 0 thus we obtain

R(t) =
n i=1 r i 2 = 6Dt (5.26) Now if the average jump distance is d nd 2 = 6Dt (5.27)
And we obtain

D = n 6t d 2 (5.28)
Finally n/t is the jump frequency Γ, thus the diffusion coefficient can be written

D = 1 6 Γd 2 (5.29)
which is equation (5.10) that we have used to estimate the diffusion coefficient. In this derivation, we assume a completely uncorrelated random walk i.e. all the jumps are completely independent however there is a lot of cases where this assumption does not hold. One example is the diffusion of atoms via a vacancy mechanism. In this case, the atom diffuses by exchanging lattice site with a neighboring vacancy. Thus this event is only possible when the atom is next to a vacancy. Just after an exchange, the atom is always next to the same vacancy and thus the probability for the reverse exchange to happen is higher than the probability for another exchange. This creates correlations between the jumps and thus r i r j = r i r j in this case. In order to quantify the importance of these correlations, the correlation factor f is defined as

f = D D rand (5.30)
where D rand is the diffusion coefficient that would be obtained if the diffusion was a completely uncorrelated random walk (i.e. D rand = Γd 2 /6). If the diffusion is well described by an uncorrelated random walk then f = 1, but if there are correlations between the jumps then the correlation factor is generally lower than one.

Constrained molecular dynamics

We have seen in this chapter that constrained MD was necessary in order to compute either static potential energy or free energy profile along the reaction coordinate. This has been done in practice using the SHAKE algorithm [28,29]. Here we present this algorithm and its use for the particular case of interest here: fixed transfer and rotation reaction coordinates.

We want to perform MD simulations under the constraint that the reaction coordinate ξ({ r}) remains fixed at a certain value ξ. The constraint can be written g({ r}) = ξ({ r})ξ = 0. What we present here is valid for any type of holonomic constraint 5 . Within the MD framework, the constraint can be kept constant by adding external forces G i that will apply on the different atoms involved in the constraint and that will ensure that the constraint is fulfilled at each time step. The constraint forces G i writes:

G i = -λ ∇ i g({ r}) (5.31)
where λ is the Lagrange multiplier associated with the constraint and whose expression needs to be determined. This is the Lagrange multiplier λ that appears in equation (5.80). In the SHAKE algorithm, the MD trajectories are integrated using the Verlet algorithm:

r i (t + δt) = 2 r i (t) -r i (t -δt) + δt 2 f i (t) m i + δt 2 G i (t) m i (5.32)
Where G i is the force associated to the constraint and f i (t) represents the other forces that applie on the atom i. One can first introduce the unconstrained positions ˆ r i (t + δt) i.e. the positions that the atoms would have at time t + δt if there was no constraint so that

r i (t + δt) = ˆ r i (t + δt) - λ(t)δt 2 m i ∇ i g({ r}) (5.33)
Now if a good guess (λ (1) (t)) for the value of the Lagrange multiplier at time t exists (for example the value at the previous time step) then one can define the new positions according to λ (1) :

r (1) i (t + δt) = ˆ r i (t + δt) - λ (1) (t)δt 2 m i ∇ i g({ r}) (5.34)
and the constrained positions can be written as

r i (t + δt) = r (1) i (t + δt) - δλ (1) (t)δt 2 m i ∇ i g({ r}) (5.35)
where δλ (1) (t) = λ(t)λ (1) (t). Now we want to ensure that the constraint remains fulfilled at time t + δt so g ( r 1 (t + δt), ..., r N (t + δt)) = 0 (5.36) which writes g r

(1) 1

δλ (1) (t)

δt 2 m 1 ∇ 1 g({ r(t)}), ..., r (1)
Nδλ (1) (t)

δt 2 m N ∇ N g({ r(t)}) = 0 (5.37)
Now performing a first order Taylor expansion around δλ (1) (t) ≈ 0 one obtains g r

(1) 1 , ..., r (1) (t + δt)}) ≈ 0 (5.38)

(1) N - N i=1 δλ (1) (t) δt 2 m i ∇ i g({ r(t)}) • ∇ i g({ r
The gradient of the constraint with respect to the centroid positions are:

∇ O 1 ,c g = d O 1 H,c d O 1 H,c (5.48)
∇ O 2 ,c g = - d O 2 H,c d O 2 H,c (5.49)
∇ H,c g = - d O 1 H,c d O 1 H,c + d O 2 H,c d O 2 H,c (5.50)
The forces that have to be applied on the centroid are given by:

G O 1 ,c = -λ ∇ O 1 ,c g (5.51) G O 2 ,c = -λ ∇ O 2 ,c g (5.52) G H,c = -λ ∇ H,c g (5.53)
and the forces that applies on the replicas are obtained using:

G O 1 ,s = -λ ∇ O 1 ,s g (5.54) G O 2 ,s = -λ ∇ O 2 ,s g (5.55) G H,s = -λ ∇ H,s g (5.56)
which writes:

G O 1 ,s = -λ 1 P d O 1 H,c d O 1 H,c = 1 P G O 1 ,c (5.57) G O 2 ,s = λ 1 P d O 2 H,c d O 2 H,c = 1 P G O 2 ,c (5.58) G H,s = λ 1 P d O 1 H,c d O 1 H,c - d O 2 H,c d O 2 H,c = 1 P G H,c (5.59)
Finally, the expression for the Lagrange multiplier is obtained using the same equation (5.39).

Rotation step

with Z = z O 1 -z H and Y = y O 1 -y H .
Then the cosine of the angle is given by (5.63)

cos θ = n • e z = -Y √ Y 2 + Z 2 (5.62)
Finally we constrain the square of the sine of the angle sin 2 θ. The forces that needs to be applied are :

G O 1 = -λ ∇ O 1 g (5.64) G H = -λ ∇ H g (5.65)
which writes:

G O 1 = -λ 2ZY 2 e z -2Z 2 Y e y (Y 2 + Z 2 ) 2 (5.66) G H = -λ -2ZY 2 e z + 2Z 2 Y e y (Y 2 + Z 2 ) 2 = -G O 1 (5.67)
Then the value for the Lagrange multiplier is obtained as usual using equation (5.39).

Free energy computation

We want to compute the free energy profile F (ξ) or more precisely ∆F (ξ) = F (ξ) -F (ξ 0 ) along the reaction coordinate ξ( r 1 , ..., r N ). We first need to define the function F (ξ) which is an "incomplete" free energy:

F (ξ) = -k B T ln Z(ξ) (5.68)
where Z(ξ) is an incomplete partition function defined as

Z(ξ) = 1 h 3N d p 1 ... d p N d r 1 ... d r N δ (ξ( r 1 , ..., r N ) -ξ) e -βE({ r},{ p}) (5.69)
where ξ denote some numerical value of the reaction coordinate defined by the expression ξ( r 1 , ..., r N ). The free energy profile ∆F (ξ) can then be expressed as:

∆F (ξ) = F (ξ) -F (ξ 0 ) = -k B T ln Z(ξ) Z(ξ 0 ) (5.70) ∆F (ξ) = -k B T ln δ(ξ( r 1 , ..., r N ) -ξ) δ(ξ( r 1 , ..., r N ) -ξ 0 ) (5.71) ∆F (ξ) = -k B T ln ρ(ξ) ρ(ξ 0 ) (5.72)
So free energy profile can be obtained from MD simulations by monitoring the value of ξ( r 1 , ..., r N ) during the simulation. Then, one can compute the probability density ρ(ξ) of the reaction coordinate to obtain ∆F (ξ) using equation (5.72). The choice of ξ 0 is arbitrary, but we generally take the value of the reaction coordinate associated with the bottom of the well (see figure 5.13) so that, at the saddle point, ∆F (ξ ‡ ) is the free energy barrier associated with the reaction.

In our case, the reaction is the transfer of the proton defined by the reaction coordinate ξ( r 1 , r 2 ) = r 2r 1 . Proton transfer is a rare event and thus very long simulations would be required in order to ensure that ρ(ξ) is correctly sampled around the saddle point. Thus, the free energy barrier can not be accessed using the probability density of ξ in our case (see figure 5.20) . In order to overcome this problem, free energy profiles are computed using the thermodynamic integration method which is based on the following relation:

∆F (ξ) = F (ξ) -F (ξ 0 ) = ξ ξ 0 ∂F ∂ξ ξ=ξ ′ dξ ′ (5.73)
If we can compute the free energy derivative at different values of ξ then we can numerically perform the integral in equation (5.73) to obtain the free energy profile.

For simple expressions of the reaction coordinate, the free energy derivative can be related to the average force that apply on the reaction coordinate [START_REF] Geneste | Landau free energy of ferroelectric crystals by thermodynamic integration[END_REF] at fixed value of ξ. As an example, if the reaction coordinate simply is the position of one particular atom, let us choose atom 1, and if we consider a one-dimensional problem for simplicity then ξ(x 1 , ..., x N ) = x 1 and the free energy derivative can be expressed as : where V is the potential energy. Thus the free energy is sometimes referred to as the potential of mean force acting on the reaction coordinate. The main idea of thermodynamic integration is thus to compute free energy derivative from equation (5.79). The average value in equation (5.79) can be calculated using constrained MD. Several MD simulations are run with the reaction coordinate constrained at different value of ξ({ r}). Then the free energy derivative is integrated along the reaction path to obtained the free energy profile (equation 5.73).

∂F ∂ξ = -k B T 1 Z(ξ) ∂Z(ξ) ∂ξ (5.74) ∂F ∂ξ = -k B T 1 Z(ξ) ∂ ∂ξ 1 h 3N dp 1 ... dp N dx 1 ... dx N δ (x 1 -ξ) e -βE (
However, for more complicated form of the reaction coordinate, the free energy derivative cannot be easily written in the simple form of equation ( 5.79) anymore. Moreover, the constrained MD algorithm not only imposes that ξ({ r}) = ξ, it also imposes that the time derivative of the reaction coordinate ξ({ r}, { ˙ r}) = 0. This is an artifact of the method that is not present in equation (5.79) [37]: the average value should be computed at fixed value of ξ({ r}) only. The blue moon ensemble method [28,37] provides a general formula that is valid for any form of the reaction coordinate and that takes into account the possible artifact coming from the use of constrained MD simulations. The formula writes:

∂F ∂ξ = Z -1/2 [-λ + k B T G] const. ξ Z -1/2 const.
Z 2 N i=1 N j=1 1 m i m j ∂ξ ∂ r i • ∂ 2 ξ ∂ r i ∂ r j • ∂ξ ∂ r j (5.82)
The free energy barriers for the transfer step presented in this chapter has been computed using this formula. The transfer reaction coordinate is ξ Figure 5.20 -top panel shows the free energy barrier as computed from thermodynamic integration with constrained MD in the blue moon ensemble (blue line and cross) and the free energy barrier as obtained from unconstrained MD simulations using equation (5.72) (red full circle). We see that the method are equivalent and give the same profile however since proton transfer is a rare event at this temperature the free energy profile at the barrier top is not accessible from unconstrained simulations.

( r 1 , r 2 ) = r 2 -r 1 with r 2 = | r 2 | = | r O 2 -r H | and r 1 = | r 1 | = | r O 1 -r H |
Z = 2 m O + 2 m H 1 - r 2 r 2 • r 1 r 1 ( 5
Finally, as we have seen previously in this chapter, in the case of PIMD we are interested in computing the free energy profile along the centroid reaction path i.e. ∆F (ξ c ) with ξ c the value of the centroid reaction coordinate ξ c (r O 1 ,c , r O 2 ,c , r H,c ) and r i,c the position of the centroid of atom i. This is done using the same method of constrained MD in the blue-moon ensemble with the centroid reaction coordinate. Figure 5.20 -bottom panel shows the barrier as obtained from unconstrained PIMD (eq. (5.72)) using both the centroid reaction coordinate probability density ρ(ξ c ) and the replicas reaction coordinate probability density ρ(ξ), and the barrier obtained from constrained PIMD in the blue moon ensemble. We see that, since we have constrained the centroid reaction coordinate ξ c , the constrained PIMD scheme is, as expected, equivalent to the free energy profile along the centroid reaction path i.e. as obtained with the centroid reaction coordinate probability density in the unconstrained PIMD scheme. 5.72) (red full circles), and as obtained from constrained MD and computed using thermodynamic integration -equation (5.73) (blue line and cross) at T = 300 K. Bottom: ξc = r2,c -r1,c as obtained from unconstrained PIMD using the probability density ρ(ξc) -equation (5.72) (red full circles), and as obtained from constrained PIMD and computed using thermodynamic integration -equation (5.73) (blue line and crosses). For comparison purposes, the free energy profile along the reaction coordinate of the replicas ξ as obtained from unconstrained PIMD using the probability density ρ(ξ) is also plotted (green triangles) Chapter 6 Proton diffusion mechanisms in GdBaCo 2 O 5.5

Introduction

In this chapter, we investigate the proton diffusion mechanisms in GdBaCo 2 O 5.5 (GBCO) whose structure is represented in figure 6.1. GBCO is a double perovskite with two different A-site cations (Gd and Ba) that alternate along the c-axis [1]. Thus the cell parameter along the c-axis, is doubled compared to the cubic perovskite structure [1]. This ordering of Gd and Ba results in three different types of planes: GdO, CoO 2 and BaO planes alternating along the c-axis. GBCO naturally contains oxygen vacancies in order to ensure the overall electrical neutrality. These oxygen vacancies are mainly located within the GdO planes [1]. Since oxygen vacancies are present, oxygen atoms can diffuse in the structure, and indeed GBCO is a good oxygen conductor [2]. Finally, GBCO exhibits a metal-insulator transition at T ≈ 360 K [1] and thus is a mixed electron/oxygen conductor at high/intermediate temperatures (800 -600 K). For this reason, GBCO has been used as a cathode material in Solid Oxide Fuel Cells (SOFC)1 and has shown good properties [3,4]. Thus, the interest for double perovskite cobaltite compounds in general, and GBCO in particular, for applications as cathode materials in hydrogen fuell cells has been growing up [2][3][4][5][6][7][8].

Oxygen diffusion mechanisms in double perovskite cobaltite compounds have been numerically studied using MD simulations [9][10][11][12]. These studies have evidenced the bidimensional nature of oxygen diffusion in these compounds which has been comfirmed experimentally by neutron diffraction [12]. It has been shown, in particular, that long range diffusion occurs through oxygen ions jumps between gadolinium and cobalt planes [9,10,12] and that BaO planes act as barriers for oxygen diffusion thus hindering long range migration of oxygen along the c-axis.

Recently, GBCO has been applied as a cathode material for Proton Conducting Fuel Cells (PCFC) [6,8] with excellent properties. This could be related to a potential diffusion of protons inside GBCO [8] indicating that this compound could incorporate water and become a mixed proton/electron conductor. This assumption has been further supported by other experimental studies [7,8] suggesting that water incorporation would be possible in double perovskite cobaltite compounds. However, recent neutron diffraction did not allow the confirmation of the presence of water in humidified NdBaCo 2 O 5+x at high temperatures [13].

Finally, DFT calculations have been performed in order to study hydration of GdBaCo inal structure considered [14]. Further experimental and theoretical studies are thus needed to investigate the possible hydration of double perovskite cobaltite compounds in order to decide whether or not these compounds can be mixed proton/electron conductors.

Mixed proton/electron conductors are interesting because they are expected to exhibit very good properties as electrode materials in PCFC. At the cathode, the reaction in PCFC is

4H + + O 2 + 4e -→ 2H 2 O (6.1)
Most cathode materials for PCFC are mixed oxygen/electron conductors. In this case, the complete electrode reaction can only occur at the triple phase boundary points near the electrolyte surface where electrons, protons and oxygen can meet [8,15,16] as represented in figure 6.2. In contrast, if the cathode material is a mixed proton/electron conductor, the cathode reaction can potentially occur on the whole cathode surface. Thus mixed proton/electron conductors are expected to be more efficient than mixed oxygen/electron conductors as cathode materials for PCFC [8,15,16].

As we just explained, the possibility to incorporate protons in GdBaCo 2 O 5.5 is still under debate. In the following, we assume that at least a small hydration of the compound would be possible and thus we study the proton diffusion mechanisms in GBCO using standard molecular dynamics2 . In the first section, we give the computational details and in particular describe the force field used to model the interatomic interactions. In a second section, we focus on the diffusion mechanisms of the proton in GBCO. Finally we discuss the results and conclude.

Computational details 127

.2: Schematic representation of the cathode and electrolyte interface for protonic conducting fuel cell.

Computational details

The interatomic potentials used in this study are those previously used to study oxygen diffusion in GdBaCo 2 O 5.5 [10]. Metal-metal and metal-oxygen interactions are described using Buckingham potentials of the form:

A exp - r ρ - C r 6 (6.2)
The interactions between the first neighboring oxygen ion O 1 of the proton and the other cations are adapted to take into account the smaller charge on the oxygen ion of the hydroxyl group, following a procedure presented by Sierka et al. [17]. The values of the parameters for the Buckingham potentials used in this study are given in The proton transfer between the two first neighboring oxygen atoms is described by the empirical valence bond (EVB) model presented in chapter 5. The EVB parameters associated to the interactions between the proton and the other ions are taken from the original article of Raiteri et al. [19] for barium zirconate and are given in chapter 5. Here, we assumed that the values of the EVB parameters (in particular λ, α and the values of the electrical charges within the OH group) are transferable and thus can be used to model the O-H interaction within GBCO. Of course, this is a strong assumption since these values have been fitted to reproduce the transfer and rotation barrier for a proton in BaZrO 3 (BZO). This assumption is based on the fact that, since both BZO and GBCO are perovskite oxides, barriers for rotation and transfer should be similar between the two compounds3 . Thus, this force field cannot claim to perfectly reproduce the potential energy surface of the proton in GBCO, nevertheless the orders of magnitude of energy barriers should be realistic. In particular, by computing rotation and transfer frequencies, we find that the rotation step is a fast process compare to the transfer step in GBCO, similarly to what we found in BZO and as it is generally accepted in the literature for perovskite oxides [22].

The size of the simulation box is 6a×6b×3c as compared to the original double perovskite cell, containing 108 Ba, 108 Gd, 216 Co and 594 O, corresponding to the formula GdBaCo 2 O 5.5 .

In the initial configuration, all the oxygen vacancies are located in GdO planes, since oxygen vacancies are experimentally found to be mainly in these planes [1]. One proton is introduced in the simulation box and the electrical neutrality is ensured by adding a uniform background charge. The simulations were performed with a time step of 0.3 fs in the NVT ensemble using a Langevin thermostat and at a pressure close to zero. The temperature was varied between 1200 and 1900 K and average quantities were obtained using equilibrium trajectories of 11 ns to 2 ns, respectively.

Proton diffusion mechanisms

Let us first note that the presence of one proton over 1026 atoms has almost no effect on the volume of the cell. As expected [10], oxygen atoms diffuse in parallel with the (a,b) plane but are blocked in the c direction. Up to 1700 K, BaO planes act as barriers for the diffusion of oxygen atoms explaining their limited displacement along the c direction, and thus the 2D nature of the diffusion. Figure 6.3 shows the mean-square displacement (MSD) of the proton as obtained from simulations at 1700 K. The proton diffuses in parallel to the (a,b) plane and is blocked along the c direction, similarly to oxygen atoms. In order to elucidate the origin of this anisotropic diffusion, one can first study the preferential locations of the proton in the material. Figure 6.4 shows the probability density map of a proton in GBCO as obtained from MD trajectories. The proton diffuses within the same GdO and CoO 2 planes which explains the bidimensional nature of proton diffusion parallel to the (a,b) plane. The proton is never found in the BaO plane: similarly to oxygen diffusion, BaO planes act as barriers for hydrogen diffusion and thus hinder long range migration of the proton along the c-axis. insight on the equilibrium locations of the proton in the material. We first note, that the distance between the two oxygen planes located on either side of the Gd plane is less than c/2 by 13%. Consequently, the distance between the two oxygen planes located on either side of the Ba plane is greater than c/2 by 13% and the proton remains between the Gd and Co planes. We also see that the hydrogen is mainly bonded to oxygen atoms of the GdO planes. More precisely, the proton is located in GdO planes during 93% and 83% of the simulation time at 1300 and 1700 K, respectively (see table 6.2). Moreover, the OH bond mostly points toward a direction that keeps the proton between the gadolinium and the cobalt planes. Comparison is made with previous experimental and theoretical results, which only exist for oxygen atoms. The calculated oxygen diffusion coefficient values are in reasonable agreement with experimental results [2,23]. The diffusion coefficient is found to be more than one order of magnitude higher for hydrogen than for oxygen atoms. The activation energy associated with proton diffusion (0.71 ± 0.03 eV) is found to be lower than the activation energy for oxygen diffusion (1.01 ± 0.03 eV). This difference of activation energy between oxygen and proton migration is typically observed in perovskite ion conducting compounds. For instance, in yttrium-doped barium zirconate, an activation energy of around 1 eV [24] is reported for oxygen ion migration, while it is around Ph.D. thesis -F. Brieuc -2016

Proton diffusion mechanisms 131

0.43 eV [22] for the proton migration. From the values of the diffusion coefficient, one can estimate the conductivity σ H of protons in GBCO using the Nernst-Einstein relation that relates the diffusion coefficient to the conductivity:

σ = q 2 nD kT (6.3)
with σ the conductivity and n the density of charge carrier of charge q. If we consider a small hydration of the compounds, e.g. 5% of the vacancies are filled by water, at the typical working temperature of PCFC i.e. T = 600 • C, we obtain that σ H ≈ 0.5σ O where σ O is the estimated conductivity for oxygen diffusion. Thus, even for a small hydration of GBCO, we expect the conductivity of the proton to be comparable to the conductivity of oxygen atoms in this compound. [9], (b) from oxygen exchange -ref. [2], (c) from conductivity relaxation -ref. [23] Further analysis of the trajectories reveals that the diffusion of the proton in GBCO follows a scheme where it is successively bonded to an oxygen atom of the gadolinium plane and then to an oxygen atom of the cobalt plane. Moreover, the long range migration of the proton is obtained through two different mechanisms (schematically represented in figure 6.7). The first one is the standard proton transfer between the two neighboring oxygen atoms combined with the rotation of H around its first neighbor oxygen atoms (mechanism (1) in figure 6.7). The limiting step is the proton transfer as evidenced by figure 6.8 showing the Arrhenius plot of the corresponding jump frequencies. The activation energies are equal to 0.63 and 0.38 eV for the proton transfer and for the rotation, respectively. In addition, the proton transfer from GdO to CoO 2 planes or reciprocally represent around 93% of the total number of transfer jumps at T = 1500 K, the other ones being within the CoO 2 planes (7%) (see table 6.3). The second mechanism consists in the migration of the OH group in which both oxygen and hydrogen atoms jump simultaneously, mainly along <101> directions of the cubic cell (mechanism (2) in figure 6.7). At 1500 K, the contributions of the two diffusion mechanisms, transfer-rotation and OH migration, are of 70% and 30%, respectively. extracted activation energies are very similar, being equal to 0.64 and 0.59 eV for the transferrotation and OH migration mechanisms, respectively. These activation energies are much smaller than the ones calculated from the diffusion coefficient values obtained using the mean square displacement. This can be explained by a spatial correlation between successive jumps leading to a non-purely random walk of the proton i.e. the diffusion cannot be described by an uncorrelated random walk.

Correlation factor

As we have seen in chapter 5 (see the complement in particular), diffusion in solids can generally be described by a series of atomic jumps on a lattice. In this case, the diffusion coefficient can be related to the atomic jump frequency. If the diffusion is well described by an uncorrelated random walk on a lattice then the diffusion coefficient can be expressed as follows

D = 1 6
Γd 2 (6.4)

with Γ the jump frequency and d the jump length. Bear in mind that this expression is valid for a three-dimensional isotropic diffusion only. We have seen in chapter 5 that the diffusion of hydrogen in BZO was well described by expression (6.4). However, this expression does not hold for every type of diffusion, for example oxygen diffusion can not be described by equation (6.4). Oxygen atoms diffuse via a vacancy mechanism in which an oxygen atom migrates by exchanging lattice site with a neighboring vacancy. In this case, the diffusion cannot be described by an uncorrelated random walk because the different atomic jumps are correlated. Indeed, just after an exchange between an oxygen atom and a vacancy, the oxygen is always next to the same vacancy and thus the probability for the reverse exchange to happen is higher than for another exchange. Thus, there is a correlation between the different atomic jumps in this case and equation (6.4) cannot be used. We have seen that a non negligible part of the proton diffusion in GBCO arises from a diffusion of the whole OH group which diffuses through a vacancy mechanism. Moreover, the presence of oxygen vacancies will sometimes block some jump directions for the proton transfer and the probability for a jump will then depend on the local arrangement of the neighboring sites, and in particular on the number of neighboring vacancies. Thus hydrogen diffusion in GBCO cannot be described by a simple uncorrelated random walk and the diffusion coefficient cannot be estimated from equation (6.4). This is coherent with the fact that the activation energies that we obtain from the jump frequencies for the two mechanisms of hydrogen migration (OH migration and proton transfer/rotation) are smaller than the one calculated from the diffusion coefficient values 134 Chapter 6. Proton diffusion in GBCO obtained using the mean square displacement.

The importance of the correlations between the jumps can be estimated by computed the correlation factor f (see the complements of chapter 5). The correlation factor can be defined as the ratio between the diffusion coefficient obtained from the MSD and the diffusion coefficient obtained from the jump frequencies (in the hypothesis of uncorrelated random walk) [25,26]:

f = D H 1 4 Γ H d 2 H + Γ OH d 2 0H (6.5)
where the factor 1/4 comes from the 2D nature of the diffusion. The frequency Γ H corresponds to the transfer-rotation combination [27]:

1 Γ H = 1 Γ rotation + 1 Γ transfer (6.6)
and Γ OH is associated to the OH migration mechanism where both oxygen and hydrogen atoms move simultaneously. The lengths of jumps d H and d OH are the proton displacements projected in the (a,b) plane. Both lengths are about a/2, since the jumps are oriented mainly along the [101] or [011] directions. In the case of a perfectly uncorrelated random walk f = 1 while, if correlations between jumps are present, f < 1. Our MD values for the correlation factor are equal to 0.28, 0.59, and 0.76 at 1300 K, 1500 K, and 1700 K, respectively. As expected, there exists strong correlations between the jumps which are related to the simultaneous diffusion of oxygen and hydrogen. Indeed, in the previous case (chapter 5) of the proton diffusion in BaZrO 3 , there was no oxygen vacancies thus the oxygen atoms were not migrating, and we obtained a correlation factor close to 1 (0.94).

Conclusion

The diffusion mechanisms of a proton in GdBaCo 2 O 5.5 have been investigated using standard (Langevin) molecular dynamics. We have found that the proton presents a 2D diffusion analogous to that of oxygen in this material. We have evaluated the proton diffusion coefficient to be more than one order of magnitude greater than the oxygen diffusion coefficient. The activation energy associated with proton diffusion that we obtain is equal to 0.71 eV which is smaller than the activation energy for oxygen diffusion. After analysis of the MD trajectories, we found that proton diffuses through two different mechanisms in this material. The main mechanism is the standard proton transfer combined with the rotation of H around its first neighbor oxygen atom. But a non negligible part of the proton diffusion comes from a second mechanism which consists in the migration of the whole OH group i.e. where both oxygen and hydrogen atoms jump simultaneously. These two mechanisms have similar activation energies of around 0.6 eV. We also evidenced strong correlations between the atomic jumps, since we obtained correlation factor values less than unity. These correlations are related to the simultaneous diffusion of hydrogen and oxygen atoms in the material. These results brings some information on the influence of a potential hydration of GBCO: if we consider for instance the filling of 5% of the oxygen vacancies, we would obtain a conductivity for protons at 600 • C that is 0.5 times that of oxygen atoms. In other words, a slight incorporation of protons in GBCO would result in a relatively high proton conductivity, at the working temperature of PCFC.

Part IV

General Conclusion 141

In a first part of this thesis, we have studied in details the QTB method. We have seen that, when dealing with highly anharmonic systems, QTB-MD simulations should be used with caution because energy is transferred between the vibrational modes, due to anharmonic couplings, leading to zero-point energy leakage (ZPEL). ZPEL is the main limitation of QTB-MD simulations and can lead the method to fail in strongly anharmonic cases. We have shown that increasing the friction coefficient in QTB-MD simulations can significantly help limiting ZPEL and even completely remove it in some cases. However, in strongly anharmonic cases, the increase of the friction coefficient is sometimes not sufficient to remove the leakage. The QTB method is thus not suitable for the study of such systems, and conventional methods such as path integral molecular dynamics can be used, in this case.

The QTB method can then be used as a thermostat for PIMD simulations. We have presented a combined QTB-PIMD method which allows to decrease the number of beads required in PIMD simulations, and thus reduces the computational cost of these simulations. The QTB includes a part of the quantum effects in the dynamics of the beads and thus fewer replicas are needed to converge towards the exact quantum results. We have seen that the gain obtained by using QTB-PIMD as compared to standard PIMD depends on the quantity of interest. For the convergence of the total energy for example, the method allows to reduce the number of replicas by at least a factor of 3, using the estimator that we propose for the kinetic energy.

In a second part of this thesis, we have studied proton conduction in perovskite materials. First we have investigated the impact of quantum effects on the diffusion of the proton in BaZrO 3 using PIMD and QTB-MD simulations. We have found that, for temperatures lower than 800 K, the diffusion coefficient obtained from QTB-MD simulations deviates from the Arrhenius law. As expected, quantum effects become significant at low temperatures and tend to increase the diffusion coefficient. However, these results have to be taken with caution since strong ZPEL was present during these simulations. In order to compare these results with PIMD, we have computed the free energy barriers associated with proton transfer using standard MD and PIMD simulations. We found that quantum effects become significant for temperatures lower than 300 K and, as expected, tend to lower the free energy barriers. We expect the impact of quantum effects to become significant on the diffusion of the proton in BaZrO 3 for temperatures lower than 300 K. And thus, we do not expect quantum effects to play a significant role at the typical working temperatures of proton conducting fuel cells (T ≈ 600 -900 K).

Finally, in the last chapter, we have investigated the proton diffusion mechanisms in the double perovskite compound GdBaCo 2 O 5.5 . We have evidenced the bidimensional nature of the proton conduction, similarly to the diffusion of oxygen in this material. Moreover, two different diffusion mechanisms for the proton exist in this material: the main one is the standard proton transfer combined with the rotation of the proton around the nearest oxygen atom, and the second mechanism is the migration of the OH group, where both the proton and the nearest oxygen diffuse simultaneously.

Further work is needed in order to complete the study on BaZrO 3 . It would be interesting to compute free energy barriers associated to the rotation step in order to evaluate the impact of quantum effects on these barriers and compare to the transfer step. Then, it would also be interesting to estimate values for the diffusion coefficient from these barriers using transition state theory. In order to achieve this, we will first need to compute values for the prefactor.

It could also be interesting to estimate the diffusion coefficient using ring-polymer molecular dynamics which allows the computation of a diffusion coefficient in the PIMD framework. One question remain concerning the calculation of free energy barriers using PIMD. Here, we have decided to compute the barrier associated with the centroid of the ring polymer because this is the quantity of interest in the Path Centroid Transition State Theory. However, this centroid free energy barrier is significantly different than the barrier computed using the replicas. The physical meaning of these two barriers and their differences remains unclear for now.

Concerning the study of GdBaCo 2 O 5.5 , the work done here could be extended by fitting the potential against different configurations obtained using density functional theory, in order to ensure that the force field accurately reproduces the potential energy surface of the proton in this material. Moreover, no quantum effects have been included in this study since there are not significant at the temperatures studied here. It could be interesting to study the impact of quantum effects at lower temperatures in this material, in particular to see what mechanism (proton transfer/rotation or OH migration) will be the most affected by the inclusion of quantum effects.

During this thesis, we have seen that several methods, such as PIMD or QTB, allow the computation of quantum time-independent average values (average energy, position distribution or pair correlation function for example). The computation of quantum time-dependent quantities however remains a computational challenge. In the path-integral framework, timedependent quantities are not directly accessible, and one need to use approximate methods such as ring-polymer molecular dynamics or centroid molecular dynamics to compute them.

The QTB method provides another framework to approximately compute quantum timedependent quantities, such as vibrational spectra or diffusion coefficients for example, with no additional computational cost compared to standard molecular dynamics simulations.

The QTB method has indeed been able to give satisfactorily results on the vibrational spectra of polyatomic molecules or high-pressure ice. However, it remains unclear to what extent the dynamics obtain with the QTB is a good approximation to the real quantum dynamics. In my opinion, this is the most fundamental question that remains on the QTB method. Thus, I think that it would be interesting to investigate this question for example by computing time-dependent quantities (such as correlation functions) on simple 1D systems using the QTB method and compare the results to other approximate methods, such as ring-polymer molecular dynamics or centroid molecular dynamics.

We have seen that the combined QTB-PIMD method reduces significantly the number of replicas needed for convergence, and thus allow a significant decrease of the computation time associated with those simulations. Nowadays, the interatomic forces are generally described using first-principle methods such as density functional theory (DFT). These methods are accurate but are also more computationally demanding than other methods such as force fields. PIMD is also computationally expensive since it requires one to simulate several replicas of the system. For this reason, PIMD and first-principle description of the forces are rarely used together. We believe that the combined QTB-PIMD will be particularly useful to make first principles PIMD studies more affordable. For future works, it would also be interesting to try to extend the combination to the case of centroid MD or ring polymer MD which allow the computation of time-independent quantities. Abstract : This thesis deals with the study of quantum effects in molecular dynamics (MD). MD is a powerful numerical method to investigate the properties of condensed matter systems. However, since the method is based on classical mechanics, quantum effects associated with the dynamics of the nuclei, such as zero-point energy or tunnelling, are not taken into account. These nuclear quantum effects can, however, play a major role in particular at low temperatures and/or in systems containing light atoms such as hydrogen. In these cases, a standard way to account for the quantum nature of the nuclei is to use path integral molecular dynamics (PIMD). This method provides exact quantum results however its high computational cost limits its range of applicability. The quantum thermal bath (QTB) method has been proposed as an alternative to PIMD. The QTB method is particularly appealing because of its computational cost that is equivalent to standard MD thus allowing to study large and complex systems. The first part of this thesis is devoted to the study of the QTB method. We have studied the behavior of the method in selected model systems in order to investigate its limitations. We have focused, in particular, on the zero-point energy leakage problem and found that increasing the friction coefficient of the QTB can significantly limit this problem. We also have developed another way to use the QTB method by combining it with PIMD. This combined QTB-PIMD method allows, in particular, to decrease the computational cost of standard PIMD simulations. In a second part, we have used these methods to study hydrogen conduction in perovskite materials. We have first investigated the impact of quantum effects on the diffusion of hydrogen in BaZrO 3 , a potential electrolyte material for hydrogen fuel cells. Since hydrogen is the lightest element, we expect quantum effects to have a significant impact on its dynamics. We find that quantum effects are indeed significant at low temperatures although their impact on the diffusion remains low at the typical working temperatures of hydrogen fuel cells. Finally, we have investigated the diffusion mechanisms of hydrogen in GdBaCo 2 O 5.5 . We evidence that the diffusion is anisotropic in this material and composed of two main diffusion mechanisms.

Université Paris-Saclay Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD 128 / 91190 Saint-Aubin, France

Figure 1 . 1 :

 11 Figure 1.1: Left: Schematic representation of the operating principle of a Proton Conducting Fuel Cell (PCFC), Right: General perovskite structure (ABO3). The B-site cations are surrounded by oxygen (O 2-) octahedra represented in light red here.

Figure 1 . 2 :

 12 Figure 1.2: Schematic illustration of proton motion in perovskite materials. The red points indicates a typical proton trajectory and the blue octahedra are the standard oxygen octahedra of the perovskite structure. The reorientation and transfer steps that result in long range diffusion of the proton are visible. Figure from reference [21].

Figure 1 . 3 : 1 . 4

 1314 Figure 1.3: Experimental diffusion coefficient for H, D and T in Nb. We clearly see the deviation from the Arrhenius law at low temperature due to quantum effects for T 250 K. Figure from reference [22].

Figure 2 . 3 :

 23 Figure 2.3: The probability density ρ(x) of finding the particle (hydrogen) at position x. Here x is in reduced units form x → x/σ with σ = m/ ω. For the PIMD calculation convergence was obtained for P ≈ 40. V(x) shows the form of the harmonic potential well.

Figure 2 . 5 :

 25 Figure 2.5: Position distribution of an hydrogen atom in the three double well potentials investigated here : C=0.1 (top panel), C=0.3 (middle panel) and C=1.0 (bottom panel). The simulations are carried out at a reduced temperature T * = kBT /V0 = 0.4. The horizontal gray line represent the average total energy of the system (computed by PIMD) to be compared to the height of the potential barrier V0 of the double-well V (y).

3 4 Quantum thermal bath for path integral molecular dynamics 4 . 1

 441 Limitations of the QTB method: zero-point energy leakage 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Coupled harmonic oscillators . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 One-dimensional chain of atoms . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Practical discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Complements of chapter 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Combining QTB and PIMD . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1 Modified power spectral density . . . . . . . . . . . . . . . . . . . . . . 4.2.2 Computation of macroscopic properties . . . . . . . . . . . . . . . . . 4.3 Results on simple one-dimensional systems . . . . . . . . . . . . . . . . . . . . 4.3.1 Harmonic oscillator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.2 Double-well potential . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Results on realistic systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.1 Ferroeletric-paraelectric phase transition in BaTiO 3 . . . . . . . . . . 4.4.2 Proton position distribution in BaZrO 3 . . . . . . . . . . . . . . . . . 4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Complements of chapter 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3 . 1 :

 31 Figure 3.1: Average energies ε1 and ε2 of the two oscillators, and average coupling energy εc computed by QTB-MD as a function of the intensity of the coupling constants c3 and c4. Top panel : cubic coupling only (c3 = 0, c4 = 0) with Ω = 0.5, Lower panel : quartic coupling only (c3 = 0, c4 = 0) with Ω = 0.25. Let us note that by symmetry εc = 0 in the cubic case.

Figure 3 . 2 :

 32 Figure 3.2: Zero-point energy leakage quantified using the deviation factor ζ (equation (3.5)) as a function of Ω. Top panel : cubic coupling only (c3 = 0, c4 = 0) and bottom panel : quartic coupling only (c3 = 0, c4 = 0))

  Figure 3.3 (b) shows the vibrational spectrum of the two oscillators for Ω = 0.2 and c 4 = 15.4 × 10 -4 , we see that many additional modes appears in the spectrum. Thus multiple resonances are likely to occurs leading to important ZPEL for Ω < 1/3.

oscillator 1 oscillator 2 Figure 3 . 3 :

 233 Figure 3.3: Vibrational spectra (logarithmic scale) of the two oscillators obtained using QTB-MD trajectories. Top panel : cubic coupling case with Ω = 0.5 and c3 = 2.4 × 10 -4 , bottom panel : quartic coupling with Ω = 0.2 and c4 = 15.4 × 10 -4 . In the two cases the results are shown for two values of γ : 4 × 10 -4 ω1 and 4 × 10 -3 ω1 to illustrates the effect of the friction coefficient on the spectrum.

Chapter 3 .

 3 Limitations of the QTB method: ZPEL

Figure 3 . 4 :

 34 Figure 3.4: Impact of the friction coefficient on the zero-point energy leakage. Top panel: cubic coupling only, Bottom panel: quartic coupling only. The values of γ are given on the graph in units of ω1.

1 Figure 3 . 5 :

 135 Figure 3.5: The computed normal mode frequencies of the OH chain as a function of the ν OO frequency (equation(3.8)). We clearly see that there is one low frequency ν2 which varies with ν OO and two very close high frequencies ν2 which are almost constant with ν OO

  7) with R the O-O distance, C 0 , α 0 the depth and the width of the potential well, respectively and R 0 the equilibrium O-O distance. The values of the parameters are set to C 0 = 3.81 eV and R 0 = 2.88 Å. The frequency associated with the O-O interactions can be computed in the harmonic approximation of V OO :

52 Chapter 3 .Figure 3 . 6 :

 52336 Figure 3.6: Effective temperature of H and O atoms computed using equation (3.9) and obtained from QTB-MD and PIMD simulations at T = 600 K. Here Ω = 0.5 and γ = 0.2 THz.

Figure 3 . 7 :

 37 Figure 3.7: Zero-point energy leakage evaluated through the deviation factor (equation (3.10)) as a function of the frequency ratio Ω and for different values of the friction coefficient γ.

54 Chapter 3 .

 543 Limitations of the QTB method:

Figure 3 . 8 :Figure 3 . 9 :

 3839 Figure 3.8: Distribution of O-H (left) and O-O(right) distances for a low and a high value of γ in the case Ω = 0.5 (top) and Ω = 0.1 (bottom)

Figure 3 . 10 :

 310 Figure 3.10: Top panel : relative deviation of the QTB-MD kinetic ( K ) and potential ( V ) energies obtained from their exact values (θ(ω0, T )/2) for an harmonic oscillator as a function of the friction coefficient γ. K vir is the kinetic energy obtained from the virial theorem (eq. 3.17). Bottom panel : Fourier transform of the kinetic energy K(ω) (eq.3.15) and potential energy V (ω) (eq.(3.16)) for different values of γ.

Figure 3 . 11 :

 311 Figure 3.11: Kinetic energy distribution of an aluminium crystal at T = 10 K as obtained from QTB-MD with a friction coefficient γ = 0.9 THz (full circles) and with γ = 10 THz (open circles). The angular frequency cut-off is ωcut = 2ωmax with ωmax/2π = 10 THz. The solid line corresponds to the exact distribution θ(ν, T ). The inset present the evolution of the slope of the energy distribution as a function of γ.

Figure 3 . 12 :

 312 Figure 3.12: Temperature evolution of the reduced polarization of BaTiO3 as obtained from QTB-MD with a friction coefficent γ = 16 THz and from PIMD with a product P × T = 1920 which corresponds to P = 16 at T = 120 K. Vertical dashed lines indicates the temperatures of the different transitions. The inset presents the evolution with γ of the different transitions temperatures obtained from QTB-MD simulations. The dashed grey lines indicates the values obtained by PIMD.

1Chapter 4 .

 4 The quantum thermal bath method is exact in the harmonic case Ph.D. thesis -F.Brieuc -2016 68 The QTB-PIMD method So the quantum thermal bath can be applied to the normal mode coordinates of the ring polymer {q k }. Using a random force with power spectral density I R = 2mγκ(ω, T ) and κ(ω, T ) obtained from equation (4.13) will ensure the correct quantum position fluctuation x 2 for an harmonic oscillator (eq. (4.3)) and thus the correct average potential energy

Finally, let us

  summarize the main steps of the QTB-PIMD method :• We first solve equation(4.11) or (4.20) to obtain the functions f

26 )

 26 can easily show that q2 0 = P ẋ2 c and we obtain 1Ph.D. thesis -F.Brieuc -2016 So if one modify the standard centroid virial estimator by replacing 1/2β by P1 2 m ẋ2 c the relation K = U is directly ensured even in the f (0)

Figure 4 . 1 :Figure 4 . 2 :

 4142 Figure 4.1: Convergence of the average potential and kinetic energy as a function of the Trotter number P for PIMD and QTB-PIMD simulations using the f (0) P and f (1) P functions referenced here as "(0)" and "(1)" respectively. The Morse potential with 1/λ 2 = 0.0015 is shown on panel (a) and (b) and with 1/λ 2 = 0.024 is shown on panel (c) and (d). The energies are normalized by the well depth D and calculations have been carried out at T = 0.02D/kB . The exact values have been obtained by numerically solving the Schrödinger equation.

( 1 )

 1 P function and P = 1 we are in the classical case and thus we find as expected that the average total energy is k B T .

Figure 4 . 3 :

 43 Figure 4.3: Kinetic and potential energy convergence with respect to the number of beads P otbained by QTB-PIMD and PIMD for an harmonic oscillator of angular frequency ω. The energies are given in units of ω. The gray line is the exact quantum result for comparison. Error bars are smaller than the symbols size and thus are not included here.

Figure 4 . 4 :

 44 Figure 4.4: Divergence factor (dP ) as a function of the number of beads for the double-well potential. Three cases are investigated : (a) C = 1.0, (b) C = 0.3 and (c) C = 0.1. The arrows indicates the smallest number of beads obtained for dP < 2%. The position distributions corresponding to the QTB-PIMD case are shown (blue solid line) for P = 16 (d), P = 7 (e) and P = 4 (f) together with the exact solution (black dashed line) and the distibution obtained with the bare QTB-MD (P = 1) (red dotted line). The V0 values are deduced from the expression of C and the distance between the well is fixed at 2a = 0.8Å . The simulations are carried out at a reduced temperature T * = kBT /V0 = 0.4.

Figure 4 .

 4 Figure 4.4 shows the position probability density ρ(y) obtained by QTB-PIMD and Standard PIMD (right panel). The simulations are carried out at a reduced temperature T * = k B T /V 0 = 0.4. The convergence of the density probability with the number of beads P is evaluated by calculating the divergence factor

Figure 4 . 5 :

 45 Figure 4.5: Convergence of the total energy with the Trotter number for the three double well potentials (C = 1.0, 0.3, 0.1) . Comparison between QTB-PIMD with f(1) P (blue points) and Standard PIMD (red square). The kinetic part of the energy is computed using the modified centroid virial estimator for QTB-PIMD and using the standard centroid virial estimator for Standard PIMD.

Figure 4 . 6 :

 46 Figure 4.6: Schematical view of the structure of barium titanate (BaTiO3) in the cubic and the tetragonal phase. The general structure of BaTiO3 is the perovskite structure with Ba 2+ cations on the A-site and Ti 4+ on the B-site, finally O 2-form the standard oxygen octahedras (indicated in ligth red on the figure) surrounding the Ti 4+ ions. Transition from cubic to tetragonal phase is associated to an atomic displacement (indicated by the arrows) that breaks the inversion symetry. The positions of the barycenter of the positive and negative charges are different which leads to the apparition of an electrical dipole P in the unit cell : the material thus becomes ferroelectric. So the cubic-tetragonal phase transition is also called the paraelectric-ferroelectric phase transition.

Figure 4 . 7 :

 47 Figure 4.7: Evolution of the reduced polarization with temperature in BaTiO3. Here we focus on the transition between the tetragonal-ferroelectric and the cubic-paraelectric phases. The polarization is obtained by PIMD (with P = 16) and QTB-PIMD using either f (0) P and f

Figure 4 . 8 :

 48 Figure 4.8: Position distribution of the proton in the Zr-O-Zr plane at T = 300 K obtained using PIMD and QTB-PIMD simulations. The top image displays an example of the distribution in the unit cell (PIMD with P = 1), where one can distinguish the eight equivalent positions for the proton. Let us note that the distribution has been symmetrized using the symmetry of the unit cell in order to obtain the eight equivalent positions since proton transfer is rare at this temperature. The other images are enlargements of the distribution.

Figure 4 . 9 :

 49 Figure 4.9: Divergence factor dP as a function of the Trotter number P . We see that PIMD is converged within an error of dP = 8% for P = 6 while QTB-PIMD is converged within the same error for P = 4. The calculations are performed using the f (0) P function. Two values of dP obtained with the f (1) P function are also added to illustrate that the convergence obtained with both f (0) P or f (1) P function is similar in this case.

P

  (u) is obtained from the function F P (u):f (0) P (u) = F P u √ P (4.42)Then we generate the stochastic forces from the power spectral densityI R (ω, T ) = 2mγκ(ω, T )Ph.D. thesis -F.Brieuc -

( 1 )

 1 P (u) cannot be generated for values of u lower than √ P sin π P . For lower values of u we arbitrarily put the function f (1) P (u) to zero which does not have a significant impact on the final result in practice. The functions f (1) P obtained for different values of P are presented in figure 4.10right panel. We see that the functions approaches u/[ √ P (P -1)] for u → ∞.

Figure 4 . 10 :

 410 Figure 4.10: Example of fP functions for some values of P . Left: f (0) P as a function of u * = u/P 3/2 , Right:

  • around the Zr-O-Zr axis and thus changes to another O-Zr-O plane, while in the second rotation (R 2 ), the proton rotates around an axis orthogonal to the Zr-O-Zr axis and remains in the same O-Zr-O plane. The second reorientation process (R 2 ) is a very fast process and thus is not limiting for the long range migration of the proton. For this reason, this rotation is generally discarded of the analysis and the proton is globally considered to be free to rotate in the same O-Zr-O plane. The proton distribution obtained from classical Langevin MD at relatively low temperature (T = 300 K) is shown in figure5.2: we see that the probability of finding the proton between the two equilibrium positions associated with rotation R 2 is high. The rotation of the proton in the same O-Zr-O plane (R 2 ) is indeed a fast process and thus is not limiting for the proton diffusion.

Figure 5 . 1 :

 51 Figure 5.1: Representation of the elementary processes for the long range diffusion of protons in barium zirconate. There is one transfer step (T) and two types of reorientations (R1 and R2).

2 Figure 5 . 2 :

 252 Figure 5.2: Proton distribution in the O-Zr-O plane calculated from standard(Langevin) MD simulations at T = 300 K. Let us note that the distribution has been symmetrized using the symmetry of the unit cell in order to obtain the eight equivalent positions since proton transfer is rare at this temperature.

Figure 5 . 3 :

 53 Figure 5.3: Schematic representation of the two states used in the EVB model for the transfer step. The thick line represent the covalent bond and the dashed line represent the hydrogen bond. Let us emphasize that ths is is only a schematic representation. In particular, the proton is generally not on the O1-O2 line in practice.

Figure 5 . 4 :

 54 Figure 5.4: Schematic representation of the potential energy within the O-H-O group along the transfer reaction coordinate ξ = r2 -r1 as obtained with the EVB model for a typical value of the O1-O2 distance (3.2 Å).The red curve is the potential energy for the proton in the |1 state and the green is the potential energy for the proton in the |2 state. Finally the blue curve is the coupling between the two states. We see that the EVB construction gives a total energy with a symmetric double well shape (dashed black line).

Figure 5 . 5 :

 55 Figure 5.5: Relaxed energy (T = 0K) per unit cell as a function of the unit cell volume for pure barium zirconate. We see that the Murnaghan equation of state fits the results and we obtain an estimation of the equilibrium unit cell volume Ω0 ≈ 73.6Å 3 and the bulk modulus B0 ≈ 213.7 GPa.

Figure 5 . 6 :

 56 Figure 5.6: Position distribution of a proton in BZO at T= 1500 K as obtained from Langevin MD. We clearly see the two types of elementary processes: rotation (red arrows) and transfer (black arrows). Here only rotation R1 is shown and we see that indeed the proton is free to rotate in the O-Zr-O plane (rotation R2). See also figure 5.1.

Figure 5 .

 5 Figure 5.6 shows the proton position distribution in BZO at high temperature (T = 1500 K).One can see the two types of elementary processes, rotation and transfer, described in section 5.2. The arrows on the figure indicate a possible diffusion pathway for the proton and we see that a combination of rotation and transfer is necessary for long range diffusion. In order to study these processes more precisely we compute the static potential energy barriers associated with these two processes. Figure5.7 shows the static potential energy barriers calculated along the reaction coordinates associated with the rotation (R 1 ) and transfer steps. These barriers are obtained by relaxing the system with the reaction coordinate constrained to a certain value. For the transfer step, we choose a standard form for the reaction coordinate ξ( r 1 , r 2 ) = r 2r 1 with r 1 being the O 1 -H distance and r 2 the O 2 -H distance. For the rotation step, we define the reaction coordinate θ as the dihedral angle between the initial O 1 -Zr-O 2 plane and the O 1 -Zr-H plane. The constraints have been included using the SHAKE algorithm[28, 29] (see the complements of the chapter for more details). The energy barrier associated with the transfer is higher than the one associated with rotation, indicating that the transfer step might be the limiting step for the diffusion of the proton as it is generally accepted[2]. These energy barrier values (0.38 eV for the transfer and 0.29 eV for the rotation) are higher than the typical values reported in the litterature[9, 17, 18]. Sundell and coworkers[17] for example report values of around 0.2 eV for the transfer step and around 0.18 eV for the rotation step. However, higher values have also been reported for example Gomez and coworkers[12] report transfer barriers of around 0.4-0.5 eV far from the dopant in Y-doped barium zirconate 2 . In their paper, Raiteri and coworkers surprisingly report a static barrier for the transfer of around 0.038 eV (We think that it is a typing error).

Figure 5 . 7 :

 57 Figure 5.7: Static potential energy associated with the two elementary processes of diffusion : transfer (left panel) and rotation (right panel).

Figure 5

 5 Figure 5.9: Jump frequencies for the rotation and transfer steps as obtained from Langevin MD. The two type of frequencies exhibit a linear behavior (Arrhenius behavior) which is expected since they are related to thermally activated processes. Frequencies for rotations are at least one order of magnitude higher than the transfer frequencies on this temperature range thus the transfer step is the limiting step.

  D 0 (cm 2 /s) 3.5 × 10 -4 2.47× 10 -3 2.77× 10 -4

Figure 5 . 10 :

 510 Figure 5.10: Diffusion coefficient as obtained from the MSD (eq. (5.9)) and from the jump frequencies (eq. (5.10)) compared to the results of Raiteri et al. modified to be consistent with the correct jump length d = a/ √ 2.

Figure 5 . 11 :

 511 Figure 5.11: Diffusion coefficient, D(T ), as obtained from the MSD using standard (Langevin) MD and QTB-MD.

Figure 5 . 12 :

 512 Figure 5.12: Bottom panel: kinetic energy distribution among the modes as obtained from QTB-MD (green line and circles) and compared to the expected distribution θ(ν, T ) (black line) zoom around the low frequency modes -Middle panel: vibrational density of state (DOS) -Top panel: zoom of the DOS as obtained from QTB-MD (green line) and from standard MD (purple line) around the low frequency region (left) and around the high frequency (OH) modes (right). The DOS are obtained from the Fourier transform of the velocity autocorrelation function normalised by kBT in the case of Standard (Langevin) MD and by θ(ν, T ) for QTB-MD.

Figure 5 . 13 :

 513 Figure 5.13: Schematic representation of a typical free energy profile ∆F (ξ) along a reaction path described by the reaction coordinate ξ({ r}). Quantities associated to the saddle point (barrier top) are represented with a ‡ superscript and quantities associated with the bottom of the well are represented with a 0 superscript.

Figure 5 . 14 :Figure 5 . 15 :

 514515 Figure 5.14: Free energy profiles as obtained using standard MD with constrained reaction coordinate ξ (top panel) and PIMD with constrained centroid reaction coordinate ξc(bottom panel).

Figure 5 . 16 :

 516 Figure 5.16: Replicas probability density ρ(ξ) with the centroid constrained at the barrier top i.e. ξc = ξ ‡ = 0.

Figure 5 . 17 :

 517 Figure 5.17: Evolution of the free energy barrier (-β∆F ‡ ) as obtained from PIMD for hydrogen and deuterium as a function of 1000/T . The Standard MD results are also plotted fro comparison purposes.

Figure 5 . 18 :

 518 Figure 5.18: Replicas probability density ρ(ξ) with the centroid constrained at the barrier top i.e. ξc = ξ ‡ = 0 for Hydrogen (left) and deuterium (right)

5 .

 5 NQE on the proton conduction in BZONow in the long time limit the Einstein relation becomes

Figure 5 . 19 :

 519 Figure 5.19: Representation of the reaction coordinate used to describe the reaction path associated with the rotation step

ξ ( 5

 5 .80) Where < ... > const. ξ means the average value computed from constrained MD simulations. λ is the Lagrange multiplier of the constrained MD algorithm (see the previous complement). And the two terms Z and G are defined as follows

  m O , m H the mass of the oxygen and the hydrogen atoms, respectively.

Figure 5 . 20 :

 520 Figure 5.20: Top: Free energy profile along the reaction coordinate ξ = r2 -r1 as obtained from unconstrained strandard MD using the probability density ρ(ξ) -equation (5.72) (red full circles), and as obtained from constrained MD and computed using thermodynamic integration -equation (5.73) (blue line and cross) at T = 300 K. Bottom: ξc = r2,c -r1,c as obtained from unconstrained PIMD using the probability density ρ(ξc) -equation (5.72) (red full circles), and as obtained from constrained PIMD and computed using thermodynamic integration -equation (5.73) (blue line and crosses). For comparison purposes, the free energy profile along the reaction coordinate of the replicas ξ as obtained from unconstrained PIMD using the probability density ρ(ξ) is also plotted (green triangles)

Figure 6 . 1 :

 61 Figure 6.1: Representation of the double perovskite structure of GBCO. Note that, in the real materials, oxygen vacancies are present and the GdO planes, in particular, are oxygen deficient. This is not represented here.

Figure 6 .Figure 6 . 3 :Figure 6 . 4 :Figure 6 . 5 :

 6636465 Figure 6.3 shows the mean-square displacement (MSD) of the proton as obtained from simulations at 1700 K. The proton diffuses in parallel to the (a,b) plane and is blocked along the c direction, similarly to oxygen atoms. In order to elucidate the origin of this anisotropic diffusion, one can first study the preferential locations of the proton in the material. Figure6.4 shows the probability density map of a proton in GBCO as obtained from MD trajectories. The proton diffuses within the same GdO and CoO 2 planes which explains the bidimensional nature of proton diffusion parallel to the (a,b) plane. The proton is never found in the BaO plane: similarly to oxygen diffusion, BaO planes act as barriers for hydrogen diffusion and thus hinder long range migration of the proton along the c-axis. Figure6.5 presents the position distribution (or density profile) of the atoms along the c-axis, in order to get more

Figure 6 .

 6 Figure 6.6 shows the diffusion coefficient in the (a,b) plane for both protons and oxygen atoms as obtained from mean square displacement. Comparison is made with previous experimental and theoretical results, which only exist for oxygen atoms. The calculated oxygen diffusion coefficient values are in reasonable agreement with experimental results[2, 23]. The diffusion coefficient is found to be more than one order of magnitude higher for hydrogen than for oxygen atoms. The activation energy associated with proton diffusion (0.71 ± 0.03 eV) is found to be lower than the activation energy for oxygen diffusion (1.01 ± 0.03 eV). This difference of activation energy between oxygen and proton migration is typically observed in perovskite ion conducting compounds. For instance, in yttrium-doped barium zirconate, an activation energy of around 1 eV[24] is reported for oxygen ion migration, while it is around

Figure 6 . 6 :

 66 Figure 6.6: Arrhenius plot of the diffusion coefficient in the (a,b) plane as obtained from MSD for both protons (blue squares) and oxygens (red full circles) atoms. (a) from MD simulations -ref.[9], (b) from oxygen exchange -ref.[2], (c) from conductivity relaxation -ref.[23]

Figure 6 .Figure 6 . 7 :Figure 6 . 8 :

 66768 Figure 6.7: Schematic representation of the different types of diffusion in GBCO. Mechanism (1) is the transfer of the proton which combined with the rotation leads to a long range migration of the proton jumping from one oxygen of the Gd-plane to an oxygen of the Co-plane. The second mechanism(2) is the migration of the whole OH group via a vacancy mechanism. Finally, the diffusion of oxygen via oxygen vacancies is also indicated.

Titre:

  Modélisation et simulation des effets quantiques en dynamique moléculaire: application à l'étude de la conduction protonique Mots clefs : dynamique moléculaire, effets quantiques, conduction protonique Résumé : Cette thèse porte sur l'étude des effets quantiques en dynamique moléculaire (DM). La DM est une méthode numérique qui permet l'étude des propriétés de la matière condensée. Cependant, la méthode étant basée sur la mécanique classique, les effets quantiques associés à la dynamique des noyaux, tels que l'énergie de point zéro ou l'effet tunnel, ne sont pas pris en compte. Ces effets quantiques nucléaires peuvent cependant jouer un rôle majeur, en particulier aux basses températures et/ou dans les systèmes contenant des atomes légers comme l'hydrogène. La dynamique moléculaire par intégrales de chemins (PIMD) est souvent utilisée, dans ce cas, pour tenir compte de la nature quantique des noyaux. Cette approche fournit des résultats quantiques exacts, mais son coût en temps de calcul élevé limite son domaine d'application. La méthode du bain thermique quantique (QTB) a été proposée comme une alternative à la PIMD. L'approche QTB est particulièrement intéressante car son coût en temps de calcul est équivalent à celui de la DM standard permettant ainsi l'étude de systèmes complexes et de plus grande taille. La première partie de cette thèse est consacrée à l'étude de la méthode QTB. Nous avons étudié le comportement de la méthode sur différents sys-tèmes modèles afin d'étudier ses limites. En particulier, le problème du "zero point energy leakage" est étudié en détail et nous montrons que l'augmentation du coefficient de friction du QTB permet de limiter ce problème. Nous avons également développé une combinaison de la méthode QTB avec la méthode PIMD. Cette méthode combinée QTB-PIMD permet de réduire le coût en temps de calcul des simulations PIMD standard. Dans une deuxième partie, nous avons utilisé ces méthodes pour étudier la conduction de l'hydrogène dans des matériaux pérovskites. Nous nous intéressons d'abord à l'impact des effets quantiques sur la diffusion de l'hydrogène dans BaZrO 3 , un matériau d'électrolyte potentiel pour piles à hydrogène. L'hydrogène étant l'élément le plus léger, un impact important des effets quantiques est attendu. Nous trouvons que les effets quantiques sont effectivement importants à basse température, mais leur impact sur la diffusion reste faible aux températures de fonctionnement typiques des piles à hydrogène. Enfin, nous avons étudié les mécanismes de diffusion de l'hydrogène dans GdBaCo 2 O 5.5 . Nous mettons en évidence une diffusion anisotrope dans ce matériau et deux mécanismes principaux de diffusion.Title : Modelling and simulation of quantum effects in molecular dynamics: application to the study of proton conduction Keywords : molecular dynamics, quantum effects, protonic conduction

Table 2

 2

	.1: Average energy obtained from QTB-MD and PIMD simulations averaged over 10 different trajec-
	tories of time duration Nsδt = 5000/γ with δt = 0.01/νmax and γ = 0.1νmin. The exact energy can easily be
	computed and is 0.25 at this temperature. The total energy is 0.5 ( ω/2).

Table 2 .

 2

	2: Average energy obtained from QTB-MD and PIMD simulations averaging over 10 different trajec-
	tories of time duration Nsδt = 5000/γ with δt = 0.01/νmax and γ = 0.1νmin.

Ph.D. thesis -F.

Brieuc -

2016 2.5. Results on one-dimensional systems 27

			Classical
		PIMD (P=30)
	ρ(x)		QTB
	probability density		V(x)
	-0.8 -0.6 -0.4 -0.2	0	0.2 0.4 0.6 0.8
		x (Å)
			Classical
		PIMD (P=40)
	ρ(x)		QTB
	probability density		V(x)
	-0.8 -0.6 -0.4 -0.2	0	0.2 0.4 0.6 0.8
		x (Å)
			Classical
		PIMD (P=50)
	ρ(x)		QTB
	probability density		V(x)
	-0.8 -0.6 -0.4 -0.2	0	0.2 0.4 0.6 0.8
		x (Å)
	Figure 2.4: Position probability density obtained using QTB-MD and PIMD for the three cases : 1/λ 2 = 0.0015
	(top panel), 1/λ 2 = 0.004 (middle panel) and 1/λ 2 = 0.024 (bottom panel). For the PIMD simulations, the
	number of replicas necessary to converge is P = 30 in the weakly anharmonic case (1/λ 2 = 0.0015), P = 40
	in the midly anharmonic case (1/λ 2 = 0.004) and P = 50 for the highly anharmonic case (1/λ 2 = 0.024).
	The simulations are carried out at a temperature T = 0.02D/kB . The form of the Morse potential V (x) is
	also shown.

Chapter 2. NQE in MD simulations

  2 the energy associated with the interaction between replicas and U = P s=1 V (x s )/P the potential energy. Now we are interested in calculating the following expression

	so that
				dx 1 ... dx P		P s=1	x s	∂φ P ∂x s	exp (-βφ P ) = -	1 β	dx 1 ... dx P	P s=1	x s	∂ ∂x s	e -βφ P (2.105)
	which can be integrated using an integration by parts
									dx s x s	∂ ∂x s	e -βφ P = x s e -βφ P +∞ -∞	-dx s e -βφ P	(2.106)
	Now if we consider that
									x s e -βφ P +∞ -∞	= 0 because	lim xs→±∞	e -βφ P = 0
	we obtain
	P s=1	x s	∂φ P ∂x s	=	1 β	1 Z P		P s=1	1 (2π ) P s=1 x s	∂φ P ∂x s	=	1 β	1 Z P	P s=1	Z P =	P β	(2.108)
									P s=1	x s	∂φ P ∂x s	=
	P s=1	x s	∂φ P ∂x s	=	1 Z P	1 (2π ) P s=1	p 2

P dp 1 ... dp P exp -β s 2m dx 1 ... dx P P s=1 x s ∂φ P ∂x s exp (-βφ P ) (2.103) One can notice that x s ∂φ P ∂x s e -βφ P = -1 β x s ∂ ∂x s e -βφ P (2.104) Ph.D. thesis -F. Brieuc -2016 P dp 1 ... dp P exp -β P s=1 p 2 s 2m dx 1 ... dx P exp (-βφ P ) (2.107) Now if we split the expression of φ P in contributions from replicas α P and from potential energy U (see eq. (2.102)) we obtain P s=1

  1 and Ba-O 2 are interchanged, same for Zr-O 1 and Zr-O 2 or O-O 1 and O-O 2 . Thus the last term U

′ represents all these other interactions (i.e. coulombic and Buckingham type interactions) with the role of O 1

Table 5 .

 5

1: Equilibrium cell parameter a0 and bulk modulus B0 values obtained from the Murnaghan fit and comparison with values obtained experimentally and by Raiteri and coworkers.

26 eV and 98 Chapter 5. NQE on the proton conduction in BZO

				T (K)
		1500 1200 1000	800	700	600
	diffusion coefficient (cm 2 /s)	1e-07 1e-06 1e-05			Raiteri et al. (6x6x6)
		0.6	0.8	1	1.2	1.4	1.6
				1000/T (K -1 )

  For the rotation step, we have chosen to define the reaction coordinate as the dihedral angle θ between the plane [ e x , d O 1 ,H ] and the plane [ e x , e z ] -see figure 5.19. With d O 1 ,H = r O 1r H . θ is given by the angle between the normal n to the plane [ e x , d O 1 H ] and e z . The normal n is obtained using the cross product between d O 1 H and e x :

	n =	d O 1 H × e x = Z e y -Y e z 1 √ Y 2 + Z 2 (Z e y -Y e z )	(5.60) (5.61)

  3N dp 1 ... dp N dx 2 ... dx N 3N dp 1 ... dp N dx 1 ... dx N

							5.75)
	∂F ∂ξ	=	1 Z(ξ)	1 h ∂E ∂ξ	e -βE	(5.76)
	∂F ∂ξ	=	1 Z(ξ)	1 h ∂E ∂ξ	e -βE δ (x 1 -ξ)	(5.77)
	∂F ∂ξ	=	∂E ∂ξ ξ			(5.78)
							be
	expressed as:	∂F ∂ξ	=	∂V ∂ξ ξ	(5.79)

Where < ... > ξ means an ensemble average with fixed value of ξ({ r}). Finally, since the reaction coordinate is a function of atomic positions only, the free energy derivative can

  the distance between H and the two first neighbor oxygens O 1 and O 2 . In this case, one finds the following expressions for Z and G[START_REF] Komeiji | Implementation of the blue moon ensemble method[END_REF]:

Table 6

 6

	.1.

Table 6 .

 6

1: Buckingham potential parameters used for GdBaCo2O5.5. Nomimal charges are used so that the charges for cations are +3, +2 and +3 for Gd, Ba and Co, respectively. O refers to oxygen ions, while O1 refers to an oxygen ion bonded to a proton. The charge of O and O1 are -2 and -1.308698 and the charge of H is 0.308698 so that the overall charge for the OH group is -1. The parameters have been taken from : Gd-O

[9]

, Ba-O, Co-O and O-O

[9, 18]

.

Table 6 . 2 :

 62 Fraction of the residence time of H in BaO, GdO and CoO2 planes.

Table 6 .

 6 3: Frequencies of H jumps corresponding to transfer from the first to the second neighboring oxygen atoms. Transfer can occur between two oxygen atoms located, in Ba and Co planes (Ba↔Co), in Gd and Co planes (Gd↔Co), or in the same Co plane (Co↔Co). Frequency for the jump of the whole OH group is also given in the last column.

	6.4. Correlation factor

Ph.D. thesis -F.

Brieuc -2016

Kröger-Vink notation : M c s -M represents the species which can be an atom (Ba,Ti, Zr...), a vacancy V, an electron e ... c is the electronic charge with a • meaning one positive charge, a ′ one negative charge and × means neutral. Finally s is the lattice site that the species occupies.Ph.D. thesis -F.Brieuc -2016

Ph.D. thesis -F.Brieuc -2016

These effects are indeed negligible as long as delocalization is not too important ie. as long as the thermal De Broglie wavelength of the particles remains small compare to the characteristic distance between these particles. It will always be the case in the following.

In particular when used in combination with first principles computation of interatomic forces.

δ n ′ ,0 = 1 if n ′ = 0, δ n ′ ,0 = 0, if n ′ = 0 Ph.D. thesis -F.Brieuc -2016

IX = 1 T limT →∞ | XT (ω)| 2 with XT (t) = X(t) for 0 < t < T ; XT (t) = 0 otherwise Ph.D. thesis -F.Brieuc -2016

ΘD = ω k B Ph.D. thesis -F.Brieuc -2016

which is possible in the case of bosons[1] but more problematic for fermions[4]

Average values can also be computed using Monte Carlo techniques which leads to methods called Path Integral Monte Carlo (PIMC) Ph.D. thesis -F.Brieuc -2016

Let us note that other type of thermostats are also used, in particular Nose-Hoover chains thermostat.Ph.D. thesis -F.Brieuc -2016

Ph.D. thesis -F. Brieuc -20162.5. Results on one-dimensional systems

This is done in practice by solving the time independent Schrödinger equation using a finite differences method that allow us to rewrite the problem as a matrix equation. Then the eigenenergies and eigenvectors are obtained through a matrix diagonalization. Finally, the effect of temperature is accounted for by occupying the resulting eigenstates with the proper Boltzmann factor.Ph.D. thesis -F.Brieuc -2016

Ph.D. thesis -F. Brieuc -20163.6. Complements of chapter 3

In this chapter and the next one, proton and hydrogen are used equivalently and both referred to the H + ion.

It is worth noting that all these values of the barriers have been obtained using Density Functional Theory within the Generalized Gradient Approximation which is known to underestimate hydrogen reaction barriers[30].Ph.D. thesis -F.Brieuc -2016

We have tested this assertion and we indeed find that, on average, Dx = Dy = Dz.Ph.D. thesis -F.Brieuc -2016

With r 1/2,c the distance between the centroid of the H atom and the centroid of the O 1/2 oxygen atoms Ph.D. thesis -F.Brieuc -2016

holonomic constraint: that can be written in the form g({ r}) = 0 i.e. only dependent on atomic positions.Ph.D. thesis -F.Brieuc -2016

Solid Oxide Fuel Cells (SOFC) here referred to hydrogen fuel cells with solid oxygen conductor electrolyte materials as opposed to Proton Conducting Fuel Cells (PCFC) for which the electrolyte is a proton conductor material.

Here, we don't include quantum effects in the study since we are working at temperatures higher than 1000 K for which quantum effects are negligible.Ph.D. thesis -F.Brieuc -2016

For example, some energy barriers values associated with the transfer/rotation step reported in the literature for different barium based perovskite materials are 0.25/0.19 eV for cubic BaTiO3[20], 0.32/0.21 eV for BaSnO3[21] while we found 0.38/0.29 eV for BaZrO3Ph.D. thesis -F.Brieuc -2016

first-principles calculations", Solid State Ionics, vol. 176, no.39-40, pp. 3035-3040, 2005. Ph.D. thesis -F.Brieuc -2016

Remerciements

The layout of this manuscript is based on the Ph.D. thesis L A T E X template created by Olivier Commowick available on http://olivier.commowick.org/thesis_template.php.

Most of the figures have been made using the following free softwares: gnuplot, VESTA and inkscape.

Chapter 5. NQE on the proton conduction in BZO

So we finally obtain the following expression for δλ (1) : δλ (1) (t) ≈ g r

(1) 1 , ..., r

(1) (1) (t + δt)}) (5.39) Because this expression for δλ (1) (t) is an approximation (and to the first order), the procedure needs to be iterated until convergence is reached. The algorithm thus computes δλ (1) (t), δλ (2) (t), ... until the constraint is fulfilled with a precision that is under a target precision.

In practice, the algorithm first check if the constraint is fulfilled with the positions r

(1) i (t+δt). If yes, then these new positions are accepted and we start the next time step, if not, then the algorithm computes δλ (1) and then new positions r

). The algorithm then check if the constraint is fulfilled with the positions r

(2) i (t + δt). If yes, these new positions are accepted and we start the next time step, if not, the iterated procedure continues until convergence is reached. At the end the Lagrange multiplier is λ = λ (1) + δλ (1) + δλ (2) + ....

Transfer step For the transfer step the reaction coordinate is

The gradients of the constraint are:

And the forces that have to be applied are given by:

While the expression for δλ (1) and then λ is obtained from equation (5.39).

In order to compute free energy barrier using PIMD we constrain the centroid reaction coordinate defined as Ph.D. thesis -F. Brieuc -2016